{ "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": 2, "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": 3, "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": 3, "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": 4, "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": [ "
" ] }, "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": 5, "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", "
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: Mon, 12 Nov 2018 Deviance: 3.0144
Time: 21:08:43 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]
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
" ], "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: Mon, 12 Nov 2018 Deviance: 3.0144\n", "Time: 21:08:43 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": 5, "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": 6, "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", "
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: Mon, 12 Nov 2018 Deviance: 18.086
Time: 21:08:47 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]
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" ], "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: Mon, 12 Nov 2018 Deviance: 18.086\n", "Time: 21:08:47 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": 6, "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": 7, "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" ] } ], "source": [ "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())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "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": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\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": 9, "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" ] } ], "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())" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "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": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\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": 12, "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd81fXZ+P/X53NWcpKQHYJsVAhbxQHKqFYJG0XaOm6po7S2d0vr73v39q7addfW2nrftva2VmurVaHF2irgAFRKlVFFVIYQNrJCyF5nfcb798dJDkRCSEJOzsj1fDwwOft6m+Rc572ut6aUUgghhBBnoMc6ACGEEPFNEoUQQog2SaIQQgjRJkkUQggh2iSJQgghRJskUQghhGiTJAohhBBtkkQhhBCiTVFPFA0NDcyaNYsjR46cdtvOnTuZN28excXF3H///ZimGe1whBBCdFBUE8WWLVu4+eabOXjwYKu3f/e73+UHP/gBq1atQinFiy++GM1whBBCdEJUE8WLL77ID3/4QwoKCk677ejRowQCAS666CIA5s2bx8qVK6MZjhBCiE5wRvPJf/rTn57xthMnTpCfnx+5nJ+fT1lZWTTDEUII0Qkxm8y2bRtN0yKXlVItLgshhIgPUe1RtKWwsJDy8vLI5YqKilaHqNqy/3AVhmEBoNGUZD6Ta7Sm/2ha8/c6OqBpWvi6Vr7qTY/RI4/R0JtuB+iOeru5uelUVjZE/4ViRNqXuJK5bZDc7dN1jezstA4/LmaJom/fvng8HjZv3sy4ceNYtmwZkydP7tBzGIZFyLA7+MpWu+6lRf7TlEQATQdd03HqGrpDw6Fr6Ho4ieha8/ddl0hsO7krwEv7Elcytw2Sv30d1e2JYuHChSxatIjRo0fzyCOP8MADD9DQ0MDIkSNZsGBBd4dzRiryn/CwGAA2gEXoM/cN91bC3Q+HruFy6DicOs6mZOLUdbQuTCBCCNGdtEQ+uGjX/vJO9Ci6V/NwlkPXcLscuCIJRP/sKFkL+fkZlJfXd1uc3U3al7iSuW2Q3O3TdY3c3PQOPy5mQ089hVLhHoltKwwznNQ0DXRNw+104Hbr4eShyyZ5IUR8kkQRA0qBpRT+kIk/FE4cDodGituJ2xnudQghRLyQRBEHlALTVDSYBhpGeJgqxY1h2bgcGqct5RJCiG4kiSLOKMC0FY0Bg+q6AE6HjjfFSYrbcXIJsBBCdCNJFHFMKTBMm9qGEA26RmqKk1S3E4cuCUMI0X0kUSQIy1Y0+Ax8AZNUjxOvx4FDJsCFEN1AEkWCsW1Fo9/AHzBJTZGEIYSIPkkUCcpWTQkj2NzDkCEpIUR0yEfRBNfcw6isC+ALmiTs7kkhRNySRJEkbFtR1xiiqi5AyIzv3epCiMQiiSLJGKZNdX2AxoAB0r8QQnQBSRRJSCmo9xnUNhooSRZCiHMkiSKJ+YMmVXVBbCVDUUKIzpNEkeQM06aqNohpSc9CCNE5kih6ANNWVNfLJLcQonMkUfQQlq2oaQjG/fkdQoj4I4miB7FtRU1jUHoWQogOkUTRw9i2orYhiClnAgsh2kkSRQ9k2Yqa+oCshhJCtIskih7KtBTV9SHsxD0yXQjRTSRR9GCGaVPnC8mWPCFEmyRR9HCBoEWj30CTwrNCiDOQRCFoDBj4Q1aswxBCxClJFAKloK4hhGHJ5LYQ4nSSKAQQPgiptjEok9tCiNNIohARphk+00JShRDiVJIoRAuBkNV0loUQQoRJohCnafQbBA2Z3BZChEmiEKdRCuoaQ1hS5kMIgSQKcQaWraj3h5DjVIUQCZ0oXnlnP3WNoViHkbQCQUv2VwghEjtR7D5Sy6MvbuG9HWWyrDNK6n0Gli37K4ToyRI6UXg9ToKGxbJ1B3j61R1U1gViHVLSsW1FXaOsghKiJ0voRHHnzOFcfGEeAAdL63nspa1s2H5cehddLGhYBGQVlBA9VkInCm+Kky9cfQFfnjaMXmluDNPm1Q0H+cOrO6iul95FV2rwhVAysS1EjxTVRLFixQpmzJjB1KlTWbx48Wm3f/LJJ9x4443MmTOHr33ta9TV1XXqdYYNyObb88cwbmg+AAdK63nspW18uLscJb2LLmFaCl9QehVC9ERRSxRlZWU8+uijLFmyhFdeeYWlS5eyd+/eFvf56U9/yqJFi1i+fDmDBw/mD3/4Q6dfL9Xj5MbPnc+CacNIT3URNCxeWruPJW/twSc7jbtEo9+QvRVC9EBRSxQbNmxg/PjxZGVl4fV6KS4uZuXKlS3uY9s2jY2NAPj9flJSUs75dYsGZLNo/hhGDMoG4JMDVfzmb9vYf6xzvRVxkm0rfEEz1mEIIbqZpqI0NvPkk0/i8/m45557APjrX//K1q1b+clPfhK5z8cff8ydd96J1+slNTWVF198kezs7Ha/xvHKxjN+wlVKsXFbKUvf3E3QsNA0mHHlYGZcNQiHntBTMzGlAbm9UvB4nLEORQjRTaL2127bNtopx6YppVpcDgQC3H///Tz77LOMGTOGZ555hnvvvZennnqq3a9RW+sjZJx5jX9Rv0z+/YZR/OXtPRyr9PHa+gNs31fBTZ+/kMw0d+ca1k1yctKoqmqMdRit8jUGyUp3cy4fMfLzMygvr++6oOJMMrcvmdsGyd0+XdfIzU3v+OOiEAsAhYWFlJeXRy6Xl5dTUFAQubx79248Hg9jxowB4Etf+hLvv/9+l8eRl5XK3dePYuLoPgB8erye//vbVvYcqeny1+opAiETw5S5CiF6iqgliiuvvJKNGzdSVVWF3+9n9erVTJ48OXL7wIEDOX78OPv37wfg7bffZvTo0VGJxenQmTFhIAuKh5HqcdAYMHn29RLe+uAwtkzOdphS4eNT5ZxtIXqGqA099e7dm3vuuYcFCxZgGAbz589nzJgxLFy4kEWLFjF69GgeeughvvOd76CUIjc3l5/97GfRCgeAooHZfHPeGP781m6OlDey5sOjHD7RwJeuuRBvioy5d0S4V+HC6ZBsIUSyi9pkdnfYtb+8zTmKMzEtmzf+dYiNnxwHIDvDw63XDeW8vLSuDrHT4nmOolmqx9npuYpkHgeG5G5fMrcNkrt9cTdHEc+cDp3ZVw3iC1efj8uhU10f5HfLtvPxnopYh5ZQZK5CiJ6hRyaKZhdfmM/d148kO8ODaSle/MdeVr73qcxbtJPMVQjRM/ToRAHQJzeNf79hNBf0zQTgnS2lPL9qF4GQbCxrj0DIxLCkDLkQyazHJwoIFxf88vQirhpVCMCuwzU88cp2KVveDkqBLyBJVYhkJomiiUPXmHnlIG6cMgSHrlFeE+CJl7dz8LiU/jibQNCSw42ESGKSKD5j3LAC7po1HG+KE1/Q5A+v7uSjPeVnf2APZiupLCtEMpNE0YpBhb34xvWjyM9KwbIVf/3HPt7efERKlrfBHzTlwCghkpQkijPI6ZXC3XNHRSa53958hL//c78MsZyBbSv8UllWiKQkiaINqR4nX54+jEuaDkTavLucP70hK6LOpDEgvQohkpEkirNw6Do3ThnCNZf0BWDv0Vp+v2IH9b5QjCOLP7at5GxtIZKQJIp20DSNay/tz41ThqBrUFrp43fLPqGyVpbPfpYvYMjJ2kIkGUkUHTBuWAH/VjysRdmPo+UNsQ4rrpimwpBehRBJRRJFBxUNyOauWcNJ9ThpDJj8/tUd7DtWG+uw4kpjwJSyHkIkEUkUnTCgdwZfmzOSzDQ3IcPmT2+UsPNgVazDihsh05JigUIkEUkUnVSQncrX5o4kLzMF01IsfnO3bMxrohT4gkaswxBCdBFJFOcgK93DV+eM5LxcL7aCv/5jH/9qOuOipwuELCypwitEUpBEcY7SU118ZfYIBhVmALB8/UHe3XIsxlHFniyVFSJ5SKLoAiluJ7fPKIrs4n7jvUNS8gPw+41OnX4nhIgvkii6iNvp4LbiYRQNyAbCJT9WvX+4RycL01YETelVCJHoJFF0IZdT59apFzJ6SA4A72w5xhvvHerRycIvS2WFSHiSKLqYQ9f50jUXcvGFeQCs21rKaxs/7bHJImRahEwppChEIpNEEQW6rnHjlPMjxQQ3bD/OivUHe2SyUAqpKitEgpNEESW6rjFvyhAuLSoA4F87ynpsspClskIkNkkUUaRrGtdPGszlw08mi1c39LxhKFkqK0Rik0QRZbqmMWfiYC5r6lls/OR4j5yz8Adkp7YQiUoSRTfQNY25kwZz6bCTcxZv/KtnrYYyLYVfDnwSIiFJougmuqZx/eQhjGtKFuu2lbJ6U8/aZ9EoZ1UIkZAkUXQjXdO4YdKQyNLZf358jDUfHo1xVN3HNBVBmasQIuFIouhm4dVQ50c25b29+Qj//LjnJItGv8xVCJFoJFHEgEPX+OI1FzBiULjcx6r3D7Nhe8+oOmuYtvQqhEgwkihixKHr3PT5CxnaPwuAVzcc5IOSEzGOqns0BsweNTcjRKKTRBFDTofOrdcNZXCfXgC8/M5+tuytiHFU0WeYFgFZASVEwpBEEWMup86C4mH0L0hHET78aOen1bEOK6qUAp9fhp+ESBTtShTPP/88DQ0N0Y6lx/K4Hdw+vYg+uV5spfjzW7vZleTJImia2DL8JERCaFei2LVrF8XFxdx///1s27at3U++YsUKZsyYwdSpU1m8ePFpt+/fv5/bbruNOXPmcNddd1FbW9v+yJNMqsfJHTOGR87g/u3ftnCkPHmTs1JIWQ8hEkS7EsWDDz7IqlWrGDVqFD/+8Y+58cYbeemllwgGg2d8TFlZGY8++ihLlizhlVdeYenSpezduzdyu1KKr3/96yxcuJDly5czfPhwnnrqqXNvUQJLT3Vx58zhZKa5CYYsnn29hLJqX6zDihp/wATZgidE3Gv3HEV6ejrTpk1j1qxZ1NTUsGTJEqZNm8aaNWtavf+GDRsYP348WVlZeL1eiouLWblyZeT2Tz75BK/Xy+TJkwG4++67ufXWW8+xOYkvK93DnTOHk+F14QuaPPPaTqrrA7EOKypMyyZkSqIQIt61K1Fs3LiR73znO0ybNo39+/fz+OOP8/e//50//elP/OAHP2j1MSdOnCA/Pz9yuaCggLKyssjlQ4cOkZeXx3333ccNN9zAD3/4Q7xe7zk2JznkZ6XyrS9ejMfloM5n8MfXS2hIwo1qSiGrn4RIAM723OnHP/4xt9xyCz/5yU/IyMiIXD9gwAC++MUvtvoY27bRTjkDUynV4rJpmrz//vu88MILjB49ml/96lf8/Oc/5+c//3m7g8/M9CbtOQc5wDe/MJbHXvyYytoAz6/ezf93yyWketr1I0sIOTlpaBpkZ3txOpJvAV5+fsbZ75SgkrltkPzt66h2vessX76clStXkpGRQXl5Oa+99hoLFixA13UWLVrU6mMKCwv54IMPIpfLy8spKCiIXM7Pz2fgwIGMHj0agFmzZp3xuc6kttZHyEjOYzZzctLITXdz0zUXsPjN3Rwuq+exv3zI7dOH43Im/ptqTk4aVVWNAJhBg1R38iRACL/RlJfXxzqMqEjmtkFyt0/XNXJz0zv+uPbc6Sc/+Qlr165teiGdzZs387Of/azNx1x55ZVs3LiRqqoq/H4/q1evjsxHAFx88cVUVVVRUlICwJo1axg5cmSHG5Dshg/KYd6U8wE4UFrP0jV7kq4X5Q9ayKS2EPGrXR/jPvroI1599VUAcnNz+fWvf83cuXPbfEzv3r255557WLBgAYZhMH/+fMaMGcPChQtZtGgRo0eP5vHHH+eBBx7A7/dTWFjIL37xi3NvURK6ZGg+voDJ6//6lB0Hq1m27gA3TBrcYigvkRmmRchUuJ3J0R4hkk27EoVhGIRCIdxuNxCeX2iP2bNnM3v27BbX/f73v498P3bsWF566aX2xtqjTRzThwa/wTtbjvFByQnSU11Mvax/rMPqEs2T2m6nO9ahCCFa0a5E8bnPfY677rqLuXPnomkar776KlOmTIl2bOIzii/vT6PfYPPuctZ+dJT0VBdXjiqMdVhdIhCySE9V6EnSSxIimbQrUfznf/4nixcv5u2338bpdHLddddx0003RTs28Rla0yl5jQGTkkPVvLbhIOmpLsacnxvr0M6ZbSsCIQtvEq3qEiJZaCqB6z3v2l+e1KuemlcFfVbItPjjazs5VNaAQ9e4fXoR5/fN7OYIz01r7XM6NXJ7pZIMfYpkXjmTzG2D5G5fVFc9vfXWW1xzzTWMGzeOSy65JPJPxIbb6WBBcREF2alYtuKF1bs5VtF6UkkkpqkwpP6TEHGnXf38X/7yl/zXf/0XI0aMSJqVNonOm+Lk9ulFPLnsE2obQ/zpjRK+NnckOb1SYh3aOfEFTTxuB4nbzxUi+bSrR9GrVy+mTp1Kv3796Nu3b+SfiK2sdA+3Ty8ixe2g3m/wzBslNAYSu9RH0LAwrOQcThQiUbUrUYwdO5Z//vOf0Y5FdELvHC8Lpg3D6dCorA3w3MpdhBJ4+Eap5g14Qoh40a6hp3/+85+88MILuFwuXC5XpG7Thx9+GO34RDsMKuzFF6+5kD+/uZvDJxr481t7+LfiYTj0xBwm9AdN0lKcslRWiDjRrkTx7LPPRjkMca5GDc5h9sRBLF93kF2Ha1j27n5umDwkIeeUbFvhD1mkyVJZIeJCu4ae+vbty7Zt23jxxRfJycnho48+kjmKODR+RCGfuzj8c/lgVzlvbz4S44g6z5+EZdWFSFTtShRPPfUUf/7zn1m5ciWBQID/+7//4/HHH492bKITrru0H5cMzQNgzYdH2bSz7CyPiE+mreSoVCHiRLsSxWuvvcbvf/97UlNTyc7O5sUXX4wUCRTxRdM0bpg8hAv7hTfgvbLuACWfVsc4qs5plF6FEHGhXYnC6XRGCgJCeLms0ynjx/HKoevcct1Q+ualoRT8+a09HD6ReDtNTdNO2p33QiSSdiWKPn36sHbtWjRNIxQK8cQTT8gcRZzzuBwsmDaMnAwPhmXzp5W7qKj1xzqsDlGAL2CQgPPxQiSVdiWK73//+zzzzDPs2rWLiy66iHfeeYfvf//70Y5NnKMMr5vbpxfh9TjxBUyefb2Eel8o1mF1SNC0MEzZpi1ELHWoKKDf78eyLNLTO15UKhp6alHAjjpUVs8fXt2JYdn0zU/jK7NG4HE5uuS5O6sj7fOmOOnlTayzKpK5sFwytw2Su32dLQrYromGZ555ptXr77jjjg6/oOh+A3pncNO1F/LC6l0cLW/kz2/t4bbioTj0xDh7OxCySEtRCbuBUIhE165EsXv37sj3oVCITZs2MWHChKgFJbre8IHZzLlqMMvWHWD34RpeefcA8xJkQ154A55Jeoor1qEI0SO1K1E89NBDLS6XlZVx//33RyUgET1XjOhNXWOIf3x0lM27yslMc3PtpYlxnKo/EC7roSXFaRVCJJZOjT307t2bo0ePdnUsohtce2k/LhmaDzRtyCs5EeOI2sdqKushhOh+HZ6jUEqxfft2cnMT//jNnii8IW8w9b4Qe47Usuzd/WSkuigamB3r0M7KFzBIdTulTyFEN2tXj2L37t2Rf3v27KFPnz488sgj0Y5NRMmpG/LsBNqQZ5oqoUuoC5Go5MzsONWVy2PPpN4X4nfLPqG6Pog3xcndc0eSl5ka1dds1tn2uV0OcjLcEOf9imReYpnMbYPkbl9Ul8fedtttba6Oee655zr8wiL2Mrxu7phRxO+WfYIvYPLM6yXcPXckGXG8Z8EwLUKmwu2M70QhRDJp19DTqFGjSElJYcGCBdx1113k5eWRlZXFrbfeyq233hrtGEUU5WWm8uVpw3A5dKrrg/xp5S6CcTxpHD4Bz5SyHkJ0o3b1KD788EOWLFmCwxHezTtp0iS++MUvUlxcHNXgRPfoX5DBzU0b8o5VNLLkrd3cVjwMpyM+N+QFQxamZSfMhkEhEl27/tKqqqoIBoORy42NjQQCgagFJbpf0cBs5k4aAsCeI7W8/M5+4nX6ylYKn5yrLUS3aVePYtasWXzpS1/iuuuuQynFG2+8wYIFC6Idm+hmlxUVUNcY4u3NR/hoTwUZXhfTrhgY67BaJedqC9F92pUovv3tbzNixAj+9a9/4fF4+O///m8uv/zyaMcmYuCaS/pS7wvx/s4TvLOllAyvm6tG94l1WKeRc7WF6D7tHuTt3bs3F154Id/5zndwuaTmTrLSNI05Vw1mxKDwBrzXNn7Klr0VMY6qdX6/gSI+h8eESCbtShR/+9vf+N73vsfTTz9NfX093/jGN3jxxRejHZuIEV3X+NI1FzKwMAOAl9buY++R2hhHdTpTynoI0S3alSheeOEFli5dSnp6Orm5ufz973/nT3/6U7RjEzHkcuosKB5GQXYqlq144c1dHC1viHVYp/EFDOlTCBFl7UoUuq63OKyoT58+kaWyInmlepzcMb2IzDQ3IcPm2TdK4u44VdNUBKWshxBR1a5EkZWVxc6dOyO7s5cvX05mZmZUAxPxITPdwx0zhpPqcdLYtHu7Ls6OU230G7EOQYik1q5Ecd999/Hd736Xffv2MXHiRH7961/zwAMPRDs2EScKslO5ffowXM7w7u1nXy/BHzRjHVaEadpJW/NLiHjQrkQRCARYtmwZL7/8Mn/84x9ZuXIlw4YNO+vjVqxYwYwZM5g6dSqLFy8+4/3Wrl3LNddc0/6oRbfrX5DBrdcNRdc0jlf5eG7VLkJmfAz5KMJzFbKlQojoaFei+I//+A8cDgfnn38+Q4cObdfy2LKyMh599FGWLFnCK6+8wtKlS9m7d+9p96uoqODhhx/ueOSi2w3tn8X8q88H4NPj9fzlrT1Ydnx8kg+aFoYZH7EIkWzalSiGDRvGihUrOHbsGDU1NZF/bdmwYQPjx48nKysLr9dLcXExK1euPO1+DzzwAN/85jc7F73odhddkMesKwcBUHKohr+t3Y8dB6U+lEKWygoRJe3a1vr222+f9iavaRo7d+4842NOnDhBfn5+5HJBQQFbt25tcZ/nnnuOESNGMHbs2I7EHJGZ6cWyY/8mFS05OWmxDqFVsyafj9I0Xlt/gI/3VpCdmcIXrx3aZin61nR1+zQNsrO9cVPMMD8/I9YhRE0ytw2Sv30d1a5EsW3btg4/sW3bLd44lFItLu/evZvVq1fz7LPPcvz48Q4/P0BtrS9pJzG74+Cic3HliAIqqn28t6OMf2w+gg58fly/dj8+Wu0zggZed+zLeiTz4TfJ3DZI7vZ19uCiNj96ff/73498X1VV1aEnLiwspLy8PHK5vLycgoKCyOWVK1dSXl7OjTfeyFe/+lVOnDjBLbfc0qHXELGjaRqzrxrE2AvCZ6e/vfkI67eVxjgq2YAnRDS0mSi2b98e+f6uu+7q0BNfeeWVbNy4kaqqKvx+P6tXr2by5MmR2xctWsSqVatYtmwZTz31FAUFBSxZsqSD4YtY0jWN+Z87n6IBWUC4LtSHu8vP8qjoMk2FIRvwhOhSbSaKU88j6OjZBL179+aee+5hwYIFXH/99cyaNYsxY8awcOHCTg1lifjk0HVuvnYog/uEx3T/9s99bNtfGdOYfAE5AU+IrtTuwdyOTlQCzJ49m9mzZ7e47ve///1p9+vXrx9r1qzp8POL+OBy6txWPIw/vraTI+WNvLhmL26nzrAB2TGJJ2TamKaNI04mtYVIdG3+Jdm2TW1tLTU1NViWFfm+PctjRc+S4nZy+/QiejcVEVz85m72H4tNxVlbKQKyp0KILqOpNsaUioqK0DSt1WGnsy2P7Q679pfLqqc4U+8L8dSKHVTWBnC7dO6cMZwBvU9fahjt9rmcGrm9UoDYjEEl88qZZG4bJHf7Orvqqc2hp5KSkk4HJHqmDK+bu2YO56nln1DTEOLZN0q4a+Zw+uZ3/JfzXJiWwrAULodMVghxrmQQV3S5rHQPd80cQS+vi0DI4o+vl1Ba2b29I6UgIDu1hegSkihEVORmpnDnrBGkpbrwB03++NpOTlR371kWwaBJHFQXESLhSaIQUVOQlcpdM0+eZfGHV3dQXtN9ycK0lRQKFKILxL7WgUhqhTle7po5nKdf3UG93+DpV3ewcNaIc6rztOtQNe9uOUZ1fZDsDA+Txp53xqW4gZCJ2+Xu9GuJ2Nm6r4KV7x2iojZAXmYK064YwJjz82IdVo8kPQoRdeflpXHXzOGkuB3U+8LJ4kS1r1PPtetQNcvXH6DOb5DicVLnN1i+/gC7DlW3ev+AYcVFdVvRMVv3VbD4zd3UNIbwpjipaQyx+M3dbN1XEevQeiRJFKJb9M1P584Zw/G4HNT5DP53yYedOn/73S3HcDh03E4HmqbhdjpwOHTe3XKs1fvbtiIgJT0Szsr3DuFw6Hhc4Z+zxxX+Oa9871CsQ+uRJFGIbtOvIJ07ZxbhcTmoqQ/y9IqOz1lU1wdxfWbHtcsRPqL1THwBQya1E0xFbQC3s+XP2e3UqagNxCiink0ShehW/QsyuHNmESmecM/i6RU7ONGBZJGd4cGwWk5QG5ZNdobnjI8xTUUwTo5tFe2Tl5lC6DMLEUKmTV5mSowi6tkkUYhu178gg+/cdEl4zsIfThZlVe2bs5g09jwsyyZkWiilCJkWlmUzaex5bT6u0R+S8uMJZNoVA7Asm6AR/jkHjfDPedoVA2IdWo8kiULExKA+vZqWzjpo8Bv8/tUdHKs4+6a8YQOymXPVYHqluggETXqluphz1eCzFiA0zPCbjUgMY87P49brhpKV5sYXMMlKc3PrdUNl1VOMtFnrKd5JrafE1dy+YxWN/PH1nfgCJiluB3fMKKJ/QXSOoXQ59ab6T9GXzPWCkrltkNzti8oJd0JE23l5aSycPYKM1KZyH6+VcKC0LiqvZZi29CqE6ARJFCLmemd7WThnBJlpboKGxbOvl5xxX8S5agyYUXleIZKZJAoRF/IyU/nqnBHk9krBsGyeX7Wbrfu6/qS8kGHJvgohOkgShYgb2RkpfHXOCApzvNhKsfTtPWwqOdHlr9PgD3X5cwqRzCRRiLiS4XXzlVkj6F+QjgJefmc/az862uEz29timgpfUIaghGgvSRQi7nhTnNw5czgX9M0EYPWmw7y28dMurdnjstBEAAAgAElEQVTU6DekBpQQ7SSJQsQlj8vBgmnDGD0kF4AN24/z4pq9mFbXLIe2bOlVCNFekihE3HI6dL70+QsYP7I3AFv3VfKnlSUEQl3zBu8LmFi29CqEOBtJFCKu6ZrG7CsHcd2l/QHYd7SOp5bvoLbx3CekbelVCNEukihE3NM0jasv6cv8z52Prmkcr/Lxu1e2c7yd9aHa4g+Y2HZy7u4XoqtIohAJ45Kh+Xx5+jA8Lge1jSGeXPYJuw/XnNNz2krRIJvwhGiTJAqRUC7sl8VX54ygV9Mu7udWlvDejrJzes5A0MKUXoUQZySJQiScPrlpfP36UZyX68VWsGzdAV7beBC7kxPTtlL4pFchxBlJohAJKTPNzcI5Ixk+MFxefP224zy3alenV0T5gyamJSughGiNJAqRsDwuB7deN5SJY/oAsPtwDU+8sr1TZ3ErBfW+kByZKkQrJFGIhKbrGjPGD+TGKUNw6BrlNQF++/J29hzp+CR30LBoCBhoWhQCFSKBSaIQSWHcsAIWzh5BetO5Fs++XtKpGlG+gCF7K4T4DEkUImkM6J3BN24YRb/8NBThGlGL39xNMNT+suJKQX2jgdFFpUKESAaSKERSyUr3sHD2SC4dlg/AjoPVPP7yNso6sDnPVorahiCmlPcQApBEIZKQy6lzw+QhXD9pMA5do6I2wG9f2c5Hu8vb/RympaiqC8jRqUIQ5USxYsUKZsyYwdSpU1m8ePFpt7/11lvMnTuXOXPm8I1vfIPa2tpohiN6EE3TuHx4b746ZyRZ6W4M0+ava/fx8jv7Mcz2DSvZtqKmIUhjwACkdyF6rqglirKyMh599FGWLFnCK6+8wtKlS9m7d2/k9oaGBn70ox/x1FNPsXz5coYNG8ZvfvObaIUjeqj+Bel8c95ohvXPAmBTyQmeeGU7ZdXtG4oKL5s1qPMZ0QxTiLgWtUSxYcMGxo8fT1ZWFl6vl+LiYlauXBm53TAMfvjDH9K7d7iE9LBhwygtLY1WOKIH86a4uG3aMIov74+uwfEqH7/9+3Y2lZxo96ooX8CkzidHqIqeKWqJ4sSJE+Tn50cuFxQUUFZ2siZPdnY21113HQCBQICnnnqKa6+9NlrhiB5O1zSmXNSXhbObhqIsm5ff2c+f39qDL9C+3oIkC9FTOaP1xLZto52yc0kp1eJys/r6ev793/+doqIibrjhhg69RmamN6kPnsnJSYt1CFEVi/bl5KQxbEguL7xRwoe7TrD9QBWHyxv58szhjBic267n0N1OsjI8OPS2d+bl52d0RchxKZnbBsnfvo6KWqIoLCzkgw8+iFwuLy+noKCgxX1OnDjBXXfdxfjx47nvvvs6/Bq1tT5CRnKud8/JSaOqqjHWYURNrNt34+TBDOqdzqsbD1LbEOSxpR8zYVQhxZf3x+10nPXxlZUOeqW5z5gs8vMzKC+v7+Ko40Mytw2Su326rpGbm97xx0UhFgCuvPJKNm7cSFVVFX6/n9WrVzN58uTI7ZZlcffddzN9+nTuv//+VnsbQkSLpmlcWlTAt24cw4De4T+cjduP85u/bePT42d/kwgaFtX1AULtXEElRCKLWo+id+/e3HPPPSxYsADDMJg/fz5jxoxh4cKFLFq0iOPHj7Njxw4sy2LVqlUAjBo1ip/+9KfRCkmI0+T2SmHh7JG88/Ex1nx4hMraAE8t/4SJY/pw7aX9cTnP/FnKtBQ19UHSUp2kpbi6MWohupemOloMJ47s2l8uQ08JKh7bt/GTUla9dzjSS8jwurhyZCF7jtRQXR8kO8PDpLHnMWxA9mmPdTo1eqW68bgdKNX54Yut+ypY+d4hKmoD5GWmMO2KAYw5P++c29YVlq/bz+pNRwgYFikuB1Mv68eciUNiHVaXk6GnVh4XhViESDi7DlWzbmspGWku0lPDvYN6n8GqTYc5WunD7XJQ5zdYvv4Auw5Vn/Z401RUNwSprg92uk7U1n0VLH5zNzWNIbwpTmoaQyx+czdb91WcU9u6wvJ1+1m+4SBBw8Kph4felm84yPJ1+2MdmugGkiiEAN7dcgyHQ8fjctIrzU1+VirNs2bBkEV5TQDTVOi6xrtbjrX6HEpBIGRRVRegosaPaakOlSxf+d6hphgcaJqGx+XA4dBZ+d6hc2/gOVq96QgaGg5dQ9P08Fc0Vm86EuvQRDeI2hyFEImkuj5Iiufkn0Pz3ETz+7ytwuU83E4d4yzDnUo1TXbX+fG4HHg9Ltyus38mq6gN4E1p+SfpdupU1AY61pgoCITM01Z46RqdPlFQJBbpUQgBZGd4Thsycjo0nA6NguxUUtzhJbMh06bOZ7DyvU/PWr68uYdRXR+gsi5A4CwFBvMyU05bRRUybfIyUzrRoq6V4nby2S1LtgpfL5KfJAohgEljz8OybEKmhVKKkGnhdjlwu51YSpGd4aFXmovmD9XvbCnlf5d+zIe7y7HPsh5EAYZpU1MfpKLGT2PQbPUx064YgGXZBI1wDEHDwrJspl0xIAot7pipl/VDobBshVJ2+CuKqZf1i3Voohs4fvSjH/0o1kF0VmW1L2l3ZqemuvH7k7cQXby1Ly8zlbzMFMoqfTT4DLLS3EwfP5ARg7Ij1+X2SmH6+IH0zvFyuKwBf8hix8Fqdh2uIT8rlewMT+T5ztQ+W0HIsPCHLBTg0DX0pomM3jleemencuREA7WNIXIyPMybPCQuVj0NG5ANSvHp8QZCliLF5WDG+AFJueopLc2DL0lLtWiahtfr7vjjEnl57L5PKwk2/cEppUAlTzHoeFw+2pUSvX01DUHe+Nchtu2vjFw3YlA20y4fQF5Warvbp2saKR4HqR4nbqdOIvw1JvPyUUju9nV2eWxCJ4qqqgYsS2ErhW0rbEWke2zbCstUWEph23b4tqamJkKLE/2N9GySpX2fHq/n9X99yuETDUB4gnfcsALmXXMhymz/oUeaBi6nA2+KE49LRyN+KxUk8xspJHf7emSiqKxswD7L0FPz8sQWycRuTiBNScVSmMpG2cRN7yRZ3kjPJJnap5Ri2/4qVr1/iOr6IABOh874kb2ZPPa8yL6M9tJ1jVSPE4/LgbuNneGxksxvpJDc7etsokj6JQvNabB5DbgD4DM135qTiVKcsXdiW+HvrabeCU23xzqhiNjTNI0x5+cyYlA2m0pO8I8Pj9LgN1i3tZT3dpQxYWRvJo5pf8KwbUWj38AXMHA6NNwuJx6njtMZ3r+QuB/tRKJK+h5FVzrZOwl/imxOKHYrvRPbUufUO0mmT9ytSeb2hQyLj/dXsepfB/EHw8NPbqfOFSN6c9WYPvTqxGQihD8NOh06Hne4p+F06DEZoErmT9yQ3O2THkU3ONk7CX+K1B2t/5lq2sneSXiYq2UysSw73DOxw/dTKPmUmETcLgfTJgxizOBsNmw/zrqtpQRCFu9uLWXjJ8cZN6yASWP6kNOrY/sjbFsRsi1ChoWmgUPTcLsdpLgc0tsQUSWJIgqa/1j18F/zZ0e6gJbJxLYVpq0wLRvDsDEtO46nMkV7pbidXHNJP64cVcjG7WWs316KL2Dy3o4y3t9ZxqjBuUwa04d+BR3/hKcUmEphBkx8ARNd13A5dNwuHZfTgcupxfWEuEgskihi5NRkojs0nA7A5UBLDQ9tZWalokwTw7QxzHDyCPdAlMyJJJgUt5OrL+nLVaML2VRygne3llLXGGLb/kq27a9kUGEGV47uw/CB2Wc9Ne9MbFsRtC2ChoWGgaZreJwOPB4HLoeO0yG9DdF5kijijFLhoa3wihdHi9PWrKahK9O2MU0bw7KxLBVOHvImEPfcLgdXje7DFSN6s21fJeu2lVJa6ePg8XoOHq8nK93NFSN6c1lRAd5zON9CEV7Z5w+Z+EMmuqbhcIDb5cTl1HHpGg5H/K2mEvFLEkUCcejhlVtudHCfHL4yLRtLKUxLYRo2hm2d7H1IAok7TofOxUPzuejCPPYdrWP99lJ2H6qhpiHEqvcP8/bmI4wanMvlIwoY2DvjnE9/tJXCNsEwwzvFm+c3nE4dt8uBU9ciE+XNv1NCnEoSRQJr/oN2OnScgMcJeMJvBM3Les2myfPI8JWSBBIvNE3jgn6ZXNAvk8raAP/65Dgf7ConaFh8vLeCj/dWUJCdyqXDCrjowrwO78c4k8j8Rsgi0FTYUNPC8Th0cDkcOJzh4armDye6JkNXPZksj41T0Vii11oCMU3VogfSXftCknl5LHS+fSHDYuu+St7bWcbR8pOP1zWNooFZXDI0n6H9s3B209DRyQSi4XY5cDl1CvIzqKvxJW3vQ5bHnk56FD2IUqdPnkP4zcBWqqkcSviOVlOvo3lVlrJpuq5pUj3ypLK8tyu5XQ4uLSrg0qICjlY08kHJCbbsrSDQVIBwx8FqUj1ORg/J4aIL8xjQOyNSVDAa1Cl7hoymEugOl5OaGn94vqNpI6CuhXsduq6hJ2kC6cmkRxGn4vFTzanvR6cOYSkVrrcV7qWES3SfbZJdehTtZ5g2Ow5W8eHucvYerW3x/zQzzc3oIbmMPj+Xfvlp5zyf0R6ttU3TwtUP0MJDoa7I0FV4f0ciJZB4/NvrKtKjEFF36h9586bD8FL98JtT835jTXNFeiiWUphmeIiruUeiAIcWnpw/7X2juTRKK6/ZU7mcOmMvyGPsBXnU+UJs3VvJx3srOFbRSG1jiHXbSlm3rZTMNDcjBucwclAOgwoz0Du51LYzmjeOoohsCoSTvyeaBrpDw6XrOBzhVVc6Grp+cmjr1BwnP/f4Ij2KOJWMn2pOfSPIy8ugoqIhfICPArSTvY9ITwXChRq1k+VSLNMmZJ0yp9KGWP5md0ePqaLGz9b9lWzdV8mJan+L27weJ8MGZFE0MJsL+2V26Ul0Xdm2U3siGqDr4NB0dF1Da+qFaE2T6ScTS/iTsYbWok5bV0nGv71m0qMQce+zf8xKqVM6JKd8ojw1o3x2W/upq7qas8lnNQ1xNFcJVip8nd703LoGqKZ3JgWWsjGtpsl9K3FK0udlpXLNJf245pJ+nKj2s+NgFZ8cqOJoRSO+oMlHeyr4aE8FuqYxsDCDYf2zGDogi97Zqd0yRNUep/ZEACwbDM5cnr05bE0L7zvXdHA2JRbd0TzEFf7naOrFNM/hxPPPMt5JohAJJzIpf7Y3u9Zqp7QqvILo1JL0VqSoY/g/NrSo1RUy7LMegdqdCrJTKcjuy+cu7kt1fZCSQ9WUfFrN/mN1WLbiQGkdB0rrWPn+ITK8Li7om8n5Tf8y0zpXpDAWTu11AmCDeYbE0rxiS9dA13WcuoajOZnoGlrTh4VTR+g0LbzyzLJVi88rzT2fZF3pdTaSKIRocmpJeucZCj5C+E3DsBQNvlDT+dbdFGA7ZWd4mDCykAkjCwmGLPYdq2XXoRp2H66htjFEvc+I9DYAcjNTGNKnF0PO68Wgwgwy0z1neYXEEFloAWBZtHa4aWufNWzdQXVt4JQ7nfyi6zRN0OtoOpGVXpoGmgrfK3K9Fh4q+2xMiUgShRAdpBQ4dY3sDA/BkI0vaGCY4R5GvL0ReNwORgzKYcSgHJRSlNcE2Hu0lr1HatlfWkvIsKmsDVBZG2BTyQkAstLdDCrsxYDCdAYUZNA7x9vpGlTxrrWfl4KWvcVTvj05NNb28FjzvIuuEZl/0TS9qXdzcoisuddzelwte7SRJz71dU4Pr9VGnXrJ6dDJPWPkZyaJQohOUgrcLh23y4MCLDu8A74xYMQ6tFZpmtY0RJXKlaMKsWybo+WNHCitY/+xOj49Xk/ItKlpCEV2hkP4LI2++Wn0y0+nX0E6I3UdTam4meeIN6fOu9gtbmn/0bjR4urkiYmy6ilOJfPKC0ju9inAm+7hyLHaprmO+OtptMayFaWVjXx6vJ6DpfUcOlFPva/1pOf1ODkvL43z8tLok+ulT24aeZkp3bokN1qSeY+Py6lTdH5+hx8nPQohupgGpKe6ye3lwbJPrtBqLpli2hZWHBZtdOhauNeQn85Vo/uglKK2McShsnqOnGjkcHkDx8obMSwbX9AMD2EdrY083uXQKchJpTDbS+8cL71zUinISqVXmlt6HwlOEoUQURMushcZ3z+lZIpSEDJtDDNcmC9cPiWOsgbhoaqsdA9Z6R7GnJ8HhJNe0FaU7K/kaEUjx8obOV7lI2hYGFZ4KOvUGlUQLpmfn5VCflYqeZmp5DV9n9PL06KMvohfkiiE6GbN+cDt1HE7ddJTXeES8bZNKGQTNEwsO756G80cuka/vHS8Tp1LhoaHMGylqKkPUlrp43iVj7JqH2VVPiprA9gKgobFkfJGjpSfPpzTKy3c88rplUJORgo5vTxkZ4T/pae6pCcSJyRRCBFjSjWfNeLA43SQgQvDsgkZFv6QiW3F98mGuqaF3+h7pTBycE7ketOyqagNUF7j50S1n/IaP5W1Acpr/YSM8DRvXWOIusYQB0pPn69yOk72aDLT3eGvaW56pbkjX1PcDkkm3UAShRBxyOXQcTnCvQ3LDh9K1byqyjDjc47js5wOncIcL4U53hbXK6Wo9xuRZbmVdQGq6oJU1weorAviD5oAmJaiojZAxal7Gj7D5dTp5XWTkeYiI9VNhtfV9M9NWoqTdK+b9FQXaSnObivNnowkUQgRx5p3obudGqCT6j5ZwsS0VLjnEQrPDyg7fnsdp9I0jV5eN728bgb36XXa7YGQSU1DiJr6IFX1QWobgtQ0hKhtDFLbEKLeF6J5saNh2lTWhZPN2aS4HaSluEhLdeL1OPGmuPCmNH8f/prqcdJo2IQCBqkeB26XI6pl3BOFJAohEsypycPt1ElPcUYSh2nbmMbJwomJeF5IittJYY7ztJ5IM9tWNAQM6pqSRq0vRH2jQb3foMEX3nne4A//s05ZPh9oOtGvsq79sWhaOMGkuJ1NX09+73E7SHGFv7pdDjyn/HO79Kav4e/dTgdOh5aww2RRTRQrVqzgiSeewDRNvvzlL3Prrbe2uH3nzp3cf//9NDY2cumll/LjH/8Yp1NylxAd0SJxoLP9aCX/+PAo1fUB8jJTmXJxXy7sm8Xf/7mXrfurMCwb21IUDchkxOBc1m09RmVtgF5pbiaNPY9hA7LZdaiad7cco7o+SHaGh0ljz+NoeQPrth4naFp4nA4mjinkmnH9zxhXa88BnHbdsAHZ7X78sAHZrP3oSDgOw8LjCsdx3WWnx6GUYtv+StZtPUZNQ4g0j5PB52WS4XVzpLyBwyfqCYSsyOl9pmUTCFotemVKgT9o4Q+e+2Y5TQO304HbqeNqSh7Nhz+5nOGhxvA5Hievc556nSN8rrmj6XuHI3y9M3J901f95FfHKWeCnFPs0dpwV1ZWxs0338zf//533G43N910E//7v//LBRdcELnPrFmzePDBB7nooou47777GDVqFLfccku7X0M23CUuaV90bN1XweI3d+NwhFdUhUwby7LJzfBQcrj25GFCTW8qbqdGXlYqHpcDywbTshg1KIeP9lagCCegoGFR1xDAHwx3UTQUhhXeWPj5S/q2mix2Hapm+foDTW9qOoZl4w+YoGmkehyR6yzLZs5Vg09LFq093rJsBhaks2V/VbjuktZ0gNYZ4jjTc4wbms/m3eWnXT/nqsFc2D+L1DQPx47XEQia7DlSw/ptxyOlNkzLxrYV/QrS8bgcBELhszeChkXQsAk2XY63dyUNcDg0CnPTeOLez3f48VH7+L5hwwbGjx9PVlYWAMXFxaxcuZJvfvObABw9epRAIMBFF10EwLx583jsscc6lCiSYRdoW6R9iS0W7XtvRxkFOd4W+xNCpkV5tZ+C7FROjaj5zSzV4wLCBe/QFB/traRXmie8ogiFruuccDvR9PCSXl3TI8fnVtaFyMrwhM8N4WTp94NldQzpm4XDcbLEd1Vt+MyM7IyU8GsrCFkWuw7XMPr83BaB7fy0mj55aTgdjkikIcPieLWf3tmpLf7f2rZi1+Faiq8Y2OL/xScHqlr9f7HrcG2r139yoIpRQ3Lp5XVjNw17vbejjH6900+7b0aKiy99/sJWfwaK8IqvkGFhmDYh08Y0bUKmhWGG55VM08Kwwod6GaYdvq75smU1bcwMX7ZU+PGWpTCaqhdbTbd1dO9NhtfVofs3i1qiOHHiBPn5J7eKFxQUsHXr1jPenp+fT1lZWYdeIzs77dwDjWOdOWAkkUj7ut59d47v9tdszX/82+Xn9Pjvd6LMxGf98ByeozA37ZyfI5lEbb2YbdstJm7UZ4qIne12IYQQ8SFqiaKwsJDy8vLI5fLycgoKCs54e0VFRYvbhRBCxIeoJYorr7ySjRs3UlVVhd/vZ/Xq1UyePDlye9++ffF4PGzevBmAZcuWtbhdCCFEfIhqmfEVK1bw5JNPYhgG8+fPZ+HChSxcuJBFixYxevRoSkpKeOCBB2hoaGDkyJE89NBDuN2JcyyjEEL0BAl9HoUQQojok+InQggh2iSJQgghRJskUQghhGiTJAohhBBtSphE8etf/5oZM2Ywc+ZMnnnmGSBcJmT27NlMnTqVRx99NMYRnruHH36Y//qv/wLCBRPnzZtHcXEx999/P6Zpxji6zrvtttuYOXMmc+fOZe7cuWzZsoUVK1YwY8YMpk6dyuLFi2Md4jlZs2YN8+bNY/r06Tz44INA8vxu/vWvf4383ObOncu4ceP47//+76RpH4SX5s+cOZOZM2fy8MMPA8nz9/fUU09RXFzM7NmzeeKJJ4BOtk0lgPfee0/ddNNNyjAM5ff71dVXX6127typpkyZog4dOqQMw1B33nmnWrt2baxD7bQNGzaoK664Qt17771KKaVmzpypPvroI6WUUt/73vfU4sWLYxlep9m2rSZOnKgMw4hcd/z4cXX11Ver6upq1djYqGbPnq327NkTwyg779ChQ2rixImqtLRUhUIhdfPNN6u1a9cm1e9ms927d6vrrrtOHTt2LGna5/P51GWXXaYqKyuVYRhq/vz5av369Unx97d+/Xo1a9YsVV9fr0zTVF/72tfUqlWrOtW2hOhRXH755Tz33HM4nU4qKyuxLIu6ujoGDhxI//79cTqdzJ49m5UrV8Y61E6pqanh0Ucf5e677wZaL5iYqG3bv38/AHfeeSdz5szhhRdeaFEw0uv1RgpGJqI333yTGTNmUFhYiMvl4tFHHyU1NTVpfjdP9aMf/Yh77rmHw4cPJ037LMvCtm38fj+maWKaJk6nMyn+/nbs2MHEiRNJT0/H4XAwadIknn/++U61LSESBYDL5eKxxx5j5syZTJgwodWigx0tKhgvfvCDH3DPPffQq1f4tK+uKJgYL+rq6pgwYQKPP/44zz77LH/5y184duxY0vzsPv30UyzL4u6772bu3LksWbIkqX43m23YsIFAIMD06dOTqn3p6el8+9vfZvr06UyZMoW+ffvicrmS4u9v5MiRrFu3jpqaGoLBIGvWrMHpdHaqbQmTKAAWLVrExo0bKS0t5eDBg0lRVPCvf/0rffr0YcKECZHrkqlg4sUXX8wvfvELMjIyyMnJYf78+Tz22GNJ0z7Lsti4cSM/+9nPWLp0KVu3buXw4cNJ075mf/nLX7jjjjuA5Pr9LCkp4W9/+xv/+Mc/ePfdd9F1nfXr1ydF+yZMmMC8efO47bbb+MpXvsK4ceMwTbNTbUuI4+T27dtHKBRi+PDhpKamMnXqVFauXInDcbJG/GeLDiaK119/nfLycubOnUttbS0+nw9N05KmYOIHH3yAYRiRRKiUom/fvm0WjEwkeXl5TJgwgZycHACuvfbapPndbBYKhdi0aRM///nPgbMX/Ewk69atY8KECeTmhs/DmDdvHn/4wx+S4u+voaGBqVOnRhL8008/Tb9+/fjggw8i92lv2xKiR3HkyBEeeOABQqEQoVCIt99+m5tuuokDBw5Euv6vvvpqQhYVfOaZZ3j11VdZtmwZixYt4pprruGhhx5KmoKJ9fX1/OIXvyAYDNLQ0MDLL7/ML3/5yzYLRiaSq6++mnXr1lFXV4dlWbz77rtMmzYtKX43m+3atYtBgwbh9YYP8xk7dmzStK+oqIgNGzbg8/lQSrFmzRouv/zypPj7O3LkCN/4xjcwTZP6+npeeukl5s+f36m2JUSPYsqUKWzdupXrr78eh8PB1KlTmTlzJjk5OXzrW98iGAwyZcoUpk2bFutQu8wjjzzSomDiggULYh1Sp1x99dVs2bKF66+/Htu2ueWWWxg3bhz33HMPCxYsiBSMHDNmTKxD7ZSxY8fyla98hVtuuQXDMLjqqqu4+eabGTJkSNL8bh4+fJjCwsLIZY/Hw89//vOkaN/EiRPZsWMH8+bNw+VyMXr0aL761a9y3XXXJfzfX1FREVOnTmXOnDlYlsXtt9/OuHHjOvXeIkUBhRBCtCkhhp6EEELEjiQKIYQQbZJEIYQQok2SKIQQQrRJEoUQQog2JcTyWCHa48EHH2TTpk1AeJNm3759SUlJAWDp0qWR7+ORUoo77riDxx57LFLKRYh4IctjRVK65ppr+PWvf83o0aNjHUq7mKbJyJEj2bRpkyQKEXekRyF6hD179vDTn/40soP69ttv54YbbmDDhg385je/IT8/n08//RSv18tXvvIVnn/+eQ4ePMj06dO599572bBhA4899hgFBQUcOHAAr9fLQw89xJAhQwiFQvziF79g8+bNWJbFyJEjuf/++0lPT2fy5MmMGzeOkpISvvvd72LbNk8//TShUIiqqipuvPFGvvWtb/G9730PgFtvvZWnn36aL3zhCzz55JMMHz4cgMmTJ/Pkk0/i9Xq54447GDBgAKWlpSxZsoQDBw7wP//zPwQCAXRdZ9GiRUyZMiWW/7tFsun6KuhCxN7VV1+ttm7dqpRSKhQKqenTp6udO3cqpZSqra1VxcXFauvWrWr9+vVqxIgRkdtuv/12dfPNN/8l0D4AAANHSURBVKtQKKQqKirU8OHDVUVFhVq/fr0qKipSmzdvVkop9fzzz6svfOELSimlfvWrX6lf/vKXyrZtpZRSDz/8sPrJT36ilFJq0qRJ6ne/+51SSinLstQtt9yiDh06pJRS6tixY6qoqEjV1NQowzDU0KFDVW1tbeRxO3bsiLSn+fLBgwfV0KFD1YcffqiUUqqqqkpNnTpVHT16VCmlVGlpqZo0aZIqLS2N0v9Z0RNJj0IkvX379nH48GHuvffeyHWhUIidO3fSr18/BgwYQFFREQD9+/cnLy8Pl8tFbm4uXq+XmpoaIFy2+ZJLLgHgC1/4Ag8++CD19fWsXbsWn8/Hu+++C4BhGC0KrY0bNw4AXdd58sknWbt2LcuWLWPv3r0opQgEAqSlpbW7PS6Xi7FjxwLw4YcfUl5ezte//vXI7bqus3v37hZlN4Q4F5IoRNKzbZusrCyWLVsWua68vJxevXqxefNm3G53i/s7na3/WZx6vW3bQPhN2bIsfvCDH3DVVVcB4aqdhmFE7tucBBoaGrjhhhsoLi5m3Lhx3Hjjjbz55puoVqYJNU1rcf2pz5eSkoKu65E4hg4dyl/+8pfI7WVlZZFqtkJ0BVkeK5LeBRdcgK7rvPbaa0D4BMFZs2ZRUlLSoefZvn07e/bsAcKrqC677DLS0tKYOHEizz//PIZhYFkW9913H7/61a9Oe/yBAwfw+/18+9vf5uqrr2bjxo2YpollWTgcDjRNi5xfnJOTw/bt24HwoUFVVVWtxnTxxRezb9++SDXQTz75hOLiYiorKzvUNiHaIj0KkfTcbjdPPPEEP/vZz/jd736HaZr8v//3/xg7diwbNmxo9/MUFBTwyCOPcPToUfLz83n44YcB+Na3vsXDDz/M9ddfH5nM/s///M/THj9ixAgmTpzI9OnTcblcFBUVMWTIEA4dOkTfvn2ZOnUqN998M7/97W/57ne/y49//GMWL17M6NGjI5Pan5WXl8djjz3GQw89RCgUQinFI488IsNOokvJ8lgh2mHDhg08/PDDLYavhOgpZOhJCCFEm6RHIYQQok3SoxBCCNEmSRRCCCHaJIlCCCFEmyRRCCGEaJMkCiGEEG2SRCGEEKJN/z9fvaNaIDS9HwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "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 }