{ "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.11 (default, Jul 27 2021, 09:42:29) [MSC v.1916 64 bit (AMD64)]\n", "uname_result(system='Windows', node='MATEIS502-R063', release='10', version='10.0.19041', machine='AMD64', processor='Intel64 Family 6 Model 142 Stepping 12, GenuineIntel')\n", "IPython 7.27.0\n", "IPython.core.release 7.27.0\n", "PIL 8.3.1\n", "PIL.Image 8.3.1\n", "PIL._version 8.3.1\n", "_csv 1.0\n", "_ctypes 1.1.0\n", "decimal 1.70\n", "argparse 1.1\n", "backcall 0.2.0\n", "bottleneck 1.3.2\n", "cffi 1.14.6\n", "colorama 0.4.4\n", "csv 1.0\n", "ctypes 1.1.0\n", "cycler 0.10.0\n", "dateutil 2.8.2\n", "decimal 1.70\n", "decorator 5.0.9\n", "defusedxml 0.7.1\n", "distutils 3.7.11\n", "entrypoints 0.3\n", "ipykernel 6.2.0\n", "ipykernel._version 6.2.0\n", "ipython_genutils 0.2.0\n", "ipython_genutils._version 0.2.0\n", "ipywidgets 7.6.4\n", "ipywidgets._version 7.6.4\n", "jedi 0.17.2\n", "joblib 1.0.1\n", "joblib.externals.cloudpickle 1.6.0\n", "joblib.externals.loky 2.9.0\n", "json 2.0.9\n", "jupyter_client 7.0.1\n", "jupyter_client._version 7.0.1\n", "jupyter_core 4.7.1\n", "jupyter_core.version 4.7.1\n", "kiwisolver 1.3.1\n", "logging 0.5.1.2\n", "matplotlib 3.4.2\n", "mkl 2.4.0\n", "numexpr 2.7.3\n", "numpy 1.20.3\n", "numpy.core 1.20.3\n", "numpy.core._multiarray_umath 3.1\n", "numpy.lib 1.20.3\n", "numpy.linalg._umath_linalg 0.1.5\n", "pandas 1.3.2\n", "parso 0.7.0\n", "patsy 0.5.1\n", "patsy.version 0.5.1\n", "pickleshare 0.7.5\n", "platform 1.0.8\n", "prompt_toolkit 3.0.17\n", "psutil 5.8.0\n", "pygments 2.10.0\n", "pyparsing 2.4.7\n", "pytz 2021.1\n", "re 2.2.1\n", "scipy 1.7.1\n", "scipy._lib._uarray 0.5.1+49.g4c3f1d7.scipy\n", "scipy._lib.decorator 4.0.5\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._fblas b'$Revision: $'\n", "scipy.linalg._flapack b'$Revision: $'\n", "scipy.linalg._flinalg b'$Revision: $'\n", "scipy.linalg._interpolative b'$Revision: $'\n", "scipy.ndimage 2.0\n", "scipy.optimize.__nnls b'$Revision: $'\n", "scipy.optimize._cobyla b'$Revision: $'\n", "scipy.optimize._lbfgsb b'$Revision: $'\n", "scipy.optimize._minpack 1.10 \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.11.2\n", "seaborn.external.husl 2.1.0\n", "six 1.16.0\n", "statsmodels 0.12.2\n", "statsmodels.__init__ 0.12.2\n", "statsmodels.api 0.12.2\n", "statsmodels.tools.web 0.12.2\n", "traitlets 5.1.0\n", "traitlets._version 5.1.0\n", "urllib.request 3.7\n", "wcwidth 0.2.5\n", "zlib 1.0\n", "zmq 22.2.1\n", "zmq.sugar 22.2.1\n", "zmq.sugar.version 22.2.1\n" ] } ], "source": [ "def print_imported_modules():\n", " import sys\n", " for name, val in sorted(sys.modules.items()):\n", " if(hasattr(val, '__version__')): \n", " print(val.__name__, val.__version__)\n", "# else:\n", "# print(val.__name__, \"(unknown version)\")\n", "def print_sys_info():\n", " import sys\n", " import platform\n", " print(sys.version)\n", " print(platform.uname())\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import seaborn as sns\n", "\n", "print_sys_info()\n", "print_imported_modules()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading and inspecting data\n", "Let's start by reading data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \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": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_url = \"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\"\n", "data = pd.read_csv(\"data_shuttle.csv\", delimiter=\",\")\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX4UlEQVR4nO3df5RcZX3H8fdnN0tI2Agx2K3NjwKSolQikuVX8ccGqw2cI6kHWtEKlhYjLbGVtsdQ22Op1p6CVSyWGiOlGHo0CqEQbVqE2uVHNZJQYyBQcEsw2YQGWANkQ9jsZr/9Y+7qZDO7O7PZO7Mzz+d1Ts7OvfeZZ7/P3p189v6YZxQRmJlZuppqXYCZmdWWg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHG5BYGkmyU9K+nREbZL0g2SuiRtlnRaXrWYmdnI8jwiuAVYPMr284D52b+lwBdzrMXMzEaQWxBExP3AT0ZpsgRYFQXrgWMkvTaveszMrLQpNfzes4HtRcvd2bpnhjeUtJTCUQPTpk1bOHfu3KoUWK7BwUGamhrzckujjs3jqj+NOrZqjevJJ598PiJeU2pbLYNAJdaVnO8iIlYCKwHa29tj48aNedZVsc7OTjo6OmpdRi4adWweV/1p1LFVa1ySfjzStlrGazdQ/Kf9HGBnjWoxM0tWLYNgLXBpdvfQWcCLEXHIaSEzM8tXbqeGJH0N6ACOldQN/AXQAhARK4B1wPlAF/AycFletZiZ2chyC4KIeN8Y2wO4Mq/vb2Zm5Wm8S/BmZlYRB4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmics1CCQtlvSEpC5JV5fYfrSkb0r6oaQtki7Lsx4zMztUbkEgqRm4ETgPOBl4n6SThzW7EngsIt4EdACflXREXjWZmdmh8jwiOAPoioinImI/sBpYMqxNADMkCWgFfgIM5FiTmZkNo4jIp2PpImBxRFyeLV8CnBkRy4razADWAq8HZgDvjYh/LdHXUmApQFtb28LVq1fnUvN49fb20traWusyctGoY/O46k+jjq1a41q0aNHDEdFeatuUHL+vSqwbnjq/BmwCzgVeB9wj6YGIeOmgJ0WsBFYCtLe3R0dHx4QXezg6OzuZbDVNlEYdm8dVfxp1bJNhXHmeGuoG5hYtzwF2DmtzGXBHFHQBWykcHZiZWZXkGQQbgPmSjs8uAF9M4TRQsW3AOwAktQEnAU/lWJOZmQ2T26mhiBiQtAy4G2gGbo6ILZKuyLavAD4F3CLpEQqnkpZHxPN51WRmZofK8xoBEbEOWDds3YqixzuBd+VZg5mZjc7vLDYzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscbkGgaTFkp6Q1CXp6hHadEjaJGmLpPvyrMfMzA41Ja+OJTUDNwLvBLqBDZLWRsRjRW2OAf4BWBwR2yT9XF71mJlZaWUdEUh64zj6PgPoioinImI/sBpYMqzN+4E7ImIbQEQ8O47vY2Zmh0ERMXYj6UHgCOAW4KsR8UIZz7mIwl/6l2fLlwBnRsSyojafB1qAXwZmAH8XEatK9LUUWArQ1ta2cPXq1WPWXE29vb20trbWuoxcNOrYPK7606hjq9a4Fi1a9HBEtJfaVtapoYh4i6T5wO8AGyU9BPxTRNwzytNUqqsS338h8A5gGvA9Sesj4slh338lsBKgvb09Ojo6yim7ajo7O5lsNU2URh2bx1V/GnVsk2FcZV8jiIgfSfpzYCNwA/BmSQI+HhF3lHhKNzC3aHkOsLNEm+cjYi+wV9L9wJuAJzEzs6oo9xrBAknXA48D5wLvjog3ZI+vH+FpG4D5ko6XdARwMbB2WJu7gLdKmiJpOnBm9j3MzKxKyj0i+HvgyxT++t83tDIidmZHCYeIiAFJy4C7gWbg5ojYIumKbPuKiHhc0r8Dm4FB4KaIePQwxmNmZhUqNwjOB/ZFxAEASU3AkRHxckTcOtKTImIdsG7YuhXDlj8DfKaiqs3MbMKU+4ayeylczB0yPVtnZmZ1rtwgODIieocWssfT8ynJzMyqqdwg2CvptKEFSQuBfaO0NzOzOlHuNYKPArdJGrr987XAe3OpyMzMqqrcN5RtkPR64CQKbxT7n4joz7UyMzOrikomnTsdOC57zpslUWo6CDMzqy9lBYGkW4HXAZuAA9nqABwEZmZ1rtwjgnbg5ChnhjozM6sr5d419Cjw83kWYmZmtVHuEcGxwGPZrKN9Qysj4oJcqjIzs6opNwiuybMIMzOrnXJvH71P0i8C8yPi3mym0OZ8SzMzs2oodxrqDwG3A1/KVs0G7sypJjMzq6JyLxZfCZwDvASFD6kB/EHzZmYNoNwg6Ms+gB4ASVM49GMnzcysDpUbBPdJ+jgwTdI7gduAb+ZXlpmZVUu5QXA18BzwCPBhCh82U/KTyczMrL6Ue9fQIIWPqvxyvuWYmVm1lTvX0FZKXBOIiBMmvCIzM6uqSuYaGnIk8BvAqye+HDMzq7ayrhFERE/Rvx0R8Xng3HxLMzOzaij31NBpRYtNFI4QZuRSkZmZVVW5p4Y+W/R4AHga+M0Jr8bMzKqu3LuGFuVdiJmZ1Ua5p4b+aLTtEfG5iSnHzMyqrZK7hk4H1mbL7wbuB7bnUZSZmVVPJR9Mc1pE7AGQdA1wW0RcnldhZmZWHeVOMTEP2F+0vB84bsKrMTOzqiv3iOBW4CFJ/0LhHcbvAVblVpWZmVVNuXcNfVrSvwFvzVZdFhE/yK8sMzOrlnJPDQFMB16KiL8DuiUdn1NNZmZWReV+VOVfAMuBP81WtQD/nFdRZmZWPeUeEbwHuADYCxARO/EUE2ZmDaHcINgfEUE2FbWko/IryczMqqncIPiGpC8Bx0j6EHAv/pAaM7OGMGYQSBLwdeB2YA1wEvCJiPhCGc9dLOkJSV2Srh6l3emSDki6qILazcxsAox5+2hEhKQ7I2IhcE+5HUtqBm4E3gl0AxskrY2Ix0q0uxa4u6LKzcxsQpR7ami9pNMr7PsMoCsinoqI/cBqYEmJdh+hcKTxbIX9m5nZBCj3ncWLgCskPU3hziFROFhYMMpzZnPwpHTdwJnFDSTNpnBH0rkUJrUrSdJSYClAW1sbnZ2dZZZdHb29vZOuponSqGPzuOpPo45tMoxr1CCQNC8itgHnjaNvlVgXw5Y/DyyPiAOFSxGlRcRKYCVAe3t7dHR0jKOc/HR2djLZapoojTo2j6v+NOrYJsO4xjoiuJPCrKM/lrQmIi6soO9uYG7R8hxg57A27cDqLASOBc6XNBARd1bwfczM7DCMFQTFf6afUGHfG4D52VQUO4CLgfcXN4iIn05TIekW4FsOATOz6horCGKEx2OKiAFJyyjcDdQM3BwRWyRdkW1fUVGlZmaWi7GC4E2SXqJwZDAteww/u1j8qtGeHBHrgHXD1pUMgIj47bIqNjOzCTVqEEREc7UKMTOz2qhkGmozM2tADgIzs8Q5CMzMEucgMDNLXDJB0NPbxw+3v0BPb1+tSzGzCvX09rGv/4BfvzlJIgju2rSDc679Dh+46fucc+13WLtpR61LMrMyDb1+tz6316/fnDR8EPT09rF8zWZe6R9kT98Ar/QP8rE1m/2XhVkdKH79Hojw6zcnDR8E3bv30dJ08DBbmpro3r2vRhWZWbn8+q2Ohg+COTOn0T84eNC6/sFB5sycVqOKzKxcfv1WR8MHwazWqVx34QKObGlixtQpHNnSxHUXLmBW69Ral2ZmYyh+/TZLfv3mpNwPpqlrF5w6m3NOPJbu3fuYM3Oaf4nM6sjQ6/eh7z3If13wFr9+c5BEEEDhLwv/ApnVp1mtU5nW0uzXcE4a/tSQmZmNzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklLtcgkLRY0hOSuiRdXWL7b0nanP37rqQ35VmPmZkdKrcgkNQM3AicB5wMvE/SycOabQXeHhELgE8BK/Oqx8zMSsvziOAMoCsinoqI/cBqYElxg4j4bkTszhbXA3NyrMfMzEpQROTTsXQRsDgiLs+WLwHOjIhlI7T/E+D1Q+2HbVsKLAVoa2tbuHr16lxqHq/e3l5aW1trXUYuGnVsHlf9adSxVWtcixYtejgi2kttm5Lj91WJdSVTR9Ii4HeBt5TaHhEryU4btbe3R0dHxwSVODE6OzuZbDVNlEYdm8dVfxp1bJNhXHkGQTcwt2h5DrBzeCNJC4CbgPMioifHeszMrIQ8rxFsAOZLOl7SEcDFwNriBpLmAXcAl0TEkznWYmZmI8jtiCAiBiQtA+4GmoGbI2KLpCuy7SuATwCzgH+QBDAw0jksMzPLR56nhoiIdcC6YetWFD2+HDjk4rBBT28f3bv3MWfmNGa1Tp2wtvWkUceVl65de9j9cj9du/ZwYtuMWpdjdSTXILDxuWvTDpav2UxLUxP9g4Ncd+ECLjh19mG3rSeNOq68fOLOR1i1fht/fMoAV11/P5eePY9PLjml1mVZnfAUE5NMT28fy9ds5pX+Qfb0DfBK/yAfW7OZnt6+w2pbTxp1XHnp2rWHVeu3HbRu1fe20bVrT40qsnrjIJhkunfvo6Xp4N3S0tRE9+59h9W2njTquPKyafsLFa03G85BMMnMmTmN/sHBg9b1Dw4yZ+a0w2pbTxp1XHk5de4xFa03G85BMMnMap3KdRcu4MiWJmZMncKRLU1cd+GCkhdLK2lbTxp1XHk5sW0Gl54976B1l549zxeMrWy+WDwJXXDqbM458diy7pippG09adRx5eWTS07h0rOO45GH13PvVWc5BKwiDoJJalbr1LL/86ukbT1p1HHl5cS2GXRPb3EIWMV8asjMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcrkEgabGkJyR1Sbq6xHZJuiHbvlnSaXnWY1apnt4+frj9BXp6+8Zsu3FrD5/79hNs3NozYX1W0rZr1x52v9xP1649Y7atRF71VlrDvv4DY/bbtWsPt2/c3rA/gzz6BZgy4T1mJDUDNwLvBLqBDZLWRsRjRc3OA+Zn/84Evph9Nau5uzbtYPmazbQ0NdE/OMh1Fy7gglNnl2z7gZvW82BXIQBu+E4Xbz1xFrdeftZh9VlJ20/c+Qir1m/jj08Z4Krr7+fSs+fxySWnjHPk+dc7nhr+4A39XHXtd0bsd+hnMKQRfwYT3e+QPI8IzgC6IuKpiNgPrAaWDGuzBFgVBeuBYyS9NseazMrS09vH8jWbeaV/kD19A7zSP8jH1mwu+dfYxq09Pw2BIQ909RxyZFBJn5W07dq156D/AAFWfW/bYf9VnFe9463hQMSI/abyM5jIfospIiass4M6li4CFkfE5dnyJcCZEbGsqM23gL+JiAez5f8AlkfExmF9LQWWZosnAU/kUvT4HQs8X+sictKoYxt1XGqZNn3KzNf+kpqamofWxeDggYHdzzwZ/fteLm7bPOPYX2g+6phD/oA5sPeFZw7seX7nePqspG3T9KNnTXnVa44DOPDyizRPPxqAgZeee3rw5RdHP081irzqHW8NQ2Mr1W/xz6BYnfwMJux3cQy/GBGvKbUht1NDgEqsG5465bQhIlYCKyeiqDxI2hgR7bWuIw+NOrZGHtfAi8823Ligccc2GX4X8zw11A3MLVqeA+wcRxszM8tRnkGwAZgv6XhJRwAXA2uHtVkLXJrdPXQW8GJEPJNjTWZmNkxup4YiYkDSMuBuoBm4OSK2SLoi274CWAecD3QBLwOX5VVPzibtaasJ0Khj87jqT6OOrebjyu1isZmZ1Qe/s9jMLHEOAjOzxDkIxkHS05IekbRJ0sZs3TWSdmTrNkk6v9Z1VkrSMZJul/Q/kh6XdLakV0u6R9KPsq8za11npUYYVyPsr5OK6t8k6SVJH633fTbKuBphn10laYukRyV9TdKRk2F/+RrBOEh6GmiPiOeL1l0D9EbE39aqrsMl6SvAAxFxU3an13Tg48BPIuJvsvmiZkbE8poWWqERxvVR6nx/FcumdNlBYYqWK6nzfTZk2Lguo473maTZwIPAyRGxT9I3KNwwczI13l8+IjAAJL0KeBvwjwARsT8iXqAwDchXsmZfAX69FvWN1yjjajTvAP43In5Mne+zYYrH1QimANMkTaHwB8lOJsH+chCMTwDflvRwNv3FkGXZLKo319vhOHAC8BzwT5J+IOkmSUcBbUPv7ci+/lwtixyHkcYF9b2/hrsY+Fr2uN73WbHicUEd77OI2AH8LbANeIbC+6a+zSTYXw6C8TknIk6jMHvqlZLeRmHm1NcBp1LYyZ+tXXnjMgU4DfhiRLwZ2AscMnV4HRppXPW+v34qO911AXBbrWuZSCXGVdf7LAuuJcDxwC8AR0n6QG2rKnAQjENE7My+Pgv8C3BGROyKiAMRMQh8mcLsq/WkG+iOiO9ny7dT+A9019CMsNnXZ2tU33iVHFcD7K9i5wH/HRG7suV632dDDhpXA+yzXwW2RsRzEdEP3AH8CpNgfzkIKiTpKEkzhh4D7wIe1cHTZ78HeLQW9Y1XRPwfsF3SSdmqdwCPUZgG5IPZug8Cd9WgvHEbaVz1vr+GeR8Hnz6p631W5KBxNcA+2wacJWm6JFH4XXycSbC/fNdQhSSdQOEoAAqnHb4aEZ+WdCuFQ9YAngY+XG/zJkk6FbgJOAJ4isJdGk3AN4B5FH6RfyMiflKrGsdjhHHdQJ3vLwBJ04HtwAkR8WK2bhb1v89KjasRXmN/CbwXGAB+AFwOtFLj/eUgMDNLnE8NmZklzkFgZpY4B4GZWeIcBGZmiXMQmJklLs8PrzeruuzWyf/IFn8eOEBhigkovPFvf00KK0FSB7A/Ir5b41IscQ4CaygR0UPhXvNJMSOspCkRMTDC5g6gFyg7CCQ1R8SBiajNbIhPDVnDk7RQ0n3ZJIF3F72dv1PS9ZLuzz6n4HRJd2Tzwv9V1ua47HMMvpJNdnZ79mansfr9a0n3AX8o6d2Svp9NenevpDZJxwFXAFdlc+u/VdItki4qqrs3+9oh6T8lfRV4RFKzpM9I2pDV9OGq/kCt4TgIrNEJ+AJwUUQsBG4GPl20fX9EvA1YQeGt/VcCbwR+OzvNBHASsDIiFgAvAb8vqWWMfo+JiLdHxGcpzEF/Vjbp3WrgYxHxdPY9r4+IUyPigTHGcQbwZxFxMvC7FGauPB04HfiQpOMr/9GYFfjUkDW6qRT+Y7+nML0LzRRmrhyyNvv6CLBlaMoCSU8Bc4EXgO0R8V9Zu38G/gD49zH6/XrR4znA17MjhiOAreMYx0MRMfS8dwELio4ejgbmj7NfMweBNTxR+A/+7BG292VfB4seDy0PvT6Gz8MSZfS7t+jxF4DPRcTa7ALxNSM8Z4DsKD2blOyIEfoT8JGIuHuEfswq4lND1uj6gNdIOhtAUoukX66wj3lDz6cwI+aDwBMV9Hs0hY9bhJ/NMgmwB5hRtPw0sDB7vARoGaG/u4Hfy05PIemXij5sx6xiDgJrdIPARcC1kn4IbKIwB3wlHgc+KGkz8GoKH3Kzv4J+rwFuk/QA8HzR+m8C7xm6WExhjv23S3qIwmf07j2kp4KbKEwR/t+SHgW+hI/u7TB49lGzUWR393wrIt5Y61rM8uIjAjOzxPmIwMwscT4iMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNL3P8Dn1j0gsgASdcAAAAASUVORK5CYII=\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": 4, "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", "
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, 24 Jan 2022 Deviance: 3.0144
Time: 16:25:21 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, 24 Jan 2022 Deviance: 3.0144\n", "Time: 16:25:21 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", "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import statsmodels.api as sm\n", "\n", "data[\"Success\"]=data.Count-data.Malfunction\n", "data[\"Intercept\"]=1\n", "\n", "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", " family=sm.families.Binomial(sm.families.links.logit())).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." ] }, { "cell_type": "code", "execution_count": 26, "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", "
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, 24 Jan 2022 Deviance: 18.086
Time: 16:57:54 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, 24 Jan 2022 Deviance: 18.086\n", "Time: 16:57:54 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", "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": 26, "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": 28, "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", "
TemperatureInterceptFrequency
030.011.0
130.511.0
231.011.0
331.511.0
432.011.0
............
11688.011.0
11788.511.0
11889.011.0
11989.511.0
12090.011.0
\n", "

121 rows × 3 columns

\n", "
" ], "text/plain": [ " 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", ".. ... ... ...\n", "116 88.0 1 1.0\n", "117 88.5 1 1.0\n", "118 89.0 1 1.0\n", "119 89.5 1 1.0\n", "120 90.0 1 1.0\n", "\n", "[121 rows x 3 columns]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAENCAYAAAAbu05nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaCElEQVR4nO3de3xU5b3v8c9MYlA0XBpTIVhB7as/C14QRUSx2ILtrrdW1FNFrZdSa221r1Pt69AttVilu7irVfe21h51i+2mN7bt0QreCmKtqMXKpSq/002RIyS0GNGYyMUwc/5YK+kQQrImzDgzD9/3P8y6zXp+mfDNM8+s9Uwqm80iIiLhSJe6ASIiUlgKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwFQn2cnMBgDPAqe7+2tdto0G7gEGAE8DV7h7e2GbKSIiSfXaYzezccAzwEd2sctPga+6+0eAFPDFwjVPRETylWQo5ovAV4DGrhvMbDiwj7s/F6+6Hzi3YK0TEZG89ToU4+7TAMysu80NQFPOchNwYB7n7weMjY/bnsdxIiJ7sipgKPBHYGvXjYnG2HuQBnLnJEgBmTyOHwv8fjfbICKypzqJaKh8B7sb7OuI/mp0GEI3QzY9aALYtKmNTCb/OWvq6vajubk17+PKkWopP6HUAaqlXPW1lnQ6xeDB+8KOIyaddivY3X2tmW0xsxPd/Q/ARcCCPJ5iO0Amk+1TsHccGwrVUn5CqQNUS7nazVq6HcLu03XsZjbfzI6NFy8AfmBmq4D9gDv61j4RESmExD12dx+R8/jUnMfLgeMK2ywREemr3R1jF5E9UDabZdOmjWzbtoUdr58ovr//PU0mk881GuWr51pS1NTszeDB9aRSqbyeV8EuInlrbX2bVCrFAQccSCr1/s5MUl2dpr09jGDvqZZsNsNbb71Ba+vb1NYOyut5NVeMiORt8+ZWamsHve+hvidJpdLU1g5m8+Y+XDVThPaISOAyme1UVekNf7FVVVWTyeR/76aCXUT6JN9xX8lfX3/G+pMrIhWvqamR88+fwogRh+ywfvbsWznggCElalXpKNhFJAj771/P/ffPLXUzyoKCXUSCNWvWTN5++23Wr3+dL3/5aurq6rjjjlvZunULAwcO4hvf+GcaGoaxatWr3HzzTQCMG3cCTzzxKPPmPcysWTM5+uhjOPXUMwCYMOFYnnlmKe+++y633jqbv/51NZlMhgsu+DynnPJPzJ//MM8//ywtLS00Nq5n7Njjufba6WSzWe666994+umnqK6u4swzp3DCCRP42te+zK9+9RDpdJo//Wkp//mfD3DLLbt/j6eCXUR2yx9WNvHMim6nLNltE44cyolHDO19R+CNNzZyySVTO5c/+cl/AmDgwIHcfPMPeO+995g27fPMnv0DhgwZwvPPL2H27FncfvsPuemm67nqqq8zbtx47rvvx72ea86cezH7KDNm3EBbWytXXHEZI0ceDsDKlSv46U9/STpdxdSpZ7N69TmsXfsaK1cu54EHfk57eztXXjmNSZNOoaGhgZdeepFjjhnLo48+wqmnnt6Hn9LOFOwiEoTuhmJmzZrZGbivv76WxsZ1TJ/+9c7tbW1tbNq0iebmZsaNGw/AmWeexfz5D/d4rqVLX2Dr1i088shDAGzZsoU1a/4KwBFHHEn//vsC0NAwjJaWt1m27EU+8YlTqKmpoaamprOdp5/+GR57bD6jRh3Biy/+kWuumV6An4SCXUR204lHJO9Vl0K/fv0A2L49Q0PDsM5Q3b59O5s2vdm5vUPuZZypVIpsNrqztr39H9/4mcls51vfuhGzwwB4881mBgwYyOOPL6CmZsfny2azVFdXk3uBS1NTI4MGDWbSpMn86Ed3smjRk4wff+JObekrXe4oInuE4cNH0NLSwvLlLwHwyCMPMXPmdfTv35/hw0fwzDOLAXjiiUc7jxk4cFBnT/zpp5/qXD9mzFh+85t5ALzxxhtcfPH5/O1vG3Z57qOOGsNTTy2kvb2dLVu2cM01V7Fx49/Ze+99OP74E/jxj3/Ipz99RsFqVY9dRPYINTU13Hjj97j99u+zbds2+vfflxkzbgBgxowb+N73buTee+9m+PCDO4/57GfP5vrrv8nFF5/HmDFjqavbH4DLLvsit9wym4su+h9kMhmuvPJqhg07sPOPRlcTJ36cVate4bLLLiCTyXLuuedz0EHDAZg06ZOsXLmcUaMOL1itqY63GSUyAljT3NzapzmJ6+tr2bjxnYI3qhRUS/kJpQ4ofC0bNqxlyJDhBXu+fBR7rpimpkauuupLzJvX8zh7IaRSWX74w39n8ODBnHfehd3u093POp1OUVe3H8DBwGtdj9FQjIhIiVx66YW4v8pZZ51T0OfVUIyISI6hQxvel946wAMP/Kwo7z7UYxcRCYyCXUT6pMSfz+0R+vozVrCLSN6qq2toa2tRuBdRNpulra2F6uqavI/VGLuI5G3w4Ho2bdpIa+tb7/u50+lwvhqvt1qqq2sYPLg+7+dVsItI3qqqqtl//9LcbarLUHunoRgRkcAo2EVEAqNgFxEJjIJdRCQwCnYRkcAo2EVEAqNgFxEJjIJdRCQwCnYRkcAo2EVEAqNgFxEJjIJdRCQwCnYRkcAo2EVEApNo2l4zmwrMAPYCbnP3O7tsHwPcDdQArwMXuvtbhW2qiIgk0WuP3cyGAbOACcBo4HIzG9llt9uB6939KMCBawvcThERSSjJUMxkYKG7v+nubcA84Jwu+1QBA+LH/YHNhWuiiIjkI8lQTAPQlLPcBBzXZZ+vA4+b2W1AGzAun0bU1e2Xz+47qK+v7fOx5Ua1lJ9Q6gDVUq6KUUuSYE8Dud9YmwI6v6TPzPYB7gUmu/sLZvZ14AHgtKSNaG5uJZPJ/0tx9RVZ5SmUWkKpA1RLueprLel0qscOcZKhmHVA7pcbDgEac5YPBza7+wvx8t3Ayfk1U0RECiVJsD8JTDKzejPrD5wNPJqz/b+BD5mZxcufAf5Y2GaKiEhSvQa7u68HrgMWAcuAufGQy3wzO9bdNwGXAL80sxXAZcClxWuyiIj0JNF17O4+F5jbZd2pOY8XAAsK2zQREekL3XkqIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGAU7CIigVGwi4gERsEuIhIYBbuISGCqk+xkZlOBGcBewG3ufmeX7QbcDQwGNgDnufumArdVREQS6LXHbmbDgFnABGA0cLmZjczZngIeAr7n7kcBLwHTi9JaERHpVZKhmMnAQnd/093bgHnAOTnbxwBt7v5ovPxd4E5ERKQkkgzFNABNOctNwHE5yx8GNpjZvcDRwKvAVQVroYiI5CVJsKeBbM5yCsh0eY6TgY+5+1IzuxG4FbgkaSPq6vZLuutO6utr+3xsuVEt5SeUOkC1lKti1JIk2NcBJ+UsDwEac5Y3AH9x96Xx8s+IhmsSa25uJZPJ9r5jF/X1tWzc+E7ex5Uj1VJ+QqkDVEu56mst6XSqxw5xkjH2J4FJZlZvZv2Bs4FHc7Y/C9Sb2VHx8hnAi3m3VERECqLXYHf39cB1wCJgGTDX3V8ws/lmdqy7bwbOAv63mb0MfAK4pohtFhGRHiS6jt3d5wJzu6w7Nefx8+z4gaqIiJSI7jwVEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCjYRUQCo2AXEQmMgl1EJDAKdhGRwCQKdjObamavmNlfzOwrPex3mpmtKVzzREQkX70Gu5kNA2YBE4DRwOVmNrKb/Q4Avg+kCtxGERHJQ5Ie+2Rgobu/6e5twDzgnG72uwe4oZCNExGR/FUn2KcBaMpZbgKOy93BzK4G/gQ815dG1NXt15fDAKivr+3zseVGtZSfUOoA1VKuilFLkmBPA9mc5RSQ6Vgws8OBs4FJwIF9aURzcyuZTLb3Hbuor69l48Z3+nLKsqNayk8odYBqKVd9rSWdTvXYIU4yFLMOGJqzPARozFk+N96+FJgPNJjZ7/NuqYiIFESSHvuTwEwzqwfaiHrnl3dsdPdvA98GMLMRwFPuflLhmyoiIkn02mN39/XAdcAiYBkw191fMLP5ZnZskdsnIiJ5StJjx93nAnO7rDu1m/1eA0YUomEiItI3uvNURCQwCnYRkcAo2EVEAqNgFxEJjIJdRCQwCnYRkcAkutxRpNiWvLyBBxevprllK3UD+jFl4qGMHzWk1M2ShPT6lRcFu5Tckpc3MGfBKra1R1MQNbdsZc6CVQAKhwqg16/8aChGSu7Bxas7Q6HDtvYMDy5eXaIWST70+pUfBbuUXHPL1rzWS3nR61d+FOxScnUD+uW1XsqLXr/yo2CXkpsy8VBqqnf8VaypTjNl4qElapHkQ69f+dGHp1JyHR+w6aqKyqTXr/wo2KUsjB81REFQwfT6lRcNxYiIBEbBLiISGAW7iEhgFOwiIoFRsIuIBEbBLiISGAW7iEhgFOwiIoFRsIuIBEbBLiISGAW7iEhgFOwiIoFRsIuIBEbBLiISGAW7iEhgFOwiIoFRsIuIBEbBLiISGAW7iEhgFOwiIoFJ9GXWZjYVmAHsBdzm7nd22f4Z4AYgBawBLnX3TQVuq4iIJNBrj93MhgGzgAnAaOByMxuZs30AcBdwmrsfBawAZhajsSIi0rskQzGTgYXu/qa7twHzgHNytu8FfMXd18fLK4CDCttMERFJKslQTAPQlLPcBBzXseDuzcCvAcxsH2A68G8FbKOIiOQhSbCngWzOcgrIdN3JzAYSBfxyd5+TTyPq6vbLZ/cd1NfX9vnYcqNayk8odYBqKVfFqCVJsK8DTspZHgI05u5gZkOBx4CFwP/MtxHNza1kMtned+yivr6WjRvfyfu4cqRayk8odYBqKVd9rSWdTvXYIU4S7E8CM82sHmgDzgYu79hoZlXAw8Av3f2mvFsoIiIF1Wuwu/t6M7sOWATUAPe4+wtmNh+4HvgQMAaoNrOOD1WXuvu0YjVaRER2LdF17O4+F5jbZd2p8cOl6EYnKTNLXt7Ag4tX09yylboB/Zgy8VCAndaNHzXkfTl3Mc6TxE8eW8XiZY1kspBOwcTRDVz0qcNK0hZ5/yQKdpFKsuTlDcxZsIpt7dFn/M0tW7nvt6+QSqdo357tXDdnwSqAgoZud+cuxnmS+Mljq1j00j8+Dstk6VxWuIdNPW0JzoOLV3cGa4ftWTpDvcO29gwPLl5d9HMX4zxJLF7WmNd6CYeCXYLT3LK1KPvuzvMV+jxJ7OpCsz5cgCYVRsEuwakb0K8o++7O8xX6PEmkU/mtl3Ao2CU4UyYeSk31jr/aVSmortox0Wqq050fqhbz3MU4TxITRzfktV7CoQ9PJTgdH1KW4qqYXZ27FFfFdHxAqqti9jypbLakA24jgDW681S1lKNQ6gDVUq4KcOfpwcBrO23f7ZaJiEhZUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiARGwS4iEhgFu4hIYKqT7GRmU4EZwF7Abe5+Z5fto4F7gAHA08AV7t5e2KaKhG/Jyxt4cPFqmlu2UjegH1MmHsozKxp5de1bnft8dPggJhzZsNN+wE7r/nvdWyxe1kgmC+kUTBzdwEWfOizRebt7vvGjhiRu9/tx7u6O3VUb9ySpbDbb4w5mNgx4BjgG2Ao8C5zv7q/k7PNnYJq7P2dm9wJL3f2uBOcfAaxpbm4lk+m5Hd2pr69l48Z38j6uHKmW8vN+17Hk5Q3MWbCKbe2ZvI+trkqRzWTZnvPfKAV097/q40fvGLDdnbcqBal0ivacJ6ypTnPxpw/bKTi7O76Y5+54Xbo7dldtLFd9/R1Lp1PU1e0HcDDwWtftSXrsk4GF7v4mgJnNA84BvhMvDwf2cffn4v3vB24AkgR7VUcj+2p3ji03qqX8vJ91PL2skUG1/Yp+nlfXbtqhrnzO+/SyRk48YuhO65IeX6hzp9OpXR7bXRvLWV9+x3KOqepue5Ie+zeBfd19Rrw8DTjO3S+Pl8cD/+ruE+LlDwPz3f0jCdo3Afh9gv1ERGRnJxGNqOwgSY89zY7vqlJAJo/tPflj3LAmYHvCY0RE9nRVwFCiDN1JkmBfRxS+HYYAjV22D+1he0+20s1fGxER6dXqXW1Icrnjk8AkM6s3s/7A2cCjHRvdfS2wxcxOjFddBCzYjcaKiMhu6DXY3X09cB2wCFgGzHX3F8xsvpkdG+92AfADM1sF7AfcUaT2iohIL3r98FRERCqL7jwVEQmMgl1EJDAKdhGRwCjYRUQCk2gSsHJhZt8hms4gC9zr7rea2WTgVmAf4Bcdd8hWAjP7PrC/u19SqXWY2SLgg8B78aovAbVUZi1nAN8G9gUed/evVeLrEt8d/tWcVQcDPwF+Q+XVciHwzXhxgbtfW4mvCYCZTQcuJbp/5xfuPqtYtVTMVTFmNhGYBZxMNMvkK8BngYeBicDrwCNEs0+W/XX0ZjYJ+DlRm78MOBVWh5mliG5QG94xm6eZ7UNl1nII0fQW44C/AQuB7wJ3U2G15DKzUUSB/gngD1RQLfF9M+uAjwBvEbX/JuBOKqgOgJwAnwC0Ab8GfgH8C0WopWKGYtx9MfDxOEA+SPRuYxDwF3dfE6//KXBu6VqZjJl9gOiP1HfjVcdRgXUAFv/7uJktN7OvUrm1nEXUY1rn7u8BnwPepTJryXUX8M/AIVReLVVEGbUvUWduL6CFyqsD4GjgMXdvcfftRDd5TqNItVRMsAO4+3tmdgNRb/13QAPRPDMdmoADS9G2PN1NdNPXpni5UusYTPQ6nAVMAq4ADqIya/kwUGVmD5nZMuBKKvd1ATp7ifu4+6+owFrc/R3gW8Aqop77a1RgHbE/AZ8ysw+Y2d7AmUS996LUUlHBDuDu3wbqgQ8RvUXr6wRkJRGPf77u7r/LWb07E6mVjLsvcffPu/vb7v4GcC/RdM4VVwvRO8DJwBeA8URDModQmbV0+BLR23+owN8xMzsSuAwYThTo26nA//MA8f/3+4GniHrrzxD9zhWllor58NTMDgP2dvdl7v6umT1I9EFq7qyQ+UxAViqfA4bGvcIPEE3BMJzKqwMzmwD0y/kjlSLqVfV1UrhS2gA86e4bAczs10RviyvudQEwsxqisdtL4lW7M1lfqXwK+J27/x3AzO4HrqUCXxMzqwX+y91vjZe/QRTyRXlNKibYiXpPN8RhkgU+QzSk8a/xHPBrgKnAfaVrYu/c/ZSOx2Z2CdGHwVcAf6mkOmKDgO+Y2QlE458XE9Xyywqs5bfAHDMbBLwDfBqYB0yvwFoAjgT+r7u3xcvPA1ZhtSwHbjazfYk+7ziDqI4LKqwOiK5MeiCeX2tfoneGXwB+XoxaKmYoxt3nE31q/BLwIvCsu/+cqEfyX0Tj7quI/jNWFHffQgXW4e6/ZcfX5D53X0Jl1vI8cDPRW+RXgLVEHzxeQoXVEjuEqJcOVObvmLs/DvyM6HdrBVHnYSYVVgeAu68gavMK4AWiq1/+QJFqqZjLHUVEJJmK6bGLiEgyCnYRkcAo2EVEAqNgFxEJjIJdRCQwlXQdu+xBzOwO4GPx4kii63w3x8vj3X1ztweWmJmNBb7g7leUui2y51KwS1ly96s7HpvZa8AF7r60dC1KbBSVMXeJBEzBLhXFzL5ANEFXGmgGvuruq+Lbzd8FjgAOAB6Kt59BdKv2NHdfGO+3GRhNNEvo48DV8QRzHwVuB+qIZha8w93vM7OT4/VtRFNAjCW6mel4ornnU0Qz9f0/orlyBprZfwBzgH9398Pjtp/csWxmM4nmpGkAlrv7hWZ2HXB2XNtrwJXuXva3y0v50Ri7VIx4Tv6LgZPc/WiicP11zi5jiOYd/xhwDdDq7icQhfL0nP3GAacQDfGMBL5kZtXEUwi4+zFE86xca2bHx8ccDpzv7kfG52kgGhIaSRTg0939deB64PfufmmCkoYDR8eh/nmiP0rHuftoYD5wT/Kfjsg/qMculeQ0oul1nzXrmAqewfH89gAPx3OpbzCzNqJZ9ABWE0241uF+d28FMLMHiL6wZSFwKHBfznPvQzSP9qtEM3KuhWhWSzObQfQH4VCi+X7e6UM9z3V8QQlwOtFc9kvj81cB/fvwnCIKdqkoVcBP3P1/AZhZmqjn3DGv/dYu+79H99pzHqeJZgusAt6Oe8vEz38A8DbRkEtrzvrTiN4F3AL8H6I5Pi7s5jxZomGaDjVdtrfmPK4CZrv7XfE5+hHNdy+SNw3FSCV5DDjfzDqmOr2C6Is+8vU5M+sXf+HBxURfr+jA5vg7NjGzDwF/Bo7p5vhTiN4d3AUsJerxV8Xb2okmqwLYCBxkZh+Mv0bwvF5qm2ZmA+Ll7xB9T6lI3hTsUjHi2f5mA0+Y2QqiaU6nuHu+M9m9S/T9pivjf//D3bcRTQU9LX7ux4FvxTPwdfUj4GQzW0n0zTirgYPjdxDPAYeY2YPu/grR1NJL4/VremjTPURTBz9nZi8TTbt7SZ51iQCa3VH2MPFVMX929++Xui0ixaIeu4hIYNRjFxEJjHrsIiKBUbCLiARGwS4iEhgFu4hIYBTsIiKBUbCLiATm/wN9qi5XvEgE7AAAAABJRU5ErkJggg==\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", "data_pred" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "The prediction from `logmodel.predict(data_pred)` does not give the expected results.\n", "\n", "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": 25, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAENCAYAAAARyyJwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABA5UlEQVR4nO3deXxU5b348c85s2chewiERUAIOyhVQRFX9k3BtogVV1rba6n2/tpapbWtWq3X1/Vq6/UKWpcKijuLNaAoKoK4E/ZddrLvme0svz8mGRIDcRIymSXf9+sVkplz5szzMJn55tm+j2KapokQQghxGmqkCyCEECK6SaAQQgjRIgkUQgghWiSBQgghRIskUAghhGiRBAohhBAtCmugqKmpYdq0aRw5cqTZsR07djBr1iwmTpzIPffcg6Zp4SyKEEKINgpboNi8eTPXXnst33777SmP/+Y3v+GPf/wjq1evxjRNXnnllXAVRQghxBkIW6B45ZVXuPfee8nOzm527OjRo3g8HkaOHAnArFmzyM/PD1dRhBBCnAFruC78wAMPnPZYUVERWVlZwdtZWVkUFhaGqyhCCCHOQEQGsw3DQFGU4G3TNJvcFkIIET3C1qJoSU5ODsXFxcHbJSUlp+yi+j77D5ehaUaT+xrCjaIEvlDUQDRUQFUUFBQUlUBgargv8OPJ7w3nEDgG0NEZsTIykigtrenYJ+1A8Vy/eK4bSP1imaoqpKUltvpxEQkUubm5OBwOvvzyS0aNGsXy5csZN25cq6/j8+n4vxMomtNDvl5DUFAIRIz6OINVUVFVBdWioKoKFjUQcFQFLBYFtf6B7R1MDCO+8zXGc/3iuW4g9etsOjRQzJ8/nwULFjBs2DAeeeQRFi5cSE1NDUOGDGHevHkdWZRTavigNzGh4ffEAO0UwUap/0dRAgHDqqpYrCoWi1IfWMBiUVGVjm+NCCFEe1JiOc34zn3FIbQoIkOpDyI2i4rNpmK1qNhUBYsltGGhrKxkiourw1zKyInn+sVz3UDqF8tUVSEjI6nVj4tI11NnYJqBQXqvoeP1B1okigIWRcFus2CzqdgtgRbIyZEVIYSIPhIoOpBpgmaaaF4NvPWBw6LgsFlx2CzYrKqEDCFE1JFAEUGmCZpmoml+6jx+LIqCwxEIGjHcIyiEiDMSKKJEsLXhDgQNq92Kx6vhtFuCs6qEECISJFBEIdMEzTCpqvVR41ZwOay4HBZsFlVmUAkhOpwEiihnGCa19a0Mp91KosuGVZUWhhCi40igiBGmCW6vhsen4XJYSXJaUVXZTkQIEX7ySRNjTBPqPBollV5qPf5IF0cI0QlIoIhRhmlSXeentMqNz28g491CiHCRQBHj/JpJeY2HylpfIPWIEEK0MwkUcaChO6qsyotfj86UJkKI2CWBIo74NYPyKi91Xg2kdSGEaCcSKOKMYQbWX1TW+jFk0YUQoh1IoIhTbq9GebUHTZdgIYQ4MxIo4phfMymr9uDxh755kxBCfJcEijhnGCaVNV5qvZpMoRVCtIkEik7ANKGm1kdVnV+GuIUQrSaBopMwgVq3n6panwQLIUSrSK6nTsbt1QDokmCXrighREikRdEJub0alXVeSVkuhAiJBIpOyuPVA8Ei0gURQkQ9CRSdmMerU+P2SxeUEKJFEig6uTq3n1qPFuliCCGimASKTs4Equt8eGVRnhDiNCRQCEwTKmt9aIaMWAghmpNAIYDACu6qWl+kiyGEiEISKESQz69T7ZbtVYUQTUmgEE3Uefz4/LL5kRDiJAkUognThKpar+xlIYQIkkAhmtEMk+o6WV8hhAiI6UCx+rPDwdxFon15fBpen0yZFULEeKAo2FfC/7yymS37SzGlq6RdmSZUuX2YkuRDiE4vpgOFzapS7fbz0nt7WPruHqrrZHpne9I0kxq3tNiE6OxiOlDcNGUQeT1TAdj2bRmPvVpAwb4SaV20ozqPH78us6CE6MxiOlCkJNqZNymPH17aD6fdQp1X4+W1e1n63h5qPbIeoD2YJlTLznhCdGphDRQrV65kypQpTJgwgSVLljQ7vm3bNmbPns2MGTP42c9+RlVVVaufQ1EUzhmQxR0/HMHAXqmB6x4o4/FXC9h1qPxMqyAILMSTXFBCdF5hCxSFhYU8+uijLF26lLfeeotly5axd+/eJuc88MADLFiwgBUrVtCnTx+eeeaZNj9fl0Q710/MY/YlfXHYLFS7/Tyfv4vl6w/g0+RD7kzV1MnAthCdVdgCxYYNGxg9ejSpqakkJCQwceJE8vPzm5xjGAa1tbUAuN1unE7nGT2noiiMysvml7OH0TsnGYBN2wv53ze3cqKs7oyu3dlpukmdpCMXolNSzDCN/D711FPU1dVx5513AvDqq69SUFDAfffdFzznm2++4eabbyYhIQGXy8Urr7xCWlpayM9xorQW/TQZTw3DZM2mg6z4eD+GYWK1qFxzRX8uOScXRVaStYmiQGaKC7vNEumiCCE6kDVcFzYMo8kHsmmaTW57PB7uuecennvuOYYPH86zzz7L7373OxYtWhTyc1RU1OHXTj8j5/y8LLqluVi2dg9l1V5eXrOLgt1FzL6kHy5H2KreLtLTEykrq410MZpx13pJSbSf8X7bWVnJFBdXt0+hokw81w2kfrFMVRUyMpJa/7gwlAWAnJwciouLg7eLi4vJzs4O3t69ezcOh4Phw4cD8OMf/5jPPvus3cvRMzuJ22cPY8TZGQBs/7acJ97YwtGS6PsQjgVur4avheAshIg/YQsUF154IRs3bqSsrAy3282aNWsYN25c8Hjv3r05ceIE+/fvB2Dt2rUMGzYsLGVx2q386LKzmTWuL1aLQlm1l/97ayubthfKmotWMk2olX22hehUwtb/0rVrV+68807mzZuH3+/nmmuuYfjw4cyfP58FCxYwbNgwHnzwQe644w5M0yQjI4O//vWv4SoOiqLwg4HZ5GYlsvS9PZRWeli+/gCHCqu56uK+2KwxvaSkQ3n9Oj7NwGaR/zMhOoOwDWZ3hJ37ilscozgdj0/jjY/2s3V/GQDdMxKYO34A6V3ObNZVe4rWMYoGTruFtGRHm8cq4rkfOJ7rBlK/WBZ1YxTRzGm3cu0V/Zl8QS8UBY6V1vHEm1vZe6Qy0kWLGV6/jlc2OBKiU+iUgQICXVEXj+jOzVMGkeC04vZqPPfODjZsPS7jFiEwzUAeKBmrECL+ddpA0aBfbgq3zxpG94wEDBNWbTjIGx/tR5NEeN9LWhVCdA6dPlAApCY5+OnMIQzrmw7Al7uKeWbVDmrckliwJaYJdTIDSoi4J4Gint1qYc4V/ZlwXk8U4GBhNU++tZXCckn90RKvJq0KIeKdBIpGFEXh0nNymTt+ADaLSnm1l6eWb2PPkYpIFy1qyViFEPFPAsUpDOmTzk9nDCY5wYbHp/P8Ozv5YmdRpIsVtWSsQoj4JoHiNHKzkvjFVUPpVj/I/cZH+3n3i8MyI+oUpFUhRHyTQNGClCQHP50+hP49UgD44KujvLZun8yIOoWG1dpCiPgjgeJ7OOwW5k3K4wd5WQB8vaeEF/J34fXJZkiNSQ4oIeKXBIoQWFSVq8f15cof9ABg79FKnl61XabPfoeMVQgRnyRQhEhRFC4/twdXj+uLosDRklqeWr6NsipPpIsWNRpaFciWqULEFQkUrXTewGx+Mn4AVotCaZWHp5Zvk21WG5FWhRDxRwJFGww6K51bpg7GabdQ7fazeOU2DhXGZ7bJtqh2+yJdBCFEO5JA0Ua9c5KZP30wSS4bbq/OP9/eIdln62maSZ1Xi3QxhBDtRALFGeiWkcjPZgwhNcmOTzN4Pn8n278ti3SxokKt248ha06EiAsSKM5QRoqTn80cSlaqC90wWfrubjbvLYl0sSJON0zcMoVYiLgggaIdpCTamT99cHAV9yvv7+XLXZLyw+31y/wnIeKABIp2kuSyceu0wfTMTsIEXv9wP59uOxHpYkWUrpn4ZQaUEDFPAkU7cjms3DxlEGd1SwZgxSffsmHr8QiXKnJMwO2TQW0hYp0EinbmsFu4cfJAzs4N5IdateEg6ws6b7Dw+nV0QzqghIhlEijCwG61cP3EvGAywX9/epCPNh+LcKkiwzBMvH4Z1BYilkmgCBObVeUnE/LI65kKQP6mQ502WLi9khNLiFgmgSKMbFaV6yYMYGCvVCAQLD4u6HzBQtNNSUEuRAyTQBFmVovK3PEDgi2Ldz491OnGLEwTPDKoLUTMkkDRARqCxYD6YPHvTw92utlQHp8uK7WFiFESKDqIzapy3fgBwQHuVRsOsml7YYRL1XEMw8QjK7WFiEkSKDpQwwB3v9wuACxff6BTreCuk5XaQsQkCRQdzGZVuX5CHmflBBblvfHhfr7pJLmhNM3EL1NlhYg5EigiwG6zcMOkgcF0H699sJdtBzpH1tk6jyb7agsRYyRQREjDCu7umYkYJry8dg+7D1dEulhh59V0NF2mygoRSyRQRJDLYeWmKQPJTgukKH9xzS72H6uKdLHCyjSR9ONCxBgJFBGW6LRx89RBpHdxoOkmL6zeyZGimkgXK6w8Hg1dWhVCxIyQAsW//vUvamri+8Mrkrok2Lll6mBSEu34/AbPvrOTY8Xx+/+tGbJVqhCxJKRAsWvXLiZOnMg999zDli1bQr74ypUrmTJlChMmTGDJkiXNju/fv5/rr7+eGTNmcMstt1BZ2Xn3nE5LdnDL1EEkumy4vRqPLfuGsipPpIsVNm4Z1BYiZoQUKO6//35Wr17N0KFD+fOf/8zs2bN57bXX8Hq9p31MYWEhjz76KEuXLuWtt95i2bJl7N27N3jcNE1+/vOfM3/+fFasWMGgQYNYtGjRmdcohmWmurhp8kCcdguVNV7++fYOqup8kS5WWPh1Q/I/CREjQh6jSEpKYtKkSUybNo2KigqWLl3KpEmTeP/99095/oYNGxg9ejSpqakkJCQwceJE8vPzg8e3bdtGQkIC48aNA+C2227juuuuO8PqxL7umYncMGkgNqtKWbWXZ9/egTtOu2nitV5CxJuQAsXGjRu54447mDRpEvv37+eJJ57gjTfe4Pnnn+ePf/zjKR9TVFREVlZW8HZ2djaFhSdTVhw6dIjMzEzuvvturr76au69914SEhLOsDrxoXdOMrfNGo5FVSgsd/N8/k58cbhQTfI/CREbrKGc9Oc//5m5c+dy3333kZycHLy/V69e/OhHPzrlYwzDQGnUCW2aZpPbmqbx2Wef8eKLLzJs2DD+53/+h4ceeoiHHnoo5MKnpibE7e5p6emJ3DhtMP9csY1DhTW8+uF+fj57OFZL/ExUS01NICHBRnKiI9JFaXdZWcnff1IMk/p1LiEFihUrVpCfn09ycjLFxcW8/fbbzJs3D1VVWbBgwSkfk5OTwxdffBG8XVxcTHZ2dvB2VlYWvXv3ZtiwYQBMmzbttNc6nYqKOvxx2s+dnp5Iv5xkpo89ixXrv2Xb/lIWvVHAjy4/GzUORoHT0xMpK6ulqkohs4sr0sVpV1lZyRQXV0e6GGEj9YtdqqqQkZHU+seFctJ9993HunXr6p9I5csvv+Svf/1ri4+58MIL2bhxI2VlZbjdbtasWRMcjwA455xzKCsrY+fOnQC8//77DBkypNUViHejB+dw5Q96AFCwr5S3NxzEjKPuGl0z8fnjM9gLES9CalF8/fXXrFq1CoCMjAwee+wxZs6c2eJjunbtyp133sm8efPw+/1cc801DB8+nPnz57NgwQKGDRvGE088wcKFC3G73eTk5PDwww+feY3i0GXn5FLr1ti47QQbt50gyWXjsnNzI12sdmECbp+G3WaPdFGEEKcRUqDw+/34fD7s9sCbWdNCm60yffp0pk+f3uS+xYsXB38eMWIEr732Wqhl7bQURWHqhb2p9fgp2FfKu18cJsll5bxBXSNdtHbh9QcGteOhS02IeBRSoLj00ku55ZZbmDlzJoqisGrVKi655JJwl000oioK11zaD7dXY8+RSt5af4AEp40hfdIjXbQzZhgmXr+Oyx7Sr6MQooOFNEbx29/+lvHjx7N27VrWrVvH+PHj+fWvfx3usonvaNhStUdWIqYJy97fw4Hj8ZFE0O3VZaW2EFFKMWN4ZHTnvuK4nvVUVlZ7ymM1bj9PrdhGaaUHp93CT2cMISc9ttagfLd+igIZXZxxMf03nmfNgNQvloV11tN7773H5ZdfzqhRozj33HODXyIyklw2bpo8kGSXDY9P57l3dlJRc/p0KrHANMETh4sKhYgHIXUK/9d//Rd33XUXgwcPbrJoTkROehcnN04ZyKIV26mq9fHsv3fysxlDSHDGbj+/x6OR6LBJF5QQUSakFkWXLl2YMGECPXr0IDc3N/glIqtbRiLXTxyARVUornDzwuqdMd0Vpxkmfk1aFUJEm5ACxYgRI/jwww/DXRbRBn27p/DDy85GAQ4V1vDy2j0xndZEBrWFiD4h9VN8+OGHvPjii9hsNmw2WzBv01dffRXu8okQDO+XQXWdj7c3HmTHwXJWfnKAmWP7xGQ3YWBPbROLGntlFyJehRQonnvuuTAXQ5ypi4Z1o7rOx0ebj/PZjiK6JNq5/NwekS5WqxmGiVfTSZA1FUJEjZC6nnJzc9myZQuvvPIK6enpfP311zJGEYUmnN+LkWdnAvDeF0f4cldRhEvUNm6PRiC5hxAiGoQUKBYtWsRLL71Efn4+Ho+Hf/zjHzzxxBPhLptoJVVRmHVJX87OTQHgzY/2s+tQeYRL1XqabuDTJFAIES1CChRvv/02ixcvxuVykZaWxiuvvBJMEiiiS2D1dn+6ZSRgmLD0vT0cKa6JdLFaxTQD+Z+EENEhpEBhtVqDCQEhMF3WapU+5GjltFu5YfJAUpPs+DWD5/N3UVrliXSxWsXt1TCl+0mIqBBSoOjWrRvr1q1DURR8Ph9PPvmkjFFEuS4Jdm6cMgiXw0qt289z/95Jjdsf6WKFLJAoMHbXhAgRT0IKFH/4wx949tln2bVrFyNHjuSjjz7iD3/4Q7jLJs5QdqqLeRPzsFoUSqs8vBBje2+7PZqsqRAiCoTUf9S1a1eef/553G43uq6TlNT6pFIiMnrnJDPniv4seXc3R4preXntHq6bkBcT6xR8mo6mG1jU2E8UKEQsCylQPPvss6e8/6abbmrXwojwGHxWOtMvPIsVn3zLzkMVrFh/gKsujv4FeaYJbp9OklMChRCRFFKg2L17d/Bnn8/H559/zpgxY8JWKNH+Rg/JobLWx4ffHOPznUWkJMXGgjyPRyPJaYt0MYTo1EIKFA8++GCT24WFhdxzzz1hKZAInwnn9aSq1sfXe0p474sjpCTaGZWXHelitUgzTDx+HafNEumiCNFptalN37VrV44ePdreZRFhpigKV49ruiBv9+GKyBYqBDKoLURktXqMwjRNtm7dSkZGRtgKJcKnYUHe4pXbOV5ax9J3d3Pr9MH0yIreCQo+TcevGXGx+50QsSikd97u3buDX3v27KFbt2488sgj4S6bCJPGC/J8MbAgr2FQWwgRGbJndpRqac/s9lJU4eap5dtwezUyujj52cwhJLk6ZuC4tfVTVYWsVCcK0d8HFc97LoPUL5a1dc/skLqerr/++hanUr7wwgutfmIReQ0L8p55e3twQd6t0wZjj8KBY8Mw8fh0XJJ+XIgOF1LX09ChQ3E6ncybN49bbrmFzMxMUlNTue6667juuuvCXUYRRg0L8hQFjhTX8tJ7e9CN6Gylye53QkRGSH+effXVVyxduhSLJfCX5sUXX8yPfvQjJk6cGNbCiY4x+Kx0ZlzUh+XrD7DrcAVvfXyAWeP6Rt2CPL+m49cNrLJSW4gOFdI7rqysDK/XG7xdW1uLxxO9g5+i9S4Y3JXLzgkkevxyVzHvfXEkwiVqLpB+PDpbO0LEs5BaFNOmTePHP/4x48ePxzRN3nnnHebNmxfusokOduUPelBV5+PLXcV88PVRkhNtjB6cE+liNeH2+kl0yDiFEB0ppHfcr371KwYPHsynn36Kw+HgL3/5C+eff364yyY6mKIoXHVxH2rq/Ow6XMHK9d+S5LIztE96pIsWpOsmft3AJmsqhOgwIb/bunbtSv/+/bnjjjuw2ST3TryyqCrXXtmfntlJmMAr7+9h/7GqSBcryDTBI2sqhOhQIQWK119/nd///vc8/fTTVFdX84tf/IJXXnkl3GUTEWK3WZg3KY/MFCeabvLiml0cLw3vmo7W8Ho1Ynf1jxCxJ6RA8eKLL7Js2TKSkpLIyMjgjTfe4Pnnnw932UQEJTpt3DRlEMkJNjw+nefe2Ul5dXRMYNAMM24XWgoRjUIKFKqqNtmsqFu3bsGpsiJ+pSU7uHHyQJx2C9V1fp6Nou1UPX4t0kUQotMIKVCkpqayY8eO4Lz6FStWkJKSEtaCiejQLSOR6+u3Uy2p9PB8/k68UTBG4PHpGNL/JESHCClQ3H333fzmN79h3759jB07lscee4yFCxeGu2wiSvTp1iW4evtocS1L3t2Npke268cwTLwxtP+3ELEspEDh8XhYvnw5b775Jv/85z/Jz88nLy/vex+3cuVKpkyZwoQJE1iyZMlpz1u3bh2XX3556KUWHW7wWelcfXFfAPYereSVD/ZiGJH9i97t1QFpVQgRbiEFiv/3//4fFouFfv36MWDAgJCmxxYWFvLoo4+ydOlS3nrrLZYtW8bevXubnVdSUsLf/va31pdcdLgfDMxm4vk9Adi6v4wVnxwgksmHAyk9JFAIEW4hBYq8vDxWrlzJsWPHqKioCH61ZMOGDYwePZrU1FQSEhKYOHEi+fn5zc5buHAht99+e5sKLzreuBHdGTu8GwCf7SiKaKoP0wS3Vwa1hQi3kFZmr127ttmHvKIo7Nix47SPKSoqIisrK3g7OzubgoKCJue88MILDB48mBEjRrSmzEGpqQnoEe7+CKf09MRIF+GUrps8CM2AT7ce54Ovj5KRlsCV5/dq9XXao36KAmlpCVG3+11WVnKkixBWUr/OJaRAsWXLllZf2DCMJtlHTdNscnv37t2sWbOG5557jhMnTrT6+gAVFXVxO5++IzYuOhNTR/eistrDjoPlvPb+HkxdZ1RedsiPb8/6aV5/VO1TEc8b34DUL5a1deOiFv8M+8Mf/hD8uaysrFUXzsnJobi4OHi7uLiY7OyTHyT5+fkUFxcze/ZsfvrTn1JUVMTcuXNb9Rwiciyqwpwr+tOnWxcA3vhoP1sPtO53pL3UeTRkUFuI8GkxUGzdujX48y233NKqC1944YVs3LiRsrIy3G43a9asYdy4ccHjCxYsYPXq1SxfvpxFixaRnZ3N0qVLW1l8EUk2q8q8iXnkZiVimrBs7R72HKno8HJouoFPk0AhRLi0GCgaz2hp7eyWrl27cueddzJv3jyuuuoqpk2bxvDhw5k/f36burJEdHLYLdw4eSDZaS50w+TFNbv59kTHJhE0Tajz+mX3OyHCJOSO3bbsdjZ9+nSmT5/e5L7Fixc3O69Hjx68//77rb6+iA4NeaEWr9hGWbWX59/Zxa3TBpGb1fq+0Lby+nR8mqQfFyIcWnxXGYZBZWUlFRUV6Loe/DmU6bGic0lJtHPz1EF0SbTj9es8+++dFJbVddjzmybUuqVVIUQ4KGYLfUoDBw5EUZRTdjt93/TYjrBzX7HMeooyRRVuFq/YRq1HI9llY/6MwWSmuJqdF476KQqkJTuxWyPbqojnWTMg9YtlbZ311GLX086dO9tcINE5Zae6uGnKIJ5etZ1qt59nVu1g/vTBpHdxhv25TRPq3H4cXRyyX4UQ7Ug6dEW7656ZyE1TBuKwWais9fHM2zuoqPF2yHN7NR2vPz5bmUJEigQKERY9s5O5YXIeNqtKebWXZ97eQVWdL+zP2zBWIesqhGg/EihE2JyV04V59XtZlFZ6eGbVdqo7IFj4NF3WVQjRjqIn74GIS/1yU/jJhDz+tXoXxRUennl7B7dOG0x6G6+361A5H28+Rnm1l7RkBxeP6E5er7Qm58hYRXwo2FdC/qZDlFR6yExxMumCXgzvlxnpYnVK0qIQYTegZyo/mTAAi6pQVO4ODHS3oWWx61A5Kz45QJXbj9NhpcrtZ8UnB9h1qLzZuTJWEdsK9pWw5N3dVNT6SHBaqaj1seTd3RTsK4l00TolCRSiQ+T1SuO68SeDxaMvfdXq/bc/3nwMi0XFbrWgKAp2qwWLReXjzceanWuaUOeRdRWxKn/TISwWFYct8Fo7bIHXOn/ToUgXrVOSQCE6zMDeacy9sj8WVeFYcW2rWxbl1d5mK69tlsBg+al4/Tp+GauISSWVnmbrYexWlZJKT4RK1LlJoBAdatBZ6cwd37gbakfIwSIt2YH/O3t1+3WDtGTHKc83TaiVVkVMykxx4vvOYlqfZpCZEv71OKI5CRSiww3qncZts4ZjURWKK9wsXrmdytrvDxYXj+iOrhv4NB3TNPFpOrpucPGI7qd9jMenSasiBk26oBe6buD1B15rrz/wWk+6oPUbZIkzJ4FCRMSwszP5yYQBWC0KJZUeFq/YdtoupAZ5vdKYcVEfurhseLwaXVw2ZlzUp9msp8ZME2rcPmlVxJjh/TK5bvwAUhPt1Hk0UhPtXDd+gMx6ipAWcz1FO8n1FLsa6rfnSAUvrt6NXzdISbRz67TBZLRz90JH54CK51xBIPWLZWHZ4U6IcOvfI5UbpwzEblOprPWxaOU2CsvbN+tsoFXRuhlWQoiTJFCIiOvTrQu3TB2E026hus7P4hXbOVJc067P4fPrePx6u15TiM5CAoWICj2zk7l12mASnVbqvBrPrNrB/mPtu1NeTQekDxEiHkmgEFGje2YiP50xhJT6zY+ee2cHO0+x6rqtNN2k1qu12/WE6CwkUIiokpXq4qczhpCR4kTTTV5cvZuvdxe32/Vr3X50Iz4nQAgRLhIoRNRJS3bw0+mD6ZaRgGGavLpuHx8XNE/T0RaGYVJdJ2nIhWgNCRQiKiUn2Jk/fTB9uiUD8M6nh3jn04MY7TCb2+PT8UjCQCFCJoFCRC2n3cqNkwcx5KxAUvKPC47z2gf70PQz/5CvqfW1S9ARojOQQCGims2qcu2V/Tl/UDYA3+wt4bl3duLxndmgtGaY1LglD5QQoZBAIaKeqirMHNuHCef1BGD/sSqeWr7tjPfhdns12bNCiBBIoBAxQVEULj0nlx9e2g9VUSgsd/N/b23l6BkszDNNOmRrViFinQQKEVPOGZDFjVMG4rRbqKrzs2jldrZ/W9bm6/k1g1qPpPcQoiUSKETMOTs3hZ/NHBLYn0IzWLJmNx8XHKOt+S1r3ZqsrRCiBRIoREzqmpbAz68aSs/sJEwC02ff+HB/m2ZEGaasrRCiJRIoRMxKctm4ddpghvfLAODL3cWt3l61gcenU+eTpIFCnIoEChHTbFaVH19+NhPO64kCHCqs4X/fbNsgd3WtD59kmBWiGQkUIuY1zIj6yYQBwX0tnlqxja9amSPKNKGi1tdsX24hOjsJFCJuDDorndtmDiWjSyCh4Gvr9rF8/YFWjVsYhklltVcGt4VoRAKFiCs56Qn84uqh5PVKBWDT9kKeXrW9VYvzNMOkvMaLYUqwEAIkUIg45HJYuX5iHpefmwsExi3+8foWdh+uCPkammZSVuVFM2QmlBBhDRQrV65kypQpTJgwgSVLljQ7/t577zFz5kxmzJjBL37xCyorK8NZHNGJqIrClT/oyQ2T8nA5ArvmPf/OTt794jBGiB/+mm5SXuWRMQvR6YUtUBQWFvLoo4+ydOlS3nrrLZYtW8bevXuDx2tqavjTn/7EokWLWLFiBXl5efz9738PV3FEJ5XXK41fzh5Gj6xETOCDr47yzNvbqQyxK0o3TCqqvfgkJ5ToxMIWKDZs2MDo0aNJTU0lISGBiRMnkp+fHzzu9/u599576dq1KwB5eXkcP348XMURnVhqkoOfzhjChUNzADhwvJrHX9/CjhBTf+iGSXmNp34bVemKEp1P2AJFUVERWVlZwdvZ2dkUFhYGb6elpTF+/HgAPB4PixYt4sorrwxXcUQnZ7WoTLvwLK6fmEeCw4rbq/GvNbtZsf4APu37106YZmCdRVWd5IUSnY81XBc2DAOlUbJ/0zSb3G5QXV3Nf/zHfzBw4ECuvvrqVj1HamoCehwPNqanJ0a6CGEVifpdlJ7I4LMz+eeKbew5XMGn2wv5trCam6cPoVdOl5CuYVpUUhLtOBynf/tkZSW3V5GjktSvcwlboMjJyeGLL74I3i4uLiY7O7vJOUVFRdxyyy2MHj2au+++u9XPUVFRh1+Lz77j9PREyspqI12MsIl0/W6YmMfHBcd49/MjnCit46Hnv+CKUT0YN7I7FvX7dzMqLlFIdFlJcNr47tlZWckUF1eHp+BRQOoXu1RVISMjqfWPC0NZALjwwgvZuHEjZWVluN1u1qxZw7hx44LHdV3ntttuY/Lkydxzzz2nbG0IES6qqnDJyFx+fvVQslKdGKbJu18c5qnlWyksr/vexzckEiyXxXmiEwhbi6Jr167ceeedzJs3D7/fzzXXXMPw4cOZP38+CxYs4MSJE2zfvh1d11m9ejUAQ4cO5YEHHghXkYRoJjczkf+YNYzVnx1m49YTHCmu5Yk3tnDlqJ5cNLzb97YufH6d0iqDLol2XHYLsg23iEeK2dYk/lFg575i6XqKUdFYvw+/OcraL4+g6YG3REYXJxcO7cq2A2WUV3tJS3Zw8Yju5PVKa/ZYRYEEp41kl+2Mui4K9pWQv+kQJZUeMlOcTLqgF8P7ZZ5RvdrLivX7WfP5ETx+HafNwoTzejBjbN9IF6vdSdfTKR4XhrIIEXN2HSrnsx2FpCTZSagfpC6t8rByw0GOl9XhsFuocvtZ8ckBdh0qb/Z404Rat5+SSjdVNV6MNvz9VbCvhCXv7qai1keC00pFrY8l7+6mYF/JGdfvTK1Yv58VG77F69exquD166zY8C0r1u+PdNFEB5BAIQTw8eZjWCwqDpuV1GQHGSnO4DG3V6e4woNhmFgsKh9vPnba62i6SbXbT0mlhxqPH80wCHX4LX/TofoyWFAUBYfNgsWikr/p0JlW74yt+fwICgoWVUFR1MB3FNZ8fiTSRRMdIGxjFELEkvJqL85G010dNgsKJ5fX6UYg95PTbsEfwgZHhmFSU+en1u3HYbWQ4LRis1paDBollR4SnE3fknarSkmlpw01al8en9ZsvEZVAveL+CctCiEgsP/2d3I6WS0KNotCVqoLuzXwVvH4dKrcfj785mhI6ctNEzx+nbJqL6WVbqrdfnynGVfLTHE2O+bTDDIbtW4ixWm38t0lS4YZuF/EPwkUQgAXj+iOrhv4NB3TNPFpOnabBbvdiolJehcHSQk2FCXw4b/6s8M89moBOw6WE+p8EM0wqXX7Ka/2UFrlwedv2i016YJe6LqB1x8og9evo+sGky7oFaZah27CeT0wMdENE9M0At8xmXBej0gXTXQA+XNACALJA2cQGKtomOE0dXRvaHRfdoqTC8b05lBhDRu3naC0ysO/Vu/i7NwUJo/uRbeM0Faamyb4NYPyGg8Ou4Uklw2bRQ3OborGWU8Ns5s6w6wn0ZxMj41S0Th9tD3Fev0Ky+t4e8NB9h4NpMZXgHMHZHHleT3p0zOtVXVTFQW7TcXltGK3qijN1npHl3iePgrxXb+2To+VFoUQbdA1LYGbpgxk1+EK3vn0EMUVbr7cXUzBvlIu+0FPzs/LajYwfTqGaeLx6Xh8Oqqq4HJYSXBYsFpUWcAnooIECiHaSFEUBvZKo3+PVL7YWcR7Xx6h1u1nzaaDfPT1ES4e3p0Lh+bgsFtCvqZRP47h9mg47BZcDmtwIF2ISJFAIcQZsqgKFwzuysj+mXyy5TjrtxzH49V594vDrN9ynIuHd2PMkFYGDNPE7dUC01ItCgkOGw6bKq0MEREyRhGlYr0P//vEc/3sLjvLP9jDp9sKg1NuExxWLhrWjdFDuuJqIT15S1RFwWpVcdhUbFZLYDyjfhZWR4rnPnyI7/rJGIUQUSLJZWPy6N6MHd6NjwuOs2lbIXVejXe/OMxHm48xekhXLhyaQ3KCvVXXNUwTn1/H59dRFD+KomC3qNjsFmwWFZtVQVUUaXGIdieBQogwSU6wM2V0by4e3o31BcfZtKMQr1/nw2+O8cmW45zTP4uLhncjO9XV6mubZmAzMI+h4/HrKICingwcdmugmyq650+JWCGBQogwS06wM3l0by4ZmcvGbSfYsPUEbq/G5zuL+HxnEXm9UrloaDf65XZp874sJmAajQKHAhZFwW63YLdZsKoKVosMiou2kUAhRAdJcFq5YlQPLh7ejS93FfPJluOUVXvZdaiCXYcqyEp1MWZoV87pn4XDFvrA96mYJmimiebRqPNowcBhs1qw21VsVhWrKoFDhEYChRAdzG6zMGZoDhcM7sr2b8vYsPUE356oprjCzYr137J602FG9s/k/EHZIa/2/j7BwOHTcPtoEjisVgWLJZAR1mqRcQ7RnAQKISJEVRWG9s1gaN8MjpXUsnHbCTbvLcHr19m0vZBN2wvpmZ3EqLwshvfLaNcEfI0DB77AfYoSWBtitajYbQ2tDgWLtDw6PZkeG6XiefooxHf9zqRudR6Nr/cU89mOQoorTqYXt1lUhvRJ59wBWfTt3gX1e7ZobS+KElgnYrdasNYPkGdnJVNRXhuRqbkdQabHNhfTLYrAX0D1Nxr9wsbh767oJBKcgfUWFw7N4cDxar7cVcTW/WX4dYNv9pbwzd4SuiTYGHF2JiPOzqRbRkKbB8BDYZqBzZg0XQNvIKeVYrNQUeFGVcGiqlgsgdaHRVWCX/EYQDqzmG5RFBVXBfc3PmUtzKY/mI3eT83ON83g6d95WKPbJx9kNr8LTDAD/3z30vWHm0ezxu0hs9ETpqUmUF5R2/T5m92ov3mq+xqXpyXNi9shpEUROo9PY8u+Ur7aXcLBwqZ/6WamOBnWL4PhfTPITnOFNWg0OF39GrquVAWsFhVrQ/BQVFQ18NesqihR3xKRFkVzMR0oSktrML67m0oMCOW9nJmZTElJ635ZzSaBzQzOtW/yCtdv22bSPMiBCUrzY6ZpYpiBPESmYWIYJpppYBgnjze+TiiviASKtimr8gRaFntKmu18l5niZEifdIb0SSc3MzFsQaMt9WsIIgpgsdSPfVgU1CgMJBIompNAEaWi+Ze14fPHCAai+qBEoze4CaZyMkiZDS22+vNSUxIoK6/FNALnGKaJUR98TE4R4GJIRwRB0zQ5XlpHwb5Stuwvpbza2+R4l0Q7A3ulMrB3Gv26p2Brx8SC4ajfdwOJzaJisSpYIxBEovm9d6Y65RiFiIyGN6pC4M0dUhPpO1KTHfg9vuDtxpcwTBNdNzEM0EwDv8/Ar+v1u6udYeHjhKIodM9MpHtmIhPP78mx0jq2HShj24FSiis8VNX6+GxHEZ/tKMJmUenTvQsDeqaS1zOVjCjYWvW7Glq/ALoRSFPS4FRBxGoNzMZSlcCHX0OXV8O1RPuSQCGiQuM3t0JgPj8WsKNCfUok3QhswenXDPx+A79uYMRwy6O9KIpCbmYiuZmJTDivJ0UVbnYeLGfnwXIOFlbj1w12H65g9+EKVhHYH/zs3BTO7pFC3+5dSHTaIl2FFn1vEEEBBVQVrEpgcF1VCf71odBo/KT+QQ1BRalvpTRcx6zvYm18/ZbK1VlIoBAxw6KqWFSwWy0oLurHTQx0IxBEGo+jNJtA0OhNfeoJBE3f9SY0yZMU6B4jOEEgmj8kslNdZKe6GDeiO7UeP3uPVAYDRa1Ho7zaG0wfApCTnkCf7l3o060LZ+Ukk+SK7sDRmNloAolhgIYOjQLJ6TQEgOBugo1ebE1VqKhwNzq50WGlIfDUt2YUBUWtDzaNL3y6521yjUaBSyUwRNgocEGguy046hjB3zkJFCImmYFx92DwgJb74NtjXNcwqG/BmPh1A79m4NN0dD16A0ei8+RUWqN+XGPfkUr2Hq3k2xOBWYMnyuo4UVbHxq0ngMCg+Fk5yfTqGvjKTHXWf2DFj5b+kDBN0L937PP7g1GolOA/J1tHjY8pKqhKYAKAUj/9WFUVVJT6Yw1dc0qzCze0poJ3tfFllEAhOoX2+CBvSHsBgdXLLnvgPr9m4PbpeL0aWhRPrlAbdVGNG9kdv2ZwuKia/ceq2H+8iiNFNWi6SUmlh5JKD1/sKgbAabfQMzuJ3KwkemYlMrQdV4iLplPcTzW9PtAE1vHRXLBldJoIEGzB1AcTh9VCenrryyivuBBnwDQDawaSXSrJLhsev47dakFVFIxobWbUs1lV+nZPoW/3FAA03eBocS3fnqji2xPVHC6soc6r4fHp7DlSyZ4jlYEHrtlNcoItOJjeLSORbhkJpCU74q7lEe0azyoM4ew2L5qSQCFEO3LaLGSludB9Prx+A003MAyzfiA+uqf9Wi0qvXOS6Z2TzCUEylpa5eFQYQ1Himo4UlzD8dI6dMOkus4fzHrbwGGz0DXdRde0BHLSE8hOD4yVJLlsHbIQUISPBAohwsCiqiQ4AuMmjded6PWBQzMa1pUEvmt+A1/9sWihKAqZKS4yU1ycOyALCLQ63JrJzv0lHC2p5VhJLYVlbvy6gdevc6iwhkOFNU2u43JYyEp1Bb5SXGSmOslMdZGe7JA9MmKEBAohwqzxuhOrRa2f9tuU4goMoGqaGVgzogeCiR5lU4CtFpXeWYkkOyycV3+fYZiUVHk4UVrLiTI3hWV1FJbXUV7lxQTc3lMHEAVISbKTkeIkPdlJehcH6V2cpCU7SEt2kOCwSkskSkigECIKmGZgsNluU7DbTrZEAkn5Aq0Nn1ePulYHBBa8NUzJHd7v5P0+TaekwkNRhZuicjclFe76gXI3mh5oSVXU+Kio8bGPqmbXtVtVUpMdpCbZSUl0kJJkJzXJQZdEOymJdrok2s94gycRGgkUQkSphlaE1RJI753osGKYJn6/iVfT8fo1DJ2oHTS3Wy3BAe/GDNOkqtZHaaWH0ioPpZUeyqq9lFd5KK3y4q1fB+HTDIrKA0HmdBw2C8kJNrok2klOsJHsCnxPctlIavjuspHgtGHpoNTs8UgChRAxIrB25GSrI9llQzcC6U78hoGuGcHuqob8WdFIVRRSkxykJjnol5vS5Jhpmnh8OuXVXipqvJRXe6ms8VFRE7hdVeujus4fnLzj9et4K/VmCRJPxeWwkOi0kei0keC0kui0kuC0kuCw4XJaSXBYcTks1GkGPo8fl8OK3apK9xcSKISIaQ37P9hRwXGyu0o3zMCYh26gaYF0J4FcWWbEUsuHQlEUXA4rLoe1WUukQWDWlY/qOh+VtX6qa31U1fmoqfNT7fZTXeejxu2n1u2ncS+d26vj9oYWVBqoCjjtVpwOS+C73RL8ctis9d8D+5A7bJbgl91mwW4L3Ge3Bn621OekikVhDRQrV67kySefRNM0brjhBq677romx3fs2ME999xDbW0tP/jBD/jzn/+M1SqxS4i2Mk0o2FdC/qZDlFR6yEp1Mnl0b4b1zeBf+TvZvK8MiyWwYHBw70B22Y8LjlNa6SHZZWP00Bz690hl16FyPt58jPJqL2nJDi4e0Z2jxTWsLziBV9NxWC2MHZ7D5aN6nrIcp3p8Xq+0097fmmt8+PWRQDn8Og5boBwTz+/V7PGGabJlXwmfFBynosZHgtNKn+5d8Pp09h+rwu3TsaoKCU4bhmlS59GC3V4nrwF1Xo06rwZ4mz1Ha6gK2OqDRkPwsFnrv+r377BZvnO74ef629b61y7wFUiMaKtPkHjyvsAe6A3HG1Zyn4mwpRkvLCzk2muv5Y033sButzNnzhz++7//m7PPPjt4zrRp07j//vsZOXIkd999N0OHDmXu3LkhP4ekGY9d8Vy/SNatYF8JS97djcWiYreq+DQDXTfISHaw83Blk3NVVcFlV8lOT8RhUwNdVSacOyCLgn2lUJ93yOfXqazx4vEZmIYBCvi1QIvl0pHduOzcpsFi16FyVnxyIJjp1a8HyjBqQBZf7i5udv+Mi/o0Cxanu0bv7CQ27y9DIfDBa9S3jq44N7dZ0DrVNdweDRQFl8NyyjLohokr0cGxE5V4vDp7jlSwYeuJ4MrmhnUxPbKScNgteHx6oPur4btfx+ePvu2ZG7a0zclI5MnfXdHqx4ftz/cNGzYwevRoUlNTAZg4cSL5+fncfvvtABw9ehSPx8PIkSMBmDVrFo8//nirAkVH7RscKVK/2BWpum3aXkh2egJ268nZQD5Np7jcHdgBr9G5DX9ipSY5mpy7cWshSYk2nHYrqqKQ6LLh1wwSXQp2i4pqCWRfNYHSKh8Zqa6Te5KYJkeKaxnQKw2bRcUwApte+XWdI8W1nNU9JTBFuD6Zn8+vs+dwBcP7ZTbZjXHnoQq6ZyViVS3Be/2aTmG5m5x0V/D/tyHb667DlUy8oHeT/4ttB8qa/V+UVgYGxjNSXE3qvO1AGUP7ZmADUhLtmBmBbq/PdhTSo2tSs//PZKeNH1/R/5SvgQFo9XnAGtbHBLId6/i1+uzHul4/FdrA79fRDAO/Fugq9Osnx5v8WiAwBQJa/ZRpzUAzA+e09u/k5IS2JXwMW6AoKioiKysreDs7O5uCgoLTHs/KyqKwsLBVz5GWduo+zHjRlg1GYkk81y9Sdbv75tERed7G7pg76oyvsbBv1vef9D3u7df2a+TUB4ozuUY8CduySMMwmgzcmKbZ5Pb3HRdCCBEdwhYocnJyKC4uDt4uLi4mOzv7tMdLSkqaHBdCCBEdwhYoLrzwQjZu3EhZWRlut5s1a9Ywbty44PHc3FwcDgdffvklAMuXL29yXAghRHQI26wnCEyPfeqpp/D7/VxzzTXMnz+f+fPns2DBAoYNG8bOnTtZuHAhNTU1DBkyhAcffBC7/btZcIQQQkRSWAOFEEKI2Cc5foUQQrRIAoUQQogWSaAQQgjRIgkUQgghWhQTgeKxxx5jypQpTJ06lWeffRYIpAiZPn06EyZM4NFHH41wCdvH3/72N+666y4gvup3/fXXM3XqVGbOnMnMmTPZvHlzXNXv/fffZ9asWUyePJn7778fiJ/X79VXXw2+bjNnzmTUqFH85S9/iZv6LV++nKlTpzJ16lT+9re/AfHz2gEsWrSIiRMnMn36dJ588kmgjfUzo9ymTZvMOXPmmH6/33S73eZll11m7tixw7zkkkvMQ4cOmX6/37z55pvNdevWRbqoZ2TDhg3mBRdcYP7ud78z3W533NTPMAxz7Nixpt/vD94XT/U7dOiQOXbsWPP48eOmz+czr732WnPdunVxU7/Gdu/ebY4fP948duxYXNSvrq7OPO+888zS0lLT7/eb11xzjbl27dq4qJtpmuYnn3xiTps2zayurjY1TTN/9rOfmcuXL29T/aK+RXH++efzwgsvYLVaKS0tRdd1qqqq6N27Nz179sRqtTJ9+nTy8/MjXdQ2q6io4NFHH+W2224DoKCgIG7qt3//fgBuvvlmZsyYwYsvvhhX9Xv33XeZMmUKOTk52Gw2Hn30UVwuV9zUr7E//elP3HnnnRw+fDgu6qfrOoZh4Ha70TQNTdNISkqKi7oBbN++nbFjx5KUlITFYuHiiy/m1VdfbVP9oj5QANhsNh5//HGmTp3KmDFjTplwsLUJBaPJH//4R+688066dOkCnDqhYqzWr6qqijFjxvDEE0/w3HPP8fLLL3Ps2LG4qd/BgwfRdZ3bbruNmTNnsnTp0rh6/Rps2LABj8fD5MmT46Z+SUlJ/OpXv2Ly5Mlccskl5Obmxk3dAIYMGcL69eupqKjA6/Xy/vvv89VXX7WpfjERKAAWLFjAxo0bOX78ON9++23cJBR89dVX6datG2PGjAneF08JE8855xwefvhhkpOTSU9P55prruHxxx+Pm/rpus7GjRv561//yrJlyygoKODw4cNxU78GL7/8MjfddBMQP7+fO3fu5PXXX+eDDz7g448/RlXVuPpsGTNmDLNmzeL666/n1ltvZdSoUWia1qb6Rf12cvv27cPn8zFo0CBcLhcTJkwgPz8fi+VkfvjvJhyMJf/+978pLi5m5syZVFZWUldXx9GjR+Omfl988QV+vz8YCE3TJDc3t8WEkbEkMzOTMWPGkJ6eDsCVV14ZV7+fAD6fj88//5yHHnoI+P6En7Fi/fr1jBkzhoyMDCCwJ84zzzwTN69dTU0NEyZMCAb4p59+mvPPP79Nr13UtyiOHDnCwoUL8fl8+Hw+1q5dy5w5czhw4ECw2b9q1aqYTSj47LPPsmrVKpYvX86CBQu4/PLLefrpp+OmftXV1Tz88MN4vV5qamp48803+fWvfx039bvssstYv349VVVV6LrOxx9/zKRJk+KmfgC7du3irLPOIiEhAYARI0bERf0GDhzIhg0bqKurwzRN3n///bipGwQ+O3/xi1+gaRrV1dW89tpr3HHHHW2qX9S3KC655BIKCgq46qqrsFgsTJgwgalTp5Kens4vf/lLvF4vl1xyCZMmTYp0UduNw+HgoYceiov6XXbZZWzevJmrrroKwzCYO3cu55xzTtzUb8SIEdx6663MnTsXv9/PRRddxLXXXkvfvn3jon4Ahw8fJicnJ3g7Xn4/x44dy/bt25k1axY2m41hw4bxy1/+kosuuijm6waBQDhhwgRmzJiBruvceOONjBo1qk2vnSQFFEII0aKo73oSQggRWRIohBBCtEgChRBCiBZJoBBCCNEiCRRCCCFaFPXTY4Vojfvvv5/PP/8cCCzWzM3Nxel0ArBs2bLgz9GmoKCA1157jb/85S+RLooQzUigEHFl4cKFwZ8vv/xyHnnkEYYNGxbBEoVm7969MZtTSMQ/CRSiU3j11Vd56aWXMAyD1NRU/vCHP9CvXz/uuusunE4nu3fvprS0lMsvv5zU1FQ++OADiouLuf/++xkzZgx33XUXDoeDnTt3UlpaykUXXcTChQux2Wzs27ePBx54gIqKCnRd5/rrr+eaa65h06ZNPPDAAyQkJFBbW8vrr7/Oww8/zObNm6mtrcU0Te6//366d+/O448/TnV1Nb///e+56qqruO+++1i1ahUAmzZtCt7++9//zjfffENRURF5eXk88sgjPPnkk6xZswbDMMjNzeXee++la9euEf4fF/FEAoWIe5999hlvvfUWS5YsweVysX79em6//XbeeecdIJCOecmSJVRUVDB27FgWLlzIyy+/zPPPP8/ixYuDeaoKCgp48cUXsdls3HzzzSxbtow5c+awYMECHn74YYYMGUJ1dTU//vGPOfvsswHYs2cP7733Hrm5uXz99dcUFRWxbNkyVFVl0aJFLF68mP/7v/9jwYIFrF69mgcffJBNmza1WJ+jR4+yatUqrFYrb731Frt37+bVV1/FarWybNkyFi5cyOLFi8P7nyo6FQkUIu6tW7eOgwcPMmfOnOB9VVVVVFRUAIE0IzabjaysLBISErj44osB6NWrV/AcgKuvvprExEQAZs6cydq1axk9ejSHDh3i7rvvDp7n8XjYvn07/fr1o1u3buTm5gKBTLopKSm8/PLLHD58mE2bNgWv1xojR47Eag28dT/44AO2bNnC7NmzAYL7KwjRniRQiLhnGAYzZ87kN7/5TfB2UVERKSkpANjt9ibnN3wIf1fjrKKmaaKqKrquk5yczPLly4PHSkpKSE5O5ptvvgkm0oNAwHrggQe46aabuOKKK+jbty8rVqxo9jyKotA4s47f729yvPE1DcMI5pqCQKbXysrKlv9DhGglmR4r4t7YsWN5++23KSoqAuCll17ihhtuaPV13nnnHXw+H16vlzfffJPLLruMPn364HQ6g4Hi+PHjTJs2ja1btzZ7/CeffMJll13G3LlzGTp0KO+99x66rgOBIKRpGgDp6ekcO3aM0tJSTNPk7bffbrFur732GjU1NUBgf/nf/va3ra6bEC2RFoWIe2PHjmX+/PncfPPNKIpCUlIS//jHP1q9IY3T6WTu3LlUVVUxceJEZs+ejaqq/O///i8PPPAATz/9NJqm8atf/YpRo0Y1G2uYM2cO//mf/8n06dPRNI2LLrooOAg9cuRInnjiCW6//Xb+8Y9/MGfOHGbPnk1WVhaXXnopW7ZsOWWZfvjDH1JYWMiPfvQjFEWhW7duwX0jhGgvkj1WiBDcdddd9O/fn1tuuSXSRRGiw0nXkxBCiBZJi0IIIUSLpEUhhBCiRRIohBBCtEgChRBCiBZJoBBCCNEiCRRCCCFaJIFCCCFEi/4/BIq0ELyZUSkAAAAASUVORK5CYII=\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, truncate=False)\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 (ipykernel)", "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.11" } }, "nbformat": 4, "nbformat_minor": 2 }