exercice 1 module 4 reproduire les résultats de l'analyse challenger dans le...

exercice 1 module 4 reproduire les résultats de l'analyse challenger dans le langage pas utilisé (python)

module4/src_Python3_challenger.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this document we reperform some of the analysis provided in \n",
+    "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n",
+    "\n",
+    "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Technical information on the computer on which the analysis is run"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "3.6.4 |Anaconda, Inc.| (default, Jan 16 2018, 18:10:19) \n",
+      "[GCC 7.2.0]\n",
+      "uname_result(system='Linux', node='3a716011d2b6', release='4.4.0-116-generic', version='#140-Ubuntu SMP Mon Feb 12 21:23:04 UTC 2018', machine='x86_64', processor='x86_64')\n",
+      "IPython 6.4.0\n",
+      "IPython.core.release 6.4.0\n",
+      "PIL 5.2.0\n",
+      "PIL.Image 5.2.0\n",
+      "PIL._version 5.2.0\n",
+      "_csv 1.0\n",
+      "_ctypes 1.1.0\n",
+      "_curses b'2.2'\n",
+      "decimal 1.70\n",
+      "argparse 1.1\n",
+      "backcall 0.1.0\n",
+      "cffi 1.11.5\n",
+      "csv 1.0\n",
+      "ctypes 1.1.0\n",
+      "cycler 0.10.0\n",
+      "dateutil 2.7.3\n",
+      "decimal 1.70\n",
+      "decorator 4.3.0\n",
+      "distutils 3.6.4\n",
+      "ipaddress 1.0\n",
+      "ipykernel 4.8.2\n",
+      "ipykernel._version 4.8.2\n",
+      "ipython_genutils 0.2.0\n",
+      "ipython_genutils._version 0.2.0\n",
+      "ipywidgets 7.2.1\n",
+      "ipywidgets._version 7.2.1\n",
+      "jedi 0.12.1\n",
+      "json 2.0.9\n",
+      "jupyter_client 5.2.3\n",
+      "jupyter_client._version 5.2.3\n",
+      "jupyter_core 4.4.0\n",
+      "jupyter_core.version 4.4.0\n",
+      "kiwisolver 1.0.1\n",
+      "logging 0.5.1.2\n",
+      "matplotlib 2.2.2\n",
+      "matplotlib.backends.backend_agg 2.2.2\n",
+      "numpy 1.13.3\n",
+      "numpy.core 1.13.3\n",
+      "numpy.core.multiarray 3.1\n",
+      "numpy.core.umath b'0.4.0'\n",
+      "numpy.lib 1.13.3\n",
+      "numpy.linalg._umath_linalg b'0.1.5'\n",
+      "numpy.matlib 1.13.3\n",
+      "optparse 1.5.3\n",
+      "pandas 0.22.0\n",
+      "_libjson 1.33\n",
+      "parso 0.3.0\n",
+      "patsy 0.5.0\n",
+      "patsy.version 0.5.0\n",
+      "pexpect 4.6.0\n",
+      "pickleshare 0.7.4\n",
+      "platform 1.0.8\n",
+      "prompt_toolkit 1.0.15\n",
+      "ptyprocess 0.6.0\n",
+      "pygments 2.2.0\n",
+      "pyparsing 2.2.0\n",
+      "pytz 2018.5\n",
+      "re 2.2.1\n",
+      "scipy 1.1.0\n",
+      "scipy._lib.decorator 4.0.5\n",
+      "scipy._lib.six 1.2.0\n",
+      "scipy.fftpack._fftpack b'$Revision: $'\n",
+      "scipy.fftpack.convolve b'$Revision: $'\n",
+      "scipy.integrate._dop b'$Revision: $'\n",
+      "scipy.integrate._ode $Id$\n",
+      "scipy.integrate._odepack  1.9 \n",
+      "scipy.integrate._quadpack  1.13 \n",
+      "scipy.integrate.lsoda b'$Revision: $'\n",
+      "scipy.integrate.vode b'$Revision: $'\n",
+      "scipy.interpolate._fitpack  1.7 \n",
+      "scipy.interpolate.dfitpack b'$Revision: $'\n",
+      "scipy.linalg 0.4.9\n",
+      "scipy.linalg._fblas b'$Revision: $'\n",
+      "scipy.linalg._flapack b'$Revision: $'\n",
+      "scipy.linalg._flinalg b'$Revision: $'\n",
+      "scipy.ndimage 2.0\n",
+      "scipy.optimize._cobyla b'$Revision: $'\n",
+      "scipy.optimize._lbfgsb b'$Revision: $'\n",
+      "scipy.optimize._minpack  1.10 \n",
+      "scipy.optimize._nnls b'$Revision: $'\n",
+      "scipy.optimize._slsqp b'$Revision: $'\n",
+      "scipy.optimize.minpack2 b'$Revision: $'\n",
+      "scipy.signal.spline 0.2\n",
+      "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n",
+      "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n",
+      "scipy.special.specfun b'$Revision: $'\n",
+      "scipy.stats.mvn b'$Revision: $'\n",
+      "scipy.stats.statlib b'$Revision: $'\n",
+      "seaborn 0.8.1\n",
+      "seaborn.external.husl 2.1.0\n",
+      "seaborn.external.six 1.10.0\n",
+      "six 1.11.0\n",
+      "statsmodels 0.9.0\n",
+      "statsmodels.__init__ 0.9.0\n",
+      "traitlets 4.3.2\n",
+      "traitlets._version 4.3.2\n",
+      "urllib.request 3.6\n",
+      "zlib 1.0\n",
+      "zmq 17.1.0\n",
+      "zmq.sugar 17.1.0\n",
+      "zmq.sugar.version 17.1.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "def print_imported_modules():\n",
+    "    import sys\n",
+    "    for name, val in sorted(sys.modules.items()):\n",
+    "        if(hasattr(val, '__version__')): \n",
+    "            print(val.__name__, val.__version__)\n",
+    "#        else:\n",
+    "#            print(val.__name__, \"(unknown version)\")\n",
+    "def print_sys_info():\n",
+    "    import sys\n",
+    "    import platform\n",
+    "    print(sys.version)\n",
+    "    print(platform.uname())\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import statsmodels.api as sm\n",
+    "import seaborn as sns\n",
+    "\n",
+    "print_sys_info()\n",
+    "print_imported_modules()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Loading and inspecting data\n",
+    "Let's start by reading data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>Date</th>\n",
+       "      <th>Count</th>\n",
+       "      <th>Temperature</th>\n",
+       "      <th>Pressure</th>\n",
+       "      <th>Malfunction</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>4/12/81</td>\n",
+       "      <td>6</td>\n",
+       "      <td>66</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>11/12/81</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>50</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>3/22/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>69</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>11/11/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>68</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>4/04/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>6/18/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>72</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>6</th>\n",
+       "      <td>8/30/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>73</td>\n",
+       "      <td>100</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>7</th>\n",
+       "      <td>11/28/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>100</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>8</th>\n",
+       "      <td>2/03/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>57</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9</th>\n",
+       "      <td>4/06/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>63</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>10</th>\n",
+       "      <td>8/30/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>11</th>\n",
+       "      <td>10/05/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>78</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>12</th>\n",
+       "      <td>11/08/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>13</th>\n",
+       "      <td>1/24/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>53</td>\n",
+       "      <td>200</td>\n",
+       "      <td>2</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>14</th>\n",
+       "      <td>4/12/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>15</th>\n",
+       "      <td>4/29/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>75</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>16</th>\n",
+       "      <td>6/17/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>17</th>\n",
+       "      <td>7/2903/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>81</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>18</th>\n",
+       "      <td>8/27/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>76</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>19</th>\n",
+       "      <td>10/03/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>79</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>20</th>\n",
+       "      <td>10/30/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>75</td>\n",
+       "      <td>200</td>\n",
+       "      <td>2</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>21</th>\n",
+       "      <td>11/26/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>76</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>22</th>\n",
+       "      <td>1/12/86</td>\n",
+       "      <td>6</td>\n",
+       "      <td>58</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "         Date  Count  Temperature  Pressure  Malfunction\n",
+       "0     4/12/81      6           66        50            0\n",
+       "1    11/12/81      6           70        50            1\n",
+       "2     3/22/82      6           69        50            0\n",
+       "3    11/11/82      6           68        50            0\n",
+       "4     4/04/83      6           67        50            0\n",
+       "5     6/18/82      6           72        50            0\n",
+       "6     8/30/83      6           73       100            0\n",
+       "7    11/28/83      6           70       100            0\n",
+       "8     2/03/84      6           57       200            1\n",
+       "9     4/06/84      6           63       200            1\n",
+       "10    8/30/84      6           70       200            1\n",
+       "11   10/05/84      6           78       200            0\n",
+       "12   11/08/84      6           67       200            0\n",
+       "13    1/24/85      6           53       200            2\n",
+       "14    4/12/85      6           67       200            0\n",
+       "15    4/29/85      6           75       200            0\n",
+       "16    6/17/85      6           70       200            0\n",
+       "17  7/2903/85      6           81       200            0\n",
+       "18    8/27/85      6           76       200            0\n",
+       "19   10/03/85      6           79       200            0\n",
+       "20   10/30/85      6           75       200            2\n",
+       "21   11/26/85      6           76       200            0\n",
+       "22    1/12/86      6           58       200            1"
+      ]
+     },
+     "execution_count": 2,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")\n",
+    "data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBNJREFUeJzt3XuQnXWd5/H3t5MACYmAxMkwBAQGhpUCJkDLReaSCFqBKpN1AYUtwXEGM1uSskRHYWZdhmGdqpVRmXFlFGRxhC2NXEbIzmaWixAdprgFiOEmTA8gdEDAGCANIemkv/vHefrxpOnLOZ1++vQ5vF9VqZznOb9++vvtp09/+rn070RmIkkSQFerC5AkTR2GgiSpZChIkkqGgiSpZChIkkqGgiSpVFkoRMTVEfFSRDwywvMREV+PiJ6IWBcRR1dViySpMVUeKfwDsHiU508BDin+LQO+WWEtkqQGVBYKmfkT4FejDFkKXJM19wB7RsQ+VdUjSRrb9BZ+7n2B5+qWe4t1LwwdGBHLqB1NMHPmzGP222+/SSmwUQMDA3R1dd7lmU7tCzq3N/tqP5PV25NPPvnLzHzXWONaGQoxzLph59zIzCuBKwG6u7tzzZo1VdbVtNWrV7Nw4cJWlzHhOrUv6Nze7Kv9TFZvEfHzRsa1Mnp7gfpf+ecDz7eoFkkSrQ2FlcA5xV1IxwOvZuZbTh1JkiZPZaePIuL7wEJgbkT0An8JzADIzG8Bq4BTgR7gDeATVdUiSWpMZaGQmWeN8XwC51X1+SVJzevMy/mSpHExFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklSqNBQiYnFEPBERPRFx4TDP7x8Rd0bEQxGxLiJOrbIeSdLoKguFiJgGXA6cAhwGnBURhw0Z9kXgusw8CjgT+Puq6pEkja3KI4VjgZ7MfCoztwIrgKVDxiTwjuLxHsDzFdYjSRpDZGY1G444HVicmecWy2cDx2Xm8rox+wC3AnsBuwMnZ+YDw2xrGbAMYN68ecesWLGikprHq6+vj9mzZ7e6jAnXqX1B5/ZmX+1nsnpbtGjRA5nZPda46RXWEMOsG5pAZwH/kJlfjYgTgGsj4vDMHNjhgzKvBK4E6O7uzoULF1ZR77itXr2aqVbTROjUvqBze7Ov9jPVeqvy9FEvsF/d8nzeenroT4DrADLzbmA3YG6FNUmSRlFlKNwPHBIRB0bELtQuJK8cMuZZ4CSAiHgPtVB4ucKaJEmjqCwUMnMbsBy4BXic2l1Gj0bEJRGxpBj2OeCTEfFT4PvAH2VVFzkkSWOq8poCmbkKWDVk3UV1jx8DTqyyBklS4/yLZklSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUqDYWIWBwRT0RET0RcOMKYj0TEYxHxaER8r8p6JEmjm97IoIg4PDMfaWbDETENuBz4ANAL3B8RKzPzsboxhwB/DpyYmRsj4jea+RySpInV6JHCtyLivoj4VETs2eDHHAv0ZOZTmbkVWAEsHTLmk8DlmbkRIDNfanDbkqQKRGY2NrD2W/0fA2cA9wHfyczbRhl/OrA4M88tls8GjsvM5XVjbgKeBE4EpgEXZ+b/G2Zby4BlAPPmzTtmxYoVjXU3Sfr6+pg9e3ary5hwndoXdG5v9tV+Jqu3RYsWPZCZ3WMOzMyG/1H7wX0asB54HPgZ8J9GGHsGcFXd8tnA/xwy5p+AHwIzgAOpnWbac7QajjnmmJxq7rzzzlaXUIlO7Suzc3uzr/YzWb0Ba7KBn/MNnT6KiCMj4rIiCN4PfCgz31M8vmyED+sF9qtbng88P8yYmzOzPzOfBp4ADmmkJknSxGv0msI3gAeB383M8zLzQYDMfB744ggfcz9wSEQcGBG7AGcCK4eMuQlYBBARc4HfAZ5qrgVJ0kRp6O4j4FRgc2ZuB4iILmC3zHwjM68d7gMyc1tELAduoXba6erMfDQiLqF2GLOyeO6DEfEYsB34fGZu2MmeJEnj1Ggo3A6cDPQVy7OAW4H3jfZBmbkKWDVk3UV1jxP4bPFPktRijZ4+2i0zBwOB4vGsakqSJLVKo6HwekQcPbgQEccAm6spSZLUKo2ePvoMcH1EDN49tA/w0WpKkiS1SkOhkJn3R8R/AA4FAvhZZvZXWpkkadI1eqQA8F7ggOJjjooIMvOaSqqSJLVEoxPiXQv8NrCW2q2jAAkYCpLUQRo9UugGDituIZUkdahG7z56BPjNKguRJLVeo0cKc4HHIuI+YMvgysxcUklVkqSWaDQULq6yCEnS1NDoLak/joh3A4dk5u0RMYvafEaSpA7S6NTZnwRuAK4oVu1LbYZTSVIHafRC83nU3h3tNYDM/DfA91OWpA7TaChsydr7LAMQEdOp/Z2CJKmDNBoKP46IvwBmRsQHgOuB/1NdWZKkVmg0FC4EXgYeBv6U2nskjPSOa5KkNtXo3UcDwLeLf5KkDtXo3EdPM8w1hMw8aMIrkiS1TDNzHw3aDTgDeOfElyNJaqWGrilk5oa6f+sz82+B91dcmyRpkjV6+ujousUuakcOcyqpSJLUMo2ePvpq3eNtwDPARya8GklSSzV699GiqguRJLVeo6ePPjva85n5tYkpR5LUSs3cffReYGWx/CHgJ8BzVRQlSWqNZt5k5+jM3AQQERcD12fmuVUVJkmafI1Oc7E/sLVueStwwIRXI0lqqUaPFK4F7ouIH1L7y+YPA9dUVpUkqSUavfvoryPin4HfL1Z9IjMfqq4sSVIrNHr6CGAW8Fpm/h3QGxEHVlSTJKlFGn07zr8ELgD+vFg1A/jfVRUlSWqNRo8UPgwsAV4HyMzncZoLSeo4jYbC1sxMiumzI2L36kqSJLVKo6FwXURcAewZEZ8Ebsc33JGkjtPo3UdfKd6b+TXgUOCizLyt0sokSZNuzCOFiJgWEbdn5m2Z+fnM/LNGAyEiFkfEExHRExEXjjLu9IjIiOgeaYwkqXpjhkJmbgfeiIg9mtlwREwDLgdOAQ4DzoqIw4YZNwf4NHBvM9uXJE28Rv+i+U3g4Yi4jeIOJIDM/PQoH3Ms0JOZTwFExApgKfDYkHH/HbgU+LNGi5YkVaPRUPi/xb9m7MuOs6j2AsfVD4iIo4D9MvOfImLEUIiIZcAygHnz5rF69eomS6lWX1/flKtpInRqX9C5vdlX+5lqvY0aChGxf2Y+m5nfHce2Y5h1WbftLuAy4I/G2lBmXglcCdDd3Z0LFy4cRznVWb16NVOtponQqX1B5/ZmX+1nqvU21jWFmwYfRMSNTW67F9ivbnk+8Hzd8hzgcGB1RDwDHA+s9GKzJLXOWKFQ/9v+QU1u+37gkIg4MCJ2Ac7k12/SQ2a+mplzM/OAzDwAuAdYkplrmvw8kqQJMlYo5AiPx5SZ24DlwC3A48B1mfloRFwSEUuaK1OSNBnGutD8uxHxGrUjhpnFY4rlzMx3jPbBmbkKWDVk3UUjjF3YUMWSpMqMGgqZOW2yCpEktV4z76cgSepwhoIkqWQoSJJKhoIkqfS2CYUNfVv46XOvsKFvS6tLkdSEDX1b2Ny/3dfuJHlbhMLNa9dz4pfv4GNX3cuJX76DlWvXt7okSQ0YfO0+/fLrvnYnSceHwoa+LVxw4zre7B9g05ZtvNk/wBduXOdvHdIUV//a3Z7pa3eSdHwo9G7czIyuHduc0dVF78bNLapIUiN87bZGx4fC/L1m0j8wsMO6/oEB5u81s0UVSWqEr93W6PhQ2Hv2rlx62pHsNqOLObtOZ7cZXVx62pHsPXvXVpcmaRT1r91pEb52J0mjb7LT1pYs2JcTD55L78bNzN9rpt9UUpsYfO3ed/dd/OuS3/O1OwneFqEAtd86/IaS2s/es3dl5oxpvn4nScefPpIkNc5QkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUqnSUIiIxRHxRET0RMSFwzz/2Yh4LCLWRcSPIuLdVdYjSRpdZaEQEdOAy4FTgMOAsyLisCHDHgK6M/NI4Abg0qrqkSSNrcojhWOBnsx8KjO3AiuApfUDMvPOzHyjWLwHmF9hPZKkMURmVrPhiNOBxZl5brF8NnBcZi4fYfw3gF9k5peGeW4ZsAxg3rx5x6xYsaKSmserr6+P2bNnt7qMCdepfUHn9mZf7Weyelu0aNEDmdk91rjpFdYQw6wbNoEi4mNAN/CHwz2fmVcCVwJ0d3fnwoULJ6jEibF69WqmWk0ToVP7gs7tzb7az1TrrcpQ6AX2q1ueDzw/dFBEnAz8V+APM3NLhfVIksZQ5TWF+4FDIuLAiNgFOBNYWT8gIo4CrgCWZOZLFdYiSWpAZaGQmduA5cAtwOPAdZn5aERcEhFLimF/A8wGro+ItRGxcoTNSZImQZWnj8jMVcCqIesuqnt8cpWfv51t6NtC78bNzN9rJnvP3nXCxraTTu2rKj0vbmLjG/30vLiJg+fNaXU5alOVhoLG5+a167ngxnXM6Oqif2CAS087kiUL9t3pse2kU/uqykU3Pcw19zzL547YxvmX/YRzTtifS5Ye0eqy1Iac5mKK2dC3hQtuXMeb/QNs2rKNN/sH+MKN69jQ99Zr8M2MbSed2ldVel7cxDX3PLvDumvufpaeFze1qCK1M0NhiunduJkZXTvulhldXfRu3LxTY9tJp/ZVlbXPvdLUemk0hsIUM3+vmfQPDOywrn9ggPl7zdypse2kU/uqyoL99mxqvTQaQ2GK2Xv2rlx62pHsNqOLObtOZ7cZXVx62pHDXmhtZmw76dS+qnLwvDmcc8L+O6w754T9vdiscfFC8xS0ZMG+nHjw3IbuvGlmbDvp1L6qcsnSIzjn+AN4+IF7uP384w0EjZuhMEXtPXvXhn8QNjO2nXRqX1U5eN4cemfNMBC0Uzx9JEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpFKloRARiyPiiYjoiYgLh3l+14j4QfH8vRFxQJX1SM3a0LeFnz73Chv6tow6bs3TG/jarU+w5ukNE7bNZsf2vLiJjW/00/PipjHHNqOqepv5/Jv7tzf8NbhhzXMd9zWocrtDTa9qwxExDbgc+ADQC9wfESsz87G6YX8CbMzMgyPiTODLwEerqklqxs1r13PBjeuY0dVF/8AAl552JEsW7PuWcR+76h7u6qmFwdfv6OH3D96ba889fqe22ezYi256mGvueZbPHbGN8y/7CeecsD+XLD1inJ1XX2+zn//T7+nn/C/f0dDXYFCnfA2q3O5wqjxSOBboycynMnMrsAJYOmTMUuC7xeMbgJMiIiqsSWrIhr4tXHDjOt7sH2DTlm282T/AF25c95bf0tY8vaEMhEH/0rNh2COGRrfZ7NieFzft8MMQ4Jq7n93p35arqnc8n3975tvya1DldkcSmVnNhiNOBxZn5rnF8tnAcZm5vG7MI8WY3mL534sxvxyyrWXAsmLxUOCJSooev7nAL8cc1X46tS8Yo7eYMXPW9L32+Z3o6po2uC4HBrZv2/jCk9m/+Y3BddPmzP2tabvvuc/Qj9/++isvbN/0y+fHs81mx3bN2mPv6e941wEA2994lWmz9gBg22svPzPwxqtjn8/aya9Bs2PH8/kH+2rka1CvTb4GE/K92IB3Z+a7xhpU2ekjYLjf+IcmUCNjyMwrgSsnoqgqRMSazOxudR0TrVP7gs7tLSLWbHv1JftqI1Pte7HK00e9wH51y/OB50caExHTgT2AX1VYkyRpFFWGwv3AIRFxYETsApwJrBwyZiXw8eLx6cAdWdX5LEnSmCo7fZSZ2yJiOXALMA24OjMfjYhLgDWZuRL4X8C1EdFD7QjhzKrqqdiUPbW1kzq1L+jc3uyr/Uyp3iq70CxJaj/+RbMkqWQoSJJKhsI4RMQzEfFwRKyNiDXFuosjYn2xbm1EnNrqOpsVEXtGxA0R8bOIeDwiToiId0bEbRHxb8X/e7W6zmaN0Fcn7K9D6+pfGxGvRcRn2n2fjdJXJ+yz8yPi0Yh4JCK+HxG7FTfj3Fvsrx8UN+a0rkavKTQvIp4Buuv/yC4iLgb6MvMrraprZ0XEd4F/ycyrim/MWcBfAL/KzP9RzF+1V2Ze0NJCmzRCX5+hzfdXvWJamfXAccB5tPk+GzSkr0/QxvssIvYF7gIOy8zNEXEdsAo4FfjHzFwREd8CfpqZ32xVnR4pCICIeAfwB9TuCCMzt2bmK+w4Fcl3gf/YmgrHZ5S+Os1JwL9n5s9p8302RH1fnWA6MLP4u6xZwAvA+6lN8wNTYH8ZCuOTwK0R8UAxBceg5RGxLiKubrdDduAg4GXgOxHxUERcFRG7A/My8wWA4v/faGWR4zBSX9De+2uoM4HvF4/bfZ/Vq+8L2nifZeZ64CvAs9TC4FXgAeCVzNxWDOsFqpnprkGGwvicmJlHA6cA50XEHwDfBH4bWEBth3+1hfWNx3TgaOCbmXkU8DrwlunO29BIfbX7/ioVp8SWANe3upaJNExfbb3PihBbChwI/BawO7WfIUO19Jy+oTAOmfl88f9LwA+BYzPzxczcnpkDwLepzRLbTnqB3sy8t1i+gdoP0xcjYh+A4v+XWlTfeA3bVwfsr3qnAA9m5ovFcrvvs0E79NUB++xk4OnMfDkz+4F/BN4H7FmcToLhpwOaVIZCkyJi94iYM/gY+CDwyOCLsPBh4JFW1DdemfkL4LmIOLRYdRLwGDtORfJx4OYWlDduI/XV7vtriLPY8RRLW++zOjv01QH77Fng+IiYFRHBr19jd1Kb5gemwP7y7qMmRcRB1I4OoHZq4nuZ+dcRcS21w9oEngH+dPC8bruIiAXAVcAuwFPU7vboAq4D9qf2TX1GZrbVpIUj9PV12nx/AUTELOA54KDMfLVYtzftv8+G66sTXmN/Re2NxLYBDwHnUruGsAJ4Z7HuY5lZ7durjVajoSBJGuTpI0lSyVCQJJUMBUlSyVCQJJUMBUlSqbJ3XpMmW3Er5o+Kxd8EtlOb4gJqf2C4tSWFjSIi/hhYVfw9hdRy3pKqjjSVZq2NiGmZuX2E5+4Clmfm2ia2N71urhxpQnn6SG8LEfHxiLivmIf/7yOiKyKmR8QrEfE3EfFgRNwSEcdFxI8j4qnB+foj4tyI+GHx/BMR8cUGt/uliLgPODYi/ioi7i/m0f9W1HyU2h9j/aD4+F0iojci9iy2fXxE3F48/lJEXBERt1Gb3G96RHyt+NzrIuLcyf+qqhMZCup4EXE4tWkR3peZC6idNj2zeHoP4NZigsOtwMXUph84A7ikbjPHFh9zNPCfI2JBA9t9MDOPzcy7gb/LzPcCRxTPLc7MHwBrgY9m5oIGTm8dBXwoM88GlgEvZeaxwHupTcy4/3i+PlI9ryno7eBkaj8419SmnGEmtSkUADZn5m3F44eBVzNzW0Q8DBxQt41bMnMjQETcBPwetdfPSNvdyq+nQwE4KSI+D+wGzKU2ZfI/N9nHzZn5ZvH4g8B7IqI+hA6hNq2FNG6Ggt4OArg6M//bDitrM1PW/3Y+AGype1z/+hh68S3H2O7mLC7YFfP4fIPa7KzrI+JL1MJhONv49RH80DGvD+npU5n5I6QJ5OkjvR3cDnwkIuZC7S6lcZxq+WDU3ut5FrU58f+1ie3OpBYyvyxm2D2t7rlNwJy65WeAY4rH9eOGugX41OCUy1F7X+OZTfYkvYVHCup4mflwMTvl7RHRBfQD/4Xm5q2/C/getTd5uXbwbqFGtpuZG6L2PtGPAD8H7q17+jvAVRGxmdp1i4uBb0fEL4D7RqnnCmqzoK4tTl29RC2spJ3iLanSGIo7ew7PzM+0uhapap4+kiSVPFKQJJU8UpAklQwFSVLJUJAklQwFSVLJUJAklf4/GKF1l7kqzyEAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "%matplotlib inline\n",
+    "pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "data[\"Frequency\"]=data.Malfunction/data.Count\n",
+    "data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n",
+    "plt.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Logistic regression\n",
+    "\n",
+    "Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<table class=\"simpletable\">\n",
+       "<caption>Generalized Linear Model Regression Results</caption>\n",
+       "<tr>\n",
+       "  <th>Dep. Variable:</th>      <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model Family:</th>       <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Link Function:</th>        <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -3.9210</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Date:</th>           <td>Wed, 24 Oct 2018</td> <th>  Deviance:          </th> <td>  3.0144</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Time:</th>               <td>11:05:55</td>     <th>  Pearson chi2:      </th>  <td>  5.00</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>No. Iterations:</th>         <td>6</td>        <th>  Covariance Type:   </th> <td>nonrobust</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>    5.0850</td> <td>    7.477</td> <td>    0.680</td> <td> 0.496</td> <td>   -9.570</td> <td>   19.740</td>\n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Temperature</th> <td>   -0.1156</td> <td>    0.115</td> <td>   -1.004</td> <td> 0.316</td> <td>   -0.341</td> <td>    0.110</td>\n",
+       "</tr>\n",
+       "</table>"
+      ],
+      "text/plain": [
+       "<class 'statsmodels.iolib.summary.Summary'>\n",
+       "\"\"\"\n",
+       "                 Generalized Linear Model Regression Results                  \n",
+       "==============================================================================\n",
+       "Dep. Variable:              Frequency   No. Observations:                   23\n",
+       "Model:                            GLM   Df Residuals:                       21\n",
+       "Model Family:                Binomial   Df Model:                            1\n",
+       "Link Function:                  logit   Scale:                          1.0000\n",
+       "Method:                          IRLS   Log-Likelihood:                -3.9210\n",
+       "Date:                Wed, 24 Oct 2018   Deviance:                       3.0144\n",
+       "Time:                        11:05:55   Pearson chi2:                     5.00\n",
+       "No. Iterations:                     6   Covariance Type:             nonrobust\n",
+       "===============================================================================\n",
+       "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
+       "-------------------------------------------------------------------------------\n",
+       "Intercept       5.0850      7.477      0.680      0.496      -9.570      19.740\n",
+       "Temperature    -0.1156      0.115     -1.004      0.316      -0.341       0.110\n",
+       "===============================================================================\n",
+       "\"\"\""
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import statsmodels.api as sm\n",
+    "\n",
+    "data[\"Success\"]=data.Count-data.Malfunction\n",
+    "data[\"Intercept\"]=1\n",
+    "\n",
+    "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n",
+    "                family=sm.families.Binomial(sm.families.links.logit)).fit()\n",
+    "\n",
+    "logmodel.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<table class=\"simpletable\">\n",
+       "<caption>Generalized Linear Model Regression Results</caption>\n",
+       "<tr>\n",
+       "  <th>Dep. Variable:</th>      <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model Family:</th>       <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Link Function:</th>        <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -23.526</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Date:</th>           <td>Wed, 24 Oct 2018</td> <th>  Deviance:          </th> <td>  18.086</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Time:</th>               <td>11:05:55</td>     <th>  Pearson chi2:      </th>  <td>  30.0</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>No. Iterations:</th>         <td>6</td>        <th>  Covariance Type:   </th> <td>nonrobust</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>    5.0850</td> <td>    3.052</td> <td>    1.666</td> <td> 0.096</td> <td>   -0.898</td> <td>   11.068</td>\n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Temperature</th> <td>   -0.1156</td> <td>    0.047</td> <td>   -2.458</td> <td> 0.014</td> <td>   -0.208</td> <td>   -0.023</td>\n",
+       "</tr>\n",
+       "</table>"
+      ],
+      "text/plain": [
+       "<class 'statsmodels.iolib.summary.Summary'>\n",
+       "\"\"\"\n",
+       "                 Generalized Linear Model Regression Results                  \n",
+       "==============================================================================\n",
+       "Dep. Variable:              Frequency   No. Observations:                   23\n",
+       "Model:                            GLM   Df Residuals:                       21\n",
+       "Model Family:                Binomial   Df Model:                            1\n",
+       "Link Function:                  logit   Scale:                          1.0000\n",
+       "Method:                          IRLS   Log-Likelihood:                -23.526\n",
+       "Date:                Wed, 24 Oct 2018   Deviance:                       18.086\n",
+       "Time:                        11:05:55   Pearson chi2:                     30.0\n",
+       "No. Iterations:                     6   Covariance Type:             nonrobust\n",
+       "===============================================================================\n",
+       "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
+       "-------------------------------------------------------------------------------\n",
+       "Intercept       5.0850      3.052      1.666      0.096      -0.898      11.068\n",
+       "Temperature    -0.1156      0.047     -2.458      0.014      -0.208      -0.023\n",
+       "===============================================================================\n",
+       "\"\"\""
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n",
+    "                family=sm.families.Binomial(sm.families.links.logit),\n",
+    "                var_weights=data['Count']).fit()\n",
+    "\n",
+    "logmodel.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\n",
+    "The Goodness of fit (Deviance) indicated for this model is $G^2=18.086$ with 21 degrees of freedom (Df Residuals).\n",
+    "\n",
+    "**I have therefore managed to fully replicate the results of the Dalal *et al.* article**."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Predicting failure probability\n",
+    "The temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOXZ//HPNUs2sgABwhI2NYDIngUQa8EqoFXcUEDEpSD2qUutlVb6WLVWuzz0+blXoYBrFakVROsjCIoLIgQEWWVHSNiXhITsyfX7YwYMMZAhmWSWXO/XK6/MOXOfc647J/nOyZkz9xFVxRhjTHhxBLoAY4wx/mfhbowxYcjC3RhjwpCFuzHGhCELd2OMCUMW7sYYE4ZqDHcRmSkiB0Rk3WmeFxF5RkS2isgaEenn/zKNMcacDV+O3F8Ghp/h+cuBFO/XROCFupdljDGmLmoMd1X9DDhyhiZXA6+qx1dAUxFp468CjTHGnD2XH9bRDthdaTrLO29v1YYiMhHP0T3R0dGp7du3r9UGKyoqcDjC4+0C60vwCZd+gPUlWNWlL5s3bz6kqi1rauePcJdq5lU7poGqTgOmAaSlpemKFStqtcHFixczePDgWi0bbKwvwSdc+gHWl2BVl76IyHe+tPPHy2AWUPkQPBnY44f1GmOMqSV/hPs84BbvVTMDgFxV/cEpGWOMMQ2nxtMyIvImMBhoISJZwCOAG0BVXwQ+AK4AtgIFwO31Vawxxhjf1BjuqjqmhucVuMtvFRljQkJpaSlZWVkUFRU1yPYSEhLYuHFjg2yrvvnSl6ioKJKTk3G73bXahj/eUDXGNEJZWVnExcXRqVMnRKq7rsK/8vLyiIuLq/ftNISa+qKqHD58mKysLDp37lyrbYTHdUXGmAZXVFREYmJigwR7YyMiJCYm1um/Igt3Y0ytWbDXn7r+bC3cjTEmDNk5d2NMyHI6nfTs2fPk9Ny5c+nUqVPgCgoiFu7GmJAVHR3N6tWrT/t8WVkZLlfjjDk7LWOMCSsvv/wyN9xwA1dddRVDhw4FYMqUKaSnp9OrVy8eeeSRk22feOIJunbtyqWXXsqYMWP429/+BsDgwYM5MTzKoUOHTv43UF5ezqRJk06ua+rUqcD3wwmMHDmSbt26MXbsWDxXiUNmZiYXXnghvXv3JiMjg7y8PIYNG3bKi9KgQYNYs2aNX38OjfMlzRjjV394bz0b9hzz6zq7t43nkasuOGObwsJC+vTpA0Dnzp2ZM2cOAEuXLmXNmjU0b96cBQsWsGXLFpYvX46qMmLECD777DOaNGnCrFmzWLVqFWVlZfTr14/U1NQzbm/GjBkkJCSQmZlJcXExgwYNOvkCsmrVKtavX0/btm0ZNGgQS5YsISMjg1GjRvHWW2+Rnp7OsWPHiI6O5pZbbuHll1/mqaeeYvPmzRQXF9OrVy8//NS+Z+FujAlZpzstc9lll9G8eXMAFixYwIIFC+jbty8A+fn5bNmyhby8PK699lpiYmIAGDFiRI3bW7BgAWvWrOHtt98GIDc3ly1bthAREUFGRgbJyckA9OnTh507d5KQkECbNm1IT08HID4+HoBrr72WQYMGMWXKFGbOnMltt91Wtx9ENSzcjTF1VtMRdkNr0qTJyceqyuTJk7nzzjtPafPUU0+d9nJDl8tFRUUFwCnXmqsqzz77LMOGDTul/eLFi4mMjDw57XQ6KSsrQ1Wr3UZMTAyXXXYZ7777LrNnz6a2I+SeiZ1zN8aEtWHDhjFz5kzy8/MByM7O5sCBA1x88cXMmTOHwsJC8vLyeO+9904u06lTJ1auXAlw8ij9xLpeeOEFSktLAdi8eTPHjx8/7ba7devGnj17yMzMBDyfTC0rKwNgwoQJ3HvvvaSnp5/8L8Of7MjdGBPWhg4dysaNGxk4cCAAsbGxvP766/Tr149Ro0bRp08fOnbsyI9+9KOTyzzwwAPceOONvPbaa1xyySUn50+YMIGdO3fSr18/VJWWLVsyd+7c0247IiKCt956i3vuuYfCwkKio6NZuHAhAKmpqcTHx3P77fU01qKqBuQrNTVVa+uTTz6p9bLBxvoSfMKlH6r125cNGzbU27qrc+zYsXpd/yOPPKJTpkyp122ccOzYMc3OztaUlBQtLy8/bbvqfsbACvUhY+20jDHGNLA33niD/v3788QTT9TbrQPttIwxxgCPPvpog23rpptu+sEbvP5mR+7GmFpTrfZ2ycYP6vqztXA3xtRKVFQUhw8ftoCvB+odzz0qKqrW67DTMsaYWklOTiYrK4uDBw82yPaKiorqFHbBxJe+nLgTU21ZuBtjasXtdtf6LkG1sXjx4pOfMg11DdEXOy1jjDFhyMLdGGPCkIW7McaEIQt3Y4wJQxbuxhgThizcjTEmDFm4G2NMGLJwN8aYMGThbowxYcjC3RhjwlDIhfuBY0V8mlVqgxUZY8wZhFy4v75sFy+tK2HCKys4mFcc6HKMMSYohVy43/eTFG7qFsHnWw8x/KnPWLhhf6BLMsaYoBNy4e5wCEM7ufnPPReRFB/FhFdX8Mi76ygqLQ90acYYEzRCLtxPSEmKY85dFzL+os68svQ7rnl+CVsP5Ae6LGOMCQohG+4AkS4nv7+yOy/dns6BvGJGPPcF767ODnRZxhgTcD6Fu4gMF5FNIrJVRB6s5vkOIvKJiKwSkTUicoX/Sz29IV1b8Z97L+KCtvH8ctZqfjdnLcVldprGGNN41RjuIuIEngcuB7oDY0Ske5VmDwGzVbUvMBr4u78LrUmbhGjevGMAP//xubyxbBc3vriU7JzChi7DGGOCgi9H7hnAVlXdrqolwCzg6iptFIj3Pk4A9vivRN+5nA4evLwbU8elsv3gca585nO+3HooEKUYY0xASU0fBhKRkcBwVZ3gnR4H9FfVuyu1aQMsAJoBTYBLVXVlNeuaCEwESEpKSp01a1atis7Pzyc2NvaMbfYdr+CZVUXsO66M6hrB0I4uRKRW26tPvvQlVIRLX8KlH2B9CVZ16cuQIUNWqmpajQ1V9YxfwA3A9ErT44Bnq7S5H/i19/FAYAPgONN6U1NTtbY++eQTn9rlFZXqHa9kasffvq/3v7Vai0rLar3N+uJrX0JBuPQlXPqhan0JVnXpC7BCa8htVfXptEwW0L7SdDI/PO0yHpjtfbFYCkQBLXxYd72KjXTx4s2p3HdpCv/+Ooux/1jGoXz7VKsxJvz5Eu6ZQIqIdBaRCDxvmM6r0mYX8BMAETkfT7gf9GehteVwCPdd2oXnb+rHuj25XP3cEjbtywt0WcYYU69qDHdVLQPuBuYDG/FcFbNeRB4TkRHeZr8G7hCRb4A3gdu8/z4EjZ/2asPsOwdSWl7ByBe+5PMtQfHaY4wx9cKn69xV9QNV7aKq56rqE955D6vqPO/jDao6SFV7q2ofVV1Qn0XXVq/kpsy9axDtmkVz+0uZvJW5K9AlGWNMvQjpT6jWRtum0fzr5wO58LwW/Pbfa3lq4WYbPtgYE3YaXbgDxEW5mXFrGiNTk3lq4RYmv7OWsvKKQJdljDF+4wp0AYHidjqYMrIXbROieObjrRzKL+a5m/oR5XYGujRjjKmzRnnkfoKIcP/Qrvzx6gtY9O0BbpmxnNzC0kCXZYwxddaow/2EcQM78czovqzafZRRU5faHZ6MMSHPwt3rqt5tmXFrOjsPH2fUVBt0zBgT2izcK7m4S0teG9+fg3nF3PjiUnYeOh7okowxplYs3KtI79ScNycOoLC0nBunLrW7OxljQpKFezV6tEtg1sQBVCiMnraUb/cdC3RJxhhzVizcT6NLUhyz7xyAy+Fg9LSvWL8nN9AlGWOMzyzcz+CclrG8decAYtxOxk5fxrpsC3hjTGiwcK9Bx8QmzJo40ALeGBNSLNx90CExhlkTB9IkwsnNM5axca+dgzfGBDcLdx91SIzhzYkDiHI5uXn6MrbstzHhjTHBy8L9LHRMbMIbd/TH6RDG/GMZ2w/aZZLGmOBk4X6WzmkZyxt39EdVGTt9GbuPFAS6JGOM+QEL91o4r1Ucr43vT0FJOWOnL2NfblGgSzLGmFNYuNdS97bxvPKzDI4cL+HmGcs4crwk0CUZY8xJFu510Kd9U6bfmsbuIwXcOnM5eUU2XLAxJjhYuNfRgHMSeeHmfmzce4wJr6ygqLQ80CUZY4yFuz9c0i2J/72xN8t3HuHuN1bZLfuMMQFn4e4nV/dpx6NXXcDCjfuZ/M5au+m2MSagGu09VOvDrRd24vDxEp5ZtIXE2EgevLxboEsyxjRSFu5+9qtLUzhyvJgXP91Gq7hIfnZR50CXZIxphCzc/UxE+MOIHhzKK+GP/9lAy7hIrurdNtBlGWMaGTvnXg+cDuGp0X1I79ic+2ev5suthwJdkjGmkbFwrydRbif/uCWNzi2acOdrK+1uTsaYBmXhXo8SYty8dHsGMZFObn8pk725hYEuyRjTSFi417N2TaN56bYM8orKuP2lTPsUqzGmQVi4N4DubeN54eZ+bD2Qzy/++TWl9iEnY0w9s3BvID9Kacmfru3J51sO8dCcdfYhJ2NMvbJLIRvQjent2X20gGc/3kqHxBjuGnJeoEsyxoQpC/cGdv9lXdh1pIAp8zfRMTGG2EAXZIwJS3ZapoGJCH+9vhdpHZtx/+xv2HrURpE0xvifT+EuIsNFZJOIbBWRB0/T5kYR2SAi60XkDf+WGV6i3E6m3ZJGm4Qonl5VZLfqM8b4XY3hLiJO4HngcqA7MEZEuldpkwJMBgap6gXAffVQa1hp3iSCmbelU14B41/J5JhdImmM8SNfjtwzgK2qul1VS4BZwNVV2twBPK+qRwFU9YB/ywxP57aM5e6+UWw/eNzGgTfG+JXUdEmeiIwEhqvqBO/0OKC/qt5dqc1cYDMwCHACj6rqh9WsayIwESApKSl11qxZtSo6Pz+f2NjweCsyPz+flUcjeWl9CZd2cHFz98hAl1Rr4bJfwqUfYH0JVnXpy5AhQ1aqalpN7Xy5WkaqmVf1FcEFpACDgWTgcxHpoao5pyykOg2YBpCWlqaDBw/2YfM/tHjxYmq7bLBZvHgxj1w5GOf7G5j+xQ4G9+vGzQM6BrqsWgmX/RIu/QDrS7BqiL74clomC2hfaToZ2FNNm3dVtVRVdwCb8IS98dHkK85nSNeWPDJvvY0iaYypM1/CPRNIEZHOIhIBjAbmVWkzFxgCICItgC7Adn8WGu6cDuGZMX05t2UTfv76SnYcOh7okowxIazGcFfVMuBuYD6wEZitqutF5DERGeFtNh84LCIbgE+ASap6uL6KDldxUW6m35KO0yGMfyWT3EK7gsYYUzs+Xeeuqh+oahdVPVdVn/DOe1hV53kfq6rer6rdVbWnqtbunVJDh8QYXrg5lV2HC7jnTbuCxhhTO/YJ1SA04JxE/nhNDz7bfJA//9+3gS7HGBOCbGyZIDUmowOb9uUx44sddGsdxw1p7WteyBhjvOzIPYg99NPzGXReIv89Zx0rvzsa6HKMMSHEwj2IuZwOnr+pH22aRnHnayvtNn3GGJ9ZuAe5pjERTL8ljaLScia+upKiUhtF0hhTMwv3EJCSFMdTo/qwbk8uv/33GruLkzGmRhbuIeLS7kk8MLQr767ew9TP7PNhxpgzs3APIb8YfC4/7dWGv374LYs32cCbxpjTs3APISLClJG96NY6nnveXMX2g/mBLskYE6Qs3ENMTISLaeNScTsd3PHqCrvJhzGmWhbuIah98xiev6kfOw8XcP9bq6mosDdYjTGnsnAPUQPPTeThK7uzcOMBnly4OdDlGGOCjA0/EMJuGdiR9XtyefbjrZzfJp4rerYJdEnGmCBhR+4hTET44zU96NuhKQ/86xu+3Xcs0CUZY4KEhXuIi3Q5efHmVGIjXdzx6gpyCkoCXZIxJghYuIeBpPgoXhyXyv7cYu5+w8aAN8ZYuIeNfh2a8fg1Pfhi6yH+YmPAG9Po2RuqYeTG9Pas25PL9C92cEG7eK7tmxzokowxAWJH7mHm91d2p3/n5jz477WsycoJdDnGmACxcA8zbqeDv4/tR4vYSO58bSUH84oDXZIxJgAs3MNQYmwkU8elcrSghF/8cyUlZfYGqzGNjYV7mOrRLoG/Xt+LzJ1HefS99YEuxxjTwOwN1TB2dZ92bNh7jKmfbueCtvGM7d8x0CUZYxqIHbmHud8M68aPu7TkkXfXs3zHkUCXY4xpIBbuYc7pEJ4Z05f2zWP4r9dXkp1jN9k2pjGwcG8EEqLd/OOWNErKKpj46goKS+wm28aEOwv3RuK8VrE8PaYPG/YeY9Lb39hNto0Jcxbujcgl3ZKYNKwr76/Zy98Xbwt0OcaYemTh3sj814/P5arebfnbgk0s2rg/0OUYY+qJhXsjIyL8z/W9uKBtPL+ctZot+/MCXZIxph5YuDdC0RFOpo1LI8rtZIKNAW9MWLJwb6TaNo1m6rhU9uYU8Yt/fk2pjQFvTFixcG/EUjs240/X9eTLbYf54/sbAl2OMcaPbPiBRm5kajKb9+cx7bPtpCTFMW6ADVFgTDiwI3fDb4d3Y0jXljw6bz1fbj0U6HKMMX7gU7iLyHAR2SQiW0XkwTO0GykiKiJp/ivR1LcTQxSc06IJ//XPr9lx6HigSzLG1FGN4S4iTuB54HKgOzBGRLpX0y4OuBdY5u8iTf2Li3Iz49Z0nA5h/MuZ5BaUBrokY0wd+HLkngFsVdXtqloCzAKurqbdH4H/AYr8WJ9pQB0SY3jx5lR2Hy3grjfsChpjQpnUNMaIiIwEhqvqBO/0OKC/qt5dqU1f4CFVvV5EFgMPqOqKatY1EZgIkJSUlDpr1qxaFZ2fn09sbGytlg02wdiXz7NKmbGuhMHtXdzaPQIR8Wm5YOxLbYRLP8D6Eqzq0pchQ4asVNUaT337crVMdX/ZJ18RRMQBPAncVtOKVHUaMA0gLS1NBw8e7MPmf2jx4sXUdtlgE4x9GQxEfPgtLyzexo96d2H8RZ19Wi4Y+1Ib4dIPsL4Eq4boiy+nZbKA9pWmk4E9labjgB7AYhHZCQwA5tmbqqFt0tCuDL+gNY//ZwMLN9gYNMaEGl/CPRNIEZHOIhIBjAbmnXhSVXNVtYWqdlLVTsBXwIjqTsuY0OFwCE+O6kOPtgncO2sV67JzA12SMeYs1BjuqloG3A3MBzYCs1V1vYg8JiIj6rtAEzjREU5m3JpG02g341/JZG+u3cXJmFDh03XuqvqBqnZR1XNV9QnvvIdVdV41bQfbUXv4aBUfxczb0zleXM7tL2WSV2SXSBoTCuwTqqZG3VrH8/zYfmw5kM9db6yySySNCQEW7sYnP+7Skieu6cFnmw/y8Lvr7DZ9xgQ5GzjM+Gx0Rgd2Hy3g+U+2kdwshruGnBfokowxp2Hhbs7Kry/rStbRQqbM30SbhCiu65cc6JKMMdWwcDdnxeEQpozszcG8Yn7z9hpaxUVxUUqLQJdljKnCzrmbsxbhcvDiuFTOaxXLz19fyfo9tb8Gfu6qbAb95WM6P/gfBv3lY+auyvZjpaa+2f4LXhbuplbio9y8fHsG8VEubnspk91HCs56HXNXZTP5nbVk5xSiQHZOIZPfWWsBESJs/wU3C3dTa60Tonh1fAYlZRXcMnM5x0rO7gqaKfM3UVhafsq8wtJypszf5M8yTT2x/RfcLNxNnZzXKo6Zt6WxN7eQJ1cUkV9c5vOye3Kq/8Tr6eab4GL7L7hZuJs6S+3YnL+P7cd3eRVMfHUFxWXlNS8EtG0afVbzTXCx/RfcLNyNX1zSLYnxPSL4ctth7pu1mvKKmk/RTBrWlWi385R50W4nk4Z1ra8yjR/Z/gtuFu7Gbwa1c/P7K7vzf+v2MfmdNTV+ivWavu3483U9adc0GgHaNY3mz9f15Jq+7RqmYFMntv+Cm13nbvxq/EWdyS0s5ZlFW4iPcvPfPz3/jHdyuqZvOwuDEGb7L3hZuBu/+9WlKRwrLGX6FzuIi3Lzy0tTAl2SMY2OhbvxOxHh4Su7k19cxpMLNxMT4eSOi88JdFnGNCoW7qZeOBzCX6/vRWFpOU98sJEot4NxAzsFuixjGg0Ld1NvnA7hyRv7UFRSzu/fXU+Ey8Go9A6BLsuYRsGuljH1KsLl4Pmx/bi4S0sefGct/16ZFeiSjGkULNxNvYtyO5k2LpULz01k0tvf8O5qG3vEmPpm4W4aRJTbyfRb0sno3JxfvbXaBpcypp5ZuJsGEx3hZOZt6fTvnMj9s1czZ5WdojGmvli4mwYVE+Fi5m3pDDgnkftnf8PsFbsDXZIxYcnC3TS46AgnM25N56LzWvCbt9fw+lffBbokY8KOhbsJiOgIJ/+4JY2fdGvFQ3PXMeOLHYEuyZiwYuFuAibK7eSFm1O5vEdr/vj+Bp5euKXGwcaMMb6xcDcBFeFy8OyYvlzfL5knF27mif9stIA3xg/sE6om4FxOB1NG9iIuysX0L3aQW1jKn6/rictpxx7G1JaFuwkKDofwyFXdSYh28/SiLRwtKOW5m/oSVeVmEMYY39ihkQkaIsKvLuvCH0ZcwKJv9zNuxjJyCkoCXZYxIcnC3QSdWy/sxDOj+/LN7lxGvriUrKMFgS7JmJBj4W6C0lW92/LKzzLYf6yI6/7+JeuycwNdkjEhxcLdBK2B5yby9s8vxOkQbpy6lI+/3R/okowJGRbuJqh1bR3H3LsG0blFEya8soJXvtwZ6JKMCQkW7iboJcVHMfvOgVzSrRWPzFvPw++uo6y8ItBlGRPUfAp3ERkuIptEZKuIPFjN8/eLyAYRWSMii0Sko/9LNY1Zk0gXU8elMfHic3h16Xfc9lImuQWlgS7LmKBVY7iLiBN4Hrgc6A6MEZHuVZqtAtJUtRfwNvA//i7UGKdD+N0V5zNlZC+W7TjMiOe/YPP+vECXZUxQ8uXIPQPYqqrbVbUEmAVcXbmBqn6iqieuV/sKSPZvmcZ874a09syaOICCknKufX4JH67bG+iSjAk6UtM4HiIyEhiuqhO80+OA/qp692naPwfsU9XHq3luIjARICkpKXXWrFm1Kjo/P5/Y2NhaLRtsrC+1d7SogmdXFbM9t4IrOru5PsWN0yF1Xq/tk+BkffEYMmTISlVNq7Ghqp7xC7gBmF5pehzw7Gna3oznyD2ypvWmpqZqbX3yySe1XjbYWF/qpqi0TCe/s0Y7/vZ9HT11qR7MK6rzOm2fBCfriwewQmvIV1X16bRMFtC+0nQysKdqIxG5FPhvYISqFvuwXmPqLNLl5E/X9uRvN/Tm611HueLpz1m67XCgyzIm4HwJ90wgRUQ6i0gEMBqYV7mBiPQFpuIJ9gP+L9OYMxuZmszcuwYRG+li7PSveGbRFsorbOhg03jVGO6qWgbcDcwHNgKzVXW9iDwmIiO8zaYAscC/RGS1iMw7zeqMqTfnt4ln3j0XMaJ3W/7fR5sZO/0r9uYWBrosYwLCpyF/VfUD4IMq8x6u9PhSP9dlTK0s3LCf5TuOALBs+xF+8r+fMjq9PfPX72dPTiFtm0YzaVhXrunbzu/bnrsqmynzN9X7dnzx0Ny1vLlsN/f1KGX85A8Y0789j1/TMyC1mMCw8dxN2Ji7KpvJ76ylsLQcAAUKS8qZuWTnyTbZOYVMfmctgF+Dt+q262s7vnho7lpe/2rXyely1ZPTFvCNhw0/YMLGlPmbTobrCdWddS8sLWfK/E31vu362I4v3ly2+6zmm/Bk4W7Cxp4c38+vZ59F27ps+2xq8pfy03x25XTzTXiycDdho23TaJ/bOh3CF1sO1fu2z6Ymf3FK9R/kOt18E54s3E3YmDSsK9FV7rnqdghu56mhFuF00LxJBDfPWMYD//qGI8frfiu/6rYd7XYyaVjXOq/7bI3p3/6s5pvwZG+omrBx4o3LqlesVDdveI/WPLNoC9M+286ijfv53RXnMzI1Ganl0e3pth2Iq2VOvGl64hy7U8SulmmELNxNWLmmb7tqA7W6eb8Z3o0Rfdry0Jx1THp7DbNX7OYPI3r4fduB8Pg1PXn8mp4sXryYbWMHB7ocEwB2WsY0at1axzP7zoH89fqebDt4nCuf/ZzXNhSTU1D3UzXGBJKFu2n0HA5hVHoHPvn1YG4e0JGPd5Xx4ymLeWnJDkrtjk8mRFm4G+OVEOPmsat78NigaHq0i+cP721g2JOf8eG6fSdGPTUmZFi4G1NF+zgHr4/vz/Rb0nA4hJ+/vpKRLy5l2XYbbdKEDgt3Y6ohIlzaPYkPf/kj/nxdT3YfKWDUtK+4ZeZy1mTlBLo8Y2pk4W7MGbicDsZkdODTSUP43RXdWJOVw4jnljD+5UwLeRPULNyN8UF0hJOJF5/L578ZwgNDu7Diu6OMeG4Jt720nMydRwJdnjE/YOFuzFmIi3Jz9yUpfPHbIUwa1pW1Wbnc8OJSbnxxKQs37KfCbhBigoSFuzG1EBfl5q4h5/HFby/h4Su7k51TyIRXV3DZk5/yxrJdFJaU17wSY+qRhbsxdRAd4eRnF3Vm8aTBPD26D1FuJ7+bs5aBf1nEX/7vW3YfKQh0iaaRsuEHjPEDt9PB1X3aMaJ3WzJ3HuWlJTuY9tk2pn62jUu6tmLsgA78uEsrnA4bmdE0DAt3Y/xIRMjo3JyMzs3Zk1PIm8t38eby3Sx6eQVtE6K4Ia09I1OTad88JtClmjBn4W5MPWnbNJpfD+3KPZeksGjjft5YvotnPt7C04u2MPCcRK5PTWZ4j9bERtqfofE/+60ypp5FuBxc3rMNl/dsQ9bRAt75Opu3V2bxwL++4aG5a7mse2tG9G7LxV1aEOly1rxCY3xg4W5MA0puFsO9P0nhnkvO4+tdR5mzKpv31+zlvW/2EBflYmj31lzeozUXpbQgym1Bb2rPwt2YABARUjs2J7Vjcx656gKWbD3Ee9/s5aMN+/j311nERrr4cdeWDO2exOCurUiIdge6ZBNiLNyNCTC308Hgrq0Y3LUVJWU9+XLbIT5ct4+FGw/wnzV7cTqE9E7NuKSbp01Kq9ha3zHKNB4W7sYEkQjX90FfUaGs2p3Dx9/uZ9HGA/zpg2/50wff0iYhih+ltOCilJZceG4iLWIjA122CUIW7sbIwl0aAAANK0lEQVQEKYdDSO3YjNSOzZg0rBt7cgr5bPNBPt18kA/X7WP2iiwAurWOY8A5iQw4J5H0Ts1ItLA3WLgbEzLaNo1mdEYHRmd0oLxCWZudy5Kth1i67TCzMnfx8pc7ATivVSxp3heF8uMVqKqdxmmELNyNCUFOh9CnfVP6tG/KXUPOo7isnLVZuSzfeYTMHUf4YO1eZmXuBuAvKz+id3JTerdvSu/kBHomJ9AqLirAPTD1zcLdmDAQ6XKS1qk5aZ2aw2CoqFC2HcznjQVfURiTxOrdOTz38RZODFrZKi6SHu0S6N4mnu5t4zm/TTwdmsfY8AhhxMLdmDDkcAgpSXH8uL2bwYN7AVBQUsb6PcdYk5XL+j25rMvO5dPNByn3Jn6U20FKqzi6JMXRJSmWlKRYzmsZR7tm0Rb6IcjC3ZhGIibCRXqn5qR3an5yXlFpOVsP5LNh7zE27ctj0748Pt9ykH9/nXWyTaTLQecWTejcogmdWjShc2ITOibG0KlFE1rGRuKw4A9KFu7GNGJRbic92iXQo13CKfNzCkrYeiCfbQfz2Xognx2HjrNpfx4fbdhPWaUbkkS6HLRvHkNys2jaN4uhXbNo2jWNPvm9RWykHfUHiIW7MeYHmsZEfH8Ov5Ky8gr25BSx4/Bxdh0pYPeRAr47fJyso4Ws2pVDbmHpKe1dDiEpPorWCVG0jo8iKT6KpPhIWsVH0iouilZxkbSIjaRpjNuu6PEzC3djjM9cTgcdEmPokFj9kMV5RaVk5xSSfbSQPblF7M0pZF9uEfuOFbFx7zE+2XSAgmruUuV2ColNIkmMjSAxNpLEJhE0r/TVLCaC746U02ZfHs1i3MRHu23snRpYuBtj/CYuyk231m66tY4/bZv84jL2HyviwLFiDuQVcSi/hEP5xRzMK+bIcc/j7QfzOXK85AcvBH9e/tnJx1FuBwnRbuKj3J7v0W7io1zERbmJ836PjXIRF+miSaSLWO9Xk0gnsZEuYiJdxLidYfuegU/hLiLDgacBJzBdVf9S5flI4FUgFTgMjFLVnf4t1ZjwNXdVNlPmb2JPTiFtm0YzaVhX/rViF0u2HTnZZtC5zbkhrcMP2gE/mLfiuyO8uWw39/UoZfzkDxjTvz2PX9PTp+1Wt75r+rbzue4T2y5XxSnyg23HRrqIbRnL2qzcGrf96FUp/KhLC44cL+HTpSvokHI+RwtK+WrbYT7dfJD9x4q9p4JiKCwtZ+uBMo4VlZJXVHbyKqCaRLudxEQ4iY448d1FtNtBTISLaLeTSLeDaLeTKLeTKLeDKNf3jyNdnucjXd7HLgeRbgcRTicRLsf3X87vv7udgmr930i9xnAXESfwPHAZkAVkisg8Vd1Qqdl44Kiqnicio4G/AqPqo2Bjws3cVdlMfmcthaWeo9TsnELue2v1D9ot2XbklLDPzilk0tvfgEKpN8iycwq5/63VVFRarlyV17/aBXBKyFa33Un/+gYESsu/X9/kd9YC/CDgq1ve39t+ZN56/nxdT67p246DiU4G92rL3FXZfPztgZPLFpVWkHW08GQ7AFWlqLSCvOJS8ovKyC/2fhWVUVBSzvGSMo4Xl3G8uJzjxWUUlJZTWFJOQUkZhaUVFJaUcTCvmELv/KLScgpLPd99fM04o3HdIxhS99WckS9H7hnAVlXdDiAis4CrgcrhfjXwqPfx28BzIiLaEC9PxoS4KfM3nQyqs3UiCCurqKYdwJvLdp8SsNVtt7Sa5CosLWfK/E0/CPfqlm+IbVe3bNV2IkK092i8VdxpiqoFVaW0XCkqK6e4tIKi0nJKyisoLq2guKyckrIKissqKCmr8MwvK6e0TCku98wrLa+gtKyCuPxd/ivqNKSm/BWRkcBwVZ3gnR4H9FfVuyu1Wedtk+Wd3uZtc6jKuiYCE72TXYFNtay7BXCoxlahwfoSfBq0HxGtz0utr3WXF+TijPn+MseSfVtX1na7lZet6/K1XLYFcOhMy1atMYjV5Xeso6q2rKmRL0fu1b3bUPUVwZc2qOo0YJoP2zxzQSIrVDWtrusJBtaX4BMu/QBPX8pyD4RNX8Jpv9R3Xxw+tMkC2leaTgb2nK6NiLiABOAIxhhjAsKXcM8EUkSks4hEAKOBeVXazANu9T4eCXxs59uNMSZwajwto6plInI3MB/PpZAzVXW9iDwGrFDVecAM4DUR2YrniH10fRaNH07tBBHrS/AJl36A9SVY1XtfanxD1RhjTOjx5bSMMcaYEGPhbowxYSjow11EokRkuYh8IyLrReQP3vmdRWSZiGwRkbe8b/YGPRFxisgqEXnfOx2q/dgpImtFZLWIrPDOay4iH3n78pGINAt0nb4QkaYi8raIfCsiG0VkYCj2RUS6evfHia9jInJfiPblV96/93Ui8qY3B0L1b+WX3n6sF5H7vPPqfZ8EfbgDxcAlqtob6AMMF5EBeIY4eFJVU4CjeIZACAW/BDZWmg7VfgAMUdU+la7XfRBY5O3LIu90KHga+FBVuwG98eyfkOuLqm7y7o8+eMZ5KgDmEGJ9EZF2wL1Amqr2wHMhx4lhTULqb0VEegB34Pmkf2/gShFJoSH2iaqGzBcQA3wN9Mfz6S6Xd/5AYH6g6/Oh/mTvjrwEeB/Ph79Crh/eWncCLarM2wS08T5uA2wKdJ0+9CMe2IH34oJQ7kuV+ocCS0KxL0A7YDfQHM8Vfe8Dw0LxbwW4Ac9giyemfw/8piH2SSgcuZ84lbEaOAB8BGwDclS1zNskC88vRLB7Cs+OPTEERyKh2Q/wfAJ5gYis9A4rAZCkqnsBvN9bBaw6350DHARe8p4umy4iTQjNvlQ2GnjT+zik+qKq2cDfgF3AXiAXWElo/q2sAy4WkUQRiQGuwPOBz3rfJyER7qparp5/NZPx/HtzfnXNGraqsyMiVwIHVLXy2Bc+DdsQpAapaj/gcuAuEbk40AXVkgvoB7ygqn2B4wT5aYuaeM9FjwD+FehaasN7/vlqoDPQFmiC5/esqqD/W1HVjXhOJ30EfAh8A5SdcSE/CYlwP0FVc4DFwACgqXeoA6h+SIRgMwgYISI7gVl4Ts08Rej1AwBV3eP9fgDPed0MYL+ItAHwfj8QuAp9lgVkqeoy7/TbeMI+FPtywuXA16q63zsdan25FNihqgdVtRR4B7iQ0P1bmaGq/VT1Yjwf8txCA+yToA93EWkpIk29j6Px7PiNwCd4hjoAz9AH7wamQt+o6mRVTVbVTnj+Zf5YVccSYv0AEJEmIhJ34jGe87vrOHUYipDoi6ruA3aLSFfvrJ/gGc465PpSyRi+PyUDodeXXcAAEYkREeH7fRJyfysAItLK+70DcB2efVPv+yToP6EqIr2AV/C8Y+4AZqvqYyJyDp4j4ObAKuBmVS0OXKW+E5HBwAOqemUo9sNb8xzvpAt4Q1WfEJFEYDbQAc8f6A2qGvQDyIlIH2A6EAFsB27H+7tG6PUlBs+bkeeoaq53XsjtF+8lz6PwnMJYBUzAc449pP5WAETkczzvr5UC96vqoobYJ0Ef7sYYY85e0J+WMcYYc/Ys3I0xJgxZuBtjTBiycDfGmDBk4W6MMWHIlxtkG9OgvJeJLfJOtgbK8QwRAJChqiUBKewMRORnwAfe6+aNCTi7FNIENRF5FMhX1b8FQS1OVS0/zXNfAHer6uqzWJ+r0lgpxviVnZYxIUVEbhXP+P6rReTvIuIQEZeI5IjIFBH5WkTmi0h/EflURLaLyBXeZSeIyBzv85tE5CEf1/u4iCwHMkTkDyKS6R2f+0XxGIVnOOq3vMtHiEhWpU9WDxCRhd7Hj4vIVBH5CM9gZS4R+X/eba8RkQkN/1M14cjC3YQM79jY1wIXegeSc/H9zdgTgAXewcxKgEfxfGz9BuCxSqvJ8C7TD7hJRPr4sN6vVTVDVZcCT6tqOtDT+9xwVX0LWA2MUs946jWdNuoLXKWq44CJeAaUywDS8QzC1qE2Px9jKrNz7iaUXIonAFd4hhwhGs9H7QEKVfUj7+O1QK6qlonIWqBTpXXMV9WjACIyF7gIz9/B6dZbwvdDLQD8REQmAVFACzxD0f7fWfbjXVUt8j4eCpwvIpVfTFLwfCTdmFqzcDehRICZqvr7U2Z6RgqsfLRcgecOXiceV/49r/omk9aw3kL1vjHlHbflOaCfqmaLyON4Qr46ZXz/n3HVNser9OkXqroIY/zITsuYULIQuFFEWoDnqppanMIYKp57psbgGTN8yVmsNxrPi8Uh76iY11d6Lg+IqzS9E8+t7qjSrqr5wC9ODGUrnvugRp9ln4z5ATtyNyFDVdd6RwtcKCIOPKPs/ZyzG9f7C+AN4FzgtRNXt/iyXlU9LCKv4Bne+DtgWaWnXwKmi0ghnvP6jwL/EJF9wPIz1DMVz8iAq72nhA7gedExpk7sUkjTaHivROmhqvcFuhZj6pudljHGmDBkR+7GGBOG7MjdGGPCkIW7McaEIQt3Y4wJQxbuxhgThizcjTEmDP1/pGenSMj5bcYAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "%matplotlib inline\n",
+    "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n",
+    "data_pred['Frequency'] = logmodel.predict(data_pred)\n",
+    "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n",
+    "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n",
+    "plt.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "hideCode": false,
+    "hidePrompt": false,
+    "scrolled": true
+   },
+   "source": [
+    "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Computing and plotting uncertainty"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9+P/XObNkmewhC0vYA4RFXBAQV6LIvrhVkSoWUWsv9ddqq3bzWnurtbetxfZ+Vdy3LooLQhSXoFAVURSJ7GsgAbJvk8w+5/z+mGRISIBJyGSWvJ+Ppsk5c+bk8zHMvOezvT+Krus6QgghxAnUUBdACCFEeJIAIYQQokMSIIQQQnRIAoQQQogOSYAQQgjRIQkQQgghOhS0APGLX/yCCy64gDlz5nT4uK7r/M///A/Tpk1j7ty5bN++PVhFEUII0QVBCxBXX301zzzzzEkf37BhA8XFxXzwwQf87ne/48EHHwxWUYQQQnRB0ALE+eefT3Jy8kkfLywsZMGCBSiKwtlnn01DQwMVFRXBKo4QQohOMobqF5eXl5Odne0/zs7Opry8nMzMzFM+z2pzoQAoCoDv5+OHzT+3ekzxHbe9TkFVfAfN3/zPEUII4ROyANFRho9A3qQbbW6qqhu7vTz+YNJcDkUBVVH8P7c+pzYfq4qCqoKqKv5rz1RGRiKVldYzvk+4iub6RXPdQOoX6TIyEjv9nJAFiOzsbMrKyvzHZWVlp209BJPe/H++73rrswFTFDAoii9gNH8ZmoOHQVUwGBQMqkwcE0JEhpAFiPz8fF555RVmz57N1q1bSUxMDGmA6A66Dh5dB+3kgUWB5mChYlAVjAYVg0HBZFBRVenmEkKEj6AFiLvvvpsvv/yS2tpaLrnkEn784x/j8XgAWLhwIZdeeinr169n2rRpxMXF8fDDDwerKGFFBzyajkfztntMVUAxGWmwuTAZVMwmVVocQoiQUSIt3fexqqagjEGEi7Q0CzU1Tf5jVVWIMaqYTQZiTIaIb2VEcz9vNNcNpH6RLqLGIERgNE3H7vJid/laHGajSqzZQKzZGPHBQggR3iRARBiXR8Pl0bDa3JhNBuJjjMSYDaEulhAiCkmAiFA64HR7cbq9GFSF+FgjcTFGVFnPIYToJjICGgW8mo7V5qaqzo7N4Ql1cYQQUUICRBTRdGiwuaiqs+N0t58lJYQQnSEBIgp5NJ1aq5P6JhdaZE1SE0KEEQkQUczu9FBV75DWhBCiSyRARDmtuTVhtblCXRQhRISRANFLNDk81DQ40E6RBkQIIVqTANGLuDwaVQ0O3B7pchJCnJ4EiF5G03RqGpzYnTIdVghxahIgeiEdqG9y0dDk6nBfDiGEAAkQvZrN6aHW6pSpsEKIDkmA6OVcHo2aegdeTQt1UYQQYUYChMCj6VQ3OHF7JEgIIY6TACGA5sFrq8xwEkIcJwFC+Ok61FqdeLzSkhBCSIAQJ9B0qJOBayEEEiBEBzyaTp3VKVNghejlJECIDrXsWieE6L0kQIiTsjk9OFyy4lqI3koChDilhiaXrJEQopeSACFOSdOhvlFShQvRG0mAEKfl8mg02mU8QojeRgKECEiT3S3rI4ToZSRAiIDoILOahOhlJECIgDndXpwuScUhRG8hAUJ0itUme0gI0VtIgBCd4tF0bLIbnRC9QsQFiD+89BV7SupCXYxerdHuRtOkFSFEtDOGugCddeBIPQeO1HP28D7MumAQCXGmUBep19F1aHK4SYw3h7ooQoggirgWhNL8/dt9Vfz1ta18u7dK+sRDwOb0SMZXIaJcxAWI+24+n+y0eMD3JvXax/t45YM9WG2y2rcn6TrYHDIWIUQ0i7gAMbR/Mv919VimTcjBoPraEzsP1fLX17fy7T5pTfQkm8MtrQgholjEBQgAg6oy9dz+LLt6HP0zLADYnV5eW7ePfxXulU+2PUTTwS4zmoSIWkENEBs2bGD69OlMmzaNFStWtHv86NGj3HTTTSxYsIC5c+eyfv36Tt0/Ky2eH84fy/SJx1sT3x2o4fGVW9lbKjOdekKTwyOtNiGiVNAChNfr5aGHHuKZZ56hoKCANWvWsG/fvjbXPPHEE8ycOZO3336bxx57jN/+9red/j0GVeHSs/vzo6vG+scmGmxunn93FwUbiyV/UJBpmo7dKaurhYhGQQsQRUVFDBo0iJycHMxmM7Nnz6awsLDNNYqi0NjYCIDVaiUzM7PLv69vuoUfXTWWi8/q65/p9Nl3ZTy5ajtVdfYu31ecns0hOZqEiEZBWwdRXl5Odna2/zgrK4uioqI21yxbtoxbb72VV155BbvdzvPPPx/QvdPSLCd9bNGs0UwYk81zq3dQ3+jkaFUT//fWNhZOH8nksX27Vpkedqr6havE5FhizYH9c8rISAxyaUInmusGUr/eJmgBoqN+aUVR2hwXFBRw1VVXsWTJErZs2cK9997LmjVrUNVTN2xqappO+XhGYgzLrh7Lm+sPsPNQLU63lxfW7GDb3krmXjgEkzF8x+bT0iynrV84arI6SE2MOe11GRmJVFZae6BEPS+a6wZSv0jXleAXtHfK7OxsysrK/Mfl5eXtupBWrlzJzJkzATjnnHNwOp3U1tZ2y++3xJr4/pUjmDtlsH8Ae/PuSp54e5t0OQWB0+2VrUmFiDJBCxDjxo2juLiYkpISXC4XBQUF5Ofnt7mmb9++bNy4EYD9+/fjdDpJS0vrtjIoisIFY7P54fwx/k+3ZTU2/v7Wd2w7WNNtv0f4yPRiIaJL0AKE0WjkgQceYOnSpcyaNYuZM2eSm5vL8uXL/YPV999/P6+99hrz5s3j7rvv5g9/+EO7bqju0D8jgWVXj2P04FQAXG6Nf3y4h/e/PCxJ57qR3SlTXoWIJooeYa/oY1VNVFU3dum5uq7z6XfHeH/TYVriwvD+yVx/+XAsseGR9C9SxyBaJFvMxMWcfGgrmvt5o7luIPWLdGE1BhGOFEXh4rP68YPZeVhifW9i+47U8//e2kZZjS3EpYsOsrJaiOjRqwJEi2H9kvmvq8cxoDlNR63VyZNvb2NHsYxLnCmXR5PFiUJEiV4ZIABSEmK4be4YzsntA/je2F75YA/rvimVfvQz5JB9q4WICr02QACYjCrXXjaMmZMH0jI2/tHmUl77eB9uj3wK7irpZhIiOvTqAAHHxyUWzxhFjMkAwNZ91TxbsINGu6SQ6AqvpuNySytCiEjX6wNEixE5KfxwwfH1EofLG3ni7W1U1Mqiuq6wSzeTEBFPAkQrWanx3LlgLAOzEoDmwetV2zh4rCHEJYs8TpesiRAi0kmAOEFCnIlbZ4/mrGHpgG/A9bmCnWzdVxXikkUWTfel3xBCRC4JEB0wGVW+lz+cS8/uB/j61P+9bh//2Xo0xCWLLLJPhBCRTQLESaiKwvSJA5l/0RD/DKf3Nh3m3S8OyT7MAXK5vZLKRIgIJgHiNCaNzuL7V47EZPD9p/q06BgrP94vmUsDoAMOl0x5FSJSSYAIQN6gVJbMziMuxjcN9tt9Vby0drdM5QyATdZECBGxJEAEaFB2IrfPG0OyxQzA3tJ6nnt3pywKOw2PV8ftkUAqRCSSANEJWanx3DF/DBkpsYBvrcTTq3dgtblCXLLwJvtECBGZJEB0UkpCDLfPG0P/Pr5Ef2U1Nla8s4NaqzPEJQtfDpdXBvaFiEASILrAEmvi1jl5DOnry69e3eBgxTvbqaqXVdcd0QGHdMUJEXEkQHRRrNnILTPzGDUwBYD6JhdPv7ODctlXokMyWC1E5JEAcQZMRpUbp41g7FDfPtpWu5unV+/gaFXk7ggXLDJYLUTkkQBxhowGlevzc/37SticHp5Zs4PSiq5tixrNbLKyWoiIIgGiGxhUhWsuG8bEvEzANyj7bMFODpdH7/62XSEJ/ISILBIguomqKMy/aAgXjMkGfInqnn93F4fKJEi0kAR+QkQWCRDdSFEU5kwZxIVjWwWJ93ZKkGhFFhYKETkkQHQzRVGYdcEgLj6rLwAutyZBohWHjEMIETEkQASBoijMmDSQS8YfDxIvvLdLxiQATdelm0mICCEBIkiU5nThLS2JljGJkgoJEg7ZjlSIiCABIohaWhIXjWsbJI708nUSTkkBLkREkAARZIqiMHPyQKY0D1y3bGFa1otXXGs6kipdiAggAaIHKIrC7AsGMWl0FuCbyfPsmh1U1Pbe3E3SzSRE+JMA0UMURWHuhYM5b2QGAE0OD88W7KC63hHikoWGQ1oQQoQ9CRA9SFUUrrp4KGcP96XlsNrcPFuwg7rG3pcqXNN06WYSIsxJgOhhanNajjFDfAn+6hpdPFuws1duOiStCCHCmwSIEDCoCtfnD2dEji9VeHW9g+cKdmJzuENcsp4l4xBChDcJECFiNKgsmjaCIX2TACivtfPCe7t61cY60s0kRHiTABFCJqPKzdNHkpOZAEBpZRNPvFmE26OFuGQ9R7qZhAhfAQWIW2+9lY8//rjTqZo3bNjA9OnTmTZtGitWrOjwmnfffZdZs2Yxe/Zs7rnnnk7dPxrEmA0snjGKrNQ4AHYfquVfhXvxar0jSEg3kxDhK6AAcf311/Piiy9yxRVXsGLFCmpra0/7HK/Xy0MPPcQzzzxDQUEBa9asYd++fW2uKS4uZsWKFfzzn/+koKCAX/7yl12rRYSLjzXyg9l5pCXFALDzUC1vrj+A1gv2TtA02WlOiHAVUIC48soreeGFF3j66aepqKhgzpw53HvvvWzbtu2kzykqKmLQoEHk5ORgNpuZPXs2hYWFba557bXXWLRoEcnJyQCkp6efQVUiW1K8mVtn55Gc4AsSW/ZW8e7GQ71igx27tCKECEvGrjzJZDIRExPDfffdx8UXX8z999/f7pry8nKys7P9x1lZWRQVFbW5pri4GIAbbrgBTdNYtmwZl1xyyWl/f1qapSvFDntpaRb+vxvO5s+vfE2Tw8Pn28rokxbPrClDQl20bnXi38+gKmSkR8ffNCMjMdRFCCqpX+8SUID44IMPeOWVV6iurubGG2+koKAAi8WCx+Phyiuv7DBAdPTJV1GUNsder5dDhw7x8ssvU1ZWxqJFi1izZg1JSUmnLE9NTfQmu+vXJ4GbZ4zk2TU7cXk03tlwADSNyaOzT//kCJCWZunw76e53JiMhhCUqPtkZCRSWRm92XqlfpGtK8EvoACxcuVKbrvtNi6++OK2TzYa+fWvf93hc7KzsykrK/Mfl5eXk5mZ2eaarKwszj77bEwmEzk5OQwZMoTi4mLOOuusztYjquRkJrLoyhG8tHY3Xk1n9afFxMeYOGtY9HbB2V3eiA8QQkSbgMYgnnrqqXbBoUV+fn6H58eNG0dxcTElJSW4XC4KCgraXXvFFVewadMmAGpqaiguLiYnJ6cz5Y9auQNS+F7+cBRAB17/eB/7jtSHulhBI7OZhAg/AQWIG2+8kfr6429OdXV1LFq06JTPMRqNPPDAAyxdupRZs2Yxc+ZMcnNzWb58uX+w+uKLLyYlJYVZs2axePFi7r33XlJTU8+gOtFl3NB05l40GACvpvPKB7s5UtkY2kIFiabJTnNChBtFD2CazPz581m1atVpz/WEY1VNVFVH55skdNxHX/h1KYVflwJgiTVyx/wx9EmOC0XxztjJxiAA4mKMJFvMPVyi7tMb+rClfpGrK2MQAbUgNE3DZju+wU1TUxNer3za6yn55/b37yXR5PDw/Lu7ojK5n9Pl6RXTeoWIFAEFiDlz5rBkyRJWrVrFqlWruPXWW5k3b16wyyaaKYrC3CmDGducAbbW6vTlbYqyrTs1HelmEiKMBDSL6Y477iAzM5N169ah6zo33HADCxYsCHbZRCuqqnDd1OE0OXZx8FgDx6ptvPLBHm6ZOQqjIXpSajlcXmLNXVqeI4ToZgGNQYST3jgG0ZrD5eHp1Ts4Vu3r8hs3NI3rL89FPWGNSbg6Xf0UICM1LmLq01pv6MOW+kWuoK2DqK6u5uWXX6akpASP53i3xvLlyzv9C8WZiTUbWTxzFE++vY26RhffHaghMf4Qsy8Y1G4hYiTSAafLS1yMtCKECLWAXoU//vGPGTZsGBdccAEGgyxmCrWkeDNLZuXx5Krt2Jy+lBxJFjOXjO8X6qJ1C4cECCHCQkCvwoaGBn73u98FuyyiE/qkxLF45kieWb0Tt1dj7abDJMabOCc3I9RFO2NOtxevpmFQo2dsRYhIFNArMDc3l/Ly8mCXRXRSTmYiC6/IRW3uWXrjkwPsLa0LbaG6id0ps5mECLWAWxDz5s3jnHPOISYmxn9exiBCb9SgVBZcPJQ3N/j2j3j1wz3cPncM/fpEdnZUm9NDQpwp1MUQolcLKEDMmTOHOXPmBLssoosmjMqkwebio82luNwaL763izvmjyEtKTbUResyTdNxurzEmGXMS4hQCShAXHXVVcEuhzhDU8/pT32ji692VWC1u3mhOUhYYiP3U7jN6ZEAIUQIBTQGUVxczMKFC/3ZWLdv387f/va3oBZMdI6iKMy7aAh5g3zJDqvqHby0djeuCN7Os2WwWggRGgEFiAcffJA777yTxETfQou8vDzWrl0b1IKJzjOoCtdfPpyBWQkAlFQ08u/CfXi1iFoL2YYMVgsROgEFCKvVyiWXXOJfiKWqKiZT5HZdRDOz0cBN00eSnuwbf9h5qJbVnx2M2CR4Nmd05ZsSIpIEFCAMBgNut9sfIMrLy1FljnrYssSa+MHMUViaZwF9ubOCT7YcDXGpukbT9KhLSihEpAh4w6Bly5ZRW1vL3/72N2688UaWLFkS7LKJM5CWFMstM0ZiNvr+xB9uLuGbPZUhLlXX2BwSIIQIhYBmMS1YsIABAwbw8ccfY7fbefTRR5kwYUKwyybOUP+MBG6cNoKX1u5C0+HN9QdIjDeROyAl1EXrFJdHw+3RMBml1SpETwo44c2ECRMkKESgETkpXHXJUN5Yf3wh3W1zx9A/whbS2Zweko2Ru9ucEJEooABxzTXXdJgpdOXKld1eINH9zhuZSX3T8YV0L0XgQjqH00NinAlVjfyMtUJEioACxH333ef/2el0UlBQQGZmZtAKJbpfpC+k05H0G0L0tIACxMSJE9scX3TRRTJIHWFaFtJZbW52Ha6lqt7By+/vZsnsPMzGyFitLAFCiJ7VpVG/xsZGSkpKurssIsgMqsINVwwnJ9O3kO5weWQtpJMpr0L0rE6PQWiaRmlpKT/4wQ+CWjARHGajgZtnjOTJVduprnf4F9LNv2hIROxIZ3N4ZM9qIXpIp8cgDAYDAwYMICsrK2iFEsHVspDuyVXbabS7+XJnBUkWM/nnDgh10U7L5dHweDWMBpnyKkSwdWkMQkS+tKRYFs8cxdOrt+Nya3y0uZSkeDMTRoX/5AObw0OSRaa8ChFsAQWIyZMnd9j9oOs6iqKwcePGbi+YCL7+fSwsmjaCF9/bjabrvP2fAyTEmxg1MDXURTslu8tDQrwJNQK6xISIZAEFiIULF1JXV8f111+Pruu88cYbZGVlMWvWrGCXTwRZ7oAUrrlsKK9/vB9Nh39+tJelc/LIyUwMddFOStfB4fQSHytjEUIEU0AduV999RX//d//zahRo8jLy+PXv/4169evp3///vTv3z/YZRRBdk5uBjMmDgTA7dF48b3dVNbZQ1yqU7M53aEughBRL6AAUVFRQU1Njf+4pqaGysrITPwmOnbx+L5MGZsN+NYbPP/uThqaXCEu1cl5vDout+wVIUQwBdRGX7x4MfPnz2fq1KkArF+/njvuuCOoBRM9S1EUZl0wCKvNzXcHqqlrdPHi2l3cNnd02E4rtTk9mE2RschPiEgU0Ct/0aJFnHfeeXz11Vfous6iRYsYOXJksMsmepiqKFw3dRhNDjcHjjZwrNrGy+/v4ZaZo8Iyk6rT5duS1CB7kwgRFAG/sgYMGMC5557LzTffLMEhihkNKt+/cgR90+MBOHisgdfW7UMLw9XWOrIlqRDBFFCAWL9+PbNnz+bHP/4xAN999x0//OEPg1owETqxZiO3zBxFWmIMANuLa3gnTLcttTncYVkuIaJBQAHi8ccfZ+XKlSQlJQEwbtw4Dh8+HNSCidBKjDfzg9l5bbYtLfy6NMSlak/TweGSVoQQwRBwF1NGRkabY7NZVrJGu/SkWH4wcxQxzQPB6745wsZtZSEuVXuyJakQwRFQgLBYLFRVVflXU2/atInExNMvpNqwYQPTp09n2rRprFix4qTXrV27lpEjR/Ldd98FWGzRU/r1sfD96SMwNG/Us/rzYr7dVxXiUrXl9moy5VWIIAhoFtM999zDbbfdRmlpKTfddBPFxcU88cQTp3yO1+vloYce4vnnnycrK4trr72W/Px8hg8f3ua6xsZGXn75ZcaPH9/1WoigGtYvmesvz+WfH+1B12Hlx/uJMxsYGUYpOWTKqxDdL6AWxPjx43nppZf405/+xNKlSykoKGDs2LGnfE5RURGDBg0iJycHs9nM7NmzKSwsbHfd8uXLWbp0KTExMV2rgegRY4ekseDioQBous4/PtzLoTJriEt1XMuUVyFE9zltC8Lr9fK9732PN954g0svvTTgG5eXl5Odne0/zsrKoqioqM01O3bsoKysjKlTp/Lcc88FfO+0NEvA10aicK3f9ClDQFV465P9uL0aL72/m3tuPJcBWZ3L2xSs+sXFmUhOCO0HjYyM8M1h1R2kfr3LaQOEwWAgNTUVp9PZqU/5HU09bJ0RVtM0HnnkER555JGA79mipqap08+JFGlplrCu34TcPlTV2PhP0THsTg+P/WsLd8wbTZ/kuICeH8z61QL2lNiQLZzLyEiksjJ8WlXdTeoX2boS/AIagxg8eDCLFi1i+vTpxMfH+88vWrTopM/Jzs6mrOz4jJfy8nIyM4/vNdDU1MSePXu4+eabAaisrOTOO+/kiSeeYNy4cZ2uiOgZiqIwY9JA7C4vm3dV0GR381zBTu6YNybkn951oMkue0UI0V0CChBNTU3k5uZy4MCBgG88btw4iouLKSkpISsri4KCAv785z/7H09MTGTTpk3+45tuuol7771XgkMEUBSFBRcNweHysO1ADXWNLp4t2Mnt88aQ0LxuIlTsTg+WOKOk3xCiG5wyQPzhD3/g/vvv55FHHuGzzz7jwgsvDPzGRiMPPPAAS5cuxev1cs0115Cbm8vy5csZO3Ysl19++RkXXoSOqip8b+pwnK7d7C2tp6rewfPv7mTpnNHExYQuuZ8ONNo9JEsrQogzpuinyFNw1VVX8dZbb7X7OZSOVTVRVd0Y6mIETbiPQZzI5fHywru7KG6e0ZSTmcCSWXnEmDuectoT9VOA9OTYHt+3ujf0YUv9IldXxiBO+QpqHTsk343oiNlo4OYZI+mf4ZuZVFLRyMsf7MbtCd2UU99YhGwoJMSZOmWAcLlc7N+/n3379rX5ueVLCPAl9/vBzFFkpfpmMh042sCrH+7G4w1dkLC7vLg9srpaiDNxyi6m/Pz8kz9RUTpc+BZs0sUUvqw2FytW76C63gFA3qBUbpyW22bAuCfrZzKopCfH9sjvgt7RRSH1i1zdPs113bp1XS6M6H0S483cOjuPp1fvoNbqZOehWv69bh/X5+f6czn1JLdXw+bwEB8bnjviCRHuZC6g6FYpCTHcOjvPP4to24EaVn4Sug2HGu2usNzsSIhIIB+tRLdLS4rl1jl5PP3ODqx2N1v3VaMqCtdcOqzL99xbWsfmXRXUWp2kJsYwYVQmuQNSTvs8TQer3S3TXiPItoPVfFp0jMo6OxkpcVx0Vl/GDkkPdbF6JWlBiKDokxzHrXNG+xfObdlbxZsbDqB1YTbc3tI63v+yhOoGJ5oO1Q1O3v+yhL2ldQE93+70hHRWlQjctoPVvLH+AOW1djQdymvtvLH+ANsOVoe6aL2SBAgRNJmpcdw6Jw9L8xjAN3sqeeW9nZ0OEpt3VXTqfEdsTtlUKBJ8WnSsU+dFcEmAEEGVlRrPrXNG+weKPy86xpvrD3RqXKDW6uzU+Y44nB4Zi4gAlXX2k5x39HBJBEiAED0gOy2epa2CxDd7Knlj/f6A37BTEztOAniy8x3RkVZEJMhI6TgrcEZKz01XFsdJgBA9oiVIJMYfH5NY+cl+vAEEiQmjMjt1/mRsTo9kBAhzF53Vt1PnRXBJgBA9Jjstnp8uPBdL88D1t/uq+Pe6vafdCS53QArTJ+aQnhSDqkB6UgzTJ+YENIupNU3TcbhkdXU4GzsknWsuHUpWahyqopCVGsc1lw6VWUwhItNcRY/ql5HAbXNG8+wa3xTYbQdq8Hj2svCKXEzGk39eyR2Q0umA0BGbwxPSbLPi9MYOSZeAECakBSF6XGZqHLfNG+1fm7DrcC2vfLAbVw/kTnJ7NcnRJESAJECIkOiTHMft80b7B5r3ltbzwru7cLiCP5Bsc8hgtRCBkAAhQiY1MZbb542hT3NCveIyK8+s2UljkFN1O1zeLi3YE6K3kQAhQirZYua2uaPpm+7b6/xoVRNPr95OfWPgaxw6S8e3LkIIcWoSIETIJcabWTpnNAOzEgDfoqin3tlO1UkWTXUHWRMhxOlJgBBhIS7GyJJZeeQOSAagrtHFk+9sp7QiOHt/eLw6LrcMVgtxKhIgRNgwmwzcNH0k44amAb7B5GfW7Ag4KV9n2aUVIcQpSYAQYcVoULk+P5fJo7MAcHk0Xlq7m2/3VXX773K4vJKfSYhTkAAhwo6qKsy9cDBXTBgAgFfTeW3dPtZ/e6RbU2Xo0CPTaoWIVBIgRFhSFIX8cwdw1cVDaNmt9P0vS3jns+Ju/dQvayKEODkJECKsnZ+Xxfenj/Sn4di0o5xXPtiDs5sGmD2aTpMjuOsuhIhUEiBE2Bs1MJXb5o72J/nbdbiWFe9031qJRpsbj1d2nBPiRBIgREQYkJHAnfPH+PcLOFZt44m3t3G0qumM760DDU2uM76PENFGAoSIGGlJsfxw/hiG9/etlWiwuXnqne1sO1hzxvd2eTSZ9irECSRAiIgSF2Nk8cyRnN+8WZDbo/GPD/dQ+HXpGedXstpcMu1ViFYkQIiIY1BVFlw8hFmTB6E0z3Aq/LqUf32bgCyaAAAfmklEQVS094xWR2s61Fgdp93ASIjeQgKEiEiKonDRWX1ZPGMUsWYDANsO1vDkqu1UN3R9g3uPV6emwSmD1kIgAUJEuBE5KfxowVh/yvCyGhv/9+Z37CnpenoOr6ZTY3Xi9kiQEL2bBAgR8fqkxHHngrGMGpgK+FJovPjeLtZ90/VxCU3TqbU6JEiIXk0ChIgKcTFGvj99BJefNwAF39TVjzaX8vLa3di6uBBO05EgIXo1CRAiaqiKwuXnDeCmGSP94xK7S+r4+5vfUVJh7dI9JUiI3iyoAWLDhg1Mnz6dadOmsWLFinaPP//888yaNYu5c+eyePFijhw5EsziiF5i1MBUll09jv59LIBvb4kV7+zgP0VHu9TldDxIyP4RoncJWoDwer089NBDPPPMMxQUFLBmzRr27dvX5pq8vDzeeOMNVq9ezfTp0/nf//3fYBVH9DJpSb79ric1pw33ajrvfXGYl9bu6tKe15oONQ3OLndXCRGJghYgioqKGDRoEDk5OZjNZmbPnk1hYWGbayZPnkxcnC91wtlnn01ZWVmwiiN6IZNRZf5FQ7g+fzgxJl+X056Sev62sqhLmxDp+FZv1zU6z3hRnhCRwBisG5eXl5Odne0/zsrKoqio6KTXr1y5kksuuSSge6elWc64fOFM6te9pk60MHZEJs+u2kbxsQasdjfPv7uL/Ak5XHXZMExGQ+dvalBJS47FYGj7GSsjI7GbSh2epH69S9ACREcbuygty15PsGrVKrZt28Yrr7wS0L1ras48QVu4SkuzSP2CwAAsmTWKjzaX8J+tx9CBdZtL2La/iuvzh9M3vfNBq6q6kdSEGH8q8oyMRCoruzYYHgmkfpGtK8EvaF1M2dnZbbqMysvLyczMbHfd559/zpNPPskTTzyB2WwOVnGEwGhQmTFpEEvm5JFs8f1bq6i18//e2sbH3xzB28k8TJqmU2N1dNveFEKEm6AFiHHjxlFcXExJSQkul4uCggLy8/PbXLNjxw4eeOABnnjiCdLT04NVFCHaGNYvmbuuPYuzhvn+zXk1nQ83l/Dkqm2U19g6dS9dhzqrU3amE1FJ0btzk98TrF+/nocffhiv18s111zDnXfeyfLlyxk7diyXX345t9xyC3v27CEjIwOAvn378uSTT57ynseqmqiqbgxWkUNOuph61tZ9VbzzWbE/1bdBVZh6bn8uGd8Po6Fzn59y+qXgsjtP2pUa6XpDF0y016+zghoggkECRGQLx/pt3V/FexsP0WA7PoU1NTGGtMQYnG4vqYkxTBiVSe6AlFPeJy3NQmODneQEMwa1c8Fl28FqPi06RmWdnYyUOC46qy9jh4RHq3rNxmI+2XKEJocHS6yRy87pz5wLBoe6WN1OAkR7QRukFiIS7C2t4z9bj2GJM6GqKg1NzuaFcU5qrU4ssUY8Xp33vywBOG2QcHk0quocxJoNxMea/APYp7LtYDVvrD/gPy6vtfuPQx0k1mwsZs1nxYBvkkmjze0/jsYgIdqSVBuiV9u8qwLwvfnFxxrJSI1HVY93ETU5PFTU2bE7PXy1szyge+qA3eWlusFBTYMDu9PT4ay+Fp8WHevU+Z70yZaOsxuc7LyILtKCEL1ardXZ5tigKhhU8IUIBa+mN2d29Q1EV9TZyWzeFzsQLo+Gy+PCaleIjzESF2No1/1UWWfv8LmVdV3f16K7nGzVeVMXVqOLyCMtCNGrpSbGtDtnUFVMRgMZqXEkxJn8551uL4+/XsS7XxzC4ercrCVN02m0u6mqc1DX6Gyz813GSQJORkpsp35HMLSuf2uWk5wX0UUChOjVJoxqvzYnPtaIJdaIqigkWcxkpMRhNvleKpqu82nRMf787618tbO803tY6/j2q6ixOqmqt2NzeLhwXHaH1150Vt9O16e7XXZO/06dF9HF8OCDDz4Y6kJ0RqPNjc3uCnUxgiYuzow9ipvv4Va/9KRYUhNjqLM6cbq8pCXFkH/uAEYOTPWfy0iJZfrEHEYPTqOkohGHy4vbo7HrcB3bD9aQlhRLenJsp+um6b5WSWKciczUOOqsThwuL5mpccyYNDDkA9Tg27EPBY5WN+HxaljiTFw5cWBUDlBbLDHYbNH73mKxtG8tn45Mcw0z4TgNtDtFev3cHo0NW4+yYevRNntEDO+fzHVXjCAxpgt5nVoxqAoxZgOxJgNm05ndq7v1hmmg0V6/zpIAEWYi/Q30dKKlfg1NLj7aXMLXuytp/QIaMySNaefndGog+2QMqkKs2UCs2RjQdNlg6w1voNFev86SWUxCdEGSxczVlw7jgrHZvP/lYfaU1AOw/WANO4prOHt4H6ae058+ZxAovJpOk8NDk8ODUVWIjTESYzKERbAQvYMECCHOQN90C7fMzOPgsQYKvznCgSP16Dps2VvFt/uqOHt4Hy47p/9JZyoFytM8C6rR7kZVwGT0BQqTQcVsUqM2vYcIrYjrYrI53JRXWNF0HV33TR/Udd9cdQ3QNZ2IqtAJoqUL5mSiuX6pqfFs/PYIH31dytGq43VUgDFD07js7P7069P9e2EogNlkIMakYjYZOp1DKlC9oQsm2uvXWRHXgoiPNZ10bnYLrTlg6Hrbn0GnZVairjcHEp3m78cDi64f389C03R/MIrkwCOCT1EURg1KZeTAFHYdqqXw61KOVtvQgW0Hath2oIYROclcdFY/hvVL6rZP/Tq+2VC+tONuVFUhxqhiMqoYDSpGo4oqLQzRBREXIAKhKgqqoftfELquNwec4z+3BKE2rZnmABMtrRrROYqikDc4jVGDUtldUscnW45wuNw3sWJPST17Surpmx7PRWf1ZdzQ9G7/xK9pOnaXF7vr+GI8o6pgMhkwtwocQpxOxHUxARHZDGwTULTWgcUXUDRdR9d0UtMsVFU1tgs+0SKau5hOVjdd1zl4zMr6b4+wt7S+zWOJcSYmjs5iYl4mifE9t2GWquBvXZgMKkaDgsFw6pZGb+iCifb6dVZUtiDCkaIoGFpefKeY3p6eHIfWQRqH1l1lHbZemgOMt/lLk1ZL2FAUhaH9khjaL4myGhufFh1j674qvJqO1e6m8OtSPtlyhLFD05iYl8Xg7MSgDzprekueKK3NeV8uKl+wMKgKJoOvxdE6gaHoPaQFEWa681OMV9Pwen0BQ9eh5T2n9diM/7um4+2B1kpvbEF0pKHJxaad5Xy5s6Jd4ruMlDgm5mVydm4fLLHhkfPIoCpkZyVRX2fDoCq+1odBiarZU9KCaE8CRJgJ9T9STdPxahoeb6uWSOtxlVY/n+5fTstbR+vLJEC05fFqfLe/mi92lFNS0XYBqEH1DXpPGJnB8AEpGEL8Kb6j+rUEC4OqoCi+1pKqKBgMiv+xSBHq116wSReTOGOqqqCqBkwB/stoGUdRUGj+H0CbT5atWzKJ8WZsjQ48Xg1NP359i+MzyaJr7OVkjAaVc0ZkcM6IDI5WNfHlznK+3VeFy63h1XS2H6xh+8EaEuJMjB+Wztm5fejXxxI2n9x9XZrekz6ugL/LymhoaXm0bX20fPBQUKQrK8xICyLM9IZPMYHW7/isMd/MMa318YlrYPTjU5hDpbtaR063l20Hqvl6dyXFZe3/W/VJjuWsYemMG5ZOVmr8Gf++QHV3609VlXYfBBSlZRzEN+6hNh+3tExU1fe8zm7pGoje8NrrLGlBiLDVMrDfmV6K1uMq/rGWDoLKic/xhtGgfozJwHkjMzlvZCZVdXb/quyWzY2q6h2s++YI6745QnZaPGOGpDFmSBpZqXFh07IIREep0nUdPF4dj/fkrRLwtUzU5gF1RfG1PFqCi6oo/haKqhzv+hKdJy2IMNMbPsWEa/1axl80zdd1ounNs8K8Gpqm4zlNEyWY4yuarnO43Mq3e6vYdrAGm6P9TLf0pFjyBqeSNyiVgVmJ3T5mEcnjR0rz/yn4goVC+4Wv6ekJ1NQ0+rpLOT6pQ2kJMs0/t7mPAgqK/1q1+XG1+RcqSvO6rDDoOpNB6igQzm+g3SGS6+fxajQ5PDicng5bGz31BurVNA4cbaBofzU7imuxO9sHi7gYIyNzUhiRk0JuTnK3zIaK5AARiGDWTwEU1dci9o3zKc2tneaWsnq8FdQScDrT6mm9iPd4l+vxljPA8MGd319EupiECJDRoJJsMZMQZ8Tm8GBzekIykG5QVXIHpJA7IIUFF2scPGb1Z5G12nxTZu1OD9/u83VNKUD/DAvDB6QwvH8SA7MSI2p2UTTQ8WVU0NDh1L1nfv5WT3OLp3Xw8K+HIrBMDV1tTUqAEKKTDKpKYrwZS5wJh9OXjtsbohFyg6oyvH8yw/snM/fCwRytamLXoVp2HarlaLUN8L05lVY2UVrZxCdbjmAyqgzOTvQv3uvXJyHkU2hFe/5ccW0+hfTsvzMJEEJ0kaooxMeaiI814dU0UpJj8ThcuD0aTo/W6f2qu6M8AzISGJCRwBUTcmiwudhbUseekjr2HanH7vR9dHV7NPaW1vvTfpiNKgOzEhmUncjgvonkZCSE3W52IjQkQAjRDQyqSqzZSHyrvn6PV/NlWW3ew7qn2xhJ8Wb/bChN0zla3cS+0nr2HanncLkVj9dXIpdHY98R33nwDbBmpcUzMCuRnMwEBmQm0Cc5todLL8KBBAghgqRlUZgl1oSm6TjdXmwOD26vdvondzNVPd66uOyc/ni8GiUVjRw42kBxWQOHyxv9e2xrOhyrtnGs2samHeWAb+rt4L5JZKbE0q+PhX59LKQnx0oa8SgnAUKIHqCqCnExRuJijDjdXprs7naJ8nqS0aAypG8SQ/omAb6ZUceqbBSXWTlcYaWkvJH6Jpf/eqfby+7Dtew+fPweJqNKdlo8fdPjyU6LJystnqzUeOJj5W0lWshfUogeFmMyEGMy4PZo2F0eHC5vj49XnMigqgxo7k6CvgDUNzqbB7cbKa1s5GiVrc2UWrfH1wo5MYdUYryJzNQ4MlLiyEzxfe+TEkdSvEkWrEUYCRBChIjJqGIymkmKB5fbi8PlxeHyhDxlSIvkhBiSE2IYMyQN8G2puv9QDUeqmjha1URZtY1jNTYaWrU0AKw2N1abm/1HGtqcN5tU+iTHkZ4UQ3pSLOnJsaQlxZKWGEOixSzdVWFIAoQQYcBsMmA2GUiymP0D2w536FsWrSmK4ntDT4pl3NDji66aHG7Ka2yU19gpq7FRUWunos7mnzXVwuXWONocXE5kNCikJMSQmhjT5ntygpmUhBiSLKag5F8SpyYBQogw09IFlYSvG8fp9uJyh2YmVCAssSaG9ktmaL9k/zld12m0u6mss1NZ56Cq3k5VnYOqBge1DU60E1YYerw6VfUOquodHf4OBUiIM5GUYCYp3kySxfc9Md7U/OX72RJrCou0FtFCAoQQYczUvIc0cSY0Xcft1nB5fMHC7dXCNiW6oijNb9rmNoEDfHmu6hqdVNc7qLU6qWlwUNPgpK7RSa3Via2D1CE6YLW7sdrdHOHk6TAUID7WSEKcCUucyfc91oQlzogl1kR8rJH4WN/PcTFG4mOMvv++okMSIISIEKqiEGM2EGM+vojN49VwuX3Bwu32hlVW2pMxqIpvDCKp47UVTpeX2kYn9Y1O6hpd1Dc6qW9y+b8amlz+Kbkn0oEmh291O7X2gMpjMqrExRhJjDdjMirEmY3ExRiIMxuJMRuINRuJNRuav1rOGfwtPZPp1Ht5RzIJEEJEsJa1Fq15vC0bNGn+XQEjaa/yGLOB7DTf1NmO6LpvTUmDzY21yYXV7qbR5sZqc9Fod7f5arJ72nVnncjt0XB7XO0G2zvDbFKJMRowmw3EGFX/mJLZqGI2qZiMLT/7vptafxlUjCf+3LwXeEcbLPWkoAaIDRs28Pvf/x5N07juuuu4/fbb2zzucrm499572b59OykpKTz22GMMGDAgmEUSImptO1jNp0XHqKyzk5ESx0Vn9WXskHSeLdjBVzsrcHs1TAaVCaMymDAqk8++O0ZlnYP0pFgm5mUyIieVnYdq+HJnBTUNDlITY5gwKpPSyka+3FGOzeUl3mxg4ugspp7T8et0b2kdm3dVUGt1+p+fOyDlpOc7c4+Pt5T6yuH0EB9jDKgcNQ0Oki1mRg1KxeH2sv1ADQ1NLmLMBjJS4og1+xIv2l0e3B4Nq82F3enF6Q4wo14zl9vXkuOE/cW7U8uOfAaDirFlq1fD8f3BW2//alBV/7avBoOKyaAwOjez078zaOm+vV4v06dP5/nnnycrK4trr72Wv/zlLwwfPtx/zauvvsru3bt56KGHKCgo4MMPP+Svf/3rae8dqemiAxHJ6bADEc31C2Xdth2s5o31B9qdt8Qa2Vlc2+acjq+fPiMlrs3580Zm8PXuyuPX6b6xgia7u3njHcWfOG7axBwuPzfHn2Za12F3SS1rN5W0K8OYIalsP1jb7vz0iTntgsTe0jre/7L9PbLT4vhuf3W785ed279dkOjoHg6XBwWIMbf9TNy6DK3Tfe8uqeX9TSX+LXVbdjQcOzSNFEsMDpcviLRMIHA0p1NpOdd6ckFLSpNQW/3n+Z1+TtBaEEVFRQwaNIicnBwAZs+eTWFhYZsAsW7dOpYtWwbA9OnTeeihh3z7G0dpf54QwfJp0bEOz+861P6NGcDewYZDn2w5QmK82X+sKApNdl9Kc0Vt28XxxfZyrr10eJvnv7F+v3/At/Xnzs27Kvz39Z/VoWh/FRNGHv9Uq+uwdV9Vh7OQOgoOAF/trGDW5MG0Hq3fsreSE99CWga+Y2PavuVt2VPJ2CG+KbuxZl/3D8DWvVW+T+An7JpubXIx/8IhHZblZHxpVjTcHi8uj2+syOXRfD/7v7y+caRW57ya7p+M4PF/6Xg8zT9rvp+9mt7mca3luLlL8UwELUCUl5eTnZ3tP87KyqKoqKjdNX37+lZtGo1GEhMTqa2tJS0t7ZT37srOSJFE6he5QlW32kZXh7NxNJ12b5bovjfqE69vcnhIO2HguKX/viU4tHy3OTzt6nqyMticXtKT22eHtdo9DByQ2u5crLn9tR5Nx2xsf97h8rbbCKfJ4SXuhECga4BCu/M2l5cRQ/v4j9OTfa2qJueOdte2XD9yWEa789EqaAGio56rE1sGgVwjhDi9v/zk0lAXoVvKEE33iAZBmwCcnZ1NWVmZ/7i8vJzMzMx21xw75msaezwerFYrKSkdD1wJIYToWUELEOPGjaO4uJiSkhJcLhcFBQXk5+e3uSY/P5+33noLgPfff5/JkydLC0IIIcJE0GYxAaxfv56HH34Yr9fLNddcw5133sny5csZO3Ysl19+OU6nk5///Ofs3LmT5ORkHnvsMf+gthBCiNAKaoAQQggRuSQJiRBCiA5JgBBCCNGhsM7F5HQ6WbRoES6Xy78y+6677qKkpIS7776b+vp6Ro8ezR//+EfMZvPpbxiGWsZnsrKyeOqpp6Kqbvn5+VgsFlRVxWAw8Oabb1JXV8dPf/pTjhw5Qv/+/fnrX/9KcnLy6W8WhhoaGvj1r3/Nnj17UBSFhx9+mCFDhkRF/Q4cOMBPf/pT/3FJSQl33XUXCxYsiIr6vfDCC7z++usoisKIESN45JFHqKioiJrX3osvvsjrr7+Orutcd9113HLLLV177elhTNM0vbGxUdd1XXe5XPq1116rb9myRb/rrrv0NWvW6Lqu67/5zW/0V199NZTFPCPPPfecfvfdd+u33367rut6VNVt6tSpenV1dZtzjz76qP7UU0/puq7rTz31lP7HP/4xFEXrFvfee6/+2muv6bqu606nU6+vr4+q+rXweDz6lClT9NLS0qioX1lZmT516lTdbrfruu57zb3xxhtR89rbvXu3Pnv2bN1ms+lut1tfvHixfvDgwS797cK6i0lRFCwWC+BbJ+HxeFAUhS+++ILp06cDcNVVV1FYWBjKYnZZWVkZn3zyCddeey3gWzgYLXU7mcLCQhYsWADAggUL+Oijj0Jcoq5pbGzkq6++8v/tzGYzSUlJUVO/1jZu3EhOTg79+/ePmvp5vV4cDgcejweHw0FGRkbUvPb279/P+PHjiYuLw2g0cv755/Phhx926W8X1gECfH/I+fPnM2XKFKZMmUJOTg5JSUkYjb7esezsbMrLy0Ncyq55+OGH+fnPf47avJVibW1t1NStxa233srVV1/Nv//9bwCqq6v9CyYzMzOpqakJZfG6rKSkhLS0NH7xi1+wYMECfvWrX2Gz2aKmfq0VFBQwZ84cIDr+fllZWSxZsoSpU6dy0UUXkZCQwJgxY6LmtTdixAg2b95MbW0tdrudDRs2UFZW1qW/XdgHCIPBwKpVq1i/fj1FRUUcONA+Y2UkLq77+OOPSUtLY+zYsae8LhLr1uKf//wnb731Fk8//TSvvvoqX331VaiL1G08Hg87duxg4cKFvP3228TFxbFixYpQF6vbuVwu1q1bx4wZM0JdlG5TX19PYWEhhYWF/Oc///G/iZ4oUl97w4YNY+nSpSxZsoSlS5cycuRIDIb2eawCEfYBokVSUhKTJk3i22+/paGhAY/Hl52xrKysXQqPSPDNN9+wbt068vPzufvuu/niiy/4/e9/HxV1a5GVlQVAeno606ZNo6ioiPT0dCoqKgCoqKg4bWLGcJWdnU12djbjx48HYMaMGezYsSNq6tdiw4YNjBkzhj59fAntoqF+n3/+OQMGDCAtLQ2TycSVV17Jli1bouq1d9111/HWW2/x6quvkpKSwqBBg7r0twvrAFFTU0NDQwMADoeDzz//nGHDhjFp0iTef/99AN566612KTwiwT333MOGDRtYt24df/nLX5g8eTJ//vOfo6JuADabjcbGRv/Pn332Gbm5ueTn5/P2228D8Pbbb3P55ZeHsphdlpGRQXZ2tr9Fu3HjRoYNGxY19WtRUFDA7Nmz/cfRUL9+/fqxdetW7HY7uq6zceNGhg8fHjWvPfB1BQIcPXqUDz74gDlz5nTpbxfWK6l37drF/fffj9frRdd1ZsyYwbJlyygpKeGnP/0p9fX15OXl8ac//Slip6MBbNq0ieeee84/zTUa6lZSUsJ//dd/Ab5xpDlz5nDnnXdSW1vLT37yE44dO0bfvn1Zvnx5xCZo3LlzJ7/61a9wu93k5OTwyCOPoGla1NTPbrdz2WWX8dFHH5GY6EvtHS1/v8cff5x3330Xo9FIXl4ev//97ykvL4+K1x7AjTfeSF1dHUajkV/84hdccMEFXfrbhXWAEEIIETph3cUkhBAidCRACCGE6JAECCGEEB2SACGEEKJDEiCEEEJ0KKyzuQpxKtdddx0ulwu3201xcTG5ubkAjB49mkceeSTEpQvM9u3bKSkpiaqVyiJ6yDRXEfFKS0u55ppr2LRpU6iL0o7H4/Hn9+nI66+/zueff85jjz3W7fcW4kzJvy4RlVauXMm//vUvvF4vSUlJ/Pa3v2Xw4MG8/vrrrF27FovFwp49e+jbty+//OUvefTRRykpKWH8+PE8+uijKIrCz372M+Li4jh8+DBlZWVMmjSJ3/zmN5hMJqxWKw8//DB79+7F6XQyZcoU7rvvPlRVZeHChUycOJEtW7YQHx/P448/7l8k6HQ6GT9+PL/97W9paGjg//7v/2hqamL+/PlMmjSJRYsWceONN/LZZ58BcOjQIf/xoUOHWLhwIddffz1ffPEFV199NfPnz+cvf/kLmzdvxuVykZeXx4MPPkhcXFyI/wIiKgQpJbkQPaakpESfOHGi//iLL77Q77jjDt3pdOq6ruuFhYX6okWLdF3X9ddee02fOHGiXlZWpuu6ri9ZskRfsGCBbrVadZfLpc+aNUv/4osvdF3X9XvuuUefP3++3tTUpLtcLv3mm2/W//GPf+i6ruv33Xefvnr1al3Xdd3r9ep33XWXvnLlSl3Xdf2GG27Qf/SjH+kej8f/eF1dnf/nu+++27+PxGuvvab/5Cc/8Ze9uLhYnzJlSofHxcXF+ogRI/S1a9f6H3/88cf9Of51XdcfeeQRffny5Wf2H1SIZtKCEFFn3bp17Nixg+uuuw7w7bPR1NTkf/y8887zJxIcPXo0DoeDhIQEAEaOHMnhw4eZNGkSALNmzSI+Ph7w5dD/5JNPWLhwIR9//DHbt2/n6aefBny5wgYOHOj/HXPnzvVn0NQ0jRUrVvDpp5+iaRp1dXVd3oUtPj7ev2dBS13tdjsFBQWAL/vqmDFjunRvIU4kAUJEHV3X+d73vseyZcs6fDwmJsb/s6qq7Y5bMnp2dN+WFNCapvHUU0/Rr1+/Dq9tCSoAq1atoqioiH/84x9YLBb+/ve/c+zYsQ6fZzAY0DTNf+x0Ok9635Yy/e53v+P888/v8H5CnAmZ5iqiTkvWypYNX7xeL9u2bevSvd577z3sdjtut5vVq1f7Wxb5+fmsWLECr9cL+DIPl5SUdHgPq9VKamoqFouF+vp6/6d9AIvFgtVq9R9nZmbicDj891qzZs1p6/rcc8/5A0ljYyP79+/vUl2FOJEECBF1Jk+ezLJly7jjjjuYN28ec+fO5ZNPPunSvc477zzuvPNO5syZQ05Ojn+L0d/85jdomsb8+fOZO3cut912G5WVlR3e46qrrqKuro45c+Zw9913t/m0f+GFF2K1Wpk3bx4PP/wwZrOZ+++/n8WLF3PTTTdhMplOWb4f/vCHDBs2jGuvvZa5c+eyaNEiDh482KW6CnEimeYqxEn87Gc/47zzzmPhwoWhLooQISEtCCGEEB2SFoQQQogOSQtCCCFEhyRACCGE6JAECCGEEB2SACGEEKJDEiCEEEJ06P8H8uC1Bmw5yFwAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "sns.set(color_codes=True)\n",
+    "plt.xlim(30,90)\n",
+    "plt.ylim(0,1)\n",
+    "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
+   ]
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Hide code",
+  "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.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}