code/chap08exmine.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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 1\n",
      "rmse xbar 0.37724940448\n",
      "rmse median 0.451699927125\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": {
    "collapsed": true
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean error biased -0.131580036087\n",
      "mean error unbiased 0.0131566245651\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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FfX1//HXIeyrssiXHVRQ0LpGbbVW6wrWlq+Cilvd\ngbrUvVChirVWcataF0RFXKpoXVqrVG2/VfurSzVaCyKKiCggKCj7loSc3x93iHdyc5MQMnfm5r6f\nj0ce5HNm7p03l0kOc2fuZ8zdERERAWgSdwAREUkONQUREamkpiAiIpXUFEREpJKagoiIVFJTEBGR\nSmoKIiJSSU1BREQqqSmIiEilpnEH2FKdO3f2vn37xh1DRCSvvPPOO8vcvUtt6+VdU+jbty8lJSVx\nxxARyStm9lld1tPbRyIiUklNQUREKqkpiIhIJTUFERGppKYgIiKVImsKZjbFzL4ys/ezLDczu93M\n5prZDDPbK6osIiJSN1FekjoVuAN4KMvyIUD/4Gs/4O7gTxGRxHF31m0ojWXb8xd9zdr1pbRq0Yy+\nPTrRrk3LyLYVWVNw93+aWd8aVhkKPOSp+4G+aWbbmFk3d18cVSYRkfLyTazfWMZnX3zNzI+/oEWz\nb38Nlsz6jPZtWrKxtJwZcxbStVN7zGDJslUxJg4bN+oo9hrUO7Lnj/PDaz2ABWnjhUEtoymY2Uhg\nJEDv3tG9GCKS/1auXs+/3p3LX16ewQ69OoeWvTnj0y16ri+/Tk4z2Kxt6xaRPn9efKLZ3ScDkwGK\ni4s95jgikgAbNpaxfmMZpWXlfL74G75ctooHnnk9tM7S5asj2Xarls0jed6aGLDvbv3o32e7SLcT\nZ1NYBPRKG/cMaiIilJaVM2POIj6at4QmRU0wg1f+PafBftG3bNGMDRvL6NqpPd/bY/vK+qo1G9ip\nX1datWxO29Yt2K5ju8plXbZtS9OmRQ2y/aSKsyk8C5xvZtNInWBeqfMJIrJ81TrO+dVDNMRbAv16\ndmbY4eELG4uKmlC8S2+aNNEV+dWJrCmY2WPAwUBnM1sIXAU0A3D3ScB04ChgLrAOOCOqLCKSfKvX\nbuD0K6Zu8eM6tGvFytXrada0iN7dOtJl27Zcesbh+qVfT1FefXRiLcsdOC+q7YtIfthYWsbVdz3P\nR58uybpO9y4d2G2nnrRr25J160sp3qUPfbp3okO7VjlMWhjy4kSziDQ+a9dv5MLfPs7yVeuyrnPL\nmOPo071TDlOJmoKI5FxtbxU9MvHMWK7wETUFEcmxBUuWc9F1j1e77IZLh7FD71pvDiYRUlMQkZxY\n+OVyLvxt9c1AbxMlh5qCiETG3bnniX/yt9dnZ13n8ZvPafTX/ucTNQURaXClZeVMefq1GpvBLjt2\n56pzj6aoSJeOJomagog0qL+8PIOpf3o96/LuXTpw4+XDadmiWQ5TSV2pKYhIg/lw3pKsDeGsYQcw\n5MBdMbMcp5ItoaYgIg1m3G1/yqhdcvrhHLDnDjGkkfpQUxCRrbZpUwXHXzI5o/7UbaNjSCNbQ01B\nROrt6xVr+NXtz1Z734Enbx0VQyLZWmoKIrLFli1fw6gJj2RdfuQBu+jcQZ5SUxCRLXLaLx9gzbqN\nWZdPvvoUOm3TNoeJpCGpKYhInQ27cFLWZfpUcuOgpiAidXLHoy9XW3/ilpH6AFojoqYgIlmVl2/i\nrffnc/MDf8tY9usLfsIuO3aPIZVESU1BRKr10adLuOLWzM8dAAw9ZHc1hEZKTUFEQt794HOuvWd6\n1uUjjtqH4UfslXW55Dc1BRGp9N6HC7I2hEE7dOPCUw+l87a6sqgxU1MQESB1r+Rr7n4+o960aRHT\nbjpbnzsoEGoKIkJ5+SZOuvz+jPq0m86hWTPd66CQqCmICCdcem9GTfMWFSZdXCxS4J57ZUZG7aHr\nz4ghiSSBjhRECpS7M/yiezLqj910Ns2b6VdDodKRgkiBGn/7nzNqe+zcSw2hwOlfX6QALViynA/n\nLQnVThhSzPGDi2NKJEmhpiBSYMrLN3HRdY+HaleeezS779QzpkSSJGoKIgXkjffmcdMDL2XU1RBk\nM51TECkQX3y1otqGoEtPJZ2agkgBWL+hlAuunZZR16WnUlWkTcHMBpvZR2Y218zGVrO8g5n9xcz+\na2azzEx7qEgEJt7/YkbtqdtG06ZVixjSSJJF1hTMrAi4ExgCDAJONLNBVVY7D/jA3XcHDgZuNrPm\nUWUSKUTuzsw5i0I1vWUk2UR5pLAvMNfd57l7KTANGFplHQfaWWqmrbbAN0B5hJlECs7Iqx4Jja+7\n+JiYkkg+iLIp9AAWpI0XBrV0dwADgS+AmcCF7l4RYSaRglJaVs43K9eGagP6do0pjeSDuE80Hwm8\nB3QH9gDuMLP2VVcys5FmVmJmJUuXLs11RpG8VFa2iRMvuy9Uu+HSYTGlkXwRZVNYBPRKG/cMaunO\nAJ72lLnAp8DOVZ/I3Se7e7G7F3fp0iWywCKNyYjLMmc+3aG3fn6kZlE2hbeB/mbWLzh5PAJ4tso6\nnwOHAphZV2AnYF6EmUQK1tTfnh53BMkDkX2i2d3Lzex84EWgCJji7rPMbHSwfBJwDTDVzGYCBoxx\n92VRZRIpFBdf/0RofMf4E2nXpmVMaSSfRDrNhbtPB6ZXqU1K+/4L4IgoM4gUmtKycj5f/E2o1q1L\nh5jSSL6J+0SziDSwqieXr79El6BK3akpiDQiZ//qoYxa/z66BFXqTk1BpJHYWFrG8lXrQrWHrz8z\npjSSr9QURBqJky6/PzQeN+ooWrfSrDGyZdQURBqBq+6oerU37DWodwxJJN+pKYjkuZlzFvH+x1+E\napoSW+pLTUEkz0248y+h8UH7DNCU2FJvagoieWr12g0Mu3BSqNa7W0d+fsohMSWSxkBNQSRPnX7F\n1IzaLWOOy30QaVTUFETy0Oq1GzJqh+8/kNStSUTqL9JpLkQkGmNveTo0fvj6M3X5qTQIHSmI5Bl3\nZ8myVaGaGoI0FDUFkTxzx6OvhMYTzvtxLDmkcVJTEMkjFRUVvPLWR6HadwZUvcutSP2pKYjkCXfn\nuIsnh2q/vuAnMaWRxkpNQSRPDL/onozaLjt2jyGJNGZqCiJ5oLpLUG8fNyKGJNLYqSmI5IGqH1S7\nZcxx9Nhum3jCSKOmpiCScO6eUevTvVMMSaQQqCmIJNzPrn40NL5nwikxJZFCoKYgkmArV69n6fLV\noVrnbdvGlEYKgZqCSEJVVFRw5vgHQzV9UE2ipqYgklDX3D09NG7TqoU+qCaRU1MQSaCVq9czY87C\nUO2O8boEVaKnpiCSQFXfNvrFWUfSvm2rmNJIIVFTEEmYioqKjNp+u/WLIYkUIjUFkYSpOr/R1N+e\nHk8QKUhqCiIJcuejr2TU2rVpmfMcUrjUFEQS5B///jA0vuvKk2JKIoVKTUEkIf7+xuzQ+JzhB9K1\nU/uY0kihUlMQSYi7p70aGg8+cJeYkkghi7QpmNlgM/vIzOaa2dgs6xxsZu+Z2Swze7W6dUQau2XL\n14TG+++5Q0xJpNA1jeqJzawIuBM4HFgIvG1mz7r7B2nrbAPcBQx298/NbLuo8ogk2agJj4TGPz/5\nkJiSSKGL8khhX2Cuu89z91JgGjC0yjonAU+7++cA7v5VhHlEEum2h/8vo9asWVEMSUSibQo9gAVp\n44VBLd0AYFsze8XM3jGzn1b3RGY20sxKzKxk6dKlEcUVicc/Sz4Ojf9ww1kxJRGJ/0RzU2Bv4EfA\nkcCvzGxA1ZXcfbK7F7t7cZcuXXKdUSQysz9ZHBr/oLg/LVs0iymNSITnFIBFQK+0cc+glm4h8LW7\nrwXWmtk/gd2BORHmEkmEiooKxt/+51Dt56foXILEK8ojhbeB/mbWz8yaAyOAZ6us82fg+2bW1Mxa\nA/sBsxEpAFWnswAwsxiSiHwrsiMFdy83s/OBF4EiYIq7zzKz0cHySe4+28xeAGYAFcB97v5+VJlE\nkmLpN6szak/dNjqGJCJhVt1NwZOsuLjYS0pK4o4hslWGXTgpNL5lzHH06d4ppjRSCMzsHXcvrm29\nuE80ixSctes3ZtTUECQp1BREcmzszU+Hxo/ddHZMSUQyqSmI5NgXS1eGxs2bRXkRoMiWUVMQyaET\nL7svNB57zuCYkohUT01BJEdKy8opLSsP1fbZtW88YUSyUFMQyZHbHgrPcXTluUfHlEQkOzUFkRx5\nc8anofHuO/WMKYlIdjU2BTObmvb9aZGnEWmknnzp3dB41PE/iCmJSM1qO1LYPe37C6MMItJYlZaV\n89jzb4Vqh+8/MKY0IjWrrSnk18edRRLojkdfCY2PH1ysOY4ksWq7QLqnmd0OWNr3ldz955ElE2kk\nXnt3bmh8wpBaZxoQiU1tTeHytO814ZDIFjr/N4+FxuNGHRVTEpG6qbEpuPuDuQoi0ti8+8HnLK7y\n6eU9B/bKsrZIMtR6SaqZnWZm75rZ2uCrJNttM0XkW9feMz00PuOY/XUuQRKvxiOF4DLUi4BLgHdJ\nnVvYC7jRzNzdH44+okj+qTol/bbtW3P0wbvFlEak7mo7UvgZcIy7v+zuK919hbv/AxgGnBd9PJH8\nNPyie0Lje399akxJRLZMbU2hvbvPr1oMau2jCCSS76reQAd0m03JH7U1hfX1XCZSkKq7k+EtY46L\nIYlI/dR2SepAM5tRTd2A7SPII5LXbrj/xdD4ynOP1l3VJK/U1hR2B7oCC6rUewFLIkkkkqfcnbdm\nzg/VNOmd5Jva3j76HbDS3T9L/wJWBstEJPD54uWh8Y2XDYspiUj91dYUurr7zKrFoNY3kkQieeqS\niU+Extv36hJTEpH6q60pbFPDslYNGUQkny1Ysrz2lUTyQG1NocTMzqlaNLOzgXeiiSSSfy667vHQ\nWJ9LkHxV24nmi4BnzOxkvm0CxUBz4Jgog4nkq7atW9CxQ5u4Y4jUS20T4n0J7G9mPwR2DcrPB59q\nFhHgw3nhC/Fuu+KEmJKIbL3ajhQAcPeXgZcjziKSl8bd9qfQeJt2rWNKIrL1ap0lVUSy27SpIjQe\nuH23mJKINAw1BZGtcPwlk0Pj8aN1Ex3Jb2oKIvU0Z/6XGbWWLZrFkESk4UTaFMxssJl9ZGZzzWxs\nDevtY2blZjY8yjwiDeXDeUv45e+eCdV+N/b4mNKINJzImoKZFQF3AkOAQcCJZjYoy3oTgZeiyiLS\nkMrLN2WcXO7epQO9u3WMKZFIw4nySGFfYK67z3P3UmAaMLSa9S4AngK+ijCLSIM5ecyU0LjTNm24\nfdyImNKINKwom0IPwrOrLgxqlcysB6kPwd0dYQ6RBuPulJdvCtUmX32qbqIjjUbcJ5pvBca4e0VN\nK5nZSDMrMbOSpUuX5iiaSKaqt9n8vY4QpJGp04fX6mkRqfsubNYzqKUrBqYF/8vqDBxlZuXuHnrD\n1t0nA5MBiouLM29tJZIDY25+OqPWfbua5owUyT9RNoW3gf5m1o9UMxgBnJS+grv32/y9mU0Fnqva\nEESS4C8vz2Du5+HTXpOuOjmmNCLRiawpuHu5mZ0PvAgUAVPcfZaZjQ6WZ97dXCShpv7p9dD4/JN+\nSJeO7WJKIxKdKI8UcPfpwPQqtWqbgbufHmUWkfpa9NWK0Lh4lz78cL+dYkojEq24TzSLJN7Pr50W\nGv/irCNjSiISPTUFkRp8vWJNRq2oSD820nhp7xapwZW/fzY0fvzmjBsRijQqagoiNViybFVo3LRp\nUUxJRHJDTUEki5JZn4XGEy85NqYkIrmjpiCSxXWT/xoa79hnu5iSiOSOmoJINV59e05oXLxLn5iS\niOSWmoJINW5/5B+h8WVnHBFTEpHcUlMQqeKqO8JXHO09qA/NmukEsxQGNQWRNM+9MoP3P/4iVLti\n1JCY0ojknpqCSJoHngnPcXTCkOKYkojEQ01BJLBy9frQeK9BvTl+sJqCFBY1BRFS910+c/yDodq4\nUUfFlEYkPmoKIsAJl94bdwSRRFBTkII37MLM2dzvu+anMSQRiZ+aghS0qh9SA7h93Ai2bd86hjQi\n8VNTkIJW9UNqV517ND1032UpYGoKUrBWrF4XGu/Yezt226lnTGlEkkFNQQrWWeMfCo2vPPdHMSUR\nSQ41BSlIi5euzKi1adUihiQiyaKmIAVn/YZSzv/NY6Hak7eOiimNSLKoKUjBOWXMlNC4eJc+mFlM\naUSSRU1BCspbM+dn1MaeMzj3QUQSSk1BCsrE+14IjR+76WwdJYikUVOQgvGvd+aGxgfstSPNmzWN\nKY1IMqkpSMH43UN/D40v/umhMSURSS41BSkIoyY8Ehr377Od3jYSqYaagjR6pWXlLFu+JlS7/pJj\nY0ojkmxqCtLonXjZfaHxtRf+b0xJRJJPTUEatdET/pBR23n7/4khiUh+UFOQRmvF6nUsXb46VHv8\n5nNiSiOSHyJtCmY22Mw+MrO5Zja2muUnm9kMM5tpZq+b2e5R5pHCUnXCu3GjjqJp06KY0ojkh8ia\ngpkVAXcCQ4BBwIlmNqjKap8CB7n7d4BrgMlR5ZHCcv9T/8qo7TWodwxJRPJLlEcK+wJz3X2eu5cC\n04Ch6Su4++vuvjwYvgloMnvZal+vWMP0f74fqt10+fCY0ojklyibQg9gQdp4YVDL5izgr9UtMLOR\nZlZiZiVLly5twIjSGN097dXQeLcBPenXs3NMaUTySyJONJvZD0k1hTHVLXf3ye5e7O7FXbp0yW04\nySsbS8v4z+wFodpV5x0dUxqR/BPlxC+LgF5p455BLcTMdgPuA4a4+9cR5pECcNLl94fG1118TExJ\nRPJTlEcKbwP9zayfmTUHRgDPpq9gZr2Bp4FT3X1OhFmkAAy7cFJGbUDfrjEkEclfkR0puHu5mZ0P\nvAgUAVPcfZaZjQ6WTwKuBDoBdwXz0JS7e3FUmaTxuvSGJzNq9/761BiSiOS3SOcNdvfpwPQqtUlp\n358NnB1lBmn8Jt73AvMXLQvVrjz3aDp2aBNTIpH8lYgTzSL19dU3qzPupnbuiQex+066ulmkPtQU\nJK/97Orw3EY/Png3Dv3uwJjSiOQ/NQXJW5dM/GNG7fRj9o8hiUjjoaYgeem9Dxfw2RfhK5g12Z3I\n1lNTkLyzYvU6rrn7+VBtwnk/1mR3Ig1ATUHyyvxFyzJmP23fthXfGVDTDCoiUldqCpI3vvx6VbWf\nR3jg2tNiSCPSOKkpSF4oLSvn3F8/mlF/8tZRMaQRabwi/fCaSEP4z+wF/GbS8xn1J24ZSfBJeBFp\nIGoKkmibNlVU2xCevHWUGoJIBPT2kSTa8Zdk3ozv9+NGqCGIRERHCpJYy1ety6g9ddvoGJKIFA4d\nKUhinf2r8KWn90w4JaYkIoVDTUESaeXq9Rm1ztu2jSGJSGFRU5DEWbe+lDPHPxiqPXrjWTGlESks\nagqSOKeOnRIaN21aRIvmzWJKI1JYdKJZEuOtmfOZeN8LGfWHrz8jhjQihUlNQWJXUVHBcRdnXnoK\n8MjEM2neTLupSK7op01iU1FRwRMvvsMfX3in2uU3XDqMVi2b5ziVSGFTU5CcKyvbxGU3PsnCL5dX\nu3zPgb244ORD6NCuVY6TiYiaguTUEy+U8PhfS7Iuv2fCKbr0VCRGagqSE7M/Wcz42/+cdfmxh+3J\n8YOLadZMN8oRiZOagjS4DRvLWLJsJbPmLmbK06/VuO6vfvYj9ti5V46SiUht1BSkwZSVbWLEZffW\nad1D9tuZ8046ONpAIrLF1BRkq3yzci1X3/lc1pPGVe02oCdXnvsjzXIqklBqCrJFyss3MW3622yq\ncJ59+b+1rt+xQxu269SOs4d9n749OqkZiCScmoLUibtzypgpbNhYVqf1b7viBHp23TbiVCLS0NQU\nJGTl6vVsKC2jrHwT78/5gsdfKGHVmswZS6sz+oQfcNA+A/QJZJE8pp/eArVhYxkff/YV36xcy8yP\nF1HUpAl/f2P2Fj3H8CP3pnmzIo49bE+9LSTSSKgpFJD1G0q587FXeeO9T7bqeW4Zczx9undsoFQi\nkiSRNgUzGwzcBhQB97n79VWWW7D8KGAdcLq7vxtlpkKzbPka3nhvHlP/9PoWPa7Ltu1YsXodZeWb\n6NujM8OP2Is9B/aiZQtNYS3SmEXWFMysCLgTOBxYCLxtZs+6+wdpqw0B+gdf+wF3B39KHXy6cBnz\nF33NN6vW8uhzb7Fj7+0oKmrCxtJy5i9aVufn2Xn7/2H9hjIO2GsHOnVow0H7DNDbQSIFKsojhX2B\nue4+D8DMpgFDgfSmMBR4yN0deNPMtjGzbu6+uKHDvPafT3jxX7MoLStv6Keul48/+wqAJk3qd5+j\nioqKjNrcz7+q8+N/OvR7DD1k93ptW0QaryibQg9gQdp4IZlHAdWt0wNo0KZQWlbOXY+9UufLKXOp\nul/uURm4fTcuO/NwtmnXOmfbFJH8khcnms1sJDASoHfv3lv8+PUbyhLZEBrK3oP60Lx5U8rKyvnJ\nIbtTFBx9dNqmDdu2b03TpppkTkTqJsqmsAhIn+msZ1Db0nVw98nAZIDi4mLf0iCtWjbj7OHf55sV\na1m5Zj2HfW/glj5FJNq0bkHXju226jn0C19EGlKUTeFtoL+Z9SP1i34EcFKVdZ4Fzg/ON+wHrIzi\nfELzZk0ZcuCuDf20IiKNTmRNwd3Lzex84EVSl6ROcfdZZjY6WD4JmE7qctS5pC5J1R3aRURiFOk5\nBXefTuoXf3ptUtr3DpwXZQYREam7+l0PKSIijZKagoiIVFJTEBGRSmoKIiJSSU1BREQqWeoCoPxh\nZkuBz7bwYZ2Bus8QFw9lbBjK2DCSnjHp+SB5Gfu4e5faVsq7plAfZlbi7sVx56iJMjYMZWwYSc+Y\n9HyQHxmro7ePRESkkpqCiIhUKpSmMDnuAHWgjA1DGRtG0jMmPR/kR8YMBXFOQURE6qZQjhRERKQO\n8r4pmNnFZjbLzN43s8fMrKWZdTSzv5nZx8Gf22Z57GAz+8jM5prZ2Bzmu9HMPjSzGWb2jJltk+Wx\n881sppm9Z2YlUeSrIeMEM1sUbPs9Mzsqy2Mjfw1ryPh4Wr75ZvZelsfm6nW8MMg3y8wuCmqJ2Rdr\nyJi0/bG6jEnbH6vLmKj9sd7cPW+/SN2681OgVTB+AjgduAEYG9TGAhOreWwR8AmwPdAc+C8wKEf5\njgCaBrWJ1eULls0HOsf0Gk4ALqvlsZG/hjVlrLLOzcCVMb6OuwLvA61JzT78d2DHpOyLtWRM0v6Y\nLWOS9sdqMyZpf9yar7w/UiD1j9LKzJqS+kf6AhgKPBgsfxD432oety8w193nuXspMC14XOT53P0l\ndy8Plr9J6o5zcaruNayLXL2GNWY0MwOOBx6LaNt1MRD4t7uvC/5tXwWOJVn7YrUZE7Y/Znsd6yLW\n13HzwoTsj/WW103B3RcBNwGfA4tJ3bntJaCrf3sHtyVA12oe3gNYkDZeGNRykS/dmcBfsz0F8Hcz\ne8dS96lucLVkvCB4S2FKlrc9In8N65AR4EDgS3f/ONtTEPHrSOp/jgeaWScza03q5lG9SMi+WEvG\ndLHuj7VkTMT+WEtGSMb+WG953RSCHWMo0A/oDrQxs1PS1/HU8Vosl1jVls/MxgHlwB+yPMX33X0P\nYAhwnpn9IIcZ7yZ1GL4HqV/ENzf0thsg42YnUvP/yiJ/Hd19Nqm3Xl4CXgDeAzZVWSe2fTHYfo0Z\nk7A/1pAxMftjHf6tY98ft0ZeNwXgMOBTd1/q7mXA08D+wJdm1g0g+POrah67iHB37xnUcpEPMzsd\nOBo4OfhlkSH4HzLu/hXwDKnD44ZWbUZ3/9LdN7l7BXBvlm3n4jXMmhEgeDvpWODxbA/O0euIu9/v\n7nu7+w+A5cAckrMv1pQxSftjtRkTtj/W9DomZn+sr3xvCp8D3zWz1sH7eIcCs4FngdOCdU4D/lzN\nY98G+ptZPzNrDowIHhd5PjMbDPwC+Im7r6vugWbWxszabf6e1MnA9xs4X00Zu6Wtc0yWbefiNcya\nMVh2GPChuy+s7oE5fB0xs+2CP3uT+sXwKMnZF7NmTNj+mC1jkvbHbP/WkKD9sd7iPtO9tV/A1cCH\npF7Yh4EWQCfg/4CPSV0Z0DFYtzswPe2xR5Hq8J8A43KYby6p9z7fC74mVc1H6lD5v8HXrKjy1ZDx\nYWAmMIPUD1a3uF7DbBmD+lRgdJV143od/x/wQbCtQ4NaYvbFGjImbX+sLmPS9seMjEnbH+v7pU80\ni4hIpXx/+0hERBqQmoKIiFRSUxARkUpqCiIiUklNQUREKqkpSKNgZuOCGStnBLNP7hfx9l4xs7y7\n/65IbZrGHUBka5nZ90h9Gncvd99oZp1JzZIpIltIRwrSGHQDlrn7RgB3X+buXwCY2ZVm9nYw9/3k\n4BPRm/+n/zszKzGz2Wa2j5k9ban7HvwmWKevpe4z8IdgnSeDCdBCzOwIM3vDzN41sz+aWdtq1ql1\ne8F6p5jZW8HRzj1mVhTU7w4eO8vMrk5bf76ZXR1se6aZ7dzAr60UGDUFaQxeAnqZ2Rwzu8vMDkpb\ndoe77+PuuwKtSB1RbFbq7sXAJFLTT5xHaq78082sU7DOTsBd7j4QWAWcm77h4KhkPHCYu+8FlACX\nZMlZ4/bMbCBwAnCApyZM2wScHDx2XPDY3YCDzGy3tOddFmz7buCy2l8ukezUFCTvufsaYG9gJLAU\neDyY4A3gh2b2bzObCRwC7JL20M3z4swEZrn74uBoYx7fTqy2wN1fC75/BPh+lc1/FxgEvGapO22d\nBvTJErW27R0a/D3eDp7rUFLTIgAcb2bvAv8J/g6D0p736eDPd4C+WbYtUic6pyCNgrtvAl4BXgka\nwGlmNg24Cyh29wVmNgFomfawjcGfFWnfbx5v/tmoOg9M1bEBf3P3E+sQs7btGfCgu/8ytAGzfqSO\nAPZx9+VmNjXL32MT+pmWraQjBcl7ZraTmfVPK+0BfMa3vziXBe/zD6/H0/cOTmQDnAT8q8ryN4ED\nzGzHIEsbMxtQj+1AauK84WkzcHY0sz5Ae2AtsNLMupKah18kEvpfhTQGbYHfW+qG8+WkZv0c6e4r\nzOxeUjMhQs8ZAAAAiElEQVSrLiE1tfKW+ojUjVCmkJoV8+70he6+NHir6jEzaxGUxxPMr78l3P0D\nMxsPvGRmTYAy4Dx3f9PM/kNqltgFwGs1PY/I1tAsqSJZmFlf4LngJLVIQdDbRyIiUklHCiIiUklH\nCiIiUklNQUREKqkpiIhIJTUFERGppKYgIiKV1BRERKTS/wcz/vXlsHXnUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbb40668>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "89.799813994632103"
      ]
     },
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(85.817601580482531, 94.016367498865662)"
      ]
     },
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4775077093404421"
      ]
     },
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse L 1.06797275063\n",
      "rmse Lm 1.81280848304\n",
      "mean error L 0.290422514039\n",
      "mean error Lm 0.39328844648\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": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean error 0.00101438588128\n",
      "mean error median 0.00089870197842\n"
     ]
    }
   ],
   "source": [
    "# Solution goes here\n",
    "iterations = 500\n",
    "sample_size = 8\n",
    "mu = 0\n",
    "sigma = 1\n",
    "means = []\n",
    "medians = []\n",
    "\n",
    "for i in range(iterations):\n",
    "    xs = [random.gauss(mu, sigma) for i in range(sample_size)]\n",
    "    xbar = np.mean(xs)\n",
    "    median = np.median(xs)\n",
    "    means.append(xbar)\n",
    "    medians.append(median)\n",
    "    \n",
    "print('mean error', MeanError(means, mu))\n",
    "print('mean error median', MeanError(medians, mu))\n",
    "print('finished')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "second experiment\n",
      "bias 0.897305218867\n",
      "unbias 1.0254916787\n",
      "finished\n"
     ]
    }
   ],
   "source": [
    "# Solution goes here\n",
    "print('second experiment')\n",
    "iterations = 500\n",
    "sample_size = 8\n",
    "mu = 0\n",
    "sigma = 1\n",
    "\n",
    "biased_est = []\n",
    "unbiased_est = []\n",
    "\n",
    "for i in range(iterations):\n",
    "    xs = [random.gauss(mu, sigma) for i in range(sample_size)]\n",
    "    bias = np.var(xs)\n",
    "    biased_est.append(bias)\n",
    "    unbias = np.var(xs, ddof=1)\n",
    "    unbiased_est.append(unbias)\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "print('bias', MeanError(biased_est, mu))\n",
    "print('unbias', MeanError(unbiased_est, mu))\n",
    "print('finished')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here\n"
   ]
  },
  {
   "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": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "experiment starting\n",
      "std error 2.47750770934\n",
      "90% confidence int = 3.20332430935\n",
      "experiment ending\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFOXV/vHvYRgWYQSEwQADQhQMRBTjCG8k7sqiRkVN\nFFdUVFTUiHkTjJqYRKNRXoNREBEJ8nMhimgMIeIS0URFAQUUECSIMiCyKbLOen5/TItTszAN09PV\nXX1/rovLOU/VVJ8Gual5qvopc3dERCRaGoTdgIiIJJ7CXUQkghTuIiIRpHAXEYkghbuISAQp3EVE\nIkjhLiISQQp3EZEIUriLiERQw7BeuE2bNt65c+ewXl5EJC3Nmzdvg7vn1rZfaOHeuXNn5s6dG9bL\ni4ikJTP7NJ79NC0jIhJBCncRkQhSuIuIRFBoc+7VKS4upqCggJ07d4bdSiQ0adKEvLw8srOzw25F\nRJIspcK9oKCAnJwcOnfujJmF3U5ac3c2btxIQUEBXbp0CbsdEUmyWqdlzGyima0zsw9r2G5m9mcz\nW25mC83sB3vbzM6dO2ndurWCPQHMjNatW+unIJEMFc+Z+yTgQWByDdsHAl1jv/oAD8X+u1cU7Imj\n30uR3XN3tu8sStrrrVy9kW07imjaOJvOHVqT06xJvb1WreHu7m+YWefd7HIGMNnLn9c328xamlk7\nd/88QT2KSAIVFZdQWFRS635rN2xm+87iJHRUvbIyZ+nKtTRq+G1Mvb1gBW1bNccaxHcvyNwPV9Kg\nQQNa7du0yrZ1G7dQFuJjRm+56hR+0KNTvR0/EXPuHYBVFeqC2FiVcDezK4ErATp1qr83lSyTJk2i\nX79+tG/fHoChQ4cyYsQIevToUafjrly5krfeeovzzz8/EW1KxBUWFbOzsPqwLi4pZeZ/FrHhq628\nMffjJHdWP1asWr+H31HK2g3h/SNVk+b7NK7X4yf1gqq7jwfGA+Tn56f9k7knTZrEIYccsivcJ0yY\nkJDjrly5kieffFLhnmF27CxiZw1n1AVrv+S9xZ+xcfM23nxvOa1bNgNg41fbktlipDVt0igpr2NA\n70O70PWAtvX6OokI99VAxwp1XmwsbT3++OP8+c9/pqioiD59+jB27Fguv/xy5s6di5lx2WWX0bFj\nR+bOncsFF1xA06ZNefvttxk4cCCjRo0iPz+f5s2bc/XVVzNjxgzatWvHH/7wB37xi1/w2WefMXr0\naE4//XRWrlzJRRddxLZt5X9BH3zwQY466ihGjhzJkiVL6NWrF5dccgnXX389I0eOZNasWRQWFnLt\ntddy1VVX8fnnn3Puuefy9ddfU1JSwkMPPcTRRx8d8u+e1GTr9kI++mQtZWVOUXEJi5av4f3Fq3Cc\nDV9u3aNjJSLUmzXd/ZljSWkZhUXFHNotr86vtbccZ+XqjRyT35XshlkAbN66g+8f2J5GjeKLLy9z\nOrbbj0bZWdVuz23VnIYNq9+WzhIR7i8Aw81sCuUXUjcnYr797BvG1bmxmjx7/7Aaty1ZsoS//vWv\nvPnmm2RnZ3PNNddwxx13sHr1aj78sPyGoa+++oqWLVvy4IMP7grzyrZt28YJJ5zAvffey6BBg7j1\n1lt5+eWXWbx4MZdccgmnn346bdu25eWXX6ZJkyZ8/PHHDB48mLlz53L33XczatQopk+fDsD48eNp\n0aIFc+bMobCwkL59+9KvXz+mTZtG//79ueWWWygtLWX79u318xsmcSsrK2P1us0s+ngNq9ZuolnT\nxsyaszQpZ9g1XZzbsq38jqlTjjmENq1yOPqIg9ivRbN670fCVWu4m9lTwHFAGzMrAH4DZAO4+zhg\nBnAKsBzYDlxaX80mw6uvvsq8efM48sgjAdixYwcDBgxgxYoVXHfddZx66qn069ev1uM0atSIAQMG\nANCzZ08aN25MdnY2PXv2ZOXKlUD5h7aGDx/O/PnzycrKYtmyZdUe66WXXmLhwoVMnToVgM2bN/Px\nxx9z5JFHctlll1FcXMyZZ55Jr169EvA7ILuzdsPXLP90HStXb6C0zGmY1YBpr7yf0NdokVP14l9p\naRlbtxfSumUz+h5+EG1aNad3z84ANGyYRat990loD5L+4rlbZnAt2x24NmEdhczdueSSS7jrrrsC\n43feeSczZ85k3LhxPP3000ycOHG3x8nOzt51K2KDBg1o3Ljxrq9LSsrnVf/0pz+x//77s2DBAsrK\nymjSpPozL3fngQceoH///lW2vfHGG/zjH/9gyJAhjBgxgosvvniP37NUr7i4lF8/+AI7C4vZtHkb\nW7cXJuS4jRtl07Nre0rLyli3cQuHd+9E9wO/wwHtW9Mut0VCXkMkpT6hWtHupk7q04knnsgZZ5zB\njTfeSNu2bdm0aRNbtmyhVatWnH322Rx88MFceOGFAOTk5LBly5a9fq3NmzeTl5dHgwYNeOyxxygt\nLa32uP379+ehhx7ihBNOIDs7m2XLltGhQwc2bNhAXl4eV1xxBYWFhbz33nsK9zpwd+Yt/ozP1mxi\n9bqvmPXu0jodr1F2Q7p1bsv3D2rPjp3FfL9rew7t1oFG2Sn7104iRP+XVdKjRw/uuOMO+vXrR1lZ\nGdnZ2dx3330MGjSIsrIygF1n9UOGDGHYsGG7LqjuqWuuuYazzz6byZMnM2DAAJo1K58HPfTQQ8nK\nyuKwww5jyJAh3HDDDaxcuZIf/OAHuDu5ubk8//zzzJo1i3vvvZfs7GyaN2/O5Mk1fc5MavL11h28\n8K8FPPfq/D36vn2bN6VNq+Yc1CmX1i2bU1hYzElHdadJo2ya79OYrCytySfhMg/pJv78/Hyv/LCO\nJUuW0L1791D6iSr9nsInBRtYuXrjrtrd+du/FlDwxZdxH+P6C0+g5b770Ll962rnxEWSxczmuXvV\nuzgq0Zm7RJK788Gy1fx27PS9+v7Du3ekSeNGHHJQe/r17U6DOD8RKZIqFO4SGc+98j6P//2dvf7+\nHx1xEEPP/lG9rvchkiwpF+7urgWvEiSsKbcw3HL/83y0Yu1u92mf24KunfffVZeVOacd25MDO+Xq\n/zmJnJQK9yZNmrBx40Yt+5sA36znXtPtlVGwbUchsxesYOxTr+92vwM75nLDxSfSoW3LJHUmEr6U\nCve8vDwKCgpYv35PFwaS6nzzJKYocXemvfI+T05/t8Z9vvfd73DTkJNpmdNUc+WSsVIq3LOzs/XU\nIKlWSUkp5970SK37DTqxFxee/j9J6EgktaVUuIvUpLZg7/7ddlx+dl+65LVJUkciqU3hLilt+44i\nLhpZ/VIPJ/T5HtcMPlbXZ0SqoXCXlLS7aZhHfneRVjUUqYXCXVLOi/9exCNT/13ttmsGH6tgF4mD\nwl1Syu7W8X9q1FAtuiUSJ/1NkZTx8NNvVDs++uZz6fidVknuRiS9KdwlJXz2+SZeenNxYOzw7h25\nddipIXUkkt4U7hK64uJSbrz76cDY1ecdy0k/zOzVLEXqQh/fk1C5O+f9vOpdMQp2kbpRuEtoNm3e\nxjk/e7jK+NTRV4XQjUi0aFpGQvHpmo2M+OMzVcafuOdyfShJJAEU7pJ0/5r9EWOemlVlfNxvLqBJ\n4+zkNyQSQQp3SRp3r3YaBsJ7ILpIVGnOXZJi/aYt1Qa7oWAXqQ86c5d698Gy1dw+5u9Vxs866XAu\n+HGfEDoSiT6Fu9Sr2QtWcO/El6qM33jxSfzoiINC6EgkMyjcpd7MencpDzzxWpXxqaOv0h0xIvVM\n4S714i/T3mL66wurjGt+XSQ5FO6ScNWt7JjbKodxt18QQjcimUl3y0hC1bRkr4JdJLl05i4JsfST\ntfxq9PNVxn86IJ9zB+aH0JFIZosr3M1sAHA/kAVMcPe7K21vATwOdIodc5S7/yXBvUoKKl/4awIl\nJaVVtt1+7Y/p2a1DCF2JSK3TMmaWBYwBBgI9gMFm1qPSbtcCi939MOA44P/MrFGCe5UU9Lux/6g2\n2H87XMEuEqZ4ztx7A8vdfQWAmU0BzgAqPlnBgRwrv7+tObAJKElwr5JiRt43jY8/XRcYO/qIrlx/\n4fE0aKDLOSJhiifcOwCrKtQFQOWPFT4IvACsAXKAc929LCEdSkpa/N/PqwT7ZWf15dRje4bUkYhU\nlKgLqv2B+cAJwIHAy2b2b3f/uuJOZnYlcCVAp06dEvTSkmyX3vIYX2/dERjLbZXDKcccElJHIlJZ\nPD87rwY6VqjzYmMVXQpM83LLgU+A71U+kLuPd/d8d8/Pzc3d254lRPc8OrNKsDcwY9ztF+hTpyIp\nJJ5wnwN0NbMusYuk51E+BVPRZ8CJAGa2P3AwsCKRjUpqeGfhJ1XGntGTk0RSTq3TMu5eYmbDgZmU\n3wo50d0Xmdmw2PZxwO+BSWb2AeWruP7S3TfUY98Sgr9MeytQ6yHWIqkrrjl3d58BzKg0Nq7C12uA\nfoltTVLJ6nVfVVkrRsEukrp0v5rUav2mLVx/55TA2E8H6FOnIqlMyw/Ibi1avoZfP1D5EgtaUkAk\nxSncpUavvbOUB5+suh77X//vihC6EZE9oXCXai1YWlBtsGs9dpH0oDl3qWJnYTG/Gzu9yvgzf7oy\nhG5EZG8o3KWKC37xaKBukdOUZ+8fpvViRNKI/rZKwOMvzK4y9ujvLw6hExGpC4W77PLaO0t57tX5\ngbEpo67QsgIiaUgXVAWA5155n8f//k5grM+hXcjOzgqpIxGpC525C9t2FFYJ9txWOfzi8v4hdSQi\ndaVwz3DuzsUjqz4RUQ+0FklvCvcM9/JbS6qM6V52kfSncM9wDz/9RqB+atTQkDoRkURSuGewv78W\nXOXxh70OpFG2rrGLRIHCPUOtWLWeSc8H12e/achJIXUjIommcM9AnxRs4H9HPRsY+/Fxh+p+dpEI\nUbhnmE2bt/Hze6dWGR8y6KgQuhGR+qJwzzDX/O7JKmNT9QxUkcjR1bMMU1xSGqh126NINOnMPYNU\nXhTsiXsuD6kTEalvCvcMUnlRsCaNs0PqRETqm8I9Q3y0Ym2gPq73wSF1IiLJoHDPELfc/3ygvu6C\n40PqRESSQeGeAT5dszFQd+7QJqRORCRZFO4ZYMQfnwnUfxwxKKRORCRZFO4RV1Lp1sfOHdrQsKEe\nwCESdQr3iDv3pkcC9Z03nBFSJyKSTAr3CFu3aUuVMd3+KJIZFO4Rdv2dUwL1o3dcHFInIpJsCveI\n+nrrjipLDbTM2SekbkQk2RTuEXXpLY8F6vt/dW5InYhIGOIKdzMbYGZLzWy5mY2sYZ/jzGy+mS0y\ns9cT26bsiaLikipjefu3CqETEQlLratCmlkWMAY4GSgA5pjZC+6+uMI+LYGxwAB3/8zM2tZXw1K7\n+ye/Gqj/fMt5IXUiImGJ58y9N7Dc3Ve4exEwBah8P935wDR3/wzA3dcltk3ZE7MXfhKoO7RtGVIn\nIhKWeMK9A7CqQl0QG6uoG9DKzGaZ2Twzq/a2DDO70szmmtnc9evX713HslvvLf4sUF92Vt+QOhGR\nMCXqgmpD4AjgVKA/cJuZdau8k7uPd/d8d8/Pzc1N0EvLN0pLy7jz4RmBsf59e4TUjYiEKZ4nMa0G\nOlao82JjFRUAG919G7DNzN4ADgOWJaRLictN9wTXkNm/9b5aakAkQ8Vz5j4H6GpmXcysEXAe8EKl\nff4G/MjMGprZPkAfYEliW5Xd+fjTL1i19svA2Nhfnx9SNyIStlrP3N29xMyGAzOBLGCiuy8ys2Gx\n7ePcfYmZvQgsBMqACe7+YX02Lt9yd0be91xg7JdDB4TUjYikgrgekO3uM4AZlcbGVarvBe5NXGsS\nrxvvfrrKWO+enZPfiIikDH1CNc25e5XpmGfvHxZSNyKSKhTuaW7KP+cG6rv1IA4RQeGe9qbOnBeo\nux6wf0idiEgqUbinsWUrvwjUJ/2we0idiEiqUbinsZv/FLxD5urzjg2pExFJNQr3NFXwRfAiasd2\n+4XUiYikIoV7mvr5PVMD9T03nRVSJyKSihTuaaryU5YaZcf1kQURyRAK9zS04cutgVq3P4pIZQr3\nNHTV7Y8Hat3+KCKVKdzTzPYdRYG66wF66JWIVKVwTzP//E9wPbabrxgYUiciksoU7mnE3Xly+ruB\nsRY5TUPqRkRSmcI9jZzzs4cD9cVn/DCkTkQk1Snc08S7H6ysMnbqMYckvxERSQsK9zTxxwkvBupJ\nfxiiR+iJSI0U7mmguDj4gaUjehxATrMmIXUjIulA4Z4GJj3/VqAeeUX/kDoRkXShcE8DL/5nUaBu\n0EB/bCKye0qJFOfugfqys/qG1ImIpBOFe4r7pGBDoD5Fd8iISBwU7inuf0c9G6jNLKRORCSdKNxT\nWOUpGRGReCncU9iQX00K1JPvvjScRkQk7SjcU9jW7YWBulnTxiF1IiLpRuGeos6+YVygvuP6M0Lq\nRETSkcI9BVU31979wHYhdCIi6UrhnoIm/212oL7rRj1GT0T2jMI9Bb3w2oJA3a2zHqMnIntG4Z5i\nKk/JaM12EdkbcYW7mQ0ws6VmttzMRu5mvyPNrMTMzklci5nluVfmB+rTju0ZUiciks5qDXczywLG\nAAOBHsBgM+tRw35/BF5KdJOZ5Inp7wTqrCz9cCUiey6e5OgNLHf3Fe5eBEwBqrsv7zrgWWBdAvvL\naCcf1T3sFkQkTcUT7h2AVRXqgtjYLmbWARgEPJS41jLPB8tWB+orf3J0SJ2ISLpL1M/8o4FfunvZ\n7nYysyvNbK6ZzV2/fn2CXjo6bh/z90CtddtFZG81jGOf1UDHCnVebKyifGBKbMXCNsApZlbi7s9X\n3MndxwPjAfLz87UqVgWlpcF/F/P2bxVSJyISBfGE+xygq5l1oTzUzwPOr7iDu3f55mszmwRMrxzs\nsntr1m8O1KNv/mlInYhIFNQa7u5eYmbDgZlAFjDR3ReZ2bDY9nG7PYDE5Wd3/TVQa912EamLeM7c\ncfcZwIxKY9WGursPqXtbmaXyB5cU6yJSV7pilwIqf3DpMa3bLiJ1pHBPAa/OXhKotW67iNSVwj0F\nrN3w9a6vjz6ia4idiEhUKNxDtmNnUaAedNLhIXUiIlGicA/Z7WOmB+pO7XR/u4jUncI9RKWlZSz/\nLLgUj26BFJFEULiH6KcjxgfqW4edGlInIhI1CveQrNu0pcrY4d07VrOniMieU7iH5OrfPhGoH7x1\ncEidiEgUKdxD8OXX26uMtcttEUInIhJVCvcQDL1tcqB+4p7LQ+pERKJK4Z5kn67ZGKj3a9GMJo2z\nQ+pGRKJK4Z5kI/74TKAe/9sLQ+pERKJM4Z5EDzzxWqDObZWj+9pFpF4o3JNo1rtLA/UDt5wXUici\nEnUK9yRZ+snaQD3kzKPIzs4KqRsRiTqFe5L8anTwqYM/Pv7QkDoRkUygcE+CkpLSQK2HX4tIfVO4\nJ8EjU/8TqO/75U9C6kREMoXCPQleeTv4pKWsLP22i0j9UsrUs7NvCD5H/PoLTwipExHJJAr3elT5\nKUsAxx7ZLYRORCTTKNzr0Zvv/zdQj7nt/JA6EZFMo3CvRw9NeT1Qf6fNviF1IiKZRuFeTyp/GvXH\nx+m+dhFJHoV7PXn6xXmB+pIzfxhSJyKSiRTu9eSLjV/v+vroI7pqgTARSSqFez34278WBOrzT+sd\nUicikqkU7vVg8t/eDtRt98sJqRMRyVQK9wRbs+6rQH3WSYeH1ImIZDKFe4Jdd+eUQH3uwPyQOhGR\nTBZXuJvZADNbambLzWxkNdsvMLOFZvaBmb1lZoclvtXUV1hUXGWsYUOt2S4iyVdruJtZFjAGGAj0\nAAabWY9Ku30CHOvuPYHfA+MT3Wg6eKbS7Y9TR18VUicikuniOXPvDSx39xXuXgRMAc6ouIO7v+Xu\nX8bK2UBeYttMD8+9Oj9Q6/ZHEQlLPOHeAVhVoS6IjdXkcuCf1W0wsyvNbK6ZzV2/fn38XaaBkfdN\nC9Tn9D8ipE5ERBJ8QdXMjqc83H9Z3XZ3H+/u+e6en5ubm8iXDpW78/Gn6wJjg085MqRuRESgYRz7\nrAY6VqjzYmMBZnYoMAEY6O4bE9Neerj5T88F6t9cc1pInYiIlIvnzH0O0NXMuphZI+A84IWKO5hZ\nJ2AacJG7L0t8m6mt8ln7oQdn5CUHEUkhtZ65u3uJmQ0HZgJZwER3X2Rmw2LbxwG/BloDY2MXEUvc\nPSNu8P7LtLcC9fDzjw+pExGRb8UzLYO7zwBmVBobV+HrocDQxLaW+tyd6a8vDIwd3+fgkLoREfmW\nPqFaB5ffNjlQX3ZW35A6EREJUrjXweYtOwL1qcf2DKkTEZEghfteWrR8TaC+68ZBIXUiIlKVwn0v\nuDu/fiBwwxDdOu8fUjciIlUp3PfCuL++EaiPOvzAkDoREamewn0vvPL2kkA94pKTQupERKR6Cvc9\ndMe4fwTqmy49WQuEiUjKUbjvofeXrArUR/XSlIyIpB6F+x5Y8t/PA/WQM48KqRMRkd1TuO+BW//8\nt0B92nG6r11EUpPCPU5vvv/fQL1fi2aaaxeRlKVwj0NRcQn3TXo5MDbmtsEhdSMiUjuFexwG/3xC\noD7s4DwaZce15pqISCgU7rV4dfaSKmO3XX1qCJ2IiMRP4b4b7s7Yp14PjE2++1LNtYtIylO478bv\nxgY/sNQyZx+aNW0cUjciIvFTuNegqLiEhcsKAmMTfn9RSN2IiOwZhXsNKl9EvevGQZqOEZG0oXCv\nxtk3jKsypiV9RSSdKNwr+cPD/6wyNua280PoRERk7yncK1ixaj3zFn8aGLvvlz/hO232DakjEZG9\no0/iANt3FHHTPc+wbtOWwPhFp/8PB7RvHVJXIiJ7T+EOXHX742zfWVRl/MwTe4XQjYhI3WV0uLs7\nV//2yWqDferoq0LoSEQkMTI63CdOe5P1XwanYu684Uy+993vhNSRiEhiZOwF1R07i5jxxoeBsQM7\n5irYRSQSMvLMfcu2nQz51aTA2LkD8/npgPxwGhIRSbCMO3P/dM2mKsEOKNhFJFIy5sy9uLiU837+\nSLXbnho1NMndiIjUr8iHe0lJKdNf/4D/98LsardPvvtSPXhDRCInrlQzswHA/UAWMMHd76603WLb\nTwG2A0Pc/b0E97pHPl+/meF3PFXj9hY5TZl4xyVJ7EhEJHlqDXczywLGACcDBcAcM3vB3RdX2G0g\n0DX2qw/wUOy/SbNjZxH//PciPlqxtsoSApVN+P3FtNp3nyR1JiKSfPGcufcGlrv7CgAzmwKcAVQM\n9zOAye7uwGwza2lm7dz980Q3/Ob7/931sOoGDWLXg90pc6/1e0855hCGnHkUWVkZdx1ZRDJMPOHe\nAVhVoS6g6ll5dft0ABIa7kXFJYx5ctauuqysLK7vG/ebC8jdLyeRrYiIpLSkXkk0syuBKwE6deq0\nx9+/Y2cxhUXFte53fJ+DOTa/Gz27ddjj1xARiYJ4wn010LFCnRcb29N9cPfxwHiA/Pz82udRKmna\nJJuh5/yITV9to6i4lItOD/4A0aCBfTtVIyKSweIJ9zlAVzPrQnlgnwdUfnrFC8Dw2Hx8H2Bzfcy3\nN8puyMCjD0n0YUVEIqfWcHf3EjMbDsyk/FbIie6+yMyGxbaPA2ZQfhvkcspvhby0/loWEZHaxDXn\n7u4zKA/wimPjKnztwLWJbU1ERPaWJqhFRCJI4S4iEkEKdxGRCFK4i4hEkMJdRCSCzONYk6VeXths\nPbD7Fb72XBtgQ4KPGbaovSe9n9QWtfcD0XtPB7t7reuphLaQubvnJvqYZjbX3SP1SKWovSe9n9QW\ntfcD0XtPZjY3nv00LSMiEkEKdxGRCIpauI8Pu4F6ELX3pPeT2qL2fiB67ymu9xPaBVUREak/UTtz\nFxERIhLuZjbRzNaZ2Ydh95IIZtbRzF4zs8VmtsjMbgi7p7owsyZm9q6ZLYi9n9+G3VMimFmWmb1v\nZtPD7iURzGylmX1gZvPjvSMjlcUe9znVzD4ysyVm9sOwe6oLMzs49mfzza+vzexnNe4fhWkZMzsG\n2Er5c1zTfsF3M2sHtHP398wsB5gHnFnpoeRpw8wMaObuW80sG/gPcIO7zw65tToxsxFAPrCvu58W\ndj91ZWYrgXx3j8Q94Wb2GPBvd59gZo2Afdz9q7D7SgQzy6L8+Rp93L3azwtF4szd3d8ANoXdR6K4\n++fu/l7s6y3AEsqfSZuWvNzWWJkd+5XWZxVmlgecCkwIuxepysxaAMcAjwK4e1FUgj3mROC/NQU7\nRCTco8zMOgOHA++E20ndxKYw5gPrgJfdPa3fDzAa+AUQ31Pa04MDr5jZvNjzjtNZF2A98JfY1NkE\nM2sWdlMJdB7w1O52ULinMDNrDjwL/Mzdvw67n7pw91J370X583V7m1naTp+Z2WnAOnefF3YvCfaj\n2J/RQODa2HRnumoI/AB4yN0PB7YBI8NtKTFiU0ynA8/sbj+Fe4qKzU0/Czzh7tPC7idRYj8avwYM\nCLuXOugLnB6bo54CnGBmj4fbUt25++rYf9cBzwG9w+2oTgqAggo/IU6lPOyjYCDwnrt/sbudFO4p\nKHYB8lFgibvfF3Y/dWVmuWbWMvZ1U+Bk4KNwu9p77n6zu+e5e2fKfzz+l7tfGHJbdWJmzWIX74lN\nX/QD0vbuM3dfC6wys4NjQycCaXlDQjUGU8uUDIS4cFgimdlTwHFAGzMrAH7j7o+G21Wd9AUuAj6I\nzVMD/Cr2LNt01A54LHaFvwHwtLtH4vbBCNkfeK78vIKGwJPu/mK4LdXZdcATsWmMFcClIfdTZ7F/\neE8Grqp13yjcCikiIkGalhERiSCFu4hIBCncRUQiSOEuIhJBCncRkQhSuEvkmdkQM2tfoZ5gZj0S\ncNzOZnZ+XY8jUh8U7pIJhgC7wt3dhyZohc3OgMJdUpLCXdKWmV0YWyd+vpk9HFucbJKZfRhbl/xG\nMzuH8mV5n4jt19TMZplZfuwYW83s3tg686+YWe/Y9hVmdnpsn85m9m8zey/266hYC3cDR8eOe2Ps\n9e81szlmttDMrop9fzszeyO234dmdnQYv1+SWfQhJklLZtYduAc4y92LzWws8AXli1+dHNunpbt/\nZWazgJ8VvwpJAAAByElEQVS7+9zY+K7azBw4xd3/aWbPAc0oX8q3B/CYu/cys32AMnffaWZdgafc\nPd/Mjosd57TYca8E2rr7HWbWGHgT+AlwFtDE3e+MfUp3n9hSziL1JhLLD0hGOhE4ApgT+8h8U+BF\n4Ltm9gDwD+ClOI5TFPs+gA+Awtg/Fh9QPu0C5evPP2hmvYBSoFsNx+oHHBr7aQGgBdAVmANMjC0G\n97y7z6/h+0USRuEu6cooP7O+OTBodgvQHxgG/BS4rJbjFPu3P76WAYUA7l5mZt/8/biR8p8KDqN8\nKnPnbnq6zt1nVtlQvnzuqcAkM7vP3SfX0pdInWjOXdLVq8A5ZtYWwMz2M7MDgAbu/ixwK98u8boF\nyKnDa7UAPnf3MsoXdMuq4bgzgatjZ+iYWbfYaosHAF+4+yOUP7kpKkvPSgrTmbukJXdfbGa3Ai+Z\nWQOgGBhB+cqG35y0fHNWPwkYZ2Y7gL15SPJY4Fkzu5jyKZxtsfGFQKmZLYi9xv2UT+W8F1u2eT1w\nJuUrlv6vmRVT/qzfi/eiB5E9oguqIiIRpGkZEZEIUriLiESQwl1EJIIU7iIiEaRwFxGJIIW7iEgE\nKdxFRCJI4S4iEkH/H3FWxUfknTiJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xc505eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution goes here\n",
    "def experiment(samples=10, lambdaa = 2, iterations=1000):\n",
    "    L = []\n",
    "    \n",
    "    print('experiment starting')\n",
    "    for i in range(iterations):\n",
    "        xs = np.random.exponential(1.0/lambdaa, samples)\n",
    "        lambdaa_hat = 1.0 / np.mean(xs)\n",
    "        L.append(lambdaa_hat)\n",
    "        \n",
    "    std_error = RMSE(L, lambdaa)\n",
    "    print('std error', stderr)\n",
    "    \n",
    "    cdf = thinkstats2.Cdf(L, label='estimates')\n",
    "    print('90% confidence int =', cdf.Percentile(90))\n",
    "    \n",
    "    thinkplot.Cdf(cdf)\n",
    "    thinkplot.Config(xlabel='estimates')\n",
    "    \n",
    "    print('experiment ending')\n",
    "experiment()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "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": 26,
   "metadata": {
    "collapsed": false
   },
   "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": 32,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE of goals 1.73778594769\n",
      "mean error 0.00442\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADjZJREFUeJzt3X+snQddx/H3h5ZFmeDUXX7YdnYxDaRBBkuzTSEYxJFu\nEro//GMVARHSLFkFDESGLvyhiyHR+INsUhucQIAuBrfYmMIYaEIMzPQOYWMbg2Ygu3WznSgQSRiV\nr3+c0+Ts7nb3ue3pPWf7vl9Jc8/znOd5+r398b7Pfc6Pm6pCktTHM2Y9gCRpfRl+SWrG8EtSM4Zf\nkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNbJz1ACs5//zza+vWrbMeQ5KeMu66665Hq2phyLZzGf6t\nW7eyuLg46zEk6Skjyb8P3dZLPZLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Izhl+SmjH8ktTM\nXL5y9+lu7w0HpnKcG6/fPZXjSOrFM35JasYz/qcBv4OQtBae8UtSM4Zfkpox/JLUjOGXpGYMvyQ1\nY/glqRnDL0nNGH5JasbwS1Izhl+SmjH8ktSM4ZekZgy/JDVj+CWpGcMvSc0YfklqZlD4k+xM8kCS\nI0muW+H+1ye5O8k9ST6f5KKh+0qS1teq4U+yAbgJuALYDuxOsn3ZZt8AfrmqfgH4I2D/GvaVJK2j\nIWf8lwBHqurBqnoMuAXYNblBVX2+qv57vHgnsHnovpKk9TUk/JuAhyaWl8brTuUtwCfXum+SPUkW\nkyweP358wFiSpNMx1Qd3k7yKUfjfvdZ9q2p/Ve2oqh0LCwvTHEuSNGHjgG2OAlsmljeP1z1OkpcA\nHwSuqKr/Wsu+kqT1M+SM/zCwLcmFSc4BrgYOTm6Q5ALgVuANVfW1tewrSVpfq57xV9WJJHuB24EN\nwM1VdW+Sa8b37wPeC/wM8FdJAE6ML9usuO9Z+lwkSQMMudRDVR0CDi1bt2/i9luBtw7dV5I0O75y\nV5KaMfyS1Izhl6RmDL8kNWP4JakZwy9JzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnD\nL0nNGH5JasbwS1Izhl+SmjH8ktSM4ZekZgy/JDVj+CWpGcMvSc0YfklqZuOsB9B82nvDgakc58br\nd0/lOJKmxzN+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Izhl+SmjH8\nktSM4ZekZgaFP8nOJA8kOZLkuhXuf1GSLyT5QZJ3Lbvvm0nuSfKlJIvTGlySdHpWfVvmJBuAm4DL\ngSXgcJKDVXXfxGbfBt4GXHWKw7yqqh4902ElSWduyBn/JcCRqnqwqh4DbgF2TW5QVceq6jDww7Mw\noyRpioaEfxPw0MTy0njdUAV8JsldSfacaqMke5IsJlk8fvz4Gg4vSVqL9Xhw9xVV9VLgCuDaJK9c\naaOq2l9VO6pqx8LCwjqMJUk9DQn/UWDLxPLm8bpBquro+OMx4DZGl44kSTMy5GfuHga2JbmQUfCv\nBn5jyMGTnAs8o6q+N779GuAPT3fYWZjWz54Ff/6spPmwavir6kSSvcDtwAbg5qq6N8k14/v3JXk+\nsAg8B/hRkncA24HzgduSnPy9Pl5Vnzo7n4okaYghZ/xU1SHg0LJ1+yZuP8LoEtBy3wUuOpMBJUnT\n5St3JakZwy9JzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Izhl+S\nmjH8ktSM4ZekZgy/JDVj+CWpGcMvSc0YfklqxvBLUjOGX5KaMfyS1Izhl6RmDL8kNWP4JakZwy9J\nzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Izhl+SmhkU/iQ7kzyQ\n5EiS61a4/0VJvpDkB0netZZ9JUnra+NqGyTZANwEXA4sAYeTHKyq+yY2+zbwNuCq09hXDey94cBU\njnPj9bunchypsyFn/JcAR6rqwap6DLgF2DW5QVUdq6rDwA/Xuq8kaX0NCf8m4KGJ5aXxuiHOZF9J\n0lkwNw/uJtmTZDHJ4vHjx2c9jiQ9bQ0J/1Fgy8Ty5vG6IQbvW1X7q2pHVe1YWFgYeHhJ0loNCf9h\nYFuSC5OcA1wNHBx4/DPZV5J0Fqz6rJ6qOpFkL3A7sAG4uaruTXLN+P59SZ4PLALPAX6U5B3A9qr6\n7kr7nq1PRpK0ulXDD1BVh4BDy9btm7j9CKPLOIP2lSTNztw8uCtJWh+GX5KaMfyS1Izhl6RmDL8k\nNWP4JakZwy9JzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Izhl+S\nmjH8ktSM4ZekZgy/JDVj+CWpGcMvSc0YfklqxvBLUjOGX5KaMfyS1Izhl6RmDL8kNWP4JakZwy9J\nzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Myj8SXYmeSDJkSTXrXB/krx/fP/dSS6euO+bSe5J\n8qUki9McXpK0dhtX2yDJBuAm4HJgCTic5GBV3Tex2RXAtvGvS4EPjD+e9KqqenRqU0uSTtuQM/5L\ngCNV9WBVPQbcAuxats0u4CM1cidwXpIXTHlWSdIUrHrGD2wCHppYXuLxZ/On2mYT8DBQwGeS/B/w\n11W1//THXd3eGw5M5Tg3Xr97KseRpHkzJPxn6hVVdTTJc4E7kny1qj63fKMke4A9ABdccME6jKWn\nsml9gQe/yKufIZd6jgJbJpY3j9cN2qaqTn48BtzG6NLRE1TV/qraUVU7FhYWhk0vSVqzIeE/DGxL\ncmGSc4CrgYPLtjkIvHH87J7LgO9U1cNJzk3ybIAk5wKvAb4yxfklSWu06qWeqjqRZC9wO7ABuLmq\n7k1yzfj+fcAh4ErgCPB94M3j3Z8H3Jbk5O/18ar61NQ/C0nSYIOu8VfVIUZxn1y3b+J2AdeusN+D\nwEVnOKMkaYp85a4kNWP4JakZwy9JzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nN\nGH5JasbwS1Izhl+SmjH8ktSM4ZekZgy/JDVj+CWpGcMvSc0YfklqxvBLUjMbZz2ANC/23nBgKse5\n8frdUzmOdLZ4xi9JzRh+SWrG8EtSM4Zfkpox/JLUjOGXpGYMvyQ1Y/glqRnDL0nNGH5JasbwS1Iz\nvlePdJZM671/wPf/0XR5xi9JzRh+SWrG8EtSM4PCn2RnkgeSHEly3Qr3J8n7x/ffneTioftKktbX\nqg/uJtkA3ARcDiwBh5McrKr7Jja7Atg2/nUp8AHg0oH7ShrIHxajaRhyxn8JcKSqHqyqx4BbgF3L\nttkFfKRG7gTOS/KCgftKktbRkKdzbgIemlheYnRWv9o2mwbuK2lG/A6ip1TVk2+Q/Dqws6reOl5+\nA3BpVe2d2OYfgfdV1b+Mlz8LvBvYutq+E8fYA+wZL74QeODMPrUndT7w6Fk8/jQ443TM+4zzPh84\n47Sc7Rl/rqoWhmw45Iz/KLBlYnnzeN2QbZ45YF8Aqmo/sH/APGcsyWJV7ViP3+t0OeN0zPuM8z4f\nOOO0zNOMQ67xHwa2JbkwyTnA1cDBZdscBN44fnbPZcB3qurhgftKktbRqmf8VXUiyV7gdmADcHNV\n3ZvkmvH9+4BDwJXAEeD7wJufbN+z8plIkgYZ9F49VXWIUdwn1+2buF3AtUP3nQPrcknpDDnjdMz7\njPM+HzjjtMzNjKs+uCtJenrxLRskqZlW4Z/3t49IsiXJPye5L8m9Sd4+65lOJcmGJP82firv3Ely\nXpJPJPlqkvuT/OKsZ1ouye+O/56/kuRAkh+bg5luTnIsyVcm1v10kjuSfH388afmcMY/Gf9d353k\ntiTnzduME/e9M0klOX8Ws0Gj8E+8fcQVwHZgd5Lts53qCU4A76yq7cBlwLVzOONJbwfun/UQT+Iv\ngU9V1YuAi5izWZNsAt4G7KiqFzN68sPVs50KgA8BO5etuw74bFVtAz47Xp6lD/HEGe8AXlxVLwG+\nBrxnvYda5kM8cUaSbAFeA3xrvQea1Cb8PAXePqKqHq6qL45vf49RrDbNdqonSrIZ+DXgg7OeZSVJ\nfhJ4JfA3AFX1WFX9z2ynWtFG4MeTbASeBfzHjOehqj4HfHvZ6l3Ah8e3Pwxcta5DLbPSjFX16ao6\nMV68k9FrhmbmFH+OAH8O/B4w0wdXO4X/VG8rMZeSbAVeBvzrbCdZ0V8w+sf7o1kPcgoXAseBvx1f\njvpgknNnPdSkqjoK/CmjM7+HGb325dOzneqUnjd+XQ7AI8DzZjnMAL8NfHLWQyyXZBdwtKq+POtZ\nOoX/KSPJTwB/D7yjqr4763kmJXktcKyq7pr1LE9iI3Ax8IGqehnwv8z+8sTjjK+T72L0RepngXOT\n/OZsp1rd+Knbc/tUwCR/wOiS6cdmPcukJM8Cfh9476xngV7hH/LWEzOX5JmMov+xqrp11vOs4OXA\n65J8k9Hlsl9J8tHZjvQES8BSVZ38bukTjL4QzJNfBb5RVcer6ofArcAvzXimU/nP8bvtMv54bMbz\nrCjJbwGvBV5f8/c89Z9n9EX+y+P/O5uBLyZ5/iyG6RT+uX/7iCRhdF36/qr6s1nPs5Kqek9Vba6q\nrYz+DP+pqubqTLWqHgEeSvLC8apXA/P2MyC+BVyW5Fnjv/dXM2cPQE84CLxpfPtNwD/McJYVJdnJ\n6PLj66rq+7OeZ7mquqeqnltVW8f/d5aAi8f/Vtddm/CPH/g5+fYR9wN/N4dvH/Fy4A2MzqK/NP51\n5ayHeor6HeBjSe4GXgr88YzneZzxdyOfAL4I3MPo/+LMX9mZ5ADwBeCFSZaSvAV4H3B5kq8z+k7l\nfXM4443As4E7xv9v9j3pQWYz49zwlbuS1EybM35J0ojhl6RmDL8kNWP4JakZwy9JzRh+SWrG8EtS\nM4Zfkpr5f+R+HilGhUhhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xca6cb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbc47860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def soccer_sim(lam=3, games=50000):\n",
    "    \"\"\"games is how many games we're simulating\"\"\"\n",
    "    \n",
    "    estimates = []\n",
    "    for i in range(games):\n",
    "        goals = SimulateGame(lam)\n",
    "        estimates.append(goals)\n",
    "        \n",
    "    print('RMSE of goals', RMSE(estimates, lam))\n",
    "    print('mean error', MeanError(estimates, lam))\n",
    "    \n",
    "    pmf = thinkstats2.Pmf(estimates)\n",
    "    thinkplot.Hist(pmf)\n",
    "    thinkplot.Show()\n",
    "    \n",
    "soccer_sim()\n",
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}