Upload New File

parent 01775af3
{
"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 Python 3 language using the pandas, statsmodels, and numpy library."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]\n",
"uname_result(system='Windows', node='MGDONDON', release='7', version='6.1.7601', machine='AMD64', processor='Intel64 Family 6 Model 94 Stepping 3, GenuineIntel')\n",
"IPython 6.5.0\n",
"IPython.core.release 6.5.0\n",
"_csv 1.0\n",
"_ctypes 1.1.0\n",
"decimal 1.70\n",
"argparse 1.1\n",
"backcall 0.1.0\n",
"colorama 0.3.9\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.7.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.4.0\n",
"ipywidgets._version 7.4.0\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.3\n",
"matplotlib.backends.backend_agg 2.2.3\n",
"numpy 1.15.4\n",
"numpy.core 1.15.4\n",
"numpy.core.multiarray 3.1\n",
"numpy.lib 1.15.4\n",
"numpy.linalg._umath_linalg b'0.1.5'\n",
"numpy.matlib 1.15.4\n",
"pandas 0.23.4\n",
"_libjson 1.33\n",
"parso 0.3.1\n",
"patsy 0.5.0\n",
"patsy.version 0.5.0\n",
"pickleshare 0.7.4\n",
"platform 1.0.8\n",
"prompt_toolkit 1.0.15\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.9.0\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.7\n",
"zlib 1.0\n",
"zmq 17.1.2\n",
"zmq.sugar 17.1.2\n",
"zmq.sugar.version 17.1.2\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/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGBNJREFUeJzt3XuQnXWd5/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>Mon, 12 Nov 2018</td> <th> Deviance: </th> <td> 3.0144</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>20:18: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: Mon, 12 Nov 2018 Deviance: 3.0144\n",
"Time: 20:18: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.0850** 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>Mon, 12 Nov 2018</td> <th> Deviance: </th> <td> 18.086</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>20:18: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: Mon, 12 Nov 2018 Deviance: 18.086\n",
"Time: 20:18: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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Temperature Intercept Frequency\n",
"0 30.0 1 1.0\n",
"1 30.5 1 1.0\n",
"2 31.0 1 1.0\n",
"3 31.5 1 1.0\n",
"4 32.0 1 1.0\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGeFJREFUeJzt3X901HV+7/Hn2wSW8LuAS8WgsC3Ftf4AEoLKXRtcBdyzIt6iyFq3bpdl72lZ67Vyj5xrV9fqOfc2nlv3bq2VKrXVowE9guwe7ga1pNt6/AEsCAINoMtqwm5RlB9xAyThff/4fhOHMMlMJjOZmQ+vxzk5zPc7n+/3+3nny7xm8pnvfMbcHRERCcs5+e6AiIhkn8JdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRAKcPdzFaa2UEze7eb+83M/q+Z7TOz7WY2LfvdFBGR3kjnlfvTwNwe7r8emBT/LAEe73u3RESkL1KGu7v/DPikhyY3Av/skTeBkWZ2XrY6KCIivVeahX2cD3yYsNwYr/tV14ZmtoTo1T1lZWUV48eP7/XBPjnunGh3LLO+FhwH1VJgQqkDVEuhGnAOjC7L7C3PPXv2fOzu56Zql41wT/b7TjqngbuvAFYAVFZW+ubNmzM6YH19PdXV1RltW2hUS+EJpQ5QLYWqL7WY2S/TaZeNq2UagcSX4OXAgSzsV0REMpSNcF8HfDO+auYK4Ii7nzEkIyIi/SflsIyZPQ9UA2PMrBG4HxgA4O5/D6wHvgbsA34DfCtXnRURkfSkDHd3X5Tifgf+LGs9EpGi0NraSmNjI8ePH++X440YMYLdu3f3y7FyLZ1aBg0aRHl5OQMGDMjoGNl4Q1VEzkKNjY0MGzaMCRMmYJb761iOHTvGsGHDcn6c/pCqFnfn0KFDNDY2MnHixIyOoekHRCQjx48fZ/To0f0S7GcbM2P06NF9+qtI4S4iGVOw505ff7cKdxGRAGnMXUSKVklJCZdeemnn8tq1a5kwYUL+OlRAFO4iUrTKysrYtm1bt/e3tbVRWnp2xpyGZUQkKE8//TQ333wzN9xwA7NnzwagpqaG6dOnc9lll3H//fd3tn344YeZPHky1157LYsWLeKRRx4BoLq6mo7pUT7++OPOvwba29tZtmxZ576eeOIJ4PPpBBYsWMBFF13EbbfdRnSVOGzatImrrrqKyy+/nKqqKo4dO8acOXNOe1KaOXMm27dvz+rv4ex8ShORrPrBj3ey68DRrO7z4nHDuf+G3++xTUtLC1OmTAFg4sSJrFmzBoA33niD7du3M2rUKDZs2MDevXt5++23cXfmzZvHz372M4YMGUJtbS1bt26lra2NadOmUVFR0ePxnnrqKUaMGMGmTZs4ceIEM2fO7HwC2bp1Kzt37mTcuHHMnDmT119/naqqKhYuXMiqVauYPn06R48epaysjG9+85s8/fTTPProo+zZs4cTJ05w2WWXZeG39jmFu4gUre6GZa677jpGjRoFwIYNG9iwYQNTp04FoLm5mb1793Ls2DFuuukmBg8eDMC8efNSHm/Dhg1s376dF198EYAjR46wd+9eBg4cSFVVFeXl5QBMmTKF/fv3M2LECM477zymT58OwPDhwwG46aabmDlzJjU1NaxcuZI77rijb7+IJBTuItJnqV5h97chQ4Z03nZ3li9fzne/+93T2jz66KPdXm5YWlrKqVOnAE671tzd+dGPfsScOXNOa19fX88XvvCFzuWSkhLa2tpw96THGDx4MNdddx0vv/wyq1evJtMZcnuiMXcRCdqcOXNYuXIlzc3NADQ1NXHw4EGuvvpq1qxZQ0tLC8eOHePHP/5x5zYTJkxgy5YtAJ2v0jv29fjjj9Pa2grAnj17+Oyzz7o99kUXXcSBAwfYtGkTEH0yta2tDYDFixdz5513Mn369M6/MrJJr9xFJGizZ89m9+7dXHnllQAMHTqUZ599lmnTprFw4UKmTJnChRdeyFe+8pXObe655x5uueUWnnnmGa655prO9YsXL2b//v1MmzYNd+fcc89l7dq13R574MCBrFq1iu9973u0tLRQVlbGq6++CkBFRQXDhw/nW9/K0VyL7p6Xn4qKCs/Uxo0bM9620KiWwhNKHe65rWXXrl0523cyR48ezen+77//fq+pqcnpMTocPXrUm5qafNKkSd7e3t5tu2S/Y2Czp5GxGpYREelnzz33HDNmzODhhx/mnHNyE8MalhERAR544IF+O9Y3vvGNM97gzTa9cheRjLkn/bpkyYK+/m4V7iKSkUGDBnHo0CEFfA54PJ/7oEGDMt6HhmVEJCPl5eU0Njby0Ucf9cvxjh8/3qewKyTp1NLxTUyZUriLSEYGDBiQ8bcEZaK+vr7zU6bFrj9q0bCMiEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISoLTC3czmmlmDme0zs3uT3H+BmW00s61mtt3Mvpb9roqISLpShruZlQCPAdcDFwOLzOziLs3uA1a7+1TgVuDvst1RERFJXzqv3KuAfe7+vrufBGqBG7u0cWB4fHsEcCB7XRQRkd6yVN9cbmYLgLnuvjhevh2Y4e5LE9qcB2wAfgsYAlzr7luS7GsJsARg7NixFbW1tRl1urm5maFDh2a0baFRLYUnlDpAtRSqvtQya9asLe5embKhu/f4A9wMPJmwfDvwoy5t7gb+Ir59JbALOKen/VZUVHimNm7cmPG2hUa1FJ5Q6nBXLYWqL7UAmz1Fbrt7WsMyjcD4hOVyzhx2+TawOn6yeAMYBIxJY98iIpID6YT7JmCSmU00s4FEb5iu69LmA+CrAGb2ZaJw/yibHRURkfSlDHd3bwOWAnXAbqKrYnaa2YNmNi9u9hfAd8zsHeB54I74zwcREcmD0nQauft6YH2Xdd9PuL0LmJndromISKb0CVURkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEAqRwFxEJkMJdRCRACncRkQAp3EVEApRWuJvZXDNrMLN9ZnZvN21uMbNdZrbTzJ7LbjdFRKQ3SlM1MLMS4DHgOqAR2GRm69x9V0KbScByYKa7f2pmX8xVh0VEJLV0XrlXAfvc/X13PwnUAjd2afMd4DF3/xTA3Q9mt5siItIb5u49NzBbAMx198Xx8u3ADHdfmtBmLbAHmAmUAA+4+0+T7GsJsARg7NixFbW1tRl1urm5maFDh2a0baFRLYUnlDpAtRSqvtQya9asLe5emapdymEZwJKs6/qMUApMAqqBcuDfzOwSdz982kbuK4AVAJWVlV5dXZ3G4c9UX19PptsWGtVSeEKpA1RLoeqPWtIZlmkExicslwMHkrR52d1b3f0XQANR2IuISB6kE+6bgElmNtHMBgK3Auu6tFkLzAIwszHA7wHvZ7OjIiKSvpTh7u5twFKgDtgNrHb3nWb2oJnNi5vVAYfMbBewEVjm7ody1WkREelZOmPuuPt6YH2Xdd9PuO3A3fGPiIjkmT6hKiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFSuIuIBEjhLiISIIW7iEiAFO4iIgFKK9zNbK6ZNZjZPjO7t4d2C8zMzawye10UEZHeShnuZlYCPAZcD1wMLDKzi5O0GwbcCbyV7U6KiEjvpPPKvQrY5+7vu/tJoBa4MUm7vwL+Gjiexf6JiEgGzN17bmC2AJjr7ovj5duBGe6+NKHNVOA+d/9DM6sH7nH3zUn2tQRYAjB27NiK2trajDrd3NzM0KFDM9q20KiWwhNKHaBaClVfapk1a9YWd0859F2axr4sybrOZwQzOwf4G+COVDty9xXACoDKykqvrq5O4/Bnqq+vJ9NtC41qKTyh1AGqpVD1Ry3pDMs0AuMTlsuBAwnLw4BLgHoz2w9cAazTm6oiIvmTTrhvAiaZ2UQzGwjcCqzruNPdj7j7GHef4O4TgDeBecmGZUREpH+kDHd3bwOWAnXAbmC1u+80swfNbF6uOygiIr2Xzpg77r4eWN9l3fe7aVvd926JiEhf6BOqIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBSutqGZFcWbu1iZq6Bg4cbmHcyDKWzZnM/Knn57tbkiadv8KlcJe8Wbu1ieUv7aCltR2ApsMtLH9pB4ACogjo/BU2DctI3tTUNXQGQ4eW1nZq6hry1CPpDZ2/wqZwl7w5cLilV+ulsOj8FTaFu+TNuJFlvVovhUXnr7Ap3CVvls2ZTNmAktPWlQ0oYdmcyXnqkfSGzl9h0xuqkjcdb7rpaovipPNX2BTuklfzp56vMChiOn+FS8MyIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBUriLiARI4S4iEiCFu4hIgBTuIiIBSivczWyumTWY2T4zuzfJ/Xeb2S4z225mr5nZhdnvqoiIpCtluJtZCfAYcD1wMbDIzC7u0mwrUOnulwEvAn+d7Y6KiEj60nnlXgXsc/f33f0kUAvcmNjA3Te6+2/ixTeB8ux2U0REesPcvecGZguAue6+OF6+HZjh7ku7af+3wK/d/aEk9y0BlgCMHTu2ora2NqNONzc3M3To0Iy2LTSqpfCEUgeolkLVl1pmzZq1xd0rU7UrTWNflmRd0mcEM/sjoBL4g2T3u/sKYAVAZWWlV1dXp3H4M9XX15PptoVGtRSeUOoA1VKo+qOWdMK9ERifsFwOHOjayMyuBf4n8AfufiI73RMRkUykM+a+CZhkZhPNbCBwK7AusYGZTQWeAOa5+8Hsd1NERHojZbi7exuwFKgDdgOr3X2nmT1oZvPiZjXAUOAFM9tmZuu62Z2IiPSDdIZlcPf1wPou676fcPvaLPdLJCNrtzZRU9fAgcMtjBtZxrI5kwHOWDd/6vn9cuxcHCcd963dwfNvfchdl7Ty7eXrWTRjPA/NvzQvfZH8SCvcRYrB2q1NLH9pBy2t7QA0HW5h2QvvgEFru3euW/7SDoCsBm+yY+fiOOm4b+0Onn3zg87ldvfOZQX82UPTD0gwauoaOsO1Q+sp7wz2Di2t7dTUNeT82Lk4Tjqef+vDXq2XMCncJRgHDrfkpG1f9pft46SjvZvPrnS3XsKkcJdgjBtZlpO2fdlfto+TjhJL9tGU7tdLmBTuEoxlcyZTNqDktHUDzjEGlJweamUDSjrfaM3lsXNxnHQsmjG+V+slTHpDVYLR8cZlPq6W6e7Y+bhapuNN044x9hIzXS1zFlK4S1DmTz0/aaD2R8h2d+x8eGj+pTw0/1Lq6+t577bqfHdH8kDDMiIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAVK4i4gESOEuIhIghbuISIBK02lkZnOBHwIlwJPu/r+63P8F4J+BCuAQsNDd92e3qyLhWru1iZq6Bg4cbmHcyDKWzZnMC5s/4PX3PulsM/N3RnFz5QVntAPOWLf5l5/w/FsfctclrXx7+XoWzRjPQ/MvTeu4yfY3f+r5afe749jt7pSY5eTYybbtro9nq5ThbmYlwGPAdUAjsMnM1rn7roRm3wY+dfffNbNbgf8NLMxFh0VCs3ZrE8tf2kFLazsATYdbuGvVtjPavf7eJ6eFfdPhFpa9+A44tJ7yznV3r9rGqYTt2t159s0PAE4L2WTHXfbCO2DQ2v75/pa/tAPgjPBMtn1/HDvZtt318WyWzrBMFbDP3d9395NALXBjlzY3Av8U334R+KqZWfa6KRKumrqGzqDqrdZ27wz2Dqe6afv8Wx+mPG7rKe8M1w4tre3U1DWcsb9k2/fHsZNt210fz2bm7j03MFsAzHX3xfHy7cAMd1+a0ObduE1jvPxe3ObjLvtaAiyJFycDmZ6NMcDHKVsVB9VSePq1joG//bsVudp3+2+OUDJ4ROfyyV/v25LpcRO37ev2GW47Bvi4p2279rGA9eX/2IXufm6qRumMuSd7Bd71GSGdNrj7CmBFGsfsuUNmm929sq/7KQSqpfCEUgdEtbQdORhMLSGdl1zXks6wTCMwPmG5HDjQXRszKwVGAJ8gIiJ5kU64bwImmdlEMxsI3Aqs69JmHfDH8e0FwL94qvEeERHJmZTDMu7eZmZLgTqiSyFXuvtOM3sQ2Ozu64CngGfMbB/RK/Zbc9lpsjC0U0BUS+EJpQ5QLYUq57WkfENVRESKjz6hKiISIIW7iEiACj7czWyQmb1tZu+Y2U4z+0G8fqKZvWVme81sVfxmb8EzsxIz22pmP4mXi7WO/Wa2w8y2mdnmeN0oM3slruUVM/utfPczHWY20sxeNLP/MLPdZnZlMdZiZpPj89Hxc9TM7irSWv57/Hh/18yej3OgWB8rfx7XsdPM7orX5fycFHy4AyeAa9z9cmAKMNfMriCa4uBv3H0S8CnRFAjF4M+B3QnLxVoHwCx3n5Jwve69wGtxLa/Fy8Xgh8BP3f0i4HKi81N0tbh7Q3w+phDN8/QbYA1FVouZnQ/cCVS6+yVEF3J0TGtSVI8VM7sE+A7RJ/0vB75uZpPoj3Pi7kXzAwwGfg7MIPp0V2m8/kqgLt/9S6P/5fGJvAb4CdGHv4qujriv+4ExXdY1AOfFt88DGvLdzzTqGA78gvjigmKupUv/ZwOvF2MtwPnAh8Aooiv6fgLMKcbHCnAz0WSLHct/CfyP/jgnxfDKvWMoYxtwEHgFeA847O5tcZNGov8Qhe5RohPbMQXHaIqzDog+gbzBzLbE00oAjHX3XwHE/34xb71L35eAj4B/jIfLnjSzIRRnLYluBZ6PbxdVLe7eBDwCfAD8CjgCbKE4HyvvAleb2WgzGwx8jegDnzk/J0UR7u7e7tGfmuVEf958OVmz/u1V75jZ14GD7p4490Va0zYUqJnuPg24HvgzM7s63x3KUCkwDXjc3acCn1HgwxapxGPR84AX8t2XTMTjzzcCE4FxwBCi/2ddFfxjxd13Ew0nvQL8FHgHaOtxoywpinDv4O6HgXrgCmBkPNUBJJ8SodDMBOaZ2X6imTWvIXolX2x1AODuB+J/DxKN61YB/2lm5wHE/x7MXw/T1gg0uvtb8fKLRGFfjLV0uB74ubv/Z7xcbLVcC/zC3T9y91bgJeAqivex8pS7T3P3q4k+5LmXfjgnBR/uZnaumY2Mb5cRnfjdwEaiqQ4gmvrg5fz0MD3uvtzdy919AtGfzP/i7rdRZHUAmNkQMxvWcZtofPddTp+GoihqcfdfAx+a2eR41VeBXRRhLQkW8fmQDBRfLR8AV5jZ4Hjq8I5zUnSPFQAz+2L87wXAfyU6Nzk/JwX/CVUzu4xorvgSoiej1e7+oJl9iegV8ChgK/BH7n4ifz1Nn5lVA/e4+9eLsY64z2vixVLgOXd/2MxGA6uBC4geoDe7e8FPIGdmU4AngYHA+8C3iP+vUXy1DCZ6M/JL7n4kXld05yW+5Hkh0RDGVmAx0Rh7UT1WAMzs34jeX2sF7nb31/rjnBR8uIuISO8V/LCMiIj0nsJdRCRACncRkQAp3EVEAqRwFxEJUDpfkC3Sr+LLxF6LF38baCeaIgCgyt1P5qVjPTCzPwHWx9fNi+SdLoWUgmZmDwDN7v5IAfSlxN3bu7nv34Gl7r6tF/srTZgrRSSrNCwjRcXM/tii+f23mdnfmdk5ZlZqZofNrMbMfm5mdWY2w8z+1czeN7OvxdsuNrM18f0NZnZfmvt9yMzeBqrM7Admtimen/vvLbKQaDrqVfH2A82sMeGT1VeY2avx7YfM7Akze4VosrJSM/s/8bG3m9ni/v+tSogU7lI04rmxbwKuiieSK+XzL2MfAWyIJzM7CTxA9LH1m4EHE3ZTFW8zDfiGmU1JY78/d/cqd38D+KG7Twcuje+b6+6rgG3AQo/mU081bDQVuMHdbweWEE0oVwVMJ5qE7YJMfj8iiTTmLsXkWqIA3BxNOUIZ0UftAVrc/ZX49g7giLu3mdkOYELCPurc/VMAM1sL/Beix0F3+z3J51MtAHzVzJYBg4AxRFPR/r9e1vGyux+Pb88GvmxmiU8mk4g+ki6SMYW7FBMDVrr7X562MpopMPHV8imib/DquJ34/7zrm0yeYr8tHr8xFc/b8rfANHdvMrOHiEI+mTY+/8u4a5vPutT0p+7+GiJZpGEZKSavAreY2RiIrqrJYAhjtkXfmTqYaM7w13ux3zKiJ4uP41kx/zDhvmPAsITl/URfdUeXdl3VAX/aMZWtRd+DWtbLmkTOoFfuUjTcfUc8W+CrZnYO0Sx7/43ezev978BzwO8Az3Rc3ZLOft39kJn9E9H0xr8E3kq4+x+BJ82shWhc/wHgH8zs18DbPfTnCaKZAbfFQ0IHiZ50RPpEl0LKWSO+EuUSd78r330RyTUNy4iIBEiv3EVEAqRX7iIiAVK4i4gESOEuIhIghbuISIAU7iIiAfr/4sNEDsfkXWAAAAAASUVORK5CYII=\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),\n",
" 'Intercept': 1})\n",
"data_pred['Frequency'] = logmodel.predict(data_pred)\n",
"print(data_pred.head())\n",
"\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": {},
"source": [
"La fonction `logmodel.predict(data_pred)` ne fonctionne pas avec les dernières versions de pandas (`Frequency = 1` pour toutes les températures).\n",
"\n",
"On peut alors utiliser le code suivant pour calculer les prédictions et tracer la courbe :"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Temperature Intercept Frequency Prob\n",
"0 30.0 1 1.0 0.834373\n",
"1 30.5 1 1.0 0.826230\n",
"2 31.0 1 1.0 0.817774\n",
"3 31.5 1 1.0 0.809002\n",
"4 32.0 1 1.0 0.799911\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3dyaTjSQsAcK+qAFEkCVhcamCWkUfBaqoKOBSlPZR21pbqvax7q3tQ/uzWrVIUar1EaRWkVorVRQVRQQEWWVHCPsaCNmT+/fHDBpDIJNkklnyeV1Xrsw5c59zvjcn+eRw5pz7mHMOERGJLZ5wFyAiIqGncBcRiUEKdxGRGKRwFxGJQQp3EZEYpHAXEYlB1Ya7mT1vZnvMbOUJ3jcze9LMNpjZcjPrH/oyRUSkJoI5cv8rMOwk718KZAa+JgB/rntZIiJSF9WGu3PuQ+DASZqMAF50fp8CzcysbagKFBGRmosLwTraA9sqTOcE5u2s3NDMJuA/uicpKSmrY8eOtdpgeXk5Hk9sfFygvkSeWOkHqC+Rqi59Wbdu3T7nXKvq2oUi3K2KeVWOaeCcmwJMAcjOznaLFy+u1QbnzZvHkCFDarVspFFfIk+s9APUl0hVl76Y2VfBtAvFn8EcoOIheAdgRwjWKyIitRSKcJ8N3BC4amYwkOucO+6UjIiINJxqT8uY2XRgCNDSzHKABwAfgHNuMvAWcBmwAcgHbq6vYkVEJDjVhrtz7rpq3nfA7SGrSESkGiUlJeTk5FBYWBjuUmqladOmrFmz5qRtEhMT6dChAz6fr1bbCMUHqiIiDSonJ4fU1FS6dOmCWVXXdES2I0eOkJqaesL3nXPs37+fnJwcunbtWqttxMZ1RSLSqBQWFpKenh6VwR4MMyM9Pb1O/zNRuItIVIrVYD+mrv1TuIuIxCCFu4hIDXm9Xvr27UuvXr24+uqryc/Pr9HybdvW/wgtCncRkRpKSkpi2bJlrFy5kvj4eCZPnvyt951zlJeXh6k6P4W7iEgdfOc732HDhg1s2bKF008/ndtuu43+/fuzbds2pk+fTu/evenVqxd33333t5b72c9+Rv/+/bnwwgvZu3dvyOvSpZAiEtUe+ucqVu84HNJ19myXxgNXnFFtu9LSUv79738zbJh/VPS1a9cybdo0nnnmGXbs2MHdd9/NkiVLaN68ORdffDGzZs1i5MiRHD16lP79+/OHP/yBhx9+mIceeoinnnoqpH3QkbuISA0VFBTQt29fsrOz6dSpE+PHjwegc+fODB48GIBFixYxZMgQWrVqRVxcHGPGjOHDDz8EwOPxcO211wIwduxY5s+fH/IadeQuIlEtmCPsUDt2zr2yJk2afP3af/N+cOrjsk4duYuI1INBgwbxwQcfsG/fPsrKypg+fTrnn38+4B/P/dVXXwXg5Zdf5txzzw359nXkLiJSD9q2bctjjz3G0KFDcc5x2WWXMWLECMB/hL9q1SqysrJo2rQpr7zySsi3r3AXEamhvLy84+Z16dKFlStXfmve9ddfz/XXX39c2507d5KamsojjzxSbzXqtIyISAxSuIuIxCCFu4hEpZpcjRKN6to/hbuIRJ3ExET2798fswF/bDz3xMTEWq9DH6iKSNTp0KEDOTk59XLbfkMoLCysNriPPYmpthTuIhJ1fD5frZ9QFAnmzZtHv3796nUbOi0jIhKDFO4iIjFI4S4iEoMU7iIiMUjhLiISgxTuIiIxSOEuIhKDFO4iIjFI4S4iEoMU7iIiMSjqwn3P4UI+yCmJ2QGDRERCIerC/aWFW5m2sphbXljM3iNF4S5HRCQiRV2433lhJtf3iOejDfsY9scPeXf17nCXJCIScaIu3D0e4+IuPv71o3PJSEvklhcX88AbKyksKQt3aSIiESPqwv2YzIxUXr/9bMaf25UXFnzFyKc/ZsOe4x9aKyLSGEVtuAMkxHn51eU9mXbzAPYcKWL4U/N5Y9n2cJclIhJ2QYW7mQ0zs7VmtsHM7qni/U5m9r6ZLTWz5WZ2WehLPbGh3Vvzrx+fyxnt0vjJjGX88vUVFJXqNI2INF7VhruZeYGngUuBnsB1ZtazUrP7gJnOuX7AaOCZUBdanbZNk5h+62B+eP6pvLxwK9dMXsD2QwUNXYaISEQI5sh9ILDBObfJOVcMzABGVGrjgLTA66bAjtCVGLw4r4d7Lu3Bs+Oy2LT3KJc/+RGfbNgXjlJERMLKqrsZyMxGAcOcc7cEpscBg5xzd1Ro0xb4D9AcaAJc5JxbUsW6JgATADIyMrJmzJhRq6Lz8vJISUk5aZtdR8t5cmkhu446ru0ez8Wd4zCzWm2vPgXTl2gRK32JlX6A+hKp6tKXoUOHLnHOZVfb0Dl30i/gamBqhelxwJ8qtbkL+Fng9VnAasBzsvVmZWW52nr//feDaneksMTd+sIi1/nuN91dryxzhSWltd5mfQm2L9EgVvoSK/1wTn2JVHXpC7DYVZPbzrmgTsvkAB0rTHfg+NMu44GZgT8WC4BEoGUQ665XKQlxTB6bxZ0XZfKPz3MY85eF7MvTXa0iEvuCCfdFQKaZdTWzePwfmM6u1GYrcCGAmZ2OP9z3hrLQ2vJ4jDsv6sbT1/dn5Y5cRjz1MWt3HQl3WSIi9aracHfOlQJ3AHOANfivilllZg+b2fBAs58Bt5rZF8B04KbAfx8ixn+d2ZaZPziLkrJyRv35Ez5aHxF/e0RE6kVQ17k7595yznVzzp3qnPt1YN79zrnZgdernXPnOOf6OOf6Ouf+U59F19aZHZox6/ZzaN88iZunLeKVRVvDXZKISL2I6jtUa6NdsyT+/sOzOPu0ltz9jxX88d11Gj5YRGJOowt3gNREH8/dmM2orA788d313PvaCkrLysNdlohIyMSFu4Bw8Xk9TBp1Ju2aJvLkexvYl1fEU9f3J9HnDXdpIiJ11iiP3I8xM+66uDuPjDiDuV/u4YbnPiO3oCTcZYmI1FmjDvdjxp3VhSdH92PptoNc++wCPeFJRKKewj3gij7teO7GAWzZf5Rrn9WgYyIS3RTuFZzXrRV/Gz+IvUeKuGbyArbsOxrukkREakXhXsmALi2YPmEwBSVlXPPsAj3dSUSiksK9Cr3aN2XGhMGUOxg9ZQFf7joc7pJERGpE4X4C3TJSmfmDwcR5PIye8imrduSGuyQRkaAp3E/ilFYpvPKDwST7vIyZupCV2xXwIhIdFO7V6JzehBkTzlLAi0hUUbgHoVN6MjMmnEWTeC9jn1vImp06By8ikU3hHqRO6clMnzCYxDgvY6cuZP1ujQkvIpFL4V4DndOb8PKtg/B6jOv+spBNe3WZpIhEJoV7DZ3SKoWXbx2Ec44xUxey7UB+uEsSETmOwr0WTmudyt/GDyK/uIwxUxeyK7cw3CWJiHyLwr2WerZL44XvD+TA0WLGPreQA0eLw12SiMjXFO510LdjM6bemM22A/nc+PxnHCnUcMEiEhkU7nU0+JR0/jy2P2t2HuaWFxZTWFIW7pJERBTuoXBBjwz+cE0fPttygDteXqpH9olI2CncQ2RE3/Y8eMUZvLtmN/e+tkIP3RaRsGq0z1CtDzee3YX9R4t5cu560lMSuOfSHuEuSUQaKYV7iP30okwOHC1i8gcbaZ2awPfP7RrukkSkEVK4h5iZ8dDwXuw7Uswj/1pNq9QErujTLtxliUgjo3Pu9cDrMf44ui8DOrfgrpnL+GTDvnCXJCKNjMK9niT6vPzlhmy6tmzCD/62RE9zEpEGpXCvR02TfUy7eSDJCV5unraInbkF4S5JRBoJhXs9a98siWk3DeRIYSk3T1uku1hFpEEo3BtAz3Zp/HlsfzbsyeO2//ucEt3kJCL1TOHeQL6T2YrffK83H63fx32vr9RNTiJSr3QpZAO6ZkBHth3M50/vbaBTejK3Dz0t3CWJSIxSuDewu77bja0H8pk0Zy2d05NJCXdBIhKTdFqmgZkZv7vqTLI7N+eumV+w4aBGkRSR0Asq3M1smJmtNbMNZnbPCdpcY2arzWyVmb0c2jJjS6LPy5QbsmnbNJEnlhbqUX0iEnLVhruZeYGngUuBnsB1ZtazUptM4F7gHOfcGcCd9VBrTGnRJJ7nbxpAWTmMf2ERh3WJpIiEUDBH7gOBDc65Tc65YmAGMKJSm1uBp51zBwGcc3tCW2ZsOrVVCnf0S2TT3qMaB15EQsqquyTPzEYBw5xztwSmxwGDnHN3VGgzC1gHnAN4gQedc29Xsa4JwASAjIyMrBkzZtSq6Ly8PFJSYuOjyLy8PJYcTGDaqmIu6hTH2J4J4S6p1mJlv8RKP0B9iVR16cvQoUOXOOeyq2sXzNUyVsW8yn8R4oBMYAjQAfjIzHo55w59ayHnpgBTALKzs92QIUOC2Pzx5s2bR22XjTTz5s3jgcuH4H1zNVPnb2ZI/x6MHdw53GXVSqzsl1jpB6gvkaoh+hLMaZkcoGOF6Q7AjiravOGcK3HObQbW4g97CdK9l53O0O6teGD2Ko0iKSJ1Fky4LwIyzayrmcUDo4HZldrMAoYCmFlLoBuwKZSFxjqvx3jyun6c2qoJP3xpCZv3HQ13SSISxaoNd+dcKXAHMAdYA8x0zq0ys4fNbHig2Rxgv5mtBt4HJjrn9tdX0bEqNdHH1BsG4PUY419YRG6BrqARkdoJ6jp359xbzrluzrlTnXO/Dsy73zk3O/DaOefucs71dM71ds7V7pNSoVN6Mn8em8XW/fn8aLquoBGR2tEdqhFo8CnpPDKyFx+u28tj//4y3OWISBTS2DIR6rqBnVi76wjPzd9MjzapXJ3dsfqFREQCdOQewe77r9M557R0/uf1lSz56mC4yxGRKKJwj2BxXg9PX9+fts0S+cHflugxfSISNIV7hGuWHM/UG7IpLCljwotLKCzRKJIiUj2FexTIzEjlj9f2ZeWOXO7+x3I9xUlEqqVwjxIX9czg5xd3541lO3j2Q90fJiInp3CPIrcNOZX/OrMtv3v7S+at1cCbInJiCvcoYmZMGnUmPdqk8aPpS9m0Ny/cJYlIhFK4R5nk+DimjMvC5/Vw64uL9ZAPEamSwj0KdWyRzNPX92fL/nzuemUZ5eX6gFVEvk3hHqXOOjWd+y/vybtr9vD4u+vCXY6IRBgNPxDFbjirM6t25PKn9zZwets0LuvdNtwliUiE0JF7FDMzHhnZi36dmvHzv3/Bl7sOh7skEYkQCvcolxDnZfLYLFIS4rj1xcUcyi8Od0kiEgEU7jEgIy2RyeOy2J1bxB0vawx4EVG4x4z+nZrz6MhezN+wj99qDHiRRk8fqMaQawZ0ZOWOXKbO38wZ7dP4Xr8O4S5JRMJER+4x5leX92RQ1xbc848VLM85FO5yRCRMFO4xxuf18MyY/rRMSeAHf1vC3iNF4S5JRMJA4R6D0lMSeHZcFgfzi7nt/5ZQXKoPWEUaG4V7jOrVvim/u+pMFm05yIP/XBXuckSkgekD1Rg2om97Vu88zLMfbOKMdmmMGdQ53CWJSAPRkXuM+8UlPTi/WyseeGMVn20+EO5yRKSBKNxjnNdjPHldPzq2SOa/X1rC9kN6yLZIY6BwbwSaJvn4yw3ZFJeWM+HFxRQU6yHbIrFO4d5InNY6hSeu68vqnYeZ+OoXesi2SIxTuDciF/TIYOIl3Xlz+U6embcx3OWISD1SuDcy/33+qVzRpx2//89a5q7ZHe5yRKSeKNwbGTPjf686kzPapfGTGctYv/tIuEsSkXqgcG+EkuK9TBmXTaLPyy0aA14kJincG6l2zZJ4dlwWOw8Vctv/fU6JxoAXiSkK90Ysq3NzfnNlbz7ZuJ9H3lwd7nJEJIQ0/EAjNyqrA+t2H2HKh5vIzEhl3GANUSASC3TkLtw9rAdDu7fiwdmr+GTDvnCXIyIhEFS4m9kwM1trZhvM7J6TtBtlZs7MskNXotS3Y0MUnNKyCf/9f5+zed/RcJckInVUbbibmRd4GrgU6AlcZ2Y9q2iXCvwYWBjqIqX+pSb6eO7GAXg9xvi/LiI3vyTcJYlIHQRz5D4Q2OCc2+ScKwZmACOqaPcI8L9AYQjrkwbUKT2ZyWOz2HYwn9tf1hU0ItHMqhtjxMxGAcOcc7cEpscBg5xzd1Ro0w+4zzl3lZnNA37unFtcxbomABMAMjIysmbMmFGrovPy8khJSanVspEmEvvyUU4Jz60sZkjHOG7sGY+ZBbVcJPalNmKlH6C+RKq69GXo0KFLnHPVnvoO5mqZqn6zv/6LYGYe4HHgpupW5JybAkwByM7OdkOGDAli88ebN28etV020kRiX4YA8W9/yZ/nbeQ7fbox/tyuQS0XiX2pjVjpB6gvkaoh+hLMaZkcoGOF6Q7AjgrTqUAvYJ6ZbQEGA7P1oWp0m3hxd4ad0YZH/7Wad1drDBqRaBNMuC8CMs2sq5nFA6OB2cfedM7lOudaOue6OOe6AJ8Cw6s6LSPRw+MxHr+2L73aNeXHM5aycntuuEsSkRqoNtydc6XAHcAcYA0w0zm3ysweNrPh9V2ghE9SvJfnbsymWZKP8S8sYmeunuIkEi2Cus7dOfeWc66bc+5U59yvA/Pud87NrqLtEB21x47WaYk8f/MAjhaVcfO0RRwp1CWSItFAd6hKtXq0SePpMf1ZvyeP219eqkskRaKAwl2Ccn63Vvx6ZC8+XLeX+99Yqcf0iUQ4DRwmQRs9sBPbDubz9Psb6dA8mduHnhbukkTkBBTuUiM/+253cg4WMGnOWto2TeTK/h3CXZKIVEHhLjXi8RiTRvVh75EifvHqclqnJnJuZstwlyUileicu9RYfJyHyeOyOK11Cj98aQmrdtT+GvhZS7dzzm/fo+s9/+Kc377HrKXbQ1ip1Dftv8ilcJdaSUv08debB5KWGMdN0xax7UB+jdcxa+l27n1tBdsPFeCA7YcKuPe1FQqIKKH9F9kU7lJrbZom8uL4gRSXlnPD859xuLhmV9BMmrOWgpKyb80rKClj0py1oSxT6on2X2RTuEudnNY6ledvymZnbgGPLy4kr6g06GV3HKr6jtcTzZfIov0X2RTuUmdZnVvwzJj+fHWknAkvLqaotKz6hYB2zZJqNF8ii/ZfZFO4S0hc0COD8b3i+WTjfu6csYyy8upP0Uy8pDtJPu+35iX5vEy8pHt9lSkhpP0X2RTuEjLntPfxq8t78u+Vu7j3teXV3sU6sl97HruyN+2bJWFA+2ZJPHZlb0b2a98wBUudaP9FNl3nLiE1/tyu5BaU8OTc9aQl+vif/zr9pE9yGtmvvcIgimn/RS6Fu4TcTy/K5HBBCVPnbyY10cdPLsoMd0kijY7CXULOzLj/8p7kFZXy+LvrSI73cut5p4S7LJFGReEu9cLjMX531ZkUlJTx67fWkOjzMO6sLuEuS6TRULhLvfF6jMev6UthcRm/emMV8XEerh3QKdxliTQKulpG6lV8nIenx/TnvG6tuOe1FfxjSU64SxJpFBTuUu8SfV6mjMvi7FPTmfjqF7yxTGOPiNQ3hbs0iESfl6k3DGBg1xb89JVlGlxKpJ4p3KXBJMV7ef6mAQzqms5dM5fx+lKdohGpLwp3aVDJ8XE8f9MABp+Szl0zv2Dm4m3hLkkkJincpcElxXt57sYBnHtaS37x6nJe+vSrcJckEnMU7hIWSfFe/nJDNhf2aM19s1by3PzN4S5JJKYo3CVsEn1e/jw2i0t7teGRN1fzxLvrqx1sTESCo3CXsIqP8/Cn6/pxVf8OPP7uOn79rzUKeJEQ0B2qEnZxXg+TRp1JamIcU+dvJreghMeu7E2cV8ceIrWlcJeI4PEYD1zRk6ZJPp6Yu56D+SU8dX0/Eis9DEJEgqNDI4kYZsZPv9uNh4afwdwvdzPuuYUcyi8Od1kiUUnhLhHnxrO78OTofnyxLZdRkxeQczA/3CWJRB2Fu0SkK/q044XvD2T34UKufOYTVm7PDXdJIlFF4S4R66xT03n1h2fj9RjXPLuA977cHe6SRKKGwl0iWvc2qcy6/Ry6tmzCLS8s5oVPtoS7JJGooHCXiJeRlsjMH5zFBT1a88DsVdz/xkpKy8rDXZZIRAsq3M1smJmtNbMNZnZPFe/fZWarzWy5mc01s86hL1UasyYJcTw7LpsJ553Ciwu+4qZpi8jNLwl3WSIRq9pwNzMv8DRwKdATuM7MelZqthTIds6dCbwK/G+oCxXxeoxfXnY6k0adycLN+xn+9HzW7T4S7rJEIlIwR+4DgQ3OuU3OuWJgBjCiYgPn3PvOuWPXq30KdAhtmSLfuDq7IzMmDCa/uIzvPf0xb6/cGe6SRCKOVTeOh5mNAoY5524JTI8DBjnn7jhB+6eAXc65R6t4bwIwASAjIyNrxowZtSo6Ly+PlJSUWi0badSX2jtYWM6flhaxKbecy7r6uCrTh9djdV6v9klkUl/8hg4dusQ5l11tQ+fcSb+Aq4GpFabHAX86Qdux+I/cE6pbb1ZWlqut999/v9bLRhr1pW4KS0rdva8td53vftONfnaB23uksM7r1D6JTOqLH7DYVZOvzrmgTsvkAB0rTHcAdlRuZGYXAf8DDHfOFQWxXpE6S4jz8pvv9eb3V/fh860HueyJj1iwcX+4yxIJu2DCfRGQaWZdzSweGA3MrtjAzPoBz+IP9j2hL1Pk5EZldWDW7eeQkhDHmKmf8uTc9ZSVa+hgabyqDXfnXClwBzAHWAPMdM6tMrOHzWx4oNkkIAX4u5ktM7PZJ1idSL05vW0as390LsP7tOP/vbOOMVM/ZWduQbjLEgmLoIb8dc69BbxVad79FV5fFOK6RGrl3dW7+WzzAQAWbjrAhX/4gNEDOjJn1W52HCqgXbMkJl7SnZH92od827OWbmfSnLX1vp1g3DdrBdMXbuPOXiWMv/ctrhvUkUdH9g5LLRIeGs9dYsaspdu597UVFJSUAeCAguIynv94y9dtth8q4N7XVgCENHgrb7u+thOM+2at4KVPt349Xebc19MK+MZDww9IzJg0Z+3X4XpMVWfdC0rKmDRnbb1vuz62E4zpC7fVaL7EJoW7xIwdh4I/v769Bm3rsu2a1BQqZSe4d+VE8yU2KdwlZrRrlhR0W6/HmL9+X71vuyY1hYrXqr6R60TzJTYp3CVmTLykO0mVnrnq8xg+77dDLd7roUWTeMY+t5Cf//0LDhyt+6P8qtp2ks/LxEu613ndNXXdoI41mi+xSR+oSsw49sFl5StWqpo3rFcbnpy7nikfbmLumt388rLTGZXVAavl0e2Jth2Oq2WOfWh67By710xXyzRCCneJKSP7ta8yUKua94thPRjetx33vb6Sia8uZ+bibTw0vFfItx0Oj47szaMjezNv3jw2jhkS7nIkDHRaRhq1Hm3SmPmDs/jdVb3ZuPcol//pI/62uohD+XU/VSMSTgp3afQ8HuPaAZ14/2dDGDu4M+9tLeX8SfOY9vFmSvTEJ4lSCneRgKbJPh4e0YuHz0miV/s0Hvrnai55/EPeXrnr2KinIlFD4S5SScdUDy+NH8TUG7LxeIwfvrSEUZMXsHCTRpuU6KFwF6mCmXFRzwze/sl3eOzK3mw7kM+1Uz7lhuc/Y3nOoXCXJ1IthbvIScR5PVw3sBMfTBzKLy/rwfKcQwx/6mPG/3WRQl4imsJdJAhJ8V4mnHcqH/1iKD+/uBuLvzrI8Kc+5qZpn7Foy4FwlydyHIW7SA2kJvq444JM5t89lImXdGdFTi5XT17ANZMX8O7q3ZTrASESIRTuIrWQmujj9qGnMf/uC7j/8p5sP1TALS8u5ruPf8DLC7dSUFxW/UpE6pHCXaQOkuK9fP/crsybOIQnRvcl0efll6+v4KzfzuW3//6SbQfyw12iNFIafkAkBHxeDyP6tmd4n3Ys2nKQaR9vZsqHG3n2w41c0L01YwZ34vxurfF6NDKjNAyFu0gImRkDu7ZgYNcW7DhUwPTPtjL9s23M/eti2jVN5OrsjozK6kDHFsnhLlVinMJdpJ60a5bEzy7uzo8uyGTumt28/NlWnnxvPU/MXc9Zp6RzVVYHhvVqQ0qCfg0l9PRTJVLP4uM8XNq7LZf2bkvOwXxe+3w7ry7J4ed//4L7Zq3guz3bMLxPO87r1pKEOG/1KxQJgsJdpAF1aJ7Mjy/M5EcXnMbnWw/y+tLtvLl8J//8YgepiXFc3LMNl/Zqw7mZLUn0Keil9hTuImFgZmR1bkFW5xY8cMUZfLxhH//8YifvrN7FPz7PISUhjvO7t+LinhkM6d6apkm+cJcsUUbhLhJmPq+HId1bM6R7a4pLe/PJxn28vXIX767Zw7+W78TrMQZ0ac4FPfxtMlun1PqJUdJ4KNxFIkh83DdBX17uWLrtEO99uZu5a/bwm7e+5DdvfUnbpol8J7Ml52a24uxT02mZkhDusiUCKdxFIpTHY2R1bk5W5+ZMvKQHOw4V8OG6vXywbi9vr9zFzMU5APRok8rgU9IZfEo6A7o0J11hLyjcRaJGu2ZJjB7YidEDO1FW7lixPZePN+xjwcb9zFi0lb9+sgWA01qnkB34o1B2tBznnE7jNEIKd5Eo5PUYfTs2o2/HZtw+9DSKSstYkZPLZ1sOsGjzAd5asZMZi7YB8Nsl79CnQzP6dGxGnw5N6d2hKa1TE8PcA6lvCneRGJAQ5yW7Swuyu7SAIVBe7ti4N4+X//MpBckZLNt2iKfeW8+xQStbpybQq31TerZNo2e7NE5vm0anFskaHiGGKNxFYpDHY2RmpHJ+Rx9DhpwJQH5xKat2HGZ5Ti6rduSycnsuH6zbS1kg8RN9HjJbp9ItI5VuGSlkZqRwWqtU2jdPUuhHIYW7SCORHB/HgC4tGNClxdfzCkvK2LAnj9U7D7N21xHW7jrCR+v38o/Pc75ukxDnoWvLJnRt2YQuLZvQNb0JndOT6dKyCa1SEvAo+COSwl2kEUv0eenVvim92jf91vxD+cVs2JPHxr15bNiTx+Z9R1m7+wjvrN5NaYUHkiTEeejYIpkOzZPo2DyZ9s2TaN8s6evvLVMSdNQfJgp3ETlOs+T4b87hV1BaVs6OQ4Vs3n+UrQfy2XYgn6/2HyXnYAFLtx4it6DkW+3jPEZFdqscAAALRElEQVRGWiJtmibSJi2RjLREMtISaJ2WQOvURFqnJtAyJYFmyT5d0RNiCncRCVqc10On9GQ6pVc9ZPGRwhK2Hypg+8ECduQWsvNQAbtyC9l1uJA1Ow/z/to95FfxlCqf10hvkkB6SjzpKQmkN4mnRYWv5snxfHWgjLa7jtA82Udakk9j71RD4S4iIZOa6KNHGx892qSdsE1eUSm7Dxey53ARe44Usi+vmH15Rew9UsSBo/7Xm/bmceBo8XF/CB777MOvXyf6PDRN8pGW6PN/T/KRlhhHaqKP1MD3lMQ4UhPiaJIQR0rgq0mCl5SEOJIT4kj2eWP2M4Ogwt3MhgFPAF5gqnPut5XeTwBeBLKA/cC1zrktoS1VJHbNWrqdSXPWsuNQAe2aJTHxku78ffFWPt544Os255zagquzOx3XDjhu3uKvDjB94Tbu7FXC+Hvf4rpBHXl0ZO+gtlvV+kb2ax903ce2XeYcXrPjtp2SEEdKqxRW5ORWu+0Hr8jkO91acuBoMR8sWEynzNM5mF/Cpxv388G6vew+XBQ4FZRMQUkZG/aUcriwhCOFpV9fBVSdJJ+X5HgvSfHHvseR5POQHB9Hks9Lgs9Dks9Los9Los9DYtw3rxPi/O8nxAVex3lI8HmI93qJj/N88+X95rvPazhX/w9SrzbczcwLPA18F8gBFpnZbOfc6grNxgMHnXOnmdlo4HfAtfVRsEismbV0O/e+toKCEv9R6vZDBdz5yrLj2n288cC3wn77oQImvvoFOCgJBNn2QwXc9coyyissV+YcL326FeBbIVvVdif+/QswKCn7Zn33vrYC4LiAr2r5UG/7gdmreOzK3ozs15696V6GnNmOWUu3896Xe75etrCknJyDBV+3A3DOUVhSzpGiEvIKS8krCnwVlpJfXMbR4lKOFpVytKiMo0Wl5JeUUVBcRn5xKQUl5RQUl7L3SBEFgfmFJWUUlPi/B/k346TG9YxnaN1Xc1LBHLkPBDY45zYBmNkMYARQMdxHAA8GXr8KPGVm5hriz5NIlJs0Z+3XQVVTx4KwovIq2gFMX7jtWwFb1XZLqkiugpIyJs1Ze1y4V7V8Q2y7qmUrtzMzkgJH461TT1BULTjnKClzFJaWUVRSTmFJGcVl5RSVlFNUWkZxaTlFpeUUl5b755eWUVLqKCrzzyspK6ektJzUvK2hK+oErLr8NbNRwDDn3C2B6XHAIOfcHRXarAy0yQlMbwy02VdpXROACYHJ7sDaWtbdEthXbavooL5EngbtR3yb07Lqa91l+bl4k7+5zLF414Yltd1uxWXrunwtl20J7DvZspVrjGB1+Rnr7JxrVV2jYI7cq/q0ofJfhGDa4JybAkwJYpsnL8hssXMuu67riQTqS+SJlX6Avy+luXtipi+xtF/quy+eINrkAB0rTHcAdpyojZnFAU2BA4iISFgEE+6LgEwz62pm8cBoYHalNrOBGwOvRwHv6Xy7iEj4VHtaxjlXamZ3AHPwXwr5vHNulZk9DCx2zs0GngP+ZmYb8B+xj67PognBqZ0Ior5EnljpB6gvkare+1LtB6oiIhJ9gjktIyIiUUbhLiISgyI+3M0s0cw+M7MvzGyVmT0UmN/VzBaa2XozeyXwYW/EMzOvmS01szcD09Hajy1mtsLMlpnZ4sC8Fmb2TqAv75hZ83DXGQwza2Zmr5rZl2a2xszOisa+mFn3wP449nXYzO6M0r78NPD7vtLMpgdyIFp/V34S6McqM7szMK/e90nEhztQBFzgnOsD9AWGmdlg/EMcPO6cywQO4h8CIRr8BFhTYTpa+wEw1DnXt8L1uvcAcwN9mRuYjgZPAG8753oAffDvn6jri3NubWB/9MU/zlM+8DpR1hczaw/8GMh2zvXCfyHHsWFNoup3xcx6Abfiv9O/D3C5mWXSEPvEORc1X0Ay8DkwCP/dXXGB+WcBc8JdXxD1dwjsyAuAN/Hf/BV1/QjUugVoWWneWqBt4HVbYG246wyiH2nAZgIXF0RzXyrVfzHwcTT2BWgPbANa4L+i703gkmj8XQGuxj/Y4rHpXwG/aIh9Eg1H7sdOZSwD9gDvABuBQ8650kCTHPw/EJHuj/h37LEhONKJzn6A/w7k/5jZksCwEgAZzrmdAIHvrcNWXfBOAfYC0wKny6aaWROisy8VjQamB15HVV+cc9uB3wNbgZ1ALrCE6PxdWQmcZ2bpZpYMXIb/hs963ydREe7OuTLn/69mB/z/vTm9qmYNW1XNmNnlwB7nXMWxL4IatiFCneOc6w9cCtxuZueFu6BaigP6A392zvUDjhLhpy2qEzgXPRz4e7hrqY3A+ecRQFegHdAE/89ZZRH/u+KcW4P/dNI7wNvAF0DpSRcKkagI92Occ4eAecBgoFlgqAOoekiESHMOMNzMtgAz8J+a+SPR1w8AnHM7At/34D+vOxDYbWZtAQLf94SvwqDlADnOuYWB6Vfxh3009uWYS4HPnXO7A9PR1peLgM3Oub3OuRLgNeBsovd35TnnXH/n3Hn4b/JcTwPsk4gPdzNrZWbNAq+T8O/4NcD7+Ic6AP/QB2+Ep8LgOOfudc51cM51wf9f5vecc2OIsn4AmFkTM0s99hr/+d2VfHsYiqjoi3NuF7DNzLoHZl2IfzjrqOtLBdfxzSkZiL6+bAUGm1mymRnf7JOo+10BMLPWge+dgCvx75t63ycRf4eqmZ0JvID/E3MPMNM597CZnYL/CLgFsBQY65wrCl+lwTOzIcDPnXOXR2M/AjW/HpiMA152zv3azNKBmUAn/L+gVzvnIn4AOTPrC0wF4oFNwM0EftaIvr4k4/8w8hTnXG5gXtTtl8Alz9fiP4WxFLgF/zn2qPpdATCzj/B/vlYC3OWcm9sQ+yTiw11ERGou4k/LiIhIzSncRURikMJdRCQGKdxFRGKQwl1EJAYF84BskQYVuExsbmCyDVCGf4gAgIHOueKwFHYSZvZ94K3AdfMiYadLISWimdmDQJ5z7vcRUIvXOVd2gvfmA3c455bVYH1xFcZKEQkpnZaRqGJmN5p/fP9lZvaMmXnMLM7MDpnZJDP73MzmmNkgM/vAzDaZ2WWBZW8xs9cD7681s/uCXO+jZvYZMNDMHjKzRYHxuSeb37X4h6N+JbB8vJnlVLizerCZvRt4/aiZPWtm7+AfrCzOzP5fYNvLzeyWhv9XlVikcJeoERgb+3vA2YGB5OL45mHsTYH/BAYzKwYexH/b+tXAwxVWMzCwTH/gejPrG8R6P3fODXTOLQCecM4NAHoH3hvmnHsFWAZc6/zjqVd32qgfcIVzbhwwAf+AcgOBAfgHYetUm38fkYp0zl2iyUX4A3Cxf8gRkvDfag9Q4Jx7J/B6BZDrnCs1sxVAlwrrmOOcOwhgZrOAc/H/HpxovcV8M9QCwIVmNhFIBFriH4r23zXsxxvOucLA64uB082s4h+TTPy3pIvUmsJdookBzzvnfvWtmf6RAiseLZfjf4LXsdcVf84rf8jkqllvgQt8MBUYt+UpoL9zbruZPYo/5KtSyjf/M67c5milPt3mnJuLSAjptIxEk3eBa8ysJfivqqnFKYyLzf/M1GT8Y4Z/XIP1JuH/Y7EvMCrmVRXeOwKkVpjegv9Rd1RqV9kc4LZjQ9ma/zmoSTXsk8hxdOQuUcM5tyIwWuC7ZubBP8reD6nZuN7zgZeBU4G/Hbu6JZj1Ouf2m9kL+Ic3/gpYWOHtacBUMyvAf17/QeAvZrYL+Owk9TyLf2TAZYFTQnvw/9ERqRNdCimNRuBKlF7OuTvDXYtIfdNpGRGRGKQjdxGRGKQjdxGRGKRwFxGJQQp3EZEYpHAXEYlBCncRkRj0/wEPe112mMXbNAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Inspiring from http://blog.yhat.com/posts/logistic-regression-and-python.html\n",
"def logit_inv(x):\n",
" return(np.exp(x)/(np.exp(x)+1))\n",
"\n",
"data_pred['Prob']=logit_inv(data_pred['Temperature'] * logmodel.params['Temperature'] + \n",
" logmodel.params['Intercept'])\n",
"print(data_pred.head())\n",
"\n",
"data_pred.plot(x=\"Temperature\",y=\"Prob\",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": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"c:\\program files\\python\\python37\\lib\\site-packages\\scipy\\stats\\stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9+P/XObMlk4TsIciOAmFXcQFkqVYIO4q0dflKXUpre1taf/f2eit2u7W1tt5ra6+1UlutCi3WVlnUgEqpsmgRlT3syBZCFrLOdpbP749JBiIQk5DJZCbv5+MRkpk558znQzLzns/2/mhKKYUQQghxAXqsCyCEEKJzk0AhhBCiWRIohBBCNEsChRBCiGZJoBBCCNEsCRRCCCGaJYFCCCFEsyRQCCGEaFbUA0VdXR0zZ87k2LFj5zy2e/du5s6dS2FhIYsWLcI0zWgXRwghRCtFNVBs3bqV2267jcOHD5/38e9+97v84Ac/YPXq1SileOmll6JZHCGEEG0Q1UDx0ksv8cMf/pC8vLxzHjt+/DiBQIDLL78cgLlz51JUVBTN4gghhGgDZzQv/tOf/vSCj506dYrc3NzI7dzcXEpLS6NZHCGEEG0Qs8Fs27bRNC1yWynV5LYQQojOIaotiubk5+dTVlYWuV1eXn7eLqrmHDpaiWHakduNYUbTwjc0TUNHAw10TUNDQ9PD9zfepzecqGsN52ta+LyGi8Uqt252dioVFXWxefIOIPWLX4lcN0js+um6RmZmSqvPi1mg6NmzJx6Phy1btjB69GiWL1/OxIkTW3WNkGERMuzPPvAzaJF/iLRqdA10XcehaegODYeu4XCEg41DB4euo2nRDSS2ndgZ4KV+8SuR6waJX7/W6vBAsWDBAhYuXMiIESN47LHHeOihh6irq2PYsGHMnz+/o4sDgIr8E+4CA7ABLOucYxuDiqZp6Do4dR2nU8fp0NE1DWdDUJFdPoQQiUKL542L9hwsa5cWRXsJ92iFA4jL4cDlCgeQtgSP3Nw0yspqo1fYGJP6xa9Erhskdv10XSM7O7XV58Ws6ykRKQUKhW2BaZn4Q+EWiNbQbeVxOXA7HLhcGromrQ4hRHyQQBFlClC2wrYVhmEDBrqu4XboeDwO3A3dVhI0hBCdlQSKGLBtRcC2CBgWuqbhcuokeRx4XA50mSIshOhkJFDEmK0UQcMiaFjoerh7KtntlFkXQohOQwJFJ2LbCn/QJBA0cXgcBAIGSS4HTqd0TQkhYkcCRSekAMuGWp9BvWbidul4PS5cLh3pmBJCdDQJFJ2crRSBkEUwZOFwaqQkuUh2y69NCNFx5B0nTijANBXVdSHqdQNvsotkjwNN2hhCiCiTHe7ikGkraupDVNQECBkWMlFKCBFNEijimGkqTtcFOV0XxLY7zwp1IURikUAR55SCQNCioiZIUFoXQogokECRICxbUVUXpMZnxLooQogEI4EigSgF9X6DylrpihJCtB8JFAkoZFhU1gQJmRIshBAXTwJFgjJtRVVtkIBx7p4aQgjRGhIoEpitFNV1QXxBUwa5hRBtJoEiwSkFtb4QdQEJFkKItpFA0QUoBXX1Ier8hgQLIUSrSaDoIhRQ5zMkWAghWk0CRReigDq/QX3QjHVRhBBxRAJFF6MU1NaHCMpsKCFEC0mg6IKUguq6kKyzEEK0iASKLqpx6qwlK7iFEJ9BAkUXZtnh/S0Uss+qEOLCJFB0cSHTpqbeAAkWQogLkEAh8AdN6gMyE0oIcX4SKAQQnjYrM6GEEOcjgUIADTOh6kOYlnRBCSGakkAhImxbUVUfwFYSLIQQZ0igEE2YpqKmPiRD20KIiLgOFK++c5Ca+lCsi5FwAiGLOr9sqSqECIvrQLH3WDWPv7SV93eVSndJO/MFDEIyuC2EIM4DhdfjJGhYLF9/iGdW7aKiJhDrIiUMpaCmPoRlSwAWoquL60Bxz4whXDEwB4DDJbU88fI2Nu44Ka2LdmLailqfdO0J0dXFdaDwJjn5wvWX8eWpg+mW4sYwbVZtPMwfVu3idK20LtpDIGRJWnIhurioBoqVK1cyffp0pkyZwpIlS855fOfOndxyyy3Mnj2br33ta9TU1LTpeQb3yeTb80YyelAuAIdKanni5e18uLcMJa2Li1bvM2R9hRBdWNQCRWlpKY8//jhLly7l1VdfZdmyZezfv7/JMT/96U9ZuHAhK1asoH///vzhD39o8/Mle5zc8rlLmT91MKnJLoKGxcvrDrD0rX34AjKD52LYSlHjCyExV4iuKWqBYuPGjYwZM4aMjAy8Xi+FhYUUFRU1Oca2berr6wHw+/0kJSVd9PMW9Mlk4byRDO2XCcDOQ5X85m/bOXiiba0VERYyLHxBCbhCdEWailLfzNNPP43P5+P+++8H4K9//Svbtm3jJz/5SeSYjz/+mHvuuQev10tycjIvvfQSmZmZLX6OkxX1F5yVo5Ri0/YSlr25l6BhoWkwfVx/pl/XD4ce10MzMaMB2d2S8HicsS6KEKIDRe0Vb9s2mqZFbiulmtwOBAIsWrSI5557jpEjR/Lss8/ywAMPsHjx4hY/R3W1j5Bx4Y13Cnql8283D+cvb+/jRIWP1zYcYseBcm79/EDSU9xtq1gHycpKobKyPtbFOEdtbYDMNA/aZx/arNzcNMrKatulTJ1RItcvkesGiV0/XdfIzk5t/XlRKAsA+fn5lJWVRW6XlZWRl5cXub137148Hg8jR44E4Etf+hL/+te/2r0cORnJ3HfTcMaP6AHAJydr+b+/bWPfsap2f66uINwFJbOghOhKohYoxo0bx6ZNm6isrMTv97NmzRomTpwYebxv376cPHmSgwcPAvD2228zYsSIqJTF6dCZPrYv8wsHk+xxUB8wee71Yt764Ci2LChrtXq/gS1bqArRZUSt66l79+7cf//9zJ8/H8MwmDdvHiNHjmTBggUsXLiQESNG8Mgjj/Cd73wHpRTZ2dn87Gc/i1ZxACjom8k3547kz2/t5VhZPWs/PM7RU3V86YaBeJOk372lbFtR4zfITPXITCghuoCoDWZ3hD0Hy5odo7gQ07J5470jbNp5EoDMNA93TB7EJTkp7V3ENuusYxSNNCAjzYPH5WjT+YncDwyJXb9Erhskdv063RhFZ+Z06My6rh9fuP5SXA6d07VBfrd8Bx/vK4910eKGAmp9IZQkJBci4XXJQNHoioG53HfTMDLTPJiW4qV/7Kfo/U9k3KKFTEvJXttCdAFdOlAA9MhO4d9uHsFlPdMBeGdrCS+s3kMgJG+ALeHzm1gysC1EQuvygQLCyQW/PK2A64bnA7DnaBVPvbpD0pa3gK0UdX4T7WIXVgghOi0JFA0cusaMcf24ZdIAHLpGWVWAp17ZweGTkvrjswRCJsE2TCoQQsQHCRSfMnpwHvfOHII3yYkvaPKHVbv5aF/ZZ5/YhSmFbJ0qRAKTQHEe/fK78Y2bhpObkYRlK/76jwO8veWYpCxvRsiwCMjWqUIkJAkUF5DVLYn75gyPDHK/veUYf//nQRm4bUadX1KRC5GIJFA0I9nj5MvTBnNlw4ZIW/aW8ac3ZEbUhZimImhKq0KIRCOB4jM4dJ1bJg3ghit7ArD/eDW/X7lL9pK+gEBAZkAJkWgkULSApmnceFVvbpk0AF2Dkgofv1u+k4pqmT77aUHTwrSke06IRCKBohVGD87j/xUObpL243hZXayL1akoBYGQdD8JkUgkULRSQZ9M7p05hGSPk/qAye9X7eLAiepYF6tT8QdMGdQWIoFIoGiDPt3T+NrsYaSnuAkZNn96o5jdhytjXaxOw7QVhgxqC5EwJFC0UV5mMl+bM4yc9CRMS7Hkzb2yMO8svqAMaguRKCRQXISMVA9fnT2MS7K92Ar++o8DvNewx0VXFzJsGdQWIkFIoLhIqckuvjJrKP3y0wBYseEw7249EeNSxZ6tFAHJ/yREQpBA0Q6S3E7uml4QWcX9xvtHJOUH4A8asq2REAlAAkU7cTsd3Fk4mII+mUA45cfqfx3t0sHCMhWGKa0KIeKdBIp25HLq3DFlICMGZAHwztYTvPH+kS4bLBRIuhMhEoAEinbm0HW+dMNArhiYA8D6bSW8tumTLhssAiELu4vWXYhEIYEiCnRd45ZJl0aSCW7ccZKVGw53yWBh24qgpB8XIq5JoIgSXdeYO2kAVxXkAfDertIuGyz8QQtkWFuIuCWBIop0TeOmCf25ZsiZYLFqY9frhjJMC8PqWnUWIpFIoIgyXdOYPb4/Vze0LDbtPNnlxiwkUaAQ8U0CRQfQNY05E/pz1eAzYxZvvNe1ZkMFgiZKup+EiEsSKDqIrmncNHEAoxuCxfrtJazZ3HXWWVi2IigrtYWISxIoOpCuadw8YUBk6uw/Pz7B2g+Px7hUHccnu98JEZckUHSw8GyoSyOL8t7ecox/ftw1goVhWoRkpbYQcUcCRQw4dI0v3nAZQ/uF032s/tdRNu5I/KyzMqgtRHySQBEjDl3n1s8PZFDvDABWbTzMB8WnYlyq6PMHTVmpLUSckUARQ06Hzh2TB9G/RzcAXnnnIFv3l8e4VNElK7WFiD8SKGLM5dSZXziY3nmpKMKbH+3+5HSsixVVvoDZZWZ7CZEIWhQoXnjhBerq6qJdli7L43Zw17QCemR7sZXiz2/tZU8CBwvTsiWrrBBxpEWBYs+ePRQWFrJo0SK2b9/e4ouvXLmS6dOnM2XKFJYsWXLO4wcPHuTOO+9k9uzZ3HvvvVRXV7e85Akm2ePk7ulDIntw//ZvWzlWlpjBWSkIBKX7SYh40aJA8fDDD7N69WqGDx/Oj3/8Y2655RZefvllgsHgBc8pLS3l8ccfZ+nSpbz66qssW7aM/fv3Rx5XSvH1r3+dBQsWsGLFCoYMGcLixYsvvkZxLDXZxT0zhpCe4iYYsnju9WJKT/tiXayoCIRMWactRJxo8RhFamoqU6dOZebMmVRVVbF06VKmTp3K2rVrz3v8xo0bGTNmDBkZGXi9XgoLCykqKoo8vnPnTrxeLxMnTgTgvvvu44477rjI6sS/jFQP98wYQprXhS9o8uxruzldG4h1sdqdrZDd74SIEy0KFJs2beI73/kOU6dO5eDBgzz55JP8/e9/509/+hM/+MEPznvOqVOnyM3NjdzOy8ujtLQ0cvvIkSPk5OTw4IMPcvPNN/PDH/4Qr9d7kdVJDLkZyXzri1fgcTmo8Rn88fVi6vxGrIvV7mT2kxDxwdmSg3784x9z++2385Of/IS0tLTI/X369OGLX/ziec+xbRvtrHwNSqkmt03T5F//+hcvvvgiI0aM4Fe/+hU///nP+fnPf97iwqene7HsxOzAyAK++YVRPPHSx1RUB3hhzV7+v9uvJNnTol9ZXEj2usnKSsGhJ2Zej9zctM8+KE4lct0g8evXWi1611mxYgVFRUWkpaVRVlbGa6+9xvz589F1nYULF573nPz8fD744IPI7bKyMvLy8iK3c3Nz6du3LyNGjABg5syZF7zWhVRX+wglaKK5rKwUslPd3HrDZSx5cy9HS2t54i8fcte0Ibic8T+rOSsrhfKKeuyQhdsV//X5tNzcNMrKamNdjKhI5LpBYtdP1zWys1Nbf15LDvrJT37CunXrGp5IZ8uWLfzsZz9r9pxx48axadMmKisr8fv9rFmzJjIeAXDFFVdQWVlJcXExAGvXrmXYsGGtrkCiG9Ivi7mTLgXgUEkty9buS6hWVMCQabJCdHYtalF89NFHrFq1CoDs7Gx+/etfM2fOnGbP6d69O/fffz/z58/HMAzmzZvHyJEjWbBgAQsXLmTEiBE8+eSTPPTQQ/j9fvLz8/nFL35x8TVKQFcOysUXMHn9vU/Ydfg0y9cf4uYJ/Zt05cWrQMgizavQiP+6CJGoWhQoDMMgFArhdruB8PhCS8yaNYtZs2Y1ue/3v/995OdRo0bx8ssvt7SsXdr4kT2o8xu8s/UEHxSfIjXZxZSre8e6WBfNthUh08bjdMS6KEKIC2hRoPjc5z7Hvffey5w5c9A0jVWrVjFp0qRol018SuE1van3G2zZW8a6j46Tmuxi3PD8WBfrovkDJklpDiSrhxCdU4sCxX/+53+yZMkS3n77bZxOJ5MnT+bWW2+NdtnEp2gNu+TVB0yKj5zmtY2HSU12MfLS7FgX7aKEDBvTsnHoiTeoLUQi0FQcZ2fbc7AsoWc9VVbWn/exkGnxx9d2c6S0Doeucde0Ai7tmd7BJbw4n65fWoqblASa+pvIM2cSuW6Q2PWL6qynt956ixtuuIHRo0dz5ZVXRr5EbLidDuYXFpCXmYxlK15cs5cT5ecPKvHCHzQkpYcQnVSLPsL98pe/5L/+678YOnRoQsy0SQTeJCd3TSvg6eU7qa4P8ac3ivnanGFkdUuKddHaxDIVhmEn5JoKIeJdi16V3bp1Y8qUKfTq1YuePXtGvkRsZaR6uGtaAUluB7V+g2ffKKY+EJ+pPhTgl9TjQnRKLQoUo0aN4p///Ge0yyLaoHuWl/lTB+N0aFRUB3i+aA+hOM2hFDSshFpMKESiaFHX0z//+U9efPFFXC4XLpcrkrfpww8/jHb5RAv0y+/GF28YyJ/f3MvRU3X8+a19/L/CwXGXQ6lxm1RvAg1qC5EIWvSKfO6556JcDHGxhvfPYtb4fqxYf5g9R6tY/u5Bbp44IO7GlHxBQwKFEJ1Mi7qeevbsyfbt23nppZfIysrio48+kjGKTmjM0Hw+d0X49/LBnjLe3nIsxiVqPctUCTvlWYh41aJAsXjxYv785z9TVFREIBDg//7v/3jyySejXTbRBpOv6sWVg3IAWPvhcTbvLv2MMzoXRbhVEWcNISESWosCxWuvvcbvf/97kpOTyczM5KWXXookCRSdi6Zp3DxxAAN7hRfgvbr+EMWfnI5xqVqncaW2EKJzaFGgcDqdkYSAEJ4u63RKP3Jn5dB1bp88iJ45KSgFf35rH0dPxc9KU1spAqH4nLklRCJqUaDo0aMH69atQ9M0QqEQTz31lIxRdHIel4P5UweTlebBsGz+VLSH8mp/rIvVYv6ArKkQorNoUaD4/ve/z7PPPsuePXu4/PLLeeedd/j+978f7bKJi5TmdXPXtAK8Hie+gMlzrxdT6wvFulgtYtqKQJyuBxEi0bQqKaDf78eyLFJTW59UKhq6alLA1jpSWssfVu3GsGx65qbwlZlD8bhiu/9DS+rncTnI6uaJy/TjiZxYLpHrBoldv7YmBWzRQMOzzz573vvvvvvuVj+h6Hh9uqdx640DeXHNHo6X1fPnt/ZxZ+GgTp/WO2RaGKbC6ZApUELEUosCxd69eyM/h0IhNm/ezNixY6NWKNH+hvTNZPZ1/Vm+/hB7j1bx6ruHmNvJF+QpFc7/lJbsinVRhOjSWhQoHnnkkSa3S0tLWbRoUVQKJKLn2qHdqakP8Y+PjrNlTxnpKW5uvKpzb6caDJqkJjtlT20hYqhNfQ/du3fn+PHj7V0W0QFuvKoXVw7KBRoW5BWfinGJmmfaimCCjkMJES9aPUahlGLHjh1kZ8f39ptdVXhBXn9qfSH2Hatm+bsHSUt2UdA3M9ZFu6BAwCTZLXtqCxErLWpR7N27N/K1b98+evTowWOPPRbtsokoOXtBnh0HC/KCpiUrtYWIIdkzu5Nqz+mxF1LrC/G75Ts5XRvEm+TkvjnDyElPjupzNmpt/dK8LlKS4mdQO5GnWCZy3SCx6xfV6bF33nlns7Njnn/++VY/sYi9NK+bu6cX8LvlO/EFTJ59vZj75gwjzev+7JM7mD9o4k1yyZC2EDHQoq6n4cOHk5SUxPz587n33nvJyckhIyODO+64gzvuuCPaZRRRlJOezJenDsbl0DldG+RPRXsIdsI8S6alMGSlthAx0aIWxYcffsjSpUtxOMKreSdMmMAXv/hFCgsLo1o40TF656VxW8OCvBPl9Sx9ay93Fg7G6ehcC/J8QROPDGoL0eFa9E5QWVlJMBiM3K6vrycQCEStUKLjFfTNZM6EAQDsO1bNK+8cpLMNX0n6cSFio0UtipkzZ/KlL32JyZMno5TijTfeYP78+dEum+hgVxfkUVMf4u0tx/hoXzlpXhdTr+0b62JF2EoRMGxSPJ2rpSNEomtRoPj2t7/N0KFDee+99/B4PPz3f/8311xzTbTLJmLghit7UusL8a/dp3hnawlpXjfXjegR62JF+APhPbVlUFuIjtPij2bdu3dn4MCBfOc738Hlip9piqJ1NE1j9nX9GdovvADvtU2fsHV/eYxLdYYMagvR8VoUKP72t7/xve99j2eeeYba2lq+8Y1v8NJLL0W7bCJGdF3jSzcMpG9+GgAvrzvA/mPVMS7VGb6gKXtqC9GBWhQoXnzxRZYtW0ZqairZ2dn8/e9/509/+lO0yyZiyOXUmV84mLzMZCxb8eKbezheVhfrYgEyqC1ER2tRoNB1vclmRT169IhMlRWJK9nj5O5pBaSnuAkZNs+9UdwptlOVPbWF6FgtChQZGRns3r07sjp7xYoVpKenR7VgonNIT/Vw9/QhJHuc1Des3q7pBNup+gOmrKcQooO0KFA8+OCDfPe73+XAgQOMHz+eX//61zz00EPRLpvoJPIyk7lr2mBczvDq7edeL8YfNGNaJtNWhExpVQjREVoUKAKBAMuXL+eVV17hj3/8I0VFRQwePPgzz1u5ciXTp09nypQpLFmy5ILHrVu3jhtuuKHlpRYdrndeGndMHoSuaZys9PH86j0xf6P2BWRQW4iO0KJA8R//8R84HA4uvfRSBg0a1KLpsaWlpTz++OMsXbqUV199lWXLlrF///5zjisvL+fRRx9tfclFhxvUO4N5118KwCcna/nLW/uw7NgNKodMi5Apg9pCRFuLAsXgwYNZuXIlJ06coKqqKvLVnI0bNzJmzBgyMjLwer0UFhZSVFR0znEPPfQQ3/zmN9tWetHhLr8sh5nj+gFQfKSKv607iB2jwQKliHkXmBBdQYtWZr/99tvnvMlrmsbu3bsveM6pU6fIzc2N3M7Ly2Pbtm1Njnn++ecZOnQoo0aNak2ZI9LTvVh24o5oZmWlxLoI5zVz4qUoTeO1DYf4eH85melJfPHGQc2moj+f9qifpkFGhheXs/Ol9cjNTYt1EaImkesGiV+/1mpRoNi+fXurL2zbdpM3DqVUk9t79+5lzZo1PPfcc5w8ebLV1weorvbJxkUxMm5oHuWnfby/q5R/bDmGDnx+dK8Wn9+e9TMa0np0Jom8+U0i1w0Su35t3bio2Y9h3//+9yM/V1ZWturC+fn5lJWVRW6XlZWRl5cXuV1UVERZWRm33HILX/3qVzl16hS33357q55DxI6macy6rh+jLgvvnf72lmNs2F4Sk7L4AgaJ264UIvaaDRQ7duyI/Hzvvfe26sLjxo1j06ZNVFZW4vf7WbNmDRMnTow8vnDhQlavXs3y5ctZvHgxeXl5LF26tJXFF7GkaxrzPncpBX0ygHBeqA/3ln3GWe3PtBQhyf8kRNQ0GyjO3o+gtXsTdO/enfvvv5/58+dz0003MXPmTEaOHMmCBQva1JUlOieHrnPbjYPo3yPcp/u3fx5g+8GKDi9Hnd+QBXhCREmLO3ZbO1AJMGvWLGbNmtXkvt///vfnHNerVy/Wrl3b6uuLzsHl1LmzcDB/fG03x8rqeWntftxOncF9MjusDIZp4wuZpHSysQohEkGzLQrbtqmurqaqqgrLsiI/t2R6rOhaktxO7ppWQPeGJIJL3tzLwRMdm3G23mfEdF2HEIlKU830KRUUFKBp2nm7nT5remxH2HOwTGY9dTK1vhCLV+6iojqA26Vzz/Qh9Ol+7lTDaNUvye0gI9UNMd7aKJFnziRy3SCx69fWWU/NttOLi4vbXCDRNaV53dw7YwiLV+ykqi7Ec28Uc++MIfTMbf0fZ1sEQhYBwybJJdmNhWgvnW+Vkoh7Gake7p0xlG5eF4GQxR9fL6akouNaR3X1sc9uK0QikUAhoiI7PYl7Zg4lJdmFP2jyx9d2c+p0x+xlYdqKgEyXFaLdSKAQUZOXkcy9M87sZfGHVbsoq+qYYOELmCDL8IRoFzKXUERVfpaXe2cM4ZlVu6j1GzyzahcLZg69qDxPe46c5t2tJzhdGyQzzcOEUZecMxXXMC1CpsLtlDzk8WrbgXKK3j9CeXWAnPQkpl7bh5GX5sS6WF2StChE1F2Sk8K9M4aQ5HZQ6wsHi1OnfW261p4jp1mx4RA1foMkj5Mav8GKDYfYc+R0k+OUgkBIMsvGq20Hylny5l6q6kN4k5xU1YdY8uZeth0oj3XRuiQJFKJD9MxN5Z7pQ/C4HNT4DP536Ydt2n/73a0ncDh03E4HmqbhdjpwOHTe3XrinGMDQSuhswsnsqL3j+Bw6Hhc4d+zxxX+PRe9fyTWReuSJFCIDtMrL5V7ZhTgcTmoqg3yzMrWj1mcrg3icjT9s3U5wlu0fpqtlLQq4lR5dQD3p1LHu5065dWBGJWoa5NAITpU77w07plRQJIn3LJ4ZuUuTrUiWGSmeTCspossDcsmM81z3uP9AQkU8SgnPemc3QtDpk1OelKMStS1SaAQHa53XhrfufXK8JiFPxwsSitbNmYxYdQlWJZNyLRQShEyLSzLZsKoS857vGkr/NKqiDtTr+2DZdkEjfDvOWiEf89Tr+0T66J1SRIoREz069GtYeqsgzq/we9X7eJE+WcvyhvcJ5PZ1/WnW7KLQNCkW7KL2df1bzYBYb3PQMlU2bgy8tIc7pg8iIwUN76ASUaKmzsmD5JZTzHSbK6nzk5yPcWvxvqdKK/nj6/vxhcwSXI7uHt6Ab3z2n8byjSvi5QkV7tf90ISOV9QItcNErt+UdnhTohouyQnhQWzhpKW3JDu47ViDpXUtPvz1AdMbMksK0SbSKAQMdc908uC2UNJT3ETNCyee734nHURF8u2FXUBkzZsqyJElyeBQnQKOenJfHX2ULK7JWFYNi+s3su2A+27U54/aJ4zk0YI8dkkUIhOIzMtia/OHkp+lhdbKZa9vY/Nxafa7fpKQa1PtkwVorUkUIhOJc3r5iszh9I7LxUFvPLOQdZ9dLzVe7bH/DHWAAAgAElEQVRfSMiwqAsY0gUlRCtIoBCdjjfJyT0zhnBZz3QA1mw+ymubPsFup2DhCxj4Q5KGXIiWkkAhOiWPy8H8qYMZMSAbgI07TvLS2v2Y1sWPMSgFNfUhTMkDJUSLSKAQnZbTofOlz1/GmGHdAdh2oII/FRW3S/4m21ZU1wVlIZ4QLSCBQnRquqYxa1w/Jl/VG4ADx2tYvGIX1e2w3alh2g0bHAkhmiOBQnR6mqZx/ZU9mfe5S9E1jZOVPn736g5OtjA/VHPq/SaWLMQTolkSKETcuHJQLl+eNhiPy0F1fYinl+9k79Gqi7qmrRS1fpkFJURzJFCIuDKwVwZfnT2Ubg2ruJ8vKub9XaUXdc1gyCIos6CEuCAJFCLu9MhO4es3DeeSbC+2guXrD/HapsPYbZzFpBTU+C9+zEOIRCWBQsSl9BQ3C2YPY0jfcHrxDdtP8vzqPW2eEWWa0gUlxIVIoBBxy+NycMfkQYwf2QOAvUereOrVHW3aixvA5zeol1lQQpxDAoWIa7quMX1MX26ZNACHrlFWFeC3r+xg37HWD3IroNYXImDIeIUQZ5NAIRLC6MF5LJg1lNSGfS2ee724TTmilIKaulDCboglRFtIoBAJo0/3NL5x83B65aagCOeIWvLm3lbPaLKVoqo+KCnJhWgggUIklIxUDwtmDeOqwbkA7Dp8midf2U5pKxfn2baiqi4oLQshkEAhEpDLqXPzxAHcNKE/Dl2jvDrAb1/dwUd7y1p1ncZgEZQxC9HFRTVQrFy5kunTpzNlyhSWLFlyzuNvvfUWc+bMYfbs2XzjG9+guro6msURXYimaVwzpDtfnT2MjFQ3hmnz13UHeOWdgxit6FKyVThY+EOyjarouqIWKEpLS3n88cdZunQpr776KsuWLWP//v2Rx+vq6vjRj37E4sWLWbFiBYMHD+Y3v/lNtIojuqjeeal8c+4IBvfOAGBz8SmeenUHpadb3hXVOMDdHokIhYhHUQsUGzduZMyYMWRkZOD1eiksLKSoqCjyuGEY/PCHP6R793AK6cGDB1NSUhKt4oguzJvk4s6pgym8pje6Bicrffz27zvYXHyqxbOiFOALmFTWBiWJoOhyohYoTp06RW5ubuR2Xl4epaVncvJkZmYyefJkAAKBAIsXL+bGG2+MVnFEF6drGpMu78mCWQ1dUZbNK+8c5M9v7cMXMFp8nZBhUVET7opC9rIQXYQzWhe2bRvtrE5dpVST241qa2v5t3/7NwoKCrj55ptb9Rzp6V6sBN6lLCsrJdZFiKpY1C8rK4XBA7J58Y1iPtxzih2HKjlaVs+XZwxhaP/sVl1LczlIT/Xgcp7/81Zublp7FLlTSuS6QeLXr7WiFijy8/P54IMPIrfLysrIy8trcsypU6e49957GTNmDA8++GCrn6O62pew0xezslKorKyPdTGiJtb1u2Vif/p1T2XVpsNU1wV5YtnHjB2eT+E1vXE7HS2+TlmZRlqKmyS3g7N7sXJz0ygrq23/gncCiVw3SOz66bpGdnZq68+LQlkAGDduHJs2baKyshK/38+aNWuYOHFi5HHLsrjvvvuYNm0aixYtOm9rQ4ho0TSNqwry+NYtI+nTPfzC2bTjJL/523Y+OdnyNwmzYQptdX0Iu5WrwIWIF1FrUXTv3p3777+f+fPnYxgG8+bNY+TIkSxYsICFCxdy8uRJdu3ahWVZrF69GoDhw4fz05/+NFpFEuIc2d2SWDBrGO98fIK1Hx6jojrA4hU7GT+yBzde1fuC3UpnUyo80B0yLLp53XjcLW+RCBEPNNXaZDidyJ6DZdL1FKc6Y/027Sxh9ftHI6k70rwuxg3LZ9+xKk7XBslM8zBh1CUM7pN5wWtoWjirbZ+eGdS0IYvttgPlFL1/hPLqADnpSUy9tg8jL81pc53a04r1B1mz+RgBwyLJ5WDK1b2YPX5ArIvV7qTr6TznRaEsQsSdPUdOs35bCWkpLlKTXQDU+gxWbz7K8QofbpeDGr/Big2H2HPk9AWvoxQEQhYV1QFq/UaruqO2HShnyZt7qaoP4U1yUlUfYsmbe9l2oPyi63exVqw/yIqNhwkaFk4dgobFio2HWbH+YKyLJjqABAohgHe3nsDh0PG4nHRLcZObkUzjqFkwZFFWFcA0Fbqu8e7WE595PQXU+w0qqgP4QmaLJtIWvX+koQwONE3D43LgcOgUvX/kYqrWLtZsPoaGhkPX0DQ9/B2NNZuPxbpoogNEbYxCiHhyujZIkufMy6FxbKIxWDSm8nA7dYxWdHdatqKmLoTfaZKa7MLjuvD4RXl1AG9S05ek26lTXh1oeUWiJBAycehNJ5zoGm3eUVDEF2lRCAFkpnkwrKYBwOnQcDo08jKTSWoYoA6ZNjU+g6L3P2lV+nLDtKmqDVJZe+EkgznpSeekNg+ZNjnpSa2sTftLcjv59JIlW4XvF4lPAoUQwIRRl2BZNiHTQilFyLRwuxy43U4spchM89AtxUXjh+p3tpbwv8s+5sO9ZS0eh1CEV3ZX1QapqA0QMCzUWZ1SU6/tg2XZBI1wGYKGhWXZTL22TxRq3DpTru6FQmHZCqXs8HcUU67uFeuiiQ7g+NGPfvSjWBeirSpO+xJ2ZXZyshu/v+WpJeJNZ6tfTnoyOelJlFb4qPMZZKS4mTamL0P7ZUbuy+6WxLQxfeme5eVoaR3+kMWuw6fZc7SK3IxkMtM8ket9Vv1sWxEIWQSCForwbJRLclLIy0zm2Kk6qutDZKV5mDtxQKeY9TS4TyYoxScn6whZiiSXg+lj+iTkrKeUFA8+X2ImgNQ0Da/X3frzZHps59QZp4+2p3ivX1VdkDfeO8L2gxWR+4b2y2TqNX3IyUhudf00DTxOB94kJy6Xg868/DSRp49CYtevrdNjpYNRiDbISPVw240DGXcyn9ff+4Sjp+rYdfg0xZ+cZvTgPObeMLBV11MKAoZFwLBwOjVSklwkuZyyB4boFCRQCHER+uancd+cYWw/WMnqfx3hdG2QzcWn+GhfOWOGdWfiqEsi6zJayjQV1XUh6h0G3mQXHqcDp0Mjftv+It5JoBDiImmaxshLsxnaL5PNxaf4x4fHqfMbrN9Wwvu7Shk7rDvjR7YhYFjhqbW6puF06iS5HbhdOk5d5qCIjiWBQoh24nTojB2Wz+hBuXx8sJLV7x3GH7R4Z2sJ7+0s5dqh3bluZA+6tXIw0VaKkGERMiw0DRwODY/LicfpwOXS0DVpbYjokkAhRDtzuxxMHduPkf0z2bjjJOu3lRAIWby7rYRNO08yenAeE0b2IKtb69dHKBXumjJNg3oMdF3D7dTxuB24nNLaENEhgUKIKElyO7nhyl6MG57Pph2lbNhRgi9g8v6uUv61u5Th/bOZMLIHvfJaPwulUWSabaihtaFpDes/HLh0DYdDAoe4eBIohIiyJLeT66/syXUj8tlcfIp3t5VQUx9i+8EKth+soF9+GuNG9GBI38xz0mS0hlJgKoUZNPEFzTPdVE4nLpeOy6E1XF+mUonWkUAhRAdxuxxcN6IH1w7tzvYDFazfXkJJhY/DJ2s5fLKWjFQ31w7tztUFeXiTWjfwfT5nd1MRIBI43E4HbpcDp67hdOiRKbgyziEuRAKFEB3M6dC5YlAulw/M4cDxGjbsKGHvkSqq6kKs/tdR3t5yjOH9s7lmaB59u6e12+6PZwKHiS8QbnFomoaugaZrODQNl8uBy6HjbGh9SPAQIIFCiJjRNI3LeqVzWa90KqoDvLfzJB/sKSNoWHy8v5yP95eTl5nMVYPzuHxgTqun134WpUAphQ1gKQzCe2loNAQOPTze4XI2phUPCzbkqJLZVl1HXKfw2HuonNCnM3GqJt/iVrynuPgsUr/zCxkW2w5U8P7uUo6XnTlf1zQK+mZw5aBcBvXOwNnBg9RnN2qyMlM4XeXD6dBxu/RIIHHq4W6s+H1HCZMUHueK6xZFRqo7nBRQnQkMTf5IVdMf1Hn6Yhvj5NmHnvn5zGORc1TT886+nvpUlLI56/4LPc/5Alucv9BE27ldDq4qyOOqgjyOl9fzQfEptu4vJ9CQgHDX4dMke5yMGJDF5QNz6NM9Db0D8nw0ec0Qnm0Vsq3IB7XGbiyHQwt3XTVM1XXooGs6uh7/AaQri+sWRUVFHXYnzB7b2tetOjtoaArbhowMLxWV9dhKhb9shW2Fv1vKxrbPnKNQcfcilBZFyxmmza7DlXy4t4z9x6ub/K7TU9yMGJDNiEuz6ZWb0m7jGc1pTd0iAUQHp+5Abxj70HUNHQ208AZI4WO0TtEikRbFueK6RdFZteUPXTvrH4cOSR4nbue53Qtnz1BRhIOKUgpbEQkqyg7vGxAOKuHbZ4JKcwU/+8cLHPlZ1xDtzuXUGXVZDqMuy6HGF2Lb/go+3l/OifJ6qutDrN9ewvrtJaSnuBnaP4th/bLol5+GfhFTbdtLZBzEBoNzd8Nr/HvWCAcJ3aHh0vVwwNPODjThFei6RqtnasmsrosnLYpOqr0+1Zz9AVOppl1faCpy/9nHnDkgfFsj3G2nGnvQFCgt3IqxlYq8GahI0GoMXHb4vobLqbOuKS2Ki1de5WfbwQq2Hajg1Gl/k8e8HieD+2RQ0DeTgb3S23Unulj97ho/TGmahq6DU9fRGwbdtYZAAoBSke4xy1KYDe8RDk1Dd2iRzaeadCk3PocWHoOprvahN+wLrmvhT+KNM8Q4z3mN5zY6+/HG+22bSAvqfOyGF4pSDQGyIXhe6PnO95zna1CefW5bWxQSKDqpeG7+nu8F09h9pgBlQ3pmMqdP+8IBBiItFdUYeOwzAchS6sz9kQs3bfV0tr/ijn4zPXXaz67Dlew8VMnx8qbPq2saffPTGNw7g0F9MuiemXxRXVRdKchHgtNZ3WS6ruM46/9P08J/3yiwUGeCS8Mxtq0wbTuylWx4VlnDBzDC/zT52+bMcXrDtGU+FTjgTGvNQkVeL40najS0xLRwN184kILL5aBnfnqr/08kUHRS8RwoWqKl9Wv6fqY1bDt6ViuGs1s64amekaCizgSms1s66nyBp0Hk5XCRXWyxfDM9XRuk+Eh4b4yDJ2rO2QUyzevisp7pXNrwlZ7SuiSFXSlQJBqXU6fg0txWnydjFKJTa/oxRjXM5W/4ZNWGT8WfPsVu6EpDaxpcbJsmXWyN7fqGD3dYDeNAlqkwTCvSvdEZZKZ5GDssn7HD8gmGLA6cqGbPkSr2Hq2iuj5Erc/go33lfLSvHIDs9CQG9OjGgEu60S8/jfRUz2c8g+hqJFCILuWcPumz/mkSfFqxTEERTsznCxgYZucJGAAet4Oh/bIY2i8LpRRlVQH2H69m/7FqDpZUEzJsKqoDVFQH2Fx8CghPO++X340++an0yUuje5b3onJQifgngUKIi6Shkex2kux2EjJsAobZKbcw1TSNvMxk8jKTGTc8H8u2OV5Wz6GSGg6eqOGTk7WETJuqulBkZTiA26nTMzeFXrmp9MpLZZiuoynVIVNxRecgYxSdlIxRxLesrBRKSmsIBC2CphUXf6eWrSipqOeTk7UcLqnlyKlaan3GeY/1epxckpPCJTkp9Mj20iM7hZz0pE4xJfdiyRjFuaRFIUQUOBw6HpcDj8sR3qHOtAkETEKW3WmDhkPXwq2G3FSuG9EDpRTV9SGOlNZy7FQ9R8vqOFFWj2HZ+IJmuAvreHXkfJdDJy8rmfxML92zvHTPSiYvI5luKW5pfcQ5CRRCRJmuaSS5HCS7HZiWwrDCQcOw7Mg6lM5I0zQyUj1kpHoYeWkOEG51BG1F8cEKjpfXc6KsnpOVPoKGhWGFu7LOzlEF4HE5yM1IIjcjmZz0ZHIafs7q5sHtdMSiaqKVJFAI0UGUCn9qd+jhoGErsG0bywJL2eHV9JbCsG2shnQtnY1D1+iVk4rXqXPloHAXhq0UVbVBSip8nKz0UXraR2mlj4rqALYKZ5s9VlbPsbJzu3O6pbjJ7uYhq1sSWWlJZHXzkJkW/kpNdklLpJOQQCFEDDSueHc0JM47e5qVpoU/uZumwrCsSNAIB5Mz60Q6C13Twm/03ZIY1j8rcr9p2ZRXByir8nPqtJ+yKj8V1QHKqv2EjHDKzJr6EDX1IQ6VnDte5XScadGkp7rD31PcdEtxR74nuR0STDqABAohOhmlwm++bpeG2xUOIGengbBsG8OyMQw7HEg6YfCAcE6m/Cwv+VneJvcrpaj1G5FpuRU1ASprgpyuDVBRE8QfDOeEMi1FeXWA8urABZ/D5dTp5nWTluIiLdlNmtfV8OUmJclJqtdNarKLlCRnh6dmTyQSKISIA41BQGtIiud06CS7iWRbNS0bSylM08YwbUzTbkh9Ettyn4+maXTzuunmddO/R7dzHg+ETKrqQlTVBqmsDVJdF6SqLkR1fZDquhC1vlAkHYZh2lTUhIPNZ0lyO0hJcpGS7MTrceJNcuFNavw5/D3Z46TesAkFDJI94S1jOyKNe2cngUKIONYYCJwOHSfgcToiwcOyw+MexlnBwz5r9XlnleR2kp/lPKcl0si2FXUBg5qGoFHtC1Fbb1DrN6jzhVee1/nDX2enLwmELAIhi4qalpdF08IBJsntbPh+5meP20GSK/zd3TDDrfHL7dIbvod/djsdOB1a3HaTRTVQrFy5kqeeegrTNPnyl7/MHXfc0eTx3bt3s2jRIurr67nqqqv48Y9/jNMpsUuIi7F1fzlF7x+hvDpATnoS08b0YeSlOTz/xm62HTyNs2FHuiF90xncN4v1W09QUR0gzetm7PB8BvbKYM+R07y79QSna4NkpnmYMOoSjpfVsX7bSYKmhcfpYPzIfG4Y3fuC5TjfNYBz7hvcJ7PF5w/uk8m6j46Fy2FYeFzhcky++txyKKXYfrCC9dtOUFUXIsXjpP8l6aR53Rwrq+PoqVoCISuy5atp2QSCVpMgqhT4gxb+oHXO9VtL08DtdOB26rgagofLqZ/5coS/Ox1n7nOefZ9Dw+nQcTT87GhoWToj9zd81898dzi0hnGwiwtQUVtwV1paym233cbf//533G43t956K//7v//LZZddFjlm5syZPPzww1x++eU8+OCDDB8+nNtvv73FzyEL7uKX1C86th0oZ8mbe3E4dNxOnZBpY1k22Wkeio9WNznW6dBIcunkZqWQ5NKxFGgorhyYy9YDFdCQVjtoWNTUBfGH7EieeaNhgH3S5ZfwuSt6nVOOPUdOs2LDoYY3NR3DsvEHTNA0kj2OyH2WZTP7uv7nBIvznW9ZNn3zUtl6sDKcWVUL5+pSwOev7HlO0LrQNUYPymXL3rJz7p99XX8G9s4gOcXDiZM1BIIm+45VsWH7yci+GGbDOpheeal4XA4CofAuf0HDImjYBBtud7Z3JQ1wODTys1N46oHPt/r8qH1837hxI2PGjCEjIwOAwsJCioqK+OY3vwnA8ePHCQQCXH755QDMnTuXJ554olWBIhFWgTZH6hffYlG/93eVkpflbbI+IWRalJ32k5eZzNklanwzyzgrCWDItNi4o5TUFBdJLie6ruFNdhEybJKTNDxOvWFvhvAn9lOnA+RmJDfJ2otSHCurZ1CfTFwOPZJevrJhUDq9myeyr0nItNh3rCqyTqOxU2zP0SouyU3F5dAj3WuGZVFaFSA/24tDO1N+y1LsOVpN4bV9m/xf7DxUed7/iz1Hq897/85DlQwfkE03rxu7odvr/V2l9Oqees6xaUkuvvT5gef9HSjCY0Yhw8IwbUIN3X4h08Iww+toTNPCsM6MKRmWjdl427IwzfA+GuGxpoYxJ0th2ArLshtmxdkN2ZRbLs3ratXxjaIWKE6dOkVu7pml4nl5eWzbtu2Cj+fm5lJaWtqq58jMTLn4gnZibdlgJJ5I/drfg/eM6fDnPJ/v3D76os5fNKD1aSY+7YdtSFXRKD875aKvkUiiNl/Mtu0mAzfqU0nEPutxIYQQnUPUAkV+fj5lZWWR22VlZeTl5V3w8fLy8iaPCyGE6ByiFijGjRvHpk2bqKysxO/3s2bNGiZOnBh5vGfPnng8HrZs2QLA8uXLmzwuhBCic4hqmvGVK1fy9NNPYxgG8+bNY8GCBSxYsICFCxcyYsQIiouLeeihh6irq2PYsGE88sgjuN2t25ZRCCFEdMX1fhRCCCGiT5KfCCGEaJYECiGEEM2SQCGEEKJZEiiEEEI0K24Cxa9//WumT5/OjBkzePbZZ4FwmpBZs2YxZcoUHn/88RiX8OI9+uij/Nd//RcQTpg4d+5cCgsLWbRoEaZpxrh0bXfnnXcyY8YM5syZw5w5c9i6dSsrV65k+vTpTJkyhSVLlsS6iBdl7dq1zJ07l2nTpvHwww8DifO3+de//jXye5szZw6jR4/mv//7vxOmfhCemj9jxgxmzJjBo48+CiTO62/x4sUUFhYya9YsnnrqKaCNdVNx4P3331e33nqrMgxD+f1+df3116vdu3erSZMmqSNHjijDMNQ999yj1q1bF+uittnGjRvVtddeqx544AGllFIzZsxQH330kVJKqe9973tqyZIlsSxem9m2rcaPH68Mw4jcd/LkSXX99der06dPq/r6ejVr1iy1b9++GJay7Y4cOaLGjx+vSkpKVCgUUrfddptat25dQv1tNtq7d6+aPHmyOnHiRMLUz+fzqauvvlpVVFQowzDUvHnz1IYNGxLi9bdhwwY1c+ZMVVtbq0zTVF/72tfU6tWr21S3uGhRXHPNNTz//PM4nU4qKiqwLIuamhr69u1L7969cTqdzJo1i6KiolgXtU2qqqp4/PHHue+++4DzJ0yM17odPHgQgHvuuYfZs2fz4osvNkkY6fV6Iwkj49Gbb77J9OnTyc/Px+Vy8fjjj5OcnJwwf5tn+9GPfsT999/P0aNHE6Z+lmVh2zZ+vx/TNDFNE6fTmRCvv127djF+/HhSU1NxOBxMmDCBF154oU11i4tAAeByuXjiiSeYMWMGY8eOPW/SwdYmFewsfvCDH3D//ffTrVt4t6/2SJjYWdTU1DB27FiefPJJnnvuOf7yl79w4sSJhPndffLJJ1iWxX333cecOXNYunRpQv1tNtq4cSOBQIBp06YlVP1SU1P59re/zbRp05g0aRI9e/bE5XIlxOtv2LBhrF+/nqqqKoLBIGvXrsXpdLapbnETKAAWLlzIpk2bKCkp4fDhwwmRVPCvf/0rPXr0YOzYsZH7Eilh4hVXXMEvfvEL0tLSyMrKYt68eTzxxBMJUz/Lsti0aRM/+9nPWLZsGdu2bePo0aMJU79Gf/nLX7j77ruBxPr7LC4u5m9/+xv/+Mc/ePfdd9F1nQ0bNiRE/caOHcvcuXO58847+cpXvsLo0aMxTbNNdYuL7eQOHDhAKBRiyJAhJCcnM2XKFIqKinA4zuSI/3TSwXjx+uuvU1ZWxpw5c6iursbn86FpWsIkTPzggw8wDCMSCJVS9OzZs9mEkfEkJyeHsWPHkpWVBcCNN96YMH+bjUKhEJs3b+bnP/858NkJP+PJ+vXrGTt2LNnZ2UC4K+YPf/hDQrz+6urqmDJlSiTAP/PMM/Tq1YsPPvggckxL6xYXLYpjx47x0EMPEQqFCIVCvP3229x6660cOnQo0vRftWpVXCYVfPbZZ1m1ahXLly9n4cKF3HDDDTzyyCMJkzCxtraWX/ziFwSDQerq6njllVf45S9/2WzCyHhy/fXXs379empqarAsi3fffZepU6cmxN9moz179tCvXz+83vBmPqNGjUqY+hUUFLBx40Z8Ph9KKdauXcs111yTEK+/Y8eO8Y1vfAPTNKmtreXll19m3rx5bapbXLQoJk2axLZt27jppptwOBxMmTKFGTNmkJWVxbe+9S2CwSCTJk1i6tSpsS5qu3nssceaJEycP39+rIvUJtdffz1bt27lpptuwrZtbr/9dkaPHs3999/P/PnzIwkjR44cGeuitsmoUaP4yle+wu23345hGFx33XXcdtttDBgwIGH+No8ePUp+fn7ktsfj4ec//3lC1G/8+PHs2rWLuXPn4nK5GDFiBF/96leZPHly3L/+CgoKmDJlCrNnz8ayLO666y5Gjx7dpvcWSQoohBCiWXHR9SSEECJ2JFAIIYRolgQKIYQQzZJAIYQQolkSKIQQQjQrLqbHCtESDz/8MJs3bwbCizR79uxJUlISAMuWLYv83Bkppbj77rt54oknIqlchOgsZHqsSEg33HADv/71rxkxYkSsi9IipmkybNgwNm/eLIFCdDrSohBdwr59+/jpT38aWUF91113cfPNN7Nx40Z+85vfkJubyyeffILX6+UrX/kKL7zwAocPH2batGk88MADbNy4kSeeeIK8vDwOHTqE1+vlkUceYcCAAYRCIX7xi1+wZcsWLMti2LBhLFq0iNTUVCZOnMjo0aMpLi7mu9/9LrZt88wzzxAKhaisrOSWW27hW9/6Ft/73vcAuOOOO3jmmWf4whe+wNNPP82QIUMAmDhxIk8//TRer5e7776bPn36UFJSwtKlSzl06BD/8z//QyAQQNd1Fi5cyKRJk2L53y0STftnQRci9q6//nq1bds2pZRSoVBITZs2Te3evVsppVR1dbUqLCxU27ZtUxs2bFBDhw6NPHbXXXep2267TYVCIVVeXq6GDBmiysvL1YYNG1RBQYHasmWLUkqpF154QX3hC19QSin1q1/9Sv3yl79Utm0rpZR69NFH1U9+8hOllFITJkxQv/vd75RSSlmWpW6//XZ15MgRpZRSJ06cUAUFBaqqqkoZhqEGDRqkqqurI+ft2rUrUp/G24cPH1aDBg1SH374oVJKqcrKSjVlyhR1/PhxpZRSJSUlasKECaqkpCRK/7OiK5IWhUh4Bw4c4OjRozzwwAOR+0KhELt376ZXr1706dOHgoICAHr37k1OTg4ul4vs7Gy8Xi9VVVVAOG3zlVdeCd6xJJkAAAJeSURBVMAXvvAFHn74YWpra1m3bh0+n493330XAMMwmiRaGz16NAC6rvP000+zbt06li9fzv79+1FKEQgESElJaXF9XC4Xo0aNAuDDDz+krKyMr3/965HHdV1n7969TdJuCHExJFCIhGfbNhkZGSxfvjxyX1lZGd26dWPLli243e4mxzud539ZnH2/bdtA+E3Zsix+8IMfcN111wHhrJ2GYUSObQwCdXV13HzzzRQWFjJ69GhuueUW3nzzTdR5hgk1TWty/9nXS0pKQtf1SDkGDRrEX/7yl8jjpaWlkWy2QrQHmR4rEt5ll12Gruu89tprQHgHwZkzZ1JcXNyq6+zYsYN9+/YB4VlUV199NSkpKYwfP54XXngBwzCwLIsHH3yQX/3qV+ecf+jQIfx+P9/+9re5/vrr2bRpE6ZpYlkWDocDTdMi+xdnZWWxY8cOILxpUGVl5XnLdMUVV3DgwIFINtCdO3dSWFhIRUVFq+omRHOkRSESntvt5qmnnuJnP/sZv/vd7zBNk3//939n1KhRbNy4scXXycvL47HHHuP48ePk5uby6KOPAvCtb32LRx99lJtuuikymP2f//mf55w/dOhQxo8fz7Rp03C5XBQUFDBgwACOHDlCz549mTJlCrfddhu//e1v+e53v8uPf/xjlixZwogRIyKD2p+Wk5PDE088wSOPPEIoFEIpxWOPPSbdTqJdyfRYIVpg48aNPProo026r4ToKqTrSQghRLOkRSGEEKJZ0qIQQgjRLAkUQgghmiWBQgghRLMkUAghhGiWBAohhBDNkkAhhBCiWf8/to/SyKjM1ykAAAAASUVORK5CYII=\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/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/5c9dbef11b4d7638b7ddf2ea71026e7bf00fcfb0/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.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment