diff --git a/module4/src_Python3_challenger.ipynb b/module4/src_Python3_challenger.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..001f33a28d4a6f477e0e4b2cc724a5e1872b1b9b --- /dev/null +++ b/module4/src_Python3_challenger.ipynb @@ -0,0 +1,895 @@ +{ + "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.7.1 (default, Dec 10 2018, 22:54:23) [MSC v.1915 64 bit (AMD64)]\n", + "uname_result(system='Windows', node='ISPCM361631', release='10', version='10.0.17134', machine='AMD64', processor='Intel64 Family 6 Model 158 Stepping 10, GenuineIntel')\n", + "IPython 7.9.0\n", + "IPython.core.release 7.9.0\n", + "PIL 6.2.1\n", + "PIL.Image 6.2.1\n", + "PIL._version 6.2.1\n", + "_csv 1.0\n", + "_ctypes 1.1.0\n", + "decimal 1.70\n", + "argparse 1.1\n", + "backcall 0.1.0\n", + "bottleneck 1.3.1\n", + "cffi 1.13.2\n", + "colorama 0.4.1\n", + "csv 1.0\n", + "ctypes 1.1.0\n", + "cycler 0.10.0\n", + "dateutil 2.8.1\n", + "decimal 1.70\n", + "decorator 4.4.1\n", + "distutils 3.7.1\n", + "ipykernel 5.1.3\n", + "ipykernel._version 5.1.3\n", + "ipython_genutils 0.2.0\n", + "ipython_genutils._version 0.2.0\n", + "ipywidgets 7.5.1\n", + "ipywidgets._version 7.5.1\n", + "jedi 0.15.1\n", + "joblib 0.14.0\n", + "joblib.externals.cloudpickle 1.2.2\n", + "joblib.externals.loky 2.6.0\n", + "json 2.0.9\n", + "jupyter_client 5.3.4\n", + "jupyter_client._version 5.3.4\n", + "jupyter_core 4.6.1\n", + "jupyter_core.version 4.6.1\n", + "kiwisolver 1.1.0\n", + "logging 0.5.1.2\n", + "matplotlib 3.0.2\n", + "matplotlib.backends.backend_agg 3.0.2\n", + "mkl_fft 1.0.15\n", + "mkl_fft._version 1.0.15\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.25.3\n", + "_libjson 1.33\n", + "parso 0.5.1\n", + "patsy 0.5.1\n", + "patsy.version 0.5.1\n", + "pickleshare 0.7.5\n", + "pkg_resources._vendor.packaging.__about__ 16.8\n", + "pkg_resources._vendor.six 1.10.0\n", + "pkg_resources._vendor.appdirs 1.4.3\n", + "pkg_resources._vendor.packaging 16.8\n", + "pkg_resources._vendor.pyparsing 2.2.1\n", + "pkg_resources._vendor.six 1.10.0\n", + "platform 1.0.8\n", + "prompt_toolkit 2.0.10\n", + "psutil 5.6.5\n", + "pygments 2.4.2\n", + "pyparsing 2.4.5\n", + "pytz 2019.3\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.13.0\n", + "statsmodels 0.10.1\n", + "statsmodels.__init__ 0.10.1\n", + "statsmodels.api 0.10.1\n", + "statsmodels.tools.web 0.10.1\n", + "traitlets 4.3.3\n", + "traitlets._version 4.3.3\n", + "urllib.request 3.7\n", + "zlib 1.0\n", + "zmq 18.1.0\n", + "zmq.sugar 18.1.0\n", + "zmq.sugar.version 18.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": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
56/18/82672500
68/30/836731000
711/28/836701000
82/03/846572001
94/06/846632001
108/30/846702001
1110/05/846782000
1211/08/846672000
131/24/856532002
144/12/856672000
154/29/856752000
166/17/856702000
177/2903/856812000
188/27/856762000
1910/03/856792000
2010/30/856752002
2111/26/856762000
221/12/866582001
\n", + "
" + ], + "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": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# le chargement par url n'a pas marché, j'ai du télécharger les données et les charger en local.\n", + "\n", + "data = pd.read_csv(\"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": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGBNJREFUeJzt3XuQnXWd5/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": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "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": 29, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\rollandm\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:7: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " import sys\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Tue, 26 May 2020 Deviance: 3.0144
Time: 09:28:42 Pearson chi2: 5.00
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
" + ], + "text/plain": [ + "\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: Tue, 26 May 2020 Deviance: 3.0144\n", + "Time: 09:28:42 Pearson chi2: 5.00\n", + "No. Iterations: 6 \n", + "Covariance Type: nonrobust \n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", + "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 29, + "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", + "# pour que l'affichage de la ligne fonctionne bien après, j'ai du inverser l'ordre des données (temperature Intercept)\n", + "# dans les régressions logistiques.\n", + "\n", + "logmodel=sm.GLM(data['Frequency'], data[['Temperature','Intercept']], \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": 31, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\rollandm\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " \n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Tue, 26 May 2020 Deviance: 18.086
Time: 09:29:10 Pearson chi2: 30.0
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
" + ], + "text/plain": [ + "\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: Tue, 26 May 2020 Deviance: 18.086\n", + "Time: 09:29:10 Pearson chi2: 30.0\n", + "No. Iterations: 6 \n", + "Covariance Type: nonrobust \n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", + "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "logmodel=sm.GLM(data['Frequency'], data[['Temperature','Intercept']], \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": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Temperature Intercept Frequency\n", + "0 30.0 1 0.834373\n", + "1 30.5 1 0.826230\n", + "2 31.0 1 0.817774\n", + "3 31.5 1 0.809002\n", + "4 32.0 1 0.799911\n", + ".. ... ... ...\n", + "116 88.0 1 0.006133\n", + "117 88.5 1 0.005791\n", + "118 89.0 1 0.005467\n", + "119 89.5 1 0.005162\n", + "120 90.0 1 0.004873\n", + "\n", + "[121 rows x 3 columns]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VOXd//H3N5N9JyQsYUchyCZrQHHBuoA+FZeqiAu4ULV1qU9bW3laq2311wXb2rqjoIgrtYjYWsGF1KosAUFWw6IsCULYk0D23L8/ZsAQAlmYMJnJ53VduTLnzJkz3zsn+czJfc65jznnEBGR0BUW6AJERKRpKehFREKcgl5EJMQp6EVEQpyCXkQkxCnoRURCXJ1Bb2bTzCzfzFYd43kzs7+Z2QYzW2Fmg/xfpoiINFZ99uhfBEYf5/mLgR6+r9uAp0+8LBER8Zc6g9459zGw5ziLXAa85LwWAslm1t5fBYqIyIkJ98M6OgBbq03n+uZ9U3NBM7sN714/0dHRgzt37uyHt2+eqqqqCAsL3UMgody+UG4bqH3Bbt26dbucc2kNeY0/gt5qmVfruArOuSnAFICMjAyXk5Pjh7dvnrKyshg5cmSgy2gyody+UG4bqH3Bzsw2N/Q1/vjYywU6VZvuCGzzw3pFRMQP/BH0c4DxvrNvhgP7nXNHdduIiEhg1Nl1Y2avASOBVDPLBR4EIgCcc88A7wKXABuAg8DNTVWsiIg0XJ1B75wbV8fzDrjTbxWJSFAoLy8nNzeXkpKSQJdyhKSkJNauXRvoMk5YdHQ0HTt2JCIi4oTX5Y+DsSLSAuXm5pKQkEDXrl0xq+2cjMAoLCwkISEh0GWcEOccu3fvJjc3l27dup3w+kL3HCQRaVIlJSW0bt26WYV8qDAzWrdu7bf/lhT0ItJoCvmm48+frYJeRCTEqY9eRIKWx+OhX79+h6dnz55N69atA1hR86SgF5GgFRMTw/Lly4+YV1hYePhxRUUF4eGKOXXdiEhIeeWVV7j66qu59NJLueiiiwCYPHkyQ4cOpX///jz44IOHl33kkUfIyMjgggsuYNy4cTz66KMAjBw5kiVLlgCwa9cuunbtCkBlZSX33Xff4XU9++yzwLfDLlx11VX06tWL66+/Hu+Z55Cdnc2ZZ57J6aefTmZmJoWFhZx99tlHfECNGDGCFStWNNnPRB91InLCfv3OatZsK/DrOnunJ/LgpX2Ou0xxcTEDBgwAoFu3brz11lsALFiwgBUrVpCSksK8efNYv349ixcvxjnHmDFj+Pjjj4mLi+P1119n2bJlVFRUMGjQIAYPHnzc95s6dSpJSUlkZ2dTWlrKiBEjDn+YLFu2jNWrV5Oens6IESP49NNPyczMZOzYsbzxxhsMHTqUgoICYmJimDhxIi+++CKPPfYY69ato7S0lP79+/vhp1Y7Bb2IBK3aum4ALrzwQlJSUgCYN28e8+bNY+DAgQAUFRWxfv16CgsLueKKK4iNjQVgzJgxdb7fvHnzWLFiBW+++SYA+/fvZ/369URGRpKZmUnHjh0BGDBgAJs2bSIpKYn27dszdOhQABITEwG4+uqr+e1vf8vkyZOZNm0aN91004n9IOqgoBeRE1bXnvfJFhcXd/ixc45JkyZx++23H7HMY489dsxTGMPDw6mqqgI44lx25xyPP/44o0aNOmL5rKwsoqKiDk97PB4qKipwztX6HrGxsVx44YW8/fbbzJw583A3UVNRH72IhLRRo0Yxbdo0ioqKAMjLyyM/P59zzjmHt956i+LiYgoLC3nnnXcOv6Zr164sXboU4PDe+6F1Pf3005SXlwOwbt06Dhw4cMz37tWrF9u2bSM7OxvwHiiuqKgAYOLEidxzzz0MHTr08H8fTUV79CIS0i666CLWrl3LGWecAUB8fDwvv/wygwYNYuzYsQwYMIAuXbpw9tlnH37NT3/6U6655hpmzJjBd77zncPzJ06cyKZNmxg0aBDOOdLS0pg9e/Yx3zsyMpI33niDu+++m+LiYmJiYvjggw+Ij49n8ODBJCYmcvPNJ2EcSOdcQL569uzpQtn8+fMDXUKTCuX2hXLbnPNf+9asWeOX9fhbQUFBo1734IMPusmTJ/u5mmPLy8tzPXr0cJWVlcdcprafMbDENTBv1XUjInKSvfTSSwwbNoxHHnnkpNz2UF03IiLAQw89dNLea/z48YwfP/6kvZ/26EWk0Zyr9fbQ4gf+/Nkq6EWkUaKjo9m9e7fCvgk433j00dHRflmfum5EpFE6duxIbm4uO3fuDHQpRygpKfFbQAbSoTtM+YOCXkQaJSIiwi93P/K3rKysw1fBipe6bkREQpyCXkQkxCnoRURCnIJeRCTEKehFREKcgl5EJMQp6EVEQpyCXkQkxCnoRURCnIJeRCTEBSzoSyoC9c4iIi1LwIJ++8Eqfv7mCvYdLAtUCSIiLULAgj4p0njz81wu+PN/+OeKbRrqVESkiQQs6FtFG3PuGkF6cgx3vbqMO15eSn5hSaDKEREJWQE9GNsnPYlZPziT+y/uxfycnVz45495e3me9u5FRPwo4GfdhHvCuOPcU/j3j86me1ocP3p9OXe++jl7DqjvXkTEH+oV9GY22sxyzGyDmd1fy/OdzWy+mS0zsxVmdklDCzklLZ437ziTn43O4P01O7joLx8z/8v8hq5GRERqqDPozcwDPAlcDPQGxplZ7xqL/RKY6ZwbCFwLPNWYYjxhxg9Hnsrbd55F67hIbn4xmwffXkVJeWVjViciItRvjz4T2OCc+8o5Vwa8DlxWYxkHJPoeJwHbTqSo3umJvH3XCG4Z0Y3pCzZz6eOfkLO98ERWKSLSYlldBz7N7CpgtHNuom/6RmCYc+6uasu0B+YBrYA44ALn3NJa1nUbcBtAWlra4JkzZ9ZZ4KpdFUxZUUZxheO6XpGM7BSOmdW7gYFSVFREfHx8oMtoMqHcvlBuG6h9we68885b6pwb0pDX1Ofm4LWlas1Ph3HAi865P5nZGcAMM+vrnKs64kXOTQGmAGRkZLiRI0fW+eYjgbGjSvnxzOVMX7OLXZ7W/P57/UiIjqhH6YGTlZVFfdoXrEK5faHcNlD7WqL6dN3kAp2qTXfk6K6ZW4GZAM65BUA0kOqPAgHSEqKYfnMmPx/di/dWb+fSxz9hzbYCf61eRCSk1Sfos4EeZtbNzCLxHmydU2OZLcD5AGZ2Gt6g3+nXQsOMH4w8hddvG05xeSWXP/Upb2Rv8edbiIiEpDqD3jlXAdwFzAXW4j27ZrWZ/cbMxvgW+wnwfTP7AngNuMk10VVPQ7um8K97ziazawo//8dKJs1aobNyRESOoz599Djn3gXerTHvV9UerwFG+Le0Y0uNj2L6LZn8aV4OT2VtZPW2Ap65YTDpyTEnqwQRkaAR8CtjG8sTZvxsdC+evXEwX+08wKWPf8Lir/cEuiwRkWYnaIP+kFF92jH7zjNJjInguucWMmPh5kCXJCLSrAR90AOc2iaB2XeO4OweqTwwexW/nL2S8sqqul8oItIChETQAyTFRPD8hKHcfm53Xl64hZteWKybmoiIEEJBD95++0kXn8bkq/qz+Os9XPHUZ3y960CgyxIRCaiQCvpDrh7SiVe/P5x9B8u44qlPdZBWRFq0kAx68J5v/9YPR5ASG8kNzy/i7eV5gS5JRCQgQjboAbqmxjHrh2cysHMyP3p9OU9lbdDdq0SkxQnpoAdIjo3kpVszGXN6On98L4cH3l5FZZXCXkRajnpdGRvsosI9PDZ2AOnJMTzzn43sKCjl8XEDiY7wBLo0EZEmF/J79IeEhRn3X9yLX4/pwwdrd3Dj1EXsP1ge6LJERJpciwn6Qyac2ZXHxw1k+dZ9XPPsArbvLwl0SSIiTarFBT3Ad/un8+LNmeTuPchVz3zGJp1rLyIhrEUGPcCIU1N57bbhHCit4KpnFuhGJiISslps0AP075jM3+84gwiPMXbKApZu1oVVIhJ6WnTQg3dAtDd/cCap8VHcOHUxn27YFeiSRET8qsUHPUCH5BjeuH04nVrFcvOL2XywZkegSxIR8RsFvU+bhGjeuH04p7VL4I6Xl/Lvld8EuiQREb9Q0FeTHBvJyxOHcXqnZO56bZnGxxGRkKCgryEhOoKXbslkaNdW3PvGct5cmhvokkREToiCvhZxUeG8cFMmZ52ayn1vfsHM7K2BLklEpNEU9McQE+nhufFDOOvUVH72jxW8vnhLoEsSEWkUBf1xREd4w/7cnmncP2slrynsRSQIKejrEB3h4dkbBzMyI41Js1aqG0dEgo6Cvh6iIzw8c8NgzumZxs9nrdABWhEJKgr6eoqO8DDlxsGHD9C+tUxhLyLBQUHfAN6wH8Lwbq35ycwv+OeKbYEuSUSkTgr6BoqJ9DD1piEM7tKKH72+nLmrtwe6JBGR41LQN0JsZDjTbhpKvw5J3PXq58zPyQ90SSIix6Sgb6SE6Aim35JJz7YJ3DFjKQs27g50SSIitVLQn4CkmAhm3DqMzimx3Do9m6Wb9wa6JBGRoyjoT1BKXCSvTBxGm4Qobn5hse5UJSLNjoLeD9okRvPyxGHER4UzftoiNu4sCnRJIiKHKej9pGOrWGZMHIZzcMPzi9hVXBXokkREgHoGvZmNNrMcM9tgZvcfY5lrzGyNma02s1f9W2ZwOCUtnhm3DuNAaQWTs0vYWVga6JJEROoOejPzAE8CFwO9gXFm1rvGMj2AScAI51wf4N4mqDUo9E5P5IWbh7K3xDFh2mL2F5cHuiQRaeHqs0efCWxwzn3lnCsDXgcuq7HM94EnnXN7AZxzLfrE8sFdUrh7YBTr8wuZOD2b4rLKQJckIi2YOeeOv4DZVcBo59xE3/SNwDDn3F3VlpkNrANGAB7gIefce7Ws6zbgNoC0tLTBM2fO9Fc7mp2ioiLWFEbz9Bel9E/zcPfAKMLDLNBl+U1RURHx8fGBLqNJhHLbQO0Lduedd95S59yQhrwmvB7L1JZONT8dwoEewEigI/BfM+vrnNt3xIucmwJMAcjIyHAjR45sSK1BJSsri599dyQdum/mF2+t4l87W/Gnq08nLETCPisri1DdfqHcNlD7WqL6BH0u0KnadEeg5mheucBC51w58LWZ5eAN/my/VBnErh/Whb0Hynh03jqSYiJ48NLemIVG2ItIcKhPH3020MPMuplZJHAtMKfGMrOB8wDMLBXoCXzlz0KD2Z3nncotI7rx4mebeCprY6DLEZEWps49eudchZndBczF2/8+zTm32sx+Ayxxzs3xPXeRma0BKoH7nHMa/MXHzPjl/5zG3oNlTJ6bQ+u4SK7N7BzoskSkhahP1w3OuXeBd2vM+1W1xw74se9LahEWZvzxqv7sOVDG/721klZxkYzq0y7QZYlIC6ArY0+iCE8YT98wiP4dk7n7tWUs/npPoEsSkRZAQX+SxUaG88JNQ+nYKoaJ07P5crsGQRORpqWgD4BWcZG8dEsmMZEeJkxbTO7eg4EuSURCmII+QDq2imX6LZkcLKtkwrTF7D1QFuiSRCREKegDqFe7RJ4fP4Ste4u5VUMliEgTUdAH2LDurfnr2AEs27qPu19bRkWlhjcWEf9S0DcDF/drz6/H9OGDtTt44O3V1DX+kIhIQ9TrPHppeuPP6Mr2/SU8lbWRdonR/OiCHoEuSURChIK+GblvVAY7Ckr5ywfraJsYpatnRcQvFPTNiJnx++/1Y2dRKf/31krSEqI4/7S2gS5LRIKc+uibmQhPGE9fP4g+6Unc+ernLNuyN9AliUiQU9A3Q3FR4Uy7aShpCVHcOn0JX+86EOiSRCSIKeibqbSEKF66ZRgAE6Yt1o3GRaTRFPTNWLfUOKZOGEJ+YQm3Ts/mQGlFoEsSkSCkoG/mBnZuxZPXDWJV3n7uevVzXVAlIg2moA8C55/Wlocv78f8nJ384q1VuqBKRBpEp1cGieuGdWb7/mL+9tEG2iVF878X9gx0SSISJBT0QeR/L+zJN/tL+OuH62mfFK0LqkSkXhT0QcTM+H9X9iO/sJRfzF5Fm8QovtNLF1SJyPGpjz7IRHjCeOr6QfRun8idryxj+dZ9gS5JRJo5BX0QOnRBVWpCJLe8mM0mXVAlIsehoA9SaQlRTL85E+cc43VBlYgch4I+iHVPi2fqTUN1QZWIHJeCPsgN6tyKJ8Z5L6j64SufU64LqkSkBgV9CLigd1seuaIf/1m3k0mzVuqCKhE5gk6vDBHjMjuz3XeOfdvEKO4b1SvQJYlIM6GgDyH3XtCD/MISnpy/kbaJ0Yw/o2ugSxKRZkBBH0LMjN9e1pedhWU8OGc1qfFRXNKvfaDLEpEAUx99iAn3hPH4uIEM7JTMvW8sZ+FXuwNdkogEmII+BMVEepg6YSidU2L5/ktL+HJ7QaBLEpEAUtCHqFZxkUy/JZO4yHAmTFtM7t6DgS5JRAJEQR/COiTHMP2WTIrLKhk/bTF7DpQFuiQRCQAFfYjLaJfA1JuGkre3mJtfWKyrZ0VaIAV9CzC0awpPXDeIVdsKuOPlpZRV6OpZkZZEQd9CXNi7Lb+7sh//Xb+Ln/z9C6qqdPWsSEtRr6A3s9FmlmNmG8zs/uMsd5WZOTMb4r8SxV+uGdKJ+y/uxTtfbOOhd1ZrqASRFqLOC6bMzAM8CVwI5ALZZjbHObemxnIJwD3AoqYoVPzjjnNPYc+BMqZ8/BUpcZHce4HuPSsS6uqzR58JbHDOfeWcKwNeBy6rZbnfAn8ESvxYnzSBSRf34qrBHXnsg/VM/2xToMsRkSZWnyEQOgBbq03nAsOqL2BmA4FOzrl/mtlPj7UiM7sNuA0gLS2NrKysBhccLIqKipp1+y5u7djYxsODc1aTt2kDZ6Y3bDSM5t6+ExHKbQO1ryWqz1+31TLvcOeumYUBfwFuqmtFzrkpwBSAjIwMN3LkyHoVGYyysrJo7u0bcXYlN7+QzdRVe8gc0I8Letf/RuPB0L7GCuW2gdrXEtWn6yYX6FRtuiOwrdp0AtAXyDKzTcBwYI4OyDZ/0REenpswhL7pifzw1c9ZsFHj4oiEovoEfTbQw8y6mVkkcC0w59CTzrn9zrlU51xX51xXYCEwxjm3pEkqFr+KjwrnhZsz6ZISy8Tp2Szfui/QJYmIn9UZ9M65CuAuYC6wFpjpnFttZr8xszFNXaA0vZS4SF6eOIzW8VFMmLaYnO2FgS5JRPyoXufRO+fedc71dM6d4px7xDfvV865ObUsO1J788GnbWI0r0wcRnREGDdMXcTXuw4EuiQR8RNdGSuHdUqJ5ZWJw6isclz/3EKNeCkSIhT0coRT2yQw49ZMikoruP75Rewo0GURIsFOQS9H6ZOexIu3ZLKrsJTrn1/ErqLSQJckIidAQS+1GtS5FVNvGkru3oPc8Pwi9mose5GgpaCXYxrevTXPjx/KV7sOcMPURew/WB7okkSkERT0clxn9Ujl2RsHs35HEeOnLWJ/8YmF/exleYz4/Ud0u/9fjPj9R8xeluenSqW50bZuPhT0UqfzMtrw1PWDWPNNAROmLaagpHFhP3tZHpNmrSRvXzEOyNtXzKRZKxUAIUjbunlR0Eu9XNC7LU9eN4hVefuZMG0xxRUNH8t+8twcissrj5hXXF7J5Lk5/ipTmglt6+ZFQS/1dlGfdjxx3SBW5u7n0eySBu/Zb9tX3KD5Ery0rZsXBb00yOi+3rDfVFDF+KkN68ZJT45p0HwJXtrWzYuCXhpsdN923DkgitXb9nPj1MX1PkB736gMYiI8R8yLifBw36iMpihTAkjbunlR0EujDGobztPXD2bttgKuf35hvc6zv3xgB353ZT86JMdgQIfkGH53ZT8uH9ih6QuWk0rbunlp2G2FRKq5oHdbpowfzG0zljLuuYW8PHEYqfFRx33N5QM76I+9hdC2bj60Ry8nZGRGG164aSibdh9g7LML2L5fY+OINDcKejlhI05NZfrNmewoKOXqZz9j6x6NeinSnCjoxS+GdW/NKxOHUVBcwVXPfMaGfN28RKS5UNCL35zeKZk3bh9OZRVc8+xCVubuD3RJIoKCXvysV7tE/n7HGcREeBj33ELdcFykGVDQi991S43jHz84k/ZJ0Ux4YTHzVm8PdEkiLZqCXppEu6RoZt5+Bqe1T+SOl5fy2uItgS5JpMVS0EuTaRUXyasTh3F2jzQmzVrJ3z5cj3MNHwxNRE6Mgl6aVFxUOM9PGMKVAzvw5/fX8YvZq6iorAp0WSItiq6MlSYX4Qnj0atPp01iNM/8ZyM79pfw+HUDiY3Ur5/IyaA9ejkpwsKM+y/uxW8v68P8nHzGTVlIfqGuohU5GRT0clLdeEZXnr1xCOt2FHHFk5+xbocurBJpagp6Oeku7N2WmbefQVllFd976jP+u35noEsSCWkKegmIfh2TmH3nCDq0iuGmF7KZsWBToEsSCVkKegmYDskxvPmDMxnZM40H3l7NA7NXUa4zckT8TkEvARUfFc6U8UO4/ZzuzFi4mQnTFtfrJiYiUn8Kegk4T5gx6ZLTePTq01myeS+XPvEJa7YVBLoskZChoJdm46rBHZl5+xlUVDq+9/RnvPPFtkCXJBISFPTSrAzolMycu0fQJz2Ru19bxm//uUb99iInSEEvzU6bhGhe/f5wbjqzK1M/+Zrrn1tEfoEurhJpLAW9NEuR4WE8NKYPf712ACvz9nPJ3z7hsw27Al2WSFCqV9Cb2WgzyzGzDWZ2fy3P/9jM1pjZCjP70My6+L9UaYkuG9CBt+8aQXJsBDdMXcTfPlxPZZVGwBRpiDqD3sw8wJPAxUBvYJyZ9a6x2DJgiHOuP/Am8Ed/FyotV8+2Cbx95wguG+AdAfOG5xexQ105IvVWnz36TGCDc+4r51wZ8DpwWfUFnHPznXMHfZMLgY7+LVNauriocP58zen88ar+LN+6j4v/+l8++nJHoMsSCQpW140gzOwqYLRzbqJv+kZgmHPurmMs/wSw3Tn3cC3P3QbcBpCWljZ45syZJ1h+81VUVER8fHygy2gygWzftqIqnv6ilK2FVZzfOZyxGZFEesxv69e2C26h3r7zzjtvqXNuSENeU58BwWv7C6r108HMbgCGAOfW9rxzbgowBSAjI8ONHDmyflUGoaysLNS+pnPlqEr++F4O0z79mi0l0Tx27QD6pCf5Zd2BbltTU/tanvp03eQCnapNdwSOupLFzC4AfgGMcc6V+qc8kdpFR3j41aW9eemWTPYXl3P5k5/yxEfrdfcqkVrUJ+izgR5m1s3MIoFrgTnVFzCzgcCzeEM+3/9litTunJ5pzL33HEb1acej89bxvWcWsCG/KNBliTQrdQa9c64CuAuYC6wFZjrnVpvZb8xsjG+xyUA88HczW25mc46xOhG/axUXyRPXDeLxcQPZvPsAl/ztvzydtVF79yI+9bppp3PuXeDdGvN+Ve3xBX6uS6TBKqsc0eEe9h0s5w/vfcmrizZz7dDOvLp4C9v2FZOeHMN9ozK4fGCHJqth9rI8Js/NOWnv1xC/nL2S1xZt5d6+5dw66V3GDevEw5f3C3RZchLo7swSEmYvy2PSrJUUl1cenrd1bzGT5+Ucns7bV8ykWSsBmiR8a9bQ1O/XEL+cvZKXF245PF3p3OFphX3o0xAIEhImz805IuSPpbi8kslzc+pczl81NOX7NcRri7Y2aL6EFgW9hIRt+4rrvWxeA5b1Rw0Nqa2pVB7jepljzZfQoqCXkJCeHFPvZQ14cv4GSurxH4A/amhIbU3FY7VfUHas+RJaFPQSEu4blUFMhOeIeRFhRkSNK2ajwsPo1zGJyXNzuOgvH/Pequ3UdXX4idQQE+HhvlEZfln/iRg3rFOD5kto0cFYCQmHDnbWPOOltnmXD+zApxt28et3VnPHy0sZ1i2FB77bm74dTuzK2mPVEOgDsfDtAddDffIeM51104LUOdZNU8nIyHA5OYE/SNVUQv0y7FBoX0VlFa9nb+Uv769jz8EyLh/QgR9f2JONKxYHfduOJxS23fGEevvMrEnGuhEJSeGeMG4Y3oUxA9J5Omsj0z75mn+t+IbzOobRb0gpreOjAl2iiF+oj15avMToCH4+uhdZ943kykEdmLe5grP/OJ9H5+awv7g80OWJnDDt0Yv4tE+K4fff60//qF18VpDME/M3MH3BJm4Z0Y1bzupGUkxEoEsUaRTt0YvUkB4fxhPXDeLde87mzFNa89cP13PW7z/iT/Ny2HOgLNDliTSYgl7kGHqnJ/LsjUN4956zOatHKo9/tIERv/+IX7+zullcBCVSX+q6EalD7/REnr5hMOt3FPL0fzby0oLNzFiwme/2b8/3z+nutxueiDQVBb1IPfVom8CfrxnAjy/syQufbuL1xVuYvXwbZ3Rvzc0junL+aW3xhOlKU2l+FPQiDdSxVSwPfLc395zfg9cWb+GlzzZx24yldEqJ4YZhXbhmSCdaxUUGukyRw9RHL9JISTER3HHuKXz8s/N46vpBtE+K4Xf//pJhv/uQH89czpJNe/w2vILIidAevcgJCveEcUm/9lzSrz052wuZsXATs5dtY9bneZzaJp5rh3bi8oEdSNUFWBIg2qMX8aOMdgk8fHk/Fv3f+fzhe/2Ijwrn4X+tZfj/+5Dvv7SE91Z9Q2mFf0fNFKmL9uhFmkBcVDhjh3Zm7NDOrNtRyD+W5jJrWR7vr9lBUkwE/9O/PWNOT2do1xQdwJUmp6AXaWI92yYw6ZLTuG9UBp9u3M3sZXm89Xkery7aQtvEKC7p157/6deeQZ1bEabQlyagoBc5ScI9YZzbM41ze6bxyBUVfLg2n3e+2MYrC7fwwqebaJMQxag+7RjVpx3DuqcQ4VHPqviHgl4kAGIjw7n09HQuPT2dwpJyPvoyn3dXfsPfl25lxsLNJEaHMzKjDeef1oZze6aRHKvTNaXxFPQiAZYQHcFlAzpw2YAOFJdV8smGXcxbvZ35OfnM+WIbYQaDOrfy/jeQkUbf9CR18UiDKOhFmpGYSA8X9m7Lhb3bUlXl+CJ3Hx99mc9/1u3kT++v40/vr6NVbARnnprKWaemckb31nRpHYvp3q9yHAp6kWYqLMwY2LkVAzu34idg1TZjAAANsElEQVQXZbCrqJRP1u/iv+t38cmGnfxrxTcAtE+K5ozurcnslkJmtxS6pcYp+OUICnqRIJEaH8XlAztw+cAOOOfYuPMAC77azYKNu/jPup3MWpZ3eLnBXZIZ0iWFQV2S6ZOeRHSNm5ZLy6KgFwlCZsapbeI5tU08Nw7vcjj4F329m6Wb9rJk817mrt4BQITH6N0+kf4dk+nfMYmywioqq5zO329BFPQiIaB68F8/rAsA+QUlLNu6j+Vb97Fsy15mfZ7LjIWbAXh48VxOa59A3w5J9G6fSO/0RHq2TdCef4hS0IuEqDaJ0YfPyweorHJ8tbOImR8spCIxnVV5+/nH0lxeKvMOyRBm0DU1jl7tEujRJoGMdgn0bBtPl9ZxOqc/yCnoRVoIT5jRo20CIzpEMHJkHwCqqhxb9x5k9bYCvtxeyJffFLB6WwH/XrWdQwNvhocZXVrHckpaPN3T4umeGke3tDi6to4jNT5SB36DgIJepAULCzO6tI6jS+s4LunX/vD84rJKNu4sYt2OQjbuLGJDfhEbdx5gfk4+5ZXfDr0cHxVO55RYurSOpXNKLJ1SYunYKoZOKbF0SI5RV1AzoaAXkaPERHro2yGJvh2OvE1iRWUVuXuL+Xr3ATbt8n5t3nOQnB2FfLg2n7LKqiOWT42PpENyDOnJMbRPiiE9OZp2SdG0T4qmbWI0bRKiiQxXt1BTU9CLSL2Fe8LomhpH19Q4yDjyuaoqR35hKVv3HmTrnoPk7S0mb5/3a92OQrJydlJcfvQQzSlxkbRJiCItIYo2CdGk+R6nxkeSGh9FanwUKXGRtIqNIFzHChpFQS8ifhEWZrRL8u6xD+2actTzzjkKiiv4pqCYb/aXkF9Qwo6CUrYXlJBfUMrOwhI25Bexq6j0iO6hQ8y8d/VKiYskJTaSVr7wbxUbSXJsJMmxESTFRLBpdyWpeftJiokgKTaC+MjwFj9khIJeRE4KMyMp1hu+vdolHnO5Qx8IO4tK2FVUxq6iUnYXlbH7QBl7DpSy90A5ew6UsXXPQVbklrH3YDllFUd2Gf0x+5Nq7+s9lpAYHUFCdLjvK4L4qHDio8NJiAonPiqcuKhw4qI83u+R3unYSA+xkd55MZEeYiM8QflfRb2C3sxGA38FPMDzzrnf13g+CngJGAzsBsY65zb5t1SRlmH2sjwmz81h275i0pNjuG9UBpcP7MD1zy3g0417Di834pQUrh7SudZla1vHks17eG3RVu7tW86tk95l3LBOPHx5vwbVcKz5DVnHL2ev5LVFW6l0Do/ZUXVU/0BYlXf0OoAj5v3hyn6M6tuO/cXl7DtYTtZn2XTL6ENBcTkFJeVkb9rDJ+t3kbevmOjwMDqmxFJcXsnGnRUUlVRQVFpBaY0PiuOJ8BgxER5iIj3ERHiI9j2ODvcQHRFGtG9edEQYUeEeoiLCiPZ9jwr3EBkeRpTvK9IT5pv2zo88PM+I8D0X4QkjIiyMCN+8xqgz6M3MAzwJXAjkAtlmNsc5t6baYrcCe51zp5rZtcAfgLGNqkikBZu9LI9Js1Ye7svO21fMpFkreXL+etbnHzhi2U837jki+A8tu2TzHv6xNO+Idfz4jeVUj7JK53h54RaAo8L+WDXUtt5Js1YCHBX2x1rH35dsOaLmhtZx39+/AONw107evmL+761VmBmXD+xA+6QYdrT2MLJvu8Pr+HjdrsPrKKmoIm9vMb+7st8RNZdXVnGgtIIDZZXe76UVHCit5EBZBcVllRwsq+RgWYXveyXFZRWUlFdxsLySEt9XcVklhaXllJRXUVJeSWmF73t51VEHqU+2+uzRZwIbnHNfAZjZ68BlQPWgvwx4yPf4TeAJMzPn3NEdbSJyTJPn5hx1wLK4vPKokD+W4vLKw3vL1R0rZl5btPWogD1WDbWtt7i8kslzc44K+mOto3rIN6aO8qqjI+VYNRyvjprLR3jCfP38tZZ3wqqqHGWVVZSWV1Fa+W34l1VUUVpRRbnvcVnFt/PLKqqoqKqirNJ5H1d6l7v7Dw1/f6sri83sKmC0c26ib/pGYJhz7q5qy6zyLZPrm97oW2ZXjXXdBtzmm+wLrGp4yUEjFdhV51LBK5TbF7C2RbY7dXBTv0flwf14Yr89bbJs+4alJ1pDc1hHtdcf3n7HW0fN9wsiGc65hIa8oD579LUdrq756VCfZXDOTQGmAJjZEufckHq8f1BS+4JXKLcNvO2r2J8f0u0L9e3X0NfUp2c/F+hUbbojsO1Yy5hZOJAE1P4/moiInFT1CfpsoIeZdTOzSOBaYE6NZeYAE3yPrwI+Uv+8iEjzUGfXjXOuwszuAubiPb1ymnNutZn9BljinJsDTAVmmNkGvHvy19bjvaecQN3BQO0LXqHcNlD7gl2D21fnwVgREQluwXeJl4iINIiCXkQkxJ2UoDezaDNbbGZfmNlqM/u1b343M1tkZuvN7A3fwd6gZGYeM1tmZv/0TYdS2zaZ2UozW37o1C4zSzGz933te9/MWgW6zsYys2Qze9PMvjSztWZ2Rqi0z8wyfNvt0FeBmd0bQu37X1+mrDKz13xZE0p/ez/ytW21md3rm9fgbXey9uhLge84504HBgCjzWw43qES/uKc6wHsxTuUQrD6EbC22nQotQ3gPOfcgGrnJ98PfOhr34e+6WD1V+A951wv4HS82zEk2uecy/FttwF4x6I6CLxFCLTPzDoA9wBDnHN98Z4scmgIlqD/2zOzvsD38Y5OcDrwXTPrQWO2nXPupH4BscDnwDC8V6+F++afAcw92fX4qU0dfT/w7wD/xHsBWUi0zVf/JiC1xrwcoL3vcXsgJ9B1NrJticDX+E5MCLX21WjTRcCnodI+oAOwFUjBewbhP4FRofK3B1yNdxDJQ9MPAD9rzLY7aX30vq6N5UA+8D6wEdjnnKvwLZKLd8MFo8fwboBDQ4q0JnTaBt6rnOeZ2VLfMBYAbZ1z3wD4vrcJWHUnpjuwE3jB1/X2vJnFETrtq+5a4DXf46Bvn3MuD3gU2AJ8A+wHlhI6f3urgHPMrLWZxQKX4L0wtcHb7qQFvXOu0nn/feyI91+R02pb7GTV4y9m9l0g3zlXfdyMeg0JEURGOOcGARcDd5rZOYEuyI/CgUHA0865gcABgrAboy6+fuoxwN8DXYu/+PqmLwO6AelAHN7f0ZqC8m/PObcWbzfU+8B7wBdAxXFfdAwn/awb59w+IAsYDiT7hkyA2odWCAYjgDFmtgl4HW/3zWOERtsAcM5t833Px9u/mwnsMLP2AL7v+YGr8ITkArnOuUW+6TfxBn+otO+Qi4HPnXM7fNOh0L4LgK+dczudc+XALOBMQutvb6pzbpBz7hy8F6OupxHb7mSddZNmZsm+xzF4N9BaYD7eIRPAO4TC2yejHn9yzk1yznV0znXF+6/xR8656wmBtgGYWZyZJRx6jLefdxVHDnsRtO1zzm0HtprZoTugno93CO6QaF814/i22wZCo31bgOFmFmtmxrfbLiT+9gDMrI3ve2fgSrzbsMHb7qRcGWtm/YHpeI+KhwEznXO/MbPuePeCU4BlwA3OudImL6iJmNlI4KfOue+GStt87XjLNxkOvOqce8TMWgMzgc54/+Cuds4F5UB2ZjYAeB6IBL4Cbsb3e0potC8W70HL7s65/b55IbH9fKdqj8XbpbEMmIi3Tz7o//YAzOy/eI/5lQM/ds592JhtpyEQRERCnK6MFREJcQp6EZEQp6AXEQlxCnoRkRCnoBcRCXH1uTm4yEnlO33sQ99kO6AS7zAFAJnOubKAFHYcZnYL8K7vvHyRZkWnV0qzZmYPAUXOuUebQS0e51zlMZ77BLjLObe8AesLrzYmi0iTUdeNBBUzm2DeexssN7OnzCzMzMLNbJ+ZTTazz81srpkNM7P/mNlXZnaJ77UTzewt3/M5ZvbLeq73YTNbDGSa2a/NLNs3Rvgz5jUW7/Dbb/heH2lmudWuBh9uZh/4Hj9sZs+a2ft4B1ILN7M/+957hZlNPPk/VQl1CnoJGr7xua8AzvQNkBfOtzeiTwLm+QZfKwMewntJ/NXAb6qtJtP3mkHAdWY2oB7r/dw5l+mcWwD81Tk3FOjne260c+4NYDkw1nnHfq+ra2kgcKlz7kbgNryD4mUCQ/EOGte5MT8fkWNRH70EkwvwhuES79AmxOC9tB+g2Dn3vu/xSmC/c67CzFYCXautY65zbi+Amc0GzsL7d3Cs9Zbx7RAQAOeb2X1ANJCKd1jcfzewHW8750p8jy8CTjOz6h8sPfBe2i7iFwp6CSYGTHPOPXDETO9IhdX3oqvw3tXs0OPqv+c1D0q5OtZb7HwHsnxjxjwBDHLO5ZnZw3gDvzYVfPsfc81lDtRo0w+dcx8i0kTUdSPB5APgGjNLBe/ZOY3o5rjIvPeIjcU7lvmnDVhvDN4Pjl2+ET2/V+25QiCh2vQmvLfuo8ZyNc0FfnhoWF3z3uM1poFtEjku7dFL0HDOrfSNVviBmYXhHdHvDho23vgnwKvAKcCMQ2fJ1Ge9zrndZjYd7zDNm4FF1Z5+AXjezIrxHgd4CHjOzLYDi49Tz7N4RyFc7us2ysf7ASTiNzq9UloM3xktfZ1z9wa6FpGTSV03IiIhTnv0IiIhTnv0IiIhTkEvIhLiFPQiIiFOQS8iEuIU9CIiIe7/A4f0FfN7o3tgAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "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)\n", + "print(data_pred)" + ] + }, + { + "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": 33, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\rollandm\\AppData\\Local\\Continuum\\anaconda3\\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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VPW9+P/XObMlk30lyI5Cwq7iAshSqRB2FGnr8pW6lNb2trT+7m29FbvbWlvvtbXXWqmtVoUWa6ssakClVFlURGUPO7KF7HtmO+d8fn9MMhAJMQmZTGbyfj4ekMxyzrw/Sea857NrSimFEEIIcQF6pAMQQgjRvUmiEEII0SpJFEIIIVoliUIIIUSrJFEIIYRolSQKIYQQrZJEIYQQolWSKIQQQrQq7Imirq6OOXPmcPLkyfMe27dvHwsWLCA/P5+lS5diGEa4wxFCCNFOYU0UO3bs4NZbb+XYsWMtPv7d736XH/7wh6xbtw6lFC+++GI4wxFCCNEBYU0UL774Ij/60Y/Izs4+77FTp07h9Xq5/PLLAViwYAEFBQXhDEcIIUQH2MN58p///OcXfKykpISsrKzQ7aysLIqLi8MZjhBCiA6IWGe2ZVlomha6rZRqdlsIIUT3ENYaRWtycnIoLS0N3S4rK2uxiao1R09WYAQsADQt+A9NQ0ML3m76qp39igZ60+2m42h8rn72GIBIrqubkZFIeXld5AIIMylf9IrlskFsl0/XNdLSEtp9XMQSRZ8+fXC5XGzfvp2xY8eyatUqJk+e3K5z+P0m/sZE0VFNSUEjmDE0QNPBpunYdA1N17DpGrquoWsaNg20xu81LbzJxLJiewV4KV/0iuWyQeyXr726PFEsXryYJUuWMGrUKB599FEefPBB6urqGDFiBIsWLerqcEIXeoWCpr8NCwzM856rNf6naRo6oNt07LqGza5jt2nYtGBCselaRGsjQgjRmbRo3rho/5HSi65RdKam5i1dB6fNht2uY7c1JpF2Jo+srCRKS2vDF2yESfmiVyyXDWK7fLqukZGR2O7jItb0FIuUCnbKWxbByYO+4P26pqHbwGm34bDrOGw6Nl1H1yPbDyKEEG0hiaILWEphGYRmnjfVPJw2HYdDx2GzYbdLk5UQonuSRBEBTTUPr2XiDZhAINhRbtNwOmy47DYM0wp7Z7kQQrSFJIpuIljrUAQMi3oCaE4b9TU+XC4bcU5bcFSWEEJEgCSKbkop8AaCNY56XSM+zk6c045Nl4QhhOhakiiigGEpahsC1HsN4l12EuLs6DKLXQjRRSRRRBHLUtR7Anh8BglxdtxxdmmSEkKEnWxcFIWsxhpGeY0Xf8BEKhdCiHCSRBHFDENRWeejqs6PKUsOCCHCRBJFlFMKPD6DihovPqldCCHCQBJFjDAtRVWdj5qGQKRDEULEGEkUMUQpqPcEqKj1YanuswaWECK6SaKIQf6ASUW1D78hyUIIcfEkUcQow1JU1foalwgRQoiOk0QRwyylqK7zUe8zpJNbCNFhkihinFJQ1+CnzivJQgjRMZIoeoBQsvAEJFkIIdpNlvDoIYLJIjh0NiHOEeFohBDRRGoUPYgC6jwBPH4j0qEIIaKIJIoeRimoqffjl9FQQog2kkTRAykFVfV+mWchhGgTSRQ9lGUpaup8mJYkCyFE6yRR9GBG4/pQlmzMLYRohSSKHi5gKKrr/CgkWQghWiaJQuALmFTX+yVVCCFaJIlCAOD1mdR5ZIlyIcT5JFGIkAZvAJ8MmxVCfIokChHSNMdCRkIJIc4liUI0Y1qK6nppghJCnCWJQpzHHzCpafDLAoJCCCDKE8Urbx+hpt4f6TBiksdn4PNLf4UQIsoTxYGT1Tz24g7e21ssk8Y6WVN/hfxchRBRnSjcLju+gMmqTUd5eu1eymu8kQ4pphiWorZB9rAQoqeL6kRx9+xhXDEkE4BjRbU8/tJOtuw+I5+CO5HXZ+CRJigherSoThTuODtfuP4yvjwjl+QEJwHDYu2WY/xp7V4qa6V20RkUUCtNUEL0aGFNFGvWrGHWrFlMnz6d5cuXn/f4nj17uPnmm5k3bx5f+9rXqKmp6dDr5PZP49sLRzN2aBYAR4tqefylXXx4oBQlF7iLZlpKZm0L0YOFLVEUFxfz2GOPsWLFCl555RVWrlzJoUOHmj3n5z//OUuWLGH16tUMGjSIP/3pTx1+vXiXnZs/dymLZuSSGO/AFzB5aeNhVrx5kAavXOQulsdnEDBlIp4QPVHYEsWWLVsYN24cqampuN1u8vPzKSgoaPYcy7Kor68HwOPxEBcXd9Gvm9c/jSULRzN8YBoAe45W8Lt/7OLI6Y7VVkTQuXtuCyF6Fk2FqW3mqaeeoqGhgfvuuw+Av//97+zcuZOf/exnoed8/PHH3H333bjdbuLj43nxxRdJS0tr82ucKa/HtFoOXynF1l1FrHzjAL6AiabBrAmDmHXdQGx6VHfNRFRygpMktzPSYQghupA9XCe2LAvtnHGVSqlmt71eL0uXLuXZZ59l9OjRPPPMM9x///0sW7asza9RXd2AP3Dh5pC8vin8x00j+dtbBzld3sCrm4+y+3AZt3x+CCkJ3ftil56eQEVFfaTDOE9NVQPpKXHoFzlmNisridLS2k6KqvuJ5fLFctkgtsun6xoZGYntPy4MsQCQk5NDaWlp6HZpaSnZ2dmh2wcOHMDlcjF69GgAvvSlL/H+++93ehyZqfHce+NIJo7qDcAnZ2r5v3/s5ODJqk5/rZ7AsBT1XiPSYQghulDYEsWECRPYunUrFRUVeDwe1q9fz+TJk0OPDxgwgDNnznDkyBEA3nrrLUaNGhWWWOw2nVnjB7AoP5d4l416r8GzrxXy5gcnsC7QdCUurMEbwG9Ix7YQPUXYmp569erFfffdx6JFiwgEAixcuJDRo0ezePFilixZwqhRo3j44Yf5zne+g1KKjIwMfvGLX4QrHADyBqTxzQWj+eubBzhZWs+GD09xoqSOL00dgjsubD+KmKMU1Db4yUi++MEHQojuL2yd2V1h/5HSVvsoLsQwLV5/9zhb95wBIC3Jxe3ThnJJZkJnh9hh3bWP4lyJbgeJcY4OHRvL7cAQ2+WL5bJBbJev2/VRdGd2m87c6wbyhesvxWHTqaz18YdVu/n4YFmkQ4sq9Z6AzK0QogfokYmiyRVDsrj3xhGkJbkwTMWL/zpEwXufSL9FGwWboAIEF/oQQsSqHp0oAHpnJPAfN43isj4pALy9o4jn1+3H65eRPW0RMEz8AUkUQsSyHp8oILi44Jdn5nHdyBwA9p+o4slXdsuy5W2gFNR7pVYhRCyTRNHIpmvMnjCQm6cMxqZrlFZ5efLl3Rw7I0t/fBa/YeI3JFEIEaskUXzK2Nxs7pkzDHecnQafwZ/W7uOjg6WffWAPphQ0eGSDIyFilSSKFgzMSeYbN44kKzUO01L8/V+HeWv7SVmyvBU+w8TXgaHKQojuTxLFBaQnx3Hv/JGhTu63tp/kn/8+gmnJxbAlUqsQInZJomhFvMvOl2fmcmXjhkjbD5Tyl9dlRNSFSK1CiNgkieIz2HSdm6cMZuqVfQA4dKqaP67ZS22DP8KRdT9KBdeBklqFELFFEkUbaJrGDVf14+Ypg9E1KCpv4A+r9lBeLcNnP80XMGXBQCFijCSKdhibm83/y89ttuzHqdK6SIfVrQRrFYbUKoSIIZIo2imvfxr3zBlGvMtOvdfgj2v3cvh0daTD6la8foOAzKsQImZIouiA/r2S+Nq8EaQkOPEHLP7yeiH7jlVEOqxuo2m2ttQqhIgNkig6KDstnq/NH0FmShyGqVj+xgGZmHcOr9+QlWWFiBGSKC5CaqKLr84bwSUZbiwFf//XYd5t3OOip1MKPD4z0mEIITqBJIqLlBjv4CtzhzMwJwmA1ZuP8c6O0xGOqnvw+Awsmc0uRNSTRNEJ4px27pyVF5rF/fp7x2XJD8CyFN6A1CqEiHaSKDqJ027jjvxc8vqnAcElP9a9f6LHJwuPV2axCxHtJFF0Iodd5/bpQxg1OB2At3ec5vX3jvfoZGGYlkzAEyLKSaLoZDZd50tTh3DFkEwANu0s4tWtn/TYZKEUsjaWEFFOEkUY6LrGzVMuDS0muGX3GdZsPtZjk4XXb2LKPuRCRC1JFGGi6xoLpgzmqrxsAN7dW9xjk4VlKXzSqS1E1JJEEUa6pnHjpEFcM+xssli7pWc2QzX4ArKrthBRShJFmOmaxryJg7i6sWaxdc+ZHtlnYRqKgOxVIURUkkTRBXRNY/6kQVyVe7bP4vV3e9ZoKAV4pFNbiKgkiaKL6JrGjZMHM7YxWWzaVcT6bT1rnoUvIJ3aQkQjSRRdSNc0bpo0ODR09t8fn2bDh6ciHFXXkU5tIaKTJIouFhwNdWloUt5b20/y7497TrLw+AKRDkEI0U6SKCLApmt8ceplDB8YXO5j3fsn2LK7Z6w6a5hKZmoLEWUkUUSITde55fNDGNovFYC1W47xQWFJhKMKP5mpLUT0kUQRQXabzu3ThjKodzIAL799hB2HyiIcVfh5/SYBqVUIETUkUUSYw66zKD+XftmJKIKbH+37pDLSYYWVZSm80lchRNRoU6J4/vnnqaurC3csPZbLaePOmXn0znBjKcVf3zzA/hhPFh6fKXtqCxEl2pQo9u/fT35+PkuXLmXXrl1tPvmaNWuYNWsW06dPZ/ny5ec9fuTIEe644w7mzZvHPffcQ3V1ddsjjzHxLjt3zRoW2oP79//YwcnS2E3OhmlJ85MQUaJNieKhhx5i3bp1jBw5kp/85CfcfPPNvPTSS/h8vgseU1xczGOPPcaKFSt45ZVXWLlyJYcOHQo9rpTi61//OosXL2b16tUMGzaMZcuWXXyJolhivIO7Zw8jJcGJz2/y7GuFFFc2RDqssFAgu98JESXa3EeRmJjIjBkzmDNnDlVVVaxYsYIZM2awYcOGFp+/ZcsWxo0bR2pqKm63m/z8fAoKCkKP79mzB7fbzeTJkwG49957uf322y+yONEvNdHF3bOHkeR20OAzeObVfVTWeiMdVlh4vQY9aGK6EFGrTYli69atfOc732HGjBkcOXKEJ554gn/+85/85S9/4Yc//GGLx5SUlJCVlRW6nZ2dTXFxcej28ePHyczM5IEHHuCmm27iRz/6EW63+yKLExuyUuP51hevwOWwUdMQ4M+vFVLnib3OX8NSBAypVQjR3dnb8qSf/OQn3HbbbfzsZz8jKSkpdH///v354he/2OIxlmWhndNbqZRqdtswDN5//31eeOEFRo0axW9+8xt++ctf8stf/rLNwaekuGN27aB04JtfGMPjL35MebWX59cf4P+77UriXW36lUWF9PQEXHYbmWnxkQ4lLLKykj77SVEqlssGsV++9mrTVWf16tUUFBSQlJREaWkpr776KosWLULXdZYsWdLiMTk5OXzwwQeh26WlpWRnZ4duZ2VlMWDAAEaNGgXAnDlzLniuC6mubsAfo0tXp6cnkJHo5Japl7H8jQOcKK7l8b99yJ0zh+GwR/+o5vT0BCoq6tF1DcMfwKbH1hCorKwkSktrIx1GWMRy2SC2y6frGhkZie0/ri1P+tnPfsbGjRsbX0hn+/bt/OIXv2j1mAkTJrB161YqKirweDysX78+1B8BcMUVV1BRUUFhYSEAGzZsYMSIEe0uQKwbNjCdBVMuBeBoUS0rNxyMqVqUZSl80vwkRLfWphrFRx99xNq1awHIyMjgt7/9LfPnz2/1mF69enHfffexaNEiAoEACxcuZPTo0SxevJglS5YwatQonnjiCR588EE8Hg85OTn86le/uvgSxaArh2bR4DV47d1P2HusklWbjnLTpEHNmvKimcdr4HbagNgojxCxpk2JIhAI4Pf7cTqdQLB/oS3mzp3L3Llzm933xz/+MfT9mDFjeOmll9oaa482cXRv6jwB3t5xmg8KS0iMdzD96n6RDqtTGKZFwFQ4bJIohOiO2pQoPve5z3HPPfcwf/58NE1j7dq1TJkyJdyxiU/Jv6Yf9Z4A2w+UsvGjUyTGO5gwMifSYV204EKBJo746O97ESIWtSlRfO9732P58uW89dZb2O12pk2bxi233BLu2MSnaI275NV7DQqPV/LqlmMkxjsYfWlGpEO7aD6fQVK8I9JhCCFaoKko3otz/5HSmB71VFFR3+JjfsPkz6/u43hxHTZd486ZeVzaJ6WLI7w4LZUvLdmFy26LUESdK5ZHzsRy2SC2yxfWUU9vvvkmU6dOZezYsVx55ZWhfyIynHYbi/LzyE6Lx7QUL6w/wOmylpNKNPF4DFkoUIhuqE1NT7/+9a/57//+b4YPHx4zI22inTvOzp0z83hq1R6q6/385fVCvjZ/BOnJcZEOrcP8hoVhWth06asQojtp0zsyOTmZ6dOn07dvX/r06RP6JyIrNdHFnTPziHPaqPUEeOb1Quq90bvUh6UUvhhtShQimrUpUYwZM4Z///vf4Y5FdECvdDeLZuRit2mUV3t5rmA//ihelbVBNjQSottpU9PTv//9b1544QUcDgcOhyO0btOHH34Y7vhEGwzMSeaLU4fw1zcOcKKkjr++eZD/l58blctimKbCb1g4Y2CZEiFiRZsSxbPPPhvmMMTFGjkonbkTB7J60zH2n6hi1TtHuGny4KjrU1IKfAFTEoUQ3Uib3o19+vRh165dvPjii6Snp/PRRx9JH0U3NG54Dp+7Ivh7+WB/KW9tPxnhiDrG52vbzH8hRNdoU6JYtmwZf/3rXykoKMDr9fJ///d/PPHEE+GOTXTAtKv6cuXQTAA2fHiKbfuKP+OI7sewVFT3swgRa9qUKF599VX++Mc/Eh8fT1paGi+++GJokUDRvWiaxk2TBzOkb3AC3iubjlL4SWWEo2o/j18ShRDdRZsShd1uDy0ICMHhsnZ77GygE2tsus5t04bSJzMBpeCvbx7kREl0zTT1BUys6F00QIiY0qZE0bt3bzZu3Iimafj9fp588knpo+jmXA4bi2bkkp7kImBa/KVgP2XVnkiH1WaWFRz9JISIvDYlih/84Ac888wz7N+/n8svv5y3336bH/zgB+GOTVykJLeTO2fm4XbZafAaPPtaIbUN/kiH1WYeryzpIUR30K5FAT0eD6ZpkpjY/kWlwqGnLgrYXseLa/nT2n0ETIs+WQl8Zc5wXI7ILr7XlvLpmkZGsgubLfqGysbywnKxXDaI7fJ1dFHANnU0PPPMMy3ef9ddd7X7BUXX698riVtuGMIL6/dzqrSev755kDvyh3b7NZUspfAZFu4oTBRCxJI2JYoDBw6Evvf7/Wzbto3x48eHLSjR+YYNSGPedYNYtekoB05U8co7R1kQBRPyPD6DhDg70q8tROS0KVE8/PDDzW4XFxezdOnSsAQkwufa4b2oqffzr49OsX1/KSkJTm64qntvpxrcJtXC3s1rP0LEsg69+3r16sWpU6c6OxbRBW64qi9XDs0CGifkFZZEOKLWBZf0iM1+KCGiRbv7KJRS7N69m4yM6N9+sycKTsgbRG2Dn4Mnq1n1zhGS4h3kDUiLdGgX5PUFSHDZgO7dTCZErGpTjeLAgQOhfwcPHqR37948+uij4Y5NhMm5E/KsKJiQZ5gKw5ROCiEiRfbM7qY6c3jshdQ2+PnDqj1U1vpwx9m5d/4IMlPiw/qaTdpbviS3g4Q4Rxgj6lyxPMQylssGsV2+sA6PveOOO1odHfPcc8+1+4VF5CW5ndw1K48/rNpDg9fgmdcKuXf+CJLczs8+uIt5/EZUJQohYkmbmp5GjhxJXFwcixYt4p577iEzM5PU1FRuv/12br/99nDHKMIoMyWeL8/IxWHTqaz18ZeC/fi64YJ8pjQ/CRExbapRfPjhh6xYsQKbLTibd9KkSXzxi18kPz8/rMGJrtEvO4lbGyfknS6rZ8WbB7gjPxd7N5roFhz9ZGC3Sa1CiK7WpitBRUUFPp8vdLu+vh6v1xu2oETXyxuQxvxJgwE4eLKal98+QnfrvvL4DLpXREL0DG2qUcyZM4cvfelLTJs2DaUUr7/+OosWLQp3bKKLXZ2XTU29n7e2n+Sjg2UkuR3MuHZApMMKMS1FQPbTFqLLtSlRfPvb32b48OG8++67uFwufvrTn3LNNdeEOzYRAVOv7ENtg5/395Xw9o4iktxOrhvVO9JhAcHmJ6/fwGnvfp3tQsSyNn8069WrF0OGDOE73/kODoe0E8cqTdOYd90ghg8MTsB7desn7DhUFuGozvL6TZQ0QAnRpdqUKP7xj3/w/e9/n6effpra2lq+8Y1v8OKLL4Y7NhEhuq7xpalDGJCTBMBLGw9z6GR1hKMKsiwlS3oI0cXalCheeOEFVq5cSWJiIhkZGfzzn//kL3/5S7hjExHksOssys8lOy0e01K88MZ+TpXWRTosALyyoZEQXapNiULX9WabFfXu3Ts0VFbErniXnbtm5pGS4MQfsHj29cJusZ2q37AwTKlVCNFV2pQoUlNT2bdvX2h29urVq0lJSQlrYKJ7SEl0cdesYcS77NQ3zt6uifB2qpaS5ichulKbEsUDDzzAd7/7XQ4fPszEiRP57W9/y4MPPhju2EQ3kZ0Wz50zc3HYg7O3n32tEI/PiGhMwdeXTm0hukKbEoXX62XVqlW8/PLL/PnPf6agoIDc3NzPPG7NmjXMmjWL6dOns3z58gs+b+PGjUydOrXtUYsu1y87idunDUXXNM5UNPDcuv34jcgt9WGYFn5DEoUQXaFNieK//uu/sNlsXHrppQwdOrRNw2OLi4t57LHHWLFiBa+88gorV67k0KFD5z2vrKyMRx55pP2Riy43tF8qC6+/FIBPztTytzcPYlqRaQJSCuo8gYi8thA9TZsSRW5uLmvWrOH06dNUVVWF/rVmy5YtjBs3jtTUVNxuN/n5+RQUFJz3vAcffJBvfvObHYtedLnLL8tkzoSBABQer+IfG49gRWipD3/AxBvofgsYChFr2jQz+6233jrvIq9pGvv27bvgMSUlJWRlZYVuZ2dns3PnzmbPee655xg+fDhjxoxpT8whKSluTCt2mx/S0xMiHUKL5ky+FKVpvLr5KB8fKiMtJY4v3jC01aXoW9IZ5bPpkJ6WgE3vfuNls7KSIh1C2MRy2SD2y9debUoUu3btaveJLctqduFQSjW7feDAAdavX8+zzz7LmTNn2n1+gOrqBtm4KEImDM+mrLKB9/YW86/tJ9GBz4/t2+bjO7N8Xk+ABFeb/pS7TCxvfhPLZYPYLl9HNy5qtenpBz/4Qej7ioqKdp04JyeH0tLS0O3S0lKys7NDtwsKCigtLeXmm2/mq1/9KiUlJdx2223teg0ROZqmMfe6gYy5LLh3+lvbT7J5V1FEYqn3BCLW/CVET9Bqoti9e3fo+3vuuaddJ54wYQJbt26loqICj8fD+vXrmTx5cujxJUuWsG7dOlatWsWyZcvIzs5mxYoV7QxfRJKuaSz83KXk9U8FgutCfXig9DOO6nyWpaj3Rna4rhCxrNVEce5+BO3dm6BXr17cd999LFq0iBtvvJE5c+YwevRoFi9e3KGmLNE92XSdW28YyqDewTbdf/z7MLuOlHd5HB6vgRWhEVhCxLo2N+y2t6MSYO7cucydO7fZfX/84x/Pe17fvn3ZsGFDu88vugeHXeeO/Fz+/Oo+TpbW8+KGQzjtOrn907osBksp6rwGyd1wv28hol2rNQrLsqiurqaqqgrTNEPft2V4rOhZ4px27pyZR6/GRQSXv3GAI6e7dsVZj8/AkFqFEJ1OU620KeXl5aFpWovNTp81PLYr7D9SKqOeupnaBj/L1uylvNqL06Fz96xh9O91/lDDcJUv3mUnNdFJpPu2Y3nkTCyXDWK7fB0d9dRq01NhYWGHAxI9U5LbyT2zh7Fs9R6q6vw8+3oh98weRp+s9v9xdoTXbxAwHNht3W9ehRDRSjYfFp0uNdHFPbOHk+x24PWb/Pm1QorKu6Z2FFzawy/7VQjRiSRRiLDISInj7jnDSYh34PEZ/PnVfZRUds1eFr6AKcuQC9GJJFGIsMlOjeee2Wf3svjT2r2UVoU/WSgFDZ6A1CqE6CTda90DEXNy0t3cM3sYT6/dS60nwNNr97J4zvCLWudp//FK3tlxmspaH2lJLiaNueS8obg+I1ircNrls1C02nm4jIL3jlNW7SUzJY4Z1/Zn9KWZkQ6rR5J3kQi7SzITuGf2MOKcNmobgsmipLKhQ+faf7yS1ZuPUuMJEOeyU+MJsHrzUfYfr2z2PKWCS3vI5kbRaefhMpa/cYCqej/uODtV9X6Wv3GAnYfLIh1ajySJQnSJPlmJ3D1rGC6HjZqGAP+74sMO7b/9zo7T2Gw6TrsNTdNw2m3YbDrv7Dh93nP9hok/IIkiGhW8dxybTcflCP6eXY7g77ngveORDq1HkkQhukzf7ETunp2Hy2GjqtbH02va32dRWevDYWv+Z+uwBbdo/TSloN4rtYpoVFbtPa/Z0GnXKav2Riiink0ShehS/bKTuHt2HnGuYM3i6TV7KWlHskhLchEwm49oCpgWaUmuFp/vlxFQUSkzJQ6/0fz35jcsMlPiIhRRzyaJQnS5ftlJfOeWK4N9Fp5gsiiuaFufxaQxl2CaFn7DRCmF3zAxTYtJYy5p8fkKqK33yzLkUWbGtf0xTQtfIPh79gWCv+cZ1/aPdGg9kiQKEREDeyc3Dp21UecJ8Me1ezld9tmT8nL7pzHvukEkxzvw+gyS4x3Mu25QqwsQGpaiTobLRpXRl2Zy+7ShpCY4afAapCY4uX3aUBn1FCGtrvXU3claT9GrqXyny+r582v7aPAaxDlt3DUrj37Znb8NpaYFZ4y7HLZOP3dLYnm9oFguG8R2+cKyw50Q4XZJZgKL5w4nKb5xuY9XCzlaVNPpr6OUNEEJ0VGSKETE9Upzs3jecFISnPgCJs++VnjevIjOYFiK2gZpghKivSRRiG4hMyWer84bTkZyHAHT4vl1B9h5uPN3yvP6Dbx+s9PPK0Qsk0Qhuo20pDi+Om84OeluLKVY+dZBthWWdOprSBOUEO0niUJ0K0luJ1+ZM5x+2Yko4OW3j7CQe89/AAAgAElEQVTxo1Pt3rO9NTIKSoj2kUQhuh13nJ27Zw/jsj4pAKzfdoJXt37SqbUAj8/A54/NEXNCdDZJFKJbcjlsLJqRy6jBGQBs2X2GFzccwjA75+KuFNR4fChZ3kOIzySJQnRbdpvOlz5/GeNG9AJg5+Fy/lJQiNdvdMr5DUPR4JOObSE+iyQK0a3pmsbcCQOZdlU/AA6fqmHZ6r1U1/s75fz1noB0bAvxGSRRiG5P0zSuv7IPCz93KbqmcaaigT+8spszbVwfqjWWpaj3dk4NRYhYJYlCRI0rh2bx5Zm5uBw2quv9PLVqDwdOVF30eT1eA9OSjm0hLkQShYgqQ/qm8tV5w0lunMX9XEEh7+0tvqhzWipYq5DhskK0TBKFiDq9MxL4+o0juSTDjaVg1aajvLr1GJbV8b4Gj884b/8DIUSQJAoRlVISnCyeN4JhA4LLi2/edYbn1u3v8IgopaBGZmwL0SJJFCJquRw2bp82lImjewNw4EQVT76yu0N7cQMEDIuaBr/MrBDiUyRRiKim6xqzxg3g5imDsekapVVefv/ybg6e7Fgnt9dnyvIeQnyKJAoRE8bmZrN47nASG/e1ePa1wg6vEdXgCVDvkyGzQjSRRCFiRv9eSXzjppH0zUpAEVwjavkbB/C1c1nxpn22GzppBrgQ0U4ShYgpqYkuFs8dwVW5WQDsPVbJEy/voridk/OaliOv9wZAei1EDyeJQsQch13npsmDuXHSIGy6Rlm1l9+/spuPDpS26zxKQW1DgFqP1CxEzxbWRLFmzRpmzZrF9OnTWb58+XmPv/nmm8yfP5958+bxjW98g+rq6nCGI3oQTdO4ZlgvvjpvBKmJTgKGxd83Hublt48QaOd8iXpPoNPWlhIiGoUtURQXF/PYY4+xYsUKXnnlFVauXMmhQ4dCj9fV1fHjH/+YZcuWsXr1anJzc/nd734XrnBED9UvO5FvLhhFbr9UALYVlvDkK7sprmxfU5THZ1BZ65N5FqJHClui2LJlC+PGjSM1NRW3201+fj4FBQWhxwOBAD/60Y/o1Su4hHRubi5FRUXhCkf0YO44B3fMyCX/mn7oGpypaOD3/9zNtsKSdo2K8gVMKmq9GBcxA1yIaBS2RFFSUkJWVlbodnZ2NsXFZ9fkSUtLY9q0aQB4vV6WLVvGDTfcEK5wRA+naxpTLu/D4rmNTVGmxctvH+Gvbx6kwRto83kMQ1FZ48UfkOU+RM9hD9eJLctCO2fWklKq2e0mtbW1/Md//Ad5eXncdNNN7XqNlBQ3Zgx/uktPT4h0CGEVifKlpyeQOziDF14v5MP9Jew+WsGJ0nq+PHsYwwdltPk8GhDndpCU4Lrgc7Kykjoh4u4plssGsV++9gpbosjJyeGDDz4I3S4tLSU7O7vZc0pKSrjnnnsYN24cDzzwQLtfo7q6IWY/2aWnJ1BRUR/pMMIm0uW7efIgBvZKZO3WY1TX+Xh85ceMH5lD/jX9cNptbTpHRSUkxjtIiLMTTB1nZWUlUVpa2/mBdwOxXDaI7fLpukZGRmL7jwtDLABMmDCBrVu3UlFRgcfjYf369UyePDn0uGma3HvvvcycOZOlS5e2WNsQIlw0TeOqvGy+dfNo+vcKvnG27j7D7/6xi0/OtO0i0TR8tkoWExQxLmw1il69enHfffexaNEiAoEACxcuZPTo0SxevJglS5Zw5swZ9u7di2marFu3DoCRI0fy85//PFwhCXGejOQ4Fs8dwdsfn2bDhycpr/aybPUeJo7uzQ1X9cNh/+zPUl6fiRHwkpzgxOloW21EiGiiqY4shtNN7D9SKk1PUao7lm/rniLWvXcitC9FktvBhBE5HDxZRWWtj7QkF5PGXEJu/7QWj9c0iHfZSYx30Cs7uUPNFzsPl1Hw3nHKqr1kpsQx49r+jL4086LK1VlWbzrC+m0n8QZM4hw2pl/dl3kTB0c6rE4nTU8tHBeGWISIOvuPV7JpZxFJCQ4S4x1AsFlp3bYTnCpvwOmwUeMJsHrzUfYfr2zxHEpBg9egvNpLTZ2P9n4E23m4jOVvHKCq3o87zk5VvZ/lbxxg5+Gyiy3eRVu96QirtxzDFzCx68Ghwqu3HGP1piORDk10AUkUQgDv7DiNzabjcthJTnCSlRof6p72+U1Kq7wYhkLXNd7ZcbrVc5mWotYToKLWgzdgtnmlqIL3jjfGYEPTNFwOGzabTsF7xy+qbJ1h/baTaGjYdA1N04Nf0Vi/7WSkQxNdIGx9FEJEk8paH3Gus2+Hpr6JpmRhKUVVnQ+nXSfQxubOgKGoqvXhsGskxDuDCaCV55dVe3HHNX9LOu06ZdXe9hQlLLx+A5vePHpdo8M7CoroIjUKIYC0JBcBs3kCsNs07DaN7LR44pzBTmq/YVHTEKDgvU/avHx5U8Ior/ZQ7zMuOEIqMyXuvH27/YZFZkpcB0rUueKcdj49ZclSwftF7JNEIQQwacwlmKaF3zBRSuE3TJwOG06nHVMp0pJcJCc4aPpQ/faOIv535cd8eKC0zUNjDVNRW++nrNpLTYMfw1TNdtKbcW1/TNPCFwjG4AuYmKbFjGv7h6HE7TP96r4oFKalUMoKfkUx/eq+kQ5NdAHbj3/84x9HOoiOKq9siNmZ2fHxTjyeti8tEW26W/kyU+LJTImjuLyBuoYAqQlOZo4bwPCBaaH7MpLjmDluAL3S3ZworsPjN9l7rJL9J6rISo0nLensLO3WyqdUcH9uj9/AH7DQdQ2bTScn3U2vtHhOltRRXe8nPcnFgsmDu8Wop9z+aaAUn5ypw28q4hw2Zo3rH5OjnhISXDQ0xOZqwZqm4XY723+cDI/tnrrj8NHOFO3lq6rz8fq7x9l1pDx03/CBacy4pj+ZqfHtLp/dphEf5yDOYTuvL6C7ieXhoxDb5evo8FhpYBSiA1ITXdx6wxAmnMnhtXc/4URJHXuPVVL4SSVjc7NZMHVIu87X1CxVr2vEOW3Eu+w47Xq7h9gKEQ6SKIS4CANykrh3/gh2Halg3fvHqaz1sa2whI8OljFuRC8mj7kkNC+jLSxL0eA18PgMnHYb7jg7TrsNWeFGRJIkCiEukqZpjL40g+ED09hWWMK/PjxFnSfApp1FvLe3mPEjejFxdPsShlLBSW3BCW4a8XF24hw27FLLEBEgiUKITmK36YwfkcPYoVl8fKSCde8ew+MzeXtHEe/uKeba4b24bnRvktvZmWhYitqGAPWagdOuExcXbJay6ZokDdElJFEI0cmcDhszxg9k9KA0tuw+w6adRXj9Ju/sLGLrnjOMzc1m0ujepCe3b36EpRTegIk3YKLrGk6bjtNlw2nTsdn0VifzCXExJFEIESZxTjtTr+zLhJE5bN1dzObdRTR4Dd7bW8z7+4oZOSiDSaN70ze7/aNQLEvhtYJJQ9PApmmN8z5sOBqH2wrRWSRRCBFmcU4711/Zh+tG5bCtsIR3dhZRU+9n15Fydh0pZ2BOEhNG9WbYgLQODY1VCgylMHwGDT4jmDh0DafdhsOh45Aah7hIkiiE6CJOh43rRvXm2uG92HW4nE27iigqb+DYmVqOnaklNdHJtcN7cXVeNu64tnd8f5pSweG2hmmAj/NqHHabhl2XGodoO0kUQnQxu03niqFZXD4kk8Onati8u4gDx6uoqvOz7v0TvLX9JCMHZXDN8GwG9Eq66N0fW6xxaBoOuw2nU8duC9Y6NA3pHBctkkQhRIRomsZlfVO4rG8K5dVe3t1zhg/2l+ILmHx8qIyPD5WRnRbPVbnZXD4ks13Da1sTShx+A48/WOPQNA2HTcfhCCYOu65hl+QhGskSHt1UtC9x8VmkfC3zB0x2Hi7nvX3FnCo9e7yuaeQNSOXKoVkM7ZeKPcyd1U3Jw2nTG5urdOy24H4UmZmxu8QFyBIeLZEahRDdiNNh46q8bK7Ky+ZUWT0fFJaw41AZ3sYFCPceqyTeZWfU4HQuH5JJ/15J6GGYtq0UKNV8ZJWmacHVc+126r0BbLbgXA5da/ynnz1WxBZJFEJ0U30yE+gzcRCzxg1g77EKPjxQyqFT1Xh8Bu/vK+H9fSWkJDgZNTiDUZdm0Dcr4aL7My6kKXFYBGeM1zYEV8bVNNDQQANdB7umY7PrOOw6uhasgeh6MMFIAolekiiE6OYcdp0xl2Uy5rJMahr87DxUzseHyjhdVk91vZ9Nu4rYtKuIlAQnwwelM2JgOgNzktC7YBVapUChQIFlgYEJgeCGThrBWgga2HSw6zZs9mDysGk6uo1QbaQpv0kyCWr6eQR3UVCAxrm/zk//nC78+UBDKdX4AUJ1eM0wSRRCRJFkt5OJo3szcXRvyqo87DxSzs7D5ZRUeqiu97N19xm27j6D22Unt38qeQPSGNI3JSI70SmCtZCmJBLAgHO2eWhqztII1kZsmh6sfdiakomGpgcft1TwhEpTodpNqJajQFmqMSEFazDBczU1l124WazpwmlaCstSKII1JtNSwePRWri4Nl18Wz5n030B08K0rNDjiuD2sZqmoTfWwmhMBBbBSZSmpbBMhakUlmWFdhXUQj8v0DQdXSPU5NiUrK2mnwvBn82nN2vXNHA6dNLTpY9CiB4jMzWeqVf2ZeqVfSmp9LD3WAV7jlZwqqyeBp/BRwfL+OhgGbqmMSAnidx+qQztn0qvtPiwNVG1R9OFHsC0IMD5W8uGwlTnXfc+U7NmMQ1serBPRdM1lHU20RiNF+TQuB6bjcpqL03Xck0HXdODF+DG4zjn/qaLtqYHk4ppBHdItJQKUw2pbVvwtkTTOhaQJAohYkB2WjzZaX343BV9qKz1UXg8uDfGkdM1mJbiaFENR4tqKHj/OEluB5f1SeHSxn8pCe3f8ayrXMyFtlmzGGCYbbvAKoLraoUykwUtXpwvdH8MkkQhRIxJS3IxfkQO40fk4PObHD5dzf7jVRw4UUV1vZ/ahkCotgGQkRLH4N7JDL4kmYE5SaQkuj7jFURPI4lCiBjmctoYPjCd4QPTUUpRWuXl0KlqDp2s5khRNf6ARXm1l/JqL9sKSwBITXQyMCeZ/jmJ9M9Oole6u9tvzyrCSxKFED2EpmmNTVTxTBiZg2lZnCqt52hRDUdO1/DJmVr8hkVVnT80MxzAadfpk5VA36xE+mYnMkLX0UIjaURPIIlCRJQW+u9Cj7fxYtSBa1bTIeeOJgmNGvmMCf9NI3ZsmoYCTGWhrGC7uBUlYzxtuk7/Xkn075XElMv7YFqKovJ6PjlTy7GiWo6X1FLbEMBvWBwtquVoUeNs5TcP4nbZuSQzgUsyE+id4aZ3RgKZKXFdMiRXdD1JFDGstYtw6AKsffr+xq8aoVEfNA7n05pGd2gaWuPQw7P3gaYah+udM4QxeLh29pyApoLt6Bhm6P5zw2j+QbXlArTlw+wFn6O0syNatNAYxfOf1sr1vqVzKxUcZukzzODomgueuXuy6Vqw1pCVyHWjeqOUorrez/HiWk6W1HOitI7TpfUETIsGnxFswjpVHTreYdPJTo8nJ81Nr3Q3vdLjyU6NJznBKbWPKNejE8W5f7qf9YbWzvvm049r59644PGh98s5F2ENDf2cJ2loxDlsuF320IQl7ZznNz01dOE950IcHP129sr76Yvw2VNo5xzPpy6eLZexMz8ou+Mc1NsjtNR1s4lLHSvUhQ6z6Rpup53MdDem38BvBPe9Nk0VGvsfLTRNIzXRRWqii9GXZgKNidBSFB4p51RZPadL6zlT0YAvYBIwg01Z565RBeBy2MhKjSMrNZ7MlHgyG79PT3bhtNsiUTTRTlGdKOKcdhw2df5FsXEST/NPjS1fSM8TmgnZpPlzW7vwnv3+7BhupfjMC3Dopc952YzUeKyA0foBnamVWZ+i/TRNw2nXcdp1kuIdmJbCMBWWUpimRcCw8BsWlhVdP2ybrtE3MxG3XefKoVlAsKmtqtZHUXkDZyoaKK5soLiigfJqL5YKTmA7WVrPydLzF0lMTnCSkewiPTmO9KQ40pNdpCUF/yXGO6Qm0k1EdaJIjHd0+zeaJhfgHk+pYPOc0974x+AIfoq2lMJvWHi8BgHTanyu6tDkskjSNS14oU+OY8Sg9ND9hmlRVu2ltMpDSaWH0ioP5dVeSqs9oVWfa+r91NT7z/Z/nMNuO1ujSUl0Br8mOElOcIa+xjltkky6QFQnCiGima4FmxjjnbbQzGDLCi7FEDBMvP5gk1W0dI5/mt2mk5PuJifd3ex+pRS1nkBoWG55jZeKGh+VtV7Ka3x4fMGatGEqyqq9lFV7L/gaDrtOsttJUoKDpHgnSW5H4z8nCXF2Et1OEuMdJMTZw740eyyTRCFEhCl1tolUtwU/HTvtOonxDgzTwjTBUBamGWy2CiYVC8sijMtEhI+maSS7nSS7nQzqnXze416/QVWdn6paHxW1PqrrfFTV+amu91Fd56e2wR9aAylgWJTXBJPNZ4lz2kiIc5AQb8ftsuOOc+COa/o++DXeZac+YOH3Boh32XA6bGFZxj3aSKIQoptSqml9InBy9tPwuYvRmVYwgQRMK1T7sBTNFpQ7r8O+mzdtxTnt5KTbz6uJNLEsRZ03QE1j0qhu8FNbH6DWE6CuITjzvM4T/Gee0zTt9QdraeU1bY9F04IJJs5pb/x69nuX00acI/jV6bDhOuef06E3fg1+77QH9yqP1maysCaKNWvW8OSTT2IYBl/+8pe5/fbbmz2+b98+li5dSn19PVdddRU/+clPsNsldwnRmnOv+6FE4jg7emjPsXI2bD9FZa2XzJR4plxxCbn90nhxw0F2HK4Irq6qawwbkELugHQ27zxNebWPpHgH40bmMKRvKvuPV/LOjtNU1vpIS3IxacwlnCqtY9POM/gME5fdxsTROUwd2++CcbZ0DuC8+3L7p7X5+Nz+aWz86GQwjoCJyxGMY9rV58ehlGLXkXI27TxNVZ2fBJedQZekkOR2crK0jhMltXj9waHMTocNw7Tw+sxmSVQp8PhMPL6LX9NJ08Bpt+G0B7ecddptOBr37nDYg/uWO+yNe5jbm38f3Nc8uD2trfF7m00P7TwYvL/xq372q82mhRZDvKjYw7UVanFxMbfeeiv//Oc/cTqd3HLLLfzv//4vl112Weg5c+bM4aGHHuLyyy/ngQceYOTIkdx2221tfo3y8rpu35ndUbG8HSNI+cJl5+Eylr9xAJstOOLKb1iYpkVGkovCE9XNnmu3acQ5dLLSE3A59NBCeFcOzWLn4XLQg8O2fQGTmjofHn9wRqGmKQwTDEvxuTG9mXJF3/Pi2H+8ktWbjzZe1HQCZrDTHk0j3mUL3WeaFvOuG3ResmjpeNO0GJCdyI4jFcGlybXg8uMK+PyVfc5LWhc6x9ihWWw/UHre/fOuG8SQfqnEJ7g4faYGr8/g4MkqNu86E5pgaZjBkWp9sxNxOWx4/Sb+QHAItC9g4Wu83d2uShpgs2nkZCTw5P2fb/fxYfv4vmXLFsaNG0dqaioA+fn5FBQU8M1vfhOAU6dO4fV6ufzyywFYsGABjz/+eLsSRazPApXyRbdIlO+9vcVkp7ubzU/wGyallR6y0+JbnDuUes4igH7DZOvuYhITHMQ57Og6uOMd+AMW8XHBIb/BfR6Cn7bPVHrJSotv3Hsi2CyEUpwsrSe3fxp2u97YSQ8VNV5QkJIcfD2lgnuEHzpZxejLMkMBKRT7T1TRJzsRu24LReo3TIqrvORkNK49FVqiXLH/RDX51w5o9rPYc7SixZ/F/hPVLd6/52gFIwdnkOx2YjU2e723t5i+vRLPe25SnIMvfX5Ii78DRXDElz9ghoZBG4aF3zAJGMFmQsMwCZgKwwgOlQ6YFkbTbdPEMBSGFbxtquDxpqkIWMF+KrPxsfYOdEhyO9r1/CZhSxQlJSVkZWWFbmdnZ7Nz584LPp6VlUVxcXG7XiMtLeHiA+3GOrIJejSR8nW+B+4e1+Wv2ZLv3Db2oo5fOjjrs5/0GX50acfPkZORcNHniCVhGy9mWVazjpuz2/G17XEhhBDdQ9gSRU5ODqWlpaHbpaWlZGdnX/DxsrKyZo8LIYToHsKWKCZMmMDWrVupqKjA4/Gwfv16Jk+eHHq8T58+uFwutm/fDsCqVauaPS6EEKJ7CNuoJwgOj33qqacIBAIsXLiQxYsXs3jxYpYsWcKoUaMoLCzkwQcfpK6ujhEjRvDwww/jdHbfbRmFEKInCmuiEEIIEf1k8RMhhBCtkkQhhBCiVZIohBBCtEoShRBCiFZFTaL47W9/y6xZs5g9ezbPPPMMEFwmZO7cuUyfPp3HHnsswhFevEceeYT//u//BoILJi5YsID8/HyWLl2KYXThbned7I477mD27NnMnz+f+fPns2PHDtasWcOsWbOYPn06y5cvj3SIF2XDhg0sWLCAmTNn8tBDDwGx87f597//PfR7mz9/PmPHjuWnP/1pzJQPgkPzZ8+ezezZs3nkkUeA2Hn/LVu2jPz8fObOncuTTz4JdLBsKgq899576pZbblGBQEB5PB51/fXXq3379qkpU6ao48ePq0AgoO6++261cePGSIfaYVu2bFHXXnutuv/++5VSSs2ePVt99NFHSimlvv/976vly5dHMrwOsyxLTZw4UQUCgdB9Z86cUddff72qrKxU9fX1au7cuergwYMRjLLjjh8/riZOnKiKioqU3+9Xt956q9q4cWNM/W02OXDggJo2bZo6ffp0zJSvoaFBXX311aq8vFwFAgG1cOFCtXnz5ph4/23evFnNmTNH1dbWKsMw1Ne+9jW1bt26DpUtKmoU11xzDc899xx2u53y8nJM06SmpoYBAwbQr18/7HY7c+fOpaCgINKhdkhVVRWPPfYY9957L9DygonRWrYjR44AcPfddzNv3jxeeOGFZgtGut3u0IKR0eiNN95g1qxZ5OTk4HA4eOyxx4iPj4+Zv81z/fjHP+a+++7jxIkTMVM+0zSxLAuPx4NhGBiGgd1uj4n33969e5k4cSKJiYnYbDYmTZrE888/36GyRUWiAHA4HDz++OPMnj2b8ePHt7joYHsXFewufvjDH3LfffeRnBzc7aszFkzsLmpqahg/fjxPPPEEzz77LH/72984ffp0zPzuPvnkE0zT5N5772X+/PmsWLEipv42m2zZsgWv18vMmTNjqnyJiYl8+9vfZubMmUyZMoU+ffrgcDhi4v03YsQINm3aRFVVFT6fjw0bNmC32ztUtqhJFABLlixh69atFBUVcezYsZhYVPDvf/87vXv3Zvz48aH7YmnBxCuuuIJf/epXJCUlkZ6ezsKFC3n88cdjpnymabJ161Z+8YtfsHLlSnbu3MmJEydipnxN/va3v3HXXXcBsfX3WVhYyD/+8Q/+9a9/8c4776DrOps3b46J8o0fP54FCxZwxx138JWvfIWxY8diGEaHyhYV28kdPnwYv9/PsGHDiI+PZ/r06RQUFGCznV0j/tOLDkaL1157jdLSUubPn091dTUNDQ1omhYzCyZ+8MEHBAKBUCJUStGnT59WF4yMJpmZmYwfP5709HQAbrjhhpj522zi9/vZtm0bv/zlL4HPXvAzmmzatInx48eTkZEBBJti/vSnP8XE+6+uro7p06eHEvzTTz9N3759+eCDD0LPaWvZoqJGcfLkSR588EH8fj9+v5+33nqLW265haNHj4aq/mvXro3KRQWfeeYZ1q5dy6pVq1iyZAlTp07l4YcfjpkFE2tra/nVr36Fz+ejrq6Ol19+mV//+tetLhgZTa6//no2bdpETU0NpmnyzjvvMGPGjJj422yyf/9+Bg4ciNsd3MxnzJgxMVO+vLw8tmzZQkNDA0opNmzYwDXXXBMT77+TJ0/yjW98A8MwqK2t5aWXXmLhwoUdKltU1CimTJnCzp07ufHGG7HZbEyfPp3Zs2eTnp7Ot771LXw+H1OmTGHGjBmRDrXTPProo80WTFy0aFGkQ+qQ66+/nh07dnDjjTdiWRa33XYbY8eO5b777mPRokWhBSNHjx4d6VA7ZMyYMXzlK1/htttuIxAIcN1113HrrbcyePDgmPnbPHHiBDk5OaHbLpeLX/7ylzFRvokTJ7J3714WLFiAw+Fg1KhRfPWrX2XatGlR//7Ly8tj+vTpzJs3D9M0ufPOOxk7dmyHri2yKKAQQohWRUXTkxBCiMiRRCGEEKJVkiiEEEK0ShKFEEKIVkmiEEII0aqoGB4rRFs89NBDbNu2DQhO0uzTpw9xcXEArFy5MvR9d6SU4q677uLxxx8PLeUiRHchw2NFTJo6dSq//e1vGTVqVKRDaRPDMBgxYgTbtm2TRCG6HalRiB7h4MGD/PznPw/NoL7zzju56aab2LJlC7/73e/Iysrik08+we1285WvfIXnn3+eY8eOMXPmTO6//362bNnC448/TnZ2NkePHsXtdvPwww8zePBg/H4/v/rVr9i+fTumaTJixAiWLl1KYmIikydPZuzYsRQWFvLd734Xy7J4+umn8fv9VFRUcPPNN/Otb32L73//+wDcfvvtPP3003zhC1/gqaeeYtiwYQBMnjyZp556CrfbzV133UX//v0pKipixYoVHD16lP/5n//B6/Wi6zpLlixhypQpkfxxi1jT+augCxF5119/vdq5c6dSSim/369mzpyp9u3bp5RSqrq6WuXn56udO3eqzZs3q+HDh4ceu/POO9Wtt96q/H6/KisrU8OGDVNlZWVq8+bNKi8vT23fvl0ppdTzzz+vvvCFLyillPrNb36jfv3rXyvLspRSSj3yyCPqZz/7mVJKqUmTJqk//OEPSimlTNNUt912mzp+/LhSSqnTp0+rvLw8VVVVpQKBgBo6dKiqrq4OHbd3795QeZpuHzt2TA0dOlR9+OGHSimlKioq1PTp09WpU6eUUkoVFRWpSZMmqaKiojD9ZEVPJDUKEfMOHz7MiRMnuP/++0P3+WMas2cAAAKLSURBVP1+9u3bR9++fenfvz95eXkA9OvXj8zMTBwOBxkZGbjdbqqqqoDgss1XXnklAF/4whd46KGHqK2tZePGjTQ0NPDOO+8AEAgEmi20NnbsWAB0Xeepp55i48aNrFq1ikOHDqGUwuv1kpCQ0ObyOBwOxowZA8CHH35IaWkpX//610OP67rOgQMHmi27IcTFkEQhYp5lWaSmprJq1arQfaWlpSQnJ7N9+3acTmez59vtLb8tzr3fsiwgeFE2TZMf/vCHXHfddUBw1c5AIBB6blMSqKur46abbiI/P5+xY8dy880388Ybb6Ba6CbUNK3Z/eeeLy4uDl3XQ3EMHTqUv/3tb6HHi4uLQ6vZCtEZZHisiHmXXXYZuq7z6quvAsEdBOfMmUNhYWG7zrN7924OHjwIBEdRXX311SQkJDBx4kSef/55AoEApmnywAMP8Jvf/Oa8448ePYrH4+Hb3/42119/PVu3bsUwDEzTxGazoWlaaP/i9PR0du/eDQQ3DaqoqGgxpiuuuILDhw+HVgPds2cP+fn5lJeXt6tsQrRGahQi5jmdTp588kl+8Ytf8Ic//AHDMPjP//xPxowZw5YtW9p8nuzsbB599FFOnTpFVlYWjzzyCADf+ta3eOSRR7jxxhtDndnf+973zjt++PDhTJw4kZkzZ+JwOMjLy2Pw4MEcP36cPn36MH36dG699VZ+//vf893vfpef/OQnLF++nFGjRoU6tT8tMzOTxx9/nIcffhi/349SikcffVSanUSnkuGxQrTBli1beOSRR5o1XwnRU0jTkxBCiFZJjUIIIUSrpEYhhBCiVZIohBBCtEoShRBCiFZJohBCCNEqSRRCCCFaJYlCCCFEq/5/LsThYimeHVsAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "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.7.1" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false + }, + "varInspector": { + "cols": { + "lenName": 16, + "lenType": 16, + "lenVar": 40 + }, + "kernels_config": { + "python": { + "delete_cmd_postfix": "", + "delete_cmd_prefix": "del ", + "library": "var_list.py", + "varRefreshCmd": "print(var_dic_list())" + }, + "r": { + "delete_cmd_postfix": ") ", + "delete_cmd_prefix": "rm(", + "library": "var_list.r", + "varRefreshCmd": "cat(var_dic_list()) " + } + }, + "types_to_exclude": [ + "module", + "function", + "builtin_function_or_method", + "instance", + "_Feature" + ], + "window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}