solutions/chap08soln.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples and Exercises from Think Stats, 2nd Edition\n",
    "\n",
    "http://thinkstats2.com\n",
    "\n",
    "Copyright 2016 Allen B. Downey\n",
    "\n",
    "MIT License: https://opensource.org/licenses/MIT\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import brfss\n",
    "\n",
    "import thinkstats2\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The estimation game\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Root mean squared error is one of several ways to summarize the average error of an estimation process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RMSE(estimates, actual):\n",
    "    \"\"\"Computes the root mean squared error of a sequence of estimates.\n",
    "\n",
    "    estimate: sequence of numbers\n",
    "    actual: actual value\n",
    "\n",
    "    returns: float RMSE\n",
    "    \"\"\"\n",
    "    e2 = [(estimate-actual)**2 for estimate in estimates]\n",
    "    mse = np.mean(e2)\n",
    "    return np.sqrt(mse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function simulates experiments where we try to estimate the mean of a population based on a sample with size `n=7`.  We run `iters=1000` experiments and collect the mean and median of each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 1\n",
      "rmse xbar 0.3666916796469523\n",
      "rmse median 0.4440447654801178\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "def Estimate1(n=7, iters=1000):\n",
    "    \"\"\"Evaluates RMSE of sample mean and median as estimators.\n",
    "\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for _ in range(n)]\n",
    "        xbar = np.mean(xs)\n",
    "        median = np.median(xs)\n",
    "        means.append(xbar)\n",
    "        medians.append(median)\n",
    "\n",
    "    print('Experiment 1')\n",
    "    print('rmse xbar', RMSE(means, mu))\n",
    "    print('rmse median', RMSE(medians, mu))\n",
    "    \n",
    "Estimate1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using $\\bar{x}$ to estimate the mean works a little better than using the median; in the long run, it minimizes RMSE.  But using the median is more robust in the presence of outliers or large errors.\n",
    "\n",
    "\n",
    "## Estimating variance\n",
    "\n",
    "The obvious way to estimate the variance of a population is to compute the variance of the sample, $S^2$, but that turns out to be a biased estimator; that is, in the long run, the average error doesn't converge to 0.\n",
    "\n",
    "The following function computes the mean error for a collection of estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MeanError(estimates, actual):\n",
    "    \"\"\"Computes the mean error of a sequence of estimates.\n",
    "\n",
    "    estimate: sequence of numbers\n",
    "    actual: actual value\n",
    "\n",
    "    returns: float mean error\n",
    "    \"\"\"\n",
    "    errors = [estimate-actual for estimate in estimates]\n",
    "    return np.mean(errors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function simulates experiments where we try to estimate the variance of a population based on a sample with size `n=7`.  We run `iters=1000` experiments and two estimates for each sample, $S^2$ and $S_{n-1}^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean error biased -0.11357007322862302\n",
      "mean error unbiased 0.03416824789993983\n"
     ]
    }
   ],
   "source": [
    "def Estimate2(n=7, iters=1000):\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    estimates1 = []\n",
    "    estimates2 = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        biased = np.var(xs)\n",
    "        unbiased = np.var(xs, ddof=1)\n",
    "        estimates1.append(biased)\n",
    "        estimates2.append(unbiased)\n",
    "\n",
    "    print('mean error biased', MeanError(estimates1, sigma**2))\n",
    "    print('mean error unbiased', MeanError(estimates2, sigma**2))\n",
    "    \n",
    "Estimate2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean error for $S^2$ is non-zero, which suggests that it is biased.  The mean error for $S_{n-1}^2$ is close to zero, and gets even smaller if we increase `iters`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The sampling distribution\n",
    "\n",
    "The following function simulates experiments where we estimate the mean of a population using $\\bar{x}$, and returns a list of estimates, one from each experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SimulateSample(mu=90, sigma=7.5, n=9, iters=1000):\n",
    "    xbars = []\n",
    "    for j in range(iters):\n",
    "        xs = np.random.normal(mu, sigma, n)\n",
    "        xbar = np.mean(xs)\n",
    "        xbars.append(xbar)\n",
    "    return xbars\n",
    "\n",
    "xbars = SimulateSample()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the \"sampling distribution of the mean\" which shows how much we should expect $\\bar{x}$ to vary from one experiment to the next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYFNXZ/vHvA8Owqwi4AQIqJqLiNiouURNc0DcRjaigKEYEEzVqjFF8jQS3/OK+RDSiKIoRNEQNJhg0LtHXiGE0ioAbAZQRZFEEBIFZnt8fXYxdvcwMMNXV3XN/rmuumXPqVPdN0TNPV3XVKXN3REREAJrFHUBERPKHioKIiNRSURARkVoqCiIiUktFQUREaqkoiIhILRUFERGppaIgIiK1VBRERKRWSdwBNlWnTp28R48ecccQESkob7311nJ371zfuIIrCj169KC8vDzuGCIiBcXMPmnIOB0+EhGRWioKIiJSS0VBRERqqSiIiEgtFQUREakVWVEws4fMbKmZzcqy3MzsbjOba2YzzWz/qLKIiEjDRHlK6njgHuDRLMuPB3oFXwcD9wXfRUQ227r1lVTX1MQdo1F9tGAp1TU1tGxRQq/u29GqZYvIniuyouDur5pZjzqGDAAe9cT9QKeb2TZmtqO7L44qk4jkj5qaGtZ8s6G2vWrNOh7/679ZsWotH87/nPZtW9G2dekmPebny1c1dsy8c9f/nk7X7TtE9vhxXrzWBViY1K4I+tKKgpmNAEYA7LzzzjkJJyKNp7KymrXrNjDt9dm8/OaHLP1ydb3rrF6zjtVr1uUgXWFp02rTCuWmirMoWIY+zzTQ3ccCYwHKysoyjhGR/LP0y9X87No/xvLcrSP+45lrJc2bccwhe7Dt1m2jfZ5IH71uFUC3pHZXYFFMWUSkEd392Ev8c8ZHDR7ftnVLANZ8sx6AU/sfQNftO9CjS0eaN9u082Hati5lq3atN2kd+VacRWEKcJGZTSLxAfNKfZ4gUpg2VFYxYcp0Vq9Zz2tvfZx1XLNmzaipqaHvPrtw7GG92Wu3nWjeXGfG55PIioKZTQSOAjqZWQXwG6AFgLv/AZgKnADMBdYCP4kqi4g0vnXrK3lrzqfMX7iMp198p86xpS1KuPWKgXTZbpscpZPNFeXZR4PrWe7AhVE9v4hEY0NlFYMvf7BBY39+5vc56qDvRJxIGlPBTZ0tIvH5bOlXXHzjpDrHnNa/jG3at+bYw3pjlul8EslnKgoiUqdv1m3ggusnsnW7Viz8fEXGMX1278r2ndpz/mlHqBAUOBUFEcno1fKPuGvCS7XtVV9/kzbm7qsH6XOCIqOiICK13J37n3yVF/71fr1j/3zXT3OQSHJNRUFEgERBGHjp/XWOufK8/uzYeWu67RDdNAsSLxUFkSZuwWfLuXPCSyxc/GXG5Qfu1YPzBh5Opw7tcpxM4qCiINJE1dTUcPO455kxa0HG5QOPO4CBx+xPixbNcxtMYqWiINIEnXnFONatr8y6/I83D4t0embJXyoKIk3Mz2+YmLUgjL7wR+y9e5ccJ5J8oqIg0oSs+WY9i5atTOs/e8AhDPjBPjEkknyjoiDShJw98uFQ+/qLB9B71x1jSiP5SEVBpAlwd0aPeTatXwVBUqkoiBS5bPMVPXHb8BjSSL7TROYiRS5TQRh0woGUlOhUU0mnPQWRIpXtCuUrz+vPQXv3yH0gKQgqCiJFatyfX0/re/L2EbrTmdRJrw6RIjS/YjnPvTYr1HfTZT9WQZB6aU9BpMiMvP0pPv5kaahv8P8cxG7dt4spkRQSFQWRInLKJX/I2D/w2P1znEQKlfYlRYqAu3PGr8ZlXKb7Hsim0J6CSIGb9fFn/Oae9AvTLhh8JP367hFDIilkKgoiBWzJF6syFoQLBx/FD/p+N4ZEUuhUFEQK1PyK5Vx+y+S0/l+cfTSHH7BbDImkGKgoiBSYqqpqrrz9aRZ8tjxtmT4/kC2loiBSYE7/5QMZ+yffeX6Ok0gx0tlHIgUk2ymnk+88HzPLcRopRtpTECkQb8/5NK3v/tFD6NShXQxppFhpT0GkAKxes44b758a6rvl8lNUEKTRqSiIFIBz/nd8qH3KMfuzS7fO8YSRoqaiIJLnvly5Jq3vjB8eFEMSaQoiLQpm1t/MPjSzuWY2MsPync3sZTP7j5nNNLMToswjUmi+XLmG4aMmhPom3npeTGmkKYisKJhZc2AMcDzQGxhsZr1Thv0aeNLd9wMGAfdGlUekEKUWBIDSFjo/RKIT5Z7CQcBcd5/n7huAScCAlDEObBX8vDWwKMI8IgVlzn8Xp/XpWgSJWpRvOboAC5PaFcDBKWNGA8+b2c+BtsDREeYRKSjX3P2XUHvirefpWgSJXJR7CplevZ7SHgyMd/euwAnABDNLy2RmI8ys3MzKly1bFkFUkfySepFaz66ddNhIciLKolABdEtqdyX98NAw4EkAd38DaAV0Sn0gdx/r7mXuXta5s07Dk+I2/d15aX23XH5KDEmkKYqyKMwAeplZTzMrJfFB8pSUMZ8C/QDMbA8SRUG7AtKk3fLQ86H2jZecpMNGkjORFQV3rwIuAqYB75M4y2i2mV1nZicGw34JDDezd4GJwDnunnqISaTJSD1s1Gf3rnx3lx1iSiNNUaQHKd19KjA1pW9U0s9zgMOizCBSCGpqajj1F2PT+kdd8D8xpJGmTFc0i8RsQ2VVxoIwcnh/HTaSnFNREInZA3/6v7S+K8/rz4F79ch9GGnydI6bSMxeevODUPuhG4aydfvWMaWRpk57CiIxmj03fJb2VSOOV0GQWKkoiMTkowVLGPX78FnaZXt2jymNSIKKgkhMrrrj6VD7yAN3jymJyLdUFERisOab9Wl9Fw/5QQxJRMJUFERybPWadZw98uFQ35/v+mlMaUTCVBREcsjd026tKZJPVBREcui0DBepTbp1eAxJRDLTdQoiOfLBvM+pSZna6093jKBZM703k/yhV6NIDqxes46r73om1HfvqDNUECTv6BUpkgOpnyO0aVXK9h23yjxYJEYqCiIRW71mXVrfuBvOjiGJSP1UFEQilOlso8duOle31pS8paIgEqFHnnkjra91q9IYkog0jIqCSESqq2t49pWZob4nbtPpp5LfVBREIjLsmkdD7SE/OpiSkuYxpRFpGBUFkYikfsB88tH7xZREpOFUFERy4Jqf6V7LUhhUFEQisGLV2lB7n+90jSmJyKZRURCJwFMvvB1qm1lMSUQ2jYqCSCNbvGwlU1+dVdvu1X27GNOIbBoVBZFG5O5cdMPEUN8xh+4RUxqRTaeiINKIrr/vb2l9/fqqKEjhUFEQaSTV1TW8+2FFqG/ynefHlEZk86goiDSS0y4L30DnyvP66wNmKTgqCiKNYPbcRWl9B+3dI/dBRLaQioJIIxj1+ymh9phrzogpiciWUVEQ2UKTnw9fk9B1+w7s0Ek30JHCpKIgsgXcnYl/+3eo7/YrT40pjciWi7QomFl/M/vQzOaa2cgsY04zszlmNtvMHo8yj0hju238P0Ltvn160ry53mtJ4Yrs9k9m1hwYAxwDVAAzzGyKu89JGtMLuAo4zN1XmJku/ZSCseab9bzxzn9DfZefe2xMaUQaR5RvaQ4C5rr7PHffAEwCBqSMGQ6McfcVAO6+NMI8Io3q7JEPh9rnnHSoTkGVghdlUegCLExqVwR9yXYHdjez181supn1z/RAZjbCzMrNrHzZsmURxRVpuFMu+UNa34++3yeGJCKNK8qikOktk6e0S4BewFHAYOBBM9smbSX3se5e5u5lnTt3bvSgIpsidW4jgN9fPSiGJCKNL8qiUAF0S2p3BVKv8KkA/uLule4+H/iQRJEQyUvuzuJlK0N9Vww7jp22S3svI1KQoiwKM4BeZtbTzEqBQcCUlDHPAN8HMLNOJA4nzYswk8gWee+jz0Ltq0Ycz8F9esaURqTxRVYU3L0KuAiYBrwPPOnus83sOjM7MRg2DfjCzOYALwO/cvcvosoksiXcnWvv/Wuor2zP7jGlEYlGZKekArj7VGBqSt+opJ8duCz4Eslb1dU1aRPeddiqTUxpRKKjq2xEGuDnN05K67t/9JAYkohES0VBpB5VVdUs+WJVqO/uqwfpymUpSnpVi9TjjXfC5z7cfuWpdNHZRlKkVBRE6nHnhBdD7e47dYwpiUj0VBRE6pA4F+JbnTu0jymJSG6oKIjU4Z7HXwm1NS22FDsVBZE6vPLvD0PtNq1LY0oikht1FgUzG5/089DI04jkkU8Wha+jPLnfvjElEcmd+vYU9kn6+ZIog4jkm8tu+lOofeaPDo4piUju1FcUUmc1FWkS/v7a7FB7t523070SpEmob5qLrmZ2N4lpsDf+XMvdL44smUiMHpj8Wqh9w8Wp94cSKU71FYVfJf1cHmUQkXyxbn1lqH3ofrvSokXzmNKI5FadRcHdH8lVEJF8ceYV40LtX55zTExJRHKv3lNSzWyomb1tZmuCr3IzOzsX4URy7fPlq+ofJFLE6txTCP74X0piauu3SXy2sD9wi5nh7o9GH1Ekdy68/vFQe9wNev8jTUt9ewoXACe7+8vuvtLdv3L3l4BTgmUiRWPR0q/S+rZpr3smSNNSX1HYyt0XpHYGfVtFEUgkLmMmvhJqP3DdWfEEEYlRfUXhm81cJlJQ3J0P5n1e2zZg263bxhdIJCb1nZK6h5nNzNBvwC4R5BGJxXnXTAi1b/rlKTElEYlXfUVhH2B7YGFKf3dgUSSJRHLsowVL+Gr12lDfrjt3jimNSLzqO3x0B7DK3T9J/gLWBstEClplZTVX3fF0qG/0hT+KKY1I/OorCj3cPe3wkbuXAz0iSSSSQ4MufyDULm1Rwt67d4kpjUj86isKrepY1roxg4jkWupd1QD+ePO5MSQRyR/1FYUZZjY8tdPMhgFvRRNJJDd+O/a5UHv8b8+hWTPdd0qatvo+aL4UeNrMzuTbIlAGlAInRxlMJGpvz/k01G7ftq4dY5Gmob4J8ZYAh5rZ94G9gu6/BVc1ixSsysrqUHvoSYfElEQkv9S3pwCAu78MvBxxFpGcGfq/40PtHx65dzxBRPKMDqBKk1NTU8P6DeF7JuizBJEE/SZIk/Pg5NdD7duuGBhTEpH8o6IgTc6018P3X+7RpVNMSUTyT6RFwcz6m9mHZjbXzEbWMW6gmbmZlUWZR2T23PDsLOefdkRMSUTyU2RFwcyaA2OA44HewGAz651hXHvgYuDNqLKIQOJitVG/nxLq+8HB34kpjUh+inJP4SBgrrvPc/cNwCRgQIZx1wM3A+sizCLC7x74e6jdplUpJSXNY0ojkp+iLApdCM+uWhH01TKz/YBu7v7XCHOI8NnSryif/Umo79Hf/SSmNCL5K8qiYBn6aiebMbNmJGZa/WW9D2Q2wszKzax82bJljRhRmoLPl6/i4hsnhfpGDu+PWaaXqEjTFmVRqAC6JbW7Er4HQ3sSV0m/YmYLgL7AlEwfNrv7WHcvc/eyzp01z700XHV1DRde/3ha/4F79ch9GJECEGVRmAH0MrOeZlYKDAJqP+Vz95Xu3snde7h7D2A6cGIwLbdIozjtsrFpfZPvPD+GJCKFIbKi4O5VwEXANOB94El3n21m15nZiVE9r8hGVVXVaX1P3DZch41E6tCguY82l7tPBaam9I3KMvaoKLNI0zPqnmdD7Vt/NVBnG4nUQ1c0S9H6cP7noXbPrrpyWaQ+KgpSlL5cuSbU/tW5x8aURKSwqChIUbp2TPjSl4P79IwpiUhhUVGQolSxZEXtz1u1a60Pl0UaSEVBis7EqTNCbR06Emk4FQUpKl+vXc/kaW+F+nrvumNMaUQKj4qCFJWhVz0cap9wxF5ZRopIJioKUjQ+XfxlWt+wUw6PIYlI4VJRkKLxi989GWr/8eZhMSURKVwqClIU3nhnXqi99+5daNWyRUxpRAqXioIUhVsffj7UHn3hj2JKIlLYVBSk4H2zbkOo3XGbtjElESl8KgpS0L5cuYYhVz4U6rvtilNjSiNS+FQUpKANHzUhra9921YxJBEpDioKUrDenDk/re+xm86NIYlI8Yj0fgoiUbp53LRQ+093jKBZM73PEdkS+g2SgvTeR5+F2n1276qCINII9FskBWfJF6sYPSZ8V7WrRvSPKY1IcVFRkIJzwXWPh9rt2rSktIWOhIo0BhUFKSip1yQAjP/tObkPIlKkVBSkoPz8xkmh9hO3DdcNdEQakYqCFJQVq9aG2iUlzWNKIlKcVBSkYKz6+ptQ+7YrBsaURKR4qShIQVi3vpKfXP1IqK/r9h1iSiNSvFQUpCCcecW4tD4dOhJpfCoKkvfumvBiWt/jt+gGOiJRUFGQvLb2mw28Wv5xqO+eXw+mZaluoCMSBRUFyWtnjQxPiz3i1O+xY+etY0ojUvxUFCRvLVr6VVrfcYfvGUMSkaZDRUHyVuqFardfeVpMSUSaDhUFyUvr1lem9XXfadsYkog0LZEWBTPrb2YfmtlcMxuZYfllZjbHzGaa2Ytm1j3KPFI4Uk9BnXzn+TElEWlaIisKZtYcGAMcD/QGBptZ75Rh/wHK3L0PMBm4Oao8Ujguv2VyWp/mNxLJjSj3FA4C5rr7PHffAEwCBiQPcPeX3X3jZDbTga4R5pECMPeTpcyvWB7qG3vtkJjSiDQ9URaFLsDCpHZF0JfNMOC5TAvMbISZlZtZ+bJlyxoxouSTxctWcuXtT4X6Tu63Lx23aRdTIpGmJ8qikGl/3zMONBsClAG3ZFru7mPdvczdyzp37tyIESWfXHTDxFC7tEUJQ07sG1MakaYpyttVVQDdktpdgUWpg8zsaOBq4Eh3Xx9hHsljz748M61v4q3nxZBEpGmLck9hBtDLzHqaWSkwCJiSPMDM9gPuB05096URZpE8tqGyivHP/CvUp7ONROIRWVFw9yrgImAa8D7wpLvPNrPrzOzEYNgtQDvgT2b2jplNyfJwUqTKZ3/C4MsfDPX1P3xPnW0kEpNI73bu7lOBqSl9o5J+PjrK55f89//Gpp9bMPzU78WQRERAVzRLTNydUy75Q1r/Pb8eHEMaEdko0j0FkWyGXfNoWt+f7/ppDElEJJn2FCTn3J2Vq8P3W77zqtNjSiMiyVQUJKc+WfQFAy+9P9R3+5Wn0m0H3W9ZJB/o8JHkzM3jpvHmzPlp/d136hhDGhHJREVBIldVVc3pv3wg47IbLzkpx2lEpC4qChK5bAVh8p3n63oEkTyjoiCRybaHMOiEAzn1uANiSCQi9VFRkEg8/Y//8Nizb6b13zvqDLbvuFUMiUSkIVQUpFG5e9rZRRtdelY/FQSRPKeiII3myb+X88Rz5RmXPXDdWWy7ddscJxKRTaWiIFvsxenvc+/Ef2ZcduL392HoSYfkOJGIbC4VBdki19z9F+b8d3HGZfePHkKnDrprmkghUVGQzfabe6ZkLAhHlPXikrP6xZBIRLaUioJssg2VVWn3QADo26cnl51zDM2ba/YUkUKloiCb5Ou16xl61cNp/T88sg8/+fGhMSQSkcakoiANdteEF3m1/OO0/pP77cuQE/vGkEhEGpuKgtTJ3XnkmTd49pWZGZePueYMduikaw9EioWKgmT1zgcLuf6+v2VdPvHW8yhtoZeQSDHRb7RkNHrMs7z30WcZl10x7DgO7tMzx4lEJBdUFCQk25lFAGed2JeT+u2b40QikksqCgIkisG8hcu5+q5n0pZ9d5cddN8DkSZCRaEJe3vOp9x4/9Q6x1x2zjEctt+uOUokInFTUWiCPlqwhDGPv0LFkhV1jhv/23No37ZVbkKJSF5QUWgiVq9Zx6vlH/PQU6/XO7a0RQmP3zJMd0UTaYJUFIrcZ0u/4uIbJ9U77sHrz6bDVm1ykEhE8pmKQpFav6GSM341rs4xHbZqwwWDj2L/3jvnKJWI5DsVhSKyfMXXnD/6sXrHHXng7gwfeDitW5XmIJWIFBIVhQJVXV3DI395g6qqGl568wMqq6rrHN+6VSnjbxxKSUnzHCUUkUKkopCnqqqqmf7ufNau28D78xZTseQr5i1cRsvSFlTX1FBVTxFIdvfVg+iy3TYRphWRYhFpUTCz/sBdQHPgQXf/XcrylsCjwAHAF8Dp7r4gykz5Zt36St54Zx7Nmhmrvl7HjFkLmD13Udbx6zdUNuhxx147hG23bqsziERkk0RWFMysOTAGOAaoAGaY2RR3n5M0bBiwwt13M7NBwE3A6VFlaixvz/mUlau/CfV9sXINnyz6gm3atw7119Q4z//rfXbt1olmzb69+cyH8z9vtDwDjzuAFiXNObhPT7bv2F6T1InIZovyr8dBwFx3nwdgZpOAAUByURgAjA5+ngzcY2bm7t7YYf7+2mxefetjampqNvsxPv5kaSzrHrLvrny5cg37frcre+62Ez27dAKgZWmJ7nImIo0qyqLQBViY1K4ADs42xt2rzGwl0BFY3phBlq/4mgcnv0ajV5pGdkRZLzZUVrN+QyVle/bgqIN2p1XLFnHHEpEmJMqikOlgdurf5YaMwcxGACMAdt5508+p/2rV2kgKwmH770aLpLN51q2vxMzYY5cd0sZWVdfQq/t2NEs6xt+ubUu227Y9LUqa69i/iOSFKItCBdAtqd0VSP0EdeOYCjMrAbYGvkx9IHcfC4wFKCsr2+S/7522bcd5Aw9nyfJVABy6BRO8NTOjR5eOOrVTRIpSlEVhBtDLzHoCnwGDgDNSxkwBhgJvAAOBl6L4PGGb9m04/nt7NfbDiogUnciKQvAZwUXANBKnpD7k7rPN7Dqg3N2nAOOACWY2l8QewqCo8oiISP0iPXfR3acCU1P6RiX9vA44NcoMIiLScDqfUUREaqkoiIhILRUFERGppaIgIiK1VBRERKSWRXBZQKTMbBnwScwxOtHIU3FEIN8z5ns+yP+M+Z4P8j9jvueDxsvY3d071zeo4IpCPjCzcncviztHXfI9Y77ng/zPmO/5IP8z5ns+yH1GHT4SEZFaKgoiIlJLRWHzjI07QAPke8Z8zwf5nzHf80H+Z8z3fJDjjPpMQUREamlPQUREaqkopDCzX5jZbDObZWYTzayVmf3RzD4M+h4ys4y3QzOzajN7J/iaksN8481sftJz75tl3aFm9nHwNTSKfHVkfC0p3yIzeybLurnYhpcE2Wab2aVB37Zm9kKwbV4wsw5Z1s3VNsyU8RYz+8DMZprZ02a2TZZ1F5jZe8E2LM9hvtFm9lnS/98JWdbtH/w+zTWzkVHkqyPjE0n5FpjZO1nWjWQbBn8/lprZrKS+jK89S7g72E4zzWz/LI95QJB1bjB+y+7Y5e76Cr5I3B50PtA6aD8JnAOcQOIucQZMBH6WZf2vY8o3HhhYz7rbAvOC7x2CnzvkKmPKmD8DZ8e0DfcCZgFtSMwS/A+gF3AzMDIYMxK4KcZtmC3jsUBJMOamTBmDZQuATjFsw9HA5fWs2xz4L7ALUAq8C/TOVcaUMbcBo3K5DYEjgP2BWUl9GV97wd+d54K/O32BN7M85r+BQ4JxzwHHb0lG7SmkKwFaW+JOcG2ARe4+1QMk/gO65lO+Bq53HPCCu3/p7iuAF4D+uc5oZu2BHwAZ9xRyYA9guruvdfcq4J/AycAA4JFgzCPASRnWzdU2zJjR3Z8P2gDTie91mG0bNsRBwFx3n+fuG4BJJLZ9TjMG76ZPI/EmL2fc/VXS7y6Z7bU3AHg0+NMzHdjGzHZMXjFob+XubwR/nx4l82u3wVQUkrj7Z8CtwKfAYmCluz+/cXlw2Ogs4O9ZHqKVmZWb2XQz26L/mM3Id2Owi3mHmbXMsHoXYGFSuyLoy2VGSPxivujuq7I8RKTbkMS7xyPMrKOZtSHxbqwbsL27Lw7+DYuB7TKsm5NtWEfGZOeSeFeYiQPPm9lblri/eS7zXRS8Dh/KcgguX7bh94Al7v5xlvWj3obJsr32GrKtugT9dY3ZJCoKSYIX8QCgJ7AT0NbMhiQNuRd41d1fy/IQO3viysMzgDvNbPNvBr1p+a4CvgscSOLQxpWZVs/Q1+innjVgGw6m7ndnkW5Dd3+fxKGXF0gU93eBqjpX+lZOtmF9Gc3s6qD9xywPcZi77w8cD1xoZkfkKN99wK7AviTeENyWYfW82IbU/zqMdBs2UEO2VaNvTxWFsKOB+e6+zN0rgaeAQwHM7DdAZ+CybCu7+6Lg+zzgFWC/XORz98XBLuZ64GESu+ipKgi/U+pKww89bXFGADPrGGT7W7aVc7ANcfdx7r6/ux9BYlf+Y2DJxl3z4PvSDKvmahtmy0jw4fYPgTODwwWZ1t24DZcCT5P59dDo+dx9ibtXu3sN8ECW582HbVgC/Bh4oo51I9+GSbK99hqyrSoIH0bc4u2pohD2KdDXzNoExxz7Ae+b2XkkjicPDl7wacysw8bDNmbWCTgMmJOjfBtfUEbieOKsDOtOA44NcnYg8aHltEbOlzVjsOxU4K+euA1rmhxtQ8xsu+D7ziT+OEwEpgAbzyYaCvwlw6q52oYZM5pZfxJ7gSe6+9os67UNPrfBzNoGGTO9HqLIl3y8++QszzsD6GVmPc2slMR92aM6yyzT/zMk3rh84O4VWdbLyTZMku21NwU4OzgLqS+JQ7GLk1cM2qvNrG/w+3Y2mV+7Dbcln1IX4xdwLfABiRfBBKAlid3O/wLvBF+jgrFlwIPBz4cC75HYTX0PGJbDfC8FzzkLeAxol5ovaJ8LzA2+fpLLbRj0vwL0TxkbxzZ8jUSxeRfoF/R1BF4k8W7yRWDbmLdhpoxzSRxj3vg6/EPQvxMwNfh5l2Cdd4HZwNU5zDch+H+bSeIP2o6p+YL2CcBHwe9UJPmyZQz6xwM/TRmbk21IojAtBipJvMsfVsdrz4AxwXZ6DyhLepx3Un6HZgXj7iG4KHlzv3RFs4iI1NLhIxERqaWiICIitVQURESkloqCiIjUUlEQEZFaKgpSFMzsakvMhjnTEjNbHhzx871iZnl9b1+RzVESdwCRLWVmh5C4ynd/d18fXPhWGnMskYKkPQUpBjsCyz0xzQfuvtyDaQrMbJSZzbDEvPpjN841H7zTv8PMXjWz983sQDMySCsFAAACh0lEQVR7yhJz2t8QjOlhifsXPBLsgUwOJlcLMbNjzewNM3vbzP5kZu0yjKn3+YJxQ8zs38Hezv1m1jzov88SEwXONrNrk8YvMLNrg+d+z8y+28jbVpoYFQUpBs8D3czsIzO718yOTFp2j7sf6O57Aa1J7FFstMET8+L8gcTUABeSmIf/nGCeJoDvAGPdvQ+wCrgg+YmDvZJfA0d7YgK1crLPj1Xn85nZHsDpJCZj2xeoBs4M1r3aExMF9gGONLM+SY+7PHju+4DL699cItmpKEjBc/evgQOAEcAy4AkzOydY/H0ze9PM3iNxH4c9k1bdOOfOe8BsT0wsuJ7EzXM2TkS20N1fD35+DDg85en7Ar2B1y1xF6+hQPcsUet7vn7Bv2NG8Fj9SEy5AHCamb0N/Cf4N/ROetyngu9vAT2yPLdIg+gzBSkK7l5NYm6lV4ICMNTMJpGY7rzM3Rea2WigVdJq64PvNUk/b2xv/N1InQcm09TFL7j74AbErO/5DHjE3a8KPYFZTxJ7AAe6+wozG5/l31GNfqdlC2lPQQqemX3HzHolde0LfMK3fziXB8f5B27Gw+8cfJANiTn4/y9l+XTgMDPbLcjSxsx234zngcRkaAOTZvfc1sy6A1sBa4CVZrY9iTn+RSKhdxVSDNoBv7fEjeyrSMwmOsLdvzKzB0gcrllAYtrmTfU+ib2O+0nMYnlf8kJ3XxYcqppo397x7tckZgHdJO4+x8x+TeKOX81IzKR5obtPN7P/kJixcx7wel2PI7IlNEuqSBZm1oPE/R/2ijmKSM7o8JGIiNTSnoKIiNTSnoKIiNRSURARkVoqCiIiUktFQUREaqkoiIhILRUFERGp9f8BQvXyf8SXzfcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cdf = thinkstats2.Cdf(xbars)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.Config(xlabel='Sample mean',\n",
    "                 ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of the sample means is close to the actual value of $\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "89.9861531227758"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(xbars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interval that contains 90% of the values in the sampling disrtribution is called a 90% confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(85.67721360886706, 94.03550577193234)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci = cdf.Percentile(5), cdf.Percentile(95)\n",
    "ci"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the RMSE of the sample means is called the standard error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4985255047938177"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stderr = RMSE(xbars, 90)\n",
    "stderr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confidence intervals and standard errors quantify the variability in the estimate due to random sampling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating rates\n",
    "\n",
    "The following function simulates experiments where we try to estimate the mean of an exponential distribution using the mean and median of a sample. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse L 1.0991868744984798\n",
      "rmse Lm 1.4874648694733599\n",
      "mean error L 0.32089631663544593\n",
      "mean error Lm 0.35349315821339045\n"
     ]
    }
   ],
   "source": [
    "def Estimate3(n=7, iters=1000):\n",
    "    lam = 2\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(iters):\n",
    "        xs = np.random.exponential(1.0/lam, n)\n",
    "        L = 1 / np.mean(xs)\n",
    "        Lm = np.log(2) / thinkstats2.Median(xs)\n",
    "        means.append(L)\n",
    "        medians.append(Lm)\n",
    "\n",
    "    print('rmse L', RMSE(means, lam))\n",
    "    print('rmse Lm', RMSE(medians, lam))\n",
    "    print('mean error L', MeanError(means, lam))\n",
    "    print('mean error Lm', MeanError(medians, lam))\n",
    "    \n",
    "Estimate3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RMSE is smaller for the sample mean than for the sample median.\n",
    "\n",
    "But neither estimator is unbiased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  In this chapter we used $\\bar{x}$ and median to estimate µ, and found that $\\bar{x}$ yields lower MSE. Also, we used $S^2$ and $S_{n-1}^2$ to estimate σ, and found that $S^2$ is biased and $S_{n-1}^2$ unbiased.\n",
    "Run similar experiments to see if $\\bar{x}$ and median are biased estimates of µ. Also check whether $S^2$ or $S_{n-1}^2$ yields a lower MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 1\n",
      "mean error xbar 0.0005094416922616375\n",
      "mean error median 0.0008580279056516486\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def Estimate4(n=7, iters=100000):\n",
    "    \"\"\"Mean error for xbar and median as estimators of population mean.\n",
    "\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        xbar = np.mean(xs)\n",
    "        median = np.median(xs)\n",
    "        means.append(xbar)\n",
    "        medians.append(median)\n",
    "\n",
    "    print('Experiment 1')\n",
    "    print('mean error xbar', MeanError(means, mu))\n",
    "    print('mean error median', MeanError(medians, mu))\n",
    "    \n",
    "Estimate4()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 2\n",
      "RMSE biased 0.5135938607123939\n",
      "RMSE unbiased 0.5763054869567986\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def Estimate5(n=7, iters=100000):\n",
    "    \"\"\"RMSE for biased and unbiased estimators of population variance.\n",
    "\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    estimates1 = []\n",
    "    estimates2 = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        biased = np.var(xs)\n",
    "        unbiased = np.var(xs, ddof=1)\n",
    "        estimates1.append(biased)\n",
    "        estimates2.append(unbiased)\n",
    "\n",
    "    print('Experiment 2')\n",
    "    print('RMSE biased', RMSE(estimates1, sigma**2))\n",
    "    print('RMSE unbiased', RMSE(estimates2, sigma**2))\n",
    "\n",
    "Estimate5()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) xbar and median yield lower mean error as m increases, so neither\n",
    "# one is obviously biased, as far as we can tell from the experiment.\n",
    "\n",
    "# 2) The biased estimator of variance yields lower RMSE than the unbiased\n",
    "# estimator, by about 10%.  And the difference holds up as m increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose you draw a sample with size n=10 from an exponential distribution with λ=2. Simulate this experiment 1000 times and plot the sampling distribution of the estimate L. Compute the standard error of the estimate and the 90% confidence interval.\n",
    "\n",
    "Repeat the experiment with a few different values of `n` and make a plot of standard error versus `n`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standard error 0.825791904537049\n",
      "confidence interval (1.311259184670713, 3.7054533946140884)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.825791904537049"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcVXX9x/HXm2EQREAR3FhVcCHDbcJcSvsVpeaSSalppoloZblkLi2m9ktbzOxXWprmloGGP/uhYmqmlaYGkhuYRog6IJtssi/z+f1xL+OcO/swZ869c9/Px2Me3u/3fufc952R+dzzPed8jyICMzMzgC5ZBzAzs+LhomBmZrVcFMzMrJaLgpmZ1XJRMDOzWi4KZmZWy0XBOiVJl0v6bf7xYEkrJFV0wOueJunJOu0VknZpp21/U9LN+cdDJYWkru207Q77GVlxc1GwdiXpEEl/l7RM0mJJT0n6QJaZIuLNiNgqIjZm8NpbRcSspsZIOkxSdQu2dVVEjG2PXJJmS/pYnW1n9jOy4tIunzLMACT1Bh4AvgTcA3QDPgSszTJXZyCpa0RsyDqHdX7eU7D2tBtARIyPiI0RsToiHomIFwEk7Srpz5LekbRI0l2Stt70zflPr9+Q9KKklZJukbS9pIckvSvpT5K2yY/dNH0yTtJcSW9L+npDoQqnWiQ9Iel7+b2YdyU9IqlfnfGnSnojn/M7hZ+qC7a9raRJkpZL+gewa8HzIWlY/vGRkmbkX3OOpAsl9QQeAnbKT9+skLRTfvproqTfSloOnFZ3SqyOLzb0/iXdJum/67Rr90Yk3QkMBu7Pv95FDfyMdsq/r8WSZko6s862Lpd0j6Q78u9luqSqBv+PsJLjomDt6TVgo6TbJR2x6Q94HQKuBnYC9gQGAZcXjDkeGE2uwBxN7g/mN4F+5P5//VrB+I8Aw4GPA5c09se7AZ8DTge2I7dHcyGApBHADcDJwI5AH2BAE9u5HliTH/vF/FdjbgHOiohewF7AnyNiJXAEMDc/fbNVRMzNjz8WmAhsDdzVyDZb/f4j4vPAm8DR+df7UQPDxgPV5H5XY4CrJH20zvPHABPy2SYBv2juda00uChYu4mI5cAhQAC/BhbmP21un39+ZkQ8GhFrI2IhcC1waMFmfh4R8yNiDvA34NmI+GdErAXuA/YtGH9FRKyMiJeAW4GTWhj31oh4LSJWk5vq2iffPwa4PyKejIh1wGX591NP/qDs8cBl+QwvA7c38ZrrgRGSekfEkoiY1kzGpyPiDxFRk8/ZkLa+/0ZJGkTu93hxRKyJiOeBm4HP1xn2ZERMzh+DuBPYe3Nf14qDi4K1q4h4JSJOi4iB5D4N7wRcByBpO0kT8lMny4HfktsDqGt+ncerG2hvVTD+rTqP38i/XkvMq/N4VZ3t7lR3mxGxCninkW30J3dcrjBDY44HjgTekPQXSQc2k/GtZp4vHNOa99+UnYDFEfFuwbbr7jEV/vy6t9eZUJYtFwVLTUT8C7iNXHGA3NRRACMjojdwCrkppc0xqM7jwcDcxga20NvAwE0NST2AbRsZuxDY0ECGBkXElIg4ltyU1R/I7aFAI3siTfTX1dj7XwlsWee5HVqx7blAX0m9CrY9pwV5rMS5KFi7kbSHpK9LGphvDyI3nfFMfkgvYAWwVNIA4Bvt8LLfkbSlpPeRO0Zw92ZubyJwtKSDJHUDrqCRwpWfOvlf4PJ8hhHAFxoaK6mbpJMl9YmI9cByYNPpn/OBbSX1aUPext7/88CRkvpK2gE4r+D75gMNXj8REW8BfweultRd0kjgDBo/rmGdiIuCtad3gQOAZyWtJFcMXgY2nRVzBbAfsAx4kNwf1M31F2Am8BhwTUQ8sjkbi4jpwFfJHUR9m9x7WkDjp9WeQ27qaR65vaJbm9j854HZ+amzs8ntKW3aoxoPzJK0VFJrpoAae/93Ai8As4FHqF8srwa+nX+9CxvY7knAUHJ7DfcB342IR1uRy0qUfJMdK0WShgKvA5Vpnr8vaStgKTA8Il5P63XMioX3FMwKSDo6PyXTE7gGeIncJ26zTs9Fway+Y8lNm8wldw3AieFdaisTnj4yM7Na3lMwM7NaJXexSb9+/WLo0KFZxzAzKynPPffcoojo39y4kisKQ4cOZerUqVnHMDMrKZKautq+lqePzMyslouCmZnVclEwM7NaLgpmZlbLRcHMzGqlVhQk/UbSAkkvN/K8JP1P/lZ/L0raL60sZmbWMmmeknobuVv03dHI80eQW0JgOLmVNX+Z/6+ZWacTEaxas65N3/va7AVsrKlhi8quDB+yHd23qGzndO9JrShExF/zK1k25ljgjvyaMs9I2lrSjhHxdlqZ0jRz5szax8OGDcswibWHzvz73LBhI6vXrk/9dRa88y4rVidXHI8IXps9n64VFam85ptvL+adpSvo02vL5gd3oKef/0+7betn3zyBgdsX3v68/WR58doAkrcSrM731SsKksYB4wAGD270xlZmZa+mpoZ/vDSb1+e8Q2XXCmbMnEsEbNm9kmde9MrfncGW3buluv0si0JDd7NqcHW+iLgJuAmgqqrKK/hZWYsI/jVrHr+c8BfmLFjKtlv3BOCdpSszTmYt1aMNf9i7VnRh9IF70rdPzxQS1XmdVLfetGqS95cdyObfX9esU1i/fiMTH3mOPz/7KjsP6Ffb//LMuaxdl5z6aWsxqOxaQbfK9P4EbKypYc3a9bx/twGo4DPg7Lnv8MG9d07lU+/qNevpu3VPduzflrubpqdmYw377DmIXj27Zx2lSVkWhUnAOZImkDvAvKxUjyeYtcbqNet4d9Va3nx7MRs2bKztf/HVOTz81PR64xcva9sf/eFDtuP9wwewYvVahuy4LX169SAIdh7Qjx369UZq8NbTVuZSKwqSxgOHAf0kVQPfBSoBIuJXwGTgSHL3l11F7qbjZp3W+vUbOfHCX7frNntv1YMTj6hivxHvHWvr26cnFRW+BMnaJs2zj05q5vkAvpLW65tlZcHid/neDQ9QWdmV3lt156XX5mzW9j44cmcOHbU7mz7YR8DuQ7enT68e7ZDWLKnkls42KyZvzVvC7x9+jhkz57Jk+apWfe9WW27BilVrGfX+obV9CxavYOiAbTnmIyMZvGNfT/FYh3NRMGuhmpoarh//F56cNpNhg/vzr1nz2rytX373ZLbr26sd05m1DxcFswY89NSrLFq6kqnT5zByz6G8Nnt+4vmWFoQvfvpgBu6wDV0rujB4x75s2b2b5/utqLkomOWtWbueky+6BYCVK98746ewIDTms4dXse+egxiyU1+26JbeMgRmaXJRsLK3cvVaTr3k1lZ9z/Gj92Pk7gPYsns3dh7Yz3P/1mm4KFjZWrRkBbfe93eeeWFWo2MO3X8XRo4Yxi4D+1FR0YVdBvaja9d01u0xKwYuClZ2nvrnf7j2tkebHHPx6YcxaPvcFbGdbUE8s6a4KFhZmLdoOU8//x9+e/+zTY67+Xunsk3vLROrpJqVExcF69R+98A/uPfRac2O+8KnDmT0gXu2aaEys87ERcE6pRderebKGx5ocszhh7yPsWMO8UFiszpcFKxTaW59ocquFYwauTOnferA1JcgNitFLgpW8uYsWMrv7n+2yZvI/NcBe/Dlkw71XoFZM1wUrKTddM/fGlxuepMd+vXmqvOO8+JxZi3komAlZ8Wqtdz90BQm//XlJsddd+kJDNohvXvZmnVGLgpWUm6598kmi8HnjhrF0YeNTPWOYmadmf/lWMn4yW2P8vd//qfB5yq7VjDhJ2d2cCKzzsdFwYrehg0bOeHrDZ9RdOSH92LMx/f3MQOzduKiYEWtscXqPn7wCM767IczSGTWubkoWNFavGwlZ152Z73+k486gE+P3jeDRGadn4uCFZ0ly1cx9jt3NPjcHT84nZ49tujgRGblw0XBisaDf3mJ3/zvU40+P/G6s3zxmVnKXBSsKPzhsee5c9IzDT4n4M4fftEFwawDuChY5u770z8bXdL6p5d8lsE79u3gRGbly0XBMvWru//Co39/JdG3/4ghXDrucO8ZmGXARcEy0di1Bzv06803zzoig0RmBtAl6wBWnhq7GO3673yug5OYWV3eU7AO95nzbqzXd+KRH+Azn9g/gzRmVpeLgnWoc6+6m5qIRN9tV51Gr57dM0pkZnW5KFiHiAiuu/MxqucvSfRfcNpoFwSzIuKiYKmLCMY0MGV09gkf5uB9d80gkZk1xgeaLXUNFYQRu+7I6INGZJDGzJqSalGQdLikVyXNlHRJA88PlvS4pH9KelHSkWnmsY730mtz6vWd9dkP872vHZtBGjNrTmrTR5IqgOuB0UA1MEXSpIiYUWfYt4F7IuKXkkYAk4GhaWWyjlVTU8Pl19+f6Lv8K0fz/t0GZJTIzJqT5p7CKGBmRMyKiHXABKDw42EAvfOP+wBzU8xjHeza2x9LtLfpvaULglmRS7MoDADeqtOuzvfVdTlwiqRqcnsJX21oQ5LGSZoqaerChQvTyGrtbN36DTz9fPLWmb++8vMZpTGzlkqzKDS0cE0UtE8CbouIgcCRwJ2S6mWKiJsioioiqvr3759CVGtP69Zv4KQLb070XTruCK9lZFYC0iwK1cCgOu2B1J8eOgO4ByAinga6A/1SzGQpW7N2fb2CAFD1viEZpDGz1kqzKEwBhkvaWVI34ERgUsGYN4GPAkjak1xR8PxQiVqzdj0nX3RLvf7brjqt48OYWZukdvZRRGyQdA7wMFAB/CYipku6EpgaEZOArwO/lnQ+uaml0yKicIrJSsDU6W9w9U0P1esff81YulX6GkmzUpHqv9aImEzuAHLdvsvqPJ4BHJxmBkvfGd++g6XvrqrXP+GaM6msrMggkZm1la9ots2yavW6BgvC+GvGuiCYlSDv19tm+c7Pk4eJqt43hEvO9F3TzEqVi4K12d+m/pvZcxYl+i4d57ummZUyTx9Zm7w1bwnX3Zm8YvmHF3w6ozRm1l5cFKxNCs806lbZlWFDtssojZm1FxcFa7W3Fy5j/jvLE33jrxmbURoza08+pmCtctWND/HcjDcSfdde/JmM0phZe/OegrXY69WL6hUEgCE7bZtBGjNLg/cUrEXWrlvPhT+emOh737CduHjsJzJKZGZpcFGwFvnBrx+u13flV4/JIImZpcnTR9YiL89MLnB778/OziiJmaXJRcGaNW/RcmpqamrbF489PMM0ZpYmFwVr1i8nPJFo77nLDtkEMbPUuShYkyKCl/+dnDrq1bN7RmnMLG0uCtakUy+5NdH+9tmfzCiJmXUEFwVr1OJlK1m1Zl2ib989BzUy2sw6AxcFa9SZl92ZaH/3y0dllMTMOoqLgjXohVer6/WN3H1gBknMrCO5KFg9NTU1XHnDA4k+L3hnVh5cFKyez5x/U72+bpW++N2sHLgoWMI9f5xar89XL5uVDxcFS7j7oWRRuOtHZ2SUxMyy4KJgtZ5+flaiPfqgPem+RWVGacwsCy4KBsC69Ru45tZHEn1nn3BoRmnMLCsuCgbANb95NNEe9f6h2QQxs0y5KBhAvTuqeSVUs/LkomBUz1+SaJ/0yVEZJTGzrLkoGOdedXeifexH9s4oiZllzUWhzK0uWPAOoLKyIoMkZlYMXBTK3DMvvJ5o3/2TMzNKYmbFINWiIOlwSa9KminpkkbGfFbSDEnTJf0uzTyWtGDxu/zid48n+rp29V6CWTlLbUEbSRXA9cBooBqYImlSRMyoM2Y4cClwcEQskbRdWnmsvi9dcVeifcDInTNKYmbFIs09hVHAzIiYFRHrgAnAsQVjzgSuj4glABGxIMU8VscT/3i1Xt+Fp4/OIImZFZM0i8IA4K067ep8X127AbtJekrSM5IaPDle0jhJUyVNXbhwYUpxy8vP70pOG0287iy6dPEhJrNyl+ZfATXQFwXtrsBw4DDgJOBmSVvX+6aImyKiKiKq+vfv3+5By01E8tfwgb2GIjX06zKzcpNmUagG6t7QdyAwt4Ex/xcR6yPideBVckXCUnTBD3+faF889hMZJTGzYpNmUZgCDJe0s6RuwInApIIxfwA+AiCpH7nppFlYaiKCN99enOjzXoKZbZJaUYiIDcA5wMPAK8A9ETFd0pWSjskPexh4R9IM4HHgGxHxTlqZDH73wD8S7YvO8F6Cmb0n1XssRsRkYHJB32V1HgdwQf7LOsAjf5+RaPs0VDOry6eblJE1a9ezYtXa2vbHDx6RYRozK0YuCmXk5ItuSbSP+9i+GSUxs2LlolAm/v3G/Hp92/XtlUESMytmLgpl4pJr70u0b7z8lIySmFkxc1EoA4uWrEi0e/bYgn7bbJVRGjMrZi4KZeCsy3+baN921RcySmJmxc5FoZObPrPwInK8xpGZNarJvw6Sbqvz2B8vS9BlP09eRD7+mrEZJTGzUtDcR8a6N+s9N80glr5dB/WnW2Wq1yuaWYlrrigUrmpqJaTwAPOXTjw0oyRmViqa+9g4UNL/kFsGe9PjWhHxtdSS2Wa78Z6/Jto7D+yXURIzKxXNFYVv1Hk8Nc0g1v6mzXiz9nGvnt0zTGJmpaLJohARt3dUEGtfby9clmh/4dgDM0piZqWk2XMTJX1B0jRJK/NfUyWd2hHhrO3++1cPJtqHjdotoyRmVkqa3FPI//E/j9zS1tPIHVvYD/ixJCLijvQjWmtVz1/CvEXLa9u9t+rhG+mYWYs0t6fwZeC4iHg8IpZFxNKI+DNwfP45K0I33fO3RPv8Uz+aURIzKzXNFYXeETG7sDPf1zuNQLb5Cq9iHrn7wIySmFmpaa4orG7jc5aR1WvWJdoXjz08oyRmVoqaOyV1T0kvNtAvYJcU8thm+sHNf0y0P7DXkIySmFkpaq4o7A1sD7xV0D8EqL/SmmXu5X8nfy0+wGxmrdHc9NFPgeUR8UbdL2BV/jkrIlNenp1of+WkwzLJYWalq7miMDQi6k0fRcRUYGgqiaxNIoIf/Do5dfRfH9wjozRmVqqaKwpNrY3Qoz2D2Ob56R2PJdp77LJDRknMrJQ1VxSmSDqzsFPSGcBz6USytnhq2sxE+8pzjskoiZmVsuYONJ8H3CfpZN4rAlVAN+C4NINZy23YsDHRvuzLR1FR4burmVnrNbcg3nzgIEkfAfbKdz+Yv6rZisT0/7ydaI/cbUBGScys1LXoNlwR8TjweMpZrI3uuv/ZRNunoZpZW3mOoRP4z1sLax/vvrMPMJtZ27kolLi35i1JtMd8fL+MkphZZ+CiUMIigvOuvjvRt88eXvzOzNrORaGEXX79/Yl2ZdcKunTxr9TM2i7VvyCSDpf0qqSZki5pYtwYSSGpKs08nUlNTU29dY7GXzM2ozRm1lmkVhQkVQDXA0cAI4CTJI1oYFwv4GvAs4XPWePOvuKuRPv04w7yWUdmttnS3FMYBcyMiFkRsQ6YABzbwLjvAT8C1qSYpdN5Z+nKRPuow0ZmlMTMOpM0i8IAkktuV+f7aknaFxgUEQ80tSFJ4yRNlTR14cKFTQ0tC3Xvvwxw6bgjMkpiZp1NmkWhobmMqH1S6kJu+e2vN7ehiLgpIqoioqp///7tGLE0Pf6PVxPtqvf5Rjpm1j7SLArVwKA67YEkb8zTi9zSGU9Img18EJjkg83Nm/jwe2sR9u3TM8MkZtbZpFkUpgDDJe0sqRtwIjBp05MRsSwi+kXE0IgYCjwDHJO/V4M1YuHidxPtzx6+f0ZJzKwzSq0oRMQG4BzgYeAV4J6ImC7pSkle17mNCs86Gn1QvRO6zMzarEUL4rVVREwGJhf0XdbI2MPSzNIZFC6Rve3Wnjoys/bly19LSOE6Rz++cExGScyss3JRKCEX/nhiot2nl++Iambty0WhRFx5Q/JSjv1H+DRUM2t/Lgol4oVXqxPtU445IKMkZtaZuSiUgGkz3ky0x3xifwbv2DejNGbWmbkolIDv35g4gYsTj/D1fWaWDheFEuTVUM0sLS4KRW5BwRXME645M6MkZlYOXBSK3Ld/9odEu7KyIqMkZlYOXBSKWEQk7pswaIdtMkxjZuXARaGIjX9wSqJ90RmfyCiJmZULF4Uidu+j0xLtnbbbOqMkZlYuXBSK1Lr1GxLtcz73kYySmFk5cVEoUg888VKi/ZEDds8oiZmVExeFInXXA89mHcHMypCLQhFavGxlon3UoSMzSmJm5cZFoQidd/U9ifZpxx2YURIzKzcuCkUmIli5em2iz8tamFlHcVEoMmdedmeiffX5x2WUxMzKkYtCkVmyfFWivdvQ7TNKYmblyEWhiEREon3JmYdnlMTMypWLQhGZu3BZor337gMzSmJm5cpFoYhce9ufEu1ulV0zSmJm5cpFoYjMnrMo6whmVuZcFIrEhg0bE+3Lv3J0RknMrJy5KBSJc6++O9Hea/hOGSUxs3LmolAk5i1anmj7gjUzy4KLQhEaO+aQrCOYWZlyUSgCT06bmWh/cO9dMkpiZuXORaEI/PT25Kmo2/TeMqMkZlbuUi0Kkg6X9KqkmZIuaeD5CyTNkPSipMckDUkzTzGqqalJtPfYZYeMkpiZpVgUJFUA1wNHACOAkySNKBj2T6AqIkYCE4EfpZWnWF1xwwOJ9kVf/ERGSczM0t1TGAXMjIhZEbEOmAAcW3dARDweEZtWgHsGKLt1HV7+99xEu0+vHhklMTNLtygMAN6q067O9zXmDOChhp6QNE7SVElTFy5c2I4Rs/XWvCWJ9g8v+HRGSczMctIsCg2daB8N9CHpFKAK+HFDz0fETRFRFRFV/fv3b8eI2Tqv4IK1YUO2yyiJmVlOmiuuVQOD6rQHAnMLB0n6GPAt4NCIWFv4fGdVuKyFmVkxSHNPYQowXNLOkroBJwKT6g6QtC9wI3BMRCxIMUvRefG1OYn2HT84PaMkZmbvSa0oRMQG4BzgYeAV4J6ImC7pSknH5If9GNgK+L2k5yVNamRznc73b5ycaPfssUVGSczM3pPqgv0RMRmYXNB3WZ3HH0vz9YtV4R3Wdh3UeY6TmFlp8xXNGXi9OnnfhG+ddWRGSczMklwUMvCnp/+VaPvaBDMrFi4KGXj4qem1j/tv0yvDJGZmSS4KHazwVNRTjj4goyRmZvW5KHSws6+4K9E+ZP9hGSUxM6vPRaEDbdiwkSXLVzU/0MwsIy4KHWjS4y8m2leff1xGSczMGuai0IHueuDZRHu3odtnlMTMrGEuCh2k8IK140fvl1ESM7PGuSh0kAkPTU20T/rkBzJKYmbWOBeFDjB3wVImPvxcok9qaGVxM7NsuSh0gK9+f0KiPXbMIRklMTNrmotCyo4/91f1+o740F4ZJDEza56LQormv7O8Xt/vfzougyRmZi3jopCiR5+akWiffcKH6dLFP3IzK17+C5WSiOC+x55P9I0+aERGaczMWsZFISVnfOeORPvko7zwnZkVPxeFlCx7d3WifdzH9skoiZlZy7kopODa2/+UaH//3E/5ugQzKwkuCu1s7boNPDVtZqJvj112yCiNmVnruCi0o5Wr1/H1ax9M9J33+Y9mlMbMrPVcFNrJ+vUbufhnD9Xr/1DV8AzSmJm1jYtCOzn/Jw/U6xt/zdgMkpiZtV3XrAN0Bs9Nf6Ne38TrzvLBZTMrOS4KmyEiOO/qe6ievyTR/9sfftEFwcxKkovCZhhz3o31+kYO35Ee3btlkMbMbPP5mEIbNbT6KcC440d1cBIzs/bjPYU2mPzXl+r1ffWkg9h9SP8M0piZtR8XhVZYsWot3/zpfcxZsDTR//1zP0XXmhUZpTIzaz8uCi3w7Iuvc8f/Pc28RfXvj3DwfsPYY5cdmDlzZgPfaWZWWlwUmvCvWfP41s/+0Ojzuw7qz/mn+oplM+s8Ui0Kkg4HfgZUADdHxA8Knt8CuAPYH3gHOCEiZqeZqTHL3l3NHx57nvUbNvLI32ewcWNNk+Mv+/JR7L37wA5KZ2bWMVIrCpIqgOuB0UA1MEXSpIioezuyM4AlETFM0onAD4ET0sjz3V9M4uV/z23wzmc1NU0XgE2GDujHmWMO8QJ3ZtZppbmnMAqYGRGzACRNAI4F6haFY4HL848nAr+QpIiI9gyyaMkKXv73XKDlBaDQhGvOpLKyoj1jmZkVnTSLwgDgrTrtaqDw9mO1YyJig6RlwLbAorqDJI0DxgEMHjy41UGWLl/VqvGnHnsgAHvusgPDh2znq5PNrGykWRQa+ktauAfQkjFExE3ATQBVVVWt3ovo13crxo45hHeWrGDwTn05aJ9dGxzXtav3BMysvKVZFKqBQXXaA4G5jYypltQV6AMsbu8gW/fakiM+tFd7bzZh2LBhqW7fOpZ/n1au0lzmYgowXNLOkroBJwKTCsZMAr6QfzwG+HN7H08wM7OWS21PIX+M4BzgYXKnpP4mIqZLuhKYGhGTgFuAOyXNJLeHcGJaeczMrHmpXqcQEZOByQV9l9V5vAb4TJoZzMys5bxKqpmZ1XJRMDOzWi4KZmZWy0XBzMxqqdTOAJW0EHgjxZfoR8EV1SWo1N+D82ev1N+D89c3JCKavRNYyRWFtEmaGhFVWefYHKX+Hpw/e6X+Hpy/7Tx9ZGZmtVwUzMyslotCfTdlHaAdlPp7cP7slfp7cP428jEFMzOr5T0FMzOr5aJgZma1XBTyJP1G0gJJL2edpS0kDZL0uKRXJE2XdG7WmVpLUndJ/5D0Qv49XJF1praQVCHpn5IeyDpLa0maLeklSc9Lmpp1ntaStLWkiZL+lf+3cGDWmVpD0u75n/2mr+WSzuvQDD6mkCPpw8AK4I6ISPeOPCmQtCOwY0RMk9QLeA74VETMaOZbi4Zy9z3tGRErJFUCTwLnRsQzGUdrFUkXAFVA74g4Kus8rSFpNlAVESV54Zek24G/RcTN+fu4bBkRS7PO1RaSKoA5wAERkeYFuwneU8iLiL+Swl3fOkpEvB0R0/KP3wVeIXcP7JIROSvyzcr8V0l9apE0EPgkcHPWWcqNpN7Ah8ndp4WIWFeqBSHvo8B/OrIggItCpyRpKLAv8Gy2SVovP/XyPLAAeDQiSu09XAdcBNRkHaSNAnhE0nOSxmUdppV2ARYCt+an726W1DPrUJuubDi6AAADq0lEQVThRGB8R7+oi0InI2kr4F7gvIhYnnWe1oqIjRGxD7l7eo+SVDJTeZKOAhZExHNZZ9kMB0fEfsARwFfy06qloiuwH/DLiNgXWAlckm2ktslPfR0D/L6jX9tFoRPJz8PfC9wVEf+bdZ7Nkd/tfwI4POMorXEwcEx+Xn4C8F+SfpttpNaJiLn5/y4A7gNGZZuoVaqB6jp7lxPJFYlSdAQwLSLmd/QLuyh0EvmDtLcAr0TEtVnnaQtJ/SVtnX/cA/gY8K9sU7VcRFwaEQMjYii5Xf8/R8QpGcdqMUk98ycpkJ92+ThQMmfjRcQ84C1Ju+e7PgqUzIkWBU4ig6kjSPkezaVE0njgMKCfpGrguxFxS7apWuVg4PPAS/k5eYBv5u+TXSp2BG7Pn3XRBbgnIkrutM4Stj1wX+7zBV2B30XEH7ON1GpfBe7KT7/MAk7POE+rSdoSGA2clcnr+5RUMzPbxNNHZmZWy0XBzMxquSiYmVktFwUzM6vlomBmZrVcFMyaIek0STvVad8saUQ7bHeopM9t7nbM2pOLglnzTgNqi0JEjG2n1WeHAi4KVlRcFKxsSTolf/+G5yXdmF+M7zZJL+fvKXC+pDHklsG+Kz+uh6QnJFXlt7FC0g/zC8j9SdKo/POzJB2THzNU0t8kTct/HZSP8APgQ/ntnp9//R9LmiLpRUmZXLxk5c0Xr1lZkrQn8CPg0xGxXtINwHzgkIgYnR+zdUQslfQEcGFETM3317YlBXBkRDwk6T6gJ7mls0cAt0fEPvkrVGsiYo2k4cD4iKiSdFh+O0fltzsO2C4i/lvSFsBTwGci4vWO+rmYeZkLK1cfBfYHpuSXdegB/BHYRdLPgQeBR1qwnXX57wN4CVibLzIvkZsegtx9IX4haR9gI7BbI9v6ODAyv3cC0AcYDrgoWIdxUbByJXKf5C9NdErfAj4BfAX4LPDFZrazPt7b3a4B1gJERI2kTf++zie3F7I3uSnbNU1k+mpEPNzK92LWbnxMwcrVY8AYSdsBSOoraQjQJSLuBb7De8suvwv02ozX6gO8HRE15BYtrGhkuw8DX8ovgY6k3Ur8JjFWgrynYGUpImZI+ja5u4x1AdYDF5BbJXTTh6VNexG3Ab+StBpoy43gbwDulfQZ4HFyN38BeBHYIOmF/Gv8jNyU07T8UugLgU+14fXM2swHms3MrJanj8zMrJaLgpmZ1XJRMDOzWi4KZmZWy0XBzMxquSiYmVktFwUzM6v1/4PlTLsEFAc8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def SimulateSample(lam=2, n=10, iters=1000):\n",
    "    \"\"\"Sampling distribution of L as an estimator of exponential parameter.\n",
    "\n",
    "    lam: parameter of an exponential distribution\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    def VertLine(x, y=1):\n",
    "        thinkplot.Plot([x, x], [0, y], color='0.8', linewidth=3)\n",
    "\n",
    "    estimates = []\n",
    "    for _ in range(iters):\n",
    "        xs = np.random.exponential(1.0/lam, n)\n",
    "        lamhat = 1.0 / np.mean(xs)\n",
    "        estimates.append(lamhat)\n",
    "\n",
    "    stderr = RMSE(estimates, lam)\n",
    "    print('standard error', stderr)\n",
    "\n",
    "    cdf = thinkstats2.Cdf(estimates)\n",
    "    ci = cdf.Percentile(5), cdf.Percentile(95)\n",
    "    print('confidence interval', ci)\n",
    "    VertLine(ci[0])\n",
    "    VertLine(ci[1])\n",
    "\n",
    "    # plot the CDF\n",
    "    thinkplot.Cdf(cdf)\n",
    "    thinkplot.Config(xlabel='estimate',\n",
    "                     ylabel='CDF',\n",
    "                     title='Sampling distribution')\n",
    "\n",
    "    return stderr\n",
    "\n",
    "SimulateSample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) With sample size 10:\n",
    "\n",
    "# standard error 0.762510819389\n",
    "# confidence interval (1.2674054394352277, 3.5377353792673705)\n",
    "\n",
    "# 2) As sample size increases, standard error and the width of\n",
    "#    the CI decrease:\n",
    "\n",
    "# 10      0.90    (1.3, 3.9)\n",
    "# 100     0.21    (1.7, 2.4)\n",
    "# 1000    0.06    (1.9, 2.1)\n",
    "\n",
    "# All three confidence intervals contain the actual value, 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In games like hockey and soccer, the time between goals is roughly exponential. So you could estimate a team’s goal-scoring rate by observing the number of goals they score in a game. This estimation process is a little different from sampling the time between goals, so let’s see how it works.\n",
    "\n",
    "Write a function that takes a goal-scoring rate, `lam`, in goals per game, and simulates a game by generating the time between goals until the total time exceeds 1 game, then returns the number of goals scored.\n",
    "\n",
    "Write another function that simulates many games, stores the estimates of `lam`, then computes their mean error and RMSE.\n",
    "\n",
    "Is this way of making an estimate biased?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SimulateGame(lam):\n",
    "    \"\"\"Simulates a game and returns the estimated goal-scoring rate.\n",
    "\n",
    "    lam: actual goal scoring rate in goals per game\n",
    "    \"\"\"\n",
    "    goals = 0\n",
    "    t = 0\n",
    "    while True:\n",
    "        time_between_goals = random.expovariate(lam)\n",
    "        t += time_between_goals\n",
    "        if t > 1:\n",
    "            break\n",
    "        goals += 1\n",
    "\n",
    "    # estimated goal-scoring rate is the actual number of goals scored\n",
    "    L = goals\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 4\n",
      "rmse L 1.4126517617587144\n",
      "mean error L -0.001049\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFDNJREFUeJzt3X+wX3V95/Hny7CA1tai3P7YhBioQQV/QL1Gu3RtK6Cx7QK7A0OoOnGXHcaOUdfubosjizWlHardnVahLYyk0BaDFqSbaWMRQXRnFEn40WCglBAUbmFLWljdEQUD7/3je+J+/XJzP/fGe/K9lzwfM3dyzud8zvm+bya5r/s5Pz4nVYUkSTN5zrgLkCQtfIaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU0HjbuA+XL44YfXihUrxl2GJC0qt9566z9V1USr37MmLFasWMHWrVvHXYYkLSpJvj6bfp6GkiQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNT1rnuBeyNZdsHHej3nReWftt+NLkiMLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTb2GRZLVSe5JsiPJudNs/7UkdyXZluSGJC8e2vZUkju6r0191ilJmllvEwkmWQJcDJwMTAFbkmyqqruGut0OTFbV40l+FfgwcGa37dtVdVxf9UmSZq/PkcUqYEdV7ayqJ4GrgFOHO1TV56vq8W71ZmBZj/VIkvZRn2GxFHhwaH2qa9ubs4HPDK0fmmRrkpuTnNZHgZKk2enzfRaZpq2m7Zi8DZgEfm6oeXlVPZTkKODGJHdW1X0j+50DnAOwfPny+alakvQMfY4spoAjhtaXAQ+NdkpyEvAB4JSqemJPe1U91P25E7gJOH5036q6tKomq2pyYmJifquXJH1Pn2GxBViZ5MgkBwNrgO+7qynJ8cAlDILikaH2w5Ic0i0fDpwADF8YlyTtR72dhqqq3UnWAdcBS4ANVbU9yXpga1VtAj4CPB/4iyQAD1TVKcDLgUuSPM0g0C4cuYtKkrQf9foO7qraDGweaTt/aPmkvez3JeCVfdYmSZo9n+CWJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJauo1LJKsTnJPkh1Jzp1m+68luSvJtiQ3JHnx0La1Se7tvtb2WackaWa9hUWSJcDFwFuAY4Czkhwz0u12YLKqXgVcDXy42/eFwAeB1wGrgA8mOayvWiVJM+tzZLEK2FFVO6vqSeAq4NThDlX1+ap6vFu9GVjWLb8ZuL6qHq2qx4DrgdU91ipJmkGfYbEUeHBofapr25uzgc/MZd8k5yTZmmTrrl27fsByJUl702dYZJq2mrZj8jZgEvjIXPatqkurarKqJicmJva5UEnSzPoMiyngiKH1ZcBDo52SnAR8ADilqp6Yy76SpP2jz7DYAqxMcmSSg4E1wKbhDkmOBy5hEBSPDG26DnhTksO6C9tv6tokSWNwUF8HrqrdSdYx+CG/BNhQVduTrAe2VtUmBqedng/8RRKAB6rqlKp6NMlvMQgcgPVV9WhftUqSZtZbWABU1WZg80jb+UPLJ82w7wZgQ3/VSZJmyye4JUlNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSk2EhSWoyLCRJTb0+wa3Fb90FG3s57kXnndXLcSX1w5GFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpKZewyLJ6iT3JNmR5Nxptr8hyW1Jdic5fWTbU0nu6L429VmnJGlmB/V14CRLgIuBk4EpYEuSTVV111C3B4B3AP9lmkN8u6qO66s+SdLszTiySHL50PLaOR57FbCjqnZW1ZPAVcCpwx2q6mtVtQ14eo7HliTtR63TUK8eWn7vHI+9FHhwaH2qa5utQ5NsTXJzktOm65DknK7P1l27ds2xPEnSbLXCon6AY+cHPN7yqpoEfgX4/SQ/9YyDVV1aVZNVNTkxMbGvdUqSGlrXLJYl+SiDH/x7lr+nqt4zw75TwBHDxwIemm1hVfVQ9+fOJDcBxwP3zXZ/SdL8aYXFfx1a3jrHY28BViY5EvgHYA2DUUJTksOAx6vqiSSHAycAH57j50uS5smMYVFVV+zrgatqd5J1wHXAEmBDVW1Psh7YWlWbkrwWuBY4DPg3ST5UVccCLwcuSfI0g1NlF47cRSVJ2o9mDIvW8w1VdUpj+2Zg80jb+UPLWxicnhrd70vAK2c6tiRp/2mdhvoZBnc0bQS+wvQXrSVJz3KtsPgJBg/VncXgesNfAxuranvfhUmSFo4Zb52tqqeq6m+qai3wemAHcFOSd++X6iRJC0Jzuo8khwC/xGB0sQL4KPDpfsuSJC0krQvcVwCvAD4DfKiqvrpfqpIkLSitkcXbgW8BRwPvTbLnCewAVVU/0mdxkqSFofWcxQHzvot1F2yc92NedN5Z835MSRqH1mmoQ4F3Ai8BtjF4sG73/ihMkrRwtEYOVwCTwJ3ALwL/vfeKJEkLTuuaxTFV9UqAJJcBt/RfkiRpoWmNLL67Z8HTT5J04GqNLF6d5JvdcoDnduveDSVJB5DW3VBL9lchkqSF64C5NVaStO8MC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktTUa1gkWZ3kniQ7kpw7zfY3JLktye4kp49sW5vk3u5rbZ91SpJm1ltYJFkCXAy8BTgGOCvJMSPdHgDeAXxiZN8XAh8EXgesAj6Y5LC+apUkzazPkcUqYEdV7ayqJ4GrgFOHO1TV16pqG/D0yL5vBq6vqker6jHgemB1j7VKkmbQZ1gsBR4cWp/q2uZt3yTnJNmaZOuuXbv2uVBJ0sz6DItM01bzuW9VXVpVk1U1OTExMafiJEmz12dYTAFHDK0vAx7aD/tKkuZZn2GxBViZ5MgkBwNrgE2z3Pc64E1JDusubL+pa5MkjUFvYVFVu4F1DH7I3w18qqq2J1mf5BSAJK9NMgWcAVySZHu376PAbzEInC3A+q5NkjQGB/V58KraDGweaTt/aHkLg1NM0+27AdjQZ32SpNnxCW5JUpNhIUlq6vU0lNSy7oKN837Mi847a96PKR3oHFlIkpoMC0lSk2EhSWoyLCRJTYaFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCbDQpLUZFhIkpoMC0lSU69hkWR1knuS7Ehy7jTbD0nyyW77V5Ks6NpXJPl2kju6rz/us05J0swO6uvASZYAFwMnA1PAliSbququoW5nA49V1UuSrAF+Fziz23ZfVR3XV32SpNnrc2SxCthRVTur6kngKuDUkT6nAld0y1cDJyZJjzVJkvZBn2GxFHhwaH2qa5u2T1XtBr4BvKjbdmSS25N8Icm/7rFOSVJDb6ehgOlGCDXLPg8Dy6vqn5O8BvjLJMdW1Te/b+fkHOAcgOXLl89DyZKk6fQ5spgCjhhaXwY8tLc+SQ4CXgA8WlVPVNU/A1TVrcB9wNGjH1BVl1bVZFVNTkxM9PAtSJKg37DYAqxMcmSSg4E1wKaRPpuAtd3y6cCNVVVJJroL5CQ5ClgJ7OyxVknSDHo7DVVVu5OsA64DlgAbqmp7kvXA1qraBFwG/FmSHcCjDAIF4A3A+iS7gaeAd1bVo33VKkmaWZ/XLKiqzcDmkbbzh5a/A5wxzX7XANf0WZskafZ8gluS1GRYSJKaDAtJUpNhIUlqMiwkSU293g0ljdu6CzbO+zEvOu+seT+mtNA5spAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU2GhSSpybCQJDUZFpKkJsNCktRkWEiSmgwLSVKTYSFJajIsJElNhoUkqcmwkCQ1GRaSpCZfqyrtI1/ZqgOJIwtJUpNhIUlqMiwkSU2GhSSpqdewSLI6yT1JdiQ5d5rthyT5ZLf9K0lWDG17f9d+T5I391mnJGlmvd0NlWQJcDFwMjAFbEmyqaruGup2NvBYVb0kyRrgd4EzkxwDrAGOBf4l8LkkR1fVU33VKy003m2lhaTPkcUqYEdV7ayqJ4GrgFNH+pwKXNEtXw2cmCRd+1VV9URV3Q/s6I4nSRqDPp+zWAo8OLQ+Bbxub32qaneSbwAv6tpvHtl3aX+lSgeWPkYt4Mjl2SxV1c+BkzOAN1fVf+zW3w6sqqp3D/XZ3vWZ6tbvYzCCWA98uar+vGu/DNhcVdeMfMY5wDnd6kuBe3r5Zr7f4cA/7YfP6cNirh0Wd/3WPh7W3vbiqppodepzZDEFHDG0vgx4aC99ppIcBLwAeHSW+1JVlwKXzmPNTUm2VtXk/vzM+bKYa4fFXb+1j4e1z58+r1lsAVYmOTLJwQwuWG8a6bMJWNstnw7cWIOhziZgTXe31JHASuCWHmuVJM2gt5FFdw1iHXAdsATYUFXbk6wHtlbVJuAy4M+S7GAwoljT7bs9yaeAu4DdwLu8E0qSxqfXiQSrajOweaTt/KHl7wBn7GXf3wZ+u8/69tF+Pe01zxZz7bC467f28bD2edLbBW5J0rOH031IkpoMizloTV+yUCU5Isnnk9ydZHuS9467prlKsiTJ7Un+aty1zEWSH01ydZK/6/7+f2bcNc1Fkvd1/2a+mmRjkkPHXdPeJNmQ5JEkXx1qe2GS65Pc2/152Dhr3Ju91P6R7t/NtiTXJvnRcdZoWMzS0PQlbwGOAc7qpiVZDHYD/7mqXg68HnjXIqp9j/cCd4+7iH3wB8DfVNXLgFeziL6HJEuB9wCTVfUKBjeqrBlvVTO6HFg90nYucENVrQRu6NYXost5Zu3XA6+oqlcBfw+8f38XNcywmL3ZTF+yIFXVw1V1W7f8fxn8wFo0T8QnWQb8EvDxcdcyF0l+BHgDg7v+qKonq+r/jLeqOTsIeG73HNTzmOZ5p4Wiqr7I4K7KYcNTCl0BnLZfi5ql6Wqvqs9W1e5u9WYGz5uNjWExe9NNX7JofuDu0c3sezzwlfFWMie/D/w68PS4C5mjo4BdwJ90p9A+nuSHxl3UbFXVPwC/BzwAPAx8o6o+O96q5uzHq+phGPzSBPzYmOvZV/8B+Mw4CzAsZi/TtC2qW8mSPB+4BvhPVfXNcdczG0l+GXikqm4ddy374CDgp4E/qqrjgW+xcE+DPEN3fv9U4EgGsz//UJK3jbeqA0+SDzA4lXzlOOswLGZvVlOQLFRJ/gWDoLiyqj497nrm4ATglCRfY3Dq741J/ny8Jc3aFDBVVXtGcVczCI/F4iTg/qraVVXfBT4N/Ksx1zRX/5jkJwG6Px8Zcz1zkmQt8MvAW2vMzzkYFrM3m+lLFqRu2vfLgLur6n+Mu565qKr3V9WyqlrB4O/8xqpaFL/dVtX/Bh5M8tKu6UQGsxIsFg8Ar0/yvO7f0Iksogv0neEphdYC/3OMtcxJktXAbwCnVNXj467HsJil7kLTnulL7gY+VVXbx1vVrJ0AvJ3Bb+V3dF+/OO6iDhDvBq5Msg04DvidMdcza92I6GrgNuBOBj8vFtRTxcOSbAS+DLw0yVSSs4ELgZOT3MvgRWwXjrPGvdlL7RcBPwxc3/2f/eOx1ugT3JKkFkcWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMix0wEjy40k+kWRnkluTfDnJv93HY60YniF0oUryjiQXjbsOLX6GhQ4I3UNlfwl8saqOqqrXMHjIb6yTs+2LbgZkab8yLHSgeCPwZFV978Gmqvp6VX0MIMmhSf4kyZ3dpH+/0LWvSPK/ktzWfT1juoskxya5pXtwaluSlSPblyS5vHsnxJ1J3te1vyTJ55L8bXfsn8rAR4b6ntn1/fnunSSfYPCAHEneNvS5l+wJkST/PsnfJ/kCgwcypR9Yr+/glhaQYxk8ibw37wKoqlcmeRnw2SRHM5hL6OSq+k4XAhuByZF93wn8QVVd2U0FM/qb/3HA0u6dEAy9xOZK4MKqurZ7qdBzgH/X9X81cDiwJckXu/6rGLzf4P4kLwfOBE6oqu8m+UPgrUmuBz4EvAb4BvB54PbZ/iVJe2NY6ICU5GLgZxmMNl7bLX8MoKr+LsnXgaOBrwMXJTkOeKprG/Vl4APdezc+XVX3jmzfCRyV5GPAXzMIoh9mECDXdp/5na6unwU2VtVTDCbB+wLwWuCbwC1VdX93zBMZBMKWwRk2nssg2F4H3FRVu7rjfXIvNUtz4mkoHSi2MzTja1W9i8EP3Imuabop6AHeB/wjg9/0J4GDRztU1SeAU4BvA9cleePI9se6/W9iMIL5+Ayft7d2GExxPtzviqo6rvt6aVX95p6PnOEY0j4xLHSguBE4NMmvDrU9b2j5i8BbAbrTT8uBe4AXAA9X1dMMJmN8xsXlJEcBO6vqowxmOX3VyPbDgedU1TXAfwN+unufyFSS07o+hyR5XlfHmd11jgkGb9q7ZZrv5wbg9CQ/1u3/wiQvZvBSq59P8qJuWvozZv9XJO2dYaEDQvcugNOAn0tyf5JbGLxm8ze6Ln8ILElyJ/BJ4B1V9UTXvjbJzQxO53zrmUfnTOCrSe4AXgb86cj2pcBN3fbL+f/vUn478J5uRtovAT8BXAtsA/6WQcD9ejfV+ej3cxdwHoNTWtsYvK/5J7u3wf0mg1Njn2Pm6zTSrDnrrCSpyZGFJKnJsJAkNRkWkqQmw0KS1GRYSJKaDAtJUpNhIUlqMiwkSU3/Dx1Fy7KRu8iXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# The following function simulates many games, then uses the\n",
    "# number of goals scored as an estimate of the true long-term\n",
    "# goal-scoring rate.\n",
    "\n",
    "def Estimate6(lam=2, m=1000000):\n",
    "\n",
    "    estimates = []\n",
    "    for i in range(m):\n",
    "        L = SimulateGame(lam)\n",
    "        estimates.append(L)\n",
    "\n",
    "    print('Experiment 4')\n",
    "    print('rmse L', RMSE(estimates, lam))\n",
    "    print('mean error L', MeanError(estimates, lam))\n",
    "    \n",
    "    pmf = thinkstats2.Pmf(estimates)\n",
    "    thinkplot.Hist(pmf)\n",
    "    thinkplot.Config(xlabel='Goals scored', ylabel='PMF')\n",
    "    \n",
    "Estimate6()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) RMSE for this way of estimating lambda is 1.4\n",
    "\n",
    "# 2) The mean error is small and decreases with m, so this estimator\n",
    "#    appears to be unbiased.\n",
    "\n",
    "# One note: If the time between goals is exponential, the distribution\n",
    "# of goals scored in a game is Poisson.\n",
    "\n",
    "# See https://en.wikipedia.org/wiki/Poisson_distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}