_1_poisson_statsmodels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple example of poisson regression with `statsmodels`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as tm\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.discrete.discrete_model import Poisson as poisson\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "plt.style.use('seaborn')\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example for generating the data is here: \n",
    "    \n",
    "https://www.theanalysisfactor.com/generalized-linear-models-in-r-part-6-poisson-regression-count-variables/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "days = [1, 2, 3, 3, 4, 4, 4, 6, 7, 8, \n",
    "        8, 8, 8, 12, 14, 15, 17, 17, 17, 18, 19, 19, 20, \n",
    "        23, 23, 23, 24, 24, 25, 26, 27, 28, 29, 34, 36, 36, \n",
    "        42, 42, 43, 43, 44, 44, 44, 44, 45, 46, 48, 48, 49, \n",
    "        49, 53, 53, 53, 54, 55, 56, 56, 58, 60, 63, 65, 67, \n",
    "        67, 68, 71, 71, 72, 72, 72, 73, 74, 74, 74, 75, 75, \n",
    "        80, 81, 81, 81, 81, 88, 88, 90, 93, 93, 94, 95, 95, \n",
    "        95, 96, 96, 97, 98, 100, 101, 102, 103, 104, 105, \n",
    "        106, 107, 108, 109, 110, 111, 112, 113, 114, 115]\n",
    "\n",
    "students = [6, 8, 12, 9, 3, 3, 11, 5, 7, 3, 8, \n",
    "            4, 6, 8, 3, 6, 3, 2, 2, 6, 3, 7, 7, 2, 2, 8, \n",
    "            3, 6, 5, 7, 6, 4, 4, 3, 3, 5, 3, 3, 3, 5, 3, \n",
    "            5, 6, 3, 3, 3, 3, 2, 3, 1, 3, 3, 5, 4, 4, 3, \n",
    "            5, 4, 3, 5, 3, 4, 2, 3, 3, 1, 3, 2, 5, 4, 3, \n",
    "            0, 3, 3, 4, 0, 3, 3, 4, 0, 2, 2, 1, 1, 2, 0, \n",
    "            2, 1, 1, 0, 0, 1, 1, 2, 2, 1, 1, 1, 1, 0, 0, \n",
    "            0, 1, 1, 0, 0, 0, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = np.array([days, students]).T\n",
    "df = pd.DataFrame(dat, columns = ['days', 'students'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following ways of calling the `statsmodels` api are equivalent in this example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 1.784900\n",
      "         Iterations: 8\n",
      "         Function evaluations: 12\n",
      "         Gradient evaluations: 12\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 1.784900\n",
      "         Iterations: 8\n",
      "         Function evaluations: 12\n",
      "         Gradient evaluations: 12\n"
     ]
    }
   ],
   "source": [
    "m1 = smf.poisson('students~days', data=df).fit(method='bfgs')\n",
    "\n",
    "m2 = poisson(endog = dat[:, 1], \n",
    "             exog = np.concatenate([np.ones(dat[:, 0].shape).reshape(-1,1),\n",
    "                                    dat[:, 0].reshape(-1,1)], axis=1)).fit(method='bfgs')\n",
    "\n",
    "m3 = sm.GLM(endog = dat[:, 1],\n",
    "            exog = np.concatenate([np.ones(dat[:, 0].shape).reshape(-1,1),\n",
    "                                   dat[:, 0].reshape(-1,1)], axis=1),\n",
    "            family = sm.families.Poisson()).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Poisson Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>students</td>     <th>  No. Observations:  </th>  <td>   109</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                <td>Poisson</td>     <th>  Df Residuals:      </th>  <td>   107</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                 <td>MLE</td>       <th>  Df Model:          </th>  <td>     1</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Sun, 12 Apr 2020</td> <th>  Pseudo R-squ.:     </th>  <td>0.2269</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>12:06:07</td>     <th>  Log-Likelihood:    </th> <td> -194.55</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>             <td>True</td>       <th>  LL-Null:           </th> <td> -251.64</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>  LLR p-value:       </th> <td>1.189e-26</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    1.9902</td> <td>    0.084</td> <td>   23.712</td> <td> 0.000</td> <td>    1.826</td> <td>    2.155</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>days</th>      <td>   -0.0175</td> <td>    0.002</td> <td>  -10.114</td> <td> 0.000</td> <td>   -0.021</td> <td>   -0.014</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                          Poisson Regression Results                          \n",
       "==============================================================================\n",
       "Dep. Variable:               students   No. Observations:                  109\n",
       "Model:                        Poisson   Df Residuals:                      107\n",
       "Method:                           MLE   Df Model:                            1\n",
       "Date:                Sun, 12 Apr 2020   Pseudo R-squ.:                  0.2269\n",
       "Time:                        12:06:07   Log-Likelihood:                -194.55\n",
       "converged:                       True   LL-Null:                       -251.64\n",
       "Covariance Type:            nonrobust   LLR p-value:                 1.189e-26\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      1.9902      0.084     23.712      0.000       1.826       2.155\n",
       "days          -0.0175      0.002    -10.114      0.000      -0.021      -0.014\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m1.summary()\n",
    "# m2.summary()\n",
    "# m3.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = m1.predict(df.days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFyCAYAAABbdsanAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxcdfX/8dfs2dMs031vmltKS0kpbaGUYikgm2xCkUXZlK+KIsoqKEXBnygIuKCoKIKiVfZNllIopSxdgbLdUKDQlkLTJM2e2X9/TCxdJslkMnfuJHk/H4880rnL53PuuXduT+7qiMViiIiIiIh9nHYHICIiIjLQqSATERERsZkKMhERERGbqSATERERsZkKMhERERGbqSATERERsZnb7gCSUVPTlPZnc5SU5FFf35ruZmUXyrH1lGPrKcfWU46tpxxbb9cc+/2Fjp7OP2CPkLndLrtD6PeUY+spx9ZTjq2nHFtPObZeb3M8YAsyERERkWyhgkxERETEZirIRERERGymgkxERETEZirIRERERGymgkxERETEZirIRERERGymgkxERER67Omnn+SWW37BjTfewEMP3c+11161c9y//33vbtPu+bk7r7/+GueddyZr167uVYwvvfQip576JbZu/aRX7WSCCjIRERHpsSeeeITTTjuDyy67iuOPP5FFi362c9y///3P3abd83N3pk3bnwkTJvY6xoMPPoQhQ4b2up1M6BOvTsoGgVCEhuYAxQU+fB498VhERLJD/qJr8D36UNcTOR2URpN/C2Hg+BNpWXR9p+NffHEZmzdv4r77FjNz5iw2b97M4sX/4L77HuXZZ5+hubmJO++8gzFjxuJwOHf7vGDBUfz5z38gEongdDrJy8vjzDO/BsCtt/6ScDjM8OEj2LZt2179rlr1Ktdddw0nn3wq5533De6668+88cbrXHfdz/jZz65jn30ms23bNvbbbxpHHnn0bvO++eZ6brrp//Hd736fcePG84tf3EBFRSXnn38hLS3N3HbbzYwaNZpt27ZxyCGHMmvWQSxbtpTXXltHeXk57777Dj/96c+TzmFPWVaQGYYxFLgemGaa5oEdw64AhgJbgRnAj03TfNeqGNIhEo2yeOkG1lXXUNcYoLTIR1Wln4XzK3A5dYBRREQGnkMOmce//vUPTjvtKwwbNhyAxYv/AcDhhx/B73//a84//8Kd0+/6+dVXX+att9Zzyy2/A+Cii77BzJmzqampYdOmTdx8868BWL582V79HnjgLL70pZPIyckFwO12c+mlV1JYWMgxxxzH3LmHEYlEOPPMU/cqyKZMmcrEiZUAlJSUMnfuYTtPZd5zz12MHDmKs88+l0CgnTPO+DKLFz/Ek08+zqGHfoGjjz6O9etfT1v+ErHyCNkhwMPA/rsMKwC+b5pmzDCMhcAvgeMtjKHXFi/dwJLVm3d+rm0M7Px8xoJKu8ISEREBoGXR9V0ezQLw+wupq2nKUERde//992hvb+eee+4CYMiQIezYUc+HH77PqFGjdk43fPiIhPOfeOIpXHzxNzn11NPZtu0zhg8fQTgcZuPGD6muNvF6fezYUd/jmIqKinfGNGFCBU1NjVx00SX8/e93cd99iznooDlMmbIfDkeP3xueFMsKMtM07zMM47A9hv1ol49OoNmq/tMhEIqwrrom4bh11ds5Zd4Enb4UERHZg9PpJBaLsWHDe0ycWLnb54qKSt58cz1nn30OAGvWrGLEiFEEgyHWrPn8Iv5PPtmSsO3Bg4cwYcJEfvGLGzjqqGMAePnlF1m16lV+/es/AHD//YsTzpuXl0drawsAn3326c7hFRWVlJaWceqppwPw1FNPUFRUzFtvvckVV1xDOBzmoou+wdy5h2EYk3qXnE7Ycg2ZYRhe4GvAt5OZvqQkz5I31fv9hV2O37q9hbqmQMJx9U3tuLwe/OX5aY+rP+kux9J7yrH1lGPrKcfWS2eOn3/+eWpqPuO//32IM844gxdeeIHW1hZefXUZxx13HIcfPp+//OV2AA4++IDdPl955ZV89NF73H33H8nPz6ehoYEjjpjHlCkVvP76Km699ecMGzYMt9vJsmXPMGfOgRQUFOzW/wUXnMuiRYu49dabAZg372AefPDf/OEPtzJkyBACgXaWLXuK8vJyamo+4+mnH+XSSy/l9NNP5fbbb6epqY7W1kY++OADGho+43vfu4hf/vKX/Oc/9xAKhfD7/QwdOoiNG6vZuLGanJwcJk+exMyZ0/B6vZbk2BGLJX+RX091HCG7yTTNGbsM8wK/B243TXNNMu3U1DSlPUi/v5Cabg7fBkIRrvnTK9Q27l2UlRXlcP3XZ+kIWReSybH0jnJsPeXYesqx9ZRj6+2aY7+/sMfnNTN6VbphGHnAHcCvTNNcYxjGKZnsv6d8HhdVlf6E46oqy1WMiYiISFpYeZflPOBsYJhhGNcANwN/B6YA4wzDAMgH7rcqhnRYOL8CiF8zVt/UTklhDlWV5TuHi4iIiPSWlRf1LwP2vGf1ZKv6s4rL6eSMBZWcMm+CnkMmIiIiltCDYZPk87gYXJJndxgiIiLSD+nJpiIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZCIiIiI2U0EmIiIiYjMVZL0QCEXYVt9KIBSxOxQRERHpw9x2B9AXRaJRFi/dwLrqGuoaA5QW+aiq9LNwfgUup2pcERER6RkVZClYvHQDS1Zv3vm5tjGw8/MZCyrtCktERET6KB3O6aFAKMK66pqE49ZVb9fpSxEREekxFWQ91NAcoK4xkHBcfVM7Dc2Jx4mIiIh0RgVZDxUX+Cgt8iUcV1KYQ3FB4nEiIiIinVFB1kM+j4uqSn/CcVWV5fg8rgxHJCIiIn2dLupPwcL5FUD8mrH6pnZKCnOoqizfOVxERESkJ1SQpcDldHLGgkpOmTeBhuYAxQU+HRkTERGRlKkg6wWfx8Xgkjy7wxAREZE+TteQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzdxWNWwYxlDgemCaaZoHdgzLAW4CtgATgZ+bplltVQwiIiIifYGVR8gOAR4GHLsM+x7wsWma/w+4BbjTwv5FRERE+gTLCjLTNO8DmvYYfCzwcsf49cA0wzCKrIpBREREpC+w7JRlJwaze5HW2DGssauZSkrycLtdaQ/G7y9Me5uyO+XYesqx9ZRj6ynH1lOOrdebHGe6INsG7BptUcewLtXXt6Y9EL+/kJqaPQ/gSTopx9ZTjq2nHFtPObaecmy9XXOcSmGW6bssHwcOAjAMYyrwummaXR4dExEREenvrLzLch5wNjDMMIxrgJuB24CbOj5XAOdb1b+IiIhIX2FZQWaa5jJgWYJR37aqTxEREZG+SA+GFREREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgS1IgFGFbfSuBUMTuUERERKSfcdsdQLaLRKMsXrqBddU11DUGKC3yUVXpZ+H8ClxO1bMiIiLSeyrIurF46QaWrN6883NtY2Dn5zMWVNoVloiIiPQjOsTThUAowrrqmoTj1lVv1+lLERERSQsVZF1oaA5Q1xhIOK6+qZ2G5sTjRERERHpCBVkXigt8lBb5Eo4rKcyhuCDxOBEREZGeUEHWBZ/HRVWlP+G4qspyfB5XhiMSERGR/kgX9Xdj4fwKIH7NWH1TOyWFOVRVlu8cLiIiItJbKsi64XI6OWNBJafMm0BDc4DiAp+OjImIiEhaqSBLks/jYnBJnt1hiIiISD+ka8hEREREbKaCTERERMRmKshEREREbKaCTERERMRmKshEREREbKaCTERERMRmKshEREREbKaCTERERMRmKshEREREbKaCTERERMRmKshEREREbKaCTERERMRmKshEREREbKaCTERERMRmKsg6EQhF2FbfSiAU6RN9ZTJeq/SHZRAREUmF245ODcO4DBgLbAcmAuebptlmRyx7ikSjLF66gXXVNdQ1Bigt8lFV6Wfh/ApczvTWr+noK5PxWqU/LIOIiEhvZLwgMwxjKHAVUG6aZtQwjIeBk4F/ZDqWRBYv3cCS1Zt3fq5tDOz8fMaCyqzrK5PxWqU/LIOIiEhv2HH4oRUIAkUdnwuAt2yIYy+BUIR11TUJx62r3p7WU2np6CuT8VqlPyyDiIhIb2X8CJlpmo0dpywXG4axFdgMbOhqnpKSPNxuV9pj8fsLd/u8dXsLdU2BhNPWN7Xj8nrwl+enpe909JXJeFO1Z4731BeWIdt1l2PpPeXYesqx9ZRj6/Umx3acstwfuAyYbppm2DCMm4EfA5d3Nk99fWva4/D7C6mpadptWCQUobTQR23j3gVCSWEOkWBor3lSlY6+MhlvKhLleE/ZvgzZLpkcS+8ox9ZTjq2nHFtv1xynUpjZccpyBFBnmma44/NWIMeGOPbi87ioqvQnHFdVWY7Pk76jdOnoK5PxWqU/LIOIiEhv2XGX5ZPAMR1HxnYAU4Dv2RBHQgvnVwDx65fqm9opKcyhqrJ85/Bs6yuT8VqlPyyDiIhIbzhisZjdMXSrpqYp7UF2d/g2EIrQ0ByguMBn+VGadPSVyXiT1dND5Nm4DNlOpyGspxxbTzm2nnJsvT1OWTp6Or8tzyHrC3weF4NL8vpMX5mM1yr9YRlERERSoaduioiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBZmIiIiIzVSQiYiIiNhMBRngaNhBzp//gKOmxu5QREREZABSQQZ4Vr5C4Q8vp3TeLLz/fdzucERERGSAUUEGBA8/kuaf/AxHUxPFX/sKBRd/C0dToyV9BUIRttW3EghFLGnfLnsuVyAUYev2lowsZ3/NqYiIDBxuuwPICk4nbf93EcHDDqfw298g959/x/viCzT95g+EDj4kLV1EolEWL93Auuoa6hoDlBb5qKr0s3B+BS5n362L91yukkIv+bleWttD1DUFKC20bjn7a05FRGTg0f9au4hM2ocd/32Wlu9fhnPLZgadeAz5P7oS2tp63fbipRtYsnoztY0BYkBtY4AlqzezeOmG3gduoz2Xq64pyKZtzfHPMWuXs7/mVEREBh4VZHvyemm98kfseOxpwhMqyLvjdkrmz8G9emXKTQZCEdZVJ75hYF319j57qq2r5dpTupezv+ZUREQGJhVknQjPmEn9sy/SeuG3cH3wPoOOO5L86xdBINDjthqaA9Q1Jp6vvqmdhuaet5kNulquPaV7OftrTkVEZGBSQdaVvDxafvpzGh58nOjI0eT9+leUHHEo7tfW9qiZ4gIfpUW+hONKCnMoLkg8Ltt1tVx7Svdy9teciojIwKSCLAmhgw+h/vkVtJ1zPu5332HQ0YeTf8N1SR8t83lcVFX6E46rqizH53GlM9yM6Wq59pTu5eyvORURkYGpxwWZYRhOwzCKrAgmm8UKCmn+xS3suP9RoiNGknfbzZQsmIt73Zqk5l84v4IFM0ZSVpSD0wFlRTksmDGShfMrLI7cWnsuV2mhj1GDCygr8lm+nP01pyIiMvA4YrFYtxMZhvFLYDtwK7AKGAb83DTNm60NL66mpqn7IHvI7y+kpqYptZmbmyn46Y/J/eufiXU8MqPl8h9CXl63swZCERqaAxQX+PrVUZw9lysQiuDyeogEQ5YvZ3/NaTJ6tR1LUpRj6ynH1lOOrbdrjv3+QkdP50/2CFmeaZo3AicDH5im6Qf26Wln/UZBAc03/oodDzxGdPQY8m7/NaWHHYRnxfJuZ/V5XAwuyet3hcOey+XzuBhWnp+R5eyvORURkYEj2YKspeP3acAfO/5dn/5w+pbQIYdS9/zLtH7ruzg//ohBJx1LwQ8uxtHYYHdoIiIi0ockW5ANMwzjMWAa8JRhGPOBfa0Lqw/Jy6Nl0fXseGIJ4X0mk3vPXymZOwvvU/+1OzIRERHpI5ItyL4B3AnMM00zAuQC11oWVR8Unj6D+mdeoOXyH+LcXkPx2Qsp/L/zcGzfbndoIiIikuWSLchOME3zQdM0N3V89gNHWhRT3+X10nrpldQ/+yKhA2aQ88B9lB4yA999iyGJmydERERkYEq2IDt41w+mad4FjEt7NP1EZNI+7HjsGZp/+v9wtLdT9K2vU7zwJJwbP7Q7NBEREclC7q5GGobxHBADJhqGMWWXUS6gx7d0DiguF20XfpvAF4+l8PJL8D73LKXzZtNy6VW0/d+3weOxO0IRERHJEl0WZMCijt8XA7ftMrwdeMOKgPqb6JixNPzrAXwP3kfBNVdQ8NMfk/PAf2i66VbCBxxod3giIiKSBbosyEzTXAZgGMYbpmnu9pgLwzAOJP6QWOmOw0Hg5FMJfuFw8q/7Ebn33sOgYxbQfs75tFx9LbGiYrsjFBERERt1d4Tsf2KGYXwdGMLn150dA8y2JKp+KlZSSvOtvyOw8AwKLvseuX/9M97HH6Xl+p8TOOFkcOgssIiIyECU7EX9TwCHAV7i147970dSEDpoDvVLV9By1Y9wNjZQ9I1zKT7tRFzvVdsdmoiIiNgg2SNkbaZpnrnrAMMwHrIgnoHD66X1kstoP/EUCq+6FO/SJZTMm03bhd+m9QeXEysotDtCERERyZBkj5A9YBjGoYZh7Hpr4LFWBDTQRMeNp+Gf99Pwt38SHT6CvN/dRsnBM/A98B89u0xERGSASLYg+w3wPBAwDCNiGEYU+KllUQ00DgfBo4+lbvlKWi69Emd9HUX/dz7FJx6D6+237I5ORERELJZsQfawaZrOjh+XaZpO4AYrAxuQcnNpvfyH1C1fSeCLx+J9eQUlhx9C/jVX4GjYYXd0IiIiYpGkCjLTNE9KMOzH6Q9HAKJjx9F49z9p+Od9REaPIe+Pv6f0oAPw/esfEI3aHZ6IiIikWVIFmWEYFYZhLDMM4wXDMPINw3jMMIyxFseWMYFQhG31rQRCkazqP3j4kdS/8CrNV1+Lo7WFou9+k0HHHYn7jdcs7zvd86Qq2+PrD5QvERH7JXuX5XXAT4AzTdNsMQzjAuKnLM9PpVPDMAzgK0AbMA9YZJrmylTa6o1INMripRtYV11DXWOA0iIfVZV+Fs6vwOVM9myuxf37fLRd/AMCp5xG/qJryHnkQQYdMY/2r55Hyw9/RKyk1Lq+0zBPqrI9vv5A+RIRyR7J7nU3mqb5LBAAME3zU6C+61kSMwzDBfwK+IlpmjcSL+pseev24qUbWLJ6M7WNAWJAbWOAJas3s3jphqzrPzpyFE1//hs77nuEyMRKcv92J6UHTSfn7r9CpOdHNlJZ9kzmK9vj6w+ULxGR7JFsQTbMMIxc4i8axzCM0cDEFPs8kPhDZb9jGMZVwPHA9hTbSll7MMy66pqE49ZVb7f89E0gFEmp/9Chh1G/dAXNi26AQJDCSy9m0Bfn416d/AHGVPpONd5UZHt8/YHyJSKSXZI9ZXk38DaQaxjGPGAwcFqKfY4BDgK+Yppmg2EYfweCwF2dzVBSkofb7Uqxu8S2bm+hrimQcFx9Uzsurwd/eX5a+0xr/9f+EL5+Dlx+OZ5//IOSYxbAeefBDTfA0KFp7zvVeP3+nj/gNpPx9QeZyvFAlkqOpWeUY+spx9brTY6TKshM03zeMIwDiBdSAC+bplmXYp+NwLumaTZ0fH6R+GuZ7upshvr61hS76lxJcS6lhT5qG/f+T6mkMIdIMERNTVPa+/2fSCjS+/49hXDL7/GcdhYFV16K+y9/Ibr437Rd/H1aL/w25Oamre9U5vH7C1PKYabi6w8ymeOBKtUcS/KUY+spx9bbNcepFGZJX7lrmmadaZqPd/zUGYZxZY97i3sVKOu4lgziR8wy/hLHHK+bqkp/wnFVleX4POk9Ircnn8eVtv5DB82h/tnlNN34K8jxkf+zn1A6p/On/afSdzrj7U62x9cfKF8iItmlyyNkhmF8SMd1YwmUAj/vaYcdxdwVwK2GYdQAfuJ3cGbcwvkVQPyamfqmdkoKc6iqLN85vE/173bTfu4FBE45lbxbbiL3T7+n6P/OJ/SnP9D8k58RPnBWr/vOZL6yPb7+QPkSEckejlgX70vsuL7rauCrwCfA8o5Rc4EppmleYnmEQE1NU9pf6rjrocVAKEJDc4DiAp8tRwas6N+58UMKfnotvkfj74BvP/FkWn70E6KjRve672TnScchcivj6w/syvFAolM91lOOraccW2+PU5aOns7f3SnLC0zT/AgYZZrmnaZpVnf83Al4upm3z/B5XAwuybPtPyMr+o+OHUfjnXdT/8hThPavIuehByg9+ADyr1+Eo6mxV31nMl/ZHl9/oHyJiNivy4LMNM32jn/uYxjG4P8NNwxjCDDNysAkPcKzD2LHk8/R+Ls/Ei0rJ+/Xv6J05jRy7rwDQiG7wxMRERGSv6j/NuBdwzDeMAzjDeAd4BbrwpK0cjoJnHo6dS+toeWqH8WfX3bVZZTMnYn3sUcSXvgvIiIimZPsy8XvAyYRv57sasAwTfMBKwMTC+Tl0XrJZdS9+hpt530d18cfUXzeWQw69gg8L6+wOzoREZEBqyePvdhmmuajHT81hmGk9B5LsV/M76f55zdTv/xVAsd+Cc/qlQw64WiKzjwV19tv2R2eiIjIgJPUg2ENw1iaYPBE4M70hiOZFJkwkca//h336pXkX78I3zNP4V3yNIEvL6TliquJjh5jd4giIiIDQrKvTtpB/DoyiN9duT/x55BJPxCeMZOGBx/Hu/QZ8q+/jpz//Avfww/Qdu4FtH7vMmJlZXaHKCIi0q8le8ryLNM0l3X8LDFN8yb60WMvBHA4CB5+JPXPLqfx9j8RHTqMvDtup/TA/ci7+UZobrY7QhERkX4r2SNkMwzD+N+/ncAw4GBLIhJ7OZ0EvryQwPEnknv3X8i75Zfk33gDuXfeQevFP6Dta+dDTo7dUYqIiPQryR4huxdYBFwH/Bg4HviBRTFJNvD5aPv6N6lb+Totl10F7QEKfnQVpbOryLnnLj3DTEREJI2SPUJ2lWma91gaiWSlWEEhrZddRdv53yDvN7eSe+cdFP7gu+T95hZaLruKwMmngktPeBcREemNZI+QFe76wTCMcwzDuK2ziaX/iZWW0XLtT6lb+Tpt516Ac8tmir79DUrmzcb38AMQjdodooiISJ+VbEE2adcPpmneBeSlPRrJetGhw2i+8VfUvbyWtjO/iuv9DRR9/RxK5h+C9/FH9dR/ERGRFHR5ytIwjA+BGFBqGMZxu4xyAW9aGZhkt+joMTTf8ltav3MJ+Tf9HN/9/6b43DMJTdmP1suuIvjFY+wOUUREpM/o7gjZYcB84KGOf38BOAI40DTNYy2NTPqE6PgJNN3+J+qXr6T9pFNwv7We4q99hUFHzINHdcRMREQkGV0WZKZpfmSa5kagBvgKsBV4EHjLMIzvWx9e3xcIRdhW30ogFLE7FEtFKg2a7vgr9cteof2Ek3Gvfx2+9KX4qcxHH0rbNWb9IZ92LkOivrM9p6nEl205zibZHp/IQJXsXZa5pmneaBjGV4APTNPczzCMP1sZWF8XiUZZvHQD66prqGsMUFrko6rSz8L5FbicSb9CtM+JTNqHpj/dRev3L6f0D7fhWryY4vO/StiYROv3LiVwwsngTnaz26XdfpBPO5chUd/TJpbjAF57b3tW5jSVfGVbjvt6PkUkc5L9FrZ2/D4N+GPHv+vTH07/sXjpBpas3kxtY4AYUNsYYMnqzSxeusHu0DIiss9kuPde6lesov30M3FteI+ib15AyZwZ+P759x4/x6w/5NPOZUjU99I1W3h2zZaszWkq+cq2HPf1fIpI5iRbkA0zDOMxYBrwlGEY84HJ1oXVtwVCEdZV1yQct656+4A6VRCZMJGmX/+eulfW0Xb2ubg2b6Lo4m/FHzB7150QCHTbRn/Ip53L0FXfdsSTjFTyla057qv5FJHMSrYg+wZwJzDPNM0IkEv8qf2SQENzgLrGxIVGfVM7Dc3dFyH9TXTMWJpvvo26VW/QesGFOGu2UXj5JZQeuB+5f7wdWls7nbc/5NPOZeiqbzviSUYq+crWHPfVfIpIZiVVkJmm2Waa5oOmaW7q+Py4aZorrQ2t7you8FFa5Es4rqQwh+KCxOMGgujwEbT87JfUrlpP67e+i7OxkYJrrqRsxlRyf3MrjuamvebpD/m0cxm66tuOeJKRSr6yNcd9NZ8iklm6ktMCPo+Lqkp/wnFVleX4PHrVUGzIEFoWXU/tmjdpueRSCAQo+OmPKT1gCnm/+BmOutqd0/aHfNq5DF31bUc8yUglX9ma476aTxHJLNeiRYvsjqFbra3BReluMz/fR2trMN3N7jR5bAltgTANzUECwTClRTnMmTqUhfMrcDoclvWbTZLKcV4eobnzaD/nPGK5uXjWrsa3dAm5f/kTjtrtRIxJxIqK+kU+rViGZLfjRH0fNGUI44cX0dgSysqcppKvbMtxX89npli9PxblOBN2zXF+vq/Hl3U5Yn3gwZ01NU1pD9LvL6SmZu/TY+kWCEVoaA5QXOAbcH+FppTjlhZy//E3cm//Da5PthBzuwl8eSGt376YiDGpX+QzncvQ0xwn6jvbc5pKfNmW42ySjfFlan88kCnH1ts1x35/YY//ylFBJpbpVY6DQXwP/Ie839yC+71qAAJHfpHWiy4hPGs2ZMERh2yg7dh6yrH1lGPrKcfW621BpmvIJDt5vQROP5P65StpuPtfhA6chat7IgQAACAASURBVO/pJyn50lEMOvYIvI89AhHdqi8iIv2DCjLJbk4nwS8ew47Hn6H+0acJHHU0ntUrKT7vLEoOPoCcv/65y0dmiIiI9AUqyKTPCM+aTeM9i6lbsZq2s8/B9ckWCq/4PmXTJ5P38+txfPaZ3SGKiIikRAWZ9DmRiZU03/xrate8Rcv3L4NYjPxf/YKyA/al4OJv4XrnbbtDFBER6REVZNJnxQYPpvXKH1G79m2abvwVkZGjyP3n3ymdN5viU0/A++zTEI3aHaaIiEi3VJBJ35efT/u5F1D/0hoa7llM8OBD8C57juKvfJmSuTPj78xsabE7ShERkU6pIJP+w+kkeNTRNDz0BPXPLqd94Rm4Nn5I4eWXUFa1D/k/+THOTR/bHaWIiMheVJBJvxSeOo2m3/yB2rVv0/KDK8DtJu+3t1J64H4UnXsWnpdehD7wDD4RERkYVJBJvxYbMoTWK66mdu3bNP7694T3nYrv8UcYdOIxlBx2MDl//XPCF5qLiIhkkgoyGRhycgicfiY7lrxA/SNP0f6lk3C9Z1J4xfcpnWpQcNkluN5cb3eUIiIyQKkgk4HF4SA8+yCa/vw36ta9TcsVVxMrLib3b3dSOn8Og45ZgO/f/4T2drsjFRGRAUQFmQxY0SFDaf3BFdStXk/D3f8icPgRuNesouiiCymbZpB/7dW4Pthgd5giIjIAqCATcbsJfvEYGv95P3Wvvkbrdy4Bl4u83/+G0tnTKf7yCfF3Z4bDdkcqIiL9lAoykV1Ex46j5UfXUbvuHRr/cCfB2QfjfeE5is87i9Lp+5L3i5/h/GSL3WGKiEg/o4JMJBGfj8DJp9LwyJPUvfAqbed/A0dLC/k3/ZzSA6ZQ9LUz8Dz3rN4EICIiaaGCLAsEQhG21bcSCEXsDqXf2jPHPcl5ZNI+NP+/m6h9w6TpV7+JPzrjv48xaOFJlM6uIve3t+GorbV6EaRDonWXynco2793mYwvmb4yFU+2rxcRqzhiNj0c0zCMXOBV4GnTNC/tatqamqa0B+n3F1JTY+/zpyLRKIuXbmBddQ11jQFKi3xUVfpZOL8Cl7Pv18rZmOOSQi/5uV5a20Op5zwWw/3aWnLuupOch+7H0dZGzOslcPyJtJ1zAeGZs8DhsHbBOmRDjjMl0fdl2sRyHMBr721Pen329HuX6Rxncr+QTF+ZiMfvL+TTzxr69f7QbgNpX2GXXXPs9xf2+D8B16JFi9IdU1J++9vf3gg0Aq3f+c53nu5q2tbW4KJ095+f76O1NZjuZnvkX8++x5LVm2kLxP8SbAtE+OCTRtoCYaaOL7M1tnTIyhwHIzS2BHuXc4eD6LDhBI8+lrZzLyA6ZCiujR/iffEFcv95D77HH8ERChIZOx7y8qxaNCA7cpwpib4vH25t4sOtTT1anz393mU6x5ncLyTTVybiyc/3ceejb/Xr/aHdBtK+wi675jg/33ddT+e35c8OwzDOBlYAH9rRfzYIhCKsq65JOG5d9XYdrk+DrnK8p1RzHhtUQts3vkX9i6vY8eDjtJ9wMq4N71FwzZWU7VdJ0flfxbP0GYhoffZGT9YldL4+s/17l8n4kukrU/G0B8NZvV5EMsGd6Q4Nw5gM7GOa5g8Nw9gvmXlKSvJwu11pj8XvL0x7m8naur2FuqZAwnH1Te24vB785fkZjir9sjXHe0pLzk88Jv5TUwN//zuOO+/E9+hD+B59CEaOhHPOgXPPhfHjU+8jATtznCk9WZfQ+fpM9XuXqRxncr+QTF9ARuIZKPtDuw2EfYXdepPjjBdkwElAu2EYVwKHAF7DML5nmuatnc1QX9+a9iDsPp8eCUUoLfRR27j3TqikMIdIMNTnz/dnc473lN6c58BZF8CZ5+Net4ace/+O78H7cF5/PVx/PcE5c2n/ylkEjjkeCgp61ZPdOc6UnqxL6Hx9pvK9y2SOM7lfSKYvICPxlBTn9vv9od0Gyr7CTntcQ9bj+TN+ytI0zRtM0/yJaZo/B14EVnZVjPVXPo+Lqkp/wnFVleX4POk/IjjQdJXjPVmSc4eD8PQZNN90K7Xrq2n87R0E58zFu2I5RRddSPmUCgq/eQHeZ5/WQ2e70ZN1CZ2vz2z/3mUyvmT6ylQ8OV53Vq8XkUyw7aJ+wzBOAU4Fxvz2t79t+853vvNmZ9P214v6J48toS0QpqE5SCAYprQohzlTh7JwfgXODN2lZ6VszHFJoY/yQbl4XA4CwUjmcu7xENl3KoHTz6T9lNOIDSrB9dFGvC+vIOf+f5P7t7/g/GQLsdJSokOHJX2XZjbkOFMSfV8OmjKE8cOLaGwJJf0d6un3LtM5zuR+IZm+MhFPfr6PsYPz+/X+0G4DaV9hl95e1G/bYy96or8+9uJ/AqEIDc0Bigt8/eovwWzOcVbkPBbDvWYVOff/G99D9+PseJZZePwEAl9eSPsppxEd1/X1ZtmU40xJtO5SWZ/JzmNXjjO5jSbTl5Xx7JrjrPhu9kMDcV+Rab197IUKMrGMctwDoRDe55/Fd99ifP99HEd7e3zwjJm0f3khgRNOJla2963/yrH1lGPrKcfWU46t19uCTE/bE8kGHg/BI75I0x1/pfatDTT++vcE530B99rVFF75A8qmTqTorNPwPXQ/tKb/JhcREbGXHXdZikgXYoVFBE4/k8DpZ+L8dCu+B++PHzl7+kl8Tz9JtKCQ4HFfov2U0+CkY+0OV0RE0kBHyESyWHToMNq+eRE7nl1O3fKVtHzvUmIlJeT86x8MOvUEGDWK/GuvxrX+DegDlx+IiEhiKshE+oiIMYnWH/6YulVvsOORJ2k7+1xobyfv97+h9PBDKDl0Frm33Yxz08d2hyoiIj2kgkykr3E6Cc0+mOabb4OtW2m4614Cx52A68MPKLjhOsoOmELxCUeTc89dOHbU2x2tiIgkQdeQifRlPh/BY44jeMxxOBp24Hv0YXz3/xvviuV4X15BwVWXEpx/BIETTyZw5NG9fjOAiIhYQwWZSD8RKx5E+1lfo/2sr+HcvAnfA/eRc/9ifE8+ju/Jx4nl5BBccBSBE04isOAoyNe7AUVEsoVOWYr0Q9GRo2j77iXUL3uFuhdepeXSK4mMGo3vsYcp+vo5lE8eT9F5Z8cfo9HcbHe4IiIDno6QifRzkUn70DppH1ovuwrXu+/ge/h+fI88hO+xh/E99nD8yNn8Iwgc9yWCR36RWFGx3SGLiAw4KshEBgqHg8g+k2ndZzKtV1wTL84eeRDfow/he+JRfE88SszrJTjvCwSOO4HgUUcTK9377QAiIpJ+KshEBqLdirOrcZnvdhwxewTfM0/he+YpYi4XoYMPIXDM8QSPOY7osOF2Ry0i0m+pIBOR+DPOjEm0/uAKXB9swPv4Y/ieeATv8mV4ly+Dqy4lNP0AAkcfR/Do44hMrARHj1/VJiIinVBBJiK7iYyvoO0736PtO9/D+ckWvP99DN/jj+J5eQWetWvghusIT6ggePRxBI4+lvABB4JT9weJiPSG9qIi0qno8BG0n38hDQ88Ru3b79P42zsIHPslXFs/Ie+3t1Jy7BGUTa2k4AcX413yFLS32x2yiEifpCNkIpKUWEkpgdO+QuC0r0BbG94Xno8fPXv6v+Te81dy7/krsbx8gofNJ/DFYwguOIpYebndYYuI9Ak6QtYLTa1B3tlYR1NrsFftBEIRttW3EghF0hRZevtJZj6rliFTuUkXO3OVyJ7baKK+k5lmr2G5uQSPOprmW39H7fr3qH/kKZouvIjg4CH4nniUou9+k7J9J5B/1Hy8N93Yq5efp/I9S1eOA6EIW7e3ZM12ncp8e86TTBuJcp6p7TZd+1WRvkZHyFIQDIe54e61bKlpJhoDpwNG+Au4+qvT8bqTT2kkGmXx0g2sq66hrjFAaZGPqko/C+dX4ErjNTmp9pPMfF1NY0fMdultrtK9THtuow4gP9eNz+Pa2ffUijI2bG7gk5qWTqeZNrEcB/Dae9s7Xy6Hg3uby1g37FjqTl7APoFtzPpgFfu+9RIVr60lb91q+MUNRIYOI7jgSIKHH0lo3mHECgp7tAzJfM/SlePd2mkKUFqYnnVl5Xexu3lKCr3k53ppbQ912kainA8vz6dyVDGvb6i1dLtN135VpK9yLVq0yO4YutXaGlyU7jbz8320pvgX2E/uWs2mbc387+/9GNDYEuT1DbV8oWpE0u3869n3WLJ6M22B+F+cbYEIH3zSSFsgzNTx6Xv+U6r9JDNfV9PMnjo85RxnKjfp0ttcpbpMnW3He26jAMFwdLe+N25torEl1OU0H25t4sOtTT1arhp3Pq/5K3l68uE8tv8xfDh4HEG3l2HbN5G7+lVyHn6A3Nt/g+elFTjraomVlBArKd3rrs1UvmfpyrFV25+V38Vu5wlGaGwJdtlGwpy3hrrdBnorP9/Hlbe/lJb9qiTWm//zJDm75jg/33ddT+fPvkMNWa6pNciWmsSvmtlS05z0YfZAKMK66pqE49ZVb0/baYFU+0lmvu6maQ+GMxqzXdKRq3QuU1fbaLoks1wAzbmFvDDpUH519CV86/v3su2Rp2n5/uWE99kX7/LnKbj2h5QefAClM6eR/8PL8Cx9BtrbU/qepSvHVq0rK7+LPZmnszZ6ut2kc7ttaA6kZb8q0pepIOuhzdvih9MTicbi45PR0BygrjGQcFx9UzsNzYnH9VSq/SQzX3fT1HcyzqqY7ZKOXKVzmbraRtMlmeXaU11LiBpjP1qvvIYdS16gdn01Tbf+jsBxJ+CorSXvz3cw6PRTKDfGUHzWQo567b/4G/cuKjr7nqUrx1atKyu/iz2Zp7M2errdpHO73bi1MS37VZG+TCfme2jk4AKcDhLuPJyO+PhkFBf4KC3yUZtgp1lSmENxga+3ofaqn2Tn62qakiIfTQ1tGYvZLunIVTqXqattNF2SWa6u5gGIDhlK+xln037G2RAM4ln5Ct4lT+Nd8hRlLz7Lt3gWuIOPykazZux01oyr4u3hk4l6PAm/Z+nabqza/qz+LiY7T2dt5PrcPdpu0rndjh1WlJb9qkhfpiNkPVSY52WEP/HOYYS/gMI8b1Lt+Dwuqir9CcdVVZbj87hSjjEd/SQzX3fT5HhTq/czlZt0SUeu0rlMXW2j6ZLMcnU2T0JeL6FDDqVl0fXUv7iK2lVv8K8TvsOqcQcwdMennLzmIW6471r+eftZ3PDYzxj8jztxbXhvtzs305Vjq9aVld/FnszTWRs93W7Sud0WF/jSsl8V6ct0UX8K5kwdwusbamluDRLj87/grv7q9B7ddTR5bAltgTANzUECwTClRTnMmTqUhfMrcKbxtTSp9pPMfF1NU1CQk3KOM5WbdOltrlJdps624z23UQdQkOumMM9DIBihtCiHWfsOJhyJ0twa6nSag6YMYfzwIhpbQkkvV0mhj/JBuXhcjp3t9HQ5Y4MGMfKYefzaM4W/7/NF3hyxL415xZSGWpjw8dt4ly4h984/kvOvf+Ay38XR3k508GAm7zMiLTm2avuz8rvY3TzJrJdE+7YR/nxmTPJ3uQ30Vn6+j/0nlKRlvyqJ6aJ+6/X2on5HLMVnA2VSTU1T2oP0+wupqWnqVRtNrUE2b2tm5ODe/QUXCEVoaA5QXOCz9OhPqv0kM1+iadKR40zlJl1SzVWqusvxnttoor6TmSaV5UrXcu4Zn3PLZrzPPYvn+aV4X3gO544dAMQcDsL77U9o3hdomTOPmsn7U1Ra1Ku+A6EILq+HSDCU1u3Pyu9id/Mk00aifZuV38Vdt+N07Vdld+nYH0vXds2x31/Y479WVJCJZZRj6w34HEciuF9fh/f5pXiWPYdn9UocoRAAsdxcQjNnE5x7GKFD5xGeOg1cPS8kBnyOM0A5tp5ybL3eFmS6qF9E+i6Xi/D0GYSnz4DvXw7NzXhfWYFn2XPxVzstew7vsucAiBYPIjRnLsG58wgdehiRiol7PftMRMQuKshEpP8oKCC44CiCC46iBXBs24b3xWV4li/Du3wZvicexffEowBEhg4jdMihhA45lOCcuUTHjLU1dBEZ2FSQiUi/FRs8mMDJpxI4+VQAnBs/xLt8GZ7lz+N98QVy7ltMzn2LAYiMGk3o4EMIzplL6JBDiY4cZWfoIjLAqCATkQEjOnYc7WPH0X72ORCN4nr3HTwvLcf74nI8L79IzuJ7yVl8LwCR0WMJHjIXvngEzn2nEx012t7gRaRfU0EmIgOT00lk8r5EJu9L+wX/Fy/Q3n4L74oX8KxYjuelFeTeew/cew9lQGTkKEKzDyZ00BxCB80hMqFC16CJSNqoIBMRgXiBNmUqbVOm0nbht+N3cL61npL1awg8sxTPKyt2O8UZLfcTOmgOwYMOJjR7DpHJ+4KelyUiKVJBJiKSiMtFeL/94fC5NJ55fvwI2nvVeF5egeeVFXheWoHv0YfwPfoQ0HEX56zZhGbPIXTQwfHHbHj1HC0RSY4KMhGRZDidRIxJRIxJtJ9zPsRiOD/aiOeVl/C8vALvyyvwPf0kvqefBCCWk0Oo6gBCsw4iPHMWoRkziQ0qsXkhRCRbqSATEUmFw0F07DgCY8cROP1MAJxbP9lZoHlWvornlZfwvrxi5yzhSfsQmnkQoZmzCM2cHX/Uhq5DExFUkImIpE102HACJ32ZwElfBsDR2IB79So8K1+OF2hrV+N+9x1y7/5LfHr/YEIzZsZ/DpxFeNr+kJtr5yKIiE1UkImIWCRWVExo/gJC8xfEB4RCuN9aj+fVl3GvWoln1av4/vsYvv8+Fp/e7SY8db94cdZRqEVHjNRRNJEBQAWZiEimeDyE959OeP/pcOG3AXBu2Yxn9Urcq+MFmnv9G3jWrYU//h6AyJChhA84kNABBxI+YAah/faHggI7l0JELKCCTETERtERIwmMGEnghJPjA9racL/xOp5VHac4V6/c7ZVPMaeTiLEPoarphKsOIDz9AMKTJoPHY+NSiEhvqSATEckmubmEZ80mPGs2bR2DnJ9swb1mFZ7Vq3CvW4Pnjddwv/MW3HsPEL+jMzx12s4iLVR1ANFx43WqU6QPUUEm3WpqDbJ5WzMjBxdQmGfdc5Uy1Y9YIxCK0NAcoLjAh8/jsjucPm+3fA4fQXD4CILHnxgfGQ7jqjbxrFuDe+0a3OvW4F67Gs+qV3fOHx00qKM4m064agbh/auIDhmq9dQNq/ZDyeTdqnUTCEXYur2FSCiCz+PSNpClMl6QGYYxAbgeWAuMBGpN0/xJpuOQ7gXDYW64ey1bapqJxsDpgBH+Aq7+6nS87vRtOpnqR6wRiUZZvHQD66prqGsMUFrko6rSz8L5Fbj05PoeSyqfbvfO1z5x5lfjw1pbO64/W437tbV41q7B+9yzeJ97dmfbzSV+Ng2ewDtl4/hs7CQK587m+JNmaT1h3X4omfVp1Xdot3abApQUeMnP9dLaHtJ3NQs5YrFYRjs0DONAYLhpmg93fH4bONs0zTWdzVNT05T2IP3+QmpqmtLdbL9y7V9Wsmlb817DRw0u4LrzZnY7f7I57m0/A1k2bMf3LqlmyerNew1fMGMkZyyotCGi9Mp0jtOZT0ddbbw4e20dnz6znGLzTcqba3ebpql0MN6ZMwjvtz/h/aYRnhY/kpZJ2bAdW7UfSmZ9WvUd6qzddPcjcbtux35/YY+vF8j44QfTNFftMcgJtGQ6DulaU2uQLTV775wAttQ009QaTMvh/Ez1I9YIhCKsq65JOG5d9XZOmTdBp0R6IN35jJWWEZp/BM1z5/PTnDnUzgtQ2lzHhM82UPHZ+1Rse5+J2z6g8Mkn8D35xM75IoOHxIuzqfsRnro/4an7ER09pt9ek2bVfiiZ9Rn/d/q/Q131nc5+JH1sPR9kGMZJwFOmab7b1XQlJXm43enfUPz+wrS32V988l4N0U6OS0Zj0BSMMn5M9/nrLsfp6mcgs3M73rq9hbqmQMJx9U3tuLwe/OX5GY4q/TKVY6vyuWu7dQWl1BXMZNWE+FEfpwP+eM6+DPnIhLVrYc0aXGvX4lryNL4lT3/eyKBBMG3a7j/77gs5OT1f0ATs3I6t2g8lsz4By9d5d/rTd9VuvdmObSvIDMP4AvAF4HvdTVtf35r2/rPhEHk2K/Q6cTpIuJNyOuLju8tfMjlORz8Dmd3bcSQUobTQR23j3jv+ksIcIsFQn19/mcyxVfnsrt1QaTk1Q4bAzEN3DnfU1uJe/zru9W/gXv8a7jdex/XCCziWLds5TczlIlIxkfC+UwhPnkpk330J7zs1fsqzB0fT7N6OrdoPJbM+gYyv83T2I5/b45Rlj+e3pSAzDONYYC5wMTDMMIwxpmm+bEcsklhhnpcR/oKE11SM8Kfv7qNM9SPW8HlcVFX6E16nUlVZrlMgPWRVPlNpN1ZWRuiw+YQOm//5wJYW3O++jfutN3G/tR73W2/ievstcsx34YH7dk4WLSsjPHlqvFDbdwrhyVOIGJPAm53fZ6v2Q8nmPdPrPJ39SPrYcVH/AcAyYHXHoHzgd6Zp3tXZPLqo3x69veso2RzrLsvUZcN2/PmdXNupb2qnpDCHqsryfnPnVqZzbFU+LVtP0SjOjz+KF2lvv7mzWHN9tHG3yWJuN5GJRkeRNnXn75jfnxXbsfV3WXae90yt80EFPvJzPbS2h6hvCvS776rdentRf8YLslSoILNXqs/l6WmO9Ryynsum7bi/PtvIrhxb+UyqTKwnR1Mjrrff3nkkzf32etzvvI2jdfdLUCKDh+Cath+tEwzCk/clMmkfwpWTIC/Psti60l+fQ+byeogEQ3oOmYVUkKUom/4j66+UY+spx9ZTjtMoGsW18QNcu5zydL/1Jq7Nm3abLOZwEB09hnClET+qZkwiMrGSSKVBrKjYpuD7Nm3H1utzj70QEZEByukkMr6CyPiKz986APi9UepfXIX7nbdwv/s2rnffwV1t4nvmKXjmqd2aiAwdRqRyEuHKSiKVk4hUGoQrJxErL8/00oiklQoyERGxV3Ex4ZmzCM+ctdtgR30drupq3NXv4qo247/fq8b7wnN4X3hut2mjpaWEKycRmWgQMQzCEw0ixiSiw4b32+enSf+igkxERLJSrKR054vWd9PcjHtDNS7zXdzvVePqKNg8K1/B+8pLu00aLSgkMnFi/KhaR5EWnlhJdMxYcOn6KckeKshERKRvKSggvP90wvtPZ7enbLW34/rg/Z1H1FzVJu73TNxvrsezbu1uTcR8PiITJhI2dr1OzSAyfkLWPp5D+jcVZCIi0j/k5Hz+0vVdhcO4PvoQl2nies/EbcZPfbrfM3G//eZuk8ZcLiKjxxCpmBi/3q1iIpEJFUQmVBAdOkynP8UyKshERKR/c7uJTJhIZMJE4LjPh0ejOLds7jiiFj/16a42cX2wAfczTwG731AQy8snPKGCyIQJRMZPIDKu42f8BGJlZSrWpFdUkImIyMDkdBIdNZrgqNFw+JG7jXLU1+F6f0P854MNuDfE/+3eUI1n/et7NRUtKiYyfjyRceN3FmmRceOJjK8gVlqqYk26pYJMRERkD7GSUsIzZhKeMXP3EdEozk+24PrwA1wfvL/L7/dxv/M2ntfW7dVWtHjQ58Xa2HFExo4nMnY80XHjiA4eomJNABVkIiIiyXM6iY4cRXTkKEJz5+0+7n/F2gfvf16sfRj/t/utN/e6sQAglptLZMzYeKE2ZhyRseOIjh0b//eo0eDzZWjBxG4qyERERNJh12Lt0MN2HxeJ4Nz6Ca6NH8YLtY7fzo0f4tr4Ie5339mruZjDQXTYcCJjxhIdMzZeuI0eQ2T0WKJjxsSPrukdlP2GCjIRERGruVyfF2uHHLr7uFgMR11d/LVSH23E9dFGnB2/XRs/xPPKSzheXrFXk7GcHCIjRxEdPSZeqI0aQ3T0aCKjRhMZNSb+9gKdDu0zVJCJiIjYyeEgVlZGuKyM8AEH7j0+EMC1ZRPOjR1F2scf4fr4I5wff4Tr4424N7yXsNlYXl5HcTYaJk4gt2xovGAbOYroqNFE/YN1hC2LqCATERHJZj7fzneAhhKMdjQ14ty0Cdemj3F9vBHnxx/j2vQxzk0f4/r4I9zmu7AECvaYL+b1EhkxkujI0URGjiQ6chSRUaOJjhgZHz5ipK5hyyAVZCIiIn1YrLAo8QNxOzgadlDeUkfDG+/i3Pwxro7izbllE65Nm3Avf77TtqP+wfFibcQoIiNGEB0+Mv55+AiiI0bqOrY0UkEmIiLSj8WKB0HFKILDxyeeoK0N15bNODdv+vz35k04t2yOPzj3rTdxJLhDFCDmdhMdNpzo8BHxgm3Y57+jI0YQGTaCmN+voi0JKshExBKBUISG5gDFBT58Hr3E2S5aD+mVaj6Tma+pNcjmbc2MHFxAYV7i92nu2U5a1m9uLq1jxtNQNoLiOQnadTlwbN9O5KOPCHywkaK6bXi3boHNm2HTx3g/+xT3qlfxvBpN2HzM4yE6dBjRocOIDB8RL+CGDSc6bBiRYSNo9w+hPr+EotLC9C1TH6SCTETSKhKNsnjpBtZV11DXGKC0yEdVpZ+F8ytw6a/kjNF6SK9U85nMfMFwmBvuXsuWmmaiMXA6YIS/gKu/Oh2v252wnZJCL/m5XlrbQ71av8m0O21iOQ7gtfeC1DWWUVI4jPyxB9I6NERdZXyaA8aXsHByPp6tW3F9+gnOTz7B+clmXFu24Nz6Cc5Pt+JeuxrPqlcTxjEMaMwrpqG4nNqCMj7LKaa6bAiFFWOYftg0GDac6NBh8VdU9dPt1xGLxeyOoVs1NU1pD9LvL6SmpindzcoulGPrZWOO711SzZLVm/cavmDGSM5YUGlDRL2TjTlORl9aD30hx6nmM5n5rv3LSjZta95rmlGDC7juvJldtpNsPJ3lONl2k9HtthWJ4KzZFj8VunUr65a9xva3P6CsuZbShcLbWAAADUJJREFU5jrKmmspa64jN9TeaRMxt5vokKFEhw4lOrjj95ChO4dFBsf/HSsvz3jhtmuO/f7CHj9vREfIRCRtAqEI66prEo5bV72dU+ZNGFCnIOyi9ZBeqeYzmfmCoQhbavYuxgC21DTT1BrE63F12k5P4ulJfKnotm+Xa+epy0Aowl0fl1Fbdsju08Ri5Abb4kVaSx1lzXWMCjVy9Gg3npptOD/divOzT3G/8TqO0OpOY4m5XET9g+OF2uBdfg8bQeDEk+PX1WUZFWQikjYNzQHqGgMJx9U3tdPQHGBwSV6Goxp4tB7SK9V8JjNfbUM70U7OAUVjsHlbM2XFOZ2205N4ehJfKtLSt8NBmy+Pzb48NpeNAuKncCd/Y/bu7UajOOrqcH72afxn22e4Ooo152ef7fztNt/B8fru7xd1NDXRdtHFvVpWK6ggE5G0KS7wUVrkozbBjrakMIfiAj3TKBO0HtIr1XwmM1+uz43TQcKizOmAkYML8HpcnbbTk3h6El8qrOo7YbtOJ7HyciLl5UT2ndL5zLFY/Dlt27bh/OxTHDt2EDp0XufT26h/XhknIrbweVxUVfoTjquqLNdpsgzRekivVPOZzHyFeV5G+Pd8ZGvcCH/8bsuu2ulJPD2JLxVW9d2rbdbhIFZUTKRiIqE5cwkeezyxwqLU2rKYa9GiRXbH0K3W1uCidLeZn++jtTWY7mZlF8qx9bIxx5PHltAWCNPQHCQQDFNalMOcqUNZOL8CZx98r1425jgZfWk99IUcp5rPZOabM3UIr2+opbk1SIzPj4xd/dXpO++Y3LOdkkIf5YNy8bgcBIKRbuPpLMfJtHvQlCGMH15EY0sopb6TzU262rXLrjnOz/dd19P5dZelWEY5tl4257i/PEsom3OcjL6wHvpSjvvqc8i6y3Ey7VryDDQL28203t5lqYJMLKMcW085tp5ybD3l2HrKsfV6W5DpGjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRERERm6kgExEREbGZCjIRkf/f3r3HyFWWcRz/7m6728vu6pZubdJKBdo+UZFaKJAGIrcaJQ0mXBKQABIkEJSGi0XAyjWA3FUS0SZCEOOFCCoqMTWgGCEN1yqI8hQoVlsbd9uu3WmXzi7T9Y/zFqfbmTK3M2935vf5pznvmTnve55z+u4z73suIiKRKSETERERiUwJmYiIiEhkSshEREREIlNCJiIiIhKZEjIRERGRyJSQiYiIiESmhExERPaQHcnRNzBEdiRXl7o2bd7xXl31rDumUvYzMzTM3/+xlczQcB1bVp2x+1VoPys5xpVst1Z118uEGJWa2RLgNKAPGHX3m2K0Q0RE/i+3axeP/P5N1qztZ+tglmndHSyc38uZJ86lrbW2v9/3qCuTpaeznamT2xnaOZJ63TGVEuPhd9/l1odfZmP/dnaNQmsLzOrtZMV5h9M+Icqf7fc1dr96uvY+ngvmTacF+PMbm0s+xpVst1Z111vL6OhoXSs0synAK8DH3T1rZo8B97v7U8W+09+fqXkje3u76O/P1HqzkkcxTp9inL5mivGPn1zLky9u2Kt8yaLZnL1kfl3qqkfdMZUS4xsefJ5/9W3f6zMfntHJTRccVVG9aZ/HpR7PQvZ1jKvZbrV1lys/xr29XS3lfj9GWrgYWO/u2bD8LLA0QjtERCTIjuRYs7a/4Lo1azfXdIpnX3WlXXdMpcQ4MzTMxv69kzGAjf3b98vpy3KOZyHFjnG1262m7hhijH3OAPLT9MFQVlRPzxQmTGireUN6e7tqvk3Zk2KcPsU4fc0Q402bd7A1ky24biCzk7b2ifROn5p6XWnXHVMpMc7sGGFXkTmhXaOQGd7FwXMqOx/TOo/LOZ6FFDvG1W63mrorVU2MYyRkfUB+i7tDWVEDA0M1b0QzTUPEohinTzFOX7PEODeSY1pXB1sG9/4D2NM1idzwSM3isK+60q47plJi3NXeSmsLBZOy1hboam+tKBZpnsflHM9Cih3jardbTd2VGDNlWfb3Y0xZrgbmmFlHWD4GeCJCO0REJOiY2MbC+b0F1y2cP52OibWbpdhXXWnXHVMpMe6a0s6s3s6Cn5nV20nXlPY0m1iRco5nIcWOcbXbrabuGOo+QubuQ2Z2CXCfmfUDr+zrgn4REamPM0+cCyTX1QxkdtLTNYmF86e/V55mXR/s7GDq5IkM7RxhIJNNte6YSonxivMOL3qX5f6qlOO5YN4B4U7HLSWfX5Vst1Z111vd77KshO6yHJ8U4/QpxulrxhhnR3Js257lA50dqY8eZEdytLVPJDc8QsfEtrrWHVMp+5kZGmZD33Zmz6h+ZKxe5/HY/Sq0n5Uc40q2W6u6S1XtXZZKyCQ1inH6FOP0KcbpU4zTpxinbzw+9kJERERE8ighExEREYlMCZmIiIhIZErIRERERCJTQiYiIiISmRIyERERkciUkImIiIhEpoRMREREJDIlZCIiIiKRjYsn9YuIiIg0Mo2QiYiIiESmhExEREQkMiVkIiIiIpEpIRMRERGJTAmZiIiISGRKyEREREQimxC7AfVmZkuA04A+YNTdb4rcpHHPzA4BbgFeBmYDW9z9ZjObBtwOrAPmAV9z9//Ea+n4Z2aTgeeA37n7cjObBNwNbCSJ8e3uvjZmG8czMzPg88A7wHHAjSR9xXXAm8BHgK+4+/ZITRz3zOwqkjhuJjlnvwhMRn1FVcxsJkk/vMDdjwxlRfsHMzsHWAjkgLfcfWWUho8jRWJ8NTAT2AQsAq5399fDurJi3FQjZGY2BfgecIW73wgcZmYnxW1VQ5gG/NTd73L3y4CzzOwI4DbgSXe/HfglSccg1bkFWJO3fDnwT3f/BvBN4IEorWoAZtYG3Avc7O53kCQKb5P0GStDjP8KXB2vleNb+IN2LbDM3W8AppL8QFZfUb1jgceBlryygv2Dmc0GlgPL3f2rwIVmNq/O7R2PCsW4E7jS3e8EHgPugspi3FQJGbAYWO/u2bD8LLA0Ynsagru/4O6P5xW1AjtIYrs6lCnWVTKzc0ni+HZe8XsxdvdXgQVm1h2heY3gSJKOdpmZXQucAvwXOAF4IXxG53F1hoBhYPc52gm8hvqKqrn7o0BmTHGx/uEzwEvuvvvJ8KuBk+vV1vGqUIzd/bq8OLYCu0fPy45xsyVkM9gzmIOhTGrEzE4FVoUh2/x4DwI9ZtZ00+S1YGYfAz7q7j8fs0rndO3MIfnR9lAYUfgUyS/cd/I6VcW3Cu4+CFwFPGJmDwEbSKaC1Veko1j/oH6jxsysHfgC8PVQVHaMmy0h6wO68pa7Q5nUgJmdQDKacEUoyo93NzDg7u/GaFsDOBXYaWbXkAybH2Vml6NzupYGgdfdfVtYfgY4FJhsZrunKBTfKpjZJ0kSsqXufj7JdWTXo74iLcX6B/UbNRSSse8CK9z9rVBcdoybLSFbDcwxs46wfAzwRMT2NAwzW0oyRHsZMNPMFpPEdnH4iGJdBXe/1d1vDtfYPAM87+7fIi/GZvYJ4C9hFELK9xxwQLiWDJIRs9eAP5BMZ4LO42rNArbmJVubgEmor0hLsf5hFXBE3g+NxcBv4zRxfAvXpq8E7nX3l8zs9LCq7Bg33cvFzezTwBlAPzCiuyyrFy7g/yPwYiiaCnwH+BVwB7AeOAS4RndOVSf8Z/8y0E4S490XQG8C5gK36S7LyoUp9xNJ+ocDgWXAh0hGcdaFsit1l2VlQrJ7H7CT5Pq8Q0kuPM+ivqIqZnYccB7wWZLRmnvCqoL9Q7gDcBHJHYBrdZfl+ysS4x+RnMf/Dh+bmncHZlkxbrqETERERGR/02xTliIiIiL7HSVkIiIiIpEpIRMRERGJTAmZiIiISGRKyEREREQi05OQRaQhhXfJ3Q90u/vxkZsjIrJPGiETkYbk7htIXhYuIrLf0wiZiDQUM1sGnA44yeuQCO9F/EUomwRscfcbwhslfkDyUMcLSF52/TDJE86/TfLw3XUk76B7zd3vru/eiEiz0AiZiDQMMzsMWAGc7O4XkzwBfreH3H25u18KLDKzo919NXAjyQvE1wGvAn9z9xXAOUC/u18NXIj6SxFJkUbIRKSRHA+86O7vhOU/kbyMPQccaGYPABngIGA+yfsrHwXuMbO5ofyp8N2ngavMrBv4GZr+FJEU6RefiDSDs0jeQXeRu18OPA+0Abj7MPAg8CXgbOAnofxl4GDg18CVwG/q32wRaRZKyESkkTxNMh05OSwfG/49ANjm7rmwfOCY760Ezgc2u/sQgJldBBzk7o8AnwOOTrHdItLk9HJxEWko4aL+M0iuB5sEnAzcCSwFtgHrgZOALcAl7v5G+N4q4NK85VOAi4HXgZnAKnf/YX33RkSahRIyEWlqZnYwsBH4vrufG7s9ItKcdFG/iDS7x0keh3FH7IaISPPSCJmIiIhIZLqoX0RERCQyJWQiIiIikSkhExEREYlMCZmIiIhIZErIRERERCJTQiYiIiIS2f8ANCqs+SxIrzkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1,1,figsize=(10, 6))\n",
    "ax.scatter(df.days, df.students)\n",
    "ax.plot(df.days, yhat, color='red', label='fitted values')\n",
    "ax.set_xlabel('days')\n",
    "ax.set_ylabel('students')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}