From 7c9760e2cd9a97dd8e75b218c7ba24c8b9a53f59 Mon Sep 17 00:00:00 2001 From: 28ec0a2489fc7a079d6e157047520de2 <28ec0a2489fc7a079d6e157047520de2@app-learninglab.inria.fr> Date: Fri, 16 Jun 2023 12:01:48 +0000 Subject: [PATCH] Upload New File --- module4/src_Python3_challenger.ipynb | 798 +++++++++++++++++++++++++++ 1 file changed, 798 insertions(+) create mode 100644 module4/src_Python3_challenger.ipynb diff --git a/module4/src_Python3_challenger.ipynb b/module4/src_Python3_challenger.ipynb new file mode 100644 index 0000000..35640a3 --- /dev/null +++ b/module4/src_Python3_challenger.ipynb @@ -0,0 +1,798 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this document we reperform some of the analysis provided in \n", + "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n", + "\n", + "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Technical information on the computer on which the analysis is run" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.6.4 |Anaconda, Inc.| (default, Mar 13 2018, 01:15:57) \n", + "[GCC 7.2.0]\n", + "uname_result(system='Linux', node='a1818438aec2', release='4.4.0-164-generic', version='#192-Ubuntu SMP Fri Sep 13 12:02:50 UTC 2019', machine='x86_64', processor='x86_64')\n", + "IPython 7.12.0\n", + "IPython.core.release 7.12.0\n", + "PIL 7.0.0\n", + "PIL.Image 7.0.0\n", + "PIL._version 7.0.0\n", + "_csv 1.0\n", + "_ctypes 1.1.0\n", + "_curses b'2.2'\n", + "decimal 1.70\n", + "argparse 1.1\n", + "backcall 0.1.0\n", + "cffi 1.13.2\n", + "csv 1.0\n", + "ctypes 1.1.0\n", + "cycler 0.10.0\n", + "dateutil 2.8.1\n", + "decimal 1.70\n", + "decorator 4.4.1\n", + "distutils 3.6.4\n", + "ipaddress 1.0\n", + "ipykernel 5.1.4\n", + "ipykernel._version 5.1.4\n", + "ipython_genutils 0.2.0\n", + "ipython_genutils._version 0.2.0\n", + "ipywidgets 7.2.1\n", + "ipywidgets._version 7.2.1\n", + "jedi 0.16.0\n", + "json 2.0.9\n", + "jupyter_client 6.0.0\n", + "jupyter_client._version 6.0.0\n", + "jupyter_core 4.6.3\n", + "jupyter_core.version 4.6.3\n", + "kiwisolver 1.1.0\n", + "logging 0.5.1.2\n", + "matplotlib 2.2.3\n", + "matplotlib.backends.backend_agg 2.2.3\n", + "numpy 1.15.2\n", + "numpy.core 1.15.2\n", + "numpy.core.multiarray 3.1\n", + "numpy.lib 1.15.2\n", + "numpy.linalg._umath_linalg b'0.1.5'\n", + "numpy.matlib 1.15.2\n", + "optparse 1.5.3\n", + "pandas 0.22.0\n", + "_libjson 1.33\n", + "parso 0.6.0\n", + "patsy 0.5.1\n", + "patsy.version 0.5.1\n", + "pexpect 4.8.0\n", + "pickleshare 0.7.5\n", + "platform 1.0.8\n", + "prompt_toolkit 3.0.3\n", + "ptyprocess 0.6.0\n", + "pygments 2.5.2\n", + "pyparsing 2.4.6\n", + "pytz 2019.3\n", + "re 2.2.1\n", + "scipy 1.1.0\n", + "scipy._lib.decorator 4.0.5\n", + "scipy._lib.six 1.2.0\n", + "scipy.fftpack._fftpack b'$Revision: $'\n", + "scipy.fftpack.convolve b'$Revision: $'\n", + "scipy.integrate._dop b'$Revision: $'\n", + "scipy.integrate._ode $Id$\n", + "scipy.integrate._odepack 1.9 \n", + "scipy.integrate._quadpack 1.13 \n", + "scipy.integrate.lsoda b'$Revision: $'\n", + "scipy.integrate.vode b'$Revision: $'\n", + "scipy.interpolate._fitpack 1.7 \n", + "scipy.interpolate.dfitpack b'$Revision: $'\n", + "scipy.linalg 0.4.9\n", + "scipy.linalg._fblas b'$Revision: $'\n", + "scipy.linalg._flapack b'$Revision: $'\n", + "scipy.linalg._flinalg b'$Revision: $'\n", + "scipy.ndimage 2.0\n", + "scipy.optimize._cobyla b'$Revision: $'\n", + "scipy.optimize._lbfgsb b'$Revision: $'\n", + "scipy.optimize._minpack 1.10 \n", + "scipy.optimize._nnls b'$Revision: $'\n", + "scipy.optimize._slsqp b'$Revision: $'\n", + "scipy.optimize.minpack2 b'$Revision: $'\n", + "scipy.signal.spline 0.2\n", + "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", + "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", + "scipy.special.specfun b'$Revision: $'\n", + "scipy.stats.mvn b'$Revision: $'\n", + "scipy.stats.statlib b'$Revision: $'\n", + "seaborn 0.8.1\n", + "seaborn.external.husl 2.1.0\n", + "seaborn.external.six 1.10.0\n", + "six 1.14.0\n", + "statsmodels 0.9.0\n", + "statsmodels.__init__ 0.9.0\n", + "traitlets 4.3.3\n", + "traitlets._version 4.3.3\n", + "urllib.request 3.6\n", + "zlib 1.0\n", + "zmq 17.1.2\n", + "zmq.sugar 17.1.2\n", + "zmq.sugar.version 17.1.2\n" + ] + } + ], + "source": [ + "def print_imported_modules():\n", + " import sys\n", + " for name, val in sorted(sys.modules.items()):\n", + " if(hasattr(val, '__version__')): \n", + " print(val.__name__, val.__version__)\n", + "# else:\n", + "# print(val.__name__, \"(unknown version)\")\n", + "def print_sys_info():\n", + " import sys\n", + " import platform\n", + " print(sys.version)\n", + " print(platform.uname())\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import statsmodels.api as sm\n", + "import seaborn as sns\n", + "\n", + "print_sys_info()\n", + "print_imported_modules()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loading and inspecting data\n", + "Let's start by reading data." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "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": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/master/data/shuttle.csv\")\n", + "data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAF9JJREFUeJzt3X2UXXV97/H3d5IACYlAg02VQAFJsVyBCOFJtDfx6Qa7JPUCBbyCl940ZUlul9y2htvVa6m1a1V8qHpFY+SiQldNVRBom14e1Ii0IASM4UHBuYBhEhogBshASGYy3/vH2bN7Mkxmzhlmz5lzeL/WmpWz9/mdne939pz5zN5nn9+JzESSJICuVhcgSZo8DAVJUslQkCSVDAVJUslQkCSVDAVJUqmyUIiIqyPiqYh4YC/3R0R8PiK6I2JDRJxQVS2SpMZUeaTwNWDxCPefAcwrvpYBX6qwFklSAyoLhcy8HfjlCEOWANdkzV3AgRHxuqrqkSSNbmoL/+9DgCfqlnuKdU8OHRgRy6gdTTB9+vQTDz300AkpsFEDAwN0dXXmyzOd2pt9tZ9O7W2i+nrkkUeeyczXjjaulaEQw6wbds6NzFwFrAJYsGBBrlu3rsq6mrZ27VoWLlzY6jIq0am92Vf76dTeJqqviPhFI+NaGbs9QP2f/HOBzS2qRZJEa0PhJuDC4iqkU4HnMvNlp44kSROnstNHEfENYCFwcET0AH8OTAPIzJXAGuA9QDfwInBRVbVIkhpTWShk5vmj3J/AJVX9/5Kk5nXeS/mSpDEzFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklSqNBQiYnFEPBwR3RFx2TD3HxAR/xARP4mIByPioirrkSSNrLJQiIgpwJXAGcAxwPkRccyQYZcAD2Xm8cBC4NMRsU9VNUmSRlblkcLJQHdmPpqZu4DVwJIhYxKYFREBzAR+CfRXWJMkaQSRmdVsOOJsYHFmLi2WLwBOyczldWNmATcBbwRmAedm5j8Ns61lwDKAOXPmnLh69epKah6r3t5eZs6c2eoyKtGpvdlX++nU3iaqr0WLFt2bmQtGGze1whpimHVDE+g/AeuBtwNvAG6NiB9m5vN7PChzFbAKYMGCBblw4cLxr/YVWLt2LZOtpvHSqb3ZV/vp1N4mW19Vnj7qAQ6tW54LbB4y5iLg+qzpBh6jdtQgSWqBKkPhHmBeRBxRvHh8HrVTRfU2Au8AiIg5wNHAoxXWJEkaQWWnjzKzPyKWAzcDU4CrM/PBiLi4uH8l8JfA1yLifmqnm1Zk5jNV1SRJGlmVrymQmWuANUPWray7vRl4d5U1SJIa5zuaJUklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVKo0FCJicUQ8HBHdEXHZXsYsjIj1EfFgRPygynokSSOb2sigiHhTZj7QzIYjYgpwJfAuoAe4JyJuysyH6sYcCHwRWJyZGyPiV5v5PyRJ46vRI4WVEXF3RHyo+EXeiJOB7sx8NDN3AauBJUPGvB+4PjM3AmTmUw1uW5JUgcjMxgZGzAN+DzgHuBv4ambeOsL4s6kdASwtli8ATsnM5XVjPgtMA/4DMAv4XGZeM8y2lgHLAObMmXPi6tWrG+tugvT29jJz5sxWl1GJTu3NvtpPp/Y2UX0tWrTo3sxcMOrAzGz4C5gCnAVsAn4K/Az4z3sZew5wVd3yBcD/HjLmC8BdwP7AwcDPgd8YqYYTTzwxJ5vvf//7rS6hMp3am321n07tbaL6AtZlA7/nG31N4TjgIuC3gVuB92bmfRHxeuBO4PphHtYDHFq3PBfYPMyYZzLzBeCFiLgdOB54pJG6JEnjq9HXFL4A3Accn5mXZOZ9AJm5GfizvTzmHmBeRBwREfsA5wE3DRlzI/C2iJgaETOAU6gdgUiSWqChIwXgPcCOzNwNEBFdwH6Z+WJmXjvcAzKzPyKWAzdTO+10dWY+GBEXF/evzMyfRsT/BTYAA9RONzV1lZMkafw0Ggq3Ae8EeovlGcAtwFtGelBmrgHWDFm3csjyJ4FPNliHJKlCjZ4+2i8zBwOB4vaMakqSJLVKo6HwQkScMLgQEScCO6opSZLUKo2ePvow8K2IGLx66HXAudWUJElqlYZCITPviYg3AkcDAfwsM/sqrUySNOEaPVIAOAk4vHjMmyOCHObdx5Kk9tXom9euBd4ArAd2F6sTMBQkqYM0eqSwADimeKu0JKlDNXr10QPAr1VZiCSp9Ro9UjgYeCgi7gZ2Dq7MzDMrqUqS1BKNhsLlVRYhSZocGr0k9QcR8evAvMy8rZi8bkq1pUmSJlpDrylExO8D3wa+XKw6BLihqqIkSa3R6AvNlwCnA88DZObPAT9PWZI6TKOhsDNrn7MMQERMpfY+BUlSB2k0FH4QEX8KTI+IdwHfAv6hurIkSa3QaChcBjwN3A/8AbXPSNjbJ65JktpUo1cfDQBfKb4kSR2q0bmPHmOY1xAy88hxr0iS1DLNzH00aD/gHOBXxr8cSVIrNfSaQmZurfvalJmfBd5ecW2SpAnW6OmjE+oWu6gdOcyqpCJJUss0evro03W3+4HHgd8d92okSS3V6NVHi6ouRJLUeo2ePvofI92fmZ8Zn3IkSa3UzNVHJwE3FcvvBW4HnqiiKElSazTzITsnZOZ2gIi4HPhWZi6tqjBJ0sRrdJqLw4Bddcu7gMPHvRpJUks1eqRwLXB3RHyH2jub3wdcU1lVkqSWaPTqo7+KiH8G3lasuigzf1xdWZKkVmj09BHADOD5zPwc0BMRR1RUkySpRRr9OM4/B1YA/7NYNQ3426qKkiS1RqNHCu8DzgReAMjMzTjNhSR1nEZDYVdmJsX02RGxf3UlSZJapdFQ+GZEfBk4MCJ+H7gNP3BHkjpOo1cffar4bObngaOBj2bmrZVWJkmacKMeKUTElIi4LTNvzcw/ycw/bjQQImJxRDwcEd0RcdkI406KiN0RcXYzxUuSxteooZCZu4EXI+KAZjYcEVOAK4EzgGOA8yPimL2M+wRwczPblySNv0bf0fwScH9E3EpxBRJAZv7hCI85GejOzEcBImI1sAR4aMi4/w5cR23CPUlSCzUaCv9UfDXjEPacRbUHOKV+QEQcQu1y17czQihExDJgGcCcOXNYu3Ztk6VUq7e3d9LVNF46tTf7aj+d2ttk62vEUIiIwzJzY2Z+fQzbjmHW5ZDlzwIrMnN3xHDDiwdlrgJWASxYsCAXLlw4hnKqs3btWiZbTeOlU3uzr/bTqb1Ntr5Ge03hhsEbEXFdk9vuAQ6tW54LbB4yZgGwOiIeB84GvhgRv9Pk/yNJGiejnT6q//P9yCa3fQ8wr5gjaRNwHvD++gGZWc6fFBFfA/4xM29AktQSo4VC7uX2qDKzPyKWU7uqaApwdWY+GBEXF/evbKpSSVLlRguF4yPieWpHDNOL2xTLmZmvGenBmbkGWDNk3bBhkJn/taGKJUmVGTEUMnPKRBUiSWq9Zj5PQZLU4QwFSVLJUJAklQwFSVLpVRMKW3t38pMnnmVr785WlyKpSVt7d7Kjb7fP3wnwqgiFG9dv4vRPfI8PXPUjTv/E97hp/aZWlySpQYPP38eefsHn7wTo+FDY2ruTFddt4KW+Abbv7OelvgE+ct0G/+KQ2kD983d3ps/fCdDxodCzbQfTuvZsc1pXFz3bdrSoIkmN8vk78To+FOYeNJ2+gYE91vUNDDD3oOktqkhSo3z+TryOD4XZM/flirOOY79pXczadyr7TeviirOOY/bMfVtdmqRR1D9/p0T4/J0AjX7ITls7c/4hnH7UwfRs28Hcg6b7AyW1kcHn79133sG/nPlWn78Ve1WEAtT+4vCHSWpPs2fuy/RpU3wOT4COP30kSWqcoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqRSpaEQEYsj4uGI6I6Iy4a5/79ExIbi618j4vgq65EkjayyUIiIKcCVwBnAMcD5EXHMkGGPAf8xM48D/hJYVVU9kqTRVXmkcDLQnZmPZuYuYDWwpH5AZv5rZm4rFu8C5lZYjyRpFJGZ1Ww44mxgcWYuLZYvAE7JzOV7Gf/HwBsHxw+5bxmwDGDOnDknrl69upKax6q3t5eZM2e2uoxKdGpv9tV+OrW3iepr0aJF92bmgtHGTa2whhhm3bAJFBGLgP8GvHW4+zNzFcWppQULFuTChQvHqcTxsXbtWiZbTeOlU3uzr/bTqb1Ntr6qDIUe4NC65bnA5qGDIuI44CrgjMzcWmE9kqRRVPmawj3AvIg4IiL2Ac4DbqofEBGHAdcDF2TmIxXWIklqQGVHCpnZHxHLgZuBKcDVmflgRFxc3L8S+CgwG/hiRAD0N3LOS5JUjSpPH5GZa4A1Q9atrLu9FHjZC8uCrb076dm2g7kHTWf2zH3HbWw76dS+qtK9ZTvbXuyje8t2jpozq9XlqE1VGgoamxvXb2LFdRuY1tVF38AAV5x1HGfOP+QVj20nndpXVT56w/1cc9dG/ujYfi79m9u58LTD+NiSY1tdltqQ01xMMlt7d7Liug281DfA9p39vNQ3wEeu28DW3p2vaGw76dS+qtK9ZTvX3LVxj3XX3LmR7i3bW1SR2pmhMMn0bNvBtK49d8u0ri56tu14RWPbSaf2VZX1Tzzb1HppJIbCJDP3oOn0DQzssa5vYIC5B01/RWPbSaf2VZX5hx7Y1HppJIbCJDN75r5ccdZx7Deti1n7TmW/aV1ccdZxw77Q2szYdtKpfVXlqDmzuPC0w/ZYd+Fph/lis8bEF5onoTPnH8LpRx3c0JU3zYxtJ53aV1U+tuRYLjz1cO6/9y5uu/RUA0FjZihMUrNn7tvwL8JmxraTTu2rKkfNmUXPjGkGgl4RTx9JkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpVGkoRMTiiHg4Iroj4rJh7o+I+Hxx/4aIOKHKeqRmbe3dyU+eeJatvTtHHbvusa185paHWffY1nHbZjNju7dsZ9uLfXRv2T7q2GZUVW+zNezo2z3qdru3bOfb657o2O9BFdsdampVG46IKcCVwLuAHuCeiLgpMx+qG3YGMK/4OgX4UvGv1HI3rt/Eius2MK2ri76BAa446zjOnH/IsGM/cNVd3NFdC4PPf6+btx01m2uXnvqKttnM2I/ecD/X3LWRPzq2n0v/5nYuPO0wPrbk2DF2Xn29Y6nhD3+zj0s/8b29bnfwezCoE78H473d4VR5pHAy0J2Zj2bmLmA1sGTImCXANVlzF3BgRLyuwpqkhmzt3cmK6zbwUt8A23f281LfAB+5bsOwf6Wte2xrGQiDfti99WVHDM1ss5mx3Vu27/HLEOCaOze+4r+Wq6p3rDXsztzrdl8t34Px3O7eRGZWs+GIs4HFmbm0WL4AOCUzl9eN+UfgrzPzjmL5u8CKzFw3ZFvLgGXF4tHAw5UUPXYHA8+0uoiKdGpvI/YV06bPmHrQ634jurqmDK7LgYHd/duefCT7drxYP3bKrINfP2X/A1/2x8zuF559cvf2ZzaPZZvNjO2accDsqa957eEAu198jikzDgCg//mnHx948bmRz2WNoKp6x1rDYG/Dbbf+e1CvTb4H4/azOIpfz8zXjjaostNHQAyzbmgCNTKGzFwFrBqPoqoQEesyc0Gr66hCp/bWyX31P/dUx/UFndvbZPtZrPL0UQ9waN3yXGDzGMZIkiZIlaFwDzAvIo6IiH2A84Cbhoy5CbiwuArpVOC5zHyywpokSSOo7PRRZvZHxHLgZmAKcHVmPhgRFxf3rwTWAO8BuoEXgYuqqqdik/bU1jjo1N7sq/10am+Tqq/KXmiWJLUf39EsSSoZCpKkkqEwBhHxeETcHxHrI2Jdse7yiNhUrFsfEe9pdZ3NiogDI+LbEfGziPhpRJwWEb8SEbdGxM+Lfw9qdZ3N2ktfnbC/jq6rf31EPB8RH273fTZCX52wzy6NiAcj4oGI+EZE7DfZ9pevKYxBRDwOLMjMZ+rWXQ70ZuanWlXXKxURXwd+mJlXFVeMzQD+FPhlZv51MX/VQZm5oqWFNmkvfX2YNt9f9YppZTZRmybmEtp8nw0a0tdFtPE+i4hDgDuAYzJzR0R8k9rFNscwifaXRwoCICJeA/wW8H8AMnNXZj5LbSqSrxfDvg78TmsqHJsR+uo07wD+X2b+gjbfZ0PU99UJpgLTI2IqtT9ONjPJ9pehMDYJ3BIR9xZTcAxaXsz2enWrDwHH4EjgaeCrEfHjiLgqIvYH5gy+d6T491dbWeQY7K0vaO/9NdR5wDeK2+2+z+rV9wVtvM8ycxPwKWAj8CS192XdwiTbX4bC2JyemSdQm+X1koj4LWozvL4BmE9th3+6hfWNxVTgBOBLmflm4AXgZdOdt6G99dXu+6tUnBI7E/hWq2sZT8P01db7rAixJcARwOuB/SPiA62t6uUMhTHIzM3Fv08B3wFOzswtmbk7MweAr1CbJbad9AA9mfmjYvnb1H6Zbhmcubb496kW1TdWw/bVAfur3hnAfZm5pVhu9302aI++OmCfvRN4LDOfzsw+4HrgLUyy/WUoNCki9o+IWYO3gXcDD8SeU36/D3igFfWNVWb+G/BERBxdrHoH8BC1qUg+WKz7IHBjC8obs7311e77a4jz2fMUS1vvszp79NUB+2wjcGpEzIiIoPaz+FMm2f7y6qMmRcSR1I4OoHZq4u8y868i4lpqh7UJPA78QbvN4xQR84GrgH2AR6ld7dEFfBM4jNoP9TmZ+cuWFTkGe+nr87T5/gKIiBnAE8CRmflcsW427b/PhuurE55jfwGcC/QDPwaWAjOZRPvLUJAklTx9JEkqGQqSpJKhIEkqGQqSpJKhIEkqVfbJa9JEKy7F/G6x+GvAbmpTXEDtDYa7WlLYCCLi94A1xfsppJbzklR1pMk0a21ETMnM3Xu57w5geWaub2J7UzOzf9wKlOp4+kivChHxwYi4u5iH/4sR0RURUyPi2Yj4ZETcFxE3R8QpEfGDiHh0cL7+iFgaEd8p7n84Iv6swe1+PCLuBk6OiL+IiHuKefRXRs251N6M9ffF4/eJiJ6IOLDY9qkRcVtx++MR8eWIuJXa5H5TI+Izxf+9ISKWTvx3VZ3IUFDHi4g3UZsW4S2ZOZ/aadPzirsPAG4pJjjcBVxObfqBc4CP1W3m5OIxJwDvj4j5DWz3vsw8OTPvBD6XmScBxxb3Lc7MvwfWA+dm5vwGTm+9GXhvZl4ALAOeysyTgZOoTcx42Fi+P1I9X1PQq8E7qf3iXFebcobp1KZQANiRmbcWt++nNp1xf0TcDxxet42bM3MbQETcALyV2vNnb9vdxb9PhwLwjoj4E2A/4GDgXuCfm+zjxsx8qbj9buA3I6I+hOZRmyZBGjNDQa8GAVydmf9rj5W1Dzqp/+t8ANhZd7v++TH0xbccZbs7snjBrpjH5wvUZmfdFBEfpxYOw+nn34/gh455YUhPH8rM7yKNI08f6dXgNuB3I+JgqF2lNIZTLe+O2mc9z6A2J/6/NLHd6dRC5pliht2z6u7bDsyqW34cOLG4XT9uqJuBDxUBNPi5xtOb7El6GY8U1PEy8/5idsrbIqIL6AMupvZRiI26A/g7ah/ycu3g1UKNbDczt0btc6IfAH4B/Kju7q8CV0XEDmqvW1wOfCUi/g24e4R6vkxtVs31xamrp6iFlfSKeEmqNIriyp43ZeaHW12LVDVPH0mSSh4pSJJKHilIkkqGgiSpZChIkkqGgiSpZChIkkr/HzHofwgP0tIHAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "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": 8, + "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", + "
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, 16 Jun 2023 Deviance: 3.0144
Time: 11:51:00 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, 16 Jun 2023 Deviance: 3.0144\n", + "Time: 11:51:00 Pearson chi2: 5.00\n", + "No. Iterations: 6 Covariance Type: nonrobust\n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", + "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 8, + "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": 9, + "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", + "
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, 16 Jun 2023 Deviance: 18.086
Time: 11:51:01 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, 16 Jun 2023 Deviance: 18.086\n", + "Time: 11:51:01 Pearson chi2: 30.0\n", + "No. Iterations: 6 Covariance Type: nonrobust\n", + "===============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "-------------------------------------------------------------------------------\n", + "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", + "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", + "===============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 9, + "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": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3dyb7QmLYISA7yA5hEXEBrYK2KiriinVBpHWp7SNVn199tE+16oNt1VZxQ3GpgisupYJa44JbQBBkX8UEkJ0kkD33748ZMGAgQzLJLPm8rivXzDlzn3O+dwY+c3LmnPuYcw4REYkunlAXICIiwadwFxGJQgp3EZEopHAXEYlCCncRkSikcBcRiUI1hruZPW1mW83s28O8bmb2sJmtMbPFZjYw+GWKiMjRCGTPfTow+givnwl09f9MBKbWvSwREamLGsPdOfcxsPMITc4FnnM+XwDpZtY6WAWKiMjRiwnCOtoC31eZzvXP23xoQzObiG/vnsTExKx27drVaoOVlZV4PNHxdYH6Ep6ipS/R0g9QX/ZbtWrVdudc85raBSPcrZp51Y5p4Jx7AngCYNCgQW7+/Pm12mB2djYjRoyo1bLhRn0JT9HSl2jpB6gv+5nZd4G0C8bHYC5QdRc8E9gUhPWKiEgtBSPc3wKu8J81czywxzn3k0MyIiLScGo8LGNmLwEjgGZmlgvcCcQCOOceA2YDZwFrgH3AVfVVrIiIBKbGcHfOXVLD6w64PmgViUhEKCsrIzc3l+Li4gbZXlpaGsuXL2+QbdW3QPqSkJBAZmYmsbGxtdpGML5QFZFGKDc3l9TUVDp06IBZdedVBFdBQQGpqan1vp2GUFNfnHPs2LGD3NxcOnbsWKttRMd5RSLS4IqLi2natGmDBHtjY2Y0bdq0Tn8VKdxFpNYU7PWnrr9bhbuISBTSMXcRiVher5c+ffocmJ41axYdOnQIXUFhROEuIhErMTGRRYsWHfb18vJyYmIaZ8zpsIyIRJXp06dz4YUXcvbZZ3PGGWcAMGXKFAYPHkzfvn258847D7S955576N69Oz/72c+45JJLeOCBBwAYMWIE+4dH2b59+4G/BioqKpg8efKBdT3++OPAj8MJjB07lh49enDZZZfhO0sccnJyOOGEE+jXrx9DhgyhoKCAUaNGHfShNHz4cBYvXhzU30Pj/EgTkaD649tLWbYpP6jr7NmmCXee3euIbYqKiujfvz8AHTt25I033gDg888/Z/HixWRkZDB37lxWr17NV199hXOOc845h48//pjk5GRmzJjBwoULKS8vZ+DAgWRlZR1xe9OmTSMtLY2cnBxKSkoYPnz4gQ+QhQsXsnTpUtq0acPw4cOZN28eQ4YM4aKLLmLmzJkMHjyY/Px8EhMTueKKK5g+fToPPvggq1atoqSkhL59+wbht/YjhbuIRKzDHZY5/fTTycjIAGDu3LnMnTuXAQMGAFBYWMjq1aspKCjgvPPOIykpCYBzzjmnxu3NnTuXxYsX8+qrrwKwZ88eVq9eTVxcHEOGDCEzMxOA/v37s2HDBtLS0mjdujWDBw8GoEmTJgCcd955DB8+nClTpvD0009z5ZVX1u0XUQ2Fu4jUWU172A0tOTn5wHPnHLfffjvXXXfdQW0efPDBw55uGBMTQ2VlJcBB55o75/j73//OqFGjDmqfnZ1NfHz8gWmv10t5eTnOuWq3kZSUxOmnn86bb77Jyy+/TG1HyD0SHXMXkag2atQonn76aQoLCwHIy8tj69atnHzyybzxxhsUFRVRUFDA22+/fWCZDh06sGDBAoADe+n71zV16lTKysoAWLVqFXv37j3stnv06MGmTZvIyckBfFemlpeXAzBhwgRuuukmBg8efOCvjGDSnruIRLUzzjiD5cuXM2zYMABSUlJ44YUXGDhwIBdddBH9+/fn2GOP5aSTTjqwzC233MK4ceN4/vnnOfXUUw/MnzBhAhs2bGDgwIE452jevDmzZs067Lbj4uKYOXMmN954I0VFRSQmJvL+++8DkJWVRZMmTbjqqnoaa9E5F5KfrKwsV1sffvhhrZcNN+pLeIqWvtRnP5YtW1Zv665Ofn5+va7/zjvvdFOmTKnXbeyXn5/v8vLyXNeuXV1FRcVh21X3OwbmuwAyVodlREQa2IsvvsjQoUO555576u3WgTosIyIC3HXXXQ22rUsvvfQnX/AGm/bcRaTWnKv2dskSBHX93SrcRaRWEhIS2LFjhwK+Hjj/eO4JCQm1XocOy4hIrWRmZpKbm8u2bdsaZHvFxcV1CrtwEkhf9t+JqbYU7iJSK7GxsbW+S1BtZGdnH7jKNNI1RF90WEZEJAop3EVEopDCXUQkCincRUSikMJdRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSiUMSF+77Sct7bUEZ5RWWoSxERCVsRF+7vLN7MP1eUMu7xz/lux+HvXSgi0phFXLiPG9SOSX3jWbO1kDMf+oSZORs15KiIyCEiLtwBjm8Tw7s3n0z/dunc+toSbnhxIXuKykJdlohI2IjIcAdok57IC9cM5dbRPZizdAtnPfQJX2/cFeqyRETCQsSGO4DHY/xqRGde/dUJeDww7rHPefLjdTpMIyKNXkDhbmajzWylma0xs9uqeT3NzN42s2/MbKmZXRX8Ug+vf7t03rnxJE47rgX3zF7OxOcX6DCNiDRqNYa7mXmBR4AzgZ7AJWbW85Bm1wPLnHP9gBHAX8wsLsi1HlFaYiyPXZ7F//yiJx+u2Mq5//iUFVvyG7IEEZGwEcie+xBgjXNunXOuFJgBnHtIGwekmpkBKcBOoDyolQbAzLj6xI7MmHg8+0orGPPIPN76ZlNDlyEiEnJW0/FpMxsLjHbOTfBPjweGOuduqNImFXgL6AGkAhc55/5VzbomAhMBWrZsmTVjxoxaFV1YWEhKSsoR2+wuqeTRRSWs2lXJWR1jGdstFo9ZrbZXnwLpS6RQX8JPtPQD1Jf9Ro4cucA5N6jGhs65I/4AFwJPVZkeD/z9kDZjgb8BBnQB1gNNjrTerKwsV1sffvhhQO1Kyirc/3tjsTv21nfcFdO+dHuKSmu9zfoSaF8igfoSfqKlH86pL/sB810Nue2cC+iwTC7Qrsp0JnDosY6rgNf9217jD/ceAay7XsXFeLh7TB/uPb8P89Zs5/xHP9NVrSLSKAQS7jlAVzPr6P+S9GJ8h2Cq2gicBmBmLYHuwLpgFloXlwxpz/PXDGV7YQnnPjKPL9btCHVJIiL1qsZwd86VAzcAc4DlwMvOuaVmNsnMJvmb/Qk4wcyWAB8AtzrnttdX0bUxrHNT3rx+OE2T4xg/7UveWJgb6pJEROpNTCCNnHOzgdmHzHusyvNNwBnBLS34jm2azOu/Gs6kFxbw25nfsHFHETed1gULwy9aRUTqIqKvUK2NtKRYnr16COcPbMvf3l/Fba8t0fDBIhJ1AtpzjzZxMR7+cmE/MtMTefg/a9hWWMI/Lh1AUlyj/HWISBRqdHvu+5kZvzujO/ec15vslVu59Mkv2bW3NNRliYgERaMN9/0uG3osUy/PYtnmfC58/HM27S4KdUkiInXW6MMdYFSvVjx71RC27Clm7NTPWLetMNQliYjUicLdb1jnpsyYeDzF5ZWMe/xzlm/WoGMiErkU7lX0bpvGy9cNI8bj4aLHP2ehbv4hIhFK4X6ILi1SeGXSMNKT4hg/7StyNuwMdUkiIkdN4V6NdhlJvHzdMFo0ieeKaV/x2ZqwuthWRKRGCvfDaJWWwMyJw2ifkcRV03P4eNW2UJckIhIwhfsRNE+N56WJx9OxWTITnpuvgBeRiKFwr0FGchwvXns8nZunMOG5+XykgBeRCKBwD0BGchwvThhK5+YpTHxuPvN0DF5EwpzCPUDHJMfxzwlD6dA0mWuezdGY8CIS1hTuRyEjOY5/XjuUzGOSuHp6Dgu+02mSIhKeFO5HqVlKPC9eO5SWTRK48ukcvs3bE+qSRER+QuFeCy1SE/jnhKE0SYxl/LQvWbmlINQliYgcROFeS23SE3nx2qHExXi4fNqXuvG2iIQVhXsdHNs0mReuGUp5RSWXPfUlW/YUh7okERFA4V5nXVum8uzVQ9i9r4zLp+mGHyISHhTuQdA3M50nrxjExp37uHJ6DntLykNdkog0cgr3IBnWuSn/uGQAS3J3M+mFBZSW66bbIhI6CvcgOqNXK+47vy+frN7Of73yDZWVLtQliUgjFRPqAqLNuMHt2LG3lPvfXUHzlHju+MVxmFmoyxKRRkbhXg8mndKJrQXFPD1vPS2bxHPdKZ1DXZKINDIK93pgZtzx855sKyjh3n+voEWTeM4bkBnqskSkEVG41xOPx/jLuH7sKCzl968upkVqAsO7NAt1WSLSSOgL1XoUH+PlsfFZdGqWwqTnF7B8c36oSxKRRkLhXs/SEmN55qrBJMfHcNUzOWzeUxTqkkSkEVC4N4A26Yk8c9VgCkvKueqZHAqKy0JdkohEOYV7AzmudRMevWwgq7cWcv2LCymr0EVOIlJ/FO4N6ORuzfnzeb35eNU2/ufNpTini5xEpH7obJkGdtHg9ny3Yx+PZq+lY7MkuoW6IBGJStpzD4FbzujOz/u05t5/r2D+Fg0yJiLBF1C4m9loM1tpZmvM7LbDtBlhZovMbKmZfRTcMqPL/nPg+7dL54nFJSzO3R3qkkQkytQY7mbmBR4BzgR6ApeYWc9D2qQDjwLnOOd6ARfWQ61RJSHWyxPjB5EaZ0x4dr5OkRSRoApkz30IsMY5t845VwrMAM49pM2lwOvOuY0AzrmtwS0zOjVPjee3WQnsK63gmunzNQ68iASN1XTGhpmNBUY75yb4p8cDQ51zN1Rp8yAQC/QCUoGHnHPPVbOuicBEgJYtW2bNmDGjVkUXFhaSkpJSq2XDTWFhIeuKEvjbghL6t/By44B4PBE6imS0vS/R0Jdo6QeoL/uNHDlygXNuUE3tAjlbprqkOfQTIQbIAk4DEoHPzewL59yqgxZy7gngCYBBgwa5ESNGBLD5n8rOzqa2y4ab7OxsbvrFCJq0Wc9dby8jp6Q1t47uEeqyaiXa3pdo6Eu09APUl6MVSLjnAu2qTGcCm6pps905txfYa2YfA/2AVUhAfnlCB1ZvLWRq9lq6NE/hgiyNIikitRfIMfccoKuZdTSzOOBi4K1D2rwJnGRmMWaWBAwFlge31OhmZtx1Ti9O6NyU219fwoLvdoa6JBGJYDWGu3OuHLgBmIMvsF92zi01s0lmNsnfZjnwLrAY+Ap4yjn3bf2VHZ1ivR4evWwgrdMTuO75BeTt1hk0IlI7AZ3n7pyb7Zzr5pzr7Jy7xz/vMefcY1XaTHHO9XTO9XbOPVhfBUe79KQ4pv1yECVllUx4VmfQiEjt6ArVMNSlRSoPXzqAlVvy+a+XdaNtETl6CvcwNbJ7C/77rON4d+kWHv7P6lCXIyIRRgOHhbFrTuzIii0FPPj+arq3TOXMPq1DXZKIRAjtuYcxM+Oe83ozsH06v3v5G5Zt0m36RCQwCvcwt/8+rGmJsVz73Hx2FJaEuiQRiQAK9wjQIjWBJ67IYnthCb/+59e6i5OI1EjhHiH6ZqZz/wV9+XL9Tv749tJQlyMiYU5fqEaQMQPasnxzPo9/vI6erdO4dGj7UJckImFKe+4R5veje3BKt+bc+da3zN+gIQpEpHoK9wjj9RgPXzyAzGOSmPTCAjZpiAIRqYbCPQKlJcXy5BVZFJdVct3zCyguqwh1SSISZhTuEapLi1T+dlF/luTt4fbXl1DTTVdEpHFRuEew03u25Hend+ONhXlM+3R9qMsRkTCicI9wN4zswuherfjz7OV8snpbqMsRkTChcI9wHo/xwLh+dGmRwg0vLmTjjn2hLklEwoDCPQqkxMfw5BW+++Ve+5zGgBcRhXvUOLZpMv+4dACrtxZwyyvf6AtWkUZO4R5FTuranNvO7MG/v93CIx+uCXU5IhJCCvcoc+1JnRjTvw1/eW8V7y/7IdTliEiIKNyjjJlx3wV96dWmCTfPXMSarYWhLklEQkDhHoUSYr08Pn4Q8TEeJj43nz1FZaEuSUQamMI9SrVNT2Tq5Vls3LmP38xYSIVusi3SqCjco9iQjhnceU4vsldu44G5K0Ndjog0II3nHuUuH9qeZZvymZq9lp6tm3B2vzahLklEGoD23KOcmfHHc3ox6NhjmPzqN3ybtyfUJYlIA1C4NwJxMR6mXp7FMUlxTHxuPtt1k22RqKdwbySap8bz5BWD2LG3lF+9sIDSct1kWySaKdwbkd5t0/i/sX3J2bCLO99aqiEKRKKYvlBtZM7t35blmwt47KO19GydyvhhHUJdkojUA+25N0KTR3Xn1B4t+OPby/hs7fZQlyMi9UDh3gh5PcZDF/enQ7Nkrv/n1xoDXiQKKdwbqdSEWJ66YhCVDiY8l0NBsYYoEIkmCvdGrEOzZB69bCBrt+3l5hmLNESBSBRRuDdyw7s0486ze/LBiq1MmaMhCkSihc6WEcYffywrt/jOoOnWMoXzB2aGuiQRqaOA9tzNbLSZrTSzNWZ22xHaDTazCjMbG7wSpb6ZGXed04thnZpy22tLWPDdzlCXJCJ1VGO4m5kXeAQ4E+gJXGJmPQ/T7n5gTrCLlPoX6/Uw9fKBtElPYOJzC8jdpTNoRCJZIHvuQ4A1zrl1zrlSYAZwbjXtbgReA7YGsT5pQOlJcUy7cjClFZVMeHa+zqARiWBW0yXo/kMso51zE/zT44GhzrkbqrRpC7wInApMA95xzr1azbomAhMBWrZsmTVjxoxaFV1YWEhKSkqtlg034diXpdsr+MuCYno383LzwHg8ZgEtF459qa1o6Uu09APUl/1Gjhy5wDk3qKZ2gXyhWt3/7EM/ER4EbnXOVdgRgsA59wTwBMCgQYPciBEjAtj8T2VnZ1PbZcNNOPZlBJCW+R1/mPUtnxS24M6zewW0XDj2pbaipS/R0g9QX45WIOGeC7SrMp0JbDqkzSBghj/YmwFnmVm5c25WUKqUBnf58ceybttenp63nk7NkjUGjUiECSTcc4CuZtYRyAMuBi6t2sA513H/czObju+wjII9wv2/nx/Hdzv2cudbS8nMSGJk9xahLklEAlTjF6rOuXLgBnxnwSwHXnbOLTWzSWY2qb4LlNDxeoyHLxnAca2bcMM/v2bZpvxQlyQiAQroPHfn3GznXDfnXGfn3D3+eY855x6rpu2V1X2ZKpEpOT6Gab8cTGpCLNc8m8OWPcWhLklEAqDhB6RGrdISmHblIPKLyrh6eg6FJeWhLklEaqBwl4D0apPGI5cNZOUPBdzw4teUV+g2fSLhTOEuARvRvQV/Orc32Su3cceb3+o2fSJhTAOHyVG5dGh7cnft49HstbRNT+SGU7uGuiQRqYbCXY7a5FHd2bKnmAfmrqJVWiJjszSKpEi4UbjLUTMz7rugL1sLSrjttcU0T43nlG7Na7WuWQvzmDJnJZt2F9EmPZHJo7ozZkDbIFcs9UXvX/jSMXeplbgY3yiS3Vqm8qsXFvDN97uPeh2zFuZx++tLyNtdhAPydhdx++tLmLUwL/gFS9Dp/QtvCneptdSEWKZfPZimKXFcPT2HLXuP7gyaKXNWUlRWcdC8orIK3REqQuj9C28Kd6mTFqkJPHvVEBzwwPxifsgP/CKnTbuLjmq+hBe9f+FN4S511ql5CtOvGkxhqeOKaV+xZ19g48C3SU88qvkSXvT+hTeFuwRF38x0bhqYwPrte7n62Rz2ldZ8FevkUd1JjPUeNC8x1svkUd3rq0wJIr1/4U3hLkHTs6mXBy/uz8KNu5j0wteUlh/5GPyYAW259/w+tE1PxIC26Ynce34fnW0RIfT+hTedCilBdVaf1tx3fl9+/9pifjtzEQ9fMgCv5/A3cBkzoK3CIILp/QtfCncJunGD25FfXMbd/1pOcryX+87vi+cIAS8iwadwl3ox4aROFBSX89AHq0mM9XLXOb040i0YRSS4FO5Sb27+WVf2lZbz5CfrSYjzctvoHgp4kQaicJd6Y2b891nHsa+0gsc/Wkd8jJffnd4t1GWJNAoKd6lXZsafzu1NaXklD3+wmliPceNpGklSpL4p3KXeeTy+gcYqKh1/eW8VMV4PvxrROdRliUQ1hbs0CK/HmHJhP8orHfe/uwKH49cjuoS6LJGopXCXBuP1GH8d1w8z+L93V+IcXD9SAS9SHxTu0qBivB7+Oq4/hm9UwfIKx02nddFZNCJBpnCXBuf1GH8Z1x+vx8Pf3l9FSXkFk0d1V8CLBJHCXULC6zGmjO1LXIyHR7PXUlxWyR2/OE4BLxIkCncJGY/H+PN5vYmP8fD0vPXsKy3nnvP6HHEsGhEJjMJdQsrMuPPsnqQmxPD3/6yhsKScv47rT1yMBiwVqQuFu4ScmfFfZ3QnJT6Ge/+9goLicqZePpCkOP3zFKkt7R5J2LjulM7cd34fPlm9jUuf/JJde0tDXZJIxFK4S1i5eEh7Hr0si2Wb87nw8c/J0/04RWpF4S5hZ3TvVjx71RB+yC/mvEfmsWxTfqhLEok4CncJS8M6N+WVScPwmDHu8c/5ZPW2UJckElEU7hK2erRqwhvXn0DmMYlc9UwOM3M2hrokkYihcJew1jotkVcmDWNY56bc+toS7n93BZWVLtRliYQ9hbuEvdSEWJ6+cjCXDGnP1Oy1/OqfC9hXWh7qskTCWkDhbmajzWylma0xs9uqef0yM1vs//nMzPoFv1RpzGK9Hv58Xm/u+EVP3lv2A2Onfs4mnUkjclg1hruZeYFHgDOBnsAlZtbzkGbrgVOcc32BPwFPBLtQETPjmhM7Mu3KwXy/cx9n//1Tvlq/M9RliYSlQPbchwBrnHPrnHOlwAzg3KoNnHOfOed2+Se/ADKDW6bIj0Z2b8Eb1w8nLTGWS5/8guc/34BzOg4vUpXV9J/CzMYCo51zE/zT44GhzrkbDtP+FqDH/vaHvDYRmAjQsmXLrBkzZtSq6MLCQlJSUmq1bLhRX2pvX5nj8cUlfLOtguFtYriiVxzx3uAMOhYt70u09APUl/1Gjhy5wDk3qMaGzrkj/gAXAk9VmR4P/P0wbUcCy4GmNa03KyvL1daHH35Y62XDjfpSNxUVle5v7610HW57x43620du/bbCoKw3Wt6XaOmHc+rLfsB8V0O+OucCOiyTC7SrMp0JbDq0kZn1BZ4CznXO7QhgvSJ15vEYN/+sG89cOZjNe4o5+++f8q/Fm0NdlkjIBRLuOUBXM+toZnHAxcBbVRuYWXvgdWC8c25V8MsUObIR3Vvwr5tOpHOLFK5/8WvumPUtxWUVoS5LJGRqDHfnXDlwAzAH3yGXl51zS81skplN8jf7H6Ap8KiZLTKz+fVWschhZB6TxMvXDePakzry/BffMeaReaz+oSDUZYmEREADZjvnZgOzD5n3WJXnE4CffIEq0tBmL9nM7CVbAFj1QwFnPfwJY/q3Zd6a7WzeU0yb9EQmj+rOmAFtg77tWQvzmDJnJZt2F9XrdgLxh1lLeOnL77m5dxnX3D6bS4a24+4xfUJSi4SG7oYgUWPWwjxuf30JRf7DMZUOXIXjlQW5B9rk7S7i9teXAAQ1eA/ddn1tJxB/mLWEF774cRyeCucOTCvgGw8NPyBRY8qclQfCdb/qTvQtKqtgypyV9b7t+thOIF768vujmi/RSeEuUeNohiMI9k1ADrftUAyRUHGYa1cON1+ik8Jdokab9MSA23oM3lyUF7QrWw+37aOpKVi8Vv2FXIebL9FJ4S5RY/Ko7iTGeg+aF+sxYg+5ajU+xkPmMUn8ZsYiJjw7n9xd++pl24mxXiaP6l7ndR+tS4a2O6r5Ep0U7hI1xgxoy73n96FteiIGtE1PZMqF/Zgytt9B8+6/oC8f3jKCP/z8OD5bu4PT//oxT3y8lrKKyqBu+97z+4TkbJm7x/Th8uPbH9hT95px+fHt9WVqI6OzZSSqjBnQttpArW7ehJM6Mbp3K+56ayl/nr2C1xbk8cdze3F8p6ZB3XYo3D2mD3eP6UN2djZrLxsR6nIkBLTnLo1a5jFJPPXLwTwxPovCknIufuILbnppITuLa78XLxIOtOcuApzRqxUndW3O1I/W8thHa3m3spLvvKu47pROJMXpv4lEHu25i/glxnn53end+OB3p9C/hZeHPljNyAeymZmzkQrdt1UijMJd5BDtMpL4df8EXpk0jNZpidz62hLOfOhj3lv2g24KIhFD4S5yGIM7ZPDGr0/g0csGUlbhuPa5+Yx59DM+Xb1dIS9hT+EucgRmxll9WvPeb0/m/gv6sC2/mMunfcm4xz9n3hqFvIQvhbtIAGK8Hi4a3J7/3DKC/z23F9/vLOKyp77kgqmf8Z8VOlwj4UfhLnIUEmK9XDGsA9mTR/Cnc3vxQ34JV0+fz1kPf8obC3PrdCGUSDAp3EVqISHWy3h/yD9wYT/KKir57cxvOPn/PuTxj9ayZ19ZqEuURk4n8IrUQazXw9isTM4f0JbsVVt54uN13PvvFTz4/mouyGrLFcM60K1laqjLlEZI4S4SBB6PcWqPlpzaoyVLN+1h+rwNvDw/lxe+2MjQjhlcfvyxnNGrJfEx3ppXJhIECneRIOvVJo0pF/bj9rOO4+X53/PCF99x40sLyUiO44KBbRk3qB1dtTcv9UzhLlJPMpLjmHRKZyae1IlP12znpa828sy8DTz5yXr6ZaYxNiuTn/dtQ0ZyXKhLlSikcBepZx6PcXK35pzcrTnbC0uYtTCPVxfkcsebS/nj28sY0b05Z/drw8+Oa0lyvP5LSnDoX5JIA2qWEs+Ekzox4aROLN+cz6yFeby5aBPvL99KQqyHU3u04MzerRnZowUpCnqpA/3rEQmR41o34bjWTbh1dA/mf7eLt7/ZxL+/3cLsJVuIi/FwUpdmnN6zJacd15LmqfGhLlcijMJdJMQ8HmNIxwyGdMzgrnN6seC7Xcxespn3lv3AByu2YraEvpnpnNajBSO6N6d3mzQ8Ht0PVY5M4S4SRrxVgv7Os3uyfHMBHyz3hfzf3l/FX99bRUZyHCd2acaJXZtxYpc9/m12AAANA0lEQVRmIbkJt4Q/hbtImDIzerZpQs82TbjxtK5sLyzh09XbyV65lU/XbOetbzYB0LFZMsd3yuD4Tk0Z0jGD1mkKe1G4i0SMZinxB+7T6pxj5Q8FfLp6O5+v3cE732zmpa++B6BdRiKDj81g4LHH4Aoqqah0eHUYp9FRuItEIDOjR6sm9GjVhAkndaK8opLlmwv4asNOvlq/g49Xb+P1hXkA3Jczh76Z6fRrl06/zDT6tkunTVoCZgr8aKZwF4kCMV4PfTLT6JOZxjUndsQ5x8ad+3jh3c8oSWnNwo27mfbpOsoqfEMTH5MUS++2afRs7Tvsc1zrJnRslkysV2MJRguFu0gUMjOObZrM8LaxjBjRG4CS8gpWbC5gce5ulm7K59tNe3hm3gZK/cMUx3k9dGqeTPdWqXRtkULXlql0aZFC+4wkhX4EUriLNBLxMV7foZl26QfmlVVUsm7bXpZvzmfFlgJWbMln/oZdvLlo04E2sV6jfUYSHZul0Kl5Msc2TaJj02TaN02idVqijueHKYW7SCMW6/XQvVUq3VsdPJBZYUk5a7YWsnZrIWu2FbJ+217Wb9/Lx6u3UVpeWWV5o216Iu0yksg8JonMYxJpm55I22MSaZOeSMvUeGK01x8SCncR+YmU+Bj6t0unf5W9fIDKSseW/GI2bN/Lhh37+H7XPjbu2Efurn3M3bSFHXtLD2rvMWiRmkCrtARaNfE9tmgST4vUBFqkxtOiSTzNUuLJSIrThVlBpnAXkYB5PEabdN9e+Qldfvr6vtJyNu0uIndXEZv3FLN5dxGb9hTzQ34xa7YVMm/tdgqKy3+ynNdjZCTH0TQ5jmYp8TRNiSMjOY6MpDgyUuI4JimO73ZU0HJzPsckxZGWGEtCrEdn/ByBwl1EgiYpLoYuLVLp0uLw49XvKy1na34JWwtK2FpQzPaCErYXlrKtoIQde0vZXljCxp372LW3lIKSgz8I7s/55MDzOK+HtKRYmiTEkJYYS5PEWJokxJKaEEPqgccYUuJ//En2//iee0mKi4na7wwCCnczGw08BHiBp5xz9x3yuvlfPwvYB1zpnPs6yLWKRK1ZC/OYMmclm3YX0SY9kcmjuvPK/I3MW7vzQJvhnTO4cFD7n7QDfjJv/nc7eenL77m5dxnX3D6bS4a24+4xfQLa7pgBbQ87P5Dl92+7wjm8Zj/ZdlJcDB2axbDo+9019uXOs3tycrfm7NxXyofzcujQrSe79pXx2drtZK/cxraCEgqKy/B6jLIKx4bte8kvLqeguOzAaZ81iY/xkBTnC/qkOC+JcV4SY70kxXlJiPU9T4jzkhDjJSHWQ0Lsj4/xMR7iY/yPsT8+j9v/4/3xebzXS2yMEev14FxgtdVFjeFuZl7gEeB0IBfIMbO3nHPLqjQ7E+jq/xkKTPU/ikgNZi3M4/bXl1BUVgFA3u4ibp656Cft5q3deVDY5+0uYvKr34CDskp3YN7vZi6isspyFc7xwhcbAQ4K2eq2e/vrS5j/3U5eW5D3k/nAQQFf3fJ12fbkV74B40Ao5+0u4o43l3Lv+X0YM6AtW5p6GdGnNbMW5vHB8q0Hli0uq+T7nUUH2u1XXFZBYUk5hcXlFBSXU1hSzt6ScvaWlrO3pIK9JeXsK61gb2k5+0p9z4tKKw487thbemC6pLyC4rJKisoqqKisezCf2TGWkSPrvJojCmTPfQiwxjm3DsDMZgDnAlXD/VzgOef7OPrCzNLNrLVzbnPQKxaJMlPmrDwQVEerur3TymraAbz05fcHBWx12y0qqziw133o/ClzVh4UntUtX5dtl1UTmoFut7p2vj1sL81SgjtccllFJcVlFZSUV1JS7nte6n/ue6ygpKyS0grf9IHH8krKKnw/nl0bg1pTdaymPw/MbCww2jk3wT89HhjqnLuhSpt3gPucc5/6pz8AbnXOzT9kXROBif7J7sDKWtbdDNhey2XDjfoSnhqsL3GtumTV17or9u3Bm5R2YLp0y5oFddluXZYPwrLNgO1HWrbqNsJcXf59Heuca15To0D23Kv7tuHQT4RA2uCcewJ4IoBtHrkgs/nOuUF1XU84UF/CU7T0xczml+/ZGvH9gOh5T6Bh+hLI1QW5QLsq05nAplq0ERGRBhJIuOcAXc2so5nFARcDbx3S5i3gCvM5Htij4+0iIqFT42EZ51y5md0AzMF3KuTTzrmlZjbJ//pjwGx8p0GuwXcq5FX1VzIQhEM7YUR9CU/R0pdo6QeoL0elxi9URUQk8mhEHxGRKKRwFxGJQmEf7maWYGZfmdk3ZrbUzP7on59hZu+Z2Wr/4zGhrjUQZuY1s4X+awMiuR8bzGyJmS0ys/n+eZHal3Qze9XMVpjZcjMbFol9MbPu/vdj/0++md0coX35rf//+7dm9pI/ByKuHwBm9ht/P5aa2c3+efXel7APd6AEONU51w/oD4z2n5FzG/CBc64r8IF/OhL8BlheZTpS+wEw0jnXv8r5upHal4eAd51zPYB++N6fiOuLc26l//3oD2ThO7nhDSKsL2bWFrgJGOSc643vRI6LibB+AJhZb+BafFf69wN+YWZdaYi+OOci5gdIAr7GN27NSqC1f35rYGWo6wug/kz/G3kq8I5/XsT1w1/rBqDZIfMiri9AE2A9/pMLIrkvh9R/BjAvEvsCtAW+BzLwndH3jr8/EdUPf50X4htscf/0HcDvG6IvkbDnvv9QxiJgK/Cec+5LoKXzn0vvf2wRyhoD9CC+N7bqEByR2A/wXYE818wW+IeVgMjsSydgG/CM/3DZU2aWTGT2paqLgZf8zyOqL865POABYCOwGd91M3OJsH74fQucbGZNzSwJ3ynj7WiAvkREuDvnKpzvT81MYIj/T52IYma/ALY65yJl7IuaDHfODcQ3Iuj1ZnZyqAuqpRhgIDDVOTcA2EsE/Ll/JP6LDc8BXgl1LbXhP/58LtARaAMkm9nloa2qdpxzy4H7gfeAd4FvgJ/eraQeRES47+ec2w1kA6OBH8ysNYD/cWsISwvEcOAcM9sAzABONbMXiLx+AOCc2+R/3IrvuO4QIrMvuUCu/69BgFfxhX0k9mW/M4GvnXM/+KcjrS8/A9Y757Y558qA14ETiLx+AOCcm+acG+icOxnYCaymAfoS9uFuZs3NLN3/PBHfG78C35AHv/Q3+yXwZmgqDIxz7nbnXKZzrgO+P5n/45y7nAjrB4CZJZtZ6v7n+I6HfksE9sU5twX43sy6+2edhm8464jrSxWX8OMhGYi8vmwEjjezJDMzfO/JciKvHwCYWQv/Y3vgfHzvTb33JeyvUDWzvsCz+L4x9wAvO+f+18yaAi8D7fH9Y7jQObfz8GsKH2Y2ArjFOfeLSOyHmXXCt7cOvsMaLzrn7onEvgCYWX/gKSAOWIdv+AwPkdmXJHxfRnZyzu3xz4u498V/yvNF+A5hLAQmAClEWD8AzOwToClQBvzOOfdBQ7wnYR/uIiJy9ML+sIyIiBw9hbuISBRSuIuIRCGFu4hIFFK4i4hEoUBukC3SoPyniX3gn2wFVOAbIgBgiHOuNCSFHYGZXQ3M9p83LxJyOhVSwpqZ3QUUOuceCINavM65isO89ilwg3Nu0VGsL8Y51yCXokvjo8MyElHM7JfmG99/kZk9amYeM4sxs91mNsXMvjazOWY21Mw+MrN1ZnaWf9kJZvaG//WVZvaHANd7t5l9hW9coz+aWY5/fO7HzOcifMNRz/QvH2dmuVWurD7ezN73P7/bzB43s/fwDVYWY2Z/9W97sZlNaPjfqkQjhbtEDP+AcecBJ/gHkovBN5QDQBow1z+YWSlwF77L1i8E/rfKaob4lxkIXGpm/QNY79fOuSHOuc+Bh5xzg4E+/tdGO+dmAouAi5xvPPWaDhsNAM52zo0HJuIbUG4IMBjfIGzta/P7EalKx9wlkvwMXwDO9w05QiK+S+0Bipxz7/mfL8E3TGy5mS0BOlRZxxzn3C4AM5sFnIjv/8Hh1lvKj0MtAJxmZpOBBKAZsAD491H2403nXLH/+RnAcWZW9cOkK75L0kVqTeEukcSAp51zdxw00ywGXwjvV4nvDl77n1f9d37ol0yuhvUWOf8XU/5xW/4BDHTO5ZnZ3fhCvjrl/PiX8aFt9h7Sp1875z5AJIh0WEYiyfvAODNrBr6zampxCOMM890zNQnfmOHzjmK9ifg+LLb7R8W8oMprBUBqlekN+G51xyHtDjUH+LX/g2T/fVATj7JPIj+hPXeJGM65Jf7RAt83Mw++UfYmAZuOYjWfAi8CnYHn95/dEsh6nXM7zOxZfMMbfwd8WeXlZ4CnzKwI33H9u4AnzWwL8NUR6nkc38iAi/yHhLbi+9ARqROdCimNhv9MlN7OuZtDXYtIfdNhGRGRKKQ9dxGRKKQ9dxGRKKRwFxGJQgp3EZEopHAXEYlCCncRkSj0/wHRUJwHFwSFegAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "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": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/conda/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", + " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9+P/XmTWTSUIWsgAJAWLYES2IoiKKAsquSBWoWhWrttpva71XvbXeaq9a76/X3tr2tqK1rqUVXFiiaAUFNxTrEgmEPRCWTEL2TGY95/z+mGQgJIFJyGSWvJ+PBybnzJmTz8dk5j2f7f1RdF3XEUIIIU5iiHQBhBBCRCcJEEIIITokAUIIIUSHJEAIIYTokAQIIYQQHZIAIYQQokNhCxAPPPAAkydPZs6cOR0+rus6//Vf/8X06dOZO3cuJSUl4SqKEEKIbghbgLjmmmt49tlnO3188+bNlJWV8e677/KrX/2KX/7yl+EqihBCiG4IW4A477zz6NevX6ePb9iwgQULFqAoCueccw4NDQ1UVlaGqzhCCCG6yBSpH+xwOMjJyQke5+Tk4HA4yMrKOuXz3B4fOgoKEPgPKErgWGk5qSgtDykKBiXwuBBCiK6JWIDoKMNHKG/kDU4fjqrGLv+81oByYjAxKErguCWgKAal5RxtvhpazhsM4Q80mZnJVHWjfrFC6he74rlu0Dfq11URCxA5OTlUVFQEjysqKk7bejgTest/2gamrqWhOh5EwGgwYDAoGE/8Z1QwGmRimBAiPkQsQEybNo2XX36Z2bNn880335CcnBzWANETdEDXdDTAr6odXqMARoOCyWQIfDUaMJsMmIwSOIQQsSVsAeKee+7h888/p7a2lksuuYS7774bv98PwOLFi5k6dSqbNm1i+vTp2Gw2HnvssXAVpVfpgF/T8XvbBhAFgsHCbDJgMRuktSGEiGpKrKX7rqxp7tYYRDQyGhQsZiMWkwGr2YjBoPSJflCpX2yK57pB36hfV0Wsi0mAqum4PH5cnsCxxWTAZrfiVzXpkhJCRJwEiCji9WvUO73U1LsxGw0kWI1YzUYJFkKIiJAAEaV8qoavWaMRH2ajAavFSIJFgoUQovdIgIgBPlXD59JocvkwGRVsVhM2i6lX1mUIIfouCRAxxq/qNDb7aGr2YTEbSbSasFqMkS6WECIOSYCIUTrg8al4fCpGg0Jiggmb1YRB0ooIIXqIBIg4oGotrQqXj0SrCXuCWbqfhBBnTAJEHNF1cLr9NLv92BJMJEmgEEKcAQkQcUgHmt1+XB4/9gQz9gSTZLQVQnSZzJmMY7oOTS4fVfVuvL6Oc0cJIURnJED0AZqmU9voweXxR7ooQogYIgGij9CBeqcXp9sX6aIIIWKEBIg+prHZR2OzN9LFEELEAAkQfZDT7ZeWhBDitCRA9FGNzT4ZkxBCnJIEiD6s3unF7ZUgIYTomASIPq6+yYvPr0W6GEKIKCQBoo/TgbomD1psbSwohOgFEiAEqqZT3yQzm4QQbUmAEEAgM2yTS2Y2CSGOkwAhgppcPknJIYQIkgAh2qhzetE0GY8QQkiAECfRNJ16p4xHCCEkQIgOeHwqzbLSWog+TwKE6FBjs0/WRwjRx0mAEB0Kro+Q8Qgh+iwJEKJTqqZT1+RBl0V0QvRJEiDEKXn9Gg0yaC1EnyQBQpyWyyuL6IToi2IuQDz87Ba+3n1Muj16WZNLBq2F6GtMkS5AVx10NHLQ0chXu6tYMGUoackJkS5Sn9Hk8pGWbI10MYQQvSTmWhAmowLA7kP1/O/KYj7+9qhkIu0lHp+KR1JxCNFnxFyA+OWyyQzJSQbA59co+vQAz6zdzrF6V4RL1jc0NctYhBB9RcwFiAH97SybO5oFU4ZiNRsBOFDRyO9XfSutiV7gUzXZqlSIPiLmAgSAQVGYNCqb/7fobApz+wGBN66iTw/wXNEO6po8ES5hfHO6fDJJQIg+IKwBYvPmzcycOZPp06ezfPnydo83NjZyxx13MG/ePGbPns1rr73WpfunJln5/lUjueaSYcHWxL4jDTy1qpiv98hMp3Dxazouj4xFCBHvwhYgVFXlkUce4dlnn6WoqIh169axZ8+eNte88sorFBQUsGbNGl566SWeeOIJvN6uLcpSFIWJI7P48bVnM2RAYGzC7VV5deMeXn1/D26vdIeEgyTzEyL+hS1AFBcXk5+fT15eHhaLhdmzZ7Nhw4Y21yiKgtPpRNd1nE4n/fr1w2Tq3szbtGQry2aP5qrzB2M0BGY6fbOnmj+89i3llU1nXB/RVqAVIcFXiHgWtnUQDoeDnJyc4HF2djbFxcVtrlm6dCl33nknU6ZMwel08tvf/haD4fQxKz3d3ulj8y8rZMLoHJ5ds42K6mZqGj08vaaEBVMLuGLSYAyK0v1K9ZJT1S+amI0GMtMTu/y8zMzkMJQmesRz/eK5bhD/9euqsAWIjvr/lZPenD/66CNGjRrFiy++yMGDB7n55puZOHEiSUlJp7x3TY3zlI8nmg3cMW8MRZ8eYGtpJZqm8/r7e9i+9xjXXnoWiQnRuz4wPd1+2vpFE6/bGxz/CUVmZjJVVY1hLFFkxXP94rlu0Dfq11Vh62LKycmhoqIieOxwOMjKympzzeuvv86MGTNQFIX8/Hxyc3PZt29fj/x8i9nI1ZcMY/EVhcE3sNKDdfzh9WLKK+P3j6C3OSVHkxBxK2wBYty4cZSVlVFeXo7X66WoqIhp06a1uWbAgAF8+umnABw7doz9+/eTm5vbs+UYlsFdC8cxMCPQFVLX5GX5mu1sKamQWU49wOvX8PllRpMQ8ShsfS0mk4mHHnqIZcuWoaoqCxcupLCwkBUrVgCwePFifvjDH/LAAw8wd+5cdF3n3nvvJT09vcfLkpGSwO3zx/LWlgN8tt2Bqums+biM8somFkwZhtkUk8tBokaTy09acujdTEKI2KDoMfYxurKmGccZ9BN+vfsYb2zeh08NZCYdmJHI0hnDoybpX6yNQbTq3y8Bk/H0gbYv9PPGa/3iuW7QN+rXVX3uo/M5hf25Y8EY0luykh6pbuaPr29j35H6CJcstjXLlFch4k6fCxAAAzLs/OiacQzPSwUCb27PFZWyZXvFaZ4pOuPy+CUPlhBxpk8GCACb1cSNM0dwyfiBAGi6zpqPynjzw32ommyM01W6jiycEyLO9NkAAWAwKFx5/mC+O+2s4D4Tn++o5Pm3S+XNrhua3fL/TIh40qcDRKtzzurPD+aNISXRDMDeww386c1tVDe4I1yy2KJquuS+EiKOSIBokZuZxJ1XH18vcazezZ/e2EZZRUOESxZbpBUhRPyQAHGCfnYLP5g3htFD0oDWwesdfLuvOsIlix2BhXMyhiNEPJAAcRKL2ciS6cOZcvYAAPyqzt/f281HxUdl5XWIXNLNJERckADRAYOicNUF+cy9aAiKAjrw1pYDFH16QKZyhsDtldQbQsQDCRCnMHlMDkunD8fcskL4k20VvLpxD35VulBORdN0PBIkhIh5EiBOY/SQdJbNHUWiNZC2qnhvNS+u3ylvgKch3UxCxD4JECHIy0rm9vljSE2yALDncD3PrNtOk6S67pTHq6Jp0h0nRCyTABGizFQbt88fS3aaDYAjx5wsX1NCXZMnwiWLTjrImgghYpwEiC5onQabnxPIinis3s3Tq0s4VueKcMmik8sj3XBCxDIJEF1ks5q4edbIYKK/eqeXp9eUcORY7KXoDjefqsmAvhAxTAJEN1hMRr43YzjjhmUA4HT7eXbddtnKtAOSBlyI2CUBoptMRgPXTTuL80YG9tl2e1X+UrSDfUckNceJJA24ELFLAsQZMBgUFkwZykVjcwDw+jReeLuUXeV1ES5Z9NB1cEsrQoiYJAHiDCmKwqzJ+Vx67iAg0O/+0js7KT1YG+GSRQ+nJPATIiZJgOgBiqIw47w8ZpyXBwTSXr/y7i52lNVEuGTRQZWV1ULEJAkQPejScwdx1fmDgZYg8c/dlOyXIAHgdMuiQiFijQSIHjZl/EBmT84HAtuYrnhvN9skXbikARciBkmACIOLxg1g7oVDgECQ+PuGPdKSAJqlFSFETJEAESaTx+Yw76IhwPGWxPY+PibhlvxMQsQUCRBhdMGYnDYtiRXv7e7TA9c6snBOiFgiASLMJo/NYc6FgTEJVdP523u7+/Q6CelmEiJ2SIDoBReOHRAcuFY1nZff3cmew/URLlVk+GSwWoiYIQGil1w0bgAzJwXWSfhVnZfe2cn+o30zLYdLupmEiAkSIHrR1HMGcfmEXCDwSfqF9aV9MsGf2+tHl/xMQkQ9CRC9bNp3BjH1nIFAIHfTX98q5Wh130oVrumBGU1CiOgmAaKXtablaE3w5/aqPFe0g8o+tumQdDMJEf0kQERAa4K/iS2pwp1uP88V7aCmwR3hkvUer182ExIi2kmAiBBFUVhw8VDOLghsOtTg9PJc0Q7q+9Ae19LNJER0kwARQQaDwqLLChiVnwZATaOH3/3jK5r7SHps6WYSIrqFFCBuvfVW3n///S7PPNm8eTMzZ85k+vTpLF++vMNrPvvsM+bPn8/s2bP53ve+16X7xwOjwcD1lxcybGAKAEeqnDz/9o4+kR5b1XR8/vivpxCxKqQAcd111/HCCy9wxRVXsHz5cmprT78ZjqqqPPLIIzz77LMUFRWxbt069uzZ0+aahoYGHn74Yf70pz9RVFTE7373u+7VIsaZTQZumDGCvKwkAA5VOXnp3Z19YkFZs0cChBDRKqQAMWPGDJ5//nmeeeYZKisrmTNnDv/+7//Otm3bOn1OcXEx+fn55OXlYbFYmD17Nhs2bGhzzdq1a5k+fToDBwamfWZkZJxBVWKb1WLkpitHMjDTDsC+Iw28unFP3Ce388iaCCGilqk7TzKbzVitVu677z6mTJnC/fff3+4ah8NBTk5O8Dg7O5vi4uI215SVleH3+7nhhhtwOp3ceOONLFiw4LQ/Pz3d3p1iR7104P9ddy7/38v/4lidi5KyGt7eWs73rhyJoiiRLl6POfn3l5ySgM3arT/FqJSZmRzpIoRNPNcN4r9+XRXSq/Ldd9/l5Zdfprq6miVLllBUVITdbsfv9zNjxowOA0RHnwpPfpNTVZWSkhKef/553G43119/PePHj2fo0KGnLE9NTfwuLEtPt3PTlSN4enUJTS4fH39zBCNwZctOdbEuPd3e7vfnbHSTlmyNUIl6VmZmMlVV8bk6Pp7rBn2jfl0VUoBYtWoVt912G1OmTGn7ZJOJBx98sMPn5OTkUFFRETx2OBxkZWW1uyYtLY3ExEQSExOZOHEipaWlpw0Q8S4jJYGbZ43kmbXbcXtVNn9zBLvNxJSzB0a6aGHh9QX2iTAY4qeVJEQ8CGkM4umnn24XHFpNmzatw/Pjxo2jrKyM8vJyvF4vRUVF7a69/PLL+eKLL/D7/bhcLoqLiykoKOhiFeLTgAw7N145ApMx8Kb59paDfLmrKsKlCg+dQH4mIUR0CSlALFmyhPr64+mp6+rqWLp06SmfYzKZeOihh1i2bBmzZs3iqquuorCwkBUrVrBixQoACgoKmDJlCvPmzWPRokVce+21DB8+/AyqE1+G5KSw+IrhtH6wfn3TXnYePP0MsljkktlMQkQdRQ9hCsn8+fNZvXr1ac/1hsqaZhxx3E/YUR/9v3ZW8tqmfQCYjQZunTOKwdmxOZjWUf1aZaQkYDbF9trNeO7Hjue6Qd+oX1eF9GrUNI3m5ubgsdPpRFXlE19vmTAiiysnBQapfarGC+t3Ulkbf8n9XNLNJERUCSlAzJkzh1tuuYXVq1ezevVqbr31VubNmxfusokTTBk/gIvHDQACKSr++lb85W1ye1VZEyFEFAlpFtPtt99OVlYWGzduRNd1rr/++pDWK4ieoygKV14wmCaXj6/3HKPe6eWvb5dy+7wxcbOGQNN0vD4Nq8UY6aIIIejCQrmrr76aq6++OpxlEadhUBSumToMp9vH7kP1VNa6ePGdndwya1TM9923cnn9EiCEiBIhBYjq6mpeeuklysvL8fuP9xP31dxJkWQyGlhyxXCeXbedw8ecHKho5B8bd7PkiuFxsY7A45U1EUJEi5ACxN13301BQQGTJ0/GaJRPd5FmtRi56aqR/Hn1NmoaPGwvq2XtJ2XMu2hIzKfkaF0TkZhgjnRRhOjzQgoQDQ0N/OpXvwp3WUQXJNnM3DxrFH9+cxtOt5/PtjvoZ7dw6bmDIl20M+byqBIghIgCIXVcFxYW4nA4wl0W0UUZKQncdNVILC3jD+9uLedfOysjXKoz51O1PpHqXIhoF3ILYt68eZx77rlYrceTqskYROTlZiaxZPpwXly/E03XeWPzPpJsZkYMTot00c5Is8dPP5Ml0sUQok8LKUDMmTOHOXPmhLssopuG56VyzdRhrPpgL5oOK97bzbK5o8nNTIp00brN7fWTnGjGEONjKkLEspAChExvjX7fGZ5Jg9PLu1vL8foDq63vmD+GjJSESBetW3Qd3B6VxIT4WOMhRCwKaQyirKyMxYsXB7OxlpSU8Pvf/z6sBRNdN/WcgVwwOhsAp8vH82+V0uTyRbhU3efySOoNISIppADxy1/+kjvvvJPk5ECyp1GjRrF+/fqwFkx0naIozLlwCKOHBMYfqhvcvLi+FK8vNvNmBQarY7PsQsSDkAJEY2Mjl1xySXCOvcFgwGyWaYjRyGBQuG5aIfkt2V4PVTlZsWE3aozubd0sacCFiJiQAoTRaMTn8wUDhMPhwGCIj9QO8chsMnDDzOH07xcYf9h5sI41H+2PyUR4bq8fLQbLLUQ8CHnDoLvuuova2lp+//vfs2TJEm655ZZwl02cgcQEMzfPGkmyLdDS21paycYvD0e4VF0XGKyWsQghIiGkKSILFiwgNzeX999/H5fLxRNPPMHEiRPDXTZxhtKSAwvplq8twevT2PCvQ/SzW5g4Muv0T44iTrek3hAiEkKeQzhx4kQJCjFoYH87S6cP54W3Awvp3vxwH8mJsbWQTtV03F4/CRaZ8ipEbwrpFbdw4cIOk8CtWrWqxwskel5hbioLpw5jZQwvpGt2S4AQoreF9Iq77777gt97PB6KiorIyoqtboq+7tzhmdTH8EI6rz8w5dVskmzCQvSWkALEpEmT2hxffPHFMkgdg6aeM5B6p5fPtjuCC+lunz+GJFts9O873X5SkyRACNFbujVXtampifLy8p4uiwgzRVGYG8ML6TxeFVWTLK9C9JYuj0FomsahQ4e4+eabw1owER6tC+n+UrSdg46m4EK6780YgTHKd3HTCYxFJCdKllchekOXxyCMRiO5ublkZ2eHrVAivMwmAzfOHMHTa0qoqnOz82Adqz/cx9WXDIv6HelcXpXkxEiXQoi+oVtjECL2JSaY+f5Vo/jz6m00Nvv4YmcVKXYLV0zMi3TRTknTdBmsFqKXhBQgLrjggg4/Weq6jqIofPrppz1eMBF+aclWvn/VSJav2Y7Hp7Lxy8MkJ1o4f3R0tw5dXgkQQvSGkALE4sWLqaur47rrrkPXdV577TWys7OZNWtWuMsnwmxAhp3vzRjO82+Xomo6az7eT3KimdFD0iNdtE55vCpIN5MQYRfSLKatW7fyn//5n4wcOZJRo0bx4IMPsmnTJgYNGsSgQYPCXUYRZgWD+nHtpQVAIPfR3zfs5kBFY4RL1TlV02XPaiF6QUgBorKykpqamuBxTU0NVVVVYSuU6H3jz+rP7Mn5APhVnRfWl+KoaY5wqTrn9koCPyHCLaQupptuuon58+dz2WWXAbBp0yZuv/32sBZM9L6Lxg2gwenlw+KjuL0qz78dWEiXmmSNdNHa8chsJiHCLqQAsXTpUiZMmMDWrVvRdZ2lS5cyYsSIcJdNRMDM8wfT5PLx1e5j1Du9PP92KT+YOybq9ob2azp+VcNklH1JhAiXkF/1ubm5qKrKmDFjwlkeEWEGReGaqcNocvnYfaieyloXL75Tyi2zR2GJsplDbq9Kkk0ChBDhEtKra9OmTcyePZu7774bgG+//ZY77rgjrAUTkWM0GFgyfTi5mXYADjqaWPHe7qhLcyHjEEKEV0gB4qmnnmLVqlWkpKQAMG7cOA4ePBjWgonIspqN3HTVyDbblr6xeV9UbVvqVwPdTEKI8Ai5fZ6Zmdnm2GKRfDjxzp5g5uZZo0ixB37XX+46xtufHYyqIOGS7UiFCJuQAoTdbufYsWPB1dSfffYZycnJp33e5s2bmTlzJtOnT2f58uWdXldcXMyoUaNYv359iMUWvaV1tbXNGhh/+Kj4KJu/ORLhUh3n8vijKmAJEU9CChA/+9nPuO222zh06BA33HAD9957b5sEfh1RVZVHHnmEZ599lqKiItatW8eePXs6vO43v/kNF198cfdqIMIuJz2RG2eOxNwyY+idz8vZWloZ4VIFaHpgsFoI0fNCmsU0fvx4XnzxRb788ksAzj333OB4RGeKi4vJz88nLy+Q/G327Nls2LCBs846q811L730EjNnzuTbb7/tTvlFL8nPSWbJ9EJeemdXcG9rm8XI2GEZkS4aLo8fmzW6puEKEQ9O+6pSVZXvfve7vPbaa0ydOjXkGzscDnJycoLH2dnZFBcXt7vmvffe44UXXuhSgEhPt4d8bSyK1vpNTrdjsph4bk0Jug6vvr+H/hl2Rg/tWpAIR/1S0xIxm6Jjymtm5um7X2NVPNcN4r9+XXXaAGE0GklLS8Pj8WC1hr6itqN+4ZMzwj766KPce++9GI1dm19fU+Ps0vWxJD3dHtX1K8hJZs5FQ1j7cRl+VedPrxVz6+xRDM4O7YUVrvq5mz2kRMFGQpmZyVRVRW8eqzMRz3WDvlG/rgqpXT5kyBCWLl3KzJkzSUw8nt9g6dKlnT4nJyeHioqK4LHD4SArK6vNNdu2beOee+4BoLa2lk2bNmEymbjiiiu6VAnRuyaPycHl8fPeF4fw+TWef7uU2+aOZkBG5Fo+bo+fZJs56jc8EiKWhBQgnE4nhYWF7Nu3L+Qbjxs3jrKyMsrLy8nOzqaoqIj/+Z//aXPNxo0bg9/ff//9XHrppRIcYsRl5w7C7VH56NtA3qa/vlXKD+aOpn+qLSLlaR2slrEIIXrOKV9Nv/71r7n//vt5/PHH+fjjj7noootCv7HJxEMPPcSyZctQVZWFCxdSWFjIihUrgMAeEyJ2KYrCVRcMxuX186+dVTS5fPylaAc/mDeGtOTIJPdzun0SIIToQYp+iknkV199NW+88Ua77yOpsqYZRxz3E0b7GMTJNE3nHxt38+2+QDr4jJQEbps3utPxgHDXLy3JitUSuZxR8dyPHc91g75Rv6465bSPE2OHLEYSHTEYFBZddhYjBqcCUN3g5rmiHTjdvoiUp8kVmZ8rRDw6ZYDwer3s3buXPXv2tPm+9Z8QACajgSVXDGfYwMDamMpaF38t2hGRNBg+VcPjk4VzQvSEU3YxTZs2rfMnKgobNmwIS6FORbqYopfHp/LXt3Zw0NEEQF5WErfMGtWmy6c36mcxGUhPSQjrz+hMPHdTxHPdoG/Ur6tOOaJ34iwjIU7Hajby/atG8pd1Ozh8zEl5ZRMvrC/l+1eNxGLuvXEBr1/D51cxR9n+FULEmuhYeiriRoLFxM2zRpKTHlgvU1bRyIvv7MTr791unyaXZHkV4kzJnEDR4xITzNwyexTPrN1OVZ2LfUcaePmdXdww88y2qd19qI4vSiupbfSQlmxl4sgsCnNTO7zW41NlS9IYtW1/NR8VH6WqzkVmqo2Lzx7A2C6mcxE9Q149IiySbGZunTMquOHQnsP1vPzuTnzdbEnsPlTHO5+XU93gQdOhusHDO5+Xs/tQXafPaZa9ImLOtv3VvLZpH45aF5oOjloXr23ax7b91ZEuWp8kAUKETUqihWVzRpPREiR2H6rnz69/i8/f9V3gvugkvXhn5yGQ5VWT6dkx5aPio106L8JLAoQIqxR7IEikpwRWV5fsq25pSXQtSNQ2erp0HkDXwe2RKa+xpKrO1cl5dy+XRIAECNEL+p0UJHYfquelLg5cd5a+43RpPZo9snAulmR2kssrMzUy05b7OgkQolekJlm5bc5ostICbwB7Dtfz4vqdeENc1DZxZFaXzrfyq3rIP0NE3sVnD+jSeRFeEiBEr+mXZOWeJROCA9f7jjTw/NuleELYMrQwN5WZk/LISLFiUCAjxcrMSXmdzmI6kQxWx46xQzNYOHUY2Wk2DIpCdpqNhVOHySymCDnlSupoJCupY1t6up2yQ7X8Zd32YL9yXlYS379qZNgysSpA/9QEjIbwfx6K59W48Vw36Bv16yppQYhel5Jo4ba5Y4KL6corm3h23fawJdrTAadbWhFCdJUECBERSTYzy+aMZlBmYBe6o9XNPLN2O/VOb1h+nsvtR9NiqrEsRMRJgBARk5hgatnPOgkITHFcvqaE6oaen9IYaEXIjCYhukIChIioQO6mUZw1qB8QWNewfHUJFTXNPf6zmj3SihCiKyRAiIizmo3ceOUIRg9JA6DR5eOZtSUcdPTsgKGuy4wmIbpCAoSICiajgcVXDOc7w/sD4PKo/GXdDnYerO3Rn+N0+yT9hhAhkgAhoobRoHDN1AIuGpcDBHaHe+mdnXy1q6rHfoauQ7PMaBIiJBIgRFQxKAqzLshn5qQ8ADQdVn6wl81fH+mxfdGdLh9+tesJA4XoayRAiKijKApTzxnEwqnDMCiBc+s/P8jaj8t6ZJBZB+qbwjOdVoh4IgFCRK0JI7L43owRmE2BP9Mt2x288s9dPbI7nU/VwrYwT4h4IQFCRLWR+WksmzMae0IgDceOA7X8Zd0OGpvPvAXgdPm6tTeFEH2FBAgR9fKykrhjwdjgxkPllU38eXUJjjNcK6ED9U5Pj41tCBFvJECImJCRksAd88eQnx1IOFbb6OHPq0tOueVoKPyqLnmahOiEBAgRM+wJgX2uzzkrsFbC41N54e1SPi2pOKNWgHQ1CdExCRAippiMBhZdVsDlE3KBwDTYtR+Xsfqj/d2euipdTUJ0TAKEiDmKonD5hFyuv/wsTMbAPNjPd1Ty3Fs7uj0zSbqahGiYcBYGAAAd9ElEQVRPAoSIWWcX9Of2eWNIsVsAKDvayP+98S2Hq5q6dT/pahKiLQkQIqYNykziR1ePJS8rkDK8rsnL02tK+LIb6Tl0oLbJI6ushWghAULEvOREC7fNHc15I7OAQHfRqg/2dmtcQtN0aholSAgBEiBEnDAZDVx9yTCunjIUY0t+js+2O1i+poTaRk+X7iVBQogACRAirpw3KpsfzBtNv5ZxiUNVTv7wejGlXUwbLkFCiDAHiM2bNzNz5kymT5/O8uXL2z2+Zs0a5s6dy9y5c7n++uspLS0NZ3FEH5GXlcxdC8dRmBvYpc7lUXlx/U7e3nKgS2/4EiREXxe2AKGqKo888gjPPvssRUVFrFu3jj179rS5Jjc3l5dffpm1a9dy55138otf/CJcxRF9jD3BzE1XjuTyCbm0JITlw+KjLF9TQk0X9rzWNJ2aBrcECdEnhS1AFBcXk5+fT15eHhaLhdmzZ7Nhw4Y213znO9+hX7/Ap7xzzjmHioqKcBVH9EEGQ2C9xM2zR5FsMwOBLqffv/YtX+85FvJ9NB1qGtz4eiCLrBCxxBSuGzscDnJycoLH2dnZFBcXd3r9qlWruOSSS0K6d3q6/YzLF82kfj1rUrqdkcP680LRdkr2VePxqby6cQ/7KxpZPGMEiQnmkO6jA9ZEC8mJZhRF6fS6zMzkHip59InnukH816+rwhYgOkpb0NmLasuWLaxatYq//e1vId27psZ5RmWLZunpdqlfmCy+/Cw+yUrinc8Pomo6W7c72HWglkWXFTBsYL+Q7lFT48RsNNAvyYLJ2L4BnpmZTFVVY08XPSrEc92gb9Svq8LWxZSTk9Omy8jhcJCVldXuutLSUh588EH+7//+j7S0tHAVRwgMisLFZw/gh1ePJSvNBkC908tf1u1g3SdlIW9E5FM1quvdNLtlwyER38IWIMaNG0dZWRnl5eV4vV6KioqYNm1am2uOHDnC3XffzX//938zdOjQcBVFiDYGZNj50dXjuHBsoAtUBz7ZVsHvV33LgYrQPkHqQEOzTwawRVwLWxeTyWTioYceYtmyZaiqysKFCyksLGTFihUALF68mD/+8Y/U1dXx8MMPA2A0Gnn99dfDVSQhgswmA3MuHMLI/DRe37SXuiYv1Q1ulq8pYfLYHKafl4fVbDztfbx+jeoGN8k2C4kJYXs5CRERih5jOY4ra5pxxHE/oYxB9L6S/dW8veUgNSesuE6ymejfLwG/qpOWbGXiyCwKc1NPeR+r2UjBkAxqqrueLHDb/mo+Kj5KVZ2LzFQbF589gLFDM7p8n3BY92kZH3x1GKfbjz3BxKXnDmLO5CGRLlaPkzGI9uQjj+jTdh+q4/2vjpBgNZGuQH2TF1XTaXL5aXI1YbMa8ak673xeDnDKIOHxqVTWNOPz+LFZQ39pbdtfzWub9gWPHbWu4HGkg8S6T8tY93EZEJhk0tTsCx7HY5AQbUmqDdGnfVFaGfw+wWIiM80WzOUEgVXYlbXNOF0+tu5wnPZ+mq5T7/RSXe/G4w1t0Puj4qNdOt+bPvjqcJfOi/giAUL0aScn8jMoCgYDGI0Ep7HqemC2087yesorQ+uC8KkatU0eahrceHynDhRVda5Ozoe+4jtcOtuAydnNjZlEbJEAIfq0tGRru3NGgwGz0UhmagIpdguty3d8fo0/vVnCqg/20tjsDen+Xr9GbaOHY/UuXJ6Od6zLTLV1cj4htEqEUZKt40WE9k7Oi/giAUL0aRNHtl+bk5hgwp5gQlEUkmxmslJtJFiOz2j6clcVT/7jGzZ9fTjkHej8aqDrqbLORZPLh6Ydnxty8dkDOnxOZ+d706XnDurSeRFfZJBa9Gmtg85flFZS2+hpmbGU1+ZcRpqNWZPzMRgU1n1chqPWhcen8s7n5Xy23cGMSYM5uyADwynSb7TSNJ0mlw+ny0eC1USi1RgciA7MYnKTmZoQNbOYWgeiP/jqMM1uP3abOW5nMYn2ZJprlInGaaA9Kdbrp2o6n2938N6/DrXpMhqUaWfmeYOZdPbALtfPZFRItJqwWowYDdHbqO8L00DjvX5dJQEiysT6G+jpxEv9XB4/H3x1mE+2VaCe0F00Ij+NaecOCu6R3VUmo4LVbMRqNmIJYaFeb+oLb6DxXr+uki4mIbrBZjVx1QX5nD86m39+Uc43e6oB2Hmglp0Hahk5OI3LJwxiUGbXAoVf1fGrfpxuPwaDQoLZiNVixGIynDKDrBDhIAFCiDOQnpLAddMKuWT8QN7dWs7Og3UAlB6spfRgLaPy07jsO4PI7WKggMB4RbPHT7PHj0IgPUggWBgxm6K3K0rEDwkQQvSAARl2brpyJDXNPl7fuJt9RxoA2HGglh0HainM7cfUcwYxdEByt1oCOoEps16/BviOty7MRixmaV2I8JAAIUQPOis3lWVzRrP/aAMb/nUoGCh2H6pn96F6BmcnMeXsgYzKT8Ng6P6bemetC7PRgMlkCGlGlRCnIwFCiDAYOiAlGCg2fX2EXeWBrqeDjiZe+ecuMlISuGhcDt8ZnnnGg9FtWxcBBoOC2WjAajZgli4p0U0SIIQIo6EDUhg6IIUjx5xs+vow2/bXoOtQ3eBmzcdlvLu1nIkjs7hgdDbpKT23clrTdDya2pLmw4dBAYs5MH5hMRs63A1PiJNJgBCiFwzsb2fxFcOpaXDz8bYK/lVaidev4faqfFR8lI+LjzJicBqTRmcxPDf1jLqfOqLp4PaquFsSCBoUMJuMmIwKJqMBo0HBaFQwKIqMZ4ggCRBC9KL0lATmXjiEKybk8kVpJZ+WVFDX5EXn+Myn1CQLE0dmMWF4Jv2S2ueK6gmaHkhP7jkp554CGI1KcKZUa/Do6YAlYkPMLZTTNJ2qqka0lmLruo6mAzptzukEsnDquh74ygnft1yntTygBZ4UvC6S4mUhWWekfm1pmk7pwVq2lDjYc7i+zWMKcFZuPyaMyGJUflpExxEUAgut6mqdGFtbHC2Bo/X7WG95yEK59mKuBWFo+aM0EL4/xuNBpTWgnBh0jp+DluBEINi0CUi6jqbpaICu6REPPCI6GQwKo4ekM3pIOtX1bj7f4eBfu6podvvROT77KcFiZOywDM4t7E9+TnKvz1LSAb/aMhDeSYJCg0JLwDBgNCqYWr4PvGaR7qsYFHMtCCAmo3xr4AgEkkBgCQSR44FI03TS0+1UHWtqG2iO3+T4/U76JlZ+idKCOD2/qrHjQC3/2lnF7kN1nPwK7We3MK4gg7MLMhjU395rb7o9UTdFaQ0UdBjklJbHFAi2Tk78Gs4gIy2I9iRARJkz/SNtDUCtAUfVdDS95esJ/1p/6a1flZZ/wQNAafnmVK/H1m66YOuq5aad/VFJgOiaeqeXr3dX8dXuY1TWtt9YKD3Zypih6Ywdls6gzKSwtiyi5XcXCBQcb50oBINH64ew1r93ON6lrCgnBpq291QUhf79k6ipdgYCVGugUuKn5SMBIg7E06eYjrrq+vdPoqqqsU13XeDawAu7NajpJ7zIY+kPNFxvorquc7S6ma93H+PbfdXUO9tvWJRitzAqP41R+WkMG5jS41NZoyVAhMup6qcAikHBQCAYKUogyCjK8YATONfSOjIoURdc+sQYhIgdrZ/CTmibBObid3Fh2ImtotZuuZNbRqqmtQSWnq1DtFAUhYH97Qzsb+fKCwZT7miieG812/ZX09gcmIrU4PTy2XYHn213YDEZOCu3HyPyUhmelxq22VB9hU5gLFED0Lr2R3ZicDn+mujgupMeCLZkWm/S2jI/8Y+85TltnqnQ7pyiKGR2qdQBEiBE1DO0tvVDcPJYT+vMtNaWiapq+FUdVddjdvKAQVHIz0kmPyeZ2Rfmc6iyiZL9NZTsr6GmZY9tr19je1kt28tqAchOs1GYm8pZuf0YkpMcdanE41mb4BJjf3ESIERcaf2EFuost9ZuLa2lRdJm0sAJj/lP6NOOJgZFYXB2MoOzk7ny/MFU1rkobUkQWF7ZFPyw6ah14ah18dG3RzEaFPKykxg2IIVhA1PIy0qWVByiQxIgRJ9mUBQMRgVC/ECtahqaRkvQaGmNqBp+Vetqz0OPUxSF7LREstMSmXrOIJrdfvYcrmNXeR27y+tpdPla6qBTdrSRsqONbPzyMCajQm5mEvk5yQzJCQQbm1XeGoQECCG6xGgw0Dr2az0pqmi6TlqqDdXjO94C0XW8XhV/BKJHYoKJswv6c3ZBf3Rdx1HrYu/hevYcqqesorElT1Ngk6KyikbKKhrZ1PLczFQbg7OTGJyVRG5WEllpib1efhF5EiCE6CEGRcFiNrb/9J0IPr+Gy+vH7VUj0lWlKAo56YnkpCdy0bgBqJrOkWNO9h2pp+xoIwccjcE8TQBVdS6q6lz8a2cVAGajgcEDkslKtTGoZbA8M9WGUVJwxDUJEEL0ArPJgNlkIdmm4/aqON0+/Grk+qSMBoW8rCTyspKYek6gteOoaeZARSMHHU0crGykpsETvN6nauw9VM/eQ8fTgZiMCllpiQxITyQnI5Hs9ESy02wk2cxRM7VTnBkJEEL0IkVRsFlN2KwmvD4Vl8eP26dGfHquQVEYkGFnQIadC8YEzjW5fByuaqK8sonDVU6OVDuDU2oh0DV15JiTI8farh1ITDCRlWYjK9VGZvBfAv2SrLKRUYyRACFEhLSuCUnWdTxeFa9PxesPrOeIBkk2MyMGpzFicBoAaWmJlJXXcuSYk8PHnFTUNHO0upnaRk+b5zW7/cFB8BOZjAr9+9nISEkgo5+VjJQE0lMSSE+xkmK3SndVFJIAIUSEGU5oVUBgppTXp+H1BTb8iZJ4gaIo9Euy0i/Jyqgh6cHzbq8fR42LippmHLXNOGoC4xdNrra5xP2qTkVNMxU1ze3ubVAUUpMtpCZZSUu2Br/2Swqc62e3yCZHESABQogoYzQYsFkNwYDh86t4WgKGt5NMqpGUYDEFF+6dqNntDw52H6t3UVXn5li9m5oGd7tWkqbr1DR42ox7nMyeYKKf3UKK3UJyYuBrSqKZ5EQLSa1fbSaMBgkkPUUChBBRLrCntBFsZjRdx+fT8PpVfH4Nn1+L2rW5iQkdBw5N02lo9nKs3k1to4eahkDQqGvyUtvoadfyaOV0+3G6/Rypbt8COZHNaiLJZibJZsZuM2FPCHyfmGDCnmAiMcFMotVEYkLgn9lokEH1TkiAECKGGBQFq8WI1XJ8DYZfDQQKf0saEZ+qReWq71YGg0JqUqAbqSNev0pdk5f6Jg/1TV7qmjzUO700OL3UO73UN3mDazg64vL4cXkCrZdQmIyBLr6kRAsWkwGbxYTNaiTBYiLBaiTB0vK9xRj8ZzWbsJoNWC2Bfb7jdcc9CRBCxDiT0dCuf17Tj+ed8qtaSx6q2MhBZTEZyUoNzILqjNen0tjso97ppcnlpbHZR2OzjyZX239Ol++0g/5+VQ8+v7vMJgNWsxGLufWrEYvJ0O6r2WQIbud68j+TseX7lq+mlu9NRkNwv/DeFtYAsXnzZh599FE0TWPRokX84Ac/aPO4rus8+uijbNq0iYSEBH79618zZsyYcBZJiLi1bX81HxUfparORWaqjYvPHsDYoRn8pWg7W3dU4lM1zEYD543K4ryRWcFrM/olcP7obIbnprL9QC2f73BQXe8mPcXKhBFZHKpq4vPtDpq9KokWI5NGZ3PZubmdlmP3oTq+KK2kttFDWrKViSOzANqdK8xNDfn5hbmpvP/VoUA5PH4SraZTlqP1HjUN7sCgen4amak29hyup/RALY3NPixmA+nJCZhNBlxePz6/TmOzF5fHj8erdimItnb3EVqjpVuMBqXlw4ASDBomoyGwc1/rcctufsaW3fxMrd8bDYwbkd3lnxm2/SBUVWXmzJn89a9/JTs7m2uvvZYnn3ySs846K3jNpk2beOmll3jmmWf45ptvePTRR1m5cuVp7x0v+yV0JJ72g+iI1C88tu2v5rVN+9qdtyeY2NGS0bWVTmB8IPOkT+gTRmQGV063qm1043T5gvsatL5dzDx/MDPOGxzImqsd3wO+9GANb2052GZdh9vrRwGslrafR2dOymsXJHYfquOdz8vb1SMn3ca3e6vbnb/0O4PaBYnO7jFmaBol+2vbnW8tx4n7Qewqr2X9Z+Vtd4DUdMYOSyc1yYrbG5hh5vGpeFq+9/m14LHXH9ie1dtyPhqs/Z/5XX5O2FoQxcXF5Ofnk5eXB8Ds2bPZsGFDmwCxYcMGFixYgKIonHPOOTQ0NFBZWUlWVla4iiVEXPqo+GiH50sPtH9DBHC5/e3OffDVYZITLW3OOV1+dD2wn8GJA7mfbKvgmksK2t3jmz3V7bq76ppUdF3HbjO3Of/t3mouGJ0dTMmODsV7j7VbD6G3XNuRrTsqmT15SJutd7/efazdjnGt155cBoCvdlUxdmhGy9hCYGzn693HAp/ET8oK3OD0Mu+ioR2WpTOaruPzH5+F5jshcLQet/mnqi1jSnpwfMmnaqgnHPtVDf+J6etbjv0t16mtiSTPcCwqbAHC4XCQk5MTPM7Ozqa4uPiU1+Tk5OBwOE4bILqzM1IskfrFtkjUr7bJ22HKbk3vYCuNljfkk693uv2kpySc9Pzj23We+LXZ7e+wnh2VQ1V1UGi3B0WDy8fAAW1bEA3N/jYD8K38mo7F1P6826tSkJ/R5lyT209CB9loq+rd9O9gXKPZqzJ8WH8AMvoFHnd61A4z2jZ7VUYUdGfrndgUtgDRUc/VyVPJQrlGCHF6T/5kaqSLAJx5OXqiHtFyj3gQthUlOTk5VFRUBI87ahmcfE1FRYV0LwkhRJQIW4AYN24cZWVllJeX4/V6KSoqYtq0aW2umTZtGm+++Sa6rvP111+TnJwsAUIIIaJE2LqYTCYTDz30EMuWLUNVVRYuXEhhYSErVqwAYPHixUydOpVNmzYxffp0bDYbjz32WLiKI4QQoovCNs1VCCFEbJOsVkIIITokAUIIIUSHojoXk8fjYenSpXi93uDK7B//+MfU1dXx05/+lMOHDzNo0CD+93//l379+kW6uN3SOj6TnZ3N008/HVd1mzZtGna7HYPBgNFo5PXXX4+r+jU0NPDggw+ya9cuFEXhscceY+jQoXFRv3379vHTn/40eFxeXs6Pf/xjFixYEBf1e/7551m5ciWKojB8+HAef/xxXC5XXNQN4IUXXmDlypXous6iRYv4/ve/363XXlS3ICwWCy+88AJr1qzhzTff5MMPP+Trr79m+fLlTJ48mXfffZfJkyezfPnySBe121588UUKCo6vSI2nukHgD3X16tW8/vrrQHzV79FHH2XKlCmsX7+e1atXU1BQEDf1GzZsGKtXrw7+7mw2G9OnT4+L+jkcDl588UVee+011q1bh6qqFBUVxUXdAHbt2sXKlStZuXIlq1ev5oMPPqCsrKxb9YvqAKEoCna7HQC/34/f70dRlGCKDoAFCxbw3nvvRbKY3VZRUcEHH3zAtddeGzwXL3XrTLzUr6mpia1btwZ/dxaLhZSUlLip34k+/fRT8vLyGDRoUNzUT1VV3G43fr8ft9tNVlZW3NRt7969jB8/HpvNhslk4rzzzuOf//xnt+oX1QECAr/I+fPnc+GFF3LhhRcyfvx4qqurg+slsrKyqKmpiXApu+exxx7j3/7t3zCcsANWvNSt1a233so111zDP/7xDyB+6ldeXk56ejoPPPAACxYs4Oc//znNzc1xU78TFRUVMWfOHCA+fn/Z2dnccsstXHbZZVx88cUkJSVx8cUXx0XdAIYPH84XX3xBbW0tLpeLzZs3U1FR0a36RX2AMBqNrF69mk2bNlFcXMyuXbsiXaQe8f7775Oens7YsWMjXZSwWbFiBW+88QbPPPMMr7zyClu3bo10kXqM3+9n+/btLF68mDfffBObzRazXRKn4vV62bhxI1deeWWki9Jj6uvr2bBhAxs2bODDDz/E5XKxevXqSBerxxQUFLBs2TJuueUWli1bxogRIzAa2+exCkXUB4hWKSkpnH/++Xz44YdkZGRQWVkJQGVlJenp6ad5dvT58ssv2bhxI9OmTeOee+5hy5Yt3HvvvXFRt1bZ2YH88xkZGUyfPp3i4uK4qV9OTg45OTmMHz8egCuvvJLt27fHTf1abd68mTFjxtC/f0syuzio3yeffEJubi7p6emYzWZmzJjBV199FRd1a7Vo0SLeeOMNXnnlFVJTU8nPz+9W/aI6QNTU1NDQ0ACA2+3mk08+YdiwYcEUHQBvvvkml19+eSSL2S0/+9nP2Lx5Mxs3buTJJ5/kggsu4De/+U1c1A2gubmZpqam4Pcff/wxhYWFcVO/zMxMcnJy2LcvsAfDp59+SkFBQdzUr1VRURGzZ88OHsdD/QYOHMg333yDy+VC1/W4/N1VVwfSox85coR3332XOXPmdKt+Ub2SurS0lPvvvx9VDeSTv/LKK7nrrruora3lJz/5CUePHmXAgAH87ne/IzW1492pYsFnn33Gc889x9NPPx03dSsvL+dHP/oREBhHmjNnDnfeeWfc1A9gx44d/PznP8fn85GXl8fjjz+OpmlxUz+Xy8Wll17Ke++9R3JyILV3vPz+nnrqKd566y1MJhOjRo3i0Ucfxel0xkXdAJYsWUJdXR0mk4kHHniAyZMnd+t3F9UBQgghROREdReTEEKIyJEAIYQQokMSIIQQQnRIAoQQQogOSYAQQgjRoajO5irEqSxatAiv14vP56OsrIzCwkIARo8ezeOPPx7h0oWmpKSE8vLyuFqpLOKHTHMVMe/QoUMsXLiQzz77LNJFacfv92Mydf45bOXKlXzyySf89re/7fF7C3Gm5K9LxKVVq1bx97//HVVVSUlJ4eGHH2bIkCGsXLmS9evXY7fb2bVrFwMGDOA//uM/eOKJJygvL2f8+PE88cQTKIrCvffei81m4+DBg1RUVHD++efzi1/8ArPZTGNjI4899hi7d+/G4/Fw4YUXct9992EwGFi8eDGTJk3iq6++IjExkaeeeiq4SNDj8TB+/HgefvhhGhoa+OMf/4jT6WT+/Pmcf/75LF26lCVLlvDxxx8DcODAgeDxgQMHWLx4Mddddx1btmzhmmuuYf78+Tz55JN88cUXeL1eRo0axS9/+UtsNluEfwMiLuhCxLjy8nJ90qRJweMtW7bot99+u+7xeHRd1/UNGzboS5cu1XVd11999VV90qRJekVFha7run7LLbfoCxYs0BsbG3Wv16vPmjVL37Jli67ruv6zn/1Mnz9/vu50OnWv16vfeOON+t/+9jdd13X9vvvu09euXavruq6rqqr/+Mc/1letWqXruq5ff/31+g9/+EPd7/cHH6+rqwt+f8899+ivvvpqsDw/+clPgmUvKyvTL7zwwg6Py8rK9OHDh+vr168PPv7UU0/pTz/9dPD48ccf13/3u9+d2f9QIVpIC0LEnY0bN7J9+3YWLVoEgK7rOJ3O4OMTJkwIJhIcPXo0brebpKQkAEaMGMHBgwc5//zzAZg1axaJiYlAIIf+Bx98wOLFi3n//fcpKSnhmWeeAQK5wgYPHhz8GXPnzg1m0NQ0jeXLl/PRRx+haRp1dXXd3qksMTGRmTNntqmry+WiqKgICGRfHTNmTLfuLcTJJECIuKPrOt/97ne56667OnzcarUGvzcYDO2O/X5/p/dVFAUIvOk//fTTDBw4sMNrW4MKwOrVqykuLuZvf/sbdrudP/zhDxw9erTD5xmNRjRNCx57PJ5O79tapl/96lecd955Hd5PiDMh01xF3GnNWulwOIBAssBt27Z1615vv/02LpcLn8/H2rVrgy2LadOmsXz5clRVBQKZh8vLyzu8R2NjI2lpadjtdurr64Of9gHsdjuNjY3B46ysLNxud/Be69atO21dn3vuuWAgaWpqYu/evd2qqxAnkwAh4s4FF1zAXXfdxe233868efOYO3cuH3zwQbfuNWHCBO68807mzJlDXl5ecIvRX/ziF2iaxvz585k7dy633XYbVVVVHd7j6quvpq6ujjlz5nDPPfe0+bR/0UUX0djYyLx583jsscewWCzcf//93HTTTdxwww2YzeZTlu+OO+6goKCAa6+9lrlz57J06VL279/frboKcTKZ5ipEJ+69914mTJjA4sWLI10UISJCWhBCCCE6JC0IIYQQHZIWhBBCiA5JgBBCCNEhCRBCCCE6JAFCCCFEhyRACCGE6ND/D2bIkFEziuh0AAAAAElFTkSuQmCC\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\"." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "celltoolbar": "Hide code", + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- 2.18.1