diff --git a/module4/data_shuttle.csv b/module4/data_shuttle.csv
new file mode 100644
index 0000000000000000000000000000000000000000..9412ab30682565afb0ca3b4bcc6b0f99f632bbf2
--- /dev/null
+++ b/module4/data_shuttle.csv
@@ -0,0 +1,24 @@
+Date,Count,Temperature,Pressure,Malfunction
+4/12/81,6,66,50,0
+11/12/81,6,70,50,1
+3/22/82,6,69,50,0
+11/11/82,6,68,50,0
+4/04/83,6,67,50,0
+6/18/82,6,72,50,0
+8/30/83,6,73,100,0
+11/28/83,6,70,100,0
+2/03/84,6,57,200,1
+4/06/84,6,63,200,1
+8/30/84,6,70,200,1
+10/05/84,6,78,200,0
+11/08/84,6,67,200,0
+1/24/85,6,53,200,2
+4/12/85,6,67,200,0
+4/29/85,6,75,200,0
+6/17/85,6,70,200,0
+7/2903/85,6,81,200,0
+8/27/85,6,76,200,0
+10/03/85,6,79,200,0
+10/30/85,6,75,200,2
+11/26/85,6,76,200,0
+1/12/86,6,58,200,1
diff --git a/module4/src_Python3_challenger.ipynb b/module4/src_Python3_challenger.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..4d9815312ca0467de0475fe78b8e0e495c58e950
--- /dev/null
+++ b/module4/src_Python3_challenger.ipynb
@@ -0,0 +1,783 @@
+{
+ "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": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "3.6.4 |Anaconda, Inc.| (default, Mar 13 2018, 01:15:57) \n",
+      "[GCC 7.2.0]\n",
+      "uname_result(system='Linux', node='e41afb9ecaf4', release='4.4.0-164-generic', version='#192-Ubuntu SMP Fri Sep 13 12:02:50 UTC 2019', machine='x86_64', processor='x86_64')\n",
+      "IPython 7.12.0\n",
+      "IPython.core.release 7.12.0\n",
+      "PIL 7.0.0\n",
+      "PIL.Image 7.0.0\n",
+      "PIL._version 7.0.0\n",
+      "_csv 1.0\n",
+      "_ctypes 1.1.0\n",
+      "_curses b'2.2'\n",
+      "decimal 1.70\n",
+      "argparse 1.1\n",
+      "backcall 0.1.0\n",
+      "cffi 1.13.2\n",
+      "csv 1.0\n",
+      "ctypes 1.1.0\n",
+      "cycler 0.10.0\n",
+      "dateutil 2.8.1\n",
+      "decimal 1.70\n",
+      "decorator 4.4.1\n",
+      "distutils 3.6.4\n",
+      "ipaddress 1.0\n",
+      "ipykernel 5.1.4\n",
+      "ipykernel._version 5.1.4\n",
+      "ipython_genutils 0.2.0\n",
+      "ipython_genutils._version 0.2.0\n",
+      "ipywidgets 7.2.1\n",
+      "ipywidgets._version 7.2.1\n",
+      "jedi 0.16.0\n",
+      "json 2.0.9\n",
+      "jupyter_client 6.0.0\n",
+      "jupyter_client._version 6.0.0\n",
+      "jupyter_core 4.6.3\n",
+      "jupyter_core.version 4.6.3\n",
+      "kiwisolver 1.1.0\n",
+      "logging 0.5.1.2\n",
+      "matplotlib 2.2.3\n",
+      "matplotlib.backends.backend_agg 2.2.3\n",
+      "numpy 1.15.2\n",
+      "numpy.core 1.15.2\n",
+      "numpy.core.multiarray 3.1\n",
+      "numpy.lib 1.15.2\n",
+      "numpy.linalg._umath_linalg b'0.1.5'\n",
+      "numpy.matlib 1.15.2\n",
+      "optparse 1.5.3\n",
+      "pandas 0.22.0\n",
+      "_libjson 1.33\n",
+      "parso 0.6.0\n",
+      "patsy 0.5.1\n",
+      "patsy.version 0.5.1\n",
+      "pexpect 4.8.0\n",
+      "pickleshare 0.7.5\n",
+      "platform 1.0.8\n",
+      "prompt_toolkit 3.0.3\n",
+      "ptyprocess 0.6.0\n",
+      "pygments 2.5.2\n",
+      "pyparsing 2.4.6\n",
+      "pytz 2019.3\n",
+      "re 2.2.1\n",
+      "scipy 1.1.0\n",
+      "scipy._lib.decorator 4.0.5\n",
+      "scipy._lib.six 1.2.0\n",
+      "scipy.fftpack._fftpack b'$Revision: $'\n",
+      "scipy.fftpack.convolve b'$Revision: $'\n",
+      "scipy.integrate._dop b'$Revision: $'\n",
+      "scipy.integrate._ode $Id$\n",
+      "scipy.integrate._odepack  1.9 \n",
+      "scipy.integrate._quadpack  1.13 \n",
+      "scipy.integrate.lsoda b'$Revision: $'\n",
+      "scipy.integrate.vode b'$Revision: $'\n",
+      "scipy.interpolate._fitpack  1.7 \n",
+      "scipy.interpolate.dfitpack b'$Revision: $'\n",
+      "scipy.linalg 0.4.9\n",
+      "scipy.linalg._fblas b'$Revision: $'\n",
+      "scipy.linalg._flapack b'$Revision: $'\n",
+      "scipy.linalg._flinalg b'$Revision: $'\n",
+      "scipy.ndimage 2.0\n",
+      "scipy.optimize._cobyla b'$Revision: $'\n",
+      "scipy.optimize._lbfgsb b'$Revision: $'\n",
+      "scipy.optimize._minpack  1.10 \n",
+      "scipy.optimize._nnls b'$Revision: $'\n",
+      "scipy.optimize._slsqp b'$Revision: $'\n",
+      "scipy.optimize.minpack2 b'$Revision: $'\n",
+      "scipy.signal.spline 0.2\n",
+      "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n",
+      "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n",
+      "scipy.special.specfun b'$Revision: $'\n",
+      "scipy.stats.mvn b'$Revision: $'\n",
+      "scipy.stats.statlib b'$Revision: $'\n",
+      "seaborn 0.8.1\n",
+      "seaborn.external.husl 2.1.0\n",
+      "seaborn.external.six 1.10.0\n",
+      "six 1.14.0\n",
+      "statsmodels 0.9.0\n",
+      "statsmodels.__init__ 0.9.0\n",
+      "traitlets 4.3.3\n",
+      "traitlets._version 4.3.3\n",
+      "urllib.request 3.6\n",
+      "zlib 1.0\n",
+      "zmq 17.1.2\n",
+      "zmq.sugar 17.1.2\n",
+      "zmq.sugar.version 17.1.2\n"
+     ]
+    }
+   ],
+   "source": [
+    "def print_imported_modules():\n",
+    "    import sys\n",
+    "    for name, val in sorted(sys.modules.items()):\n",
+    "        if(hasattr(val, '__version__')): \n",
+    "            print(val.__name__, val.__version__)\n",
+    "#        else:\n",
+    "#            print(val.__name__, \"(unknown version)\")\n",
+    "def print_sys_info():\n",
+    "    import sys\n",
+    "    import platform\n",
+    "    print(sys.version)\n",
+    "    print(platform.uname())\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import statsmodels.api as sm\n",
+    "import seaborn as sns\n",
+    "\n",
+    "print_sys_info()\n",
+    "print_imported_modules()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Loading and inspecting data\n",
+    "Let's start by reading data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>Date</th>\n",
+       "      <th>Count</th>\n",
+       "      <th>Temperature</th>\n",
+       "      <th>Pressure</th>\n",
+       "      <th>Malfunction</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>4/12/81</td>\n",
+       "      <td>6</td>\n",
+       "      <td>66</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>11/12/81</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>50</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>3/22/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>69</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>11/11/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>68</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>4/04/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>6/18/82</td>\n",
+       "      <td>6</td>\n",
+       "      <td>72</td>\n",
+       "      <td>50</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>6</th>\n",
+       "      <td>8/30/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>73</td>\n",
+       "      <td>100</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>7</th>\n",
+       "      <td>11/28/83</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>100</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>8</th>\n",
+       "      <td>2/03/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>57</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9</th>\n",
+       "      <td>4/06/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>63</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>10</th>\n",
+       "      <td>8/30/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>11</th>\n",
+       "      <td>10/05/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>78</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>12</th>\n",
+       "      <td>11/08/84</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>13</th>\n",
+       "      <td>1/24/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>53</td>\n",
+       "      <td>200</td>\n",
+       "      <td>2</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>14</th>\n",
+       "      <td>4/12/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>67</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>15</th>\n",
+       "      <td>4/29/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>75</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>16</th>\n",
+       "      <td>6/17/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>70</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>17</th>\n",
+       "      <td>7/2903/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>81</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>18</th>\n",
+       "      <td>8/27/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>76</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>19</th>\n",
+       "      <td>10/03/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>79</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>20</th>\n",
+       "      <td>10/30/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>75</td>\n",
+       "      <td>200</td>\n",
+       "      <td>2</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>21</th>\n",
+       "      <td>11/26/85</td>\n",
+       "      <td>6</td>\n",
+       "      <td>76</td>\n",
+       "      <td>200</td>\n",
+       "      <td>0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>22</th>\n",
+       "      <td>1/12/86</td>\n",
+       "      <td>6</td>\n",
+       "      <td>58</td>\n",
+       "      <td>200</td>\n",
+       "      <td>1</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "         Date  Count  Temperature  Pressure  Malfunction\n",
+       "0     4/12/81      6           66        50            0\n",
+       "1    11/12/81      6           70        50            1\n",
+       "2     3/22/82      6           69        50            0\n",
+       "3    11/11/82      6           68        50            0\n",
+       "4     4/04/83      6           67        50            0\n",
+       "5     6/18/82      6           72        50            0\n",
+       "6     8/30/83      6           73       100            0\n",
+       "7    11/28/83      6           70       100            0\n",
+       "8     2/03/84      6           57       200            1\n",
+       "9     4/06/84      6           63       200            1\n",
+       "10    8/30/84      6           70       200            1\n",
+       "11   10/05/84      6           78       200            0\n",
+       "12   11/08/84      6           67       200            0\n",
+       "13    1/24/85      6           53       200            2\n",
+       "14    4/12/85      6           67       200            0\n",
+       "15    4/29/85      6           75       200            0\n",
+       "16    6/17/85      6           70       200            0\n",
+       "17  7/2903/85      6           81       200            0\n",
+       "18    8/27/85      6           76       200            0\n",
+       "19   10/03/85      6           79       200            0\n",
+       "20   10/30/85      6           75       200            2\n",
+       "21   11/26/85      6           76       200            0\n",
+       "22    1/12/86      6           58       200            1"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data = pd.read_csv(\"data_shuttle.csv\")\n",
+    "data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAF9JJREFUeJzt3X2UXXV97/H3d5IACYlAg02VQAFJsVyBCOFJtDfx6Qa7JPUCBbyCl940ZUlul9y2htvVa6m1a1V8qHpFY+SiQldNVRBom14e1Ii0IASM4UHBuYBhEhogBshASGYy3/vH2bN7Mkxmzhlmz5lzeL/WmpWz9/mdne939pz5zN5nn9+JzESSJICuVhcgSZo8DAVJUslQkCSVDAVJUslQkCSVDAVJUqmyUIiIqyPiqYh4YC/3R0R8PiK6I2JDRJxQVS2SpMZUeaTwNWDxCPefAcwrvpYBX6qwFklSAyoLhcy8HfjlCEOWANdkzV3AgRHxuqrqkSSNbmoL/+9DgCfqlnuKdU8OHRgRy6gdTTB9+vQTDz300AkpsFEDAwN0dXXmyzOd2pt9tZ9O7W2i+nrkkUeeyczXjjaulaEQw6wbds6NzFwFrAJYsGBBrlu3rsq6mrZ27VoWLlzY6jIq0am92Vf76dTeJqqviPhFI+NaGbs9QP2f/HOBzS2qRZJEa0PhJuDC4iqkU4HnMvNlp44kSROnstNHEfENYCFwcET0AH8OTAPIzJXAGuA9QDfwInBRVbVIkhpTWShk5vmj3J/AJVX9/5Kk5nXeS/mSpDEzFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklSqNBQiYnFEPBwR3RFx2TD3HxAR/xARP4mIByPioirrkSSNrLJQiIgpwJXAGcAxwPkRccyQYZcAD2Xm8cBC4NMRsU9VNUmSRlblkcLJQHdmPpqZu4DVwJIhYxKYFREBzAR+CfRXWJMkaQSRmdVsOOJsYHFmLi2WLwBOyczldWNmATcBbwRmAedm5j8Ns61lwDKAOXPmnLh69epKah6r3t5eZs6c2eoyKtGpvdlX++nU3iaqr0WLFt2bmQtGGze1whpimHVDE+g/AeuBtwNvAG6NiB9m5vN7PChzFbAKYMGCBblw4cLxr/YVWLt2LZOtpvHSqb3ZV/vp1N4mW19Vnj7qAQ6tW54LbB4y5iLg+qzpBh6jdtQgSWqBKkPhHmBeRBxRvHh8HrVTRfU2Au8AiIg5wNHAoxXWJEkaQWWnjzKzPyKWAzcDU4CrM/PBiLi4uH8l8JfA1yLifmqnm1Zk5jNV1SRJGlmVrymQmWuANUPWray7vRl4d5U1SJIa5zuaJUklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVKo0FCJicUQ8HBHdEXHZXsYsjIj1EfFgRPygynokSSOb2sigiHhTZj7QzIYjYgpwJfAuoAe4JyJuysyH6sYcCHwRWJyZGyPiV5v5PyRJ46vRI4WVEXF3RHyo+EXeiJOB7sx8NDN3AauBJUPGvB+4PjM3AmTmUw1uW5JUgcjMxgZGzAN+DzgHuBv4ambeOsL4s6kdASwtli8ATsnM5XVjPgtMA/4DMAv4XGZeM8y2lgHLAObMmXPi6tWrG+tugvT29jJz5sxWl1GJTu3NvtpPp/Y2UX0tWrTo3sxcMOrAzGz4C5gCnAVsAn4K/Az4z3sZew5wVd3yBcD/HjLmC8BdwP7AwcDPgd8YqYYTTzwxJ5vvf//7rS6hMp3am321n07tbaL6AtZlA7/nG31N4TjgIuC3gVuB92bmfRHxeuBO4PphHtYDHFq3PBfYPMyYZzLzBeCFiLgdOB54pJG6JEnjq9HXFL4A3Accn5mXZOZ9AJm5GfizvTzmHmBeRBwREfsA5wE3DRlzI/C2iJgaETOAU6gdgUiSWqChIwXgPcCOzNwNEBFdwH6Z+WJmXjvcAzKzPyKWAzdTO+10dWY+GBEXF/evzMyfRsT/BTYAA9RONzV1lZMkafw0Ggq3Ae8EeovlGcAtwFtGelBmrgHWDFm3csjyJ4FPNliHJKlCjZ4+2i8zBwOB4vaMakqSJLVKo6HwQkScMLgQEScCO6opSZLUKo2ePvow8K2IGLx66HXAudWUJElqlYZCITPviYg3AkcDAfwsM/sqrUySNOEaPVIAOAk4vHjMmyOCHObdx5Kk9tXom9euBd4ArAd2F6sTMBQkqYM0eqSwADimeKu0JKlDNXr10QPAr1VZiCSp9Ro9UjgYeCgi7gZ2Dq7MzDMrqUqS1BKNhsLlVRYhSZocGr0k9QcR8evAvMy8rZi8bkq1pUmSJlpDrylExO8D3wa+XKw6BLihqqIkSa3R6AvNlwCnA88DZObPAT9PWZI6TKOhsDNrn7MMQERMpfY+BUlSB2k0FH4QEX8KTI+IdwHfAv6hurIkSa3QaChcBjwN3A/8AbXPSNjbJ65JktpUo1cfDQBfKb4kSR2q0bmPHmOY1xAy88hxr0iS1DLNzH00aD/gHOBXxr8cSVIrNfSaQmZurfvalJmfBd5ecW2SpAnW6OmjE+oWu6gdOcyqpCJJUss0evro03W3+4HHgd8d92okSS3V6NVHi6ouRJLUeo2ePvofI92fmZ8Zn3IkSa3UzNVHJwE3FcvvBW4HnqiiKElSazTzITsnZOZ2gIi4HPhWZi6tqjBJ0sRrdJqLw4Bddcu7gMPHvRpJUks1eqRwLXB3RHyH2jub3wdcU1lVkqSWaPTqo7+KiH8G3lasuigzf1xdWZKkVmj09BHADOD5zPwc0BMRR1RUkySpRRr9OM4/B1YA/7NYNQ3426qKkiS1RqNHCu8DzgReAMjMzTjNhSR1nEZDYVdmJsX02RGxf3UlSZJapdFQ+GZEfBk4MCJ+H7gNP3BHkjpOo1cffar4bObngaOBj2bmrZVWJkmacKMeKUTElIi4LTNvzcw/ycw/bjQQImJxRDwcEd0RcdkI406KiN0RcXYzxUuSxteooZCZu4EXI+KAZjYcEVOAK4EzgGOA8yPimL2M+wRwczPblySNv0bf0fwScH9E3EpxBRJAZv7hCI85GejOzEcBImI1sAR4aMi4/w5cR23CPUlSCzUaCv9UfDXjEPacRbUHOKV+QEQcQu1y17czQihExDJgGcCcOXNYu3Ztk6VUq7e3d9LVNF46tTf7aj+d2ttk62vEUIiIwzJzY2Z+fQzbjmHW5ZDlzwIrMnN3xHDDiwdlrgJWASxYsCAXLlw4hnKqs3btWiZbTeOlU3uzr/bTqb1Ntr5Ge03hhsEbEXFdk9vuAQ6tW54LbB4yZgGwOiIeB84GvhgRv9Pk/yNJGiejnT6q//P9yCa3fQ8wr5gjaRNwHvD++gGZWc6fFBFfA/4xM29AktQSo4VC7uX2qDKzPyKWU7uqaApwdWY+GBEXF/evbKpSSVLlRguF4yPieWpHDNOL2xTLmZmvGenBmbkGWDNk3bBhkJn/taGKJUmVGTEUMnPKRBUiSWq9Zj5PQZLU4QwFSVLJUJAklQwFSVLpVRMKW3t38pMnnmVr785WlyKpSVt7d7Kjb7fP3wnwqgiFG9dv4vRPfI8PXPUjTv/E97hp/aZWlySpQYPP38eefsHn7wTo+FDY2ruTFddt4KW+Abbv7OelvgE+ct0G/+KQ2kD983d3ps/fCdDxodCzbQfTuvZsc1pXFz3bdrSoIkmN8vk78To+FOYeNJ2+gYE91vUNDDD3oOktqkhSo3z+TryOD4XZM/flirOOY79pXczadyr7TeviirOOY/bMfVtdmqRR1D9/p0T4/J0AjX7ITls7c/4hnH7UwfRs28Hcg6b7AyW1kcHn79133sG/nPlWn78Ve1WEAtT+4vCHSWpPs2fuy/RpU3wOT4COP30kSWqcoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqRSpaEQEYsj4uGI6I6Iy4a5/79ExIbi618j4vgq65EkjayyUIiIKcCVwBnAMcD5EXHMkGGPAf8xM48D/hJYVVU9kqTRVXmkcDLQnZmPZuYuYDWwpH5AZv5rZm4rFu8C5lZYjyRpFJGZ1Ww44mxgcWYuLZYvAE7JzOV7Gf/HwBsHxw+5bxmwDGDOnDknrl69upKax6q3t5eZM2e2uoxKdGpv9tV+OrW3iepr0aJF92bmgtHGTa2whhhm3bAJFBGLgP8GvHW4+zNzFcWppQULFuTChQvHqcTxsXbtWiZbTeOlU3uzr/bTqb1Ntr6qDIUe4NC65bnA5qGDIuI44CrgjMzcWmE9kqRRVPmawj3AvIg4IiL2Ac4DbqofEBGHAdcDF2TmIxXWIklqQGVHCpnZHxHLgZuBKcDVmflgRFxc3L8S+CgwG/hiRAD0N3LOS5JUjSpPH5GZa4A1Q9atrLu9FHjZC8uCrb076dm2g7kHTWf2zH3HbWw76dS+qtK9ZTvbXuyje8t2jpozq9XlqE1VGgoamxvXb2LFdRuY1tVF38AAV5x1HGfOP+QVj20nndpXVT56w/1cc9dG/ujYfi79m9u58LTD+NiSY1tdltqQ01xMMlt7d7Liug281DfA9p39vNQ3wEeu28DW3p2vaGw76dS+qtK9ZTvX3LVxj3XX3LmR7i3bW1SR2pmhMMn0bNvBtK49d8u0ri56tu14RWPbSaf2VZX1Tzzb1HppJIbCJDP3oOn0DQzssa5vYIC5B01/RWPbSaf2VZX5hx7Y1HppJIbCJDN75r5ccdZx7Deti1n7TmW/aV1ccdZxw77Q2szYdtKpfVXlqDmzuPC0w/ZYd+Fph/lis8bEF5onoTPnH8LpRx3c0JU3zYxtJ53aV1U+tuRYLjz1cO6/9y5uu/RUA0FjZihMUrNn7tvwL8JmxraTTu2rKkfNmUXPjGkGgl4RTx9JkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpVGkoRMTiiHg4Iroj4rJh7o+I+Hxx/4aIOKHKeqRmbe3dyU+eeJatvTtHHbvusa185paHWffY1nHbZjNju7dsZ9uLfXRv2T7q2GZUVW+zNezo2z3qdru3bOfb657o2O9BFdsdampVG46IKcCVwLuAHuCeiLgpMx+qG3YGMK/4OgX4UvGv1HI3rt/Eius2MK2ri76BAa446zjOnH/IsGM/cNVd3NFdC4PPf6+btx01m2uXnvqKttnM2I/ecD/X3LWRPzq2n0v/5nYuPO0wPrbk2DF2Xn29Y6nhD3+zj0s/8b29bnfwezCoE78H473d4VR5pHAy0J2Zj2bmLmA1sGTImCXANVlzF3BgRLyuwpqkhmzt3cmK6zbwUt8A23f281LfAB+5bsOwf6Wte2xrGQiDfti99WVHDM1ss5mx3Vu27/HLEOCaOze+4r+Wq6p3rDXsztzrdl8t34Px3O7eRGZWs+GIs4HFmbm0WL4AOCUzl9eN+UfgrzPzjmL5u8CKzFw3ZFvLgGXF4tHAw5UUPXYHA8+0uoiKdGpvI/YV06bPmHrQ634jurqmDK7LgYHd/duefCT7drxYP3bKrINfP2X/A1/2x8zuF559cvf2ZzaPZZvNjO2accDsqa957eEAu198jikzDgCg//mnHx948bmRz2WNoKp6x1rDYG/Dbbf+e1CvTb4H4/azOIpfz8zXjjaostNHQAyzbmgCNTKGzFwFrBqPoqoQEesyc0Gr66hCp/bWyX31P/dUx/UFndvbZPtZrPL0UQ9waN3yXGDzGMZIkiZIlaFwDzAvIo6IiH2A84Cbhoy5CbiwuArpVOC5zHyywpokSSOo7PRRZvZHxHLgZmAKcHVmPhgRFxf3rwTWAO8BuoEXgYuqqqdik/bU1jjo1N7sq/10am+Tqq/KXmiWJLUf39EsSSoZCpKkkqEwBhHxeETcHxHrI2Jdse7yiNhUrFsfEe9pdZ3NiogDI+LbEfGziPhpRJwWEb8SEbdGxM+Lfw9qdZ3N2ktfnbC/jq6rf31EPB8RH273fTZCX52wzy6NiAcj4oGI+EZE7DfZ9pevKYxBRDwOLMjMZ+rWXQ70ZuanWlXXKxURXwd+mJlXFVeMzQD+FPhlZv51MX/VQZm5oqWFNmkvfX2YNt9f9YppZTZRmybmEtp8nw0a0tdFtPE+i4hDgDuAYzJzR0R8k9rFNscwifaXRwoCICJeA/wW8H8AMnNXZj5LbSqSrxfDvg78TmsqHJsR+uo07wD+X2b+gjbfZ0PU99UJpgLTI2IqtT9ONjPJ9pehMDYJ3BIR9xZTcAxaXsz2enWrDwHH4EjgaeCrEfHjiLgqIvYH5gy+d6T491dbWeQY7K0vaO/9NdR5wDeK2+2+z+rV9wVtvM8ycxPwKWAj8CS192XdwiTbX4bC2JyemSdQm+X1koj4LWozvL4BmE9th3+6hfWNxVTgBOBLmflm4AXgZdOdt6G99dXu+6tUnBI7E/hWq2sZT8P01db7rAixJcARwOuB/SPiA62t6uUMhTHIzM3Fv08B3wFOzswtmbk7MweAr1CbJbad9AA9mfmjYvnb1H6Zbhmcubb496kW1TdWw/bVAfur3hnAfZm5pVhu9302aI++OmCfvRN4LDOfzsw+4HrgLUyy/WUoNCki9o+IWYO3gXcDD8SeU36/D3igFfWNVWb+G/BERBxdrHoH8BC1qUg+WKz7IHBjC8obs7311e77a4jz2fMUS1vvszp79NUB+2wjcGpEzIiIoPaz+FMm2f7y6qMmRcSR1I4OoHZq4u8y868i4lpqh7UJPA78QbvN4xQR84GrgH2AR6ld7dEFfBM4jNoP9TmZ+cuWFTkGe+nr87T5/gKIiBnAE8CRmflcsW427b/PhuurE55jfwGcC/QDPwaWAjOZRPvLUJAklTx9JEkqGQqSpJKhIEkqGQqSpJKhIEkqVfbJa9JEKy7F/G6x+GvAbmpTXEDtDYa7WlLYCCLi94A1xfsppJbzklR1pMk0a21ETMnM3Xu57w5geWaub2J7UzOzf9wKlOp4+kivChHxwYi4u5iH/4sR0RURUyPi2Yj4ZETcFxE3R8QpEfGDiHh0cL7+iFgaEd8p7n84Iv6swe1+PCLuBk6OiL+IiHuKefRXRs251N6M9ffF4/eJiJ6IOLDY9qkRcVtx++MR8eWIuJXa5H5TI+Izxf+9ISKWTvx3VZ3IUFDHi4g3UZsW4S2ZOZ/aadPzirsPAG4pJjjcBVxObfqBc4CP1W3m5OIxJwDvj4j5DWz3vsw8OTPvBD6XmScBxxb3Lc7MvwfWA+dm5vwGTm+9GXhvZl4ALAOeysyTgZOoTcx42Fi+P1I9X1PQq8E7qf3iXFebcobp1KZQANiRmbcWt++nNp1xf0TcDxxet42bM3MbQETcALyV2vNnb9vdxb9PhwLwjoj4E2A/4GDgXuCfm+zjxsx8qbj9buA3I6I+hOZRmyZBGjNDQa8GAVydmf9rj5W1Dzqp/+t8ANhZd7v++TH0xbccZbs7snjBrpjH5wvUZmfdFBEfpxYOw+nn34/gh455YUhPH8rM7yKNI08f6dXgNuB3I+JgqF2lNIZTLe+O2mc9z6A2J/6/NLHd6dRC5pliht2z6u7bDsyqW34cOLG4XT9uqJuBDxUBNPi5xtOb7El6GY8U1PEy8/5idsrbIqIL6AMupvZRiI26A/g7ah/ycu3g1UKNbDczt0btc6IfAH4B/Kju7q8CV0XEDmqvW1wOfCUi/g24e4R6vkxtVs31xamrp6iFlfSKeEmqNIriyp43ZeaHW12LVDVPH0mSSh4pSJJKHilIkkqGgiSpZChIkkqGgiSpZChIkkr/HzHofwgP0tIHAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "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": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<table class=\"simpletable\">\n",
+       "<caption>Generalized Linear Model Regression Results</caption>\n",
+       "<tr>\n",
+       "  <th>Dep. Variable:</th>      <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model Family:</th>       <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Link Function:</th>        <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -3.9210</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Date:</th>           <td>Sun, 26 Apr 2020</td> <th>  Deviance:          </th> <td>  3.0144</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Time:</th>               <td>00:12:48</td>     <th>  Pearson chi2:      </th>  <td>  5.00</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>No. Iterations:</th>         <td>6</td>        <th>  Covariance Type:   </th> <td>nonrobust</td>\n",
+       "</tr>\n",
+       "</table>\n",
+       "<table class=\"simpletable\">\n",
+       "<tr>\n",
+       "       <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Intercept</th>   <td>    5.0850</td> <td>    7.477</td> <td>    0.680</td> <td> 0.496</td> <td>   -9.570</td> <td>   19.740</td>\n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Temperature</th> <td>   -0.1156</td> <td>    0.115</td> <td>   -1.004</td> <td> 0.316</td> <td>   -0.341</td> <td>    0.110</td>\n",
+       "</tr>\n",
+       "</table>"
+      ],
+      "text/plain": [
+       "<class 'statsmodels.iolib.summary.Summary'>\n",
+       "\"\"\"\n",
+       "                 Generalized Linear Model Regression Results                  \n",
+       "==============================================================================\n",
+       "Dep. Variable:              Frequency   No. Observations:                   23\n",
+       "Model:                            GLM   Df Residuals:                       21\n",
+       "Model Family:                Binomial   Df Model:                            1\n",
+       "Link Function:                  logit   Scale:                          1.0000\n",
+       "Method:                          IRLS   Log-Likelihood:                -3.9210\n",
+       "Date:                Sun, 26 Apr 2020   Deviance:                       3.0144\n",
+       "Time:                        00:12:48   Pearson chi2:                     5.00\n",
+       "No. Iterations:                     6   Covariance Type:             nonrobust\n",
+       "===============================================================================\n",
+       "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
+       "-------------------------------------------------------------------------------\n",
+       "Intercept       5.0850      7.477      0.680      0.496      -9.570      19.740\n",
+       "Temperature    -0.1156      0.115     -1.004      0.316      -0.341       0.110\n",
+       "===============================================================================\n",
+       "\"\"\""
+      ]
+     },
+     "execution_count": 13,
+     "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": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<table class=\"simpletable\">\n",
+       "<caption>Generalized Linear Model Regression Results</caption>\n",
+       "<tr>\n",
+       "  <th>Dep. Variable:</th>      <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Model Family:</th>       <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Link Function:</th>        <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -23.526</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Date:</th>           <td>Sun, 26 Apr 2020</td> <th>  Deviance:          </th> <td>  18.086</td> \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Time:</th>               <td>00:13:21</td>     <th>  Pearson chi2:      </th>  <td>  30.0</td>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>No. Iterations:</th>         <td>6</td>        <th>  Covariance Type:   </th> <td>nonrobust</td>\n",
+       "</tr>\n",
+       "</table>\n",
+       "<table class=\"simpletable\">\n",
+       "<tr>\n",
+       "       <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Intercept</th>   <td>    5.0850</td> <td>    3.052</td> <td>    1.666</td> <td> 0.096</td> <td>   -0.898</td> <td>   11.068</td>\n",
+       "</tr>\n",
+       "<tr>\n",
+       "  <th>Temperature</th> <td>   -0.1156</td> <td>    0.047</td> <td>   -2.458</td> <td> 0.014</td> <td>   -0.208</td> <td>   -0.023</td>\n",
+       "</tr>\n",
+       "</table>"
+      ],
+      "text/plain": [
+       "<class 'statsmodels.iolib.summary.Summary'>\n",
+       "\"\"\"\n",
+       "                 Generalized Linear Model Regression Results                  \n",
+       "==============================================================================\n",
+       "Dep. Variable:              Frequency   No. Observations:                   23\n",
+       "Model:                            GLM   Df Residuals:                       21\n",
+       "Model Family:                Binomial   Df Model:                            1\n",
+       "Link Function:                  logit   Scale:                          1.0000\n",
+       "Method:                          IRLS   Log-Likelihood:                -23.526\n",
+       "Date:                Sun, 26 Apr 2020   Deviance:                       18.086\n",
+       "Time:                        00:13:21   Pearson chi2:                     30.0\n",
+       "No. Iterations:                     6   Covariance Type:             nonrobust\n",
+       "===============================================================================\n",
+       "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
+       "-------------------------------------------------------------------------------\n",
+       "Intercept       5.0850      3.052      1.666      0.096      -0.898      11.068\n",
+       "Temperature    -0.1156      0.047     -2.458      0.014      -0.208      -0.023\n",
+       "===============================================================================\n",
+       "\"\"\""
+      ]
+     },
+     "execution_count": 14,
+     "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": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3dyb7QmLYISA7yA5hEXEBrYK2KiriinVBpHWp7SNVn199tE+16oNt1VZxQ3GpgisupYJa44JbQBBkX8UEkJ0kkD33748ZMGAgQzLJLPm8rivXzDlzn3O+dwY+c3LmnPuYcw4REYkunlAXICIiwadwFxGJQgp3EZEopHAXEYlCCncRkSikcBcRiUI1hruZPW1mW83s28O8bmb2sJmtMbPFZjYw+GWKiMjRCGTPfTow+givnwl09f9MBKbWvSwREamLGsPdOfcxsPMITc4FnnM+XwDpZtY6WAWKiMjRiwnCOtoC31eZzvXP23xoQzObiG/vnsTExKx27drVaoOVlZV4PNHxdYH6Ep6ipS/R0g9QX/ZbtWrVdudc85raBSPcrZp51Y5p4Jx7AngCYNCgQW7+/Pm12mB2djYjRoyo1bLhRn0JT9HSl2jpB6gv+5nZd4G0C8bHYC5QdRc8E9gUhPWKiEgtBSPc3wKu8J81czywxzn3k0MyIiLScGo8LGNmLwEjgGZmlgvcCcQCOOceA2YDZwFrgH3AVfVVrIiIBKbGcHfOXVLD6w64PmgViUhEKCsrIzc3l+Li4gbZXlpaGsuXL2+QbdW3QPqSkJBAZmYmsbGxtdpGML5QFZFGKDc3l9TUVDp06IBZdedVBFdBQQGpqan1vp2GUFNfnHPs2LGD3NxcOnbsWKttRMd5RSLS4IqLi2natGmDBHtjY2Y0bdq0Tn8VKdxFpNYU7PWnrr9bhbuISBTSMXcRiVher5c+ffocmJ41axYdOnQIXUFhROEuIhErMTGRRYsWHfb18vJyYmIaZ8zpsIyIRJXp06dz4YUXcvbZZ3PGGWcAMGXKFAYPHkzfvn258847D7S955576N69Oz/72c+45JJLeOCBBwAYMWIE+4dH2b59+4G/BioqKpg8efKBdT3++OPAj8MJjB07lh49enDZZZfhO0sccnJyOOGEE+jXrx9DhgyhoKCAUaNGHfShNHz4cBYvXhzU30Pj/EgTkaD649tLWbYpP6jr7NmmCXee3euIbYqKiujfvz8AHTt25I033gDg888/Z/HixWRkZDB37lxWr17NV199hXOOc845h48//pjk5GRmzJjBwoULKS8vZ+DAgWRlZR1xe9OmTSMtLY2cnBxKSkoYPnz4gQ+QhQsXsnTpUtq0acPw4cOZN28eQ4YM4aKLLmLmzJkMHjyY/Px8EhMTueKKK5g+fToPPvggq1atoqSkhL59+wbht/YjhbuIRKzDHZY5/fTTycjIAGDu3LnMnTuXAQMGAFBYWMjq1aspKCjgvPPOIykpCYBzzjmnxu3NnTuXxYsX8+qrrwKwZ88eVq9eTVxcHEOGDCEzMxOA/v37s2HDBtLS0mjdujWDBw8GoEmTJgCcd955DB8+nClTpvD0009z5ZVX1u0XUQ2Fu4jUWU172A0tOTn5wHPnHLfffjvXXXfdQW0efPDBw55uGBMTQ2VlJcBB55o75/j73//OqFGjDmqfnZ1NfHz8gWmv10t5eTnOuWq3kZSUxOmnn86bb77Jyy+/TG1HyD0SHXMXkag2atQonn76aQoLCwHIy8tj69atnHzyybzxxhsUFRVRUFDA22+/fWCZDh06sGDBAoADe+n71zV16lTKysoAWLVqFXv37j3stnv06MGmTZvIyckBfFemlpeXAzBhwgRuuukmBg8efOCvjGDSnruIRLUzzjiD5cuXM2zYMABSUlJ44YUXGDhwIBdddBH9+/fn2GOP5aSTTjqwzC233MK4ceN4/vnnOfXUUw/MnzBhAhs2bGDgwIE452jevDmzZs067Lbj4uKYOXMmN954I0VFRSQmJvL+++8DkJWVRZMmTbjqqnoaa9E5F5KfrKwsV1sffvhhrZcNN+pLeIqWvtRnP5YtW1Zv665Ofn5+va7/zjvvdFOmTKnXbeyXn5/v8vLyXNeuXV1FRcVh21X3OwbmuwAyVodlREQa2IsvvsjQoUO555576u3WgTosIyIC3HXXXQ22rUsvvfQnX/AGm/bcRaTWnKv2dskSBHX93SrcRaRWEhIS2LFjhwK+Hjj/eO4JCQm1XocOy4hIrWRmZpKbm8u2bdsaZHvFxcV1CrtwEkhf9t+JqbYU7iJSK7GxsbW+S1BtZGdnH7jKNNI1RF90WEZEJAop3EVEopDCXUQkCincRUSikMJdRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSiUMSF+77Sct7bUEZ5RWWoSxERCVsRF+7vLN7MP1eUMu7xz/lux+HvXSgi0phFXLiPG9SOSX3jWbO1kDMf+oSZORs15KiIyCEiLtwBjm8Tw7s3n0z/dunc+toSbnhxIXuKykJdlohI2IjIcAdok57IC9cM5dbRPZizdAtnPfQJX2/cFeqyRETCQsSGO4DHY/xqRGde/dUJeDww7rHPefLjdTpMIyKNXkDhbmajzWylma0xs9uqeT3NzN42s2/MbKmZXRX8Ug+vf7t03rnxJE47rgX3zF7OxOcX6DCNiDRqNYa7mXmBR4AzgZ7AJWbW85Bm1wPLnHP9gBHAX8wsLsi1HlFaYiyPXZ7F//yiJx+u2Mq5//iUFVvyG7IEEZGwEcie+xBgjXNunXOuFJgBnHtIGwekmpkBKcBOoDyolQbAzLj6xI7MmHg8+0orGPPIPN76ZlNDlyEiEnJW0/FpMxsLjHbOTfBPjweGOuduqNImFXgL6AGkAhc55/5VzbomAhMBWrZsmTVjxoxaFV1YWEhKSsoR2+wuqeTRRSWs2lXJWR1jGdstFo9ZrbZXnwLpS6RQX8JPtPQD1Jf9Ro4cucA5N6jGhs65I/4AFwJPVZkeD/z9kDZjgb8BBnQB1gNNjrTerKwsV1sffvhhQO1Kyirc/3tjsTv21nfcFdO+dHuKSmu9zfoSaF8igfoSfqKlH86pL/sB810Nue2cC+iwTC7Qrsp0JnDosY6rgNf9217jD/ceAay7XsXFeLh7TB/uPb8P89Zs5/xHP9NVrSLSKAQS7jlAVzPr6P+S9GJ8h2Cq2gicBmBmLYHuwLpgFloXlwxpz/PXDGV7YQnnPjKPL9btCHVJIiL1qsZwd86VAzcAc4DlwMvOuaVmNsnMJvmb/Qk4wcyWAB8AtzrnttdX0bUxrHNT3rx+OE2T4xg/7UveWJgb6pJEROpNTCCNnHOzgdmHzHusyvNNwBnBLS34jm2azOu/Gs6kFxbw25nfsHFHETed1gULwy9aRUTqIqKvUK2NtKRYnr16COcPbMvf3l/Fba8t0fDBIhJ1AtpzjzZxMR7+cmE/MtMTefg/a9hWWMI/Lh1AUlyj/HWISBRqdHvu+5kZvzujO/ec15vslVu59Mkv2bW3NNRliYgERaMN9/0uG3osUy/PYtnmfC58/HM27S4KdUkiInXW6MMdYFSvVjx71RC27Clm7NTPWLetMNQliYjUicLdb1jnpsyYeDzF5ZWMe/xzlm/WoGMiErkU7lX0bpvGy9cNI8bj4aLHP2ehbv4hIhFK4X6ILi1SeGXSMNKT4hg/7StyNuwMdUkiIkdN4V6NdhlJvHzdMFo0ieeKaV/x2ZqwuthWRKRGCvfDaJWWwMyJw2ifkcRV03P4eNW2UJckIhIwhfsRNE+N56WJx9OxWTITnpuvgBeRiKFwr0FGchwvXns8nZunMOG5+XykgBeRCKBwD0BGchwvThhK5+YpTHxuPvN0DF5EwpzCPUDHJMfxzwlD6dA0mWuezdGY8CIS1hTuRyEjOY5/XjuUzGOSuHp6Dgu+02mSIhKeFO5HqVlKPC9eO5SWTRK48ukcvs3bE+qSRER+QuFeCy1SE/jnhKE0SYxl/LQvWbmlINQliYgcROFeS23SE3nx2qHExXi4fNqXuvG2iIQVhXsdHNs0mReuGUp5RSWXPfUlW/YUh7okERFA4V5nXVum8uzVQ9i9r4zLp+mGHyISHhTuQdA3M50nrxjExp37uHJ6DntLykNdkog0cgr3IBnWuSn/uGQAS3J3M+mFBZSW66bbIhI6CvcgOqNXK+47vy+frN7Of73yDZWVLtQliUgjFRPqAqLNuMHt2LG3lPvfXUHzlHju+MVxmFmoyxKRRkbhXg8mndKJrQXFPD1vPS2bxHPdKZ1DXZKINDIK93pgZtzx855sKyjh3n+voEWTeM4bkBnqskSkEVG41xOPx/jLuH7sKCzl968upkVqAsO7NAt1WSLSSOgL1XoUH+PlsfFZdGqWwqTnF7B8c36oSxKRRkLhXs/SEmN55qrBJMfHcNUzOWzeUxTqkkSkEVC4N4A26Yk8c9VgCkvKueqZHAqKy0JdkohEOYV7AzmudRMevWwgq7cWcv2LCymr0EVOIlJ/FO4N6ORuzfnzeb35eNU2/ufNpTini5xEpH7obJkGdtHg9ny3Yx+PZq+lY7MkuoW6IBGJStpzD4FbzujOz/u05t5/r2D+Fg0yJiLBF1C4m9loM1tpZmvM7LbDtBlhZovMbKmZfRTcMqPL/nPg+7dL54nFJSzO3R3qkkQkytQY7mbmBR4BzgR6ApeYWc9D2qQDjwLnOOd6ARfWQ61RJSHWyxPjB5EaZ0x4dr5OkRSRoApkz30IsMY5t845VwrMAM49pM2lwOvOuY0AzrmtwS0zOjVPjee3WQnsK63gmunzNQ68iASN1XTGhpmNBUY75yb4p8cDQ51zN1Rp8yAQC/QCUoGHnHPPVbOuicBEgJYtW2bNmDGjVkUXFhaSkpJSq2XDTWFhIeuKEvjbghL6t/By44B4PBE6imS0vS/R0Jdo6QeoL/uNHDlygXNuUE3tAjlbprqkOfQTIQbIAk4DEoHPzewL59yqgxZy7gngCYBBgwa5ESNGBLD5n8rOzqa2y4ab7OxsbvrFCJq0Wc9dby8jp6Q1t47uEeqyaiXa3pdo6Eu09APUl6MVSLjnAu2qTGcCm6pps905txfYa2YfA/2AVUhAfnlCB1ZvLWRq9lq6NE/hgiyNIikitRfIMfccoKuZdTSzOOBi4K1D2rwJnGRmMWaWBAwFlge31OhmZtx1Ti9O6NyU219fwoLvdoa6JBGJYDWGu3OuHLgBmIMvsF92zi01s0lmNsnfZjnwLrAY+Ap4yjn3bf2VHZ1ivR4evWwgrdMTuO75BeTt1hk0IlI7AZ3n7pyb7Zzr5pzr7Jy7xz/vMefcY1XaTHHO9XTO9XbOPVhfBUe79KQ4pv1yECVllUx4VmfQiEjt6ArVMNSlRSoPXzqAlVvy+a+XdaNtETl6CvcwNbJ7C/77rON4d+kWHv7P6lCXIyIRRgOHhbFrTuzIii0FPPj+arq3TOXMPq1DXZKIRAjtuYcxM+Oe83ozsH06v3v5G5Zt0m36RCQwCvcwt/8+rGmJsVz73Hx2FJaEuiQRiQAK9wjQIjWBJ67IYnthCb/+59e6i5OI1EjhHiH6ZqZz/wV9+XL9Tv749tJQlyMiYU5fqEaQMQPasnxzPo9/vI6erdO4dGj7UJckImFKe+4R5veje3BKt+bc+da3zN+gIQpEpHoK9wjj9RgPXzyAzGOSmPTCAjZpiAIRqYbCPQKlJcXy5BVZFJdVct3zCyguqwh1SSISZhTuEapLi1T+dlF/luTt4fbXl1DTTVdEpHFRuEew03u25Hend+ONhXlM+3R9qMsRkTCicI9wN4zswuherfjz7OV8snpbqMsRkTChcI9wHo/xwLh+dGmRwg0vLmTjjn2hLklEwoDCPQqkxMfw5BW+++Ve+5zGgBcRhXvUOLZpMv+4dACrtxZwyyvf6AtWkUZO4R5FTuranNvO7MG/v93CIx+uCXU5IhJCCvcoc+1JnRjTvw1/eW8V7y/7IdTliEiIKNyjjJlx3wV96dWmCTfPXMSarYWhLklEQkDhHoUSYr08Pn4Q8TEeJj43nz1FZaEuSUQamMI9SrVNT2Tq5Vls3LmP38xYSIVusi3SqCjco9iQjhnceU4vsldu44G5K0Ndjog0II3nHuUuH9qeZZvymZq9lp6tm3B2vzahLklEGoD23KOcmfHHc3ox6NhjmPzqN3ybtyfUJYlIA1C4NwJxMR6mXp7FMUlxTHxuPtt1k22RqKdwbySap8bz5BWD2LG3lF+9sIDSct1kWySaKdwbkd5t0/i/sX3J2bCLO99aqiEKRKKYvlBtZM7t35blmwt47KO19GydyvhhHUJdkojUA+25N0KTR3Xn1B4t+OPby/hs7fZQlyMi9UDh3gh5PcZDF/enQ7Nkrv/n1xoDXiQKKdwbqdSEWJ66YhCVDiY8l0NBsYYoEIkmCvdGrEOzZB69bCBrt+3l5hmLNESBSBRRuDdyw7s0486ze/LBiq1MmaMhCkSihc6WEcYffywrt/jOoOnWMoXzB2aGuiQRqaOA9tzNbLSZrTSzNWZ22xHaDTazCjMbG7wSpb6ZGXed04thnZpy22tLWPDdzlCXJCJ1VGO4m5kXeAQ4E+gJXGJmPQ/T7n5gTrCLlPoX6/Uw9fKBtElPYOJzC8jdpTNoRCJZIHvuQ4A1zrl1zrlSYAZwbjXtbgReA7YGsT5pQOlJcUy7cjClFZVMeHa+zqARiWBW0yXo/kMso51zE/zT44GhzrkbqrRpC7wInApMA95xzr1azbomAhMBWrZsmTVjxoxaFV1YWEhKSkqtlg034diXpdsr+MuCYno383LzwHg8ZgEtF459qa1o6Uu09APUl/1Gjhy5wDk3qKZ2gXyhWt3/7EM/ER4EbnXOVdgRgsA59wTwBMCgQYPciBEjAtj8T2VnZ1PbZcNNOPZlBJCW+R1/mPUtnxS24M6zewW0XDj2pbaipS/R0g9QX45WIOGeC7SrMp0JbDqkzSBghj/YmwFnmVm5c25WUKqUBnf58ceybttenp63nk7NkjUGjUiECSTcc4CuZtYRyAMuBi6t2sA513H/czObju+wjII9wv2/nx/Hdzv2cudbS8nMSGJk9xahLklEAlTjF6rOuXLgBnxnwSwHXnbOLTWzSWY2qb4LlNDxeoyHLxnAca2bcMM/v2bZpvxQlyQiAQroPHfn3GznXDfnXGfn3D3+eY855x6rpu2V1X2ZKpEpOT6Gab8cTGpCLNc8m8OWPcWhLklEAqDhB6RGrdISmHblIPKLyrh6eg6FJeWhLklEaqBwl4D0apPGI5cNZOUPBdzw4teUV+g2fSLhTOEuARvRvQV/Orc32Su3cceb3+o2fSJhTAOHyVG5dGh7cnft49HstbRNT+SGU7uGuiQRqYbCXY7a5FHd2bKnmAfmrqJVWiJjszSKpEi4UbjLUTMz7rugL1sLSrjttcU0T43nlG7Na7WuWQvzmDJnJZt2F9EmPZHJo7ozZkDbIFcs9UXvX/jSMXeplbgY3yiS3Vqm8qsXFvDN97uPeh2zFuZx++tLyNtdhAPydhdx++tLmLUwL/gFS9Dp/QtvCneptdSEWKZfPZimKXFcPT2HLXuP7gyaKXNWUlRWcdC8orIK3REqQuj9C28Kd6mTFqkJPHvVEBzwwPxifsgP/CKnTbuLjmq+hBe9f+FN4S511ql5CtOvGkxhqeOKaV+xZ19g48C3SU88qvkSXvT+hTeFuwRF38x0bhqYwPrte7n62Rz2ldZ8FevkUd1JjPUeNC8x1svkUd3rq0wJIr1/4U3hLkHTs6mXBy/uz8KNu5j0wteUlh/5GPyYAW259/w+tE1PxIC26Ynce34fnW0RIfT+hTedCilBdVaf1tx3fl9+/9pifjtzEQ9fMgCv5/A3cBkzoK3CIILp/QtfCncJunGD25FfXMbd/1pOcryX+87vi+cIAS8iwadwl3ox4aROFBSX89AHq0mM9XLXOb040i0YRSS4FO5Sb27+WVf2lZbz5CfrSYjzctvoHgp4kQaicJd6Y2b891nHsa+0gsc/Wkd8jJffnd4t1GWJNAoKd6lXZsafzu1NaXklD3+wmliPceNpGklSpL4p3KXeeTy+gcYqKh1/eW8VMV4PvxrROdRliUQ1hbs0CK/HmHJhP8orHfe/uwKH49cjuoS6LJGopXCXBuP1GH8d1w8z+L93V+IcXD9SAS9SHxTu0qBivB7+Oq4/hm9UwfIKx02nddFZNCJBpnCXBuf1GH8Z1x+vx8Pf3l9FSXkFk0d1V8CLBJHCXULC6zGmjO1LXIyHR7PXUlxWyR2/OE4BLxIkCncJGY/H+PN5vYmP8fD0vPXsKy3nnvP6HHEsGhEJjMJdQsrMuPPsnqQmxPD3/6yhsKScv47rT1yMBiwVqQuFu4ScmfFfZ3QnJT6Ge/+9goLicqZePpCkOP3zFKkt7R5J2LjulM7cd34fPlm9jUuf/JJde0tDXZJIxFK4S1i5eEh7Hr0si2Wb87nw8c/J0/04RWpF4S5hZ3TvVjx71RB+yC/mvEfmsWxTfqhLEok4CncJS8M6N+WVScPwmDHu8c/5ZPW2UJckElEU7hK2erRqwhvXn0DmMYlc9UwOM3M2hrokkYihcJew1jotkVcmDWNY56bc+toS7n93BZWVLtRliYQ9hbuEvdSEWJ6+cjCXDGnP1Oy1/OqfC9hXWh7qskTCWkDhbmajzWylma0xs9uqef0yM1vs//nMzPoFv1RpzGK9Hv58Xm/u+EVP3lv2A2Onfs4mnUkjclg1hruZeYFHgDOBnsAlZtbzkGbrgVOcc32BPwFPBLtQETPjmhM7Mu3KwXy/cx9n//1Tvlq/M9RliYSlQPbchwBrnHPrnHOlwAzg3KoNnHOfOed2+Se/ADKDW6bIj0Z2b8Eb1w8nLTGWS5/8guc/34BzOg4vUpXV9J/CzMYCo51zE/zT44GhzrkbDtP+FqDH/vaHvDYRmAjQsmXLrBkzZtSq6MLCQlJSUmq1bLhRX2pvX5nj8cUlfLOtguFtYriiVxzx3uAMOhYt70u09APUl/1Gjhy5wDk3qMaGzrkj/gAXAk9VmR4P/P0wbUcCy4GmNa03KyvL1daHH35Y62XDjfpSNxUVle5v7610HW57x43620du/bbCoKw3Wt6XaOmHc+rLfsB8V0O+OucCOiyTC7SrMp0JbDq0kZn1BZ4CznXO7QhgvSJ15vEYN/+sG89cOZjNe4o5+++f8q/Fm0NdlkjIBRLuOUBXM+toZnHAxcBbVRuYWXvgdWC8c25V8MsUObIR3Vvwr5tOpHOLFK5/8WvumPUtxWUVoS5LJGRqDHfnXDlwAzAH3yGXl51zS81skplN8jf7H6Ap8KiZLTKz+fVWschhZB6TxMvXDePakzry/BffMeaReaz+oSDUZYmEREADZjvnZgOzD5n3WJXnE4CffIEq0tBmL9nM7CVbAFj1QwFnPfwJY/q3Zd6a7WzeU0yb9EQmj+rOmAFtg77tWQvzmDJnJZt2F9XrdgLxh1lLeOnL77m5dxnX3D6bS4a24+4xfUJSi4SG7oYgUWPWwjxuf30JRf7DMZUOXIXjlQW5B9rk7S7i9teXAAQ1eA/ddn1tJxB/mLWEF774cRyeCucOTCvgGw8NPyBRY8qclQfCdb/qTvQtKqtgypyV9b7t+thOIF768vujmi/RSeEuUeNohiMI9k1ADrftUAyRUHGYa1cON1+ik8Jdokab9MSA23oM3lyUF7QrWw+37aOpKVi8Vv2FXIebL9FJ4S5RY/Ko7iTGeg+aF+sxYg+5ajU+xkPmMUn8ZsYiJjw7n9xd++pl24mxXiaP6l7ndR+tS4a2O6r5Ep0U7hI1xgxoy73n96FteiIGtE1PZMqF/Zgytt9B8+6/oC8f3jKCP/z8OD5bu4PT//oxT3y8lrKKyqBu+97z+4TkbJm7x/Th8uPbH9hT95px+fHt9WVqI6OzZSSqjBnQttpArW7ehJM6Mbp3K+56ayl/nr2C1xbk8cdze3F8p6ZB3XYo3D2mD3eP6UN2djZrLxsR6nIkBLTnLo1a5jFJPPXLwTwxPovCknIufuILbnppITuLa78XLxIOtOcuApzRqxUndW3O1I/W8thHa3m3spLvvKu47pROJMXpv4lEHu25i/glxnn53end+OB3p9C/hZeHPljNyAeymZmzkQrdt1UijMJd5BDtMpL4df8EXpk0jNZpidz62hLOfOhj3lv2g24KIhFD4S5yGIM7ZPDGr0/g0csGUlbhuPa5+Yx59DM+Xb1dIS9hT+EucgRmxll9WvPeb0/m/gv6sC2/mMunfcm4xz9n3hqFvIQvhbtIAGK8Hi4a3J7/3DKC/z23F9/vLOKyp77kgqmf8Z8VOlwj4UfhLnIUEmK9XDGsA9mTR/Cnc3vxQ34JV0+fz1kPf8obC3PrdCGUSDAp3EVqISHWy3h/yD9wYT/KKir57cxvOPn/PuTxj9ayZ19ZqEuURk4n8IrUQazXw9isTM4f0JbsVVt54uN13PvvFTz4/mouyGrLFcM60K1laqjLlEZI4S4SBB6PcWqPlpzaoyVLN+1h+rwNvDw/lxe+2MjQjhlcfvyxnNGrJfEx3ppXJhIECneRIOvVJo0pF/bj9rOO4+X53/PCF99x40sLyUiO44KBbRk3qB1dtTcv9UzhLlJPMpLjmHRKZyae1IlP12znpa828sy8DTz5yXr6ZaYxNiuTn/dtQ0ZyXKhLlSikcBepZx6PcXK35pzcrTnbC0uYtTCPVxfkcsebS/nj28sY0b05Z/drw8+Oa0lyvP5LSnDoX5JIA2qWEs+Ekzox4aROLN+cz6yFeby5aBPvL99KQqyHU3u04MzerRnZowUpCnqpA/3rEQmR41o34bjWTbh1dA/mf7eLt7/ZxL+/3cLsJVuIi/FwUpdmnN6zJacd15LmqfGhLlcijMJdJMQ8HmNIxwyGdMzgrnN6seC7Xcxespn3lv3AByu2YraEvpnpnNajBSO6N6d3mzQ8Ht0PVY5M4S4SRrxVgv7Os3uyfHMBHyz3hfzf3l/FX99bRUZyHCd2acaJXZtxYpc9/m12AAANA0lEQVRmIbkJt4Q/hbtImDIzerZpQs82TbjxtK5sLyzh09XbyV65lU/XbOetbzYB0LFZMsd3yuD4Tk0Z0jGD1mkKe1G4i0SMZinxB+7T6pxj5Q8FfLp6O5+v3cE732zmpa++B6BdRiKDj81g4LHH4Aoqqah0eHUYp9FRuItEIDOjR6sm9GjVhAkndaK8opLlmwv4asNOvlq/g49Xb+P1hXkA3Jczh76Z6fRrl06/zDT6tkunTVoCZgr8aKZwF4kCMV4PfTLT6JOZxjUndsQ5x8ad+3jh3c8oSWnNwo27mfbpOsoqfEMTH5MUS++2afRs7Tvsc1zrJnRslkysV2MJRguFu0gUMjOObZrM8LaxjBjRG4CS8gpWbC5gce5ulm7K59tNe3hm3gZK/cMUx3k9dGqeTPdWqXRtkULXlql0aZFC+4wkhX4EUriLNBLxMV7foZl26QfmlVVUsm7bXpZvzmfFlgJWbMln/oZdvLlo04E2sV6jfUYSHZul0Kl5Msc2TaJj02TaN02idVqijueHKYW7SCMW6/XQvVUq3VsdPJBZYUk5a7YWsnZrIWu2FbJ+217Wb9/Lx6u3UVpeWWV5o216Iu0yksg8JonMYxJpm55I22MSaZOeSMvUeGK01x8SCncR+YmU+Bj6t0unf5W9fIDKSseW/GI2bN/Lhh37+H7XPjbu2Efurn3M3bSFHXtLD2rvMWiRmkCrtARaNfE9tmgST4vUBFqkxtOiSTzNUuLJSIrThVlBpnAXkYB5PEabdN9e+Qldfvr6vtJyNu0uIndXEZv3FLN5dxGb9hTzQ34xa7YVMm/tdgqKy3+ynNdjZCTH0TQ5jmYp8TRNiSMjOY6MpDgyUuI4JimO73ZU0HJzPsckxZGWGEtCrEdn/ByBwl1EgiYpLoYuLVLp0uLw49XvKy1na34JWwtK2FpQzPaCErYXlrKtoIQde0vZXljCxp372LW3lIKSgz8I7s/55MDzOK+HtKRYmiTEkJYYS5PEWJokxJKaEEPqgccYUuJ//En2//iee0mKi4na7wwCCnczGw08BHiBp5xz9x3yuvlfPwvYB1zpnPs6yLWKRK1ZC/OYMmclm3YX0SY9kcmjuvPK/I3MW7vzQJvhnTO4cFD7n7QDfjJv/nc7eenL77m5dxnX3D6bS4a24+4xfQLa7pgBbQ87P5Dl92+7wjm8Zj/ZdlJcDB2axbDo+9019uXOs3tycrfm7NxXyofzcujQrSe79pXx2drtZK/cxraCEgqKy/B6jLIKx4bte8kvLqeguOzAaZ81iY/xkBTnC/qkOC+JcV4SY70kxXlJiPU9T4jzkhDjJSHWQ0Lsj4/xMR7iY/yPsT8+j9v/4/3xebzXS2yMEev14FxgtdVFjeFuZl7gEeB0IBfIMbO3nHPLqjQ7E+jq/xkKTPU/ikgNZi3M4/bXl1BUVgFA3u4ibp656Cft5q3deVDY5+0uYvKr34CDskp3YN7vZi6isspyFc7xwhcbAQ4K2eq2e/vrS5j/3U5eW5D3k/nAQQFf3fJ12fbkV74B40Ao5+0u4o43l3Lv+X0YM6AtW5p6GdGnNbMW5vHB8q0Hli0uq+T7nUUH2u1XXFZBYUk5hcXlFBSXU1hSzt6ScvaWlrO3pIK9JeXsK61gb2k5+0p9z4tKKw487thbemC6pLyC4rJKisoqqKisezCf2TGWkSPrvJojCmTPfQiwxjm3DsDMZgDnAlXD/VzgOef7OPrCzNLNrLVzbnPQKxaJMlPmrDwQVEerur3TymraAbz05fcHBWx12y0qqziw133o/ClzVh4UntUtX5dtl1UTmoFut7p2vj1sL81SgjtccllFJcVlFZSUV1JS7nte6n/ue6ygpKyS0grf9IHH8krKKnw/nl0bg1pTdaymPw/MbCww2jk3wT89HhjqnLuhSpt3gPucc5/6pz8AbnXOzT9kXROBif7J7sDKWtbdDNhey2XDjfoSnhqsL3GtumTV17or9u3Bm5R2YLp0y5oFddluXZYPwrLNgO1HWrbqNsJcXf59Heuca15To0D23Kv7tuHQT4RA2uCcewJ4IoBtHrkgs/nOuUF1XU84UF/CU7T0xczml+/ZGvH9gOh5T6Bh+hLI1QW5QLsq05nAplq0ERGRBhJIuOcAXc2so5nFARcDbx3S5i3gCvM5Htij4+0iIqFT42EZ51y5md0AzMF3KuTTzrmlZjbJ//pjwGx8p0GuwXcq5FX1VzIQhEM7YUR9CU/R0pdo6QeoL0elxi9URUQk8mhEHxGRKKRwFxGJQmEf7maWYGZfmdk3ZrbUzP7on59hZu+Z2Wr/4zGhrjUQZuY1s4X+awMiuR8bzGyJmS0ys/n+eZHal3Qze9XMVpjZcjMbFol9MbPu/vdj/0++md0coX35rf//+7dm9pI/ByKuHwBm9ht/P5aa2c3+efXel7APd6AEONU51w/oD4z2n5FzG/CBc64r8IF/OhL8BlheZTpS+wEw0jnXv8r5upHal4eAd51zPYB++N6fiOuLc26l//3oD2ThO7nhDSKsL2bWFrgJGOSc643vRI6LibB+AJhZb+BafFf69wN+YWZdaYi+OOci5gdIAr7GN27NSqC1f35rYGWo6wug/kz/G3kq8I5/XsT1w1/rBqDZIfMiri9AE2A9/pMLIrkvh9R/BjAvEvsCtAW+BzLwndH3jr8/EdUPf50X4htscf/0HcDvG6IvkbDnvv9QxiJgK/Cec+5LoKXzn0vvf2wRyhoD9CC+N7bqEByR2A/wXYE818wW+IeVgMjsSydgG/CM/3DZU2aWTGT2paqLgZf8zyOqL865POABYCOwGd91M3OJsH74fQucbGZNzSwJ3ynj7WiAvkREuDvnKpzvT81MYIj/T52IYma/ALY65yJl7IuaDHfODcQ3Iuj1ZnZyqAuqpRhgIDDVOTcA2EsE/Ll/JP6LDc8BXgl1LbXhP/58LtARaAMkm9nloa2qdpxzy4H7gfeAd4FvgJ/eraQeRES47+ec2w1kA6OBH8ysNYD/cWsISwvEcOAcM9sAzABONbMXiLx+AOCc2+R/3IrvuO4QIrMvuUCu/69BgFfxhX0k9mW/M4GvnXM/+KcjrS8/A9Y757Y558qA14ETiLx+AOCcm+acG+icOxnYCaymAfoS9uFuZs3NLN3/PBHfG78C35AHv/Q3+yXwZmgqDIxz7nbnXKZzrgO+P5n/45y7nAjrB4CZJZtZ6v7n+I6HfksE9sU5twX43sy6+2edhm8464jrSxWX8OMhGYi8vmwEjjezJDMzfO/JciKvHwCYWQv/Y3vgfHzvTb33JeyvUDWzvsCz+L4x9wAvO+f+18yaAi8D7fH9Y7jQObfz8GsKH2Y2ArjFOfeLSOyHmXXCt7cOvsMaLzrn7onEvgCYWX/gKSAOWIdv+AwPkdmXJHxfRnZyzu3xz4u498V/yvNF+A5hLAQmAClEWD8AzOwToClQBvzOOfdBQ7wnYR/uIiJy9ML+sIyIiBw9hbuISBRSuIuIRCGFu4hIFFK4i4hEoUBukC3SoPyniX3gn2wFVOAbIgBgiHOuNCSFHYGZXQ3M9p83LxJyOhVSwpqZ3QUUOuceCINavM65isO89ilwg3Nu0VGsL8Y51yCXokvjo8MyElHM7JfmG99/kZk9amYeM4sxs91mNsXMvjazOWY21Mw+MrN1ZnaWf9kJZvaG//WVZvaHANd7t5l9hW9coz+aWY5/fO7HzOcifMNRz/QvH2dmuVWurD7ezN73P7/bzB43s/fwDVYWY2Z/9W97sZlNaPjfqkQjhbtEDP+AcecBJ/gHkovBN5QDQBow1z+YWSlwF77L1i8E/rfKaob4lxkIXGpm/QNY79fOuSHOuc+Bh5xzg4E+/tdGO+dmAouAi5xvPPWaDhsNAM52zo0HJuIbUG4IMBjfIGzta/P7EalKx9wlkvwMXwDO9w05QiK+S+0Bipxz7/mfL8E3TGy5mS0BOlRZxxzn3C4AM5sFnIjv/8Hh1lvKj0MtAJxmZpOBBKAZsAD491H2403nXLH/+RnAcWZW9cOkK75L0kVqTeEukcSAp51zdxw00ywGXwjvV4nvDl77n1f9d37ol0yuhvUWOf8XU/5xW/4BDHTO5ZnZ3fhCvjrl/PiX8aFt9h7Sp1875z5AJIh0WEYiyfvAODNrBr6zampxCOMM890zNQnfmOHzjmK9ifg+LLb7R8W8oMprBUBqlekN+G51xyHtDjUH+LX/g2T/fVATj7JPIj+hPXeJGM65Jf7RAt83Mw++UfYmAZuOYjWfAi8CnYHn95/dEsh6nXM7zOxZfMMbfwd8WeXlZ4CnzKwI33H9u4AnzWwL8NUR6nkc38iAi/yHhLbi+9ARqROdCimNhv9MlN7OuZtDXYtIfdNhGRGRKKQ9dxGRKKQ9dxGRKKRwFxGJQgp3EZEopHAXEYlCCncRkSj0/wHRUJwHFwSFegAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "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": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9+P/XmS2ZTBKykAVICIthFcGCKCqiKIvsilSBqlWxaqv9ttZ71Vvrrfaq9f567a1tbyta61pawYUlilZQcEOhLpFA2ANhySRkz2S2s/z+mGQgEGASMpkl7+fjkUfmnDlz5vMhzHnP+Szvj2IYhoEQQghxAlOkCyCEECI6SYAQQgjRLgkQQggh2iUBQgghRLskQAghhGiXBAghhBDtCluAePDBB5kwYQKzZs1q93nDMPiv//ovpkyZwuzZsykpKQlXUYQQQnRC2ALEtddey3PPPXfK5zdu3EhZWRnvvfcev/rVr/jlL38ZrqIIIYTohLAFiAsuuIBevXqd8vl169Yxb948FEVhzJgxNDQ0UFlZGa7iCCGE6KCI9UE4nU5yc3OD27m5uTidzjO+zutT8fo1fH4Nv6rhV3VUTUfTdDTdQNcNZHK4EEKcPUuk3ri9i7iiKGd8XX2TD2dVY0jvoSigtDwwtZxfUU7+bVIUTC3bJkXBZGp5bApsd6esrBSqQqxfLJL6xa54rhv0jPp1VMQCRG5uLhUVFcHtiooKsrOzu/Q9DAOMlgd6YE+Hz6EAiknB3BIwzKZjv80mBbNZwWySwWBCiPgTsQAxefJkXnnlFWbOnMk333xDSkpKlweIrmAAhm6gY4DW/jEKtAQLE2azgsVkwmJWsFhM3X4HIoQQXSVsAeLee+/liy++oLa2lssuu4x77rkHVVUBWLhwIZMmTWLDhg1MmTIFu93O448/Hq6ihJ0BqLqBqmvgb/uc2aRgMZuwWo79SNAQQsQCJdbSfVfWNIfcBxGNFAgEDKuJBIsZm9XUpu+lJ7SDSv1iUzzXDXpG/ToqYk1MPZUB+DUdv6bTjIoC2KxmEqwmEm3y5xBCRA+5IkWYAXj9Gl6/RkOzH3OClWaPSmKCWZqihBARJQEiygQChY9GNyTaLCQlWLBaZJSUEKL7SYCIUoYBbq+K26tis5hw2K0kWM2RLpYQogeRABEDfKqOr9GL1Wwi2W4lwSaBQggRftJ2EUP8mk5tk5faRi+qpke6OEKIOCd3EDHI69fw1WvYEywkJ1mlM1sIERZyBxGjDKDZq3K03oPbq0a6OEKIOCQBIsbpukG9y0dtoxdNl2YnIUTXkQARJ7x+jeoGLz7/KRJGCSFEB0mAiCO6blDb6KXZ4z/zwUIIcQYSIOKMATQ0+6lv8srCSUKIsyIBIk65fRp1TT4JEkKITpMAEce8fo3aRi+6BAkhRCdIgIhzPlWnpsGDrkuQEEJ0jASIHkDVDGoaJUgIITpGAkQPoWqGNDcJITpEAkQP4td06iRICCFCJAGih/GpOvVNvkgXQwgRAyRA9EBev0aDS4KEEOL0JED0UM1eVWZcCyFOSwJED9bQ7Mfrk9xNQoj2SYDo4epcsviQEKJ9EiB6OMOAeknJIYRohwQIgV/TaXRLf4QQoi0JEAKAZo8q/RFCiDYkQIigepesSieEOEYChAjSDWhwSVOTECJAAoRow+vXcHvVSBdDCBEFYi5APPLcJr7edVRG3YRRY7NPMr8KIbBEugAddcDZyAFnI1/tqmLexIGkpyRGukhxRzeg3uUjPSUh0kURQkRQzN1BWMwKALsO1vO/y4v55Nsjkp00DKSpSQgRcwHil0smMCA3BQC/qlP02X6eXb2No/XuCJcs/jQ2+2RUkxA9WMwFiD69HSyZPYJ5EweSYDUDsL+ikd+v+FbuJrqY3jLLWgjRM8VcgAAwKQrjh+fw/xacR2FeLyAwG7jos/08X7SduiZvhEsYP3yqTpPMshaiRwprgNi4cSPTpk1jypQpLF269KTnGxsbufPOO5kzZw4zZ87k9ddf79D505IT+P7Vw7j2skHBu4m9hxt4ekUxX++WkU5dpcntx+eXWdZC9DRhCxCapvHoo4/y3HPPUVRUxJo1a9i9e3ebY1599VUGDx7MqlWrePnll3nyySfx+TrWpKEoCuOGZfPj685jQJ9A34THp/Ha+t289sFuPD7paO0K9S6fNN8J0cOELUAUFxdTUFBAfn4+NpuNmTNnsm7dujbHKIqCy+XCMAxcLhe9evXCYuncyNv0lASWzBzB1Rf2x2wKjHT6Znc1f3j9W8orm866Pj2dphu4pKlJiB4lbPMgnE4nubm5we2cnByKi4vbHLN48WLuuusuJk6ciMvl4re//S0m05ljVkaG45TPzb2ikLEjcnlu1VYqqpupafTyzKoS5k0azFXj+2NSlM5Xqpucrn6RpABp6UlYLWf3vSIrK6VrChSl4rl+8Vw3iP/6dVTYAkR77f/KCRfnjz/+mOHDh/PSSy9x4MABbrnlFsaNG0dycvJpz11T4zrt80lWE3fOGUnRZ/vZXFqJrhu88cFutu05ynWXn0NSYvTOD8zIcJyxfpHU2OAmI7XzkxOzslKoqmrswhJFl3iuXzzXDXpG/ToqbE1Mubm5VFRUBLedTifZ2dltjnnjjTeYOnUqiqJQUFBAXl4ee/fu7ZL3t1nNXHPZIBZeVRjswC49UMcf3iimvDJ+/xOEm0/VZQKdED1E2ALEqFGjKCsro7y8HJ/PR1FREZMnT25zTJ8+ffjss88AOHr0KPv27SMvL69ryzEok7vnj6JvZhIAdU0+lq7axqaSChnl1EmNbr90WAvRA4StrcVisfDwww+zZMkSNE1j/vz5FBYWsmzZMgAWLlzID3/4Qx588EFmz56NYRjcd999ZGRkdHlZMlMTuWPuuby9aT+fb3Oi6QarPimjvLKJeRMHnXWbek+j6wbNHpVkuzXSRRFChJFixNjX6MqaZpxn0U749a6jvLlxL34tkEKib2YSi6cOiZqkf9HeB9HKpEDvNHuHO/17QjtvvNYvnusGPaN+HdXjvjqPKezNnfNGktGSqfRwdTN/fGMrew/XR7hksUU3kL4IIeJcjwsQAH0yHfzo2lEMyU8DoNmr8nxRKZu2VZzhleJ4Lo8q/ThCxLEeGSAA7AkWbpo2lMtG9wVANwxWfVzGWx/tlQymIdJ1A7dXUnAIEa96bIAAMJkUpl/Yn+9OPie4zsQX2yt54Z1SaT4Jkcsjs6uFiFc9OkC0GnNOb34wZySpSYFROXsONfCnt7ZS3eCJcMmin6YbEkyFiFMSIFrkZSVz1zXH5kscrffwpze3UlbREOGSRT/J0SREfJIAcZxeDhs/mDOSEQPSgdbO6+18u7c6wiWLbqrcRQgRlyRAnMBmNbNoyhAmntcHAFUz+Pv7u/i4+IiM2DkN6YsQIv5IgGiHSVG4+qICZl8yAEUBA3h7036KPtsvKSZOQdUMWXtDiDgjAeI0JozMZfGUIVjNgX+mT7dW8Nr63aiaDINtjyxNKkR8kQBxBiMGZLBk9nCSEgJpq4r3VPPS2h14fTL+/0RyFyFEfJEAEYL87BTumDuStGQbALsP1fPsmm3yjbkd8m8iRPyQABGirDQ7d8w9l5x0OwCHj7pYuqqEuiZvhEsWXVTNkLsrIeKEBIgOaB0GW5AbyIp4tN7DMytLOFrnjnDJoouMaBIiPkiA6CB7goVbZgwLJvqrd/l4ZlUJh49Gf4ru7uJTdXx+uYsQItZJgOgEm8XM96YOYdSgTCCQ1fS5NdtkKdPjuDzSWS1ErJMA0UkWs4nrJ5/DBcMC62x7fBp/KdrO3sOSmgPA69fwqzIcWIhYJgHiLJhMCvMmDuSSc3MB8Pl1XnynlJ3ldREuWXSQvgghYpsEiLOkKAozJhRw+fn9APBrOi+/u4PSA7URLlnkeXyaTCoUIoZJgOgCiqIw9YJ8pl6QDwRSYL/63k62l9VEuGSR1yx9EULELAkQXejy8/tx9YX9gZYg8c9dlOzr2UHC7ZNlSYWIVRIgutjE0X2ZOaEACCxjuuz9XWztwenCDSPQ1CSEiD0SIMLgklF9mH3xACAQJP6+bnePvpOQtSKEiE0SIMJkwrm5zLlkAHDsTmJbD+2T8Km6dFYLEYMkQITRRSNz29xJLHt/V4/tuG6WuwghYo4EiDCbcG4usy4O9ElousHf3t/VI+dJeLzSWS1ErJEA0Q0uPrdPsONa0w1eeW8Huw/VR7hU3UuXzmohYo4EiG5yyag+TBsfmCehagYvv7uDfUd6VloO6awWIrZIgOhGk8b048qxeQD4VZ0X15b2qAR/PlVH06WZSYhYIQGim03+Tj8mjekLBHI3/fXtUo5U95xU4R65ixAiZkiA6GataTlaE/x5fBrPF22nsocsOiRrVgsROyRAREBrgr9xLanCXR6V54u2U9PgiXDJws/r09BlNJMQMUECRIQoisK8Swdy3uDAokMNLh/PF22nPs7XuDZA1qwWIkZIgIggk0lhwRWDGV6QDkBNo5ff/eOruM+AKsNdhYgNIQWI2267jQ8++KDDE502btzItGnTmDJlCkuXLm33mM8//5y5c+cyc+ZMvve973Xo/PHAbDJxw5WFDOqbCsDhKhcvvLM9rr9l+/zSzCRELAgpQFx//fW8+OKLXHXVVSxdupTa2jMvhqNpGo8++ijPPfccRUVFrFmzht27d7c5pqGhgUceeYQ//elPFBUV8bvf/a5ztYhxVouJG6cOJT87GYCDVS5efm9H3C7ZKc1MQsSGkALE1KlTeeGFF3j22WeprKxk1qxZ/Pu//ztbt2495WuKi4spKCggPz8fm83GzJkzWbduXZtjVq9ezZQpU+jbNzDsMzMz8yyqEtsSbGZunj6MvlkOAPYebuC19bvR43TegNcvAUKIaGfpzIusVisJCQncf//9TJw4kQceeOCkY5xOJ7m5ucHtnJwciouL2xxTVlaGqqrceOONuFwubrrpJubNm3fG98/IcHSm2FEvA/h/15/P//fKvzha56akrIZ3NpfzvenDUBQl0sXrMhkZDhQFemc64qperbKyUiJdhLCJ57pB/Nevo0IKEO+99x6vvPIK1dXVLFq0iKKiIhwOB6qqMnXq1HYDRHv9FSdeDDRNo6SkhBdeeAGPx8MNN9zA6NGjGThw4GnLU1MTvxPLMjIc3Dx9KM+sLKHJ7eeTbw5jBqa3rFQX6zIyHMG/n9/tJymxU99RolZWVgpVVfE5Oz6e6wY9o34dFdKnc8WKFdx+++1MnDix7YstFh566KF2X5Obm0tFRUVw2+l0kp2dfdIx6enpJCUlkZSUxLhx4ygtLT1jgIh3mamJ3DJjGM+u3obHp7Hxm8M47BYmntc30kXrUk0eP4kJZkxxeBchRDwIqQ/imWeeOSk4tJo8eXK7+0eNGkVZWRnl5eX4fD6KiopOOvbKK69ky5YtqKqK2+2muLiYwYMHd7AK8alPpoObpg/FYg5cPN/ZdIAvd1ZFuFRdS9eNuB/SK0QsCylALFq0iPr6Y+mp6+rqWLx48WlfY7FYePjhh1myZAkzZszg6quvprCwkGXLlrFs2TIABg8ezMSJE5kzZw4LFizguuuuY8iQIWdRnfgyIDeVhVcNwdTyBfuNDXvYceDMI8hiicvtR9Pjc7SWELFOMUKY3DB37lxWrlx5xn3dobKmGWcctxMe30bf6l87Knl9w14ArGYTt80aTv+c2OxMa69+9gQLvRy2CJWoa8VzO3Y81w16Rv06KqQ7CF3XaW5uDm67XC40TYYpdpexQ7OZPj7QSe3XdF5cu4PK2vhJ7uf2qnE750OIWBZSgJg1axa33norK1euZOXKldx2223MmTMn3GUTx5k4ug+XjuoDBC6of307vvI2Nbn9kS6CEOIEIY1iuuOOO8jOzmb9+vUYhsENN9wQ0nwF0XUURWH6Rf1pcvv5evdR6l0+/vpOKXfMGYk9IfaHinr9GqqmYzFLejAhokXIV5ZrrrmGa665JpxlEWdgUhSunTQIl8fProP1VNa6eendHdw6YzhWS+xfWJu9KqlJ8dEXIUQ8CClAVFdX8/LLL1NeXo6qHhuW2FNzJ0WSxWxi0VVDeG7NNg4ddbG/opF/rN/FoquGYDLF9nwCj1clxW6Ny9nVQsSikALEPffcw+DBg5kwYQJmszncZRJnkGAzc/PVw/jzyq3UNHjZVlbL6k/LmHPJgJi+uOoGuL1a3M2uFiJWhfRJbGho4Fe/+lW4yyI6INlu5ZYZw/nzW1txeVQ+3+akl8PG5ef3i3TRzkqzN/7SbwgRq0JquC4sLMTpdIa7LKKDMlMTufnqYdha+h/e21zOv3ZURrhUZ0fVDHyS6VWIqBDyHcScOXM4//zzSUhICO6XPojIy8tKZtGUIby0dge6YfDmxr0k260M7Z8e6aJ1msujYrNKU6YQkRZSgJg1axazZs0Kd1lEJw3JT+PaSYNY8eEedAOWvb+LJbNHkJeVHOmidYrXr6HpOmZT7I/MEiKWhRQgZHhr9PvOkCwaXD7e21yOTw3Mtr5z7kgyUxMjXbROcXs1ku0SIISIpJA+gWVlZSxcuDCYjbWkpITf//73YS2Y6LhJY/py0YgcIJAE74W3S2N2hnKzV7K8ChFpIQWIX/7yl9x1112kpASSPQ0fPpy1a9eGtWCi4xRFYdbFAxgxIND/UN3g4aW1pTHZ6avrBh6fBAkhIimkANHY2Mhll10WHGNvMpmwWq1hLZjoHJNJ4frJhRS0ZHs9WOVi2bpdaDG4trXbG3uBTYh4ElKAMJvN+P3+YIBwOp2YpAMxalktJm6cNoTevQL9DzsO1LHq433tLgMbzVrzMwkhIiPkBYPuvvtuamtr+f3vf8+iRYu49dZbw102cRaSEq3cMmMYKfbAnd7m0krWf3kowqXqOLf0RQgRMSGNYpo3bx55eXl88MEHuN1unnzyScaNGxfusomzlJ4SmEi3dHUJPr/Oun8dpJfDxrhh2Wd+cZRw+zSS7UZMpxARIlaFnNNg3LhxEhRiUN/eDhZPGcKL7wQm0r310V5SkmJnIp2uG3j9Gok2Sb8hRHcL6VM3f/78dr/BrVixossLJLpeYV4a8ycNYnmMTqRzeyVACBEJIX3q7r///uBjr9dLUVER2dmx00wh4PwhWdTH6EQ6mVktRGSEFCDGjx/fZvvSSy+VTuoYNGlMX+pdPj7f5gxOpLtj7kiS7dE/ZFlmVgvR/Tr1iWtqaqK8vLyryyLCTFEUZsfoRDoZzSRE9+twH4Su6xw8eJBbbrklrAUT4dE6ke4vRds44GwKTqT73tShmKN4RTqtpbM6QbK8CtFtOtwHYTabycvLIycnJ2yFEuFltZi4adpQnllVQlWdhx0H6lj50V6uuWxQVA8ndXtVCRBCdKNO9UGI2JeUaOX7Vw/nzyu30tjsZ8uOKlIdNq4alx/pop2S16eh60bMr70tRKwIKUBcdNFF7X6zNIzABKbPPvusywsmwi89JYHvXz2Mpau24fVrrP/yEClJNi4cEZ13hwbg8akkJUZ/p7oQ8SCkALFw4ULq6uq4/vrrMQyD119/nZycHGbMmBHu8okw65Pp4HtTh/DCO6VousGqT/aRkmRlxICMSBetXc1eCRBCdJeQRjFt3ryZ//zP/2TYsGEMHz6chx56iA0bNtCvXz/69esX7jKKMBvcrxfXXT4YAMOAv6/bxf6KxgiXqn2qFuisFkKEX0gBorKykpqamuB2TU0NVVVVYSuU6H6jz+nNzAkFQOAi/OLaUpw1zREuVfuaPTLkVYjuEFIT080338zcuXO54oorANiwYQN33HFHWAsmut8lo/rQ4PLxUfERPD6NF94JTKRLS06IdNHaaE0DbjHLxDkhwimkALF48WLGjh3L5s2bMQyDxYsXM3To0HCXTUTAtAv70+T289Wuo9S7fLzwTik/mD2SpMToyoXU7FFJddgiXQwh4lrIn/q8vDw0TWPkyJHhLI+IMJOicO2kQTS5/ew6WE9lrZuX3i3l1pnDsVmiZw6C26eSnGTFFMXzNoSIdSHdo2/YsIGZM2dyzz33APDtt99y5513hrVgInLMJhOLpgwhL8sBwAFnE8ve34WmR8/qboYh6TeECLeQAsTTTz/NihUrSE1NBWDUqFEcOHAgrAUTkZVgNXPz1cPaLFv65sa9UbVsqXRWCxFeIffyZWVltdm22aT9N945Eq3cMmN4sK3/y51HeefzA1ETJDTdwOOTICFEuIQUIBwOB0ePHg3Opv78889JSUk54+s2btzItGnTmDJlCkuXLj3lccXFxQwfPpy1a9eGWGzRXVpnW9sTAv0PHxcfYeM3hyNcqmPkLkKI8AkpQPzsZz/j9ttv5+DBg9x4443cd999bRL4tUfTNB599FGee+45ioqKWLNmDbt37273uN/85jdceumlnauBCLvcjCRumjYMa8uw0ne/KGdzaWWESxXgU3X8avT0jQgRT0IaxTR69GheeuklvvzySwDOP//8YH/EqRQXF1NQUEB+fiD528yZM1m3bh3nnHNOm+Nefvllpk2bxrffftuZ8otuUpCbwqIphbz87s7g2tZ2m5lzB2VGumg0e/z0irK5GkLEgzMGCE3T+O53v8vrr7/OpEmTQj6x0+kkNzc3uJ2Tk0NxcfFJx7z//vu8+OKLHQoQGRmOkI+NRdFavwkZDiw2C8+vKsEw4LUPdtM708GIgR0LEl1dPwXIyHREzXoWWVlnbn6NVfFcN4j/+nXUGQOE2WwmPT0dr9dLQkLo39La68g8MSPsY489xn333YfZ3LHx9TU1rg4dH0syMhxRXb/BuSnMumQAqz8pQ9UM/vR6MbfNHE7/nNA+WOGqn9fti4qlU7OyUqiqis48VmcrnusGPaN+HRVSE9OAAQNYvHgx06ZNIykpKbh/8eLFp3xNbm4uFRUVwW2n00l2dnabY7Zu3cq9994LQG1tLRs2bMBisXDVVVd1qBKie00YmYvbq/L+loP4VZ0X3inl9tkj6JMZuTufZo8fR6Ilqhc8EiLWhBQgXC4XhYWF7N27N+QTjxo1irKyMsrLy8nJyaGoqIj/+Z//aXPM+vXrg48feOABLr/8cgkOMeKK8/vh8Wp8/G0gb9Nf3y7lB7NH0DvNHpHy6AZ4fBr2hOhKCSJELDvtp+nXv/41DzzwAE888QSffPIJl1xySegntlh4+OGHWbJkCZqmMX/+fAoLC1m2bBkQWGNCxC5FUbj6ov64fSr/2lFFk9vPX4q284M5I0lPiUyHscvjlwAhRBdSjNPMerrmmmt48803T3ocSZU1zTjjuJ0w2vsgTqTrBv9Yv4tv9wbSwWemJnL7nBGkJrU/kTLc9ctIScAWwXWr47kdO57rBj2jfh112nkQx8eOaJk9K6KLyaSw4IpzGNo/DYDqBg/PF23H5fFHpDwumTgnRJc5bYDw+Xzs2bOH3bt3t3nc+iMEgMVsYtFVQxjUNzA3prLWzV+LtkckmV7rWhFCiLN32iamyZMnn/qFisK6devCUqjTkSam6OX1a/z17e0ccDYBkJ+dzK0zhpNgO9bk0x31S0q0nLKJK9ziuZkinusGPaN+HXXaHr3jRxkJcSYJVjPfv3oYf1mznUNHXZRXNvHi2lK+f/Wwbu0XcHtVku2yVoQQZ0vWbBRdKtFm4ZYZw8jNCMyXKato5KV3d+BTtW4rg2GAR9aKEOKsyZhA0eWSEq3cOnM4z67eRlWdm72HG3jl3Z3cOO3slqnddbCOLaWV1DZ6SU9JYNywbArz0to91uVRsSfIxLlYtHVfNR8XH6Gqzk1Wmp1Lz+vDuR1M5yK6htxBiLBItlu5bdbw4IJDuw/V88p7O/B38k5i18E63v2inOoGL7oB1Q1e3v2inF0H69o9XtMNGdEUg7buq+b1DXtx1rrRDXDWunl9w1627quOdNF6JAkQImxSk2wsmTWCzJYgsetgPX9+49tOpefecor04qfaD4GJc9G0TKo4s4+Lj3RovwgvCRAirFIdgSCRkRqYXV2yt7rlTqJjF+7aRm+H9kOgL6KxOTLzMUTnVNW5T7Hf080lESABQnSDXicEiV0H63m5gx3Xp0rfcaa0Hh6fhtfffR3k4uxknSKXV1ZaYjeXRIAECNFN0pITuH3WCLLTAxeA3YfqeWntDnwhXrzHDcvu0P7jNbp8kgkgRlx6Xp8O7RfhJQFCdJteyQncu2hssON67+EGXninFK/vzEGiMC+NaePzyUxNwKRAZmoC08bnn3IU0/FU3YjIrG7RcecOzGT+pEHkpNsxKQo56XbmTxoko5gi5LQzqaORzKSObRkZDsoO1vKXNduC7cr52cl8/+phYc3EajIpZPVKDPuw13iejRvPdYOeUb+OkjsI0e1Sk2zcPntkcDJdeWUTz63ZRpM7fB3Kum7QLHcRQnSIBAgREcl2K0tmjaBfVmAVuiPVzTy7ehv1Ll/Y3tPl9qPH1g2zEBElAUJETFKipWU962QgMMRx6aoSqhvCM6RRN6BZJs8JETIJECKiArmbhnNOv15AYF7D0pUlVNQ0h+X9XB4/ui53EUKEQgKEiLgEq5mbpg9lxIB0ABrdfp5dXcIBZ9d3GBoGEVvMSIhYIwFCRAWL2cTCq4bwnSG9AXB7Nf6yZjs7DtR2+Xs1e1RJwSFECCRAiKhhNilcO2kwl4zKBcCv6bz87g6+2lnVpe9jICk4hAiFBAgRVUyKwoyLCpg2Ph8IdCwv/3APG78+3KWzoT0+rdOZZYXoKSRAiKijKAqTxvRj/qRBmFrmta394gCrPynr0g7mBpfcRQhxOhIgRNQaOzSb700ditUS+G+6aZuTV/+5s8tWp/NruqTgEOI0JECIqDasIJ0ls0bgSAyk4di+v5a/rNlOY3PXTKhrlMlzQpySBAgR9fKzk7lz3rnBhYfKK5v488oSnF0wV0LXDeqbwjd7W4hYJgFCxITM1ETunDuSgpxAwrHaRi9/XllyyiVHO8Lr18KaB0qIWCUBQsQMR2Jgnesx5wTmSnj9Gi++U8pnJRVnPcKpye3H45P+CCGOJwFCxBSL2cSCKwZz5dg8IDAMdvUnZaz8eB+qdnaT3+pdvrM+hxC4wXT7AAAeVklEQVTxRAKEiDmKonDl2DxuuPIcLObAONgvtlfy/Nvbz6qpyDCgptErQUKIFhIgRMw6b3Bv7pgzklSHDYCyI43835vfcqiqqdPn1HWD2kavpOIQAgkQIsb1y0rmR9ecS352IGV4XZOPZ1aV8OVZpOfQdIPaBgkSQkiAEDEvJcnG7bNHcMGwbABUzWDFh3vOql9CbQkSkhpc9GQSIERcsJhNXHPZIK6ZOBBzS36Oz7c5WbqqhNpGb6fOqeoGNY0emUgneiwJECKuXDA8hx/MGUGvln6Jg1Uu/vBGMaWdTBuuagZ1jV4JEqJHCmuA2LhxI9OmTWPKlCksXbr0pOdXrVrF7NmzmT17NjfccAOlpaXhLI7oIfKzU7h7/igK8wKr1Lm9Gi+t3cE7m/Z3qsnJp+rUNXq7NJusELEgbAFC0zQeffRRnnvuOYqKilizZg27d+9uc0xeXh6vvPIKq1ev5q677uIXv/hFuIojehhHopWbpw/jyrF5tCSE5aPiIyxdVUJNJ9a89qk6dU0+CRKiRwlbgCguLqagoID8/HxsNhszZ85k3bp1bY75zne+Q69egW95Y8aMoaKiIlzFET2QyRSYL3HLzOGk2K1AoMnp969/y9e7j3b4fF6/Rr1L8jaJnsMSrhM7nU5yc3OD2zk5ORQXF5/y+BUrVnDZZZeFdO6MDMdZly+aSf261vgMB8MG9ebFom2U7K3G69d4bf1u9lU0snDqUJISrR06nyXRQnpK4imfz8pKOdsiR614rhvEf/06KmwBor1bcUVR2jkSNm3axIoVK/jb3/4W0rlralxnVbZolpHhkPqFycIrz+HT7GTe/eIAmm6weZuTnftrWXDFYAb17RXyeWqA6gRXsCP8eFlZKVRVNXZhqaNHPNcNekb9OipsTUy5ubltmoycTifZ2dknHVdaWspDDz3E//3f/5Genh6u4giBSVG49Lw+/PCac8lOtwOB/Et/WbOdNZ+WdWghIrdXpVZGN4k4F7YAMWrUKMrKyigvL8fn81FUVMTkyZPbHHP48GHuuece/vu//5uBAweGqyhCtNEn08GPrhnFxecGmkAN4NOtFfx+xbfsrwj9G6TXr1HT4JHcTSJuha2JyWKx8PDDD7NkyRI0TWP+/PkUFhaybNkyABYuXMgf//hH6urqeOSRRwAwm8288cYb4SqSEEFWi4lZFw9gWEE6b2zYQ12Tj+oGD0tXlTDh3FymXJBPgtV8xvOomkFNg4eUJBv2hLB9nISICMWIsXF7lTXNOOO4nVD6ILpfyb5q3tl0gJrjZlwn2y307pWIqhmkpyQwblg2hXlppz2PzWJiUEEGdbUdX+lu675qPi4+QlWdm6w0O5ee14dzB2Z2+DzhsOazMj786hAuj4oj0cLl5/dj1oQBkS5Wl5M+iJPJVx7Ro+06WMcHXx0mMcFChgL1TT403aDJrdLkbsKeYMavGbz7RTnAaYOET9WprHXjdvlIsVsxmdoflHGirfuqeX3D3uC2s9Yd3I50kFjzWRlrPikDAoNMmpr9we14DBKiLUm1IXq0LaWVwceJNgtZ6fZgLicIzMKurG3G5fazebszpHO6vSpH6900e0Jbm+Lj4iMd2t+dPvzqUIf2i/giAUL0aCcm8jMpCiYTmM2BBIAQWEio3uVjR3k95ZWhNUHoBjQ0+zla78brO/3oqKo69yn2d3zGd1c71QJMLlnDu0eQACF6tPSUhJP2mU0mrGYzWWmJpDpstE7f8as6f3qrhBUf7qGxObQZ1apmUNvkpbreg9fffqDISrOfYv+pJ+N1l2R7+5MIHafYL+KLBAjRo40bdvLcnKREC45EC4qikGy3kp1mJ9F2bETTlzureOof37Dh60P41dCGuPo1ndpGLzUNJweKS8/r0+5rTrW/O11+fr8O7RfxRTqpRY/W2um8pbSS2kZvy4il/Db7MtPtzJhQgMmksOaTMpy1brx+jXe/KOfzbU6mju/PeYMzMZ0iU8DxfKqOr9GLzWLCkWglwWYOdkQHRjF5yEpLjJpRTK0d0R9+dYhmj4rDbo3bUUziZDLMNcpE4zDQrhTr9dN0gy+2OXn/Xwdxe9Xg/n5ZDqZd0J/x5/XtUP3MJgV7goWkBEvIo54ipScMA433+nWUBIgoE+sX0DOJl/q5vSoffnWIT7dWoB23LOnQgnQmn98vuEZ2qBQCk/dsVjMJVjNWS/S1/vaEC2i816+jpIlJiE6wJ1i4+qICLhyRwz+3lPPN7moAduyvZcf+Wob1T+fKsf3olxVaoDBoaX5SdZrcfkwmhUSbGbvNEpXBQvQMEiCEOAsZqYlcP7mQy0b35b3N5ew4UAdA6YFaSg/UMrwgnSu+04+8EANFK103aPaoNHtULCaFBFvgzsIWQvoPIbqKBAghukCfTAc3Tx9GTbOfN9bvYu/hBgC2769l+/5aCvN6MWlMPwb2STll2vtTUXUD1aPi8qiYFEiwmkm0WbBZTR0+lxAdIQFCiC50Tl4aS2aNYN+RBtb962AwUOw6WM+ug/X0z0lm4nl9GV6Q3qlOad0At0/D7dNQWoKFzWLCajFhMUvAEF1LAoQQYTCwT2owUGz4+jA7ywNNTwecTbz6z51kpiZyyahcvjMkq9PNRoYBHp+Gp2WmtgLYrGYSbWYSbOaQht0KcToSIIQIo4F9UhnYJ5XDR11s+PoQW/fVYBhQ3eBh1SdlvLe5nHHDsrloRA4ZqWc3c9ogsEaF16+huAKjolr7LlrThgjRERIghOgGfXs7WHjVEGoaPHyytYJ/lVbiU3U8Po2Pi4/wSfERhvZPZ/yIbIbkpZ31nIjjR0U1EhgVZTWbsJgVLGZTy48iTVLitCRACNGNMlITmX3xAK4am8eW0ko+K6mgrsmHwbGRT2nJNsYNy2bskCx6JZ+cK6ozdN3Aq2t4T8ixZzEpmM0mTAqYTAomk4LFZMLcEkhEzyYBQogIsCdYmDi6L5eM6kPpgVo2lTjZfagegLomH+9vOci6LQc5J68XY4dmM7wgPSzzIVTdQNXbTyKoEMhoa7GYsJoVfH4N3TCkb6MHkQAhRASZTAojBmQwYkAG1fUevtju5F87q2j2qBgcG/2UaDNz7qBMzi/sTUFuSrdcpA0CSQb9mo4bsNS5qal1YzYF7i7MLXccJkXBbFIwmwO/pdkqfkiAECJKZPZK5OqLCphyQT7b99fyrx1V7DpYFxyttKW0ki2llfRy2Bg1OJPzBmfSr7ej2y/Imm6gneKuAwj2d5jNgeYqkwkUWsqoBO5MFEVBUZCAEuUkQAgRZSxmE6MGZTJqUCb1Lh9f76riq11HqawNLCxU7/LxcfERPi4+QkZKAiMHZnDuoAz6ZSVHRfNPa38HIawppABmc2tAMaEoLYs2tSzcpChyVxJJEiCEiGK9HDYmjenHZaP7cqS6ma93HeXbvdXUuwILFtU0evmo+AgfFR8h1WFjeEE6wwvSGdQ3NSY6mQ0Ciyqpmgac+q5EARSTgoljdx+nChqBZq/A84HBYEpw0SfdMGhNTxo8pqWZTNX0mOtj0Q0DXQ/8aLpxbNsIBGojuG1Isj4h4pWiKPTt7aBvbwfTL+pPubOJ4j3VbN1XTWNz4Kt6g8vH59ucfL7Nic1i4py8XgzNT2NIflqXjYaKFAMwdAM9uBWG97CYqal1tzSBHReIjgswcCwwmZS2wcowjEDJDAIX5pbfAMfnzDaOK//x524NaEpLwIKWgKYfO5fe8m9gHBfowkkChBAxxqQoFOSmUJCbwsyLCzhY2UTJvhpK9tVQ07LGtk/V2VZWy7ayWgBy0u0U5qVxTl4vBuSmSNK/0zBouaAHr8DhvBJH92oLEiCEiGEmRaF/Tgr9c1KYfmF/KuvclLYkCCyvbApe45y1bpy1bj7+9ghmk0J+TjKD+qQyqG8q+dkpklJctEsChBBxQlEUctKTyElPYtKYfjR7VHYfqmNneR27yutpdAeaojTdoOxII2VHGln/5SEsZoW8rGQKclMYkBsINvYEuTQICRBCxK2kRAvnDe7NeYN7YxgGzlo3ew7Vs/tgPWUVjXj9gU5hVTMoq2ikrKKRDS2vzUqz0z8nmf7ZyeRlJ5OdnhS5ioiIkQAhukVrH19gGGPggdK6Xzm5s6712NbXGbSOymjprCP0jrpj73GMHt1Nv11OURRyM5LIzUjiklF90HSDw0dd7D1cT9mRRvY7G4NZYQGq6txU1bn5144qAKxmE/37pJCdZqdfS2d5Vpodc5Svoy3OjgSIONfm4th6UW7ngqyccAwt49Hh2IW67YVbOW7SU+AFx8534vPHLiJZvZOxdvHwi8DQxWPB4vj3P9WQRd0w0DQdVQsMAWypWcuIlNYRJUpwXL6itMwsVnVUTUdV9WMjS1reW4+h5d3NJoX87GTys5OZNCZQdmdNM/srGjngbOJAZSM1Dd7g8X5NZ8/BevYcrA/us5gVstOT6JORRG5mEjkZSeSk20m2W2XeQpyQABEm7V6Yg/vajuM+/uLsSLTiSbS0mW16/MXr2PGtF8D2nzvdOPF4Y2q9qnfwNSaLGWsHPgEKgQV6Ek4zAigrKxmL3hJENANN1/GrgXQV0Rw/TIpCn0wHfTIdXDQysK/J7edQVRPllU0cqnJxuNoVHFILgaapw0ddHD7qanOupEQL2el2stPsZAV/EumVnBBTcwxEDAYIe6KFpNYONKXNr8DjE/4Dthm/3N6jE/6/Bi+0nPxCpc0xoX1j7qi0lAT8Hl+nXy8iz2RSsJnM2Kxt9/tVHU3XW5rJAt/afX4Nn6pHpqBnkGy3MrR/OkP7pwOQnp5EWXkth4+6OHTURUVNM0eqm6lt9LZ5XbNHDXaCH89iVujdy05maiKZvRLITE0kIzWRjNQEUh0J0lwVhWIuQKQk2fA4bJEuhhAdZrWYsHLCcFK7FV038Pg0fKrWkjE1MuU7E0VR6JWcQK/kBIYPyAju9/hUnDVuKmqacdY246wJ9F80udvm2lA1g4qaZipqmk86t0lRSEuxkZacQHpKQvB3r+TAvl4OW0zMDI83MRcghIg3JpNCUqKFpJaPo1/V8akaXp+GX9WjfCoVJNoswYl7x2v2qMHO7qP1bqrqPByt91DT4EE7IQrqhkFNg7dNv8eJHIkWejlspDpspCQFfqcmWUlJspHc+ttuwWySQNJVJEAIEWWsFhNWiwlHohXdMPD7AwHDr+oxETBaJSW2Hzh03aCh2cfReg+1jV5qGgJBo67JR22j96Q7j1Yuj4rLo3K4+uQ7kOPZEywk260k26047BYciYHHSYkWHIkWkhKtJCVYAkE50YLVbOox/XUdJQFCiChmUpTAutK2QMe4YRgtye0CHd9qy6iqaG2Wao/JpJCWHGhGao9P1ahr8lHf5KW+yUddk5d6l48Gl496l4/6Jl9wDkd73F4Vtzdw9xIKi1kJBJUkGzaLCbvNgj3BTKLNQmKCmURby2ObOfiTYLWQYA2s+W2zmM96idhoJQFCiBiiKApWi4LVYsJ+3H5dN4Ijp1Q9EDi0lgyfscZmMZOdFhgFdSo+v0Zjs596l48mt4/GZj+NzX6a3G1/XG7/Gf8NVM0Ivr6zrBYTCVYzNmvrbzM2i+mk31aLCZvFHLxLPP7HYm553PLb0vLY0rK2RiRGgIU1QGzcuJHHHnsMXddZsGABP/jBD9o8bxgGjz32GBs2bCAxMZFf//rXjBw5MpxFEiIumUwKO/fX8XHxEarq3GSl2bn0vD6MHJDBc2u2saW0Er9mYDUrjCnszZjCLL7Y7qS63kN6SgLjhmVTmJfGroN1bCmtpLbRG9x/sKqJL7Y5afZpJNnMjB+RwxXn552yLO2dAzhpX2FeWsivL8xL44OvDgbK4VVJSrCcthyt56hp8AQ61QvSyUqzs/tQPaX7a2ls9mOzmshIScRqMeH2qfhVg8ZmH26vitendagpr7X5j9BuWjqldSU/S8t64a3rhreuK966QFPryn5mU+BYc8vzo4bmdPg9FcMIz+hsTdOYNm0af/3rX8nJyeG6667jqaee4pxzzgkes2HDBl5++WWeffZZvvnmGx577DGWL19+xnNXVTWe8ZhYlZWVIvWLYZGq39Z91by+Ye9J+x2JFra3ZHRtZRDoH8hq+Ybeegk4f0gWX+6sCiYYNYB6l5dmt79lsqASPHby2DyuOL8fbZKeErgwv/tFeZv38/jUwBwSW9vvo9PG558UJNp7PUBuhp1v91SftP/y7/Q7KUic6hwjB6ZTsq/2pP2t5cjIcFBTE5jTsbO8lrWflwcnYepGIN34uYMySEtOwOPT8Ppbfloe+1U9uO1T9cBPy/5osPp/5nb4NWG7gyguLqagoID8/HwAZs6cybp169oEiHXr1jFv3jwURWHMmDE0NDRQWVlJdnZ2uIolRFz6uPhIu/tL9598QQRwe9Tg49YO2o++OUxKkq3NxKJmt4phBNYoOL4jd3NpJddPLgxut85kf3PjHizmY8cZBPoEAJIS204M+Wb3Ub4zJKvl9YF9X+86ykmDkAzaDQ4Am7dXcvWFBceVA77aWdXuvMnN2ytx2K0n7f9qZxUjB2QEmoZastp+veto4Js4bU9U7/Ix55KB7ZblVF+1dcNAbQkWvpY7DZ967PFJP5rWMmO/pa+pZaKldty2qumo+rFsAFrLttpynKa3TtQ8u+//YQsQTqeT3Nzc4HZOTg7FxcWnPSY3Nxen03nGANGZlZFiidQvtkWifrVNvnZTdutGO5PMjcCF+8TjXR6VjNTEE14fuMAcm/Uf+N3sUdutZ53Lf9JaE7oOKAQ72ls1eVQG5GectC/RdvJlSdUNbJaTZ7B7fBqFA3u3rYdXazcb7dF6D73b6ddo9mkMHRwIVK3Pu7zb2j2H+7hje4KwBYj2Wq5OHEoWyjFCiDN76ieTIl0E4OzL0RX1iJZzxIOwzSjJzc2loqIiuN3encGJx1RUVEjzkhBCRImwBYhRo0ZRVlZGeXk5Pp+PoqIiJk+e3OaYyZMn89Zbb2EYBl9//TUpKSkSIIQQIkqErYnJYrHw8MMPs2TJEjRNY/78+RQWFrJs2TIAFi5cyKRJk9iwYQNTpkzBbrfz+OOPh6s4QgghOihsw1yFEELENslqJYQQol0SIIQQQrQrqnMxeb1eFi9ejM/nC87M/vGPf0xdXR0//elPOXToEP369eN///d/6dWrV6SL2ymt/TM5OTk888wzcVW3yZMn43A4MJlMmM1m3njjjbiqX0NDAw899BA7d+5EURQef/xxBg4cGBf127t3Lz/96U+D2+Xl5fz4xz9m3rx5cVG/F154geXLl6MoCkOGDOGJJ57A7XbHRd0AXnzxRZYvX45hGCxYsIDvf//7nfrsRfUdhM1m48UXX2TVqlW89dZbfPTRR3z99dcsXbqUCRMm8N577zFhwgSWLl0a6aJ22ksvvcTgwYOD2/FUNwj8R125ciVvvPEGEF/1e+yxx5g4cSJr165l5cqVDB48OG7qN2jQIFauXBn829ntdqZMmRIX9XM6nbz00ku8/vrrrFmzBk3TKCoqiou6AezcuZPly5ezfPlyVq5cyYcffkhZWVmn6hfVAUJRFBwOBwCqqqKqKoqiBFN0AMybN4/3338/ksXstIqKCj788EOuu+664L54qdupxEv9mpqa2Lx5c/BvZ7PZSE1NjZv6He+zzz4jPz+ffv36xU39NE3D4/Ggqioej4fs7Oy4qduePXsYPXo0drsdi8XCBRdcwD//+c9O1S+qAwQE/pBz587l4osv5uKLL2b06NFUV1cH50tkZ2dTU1MT4VJ2zuOPP86//du/YTou+Uy81K3VbbfdxrXXXss//vEPIH7qV15eTkZGBg8++CDz5s3j5z//Oc3NzXFTv+MVFRUxa9YsID7+fjk5Odx6661cccUVXHrppSQnJ3PppZfGRd0AhgwZwpYtW6itrcXtdrNx40YqKio6Vb+oDxBms5mVK1eyYcMGiouL2blzZ6SL1CU++OADMjIyOPfccyNdlLBZtmwZb775Js8++yyvvvoqmzdvjnSRuoyqqmzbto2FCxfy1ltvYbfbY7ZJ4nR8Ph/r169n+vTpkS5Kl6mvr2fdunWsW7eOjz76CLfbzcqVKyNdrC4zePBglixZwq233sqSJUsYOnQoZvPJeaxCEfUBolVqaioXXnghH330EZmZmVRWVgJQWVlJRkbGGV4dfb788kvWr1/P5MmTuffee9m0aRP33XdfXNStVU5OIP98ZmYmU6ZMobi4OG7ql5ubS25uLqNHjwZg+vTpbNu2LW7q12rjxo2MHDmS3r0DCfHioX6ffvopeXl5ZGRkYLVamTp1Kl999VVc1K3VggULePPNN3n11VdJS0ujoKCgU/WL6gBRU1NDQ0MDAB6Ph08//ZRBgwYFU3QAvPXWW1x55ZWRLGan/OxnP2Pjxo2sX7+ep556iosuuojf/OY3cVE3gObmZpqamoKPP/nkEwoLC+OmfllZWeTm5rJ3b2ANhs8++4zBgwfHTf1aFRUVMXPmzOB2PNSvb9++fPPNN7jdbgzDiMu/XXV1ID364cOHee+995g1a1an6hfVM6lLS0t54IEH0DQNwzCYPn06d999N7W1tfzkJz/hyJEj9OnTh9/97nekpbW/OlUs+Pzzz3n++ed55pln4qZu5eXl/OhHPwIC/UizZs3irrvuipv6AWzfvp2f//zn+P1+8vPzeeKJJ9B1PW7q53a7ufzyy3n//fdJSQmk9o6Xv9/TTz/N22+/jcViYfjw4Tz22GO4XK64qBvAokWLqKurw2Kx8OCDDzJhwoRO/e2iOkAIIYSInKhuYhJCCBE5EiCEEEK0SwKEEEKIdkmAEEII0S4JEEIIIdoV1dlchTidBQsW4PP58Pv9lJWVUVhYCMCIESN44oknIly60JSUlFBeXh5XM5VF/JBhriLmHTx4kPnz5/P5559HuignUVUVi+XU38OWL1/Op59+ym9/+9suP7cQZ0v+d4m4tGLFCv7+97+jaRqpqak88sgjDBgwgOXLl7N27VocDgc7d+6kT58+/Md//AdPPvkk5eXljB49mieffBJFUbjvvvuw2+0cOHCAiooKLrzwQn7xi19gtVppbGzk8ccfZ9euXXi9Xi6++GLuv/9+TCYTCxcuZPz48Xz11VckJSXx9NNPBycJer1eRo8ezSOPPEJDQwN//OMfcblczJ07lwsvvJDFixezaNEiPvnkEwD2798f3N6/fz8LFy7k+uuvZ9OmTVx77bXMnTuXp556ii1btuDz+Rg+fDi//OUvsdvtEf4LiLhgCBHjysvLjfHjxwe3N23aZNxxxx2G1+s1DMMw1q1bZyxevNgwDMN47bXXjPHjxxsVFRWGYRjGrbfeasybN89obGw0fD6fMWPGDGPTpk2GYRjGz372M2Pu3LmGy+UyfD6fcdNNNxl/+9vfDMMwjPvvv99YvXq1YRiGoWma8eMf/9hYsWKFYRiGccMNNxg//OEPDVVVg8/X1dUFH997773Ga6+9FizPT37yk2DZy8rKjIsvvrjd7bKyMmPIkCHG2rVrg88//fTTxjPPPBPcfuKJJ4zf/e53Z/cPKkQLuYMQcWf9+vVs27aNBQsWAGAYBi6XK/j82LFjg4kER4wYgcfjITk5GYChQ4dy4MABLrzwQgBmzJhBUlISEMih/+GHH7Jw4UI++OADSkpKePbZZ4FArrD+/fsH32P27NnBDJq6rrN06VI+/vhjdF2nrq6u0yuVJSUlMW3atDZ1dbvdFBUVAYHsqyNHjuzUuYU4kQQIEXcMw+C73/0ud999d7vPJyQkBB+bTKaTtlVVPeV5FUUBAhf9Z555hr59+7Z7bGtQAVi5ciXFxcX87W9/w+Fw8Ic//IEjR460+zqz2Yyu68Ftr9d7yvO2lulXv/oVF1xwQbvnE+JsyDBXEXdas1Y6nU4gkCxw69atnTrXO++8g9vtxu/3s3r16uCdxeTJk1m6dCmapgGBzMPl5eXtnqOxsZH09HQcDgf19fXBb/sADoeDxsbG4HZ2djYejyd4rjVr1pyxrs8//3wwkDQ1NbFnz55O1VWIE0mAEHHnoosu4u677+aOO+5gzpw5zJ49mw8//LBT5xo7dix33XUXs2bNIj8/P7jE6C9+8Qt0XWfu3LnMnj2b22+/naqqqnbPcc0111BXV8esWbO4995723zbv+SSS2hsbGTOnDk8/vjj2Gw2HnjgAW6++WZuvPFGrFbract35513MnjwYK677jpmz57N4sWL2bdvX6fqKsSJZJirEKdw3333MXbsWBYuXBjpoggREXIHIYQQol1yByGEEKJdcgchhBCiXRIghBBCtEsChBBCiHZJgBBCCNEuCRBCCCHa9f8DD1fSWw4ItgAAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "sns.set(color_codes=True)\n",
+    "plt.xlim(30,90)\n",
+    "plt.ylim(0,1)\n",
+    "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
+   ]
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Hide code",
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}