diff --git a/module4/install.sh b/module4/install.sh new file mode 100755 index 0000000000000000000000000000000000000000..422564265c97834d0e24dccfef5b736f6e5765da --- /dev/null +++ b/module4/install.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +python3.9 -m venv venv +source venv/bin/activate +pip install -r requirements.txt diff --git a/module4/repro.ipynb b/module4/repro.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..83b8e6b34c8d228b8790d3b34fc3ce06a784672a --- /dev/null +++ b/module4/repro.ipynb @@ -0,0 +1,801 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "considerable-stanford", + "metadata": {}, + "source": [ + "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure\n", + "(copy with some modifications from https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/src_Python3_challenger__1_.ipynb graphs doesnt give the same output)\n", + "\n", + "---\n", + "\n", + "In this document we reperform some of the analysis provided in *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 \n", + "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. \n", + "\n", + "Our goal is to reproduce the computation behind these values and the Figure 4 of this article,\n", + "possibly in a nicer looking way.\n", + "\n", + "## Technical information on the computer on which the analysis is run\n", + "\n", + "We will be using the python3.9 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "suffering-night", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.9.0+ (default, Oct 19 2020, 09:51:18) \n", + "[GCC 10.2.0]\n", + "uname_result(system='Linux', node='dell', release='5.8.0-44-generic', version='#50-Ubuntu SMP Tue Feb 9 06:29:41 UTC 2021', machine='x86_64')\n", + "IPython 7.21.0\n", + "IPython.core.release 7.21.0\n", + "PIL 8.1.2\n", + "PIL.Image 8.1.2\n", + "PIL._version 8.1.2\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.2.0\n", + "cffi 1.14.5\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.9.0+\n", + "ipykernel 5.5.0\n", + "ipykernel._version 5.5.0\n", + "ipython_genutils 0.2.0\n", + "ipython_genutils._version 0.2.0\n", + "ipywidgets 7.6.3\n", + "ipywidgets._version 7.6.3\n", + "jedi 0.18.0\n", + "json 2.0.9\n", + "jupyter_client 6.1.11\n", + "jupyter_client._version 6.1.11\n", + "jupyter_core 4.7.1\n", + "jupyter_core.version 4.7.1\n", + "kiwisolver 1.3.1\n", + "logging 0.5.1.2\n", + "matplotlib 3.3.4\n", + "numpy 1.20.1\n", + "numpy.core 1.20.1\n", + "numpy.core._multiarray_umath 3.1\n", + "numpy.lib 1.20.1\n", + "numpy.linalg._umath_linalg 0.1.5\n", + "pandas 1.2.3\n", + "pandas.compat.numpy.function 1.20.1\n", + "parso 0.8.1\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.17\n", + "ptyprocess 0.7.0\n", + "pygments 2.8.1\n", + "pyparsing 2.4.7\n", + "pytz 2021.1\n", + "re 2.2.1\n", + "scipy 1.6.1\n", + "scipy._lib._uarray 0.5.1+49.g4c3f1d7.scipy\n", + "scipy._lib.decorator 4.0.5\n", + "scipy.integrate._dop b'$Revision: $'\n", + "scipy.integrate._ode $Id$\n", + "scipy.integrate._odepack 1.9 \n", + "scipy.integrate._quadpack 1.13 \n", + "scipy.integrate.lsoda b'$Revision: $'\n", + "scipy.integrate.vode b'$Revision: $'\n", + "scipy.interpolate._fitpack 1.7 \n", + "scipy.interpolate.dfitpack b'$Revision: $'\n", + "scipy.linalg._fblas b'$Revision: $'\n", + "scipy.linalg._flapack b'$Revision: $'\n", + "scipy.linalg._flinalg b'$Revision: $'\n", + "scipy.linalg._interpolative b'$Revision: $'\n", + "scipy.ndimage 2.0\n", + "scipy.optimize.__nnls b'$Revision: $'\n", + "scipy.optimize._cobyla b'$Revision: $'\n", + "scipy.optimize._lbfgsb b'$Revision: $'\n", + "scipy.optimize._minpack 1.10 \n", + "scipy.optimize._slsqp b'$Revision: $'\n", + "scipy.optimize.minpack2 b'$Revision: $'\n", + "scipy.signal.spline 0.2\n", + "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", + "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", + "scipy.special.specfun b'$Revision: $'\n", + "scipy.stats.mvn b'$Revision: $'\n", + "scipy.stats.statlib b'$Revision: $'\n", + "seaborn 0.11.1\n", + "seaborn.external.husl 2.1.0\n", + "six 1.15.0\n", + "statsmodels 0.12.2\n", + "statsmodels.__init__ 0.12.2\n", + "statsmodels.api 0.12.2\n", + "statsmodels.tools.web 0.12.2\n", + "traitlets 5.0.5\n", + "traitlets._version 5.0.5\n", + "urllib.request 3.9\n", + "wcwidth 0.2.5\n", + "zlib 1.0\n", + "zmq 22.0.3\n", + "zmq.sugar 22.0.3\n", + "zmq.sugar.version 22.0.3\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", + "id": "purple-budapest", + "metadata": {}, + "source": [ + "## Loading and inspecting data\n", + "---\n", + "Let's start by reading data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "aquatic-starter", + "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/raw/master/data/shuttle.csv\")\n", + "data" + ] + }, + { + "cell_type": "markdown", + "id": "steady-seven", + "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, + "id": "diagnostic-optimum", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYAUlEQVR4nO3dfXRc9X3n8fdHtvADdsCxswrFJpDisKHBS0A8eGkauXlYkz2xm4VNoYeQZAPO7uJtQnIaaDaHsmxzzkKb0E2WJnFZmkDbiAenRG3dklCihHQD2A7GgKmJCsTIEAOKMRYYW7K++8dcJSMxku5IujOa+X1e5/h47r2/ufP9zdXoo/swv6uIwMzM0tVS7wLMzKy+HARmZolzEJiZJc5BYGaWOAeBmVniHARmZokrLAgk3STpOUmPjLFckr4kqUfSdkmnFVWLmZmNrcg9gq8Dq8dZfi6wPPu3DvhKgbWYmdkYCguCiPgB8PNxmqwFbo6S+4CjJR1TVD1mZlbZ7Dq+9rHA02XTvdm8Z0c3lLSO0l4D8+bNO33ZsmU1KTCvoaEhWlqa73RLs/YLmrdv7lfjqVXfHn/88Rci4g2VltUzCHKLiA3ABoD29vbYsmVLnSsaqbu7m46OjnqXMe2atV/QvH1zvxpPrfom6adjLatnxO4Gyv+0X5rNMzOzGqpnEHQBF2dXD50N7IuI1xwWMjOzYhV2aEjSN4EOYImkXuAPgFaAiPgqsAl4H9ADvAJ8tKhazMxsbIUFQURcOMHyAC4r6vXNzCyf5jwNb2ZmuTkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcoUEgabWknZJ6JF1ZYflxkr4n6UFJ2yW9r8h6zMzstQoLAkmzgBuAc4GTgQslnTyq2eeA2yLi7cAFwJ8WVY+ZmVVW5B7BmUBPRDwREYeATmDtqDYBvC57fBTwTIH1mJlZBYqIYlYsnQ+sjohLsukPAWdFxPqyNscA3wEWAUcC746IrRXWtQ5YB9DW1nZ6Z2dnITVPVn9/PwsWLKh3GdOuWfsFzds396vx1Kpvq1at2hoR7ZWWzS781cd3IfD1iPiCpJXALZLeFhFD5Y0iYgOwAaC9vT06OjpqX+k4uru7mWk1TYdm7Rc0b9/cr8YzE/pW5KGh3cCysuml2bxyHwNuA4iIHwFzgSUF1mRmZqMUGQSbgeWSTpB0BKWTwV2j2uwC3gUg6a2UguD5AmsyM7NRCguCiBgE1gN3AY9RujroUUnXSFqTNfs0cKmkh4BvAh+Jok5amJlZRYWeI4iITcCmUfOuKnu8AzinyBrMzGx8/maxmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZokrNAgkrZa0U1KPpCvHaPNBSTskPSrpr4qsx8zMXmt2USuWNAu4AXgP0AtsltQVETvK2iwHfh84JyL2SvpXRdVjZmaV5dojkHTKJNZ9JtATEU9ExCGgE1g7qs2lwA0RsRcgIp6bxOuYmdkUKCImbiTdC8wBvg78ZUTsy/Gc84HVEXFJNv0h4KyIWF/W5k7gceAcYBZwdUT8Q4V1rQPWAbS1tZ3e2dk5Yc211N/fz4IFC+pdxrRr1n5B8/bN/Wo8terbqlWrtkZEe6VluQ4NRcQ7ssM4/wnYKukB4M8j4rtTrG02sBzoAJYCP5B0SkS8OOr1NwAbANrb26Ojo2OKLzu9uru7mWk1TYdm7Rc0b9/cr8YzE/qW+2RxRPwE+BxwBfBO4EuS/lnSfxjjKbuBZWXTS7N55XqBrogYiIgnKe0dLM9bk5mZTV3ecwQrJF0PPAb8JvD+iHhr9vj6MZ62GVgu6QRJRwAXAF2j2txJaW8ASUuAtwBPVNkHMzObgrxXDX0ZuBH4bEQcGJ4ZEc9I+lylJ0TEoKT1wF2Ujv/fFBGPSroG2BIRXdmy90raARwGfi8i+qbQHzMzq1LeIPj3wIGIOAwgqQWYGxGvRMQtYz0pIjYBm0bNu6rscQCfyv6ZmVkd5D1HcDcwr2x6fjbPzMwaXN4gmBsR/cMT2eP5xZRkZma1lDcIXpZ02vCEpNOBA+O0NzOzBpH3HMEngdslPQMIeCPw20UVZWZmtZP3C2WbJf1r4KRs1s6IGCiuLDMzq5VqBp07Azg+e85pkoiImwupyszMaiZXEEi6BfhVYBul6/0BAnAQmJk1uLx7BO3AyZFnhDozM2soea8aeoTSCWIzM2syefcIlgA7slFHDw7PjIg1hVRlZmY1kzcIri6yCDMzq5+8l49+X9KbgOURcbek+ZQGkjMzswaXdxjqS4E7gK9ls46lNIS0mZk1uLwniy+jdDvJl+AXN6nxjebNzJpA3iA4mN2AHgBJsyl9j8DMzBpc3iD4vqTPAvMkvQe4Hfib4soyM7NayRsEVwLPAw8DH6d0s5mKdyYzM7PGkveqoSHgz7J/ZmbWRPKONfQkFc4JRMSbp70iMzOrqWrGGho2F/iPwOunvxwzM6u1XOcIIqKv7N/uiPgTSje0NzOzBpf30NBpZZMtlPYQqrmXgZmZzVB5f5l/oezxIPAU8MFpr8bMzGou71VDq4ouxMzM6iPvoaFPjbc8Ir44PeWYmVmtVXPV0BlAVzb9fuAB4CdFFGVmZrWTNwiWAqdFxH4ASVcDfxcRFxVVmJmZ1UbeISbagENl04eyeWZm1uDy7hHcDDwg6a+z6d8CvlFIRWZmVlN5rxr6vKS/B96RzfpoRDxYXFlmZlYreQ8NAcwHXoqI/w30SjqhoJrMzKyG8t6q8g+AK4Dfz2a1An9RVFFmZlY7efcIPgCsAV4GiIhngIVFFWVmZrWTNwgORUSQDUUt6cjiSjIzs1rKGwS3SfoacLSkS4G78U1qzMyawoRBIEnArcAdwEbgJOCqiPhyjueulrRTUo+kK8dpd56kkNQ+VhszMyvGhJePRkRI2hQRpwDfzbtiSbOAG4D3AL3AZkldEbFjVLuFwCeA+6uq3MzMpkXeQ0M/lnRGles+E+iJiCci4hDQCayt0O5/AtcCr1a5fjMzmwZ5v1l8FnCRpKcoXTkkSjsLK8Z5zrHA02XTvdl6fiG74c2yiPg7Sb831ookrQPWAbS1tdHd3Z2z7Nro7++fcTVNh2btFzRv39yvxjMT+jZuEEg6LiJ2Af9uul9YUgvwReAjE7WNiA3ABoD29vbo6OiY7nKmpLu7m5lW03Ro1n5B8/bN/Wo8M6FvE+0R3Elp1NGfStoYEedVse7dwLKy6aXZvGELgbcB3aXz0bwR6JK0JiK2VPE6ZmY2BROdI1DZ4zdXue7NwHJJJ0g6AriAX97PgIjYFxFLIuL4iDgeuA9wCJiZ1dhEQRBjPJ5QRAwC64G7gMeA2yLiUUnXSFpTXZlmZlaUiQ4N/RtJL1HaM5iXPYZfnix+3XhPjohNwKZR864ao21HrorNzGxajRsEETGrVoWYmVl9VDMMtZmZNSEHgZlZ4hwEZmaJcxCYmSUumSDo6z/IQ0+/SF//wXqXYmZV6Os/yIGBw/7sFiiJIPj2tt2cc+09XHTj/Zxz7T10bds98ZPMrO6GP7tPPv+yP7sFavog6Os/yBUbt/PqwBD7Dw7y6sAQn9m43X9dmM1w5Z/dwxH+7Bao6YOgd+8BWltGdrO1pYXevQfqVJGZ5eHPbu00fRAsXTSPgaGhEfMGhoZYumhenSoyszz82a2dpg+CxQvmcN15K5jb2sLCObOZ29rCdeetYPGCOfUuzczGUf7ZnSX5s1ugvDemaWhrTj2Wc05cQu/eAyxdNM8/SGYNYviz+8CPfsg/rfl1f3YLkkQQQOmvC/8QmTWexQvmMK91lj+/BWr6Q0NmZjY+B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpa4QoNA0mpJOyX1SLqywvJPSdohabukf5T0piLrMTOz1yosCCTNAm4AzgVOBi6UdPKoZg8C7RGxArgDuK6oeszMrLIi9wjOBHoi4omIOAR0AmvLG0TE9yLilWzyPmBpgfWYmVkFiohiViydD6yOiEuy6Q8BZ0XE+jHa/x/gZxHxhxWWrQPWAbS1tZ3e2dlZSM2T1d/fz4IFC+pdxrRr1n5B8/bN/Wo8terbqlWrtkZEe6Vlswt/9RwkXQS0A++stDwiNgAbANrb26Ojo6N2xeXQ3d3NTKtpOjRrv6B5++Z+NZ6Z0Lcig2A3sKxsemk2bwRJ7wb+O/DOiDhYYD1mZlZBkecINgPLJZ0g6QjgAqCrvIGktwNfA9ZExHMF1mJmZmMoLAgiYhBYD9wFPAbcFhGPSrpG0pqs2R8BC4DbJW2T1DXG6szMrCCFniOIiE3AplHzrip7/O4iX7+R9fUfpHfvAZYumsfiBXOmrW0jadZ+FaVnz372vjJAz579nNi2sN7lWAOZESeLbaRvb9vNFRu309rSwsDQENedt4I1px475baNpFn7VZSr7nyYm+/bxadPGeTy63/AxSuP45q1p9S7LGsQHmJihunrP8gVG7fz6sAQ+w8O8urAEJ/ZuJ2+/teeR6+mbSNp1n4VpWfPfm6+b9eIeTf/aBc9e/bXqSJrNA6CGaZ37wFaW0ZultaWFnr3HphS20bSrP0qyranX6xqvtloDoIZZumieQwMDY2YNzA0xNJF86bUtpE0a7+Kcuqyo6uabzaag2CGWbxgDtedt4K5rS0snDObua0tXHfeioonS6tp20iatV9FObFtIRevPG7EvItXHucTxpabTxbPQGtOPZZzTlyS64qZato2kmbtV1GuWXsKF599PA9vvY+7Lz/bIWBVcRDMUIsXzMn9y6+ato2kWftVlBPbFtI7v9UhYFXzoSEzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxhQaBpNWSdkrqkXRlheVzJN2aLb9f0vFF1mNWrb7+gzz09Iv09R8ct92WJ/v44nd2suXJvmlbZ7Vte/bsZ+8rA/Ts2T9h22oUVW81r39g4HDu9+COLU833XtQ5HoBZk/7GjOSZgE3AO8BeoHNkroiYkdZs48BeyPiREkXANcCv11UTWbV+Pa23VyxcTutLS0MDA1x3XkrWHPqsa9pd9GN9/HDnlIAfOmeHt5x4mJuueTsKa2z2rZX3fkwN9+3i0+fMsjl1/+Ai1cexzVrT5lkz4uvt9rX/923DnD5tffkeg+GNct7UOR6hxW5R3Am0BMRT0TEIaATWDuqzVrgG9njO4B3SVKBNZnl0td/kCs2bufVgSH2Hxzk1YEhPrNx+2v+GtvyZN8vQmDYvT19FfcM8q6z2rY9e/aP+AUIcPOPdk35r+Ki6p3M6x+OSPI9KHK95RQR07ayESuWzgdWR8Ql2fSHgLMiYn1Zm0eyNr3Z9L9kbV4Yta51wLps8iRgZyFFT94S4IUJWzWeZu0XTNA3tc6bP3vRMW9RS8us4XkxNHR4cO+zj8fAgVeG581auORXZh159DGjn3/45RefPbz/hWcms85q27bMP2rx7Ne94XiAw6/sY9b8owAYfOn5p4Ze2TfxsaopvgfVtp3M6w/3K897UK5B3oNp+VnM4U0R8YZKCwo7NDSdImIDsKHedYxF0paIaK93HdOtWfsFzds3SVsG9z3nfjWQmfCzWOShod3AsrLppdm8im0kzQaOAiad3mZmVr0ig2AzsFzSCZKOAC4Auka16QI+nD0+H7gnijpWZWZmFRV2aCgiBiWtB+4CZgE3RcSjkq4BtkREF/B/gVsk9QA/pxQWjWjGHraaombtFzRv39yvxlP3vhV2stjMzBqDv1lsZpY4B4GZWeIcBFWS9JSkhyVtk7Qlm3e1pN3ZvG2S3lfvOidD0tGS7pD0z5Iek7RS0uslfVfST7L/F9W7zmqN0a+G32aSTiqrf5uklyR9stG32Tj9aoZtdrmkRyU9IumbkuZmF9Tcnw21c2t2cU1t6/I5gupIegpoL//Sm6Srgf6I+ON61TUdJH0DuDcibsx+GOcDnwV+HhH/KxsvalFEXFHXQqs0Rr8+SRNss2HZkC67gbOAy2jwbTZsVL8+SgNvM0nHAj8ETo6IA5JuAzYB7wO+FRGdkr4KPBQRX6llbd4jMAAkHQX8BqUruYiIQxHxIiOHAfkG8Fv1qG+yxulXs3kX8C8R8VMafJuNUt6vZjAbmJd9b2o+8Czwm5SG2IE6bS8HQfUC+I6krdnQF8PWS9ou6aZG2xXPnAA8D/y5pAcl3SjpSKAtIp7N2vwMaKtbhZMzVr+g8bdZuQuAb2aPG32blSvvFzTwNouI3cAfA7soBcA+YCvwYkQMZs16gekbTS4nB0H1fj0iTgPOBS6T9BvAV4BfBU6ltIG/UL/yJm02cBrwlYh4O/AyMGLo8OzLfo12LHGsfjXDNgMgO9y1Brh99LIG3WZAxX419DbLgmstpT9OfgU4Elhd16IyDoIqZalORDwH/DVwZkTsiYjDETEE/BmlkVcbTS/QGxH3Z9N3UPoFukfSMQDZ/8/Vqb7JqtivJtlmw84FfhwRe7LpRt9mw0b0qwm22buBJyPi+YgYAL4FnAMcnR0qgspD8RTOQVAFSUdKWjj8GHgv8Mjwhy7zAeCRetQ3FRHxM+BpSSdls94F7GDkMCAfBr5dh/Imbax+NcM2K3MhIw+fNPQ2KzOiX02wzXYBZ0uaL0n88jP2PUpD7ECdtpevGqqCpDdT2guA0iGHv4qIz0u6hdLuagBPAR8vO0bbMCSdCtwIHAE8QekqjRbgNuA44KfAByPi5/WqcTLG6NeXaI5tdiSlXzBvjoh92bzFNP42q9Svhv+cSfoflG6+NQg8CFxC6ZxAJ/D6bN5FETH9tyEbry4HgZlZ2nxoyMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ1x83qzvLJLJ/8xm3wjcJjSEBNQ+vLfoboUVoGkDuBQRPy/OpdiiXMQWFOJiD5K15rPiFFhJc0uG0dmtA6gH8gdBBOsz2xSfGjImp6k0yV9Pxso8K6y4Re6JV0vaUt2n4IzJH0rG8f/D7M2x2f3MfjLrM0dkubnWO+fqHS/ik9Ien823vyDku6W1CbpeOA/A5dnY+u/Q9LXJZ1fVnd/9n+HpHsldVH6VvQsSX8kaXM2ANvHa/qGWtNxEFizE/Bl4PyIOB24Cfh82fJDEdEOfJXSV/svA94GfCQ7zARwEvCnEfFW4CXgv0pqnWC9R0REe0R8gdIY9Gdng951Ap+JiKey17w+Ik6NiHsn6MdpwCci4i3Ax4B9EXEGcAZwqaQTqn9rzEp8aMia3RxKv9i/WxrehVmURq4c1pX9/zDw6PCQBZKeAJYBLwJPR8Q/Ze3+Avhd4B8mWO+tZY+XArdmewxHAE9Ooh8PRMTw894LrCjbezgKWD7J9Zo5CKzpidIv+JVjLB8e02Wo7PHw9PDnY/Q4LJFjvS+XPf4y8MWI6MpOEF89xnMGyfbSJbVQCo1K6xPw3yLirjHWY1YVHxqyZncQeIOklQCSWiX9WpXrOG74+cDvUDrUs7OK9R7FL4cW/nDZ/P3AwrLpp4DTs8drgNYx1ncX8F+yw1NIekvZzXbMquYgsGY3RGmI32slPQRsA/5tlevYSekmRI8Biyjd5OZQFeu9Grhd0lbghbL5fwN8YPhkMaUx9t+ZrW8lI/cCyt1IafjiH0t6BPga3ru3KfDoo2bjyK7u+duIeFu9azErivcIzMwS5z0CM7PEeY/AzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxx/x8am+M/AL5ZigAAAABJRU5ErkJggg==\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", + "id": "focal-screening", + "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, + "id": "tribal-sleep", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + ":7: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " family=sm.families.Binomial(sm.families.links.logit)).fit()\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: Mon, 15 Mar 2021 Deviance: 3.0144
Time: 14:12:16 Pearson chi2: 5.00
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
" + ], + "text/plain": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: Frequency No. Observations: 23\n", + "Model: GLM Df Residuals: 21\n", + "Model Family: Binomial Df Model: 1\n", + "Link Function: logit Scale: 1.0000\n", + "Method: IRLS Log-Likelihood: -3.9210\n", + "Date: Mon, 15 Mar 2021 Deviance: 3.0144\n", + "Time: 14:12:16 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", + "id": "tested-jungle", + "metadata": {}, + "source": [ + " The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0850$ 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, + "id": "administrative-topic", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + ":2: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", + "Use an instance of a link class instead.\n", + " family=sm.families.Binomial(sm.families.links.logit),\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: Mon, 15 Mar 2021 Deviance: 18.086
Time: 14:12:16 Pearson chi2: 30.0
No. Iterations: 6
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" + ], + "text/plain": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: Frequency No. Observations: 23\n", + "Model: GLM Df Residuals: 21\n", + "Model Family: Binomial Df Model: 1\n", + "Link Function: logit Scale: 1.0000\n", + "Method: IRLS Log-Likelihood: -23.526\n", + "Date: Mon, 15 Mar 2021 Deviance: 18.086\n", + "Time: 14:12:16 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", + "id": "martial-colleague", + "metadata": {}, + "source": [ + "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. 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.*\n", + "\n", + "## Predicting failure probability\n", + "---\n", + "\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, + "id": "becoming-plaintiff", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZpklEQVR4nO3df5RV5X3v8ffXAWT4IRhUKg4KbRBrUX7M8CtYMyRRMEkRb6lILIlZIaT3hkRrpUtW02itrnVzx1vttdZKlZrWqwOy7ASzWBkSy9TUVh0IIAIdQDORGU1RDD/GDDIM3/6x9xkPw8ycM2fOmXPOw+e11qw5e59n7/18ZzMf9jxnn+eYuyMiIsXvnHx3QEREskOBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiJSBbmZrzOygmb3RzfNmZv/PzPab2etmNi373RQRkVTSuUJ/Cpjfw/M3ABPir+XAY33vloiI9FbKQHf3l4APemhyI/CPHnkFGGlmF2ergyIikp4BWdjHJcCBpOWmeN27nRua2XKiq3hKS0vLx44d2+uDfXDc+ajdscz6WnAcVEuBCaUOUC2FauA5MKo0s5cw9+7d+767X9jVc9kI9LS5+2pgNUBFRYVv2bIlo/3U1dVRWVmZxZ7lj2opPKHUAaqlUPWlFjP7RXfPZeMul2Yg+VK7LF4nIiL9KBuBvgH4cny3yyzgiLufMdwiIiK5lXLIxcyeBSqBC8ysCbgHGAjg7n8HbAQ+D+wHfg18NVedFRGR7qUMdHdfkuJ5B76ZtR6JSNFoa2ujqamJ48eP5/xYI0aMYM+ePTk/Tn9Ip5bBgwdTVlbGwIED095vv74oKiJhaWpqYvjw4YwbNw6z3N6DcuzYMYYPH57TY/SXVLW4O4cOHaKpqYnx48envV+99V9EMnb8+HFGjRqV8zA/25gZo0aN6vVfPgp0EekThXluZPJzVaCLiARCY+giUtRKSkq46qqrOpZramoYN25c/jqURwp0ESlqpaWlbN++vcvn3B1355xzzo7BiLOjShE5azQ2NjJx4kS+/OUvM2nSJA4cOEBVVRXTp0/n6quv5p577ulo+8ADD3D55ZdzzTXXsGTJEh588EEAKisrSUxN8v7773dc8be3t7Ny5cqOfT3++OPAx2/lX7RoEVdccQW33nor0R3dUF9fz6c+9SkmT57MjBkzOHbsGPPnzz/tP6FrrrmGHTt29Ll2XaGLSFb8xQu72P3O0azu88ox53HP7/1Oj21aW1uZMmUKAOPHj+ehhx5i3759fP/732fWrFls2rSJffv28dprr+HuLFiwgJdeeomhQ4dSXV3N9u3bOXnyJNOmTaO8vLzHYz355JOMGDGC+vp6PvroI+bMmcP1118PwLZt29i1axdjxoxhzpw5vPzyy8yYMYPFixezdu1apk+fztGjRyktLWXp0qU89dRTPPzww+zdu5fjx48zefLkPv+8FOgiUtQ6D7k0NjZy2WWXMWvWLAA2bdrEpk2bmDp1KgAtLS3s27ePY8eOcdNNNzFkyBAAFixYkPJYmzZt4vXXX2f9+vUAHDlyhH379jFo0CBmzJhBWVkZAFOmTKGxsZERI0Zw8cUXM336dADOO+88AG666SbmzJlDVVUVa9as4bbbbsvKz0KBLiJZkepKuj8NHTq047G7s2rVKr7xjW+c1ubhhx/udvsBAwZw6tQpgNPuBXd3HnnkEebNm3da+7q6Os4999yO5ZKSEk6ePNnt/ocMGcJ1113HD37wA9atW8fWrVvTqisVjaGLSNDmzZvHmjVraGlpAaC5uZmDBw9y7bXXUlNTQ2trK8eOHeOFF17o2GbcuHEdIZu4Gk/s67HHHqOtrQ2AvXv38uGHH3Z77IkTJ/Luu+9SX18PRO8QTQT9smXL+Pa3v8306dM5//zzs1KrrtBFJGjXX389e/bsYfbs2QAMGzaMp59+mmnTprF48WImT57MRRdd1DEsAnDXXXdx8803s3r1ar7whS90rF+2bBmNjY1MmzYNd+fCCy+kpqam22MPGjSItWvX8q1vfYvW1lZKS0v5yU9+AkB5eTnnnXceX/1qFuczTNzW099f5eXlnqnNmzdnvG2hUS2FJ5Q63HNfy+7du3O6/2RHjx7N6f7vuecer6qqyukxEo4ePerNzc0+YcIEb29v77ZdVz9fYIt3k6sachER6WfPPPMMM2fO5IEHHsjqPfIachERAe69995+O9aXvvSlM16kzQZdoYtIn3j8BhrJrkx+rgp0EcnY4MGDOXTokEI9yzyeD33w4MG92k5DLiKSsbKyMpqamnjvvfdyfqzjx4/3OuAKVTq1JD6xqDcU6CKSsYEDB/bqE3X6oq6uruPdnsUuV7VoyEVEJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQlEWoFuZvPNrMHM9pvZ3V08f6mZbTazbWb2upl9PvtdFRGRnqQMdDMrAR4FbgCuBJaY2ZWdmn0HWOfuU4FbgL/NdkdFRKRn6VyhzwD2u/tb7n4CqAZu7NTGgfPixyOAd7LXRRERSYel+rRuM1sEzHf3ZfHyUmCmu69IanMxsAk4HxgKfM7dt3axr+XAcoDRo0eXV1dXZ9TplpYWhg0bltG2hUa1FJ5Q6gDVUqj6UsvcuXO3untFl0+6e49fwCLgiaTlpcDfdGpzJ/An8ePZwG7gnJ72W15e7pnavHlzxtsWGtVSeEKpw121FKq+1AJs8W5yNZ0hl2ZgbNJyWbwu2deAdfF/EP8BDAYuSGPfIiKSJekEej0wwczGm9kgohc9N3Rq8zbwWQAz+22iQH8vmx0VEZGepQx0dz8JrABqgT1Ed7PsMrP7zGxB3OxPgK+b2Q7gWeC2+E8DERHpJwPSaeTuG4GNndZ9N+nxbmBOdrsmIiK9oXeKiogEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhKItALdzOabWYOZ7Tezu7tpc7OZ7TazXWb2THa7KSIiqQxI1cDMSoBHgeuAJqDezDa4++6kNhOAVcAcd/+VmV2Uqw6LiEjX0rlCnwHsd/e33P0EUA3c2KnN14FH3f1XAO5+MLvdFBGRVMzde25gtgiY7+7L4uWlwEx3X5HUpgbYC8wBSoB73f1HXexrObAcYPTo0eXV1dUZdbqlpYVhw4ZltG2hUS2FJ5Q6QLUUqr7UMnfu3K3uXtHVcymHXNI0AJgAVAJlwEtmdpW7H05u5O6rgdUAFRUVXllZmdHB6urqyHTbQqNaCk8odYBqKVS5qiWdIZdmYGzSclm8LlkTsMHd29z950RX6xOy00UREUlHOoFeD0wws/FmNgi4BdjQqU0N0dU5ZnYBcDnwVva6KSIiqaQMdHc/CawAaoE9wDp332Vm95nZgrhZLXDIzHYDm4GV7n4oV50WEZEzpTWG7u4bgY2d1n036bEDd8ZfIiKSB3qnqIhIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiAQirUA3s/lm1mBm+83s7h7a/b6ZuZlVZK+LIiKSjpSBbmYlwKPADcCVwBIzu7KLdsOB24FXs91JERFJLZ0r9BnAfnd/y91PANXAjV20+0vge8DxLPZPRETSZO7ecwOzRcB8d18WLy8FZrr7iqQ204A/c/ffN7M64C5339LFvpYDywFGjx5dXl1dnVGnW1paGDZsWEbbFhrVUnhCqQNUS6HqSy1z587d6u5dDmsP6FOvADM7B/gr4LZUbd19NbAaoKKiwisrKzM6Zl1dHZluW2hUS+EJpQ5QLYUqV7WkM+TSDIxNWi6L1yUMByYBdWbWCMwCNuiFURGR/pVOoNcDE8xsvJkNAm4BNiSedPcj7n6Bu49z93HAK8CCroZcREQkd1IGurufBFYAtcAeYJ277zKz+8xsQa47KCIi6UlrDN3dNwIbO637bjdtK/veLRER6S29U1REJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBB9fqeoSF/UbGumqraBdw63MmZkKSvnTWTh1Evy3S1Jk85fYVGgS97UbGtm1fM7aW1rB6D5cCurnt8JoFAoAjp/hUdDLpI3VbUNHWGQ0NrWTlVtQ556JL2h81d4FOiSN+8cbu3VeiksOn+FR4EueTNmZGmv1kth0fkrPAp0yZuV8yZSOrDktHWlA0tYOW9innokvaHzV3j0oqjkTeKFM90lUZx0/gqPAl3yauHUSxQARUznr7BoyEVEJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBBpBbqZzTezBjPbb2Z3d/H8nWa228xeN7MXzeyy7HdVRER6kjLQzawEeBS4AbgSWGJmV3Zqtg2ocPergfXA/8l2R0VEpGfpXKHPAPa7+1vufgKoBm5MbuDum9391/HiK0BZdrspIiKpmLv33MBsETDf3ZfFy0uBme6+opv2fwP80t3v7+K55cBygNGjR5dXV1dn1OmWlhaGDRuW0baFRrUUnlDqANVSqPpSy9y5c7e6e0VXzw3oU686MbM/BCqAT3f1vLuvBlYDVFRUeGVlZUbHqaurI9NtC41qKTyh1AGqpVDlqpZ0Ar0ZGJu0XBavO42ZfQ74M+DT7v5RdronIiLpSmcMvR6YYGbjzWwQcAuwIbmBmU0FHgcWuPvB7HdTRERSSRno7n4SWAHUAnuAde6+y8zuM7MFcbMqYBjwnJltN7MN3exORERyJK0xdHffCGzstO67SY8/l+V+iWSkZlszVbUNvHO4lTEjS1k5byLAGesWTr2kX46di+Ok4zs1O3n21QPcMamNr63ayJKZY7l/4VV56Yv0n6y+KCqSTzXbmln1/E5a29oBaD7cysrndoBBW7t3rFv1/E6ArIZtV8fOxXHS8Z2anTz9ytsdy+3uHcsK9bDprf8SjKraho5ATWg75R1hntDa1k5VbUPOj52L46Tj2VcP9Gq9hEOBLsF453BrTtr2ZX/ZPk462rt5b0l36yUcCnQJxpiRpTlp25f9Zfs46Sgx69V6CYcCXYKxct5ESgeWnLZu4DnGwJLTg6x0YEnHi6W5PHYujpOOJTPH9mq9hEMvikowEi8+5uMul+6OnY+7XBIvfCbGzEvMdJfLWUKBLkFZOPWSLkO0P4K1u2Pnw/0Lr+L+hVdRV1fHm7dW5rs70k805CIiEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARiQDqNzGw+8NdACfCEu//vTs+fC/wjUA4cAha7e2N2uyoSrpptzVTVNvDO4VbGjCxl5byJPLflbV5+84OONnN+6xP8QcWlZ7QDzli35Rcf8OyrB7hjUhtfW7WRJTPHcv/Cq9I6blf7Wzj1krT7nTh2uzslZjk5dlfbdtfHs0nKQDezEuBR4DqgCag3sw3uvjup2deAX7n7J83sFuB7wOJcdFgkNDXbmln1/E5a29oBaD7cyh1rt5/R7uU3Pzgt4JsPt7Jy/Q5waDvlHevuXLudU0nbtbvz9CtvA5wWrF0dd+VzO8Cgrf3j/a16fifAGYHZ1fb9ceyutu2uj2ebdIZcZgD73f0tdz8BVAM3dmpzI/D9+PF64LNmZtnrpki4qmobOsKpt9ravSPME0510/bZVw+kPG7bKe8I1ITWtnaqahvO2F9X2/fHsbvatrs+nm3M3XtuYLYImO/uy+LlpcBMd1+R1OaNuE1TvPxm3Ob9TvtaDiyPFycCmZ6BC4D3U7YqDqql8PRrHYN+45Pludp3+6+PUDJkRMfyiV/u35rpcZO37ev2GW57AfB+T9t27mMB68u/scvc/cKunkhrDD1b3H01sLqv+zGzLe5ekYUu5Z1qKTyh1AFRLSePHAymlpDOSy5qSWfIpRkYm7RcFq/rso2ZDQBGEL04KiIi/SSdQK8HJpjZeDMbBNwCbOjUZgPwlfjxIuBfPNVYjoiIZFXKIRd3P2lmK4BaotsW17j7LjO7D9ji7huAJ4F/MrP9wAdEoZ9LfR62KSCqpfCEUgeolkKVk1pSvigqIiLFQe8UFREJhAJdRCQQBR/oZjbYzF4zsx1mtsvM/iJeP97MXjWz/Wa2Nn7BtuCZWYmZbTOzH8bLxVpHo5ntNLPtZrYlXvcJM/uxme2Lv5+f736mw8xGmtl6M/tPM9tjZrOLsRYzmxifj8TXUTO7o0hr+eP49/0NM3s2zoFi/V25Pa5jl5ndEa/LyTkp+EAHPgI+4+6TgSnAfDObRTS9wEPu/kngV0TTDxSD24E9ScvFWgfAXHefknQ/7d3Ai+4+AXgxXi4Gfw38yN2vACYTnZ+iq8XdG+LzMYVoXqVfA/9MkdViZpcA3wYq3H0S0c0YiSlFiup3xcwmAV8nesf9ZOCLZvZJcnVO3L1ovoAhwM+AmUTvshoQr58N1Oa7f2n0vyw+eZ8BfghYMdYR97URuKDTugbg4vjxxUBDvvuZRh0jgJ8T3yBQzLV06v/1wMvFWAtwCXAA+ATRnXg/BOYV4+8K8AfAk0nLfw78aa7OSTFcoSeGKbYDB4EfA28Ch939ZNykiegfQaF7mOhkJqa8GEVx1gHgwCYz2xpP6QAw2t3fjR//Ehidn671ynjgPeAf4qGwJ8xsKMVZS7JbgGfjx0VVi7s3Aw8CbwPvAkeArRTn78obwO+a2SgzGwJ8nuhNmDk5J0UR6O7e7tGfkWVEf7pckd8e9Z6ZfRE46O7FMtdEKte4+zTgBuCbZnZt8pMeXXoUwz2xA4BpwGPuPhX4kE5//hZRLQDEY8sLgOc6P1cMtcTjyTcS/Wc7BhgKzM9rpzLk7nuIhoo2AT8CtgPtndpk7ZwURaAnuPthYDPRn1sj42kGoOvpCArNHGCBmTUSzVj5GaKx22KrA+i4isLdDxKN084A/svMLgaIvx/MXw/T1gQ0ufur8fJ6ooAvxloSbgB+5u7/FS8XWy2fA37u7u+5exvwPNHvT7H+rjzp7uXufi3R2P9ecnROCj7QzexCMxsZPy4lmpd9D1GwL4qbfQX4QV46mCZ3X+XuZe4+jujP4X9x91spsjoAzGyomQ1PPCYar32D06eAKIpa3P2XwAEzmxiv+iywmyKsJckSPh5ugeKr5W1glpkNMTPj43NSdL8rAGZ2Ufz9UuB/AM+Qo3NS8O8UNbOrieZaLyH6D2idu99nZr9JdKX7CWAb8Ifu/lH+epo+M6sE7nL3LxZjHXGf/zleHAA84+4PmNkoYB1wKfAL4GZ3/6Cb3RQMM5sCPAEMAt4Cvkr8b43iq2UoUSD+prsfidcV3XmJb09eDJwk+r1YRjRmXlS/KwBm9lOi18vagDvd/cVcnZOCD3QREUlPwQ+5iIhIehToIiKBUKCLiARCgS4iEggFuohIIPr1Q6JF0hXf1vVivPgbRO+uey9enuHuJ/LSsS7Et6GecPd/z3NX5CynQJeC5O6HiGbXxMzuBVrc/cF89cfMBiTNI9JZJdACpB3oKfYnkhENuUjRMLNyM/vXeEKw2qS3TteZ2UNmtiWez3y6mT0fzzV9f9xmXDzf+f+P26yPJ0tKtd+HLZrv/XYz+714Pu5tZvYTMxttZuOAPwL+OJ6D/HfN7CkzW5TU75b4e6WZ/dTMNgC740nnqsys3sxeN7Nv9OsPVIKjQJdiYcAjwCJ3LwfWAA8kPX/Co3nZ/47obdTfBCYBt8XDNwATgb91998GjgL/y8wGptjvIHevcPf/C/wbMCuexKsa+FN3b4yP+ZBHc5H/NEUd04Db3f1yovm8j7j7dGA68HUzG9/7H41IREMuUizOJQroH0fTe1BCNLVqwob4+05gV2JqUjN7i2i60sPAAXd/OW73NNGHKPwoxX7XJj0uA9bGV/CDiOZR763X3D2x3fXA1UlX8yOACRnuV0SBLkXDiIJ6djfPJ+b0OJX0OLGc+HfeeZ4LT2O/HyY9fgT4K3ffEL8Qem8325wk/uvXzM4hCv+u9mfAt9y9tpv9iPSKhlykWHwEXGhmswHMbKCZ/U4v93FpYnvgS0RDKA292O8IPp6y9StJ648Bw5OWG4k+Ag6ieckHdrO/WuB/xsM+mNnl8eRaIhlRoEuxOEU0der3zGwH0QcFfKqX+2gg+jCOPcD5RB9qcaIX+70XeM7MthJ9HFrCC8BNiRdFgb8HPh3vbzanX5Une4JoWtifmdkbwOPor2bpA822KGeF+G6UH3r0ocMiQdIVuohIIHSFLiISCF2hi4gEQoEuIhIIBbqISCAU6CIigVCgi4gE4r8Bc2QgqRDpwXwAAAAASUVORK5CYII=\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", + "id": "breathing-junction", + "metadata": {}, + "source": [ + "This figure is very similar to the Figure 4 of Dalal et al. *I haven't managed to replicate the Figure 4 of the Dalal et al. article. The curve is missing*.\n", + "\n", + "## Computing and plotting uncertainty\n", + "---\n", + "Following the documentation of Seaborn, I use regplot." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "loaded-collapse", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxjUlEQVR4nO3deXgUZZ4H8G9V9Zn7vrhBwcipgHEUZQxIMhIQdRGXWd31gFlFUWfcR0RHDhUH3VURr1FnxnFgdZf1gqiAqCggcigSIBwaAiGkc3Xuo6+qd/9o0iYEpFPk6O58P8/DA92pqv69dJJvV71vva8khBAgIiLqILmnCyAiouDEACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLSpVsCZPny5cjMzMSwYcNw5MiRM26jqiqWLFmCyZMn49prr8WaNWu6ozQiItKpWwJk0qRJWL16Nfr06XPWbdatW4eioiJs3LgR//M//4OVK1eiuLi4O8ojIiIduiVAxo0bh9TU1F/c5pNPPsHMmTMhyzLi4uIwefJkrF+/vjvKIyIiHQKmD8RmsyEtLc33ODU1FaWlpT1YERER/ZKACRAiIgouhp4uoEVqaipKSkowatQoAO3PSPxVXd0ITQvN6b3i4yNgtzf0dBldJpTbF8ptA9i+YCbLEmJjw3XtGzABkp2djTVr1mDKlCmoqanBpk2bsHr16g4fR9NEyAYIgJBuGxDa7QvltgFsX2/ULZewnnzySVx99dUoLS3F7bffjqlTpwIA5syZg3379gEArr/+evTt2xdTpkzBzTffjHnz5qFfv37dUR4REekghdp07nZ7Q8h+UkhMjERFRX1Pl9FlQrl9odw2gO0LZrIsIT4+Qt++nVwLERH1EgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6GLrrhQoLC7FgwQLU1NQgJiYGy5cvx8CBA9tsY7fb8cgjj8Bms8Hj8SAjIwOPPfYYDIZuK5OIiPzUbWcgixYtwuzZs7FhwwbMnj0bjz/+eLttXnvtNQwZMgTr1q3D2rVrceDAAWzcuLG7SiQiog7olgCx2+3Iz89HTk4OACAnJwf5+fmoqqpqs50kSWhsbISmaXC5XHC73UhOTu6OEomIqIO65dqQzWZDcnIyFEUBACiKgqSkJNhsNsTFxfm2u+eee3DfffdhwoQJaG5uxm9/+1uMHTu2Q68VHx/RqbUHmsTEyJ4uoUuFcvtCuW0A29cbBVTnwvr16zFs2DD8/e9/R2NjI+bMmYP169cjOzvb72PY7Q3QNNGFVfacxMRIVFTU93QZXSaU2xfKbQPYvmAmy5LuD97dcgkrNTUVZWVlUFUVAKCqKsrLy5Gamtpmu1WrVmH69OmQZRmRkZHIzMzEjh07uqNEIiLqoG4JkPj4eKSnpyM3NxcAkJubi/T09DaXrwCgb9+++PrrrwEALpcL27dvx4UXXtgdJRIRUQd12yisxYsXY9WqVcjKysKqVauwZMkSAMCcOXOwb98+AMDChQvx3XffYdq0aZgxYwYGDhyIm2++ubtKJCKiDpCEECHVYcA+kOAVyu0L5bYBbF8wC/g+ECIiCj0MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0oUBQkREujBAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHTxO0A2bdoEj8fTlbUQEVEQ8TtAXnzxRUyYMAFLly7F3r17u7ImIiIKAn4HyNq1a/HWW2/BbDbjvvvuQ1ZWFl555RUUFxf7tX9hYSFmzZqFrKwszJo1C8eOHTvjdp988gmmTZuGnJwcTJs2DZWVlf6WSERE3UgSQoiO7iSEwPbt2/GnP/0JP/74Iy699FLMmjULOTk5kOUzZ9Jtt92Gm266Cddffz0++ugjvPfee3j77bfbbLNv3z48/PDD+Pvf/47ExETU19fDZDLBbDb7XZvd3gBN63CTgkJiYiQqKup7uowuE8rtC+W2AWxfMJNlCfHxEfr27egORUVFePnll7F48WI4nU7Mnz8fM2fOxOrVqzF//vwz7mO325Gfn4+cnBwAQE5ODvLz81FVVdVmu7feegt33HEHEhMTAQCRkZEdCg8iIuo+Bn83XL16NT766CMcP34cv/nNb/DMM89gzJgxvq9nZWXhiiuuOOO+NpsNycnJUBQFAKAoCpKSkmCz2RAXF+fbrqCgAH379sVvf/tbNDU14dprr8Xdd98NSZJ0No+IiLqK3wHy9ddf4/bbb8ekSZNgMpnafd1qtWLlypXnVYyqqjh8+DD+9re/weVy4a677kJaWhpmzJjh9zH0nooFi8TEyJ4uoUuFcvtCuW0A29cb+R0gL774ImRZhtFo9D3ndrshhPAFyoQJE864b2pqKsrKyqCqKhRFgaqqKC8vR2pqapvt0tLSkJ2dDZPJBJPJhEmTJiEvL69DAcI+kOAVyu0L5bYBbF8w65Y+kDvuuAMHDhxo89yBAwdw5513nnPf+Ph4pKenIzc3FwCQm5uL9PT0NpevAG/fyNatWyGEgNvtxrfffouLLrrI3xKJiKgb+R0ghw8fxujRo9s8N2rUKBw6dMiv/RcvXoxVq1YhKysLq1atwpIlSwAAc+bMwb59+wAAU6dORXx8PK677jrMmDEDF1xwAf7pn/7J3xKJiKgb+X0JKyoqCpWVlb4RUgBQWVkJq9Xq1/5DhgzBmjVr2j3/xhtv+P4tyzIeeeQRPPLII/6WRUREPcTvM5ApU6bgD3/4A44cOYLm5mYcPnwYDz/8MH7zm990ZX1ERBSg/A6QBx98EEOGDMHMmTN9Nw4OGjQIv//977uyPiIiClAdvhNdCIHq6mrExsYG5P0ZHIUVvEK5faHcNoDtC2bnMwrL7z4QAKivr0dhYSEaGxvbPP+rX/1K14sTEVHw8jtA3n//fSxduhRhYWGwWCy+5yVJwueff94lxRERUeDyO0Cef/55rFixAhMnTuzKeoiIKEj43YmuqupZ7zQnIqLex+8AmTNnDl599VVomtaV9RARUZDw+xLWW2+9hcrKSrz55puIiYlp87XNmzd3cllERBTo/A6QZ599tivrICKiION3gFx22WVdWQcREQUZv/tAXC4Xnn/+eUyaNAljx44FAGzduhWrVq3qsuKIiChw+R0gy5Ytw5EjR/Cf//mfvjvQL7zwQrzzzjtdVhwREQUuvy9hbdq0CRs3bkRYWBhk2Zs7ycnJKCsr67LiiIgocPl9BmI0GqGqapvnqqqq2o3IIiKi3sHvAMnOzsbDDz+MEydOAADKy8uxdOlSTJ06tcuKIyKiwNWh6dz79u2L6dOno66uDllZWUhKSsK8efO6sj4iIgpQHZ7OHfBeuuJ07t0vlKeUBkK7faHcNoDtC2bdMp17y6WrFq2ndO/Xr5+uFyciouDld4Bce+21kCQJrU9YWs5ADh482PmVERFRQPM7QA4dOtTmcUVFBV566SWMGzeu04siIqLA53cn+ukSExPx6KOP4rnnnuvMeoiIKEjoDhAAOHr0KJqbmzurFiIiCiJ+X8KaPXt2m1FXzc3N+OmnnziMl4iol/I7QGbOnNnmsdVqxUUXXYSBAwd2dk1ERBQE/A6QG264oSvrICKiION3gKxYscKv7e6//37dxRARUfDwO0COHz+OjRs3YsSIEejTpw9KSkqwb98+TJkyBWazuStrJCKiAOR3gAgh8F//9V/IysryPbdx40asX78eTz/9dJcUR0REgcvvYbxff/01Jk+e3Oa5zMxMfPXVV51eFBERBT6/A2TAgAFYvXp1m+feeecd9O/fv9OLIiKiwOf3Jawnn3wS9957L958803fSoQGgwErV67syvqIiChA+R0gF198MTZs2IC9e/eivLwciYmJGDNmDIxGY1fWR0REAUr3VCbjx4+H2+1GU1NTZ9ZDRERBwu8zkMOHD+Puu++GyWRCWVkZrrvuOuzatQsffPABXnjhhS4skagjBIDAW+iMKBT5fQayePFizJ8/H+vXr4fB4M2d8ePH47vvvvNr/8LCQsyaNQtZWVmYNWsWjh07dtZtjx49itGjR2P58uX+lkcEIQCXJzRXoyQKRH4HyE8//YTrr78ewM8LSYWFhcHpdPq1/6JFizB79mxs2LABs2fPxuOPP37G7VRVxaJFi9oNGSbyh9Ot9nQJRL2G3wHSp08f7N+/v81zeXl5fg3jtdvtyM/PR05ODgAgJycH+fn5qKqqarft66+/jl//+tecpJF0cblVaIJnIUTdwe8+kPvvvx+/+93vcMstt8DtduPPf/4z3n33XTzxxBPn3NdmsyE5ORmKogAAFEVBUlISbDYb4uLifNsdOnQIW7duxdtvv41XXnlFR3Oge3H4YJGYGNnTJXSp82mfqgl4ZAlhViMiwwNveh2+d8Et1Nunh98Bcs011+DNN9/E//7v/2L8+PE4efIkVq5ciREjRnRKIW63G3/84x/x9NNP+4JGD7u9AZoWmp9AExMjUVFR39NldJnzbZ8QQHVtM2prgfgoa0B1pfO9C26h3D5ZlnR/8PYrQFRVRVZWFj755BMsXry4wy+SmpqKsrIyqKoKRVGgqirKy8uRmprq26aiogJFRUWYO3cuAKCurg5CCDQ0NPh1lkPUwuMRcLlVmI36P4gQ0bn5FSCKokBRFDidTphMpg6/SHx8PNLT05Gbm4vrr78eubm5SE9Pb3P5Ki0tDTt27PA9XrlyJZqamvDwww93+PWImhweWEwK2B1C1HX87kS/7bbb8MADD2Dnzp0oKirCiRMnfH/8sXjxYqxatQpZWVlYtWoVlixZAgCYM2cO9u3bp696orNweVQ43VpPl0EU0iQhfvkzWkVFBRITE3HRRRd5d5AktN5FkiQcPHiwa6vsAPaBBK/O6AOpqG32vf9WswExEaaAOAvhexfcQrl959MHcs4zkJb1Pw4dOoRDhw4hMzPT9+9Dhw4FVHgQteZweeBWeRZC1FXOGSCnn6Ds2rWry4oh6kxCePtCiKhrnDNAWu46b3GOK15EAaXZ6YEnRC9pEvW0c47CUlUV3377rS84Tn8MAL/61a+6rkKi8yAE0NjsDpi+EKJQcs4AiY+Px8KFC32PY2Ji2jyWJAmff/5511RH1AkcLg9cHgOMiu7VC4joDM4ZIF988UV31EHUZYQAGprciI0MvOlNiIKZ31OZEAWqvIJKrN9RBIdLRbjFgHEXJWFY/9g22zjdKppdHlhN/JYPZi3vdWWtAwnRFmRn9MeoIQk9XVavxXN6Cmp5BZVY/dkR1DS6YLUoqGt2Y+22Qhwuqm63bWOTmzP1BrHW73WYxYCaRhdWf3YEeQWVPV1ar8UAoaC2fkcRFEWG2ahAkiSYDAoURcaWvSXttvVoAg3N7h6okjrD6e+12eh9r9fvKOrp0notBggFtcpaB0yGtt/GRkVGdf2ZFzprdnrg8vDmwmB0pvfaZJBRWevooYqIAUJBLSHa0i4Q3Kp21g5zIYCmZjekQJrrnfxypvfa5dGQEG3poYqIAUJBLTujP1RVg9OtQggBl0eFqmq4anTaWfdxcqLFoHT6e+10e9/r7Ixzr4pKXYNDUiiotYzAWb+jCM0OFVFWIzIv6dNuFFZrLTcXmiJNQEAtO0W/pPV7zVFYgYEBQkFv1JAEjBqS0G423l/i8qhwuQVMRgZIMGl5rykw8BIW9UpCAA0OjsgiOh8MEOq1XG4VDrfa02UQBS0GCPVqjbwvhEg3Bgj1ah6PBhdHZBHpwgChXk0AaHLwvhAiPRgg1OvxvhAifRgg1OsJATQ0u8FpFok6hgFCBO+IrGYX108n6ggGCNEpDY1uuFVeyiLyFwOEQoK91oFNu0+c13TtmhCobXByzRAiP3EqEwoJa7cVYkueDbGRZvxr9kVIirXqOo5HFahvciM63AjOk0X0y3gGQiHhypGpMBm964C89tF+/FRcq/tYDqcHDo7KIjonBgiFhKH9YvDw7EsRGWaEw6XirU8PYkd+ma5jCQD1jS6ofkzKSNSbMUAoZAxMicK8G0ciNT4MmgA+2lqIdd8c0xUEqiZQ3+wCOLiX6KwYIBRSYiLMmDt9ONIHeNcD2b6/FH//9BCanR0foutwqmhycbJForNhgFDIMRsV/Pbaobj61KqEP52sxSsf7EdZVVOHj9XQ5IamsT+E6EwYIBSSZFlCdkZ/3HzNBTAoEux1Drz60X7sL6zq0HE0TaCuiTP2Ep0JA4RC2pgLEzB3+nBEh5vgcmv478+OYMPOIr9WLWzhcKmob+aEi0SnY4BQyOubGIF5N47EoNQoAMBXP5Tgb58e7NBNh00ONxp19KMQhTIGCPUKEVYj7piajgkjUwEABSfr8NL7+3C8tN6v/YXwDu1tcnp4JkJ0SrcFSGFhIWbNmoWsrCzMmjULx44da7fNyy+/jKlTp2LatGm48cYbsWXLlu4qj3oBRZZw3a8G4J8nXwiTUUZdowtvrMvH13tL/Jq+RAigvsnFPhGiU7otQBYtWoTZs2djw4YNmD17Nh5//PF224waNQr/93//h3Xr1mHZsmV48MEH4XA4uqtE6iVGDo7HvBtGIjnWCk0IrN9RhH+sP+zXJS0hvMvgVtU7OTqLer1uCRC73Y78/Hzk5OQAAHJycpCfn4+qqrYjYq666ipYrd45jIYNGwYhBGpqarqjROplEmOsuPuGERg3LBEAcPhEDV56Lw8FJf5NgeJyq6iqc8LlYYhQ79UtAWKz2ZCcnAxFUQAAiqIgKSkJNpvtrPt8+OGH6N+/P1JSUrqjROqFTAYFN04cgpszL/Be0mpy46+5B7FhZxE8fkzr7tEEqusd7FynXisgZ+PduXMnVqxYgb/+9a8d3jc+PqILKgociYmRPV1Clzqf9qmagEeW0NHZ2DMvC8eICxPxl7UHcNxWh69+KEGhrR53TB+OlPhwv44hDDKiI8wwG5WzbsP3LriFevv0kITo+sUP7HY7srKysGPHDiiKAlVVkZGRgY0bNyIuLq7Ntnv27MEDDzyAV155BcOHD9fxWg0dGuMfTBITI1FR4d+ooWB0vu0TAqiobdb9/quahk27i/H1DyUQAAyK92bEy4enQPZj6JUsS4gKN8FqUtqFGN+74BbK7ZNlSfcH7265hBUfH4/09HTk5uYCAHJzc5Gent4uPPLy8vDggw/ixRdf1BUeROdDkWVkXdYfd027GDERJnhUgdxvjuOvHx9Edf25B3NomndBqpoGFzvYqVfoljMQACgoKMCCBQtQV1eHqKgoLF++HIMHD8acOXMwf/58jBw5EjfddBNOnjyJ5ORk337PPPMMhg0b5vfr8AwkePX0GUhrDpcHH39zHN8dqQAAmIwysjP647L0ZL/PRsIsBoRbDJAg8b0LcqHcvvM5A+m2AOkuDJDgFUgB0mLT7hP46ocS35TwybFWXDEiBXt/qkR1vROxkWZcNToNw/rHnnF/gyIhMsyEvmkxutuWV1CJ9TuKUFnrQEK0BdkZ/TFqSILuNnWmtVuPYuOuYjjcKixGBVPG98X0CYN7uqxOF8o/ewF/CYsoGB0uqsaeHysQHWGC1eztHC+rbsYHWwpRWt0Ms0lBXbMba7cV4nBR9RmP4VEFahqcsNc6dK21nldQidWfHUFNowthFgNqGl1Y/dkR5BVUnlfbOsParUex9ptjcLpVGGTA6Vax9ptjWLv1aE+XRt2EAUJ0Flv2lkBRZFhMBsRGWhAXZfZ9rcnhQWWtAxCAosjYsrfkrMcRwntJrKrWgSanB6IDi1St31EERZFhNiqQJAlmowJFkbF+R9F5ta0zbNxVDAkSFFmCJMnevyFh467ini6NuklADuMlCgTV9U5YzD//iFhMBkhw+n79e1SByloHrGYFbve5F57yaAJ1jS40NUuwWgywmg3n7E+prHUgzNL2x9RkkL3h1cMcLg8UuW39suR9nnoHnoEQnUVspBnu024oNCgSjIqEhGgLjIr3x6fZqaK+yY1v80v96n/xaAL1TW5U1jrQ4HD/4qWthGhLu7vdXR4NCdEWHS3qXBaTAac3VxPe56l3YIAQncVVo9OgqhpcHhVCCLg8KkxGBSaTAZCA+Ggzwq0GSPD+4ly79Rhe/mAfjpbU+XV8TRNoaHLDfipIVE1rN9NvdkZ/qKoGp9tbg9OtQlU1ZGf07/wGd9CU8X0hIKBqAkJo3r8hMGV8354ujboJPyoQncWw/rGYDm9fSMuIq6mXDwBaPZccY8W4y5NQcLIO3x2pgM3ehDdz8zF8YByyM/oj3o8zBfVUkDQ1e2A2KbCYFRgVb59Cy2irQByF1TLaqjeMwqIz4zDeIBLKQwmBwBzG2xHF5Q3I3X4MRWUNALzTx2dcnIxrLu2DfmkxqKpq9PtYsizBqMgwmxSYDDIMigQgcBci4fdm8DqfYbw8AyHqJH2TIvC76cORV2DHhp1FqGlw4Zv9pfjucAWyLh+ASy6I/8W5slrTNAGnpsLpViFJgEGREW41wmKUEchBQr0LA4SoE0mShNEXJODigXHYfqAUm/echMOlYu2Wo/h89wn8ekwaLktPhtHgf/ejEIDbo6G23gmHSYHFbIBBlmBQ2IVJPYsBQtQFjAYZV49Ow7hhSfh670lsP1CGxmY3Pt5+HFvybPj1mDSMuyipQyEgADhcKhwu71mJIkswGRVYjApMfp7ZEHUm9oEEkVC+DgsEfx/IL5GMCj744kd8d7jCNy1KVLgJV49OxbiLkmAynF8AyLIEi6nnwoTfm8GLfSBEAS420oIZVw3GxDFp+HJPCb4/XIG6RhdyvzmOL/eU4MoRKci4OBlWs74fSU0TaHJ40OTweMPEqMBkUmBUvJe6QutjIgUKBghRN4qNtODGqwfjmkvS8NUPJfjucAUam93YuMs7aeP4i5JwxcgUxESYz32ws9A0gSanB01OD2RJgqIAJqMBJqMMg+wd0cVAoc7AACHqAS1nJJmX9sW2fTbsOFgGp1vF1n02fLPfhhGD43HlyBT0Szq/VfA0IaB5ALfHjcZmtAkUs0GGokhQZHbGkz4MEKIeFBVuwm8uH4BfX9IHO/LLsH1/Keqb3cgrsCOvwI6+ieH41YgUjBwc3ymjrtoECvBzZ7xBgdEoe29gVGQOFCa/MECIAoDVbMCvL+mDCaNSkVdgx7Z9NtjsTSiuaMSaLwvw8fbjGDcsEePTkxEf1XnzYAnhnRTSo3oApzdQZEmC0SDDZFRgULyXvBQda81T6GOAEAUQgyLj0qGJuOTCBBwrrcf2A6XIL6xCk8ODr/fa8PVeG4b0icK4YUm4eGBch+4n8YcQgCoE1JbhwgAkWYIiA0aDAqNB9oWKLDFUejsGCFEAkiQJg1KjMCg1CnWNLuw6VI7dh8pR2+hCwck6FJysg8WkYPQFCbh0aAL6JkZA8mOp3Y4SAIQmoGmA2+M5VZu3PkUGDLICg0GCpdEFp0f1rQ8iyxJkCQyYEMcAIQpwUeEmTBrbF9dc0gc/Ftdg16FyHDpeA4dLxY78MuzIL0NCtAWjL0jA6AvikRBt7dJ6hACEOBUq8AAuwGhxobrO6Q0XSIAEyDJgkL19KooseTvspZZwkSAxYIIeA4QoSMiyhGH9YzGsfyzqm1zY+5Md3x+pQGlVEyprHfj8u2J8/l0x+iSEY+SQeIwcHIfYyO5dN0QIeFdcFICmAR6oAH5ebKvl7EWWToXLqYkiFUmGLAMKL40FFQYIURCKDDNhwqhUTBiVCpu9ET/8WIm9BXbUNbpwsrIRJysbsX5HEfokhGP4oDhcPDAOSbFde2biD9/ZCwCPqgLuM4TLqTMX+dSlMOnU2UrLkrktZy7e7U8NTW61MiLDp/swQIiCXGp8OFLjw5GV0R/HS+uRV2DHgcIqNDS7fWGycdcJJERbkD4gFhcNiEX/5Mh2y9H2tNaXxjw4+xLBp3f1SJIECd5LZookQ1Yk71Bk75U0byidtqMswRdOLf8NDJ6OY4AQhQi5Vcf7tCsG4lhpPfYX2nHwWDVqG12orHVgS54NW/JssJoVXNg3BkP7xeCCvtGICjP1dPl+O/0Xfct0fqoGuKEC7jPvJ8E7KEA69aDlbEZq6auRT4WPLLVaqlWCJAMutwoB0eYMiBggRCFJliUMTovC4DRvmJysbMTBY9U4VFQNm70JzU7Vd7MiAKTEheGCvtEYkhaFgalRfq9bEkxE679b9dXgDH01rUkSAEVBdY3j5zOdU4MDIqwG76CBXooBQiGl9/4on50kSeibGIG+iRG4dnw/1DQ4cbioBj8W1+Cnk7VwuTWUVjWhtKoJW/NskCUJfZPCMTgtGoNSI9E/KRJmU+gFir9ES85obc90ZElDuMXQ7pJab8IAoZAhSQLREWbvp0tJ+K6pt/wC8P1bE9A0AY/QoGmndm75RIqfH581jUTbSyG+p4PkskZMhBkZFycj4+JkeFQNJ8ob8GNxLQpO1qK4ogGaECgqa0BRWQM27/H2F6QmhGNAciQGpESif3IkosOD55IXdR0GCIUQCSaDfx8HWz41aq1D5lQCeMOm7Xbeo3u3b611P7QmvHNNqaoGVRVwqSo01ftcoDIosq/fBOP7weHyoNBWj6MltSi01cNmb4QmgJMVjThZ0Yhv9pcC8N6b0i8xAv2SItAnKRzWcP2zB1PwYoBQr+QLCHgv8bQ7nTgfp/oPJAnwqBocbg0BNuDprCwmA9IHxCJ9QCwAwOHy4HhpvfdPWT2KyxvhVjXUNbpwoLEKB45VndrzIBKiLUhLCEdafDhSE8KQGh+OCKux5xpDXY4BQtRFhPB2toabZcTGhsHlcHuHqQpAUwXcqgpVEwF96ctiMvhuXgQAVdNQam9CUXkDissbcKK8AfZaBwSAyloHKmsdvo55AIi0GpESH4bkuDAkx1qRHBuGxFhrSHbS90YMEKJuYFBkhJ1htUGPpsHhUqF6hHcSQ1WDJgI3VBRZRp/ECPRJjACGe59zuDyod2o4VFgJW2UTTlY2orK2GUIA9c1u1BfX4sfi2jbHiYkwITHGisQYKxJiLEiItiIh2oKocBPk3twrHWQYIEQ9yCDLiLB47zqQJG8/itujwe1RoQlAVTV4PBrUAA4Vi8mAtJRwJEb+3LHu9mgoq25Cqd07uqusugmlVc1obPbepFHT4EJNg6tdsBgVGXFRZsRFWRAfZUFslBlxkWbERlkQG2Hu9NmH6fwwQIgChBDeXhiTQYap1S9KAe+lI4/6cwe9W9XgUbWADRWjQfYNHW6t0eFGeXUzyqubUVnTjIraZlTUOFBT74QA4FY1lFU3o6y6+YzHjbAaERNhQkykGTHhZkRHmBAdYUZ0uAnR4SZEWI2Qg6XDKQQwQIgCnATvmYpBRpsOerdHg9OjwelS4fFoAT3aq0W4xYhBqUbvqK9W3B4NVfUO2GsdsNd5/66ud6KqzomaBifUU8PfGprdaGh2o7ii8YzHlyRvv0tkuAmRVhMiw4yn/njDJTLMiAirEeFWI0wGuUumwO9NGCBEQUgInFrYSUaExeBdVVDT4HJpcHk8UDXv/S6BHyleRoOM5NgwJMeGtfuapgnUNblQXe9Edb03UGoaXKhtcKK20YXaBhecpyZlFAKoa3KjrskN4Mwh43tNRUaYxYBwqxHhFgPCLUZYLQaEmQ0Ia/W31WyAKstwOT0wmxT20bTCACEKct7RXhIUWYHZoECSjFA1cWqpWg1utwa3qkLTTt3z0tMFd5AsS4iJMCMmwoxBqWfexuHyoK7RjbpGF+qaXKhvcqGu0Y36Zhfqm9xoaPL+2+XWfPu4Vc0bQI0uv2uRAJhNCiwmBeEWI26YOBiXXph4ni0MXt0WIIWFhViwYAFqamoQExOD5cuXY+DAgW22UVUVTz75JLZs2QJJkjB37lzMnDmzu0okCglCeCdWNBkkb1+KGcg/ZsfmH0pQ1+RGXKQZ4y9KwqCUKPzf5h+xt6DKN5x4xMAYjBiSgC17S1Bd70RspBlXjU7DsP6xOFxU3e75kxUN2JpXCqdHhdmgYMKoFGSO7XfGus60/9mO2zJs2N9jfLPP5q3DrcJs9NYx7cqB7fZ3eVTs+8mO7QdKUdvoQpjZgP7JEXC4VRwvrYfTpUKWJZiNCjyagMPpaXPzqADgOLXcb02DC3/7+CA2JZ1AdkZ/jBqS0AnvXnCRhOieC6e33XYbbrrpJlx//fX46KOP8N577+Htt99us82HH36IdevW4Y033kBNTQ1mzJiB//7v/0bfvn39fh27vcE3Z02oSUyMREVFfU+X0WVCuX092ba8gkqs/uwIFMXbOe/2aAAEkqKtOFpaD0WRfSsEyrIEs0FGTKQFsiTgcGtwu1WMGBSHPT9VQpIkGGQJTo+G+gYHml3eT/QShPeyGYBJl/ZpFyKHi6qxdlshFEWGUZHhVjWoqoaxQxPx3ZGKds9Pv3JQuxA52zEGJEVg79Eq70SHp0aydaSOZocHkCRYzUq7Gob2i4HTrcJkMaG0vB7NTg8KTtZi56EymAwKkmKtqGtyQ1U1/PbaoUEZIrIsIT4+4twbnkG3nIHY7Xbk5+fjb3/7GwAgJycHTzzxBKqqqhAXF+fb7pNPPsHMmTMhyzLi4uIwefJkrF+/HnfddZffrxXqIzDYvuDVU23bkV+GpLgwmAw/37zn8qiwVTfDaja0uf9eAHB5NJgM3gWdzGYJHlXDj8W1SI0P994AeGrhp4pqEyTJ23+hSDIEBFRNQ2WtE9HhZggIb8e+AI6V1WFwnxgoivfVhPBOkV5U1oh+yZEwKgoEvGdBblXFoaIajBwS75tCVwDIP1aN5PhT7RAttaqwVTUjOdba5v9X0wQOn6hFVsaANv8XBwqr2v1f2Gu9I77iWy0F7PKoOFBYhRGD42EyKoiJCYP51MC4vT9Vol9yJCKsJljNCsKtHrg83uWFxwTh5azz+b7slgCx2WxITk6GonjfNEVRkJSUBJvN1iZAbDYb0tLSfI9TU1NRWlraodeKjQ3vnKIDlN5PCsEilNvXU21beMflPfK6rT30L5ed9zEeH3L+v5wXnccxUuLDz/sYoYZ35RARkS7dEiCpqakoKyuDqnqH2qmqivLycqSmprbbrqSkxPfYZrMhJSWlO0okIqIO6pYAiY+PR3p6OnJzcwEAubm5SE9Pb3P5CgCys7OxZs0aaJqGqqoqbNq0CVlZWd1RIhERdVC3jcIqKCjAggULUFdXh6ioKCxfvhyDBw/GnDlzMH/+fIwcORKqqmLp0qXYtm0bAGDOnDmYNWtWd5RHREQd1G0BQkREoYWd6EREpAsDhIiIdGGAEBGRLgwQIiLSJWhn473nnntQXFwMWZYRFhaGP/7xj0hPT/dr0sZg8dJLL2HlypVYt24dhg4dih9++AGPP/44nE4n+vTpg2effRbx8fE9XaYumZmZMJlMMJvNAICHHnoIV111VUi00el0YtmyZdi+fTvMZjPGjBmDJ554IiS+N4uLizFv3jzf4/r6ejQ0NGDnzp0h0T4A+PLLL7FixQrvzMVC4N5778WUKVNCpn2bN2/GihUr4PF4EB0djaeffhr9+vXT1z4RpOrq6nz//uyzz8SMGTOEEELceuut4sMPPxRCCPHhhx+KW2+9tUfqO1/79+8Xd955p7jmmmvE4cOHhaqqYvLkyWLXrl1CCCFefvllsWDBgh6uUr+WdrUWKm184oknxFNPPSU0TRNCCFFRUSGECJ3vzdaefPJJsWTJEiFEaLRP0zQxbtw43/fmwYMHxZgxY4SqqiHRvpqaGnHZZZeJo0ePCiG87bjjjjuEEPrev6ANkNY++OADccMNN4jKykoxduxY4fF4hBBCeDweMXbsWGG323u4wo5xOp3i5ptvFidOnPD9ot27d6+YOnWqbxu73S7GjBnTg1WenzMFSCi0saGhQYwdO1Y0NDS0eT5UvjdbczqdIiMjQ+zfvz9k2qdpmrjsssvE7t27hRBC7Ny5U0yZMiVk2rd3715x3XXX+R5XV1eLoUOH6m5f0F7CAoBHH30U27ZtgxACb775pt+TNga6FStWYPr06W2msT99osm4uDhomuY73QxGDz30EIQQGDt2LH7/+9+HRBtPnDiBmJgYvPTSS9ixYwfCw8Nx//33w2KxhMT3ZmtffPEFkpOTMXz4cOzfvz8k2idJEl544QXcc889CAsLQ2NjI15//fWQ+d0yaNAgVFZWIi8vD6NGjcK6desA+D/h7emCuhP9qaeewubNm/Hggw/imWee6elyOsWePXuwf/9+zJ49u6dL6VKrV6/G2rVr8d5770EIgaVLl/Z0SZ1CVVWcOHECF198Md5//3089NBDuO+++9DU1NTTpXW69957DzfddFNPl9GpPB4P/vznP+OVV17Bl19+iVdffRUPPPBAyLx/kZGReP755/H000/jxhtvhN1uR1RUlO72BXWAtJgxYwZ27NiBlJQUvyZtDGS7du1CQUEBJk2ahMzMTJSWluLOO+/E8ePH20w0WVVVBVmWg+aT+ela3hOTyYTZs2fj+++/bzeZZjC2MTU1FQaDATk5OQCA0aNHIzY2FhaLJei/N1srKyvDrl27MG3aNAD+T5ga6A4ePIjy8nKMHTsWADB27FhYrVaYzeaQaB8AXHHFFXjnnXfw/vvv41/+5V/gcDjQp08fXe0LygBpbGyEzWbzPf7iiy8QHR3t96SNgWzu3LnYunUrvvjiC3zxxRdISUnBX/7yF9x1111wOBzYvXs3AODdd99FdnZ2D1erT1NTE+rrvavzCSHwySefID09HSNGjAj6NsbFxSEjI8M3n1thYSHsdjsGDhwY9N+brX3wwQeYOHEiYmO9qwaGws8eAKSkpKC0tBRHjx4F4J3Dz263Y8CAASHRPgCoqKgAAGiahueeew633HIL+vTpo6t9QTkXVmVlJe655x40NzdDlmVER0fj4YcfxvDhw886aWOwyszMxGuvvYahQ4fi+++/x6JFi9oMcU1ICL4lNE+cOIH77rsPqqpC0zQMGTIEjz32GJKSkkKijSdOnMDChQtRU1MDg8GABx54ABMnTgyp782srCw8+uijuPrqq33PhUr71q5dizfeeAOS5F2pb/78+Zg8eXLItO/RRx/F999/D7fbjSuvvBILFy6E2WzW1b6gDBAiIup5QXkJi4iIeh4DhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEiXoJ4Li+h0l1xyie/fzc3NMJlMvvl9lixZgunTp/dUabplZmbiySefxBVXXNHTpRC1wQChkLJnzx7fv4PhF6/H44HB0LU/ht3xGtQ78RIW9QqapuH111/H5MmTkZGRgfvvvx81NTUAvIskDRs2DO+99x4mTpyI8ePH45133kFeXh6mTZuGcePGtZns8f3338ctt9yCpUuXYuzYscjOzsb27dt9X6+vr8fChQsxYcIEXHXVVXj++ed9cwy17Lts2TJkZGRg5cqVKCoqwm233YaMjAxkZGTgD3/4A+rq6gAA//Ef/4GSkhL8+7//Oy655BK88cYb2LFjR5s7wAFvWH7zzTcAgJUrV2L+/Pl46KGHcOmll+KDDz74xZqI9GKAUK/wj3/8A5s2bcKqVauwZcsWREdHt5sBeO/evdi4cSOef/55LFu2DK+99hreeustfPzxx/j000+xc+dO37Z5eXno378/vv32W8yfPx/33nuvL5AWLFgAg8GAjRs34sMPP8S2bduwZs2aNvv269cP27Ztw9133w0hBH73u99hy5Yt+PTTT1FaWoqVK1cCAJ599lmkpaXhtddew549ezBnzhy/2vv5558jOzsbu3fvxrRp085ZE5EeDBDqFd599108+OCDSElJgclkwr333osNGzbA4/H4tpk3bx7MZjMmTJiAsLAw5OTkID4+HsnJyRg3bhzy8/N928bFxeFf//VfYTQacd1112HQoEHYvHkzKisr8dVXX2HhwoUICwtDfHw8/u3f/g0ff/yxb9+kpCTceuutMBgMsFgsGDBgAK688kqYTCbExcXh9ttvx65du86rvWPGjMHkyZMhyzIaGhrOWRORHrwwSr1CSUkJ5s2bB1n++TOTLMuw2+2+x63XXjebze0et14zITk52TfZHgCkpaWhvLwcJSUl8Hg8mDBhgu9rmqa1mRY7JSWlTW2VlZV46qmnsHv3bjQ2NkIIgaioqPNqb+vX8KcmIj0YINQrpKSkYNmyZb51HlorLi7u8PHKysoghPCFiM1mQ2Zmpu8M59tvvz1rx3Xr4AGA5557DpIkYd26dYiJicGmTZt+cYEtq9UKh8Phe6yqKqqqqs76Gv7URKQHL2FRr/DP//zPeOGFF3Dy5EkA3sWqNm3apPt4VVVVePvtt+F2u/Hpp5+ioKAAEydORFJSEq688kr86U9/QkNDAzRNQ1FRUZv+k9M1NjYiLCwMkZGRKCsrw5tvvtnm6wkJCThx4oTv8aBBg+B0OrF582a43W68+uqrcLlcZz2+npqI/MEAoV7htttuQ2ZmJu644w5ccskluPnmm5GXl6f7eKNGjcLx48dx+eWX44UXXsCLL77oW1zpmWeegdvtxnXXXYfx48dj/vz5vkV8zuTee+9Ffn4+xo0bh7lz52LKlCltvj537ly8+uqrGDduHP7yl78gMjISixYtwmOPPYarr74aVqu13WWx03W0JiJ/cD0Qog56//33sWbNGrzzzjs9XQpRj+IZCBER6cIAISIiXXgJi4iIdOEZCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItLl/wFbCrmpbWMOWwAAAABJRU5ErkJggg==\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", + "id": "vanilla-baseline", + "metadata": {}, + "source": [ + "*I think I haven't 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, 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\"." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "innovative-compensation", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "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.9.0+" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/module4/requirements.txt b/module4/requirements.txt new file mode 100644 index 0000000000000000000000000000000000000000..44b08df60a73bb0516441ad8506e4462a5f05068 --- /dev/null +++ b/module4/requirements.txt @@ -0,0 +1,67 @@ +argon2-cffi==20.1.0 +async-generator==1.10 +attrs==20.3.0 +backcall==0.2.0 +bleach==3.3.0 +cffi==1.14.5 +cycler==0.10.0 +decorator==4.4.2 +defusedxml==0.7.1 +entrypoints==0.3 +ipykernel==5.5.0 +ipython==7.21.0 +ipython-genutils==0.2.0 +ipywidgets==7.6.3 +jedi==0.18.0 +Jinja2==2.11.3 +jsonschema==3.2.0 +jupyter==1.0.0 +jupyter-client==6.1.11 +jupyter-console==6.2.0 +jupyter-core==4.7.1 +jupyterlab-pygments==0.1.2 +jupyterlab-widgets==1.0.0 +kiwisolver==1.3.1 +MarkupSafe==1.1.1 +matplotlib==3.3.4 +mistune==0.8.4 +mpmath==1.2.1 +nbclient==0.5.3 +nbconvert==6.0.7 +nbformat==5.1.2 +nest-asyncio==1.5.1 +nose==1.3.7 +notebook==6.2.0 +numpy==1.20.1 +packaging==20.9 +pandas==1.2.3 +pandocfilters==1.4.3 +parso==0.8.1 +patsy==0.5.1 +pexpect==4.8.0 +pickleshare==0.7.5 +Pillow==8.1.2 +prometheus-client==0.9.0 +prompt-toolkit==3.0.17 +ptyprocess==0.7.0 +pycparser==2.20 +Pygments==2.8.1 +pyparsing==2.4.7 +pyrsistent==0.17.3 +python-dateutil==2.8.1 +pytz==2021.1 +pyzmq==22.0.3 +qtconsole==5.0.2 +QtPy==1.9.0 +scipy==1.6.1 +Send2Trash==1.5.0 +six==1.15.0 +statsmodels==0.12.2 +sympy==1.7.1 +terminado==0.9.2 +testpath==0.4.4 +tornado==6.1 +traitlets==5.0.5 +wcwidth==0.2.5 +webencodings==0.5.1 +widgetsnbextension==3.5.1