{ "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": 2, "metadata": { "trusted": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.11.15 (main, Mar 11 2026, 17:14:47) [Clang 20.1.8 ]\n", "uname_result(system='Darwin', node='MacBook-Pro-de-Marine.local', release='21.6.0', version='Darwin Kernel Version 21.6.0: Sat Jun 18 17:07:28 PDT 2022; root:xnu-8020.140.41~1/RELEASE_ARM64_T8110', machine='arm64')\n", "IPython 9.10.0\n", "IPython.core.release 9.10.0\n", "IPython.external.pickleshare 0.7.5\n", "PIL 12.3.0\n", "PIL.Image 12.3.0\n", "PIL._deprecate 12.3.0\n", "PIL._version 12.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", "appnope 0.1.4\n", "argparse 1.1\n", "charset_normalizer 3.4.4\n", "charset_normalizer.version 3.4.4\n", "comm 0.2.3\n", "csv 1.0\n", "ctypes 1.1.0\n", "ctypes.macholib 1.0\n", "cycler 0.12.1\n", "dateutil 2.9.0.post0\n", "dateutil._version 2.9.0.post0\n", "debugpy 1.8.16\n", "debugpy.public_api 1.8.16\n", "decimal 1.70\n", "decorator 5.2.1\n", "defusedxml 0.7.1\n", "executing 2.2.1\n", "executing.version 2.2.1\n", "fontTools 4.63.0\n", "http.server 0.6\n", "ipaddress 1.0\n", "ipykernel 7.2.0\n", "ipykernel._version 7.2.0\n", "ipywidgets 8.1.7\n", "ipywidgets._version 8.1.7\n", "jedi 0.19.2\n", "json 2.0.9\n", "jupyter_client 8.8.0\n", "jupyter_client._version 8.8.0\n", "jupyter_core 5.9.1\n", "jupyter_core.version 5.9.1\n", "kiwisolver 1.5.0\n", "kiwisolver._cext 1.5.0\n", "logging 0.5.1.2\n", "matplotlib 3.11.1\n", "matplotlib_inline 0.2.1\n", "numpy 2.4.6\n", "numpy._core 2.4.6\n", "numpy._core._multiarray_umath 3.1\n", "numpy.core 2.4.6\n", "numpy.f2py \n", "numpy.f2py.auxfuncs \n", "numpy.f2py.capi_maps \n", "numpy.f2py.cb_rules \n", "numpy.f2py.cfuncs \n", "numpy.f2py.common_rules \n", "numpy.f2py.crackfortran \n", "numpy.f2py.f2py2e \n", "numpy.f2py.f90mod_rules 1.27 \n", "numpy.f2py.rules \n", "numpy.f2py.use_rules 1.3 \n", "numpy.linalg._umath_linalg 0.1.5\n", "numpy.version 2.4.6\n", "packaging 25.0\n", "pandas 3.0.5\n", "pandas._version_meson 3.0.5\n", "parso 0.8.5\n", "patsy 1.0.2\n", "patsy.version 1.0.2\n", "platform 1.0.8\n", "platformdirs 4.9.4\n", "platformdirs.version 4.9.4\n", "prompt_toolkit 3.0.52\n", "psutil 7.0.0\n", "pure_eval 0.2.3\n", "pure_eval.version 0.2.3\n", "pydevd 3.2.3\n", "pygments 2.19.2\n", "pyparsing 3.3.2\n", "re 2.2.1\n", "scipy 1.17.1\n", "scipy._lib._uarray 0.8.8.dev0+aa94c5a4.scipy\n", "scipy._lib.array_api_compat 1.13.0\n", "scipy._lib.array_api_compat.numpy 2.4.6\n", "scipy._lib.array_api_extra 0.9.1\n", "scipy.interpolate._dfitpack 2.4.2\n", "scipy.linalg._fblas 2.4.2\n", "scipy.linalg._flapack 2.4.2\n", "six 1.17.0\n", "socketserver 0.4\n", "stack_data 0.6.3\n", "stack_data.version 0.6.3\n", "statsmodels 0.14.6\n", "statsmodels.__init__ 0.14.6\n", "statsmodels._version 0.14.6\n", "statsmodels.api 0.14.6\n", "traitlets 5.14.3\n", "traitlets._version 5.14.3\n", "urllib.request 3.11\n", "wcwidth 0.2.14\n", "xmlrpc.client 3.11\n", "zlib 1.0\n", "zmq 27.1.0\n", "zmq.sugar 27.1.0\n", "zmq.sugar.version 27.1.0\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/9m/6k7qyp3d5mz0cz8dg3pz4ffh0000gn/T/ipykernel_1147/3691746940.py:4: DeprecationWarning: numpy.core is deprecated and has been renamed to numpy._core. The numpy._core namespace contains private NumPy internals and its use is discouraged, as NumPy internals can change without warning in any release. In practice, most real-world usage of numpy.core is to access functionality in the public NumPy API. If that is the case, use the public NumPy API. If not, you are using NumPy internals. If you would still like to access an internal attribute, use numpy._core.__version__.\n", " if(hasattr(val, '__version__')):\n", "/var/folders/9m/6k7qyp3d5mz0cz8dg3pz4ffh0000gn/T/ipykernel_1147/3691746940.py:5: DeprecationWarning: numpy.core is deprecated and has been renamed to numpy._core. The numpy._core namespace contains private NumPy internals and its use is discouraged, as NumPy internals can change without warning in any release. In practice, most real-world usage of numpy.core is to access functionality in the public NumPy API. If that is the case, use the public NumPy API. If not, you are using NumPy internals. If you would still like to access an internal attribute, use numpy._core.__version__.\n", " print(val.__name__, val.__version__)\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", "sns = None # seaborn not available in this environment\n", "\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": 3, "metadata": { "trusted": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Current working directory: /Users/marine/Desktop/Module4_Challenger\n", "Files here: ['data_shuttle.csv', 'src_Python3_challenger__1___1_.ipynb']\n" ] } ], "source": [ "import os\n", "print(\"Current working directory:\", os.getcwd())\n", "print(\"Files here:\", os.listdir(\".\")[:50])\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true, "trusted": true }, "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", "\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": { "trusted": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG1CAYAAAAFuNXgAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMRFJREFUeJzt3QmczfX+x/HPmBnLyD5jX7JFIrLLrri69rXcFBIlSpe64fav3ES3dEsqShKuqBsh0bXvEcIVZU8iYSwzZsbMmDn/x+d77znNmYUZs/zO+c7r+Xicxzi/8/ud8ztfvznnPd81wOVyuQQAAMASeZw+AQAAgKxEuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWCXI6ROIjIyU5cuXS8GCBaVz587pOubXX3+VzZs3S/78+aVNmzZSqFChbD9PAADgHwKcmqE4ISFBRowYIYsXL5Z8+fJJaGio7Ny584bHzZs3T4YOHSrNmjWTixcvyi+//CIrVqyQ+vXr58h5AwAA3+ZYs5RmqjvvvFMOHTok3bt3T9cxZ8+eNcFm4sSJsnr1atm1a5epuRkwYEC2ny8AAPAPjoWboKAgGTZsWIaalLSWR0PRkCFDPNueeuop+f7772Xfvn3ZdKYAAMCfON7nJiM0xFSuXFlCQkI822rXru15rE6dOimOiY2NNTe3xMREuXDhgpQoUUICAgJy6MwBAEBmaOWG9tMtW7as5MmTx55wc/nyZSlatKjXtiJFikhgYKB5LDWTJk2S8ePH59AZAgCA7HTy5EkpX768PeGmQIECEhUV5bUtJibGdE7Wx1IzduxYGTVqlOe+hqCKFSvK8ePHrRtlFR8fL+vWrZO2bdtKcHCw06fjdyg/ytBpXIOUoS+I99HvEq210dab9Hx3+1W4qVq1qvzrX/8yTUvuKikNKe7HUqMjsfSWXPHixaVw4cJi2wWpTXba5OZLF6S/oPwoQ6dxDVKGviDeR79L3OeSni4lPj2JX1xcnHz88cdy5MgRc79Tp05m+PeqVas8+8yfP1/CwsKkSZMmDp4pAADwFY7W3CxZssSElQMHDkh4eLgJMqp///5mNFV0dLQMGjRIZs2aJdWqVZNatWqZuXH08aefftp0DH777bfNcb6ULgEAQC4NNzt27DCT8GnPZ72tX7/ebO/Xr58JN3nz5jVz2GiwcdMw07p1a1m7dq1pbtq0aZM0bdrUwXcBAAB8iaPhZsKECdd9XNv83LU5SfXq1cvcAAAA/KrPDQAAQEYRbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYJcjpE/jhhx9k/fr1kj9/frnvvvukdOnSNzxm06ZNcuDAAcmbN680btxY7rjjjhw5VwAA4Pscrbl59913pUGDBrJhwwZZsGCB3HbbbbJx48Y094+Pj5c//OEP0qdPH/nuu+9k9erV0qhRIxkzZkyOnjcAAPBdjtXc/PLLLzJq1CiZNm2aPPLII2bb4MGDze3QoUMSEBCQ4pi1a9fKypUrTW1PzZo1zbZZs2aZYzTgFC1aNMffBwAA8C2O1dwsWbLENCs9+OCDnm2PP/64HDlyRHbv3p3qMfny5TM/CxYs6NkWEhJinidPHroPAQAAB2tutPalUqVKnsCi3LUx+lj9+vVTHNOmTRtTQ9OjRw/p3LmzREVFyYoVK0ztTeHChVN9ndjYWHNzi4iI8DRx6c0m7vdj2/vKKZQfZeg0rkHK0BfE++h3SUbOx7FwExkZmaIZqVChQhIYGGgeS01iYqL5eeHCBTl16pQJN9HR0XLt2rU0X2fSpEkyfvz4FNu1eUtrfWy0atUqp0/Br1F+lKHTuAYpQ1+wyse+S/T73ufDjQYLdy2K25UrVyQhIcGr2Smpf/7zn/LWW2/Jjz/+aGp91NKlS6Vnz55m1FSNGjVSHDN27FjTt8dNX7NChQrSoUOHNGt7/JWmWr0Y27dvL8HBwU6fjt+h/ChDp3ENUoa+IN5Hv0uSZwafDDcaRObPn29qXYKC/nsaR48eNT+rV6+e6jG7du0yj7mDjbrnnntMINLRU6mFG232Str05ab/Yb70n5aVbH5vOYHyowydxjVIGfqCYB/7LsnIuTjWC7dLly6mpkY7Frt9/PHHUq5cOTO8W2lfmcmTJ8v+/fvN/apVq8rx48fl3LlznmO2b99uflarVi3H3wMAAPA9jtXcaFD561//KoMGDZLNmzebfjQ6183nn39u+t2omJgYefbZZyU0NNRM1KdDvmfOnCnNmjWTfv36mT432plYR1y5AxEAAMjdHJ2hWDv6tmvXzsxfExYWJnv37vWMmFI6a/Ho0aOldu3a5r72xdGmqS+++MKMqCpZsqQsXLjQPAcAAIBPLL/QunVrc0uNhhttlkpK++foDMUAAACpYeY7AABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsEOX0CGzdulDVr1kj+/PmlR48eUrNmzRseExERIQsXLpTjx4/LXXfdJd27d5eAgIAcOV8AAODbHK25eemll6Rz585y+fJl+fHHH6Vu3bry5ZdfXveYvXv3mgA0e/ZsyZs3rwk5Dz/8cI6dMwAA8G2O1dwcPXpUJkyYIJ9++qn06tXLbCtevLgMGzZM/vjHP0pgYGCKY65duya9e/eWNm3ayCeffOLZfvjw4Rw9dwAA4Lscq7nRGppbbrlFunXr5tk2cOBAOXXqlOzcuTPVY77++ms5cuSIvPjii17bq1evnu3nCwAA/INjNTcHDx6UihUrSlDQ76dQtWpV8/PQoUPSpEmTFMfs2LFDwsLCpGjRovLGG29IdHS06XOjTVtpiY2NNbek/XVUfHy8udnE/X5se185hfKjDJ3GNUgZ+oJ4H/0uycj5ZDjcbN++XYKDg6V+/fqSGRpMChcu7LVNa3K0OSoqKirVY7RvjtJmqfvuu890Qh4yZIg0bNhQli5dmmqn4kmTJsn48eNTbF+5cqWEhISIjVatWuX0Kfg1yo8ydBrXIGXoC1b52HeJ5oZsCzf79u0zgUI7/w4ePFgefPBB01cmozTIXLp0yWtbZGSkJCQkSKFChdI85ty5czJr1izp1KmT2aavX7t2bdNkpYEnubFjx8qoUaO8am4qVKggHTp0SBGu/J2mWr0Y27dvbwIoKD+uQf/C7zBl6AviffS7xN3yki3h5tFHH5VWrVrJRx99ZGpFnn32WdNv5pFHHjEFkSdP+rrx1KpVS+bMmWOajPLly2e26Ygpdfvtt6d6zB133GF+NmrUyGtbgQIF5NixY6keo8/tfv6k9D/Ml/7TspLN7y0nUH6UodO4BilDXxDsY98lGTmXm+pQfNttt8mrr74qJ0+elM8//9yMYurSpYvceuut8sILL8iJEydu+BwaiOLi4mTevHmebdOnT5dq1aqZfjTq6tWr8swzz3g6GGvNjNa2rF271nPMtm3bJCYmRu68886beSsAAMAymepQrP1jtHlHQ8jp06dN0Jg5c6ZMnDjR1ORMnTo11VoTVb58edMpeMSIEab668KFC/LNN9/IsmXLPH1n9Hl1H2120n412pF4xowZpvZI+9honxsNV0899ZS0bNkyM28FAADk9qHgOpneyJEjpWzZsvL4449L06ZNZf/+/WYot846vGHDBtM35no02Ozatcs0c91///1mlJT+202bm15//XWvZqi+ffuafj9t27aV5s2by+bNm2XKlCk3+zYAAEBur7nRMKGhZvfu3dKuXTt55513pGfPnma2YLe7775bBg0aZCbquxHtX5NWHxut9dFmqeQqVapkOjUDAABkOtzoek4dO3Y0zUGVK1dOcz8NN742Rh4AANgvw+HmoYceStd+pUqVupnzAQAAyPk+N7pY5Z49e7y26bIIupglAACAX4Ub7Uejw711npqkdOkEHcqtSyQAAAD4TbjROWa0w3DSDsRKh2/fe++9ZlkDAAAAvwk3OoJJJ+9Lzc8//2zmvgEAAPCbcKNLLKxfv940QSUmJnq2L1iwQObPn+9Z8wkAAMAvwk2ZMmXk/fffN7MCh4WFSb169czIKF3AUteaqlOnTvacKQAAQHYtvzBgwABp3bq1LFmyRM6cOSOhoaHSuXNnqVGjxs08HQAAgPNrS+kimTpTMQAAgN+HG5fLZRbJ1NmKdWXvpNyLXAIAAPhFuNElFbRTsS6OqX1ugoODvR7XRTQJNwAAwG/CzbJly8yQ7xMnTkiFChWy56wAAAByarTU2bNnzcKZBBsAAGBFuKlbt67s378/e84GAAAgp5ulKlWqJFevXpUhQ4ZI3759pVChQl6Ply1bVipWrJjZ8wIAAMiZcDNv3jz59ttvze3DDz9M8fjo0aNl8uTJN3c2AAAAOR1udG4bHRGVluQLagIAAPh0uNGh38mHfwMAAPj9DMWrV6+WrVu3mrWlunbtKr/++qtcunRJbr/99qw9QwAAgOwcLaW0Wap3794ye/ZsM5mfCgwMlG7duklUVNTNPCUAAIAz4Wbnzp2ydOlSOXjwoDzxxBOe7SVLlpSmTZvKZ599ljVnBgAAkBPhZteuXdKlSxcpVaqUBAQEeD1WtWpVOXz48M2cBwAAgDPhRkdDXb58OdXHdHK/0NDQrDgvAACAnAk3HTp0kK+//tqsCu6uudFVwt9//3354osvTK0OAACA34yWKleunLzzzjvSpk0bKVCggKnJmTlzpkRERMiUKVOkevXq2XOmAAAA2TUUvH///ibcLF68WE6dOiUlSpQwNTY1atS4macDAABwfp6b8uXLy4gRI7LuTAAAAJwIN8eOHZMDBw6k+biOmGIiPwAA4DfhRue4GTNmjNe2+Ph4SUxMNBP5PfvsszJp0qSsPEcAAIDsGy319NNPy9WrV71uMTExZqSU9rkZN25cRp8SAADA2eUXktMRU927d5fWrVvL559/nhVPCQAA4Fy4cStevLicOHEiK58SAAAge/vc6Mrf58+f99qWkJAg+/btkw8++EDee++9jD4lAACAc+Hmww8/NJ2GkwsODjarhffq1Surzg0AACD7w83QoUOld+/e3k8SFCSlS5c2PwEAAJyU4TRSuHBhcwMAAMgVk/glxYR+AADA58PNsmXLTJ+buLi4/z5BUJBcu3bN0+9Gh4W7jRw5Ul555ZWsPF8AAICsHQo+ZMgQM1nfiy++KGfOnDGzE4eHh8ubb74pFSpUkJ9//lmuXLlibgQbAADg8+Fm9erVUrlyZXnppZekVKlSnvltdObie++9VxYtWpQd5wkAAJA94UYn6UurQ7FuZxI/AADgV+GmYcOGsnDhQtP3JqlNmzaZOXD0cQAAAL/pUNy0aVN57rnnpEePHlKyZEkpW7asnDt3Tk6ePCkjRoyQbt26Zc+ZAgAApMNNzbqnnYkffvhh+fe//y2nTp0yE/i1a9dObr/99pt5OgAAgCxz01MKa6diXW4BAADAinCjo6a2bt0q9erVk65du8qvv/5qFtWk9gYAAPhVh2KlNTa6vtTs2bNl48aNZltgYKDpbxMVFZXV5wgAAJB94Wbnzp2ydOlSOXjwoDzxxBOe7dq5WDsbf/bZZxl9SgAAAOfCza5du6RLly5mAr+AgIAUa0kdPnw4684OAAAgu8ONrh11+fLlVB/bv3+/hIaGZvQpAQAAnAs3HTp0kK+//lq2bdvmqblxuVzy/vvvyxdffGFqdQAAAPxmtFS5cuXknXfekTZt2kiBAgVMTc7MmTMlIiJCpkyZItWrV8+eMwUAAMiuoeD9+/c34Wbx4sVmEr8SJUqYGhtdLRwAAMCvwo2OhgoPD5dhw4aZ5RYAAAD8us+NTtS3d+/e7DkbAACAnA43uobUunXr0hwxBQAA4FfNUrrMgnYi1v41nTp1krCwMK/HW7VqJX/84x+z8hwBAACyL9ycOXPGTOCntxMnTphbUlWqVMnoUwIAAOR8uDl58qTExcVJnz59zA0AAMCv+9x8+umnMm3aNM/9OXPmyJtvvpld5wUAAJBz89yos2fPmiYqAAAAX3LT4SarLFq0SNasWSP58+c3zV26snh66JIP48aNMx2ctQapWLFi2X6uAADAwqHgWempp56Sxx57TMqUKWPut2zZUubNm5euY//xj3/I3LlzZfbs2RIVFZXNZwoAAKysuVmwYIFZMFPpsgvx8fGe+279+vWT4cOH3/C5Dhw4YNaoWrFihfzhD38w20JCQuTPf/6z9O3bV4KDg9M8dufOnWYdq9dff13+9Kc/ZeQtAAAAy6U73Nx5553SsWNHz/2aNWumul/JkiXT9XzLly83TUnt27f3CkYTJkyQ7du3S4sWLVI9LjIy0uynq5BfLwABAIDcKd3hpkOHDuaWVY4cOSIVKlSQPHl+bxmrXLmy+Xn06NE0w402Y2nIuu+++2T16tU3fJ3Y2Fhzc9PVy5XWOunNJu73Y9v7yimUH2XoNK5BytAXxPvod0lGzsexDsUxMTFSsGBBr20FChSQwMBA81hqPvroI7OulTZLpdekSZNk/PjxKbavXLnSNIPZaNWqVU6fgl+j/ChDp3ENUoa+YJWPfZdER0f7frgpXLiwWYQzKV2vKiEhQYoUKZLqMa+88orpfKwrkisdKaVGjRol3bt3T7X/zdixY83jSWtutMZIa6H0HGyiqVYvRm3qo8mO8uMa9D/8DlOGviDeR79L3C0vPh1u6tSpIzNnzjRJzF2D8v3335uftWvXTvWYv//973LlyhWvTslaA9OsWTOpVq1aqsfky5fP3JLT/zBf+k/LSja/t5xA+VGGTuMapAx9QbCPfZdk5FwcGwrerVs3CQgIkBkzZni2vf322ybYaPBRGnwGDhwomzdvNvd79+5t7rtv7j5AOj9O48aNHXonAADAlzhWc6MLb+qIJ+0gvGzZMrl48aJZv0qHhrvpWlY6j02bNm3S7GAMAADgMzMU9+/fX9q1aydbtmwxTUdt27aVQoUKeR7X5qpZs2alGWzuuOMO83jx4sVz8KwBAIAvc3z5hbJly6a5ynjevHlN81NatHPx9R4HAAC5j6PLLwAAAGQ1wg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBustCxc1dk3cGzcvx8VFY+LQAgnX763+fvifBoyiwXC3L6BGxwKTpOnpq/RzYePufZ1qp6mEztd5cUCQl29NwAIDd9Dm8/dlZeayzSaeomaVKlJJ/DuRQ1N1lAf6G2HDnvtU3vPzl/d1Y8PQCAz2FkAOEmC5qitMYmweXy2q73dTtNVACQvfgcRnKEm0w6ceH67bo/hdP/BgCyE5/DSI5wk0mViodc9/FbSxTM7EsAAPgcRgYQbjKpStgtpvNwYECA13a9r9srhxJuACA78TmM5Ag3WUBHRTWvFuq1Te/rdgBA9uNzGEkxFDwL6HDvOYMbm87D2sdGm6KosQGAnP8cPnLmsuzfvl6+erKlVCtdhP+CXIpwk4U00BBqAMA5lUqEyP7//UTuRbMUAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFWCnHzxhIQEmTlzpqxZs0by588vffv2lU6dOl33mH379sknn3wix44dkwoVKsgjjzwitWrVyrFzBgAAvs3RmptBgwbJyy+/LG3btpWaNWtKr169ZNq0aWnuP2vWLHnooYekSJEi0rNnT7l69arUrVtXvv766xw9bwAA4Lscq7nZvXu3zJ07VzZu3CgtW7b0bB83bpwJPVqTk5zW6gwcOFACAgLM/fvvv1/Onj0rEydOlI4dO+bo+QMAAN/kWM2N1raEhYVJixYtPNt69+4tly5dkm3btqV6TMmSJT3Bxq106dISERGR7ecLAAD8g2M1N8ePH5dy5cp5hRXtQ+N+rE2bNjd8jt9++830vxk6dGia+8TGxpqbmzsIxcfHm5tN3O/HtveVUyg/ytBpXIOUoS+I99Hvkoycj2PhJi4uTkJCQry25cuXTwIDA81jNxIdHW363ZQvX16ef/75NPebNGmSjB8/PsX2lStXpnh9W6xatcrpU/BrlB9l6DSuQcrQF6zyse8S/d73+XBTtGhRuXDhgtc2bZLSEVT62PXExMRI165dzf7r1q27bkgZO3asjBo1yqvmRmuIOnToIIULFxabaKrVi7F9+/YSHBzs9On4HcqPMnQa1yBl6AviffS7JCNdUBwLNzrKSUdGRUZGSqFChcy2PXv2eB5Li46Q6tatm5w+fdoEG+2Hcz1aG6S35PQ/zJf+07KSze8tJ1B+lKHTuAYpQ18Q7GPfJRk5F8c6FGtA0RFRb731lrmfmJgokydPlkaNGplh4SoqKko6d+5s5sFR2ndGjzt16pQJNqVKlXLq9AEAgI9yrOamePHiMm/ePOnfv78sWrRILl++bLYvX77cq2rsq6++MqOo1GuvvWb6yjRp0kQGDx7s2e+WW26RBQsWOPAuAACAr3F0hmKtlTl58qTs3LnTNB1prU3SaicNLV9++aXUq1fP3O/Tp4/cddddKZ7Hl6rNAABALg43Svvb6AzFqQkKCjIByE2bq9xNVgAAAKlh4UwAAGAVwg2yzLFzV2TdwbNy/HyUI8eDMvR3Ww6fMz+/OXre6VMB/JrjzVLwf5ei4+Sp+Xtk4/8+mFWr6mEytd9dUiQkONuPB2Xo706ER0n3d7dIdGycvNZYZMjcXRKSL68sHd5CKpSwc7JRIDtRc4NM02Cy5Yj3X5p6/8n5u3PkeFCG/k6DzcVo76nl9X7Xdzc7dk6APyPcIFO0KUlrXBJcLq/tel+336iJKbPHgzL0dxsOnk0RbNx0+6YkNZoA0odwg0w5ceH6a338FB6VrceDMvR3e365dN3Hv/v5Yo6dC2ALwg0ypVLx6/cHuLVEwWw9HpShv6tX/vpr6dWvWCzHzgWwBeEGmVIl7BbT+TcwIMBru97X7ZVDC2br8aAM/V3rGiWlWBod53V7y+phOX5OgL8j3CDTdFRT82qhXtv0vm7PieNBGfo7HRWVPODofd0OIOMYCo5M0+HacwY3Np1/tY+MNiVlpMYls8eDMvR3Otx79wsdZOOPv8rFg9/KjIcaSKuaZZw+LcBvEW6QZTSQZCaUZPZ4UIb+rlnVUFl+8L8/Adw8mqUAAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVR8PN1atXZeLEiXLPPfdIp06dZM6cOdlyDAAAyD2CnHzxBx54QH744Qd55ZVX5OLFizJ8+HA5ffq0jBkzJkuPAeAfjp27IicuRMutJQpK5dCCGT7+029/lm+Oh0vzqqHSp2GFHH/9zB6/5fA58/Obo+elVc0y4gSnyyCzfjofZX6eCI+WaqWLZPj4DQfPyp5fLkn9isWkZfUwyWn+Xv6+cg6OhZtvvvlGlixZIjt27JCGDRuabVFRUfL888/Lk08+KQULFsySYwD4vkvRcfLU/D2y8X9f7qpV9TCZ2u8uKRISfMPj9/1ySXq8t1WuJbrM/cW7T8vYRftk6fDmUqtckWx//cwefyI8Srq/u0WiY+PktcYiQ+bukpB8eWXp8BZSoUSI5ASnyyCz3K+//dhZU4adpm6SJlVKZvj/4GJ0vGdbsZDgHPs/8Pfy95VzcLxZas2aNVK6dGlPSFHdunUzYWXbtm1ZdgwA36cfiFuOnPfapvefnL87XccnDTZuer/ru1ty5PUze3zyL1Wl97u+u1lyitNlkFn+/n/g7+XvK+fgeM3NiRMnpGzZsl7bypUr53ksq46JjY01N7fLly+bnxcuXJD4eO8L2d/p+4mOjpbw8HAJDs7ZlGwDys+ZMvw5PEq2/vCTBKTygbT1hyuy+3AZqVg87b+cl+45JRJ3Jc0Ps9lr90nnumWz7fUze/y2o+clMuKSOTYo0SXR0YkSFJ9HEhIDJDJCZMXOQ9K4cgnJTk6XQWZ5vX4e7zLM6P9Bcjnxf+Br5R/vwO9xekRGRpqfLpf3HzKpcjlkwIABrmbNmnltS0xMdOXJk8c1ffr0LDvmxRdf1FLgRhlwDXANcA1wDXANiP+XwcmTJ2+YMRyruSlevLipPUnq0qVLkpiYKCVKlMiyY8aOHSujRo3y3Nd99Tl0/4AAzZj2iIiIkAoVKsjJkyelcOHCTp+O36H8KEOncQ1Shr4gwke/S7TGRmtvkrfgpMaxcFO/fn2ZOnWqGfFUrFgxs2379u3m51133ZVlx+TLl8/ckipatKjYTC9GX7og/Q3lRxk6jWuQMvQFhX3wu6RIkSK+3aFYOwJryNAh3SouLk5effVVad26tVStWtVsu3LlijRt2lS++uqrdB8DAAByN8fCTaFChWThwoXyySefmOovHQWlnX2TTsp37do1UzNz7ty5dB8DAAByN0cn8WvVqpX8/PPPZlI+bTq67bbbvB7XMKNz2yStlbnRMbmZlseLL76YohkOlB/XoH/gd5gy9AX5LPguCdBexU6fBAAAQFZh4UwAAGAVwg0AALAK4QYAAFjF0Q7FyLi33npLFixY4LVNR47961//ytA+ud3OnTtl2rRpcvz4cTNH0rhx47wmgrx69ar84x//MOuZ5c+fX+6//355+OGHHT1nX6JLmmj5rVy50kzPPmDAAOnZs6fn8ffeey/FKMaSJUvK0qVLJbf797//bTprpkbLzD1IQqe6ePvtt83+QUFBpnwfffRR6yYfvRkPPPCA/PTTTym2N2/eXN544w3z748//limT5+eYpDKqlWrcuw8fdmpU6dkypQpsnfvXsmTJ4+ZR27kyJHm9zTpiGWdW2758uUSGBhopmN57LHHzP6+jnDjZ/QXWj/oJk+e7NmmX74Z3Sc3++yzz+Shhx6SZ555xvz88ccfZejQoWaagaQfnjoiT+dU0kkjhw8fLqdPn5YxY8ZIbqfBr127dmbB2ueff95MqDljxgy55ZZbpEOHDmYfHdGos4Hrl7ObP4+8yEoNGjQwf4AkNWHCBNm1a5dUqVLFs23QoEGyZcsWee2110xZ60zruoae7pvb6czzMTExnvs6XUjXrl3N763bL7/8YuZK+/DDDz3b9HMRYmb5vfvuu02QHj16tAkxEydONJ+B+/bt86wnpUFG/4B5/fXXTdjWa/Do0aNe3y0+KyPrQcF5I0eOdHXq1CnT++RWly5dchUuXNj10ksveW2PiYnx/Hvr1q1m/ZIdO3Z4tr355puuggULuq5cueLK7V5++WVXsWLFXOfOnUuzDJ977jnXPffc48DZ+R8tNy3PMWPGeLb95z//MdfgunXrPNtmzJjhyps3r+vChQsOnanvmjx5sitfvnyu8+fPe12nTZo0cfS8fNWGDRvM9XXkyBHPNv2802179+419w8dOmTur1ixwrPP3LlzXUFBQa4zZ864fJ3v1y0hhT179pi/nLt3724Stf4lfTP75EbaLKLrpgwbNsxre9KaLW2K0gkiGzZs6Nmm1bH61/O2bdskt5s7d6706dNHQkNDvbYnrx3cv3+/3HPPPabsdCbxpH9p43eLFi0ya+QNHjzY6xrUae91Xi83LUf963njxo0UXzIfffSRabZLvsag1jK0b99eunTpIi+//LKpyYFIjRo1TE3r1q1bPcWxefNm8ztduXJlzzWov9P33nuv1zWotTzr16/3+WKkjs7P6MXWv39/adu2rfz222+m2USbWfQidVclpmef3OrAgQNSsWJF0xSlVa7R0dGmz402Ubm/rLXqP/nCbOXKlfM8lpvpB9uhQ4dkxIgRpklKm01KlSol/fr1Mx98SZugHnzwQfPBeP78eVPlPX/+fPn2229pnkpm5syZ5ne1WrVqnm16nWnATtq3ISwszPz+5vZrMDmd6FV/r9955x2v7VpWf/rTn6Rjx44mPP7973+Xf/7zn/Ldd99JwYIFJTcrVaqUCSi9e/c2v5sJCQnmWtu0aZPpl6T0OtP+N0mb8vQxDUV+cQ06XXWEjLl69arX/Z9++slUVc+aNStD++RWTz/9tGleqlevnmvhwoWur776ytWyZUtX1apVXZGRkWafAQMGuJo1a+Z1XGJioitPnjyu6dOnu3IzbZbTj42iRYu6xo4d61qzZo3rjTfecOXPn981derUNK/BU6dOuQoUKOCaNm2aA2ftu44dO+YKCAhwzZ8/32v78OHDXXXq1Emxf6FChUwTDH43ePBgV7Vq1czvaFLJr0FtRi1SpIjr9ddfz/XFd+XKFVejRo1M07F+Bi5ZssR19913u9q1a+eKj4835TN69GhX9erVU5RVWFiYa8KECT5fhtTc+JnknTIrVaok1atXNz3eM7JPblW8eHHTvKSdDLVjp2rcuLH5C0VHBPTt29fsc+HCBa/j9C8/7SCbvNo7twkJCTE1g7qgrf7Fp7T58+TJk6bzsNbopHYNak1YrVq1uAZTaU7R661Hjx5e21O7BrVJSjuC5vZrMCltZvr0009NLWLyUWTJr0Gtma1Xrx7XoIipwdKOw9oRW2tiVLNmzUyNjjaTpvU5qAsa6DZ/uAbpc+Pn9Av37Nmz161mTc8+uYW7H03SZif9JdYPQh0VpXRIpLbVu+8rXcBVaRNWbqZfIFo+yZvtypQp41VeyemHojaRcg3+TpsCdLiyTjGQ/ItYy1hH5+nNTZv0VG6/BpPSYKP9CQcOHJiu/c+cOcM1KGICivbpcgcbpYElb968nkCj12B4eLiZLsNNR/TpdesX16DTVUfIGB0B4B6Vcu3aNdM0EBgY6NqzZ0+G9smt4uLiXFWqVHH99a9/9Wz74IMPXMHBwa4DBw6Y+xEREa7Q0FBTLatiY2NdrVu3Nje4XLNnzzZV09rcqS5evOiqW7eu64EHHvAUzyuvvOKKjo42/05ISDCj07RZb/v27RTh/2hzgH4E79+/P0WZ6O9vuXLlXI899pi5r00FHTt2dDVo0IDyS0Kbj3v27Jlqmbz66quepmZtstLmKC3vtWvX5voy3LRpkymLpM2h7733nvkddY+W0s+9W2+91TVw4EDP73G3bt1Mc2nyJkBfRLOUn9FaGJ2QT6sP9S9hTd5ajVi3bt0M7ZNbaSfDL774wnSk01E/+hez/qUya9Ysuf322z2d5nS+B50zQ/8y1GYsbdpbsmSJ06fvE7SmQTtw1q5d24ys0HmVdPI0nezLTTsnaplpc5/WGmpTlnZq1yZA/N6RWMtNm+uS0/LSa1BHpWlndp00UTsYL168mOL7H52HSjsTr1ixItUy0TLU5nh384rWOs6bN8903s7tWrRoYSY7HDJkiLzwwgumNkbLSCfmvPPOO80+Wouj12CvXr1Mzazuo2Wp16A/TCTJquB+KD4+Xg4fPixFihQxzQOpXWjp2Sc302YSHTGlX8I6cVpqo8h0ZJB+gGoAcs8ai99pM5QGm/Lly5uRPGmNrNKwqPtwDXrTKn4Nf/qHSFr0C0WvQR2xosN3KcPfaWg+duyYCcxpzZir5aefg9pXTK9Bf5hZNyfFxsaaZictl1tvvdUEmuT0j2W9BnWfmjVr+s01SLgBAABWIcYCAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKzCDMUAUrVnzx4z0eH1dO/e3cwEa6Pdu3ebhRlbtmzp9KkAyCDCDYBU7d+/X7788kvPfV22QmcodS9Toe677z5rw83s2bPlyJEjhBvADzFDMYB0CQ0NlWeeeUbGjBmTIgTpFPe6BpKuFqxLBbjt3LnTTPFer149UxOkSzZoTYguC6K1Ilu2bJHAwEC5++67zRT5qR2nNSiXL1+WVq1amaUckkvv62/dutWsIN2lSxc5ePCgeV6lqyPXqVPHaxkEXTtL36uuIv2Xv/zFbNPX16nqdTkO9+ryat++fXLu3Dlp167ddV/TvSzKtm3b5NKlSyYkVqtWjasPyAbU3AC4KdHR0dKvXz/zZd6gQQNTy6Ff/Frbo+v4qA8//FA2bdpk9tVaH12L6vz582bRvpdeesls0+M04OzYscMs8uo+ToNPTEyMWfNGQ4YGiFWrVnkW9kvv62vA0LChz6NhQoOGhiH3IpQaNDZv3izPPfec/N///Z/ZpuHn6NGjJoC599O1nd58800T8pKGG11cVY93h5u0XlNDULdu3Uyw0yCliz5qs94HH3zgN+v1AP6CcAPgpowdO9asmK6LF2qo0MVI+/fvL08//bR8/vnnnv00JOgikXfccYdZTFNrLIYNG2ZqTnRBUg0BVatWlTlz5sgTTzzhOe7777+X5cuXm6YvXbxPV2nXxzVIZOT1tWZHg0TSFck7d+5sbklfSwOLrhav59ejRw/ZsGGDCUwLFizIcNkkf0193/qcev5aI6TCw8Olbt26ZqVqPW8AWYdwAyDDNGxon5SBAweamhINFnorU6aMzJo1y2vf5s2bm2BjPnCCgkyIiIiI8Ky0rsFEm5N0BfGkateubYKN0hWJn332WRMWTp06ZV4nva/fpEkTr2DjprUy3333nfz2229m9egSJUqYWqCkfYpuVvLX3LhxowlhpUuXNsHLfb5aq7Nu3TrCDZDFCDcAMkybcrQfjDa1aJNRUu3btzfhRwOJKlasmNfjGmZS26Z9U5LSJp2kKleubH6eOHFCChQokO7X18CTnAaihx9+2LyG3vT1tZnr7NmzWXI1JH9NbY7TYLds2TKv7Rp2tLkLQNYi3ADIsIIFC5p+MgMGDDAhITto5+PU7mufl4y8fmr9WUaOHCnjxo0ztUFu2k9Ha1OuRwOTBqekkoey1F5TOy1rZ+L333/f9LkBkL2YxA9AhmlNR5s2bUxn2ORf9tpslBW0g3HS51q0aJGpEdEanMy+vjZFJa0x0RFMyY/Tzs3Jg4uOyNJ+OG4ahrRvzo20bt3aDJnXcJOUNofpuQDIWtTcALgpU6dOlbZt25ovbu3sqzUT2n+kZMmSMmPGjEyXqoYLbWIaPny4CR6TJ082YSY4ODjTr6+jlkaPHi2nT582NUJTpkzxjNRy075Buv3dd981/XF0KPiDDz5oXvfPf/6z6UekI6m0maxWrVrXfb2wsDBznHak1r5FzZo1M+9p4cKF8re//c2cD4CsQ7gBkC49e/b0+hLXjrc6Kujjjz82tSwaAIYOHSqdOnXy7NOoUSOJjIxM0dk2eW1LixYtUsxho8FlxIgRpn+M9q9ZsmSJp4NxZl5fffTRRzJ9+nRTY6P9f/Q1NGjo0HQ3DRwabHR0lnaA1poeDTw6tH3+/Pnyn//8Rx577DHTLJZ0Jue0XvPRRx81j+noK+1gXKVKFTOMPOlrAsgaTOIHwOc8/vjjZj6cpEO6ASC96HMDAACsQrMUAJ+TVtMOAKQHzVIAAMAqNEsBAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAALHJ/wPPVMov7Mf/dAAAAABJRU5ErkJggg==", "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": 6, "metadata": { "trusted": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/miniconda3/envs/python/lib/python3.11/site-packages/statsmodels/genmod/families/links.py:13: FutureWarning: The logit link alias is deprecated. Use Logit instead. The logit link alias will be removed after the 0.15.0 release.\n", " warnings.warn(\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Wed, 29 Jul 2026 Deviance: 3.0144
Time: 20:33:22 Pearson chi2: 5.00
No. Iterations: 6 Pseudo R-squ. (CS): 0.04355
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
" ], "text/latex": [ "\\begin{center}\n", "\\begin{tabular}{lclc}\n", "\\toprule\n", "\\textbf{Dep. Variable:} & Frequency & \\textbf{ No. Observations: } & 23 \\\\\n", "\\textbf{Model:} & GLM & \\textbf{ Df Residuals: } & 21 \\\\\n", "\\textbf{Model Family:} & Binomial & \\textbf{ Df Model: } & 1 \\\\\n", "\\textbf{Link Function:} & logit & \\textbf{ Scale: } & 1.0000 \\\\\n", "\\textbf{Method:} & IRLS & \\textbf{ Log-Likelihood: } & -3.9210 \\\\\n", "\\textbf{Date:} & Wed, 29 Jul 2026 & \\textbf{ Deviance: } & 3.0144 \\\\\n", "\\textbf{Time:} & 20:33:22 & \\textbf{ Pearson chi2: } & 5.00 \\\\\n", "\\textbf{No. Iterations:} & 6 & \\textbf{ Pseudo R-squ. (CS):} & 0.04355 \\\\\n", "\\textbf{Covariance Type:} & nonrobust & \\textbf{ } & \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "\\begin{tabular}{lcccccc}\n", " & \\textbf{coef} & \\textbf{std err} & \\textbf{z} & \\textbf{P$> |$z$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n", "\\midrule\n", "\\textbf{Intercept} & 5.0850 & 7.477 & 0.680 & 0.496 & -9.570 & 19.740 \\\\\n", "\\textbf{Temperature} & -0.1156 & 0.115 & -1.004 & 0.316 & -0.341 & 0.110 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "%\\caption{Generalized Linear Model Regression Results}\n", "\\end{center}" ], "text/plain": [ "\n", "\"\"\"\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: Frequency No. Observations: 23\n", "Model: GLM Df Residuals: 21\n", "Model Family: Binomial Df Model: 1\n", "Link Function: logit Scale: 1.0000\n", "Method: IRLS Log-Likelihood: -3.9210\n", "Date: Wed, 29 Jul 2026 Deviance: 3.0144\n", "Time: 20:33:22 Pearson chi2: 5.00\n", "No. Iterations: 6 Pseudo R-squ. (CS): 0.04355\n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import statsmodels.api as sm\n", "\n", "data[\"Success\"]=data.Count-data.Malfunction\n", "data[\"Intercept\"]=1\n", "\n", "logmodel = sm.GLM(\n", " data['Frequency'],\n", " data[['Intercept','Temperature']],\n", " family=sm.families.Binomial(link=sm.families.links.logit())\n", ").fit()\n", "\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": 7, "metadata": { "trusted": true }, "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", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Wed, 29 Jul 2026 Deviance: 18.086
Time: 20:33:27 Pearson chi2: 30.0
No. Iterations: 6 Pseudo R-squ. (CS): 0.2344
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" ], "text/latex": [ "\\begin{center}\n", "\\begin{tabular}{lclc}\n", "\\toprule\n", "\\textbf{Dep. Variable:} & Frequency & \\textbf{ No. Observations: } & 23 \\\\\n", "\\textbf{Model:} & GLM & \\textbf{ Df Residuals: } & 21 \\\\\n", "\\textbf{Model Family:} & Binomial & \\textbf{ Df Model: } & 1 \\\\\n", "\\textbf{Link Function:} & Logit & \\textbf{ Scale: } & 1.0000 \\\\\n", "\\textbf{Method:} & IRLS & \\textbf{ Log-Likelihood: } & -23.526 \\\\\n", "\\textbf{Date:} & Wed, 29 Jul 2026 & \\textbf{ Deviance: } & 18.086 \\\\\n", "\\textbf{Time:} & 20:33:27 & \\textbf{ Pearson chi2: } & 30.0 \\\\\n", "\\textbf{No. Iterations:} & 6 & \\textbf{ Pseudo R-squ. (CS):} & 0.2344 \\\\\n", "\\textbf{Covariance Type:} & nonrobust & \\textbf{ } & \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "\\begin{tabular}{lcccccc}\n", " & \\textbf{coef} & \\textbf{std err} & \\textbf{z} & \\textbf{P$> |$z$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n", "\\midrule\n", "\\textbf{Intercept} & 5.0850 & 3.052 & 1.666 & 0.096 & -0.898 & 11.068 \\\\\n", "\\textbf{Temperature} & -0.1156 & 0.047 & -2.458 & 0.014 & -0.208 & -0.023 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "%\\caption{Generalized Linear Model Regression Results}\n", "\\end{center}" ], "text/plain": [ "\n", "\"\"\"\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: Frequency No. Observations: 23\n", "Model: GLM Df Residuals: 21\n", "Model Family: Binomial Df Model: 1\n", "Link Function: Logit Scale: 1.0000\n", "Method: IRLS Log-Likelihood: -23.526\n", "Date: Wed, 29 Jul 2026 Deviance: 18.086\n", "Time: 20:33:27 Pearson chi2: 30.0\n", "No. Iterations: 6 Pseudo R-squ. (CS): 0.2344\n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logmodel = sm.GLM(\n", " data['Frequency'],\n", " data[['Intercept','Temperature']],\n", " family=sm.families.Binomial(link=sm.families.links.Logit()),\n", " var_weights=data['Count']\n", ").fit()\n", "\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": 8, "metadata": { "trusted": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG1CAYAAAAr/fRyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMq9JREFUeJzt3Qd4VFX6x/E3JJAQei/SQu8qAitdhQACkWKhiOyuoiiyLqKroiKCBRuKivK3YAPRBSlKkaYgIE3RqIiiIARBykpoEkIC3P/znvXOzqRPMsNJMt/P88wT5s5tOdyZ+eXcU8Icx3EEAADAkiK2DgwAAKAIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsisjNRt988418/fXX0rlzZ6lbt2626+u4ahs2bJCEhARp0KCBtG7dOjeHBQAAoR5GNm7cKHfddZccP35ctm7dKjNmzMg2jCQnJ0tcXJxZv23btrJu3Trp1q2bzJo1S8LDw/N6/gAAIJRu05w+fVqefPJJ+e6773K8zTPPPGOCSHx8vHz44Ycm0CxatEjeeuut3JwvAAAI5TDSpUsX6dixo18H0BqQgQMHSpUqVcxzvU3Tq1cvsxwAACBXbUZyKjU1VbZv3y533nmnz/LmzZvLyy+/nGUNjD5c586dk8TERKlQoYKEhYXxvwYAQAGgbUZPnDgh1atXlyJFitgJI3/88YcJEmXLlvVZXr58eTl27Fim202aNEkmTJgQzFMDAADnya+//io1atSwE0aKFy9ufp48edJnuaYk97WMjB07VsaMGeN5rsGlVq1asmvXLilVqlTg0lrSaflszRrp0rmzRBQNalEUeGdSz1BWlBfXlmW8DymrYF5XPbpeJsWKFQvovvX7PiYmJtvv7qB+A0dFRZmqmd27d/ss1+f16tXLdLvIyEjzSEtrVEqXLh2w8yuTmirlSkVLjWqVpWjRogHbb2Gkt9woK8qLa4v3YUHBZ5b/ZVWxYsWAfxe6+8uuiUXABz3TXjPejVP79Okjc+fOlTNnznhqSRYuXGiWAwAA+FUzsn//flm2bJnn+dq1a03IaNiwobRv394sW7BggUyZMkWGDBlino8bN86ML9K7d2/p2bOnzJkzx1TXjB49mtIHAAD+1YwcOXJEVq9ebR5//etfTY8X/bf2mHFddNFFniCitMGKO1rrjz/+KAMGDJAvvvgiXaNWAAAQmvyqGWnatGm2g5X169fPPLzpGCMPPPBA7s4QAFBgnT171rRJON/0mBEREWYUcD0HBKestE1IIEZTpwsJACDgtMfigQMH5OjRo9aOX7VqVdOllPGpgltWeqdDt89LORNGAAAB5waRypUrS3R09HkPBDrGlY51VbJkySwH24Lkuqw0xCQlJcmhQ4fM82rVquW6OAkjAICA0qp+N4joyNm2vmBTUlLMEBOEkeCVlTtmmAYS/f/O7S0b4iIAIKDcNiJaI4LCL/rP/+e8tA0ijAAAgoK2GqEhLAC34AgjAADAKtqMAADwZ1uXn3/+OcN2EbVr16aMgogwAgDAnwN7NmnSxAzWWaJECZ/BPN9//33KKIi4TQMAgJdp06aZEcPdhwYRHRBM/609T06dOiU7d+70mZFeG2/qJLDa1TUziYmJnm6w2ttIZ6J3HTx4UH777bd04ch7nZwc6/fff5c9e/aYf+t5pt1n2pqghIQE87t5z7LrPaq6S39X9/cPBsIIAADZ2Lp1q6k10dHEdTyNXr16malO1OTJk81I4506dTI/+/fvL4cPH/Zsq3O43Xjjjabrq9ay6Kz1d9xxh8TFxXnWGT9+vFnm7e233/ZZJyfH0rnh+vbtK1dffbXExMSYkdMbN24sO3bsSLcfnaW3TZs2Ur16dXN8DRoaknSbjRs3+qz/1FNPmWMFq5s0YQQAEHRmgKyUM+f1cSrlrDmuv/bu3etTM6K1Bd4z0+/bt8/UHnTs2NFMkfL888/L+vXrzQimOtib1lyMGjXKp6Zl8eLFJtBoTYWGjNmzZ/t9Xm/l4FjuOcbGxprXtcblggsukPvvv9/z+muvvSYPPfSQ/Pvf/zY1NbpOuXLlzHxzNWvWlB49esgbb7zh83+n56yBKlhoMwIACLpTqWel6UP/m/X9fNn6cKyU9HMgrscee8zUMLieffZZU6uhJk6c6NOe5LnnnjOTw+rcLtr4Vb+4Bw4cKLfccoupadCahFdffVVuvfVWU0OhNMRcd9118tVXX/l1Xs/l4Fiqfv365ngqMjJSrr32WlMT4tLfbcSIEdK9e3fP/DKjR4/2DGA2fPhw+fvf/27W0zFEVq5caQLYsGHDJFgIIwAAeNGajD59+viUyZdffml+1q1b12f5tm3bzG2SBQsW+CzX3jfa5kNHoNX2Jc2aNfN5vXnz5n6HkW05OJbSmhBvOsy7d+2O1uqMHTs20+PorSENJnPnzpUbbrhB3nzzTendu7e5LRQshBEAQNAVLxou2yb2OG8lrTUFJ46fMMcNpLTDnWutgt4CGTlyZKbb6DDr3o1EVdrGp2EZDBymbU38PVZOFCtWLMuGtnqcv/71r+ZWjQaT+fPnm1s6wUSbEQBA0OmXbXSxiPP6KF4sPOijwHbo0EE++OCDdMu9h0a/+OKL5bPPPvN5ffXq1T7PK1WqJPv370/X9sPfY+VE+/btZcmSJenCm3f7Gr1Vs2bNGnn00UelTJkypsFuMBFGAADIpUmTJsmmTZvM7QwNHGvXrjXLBgwY4FlHe+DMmDFDnn76adm8ebPcc8898vnnn/vsp1evXrJhwwZ56aWXZMuWLfL444+nCx45OVZOaMBYtmyZafiqjWG1ce1VV11lugK7GjRoYHrsaFsTbSui7VSCiTACAIC2W4iIkEaNGpk2FmlpGwp9Le1tmtatW5v2JHrr484775QHH3zQ3AJ55513POtcccUVZqySRYsWyT//+U9z+0Ubwnq79NJL5d133zUBRLv4ahsPbbDq3UYlJ8fSGpZatWr57FtrNjRceB9Lg4+Oe6K3fF5//XUZN25cuokN3QarwexF46LNCAAAIlK2bFnTlTcj2gA1s9d0/JHp06dnWYZae+FdgzF16tR06wwePNg8vN12221+HUvDTlra7iPteCU63smsWbM8t2iOHz+ebju9TaS3htxeQMFEGAEAAB46Fsp3331nGrBqbc35QBgBAOA800HG0nYTzi+efPJJ0yblvvvuM6O5ng+EEQAAzrPrr7/ePPKj559//rwfkwasAADAKsIIAACwijACAAiK3ExSh9D8fyaMAAACSocTV1kNOY7CI+nP/2f3/z03aMAKAAgoHRhMx+zQ6emVDqYV7GHZ09KxM1JSUsycMO5stghsWWmNiAYR/X/W/++0A8L5gzACAAi4qlWrmp9uIDnf9ItShzfXkVPPdxAqaJw8lpUGEff/O7cIIwCAgNMvtWrVqknlypX9nsgtEPSYOtFb586d83T7IBSk5qGsdP281Ii4CCMAgKDRL6pAfFnl5rg6B0xUVBRhpACUFTfSAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFUR/m6QmpoqH3/8sSQkJEiDBg2ke/fuUqRI1pnm6NGjsnz5cjl06JBUr15devToISVKlMjLeQMAgFCsGTlx4oS0b99e7r77bomPj5dbbrlFevbsKSkpKZlus27dOqlVq5ZMmzZNfvzxR3niiSekTp06snXr1kCcPwAACKWakUmTJsnBgwfl22+/lbJly8q+ffukadOm8uqrr8qoUaMy3Oaxxx4zAWbp0qXm+blz5+SSSy6RyZMny5tvvhmY3wIAAIRGzcicOXNk4MCBJoioCy64QPr06SOzZ8/OdJvIyEifWzJ6S6d48eLmAQAAkOOaEb0Vs3PnTmnUqJHP8saNG5v2IJl55pln5IYbbpDrr7/e1KJ88cUXUrJkSRk/fnym25w+fdo8XMePH/e0V9FHoLj7CuQ+CyvKivLi2rKP9yFlVdCuq5zuM8dh5OTJk+I4jqdWxKXPtS1JZs6ePSthYWGyZ88eKVOmjLm1U65cObM8q9tBEyZMSLdcQ090dLQE2ooVKwK+z8KKsqK8uLbs431IWRWU6yopKSmwYcQNAW4thevYsWNZ9owZPHiwNGnSRN59913zXANNbGysafy6aNGiDLcZO3asjBkzxvNcj1mzZk3Tc6d06dISyMSmha/nU7Ro0YDttzCirCgvri37eB9SVgXtukqbGfIcRrTtR+3atc2tGm/6XLv4ZkRrP7TXjXfjVq0lufzyy+WFF17I8lj6SEsLKRihIVj7LYwoK8qLa8s+3oeUVUG5rnK6P78asPbt29c0Yk1OTjbPExMTZeHChdK/f3/POuvXr5epU6eaf4eHh5tuvJs2bfLZjz6vX7++P4cGAACFlF9dex988EFZvHixdOnSRbp16yYfffSRqS3xrvnQdh1TpkzxLNMGrEOGDDHBpWXLlrJx40ZZvXq1GTgNAADArzBSqVIl+frrr2XWrFmmQepdd90lgwYNkqioKM86OqaId+PUAQMGyA8//GBCjI5REhcXJ6+//rpUq1aN0gcAAP4PB1+qVCkZMWJEpq9rI1N9eIuJicl0UDQAABDamCgPAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWRfi7wYkTJ+S9996ThIQEadCggQwaNEiioqKy3S4+Pl6WLl0qYWFhMmDAALMtAACAXzUjv//+u7Rq1UqmT59unk+ePFk6dOggJ0+ezHK7sWPHSqdOnWTv3r2SnJwsAwcOlLVr11L6AADAv5qRRx991NRsrF69WooXLy533XWXNGzYUKZOnSr33ntvhtssWLBAnnrqKRM+2rdvb5bdd999cvDgQYofAAD4VzOiweLaa681QUSVL19e4uLizPLMaFDp3r27J4ioyMhIqVWrFsUPAAByXjNy+vRp006kXr16Psv1+UcffZTpdl988YWpNfnkk09MjUrlypWlT58+EhMTk+Wx9OE6fvy4+ZmammoegeLuK5D7LKwoK8qLa8s+3oeUVUG7rnK6zzDHcZycrHjkyBFTEzJnzhy55pprPMtffPFF+de//mXagqSluy5SpIi0bNlSSpQoIT169JBt27aZ8KL70VCSkYcfflgmTJiQbvmsWbMkOjo6R78YAACwKykpSYYMGSLHjh2T0qVL571mRMOEthc5evSoz3J9XqpUqQy30fV1u5SUFNmyZYtERPz3cDfffLNpb5JZGNEGr2PGjPGpGalZs6a53ZPVL5ObxLZixQqJjY2VokWLBmy/hRFlRXlxbdnH+5CyKmjXlXtnIzs5DiPFihUzt2S2b9/us/zHH3+UJk2aZLpds2bNzHZuEFFt2rSRt99+29ScaGBJS9uU6CMtLaRghIZg7bcwoqwoL64t+3gfUlYF5brK6f78asCqjVdnz55tqlvUvn37ZNGiRWa5a/ny5TJu3DjPcx2HZMOGDXLq1CnPsk8//VRatGiRYRABAAChxa8wol1yK1WqZGo2brrpJmnXrp20bdtWRowY4Vln/fr1ph2J6/bbb5dGjRrJRRddJCNHjpQuXbqYMDJt2rTA/iYAAKDwjzOi7TW0lmPJkiWyZ88eUyOi7Ti0kapLn1esWNHn9s7HH38sK1euNLd4unbtKt26dZMyZcoE9jcBAAChMRy83v/p27dvpq/reCLeY4oovR2jDWP0AQAA4I2J8gAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYFWEvxscPHhQpk+fLgkJCdKgQQO5+eabpUyZMjna9ttvv5Vnn31W2rRpI7fffntuzhcAAIRyzci+ffvk4osvls8++0waN24s8+bNk7Zt28qxY8ey3fbkyZMyaNAgWb58uaxatSov5wwAAEI1jEycOFEqVKggixcvljvvvFNWrFghx48flylTpmS77ahRo+TKK6+UVq1a5eV8AQBAKIcRDSFXX321RET89+5OiRIlJC4uThYuXJjldrNmzZKvvvpKHn/88bydLQAACN02I8nJyeY2TZ06dXyW6/M5c+Zkut2OHTtk9OjR8umnn0pkZGSOjnX69GnzcGnti0pNTTWPQHH3Fch9FlaUFeXFtWUf70PKqqBdVzndZ47DyKlTpzy1Id5KlSrleS2tlJQU007kwQcflObNm+f0UDJp0iSZMGFCuuXa3iQ6OloCTW83gbIKBq4tyorryi7eg3bLKikpKbBhpGTJklKkSBE5evSoz/LExMRMe9MsW7ZMvvvuO9myZYv87W9/M8u++eYbCQ8PN8+feeYZqVixYrrtxo4dK2PGjPGpGalZs6Z0795dSpcuLYFMbFr4sbGxUrRo0YDttzCirCgvri37eB9SVgXtunLvbAQsjOgJNmrUSL7//nuf5Vu3bs201uPCCy+UV155xWdZfHy82ddll10mUVFRGW6nt3MyuqWj2wUjNARrv4URZUV5cW3Zx/uQsioo11VO9+fXOCODBw+Wl19+2dRcVKlSRX7++WdZsmSJvPDCC551FixYYG6n6Hq1atXy1Ii4PvjgAxNC0i4HAAChya/eNHfffbc0a9bMjDXSt29fadeunfTp08cnWGjNh/aeAQAACHjNSPHixc19pfXr18uePXtk3Lhx0rp1a591+vXrJw0bNsx0Hzo+ibYZAQAA8DuMqLCwMOnQoYN5ZOSiiy4yj8x07dqVkgcAAB5MlAcAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAArWrL0A7Dh7zpHNuxLl0IlkqVwqStrGlJfwImH8d8AarkkECmEEKACWbt0vExZuk/3Hkj3LqpWJkvFxTaVn82pWzw2hiWsSgcRtGqAAfOjfNvMrnyCiDhxLNsv1dYBrEgUZYQTI59XgWiPiZPCau0xf1/UArkkUVIQRIB/TNiJpa0S8aQTR13U9gGsSBRVhBMjHtLFqINcD8oprEsFAGAHyMe01E8j1gLzimkQwEEaAfEy772qvmcw68OpyfV3XA7gmUVARRoB8TMcR0e67Km0gcZ/r64w3Aq5JFGSEESCf03FEpg1tJVXL+N6K0ee6nHFGwDWJgo5Bz4ACQANHbNOqjMCKfINrEoFEGAEKCL0V065eBdunAXhwTSJQuE0DAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAMAqwggAALCKMAIAAKwijAAAAKsIIwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsivB3g19++UWmTp0qCQkJ0qBBAxk9erRUrVo10/VTU1Pl/fffl88++8z8u02bNjJ8+HCJiorK67kDAIBQqxnZuXOntG7dWn777Tfp16+ffPPNN+b5f/7zn0y3ad++vaxcuVIuvfRSueyyy+SVV16Rjh07SnJyciDOHwAAhFLNyMSJE6Vu3bry3nvvSVhYmAwcOFDq168vkydPlieeeCLDbRYtWiRVqlTxPO/evbvUqFFDlixZIgMGDMj7bwAAAEKnZmTp0qXSv39/E0RUsWLFJC4uzizPjHcQURUrVpSiRYvK8ePHc3vOAAAgFGtGkpKS5NChQ6ZWw1vNmjVl165dOT7giy++KEWKFJGuXbtmus7p06fNw+UGF21zoo9AcfcVyH0WVpQV5cW1ZR/vQ8qqoF1XOd1njsNISkqK+RkdHe2zXJ+7r2VnxYoVcv/998uUKVNMiMnMpEmTZMKECemWL1++PN3xA0HPC5RVMHBtUVZcV3bxHrRbVlqREdAwUrJkSQkPD5fExESf5YcPH5ayZctmu732ptFGr+PHj5eRI0dmue7YsWNlzJgxPjUjGl60vUnp0qUlkIlNCz82NtbcOgJlxbV1/vE+pKy4rgrvezCnTTJyHEYiIiKkWbNmpgeNt/j4eLnwwguz3HbNmjXSu3dvue++++SBBx7I9liRkZHmkZYWUjBCQ7D2WxhRVpQX15Z9vA8pq4JyXeV0f341YB02bJjMnj1bdu/ebZ5rMFm2bJlZ7po1a5YMGjTI83zdunXSq1cvuffee2XcuHH+HA4AAIQAv7r23nHHHbJ582Zp2bKltGjRwtSK3HjjjTJ48GDPOj/99JNP75q+ffua3jebNm2SPn36eJYPGTLEPAAAQGiL8Le65d///rds375d9uzZY8YYiYmJ8VlHA0a7du08z2fOnClnz55Nt6+GDRvm5bwBAECoDgevGjVqZB4Z0ZDhHTSuvPLK3J8dAAAo9HIVRgCEjrPnHNm8K1EOnUiWyqWipG1MeQkvEpbj1/PjORdEKWfOycwNu6WCiMzYsFuGtq8nxSKY6xSFA2EEQKaWbt0vExZuk/3H/jeXVLUyUTI+rqn0bF4t29dtyI/nlFeTlmyT19bukqJFHHmqrciTy7bLox//JDd3ipGxvZraPj0gz4jVADL9Ur9t5lc+X+rqwLFks1y/ILN6XbfPb+ds45zySsv5lTW75Jzju1yf63J9HSjoCCMAMrzNobULab7/DOfPh/6lntnrSrfX/eSXc7ZxToG4NaPlnBV9XdcDCjLCCIB0tL1F2tqFtLL6TteXdHvdT345ZxvnlFfaNiS77KSv63pAQUYYAZCONvzMT/sJ5LHO5znlVUJiUkDXA/IrwgiAdLQHSn7aTyCPdT7PKa9ql48O6HpAfkUYAZCOdoXVHihZdYbVnrKZva7LdXvdT345ZxvnlFc3tKtjyjkr+rquBxRkhBEA6eiYHNoVVqX9Lgz786HdSjN7Xen253Nsj+zO2cY55ZWOI+KWc2b0dcYbQUFHGAGQIR2TY9rQVlK1jO9tDX2uy3V8i6xetzGmR3bnXBDHGdFyHtE5Jl0NiT7X5YwzgsKAQc8AZEq/vGObVs10NNPsXrchP55TXmnguKt7Y5m5fqfIkW1yb49GjMCKQoUwAiBL+iXerl6FXL9uQ348p7zSWzHaNmTJkm3mZ1GGgkchwm0aAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWEUYAAIBVhBEAAGAVYQQAAFhFGAEAAFYRRgAAgFWEEQAAYBVhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYFeHvBvHx8fLss89KQkKCNGjQQO677z6pX79+wLcBgPzo7DlHNu9KlEMnkqVyqShpG1NewouEmddOpZyVx5dsk92Hk6ROhWi5v1dTKV4sPEfbZvWaSjlzTmZu2C0VRGTGht0ytH09KRaRs78n83Lc3O7XPWc914TEJKldPlpuaFcnIOcc7PNGPg8jW7dulY4dO8qwYcNMoHjnnXfk0ksvNWGjRo0aAdsGAPKjpVv3y4SF22T/sWTPsmplomR8XFOZ+9VeWbHtkGf52p9FZmzcI7FNK8trw9pkua3K7LWezavJpCXb5LW1u6RoEUeeaivy5LLt8ujHP8nNnWJkbK+muT7n7I6b2/16n/M553/bPLbkhzyfs+47u9fzsm/YEeY4jtelkrWBAwfKvn37ZN26deb52bNnpVGjRhIXFyfPPfdcwLZJ6/jx41KmTBk5duyYlC5dWgIlNTVVlixZIr169ZKiRYsGbL+FEWVFeYX6taVfYrfN/ErSfmDq39PZfYi2rFFavtt73K9t3b/TuzWt7Ak5keEaRs7KPZvD5fTZ/64xonPmX+65OWf3uNOGtsr0yzmr/aY954zk9pzVLZ1j5NU1uzJ93T3vjK6r7Pad1e9cmKUG8T2Y0+9vv9qMfPLJJyZEuMLDw6V3796ycuXKgG4DAPmJVuvrX9MZfYHn5K+5bzMIItlt6/z5yOpLXWntg94OCdQ5u6/ptroPf/cbzHN2/tzWCcJ5Z7Ut8tFtmpMnT8rhw4elevXqPsv1ubYFCdQ26vTp0+bh0kSlEhMTTYILFN1XUlKSOcf8/BdZfkBZUV6hfG1tSTgi/zl82P9GdgEWcc6RpKRzEpFaRM6e+18bh9dWfCOD2tYK6Dn/5/BJ+SR+p1xSu1xA9xuIcy6Sg/NuWb2kz3WVk31n9jsXdqlBfA+eOHHC/MzuJkyOryc3BERGRvosL168eKYBITfbqEmTJsmECRPSLY+Jicnp6QJAoTQkg2WjJouMCsKxekyWoAnWOef1vIP5O4eyEydOmNs1eQ4jpUqVMolJaye8aZKqUKFCwLZRY8eOlTFjxnienzt3zuxDtwkLCwvovayaNWvKr7/+GtC2KIURZUV5cW3Zx/uQsipo15XWiGgQSXuHJNdhRNt6tGzZUr744gu59dZbPcs3bdokF198ccC2cWtS0tamlC1bVoJFC58wQllxbdnF+5Cy4roqnO/BrGpEctWA9aabbpIPPvhAtm3bZp5rDxltoKrLXa+//rrExsb6tQ0AAAhdfrVB0toNHTdEazXq1KljGqHqLZV+/fp51tm7d6+pCfFnGwAAELr8CiPaXuOll16Shx9+2IQODRflyvm2Oh4+fLj06dPHr21s0VtB48ePT3dLCJQV1xbvw/yIzyzKqrBeV34NegYAABBoTJQHAACsIowAAACrCCMAAMAq26MbnxcbN26UN954Q3bs2CHVqlWTwYMH+zSyVb///rsZ+XXLli1SsWJFueWWW6R79+4SqnQ4/uuuu04OHjxoJlAqX76857VTp07J5MmTZdWqVWY03UGDBsnQoUMl1OgcSzqAn7e///3vMmLECJ8B+1555RX58MMPzSSROhHVP/7xD4mICIm3ng+9lnRyTO1tp+/Du+66K914Q7Nnz5aZM2fKH3/8IZ06dZJ77rlHSpQoIaGkS5cuPtNhuK655hq5++67Pc9Xr14tL7/8shw6dMiM53T//fdL1apVJdTMnTvXXDd6fenvr5/vffv29VlHe3Q+88wzsmvXLqlXr565rho3biyhZtWqVea7UAc3q1u3rnkPNmvWzGcd7Wjy+OOPm+E49H2qn1ft27cP/sk5hdz8+fOdbt26OW+88Ybz6aefOk899ZRTrFgx56WXXvKsk5yc7DRt2tTp0qWLs3DhQufRRx91wsPDnQ8//NAJVSNHjnSaN29u5qbav3+/z2u9e/d2Gjdu7MydO9f5v//7Pyc6Otp5+umnnVBTpUoV5+GHH3Y2bNjgefz6668+64wZM8apWLGi88477zjvvfeeU716deemm25yQs2OHTucqlWrOn379nU+/vhj875s166dc+TIEc86L7/8shMVFeW8+OKLzrx585wWLVo4l112mXPu3DknlGzatMnnmnrllVfM+1DLzLVixQonIiLCeeihh5xFixY5sbGxTr169ZwTJ044oUQ/xyMjI53nn3/eWbVqlfkc0nKZPn26Z53t27c7pUqVcoYPH+4sWbLEueGGG5yyZcs6v/zyixNKZsyY4RQtWtSZNGmSs3LlSueBBx5wSpcu7Wzbts2zTmJiolOjRg0nLi7OWbx4sfOvf/3LbPP5558H/fwKfRg5efJkumXDhg0zwcP16quvmgva+4NRvzD0yzgU6RdBs2bNzIdf2jCyZs0as+zrr7/2LNMPAH2zJyUlOaEWRvQNnhktNw21GkJcGnbDwsKcn3/+2QklPXr0cDp06OCcPXvWsyw1NdU5c+aM598VKlRwHnnkEc/rP/zwg7nWli5d6oSyUaNGmSCnZeRq27atc/3113ueawgpWbKkM2XKFCeUaAgbNGiQzzINvPpl6v1537p1a89zDbdNmjRxbr31VieUtGjRIt3vrGU1cOBAz/OJEyc6lSpVck6fPu3zx6f+QR9shb7NSHR0tM9zHSP/66+/lgsvvNCzTEeE1Sph7yHntZpPq/a06i+UaPXdyJEj5d1335WoqKh0r2tZ1ahRQy666CKfstJy3bx5s4SaF154QS677DIZNmyYuZ3lTavR9daM9y3BHj16SLFixUw5hgq9jbBs2TJz+6pIkf995OitKp0yQsXHx5tbXnFxcZ7XtRq9YcOGsnLlSglVertG34t6+8+9tafvNb3V5V1WJUuWlCuuuCLkyqp169by7bffmtt67gzv+rndtm1bzzr6Xks79pWWXaiV1eHDh9PND3PBBRfI8uXLfcrK/Yzy/nz/7LPPspzcNhAKfRjxvr9/ySWXmMLX+19PP/205zUdFTbtf5L7XF8LFfrFOWTIEHNf2juseaOs/qd+/frmS+LBBx8096GvvvpqeeKJJ3zKSud50C8Kl04cWalSpZC6rn744QfzU+8/6zQQl19+ufn5zTffeNZxyyOj92EolVVa8+bNk6NHj/pMn7Fnzx4z+RhlJfLII4/IlVdeaSZ50/ZHtWvXNm3YtP2M+5m2b98+ykr+2xbpvffe80xcq+WyYMECOXLkiJkoL6vPdw0i+/fvD+q1HhFKF60W+IYNG+Sxxx6TNm3aeN7gWtBpR57Thpnua6HCHYHPe8bktDIqK7cGJZTKyv0rwi2Lbt26mcmgtGGc1ixpCMmorNxrK5TKKjk52fzU99u9995rAu/8+fPNX7Xr168370W3PDJ6H4ZSWaU1ffp0U+OhYddFWfk2XtX50MaNGyetWrUyNUb6+a41I1dddVWWZXXmzBkT6gI5E3x+9txzz5nGvRrYYmJiTI2l1hjpNaZlYfu7MGTCiF6oSqvUU1JSzJeGG0a0p4ibFl1uL4kKFSpIqHj77bfNba127dp5qjzdWqW//e1vplW1lpVWqXtzyy6UykqlfdNqINE3rNYE/OUvfzFlpX91pP3A02srlMrK7YmlPdQ0qKmuXbua2bunTZtmwoi7jl5L3rOGalm1aNFCQpH2/Pj000/NX7PevMvKW6hdV0p7g+j8Z/pTaa3bb7/9Zmp3NYzoH0r6ZZpRWWk5hkoQUVWqVDHX04EDB8yjQYMGpjeWfua7TRSy+i707lEZDCFzm8abVhdr1aebBjWofPnllz7r6AdlqVKlfP4iKew++ugjE0imTJliHm4XVa1VcrvKaVn9/PPPnqDilpXybkcSivQNrtyuqFpWeo15h7edO3eaN3vaLq2FmXYd1C+FtNW/+j7UsOZeO9qexHuSzaSkJM8km6FIu2BquOjfv7/Pcr0locMPeJeV0jZboVZWev3orXdvep15f6G6NSbe9DMr1MrKpd2f9f2mn1MLFy40fxi4bbkyK6vzMqecU8i9/vrrTkJCguf5b7/95lx44YWmdb9LuzZpr4c333zTPD906JATExPj3H777U4o0y6YaXvTHD161Clfvrxz7733erpFd+zY0enatasTStatW+fTy+PAgQNOmzZtTIt1tyuq/tQeWQMGDPD0ItFuhXptpaSkOKHk5ptvdjp16uT88ccf5vn3339ven9MnTrVs85VV11lytDtAafdprWXlr4fQ41eL9rFUruGZ+See+4xr+vnmXr//fdNL60tW7Y4oaRnz57OJZdc4ukJefjwYfMe7Nevn2cdHdahRIkSTnx8vKfrtPae9O7lFgq+/PJLZ/369Z7n2oVeu+3qcu/eknodaRdotWvXLtO7xruXW7AU+jCi/fG1m6p+AehYIjqOwbXXXmu+PLzpOBD64ah99YsXL+706tXL88EZqjIKI0rHa9FurTVr1jT99Vu1apVufI3Cbu/evU7//v1Nd1QNHHpd6TWze/dun/W0e2qjRo3MG1q7Z9apU8f56quvnFCjXU+1i2C5cuXMl4WOTXPHHXf4dPU9ePCgGXtExz6oXbu2GZ/F/VAMNfp763vPewwIb6dOnTIhV6+7+vXrm/LUMX9Cjf6hqX8MaWht2bKl+Qy//PLLnX379nnW0T8KRo8ebcaXatiwofnp/jEVSvbv32/G7dHrRT+7a9WqZT7j03r22WfNd6CWlYY27UJ+Pv54CplZe7UFunb/qlWrlk/vBm9aLbx9+3ZTNarrhTq9FaPtH7QXkvYC8ea2jdD7sXrvMVTp7T7tDq3XizZgzYi+xX788UczGmuTJk18ureGGh3dUavQtQGd3gbNiN7K0veqlpV3F8NQoteUNjDU915WtEeErqfvwcw+10KBDsGgt0n11l/lypUzXEdH2dbvAW3AGWpta9K2RdIu440aNcq0zYx29tD3od7S0TI9H0ImjAAAgPwpdP9EAwAA+QJhBAAAWEUYAQAAVhFGAACAVYQRAABgFWEEAABYRRgBAABWhcxEeUAomD17thlcLTM64JM7EWJh9P7775vJ0nRSMAAFB2EEKEQ+/PBDOXv2rPl3QkKCbNy4Ua699lrPqK/t27cv1GFEp0hfsWIFYQQoYAgjQCHy7rvvev49c+ZME0beeecdM2uu97QHGzZsMENCt2zZUmrUqOF5TWtVtHalW7ducvLkSTNrbqVKlaRt27bmdZ0uQYe21+HHmzZtmuF2OpT7999/bwJB69at051jTo+vM7Lq8fU4OnS1ziqt22qw0m101lWdjsA1f/5883P16tVm6G8dbv6KK64wAa13794+w8/PnTtXLr30UjPja1bHVDolvc7qrcP966ymmQ1jDyAPgj77DQArZsyYYSZb00nVXJ988olTuXJlMyGdTuynEx3q7LguXVe36dy5s5koSye300nYhg4daia2a9y4sdlOJ9J6+umn022ns6jqJHf6Uye8GzRokGcWY3+Or6/VrVvXufrqq52FCxea10aOHOkMHDjQueaaa8zklzrp4NatW31mBtZtdTIwXe+f//ynmcBRl+mEhd50Ftf58+dne8zx48ebc9TfRydk03PX3wFAYBFGgBAJIzq9un6xzpkzx7POjh07zEyna9eu9fli7tu3r3PmzBmzTL+0ddl1113nmWVX960hJTU11Wc7ncHZne16+/btJrTMnj3b7+NfccUVzunTp7P8/W677TYTErzptjpTt8ufMJL2mPPmzTMzLXvPSK0z41avXj3bcwPgH27TACFiwYIFnsatc+bMMT/1+1tnHF61apV07NjRs+5NN90k4eHh5t9uG5Phw4d72p7oMr1lorPw1qlTx7PdyJEjpUSJEubfDRs2lP79+5tbINpuxZ/jjxgxIsMZe3fs2CE//fSTmVVUb5ts3rw5YOWT9phvvvmmNG/e3Nzq+vMPN/O63rbR21UtWrQI2LGBUEcYAULE7t27TZj44IMPfJbrl6p3uw1Vrlw5z78jIyMzXZacnOyznXcwUTExMbJ8+XK/j5922nJtlKuNU5ctW2bar+i5aLuQxMRE85obnPIi7TH1fDWApD3fgQMHZjr1OoDcIYwAIaJ06dLmS1S7vwaLNgBN+7xixYp+Hz/tl/3ixYtNqNGaEW1QqzQkaI3Kf+/OZMytyfHu7qz/TklJyfaYer5au/PGG29ke74A8oZBz4AQ0aNHDxMO9LaJN63dOHz4cECOobdiXPqFv2jRIunQoUOej3/gwAEpX768J4iotDUWSm8RedfWVK5cWSIiIkyIcX3++eeSmpqa7e/Ss2dPmTdvnhw6dMhn+b59+7LdFoB/qBkBQoTeDhk3bpwMGzbMdK3V9hC//PKL6eaq3YArVKiQ52PobRRtW6K3UmbNmmWW/eMf/8jz8WNjY+XOO++UG2+80YQbHUtEj5WWdiV+4YUX5OjRo+ZWjnbpHTJkiIwePdqECG1r8tZbb+Xots6YMWNMjUybNm1MW5iyZcvKli1bzLl/9913eSonAL6oGQEKKW2/oe0bvL94J06caL7I9ZbEunXrzJgZK1eu9IwHouvqNt41ENpoU5dpzYQrOjraLNNbGd40gGjo0IalGho2bdrks05uju+2PdGGpNpode3atSYgLF261Kzr3opxx1nRIKRBRY+jXnvtNRNGdKwQHdtElw8dOtTTTiWzY+rvuGbNGnnkkUdk586dEh8fL3/5y1/MfgAEVph2qQnwPgGEGL01ogOQaVDw7hUDADlBzQgAALCKMAIgzzK71QEAOcFtGgAAYBU1IwAAwCrCCAAAsIowAgAArCKMAAAAqwgjAADAKsIIAACwijACAACsIowAAACrCCMAAEBs+n+tfZjksK4FeQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n", "data_pred['Frequency'] = logmodel.predict(data_pred)\n", "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing and plotting uncertainty" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "trusted": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGiCAYAAADEJZ3cAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHb9JREFUeJzt3QeQFfXhwPHfKYKIFAVRkCaKEmMJYpkoIGJAoo5iRGHsCYplRCxEBkf8jxpDBjGggmIDEdSARhEFWwwaKyoqNsASCwEHFJDeef/5beaed8DBHZy5n97nM/Pm2L3du2Xn7Xvf2/YKcrlcLgAAJGy7il4AAIAtESwAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAD+vYFmzZk145JFHQocOHUKdOnXCuHHjSjXfE088EQ4//PCw++67hzZt2oQXX3xxa5cXAKiEyhQsQ4cODWPHjg3XXHNNWLRoUVi9evUW53nppZdC165dw1lnnRXeeOONLHY6d+4cPvzww21ZbgCgEikoy2cJxUkLCgr+O2NBQRg9enQWIptz/PHHh+222y489dRT+XGtWrXKHiNGjNiWZQcAKoky7WEpjJWyeOWVV8Kxxx5bbFzHjh2z8QAApVEl/IgWL14clixZkp27UlT9+vXDnDlzSpxv1apV2aPQ+vXrw4IFC0LdunW3KpoAgP+9eGQmdkDDhg2zoy3JBkuhDReySpUq2X+iJAMGDAjXX3/9/2DJAIAf26xZs0KjRo3SDZaaNWuG6tWrh++++67Y+Hnz5mV7WUrSr1+/cOWVV+aH4wm+TZo0yf7DtWrV+jEXGQAoxyMtjRs3znpgW/2owRIP3xx22GHh5ZdfDpdeemmxK4eOOOKIEuerVq1a9thQjBXBAgA/LeVxOke53zjuL3/5S2jatGl+uHfv3mH8+PFhwoQJYe3atWHkyJFhypQpoVevXuX9qwGAn6kyBUvcMxJvGBcf0YUXXpj9+7LLLstPs3LlyuwQTqHf/e53YdCgQaFHjx5hxx13DP379w+jRo0KRx11VHn+PwCAn7Ey3Ycl7iFZunTpRuPj4Zt4rkphsMQrfGrXrr3RdCtWrMhPV9ZjYPHnxRBySAgAfhrK8/27TOewxKt7CveulCTuRYmPTdmaWAEA8OGHAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsA8PMLlpEjR4bmzZuHKlWqhJYtW4bx48dvdvrly5eHXr16hYYNG4YddtghNGnSJPTr1y+sXbt2W5YbAKhEyhQskyZNCj179gx/+tOfwsKFC8Mll1wSTjvttPDWW2+VOE/fvn2zqHn66afDihUrwkMPPRTuvPPOMHDgwPJYfgCgEijI5XK50k587LHHhjp16oS///3v+XFHHHFEaNGiRRgzZswm52nbtm3Yb7/9wr333psfd8IJJ4QaNWqEcePGler3Ll68ONSuXTssWrQo1KpVq7SLCwBUoPJ8/y71HpbYNW+88UY4+uiji43v0KFDeO2110qc78wzzwzPPPNMthdm2bJl4Z///Gd4/fXXwxlnnLFNCw4AVB5VSjvhkiVLsvNRdtttt2Lj69evH+bOnVvifBdddFGYOXNmOPzww7Ph7bbbLtx0002hS5cuJc6zatWq7FG00ACAymubrxJav359KCgoKPH7ffr0yc5hefvtt7MImTx5chg0aFAYPHhwifMMGDAg24VU+GjcuPG2LiYAUBmCpWbNmtl5J/PmzSs2/ttvvw177LFHiTEzbNiwcMUVV4TWrVuHqlWrhnbt2mUn7t56660l/q54FVE83lX4mDVrVln+TwBAZQ2WuBflyCOPzPaQFPXCCy9k44tGyrp16/LzxMufNzyvN04TL3EuSbVq1bKTc4o+AIDKq0yHhOLhnYkTJ4YRI0Zke1ZuvvnmMG3atGwPSqEbbrgh1K1bNx8sp556arjlllvCq6++mp0H8/zzz4e77747Gw8AUK4n3UadOnUKo0aNCjfeeGN2M7h4OXM8P6VVq1b5aeJJtXGvSqGhQ4eG66+/Ppx77rnZybl77rln6N27d3bYBwCg3O/DUlHchwUAfnoq5D4sAAAVRbAAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAD/PYFm/fn1YuHBhyOVyZZpv1apVYcWKFVvzKwGASqzMwXLzzTeHunXrhoYNG4bdd989jBgxYovzzJgxI3Tq1CnUrl07NGrUKHTr1i3Mnz9/a5cZAKhkyhQs48aNC/379w9jx44Ny5cvD0OGDAkXXHBBeOmll0qcZ86cOaFt27ahcePG4dtvvw3fffdd6N69e3j77bfLY/kBgEqgIFeG4zoxPOIekocffrjYuAYNGmQxsymXXnppmDhxYvjkk0/CDjvssFULuXjx4mzvzKJFi0KtWrW26mcAAP9b5fn+vV1ZzluJe0XatGlTbHy7du3Cm2++WeJ8MVZOPvnkLFbiAsefAwDwowTL0qVLw8qVK0O9evWKjd9tt93CvHnzSpxv1qxZoaCgILRq1So7LLTzzjuHM888MyxYsGCzJ+fGKiv6AAAqr1IHS4yOaO3atcXGx+Htt99+s/MOHz483HrrrVl4xBNw33nnnXDxxReXOP2AAQOyXUiFjxg6AEDlVepgqVmzZnb8ae7cucXGx+F4xVBJ9txzz3D88cdnh46iJk2ahEsuuSRMmjSpxMui+/Xrlx0+KnzEvTQAQOVVpquE4gm2zz//fLFxzz77bD5GCg8dFT1EdMwxx4Rly5YVmydOU7169fxemw1Vq1Yti6OiDwCg8ipTsMQ9H//4xz/CwIEDw/Tp00Pfvn3DZ599Fq666qr8NIMGDQr77rtvfvjqq68Or7zySrjtttvCv//97/Dkk0+GwYMHZ5dDAwCUe7AcddRRYcKECVl0dO7cOUydOjXb49KyZcv8NPGk2nhDuUL7779/FjnxaqEOHTqEP//5z+Haa68NN9xwQ1l+NQBQiZXpPiwVxX1YAOCnp0LuwwIAUFEECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwDw8wyWpUuXhk8//TSsWLGiTPOtWbMmfPjhh2HWrFlb82sBgEqqzMHSt2/fUK9evdC+fftQt27dcPPNN5dp3gMPPDBcccUVZf21AEAlVqZgGTFiRBg2bFh49dVXw+zZs8P48eNDv379wtNPP73FeSdNmhSeffbZ0KFDh21ZXgCgEipTsNx1112ha9euoXXr1tlwp06dsj0tcfzmzJkzJ/Ts2TOMGTMmVK9efduWGACodEodLOvXrw/vvfdeOOKII4qNP/LII8PUqVM3O99ZZ50VevfuHVq1arVtSwsAVEpVSjvhkiVLwurVq7PzVoqKw999912J8910003Z16uuuqrUC7Vq1arsUWjx4sWlnhcAqMR7WKpU+W/bxGgpKobFDjvssMl5pk2bFgYMGBCuvvrq8PHHH2dXCMXwiQES/x2vGtqUOE/t2rXzj8aNG5ftfwUAVM49LDVq1Ai77LJL+Oabb4qNj8MlBcX8+fND8+bNQ58+ffLjvv7661BQUBC6d+8ennvuudCwYcON5osn8l555ZX54Rg4ogUAKq+CXC6XK+3Ep556aliwYEGYPHlyNhxnbdmyZejYsWMYOnRoNm7evHlZqPziF7/Y5M848cQTw4477hgeffTRUi9kDJa4p2XRokWhVq1apZ4PAKg45fn+XaarhK699trw+uuvZ4d4Xn755XDBBReEuXPnFtuDcscdd4Rf//rX27RQAABbHSzxKp+4dyXe5fbyyy/P7ngbw6VZs2b5aerXrx/233//En9G06ZNQ5MmTcryawGASq5Mh4QqikNCAPDTU2GHhAAAKoJgAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQB+nsEye/bs8Prrr4d58+aVavp169aFGTNmhA8++CCsWLFia34lAFCJlSlY1q9fH3r06BH22WefcNFFF4UmTZqEvn37bnaeoUOHhmbNmoUuXbqE7t27h4YNG4b77rtvW5cbAKhEqpRl4mHDhoXHHnssTJs2Ley7775hypQpoW3btuGwww4LXbt23eQ8ixcvDm+99VbYY489suGRI0eG888/Pxx66KHh4IMPLp//BQDws1aQy+VypZ34kEMOyULj7rvvzo874YQTsq8TJ04s1c+Iv27HHXcMt99+e+jZs2ep5onRU7t27bBo0aJQq1at0i4uAFCByvP9u9SHhOJ5KPEclNatWxcbH4fffffdUv/C9957L6xevTq0aNGixGlWrVqV/SeLPgCAyqvUwbJkyZKwdu3asOuuuxYbX69evbBw4cJS/Yxly5aF3//+96Fdu3ahffv2JU43YMCArMgKH40bNy7tYgIAlTlYqlatmn3d8Cqf5cuX57+3OStXrsxOvI17Vx555JFQUFBQ4rT9+vXLdh8VPmbNmlXaxQQAKvNJtzvttFO2NyVe0lxUHG7atOlm542HeGKsxPB48cUXQ/369Tc7fbVq1bIHAECZL2vu2LFjmDBhQn44HiJ66qmnsvGFvvzyy/DKK69sFCtx/OTJk/NXCwEA/CiXNV933XXZJczxXiwnnnhiGDNmTHZIqE+fPvlp7r///jBkyJDw/fffZ8Onn356FjD33HNPmDlzZvaI4r1Z4gMAoFyDpWXLltm9VwYPHhyGDx+eXekThxs0aJCfJkZImzZt8sPxnJd4JVGcvqjzzjsvewAAlOt9WCqK+7AAwE9PhdyHBQCgoggWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5VbZmpvfffz989dVXoUWLFqFly5Y/2jwAAGXew7J69erQpUuXcMwxx4QhQ4aEww8/PPTo0SPkcrlynQcAYKv3sMTgePXVV8O0adNCo0aNwkcffRQOPfTQ0L59+3D22WeX2zwAAFu9h2X06NGhW7duWXhEv/zlL8Nvf/vbbHx5zgMAsFV7WNauXRumT58eevXqVWz8QQcdFIYPH15u80SrVq3KHoUWLVqUfV28eHFpFxcAqGCF79vlcRpIqYNl6dKlYd26dWGXXXYpNr5u3brh+++/L7d5ogEDBoTrr79+o/GNGzcu7eICAImYP39+qF279v8mWKpVq5Z9Xb58+UZRsuOOO5bbPFG/fv3ClVdemR+OcdO0adPw9ddfb/N/uLKXboy+WbNmhVq1alX04vykWZfWZWo8J63LFMUjJE2aNAm77rrrNv+sUgdL9erVwx577JFFQ1FxuHnz5uU2T2HoFMZOUTFWvNFuu7gOrcfyYV2WH+vSekyN52T52W67bb/tW5l+QjxZ9rHHHgvr16/PhleuXBmefPLJbHyhjz/+OEyYMKFM8wAAlFuwXHfddeE///lPOPXUU8M999wTTjjhhFClSpVih2/GjRsXzjnnnDLNAwBQbsHSrFmz8M4772R3qn3xxRdD27Ztw1tvvZWdRFto//33DyeffHKZ5tmSeHjo//7v/zZ5mIjSsx7Lj3VpXabGc9K6/Lk/LwtybjkLACTOhx8CAMkTLABA8gQLAPDz+vDDH9vChQvD22+/HQoKCrLb99evX3+jaeIt++OHKcabz8VPfo73eaFkzz33XFiwYEF2ldYOO+xQ7Htz584NU6ZMCTVq1Aht2rRxUvMGPv300zB16tSN1mn37t03Gvf++++Hzz//POy1117hV7/6ladkCebMmZOt0wYNGmQfgmr7LpuZM2eGd999d5Pfixc7xHtfFbJ9b1m8E3t8z/nmm2/C7rvvHg477LDsKtYN2b637LvvvgvvvfdeWL16dTjyyCNDnTp1yv/9O5eIa665Jte4ceNc586dc23bts1Vr149d9NNNxWb5tNPP801a9Yst99+++XatWuX22mnnXL33XdfhS1z6p566qlctWrV4gc45BYuXFjse6NGjcrWX1yPLVu2zDVt2jQ3c+bMClvWFN1+++25mjVr5rp161bsUdTatWuzcXXq1Ml16tQpt+uuu+ZOOeWU3OrVqytsuVO0bt263BVXXJE954477rjc0UcfnX1duXJlfhrb95Y9+eSTGz0f42vizjvvnFu2bFl+Otv3ln3++ee5ffbZJ9e8efPcySefnNt7772zdTljxoz8NLbv0rnjjjtyNWrUyLVv3z571K5dOzdhwoRi05TH9p1MsDz44IPFXuQfeuih7I02PqkKxRe5+KYQn0TRnXfematatWruq6++qpBlTtns2bNzjRo1yt14440bBcusWbOykInrL4rrM67XNm3aVOASpxkscePanLgOa9Wqlfvss8+y4S+++CKLlyFDhvyPlvKnIf7xEV/EPvroo/y4F154IbdkyZL8sO1760KwSZMmufPPPz8/zvZdOuecc07uoIMOyq1Zsyb/Oti6devc6aefnp/G9r1l8Q/d7bffPjd69Oj8uAceeCDb3ufPn1+u23cywbKhqVOnZm+006ZNy2+EcXjixIn5aWLgxDeHQYMGVeCSpvkiFis3vmk+/vjjGwXL4MGDszfZVatW5cdNmjQpmy6+4fJDsMS/COK6ee6553LffPPNRqvmqKOOyp199tnFxv3hD3/IHXrooVZjke10l112yV133XUlrhPb99Z55plnsu12ypQptu8yinunjj/++GLj4t7RLl265Idt31s2ZsyY7DlYdG/pihUrsnEjRowo1+17u9SOz/7tb38Lt912W3a33N69e2fnskQffPBB9vWAAw7ITx/Pydhvv/3y3+O//vSnP4WqVauGyy67bJOrJK6vuN7iNIUOPPDA7OuHH35oNRbx7bffhr/+9a/ZjY/iB3D2799/o3VZ9DlZuC49J4uvo3h+WqdOnbJ/P/HEExutH9v31rnvvvuy18h4PoDtu2ziNv3ZZ5+Fyy+/PIwaNSq7+3p8/Yuvn0XXpe178xo2bJh9nTFjRrH38mjatGnlun1XSe0kx8cffzzMnj07LFmyJLRq1arYJz5GG37iY7xjbvw0Z/7rlVdeCXfccUd2Yl48eXlT4rrc1HqMrMsfHHHEEeHLL78M9erVy4afffbZcPzxx2dvEKeddlrcO5l9Qu6m1mU8uSx+btbmPpW8spg3b1729fbbb89ewPbZZ5/w2muvhdatW2fxEk8UtX2X3fz587P1d8sttxQbb/suncKTbONn333xxRdZrMThOD6yfZfO0UcfHdq3bx9OOeWUcMUVV2Tr7a677speNwu36/LavpPaw3LiiSeGsWPHZm+6d955ZzjvvPOyF7ao8La+8eziouKwN4UfnH/++aFz587hpZdeyvZWxXUZxQ+gLCzZuC43tR4j6/IH8cWrMFai4447LruaKn54ZxSDMO6l2tS6LPwePzynYsB99NFH2fqLH5Iao3rgwIG27600evTo7BNwzzrrrGLjbd+lE99fvvrqq2zPQAy/uFcgXuly9tln277LID4H49Wo11xzTbZ9x/UYn5u77bZb2Hnnncv1/TupYCkq/iUba6zwDXfvvffOvn799dfFpovDzZs3r5BlTFG7du2yN4bx48dnj3jZcjRx4sQwffr0/LrccD3GDTeyLrf8cfPxMFGhktZl/Ayt8vg49Z+Dwm03/gVWuE7iX7Ex/govG7d9b93hoLinb8PLR23fpRP/qOvSpUv+D4t4OXMcnjx58hbXpe27uHh4J/6xPHz48DBs2LDs9g6ffPJJ/hYP5bV9J/GKGncXrVixoti4uIsuHvdu3Lhx/ryA+O9HHnkkP018M4677OMnQPNfd999d7ZnpfBx1VVX5V/cTj/99HwMzpo1Kx8zUdyzFY9FuofID+K9GTY8tBEDOu55KRTXZfzrLN57IFqzZk12WNNz8gdxuz344IOLHeNev3599oJm+946b775ZnYI44ILLtjoe7bv0onPvcI/4grF4UaNGhVbl7bvLYt7pooaMmRIqF27dujatWv5vn/nEjB9+vTcAQcckOvfv39u5MiR2aW48Z4s8Vrtopc6P/roo7kqVark+vTpk13pEqfZ8L4YFLepq4SiM844I7vsOa7HP/7xj9l6HTt2rNVXRMeOHXM9evTIDR8+PDdw4MDcXnvtlV0GuWDBgvw0c+fOzZ6Hv/nNb7J7EcR7izRo0CC7rJwfvPTSS9k9beJzLd574aSTTsrVrVu32FVptu/S69mz52Yvubd9b9nDDz+cve5dfPHF2dUsvXr1yobvv/9+23cZxSsle/funb1/x9fMeI+VeHVlUeWxfSfzac3xrozxTO34V1g8ESfuLj7ppJM2OnH0jTfeCA8++GB27Ktt27bh3HPPDdtvv32FLXfq4l0cBw0aFEaMGBF22mmnYn/hPvDAA+Ff//pXNv6MM87I7k5I8btgxr1U8c6McbfxIYcckq2nDe+EGfe8xHOuCu90e/HFF7sD8ybE49v3339/drJoixYtsr0DRc8RimzfpXPhhReGDh06hG7dum3y+7bv0ol3Zh03bly2NzXedTXuEYgng9u+yybuWb733nuzPX9xT0o8P2hTh3q2dftOJlgAAJI+hwUAYHMECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAAAhdf8PFNxL2/4KN+gAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#sns.set(color_codes=True)\n", "plt.xlim(30,90)\n", "plt.ylim(0,1)\n", "#sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "trusted": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAG1CAYAAAARLUsBAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVXpJREFUeJzt3Qd4U1UbB/B/6WIXyl6VLTJlL0GG7A3KHrIUBD5QRIYoAkIZDmQpQ0BApoCg7L13BWQIyN6z0EJpgTbf855r0qR0pG3am/H/Pc99ktzcJKe3SfP2nPe8x81gMBhARERE5EJS6N0AIiIiouTGAIiIiIhcDgMgIiIicjkMgIiIiMjlMAAiIiIil8MAiIiIiFwOAyAiIiJyOQyAiIiIyOV4wA4cOHAA//zzD+rXr4/s2bPHefzLly+xe/du3LlzByVKlECxYsWSpZ1ERETkHHQNgDZs2IAhQ4YgPDwcJ0+exPbt2+MMgB4+fIg6deqoy+LFi2Pnzp3o1q0bJk2alGztJiIiIsemawAkgc+8efOQOXNm5MmTx6rHDBs2DM+ePcPff/+NtGnT4uDBg6hcuTIaNGiAevXqJXmbiYiIyPHpmgPUqFEjvPnmm1YfL8uWLVmyRPX4SPAjKlasqLZFixYlYUuJiIjImdhFDpC1rl27hsePH7+S8yNDYQEBATE+LiwsTG1GERERaggtU6ZMcHNzS9I2ExERkW1IR0hwcDBy5syJFClSuE4AJMGPyJAhg8V+X19f033R8ff3x8iRI5O8fURERJQ8HSK5c+d2nQAoVapU6vLp06cW+yUaNN4XnaFDh+KTTz4x3ZZgyc/PT53A9OnTJ2GLiYiIyFaCgoJUznC6dOkS/VwOFQBJ0OLp6YnLly9b7JfbBQoUiPFx3t7eaotKgh8GQERERI7FFukrdl8Icd++fVi9erW67uXlpabAL1261HS/1ALatm0bGjdurGMriYiIyJHo2gP077//Ys+ePQgMDDTVBZLeHJkZZpwdNmfOHFUosVmzZur2+PHjUbVqVbRu3RqVKlXC3Llz1bFdunTR80chIiIiB6JrAHT79m3s2LFDXZcARm7LJknOxgBIgh3z4ogy4+v48eMqMDp79iw+/PBD9OjRQw2NEREREVnDzSBzylwwicrHx0clQzMHiIjI8Ugh3RcvXujdDLIx6cxwd3dPlu9vh0qCJiIi1yb/s8tIwaNHj/RuCiURGQWSkZ+krtPHAIiIiByGMfjJmjUrUqdOzWK2ThbchoSE4O7du+p2jhw5kvT1GAAREZHDDHsZgx+p5E/OJ9V/Nf0kCJLfc2zDYU4/DZ6IiEgYc36k54ecV+r/fr9JnePFAIiIiBwK13B0bm7JtEYnAyAiIiI7IbObLly4AGfy9OlTnDt3DvaGARAREVEykWGdK1eu4ObNm9Hev2rVKrz99ttO9fvYvXs3SpYsCXvDAIiIiCiJPXnyBB999JFK3q5WrRpKlSqlVjOfPn06z71OGAAREZFLCY8wYP+FB1h97Ia6lNtJ3etTr1497Nq1Sy3tdPXqVdy7dw8zZszA0KFDMXz48Ggf9/z5c9VbFBEREe398hz379+P9j55jLyOFA6MbZgtKCgI58+fV8NU//zzzyuJx3Jb9oeGhlr13EbBwcEx9nLZCwZARETkMjacvIW3xm9Du1kH0H/JMXUpt2V/UpE1K/fv34+FCxeiaNGipv2NGjXCyJEjMW7cOIu8Hwk6evbsqYoBSk9Rvnz5cPDgQdP9p0+fRokSJVC4cGF1Wbp0aRw7dsx0/+LFi1XvkqyX6efnp3qcLl68aDHMVqVKFfTq1Uvd36RJExVMVahQAStXrrRo+4oVK1CuXDlVgsCa5xZffvml6umSx8nP8Mcff8AeMQAiIiKXIEFO74UBuPU4sjdD3H4cqvYnVRD0+++/o2zZsqY1Ls3JWpYSXJgHCVID59mzZ+ry4cOHqF+/Ptq1a2fqnfnss8/U80nvz61bt7BgwQKcOHFC3Se9TL1798by5ctVD4wcI2totm3bVhUaNH8NT09P9fzSw5M3b160aNECv/76q0X75LbsT5MmjVXPvX79ekyYMAHbtm1Tx+zcuRO//fYb7BEDICIicnoyzDXyj9OIbrDLuE/uT4rhMBnGKlCgQLT3pU2bFtmyZcPly5ct9k+cOBEeHh5IkSIFxo8fj+vXr2Pz5s3qPhl6ypw5s6lIoAQhnTt3Vtd/+OEHNGjQQD2nDG1J74wEKIcPH1bPYZQiRQr4+/urS6OOHTtiw4YNePDggbotlxs3blT7rX3umTNnolWrVnjrrbfU7ddffx19+vSBPWIlaCIicnqHLj18pefHnIQ9cr8cV7mAbatMS0+LLPEQE7nPy8vLdDtjxowWy0DI2lgy7PTvv/+ahpjatGmjAqI6deqgcePGqFGjhrrv1KlTqlfnr7/+sngNCUSkxyZPnjzqdtasWVXwZa5WrVoqsJIeHhkeW7ZsmRrKeuedd6x+bmmj9FaZkwDNHjEAIiIip3c3ONSmx8WH5Ons3btXDRNFLfInPT+SMCzHGJknHJsHScZlIiQgkeElGZLaunWr6nHp0KEDJk+erIIt6bH57rvvYm2TezRLTMg+6dGRYS8JgORSbhuPtea5pY1R2x9b8KcnDoEREZHTy5oupU2Pi4/u3bvj0qVLWLRo0Sv3ff3116qXpXnz5qZ9kv8jw0pGZ86cUTk7xhwiyQXy9vZWvT+SQC3bkiVL1H0y9LRmzRo1g8yctctKdOjQQQVrO3bswL59+0zDX9Y+tyRkS96POXkue8QeICIicnoV8vkih09KlfAcXZaP9Mtk90mpjrM1GZ6Sqe4ys+vatWtqSrz0ivzyyy9qVpUkSadLl850vPS4dOnSBd9//73KAxo4cCDq1q2L8uXLm3qAmjZtiqpVq6rgQ2aXVa5cWd03bNgwNctLhsUGDRqE9OnTq2BKXss8qIqJJFfLkFanTp3Updw2sua5P/30U9WbJZetW7dWvVRRE6vtBXuAiIjI6bmncMOIJtoU9KgrTRlvy/1yXFIYPXq0CnQkf6Zr164YMGCACm7ktvTkGPn4+KieFkl8njp1qgp+JNAx9vAIyc25ffs2Pv74Y1VHSO6fP3++uk/ycI4ePaqm23/++efqGMnLMZ/eLq9RsGDBGNsqM71k1pdcmrPmuQsVKqRyk/7++2+V/CwzzKRtEkzZGzeD+bw4FyEZ9PIGkGJQEsESEZH9k9wSGUqSujgpUyZsqEqmustsL/OEaOkZkuCnfvHIxGOyz9+zLb+/OQRGREQuQ4KcOkWzq9lekvAsOT8y7JVUPT9kvxgAERGRS5Fgx9ZT3cnxMAeIiIiIXA4DICIiInI5DICIiIjI5TAAIiIiIpfDAIiIiIhcDgMgIiIicjkMgIiIiMjlMAAiIiJyAiNGjMAHH3yg2+PFjBkz1JplUsV506ZN6NevH4YMGQJ7xACIiIgoCcnq6nnz5lXrdyWlBw8eqFXjrfH555/jo48+SvDjoyOr1sv6X6NGjcL27dvVmmb37t3D/fv3Y31dvbASNBERURJ69uwZrly5gpcvXybpeZbAw9rXePDggUVgEt/HR0cWQPX19UWDBg1M+6ZMmYIUKVLE+rp6YQ8QERGRzhYtWoTq1aurVdMbN26M/fv3v3LM5MmTUaZMGVSoUAGDBw/GyJEj0aVLF9P9P/zwA4YPH266LSu1d+zYESVKlECNGjUwc+ZMyPrn3377LX799VesX79e9UzJJr1UUR8vpCenSZMmKFKkCJo2bYqAgIBo2y/Bk/T+SIBjfE7ZZDhs4sSJ6piYXlcv7AEiIiLHZDAAISH6vHbq1ICbm82Cnx49emDSpEmoWLGiui0By/Hjx1XgIWbNmqWGj6ZNm4ZSpUqpYOabb75BnTp1YhzCkkCqSpUqKuh48uSJuly3bh26d++uApnAwEBMnz5dHZs9e3YsXrzY4vGrVq1C27Zt1et+/fXXuH79usrnkdyeqPr27atWZx8zZgx27Nhh2i85QMYen5heVy8MgIiIyDFJ8JM2rT6v/eQJkCaNTZ7qq6++woABA0wJyBLgHDhwAGPHjsX8+fPVPn9/fwwaNAidO3dWt6dOnWoRaEQVGhqKs2fPYuHChShZsqTaJ8FQeHg43N3dkS5dOoSFhalemJhIL1OvXr3w5ZdfmtplPrxlToa+MmfOrJ7b/DnTmJ2jDBkyWPW6yYVDYERERDoJCQnB+fPnVY+PuVq1aqkeIPH06VNcunQJ1apVM93v5uaGqlWrxvi8KVOmVMFK+/bt1RDUkSNH1PCXBCjWePDggWpX/fr1Lfab5/M4OvYAERGRY5JhKOmJ0eu1bUB6aoS3t7fFfrltvC+2Y2KzZs0arFy5Ehs3blT5Q6lSpVK3ixcvHme7Xrx4oS7lMc7KeUI5IiJyLZKDI0Msemw2yv+RoaNMmTLhxIkTFvul96dQoULqutyfMWNGnDx50uKYqLej8vDwQOvWrfHzzz/j8uXL8PPzw7hx49R97u7uqkcoJlmzZlVtO3z4MGwprtdNTgyAiIiIdCQJxJLQLENOQnpsVqxYoRKIjXr27IkJEybg6tWr6rbcH1sO0M2bN1XCsjGpWYbagoKCVFAjcubMiYsXL6qcoOjIUFf//v3Vaxpnasnjv/jii0T9rHG9bnLiEBgREVEyqFy58is5OH/99ReGDh2qAhaZrp46dWpERESopOd69eqZjpNE5HPnzqkKy5JMXKBAATUtXY6NTpYsWVSvkSRAS49LcHCwmjEmU+dF586dVYK1JC77+PioGWJRyZR4GQqTXCJpt2zGxydUdK8bWy5TUnIz2EtfVDKSKFZO/OPHj9W0PSIisn+SCyPJwBIESJKvI7U7pirQMixlTCyWgokyRVyGn2T4KjqPHj1S96VNm1ZNc8+TJw9+/PFHdd/Dhw9VIUN5vDmpxizfeV5eXq8838OHD9V3okxHl16i6B4v+2Qqe7Zs2VTydUwkWVueT9pk/try88kwXkyvG/V3Gdvv2Zbf3+wBIiIiSkLyJW7NtG9JOI4p6Vh6fyQfp02bNioAWr58uSoouGXLFtMxxuGt6HqDYuLr62t6XExBpbyeNfV6ZMq7+bT32F7b/HX1whwgIiIiOye9KlKAUAIKqaUjeUNSDLFmzZp6N81hsQeIiIjIzknP0C+//KJyciSfR+/eE2fAHiAiIiIH4enpyeDHRhgAERERkcthAERERA7FBScvuxRDMv1+GQAREZHDDP8Ima5Nzivkv9+v8fedVJgETUREDkEK8UkRQGN1YykaGFtdGnK8np+QkBD1+5Xfs7ULtyYUAyAiInIYxno0xiCInE+GDBmsqjuUWAyAiIjIYUiPT44cOVS1YuOK5eQ8PD09k7znx4gBEBERORzj2lRECcUkaCIiInI5DICIiIjI5TAAIiIiIpfDAIiIiIhcjmsHQL/9pncLiIiISAeuHQB17w588IGUndS7JURERJSMXDsAErNmARUqAKdP690SIiIiSiauHQCtXg1kywacOgWUKwfMmSO1uPVuFRERESUx1w6AatQAjh8H6tQBnj3ThsQ6dgSCg/VuGRERESUh1w6AhPQAbdgA+PtLaVFg0SKgbFng5Em9W0ZERERJhAGQOgspgCFDgJ07gTx5gPPngYoVgaVLk+q8ExERkY4YAJmrWhUICABq19ZmhrVtC3z6KfDypW6/ICIiIrI9N4NB36zfkJAQrF27Fnfu3EGJEiXw9ttvx/mY27dvY9u2bQgMDISfnx/q1asHLy8vq18zKCgIPj4+ePz4MdKnT//qARLwDB8OjB+v3a5ZU+sNypIlXj8bERER2U6c39+O0gN069YtlCpVCmPGjMHRo0fx3nvvoaMkIcdixYoVyJcvH3799VecPn0aQ4cOxeuvv47r16/brmEeHsC4ccDy5UCaNMD27Vpe0OHDtnsNIiIics0eoM6dO+PkyZPYv38/vL291XUJiCTIad68ebSPkV6iihUrYvbs2ep2aGgoChQogG7dumH06NG2jyClPlCLFsC5c4C3NzB9OtCtW/x/WCIiIkoUp+gBCg8Px8qVK9GlSxcV/IjixYvjrbfewrJly2J8nBybRnpl/uPp6amGv1KlSpU0DS1aFDh0CGjaFAgL06bK9+3LvCAiIiIHplsAdO3aNTx9+lQNX5krUqQIzpw5E+PjZsyYgV27dqFnz574+uuvUb9+fZQrVw79+vWL8TFhYWEqajTf4sXHB1i1CpAeJjc3YNo0oFEj4PHj+D0PERERuXYAFPxfscEMGTJY7Jfbxvti6jmS7erVq7h586ZKiBYRERExPsbf3191mRm3PDLVPSFT5SUxesUKIHVqYNMmoHJl4OLF+D8XERERuWYAlFqCiP/G88zJuJ75EJe5ly9fokWLFqrXZ+PGjZg+fToCAgJUj9HgwYNjfC1JlJbnNW7S+5Rgkg+0Zw+QKxcgPVWyjtju3Ql/PiIiInKdAOi1115T+TwXLlyw2C+3CxUqFOOsMen1eeeddyxygKpVq4bDsczQkteRZCnzLVFKl9bygmRm2IMHWt2gX35J3HMSERGR8wdAHh4eaNSokZrObhy+unz5Mnbu3Kl6eYw2bdqEefPmqes5c+ZUyc4HDx403S+T2I4cOYKCBQsm7w+QMyewaxfQqhXw4gXw/vvS1SRjccnbDiIiInKsafDS21OlShU1tb1ChQpYsmSJ6v1Zt24d3GVdLgA9evTAgQMH1BR5MW3aNHzyySdo164d8ufPjy1btqj79uzZg6IyYyuZp9GpgOeLL4CxY7XbLVsCCxcCSTUrjYiIyEUF2fD7W/dK0Pfu3cPixYtNlaClGKIx+BHLly9XCc8DBw407Tt16hQ2b96MBw8eqKG0d99995Vk6mQLgIwWLJBoDXj+HKhSBfjjD8DX1zbPTURERHCqAEgPSRIACRkSk3pBMj2+SBFtlfnXXrPd8xMREbmwIGcohOiUqlfXZojlzg3884/WE3TihN6tIiIioigYANla8eLAvn1AsWLAzZtAtWrAjh02fxkiIiJKOAZASUEKLUptIOkRkjpH9eppq8kTERGRXWAAlFQyZgQ2btSmyUtidNu2wKRJSfZyREREZD0GQEkpZUqt50cWTxUff6zVCnK9vHMiIiK7wgAoqcmU/smTgXHjtNtyKQERCyYSERHphgFQcpAV5GWtshkztOvTpwNdu8riZsny8kRERGSJAVBy+uADrUq09ArNnw+0aQOEhSVrE4iIiIgBUPJr3x5YsQLw8gJWrgSaNQNCQvheJCIiSkbsAdKDBD1r1wKpU2szxWSavFSPJiIiomTBAEgv77wDbN4M+Pho1aNr1wbu39etOURERK6EAZCeZKmM7duBzJmBo0eBGjWAO3d0bRIREZErYACkt9KltarRuXLJMvdArVoMgoiIiJIYAyB7ICvHy3phsojq6dNAzZrA7dt6t4qIiMhpMQCyFwULRgZBZ84wCCIiIkpCDIDsSYECWhAki6n+848WBN26pXeriIiInA4DIHsMgiQxmkEQERFRkmEAZO89QWfPsieIiIjIxhgA2av8+bUgyM9PC4JkivzNm3q3ioiIyCkwAHKUIOjcOa144r17ereKiIjI4TEAsnf58lnODqtbFwgM1LtVREREDo0BkKMEQVu2AFmzAseOAQ0aAMHBereKiIjIYTEAchSvv64FQb6+wMGDQJMmXEWeiIgogRgAOZISJbTV49OnB3buBFq1AsLC9G4VERGRw2EA5GjKlQPWrQNSpwY2bADatQNevtS7VURERA6FAZAjqloVWLMG8PYGVq0CunQBwsP1bhUREZHDYADkqGrXBlasADw8gEWLgF69AINB71YRERE5BAZAjqxRIy34SZECmD0bGDpU7xYRERE5BAZAju6994BZs7Tr48cD336rd4uIiIjsHgMgZ9Ctmxb8iE8/BebP17tFREREdo0BkLMYNAgYODAyIPrjD71bREREZLcYADkLNzdgwoTIGWGtWwN79ujdKiIiIrvEAMiZSDK05AM1bgyEhmqXJ07o3SoiIiK7wwDI2Xh6AkuXAm+9BTx+DNSrB1y6pHeriIiI7AoDIGckVaIlB0iWzrh9G6hTB7hzR+9WERER2Q0GQM4qQwZt3bC8eYELF7SaQU+e6N0qIiIiu8AAyJnlyAFs2gRkzgwcPQq0acN1w4iIiBgAuYBChbThsFSptEVUP/qIS2YQEZHLYw+QK6hUCVi8OHKW2NixereIiIhIVwyAXEWzZsDkydr14cNZLZqIiFwaAyBX0qcP8Nln2vXu3YHNm/VuERERkS4YALkaf3+gXTstGbpVK+D4cb1bRERE5BgB0M2bN23fEkoekgc0dy5QowYQHAw0bAhcu8azT0RELiVBAZCfnx8aN26MlStX4sWLF7ZvFSUtb29g1SqgWDGJZoEGDbSq0URERC4iQQHQ5s2bkTFjRnTs2BG5cuXCwIEDcerUKdu3jpK2UOL69VqtIPndsUYQERG5kAQFQDVr1sSCBQtw+/ZtjB49Gnv27EHx4sVRsWJFzJgxA0FBQbZvKdlenjxajSBZOkOqRv/vf6wRRERELiFRSdDp06fHhx9+iAMHDmDy5Mk4duwYevXqhRw5cqBfv3549OiR7VpKSaNsWeDXXwE3N+DHHyOnyhMRETmxRAVAV69exahRo1CgQAEMGTIE7dq1w+7du1Xv0MGDB9GiRQvbtZSSTvPmwMSJ2vWPP9Z6hYiIiJyYm8FgMMT3QUuWLMGcOXOwdetWvPnmm+jRowc6dOigeoSMHj9+jMyZM9tlkrQM0fn4+Kg2mrfZpcnboFcvYOZMIE0aYM8e4M039W4VERFRknx/J6gHSIa5ChYsiCNHjuDo0aPo3bv3Kw2RBo4YMSJRjaNkJENgU6cC77wDPH0KNG6szRAjIiJyQgnqAXr27BlSyeKaDoo9QLGQvK0qVYAzZ4AyZYBdu7QeISIiIlfvAfL09FQ5PlHJvpdSYZgce3r82rVAlixAQADQoQMQHq53q4iIiGwqQQGQTH3fKNOmo9iwYQPGjBlji3aRnvLlA37/XSuYuHo1MGwYfx9ERORUEjQEJsUPJf9Hprubu3XrFipVqoQrV67AnnEIzEqLFwPt22vX588HOnVKwt8KERGRnQ+ByQtHFzdFRETg/v37iWoQ2RFZNNXY+9Ozp4xx6t0iIiIim0hQAFS5cmVMNNaNMTN+/HjVA0ROZPRooFkzICxMqxd0/breLSIiIko0j4Q8yN/fH2+//TZ27dqFatWqqd4gKYB49uxZ7NixI/GtIvtaPX7BAm1m2MmTWhAkM8Nk+QwiIiJX6gEqV66cqv8jl3v37sX+/ftRvnx5BAQEqEtyMunSAWvWAJkzA0ePAt27c80wIiJyvSRoR8ck6ATauVMrlCilDr7+Gvj8c9v+YoiIiOw5CdooNDRUrQgfdSMn9fbbwLRp2vXhw7Up8kRERA4oQQHQyZMn1VBX6tSp1VT4qBs5sQ8+APr21a5LkcS//9a7RURERMmTBN2zZ0/4+fnhu+++Q8aMGRPyFOTIvv9eWypj61agaVPg0CGtcjQREZEz5wClSZMG169ft0nwc+/ePbW6/J07d1CiRAm8++67cHd3j/Nx+/btw7Zt21QvVJs2bVRxRmsxB8gGHj4EKlQALlwAqlcHNm8GvLxs8cxERET2mQMkvT9PnjxBYl24cEEFPb///jvCw8MxdOhQNGzYUF2PiRRb7Nq1K5o2bYrAwEC1yWP++eefRLeH4sHXV5sZJjPEZFp8v36cGUZERM7dAzRz5kysW7dOXWbNmjXBLy69PZI0LfWEUqRIgcuXL6Nw4cKYO3cuOkh+STSmTp2KwYMH49ixYyhUqJApIpQV6rNly2bV67IHyIbWrQMaN9aCH0mQ/ugjWz47ERFRknx/JygAypAhg3pxIQ1wc3OzuP/Ro0dxPoesGp8uXTqVR9S7d2/T/jp16qjnX758ebSPK1KkiCq+OGvWLCQUAyAbk6rgn30GeHgA27YB1arZ+hWIiIhgy+/vBCVBz549O9G/BlkwVabRFyhQwGK/3JbCitEJDg5W1aaHDRumhs2kGGPOnDnRsmXLWHt/wsLC1GZ+AsmGPv0UCAgAliyRbj2tWGLu3DzFRERktzwSOnSVWCEhIeoyagQnkd3Tp0+jfYyx1+nbb79F9uzZUbVqVfz5559qSGzLli2oIEm5MSzdMXLkyES3mWIgPYASFJ8+DZw4AbRqpRVNTJmSp4yIiOxSisQmMW+W2T8JkDZt2miHy+S2DI1Fx7jf19cXGzduxJdffom1a9eqQEgSqGMi90nwZNyuXbuWoDZTLNKkAX7/XUuOlmnxffowKZqIiJwrAJKp6zVr1kTBggVRt25d0/4GDRpg+/btVs8kk+n0MqRlTmZzvfHGG9E+RnqHZMgr6npjcvvixYsxvpa3t7fqaTLfKAnky6cNg8kCqnPmAD/+yNNMRETOEwANHDgQmTJlUrV7zH322Wf4WtaIsoLU+pHcnV9++cWUnyMVpvfs2YP33nvPdJwkQ8uQl1Hbtm2xc+dONR1eSA63rEBfsmTJhPwoZGt16gDjxmnX+/cHdu/mOSYiIruToFlgknAs09Bl2QuZAWZ8Chlekn3G/J643Lp1S83okuGwMmXKqHwe6VFauHCh6ZgePXrgwIEDKjgyvkbt2rVVraDKlSvj0KFDqkdKiiJGTaiOCWeBJTF5P7Rvr/UGSZkEJkUTEZEzTIOX6ssSvEgjpH6PsTdG6viUKlXKlKxsDQmWJPAxVoKuUaOGxf2bNm3CzZs38f7775v2vXjxAhs2bFAzyWQoTYKmlPFIuGUAlAwkkb1KFS0pWpLTmRRNRESOHgDVqlULrVu3Rq9evdRQlvTGyNa9e3fcv39fBTT2jAFQMpG8rHLlgMBAoFs3baZYlJpRREREDlMHaMKECWoYSqaeS/zUr18/df3GjRsqh4dIyZ8fWLoUqF9fS4ouW5aVoomIyHGToMuVK4cjR46oXCCZgSU5OtWrV1eFCZmMTBaYFE1ERHYoQUNgjo5DYMmMSdFEROQMq8HLOl6xbUTRVoqWUgV37wJS5uD5c54kIiLSTYICIE9Pz1g3omgrRa9cKSvpAgcOAJ98wpNERES6SVASdNRqzzIN/vz58xg9ejQ+4RcbxUTqNEmNp8aNgWnTgIoVgU6deL6IiMixc4BkBtjnn3+uKjXbM+YA6WzECGDUKCBVKmD/fqBUKb1bREREDkD3HKCYlC5dGn/99Zctn5Kc0ZdfalPjnz0DWrbU6gQRERElI5sGQLKuV5YsWWz5lOSM3N2BX38F8ubViiV27izjqHq3ioiIXEiCcoBkFfioAgMDERwcjLlz59qiXeTsfH2BFSu05TKkcvjYscDw4Xq3ioiIXESCcoB++umnV/ZlzJgRlSpVwmuvvQZ7xxwgOyIBsyyTIVPl168H6tXTu0VERGSndF8LzNExALIzH34IzJyp9QrJyvEyNEZERORohRBZFJHiZfJkoHx54OFDoFUrIDSUJ5CIiByvECKLIlK8eHsDv/0GZMoEBAQAffvyBBIRkf0lQU+cOBFjxoxB//791WKo4vDhw/jhhx8wfPhwVJQCd0Tx4ecHLFmi5QD9/DNQqRLQowfPIRERJYkE5QBVrlwZo0aNQh1Z6dvMpk2bMGLECOyX4nZ2jDlAdszfHxg2DPDyksqa2tAYERER7CAJOm3atLh58+YrLy4Nyp07t5oOb88YANkxqQckxRFXr9Z6hSQpOnNmvVtFRER2QPck6KxZs2L+/Pmv7Jd92bJlS1SDyMWlSCEVNaXYFHD1KtCuHRAerneriIjIySQoB2j8+PFo164dVq5cqXKApBPpyJEj2LVrF5YuXWr7VpJr8fHRVo6XPKAtW7SlM8aM0btVRETkRBLUA/Tee+/h+PHjKFSoEPbt26dyfuT6iRMn0EqmMRMlVokSwOzZ2nWpEi1DYkRERDbCQoiJHEOkJNa/v1YnSH5PR44AhQrxlBMRuaggvXOAjC5cuIDNmzcnqgFEsZo4EahaVd71WnL006c8YURElGgJCoDu3buHmjVrqkVR69ata9rfoEEDbN++PfGtIjKS6fDLlwPZswMnTwIffAC43uotRERkDwHQwIEDkSlTJty5c8di/2effYavv/7aVm0j0uTIASxbBri7A4sWAdOm8cwQEVHy5wDJVPdjx44hR44ccHNzU7PAhIzJyb6QkBDYM9YBclDffSfRt6zFAuzcKRU59W4RERG5Ug6QFDpMnTq1ui4BkFFgYKBaB4woSXz8sUxBBF68AN59F4jSA0lERJSkAVClSpWwePFiiwAoPDwcX331FapVq5aQpySKm7zXZJ2wIkWAmzeBtm2Bly955oiIKHkKIU6YMAG1a9fGli1b1PBXv3791PUbN25gj6zfRJRU0qXTiiRWqADs2AF8/rlU5uT5JiKipO8BKleunKr8LLlAUgn6wIEDqF69Oo4ePYqSJUsm5CmJrPfGG8CcOdr1CRO0gIiIiCipk6DHjRuHIUOGwFExCdpJSEK0JEZLr9Dhw8Drr+vdIiIicubV4L29vfH06VN4eCRoBE13DICchCRD164N7N4NFCsGHDwIpEmjd6uIiMhZZ4GVKlVKrf9FpCuZcSj1gaRO0KlTQM+eLJJIRERWSVAXTuvWrdG2bVsMGjQIRYsWhZdU6zVTo0aNhDwtUfxJhWgJgmrWBGRmotQG6tePZ5KIiGw/BGZe+yc6CXjKZMUhMCc0aZJWJ0iGZWV2mKwfRkRETiVI7yGwFy9exLoR6bJqfJs2Wl2g1q1ZJJGIiGwTAGXOnNl0fcCAASoBOqaNKNlJr+Ts2doUeSmSaAyGiIiIEhMAyayvZ8+eqevTuBgl2aO0abWaQHIpa4UNHap3i4iIyE5Z3V0jBQ9btGiBsmXLqtvDhw+P8ViuCE+6kWUy5s7V1gz75htZtwVo1Yq/ECIiSlgAtGDBAvj7++Og1FoBuOQF2S9ZKPXTT7UAqGtXoHhxFkkkIqLEzwJLmzYtnjx5AkfFWWAuQPJ/3nlHGworWlQrkihDY0RE5LB0nwXmyMEPuQhJxl+yRCuSePo00KMHiyQSEVHiAiAihymSuHy5FgwtXQpMnqx3i4iIyE4wACLnJgURv/1Wuy55QXv26N0iIiKyAwyAyPnJ0hjt2kUWSbx9W+8WERGRzhgAkWsUSZw5U1sx/tYtrUgiK5YTEbm0BAdA165dw7fffos+ffqY9m3evJlLYZB9F0lMlw7YtYtFEomIXFyCAqADBw6gWLFiWLFiBaZPn27av3btWsyU/7SJ7FHhwsC8edp1yQv67Te9W0RERI4UAA0aNEgVRdy3b5/F/h49emDq1Km2ahuR7bVsCXz2mXZdiiSeOcOzTETkghJUCDFdunS4deuWKoiYIkUKREREmNYLy5gxI54/fw57xkKILk6SoevWBbZv1xZPPXSIRRKJiByA7oUQvb298fDhw1f2nzhxAlmzZk1Ug4iSnNQFWrwYyJlT6wHq3p1FEomIXEyCAqDmzZtj6NChCAsLg5vMsAFw6tQpfPDBB2gpQwxE9i5bNi0HyNMTWLYMmDRJ7xYREZG9B0ATJ07ExYsXkSlTJjX8lSdPHhQvXlwNjXEleHIYlSsD332nXR80CNi9W+8WERGRPecACQl8tm7diiNHjqjrZcqUQb169VROkL1jDhCZyNu/Y0dg0SJt6YyAAG39MCIicurv7wQFQJUrV8b+/fvhqBgAkYWnT4FKlYCTJ4Fq1YCtW7WhMSIisiu6J0GfPn0awcHBiXphIruRJo1WJFE+TDIMNniw3i0iIqIklqAAqHHjxpg/f77tW0Okl0KFgF9+0a5//72WGE1ERE7LIyEPcnd3R9++fbF8+XIULVoUXl5eFvdP4owackTNmwNDhgDjxgHdugElSmh1goiIyOkkKACSGkCNGjVS169evWrrNhHpZ/RorTDitm1aQCTXfXz4GyEicjIJngXmyJgETbG6dw8oW1ZW/JXxXmD1asABZjcSETm7IL2ToImcWpYsWlK0tzfw55/AyJF6t4iIiOxhCKyj1E2JxcKFCxPaHiL7UK4cMHMm0KULMGoUULq0NiRGREROIUE9QB4eHhabFD+8cOECfv31V4SGhtq+lUR66NwZ+N//tOudOnHleCIiV+8BmjdvXrT7J0yYgEuXLsXruS5fvoy5c+fizp07KFGiBLp3746UKVNa9dgdO3aottStWxft27eP1+sSWeWbb4Djx4GdO5kUTUTkRGyaAySLof7xxx9WH3/y5EmUKlVKFVYsVKgQfvrpJ9SoUQMvXryI87F3797F+++/jz///BOHZKYOUVIwLpaaJw9w7py2bEZEBM81EZGDs2kAdOPGDYSEhFh9/ODBg1GhQgVVT2jgwIHYsmULjh07FmeRRZm41qlTJ/Tv3x+5c+e2QcuJYpE1K7BqFSA9k5IU/dVXPF1ERK44BBbdiu+BgYFYunQpmjVrZtVzPH/+HJs3b8b06dNN+7Jly4ZatWqpXiQZCouJDLXJAqwDBgzAL8bqvURJSabFz5ql5QJJrSBJim7RgueciMiVAiAZdooqY8aM+Oijj1RQYg0poChDXXnz5rXYL7d3y3pMMTh48CC+//57BAQEwM3NzarXCgsLU5t5HQGieJPhr6NHpdS5liB98CBQtChPJBGRqwRABw4cSPQLP3v2TF2mkYUozaRLl850X1RS+Khdu3aq1yhnzpxWv5a/vz9GspYL2cLEiVpS9PbtkUnRGTLw3BIRORjdCiEaKzg+evTolWU2pMpjdBYtWqTuX7NmjUqAlk16kjZt2qSuy7BYdIYOHaqCJ+N2TSr8EiWEhwewdCng5wecPw906ACEh/NcEhE5aw9Qjx49rH7S2bNnx3lMnjx5VKBz6tQpNGjQwGJmWPHixaN9TM2aNV9ZaHXbtm0qEVpmj8U0JObt7a02IptVipak6KpVgXXrgC+/BMaM4cklInLGAOjJkyc2fWEpntimTRvMmTMHvXr1Qtq0aVV+j2xfmc2ykRpBUmRREq+LFCmiNnMSEMmK9NIDRJRsypSRSF/LCxo7Vls5vm1b/gKIiJwtAFqyZInNX1xyc9555x1VAFE2KWzYr18/1KtXz3TM3r17Vc5RdDPPiHQlw1+SDyR5QV27AgULaktoEBGR3dN9NfiXL19i165dpkrQUYe/9u3bh3v37sU4vV7ygbJnz67qCVmLq8GTzUj+j7w3164FJDH/8GHtkoiIbM6W399WB0DGHCDJ74krH8iaHCA9MQAiG7+hgMqVgdOnAQnEd+wAUqXiSSYisuPv7wTlANk6H4jIocmHcM0aLfiRafE9ewILFgBW1qkiIiIXHALTA3uAKEls2wbUrasNi/n7A0OG8EQTEdnp97dudYCInE6tWsCUKdr1YcOAeCwMTEREDlAJ2riWlyxcKoUIJZHZXFtOByZX1bs38PffwI8/Au3bA/v3AzHUtSIiIgcbAjtz5gyaNm2qKirLGluynMXTp0/VfZkyZcL9+/dhzzgERknqxQtASjnIchmy1p3MDMucmSediMjRh8BkwdNGjRqZkqHlUio6ly9fHoMHD05Ug4gcnqcnsHw5kD8/cPky8O670mWqd6uIiCixPUC+vr44d+4cMmfOrCo6Sy+Qp6cnTp8+rQKjS5cuwZ6xB4iSxalT2vT44GDggw+An37izDAiIkfuAQoMDFTBj5DLmzdvqut+fn64detWohpE5DSKFQMWL9aCnpkzgalT9W4RERHZahZYpUqV1DIVsojpiBEjULhw4cQ+JZHzaNQIGDdOuz5ggLZ4KhEROWYA1L9/f9P18ePHqzW8ZBmLBQsWYIpxGjARaQYNArp1AyIigDZtgBMneGaIiBwpB+jTTz/FN998Y7p9/fp15M6dW12XmV+SGyQ5QfaOOUCU7CQJun59bWZYnjzAwYNAjhz8RRAR2ftaYOpgNzeYHx71tqNgAES6CAzUkqLPngXKlgV27gTSpOEvg4jIUZKgiSgBMmbUVo2XCQRHjwKdOmnDYkRElOwYABElpwIFgN9/B7y8gFWrkm29sPAIA/ZfeIDVx26oS7lN5Oj4vqZkXQrj0aNHsd4WGTJkSFSjiJxa1arAvHnaUhkTJwKFCmkryCeRDSdvYeQfp3HrcahpXw6flBjRpCjqF2ceEjkmvq8pseKdA2QNe88LYg4Q2YXRo4EvvwTc3YENG4B33kmSL4neCwMQ9RNp/CT/2LEMgyByOHxfu64gG+YAxasHaJV02RORbQwfDpw/DyxYALRqBezZA5QoYdPhAen5ie7fEcN/QZDcX6dodrinsO6fGyK98X1NthKvAKh58+Y2e2Eilyc9qrNmAVeuALt2AQ0baqvH/1daIrEOXXpoMewVXRAk98txlQtkcvlfBzkGvq/JVpgETaQnb28tGfqNN6SwlhYEPX5sk6e+Gxxq0+OI7AHf12QrDICI9ObrC6xfD2TPDvz9N9CypU1Wj8+aLqVNjyOyB3xfk60wACKyB6+9pgVBadMC27ZpS2ckcjJBhXy+arZXTNk9sl/ul+OIHAXf12QrDICI7MWbbwIrVgAeHsCvvwKff56op5PEZpnqLqIGQcbbcj8ToMmR8H1NtsIAiMie1K0LzJ6tXff3B378MVFPJ3V+ZKp7dh/LYS65zSnw5Kj4vqZkrwPkLFgHiOze118DX3wByOLCK1cCzZoleuqwzJ6RBFLJoZBhBPb8kKPj+9r1BOm1GKqzYABEdk8+lh9+qE2TT5UK2LpVW0iViMiFBXExVCIXqBE0fbo2Lf7ZM6BxY+D0ab1bRUTkNJgDRGSvJBl62TKgYkXg4UMtP0iKJhIRUaIxACKyZ2nSAGvXaoUSb9zQgqB79/RuFRGRw2MARGTvMmUCNm0C/PyAc+eABg2A4GC9W0VE5NAYABE5AlkfbPNmIHNm4OhRWZgPCOUSFkRECcUAiMhRFC4MbNgQWS26QwcgPFzvVhEROSQGQESOpGxZYPVqwMtLqw/Uu3eil8wgInJFDICIHE2tWsDixVqRRKkTlMglM4iIXBEDICJHJCvG//RT5JIZshERkdUYABE5qp49gfHjtevDhgE//KB3i4iIHAYDICJH9tlnwIgR2vUBA7QhMSIiihMDICJHJwHQp59q12X9sIUL9W4REZHdYwBE5Azrhk2YAHz0kTYjrEsX4Lff9G4VEZFdYwBE5CxB0JQpQNeuQEQE0K6dtoQGERFFiwEQkbMwTotv2xZ4+RJo1QrYskXvVhER2SUGQETOxN0dmD8faNYMCAsDmjYFdu3Su1VERHaHARCRs/H0BJYuBerVA549Axo2ZBBERBQFAyAiZ+TtDaxaBdSpAzx9qq0gv3On3q0iIrIbDICInFWqVNq6YXXrAiEhWk/Qjh16t4qIyC4wACJyhSCofv3IIEhWkicicnEMgIicXcqU2nCYDINJTlDjxsDWrXq3iohIVwyAiFwpCGrUKDII4hR5InJhDICIXCkxesUKLfgJDQWaNAE2bdK7VUREumAARORqQZAskyH1gSQIkss//tC7VUREyY4BEJErBkHLlwMtWmjFEuVy0SK9W0VElKwYABG5Ii8vYNkyoFMnIDwc6NgRmDFD71YRESUbBkBErsrDA5g3L3IV+V69tFXliYhcAAMgIldfQHXqVGDoUO324MHA8OFaQERE5MQYABG5Ojc3YOxYYNw47faYMUD//kBEhN4tIyJKMgyAiCiy92f6dC0gmjIF6NYNePmSZ4eInBIDICKK1Ls3sGAB4O4O/PIL0KqVtoQGEZGTYQBERJY6dNAKJkr16DVrgNq1gfv3eZaIyKkwACKiVzVrpi2VkTEjcOAAUKUKcOkSzxQROQ0GQEQUvapVgb17AT8/4Px5oHJlICCAZ4uInAIDICKK2RtvAPv3A6VKAXfuAG+/DWzcyDNGRA6PARARxS5nTmDnTi0X6MkTbTHV+fN51ojIoTEAIqK4+fgA69YB7dtrU+O7dNHqBbFgIhE5KA+9G3DixAnMmDEDd+7cQYkSJdC/f39kyJAhxuNDQkKwYMEC7Nu3Dx4eHnjrrbfQuXNnuMu0XSJK2vXDZIp8rlzAxIlaxehz54CZM7UFVomIHIiuPUCHDh1CxYoVERERgSZNmmDjxo0qoHn27Fm0x8txxYoVw7Fjx1C7dm2UL18eI0eOVI+V+4goGZbOkPXCZPkM+adDhsJkaOzePZ56InIobgaDfn3YEsSkTZsWq1evVrcfPXqE3LlzY9y4cejbt+8rx0tT7927h6xZs74SRB08eBAVKlSw6nWDgoLg4+ODx48fI3369Db8iYhcyKZNQOvWwOPHQN68wB9/AMWL690qInJiQTb8/tatByg0NBQ7d+5Ey5YtTftk6EuCog0bNkT7GDc3N4vgR2TPnt10UogoGdWtq80QK1AAuHxZqxUkeUJERA5AtwDo6tWrCA8PVz0+5vLkyYNL8Si49t1338HX11f1AsUkLCxMBUjmGxHZaJr8wYNA9epAcDDQpAkwaRKTo4nI7ukWAD1//lxdpk6d2mK/3DbeFxdJhp46dSp+/vlnpEuXLsbj/P39VZeZcZMgi4hsJFMmYPNmbfFUycX7+GOgVy/gxQueYiKyW7oFQMaZXg8fPrTY/+DBg1hngRktX74c3bt3x6xZs9C8efNYjx06dKgaLzRu165dS2TrieiVGWKzZwPffKOtJi8zwyQ5+vZtnigisku6BUC5cuVCpkyZcPz4cYv9MsOrlFSdjcWKFSvQsWNH/PTTT+jatWucr+Xt7a2Spcw3IrIxCXwGDtQWUJXP2O7dQJkywL59PNVEZHd0C4AkoblTp06YPXu26vURmzdvRkBAgKrrYzRlyhT06dPHdHvVqlVo3749fvzxR3STLncisi9SKfrwYaBoUeDWLaBGDWDaNOYFEZFd0XUa/JMnT9C0aVNVDLFw4cL466+/MGzYMHzxxRemY3r06IEDBw7g5MmTKnk5c+bMKo8natLz//73P9SVWSlW4DR4omQgy2bIPynLl2u35R+bn34CUqXi6SeiBLHl97euAZCRBEBSCVqKHOaUdYei3Cf1gapXr44XL16oYonRkWEza5ObGQARJRP58/Ldd8DgwUB4OPDmm8DKlUC+fPwVEFG8OV0AlNwYABEls+3bgTZttIrRGTMCixYB9evz10BErlcIkYhcSM2aCD98BMGlygCBgUCDBogYMsRiqnx4hAH7LzzA6mM31KXctjeO0EZbeP4yAj/vvogvV59Ul3KbyNmwB4gzwoiS3IaTtzDyj9N48CAYX2ybhU5/aRWjH5Usiwyrf8OGJ97q/luPQ02PyeGTEiOaFEX94jns6mew5zbagv+605i1+xLMY7sUbkDPavkwtGFRPZtGBA6BJRKHwIiSN3DovTAA5n0l9c/uxYT1k5E+7CmepUmPT+r0wfrXq1o8zu2/yx87ltE9wIjuZ7C3Ntoq+JmxK+ZK/B9WZxBE+uIQGBE5BBkikl6TqIHDhteromHXyQjI+TpSPQ3Cj7/7Y/Sm6fB+EWY6xvgYebyeQ00x/Qz21EZbkGEu6fmJjdzP4TByFswBIqIkc+jSQ4shI3PXfbKhdfvx+LHiu+q2DIv9vmAgCtyPrNQuIYU8Xp7HHn8Ge2mjLSzYf9li2Cs6cr8cR+QMGAARUZK5Gxxz4CBeuntgfI330an1KNxLnQFv3LuMP+YPQAfJETKboBrX8yQla19bzzbawpWHITY9jsjeMQAioiSTNV1Kq47bna8MGnadgt2vvYnUL8IwZtN0zFv+FbIF34/X8yQFa19bzzbawmu+qW16HJG9YwBEREmmQj5fNVPKLbY/Qm5aMvG9tBnRuc0ojKrVE2Hunqhx6Sg2zumLTpf3q+ex159B9sv9erbRFjpVzqt+F7GR++U4ImfAAIiIkox7Cjc1TVxE/W51+2+T6dXG2wa3FJhTvhkavf8DTmQviAyhTzB66Ri4t2sL/LdmoL39DELul+McmZdHCtPvIiZyvxxH5Az4TiaiJCXTw2WaeHYfyyEiuS37pbZM1Pv/zeyHj/pMxfmPBgLu7sCyZUDx4sA6rX6Qvf0MzjAFXsjvQqa6R43l5DanwJOzYSFEFkIkShYyTVxmSkmysOTLyJCRea9JjPcfOaItpHrmjHZg9+7AxInakhp29jM4C5nqLrO9JOFZcn5k2Is9P2QPWAjRjk4gESWDZ8+Azz8HJk3SZodlzw5MmQK0agW4OV8AQkTRYyFEInItqVJpq8rv3Am8/jpw+zbw3ntA8+bA9et6t46IHBBzgIjIcVSrBhw7BnzxBeDpCaxZAxQtCkybBkRwwU4ish4DICJyLClTAqNGAX/9BVSqBAQHA337asHR6dN6t46IHAQDICJyTMWKAXv2aLlAadMC+/YBb74JDBkCPHmid+uIyM4xACIixyVT5KX3R3p+mjQBXrwAxo8HihQBFi+2WE6DiMgcAyAicnx58gCrV2s5QfnzAzduAO3bAzVqACdO6N06IrJDDICIyDnIdHjpBTp1Chg9Wps5tmsXULo00K8fEBiodwuJyI4wACIi50uSHj4c+Ocf4N13tdlhU6cChQsDM2YAL1/q3UIisgMMgIjIOfn5AcuXA1u2aFPl798HevUCSpTQhsqYH0Tk0hgAEZFzq11bqx30ww9Apkxaz1CzZsDbbwMHD+rdOiLSCQMgInJ+UjTxf/8DLlwAhg7Vhsl279bqCLVuDfz7r94tJKJkxgCIiFyHjw8wdixw/jzQtauWOC3DZG+8oSVKyxIbROQSGAARkevJnRuYMwc4fhxo2FBLjJZE6Xz5gIEDgTt39G4hESUxBkBE5LokIXrtWmDrVqByZSA0VFt0VQKhQYOAe/f0biERJREGQEREtWoBe/cC69cDFSoAz54B33yjBUKytIbMICMip8IAiIhISD5Q/frAgQNar1C5csDTp9rSGhIIDR4M3LrFc0XkJBgAERFFDYQkL+jQIa1eUJky2uKqEyYAefMCPXsC587xnBE5OAZARESxLa1x5IgWCFWtCjx/DsyerS222qqVFiQRkUNiAEREZE0gtGePtjVtqlWRXrkSqFgRqFlTyx2SJTeIyGEwACIispb0Asmq87Lg6vvvAx4ewI4d2pCZLLchU+mDg3k+iRwAAyAioviSYGfuXODSJeCTT4D06YGzZ7ViirlyaVWnmSdEZNcYABERJaag4rffAteva70/khskPUBTpgCvvw40aACsW8fhMSI7xACIiCix0qUD+vQBTp8GNm0CGjfWcoc2bAAaNQIKFgTGjAFu3OC5JrITDICIiGxFgp46dYA//tDWG5PhMVl/TIbKhg8H/Py0hGqZVSbLbxCRbhgAERElhQIFtOGxmzeB+fOB6tW1obA//wSaNdOCoWHDuBI9kU7cDAaZz+lagoKC4OPjg8ePHyO9JC8SESUHSZSWRVjnzQPu3o3cL+uQdewItG4NZM7M3wVRMnx/MwBiAEREyU0KKkpP0KxZWs6QsYaQTKuX5TgkGJKhstSp+bshMsMAKJHYA0REdkPWF1uyBFi4EAgIsEysbtkSaNcOarFWT089W0lkFxgA2dEJJCKymTNngF9/1YKhK1ci92fIoOUNvfuulmTt7c2TTi4piENg9nMCiYhsTobE9u0DFi3Slty4cyfyPvmbJcNjEgzVqwekSsVfALmMIAZA9nMCiYiSVHg4sHcv8NtvwIoV2qwyI8kRkh4hCYik3lD27PxlkFMLYgBkPyeQiChZe4YOHNACIQmIrl61vL98eS0YkkKMb76p1SUiciJBDIDs5wQSEelCKpgcO6YVXZQZZYcPv7pMh8woq1sXqF0b8PXlL4ocXhADIPs5gUREdjObTNYdk4Bo82YgJCTyPukJKldOC4Zkq1QJ8PLSs7VECcIAKJEYABGRUwsNBXbs0GoMyXbqlOX9adIANWpEbjJcJjWIiOwcAyA7OoFERHZPFmHdskULhqR36N49y/ul5lC1asDbb2sBUZkyDIjILjEAsqMTSETkcInUx48D27drvUS7dgGPH1sekzatNkwmS3TIJtczZtSrxUQmDIASiQEQEZHZNPsTJ7RgyBgQPXr06ukpUiQyIJKtaFEgBdfTpuTFAMiOTiARkdMFRH//DezfH7n9+++rx8nfzooVteTqsmW1YbO8eTn1npIUAyA7OoFERE5Pcoak/pAEQ3J56BDw9Omrx8mSHRIIyWYMigoWZE8R2QwDIDs6gURELuflS+DkSS0YkgVcZZNeI1nlPirJJypdGihZEiheXNuKFWNOESUIA6BEYgBERGRjEvycPq0FQ0ePapeSbP3sWfTH58oVGRAZtzfe0KboE8WAAVAiMQAiIkqmnqKzZ4G//tJ6jKSXSC6jLuFhzs8PeP11bStcOPJ6njwcSiMwAEokBkBERDqSaffSWyTBkHGT4ChqfSJzKVMChQpFBkYFCgD58mmbLPvBQo4uIYhLYdjPCSQiIhuRAEh6jM6d0y6N12UW2osXMT9Ogh/pOcqfPzIoMr+eOTNnpzmJIBt+f7P2ORER2YcsWbTtrbdeHUq7fDkyMJLLS5eAixeBK1e0/CO5Llt0JBH7tde0niIZSpPLqNf5z7DLcTMYZElh18IeICIiJ6psffOmFvwYgyK5NF6X+6why4GYB0Sy5cgBZM8euWXLBqRKldQ/EcWCQ2CJxACIiMiFFoaV3qNr17Tt+nVtM78eXeXrmPj4aIGQeWBk3LJm1YbbjJv0Krm5JeVP53KCOARGRERkBUmelmU8ZIvJkyeRwZB5cHTnDnD7duQWFqYlcMsmw3BxkdykTJksgyLzzfw+X1+tkKQEWEzoThbMASIiItcmOUJxBUmSLRIUZBkQRd0kifvBA+D+fS2oktwlCaJkiw8ZjpNgKOomC9JGt0+ON9+klhJ7nuw/ANq1axemTZuGO3fuoESJEhg2bBhyyLirjR9DROQIwiMMOHTpIe4GhyJrupSokM8X7ikih1GePQ/H2HWncflBCPJmSo1hDYsilZd7vJ4jrvufv4zAgv2XceVhCF7zTY1OlfPCyyN+C58mtg22eI3E/hzRPr/00MhU/LjaIENvxmAopu3BAxju38fzW3fgHvQYHiH/LS8SHKxt0hOVELJIrQR1UQOj2La0abXAKXXqmDdPTzgTXZOgt2/fjrp162Lw4MGoXLkypkyZgnPnzuH48eNIJ78QGz0mKuYAEZE92nDyFkb+cRq3Hoea9uXwSYkRTYqifvEc6Dn/MDafvvvK4+oUzYpZnctb9Rxx3e+/7jRm7b6ECLNvBvk+71ktH4Y2LGqTnyOu+23xGon9OaxpY2J/jqiP9wh/iYLe4RhWJTuqZ/XScpNkCwyMvG6+GfdLz5QxaErKr3QPj9gDJPNNksVl+FE2b2/rLmO6TwKv/3q0nCYJumrVqsiTJw+WLFmiboeEhKienC+//BIDBw602WOiYgBERPZGvgx7LwxA1D/Ixv6MErnT48T1oBgfL0FQqzK5Y32OD6rnw8xdl2K8/52iWaMNsIw+rB538BDXzxFXG37sWCbO4CGu10jszxHX80sbRVzHxPZzWPMa1gaDJvJ1HhISGQxF3cwDpeAomwzZyWOj22SmnZ4k+PkvIAry8oLP3buOHQBJ4CI9NvPnz0eHDh1M+1u1aqXuW79+vU0eEx0GQERkT2QY5a3x2yx6EhIia1pP3H0SfcFAt/++R8x7ROJLelD+Gd0gxmEka36OFLG0QdqY3Scl9gyuFeNwmC3OVWw/R1zPL63Klt5bXbsdFJqgn8Oa14jrPCQbg0GrsxRTcBTbJknjsslwoGzG69ZcRrewrnx/y0Q8VUzcgQshXrt2DREREciZM6fFfrm9detWmz1GhIWFqc1ITpwxECIi0tuhiw9x4+7DRD/P7cg/c0lC+gFmbvkbnavkS/DPEVdfwo27Idh+4goq5PdN8Gsk5uew5vlv3guJ8zVi+zmseY24zkOyc3ePzBdKDtLrJEGQBERml0GSV1W3LmzRd6NbAPTiv7Lm3tKtZSZVqlSm+2zxGOHv74+RI0e+sl+G0oiIyHr9JgH9kviE1ZmUxC/gID9HcpwHR/XgwQOVC+SQAZCv1DwA8PDhw1d+qExSG8FGjxFDhw7FJ598Yrr96NEjvPbaa7h69WqiT6Arkx40CSKlZ45rqvFc2gu+L3ke7Q3fk7YjIzh+fn6meMAhAyAZtsqePTsOHz6Mxo0bm/YfPHgQ1apVs9ljjD1GUXuNhAQ//OJOPDmHPI+2wXNpOzyXPI/2hu9J20khU/0T+xzQUbdu3TB79mzcuHFD3V6xYgVOnz6t9psPX5knPFvzGCIiIiK7LYQoU9fPnz+PggULqqGU69evY+rUqShfXqtnIS5cuKBq/MTnMURERER2GwDJsNSyZctw8+ZNVdVZgpqoxQylyvMTqU8Qj8dY87ojRoyIdliMeB71wPckz6W94XuS59LZ35e6FkIkIiIi0oOuOUBEREREemAARERERC6HARARERG5HF2ToJPS3bt3MW3aNFUjSJKlpE5Qr169kDZtWtMxkv4kU+pXrlyJly9fol69eujfvz88ZeVZesW9e/fQunVruLm5Ydu2bRb3SXFKKVlw5MgRVaCqZ8+eqF+/Ps+imenTp6t17MxlzZoVa9assdh36NAh/PDDD2qG4xtvvKEKeUrhTrJ06dIlTJo0CX///Tfy58+vzlOBAgX4+Y4HOX/GhaXNubu7Y+/evfx8x4N8nyxcuBCrV6/G/fv3kTt3brz//vt45513+PmOJ/k++fbbb3H06FFV76dly5bo3r27+u6x6fe3wQm9fPnS8MYbbxj8/f0NmzdvNixevFjdrlatmsVxw4YNM2TMmNEwd+5cw9KlSw1+fn6GDh066NZuexYREWGoX7++oXjx4oY0adJY3BcWFmYoUaKEOr9r1qwxjB071uDu7m5YuXKlbu21R4MHDzaUL1/esH//ftMWEBBgcczhw4cN3t7eho8//tiwbt06Q4sWLQw5cuQw3L17V7d22yM5T+nTpzd07drV9BmvWLGixTH8fMft8uXLFu9H2QoUKGCoV6+e6Rh+vq0zcuRIQ7p06Qw//fSTYfv27YavvvrKkCJFCsPvv/9u8b7l5zt2jx8/NuTPn99Qq1Yt9TdQvpsLFixoGDhwoM0/304ZAImQkBCL23IiJd67deuWuv3gwQODp6enOnlG8odUjjl58mSyt9feTZgwQf1RnDx58isB0Jw5cwxeXl7qnBp9+OGHhiJFiujQUvsOgGrXrh3rMQ0bNrT48nn+/LkKgIYPH54MLXQcEoi/++67FvuePXtmus7Pd8LI3z75G/jbb7+Z9vHzbZ3SpUsb+vbta7GvatWqKkg34uc7bj///LP6PpFAyEgCSvmn+vr16zb9fDttDpAskGpu9+7dqkvSuGbYrl271AKqTZo0MR1Ts2ZNVVNoy5Ytyd5eeyZLj3z//feYO3euRRek0datW1GlShWLtVmaNWuGf/75x1SxmzSnTp1C7dq11fkZN24cnj17Zjo1ERER2L59u8V7UrpzGzRowPekmb/++gsnT57ERx99ZPG2Spkypek6P98J8/PPP6th2aZNm/LzHU/lypVTQzZhYWHqttSp+/fff1GhQgV1m59v68janpKqYr68Uq5cuRAeHq6+a2z5+XbaAMg4vl2xYkUV+KxduxY7duwwjQ9euXIFXl5eFouoyrh3tmzZ1H0UuYhfu3btVP5Kjhw5oj0tcr5knTZzxts8l5EkF02WdRk0aBDee+89lQ9UqVIl0x9MybGSgCi6c8nzGEmWvhHyR1LOpwSUEgxJ1Xjz9yQ/3/Hz/PlzLFiwQOWtmOdR8PNtnSlTpqBYsWLq81q6dGkULlwYn376qco95efbetWrV0dgYKApN01GquT7R1y+fNmmn2+nTYIWLVq0UNH3xYsXMXbsWPTr1w9//vmnSqqS6DG6SpLScyT3kebDDz9UXzDNmzeP8ZREdy6NPXA8l5ZVzc3PU61atVQlc+lZkz+SxnMV3bnkeYwUGhqqeiI7deqEIUOGqP8O5RzKl86xY8dUQjQ/3/EnyfiSvNujRw9+vhNg3rx5KiF31KhRavKCjDp8/fXX6jtIvtT5+baOdFp899136n341Vdf4enTp2oSkwQ3kuwsbPX5duoASGbOyCbDM7JWWJEiRVQXWp06ddRwTXBwsDpZ5v/tSPebeVTpyqRnQqLwkiVLqp4KY7eu9FLI7cGDB6sgU86lZO2bk/MoeC4jRf3Ayn+KRYsWNa11lzFjRvXFHt255HmMJO83+a9QZn116dLFFEzmzZtXfQnJFxA/3wkb/qpRowYKFSpksZ+f77jJ8MzAgQPVF3afPn1M70lZt1KC9H379vHzHQ8DBgxQ/3xL54X09MqwrAxvyT87wlafb6ceAjNnHL6R/3BEmTJl1KVM2zaSace3bt1S/0kSVBfj/v37MWPGDDWcKFurVq3UF7lcr1q1qulcyti3+aoqUn4gTZo0r/wxpUhyviSglPMk5FK6zSXnypycS74nI5UtW1Zdmg8VSvd3lixZVNc5P9/xd+3aNWzatEmVr4iKn2/r/lkMCQkxfUEbyXvU+A8NP9/xI705MqQonRjr169XOVQyGmHT72+DE9q3b59h/fr1FjNpZApdqlSpDDdu3DDtL1u2rKFRo0Zq2rzo2bOnIVeuXBazScjSlClTXpkFdvbsWYOHh4dh1qxZ6vb9+/fVVFqZCUaRxowZY5qdGB4ebpome/DgQYvZdpkyZTL8+++/6vamTZvUMTLDgSI1aNBAzQKTz7bYtWuXmiViPuWYn+/4TeH29fU1hIaGvnIfP9/WKVeunJq6/eTJE3X75s2bhtdee83Qo0cPfr7jSUoJGD/bUqpBvk+6detmcYwtPt9OGQDJG0/+OMoXSalSpVStAKlTs2XLFovjzp8/byhWrJg6LmfOnIY8efJYfBmRdQGQWLhwoaqBIW/U1KlTq6ncQUFBPIVmpC5VlixZ1HtOLuX9Zj7dWMiH+f3331e1QgoXLmxImTKlYdy4cTyPUdy+fdtQpUoVQ9asWdX5lPfk6NGj+flOYI2vvHnzGvr37x/jMfx8x+306dPqSzlDhgyGkiVLqn+4Zdr7w4cP+fmOp++//159tuU8yveJ/DMt9ahs/f3t1KvBywwmqRabPXt2lUAVHfnxz549q5KrJHFNutIpZjJkc/XqVZVTFZV0AZ87d06Nz/r5+fE0RkPeZ3KOZDxbZidGV1bAeJ5v3rypEnp9fHx4LmMgn29JkpTzlDp1an6+Ezj7KyAgQA2/mpeyiIqfb+vI51ZmdJqXXYmKn++4yedayghIbl9MfwMT+/3t1AEQERERkUsnQRMREREZMQAiIiIil8MAiIiIiFwOAyAiIiJyOQyAiIiIyOUwACIiIiKXwwCIiIiIXI5TL4ZKRIlbdf3333+P9RhZYPjNN990ytMsi/6uXr0ajRo1UoUrici5MAAiohgXeDQPgM6cOYN//vkHLVq0MO1r0qSJ0wZAsrJ0u3bt1M8tgR4RORdWgiYiq4wbNw7ffPMN7t+/b7H/8ePHOHDggLouwZD5sjNPnjzBn3/+qQIlWSJAAihZ3blkyZLq/uPHj+Py5csoXrw4ChQoEO3jZJVnWT4kX7586riorH39ixcvqueRZVxkqZYlS5aoY6R8vrRJVpH29PRU+2Tl6Z9++gl9+vTBt99+q1b1zpo1K0qVKoXNmzfj3XffhYeH9v/jixcvsGLFCtSrVw8ZM2aM9TWNy3ecOHECmTNnVqtay6rXRJT82ANERAm2bNkyfPjhhyowSJkypQpEJFDq1auXuv/27duqF+Xtt99GYGAgcuXKhS1btmDAgAEq8Dl//rwKWHbs2IGff/4ZHTp0sHhcgwYN1HpAEvzs3r0bvXv3VgFJfF+/cePG6nlKlCihghkJRoy9WxLASCAmgciGDRtUG8PDw7Fx40Z1/9atW9UQWNGiRdV6Y8bnS5s2rWnNItl3+PBhlCtXLsbXlLWhJKBavnw5KlWqpNaDkmBy1apVTtuLRmTX4rtKKxG57mr2svKy+WrMsgr77t27TfsOHDigVrA/d+6c6Rj5M9OrVy/TMVOmTFH7BgwYYPHc+fLls3huOaZevXqGFy9eqH379u0zuLm5qcv4vn67du0M4eHhMf5scl/z5s0NPXr0MO27du2aeuyZM2dM+/bv36/2BQcHm/YFBgaqfYcPH471NSdPnmwoUqSIOt5o+PDhhlKlSsV57onI9tgDREQJsmjRIjUsJD0e0qthXFdZVm7es2cPChUqZDr2gw8+MF2vXLlytPuGDRumVib38vIy7f/4449NQ01yTLVq1VSvj1yPz+v37dsXKVK8Oun11KlTaphKenHkuYxDabYQ9TXnzp2reqqkB0zaKpu0VXqfpHdMhs+IKPkwACKiBJEhLEmU/u233yz216hRA5kyZbLYZ/7l7u3tHe0+CQiiBkB58+a1eB4ZCrty5Uq8Xz9HjhwWt0NCQtQQleTiyLCVBCJXr17F3bt3YStRX1PaK8Fc1Pa2adNGzThjAESUvBgAEVGCpE+fXgUaxmTipCA9I1FvG5Oc4/P6bm5uFrfnzZunkpEl6JG8HjFp0iSVPxQbY4+OJEmblwuw5jWlvfXr18eoUaPibC8RJT0WQiSiBJEv87///ht79+612B8UFKRmQtmC+TR8CX62b9+OqlWrJvr1ZdgsT548puBHep9kJpc5Y5KzeYAjCdJCkpuNpE3WkPbOnz9f9T6Zu3HjhlWPJyLbYg8QESWIfKF369ZNzdTq168f8ufPr6a5S9AiM6eMAURiSN6MDHMVLlwYM2fOVK9hnCmWmNeXKepjx47FJ598omZ3SfBz8uRJ0/CcyJAhg5qa7+/vj2bNmiF79uyoVauW2rp27apyfCR4WbhwoVU/y+jRo9VstwoVKqh2y6y1/fv3q2BMptYTUfJiDxARWUUChZYtW1rsk6nrS5cuRXBwMPbt26fyXiSR2FjzRqaPS45LmjRpLAIL2Wde/0Zq4sg+Yx0eI5mKLoGHTDGX1965c6cpKTqhry+kLs+uXbtUcCWPa9q0qcrNad68ucVx69evVzWC1q5dqxKrhVSHlmnu8joSxEhQI6/h6+sb62tmyZIFAQEBKliT5GspsNiwYUM19Z6Ikh8LIRKR3ZEhJpnFJXk6UROhiYhsgT1ARERE5HIYABGR3YlpGImIyFY4BEZEREQuhz1ARERE5HIYABEREZHLYQBERERELocBEBEREbkcBkBERETkchgAERERkcthAEREREQuhwEQERERuRwGQERERARX83+WV4Dbsssi8QAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "\n", "# Scatter plot of observed frequencies\n", "plt.scatter(data['Temperature'], data['Frequency'], label=\"Observed\")\n", "\n", "# Logistic curve from the fitted model\n", "x = np.linspace(30, 90, 200)\n", "X = np.column_stack([np.ones_like(x), x])\n", "y = logmodel.predict(X)\n", "\n", "plt.plot(x, y, color='red', label=\"Logistic fit\")\n", "\n", "plt.xlabel(\"Temperature\")\n", "plt.ylabel(\"Failure frequency\")\n", "plt.xlim(30, 90)\n", "plt.ylim(0, 1)\n", "plt.legend()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "trusted": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== SYSTEM ===\n", "OS: macOS-12.5-arm64-arm-64bit\n", "Python: 3.11.15\n", "\n", "=== LIB VERSIONS ===\n", "numpy: 2.4.6\n", "pandas: 3.0.5\n", "matplotlib: 3.11.1\n", "statsmodels: 0.14.6\n", "\n", "=== MODEL RESULTS (GLM Binomial Logit) ===\n", "params:\n", " Intercept 5.084977\n", "Temperature -0.115601\n", "dtype: float64\n", "\n", "std errors:\n", " Intercept 3.052484\n", "Temperature 0.047024\n", "dtype: float64\n", "\n", "AIC: 51.05194634396779\n", "Deviance: 18.08632674249745\n", "Null deviance: 24.230361814971374\n" ] } ], "source": [ "import sys, platform\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "import statsmodels\n", "\n", "print(\"=== SYSTEM ===\")\n", "print(\"OS:\", platform.platform())\n", "print(\"Python:\", sys.version.split()[0])\n", "\n", "print(\"\\n=== LIB VERSIONS ===\")\n", "print(\"numpy:\", np.__version__)\n", "print(\"pandas:\", pd.__version__)\n", "print(\"matplotlib:\", matplotlib.__version__)\n", "print(\"statsmodels:\", statsmodels.__version__)\n", "\n", "print(\"\\n=== MODEL RESULTS (GLM Binomial Logit) ===\")\n", "print(\"params:\\n\", logmodel.params)\n", "print(\"\\nstd errors:\\n\", logmodel.bse)\n", "print(\"\\nAIC:\", logmodel.aic)\n", "print(\"Deviance:\", logmodel.deviance)\n", "print(\"Null deviance:\", logmodel.null_deviance)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "trusted": true }, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Hide code", "kernelspec": { "display_name": "python", "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.11.15" } }, "nbformat": 4, "nbformat_minor": 4 }