diff --git a/module4/challenger.ipynb b/module4/challenger.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..e2e81fc6bee46203d8fc583e79e507c14c57a1d4 --- /dev/null +++ b/module4/challenger.ipynb @@ -0,0 +1,802 @@ +{ + "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.3 (v3.7.3:ef4ec6ed12, Mar 25 2019, 22:22:05) [MSC v.1916 64 bit (AMD64)]\n", + "uname_result(system='Windows', node='Travelmate-P253', release='7', version='6.1.7601', machine='AMD64', processor='Intel64 Family 6 Model 58 Stepping 9, GenuineIntel')\n", + "IPython 7.18.1\n", + "IPython.core.release 7.18.1\n", + "html 6.1.4\n", + "_csv 1.0\n", + "_ctypes 1.1.0\n", + "decimal 1.70\n", + "argparse 1.1\n", + "backcall 0.2.0\n", + "colorama 0.4.3\n", + "csv 1.0\n", + "ctypes 1.1.0\n", + "cycler 0.10.0\n", + "dateutil 2.8.0\n", + "decimal 1.70\n", + "decorator 4.4.2\n", + "distutils 3.7.3\n", + "ipykernel 5.3.4\n", + "ipykernel._version 5.3.4\n", + "ipython_genutils 0.2.0\n", + "ipython_genutils._version 0.2.0\n", + "jedi 0.17.2\n", + "json 2.0.9\n", + "jupyter_client 6.1.7\n", + "jupyter_client._version 6.1.7\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.1\n", + "matplotlib.backends.backend_agg 3.1.1\n", + "numpy 1.16.4\n", + "numpy.core 1.16.4\n", + "numpy.core._multiarray_umath 3.1\n", + "numpy.lib 1.16.4\n", + "numpy.linalg._umath_linalg b'0.1.5'\n", + "pandas 0.25.0\n", + "_libjson 1.33\n", + "parso 0.7.1\n", + "patsy 0.5.1\n", + "patsy.version 0.5.1\n", + "pickleshare 0.7.5\n", + "platform 1.0.8\n", + "prompt_toolkit 3.0.7\n", + "pygments 2.7.1\n", + "pyparsing 2.4.1\n", + "pytz 2019.1\n", + "re 2.2.1\n", + "scipy 1.3.0\n", + "scipy._lib.decorator 4.0.5\n", + "scipy._lib.six 1.2.0\n", + "scipy.fftpack._fftpack b'$Revision: $'\n", + "scipy.fftpack.convolve b'$Revision: $'\n", + "scipy.integrate._dop b'$Revision: $'\n", + "scipy.integrate._ode $Id$\n", + "scipy.integrate._odepack 1.9 \n", + "scipy.integrate._quadpack 1.13 \n", + "scipy.integrate.lsoda b'$Revision: $'\n", + "scipy.integrate.vode b'$Revision: $'\n", + "scipy.interpolate._fitpack 1.7 \n", + "scipy.interpolate.dfitpack b'$Revision: $'\n", + "scipy.linalg 0.4.9\n", + "scipy.linalg._fblas b'$Revision: $'\n", + "scipy.linalg._flapack b'$Revision: $'\n", + "scipy.linalg._flinalg b'$Revision: $'\n", + "scipy.ndimage 2.0\n", + "scipy.optimize._cobyla b'$Revision: $'\n", + "scipy.optimize._lbfgsb b'$Revision: $'\n", + "scipy.optimize._minpack 1.10 \n", + "scipy.optimize._nnls b'$Revision: $'\n", + "scipy.optimize._slsqp b'$Revision: $'\n", + "scipy.optimize.minpack2 b'$Revision: $'\n", + "scipy.signal.spline 0.2\n", + "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", + "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", + "scipy.special.specfun b'$Revision: $'\n", + "scipy.stats.mvn b'$Revision: $'\n", + "scipy.stats.statlib b'$Revision: $'\n", + "seaborn 0.9.0\n", + "seaborn.external.husl 2.1.0\n", + "seaborn.external.six 1.10.0\n", + "six 1.12.0\n", + "statsmodels 0.12.0\n", + "statsmodels.__init__ 0.12.0\n", + "statsmodels.api 0.12.0\n", + "statsmodels.tools.web 0.12.0\n", + "traitlets 5.0.4\n", + "traitlets._version 5.0.4\n", + "urllib.request 3.7\n", + "wcwidth 0.2.5\n", + "zlib 1.0\n", + "zmq 19.0.2\n", + "zmq.sugar 19.0.2\n", + "zmq.sugar.version 19.0.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": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
56/18/82672500
68/30/836731000
711/28/836701000
82/03/846572001
94/06/846632001
108/30/846702001
1110/05/846782000
1211/08/846672000
131/24/856532002
144/12/856672000
154/29/856752000
166/17/856702000
177/2903/856812000
188/27/856762000
1910/03/856792000
2010/30/856752002
2111/26/856762000
221/12/866582001
\n", + "
" + ], + "text/plain": [ + " Date Count Temperature Pressure Malfunction\n", + "0 4/12/81 6 66 50 0\n", + "1 11/12/81 6 70 50 1\n", + "2 3/22/82 6 69 50 0\n", + "3 11/11/82 6 68 50 0\n", + "4 4/04/83 6 67 50 0\n", + "5 6/18/82 6 72 50 0\n", + "6 8/30/83 6 73 100 0\n", + "7 11/28/83 6 70 100 0\n", + "8 2/03/84 6 57 200 1\n", + "9 4/06/84 6 63 200 1\n", + "10 8/30/84 6 70 200 1\n", + "11 10/05/84 6 78 200 0\n", + "12 11/08/84 6 67 200 0\n", + "13 1/24/85 6 53 200 2\n", + "14 4/12/85 6 67 200 0\n", + "15 4/29/85 6 75 200 0\n", + "16 6/17/85 6 70 200 0\n", + "17 7/2903/85 6 81 200 0\n", + "18 8/27/85 6 76 200 0\n", + "19 10/03/85 6 79 200 0\n", + "20 10/30/85 6 75 200 2\n", + "21 11/26/85 6 76 200 0\n", + "22 1/12/86 6 58 200 1" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")\n", + "data = pd.read_csv(\"C:/Users/José/Desktop/Cours/Moocs Plateforme FUN/Recherche reproductible/Semaine 4/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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYMklEQVR4nO3df5TU9X3v8edrAWERIgRSagWiVuqNRy3KBjX2B0STg54TaK6a4D3RNC2h90ROTkyaaHtzreWm59zYJLa5sYnEaxrtSYhKo9xeev0RJak9/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": [ + "c:\\program files\\python37\\lib\\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: Fri, 09 Oct 2020 Deviance: 3.0144
Time: 17:52:58 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: Fri, 09 Oct 2020 Deviance: 3.0144\n", + "Time: 17:52:58 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": [ + "c:\\program files\\python37\\lib\\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: Fri, 09 Oct 2020 Deviance: 18.086
Time: 17:52:58 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: Fri, 09 Oct 2020 Deviance: 18.086\n", + "Time: 17:52:58 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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZd0lEQVR4nO3df3TV9Z3n8efLACX8EAekWyFY6AyNdf3BDwkqqxOtCvZUtLsisk4du1K6Z8Y6blf2yFmnWqvn7E7cHTsdx5GpjqO2CLqKtIdpUMdsZzyKYEEQ2AhaKoF2UBQEGzSJ7/3jexNDTMhNcvPjfng9zsnhfr/3c7/f9+d+ua9887nf+7mKCMzMLF3H9XcBZmbWuxz0ZmaJc9CbmSXOQW9mljgHvZlZ4hz0ZmaJ6zToJT0oaa+k1zq4X5L+StIOSZskTSt8mWZm1l35nNE/BMw5yv2XApNzP4uA+3pelpmZFUqnQR8RvwDePUqTy4GHI/MScIKkkwpVoJmZ9cygAmxjPLCr1XJdbt1v2jaUtIjsrJ+SoSOmDx/zuQLsfmAKQP1dRC9KuX8p9w3cv2L3/u4d70TE2K48phBB395z2u68ChGxFFgKUF5eHrW1tQXY/cBUU1NDZWVlf5fRa1LuX8p9A/ev2En6dVcfU4irbuqACa2Wy4A9BdiumZkVQCGCfhVwbe7qm7OBAxHxqWEbMzPrH50O3UhaBlQCJ0qqA24DBgNExN8Cq4GvADuA3wHf6K1izcys6zoN+ohY0Mn9AfxpwSoys6LQ0NBAXV0dhw8f7u9SjjBq1Ci2bdvW32X02NChQykrK2Pw4ME93lYh3ow1s2NQXV0dI0eOZOLEiUgD5zqXgwcPMnLkyP4uo0cign379lFXV8ekSZN6vD1PgWBm3XL48GHGjBkzoEI+FZIYM2ZMwf5actCbWbc55HtPIZ9bB72ZWeI8Rm9mRaukpITTTz+9ZXnlypWMGTOmHysamBz0Zla0SktL2bhx4xHrDh482HK7sbGRQYMccx66MbOk/PjHP2bevHlcdtllXHLJJQBUVVUxY8YMzjjjDG677baWtnfddRfl5eVcdNFFLFiwgLvvvhuAyspK1q9fD8A777zDxIkTAWhqamLx4sUt27r//vuBT6ZduPLKKznllFO45ppryK48h3Xr1nHuuedy5plnUlFRwcGDBznvvPOO+AU1a9YsNm3a1GvPiX/VmVmPfe+nW9i65/2CbvPUccdz22X/9qht6uvrmTJlCgCTJk3iqaeeAuDFF19k06ZNjB49mjVr1rB9+3ZefvllIoK5c+fyi1/8guHDh/PYY4+xYcMGGhsbmTZtGtOnTz/q/h544AFGjRrFunXr+PDDD5k1a1bLL5MNGzawZcsWxo0bx6xZs3jhhReoqKhg/vz5LF++nBkzZvD+++9TWlrKwoULeeihh7jnnnt4/fXX+fDDDznjjDMK8Ky1z0FvZkWrvaEbgIsvvpjRo0cDsGbNGtasWcPUqVMBOHToENu3b+fgwYN87WtfY9iwYQDMnTu30/2tWbOGTZs28cQTTwBw4MABtm/fzpAhQ6ioqKCsrAyAKVOmsHPnTkaNGsVJJ53EjBkzADj++OMBmDdvHt///vepqqriwQcf5LrrruvZE9EJB72Z9VhnZ959bfjw4S23I4IlS5bwrW9964g299xzT4eXMA4aNIiPP/4Y4Ihr2SOCH/7wh8yePfuI9jU1NXzmM59pWS4pKaGxsZGIaHcfw4YN4+KLL+bpp59mxYoVLcNEvcVj9GaWtNmzZ/Pggw9y6NAhAHbv3s3evXs5//zzeeqpp6ivr+fgwYP89Kc/bXnMxIkTeeWVVwBazt6bt3XffffR0NAAwOuvv84HH3zQ4b5POeUU9uzZw7p164DsjeLGxkYAFi5cyI033siMGTNa/vroLT6jN7OkXXLJJWzbto1zzjkHgBEjRvDoo48ybdo05s+fz5QpU/j85z/Peeed1/KYm2++mauuuopHHnmECy+8sGX9woUL2blzJ9OmTSMiGDt2LCtXruxw30OGDGH58uV8+9vfpr6+ntLSUp599llGjBjB9OnTOf744/nGN/pgHsiI6JefL37xi5Gy559/vr9L6FUp9y/lvkUUrn9bt24tyHYK7f333+/W42677baoqqoqcDUd2717d0yePDmampo6bNPecwysjy7mrYduzMz62MMPP8zMmTO56667OO643o9hD92YmQG33357n+3r2muv5dprr+2z/fmM3sy6LaLdr4e2Aijkc+ugN7NuGTp0KPv27XPY94LIzUc/dOjQgmzPQzdm1i1lZWXU1dXx9ttv93cpRzh8+HDBArI/NX/DVCE46M2sWwYPHlyQbz8qtJqampZPwVrGQzdmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVni8gp6SXMk1UraIemWdu4/WdLzkjZI2iTpK4Uv1czMuqPToJdUAtwLXAqcCiyQdGqbZrcCKyJiKnA18DeFLtTMzLonnzP6CmBHRLwZER8BjwGXt2kTwPG526OAPYUr0czMekKdfYO7pCuBORGxMLf8dWBmRNzQqs1JwBrg94DhwEUR8Uo721oELAIYO3bs9BUrVhSqHwPOoUOHGDFiRH+X0WtS7l/KfQP3r9hdcMEFr0TEWV15TD5fDq521rX97bAAeCgi/pekc4BHJJ0WER8f8aCIpcBSgPLy8qisrOxKrUWlpqYG9684pdw3cP+ORfkM3dQBE1otl/HpoZnrgRUAEfEiMBQ4sRAFmplZz+QT9OuAyZImSRpC9mbrqjZt3gK+DCDpS2RB/3YhCzUzs+7pNOgjohG4AagGtpFdXbNF0h2S5uaa/Vfgm5JeBZYB10Vng/9mZtYn8hmjJyJWA6vbrPtuq9tbgVmFLc3MzArBn4w1M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHEOejOzxDnozcwS56A3M0ucg97MLHF5Bb2kOZJqJe2QdEsHba6StFXSFkk/KWyZZmbWXYM6ayCpBLgXuBioA9ZJWhURW1u1mQwsAWZFxHuSPttbBZuZWdfkc0ZfAeyIiDcj4iPgMeDyNm2+CdwbEe8BRMTewpZpZmbd1ekZPTAe2NVquQ6Y2abNFwEkvQCUALdHxM/bbkjSImARwNixY6mpqelGycXh0KFD7l+RSrlv4P4di/IJerWzLtrZzmSgEigD/lnSaRGx/4gHRSwFlgKUl5dHZWVlV+stGjU1Nbh/xSnlvoH7dyzKZ+imDpjQarkM2NNOm6cjoiEifgXUkgW/mZn1s3yCfh0wWdIkSUOAq4FVbdqsBC4AkHQi2VDOm4Us1MzMuqfToI+IRuAGoBrYBqyIiC2S7pA0N9esGtgnaSvwPLA4Ivb1VtFmZpa/fMboiYjVwOo2677b6nYA38n9mJnZAOJPxpqZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVni8gp6SXMk1UraIemWo7S7UlJIOqtwJZqZWU90GvSSSoB7gUuBU4EFkk5tp91I4EZgbaGLNDOz7svnjL4C2BERb0bER8BjwOXttPs+8BfA4QLWZ2ZmPTQojzbjgV2tluuAma0bSJoKTIiIn0m6uaMNSVoELAIYO3YsNTU1XS64WBw6dMj9K1Ip9w3cv2NRPkGvdtZFy53SccBfAtd1tqGIWAosBSgvL4/Kysq8iixGNTU1uH/FKeW+gft3LMpn6KYOmNBquQzY02p5JHAaUCNpJ3A2sMpvyJqZDQz5BP06YLKkSZKGAFcDq5rvjIgDEXFiREyMiInAS8DciFjfKxWbmVmXdBr0EdEI3ABUA9uAFRGxRdIdkub2doFmZtYz+YzRExGrgdVt1n23g7aVPS/LzMwKxZ+MNTNLnIPezCxxDnozs8Q56M3MEuegNzNLXF5X3ZgVysoNu6mqrmXP/nrGnVDK4tnlXDF1fH+XZb3Ax3rgcNBbn1m5YTdLntxMfUMTALv317Pkyc0ADoDE+FgPLB66sT5TVV3b8sJvVt/QRFV1bT9VZL3Fx3pgcdBbn9mzv75L6614+VgPLA566zPjTijt0norXj7WA4uD3vrM4tnllA4uOWJd6eASFs8u76eKrLf4WA8sfjPW+kzzm3C+EiN9PtYDi4Pe+tQVU8f7xX6M8LEeODx0Y2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWuLyCXtIcSbWSdki6pZ37vyNpq6RNkp6T9PnCl2pmZt3RadBLKgHuBS4FTgUWSDq1TbMNwFkRcQbwBPAXhS7UzMy6J58z+gpgR0S8GREfAY8Bl7duEBHPR8TvcosvAWWFLdPMzLprUB5txgO7Wi3XATOP0v564B/bu0PSImARwNixY6mpqcmvyiJ06NAh969Ipdw3cP+ORfkEvdpZF+02lP4IOAv4w/buj4ilwFKA8vLyqKyszK/KIlRTU4P7V5xS7hu4f8eifIK+DpjQarkM2NO2kaSLgP8O/GFEfFiY8szMrKfyGaNfB0yWNEnSEOBqYFXrBpKmAvcDcyNib+HLNDOz7uo06COiEbgBqAa2ASsiYoukOyTNzTWrAkYAj0vaKGlVB5szM7M+ls/QDRGxGljdZt13W92+qMB1mXXZyg27qaquZc/+esadUMri2eUAn1p3xdTxfVpDb+6vK25duZlla3dx02kNXL9kNQtmTuDOK07v77KsD+QV9GYD3coNu1ny5GbqG5oA2L2/nsWPvwqChqZoWbfkyc0AvRK+7dXQm/vriltXbubRl95qWW6KaFl22KfPUyBYEqqqa1sCtlnDx9ES8s3qG5qoqq7tsxp6c39dsWztri6tt7Q46C0Je/bX90rbQtTQW/vriqZo94roDtdbWhz0loRxJ5T2SttC1NBb++uKErX3cZiO11taHPSWhMWzyykdXHLEusHHicElRwZZ6eCSljdp+6KG3txfVyyYOaFL6y0tfjPWktD8Zmd/XnXTUQ39/UYsfPKGa/OYfInkq26OIQ56S8YVU8e3G6p9GbQd1TAQ3HnF6dx5xenU1NTwxjWV/V2O9SEP3ZiZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZolz0JuZJc5Bb2aWOAe9mVniHPRmZokblE8jSXOAHwAlwI8i4n+0uf8zwMPAdGAfMD8idha2VLNjw8oNu6mqrmXP/nrGnVDK4tnlXDF1PNf83Yu88Ma7Le1m/f5o5p11crtt29vG+l+/y7K1u7jptAauX7KaBTMncOcVp3epho7Wd2Ubt67czLK1u2iKoETqch1A3jUcrY5jSadBL6kEuBe4GKgD1klaFRFbWzW7HngvIv5A0tXA/wTm90bBZilbuWE3S57cTH1DEwC799ez5MnN3Pv8drbv/eCIti+88e4Rwd/cdv2v3+X/vLL7iG18Z/lGPm712KYIHn3pLYBPhWxHNbS33SVPbgb4VHB2tI3H1791RM1drWPx46+CoKEpOq3haHV01D5V+QzdVAA7IuLNiPgIeAy4vE2by4F/yN1+AviyJBWuTLNjQ1V1bUsoNatvaPpUyHekvqGJZWt3fWobH3fQftnaXXnX0N526xuaqKquzXsbrUO+O3U0fBwtId9ZDUero6P2qVJEHL2BdCUwJyIW5pa/DsyMiBtatXkt16Yut/xGrs07bba1CFiUWzwNeK1QHRmATgTe6bRV8Uq5f/3WtyGf+4Ppvb2Ppt8doGTYqJblj36745We1jAQttHq8S3H72jbaLu/IlIeESO78oB8xujbOzNv+9shnzZExFJgKYCk9RFxVh77L0ruX/FKuW+Q9a/xwN6k+5f68evqY/IZuqkDJrRaLgP2dNRG0iBgFND+32hmZtan8gn6dcBkSZMkDQGuBla1abMK+OPc7SuBf4rOxoTMzKxPdDp0ExGNkm4Aqskur3wwIrZIugNYHxGrgAeARyTtIDuTvzqPfS/tQd3FwP0rXin3Ddy/Ytfl/nX6ZqyZmRU3fzLWzCxxDnozs8T1SdBLGirpZUmvStoi6Xu59ZMkrZW0XdLy3Ju9RUlSiaQNkn6WW06pbzslbZa0sfnSLkmjJT2T698zkn6vv+vsLkknSHpC0v+TtE3SOan0T1J57rg1/7wv6aaE+vdfcpnymqRluaxJ6bX3Z7m+bZF0U25dl49dX53RfwhcGBFnAlOAOZLOJpsq4S8jYjLwHtlUCsXqz4BtrZZT6hvABRExpdX1ybcAz+X691xuuVj9APh5RJwCnEl2HJPoX0TU5o7bFLK5qH4HPEUC/ZM0HrgROCsiTiO7WKR5Cpaif+1JOg34JtnsBGcCX5U0me4cu4jo0x9gGPBLYCbZp9cG5dafA1T3dT0F6lNZ7gm/EPgZ2QfIkuhbrv6dwIlt1tUCJ+VunwTU9ned3ezb8cCvyF2YkFr/2vTpEuCFVPoHjAd2AaPJriD8GTA7ldceMI9sEsnm5T8H/lt3jl2fjdHnhjY2AnuBZ4A3gP0R0ZhrUkd24IrRPWQHoHlKkTGk0zfIPuW8RtIruWksAP5NRPwGIPfvZ/utup75AvA28Pe5obcfSRpOOv1r7WpgWe520fcvInYDdwNvAb8BDgCvkM5r7zXgfEljJA0DvkL2wdQuH7s+C/qIaIrsz8cysj9FvtRes76qp1AkfRXYGxGt583Ia0qIIjIrIqYBlwJ/Kun8/i6ogAYB04D7ImIq8AFFOIzRmdw49Vzg8f6upVByY9OXA5OAccBwsv+jbRXlay8itpENQz0D/Bx4FWg86oM60OdX3UTEfqAGOBs4ITdlArQ/tUIxmAXMlbSTbGbPC8nO8FPoGwARsSf3716y8d0K4F8lnQSQ+3dv/1XYI3VAXUSszS0/QRb8qfSv2aXALyPiX3PLKfTvIuBXEfF2RDQATwLnktZr74GImBYR55N9GHU73Th2fXXVzVhJJ+Rul5IdoG3A82RTJkA2hcLTfVFPIUXEkogoi4iJZH8a/1NEXEMCfQOQNFzSyObbZOO8r3HktBdF27+I+C2wS1J5btWXga0k0r9WFvDJsA2k0b+3gLMlDctNi9587JJ47QFI+mzu35OBf092DLt87Prkk7GSziCbr76E7JfLioi4Q9IXyM6CRwMbgD+KiA97vaBeIqkSuDkivppK33L9eCq3OAj4SUTcJWkMsAI4mewFNy8iinIiO0lTgB8BQ4A3gW+Q+39KGv0bRvam5Rci4kBuXRLHL3ep9nyyIY0NwEKyMfmif+0BSPpnsvf8GoDvRMRz3Tl2ngLBzCxx/mSsmVniHPRmZolz0JuZJc5Bb2aWOAe9mVni8vlycLM+lbt87Lnc4ueAJrJpCgAqIuKjfinsKCT9J2B17rp8swHFl1fagCbpduBQRNw9AGopiYimDu77F+CGiNjYhe0NajUni1mv8dCNFRVJf6zsuw02SvobScdJGiRpv6QqSb+UVC1ppqT/K+lNSV/JPXahpKdy99dKujXP7d4p6WWgQtL3JK3LzRH+t8rMJ5t+e3nu8UMk1bX6NPjZkp7N3b5T0v2SniGbSG2QpP+d2/cmSQv7/lm11DnorWjk5uf+GnBuboK8QXzyRfSjgDW5ydc+Am4n+0j8POCOVpupyD1mGvAfJU3JY7u/jIiKiHgR+EFEzABOz903JyKWAxuB+ZHN/d7Z0NJU4LKI+DqwiGxSvApgBtmkcSd35/kx64jH6K2YXEQWhuuzqU0oJftoP0B9RDyTu70ZOBARjZI2AxNbbaM6It4DkLQS+Hdkr4OOtvsRn0wBAfBlSYuBocCJZNPi/mMX+/F0RBzO3b4E+JKk1r9YJpN9tN2sIBz0VkwEPBgRf37EymymwtZn0R+TfatZ8+3W/8/bvikVnWy3PnJvZOXmjPlrYFpE7JZ0J1ngt6eRT/5ibtvmgzZ9+pOIeA6zXuKhGysmzwJXSToRsqtzujHMcYmy74gdRjaX+Qtd2G4p2S+Od3Izev6HVvcdBEa2Wt5J9tV9tGnXVjXwJ83T6ir7jtfSLvbJ7Kh8Rm9FIyI252YrfFbScWQz+v1nujbf+L8APwF+H3ik+SqZfLYbEfsk/QPZNM2/Bta2uvvvgR9Jqid7H+B24O8k/RZ4+Sj13E82C+HG3LDRXrJfQGYF48sr7ZiRu6LltIi4qb9rMetLHroxM0ucz+jNzBLnM3ozs8Q56M3MEuegNzNLnIPezCxxDnozs8T9f6DkuZ+X6fjJAAAAAElFTkSuQmCC\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU1f34/9e9d5ZkkpA9BNlRSdhVXABZKlXCjiJtXb5Sl9Laflpaf59P66dqd1tr6+djaz/WSm21KrRYW2VRAy6lyqIiKnsIq2wh+57Z7r3n98ckAwjEJGQymcn7+XiEzHZn3odk5p1zzznvoymlFEIIIcQ56NEOQAghRPcmiUIIIUSrJFEIIYRolSQKIYQQrZJEIYQQolWSKIQQQrQq4omioaGB2bNnc/To0TPu2717N/Pnz6egoID7778f0zQjHY4QQoh2imii2Lp1KzfffDOHDh066/3f/e53+eEPf8iaNWtQSvHCCy9EMhwhhBAdENFE8cILL/CjH/2InJycM+47duwYPp+PSy65BID58+dTWFgYyXCEEEJ0gCOST/7zn//8nPeVlZWRnZ0dvp6dnU1paWkkwxFCCNEBURvMtm0bTdPC15VSp10XQgjRPUS0R9Ga3NxcysvLw9crKirOeoqqNZ8crca0bNA0dEDTNDQNNDQ0PXQdLZQNNV1DU4Suh24GTUPTNPRT8lN3qXyVmZlMZWVDtMOIGGlf7IrntkF8t0/XNdLTk9p9XNQSRd++fXG73WzZsoWxY8eyYsUKJk+e3K7n8AVMAkG7Xcdo4X8I92B0LXTZ0DR0XUM3NPSWy82JpOWypnVdMrHtbpK1IkTaF7viuW0Q/+1rry5PFIsWLWLx4sWMGjWKRx55hAceeICGhgZGjBjBwoULI/76KvxP6HQXgN18T/Acx2jNiUQDdEPDqesYDh1D10776i69ESGE6ExaLJcZ33OgvN09is7WkkR0DRyGjtOh42hJIpqOrnesB5KdnUJ5eX3nB9xNSPtiVzy3DeK7fbqukZmZ3O7jonbqKV4oFeqZ2IBpWfgCFnAygTgNHaezOYHoOg5Deh5CiNgiiSJCWhKI37bwB0PJQ9c0dAPcTgduRyiBaMhMLyFE9yaJogvZSmGbYJpBGgl1A91OgwS3gcshSUMI0T1Joogi21Z4/SZev3kyabgM3E4j2qEJIUSYJIpu4tSk4dA1HG4npqVwOmRMQwgRXZIouiHTVjT5TarrvDgdBkkJDlxOQ05MCSGiQhJFN6YUBIIWgaCFw9DwJDhJdBsyliGE6FKycVGMMC1FXWOAylof3oApp6OEEF1GEkWMMS1FbUOAynovgaCN1FEUQkSaJIoYZZqK6gYftY0BFNK9EEJEjiSKGKYUNPlMKut8BMzoljIRQsQvSRRxwDQV1fU+Gv0mSO9CCNHJJFHECaWgvjFAbWNQTkUJITqVJIo44/WbVNf5sW05FSWE6BySKOJQwLSpqvdjWtKzEEKcP0kUccq0QuMWMsgthDhfkijimGUraur9BJrLnAshREdIoohztlLUNASivhOgECJ2SaLoAULJwi+noYQQHSKJooeQZCGE6ChJFD2IbSvqGvzYSpKFEKLtJFH0MKatqK4PYEv5WSFEG0mi6IGCpk1dY0DWbwsh2kQSRQ/lC1g0eINSplwI8ZkkUfRgTd4g3oCssRBCtE4SRQ+mgLrGAKYtJ6GEEOcmiaKHa5kJJRVnhRDnIolCEDBt6ptkvEIIcXaSKAQAXp9Jk9+MdhhCiG5IEoUAQuMV9Y1BgpYsxhNCnC6mE8XLbx+grjEQ7TDihq0UtQ1+LBncFkKcIqYTRfHRWh59YSvv7SqVlcadxLRUaDGe/HcKIZrFdKLwuB34gxYr1h/kqdW7qKzzRTukuOAPWjT4ZHBbCBES04nizlnDuPTiLAAOldTz2Ivb2LjjhPQuOkGTTxbjCSFCYjpReBIcfOGai/jy9Dx6JbkImjarNx7iT6t3UV0vvYvzoRTUNwak0qwQIrKJYtWqVcycOZNp06axdOnSM+7fuXMnN954I3PnzuVrX/sadXV1HXqdvAHpfHvBaMYOzQbgYEk9j724nQ+Ly1HSu+gwy1bUNgRlKZ4QPVzEEkVpaSmPPvooy5Yt4+WXX2b58uXs27fvtMf8/Oc/Z/HixaxcuZLBgwfzpz/9qcOvl+h2cOPnLmTh9DySE534gxYvrtvPsjf20uQLnm9zeix/0KLRK/9/QvRkEUsUGzduZNy4caSlpeHxeCgoKKCwsPC0x9i2TWNjIwBer5eEhITzft38AeksXjCa4YPSAdh5sIrf/WM7B453rLcioNEXlJ3xhOjBNBWhczNPPvkkTU1N3HPPPQD8/e9/Z9u2bfzsZz8LP+bjjz/mzjvvxOPxkJiYyAsvvEB6enqbX+NEZeM55/wrpdi0vYTlrxfjD1poGsycMJiZVw/C0GN6aCYqDB2y0zwYhvzfCdHTOCL1xLZto50yv1Ipddp1n8/H/fffzzPPPMPo0aN5+umnuffee1myZEmbX6O2tolA8Nx/6eb3S+U/bhjJ397cy/HKJl7ZcJAd+yu46fMXk5rk6ljDukhGRhJVVY3RDuM0jQ1+UpNcnbLGIjs7hfLy+vN/om4qntsXz22D+G6frmtkZia3/7gIxAJAbm4u5eXl4evl5eXk5OSErxcXF+N2uxk9ejQAX/rSl3j//fc7PY6stETuvn4kE0f1AeCTE/X83z+2sfdoTae/Vrzz+k38ATkFJURPE7FEMWHCBDZt2kRVVRVer5e1a9cyefLk8P0DBw7kxIkTHDhwAIA333yTUaNGRSQWh6Ezc/xAFhbkkeg2aPSZPPNqEW98cARbylW0mVJQ1+SXdSpC9DARO/XUu3dv7rnnHhYuXEgwGGTBggWMHj2aRYsWsXjxYkaNGsVDDz3Ed77zHZRSZGZm8otf/CJS4QCQPzCdb84fzV/fKOZoeSNvfXiMI2UNfGnqxXgSIvZfEVdMS9HgDdLL071P3QkhOk/EBrO7wp4D5a2OUZyLadm89u5hNu08AUB6iptbrxvKBVlJnR1ih3XHMYoWmgbpKQm4HB3vkMbzeWCI7/bFc9sgvtvX7cYoujOHoTPn6kF84ZoLcRo61fV+/rBiBx/vrYh2aDFBKWjwykI8IXqKHpkoWlx6cTZ3Xz+C9BQ3pqV44V/7KHzvExm3aINA0JKNjoToIXp0ogDok5nEf9wwiov6pgLw9tYSnluzB19APgQ/S6M3iG3LLCgh4l2PTxQQKi745Rn5XD0yF4A9R2p44uUdUrb8M9i2os4r5ciFiHeSKJoZusasCYO4ccoQDF2jvMbHEy/t4NAJKf3RGn/AkrUVQsQ5SRSfMjYvh7tmD8OT4KDJb/Kn1bv5aG/5Zx/YQykF9d6ADGwLEcckUZzFoNxefOP6kWSnJWDZir//az9vbjkqJcvPIWjaeGVMR4i4JYniHDJ6JXD3vJHhQe43txzln/8+gCWDt2fV2BSUFdtCxClJFK1IdDv48ow8LmveEGlLcTl/eU1mRJ2NZSsaffL/IkQ8kkTxGQxd58YpQ5h6WV8A9h2r5Y+rdlHfFIhyZN2P12dKj0uIOCSJog00TePay/tz45Qh6BqUVDbxhxU7qayV6bOnspXC14GSKkKI7k0SRTuMzcvh/xXknVb241h5Q7TD6la8ftk2VYh4I4minfIHpHPX7GEkuh00+kz+uHoX+4/XRjusbsOylGybKkSckUTRAQN6p/C1uSNITXIRCNr85bUidh+qinZY3YJS4A9a0Q5DCNGJJFF0UE56Il+bN4Ks1ARMS7H09WJZmNfML8UChYgrkijOQ1qym6/OHcEFmR5sBX//137ebd7joiczbUVAehVCxA1JFOcpOdHJV+YMZ1BuCgArNxzina3HoxxV9Hn9kiiEiBeSKDpBgsvB7TPzw6u4X3vvcI8v+eE3LVmpLUSckETRSVwOg9sK8sgfkA6ESn6sef9Ij00Wtq1kUFuIOCGJohM5HTq3TruYUUMyAHh763Fee+9wj00WjT7ZLlWIeCCJopMZus6Xpl7MpRdnAbB+WwmvbPqkRyYL05RehRDxQBJFBOi6xo1TLgwXE9y44wSrNhzqkcmi0SsrtYWIdZIoIkTXNeZPGcLl+TkAvLurtEcmC9O0ZaqsEDFOEkUE6ZrG9ZMGc+Wwk8li9caedRpKAQ0+s/mSECIWSaKIMF3TmDtxMFc09yw27TzR48YsgqZFwOw57RUi3kii6AK6pjFv0mAuzzs5ZvHauz1nNpRS0OQNomnRjkQI0RGSKLqIrmlcP3kIY5uTxfrtJazd3HPWWfhNi6AlVWWFiEWSKLqQrmncMGlIeOrsvz8+zlsfHotyVF1DKfAFZFBbiFgkiaKLhWZDXRhelPfmlqP8++OekSx8sqe2EDFJEkUUGLrGF6dexPBBoXIfa94/wsYd8V911pSyHkLEJEkUUWLoOjd9/mKG9k8DYPXGQ3xQVBblqCLP6zNlUFuIGCOJIoochs6t1w1lcJ9eALz09gG27quIclSRFTBtTBnUFiKmSKKIMqdDZ2FBHv1zklGENj/a/Ul1tMOKGFspGdQWIsa0KVE899xzNDQ0RDqWHsvtMrh9Rj59Mj3YSvHXN4rZE8fJwus3se2eMS1YiHjQpkSxZ88eCgoKuP/++9m+fXubn3zVqlXMnDmTadOmsXTp0jPuP3DgALfddhtz587lrrvuora2tu2Rx5lEt4M7Zg4L78H9+39s5Wh5fCZn01J4AzIDSohY0aZE8eCDD7JmzRpGjhzJT37yE2688UZefPFF/H7/OY8pLS3l0UcfZdmyZbz88sssX76cffv2he9XSvH1r3+dRYsWsXLlSoYNG8aSJUvOv0UxLDnRyZ2zhpGa5MIfsHjm1SJKq5uiHVZENHllUFuIWNHmMYrk5GSmT5/O7NmzqampYdmyZUyfPp233nrrrI/fuHEj48aNIy0tDY/HQ0FBAYWFheH7d+7cicfjYfLkyQDcfffd3HrrrefZnNiXluzmzlnDSPE4afKbPP3KbqrrfdEOq9MFZaW2EDGjTYli06ZNfOc732H69OkcOHCAxx9/nH/+85/85S9/4Yc//OFZjykrKyM7Ozt8PScnh9LS0vD1w4cPk5WVxX333ccNN9zAj370Izwez3k2Jz5kpyXyrS9eittpUNcU5M+vFtEQZ/s6KMDrl0FtIWKBoy0P+slPfsItt9zCz372M1JSUsK3DxgwgC9+8YtnPca2bbRTzi0opU67bpom77//Ps8//zyjRo3iN7/5Db/85S/55S9/2ebgU1M9WHE6KJoBfPMLY3jshY+prPXx3Npi/r9bLiPR3aYfWUxI8LhIT/fgMOJz8l12dspnPyhGxXPbIP7b115t+tRZuXIlhYWFpKSkUF5eziuvvMLChQvRdZ3Fixef9Zjc3Fw++OCD8PXy8nJycnLC17Ozsxk4cCCjRo0CYPbs2ed8rnOprW0iEIzP0xcZGUlkJru4aepFLH29mCOl9Tz2tw+5fcYwnI7Y/2DNyEiisrIR0x8k0RU/ya9FdnYK5eX10Q4jIuK5bRDf7dN1jczM5PYf15YH/exnP2PdunXNL6SzZcsWfvGLX7R6zIQJE9i0aRNVVVV4vV7Wrl0bHo8AuPTSS6mqqqKoqAiAt956ixEjRrS7AfFu2KAM5k+5EICDJfUsf2tvXPWimqT+kxDdXpv+lPvoo49YvXo1AJmZmfz2t79l3rx5rR7Tu3dv7rnnHhYuXEgwGGTBggWMHj2aRYsWsXjxYkaNGsXjjz/OAw88gNfrJTc3l1/96lfn36I4dNnQbJp8Jq+++wm7DlWzYv1Bbpg0+LRTebHKtGwCpo0rDnpJQsSrNiWKYDBIIBDA5XIBofGFtpgzZw5z5sw57bY//vGP4ctjxozhxRdfbGusPdrE0X1o8AZ5e+txPigqIznRybQr+kc7rPOmVGgBntvpoodszSFEzGlTovjc5z7HXXfdxbx589A0jdWrVzNlypRIxyY+peDK/jR6g2wpLmfdR8dITnQyYWRutMM6b/6ghWkpDD32e0hCxKM2JYrvfe97LF26lDfffBOHw8F1113HTTfdFOnYxKdozbvkNfpMig5X88rGQyQnOhl9YWa0QzsvdnP5cU8czegSIp5oKob34txzoDyuZz1VVTWe9b6AafHnV3ZzuLQBQ9e4fUY+F/ZN7eIIz8+n2+dwaGT1SoxiRJ0rnmfOxHPbIL7bF9FZT2+88QZTp05l7NixXHbZZeEvER0uh8HCgnxy0hOxbMXza4s5XnH2pBIrLEsRMOMz6QsR69rU1//1r3/Nf//3fzN8+PC4mGkTDzwJDm6fkc+TK3ZS2xjgL68V8bV5I8jolRDt0DoktKe2icvhinYoQohPaVOPolevXkybNo1+/frRt2/f8JeIrrRkN7fPyCfBZVDvDfL0a0U0+mK31IcvYGHH7plQIeJWmxLFmDFj+Pe//x3pWEQH9M7wsHB6Hg5Do7LWx7OFewjE6L7Utq3wxWjsQsSzNp16+ve//83zzz+P0+nE6XSG6zZ9+OGHkY5PtMGg3F58cerF/PX1Yo6UNfDXN/by/wryYnK6qddn4nEZQOzFLkS8alOieOaZZyIchjhfIwdnMGfiIFauP8SeIzWseOcAN0weEnNjSqZlE7QUTiO24hYinrXp1FPfvn3Zvn07L7zwAhkZGXz00UcyRtENjRuey+cuDf1cPthTzptbjkY5ovYLDWrL6SchupM2JYolS5bw17/+lcLCQnw+H//3f//H448/HunYRAdcd3k/LhuaBcBbHx5j8+7Szzii+/H7TWRIW4juo02J4pVXXuGPf/wjiYmJpKen88ILL4SLBIruRdM0bpg8hIv7hRbgvbz+IEWfVEc5qvYxbUUwThdSChGL2pQoHA5HuCAghKbLOhxSbqG7MnSdW64bSt+sJJSCv76xlyNlsbXS1BeU8uNCdBdtShR9+vRh3bp1aJpGIBDgiSeekDGKbs7tNFg4PY+MFDdBy+YvhXuoqPVGO6w2kzUVQnQfbUoUP/jBD3j66afZs2cPl1xyCW+//TY/+MEPIh2bOE8pHhe3z8jH43bQ5DN55tUi6psC0Q6rTWxbSnoI0V20qyig1+vFsiySk9tfVCoSempRwPY6XFrPn1bvJmjZ9M1O4iuzh+N2Gp3y3B3VlvYlOA3Se7ljcp+KeC4sF89tg/huX0eLArZpoOHpp58+6+133HFHu19QdL0BvVO46dqLeX7tHo6VN/LXN/ZyW8FQDL177yoXMG1My+72cQoR79qUKIqLi8OXA4EAmzdvZvz48RELSnS+YQPTmXv1YFasP0jxkRpefucg87v5gjxbKfxBG49bEoUQ0dSmRPHQQw+ddr20tJT7778/IgGJyLlqeG/qGgP866NjbNlTTmqSi2sv797bqXr9Jh63lPQQIpo69Kda7969OXbsWGfHIrrAtZf347Kh2UDzgryisihH1LqWkh5CiOhp9xiFUoodO3aQmRnb22/2VKEFeYOpbwqw92gtK945QEqik/yB6dEO7axaSno4E+X0kxDR0qZ3X3Fxcfhr79699OnTh0ceeSTSsYkIOXVBnh0DC/L8fjMmZz4JES9kz+xuqjOnx55LfVOAP6zYSXW9H0+Cg7vnjSArtWv2rW5v+zJS3LiiPKW3PeJ5imU8tw3iu30RnR572223tTo75tlnn233C4voS/G4uGNmPn9YsZMmn8nTrxZx97wRpHi633ak3oAVU4lCiHjSplNPI0eOJCEhgYULF3LXXXeRlZVFWloat956K7feemukYxQRlJWayJen5+E0dKrr/fylcA/+bljm2x+Ukh5CREubehQffvghy5YtwzBCf9FNmjSJL37xixQUFEQ0ONE1+uekcHPzgrzjFY0se6OY2wrycBjdZwDZthX+oEWiS4pRCtHV2vRJUFVVhd/vD19vbGzE5/NFLCjR9fIHpjNv0hAA9h6t5aW3D9Ddhq98fotuvD5QiLjVpj/PZs+ezZe+9CWuu+46lFK89tprLFy4MNKxiS52RX4OdY0B3txylI/2VpDicTL9qoHRDissYFoELRuHlPQQoku1KVF8+9vfZvjw4bz77ru43W5++tOfcuWVV0Y6NhEFUy/rS31TgPd3l/H21hJSPC6uHtUn2mEBoTUV/oCFI0EShRBdqc3vuN69e3PxxRfzne98B6fTGcmYRBRpmsbcqwczfFBoAd4rmz5h676KKEd1kjcgGxoJ0dXalCj+8Y9/8P3vf5+nnnqK+vp6vvGNb/DCCy9EOjYRJbqu8aWpFzMwNwWAF9ftZ9/R2ihHFWJZsk+FEF2tTYni+eefZ/ny5SQnJ5OZmck///lP/vKXv0Q6NhFFTofOwoI8ctITsWzF86/v4Vh5Q7TDai7pIb0KIbpSmxKFruunbVbUp0+f8FRZEb8S3Q7umJFPapKLQNDmmdeKusV2qr6AhaJ7zcgSIp61KVGkpaWxe/fu8OrslStXkpqaGtHARPeQmuzmjpnDSHQ7aGxevV0X5e1UQ2sq5PSTEF2lTYnivvvu47vf/S779+9n4sSJ/Pa3v+WBBx6IdGyim8hJT+T2GXk4HaHV28+8WoTXH93TPz6fKWsqhOgibUoUPp+PFStW8NJLL/HnP/+ZwsJC8vLyPvO4VatWMXPmTKZNm8bSpUvP+bh169YxderUtkctulz/nBRuvW4ouqZxoqqJZ9fsIWBGr9RHwLSxLOlVCNEV2pQo/uu//gvDMLjwwgsZOnRom6bHlpaW8uijj7Js2TJefvllli9fzr59+854XEVFBQ8//HD7Ixddbmj/NBZccyEAn5yo529v7MWyo/NhbStFg08GtYXoCm1KFHl5eaxatYrjx49TU1MT/mrNxo0bGTduHGlpaXg8HgoKCigsLDzjcQ888ADf/OY3Oxa96HKXXJTF7AmDACg6XMM/1h2IWrE+r98kKL0KISKuTSuz33zzzTM+5DVNY/fu3ec8pqysjOzs7PD1nJwctm3bdtpjnn32WYYPH86YMWPaE3NYaqoHy47f2S8ZGUnRDuGsZk++EKVpvLLhIB/vqyA9NYEvXju01VL0Z9MZ7XM6DLLSEtr92l0hOzsl2iFETDy3DeK/fe3VpkSxffv2dj+xbdunvXmVUqddLy4uZu3atTzzzDOcOHGi3c8PUFvbJBsXRcmE4TlUVDfx3q5S/rXlKDrw+bH92nx8Z7VPAwJef7fbqyKeN7+J57ZBfLevoxsXtXrq6Qc/+EH4clVVVbueODc3l/Ly8vD18vJycnJywtcLCwspLy/nxhtv5Ktf/SplZWXccsst7XoNET2apjHn6kGMuSi0d/qbW46yYXtJl8ehgHpvQFZVCBFBrSaKHTt2hC/fdddd7XriCRMmsGnTJqqqqvB6vaxdu5bJkyeH71+8eDFr1qxhxYoVLFmyhJycHJYtW9bO8EU06ZrGgs9dSP6ANCBUF+rD4vLPOKrzBU0lq7WFiKBWE8Wp+xG0d2+C3r17c88997Bw4UKuv/56Zs+ezejRo1m0aFGHTmWJ7snQdW6+diiD+4TO6f7j3/vZfqCyy+No9AW7/DWF6CnavF1YRwYL58yZw5w5c0677Y9//OMZj+vXrx9vvfVWu59fdA9Oh85tBXn8+ZXdHC1v5IW39uFy6OQNSO+yGCxTEQjauJxSglyIztbqu8q2bWpra6mpqcGyrPDltkyPFT1LgsvB7TPy6d1cRHDp68UcON51FWcV0OQPymptISJAU62cU8rPz0fTtLOedvqs6bFdYc+Bcpn11M3UNwVYsmoXlbU+XE6dO2cOY0DvM6caRqJ9uqaRmerG6AY74MXzzJl4bhvEd/s6Ouup1VNPRUVFHQ5I9EwpHhd3zRrGkpU7qWkI8MxrRdw1axh9s9v/y9letlJ4AxbJsgOeEJ1K3lGi06Ulu7lr1nB6eZz4AhZ/frWIksqu6R35pKyHEJ1OEoWIiMzUBO6cPZykRCdev8mfX9lNWXXk97IwbYUvGL1ihULEI0kUImJy0hK5a9bJvSz+tHoX5TWRTxZNPhNkCZ4QnabN02OF6IjcDA93zRrGU6t3Ue8N8tTqXSyaPfy86jztOVzNO1uPU13vJz3FzaQxF5w2FTcYtPAHbdzdrKyHaJ9t+ysofO8wFbU+slITmH7VAEZfmBXtsHok6VGIiLsgK4m7Zg0jwWVQ3xRKFmXVTR16rj2Hq1m54SB13iAJbgd13iArNxxkz+Hq8GMU0OANEKWitqITbNtfwdLXi6lpDOBJcFDTGGDp68Vs218R7dB6JEkUokv0zU7mzpnDcDsN6pqC/O+yDzu0//Y7W49jGDouh4GmabgcBoah887W46c9LmgqvEEZ2I5Vhe8dxjB03M7Qz9ntDP2cC987HO3QeiRJFKLL9MtJ5s5Z+bidBjX1fp5a1f4xi+p6P07j9F9bpxHaovXTGpuCUdsrQ5yfilofLsfpP2eXQ6ei1heliHo2SRSiS/XPSeHOWfkkuEM9i6dW7aKsHckiPcV9xmZFQcsmPcV9xmMtW9Eo02VjUlZqAgHz9J9zwLTJSk2IUkQ9myQK0eX656TwnZsuC41ZeEPJorSqbWMWk8ZcgGXZBEwLpRQB08KybCaNueCsj/f6zKht1yo6bvpVA7AsG38w9HP2B0M/5+lXDYh2aD2SJAoRFYP69GqeOmvQ4A3yx9W7OF7x2Yvy8gakM/fqwfRKdOLzm/RKdDL36sHnLEBoK0W9V2pAxZrRF2Zx63VDSUty0eQzSUtycet1Q2XWU5S0Wuupu5NaT7GrpX3HKxr586u7afKZJLgM7piZT/+czt2GUtMgPdndpbvgxXO9oHhuG8R3+yKyw50QkXZBVhKL5gwnJbG53McrRRwsqevU11AK6ryBTn1OIXoSSRQi6nqne1g0dzipSS78QYtnXi06bV1EZzBNJZsbCdFBkihEt5CVmshX5w4ns1cCQcvmuTXFbNvfuTvlNXplYFuIjpBEIbqN9JQEvjp3OLkZHmylWP7mXjYXlXXa89tKUd8UROpACdE+kihEt5LicfGV2cPpn5OMAl56+wDrPjrW7j3bz8UXsPAGpC+f120AACAASURBVLqsEO0hiUJ0O54EB3fOGsZFfVMBWLv5CK9s+qTTVlnXNwWx5RSUEG0miUJ0S26nwcLpeYwakgnAxh0neOGtfZjW+X/A27aiTk5BCdFmkihEt+UwdL70+YsYN6I3ANv2V/KXwiJ8gfMvy+ELWDTJKSgh2kQShejWdE1jzoRBXHd5fwD2H6tjycpd1Dae/7qIhqYgpi29CiE+iyQK0e1pmsY1l/VlwecuRNc0TlQ18YeXd3CijfWhzsW2FXWdkHCEiHeSKETMuGxoNl+ekYfbaVDbGODJFTspPlJzXs8ZCFrUe2UhnhCtkUQhYsrF/dL46tzh9Gpexf1sYRHv7So9r+ds8gUJBGW8QohzkUQhYk6fzCS+fv1ILsj0YCtYsf4gr2w6hN3B8QaloK4xgCXjFUKclSQKEZNSk1wsmjuCYQND5cU3bD/Bs2v2dHhGlGkrahsDsiOeEGchiULELLfT4NbrhjJxdB8Aio/U8MTLOzq0FzeExitqGyRZCPFpkihETNN1jZnjBnLjlCEYukZ5jY/fv7SDvUc7Nsjtb04WShbjCREmiULEhbF5OSyaM5zk5n0tnnm1qMM1ovxBi5oGmTYrRAtJFCJuDOidwjduGEm/7CQUoRpRS18vxt+BFdj+gEVtY0D6FUIgiULEmbRkN4vmjODyvGwAdh2q5vGXtlPagcV5Xr9JfVMAqQklejpJFCLuOB06N0wewvWTBmPoGhW1Pn7/8g4+Ki5v93M1+UwafSaaFoFAhYgREU0Uq1atYubMmUybNo2lS5eecf8bb7zBvHnzmDt3Lt/4xjeora2NZDiiB9E0jSuH9earc0eQluwiaNr8fd1+Xnr7AEGzfRVoG7xBGn3nX4hQiFgVsURRWlrKo48+yrJly3j55ZdZvnw5+/btC9/f0NDAj3/8Y5YsWcLKlSvJy8vjd7/7XaTCET1U/5xkvjl/FHn90wDYXFTGEy/voLS67aeilIL6pgBNnVC1VohYFLFEsXHjRsaNG0daWhoej4eCggIKCwvD9weDQX70ox/Ru3eohHReXh4lJSWRCkf0YJ4EJ7dNz6Pgyv7oGpyoauL3/9zB5qKyNs+KUgrqGwN4JVmIHihiiaKsrIzs7Ozw9ZycHEpLT9bkSU9P57rrrgPA5/OxZMkSrr322kiFI3o4XdOYcklfFs1pPhVl2bz09gH++sZemnxtKwrYUurDG5AxC9GzOCL1xLZto53yblJKnXa9RX19Pf/xH/9Bfn4+N9xwQ7teIzXVE9f1eTIykqIdQkRFo30ZGUnkDcnk+deK+HBPGTsOVnGkvJEvzxrG8MGZbX4ew+WgV7IbQz93xsjOTumMkLuleG4bxH/72itiiSI3N5cPPvggfL28vJycnJzTHlNWVsZdd93FuHHjuO+++9r9GrW1TQSC8bn3cUZGElVVjdEOI2Ki3b4bJw9mUO9kVm86RG2Dn8eWf8z4kbkUXNkfl8P4zOOrAJfToFeSC8dZkkV2dgrl5fWdH3g3EM9tg/hun65rZGYmt/+4CMQCwIQJE9i0aRNVVVV4vV7Wrl3L5MmTw/dblsXdd9/NjBkzuP/++8/a2xAiUjRN4/L8HL5142gG9A69cTbtOMHv/rGdT0607UMiELSoqvPR6Del5IeIaxHrUfTu3Zt77rmHhQsXEgwGWbBgAaNHj2bRokUsXryYEydOsGvXLizLYs2aNQCMHDmSn//855EKSYgzZPZKYNGcEbz98XHe+vAolbU+lqzcycTRfbj28v44Ha3/LWXbKjTI7dXwJDrxuCP2lhIiajTVkWI43cSeA+Vy6ilGdcf2bdpZwpr3jhBoXmeR4nEyYUQue4/WUF3vJz3FzaQxF5A3IP2cz+EwNJI9Lvr1SaWioqHdMWzbX0Hhe4epqPWRlZrA9KsGMPrCrA63qTOtXH+AtZuP4gtaJDgNpl3Rj7kTh0Q7rE4np57OclwEYhEi5uw5XM36bSWkJDlJTnQCUN8UZM3mIxyrbMLlNKjzBlm54SB7Dlef83lMS1FT76e82tv8R0zb/w7btr+Cpa8XU9MYwJPgoKYxwNLXi9m2v+J8m3feVq4/wMqNh/AHLRx6qHDiyo2HWLn+QLRDE11AEoUQwDtbj2MYOm6ng15JLrLTEmkZNfMHLMprfJimQtc13tl6/DOfL2jZVDf4qKzz4w2YbdrjovC9w80xGGiahttpYBg6he8dPs/Wnb+1m4+ioWHoGpqmh76jsXbz0WiHJrqAnFAVAqiu95NwyvhCy9hES7KwlaKmwY/LoRNs4+lOpSBo2tQ2BNB1jUS3A4/bwNDP/vdZRa0PT8Lpb0mXQ6ei1tf+BnUyX8A8YyqwrtHhHQVFbJEehRBAeoqboHV6AnAYGg5DIyc9kQRXaMpswLSpawpS+N4n7SpfbtuKRm+Qylo/dU0BTEudsWgvKzUhPD7SImDaZKUmdKxRnSjB5eDTS5ZsFbpdxD9JFEIAk8ZcgGXZBEwLpRQB08LlNHC5HFhKkZ7ipleSk5Y/qt/eWsL/Lv+YD4vL27V1qq0UTT6Tyjov1XX+5sQQOn76VQOwLBt/MBSDP2hhWTbTrxoQgRa3z7Qr+qFQWLZCKTv0HcW0K/pFOzTRBYwf//jHP452EB1VWd0UtyuzExNdeL1tKy0Ri7pb+7JSE8lKTaC0somGpiBpSS5mjBvI8EHp4dsyeyUwY9xAemd4OFLagDdgsetQNXuO1JCdlkh6ijv8fG1pn2krfAGTQPPYR9+sZHqnJ3K0rIHaxgAZKW7mTx7SLWY95Q1IB6X45EQDAUuR4DSYOW5AXM56Skpy09QUnzscapqGx+Nq/3EyPbZ76o7TRztTrLevpsHPa+8eZvuByvBtwwelM/3KAWSlJXaofU6HTlKiMzSY3dkBd6J4nj4K8d2+jk6PlROMQnRAWrKbm6+9mAkncnn13U84UtbArkPVFH1Szdi8HOZPvbjdzxk0bWrq/TgMjcQEJwlOHYehE7t/yol4IYlCiPMwMDeFu+eNYPuBKta8f5jqej+bi8r4aG8F40b0ZvKYC8LrMtrKtEKrvRs1DZdDJyHBgdPQcRiaJA0RFZIohDhPmqYx+sJMhg9KZ3NRGf/68BgN3iDrt5Xw3q5Sxo/ozcTR7U8YtlL4gha+oIWuaTgdOm63gcuh4zjHFFshIkEShRCdxGHojB+Ry9ih2Xx8oIo17x7C67d4e2sJ7+4s5arhvbl6dB96dWAw0W6eBeUPWmgaGIaG2+nA7dBxOEIL4KS3ISJFEoUQnczlNJg+fhCjB6ezcccJ1m8rwReweGdbCZt2nmBsXg6TRvcho1fH1kcoBaapMM0gjYQGKF2Gjstl4DB0nA4NXZPEITqPJAohIiTB5WDqZf2YMDKXTTtK2bCjhCafyXu7Snl/dykjB2cyaXQf+uW0fxbKqWxb4bNDp6g0QNO15oQR+nLoGoahd+uZVKJ7k0QhRIQluBxcc1lfrh6Vy+aiMt7ZVkJdY4DtByrZfqCSQbkpTBjVh2ED01vdMa8tFKBsRcC2CARDK8c1DQxNw+UycDkNnIYms6lEu0iiEKKLuJwGV4/qw1XDe7N9fyXrt5dQUtnEoRP1HDpRT1qyi6uG9+aK/Bw8Ce0b+G6NUmAqhekzafKZ6JqGboDb6Qj3OByGjqYhyUOclSQKIbqYw9C5dGg2l1ycxf5jdWzYUULx4RpqGgKsef8Ib245ysjBmVw5PIeBvVM6ffdHWylsE0wztHJc00IztwxDw+UwcBgahq43f5exDiGJQoio0TSNi/qlclG/VCprfby78wQf7CnHH7T4eF8FH++rICc9kcvzcrjk4qx2T69tK6VAKYVtq3Bl3Jbkoevg1A2cTj00LdeQabk9kZTw6KZivcTFZ5H2nV0gaLFtfyXv7S7lWPnJ43VNI39gGpcNzWZo/7SofGBrGhi6Rm5OLxoafBh6qMeha1pcnbaSEh5nkh6FEN2Iy2lweX4Ol+fncKyikQ+Kyti6rwJfcwHCXYeqSXQ7GDUkg0suzmJA7xT0Tj41dS5KhVaNN/lNaur9oV4HoSShGxoOTUc3QgsDDT1Uor0ltnhJIj2VJAohuqm+WUn0nTiYmeMGsutQFR8Wl7PvWC1ev8n7u8t4f3cZqUkuRg3JZNSFmfTLTur08YzWKAUKBQosWxHkZO/+1FNXRvOOeIYRmqaraxq6FvrrVtMApTU/XhJKdyWJQohuzunQGXNRFmMuyqKuKcC2fZV8vK+C4xWN1DYGWL+9hPXbS0hNcjF8cAYjBmUwKDcF/Tyn2p6Pk+MeYHLmBk8tvZGWxR26Bg5dx9FcCFEnNDMrlFS08DHnei0RWZIohIghvTwuJo7uw8TRfaio8bLtQCXb9ldSVu2ltjHAph0n2LTjBB63g7wBaeQPTOfifqndbie6U3sjADZgWhYETyaVU5OJBqHvn0oWmqZhaC3jJJ++85SeTfNaknAvpvkBLT2ZllNo4uy612+PEKLNstISmXpZP6Ze1o+yai+7DlWx82AVxyoaafKbfLS3go/2VqBrGgNzU8jrn8bQAWn0Tk/s0lNUHfXpZHKOR9Ge7a9ObbbW0p1pTkSaFlrVrgyD2oZA+HpLsqG5h2Q37+7Xkpw0vbnno2uhLUObE5PWPMgfejwoO/Rap96nlDqtR6Q1Z0RNO3n50x3DaPSgJFEIEQdy0hPJSe/L5y7tS3W9n6LDob0xDhyvw7IVB0vqOFhSR+H7h0nxOLmobyoXNn+lJrW/SGGsOvVDVrVkoFM/eC1FwLTxBswOv0a499OcgVomlir1qftOjeG040+ekmt5Pk1v6fU0J6GWZHJKItPCCe9k7yjcXqWwVfOspw60SRKFEHEmPcXN+BG5jB+Riz9gsf94LXsO11B8pIbaxgD1TcFwbwMgMzWBIX16MeSCXgzKTSE12f0ZryBao5r/OVtv6LT7znn8WXpRdvjo8+J0dGxatSSKDtDOuNAG6qwXhYgot8tg+KAMhg/KQClFeY2Pfcdq2Xe0lgMltQSCNpW1PiprfWwuKgMgLdnFoNxeDMhNZkBOCr0zPOddg0rEtphOFC6ncXJWxKndrzP6XWf/cD75ga+d0iVs6R6ePtNCU6c84ORh58wVp75eq+cU1dmvpCW7UObJ7u/ZnkOp5nOf6pQHnHISVgv/czLOU+NV57h8ao/80+dQP+2sUxrDx5/9/LJS4HYaJDYPsKpzZNFPL6U8/VTBye68Otdp7FP+cpOZMaH3RegUVSITRuZi2TbHyhs5WFLHgeN1fHKinoBpU9MQCK8MB3A5dPpmJ9EvO5l+OcmM0HU0pWJinEN0jphemV1V1YBlxWz4rYrn1aFw7vZ17LNHw1Ynk4imheoZqeaEYiswLZtg0CZgWtgtD/2MX32tlS6jovUEGosrzy1bUVLZyCcn6jlUUs/hsnrqm84+VOxxO7ggK4kLspLok+mhT2YSWakJUZ2S21li8WfXVk6HTv6F2e0+LqZ7FLGb4sS5dOxnqs44HaiHu3saBuA0dBJdJ18j9EF/thknpzrzQ6+lB2XZCsu2MU1F0LKxLLtNPbDuzNC1UK8hO5mrR/VBKUVtY4DDpfUcLWvkSHkDx8sbCVo2TX4zdArrWG34eKehk5ORSG66h94ZHnpnJJKTlkivJJf0PmJcTCcKITri5Pz8Dn54Ndc8Ah1cJ5OH3VxYz7QViS4HDkPDsmM3cWiaRlqym7RkN6MvzAJCCdJvK4oOVHKsopHj5Y2cqGrCH7QIWqFTWafWqILQacbstASy0xLJSk0kq/lyRi83LocRjaaJdpJEIcR5akkEoX0eNBwGZKQmYAWCWJZNwLIJBEKnveKi15GVjMehc9nQ0CkMWylq6v2UVDZxoqqJ0uomSquaqKz1YSvwBy2OljdytPzM0zm9klxk9nKT0SuBjJQEMnq5SU8JfSUnOqUn0k1IohAiggxDJ/GU014tvQ7bVljNly1LnRxjAWwFtm2HxlKaeyrdma5poQ/6XgmMGJwRvt20bCpqfZTXeCmr9lJe46Wy1kd5rTdc9bmuMUBdY4CDJWeOVzmMkz2a1GRX6HuSi15JrvD3BJchyaQLSKIQogu19Do4xxmXUz/zbBtsZROqbGERCNoELTtmeiQOQyc3w0Nuhue025VS1HuD4Wm5lXU+qur8VNf7qKzz4/WHZvuZlqKi1kdFre+cr+F06PTyuEhJcpKS6CLF42z+cpGU4CDZ4yI50UlSgkP20jgPkiiE6EY+PbgeqrwKLqdOcuLJ3oZlExpMtxRB08a07DMWcnXXZKJpGr08Lnp5XAzu0+uM+30Bk5qGADX1fqrq/dQ2+KlpCFDb6Ke2IUB9UyA8cy1o2lTWhZLNZ0lwGSQlOElKdOBxO/AkOPEktFwOfU90O2gM2gR8QRLdRngKfk8niUKIGNFSAsLQQ8mDUGWh5unALWMfqnlgnZidlZXgcpCb4TijJ9LCthUNviB1zUmjtilAfWOQem+QhqbQyvMGb+jLsk821hew8AUsKuvaHoumhRJMgsvR/P3kZbfLIMEZ+u5yGrhP+XI59ebvocstW8zG6mmyiCaKVatW8cQTT2CaJl/+8pe59dZbT7t/9+7d3H///TQ2NnL55Zfzk5/8BIdDcpcQ7dGSQLTmSnLb9ldQ+N5hKmt99M5IpOCqAYwYlMnza4vYtr8KwwjtDzFsQCp5AzPYsO04VXV+UjxOxo3I5eL+aRR9Us07W49TXe8nPcXNpDEXcKy8gfXbTuA3LdwOg4mjc5k6tv8549pz+MznAM64LW9AepuPzxuQzrqPjobiCFq4naE4rrvizDiUUmw/UMn6bcepaQiQ5HYw+IJUUjwujpY3cKSsHl/AwtA1XE4D07Lx+a0zFst6/RZe/5ml0ttL08DlMHA59OatZQ2cDv3klxH67jBO3uY49TZDw2HoGM2XDSN0uyN8e/N3/eR3o3n/8/NdWR+xBXelpaXcfPPN/POf/8TlcnHTTTfxv//7v1x00UXhx8yePZsHH3yQSy65hPvuu4+RI0dyyy23tPk1KisbsO1u/udRB/XUBXfxIlrt27a/gqWvF2MYoT2uA2aoN5GZ4qboSO1pj3UYOokujaz0JBKcOqYd6pFcnpfNx/sqcRgGhg5+06K23o8vYNO8Xh/TVFgKJo/uw5RL+57RS9lzuJqVGw42f6jpBC0br88ETSPRbYRvsyybuVcPPiNZnO14y7IZmJPM1gNVaISqsNrNq/I/f1nfM5LWuZ5j7NBsthSXn3H73KsHc3H/NBKT3Bw/UYfPb7L3aA0btp8IV3wwLRvbVvTLScbtNPAFLAJBC3/Qwh+08Tdf726fShpgGBq5mUk8ce/n2318xP5837hxI+PGjSMtLQ2AgoICCgsL+eY3vwnAsWPH8Pl8XHLJJQDMnz+fxx57rF2JIh5WgbZG2hfbotG+93aVkpPhOW19QsC0KK/2kpOeeNYSLmmnFAEMmBabdpSRnOTE6QjtRudyGjR6TVIcDtwuo7mUdujxx6uayEpLDM3O4uSpreOVTeQPSsfQ9dDYiVJU1vpQQFqKO1x6Jmha7D9Wy+iLsk4LqvhIDf1yknEYRniFfcC0KKv10SfTc/L/ViksW1F8pJaCqwae9n+x82DVWf8v9hypPevtOw9WMXJIJr08Luzm017v7SqlX+/kMx6bkuDkS5+/+Kw/A0VoxlcgaBE0bQKmjWmGpkcHm08FmqZF0FKYpk3QDE1SMFuuWxamGVqPY5o2lgodb1mKoK2wLBur+b72zohL8Tjb9fgWEUsUZWVlZGefXCqek5PDtm3bznl/dnY2paWl7XqN9PSk8w+0G+vIJuixRNrX+e67c1yXv+bZLL7psvM6/r4hWecdw486UKqiRW5m0nk/RzyJ2Hwx27ZPG7hRnyoi9ln3CyGE6B4ilihyc3MpLy8PXy8vLycnJ+ec91dUVJx2vxBCiO4hYoliwoQJbNq0iaqqKrxeL2vXrmXy5Mnh+/v27Yvb7WbLli0ArFix4rT7hRBCdA8RLTO+atUqnnzySYLBIAsWLGDRokUsWrSIxYsXM2rUKIqKinjggQdoaGhgxIgRPPTQQ7hcPWdbRiGEiAUxvR+FEEKIyJPiJ0IIIVoliUIIIUSrJFEIIYRolSQKIYQQrYqZRPHb3/6WmTNnMmvWLJ5++mkgVCZkzpw5TJs2jUcffTTKEZ6/hx9+mP/+7/8GQgUT58+fT0FBAffffz+maUY5uo677bbbmDVrFvPmzWPevHls3bqVVatWMXPmTKZNm8bSpUujHeJ5eeutt5g/fz4zZszgwQcfBOLnd/Pvf/97+Oc2b948xo4dy09/+tO4aR+EpubPmjWLWbNm8fDDDwPx8/5bsmQJBQUFzJkzhyeeeALoYNtUDHjvvffUTTfdpILBoPJ6veqaa65Ru3fvVlOmTFGHDx9WwWBQ3XnnnWrdunXRDrXDNm7cqK666ip17733KqWUmjVrlvroo4+UUkp9//vfV0uXLo1meB1m27aaOHGiCgaD4dtOnDihrrnmGlVdXa0aGxvVnDlz1N69e6MYZccdPnxYTZw4UZWUlKhAIKBuvvlmtW7durj63WxRXFysrrvuOnX8+PG4aV9TU5O64oorVGVlpQoGg2rBggVqw4YNcfH+27Bhg5o9e7aqr69Xpmmqr33ta2rNmjUdaltM9CiuvPJKnn32WRwOB5WVlViWRV1dHQMHDqR///44HA7mzJlDYWFhtEPtkJqaGh599FHuvvtu4OwFE2O1bQcOHADgzjvvZO7cuTz//POnFYz0eDzhgpGx6PXXX2fmzJnk5ubidDp59NFHSUxMjJvfzVP9+Mc/5p577uHIkSNx0z7LsrBtG6/Xi2mamKaJw+GIi/ffrl27mDhxIsnJyRiGwaRJk3juuec61LaYSBQATqeTxx57jFmzZjF+/PizFh1sb1HB7uKHP/wh99xzD716hXb76oyCid1FXV0d48eP5/HHH+eZZ57hb3/7G8ePH4+bn90nn3yCZVncfffdzJs3j2XLlsXV72aLjRs34vP5mDFjRly1Lzk5mW9/+9vMmDGDKVOm0LdvX5xOZ1y8/0aMGMH69eupqanB7/fz1ltv4XA4OtS2mEkUAIsXL2bTpk2UlJRw6NChuCgq+Pe//50+ffowfvz48G3xVDDx0ksv5Ve/+hUpKSlkZGSwYMECHnvssbhpn2VZbNq0iV/84hcsX76cbdu2ceTIkbhpX4u//e1v3HHHHUB8/X4WFRXxj3/8g3/961+888476LrOhg0b4qJ948ePZ/78+dx222185StfYezYsZim2aG2xcR2cvv37ycQCDBs2DASExOZNm0ahYWFGMbJGvGfLjoYK1599VXKy8uZN28etbW1NDU1oWla3BRM/OCDDwgGg+FEqJSib9++rRaMjCVZWVmMHz+ejIwMAK699tq4+d1sEQgE2Lx5M7/85S+Bzy74GUvWr1/P+PHjyczMBEKnYv70pz/FxfuvoaGBadOmhRP8U089Rb9+/fjggw/Cj2lr22KiR3H06FEeeOABAoEAgUCAN998k5tuuomDBw+Gu/6rV6+OyaKCTz/9NKtXr2bFihUsXryYqVOn8tBDD8VNwcT6+np+9atf4ff7aWho4KWXXuLXv/51qwUjY8k111zD+vXrqaurw7Is3nnnHaZPnx4Xv5st9uzZw6BBg/B4Qpv5jBkzJm7al5+fz8aNG2lqakIpxVtvvcWVV14ZF++/o0eP8o1vfAPTNKmvr+fFF19kwYIFHWpbTPQopkyZwrZt27j++usxDINp06Yxa9YsMjIy+Na3voXf72fKlClMnz492qF2mkceeeS0gokLFy6Mdkgdcs0117B161auv/56bNvmlltuYezYsdxzzz0sXLgwXDBy9OjR0Q61Q8aMGcNXvvIVbrnlFoLBIFdffTU333wzQ4YMiZvfzSNHjpCbmxu+7na7+eUvfxkX7Zs4cSK7du1i/vz5OJ1ORo0axVe/+lWuu+66mH//5efnM23aNObOnYtlWdx+++2MHTu2Q58tUhRQCCFEq2Li1JMQQojokUQhhBCiVZIohBBCtEoShRBCiFZJohBCCNGqmJgeK0RbPPjgg2zevBkILdLs27cvCQkJACxfvjx8uTtSSnHHHXfw2GOPhUu5CNFdyPRYEZemTp3Kb3/7W0aNGhXtUNrENE1GjBjB5s2bJVGIbkd6FKJH2Lt3Lz//+c/DK6hvv/12brjhBjZu3Mjvfvc7srOz+eSTT/B4PHzlK1/hueee49ChQ8yYMYN7772XjRs38thjj5GTk8PBgwfxeDw89NBDDBkyhEAgwK9+9Su2bNmCZVmMGDGC+++/n+TkZCZPnszYsWMpKiriu9/9LrZt89RTTxEIBKiqquLGG2/kW9/6Ft///vcBuPXWW3nqqaf4whe+wJNPPsmwYcMAmDx5Mk8++SQej4c77riDAQMGUFJSwrJlyzh48CD/8z//g8/nQ9d1Fi9ezJQpU6L53y3iTedXQRci+q655hq1bds2pZRSgUBAzZgxQ+3evVsppVRtba0qKChQ27ZtUxs2bFDDhw8P33f77berm2++WQUCAVVRUaGGDRumKioq1IYNG1R+fr7asmWLUkqp5557Tn3hC19QSin1m9/8Rv36179Wtm0rpZR6+OGH1c9+9jOllFKTJk1Sf/jDH5RSSlmWpW655RZ1+PBhpZRSx48fV/n5+aqmpkYFg0E1dOhQVVtbGz5u165d4fa0XD906JAaOnSo+vDDD5VSSlVVValp06apY8eOKaWUKikpUZMmTVIlJSUR+p8VPZH0KETc279/P0eOHOHee+8N3xYIBNi9ezf9+vVjwIAB5OfnA9C/f3+ysrJwOp1kZmbi8XioqakBQmWbL7vsMgC+8IUv8OCDD1JfX8+6detoamrinXfeASAYDJ5WaG3s2LEA6LrOk08+ybp1rDlRwwAAAjRJREFU61ixYgX79u1DKYXP5yMpKanN7XE6nYwZMwaADz/8kPLycr7+9a+H79d1neLi4tPKbghxPiRRiLhn2zZpaWmsWLEifFt5eTm9evViy5YtuFyu0x7vcJz9bXHq7bZtA6EPZcuy+OEPf8jVV18NhKp2BoPB8GNbkkBDQwM33HADBQUFjB07lhtvvJHXX38ddZZhQk3TTrv91OdLSEhA1/VwHEOHDuVvf/tb+P7S0tJwNVshOoNMjxVx76KLLkLXdV555RUgtIPg7NmzKSoqatfz7Nixg7179wKhWVRXXHEFSUlJTJw4keeee45gMIhlWdx333385je/OeP4gwcP4vV6+fa3v80111zDpk2bME0Ty7IwDANN08L7F2dkZLBjxw4gtGlQVVXVWWO69NJL2b9/f7ga6M6dOykoKKCysrJdbROiNdKjEHHP5XLxxBNP8Itf/II//OEPmKbJf/7nfzJmzBg2btzY5ufJycnhkUce4dixY2RnZ/Pwww8D8K1vfYuHH36Y66+/PjyY/b3vfe+M44cPH87EiROZMWMGTqeT/Px8hgwZwuHDh+nbty/Tpk3j5ptv5ve//z3f/e53+clPfsLSpUsZNWpUeFD707Kysnjsscd46KGHCAQCKKV45JFH5LST6FQyPVaINti4cSMPP/zwaaevhOgp5NSTEEKIVkmPQgghRKukRyGEEKJVkiiEEEK0ShKFEEKIVkmiEEII0SpJFEIIIVoliUIIIUSr/n+qkd4OsQ4/eAAAAABJRU5ErkJggg==\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.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}