{ "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 Python 3 language using the pandas, statsmodels, and numpy library." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.8.13 | packaged by conda-forge | (default, Mar 25 2022, 06:04:18) \n", "[GCC 10.3.0]\n", "uname_result(system='Linux', node='felic-ThinkPad-E14-Gen-2', release='5.15.0-91-generic', version='#101~20.04.1-Ubuntu SMP Thu Nov 16 14:22:28 UTC 2023', machine='x86_64', processor='x86_64')\n", "IPython 8.6.0\n", "IPython.core.release 8.6.0\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", "_curses b'2.2'\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", "cffi 1.15.1\n", "csv 1.0\n", "ctypes 1.1.0\n", "cycler 0.10.0\n", "dateutil 2.8.2\n", "debugpy 1.6.3\n", "debugpy.public_api 1.6.3\n", "decimal 1.70\n", "decorator 5.1.1\n", "defusedxml 0.7.1\n", "distutils 3.8.13\n", "entrypoints 0.4\n", "executing 1.2.0\n", "executing.version 1.2.0\n", "http.server 0.6\n", "ipykernel 6.17.0\n", "ipykernel._version 6.17.0\n", "ipywidgets 8.0.2\n", "ipywidgets._version 8.0.2\n", "jedi 0.18.1\n", "joblib 1.2.0\n", "joblib.externals.cloudpickle 2.2.0\n", "joblib.externals.loky 3.3.0\n", "json 2.0.9\n", "jupyter_client 7.4.4\n", "jupyter_client._version 7.4.4\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.6.2\n", "matplotlib._version 3.6.2\n", "matplotlib_inline 0.1.6\n", "numpy 1.23.4\n", "numpy.core 1.23.4\n", "numpy.core._multiarray_umath 3.1\n", "numpy.lib 1.23.4\n", "numpy.linalg._umath_linalg 0.1.5\n", "numpy.version 1.23.4\n", "packaging 21.3\n", "packaging.__about__ 21.3\n", "pandas 2.0.3\n", "parso 0.8.3\n", "patsy 0.5.6\n", "patsy.version 0.5.6\n", "pexpect 4.8.0\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.32\n", "psutil 5.9.3\n", "ptyprocess 0.7.0\n", "pure_eval 0.2.2\n", "pure_eval.version 0.2.2\n", "pydevd 2.8.0\n", "pygments 2.13.0\n", "pyparsing 3.0.9\n", "pytz 2023.3.post1\n", "re 2.2.1\n", "scipy 1.9.3\n", "scipy._lib._uarray 0.8.8.dev0+aa94c5a4.scipy\n", "scipy._lib.decorator 4.0.5\n", "scipy.integrate._dop 1.20.3\n", "scipy.integrate._lsoda 1.20.3\n", "scipy.integrate._vode 1.20.3\n", "scipy.interpolate.dfitpack 1.20.3\n", "scipy.linalg._fblas 1.20.3\n", "scipy.linalg._flapack 1.20.3\n", "scipy.linalg._flinalg 1.20.3\n", "scipy.linalg._interpolative 1.20.3\n", "scipy.optimize.__nnls 1.20.3\n", "scipy.optimize._cobyla 1.20.3\n", "scipy.optimize._lbfgsb 1.20.3\n", "scipy.optimize._minpack2 1.20.3\n", "scipy.optimize._slsqp 1.20.3\n", "scipy.sparse.linalg._eigen.arpack._arpack 1.20.3\n", "scipy.sparse.linalg._isolve._iterative 1.20.3\n", "scipy.special._specfun 1.20.3\n", "scipy.stats._mvn 1.20.3\n", "scipy.stats._statlib 1.20.3\n", "seaborn 0.12.1\n", "seaborn.external.appdirs 1.4.4\n", "seaborn.external.husl 2.1.0\n", "setuptools 65.5.0\n", "distutils 3.8.13\n", "setuptools._vendor.more_itertools 8.8.0\n", "setuptools._vendor.ordered_set 3.1\n", "setuptools._vendor.packaging 21.3\n", "setuptools._vendor.packaging.__about__ 21.3\n", "setuptools._vendor.pyparsing 3.0.9\n", "setuptools._vendor.more_itertools 8.8.0\n", "setuptools._vendor.ordered_set 3.1\n", "setuptools._vendor.packaging 21.3\n", "setuptools._vendor.pyparsing 3.0.9\n", "setuptools.version 65.5.0\n", "six 1.16.0\n", "socketserver 0.4\n", "stack_data 0.6.0\n", "stack_data.version 0.6.0\n", "statsmodels 0.14.1\n", "statsmodels.__init__ 0.14.1\n", "statsmodels._version 0.14.1\n", "statsmodels.api 0.14.1\n", "statsmodels.tools.web 0.14.1\n", "traitlets 5.5.0\n", "traitlets._version 5.5.0\n", "urllib.request 3.8\n", "wcwidth 0.2.5\n", "xmlrpc.client 3.8\n", "zlib 1.0\n", "zmq 24.0.1\n", "zmq.sugar 24.0.1\n", "zmq.sugar.version 24.0.1\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": 10, "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": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = pd.read_csv(\"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": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxdElEQVR4nO3de1xVdb7/8fcGuUjeRUBNhdTU0rwmP1K7nFSsxjRPjanlJfM8LB0dyTKalMxTlE2M1VieLl6avOUcj3XGspBkMmXGvKA1Y6h4oQwQbyGQsIP1+6Nxn7agwl5b9ubr6/l47Eeu71rfvT77Ay7frbX23g7LsiwBAAAYIsDXBQAAAHgT4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGMWn4ebzzz/X0KFD1apVKzkcDq1bt+6Sc9LT09WrVy+FhISoQ4cOWrp06WWvEwAA1B0+DTfFxcXq3r27Fi5cWK3tDx06pLvuuku33XabMjMz9dvf/lYPP/ywPvnkk8tcKQAAqCsc/vLFmQ6HQ//zP/+j4cOHX3CbWbNmaf369fr6669dY/fff79Onz6tDRs21EKVAADA39XzdQE1kZGRoYEDB7qNxcfH67e//e0F55SWlqq0tNS1XFFRoZMnT6p58+ZyOByXq1QAAOBFlmXpzJkzatWqlQICLn7hqU6Fm7y8PEVGRrqNRUZGqrCwUD/++KPq169faU5ycrLmzp1bWyUCAIDL6Ntvv9XVV1990W3qVLjxRGJiohISElzLP/zwg9q2batDhw6pYcOGPqzM+5xOpzZt2qTbbrtNQUFBvi6nTqKH9tA/++ihPfTPPn/t4ZkzZxQTE1Otf7vrVLiJiopSfn6+21h+fr4aNWpU5VkbSQoJCVFISEil8WbNmqlRo0aXpU5fcTqdCgsLU/Pmzf3qF7IuoYf20D/76KE99M8+f+3huVqqc0tJnfqcm7i4OKWlpbmNpaamKi4uzkcVAQAAf+PTcFNUVKTMzExlZmZK+vmt3pmZmcrJyZH08yWlsWPHurafPHmyDh48qCeeeELffPONXn/9db3//vuaMWOGL8oHAAB+yKfhZvv27erZs6d69uwpSUpISFDPnj01Z84cSVJubq4r6EhSTEyM1q9fr9TUVHXv3l0vv/yy3n77bcXHx/ukfgAA4H98es/Nrbfeqot9zE5Vnz586623ateuXZexKgAAUJfVqXtuAAAALoVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUn4ebhQsXKjo6WqGhoYqNjdW2bdsuuv2CBQvUqVMn1a9fX23atNGMGTN09uzZWqoWAAD4O5+Gm9WrVyshIUFJSUnauXOnunfvrvj4eB07dqzK7VesWKEnn3xSSUlJ2rt3r9555x2tXr1aTz31VC1XDgAA/JVPw01KSoomTZqkCRMm6LrrrtOiRYsUFhamxYsXV7n91q1b1a9fP40ePVrR0dEaPHiwRo0adcmzPQAA4MpRz1c7Lisr044dO5SYmOgaCwgI0MCBA5WRkVHlnJtuuknvvfeetm3bpr59++rgwYP66KOP9OCDD15wP6WlpSotLXUtFxYWSpKcTqecTqeXXo1/OPd6THtdtYke2kP/7KOH9tA/+/y1hzWpx2fh5vjx4yovL1dkZKTbeGRkpL755psq54wePVrHjx9X//79ZVmWfvrpJ02ePPmil6WSk5M1d+7cSuOffvqpwsLC7L0IP5WamurrEuo8emgP/bOPHtpD/+zztx6WlJRUe1ufhRtPpKen6/nnn9frr7+u2NhYHThwQNOnT9e8efM0e/bsKuckJiYqISHBtVxYWKg2bdpo8ODBatSoUW2VXiucTqdSU1M1aNAgBQUF+bqcOoke2kP/7KOH9tA/+/y1h+euvFSHz8JNeHi4AgMDlZ+f7zaen5+vqKioKufMnj1bDz74oB5++GFJUrdu3VRcXKz/+I//0O9+9zsFBFS+hSgkJEQhISGVxoOCgvzqh+ZNJr+22kIP7aF/9tFDe+ifff7Ww5rU4rMbioODg9W7d2+lpaW5xioqKpSWlqa4uLgq55SUlFQKMIGBgZIky7IuX7EAAKDO8OllqYSEBI0bN059+vRR3759tWDBAhUXF2vChAmSpLFjx6p169ZKTk6WJA0dOlQpKSnq2bOn67LU7NmzNXToUFfIAQAAVzafhpuRI0eqoKBAc+bMUV5ennr06KENGza4bjLOyclxO1Pz9NNPy+Fw6Omnn9bRo0fVokULDR06VM8995yvXgIAAPAzPr+heOrUqZo6dWqV69LT092W69Wrp6SkJCUlJdVCZQAAoC7y+dcvAAAAeBPhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAoPg83CxcuVHR0tEJDQxUbG6tt27ZddPvTp09rypQpatmypUJCQnTttdfqo48+qqVqAQCAv6vny52vXr1aCQkJWrRokWJjY7VgwQLFx8crKytLERERlbYvKyvToEGDFBERoT//+c9q3bq1jhw5oiZNmtR+8QAAwC/5NNykpKRo0qRJmjBhgiRp0aJFWr9+vRYvXqwnn3yy0vaLFy/WyZMntXXrVgUFBUmSoqOja7NkAADg53wWbsrKyrRjxw4lJia6xgICAjRw4EBlZGRUOefDDz9UXFycpkyZog8++EAtWrTQ6NGjNWvWLAUGBlY5p7S0VKWlpa7lwsJCSZLT6ZTT6fTiK/K9c6/HtNdVm+ihPfTPPnpoD/2zz197WJN6fBZujh8/rvLyckVGRrqNR0ZG6ptvvqlyzsGDB/XZZ59pzJgx+uijj3TgwAE9+uijcjqdSkpKqnJOcnKy5s6dW2n8008/VVhYmP0X4odSU1N9XUKdRw/toX/20UN76J99/tbDkpKSam/r08tSNVVRUaGIiAi9+eabCgwMVO/evXX06FG99NJLFww3iYmJSkhIcC0XFhaqTZs2Gjx4sBo1alRbpdcKp9Op1NRUDRo0yHXZDjVDD+2hf/bRQ3von33+2sNzV16qw2fhJjw8XIGBgcrPz3cbz8/PV1RUVJVzWrZsqaCgILdLUF26dFFeXp7KysoUHBxcaU5ISIhCQkIqjQcFBfnVD82bTH5ttYUe2kP/7KOH9tA/+/ythzWpxWdvBQ8ODlbv3r2VlpbmGquoqFBaWpri4uKqnNOvXz8dOHBAFRUVrrF9+/apZcuWVQYbAABw5fEo3Bw8eNArO09ISNBbb72lZcuWae/evXrkkUdUXFzsevfU2LFj3W44fuSRR3Ty5ElNnz5d+/bt0/r16/X8889rypQpXqkHAADUfR5dlurQoYNuueUWTZw4Uffee69CQ0M92vnIkSNVUFCgOXPmKC8vTz169NCGDRtcNxnn5OQoIOD/8lebNm30ySefaMaMGbrhhhvUunVrTZ8+XbNmzfJo/wAAwDwehZudO3dqyZIlSkhI0NSpUzVy5EhNnDhRffv2rfFzTZ06VVOnTq1yXXp6eqWxuLg4/e1vf6vxfgAAwJXBo8tSPXr00CuvvKLvv/9eixcvVm5urvr376+uXbsqJSVFBQUF3q4TAACgWmzdUFyvXj2NGDFCa9as0YsvvqgDBw5o5syZatOmjcaOHavc3Fxv1QkAAFAttsLN9u3b9eijj6ply5ZKSUnRzJkzlZ2drdTUVH3//fcaNmyYt+oEAACoFo/uuUlJSdGSJUuUlZWlO++8U++++67uvPNO182/MTExWrp0Kd/7BAAAap1H4eaNN97QQw89pPHjx6tly5ZVbhMREaF33nnHVnEAAAA15VG42b9//yW3CQ4O1rhx4zx5egAAAI95dM/NkiVLtGbNmkrja9as0bJly2wXBQAA4CmPwk1ycrLCw8MrjUdEROj555+3XRQAAICnPAo3OTk5iomJqTTerl075eTk2C4KAADAUx6Fm4iICO3Zs6fS+O7du9W8eXPbRQEAAHjKo3AzatQoTZs2TZs2bVJ5ebnKy8v12Wefafr06br//vu9XSMAAEC1efRuqXnz5unw4cO6/fbbVa/ez09RUVGhsWPHcs8NAADwKY/CTXBwsFavXq158+Zp9+7dql+/vrp166Z27dp5uz4AAIAa8SjcnHPttdfq2muv9VYtAAAAtnkUbsrLy7V06VKlpaXp2LFjqqiocFv/2WefeaU4AACAmvIo3EyfPl1Lly7VXXfdpa5du8rhcHi7LgAAAI94FG5WrVql999/X3feeae36wEAALDFo7eCBwcHq0OHDt6uBQAAwDaPws1jjz2mV155RZZlebseAAAAWzy6LPXFF19o06ZN+vjjj3X99dcrKCjIbf3atWu9UhwAAEBNeRRumjRponvuucfbtQAAANjmUbhZsmSJt+sAAADwCo/uuZGkn376SRs3btR//dd/6cyZM5Kk77//XkVFRV4rDgAAoKY8OnNz5MgRDRkyRDk5OSotLdWgQYPUsGFDvfjiiyotLdWiRYu8XScAAEC1eHTmZvr06erTp49OnTql+vXru8bvuecepaWlea04AACAmvLozM3mzZu1detWBQcHu41HR0fr6NGjXikMAADAEx6duamoqFB5eXml8e+++04NGza0XRQAAICnPAo3gwcP1oIFC1zLDodDRUVFSkpK4isZAACAT3l0Werll19WfHy8rrvuOp09e1ajR4/W/v37FR4erpUrV3q7RgAAgGrzKNxcffXV2r17t1atWqU9e/aoqKhIEydO1JgxY9xuMAYAAKhtHoUbSapXr54eeOABb9YCAABgm0fh5t13373o+rFjx3pUDAAAgF0ehZvp06e7LTudTpWUlCg4OFhhYWGEGwAA4DMevVvq1KlTbo+ioiJlZWWpf//+3FAMAAB8yuPvljpfx44d9cILL1Q6qwMAAFCbvBZupJ9vMv7++++9+ZQAAAA14tE9Nx9++KHbsmVZys3N1R//+Ef169fPK4UBAAB4wqNwM3z4cLdlh8OhFi1a6N/+7d/08ssve6MuAAAAj3gUbioqKrxdBwAAgFd49Z4bAAAAX/PozE1CQkK1t01JSfFkFwAAAB7xKNzs2rVLu3btktPpVKdOnSRJ+/btU2BgoHr16uXazuFweKdKAACAavIo3AwdOlQNGzbUsmXL1LRpU0k/f7DfhAkTNGDAAD322GNeLRIAAKC6PLrn5uWXX1ZycrIr2EhS06ZN9Z//+Z+8WwoAAPiUR+GmsLBQBQUFlcYLCgp05swZ20UBAAB4yqNwc88992jChAlau3atvvvuO3333Xf67//+b02cOFEjRozwdo0AAADV5tE9N4sWLdLMmTM1evRoOZ3On5+oXj1NnDhRL730klcLBAAAqAmPwk1YWJhef/11vfTSS8rOzpYktW/fXldddZVXiwMAAKgpWx/il5ubq9zcXHXs2FFXXXWVLMvyVl0AAAAe8SjcnDhxQrfffruuvfZa3XnnncrNzZUkTZw4kbeBAwAAn/Io3MyYMUNBQUHKyclRWFiYa3zkyJHasGGD14oDAACoKY/uufn000/1ySef6Oqrr3Yb79ixo44cOeKVwgAAADzh0Zmb4uJitzM255w8eVIhISG2iwIAAPCUR+FmwIABevfdd13LDodDFRUVmj9/vm677TavFQcAAFBTHl2Wmj9/vm6//XZt375dZWVleuKJJ/SPf/xDJ0+e1JYtW7xdIwAAQLV5dOama9eu2rdvn/r3769hw4apuLhYI0aM0K5du9S+fXtv1wgAAFBtNT5z43Q6NWTIEC1atEi/+93vLkdNAAAAHqvxmZugoCDt2bPnctQCAABgm0eXpR544AG988473q4FAADANo9uKP7pp5+0ePFibdy4Ub179670nVIpKSleKQ4AAKCmahRuDh48qOjoaH399dfq1auXJGnfvn1u2zgcDu9VBwAAUEM1CjcdO3ZUbm6uNm3aJOnnr1t49dVXFRkZeVmKAwAAqKka3XNz/rd+f/zxxyouLvZqQQAAAHZ4dEPxOeeHHU8tXLhQ0dHRCg0NVWxsrLZt21ateatWrZLD4dDw4cO9UgcAAKj7ahRuHA5HpXtq7N5js3r1aiUkJCgpKUk7d+5U9+7dFR8fr2PHjl103uHDhzVz5kwNGDDA1v4BAIBZanTPjWVZGj9+vOvLMc+ePavJkydXerfU2rVrq/2cKSkpmjRpkiZMmCBJWrRokdavX6/FixfrySefrHJOeXm5xowZo7lz52rz5s06ffp0TV4GAAAwWI3Czbhx49yWH3jgAVs7Lysr044dO5SYmOgaCwgI0MCBA5WRkXHBec8++6wiIiI0ceJEbd68+aL7KC0tVWlpqWu5sLBQ0s+ftOx0Om3V72/OvR7TXldtoof20D/76KE99M8+f+1hTeqpUbhZsmRJjYu5mOPHj6u8vLzSu60iIyP1zTffVDnniy++0DvvvKPMzMxq7SM5OVlz586tNP7pp58qLCysxjXXBampqb4uoc6jh/bQP/vooT30zz5/62FJSUm1t/XoQ/x85cyZM3rwwQf11ltvKTw8vFpzEhMTlZCQ4FouLCxUmzZtNHjwYDVq1OhyleoTTqdTqampGjRokIKCgnxdTp1ED+2hf/bRQ3von33+2sNzV16qw6fhJjw8XIGBgcrPz3cbz8/PV1RUVKXts7OzdfjwYQ0dOtQ1VlFRIUmqV6+esrKyKn0reUhIiOseoV8KCgryqx+aN5n82moLPbSH/tlHD+2hf/b5Ww9rUoutt4LbFRwcrN69eystLc01VlFRobS0NMXFxVXavnPnzvrqq6+UmZnpetx999267bbblJmZqTZt2tRm+QAAwA/5/LJUQkKCxo0bpz59+qhv375asGCBiouLXe+eGjt2rFq3bq3k5GSFhoaqa9eubvObNGkiSZXGAQDAlcnn4WbkyJEqKCjQnDlzlJeXpx49emjDhg2um4xzcnIUEODTE0wAAKAO8Xm4kaSpU6dq6tSpVa5LT0+/6NylS5d6vyAAAFBncUoEAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGy86WFCkTVnHdOh4sa9LAYAr0uF/HX+PnCjxcSXwpXq+LsAEp0vKNG1lpj7fX+Aau7ljC702qqcahwX5sDIAuDKcOw7//eAxze8r3fXaZsVeE8Fx+ArFmRsvmLYyU1sOHHcb23LguH6zcpePKgKAKwvHYfwS4camgwVF+nx/gcoty2283LL0+f4CLlEBwGXGcRjnI9zYdOTkxa/rHj7BXyoAuJw4DuN8hBub2jULu+j66OZX1VIlAHBl4jiM8xFubLqmRQPd3LGFAh0Ot/FAh0M3d2yhmHD+UgHA5cRxGOcj3HjBa6N6ql+HcLexfh3C9dqonj6qCACuLByH8Uu8FdwLGocF6d2JfXXoeLEOnyhWdPOr+D8FAKhF547DB/J+0D/+nq71vxmgDlGNfV0WfIRw40Ux4YQaAPClds3D9I9//RdXLi5LAQAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIziF+Fm4cKFio6OVmhoqGJjY7Vt27YLbvvWW29pwIABatq0qZo2baqBAwdedHsAAHBl8Xm4Wb16tRISEpSUlKSdO3eqe/fuio+P17Fjx6rcPj09XaNGjdKmTZuUkZGhNm3aaPDgwTp69GgtVw4AAPyRz8NNSkqKJk2apAkTJui6667TokWLFBYWpsWLF1e5/fLly/Xoo4+qR48e6ty5s95++21VVFQoLS2tlisHAAD+qJ4vd15WVqYdO3YoMTHRNRYQEKCBAwcqIyOjWs9RUlIip9OpZs2aVbm+tLRUpaWlruXCwkJJktPplNPptFG9/zn3ekx7XbWJHtpD/+yjh/bQP/v8tYc1qcen4eb48eMqLy9XZGSk23hkZKS++eabaj3HrFmz1KpVKw0cOLDK9cnJyZo7d26l8U8//VRhYWE1L7oOSE1N9XUJdR49tIf+2UcP7aF/9vlbD0tKSqq9rU/DjV0vvPCCVq1apfT0dIWGhla5TWJiohISElzLhYWFrvt0GjVqVFul1gqn06nU1FQNGjRIQUFBvi6nTqKH9tA/++ihPfTPPn/t4bkrL9Xh03ATHh6uwMBA5efnu43n5+crKirqonN///vf64UXXtDGjRt1ww03XHC7kJAQhYSEVBoPCgryqx+aN5n82moLPbSH/tlHD+2hf/b5Ww9rUotPbygODg5W79693W4GPndzcFxc3AXnzZ8/X/PmzdOGDRvUp0+f2igVAADUET6/LJWQkKBx48apT58+6tu3rxYsWKDi4mJNmDBBkjR27Fi1bt1aycnJkqQXX3xRc+bM0YoVKxQdHa28vDxJUoMGDdSgQQOfvQ4AAOAffB5uRo4cqYKCAs2ZM0d5eXnq0aOHNmzY4LrJOCcnRwEB/3eC6Y033lBZWZnuvfdet+dJSkrSM888U5ulAwAAP+TzcCNJU6dO1dSpU6tcl56e7rZ8+PDhy18QAACos3z+IX4AAADeRLiB1xwsKNKmrGM6dLzYJ/NBD02wZX+BJCkj+7iPKwHqLr+4LIW67XRJmaatzNTn/zooS9LNHVvotVE91Tjs0m/dszsf9NAER04Ua/jCLSopLdP8vtKkP+1QWEiwPpzSX22am/mBo8Dlwpkb2DZtZaa2HHD/v8wtB47rNyt31cp80EMTDF+4RadK3D9e/lSJU3cv/MJHFQF1F+EGthwsKNLn+wtUbllu4+WWpc/3F1zy8ojd+aCHJvhr1rFKweacUyVObf7FGTkAl0a4gS1HTl78uz4On7j4P6x254MemiDzu9MXXb8z51TtFAIYgnADW9o1u/i9ANHNr7qs80EPTdDj6iYXXd+rbdPaKQQwBOEGtlzTooFu7thCgQ6H23igw6GbO7ZQTPjF/2G1Ox/00AS3dIpQ0wvc+N00LEgDOrao5YqAuo1wA9teG9VT/TqEu4316xCu10b1rJX5oIcm+HBK/0oBp2lYkD6c0t9HFQF1F28Fh22Nw4L07sS+OnS8WIdPFCu6+VU1Oltgdz7ooQnaNA/TrjmD9fk3uTqVtU1vPdhbN3du6euygDqJcAOviQm39w+q3fmghyaIax+uj7J+/i8Az3BZCgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAo/hFuFm4cKGio6MVGhqq2NhYbdu27aLbr1mzRp07d1ZoaKi6deumjz76qJYqBXA5HSwo0qasYzp0vNij+au35ei3q3dpzfZvfbJ/u/Mlacv+AklSRvZxj5/DDn/ogR2H/7XfIydKPJr/16xjeiVtnzb/6+dQ27zRP1//DHy9f0mq57M9/8vq1auVkJCgRYsWKTY2VgsWLFB8fLyysrIUERFRafutW7dq1KhRSk5O1q9+9SutWLFCw4cP186dO9W1a1cfvAIAdp0uKdO0lZn6/Bf/oNzcsYVeG9VTjcOCLjn/q+9O657Xt+qnCkuStG7X90pc+5U+nNJP17VufNn3b3e+JB05UazhC7eopLRM8/tKk/60Q2EhwfpwSn+1aR5Wreewwx96YMe5/f/94DHN7yvd9dpmxV4TUe39n+v/qRKna6xpWFCd6b+3nsMOX+//l3x+5iYlJUWTJk3ShAkTdN1112nRokUKCwvT4sWLq9z+lVde0ZAhQ/T444+rS5cumjdvnnr16qU//vGPtVw5AG+ZtjJTWw64n6nYcuC4frNyV7Xm/zLYnPNThaW7F26plf3bnS+p0j+sknSqxKm7F35R7eewwx96YIfd/df1/nvrOezw9f5/yadnbsrKyrRjxw4lJia6xgICAjRw4EBlZGRUOScjI0MJCQluY/Hx8Vq3bl2V25eWlqq0tNS1/MMPP0iSTp48KafTWeWcusrpdKqkpEQnTpxQUFDtpmRT0EN7POlfzolibd17WA5VPiBt3VukXftbqm2zC/+f84eZR6WyogsezJZ99pV+1b3VZdu/3fmS9Lfs4zpTeFr1JNWrsFRSUqF6zgCVVzh0plD6ePs+9Y1pftHnsMMfemCH2/4D3PtXnf3/sv/nqwv999ZznOOLv8fVcebMGUmSZVmX2PLnjXzm6NGjliRr69atbuOPP/641bdv3yrnBAUFWStWrHAbW7hwoRUREVHl9klJSZYkHjx48ODBg4cBj2+//faS+cLn99xcbomJiW5neioqKnTy5Ek1b95cDofDh5V5X2Fhodq0aaNvv/1WjRo18nU5dRI9tIf+2UcP7aF/9vlrDy3L0pkzZ9Sq1YXPxJ7j03ATHh6uwMBA5efnu43n5+crKiqqyjlRUVE12j4kJEQhISFuY02aNPG86DqgUaNGfvULWRfRQ3von3300B76Z58/9rBx48bV2s6nNxQHBwerd+/eSktLc41VVFQoLS1NcXFxVc6Ji4tz216SUlNTL7g9AAC4svj8slRCQoLGjRunPn36qG/fvlqwYIGKi4s1YcIESdLYsWPVunVrJScnS5KmT5+uW265RS+//LLuuusurVq1Stu3b9ebb77py5cBAAD8hM/DzciRI1VQUKA5c+YoLy9PPXr00IYNGxQZGSlJysnJUUDA/51guummm7RixQo9/fTTeuqpp9SxY0etW7eOz7jRz5fgkpKSKl2GQ/XRQ3von3300B76Z58JPXRYVnXeUwUAAFA3+PxD/AAAALyJcAMAAIxCuAEAAEYh3AAAAKMQbuqYZ555Rg6Hw+3RuXNn1/pbb7210vrJkyf7sGL/dPToUT3wwANq3ry56tevr27dumn79u2u9ZZlac6cOWrZsqXq16+vgQMHav/+/T6s2L9cqn/jx4+v9Hs4ZMgQH1bsX6Kjoyv1x+FwaMqUKZKks2fPasqUKWrevLkaNGigf//3f6/04aVXskv1j+PgpZWXl2v27NmKiYlR/fr11b59e82bN8/te5vq8nHQ528FR81df/312rhxo2u5Xj33H+OkSZP07LPPupbDwi7fF9bVRadOnVK/fv1022236eOPP1aLFi20f/9+NW3a1LXN/Pnz9eqrr2rZsmWKiYnR7NmzFR8fr3/+858KDQ31YfW+V53+SdKQIUO0ZMkS13Jdflupt3355ZcqLy93LX/99dcaNGiQ7rvvPknSjBkztH79eq1Zs0aNGzfW1KlTNWLECG3ZUr1vOTfdpfoncRy8lBdffFFvvPGGli1bpuuvv17bt2/XhAkT1LhxY02bNk1SHT8OXvLbp+BXkpKSrO7du19w/S233GJNnz691uqpi2bNmmX179//gusrKiqsqKgo66WXXnKNnT592goJCbFWrlxZGyX6tUv1z7Isa9y4cdawYcNqpyADTJ8+3Wrfvr1VUVFhnT592goKCrLWrFnjWr93715LkpWRkeHDKv3XL/tnWRwHq+Ouu+6yHnroIbexESNGWGPGjLEsq+4fB7ksVQft379frVq10jXXXKMxY8YoJyfHbf3y5csVHh6url27KjExUSUlJT6q1D99+OGH6tOnj+677z5FRESoZ8+eeuutt1zrDx06pLy8PA0cONA11rhxY8XGxiojI8MXJfuVS/XvnPT0dEVERKhTp0565JFHdOLECR9U6//Kysr03nvv6aGHHpLD4dCOHTvkdDrdfv86d+6stm3b8vtXhfP7dw7HwYu76aablJaWpn379kmSdu/erS+++EJ33HGHpLp/HOSyVB0TGxurpUuXqlOnTsrNzdXcuXM1YMAAff3112rYsKFGjx6tdu3aqVWrVtqzZ49mzZqlrKwsrV271tel+42DBw/qjTfeUEJCgp566il9+eWXmjZtmoKDgzVu3Djl5eVJkutTss+JjIx0rbuSXap/0s+XpEaMGKGYmBhlZ2frqaee0h133KGMjAwFBgb6+BX4l3Xr1un06dMaP368JCkvL0/BwcGVvuCX37+qnd8/SRwHq+HJJ59UYWGhOnfurMDAQJWXl+u5557TmDFjJKnuHwd9feoI9pw6dcpq1KiR9fbbb1e5Pi0tzZJkHThwoJYr819BQUFWXFyc29hvfvMb6//9v/9nWZZlbdmyxZJkff/9927b3Hfffdavf/3rWqvTX12qf1XJzs62JFkbN2683OXVOYMHD7Z+9atfuZaXL19uBQcHV9ruxhtvtJ544onaLK1OOL9/VeE4WNnKlSutq6++2lq5cqW1Z88e691337WaNWtmLV261LKsun8c5LJUHdekSRNde+21OnDgQJXrY2NjJemC669ELVu21HXXXec21qVLF9flvaioKEmq9O6U/Px817or2aX6V5VrrrlG4eHh/B6e58iRI9q4caMefvhh11hUVJTKysp0+vRpt235/ausqv5VheNgZY8//riefPJJ3X///erWrZsefPBBzZgxw/Ul1XX9OEi4qeOKioqUnZ2tli1bVrk+MzNTki64/krUr18/ZWVluY3t27dP7dq1kyTFxMQoKipKaWlprvWFhYX6+9//rri4uFqt1R9dqn9V+e6773TixAl+D8+zZMkSRURE6K677nKN9e7dW0FBQW6/f1lZWcrJyeH37zxV9a8qHAcrKykpcftSakkKDAxURUWFJAOOg74+dYSaeeyxx6z09HTr0KFD1pYtW6yBAwda4eHh1rFjx6wDBw5Yzz77rLV9+3br0KFD1gcffGBdc8011s033+zrsv3Ktm3brHr16lnPPfectX//fmv58uVWWFiY9d5777m2eeGFF6wmTZpYH3zwgbVnzx5r2LBhVkxMjPXjjz/6sHL/cKn+nTlzxpo5c6aVkZFhHTp0yNq4caPVq1cvq2PHjtbZs2d9XL3/KC8vt9q2bWvNmjWr0rrJkydbbdu2tT777DNr+/btVlxcXKVLgVe6C/WP42D1jBs3zmrdurX1l7/8xTp06JC1du1aKzw83O3SZ10+DhJu6piRI0daLVu2tIKDg63WrVtbI0eOdF1HzsnJsW6++WarWbNmVkhIiNWhQwfr8ccft3744QcfV+1//vd//9fq2rWrFRISYnXu3Nl688033dZXVFRYs2fPtiIjI62QkBDr9ttvt7KysnxUrf+5WP9KSkqswYMHWy1atLCCgoKsdu3aWZMmTbLy8vJ8WLH/+eSTTyxJVf5e/fjjj9ajjz5qNW3a1AoLC7PuueceKzc31wdV+q8L9Y/jYPUUFhZa06dPt9q2bWuFhoZa11xzjfW73/3OKi0tdW1Tl4+DDsv6xccRAgAA1HHccwMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwA+CiHA7HRR/PPPOMr0v0uujoaC1YsMDXZQDwUD1fFwDAv+Xm5rr+vHr1as2ZM8ftizMbNGjgi7JqzLIslZeXq1692jvslZWVKTg4uNb2B+BnnLkBcFFRUVGuR+PGjeVwONzGVq1apS5duig0NFSdO3fW66+/7pp7+PBhORwOvf/++xowYIDq16+vG2+8Ufv27dOXX36pPn36qEGDBrrjjjtUUFDgmjd+/HgNHz5cc+fOVYsWLdSoUSNNnjxZZWVlrm0qKiqUnJysmJgY1a9fX927d9ef//xn1/r09HQ5HA59/PHH6t27t0JCQvTFF18oOztbw4YNU2RkpBo0aKAbb7xRGzdudM279dZbdeTIEc2YMcN1dkqSnnnmGfXo0cOtNwsWLFB0dHSlup977jm1atVKnTp1kiR9++23+vWvf60mTZqoWbNmGjZsmA4fPuyNHw+AKhBuAHhs+fLlmjNnjp577jnt3btXzz//vGbPnq1ly5a5bZeUlKSnn35aO3fuVL169TR69Gg98cQTeuWVV7R582YdOHBAc+bMcZuTlpamvXv3Kj09XStXrtTatWs1d+5c1/rk5GS9++67WrRokf7xj39oxowZeuCBB/TXv/7V7XmefPJJvfDCC9q7d69uuOEGFRUV6c4771RaWpp27dqlIUOGaOjQocrJyZEkrV27VldffbWeffZZ5ebmup25qo60tDRlZWUpNTVVf/nLX+R0OhUfH6+GDRtq8+bN2rJlixo0aKAhQ4a4hTUAXuTjL+4EUIcsWbLEaty4sWu5ffv21ooVK9y2mTdvnhUXF2dZlmUdOnTIkmS9/fbbrvUrV660JFlpaWmuseTkZKtTp06u5XHjxlnNmjWziouLXWNvvPGG1aBBA6u8vNw6e/asFRYWZm3dutVt3xMnTrRGjRplWZZlbdq0yZJkrVu37pKv6/rrr7dee+0113K7du2sP/zhD27bJCUlWd27d3cb+8Mf/mC1a9fOre7IyEi3b1b+05/+ZHXq1MmqqKhwjZWWllr169e3Pvnkk0vWBqDmuOcGgEeKi4uVnZ2tiRMnatKkSa7xn376SY0bN3bb9oYbbnD9OTIyUpLUrVs3t7Fjx465zenevbvCwsJcy3FxcSoqKtK3336roqIilZSUaNCgQW5zysrK1LNnT7exPn36uC0XFRXpmWee0fr165Wbm6uffvpJP/74o+vMjV3dunVzu89m9+7dOnDggBo2bOi23dmzZ5Wdne2VfQJwR7gB4JGioiJJ0ltvvaXY2Fi3dYGBgW7LQUFBrj+fu4fl/LGKiooa73v9+vVq3bq127qQkBC35auuuspteebMmUpNTdXvf/97dejQQfXr19e99957yUtEAQEBsizLbczpdFba7vz9FRUVqXfv3lq+fHmlbVu0aHHRfQLwDOEGgEciIyPVqlUrHTx4UGPGjPH68+/evVs//vij6tevL0n629/+pgYNGqhNmzZq1qyZQkJClJOTo1tuuaVGz7tlyxaNHz9e99xzj6Sfw8f5N/cGBwervLzcbaxFixbKy8uTZVmugJaZmXnJ/fXq1UurV69WRESEGjVqVKNaAXiGG4oBeGzu3LlKTk7Wq6++qn379umrr77SkiVLlJKSYvu5y8rKNHHiRP3zn//URx99pKSkJE2dOlUBAQFq2LChZs6cqRkzZmjZsmXKzs7Wzp079dprr1W6mfl8HTt21Nq1a5WZmandu3dr9OjRlc4aRUdH6/PPP9fRo0d1/PhxST+/i6qgoEDz589Xdna2Fi5cqI8//viSr2PMmDEKDw/XsGHDtHnzZh06dEjp6emaNm2avvvuO88bBOCCCDcAPPbwww/r7bff1pIlS9StWzfdcsstWrp0qWJiYmw/9+23366OHTvq5ptv1siRI3X33Xe7fWDgvHnzNHv2bCUnJ6tLly4aMmSI1q9ff8l9p6SkqGnTprrppps0dOhQxcfHq1evXm7bPPvsszp8+LDat2/vunTUpUsXvf7661q4cKG6d++ubdu2aebMmZd8HWFhYfr888/Vtm1bjRgxQl26dNHEiRN19uxZzuQAl4nDOv8iMgD42Pjx43X69GmtW7fO16UAqIM4cwMAAIxCuAEAAEbhshQAADAKZ24AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFH+PxzP5o/D94vEAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "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": 12, "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "Calling Family(..) with a link class is not allowed. Use an instance of a link class instead.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [12], line 7\u001b[0m\n\u001b[1;32m 3\u001b[0m data[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSuccess\u001b[39m\u001b[38;5;124m\"\u001b[39m]\u001b[38;5;241m=\u001b[39mdata\u001b[38;5;241m.\u001b[39mCount\u001b[38;5;241m-\u001b[39mdata\u001b[38;5;241m.\u001b[39mMalfunction\n\u001b[1;32m 4\u001b[0m data[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIntercept\u001b[39m\u001b[38;5;124m\"\u001b[39m]\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m\n\u001b[1;32m 6\u001b[0m logmodel\u001b[38;5;241m=\u001b[39msm\u001b[38;5;241m.\u001b[39mGLM(data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m'\u001b[39m], data[[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mIntercept\u001b[39m\u001b[38;5;124m'\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m]], \n\u001b[0;32m----> 7\u001b[0m family\u001b[38;5;241m=\u001b[39m\u001b[43msm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfamilies\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mBinomial\u001b[49m\u001b[43m(\u001b[49m\u001b[43msm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfamilies\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinks\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlogit\u001b[49m\u001b[43m)\u001b[49m)\u001b[38;5;241m.\u001b[39mfit()\n\u001b[1;32m 9\u001b[0m logmodel\u001b[38;5;241m.\u001b[39msummary()\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/statsmodels/genmod/families/family.py:932\u001b[0m, in \u001b[0;36mBinomial.__init__\u001b[0;34m(self, link, check_link)\u001b[0m\n\u001b[1;32m 929\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m 930\u001b[0m \u001b[38;5;66;03m# overwritten by initialize if needed but always used to initialize\u001b[39;00m\n\u001b[1;32m 931\u001b[0m \u001b[38;5;66;03m# variance since endog is assumed/forced to be (0,1)\u001b[39;00m\n\u001b[0;32m--> 932\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mBinomial\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 933\u001b[0m \u001b[43m \u001b[49m\u001b[43mlink\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlink\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 934\u001b[0m \u001b[43m \u001b[49m\u001b[43mvariance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mV\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mBinomial\u001b[49m\u001b[43m(\u001b[49m\u001b[43mn\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mn\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 935\u001b[0m \u001b[43m \u001b[49m\u001b[43mcheck_link\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcheck_link\u001b[49m\n\u001b[1;32m 936\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/statsmodels/genmod/families/family.py:94\u001b[0m, in \u001b[0;36mFamily.__init__\u001b[0;34m(self, link, variance, check_link)\u001b[0m\n\u001b[1;32m 89\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m inspect\u001b[38;5;241m.\u001b[39misclass(link):\n\u001b[1;32m 90\u001b[0m warnmssg \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 91\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCalling Family(..) with a link class is not allowed. Use an \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 92\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minstance of a link class instead.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 93\u001b[0m )\n\u001b[0;32m---> 94\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(warnmssg)\n\u001b[1;32m 96\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlink \u001b[38;5;241m=\u001b[39m link\n\u001b[1;32m 97\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvariance \u001b[38;5;241m=\u001b[39m variance\n", "\u001b[0;31mTypeError\u001b[0m: Calling Family(..) with a link class is not allowed. Use an instance of a link class instead." ] } ], "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.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": 13, "metadata": {}, "outputs": [ { "ename": "TypeError", "evalue": "Calling Family(..) with a link class is not allowed. Use an instance of a link class instead.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [13], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m logmodel\u001b[38;5;241m=\u001b[39msm\u001b[38;5;241m.\u001b[39mGLM(data[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m'\u001b[39m], data[[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mIntercept\u001b[39m\u001b[38;5;124m'\u001b[39m,\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m]], \n\u001b[0;32m----> 2\u001b[0m family\u001b[38;5;241m=\u001b[39m\u001b[43msm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfamilies\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mBinomial\u001b[49m\u001b[43m(\u001b[49m\u001b[43msm\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfamilies\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlinks\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlogit\u001b[49m\u001b[43m)\u001b[49m,\n\u001b[1;32m 3\u001b[0m var_weights\u001b[38;5;241m=\u001b[39mdata[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mCount\u001b[39m\u001b[38;5;124m'\u001b[39m])\u001b[38;5;241m.\u001b[39mfit()\n\u001b[1;32m 5\u001b[0m logmodel\u001b[38;5;241m.\u001b[39msummary()\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/statsmodels/genmod/families/family.py:932\u001b[0m, in \u001b[0;36mBinomial.__init__\u001b[0;34m(self, link, check_link)\u001b[0m\n\u001b[1;32m 929\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mn \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m 930\u001b[0m \u001b[38;5;66;03m# overwritten by initialize if needed but always used to initialize\u001b[39;00m\n\u001b[1;32m 931\u001b[0m \u001b[38;5;66;03m# variance since endog is assumed/forced to be (0,1)\u001b[39;00m\n\u001b[0;32m--> 932\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mBinomial\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\n\u001b[1;32m 933\u001b[0m \u001b[43m \u001b[49m\u001b[43mlink\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mlink\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 934\u001b[0m \u001b[43m \u001b[49m\u001b[43mvariance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mV\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mBinomial\u001b[49m\u001b[43m(\u001b[49m\u001b[43mn\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mn\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 935\u001b[0m \u001b[43m \u001b[49m\u001b[43mcheck_link\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcheck_link\u001b[49m\n\u001b[1;32m 936\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/statsmodels/genmod/families/family.py:94\u001b[0m, in \u001b[0;36mFamily.__init__\u001b[0;34m(self, link, variance, check_link)\u001b[0m\n\u001b[1;32m 89\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m inspect\u001b[38;5;241m.\u001b[39misclass(link):\n\u001b[1;32m 90\u001b[0m warnmssg \u001b[38;5;241m=\u001b[39m (\n\u001b[1;32m 91\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mCalling Family(..) with a link class is not allowed. Use an \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 92\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minstance of a link class instead.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 93\u001b[0m )\n\u001b[0;32m---> 94\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(warnmssg)\n\u001b[1;32m 96\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mlink \u001b[38;5;241m=\u001b[39m link\n\u001b[1;32m 97\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mvariance \u001b[38;5;241m=\u001b[39m variance\n", "\u001b[0;31mTypeError\u001b[0m: Calling Family(..) with a link class is not allowed. Use an instance of a link class instead." ] } ], "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": 14, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'logmodel' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [14], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m data_pred \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mDataFrame({\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m: np\u001b[38;5;241m.\u001b[39mlinspace(start\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m30\u001b[39m, stop\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m90\u001b[39m, num\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m121\u001b[39m),\n\u001b[1;32m 2\u001b[0m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mIntercept\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;241m1\u001b[39m})\n\u001b[0;32m----> 3\u001b[0m data_pred[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[43mlogmodel\u001b[49m\u001b[38;5;241m.\u001b[39mpredict(data_pred)\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28mprint\u001b[39m(data_pred\u001b[38;5;241m.\u001b[39mhead())\n", "\u001b[0;31mNameError\u001b[0m: name 'logmodel' is not defined" ] } ], "source": [ "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121),\n", " 'Intercept': 1})\n", "data_pred['Frequency'] = logmodel.predict(data_pred)\n", "print(data_pred.head())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "ename": "KeyError", "evalue": "'Frequency'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/indexes/base.py:3653\u001b[0m, in \u001b[0;36mIndex.get_loc\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3652\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m-> 3653\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_engine\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_loc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcasted_key\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3654\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m err:\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/_libs/index.pyx:147\u001b[0m, in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/_libs/index.pyx:176\u001b[0m, in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32mpandas/_libs/hashtable_class_helper.pxi:7080\u001b[0m, in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32mpandas/_libs/hashtable_class_helper.pxi:7088\u001b[0m, in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n", "\u001b[0;31mKeyError\u001b[0m: 'Frequency'", "\nThe above exception was the direct cause of the following exception:\n", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [8], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m get_ipython()\u001b[38;5;241m.\u001b[39mrun_line_magic(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmatplotlib\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124minline\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 2\u001b[0m \u001b[43mdata_pred\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mTemperature\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43my\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mFrequency\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mkind\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mline\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mylim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3\u001b[0m plt\u001b[38;5;241m.\u001b[39mscatter(x\u001b[38;5;241m=\u001b[39mdata[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m\"\u001b[39m],y\u001b[38;5;241m=\u001b[39mdata[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m\"\u001b[39m])\n\u001b[1;32m 4\u001b[0m plt\u001b[38;5;241m.\u001b[39mgrid(\u001b[38;5;28;01mTrue\u001b[39;00m)\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/plotting/_core.py:961\u001b[0m, in \u001b[0;36mPlotAccessor.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 958\u001b[0m \u001b[38;5;28;01mpass\u001b[39;00m\n\u001b[1;32m 960\u001b[0m \u001b[38;5;66;03m# don't overwrite\u001b[39;00m\n\u001b[0;32m--> 961\u001b[0m data \u001b[38;5;241m=\u001b[39m \u001b[43mdata\u001b[49m\u001b[43m[\u001b[49m\u001b[43my\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mcopy()\n\u001b[1;32m 963\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(data, ABCSeries):\n\u001b[1;32m 964\u001b[0m label_name \u001b[38;5;241m=\u001b[39m label_kw \u001b[38;5;129;01mor\u001b[39;00m y\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/frame.py:3761\u001b[0m, in \u001b[0;36mDataFrame.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3759\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcolumns\u001b[38;5;241m.\u001b[39mnlevels \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[1;32m 3760\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_getitem_multilevel(key)\n\u001b[0;32m-> 3761\u001b[0m indexer \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcolumns\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_loc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3762\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m is_integer(indexer):\n\u001b[1;32m 3763\u001b[0m indexer \u001b[38;5;241m=\u001b[39m [indexer]\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/indexes/base.py:3655\u001b[0m, in \u001b[0;36mIndex.get_loc\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3653\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_engine\u001b[38;5;241m.\u001b[39mget_loc(casted_key)\n\u001b[1;32m 3654\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m err:\n\u001b[0;32m-> 3655\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(key) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01merr\u001b[39;00m\n\u001b[1;32m 3656\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m:\n\u001b[1;32m 3657\u001b[0m \u001b[38;5;66;03m# If we have a listlike key, _check_indexing_error will raise\u001b[39;00m\n\u001b[1;32m 3658\u001b[0m \u001b[38;5;66;03m# InvalidIndexError. Otherwise we fall through and re-raise\u001b[39;00m\n\u001b[1;32m 3659\u001b[0m \u001b[38;5;66;03m# the TypeError.\u001b[39;00m\n\u001b[1;32m 3660\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_check_indexing_error(key)\n", "\u001b[0;31mKeyError\u001b[0m: 'Frequency'" ] } ], "source": [ "%matplotlib inline\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": {}, "source": [ "La fonction `logmodel.predict(data_pred)` ne fonctionne pas avec les dernières versions de pandas (`Frequency = 1` pour toutes les températures).\n", "\n", "On peut alors utiliser le code suivant pour calculer les prédictions et tracer la courbe :" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'logmodel' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [15], line 5\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mlogit_inv\u001b[39m(x):\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m(np\u001b[38;5;241m.\u001b[39mexp(x)\u001b[38;5;241m/\u001b[39m(np\u001b[38;5;241m.\u001b[39mexp(x)\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m))\n\u001b[0;32m----> 5\u001b[0m data_pred[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mProb\u001b[39m\u001b[38;5;124m'\u001b[39m]\u001b[38;5;241m=\u001b[39mlogit_inv(data_pred[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m*\u001b[39m \u001b[43mlogmodel\u001b[49m\u001b[38;5;241m.\u001b[39mparams[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m+\u001b[39m \n\u001b[1;32m 6\u001b[0m logmodel\u001b[38;5;241m.\u001b[39mparams[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mIntercept\u001b[39m\u001b[38;5;124m'\u001b[39m])\n\u001b[1;32m 7\u001b[0m \u001b[38;5;28mprint\u001b[39m(data_pred\u001b[38;5;241m.\u001b[39mhead())\n", "\u001b[0;31mNameError\u001b[0m: name 'logmodel' is not defined" ] } ], "source": [ "# Inspiring from http://blog.yhat.com/posts/logistic-regression-and-python.html\n", "def logit_inv(x):\n", " return(np.exp(x)/(np.exp(x)+1))\n", "\n", "data_pred['Prob']=logit_inv(data_pred['Temperature'] * logmodel.params['Temperature'] + \n", " logmodel.params['Intercept'])\n", "print(data_pred.head())" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "ename": "KeyError", "evalue": "'Prob'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/indexes/base.py:3653\u001b[0m, in \u001b[0;36mIndex.get_loc\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3652\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m-> 3653\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_engine\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_loc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcasted_key\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3654\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m err:\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/_libs/index.pyx:147\u001b[0m, in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/_libs/index.pyx:176\u001b[0m, in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32mpandas/_libs/hashtable_class_helper.pxi:7080\u001b[0m, in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n", "File \u001b[0;32mpandas/_libs/hashtable_class_helper.pxi:7088\u001b[0m, in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n", "\u001b[0;31mKeyError\u001b[0m: 'Prob'", "\nThe above exception was the direct cause of the following exception:\n", "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn [16], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m get_ipython()\u001b[38;5;241m.\u001b[39mrun_line_magic(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmatplotlib\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124minline\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m----> 2\u001b[0m \u001b[43mdata_pred\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mTemperature\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43my\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mProb\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mkind\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mline\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43mylim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3\u001b[0m plt\u001b[38;5;241m.\u001b[39mscatter(x\u001b[38;5;241m=\u001b[39mdata[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m\"\u001b[39m],y\u001b[38;5;241m=\u001b[39mdata[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m\"\u001b[39m])\n\u001b[1;32m 4\u001b[0m plt\u001b[38;5;241m.\u001b[39mgrid(\u001b[38;5;28;01mTrue\u001b[39;00m)\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/plotting/_core.py:961\u001b[0m, in \u001b[0;36mPlotAccessor.__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 958\u001b[0m \u001b[38;5;28;01mpass\u001b[39;00m\n\u001b[1;32m 960\u001b[0m \u001b[38;5;66;03m# don't overwrite\u001b[39;00m\n\u001b[0;32m--> 961\u001b[0m data \u001b[38;5;241m=\u001b[39m \u001b[43mdata\u001b[49m\u001b[43m[\u001b[49m\u001b[43my\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mcopy()\n\u001b[1;32m 963\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(data, ABCSeries):\n\u001b[1;32m 964\u001b[0m label_name \u001b[38;5;241m=\u001b[39m label_kw \u001b[38;5;129;01mor\u001b[39;00m y\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/frame.py:3761\u001b[0m, in \u001b[0;36mDataFrame.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3759\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcolumns\u001b[38;5;241m.\u001b[39mnlevels \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[1;32m 3760\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_getitem_multilevel(key)\n\u001b[0;32m-> 3761\u001b[0m indexer \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcolumns\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_loc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 3762\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m is_integer(indexer):\n\u001b[1;32m 3763\u001b[0m indexer \u001b[38;5;241m=\u001b[39m [indexer]\n", "File \u001b[0;32m~/Programs/miniconda3/envs/cosmicfish/lib/python3.8/site-packages/pandas/core/indexes/base.py:3655\u001b[0m, in \u001b[0;36mIndex.get_loc\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 3653\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_engine\u001b[38;5;241m.\u001b[39mget_loc(casted_key)\n\u001b[1;32m 3654\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m err:\n\u001b[0;32m-> 3655\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m(key) \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01merr\u001b[39;00m\n\u001b[1;32m 3656\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m:\n\u001b[1;32m 3657\u001b[0m \u001b[38;5;66;03m# If we have a listlike key, _check_indexing_error will raise\u001b[39;00m\n\u001b[1;32m 3658\u001b[0m \u001b[38;5;66;03m# InvalidIndexError. Otherwise we fall through and re-raise\u001b[39;00m\n\u001b[1;32m 3659\u001b[0m \u001b[38;5;66;03m# the TypeError.\u001b[39;00m\n\u001b[1;32m 3660\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_check_indexing_error(key)\n", "\u001b[0;31mKeyError\u001b[0m: 'Prob'" ] } ], "source": [ "%matplotlib inline\n", "data_pred.plot(x=\"Temperature\",y=\"Prob\",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": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAG/CAYAAACnlPiuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWlklEQVR4nO3deXhU1cEG8PfeOzOZTJLJAoGwQ4AECLsKhCAirgii4kbVFhUrtKlWtF9F/UrlQ1uli1VwX1q0itWCK4JiVUA2tSDIIksCYQmEkG0myaz3nu+PSYYMyYWQuclMkvf3PBHmrmcOk+T1nHPPkYQQAkRERERUjxzpAhARERFFKwYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRjqgKSgUFBZg3bx6uueYaDBo0CFOmTGnUeUIIvPTSS5gwYQKGDh2Km2++Gd9//33zFpaIiIjavKgKSvv27cOaNWvQq1cv9O3bt9Hnvfzyy3jmmWdw++2348UXX0RqairuvPNOHD58uBlLS0RERG2dFE1rvWmaBlkOZLe5c+dix44d+Pjjj894jsfjwdixY3Hrrbfi/vvvBwB4vV5ceeWVGD9+PB599NHmLjYRERG1UVHVolQbks7Fli1bUFlZiUmTJgW3WSwWXHbZZVi7dq2RxSMiIqJ2JqqCUlPk5+cDANLT00O29+3bF4WFhXC73ZEoFhEREbUBrT4oORwOWCwWxMTEhGy32+0QQqCioiJCJSMiIqLWrtUHpeYURcO3iIiIKAJMkS5AuOx2O7xeLzweT0irksPhgCRJSExMbPK1JUmCw+GCqmpGFLXdUhQZdnss6zJMrEfjsC6Nw7o0BuvROImJsU0a86yn1Qel2rFJBw4cwIABA4Lb8/Pz0bVrV1it1rCur6oa/H5+aI3AujQG69E4rEvjsC6NwXoMn9GdQa2+623kyJGIj4/HypUrg9t8Ph8+++wzjB8/PoIlIyIiotYuqlqUXC4X1qxZAwA4evQoKisrsWrVKgDAqFGjkJKSghkzZqCwsBCrV68GAMTExGDWrFlYtGgRUlJSkJGRgaVLl6K8vBwzZ86M2HshIiKi1i+qglJJSQl+/etfh2yrff36669j9OjR0DQNqqqGHPPzn/8cQgi89tprKC0txcCBA/Hqq6+iR48eLVZ2IiIianuiambuaFRWVsX+4jCZTDKSk+NYl2FiPRqHdWkc1qUxWI/GSUmJg6IYN7Ko1Y9RIiIiImouDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0hF1QSkvLw933HEHhg8fjpycHCxcuBBer/es55WVlWHevHmYMGEChg8fjilTpmDp0qUtUGIiIiJqq0yRLkBdFRUVmDFjBnr37o1FixahqKgITzzxBNxuN+bNm3fGc3/9618jPz8f999/P7p06YK1a9fi0UcfhaIouOmmm1roHRAREVFbElVB6e2330ZVVRUWL16MpKQkAICqqpg/fz5mzZqFzp07N3hecXExNm/ejD/+8Y+YNm0aACA7Oxs//PADVqxYwaBERERETRJVXW9r165FdnZ2MCQBwKRJk6BpGtavX697nt/vBwAkJCSEbI+Pj4cQolnKSkRERG1fVAWl/Px8pKenh2yz2+1ITU1Ffn6+7nldunTBuHHj8MILL2D//v2orKzEJ598gvXr1+PWW29t7mITERFRGxVVXW8OhwN2u73e9sTERFRUVJzx3EWLFmHOnDmYPHkyAEBRFPzv//4vrrjiirDKpChRlSVbpdo6ZF2Gh/VoHNalcViXxmA9GkeSjL1eVAWlphJC4KGHHsLBgwfxl7/8BampqdiwYQP+8Ic/IDExMRiemsJujzWwpO0b69IYrEfjsC6Nw7o0Busx+kRVULLb7XA6nfW2V1RUIDExUfe8r776CqtWrcKHH36IzMxMAMDo0aNRUlKCJ554Iqyg5HC4oKpak8+nwP8h2e2xrMswsR6Nw7o0DuvSGKxH4yQmxkKWjWuZi6qglJ6eXm8sktPpRHFxcb2xS3Xt378fiqIgIyMjZPvAgQPx7rvvwuVyITa2aSldVTX4/fzQGoF1aQzWo3FYl8ZhXRqD9Rg+o5/hiqrO0PHjx2PDhg1wOBzBbatWrYIsy8jJydE9r1u3blBVFXv27AnZvnPnTnTo0KHJIYmIiIjat6gKStOnT0dcXBxyc3Px9ddfY9myZVi4cCGmT58eMofSjBkzcNlllwVfjx8/Hl27dsW9996LDz74ABs3bsSf/vQnvPfee7jtttsi8VaIiIioDYiqrrfExEQsWbIECxYsQG5uLuLi4nDDDTdgzpw5IcdpmgZVVYOv4+Pj8Y9//ANPPfUU/vznP8PpdKJ79+6YO3cugxIRERE1mSQ4I+MZlZVVsb84TCaTjOTkONZlmFiPxmFdGod1aQzWo3FSUuIMnWYhqrreiIiIiKIJgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdDAoEREREelgUCIiIiLSwaBEREREpINBiYiIiEgHgxIRERGRDgYlIiIiIh0MSkREREQ6GJSIiIiIdERdUMrLy8Mdd9yB4cOHIycnBwsXLoTX623UuUVFRXjwwQcxZswYDB06FJMmTcKHH37YzCUmIiKitsoU6QLUVVFRgRkzZqB3795YtGgRioqK8MQTT8DtdmPevHlnPPfEiRO4+eab0adPHyxYsADx8fHYt29fo0MWERER0emiKii9/fbbqKqqwuLFi5GUlAQAUFUV8+fPx6xZs9C5c2fdc//0pz8hLS0Nr7zyChRFAQBkZ2e3RLGJiIiojQqr6+3EiRNGlQMAsHbtWmRnZwdDEgBMmjQJmqZh/fr1uudVVlZi5cqVuOWWW4IhiYiIiChcYbUoTZgwAWPGjMHUqVNx+eWXw2azhVWY/Px8XH/99SHb7HY7UlNTkZ+fr3vezp074fP5YDKZcNttt2Hr1q1ISkrCtddei/vuuw9ms7nJZVKUqBvG1erU1iHrMjysR+OwLo3DujQG69E4kmTs9cIKSvfeey8+/vhjzJ07F/Pnz8cll1yCqVOnYty4cZDlc//HdjgcsNvt9bYnJiaioqJC97yTJ08CAP73f/8XN910E371q19h+/bteOaZZyDLMh544IFzLkstuz22yedSKNalMViPxmFdGod1aQzWY/QJKyjNnj0bs2fPxq5du/DRRx9hxYoV+Pjjj9GhQwdMnjwZV199NYYMGWJUWXVpmgYAGDt2LObOnQsAGDNmDKqqqvDaa68hNzcXVqu1Sdd2OFxQVc2wsrZHiiLDbo9lXYaJ9Wgc1qVxWJfGYD0aJzExtkmNNXoMGcw9aNAgDBo0CL/97W+xadMmfPTRR1i+fDneeOMN9OnTB1OnTsXUqVPRtWvXM17HbrfD6XTW215RUYHExMQzngcEwlFd2dnZeOGFF1BQUIDMzMwmvDNAVTX4/fzQGoF1aQzWo3FYl8ZhXRqD9Rg+IYy9nqGdoZIk4bzzzsNFF12EYcOGQQiBgoICLF68GJdeeinuvffeMw4AT09PrzcWyel0ori4GOnp6brn9evX74zl8ng85/ZGiIiIiGBgUNq0aRMeeeQR5OTk4L777sPJkyfx4IMPYs2aNVi3bh0eeOABbNq0Cb/97W91rzF+/Hhs2LABDocjuG3VqlWQZRk5OTm653Xr1g0ZGRnYsGFDyPYNGzbAarWeNUgRERERNSSsrrcff/wRH374IVasWIETJ06gY8eOuOGGG3DttdfW6+qaOXMmYmJi8OSTT+peb/r06XjjjTeQm5uLWbNmoaioCAsXLsT06dND5lCaMWMGCgsLsXr16uC2OXPm4Je//CUef/xxTJgwAT/88ANee+01zJw5M+yn8YiIiKh9CisoXXvttbBarbjkkktw7bXXIicn54wDqPr164fhw4fr7k9MTMSSJUuwYMEC5ObmIi4uDjfccAPmzJkTcpymaVBVNWTbxIkT8de//hXPPfccli5dik6dOuGee+7B3XffHc5bJCIionZMEqLpw56WL1+OK664AnFxcUaWKaqUlVVxYF2YTCYZyclxrMswsR6Nw7o0DuvSGKxH46SkxBk6H1VYLUrTpk0zqhxEREREUSesyPX6669j5syZuvvvuusuvPXWW+HcgoiIiChiwgpK//73v9G3b1/d/f369cM777wTzi2IiIiIIiasoHT48OEzBqX09HQcOnQonFsQERERRUxYQclsNqO4uFh3/4kTJwydRpyIiIioJYWVYoYNG4b33nsPlZWV9fY5nU4sX74cw4YNC+cWRERERBET1lNvv/rVr3Dbbbfh2muvxYwZM4IzYO/btw9LlixBcXEx/vKXvxhSUCIiIqKWFlZQGjZsGF544QXMmzcPjz/+OCRJAgAIIdC9e3c8//zzGDFihCEFJSIiImppYQUlAMjJycHq1auxa9eu4MDtnj17IisrKxiciIiIiFqjsIMSAMiyjMGDB2Pw4MFGXI6IiIgoKhgSlPbv34/Dhw+joqKiwf3XXnutEbchIiIialFhBaVDhw7hf/7nf7B9+3boLRknSRKDEhEREbVKYQWlefPmYe/evXj44Ydx/vnnw263G1UuIiIioogLKyht2bIFs2bNwk9/+lOjykNEREQUNcKacDI5ORkJCQlGlYWIiIgoqoQVlKZPn44PP/wQqqoaVR4iIiKiqBFW11vv3r2haRquueYaXH/99UhLS4OiKPWOu/zyy8O5DREREVFEhBWU5syZE/z7k08+2eAxkiRh9+7d4dyGiIiIKCLCCkqvv/66UeUgIiIiijphBaVRo0YZVQ4iIiKiqGPIzNxerxc7d+5ESUkJRo4ciZSUFCMuS0RERBRRYT31BgS638aNG4dbbrkF99xzD/bs2QMAKC0txejRo/Hvf/877EISERERRUJYQWnZsmX4wx/+gAsvvBCPP/54yDImKSkpGDNmDD755JOwC0lEREQUCWEFpb///e+45JJL8Je//AUXX3xxvf1ZWVnYt29fOLcgIiIiipiwglJBQQHGjx+vuz8pKQnl5eXh3IKIiIgoYsIKSna7HWVlZbr79+/fj9TU1HBuQURERBQxYQWl8ePH45133oHD4ai3b9++fXj33XcxceLEcG5BREREFDFhTQ9w33334aabbsKUKVNw8cUXQ5IkvP/++1i2bBk+++wzpKam4pe//KVRZSUiIiJqUWG1KHXu3BnLly/HhRdeiJUrV0IIgQ8++ABffvklJk+ejHfeeYdzKhEREVGrJYm6z/SHqbS0FJqmISUlBbIc9hRNUaGsrAp+vxbpYrRqJpOM5OQ41mWYWI/GYV0ah3VpDNajcVJS4qAoxmUQQ2bmrsXWIyIiImpLwgpKixcvPusxkiQhNzc3nNsQERERRUSzBSVJkiCEYFAiIiKiViusoPTjjz/W26ZpGo4ePYq33noL3377LV5++eVwbkFEREQUMYaPuJZlGT169MCDDz6IXr164bHHHjP6FkREREQtolkfTbvggguwZs2a5rwFERERUbNp1qC0Y8eONjNNABEREbU/YY1Rev/99xvc7nA48N133+Gzzz7DjTfeGM4tiIiIiCImrKA0d+5c3X3Jycm4++67+cQbERERtVphBaX//Oc/9bZJkgS73Y74+PhwLk1EREQUcWEFpW7duhlVDiIiIqKow5HWRERERDrCalEaMGAAJEk6p3MkScKuXbvCuS0RERFRiwgrKOXm5uLzzz/H/v37MW7cOPTp0wcAkJ+fj/Xr16N///649NJLDSkoERERUUsLKyh16tQJJSUl+Oijj5Cenh6yLy8vDzNmzECnTp1w0003hVVIIiIiokgIa4zSq6++ittuu61eSAKAvn374tZbb8Urr7wSzi2IiIiIIiasoHT8+HGYTPqNUiaTCcePHw/nFkREREQRE1ZQ6t+/P9566y0UFRXV23f8+HEsXboUGRkZ4dyCiIiIKGLCGqP00EMP4a677sIVV1yBSy+9FL169QIAHDx4EP/5z38ghMDChQsNKSgRERFRSwsrKJ1//vl455138PTTT+Pzzz+H2+0GAFitVowbNw733HMPMjMzDSkoERERUUsLKygBQEZGBp599llomobS0lIAQEpKCmSZc1kSERFR6xZ2UKolyzJiYmJgs9kYkoiIiKhNCDvR/PDDD5g5cyaGDRuG0aNH45tvvgEAlJaW4he/+AU2b94cdiGJiIiIIiGsoLRlyxbccsstKCgowNSpU6FpWnBfSkoKKisr8a9//SvsQhIRERFFQlhB6amnnkLfvn3xySefYM6cOfX2jx49Gtu2bQvnFkREREQRE1ZQ+uGHHzBt2jRYLJYGF8ft3LkzTp48Gc4tiIiIiCImrKBkMplCuttOV1RUBJvNFs4tiIiIiCImrKA0bNgwfPrppw3uq66uxvLly3HBBReEcwsiIiKiiAkrKN17773YsWMH7r77bqxduxYAsGfPHrz77ruYNm0aSktL8ctf/tKQghJRy2igF52IqN2ShBAinAts3LgRjz76KAoKCkK29+zZE4899hhGjRoVVgEjraysCn6/fvcinZ3JJCM5OY51GaaWqkcBQJaA8H4yRDd+Jo3DujQG69E4KSlxUBTj5nNs8oSTQghUVVVh5MiR+PTTT7F7924cPHgQQgj06NEDgwcPbnCANxFFN7+qwaRIkMDvXyKiJkcun8+HUaNG4fXXXwcADBw4EJMmTcJVV12FIUOGNDkk5eXl4Y477sDw4cORk5ODhQsXwuv1ntM1/vGPfyAzMxOzZs1qUhmI2juPV410EYiIokKTW5QsFgs6duwIi8ViWGEqKiowY8YM9O7dG4sWLUJRURGeeOIJuN1uzJs3r1HXKC4uxrPPPosOHToYVi6i9sbt02CzStC0Ntz/RkTUCGGt9Xbdddfhgw8+wE9+8hNDAtPbb7+NqqoqLF68GElJSQAAVVUxf/58zJo1C507dz7rNf70pz9h4sSJKCwsDLs8RO2Vqmrw+TQoCrvfiKh9C2u0U2ZmJrxeL6ZMmYLnn38eH374IT777LN6X421du1aZGdnB0MSAEyaNAmapmH9+vVnPf+7777D559/jgceeKApb4eIaqiagMvr5xNwRNTuhdWidP/99wf//vTTTzd4jCRJ2L17d6Oul5+fj+uvvz5km91uR2pqKvLz8894rqqqWLBgAWbPno1OnTo16n6NYeTI+faqtg5Zl+FpqXrUACiyBL8mILfRQd38TBqHdWkM1qNxjP4fvHMOSn/9619x1VVXYcCAAcGB3EZxOByw2+31ticmJqKiouKM57711ltwuVy4/fbbDS2T3R5r6PXaM9alMZq7Hl0eH1RJghCAyWyBPd64cYjRhp9J47AujcF6jD7nHJReeukl9O/fHwMGDMCoUaNQVlaGsWPH4rXXXkN2dnZzlPGsSkpK8Mwzz+DJJ580dHA5ADgcLqgq57QIh6LIsNtjWZdhaql69Po1OBxuaEKgusoD1WdttntFCj+TxmFdGoP1aJzExFjIchTMo1RXmHNWBtntdjidznrbKyoqkJiYqHve008/jczMTJx//vlwOBwAAL/fD7/fD4fDAZvNBpOpaW9VVTVO/mUQ1qUxmrseVVWDqgpoQkDT/Kis9sJqMeRHRdThZ9I4rEtjsB7DZ/RkuVH10y89Pb3eWCSn04ni4mKkp6frnnfgwAF8++23Da4rd8EFF+Dll1/G+PHjDS8vUVsnBFDt8bfZoEREdDZR9dNv/PjxeOGFF0LGKq1atQqyLCMnJ0f3vIcffjjYklTrD3/4A6xWK+6//35kZmY2a7mJ2jKfX4Pbp8JqViJdFCKiFtekoHT06FHs3LkTAIJdZQUFBQ0OxAaArKysRl13+vTpeOONN5Cbm4tZs2ahqKgICxcuxPTp00PmUJoxYwYKCwuxevVqAIFZwU9nt9ths9kwevToc3pvRBRKCKDa7YPVogQWgiMiakeaFJSefvrpetMBzJ8/v95xQohzmh4gMTERS5YswYIFC5Cbm4u4uDjccMMNmDNnTshxmqZBVbnEAlFL8fk1uL1sVSKi9kcS5zgS+7333jvnm1x33XXnfE604ErO4eOq2MZoqXr0qRrKHB5op/1oiDHLSLZb20SrEj+TxmFdGoP1aJyUlDhD56M65xal1hx6iKjpfH4Bj09FjImtSkTUfnAKUCJqFE0IVLt8XNaEiNoVBiUiarTaViUiovaCQYmIGk0TAlVuLpZLRO0HgxIRnROfX2OrEhG1GwxKRHRONI2tSkTUfjAoEdE58/k1ePkIMxG1AwxKRHTONK32CTg2KxFR28agRERN4vFr8Po5VomI2jYGJSJqEk0TqHb72apERG0agxIRNZnXp3K5BSJq0xiUiKjJVE2g2sOxSkTUdp3zWm9E1DZpQuBQkRMujwqPz49OyTbIjQhAbp8Km6pBkRmWKPrUfq4rq32It5nRs3NCoz7XRLUYlIgIuw+WYsWmAhSVVqNDYiwqKt2It1kwYVhX9O2edMZzVVWg2u2DPS4GQoiWKTBRI9R+ro+XVkNVBRRFQlqKDZPH9MLA3imRLh61Eux6I2rndh8sxZJP9+BIcSVizAriYs2wmBUcL3Xhva8PIO9I+Vmv4fap8Kscq0TR4/TPtT3eghizgiPFVVjy6R7sPlga6SJSK8GgRNSOaUJgxaYCuL1+JMXHwGJWIEuAyaTAbjPD49Pw1bZCaGdpKVJVgWoPn4Cj6NDw51qCxawgKd4Ct1fFik0FZ/1cEwEMSkTt2qEiJ46XViPOaq4XciRJgi1GwclyF46drDrrtdweP3ycV4miwNk+13FWE46XVuNQkTNCJaTWhEGJqB2rrPZBVQVMpoZ/FCiKDFUDqt3+s15L1QQqXT6uAUcRd7bPtckkQ1UFKqt9LVwyao0YlIjasXibGYoi6c6FpKoaFBmwWRv33IfXp8HlZasSRdbZPtd+vwZFkRBvM7dwyag1YlAiasd6dk5AWooNVW5/vSfWhBCo9qjomBSLLh3jGnU9TQSegCOKpLN9rqvcfqSl2NCzc0KESkitCYMSUTsmSxImj+kFq0VBeaUXXp8KTQB+vwpHtQ8xZhkThnU9p3lnfH4Nbh9blShyGv5cC3h9KsorvbBaFEwe04vzKVGjMCgRtXMDe6dgxhWZ6J4aB49PRZXLB69PRVpKLK4b1+es8yidTojAwG6ZE1BSBJ3+uXZUeuHxqeieGocZV2RyHiVqNE44SUQY2DsFmb2SmzQzd0O8fg1enwqTwv8Xo8ip+7nmzNzUVAxKRAQg0F3RO80On6qhzOEJa44ZTRNwefywx1nAqWookmo/10RNxf/dI6Jm4fap4GTdRNTaMSgRUbPQVAGP1895lYioVWNQIqJmIQC4vGefqJKIKJoxKBFRs/GrAl6dSf+IiFoDBiUiajaaJuD2sPuNiFovBiUialYen8ZV2omo1WJQIqJmpWoa3Fz/jYhaKQYlImpWQgAuDwd1E1HrxKBERM3O59f4BBwRtUoMSkTU7IQAqt1+BCYNICJqPRiUiKhF+FUNLo/KJ+CIqFVhUCKiFlHbqsQn4IioNWFQIqIW41M1VLr8kNisREStBIMSEbUot8cPr5/TBRBR68CgREQtStUEKl0+gI1KRNQKMCgRUYvz+lS43FzahIiiH4MSEbU4IYAqtx9+lQO7iSi6MSgRUUT4VQ2V1V62KhFRVDNFugBEFB00IbByUwFKnR5kD0pDvM3c7Pd0+1RUe/ywxZjAWQOIKBoxKBERAGDbvpNYtiY/+PeZUwYhOSGmWe8pBFDt8sNqUSBxdDcRRSF2vRERACAx/lQoKnV68PJHO1HqcDf7fX2qhioXB3YTUXRiUCIiAEB6VzumX9I/+Lq80ouXP9qFkormD0surx8+DuwmoijEoEREQZdf0AM3X9Iv+LqiyouXP9qJE+WuZr2vqgpUuXycsZuIog6DEhGFmDiyO64d1yf42lHtwysf7cLx0upmva/Hp3LGbiKKOgxKRFTPmKw0XDc+PTi8utLlw8sf7cLR4spmu6em1bYqNdstiIjOGYMSETXoggGdcP2EvsHg4vL48crHu1Fw3Nls9/T6NLh9bFUioujBoEREukZmpOLmif0h16Qlj0/Fa5/sxr4j5c1yP00IVHMdOCKKIgxKRHRGQ/t2wG2XZ8CkBNKLz6/h9VV7sONAabPcz+vXUM114IgoSjAoEdFZDeiVjJ9dOQAWU+BHhqoJLP18L7778YTh9xICqHL54PFphl+biOhcMSgRUaP065aImVMGIjZGARAINMvX5mPdtkLD76VqAs5qLzSua0JEEcagRESN1qNTAn5+dRYS6qwDt3LzIazcVABhcKjx+TU4qryGXpOI6FwxKBHROUlLsWHW1Cyk2E8tebJu+zEsW5MHVTO2u8zjU+HkRJREFEEMSkR0zlLsVsyamoUuHWzBbVv2nsQbn+6Fx8DH+4UITEvg8nB+JSKKDAYlImqSBJsFP796EPp0sQe37T1cjlc+3oVKl8+w+2iagNPlg9evMSwRUYtjUCKiJrNaTLh90gBk9UkJbjtaXIUXPtiBkwauD6eqAhWVXvj8fBKOiFoWgxIRhcVskvGTS/pjTFbn4LZShwcvfLDT0Fm8/aqGikovVJVPwhFRy2FQIqKwybKEq8f2xpWjega3VXv8eHXFLvyQX2LYfXyqhopqDzRmJSJqIVEXlPLy8nDHHXdg+PDhyMnJwcKFC+H1nvkR4RMnTmDhwoW45pprMGLECIwfPx4PPPAAjh492kKlJiJJkjB+eFfcPLEfFDkwmMivCiz9fB/WfH/UsOkDvD4NFVUeCDAtEVHzi6qgVFFRgRkzZsDn82HRokWYM2cO3nnnHTzxxBNnPG/nzp1YvXo1Jk2ahOeeew5z587F3r17ceONN6K0tHmWWSCihg3r1xF3XHVqYkoA+PSbw3hvbT78qjFjjDxeFY4qPglHRM3PFOkC1PX222+jqqoKixcvRlJSEgBAVVXMnz8fs2bNQufOnRs877zzzsPKlSthMp16OyNHjsSECRPw/vvv484772yJ4hNRjfSudsy+ZjCWrPwRpU4PAOC7PcUocXhw62UZsFnD/9Hj8apwefyIjTGBE3gTUXOJqhaltWvXIjs7OxiSAGDSpEnQNA3r16/XPc9ut4eEJABIS0tDSkoKTpwwfi0qIjq71KRYzL52MHp1TghuO3DMgeff34ETBjwRpwmBSpcfKgcsEVEziqoWpfz8fFx//fUh2+x2O1JTU5Gfn39O1zpw4ABKSkrQt2/fsMqkKFGVJVul2jpkXYanpepRA6AoEiQD8kdivAV3XzMI//4qD1v3ngQAlDjceOH9Hbjlsv7I7Jkc1vUFBFwePxLjY85+cB38TBqHdWkM1qNxjO6Sj6qg5HA4YLfb621PTExERUVFo68jhMBjjz2GTp06YfLkyWGVyW6PDet8OoV1aYzmrkeXxwdVkgztzrr7uqFYufEgPlwb+B8et1fF3z/5EdMm9MOlo3qGtUSJBAAmBckJ1nM+l59J47AujcF6jD5RFZSMsmjRImzatAmvvPIKbDbb2U84A4fDBdWgAajtlaLIsNtjWZdhaql69Po1OBxuaAYP/MnJ6ozEWBPe/s9++PwahACWfbkf+UfLcf1F6TCblLNfRIfDIaGy0o2EWHOjAh4/k8ZhXRqD9WicxMRYyLJxLXNRFZTsdjuczvoT1FVUVCAxMbFR13jnnXfw7LPP4vHHH0d2dnbYZVJVDX7OBmwI1qUxmrseVVWDqgrDgxIADOyVgtnXZOGNT/egvDIw7cfWvSdRVOrCrZdlIDnh3LrQaqkQcFR6oakCcdbGD+7mZ9I4rEtjsB7DZ/SPrqjqDE1PT683FsnpdKK4uBjp6elnPX/16tV49NFHce+99+KGG25ormISURi6dIjDL68bgj5dTg3yLjxZhWff+wH7jza+i/10mibgrPaiosoLcI4lIjJIVAWl8ePHY8OGDXA4HMFtq1atgizLyMnJOeO5mzdvxv33348bb7wRubm5zV1UIgpDfKwZd04eGLLsSbXbj79/sjusySmFCFynzOmBT9XCGvtERAREWVCaPn064uLikJubi6+//hrLli3DwoULMX369JA5lGbMmIHLLrss+DovLw+5ubno3bs3rrnmGnz//ffBr0OHDkXirRDRWSiyjKk5fXDDhL4wKYFAI0Rgcso3V++F2+tv8rU9Pg3lTg+c1Wee1Z+I6GyiaoxSYmIilixZggULFiA3NxdxcXG44YYbMGfOnJDjNE2DqqrB19u2bYPT6YTT6cRPfvKTkGOvu+66s87sTUSRMzIjFZ1TbHjzs1PjlnYdLMPi5T/glksz0LVjXJOuq2oCVW4ffH4VcbEWxJgVw5ZRIaL2QxL8yXFGZWVVHFgXJpNJRnJyHOsyTC1Vjz5VQ5nD0yyDuc+k2u3Dv77Yj31HTo1TMikSpoztjQsGdAqrG02WJcRaFMTbzJAg8TNpINalMViPxklJiTN0PqqoalEiovbLZjXjp1dk4uMNB/HN7sCM+n5V4P11B5Bf6MC1F/aBxazg2MkqVLv9sFlN6NIxDnIjApSmCVS5/fD6NcTFmhFvauZJO4XAoSInKqt9iLeZ0bNzQqPK2dr4NQ2bfyhCtUeDLUbGeZmpMBn4WDZRNGBQIqKokHekHF9tK8TJchesFgUerxp8dm17XgkOHnPAHmeBs9oLVQMUGeiYFIsJw7qib/ekRt3D59fgqPRCVTXExZ/7BJWNsftgKVZsKsDx0mqoqoCiSEhLsWHymF4Y2DulWe4ZCas2F2DFxgK4PH4IBCb+/GeMCZOze+HK0b0iXTwiwzD6E1HE5R0px3tfH8Dx0mpYzAqSEmKQbI+BLJ9qhXFU+3CkuAp+VSAu1gSLWcHxUhfe+/oA8o6UN/pemhCo9vhRUuFGldsHQBi25MHug6VY8ukeHCmuRIxZgT0+MDbqSHEVlny6B7sPlhpzowhbtbkAy9bko8rthyxLMCkSZFlClduPZWvysWpzQaSLSGQYBiUiiihNCHy1rRAenwq7zQKzSYEsSbBaTOiUZIVZCU0xVW4/yp1eyLIMu80Mj0/DV9sKz3lMlSYEnFVelDg8qPb4ISDCGgelCYEVmwrg9vqRFB8DiznwPixmBUnxFri9KlZsKmjxsV9G82saVmwsgKoJmGsCkizJkGUJZkWCqgms2FgAv8ZxNtQ2MCgRUUQdO1mFk+Uu2GJM9YKKLMuIjzUDCF3o0uNTUVxWDY9PhS1GwclyF46drDrnewvUdMdVeXGywg1HlQdaE3PMoSInjpdWI85qrvc+JElCnNWE46XVOFRUf/WB1uSbXUVwefwwyVKD79MkS3B5/PhmV1GESkhkLAYlIoqoarc/MOZI7ymVmt/FiTYzYiyn1oPTBFDq8KDS5YdfFah2N33eJSEAVQ0M+C5zuuHxq+fculRZ7YOqCph0BoqbTDJUVaCy2tfkckaDUoc7MHZMr3qkQAAtdbhbrlBEzYhBiYgiymY1QZGhvxCoCPxOliQJKQkxSIyzhOx2efxwefxwVBkzuaTPr6Gi0ovymsCkGwhOE28zQ1Ek3Ue7/X4NiiIh3mY2pJyRkmK3BqpEr+Wt5t8rxd48g+WJWhqDEhFFVJeOceiYFItqj1pvQkghBLx+DVaLAm9NkIqLNSM1KRbmOi1QmgCWr8vHp98cgt+Aldc1TcDlVVHh9KLM4YaqnX3Ad8/OCUhLsaHK7W/wfVS5/UhLsaFn5wSdK7QOowZ1RmyMCX5NNPg+/ZpAbIwJowZ11rkCUevCoEREESVLEiYM64oYswxHdWAmbU0I+PwqHNU+xFgUTBjRFTFmJbhfUSQkxpthrtPNJQSw5vtCPLv8BxwtrjSkbJoQ8Pg0lDrdqPb4gTMM+JYlCZPH9ILVoqC80guvL/A+vD4V5ZVeWC0KJo/p1ernUzLJMiZn94IiS/CpApomoAkNmibgUwUUWcLk7F6cT4naDM7MfRacJTV8nHHWGG19Zu668yg1NE+S3v6BPZOwcecJlNQZEyNLwPjh3TBxZDeYGhj7pCgSEhNtqKgIzHXUGJIUWJ/OYpJhsSiwmGQoslyvVaU9z6MUy3mUmow/J41j9MzcDEpnwQ9t+PgDwBhtPSgBgRacM828rbff61Ox6ptD2LQz9Emr1KRYXH9Rer3urqYEpbokKbAsitWsICHOUm+8Tnuamfu/e4o5M7cB+HPSOAxKLYwf2vDxB4Ax2kNQCld+YQWWrclHmdMT3CYBGJOVhssv6BF8ai7coBS8thRoRbE3EJbaC35/G4P1aByjgxKjPxG1GeldE/HrG4Zi7OC04MNqAsDGncfxt3e3YXdBmaH3EyLw1F1FpbdmwHfbazUiau8YlIioTbGYFUwZ2xt3T81CalJscHtFlRdvfLoH//xsD8orPWe4wrmpDUulzsCElaoqQpZeIaLWjUGJiNqkXmkJuOf6IZg4shuUOsFl18Ey/GXp9/hsc4H+3E1NUDthZYnTjRKHG+6ap97YykTUujEoEVGbZVJkXHp+D9xz/VD0Tjs1oNvr17D8y/14+t3tyC+sMPSemibg8aqoqPSgxOFGuTMQmoxcfJeIWg6DEhG1eZ2SY3HX1YMwbXw6bDGm4PaiMhde+Xg33v7PPkO744BTy6K4akLTyQoPHFVeuL3+4ASWDE5E0c909kOIiFo/WZJw/oBOGNQ7GZ99exjf7j4RfFBte14JdheU4aLhXXHh0K4hE1kaQQjAr2rwq1pNQJKgyBIsJhlmswKzIsOkSMFjiSh6MCgRUbtis5px/YS+uPiCnnhz5W4cKa4CEFjj7fPvjuC7H0/gilE9MbRvh2YZXyREYKkPTRPw+TVIHj9kSYKiSLCYFFjMMswmGbJUfzJLImp57HojonapT9dE5F4/BNPGpyPOeur/GcsrvfjXF/vxwgc7UXDc2ezlEAJQNQGvT0Oly4cypwcnK9wo49gmoqjAFiUiardqu+MGp6fgi/8excadx6FqgVacwycq8eKHO5HVJwVXXNADHetMNdCcasc2qaoKj09FlSwjxqIgxqLAbJIhgd1zRC2JQYmI2j2rxYSrsnth1KBOWLX5EHYdPDUx5c4Dpdh9sBTnD+iEied1h91mabFyCRGYKd3n0lDt8UORJcSYFVjMgbXmJImhiai5MSgREdXomBiL2y7PRH6hAys3F+BozfglTQDf7D6BrXtPYuyQNIwf1hWxMS3741PTTo1rkt3+wFpzlkBoMpskyJIMQDA4ERmMQYmI6DTpXe34xbWD8UNeCVZ/exilNWvH+VQNa74vxOZdRRg3tAvGDk6D1dLyP0Y1IaCpApWuQGiSZMAkSzWhiYPBiYzEoERE1ABZkjCsX0dk9UnBN7tP4MutR1Hl8gEA3F4Vn393BOt/OI4Lh3ZBdlZacMHdlqYJAaiBcU0enwZZkiDLEswmGRazDEmSIEtSsNWJ4Yno3DAoERGdgUmRMXZwGs7LTMWGH45j3fZCuL0qgMAab599exjrth/DuCFdkD24c0RamOqqbW3yqxpcnsCklhIC4clikmGxBMY3mZRAaGJuIjozBiUiokaIMSu4eGQ3jMnqjK+3H8P6Hcfg9QXWinN5/Fj93WGs216I7MFpGDs4DXFWc4RLHCAEIFAnPHkD8zaZFAkxFlMwNEmSxNYmogYwKBERnYPYGBMuu6AHcoak4evtx7Bh5/FgYHJ7VXy55Si+3n4MowZ0Qs7QLkiKj4lwiUMJAahCQNUEPD4vZPlU11xgjJMCsyIBgYkI2OJE7R6DEhFRE9isZlw+qifGDe2C9T8cx4Ydx+HxBbrkfH4N63ccx8adRRjevwPGDe2KtBRbhEvcME0T0CDgVwNBr3aMkyJLUBQZZkWCySRDqQlUDFDU3jAoERGFwWY147ILemDc0C7YtLMI63ccQ7XbDyAwXmjL3pPYsvck+ndPxLihXdCvW2KzLI1ihLqtTYFh62pwbbrasU6KEhgobqpZn07SpJpFfqXgDOIMUdSWMCgRERkgNsaEi0d2Q86QNHy35wS+3n4M5ZXe4P59Ryqw70gFOiXHImdwGob17wiLKTJPyp2L2rXpal7BrwIe76kAZVIk+ISEqko3JAAmRYJikgMtUzXhSZHZEkWtF4MSEZGBLGYFYwd3wehBnfFDXinWbS/EsZLq4P4TZS68t+4AVn1zCBcM6ITRgzojOcEawRI3TW2AUqVAy5nXr0JVAymotvUJ0qm/m2q68EyKDLMS6MoLtKwxPFF0Y1AiImoGiixjeP+OGNavA/IKHVi//Rj2HC4P7nd5VKzddgzrth1DZs8kjB7UGf17JNWMA2rdap+0QzAABVqiUKclSpbqdOPJgXFRkiRBhgRJDp3WQNS5Tu31iVoKgxIRUTOSJAn9uiWiX7dEFJe7sGHHcWzdWwyvP/CknADw46Fy/HioHMkJMbhgQCeMzExt0TXlWlJtS5R2WjceUKcVKvAi8BqALAcmAFXkUwPNZVmu+Xv9eyjyqYk1GaooXJLgxBlnVFZWBX/NDzRqGpNJRnJyHOsyTC1Vjz5VQ5nDE5jxuY1SFAmJiTZUVFQHu4taktvrx3/3FGPTriKUVLjr7ZclYECvZJyf2Qn9eyTVjPGJTpGoy4a69kRNa1OgFQowyYHB5nJNN1+gq682RAF1mrtQE81qnNrekt8C/DlpnJSUOCgNJegmYosSEVELs1pMyBnSBdmD05B3tAKbdhbhx0NlwV/MmgB2HSzDroNlSLCZMaJ/KkZmpKJTcmxkCx4lGuraq7MXUAEvAmGjNkjVfXovMHdUIDRpmoAaaOaqGXwuQarZH3zar7Y1SwqEMEmSgq1d7A5s+xiUiIgiRJYk9O+ehP7dk1BR6cF3e4rx3Y8nUFF16mk5Z7UPa7cVYu22QvToFI/h/TtiaN8OUTPzd7QLDVWnpxm1UddoqAVLllHTUiUHw5ciS1AkGbKCOnNOBe7LINV6MSgREUWBxPgYXHJed1w8ohv2HSnHf/cUY3dBGVTt1G/YwycqcfhEJVZsKEBGj0QM69cRA3slw2KO/mkGWrMGW7CCGetU2Ko355RcG6AC46mCYSrYOiU10A1I0YZBiYjqkWUJkkDo0I06arsdRN1uB3Hqdb1jpdO3oeZ8cGHW08iyhMyeycjsmYwqtw/f7zuJrXuLUVhnigFNiOAAcLNJxsBeyRjatwP6d0+C2WTc2Aw6Nw3NORUQ+EswSOFUoDLJgXFUZrMMc7UXXp8GoQlINQPY6z71J0n8fokEDuY+Cw6sCx8HKRqjJeux7kBu6fRxrg2EJ00LnCOEqAk/gQHJwKnZmoPdF3WuqWqBWaBVVYPfr8GnBq6hNfMvg0gP5m6KYyVV+H7fSXy//ySc1b4Gj4kxKxjYKxmD01NaLDS1xrqMRrX16HBUQ9NCx1aZZCkw4koASQmWOuOjqCEczE1EzU53Lh+dzYoMKE344S0rEswKIFkUoObJJVUVULXa4KTBr9YEJw3BINYedekQhy4d4nDFqJ7IL3Rg2/6T2HGgNLi+HAB4fCq+3x8IUxaTjIweSRjUJwUDeibBauGP+9ZAiMD6ezWvAAC1sdhk4C9/ajx+5xBRxNUdpxEYw6EgxhwIT4AItjz5a4KT369CrdOKFTJ8pI2TZQn9uieiX/dETB3XB/uOlGN7Xgl+LCgLzs0EAF6/hh0HSrHjQCkUWUKfLnYM7JWMAb2SkZwQE8F3QNS6MCgRUVSqG55kSappfZJruu0s0IQGTQNUTYOqicBj3qqAz69Bre0GbOPpyWySMah3Cgb1ToHXr2Lv4QrsyC/Bj4fK4PWdCk2qJrD/aAX2H63ARxsOIi3FhsyeScjsmYQenRKiep4mokhjUCKiVqU2QNU+VaTIgSe+givXA1D9Aj5Ng1/V4PGp8Pu1Nh+aLCYFg/ukYHCfFPj8GvKOVmDnwVLsLihDtdsfcuzx0mocL63Gmu8LYbUo6Nc9ERndk9C/eyIS49naRFQXgxIRtQl1g5CiSFAUBZJFQbzVjGqPH9VuP/xq+3iYwGySMaCmm03TBA6dcGL3wTLsLijDydNmAnd7VezIL8WO/FIAQGpSLPp1T0T/bono08WOGAunHqD2jUGJiNqs2vAUZzXBajHBWe2F2+s/80ltjCxL6J1mR+80OyaN6YWTFS78WFCOvYfLceCYI2SeJgAoLnehuNyFjTuOQ5aA7p3ikd7FjvSuieiZFg+LicGJ2hcGJSJq82qnK0iKt6DaI8Pj06D3YF9b1zExFuOGxmLc0C7w+FTkFzqw93A59h0pR6nDE3KsJoBDRZU4VFSJr74vhCJL6J4ajz5dEtC7ix09O8cjLpYzhFPbxqBERO2GEIAtxoSEOBnWWAt8Hh+qXX6oWtsfw9SQ2nmXBvZKBgCUONzYfyQw6Du/sAIuT+gSH6omUFDkREGRE/i+EJIEdOlgQ/+eKeiaEovuqfFIirdAaq8plNokBiUialdqZxGPjTEjMc6CWIsCj09DtdsHn09rN9MMNKSD3YoOg6wYPagzNE3gWEkV8godyC90oOC4M2TOJiBQl4Unq1F48tSs4Qk2M3p2SkCPzvHo0Ske3TrGcYkVatUYlIio3RIiMPux1azAapZR5faj2uPnDNMIjG3qlhqPbqnxGD+sK1RN4NjJKhw45sCBYw4UFDnrtTgBgUV8dx4sxc6DgcHhkgR0TrahW2pc4KtjPNJSbFxqhVoNBiUiIgCAhPhYM2IsJlS5fPCeWqgruIZXe+yeq6XIErp3ikf3TvG4cFhXaELgRJkLh084UVjiwv7DZSg5bYwTEKi72ukI/runGEBgXqzOKbHo0iEOXTvaamYdt3H2cIpK/FQSEdUQIrCuVlJ8DDShBde2UzURXIvO51ehaqLdhydZkpCWEmgpql3rraLSi8NFThw6UYnDJypxpLgyZOLLWpoQOFZSjWMl1diy99T25IQYpKXY0KWDDZ1TAl8d7FZOiEkRxaBERHQaIQITWtYuX2eSJZgsMmIlQJIsUDUNmgqoomZNOr8Gr6pBU9vvWnQAEB9rxsDeKRjYOwVAYM2y4goXjhZX4UhxJY4WV+FYSRX8Ol2bZU4Pypwe7C4oC24zKRJSk2LRKTkWnZJs6JQci9TkWHSwx0CR2X1HzY9BiYiokWpbkSRIUBRAwak16VRNwO31w+31w+8PLOTb3smyhM7JNnROtmFkRiqAQOtccbkLhSerAl8lVTh2srreQPFafvVU6xNQEtyuyBJS7FakJlmRmhSLjolWdEyMRYdEK+KsJj55R4ZhUCIiCkPtkiqyFJjY0mY1wetT4XL7A61MWvvtnmuIIge67NJSToUnIQTKnJ7gWKbjJYE/Sxxu3bqrDVzF5S4AZSH7rBYFHROtSLFbA0/yJVqRYo9BSoIV8TYzZIYoOgcMSkREBqn9pW4xKYhJMEETGjw+DT6fCq//1OK9FEqSAq1DKXYrBtV02wGAz6+huNyForJqnChzBb7KXSg9Q4ACAsuyHCmuwpHiqnr7TIqE5IQYJCdYa/4MfCXFxyAp3oL4WDNboygEgxIRUTOo7aKzmhXE1qyX5vcLeHx+uH01A8I1sIvuDMwmGV07xqFrx7iQ7X5Vw8kKd7BFqaTm7ycr3HB7G+7CO3WuQHG5G8Xl7gb3mxQJiXExSIy3ICneEvx7YpwF9rjAn7Ex7NprTxiUiIiaWW0WUhQJcSYz4mLNgSfpVAGvX4XHp0JVBYTWvgeDN5ZJkYPdd3UJIVDt8eNkuRsljsBXqcONUocHJQ43qt1nX+fPr4rgufr3l2C3BYJTgs0Mu82CBJsF8TYzEmxmJNgC22NjTOzmawMYlIiIWlBtaJIlCRaThBizjIRYM/z+QGjy+lT4VAEhAgPC2eDUeJIkIc5qRlyaGb3SEurtd3v9wSfrSh0elFV6UF7zuszp0R1Qfjq/KlDq9KDUWX/eqLpkSUJ8rAnxsYFwHF/zFRdrRpy1ZrvVDJvVBHu8BYL/2FGJQYmIKILqtjbFKibYrGYICKiqgKpp8NfM3VT7JB3DU9NZLSZ06WBClw5xDe53efyoqPKivNKD8koPKiq9ga8qLyqqPHBUeXWnNmiIJgQc1T44qn2NOt6kyLBZFdhiAuHJFhN4OCC25k+7zYLswZ1ht8U0ugwUPgYlIqIoUtuqoMgSFFmBxQRIVjOAQFedXzs1d5NfrQ1ODE9GiI0JhJLTu/RqCSHg8qhwVHvhqKr5qvbCWe2Ds86flS7fOQWqWn5Vg6NKg6NKP1it/u4w/nj3GJhNXD+vpURdUMrLy8Njjz2GrVu3Ii4uDtdccw3uu+8+WCyWM54nhMDLL7+Mt956C6WlpRg4cCAeeughDB8+vGUKTkTUTE4PT1azAk0AR4srUeXyIcaioGNSLFT1VHjy+lV8/X0hTpa7kGy34sLhXWFWTv1y1URg7bZqtx82qwldOsaFjKc5234AUDUN2/JK4fZqsFpkDO6Tck6TQBpRhnDvoWoatu8/ifJKL5LiLRjar6Pue5AkKdDSYw0NU6ffI62DDV6fhkqXLxieCk9W1bRIBRZernL7UeXyocrth8tz9rFTtUodHuQdrUBGz2SOf2ohkoiiTtGKigpMnjwZvXv3xqxZs1BUVIQnnngCU6dOxbx588547ksvvYRnnnkGv/nNb5CZmYk333wTGzZswAcffIAePXo0uUxlZVXw++tPwU+NZzLJSE6OY12GifVonNZel7sPlmLFpgIcL62GqgpYzDK6pcZj0uieyOyZjH9/tR9b950MtjIFpiXQkNUnBZNG90L+0Qqs2V6I4jIXfKqABIGOSbGYMKwr+nZPQt6Rcny1LRCyVA1QZITsB4B1247iq62FcHvV2pVeYLUomDCiKy4c1u2s7+Fs92hMGcK9R7jvwaj3oWqBpyFlswlFxU44q/1wuX2o9vhR7fbjeGk1Ck9WQdUEEmLNUIVAWooNk8f0Cs6CTqekpMRBUYybtT2qgtKLL76IF154AV9++SWSkpIAAP/6178wf/58fPnll+jcuXOD53k8HowdOxa33nor7r//fgCA1+vFlVdeifHjx+PRRx9tcpla6w/SaNLafylFC9ajcVpzXe4+WIoln+6B2+tHnNUMk0mG36+hyu2H1aIgLSUW2/NKIUmBMS8mRYbZJEORJZgUGQN6JaHU6YHXpyLWYoKsSFBVDdVuFSZFwtD0DvhubzHcXj9iTAokWYLfr6LKrcKsSLh6bC8UllRh1ebD0ETgF78kBcZaqRogS8AVo3qcMWjkHSnHe18fgMenwhZjgqLIgTJ4VMSYZYwe2Ambd5/Q3X/duD5nDUtnu0ffrnb8d09x4D1ICKQkAaiice/B6PehKFJwzTy1Trdd3Xukpdhgs5pwstwd/PeecUUmw9JpjA5KUbVQztq1a5GdnR0MSQAwadIkaJqG9evX6563ZcsWVFZWYtKkScFtFosFl112GdauXducRSYiajGaEFixqQBurx9J8TGwmJXA03NmBUnxFlS5vNieVwogEFz8fg1ujx/OKi/KnB4Ul7uwbtsxHCwsh68mXFW5/PD6NJgUCZUuH9bvOA5VVdEpKRYpiVakJMSgU7INvdMSYI+Pwc6DpfjxUAXSOsShe2oc0jrEoXNyHFKTY5GaZIU9PgZb9p4EJAFZlk59SVIgUEFg3Y5j8KsakuJiYDEF3oPZpMBuM8PjVfHV1kJ4fCrsNgvMp+/3afhqW+EZ55/ShMBX2/Sv4fb48V1NSDLJCJZPliWYZEATwFdbC6Fq+iH6bPdojvdhUkL/vd1eFSs2FXAurmYWVWOU8vPzcf3114dss9vtSE1NRX5+/hnPA4D09PSQ7X379sWSJUvgdrthtVqbVKbExFgOkgxTbTc66zI8rEfjtNa69KsaHrjtfEhAgxMeqppo1CPmkiRBkeufL4SAqolgcKhPQBP1Fw1u6Mi4WHPNOngImRtKVTU8fPtoSJJU24hTrwxaTRnqHiFCjgHiYk26Y4lUVcMjd9TcQ6pzrqhzD1G74fTSi+AfcbFmWMwND5r2qxoevmO07r+FVud9NFDVNesGAvE2M0w1rR+yHFgAOOQet48OfF4l1IRNKTi7uxCBebcS4yzBa1CgHo0UVUHJ4XDAbrfX256YmIiKiooznmexWBATE/rIpN1uhxACFRUVTQ5KMlenNgzr0hisR+O0trqsfZKqNgCEo6HzazOWpLtfghBa4Jc8BKRTmeLUMTV/apposPtDVWsDROh5Us1/awPEqY1Snf215RCQpUCXYkOEQGjYO+29aBoATdQEnNPfaCCcCRGoZ70AIoSAohsoAVlIUKFBkWXdfytNBIJU3XpSlNDB7IpS/x6nXko1wVgytKuJQkVVUCIiIn0xFgWdLA0/ut5axFgUpDbze4ixKIixxJ79wLDuYUKMpXl/hSoWE6zNfA86u6iKoHa7HU6ns972iooKJCYmnvE8r9cLjyd0llSHwwFJks54LhEREZGeqApK6enp9cYiOZ1OFBcX1xt/dPp5AHDgwIGQ7fn5+ejatWuTu92IiIiofYuqoDR+/Hhs2LABDocjuG3VqlWQZRk5OTm6540cORLx8fFYuXJlcJvP58Nnn32G8ePHN2uZiYiIqO2Kqs7P6dOn44033kBubm5wwsmFCxdi+vTpIXMozZgxA4WFhVi9ejUAICYmBrNmzcKiRYuQkpKCjIwMLF26FOXl5Zg5c2ak3g4RERG1clEVlBITE7FkyRIsWLAAubm5iIuLww033IA5c+aEHKdpGlQ1dJXnn//85xBC4LXXXgsuYfLqq6+GNSs3ERERtW9RNTM3ERERUTSJqjFKRERERNGEQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHS026C0Zs0a3HbbbRgzZgwGDx6MSy65BH/84x/rrTX3xRdfYOrUqRgyZAiuuOIKLFu2LEIlbj2qqqowfvx4ZGZm4ocffgjZ9+677+KKK67AkCFDMHXqVHz55ZcRKmV0Wr58OTIzM+t9/fnPfw45jvXYOO+99x6uvfZaDBkyBKNHj8Zdd90Ft9sd3M/v77P76U9/2uBnMjMzEytWrAgex89k4/znP//BjTfeiBEjRmDcuHH49a9/jcOHD9c7jvV5Zl9++SWuu+46DB48GBdddBGeeeaZevMrAgZ9j4t26v333xdPPvmkWLVqldi0aZN44403xKhRo8Qdd9wRPObbb78VAwcOFL/73e/Exo0bxVNPPSUyMzPFypUrI1jy6Ldw4UIxduxYkZGRIbZv3x7c/vHHH4vMzEzx1FNPiY0bN4rf/e53YtCgQWLr1q2RK2yUWbZsmcjIyBBr164VW7duDX4VFhYGj2E9Ns5zzz0nRowYIV588UWxefNmsWrVKvH73/9eVFZWCiH4/d1Y+/btC/ksbt26Vdx3331i0KBBoqSkRAjBz2Rjbdq0SQwYMEDMnTtXrF+/XqxYsUJcfvnl4tJLLxUulyt4HOvzzLZu3SoGDBggHnjgAbF27Vrx2muviaFDh4onnngi5DijvsfbbVBqyL/+9S+RkZEhjh8/LoQQ4s477xQ333xzyDH333+/mDRpUiSK1yrs379fDB8+XCxdurReULr88svF/fffH3L8zTffLO66666WLmbUqg1Ktb+AGsJ6PLu8vDwxaNAg8dVXX+kew+/vpps4caL4+c9/HnzNz2Tj/O53vxMTJ04UmqYFt23cuFFkZGSIb7/9NriN9Xlmd955p7juuutCtr366qsiKytLFBcXhxxnxPd4u+16a0hSUhKAwIK6Xq8XmzdvxpVXXhlyzFVXXYW8vDwcOXIkAiWMfo899himT5+OPn36hGw/fPgwDh48iEmTJoVsv+qqq7Bx40Z4vd6WLGarxXpsnOXLl6N79+646KKLGtzP7++m27JlC44cOYKrr74aAD+T58Lv9yMuLg6SJAW3JSQkAABEzSIZrM+z2717N3JyckK2jRs3Dj6fD19//TUAY7/H231QUlUVHo8HO3fuxLPPPouJEyeie/fuOHToEHw+H9LT00OO79u3LwAgPz8/EsWNaqtWrcLevXuRm5tbb19tfZ0eoPr27Qufz9dgH317NmXKFAwcOBCXXHIJXnzxxWDfO+uxcbZt24aMjAw899xzyM7OxuDBgzF9+nRs27YNAPj9HYaPP/4YNpsNl1xyCQB+Js/FtGnTkJeXhzfffBNOpxOHDx/GX//6VwwaNAgjR44EwPpsDI/HA4vFErKt9nVeXh4AY7/Ho2pR3Ei4+OKLUVRUBAC48MIL8Ze//AUAUFFRAQCw2+0hx9e+rt1PAS6XC0888QTmzJmD+Pj4evtZn42TmpqKe+65B8OGDYMkSfjiiy/wt7/9DUVFRZg3bx7rsZGKi4uxY8cO7N27F7///e8RGxuLF154AXfeeSc+++wz1mMT+f1+rFy5EhMnToTNZgPA7+1zcf7552Px4sV44IEH8H//938AgIEDB+KVV16BoigAWJ+N0atXL2zfvj1k2/fffw/gVP0YWY/tPii99NJLcLlc2L9/P55//nnMnj0bf//73yNdrFbn+eefR4cOHXD99ddHuiit2oUXXogLL7ww+HrcuHGIiYnBkiVLMHv27AiWrHURQqC6uhpPP/00BgwYAAAYNmwYJk6ciH/+858YN25chEvYOq1fvx6lpaWYMmVKpIvSKm3ZsgW//e1vcdNNN2HChAkoLy/Hc889h7vvvhtvvfUWrFZrpIvYKtxyyy145JFHsGTJElxzzTXYv38//va3vwXDptHafdfbgAEDMGLECNx444147rnnsHnzZqxevRqJiYkAUG+6AIfDAQDB/QQcPXoUr732Gu699144nU44HA5UV1cDAKqrq1FVVcX6DMOkSZOgqip2797Nemwku92OpKSkYEgCAmMQBw0ahP3797Mem+jjjz9GUlJSSNBkXTbeY489hjFjxmDu3LkYM2YMrrzySrz00kvYtWsXPvjgAwCsz8aYNm0aZsyYgYULF2L06NG4/fbbMX36dCQmJqJTp04AjK3Hdh+U6srMzITZbMahQ4fQs2dPmM3mev2Yta9P7/dsz44cOQKfz4e7774bF1xwAS644IJg68fPfvYz3HHHHcH6aqg+zWYzevTo0eLlbo1Yj43Tr18/3X0ej4ff303gdrvx+eef48orr4TZbA5u52ey8fLy8kLCOwCkpaUhOTkZhw4dAsD6bAxZlvHwww9j06ZN+OCDD7BhwwbcdNNNKC0txbBhwwDA0O9xBqU6tm3bBp/Ph+7du8NisWD06NH49NNPQ4755JNP0LdvX3Tv3j1CpYw+AwcOxOuvvx7y9dBDDwEA5s+fj9///vfo0aMHevfujVWrVoWc+8knnyA7O7vewDw65ZNPPoGiKBg0aBDrsZEuvvhilJeXY/fu3cFtZWVl2LlzJ7Kysvj93QRffPEFqqurg0+71eJnsvG6du2KXbt2hWw7evQoysrK0K1bNwCsz3ORkJCAAQMGwG6344033kD37t0xduxYADD0e7zdjlH61a9+hcGDByMzMxNWqxU//vgjXn31VWRmZuLSSy8FAPziF7/Az372Mzz66KOYNGkSNm/ejI8//hhPPfVUhEsfXex2O0aPHt3gvqysLGRlZQEA7rnnHvzmN79Bz549MXr0aHzyySfYvn07/vnPf7ZkcaPazJkzMXr0aGRmZgIIzOL7zjvv4Gc/+xlSU1MBsB4b49JLL8WQIUNw7733Ys6cOYiJicFLL70Ei8WCW265BQC/v8/VRx99hK5du+K8886rt4+fycaZPn06/vCHP+Cxxx7DxIkTUV5eHhzfWXc6ANbnmW3fvh3ffPMNBg4cCLfbjS+++AIffPABXn755ZBxSoZ9jzd1wqfW7sUXXxTXXHONGDFihBg+fLiYPHmy+Nvf/iacTmfIcZ9//rmYMmWKyMrKEpdddpl49913I1Ti1mXTpk31JpwUQoh33nlHXHbZZSIrK0tMmTJFfPHFFxEqYXRasGCBuPzyy8XQoUPF4MGDxZQpU8SSJUtCJqgTgvXYGCUlJeI3v/mNOO+888TQoUPFnXfeKfbt2xdyDL+/G6e8vFxkZWWJhQsX6h7Dz+TZaZom3nrrLXH11VeL4cOHi5ycHJGbmyv2799f71jWp75du3aJG2+8UQwfPlwMHz5czJgxQ2zZsqXBY434HpeEqJnlioiIiIhCcIwSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIR7tdwoSIGq92SZWzef3113WXs2nr3nzzTcTGxmLatGmRLgoRGYhBiYjOauHChSGvP/jgA6xfv77e9r59+7ZksaLK0qVLkZyczKBE1MYwKBHRWV1zzTUhr7dt24b169fX295WCCHg8XhgtVpZDqJ2jmOUiMgQmqbhH//4ByZPnowhQ4Zg7NixmDdvHioqKkKOmzhxImbNmoXNmzdj2rRpGDp0KK6++mps3rwZAPDZZ5/h6quvxpAhQzBt2jTs2rUr5Py5c+dixIgROHz4MGbOnInhw4dj3LhxWLx4MU5fuvJcy7Ru3bpgmd5++20AwLJly/Czn/0M2dnZGDx4MK666iq89dZb9c7ft28fvvnmG2RmZiIzMxM//elPAQCLFi1qsOty+fLlyMzMxJEjRxpVDofDgccffxwXXXQRBg8ejMsuuwwvvfQSNE1r9L8REZ07tigRkSHmzZuH9957D9OmTcNPf/pTHDlyBG+++SZ27dqFpUuXwmw2B48tKCjAAw88gOnTp2Pq1Kl47bXXMHv2bMyfPx9PPfUUfvKTnwAAXnrpJdx3331YtWoVZPnU/9epqoq77roLw4YNw//8z/9g3bp1WLRoEVRVxa9//esmlenAgQN44IEHcPPNN+Omm25Cnz59AAS61Pr374+JEyfCZDLhyy+/xPz58yGEwK233goAePjhh7FgwQLYbDbMnj0bANCxY8cm1WND5XC5XLjttttQVFSE6dOno0uXLti6dSv++te/ori4GI888kiT7kVEjSCIiM7R/PnzRUZGRvD1t99+KzIyMsSHH34YctzatWvrbb/44otFRkaG2LJlS3DbunXrREZGhhg6dKg4evRocPvbb78tMjIyxKZNm4LbHnzwQZGRkSEWLFgQ3KZpmrj77rtFVlaWKCkpaXKZ1q5dW++9ulyuetvuvPNOcckll4Rsmzx5srjtttvqHfvMM8+E1FWtZcuWiYyMDHH48OGzluPZZ58Vw4cPFwcOHAjZ/uc//1kMHDhQFBYW1rs+ERmDXW9EFLZVq1YhISEBOTk5KC0tDX5lZWXBZrMFu9Vq9evXDyNGjAi+HjZsGABgzJgx6Nq1a73thw8frnfP2tYcAJAkCbfeeit8Ph82btzYpDJ1794dF154Yb371B0f5HQ6UVpailGjRuHw4cNwOp2NrqPGaqgcq1atwnnnnQe73R7yXsaOHQtVVfHtt98aXg4iCmDXGxGFraCgAE6nE9nZ2Q3uLykpCXndpUuXkNcJCQkAgLS0tJDt8fHxAALjc+qSZRk9evQI2VbbVXb06NEmlal79+4NHvff//4XixYtwvfffw+XyxWyz+l0BstulIbKUVBQgD179ui+l9LSUkPLQESnMCgRUdg0TUOHDh3w5z//ucH9KSkpIa8VRWnwOL3t4rRB2s1RpoaeLDt06BBuv/12pKenY+7cuejSpQvMZjPWrFmDf/zjH40aSC1JUoPbVVVtcHtD5dA0DTk5ObjrrrsaPKd3795nLQcRNQ2DEhGFrWfPnti4cSNGjhzZIo+ya5qGw4cPB1uRgMAgaADo1q2bYWX64osv4PV68fzzz4d0CZ7ebQfoByK73Q4g0CpW+3cAKCwsbHQ5evbsierqaowdO7bR5xCRMThGiYjCNmnSJKiqiueee67ePr/fX6/rzAhvvvlm8O9CCLz55pswm83B7ikjylTbwlW3RcvpdGLZsmX1jo2NjW3wmj179gSAkHFE1dXVeP/99896/1qTJk3C1q1bsW7dunr7HA4H/H5/o69FROeGLUpEFLZRo0bh5ptvxosvvojdu3cjJycHZrMZBw8exKpVq/DII4/gyiuvNOx+MTExWLduHR588EEMHToU69atw1dffYXZs2cHu9SMKFPtObNnz8b06dNRVVWFd999Fx06dEBxcXHIsVlZWVi6dCmee+459OrVCykpKcjOzkZOTg66du2KRx55BPn5+VAUBcuWLUNycnKjW5VmzpyJL774ArNnz8Z1112HrKwsuFwu7N27F59++in+85//1OtKJCJjMCgRkSH+7//+D4MHD8bbb7+Np556CoqioFu3bpg6dSpGjhxp6L0URcErr7yCRx99FH/6058QFxeHX/3qV8jNzTW0TOnp6XjmmWfwt7/9DU8++SQ6duyIn/zkJ0hJScHDDz8ccmxubi4KCwvxyiuvoKqqCqNGjUJ2djbMZjMWL16M+fPn4+mnn0ZqaipmzJgBu92Ohx56qFHvNzY2Fm+88QZefPFFrFq1Cu+//z7i4+PRu3dv3HPPPYYPKCeiUyTRlFGSREQRMnfuXHz66afYunVrpItCRO0AxygRERER6WBQIiIiItLBoERERESkg2OUiIiIiHSwRYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIx/8DdmKQfBcfotQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "sns.set(color_codes=True)\n", "plt.xlim(30,90)\n", "plt.ylim(0,1)\n", "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] } ], "metadata": { "celltoolbar": "Hide code", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.13" } }, "nbformat": 4, "nbformat_minor": 2 }