{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap confidence intervals\n",
    "> A Summary of lecture \"Statistical Thinking in Python (Part 2)\", via datacamp\n",
    "\n",
    "- toc: true \n",
    "- badges: true\n",
    "- comments: true\n",
    "- author: Chanseok Kang\n",
    "- categories: [Python, Datacamp, Data Science, Statistics]\n",
    "- image: images/bootstrap-reg.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ecdf(data):\n",
    "    \"\"\"Compute ECDF for a one-dimensional array of measurements.\"\"\"\n",
    "    # Number of data points: n\n",
    "    n = len(data)\n",
    "\n",
    "    # x-data for the ECDF: x\n",
    "    x = np.sort(data)\n",
    "\n",
    "    # y-data for the ECDF: y\n",
    "    y = np.arange(1, n + 1) / n\n",
    "\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating bootstrap replicates\n",
    "### Bootstrapping\n",
    "- The use of resampled data to perform statistical inference\n",
    "- Bootstrap sample : A resampled array of the data\n",
    "- Bootstrap replicate : A statistic computed from a resampled array\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing bootstrap samples\n",
    "In this exercise, you will generate bootstrap samples from the set of annual rainfall data measured at the Sheffield Weather Station in the UK from 1883 to 2015. The data are stored in the NumPy array ```rainfall``` in units of millimeters (mm). By graphically displaying the bootstrap samples with an ECDF, you can get a feel for how bootstrap sampling allows probabilistic descriptions of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_fwf('./dataset/sheffield_weather_station.csv', skiprows=8)\n",
    "rainfall = df['rain'].astype('float')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de1zUdd7//8fMcBwBFeSggKgo4AE8l2mZtiWlmFl2srK2vdy6ts1f7nVrazevba92223brpsddm9d1bdru7XpVXZYzQ6K1ZpbWh7ylAdQEBSEYWAMGBhg+Mzn98fAhxkYEMUZZobX/R/nc+Az77ej8+T9eX/e77dOVVUVIYQQA56+vwsghBDCP0ggCCGEACQQhBBCtJFAEEIIAUggCCGEaCOBIIQQAvBBIFitVvLy8igrK+ty7NixY9x8883k5ubyxBNP0Nra6u3iCCGE6EaINy9+8OBB1qxZQ0lJicfjjz76KL///e+ZMmUKv/71r9mwYQPLly/v9fXPnWvA4XAfRhEXF0VNjbUvxfZrwV4/CP46Bnv9IPjrGKj10+t1DB06qNvjXg2EDRs28OSTT/LLX/6yy7Hy8nKampqYMmUKADfffDMvvvjiBQWCw6F2CYT2/cEs2OsHwV/HYKtfbW0tJtNZTp0qwm5v4cSZc1iawkkeGsqD993e5fz/27iF4ooGxgwfxJ03Xd+n9/7f//sHJaYmRiVGcP+dS/t0rQvhz5/hG29vpKSqiVEJEdx3x029/jmvBsLTTz/d7bGqqiri4+O17fj4eEwmkzeLI8SA0dLSQkVFGWazmUGDjIwYMZLBgwcDUFBwlMrKCsrKztDc3ERNgwF7aBzpw6O4fckC/v7uZk6ZbIxOjOSeWxdr11z3wScUVzQyZriRu25eqO2vra3l+f/9B2Z7HDoc6HV66pVRAJRXqlS/vJF5cy7HoaooioPdB45SWBWCyhCKa1VK/l8+ORMzcKigqiqo4FBVnC/b/mw7prYdA+efRwqKqKiPAWIoK4WTf/mU8ePGoH1Vt52rum9qe9q3O+1GRe203UFVISIihKYmezc/236ee2B0fS+107an9/J8jc7HXXefOWvinC0aiOZMiQpvb+x1KHg1EHricDjQ6XTatqqqbtu9ERcX5XF/fHx0n8rm74K9fhD8dexr/RRFQVEUbDYbxcXFNDQ0MHLkSCIjIykqKmLr1q1YW41U24eiA9KGHuSRB5dTWlrK25/s5lzrEFR1MDZHEq2EAlD0g8rBU+9jahqMipFSq8rZ/93KrJmT2bnnMCeqwoFwTtXDiVe2MT5zDIqisnPfMRodqW0l6/x/WEdxbQzFnxxz2dfxtaOi40S1nhNfFnepo04HOp0Ovc55Hb0OdHodOjr2Nza1X0sHqFRZQ7CdMLft0bkVSedyXdc9uk7H23d02e96rNO1Lvi9zlMGT9+F3b5Xp59taG6PBx0qcLq6udf/3votEJKSkjCbzdp2dXU1CQkJF3SNmhprl2ZbfHw0ZnP9JSmjPwr2+kHw1/FC6mez2WhpaSIsLILIyEhKS4upqnK2pKuqKiktPYW11Ui9EkW0YRc2JQKzPY4Wx3haCdOuY66Cx1/8jB/qbFi7+fJW0VHZNNRtu6BKT8HHRwEDzt9DnV+8Zed0/HC4AoNeh80R4eF6HecONZzj8ZWLMOh1GPQ6Nny4lW9OG1HRoUNlTloDK25ztkT0Op0WBL3xP2+8zZ7KYW3vpjIjqYYH77ujVz/bF/78b/SNtzfyr5Io7e9k5LBwrax6va7bX6ShHwMhOTmZ8PBw9u3bx/Tp09m0aRNz587tr+II0W8URaGkpIhz52qIj08kLi6B48ePUllZhsViwWqtA8BojKaqTqFeiSKEVn5ojcaqTEBp+w2/K/cv6PIf9MAgl2Oqy5/OL4+kMBOVLYnal3VW9Gn+v5/ey1vvfcjXpdHal8yctHp+fOe1APzX2jcobU7F/cYFbdsqGYkO4odEantvW3Q1Te9uoqreQEK0wi0LlxBiuLgHHu9cegMNb73HDy3hDAlr5s6lyy7qOsHkvjtuAn/sQ/Bk5cqVrFq1iuzsbJ577jnWrFmD1Wpl4sSJrFixwtfFEcJnFEWhutpMba0JRQklLCyMb77ZQWHhcbfzQkMjOGfTU20fSpMyjBbHCHQ6lehGKzWtcagebsv0zPlFnT30DLWNuk5f3s4/o/RWUiIqGD86gWOnTra1OKwsnDebsFADtyy8hrq3/0FNYwhxxlZuXtjRebs8bxbrP/oGsz0OPQ4iDU0Y9TZaCSElLowHVtzmVprBgwez4tYlNDZaMRqjtL6NizF48GD+7e5ll+RaweRCQsCVLpCnv5ZbRsEpkOuoKAoOh4Jeb8Bms3HunJn6eitHjx6mqqpSOy8ychA2WwPm5ljONifSioEQnUJcqAVFNWC2x3u4eu+++D1JCz9D3twJxMTEsP6jb9r6EFRa1HCGhvzA1JEORo0ay8yZszhy5BAVFeUMH57MxIk52jVqa2u7/eI9ceI4Z86UMnjwUDIzJzB69PCA/Qx7I1D/jZ7vlpEEQoAJ9vpB4NWx/T5/dbWZgoJj2O0t1NfXU1f3g3ZOmS2JGvtQdDgI1Tu/+BuVCKpbXb/4e/uFr6O7L/+kUBN1SjStDj0KBiL1zaREVDD/iinMnDkLcH55m0yVREYaiY2NJSpqsNsTf5dCoH2GFypQ6+e3fQhCBCpFUTh3zoLN1khrayslJaewWMyYTBXaOdZWIxXNaTQ7wmlVDW4dvC0OaGju7j/l+b70O271ROmtNDnCUTEQqbeRElHJqKRBVFefAGDEiBRGjx7HpEl5GAwG7QrjxmUxblzWRddfBC8JBCF6SVEUmpps7Nz5JSdPnkBVHW7HixtSOacMRgc4MHT66c6//aud9qsuZ7Y/Ba8SqbOhqAb0OpXEMDORhibt/v6QCDutrXbi4uLJyprIpEl5KIpCa2sLISFhhIWFIcSFkEAQohNFUbDbWwBobm6hqcnZEjh3rppDhw7xww81bS2AeJod4ehQaVKdv6kDHjp9Xandvg6hleTwSoYPi+LsuVaM/EBUSCMAoaFhXHbZbJKSRlBX9wNhYWEMGhRFRISRqKiO1obBYJAgEBdNAkEIF2azmcrKMlpbFQBOnz5FVVUldrtdO+dwXQbNRHb6ye5CoOttnzBdC8PDTEQamrSBY3Gh54gKaeS66xYyfHgKTU2N2O12WlqaMBhCiItLIDLS+Z6JiUmXoKZCdCWBIESboqJCvvnmK6xWK4rSMfOuuTkWsz0OVQWbGo5zkFbvWwGROhsGnYKKnmGhFtLjW8nMnMCQIcMoKTlBfX09iYnjmTRpivb0jutv/UL4igSCGPAUReH48SN8+eVn2j5zcyyVLfE0q6HQY3+A547fSJ2NqJBG7Tf/8PBwIiIiGDYsgYULr8dmUzEYDIwaNQq93uDW6StEf5FAEAOWzWajsrKMoqIT2uCwMlsSlfZhuIdAz7eDwmhCr1MJ0SlEGpqICz1HZlo8SUkpxMZOZsSIFPR6A6qqEBISRmxsrPbIogSB8CcSCGLAKS0tpqLiLCZTBeXlZ7T9B+uysBPetuUpBLq2BpJCTaREVhISEkJMzBDi4xOZMOFKhg9P9k7hhfAiCQQxYCiKwu7d37B//7favvb+gUZHe98AeHoUtF0YLRh0reh1KsNCLWSnRTB9+lKGDYsnIiJSfuMXAU0CQQS19lHEiqJy4MAejh8/oh07Wj+WRtV19ShPfQMKBmgb+FWhPQaq1xsYPjyZH/3oeukAFkFDAkEELZvNxoED31FXV0tx8QlUVcHaauRE4ygUQnAGQHcDxiAcG9kxhdp2SspIJk2aTGhoCDpdCEOGDJUwEEFFAkEEnfZZRY8dO8zRo4e1/cUNqViU2Lat7oMAVJJCq0iJrCQ2dhhDhgwhJSWN5OSRxMTEyG0hEbQkEERQaWlpYfv2bZw8WQCgrRpWYx+qjSTuvo9AJVpfT2bUKQByc/NITnYuJhMaGiZBIIKeBIIIGmazmX/9659UVpY5t5tj2+b+b+epj0BFj4OE0GpSIp3TU+v1YSxdukxGBIsBRwJBBIWKinK2bduC1VqLtdVIiS2ZJtXYdrS7gWQKM2K+d7vOuHFZXHHFXOkbEAOSBIIIaIqiUFlZzhdf5GO11vXQKnDvI3C9NQQwYkQql18+h4SERLk1JAYsCQQRsGw2G8eOHWb37l04HEqnx0g93x7qHATz5+cSHR0jTwwJgQSCCFBWq5Vdu7Zz4oTzsVDPo4w7WgV67EyLOep2jdzcPNLTM3xQWiECgwSCCDgnThxn586vaGioA1zDwHOroP0RUoCwsAhSUlLJzp5McvJIn5ZbCH8ngSACRkVFOYcP7+fkyY7BYvvrxqMQhqcw6NwqmDNnHomJIxg8eLC2toAQooMEgggIpaXFfPbZpzQ3N2v7DtZltYWBK2cYuI4yTksbw5Qp06RFIMR5SCAIv2cyVfL559u0MOg6RbV7n4FrGEyePJ2ZM6+QZSWF6AUJBOHXystPs3nzRhwO5wpm53uSyDUMZs6czcyZs3xYWiECmwSC8FtnzpSwefMH2vZ3dRNwENq21XMYXHfdQsaNy/JRSYUIDhIIwi/t27ebb7/9CmifgmIEXdcy7nisNNZgYcwg52I3l112JWPGjPNdYYUIEhIIwu9s2LCBY8eOAVBgHU29I6btSNfxBToUMo3F2joFV199LVlZE2W0sRAXQQJB+A1FUdix459aGByuy6CZ9sdDu05X3XnU8dVXX8vEiTm+KawQQUgCQfiFlpYWPv98C6dOnQS66y/oaBmkhZcRH24BIDw8klmzriIra7wPSyxE8JFAEP3OZrOxadP7WCxVAOytm0T3/QUO0sLLtTBISRnF5ZfPYdiwYXKbSIg+kkAQ/cpsNvPxx/+gsdEKuIaBK2cYhNLM5Jjj2t6RI0dx1VU/YvDgwT4qrRDBTQJB9IuWlhYKC4+xY8fn2r6uLYOOW0TR+jq3/oKJE6cwZcp0CQMhLiEJBOFzNpuNt9/+OzabVdvXfRionfoLIpg9+2oyMrLkFpEQl5jemxffvHkzCxcuZMGCBaxbt67L8SNHjnDLLbdw44038sADD1BXV+fN4gg/YLVaefPN188TBu0UZsQc0sIgPj6RZcvuYvx4eaxUCG/wWiCYTCbWrl3L+vXr2bhxI++88w4nT550O+fpp59m1apVfPjhh4wePZrXX3/dW8URfsBms/HZZ1tQlBZtX08dyFnGYm1Pevo4liy5VW4RCeFFXguEnTt3MmvWLIYMGYLRaCQ3N5ctW7a4neNwOGhoaACcXxYRERHeKo7wA/v2fcvZs6cB54CzvXXZdL1N5AyDpFCTNths6tTLyM1dLBPUCeFlXutDqKqqIj4+XttOSEjg0KFDbuc8/vjj3H///fzhD38gMjKSDRs2eKs4op8dOXKIQ4e+A1wnqPO8oI1rn4FMUCeE73gtEBwOBzpdx394VVXdtpuamnjiiSd44403yMnJ4W9/+xuPPfYYr776aq/fIy7O8xq48fHRF1/wABBo9Ttw4ABffvkZANZWo8tspe3aw0BhRsz32t6JEyeycOF1vimkjwXaZ3gxgr2OwVg/rwVCUlISe/fu1bbNZjMJCQnadmFhIeHh4eTkOKcauP3223nhhRcu6D1qaqw4HKrbvvj4aMzm+j6U3L8FWv327PmGPXt2atvHG9PbXrk/TaSjlekxR7TzoqKiufLKawOqrr0VaJ/hxQj2OgZq/fR6Xbe/SIMX+xBmz57Nrl27sFgs2Gw28vPzmTt3rnY8LS2NyspKioudHYeff/452dnZ3iqO8DGbzcY33+x0C4P9deNx/pNzDwMDLW5hEB0dzZIlt8mTREL4mNdaCImJiaxevZoVK1Zgt9tZtmwZOTk5rFy5klWrVpGdnc0f//hHHnnkEVRVJS4ujj/84Q/eKo7wIZvNxo4dn1NU1LH28eG6DA/LXQIoTI05pm2NGJHKnXfehs2mejhXCOFNOlVVA/Z/ntwy8j9ms5n9+3dz8mSBtu9gXRZ2wtu23FsHaeFntA7k2Nhh3HHHCr+vY18Fe/0g+OsYqPU73y0jGaksLpny8tNs3/45tbXntH0dYdD1iSKjrkELg2HDElm8+GYfllYI0ZkEgrgkbDYbX3yxjfr6Wm3f4boMl5ZBu46J6iZEOwcqZmVN4oorriIyMhIhRP+RQBCXxJdfftYlDLoubtN11tKEhCSuuWaBD0sqhOiOBILos1OnTlJcfELbLm5I7TYMwrGRHdPR2Xz11cE5zkCIQOTVye3EwPDdd3vdti3K0LZXXVsGrmEwefJ0t9HsQoj+JYEg+sRsNmMyndW2nZPV6fA01sB1cZvU1JFMnjzddwUVQpyX3DISfXLkyAHtdU/TWLuONcjIyGLWrLlERXX/+JsQwvckEMRFs1gsHD16GHB2Ine32pnrNNYpKWkSBkL4KQkEcdG++2434JywrqMT2ZWDLGORNo31uHFZXHGFhIEQ/koCQVwUi8VCYeFRAE40jmrb6946cA0DmcZaCP8nncriovzzn/naa8Xt94qOJ4rawyA1NU3CQIgAIIEgLlhRUaH2ZNHBuizcnyoCULUnikJCwrn22oU+L6MQ4sJJIIgLtnXrx9pru9sMpu2tg/Y1k3XcfvvdMiWFEAFCAkFckF27vqL9i9+5vkH3rYM77riXwYMH+7qIQoiLJIEgeq22tpb9+51PFpmbYz2ub2CgFYD09ExiY2N9Wj4hRN9IIIhe+/zzLdrr0uaUtlfuTxaNM5YAMG3aZb4rmBDikpBAEL1WWVkOwNH6sXianiIcG1EhjYwZM07mKBIiAEkgiF757LNPAecgtEZ1kIczFG3iussum+PDkgkhLhUJBNErhYXOuYiON6a37XFvHSSFVgMQFhYhfQdCBCgJBHFeVqsVcK5z4Pwn4x4GeuykRFYCMH78JN8XUAhxSUggiPN6883XANd1DlwpTIs5qm1JIAgRuCQQRI8+/ngjoFJmS8JTR3KsoWPZzPnzF8jtIiECmASC6FZRUSGlpc6pqyvtCR7OUBgz6AwAaWljpHUgRICTQBAeKYrC1q0fdXO0a+vgyivn+6BUQghvkkAQHm3YsE573XXcAYCqtQ4mT54hU1QIEQQkEEQXmza9x7lzzsdIy2xJ3Yw7cABgMITKrSIhgoQEgnCzZ883lJefBpyD0CrtiW1HPI87uPbaXOlIFiJISCAITW1tLXv27NS2CxrHtL1yDwMDLaREVpKcnEZ6eoZvCymE8BpZQlMAziUx3377DW27wDoaFYPLGWrbnwpTY5yjlmfMmOmz8gkhvE9aCAKLxcI777ypbZubY6l3xLRtuXYkO5gR8z0AISFhJCeP9F0hhRBeJ4Eg+Oc/81FVh7bd3dTW0Xqrdk5GRqaPSieE8BUJhAHOZKrU1kcG2Fs3CU8jknW0khl1yvlaZ2DixCm+LagQwuskEAa4L77I115/VzcBMNA5DEBheswR7bxrr82V9Q6ECEJeDYTNmzezcOFCFixYwLp167ocLy4u5p577uHGG2/kJz/5CbW1tR6uIrzlxInj2niD7+om4CAU9z4DAEXrNwC47rqFjBuX5btCCiF8xmuBYDKZWLt2LevXr2fjxo288847nDx5Ujuuqir//u//zsqVK/nwww8ZP348r776qreKIzqxWCxs2/YJ4ByJ7AwDV87WQVp4x+2ktLR0CQMhgpjXAmHnzp3MmjWLIUOGYDQayc3NZcuWjjV5jxw5gtFoZO7cuQA8+OCD3HXXXd4qjuikqKhAe92oGttede5EriM+3KKdN2XKVB+VTgjRH7wWCFVVVW73mRMSEjCZTNr26dOnGTZsGL/+9a9ZunQpTz75JEaj0dOlhBfs2bMLgMN1GXjqRI7W12mdyAATJmSTlJTs20IKIXzKawPTHA4HOl3H/WhVVd22W1tb2b17N2+99RbZ2dk8//zzPPPMMzzzzDO9fo+4uCiP++Pjoy++4AGgr/U7darji76ZCJcjHSuguYbB2LFjycu7gcjIyD6974WQzzDwBXsdg7F+XguEpKQk9u7dq22bzWYSEjrm1I+PjyctLY3s7GwA8vLyWLVq1QW9R02NFYdDddsXHx+N2Vzfh5L7t0tRv7fffhdof6qo6yymriugRUXFMGPGlVitrVitvvl7lc8w8AV7HQO1fnq9rttfpMGLt4xmz57Nrl27sFgs2Gw28vPztf4CgKlTp2KxWDh+/DgAX3zxBRMnTvRWcUQbRVFoabFhbo7t1JHsDNZQWrQ9oaHhzJ9/rUxeJ8QA4bUWQmJiIqtXr2bFihXY7XaWLVtGTk4OK1euZNWqVWRnZ/PXv/6VNWvWYLPZSEpK4tlnn/VWcUSb9r6D083t/QHuU1NMjnEG9KBB0cyb9yNSU0f5tHxCiP6jU1VVPf9p/kluGV2Y2tpa1q17HYC9dTl03C7q2pG8ePHN/RYG8hkGvmCvY6DWr99uGQn/s2vXl0D3K6C1h8GYMeOkZSDEACSBMEDU1tZSXOwcGNgx7qCDHkV7PXWqTGstxEDUYyAsXbpUe71jxw6vF0Z4h6IobNrkfLLI3BxL55YBQIbR2TqYMGEKiYlJPi6hEMIf9BgIrt0La9eu9XphhHecPn0Kq7UO8DS1NYCDqJBGAKZOne7bwgkh/EaPgdB5YJkITJ9++iHQ06hk5zoHU6dexuDBg31fQCGEX+h1H4JrOIjA0T69tbXVSDOeRhorWmdyZuYEH5ZMCOFvehyHUFdXx7Zt21BVlfr6evLz892OL1iwwKuFE31TVFTI8ePOqasLGse07XVvHbTPZpqcnCoD0IQY4HoMhBEjRvDmm861docPH87f//537ZhOp5NA8HP5+c7ZZa2tRlQMLkfaRyU3a7OZTpkyzdfFE0L4mR4DwTUARGCx2WyoaivgqXUAoGqjksPCIkhIGOHbAgoh/M55p65oaGjgo48+orCwkIiICDIzM7n++usJCwvzRfnERdqx43MAymxJHlsH4TRpe6666mqfzmQqhPBPPXYql5aWsmjRIvLz8wkPDwfgvffe4/rrr6e8vNwnBRQXp6ioEGurkUp7Ytse99ZBdkwhAEOHDmPkyDFdfl4IMfD02EJ48cUXWb16NUuWLHHb/+677/Lcc8/J2AQ/ZTJVAnCicVTbHveOZKOuUTt3woRJ0joQQgDnaSEUFhZ2CQOAW2+91W2RFeFf3n9/PQCKW963jyNRmBDtnMIiOnowo0al+7ZwQgi/1WMgGAyGbo/JuAT/tGvXV4Cz78DTBHYzYr7XtrKyJspANCGEptcjlUVg2L9/NwCV9ngPRztGm+v1BtLTM3xUKiFEIOixD6GyspLf//73Ho+ZTCavFEhcPIvF4rLVdQK79ikqADIzJ8pANCGEmx4D4a677ur22PLlyy95YUTf5Od/AnS/VnL7FBUA48aN9WnZhBD+r8dA+PnPf95lX0tLi4xB8FMWSxVltqROayU7ua53MHbsOBmIJoToosc+hJaWFh577DG2bdum7Xv44Yf51a9+RWtrq9cLJ3qv/XZRpT2hbY/7o6bt6x0YjdFMnjxTQl0I0UWPgfDiiy9itVqZNq1jnpunnnqK2tpaXnrpJa8XTvTe7t1fedjbsVZyx3oH02QBHCGERz0Gwvbt2/nv//5v4uLitH2JiYk8++yzfPbZZ14vnOi94uKT510rOTp6COnpmf1RPCFEAOgxEEJDQ4mIiOiyPyoqSm45+JFTp3q3VnJW1gSioqJ8Vi4hRGDpMRD0ej1Wq7XLfqvVKn0IfqK2tpbdu7/xcMS970Cn08u4AyFEj3oMhLy8PNasWUNjY8fcN42NjaxZs0bWQvADZrOZr776gpqaqm5HJrf3HYwenS7jDoQQPeoxEO69916io6OZM2cOt912G8uWLWPOnDnExMTw0EMP+aqMwoOWlha++243paXOFkDH00WuHNqrxMThPiqZECJQ9TgOQa/X87vf/Y4HH3yQI0eOoNfrycnJISHB05eP8KWCgqMUFRV0c9R5uyjWUAtAWFg4KSlpPiqZECJQ9RgIZ8+eZcSIESQnJ5OcnOx2bMeOHcydO9erhROemc1m9u379jxnqYwZdAaAcePGEx/vaW4jIYTo0OMtI9fbQg8//LDbMVkLof8cP36YxsYGbftwXQZd+w+cwsONTJiQ7bvCCSECVo+BoKods2OeOXOm22PCdyorKzl8+IDbvma6PhpswPkUWE7OFGkdCCF6pdfTX3eeClumxu4fe/bs6eaI+1QV44wlAGRkjPd6mYQQwaHXLQTR/2pra/n+++/d9u2tm0TXW0XOx01HjkyTBXCEEL3WY6eyw+GgtrYWVVVRFEV7DaAoSk8/KrygrOw0LS0t2rYzDAx0bh2E0wTAoEHRvi2gECKg9RgIhYWFzJo1SwuByy+/XDsmt4x8r7S0SHvtnLfINQzaKWTHFAIQGhrps7IJIQJfj4Fw/PjxPl188+bNvPzyy7S2tnLvvfd2u+DO9u3beeqpp/jiiy/69H7BzGw2U1JSDIC11UijOqjTGc7QzjIWa3uSk2VWUyFE7/UYCH1hMplYu3YtH3zwAWFhYdxxxx1cfvnljB3rvlJXdXU1f/rTn7xVjKBRVlaivT7emN72yv1WUVKoSZuqYsiQOJKSUnxXQCFEwOuxU7kvdu7cyaxZsxgyZAhGo5Hc3Fy2bNnS5bw1a9Z4XJlNuKuuNgPtt4r0dA4DAy2kRFZq5yclDScyUm4ZCSF6z2uBUFVV5fb8e0JCAiaTye2cN998kwkTJjB58mRvFSNoVFU5/+48TXENClNjjrntSU1N9UGphBDBxGu3jBwOh1vHs6qqbtuFhYXk5+fzxhtvUFlZ6ekS5xUX53lu//j44Hu6prb2nIcRyV37DcAZvjk5gb32QTB+hq6CvX4Q/HUMxvp5LRCSkpLYu3evtm02m90mxduyZQtms5lbbrkFu91OVVUVy5cvZ/369b1+j5oaKw6H+1iJ+PhozOb6vlfAj5jNzttFnkYk61C0foN2iYkp2GwqNltg/j0E42foKtjrB8Ffx0Ctn16v6/YXafDiLaPZs2eza9cuLC01jcMAABg4SURBVBYLNpuN/Px8t8nwVq1axdatW9m0aROvvvoqCQkJFxQGA8nOnds97HUGYWan1kF4eCQZGVneLpIQIgh5LRASExNZvXo1K1as4KabbiIvL4+cnBxWrlzJ4cOHvfW2Qcdms1Fefqbb9ZJdWwehoeHMnn0ViYnyuKkQ4sJ57ZYRwOLFi1m8eLHbvtdee63LeSkpKTIGoRtnz5YB518vGSAjI0vmLhJCXDSvtRDEpfH99wc87HVfLxkgKiqG7OypGAwGH5VMCBFsJBD8WG1tLeXlZyiwjuZ8t4vS0zNkzWQhRJ9IIPixfft2A1Dv8PR4W8d6yRERg+RWkRCizyQQ/Njx45463523i9LCz2p7pk6dJovgCCH6TALBT5044ZxYsLuni+LDLQDExSUybpy0DoQQfSeB4Ke2bfsEOP/TRTk5OQE9IlkI4T8kEPxQUZFzPYOOzuR27k8XxcUlMGrUWIQQ4lKQQPAzLS0tbN36EeDamewaCh1TVYwaNUZmNBVCXDISCH5mz55vALqdyC4ptBqAsLBwxowZ5/sCCiGClgSCnzl40DkhoKeJ7EDR1jwYMyZDniwSQlxSEgh+pKKiHMDDQLSu01xnZmb6tnBCiKAngeAnFEVh48Z3ge4GonX0HQwfnkJSUrIPSyeEGAgkEPzEoUP7UVUHB+uy8NQ6cB2IlpQ0QuYsEkJcchIIfsBms7Fr1w4A7IS5HHGGgR67NhANIDl5hC+LJ4QYICQQ/MDOndsBKG5Ixf0RUwCVaTFHta2xY8eSlJTio5IJIQYSCYR+Vl5+moKCYwBYlKFte11DweF2/uTJkwkLC0MIIS41CYR+tnevc0bTjtaBe99BrKFWOzc5OY309HTfFlAIMWBIIPSz8vLTgGvrwJXCmEFntK1p06bLyGQhhNdIIPSjlpYWAL6rm0BPo5IBRo4cTWrqKJ+WTwgxsEgg9KODB7/D3ByLg1CXvc4w0NGqjUoGSE+XaSqEEN4lgdBPrFYre/bs5HRz+wAz947k6TFHtK2hQ+MYMSLVp+UTQgw8Egj9ZOfO7QCoHqa3jtZb3c6dPfsqBg8e7JNyCSEGLgmEflBaWszJk4XdHFXJjDqlbaWnZ5CWNsY3BRNCDGgSCP3g4483AniYpqKradMu902hhBADngSCj7Wvhgaep6kIp0nbM3365TLFtRDCZyQQfGzPnl0A7K8bT9fWgUp2jDMwEhOHk5MzzeflE0IMXBIIPmQyVWKx1FBmS0Kh6/QTehQAdDoDs2fPlUFoQgifkkDwoX37nMtjVtoT2va4D0TLMDo7k+fP/xHDh8t6B0II35JA8JEzZ0ooKSnG2mqk820iAKOugaiQRgyGEMaMyeiXMgohBjYJBB8wmSrZvPkDAAobR7ftdQ+FCdEnARg7NkNmMxVC9AsJBC+zWCxs3fqxtu3AdaWz9tZBo7YnJ2e6r4omhBBuJBC8qKWlhT17vsZqdU5hfbguA09PFrW3DhITR8hjpkKIfiOB4EUVFeUUFZ0AwNwcSzORdB2E1rEAzpQp8pipEKL/eDUQNm/ezMKFC1mwYAHr1q3rcvyzzz5jyZIl3HjjjfzsZz+jtrbWw1UC1+HD32mvS5s7L3vpvF2UFn4WcM5mOmqULH4jhOg/XgsEk8nE2rVrWb9+PRs3buSdd97h5MmT2nGr1cpvf/tbXn31VT788EMyMzN56aWXvFUcn7NYLJw+XQp0v95BKM3Eh1vQ6fTMmjUXg8Hg8VpCCOELXguEnTt3MmvWLIYMGYLRaCQ3N5ctW7Zox+12O08++SSJiYkAZGZmUlFR4a3i+JTNZmPLlk2Ac2lM9/UO2ilMjjkOwNy518hspkKIfhfirQtXVVW5dZAmJCRw6NAhbXvo0KFcd911ADQ1NfHqq69yzz33XNB7xMVFedwfHx99ESW+dA4cKOKHH84BrktjurcO2m8VAcybN+eCrt/f9fOFYK9jsNcPgr+OwVg/rwWCw+FAp+voQFVV1W27XX19PQ899BBZWVksXbr0gt6jpsaKw6G67YuPj8Zsrr+4Ql8iX3+9E4AC62i6G4QWH24BYMyYsRdUXn+on7cFex2DvX4Q/HUM1Prp9bpuf5EGL94ySkpKwmw2a9tms5mEhAS3c6qqqli+fDmZmZk8/fTT3iqKT1mtVqqrnfWud7T/BuEaCor2mCnAxIk5viucEEL0wGuBMHv2bHbt2oXFYsFms5Gfn8/cuXO144qi8OCDD3LDDTfwxBNPeGw9BKJ9+3YDuExR4X6rKCm0Wjs3O3saqamjfFo+IYTojtduGSUmJrJ69WpWrFiB3W5n2bJl5OTksHLlSlatWkVlZSVHjx5FURS2bt0KwKRJkwK6pWAyVXLkyAEAjjd6eoRUJSWyEoARI0Yye/ZVPiydEEL0zGuBALB48WIWL17stu+1114DIDs7m+PHj3vz7X3KZrOxbZtziooyWxLujS/3KSoiIozMnXuNPGYqhPArMlL5Ejl9uoS6OufAukp7+9NVrrfBHFrfwfz51xIbG+vbAgohxHlIIFwi3333rctW1yeLovVWAMaPn8Lo0WN9VzAhhOglCYRLwGq1cu6c8zHSvXWT8DSBXWaUc/GbmTMv83n5hBCiNyQQLoF//esLAI7WjwUMdJ2iogWAmTOvJCqq+2eAhRCiP0kg9NGZMyWcOnUSc3MsjeogD2c4p6iIjBzEpEnZPi+fEEL0lgRCH+3d6xx3cLq5fQ1kz+MOrr02l8jISB+XTgghek8CoQ8URaGiogwA1UNHsoEWUiIriY9PkgFoQgi/J4HQB99+65yzyNl30LUjeWrMMQAWLFjk87IJIcSFkkC4SGazmQMH9mBtNXbTd+BcCe3WW++Rqa2FEAFBAuEiWCwW3n337wAUNI5p29u17+CKK+bKGslCiIAhgXCBbDYbb7/9BuBcJ1nFdfoJZxjosZMSWcnUqTN8X0AhhLhIEggX6OOPN2qvO9ZJdu87mBZzlCFD4nxaLiGE6CsJhAtQUHCUqirnMp9dRyS7T2B3/fWLu15ACCH8mARCL5lMlXz+uXNNaGcYuI5Ibudc/CYzc4JMXieECDgSCL2gKArvv78egO/qJtA1DJytgyxjMQCzZ1/t2wIKIcQlIIHQC++993+Ac41kB6GdjjrDIC38DFEhjYwePUZGJAshApIEwnl88skmamqqKLMlUe+Iadvr3m+QFn6G+HDnbKeXXz6360WEECIASCD04PPPt1BSUoS11UilPbFtr3sYROvrtDAYNy5L+g6EEAFLAqEbJ04cp6DgKOC6PrJ7GOixa+scREQYueyyOT4upRBCXDpeXVM5UJWWFrNt2ydA+xNFejqHAShMi3EGRmRkFNdfv0imqBBCBDRpIXRSWlrMJ598CLg+XtqZwoyY77Wt669fxPDhyR7OE0KIwCGB4MJisbB16xZU1eGy+hl07UQ+q/3M1KmXSxgIIYKCBIKLPXu+prW1qdMMpp0nrTNpncjDhsUzffpM3xdUCCG8QAKhjclUSVHRCaD7TuSkUBMpkZUAxMbGMX/+9YSFhfm4pEII4R3Sqdzm4483Ad13IrevfgaQlTWRKVNmyiOmQoigIoEA7N+/l6amhh7nKGpf/Wz27LlkZ0/FYPDU2SyEEIFrwAfCnj3fsGfPTg9hoGrntM9RdNVVPyI7e7LPyyiEEL4woANh//697Nmzk4N1WeBhoRtwkGUsIiqkkYyMLCZMmNQPpRRCCN8YsIFQUVHOrl072F83HoX2jmHXW0UdYw3i4pKYOXOO3CYSQgS1ARkI5eWnWff+Jxxv9NRn4D7WYPTodGbMmC2jkIUQQW/ABUJRUSHvf/wFxxvHte3pGgbh2IgPt3DDDUtIShoh01kLIQaEARMIVquVw4cP8q89hz2MM4COJTAbmBB9knnzrmP06PQu1xFCiGAV9IFgs9koLS3mnU/3cLY5ETudWwYdTxO1DzzLyBjP2LGZPi+rEEL0J68GwubNm3n55ZdpbW3l3nvv5a677nI7fuzYMZ544gkaGhqYMWMG//Vf/0VISN+LtHXrR5hMlZyqCaGyJZ5WNQSF1LajXVsFoJIWXkZ8uIXc3DxSU0fJCGQhxIDjtakrTCYTa9euZf369WzcuJF33nmHkydPup3z6KOP8pvf/IatW7eiqiobNmzo8/s++p9/ZtMBle1n0yhtTqVZjUAhBGcQeAoDhRkxh5iZOZSf/ewXpKdnSBgIIQYkrwXCzp07mTVrFkOGDMFoNJKbm8uWLVu04+Xl5TQ1NTFlyhQAbr75ZrfjF+OVN9/leOM4mlQjHU8PeQqCjgVuZsR8z3XXLWThwiV9em8hhAh0XrtlVFVVRXx8vLadkJDAoUOHuj0eHx+PyWS6oPeIi4ty2y6raWl71bV/wJ1KrOEcYwad4cc//jEjR468oPftb/Hx0f1dBK8L9joGe/0g+OsYjPXzWiA4HA50uo7fzFVVdds+3/HeqKmx4nB0fOmnDgujvNw9BPQohKAADlSdnriQc2SPaGXChIlkZOQRGRmJ2Vx/gbXrP/Hx0QFV3osR7HUM9vpB8NcxUOun1+u6/CLtymuBkJSUxN69e7Vts9lMQkKC23Gz2axtV1dXux2/GD+951bC3v6A3aWhtKjhDDHUkjG4kkGDopg0aQpJScMZNCiKqKju/0KEEGKg8logzJ49m5deegmLxUJkZCT5+fn87ne/044nJycTHh7Ovn37mD59Ops2bWLu3LkX9B56fdcWxaMP30tNjbXP5fdnnuodbIK9jsFePwj+OgZi/c5XZp2qqt3daO+zzZs388orr2C321m2bBkrV65k5cqVrFq1iuzsbI4fP86aNWuwWq1MnDiRP/7xj/KEjxBC9BOvBoIQQojAIUtoCiGEACQQhBBCtJFAEEIIAUggCCGEaCOBIIQQApBAEEII0UYCQQghBCCBIIQQok3QBMLmzZtZuHAhCxYsYN26df1dnEvmnnvuYdGiRSxZsoQlS5Zw8ODBoKir1WolLy+PsrIywDld+uLFi1mwYAFr167Vzjt27Bg333wzubm5PPHEE7S2tvZXkS9I5/r96le/YsGCBdrnuG3bNqD7evu7v/zlLyxatIhFixbx7LPPAsH3GXqqY7B9jl2oQaCyslKdP3++eu7cObWhoUFdvHixeuLEif4uVp85HA71yiuvVO12u7YvGOp64MABNS8vT504caJ65swZ1WazqVdffbV6+vRp1W63q/fff7+6fft2VVVVddGiRer+/ftVVVXVX/3qV+q6dev6s+i90rl+qqqqeXl5qslkcjuvp3r7s6+//lq9/fbb1ebmZrWlpUVdsWKFunnz5qD6DD3VMT8/P6g+R0+CooVwvsV4AlVxcTEA999/PzfeeCNvvfVWUNR1w4YNPPnkk9rstocOHSItLY3U1FRCQkJYvHgxW7Zs8coiSr7QuX42m42zZ8/y61//msWLF/Piiy/icDi6rbe/i4+P5/HHHycsLIzQ0FDS09MpKSkJqs/QUx3Pnj0bVJ+jJ15dU9lXzrcYT6Cqq6vjiiuu4D//8z+x2+2sWLGCG264IeDr+vTTT7tte/r8TCbTJVlEqT90rl91dTWzZs3iySefJDo6mgceeID33nsPo9Hosd7+bty4cdrrkpISPv30U+6+++6g+gw91XHdunXs3r07aD5HT4KihXApFtvxR1OnTuXZZ58lOjqa2NhYli1bxosvvhh0de3u8wuWzzU1NZW//vWvJCQkEBkZyT333MOXX34Z8PU7ceIE999/P7/85S9JTU0Nys/QtY5jxowJys/RVVAEQufFdjovxhOo9u7dy65du7RtVVVJTk4Ourp29/l5YxGl/lBQUMDWrVu1bVVVCQkJCeh/t/v27eO+++7jP/7jP1i6dGlQfoad6xiMn2NnQREIs2fPZteuXVgsFmw2G/n5+Re82I4/qq+v59lnn6W5uRmr1co//vEP/vznPwddXSdPnsypU6coLS1FURQ++ugj5s6d67aIEnBRiyj5A1VV+cMf/kBtbS12u5133nmH6667rtt6+7uKigoeeughnnvuORYtWgQE32foqY7B9jl6EhR9CImJiaxevZoVK1Zoi/Hk5OT0d7H6bP78+Rw8eJCbbroJh8PB8uXLmT59etDVNTw8nGeeeYaHH36Y5uZmrr76aq6//noAnnvuObdFlFasWNHPpb1wWVlZ/PSnP+XOO++ktbWVBQsWkJeXB9Btvf3Z66+/TnNzM88884y274477giqz7C7OgbT5+iJLJAjhBACCJJbRkIIIfpOAkEIIQQggSCEEKKNBIIQQghAAkEIIUQbCQQxYFxzzTUcPnz4kl/38OHDrFq16rznffDBB8ybN4+f/OQnPZ6XmZmJxWLhgw8+4IEHHuj2vEcffZTCwsILLu/5VFRU8POf/xyHw3HJry38mwSCEH2UnZ3Niy++eN7zNm7cyOrVq3n99df7/J6ffPIJ0dHRZGRk9PlanQ0fPpysrCzWr19/ya8t/FtQDEwTwWHNmjXExcWxevVqwDmqNT8/n7/+9a988cUXvPzyy9jtdiIiInjssceYOnUq1dXV/OY3v6Gmpgaz2UxycjLPP/88cXFxXHPNNeTk5FBQUMAvfvGLXr2Pq84/HxISwiuvvEJLSwsWi4WbbrqJRx55hG+//Zbf/e53fPTRRzz++ONERUVRUFBAZWUlmZmZ/OlPf+KFF17g8OHDlJWVce7cOa6++mqeeuopGhoaMJvNZGVl8fzzzxMeHt6rv6uXXnqJF154AXC2PPLz83E4HJw9e5bExERuu+023nrrLUpKSvjxj3/M/fff3+vzAG699VaWLVvGbbfdRlhYWJ8/WxEg+mHKbSE8Onr0qDpnzhxt/Yfly5erO3bsUE+dOqXm5eWpFotFVVVVLSwsVOfMmaM2NDSob7zxhvrKK6+oqupcP+Lf/u3f1Ndff11VVVWdP3+++pe//EW7/vz589VDhw51+z6duf68w+FQ7777bvXUqVOqqjrXpRg/frxaU1OjfvPNN+qiRYtUVVXVxx57zG0e/Ztuukl97733VFVV1bvvvlv99NNPVVVV1WeeeUbduHGjqqqq2tLSoubl5albtmxRVVVVMzIy1JqaGvX9999Xf/rTn3YpV0FBgTp//nxt+/3331enT5+unj17VlUURV24cKH68MMPq4qiqMeOHVOzs7NVRVF6fV67vLw8ddeuXb3/AEXAkxaC8Bvjx48nJSWF7du3M3r0aKqqqrjyyitZv349VVVV3Hfffdq5Op2O06dPc++997J3717+9re/UVJSwokTJ5g8ebJ23owZM3r9Pp60/7xOp+N//ud/2L59Ox999BFFRUWoqorNZuvyM1dddZX2W3VGRga1tbVdznn00Uf5+uuvee211ygpKaGqqorGxsZe/T0VFxczcuRIt33Z2dkMHz4cgJSUFK688kr0ej2pqak0Nzdr5ezNeYMGDdKOnzp1ilmzZvWqXCLwSSAIv3LXXXfx/vvvM2rUKG677TZtCuUrrriC559/XjuvoqKChIQE/vznP3Po0CFuueUWLr/8clpbW1FdZmMxGo29fh9P2n++sbGRpUuXcu211zJjxgxuueUWPvvsM7f3ahcREaG91ul0Hs/5xS9+gaIo3HDDDcybN4+KigqP53nS/nfiqvNtnZAQz/+1e3seQGhoKAaDoVdlEsFBOpWFX8nNzeXYsWNs3bqVW265BYArrriCr7/+mqKiIgC+/PJLbrzxRpqamvjqq6+49957uemmm4iLi2Pnzp0oinJR79OT0tJSrFYrjzzyCNdccw3ffvstLS0tF/0kzldffcVDDz3EwoULATh48GCvyg0wevRozpw5c1HveyHKysoYM2aM199H+A9pIQi/EhYWRm5uLtXV1cTGxgIwduxYnnrqKX7xi19oc9C//PLLDBo0iIceeohnn32WF154gdDQUKZNm8bp06cv6n16kpmZybx587jhhhsICwsjIyODsWPHUlpaelGdrqtXr+ahhx7CaDQSFRXFzJkze1VucN6GCg8Pp6ioiPT09At+796orq6mpqaGadOmeeX6wj/JbKfCrzQ2NnL33Xfzm9/8RluHN5Dfx1s2b97Mvn37+O1vf+uV67/00kvExsZy1113eeX6wj/JLSPhN/71r38xb948rrrqKq9+Sfvqfbxp8eLF1NbWUlBQcMmvXVFRwZEjR7jjjjsu+bWFf5MWghBCCEBaCEIIIdpIIAghhAAkEIQQQrSRQBBCCAFIIAghhGgjgSCEEAKA/x8nfVDbNQXaqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _ in range(50):\n",
    "    # Generate bootstrap samle: bs_sample\n",
    "    bs_sample = np.random.choice(rainfall, size=len(rainfall))\n",
    "    \n",
    "    # Compute and plot ECDF from bootstrap sample\n",
    "    x, y = ecdf(bs_sample)\n",
    "    _ = plt.plot(x, y, marker='.', linestyle='none', color='gray', alpha=0.1)\n",
    "    \n",
    "# Compute and plot ECDF from original data\n",
    "x, y = ecdf(rainfall)\n",
    "_ = plt.plot(x, y, marker='.')\n",
    "\n",
    "# Make margins and label axes\n",
    "plt.margins(0.02)\n",
    "_ = plt.xlabel('yearly rainfall (mm)')\n",
    "_ = plt.ylabel('ECDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrap confidence intervals\n",
    "- Confidence interval of a statistics\n",
    "    - If we repeated measurements over and over again, ```p%``` of the observed values would lie within the ```p%``` confidence interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating many bootstrap replicates\n",
    "The function ```bootstrap_replicate_1d()``` from the video is available in your namespace. Now you'll write another function, ```draw_bs_reps(data, func, size=1)```, which generates many bootstrap replicates from the data set. This function will come in handy for you again and again as you compute confidence intervals and later when you do hypothesis tests.\n",
    "\n",
    "For your reference, the ```bootstrap_replicate_1d()``` function is provided below:\n",
    "```python\n",
    "def bootstrap_replicate_1d(data, func):\n",
    "    \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n",
    "    bs_sample = np.random.choice(data, len(data))\n",
    "    return func(bs_sample)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_replicate_1d(data, func):\n",
    "    \"\"\"Generate bootstrap replicate of 1D data.\"\"\"\n",
    "    bs_sample = np.random.choice(data, len(data))\n",
    "    return func(bs_sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_bs_reps(data, func, size=1):\n",
    "    \"\"\"Draw bootstrap replicates.\"\"\"\n",
    "    \n",
    "    # Initialize array of replicates: bs_replicates\n",
    "    bs_replicates = np.empty(size)\n",
    "    \n",
    "    # Generate replicates\n",
    "    for i in range(size):\n",
    "        bs_replicates[i] = bootstrap_replicate_1d(data, func)\n",
    "        \n",
    "    return bs_replicates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap replicates of the mean and the SEM\n",
    "In this exercise, you will compute a bootstrap estimate of the probability density function of the mean annual rainfall at the Sheffield Weather Station. Remember, we are estimating the mean annual rainfall we would get if the Sheffield Weather Station could repeat all of the measurements from 1883 to 2015 over and over again. This is a probabilistic estimate of the mean. You will plot the PDF as a histogram, and you will see that it is Normal.\n",
    "\n",
    "In fact, it can be shown theoretically that under not-too-restrictive conditions, the value of the mean will always be Normally distributed. (This does not hold in general, just for the mean and a few other statistics.) The standard deviation of this distribution, called the standard error of the mean, or SEM, is given by the standard deviation of the data divided by the square root of the number of data points. I.e., for a data set, ```sem = np.std(data) / np.sqrt(len(data))```. Using hacker statistics, you get this same result without the need to derive it, but you will verify this result from your bootstrap replicates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9488593574676786\n",
      "0.956034708235844\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAeI0lEQVR4nO3deXgU9eHH8U8SQiAcIpjDCmp9aLhPtSAIiAQQJAWEKnJFLSAVpEbLUZsiBASkQhRtaQSrLRIJKkeDj+EQ8RFIAZESReRQDEFCsiQIBHLufn9/uNkfkVyQTHaTvF/Pw/NkdiYzH7K7+WRmdr7jZYwxAgDUet7uDgAA8AwUAgBAEoUAAHCiEAAAkigEAIAThQAAkEQhAACc6rg7QEWcO3dJDkfFL6No1qyhMjKyKiGRNchXMeSrGPJVjCfl8/b20o03NihxfrUuBIfDVEohFK7Lk5GvYshXMeSrGE/PV4hDRgAASRQCAMCJQgAASKIQAABOFAIAQBKFAABwohAAAJKq+XUIgKdq1Li+6vld/fbKyS3QxQvZbkgElI1CACxQz6+Owp7beNXj8UuG6qIb8gDlwSEjAIAkCgEA4EQhAAAkUQgAACcKAQAgiUIAADhRCAAASRQCAMCJQgAASKIQAABOFAIAQBKFAABwohAAAJIY7RQoU0lDWUsMZ42axdJCiI+P1/Lly1VQUKDw8HCNGTOm2OV27NihqKgobd++3co4wHUpaShrieGsUbNYVghpaWmKjo7WunXrVLduXY0aNUrdunVTy5Ytiyx39uxZvfTSS1bFACyVl29XQEAjd8cAKoVlhbB79251795dTZo0kSQNHDhQCQkJmjp1apHlIiMjNXXqVC1ZssSqKIBl6vr6lHgjnGvBHdbgCSwrhPT0dAUEBLimAwMDlZSUVGSZf//732rbtq06dep0Xdto1qxhhTJeydP/yiNfxXhSvuKylHaHtXoekN2Tfn7FIV/lsKwQHA6HvLy8XNPGmCLTR48e1ZYtW/T222/rzJkz17WNjIwsORymwlkDAhrJZvPcI8Hkq5iK5qvsN/PPs5S1fnf/bGv682s1T8rn7e1V6h/Sln3sNDg4WDabzTVts9kUGBjomk5ISJDNZtOIESM0adIkpaena/To0VbFAQCUwbJC6NGjhxITE5WZmans7Gxt2bJFvXv3ds2fNm2aNm/erI0bN+qNN95QYGCgYmNjrYoDACiDZYeMgoKCFBERofHjxys/P18jR45Ux44dNXHiRE2bNk0dOnSwatPAdSntegOgNrD01R8WFqawsLAij61YseKq5Zo3b841CHC70k7sArUBQ1cAACQxdAVQpbiQDZ6MQgCqUGVdyAZYgUNGAABJFAIAwIlCAABIohAAAE6cVAY8WEmfSmIUVFiBQgA8WGmfSvKM4dJQk3DICAAgiUIAADhRCAAASRQCAMCJQgAASKIQAABOFAIAQBKFAABwohAAAJIoBACAE4UAAJBEIQAAnCgEAIAkCgEA4MTw16h1GjWur3p+vPSBn+NdgVqnnl+dEu8xANRmHDICAEiiEAAAThQCAEAS5xCAaikv366AgEZXPZ6TW6CLF7LdkAg1AYUAVEN1fX1KPDF+0Q15UDNwyAgAIIlCAAA4ccgINdaVF6AVd7wdQFEUAmosLkADrg2HjAAAkigEAIAThQAAkEQhAACcLC2E+Ph4DR48WAMGDNDq1auvmr9161aFhYXpwQcf1KxZs5SXl2dlHABAKSwrhLS0NEVHRys2NlYbNmxQXFycjh8/7pp/+fJlRUVF6a233tKHH36o3NxcrV+/3qo4AIAyWFYIu3fvVvfu3dWkSRP5+/tr4MCBSkhIcM339/fX9u3bddNNNyk7O1sZGRlq3LixVXEAAGWw7DqE9PR0BQQEuKYDAwOVlJRUZBlfX199+umnmjFjhgIDA3Xvvfde0zaaNWtYKVklz79wiXwoLyueC09/fslXOSwrBIfDIS8vL9e0MabIdKE+ffpoz549Wrp0qebMmaMlS5aUexsZGVlyOEyFswYENJLN5rlDgpHv+lSXN2Flq+znwlOf30LkKz9vb69S/5C27JBRcHCwbDaba9pmsykwMNA1/eOPP2rnzp2u6bCwMB05csSqOACAMlhWCD169FBiYqIyMzOVnZ2tLVu2qHfv3q75xhhNnz5dp0+fliQlJCSoa9euVsUBAJTBskNGQUFBioiI0Pjx45Wfn6+RI0eqY8eOmjhxoqZNm6YOHTpo3rx5evLJJ+Xl5aWWLVtq7ty5VsUBAJTB0sHtwsLCFBYWVuSxFStWuL4ODQ1VaGiolREAAOXElcoAAEkMfw3UKNxrGRVBIQA1CPdaRkVQCKj2rrwzGoDrx7sI1R53RgMqByeVAQCSKAQAgBOFAACQRCEAAJwoBACAJAoBAOBEIQAAJFEIAAAnCgEAIIlCAAA4UQgAAEmMZQTUCgyLjfIotRAOHTqkdu3aVVUWABZhWGyUR6mHjCIjI11f//3vf7c8DADAfUotBGOM6+utW7daHgYA4D6lFoKXl5fr6yvLAQBQ85T7U0ZXlgMAoOYp9aTymTNnNH/+/Ku+LnTlOQbAatwqE7BWqe+uMWPGFPs14A7cKhOwVqmFMHXq1KrKAQBwszL3v5OSkvSvf/1LR48eVb169RQSEqLw8HCFhIRURT4AQBUp9aRyYmKinnrqKYWEhOiPf/yjpk6dqptvvllPPPGE9u7dW1UZAQBVoNQ9hJiYGK1cuVKtW7d2PdanTx/17t1bS5Ys0a9//WvLAwIAqkapewgZGRlFyqBQx44ddfEiF7wDQE1SaiH4+PiUOI8L1QCgZin1kBEXowE1G6Og4kqlFsJ3332nsLCwYuelpKRYEghA1WEUVFyp1EJYsWKFjh49qoCAAOXm5io4OLiqcgEAqliphZCSkqJly5bptttu08mTJ7VkyRLde++9VZUNAFCFSi2EVatWKT4+XkFBQTpw4ICio6MpBACoococ7TQoKEiS1KVLF507d87yQAAA9yj3/RCk0j+GCgCo3sp9PwSJj6ECQE1W6jmEI0eOqGvXrq7pnJwcde3aVcYYeXl56Ysvvih15fHx8Vq+fLkKCgoUHh5+1RDa27Zt02uvvSZjjJo3b66FCxfqhhtuqMB/BwBwvUothIrcRzktLU3R0dFat26d6tatq1GjRqlbt25q2bKlJCkrK0tz5szRBx98oKCgIL366qt67bXXuOkOALhJqYVwyy23XPeKd+/ere7du6tJkyaSpIEDByohIcF1j4X8/Hy98MILrpPWrVq1Unx8/HVvDwBQMZbdjzA9PV0BAQGu6cDAQCUlJbmmb7zxRvXv31/ST4ei3njjDY0bN+6attGsWcPKCSsVe/m+JyEfqtqVz6mnP7/kqxyWFYLD4ShyErrwvMPPXbx4UVOmTFHr1q01fPjwa9pGRkaWHI6KD7IXENBINpvnXqhPvv/fDqpO4XPK669iPCmft7dXqX9IX9OnjK5FcHCwbDaba9pmsykwMLDIMunp6Ro9erRatWqlF1980aooAIBysKwQevToocTERGVmZio7O1tbtmxR7969XfPtdrsmT56sQYMG6c9//jMfaQUAN7PskFFQUJAiIiI0fvx45efna+TIkerYsaMmTpyoadOm6cyZM/r6669lt9u1efNmSVL79u3ZU4AaNa6ven6WvTQBlMDSd11YWNhVw2evWLFCktShQwd98803Vm4e1VQ9vzolDskMwDqWHTICAFQvFAIAQBKFAABwohAAAJIoBACAE4UAAJBEIQAAnCgEAIAkCgEA4EQhAAAkWTx0BYDqKS/fXuz9EHJyC3TxQra7YsFiFAKAq9T19SlxPCnPGNkfVuCQEQBAEoUAAHCiEAAAkigEAIAThQAAkEQhAACc+Ngp3IZ7JwOehXcj3IZ7JwOehUNGAABJFAIAwIlCAABIohAAAE4UAgBAEoUAAHCiEAAAkigEAIAThQAAkEQhAACcGLoClmPMIqB64F0KyzFmEVA9UAioNOwJ1Hx5+XYFBDS66vGc3AJdvJDthkSoTLx7UWnYE6j56vr6lPgcX3RDHlQuTioDACRRCAAAJwoBACCJQgAAOFlaCPHx8Ro8eLAGDBig1atXl7jcjBkztG7dOiujAADKYFkhpKWlKTo6WrGxsdqwYYPi4uJ0/Pjxq5aZPHmyNm/ebFUMAEA5WVYIu3fvVvfu3dWkSRP5+/tr4MCBSkhIKLJMfHy8+vXrp0GDBlkVAwBQTpZdh5Cenq6AgADXdGBgoJKSkoosM2HCBEnS/v37r2sbzZo1vP6AP1PcxTaehHzwdO58DXj668/T8xWyrBAcDoe8vLxc08aYItOVISMjSw6HqfB6AgIayWbz3Mtqqku+6vKihzXc9RqtLu8PT+Dt7VXqH9KWHTIKDg6WzWZzTdtsNgUGBlq1OQBABVlWCD169FBiYqIyMzOVnZ2tLVu2qHfv3lZtDgBQQZYVQlBQkCIiIjR+/HgNGzZMQ4YMUceOHTVx4kR9+eWXVm0WAHCdLB3cLiwsTGFhYUUeW7FixVXLLVq0yMoYACxW0iioEiOhVieMdgqgwkoaBVViJNTqhEIAYCnuoVB9UAgALMU9FKoPBrcDAEhiDwHXobhbZXJRGlD9UQi4ZtwqE6iZOGQEAJBEIQAAnCgEAIAkCgEA4EQhAAAkUQgAACcKAQAgiUIAADhxYRqKVdzVyABqNt7xKFZJVyNLXJEM1FQUAgC3YFhsz0MhAHALhsX2PJxUBgBIohAAAE4UAgBAEoUAAHDipDIAj8Knj9yHQqjluAANnoZPH7kPvwlqOW6HCaAQhVBLsCcAoCz8hqgl2BMAUBY+ZQQAkMQeAoBqoqRPH+Xl292QpmaiEABUC6V9+giVg0NGAABJFAIAwIlDRgBqpJI+as0VzyWjEABUayWdbJbEFc/XiEIAUK1xsrnyUAg1DFckA7he/Oaopkr7xc9fS0DJGE21ZBSCh2vUuL4kFfsC5hc/cO1KOsT0waIhtb4oLC2E+Ph4LV++XAUFBQoPD9eYMWOKzD98+LD+/Oc/69KlS7rrrrs0d+5c1alDR12JMYiAqsGw2xZeh5CWlqbo6GjFxsZqw4YNiouL0/Hjx4ssM336dM2ePVubN2+WMUZr1661Ko7Ha9S4vgICGl31D4B7FR5i+vm/wr33msSyP8d3796t7t27q0mTJpKkgQMHKiEhQVOnTpUk/fDDD8rJyVHnzp0lSQ899JCWLVum0aNHl3sb3t5elZa3MtclSQ0b1pNfMcf4c/Ps8qvrU+z3/G7+lqseezNygAJvLP6F567H3bltT3vcEzN52uOemOlaHq/r61Pie/NSMb83invvF/5xV9L7Pze3QFlZOcVmqkxl/Z7zMsYYKzYcExOjy5cvKyIiQpL03nvvKSkpSfPmzZMkHThwQIsXL9a7774rSUpOTtakSZO0efNmK+IAAMpg2SEjh8MhL6//byNjTJHpsuYDAKqWZYUQHBwsm83mmrbZbAoMDCxx/tmzZ4vMBwBULcsKoUePHkpMTFRmZqays7O1ZcsW9e7d2zX/lltukZ+fn/bv3y9J2rhxY5H5AICqZdk5BOmnj53GxMQoPz9fI0eO1MSJEzVx4kRNmzZNHTp00DfffKPIyEhlZWWpXbt2WrhwoerWrWtVHABAKSwtBABA9cH9EAAAkigEAIAThQAAkEQhAACcat1Ictu3b9frr7+u7Oxs9ezZU5GRkYqNjdXq1atljFGfPn00Y8YMt10kV1y+Qu+88442b96sVatWuSVbSfn+9Kc/af/+/apf/6fL/qdOnar+/ft7TL4DBw5o4cKFunTpklq1aqVFixa57dNsP8/Xq1cvLV261DU/LS1NnTp1UkxMjEfki4yM1M6dO7V48WI5HA61bdtW8+fP95ifX2RkpNatW6eVK1fKx8dH3bp106xZs9wySOZ7772nd955xzV96tQpDR06VKGhoVq4cKFyc3M1aNAg1+gNHsnUIidPnjT33nuvSU1NNXl5eebRRx81q1evNv379zeXLl0yBQUF5pFHHjGfffaZx+TbsWOHMcaYY8eOmV69epmxY8e6JVtp+YYMGWLS0tLclqu0fFu3bjU9e/Y0hw8fNsYYExERYVavXu0x+QqfX2OMSU9PN/369TMnTpzwqHy9e/c2x48fN8YY8/TTT5u1a9d6TL63337b9OrVy/X6e+GFF8w///lPt+S70tGjR03//v3N6dOnTZ8+fczJkydNfn6+eeKJJ4o8556mVh0y2rp1qwYPHqzg4GD5+voqOjpagwcP1ocffih/f39duHBBWVlZaty4scfk69Spk/Ly8jR79mxNmzbNLblKy9e6dWudPn1azz//vMLCwrRs2TI5HA6PyWe329W5c2e1bt1akhQZGem2vZeSnt9Cixcv1qhRo3T77bd7VD673a6srCzZ7Xbl5ubKz8/PY/IFBASoc+fOrlEO+vbtq23btrkl35XmzJmjiIgIpaSk6LbbblOLFi1Up04dhYWFKSEhwd3xSlSrCiE5OVl2u12TJ0/W0KFDFRsbqxtuuEG+vr5au3atQkNDFRAQ4Prl4Sn5lixZohEjRqhFixZuyVVavtzcXHXv3l0LFizQ2rVr9fnnn+v999/3mHzJycny9/dXRESEhg4dqtdee81thV/S8ytJ33//vfbu3avx48e7JVtp+ebMmaNx48apV69eOnfunB544AGPyde6dWsdPHhQqampstvtSkhI0NmzZ92Sr9Du3buVk5OjQYMGKT09XQEBAa55gYGBSktLc2O60tWqQrDb7UpMTNSCBQsUFxenpKQkrV+/XpL08MMPa8+ePbrpppv0+uuve0y+9957T6mpqRoxYoRbMpWV7/PPP9ff/vY3BQYGqn79+ho3bpw+/fRTj8lnt9u1c+dOPfvss1q3bp2ys7P1xhtveEy+wtdfXFycRo8e7dYr9YvLt2LFCr388svatGmTdu7cqU6dOmnhwoUek+9///ufnnvuOf3+97/XmDFj1KpVK/n6+rolX6E1a9bo8ccfl1T9BvGsVYVw00036Z577lHTpk1Vr149hYaGat++fa7xlOrUqaMHH3xQR44c8Zh8Bw4c0LFjxzR06FBFRkbqq6++0jPPPOMx+davX19kyHJjjNvueldcvuXLl6tTp05q0aKFfHx8NGjQICUlJXlMvsIsH3/8sQYPHuyWXKXlW79+vUJCQnTrrbfK29tbDz/8sPbu3esx+fbt26eOHTtqw4YNWrNmjYKCgty6J52Xl6d9+/bp/vvvl1T2IJ+eplYVQt++fbVz505duHBBdrtdn332mVq3bq3p06frwoULMsZo8+bNuvPOOz0mX9euXfXRRx9p48aNmj9/vtq3b69XXnnFY/KFhoZqwYIFOn/+vPLz8xUXF+e2Y/TF5Zs0aZIOHTqk1NRUSdInn3yidu3aeUy+du3aKTMzUzk5OW4/JFhcvrFjxyopKcl1GObjjz9Whw4dPCbfr371Kz322GPKyspSXl6e3nnnHbcW65EjR3T77bfL399fktSpUyedOHHCdbhr06ZNHj2IZ6362GmnTp00YcIEjR49Wvn5+erZs6fGjRsnPz8/jRo1Sj4+Prrrrrtcu3uekM8TDhUVKunnV6dOHT366KMqKCjQgAEDNGTIEI/J99RTT6l9+/aaPHmycnNz1aZNG82cOdNj8o0YMUJfffWVgoOD3ZKprHyPPvqo/P39NX78ePn4+Oi2225TVFSUx+R77LHH1KhRIz3yyCMqKCjQkCFDFBYW5pZ8kpSSklLkufTz89OiRYv09NNPKzc3V3369HHbOZjyYHA7AICkWnbICABQMgoBACCJQgAAOFEIAABJFAIAwIlCACrJrFmz9Oabb1ZoHUOHDtWFCxdKXebw4cMKDQ3VQw89pFOnTpW43Lhx45SQkKBTp06pS5cuJS4XGxuruLi4685cErvdrieffFIZGRmVvm5Yg0IAPMjGjRvLHGvp448/Vrdu3bRu3To1b968Qtv74YcftH79ej388MMVWk9xfHx8NGHCBM2dO7fS1w1r1KoL01D59uzZo6VLl+rmm2/WiRMnVL9+fU2aNEmrVq3SiRMnNGDAAD3//POSfhrLfvny5crPz1e9evU0c+ZMdenSRWfPntXs2bOVkZEhm82mW265Ra+88oqaNWum+++/X8OHD1diYqJSU1M1dOjQYofu+OSTTxQTE6O8vDxlZmZq2LBheuaZZ7Rnzx5FR0erRYsWOnbsmAoKCjR37lzdeeedmjVrlho2bKgjR47ozJkzatWqlV566SU1aNBArVq1UmJiopo2bSpJrukmTZpowYIFOnjwoC5duiRjjObPn1/q1e2zZs3Sjz/+qJSUFN13330aOXKkoqKidOnSJdlsNrVu3VqvvPKK/Pz8XNvZsWOHtm7dKm9vbyUnJ6tevXp66aWXdOjQIb377ruy2+3KycnRvHnzNGfOHCUnJ+vHH39UgwYN9PLLL+uOO+4o1/MXExOjoUOHysvLS6dOnVJ4eLh69uypr776Sna7XdOmTVNcXJy+++47tW/fXkuXLtXp06fLtZy3t7fuvvtuvfDCCzp8+LDatGlzHa8wVCl3jr2N6u+///2vadOmjTl06JAxxpjf/e535pFHHjG5ubkmIyPDtGvXzpw5c8acOHHCDBkyxGRmZhpjfhovvmfPnubSpUvm7bffNjExMcYYYxwOh5kwYYJ58803jTHG9O3b1yxatMgYY8yZM2dMhw4dzMmTJ4tkcDgcZuzYsa77CJw5c8a0adPGZGRkuPJ9/fXXxhhj3nzzTTNmzBhjjDEzZ850Zc3LyzPDhg0z77//vjHGmJCQEJORkeHaRuH0F198YZ5++mljt9uNMcbExMSYJ5980rW+lStXXvUzmjlzpgkPD3dNL1q0yGzYsMEYY0xeXp4ZMmSISUhIKLKdDz74wNx5550mNTXVGGNMVFSUmTFjhjHGmGXLlpm5c+caY4z56KOPzLx581zr/stf/mKioqKMMcaMHTvWfPTRRyYlJcV07tz5qlwOh8N069bNpKSkGGOMSUlJMSEhIWbbtm3GGGNmz55t+vbtay5evGhycnJMz549zf79+8u9XKF58+aZV1999artw/Owh4AKa968udq2bStJuvXWW9WoUSPVrVtXTZs2VYMGDXT+/Hnt27dP6enpeuyxx1zf5+XlpZMnTyo8PFyff/653nrrLX3//fc6duxYkfsE9OvXT5IUFBSkZs2a6fz580XG/fHy8tI//vEP7dixQ5s2bdK3334rY4yys7MlSb/4xS9cf522bdvWNcKoJPXq1cs1wmhISIjOnz9f6v+1S5cuuuGGG7RmzRqlpKRoz549atCgQZk/oyv3IKZPn65du3ZpxYoV+v7775Wenq7Lly9f9T3t2rVzDYPQtm1bbd269aplHnjgAbVo0UKrVq1ScnKy9u7dW+r5giudO3dOFy9eLHLYydfX1zUw26233qouXbqoYcOGkn4auvn8+fMKDAws13KFmjdvroMHD5YrE9yLQkCF/XzI5uJGO3U4HLrnnnuKDMyXmpqqwMBA/fWvf1VSUpJGjBihbt26qaCgQOaKEVWuvCGLl5dXkXmSdPnyZQ0fPlyhoaG66667NGLECG3bts21XL169Ur8/tLmFcrLy3N9vWPHDr344ot6/PHH1a9fP91xxx36z3/+U/IPx6lwsDNJevbZZ2W32zVo0CDdd999Sk1NLXa75ckWGxurtWvXasyYMQoLC1OTJk1KPdF8pcJ1OhwOeXv/dDrR19e3yPDMJQ0lXd7lpJ9eD4Xrh2fjWUKVuOeee7Rr1y59++23kqRPP/1Uv/nNb5STk6OdO3cqPDxcw4YNU7NmzbR7927Z7fZyrzs5OVlZWVl65plndP/992vPnj3Ky8ur0J3bmjZtqi+//FKStGnTJtfju3btUt++fTV69Gi1b99e27Ztu6askrRz505NmTLFNSrnwYMHr3kdV65r+PDh+u1vf6tf/vKX2r59e7nXdeONN6px48b64Ycfrmvb5XXq1Klyn9OAe7GHgCrRsmVLRUVF6dlnn3XdM2H58uVq0KCBpkyZosWLF+vVV1+Vr6+vunbtqpMnT5Z73a1atdJ9992nQYMGqW7dugoJCVHLli2VnJx83TeciYyMVFRUlBo3bqwePXq47no1atQoPffccwoLC1NBQYF69uypLVu2XFP5REREaMqUKfL391fDhg119913X9P/90pPPPGEZs+e7bpLXefOnXX06NFyf/+AAQP02WefafTo0de1/fLYtWuX24Zsx7VhtFOgFktJSdEf/vAHffDBB5bcyWvPnj1avXq1li1bVunrRuXjkBFQi7Vo0ULDhg3TmjVrKn3ddrtdK1euVGRkZKWvG9ZgDwEAIIk9BACAE4UAAJBEIQAAnCgEAIAkCgEA4EQhAAAkSf8HS2m8y5IvbacAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Take 10,000 bootstrap replicates of the mean: bs_replicates\n",
    "bs_replicates = draw_bs_reps(rainfall, np.mean, size=10000)\n",
    "\n",
    "# Compute and print SEM\n",
    "sem = np.std(rainfall) / np.sqrt(len(rainfall))\n",
    "print(sem)\n",
    "\n",
    "# Compute and print standard deviation of bootstrap replicates\n",
    "bs_std = np.std(bs_replicates)\n",
    "print(bs_std)\n",
    "\n",
    "# Make a histogram of the results\n",
    "_ = plt.hist(bs_replicates, bins=50, density=True)\n",
    "_ = plt.xlabel('mean annual rainfall (mm)')\n",
    "_ = plt.ylabel('PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence intervals of rainfall data\n",
    "A confidence interval gives upper and lower bounds on the range of parameter values you might expect to get if we repeat our measurements. For named distributions, you can compute them analytically or look them up, but one of the many beautiful properties of the bootstrap method is that you can take percentiles of your bootstrap replicates to get your confidence interval. Conveniently, you can use the ```np.percentile()``` function.\n",
    "\n",
    "Use the bootstrap replicates you just generated to compute the 95% confidence interval. That is, give the 2.5th and 97.5th percentile of your bootstrap replicates stored as ```bs_replicates```. What is the 95% confidence interval?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([64.86782928, 68.62363296])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.percentile(bs_replicates, [2.5, 97.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap replicates of other statistics\n",
    "We saw in a previous exercise that the mean is Normally distributed. This does not necessarily hold for other statistics, but no worry: as hackers, we can always take bootstrap replicates! In this exercise, you'll generate bootstrap replicates for the variance of the annual rainfall at the Sheffield Weather Station and plot the histogram of the replicates.\n",
    "\n",
    "Here, you will make use of the draw_bs_reps() function you defined a few exercises ago."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEJCAYAAACUk1DVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfVBU9f4H8PeCICAQKbvQLWtul9TkwacMRCVNhVQ2tBgjH6BU0NKo1VATRhLUtFRMK0OaOz1oXQ1RL3UDLMsmIBvTCfCGVlwVS2CFElAelt3v7w+X/bkCy4McdhferxlnOHvOnvPZr7v73vM953yPTAghQEREfZ6NuQsgIiLLwEAgIiIADAQiItJjIBAREQAGAhER6TEQiIgIAAOBiIj0+pm7gNvx55/XoNP13csoBg1yRmVlrbnLsAhsC2NsD2NsjxtsbGS4884Bbc636kDQ6USfDgQAff7134xtYYztYYzt0T52GREREQAGAhER6TEQiIgIAAOBiIj0GAhERASAgUBERHoMBCIiAiDxdQiZmZnYvXs3mpqaEBUVhXnz5hnNLykpQWJiIq5evQq5XI7t27fjjjvukLIkIrNycXWEQ/+WH7v6hibUVNeZoSKi/ydZIJSXlyMlJQUZGRmwt7dHREQE/P394eXlBQAQQuC5555DfHw8goKCsHXrVuzZswdxcXFSlURkdg79+0G58kiLxzO3haHGDPUQ3UyyLqO8vDwEBATAzc0NTk5OCAkJQVZWlmH+mTNn4OTkhKCgIADA0qVLW+xBEBFRz5EsECoqKiCXyw3TCoUC5eXlhumLFy/C3d0da9euxezZs5GYmAgnJyepyiEionZI1mWk0+kgk8kM00IIo+mmpib88MMP2Lt3L3x9fbFjxw5s3rwZmzdv7vA2Bg1y7taarZFc7mLuEizG7bZFo0YLezvbDj/e3br7/5LvDWNsj/ZJFgienp44efKkYVqtVkOhUBim5XI57rvvPvj6+gIAQkNDERsb26ltVFbW9ukBq+RyF6jV7HkGuqct5HKXNvv321p3Zw8Sm/pS6s7/S743jLE9brCxkZn8IS1ZIAQGBmLXrl2oqqqCo6MjcnJykJycbJg/atQoVFVVobi4GMOGDcOxY8fg7e0tVTlEkuBBYupNJAsEDw8PqFQqREZGQqPRIDw8HH5+foiOjkZsbCx8fX3x9ttvIyEhAXV1dfD09MTrr78uVTlEPapRo2UXBVkdSa9DUCqVUCqVRo+lpaUZ/h4xYgTS09OlLIHILOztbNvccyCyVFZ9gxwiU7rrIjD+2qe+goFAvVZ39e+39Wu/eV1EvQXHMiIiIgDcQ6A+iF1ARK1jIFCfwwO+RK1jlxEREQFgIBARkR4DgYiIADAQiIhIj4FAREQAGAhERKTHQCAiIgAMBCIi0mMgEBERAAYCERHpMRCIiAgAA4GIiPQYCEREBICBQEREegwEIiICwEAgIiI9BgIREQFgIBARkR4DgYiIADAQiIhIT9JAyMzMxIwZMxAcHIx9+/a1mP/WW29h8uTJCAsLQ1hYWKvLEBFRz+gn1YrLy8uRkpKCjIwM2NvbIyIiAv7+/vDy8jIsU1RUhO3bt2PUqFFSlUFERB0k2R5CXl4eAgIC4ObmBicnJ4SEhCArK8tomaKiIqSmpkKpVCIpKQkNDQ1SlUNERO2QLBAqKiogl8sN0wqFAuXl5Ybpa9eu4cEHH0RcXBwOHTqE6upqvPPOO1KVQ0RE7ZCsy0in00EmkxmmhRBG0wMGDEBaWppheuHChVi7di1UKlWHtzFokHP3FGvF5HIXc5dgMay9Lbq7fmtvj+7G9mifZIHg6emJkydPGqbVajUUCoVh+o8//kBeXh7Cw8MB3AiMfv06V05lZS10OtE9BVshudwFanWNucuwCK21hbV9AXTn/yXfG8bYHjfY2MhM/pCWrMsoMDAQ+fn5qKqqQl1dHXJychAUFGSY7+DggDfeeAOlpaUQQmDfvn2YNm2aVOUQEVE7JAsEDw8PqFQqREZGYtasWQgNDYWfnx+io6NRWFiIgQMHIikpCc899xwee+wxCCHw7LPPSlUOERG1Q7IuIwBQKpVQKpVGj9183CAkJAQhISFSlkC9iIurIxz6t3zLNjRqAVhfF9HNGjXaVuuvb2hCTXWdGSqivkjSQCDqTg79+0G58kiLxzO3hbX5uLWwt7Nt8zWw55t6CoeuICIiAAwEIiLSY5cRkQXjsQXqSQwEIgvGYwvUk9hlREREABgIRESkx0AgIiIADAQiItJjIBAREQAGAhER6TEQiIgIAAOBiIj0GAhERASAgUBERHoMBCIiAsBAICIiPQYCEREB4GinZIHaulUmEUmLnzqyOKZulUlE0mGXERERAWAgEBGRHgOBiIgAMBCIiEiPgUBERAAkDoTMzEzMmDEDwcHB2LdvX5vLffPNN3j00UelLIWIiNoh2Wmn5eXlSElJQUZGBuzt7REREQF/f394eXkZLXflyhVs2bJFqjKIiKiDJNtDyMvLQ0BAANzc3ODk5ISQkBBkZWW1WC4hIQHLly+XqgyiXqlRo4Vc7tLin4uro7lLIysm2R5CRUUF5HK5YVqhUKCgoMBomQ8//BDDhw/HiBEjpCqDqFeyt7Nt8+K9GjPUQ72DZIGg0+kgk8kM00IIo+lz584hJycH77//PsrKyrq0jUGDnG+7Tmsnl7uYuwSyMM3vCb43jLE92idZIHh6euLkyZOGabVaDYVCYZjOysqCWq3Gk08+CY1Gg4qKCsydOxcff/xxh7dRWVkLnU50a93WRC53gVptvb8HOWaRNNTqGqt/b3Q3tscNNjYykz+kJfs0BgYGYteuXaiqqoKjoyNycnKQnJxsmB8bG4vY2FgAwKVLlxAZGdmpMCDrxzGLiCyLZAeVPTw8oFKpEBkZiVmzZiE0NBR+fn6Ijo5GYWGhVJslIqIuknR/XalUQqlUGj2WlpbWYrl77rkHx44dk7IUIiJqB69UJiIiAAwEIiLS4ykeRL1I8wVrgPFplvUNTaiprjNXWWQlGAhEvQgvWKPbwUAgyfF6AyLrwE8pSY7XGxBZBwYCUR9w87GFm/HYAt2MgUDUB/DYAnUETzslIiIADAQiItJjIBAREQAGAhER6ZkMhDNnzvRUHUREZGYmAyEhIcHw9zvvvCN5MUREZD4mA0GI/78b2dGjRyUvhoiIzMdkINx6T2QiIuq9OnxQ+eZwICKi3sfklcplZWXYsGFDi7+b3XyMgYiIrJvJQJg3b16rfxMRUe9jMhCWL1/eU3UQEZGZtTu4XUFBAT744AOcO3cODg4OGDJkCKKiojBkyJCeqI+IJMRRUOlmJgMhPz8fcXFxWLBgAR5//HEAQGFhIRYuXIjt27fj4Ycf7pEiiUgaHAWVbmYyEFJTU/Hee+9h2LBhhsceeeQRBAUFYdu2bQwEIqJexORpp5WVlUZh0MzPzw81Nfz9QETUm5gMBFtb2zbn8UI1IqLepcNXKndFZmYmZsyYgeDgYOzbt6/F/KNHj0KpVGLmzJlYs2YNGhsbb2t7RETUdSaPIZSUlECpVLY6r7S01OSKy8vLkZKSgoyMDNjb2yMiIgL+/v7w8vICAFy/fh1JSUk4dOgQ3N3doVKpcOjQITz11FNdfClERHQ7TAZCWloazp07B7lcjoaGBnh6enZ4xXl5eQgICICbmxsAICQkBFlZWYZrG5ycnHDs2DHY2dmhrq4OlZWVcHV1vY2XQkREt8NkIJSWlmLnzp247777cPHiRWzbtg0TJkzo0IorKiogl8sN0wqFAgUFBUbL2NnZ4fjx41i1ahUUCkWH101ERN3PZCB89NFHyMzMhIeHB06fPo2UlJQOf2nrdLoWo6W2dkzikUcewYkTJ7B9+3a8+uqr2LZtW4eLHzTIucPL9latXVRE1B1623urt70eKbR7pbKHhwcAYNSoUfjzzz87vGJPT0+cPHnSMK1Wq6FQKAzTf/31F4qKigwBo1QqoVKpOrx+AKisrIVO13fPdpLLXaBWW/7pv/wgWidreG91lLV8VqRmYyMz+UO6U2cZmToN9VaBgYHIz89HVVUV6urqkJOTg6CgIMN8IQTi4uLwxx9/AACysrIwevToDq+fiIi6V7t7CDfrzGmoHh4eUKlUiIyMhEajQXh4OPz8/BAdHY3Y2Fj4+voiOTkZS5YsgUwmg5eXF9avX9/pF0BERN3DZCCcPXvW6Fd7fX09Ro8ebTgecOrUKZMrVyqVLU5bTUtLM/w9depUTJ06tSt1ExFRNzMZCLyPMhFR32EyEO6+++6eqoOIiMysw/dUJiKi3o2BQEREADp5lhER9Q28k1rfxEAgohZ4J7W+iV1GREQEgIFARER6DAQiIgLAQCAiIj0GAhERAWAgEBGRHgOBiIgAMBCIiEiPF6ZRt3FxdYRDf76liKwVP73UbRz692vz6lYisnzsMiIiIgAMBCIi0mMgEBERAAYCERHpMRCIiAgAA4GIiPQYCEREBIDXIRBRJ/DWmr0bA4GIOoy31uzdJO0yyszMxIwZMxAcHIx9+/a1mP/ll18iLCwMjz/+OJ5//nlcvXpVynKIiMgEyQKhvLwcKSkp+Pjjj3H48GHs378fv/76q2F+bW0tXn31VezZswf//ve/MXToUOzatUuqcqgbubg6Qi53afGPiKybZF1GeXl5CAgIgJubGwAgJCQEWVlZWL58OQBAo9EgMTERHh4eAIChQ4ciMzNTqnKoG3HMIqLeSbI9hIqKCsjlcsO0QqFAeXm5YfrOO+/EtGnTAAD19fXYs2cPpk6dKlU5RETUDsn2EHQ6HWQymWFaCGE03aympgbLli3DsGHDMHv27E5tY9Ag59uu09qxq4YshaW/Fy29PksgWSB4enri5MmThmm1Wg2FQmG0TEVFBRYtWoSAgACsXbu209uorKyFTiduu1ZrJZe7QK3u+XM7+MGi1pjjvdhR5vqsWBobG5nJH9KSdRkFBgYiPz8fVVVVqKurQ05ODoKCggzztVotli5diunTpyM+Pr7VvQciIuo5ku0heHh4QKVSITIyEhqNBuHh4fDz80N0dDRiY2NRVlaG//73v9BqtcjOzgYA+Pj4YOPGjVKVREREJkh6YZpSqYRSqTR6LC0tDQDg6+uL4uJiKTdPRESdwLGMiIgIAIeuoDa4uDrCoT/fHkR9CT/x1Kq2Lj4DeAEaUW/FLiMiIgLAQCAiIj0GAhERAWAgEBGRHgOBiIgA8CwjIuoGbd1aE+DtNa0JA4GIbltbt9YEeHtNa8IuIyIiAsBAICIiPXYZ9XEcooKImvGboI/j/ZGJqBm7jIiICAADgYiI9BgIREQEgIFARER6DAQiIgLAQCAiIj0GAhERAeB1CEQksbYGvuOgd5aHgUBEkmpr4DsOemd52GVEREQAGAhERKQnaZdRZmYmdu/ejaamJkRFRWHevHmtLrdq1SoEBATgiSeekLKcPo2D2JGl4bEFyyPZN0R5eTlSUlKQkZEBe3t7REREwN/fH15eXkbLJCYmIj8/HwEBAVKVQuAgdmR5eGzB8kjWZZSXl4eAgAC4ubnByckJISEhyMrKMlomMzMTU6ZMwfTp06Uqg4iIOkiyPYSKigrI5XLDtEKhQEFBgdEyixcvBgD8+OOPXdrGoEHOXS+wl2jrPrZE1kyK9zU/K+2TLBB0Oh1kMplhWghhNN0dKitrodOJbl2nNZHLXaBWd2znmh8GsiYdfV93VGc+K72ZjY3M5A9pybqMPD09oVarDdNqtRoKhUKqzRER0W2SLBACAwORn5+Pqqoq1NXVIScnB0FBQVJtjoiIbpNkgeDh4QGVSoXIyEjMmjULoaGh8PPzQ3R0NAoLC6XaLBERdZGkJ6YrlUoolUqjx9LS0lost3nzZinLICKiDuCVSr0ML0Ajoq7iN0cvwwvQiKirGAhEZFE4pIX5MBCIyKK0NaTFwc2hDAqJMRCIyCpw7CPpcfhrIiICwEAgIiI9BgIREQFgIBARkR4DgYiIADAQiIhIj6edWikXV0cAvM8BEXUfBoKV4hAVRNTd2GVEREQAuIdg8Th6KRH1FH7TWDh2DRF1za0/ppqPt3Hso7YxEIioVzL1Y4pjH7WOgUBEVq2t4bKp8xgIFoLHCoi6xtQoqNQ5/AayEDxWQETmxtNOiYgIAPcQehy7hojMi7fobBu/mXoYu4aIzIu36GwbA4GICLxFJ8BAICIyqS91MUkaCJmZmdi9ezeampoQFRWFefPmGc3/+eefER8fj2vXruGhhx7C+vXr0a+fdWVUW8cEGhq16G9va4aKiKg79aUuJsm+fcvLy5GSkoKMjAzY29sjIiIC/v7+8PLyMiwTFxeHDRs2YOTIkVi7di0OHDiAuXPnSlWSJEwdE+CxAqLeqzd2MUkWCHl5eQgICICbmxsAICQkBFlZWVi+fDkA4Pfff0d9fT1GjhwJAHjiiSewc+fOTgWCjY2sS7U5Ozugfxtn+rT1y97UL37FnY596nFLrIltYbmPW2JNUj5u6srphoYm1NbWtzqvJ7T3nSkTQggpNpyamorr169DpVIBAD799FMUFBQgOTkZAHD69Gm8/vrr+OSTTwAAFy5cQExMDLKzs6Uoh4iI2iHZhWk6nQ4y2f+nkRDCaLq9+URE1LMkCwRPT0+o1WrDtFqthkKhaHP+lStXjOYTEVHPkiwQAgMDkZ+fj6qqKtTV1SEnJwdBQUGG+XfffTf69++PH3/8EQBw5MgRo/lERNSzJDuGANw47TQ1NRUajQbh4eGIjo5GdHQ0YmNj4evri+LiYiQkJKC2thbe3t547bXXYG9vL1U5RERkgqSBQERE1oOjnRIREQAGAhER6TEQiIgIAAOBiIj0GAhWpLa2FqGhobh06RIA4OOPP8bMmTMxY8YMbNmyBX3p/IBb26LZ3r17sWDBAjNVZT63tscrr7yC4OBghIWFISwsDEePHjVzhT3r1vY4ffo05syZg5kzZ2LFihVobGw0c4WWiYFgJX766Sc8/fTTOH/+PACgtLQU77//Pj799FNkZmbi9OnTyM3NNW+RPeTWtmj266+/Ys+ePeYpyoxaa4+ioiLs3bsXR44cwZEjRzBt2jTzFdjDbm2P2tpavPDCC0hKSsLnn38OAEhPTzdjhZaLgWAlDhw4gMTERMPV3IMHD8bnn38OJycnVFdXo7a2Fq6urmausmfc2hYA0NjYiHXr1iE2NtaMlZnHre1RV1eHP/74A2vXroVSqcTOnTuh0+nMXGXPubU9cnNzMXLkSAwbNgwAkJCQ0KcCsjOs6+YDfdjGjRtbPGZnZ4cDBw5gy5Yt8PPzM7zhe7vW2mLbtm148skncc8995ihIvO6tT2uXLmCgIAAJCYmwsXFBUuWLEF6ejrmzJljpgp71q3tceHCBTg5OUGlUqGkpASjR4/GmjVrzFSdZeMegpWbM2cOTpw4AXd3d7z11lvmLscscnNzcfnyZTz55JPmLsUiDB48GG+//TYUCgUcHR2xYMECHD9+3NxlmY1Wq8V3332HFStWICMjA3V1dX2ya7EjGAhW6vLly4ZxoPr164eZM2fi7NmzZq7KPD777DP88ssvCAsLQ0JCAoqKivDSSy+ZuyyzOXv2rNEw8kIIq7sTYXdyd3fHiBEjMHjwYNja2mL69OkoKCgwd1kWiYFgpWpqahAXF4fq6moIIZCdnY0xY8aYuyyzeO211/DFF1/gyJEj2LBhA3x8fLBjxw5zl2U2Qghs2rQJV69ehUajwf79+/t0n/mECRNw5swZXL58GQDw9ddfw9vb28xVWaa++7PByg0ZMgQxMTGIiIiAra0tHnroITz77LPmLosswLBhwxATE4Onn34aTU1NCA4ORmhoqLnLMpu77roLSUlJWLp0KRoaGvDggw9i9erV5i7LInFwOyIiAsAuIyIi0mMgEBERAAYCERHpMRCIiAgAA4GIiPQYCNRhb775Jg4fPmzuMlqVkZGBSZMmYdGiReYupVVDhw5FVVVVl5//1VdfYcOGDe0ut3v3bkyaNAmvvPJKm8tcunQJo0aNAgDs2rULSUlJrS6n1WqxZMkSXLlypWtFS6S5rsrKSnOX0uvwOgTqsBdffNHcJbTp8OHDUKlUCAsLM3cpkpgyZQqmTJnS7nLp6enYunUrHnroodve5j//+U88/PDDcHd3v+11dSdbW1ssXrwY69evx86dO81dTq/CQOgjVq5cCW9vbyxcuBDAjXsp/PDDD9i+fTs2bdqEn376CdeuXYMQAhs2bMCYMWOwZs0a/PXXXygtLcWkSZNQWVmJBx54AIsWLUJ6ejr2798PjUaDq1evIjo6GnPnzkVGRgaOHj0KGxsbXLhwAQ4ODtiyZQv+8Y9/QK1WIzExESUlJbCxsUFERAQiIyNRU1ODjRs34ty5c9BoNBg3bhxWrVrVYriFmpoarF+/HsXFxZDJZJg4cSJWrFiB119/HYWFhbh06RL+/PNPPPPMM4bn6HQ6k6/P2dkZZ8+eRVlZGYYOHYotW7ZgwIAB8PX1RUxMDHJzc1FRUYHFixcbXl92djZSU1MBwGj6f//7H5KSknDt2jWo1WoMGzYMO3bsQP/+/dv8f/Hx8cGUKVNQXFyMrVu34uzZs222a/N2FixYgJEjR+LUqVO4fPkyxo0bh+TkZKxYsQLl5eWIj4/Hiy++iL/97W9444030NjYCLVajcDAQGzatKlD75e6ujp88MEHyMzMBAD89ttviI+PR2NjI4QQCA8Px7x581BTU4OEhAQUFxdDoVDgrrvuwuDBg/HCCy+YXP/XX3+NHTt2QKfTwcnJCevXr4ezszOioqIwfvx4FBUVQavVIjY2Fvv370dJSQl8fHywfft22NjYYOzYsUhMTMTPP/+MBx98sEOviTpAUJ+Qn58vQkNDDdPh4eEiNzdXnDp1SrzwwgtCq9UKIYRITU0VS5YsEUIIsXr1ahEVFWV4zurVq8V7770namtrxZw5c0RVVZUQQojTp0+LkSNHCiGEOHjwoBgzZoy4fPmyEEKIpKQksWrVKiGEEMuWLRNbtmwRQghRXV0tZs6cKc6fPy/WrFkjPvzwQyGEEE1NTeLll18We/bsafEaVq1aJZKTk4VOpxMNDQ1i4cKFIjU1VQghxPz588UXX3zR4jntvb6nnnpKNDQ0iMbGRjFr1iyRnp4uhBBiyJAh4qOPPhJCCFFYWCh8fHxEfX29OHjwoIiJiTGs/+bpzZs3i8OHDwshhGhsbBShoaEiKyvLsL7KysoW9Q0ZMkQcOnRICCHabdfm7cyfP1/ExsYKrVYrampqxIQJE0R+fr4QQojJkyeLgoICIYQQKpVKfP/994Z1+/v7i8LCQlFaWmpY786dO8X69etb1HXs2DExf/58w/Qrr7xiaOuKigrx0ksvCa1WK5KTk0VcXJzQ6XRCrVaLiRMnip07d7ZY383UarUYM2aMOHPmjBBCiOzsbLFo0SJRWloqhgwZIr788kshhBDr1q0TkydPFjU1NaK+vl6MHz9e/Pjjj4b1JCcnizfffNPktqhzuIfQR/j7+6OhoQGFhYVwdHREVVUVxo0bB5lMhjvuuAP/+te/UFpaihMnTmDAgAGG57U2PtKAAQPw7rvv4vjx4zh//jyKi4tx/fp1w3xvb294enoCAIYPH264W1deXh7i4uIAAC4uLvjss88AAN988w0KCwsNNy2pr69v9TV8++23+OSTTyCTyWBvb4+IiAh88MEHiImJafN1jxo1yuTrmzhxIuzt7QHcGA7k6tWrhnnNXTTe3t5obGw0eo2tiYuLQ25uLtLS0nD+/HlUVFS0+xwAhu6d9tr1ZpMnT4aNjQ2cnZ1x3333GdXdbPPmzfj222/x7rvvoqSkBA0NDbh+/Trc3NzaramkpAT33nuvYXratGlYvXo1CgoKMG7cOCQkJMDGxgbff/894uPjIZPJ4O7ujpCQkHbXferUKTzwwAMYPnw4ACA4OBjBwcG4dOkS7Ozs8OijjwIA7r33XowaNQrOzs4AAIVCYfQ677nnHvz000/tbo86jgeV+wiZTIbw8HAcOXIEBw8eRHh4OGQyGb755hssWbIEwI0vwKefftroeU5OTi3WVVZWhlmzZuH333/HmDFjWows6uDgYLRdoR8dpV+/fpDJZIZ5paWlqK2thU6nw5tvvmm4u9enn36KdevWtdiuTqczer5Op0NTU5PJ193e62urVgCGrp7mbQohWiyj0WgMf69YsQIHDhzA3XffjWeeeQbe3t4duq1pcxu3164drbvZ/Pnzcfz4cdx///1YtmwZFApFh2+zKpPJjG6qM3nyZGRnZ2P69On4+eefoVQqUVZWhv79+xut087Ort1129raGv0/CiFQXFxseP7N80ytr1+/frCx4VdYd2Jr9iGzZ8/GsWPHkJ2djSeeeALAjXsJTJ48GXPnzoWPjw++/PJLaLVak+spKirCwIED8fzzz2PChAn4+uuvAaDd540bNw4HDx4EcON4QFRUFM6fP48JEybg/fffhxACjY2NeO6557B3794Wz58wYQL27t1rWO7AgQMIDAw0uc2uvD5TBg4ciF9++QUNDQ3QaDRGw0x/9913WLZsGWbMmAHgxq0cO7OtrrZra6qrq1FYWIiXX34ZwcHBKCsrw8WLFzt857S///3vKC0tNUyvXLkS//nPfzBz5kwkJibC2dkZFy9exKRJk3DgwAFotVrU1NTgq6++anfdI0aMwG+//YZffvkFwI0zqJr3HDvj0qVLuP/++zv9PGobu4z6ELlcjuHDh6OpqQkeHh4AgIiICKxcuRJKpRJNTU0YP348cnJyTH5xjB8/Hunp6Xjssccgk8nw8MMPY+DAgbhw4YLJ7a9btw6vvvoqlEolhBBYsmQJfHx8EB8fj40bN0KpVEKj0SAwMBCLFy9u8fyEhARs2LDBsNzEiROxdOlSk9vsyuszZfz48Rg7diymT58OuVwOf39/w30oVCoVli1bBicnJzg7O2Ps2LG4ePFip9bdlXZtjaurK2JiYjB79mw4OTnBw8MDo0ePxoULFzB48OB2nx8YGIj4+HhUV1fD1dUVzz//POLj47F//37Y2tpi6tSpGDt2LEaMGMOCMzYAAAC1SURBVIFNmzbh8ccfh6urK+RyuWEdn3zyCYqKilrcwczd3R1bt27F6tWrodVq4ezsjJSUlE6/xtzc3D49zLkUONopEbXq3Xffha2tLaKjozv8nKSkJNx5553tnmV0u06cOIF9+/bxtNNuxi4jImrVwoUL8f3330OtVpu7FCNarRbvvfceEhISzF1Kr8M9BCIiAsA9BCIi0mMgEBERAAYCERHpMRCIiAgAA4GIiPQYCEREBAD4PxJXDpJ0yiLtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate 10,000 bootstrap replicates of the variance: bs_replicates\n",
    "bs_replicates = draw_bs_reps(rainfall, np.var, size=10000)\n",
    "\n",
    "# Put the variance in units of square centimeters\n",
    "bs_replicates /= 100\n",
    "\n",
    "# Make a histogram of the results\n",
    "_ = plt.hist(bs_replicates, bins=50, density=True)\n",
    "_ = plt.xlabel('variance of annual rainfall (sq. cm)')\n",
    "_ = plt.ylabel('PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confidence interval on the rate of no-hitters\n",
    "Consider again the inter-no-hitter intervals for the modern era of baseball. Generate 10,000 bootstrap replicates of the optimal parameter $\\tau$. Plot a histogram of your replicates and report a 95% confidence interval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "nohitter_times = np.array([ 843, 1613, 1101,  215,  684,  814,  278,  324,  161,  219,  545,\n",
    "        715,  966,  624,   29,  450,  107,   20,   91, 1325,  124, 1468,\n",
    "        104, 1309,  429,   62, 1878, 1104,  123,  251,   93,  188,  983,\n",
    "        166,   96,  702,   23,  524,   26,  299,   59,   39,   12,    2,\n",
    "        308, 1114,  813,  887,  645, 2088,   42, 2090,   11,  886, 1665,\n",
    "       1084, 2900, 2432,  750, 4021, 1070, 1765, 1322,   26,  548, 1525,\n",
    "         77, 2181, 2752,  127, 2147,  211,   41, 1575,  151,  479,  697,\n",
    "        557, 2267,  542,  392,   73,  603,  233,  255,  528,  397, 1529,\n",
    "       1023, 1194,  462,  583,   37,  943,  996,  480, 1497,  717,  224,\n",
    "        219, 1531,  498,   44,  288,  267,  600,   52,  269, 1086,  386,\n",
    "        176, 2199,  216,   54,  675, 1243,  463,  650,  171,  327,  110,\n",
    "        774,  509,    8,  197,  136,   12, 1124,   64,  380,  811,  232,\n",
    "        192,  731,  715,  226,  605,  539, 1491,  323,  240,  179,  702,\n",
    "        156,   82, 1397,  354,  778,  603, 1001,  385,  986,  203,  149,\n",
    "        576,  445,  180, 1403,  252,  675, 1351, 2983, 1568,   45,  899,\n",
    "       3260, 1025,   31,  100, 2055, 4043,   79,  238, 3931, 2351,  595,\n",
    "        110,  215,    0,  563,  206,  660,  242,  577,  179,  157,  192,\n",
    "        192, 1848,  792, 1693,   55,  388,  225, 1134, 1172, 1555,   31,\n",
    "       1582, 1044,  378, 1687, 2915,  280,  765, 2819,  511, 1521,  745,\n",
    "       2491,  580, 2072, 6450,  578,  745, 1075, 1103, 1549, 1520,  138,\n",
    "       1202,  296,  277,  351,  391,  950,  459,   62, 1056, 1128,  139,\n",
    "        420,   87,   71,  814,  603, 1349,  162, 1027,  783,  326,  101,\n",
    "        876,  381,  905,  156,  419,  239,  119,  129,  467])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% confidence interval = [663.32300797 871.85756972] games\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df1iVdZ7/8SdwEFEww85BL3Rad6phVzQtZlbJjlc75VGSNNamTYp2nMHGtnWjhk3T1aEmrYmkvBwZtauuK0c2CB0Yij1Yuf1QaAL7oU1UY9uUSR0OSAkIdODc3z/6erYjv/TWG476elzXua5z35/7PrzvD+fw4tw/PneYYRgGIiIipyh8qAsQEZGzkwJERERMUYCIiIgpChARETFFASIiIqYoQERExBRLA6S8vJzU1FRmz57N9u3be7TX1dWRnp6Oy+Vi5cqVdHV1AVBfX09GRgZz5sxh6dKltLW1AfD111+TlZXFDTfcwMKFC6mrq7OyfBER6YdlAeLxeMjPz6ewsJDS0lKKioo4ePBg0DI5OTmsXr2ayspKDMOguLgYgNzcXBYtWoTb7SYpKYlNmzYB8PTTT3PZZZfxxz/+kTvvvJMHHnjAqvJFRGQAlgVIVVUV06dPZ/To0YwYMQKXy4Xb7Q60Hz58mI6ODqZOnQpAeno6brcbn89HTU0NLpcraD6A3+8PfBtpb29n+PDhVpUvIiIDsFn1wg0NDdjt9sC0w+Fg//79fbbb7XY8Hg/Nzc3ExMRgs9mC5gMsXryYm2++mZkzZ9LW1sZTTz11SjU1N7fh9w/uhfdjxsTQ1NQ6qD/zZIVybRDa9ak2c1SbOUNVW3h4GBdeOLLPdssCxO/3ExYWFpg2DCNouq/2E5cDAtMPPvggGRkZZGZm8vbbb5Odnc0LL7zAyJF9b+B39dcRVhozJmZIfu7JCOXaILTrU23mqDZzQrE2ywJk7Nix1NbWBqa9Xi8OhyOo3ev1BqYbGxtxOBzExcXR0tJCd3c3ERERQeu9/PLLgeMe06ZNY8yYMXz88cdMmTLlpGpqamod9G8gdnssXm/LoP7MkxXKtUFo16fazFFt5gxVbeHhYf0Gl2XHQFJSUqiurubIkSO0t7eza9cunE5noD0hIYGoqCj27dsHQFlZGU6nk8jISJKTk6moqACgtLQ0sF5iYiIvvfQSAH/9619paGhg4sSJVm2CiIj0w7IAiY+PJzs7m8zMTBYsWMC8efOYMmUKWVlZHDhwAIC8vDzWrVvHnDlzOHbsGJmZmQCsWbOG4uJiUlNTqa2t5e677wbg4YcfZseOHcybN4977rmHRx55hNjYWKs2QURE+hF2Pg3nrl1YwUK5Ngjt+lSbOarNnPNuF5aIiJzbFCAiImKKAkREREyx7DRekfNB7Khohkf1/Bh1dHbRcrR9CCoSGTwKEJHTMDzKRtq9ZT3mlz82n9A8HCty5mgXloiImKJvICIW+MbXjd3e8xol7dqSc4kCRMQCwyIjtGtLznnahSUiIqYoQERExBQFiIiImKIAERERUxQgIiJiigJERERMUYCIiIgpChARETFFASIiIqYoQERExBRLhzIpLy+noKCArq4ubr/9djIyMoLa6+rqWLlyJW1tbSQnJ5Obm4vNZqO+vp6cnByampqYOHEieXl5jBw5kvT0dLq7uwHo6Ojg0KFDvPbaa1x00UVWboaIiPTCsm8gHo+H/Px8CgsLKS0tpaioiIMHDwYtk5OTw+rVq6msrMQwDIqLiwHIzc1l0aJFuN1ukpKS2LRpEwA7d+6krKyMsrIyLr/8cpYtW6bwkEEROyoauz028AB6HSxR5HxiWYBUVVUxffp0Ro8ezYgRI3C5XLjd7kD74cOH6ejoYOrUqQCkp6fjdrvx+XzU1NTgcrmC5n9XdXU1H3zwAVlZWVaVLxLk+H0/TnyInM8sC5CGhgbsdntg2uFw4PF4+my32+14PB6am5uJiYnBZrMFzf+uDRs2kJ2dTUREhFXli4jIACw7BuL3+wkLCwtMG4YRNN1X+4nLAUHTf/nLX2hubuaaa6455ZrGjIk55XXOhFDe1RHKtUHo13eq+rpPyDe+boZFnrl/iEK531SbOaFYm2UBMnbsWGprawPTXq8Xh8MR1O71egPTjY2NOBwO4uLiaGlpobu7m4iIiB7rvfTSS6SmppqqqampFb/fMLWuWXZ7LF5vaN4BIpRrg9Cq70x9ePu7T8iZ2tZQ6rcTqTZzhqq28PCwfv/xtmwXVkpKCtXV1Rw5coT29nZ27dqF0+kMtCckJBAVFcW+ffsAKCsrw+l0EhkZSXJyMhUVFQCUlpYGrffOO++QnJxsVdkiInKSLPsGEh8fT3Z2NpmZmfh8PhYuXMiUKVPIyspi2bJlTJ48mby8PFatWkVrayuTJk0iMzMTgDVr1rB8+XIKCgoYN24c69evD7zuoUOHiI+Pt6pskSGhW+DK2cjS60DS0tJIS0sLmrd169bA88TEREpKSnqsl5CQwLZt23p9zePfTETOJboFrpyNdCW6iIiYYuk3EJFQFDsqmuFRPd/62l0kcmoUIHLeOX5R4Il2PDwvJE+VFAlVChCR/6+v4xDw7bEIEQmmABEJYTo7S0KZAkQkhOnsLAllOgtLRERMUYCIiIgpChARETFFASIiIqYoQERExBQFiIiImKIAERERUxQgIiJiigJERERMUYCIiIgpChARETFFASIiIqZYGiDl5eWkpqYye/Zstm/f3qO9rq6O9PR0XC4XK1eupKurC4D6+noyMjKYM2cOS5cupa2tDYDW1lbuvfdeFixYwIIFC/jzn/9sZfkiItIPywLE4/GQn59PYWEhpaWlFBUVcfDgwaBlcnJyWL16NZWVlRiGQXFxMQC5ubksWrQIt9tNUlISmzZtAmDdunWMGzeO0tJS7rnnHn71q19ZVb6IiAzAsgCpqqpi+vTpjB49mhEjRuByuXC73YH2w4cP09HRwdSpUwFIT0/H7Xbj8/moqanB5XIFzTcMg127drFkyRIAnE4na9eutap8EREZgGUB0tDQgN1uD0w7HA48Hk+f7Xa7HY/HQ3NzMzExMdhstqD5TU1NDBs2jMLCQm6++WYyMzPp7u62qnwRERmAZTeU8vv9hIWFBaYNwwia7qv9xOUAwsLC6O7uprGxkdjYWIqKiti7dy//+q//yssvv3zSNY0ZE3MaW2ReKN9nO5Rrg9Cvbyj11zeh3G+qzZxQrM2yABk7diy1tbWBaa/Xi8PhCGr3er2B6cbGRhwOB3FxcbS0tNDd3U1ERERgvQsvvBCbzca8efMAuOqqqzh27BhNTU2MGTPmpGpqamrF7zfO0BaeHLs9Fq83NO8dF8q1gXX1heIH0Yy++iaUf6+qzZyhqi08PKzff7wt24WVkpJCdXU1R44cob29nV27duF0OgPtCQkJREVFsW/fPgDKyspwOp1ERkaSnJxMRUUFAKWlpTidToYNG0ZKSgovvPACAO+88w7R0dFceOGFVm2CSMg6fq/0Ex+xo6KHujQ5j1j2DSQ+Pp7s7GwyMzPx+XwsXLiQKVOmkJWVxbJly5g8eTJ5eXmsWrWK1tZWJk2aRGZmJgBr1qxh+fLlFBQUMG7cONavXw/AQw89xOrVqyksLMRms5Gfn094uC5lOd/FjopmeFTPt3LnN91EDYsYgoqs19+90kUGi2UBApCWlkZaWlrQvK1btwaeJyYmUlJS0mO9hIQEtm3b1mO+w+Hgd7/73ZkvVM5qw6Nsff4x1R9ZEevo33cRETFFASIiIqYoQERExBQFiIiImKIAERERUxQgIiJiigJERERMUYCIiIgpChARETHF0ivRRWRwfePrZlhkRI8BIzs6u2g52j5EVcm5SgEicg7pb4ys0BxnVs5m2oUlIiKmKEBERMQUBYiIiJiiABEREVMUICIiYooCRERETFGAiIiIKZYGSHl5OampqcyePZvt27f3aK+rqyM9PR2Xy8XKlSvp6uoCoL6+noyMDObMmcPSpUtpa2sD4M033+Qf/uEfmD9/PvPnz2fFihVWli8iIv2wLEA8Hg/5+fkUFhZSWlpKUVERBw8eDFomJyeH1atXU1lZiWEYFBcXA5Cbm8uiRYtwu90kJSWxadMmAN577z0WL15MWVkZZWVlrFu3zqryRURkAJYFSFVVFdOnT2f06NGMGDECl8uF2+0OtB8+fJiOjg6mTp0KQHp6Om63G5/PR01NDS6XK2g+wIEDB9izZw9paWn84he/4IsvvrCqfBERGYBlAdLQ0IDdbg9MOxwOPB5Pn+12ux2Px0NzczMxMTHYbLag+QCxsbHcdtttlJeXM2vWLLKzs60qX0REBmDZWFh+v5+wsLDAtGEYQdN9tZ+4HBCYfuCBBwLzbrnlFh577DFaWlqIjQ0eOK4vY8bEmNqW03XiwHahJJRrg9Cv72wSKn0ZKnX0RrWdGssCZOzYsdTW1gamvV4vDocjqN3r9QamGxsbcTgcxMXF0dLSQnd3NxEREYH1/H4/mzdvZsmSJURERATW++7zgTQ1teL3G6e5ZafGbo/F6w3NYexCuTboWV/sqGiGR2n8T7NC4Xcdyu851dZTeHhYv/94W7YLKyUlherqao4cOUJ7ezu7du3C6XQG2hMSEoiKimLfvn0AlJWV4XQ6iYyMJDk5mYqKCgBKS0txOp2Eh4fz4osvUllZGZh/+eWXM2LECKs2QULM8CgbafeW9XiIyNCwLEDi4+PJzs4mMzOTBQsWMG/ePKZMmUJWVhYHDhwAIC8vj3Xr1jFnzhyOHTtGZmYmAGvWrKG4uJjU1FRqa2u5++67AXjkkUd45plnuP7669mxYwe//vWvrSpfREQGYOn+gLS0NNLS0oLmbd26NfA8MTGRkpKSHuslJCSwbdu2HvMvvfRSnn322TNfqIiInDJdiS4iIqYoQERExBQFiIiImKIAERERU3RSvch54Btfd68XonV0dtFytH0IKpJzgQJE5DwwLDKi12tmyh+bT2heOidnA+3CEhERUxQgIiJiigJERERMUYCIiIgpChARETGl3wD585//PFh1iIjIWabfAFm1alXg+fH7kouIiMAAAWIY/3fzpRdffNHyYkRE5OzRb4CceMtZERGR4076IPqJ9ykXEZHzW79DmXz55ZeBu/599/lx3z1GIiJnH42RJaej3wDJyMjo9bmInBs0Rpacjn4D5K677jqtFy8vL6egoICuri5uv/32HiFUV1fHypUraWtrIzk5mdzcXGw2G/X19eTk5NDU1MTEiRPJy8tj5MiRgfW+/PJLbrjhBnbu3Mn48eNPq0YRETFnwGMg+/fv59577yUtLY2bbrqJlStX8tFHHw34wh6Ph/z8fAoLCyktLaWoqIiDBw8GLZOTk8Pq1auprKzEMAyKi4sByM3NZdGiRbjdbpKSkoJOIfb7/axcuRKfz3eq2yoiImdQvwFSXV3NnXfeyWWXXcYvf/lL7rrrLsaNG8fixYt58803+33hqqoqpk+fzujRoxkxYgQulwu32x1oP3z4MB0dHUydOhWA9PR03G43Pp+PmpoaXC5X0PzjnnzySVJSUrjwwgtNb7SIiJy+fndhbd68mSeffJLExMTAvFmzZuF0Onnsscf40Y9+1Oe6DQ0N2O32wLTD4WD//v19ttvtdjweD83NzcTExGCz2YLmA7z33nu88cYbPPnkk2zfvv0UN1VERM6kfgOkqakpKDyOmzJlCi0t/R9i8/v9Pa4j+e50X+0nLgffnkLc3t5Obm4uTzzxBOHh5obwGjMmxtR6p6u3s1xCRSjW9o2vm2GREUBo1ne+sKrvQ/l3qtpOTb8BEhER0WfbQBcWjh07ltra2sC01+vF4XAEtXu93sB0Y2MjDoeDuLg4Wlpa6O7uJiIiIrBebW0tTU1NLF26FPj2G8ySJUvYuHEjf/u3f9v/Vv5/TU2t+P2De0Gk3R6L1xua57OEam12e2yfZwbJ4LHivRGq7zlQbb0JDw/r9x/vk74S/VSlpKRQXV3NkSNHaG9vZ9euXTidzkB7QkICUVFR7Nu3D4CysjKcTieRkZEkJydTUVEBQGlpKU6nk6uvvprdu3dTVlZGWVkZDoeDLVu2nHR4iIjImdXvN5D//d//JS0trde2Q4cO9fvC8fHxZGdnk5mZic/nY+HChUyZMoWsrCyWLVvG5MmTycvLY9WqVbS2tjJp0iQyMzMBWLNmDcuXL6egoIBx48axfv16k5snIiJW6TdAtm7dykcffYTdbqezs5OxY8ee0ounpaX1CKCtW7cGnicmJlJSUtJjvYSEBLZt29bva+/evfuUahERkTOr3wA5dOgQGzZs4OKLL+azzz7jscceY+bMmYNVm4iIhLB+A2Tbtm2Ul5cTHx/P22+/TX5+vgJERESAk7gSPT4+HoBp06bR3NxseUEiInJ2OKWzsPo7rVdERM4vp3RFnu4JIiIix/V7DOTDDz/kiiuuCEx3dHRwxRVXBK4Wf+uttywvUEREQlO/AaL7oIuISF/6DZCEhITBqkNEQojuVCgno98AEZHzk+5UKCfD3LC2IiJy3lOAiIiIKQoQERExRQEiIiKmKEBERMQUBYiIiJiiABEREVMUICIiYoouJJQhEzsqmuFReguKnK0s/fSWl5dTUFBAV1cXt99+OxkZGUHtdXV1rFy5kra2NpKTk8nNzcVms1FfX09OTg5NTU1MnDiRvLw8Ro4cycGDB1m1ahXHjh3jggsu4OGHH9ZwK2ex4VG2Pq92FpHQZ9kuLI/HQ35+PoWFhZSWllJUVMTBgweDlsnJyWH16tVUVlZiGAbFxcUA5ObmsmjRItxuN0lJSWzatCkw/8477+SPf/wjqamprF+/3qryRaQXx8fI6u0ROyp6qMuTQWbZN5CqqiqmT5/O6NGjAXC5XLjdbu666y4ADh8+TEdHB1OnTgUgPT2dDRs2cNNNN1FTU8Nvf/vbwPxbb72VnJwcnn76aWw2G36/n/r6ekaNGmVV+SLSi77GyAKNk3U+sixAGhoasNvtgWmHw8H+/fv7bLfb7Xg8Hpqbm4mJicFmswXNB7DZbBw9epTU1FQ6OjrYtm2bVeWLiMgALAsQv98fdAfD4zehGqj9xOUg+E6Io0aNYs+ePbz22mssXbqUl19++aRvtTtmTIzZzTktvQ2LHSpCuTY5+5zM+ymU33Oq7dRYFiBjx46ltrY2MO31enE4HEHtXq83MN3Y2IjD4SAuLo6Wlha6u7uJiIgIWq+iooK5c+cSFhaG0+mko6ODr7/+mri4uJOqqampFb/fOENbeHLs9li83tD8Yj/UtYXiB0JOz0Dvp6F+z/VHtfUUHh7W7z/elh1ET0lJobq6miNHjtDe3s6uXbtwOp2B9oSEBKKioti3bx8AZWVlOJ1OIiMjSU5OpqKiAoDS0tLAek899VTgLolvvPEGF1544UmHh4iInFmWBUh8fDzZ2dlkZmayYMEC5s2bx5QpU8jKyuLAgQMA5OXlsW7dOubMmcOxY8fIzMwEYM2aNRQXF5OamkptbS133303AA8//DBPP/008+fPZ+PGjWzYsMGq8kVEZACWXgeSlpZGWlpa0LytW7cGnicmJlJSUtJjvYSEhF4PkF9yySX813/915kvVERETpmGMhEREVMUICIiYooCRERETNFIdmI5DZoocm7Sp1osp0ETRc5N2oUlIiKmKEBERMQU7cISkTPi+FDvJ+ro7KLlaPsQVCRWU4CIyBnR11DvGub93KVdWCIiYooCRERETFGAiIiIKQoQERExRQEiIiKm6CwsEbHUiaf3Hn+u03vPfgoQEbGUTu89d2kXloiImKIAERERUywNkPLyclJTU5k9ezbbt2/v0V5XV0d6ejoul4uVK1fS1dUFQH19PRkZGcyZM4elS5fS1tYGwMcff0xGRgbz58/n5ptvpq6uzsryRUSkH5YFiMfjIT8/n8LCQkpLSykqKuLgwYNBy+Tk5LB69WoqKysxDIPi4mIAcnNzWbRoEW63m6SkJDZt2gTAqlWryMrKoqysjLvvvpv77rvPqvJFRGQAlgVIVVUV06dPZ/To0YwYMQKXy4Xb7Q60Hz58mI6ODqZOnQpAeno6brcbn89HTU0NLpcraD7ATTfdxNVXXw3AD37wA7744guryhcRkQFYFiANDQ3Y7fbAtMPhwOPx9Nlut9vxeDw0NzcTExODzWYLmg/fhklERAQAGzZs4Nprr7WqfBERGYBlp/H6/X7CwsIC04ZhBE331X7ickCP5X7zm9/w7rvv8swzz5xSTWPGxJzqZpwRvQ1xHSpCuTY594Xa+y/U6vmuUKzNsgAZO3YstbW1gWmv14vD4Qhq93q9genGxkYcDgdxcXG0tLTQ3d1NRERE0HpdXV3cd999eDwennnmGWJjT61Dm5pa8fuN09yyU2O3x+L1hubZ7oNVWyi+8SU0hNJnQ5/VnsLDw/r9x9uyXVgpKSlUV1dz5MgR2tvb2bVrF06nM9CekJBAVFQU+/btA6CsrAyn00lkZCTJyclUVFQAUFpaGljvkUceobW1laeeeuqUw0NERM4sy76BxMfHk52dTWZmJj6fj4ULFzJlyhSysrJYtmwZkydPJi8vj1WrVtHa2sqkSZPIzMwEYM2aNSxfvpyCggLGjRvH+vXrOXLkCNu3b2f8+PHcdNNNgZ9TVtbzClcREbGepUOZpKWlkZaWFjRv69atgeeJiYmUlJT0WC8hIYFt27b1mP/++++f+SJFRMQUjYUlZ0zsqGiGR+ktJXK+0KddzpjhUbY+B80TkXOPxsISERFTFCAiImKKAkRERExRgIiIiCk6iC4iQ+LEW90ep1vdnj0UICIyJHSr27OfdmGJiIgpChARETFFu7BEJKTo2MjZQwEiIiFFx0bOHgoQOWUa80pEQAEiJmjMKxEBHUQXERGTFCAiImKKAkREREyxNEDKy8tJTU1l9uzZbN++vUd7XV0d6enpuFwuVq5cSVdXFwD19fVkZGQwZ84cli5dSltbW9B6zz33HMuXL7eydOHbg+V2e2yPh4gIWBggHo+H/Px8CgsLKS0tpaioiIMHDwYtk5OTw+rVq6msrMQwDIqLiwHIzc1l0aJFuN1ukpKS2LRpEwCdnZ3k5eWxdu1aq8qW7zh+sPzEh4gIWBggVVVVTJ8+ndGjRzNixAhcLhdutzvQfvjwYTo6Opg6dSoA6enpuN1ufD4fNTU1uFyuoPkANTU1+P1+cnJyrCpbREROkmUB0tDQgN1uD0w7HA48Hk+f7Xa7HY/HQ3NzMzExMdhstqD5ADNnzuQ//uM/GD58uFVli4jISbLsOhC/309YWFhg2jCMoOm+2k9cDugxbdaYMTFn5HVOVSgfNwjl2kROZPX7NZQ/D6FYm2UBMnbsWGprawPTXq8Xh8MR1O71egPTjY2NOBwO4uLiaGlpobu7m4iIiB7rnY6mplb8fuOMvNbJsttj8XpDcwCGgWoLxTesnN+s/CydzZ9Vq4SHh/X7j7dlu7BSUlKorq7myJEjtLe3s2vXLpxOZ6A9ISGBqKgo9u3bB0BZWRlOp5PIyEiSk5OpqKgAoLS0NGg9ETk/HR9k8cRH7KjooS7tvGXZN5D4+Hiys7PJzMzE5/OxcOFCpkyZQlZWFsuWLWPy5Mnk5eWxatUqWltbmTRpEpmZmQCsWbOG5cuXU1BQwLhx41i/fr1VZYrIWUKDLIYeS8fCSktLIy0tLWje1q1bA88TExMpKSnpsV5CQgLbtm3r83XT09NJT08/c4WexzQwooiYpb8c57m+BkYEDY4oIv3TUCYiImKKvoGIyFlNdzAcOgqQ84SOdci5SgfXh47+opwndBMoETnTdAxERERM0TcQETkn6diI9RQgInJO0rER62kXloiImKIAERERUxQgIiJiigJERERMUYCIiIgpOgtLRKQffY3ioNOBFSBnrb7e1J3fdBM1LGIIKhI5O/R1fcg3vu5el+9vFIfz/XRgBchZqr83tYYsEelbX9eH7Hh4nm7jfIoUICFOgyCKDI7+LjyU3ukvU4jTIIgioamvXWFw/hwfsTRAysvLKSgooKuri9tvv52MjIyg9rq6OlauXElbWxvJycnk5uZis9mor68nJyeHpqYmJk6cSF5eHiNHjuTo0aP88pe/5NChQ8TFxfH4449jt9ut3IRBo28aImeXvr6xwPlzfMSy03g9Hg/5+fkUFhZSWlpKUVERBw8eDFomJyeH1atXU1lZiWEYFBcXA5Cbm8uiRYtwu90kJSWxadMmAB5//HGSk5P57//+b2666SYeeughq8ofdMe/aZz4EBEJVZYFSFVVFdOnT2f06NGMGDECl8uF2+0OtB8+fJiOjg6mTp0KQHp6Om63G5/PR01NDS6XK2g+wCuvvEJaWhoA8+bN47XXXsPn81m1CSIigyJ2VDR2e2yPR+yo6KEurV+W7TNpaGgI2r3kcDjYv39/n+12ux2Px0NzczMxMTHYbLag+SeuY7PZiImJ4ciRI8THx59UTeHhYae9XScrJmY4UVHHt+H/9pP2d5qt48Le3yxDNT8Ua1JfhO78UKxpqOb3dXykv8//z369q8e8gvt+HHidoL8jnV20tnb0+jpn0oB/Mw2LbNq0ycjPzw9MFxUVGf/5n/8ZmK6trTVuueWWwPQnn3xiuFwu48svvzScTmdgvs/nM5KSkgzDMIxJkyYZPp8v0DZz5kyjoaHBqk0QEZF+WLYLa+zYsXi93sC01+vF4XD02d7Y2IjD4SAuLo6Wlha6u7t7rOdwOGhsbASgq6uLtrY2Ro8ebdUmiIhIPywLkJSUFKqrqzly5Ajt7e3s2rULp9MZaE9ISCAqKop9+/YBUFZWhtPpJDIykuTkZCoqKgAoLS0NrDdr1ixKS0sBqKioIDk5mcjISKs2QURE+hFmGIZh1YuXl5ezefNmfD4fCxcuJCsri6ysLJYtW8bkyZP54IMPWLVqFa2trUyaNIl169YxbNgwDh8+zPLly2lqamLcuHGsX7+eCy64gK+++orly5dz6NAhYmNjycvLY/z48VaVLyIi/bA0QERE5Nyl4dxFRMQUBYiIiJiiABEREVMUICIiYooCRERETNHwr6dp9+7dbNy4ke4mVocAAApISURBVPb2dq666ipWrVpFVVUV69ato7Ozk7lz55KdnQ30PfrwYNe3YsUK9u3bR3T0t8Mw3HXXXVx33XV91m2F5557jt///veB6c8//5z58+dz7bXXDnnf9VVbe3v7kPcbfHvN1JYtWwBwOp3cd999pzyy9WDWtnHjRnbs2MGoUaMA+MlPfkJGRsagfx62bNnCjh07GDZsGKmpqSxdujRk+q232kKl3/o1tBfCn90+++wzY+bMmcYXX3xhfPPNN8Ytt9xivPLKK8asWbOMzz77zPD5fMbixYuNV155xTAMw7j++uuNt99+2zAMw1ixYoWxffv2Ialv3rx5hsfjCVq2vb29z7qt9tFHHxnXXXedUV9fHzJ9d2JtTU1NIdFvx44dM374wx8aTU1Nhs/nMxYuXGjs3bu3z/5ZsmSJ8fzzzxuGYRgbN240fvOb3wx6bXfccYfx1ltv9Vh+MH+ne/fuNebNm2e0tLQYXV1dxh133GFUVlaGRL/1VVso9NtAtAvrNLz44oukpqYyduxYIiMjyc/PJzo6mosvvpgJEyZgs9lIS0vD7Xb3OfrwYNeXmJhIfX09999/P2lpaWzYsAG/38/+/ft7rXsw/OpXvyI7O5tDhw6FTN+dWFt0dHRI9Ft3dzd+v5/29na6urro6urCZrOd8sjWg1VbVFQU7733Hps3byYtLY0HHniAzs7OQf+dvv/++8ycOZOYmBgiIiK4+uqr2bZtW0j0W2+1vfTSSyHRbwNRgJyGTz/9lO7ubn7xi18wf/58CgsLex2F2OPx9Dn68GDX19nZyfTp01m7di3FxcXU1tZSUlLSZ91Wq6qqoqOjg7lz54ZU351YW2NjY0j0W0xMDP/+7//O3LlzmTVrFgkJCURGRp7yyNaDVdsPfvAD/u7v/o6cnBz+8Ic/cPToUTZt2jTov9NJkyaxZ88evvrqKzo7O9m9ezc2my0k+q232r788suQ6LeBKEBOQ3d3N9XV1axdu5aioiL279/PoUOHCAv7vyGQDcMgLCwMv9/f6/zBrq+2tpbf/va3OBwOoqOjue2223j11VeHpD6AZ599lp/+9KcAfdYQCrVNmDAhJPrtgw8+YMeOHfzP//wPr7/+OuHh4ezdu7fXGnqrZbBre/bZZ9m6dSvf//73sdlsLF68eEj6bcaMGaSnp3Pbbbfx85//nCuvvJKurq6Q6LfeaouOjg6JfhuIAuQ0XHTRRcyYMYO4uDiGDx/OtddeS1VVVa+jEPc1+vBg1/eHP/yBysrKwDKGYWCz2QYcPdkK33zzDTU1NfzjP/4j0PcIzkPRdyfW9uGHH4ZEv+3Zs4cZM2YwZswYhg0bRnp6On/6059OeWTrwart1VdfpaSkJLBMX/1m9e+0tbWV2bNnU15ezrZt2xg2bBjjx48PiX7rrba4uLiQ6LeBKEBOwzXXXMOePXs4evQo3d3dvP7668yZM4dPPvkksPvo+eefx+l09jn68GDXd+2117J27Vq+/vprfD4fRUVFXHfddVx++eW91m2lDz/8kL/5m79hxIgRAH3WMBR9d2JthmGERL8lJiZSVVXFsWPHMAyD3bt386Mf/eiUR7YerNouueQSHn30UQ4dOoRhGGzfvp3rrrtu0H+nn3/+OXfeeSddXV20tLRQUlLCwoULQ6LfeqvtxhtvDIl+G9CgHrI/Bz333HPG9ddfb8yePdvIzc01uru7jaqqKiMtLc2YPXu28dBDDxl+v98wDMOoq6sz/umf/slwuVzGPffcY3R2dg5Jfb///e+NuXPnGtddd53x6KOPBpbtq26rvPDCC8bdd98dNC9U+q632kKl3zZv3my4XC5j3rx5xooVK4yOjo4+++fzzz83br31VmPu3LnG4sWLja+++mrQa3O73YH34PLlywO1DfbvdOPGjcbcuXON2bNnG4WFhf3WMNj91lttodJv/dFovCIiYop2YYmIiCkKEBERMUUBIiIipihARETEFAWIiIiYogARERFTFCAiImKKAkTkDHn00UfZs2fPUJfRQ2trKz//+c/p6OgY6lLkHKMAEenH3LlzmTZtGklJSSQlJTFt2jSmTZvGxx9/HLTcO++8w8cff8zMmTOHqNK+xcTEMG/ePJ544omhLkXOMboSXeQk3H///UyYMIGlS5f22v6zn/2MW2+9lWuuuQb49g5zJSUljBw5kuTkZF5++WV2796N3+9n7dq1vPvuu7S1tWEYBr/+9a/p6upi/fr1jBs3jk8++YTo6GiWLFnCtm3b+OSTT5g9ezb3338/8O1dJgsKCvD5fAwfPpz77ruPadOm0dbWxooVK/j0008JDw9n0qRJPPDAA4SHh9PZ2cmPf/xjSktLueiiiwat3+Tcpm8gIifhww8/5NJLL+217ejRo+zbt4+rrroKgNdff52dO3dSUlLCzp07aWtrCyz77rvv0tDQQFFRERUVFdx4441s3boVgAMHDrBkyRLKysqIiYlhy5YtbN68mZ07d1JYWIjH4+Gvf/0r+fn5bNmyhdLSUh588EH+7d/+jWPHjvHiiy/S1tZGWVlZYCTXQ4cOARAVFUVSUhKvvvqqld0k5xndE11kAH6/n48//pjLLrus1/ZPP/0Uu93OsGHDAHj11VeZM2dO4F7WGRkZvPHGGwBMmzaNCy64gGeffZZDhw7xpz/9KXCv7fHjx/P3f//3AHzve98jNjY2MLT3yJEj+frrr6mpqaGhoYF/+Zd/Cfz8sLAwPvvsM6688kry8/O57bbbSElJ4fbbb+fiiy8OLDd+/Hg++eSTM94/cv7SNxCRAdTX1+P3+5kwYUKv7cdvenWczWbju3uGIyIiAs9feeUV7rjjDgB+/OMfc8sttwTajgfQd1/nRH6/nxkzZlBWVhZ4FBcXc+mllzJhwgRefPFFlixZQmtrKz/96U/ZvXt3YN3IyMigWkROlwJEZACtra1ER0fj8/l6bf/e975HU1MTnZ2dAMyaNYtdu3bR0tICEHRjoL1793LNNdewaNEikpKSeOmllwI3LjoZM2bMYO/evYGD+K+++io33HADHR0dFBYWsmLFCmbOnElOTg4zZ87k/fffD6z7+eefM3HixFPefpG+KEBEBvD973+fxMREfvjDH/Y4+wpg1KhRXHnllYHdVDNmzOAnP/kJN998M+np6bS0tBAdHQ3AP//zP/Pmm2+SlpbGjTfeyIQJE/j888+DvsH055JLLuGBBx7gnnvu4YYbbuCJJ56goKCAkSNHsmDBArq7u0lNTQ383Ntuuw349g6L77zzTuAOiyJngs7CEjkD3nrrLX73u9+xZcsWDhw4wNtvv01mZiYATz/9NO+++y6PP/74kNW3c+dO/vKXv3DfffcNWQ1y7tFBdJEz4IorrmDixIm89tprXHHFFWzdupXi4mLCwsIYN24cDz744JDV1tbWxvPPP8/GjRuHrAY5N+kbiIiImKJjICIiYooCRERETFGAiIiIKQoQERExRQEiIiKmKEBERMQUBYiIiJiiABEREVP+H6w135SnlHcLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Draw bootstrap replicates of the mean no-hitter time (equal to tau): bs_replicates\n",
    "bs_replicates = draw_bs_reps(nohitter_times, np.mean, 10000)\n",
    "\n",
    "# Compute the 95% confidence interval: conf_int\n",
    "conf_int = np.percentile(bs_replicates, [2.5, 97.5])\n",
    "\n",
    "# Print the confidence interval\n",
    "print('95% confidence interval =', conf_int, 'games')\n",
    "\n",
    "# Plot the histogram of the replicates\n",
    "_ = plt.hist(bs_replicates, bins=50, density=True)\n",
    "_ = plt.xlabel(r'$\\tau$ (games)')\n",
    "_ = plt.ylabel('PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pairs bootstrap\n",
    "- Nonparametric inference\n",
    "    - Make no assumptions about the model or probability distribution underlying the data\n",
    "- Pairs bootstrap for linear regression\n",
    "    - Resample data in pairs\n",
    "    - Compute slope and intercept from resampled data\n",
    "    - Each slope and intercept is a bootstrap replicate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A function to do pairs bootstrap\n",
    "As discussed in the video, pairs bootstrap involves resampling pairs of data. Each collection of pairs fit with a line, in this case using ```np.polyfit()```. We do this again and again, getting bootstrap replicates of the parameter values. To have a useful tool for doing pairs bootstrap, you will write a function to perform pairs bootstrap on a set of ```x,y``` data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_bs_pairs_linreg(x, y, size=1):\n",
    "    \"\"\"Perform pairs bootstrap for linear regression.\"\"\"\n",
    "\n",
    "    # Set up array of indices to sample from: inds\n",
    "    inds = np.arange(len(x))\n",
    "\n",
    "    # Initialize replicates: bs_slope_reps, bs_intercept_reps\n",
    "    bs_slope_reps = np.empty(size)\n",
    "    bs_intercept_reps = np.empty(size)\n",
    "\n",
    "    # Generate replicates\n",
    "    for i in range(size):\n",
    "        bs_inds = np.random.choice(inds, size=len(inds))\n",
    "        bs_x, bs_y = x[bs_inds], y[bs_inds]\n",
    "        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, deg=1)\n",
    "\n",
    "    return bs_slope_reps, bs_intercept_reps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pairs bootstrap of literacy/fertility data\n",
    "Using the function you just wrote, perform pairs bootstrap to plot a histogram describing the estimate of the slope from the illiteracy/fertility data. Also report the 95% confidence interval of the slope. The data is available to you in the NumPy arrays ```illiteracy``` and ```fertility```.\n",
    "\n",
    "As a reminder, ```draw_bs_pairs_linreg()``` has a function signature of ```draw_bs_pairs_linreg(x, y, size=1)```, and it returns two values: ```bs_slope_reps``` and ```bs_intercept_reps```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./dataset/female_literacy_fertility.csv')\n",
    "fertility = np.array(df['fertility'])\n",
    "illiteracy = np.array(100 - df['female literacy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.04409522 0.05545012]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEJCAYAAABc/7oDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de1hUdf4H8PdwFRfKzZ0RU3/u1qZsllJSpkuQlQKDiI224b1VIto0s6JNJV1bXe3i4rOPi9bmVo+3DW9cXBrX1iJNV4xar5SlQuYFhovJRW4z398fPpwFYUYGvnPmAO/X8/A8zJnDOW++nMNnzvec8z06IYQAERGRBB7uDkBERF0HiwoREUnDokJERNKwqBARkTQsKkREJA2LChERScOiQkRE0ni5O4As5eVVsNk6dstN797+KC2tlJRILi1nA7Sdj9naR8vZAG3n6wzZPDx0+OlPfyJ9+V2mqNhsosNFpXE5WqXlbIC28zFb+2g5G6DtfN01G7u/iIhIGhYVIiKShkWFiIikYVEhIiJpWFSIiEgaFhUiIpKGRYWIiKTpMvepEGlVwE1+6OHbclerqW1wev6KK1el5yOSiUWFyMV6+Hoh5sWMFtOzVsU6PX+F9HREcrH7i4iIpGFRISIiaVhUiIhIGhYVIiKShkWFiIikYVEhIiJpWFSIiEgaFhUiIpLGpUWlsrIS48aNww8//NBs+saNGzF9+nTl9YULFzB16lRERkbimWeeQVVVlStjERGRi7isqBw5cgSTJ09GQUFBs+nfffcd3nnnnWbTli5diilTpsBsNuOuu+5Camqqq2IREZELuayopKWlYcmSJTAYDMq0uro6LF68GM8995wyrb6+HocPH0ZERAQAwGQywWw2uyoWERG5kMvG/lq+fHmLaatWrcLEiRPRv39/ZVp5eTn8/f3h5XUtil6vR1FRkdPr693bv/1hm9DrA6QsxxW0nA3Qdr6ukk3t30PL7QZoO193zabagJKff/45Ll68iAULFuDQoUPKdCEEdDpds3mvf90WpaWVsNlEhzLq9QGwWLQ5ZJ+WswHaziczm70RhAH7owjfaAe+Ppuj+dVsYy3/TQFt5+sM2Tw8dNI+jDelWlHZtWsXvv32W8TGxqK6uholJSV4/vnn8eabb6KiogJWqxWenp6wWCzNusyItMTeCMIARxEmAlQsKitWrFC+P3ToENasWYPVq1cDAEJCQpCdnY2YmBikp6cjLCxMrVhERCSRJu5TWbJkCdLS0mA0GvHFF1/g+eefd3ckIiJqB5cfqezdu7fFtBEjRmDEiBHK6379+mHDhg2ujkLUZo7OnchSV2+Fj7enpk/oEjmLT34kaoWzT2tsDx9vT5evg0htmuj+IiKiroFFhYiIpGFRISIiaVhUiIhIGhYVIiKShkWFiIikYVEhIiJpWFSIiEgaFhUiIpKGRYWIiKRhUSEiImlYVIiISBoWFSIikoZFhYiIpGFRISIiaVhUiIhIGj6ki6iTqKu3tvqUyJraBlRcueqGREQtubyoVFZWIi4uDuvWrUP//v3x4YcfYsOGDdDpdLjrrruwdOlS+Pj4ID8/H4sWLUJVVRVCQkKwdOlSeHmx5hE1cvSkyAo35CFqjUu7v44cOYLJkyejoKAAAHD27FmsX78e//jHP5CZmQmbzYbNmzcDAJKSkrB48WLs3r0bQgikpaW5MhoREbmAS4tKWloalixZAoPBAADw8fHBkiVL4O/vD51Oh0GDBuHChQs4f/48ampqEBwcDAAwmUwwm82ujEYkXWP31PVfRN2JS/uXli9f3ux1v3790K9fPwBAWVkZNm3ahBUrVqC4uBh6vV6ZT6/Xo6ioyJXRiKRz1D1F1F245aRFUVER4uPjMXHiRIwYMQJ5eXnQ6XTK+0KIZq/bondvfynZtPzJUsvZAG3n03I2GVz1+2m93bScr7tmU72onD59GvHx8Zg+fTpmzZoFAAgMDITFYlHmKSkpUbrM2qq0tBI2m+hQNr0+ABaLNk95ajkboO187cmm5X8IrXFF22v5bwpoO19nyObhoZP2YbwpVe9TqaysxOzZszFv3jyloADXusV8fX2Rl5cHAMjIyEBYWJia0YiISAJVj1S2bduGkpISvPfee3jvvfcAAA8//DDmzZuHt956C8nJyaisrMSQIUMwY8YMNaMREZEEqhSVvXv3AgCefPJJPPnkk63OExQUhG3btqkRh4iIXITDtBARkTQsKkREJA2LChERScOiQkRE0nDERurWAm7yQw/f7rUb2PudOdoxydC99iai6/Tw9ep2Q6s4+p21ebsedSbs/iIiImlYVIiISBp2f1G30B3PnRC5A/cy6ha647kTIndg9xcREUnDIxUiAvC/J1c2pdcHoLbOCl8fzxbz8xJkag2LChEBcPzkSl6CTG3F7i8iIpKGRYWIiKRhUSEiImlYVIiISBoWFSIikoZXfxF1cq1dCgzwkl9yD5cXlcrKSsTFxWHdunXo378/Dhw4gBUrVqC2thZRUVGYP38+ACA/Px+LFi1CVVUVQkJCsHTpUnh5seYR3YijS4F5yS+pzaXdX0eOHMHkyZNRUFAAAKipqcHChQuRmpqK7OxsHD9+HDk5OQCApKQkLF68GLt374YQAmlpaa6MRkRELuDSopKWloYlS5bAYDAAAI4ePYqBAwdiwIAB8PLyQkxMDMxmM86fP4+amhoEBwcDAEwmE8xmsyujERGRC7i0f2n58uXNXhcXF0Ov1yuvDQYDioqKWkzX6/UoKipyZTSiLs/euRYiV1L1pIXNZoNOp1NeCyGg0+nsTndG797+UjJqeSfUcjZA+/m6G0fnWmRx99/c3et3pLtmU7WoBAYGwmKxKK8tFgsMBkOL6SUlJUqXWVuVllbCZhMdyqfXB8Bi0eapTS1nA7SdT8s7d2fnzr+51rc5rWfz8NBJ+zDelKr3qQwbNgxnz55FYWEhrFYrdu3ahbCwMPTr1w++vr7Iy8sDAGRkZCAsLEzNaEREJIGqRyq+vr5YuXIl5s6di9raWoSHhyMyMhIA8NZbbyE5ORmVlZUYMmQIZsyYoWY0IiKSQJWisnfvXuX7kSNHIjMzs8U8QUFB2LZtmxpxiIjIRThMCxERScNb1qlLCbjJDz18uVkTuQv3PupSevh6ufwyWiKyj91fREQkDYsKERFJw6JCRETSsKgQEZE0LCpERCQNiwoREUnDokJERNKwqBARkTQsKkREJA2LChERScOiQkRE0rCoEBGRNA6LyokTJ9TKQUREXYDDopKcnKx8n5qa6vIwRETUuTksKkII5fs9e/a4PAwREXVuDouKTqdTvm9aYIiIiFrT5hP1TQtMR2VkZCA6OhrR0dF4/fXXAQD5+fkwmUyIiIjAokWL0NDQIG19RESkDodPfrx06RKWLVvW4vtGTc+5tNXVq1exfPlymM1m3HTTTZg8eTIOHDiAP/3pT1i2bBmCg4OxcOFCpKWlYcqUKU4vn4iI3MfhkcrUqVPRq1cv9OrVq9n3jV/tYbVaYbPZcPXqVTQ0NKChoQFeXl6oqalBcHAwAMBkMsFsNrdr+URE5D4Oj1TmzJkjfYX+/v6YN28eoqKi4Ofnh/vuuw/e3t7Q6/XKPHq9HkVFRdLXTUREruWwqADA0aNH8cEHH+DUqVPo0aMHBg0ahJkzZ2LQoEHtWuHXX3+N7du345NPPkFAQABeeuklfP755y0uCnD2HE7v3v7tynM9vT5AynJcQcvZANfkq6u3wsfbs83TSV3u3ibdvX5Hums2h0Xl4MGDSEpKwvTp0zF+/HgAwLFjxzBr1iz8+c9/xv333+/0Cvfv34+RI0eid+/eAK51da1fvx4Wi0WZp6SkBAaDwanllpZWwmbr2BVqen0ALJaKDi3DVbScDXBdPr0+ADEvZrSYnrUqttX1aXlH7orcuU1qeZ/oDNk8PHTSPow35bCovP3223j33XcRFBSkTAsPD0dYWBhWrVrVrqISFBSEN998E9XV1fDz88PevXtx//33Y/fu3cjLy8Pw4cORkZGBsLAw538bIiJyK4dFpbS0tFlBaTR06FBUVLSvCoeGhuLkyZMwmUzw9vbG3XffjYSEBIwZMwbJycmorKzEkCFDMGPGjHYtn4iI3MdhUfH0tN9n3ZGbIRMSEpCQkNBsWlBQELZt29buZRIRkfu1+Y56IiKiG3F4pHLmzBnExMS0+t65c+dcEoiIiDovh0Xlb3/7G06dOgW9Xo/a2loEBgaqlYu6mYCb/NDD94ZXuBORxjnci8+dO4e//OUvGDhwIL7//nusWrUKoaGhamWjbqSHr5fdS4eJqPNwWFQ2bNiArKws9OnTB1999RVSUlJYVIiIyK4bjlLcp08fAMA999yD8vJylwciIqLOy6mrvxxdYkxERNTm56kAvMSYiIgcc3hO5ZtvvsG9996rvK6pqcG9996rDPj45ZdfujwgERF1Hg6LCp9LT1pVV2/l4JFEGuSwqPTr10+tHERO8fH25CXIRBrk1DkVIiIiR3gLMxG1i70uyJraBlRcueqGRKQFLCpE1C6OuiC1+XgqUgO7v4iISBoWFSIikoZFhYiIpGFRISIiaVhUiIhIGrcUlb1798JkMiEqKgrLli0DABw4cAAxMTEYO3YsUlJS3BGLiIg6SPWicu7cOSxZsgSpqanIzMzEyZMnkZOTg4ULFyI1NRXZ2dk4fvw4cnJy1I5GREQdpHpR2bNnD4xGIwIDA+Ht7Y2UlBT4+flh4MCBGDBgALy8vBATEwOz2ax2NCIi6iDVb34sLCyEt7c3EhMTcfHiRTz00EO44447oNfrlXkMBgOKioqcWm7v3v5S8ml5kEItZwO0n4/Uo9a2oOVtrrtmU72oWK1WfPHFF9iwYQN69uyJZ555Bj169Gj2rJbGofWdUVpaCZtNdCibXh8Ai0Wb9wJrORvQ8Xxa3gHJeWpsq1reJzpDNg8PnbQP402pXlR+9rOfYeTIkbjlllsAAI8++ijMZnOzp0paLBYYDAa1oxERUQepfk5l9OjR2L9/P65cuQKr1Yp9+/YhMjISZ8+eRWFhIaxWK3bt2oWwsDC1oxERUQepfqQybNgwxMfHY8qUKaivr8evf/1rTJ48Gbfddhvmzp2L2tpahIeHIzIyUu1oRETUQW4ZpXjSpEmYNGlSs2kjR45EZmamO+IQEZEkvKOeiIikYVEhIiJpWFSIiEgaFhUiIpKGRYWIiKRhUSEiImlYVIiISBoWFSIikoZFhYiIpHHLHfXUfQXc5IcevtzsiLoq7t2kqh6+Xoh5MaPF9KxVsW5IQ0SysfuLiIikYVEhIiJp2P1FLsFzJ91XXb3V7pM8a2obUHHlqsqJSE3c68kleO6k+/Lx9mz1bw9c+/tr8yG7JAu7v4iISBoWFSIikoZFhYiIpHFrUXn99dfxyiuvAADy8/NhMpkQERGBRYsWoaGhwZ3RiIioHdxWVA4ePIidO3cqr5OSkrB48WLs3r0bQgikpaW5KxoREbWTW4rK5cuXkZKSgsTERADA+fPnUVNTg+DgYACAyWSC2Wx2RzQicqHGy42v/wq4yc/d0UgSt1xSvHjxYsyfPx8XL14EABQXF0Ov1yvv6/V6FBUVuSMaEbmQvcuNealx16F6Udm6dSv69u2LkSNHYseOHQAAm80GnU6nzCOEaPa6LXr39peSz95NW1qg5WyA9vORtrVn+9HyNtdds6leVLKzs2GxWBAbG4sff/wR1dXV0Ol0sFgsyjwlJSUwGAxOLbe0tBI2m+hQNr0+ABaLNj8vaTkb0DKflnco0iZnt28t7xOdIZuHh07ah/GmVC8q7733nvL9jh07kJubixUrVmDcuHHIy8vD8OHDkZGRgbCwMLWjERFRB2lmmJa33noLycnJqKysxJAhQzBjxgx3RyIiIie5taiYTCaYTCYAQFBQELZt2+bOOERE1EG8o56IiKRhUSEiImlYVIiISBoWFSIikkYzV39R59T0CY+8N4Vks/cE0ZpaDjirVSwq1CF8wiO5ErevzofdX0REJA2LChERScOiQkRE0rCoEBGRNCwqREQkDa/+IiK3a3wiJHV+LCpE5HaOnghJnQu7v4iISBoWFSIikobdX3RD9obKICK6Hv9T0A3ZGyoDYJ83ETXH7i8iIpKGRypE1OnU1Vvh4+3Z4jLkmtoGVFy56qZUBLipqKxZswYfffQRACA8PBwvv/wyDhw4gBUrVqC2thZRUVGYP3++O6IRUSfg6BLkCjfkof9RvfvrwIED2L9/P3bu3In09HScOHECu3btwsKFC5Gamors7GwcP34cOTk5akcjIqIOUr2o6PV6vPLKK/Dx8YG3tzduv/12FBQUYODAgRgwYAC8vLwQExMDs9msdjQiIuog1YvKHXfcgeDgYABAQUEBPvroI+h0Ouj1emUeg8GAoqIitaMREVEHue1E/bfffounn34aL7/8Mjw9PVFQUKC8J4SATqdzanm9e/tLyaXl8Ye0nI1IK7Syn2glR2tcmc0tRSUvLw/PPfccFi5ciOjoaOTm5sJisSjvWywWGAwGp5ZZWloJm010KJdeHwCLRZun+dyZTcs7B9H1tLAPd4b/JR4eOmkfxptSvahcvHgRzz77LFJSUjBy5EgAwLBhw3D27FkUFhaif//+2LVrFyZOnKh2NCLqouyNCsFLkOVTvaisX78etbW1WLlypTItLi4OK1euxNy5c1FbW4vw8HBERkaqHY2Iuih7o0LwEmT5VC8qycnJSE5ObvW9zMxMldMQEZFMHKaFiIik4TAtXRj7kam74RMk3Y9FpQtjPzJ1N3yCpPux+4uIiKThkUoX4OxDtOx1EdTWWeHr4ykzGpGm2dsX2EXcfiwqXYCjbq7WOOoiYNcBdScc7Vg+dn8REZE0LCpERCQNiwoREUnDokJERNKwqBARkTQsKkREJA0vKVYBh0sh6lx4/0r7saiogMOlEHUuvH+l/dj9RURE0vBIxY2cPcR2djgWIlJHa/umXh/QLbvL+B/KjZw9xHZ2OBYiUge7uP+H3V9ERCQNiwoREUmjqe6vrKwsrF27Fg0NDZg5cyamTp3q8nVe3xfa9ByHu/pD+fQ6Im1ydt90NL+9R004O11r5200U1SKioqQkpKCHTt2wMfHB3FxcRgxYgR++ctfunS99vpCAff1h/LpdUTa5Oy+aW/+xp9x9hEUneG8jWaKyoEDB/DAAw+gV69eAICIiAiYzWbMmTOnTT/v4aFr97oNP/VzyXLbso7OPl2LmdgW2p2uxUxdoS2c/T/l4aGT9r/tejohhHDJkp309ttvo7q6GvPnzwcAbN26FUePHsUf//hHNycjIqK20syJepvNBp3uf5VTCNHsNRERaZ9mikpgYCAsFovy2mKxwGAwuDERERE5SzNFZdSoUTh48CDKyspw9epV/Otf/0JYWJi7YxERkRM0c6K+T58+mD9/PmbMmIH6+npMmjQJQ4cOdXcsIiJygmZO1BMRUeenme4vIiLq/FhUiIhIGhYVIiKShkWFiIik6XJFJSsrC0ajEWPHjsWmTZtavJ+fnw+TyYSIiAgsWrQIDQ0Nzd4/efIk7rrrLuV1XV0dkpKSEBUVhcceewynT58GcO3mzNdffx2RkZEwGo3Iy8tTPVtxcTFmz56N2NhYPPbYYzh48CAAoL6+Hvfeey9iY2OVL6vVqmq28+fP45577lHWP3v2bAD22/NGZOdLTExUssXExGDw4ME4duyYqm23c+dOhIaGKutJSUkBAFy4cAFTp05FZGQknnnmGVRVVQEArly5goSEBERFRWHq1KnN7utSK1teXh4mTZqE2NhYzJw5E+fPnwcA5ObmYsSIEcr8CxYsUD2bs+2pZr7S0tJm29TDDz+Me+65R/W2Ky4uRkJCAiZMmIC4uDj88MMPAOxvW+3aX0UXcunSJTF69GhRXl4uqqqqRExMjPj222+bzRMdHS2++uorIYQQCxYsEJs2bVLeq66uFnFxcWLQoEHKtHfffVe8+uqrQgghcnNzxeOPPy6EEOKjjz4STz31lLBareLMmTNizJgxor6+XtVsL774oti4caMQQojTp0+LUaNGiYaGBnHs2DExa9Yst7ab2WxW2q0pe+2pdr6mVq9eLZKTk4UQQtW2e+2110RWVlaLZSYkJIhdu3YJIYRYs2aNeOONN4QQQixdulS8/fbbQgghdu7cKebNm6d6ttGjR4v8/HwhhBBbt24ViYmJQggh1q9fL9atW+cwj6uzOdueaudrZLVaxbRp00RmZqYQQt22mzlzpti8ebMQQojNmzcr25C9bas9+2uXOlJpOihlz549lUEpG50/fx41NTUIDg4GAJhMpmbvr1y5EjNnzmy2zE8//RTjx48HANx3330oKyvDhQsXkJOTA6PRCA8PD/ziF79A37598dVXX6mabcyYMRg3bhwAYODAgaitrUV1dTWOHTuGsrIymEwm/OY3v0Fubq7q7Xbs2DGcOnUKsbGxmDFjBr755huH7al2vkZnzpxBeno6fv/73yu51Wq7Y8eOYefOnYiJicFLL72EH3/8EfX19Th8+DAiIiJazP/pp58iJiYGADBu3Dh89tlnqK+vVy1bXV0d5s2bh6CgIADA4MGDcfHiRWX+/fv3IyYmBomJicp0tbK1pz3Vztdo+/bt8PPzU/6WarVdWVkZvv76a8TFxQEAJk6ciOeffx6A/W2rPftrlyoqxcXF0Ov1ymuDwYCioiK77+v1euX9f//736ipqUFkZKTDZer1ely6dAnFxcXNhpFpnK5mtoiICNx8880AgPXr1+NXv/oVAgICoNPp8Mgjj+DDDz/EH/7wB8yfPx9lZWWqZvP19cX48eOxc+dOzJ49G88++yzq6urstqcjrsjXKDU1FbNnz4a/vz8AqNp2er0ev/vd75CZmYm+ffvitddeQ3l5Ofz9/eHl5dVi/qbL8vLygr+/v6rZfHx8EBt7bYh3m82GNWvW4NFHHwUABAQEYPr06cjKykJ4eLgyMKxa2drTnmrnAwCr1Yp169bhxRdfVKap1Xbnzp3DrbfeipUrV2LixIl47rnn4O3t3eJnmm5b7dlfu1RRudGglPbet1gsWLt2LV599dUWy7x+GUIIeHh4tLosDw/7zemKbI3ef/99fPjhh3jjjTcAAHFxcZgzZw68vb1x5513YujQofjyyy9VzTZ37lxMmTIFHh4eCA8PR8+ePXHmzBm77emIq9ruxx9/xOeff47HH39cmaZW2wHAX//6VwwfPhw6nQ7x8fHYt29fqwOp2htY1VXbnL1sjerq6vDSSy+hoaEBTz/9NADgtddew9ixYwEAkydPxnfffYeKCvtP+XBFto62p6vzAcC+ffvw85//HIMHD1amqdV2DQ0NOHnyJB544AFs374djzzyCF555ZVW19G4bbVnf+1SReVGg1Je/35JSQkMBgM+/fRTXL58GVOnTlU+icXGxqKyshJ9+vRBcXFxi58JDAxsdbqa2QDgjTfewNatW7Fp0yb07dsXAJCeno7vv/9eWZYQQvlEola2DRs2oLy8vFkGLy8vu+3piKvaLicnB2FhYfD19VV+Vq22q6iowPvvv99sPZ6enrjllltQUVGhXBzQdHkGgwElJSUAgIaGBlRVVSnPH1IjGwBUVVUhPj4eDQ0NWLt2Lby9vWGz2bB27doWFzQ0/owa2drTno64ou0A4OOPP4bRaFReq9l2er0eP/nJTzB69GgA17q5jh49CsD+ttWe/bVLFZUbDUrZr18/+Pr6KldqZWRkICwsDI8//jg+/vhjZGRkICMjQ3nP398f4eHhyrQvvvgCvr6+uPXWWxEWFoasrCxYrVYUFhaioKAAd999t6rZ3n//fRw6dAhbtmxBYGCgsqxvvvkGf//73wFcO2eQn5+P4cOHq5rt8OHD2LZtG4BrV7fYbDbcdtttdtvTEVfkA4D//ve/CAkJabYutdquZ8+eePfdd3HkyBEAwMaNGzFmzBh4e3sjJCQE2dnZAK4VucblhYeHIz09HQCQnZ2NkJAQhwVPdjYASEpKwsCBA7F69Wr4+PgAADw8PLBnzx7s3r1byTxs2DD07NlTtWztaU9HXNF2QMttTs22+7//+z8EBgYiJycHAPDJJ59gyJAhAOxvW+3ZX7vU1V9CCJGZmSmio6PF2LFjxTvvvCOEECI+Pl4cPXpUCCFEfn6+mDhxooiIiBAvvPCCqK2tbbGMplcJ1dTUiJdfflkYjUYxYcIEcfz4cSGEEDabTaxcuVIYjUZhNBrFvn37VM1ms9lESEiIeOihh8T48eOVr0uXLomKigoxd+5cER0dLcaNGycOHjyoertdunRJPPnkkyI6OlqYTCbliiF77al2vsafz8nJaTZNzbY7fPiwmDBhgoiMjBSJiYniypUrQgghfvjhBzFt2jQRFRUlZs2aJS5fviyEEKK8vFw8/fTTwmg0iieeeEKcO3dO1WwnTpwQgwYNEkajUdne4uPjhRBCnDp1SjzxxBPCaDSKadOmiQsXLqjebs62p9r5hBBi6NChoqamptl61Gy706dPi2nTpono6GjxxBNPiLNnzwoh7G9b7dlfOaAkERFJ06W6v4iIyL1YVIiISBoWFSIikoZFhYiIpGFRISIiaVhUiDrg0KFDyvhrRMSiQkREEnm5OwBRZ1FVVYUFCxagsLAQHh4eGDJkCKKjo5X3KyoqsHTpUnz99dfQ6XR48MEH8cILL8DLywt33nknnnrqKezbtw/V1dV44YUXlPGetm7dii1btsBms6FXr1549dVXcfvtt7vr1yTqEB6pELXRnj17UFVVhYyMDGUImsaHHAHAsmXL0KtXL2RlZWH79u3NhnyxWq3w8/PDjh07sHr1aixcuBBlZWXIzc1Feno6Nm3ahPT0dMTHx2POnDlu+f2IZOCRClEbDR8+HCkpKZg+fTpGjRqFmTNnNht6/rPPPsOWLVug0+ng4+ODuLg4fPDBB0hISAAATJs2DQAQFBSEQYMG4fDhwzhy5AgKCwuVZ1wA157Cd/nyZYeDRRJpFYsKURsNGDAAe/bswaFDh/Cf//wHv/3tb5s9K+P6IcdtNluzxxo3HXnWZrPB09MTNpsNsQSeEUcAAAEQSURBVLGxSEpKUqYXFxcrz8kh6mzY/UXURps3b8aCBQsQGhqKpKQkhIaG4uTJk8r7oaGh2LhxI4QQqKurQ1paGkaNGqW83zgK7IkTJ3D27Fncd999CA0NxT//+U9lePEtW7bYfUolUWfAIxWiNpowYQJyc3NhNBrh5+eHvn37YvDgwcqjXJOTk7Fs2TLExMSgvr4eDz74IBITE5Wf//LLL5GWlgabzYaUlBTcfPPNCA0NxVNPPYVZs2ZBp9PB398fa9asadODpIi0iKMUE6lg8ODBOHjwIG655RZ3RyFyKXZ/ERGRNDxSISIiaXikQkRE0rCoEBGRNCwqREQkDYsKERFJw6JCRETSsKgQEZE0/w/YBpd9LGx4VQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate replicates of slope and intercept using pairs bootstrap\n",
    "bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(illiteracy, fertility, 1000)\n",
    "\n",
    "# Compute and print 95% CI for slope\n",
    "print(np.percentile(bs_slope_reps, [2.5, 97.5]))\n",
    "\n",
    "# Plot the histogram\n",
    "_ = plt.hist(bs_slope_reps, bins=50, density=True)\n",
    "_ = plt.xlabel('slope')\n",
    "_ = plt.ylabel('PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting bootstrap regressions\n",
    "A nice way to visualize the variability we might expect in a linear regression is to plot the line you would get from each bootstrap replicate of the slope and intercept. Do this for the first 100 of your bootstrap replicates of the slope and intercept.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEJCAYAAAB4yveGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9aYwc13n3+zunqpfp7pmenpleZuXsC0mJlrVYix3JtmzZThw7iZ3ERhYgyIcECQIHSF7nOnY2IDf3RS4QIIiB92MAB9l8ncRx4jh2JMvWvnEnhzPkcPZ9unt6X6rqnPvhkKIsyRJJcUiKrB8gQFJPz6kqUf8+/T/P83+E1lrj4+Pj43NbIG/0Bfj4+Pj4XD980ffx8fG5jfBF38fHx+c2whd9Hx8fn9sIX/R9fHx8biN80ffx8fG5jfBF38fHx+c2wr7RF/DjyOcrKPWjLQSdnTGy2fINuqKbg9v9Gdzu9w/+M/Dv/83vX0pBIhF92/fftKKvlH6D6F/897c7t/szuN3vH/xn4N//1d+/b+/4+Pj4vJtwnHf0dl/0fXx8fN4FiHIJOXcO8vl39HtuWnvHx8fHxwdoNJAb64h6DRUKg+e9I+X2Rd/Hx8fnZsTzkFubsFuASAvUGkitIZmEfO2qf60v+j4+Pj43E1ojslnk9hY6FALPQRSbeBMT0NIC9juTbV/0fXx8fG4SRLGA3NxAC4kWIAu7uEMj0N5udv7Hj0IyDt1DV72GL/o+Pj4+N5paDbm+hnBdtBDIwi5efz9uKgNaI8+exTryMjSa8IlH39FSvuj7+Pj43CgcB7m5gSgVUYEgslxCd3bhvvceEAK5tor10vOIYgnV040KtkChAJ29V73knon+17/+df7u7/7u1X9eWVnhU5/6FH/0R3+0V0v6+Pj4vCXnVgvMLOWZGEgw2hu/cReiFGJ7G5HdgVAQ6g2EtHAP3QWWBbks9osvIDfWUakU3ugY9tEjWFtb8HM//Y6W3jPR/+xnP8tnP/tZAM6ePctv/dZv8du//dt7tZyPj4/PW3JutcBf/sMRXE9hW5Lf/9xdN0T4xW4eubGOtgPQbCKUh7f/AIRCUCphH34Za/48KpHAm5zCOnUK+5WXUB2dqKkpU73zDrgu9s6f/Mmf8Lu/+7t0dHRcj+V8fHx83sDMUh7XU2gNnqeYWcpfX9GvVJBrq6AVOC6i3sAbHoG2NqjXsV56AWv6NDocxp2cQs7NETjxHVRrAm90DB2OQDQKsdg7uow9F/1nn32Wer3Oxz/+8b1eysfHx+fHMjGQwLYknqewLMnEQOL6LNxsItfXEdUyGpC1Gt6+fejOJLgu8sQxrJPHQYEaGoHNdQKPfxcViuAOjiAiLeiWKLqzE5nLwvY2RDuv+nKE1npPk4t+53d+h49+9KP81E/91F4u4+Pj4/O2nFnIcWJuhztGupgc3GPnwfNgYwOyWePTl0rQ3Q379oFSMDsLr7wC1SoMDJjXjx41P9vba3b04TCk0+bwdm4O1tfhvvvg05++6svaU9FvNps8/PDDPP7440QikSt6bzZbfkOSXDLZyvZ26Vpe4ruO2/0Z3O73D/4zuOnv/2Jz1eYG2raQ5TKqLY4aGQUpkavLWC+/jNjdRaUzaDT28WMIz0FlMqjWOMIOoLq7Td7O/AJsbyK1pvmBh+l84G627TdaPFIKOjvf3vrZU3tnZmaGwcHBKxZ8Hx8fn3cjolxCrqygtYJGA0EI945DEAzC9jaBw6/Axiq0d+L1dGOfOAGVKjrThYrFIRyBVBLlNLFOnkTsbCBchfPAQ3ihMPbieRgfhPTV+/p7KvrLy8tkMpm9XMLHx8fnxlOvI1dXEI06OE0EAm9i0hy8FgrYzz+HXJhDR6KofUPYxw8ji0VUMoPIZNDhFnQqiVIa6+wMcmcHWavh3nMfuq0Na3EeLSzE9BlItMLH9131pe6p6H/iE5/gE5/4xF4u4ePj43PjcF3TSVsomNp7p4k7OAydnVCrYb/4PHJ2Fm1ZeH37sGamCUyfwutMwr5RCAdR6W6UAPv8eWR2G1mu4B48hJdOYy2cR+/uwswMttOg+eDD8Mj739El+x25Pj4+PleK1ojtbePbA6Jew+vtQ/f0mmqd40eR06eQjotKZ5BLCwSf+gFeRwJv3yCEWtCpFCoYxD4/h8zlkcUc3vh+mgMD2PPn0QvziLMz2LUqjYc+gHI9gk98B3QdPvnZq750X/R9fHx8rgCxm0euraKVRjRq6I4k7sE7TEbOzBmsU8ehVEGnUuhsFvv5Z1GxKO7AAIRb0MkkqqUFe/6C2BfyeIMjNO+7D3tpAbkwj5g7i10u07j/QTwg+Pj/oBF4qTS0tr6j6/dF38fHx+dyqFaxlpeg2YBGExGN4r7nbrAs5NIi8sQxZC6Lau+AFg/r6CtgB/HS3RCNors6UbFW7IV5ZD6PLORQ3X00770Pe3kRubCAmDtHsFigce99eLZF8Mnvo6XAy2SQ21lkuQzvsMnVF30fHx+ft8JxzCFtsQjNJgQDeAcOQksLYmMD6+Rx5OoKqq0V1dqGdfoUWgh0eye6rRU6OvHaE9gLc8i5OURhF53ooPmhR7GXVrDmz8PCPMF8nuZ778YLhQn+4Em0BV46g8zvIEtl1PgYOtMNQ1cfqwy+6Pv4+Pi8OUohNjeQO9smI0cI3NExk22fz2G/9CJycR4dDODFE1hnzyBcF1pb0dEYpNJ4XV3YC/PY8/OI4i66tY3mQx/A3tzEOjcHa+vYO1s4hw7RnDpI4AdPoFF46R7k7g6yUkaNjKO7u1G9/aYiqK8LCo2rvi1f9H18fHxeh8hmkasr4HkIxzGxCekMVCrYLzyPOH8OoRQq3o41dxZZq0Eshk50olNJvFQae2kJe2kRUSxASxjnvvuxcjlzcLu9g9zawJk6gBobN2IvwOtKIwtZRLWMGh5DpzOo/gG8ySlEpYK1cB5iNoiWq743X/R9fHx8LlIuYy0toOsNhOOg0mlTbeM4WEePIGfPQL2JjrdhLc0j5+dRbTHoSqJSKXQ6jVxbxT78ClRrCAvcu+9BFopYS4uI/A7WxibuyDjuyAjBH3wfrTVeuhu5m0XUq6ihEXR3jxH7qQOIUgnr/Dm04yBWVsBWMHXXVd+iL/o+Pj4+zSZycQFRKoHThHg77tg4APLMNNbMaSgU0NE2rNIOcmkBLxpFJBKo7m50uhu5tYk8fgzKFQQK78AdiHodubSIKBax1tdwB4Zovu8B49krbQ5oi7uIegU1MmYatAaG8PYfQBQKWOdmjdivrmAvzqPb2qHnwXd0q77o+/j43L54njmkze5AvQGxqDmkDQSQiwvI6dOI7DY61IJsOlgrp9DhCCreju7O4HX3ILe3kKdOQLWCcBzcyf1Iz0UsryBqVaz1Fdzefpr3vo/g0z9ET58y7yvsIqo11OAQOpVGDQziHTiI2M1jnZ0xMQ4rK1hLC9DRifPQT0CszUQxvwN80ffx8bn9uNBcZa2voes1CAbxDhzkXFEx8+QM+2sbTJTW0HYA7XpYS7MQDKCiMVS6B93Xi8zuYJ0+ga7VEfU67ugYwpLI9VXEhThl1dlJ8+77CD7zA/TsjKnGKRehXkMN7kOnu1H7Luzs83ns2Rl0vY5cWUIsrUCqC/eRD0M4jCgVkWdnINECqYGrvnVf9H18fG4rRLGAXFyEWg3QqJERdFeKuVOL/L//OYerNN8Skj/oDDC+eh6EhQ6F0T3dqJ4+RKmAdWYaXW8iqiXU0BBeMIRc30R4DnJjA90ex7nzEMHnn8U6P4+XSSNLRXAaqH5zKKwGB40FlM1eEPsaYmEBubEOXV24jz4KARu5u4tYXkJHI7jveS/0doNz9ffvi76Pj8/tQb2OtTAPxYKpvOntw+sfgFIR+4XnmD26jauiKASuhjM7DSYsG9JpvH1DiHIZ69wM2mkiyhV0dy9efx9yYw2UQm5soqMRnIkJgkdeJrCwgEqnEOUKNJvo/gGzsx8awtt/0Ij9zDS6WkUszmNtbKIzadxHPwqAyOWQu3l0vB3nrruxclms8+egpwPk1ScX+6Lv4+Nza+O6yOUlM4jcc1BdKbyRUWg0sI68glych1KZKRnApgVXC2yhmWq3ccfvRdQqWHNn0bUa1OqQ6sJLdWNtbqCKRcT2BgSDOKOjBI++QnB5yRzQlspo14P+AXQ6jTcyaqpxtrexZ8+gyxXE4nmszW10dwbnsY8jXMeUixby6GQS5577sNbXsOfnjOc/MQltUSi7V/04fNH38fG5NdEasb6GXF8zYt0exx2/C4TAOjONPHcOCnlQCiu7w2SjwZfY4lSkh8mBBCOhDuTieag1oF5FxDvwMhkC65t45SrkdpBa4/UPEjh5jODaKl4yg6xWwHHRvb3Q3YM7Moo3ud+I/fQpI/ZL57G2c6ieXpyPPoZwHZPWWa2g0hncsXuxVlaw5+fwhobRXUnk9CkCh1+G998Hyf6rfiy+6Pv4+NxyiFwOubQEpV2IteLd9V4Ih5Hzc8izs4idHXAVsphHVMpmPm1HB0NjYww3HcTOMqJcQTTq6LY4qiuJtbWFbNTQpSKy2UBlurFnppGbG+aAtlYFz0X39kCmF3d0FG9iCrm1iT19ygxGP38eCgVUXw/Oo3cjmg3k2qrpCejpwRubwFpdwlpaxhseQbe2Yh07jHzuOVR/L+5dd0N/P9Sv/tn4ou/j43PrUKlgzZ9H7OYgEMLdfxA6OhEry1izs8jNNVNtU6kgSyVwXVSiHW9qP8JTyNU1RK2CqNfQ0RheZw/WThbRaKDrFTP6MJPBPnsOub2Fl0oh6zXQoFPd0NuDOzaBNzZxSexLJeT5OSiXUX39qLvvQdRryNVlhAJvYAAdi2EvL4Hr4U0dQEsL6/DLWLks7vAI3j33QjiE7ug0KZv1qx8X6Yu+j4/Pux/HMREFW9sINN7QCLq3D7G1hfXDJxEb64hSGSolrEoF7bjoWBTnwD1IFHJ9HVEuIxoNCIWMZ5/LIndyqEYDWcijUimstTVkdhuvM4MIhhBKQSqD7unFG5/EGxtHbm5inzmNzueR83NQr6MGBlF334Os1ZDLCwjLxh0aRlgW1uoqynVx33MXutk0g1cqFZyRMfTwCDoSRXd0gJSgtQl9ewf4ou9zzTm3WmBmKc/EQILR3viNvhyfWxmlEMvLpsHKdUwo2fAworCL/fwziLUV2C0iymWoVqFeR4dCeHfdjQgFTGNVqYxoNiEQQHV2Ind3kfkcSrlY2RxeRwJZrSKnT6OSGRO+JjUik8Lr7kNNTOINj1wS+1zOiL3TRA0Oo3r7zPsXlyAcwhubRDhNrLV1dDqFc9/9iHwO+8knQIM3No5qa4NYDNUWByHAdc1BtPIgPviOHpkv+j7XlHOrBf7yH47gegrbkvz+5+7yhd9nTxCbm1jLi+hSCVJp3IlJU5HzykuItRVkNo8o7aKrNajXwA6i3vMeCIcgm0WslxGNOiIQQLW3GvEvFVDCwt7awGtrQ1Qr2DvbqGQaAmEQGpJdxqaZmMQbGjZiP30Kncubg1/PRQ0Oo1NpRKOBXF6CSATv4B2m7HNlCdXbj/vgQ8jlZQL/9Z/ocAvexBS0thqhj0bNTdZqpmxTWuiODnQgAKEQcPW7fV/0fa4pM0t5XE+hNXieYmYp74u+z7WlWMQ6cgTyOXSiA+/+B0FKrNOnkMuLsJ1FFrLohgO1KgKJOngIHY0icjnYXEc0GwgrgBdpxaoUEa5C2Tb2+hpeJArVMnZ2G68zBcEwSAGdnSZtc2ICNTiE2NzEPn0and1BLi6glYcaGkV3dSDqDeTaGrqtDe/QXcjtbaylBbzBIdyJCeTsDIFv/iu6qwPvPe+FSAsq0XFB0E0DmSiX0eEwKp0x1s41whd9n2vKxEAC25J4nsKyJBMDiRt9ST63Co0G1txZqBcRrsB97z0QjSLn57Dm59GbW4jdLLLeQFfLSMCdOgjxNkQ+h1jeQTQdhCXQoQiUK1iuiwq3YK2uokMBdK2Knd0xYh8Kg22j29vRAwOo8QnUwD5zTjB9Gr29hVxcRFsWangY3Z5A1GrIrS3TUHXoPdjr61gL83ijY+iDd2AdO4L1wnN43b14992Pjl7w6y3LWFXZrKkYam1F9fTuyWP0Rd/nmjLaG+f3P3eX7+n7XDs8D7k4b2wSKeH+u3FlBLG8jPXyS8idLdjeNpn21QrCU3jjY+hkGpHbgcUlRLMO0gLLNj9jOXixKPbyEjJgQ7mCzFZRXUm8QBht2RCPG7GfmDRiv7lpxH5rC7m0hA4GUKOj6NY2RLWG3NpGd3XhDo1grSxjLy/jjU+go1Gsl57H2t7G6x/AffD96FgM3Z4wfr3jmEEtWqPaE+jOzj19nL7o+1xzRnvjvtj7vHO0RqytIufnkc0G3tAIat8+cMvYT/8QsbWB2NxEVMqIeh2aLt7wCF5vN2K3gFheQjTqZhctJaJWBSuI15HAXlzEymehVEE2qibSOBIGaaE729D9g6ipKVRvP2Jr06Rtrq9jra2iAwHU2Bg61oqolZG7eVQyhTc6ir20iLW+hnfwIFqD/ezTyGoZd2gUd2wCHY+jYxcGm1eryMIu2rJRnV1gXx859kXfx8fnpmNuepnZE/NM2nVGpvpwxicQu6Yih0YZcWYOK5dFex6iVkfv24e7bx+yWDDfCOpNsC0QIIsltB3E60piL8xhZTfNzt5x0B0JVCwC0ka3xdADF8Q+04PY3kZOnzJiv7GKDkfxRschFoVyCVkq4GV60bEI9soKWBbO3fdCqUTg+08gPBdnYgrV1YVqT0A4bD7ICruISgXd0oLKdJvd/nVkT0X/iSee4G/+5m+o1Wo89NBDfPnLX97L5Xx8fN7tVKvMv3CSv3yxhKsFlh3kf93RzvjLLyN2NpCra1DZRe5W0c0mdKdx3ztiSiKXlsyOPxwCKZD5XXQohJtKmx14dst0xTab6M4OlAYCQeOr7xtETUyh0hlEdgd55jRyYwW5voluieCNjENrDFEqocslVF+/OTxeX0cFUjgPvh+xukLgv/7DxDTv34+OG6tG27bx63d2EM0Gqi2Ovhq/vtlE5PMQfmePeM9Ef3l5mT/+4z/m61//Op2dnfzqr/4qP/jBD3j44Yf3akkfH593K66Lde4scnGBM7tRHAQawFXMPnucSbWKXFlCCw3ag3gcb2wcUa8hNjagVkMEAmhLYu3sQCiM29uDvbSEtbkJ9TrSaaI7OlBCQOiC2A8MoyYnUV1JRC6LnDmNXFlB7mRRsYiZUxuLQrkC1SrewBC4DnJzE9XXj/OBh5GzMwT/7V9QrVG8996DjrejEwlz/tBsIjc3jF+f6ECHuq7suWiNKJk+A20HIPDOJXvPRP973/sen/jEJ8hkMgD81V/9FaEL5Ug+Pj4+gBG1pSWss7NgW7h3381YSRBYPIGrNbbQHFg5hXDziHod3d4B73svansXsbWJqFbRgRBaYMYVtrTg9Q8gl5axTpw0pZmuaypwrFZ0IIiORdGDI6Yap6MDmc8jZ6aRq6vI3A4qEkEND6FjbVAqQrOBNzqCLFeQm+t4Q8O4k1NYR18h+I1X8Lq6cB940KzRduEsq1JBFgtoO4DqSppzhSvBcRD5PMJz0S0taDuA8FxUJHLzxjAsLi4SCAT4jd/4DdbX13nkkUf4whe+sFfL+fj4vMsQW1tYZ04jmg28iUlUdy9yaZHJ+fP8QXqXM6tF9js7jDey6NY23DsPITwFW1uI9S1UKIwQErm5gWgJ4w0MIjZWsY4dNQe4ykPHE+ZwNhhEt0bR+4aN2CcSyN1ds/7aGqK4i4614Q0OQrQVKhVQHt7EJDKbxVpbM5U4+w9gP/c08vln8Hr7cR9+xNTXt7SYD7DdvPkgikSuyq8XpaKxkOwAuqUFUS5Bo3mpMetaPHettb4mv+l1fPnLX+bIkSN87WtfIxKJ8Ju/+Zt88pOf5Gd/9mf3YjkfH593C8UiHDsGuRyMjsL4OGxswPQ0ZLNw6pSJTHAcs6s9cMBkzhSL5q9QyPzz+roR2/5+WFmBtTVoNMwaHR2mGkZKiMdhZAQmJyGRgHze/OzKChQKZuZsVxdEImaaVmsr9PWZa1IK7rwTgkF44gnz82NjMDxs3hMIgOeZ63YcaG+/1E17ubiueRbNpllbKfOh09JirvcaNmbBHu70u7q6eOCBB+jo6ADg0Ucf5fjx45ct+tlsGaV+9PMomWxle/vqv9bcCtzuz+BK7/9WzAF61/4ZaDRMjv3yIqqvH+/uhxCFAtZ/fBdyOcS5GeztbRTC7NzH9htbY3UbUSwYH14r2rdXKGoL1TeA2FhHfv9pRLOOFkBbwuyIXYUOh1CjY+jhMVRbHLmTR5ycQayvQ6UKsQg6kYZoFFGooa0WvHQ/9voaurqIe8cdUK0R+Lf/RDQauFMH0JPvuRR+tllA5nMghNntB1ugqqB6mf9tymVkqWgiFlpbEeUGIle5EMMQBwVkK29424/77y+loLMz9rbL7pnof/CDH+SLX/wixWKRaDTKU089xYc//OG9Ws7H5w34OUA3CUohz5/DmplGdXThPPJhhNPEPnYEtrcQc+ex1pYRgA5FUKNDEI6iq1VzCBoIogMBrLVldDACBybQc4smY6dWAynR7Ql0wOzsdSxumqZGRlGxVmRhF+vUSeTmOlRr6GgI3dcLkag5J4jGcHt7sdfXsAoFnHvvR2ysE/jOtwGBd+chdCZj/HohjFgXC+hgCJVMXZlf73nGq2820NGo+TAqFqBUuqYWzluxZ6J/6NAhfv3Xf53Pf/7zOI7DQw89xM/93M/t1XI+Pm/AzwG68YjVFawTx804wfsfQoRCWGdnEZsbiIUFrPlz5gdbWnB7+iHRDs06cn2Fs6qNU7qXA9uLjNs1vOEJxM42PP88Vr4ElkB3JMAKmb+Px83OfnAI1dpmPPvTJ5Gbm9Aw6Zq6rw9awohmE5VIoEJhrM1NREsLzkM/gTwzTeCb/4JuCePd+z4j6tHoJb++VjNi3d1zZX79hUYspETF26EmERcsnOtdq7+ndfqf+cxn+MxnPrOXS/j4/Fj8HKAbh8jnzIFqrYZ38A5UV9Jk5KyuwvIy1uy08cIjLXiZHuhKmvLGtXWQmlni/EVtHy6Cb4aTfNGeZerwS8hqFdqi6M4OsIMgMeI9Mo4eHERHIojdItbx48jtTXBcRDCI192NjkSRjoOXTJumrZ0ddE8vzk88jPXKSwT/v39EdXThPfxBU3ETDJpI460tUzkTb0cnOi7/IShldvWNujnY7exC5PMmtvlKf9c1xO/I9bll8XOAbgC1Gtaxo8jtTbyJKdS+QeTyEoFnn4KNTexTR/FqDYi24vWkIJ0xg8tX1wAPHQgiNjeZtltxgxKFwNWKsxsV9jeaRozjUag5qEQHamQEvW/QJGgWilhzxxG5LDQcRMDG605DSxQ8zwwnbzTNEJQLsQj2C89hv/ISXv8Azkc+bnJvLMvU9W+so4W8ctulVjO7eiGMuIfDxsJpOpeatW4gvuj73NS804NYPwfoOuG6yOnTWHPnUEODOI8+htzewn7macjtYB99BYpldHsc+rvxerpBK+TWuqlWCQRgO4vwNKqvn8myi93wcIXE1orJcAMdShvFSqVwU32ogQEItyCKReS5o8jdPDQaiGAIL5M21TgaVHc3slRC5rJ445O4o2MEnn8W8cxTeBMTOD/106+Gn71aMhkKo1Lpy6+cUcrEK9RqJg45mUIUdpG57FWXb+4Vvuj73LS82UFsMtl6oy/L57VobWybkydQyRTOox9FlEtYL7+EyO9gHT6CzG6j4+3Q14/X12fCz7a2Ea6DDgWR2ztopVH9fVBrYJ8+wVRhlz+MzXG6c4gpUWZcVlAdXaiRUbhjAlVTiFIJMXsGWSyB00DYAbyePpNxY1modBq5W0Dm8ngHDqIdh8DzzyFqZdyDd6J/4mETfqaUiV6u19Gx2JVFGtfrZlevtdnVx1rNUPatTZOY2bG3iZlXgy/6Pjctb3YQ+8B7+m70ZflcQGxtYR9+EW0FcN7/MEIp7BPHIbeDPHYMe2UJt70D1duP6smAHURkdxCOgw6GEOUi5Bp4vf3QdLBPnUKUCmjLRqczjAUDjKlVVFcn7sgDRozDYWg0kMdOI8tFaDoQDKDTPahwGEIhdGcnMpdHloq4dx5C7mxjP/4/CBTu3feg+veZ3+M4JqVTeVcm0BdD06pV842gK3lpwpVl3xQWzltx816Zz22PfxD7ztizHoVSCfvlFxHlMu4990Ikgn3+HGxtIU+exJ6bwYu14+0bglQSHYogdnNQ2TKNVdUqcmsbr7cPtMI6cwa5mzORxekMBIKARiVTqJERVKYHwmHE7i5i5gxYHjJXgEAA3dsLARsdjaFa27DyOXStjnPvfcjz5wn8xzfRLRHcBx80IWe2bQR6fc3Ux1+JQDebpi5fKROalokjdvPIzQ1T0XMTWThvhS/6PjctV3sQeys2ZF0pe9Kj0GhgHTmMXF/FO3gI1Z1BLswjT6wiZmexp0+iwxG8gWFUe9xMrCqUYGsLHQojXBeZ3cHL9EG8HevcLDKfhUDYlEDaAUCjUhkj9qkMIhREFvIwcwpZq5tM/K6EqbO3bHQ8gRcOYu/uAuDcdz/y1EmC3/hnVDKN+5HH0MmUsZReO4LwcgX6tYFngaDZ1XueEf9S0Vg6N6GF81b4ou9zU3OlB7F+Q5bhmvYoKIU8fQrrzDTeyAjORz6GXFsl8MLzMH8e+/DL6FAQ1dOPam+DRIfZzS8soiItCAT2+ipuOoPX24t9bg6xmwUrgOrpuxCXIFDJbtTIMCqZRoQCyHweMbsMlTqyXjERxwNDkIii7SgKsEpFRCKBMzqBffRl7G/8M2pgH86nfsYczr62bLKt7fL9esdB5HLG+mltQ/f0Isols6u/KP5XGqJ2k+CLvs8thd+QZbhW1phYXsR+6UVUKonz0Y8hc1nsl55HLK8QePE5FJhs+UgMlUpCtYo4v4AIh9DBIPbqCl5nCrenF3t+HpHbBstGZXrBlkbsu3tQg0NG7G0bmc8iZ9bQtQayVkZHInhDQ4CEziS0BWFtBwaHcMfGsF94HvulF/Cmpmj+/KDflpUAACAASURBVOdM1Y7jXIo0voIRhK/u6i9681KaD41iwRzyXmlT1k2IL/o+txT+OYDhHfco5LIEnnkaHQzgPPJhZK2KdfQwYn2d4HNPQ71uKmUiLajuXmjUkefnEKEQKhJGbqzjJRK4vf3YiwuwvY2wLVSPGT6CLVGpHhNh3JWEQACxs4O1voFuNhDVIrTE8IZHAdCpFKJaNbk2h+7Bi6ewX3wRUS3h3vke1Ec/Zso+q1Xj11/JCELXNbt6z70k7I0GMrtz3ebWXk980ffZM26Et+43ZF3iqnoUajXs555Blgo49z0A0no1fjjw/DPIQhGnrxcRiaG7M2jHRc6dQwRD6GgrYmMVEWszWTaLC+Zw1w7gdfeiLQm2herpQw0Ooju7jEWytYW1s4WuN6FagUgEb2gMISUqmUSWy4hKFXdqP6JUhCefJFCs4dzzPvTwsKmvLxaufAThawPPLjRgiVIRubZqcnXexRbOW+GLvs+ecCO9db8h6ypwXeThl7EXF3EPvQev6y6shfOwvIL94rNYm9s0+3qQySRkelBKI+cXkLbJkpGrK9ASxk11Y68uImdPI6WN19OPJ4WxdAYGUPv2oRMdCARiewtrO4tq1qFaQURiqOERsAOojnZksYhs1HEP3om1sUbwu/9l6up/8iM0g23Gr8/lEE7zVd/9bXlt4FksZj4gLvr+zQa6tfXK6vTfhfii77MnvNZbd13FN58+z6feP+yL8c2G1sjp09inTuCOjOM8+hHE8jL2M88gj7yEtbyEm0yjx0YhmUJbAcTSAhKBSsQRK+vIShU33Y29uow9P4cEvEwfnjA2jrdvEL1vH7o9YcR+axMrl0fXquhaFRGPo0bGIBhEt7YiKxWEAveue7DOzhD4z39HZzI0f/KnobMT4iHk7OJrRhBexkS+C5OskNK8Jxg0Fs7mBsDVjTJ8l+KLvs+ecNFbd12FBk7P55ldPnLbVtPcjIi1Nexnn0KlUjQ/8hhyYwPrxeexDh9Gzp9Dtxsx1p0JdDCCXFlCuA5eMoVcW0GuVvA6OwhsrCPnXzSHsplulLDAti+I/QC6rR2BRmysYxUK6GoFXW+gE3F0dwaiMXSwxVTohMK4o2PIk8cJfuubqMFBnF/4vDmcrVSQ62ugE2+wXt7USnx94NkF20cUC4idbWPhXEnUwi2CL/o+e8JFb/2bT5/n9Hweze1dTXNTUSgQeOr7YFk4D38QeaHZyjpxHGv6NCrWihocQscTqHgbcmkR0Wigkt3IrTWshQUj9ltbyGPHQCjODBxgOtjFlCowvK8L0d+HirUhPA+xvoJVLKGrVWg2zdDw7h4zJcqyTcxxPI47MoJ15BUCp07gTUzR/KVfAds23a+7+UvCnWqD1wwReYOV+LP7GY94lwLPOjuNrbOzg3CaV1a6eQvii77PnjHaG+dT7x9mdvnIbV9Nc1PQaGA98xQyn8N53wNIrbFOnkSePkHgxHG8QACvfwCdaEd1dGAtLGCvreL19iO2NrDnz+J1diGrFeyTJxFovHSas6Ek/3fr/SYcTcAXe2zGwgq5uowsldG1CjgNdEcXuq0XnWhHewrpeajeHnQoROCVF+FkE/fOO3F/8pOX/HrXMd2vPT/+z80bynTndxj94KTZwV9Iy7w43UoHg9fxgd+c+KLvs+c8dDADwIN3dPu7/BuBUuaQ9tw53Lvuwt1/AOvsWeSZUwSOHUUBXnc3tCfwupLIxUXslVW8/kGE52GfO4PX3glaY0+fRgiNl0yhrACEgpzKTOGWL8Yga84s7TIZ3GK26DHtRJmMtzDWG0bHE2jXQbguanAYrTysl15C2BbNe+83c2cv+uyXK9L1OlNRj4AUuEqbjcVEj6m3r1SuPC3zNsAX/avEb/V/e17/tfvBO7pv9CXddsjZM9iHX8EdG8f50KPIxfNYzz6Dffww1JuoZBeqvROV7EKurhI4ehR3cAglNPbcNKrVRA7bs7OgPVQqjbAD6FAQb2wM0t1MuGHsUxpXa2w0+6vrnK14/IUzgYvA3oUvdjYZQ+GNTWDlc9jPPY1ojdN89COQziDKJcTFUsm3E+mLU6wuBJ4NHxzi9zq7mFnIMRkXjFlVtHXrV+FcLb7oXwV+q//l4XfH3kA2Ngj88El0KoXzoY8g11ax/vvbyOPHTG16Zye6twOdSsL6OoGTx015ZYfGnj2DamtFI7Hm58Dz8Lq6kMGQ2TmPjqAz3aiWFmSlynh5jf+r1WW6arFflhhJRvj3ZhJ3W5jdv4LTLRlGwxWCP3wSlU7hfOrnoLXVRBqvrZrAsrcT6UbDZOY3I+hA8JLlU6sxblUZH20xg1V8C+ct8UX/KvDF7PLwu2NvAKUSgScfBwTuBx5G5LLYTzyOdewwMp+9UJEzgupMInayWMeO4fYP4CU6sc6eRbVG0baFtbgEnoPbkcQK2RCO4I2OoTMZVKgFWS4hcyuIag1KRcaFYDSdRMe7US1RJncd7KzA1WALzYG1aYhmaP7C58zhbC6HqJTfPrBM60uNVxcbpjLtsFU0B7yViglQ8y2cy8YX/avAF7PLw++OvY64LvZTTyJ2dnDueR/SaWI98xT28SNmoEeszaRfproQ+Rz29Em87j687h7s8+dRLS1o20KuriE8B6+jExkKIcIteOPj6HQGZQeRpSJyZ8ekXRbLCEuiutKoRBxtB5Ceh26NMTKS5ot6mpkdh4mxNIMf+gU8xzGZ8z9mBOGPWKapyBsCzy7eJ5ubyM3dy2/I8vkRfNG/Cnwxu3xu5e7Ym+VcRx5+GXtmGufQexFDIwSOH8U6dgS5voEOBPAGhsyYwmoNa2YGlUrhZXqxlhZQtgQEcn0NoRReoh0ZCkNLBG9sApVJISwb8nlktYKoN6FYQIRsZuL9TMs2JkOKcQt0ZxdePI514gRiYZ6ROw4xtP+A8es3N97yUPW1lmlACn7/40OMTPZdyrqvVpGFXbS0YGIfSkau70O+hfBF/yq5lcXM5+25Gc51xNxZAi++gDs0gvv+h7FOnsQ6fhhraQktJF5PD7q93ZQtzp9HJZN43T1YSwtoBCgXuV1CuC5eIoEMBMzA8vFxVCaNwDLDT6o1RLMBhSIiZKPSaWYCnfzvxYixb3YE/+uBOJPnZ8yHzN33oAb2IfJ5E372dtEGrsvM6ZVXG/lcpTlTghHLunRg+9pMnUAAqF+vx3zL4Yu+j89VsBfnOpf9zWF7m8A//xs6nsD5wMNYM2eQX/9H5MKC6ZhNZ9CJDmg0kUtLkGjH7R3Anp9Da41wXUSlbEonW+MQDkGsDW9iHJXsQiAROzloNKHZMBOroi3o7jReWztoxZlCCEeLV0V65uwWwz/1YUh0mBmxmxtvm04pyiUzhFxaTIxnsI9tv2qZTraC3Fh/2xp9nyvHF30fn6vgWp/rXNY3h2qVwP98D1qDuPc/gFxawP7mv2DNnYNGA9XRBV3t0PSQa6uIaASntw97cQ6raewbUSkiNbjRGDIYRrfGUBOTqGQS4SnEzg40PXAaZuB3JIru7caLxUFrCLegujNMLuUIrHu4WmDZkrFH7kZ6Ljqff+sRhJ73akjaq4FnQjAK/P6nJ5k9v81Ed4yRiV7U6zx/n2vDnor+L//yL5PL5bAv/AH4sz/7Mw4dOrSXS/r4XBeu9bnOW35z8DxzSLu1hXvoLtB17O98G2vmDKJWQSc60N09oEGub0FLFC+dwVpeJrC2CgpErYoE3GgUHQhBexxv7ILYO03EuhkQri+IvY7GUD096EgrQiuIRVHpNHJzE+vkCUb7B/i9zx5gZjHHVIfNSLuFSiR/fKTxmwWeganOyecQtRpjXRFGRw6+64eU3OzsmehrrVlYWOD73//+q6Lv43OzcjWHstfyXOfHfXOwjryCNX0K98BBRKwN+4Vn4NwM9nYO3R7HSw0htIbdHNIO4qUz2CvL2KvLILTx87XAjbSggyFo70CNjRmxrzcQa6sIpdCNJqKyi460onoH0C0hQEC8DS+ZRK6sYE+fxh0ZRX3gYUSpxHizwdi9vehYK/rNbsrzjDX0usAz4EfHEcbbjR3lc13YMzU+f/48AL/2a7/G7u4uP//zP88v/dIv7dVyPj5Xzc1wKPuGbw5OjuA//Duqrx/vwJ1Yp04gT51E5rKQ7jIHpUojiiWELfCSKeyVFey1FUBAvYYQEhWNoCwb3ZFEjwyjMhlEuYxcWUFrbcS3UIS2KKp7AIIBCAYhkUAlOpCLC1j5PN7+A7jDI8jCLjKfe+so4lrNWEOvDTy7yIUd/6vjCP0N4XVnz554sVjkgQce4Ctf+QqO4/Arv/IrDA0N8dBDD+3Vkj4+V8XN0mw32htnNOwSfOI7qFjMiP3sGcTj38Pa3ORMvJ/T/RPcG67Sk1tGS4nqaMfe2CSwdtJ47rUG2pbotlZAojpTMDJiqnFKJeTSAhrJbM1ipmoz2WIxuq8dkGb0YXsCWluRS0tYpRLewTtQ6W5ksYAoFn/8CEKlTKVNvW4qbV5bmvna116/4/e57git9Zt+M7vW/O3f/i1ra2t86Utfuh7L+fhcNmcWcvzh/3kG11XYtuTPf+MhJgevs91Qr8O3vw2NBoyOwuIiHD4MKysQCnEmsY8/9PbjaoktFH/eco7J3AJsb4NS0GyaUsZIxIhtVxdMTUEmA/k8ZLMmf77R4ExB84e1MZOLI+DPJ5tMjqXN60tLEIvB3XdDNArlsvmdicSbC3W9Drmc+ftEAlpaLr3WbJp1tX7jaz43jD3b6b/88ss4jsMDDzwAGI//Srz9bLaMUj/6eZRMtrL9mhzt25Hb/Rnsxf13RgP83i9eslY6o4Hr94yVwn7qSeTaGu7gCLLcwPr6vyLn5yEYMmWPGl6qh3FtgRICVwteXirRv70ETQdawmZ6lLRQ8U700AgqnUbkdhAvHzevV0tYpSqqvY0XA/24tQu5OBpeqIQpHVlkWrcxdtcdjKWiiGwZ1ZQQi4MH7JQvXbPWJgLhQuCZTiTMB03ZhXLp0uxZO2BeC9iXXrsG+P8PvPn9Syno7Iy97fv3TPRLpRJ//dd/zT/+4z/iOA7/+q//yp/+6Z/u1XLvem6W7s7blRvRbCePHsU+eRy3txfSaQIvPY+YPYOwLFQqCa4HykMC+0UZG42rFbZWHNhdAMtCJVoQSFQqhR4ZQ6XSyK1NrHMzqGgMPAexXTQHsj19oDVTYYldNkFottC0hiz+4nw7rtIE1uf5vc8cZHSo540XfDHwTCnj1be/pkz1tVOqolHfwrmJ2TPR/+AHP8ixY8f49Kc/jVKKz3/+89x11117tdy7mpvhINHnOrK4QPDZp9Ftraj+fuyZaawTx9HSQiUSSNcDpRFCgBWEapmJlaN82V7kdCjF3bJAr10HEUan0qjBIVQqg7W1jjU7jYonwHGRm9tmjXgrYEFHApVMMlqp8AeUOB3sYnyyj9mVAq7aMfOMlWZmo8LoUNJc65sFnr1mTCHNJjKfM/Nq36YZy+fmYE+Pzr/whS/whS98YS+XuCW4WQ4SffaYfJ7gE99DeR6quweWFwh89ztopfESHchmE6kFBMOAhkYVsbyMaNYQrsd4RDEqKrS2BChmMriDw3jpNPbKCtbZabxEF5bjmuiD9gS6zUJLGzq6jFhXSsjsDt6+IYYOZRh2HHQkAvE49oncj5aLNptm5/76wLMLiHIJUSyiA8E3fhD43NT49VI3AX5q5y1Oo4H9+HcRuSyqtx+xvY31H/+GqDXwOhPIag3Lc9EtYVTTQdaqyJ1NRLWGcF10LIbXHkVqD51MwR37ccOt2GurBM6cwe1KYe3uYm2uouMJ0AIdDKHTSXRHJ6JQQBZ28YZHceNxhFboSBQdjQIw2s6r5aKTXUHGRBldaL6xpPK1Fk4s5g8peZfii/5NgJ/aeYvieVjPPoN1dhbV3YuOxbH+61tYpQpuZxcWEstVJpCs6SIqNaztTUS9hnAa6LY4XjiMVEBnB+7gCF6qC7Kb2CsbeF1JI/Zry6am3nUgFEKlMuh4O6K0iyyX8UbHIRYFIUz37uuHjDgOY4EG44NBcy2x1h99/aKX71s4twS+6N8k+KmdtxbW8aNYL72I6uzA6+0l8NSTxlrp6IKIQioH3RoDpwmlIlZuF1EpIlwP1dZm8uw9DxJJ3KFhvGQX9uI89rk8jA4h1rax1tdQ8Tii2UQEQ7iDg4jWVigUkU4Td2IKEQqbXX9HxxsijV8NPLNsdCLxhnx7USqa19/My/d51+KLvo/PtWRxkeD3H4eAhTcwiP38U1grK7hdKWQshnBc3NZWhOvAbgFrNwelirFxWmN4nXGk8kwy5r5hZsJpZjbK7C9uMpzsxF5ahKUlvPY4otGAQAg1PAKRKKJUBC3wJqcQoRDEYqi2120kXNdYNK8LPHsVz/vR130L55bDF30fn2tBLkfgv7+NrBg7xXrxOUI//AFeRwe6PYGsNfDa4ohm40KUQQEqZYTjoGIxVHvKiH28DXdoCC+eYH4py//2WnFpw3YVX1qfZiSeAOlBJIoaH0e3RBClsjlQndpvdvbtCfTrG6Euxh9Iy+zqX2/x1OvGwhHCWDih0PV7dj7XFV/0fa4Lt2wfQq2G/d3/Mrv5yf3I08cJ/sPXcONxvHg7oumiYhGk5SJy21DcRZSqCK+JirWhMhmkUuhYG+7wEF6sjcDyIjSanIoO4BYwDVlCclq1MhpTcOcUquohKmV0IISeGkCHW944gvDirr3Z+LG186JYQJTLbznVyufWwhd9nz3nluxD8Dys738f+/QJvNFxdChI6J//Hrc1hteVRNRqEIuY5qpsFspl46E3m6hoFNXVbdIxYzHcgUFUayv26go4Ci/djbUwz1TUxhZtpiFLwsTBQbw2aebEtoRRAwPoaOxSR+xFLowWRMoLu/auN1y7yOUQrvP2U618bjneVvS/9rWv8TM/8zPEYm/f3uvj82bcan0I8pWXCDzzNF5vH15nJ4FvfgMVCZsI4moNhURGooidbSiXEJUysuniRSKofd0IyzI7674BdCSCtbUFQuClurEW5lCtrahYlIlGlj/o3uZ0tIepiMd42EFFOmH/MMq10K1tly7q9YFn6cwbd+2vtXDerIrH57bgbUV/ZmaGxx57jEceeYRf/MVf5I477rge1+XzLuOt7JubrQ/haq0mcW6W4He/g4rHcbp7CT3+PZQEL5XCKpVQSqGjEaydbahWjdjXm6hIC95QP0iJDgZRffvQoSBWNosXbsFLdmEtLaNaI6hwC6JRRY2M4+0bZrRRZ8xroDq78DLdJgWzP4m+mL1yUcjBRCN0vK6c8rUdtWHfwvG5zJTNcrnMt771Lb7xjW+gteZzn/scn/zkJwnt4WGPH7j25tyMz+By7Jtr5em/0/u/KqtpfZ3gt78FSuF2pQg9+QQ0amaSVLmCkhYEQ1g7G+hqFcoVpNNEt0TQ6W50KABConp60aEQslRkJjHAdCPEgdw8oy0uwnFAC9TUFGrfIKJcBq1Q6W50MmVKLi80SiW7YuycXULUaj8aePZaLlbpuA6qtc0kZ94i3Iz/D1xPrkvgWiwW42Mf+xiNRoOvfe1r/P3f/z1f/epX+cpXvsKHPvShK79qn1uKy7FvbpY+hCuymnZ3CXz7P5D5HG5vL4Fnnib84vN43b1Qt6FUxou0YG1tQqUC9TqiVkO0tKAGR9AtYdAeKtONsILmQLWrk1kV4f/JdeMi+Waogy85xxk+MIXqH0CWCohSEdXbj04mzUSpi4evF5uknOgFsX+T+OdaDbmb//FVOj63PW8r+s899xz/9E//xHPPPcdjjz3GV7/6VSYnJ1laWuLzn/+8L/o+19W+ObOQ4/njq1f9jeGyrrVWw/7ef2MvLeL09GKtrxI8/DK6uwe3K4ksF1EtLcj8Jtb6KqJRh4YDoTB6eAQ3FkM4TVRHB8K2Eer/b+/N4+sor/v/9zMz90pX0tUuW4tXWZsXecHY2MZgQ8ABTOrUJQGTlrRpoLQEWvpqwCGEtmEJUPolIfx+7Sv50uSbfAk0NBAChJ2YYMvG4E3eZFmytVq29vWuM/N8/xhL3iRbtnVlXd3n/Y915947c57xzGfOPc95zhFY2ZloXV0YlZXsSyjEFMfLJKOxp/gK8pP9iJ4e7MlTkdnZyP5Vsf1ljE8ueJadCid7eieXOo6PHzyer1Ac55yi/6//+q/cfvvtPProo3i9J5ZnT5kyha9+9asRNU4RHYxWGYmqxi6eeXkHYfPCs4DOamswiLF5I/qe3ViZWci2ZuK2f46Vm4PMzUV09aB53IieHoyGegiH0YJBZHw8cto0rOOLrmRiIhhpCMPAystDa23BOLAPyx2HDJvM1np4XesvbQwz40PYhUWOWPeHTE8ueJacckbBM8AJ4bS3Iyzz+GdUzSbFuTmn6P/FX/wFX/va107Z9pOf/IS77rqL++67L2KGKaKL0QjfHKjrwDQvPgvoDFvDYfSdOzG2fYqVkIjs6sS97TOsiTmY06ejt7Zix8Uh/H2Io41UujLYn1bKrL4jFOboWGlpiKAP4uKQ3mSEJwErJw+t+QjGnnJsoSNCAURyCuHrvsh0l5v1h1vYb3oomj2N6bOnYuv6iUnXXmex1ZA9ZH0+p5Kmpp+Zm69QnIMhRf+ll14iEAjw85//nFAoNLA9HA7z8ssvc9ddd42KgQpFP8VT0jAMDdMcoTBSOIxWVYm+ZbPTcrC7B/dnW7HS0jHzCzBaW7FDQfD7MBrqkWGLqoRMHpv1VUxNx0DyUMdGimQfeFMgJRUrNwetvhGjfAcSiQiEEFMmE1x+FZoNWmszMjOL6dctYVrq8RaE4TCirQ1hmUPnzUvppGT6/TB5gmpSorhghhR9wzCorKwkEAhQWVk5sF3XddavXz8qxikUJ1OQl8Ljd195UTF9AEwTrb4W/dMtjmfd2YVRWYGVnIJZUIB+rBm7rR36+tDbW5GW5QizYbBnylxMTccWGqa02e/NozDVhzUhG622BuPzz7GRzkrcggKCV65A6+tGO3YMOWUK5oprBzJpBmrSH28rOKjHHg47IRzbclIy09Ih3XtqTF+hOA+GFP2vfOUrfOUrX+GDDz7guuuuG02bFIohKZmWTkbiBYYzTBPtaBParh3oDfXIri6M6oPY8R7CBUUYrcegsQERCGC0tCKREAojDA2ZOwlr2nRKenowpI2JE48vmpaBPNKG8fkWMC00aWGXzie4eAlayzH0xjqsgkLsxUvB7Xbi8C0tJwqa5eQO7rH318rRjaHDPArFBTDklfTTn/6UO++8k82bN7Nly5Yz3n/44YcjaphCMWKYJqK1BW3fXoyqSmRHO9rhGqShY02fgWjrcB4Cfh9663HPXloIBHLSJKz8fERrK1pDPUUZGaxPbWC/kcmstkMUbz+ECAfAcGMuXkJw3gL0ulq0hjqs4lnIqVOdksS9vWhtrWdPpZQS0dHurKpNSFAhHEVEGFL0+zN10tJURoBibDHshV6W5Yh1bQ3a3nJEUxOisRFNAzllEvT50OobwN+D1tKGtEIghaOzObmEi0rQmlvQq6uxJ2Zjzp6NTEymuPogM49+BoEAJCUT/sIN2IWF6IcOoTc1YC1YgJyYfaLT1FkKngFnhnBOX1WrUIwgw1qReylQK3IHJ9bPQVtfmIf+Y9PZV9QeF3vR1oq+aydaVSWipQUtHMbKzQUzhNbWjujtderjmBZCgBRAdg7hollobS1ozcew8/Kwi2diezzo+yvQjjaCGUKmZhC+9npkXi56VRUyPQNzwWWQknJGwTOGWrne24vW0+2EcE5acXsuYv0aUOOP0IrcBQsWIM7y03L79u3DNFGhGDl2V7cOvaLWspwsmO4utL270ct3obW0ghXCzsrC1HS0Y01oXV1o7R0cSJrI3pwrmN1dS2GCjVkyC73lGMbendi5UzCvW4UdH49esRdXTa3jieflEf7CjeDxOG0KfT7CN90Mbrfj1R9pdAqeDeXVn1wY7Wzev0IRIYYU/TfffHM07VAohkXpjMwzV9T2i73fhzh4AKOsDK21BWwTmZKGHedBazmK3tGG1tGFDIc5kDqJR+f+uZN6KSTf7dxE8Y5tyKlTCa/+E6RmoO/djetwNZoQWIVFhK69HhEOo7e3Y6XPILRmLYRCJwqena1/bCiE1tHu9JlVIRzFJWRI0a+rq2Pp0qW89957g76fl6dqcCtGn5Jp6SdW1E5KodAVRBztgfo6XB+9j954BKkLp7l3QoIT4qmpRu/ohFAQWzPQNNibOwtTN7ARmNJmX1wW0//sK2ALjPLtaNXV4DKw5pQSvPY6tPY2tJYWrNlzMK9agejqRDt29JyVKwfSMl1u1WdWMSYYUvTfeustli5dyi9/+csz3hNCsGrVqogaplAMRUGOl0J3CGEF4Fgzxnu/xzhcg9R1p/VgUhJaZxfaob1one3IcAhpuBC6AZkZhOYtYGa7H8MyHU9fExQV5GB89gnaoWpITMJatpTwkuVox5rQm48Rnr8QMjPROjucWH9yyuAFz+DEBG4wcPa0TIXiEnDOidzy8nLmzp17yraysjKWLVs2rAM89dRTdHR08OSTT56XYWoid3Bi+hzYNm1HWvl03zGK03RmbXkPvWIflkuHOA/EeyAQQqs/jN7eCkET6TIQmoadlo41dwFay1G02sNYBcXsn3sVFb0asw99Tknl58i0NKxlKzBnzkJvPort9WJddjlC05xiZu44p4zxUN56fxVMKZ0J3Pj4iJyGmL4GUOOP2ETuvn37kFLy4IMP8u///u/0PxtM0+Rf/uVfhgz7nMzmzZt57bXXWLly5Tk/q4g+Rq3vrW0j2tupbuzi6ffrMU0Lw7Z4qKONovgERLwHiUSvPoh29CiYpuPZx7mxU1I4MP9q9vfozNm3ixk5SQTv+GsIB5n1h7eYU1uPlZND+JZ12JMnozc3IyyL8DXXIYJBtGDwnMXMUnI8xQAAIABJREFURE83oqfnRBVMFcJRjGHOWntn06ZNNDc3861vfevEFwyD66+//pw77uzs5Nlnn+Xuu++moqJiZKxVjBlGpe+tlM4EbTiEHRfPgbJyzLDXKYEgNPYn5FIgGzAOV6EdaQIrjG240OLisFNTsWfPpbpX8oQ9C9Or89qc2TyY08bs372G1nQEc1o+gb/8K0RyClpPN3gSCF/7BUQggAiHz74S9viDSISCqs+sIqoYUvQfffRRAL773e/y+OOPn/eOH3nkEe6//36ampouyLChfqZkZXkH3R5LjIVz8MrHhwibNuCkTja0+Vg6f9LI7FxKaGuDUBAmZ8Fbb8H777NYT+P1xKVOo3AhubztICmVnzuNwl0u8CRCaiqUlkJnJ5Rvo2ru6hMTtrZNbdkeFk/IgL+43fHIQyEoKYQJE5zjJieD9yznNxh0bBMC8nOHzsGPMGPhGriUqPFf+PjPuRpk27Zt573TV155hZycHJYuXcqrr756QYadLaY/amGFMchYiGdWNXbx3tbagdeaJpiUkXDxdkk54D3b3mSMjRsx/vCe06BEaOSJVh5PLOfzVovSQzuY2nKYvng3muHC9qYhS2YifX0YZVswp+Vj/vXfMrWhE6PFxETDEDB15eV0ugPIniAyPx+Zlu4UPBMecBkQAAJnjmOg5LE7zllIpWnQHQJCZ3w20oyFa+BSosYf4XaJeXl5bN++nfnz56MNsxvP73//e1paWlizZg1dXV34fD6eeOIJHnrooWF9/2yMSlhhDHOxnaNGggN1HQMPZAEsL80B4K3NNRdm10k1Z+yUVLR9+4h/503o6UUzdGwkdlIi+pFGSvbuYlqvD8vtRibEIZNTsQoLkQE/xp5dWFOn4//be6HxCPH/+z8o9fl5aMkX2Zs3i5lxAQpSE7CmlkJKMtLrRXqTh7bLsk6UUUhOViEcxbjgnKJfXV3N7bffjmEYuN1upJQIIc66IvdnP/vZwN+vvvoqW7duHRHBh1N7nJqWzesbD7FmeX5MCP9IdI4aCU5vOTgl23thD+KTxT45BXGomvgXfoJob0NqOmBheRLQ2n24d25EBP2QlICdkIhMTkROL0KaYYz9+7CmTsP/t/ciKiuJ/4/nEeEg4SuvxCqZzQwzTH6qjj19NlZKupOBc7besYGAk4UjhLPgKi5zxM7dxVLV2MWG8iYmZSTExDWvGHnOKfovvvjiaNgxbPoFp1/49x7uoLJ+R0x4/CPVOepiOb3l4Hk1G4dTGoLYySmIpibifvV/0Y40YGu6U3Y4PgHNtDHKysDnw07wIBITIC0NmTcFKSVGdSXW5En4/+5etN3leJ77X0gpCF99NTK/AMwwMjkVe3q+02Q8OWXofHkpnSyc3l5k3NkXXF0q+n/l9j9sY+GaV4w8wwrvvPPOO+zfv5+7776bDz/8kJtvvnnYB1i7di1r1669KCNPpl9wXt94iL2HneXvl1IAR5MR7xx1EZzecnBYjdFPFntvMqKnGffrr6IfqsaWIAJ+dHc8tm7g2v4Z9PZiJyQikryIpETk1OmQ4EarqEJOzsN/59+h7/wMz7P/hnTHEbzmCzB5CoSCyMwM7ILT+s4OhmU58whmeMxn4Zz8cCVGrnnFyHNO0f/JT37Cpk2bOHr0KH/5l3/J888/T21tLffcc89o2DcoBXkprFmeT2X9jnMLzThixDpHjTDDaYwuOjsQPh92khfh68P1wXvo+/YhLRN6etA1A+n2oO3dhdHZgfQkIFJSER4PMm8ytkvHqKuD/Cn4v/FNjM8+xfPjZ5CeRILX3gA5WQhbYOflYZXMQmZknn0VrN/vVMIUAjstffD69mOM08NqsXDNK0aec67IXbNmDa+88gpf/epX+e1vf0t3dze33norb7/9dkQNG86K3FjM4om2zIUBsU9MQmttRt+xHW3PbqczVG8PmmVjJSag79uHaGtBJiSAyw3xcdh5k5CeBIz6OqzsHELXrSJ91+f4N2/FTvZiLluOTE2HOANZNAtzTil4PEMb0994vK8PGR+PTE0bcyGcc1HV2EVDmy+mY/rRdg+MNBHP3umfwO0nOTkZY4y0bjs9xBCNjNcHl+jqRPT1OeGZgB9913b0nTvR2tucrJygHyspBe1wJUbrUYhPQHq9EJeAPSnveOniRqyJLvy33Ibr80/x/Mf/DzlZhL94A3ZiEqSlYc2dh11UcvZVsKbphHAsE9ubjBzDIZxzUZCXwtL5k2Ja9BQXxznVOycnhw0bNiCEIBQK8cILL6gKmyPEeEw/Fd1dTkmChAQI+NEr9qOX70RrakL0dEFvH3ZyCjS24C7f5dSe93rRXHHYkyYhk5Ixjh7Bcrvw33Qz7l278PziBez0TEKrVuGZmIHtTcNceAVyypSzG3O8mYnUdGR6+uCNxxWKGGNI0f/hD3/IP/zDP/D1r3+d559/ngMHDjB//nzmzZvHM888M5o2jlvOO+tlDHOq2AfQamrQ9+xCq6lF9HZDWxsy2YvW1Y177x6nJLHXi9B17NzJ2GkpGMeOYbndBK9ajlFRiefXv8ZOSye08lqENwV7xgxYvYqwdRbxltL5leHznb2ZiUIRo5y1icq6det49NFH+cUvfoHP50MIgedsMVPFeTEeJuZET7dTLz4+3hH7pib0in1oByqgp9spQ+zxoAUC6AcqIC4OOzERNB1r8mRkagqulhasQJDwwivQDlUR9/u3sZNTCF11FSIrC2vmHKzLF0FCAqR7YbDQxskhnHMUSFMoYpkhRf/KK69k5cqVSClZunTpwPb+xVn79+8fFQPHM8PJehmrDIi92+2IfWsLWnUVevkuZHe300rQFQ/BEMaBA0hPnDPJqutYk6cg01Ix2tqw/EHCc+ai1dXh2vAednIqweVXI3JyMedfhj2n9OwLqVQIR6E4L86ZvfO1r33tkizQUvX0B+dSn4MTnaAMRE8voq8XefQork+3QHcHRo1Tk8e2bfTaQ0hXHAJA07CmTUF6U9C6erDj4pB5eWgtLeg1NdgpiZiXLYKp0xyxLywetFF4VpaXlubugXx/6fE4WTgxFMK51NfApUaNP8LZO2NtRa7iEtHbi9bdhdR1J8e9rQ/Z3Y1r80bsY63ojbVOs29NR6urQTMMp1m4DdbkyZCWhujthT4f1qRJaJ2dGJ+WIdPSCF5/PRQUYs2Ziz2jYGgBD4fh6FG05i6nz+xQnasUCsWQjI3cS8XYpV/sNc0J4/h8SNPE9cknyCM1iLoGXN2dSHcc4kgjQtdB1xG2jTV5CqSmowV82H29WNk5iK5uXNs+Q2ZmEVyzFpmfjzV7LnLqtKFt6Os7/sAxoHgqtp44asNXKMYbSvQVg9PX58TKwWkD6PchAf3TLWiHqxB1tejNrU5Mv60dTdqAgHAYa/JUSE9H+H3g68XKzILODlw7tiGzswncejtMnYI19zJkTs7gx7dtJ4QTCCATEk5k4YyRNSIKRbSi7qBRIKoWYB0XeywLwmG0cAipG2i7d6Hv2YNsaEA/0uCEeXq70CRO8xHTxMrJg8xMRCgE3T3YE7IQbe0Ye3YhJ00lePtfIKdOx1y4ENIzBj9+KOSUM7Ytp8LlUJ9TKBQXhBL9CBM1C7B8PqecsGVC2HTaBbrdcOAArh3boK4Bo6YaqQlk0I9m2xAKg2Vj5WRDVjYiHARfHzI5GdHZibF/H+b0fEJf+0vs6fmYCxY6na1Oo6qxiwMHjzEzTWdGbvLZ2xQqFIqLYtzeWZH2roe7/7LdTae0FRxzC7D6xT4cAstGmCZ2XDyipgbXljJEfS3GwUqkZWLbNpplQdiEUBA7OxuZNQHNtJzyx/FxCL8fvaaOcGEh5k1fwp5RiLnwcifH/nRsm+qKBp55q5qwLU88FJXgKxQRI+ruruGIbaS96+Huv6qxi092n+gRrGli7CzA8vvROtqdHrHgiL3bDa3NuD/+GHG4GqOiAswQttDRLRPLDEMghMzKxM7KQkjQfAHsuDhEyIfW04M5I5/w5YuwS2ZizV84eFnjUMg5tpRUtIYI23JcrEpWKKKBqBL9ipr2YYltpMsbDHf/g7UVvOSC5vc7Rc+CIdAE2BbScEEggOvtN9AOVjkrant7sHQXugACfmTAD+lp2FOmga6hhUJOb9lwEC0YxMzPx154OfacUqw5cwddUHUix9+NnZkFuk6x1YWxtSGqVyUrFNFEVIn+7urWYYltpMsbDLb/wX6BnP65ZaVDZKqMBoGAI/Y+//EMGIk0XAhfCNeHH6DvLUevqEC0t2G745Bx8RDwgS8AySnYeZPA7YJQEGFrYFpOWuaMGVgLF2OXzsOaNfvM7BrbdiZmgwFkUtIZTUqieVWyQhGNRJXol87IHJaYDyYkIxnjP33/wKC/QIYraBGdf+gX+74+p069oSMNHWFKjA0fYuzYhrZvH1prK7bbBUlepP94rXtPIva06Qi3C2nbEA4hLNupm5Ofj7VoCdaChdgFhWeKfTDozBVI6WThZAydhTMeSmQrFNFCVIl+ybT0YXuFJwtJJGL8J+//rc01Q/4COZegRWz+IRhEa21xxN4dB2430jAQto3++ae4ysrQ9pajHTuKbcQhk5ORfh9aWytafBwyOxuZ5MVGIGwbEbaQnnis6VOwLl+MuWixs6DqNLEXPd1OtU133EAIR6FQjB2iSvThwrzCSMf4LyacNOLZPf1i39MDngSIj3eE2bbR9+zB+PgjxJ5yjIYGpKEhU9LA70O0taIZLmRWJjIlDakbCMtEC4exExKxC6diLlyMtXgxMm/yqWJvWU4IJxQc831mFYpYJ+pE/0KIdIz/QuPSI5rdEwqhNR+jqradii5JcbaHwhQXUoJWfRDjo/fRdu3AqKkBBDIzHdnrR2trdV57U5ETs7B1HS0UBjOInZiAXTwTu9+zz8071XMPBJwQjhBOCCcu88JsVygUo0ZMiP5oTBZe6C+Qi87uCYXQWpqhu4cqv8bTZe2ELYlrr+CByzzM3PoB2uefotccRlgWdkamI9bNzWBJpDcJO28K0qWh9fmc1bSJXsyCAuTlizEXXYGcOPEUsRfdXYjeXmRcPPaEiVHXZ1ahiGViQvRhbE4WXlR2z3HPnu5uSEpCJidTUe8IvgRMy+bgax8wb/NLyHAYmZkF4TCirQXCphOvn5YPukD09SGCIJOSMAuLsRdfgbVg4alib1lOk5JwCJmcrEI4CkWUEjOiPxa5oF8g4TBaQz10dUFKCjIlFaEJ0DRKkixc2JgSDMukdN8n2N5kp4VgW5tTWiEpCWvmbIQUCF8vArCTvJiz5mAvWow1dz5ywoQTYn9yCCct3SnNoFAoopaIiv6PfvQj3n33XYQQ3HLLLfzVX/1VJA93XoyVImjD/QVSVdtG5f4GrshxkZmZDOkZTtl5XYfubvQNHzLngw94pC3I3vSpzO47SqHZiehod3LqvUmYc+YgTAvR04smJHZyKuE5pVgLF2GXznM8e0070We2r0+FcBSKcUbERH/r1q1s2bKF3/3ud5imyU033cSKFSvIz8+P1CGHzUikSY7aQ8M0ObS7mmfebyRsw+/2aDxwfRwz8pKhtxdj48e43nsX7cA+6OmlODWVQl8ToqMDpI1ITMKclo/QdERXNxrSaTY+dwHWZQuxZ805Ifam6fwiMMPY3mRn4lahUIwrIib6ixcv5he/+AWGYXDs2DEsyyJhsKJbl4CTUzhN06Zsd9N5CfeoVM40TbSmI1Qfaua1GpOQfXyzZVNxpIfiys8xfv8G+v490NODSEzGTk6Gri7nF0BCAtakKYi4OERvL5ptY2dmELrscqz5l2EXlSCzsx2x9/ud2vlCQ6alqRCOQjGOiWh4x+Vy8dxzz/Ff//Vf3HDDDUycOHHY3x2q12NWlvei7VoyN4/XNx7GPD7puWnPUVZfNYOSacNrv7ehvAnrpLz/hjYfS+dPumi7ADBNaGiA9nYqrESe+qx3II9fAIaQXPHBy6R++oGTix8fD8mO14+mgTcJcnKcbd3d0NcNebmwbBnMnQtz5kB2ttOQpLMT+nog0QOTCqImhDMS10C0E+vnQI3/wscf8Ync++67jzvvvJO7776bX//619x6663D+l4kG6NnJLpYXprDhp1HAEe4t5Q3kpHoGtb3J2UkoOsaHM+6mZSRcPF2WRaioQG9swNr4kSIT+XT7U2ELRsJCCRz/Ee5Y8drTK3ehd/tAt0FPX1OTRyPB5mWgZ2SgtbXi2huxc7NJbx4CbJkJlbhcc/ethH7DyMsEzs5BRJTwAba+i7O/lEi1ptigzoHavwRbox+oVRXVxMKhZg5cyYej4dVq1Zx4MCBSB3uvFlWmsOmPUcvaMHWiOb9Wxaivh69qxNr4kSsSZMQwSB2vIuSiR5cQmLaEsM2ufWPv6CkuwG/S3dKIse5EYkJ2KmpWClpaL4+9K4e7El5hJZdiZxRiFVU4rQk9PvRjh1FajoyPR3pGt4DTqFQjC8iJvoNDQ0899xzvPTSSwB8+OGH/Nmf/VmkDnfeXKxwX3Tef7/Yd3Zg5eRgeY+LveFCBAJoO3Yw643XeKSykb2JOcxqPkhRy2HQBVJ3oXnc2MlerNR08Pkwujuxpk0neNUK5KQpWCUzkROzEd1daE1HkB7PiT6z58FYyXJSKBQjQ8REf8WKFZSXl/PlL38ZXddZtWoVq1evjtThLohLsmDLshD1deidnVg5uRzwaxzY2UzR9AyKEiXarl2433gVY/Mm6OygUAiK7HLQQMZ7QBfg8WAlp0E4jN7bjZVfQODqayA7G6tkFjIzE9HZ6RRTS0lFpl5YaYeoafXI2H84jXX7FLFDRGP69957L/fee28kDxE92DZabQ1aVxdWbh5WairVtW08/X49pilx7Wrlu50bmfnHN6GrCyklwpIIQ4ckD8KykXEeyM5E+oOIgB+7oIDA1ddCZibW7DlIbzJaTzd0djohnItsOxjpQnUjxVh/OI11+xSxxbhfkXvJPSzbRqupQevqwJo8FSstDeHzYXs87D/mxzKPl00wLQ4cbqOkowPNNCHOBYnJIG2nQ1WaF4mAQABmFBK8egWkpmPNnA0JHkQgAGb4gkI4QxHpQnUjxVh/OI11+xSxxbgW/UvqYdk22uHDaN2djthnpCP6+rDjvKBpGPv2UvrRu7zhvQJT6Bi2xaz6vQgNZFoaWBYAMikZ3C6ElNgzS2D1jYRwYxUWQ3w8QtrY8R5k+tBNSi6UaOlqNdYfTmPdPkVsMa5F/5J4WP1i39PliH1WptNUJC4ZqesYu8sx/ud/MHZ8xixfH494ytiXmc+stkMU2p2QlASm7aRgeuIQmo45pxRz2VWQmAQL5mCHBBiuEQnhnIuxWKjudMb6w2ms26eILca16I+qh2XbaIcOofV2Y06egpyQdbwJuBfpcqN/vhXjrd9h7NyBDASQHe2IUIiieB+FoTaElIiwiR0fh0yKQzPcmHPnYi1aCklJ2LmTICMdcnKoPNzDgepOisNxp3QHi2VRieTDaSTObTQ8PBWxwbgW/VHxsE4W+ylTkRMnoHV3I91upG6gf/Ix+oYPce/ejeXrg9ZWtGAQ6fHAhGwkoAX9SJcbOyUZERePOW8+ocsXU2V52CfSKMnNJr8wFzweKjrD/NvLO08JWcHgPXoVF4+ahFWMN8a16EMEPSzbRjtUTXV9O/tsLyVTcij0+5GGgbRtjPffQ3y2Bffu3djdncjWNvRw0KlamZPrFEML+BAuN7Y3FZGQSHjefOzLFiAMFweTcnhqR5CwHcQ4XM+3MzIpyPOwu7r1jJAVoCYKI4SahFWMN8a96I84loV2uBqtp5eK+Cz+bZsf0+zDtf0Y3746i5LaTxDl5bj27nLKGre0owd8zuKoiTkgQPT2gsuN9CYjvF7C8xdhz5mDECCnTscqLmHfng7C9qEzxKZ0RuagISs1URgZ1CSsYryhRH+4WBbaoSq03j7MadOxc/Ko3FiNadonOlV9sIXSba+jtbZCZxvC53M8+9xJSNtG6+1xJmC9XkjNcDz7khKEEJCfT7h0nlNADSieCkbZmWJTMi190JCVmiiMDGN9EjbW53IU509Uin4kLvQh92maaIerEX0+DsRlcaBbMLOmhYJUgxKvPFEbR1rM3fIuoqEa0d0NcfHI3ElIy0Lr6kAaLvB6qZg+jz1FiynOiqM4QSJnzsacN/+UHrT9tqy7rpA+f/gMmwYLWamJwsgxVs+tmm9QXAhRJ/oXc6EPJeyD7nNiItqhKoTPjzV1GlVaiGde3YdpSd7QYX1OFzPr9/DdyiNsjHPEnfY2tHDIaT5iS+jqRGjC6Sk7cSIV81byhHs+JgKjS/BPNy88w3Z1IyuGi5pvUFwIUSf6F3qhn01MT99n5fZKiqcYWFOng8uFXldDZWXvKU3HKz/dy5ydb0JSDh9f+SeYmsEfpy3mkc3/m6KGvc6q2IRErLxJ2LPnYKdnsc81mXC35uxDMqjt6kZWDBc136C4EKJO9IunpKFpAtuSCE0M+0I/m5gO3DymjSGgqDAHKzsJo7YG6XZheZOZ2X0AF15MKTBsk5LDOxHNx6jIWYSpGdiajmnDXm8eRXFV2JOnYc0sQaZnIidOwFq0hALPBIyXd571JlU3smK4jPX5BsXYJOpEH5wOUvL4v8NlKDGtqmmlcvdh/rwknu6EFEomxFHsP4bd4sPMzsb18UcYVYeYtXsnj/TY7JlQxJyqbRS1VCHT0pjd3YhhW5iAYVuUuAKEv7AKmZGFzMvDunwR9owCEIICzj3hqm5kxfkwVucbFGOXqBP9A3UdWMc7atm2HHb44wwxzfJweEs5z/yxlbAElyZ48HKdfG8CZnwe7j98iDhUzcHqY2wycpEZi1l59CPWlv03JHkhJQVMC0IBVtR+hu1N4qr4XmbMn4E9eQrWvAXYhUVnFD8bzk2qbmSFQhEpok70Lyb8UZCXQkFmPHrtYWgLsL9TELadXw2mLdkXjKN40ydoByugYj+HegSPLb0LU3Myaz7Omcc///H/o6j9MAQCVE4o4vsr7sHUdAwhuXJSD+FFs5BFxSNW6XI8E8vphrE8dsWlJepEH+DKOdmA0/Jw2DdMIIBeXwu9fU6efEoqM/09uHThpFxiU7r1HbTKrRh1jdiJSezNKMUS2oCAm5rOvvTpFLUcomLWUl4pWUVYN5AITGBv/gKmFU8fkTGOd1GI5SylWB674tITVaJfUdN+ys2yrDTn3F/y+9HraxEdHcjkFOwJExC2jXTHkZ9p8Z20I+xv6mPW4Z3MOrQdOzEJrDB69UFmBzX0ki9iHhd9w7YoEd28e9Nf8/OkUiwEIBAC51fH1PQRGWcsiEK0ZimNxMM4WseuGB9EleifXHfGtGxe33iINcvzB79hfD70uhpESysyKwtrRgEiHHa8dl8frnffQT9UTfHhKmYfO+Z4/+ho1dVgGMjEBAqDrfzLx8/x8ZTFyMRElif1YV92OT8zi7HRBg41a1ra0HZcALEgCtGYpTRSD+NoHLti/BBVot9fd6ZfEPce7qCyfsdApckDdR0UZ8ZR5Dt2vOF4HtZlC50WgqEQdHbg+uQT9JrDUF+D3nwM6UkAXUccrATbBnccMikRLRiEHh+FyX7yJ/QiM9zYE6fyX0Yxdq1/wCZNMKKCD7EhCtGYpTRSD+NoHLti/BBVot9fd+b1jYfYe9ipLmlZNmW7myjb04RpSgwNHvhCLvlLr0TraEe0tSI6O3D98WNEfS00NaI3t2J73GAYiOqDThaOKw6ZnICwLURPD/bkSVgLr0B6E7En5GAtvBx71hzk+wfhJNGfV5A54jdtrIhCtGUpjeTDONrGrhg/RJXog3OzLCyewP7aTqQtMTSBdvQopimxcVa6VrSbzGhsQLS14tr4R8SRRmg+htHaih3nApfhhHECAYiLozJ/Pnsy8plTv4ciV5DwjTdz0J3J3uQ8ikunk798ARzvULWsNIeNu5swLYmhC25cMjVi41SiMLaIlYexYnwTdaJf1djFSx8cRNoSTcDXCl3kFk/hk7cOgWVjaIKZVgfu37yJONqIaOtCa2tFagJpuBFV1QhfL9IdD2kZVGRN57ElTlrmq7Nu4CF7NzItlR/0TiPsA+PzPr5d2DdwgxfkpfDA7ZdRtrvpEp8JxaVAPYwV0U50iX4gQOX2g1jHyxlbEjZ0urndNHlwWRoVh1uZ07ibog9qEL3diNZmkCANDe3QIURvD9LlRqalY8d70Lu72D+xGFM3sIWGiWRP8RXI7GzCm2rPGrvdtOcopmWzac/RiGXXjPe0TYVCMfpEl+jX1VGS6UbTnDlXgMNHe3nqrR4eZjdr22oRvT2ItlawLBCg1dYgujqRuo70JmMlezF6e9GEILT6T0iIn4wQTkkH3dApWjoHAGNL/ZCx29HIromFtE2FQjH6RFT0n3/+ed5++20AVqxYwQMPPHBxO5w0iXxvFnMPbWd7fe/AZlPC/i4oqa8H00RKC72uHtHZdrzapRc7NQU9EEAHQtetgqQkDsZP4JeyBNsGTROsu65wyKYkJ3vdo5FdEwtpmwqFYvSJmOiXlZWxceNGXnvtNYQQfPOb3+T999/n+uuvv/Cddnai79yHaG0FGe8IupSAxNvdih0OYDQ2Om0KLQuZlIyVnoEe8KPZFqElyyA9C82tE1pxLbvNiYQ31yEBKSV9/vDAoU6O3Q7mdUd6Qi8W0jYVCsXoEzHRz8rKYv369bjdbgBmzJjBkSNHLm6nn36KsWETqSIftNyT3hD8n7yrqPXBCtFLUbwfa8JEx7MPhzDnX4adm4suILRsBdaSJeDxUNzYhbG14ZzCOpjXvXrptIh63ipTRKFQRAIhpZSRPkhNTQ3r1q3jpZdeYtq0aRe+o5/9DDZtosJn8J3sVZj9q2KPe/xCSlzS4vHtP6ekvYaKeVexO282pbKDkhUL4JprIDHxlF1W1LSzu7qV0hmZlEzHnDc0AAAPJElEQVQbvIxCRU073/3PTZimjWFoPH73lUN+VqFQKMYyEZ/IPXjwIH/zN3/DAw88cF6C39bWi22f+jzKqqqir+kYftMDE2xH7AdCPCA1DdOWbJsyj94FC/mBPg/T0tCNqXx70nwKfDb4ek7ZZ0aii5VznRo+LS0n3js9c+afbjvhdWckuk757GiSleW9ZMceC8T6+EGdAzX+wcevaYKMjKRzfj+ior9t2zbuu+8+HnroIVavXn1R+6pq7GJDWwJFrUH2J6Y75Y6FNuDha9JGahqGgJJUwd6JRYTbded5cJ4ToUNlzqgQi0KhiHYiJvpNTU3cc889PPvssyxduvSi9lXV2MXTv9qOmTAL15JiltR+Bhz38IVAIllhNpFlhyjOTWTKLX9DqFdgvLTjgiZCRzJzJtpz7aPdfoVCcSoRE/0XXniBYDDIk08+ObDttttuY926dee9r7LjZQ8QGmFdsG2KU2DtRGhHMDUjnqvW3QzpGQAUpJy7NeFQFE9JQ9cEpiXRjvfhvRDxi/Zc+2i3X6FQnEnERP/hhx/m4YcfHpF9dfWFTnltndyV6vjf3ZddMSD4/VxMSEae9G9DSy8vfXDwvMXv9F8MZbubIu41j6RnrtYKKBTjj6hYkZuS6D7ltUeD4Em63++NjxQH6joGJpGlLflk1xHCprME+HzE7+Rce6EJNu5uwrJlxLzmkfbM1VoBhWL8ERWif6KypQ0COkX8Ke8X5CZzoM4ptTwSQnqy2GmaoPbYiZny83nAnJxr39YV4ONdRyLqNY+0Z67WCigU44+oEP3+ypb/88EBKpt6z3i/qrGLg41dI+ZBDybW4NTnWX4+fXk5EWKqauxi056jEfWaI+GZq6wlhWJ8ERWi38/hFt8Z2/rnciUj60EPJdbD6ss7xP4i7TUrz1yhUJyLqBH9A3UdWJY98HpBYSb5ucn4AibvfVaPLWVEPOiRFNLR8JqVZ65QKM5G1Ih+8ZQ0DEPDNB2Pu79j1b+9tAPLlmdUyYSRy2RRQqpQKMYLUSP6BXkp3LmmlA3b6lhYPAGA1zcewjyeVcNpVTJVjrlCoVCcSdSIflVjFz99fTdh0+ZAXefxGL6TVikEZ4R2VI65QqFQnEnUiP6Bug7CpiPipnWiEJsAZk1LY2HxhIG0TYC2rgCaJpB2ZGL9CoVCEY1EjegnelycXARa0wDpePgLiycMrJjVNIEALFuia4Kr5uWy7DzTLBUKhWK8EjWi3+cPI3BSMwVQkJtCZ1+QhUUTqDvac9KK2RNPBtuWZKTEK8FXKBSK40SN6BdPSUPXhSPqAiobugB4+9M6Ti7FczIXUyxNoVAoxiNRI/rgePmy/4+Ttw/S+6t/9WxDSy8vvleJbUsMQ2XxKBSK2Ea71AYMl7LdTaeEboZCCNAEGIbGlGwv//e9SixbIgHzeBbP6VQ1dvHW5hqqGrtG3nCFQqEYQ0SNp396eeXB0HXBqssnU9fcw8LiCfT5w6e0XNTEmcXSVD6/QqGIJaJG9E8vr3wyAlgxP5cp2d6BLJ7K+i7WXVeI6/gqXqEJvraq6AxBV/n8CoUiloga0V9WmsMn5Uc4qfzOAPMLM7njhhLe2lxzioD3+cPnrJszWGVKNfGrUCjGK1Ej+gBC9CdtnkrpDKdj1mACfra6Of3ivu66Qvr84YHQjwr3KBSK8UrUiP7J3axOp7/mTn9FzLLdTQPvDeW1DxXLP/3Xggr3KBSK8UTUiH5/lc3+Ugz96PqZk7Ob9hzFtGw+2d00sDr3dK99qFi+ahGoUCjGM1Ej+gV5KTx+95VsKW8k0eOi7qjTwvD0Egsni7ltyYFg0Ole+1DirhqRKBSK8UzUiD5AybR0MhJdZ/3M6c3IBU45htO99rOJu6qfr1AoxitRJfrD4XQxB4b02pW4KxSKWCPiot/b28ttt93Gf/7nfzJp0qRIHw44U8yVsCsUCoVDRMsw7Nq1i3Xr1lFTUxPJwygUCoVimERU9H/961/zz//8z0yYMCGSh1EoFArFMIloeOfxxx+P5O5HBbU6V6FQjCfG7ERuRkbSoNuzsryjZkNFTTvPvLwD07QxDI3H776Skmnpo3b8oRjNczAWifXxgzoHavwXPv4xK/ptbb1nrMDNyvLS0tIzajZsKW880ZfXtNlS3njOlNFIM9rnYKwR6+MHdQ7U+Acfv6aJIZ3lUz4XCaPGC/05/5pArc5VKBTjgjHr6Y8F1OpchUIx3hgV0f/oo4/O+zuaNnjj26G2R4qiyakUTU4d1WOei9E+B2ONWB8/qHOgxn/m+Id7ToSUg3WYVSgUCsV4RMX0FQqFIoZQoq9QKBQxhBJ9hUKhiCGU6CsUCkUMoURfoVAoYggl+gqFQhFDKNFXKBSKGEKJvkKhUMQQSvQVCoUihoga0X/jjTe46aabWLVqFS+++OKlNmdUeP7551m9ejWrV6/m6aefBqCsrIwvfelLrFq1imefffYSWzg6PPXUU6xfvx6A/fv3s3btWr74xS/y3e9+F9M0L7F1keWjjz5i7dq13HjjjTz22GNAbF0Dr7/++sA98NRTTwGxcQ309vZy880309DQAAz9f35B50JGAUePHpXXXHON7OjokH19ffJLX/qSPHjw4KU2K6Js2rRJ3nrrrTIYDMpQKCTvuOMO+cYbb8gVK1bIuro6GQ6H5Te+8Q25YcOGS21qRCkrK5NXXHGFfPDBB6WUUq5evVru2LFDSinld77zHfniiy9eSvMiSl1dnVy+fLlsamqSoVBIrlu3Tm7YsCFmrgGfzycXLVok29raZDgclrfccovctGnTuL8Gdu7cKW+++WY5e/ZsWV9fL/1+/5D/5xdyLqLC0y8rK2PJkiWkpqaSkJDAF7/4Rd55551LbVZEycrKYv369bjdblwuFzNmzKCmpoapU6cyefJkDMPgS1/60rg+D52dnTz77LPcfffdADQ2NhIIBJg/fz4Aa9euHdfjf//997npppvIzs7G5XLx7LPP4vF4YuYasCwL27bx+/2YpolpmhiGMe6vgdPbzJaXlw/6f36h90NUlFZubm4mKytr4PWECRMoLy+/hBZFnsLCwoG/a2pqePvtt/nzP//zM87DsWPHLoV5o8IjjzzC/fffT1NTE3DmdZCVlTWux19bW4vL5eLuu++mqamJlStXUlhYGDPXQFJSEn//93/PjTfeiMfjYdGiRbhcrnF/DZzeZnYw/Tt27NgF3w9R4enbto0QJ8qGSilPeT2eOXjwIN/4xjd44IEHmDx5csych1deeYWcnByWLl06sC3WrgPLsti8eTNPPPEE//3f/015eTn19fUxcw4qKir4zW9+wx/+8Ac++eQTNE1j06ZNMTP+foa67i/0fogKTz87O5vPP/984HVLS8vAT5/xzLZt27jvvvt46KGHWL16NVu3bqWlpWXg/fF8Hn7/+9/T0tLCmjVr6OrqwufzIYQ4Zfytra3jdvwAmZmZLF26lPR0py/zddddxzvvvIOu6wOfGc/XwMaNG1m6dCkZGRmAE7544YUXYuoaAEf/BrvvT98+3HMRFZ7+smXL2Lx5M+3t7fj9ft577z2uvvrqS21WRGlqauKee+7hmWeeYfXq1QDMmzePw4cPU1tbi2VZvPnmm+P2PPzsZz/jzTff5PXXX+e+++7j2muv5Qc/+AFxcXFs27YNcDI7xuv4Aa655ho2btxId3c3lmXxySefcMMNN8TMNVBSUkJZWRk+nw8pJR999BGLFy+OqWsAhr7v8/LyLuhcRIWnP3HiRO6//37uuOMOwuEwt9xyC3Pnzr3UZkWUF154gWAwyJNPPjmw7bbbbuPJJ5/k3nvvJRgMsmLFCm644YZLaOXo88wzz/Dwww/T29vL7NmzueOOOy61SRFj3rx5fPOb3+T2228nHA5z5ZVXsm7dOvLz82PiGli+fDn79u1j7dq1uFwuSktLueuuu7j++utj5hoAiIuLG/K+v5D7QXXOUigUihgiKsI7CoVCoRgZlOgrFApFDKFEX6FQKGIIJfoKhUIRQyjRVygUihhCib5iXLN7927uu+8+1q9fzwsvvABAcXEx7e3tfPjhhwOVKzds2MCPfvSjS2mqQjEqKNFXjGtKS0t57rnnBn3vC1/4Ag8//DDgPBy6urpG0zSF4pIQFYuzFIoL5dNPP+XRRx9lzpw5Z7z36quv8u677/J3f/d3vPzyy1iWhdfr5f777+eVV17hpZdewrZtUlNT+d73vseMGTNYv349nZ2d1NfXs3LlSm655Ra+//3v09fXR0tLCyUlJfzwhz8kLi6OXbt28dhjj+H3+3G5XDzwwAO0tLTwq1/9ipdffhmAI0eO8NWvfpWPPvoIt9s92qdHEYMo0VfEPPPmzeO2226jo6OD+++/n61bt/Lb3/6WF198EY/Hw8aNG/nWt77F22+/DUAgEOCtt94CnAYvX/7yl1mzZg3hcJi1a9eyYcMGrr32Wu655x4ee+wxVq5cyZ49e/jOd77Db37zG5588kkOHjxIYWEhr7zyCn/6p3+qBF8xaijRVyhOY8OGDdTW1nLbbbcNbOvu7qazsxOAhQsXDmz/9re/zaZNm/jpT39KTU0Nzc3N+Hw+Kisr0TSNlStXAjBnzhzeeOMNAL7yla/wyiuv8OCDD/Laa6/xy1/+cvQGp4h5lOgrFKdh2zZr1qzh29/+9sDr5uZmUlJSAEhISBj47D/+4z9iWRY33ngjK1eupKmpCSkluq6fUea2srKS/Px8brvtNm655RYWL15MYWEhkydPHr3BKWIeNZGrUAC6rg/0F12+fDlvvfUWzc3NALz00kt8/etfH/R7Gzdu5J577uGmm24CYNeuXViWRX5+PkIINm3aBMDevXv5+te/jm3b5OTkMH/+fJ544gnWrVs3CqNTKE6gPH2FAliyZAn/9E//xKOPPsr3vvc97rzzTr7xjW8ghCApKYnnn39+0AYV999/P/fccw8JCQkkJSWxaNEi6urqcLvd/PjHP+aJJ57g6aefxuVy8eMf/3ggdr927VoeffRRVqxYMdpDVcQ4qsqmQjHK2LbN97//fXJzc7nrrrsutTmKGEOFdxSKUaS3t5crrriCpqamcV8HXjE2UZ6+QqFQxBDK01coFIoYQom+QqFQxBBK9BUKhSKGUKKvUCgUMYQSfYVCoYghlOgrFApFDPH/AEgXw+qk0GT8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate array of x-values for bootstrap lines: x\n",
    "x = np.array([0, 100])\n",
    "\n",
    "# Plot the bootstrap lines\n",
    "for i in range(100):\n",
    "    _ = plt.plot(x, \n",
    "                 bs_slope_reps[i] * x + bs_intercept_reps[i],\n",
    "                 linewidth=0.5, alpha=0.2, color='red')\n",
    "\n",
    "# Plot the data\n",
    "_ = plt.plot(illiteracy, fertility, marker='.', linestyle='none')\n",
    "\n",
    "# Label axes, set the margins, and show the plot\n",
    "_ = plt.xlabel('illiteracy')\n",
    "_ = plt.ylabel('fertility')\n",
    "plt.margins(0.02)\n",
    "plt.savefig('../images/bootstrap-reg.png')"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}