diff --git a/module4/exo/challenger.ipynb b/module4/exo/challenger.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..d305d445dd5c921fee32756dde663a9fdb407c32 --- /dev/null +++ b/module4/exo/challenger.ipynb @@ -0,0 +1,804 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this document we reperform some of the analysis provided in \n", + "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n", + "\n", + "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Technical information on the computer on which the analysis is run" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.7.7 (default, Mar 23 2020, 22:36:06) \n", + "[GCC 7.3.0]\n", + "uname_result(system='Linux', node='BI-platform', release='4.15.0-65-generic', version='#74-Ubuntu SMP Tue Sep 17 17:06:04 UTC 2019', machine='x86_64', processor='x86_64')\n", + "IPython 7.13.0\n", + "IPython.core.release 7.13.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", + "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.2\n", + "distutils 3.7.7\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.5.1\n", + "ipywidgets._version 7.5.1\n", + "jedi 0.16.0\n", + "json 2.0.9\n", + "jupyter_client 6.1.2\n", + "jupyter_client._version 6.1.2\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 3.1.3\n", + "matplotlib.backends.backend_agg 3.1.3\n", + "mkl 2.3.0\n", + "numpy 1.18.1\n", + "numpy.core 1.18.1\n", + "numpy.core._multiarray_umath 3.1\n", + "numpy.lib 1.18.1\n", + "numpy.linalg._umath_linalg b'0.1.5'\n", + "pandas 1.0.3\n", + "parso 0.6.2\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.4\n", + "ptyprocess 0.6.0\n", + "pygments 2.6.1\n", + "pyparsing 2.4.6\n", + "pytz 2019.3\n", + "re 2.2.1\n", + "scipy 1.4.1\n", + "scipy._lib._uarray 0.5.1+5.ga864a57.scipy\n", + "scipy._lib.decorator 4.0.5\n", + "scipy._lib.six 1.2.0\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.10.0\n", + "seaborn.external.husl 2.1.0\n", + "six 1.14.0\n", + "statsmodels 0.11.1\n", + "statsmodels.__init__ 0.11.1\n", + "statsmodels.api 0.11.1\n", + "statsmodels.tools.web 0.11.1\n", + "traitlets 4.3.3\n", + "traitlets._version 4.3.3\n", + "urllib.request 3.7\n", + "zlib 1.0\n", + "zmq 17.1.2\n", + "zmq.sugar 17.1.2\n", + "zmq.sugar.version 17.1.2\n" + ] + } + ], + "source": [ + "def print_imported_modules():\n", + " import sys\n", + " for name, val in sorted(sys.modules.items()):\n", + " if(hasattr(val, '__version__')): \n", + " print(val.__name__, val.__version__)\n", + "# else:\n", + "# print(val.__name__, \"(unknown version)\")\n", + "def print_sys_info():\n", + " import sys\n", + " import platform\n", + " print(sys.version)\n", + " print(platform.uname())\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import statsmodels.api as sm\n", + "import seaborn as sns\n", + "\n", + "print_sys_info()\n", + "print_imported_modules()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loading and inspecting data\n", + "Let's start by reading data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "scrolled": true + }, + "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/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": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYMklEQVR4nO3df5TU9X3v8edrAWERIgRSagWiVuqNRy3KBjX2B0STg54TaK6a4D3RNC2h90ROTkyaaHtzreWm59zYJLa5sYnEaxrtSYhKo9xeev0RJak9/gCVgGKwWzW4YFA3qKwg7LLv+8d8txmG2eU7y35ndubzepyzZ+f7/X7mO+/PfHf2Nd8f8xlFBGZmlq62RhdgZmaN5SAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcYUEg6VZJr0h6epDlkvR1SZ2SNks6u6hazMxscEXuEfw9sGiI5RcBc7Kf5cA3C6zFzMwGUVgQRMRPgF8O0WQJcFuUPApMkXR8UfWYmVl1Yxv42CcAL5VNd2XzXq5sKGk5pb0G2tvb582aNasuBebV399PW1vrnW5p1X5B6/bN/Wo+9erbc88991pEvKvaskYGgarMqzreRUSsAlYBdHR0xMaNG4usq2br169nwYIFjS5jxLVqv6B1++Z+NZ969U3Szwdb1siI7QLK39rPBHY2qBYzs2Q1MgjWAldmVw+dC7wREYcdFjIzs2IVdmhI0veBBcB0SV3AXwDjACLiW8A64GKgE9gLfKKoWszMbHCFBUFEXH6E5QFcVdTjm5lZPq15Gt7MzHJzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWuEKDQNIiSdskdUq6tsry2ZIekvSUpM2SLi6yHjMzO1xhQSBpDHATcBFwGnC5pNMqmn0RuCMizgKWAn9XVD1mZlZdkXsE84HOiHg+Ig4Aq4ElFW0CeEd2+zhgZ4H1mJlZFYqIYlYsXQosiohl2fQVwDkRsaKszfHAfcBU4Fjgwoh4osq6lgPLAWbMmDFv9erVhdQ8XD09PUyaNKnRZYy4Vu0XtG7f3K/mU6++LVy48ImI6Ki2bGyBj6sq8ypT53Lg7yPiq5LOA26XdHpE9B9yp4hVwCqAjo6OWLBgQRH1Dtv69esZbTWNhFbtF7Ru39yv5jMa+lbkoaEuYFbZ9EwOP/Tzx8AdABHxCDABmF5gTWZmVqHIINgAzJF0kqRjKJ0MXlvRZjtwAYCk91AKglcLrMnMzCoUFgQR0QesAO4FnqV0ddAzklZKWpw1+xzwSUk/Bb4P/GEUddLCzMyqKvIcARGxDlhXMe+6sttbgfOLrMHMzIbmTxabmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmlrhCg0DSIknbJHVKunaQNh+RtFXSM5K+V2Q9ZmZ2uLF5Gkk6PSKermXFksYANwEfALqADZLWRsTWsjZzgD8Dzo+I3ZJ+rZbHMDOzo5d3j+Bbkh6X9ClJU3LeZz7QGRHPR8QBYDWwpKLNJ4GbImI3QES8knPdZmY2QhQR+RqW3r3/EXAZ8DjwnYi4f4j2lwKLImJZNn0FcE5ErChrczfwHHA+MAa4PiL+X5V1LQeWA8yYMWPe6tWr8/WuTnp6epg0aVKjyxhxrdovaN2+uV/Np159W7hw4RMR0VF1YUTk/qH0z/oSYAfwLPAz4D8P0vYy4Jay6SuA/1XR5p+AHwLjgJMoHUKaMlQN8+bNi9HmoYceanQJhWjVfkW0bt/cr+ZTr74BG2OQ/6u5Dg1JOlPSjdk///cDH4qI92S3bxzkbl3ArLLpmcDOKm3uiYjeiHgB2AbMyVOTmZmNjLznCL4BPAn8dkRcFRFPAkTETuCLg9xnAzBH0kmSjgGWAmsr2twNLASQNB34LeD52rpgZmZHI9dVQ8DFwL6IOAggqQ2YEBF7I+L2aneIiD5JK4B7KR1SujUinpG0ktIuytps2QclbQUOAp+PiO6j7JOZmdUgbxA8AFwI9GTTE4H7gPcNdaeIWAesq5h3XdntAD6b/ZiZWQPkPTQ0ISIGQoDs9sRiSjIzs3rKGwRvSTp7YELSPGBfMSWZmVk95T009BngTkkDV/0cD3y0mJLMzKyecgVBRGyQ9J+AUwEBP4uI3kIrMzOzusi7RwDwXuDE7D5nSSIibiukKjMzq5u8g87dDvwmsInSZZ4AATgIzMyaXN49gg7gtOxyTzMzayF5rxp6Gvj1IgsxM7PGyLtHMB3YKulxYP/AzIhYXEhVZmZWN3mD4PoiizAzs8bJe/nojyW9G5gTEQ9Imkhp/CAzM2tyeYeh/iRwF3BzNusESiOHmplZk8t7svgqSt8i9iZARPwb4O8XNjNrAXmDYH+UvncYAEljKX2OwMzMmlzeIPixpD8H2iV9ALgT+D/FlWVmZvWSNwiuBV4FtgB/Quk7Bgb7ZjIzM2siea8a6ge+nf2YmVkLyTvW0AtUOScQESePeEVmZlZXtYw1NGACcBnwzpEvx8zM6i3XOYKI6C772RERfwO8v+DazMysDvIeGjq7bLKN0h7C5EIqMjOzusp7aOirZbf7gBeBj4x4NWZmVnd5rxpaWHQhZmbWGHkPDX12qOUR8bWRKcfMzOqtlquG3guszaY/BPwEeKmIoszMrH5q+WKasyNiD4Ck64E7I2JZUYWZmVl95B1iYjZwoGz6AHDiiFdjZmZ1l3eP4HbgcUk/pPQJ4w8DtxVWlZmZ1U3eq4b+StI/A7+bzfpERDxVXFlmZlYveQ8NAUwE3oyIvwW6JJ1UUE1mZlZHeb+q8i+Aa4A/y2aNA/6hqKLMzKx+8u4RfBhYDLwFEBE78RATZmYtIW8QHIiIIBuKWtKxxZVkZmb1lDcI7pB0MzBF0ieBB/CX1JiZtYS8Vw19Jfuu4jeBU4HrIuL+QiszM7O6OOIegaQxkh6IiPsj4vMR8ad5Q0DSIknbJHVKunaIdpdKCkkdg7UxM7NiHDEIIuIgsFfScbWsWNIY4CbgIuA04HJJp1VpNxn4NPBYLes3M7ORkfeTxW8DWyTdT3blEEBEfHqI+8wHOiPieQBJq4ElwNaKdv8DuAH407xFm5nZyMkbBP83+6nFCRw6OmkXcE55A0lnAbMi4p8kDRoEkpYDywFmzJjB+vXrayylWD09PaOuppHQqv2C1u2b+9V8RkPfhgwCSbMjYntEfHcY61aVeVG27jbgRuAPj7SiiFgFrALo6OiIBQsWDKOc4qxfv57RVtNIaNV+Qev2zf1qPqOhb0c6R3D3wA1Ja2pcdxcwq2x6JrCzbHoycDqwXtKLwLnAWp8wNjOrryMFQfm7+pNrXPcGYI6kkyQdAyzlV19sQ0S8ERHTI+LEiDgReBRYHBEba3wcMzM7CkcKghjk9hFFRB+wArgXeBa4IyKekbRS0uLayjQzs6Ic6WTxb0t6k9KeQXt2m2w6IuIdQ905ItYB6yrmXTdI2wW5KjYzsxE1ZBBExJh6FWJmZo1Ry/cRmJlZC3IQmJklzkFgZpY4B4GZWeKSCYLunv389KXX6e7Z3+hSzKwG3T372dd70K/dAiURBPds2sH5X36Qj93yGOd/+UHWbtrR6JLMLIeB1+4Lr77l126BWj4Iunv2c82azbzd28+e/X283dvPF9Zs9rsLs1Gu/LV7MMKv3QK1fBB07d7HuLZDuzmurY2u3fsaVJGZ5eHXbv20fBDMnNpOb3//IfN6+/uZObW9QRWZWR5+7dZPywfBtEnjueGSM5kwro3J48cyYVwbN1xyJtMmjW90aWY2hPLX7hjJr90C5f1imqa2eO4JnH/KdLp272Pm1Hb/IZk1iYHX7uOPPMy/Lv4dv3YLkkQQQOndhf+IzJrPtEnjaR83xq/fArX8oSEzMxuag8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcoUEgaZGkbZI6JV1bZflnJW2VtFnSjyS9u8h6zMzscIUFgaQxwE3ARcBpwOWSTqto9hTQERFnAncBNxRVj5mZVVfkHsF8oDMino+IA8BqYEl5g4h4KCL2ZpOPAjMLrMfMzKpQRBSzYulSYFFELMumrwDOiYgVg7T/BvCLiPhSlWXLgeUAM2bMmLd69epCah6unp4eJk2a1OgyRlyr9gtat2/uV/OpV98WLlz4RER0VFs2tsDHVZV5VVNH0seADuD3qy2PiFXAKoCOjo5YsGDBCJU4MtavX89oq2kktGq/oHX75n41n9HQtyKDoAuYVTY9E9hZ2UjShcB/A34/IvYXWI+ZmVVR5DmCDcAcSSdJOgZYCqwtbyDpLOBmYHFEvFJgLWZmNojCgiAi+oAVwL3As8AdEfGMpJWSFmfN/hqYBNwpaZOktYOszszMClLkoSEiYh2wrmLedWW3Lyzy8ZtZd89+unbvY+bUdqZNGj9ibZtJq/arKJ279rB7by+du/ZwyozJjS7HmkihQWDDc8+mHVyzZjPj2tro7e/nhkvOZPHcE466bTNp1X4V5bq7t3Dbo9v53Bl9XH3jT7jyvNmsXHJGo8uyJuEhJkaZ7p79XLNmM2/39rNnfx9v9/bzhTWb6e45/Dx6LW2bSav2qyidu/Zw26PbD5l32yPb6dy1p0EVWbNxEIwyXbv3Ma7t0M0yrq2Nrt37jqptM2nVfhVl00uv1zTfrJKDYJSZObWd3v7+Q+b19vczc2r7UbVtJq3ar6LMnTWlpvlmlRwEo8y0SeO54ZIzmTCujcnjxzJhXBs3XHJm1ZOltbRtJq3ar6KcMmMyV543+5B5V5432yeMLTefLB6FFs89gfNPmZ7ripla2jaTVu1XUVYuOYMrzz2RLU88ygNXn+sQsJo4CEapaZPG5/7nV0vbZtKq/SrKKTMm0zVxnEPAauZDQ2ZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIKDQJJiyRtk9Qp6doqy8dL+kG2/DFJJxZZj1mtunv289OXXqe7Z/+Q7Ta+0M3X7tvGxhe6R2ydtbbt3LWH3Xt76dy154hta1FUvbU8/r7eg7mfg7s2vtRyz0GR6wUYO+JrzEgaA9wEfADoAjZIWhsRW8ua/TGwOyJOkbQU+DLw0aJqMqvFPZt2cM2azYxra6O3v58bLjmTxXNPOKzdx255lIc7SwHw9Qc7+d1TpnH7snOPap21tr3u7i3c9uh2PndGH1ff+BOuPG82K5ecMcyeF19vrY//6ff0cvWXH8z1HAxoleegyPUOKHKPYD7QGRHPR8QBYDWwpKLNEuC72e27gAskqcCazHLp7tnPNWs283ZvP3v29/F2bz9fWLP5sHdjG1/o/o8QGPAvnd1V9wzyrrPWtp279hzyDxDgtke2H/W74qLqHc7jH4xI8jkocr3lFBEjtrJDVixdCiyKiGXZ9BXAORGxoqzN01mbrmz637M2r1WsazmwPJs8FdhWSNHDNx147Yitmk+r9guO0DeNa584durxv6W2tjED86K//2Df7pefi959ewfmjZk8/TfGHDvl+Mr7H3zr9ZcP7nlt53DWWWvbtonHTRv7jnedCHBw7xuMmXgcAH1vvvpi/943jnys6iifg1rbDufxB/qV5zko1yTPwYj8Lebw7oh4V7UFhR0aAqq9s69MnTxtiIhVwKqRKKoIkjZGREej6xhprdovaN2+SdrY98Yr7lcTGQ1/i0UeGuoCZpVNzwR2DtZG0ljgOOCXBdZkZmYVigyCDcAcSSdJOgZYCqytaLMW+Hh2+1LgwSjqWJWZmVVV2KGhiOiTtAK4FxgD3BoRz0haCWyMiLXA/wZul9RJaU9gaVH1FGzUHrY6Sq3aL2jdvrlfzafhfSvsZLGZmTUHf7LYzCxxDgIzs8Q5CIZB0ouStkjaJGljNu96STuyeZskXdzoOmslaYqkuyT9TNKzks6T9E5J90v6t+z31EbXWatB+tUK2+vUsvo3SXpT0meafZsN0a9W2GZXS3pG0tOSvi9pQnZBzWPZ9vpBdnFNfevyOYLaSXoR6Cj/4Juk64GeiPhKo+o6WpK+C/xLRNyS/TFOBP4c+GVE/M9svKipEXFNQwut0SD9+gxNvr3KZUO67ADOAa6iybfZgIp+fYIm3maSTgAeBk6LiH2S7gDWARcD/xgRqyV9C/hpRHyznrV5j8AAkPQO4PcoXclFRByIiNc5dBiQ7wJ/0JgKh2eIfrWaC4B/j4if0+TbrEJ5v1rBWKA9+9zUROBl4P2UhtiBBm0vB8HwBHCfpCey4S8GrJC0WdKtzbY7DpwMvAp8R9JTkm6RdCwwIyJeBsh+/1ojixyGwfoFzb29Ki0Fvp/dbvZtVq68X9DE2ywidgBfAbZTCoA3gCeA1yOiL2vWBYzcaHI5OQiG5/yIOBu4CLhK0u8B3wR+E5hLaSN/tYH1DcdY4GzgmxFxFvAWcNjQ4U1osH41+/b6D9nhrsXAnY2uZSRV6VdTb7MsuJYAJwG/ARxL6X9Ipbofr3cQDENE7Mx+vwL8EJgfEbsi4mBE9APfpjT6ajPpAroi4rFs+i5K/0B3SToeIPv9SoPqG66q/WqB7VXuIuDJiNiVTTf7NhtwSL9aYJtdCLwQEa9GRC/wj8D7gCnZoSKoPhRP4RwENZJ0rKTJA7eBDwJPD7zwMh8Gnm5EfcMVEb8AXpJ0ajbrAmArhw4D8nHgngaUN2yD9avZt1eFyzn08ElTb7Myh/SrBbbZduBcSRMliV+9xh6iNMQONGh7+aqhGkk6mdJeAJQOO3wvIv5K0u2UdlkDeBH4k4HjtM1C0lzgFuAY4HlKV2m0AXcAsyn9IV8WEU01MOAg/fo6Tb69ACRNBF4CTo6IN7J502j+bVatX63wGvtLSl++1Qc8BSyjdE5gNfDObN7HImLkv4ZsqLocBGZmafOhITOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxBX55fVmdZVdNvmjbPLXgYOUhpeA0of+DjSksCFI+iNgXfZ5B7OG8OWj1pJG02iwksZExMFBlj0MrIiITTWsb2zZ2DRmR82HhiwJkj4u6fFsHPu/k9Qmaayk1yX9taQnJd0r6RxJP5b0/MB495KWSfphtnybpC/mXO+XJD0OzJf0l5I2ZOPQf0slH6X0AakfZPc/RlKXpCnZus+V9EB2+0uSbpZ0P6UB9MZK+lr22JslLav/s2qtwkFgLU/S6ZSGJHhfRMyldEh0abb4OOC+bBDBA8D1lD76fxmwsmw187P7nA38F0lzc6z3yYiYHxGPAH8bEe8FzsiWLYqIHwCbgI9GxNwch67OAj4UEVcAy4FXImI+8F5Kgx/OHs7zY+ZzBJaCCyn9s9xYGuKFdkrDFwDsi4j7s9tbgDciok/SFuDEsnXcGxG7ASTdDfwOpdfPYOs9wK+GIgG4QNLngQnAdErDD/9zjf24JyLezm5/EHiPpPLgmUNpSAmzmjgILAUCbo2I/37IzNKIj+XvwvuB/WW3y18flSfT4gjr3RfZCbhs3JxvUBr1dIekL1EKhGr6+NWeemWbtyr69KmI+BFmR8mHhiwFDwAfkTQdSlcXDeMwygdV+u7jiZTGlP/XGtbbTilYXstGrr2kbNkeYHLZ9IvAvOx2ebtK9wKfGhi+WKXv+W2vsU9mgPcILAERsSUb9fEBSW1AL/BfqW3c94eB71H6YpTbB67yybPeiOhW6XuTnwZ+DjxWtvg7wC2S9lE6D3E98G1JvwAeH6KemymNLropOyz1CqWAMquZLx81O4LsipzTI+Izja7FrAg+NGRmljjvEZiZJc57BGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmifv/4gtsoMkmdNcAAAAASUVORK5CYII=\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": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mag/miniconda3/envs/mooc-rr-jupyter/lib/python3.7/site-packages/ipykernel_launcher.py:7: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " import sys\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Wed, 22 Apr 2020 Deviance: 3.0144
Time: 11:29:55 Pearson chi2: 5.00
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
" + ], + "text/plain": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: Frequency No. Observations: 23\n", + "Model: GLM Df Residuals: 21\n", + "Model Family: Binomial Df Model: 1\n", + "Link Function: logit Scale: 1.0000\n", + "Method: IRLS Log-Likelihood: -3.9210\n", + "Date: Wed, 22 Apr 2020 Deviance: 3.0144\n", + "Time: 11:29:55 Pearson chi2: 5.00\n", + "No. Iterations: 6 \n", + "Covariance Type: nonrobust \n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", + "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import statsmodels.api as sm\n", + "\n", + "data[\"Success\"]=data.Count-data.Malfunction\n", + "data[\"Intercept\"]=1\n", + "\n", + "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", + " family=sm.families.Binomial(sm.families.links.logit)).fit()\n", + "\n", + "logmodel.summary()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mag/miniconda3/envs/mooc-rr-jupyter/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " \n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Wed, 22 Apr 2020 Deviance: 18.086
Time: 11:29:55 Pearson chi2: 30.0
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" + ], + "text/plain": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: Frequency No. Observations: 23\n", + "Model: GLM Df Residuals: 21\n", + "Model Family: Binomial Df Model: 1\n", + "Link Function: logit Scale: 1.0000\n", + "Method: IRLS Log-Likelihood: -23.526\n", + "Date: Wed, 22 Apr 2020 Deviance: 18.086\n", + "Time: 11:29:55 Pearson chi2: 30.0\n", + "No. Iterations: 6 \n", + "Covariance Type: nonrobust \n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", + "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 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", + " 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": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAZ1klEQVR4nO3df3BVdX7/8eebBJbwu4BLxaCwLcW1/gASgprv2uAq4M6K2KLIWrdul2W/07LWr5WOTO3qWp1pG6d1u7VWqtRWRwM6guwO3aCWdFvHH4ENgkAD6LKasFsU5UfcAEl4f/84J/ESbpKbm3uTez+8HjMZ7jn3c875vHPI69587jmfmLsjIiL5b9BAd0BERDJDgS4iEggFuohIIBToIiKBUKCLiARCgS4iEogeA93MVpvZQTN7p4vnzcz+3sz2mdl2M5uZ+W6KiEhPUnmH/hQwv5vnrwOmxl/LgMf63i0REemtHgPd3X8CfNxNkxuAf/PIG8AYMzs3Ux0UEZHUFGZgH+cBHyQsN8TrftG5oZktI3oXT1FRUcmkSZN6fbCPjzsn2hxLr685x0G15JhQ6gDVkqsGD4JxRel9hLlnz56P3P2cZM9lItCTfY+Tzifg7quAVQClpaW+ZcuWtA5YU1NDRUVFWtvmGtWSe0KpA1RLrupLLWb2866ey8RVLg1A4lvtYuBABvYrIiK9kIlA3wB8Pb7a5XLgiLufMdwiIiLZ1eOQi5k9B1QA482sAbgPGAzg7v8EbAS+AuwDfgV8I1udFRGRrvUY6O6+pIfnHfjjjPVIRPJCS0sLDQ0NHD9+vF+ON3r0aHbv3t0vx8q2VGoZOnQoxcXFDB48OOX9ZuJDURE5CzU0NDBy5EgmT56MWfavPzl27BgjR47M+nH6Q0+1uDuHDh2ioaGBKVOmpLxf3fovImk5fvw448aN65cwP9uYGePGjev1bz8KdBFJm8I8e9L53irQRUQCoTF0EclbBQUFXHLJJR3L69evZ/LkyQPXoQGmQBeRvFVUVMS2bdu6fL61tZXCwrMn5jTkIiJBeeqpp7jpppu4/vrrmTt3LgCVlZXMmjWLSy+9lPvuu6+j7UMPPcS0adO45pprWLJkCQ8//DAAFRUVtE9N8tFHH3W8629ra2PFihUd+3r88ceBz27lX7RoERdeeCG33nor0RXdUFtby5VXXslll11GWVkZx44dY968eae9EJWXl7N9+/Y+1372vHSJSNZ874c72XXgaEb3edHEUdx3/W9326a5uZnp06cDMGXKFNatWwfA66+/zvbt2xk7diybNm1i7969vPXWW7g7CxYs4Cc/+QnDhw+nqqqKuro6WltbmTlzJiUlJd0e78knn2T06NHU1tZy4sQJysvLO1406urq2LlzJxMnTqS8vJzXXnuNsrIyFi9ezJo1a5g1axZHjx6lqKiIr3/96zz11FM88sgj7NmzhxMnTnDppZf2+XumQBeRvNXVkMu1117L2LFjAdi0aRObNm1ixowZADQ1NbF3716OHTvGjTfeyLBhwwBYsGBBj8fbtGkT27dv54UXXgDgyJEj7N27lyFDhlBWVkZxcTEA06dPZ//+/YwePZpzzz2XWbNmATBq1CgAbrzxRsrLy6msrGT16tXcfvvtfftGxBToItJnPb2T7m/Dhw/veOzurFy5km9/+9untXnkkUe6vDSwsLCQU6dOAZx2Lbi784Mf/IB58+ad1r6mpobPfe5zHcsFBQW0trbi7kmPMWzYMK699lpeeukl1q5dS7ozz3amMXQRCdq8efNYvXo1TU1NADQ2NnLw4EGuuuoq1q1bR3NzM8eOHeOHP/xhxzaTJ09m69atAB3vxtv39dhjj9HS0gLAnj17+PTTT7s89oUXXsiBAweora0FojtEW1tbAVi6dCl33HEHs2bN6vhtoq/0Dl1EgjZ37lx2797NFVdcAcCIESN45plnmDlzJosXL2b69OlccMEFfOlLX+rY5u677+bmm2/m6aef5uqrr+5Yv3TpUvbv38/MmTNxd8455xzWr1/f5bGHDBnCmjVr+M53vkNzczNFRUW88sorAJSUlDBq1Ci+8Y0Mzmfo7gPyVVJS4unavHlz2tvmGtWSe0Kpwz27tezatStr+07m6NGjWd3/fffd55WVlVk9RrujR496Y2OjT5061dva2rpsl+x7DGzxLnJVQy4iIv3s2WefZfbs2Tz00EMMGpS5GNaQi4gIcP/99/fbsb72ta+d8SFtJugduoikzT3pnw+WDEjne6tAF5G0DB06lEOHDinUs8Dj+dCHDh3aq+005CIiaSkuLqahoYEPP/ywX453/PjxXgdcrkqllva/WNQbCnQRScvgwYN79dd0+qqmpqbjbs98l61aNOQiIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEIqVAN7P5ZlZvZvvM7J4kz59vZpvNrM7MtpvZVzLfVRER6U6PgW5mBcCjwHXARcASM7uoU7N7gbXuPgO4BfjHTHdURES6l8o79DJgn7u/5+4ngSrghk5tHBgVPx4NHMhcF0VEJBXW01/sNrNFwHx3Xxov3wbMdvflCW3OBTYBvwYMB65x961J9rUMWAYwYcKEkqqqqrQ63dTUxIgRI9LaNteoltwTSh2gWnJVX2qZM2fOVncvTfqku3f7BdwEPJGwfBvwg05t7gL+NH58BbALGNTdfktKSjxdmzdvTnvbXKNack8odbirllzVl1qALd5FrqYy5NIATEpYLubMIZVvAmvjF4jXgaHA+BT2LSIiGZJKoNcCU81sipkNIfrQc0OnNu8DXwYwsy8SBfqHmeyoiIh0r8dAd/dWYDlQDewmupplp5k9YGYL4mZ/CnzLzN4GngNuj381EBGRflKYSiN33whs7LTuuwmPdwHlme2aiIj0hu4UFREJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQKQW6mc03s3oz22dm93TR5mYz22VmO83s2cx2U0REelLYUwMzKwAeBa4FGoBaM9vg7rsS2kwFVgLl7v6JmX0+Wx0WEZHkUnmHXgbsc/f33P0kUAXc0KnNt4BH3f0TAHc/mNluiohIT8zdu29gtgiY7+5L4+XbgNnuvjyhzXpgD1AOFAD3u/uPk+xrGbAMYMKECSVVVVVpdbqpqYkRI0aktW2uUS25J5Q6QLXkqr7UMmfOnK3uXprsuR6HXABLsq7zq0AhMBWoAIqB/zKzi9398Gkbua8CVgGUlpZ6RUVFCoc/U01NDelum2tUS+4JpQ5QLbkqW7WkMuTSAExKWC4GDiRp85K7t7j7z4B6ooAXEZF+kkqg1wJTzWyKmQ0BbgE2dGqzHpgDYGbjgd8C3stkR0VEpHs9Brq7twLLgWpgN7DW3Xea2QNmtiBuVg0cMrNdwGZghbsfylanRUTkTKmMoePuG4GNndZ9N+GxA3fFXyIiMgB0p6iISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEIqVAN7P5ZlZvZvvM7J5u2i0yMzez0sx1UUREUtFjoJtZAfAocB1wEbDEzC5K0m4kcAfwZqY7KSIiPUvlHXoZsM/d33P3k0AVcEOSdn8J/A1wPIP9ExGRFJm7d9/AbBEw392Xxsu3AbPdfXlCmxnAve7+e2ZWA9zt7luS7GsZsAxgwoQJJVVVVWl1uqmpiREjRqS1ba5RLbknlDpAteSqvtQyZ86cre6edFi7MIXtLcm6jlcBMxsE/B1we087cvdVwCqA0tJSr6ioSOHwZ6qpqSHdbXONask9odQBqiVXZauWVIZcGoBJCcvFwIGE5ZHAxUCNme0HLgc26INREZH+lUqg1wJTzWyKmQ0BbgE2tD/p7kfcfby7T3b3ycAbwIJkQy4iIpI9PQa6u7cCy4FqYDew1t13mtkDZrYg2x0UEZHUpDKGjrtvBDZ2WvfdLtpW9L1bIiLSW7pTVEQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkECld5SKSLevrGqmsrufA4WYmjilixbxpLJxx3kB3S1Kk85dbFOgyYNbXNbLyxR00t7QB0Hi4mZUv7gBQKOQBnb/coyEXGTCV1fUdYdCuuaWNyur6AeqR9IbOX+5RoMuAOXC4uVfrJbfo/OUeBboMmIljinq1XnKLzl/uUaDLgFkxbxpFgwtOW1c0uIAV86YNUI+kN3T+co8+FJUB0/7Bma6SyE86f7lHgS4DauGM8xQAeUznL7doyEVEJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBApBbqZzTezejPbZ2b3JHn+LjPbZWbbzexVM7sg810VEZHu9BjoZlYAPApcB1wELDGzizo1qwNK3f1S4AXgbzLdURER6V4q79DLgH3u/p67nwSqgBsSG7j7Znf/Vbz4BlCc2W6KiEhPzN27b2C2CJjv7kvj5duA2e6+vIv2/wD80t0fTPLcMmAZwIQJE0qqqqrS6nRTUxMjRoxIa9tco1pyTyh1gGrJVX2pZc6cOVvdvTTZc4UpbG9J1iV9FTCz3wdKgd9J9ry7rwJWAZSWlnpFRUUKhz9TTU0N6W6ba1RL7gmlDlAtuSpbtaQS6A3ApITlYuBA50Zmdg3w58DvuPuJzHRPRERSlcoYei0w1cymmNkQ4BZgQ2IDM5sBPA4scPeDme+miIj0pMdAd/dWYDlQDewG1rr7TjN7wMwWxM0qgRHA82a2zcw2dLE7ERHJklSGXHD3jcDGTuu+m/D4mgz3SyQt6+saqayu58DhZiaOKWLFvGkAZ6xbOOO8fjl2No6TinvX7+C5Nz/gzotb+ObKjSyZPYkHF14yIH2R/pNSoIvkg/V1jax8cQfNLW0ANB5uZsXzb4NBS5t3rFv54g6AjIZtsmNn4zipuHf9Dp554/2O5Tb3jmWFeth0678Eo7K6viNQ27Wc8o4wb9fc0kZldX3Wj52N46TiuTc/6NV6CYcCXYJx4HBzVtr2ZX+ZPk4q2rq4t6Sr9RIOBboEY+KYoqy07cv+Mn2cVBRYsltHul4v4VCgSzBWzJtG0eCC09YNHmQMLjg9yIoGF3R8WJrNY2fjOKlYMntSr9ZLOPShqASj/cPHgbjKpatjD8RVLu0ffLaPmReY6SqXs4QCXYKycMZ5SUO0P4K1q2MPhAcXXsKDCy+hpqaGd2+tGOjuSD/RkIuISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhKIwlQamdl84PtAAfCEu/9Vp+c/B/wbUAIcAha7+/7MdlUkXOvrGqmsrufA4WYmjilixbxpPL/lfV579+OONuW/MZabSs8/ox1wxrotP/+Y5978gDsvbuGbKzeyZPYkHlx4SUrHTba/hTPOS7nf7cduc6fALCvHTrZtV308m/QY6GZWADwKXAs0ALVmtsHddyU0+ybwibv/ppndAvw1sDgbHRYJzfq6Rla+uIPmljYAGg83c+eabWe0e+3dj08L+MbDzax44W1waDnlHevuWrONUwnbtbnzzBvvA5wWrMmOu+L5t8Ggpe2z/a18cQfAGYGZbPv+OHaybbvq49kmlSGXMmCfu7/n7ieBKuCGTm1uAP41fvwC8GUzs8x1UyRcldX1HeHUWy1t3hHm7U510fa5Nz/o8bgtp7wjUNs1t7RRWV1/xv6Sbd8fx062bVd9PNuYu3ffwGwRMN/dl8bLtwGz3X15Qpt34jYN8fK7cZuPOu1rGbAsXpwGpHsGxgMf9dgqP6iW3NOvdQz59d8syda+2351hIJhozuWT/5y39Z0j5u4bV+3T3Pb8cBH3W3buY85rC//xy5w93OSPZHKGHqyd9qdXwVSaYO7rwJWpXDM7jtktsXdS/u6n1ygWnJPKHVAVEvrkYPB1BLSeclGLakMuTQAkxKWi4EDXbUxs0JgNPAxIiLSb1IJ9FpgqplNMbMhwC3Ahk5tNgB/ED9eBPyH9zSWIyIiGdXjkIu7t5rZcqCa6LLF1e6+08weALa4+wbgSeBpM9tH9M78lmx2mgwM2+QQ1ZJ7QqkDVEuuykotPX4oKiIi+UF3ioqIBEKBLiISiJwPdDMbamZvmdnbZrbTzL4Xr59iZm+a2V4zWxN/YJvzzKzAzOrM7Efxcr7Wsd/MdpjZNjPbEq8ba2Yvx7W8bGa/NtD9TIWZjTGzF8zsf8xst5ldkY+1mNm0+Hy0fx01szvztJb/F/+8v2Nmz8U5kK8/K38S17HTzO6M12XlnOR8oAMngKvd/TJgOjDfzC4nml7g79x9KvAJ0fQD+eBPgN0Jy/laB8Acd5+ecD3tPcCrcS2vxsv54PvAj939QuAyovOTd7W4e318PqYTzav0K2AdeVaLmZ0H3AGUuvvFRBdjtE8pklc/K2Z2MfAtojvuLwO+amZTydY5cfe8+QKGAT8FZhPdZVUYr78CqB7o/qXQ/+L45F0N/Ijohqy8qyPu635gfKd19cC58eNzgfqB7mcKdYwCfkZ8gUA+19Kp/3OB1/KxFuA84ANgLNGVeD8C5uXjzwpwE9GEhu3LfwH8WbbOST68Q28fptgGHAReBt4FDrt7a9ykgeg/Qa57hOhktk95MY78rAOiO4E3mdnWeEoHgAnu/guA+N/PD1jvUvcF4EPgX+KhsCfMbDj5WUuiW4Dn4sd5VYu7NwIPA+8DvwCOAFvJz5+Vd4CrzGycmQ0DvkJ0E2ZWzkleBLq7t3n0a2Qx0a8uX0zWrH971Ttm9lXgoLsnzjWR0pQJOarc3WcC1wF/bGZXDXSH0lQIzAQec/cZwKfk+JBET+Kx5QXA8wPdl3TE48k3AFOAicBwov9nneX8z4q77yYaKnoZ+DHwNtDa7UZ9kBeB3s7dDwM1wOXAmHiaAUg+HUGuKQcWmNl+ohkrryZ6x55vdQDg7gfifw8SjdOWAf9rZucCxP8eHLgepqwBaHD3N+PlF4gCPh9raXcd8FN3/994Od9quQb4mbt/6O4twIvAleTvz8qT7j7T3a8iuvFyL1k6Jzkf6GZ2jpmNiR8XEZ3s3cBmomkGIJp24KWB6WFq3H2luxe7+2SiX4f/w91vJc/qADCz4WY2sv0x0XjtO5w+BURe1OLuvwQ+MLNp8aovA7vIw1oSLOGz4RbIv1reBy43s2HxNNzt5yTvflYAzOzz8b/nA79LdG6yck5y/k5RM7uUaK71AqIXoLXu/oCZfYHone5YoA74fXc/MXA9TZ2ZVQB3u/tX87GOuM/r4sVC4Fl3f8jMxgFrgfOJfihvcvecn6TNzKYDTwBDgPeAbxD/XyP/ahlG9IHiF9z9SLwu785LfHnyYqLhiTpgKdGYeV79rACY2X8RfV7WAtzl7q9m65zkfKCLiEhqcn7IRUREUqNAFxEJhAJdRCQQCnQRkUAo0EVEApHKH4kW6VfxJV2vxou/DrQR3Z4PUObuJwekY90wsz8ENsbXtYsMCF22KDnNzO4Hmtz94RzoS4G7t3Xx3H8Dy919Wy/2V5gwN4lIn2nIRfKKmf2BRfPjbzOzfzSzQWZWaGaHzazSzH5qZtVmNtvM/tPM3jOzr8TbLjWzdfHz9WZ2b4r7fdDM3gLKzOx7ZlYbz2/9TxZZTDS185p4+yFm1pBwh/PlZvZK/PhBM3vczF4mmhCs0Mz+Nj72djNb2v/fVQmFAl3yRjy39I3AlfFkbYV89gfJRwOb4gnDTgL3E90yfhPwQMJuyuJtZgJfM7PpKez3p+5e5u6vA99391nAJfFz8919DbANWOzRfOQ9DQnNAK5399uAZUSTtpUBs4gmOjs/ne+PiMbQJZ9cQxR6W6IpPigius0doNndX44f7wCOuHurme0AJifso9rdPwEws/XA/yH6Oehqvyf5bJoDgC+b2QpgKDCeaFrXf+9lHS+5+/H48Vzgi2aW+AIyleh2cJFeUaBLPjFgtbv/xWkroxn4Et8VnyL6S1ftjxP/n3f+0Mh72G+zxx80xfOk/AMw090bzexBomBPppXPfgPu3ObTTjX9kbu/ikgfachF8skrwM1mNh6iq2HSGJ6Ya9HfEB1GNOf2a73YbxHRC8RH8WyTv5fw3DFgZMLyfqI/A0endp1VA3/UPi2sRX8XtKiXNYkAeocuecTdd8Sz8L1iZoOIZq/7v/RuXuz/Bp4FfgN4uv2qlFT26+6HzOxfiaYK/jnwZsLT/wI8YWbNROP09wP/bGa/BN7qpj+PE824ty0e7jlI9EIj0mu6bFHOGvEVJBe7+50D3ReRbNCQi4hIIPQOXUQkEHqHLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiP8PfJgpKZPR3KwAAAAASUVORK5CYII=\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": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAENCAYAAAARyyJwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxTVd4/8M+92dp0oVvaYlkEFRBoi1ZlkWVwpGVHAR3Bnx03HMdnxOH1PK6o4wqiPoMyo46oLx0VRhh1WHSmoPLgCO0oi1L2pYBl6ZLuW5J7c+/5/ZE2WgshLY1N0s/79dI2uUn6PWT55Jx77rmSEEKAiIjoLOSuLoCIiIIbg4KIiHxiUBARkU8MCiIi8olBQUREPjEoiIjIp4AHRUNDA6ZOnYqTJ0+22bZ//37MnDkTOTk5WLhwIdxud6DLISKidgpoUOzatQtz5szB8ePHz7j9/vvvx+OPP44NGzZACIHVq1cHshwiIuqAgAbF6tWr8Yc//AHJyclttp06dQpOpxPDhg0DAMycORN5eXmBLIeIiDrAGMgHf/bZZ8+6rby8HDabzXvZZrOhrKwskOUQEVEHdNnObF3XIUmS97IQotVlIiIKDgHtUfiSmpoKu93uvVxRUXHGISpfqqsboevhuVRVYmI0KisburqMgGH7Qlc4tw0I7/bJsoT4+Kh236/LgiItLQ0WiwU7duxAVlYW1q5di7Fjx7brMXRdhG1QAAjrtgFsXygL57YB4d++9vrZh57mzZuH3bt3AwBefPFFLF68GBMnTkRTUxNyc3N/7nKIiOgcpFBeZryysiFsk99mi4HdXt/VZQQM2xe6wrltQHi3T5YlJCZGt/9+AaiFiIjCCIOCiIh8YlAQEZFPDAoiIvKJQUFERD4xKIiIyCcGBRER+cSgICIinxgURETkE4OCiIh8YlAQEZFPDAoiIvKJQUFERD4xKIiIyCcGBRER+cSgICIinxgURETkE4OCiIh8YlAQEZFPDAoiIvKJQUFERD4xKIiIyCcGBRER+cSgICIinxgURETkE4OCiIh8YlAQEZFPDAoiIvKJQUFERD4xKIiIyCcGBRER+cSgICIinxgURETkE4OCiIh8YlAQEZFPAQ2K9evXY/LkycjOzsaKFSvabN+7dy9mzZqF6dOn4ze/+Q3q6uoCWQ4REXVAwIKirKwMS5cuxcqVK7FmzRqsWrUKR44caXWbZ599FvPnz8e6devQr18/vPXWW4Eqh4iIOihgQZGfn48RI0YgLi4OVqsVOTk5yMvLa3UbXdfR2NgIAHA4HIiIiAhUOURE1EHGQD1weXk5bDab93JycjIKCwtb3eahhx7C7bffjkWLFiEyMhKrV69u199ITIzulFqDlc0W09UlBBTbF7rCuW1A+LevvQIWFLquQ5Ik72UhRKvLTqcTCxcuxDvvvIOMjAy8/fbbePDBB7F8+XK//0ZlZQN0XXRq3cHCZouB3V7f1WUEDNsXusK5bUB4t0+WpQ59wQ7Y0FNqairsdrv3st1uR3JysvfyoUOHYLFYkJGRAQD41a9+hW+++SZQ5RARUQcFLChGjRqFgoICVFVVweFwYOPGjRg7dqx3e9++fVFaWoqjR48CAL744gukp6cHqhwiIuqggA09paSkYMGCBcjNzYWqqpg9ezYyMjIwb948zJ8/H+np6Vi8eDF+//vfQwiBxMRELFq0KFDlEBFRB0lCiJAd5Oc+itDF9oWucG4bEN7tC7p9FEREFB4YFERE5BODgoiIfGJQEBGRTwwKIiLyiUFBREQ+MSiIiMgnBgUREfnEoCAiIp8YFERE5BODgoiIfGJQEBGRTwwKIiLyiUFBREQ+MSiIiMgnBgUREfnEoCAiIp8YFERE5BODgoiIfGJQEBGRTwwKIiLyiUFBREQ+MSiIiMgnBgUREfnEoCAiIp8YFERE5BODgoiIfGJQEBGRTwwKIiLyiUFBREQ+MSiIiMgnBgUREfnEoCAiIp/8Cor33nsPDQ0Nga6FiIiCkF9BcfDgQeTk5GDhwoXYvXu33w++fv16TJ48GdnZ2VixYkWb7UePHsUtt9yC6dOn44477kBtba3/lRMR0c/Cr6B45plnsGHDBgwdOhRPPvkkZs2ahQ8//BAul+us9ykrK8PSpUuxcuVKrFmzBqtWrcKRI0e824UQ+O1vf4t58+Zh3bp1uPTSS7F8+fLzbxEREXUqv/dRREdHY+LEiZg6dSpqamqwcuVKTJw4EZs2bTrj7fPz8zFixAjExcXBarUiJycHeXl53u179+6F1WrF2LFjAQB33303br755vNsDhERdTajPzcqKCjAqlWrUFBQgJycHLzyyisYNGgQiouLMXfuXFxzzTVt7lNeXg6bzea9nJycjMLCQu/l4uJiJCUl4ZFHHsH+/fvRv39/PPbYY53QJCIi6kx+BcWTTz6JuXPn4umnn0ZMTIz3+j59+uDGG2884310XYckSd7LQohWl91uN7755hu8//77SE9Px0svvYTnnnsOzz33nN/FJyZG+33bUGSzxZz7RiGM7Qtd4dw2IPzb115+BcW6deuQl5eHmJgY2O12fPrpp8jNzYUsy5g/f/4Z75Oamort27d7L9vtdiQnJ3sv22w29O3bF+np6QCAqVOnnvWxzqaysgG6Ltp1n1Bhs8XAbq/v6jIChu0LXeHcNiC82yfLUoe+YPu1j+Lpp5/G5s2bm/+QjB07dmDRokU+7zNq1CgUFBSgqqoKDocDGzdu9O6PAIDLLrsMVVVVOHDgAABg06ZNGDJkSLsbQEREgeVXj+Lbb7/FJ598AgBITEzEyy+/jBkzZvi8T0pKChYsWIDc3FyoqorZs2cjIyMD8+bNw/z585Geno5XXnkFjz76KBwOB1JTU/H888+ff4uIiKhT+RUUqqpCURSYzWYAnv0L/pg2bRqmTZvW6ro33njD+3tmZiY+/PBDf2slIqIu4FdQ/OIXv8Add9yBGTNmQJIkfPLJJxg3blygayMioiDgV1A88MADWLFiBb744gsYjUZMmDABN910U6BrIyKiIOBXUBgMBuTm5iI3NzfQ9RARUZDxKyg+//xzLFq0CLW1tRDih+moO3fuDFhhREQUHPwKihdeeAEPPfQQBg8e3OqgOSIiCn9+BUVsbCyys7MDXQsREQUhvw64y8zMxJdffhnoWoiIKAj51aP48ssv8f7778NkMsFkMnnXbeI+CiKi8OdXULzzzjsBLoOIiIKVX0NPaWlp2L17N1avXo2EhAR8++23SEtLC3RtREQUBPwKiuXLl+Nvf/sb8vLy4HQ68ec//xmvvPJKoGsjIqIg4FdQfPrpp3jjjTcQGRmJ+Ph4rF692rtIIBERhTe/gsJoNHoXBAQ802WNRr92bxARUYjz69O+Z8+e2Lx5MyRJgqIoeOutt7iPgoiom/ArKB577DE88MADOHjwIIYNG4bMzEy8+OKLga6NiIiCgF9BkZKSgr/+9a9wOBzQNA3R0eF9rmoiIvqBX0Hx9ttvn/H62267rVOLISKi4ONXUBw6dMj7u6Io2LZtG0aOHBmwooiIKHj4FRSLFy9udbmsrAwLFy4MSEFERBRc/Joe+1MpKSk4depUZ9dCRERBqN37KIQQ2LNnDxITEwNWFBERBY9276MAPMdVPPDAAwEpiIiIgkuH9lEQEVH34VdQ3HLLLT5Pgfruu+92WkFERBRc/AqKoUOHoqioCDfeeCNMJhPWrl0Lt9uNKVOmBLo+IiLqYn4Fxc6dO7Fy5UoYDAYAwJgxY3DjjTciJycnoMUREVHX82t6bFVVFVwul/dyY2MjnE5nwIoiIqLg4VePYurUqfjVr36FCRMmQAiBf/3rX8jNzQ10bUREFAT8Cor77rsPgwcPxn/+8x9YLBY89dRTuOqqqwJdGxERBQG/j8xOSUnBJZdcgt///vcwmUyBrInonCQJEBBdXQZRt+BXUHz00Ud4+OGH8eabb6K+vh733HMPVq9eHejaiHxS3XpXl0DULfgVFO+//z5WrVqF6OhoJCYm4uOPP8Zf//rXQNdG5JNT0eHj8B4i6iR+BYUsy61OVtSzZ0/vVFmiruLWdGg6h5+IAs2voIiLi8P+/fu9R2evW7cOPXr0CGhhROei6wIKh5+IAs6vWU+PPPII7rvvPhQXF2P06NGwWCx49dVXA10b0TkIOJxuRJoNEOxYEAWMX0HhdDqxdu1aHD9+HJqmoV+/fn7NfFq/fj1ee+01uN1u/PrXv8bNN998xttt3rwZTz31FDZt2tS+6qnbU9wa3JoOg9yhU6sQkR/8enf9z//8DwwGAy666CIMGDDAr5AoKyvD0qVLsXLlSqxZswarVq3CkSNH2tyuoqICS5YsaX/lRACEAJyK1tVlEIU1v4Ji4MCBWL9+PU6fPo2amhrvf77k5+djxIgRiIuLg9VqRU5ODvLy8trc7tFHH8Xvfve7jlVPBMDhcvOICqIA8mvo6YsvvmjzIS9JEvbv33/W+5SXl8Nms3kvJycno7CwsNVt3n33XQwePBiZmZntqdkrMTH63DcKYTZbTFeXEFDn2z43JGjNCREVY0FURHAdCBrOz184tw0I//a1l19BsXv37nY/sK7rrc5hIYRodfnQoUPYuHEj3nnnHZSWlrb78QGgsrIBephOj7TZYmC313d1GQFzvu2TJKC6xgF3c1I01juREGsJmp3a4fz8hXPbgPBunyxLHfqC7XPo6bHHHvP+XlVV1a4HTk1Nhd1u91622+1ITk72Xs7Ly4PdbsesWbNw1113oby8HHPnzm3X3yBqobg1uFROlSUKBJ9BsWfPHu/vd9xxR7seeNSoUSgoKEBVVRUcDgc2btyIsWPHerfPnz8fGzZswNq1a7F8+XIkJydj5cqV7SyfyEMIoMmp8khtogDwGRTiR/140c4+fUpKChYsWIDc3Fxcd911mDp1KjIyMjBv3rwODWURnYtL1TgDiigA/NpHAcDnObPPZtq0aZg2bVqr69544402t+vVqxePoaDzJgRQ36jA1CMCMrsWRJ3GZ1Douo7a2loIIaBpmvf3FnFxcQEvkKg93LpAg0NFrNXc1aUQhQ2fQXHo0CGMGDHCGw7Dhw/3bjvX9FiiruJwuRFhNsJs5NHaRJ3BZ1AcOHDg56qDqNMIAdQ3uZAQE8md20SdgF+5KCypbgGH6u7qMojCAoOCwlZTkwo9WI7AIwphDAoKW25dwOFir4LofPk9PZYoWBQWVeDL706jstYJo0HCmMwLMLBP/Blv2+h0I9Ji5HTZEFRYVIG8r4tRUetEUo8ITBzeBxkXJXV1Wd0SexQUUgqLKrDis0OobVJgsRhQ51CxbusxHCyuPuPtdfYqQlLL81zTqMAaYURNo4IVnx1CYVFFV5fWLTEoKKTkfV0Mg0GGxWSABAlmowEGg4yvdp0+630anW6eWzvEtHqeJQkWk+d5zvu6uKtL65YYFBRSKmqdbY6PMBlkVNe7znofXReodygAz1oRMs70PJuNMipqnV1UUffGoKCQktQjAoq79SqxqqYjPsbi835Ol4YmrgMVMs70PCtuHUk9Irqoou6NQUEhZeLwPtA0HS5Vg4CA4tagaTrGZF5wzvs2NKlwcwgqJLR6noWAS/U8zxOH9+nq0rolznqikNIy66Vl1lNspMnnrKcf03WBukYF8dEWHrEd5FqeZ856Cg4MCgo5GRclIfPiJNh/dIY7fymqBofqhtXMl36wy7goicEQJDj0RN1Ok0Pt6hKIQgqDgrodtybQpPDYCiJ/MSioW2pqYq+CyF8MCuqW3Dp7FUT+YlBQt8V9FUT+YVBQt+XWBBzsVRCdE4OCujWHSwOX9iDyjUFB3Zrq1qC281gMou6GQUHdmhCAk2tAEfnEoKCQU9uo4JP84yivdnTK4zlcXIacyBcGBYWcf393Ch99eRR//ng3Dp+sOe/H03WBRidnQBGdDYOCQs7wwSmIizbDqWj4678O4Jv9Zef9mA6XG6qmn/uGRN0Qg4JCTnK8FY/9+gr0TLRCF8Car47h04Lj0M9j+EgIoJ5HaxOdEYOCQlJCbAR+e91QDOoTBwDYursU7204COd5HBehqBrqmhQuQU70EwwKClkRZgP+X/ZAjMnoCQA4eKIGr63Zi8rzOF2mw+XmLCiin2BQUEiTZQmTRvTFzLH9YZAl2GsceHVNx3dye4agFAgehEfkxaCgsHDFoGTcOXUwoiJNcLg0vPOvA/j3d6chRPs/8N2aQIODS3sQtWBQUNjomxqD/7p+KNJsURACyPumGH/7/DBcHRhKanKqUFTOgiICGBQUZuKiLbhr2hBkDbABAPYcq8Kra3ajrLqpXY8jBFDX5ILegR4JUbhhUFDYMRllzBzXHzNG92veb+HEa//Yg++OVLTrcdyaQE29C7rOngV1bwENivXr12Py5MnIzs7GihUr2mz//PPPMWPGDEyfPh333HMPamtrA1kOdSOSJGH44BTcNX0wekSZobh1rN50BP/491Gobv8/+BW3jqp6F5R23Ico3AQsKMrKyrB06VKsXLkSa9aswapVq3DkyBHv9oaGBjzxxBNYvnw51q1bh4EDB+JPf/pToMqhbqp3cgzunZWOAb09x1tsO1CO19bsQXmN/+tEtfQsPGfE41AUdT8BC4r8/HyMGDECcXFxsFqtyMnJQV5enne7qqr4wx/+gJSUFADAwIEDUVJSEqhyqBuzRpiQO3Egcq7qDVkCSqua8MrHu7H9QLnfs6J0IVDfoKCuSeXUWep2AhYU5eXlsNls3svJyckoK/thTZ74+HhMmDABAOB0OrF8+XJce+21gSqHujlZkjBuWBrmTRuCuGgzVLeOj/99FH/74jCanP5NhRUAmpxuVNU5uS4UdSvGQD2wruuQfrQWghCi1eUW9fX1+K//+i8MGjQI119/fbv+RmJi9HnXGcxstpiuLiGgzrd9bkho7zmHEhKiMKB/Ilb86wB2HizHnqNVOGlvxK1TBmPQhQnteixzpAmxUeYzvq6B8H7+wrltQPi3r70CFhSpqanYvn2797LdbkdycnKr25SXl+OOO+7AiBEj8Mgjj7T7b1RWNpzXQnDBzGaLgd1e39VlBMz5tk+SgOoaB9wdPDvdrLH9cGFKND7JP46aehde+uBbjBqaipyr+sBk9K+jXQnAYjIgxmqG0dA6LML5+QvntgHh3T5Zljr0BTtgQ0+jRo1CQUEBqqqq4HA4sHHjRowdO9a7XdM03H333Zg0aRIWLlx41m9lRIEgSRKuGJSM381KR+9kzxsnf08p/vRRIU6U+/8h4VI1VNY5UO9QecwFha2A9ShSUlKwYMEC5ObmQlVVzJ49GxkZGZg3bx7mz5+P0tJS7Nu3D5qmYcOGDQCAoUOH4tlnnw1USURtJPWIxF3Th+DL705h045TqKh14i9r92J0ek9ce0Vvv3oXQgCNDhUuxY3YKAvMfvZIiEKFJDqyGE6Q4NBT6OqMoSf7eQw9nUnBnhJs+OaE95iJWKsJo4am4tCJGlTXuxAfY8GYzAswsE+8z7qsESb0SYtDdVVju2soLKpA3tfFqKh1IqlHBCYO74OMi5I63KbOtG7LUWzcdhJOVUOEyYDsK3th+uj+XV1Wpwvn917QDT0RhZKDxdXYsrsEMVEmREd6Otp1TSryvjmB05VNMJsMqHOoWLf1GA4WV5/1cVp6F/bqJjhVDe057qKwqAIrPjuEmkYF1ggjahoVrPjsEAqL2ndEeSCs23IU6/KPw6VqMMqeIbd1+cexbsvRri6NfgYMCiIAX+06DYNBhsVkRGyUBba4SLTsNXMqGuw1TmiagCxL+GrX6XM+ni6A2noXahoUuP1cAiTv6+LmGgyQJAkWkwEGg4y8r4vPo2WdY+O2k5AgwSBLkCTZ8xMSNm472dWl0c8gYPsoiEJJdb0LEZYf3g4t+yZawkIXAtX1LphNMlQ/V5UV8ISMouqIijQiKsLk8/YVtU5YI1q/Jc1GGRXncSKmzuJU3DDIrSecyBLO64yCFDrYoyACEB9jaXMQndEgwWiQYIuPRITZAABQVN0zJPV1MVyqf8uX60KgvklFRa0DTYr7rEd2J/WIaLOmlOLWkdQjogMt6lwRZiN+ujtQF57rKfwxKIgAjMm8AJqmQ3FrEEJAcWswmwwwm43QhUB8jAWxUSa0fKn+967T+OOq7/DtIbvf02LdmkBdg4LKGieaXG0DY+LwPtA0HS7VU4NL1aBpOiYO79PZzW237Ct7QUBA0wWE0D0/IZB9Za+uLo1+BoYnnnjiia4uoqMcDgWhO2fLt6goC5qalK4uI2DOt32S5FlOo7MmvSX1iERSjwiUVTahoUlFXJQZk0b0xeAL473XJcZGYNKIvkiOt+JEeT2cioZ9x6tx6EQNbPGRiIu2eB8vMtIMh0M949/ShWdnsMulQTbI3mGulAQrUuIjcbK8AbWNChJiLJg5tn9QzHoa2CceEALflzZA0QQiTAZMHtEnLGc9hfN7T5IkWK3m9t+P02ODUzhP0QOCc3pse1TXu5D39ffYfbTKe92QCxOQfVVv2OIikZAQhSo/psdKACwWA6IjzDAZpZD44sPXZujq6PRYDjASdUB8jAVzrh2AkaV1+GfB9zhpb8Te41XY/30VsgYmY9YvB/j1OAKA06XBpThgMRlgjTDBbJTww250oq7HoCA6DxemxuLu64Zid1ElNm47gep6F7YdKMd3RyowfHAKxmZegOhI37OdAM/xF05F8xynYJARGWGExSjDaJBDopdB4Y1BQXSeZElC5sVJGNIvAdv2l2PTt6fQ6FCxpbAE3+wrw8ihqRid0fOc02MBT2Cobh1qgwJZkmAyNoeGycA+BnUZBgVRJzEaZIwcmorLB9rwXVEVNvznOJyKhi+/O42CPaUYMSQFV6f3RIyfOxP15plPLlWDySghwmKCSZZhNHoOfGNPg34uDAqiTmYxGTBp1IXI7B+PrbtLsXV3CZyKhn/vKkH+nlJkDUzG2MyeiI/x//gI1S2guj0zcWRZgtkow2oxwWSS2dOggGNQEAVIhNmIX2b1wtXpqcjfU4r83aVocrnx9b4ybNtfhqH9EzEmoyfSbO2bhaLrwrM/Q9FgMEqIMBthbt6fwZ4GBQKDgijAIsxGXHN5L1yd3hPbD5RjS2EJahsVFBZVorCoEhf2jMHVQ3vi0r7xkGX/+wcCgNst0OD2HK8hyxJMBhlmswEmgwyTUYIsMTjo/DEoiH4mFpMBV6f3xPDBKSgsqsSWwhKUVjXheEk9jpfUIz7GguGDU5A10ObXju+f0nUBl+7ZpyEBkJqHqCxmA8wGGQYDF2KgjmFQEP3MjAYZlw+w4bJLknDkVC3yd5fiYPM5L/K+Lsbn208gvX8irro0BX1Sojt09kcBQDQPUTkVDZIEGAwSLCYjLCYZBlmG0cDeBvmHQUHURSRJwiW94nBJrzhU1Drwn71l2HHQDpeq4dvDFfj2cAWS4yNxxcBkDLskya/jMc5GCM8wldutotHhmdIrGwCLyQiT0RMaBkmGLIPhQW1wCY8gFc7LCAChv4THufi7hMdPuVQNu45U4Jv95Thd8cP9DbKEgX3icPkAGwb0joOxk4eRJMkTXEaDDJPBMwXXaPAEyE/3c/C1Gbq4hAdRGLCYDLjq0hRcdWkKTlU0YvuBcuw6UuFdgHDf8WpYLUakX5SIYRcnoXdKNOQODE39lBDwrJqra1Cal09vFR5GGWaTDKMsQdcFJIk9j+6EPYogFc7fagD2KNpDdevYe7wKOw/aUXS6ttUHdI8oM9IvSkR6/wT0snVsf0Z7SJKnbY11TpiaZ1cZ5JYz34VHeITze489CqIwZTLKGHZxEoZdnIS6RgW7iiqw63AFTlc2obZRwZbCEmwpLEGPKDOG9EvA4AsT0Dc1ps0Z6TqDp+cBOFUNzpbZVZLUvLPc0+MwGCQYDDJkSIDkOROe4QxDWBQ6GBREISQ2yowxGRdgTMYFsNc4UFhUid1HK1Fe7UBto+I5sG9PKawWIwb1jcOgPvG4pFccLM1n6OtsAp4hKwhA0zX89CwOLR0cSZIgS54ZX0ajpxciN593W5Y833RbRgdkWfLeHgiPXkqoY1AQhShbXCR+mdULv8zqhfJqB/Ycq8S+Y1U4XdmEJpcbOw9VYOehChhkCRf2jMGAXnG4pHccUuIjAz5E1aLlQ14IAR2AW9MA5YdTyErN/5Mgec/411KbJ0Bkz8/m3ogsS5BlAM29GBk//C4Jz+8twdO8ybsPh4HTcQwKojCQHB+Ja+J74ZrLe6G63oX931fjwPfVOFZSB00XKDpVh6JTdfjX18WItZpwUVoPXJzWA/3TeqBHVPvPeNZZRPP/fnxa2JbdpjoAaP6dl/zHpObgaVkESwIgGyQYJdnTW5Gbh8uab2uQmns3zTvv3ZoOAdEmWOSWVOuGGBQUskxGGQb5h7euDkAXOnQdbT58gLYfHpBaPlQ83zZ18aP7CSBUv4DGx1gwamgqRg1NhVNx48ipOhw+UYNDJ2pQ26igrkn1HqcBAEk9ItD/glhc2DMW/VJj0ONHp3QNReJHz2ELTRdQPdFzRt4hMkjQZBnVNc422+OjLTAaGBREIUMItDpH9Y/pwrOcxY8n9EnNyfDjoYifjr7ouucDRtcFdPGj3/UfIkfXBFRdh6YJb6gE85BGhNmIof0SMLRfAoQQsNc4cfhkDY6cqsWxkjooqo6KWicqap34Zn85AE/Q9E2JQd/UGPRJiUZKvLVda1CFIu8QWfPz+dPZlD/TSF3QYlBQyDrbB7QENM/48f3u/un9W4Ys5HN8a2yZBqrpojlMWoeKw+X2vxE/I0mSkBwfieT4SFyd3hOaruOUvRFHT9fhWEkdvi+th+LWUV3vQnW9C98d8fQ4zEYZabYo9E6ORi9bNIYYPIdv/1z7Oajr8TiKIBXOc7mB8G6fgIAl0oLKqkZomg7FrbdsgB7EbzdN11FS0YTjpfX4vrQexWX1qHeoZ7ytNcKItKQo9EyMwgVJVvRMjEJibERY9DzOdAyMJAGJsRGdfkT8z43HURAFCQkSekRboDgUSJKn5wF4dtK6VM+wVcv3M8WtBc1BgwZZRq/kaPRKjsbojJ4QQqCmwan4u7wAABFRSURBVIXisgacLG/ACXsDSiqaoGo6mpxuHD5Zi8Mna733NxlkJCdEIjXeipQEK5LjI5ESH4nYKDN7HyGOQUEUQEL8sE8EkgSrpfU3UgFAVTW43DoUVYPWPHwVDCRJQnxMBOJjIpB5cRIAT+gpOrCvyI7TFU04XdGI0qomuFQNquYZyjplb/1t3GIywBYXAVtcJGxxkUjqEYGkuEgkxFpgNgbm+A7qXAwKoi4kATCbDDCbDJCsJs8HsVuHomhQ3J7gCKbRKoMsIS0pCpFGCVkDPdfpQqCm3oWSyiaUVjWhrMrzs6rOCV14Fjo8aW/ESXvbJU1io8xIjLUgITYCCTERSIi1ID7G8190pIk9kSDBoCAKEi29jwiTAREmAyQJUDUdLlWH0+X2zO8PotBoIUuS54M+NgJD+iV4r3drnhlV5dVNsNc4Ya9xoKLGgYo6JxTVs9+mrlFBXaOCYyVt91cZDRLioi2Ii7agR7TZ8zPKjB7RZsRGmdEjygyLycAw+RkwKIiClBCAUZZhtMiIjjDC7dbh0nS4VT0oexs/ZTTISE2wIjXB2up6IQTqm1RU1jlRWetEZZ0TVXVOVNW7UF3nQlPzrDG3JrxTd8/GbJQRE2VGjNWEmEgzYq0mRFtNiI40IcZqRnSkCVGRJkRFGEN+R3RXYlAQhQAhPIvuWQ0y0HwgtVvXobp1qGrz/g0R3MHRQpIkxEZ5egX9esa22e5U3Kiud6Gm3oWaBgU1DZ6ftY0u1DYoqG9S0LIbR3HrnrDxESYtIswGREWYYI0wIirCCGvz71aL0fsz0mJEg6JDdamItBjYY2nGoCAKUUZZhtEsI7IlODQdmhBwawKaW4eqadB0zzf4UAiQFhFmI3omGtEzMeqM23VdoMGheoatmjz/1Tepzf8paGhSUe9Q0ehQvTPOAHhPC1tZ538tkuQJmEiLEdlX9kb2lX3Ot3khKaBBsX79erz22mtwu9349a9/jZtvvrnV9v3792PhwoVobGzEFVdcgSeffBJGI7OLqCOMBhlGAIeKK/D59pOobVRgi4vALy7rhUt69cBHm4uwq6gSuhDQNYFBfXpgcL9E5O8uQWWdC7FRJowa2hMD+8TjYHE1vtp1GtX1LsTHWDAm8wKcsjdgS2EpXG4NFqMBozNScU1W77PWc6bHANDmuoF94v2+/8A+8dj87UlPHaoGi8lTx4Qr2tYhhMDuo5XYUngaNQ0KoixG9LugB2KsZpy0N+BEWT2cqgaDLMFkNEDTdThdWqulW4QAHC4NDpeGNV8dw3eHKzBxeB9kXJR0Pk9VyAnYAXdlZWWYM2cOPv74Y5jNZtx000344x//iIsvvth7m6lTp+KZZ57BsGHD8Mgjj2Do0KGYO3eu33+DB9yFLrYvMAqLKrDis0MwGGSYjTIUtw5N05EYY8Ghk7WQ5R+W+DbIEsxGCYlxkbCYDJ5v3wLIvDgRO5vXgZIlz/BOfaMLDpcOIQQkCKiagKYLjL/sAoy/vO2H9MHiaqzbegyG5lOrqpoOh9MNSBIiLQbvdZqmY/rV/dqExZnur2k6+iZHY9fRKs9CfxKal1oBfnl5WpvQOttjZA2wYcche5vrp1/dD5f0jkOk1YLTZXVwuNw4fKIG+XtKPftbEiNR26hC03TcPGFASIZF0B1wl5+fjxEjRiAuLg4AkJOTg7y8PPzud78DAJw6dQpOpxPDhg0DAMycORPLli1rV1CEw1GgvrB9oa0r2vf1vjIkJ1hbHZ+guDXYqx1IiotstahJy1es6OaxK1Pzbf+ztwyx0WZYzEbIzYsm6gKIiZJgNkgwGGXIzfevrFeQGBfp+cImmtfIEgInyhtwce84GCTZs34SgKpaJ4QQiI+N8AyFCUDRNBw8UYP0ixJbtWP/99XomRQFo8HgrVZxayitdiAlPrLVv62uCxw8UYuc4X1bPcbeY1Vn/Lc4eKL2jNfvPVaFof0TERtlht68A/6bfWXolRKNSLMJsVEmWCNUKG4NX+8rw7BLbO1+frpaR1+TAQuK8vJy2Gw//EMmJyejsLDwrNttNhvKysra9Tfi4888hhkuOpL8oYTt63yP3D7iZ/+bZ7Lg5ivO6/6P9T//D+E/XNTxx0ht3j9yPo8RTgI2X0zX9VazBcRPFhE713YiIgoOAQuK1NRU2O1272W73Y7k5OSzbq+oqGi1nYiIgkPAgmLUqFEoKChAVVUVHA4HNm7ciLFjx3q3p6WlwWKxYMeOHQCAtWvXttpORETBIaDLjK9fvx6vv/46VFXF7NmzMW/ePMybNw/z589Heno6Dhw4gEcffRQNDQ0YMmQIFi9eDLO5607LSEREbYX0+SiIiCjwuPgJERH5xKAgIiKfGBREROQTg4KIiHwKmaB4+eWXMXnyZEyZMgVvv/02AM8yIdOmTUN2djaWLl3axRWevyVLluChhx4C4FkwcebMmcjJycHChQvhdru7uLqOu+WWWzBlyhTMmDEDM2bMwK5du7B+/XpMnjwZ2dnZWLFiRVeXeF42bdqEmTNnYtKkSXjmmWcAhM9r8+9//7v3eZsxYwaysrLw1FNPhU37AM/U/ClTpmDKlClYsmQJgPB5/y1fvhw5OTmYNm0aXnvtNQAdbJsIAV9//bW46aabhKqqwuFwiPHjx4v9+/eLcePGieLiYqGqqrj99tvF5s2bu7rUDsvPzxfDhw8XDz74oBBCiClTpohvv/1WCCHEww8/LFasWNGV5XWYruti9OjRQlVV73WlpaVi/Pjxorq6WjQ2Nopp06aJw4cPd2GVHVdcXCxGjx4tSkpKhKIoYs6cOWLz5s1h9dpscejQITFhwgRx+vTpsGlfU1OTuPLKK0VlZaVQVVXMnj1bbN26NSzef1u3bhVTp04V9fX1wu12i9/85jdiw4YNHWpbSPQorrrqKrz77rswGo2orKyEpmmoq6tD37590bt3bxiNRkybNg15eXldXWqH1NTUYOnSpbj77rsBnHnBxFBt29GjRwEAt99+O6ZPn47333+/1YKRVqvVu2BkKPrss88wefJkpKamwmQyYenSpYiMjAyb1+aPPfHEE1iwYAFOnDgRNu3TNA26rsPhcMDtdsPtdsNoNIbF+2/fvn0YPXo0oqOjYTAYMGbMGLz33nsdaltIBAUAmEwmLFu2DFOmTMHIkSPPuOhgexcVDBaPP/44FixYgNhYz9m+OmPBxGBRV1eHkSNH4pVXXsE777yDDz74AKdPnw6b5+7777+Hpmm4++67MWPGDKxcuTKsXpst8vPz4XQ6MWnSpLBqX3R0NO677z5MmjQJ48aNQ1paGkwmU1i8/4YMGYItW7agpqYGLpcLmzZtgtFo7FDbQiYoAGD+/PkoKChASUkJjh8/HhaLCv79739Hz549MXLkSO914bRg4mWXXYbnn38eMTExSEhIwOzZs7Fs2bKwaZ+maSgoKMCiRYuwatUqFBYW4sSJE2HTvhYffPABbrvtNgDh9fo8cOAAPvroI/zf//0fvvrqK8iyjK1bt4ZF+0aOHImZM2filltuwZ133omsrCy43e4OtS0kTidXVFQERVFw6aWXIjIyEtnZ2cjLy4PB8MN68j9ddDBU/POf/4TdbseMGTNQW1uLpqYmSJIUNgsmbt++HaqqeoNQCIG0tDSfC0aGkqSkJIwcORIJCQkAgGuvvTZsXpstFEXBtm3b8NxzzwE494KfoWTLli0YOXIkEhM958OYOXMm3nrrrbB4/zU0NCA7O9sb8G+++SZ69eqF7du3e2/jb9tCokdx8uRJPProo1AUBYqi4IsvvsBNN92EY8eOebv+n3zySUguKvj222/jk08+wdq1azF//nxcc801WLx4cdgsmFhfX4/nn38eLpcLDQ0N+Mc//oEXXnjB54KRoWT8+PHYsmUL6urqoGkavvrqK0ycODEsXpstDh48iAsvvBBWq+dkPpmZmWHTvkGDBiE/Px9NTU0QQmDTpk246qqrwuL9d/LkSdxzzz1wu92or6/Hhx9+iNmzZ3eobSHRoxg3bhwKCwtx3XXXwWAwIDs7G1OmTEFCQgLuvfdeuFwujBs3DhMnTuzqUjvNiy++2GrBxNzc3K4uqUPGjx+PXbt24brrroOu65g7dy6ysrKwYMEC5ObmeheMzMjI6OpSOyQzMxN33nkn5s6dC1VVcfXVV2POnDno379/2Lw2T5w4gdTUVO9li8WC5557LizaN3r0aOzbtw8zZ86EyWRCeno67rrrLkyYMCHk33+DBg1CdnY2pk+fDk3TcOuttyIrK6tDny1cFJCIiHwKiaEnIiLqOgwKIiLyiUFBREQ+MSiIiMgnBgUREfkUEtNjifzxzDPPYNu2bQA8B2mmpaUhIiICALBq1Srv78FICIHbbrsNy5Yt8y7lQhQsOD2WwtI111yDl19+Genp6V1dil/cbjeGDBmCbdu2MSgo6LBHQd3C4cOH8eyzz3qPoL711ltx/fXXIz8/H3/6059gs9nw/fffw2q14s4778R7772H48ePY9KkSXjwwQeRn5+PZcuWITk5GceOHYPVasXixYvRv39/KIqC559/Hjt27ICmaRgyZAgWLlyI6OhojB07FllZWThw4ADuv/9+6LqON998E4qioKqqCrNmzcK9996Lhx9+GABw8803480338QNN9yA119/HZdeeikAYOzYsXj99ddhtVpx2223oU+fPigpKcHKlStx7Ngx/O///i+cTidkWcb8+fMxbty4rvznpnDT+augE3W98ePHi8LCQiGEEIqiiEmTJon9+/cLIYSora0VOTk5orCwUGzdulUMHjzYu+3WW28Vc+bMEYqiiIqKCnHppZeKiooKsXXrVjFo0CCxY8cOIYQQ7733nrjhhhuEEEK89NJL4oUXXhC6rgshhFiyZIl4+umnhRBCjBkzRvzlL38RQgihaZqYO3euKC4uFkIIcfr0aTFo0CBRU1MjVFUVAwYMELW1td777du3z9uelsvHjx8XAwYMEDt37hRCCFFVVSWys7PFqVOnhBBClJSUiDFjxoiSkpIA/ctSd8QeBYW9oqIinDhxAg8++KD3OkVRsH//fvTq1Qt9+vTBoEGDAAC9e/dGUlISTCYTEhMTYbVaUVNTA8CzbPPll18OALjhhhvwzDPPoL6+Hps3b0ZTUxO++uorAICqqq0WWsvKygIAyLKM119/HZs3b8batWtx5MgRCCHgdDoRFRXld3tMJhMyMzMBADt37oTdbsdvf/tb73ZZlnHo0KFWy24QnQ8GBYU9XdcRFxeHtWvXeq+z2+2IjY3Fjh07YDabW93eaDzz2+LH1+u6DsDzoaxpGh5//HFcffXVADyrdqqq6r1tSwg0NDTg+uuvR05ODrKysjBr1ix89tlnEGfYTShJUqvrf/x4ERERkGXZW8eAAQPwwQcfeLeXlZV5V7Ml6gycHkth7+KLL4Ysy/j0008BeM4gOHXqVBw4cKBdj7Nnzx4cPnwYgGcW1ZVXXomoqCiMHj0a7733HlRVhaZpeOSRR/DSSy+1uf+xY8fgcDhw3333Yfz48SgoKIDb7YamaTAYDJAkyXv+4oSEBOzZsweA56RBVVVVZ6zpsssuQ1FRkXc10L179yInJweVlZXtahuRL+xRUNgzm8147bXXsGjRIvzlL3+B2+3Gf//3fyMzMxP5+fl+P05ycjJefPFFnDp1CjabDUuWLAEA3HvvvViyZAmuu+46787sBx54oM39Bw8ejNGjR2PSpEkwmUwYNGgQ+vfvj+LiYqSlpSE7Oxtz5szBq6++ivvvvx9PPvkkVqxYgfT0dO9O7Z9KSkrCsmXLsHjxYiiKAiEEXnzxRQ47Uafi9FgiP+Tn52PJkiWthq+IugsOPRERkU/sURARkU/sURARkU8MCiIi8olBQUREPjEoiIjIJwYFERH5xKAgIiKf/j8hUuX1zVFMwgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sns.set(color_codes=True)\n", + "plt.xlim(30,90)\n", + "plt.ylim(0,1)\n", + "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." + ] + } + ], + "metadata": { + "celltoolbar": "Hide code", + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/module4/exo/challenger.pdf b/module4/exo/challenger.pdf new file mode 100644 index 0000000000000000000000000000000000000000..ee646b2c43bc53934bb641dda59f032786bc0a11 Binary files /dev/null and b/module4/exo/challenger.pdf differ