{ "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.9.12 (main, Apr 5 2022, 06:56:58) \n", "[GCC 7.5.0]\n", "uname_result(system='Linux', node='CyberBandit', release='5.15.0-83-generic', version='#92~20.04.1-Ubuntu SMP Mon Aug 21 14:00:49 UTC 2023', machine='x86_64')\n", "IPython 8.2.0\n", "IPython.core.release 8.2.0\n", "PIL 9.0.1\n", "PIL.Image 9.0.1\n", "PIL._version 9.0.1\n", "_csv 1.0\n", "_ctypes 1.1.0\n", "_curses b'2.2'\n", "decimal 1.70\n", "_pydev_bundle.fsnotify 0.1.5\n", "_pydevd_frame_eval.vendored.bytecode 0.13.0.dev\n", "argparse 1.1\n", "backcall 0.2.0\n", "bottleneck 1.3.4\n", "cffi 1.15.0\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", "debugpy 1.5.1\n", "decimal 1.70\n", "decorator 5.1.1\n", "defusedxml 0.7.1\n", "distutils 3.9.12\n", "executing 0.8.3\n", "executing.version 0.8.3\n", "http.server 0.6\n", "ipykernel 6.9.1\n", "ipykernel._version 6.9.1\n", "ipython_genutils 0.2.0\n", "ipython_genutils._version 0.2.0\n", "ipywidgets 7.6.5\n", "ipywidgets._version 7.6.5\n", "jedi 0.18.1\n", "joblib 1.1.0\n", "joblib.externals.cloudpickle 2.0.0\n", "joblib.externals.loky 3.0.0\n", "json 2.0.9\n", "jupyter_client 6.1.12\n", "jupyter_client._version 6.1.12\n", "jupyter_core 4.9.2\n", "jupyter_core.version 4.9.2\n", "kiwisolver 1.3.2\n", "logging 0.5.1.2\n", "matplotlib 3.5.1\n", "numexpr 2.8.1\n", "numpy 1.22.4\n", "numpy.core 1.22.4\n", "numpy.core._multiarray_umath 3.1\n", "numpy.lib 1.22.4\n", "numpy.linalg._umath_linalg 0.1.5\n", "numpy.version 1.22.4\n", "packaging 21.3\n", "packaging.__about__ 21.3\n", "pandas 1.4.3\n", "parso 0.8.3\n", "patsy 0.5.2\n", "patsy.version 0.5.2\n", "pexpect 4.8.0\n", "pickleshare 0.7.5\n", "pkg_resources._vendor.appdirs 1.4.3\n", "pkg_resources._vendor.more_itertools 8.12.0\n", "pkg_resources._vendor.packaging 21.3\n", "pkg_resources._vendor.packaging.__about__ 21.3\n", "pkg_resources._vendor.pyparsing 3.0.8\n", "pkg_resources._vendor.appdirs 1.4.3\n", "pkg_resources._vendor.more_itertools 8.12.0\n", "pkg_resources._vendor.packaging 21.3\n", "pkg_resources._vendor.pyparsing 3.0.8\n", "platform 1.0.8\n", "prompt_toolkit 3.0.20\n", "psutil 5.8.0\n", "ptyprocess 0.7.0\n", "pure_eval 0.2.2\n", "pure_eval.version 0.2.2\n", "pydevd 2.6.0\n", "pygments 2.11.2\n", "pyparsing 3.0.4\n", "pytz 2022.1\n", "re 2.2.1\n", "scipy 1.8.1\n", "scipy._lib._uarray 0.8.2+14.gaf53966.scipy\n", "scipy._lib.decorator 4.0.5\n", "scipy.integrate._dop b'$Revision: $'\n", "scipy.integrate._lsoda b'$Revision: $'\n", "scipy.integrate._ode $Id$\n", "scipy.integrate._odepack 1.9 \n", "scipy.integrate._quadpack 1.13 \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.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._minpack2 b'$Revision: $'\n", "scipy.optimize._slsqp 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", "setuptools 63.2.0\n", "distutils 3.9.12\n", "setuptools._vendor.more_itertools 8.8.0\n", "setuptools._vendor.ordered_set 3.1\n", "setuptools._vendor.packaging 21.3\n", "setuptools._vendor.packaging.__about__ 21.3\n", "setuptools._vendor.more_itertools 8.8.0\n", "setuptools._vendor.ordered_set 3.1\n", "setuptools._vendor.packaging 21.3\n", "setuptools.version 63.2.0\n", "six 1.16.0\n", "socketserver 0.4\n", "stack_data 0.2.0\n", "stack_data.version 0.2.0\n", "statsmodels 0.13.2\n", "statsmodels.__init__ 0.13.2\n", "statsmodels.api 0.13.2\n", "statsmodels.tools.web 0.13.2\n", "traitlets 5.1.1\n", "traitlets._version 5.1.1\n", "urllib.request 3.9\n", "wcwidth 0.2.5\n", "xmlrpc.client 3.9\n", "zlib 1.0\n", "zmq 22.3.0\n", "zmq.sugar 22.3.0\n", "zmq.sugar.version 22.3.0\n" ] } ], "source": [ "def print_imported_modules():\n", " import sys\n", " for name, val in sorted(sys.modules.items()):\n", " if(hasattr(val, '__version__')): \n", " print(val.__name__, val.__version__)\n", "# else:\n", "# print(val.__name__, \"(unknown version)\")\n", "def print_sys_info():\n", " import sys\n", " import platform\n", " print(sys.version)\n", " print(platform.uname())\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import seaborn as sns\n", "\n", "print_sys_info()\n", "print_imported_modules()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading and inspecting data\n", "Let's start by reading data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \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 = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")\n", "data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv\")\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX4UlEQVR4nO3df5RcZX3H8fdnN0tI2Agx2K3NjwKSolQikuVX8ccGqw2cI6kHWtEKlhYjLbGVtsdQ22Op1p6CVSyWGiOlGHo0CqEQbVqE2uVHNZJQYyBQcEsw2YQGWANkQ9jsZr/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: Wed, 27 Sep 2023 Deviance: 3.0144
Time: 23:55:55 Pearson chi2: 5.00
No. Iterations: 6 Pseudo R-squ. (CS): 0.04355
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: Wed, 27 Sep 2023 Deviance: 3.0144\n", "Time: 23:55:55 Pearson chi2: 5.00\n", "No. Iterations: 6 Pseudo R-squ. (CS): 0.04355\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", " family=sm.families.Binomial(sm.families.links.logit())).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\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: Wed, 27 Sep 2023 Deviance: 18.086
Time: 23:55:55 Pearson chi2: 30.0
No. Iterations: 6 Pseudo R-squ. (CS): 0.2344
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: Wed, 27 Sep 2023 Deviance: 18.086\n", "Time: 23:55:55 Pearson chi2: 30.0\n", "No. Iterations: 6 Pseudo R-squ. (CS): 0.2344\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": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", "# family=sm.families.Binomial(sm.families.links.logit),\n", " family=sm.families.Binomial(sm.families.links.logit()),\n", " var_weights=data['Count']).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\n", "The Goodness of fit (Deviance) indicated for this model is $G^2=18.086$ with 21 degrees of freedom (Df Residuals).\n", "\n", "**I have therefore managed to fully replicate the results of the Dalal *et al.* article**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predicting failure probability\n", "The temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZh0lEQVR4nO3df3RV5Z3v8ffXEEr4IYxoGSEodC4XapVfCUHL1QZbATsdxbtoFR2trjLUabG1s2RG1p17de7orOmKc6XXpVKuUqd1aVCHIjqsBttpaseOChQEkQlQSyWhdxC5/IgNmITv/WPvkx7CCWfncA455+HzWisr2Xs/ez/PN0c/7Dxn733M3RERkdJ3Tl8PQERE8kOBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiKyBbmYrzGyfmb3dw3Yzs/9tZrvMbIuZTc3/MEVEJJskZ+hPAXNOsf1aYFz8tRB4/PSHJSIivZU10N39VeDAKZpcD3zfI68Dw8zswnwNUEREkumXh2OMAvakLTfH637bvaGZLSQ6i6eioqJq9OjRve7swFHnWKdjuY216DioliIUSi2h1AFh1VJ+DgyvyO0tzB07dux39wsybctHoGf6HWd8noC7LweWA1RXV/uGDRty6rCxsZHa2tqc9i02qqU4hVJLKHWAakkxs9/0tC0fV7k0A+mn2pXA3jwcV0REeiEfgb4GuC2+2uVy4JC7nzTdIiIihZV1ysXMngVqgfPNrBm4DygHcPdlwFrg88Au4HfAHYUarIiI9CxroLv7/CzbHfh63kYkIiWjvb2d5uZmjh49WvC+hg4dyvbt2wvez5mQpJYBAwZQWVlJeXl54uPm401RETlLNTc3M2TIEMaMGYNZYa9BOXLkCEOGDCloH2dKtlrcnQ8++IDm5mbGjh2b+Li69V9Ecnb06FGGDx9e8DA/25gZw4cP7/VfPgp0ETktCvPCyOX3qkAXEQmE5tBFpKSVlZVx2WWXdS2vXr2aMWPG9N2A+pACXURKWkVFBZs3b864zd1xd8455+yYjDg7qhSRs8bu3bv55Cc/yde+9jWmTp3Knj17qKurY9q0aUycOJH77ruvq+2DDz7I+PHj+dznPsf8+fN56KGHAKitrSX1aJL9+/d3nfF3dnayePHirmN997vfBX5/K/+8efOYMGECt9xyC9EV3bB+/Xo+/elPM2nSJGpqajhy5AizZ88+4R+hGTNmsGXLltOuXWfoIpIXf/PSNt7Zezivx7xk5Lnc9yefOmWbtrY2Jk+eDMDYsWN5+OGHaWpq4nvf+x6PPfYY69atY+fOnbz55pu4O9dddx2vvvoqgwYNor6+nk2bNtHR0cHUqVOpqqo6ZV9PPvkkQ4cOZf369Rw7dowZM2Ywa9YsADZt2sS2bdsYOXIkM2bM4LXXXqOmpoYbb7yRlStXMm3aNA4fPkxFRQW33XYbTz31FEuXLmXHjh0cO3aMiRMnnvbvS4EuIiWt+5TL7t27ufjii7n88ssBWLduHevWrWPKlCkAtLa2snPnTo4cOcINN9zAwIEDAbjuuuuy9rVu3Tq2bNnCCy+8AMChQ4fYuXMn/fv3p6amhsrKSgAmT57M7t27GTp0KBdeeCHTpk0D4NxzzwXghhtuYMaMGdTV1bFixQpuv/32vPwuFOgikhfZzqTPpEGDBnX97O4sWbKEr371qye0Wbp0aY+XBvbr14/jx48DnHAtuLvzyCOPMHv27BPaNzY28rGPfaxruaysjI6ODtw9Yx8DBw7kmmuu4cUXX+S5554j1yfPdqc5dBEJ2uzZs1mxYgWtra0AtLS0sG/fPq666ip++MMf0tbWxpEjR3jppZe69hkzZgwbN24E6DobTx3r8ccfp729HYAdO3bw4Ycf9tj3hAkT2Lt3L+vXrweiO0Q7OjoAWLBgAd/4xjeYNm0a5513Xl5q1Rm6iARt1qxZbN++nSuuuAKAwYMH8/TTTzN16lRuvPFGJk+ezMUXX8yVV17Ztc8999zDl770JX7wgx9w9dVXd61fsGABu3fvZurUqbg7F1xwAatXr+6x7/79+7Ny5Uruuusu2traqKio4Mc//jEAVVVVnHvuudxxRx6fZ5i6rOdMf1VVVXmufvrTn+a8b7FRLcUplFoKXcc777xT0OOnO3z4cEGPf99993ldXV1B+0g5fPiwt7S0+Lhx47yzs7PHdpl+v8AG7yFXNeUiInKGPfPMM0yfPp0HH3wwr9fIa8pFRAS4//77z1hfN99880lv0uaDztBF5LS4Z/wIYTlNufxeFegikrMBAwbwwQcfKNTzzOPnoQ8YMKBX+2nKRURyVllZSXNzM++//37B+zp69GivA65YJakl9YlFvaFAF5GclZeX9+oTdU5HY2Nj192epa5QtWjKRUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCUSiQDezOWbWZGa7zOzeDNuHmtlLZvaWmW0zszvyP1QRETmVrIFuZmXAo8C1wCXAfDO7pFuzrwPvuPskoBb4BzPrn+exiojIKSQ5Q68Bdrn7u+7+EVAPXN+tjQNDzMyAwcABoCOvIxURkVOybJ/WbWbzgDnuviBevhWY7u6L0toMAdYAE4AhwI3u/s8ZjrUQWAgwYsSIqvr6+pwG3drayuDBg3Pat9ioluIUSi2h1AGqJWXmzJkb3b0607YkHxJtGdZ1/1dgNrAZuBr4I+AVM/u5ux8+YSf35cBygOrqaq+trU3Q/ckaGxvJdd9io1qKUyi1hFIHqJYkkky5NAOj05Yrgb3d2twBrPLILuDXRGfrIiJyhiQJ9PXAODMbG7/ReRPR9Eq694DPApjZCGA88G4+ByoiIqeWdcrF3TvMbBHQAJQBK9x9m5ndGW9fBvwt8JSZbSWaovkrd99fwHGLiEg3SebQcfe1wNpu65al/bwXmJXfoYmISG/oTlERkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEApEo0M1sjpk1mdkuM7u3hza1ZrbZzLaZ2c/yO0wREcmmX7YGZlYGPApcAzQD681sjbu/k9ZmGPAYMMfd3zOzjxdovCIi0oMkZ+g1wC53f9fdPwLqgeu7tbkZWOXu7wG4+778DlNERLIxdz91A7N5RGfeC+LlW4Hp7r4orc1SoBz4FDAE+I67fz/DsRYCCwFGjBhRVV9fn9OgW1tbGTx4cE77FhvVUpxCqSWUOkC1pMycOXOju1dn2pZ1ygWwDOu6/yvQD6gCPgtUAP9mZq+7+44TdnJfDiwHqK6u9tra2gTdn6yxsZFc9y02qqU4hVJLKHWAakkiSaA3A6PTliuBvRna7Hf3D4EPzexVYBKwAxEROSOSzKGvB8aZ2Vgz6w/cBKzp1uZF4Eoz62dmA4HpwPb8DlVERE4l6xm6u3eY2SKgASgDVrj7NjO7M96+zN23m9mPgC3AceAJd3+7kAMXEZETJZlywd3XAmu7rVvWbbkOqMvf0EREpDd0p6iISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEIlGgm9kcM2sys11mdu8p2k0zs04zm5e/IYqISBJZA93MyoBHgWuBS4D5ZnZJD+2+DTTke5AiIpJdkjP0GmCXu7/r7h8B9cD1GdrdBfwTsC+P4xMRkYTM3U/dIJo+mePuC+LlW4Hp7r4orc0o4BngauBJ4GV3fyHDsRYCCwFGjBhRVV9fn9OgW1tbGTx4cE77FhvVUpxCqSWUOkC1pMycOXOju1dn2tYvwf6WYV33fwWWAn/l7p1mmZrHO7kvB5YDVFdXe21tbYLuT9bY2Eiu+xYb1VKcQqkllDpAtSSRJNCbgdFpy5XA3m5tqoH6OMzPBz5vZh3uvjofgxQRkeySBPp6YJyZjQVagJuAm9MbuPvY1M9m9hTRlMvq/A1TRESyyRro7t5hZouIrl4pA1a4+zYzuzPevqzAYxQRkQSSnKHj7muBtd3WZQxyd7/99IclIiK9pTtFRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCkegqF5FCWb2phbqGJvYebGPksAoWzx7P3Cmj+npYkpBev+KiQJc+s3pTC0tWbaWtvROAloNtLFm1FUChUAL0+hUfTblIn6lraOoKg5S29k7qGpr6aETSG3r9io8CXfrM3oNtvVovxUWvX/FRoEufGTmsolfrpbjo9Ss+CnTpM4tnj6eivOyEdRXlZSyePb6PRiS9odev+OhNUekzqTfOdJVEadLrV3wU6NKn5k4ZpQAoYXr9ioumXEREAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEApEo0M1sjpk1mdkuM7s3w/ZbzGxL/PULM5uU/6GKiMipZA10MysDHgWuBS4B5pvZJd2a/Rr4jLtPBP4WWJ7vgYqIyKklOUOvAXa5+7vu/hFQD1yf3sDdf+Hu/y9efB2ozO8wRUQkG3P3UzcwmwfMcfcF8fKtwHR3X9RD+3uACan23bYtBBYCjBgxoqq+vj6nQbe2tjJ48OCc9i02qqU4hVJLKHWAakmZOXPmRnevzrStX4L9LcO6jP8KmNlM4CvAf8m03d2XE0/HVFdXe21tbYLuT9bY2Eiu+xYb1VKcQqkllDpAtSSRJNCbgdFpy5XA3u6NzGwi8ARwrbt/kJ/hiYhIUknm0NcD48xsrJn1B24C1qQ3MLOLgFXAre6+I//DFBGRbLKeobt7h5ktAhqAMmCFu28zszvj7cuA/wEMBx4zM4COnuZ4RESkMJJMueDua4G13dYtS/t5AXDSm6AiZ9rqTS3UNTSx92AbI4dVsHj2eICT1s2dMuqM9F2IfpL469VbefaNPdx9aTtfWbKW+dNH88Dcy/pkLHLmJAp0kVKwelMLS1Ztpa29E4CWg20sfv4tMGjv9K51S1ZtBchr2GbquxD9JPHXq7fy9OvvdS13unctK9TDplv/JRh1DU1dgZrSfty7wjylrb2TuoamgvddiH6SePaNPb1aL+FQoEsw9h5sK0jb0zlevvtJorOHe0t6Wi/hUKBLMEYOqyhI29M5Xr77SaLMMt060vN6CYcCXYKxePZ4KsrLTlhXfo5RXnZikFWUl3W9WVrIvgvRTxLzp4/u1XoJh94UlWCk3nzsi6tceuq7L65ySb3xmZozLzPTVS5nCQW6BGXulFEZQ/RMBGtPffeFB+ZexgNzL6OxsZFf3VLb18ORM0RTLiIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCD6JWlkZnOA7wBlwBPu/vfdtlu8/fPA74Db3f2XeR6rSLBWb2qhrqGJvQfbGDmsgsWzx/P8hvd47VcHutrM+KPz+GL1RSe1A05at+E3B3j2jT3cfWk7X1mylvnTR/PA3MsS9Tt3yqge1yfZP9V3pztlZr3qO1MtSfvN1O5skzXQzawMeBS4BmgG1pvZGnd/J63ZtcC4+Gs68Hj8XUSyWL2phSWrttLW3glAy8E27l65+aR2r/3qwAkB33KwjcUvvAUO7ce9a91frNzM8bT9Ot15+vX3AE4I1kz9Llm1lQ2/OcA/bWw5aT1wQmhm2v90+l78/Ftg0N75+1qS9pup3dkoyZRLDbDL3d9194+AeuD6bm2uB77vkdeBYWZ2YZ7HKhKkuoamrnDqrfZO7wrzlOM9tH32jT1Z+21r7+TZN/ZkXF/X0JR1/9Ppu/24d4V5b/vN1O5slGTKZRSQ/mo0c/LZd6Y2o4Dfpjcys4XAwnix1cxyfQXOB/bnuG+xUS3F6YzV0v8P/1NVoY79rd8domzg0K5l+/s/3phrv78FbMmunPfPte+0fs8H9ve0b/fxFbnT+e/r4p42JAl0y7DOc2iDuy8Hlifo89QDMtvg7tWne5xioFqKUyi1mNmGjkP7Sr4OCOc1gcLVkmTKpRkYnbZcCezNoY2IiBRQkkBfD4wzs7Fm1h+4CVjTrc0a4DaLXA4ccvffdj+QiIgUTtYpF3fvMLNFQAPRZYsr3H2bmd0Zb18GrCW6ZHEX0WWLdxRuyEAepm2KiGopTqHUEkodoFqyMveTprpFRKQE6U5REZFAKNBFRAJR9IFuZgPM7E0ze8vMtpnZ38TrzzOzV8xsZ/z9D/p6rEmYWZmZbTKzl+PlUq1jt5ltNbPNZrYhXleqtQwzsxfM7N/NbLuZXVGKtZjZ+Pj1SH0dNrO7S7SWb8X/v79tZs/GOVBydQCY2TfjOraZ2d3xuoLUUvSBDhwDrnb3ScBkYE58Jc29wE/cfRzwk3i5FHwT2J62XKp1AMx098lp19OWai3fAX7k7hOASUSvT8nV4u5N8esxGagiukDhh5RYLWY2CvgGUO3ulxJdjHETJVYHgJldCvwZ0R33k4AvmNk4ClWLu5fMFzAQ+CXRnapNwIXx+guBpr4eX4LxV8Yv3tXAy/G6kqsjHutu4Pxu60quFuBc4NfEFwiUci3dxj8LeK0Ua+H3d56fR3Ql3stxPSVVRzzOLxI90DC1/N+BvyxULaVwhp6aptgM7ANecfc3gBEeX+sef/94Hw4xqaVEL2b6Iy9KsQ6I7gReZ2Yb40c6QGnW8gngfeB78VTYE2Y2iNKsJd1NwLPxzyVVi7u3AA8B7xHd0X/I3ddRYnXE3gauMrPhZjaQ6PLu0RSolpIIdHfv9OjPyEqgJv4zpqSY2ReAfe5eKs+ayGaGu08letLm183sqr4eUI76AVOBx919CvAhJfCn/KnENwBeBzzf12PJRTyffD0wFhgJDDKzP+3bUeXG3bcD3wZeAX4EvAV0FKq/kgj0FHc/CDQCc4D/SD3RMf6+r+9GlsgM4Doz2030xMqrzexpSq8OANx9b/x9H9E8bQ2lWUsz0Bz/1QfwAlHAl2ItKdcCv3T3/4iXS62WzwG/dvf33b0dWAV8mtKrAwB3f9Ldp7r7VcABYCcFqqXoA93MLjCzYfHPFUQv9r8TPW7gy3GzLwMv9skAE3L3Je5e6e5jiP4c/hd3/1NKrA4AMxtkZkNSPxPNb75NCdbi7v8X2GNm4+NVnwXeoQRrSTOf30+3QOnV8h5wuZkNNDMjek22U3p1AGBmH4+/XwT8V6LXpiC1FP2domY2EfhHone6zwGec/f/aWbDgeeAi4j+A/iiux/o+UjFw8xqgXvc/QulWIeZfYLorByiKYtn3P3BUqwFwMwmA08A/YF3iR5dcQ6lWctAojcUP+Huh+J1Jfe6xJcn30g0PbEJWAAMpsTqADCznwPDgXbgL9z9J4V6TYo+0EVEJJmin3IREZFkFOgiIoFQoIuIBEKBLiISCAW6iEggknxItMgZF1/W9ZN48Q+BTqJb9AFq3P2jPhlYBvFlqB+5+y/6eChyllOgS1Fy9w+Inq6Jmd0PtLr7Q301HjPr5+493bJdC7QCiQPdzMrcvTMfYxNJ0ZSLlAwzqzKzn8UPBGtIu3W60cweNrNX4+eZTzOzVfGzph+I24yJn3f+j2a2JX7++cAEx/07M/sZ8E0z+xMzeyN+iNePzWyEmY0B7gS+ZdEzyK80s6fMbF7auFvj77Vm9lMzewbYGj90rs7M1sdj+uoZ/YVKcBToUioMeASY5+5VwArgwbTtH8XPylhGdBv114FLgdvj6RuA8cByd58IHAa+ZmblWY47zN0/4+7/APwrcHn8EK964C/dfXfc58MePYv851nqqAH+m7tfAnyF6EmC04BpwJ+Z2dje/2pEIppykVLxMaKAfiV6vAdlRI9WTVkTf98KbEs9mtTM3iV6XOlBYI+7vxa3e5roQxR+lOW4K9N+rgRWxmfw/Ymeo95bb7p7ar9ZwMS0s/mhwLgcjyuiQJeSYURBfUUP24/F34+n/ZxaTv133v05F57guB+m/fwI8L/cfU38Ruj9PezTQfzXb/xwqf49HM+Au9y9oYfjiPSKplykVBwDLjCzKwDMrNzMPtXLY1yU2p/oiYT/SvTJMUmPOxRoiX/+ctr6I8CQtOXdRB8BB9Fzvct7OF4D8OfxtA9m9p/jp1eK5ESBLqXiODAP+LaZvQVsJnpGdm9sB75sZluIPt7s8fjyx6THvR94Pn563v609S8BN6TeFAX+D/AZM3uT6OMSPzzpSJEniB7V+0szexv4LvqrWU6DnrYoZ4X4apSXPfrQYZEg6QxdRCQQOkMXEQmEztBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRALx/wGco9P8s4GD7wAAAABJRU5ErkJggg==\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)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing and plotting uncertainty" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxlUlEQVR4nO3deXgUVb4+8Leqes2edHZ2UCAiiwLGURQFJAwGRL2Ig1fnqoNzFUWZ8T4oOgrK4OAsijguozPj+IPxzsN1g6iA6DCAYkBlCBAWhUBYOltn7/RadX5/NGkT1u7K2p338zwC3amq/h47ydt1zqlTkhBCgIiIKExyVxdARESRiQFCRES6MECIiEgXBggREenCACEiIl0YIEREpEunBMiyZcswYcIEDBkyBAcPHjzrNqqqYvHixZg0aRJuuOEGrF69ujNKIyIinTolQCZOnIhVq1ahV69e59xm7dq1KC0txYYNG/CPf/wDK1aswPHjxzujPCIi0qFTAmTMmDHIyso67zYff/wxZs6cCVmWkZKSgkmTJmHdunWdUR4REenQbcZA7HY7srOzg4+zsrJQVlbWhRUREdH5dJsAISKiyGLo6gKaZWVl4eTJkxgxYgSAM89IQlVT44SmRefyXjZbHByOxq4uo8NEc/uiuW0A2xfJZFlCcnKsrn27TYBMmTIFq1evxuTJk1FbW4uNGzdi1apVYR9H00TUBgiAqG4bEN3ti+a2AWxfT9QpXVhLlizBtddei7KyMtx999248cYbAQBz5szB7t27AQA33XQTevfujcmTJ+O2227D3Llz0adPn84oj4iIdJCibTl3h6Mxaj8ppKXFo7KyoavL6DDR3L5obhvA9kUyWZZgs8Xp27edayEioh6CAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBdDZ71QSUkJHnvsMdTW1iIpKQnLli1D//79W23jcDjw+OOPw263w+fz4corr8STTz4Jg6HTyiQiohB12hnI008/jdmzZ2P9+vWYPXs2nnrqqTO2ee211zBo0CCsXbsWa9euxd69e7Fhw4bOKpGIiMLQKQHicDhQXFyM/Px8AEB+fj6Ki4tRXV3dajtJkuB0OqFpGrxeL3w+HzIyMjqjRCIiClOn9A3Z7XZkZGRAURQAgKIoSE9Ph91uR0pKSnC7Bx54AA899BDGjRsHl8uFO+64A6NHjw7rtWy2uHatvbtJS4vv6hI6VDS3L5rbBrB9PVG3GlxYt24dhgwZgr/97W9wOp2YM2cO1q1bhylTpoR8DIejEZomOrDKrpOWFo/KyoauLqPDRHP7orltANsXyWRZ0v3Bu1O6sLKyslBeXg5VVQEAqqqioqICWVlZrbZbuXIlpk+fDlmWER8fjwkTJqCwsLAzSiQiojB1SoDYbDbk5OSgoKAAAFBQUICcnJxW3VcA0Lt3b2zevBkA4PV6sW3bNlx88cWdUSIREYWp02ZhLVq0CCtXrkReXh5WrlyJxYsXAwDmzJmD3bt3AwAWLlyIb775BtOmTcOMGTPQv39/3HbbbZ1VIhERhUESQkTVgAHHQCJXNLcvmtsGsH2RrNuPgRARUfRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKRLyAHy2Wefwe/3d2QtREQUQUIOkOXLl2PcuHF45plnsGvXro6siYiIIkDIAbJmzRq89dZbMJvNeOihh5CXl4dXXnkFx48fD2n/kpISzJo1C3l5eZg1axaOHDly1u0+/vhjTJs2Dfn5+Zg2bRqqqqpCLZGIiDqRJIQQ4e4khMC2bdvwm9/8Bt999x0uv/xyzJo1C/n5+ZDls2fSXXfdhVtvvRU33XQTPvzwQ7z77rt4++23W22ze/duLFiwAH/729+QlpaGhoYGmEwmmM3mkGtzOBqhaWE3KSKkpcWjsrKhq8voMNHcvmhuG8D2RTJZlmCzxenbN9wdSktL8cc//hGLFi2Cx+PBvHnzMHPmTKxatQrz5s076z4OhwPFxcXIz88HAOTn56O4uBjV1dWttnvrrbdwzz33IC0tDQAQHx8fVngQEVHnMYS64apVq/Dhhx/i6NGj+PGPf4znn38eo0aNCn49Ly8PV1111Vn3tdvtyMjIgKIoAABFUZCeng673Y6UlJTgdocOHULv3r1xxx13oKmpCTfccAPuv/9+SJKks3lERNRRQg6QzZs34+6778bEiRNhMpnO+LrVasWKFSvaVIyqqjhw4AD++te/wuv14mc/+xmys7MxY8aMkI+h91QsUqSlxXd1CR0qmtsXzW0D2L6eKOQAeemllyDLMoxGY/A5n88HIUQwUMaNG3fWfbOyslBeXg5VVaEoClRVRUVFBbKyslptl52djSlTpsBkMsFkMmHixIkoKioKK0A4BhK5orl90dw2gO2LZJ0yBnLPPfdg7969rZ7bu3cv7r333gvua7PZkJOTg4KCAgBAQUEBcnJyWnVfAYGxka1bt0IIAZ/Ph6+++gpDhw4NtUQiIupEIQfIgQMHMHLkyFbPjRgxAvv37w9p/0WLFmHlypXIy8vDypUrsXjxYgDAnDlzsHv3bgDAjTfeCJvNhqlTp2LGjBm46KKL8B//8R+hlkhERJ0o5C6shIQEVFVVBWdIAUBVVRWsVmtI+w8aNAirV68+4/k33ngj+G9ZlvH444/j8ccfD7UsIiLqIiGfgUyePBm//OUvcfDgQbhcLhw4cAALFizAj3/8446sj4iIuqmQA2T+/PkYNGgQZs6cGbxwcMCAAfjFL37RkfUREVE3FfaV6EII1NTUIDk5uVten8FZWJErmtsXzW0D2L5I1pZZWCGPgQBAQ0MDSkpK4HQ6Wz3/ox/9SNeLExFR5Ao5QN577z0888wziImJgcViCT4vSRI+++yzDimOiIi6r5AD5IUXXsDy5csxfvz4jqyHiIgiRMiD6KqqnvNKcyIi6nlCDpA5c+bg1VdfhaZpHVkPERFFiJC7sN566y1UVVXhzTffRFJSUquvbdq0qZ3LIiKi7i7kAPntb3/bkXUQEVGECTlArrjiio6sg4iIIkzIYyBerxcvvPACJk6ciNGjRwMAtm7dipUrV3ZYcURE1H2FHCBLly7FwYMH8bvf/S54BfrFF1+Md955p8OKIyKi7ivkLqyNGzdiw4YNiImJgSwHcicjIwPl5eUdVhwREXVfIZ+BGI1GqKra6rnq6uozZmQREVHPEHKATJkyBQsWLMCxY8cAABUVFXjmmWdw4403dlhxRETUfYW1nHuvXr0wffp01NfXIy8vD+np6Zg7d25H1kdERN1U2Mu5A4GuKy7n3vmieUlpILrbF81tA9i+SNYpy7k3d101a7mke58+fXS9OBERRa6QA+SGG26AJEloecLSfAayb9++9q+MiIi6tZADZP/+/a0eV1ZW4uWXX8aYMWPavSgiIur+Qh5EP11aWhqeeOIJ/OEPf2jPeoiIKELoDhAAOHz4MFwuV3vVQkREESTkLqzZs2e3mnXlcrnw/fffcxovEVEPFXKAzJw5s9Vjq9WKoUOHon///u1dE5EuAkD3m1hOFL1CDpCbb765I+sgajsB+DQNRqVNPbNEFKKQA2T58uUhbffwww/rLoaorXw+lQFC1ElCDpCjR49iw4YNuPTSS9GrVy+cPHkSu3fvxuTJk2E2mzuyRqKQef0aYiUg/PUViChcIQeIEAK///3vkZeXF3xuw4YNWLduHZ577rkOKY4oXH5VwKdqMMg8CyHqaCH/lG3evBmTJk1q9dzEiRPxr3/9q92LItJLEwJen9bVZRD1CCEHSL9+/bBq1apWz/39739H3759270oorZwe1V0w3U+iaJOyF1YS5YswYMPPog333wzeCdCg8GAFStWdGR9RGHz+zX4VQ0Ku7GIOlTIAXLJJZdg/fr12LVrFyoqKpCWloZRo0bBaDR2ZH1EYdOEgNevwWpigBB1JN0/YWPHjoXP50NTU1N71kPULjwedmMRdbSQz0AOHDiA+++/HyaTCeXl5Zg6dSp27NiB999/Hy+++GIHlkgUPi+7sYg6XMg/XYsWLcK8efOwbt06GAyB3Bk7diy++eabkPYvKSnBrFmzkJeXh1mzZuHIkSPn3Pbw4cMYOXIkli1bFmp5RK1oQsDD2VhEHSrkAPn+++9x0003AfjhRlIxMTHweDwh7f/0009j9uzZWL9+PWbPno2nnnrqrNupqoqnn376jCnDROFq8vi6ugSiqBZygPTq1Qt79uxp9VxRUVFI03gdDgeKi4uRn58PAMjPz0dxcTGqq6vP2PZPf/oTrrvuOi7SSG2mqoHBdCLqGCGPgTz88MP4+c9/jttvvx0+nw+vv/46/vd//xfPPvvsBfe12+3IyMiAoigAAEVRkJ6eDrvdjpSUlOB2+/fvx9atW/H222/jlVde0dEc6L45fKRIS4vv6hI6VFvap2oCfllqtYyJxWiALcnSDpW1Hd+7yBbt7dMj5AC5/vrr8cYbb2D16tUYO3YsTpw4gRUrVuDSSy9tl0J8Ph9+9atf4bnnngsGjR4ORyM0LToXQkpLi0dlZUNXl9Fh2to+IYCaOler91+WJPg8ni4fTOd7F9miuX2yLOn+4B1SgKiqiry8PHz88cdYtGhR2C+SlZWF8vJyqKoKRVGgqioqKiqQlZUV3KayshKlpaW47777AAD19fUQQqCxsTGksxyis9GEQJNHRbyVs7GI2ltIAaIoChRFgcfjgclkCvtFbDYbcnJyUFBQgJtuugkFBQXIyclp1X2VnZ2NwsLC4OMVK1agqakJCxYsCPv1iFpyefyItRgg88IQonYV8seyu+66C4888gi2b9+O0tJSHDt2LPhfKBYtWoSVK1ciLy8PK1euxOLFiwEAc+bMwe7du/VVTxQCTRNw+9SuLoMo6khCnP/OCZWVlUhLS8PQoUMDO0gSWu4iSRL27dvXsVWGgWMgkas9xkAqTxsDaWY0SLAlWNtSXpvwvYts0dy+toyBXPAMpPn+H/v378f+/fsxYcKE4L/379/frcKD6Fz8nNJL1O4uGCCnn6Ds2LGjw4oh6ihCAG6vv6vLIIoqFwwQ6bSBxwv0eBF1W26PCjVKuzeJusIFZ2GpqoqvvvoqGBynPwaAH/3oRx1XIVE70YSA2+tHrIW3ICBqDxcMEJvNhoULFwYfJyUltXosSRI+++yzjqmOqJ253AwQovZywQD5/PPPO6MOok7h1wRcXj+sppAXYSCic+BPEUW8okNVWFdYCrdXRazFgDFD0zGkb/I5t3e6fAyQCNX8XlfVuZGaaMGU3L4YMSi1q8vqsbi+A0W0okNVWPXpQdQ6vbBaFNS7fFjzRQkOlNaccx+/KtDk4YysSNPyvY6xGFDr9GLVpwdRdKiqq0vrsRggFNHWFZZCUWSYjQokSYLJoEBRZGzZdfK8+zW5fBDgjKxIcvp7bTYG3ut1haVdXVqPxQChiFZV54bJ0Prb2KjIqGk4/43O/JpAk5tnIZHkbO+1ySCjqs7dRRURA4QiWmqi5YwrzH2qhuR48wX3dbr90HhdU8Q423vt9WtITewe93vpiRggFNGm5PaFqmrw+FQIIeD1q1BVDdeMzL7gvpom4ORZSMQ4/b32+ALv9ZTcC98VlToGp6JQRGuegbOusBQut4oEqxETLut13llYLTW5AzOyDAqXeu/uWr7XnIXVPTBAKOKNGJSKEYNSz7sa77kIATQ0eZEcbwLAEOnumt9r6h7YhUU9nsenwu3jSr1E4WKAEAFodHo5oE4UJgYIEQLTejmgThQeBgjRKU1uH/wqz0KIQsUAITqleUCdiELDACFqweNT4eKdC4lCwgAhOk1jk48D6kQhYIAQnUbVBBpdPki8LITovBggRGfh8vjh8fLaEKLzYYBQVHDUubHx62NodPna5XhCAPVNHnZlEZ0HlzKhqLDmixJsKbIjOd6Mn04ZivRka5uP6VcD14bEW3kPdaKz4RkIRYWrh2fBZAzcB+S1D/fg++N17XJcXhtCdG4MEIoKg/skYcHsyxEfY4Tbq+KtT/ahsLi8zcfltSFE58YAoajRPzMBc28ZjixbDDQBfLi1BGu/PAI1jNV5z8bjU9HEa0OIzsAAoaiSFGfGfdOHIadf4H4g2/aU4W+f7IfL07YAaGzyQdM4K4uoJQYIRR2zUcEdNwzGtafuSvj9iTq88v4elFc36T6mpgnUO33gaAjRDxggFJVkWcKU3L647fqLYFAkOOrdePXDPdhTUq37mG6fygsMiVpggFBUG3VxKu6bPgyJsSZ4fRr+/ulBrN9eGtZdC1tqcvvg8qrtXCVRZGKAUNTrnRaHubcMx4CsBADAv/59En/9ZJ+uiw6FAOobvfD6OR5CxAChHiHOasQ9N+Zg3PAsAMChE/V4+b3dOFrWEPaxNCFQ2+iB26eyO4t6tE4LkJKSEsyaNQt5eXmYNWsWjhw5csY2f/zjH3HjjTdi+vTpuOWWW7Bly5bOKo96AEWWMPVH/fCTSRfDZJRR7/TijbXF2LzrZNhLlmiaQF2jBw0cE6EerNMC5Omnn8bs2bOxfv16zJ49G0899dQZ24wYMQL/93//hzVr1mDp0qWYP38+3G53Z5VIPcTwgTbMvXk4MpKt0ITAusJS/L91B8Lu0hICcDb5OLBOPVanBIjD4UBxcTHy8/MBAPn5+SguLkZ1desZMddccw2s1sAaRkOGDIEQArW1tZ1RIvUwaUlW3H/zpRgzJA0AcOBYLV5+twiHToa3BIpA4BqRQPhwki/1LJ0SIHa7HRkZGVAUBQCgKArS09Nht9vPuc8HH3yAvn37IjMzszNKpB7IZFBwy/hBuG3CRYEurSYf/lKwD+u3l8Kvhj5I3hwi9U0+CIYI9SDdcjXe7du3Y/ny5fjLX/4S9r42W1wHVNR9pKXFd3UJHaot7VM1Ab8sIdwV2CdcEYtLL07Dn9fsxVF7Pf7175MosTfgnunDkGmLDetYmgwkxJoRYzlzBV++d5Et2tunhyREx9/wwOFwIC8vD4WFhVAUBaqqIjc3Fxs2bEBKSkqrbXfu3IlHHnkEr7zyCoYNG6bjtRp1z/Hv7tLS4lFZGf6soUjR1vYJAVTWuXS//6qmYePXx7H53ychABiUwMWIVw7LhBzGIIckAVazAfExRkgI7Mf3LrJFc/tkWdL9wbtTurBsNhtycnJQUFAAACgoKEBOTs4Z4VFUVIT58+fjpZde0hUeRG2hyDLyruiLn027BElxJvhVgYIvj+IvH+1DTUPokzmEAJrcfjjq3fD6VHBshKJVp5yBAMChQ4fw2GOPob6+HgkJCVi2bBkGDhyIOXPmYN68eRg+fDhuvfVWnDhxAhkZGcH9nn/+eQwZMiTk1+EZSOTq6jOQltxePz768ii+OVgJADAZZUzJ7YsrcjLCPhsxGxT07ZWE+npXm+vqrvi9GbnacgbSaQHSWRggkas7BUizjV8fw7/+fTK4JHxGshVXXZqJXd9XoabBg+R4M64ZmY0hfZPPexxbSizcLi9iLYawAggAig5VYV1hKarq3EhNtGBKbl+MGJSqu03tac3Ww9iw4zjcPhUWo4LJY3tj+riBXV1Wu4vmn71u34VFFIkOlNZg53eVSIwzwWoOzCAsr3Hh/S0lKKtxwWxSUO/yYc0XJThQWnPeYwkATpcPjjo3msJYWr7oUBVWfXoQtU4vYiwG1Dq9WPXpQRQdqmpL09rFmq2HsebLI/D4VBjkwH1T1nx5BGu2Hu7q0qiTMECIzmHLrpNQFBkWkwHJ8RakJJiDX2ty+1FV5wYEoCgytuw6GdIxVU2g3ulFVZ0Lbp96wdGRdYWlUBQZZqMCSZJgNipQFBnrCkvb0LL2sWHHcUiQoMgSJEkO/A0JG3Yc7+rSqJN0y2m8RN1BTYMHFvMPPyIWkwESPMFf+n5VoKrODatZgc8X3gq9flWgtsEDgyIhxmJEjPnsP4pVdW7EWFp/zWSQA+HVxdxePxS5dXecLAWep56BZyBE55Acb4bvtAsKDYoEoyIhNdECoxL48XF5VDQ0+fBVcVnY4y9+NXBGEpixpZ2xJEpqouWMlX+9fg2piZbwG9TOLCYDTm+uJgLPU8/AACE6h2tGZkNVNXj9KoQQ8PpVmIwKTCYDIAG2RDNirQZICPziXLP1CP74/m4cPlkf9mv5/BpqGt2obfS2unXulNy+UFUNHl+gBo9PhapqmJLbtx1bqs/ksb0hIKBqAkJogb8hMHls764ujToJPyoQncOQvsmYjsBYSPOMqxuv7Ae0eC4jyYoxV6bj0Il6fHOwEnZHE94sKMaw/imYktsXtjDOFIQAXB4/PD4VMRYDzEYFIwalQpKAT77qfrOwmmdb9YRZWHR2nMYbQaJ5KiHQPafxhuN4RSMKth1BaXkjgMDy8bmXZOD6y3uhT3YSqqudYR1PkgBFkRBjNsJiUsKe/tuZ+L0ZuTiNl6gb6J0eh59PH4ZZEy5CUpwJqibw5Z4y/O6df+PjL0rgCXOgXQjA7w+MkVTXueF0+6BqZ46TEHUVdmERtSNJkjDyolRc0j8F2/aWYdPOE3B7VazZchiffX0M143KxhU5GTAawvvs5tcEGpp8cLr8MBllWC0GGBW5W5+VUPRjgBB1AKNBxrUjszFmSDo27zqBbXvL4XT58NG2o9hSZMd1o7IxZmg6DEp4QaIJAbdXhdurQpYlmA0KLBYFJoMcXLiRqLNwDCSCRHM/LBD5YyDnIxkVvP/5d/jmQGVwWZSEWBOuHZmFMUPTYTIobTq+LEuwmBRYjAqMBqXTu7n4vRm52jIGwjMQok6QHG/BjGsGYvyobPxz50l8e6AS9U4vCr48in/uPImrL81E7iUZsJ7jgsIL0TSBJrcfTW4/DLIEs9kAs0GG0cgzE+o4DBCiTpQcb8Et1w7E9Zdl41//PolvDlTC6fJhw47Aoo1jh6bjquGZSIozX/hg5+DXBPwuH5xAsJvLZJJhNMgwyJw3Q+2HAULUBZrPSCZc3htf7LajcF85PD4VW3fb8eUeOy4daMPVwzPRJ71td8HTNAGX1w+X99S0YFmCyajAbFJgkCUYFDnsOzgSNWOAEHWhhFgTfnxlP1x3WS8UFpdj254yNLh8KDrkQNEhB3qnxeJHl2Zi+EBb2APupxMisHSKXw10dcmSBFkBzEYDTEYZBlk6tbYVu7woNBxEjyDRPJAHRPcgekpKbEgXEvpVDUWHHPhitx12R1Pw+RiLAWOGpGFsTgZsCR2zDlZzoJgMCowGGYosw6AEQuVCvyX4vRm5OIhOFCUMiozLB6fhsotTcaSsAdv2lqG4pBpNbj8277Jj8y47BvVKwJgh6bikf0rY15OcjyYEND/g9wdW05UQuK5FlgGjosBolGFQZBgNEmTpwqFC0Y8BQtQNSZKEAVkJGJCVgHqnFzv2V+Dr/RWoc3px6EQ9Dp2oh8WkYORFqbh8cCp6p8VBaue5uwKAEAKaCvjVH8ZRJEmCUZFhMikwngoU6pkYIETdXEKsCRNH98b1l/XCd8drsWN/BfYfrYXbq6KwuByFxeVITbRg5EWpGHmRDamJ1g6rRYhAqHg0FR6fGjhLkSVIBgMa3b5TM70kyLIMWQLPUqIcA4QoQsiyhCF9kzGkbzIamrzY9b0D3x6sRFl1E6rq3Pjsm+P47Jvj6JUai+GDbBg+MAXJ8R173xABQGgCXr+GxiYfgMBZiixJUJTA1GFZCix3DwHIcuAOjgZZgqKwKyzSMUCIIlB8jAnjRmRh3Igs2B1O/Pu7Kuw65EC904sTVU6cqHJiXWEpeqXGYtiAFFzSPwXpyR13ZtKSEIAqBFRNhfccC0g2d4XJMmCQAwP2khx4rjmA5FMbKnLzbXN5RtPdMECIIlyWLRZZtljk5fbF0bIGFB1yYG9JNRpdvmCYbNhxDKmJFuT0S8bQfsnomxF/xu1oO1NzV5imAX6oAM4RNPghVIyKDMUQuPd689mLLAXOzCQwYLoCA4QoSsgtBt6nXdUfR8oasKfEgX1HalDn9KKqzo0tRXZsKbLDalZwce8kDO6ThIt6JyIhxtTV5Z9V80A+BKBqKtDijKZ5zoAkSafGYgJnM4ZTIdMcLIocCJnmlYsZMu2HAUIUhWRZwsDsBAzMDoTJiSon9h2pwf7SGtgdTXB51ODFigCQmRKDi3onYlB2AvpnJcBsbNvijp2hOQiCl7I1n814W4TMqT+kU2crBlmGLLcOEknGqbMaOXhWA/zQzSZJCF5bdPpEt54eRgwQih4SYDLK8Ps1aFrgugYK/BLsnRaH3mlxuGFsH9Q2enCgtBbfHa/F9yfq4PVpKKtuQll1E7YW2SFLEnqnx2JgdiIGZMWjb3o8zKbuHyhnI079IYSABsCvnv+mXlLwj+bHEiABfklGba3rVBgFnpdlCYmxph59TxYGCEUNCUBSbGARQlUTUFUBv6bB79fg1wSE0HCqNyT4i0KWArOEZEmCJAc+fTZ/6pQQ2FYL9tcLCC3wbxUCQhPB4wkhIubTaFKcGbmXZCD3kgz4VQ3HKhrx3fE6HDpRh+OVjdCEQGl5I0rLG7FpJyBLQFZqLPplxKNfZjz6ZsQjMbZ7dnm1lQj+0fxYnAogBJfhb/5KTw6OZgwQikrNM3dMkAHzmV0PLYX7i7/lsZr39asavGogoDRNQPVr8GsauuGqKq0YFDk4boKxfeD2+lFib8Dhk3UosTfA7nBCE8CJSidOVDrx5Z4yAIFrU/qkxaFPehx6pcfCGqt/9WCKXAwQ6hHa8+zgbMcyKHKrxQ6bQ0Y79ZE2PsEMr9sLSZKgqhq8PvWHcBGnPum2c516WEwG5PRLRk6/ZACA2+vH0bKGwH/lDThe4YRP1VDv9GKvsxp7j1Sf2nMfUhMtyE6NRbYtFlmpMciyxSLOauy6xlCHY4AQdYDgAO2pP61mI+IsZ/4y1YJdXwKaAHyqBq9HhdevdYsxHIvJELx4EQBUTUOZowmlFY04XtGIYxWNcNS5IQBU1blRVecODswDQLzViExbDDJSYpCRbEVGcgzSkq0RMUhPF8YAIepCcvNgCyQoCFzrEGs2wK8KuL1+uNx+qN1ofEWRZfRKi0OvtDhgWOA5t9ePBo+G/SVVsFc14USVE1V1LggBNLh8aDheh++O17U6TlKcCWlJVqQlWZGaZEFqohWpiRYk9PBB6UjDACHqZoQIjOHEWoyIsRjhVzU43T54vGq3CZKWLCYDsjNjkRb/w8C6z6+hvKYJZY7A7K7ymiaUVbvgdAWWO6lt9KK20XtGsBgVGSkJZqQkWGBLsCA5wYyUeDOSEyxIjjO36+rD1HYMEKJuTELgl2pynBk+vwZVCHh9GjxeP1St+5yZnM5okINTh1tyun2oqHGhosaFqloXKutcqKx1o7bBA4FAF155jQvlNa6zHjfOakRSnAlJ8WYkxZqRGGdCYpwZibEmJMaaEGc1Bq/zoI7HACGKAEKcGqgHYDYoiLcaoWoaVFXAp2rw+DT4Va1b3kyrpViLEQOyjIFZXy34/BqqG9xw1LnhqA/8XdPgQXW9B7WNnuAU2kaXD40uH45Xnv3mXJIUGHeJjzUh3mpCfIzx1H+BcImPMSLOakSs1QiTQW73JfB7GgYIUYRSZBmKDJiMCuKsgesU/KqAJsSp62ACAeNXtRaD9d2T0SAjIzkGGckxZ3xN0wTqm7yoafCgpiEQKLWNXtQ1elDn9KKu0QvPqSVOhADqm3yob/IBOP8dII2KjBiLAbFWI2ItBsRajLBaDIgxGxDT4m+r2QBVluH1+GE2KRyjaYEBQhQFhAgMyJtOu7lT8wKDgbMVwC8CF1b6ToVLJFwAKcsSkuLMSIozY0DW2bdxe/2od/pQ7/SivsmLhiYv6p0+NLi8aGjyobEp8G+vTwvu41O1QAA5vSHXIgEwmxRYTApiLUbcPH4gLr84rY0tjFydFiAlJSV47LHHUFtbi6SkJCxbtgz9+/dvtY2qqliyZAm2bNkCSZJw3333YebMmZ1VIlHUaQ6HvSXVWFdYCke9GxnJMZiS2xeX9E/B3z89gOIjtYHl11UNQ/smYmi/FGwtssNR50ZCrBHjRmRjSN9kHCitwZZdJ1HT4EFyvBnXjMzGicpGbC0qg8evwmxQMG5EJiaM7nPWWs62/7mO2zxtONRjfLnbHqjDp8JsDNQx7er+Z+zv9avY/b0D2/aWoc7pRYzZgL4ZcXD7VBwta4DHq0KWJZiNCvyagNvjb3UxqADg9qpwe1XUNnrx14/2YWP6MUzJ7YsRg1Lb+G5FHkmIzvn8cdddd+HWW2/FTTfdhA8//BDvvvsu3n777VbbfPDBB1i7di3eeOMN1NbWYsaMGfj73/+O3r17h/w6Dkdjt+8H1istLR6VlQ1dXUaHieb2dWXbig5VYdWnB6EoMkwGGV6/BlXVYIs3Y/+xwCwoWQ7cptagSLAYZaSmxMBsUALdYarApQOS8e13VaeWf5Hg9Wmob3TD5dUAISBJAj6/gCqAiZf3OiNEDpTWYM0XJYGbTCnyqTMgDaMHp+Gbg5VnPD/96gFnhMi5jtEvPQ67DldDAoI3rxIIvQ6X2w9IEqxm5YwaBvdJgsenwmQxoayiAS6PH4dO1GH7/nKYDArSk62ob/JBVTXcccPgiAwRWZZgs8VdeMOz6JQzEIfDgeLiYvz1r38FAOTn5+PZZ59FdXU1UlJSgtt9/PHHmDlzJmRZRkpKCiZNmoR169bhZz/7WcivFe0zMNi+yNVVbSssLkd6SgxMhh8u3vP6VVTWuJCebG25dmBwGagYc+Cix8AeKr4+WIXEWBMsJuWH+3MYFcgILGDZsm01jT4kJ1gC3WNa4Cr70opGXNQnCYokQ5y67t7rU1Fa4USfjHgYFSXwvAB8qor9pbUYPsh26ir9gOIj1ciwxcCkKMHnvH4V9moXMpKtrWrQNIEDx+qQl9uv1f+LvSXVZ/y/cNQFZnzZWtwK2OtXsbekGpcOtMFkVJCUFAPzqRnEu76vQp+MeMRZTbCaFcRa/fD6A7cXHhWB3Vlt+b7slACx2+3IyMiAogTeNEVRkJ6eDrvd3ipA7HY7srOzg4+zsrJQVlYW1mslJ8e2T9HdlN5PCpEimtvXVW1beM+VXfK6Lf3ijjFtPsZTg9r+y/npNhwj0xbb5mNEG16VQ0REunRKgGRlZaG8vBzqqbX4VVVFRUUFsrKyztju5MmTwcd2ux2ZmZmdUSIREYWpUwLEZrMhJycHBQUFAICCggLk5OS06r4CgClTpmD16tXQNA3V1dXYuHEj8vLyOqNEIiIKU6fNwjp06BAee+wx1NfXIyEhAcuWLcPAgQMxZ84czJs3D8OHD4eqqnjmmWfwxRdfAADmzJmDWbNmdUZ5REQUpk4LECIiii4cRCciIl0YIEREpAsDhIiIdGGAEBGRLhG7Gu8DDzyA48ePQ5ZlxMTE4Fe/+hVycnJCWrQxUrz88stYsWIF1q5di8GDB0dN2yZMmACTyQSz2QwAePTRR3HNNddETfs8Hg+WLl2Kbdu2wWw2Y9SoUXj22Wejon3Hjx/H3Llzg48bGhrQ2NiI7du3R0X7AOCf//wnli9fDiEENE3DQw89hMmTJ0dN+zZt2oTly5fD7/cjMTERzz33HPr06aOvfSJC1dfXB//96aefihkzZgghhLjzzjvFBx98IIQQ4oMPPhB33nlnl9TXVnv27BH33nuvuO6668SBAweEENHTtuuvvz7YppaipX3PPvus+PWvfy00TRNCCFFZWSmEiJ72tbRkyRKxePFiIUR0tE/TNDFmzJjg9+e+ffvEqFGjhKqqUdG+2tpaccUVV4jDhw8LIQLtuOeee4QQ+t6/iA2Qlt5//31x8803i6qqKjF69Gjh9/uFEEL4/X4xevRo4XA4urjC8Hg8HnHbbbeJ0tLS4C/baGmbEGcPkGhpX2Njoxg9erRobGxs9Xy0tK8lj8cjcnNzxZ49e6KmfZqmiSuuuEJ8/fXXQgghtm/fLiZPnhw17du1a5eYOnVq8HFNTY0YPHiw7vZFbBcWADzxxBP44osvIITAm2++GfKijd3d8uXLMX36dPTp88NS1NHStmaPPvoohBAYPXo0fvGLX0RN+44dO4akpCS8/PLLKCwsRGxsLB5++GFYLJaoaF9Ln3/+OTIyMjBs2DDs2bMnKtonSRJefPFFPPDAA4iJiYHT6cTrr78eNd+fAwYMQFVVFYqKijBixAisXbsWgP7fLxE9iP7rX/8amzZtwvz58/H88893dTntYufOndi9ezdmz57d1aV0mFWrVmHNmjV49913IYTAM88809UltRu/349jx47hkksuwXvvvYdHH30UDz30EJqamrq6tHb37rvv4tZbb+3qMtqV3+/H66+/jldeeQX//Oc/8eqrr2L+/PlR8/7Fx8fjhRdewHPPPYdbbrkFDocDCQkJutsX0QHSbMaMGSgsLERmZmZIizZ2Zzt27MDhw4cxceJETJgwAWVlZbj33ntRWloa8W1r1lyzyWTC7Nmz8e2334a84GZ3l52dDYPBgPz8fADAyJEjkZycDIvFEhXta1ZeXo4dO3Zg2rRpAEJfMLW727dvHyoqKjB69GgAwOjRo2G1WmE2m6OifQBw1VVX4Z133sF7772H//zP/4Tb7UavXr10tS8iA8TpdMJutwcff/7550hMTAx50cbu7L777sPWrVvx+eef4/PPP0dmZib+/Oc/Y+rUqRHfNgBoampCQ0PgznxCCHz88cfIycmJivcOAFJSUpCbmxtcz62kpAQOhwP9+/ePivY1e//99zF+/HgkJwfuGhgt719mZibKyspw+PBhAIE1/KqqqtCvX7+oaB8AVFZWAgA0TcMf/vAH3H777ejVq5eu9kXkWlhVVVV44IEH4HK5IMsyEhMTsWDBAgwbNuycizZGqgkTJuC1117D4MGDo6Jtx44dw0MPPQRVVaFpGgYNGoQnn3wS6enpUdE+INDGhQsXora2FgaDAY888gjGjx8fNe0DgLy8PDzxxBO49tprg89FS/vWrFmDN954A5IUuFPfvHnzMGnSpKhp3xNPPIFvv/0WPp8PV199NRYuXAiz2ayrfREZIERE1PUisguLiIi6HgOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISJeIXguL6HSXXXZZ8N8ulwsmkym4vs/ixYsxffr0ripNtwkTJmDJkiW46qqruroUolYYIBRVdu7cGfx3JPzi9fv9MBg69sewM16DeiZ2YVGPoGka/vSnP2HSpEnIzc3Fww8/jNraWgCBmyQNGTIE7777LsaPH4+xY8finXfeQVFREaZNm4YxY8a0WvDxvffew+23345nn30Wo0ePxpQpU7Bt27bg1xsaGrBw4UKMGzcO11xzDV544YXgGkPN+y5duhRXXHEFVqxYgdLSUtx1113Izc1Fbm4ufvnLX6K+vh4A8D//8z84efIk/vu//xuXXXYZ3njjDRQWFra6AhwIhOWXX34JAFixYgXmzZuHRx99FJdffjnef//989ZEpBcDhHqEt99+Gxs3bsTKlSuxZcsWJCYmnrEK8K5du7Bhwwa88MILWLp0KV577TW89dZb+Oijj/DJJ59g+/btwW2LiorQp08ffPXVV5g3bx4efPDBYCAtWLAABoMBGzZswAcffIAvvvgCq1evPmPfL7/8Evfffz+EEPj5z3+OLVu24JNPPkFZWRlWrFgBAPjtb3+L7OxsvPbaa9i5cyfmzJkTUns/++wzTJkyBV9//TWmTZt2wZqI9GCAUI/wj3/8A/Pnz0dmZiZMJhMefPBBrF+/Hn6/P7jN3LlzYTabMW7cOMTExCA/Px82mw0ZGRkYM2YMiouLg9umpKTgpz/9KYxGI6ZOnYoBAwZg06ZNqKqqwubNm7Fw4ULExMTAZrPhv/7rv/DRRx8F901PT8edd94Jg8EAi8WCfv364eqrr4bJZEJKSgruvvtu7Nixo03tHTVqFCZNmgRZltHY2HjBmoj0YMco9QgnT57E3LlzIcs/fGaSZRkOhyP42GazBf9tNpvPeNzyngkZGRnBxfaAwDLuFRUVOHnyJPx+P8aNGxf8mqZprZbFzszMbFWbw+HAkiVL8PXXX8PpdEIIgYSEhDa1t+VrhFITkR4MEOoRMjMzsXTp0uB9Hlo6fvx42McrLy+HECIYIna7HRMmTAie4Xz11VfnHLhuGTwA8Pvf/x6SJGHNmjVITk7Gxo0bz3uTLavVCrfbHXysqiqqq6vP+Rqh1ESkB7uwqEf4yU9+ghdffBEnTpwAAFRXV2Pjxo26j1ddXY23334bPp8Pn3zyCQ4dOoTx48cjPT0dV199NX7zm9+gsbERmqahtLS01fjJ6ZxOJ2JiYpCQkIDy8nK8+eabrb6empqKY8eOBR8PGDAAHo8HmzZtgs/nw6uvvgqv13vO4+upiSgUDBDqEe666y5MmDAB99xzDy677DLcdtttKCoq0n28ESNG4OjRo7jyyivx4osv4qWXXgreXOn555+Hz+fD1KlTMXbsWMybNy94E5+zefDBB1FcXIwxY8bgvvvuw+TJk1t9/b777sOrr76KMWPG4M9//jPi4+Px9NNP48knn8S1114Lq9V6RrfY6cKtiSgUvB8IUZjee+89rF69Gu+8805Xl0LUpXgGQkREujBAiIhIF3ZhERGRLjwDISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLr8fyTBBLPSycNMAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.set(color_codes=True)\n", "plt.xlim(30,90)\n", "plt.ylim(0,1)\n", "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] } ], "metadata": { "celltoolbar": "Hide code", "kernelspec": { "display_name": "Python 3 (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.9.12" } }, "nbformat": 4, "nbformat_minor": 2 }