From 222f5bd6e7b3e3ce1f5af49fa223c2ab453fcdc9 Mon Sep 17 00:00:00 2001 From: fac8f4cc02bfb4e3dfb02802cea5d62a Date: Mon, 22 Jan 2024 17:49:09 +0000 Subject: [PATCH] Upload New File --- ..._Python3_challenger_mooc_environment.ipynb | 679 ++++++++++++++++++ 1 file changed, 679 insertions(+) create mode 100644 module4/src_Python3_challenger_mooc_environment.ipynb diff --git a/module4/src_Python3_challenger_mooc_environment.ipynb b/module4/src_Python3_challenger_mooc_environment.ipynb new file mode 100644 index 0000000..d19e2db --- /dev/null +++ b/module4/src_Python3_challenger_mooc_environment.ipynb @@ -0,0 +1,679 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this document we reperform some of the analysis provided in \n", + "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n", + "\n", + "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Technical information on the computer on which the analysis is run" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.12.1 | packaged by conda-forge | (main, Dec 23 2023, 08:03:24) [GCC 12.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')\n", + "IPython 8.20.0\n", + "IPython.core.release 8.20.0\n", + "PIL 10.2.0\n", + "PIL.Image 10.2.0\n", + "PIL._deprecate 10.2.0\n", + "PIL._version 10.2.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", + "cffi 1.16.0\n", + "comm 0.2.1\n", + "csv 1.0\n", + "ctypes 1.1.0\n", + "cycler 0.12.1\n", + "dateutil 2.8.2\n", + "debugpy 1.8.0\n", + "debugpy.public_api 1.8.0\n", + "decimal 1.70\n", + "decorator 5.1.1\n", + "defusedxml 0.7.1\n", + "executing 2.0.1\n", + "executing.version 2.0.1\n", + "http.server 0.6\n", + "ipaddress 1.0\n", + "ipykernel 6.29.0\n", + "ipykernel._version 6.29.0\n", + "ipywidgets 8.1.1\n", + "ipywidgets._version 8.1.1\n", + "jedi 0.19.1\n", + "json 2.0.9\n", + "jupyter_client 8.6.0\n", + "jupyter_client._version 8.6.0\n", + "jupyter_core 5.7.1\n", + "jupyter_core.version 5.7.1\n", + "kiwisolver 1.4.5\n", + "kiwisolver._cext 1.4.5\n", + "logging 0.5.1.2\n", + "matplotlib 3.8.2\n", + "matplotlib._version 3.8.2\n", + "numpy 1.26.3\n", + "numpy.core 1.26.3\n", + "numpy.core._multiarray_umath 3.1\n", + "numpy.linalg._umath_linalg 0.1.5\n", + "numpy.version 1.26.3\n", + "packaging 23.2\n", + "pandas 2.2.0\n", + "pandas._version_meson 2.2.0\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", + "platform 1.0.8\n", + "platformdirs 4.1.0\n", + "platformdirs.version 4.1.0\n", + "prompt_toolkit 3.0.42\n", + "psutil 5.9.8\n", + "ptyprocess 0.7.0\n", + "pure_eval 0.2.2\n", + "pure_eval.version 0.2.2\n", + "pydevd 2.9.5\n", + "pygments 2.17.2\n", + "pyparsing 3.1.1\n", + "pytz 2023.3.post1\n", + "re 2.2.1\n", + "scipy 1.12.0\n", + "scipy._lib._uarray 0.8.8.dev0+aa94c5a4.scipy\n", + "scipy._lib.array_api_compat.array_api_compat 1.4\n", + "scipy._lib.array_api_compat.array_api_compat.numpy 1.26.3\n", + "scipy._lib.decorator 4.0.5\n", + "scipy.integrate._dop 1.26.3\n", + "scipy.integrate._lsoda 1.26.3\n", + "scipy.integrate._vode 1.26.3\n", + "scipy.interpolate.dfitpack 1.26.3\n", + "scipy.linalg._fblas 1.26.3\n", + "scipy.linalg._flapack 1.26.3\n", + "scipy.linalg._flinalg 1.26.3\n", + "scipy.linalg._interpolative 1.26.3\n", + "scipy.optimize._cobyla 1.26.3\n", + "scipy.optimize._lbfgsb 1.26.3\n", + "scipy.optimize._minpack2 1.26.3\n", + "scipy.optimize._slsqp 1.26.3\n", + "scipy.sparse.linalg._eigen.arpack._arpack 1.26.3\n", + "scipy.special._specfun 1.26.3\n", + "scipy.stats._mvn 1.26.3\n", + "seaborn 0.13.1\n", + "seaborn.external.appdirs 1.4.4\n", + "seaborn.external.husl 2.1.0\n", + "six 1.16.0\n", + "socketserver 0.4\n", + "stack_data 0.6.2\n", + "stack_data.version 0.6.2\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.14.1\n", + "traitlets._version 5.14.1\n", + "urllib.request 3.12\n", + "wcwidth 0.2.13\n", + "xmlrpc.client 3.12\n", + "zlib 1.0\n", + "zmq 25.1.2\n", + "zmq.sugar 25.1.2\n", + "zmq.sugar.version 25.1.2\n" + ] + } + ], + "source": [ + "def print_imported_modules():\n", + " import sys\n", + " for name, val in sorted(sys.modules.items()):\n", + " if(hasattr(val, '__version__')): \n", + " print(val.__name__, val.__version__)\n", + "# else:\n", + "# print(val.__name__, \"(unknown version)\")\n", + "def print_sys_info():\n", + " import sys\n", + " import platform\n", + " print(sys.version)\n", + " print(platform.uname())\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import statsmodels.api as sm\n", + "import seaborn as sns\n", + "\n", + "print_sys_info()\n", + "print_imported_modules()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Loading and inspecting data\n", + "Let's start by reading data." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "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": 4, + "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": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwQ0lEQVR4nO3de1yUdf7//+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==", + "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": 7, + "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[7], 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/rech_repro/lib/python3.12/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/rech_repro/lib/python3.12/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.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "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[8], 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/rech_repro/lib/python3.12/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/rech_repro/lib/python3.12/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": 9, + "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[9], line 3\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[1;32m 2\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), \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 data_pred\u001b[38;5;241m.\u001b[39mplot(x\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m\"\u001b[39m,y\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFrequency\u001b[39m\u001b[38;5;124m\"\u001b[39m,kind\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mline\u001b[39m\u001b[38;5;124m\"\u001b[39m,ylim\u001b[38;5;241m=\u001b[39m[\u001b[38;5;241m0\u001b[39m,\u001b[38;5;241m1\u001b[39m])\n\u001b[1;32m 5\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[0;31mNameError\u001b[0m: name 'logmodel' is not defined" + ] + } + ], + "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": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAG8CAYAAADU/gGwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABUsUlEQVR4nO3deXwU9f0/8NfM7JVrcxBuCEiQhPtQhCAgSAWBAJ6VVpAqighCq0WL9LBVWlp/fIuFooJXUVGkVbGCQvEAL9BWuS8lHOEMuTfJ3jOf3x+bLLMmgWQzyeZ4PR8PlMz5mTdL8uIzn/mMJIQQICIiIiIAgBzpBhARERE1JgxHRERERDoMR0REREQ6DEdEREREOgxHRERERDoMR0REREQ6DEdEREREOgxHRERERDoMR0REREQ6EQ9HJ0+exO9+9ztMmTIFvXr1QmZmZo33feedd3DjjTeib9++yMzMxAcffFCPLSUiIqKWIOLh6Pvvv8f27dvRpUsXpKam1ni/zZs3Y+HChbjhhhvw/PPPY+jQoXjooYfw+eef12NriYiIqLmTIv1uNU3TIMuBjLZw4ULs378fGzduvOx+48ePR48ePfC3v/0tuGzmzJkoKSnB+vXr6629RERE1LxFvOeoIhjVxqlTp3Ds2LFKt+AyMzOxd+9eFBQUGNU8IiIiamEiHo7CcezYMQBAt27dQpanpqZCCBFcT0RERFRbTTIcFRcXAwDsdnvI8vj4+JD1RERERLXVJMNRBUmSQr6uGD71w+W1EeEhWERERBRhpkg3IBz6HqLk5OTgcofDAaByj1JtSJIEh8MFVdXq1sgWTlFk2O1RrGUdsY7GYS2Nw1oag3U0Tnx8VFhjmKvTJMNRxVijY8eOhTz+n5WVBUmSKo1Fqi1V1eD384NqBNbSGKyjcVhL47CWxmAd687omz5N8rZa586d0a1bN7z//vshyzdu3Ih+/fohKSkpQi0jIiKipi7iPUculwvbt28HAJw5cwalpaXYvHkzAOCaa65BUlISFi1ahA0bNuDgwYPB/ebPn4+HHnoIKSkpGDZsGD766CN88cUXeOGFFyJyHURERNQ8RDwc5efn4+c//3nIsoqvX3nlFQwZMgSapkFV1ZBtxo8fD7fbjeeeew4vvvgiunTpgmXLlmH48OEN1nYiIiJqfiI+Q3ZjVFhYxvu/dWQyyUhMjGEt64h1NA5raRzW0hiso3GSkmKgKMaNFGqSY46IiIiI6gvDEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkU6jCEfHjx/HzJkzMWDAAGRkZGDx4sVwu92X3c/pdGLp0qX40Y9+hP79+2Ps2LFYsWIFvF5vA7SaiIiImiNTpBvgcDgwY8YMdOjQAcuXL0dBQQGWLFmCoqIiLF269JL7/v73v8eHH36Ihx56CFdeeSX27t2L5cuXo7i4GL/5zW8a6AqIiIioOYl4OFq3bh0cDgc2bNiApKQkAICiKFiwYAEeeOABpKamVrmf3+/H5s2bce+992L69OkAgKFDh+Ls2bN4//33GY6IiIgoLBG/rfbpp58iIyMjGIwAYNy4cbBYLNi+fXu1+wkhoKoq4uLiQpbb7XYIIeqtvURERNS8RTwcZWVlVeodslgsSElJQVZWVrX7mc1m3HLLLXj11VexZ88elJWVYefOnVi/fj3uvPPO+m42ERERNVMRv63mcDhgt9srLbfb7SguLr7kvr///e/x+OOP48c//nFw2fTp0/Hggw/WqU2KEvHM2ORV1JC1rBvW0TispXFYS2OwjsaRJGOPF/FwVB0hBKTLXO3SpUuxbds2PPnkk7jiiitw4MABLF++HHa7HfPnzw/73HZ7VNj7UijW0hiso3FYS+OwlsZgHRufiIcju90Oh8NRaXlJSUm1g7EB4LvvvsNLL72EZ555BmPGjAEADB48GJIk4amnnsKdd96JVq1ahdUmh8MFVdXC2pcCFEWG3R7FWtYR62gc1tI4rKUxWEfjxMdHQZaN64GLeDhKTU2tNLbI6/UiOzsbt956a7X7HT16FADQs2fPkOU9e/aE3+/HmTNnwg5HqqrB7+cH1QispTFYR+OwlsZhLY3BOtad0c9hRfxG58iRI7Fz504UFhYGl23duhVerxfXXXddtft17NgRAHDgwIGQ5fv37wcAdOrUqR5aS0RERM1dxHuOpk6ditdeew1z5szBnDlzkJ+fjz//+c+YNGlSyG21RYsWYcOGDTh48CAAoE+fPujXrx8ef/xx5OXl4YorrsC+ffvwzDPPYMKECSFTAxARERHVVMTDkd1ux5o1a7B48WLMmzcPNpsNmZmZWLBgQch2mqZBVdXg14qi4LnnnsPf/vY3PP/888jLy0P79u0xbdo0zJ49u6Evg4iIiJoJSXDGxEoKC8t4/7eOTCYZiYkxrGUdsY7GYS2Nw1oag3U0TlJSjKFTIkR8zBERERFRY8JwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTDcERERESkw3BEREREpNMowtHx48cxc+ZMDBgwABkZGVi8eDHcbneN9i0qKsLvf/97DB8+HH379sW4ceOwbt26em4xERERNVemSDfA4XBgxowZ6NChA5YvX46CggIsWbIERUVFWLp06SX3LSsrw/Tp02G1WrFo0SK0atUKJ0+ehM/na6DWExERUXMT8XC0bt06OBwObNiwAUlJSQAARVGwYMECPPDAA0hNTa1231WrVsHtduOf//wnbDYbAGDIkCEN0m4iIiJqnsK+rZabm2tIAz799FNkZGQEgxEAjBs3DhaLBdu3b7/kvm+99RZuu+22YDAiIiIiqquwe45Gjx6NsWPH4s4778RVV10VdgOysrJw6623hiyzWCxISUlBVlZWtfudOnUKeXl5sNvtuP/++/HFF18gJiYGEyZMwK9+9as6BSZFaRRDsZq0ihqylnXDOhqHtTQOa2kM1tE4kmTs8cIOR7Nnz8b69evxwQcfoEePHpg2bRomTZpU61DicDhgt9srLbfb7SguLq52v7y8PADAU089hRtvvBHPP/88jh49ir/+9a/w+XxYvHhx7S4o5NxRYe9LoVhLY7COxmEtjcNaGoN1bHzCDkcPPvggHnjgAWzZsgVr167Fb3/7WyxduhS33HILfvKTnyAlJaVODRNCQLpEFNQ0DQCQmpqKJUuWAAAyMjLg9/vx1FNP4ec//zlat24d1rkdDhdUVQtrXwpQFBl2exRrWUeso3FYS+OwlsZgHY0THx8FWTauB65OA7IVRcGECRMwYcIEHD58GGvXrsW6deuwZs0ajBgxAtOmTcOIESMueQy73Q6Hw1FpeUlJySUHYyckJAAAhg4dGrJ86NCh0DQNWVlZYYcjVdXg9/ODagTW0hiso3FYS+OwlsZgHetOCGOPZ1jMSktLw8iRI3HllVdC0zTs2LEDs2bNwi233ILjx49Xu19qamqlsUVerxfZ2dmXDEedO3eG2WyutFyUV8jIBElEREQtR50TREFBAVatWoUxY8Zg/vz5UBQFy5YtwzfffIOVK1eirKwMjz32WLX7jxw5Ejt37kRhYWFw2datW+H1enHddddVu5/FYsG1116LHTt2hCzfsWMHTCYTunfvXtdLIyIiohYo7Ntqe/bswdq1a7F582YIITBhwgTcdddd6N27d3Cb66+/HoqiYO7cudUeZ+rUqXjttdcwZ84czJkzB/n5+fjzn/+MSZMmhfQcLVq0CBs2bMDBgweDy+bOnYuf/vSnePTRRzF58mQcPXoUK1aswJ133hkyNQARERFRTYUdju644w4kJydj1qxZ+MlPfoJWrVpVuV2nTp0wcODAao9jt9uxZs0aLF68GPPmzYPNZkNmZiYWLFgQsp2maVBVNWRZv379sGrVKvzf//0fZs+ejYSEBEybNg0///nPw70sIiIiauEkIcIbxvTuu+9iwoQJVY77aeoKC8s4OK6OTCYZiYkxrGUdsY7GYS2Nw1oag3U0TlJSjKHzRYXdczRlyhTDGkFERETUWIQds1avXo0nn3yyynVPPvkkXnzxxbAbRURERBQpYYejDRs24Morr6xyXXp6OjZs2BDuoYmIiIgiJuxwdPbsWXTt2rXKdSkpKTh9+nS4hyYiIiKKmLDDkclkQkFBQZXr8vPzL/nqDyIiIqLGKuxw1KdPH6xfv77KdevXr0efPn3CbhQRERFRpIT9tNo999yD+++/H9OnT8dPfvITtG3bFjk5OXjjjTfwv//9D6tXrzaynUREREQNIuxwNHLkSDzxxBP4y1/+gocffhiSJEEIgbi4ODz55JOXfeEsERERUWMUdjgCgNtvvx0TJ07Erl27UFBQgKSkJAwcOBDR0dFGtY+IiIioQdUpHAFAdHQ0rr32WiPaQkRERBRxdQpHQgjs27cPZ86cgcfjqbT+pptuqsvhiYiIiBpc2OHo+PHjeOCBB3Dy5ElU9Xo2SZIYjoiIiKjJCTscPfHEE/B6vVi2bBnS0tJgsViMbBcRERFRRIQdjvbu3Ysnn3wSN954o5HtISIiIoqosCeBjI6ORmxsrJFtISIiIoq4sMPRLbfcgo0bNxrZFiIiIqKIC/u2Wo8ePbBp0ybMnj0b119/PRISEiptM3bs2Lq0jYiIiKjBhR2OfvnLXwIATp8+jW3btlVaL0kSDh06FHbDiIiIiCIh7HD0yiuvGNkOIiIiokYh7HB0zTXXGNkOIiIiokahzq8PKSkpwe7du1FYWIjrrrsO8fHxRrSLiIiIKCLqFI5WrlyJ559/Hm63G5Ik4V//+hfi4+MxY8YMXHvttZg1a5ZR7SQiIiJqEGE/yr927VqsXLkSt912G1atWhXyCpHRo0dXOUibiIiIqLELu+do7dq1+NnPfoZHH30UqqqGrOvSpQtOnjxZ58YRERERNbSwe45OnTqFESNGVLkuJiYGDocj7EYRERERRUrY4SguLg55eXlVrjtz5gxatWoVdqOIiIiIIiXscJSRkYEXXngBTqczuEySJPj9frzxxhsYPny4IQ0kIiIiakhhjzmaP38+brvtNkycOBE/+tGPIEkSXnvtNRw6dAhnz57F008/bWAziYiIiBpG2D1HXbp0wRtvvIFu3brhjTfegBAC7777LhITE/H666+jQ4cORraTiIiIqEHUaZ6j7t2748UXX4TX60VhYSHi4+Nhs9mMahsRERFRg6vzDNkAYLFY0LZtWyMORURERBRRYYejv//975dcL0kS5s6dG+7hiYiIiCKC4YiIiIhIJ+xwdPjw4UrLioqK8OGHH2LNmjVYvXp1nRpGREREFAlhP61WlYSEBNx2222YNGkSFi9ebOShiYiIiBqEoeGoQt++fbFjx476ODQRERFRvaqXcHTkyBFER0fXx6GJiIiI6lXYY442bNhQaZnX68WRI0fw1ltvYfLkyXVpFxEREVFEhB2OFi5cWOVyq9WKyZMn49FHHw27UURERESREnY4+uijjyots1qtSE5OrlODiIiIiCIp7HDUsWNHI9tBRERE1CjUy4BsIiIioqYq7J6j9PR0SJJUo20lScLBgwfDPRURERFRgwk7HM2dOxfvvPMOysrKcP311yM5ORm5ubn45JNPEBMTg1tuucXIdhIRERE1iLDDUUxMDJKTk/Hee+8hJiYmuLy0tBR33303bDYb7r33XkMaSURERNRQwh5z9Prrr+Pee+8NCUYAEBsbi3vvvRevv/56nRtHRERE1NDCDkc5OTlQFKXKdYqiIC8vL+xGEREREUVK2OEoNTUV//jHP+Dz+UKWe71evPzyy+jWrVudG0dERETU0MIec/SLX/wCc+fOxY9+9CPccMMNaN26NXJzc7F161bk5eVh5cqVRraTiOqRJEkQQkS6GUREjULY4WjUqFF44YUXsGzZMrz++uvQNA2SJKFfv35YsmQJhg0bZmQ7iageqZoGuYZTcxARNXdhhyMAyMjIQEZGBlwuFxwOB+x2O6KiooxqGxE1ECEAt1+FzVz1OEIiopbEkBmyKyaDNJvNRhyOiCLA7fGDnUdERHUMRzt37sQdd9yBQYMGYfTo0Thy5AgA4A9/+AP+85//GNJAImoYqibg8WmRbgYRUcSFHY527NiBmTNnwuPx4J577oGmXfymmpiYiLffftuQBhJRw9A0AZfHX+PXAhERNVdhh6Ply5dj5MiR2LBhA37xi1+ErEtPT8fhw4fr2jYiakACgNenwqey94iIWraww9GhQ4cwdepUAKj0L82kpCTk5+fXrWVE1OBUTcDj5dgjImrZwg5HiqJUmgCyQn5+fqXXihBR0+D2quCUR0TUkoUdjvr27Yt///vfVa7bsmULBgwYEO6hiSiCVFXA41cj3QwioogJOxzNmjULW7duxdy5c/Hxxx9DkiTs2bMHTzzxBLZs2YJ7773XyHYSUQPRhIDbzYHZRNRySaIO7wx499138ac//QnFxcXBZXa7Hb/5zW8wefJkQxoYCYWFZfD7OSi1LkwmGYmJMaxlHTVUHVVNIN/hhqYFvh0oioSkOBsUufkEJH4mjcNaGoN1NE5SUgwUxZCpGwGEOUO2qqrIzs7G6NGjMW7cOOzatQt5eXlITEzEoEGDEB0dbVgDiajhqWpgYHZMlJnjj4ioxQkrZgkhMHHiROzatQs2mw0ZGRmYNGkShg8fHlYwOn78OGbOnIkBAwYgIyMDixcvhtvtrtUxtm7dirS0NGRmZtb6/ERUmcuj8mW0RNQihdVzZDKZkJycbMg3TofDgRkzZqBDhw5Yvnw5CgoKsGTJEhQVFWHp0qU1Oobb7caSJUuQnJxc5/YQUYBf0+D2aXzfGhG1OGG/eHbixInYsGEDRo0aVacGrFu3Dg6HAxs2bEBSUhKAwDQBCxYswAMPPIDU1NTLHmPVqlXo0KEDOnXqhP3799epPUQUIATgcvsQZVF4a42IWpSww1F6ejref/993HXXXRg7dixat25d6emWsWPHXvY4n376KTIyMoLBCADGjRuHRYsWYfv27ZcNR9nZ2Xj55Zexbt06/OMf/wjrWoioaj5VwOvXYDZwoCMRUWMXdjj61a9+BQDIycnB119/XWm9JEk4dOjQZY+TlZWFW2+9NWSZxWJBSkoKsrKyLrv/H//4R0yZMgXp6ek1bPnlGTnivaWqqCFrWTcNVUdJ06AoUpUzY3v9GqJtpibfe8TPpHFYS2OwjsYxeuaRWoWjp556CnfddRfatWuHV155BUDgyTVFCX9MgsPhgN1ur7TcbreHTBFQlY8//hi7du3C5s2bwz5/Vez2KEOP15Kxlsao7zp6fCq8QqoyAMmShJjYKJhNzeMbOD+TxmEtjcE6Nj61Ckcvv/wybrzxRrRr1w7XXHMNVFVFnz598K9//Qu9e/c2tGFCiEtOQufxePCnP/0J8+bNC7klZwSHwwWVL9+sE0WRYbdHsZZ11FB1VDUNDt08R3qSBAjVjyhL2B3NjQI/k8ZhLY3BOhonPj4KshyheY6qejqtrk+s2e12OByOSstLSkouOd5ozZo1kGUZEydODO7v8/mgaRocDgdsNhssFktYbVJVjRNyGYS1NEZ911HVBFRVVBmOAKDU6YNFUZrFo/38TBqHtTQG61h3Rn9rivg/BVNTUyuNLfJ6vcjOzq40Fknv2LFjOHnyJDIyMiqtGzx4MH7/+9/jJz/5ieHtJWqJ/H4NflVrVjNmExFVJ+LhaOTIkXj22WdRWFiIxMREAIEJHb1eL6677rpq97vvvvtw8803hyxbvXo1jh8/jiVLlqBr16712WyiFkXVOGM2EbUctQ5Hx44dCw7AVlU1uKwqNRmHNHXqVLz22muYM2cO5syZg/z8fPz5z3/GpEmTQm6rLVq0CBs2bMDBgwcBBHqcfnjb7Z133kFOTg6GDBlS28siostw+1TERJkj3QwionpX63D02GOPVVr26KOPhnxdMZi6Jo/y2+12rFmzBosXL8a8efNgs9mQmZmJBQsWhGynaVowjBFRw/NzziMiaiEkUYsRlu+8806tDv7D215NBd+QXHd827QxGqqOqiaQX83TanoxNhPsMZYmeWuNn0njsJbGYB2Nk5QUY+h8UbXqOWqqYYeIjOHxadCEgAQOzCai5ov940RUY6qmwe3l7W0iat4YjoioxoQA3B6/4VP1ExE1JgxHRFQrFS+jJSJqrhiOiKhWNE3A41XZe0REzRbDERHVmserNskn1oiIaoLhiIhqTdV4a42Imi+GIyKqNU0EXifCW2tE1BwxHBFRWDw+DbWYQ5aIqMlgOCKisKiaBg9vrRFRM8RwRERhCcx5pELivTUiamYYjogobF6fyndCEVGzU6t3qxFR86IJgeycEri9KjxeFa0ToyDXoidI1QScHh/sMVaOP6JGo+JzXer0ITbajJS2cbX6XBMxHBG1UIdOFGDTzpM4X+BEK7sNpS4vom1mjOrfAamdEmp8HLdPRbSqQZH5w4ciT/+5VlUBRZHQLikaE4d2Qc+uSZFuHjURvK1G1AIdOlGANVuO4HRuKaxmBTHRZljMCs4XuPDO58eRdbqoxsdSVQEn37dGjcAPP9f2WAusZgWnc8uwZssRHDpREOkmUhPBcETUwmhCYNPOk3B7/UiItcJiViBLEswmBfZoMzw+Ddv2nIVWi9tkbq8fKoceUQRV97m2mBUkxFrg9qrYtPNkrT7X1HIxHBG1MNk5JThf4ESMzVzpSTNJkhBtVZBX5MK5vLIaH1PTBDw+v9FNJaqxy32uY2wmnC9wIjunJEItpKaE4YiohSl1+qCqAiZT1X/9FUWGqgFOd83DTuCxfj9kjjuiCLnc59pkkqGqAqVOXwO3jJoihiOiFiY22gxFkap9BF9VNSgyEG2r3fMaPlXA5+O9NYqMy32u/X4NiiIhNtrcwC2jpojhiKiFSWkbh3ZJ0Shz+ys9fi+EgNOjIjkhCu2TY2p1XE0TcPs4MJsi43Kf6zK3H+2SopHSNi5CLaSmhOGIqIWRJQkTh3aBzaKgqNQLr0+FJgR8fhUOpw9Ws4xR/TuENS+M26uCw10pEqr7XHt9KopKvbBZFEwc2oXzHVGNMBwRtUA9uyZhxrg0dGodA49PRZnTB69PRbukKNw8/IpazXOkp6oCXt5aowj54efaUeqFx6eiU+sYzBiXxnmOqMY4CSRRC9WzaxLSuiTWaYbsH9KEgNPtgzXOCj4xTZGg/1xzhmwKF8MRUQsmSxK6trND1QTyHW5oWt0TjdevweNTYTEpBrSQqPYqPtdE4eJtNSIylKYJlLl8AP+hTkRNFMMRERnO5xdwe9VIN4OIKCwMR0RkuIqxR0RETRHDERHVC59fg4svpCWiJojhiIjqhRCBV5DwRZ9E1NQwHBFRvfGpGspc7D0ioqaF4YiI6pXL64dfZe8RETUdDEdEVK9UVaDU5YPE7iMiaiIYjoio3nl8Kjw+PtpPRE0DwxER1TtNEyhzc2JIImoaGI6IqEF4fSonhiSiJoHhiIgaRODRfk4MSUSNH8MRETUYn1+Dm2OPiKiRYzgiogbD3iMiagoYjoioQbH3iIgaO4YjImpQ7D0iosaO4YiIGpzPr8Ht9Ue6GUREVWI4IqIGJwRQ5mE4IqLGyRTpBhBRZG3bfQYnzpVgcM82aGW3Ndh5/X4NLq8f0VYTBF+9RkSNCMMRUQuWnVOCVzYfAQD87/AF3DOxJzokxzTIuQNjj/yIsijg1NlE1JjwthpRCxYbZYZJCQQTp8ePFzcdxJnc0gY7v1/V4PSo4DtpiagxYTgiasGS7DbMmtQbihxIJy6Pihc3HcKpCw0TkCp6j1StQU5HRFQjDEdELdzV6W0wa3JvyOUBye1V8dKmQzh5vqRBzu9XNZS5fZDYfUREjQTDERFhUI/WuPOGHsEeJI9PxcvvH0LW2eIGOb/b64ffz+4jImocGI6ICADQ+4okTBvbIzgGyevXsOaDw/juVFG9n1tVBUpdXo49IqJGgeGIiILSUhIxfVwazErgW4NfFXh1yxHsP15Q7+f2+DS4vHytCBFFHsMREYW4slMCfjYhHRZz4NuDqgm88eF3+Pa73Ho9ryYEytw+AJz0iIgii+GIiCq5or0dMyf2gs2iAAg8VfavbVn4cv/5ej2v36eh1OXn4GwiiiiGIyKqUuc2sbhvUi/ERpmDyzZ+eQIffXMaop6mtBYIzLfk8fk5/oiIIobhiIiq1b5VDGZN7oWEWEtw2UffnMbGL09Cq6eApGkCJWU++FTeXiOiyGA4IqJLSo6Pwv2Te6N1wsX3ru04cB7//OQo/PU0e6NP1VBS5gHHHxFRJDAcEdFlxcdaMWtyb3RqffG9a3uO5uPVLUfg8dXPE2Yen4YSp4+314iowTEcEVGNxNjMmDmxF7p3jA8u+/50MV547yBKnN56OafLo8Ll4fgjImpYDEdEVGNWi4K7bkxD326tgsvO5JVh1b8PIL/Ybfj5NCFQ4vLB7+ftNSJqOAxHRFQrJkXGHWO6I6NPu+CyAocHz767H6cuGP8+NlUVcDg5/oiIGg7DERHVmixJyMzoghuHpASXOd1+vPDeIRw8Yfxs2h6fBgfHHxFRA2E4IqKwSJKEkf074MejuwdfWOtTNaz9z3f4Yt85w8/n8vhR5ub4IyKqf40iHB0/fhwzZ87EgAEDkJGRgcWLF8PtvvT4hdLSUqxYsQK33347rr76agwdOhQzZ87EgQMHGqjVRAQAA65Mxs8mpF+cTRvAph0n8d4XJ6Bpxt0KEwIoc/ng9dfP9AFERBUiHo4cDgdmzJiBsrIyLF++HL/61a/w3nvv4Te/+c0l9zt79izefPNNDBs2DMuWLcOSJUugaRqmTp3KgETUwFI7xGPW5N4hk0XuOHAer/7nCDwGvkxW1QQcZV4YmLmIiCoxRboB69atg8PhwIYNG5CUlAQAUBQFCxYswAMPPIDU1NQq9+vUqRO2bt2KqKio4LJhw4ZhzJgxeO2117BkyZIGaT8RBbRLisbsm/rg1S1HcCa3DABwJLsIq/59ANPHpSExzmrIeXx+DY4yT3kQ4z02IjJexHuOPv30U2RkZASDEQCMGzcOFosF27dvr3a/6OjokGAEAFarFampqbhw4UK9tZeIqmePtuC+zF7o1TUxuOx8gRPPbNiPk+eNe5LN41VRWOKB16/xJbVEZLiI9xxlZWXh1ltvDVlmsViQkpKCrKysWh3L6XTi0KFDmDJlSp3apCgRz4xNXkUNWcu6aag6SpoGRZEMGewcpZgw/cY0bPkqG9t2nQUQGCv0wsaDuPW6brgqvU3dTwLArwmUuLxQNRPioi2X3Z6fSeOwlsZgHY1j9L+RIh6OHA4H7HZ7peV2ux3FxcW1OtbTTz8Nl8uFadOm1alNdnvU5TeiGmEtjVHfdfT4VHiFBCPfJTt1XE+ktI/H2s2HoWoCqiaw/pMsFJT6cPPoVCiyMT8QJAAwKUiItdaoF4mfSeOwlsZgHRufiIej6gghatVd/t5772HNmjX43e9+hy5dutTp3A6HC2o9vVCzpVAUGXZ7FGtZRw1VR1XT4HC4DX26DAB6d0nAfZN74dXNR1Dm9gMAPvxvNk6eK8ZPb7gS0TazIecpdgClJRbERVsgqkl4/Ewah7U0ButonPj4KMgG/YMLaAThyG63w+FwVFpeUlJS7WDsH/riiy/w2GOPYebMmbjzzjvr3CZV1eDn48KGYC2NUd91VDUBVRWGhyMASGkThzk398Fr//kO5/KdAALvZFvxr324c2wPtG8Vc5kj1Iyj1AtNE4ixmasNSAA/k0ZiLY3BOtadkb3eQCMYkJ2amlppbJHX60V2dnaNwtHevXvx4IMP4sYbb8QjjzxSX80kojpIjLPh/sm90afbxQcvCko8eO7dA9hzNM+Qc2hCoNTlg8vLiSKJqG4iHo5GjhyJnTt3orCwMLhs69at8Hq9uO666y65b1ZWFu677z4MGjQIS5Ys4VMrRI2YxazgJ2OuxNjBnYMP4Pv8Gt78+Cg2fnkCqlb3fzlrmkCJ0wuXV+X3AyIKW8TD0dSpUxEXF4c5c+bgs88+w4YNG/Dkk09i0qRJIT1HixYtQq9evYJf5+fnY+bMmTCbzbj33ntx4MAB7N69G7t378bBgwcjcSlEdBmSJGHUwI6YMf7ijNoA8OX+83jhvUNwlHnrfA5VDUwUWeKs+7GIqGVqFGOO1qxZg8WLF2PevHmw2WzIzMzEggULQrbTNA2qenGm3aNHj+LcucD7m372s5+FbNuxY0d8/PHH9d52IgpPj84JePCWvli79eI4pJM5JVjx9j7ccX13dO8YX6fja5pAmdsHv6rBHmOFzE4kIqoFSVxq5GILVVhYxsFxdWQyyUhMjGEt66ih6qhqAvn18LTa5Xj9Kt797Dh2fX9x3JEkAdcP6oTRAztCNiDVmE0y4qItiIkyISGBn0kj8O+3MVhH4yQlxRg6X1TEe46IqOWymBTccl03xMdY8Nnec1A1ASGAj745jRPnHfjx6O6Ii7ZAEwLn8srgdPsRbTOhfXIM5BqOKfL5NRSXeqBBwG6vv/CnCYHsnBKUOn2IjTYjpW1cjdvYlPg1DV/ty4HToyHaKuOqtNYwGfgINVFjwHBERBGTdboI2/acRV6RCxazAo/XH3ypbNYZB1a8tQ/D+7TD92eLkVfkgqoBigwkJ0RhVP8OSO2UUKPzqJpAqdOLAocbUj10lh86UYBNO0/ifIETqiqgKBLaJUVj4tAu6Nk16fIHaCI2f3USm3achMvjh0BgAs7XrCZMzOiCG4fUbX45osaEcZ+IIiLrdBHe+fw4zhc4YTEriI+1INFuhaK7lVbq8mHzf08hO6cEZpOM2GgzLGYF5wtceOfz48g6XVTj8wkRmAm80OGBX9UMe9z/0IkCrNlyBKdzS2E1K7DHWmA1KzidW4Y1W47g0IkCY04UYZu/Oom3th9DmdsPWZZgUiTIsoQytx9vbT+GzV+djHQTiQzDcEREDU4TAtv2nIXHp8IebYHZpECWJFjNJrROsMFilqHPLj6/QHGZD5omYDYpsEeb4fFp2LbnLLRa9gR5/SqKSise96/7dWzaeRJurx8JsVZYzIHrsJgVJMRa4Paq2LTzZK3b2Nj4NQ2bdpyEqgmYy0ORLMmQZQlmRYKqCWzacRJ+A6ZjIGoMGI6IqMGdyytDXpEL0VZTpfmIZFlGfIwFVqsSEl58fg25hS443T4AQLRVQV6RC+fyymp9fr+qwVHmRVGpF5oI/6WV2TklOF/gRIzNXOk6JElCjM2E8wVOZOeUhHeCRuLrgzlwefwwyVKV12mSJbg8fnx9MCdCLSQyFsMRETU4p9sfGD9UzdMliiIDGmAxyYiyXpwPSQAoKvWisMQDSZKgaoFjhUPTBFwePwodbrg84c2qXer0QVUFTKaqr8NkkqGqAqVOX1htbCwKHG4IAKiuRlLgz6bA4W64RhHVI4YjImpw0TYTFBnVvmxTVTUoJgkmk4LYKDMS46wh4cXtVZFb5IKmaYi21e25Ep+qobgsELgCY5FqnpJio81QFKnax7D9fg2KIiE22pgX7EZKkt0WyEXV3R0UgdyUZLc1XKOI6hHDERE1uPbJMUhOiILTo1Z6SawQAk6PiraJUWiXFNjGZlHQJiEKFrOs2w5welTsPJgDtze83iP9sdxeFYUlHhSV1jwkpbSNQ7ukaJS5/VVeR5nbj3ZJ0UhpG1en9kXaNb3aIspqgl8TVV6nXxOIsppwTa+2EWohkbEYjoiowcmShFH9O8BqluFw+uDzq9CEgM+vwuH0wWqWMXpAx5BtNCGQGGdFzA96ir45kovl/9qL72vx5Fp11IpbbSUeFJd5oGrikiFJliRMHNoFNouColIvvL7AdXh9gUHfNouCiUO7NPn5jkyyjIkZXaDIEnyqgKYJaEKDpgn4VAFFljAxowvnO6JmgzNkV4GzldYdZ341RnOfIVs/z1F1cxhVtY09xgKXR0VecegYl8HpbTB+aApslsq32hRFQnx8NIqLA3MR1YSiSLCZFURbzTCZ5Eq9JhVa8jxHUZznKGz8Pmkco2fIZjiqAj+odce/9MZo7uEIQI1mv65qGyEEtu8+i4+/ORPyqLw9xoKbRlyB9JTEkGOEE46C+8oSrGYlMMZIllDVd82WNEP2N0dyOUO2Afh90jgMRw2AH9S64196Y7SEcFRX5/LL8K9tWcEX2Fbol9oKmcO6IjYqMBi6LuGogqJIiLGaygeBN7/gU1P8+20M1tE4Rocjxn0iatLat4rBnJv74IarO4fMrr03Kx/L1u/BN0cuVHs7rLZUVaDE5UN+iQduX90nkSSixonhiIiaPEWWMXpQRzx4S190bhMbXO7yBF5t8fzGg8gpdF7iCDUnBODzaXCUelFQ4oHXz5BE1NwwHBFRs9E2KRr3T+6NScO6hjz2f+JcCZ5evxfvbDsKr0815FyaEPB4A0+lFZT3JAHhz7ZNRI0HwxERNSuyLCGjTzs8dHt/9Op6cVC2pgls2XkS/7duN/YfyzfsVpumBUJScakH+cVulDh9wRfbMigRNU0MR0TULMXHWjFtbBqmj+2BhFhLcHlRqRevf/g9Xn7/MC4Uugw7nxCB2bZLXT4UlHiQ73CjzOW77FxJRNT4MBwRUbPWs2sSfnF7f4wa2CFkwPbRM8VY/q892PjlCbg8dZth+4c0TcDr0+Bw+lBQ4kZxWWDWbVlmSCJqChiOiKjZs5gVjB/aBb+dOQRXdo4PLtcE8OX+8/i/dbux88B5qPUwlYGqCjjd/mBvUmBskuAtN6JGjOGIiFqMdq1iMHNiT9x5Qw8kxlmDy50eP/79xQks/9deHM4uNGw8kp5+bFJesYdjk4gasbq9zpqIqImRJAm9r0hCj84J+GLfOWzbdQbe8gn4cotceGXzEXTrYMeNQ1LQqXXsZY5We0IAflVDqUuD0+OHSZFgs5hgNSswKVJwGyKKHIYjImqRzCYZowZ2xKC01tj631P49kguKjLJsbMOPPPOfvTt1gpjB3dGq3hbvbRB0wS8moDX54UsSzArEqzlQclc/i43BiWihsdwREQtmj3agluvS8WwPu3wwc5sHD1THFy371g+DhzPx9XpbTB6UCfEx1gucaS60TQBjybg8XmhyBJMFUHJpMBkYo8SUUNiOCIiQuA1JHdPSMf3p4ux5evs4LvaNAF8fegCvv0uF0N6tcV1AzoG39dWX1RNQC0PSsEeJbMCq8UEkyLXy5goIrqI4YiIqJwkSejROQHdO8Vjz9E8fPi/0ygs8QAA/KrAF/vO4+tDF5DRuy2G9+tQ7yEJ0PcoaVDcfphMMqIsCqwWBbLEoERUHxiOiIh+QJYkDLyyNfp2a4X/HbmAT749gxKnDwDg82v4dM857DiQg6G92mJ4v/aIi66/2216qiagelV4fSpklwSLSYHFLMOsyDCZZEjgrTciIzAcERFVw6TIGNqrHa7q0QZfHczBtt1n4HQHJoz0+TV8tvccdhw4j8E922JEv/ZIiLVe5ojGECIwf5JL9cPtASRZgiIHbr1ZzArMJgmyJAPggG6icDAcERFdhtkkY3i/9hjcMxCSPt1zNhiS/KrAjv3n8fXBHAzonowRAzqgTUJUg7VNABCagKYJ+PwaZLcfkhwIdlaTDLNZgVmRIUkSb8ER1RDDERFRDVnNCkb274Chvdriq0M5+GzPOZS6ArfbVE3gm+9y8e13uejZNREj+nVAl3ZxDd5GTQhABVRVhcerQpb9kCUJZpMMi0kOTDgpSTDLMhSTBAkMTUQ/xHBERFRLFrOCEf06YGivdvjf4Qv4bO9ZFJV6AQR6cg6eKMTBE4Xo3CYWw/u1R6+uSSHvdWtImiagQcCvanB5AAkApMDgc1mSYDHJMJtlWBQFiiKVvySXt+OoZWM4IiIKk9kkI6NPO1zTqw32HM3Hp3vO4kKhK7j+1IVSvPHh90iItSCjTztcndYGUdbIftsV5f8R4mJokrwXw5JJCQQmk0mGIstQZIm35KjFYTgiIqojRZYxqEdrDLgyGUeyi/D53rM4fq4kuL6o1IsPdmbjw/+dxsArk5HRux3aJkVHsMWhREhYAtxeFZIUeGpPkiQoFfMsmRRIcuDWnRB8eS41XwxHREQGkSUJPbskomeXRJy+UIrP953D/mMFgXFACDzh9vWhC/j60AVc0T4OQ3q1Q6+uiTApje8d4EIAqhBAeWAKjF+SIEmAIkvwCQkelxeyBJhN5bfkOH6JmgmGIyKietCpTSymjrkSxUM82HkwB/89dAFOjz+4/vi5Ehw/V4KYKDOuTmuNweltkGSvn3e4GUXTLgYfVRMoc/ugaSLQuyRJMJlkmE0yFEUq73UKbKdpAhIAufw2nSwHettExT0+cH4malwYjoiI6lF8rBXjrknB9YM6Yc/RPOw8cB5ny19NAgBlLh+27z6L7bvPIrWjHVentUGvrkkwmxpfb1JV9LfkfBWDviVAggRICL48V7+s4padLAdClaxI5aEp0PskyxfXS+AAcWp4DEdERA3AbJJxdXobXJXWGqculGLngRzsP54Pv3rxp37WGQeyzjgQZVXQPzUZg9Jao2NyTPkTZE2HEICAqOgUqnKZitC0E3yKriJAAcHxTmYl0BtlUgI9T4En/wKhqeLYREZiOCIiakCSJCGlbRxS2sYh090F336Xh/8ezkFukTu4jcujYufBHOw8mIM2iVEYeGUyBnRPRnwDzcAdCcGn6EJCVfl4J6jBsCTJgQBlViTIigxFCswQLsuBW3kVt/Mqep6ECDkDgxTVCMMREVGERNvMGN6vPa7t2w4nc0rwv8O52HcsHz6/FtzmQqELW74+hf98fQpXdLCjf/dk9LkiKeJTAjQ0gcAtOqiBr/wqUP5F5V4n3e8VWYIslY93Kg9SklR+G688SFVsX9FBV9FT98NgFdIehqxmrWX97SIiaoQkSULXdnZ0bWfHpGFdse9YPr79Lhcnzl+cDkAAOHbWgWNnHfj358dxZacE9OveCj1TEmG1KJFrfCNQXa8TAFwcAl9NkNJ/XaF8uSxLkMv/L+km8ZRxMWRJuluAwV6r4Fgp3SF146+EEBAAVE2D16cGpkaAgKYh+GRjRWATAEyKVFU+o3rEcERE1IhYLQquTm+Dq9PbIN/hxu7v87Dr+1wUODzBbVRN4HB2IQ5nF8KkSOjROQF9urVCekoCbBZ+W7+UqoNUcE2tSCF5KnSslCxXBC9JF4oqzhT4jVw+JYLD4YKqlQ86F8EDBrdJjLVGbIb1lop/i4iIGqlWdhvGXNUJ1w/qiFMXSrH7aB72ZeWjzH2xP8SviuDrShRZQveO8eh1RRJ6dklEbJQ5gq1v/vS31ir1Wqk1O4YmAmFXVX8Qzsq/ZCSKDIYjIqJGTj+Ie2JGVxw/68DerDwcOFEAl+fiT2FVEzhyqghHThVhA4CUtnHo2TUR6V0S0Tre1uSeeiOKFIYjIqImRJEldO8Uj+6d4jFlxBXIOuPA/mP5OHiiMGSSSQHgZE4JTuaUYPNX2UiyW5Gekoi0lARc0d7eKGflJmosGI6IiJooRZbRo3MCenROwJQRAifOOXDgRAEOnShEcZk3ZNsChwdf7j+PL/efh9kkI7WDHVd2TkCPTglIslvZq0Skw3BERNQMKLKE1I7xSO0Yj0nDuuJsXhkOnSzEoZOFOKebkRsIvOPtcHYRDmcXAQAS46zo3jHQG5XawY5oG8cqUcvGcERE1MxIkoSOrWPRsXUsfnR1ZxSXenA4uwjfnSrC0TPFIfMoAUBhiQf/PXwB/z18ARKA9q2i0a1DPLp1sKNr+zg+AUctDj/xRBRxFXPNcGK9+hEfa8WQXm0xpFdb+PwaTp4vwXenA2HpQqErZFsB4Gy+E2fznfh83zlIEtChVQyuaB8ISl3axcEeY4nMhRA1EIYjIqoVKfifizMJB2cW1j14HJx5WJLK5365OAuxJJcvK9+wYiK9wCzIFyfJE+VvdNdQ/nLTig6PipeXlu8vly+Ty08aHD5T8VLU8jefVrxSQtNa7otMzSY5OKB7wtAuKC7z4ujpovL3uhWjxOUL2V4I4ExeGc7kleHzfecAAG0So3BlSiI6JEWhc+tYtOKTcNTMMBwREYBAsJBkhMwWHHw7evmb0/UzA8vBie5CX8EQDCmoqicodEF1ASX056xUxbKq9r9U2pFgMsmIs0fBJGlwu/3w+jSoQoQEJaE7hhCApjX/BBUfY8FVaW1wVVobCCFwociFY2ccyDpbjBPnSkKegKtwodAV0uMUbTWhc9tYdG4T+NWpdWyLe70JNS/89BIRFFlCQpwFEkIDz8Up6EJDwyWJi69ACFdVoaduhxQQQkCRJZhkGTaLCVFWAOXXqX+tQ8Ul+1UBn0+FT9XgVwO9T829x0mSJLRNjEbbxGhk9GkHTQhcKHTh+FkHjp9z4MT5EpT+oGcJAJweP45kF+FI+QBvAEiOt6FT61h0bB2Djq1j0L5VDKzmlv2aE2o6GI6ICABgkkPnvREV71lopiquTx929LeGzIoEiynQlSZJgbDk9fnh8qrw+zWoLaBXSZYktEuKRrukQFgSQqDA4UH2hRKcK3Dhu+zCSmOWKuQVu5FX7Mbuo3kAApkzOSEKHZKj0aFVDNonx6B9q2jE8Mk4aoQYjoiIqqEPULKE8h4nM/x+DS6vHx6vCr+mNeveJD1JktAq3oY2SVGIj49GcbETpU4fTl0oDfnlquJWnACQW+RCbpELe47mB5fbYyxolxSN9q2i0bY8iCXH2zhJJUUUwxERUS0IIaAoEuKizYiNMsHt0+D2+OFTtWZ/260qUVZTcCJKIFCfghIPTl8oxZm8MpzOLcXZvDJ4fVqV+zvKvHCUefHdqaLgMrk8hLVNjEKbxCi0SYxGm8QohiZqMAxHRERhCIQgCTazgiiLCaqmweNTy4OSaBGDuasiSRJa2W1oZbehf/dkAIExaPnFbpzNK8PZvDKcy3fibH4ZnO7KPUwV21f0MuG4/thAUpwNyQk2tE4IhKXk+CgkJ9gQF2XmE3NkGIYjIqI6EkJAliREWUyItprg9WuBwdx+LRiU6jpIvSmTJQmtE6LQOiEqGJiEEHA4fTifX4bzBU6cy3fiQmEgEFU3nksIIN/hRr7DHTL4GwAsZjkYylrF25BktyHJbkUruw32aAtkmcGJao7hiIjIQEIAZkUODubWhAafX8DrU+HxqlBbeFCqIEkS4mMsiI+xIC0lMbhc1TTkFbuRU+DChcJAYLpQ5EJ+sfuSg+C9Pg3n8p2VXpUClD+NGWtFYlzor4RYKxLirIiLMjM8UQiGIyKielAxmFuCBItJgtUsIzbKDK8/EJI8Pq3F9yhVRZHl4HQCQKvgclXTUODwILfIhbwid+D/xW7kFruqvT13cV8R7HGq+pwS7DEWxMdakBBjRXxsILTFx1oRH2OBPcaCaJspOH8XNX8MR0REDaAiA1lMCqxmBZoQ8KsCfr8Gv6rB69daxFxK4VJkOXhr7odcHj/yit0oKA9A+cVuFDg8KHC4K834XRVVEygs8aCwxAOgpJrzBwbh22MsiIu2BH5f/v+Kr2OjzIixsReqOWA4IiJqYIG5JiWYFQlmRQ7O/q1qAn5Vg88fCEtq+eSTQhPNeMapuouymoKzc/+Q16+i0BEIPgUlHhSWuFFU4kVhiRuFpd4qpx2oiqoJFJV6UVTqveR2kgTE2AJBqeJXTJQpGJxiosyIsZkQE2UOBCwm4UaJ4YiIKMIqfj7KkhTsWaoYr+RXywOTTwtOF6AJ9i7VlMWkoG1SYA6lqni8KorKPCgu9aKoNPD/4jIPHGU+FJV64CjzwuuvehqCqggBlLp8Vc4kXhVZlhBjMyHKakK0LTCgP9pmDvzfakJMlAnX9GqL1vGVe8yo/jAcERE1MvrxShW9S9HWwLqK3iW/P9DD5NeFJcHQVGtWi4K2looxTpUJIeDxqSgu86KkzAeH04sSpzfwtdOHEmfg/6VOH3xqzUNUBU0T5cepPkxt2nESS+7PQHyMpdbHp/A0inB0/PhxLF68GN988w2ioqIwceJELFiwADab7bL7vvPOO1i1ahXOnDmDLl26YO7cuRg/fnwDtJqIqOGI8nfWZeeUoMzlQ1yMBZ3bxEGWAoFJ1QQ8fhVf7jmH/BI3Ym1mXN2zDUyyHAhOCGxzLi8wv1C0zYT2yTGVBhlr4tLbqJqGPVkFcHs12Cwy+lyRBEWu+cSMlzt+TbepyzlUTcPeo3koKvUiIdaCft2Tq70GSZJgs5hgs5jQ9uJDdZXO0a5VNPx+gVKXDyWuQGA6m1uK4jIv/KqAgECZ248ylx9lbh9cbn+Nb5W6vSoOZxdgcHpbDgpvIBEPRw6HAzNmzECHDh2wfPlyFBQUYMmSJSgqKsLSpUsvue/mzZuxcOFCzJo1C9deey0+/PBDPPTQQ4iLi8Pw4cMb6AqIiOrfoRMF2LTzJM4XOKGqgVm62yVFY+LQLujZNQmvbj2MHQfOQ5JkmBQZZpOMz/efx5BebTBp2BU4eroYX+w7h9xiN7w+FRAC8bFWDOvdFl3a2wEAWaeLsH3POVwodELVAEUOvA9tVP8OSO2UgM/2nMG2XWfh9qoQCLwv7d8WBaMGdsCI/h0vew1Zp4uwbc9Z5BW5qjx+TbepyzmquoaNX56s8TXU5BxFJW58fSin0voxAzsGr0HTBLx+FbLZhPO5JSh1+uB0+wO/PH6cLyjDmdwy+FUNcdEW/POTLHy251zwz5vqlyREZDthV69ejWeeeQYff/wxkpICf+DvvfceFixYgPfffx+pqanV7jt+/Hj06NEDf/vb34LLZs6ciZKSEqxfvz7sNhUWlsFfi3vMVJnJJCMxMYa1rCPW0ThNuZaHThRgzZYjcHv9iLGZYTLJ8Ps1lLn9sFkUtEuKwt6sgmr37989CblFHqha4Aet1aJAUzV4fBqibCZMGtoVmhDYtPMkPH4VURYTZFmCz6fC6fFDkSV0bhOL3d/nwacKSBAQQpTf4hOAELhhcOdLhous00V45/Pj8PhURFtNUBQZqqrB6VFhNcu4efgVAHDZbS4VkC53jtQOdnxzJBeaABQJgWQkALX83Xnjrrn0NdTkHEN6tsFXhy7U6BoURQq+o05VRZXnaB0fhXi7FRfyncE/7xnj0hiQfiApKQaKga+WifhLaj799FNkZGQEgxEAjBs3DhaLBdu3b692v1OnTuHYsWPIzMwMWZ6ZmYm9e/eioKD6bxRERE1FRWhxe/1IiLXCYlYCA7fNChJiLShzeUOCkaT7VWHP0QIUl7pgNSvw+bVAL4VHhV/VcCqnFBt3HMcHX53A+YIyQFycsBIIPHnl9ak4cLwAVosJibEWJMXbkJwQjTZJ0ejQKhptkmJw5FQxbFYFUVYFVosCiznQg2VSZMiyhP8eyYXJFJjDKC7agihr4ImtVnYrNAF8uvcsPt13Dl6/ivhoC8ymwHWaTQrs0WZ4fBq27Tlb7bxQmhDYtucsPD4V9ir2d3v8+F95MDLJgYHQsiRBliWYZEATwLZdZ6Fq1Qfny53D41Wxbdcl1l/mGqo6h8kkh/x5u70qNu08yfmx6lnEb6tlZWXh1ltvDVlmsViQkpKCrKysavc7duwYAKBbt24hy1NTUyGEwLFjx0ICV23Ex0dxUGMdVdwWZy3rhnU0TlOtpV/V8MtpVwcCTxXjTVQt0ItzOZIkQali/h39IG5Jqvoc+skqg2ur2C4u2lz+pF3la3j4zquqvQZR3g5Usx669TE2c8h1VFy5qmn4zT1Dqt0/cA3l7Zd0O/7gODE2EyxVXEPFdfz67iGAFJiKoao2qpoIBq+q1gsAsVHm4At0ZVkKmbvJr2pY9LMh5X8WgXpIkhR8V1/FMeJjLHwJr47Rc0tFPBw5HA7Y7fZKy+12O4qLi6vdr2LdD/eNj48PWR8OuRaDC+nSWEtjsI7GaWq19KsXQ0Ndx+JWvb8UDD5ybc6hC2QVv9M0UeWtDb9aHhokCVI1w5ArXg1Ssbuk+2/gd4F2SrIEcxXhxeMr30qqKrYgcJ9E0wWwShsFQqIkSTCZ5ErhCSh/h15F8JGCu+nWSwA0KLJ8yVrLshQSbGTl4saaCIwn+2G4uvilVB4UJUNvI1GoiIej6gghavSG5R9uc7l/fRARNSVWi4I2lqofM28qGuIarGYF1ipmzw5bFT9CrBYTrJb6/bFZ8WQcRVbEY6fdbofD4ai0vKSkpMoepQrV9RBVHOtS+xIRERFVJ+LhKDU1tdLYIq/Xi+zs7Es+qVYx1qhi7FGFrKwsSJJUaSwSERERUU1EPByNHDkSO3fuRGFhYXDZ1q1b4fV6cd1111W7X+fOndGtWze8//77Ics3btyIfv36hT0Ym4iIiFq2iIejqVOnIi4uDnPmzMFnn32GDRs24Mknn8SkSZNCeo4WLVqEXr16hew7f/58fPDBB1i2bBm++uor/OlPf8IXX3yB+fPnN/RlEBERUTMR8VFfdrsda9asweLFizFv3jzYbDZkZmZiwYIFIdtpmgZVVUOWjR8/Hm63G8899xxefPFFdOnSBcuWLePs2ERERBS2iM+QTURERNSYRPy2GhEREVFjwnBEREREpMNwRERERKTDcERERESkw3BEREREpMNwRERERKTT4sLRZ599hmnTpmHo0KHo06cPxowZgyVLlqCkpCRku+3bt+Omm25C3759ccMNN2Dt2rURanHTUFZWhpEjRyItLQ379u0LWcdaXtrbb7+NtLS0Sr+WLl0ash3rWDP//Oc/MXnyZPTt2xcZGRmYPXt2yHrW8fKmT59e5WcyLS0NmzZtCm7HWtbMhx9+iNtvvx2DBg3CsGHD8OCDD1Z69RXAel7OJ598gptvvhl9+vTBddddh+XLl1ea/xAwpo4RnwSyoRUXF2PgwIGYMWMG7HY7vv/+e6xYsQLff/89XnrpJQDArl27MGfOHEyZMgULFy7Et99+i8WLF8NiseD222+P8BU0Ts8880yVH1LWsuZeeOEFxMXFBb9u27Zt8PesY82sWLEC//jHPzB79mz0798fxcXF+Oyzz4LrWceaefzxx1FaWhqybM2aNfjPf/6DjIwMAKxlTX355Zd48MEHMXnyZPziF7+Aw+HA3//+d9x9993YtGkTYmNjAbCel7N7927MmTMHEyZMwMMPP4ysrCwsW7YMLpcLv/rVr4LbGVZHQeLNN98UPXr0EOfPnxdCCDFz5kxx2223hWzzm9/8Rlx77bVCVdVINLFRO3r0qBgwYIB44403RI8ePcTevXuD61jLy3vrrbdEjx49RH5+frXbsI6Xd/ToUdGzZ0/x2WefVbsN6xi+66+/Xtx3333Br1nLmlm0aJEYPXq00DQtuGzPnj2iR48eYtu2bcFlrOel3XPPPeLmm28OWfbCCy+I3r17i9zc3OAyo+rY4m6rVSUhIQEA4Pf74fV6sXPnTkycODFkm0mTJiE3NxcHDx6MQAsbtz/+8Y+YOnUqrrjiipDlrKUxWMeaefvtt9G5c+dqXx/EOobv22+/xenTpzFp0iQArGVt+P1+xMTEQJKk4DJ9DzHAetbEoUOHKv3dHjFiBHw+Hz7//HMAxtaxxYYjVVXh8Xhw4MABrFy5EqNHj0bHjh2RnZ0Nn8+Hbt26hWzfvXt3AEBWVlYkmttobd68GYcPH8bcuXMrrWMtayczMxM9e/bEmDFjsGrVquBtStaxZvbs2YMePXpg5cqVyMjIQJ8+fTBt2jQcOnQIAOtYFxs3bkRUVBTGjBkDgLWsjdtuuw3Hjh3Dq6++CofDgdOnT+Mvf/kLUlNTg7coWc/L83g8MJvNIcssFguAi/Uxso4tbsxRhdGjRyMnJwdAIH3+9a9/BRAYkwQEXoirV/F1xXoCXC4X/vznP+Phhx8O3jfXYy1rpnXr1pg3bx769+8PSZLw8ccf4+mnn0ZOTg5+97vfsY41lJubiwMHDuD777/HH/7wB5jN5uDYjv/85z+sY5j8fj82b96MMWPGIDo6GgD/btfG4MGD8fe//x2//OUvsXjxYgCBH9YvvfRS8Ic763l5Xbt2xd69e0OW7d69G8DF+hhZxxYbjlavXg2n04mjR4/imWeewezZs/Hyyy8H1+u7QPWqW94SPfvss2jVqhVuueWWS27HWl7aiBEjMGLEiODXw4cPh9VqxZo1a0KetGIdL00IAafTiRUrVuDKK68EAPTu3RtjxozBm2++iUGDBgFgHWvriy++QH5+PjIzMyutYy0v79tvv8UjjzyCW2+9Fddffz1KS0vx3HPP4b777sMbb7wR8g9L1rN6d955Jx577DGsWbMGU6ZMwdGjR/H0009DUZRK9TGiji32tlp6ejoGDRqEH//4x/j73/+Or776Clu3bkV8fDyAygnT4XAAqJxIW6ozZ87gpZdewvz581FaWgqHwwGn0wkAcDqdKCsrYy3rYPz48VBVFYcOHWIdayg+Ph7JycnBYAQAbdq0Qbdu3XD06FHWMUwbN25EQkJCyHgP1rLmFi9ejKFDh+LXv/41MjIycMMNN+D555/H8ePH8c9//hMA61kTN998M+6++2489dRTGDJkCH72s59h6tSpiI+PR+vWrQEYW8cWG470evbsCUVRkJ2djZSUFJjN5kpzUBw9ehQAkJqaGokmNjqnT5+Gz+fDrFmzMHjwYAwePDjYy3HXXXfh7rvvZi0NwjrWTHV1EEJAlmXWMQxutxsfffQRbrzxxpDxHqxlzWVlZSE9PT1kWVJSEtq0aYPs7GwArGdNSJKEhQsXYufOnXj33Xfx5Zdf4sc//jEKCgrQv39/AMbWkeEIgXkRVFVFp06dYLFYMHToUHzwwQch22zcuBGtW7dGr169ItTKxqVnz5545ZVXQn499thjAIA//OEPePzxx1nLOnj//fehKAp69erFOtbQqFGjkJeXh++++y64LCcnB8eOHUNaWhrrGIaPP/4YZWVlwafUKrCWNdehQwccOHAgZFlubi4uXLiAjh07AmA9ayMuLg7p6emw2+149dVX0bFjRwwbNgyAwXWsy7wDTdHcuXPFs88+Kz7++GPx5ZdfipdeekkMGzZMTJo0SXg8HiGEEN9++63o1auX+PWvfy127twpnnnmGZGeni7Wr18f4dY3bjt37qw0zxFreXn33HOPWL16tdi2bZvYtm2b+O1vfyvS0tLEH//4x+A2rOPl+f1+cfPNN4uxY8eKTZs2ia1bt4qbbrpJjBgxQpSVlQkhWMfamj17thg1alTIHD0VWMuaefXVV0WPHj3EH/7wB/H555+L999/X0yZMkUMHjxY5OTkBLdjPS9tz5494vnnnxeff/65+PDDD8WiRYtE7969xZdffhmynVF1bHHhaNWqVWLKlCli4MCBYsCAAWLixIni6aefFiUlJSHbbdu2TUyePFn07t1bjBkzRrz22msRanHTUVU4EoK1vJwnn3xSjB07VvTr10/06dNHZGZmijVr1lT6gcQ6Xl5eXp54+OGHxVVXXSX69+8v7r33XpGVlRWyDetYM0VFRaJ3797iqaeeqnYb1vLyNE0T69atE5MnTxYDBgwQw4YNE/fff784fPhwpW1Zz+odPHhQ3H777WLAgAFiwIABYsaMGeLbb7+tclsj6igJIYSRXV5ERERETRnHHBERERHpMBwRERER6TAcEREREekwHBERERHpMBwRERER6TAcEREREekwHBERERHpmCLdACJqnNLS0mq03SuvvIIhQ4bUc2sap+3bt2Pv3r2YN29epJtCRAbiJJBEVKXdu3eHfP3MM8/gq6++wpo1a0KWd+/eHbGxsQ3YssbjiSeewNq1a3HkyJFIN4WIDMSeIyKq0oABA0K+TkpKgizLlZY3Jy6XC1FRUZFuRqNpB1FLxTFHRBQ2r9eLZ555BjfeeCP69OmDoUOH4rHHHkNBQUHIdtdffz3uv/9+fPLJJ7jpppvQr18/jB8/Hp988gkA4O2338b48eMxYMAA3Hbbbdi3b1/I/gsXLsTAgQPx/fffY8aMGRgwYACGDh2KJ554Ai6XK2RbIQTWrl2LKVOmoF+/fhg8eDDmz5+PU6dOhWw3ffp0ZGZm4r///S+mTp2K/v37Y9GiRQCA999/H/fccw+GDx8ebOvSpUvhdDpD2rR27VoAgVuQFb9Onz6N06dPIy0tDW+//XalmqWlpWHFihXBr1esWIG0tDQcOHAA8+fPx+DBg3HDDTfU6lqIyFjsOSKisGiahjlz5uCbb77BzJkzMWjQIJw5cwYrVqzA3r178dZbb8FmswW3P3z4MP76179i9uzZiI2NxcqVKzFv3jzMmjULO3bswMMPPwxJkvD//t//w+zZs/HRRx+F7O/z+TBr1izccccdmDVrFnbt2oVnn30WZ8+exXPPPRfc7ne/+x3eeecdTJ8+HQsWLEBxcTFWrlyJqVOn4t1330VycnJw29zcXDzyyCO499578dBDD0GWA/9ePHHiBEaOHIkZM2YgKioKx44dw/PPP4+9e/filVdeAQDMmTMHTqcTW7ZswZtvvhk8Zps2bXDhwoVa13PevHmYMGECpk6dGgxhtbkWIjIOwxERheWDDz7AZ599hhUrVmDs2LHB5enp6bjtttvw9ttv46c//WlweVFREdavX4+2bdsCANq2bYspU6Zg/fr12Lp1a8htpLlz5+LLL7/E9ddfH1zm8/lw991346677gIAXHvttTCZTFi2bBm++eYbXHXVVdi9ezfWr1+PhQsX4u677w7ue/XVV2PcuHF4+eWX8cgjj4S06emnn0ZGRkbItc2ZMyf4eyEEBg0ahNTUVEybNg2HDx9Geno6UlJSguHEiFuNN910E+bPnx/8urbXQkTG4W01IgrLJ598ArvdjtGjR8Pv9wd/9ezZE61bt8bXX38dsn3Pnj2DwQgAunXrBgAYMmRISDBKTU0FAJw9e7bSOSdNmhTydWZmJgDgq6++CrZJkiRMnjw5pE3JyclIT0+v1Kb4+PhKwQgATp06hV/+8pe49tpr0bNnT/Tu3RvTpk0DABw7dqxmBaolfcAM51qIyDjsOSKisOTn58PhcKBPnz5Vri8sLAz5Oj4+PuRri8VS5XKz2QwA8Hg8IctNJhMSExNDlrVu3RpAoAeook1CCAwbNqzKNnXu3LnK/fXKysrw05/+FFarFb/4xS/QtWtX2Gw2nD9/Hg8++CDcbneVx66rNm3ahHxd22shIuMwHBFRWBITE5GQkIAXXnihyvUxMTGGns/v96OwsDAkIOXm5gIAEhISgm2SJAlr164Nhi+9Hy6TJKnSNjt37sSFCxfw6quv4pprrgkuLykpqXFbrVYrgMCAdb0fBsZLqe21EJFxGI6IKCyjRo3Cpk2boGka+vfv3yDnfO+994JjjgBg48aNABAMMaNGjcLq1auRk5ODCRMmhHWOisD0w/Cxbt26SttWbON2u0MGjycnJ8NqtVaa/+ijjz6qcTuMuBYiCg/DERGFZeLEiXjvvfcwa9YsTJ8+Hf369YPZbMb58+fx1VdfYcyYMcFH0o1gNpvx8ssvw+l0om/fvsGn1UaOHImrr74aAHDVVVfhjjvuwKJFi7B//34MHjwYUVFRyM3NxTfffIMePXqEDBKvysCBAxEfH4/HH38cDz74IEwmE957770qJ3rs0aMHAOD555/HyJEjIcsy0tLSYLFYMHnyZLz11ltISUlBeno69u7dGwxzNWHEtRBReBiOiCgsiqLg2WefxSuvvIJ3330Xq1evhqIoaNeuHQYPHhwMDkYxm8147rnnsHjxYjz77LOw2Wy4/fbb8eijj4Zs98QTT6B///5488038cYbb0DTNLRp0waDBg1Cv379LnuexMRErFq1Cn/5y1/wyCOPICoqCmPGjMGyZctw8803h2ybmZmJb7/9Fq+//jpWrlwJIQQ++ugjdOrUCQsXLgQAvPDCC3A6nRgyZAiee+65kCfwLqeu10JE4eHrQ4io0Vu4cCG2bNmCXbt2RbopRNQC8FF+IiIiIh2GIyIiIiId3lYjIiIi0mHPEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkQ7DEREREZEOwxERERGRDsMRERERkc7/B/PUzXfwXqNXAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sns.set(color_codes=True)\n", + "plt.xlim(30,90)\n", + "plt.ylim(0,1)\n", + "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." + ] + } + ], + "metadata": { + "celltoolbar": "Hide code", + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- 2.18.1