Upload New File

parent 3f2dce47
{
"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": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.7.16 (default, Jan 17 2023, 16:06:28) [MSC v.1916 64 bit (AMD64)]\n",
"uname_result(system='Windows', node='CLT-PHAC2105-P', release='10', version='10.0.19041', machine='AMD64', processor='Intel64 Family 6 Model 140 Stepping 1, GenuineIntel')\n",
"3.7.16\n",
"IPython 7.31.1\n",
"IPython.core.release 7.31.1\n",
"PIL 9.3.0\n",
"PIL.Image 9.3.0\n",
"PIL._deprecate 9.3.0\n",
"PIL._version 9.3.0\n",
"_csv 1.0\n",
"_ctypes 1.1.0\n",
"decimal 1.70\n",
"_pydev_bundle.fsnotify 0.1.5\n",
"_pydevd_frame_eval.vendored.bytecode 0.13.0.dev\n",
"argparse 1.1\n",
"backcall 0.2.0\n",
"bottleneck 1.3.5\n",
"cffi 1.15.1\n",
"colorama 0.4.6\n",
"csv 1.0\n",
"ctypes 1.1.0\n",
"cycler 0.10.0\n",
"dateutil 2.8.2\n",
"debugpy 1.5.1\n",
"decimal 1.70\n",
"decorator 5.1.1\n",
"defusedxml 0.7.1\n",
"entrypoints 0.4\n",
"http.server 0.6\n",
"ipykernel 6.15.2\n",
"ipykernel._version 6.15.2\n",
"ipython_genutils 0.2.0\n",
"ipython_genutils._version 0.2.0\n",
"ipywidgets 7.6.5\n",
"ipywidgets._version 7.6.5\n",
"jedi 0.18.1\n",
"json 2.0.9\n",
"jupyter_client 7.4.9\n",
"jupyter_client._version 7.4.9\n",
"jupyter_core 4.11.2\n",
"jupyter_core.version 4.11.2\n",
"kiwisolver 1.4.4\n",
"kiwisolver._cext 1.4.4\n",
"logging 0.5.1.2\n",
"matplotlib 3.5.3\n",
"matplotlib_inline 0.1.6\n",
"mkl 2.4.0\n",
"numexpr 2.8.4\n",
"numpy 1.21.5\n",
"numpy.core 1.21.5\n",
"numpy.core._multiarray_umath 3.1\n",
"numpy.lib 1.21.5\n",
"numpy.linalg._umath_linalg 0.1.5\n",
"packaging 22.0\n",
"packaging.__about__ 22.0\n",
"pandas 1.3.5\n",
"parso 0.8.3\n",
"patsy 0.5.6\n",
"patsy.version 0.5.6\n",
"pickleshare 0.7.5\n",
"pkg_resources._vendor.appdirs 1.4.3\n",
"pkg_resources._vendor.more_itertools 8.12.0\n",
"pkg_resources._vendor.packaging 21.3\n",
"pkg_resources._vendor.packaging.__about__ 21.3\n",
"pkg_resources._vendor.pyparsing 3.0.9\n",
"pkg_resources._vendor.appdirs 1.4.3\n",
"pkg_resources._vendor.more_itertools 8.12.0\n",
"pkg_resources._vendor.packaging 21.3\n",
"pkg_resources._vendor.pyparsing 3.0.9\n",
"platform 1.0.8\n",
"prompt_toolkit 3.0.36\n",
"psutil 5.9.0\n",
"pydevd 2.6.0\n",
"pygments 2.11.2\n",
"pyparsing 3.0.9\n",
"pytz 2022.7\n",
"re 2.2.1\n",
"scipy 1.7.3\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.12.2\n",
"seaborn.external.appdirs 1.4.4\n",
"seaborn.external.husl 2.1.0\n",
"six 1.16.0\n",
"socketserver 0.4\n",
"statsmodels 0.13.5\n",
"statsmodels.__init__ 0.13.5\n",
"statsmodels._version 0.13.5\n",
"statsmodels.api 0.13.5\n",
"statsmodels.tools.web 0.13.5\n",
"traitlets 5.7.1\n",
"traitlets._version 5.7.1\n",
"urllib.request 3.7\n",
"wcwidth 0.2.5\n",
"xmlrpc.client 3.7\n",
"zlib 1.0\n",
"zmq 23.2.0\n",
"zmq.sugar 23.2.0\n",
"zmq.sugar.version 23.2.0\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",
" print(platform.python_version())\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": [
"Main libraries :\n",
"- statsmodels 0.13.5\n",
"- numpy 1.21.5\n",
"- pandas 1.3.5\n",
"- matplotlib 3.5.3\n",
"- seaborn 0.12.2\n",
"\n",
"Python version : 3.7.16\n",
"\n",
"platform : Jupyter\n",
"\n",
"OS : Windows x64 (build 19045)\n",
"\n",
"## Loading and inspecting data\n",
"Let's start by reading data."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Date</th>\n",
" <th>Count</th>\n",
" <th>Temperature</th>\n",
" <th>Pressure</th>\n",
" <th>Malfunction</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>4/12/81</td>\n",
" <td>6</td>\n",
" <td>66</td>\n",
" <td>50</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>11/12/81</td>\n",
" <td>6</td>\n",
" <td>70</td>\n",
" <td>50</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>3/22/82</td>\n",
" <td>6</td>\n",
" <td>69</td>\n",
" <td>50</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>11/11/82</td>\n",
" <td>6</td>\n",
" <td>68</td>\n",
" <td>50</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>4/04/83</td>\n",
" <td>6</td>\n",
" <td>67</td>\n",
" <td>50</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>6/18/82</td>\n",
" <td>6</td>\n",
" <td>72</td>\n",
" <td>50</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>8/30/83</td>\n",
" <td>6</td>\n",
" <td>73</td>\n",
" <td>100</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>11/28/83</td>\n",
" <td>6</td>\n",
" <td>70</td>\n",
" <td>100</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>2/03/84</td>\n",
" <td>6</td>\n",
" <td>57</td>\n",
" <td>200</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>4/06/84</td>\n",
" <td>6</td>\n",
" <td>63</td>\n",
" <td>200</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>8/30/84</td>\n",
" <td>6</td>\n",
" <td>70</td>\n",
" <td>200</td>\n",
" <td>1</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>10/05/84</td>\n",
" <td>6</td>\n",
" <td>78</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>11/08/84</td>\n",
" <td>6</td>\n",
" <td>67</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>1/24/85</td>\n",
" <td>6</td>\n",
" <td>53</td>\n",
" <td>200</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>4/12/85</td>\n",
" <td>6</td>\n",
" <td>67</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>15</th>\n",
" <td>4/29/85</td>\n",
" <td>6</td>\n",
" <td>75</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>16</th>\n",
" <td>6/17/85</td>\n",
" <td>6</td>\n",
" <td>70</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>17</th>\n",
" <td>7/2903/85</td>\n",
" <td>6</td>\n",
" <td>81</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>8/27/85</td>\n",
" <td>6</td>\n",
" <td>76</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>10/03/85</td>\n",
" <td>6</td>\n",
" <td>79</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>10/30/85</td>\n",
" <td>6</td>\n",
" <td>75</td>\n",
" <td>200</td>\n",
" <td>2</td>\n",
" </tr>\n",
" <tr>\n",
" <th>21</th>\n",
" <td>11/26/85</td>\n",
" <td>6</td>\n",
" <td>76</td>\n",
" <td>200</td>\n",
" <td>0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>22</th>\n",
" <td>1/12/86</td>\n",
" <td>6</td>\n",
" <td>58</td>\n",
" <td>200</td>\n",
" <td>1</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Date Count Temperature Pressure Malfunction\n",
"0 4/12/81 6 66 50 0\n",
"1 11/12/81 6 70 50 1\n",
"2 3/22/82 6 69 50 0\n",
"3 11/11/82 6 68 50 0\n",
"4 4/04/83 6 67 50 0\n",
"5 6/18/82 6 72 50 0\n",
"6 8/30/83 6 73 100 0\n",
"7 11/28/83 6 70 100 0\n",
"8 2/03/84 6 57 200 1\n",
"9 4/06/84 6 63 200 1\n",
"10 8/30/84 6 70 200 1\n",
"11 10/05/84 6 78 200 0\n",
"12 11/08/84 6 67 200 0\n",
"13 1/24/85 6 53 200 2\n",
"14 4/12/85 6 67 200 0\n",
"15 4/29/85 6 75 200 0\n",
"16 6/17/85 6 70 200 0\n",
"17 7/2903/85 6 81 200 0\n",
"18 8/27/85 6 76 200 0\n",
"19 10/03/85 6 79 200 0\n",
"20 10/30/85 6 75 200 2\n",
"21 11/26/85 6 76 200 0\n",
"22 1/12/86 6 58 200 1"
]
},
"execution_count": 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\") # Modify \"blob to raw\"\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": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwQ0lEQVR4nO3de1yUdf7//+fIUVyhREVMRExNS7cSt8JTWiuu2sHcLdM1PGAbN01UtNJ1Nw9ZdlhZO3moVDLLjx2sj/WjlMrU0G3TxE7+zDxhCkuiBcpHGOH6/uE6t0ZQYK6RGd497rcbt7re1/uaeV0vLsen13XNjMOyLEsAAACGaODrAgAAALyJcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjOLTcLNp0ybdeuutatmypRwOh955551qt9m4caPi4+MVGhqqtm3bavHixRe/UAAAUG/4NNycPHlSV199tZ577rkazd+/f78GDhyoXr16aceOHfrrX/+q1NRUvfXWWxe5UgAAUF84/OWLMx0Oh95++20NHjz4vHMeeughrV27Vrt27XKNpaSkaOfOndq6dWsdVAkAAPxdoK8LqI2tW7cqMTHRbax///5aunSpnE6ngoKCKm1TWlqq0tJS13JFRYWOHTumyMhIORyOi14zAACwz7IsFRcXq2XLlmrQ4MIXnupVuMnPz1dUVJTbWFRUlE6fPq2jR48qOjq60jbz5s3T7Nmz66pEAABwER06dEitWrW64Jx6FW4kVTrbcvaq2vnOwkyfPl1paWmu5Z9//lmtW7fW/v371bhx44tXqA84nU5t2LBBffv2rfIsFqpHD+2hf/bRQ3von33+2sPi4mLFxcXV6O/uehVuWrRoofz8fLexgoICBQYGKjIyssptQkJCFBISUmm8SZMmCg8Pvyh1+orT6VRYWJgiIyP96oCsT+ihPfTPPnpoD/2zz197eLaWmtxSUq8+5yYhIUFZWVluY+vXr1e3bt386hcAAAB8x6fh5sSJE8rJyVFOTo6kM2/1zsnJUW5urqQzl5SSkpJc81NSUnTw4EGlpaVp165dWrZsmZYuXaqpU6f6onwAAOCHfHpZatu2berbt69r+ey9MSNHjlRGRoby8vJcQUeS4uLilJmZqcmTJ+v5559Xy5Yt9cwzz+iPf/xjndcOAAD8k0/DTZ8+fXShj9nJyMioNHbjjTfqiy++uIhVAQCA+qxe3XMDAABQHcINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFF8Hm4WLlyouLg4hYaGKj4+Xps3b77g/FdffVVXX321wsLCFB0drdGjR6uwsLCOqgUAAP7Op+Fm9erVmjRpkmbMmKEdO3aoV69eGjBggHJzc6uc/+mnnyopKUnJycn65ptv9MYbb+jzzz/X2LFj67hyAADgr3wabtLT05WcnKyxY8eqU6dOWrBggWJiYrRo0aIq5//rX/9SmzZtlJqaqri4OPXs2VP33Xeftm3bVseVAwAAfxXoqycuKyvT9u3bNW3aNLfxxMREbdmypcptunfvrhkzZigzM1MDBgxQQUGB3nzzTQ0aNOi8z1NaWqrS0lLXclFRkSTJ6XTK6XR6YU/8x9n9MW2/6hI9tIf+2UcP7aF/9vlrD2tTj8OyLOsi1nJeR44c0WWXXabs7Gx1797dNf7YY4/p5Zdf1u7du6vc7s0339To0aN16tQpnT59WrfddpvefPNNBQUFVTl/1qxZmj17dqXx1157TWFhYd7ZGQAAcFGVlJRo+PDh+vnnnxUeHn7BuT47c3OWw+FwW7Ysq9LYWd9++61SU1P18MMPq3///srLy9MDDzyglJQULV26tMptpk+frrS0NNdyUVGRYmJilJiYWG1z6hun06msrCz169fvvGEPF0YP7aF/9tFDe+ifff7aw7NXXmrCZ+GmadOmCggIUH5+vtt4QUGBoqKiqtxm3rx56tGjhx544AFJ0m9/+1s1atRIvXr10ty5cxUdHV1pm5CQEIWEhFQaDwoK8qtfmjeZvG91hR7aQ//so4f20D/7/K2HtanFZzcUBwcHKz4+XllZWW7jWVlZbpepfqmkpEQNGriXHBAQIOnMGR8AAACfvlsqLS1NL730kpYtW6Zdu3Zp8uTJys3NVUpKiqQzl5SSkpJc82+99VatWbNGixYt0r59+5Sdna3U1FRdd911atmypa92AwAA+BGf3nMzdOhQFRYWas6cOcrLy1Pnzp2VmZmp2NhYSVJeXp7bZ96MGjVKxcXFeu655zRlyhRdcskluummm/TEE0/4ahcAAICf8fkNxePGjdO4ceOqXJeRkVFpbMKECZowYcJFrgoAANRXPv/6BQAAAG8i3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYxefhZuHChYqLi1NoaKji4+O1efPmC84vLS3VjBkzFBsbq5CQEF1++eVatmxZHVULAAD8XaAvn3z16tWaNGmSFi5cqB49emjJkiUaMGCAvv32W7Vu3brKbe666y795z//0dKlS9WuXTsVFBTo9OnTdVw5AADwVz4NN+np6UpOTtbYsWMlSQsWLNC6deu0aNEizZs3r9L8Dz74QBs3btS+ffvUpEkTSVKbNm3qsmQAAODnfBZuysrKtH37dk2bNs1tPDExUVu2bKlym7Vr16pbt2568skn9corr6hRo0a67bbb9Mgjj6hhw4ZVblNaWqrS0lLXclFRkSTJ6XTK6XR6aW/8w9n9MW2/6hI9tIf+2UcP7aF/9vlrD2tTj8/CzdGjR1VeXq6oqCi38aioKOXn51e5zb59+/Tpp58qNDRUb7/9to4ePapx48bp2LFj573vZt68eZo9e3al8fXr1yssLMz+jvihrKwsX5dQ79FDe+ifffTQHvpnn7/1sKSkpMZzfXpZSpIcDofbsmVZlcbOqqiokMPh0KuvvqqIiAhJZy5t/elPf9Lzzz9f5dmb6dOnKy0tzbVcVFSkmJgYJSYmKjw83It74ntOp1NZWVnq16+fgoKCfF1OvUQP7aF/9tFDe+ifff7aw7NXXmrCZ+GmadOmCggIqHSWpqCgoNLZnLOio6N12WWXuYKNJHXq1EmWZemHH35Q+/btK20TEhKikJCQSuNBQUF+9UvzJpP3ra7QQ3von3300B76Z5+/9bA2tfjsreDBwcGKj4+vdNorKytL3bt3r3KbHj166MiRIzpx4oRr7LvvvlODBg3UqlWri1ovAACoHzwKN/v37/fKk6elpemll17SsmXLtGvXLk2ePFm5ublKSUmRdOaSUlJSkmv+8OHDFRkZqdGjR+vbb7/Vpk2b9MADD2jMmDHnvaEYAAD8unh0Wapdu3bq3bu3kpOT9ac//UmhoaEePfnQoUNVWFioOXPmKC8vT507d1ZmZqZiY2MlSXl5ecrNzXXN/81vfqOsrCxNmDBB3bp1U2RkpO666y7NnTvXo+cHAADm8Sjc7Ny5U8uWLdOUKVN0//33a+jQoUpOTtZ1111X68caN26cxo0bV+W6jIyMSmMdO3b0uzu4AQCA//DoslTnzp2Vnp6uw4cPa/ny5crPz1fPnj111VVXKT09XT/++KO36wQAAKgRWzcUBwYG6o477tDrr7+uJ554Qnv37tXUqVPVqlUrJSUlKS8vz1t1AgAA1IitcLNt2zaNGzdO0dHRSk9P19SpU7V37159/PHHOnz4sG6//XZv1QkAAFAjHt1zk56eruXLl2v37t0aOHCgVqxYoYEDB6pBgzNZKS4uTkuWLFHHjh29WiwAAEB1PAo3ixYt0pgxYzR69Gi1aNGiyjmtW7fW0qVLbRUHAABQWx6Fmz179lQ7Jzg4WCNHjvTk4QEAADzm0T03y5cv1xtvvFFp/I033tDLL79suygAAABPeRRuHn/8cTVt2rTSePPmzfXYY4/ZLgoAAMBTHoWbgwcPKi4urtJ4bGys2ycKAwAA1DWPwk3z5s315ZdfVhrfuXOnIiMjbRcFAADgKY/Czd13363U1FRt2LBB5eXlKi8v18cff6yJEyfq7rvv9naNAAAANebRu6Xmzp2rgwcP6uabb1Zg4JmHqKioUFJSEvfcAAAAn/Io3AQHB2v16tV65JFHtHPnTjVs2FBdunRxfZs3AACAr3gUbs7q0KGDOnTo4K1aAAAAbPMo3JSXlysjI0MfffSRCgoKVFFR4bb+448/9kpxAAAAteVRuJk4caIyMjI0aNAgde7cWQ6Hw9t1AQAAeMSjcPM///M/ev311zVw4EBv1wMAAGCLR28FDw4OVrt27bxdCwAAgG0ehZspU6bo6aeflmVZ3q4HAADAFo8uS3366afasGGD3n//fV111VUKCgpyW79mzRqvFAcAAFBbHoWbSy65RHfccYe3awEAALDNo3CzfPlyb9cBAADgFR7dcyNJp0+f1ocffqglS5aouLhYknTkyBGdOHHCa8UBAADUlkdnbg4ePKg//OEPys3NVWlpqfr166fGjRvrySef1KlTp7R48WJv1wkAAFAjHp25mThxorp166bjx4+rYcOGrvE77rhDH330kdeKAwAAqC2P3y2VnZ2t4OBgt/HY2FgdPnzYK4UBAAB4wqMzNxUVFSovL680/sMPP6hx48a2iwIAAPCUR+GmX79+WrBggWvZ4XDoxIkTmjlzJl/JAAAAfMqjy1L//Oc/1bdvX1155ZU6deqUhg8frj179qhp06ZatWqVt2sEAACoMY/CTcuWLZWTk6NVq1bpiy++UEVFhZKTk/XnP//Z7QZjAACAuuZRuJGkhg0basyYMRozZow36wEAALDFo3CzYsWKC65PSkryqBgAAAC7PAo3EydOdFt2Op0qKSlRcHCwwsLCCDcAAMBnPHq31PHjx91+Tpw4od27d6tnz57cUAwAAHzK4++WOlf79u31+OOPVzqrAwAAUJe8Fm4kKSAgQEeOHPHmQwIAANSKR/fcrF271m3Zsizl5eXpueeeU48ePbxSGAAAgCc8CjeDBw92W3Y4HGrWrJluuukmzZ8/3xt1AQAAeMSjcFNRUeHtOgAAALzCq/fcAAAA+JpHZ27S0tJqPDc9Pd2TpwAAAPCIR+Fmx44d+uKLL3T69GldccUVkqTvvvtOAQEB6tq1q2uew+HwTpUAAAA15FG4ufXWW9W4cWO9/PLLuvTSSyWd+WC/0aNHq1evXpoyZYpXiwQAAKgpj+65mT9/vubNm+cKNpJ06aWXau7cubxbCgAA+JRH4aaoqEj/+c9/Ko0XFBSouLjYdlEAAACe8ijc3HHHHRo9erTefPNN/fDDD/rhhx/05ptvKjk5WUOGDPF2jQAAADXm0T03ixcv1tSpUzVixAg5nc4zDxQYqOTkZD311FNeLRAAAKA2PAo3YWFhWrhwoZ566int3btXlmWpXbt2atSokbfrAwAAqBVbH+KXl5envLw8dejQQY0aNZJlWd6qCwAAwCMehZvCwkLdfPPN6tChgwYOHKi8vDxJ0tixY3kbOAAA8CmPws3kyZMVFBSk3NxchYWFucaHDh2qDz74wGvFAQAA1JZH99ysX79e69atU6tWrdzG27dvr4MHD3qlMAAAAE94dObm5MmTbmdszjp69KhCQkJsFwUAAOApj8JN7969tWLFCteyw+FQRUWFnnrqKfXt29drxQEAANSWR5elnnrqKfXp00fbtm1TWVmZHnzwQX3zzTc6duyYsrOzvV0jAABAjXl05ubKK6/Ul19+qeuuu079+vXTyZMnNWTIEO3YsUOXX365t2sEAACosVqfuXE6nUpMTNSSJUs0e/bsi1ETAACAx2p95iYoKEhff/21HA7HxagHAADAFo8uSyUlJWnp0qXergUAAMA2j24oLisr00svvaSsrCx169at0ndKpaene6U4AACA2qpVuNm3b5/atGmjr7/+Wl27dpUkfffdd25zuFwFAAB8qVbhpn379srLy9OGDRsknfm6hWeeeUZRUVEXpTgAAIDaqtU9N+d+6/f777+vkydPerUgAAAAOzy6ofisc8OOJxYuXKi4uDiFhoYqPj5emzdvrtF22dnZCgwM1DXXXGO7BgAAYI5ahRuHw1Hpnho799isXr1akyZN0owZM7Rjxw716tVLAwYMUG5u7gW3+/nnn5WUlKSbb77Z4+cGAABmqtU9N5ZladSoUa4vxzx16pRSUlIqvVtqzZo1NXq89PR0JScna+zYsZKkBQsWaN26dVq0aJHmzZt33u3uu+8+DR8+XAEBAXrnnXdqswsAAMBwtQo3I0eOdFseMWKEx09cVlam7du3a9q0aW7jiYmJ2rJly3m3W758ufbu3auVK1dq7ty51T5PaWmpSktLXctFRUWSznzSstPp9LB6/3R2f0zbr7pED+2hf/bRQ3von33+2sPa1FOrcLN8+fJaF3M+R48eVXl5eaV3WkVFRSk/P7/Kbfbs2aNp06Zp8+bNCgysWenz5s2r8msi1q9fr7CwsNoXXg9kZWX5uoR6jx7aQ//so4f20D/7/K2HJSUlNZ7r0Yf4edO59+xYllXlfTzl5eUaPny4Zs+erQ4dOtT48adPn660tDTXclFRkWJiYpSYmKjw8HDPC/dDTqdTWVlZ6tevn4KCgnxdTr1ED+2hf/bRQ3von33+2sOzV15qwmfhpmnTpgoICKh0lqagoKDKz80pLi7Wtm3btGPHDt1///2SpIqKClmWpcDAQK1fv1433XRTpe1CQkJc9wj9UlBQkF/90rzJ5H2rK/TQHvpnHz20h/7Z5289rE0ttt4KbkdwcLDi4+MrnfbKyspS9+7dK80PDw/XV199pZycHNdPSkqKrrjiCuXk5Oj666+vq9IBAIAf8+llqbS0NN1zzz3q1q2bEhIS9MILLyg3N1cpKSmSzlxSOnz4sFasWKEGDRqoc+fObts3b95coaGhlcYBAMCvl0/DzdChQ1VYWKg5c+YoLy9PnTt3VmZmpmJjYyVJeXl51X7mDQAAwC/5/IbicePGady4cVWuy8jIuOC2s2bN0qxZs7xfFAAAqLd8ds8NAADAxUC4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwo0X7fvxhDbsLtD+oyd9XQoA/Cod+O/r78HCEh9XAl8K9HUBJvippEypq3K0ac+PrrHe7Zvp2WHXKiIsyIeVAcCvw9nX4c/2FejJ66RBz27W9W2b8zr8K8WZGy9IXZWj7O+Puo1lf39UE1bt8FFFAPDrwuswfolwY9O+H09o054fVW5ZbuPllqVNe37kEhUAXGS8DuNchBubDh678HXdA4X8oQKAi4nXYZyLcGNTbJOwC65vE9mojioBgF8nXodxLsKNTW2b/Ua92zdTgMPhNh7gcKh3+2aKa8ofKgC4mHgdxrkIN17w7LBr1aNdU7exHu2a6tlh1/qoIgD4deF1GL/EW8G9ICIsSCuSr9P+oyd1oPCk2kQ24l8KAFCHzr4Of5//s7757BP9fxN6qV2LCF+XBR8h3HhRXFNCDQD4UmxkmL7573/x68VlKQAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRfB5uFi5cqLi4OIWGhio+Pl6bN28+79w1a9aoX79+atasmcLDw5WQkKB169bVYbUAAMDf+TTcrF69WpMmTdKMGTO0Y8cO9erVSwMGDFBubm6V8zdt2qR+/fopMzNT27dvV9++fXXrrbdqx44ddVw5AADwVz4NN+np6UpOTtbYsWPVqVMnLViwQDExMVq0aFGV8xcsWKAHH3xQv/vd79S+fXs99thjat++vd599906rhwAAPirQF89cVlZmbZv365p06a5jScmJmrLli01eoyKigoVFxerSZMm551TWlqq0tJS13JRUZEkyel0yul0elC5/zq7P6btV12ih/bQP/vooT30zz5/7WFt6vFZuDl69KjKy8sVFRXlNh4VFaX8/PwaPcb8+fN18uRJ3XXXXeedM2/ePM2ePbvS+Pr16xUWFla7ouuJrKwsX5dQ79FDe+ifffTQHvpnn7/1sKSkpMZzfRZuznI4HG7LlmVVGqvKqlWrNGvWLP3v//6vmjdvft5506dPV1pammu5qKhIMTExSkxMVHh4uOeF+yGn06msrCz169dPQUFBvi6nXqKH9tA/++ihPfTPPn/t4dkrLzXhs3DTtGlTBQQEVDpLU1BQUOlszrlWr16t5ORkvfHGG/r9739/wbkhISEKCQmpNB4UFORXvzRvMnnf6go9tIf+2UcP7aF/9vlbD2tTi89uKA4ODlZ8fHyl015ZWVnq3r37ebdbtWqVRo0apddee02DBg262GUCAIB6xqeXpdLS0nTPPfeoW7duSkhI0AsvvKDc3FylpKRIOnNJ6fDhw1qxYoWkM8EmKSlJTz/9tG644QbXWZ+GDRsqIiLCZ/sBAAD8h0/DzdChQ1VYWKg5c+YoLy9PnTt3VmZmpmJjYyVJeXl5bp95s2TJEp0+fVrjx4/X+PHjXeMjR45URkZGXZcPAAD8kM9vKB43bpzGjRtX5bpzA8snn3xy8QsCAAD1ms+/fgEAAMCbCDfwmn0/ntCG3QXaf/SkT7YHPTRB9p4fJUlb9x71cSVA/eXzy1Ko/34qKVPqqhxt+u+LsiT1bt9Mzw67VhFh1b91z+72oIcmOFh4UoOfz1ZJaZmevE6695XtCgsJ1trxPRUTaeYHjgIXC2duYFvqqhxlf+/+r8zs749qwqqafaGp3e1BD00w+PlsHS9x/3j54yVO3fb8pz6qCKi/CDewZd+PJ7Rpz48qtyy38XLL0qY9P1Z7ecTu9qCHJti4u6BSsDnreIlTm39xRg5A9Qg3sOXgsQt/18eBwgv/xWp3e9BDE+T88NMF13+Re7xuCgEMQbiBLbFNLnwvQJvIRhd1e9BDE1zT6pILru/a+tK6KQQwBOEGtrRt9hv1bt9MAed82WmAw6He7ZsprumF/2K1uz3ooQluvKK5Lj3Pjd+XhgWpV/tmdVwRUL8RbmDbs8OuVY92Td3GerRrqmeHXVsn24MemmDt+J6VAs6lYUFaO76njyoC6i/eCg7bIsKCtCL5Ou0/elIHCk+qTWSjWp0tsLs96KEJYiLDtOPhRG36//N0fPe/9eI98erdMdrXZQH1EuEGXhPX1N5fqHa3Bz00QcLlTZW5+8x/AXiGy1IAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABjF5+Fm4cKFiouLU2hoqOLj47V58+YLzt+4caPi4+MVGhqqtm3bavHixXVUKYCLbd+PJ7Rhd4H2Hz3p0far/52rSat36I1th3xWg93ts/f8KEnauveoR9vb5ev9t+vAf5/3YGGJR9tv3F2gpz/6Tpv/+3vwhfr+O/D180tSoM+eWdLq1as1adIkLVy4UD169NCSJUs0YMAAffvtt2rdunWl+fv379fAgQN17733auXKlcrOzta4cePUrFkz/fGPf/TBHgDwhp9KypS6KkebfvEXSu/2zfTssGsVERZU7fZf/fCT7li4RacrLEnSOzuOaPqar7R2fA9deVlEndRgd/uDhSc1+PlslZSW6cnrpHtf2a6wkGCtHd9TMZFhNdoHO3y9/3adff7P9hXoyeukQc9u1vVtm9e6/8dLnK6xS8OC6qz/kjm/A189/y/59MxNenq6kpOTNXbsWHXq1EkLFixQTEyMFi1aVOX8xYsXq3Xr1lqwYIE6deqksWPHasyYMfrHP/5Rx5UD8KbUVTnK/t79TEX290c1YdWOGm3/y2Bz1ukKS7c9n11nNdjd/ty/WCXpeIlTtz3/aY22t8vX+29Xfe+/xO/Am3x25qasrEzbt2/XtGnT3MYTExO1ZcuWKrfZunWrEhMT3cb69++vpUuXyul0KiiocjIsLS1VaWmpa/nnn3+WJB07dkxOp7PS/PrM6XSqpKREhYWFVfYC1aOH9njSv9zCk9qy64AcqvyCtGXXCe3YE63WTc7/L+e1OYelshPnfTF7+eOvdMvVLS9qDXa3/9feoyou+kmBkgIrLJWUVCjQ2UDlFQ4VF0nvb/tO18VFXnAf7PD1/tvl9vwN3PtX2/6fqy76L/nX78AXf45rori4WJJkWVY1M89M8onDhw9bkqzs7Gy38UcffdTq0KFDldu0b9/eevTRR93GsrOzLUnWkSNHqtxm5syZliR++OGHH3744ceAn0OHDlWbMXx6z40kORwOt2XLsiqNVTe/qvGzpk+frrS0NNdyRUWFjh07psjIyAs+T31UVFSkmJgYHTp0SOHh4b4up16ih/bQP/vooT30zz5/7aFlWSouLlbLlhc+Eyv58LJU06ZNFRAQoPz8fLfxgoICRUVFVblNixYtqpwfGBioyMiqTxmGhIQoJCTEbeySSy7xvPB6IDw83K8OyPqIHtpD/+yjh/bQP/v8sYcRERE1muezG4qDg4MVHx+vrKwst/GsrCx17969ym0SEhIqzV+/fr26devG/REAAECSj98tlZaWppdeeknLli3Trl27NHnyZOXm5iolJUXSmUtKSUlJrvkpKSk6ePCg0tLStGvXLi1btkxLly7V1KlTfbULAADAz/j0npuhQ4eqsLBQc+bMUV5enjp37qzMzEzFxsZKkvLy8pSbm+uaHxcXp8zMTE2ePFnPP/+8WrZsqWeeeYbPuPmvkJAQzZw5s9JlONQcPbSH/tlHD+2hf/aZ0EOHZdXkPVUAAAD1g8+/fgEAAMCbCDcAAMAohBsAAGAUwg0AADAK4aaemTVrlhwOh9tPixYtXOtHjRpVaf0NN9zgw4r90+HDhzVixAhFRkYqLCxM11xzjbZv3+5ab1mWZs2apZYtW6phw4bq06ePvvnmGx9W7F+q6x/H4YW1adOmUn8cDofGjx8vieOvOtX1j+OveqdPn9bf/vY3xcXFqWHDhmrbtq3mzJmjiooK15z6fBz6/OsXUHtXXXWVPvzwQ9dyQECA2/o//OEPWr58uWs5ODi4zmqrD44fP64ePXqob9++ev/999W8eXPt3bvX7ZOrn3zySaWnpysjI0MdOnTQ3Llz1a9fP+3evVuNGzf2XfF+oCb9kzgOL+Tzzz9XeXm5a/nrr79Wv379dOedd0ri+KtOdf2TOP6q88QTT2jx4sV6+eWXddVVV2nbtm0aPXq0IiIiNHHiREn1/Dis9tun4FdmzpxpXX311eddP3LkSOv222+vs3rqo4ceesjq2bPneddXVFRYLVq0sB5//HHX2KlTp6yIiAhr8eLFdVGiX6uuf5bFcVhbEydOtC6//HKroqKC488Dv+yfZXH81cSgQYOsMWPGuI0NGTLEGjFihGVZ9f91kMtS9dCePXvUsmVLxcXF6e6779a+ffvc1n/yySdq3ry5OnTooHvvvVcFBQU+qtQ/rV27Vt26ddOdd96p5s2b69prr9WLL77oWr9//37l5+crMTHRNRYSEqIbb7xRW7Zs8UXJfqW6/p3FcVgzZWVlWrlypcaMGSOHw8HxV0vn9u8sjr8L69mzpz766CN99913kqSdO3fq008/1cCBAyXV/9dBwk09c/3112vFihVat26dXnzxReXn56t79+4qLCyUJA0YMECvvvqqPv74Y82fP1+ff/65brrpJpWWlvq4cv+xb98+LVq0SO3bt9e6deuUkpKi1NRUrVixQpJcX8567he4RkVFVfri1l+j6voncRzWxjvvvKOffvpJo0aNksTxV1vn9k/i+KuJhx56SMOGDVPHjh0VFBSka6+9VpMmTdKwYcMkGXAc+vrUEew5ceKEFRUVZc2fP7/K9UeOHLGCgoKst956q44r819BQUFWQkKC29iECROsG264wbIsy8rOzrYkWUeOHHGbM3bsWKt///51Vqe/qq5/VeE4PL/ExETrlltucS1z/NXOuf2rCsdfZatWrbJatWplrVq1yvryyy+tFStWWE2aNLEyMjIsy6r/xyFnbuq5Ro0aqUuXLtqzZ0+V66OjoxUbG3ve9b9G0dHRuvLKK93GOnXq5Poes7PvPjv3XycFBQWV/hXza1Rd/863DcdhZQcPHtSHH36osWPHusY4/mquqv5VheOvsgceeEDTpk3T3XffrS5duuiee+7R5MmTNW/ePEn1/zgk3NRzpaWl2rVrl6Kjo6tcX1hYqEOHDp13/a9Rjx49tHv3brex7777zvWFrXFxcWrRooWysrJc68vKyrRx40Z17969Tmv1R9X1ryoch1Vbvny5mjdvrkGDBrnGOP5qrqr+VYXjr7KSkhI1aOAeAQICAlxvBa/3x6GvTx2hdqZMmWJ98skn1r59+6x//etf1i233GI1btzYOnDggFVcXGxNmTLF2rJli7V//35rw4YNVkJCgnXZZZdZRUVFvi7db/z73/+2AgMDrUcffdTas2eP9eqrr1phYWHWypUrXXMef/xxKyIiwlqzZo311VdfWcOGDbOio6Ppo1V9/zgOa6a8vNxq3bq19dBDD1Vax/FXvfP1j+OvZkaOHGlddtll1nvvvWft37/fWrNmjdW0aVPrwQcfdM2pz8ch4aaeGTp0qBUdHW0FBQVZLVu2tIYMGWJ98803lmVZVklJiZWYmGg1a9bMCgoKslq3bm2NHDnSys3N9XHV/ufdd9+1OnfubIWEhFgdO3a0XnjhBbf1FRUV1syZM60WLVpYISEhVu/eva2vvvrKR9X6nwv1j+OwZtatW2dJsnbv3l1pHcdf9c7XP46/mikqKrImTpxotW7d2goNDbXatm1rzZgxwyotLXXNqc/HocOyLMvXZ48AAAC8hXtuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AXJDD4bjgz6hRo3xdotf16dNHkyZN8nUZADwU6OsCAPi3vLw81/+vXr1aDz/8sNsXZzZs2NAXZXnE6XQqKCjI2OcDcAZnbgBcUIsWLVw/ERERcjgcbmObNm1SfHy8QkND1bZtW82ePVunT592be9wOLRkyRLdcsstCgsLU6dOnbR161Z9//336tOnjxo1aqSEhATt3bvXtc2sWbN0zTXXaMmSJYqJiVFYWJjuvPNO/fTTT261LV++XJ06dVJoaKg6duyohQsXutYdOHBADodDr7/+uvr06aPQ0FCtXLlShYWFGjZsmFq1aqWwsDB16dJFq1atcm03atQobdy4UU8//bTr7NSBAweUkZGhSy65xO3533nnHTkcjkp1L1u2TG3btlVISIgsy9LPP/+sv/zlL2revLnCw8N10003aefOnV76DQE4F+EGgMfWrVunESNGKDU1Vd9++62WLFmijIwMPfroo27zHnnkESUlJSknJ0cdO3bU8OHDdd9992n69Onatm2bJOn+++932+b777/X66+/rnfffVcffPCBcnJyNH78eNf6F198UTNmzNCjjz6qXbt26bHHHtPf//53vfzyy26P89BDDyk1NVW7du1S//79derUKcXHx+u9997T119/rb/85S+655579Nlnn0mSnn76aSUkJOjee+9VXl6e8vLyFBMTU+OenK37rbfeUk5OjiRp0KBBys/PV2ZmprZv366uXbvq5ptv1rFjx2r8uABqwcdf3AmgHlm+fLkVERHhWu7Vq5f12GOPuc155ZVXrOjoaNeyJOtvf/uba3nr1q2WJGvp0qWusVWrVlmhoaGu5ZkzZ1oBAQHWoUOHXGPvv/++1aBBAysvL8+yLMuKiYmxXnvtNbfnfuSRR6yEhATLsixr//79liRrwYIF1e7XwIEDrSlTpriWb7zxRmvixIkX3HfLsqy3337b+uXL6MyZM62goCCroKDANfbRRx9Z4eHh1qlTp9y2vfzyy60lS5ZUWxuA2uOeGwAe2759uz7//HO3MzXl5eU6deqUSkpKFBYWJkn67W9/61ofFRUlSerSpYvb2KlTp1RUVKTw8HBJUuvWrdWqVSvXnISEBFVUVGj37t0KCAjQoUOHlJycrHvvvdc15/Tp04qIiHCrsVu3bm7L5eXlevzxx7V69WodPnxYpaWlKi0tVaNGjey2Q5IUGxurZs2auZa3b9+uEydOKDIy0m3e//3f/7ldigPgPYQbAB6rqKjQ7NmzNWTIkErrQkNDXf//y5tqz96jUtVYRUXFeZ/r7ByHw+Ga9+KLL+r66693mxcQEOC2fG5omT9/vv75z39qwYIF6tKlixo1aqRJkyaprKzs/DsqqUGDBrIsy23M6XRWmnfu81VUVCg6OlqffPJJpbnn3sMDwDsINwA81rVrV+3evVvt2rXz+mPn5ubqyJEjatmypSRp69atatCggTp06KCoqChddtll2rdvn/785z/X6nE3b96s22+/XSNGjJB0Jnzs2bNHnTp1cs0JDg5WeXm523bNmjVTcXGxTp486QowZ++puZCuXbsqPz9fgYGBatOmTa1qBeAZwg0Ajz388MO65ZZbFBMTozvvvFMNGjTQl19+qa+++kpz58619dihoaEaOXKk/vGPf6ioqEipqam666671KJFC0ln3pmUmpqq8PBwDRgwQKWlpdq2bZuOHz+utLS08z5uu3bt9NZbb2nLli269NJLlZ6ervz8fLdw06ZNG3322Wc6cOCAfvOb36hJkya6/vrrFRYWpr/+9a+aMGGC/v3vfysjI6Pa/fj973+vhIQEDR48WE888YSuuOIKHTlyRJmZmRo8eHCly2YA7OPdUgA81r9/f7333nvKysrS7373O91www1KT09XbGys7cdu166dhgwZooEDByoxMVGdO3d2e6v32LFj9dJLLykjI0NdunTRjTfeqIyMDMXFxV3wcf/+97+ra9eu6t+/v/r06aMWLVpo8ODBbnOmTp2qgIAAXXnllWrWrJlyc3PVpEkTrVy5UpmZma63j8+aNava/XA4HMrMzFTv3r01ZswYdejQQXfffbcOHDjguv8IgHc5rHMvIgOAj82aNUvvvPNOjS77AMC5OHMDAACMQrgBAABG4bIUAAAwCmduAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBR/h9DigUvvZMUYwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n",
"import matplotlib.pyplot as plt\n",
"\n",
"data[\"Frequency\"]=data.Malfunction/data.Count\n",
"data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Logistic regression\n",
"\n",
"Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Generalized Linear Model Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>Frequency</td> <th> No. Observations: </th> <td> 23</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>GLM</td> <th> Df Residuals: </th> <td> 21</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model Family:</th> <td>Binomial</td> <th> Df Model: </th> <td> 1</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Link Function:</th> <td>logit</td> <th> Scale: </th> <td> 1.0000</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>IRLS</td> <th> Log-Likelihood: </th> <td> -3.9210</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Thu, 22 Aug 2024</td> <th> Deviance: </th> <td> 3.0144</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:07:03</td> <th> Pearson chi2: </th> <td> 5.00</td> \n",
"</tr>\n",
"<tr>\n",
" <th>No. Iterations:</th> <td>6</td> <th> Pseudo R-squ. (CS):</th> <td>0.04355</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 5.0850</td> <td> 7.477</td> <td> 0.680</td> <td> 0.496</td> <td> -9.570</td> <td> 19.740</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Temperature</th> <td> -0.1156</td> <td> 0.115</td> <td> -1.004</td> <td> 0.316</td> <td> -0.341</td> <td> 0.110</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Generalized Linear Model Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Frequency No. Observations: 23\n",
"Model: GLM Df Residuals: 21\n",
"Model Family: Binomial Df Model: 1\n",
"Link Function: logit Scale: 1.0000\n",
"Method: IRLS Log-Likelihood: -3.9210\n",
"Date: Thu, 22 Aug 2024 Deviance: 3.0144\n",
"Time: 09:07:03 Pearson chi2: 5.00\n",
"No. Iterations: 6 Pseudo R-squ. (CS): 0.04355\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": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import statsmodels.api as sm\n",
"\n",
"data[\"Success\"]=data.Count-data.Malfunction\n",
"data[\"Intercept\"]=1\n",
"\n",
"logmodel=sm.GLM(data['Frequency'],\n",
" data[['Intercept','Temperature']],\n",
" family=sm.families.Binomial(sm.families.links.logit())).fit() # Added \"()\" after \"logit\"\n",
"\n",
"logmodel.summary()"
]
},
{
"cell_type": "markdown",
"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": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>Generalized Linear Model Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>Frequency</td> <th> No. Observations: </th> <td> 23</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>GLM</td> <th> Df Residuals: </th> <td> 21</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Model Family:</th> <td>Binomial</td> <th> Df Model: </th> <td> 1</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Link Function:</th> <td>logit</td> <th> Scale: </th> <td> 1.0000</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>IRLS</td> <th> Log-Likelihood: </th> <td> -23.526</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Thu, 22 Aug 2024</td> <th> Deviance: </th> <td> 18.086</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:09:20</td> <th> Pearson chi2: </th> <td> 30.0</td> \n",
"</tr>\n",
"<tr>\n",
" <th>No. Iterations:</th> <td>6</td> <th> Pseudo R-squ. (CS):</th> <td>0.2344</td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>nonrobust</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 5.0850</td> <td> 3.052</td> <td> 1.666</td> <td> 0.096</td> <td> -0.898</td> <td> 11.068</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Temperature</th> <td> -0.1156</td> <td> 0.047</td> <td> -2.458</td> <td> 0.014</td> <td> -0.208</td> <td> -0.023</td>\n",
"</tr>\n",
"</table>"
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" Generalized Linear Model Regression Results \n",
"==============================================================================\n",
"Dep. Variable: Frequency No. Observations: 23\n",
"Model: GLM Df Residuals: 21\n",
"Model Family: Binomial Df Model: 1\n",
"Link Function: logit Scale: 1.0000\n",
"Method: IRLS Log-Likelihood: -23.526\n",
"Date: Thu, 22 Aug 2024 Deviance: 18.086\n",
"Time: 09:09:20 Pearson chi2: 30.0\n",
"No. Iterations: 6 Pseudo R-squ. (CS): 0.2344\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": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n",
" family=sm.families.Binomial(sm.families.links.logit()), # Added \"()\" after \"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": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzZUlEQVR4nO3deXBUVf7+8aezE0KCBMiiAcIii+BCGDQoiwtBQNRyRlBkE3DMgCLGDWSURQUdlUFLiYhsEhdUlK84EchYsggoEogDhgKUYFA7k5+gJBqTdNL39weTlqazdGfhQHi/qrqKe+65957+5EIe7mqzLMsSAACAIX6mBwAAAM5thBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABglM9hZPPmzRo2bJhiY2Nls9m0Zs2aGpfZtGmTEhISFBISovbt2+uVV16pzVgBAEAj5HMY+e2333TJJZfopZde8qp/Tk6OhgwZor59+2r37t169NFHNWXKFK1evdrnwQIAgMbHVpcX5dlsNn3wwQe6+eabq+zzyCOP6MMPP9S+fftcbcnJyfrqq6+0ffv22m4aAAA0EgENvYHt27crKSnJrW3QoEFasmSJHA6HAgMDPZYpKSlRSUmJa9rpdOrYsWOKjIyUzWZr6CEDAIB6YFmWCgsLFRsbKz+/qk/GNHgYycvLU1RUlFtbVFSUysrK9NNPPykmJsZjmXnz5mn27NkNPTQAAHAaHDlyRBdccEGV8xs8jEjyOJpRcWaoqqMc06dPV0pKimv6+PHjatOmjXJyctSsWbN6GZNlWSosKtGmzZvVv18/BQSellKc1cocZdTLS9TKe9TKe9TKe9TKexW1GnTtAAUFBdXrugsLCxUfH1/j7+4G/wlFR0crLy/PrS0/P18BAQGKjIysdJng4GAFBwd7tLdo0ULh4eH1NrYIh0PnNQvVBTGtKz1dBHcO6uU1auU9auU9auU9auW9ilq1bNmy3mtVsb6aLrFo8OeMJCYmKiMjw61tw4YN6tWrFzsIAADwPYz8+uuvysrKUlZWlqQTt+5mZWUpNzdX0olTLGPGjHH1T05O1nfffaeUlBTt27dPS5cu1ZIlS/Tggw/WzzcAAABnNZ9P0+zcuVNXX321a7ri2o6xY8dq+fLlstvtrmAiSfHx8UpPT9f999+vl19+WbGxsXrxxRf15z//uR6GDwAAznY+h5EBAwaoukeTLF++3KOtf//+2rVrl6+bAgA0EuXl5XI4HKdtew6HQwEBASouLlZ5eflp2+7ZqC61CgwMlL+/f53HwCXGAIAGY1mW8vLy9Msvv5z27UZHR+vIkSM8n6oGda1V8+bNFR0dXac6E0YAAA2mIoi0bt1aoaGhpy0YOJ1O/frrrwoLC6v2YVuofa0sy1JRUZHy8/MlqdLnhnmLMAIAaBDl5eWuIFLVoxwaitPpVGlpqUJCQggjNahLrZo0aSLpxCM7WrduXetTNvyEAAANouIakdDQUMMjQUOq+PnW5ZogwggAoEFxzUbjVh8/X8IIAAAwijACAACMIowAAHCKcePGyWazeXy++eYb00NrlLibBgCASlx//fVatmyZW1urVq3cpktLS+v9TbfnIo6MAABQieDgYEVHR7t9rr32Wt1zzz1KSUlRy5YtNXDgQElSdna2hgwZorCwMEVFRWn06NH66aefXOv67bffNGbMGIWFhSkmJkbPP/+8BgwYoKlTp7r62Gw2rVmzxm0MzZs3d3uy+Q8//KARI0bovPPOU2RkpG666SYdPnzYNX/cuHG6+eab9dxzzykmJkaRkZGaPHmy250uJSUlevjhhxUXF6fg4GB17txZK1eulGVZ6tixo5577jm3Mezdu1d+fn769ttv617UKhBGAACnjWVZKiotOy2f30vLXX+u7jUmvlqxYoUCAgK0detWLVq0SHa7Xf3799ell16qnTt3at26dfrvf/+r4cOHu5Z56KGH9Omnn+qDDz7Qhg0btHHjRmVmZvq03aKiIl199dUKCwvT5s2b9dlnnyksLEzXX3+9SktLXf0+/fRTffvtt/r000+1YsUKLV++3C3QjBkzRm+//bZefPFF7du3TwsXLlTTpk1ls9k0fvx4j6NBS5cuVd++fdWhQ4faFcwLnKYBAJw2vzvK1e3x9ad9u9lzBik0yLdfeR999JHCwsJc04MHD5YkdezYUf/4xz9c7Y8//rh69uypuXPnutqWLl2quLg4HThwQLGxsVqyZIlef/1115GUFStW6IILLvBpPG+//bb8/Pz02muvuW6nXbZsmZo3b66NGzcqKSlJknTeeefppZdekr+/v7p06aKhQ4fqk08+0V133aUDBw7onXfeUUZGhq677jpJUrt27VRQUCBJuvPOO/X4449rx44d6t27txwOh9LS0vTss8/6NFZfEUYAAKjE1VdfrdTUVNd006ZNdfvtt6tXr15u/TIzM/Xpp5+6BZcK3377rX7//XeVlpYqMTHR1d6iRQt17tzZp/FkZmbqm2++UbNmzdzai4uL3U6hXHTRRW5PQo2JidGePXskSVlZWfL391f//v0r3UZMTIyGDh2qpUuXqnfv3vroo49UXFysW2+91aex+oowAgA4bZoE+it7zqAG347T6VRhQaGahTeTn5+fmgT6/pjypk2bqmPHjpW2n7qtYcOG6ZlnnvHoGxMTo4MHD3q1PZvN5nE66eRrPZxOpxISEvTGG294LHvyhbWBgYEe63U6nZL+eHx7dSZOnKjRo0frn//8p5YtW6YRI0Y0+FN0CSMAgNPGZrP5fLqkNpxOp8qC/BUaFNDg76bp2bOnVq9erXbt2ikgwPO7dezYUYGBgfr888/Vpk0bSdLPP/+sAwcOuB2haNWqlex2u2v64MGDKioqctvOqlWr1Lp1a4WHh9dqrD169JDT6dSmTZtcp2lONWTIEDVt2lSpqan6+OOPtXnz5lptyxdcwAoAQB1MnjxZx44d0+23364dO3bo0KFD2rBhg8aPH6/y8nKFhYVpwoQJeuihh/TJJ59o7969GjdunEdIuuaaa/TSSy9p165d2rlzp5KTk92Octxxxx1q2bKlbrrpJm3ZskU5OTnatGmT7rvvPn3//fdejbVdu3YaO3asxo8frzVr1ignJ0cbN27UBx984Orj7++vcePGafr06erYsaPb6aWGQhgBAKAOYmNjtXXrVpWXl2vQoEHq3r277rvvPkVERLgCx7PPPqt+/frpxhtv1HXXXaerrrpKCQkJbut5/vnnFRcXp379+mnkyJF68MEH3U6PhIaGavPmzWrTpo1uueUWde3aVePHj9fvv//u05GS1NRU/eUvf9GkSZPUpUsX3X333W5HYCRpwoQJKi0t1fjx4+tQGe9xmgYAgFOcfCvsyTZu3Fhpe6dOnfT+++9Xub6wsDCtXLlSK1eudLX961//cusTGxur9evd7zT65Zdf3Kajo6O1YsUKn8a9YMECt+mQkBDNnz9f8+fPl3TilFbF3TQV7Ha7AgICNGbMmCq3VZ8IIwAAQNKJB6IdOXJEjz32mIYPH66oqKjTsl1O0wAAAEnSW2+9pc6dO+v48eNuz1JpaBwZAQDAgKpO+Zg0btw4jRs37rRvlyMjAADAKMIIAKBB1ed7YXDmqY+fL2EEANAgKp6Rcepto2hcKn6+pz751RdcMwIAaBD+/v5q3ry58vPzJZ14TkbFC94amtPpVGlpqYqLixv8Caxnu9rWyrIsFRUVKT8/X82bN3d7H46vCCMAgAYTHR0tSa5AcrpYlqXff/9dTZo0OW0B6GxV11o1b97c9XOuLcIIAKDB2Gw2xcTEqHXr1m4vfWtoDodDmzdvVr9+/ep0+uBcUJdaBQYG1umISAXCCACgwfn7+9fLLy1ftldWVqaQkBDCSA3OhFpxIg0AABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRtQojCxcuVHx8vEJCQpSQkKAtW7ZU2/+NN97QJZdcotDQUMXExOjOO+/U0aNHazVgAADQuPgcRlatWqWpU6dqxowZ2r17t/r27avBgwcrNze30v6fffaZxowZowkTJujrr7/Wu+++qy+//FITJ06s8+ABAMDZz+cwMn/+fE2YMEETJ05U165dtWDBAsXFxSk1NbXS/p9//rnatWunKVOmKD4+XldddZXuvvtu7dy5s86DBwAAZ78AXzqXlpYqMzNT06ZNc2tPSkrStm3bKl2mT58+mjFjhtLT0zV48GDl5+frvffe09ChQ6vcTklJiUpKSlzTBQUFkiSHwyGHw+HLkKtVsa76XGdjRr28R628R628R628R62815C18nadNsuyLG9X+uOPP+r888/X1q1b1adPH1f73LlztWLFCu3fv7/S5d577z3deeedKi4uVllZmW688Ua99957CgwMrLT/rFmzNHv2bI/2N998U6Ghod4OFwAAGFRUVKSRI0fq+PHjCg8Pr7KfT0dGKthsNrdpy7I82ipkZ2drypQpevzxxzVo0CDZ7XY99NBDSk5O1pIlSypdZvr06UpJSXFNFxQUKC4uTklJSdV+GV85HA5lZGRo4MCBVQYj/IF6eY9aeY9aeY9aeY9aea8ha1VxZqMmPoWRli1byt/fX3l5eW7t+fn5ioqKqnSZefPm6corr9RDDz0kSbr44ovVtGlT9e3bV08++aRiYmI8lgkODlZwcLBHe2BgYIPsVA213saKenmPWnmPWnmPWnmPWnmvIWrl7fp8uoA1KChICQkJysjIcGvPyMhwO21zsqKiIvn5uW/G399f0okjKgAA4Nzm8900KSkpeu2117R06VLt27dP999/v3Jzc5WcnCzpxCmWMWPGuPoPGzZM77//vlJTU3Xo0CFt3bpVU6ZMUe/evRUbG1t/3wQAAJyVfL5mZMSIETp69KjmzJkju92u7t27Kz09XW3btpUk2e12t2eOjBs3ToWFhXrppZf0wAMPqHnz5rrmmmv0zDPP1N+3AAAAZ61aXcA6adIkTZo0qdJ5y5cv92i79957de+999ZmUwAAoJHj3TQAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIyqVRhZuHCh4uPjFRISooSEBG3ZsqXa/iUlJZoxY4batm2r4OBgdejQQUuXLq3VgAEAQOMS4OsCq1at0tSpU7Vw4UJdeeWVWrRokQYPHqzs7Gy1adOm0mWGDx+u//73v1qyZIk6duyo/Px8lZWV1XnwAADg7OdzGJk/f74mTJigiRMnSpIWLFig9evXKzU1VfPmzfPov27dOm3atEmHDh1SixYtJEnt2rWr26gBAECj4VMYKS0tVWZmpqZNm+bWnpSUpG3btlW6zIcffqhevXrpH//4h1auXKmmTZvqxhtv1BNPPKEmTZpUukxJSYlKSkpc0wUFBZIkh8Mhh8Phy5CrVbGu+lxnY0a9vEetvEetvEetvEetvNeQtfJ2nT6FkZ9++knl5eWKiopya4+KilJeXl6lyxw6dEifffaZQkJC9MEHH+inn37SpEmTdOzYsSqvG5k3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/Xw+TSNJNpvNbdqyLI+2Ck6nUzabTW+88YYiIiIknTjV85e//EUvv/xypUdHpk+frpSUFNd0QUGB4uLilJSUpPDw8NoMuVIOh0MZGRkaOHCgAgMD6229jRX18h618h618h618h618l5D1qrizEZNfAojLVu2lL+/v8dRkPz8fI+jJRViYmJ0/vnnu4KIJHXt2lWWZen7779Xp06dPJYJDg5WcHCwR3tgYGCD7FQNtd7Ginp5j1p5j1p5j1p5j1p5ryFq5e36fLq1NygoSAkJCR6HcjIyMtSnT59Kl7nyyiv1448/6tdff3W1HThwQH5+frrgggt82TwAAGiEfH7OSEpKil577TUtXbpU+/bt0/3336/c3FwlJydLOnGKZcyYMa7+I0eOVGRkpO68805lZ2dr8+bNeuihhzR+/PgqL2AFAADnDp+vGRkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7Xbm5ua7+YWFhysjI0L333qtevXopMjJSw4cP15NPPll/3wIAAJy1anUB66RJkzRp0qRK5y1fvtyjrUuXLlzRDAAAKsW7aQAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGFWrMLJw4ULFx8crJCRECQkJ2rJli1fLbd26VQEBAbr00ktrs1kAANAI+RxGVq1apalTp2rGjBnavXu3+vbtq8GDBys3N7fa5Y4fP64xY8bo2muvrfVgAQBA4+NzGJk/f74mTJigiRMnqmvXrlqwYIHi4uKUmppa7XJ33323Ro4cqcTExFoPFgAAND4BvnQuLS1VZmampk2b5taelJSkbdu2VbncsmXL9O233yotLU1PPvlkjdspKSlRSUmJa7qgoECS5HA45HA4fBlytSrWVZ/rbMyol/eolfeolfeolfeolfcaslbertOnMPLTTz+pvLxcUVFRbu1RUVHKy8urdJmDBw9q2rRp2rJliwICvNvcvHnzNHv2bI/2DRs2KDQ01JcheyUjI6Pe19mYUS/vUSvvUSvvUSvvUSvvNUStioqKvOrnUxipYLPZ3KYty/Jok6Ty8nKNHDlSs2fP1oUXXuj1+qdPn66UlBTXdEFBgeLi4pSUlKTw8PDaDLlSDodDGRkZGjhwoAIDA+ttvY0V9fIetfIetfIetfIetfJeQ9aq4sxGTXwKIy1btpS/v7/HUZD8/HyPoyWSVFhYqJ07d2r37t265557JElOp1OWZSkgIEAbNmzQNddc47FccHCwgoODPdoDAwMbZKdqqPU2VtTLe9TKe9TKe9TKe9TKew1RK2/X59MFrEFBQUpISPA4lJORkaE+ffp49A8PD9eePXuUlZXl+iQnJ6tz587KysrS5Zdf7svmAQBAI+TzaZqUlBSNHj1avXr1UmJiol599VXl5uYqOTlZ0olTLD/88INef/11+fn5qXv37m7Lt27dWiEhIR7tAADg3ORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7vcZnjgAAAFSo1QWskyZN0qRJkyqdt3z58mqXnTVrlmbNmlWbzQIAgEaId9MAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjavXWXgCnX7nT0o6cY8ovLFbrZiHqHd9C/n4208PCOY79EvWBMAKcBdbttWv22mzZjxe72mIiQjRzWDdd3z3G4MhwLmO/RH3hNA1whlu3166/pe1y+wdfkvKOF+tvabu0bq/d0MhwLmO/RH0ijABnsHKnpdlrs2VVMq+ibfbabJU7K+sBNAz2S9Q3wghwBtuRc8zjf54nsyTZjxdrR86x0zconPPYL1HfCCPAGSy/sOp/8GvTD6gP7Jeob4QR4AzWullIvfYD6gP7JeobYQQ4g/WOb6GYiBBVdaOkTSfuXugd3+J0DgvnOPZL1DfCCHAG8/ezaeawbpLk8Q9/xfTMYd14rgNOK/ZL1DfCCHCGu757jFJH9VR0hPsh7+iIEKWO6snzHGAE+yXqEw89A84C13eP0cBu0TzpEmcU9kvUF8IIcJbw97MpsUOk6WEAbtgvUR84TQMAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKNqFUYWLlyo+Ph4hYSEKCEhQVu2bKmy7/vvv6+BAweqVatWCg8PV2JiotavX1/rAQMAgMbF5zCyatUqTZ06VTNmzNDu3bvVt29fDR48WLm5uZX237x5swYOHKj09HRlZmbq6quv1rBhw7R79+46Dx4AAJz9fA4j8+fP14QJEzRx4kR17dpVCxYsUFxcnFJTUyvtv2DBAj388MP605/+pE6dOmnu3Lnq1KmT1q5dW+fBAwCAs1+AL51LS0uVmZmpadOmubUnJSVp27ZtXq3D6XSqsLBQLVq0qLJPSUmJSkpKXNMFBQWSJIfDIYfD4cuQq1WxrvpcZ2NGvbxHrbxHrbxHrbxHrbzXkLXydp0+hZGffvpJ5eXlioqKcmuPiopSXl6eV+t4/vnn9dtvv2n48OFV9pk3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/XwKIxVsNpvbtGVZHm2VeeuttzRr1iz93//9n1q3bl1lv+nTpyslJcU1XVBQoLi4OCUlJSk8PLw2Q66Uw+FQRkaGBg4cqMDAwHpbb2NFvbxHrbxHrbxHrbxHrbzXkLWqOLNRE5/CSMuWLeXv7+9xFCQ/P9/jaMmpVq1apQkTJujdd9/VddddV23f4OBgBQcHe7QHBgY2yE7VUOttrKiX96iV96iV96iV96iV9xqiVt6uz6cLWIOCgpSQkOBxKCcjI0N9+vSpcrm33npL48aN05tvvqmhQ4f6skkAANDI+XyaJiUlRaNHj1avXr2UmJioV199Vbm5uUpOTpZ04hTLDz/8oNdff13SiSAyZswYvfDCC7riiitcR1WaNGmiiIiIevwqAADgbORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa73e2ZI4sWLVJZWZkmT56syZMnu9rHjh2r5cuX1/0bAACAs1qtLmCdNGmSJk2aVOm8UwPGxo0ba7MJAABwjqhVGAFw7ih3WtqRc0z5hcVq3SxEveNbyN/P5vV8E87EMdVVaZlTadsPK1LSyu2HNapPBwUF8HoxNA6EEQBVWrfXrtlrs2U/Xuxqi4kI0cxh3XR995ga55twJo6prualZ2vxlhwF+ln6R2/pmfX79eTHB3RX33hNH9LN9PCAOiNWA6jUur12/S1tl9svdUnKO16sv6Xt0rz07Grnr9trP53DlVTzmE2Mqa7mpWdr0eYcOS33dqclLdqco3np2WYGBtQjwggAD+VOS7PXZsuqZJ71v8/iLTlVzpek2WuzVX7qb9AGVNOYTYyprkrLnFq8JafaPou35Ki0zHmaRgQ0DMIIAA87co55HF04VXW/0y1J9uPF2pFzrH4HVo2axmxiTHW1cvvhaussnfg5rNx++LSMB2gohBEAHvILqw8ip3s99bmt0zmmuvrumHfv9fC2H3CmIowA8NC6WcgZtZ763NbpHFNdtW3h3YtBve0HnKkIIwA89I5voZiIEFV3M6yfTVXOt+nEHSy941s0wOgqV9OYTYyprkYntlNNdyT72U70A85mhBEAHvz9bJo57MQto6f+LrT973NX3/gq50vSzGHdTuuzPWoas4kx1VVQgJ+rzlW5q288zxvBWY89GEClru8eo9RRPRUd4X5aIzoiRKmjemr6kG7VzjfxTI+axnw2Pmdk+pBuurtfvMcREj+bdHc/njOCxoGHngGo0vXdYzSwW3SVTzOtaf6ZOOaz0fQh3fRAUhelbftW+jlbjwzqzBNY0agQRgBUy9/PpsQOkbWeb8KZOKa6Cgrw0+jEdkpPz9boxHYKJIigEWFvBgAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUQGmBwAAZ5Nyp6UdOceUX1is1s1C1Du+hfz9bJKk30vLNTc9W4ePFqldZKgeHdJNTYL8vVq2unmSVFrmVNr2w4qUtHL7YY3q00FBAd79f7Kmddc0v7brLi1zauX2w/ruWJHatgjV6MR2Z/yYYUatwsjChQv17LPPym6366KLLtKCBQvUt2/fKvtv2rRJKSkp+vrrrxUbG6uHH35YycnJtR40AJiwbq9ds9dmy3682NUWExGimcO6afWu75WRne9q33JQWvl5rgZ2a63FY/5U7bKSqpx3ffcYzUvP1uItOQr0s/SP3tIz6/fryY8P6K6+8Zo+pFutx3x995ga59d23btzf9biLTlyWn/0fyp93xk95pqWRcPxOYysWrVKU6dO1cKFC3XllVdq0aJFGjx4sLKzs9WmTRuP/jk5ORoyZIjuuusupaWlaevWrZo0aZJatWqlP//5z/XyJQCgoa3ba9ff0nbJOqU973ixktN2VblcRna+bnxpi/Z8X+DTsnnHi/W3tF26rltrt5BTwWlJizbnSFKVv9yrG/Pf0nbpr/3i9ermnCrnp47qWeUv6NrU40wec03LomH5fM3I/PnzNWHCBE2cOFFdu3bVggULFBcXp9TU1Er7v/LKK2rTpo0WLFigrl27auLEiRo/fryee+65Og8eAE6Hcqel2WuzPX6JSaq07VT/qSSI1LSs9b9PZUHkZIu35Ki0zOnRXtOYrf8tW924Zq/NVrnTs0dd63GmjrmqZdHwfDoyUlpaqszMTE2bNs2tPSkpSdu2bat0me3btyspKcmtbdCgQVqyZIkcDocCAwM9likpKVFJSYlr+vjx45KkY8eOyeFw+DLkajkcDhUVFeno0aOVjgPuqJf3qJX3zoZaZX73s/7f0aPGL7ILcFoqKnIqwOGncucf1zgszvhKt/V2PzLt7Zir+x/p/zv6mz7J+lYJbc+r1bqr09Bjvjg2zG2/8mbdVX3fxq4h/w4WFhZKkiyrhpBn+eCHH36wJFlbt251a3/qqaesCy+8sNJlOnXqZD311FNubVu3brUkWT/++GOly8ycObMiBPPhw4cPHz58zvLPkSNHqs0XtQq2Npv7VceWZXm01dS/svYK06dPV0pKimva6XTq2LFjioyMrHY7viooKFBcXJyOHDmi8PDweltvY0W9vEetvEetvEetvEetvNeQtbIsS4WFhYqNja22n09hpGXLlvL391deXp5be35+vqKioipdJjo6utL+AQEBioyMrHSZ4OBgBQcHu7U1b97cl6H6JDw8nJ3VB9TLe9TKe9TKe9TKe9TKew1Vq4iIiBr7+HQBa1BQkBISEpSRkeHWnpGRoT59+lS6TGJiokf/DRs2qFevXmfs+WEAAHD6+Hw3TUpKil577TUtXbpU+/bt0/3336/c3FzXc0OmT5+uMWPGuPonJyfru+++U0pKivbt26elS5dqyZIlevDBB+vvWwAAgLOWz9eMjBgxQkePHtWcOXNkt9vVvXt3paenq23btpIku92u3NxcV//4+Hilp6fr/vvv18svv6zY2Fi9+OKLZ8QzRoKDgzVz5kyPU0KoHPXyHrXyHrXyHrXyHrXy3plQK5tl1XS/DQAAQMPhRXkAAMAowggAADCKMAIAAIwijAAAAKPOiTCSmpqqiy++2PVAl8TERH388ceu+ZZladasWYqNjVWTJk00YMAAff311wZHfGaYN2+ebDabpk6d6mqjVn+YNWuWbDab2yc6Oto1n1q5++GHHzRq1ChFRkYqNDRUl156qTIzM13zqdcJ7dq189ivbDabJk+eLIk6naysrEx///vfFR8fryZNmqh9+/aaM2eOnM4/XsJHvf5QWFioqVOnqm3btmrSpIn69OmjL7/80jXfaK1qeB1No/Dhhx9a//rXv6z9+/db+/fvtx599FErMDDQ2rt3r2VZlvX0009bzZo1s1avXm3t2bPHGjFihBUTE2MVFBQYHrk5O3bssNq1a2ddfPHF1n333edqp1Z/mDlzpnXRRRdZdrvd9cnPz3fNp1Z/OHbsmNW2bVtr3Lhx1hdffGHl5ORY//73v61vvvnG1Yd6nZCfn++2T2VkZFiSrE8//dSyLOp0sieffNKKjIy0PvroIysnJ8d69913rbCwMGvBggWuPtTrD8OHD7e6detmbdq0yTp48KA1c+ZMKzw83Pr+++8tyzJbq3MijFTmvPPOs1577TXL6XRa0dHR1tNPP+2aV1xcbEVERFivvPKKwRGaU1hYaHXq1MnKyMiw+vfv7woj1MrdzJkzrUsuuaTSedTK3SOPPGJdddVVVc6nXlW77777rA4dOlhOp5M6nWLo0KHW+PHj3dpuueUWa9SoUZZlsV+drKioyPL397c++ugjt/ZLLrnEmjFjhvFanROnaU5WXl6ut99+W7/99psSExOVk5OjvLw8JSUlufoEBwerf//+2rZtm8GRmjN58mQNHTpU1113nVs7tfJ08OBBxcbGKj4+XrfddpsOHTokiVqd6sMPP1SvXr106623qnXr1rrsssu0ePFi13zqVbnS0lKlpaVp/Pjxstls1OkUV111lT755BMdOHBAkvTVV1/ps88+05AhQySxX52srKxM5eXlCgkJcWtv0qSJPvvsM+O1OmfCyJ49exQWFqbg4GAlJyfrgw8+ULdu3Vwv8Tv1RX9RUVEeL/g7F7z99tvatWuX5s2b5zGPWrm7/PLL9frrr2v9+vVavHix8vLy1KdPHx09epRaneLQoUNKTU1Vp06dtH79eiUnJ2vKlCl6/fXXJbFvVWXNmjX65ZdfNG7cOEnU6VSPPPKIbr/9dnXp0kWBgYG67LLLNHXqVN1+++2SqNfJmjVrpsTERD3xxBP68ccfVV5errS0NH3xxRey2+3Ga+Xz4+DPVp07d1ZWVpZ++eUXrV69WmPHjtWmTZtc8202m1t/y7I82hq7I0eO6L777tOGDRs80vPJqNUJgwcPdv25R48eSkxMVIcOHbRixQpdccUVkqhVBafTqV69emnu3LmSpMsuu0xff/21UlNT3d5lRb3cLVmyRIMHD/Z4/Tp1OmHVqlVKS0vTm2++qYsuukhZWVmaOnWqYmNjNXbsWFc/6nXCypUrNX78eJ1//vny9/dXz549NXLkSO3atcvVx1StzpkjI0FBQerYsaN69eqlefPm6ZJLLtELL7zguvvh1OSXn5/vkRAbu8zMTOXn5yshIUEBAQEKCAjQpk2b9OKLLyogIMBVD2pVuaZNm6pHjx46ePAg+9UpYmJi1K1bN7e2rl27ut5jRb08fffdd/r3v/+tiRMnutqok7uHHnpI06ZN02233aYePXpo9OjRuv/++11HdqmXuw4dOmjTpk369ddfdeTIEe3YsUMOh0Px8fHGa3XOhJFTWZalkpIS1w8hIyPDNa+0tFSbNm1Snz59DI7w9Lv22mu1Z88eZWVluT69evXSHXfcoaysLLVv355aVaOkpET79u1TTEwM+9UprrzySu3fv9+t7cCBA64XbFIvT8uWLVPr1q01dOhQVxt1cldUVCQ/P/dfY/7+/q5be6lX5Zo2baqYmBj9/PPPWr9+vW666SbztWrwS2TPANOnT7c2b95s5eTkWP/5z3+sRx991PLz87M2bNhgWdaJ25kiIiKs999/39qzZ491++23n7O3fp3q5LtpLItaneyBBx6wNm7caB06dMj6/PPPrRtuuMFq1qyZdfjwYcuyqNXJduzYYQUEBFhPPfWUdfDgQeuNN96wQkNDrbS0NFcf6vWH8vJyq02bNtYjjzziMY86/WHs2LHW+eef77q19/3337datmxpPfzww64+1OsP69atsz7++GPr0KFD1oYNG6xLLrnE6t27t1VaWmpZltlanRNhZPz48Vbbtm2toKAgq1WrVta1117rCiKWdeL2r5kzZ1rR0dFWcHCw1a9fP2vPnj0GR3zmODWMUKs/VNyDHxgYaMXGxlq33HKL9fXXX7vmUyt3a9eutbp3724FBwdbXbp0sV599VW3+dTrD+vXr7ckWfv37/eYR53+UFBQYN13331WmzZtrJCQEKt9+/bWjBkzrJKSElcf6vWHVatWWe3bt7eCgoKs6Ohoa/LkydYvv/zimm+yVjbLsqyGP/4CAABQuXP2mhEAAHBmIIwAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAjRCNput2k/FK+kbkwEDBmjq1KmmhwGgFgJMDwBA/bPb7a4/r1q1So8//rjbi+qaNGliYli14nA4FBgY2Gi3B4AjI0CjFB0d7fpERETIZrO5tW3evFkJCQkKCQlR+/btNXv2bJWVlbmWt9lsWrRokW644QaFhoaqa9eu2r59u7755hsNGDBATZs2VWJior799lvXMrNmzdKll16qRYsWKS4uTqGhobr11lv1yy+/uI1t2bJl6tq1q0JCQtSlSxctXLjQNe/w4cOy2Wx65513NGDAAIWEhCgtLU1Hjx7V7bffrgsuuEChoaHq0aOH3nrrLddy48aN06ZNm/TCCy+4jv4cPnxYy5cvV/Pmzd22v2bNGtlsNo9xL126VO3bt1dwcLAsy9Lx48f117/+Va1bt1Z4eLiuueYaffXVV/X0EwJwMsIIcI5Zv369Ro0apSlTpig7O1uLFi3S8uXL9dRTT7n1e+KJJzRmzBhlZWWpS5cuGjlypO6++25Nnz5dO3fulCTdc889bst88803euedd7R27VqtW7dOWVlZmjx5smv+4sWLNWPGDD311FPat2+f5s6dq8cee0wrVqxwW88jjzyiKVOmaN++fRo0aJCKi4uVkJCgjz76SHv37tVf//pXjR49Wl988YUk6YUXXlBiYqLuuusu2e122e12xcXFeV2TinGvXr1aWVlZkqShQ4cqLy9P6enpyszMVM+ePXXttdfq2LFjXq8XgJdOy+v4ABizbNkyKyIiwjXdt29fa+7cuW59Vq5cacXExLimJVl///vfXdPbt2+3JFlLlixxtb311ltWSEiIa3rmzJmWv7+/deTIEVfbxx9/bPn5+Vl2u92yLMuKi4uz3nzzTbdtP/HEE1ZiYqJlWZaVk5NjSbIWLFhQ4/caMmSI9cADD7imT33DdGXf3bIs64MPPrBO/qdv5syZVmBgoJWfn+9q++STT6zw8HCruLjYbdkOHTpYixYtqnFsAHzDNSPAOSYzM1Nffvml25GQ8vJyFRcXq6ioSKGhoZKkiy++2DU/KipKktSjRw+3tuLiYhUUFCg8PFyS1KZNG11wwQWuPomJiXI6ndq/f7/8/f115MgRTZgwQXfddZerT1lZmSIiItzG2KtXL7fp8vJyPf3001q1apV++OEHlZSUqKSkRE2bNq1rOSRJbdu2VatWrVzTmZmZ+vXXXxUZGenW7/fff3c7NQWgfhBGgHOM0+nU7Nmzdcstt3jMCwkJcf355Is4K66xqKzN6XRWua2KPjabzdVv8eLFuvzyy936+fv7u02fGjKef/55/fOf/9SCBQvUo0cPNW3aVFOnTlVpaWnVX1SSn5+fLMtya3M4HB79Tt2e0+lUTEyMNm7c6NH31GtQANQdYQQ4x/Ts2VP79+9Xx44d633dubm5+vHHHxUbGytJ2r59u/z8/HThhRcqKipK559/vg4dOqQ77rjDp/Vu2bJFN910k0aNGiXpRFg4ePCgunbt6uoTFBSk8vJyt+VatWqlwsJC/fbbb67AUXFNSHV69uypvLw8BQQEqF27dj6NFYDvCCPAOebxxx/XDTfcoLi4ON16663y8/PTf/7zH+3Zs0dPPvlkndYdEhKisWPH6rnnnlNBQYGmTJmi4cOHKzo6WtKJO1emTJmi8PBwDR48WCUlJdq5c6d+/vlnpaSkVLnejh07avXq1dq2bZvOO+88zZ8/X3l5eW5hpF27dvriiy90+PBhhYWFqUWLFrr88ssVGhqqRx99VPfee6927Nih5cuX1/g9rrvuOiUmJurmm2/WM888o86dO+vHH39Uenq6br75Zo/TSADqhrtpgHPMoEGD9NFHHykjI0N/+tOfdMUVV2j+/Plq27ZtndfdsWNH3XLLLRoyZIiSkpLUvXt3t1t3J06cqNdee03Lly9Xjx491L9/fy1fvlzx8fHVrvexxx5Tz549NWjQIA0YMEDR0dG6+eab3fo8+OCD8vf3V7du3dSqVSvl5uaqRYsWSktLU3p6uut24FmzZtX4PWw2m9LT09WvXz+NHz9eF154oW677TYdPnzYdf0MgPpjs049oQoAtTBr1iytWbPGq9MgAHAyjowAAACjCCMAAMAoTtMAAACjODICAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjPr/yAWkb36Wd7EAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n",
"data_pred['Frequency'] = logmodel.predict(data_pred)\n",
"data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n",
"plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"scrolled": true
},
"source": [
"This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing and plotting uncertainty"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG6CAYAAAALTELXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABThUlEQVR4nO3deXxU1f0//te9s2ayJ5CFHdlCIASQVUAQRf2BReWDrViQakVbF4q0gv26QNWKrSAtAuICYl0qKopVEUWtdWMXQYSwCYFA9j2Z/d7z+2OSIUMCJDM3mYT7ej4ekeTeO3PPvDMmr5xz7rmSEEKAiIiISKfkcDeAiIiIKJwYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI11pVGFq5ciVmzJhx3mNKS0vxxz/+EUOHDsXQoUPxyCOPwG63t1ALiYiI6GLTasLQ2rVrsWzZsgseN3v2bJw8edJ//Lfffou//OUvLdBCIiIiuhgZw92A/Px8PPTQQ9i1axe6d+9+3mN3796N7du3Y+PGjejRowcA4LHHHsMdd9yBuXPnIjk5uSWaTERERBeRsPcM/fTTT4iNjcV//vMfZGZmnvfYnTt3on379v4gBADDhg2DJEnYtWtXczeViIiILkJh7xkaP348xo8f36hj8/PzkZqaGrDNbDYjLi4Oubm5zdE8IiIiusiFvWeoKRwOB8xmc73tFosFLpcr6OcVQoTSLCIiImrDwt4z1BRWqxVut7vedpfLBZvNFvTzSpKEigoHFEUNpXm6ZzDIiImJYC1DxDpqh7XUDmupDdZRO7GxEZBlbfp02lQYSklJwWeffRawze12o6ysLOTJ04qiwuvlG1MLrKU2WEftsJbaYS21wTqGTstBnTY1TDZ06FDk5eUhOzvbv23btm0AgMGDB4erWURERNSGteowpCgKCgsL4XQ6AQCZmZkYPHgw7r//fuzduxdbt27FggULcMMNN/CyeiIiIgpKqw5Dubm5GD16NDZu3AjAN7dn+fLl6NSpE2bOnIk5c+bg8ssvx8KFC8PbUCIiImqzJMFLqQAApaXVHL8NkdEoIz4+krUMEeuoHdZSO6ylNlhH7SQkRMJg0KZPp1X3DBERERE1N4YhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0rWwhyFVVbFs2TKMGTMGmZmZuP3225GdnX3O4wsLCzF37lwMHz4cw4cPxx/+8Afk5eW1YIuJiIjoYhL2MLRy5Uq8+eabeOKJJ7Bu3TpIkoRZs2bB7XY3ePz999+P3NxcvPzyy3j55ZeRl5eHu+++u4VbTURERBeLsIYht9uNNWvW4L777sPYsWORlpaGpUuXIj8/H5s3b653fEVFBXbs2IFZs2YhPT0d6enpuPPOO/HTTz+htLQ0DK+AiIiI2jpjOE+elZWF6upqjBgxwr8tJiYG6enp2LFjByZNmhRwvMVigc1mw4YNGzBs2DAAwPvvv49u3bohNjY2pLYYDGHvJGvzamvIWoaGddQOa6kd1lIbrKN2JEm75wprGKqd65OamhqwPSkpCbm5ufWOt1gs+Otf/4rHHnsMQ4YMgSRJaN++PV577TXIcmhvrJiYiJAeT2ewltpgHbXDWmqHtdQG69i6hDUMORwOAIDZbA7YbrFYUF5eXu94IQQOHjyIQYMG4Y477oCiKFi6dCnuuece/Pvf/0ZUVFTQbamocEBR1KAfT76/dGJiIljLELGO2mEttcNaaoN11E5sbETIHSG1whqGrFYrAN/codrPAcDlciEion5q/uijj/DGG2/gv//9rz/4rFq1CldccQXWr1+PmTNnBt0WRVHh9fKNqQXWUhuso3ZYS+2wltpgHUMnhHbPFdZBy9rhsYKCgoDtBQUFSElJqXf8rl270L1794AeoNjYWHTv3h3Hjx9v1rYSERHRxSmsYSgtLQ1RUVHYtm2bf1tFRQX279+PIUOG1Ds+NTUV2dnZcLlc/m0OhwM5OTno2rVri7SZiIiILi5hDUNmsxnTp0/H4sWL8fnnnyMrKwv3338/UlJSMGHCBCiKgsLCQjidTgDADTfcAACYM2cOsrKy/MebzWZMmTIljK+EiIiI2qqwX9s3e/ZsTJ06FQ8//DCmTZsGg8GA1atXw2w2Izc3F6NHj8bGjRsB+K4ye+ONNyCEwMyZM3HbbbfBZDLh3//+N2JiYsL8SoiIiKgtkoTQcgpS21VaWs3JbCEyGmXEx0eyliFiHbXDWmqHtdQG66idhIRIzdZrCnvPEBEREVE4MQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRroU9DKmqimXLlmHMmDHIzMzE7bffjuzs7HMe7/F4sGTJEowZMwYDBw7E9OnTceDAgRZsMREREV1Mwh6GVq5ciTfffBNPPPEE1q1bB0mSMGvWLLjd7gaPX7hwId555x08/vjjWL9+PeLi4jBr1ixUVla2cMuJiIjoYhDWMOR2u7FmzRrcd999GDt2LNLS0rB06VLk5+dj8+bN9Y4/efIk3nnnHSxatAjjxo1Djx498OSTT8JsNmPfvn1heAVERETU1hnDefKsrCxUV1djxIgR/m0xMTFIT0/Hjh07MGnSpIDjv/nmG8TExODyyy8POP6LL74IuS0GQ9g7ydq82hqylqFhHbXDWmqHtdQG66gdSdLuucIahvLy8gAAqampAduTkpKQm5tb7/jjx4+jc+fO+PTTT/HCCy8gPz8f6enpePDBB9GjR4+Q2hITExHS4+kM1lIbrKN2WEvtsJbaYB1bl7CGIYfDAQAwm80B2y0WC8rLy+sdX1VVhRMnTmDlypWYN28eYmJi8Nxzz+GWW27Bxo0bkZiYGHRbKiocUBQ16MeT7y+dmJgI1jJErKN2WEvtsJbaYB21ExsbAVnWpoctqDC0YsUKTJkypV6PTlNZrVYAvrlDtZ8DgMvlQkRE/dRsMplQWVmJpUuX+nuCli5dirFjx+K9997DHXfcEXRbFEWF18s3phZYS22wjtphLbXDWmqDdQydENo9V1CR6pVXXsGVV16J2267DR988AFcLldQJ68NUwUFBQHbCwoKkJKSUu/4lJQUGI3GgCExq9WKzp07IycnJ6g2EBERkb4FFYa++eYbLF68GCaTCQ8++CBGjRqFRx99FLt3727S86SlpSEqKgrbtm3zb6uoqMD+/fsxZMiQescPGTIEXq8XP/74o3+b0+nEyZMn0bVr12BeChEREelcUMNkZrMZEydOxMSJE1FQUID//Oc/+Pjjj/H222+jW7dumDJlCqZMmXLBOTxmsxnTp0/H4sWLkZCQgI4dO+Lpp59GSkoKJkyYAEVRUFJSgujoaFitVgwZMgSXXXYZ5s+fj8ceewxxcXFYtmwZDAYDrr/++qAKQERERPoW8syjpKQk3Hrrrfjd736HIUOG4NixY3jmmWcwduxYPPLII6iqqjrv42fPno2pU6fi4YcfxrRp02AwGLB69WqYzWbk5uZi9OjR2Lhxo//4Z599FsOGDcO9996LqVOnoqqqCv/617+QkJAQ6kshIiIiHZKECH4K0vbt2/H+++/jk08+gd1ux4gRIzB16lSMHTsW//vf//DYY48hIyMDL774opZtbhalpdWczBYio1FGfHwkaxki1lE7rKV2WEttsI7aSUiI1Gy9pqCGyZYuXYoPPvgAubm5SE1NxW9+8xtMmTIFHTp08B8zceJEHDx4EP/61780aSgRERFRcwgqDL388su46qqr8Pjjj+Oyyy6DdI5lIDMyMjBnzpxQ2kdERETUrIIKQ19//TViY2NRWFjoD0Ll5eXIzc1FWlqa/7irrrpKm1YSERERNZOgBttkWcZtt92GGTNm+Lft2bMHN9xwA+6++27/ytJERERErV1QYejpp5/G4cOHMXfuXP+2ESNGYOXKldi3bx+WLVumWQOJiIiImlNQYeiLL77A/PnzcfXVV/u3mc1mjB8/HnPnzsXHH3+sWQOJiIiImlNQYai6uhoxMTEN7ktMTERpaWlIjSIiIiJqKUGFoX79+mH9+vUN7nv33XfRp0+fkBpFRERE1FKCuprs97//PWbNmoUpU6ZgwoQJSExMRElJCT7//HP89NNPWLVqldbtJCIiImoWQYWhUaNG4bnnnsOyZcuwbNkyCCEgSRL69u2LlStX4vLLL9e6nURERETNIqgwBABjx47F2LFj4XK5UFZWhujoaNhsNi3bRkRERNTsgg5DgG+hRYfDAVVVUVZWhrKyMv++urfmICIiImqtggpDx48fx4MPPog9e/ac85gDBw4E3SgiIiKilhJUGHr88cdx/Phx3HvvvUhJSYEsa3PXWCIiIqKWFlQY2rlzJ/7617/iuuuu07o9RERERC0qqC6dqKgoxMbGat0WIiIiohYXVBi6/vrr8frrr0MIoXV7iIiIiFpUUMNkERER2LVrFyZMmICMjAxYrdaA/ZIk4cknn9SkgURERETNKagw9N577yE6OhqqqjZ4RZkkSSE3jIiIiKglBBWGvvjiC63bQURERBQWIV0Tr6oqsrKy8NVXX6Gqqipg0UUiIiKitiDoFajff/99LFmyBAUFBZAkCe+88w6effZZmEwmLFmyBGazWct2EhERETWLoHqGNm7ciPnz52PEiBFYunSp/6qyq6++Gl999RVWrlypaSOJiIiImktQPUOrVq3CzTffjIULF0JRFP/2KVOmoLi4GG+99RbmzJmjVRuJiIiImk1QPUPHjh3DhAkTGtyXmZmJ/Pz8kBpFRERE1FKCCkOJiYk4evRog/uOHj2KxMTEkBpFRERE1FKCCkMTJ07EsmXLsGnTJrjdbgC+tYX27duHlStX4tprr9W0kURERETNJag5Q3PmzMGhQ4cwZ84c/x3rZ8yYAbvdjiFDhuAPf/iDpo0kIiIiai5BhSGz2YyXXnoJ3377LbZs2YLy8nJER0dj2LBhGDt2LFegJiIiojYj6HWGAGDUqFEYNWqUVm0hIiIianFBhaHly5df8Jh77703mKcmIiIialGah6GoqCgkJSUxDBEREVGbEFQYysrKqrfNbrdj165dWLhwIR555JGQG0ZERETUEkK6UWtdNpsNY8aMwT333IO///3vWj0tERERUbPSLAzVSk1NPeeCjEREREStTUhXk9UlhEBubi5efPFFdOzYUaunJSIiImpWQYWhtLS0c64lJITgMBkRERG1GUGFoXvuuafBMBQVFYVx48ahW7duobaLiIiIqEUEFYbuu+8+rdtBREREFBZBhaHTp0836fgOHToEcxoiIiKiZhdUGBo/fnyT7j924MCBYE5DRERE1OyCCkP/+Mc/sGDBAvTr1w+TJ09GcnIySktL8cUXX+Djjz/G73//e15RRkRERG1CUGFow4YNGD9+PBYtWhSwfeLEiUhMTMT333/P23EQERFRmxDUootbt27Fdddd1+C+yy+/HLt27QqpUUREREQtJagwFB8fjx9++KHBfd9++y2Sk5NDaRMRERFRiwlqmGzq1Kl47rnn4HA4MH78eCQkJKCoqAgbN27Em2++iUcffVTrdhIRERE1i6DC0N13343KykqsXbsWq1evBuBbeToiIgJz587FzTffrGkjiYiIiJpLUGFIkiQ8+OCDuPvuu/HDDz+gvLwc8fHxGDhwIKKiorRuIxEREVGzCelGrVFRUUhKSgIADBw4EF6vV5NGEREREbWUoMPQ+++/jyVLlqCwsBCSJOHtt9/Gs88+C5PJhCVLlsBsNmvZTiIiIqJmEdTVZBs3bsT8+fMxYsQIPPPMM1BVFQBw9dVX46uvvsLKlSs1bSQRERFRcwmqZ2jVqlW4+eabsXDhQiiK4t8+ZcoUFBcX46233sKcOXO0aiMRERFRswmqZ+jYsWOYMGFCg/syMzORn58fUqOIiIiIWkpQYSgxMRFHjx5tcN/Ro0eRmJgYUqOIiIiIWkpQYWjixIlYtmwZNm3aBLfbDcB3uf2+ffuwcuVKXHvttZo2koiIiKi5BDVnaM6cOTh06BDmzJkDWfblqRkzZsBut2PIkCH4wx/+oGkjiYiIiJpLUGHIbDbjpZdewrfffoutW7eirKwM0dHRGDZsGMaOHQtJkrRuJxEREVGzCCoM/e53v8Ott96KUaNGYdSoUVq3iYiIiKjFBDVnaMeOHTAYDFq3hYiIiKjFBRWGRo0ahbfffhsul0vr9hBRC5FlCRzRJiIKcpjMYrHg448/xubNm9GpU6d6l9JLkoRXXnlFkwYSUfNwexWYDEH9PUREdFEJKgzl5eVh0KBB/q+FEAH7z/6aiFofIQC3V2UgIiLda3QY+uCDDzBmzBjExcXh1Vdf1awBqqpi+fLlePvtt1FRUYFLL70UCxYsQNeuXRvVpj/96U/4/PPP0alTJ83aRKQXLo8Cs1EG/34hIj1r9J+E8+bNw4kTJwK2rVq1CkVFRSE1YOXKlXjzzTfxxBNPYN26dZAkCbNmzfIv5ngup06dwl/+8peQzk2kdx6PyiBERLrX6DB09tCXoij45z//GdJ9yNxuN9asWYP77rsPY8eORVpaGpYuXYr8/Hxs3rz5nI9TVRUPPPAA+vXrF/S5iQhQVAG3V7nwgUREF7Gg5gzVCnVuUFZWFqqrqzFixAj/tpiYGKSnp2PHjh2YNGlSg49btWoVPB4P7r33XmzdujWkNtQycN5EyGpryFqGpqXqqAKQJMCrCERGXJxDZXxPaoe11AbrqB0tr4YNKQyFKi8vDwCQmpoasD0pKQm5ubkNPmbv3r1Ys2YN3nnnnZB6pc4WExOh2XPpHWupjeauo8PlgSJJkCAhMioCJuPF+8OZ70ntsJbaYB1bl7CGIYfDAcB3e4+6LBYLysvL6x1vt9vxpz/9CX/605/QrVs3TcNQRYUDiqJq9nx6ZDDIiImJYC1D1FJ1dHtVVFQ4ISAgFC8izGH9cdAs+J7UDmupDdZRO7GxEf77o4Yq5J9+odyHzGq1AvDNHar9HABcLhciIuqn5ieeeALdunXDzTffHPQ5z0VRVHi9fGNqgbXURnPXUVFUKIqAKgSq7B6YDRfnUBnA96SWWEttsI6h0/LnVZPC0D333FOvF+d3v/sdTCZTwDZJkvDZZ59d8Plqh8cKCgrQpUsX//aCggKkpaXVO379+vUwm83+NY4UxTfx87rrrsPkyZPx2GOPNeXlEFENj1eFx6vCyHkMRKRDjQ5DN954o+YnT0tLQ1RUFLZt2+YPQxUVFdi/fz+mT59e7/hPP/004Os9e/bggQcewAsvvIAePXpo3j4ivVBVAadbQbTt4u0dIiI6l0aHoUWLFml+crPZjOnTp2Px4sVISEhAx44d8fTTTyMlJQUTJkyAoigoKSlBdHQ0rFZrvYUYaydgd+jQod4tQYioaZxuBVERRgC8YRkR6UvY+8Rnz56NqVOn4uGHH8a0adNgMBiwevVqmM1m5ObmYvTo0di4cWO4m0l00VNUFU4P5zAQkf5IgjcSAwCUllZzMluIjEYZ8fGRrGWIWqqOHkVFaYULap0fARazAQnR1ovm/oJ8T2qHtdQG66idhIRIzdZrCnvPEBG1Hl4vr3AhIv1hGCIiP0UVcLi9mq7sSkTU2jEMEVEAp1u5aIbJiIgag2GIiAJwIjUR6Q3DEBEFEAJwOD28wp6IdINhiIjqcXtVOFycO0RE+sAwRET1CAHYnd6Ay+6JiC5WDENE1CCvV4Xdyd4hIrr4MQwRUYMEAIdLYe8QEV30GIaI6JwUVYWLV5YR0UWOYYiIzkkIwOXyQuJYGRFdxBiGiOi83IoKReVQGRFdvIzhbgARtSxVCJzIr4TDpcDl8SIp3gb5PD0/qirg8ngRYeaPC2q9at/XVXYPomwmdEmOPu/7mqgu/nQj0pEDx0vw0dZs5JfYkRgbgfIqJ6JsZozL7IAeneIafIwQvlt02Cwm3qaDWqXa93VeiR2KImAwSEhJsGHSiK7o2y0h3M2jNoDDZEQ6ceB4CV755CByCqtgMRkQGWGC2WRAXokD731zDEdzys75WK9XhVfhRGpqfc5+X8dEmWExGZBTWI1XPjmIA8dLwt1EagMYhoh0QBUCH23NhtPtRVyUBWaTAbIEGI0GxNhMcHlUfLnn9Dkvo1dUwRWpqdVp+H0twWwyIC7KDKdbwUdbs7k8BF0QwxCRDpzIr0ReiR2RVlO9K8MkSYLNYkBRmQO5RdXnfA6nm2sOUetyofd1pNWIvBI7TuRXhqmF1FYwDBHpQJXdA0URMBob/l/eYJChqL5bcJyLoqhwcEVqakUu9L42GmUoikCV3dPCLaO2hmGISAeibCYYDBK83obn/SiKCoMM2KznvqaidkVqTqKm1uJC72uvV4XBICHKZmrhllFbwzBEpANdkqORkmBDtdNbL8wIIWB3KWgXF4HUdpHnfR6vqsLhVtg7RK3Chd7X1U4vUhJs6JIcHaYWUlvBMESkA7IkYdKIrrCaDSircsPtUaAKwOtVUGH3wGKSMS6zwwXXZRECsDu84IVl1Bo0/L4WcHsUlFW5YTUbMGlEV643RBfEMESkE327JWDmNX3QqX0kXB4F1Q4P3B4FKQkRuHF093OuM3Q2j6Ki0u5m7xC1Cme/ryuq3HB5FHRqH4mZ1/ThOkPUKFx0kUhH+nZLQJ+u8U1agbohLrcCh1uB1WRoppYSNV7d9zVXoKZgMAwR6YwsSeiWEgOPoqK0whXU5fKqEKh2eGA1yQD4C4fCr/Z9TRQMDpMRUVC8igqnWwl3M4iIQsYwRERBEQJclZqILgoMQ0QUNI9XwOXhpWVE1LYxDBFR0FQh4HB66t0KgYioLWEYIqKQuL3qOVcAJiJqCxiGiCgkiirgcHPuEBG1XQxDRBQy3tGeiNoyhiEiCpmiqHC4eM8yImqbGIaIKGS+O9p7wc4hImqLGIaISBNeRYXTw0UYiajtYRgiIk0IgZrL7MPdEiKipmEYIiLNeBQuwkhEbQ/DEBFpRlW5CCMRtT0MQ0SkKS7CSERtDcMQEWlKUQXsXISRiNoQhiEi0pzT7YVX4XX2RNQ2MAwRkeYURaDC7g53M4iIGoVhiIiahdujoJqX2hNRG2AMdwOIqGV5FRXr/3cUJZUujM3sgIQYa7OcRwjA7vTCbDLAZODfXUTUevEnFJHO/HC4CJ9sP4kdBwqw6v2fkF9qb7ZzKapAlcMDsHeIiFoxhiEinUlOsPmzSZXDg5c+PIC8kuYLRG6PAqfL22zPT0QUKoYhIp3pnBSF2yf19QeiaocHL324H7nF1c1yPiGAaqcXAry6jIhaJ4YhIh0alZGK2yb19U9utju9eOnD/cgprGqW83m9KqodXq5MTUStEsMQkU4NT0/GzeN7Qa7JJw6XgtUfHkB2XqXm5xIA7C4vXB4OlxFR68MwRKRjmT3b4eYre0Gu6bFxeRS8vPEAjp4u1/xcqipQWe2BytEyImplGIaIdK7/JYn49YReMNR0Ebm9Kl75OAsHT5Rqfi6PoqKi2sWry4ioVWEYIiL07ZaAW6/t418PyKsIvPrJIew9Wqz5uVweBZV2LsZIRK0HwxARAQB6dYrDbyamwWIyAABUIbDu88PYkVWg6XmEABxOL6qdnFBNRK0DwxAR+XVPjcFvr+uLCItvcXoB4L2vfsZXe05reh5V+BZjdPDu9kTUCjAMEVGATu2jcOcv0hFtM/m3bdp2Apu2ZUMI7WY/q6pApd0Nl0dlICKisGIYIqJ6khNsuGtyPyTEWPzbvtqTi3e/+hmKhpeDKYpARbUbHq+q2XMSETUVwxARNSghxoq7JvdDaqLNv23XwUK8/ukhuL2KZufxKirKq9xQFF5zT0ThwTBEROcUbTPjjuvS0S012r8t60Qp1nx0AHanR7PzeBQVZdUuTXudiIgai2GIiM4rwmLEbf9fX6R3i/dvO5Ffhef/8xNKK12ancfjVVFWxUBERC2PYYiILshklDHtqt4Ympbk31ZY5sSqDftwqki7G7x6vOwhIqKWxzBERI1ikCXcMKY7rhrSyb+t0uHBi//5SdPVqj0eFeXVLqgaXrlGRHQ+DENE1GiSJGH84E6Ycvkl/hu8ur0qXv3kILbtz9fsPG4Ph8yIqOUwDBFRkw1JS8Kt16bBbPL9CFEF8P43x/Dx1mzNenTcNT1EDERE1NzCHoZUVcWyZcswZswYZGZm4vbbb0d2dvY5jz98+DDuvPNODB8+HCNHjsTs2bNx+rS2q+MS0YX17hyHO3/RDzF1Fmf8em8u3th8CG6PNpfe+3uIeNk9ETWjsIehlStX4s0338QTTzyBdevWQZIkzJo1C263u96xpaWluO222xAZGYnXXnsNL774IkpLS3HHHXfA5dLuqhYiapwO7SLx+xv6B6xFtP94KV74YD/Kq+v/PxyM2qvMPApXqiai5hHWMOR2u7FmzRrcd999GDt2LNLS0rB06VLk5+dj8+bN9Y7/7LPP4HA48NRTT6FXr17o378/nn76aRw9ehTff/99GF4BEcVGWXDnL/qhT+c4/7bTRdVY+d6PyCmo0uQcnpqFGZ0ehYGIiDRnDOfJs7KyUF1djREjRvi3xcTEID09HTt27MCkSZMCjh85ciRWrFgBi8Vy9lOhvLw8pLYYDGHvJGvzamvIWoampeqoAjAYJEgajEDZIoz4zcQ0fLQlG9/szQUAVNo9eOGDn3DTFT0xsFe7kM8h4Lu5qyRJiLQa0ZipSXxPaoe11AbrqB0t/zAKaxjKy8sDAKSmpgZsT0pKQm5ubr3jO3XqhE6dOgVse/7552GxWDB06NCQ2hITExHS4+kM1lIbzV1Hh8sDRZIaFSoaa/rEdHRJjcGbmw9BVQW8isC/PzuMkio3rr+8B2Q59J9ekgQIgwExNnOjf6HwPakd1lIbrGPrEtYw5HA4AABmszlgu8ViaVRPz7/+9S+88cYb+POf/4zExMSQ2lJR4YCi8GaRoTAYZMTERLCWIWqpOrq9KioqnJqv55N5SQKiruuL1z45BLvLCwD4ZGs2sk+X4+areiHCEvqPnTIAxSYDoiPMsJjlcwY6vie1w1pqg3XUTmxsBGRZmx62sIYhq9UKwDd3qPZzAHC5XIiIOHdqFkLgn//8J5577jncdddd+M1vfhNyWxRFhZd3ztYEa6mN5q6joqhQFNEsixt2S4nB72/sj1c/OYiCUt8fPVknyvDs+h8xfUJvJCfYLvAMF+ZQvHC7FURFmBAZYYJ6nkvw+Z7UDmupDdYxdFr+6ArroGXt8FhBQUHA9oKCAqSkpDT4GI/HgwceeACrVq3CvHnzMHfu3GZvJxE1XWKMFb+/vn/APc2Ky514bsM+/PhzsSbnUFSBSrsHFdVugBOriShIYQ1DaWlpiIqKwrZt2/zbKioqsH//fgwZMqTBx8ybNw+bNm3CkiVL8Nvf/ralmkpEQbCYDbhlQm9ceemZuX5ur4p/f3YYH2/N1mRBRVUIVDs9qGQgIqIghXWYzGw2Y/r06Vi8eDESEhLQsWNHPP3000hJScGECROgKApKSkoQHR0Nq9WKd999Fxs3bsS8efMwbNgwFBYW+p+r9hgial1kScKVl3ZCx3aRWPfFEbhqFmT8em8ucgqrcPOVvRBtM1/gWc5PCMDu8kIAiIk0A1yjkYiaIOzX9s2ePRtTp07Fww8/jGnTpsFgMGD16tUwm83Izc3F6NGjsXHjRgDAhx9+CAD4+9//jtGjRwd81B5DRK1TWtd43DOlP5Ljz8wHPJZbieXrf8TPpytCfn4hAIfLi/IqNwTTEBE1gSQEbw0NAKWl1ZzMFiKjUUZ8fCRrGaKWqqNHUVFa0fJ3h3d7FLz39c/Yc+TMvCFJAiYM6YzLB3aArMHiIWaTjGibGTarEXFxfE9qgf9/a4N11E5CQqRm6zWFdZiMiPTHbDJg6rgeiI0045u9uVCFr1fn0x0ncSy3Ajdd0RM2qxG5RdWwO72wWY1IbRfZpJDk9qgoq3TBq6iIim6+9VxUIXAivxJVdg+ibCZ0SY7WJMy1Nl5VxbYf82F3qbBZZFzapz2MGl3STNQaMAwRUYs6mlOGL/ecRlGZAxazAU634r9E9nBOOf7x1h7ERpthd3igqIBBBtrFRWBcZgf06BTX6PMoqm/F6qIyB4RXgcmg7QKTB46X4KOt2cgrsUNRBAwGCSkJNkwa0RV9uyVod6Iw27QtGx9tyYajZk6WBOA1ixGTRnbFtcO7hrt5RJpgtCeiFnM0pwzvfXMMeSV2mE0GxEZZkBBtgaHOytR2lxe5RXZ4vCoiI4wwmwzIK3HgvW+O4WhOWZPPqagC5VUulFW5oQqhyRL+B46X4JVPDiKnsAoWkwExUWZYTAbkFFbjlU8O4sDxktBP0gps2paN9f/7GdVOL2RZgtEgQZYlVDu9WP+/n7FpW3a4m0ikCYYhImoRqhD4cs9puDwKYmxmmIwGyJIEi9mI9nFWmI2BP47sLgUlFS5IkoQYmwkuj4ov95wOao6TKgQcLi9KKlxwuLwhBSJVCHy0NRtOtxdxURaYTb7XYTYZEBdlhtOt4KOt2S0+F0trXlXFR1t8yx+YakKQLMmQZQkmgwRFFfhoSza8Kue9UNvHMERELSK3qBpFZQ7YLEZIZ6URWZYRafWN2te9fZnHq6KwzAGHy4sIs4yiMgdyi6qDboNXUVFR7UFZlRuKKuq1ozFO5Fcir8SOSKup3uNrbyKbV2LHifzKoNvZGmzfnw+HywujLDX4Oo2yBIfLi+3788PUQiLtMAwRUYuwO72+OUDnuvqj5vdtbKQp4P5lQgBlVW5UOrzwKgJ2pzekdvh7iSqdqKh2NTkUVdk9UBQBo7Hh12E0ylAUgSq7J6R2hltJhdO3QMG5SiP5lnMqqXC2XKOImgnDEBG1CJvVCIOMc9+cUtT+3pUQH21BfLQlYDjL5VZgd3lRUObQpD2KIlDtPBOKvI0MRVE2EwwG6ZyXRXu9KgwGCVE2kybtDJeEGKvv+3Gu0b6a71dCDBe7pbaPYYiIWkRqu0i0i4uA3aXg7OXNhBBwe1VYzQa4FRVCCERYjEiKi4DZJNc5DvhoSzbe++pnuNyKJu2qDUWljewp6pIcjZQEG6qd3gZfR7XTi5QEG7okR2vSvnAZlp6MCIsRXlU0+Dq9qu97NCw9OUwtJNIOwxARtQhZkjAuswMsJhkVdg88XgWqEPB4FVTYPbCYDRg3qAMsJoN/vyT7Jk+fPbl6R1YB/vnOHhzJKdesfWf3FClKw6FIliRMGtEVVrMBZVVuuD2+1+H2KCircsNqNmDSiK5tfr0hoyxj0siuMMgSPIqAqgqoQoWqCngUAYMsYdLIrlxviC4KfBcTUYvp0SkON47ujpSECLg9CqrsHrg9ClISInDj6O4Yk9mx3n6PV0Xn5ChMGdMdHdpF+p+rrMqNNRsP4L2vfobTHdo8orrqhqKyKhc8ilrv6rO+3RIw85o+6NQ+Ei6PgooqN1weBZ3aR2LmNX0umnWGrh3eFf839hJEWo1QVQFvTSiKtBrxf2Mv4TpDdNHg7ThqcGn00HGZeW1c7LfjAHyTmM+3wvS59iuqii93n8Z/vz8V0O7YSDNuGNMdfbrEB5zHYJAQG2tDeblvYcRgyLIEk1GG1WyAxWSAQZb9w0Z6WoF618FCrkCtAf6c1I6Wt+NgGKrBN2bo+D+5NvQQhkJ1uqga6/93FLnF9oDtA3ok4rrLuiEqwjd5WYswVEuSfMHIYjQgwmqE2ShruqJ1a8f/v7XBOmpHyzDEaE9EbU6HdpG4+8b+uGpIp4DVq/ceLcbSt/Zg18GCepN+QyWEbwjN7vKitNKFkgonnB7fJO6LsDOISFcYhoioTTLIMsYP7oR7p2Sgc1KUf7vD5btVxIsf7kd+if08zxA8VRVweVSUV7lQXO5EtcNTc6sPpiKitohhiIjatOQEG+6a3A+/uKxbwGX4x3Mr8Y+39+K9L4/A7dHmMvyzCeEbbqywe1Bc4UT5OSZcE1HrxjBERG2eLEsY2T8F99+UifRuZyZRq6rAJ1uzsfjNH/Djz8WaD53VVXcIrbSyNhQxFRG1BQxDRHTRiI2yYPrVfTDj6t6IizL7t5dXufHvzw7j5Y1ZyC9tnqGzWqoq4HQrKK10oazKBbeXPUVErR3DEBFddPp2S8CcX2biikEdAyZYHzlVjmff2YsPvjsOh0u7tYkaoqq+e6CVVdVMtnZ7AQgGI6JWiGGIiC5KZqMB147ogkd+Oxy9Osf6t6sC2LIvD4vf/AFb9uVBUZv38mb/ZOtqN4rLXai0e+CpuT8bh9GIWgeGISK6qKUkRuK3k/pi+tW9kRBt8W93uLz44LvjWPbOXhw4XtKs84mAM5OtqxwelFa6UFTuW+Ha4fZe8H5oRNS8jOFuABFRc5MkCendEtCrUxy+/TEXX/5wCm6Pr3emsMyJVz89hO6p0bh2eNeAy/Sbi6oKqBDwKiqc7ppVrg2+Va5NRgOMBl8w0tOijkThxDBERLphMsoYN6gjLu3THpt35mBXVgFq88ax3Eo8t2Ef+ndPwIShndE+LqJF2lS7mKOiKHB5FEiSBIMswWyUYTIZYDHJkCW52XuuiPSMYYiIdCfaZsaUyy/ByH7J2LTtBA7nlPv37TtWgv3HSzC4TxLGD+6IuCjLeZ5JW0IAQvhuhurxqpBcXsiyBKvJAKtFf7cAIWopDENEpFupiZG4bWJfHM4pwyfbTuB0zb3OVAHszCrAD4cLMbRvMsYN7IBom/kCz6a92l6jasULh1uBySDBajbCbOJQGpGWGIaISPd6dYpDj46x+PFoMTbvPImSChcAwKsIbNmXh50HCjA8PRljMlPDEoqAmqvSVAG3xw1ZlmAw+HqMLGYjTP4eI8FwRBQEhiEiIgCyJCGzZzv0vyQBO7MK8d/dp1BR7Qbguwrsmx9zsW1/Pob1TcKYzA6IiQxPKBIAFFVAUQXcHhUGp28ozSBLMBhkmIyyf54RwxFR4zAMERHVYZBlDE9PxuDe7bH9QD7+98NpVDk8AHyh6Nt9edh2IB+X9knCmAGpSIixhrW9tcHI10IFkuS7Os0oSzCZDDAZZRglGbIBDEhE58AwRETUAJNRxqiMVAztm4Tt+wvw1Z4zocirCGzbn48dB/IxoEc7XD6wA1ISbGFusc+Zq9N8iz1Kkm9pAVnyDa2ZjTLMNSFJAuccEQEMQ0RE52U2GjB6QCqGpSdhx4ECfL3nNCrsvlCkCuCHI0X44UgR+nSJw5gBHdA9NbpVLaDov0INAl4FcLkVyLLXd/m+yQCTQYLRKMNkkP3HE+kNwxARUSOYjQaMykjFsL7J+P5QIb7ecxollS7//oMnynDwRBk6tovEqAGpyLgkAQa5dS7yr6p1Lt+v6TkyyBIsJgPMJhlGg+9DCA6pkT4wDBERNYHJ6JtTNCQtCT/+XIyv95xGbs0l+QBwqqgab31xBJu2mTEiPRlD+yYh0moKY4vPr97aRk74J2QbDb4eI4NBhiT7jvOFJ/Yg0cWFYYiIKAgGWcLAnu2Q2SMRh3PK8c3eXBw5dWbxxopqNz7dcRJffJ+DzJ7tMKJfCjq2iwxjixun7pyj2luWyJIEo1GCFzIcdhdkAAajDKMsQ5YBCRJkWfKvks2gRG0NwxARUQgkSULvznHo3TkOp4uq8e2Pudh7tBiK6ksEXkVg18FC7DpYiM5JURiRnoz+lyTCZGydQ2gNUYWAovquXHO4vFAU4R9ekyRfGJJkwFjTm2QwyL5L/WVfSJJ8U7UZkqjVYhgiItJIh3aRuOmKnrh2eBds25+P7QcK/FegAcDJgiqcLKjCh1uyMbh3Owztm4ykFroHmtZqh9dqvgIUBFzeL0HyByZjzSRtY81aSLVBCTgz0dw39MbAROHBMEREpLFomxlXDemMcYM64sefi7H1p3ycLKjy73e4vPj2xzx8+2MeuqZEY0if9si4JBFmkyGMrdaOEICA8K0QWXMVG9xKQG+SLPnWQoLkm9At4JukbjKdWRfJIJ9ZWbv2eYmaA8MQEVEzMRpkDOrVHoN6tcepomps+ykPe44Ww+NV/cdk51UiO68SH3x3HAMuScTgPu3RNbl1XZ6vlbq9SQpqF4o8w+3xTeA+OzDJBhlyzdeyoWb4rWbtJEi1E8ADn6t2ojeE7/lQ8/jadtS0iAGLADAMERG1iI7tIjFlbA9MHNkVuw8XYWdWQcBVaG6Pip0HC7HzYCESYiwY1Ks9BvZqh8Qwr3Dd0hoOTIp//9lzlQBfL9TZoUby/+fMcWcej5pFKGX/vCaj4UzAkiTfPKmzA1bt4+vOg6ptM7VtDENERC3IajZiZL8UjEhPxqmiauzMKsCeI8Vwec78wi+pcOHzXTn4fFcOOidFIbNnO2RckhC2m8S2JvXmKjXuUefY7qt5YwOW/1hIkGX47wdXOwfKIPmWIJCl2qvrzjymtu0X6vCreyy1HIYhIqIwkCQJndpHoVP7KEwc2RX7j5Xi+0OFOHqqPOBXd+2k64+2HEePDrEY0CMR6d3iYWvFaxe1NU0PWKJOZ5Xi74UKuLrO34MEKDUTw2tX/RYGA+wuL1RV+I9XhfCv96QKIDqC39+WxDBERBRmZqMBA3u1w8Be7VBW5cKeI0XYfbgIBaUO/zFCAEdOlePIqXJs+FpCj44x6H9JIvp2jUcUf3GGVe0c7wsFKi8Ar6LCZDGhvMoFVT3ruJo550aDjKgIY8DwHjUvhiEinZIkwGKqWeumzl+ztTsl/3yLMxdAq6j5K1qt+SsWvrkVomZuhX9oQZx7mIHOLy7KgrEDO+LyzA7ILbZjz5Ei7D1ajPJqt/8YVQgczinH4ZxybJCAbikx6Nc9HundEhAXZQlj66kp+P9H68EwRKRTJoOMuOjaybnn/6l89tyHmq/824QQNd38gFABFQKKGtjtr6oqIHyBCsK3rc5V0/7wxLVmfCRJQod2kejQLhLXDO+CE/mV2Hu0GD/9XILKOmsXCQEcy63AsdwKfPhdNjok2pDWNR59uyWgQ6LtorwqjUhrDENEOlV3/ZamPcb/VcA2f8+SATBAQu2SOQ0FqLrba8NPbQ+TIlQoioBXFVAVFYqoCVZqzbFNavHFQZYkdEuJQbeUGFw3shuy8yux7+cS/HSsGBX2wAvUTxfbcbrYji++P4UYmwm9u8SjT+c49OwYC4v54ljHiEhrDENE1KwaClD1t9esAWMADDAAxtqw5AtYiqrCqwgoiurrZRICqiJ8QUkRAT1TFztZltA9NQbdU2Mw6bKuyCmowv7jJfjpeCmKy50Bx1bYPdiZVYCdWQUwyBK6pkSjd6c49Ooci+QEm3/dHSK9YxgiolaptudKCF+vk8kgwWTwzXE68ztcgoAvENUGJo9Xhduj+O8NdjGTJQldkqPRJTka1wzrgsJyJw4cL0FWdhlOFFQGhENFFfj5dAV+Pl2BTduByAgTenaMQc+OsejZMRaxnGtEOsYwRERtTt0VhIGa9V5kA8xGQLJK8Coq7C4vvIp6wXVdLhaSJCEpLgJJAzti7MCOqHZ6cOhEGQ6eLMPhnDI4XErA8dUOD/YcKcaeI8UAgMRYKy5JjcElHXwfXNOI9IRhiIguKkIIGGQJMTYTDEYZNpsFwqug2u6GRxH1L2e+SEVaTRjUuz0G9W4PRRU4VViFgyfLcCSnHDmFVfWGFIvLnSgud2JHVgEAoF2s1T8c1zUlGomx7DmiixfDEBFdlGov8beYjbBZjLCaDPB4VFQ7PXDpZBitlkE+M5w2YUhn2J1eHD1djqOnynEkpxwlla56jykqd6KoTjiKjTSjV5c4dEi0oXP7KKQk2mCQ5ZZ+KUTNgmGIiHRBVQUMBgmxURZ4vAoqHR64PYouJl2fzWY1IuOSRGRckggAKK104efT5f45RXXXNKpVXu3GzgMF/q9NBhkdkyLRJSkKnZKi0TkpCrGRHFqjtolhiIh0RQgBo0FGfLQF1Q4P7C4vFEWHiaiO+GgLLu2ThEv7JEEIgdJKV83aRZXIzqtEcYWz3mM8iorjuZU4nlsJIBcAEG0zoVP7KHRsH4mONWskce4RtQUMQ0SkTwKIijDBajbC4fLC6Vb86xjpZV5RQyRJQkKMFQkxVlzaJwkAUGF342RBFfJKHTiUXYrTRdUNDjNW2j04kF2KA9ml/m0xkWZ0SLQhtV0kOiRGIjXRhvhoCxeDpFaFYYiIdKv25pnRNhOiIky+9YuEgNutwOFW/Kto1x6rVzE2Mwb0SMSYWBvKy+1wuhScLqquuYlsJU4WVKGsqv7QGgBUVLtRUe1G1oky/zaLyYCUBBtSEm1ISbAhOSECyfE2RFj4K4nCg+88ItK92qAjSxJkSYIpQobNaoJXVaGqvluJuNwK3IrKW4YAMBlldE2JRteUaACpAIAqhwc5hVU4VViNnMIqnC6sDrhtSF0uj4Ls/Epk51cGbI+JNCM5PsK3REB8BNrXfG6z8ka01LwYhoiIzlIbdIyyDMgAICPCYoIqVCiqgKoCHq/C+UZ1REWYkNYlHmld4v3bKuxunC6qxumiauQW2ZFbXN3glWv+42t6kQ7nlAdsj7Qa0S4uAu3jItAu1or2sVYkxkYgIcYCo4FXtFHoGIaIiBpBCAEJEoyyBMiAxSQjwmJEtcMDp0dhKGpAjM2MmC7mgIDkdHuRX+JAbkk18ortyC9xIL/UDqdbOefzVDu9qM7zTeauS5KAuCgLEmOsSIixIDHWisQYK+KjLUiItvJebNRoDENEREEQwjesFhNpQaSqwuHy+tcvUlV9D6Odj9VsrDPE5iOEQEW1G/mlDhSUOlBQakdBme/z84UkIXzLApRWuoBT9ffbrEYkRFsQX/MRF21BfJQFsVG+fxmWqBbDEBFRCIQQkCUJUREmRNvMvvujKQo8HhUuj+K/gaze5xmdjyT51n+KjbKgd+c4/3YhBKocHhSWOVBY5kRRuQNFZb7FIEsrnbjQRX92pxd2pxc5hdUN7reaDYiLsiA2yozYSDPioiyIifR9Xvuv2cTApAcMQ0REGqgNPLIEWIwGWE0GRMN3hZpXEVCUmhvJKqr/KjWGpPOTJAnRNjOibWZc0iE2YJ+iqiitcKG4wlnz4UJJhRMlFU6UVrrgbcSwpdOtIK/EjrwS+zmPsZgMiIk0I9pmQozN92+0/18TomxmREeYYDUbuFxAG8YwRETUDGoDjgQJJoMEk0GuuWmsBMB3Cb+iCni9KjxeFS6P6ru0X8drHDWFQZbRLi4C7eIi6u1ThUCl3YPSSidKK1worXL5h9PKKl0or3Y3+nYsLo9S0zPluEB7fL2DtR+RESZERRgRafV9Hmn1fR4daYIlwuxfsoFaB4YhIqIW4vv95/slWDsZ22iWEWHxbfUqKjweFW6P4ruprBDsOQqCLEmIrRnm6pZSf79aM/xWXuVCWZUb5VVulFf5QlJ5zRVtlXb3BYfh6lJU4X98YxhkCbaagBRhMcJm9d1Dz2Y1ItpmxpgBKWgfZ2t8AygkDENERGFW91J+o8W3xpEQAooi4FFVX++RokJRBASE/ya0/s+pSWRJ8l3pZjOjc1LDx6iqQJXTg8qacFRh96DC7kal3YPKmn+rHB5U2T1Qg/gmKKqoea6G12L6as9p/O13I2HhnKUWwTBERNTK1A6hGAwSJFlGYakddqcXkTYTUhMjARVQalfL9irY8mMuisqdiI4wYWh6Mgyy7O9RUoVAblE17E4vbFYjUttFQq4zt+VC+wHf/Jw9R0vgdKuwmmX0757QpDvWa9GGUM+hqCr2HilCWZUbcVFmDOjZ7ryvQZbPBKaO7c99DgBwuryorAlLpwp9q3ErqoAk+SZxVzk8sDu9qHZ64HCd++q4uqocHhzPrUCvznFNrgU1nSTCPHCpqiqWL1+Ot99+GxUVFbj00kuxYMECdO3atcHjS0tL8cQTT+Crr74CAFx77bX485//DJsttO7E0tJqeL1qSM+hd0ajjPj4SNYyRKyjdtp6LQ8cL8FHW7ORV2KHoggYDBJSEmyYNKIr+nZLwGufZmHLT/mQJAlGgwyjQYbJKGFY3yRMHNkNR0+V49u9ecgvs8PjUSGEQGyUGaP6paB7h1gcySnDl3tOo6jMAUUFDDLQLi4C4zI7oEenOADA13tO4cvdp333boNvxpPVbMC4QR0wJrPjBV/D0Quc40L7G+NCzxHqa9D0dUiA0WxCbkEFquy+kGR3eWF3enC6qBrH8yqhKCqibWZ4FDXg+02BEhIiYdBo0c2wh6Hly5fjjTfewKJFi5CcnIynn34aJ0+exIcffgizuf7djmfMmAGXy4UFCxagoqICDz30EIYOHYq//e1vIbWjrf6wbE3a+i+e1oJ11E5bruWB4yV45ZODcLq9iLSaYDTK8HpVVDu9sJoNSEmIwN6jJed8/KBeiSgqd8GjKIiJtMBskgEVcHlVmI0yBvduh12HiuBye2E2GSDLErxeFXaXAqNBwrVDOyO3pBqf78yBr3S++UuK6rs6TgJwzbDO5w0TR3PK8N43x+DyKLBZjDAYZCiK7xwWk4zhfZOw7UDBOfffOLr7BQPRhc7Ro0MMdh0shCoAg4Ta+etQBCBLF34NWr8Og0FCbM093uou1Fn3HCkJNtisRhSVOf3f75nX9GEgOouWYSis65i73W6sWbMG9913H8aOHYu0tDQsXboU+fn52Lx5c73jd+/eje3bt2PRokXo168fRo4cicceewzvv/8+8vPzw/AKiIi0pwqBj7Zmw+n2Ii7K4gsrkgSzyYC4KDOqHe6AICTV+ai1+3AxSirsMBsNcLkVVFZ7UOnwwOX24nRRFT787jiKyuwwGmT/MIzJKCM20gRVFdh9uAA/HC5GVKQF7WMtaB8fgeSESKS2i0Tn9pFITozEgRNlMJlkWMwGmIyyv3fKN7wHbD9YAINBQlJcBKIizLCYDIiwGBEfZYbHq+J/P5yGy6MgxmaGyeh7jSajATE2E1weFV/uOX3e+TiqEPhyz7mfw+nyYmdNEDLKvqEvWZIgyxKMMqAK4Mvdp6Go5w7KFzqHy63gy93avg6jIfD77XQr+GhrdlBzk6hxwtoztHfvXtx0003YtGkTunfv7t8+bdo09OnTBwsXLgw4/sUXX8Qrr7yCb775xr/N7XYjMzMTS5YswcSJE4Nui6qqnIgYIkkCZFlmLUPEOmqnrdbSq6gor3b7Ak4D80Vq1ym6EEmSYJDrP762h6c2HDS0XxVn5i6dOUTC2UdH2UwNTvL1KiqqHB5IUv3HAPAvIyDLUoOvsXbOU2SEyXcLlNrtdY5RVBXVNedoiKr6Xod/RYOzn6Dmy0ir8ZyLKyr+11H3SQLbecHXAd+922rvoybLUsASCv5awfc9kyTfv7XH1D5HbKSZ92Kr41w1D0ZYJ1Dn5eUBAFJTUwO2JyUlITc3t97x+fn59Y41m82Ii4tr8PimkJswGZDOj7XUBuuonbZWy9oFA2t/MYaiocefWQOp4f2ABIg6vSXizCfirE2qKhocqvAqNYFLOncbVAFIApD8M3nq/iNBwBcyTGcFldqQJrw1R54jcEEGaq+P9//SDDjQF7gkSYLReNZrEGfO5Q+NDb0OVQKgwiDLkGTUC1uABFUIGGQpIMjIhjoBr2b/2cG0bgj1vWZJs2EhChTWMORw+BaxOntukMViQXl5eYPHNzSPyGKxwOU6952QiYjaEovZgCRz215jpjlfQ22wsZgMsDSw6KI2J/H9YzEbYTE376/KljgHnV9YI6bVagXgG+qqy+VyISKi/hvcarXWO7b2+FCvJiMiIiJ9CmsYqh3yKigoCNheUFCAlJT6y4ampKTUO9btdqOsrAzJycnN11AiIiK6aIU1DKWlpSEqKgrbtm3zb6uoqMD+/fsxZMiQescPHToUeXl5yM7O9m+rfezgwYObv8FERER00QnrIKXZbMb06dOxePFiJCQkoGPHjnj66aeRkpKCCRMmQFEUlJSUIDo6GlarFZmZmRg8eDDuv/9+LFy4EHa7HQsWLMANN9zAniEiIiIKStgXXVQUBc888wzeffddOJ1ODB06FI8++ig6deqEnJwcXHnllVi0aBGmTJkCACguLsZf/vIXfP3117BYLP4VqC0WSzhfBhEREbVRYQ9DREREROHEBQuIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXdBOGiouL8cADD2DEiBEYNGgQ7rzzThw5csS//8CBA5g+fToGDhyIcePGYfXq1WFsbdtw7NgxDBo0CO+++65/G+vYNKdOnUKfPn3qfbz99tsAWM+m2LBhAyZOnIiMjAxMmjQJH3/8sX8f69g427Zta/D92KdPH1x55ZUAWMvG8ng8WLp0KcaNG4dBgwbhlltuwffff+/fzzo2XnV1NR5//HGMHTsWl156Ke6++26cOHHCv1+TWgqduOmmm8SvfvUrsXfvXnHkyBFx3333iVGjRgm73S5KSkrE8OHDxUMPPSSOHDki3nnnHZGRkSHeeeedcDe71XK73WLKlCmid+/eYv369UIIwToG4fPPPxcZGRkiPz9fFBQU+D8cDgfr2QQbNmwQffv2FWvXrhXHjx8Xy5cvF2lpaeL7779nHZvA5XIFvA8LCgrEN998I9LT08Vbb73FWjbBP//5TzFq1Cjx9ddfi+PHj4uHHnpIDB48WOTl5bGOTXTHHXeIMWPGiC+++EIcOXJEPPzww+Kyyy4TJSUlmtVSF2GopKRE3H///eLQoUP+bQcOHBC9e/cWe/bsEatWrRJjxowRHo/Hv3/JkiXimmuuCUdz24QlS5aIGTNmBIQh1rHpnnvuOTF58uQG97GejaOqqrjiiivEU089FbD99ttvF6tWrWIdQ+B2u8WkSZPEnDlzhBB8TzbF5MmTxaJFi/xfV1ZWit69e4tNmzaxjk1Q+7v6yy+/9G9TFEVcffXVYvny5ZrVUhfDZPHx8XjmmWfQq1cvAEBRURFWr16NlJQU9OzZEzt37sTQoUNhNJ65VduIESNw7NgxFBcXh6vZrdaOHTuwbt06/O1vfwvYzjo23cGDB9GzZ88G97GejfPzzz/j1KlT+MUvfhGwffXq1bjrrrtYxxC8/vrryM3NxZ///GcAfE82RVxcHP773/8iJycHiqJg3bp1MJvN6Nu3L+vYBMeOHQOAgJu3y7KMtLQ07NixQ7Na6iIM1fXII49g1KhR2LRpE/7617/CZrMhLy8PKSkpAcclJSUBAE6fPh2OZrZaFRUVmDdvHh5++GGkpqYG7GMdm+7QoUMoLi7GLbfcgssuuwzTpk3D119/DYD1bKzjx48DAOx2O377299i5MiRuOmmm/DFF18AYB2D5XK5sGrVKsycOdNfL9ay8R566CEYjUZceeWVyMjIwNKlS/GPf/wDXbp0YR2boH379gB87726Tp06heLiYs1qqbswNHPmTKxfvx6TJ0/GPffcg59++glOpxNmsznguNobv7pcrnA0s9VauHAhBg4cWO+vcACsYxO53W4cP34cVVVVmDNnDl544QVkZGRg1qxZ2LJlC+vZSFVVVQCA+fPn47rrrsOaNWswatQo3H333axjCN5//324XC7MmDHDv421bLyjR48iJiYGK1aswLp16zBlyhTMnz8fWVlZrGMTZGZmokePHliwYAFyc3Phdruxdu1aHDhwAG63W7NaGi98yMWldkji8ccfxw8//IDXXnsNVqsVbrc74LjaItpsthZvY2u1YcMG7Ny5Ex988EGD+1nHpjGbzdixYweMRqP/f+b+/fvj6NGjWL16NevZSCaTCQDw29/+FjfeeCMAoG/fvti/fz9efvll1jFIGzZswNVXX434+Hj/NtaycU6dOoUHHngAa9eu9Q/vZGRk4MiRI3j22WdZxyYwmUxYsWIFHnzwQYwbNw5GoxHjxo3D1KlTsW/fPrjdbk1qqYueoeLiYnz44YdQFMW/TZZl9OjRAwUFBUhJSUFBQUHAY2q/Tk5ObtG2tmbr169HcXGx/1LRQYMGAQAWLFiASZMmsY5BsNls9f6q6d27N/Lz81nPRqrtIu/du3fA9p49eyInJ4d1DEJJSQl2796NiRMnBmxnLRtn79698Hg8yMjICNiemZmJ48ePs45N1L17d6xbtw7bt2/Hli1bsGLFCpSVlaFbt26a1VIXYaigoAB//OMfsX37dv82j8eD/fv3o0ePHhg6dCh27doVEJa2bNmC7t27IzExMRxNbpUWL16MjRs3YsOGDf4PAJg9ezZeeOEF1rGJsrKyMGjQIOzcuTNg+759+9CzZ0/Ws5HS09MRGRmJPXv2BGw/dOgQunTpwjoG4fvvv4ckSRg2bFjAdtaycWrnUx48eDBg+6FDh9C1a1fWsQmqqqowffp07Nu3D7GxsYiJiUFlZSW+++47jBkzRrtaanb9Wyumqqq4/fbbxTXXXCN27NghDh48KO6//34xdOhQcerUKVFUVCSGDh0q5s+fLw4fPizWr18vMjIyxLvvvhvuprd6dS+tZx2bRlEUcdNNN4nrrrtO7NixQxw5ckQ8+eSTon///iIrK4v1bIIVK1aIQYMGiQ8++EBkZ2eLlStXirS0NLF161bWMQjPPvusuPrqq+ttZy0bR1EUccstt4hrr71WbNmyRRw7dkwsXbpU9O3bV+zevZt1bKLp06eLadOmiaysLHHgwAFxyy23iMmTJwuPx6NZLXURhoQQoqKiQixYsECMGjVKDBgwQNx+++0B6w7t2bNH/PKXvxT9+/cXV1xxhXj11VfD2Nq2o24YEoJ1bKri4mLx5z//WYwaNUpkZGSIX/3qV2LHjh3+/axn461Zs0aMHz9e9OvXT0yePFls3rzZv491bJoFCxaIX/7ylw3uYy0bp6ysTCxcuFCMGzdODBo0SPzqV78S27Zt8+9nHRsvPz9f3HfffWLIkCFi2LBhYv78+aK4uNi/X4taSkII0Qw9W0RERERtgi7mDBERERGdC8MQERER6RrDEBEREekawxARERHpGsMQERER6RrDEBEREekawxARERHpGsMQEVEjcVk2oouT7u5aT0RN9+CDD+K999477zEdO3bEF1980UItanmff/45PvnkE/z9738Pd1OISGNcgZqILujEiRMoKSnxf71y5Urs378fy5cv928zm81IT08PR/NaxIwZMwAAr776aphbQkRaY88QEV1Qly5d0KVLF//XCQkJMJvNGDhwYPgaRUSkEc4ZIiJNHDp0CHfddRcGDx6MwYMH45577sHJkyf9+7dt24Y+ffpgy5YtmDFjBgYMGIBx48bh7bffRkFBAe69914MGjQIY8eOxdq1a+s97ptvvsGvf/1rDBgwABMmTMBrr70WcH5VVfHCCy9gwoQJ6N+/P6655pp6vTgzZszAn/70J8yePRuDBw/GnXfeCQDIycnBvHnzMHr0aPTr1w8jR47EvHnzUFpa6n/c9u3bsX37dvTp0wfbtm3zt2vbtm31zlHbiwQA48ePx5NPPomZM2di8ODBePTRRwEAZWVlePTRR3HZZZchIyMDv/zlL7Fly5bQvxFE1GQMQ0QUsmPHjuHmm29GcXExnnrqKfz1r3/FyZMnMW3aNBQXFwccO3fuXIwfPx6rVq1Ct27dsGDBAtx6663o3bs3li1bhn79+mHRokXYu3dvwOPuv/9+pKenY8WKFRg1ahQef/zxgLCzcOFCLFu2DJMnT8aqVatw7bXX4sknn8SKFSsCnufjjz+GyWTCihUrcOutt8LhcODWW2/F0aNHsWDBAqxevRrTp0/Hhx9+iGeeeQYAsGDBAqSnpyM9PR3r1q1Dv379mlSf119/HX369MGzzz6L66+/Hi6XCzNnzsTnn3+O+++/H8uXL0dKSgruuOMOBiKiMOAwGRGFbPny5bBarVi7di2ioqIAACNHjsRVV12Fl156CfPnz/cf+3//93+47bbbAAA2mw2/+tWvMGDAAMyePRsA0L9/f3z++ef4/vvvMWDAAP/jrrrqKjz00EMAgDFjxqCgoADPPfccfv3rXyM7OxtvvfUW5s6d6+/tGT16NCRJwvPPP49bbrkF8fHxAABZlvH444/DZrMBAA4cOICUlBQ89dRT/qHAESNG4Mcff8T27dsBAD179vS/rmCGBpOSkvDggw9Cln1/f7711lvIysrCW2+9hczMTADA5ZdfjhkzZmDx4sVYv359k89BRMFjzxARhWzr1q0YPnw4rFYrvF4vvF4voqKiMGTIEHz33XcBxw4aNMj/ebt27QDAHwgA+ENLZWVlwOOuv/76gK+vvvpqFBcX49ixY9i6dSuEEBg/frz//F6vF+PHj4fL5cKuXbv8j+vUqZM/CAFA37598cYbb6BTp044efIkvv76a6xZswY///wzPB5PiJXx6dGjhz8IAcCWLVvQvn179OvXz99WRVFwxRVXYN++fSgvL9fkvETUOOwZIqKQlZWVYePGjdi4cWO9fQkJCQFf1/aw1BUREXHBcyQlJQV8nZiYCACoqKhAWVkZAGDSpEkNPjY/P9//eW0Aq+vll1/G888/j9LSUrRr1w79+vVDREREvUAWrLPPWVZWhsLCwnMOtxUWFiI2NlaTcxPRhTEMEVHIoqOjcdlll/mHv+oyGrX5MVMbeGrVzkVKTExETEwMAOCVV15BZGRkvcd26NDhnM/7wQcf4KmnnsIf//hHTJ061R/e/vCHP+DHH3885+MkSQLgm7hdV3V1dYNtqCs6OhrdunXD4sWLG9zfqVOn8z6eiLTFYTIiCtmwYcNw5MgR9O3bFxkZGcjIyED//v2xdu1abN68WZNznL2g46ZNm9CxY0d06dIFQ4cOBQCUlpb6z5+RkYGysjL84x//qBek6tq1axeio6Nx5513+oNQdXU1du3aFRB06g5zAWd6uHJzc/3bysvLcfTo0Qu+lmHDhiE3NxeJiYkB7d2yZQteeuklGAyGCz4HEWmHPUNEFLK7774bN998M+666y5MmzYNFosF69atw2effYZly5Zpco61a9fCarVi4MCB+PTTT/Hf//4XS5YsAQD07t0bkydPxiOPPIJTp06hf//+OHbsGJYuXYpOnTqhW7du53zeAQMG4N///jeeeuopXHHFFSgoKMDq1atRVFQUMFQVExOD3bt3Y8uWLUhPT0efPn2QmpqK5cuXIzo6GrIs44UXXmjUkN+UKVPw2muv4bbbbsPvfvc7pKam4rvvvsOLL76I6dOnw2QyhVwvImo8hiEiCllaWhpef/11LF26FPPmzYMQAr1798aKFStw5ZVXanKO//f//h/ee+89PP/887jkkkuwbNkyXHPNNf79ixYtwvPPP48333wTeXl5SExMxMSJEzFnzpzz9rTceOONyMnJwfr16/HGG28gOTkZY8eOxS233IJHHnkER44cQc+ePfHrX/8a+/btw6xZs7Bo0SL84he/wLJly/Dkk09i7ty5aNeuHWbOnImff/4Zx44dO+9rsdlseP3117FkyRI8/fTTqKysRMeOHfHHP/4Rt99+uyb1IqLG4+04iKhV27ZtG2699Vb861//wvDhw8PdHCK6CHHOEBEREekawxARERHpGofJiIiISNfYM0RERES6xjBEREREusYwRERERLrGMERERES6xjBEREREusYwRERERLrGMERERES6xjBEREREuvb/A03Thu9NOyiGAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.set(color_codes=True)\n",
"plt.xlim(30,90)\n",
"plt.ylim(0,1)\n",
"sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
]
}
],
"metadata": {
"celltoolbar": "Hide code",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.16"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment