{
  "metadata": {
    "kernelspec": {
      "name": "python",
      "display_name": "Python (Pyodide)",
      "language": "python"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "python",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8"
    },
    "celltoolbar": "Hide code"
  },
  "nbformat_minor": 4,
  "nbformat": 4,
  "cells": [
    {
      "cell_type": "markdown",
      "source": "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure",
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "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\nOn 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.",
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": "## Technical information on the computer on which the analysis is run",
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries.",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "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)\")\ndef print_sys_info():\n    import sys\n    import platform\n    print(sys.version)\n    print(platform.uname())\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport statsmodels.api as sm\nsns = None  # seaborn not available in this environment\n\n\nprint_sys_info()\nprint_imported_modules()",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "3.13.2 (main, Oct 20 2025, 18:07:39) [Clang 21.0.0git (https:/github.com/llvm/llvm-project 2f05451198e2f222ec66cec489\nuname_result(system='Emscripten', node='emscripten', release='4.0.9', version='#1', machine='wasm32')\nIPython 9.0.2\nIPython.core.release 9.0.2\nIPython.external.pickleshare 0.7.5\nPIL 11.3.0\nPIL.Image 11.3.0\nPIL._deprecate 11.3.0\nPIL._version 11.3.0\n_ctypes 1.1.0\n_decimal 1.70\nargparse 1.1\ncomm 0.2.3\ncsv 1.0\nctypes 1.1.0\ncycler 0.12.1\ndateutil 2.9.0.post0\ndateutil._version 2.9.0.post0\ndecimal 1.70\ndecorator 5.2.1\nexecuting 2.2.0\nexecuting.version 2.2.0\nipaddress 1.0\njedi 0.19.2\njson 2.0.9\nkiwisolver 1.4.8\nkiwisolver._cext 1.4.8\nlogging 0.5.1.2\nmatplotlib 3.8.4\nmatplotlib._version 3.8.4\nmicropip 0.11.0\nmicropip._vendored.packaging.src.packaging 24.2\nmicropip._version 0.11.0\nnumpy 2.2.5\nnumpy._core 2.2.5\nnumpy._core._multiarray_umath 3.1\n"
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": "<ipython-input-1-4e2079048fa1>: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<ipython-input-1-4e2079048fa1>: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"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "numpy.core 2.2.5\nnumpy.f2py <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.auxfuncs <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.capi_maps <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.cb_rules <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.cfuncs <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.common_rules <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.crackfortran <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.f2py2e <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.f90mod_rules  1.27 \nnumpy.f2py.rules <module 'numpy.f2py.__version__' from '/lib/python3.13/site-packages/numpy/f2py/__version__.py'>\nnumpy.f2py.use_rules  1.3 \nnumpy.linalg._umath_linalg 0.1.5\nnumpy.version 2.2.5\npackaging 24.2\npandas 2.3.2\npandas._version_meson 2.3.2\nparso 0.8.4\npatsy 1.0.1\npatsy.version 1.0.1\npiplite 0.7.0\nplatform 1.0.8\nprompt_toolkit 3.0.50\npure_eval 0.2.3\npure_eval.version 0.2.3\npygments 2.19.1\npyodide 0.29.0\npyodide_kernel 0.7.0\npyparsing 3.2.1\npytz 2025.2\nre 2.2.1\nscipy 1.14.1\nscipy._lib._uarray 0.8.8.dev0+aa94c5a4.scipy\nscipy._lib.array_api_compat 1.5.1\nscipy._lib.array_api_compat.numpy 2.2.5\nscipy._lib.decorator 4.0.5\nscipy.integrate._dop 2.2.5\nscipy.integrate._lsoda 2.2.5\nscipy.integrate._vode 2.2.5\nscipy.interpolate._dfitpack 2.2.5\nscipy.linalg._fblas 2.2.5\nscipy.linalg._flapack 2.2.5\nscipy.optimize._cobyla 2.2.5\nscipy.optimize._lbfgsb 2.2.5\nscipy.optimize._slsqp 2.2.5\nscipy.sparse.linalg._eigen.arpack._arpack 2.2.5\nscipy.sparse.linalg._propack._cpropack 2.2.5\nscipy.sparse.linalg._propack._dpropack 2.2.5\nscipy.sparse.linalg._propack._spropack 2.2.5\nscipy.sparse.linalg._propack._zpropack 2.2.5\nscipy.stats._mvn 2.2.5\nsix 1.17.0\nsocketserver 0.4\nstack_data 0.6.3\nstack_data.version 0.6.3\nstatsmodels 0.14.4\nstatsmodels.__init__ 0.14.4\nstatsmodels._version 0.14.4\nstatsmodels.api 0.14.4\ntraitlets 5.14.3\ntraitlets._version 5.14.3\nurllib.request 3.13\nwcwidth 0.2.13\nzlib 1.0\n"
        }
      ],
      "execution_count": 1
    },
    {
      "cell_type": "markdown",
      "source": "## Loading and inspecting data\nLet's start by reading data.",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "import os\nprint(\"Current working directory:\", os.getcwd())\nprint(\"Files here:\", os.listdir(\".\")[:50])\n",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "Current working directory: /drive\nFiles here: ['Untitled.ipynb', 'data_shuttle.csv', 'src_Python3_challenger__1_.ipynb', 'README.md', 'data', 'notebooks']\n"
        }
      ],
      "execution_count": 7
    },
    {
      "cell_type": "code",
      "source": "data = pd.read_csv(\"data_shuttle.csv\")\n\ndata",
      "metadata": {
        "trusted": true,
        "scrolled": true
      },
      "outputs": [
        {
          "execution_count": 8,
          "output_type": "execute_result",
          "data": {
            "text/plain": "         Date  Count  Temperature  Pressure  Malfunction\n0     4/12/81      6           66        50            0\n1    11/12/81      6           70        50            1\n2     3/22/82      6           69        50            0\n3    11/11/82      6           68        50            0\n4     4/04/83      6           67        50            0\n5     6/18/82      6           72        50            0\n6     8/30/83      6           73       100            0\n7    11/28/83      6           70       100            0\n8     2/03/84      6           57       200            1\n9     4/06/84      6           63       200            1\n10    8/30/84      6           70       200            1\n11   10/05/84      6           78       200            0\n12   11/08/84      6           67       200            0\n13    1/24/85      6           53       200            2\n14    4/12/85      6           67       200            0\n15    4/29/85      6           75       200            0\n16    6/17/85      6           70       200            0\n17  7/2903/85      6           81       200            0\n18    8/27/85      6           76       200            0\n19   10/03/85      6           79       200            0\n20   10/30/85      6           75       200            2\n21   11/26/85      6           76       200            0\n22    1/12/86      6           58       200            1",
            "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Date</th>\n      <th>Count</th>\n      <th>Temperature</th>\n      <th>Pressure</th>\n      <th>Malfunction</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>0</th>\n      <td>4/12/81</td>\n      <td>6</td>\n      <td>66</td>\n      <td>50</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>1</th>\n      <td>11/12/81</td>\n      <td>6</td>\n      <td>70</td>\n      <td>50</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>2</th>\n      <td>3/22/82</td>\n      <td>6</td>\n      <td>69</td>\n      <td>50</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>3</th>\n      <td>11/11/82</td>\n      <td>6</td>\n      <td>68</td>\n      <td>50</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>4</th>\n      <td>4/04/83</td>\n      <td>6</td>\n      <td>67</td>\n      <td>50</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>5</th>\n      <td>6/18/82</td>\n      <td>6</td>\n      <td>72</td>\n      <td>50</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>6</th>\n      <td>8/30/83</td>\n      <td>6</td>\n      <td>73</td>\n      <td>100</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>7</th>\n      <td>11/28/83</td>\n      <td>6</td>\n      <td>70</td>\n      <td>100</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>8</th>\n      <td>2/03/84</td>\n      <td>6</td>\n      <td>57</td>\n      <td>200</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>9</th>\n      <td>4/06/84</td>\n      <td>6</td>\n      <td>63</td>\n      <td>200</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>10</th>\n      <td>8/30/84</td>\n      <td>6</td>\n      <td>70</td>\n      <td>200</td>\n      <td>1</td>\n    </tr>\n    <tr>\n      <th>11</th>\n      <td>10/05/84</td>\n      <td>6</td>\n      <td>78</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>12</th>\n      <td>11/08/84</td>\n      <td>6</td>\n      <td>67</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>13</th>\n      <td>1/24/85</td>\n      <td>6</td>\n      <td>53</td>\n      <td>200</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>14</th>\n      <td>4/12/85</td>\n      <td>6</td>\n      <td>67</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>15</th>\n      <td>4/29/85</td>\n      <td>6</td>\n      <td>75</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>16</th>\n      <td>6/17/85</td>\n      <td>6</td>\n      <td>70</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>17</th>\n      <td>7/2903/85</td>\n      <td>6</td>\n      <td>81</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>18</th>\n      <td>8/27/85</td>\n      <td>6</td>\n      <td>76</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>19</th>\n      <td>10/03/85</td>\n      <td>6</td>\n      <td>79</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>20</th>\n      <td>10/30/85</td>\n      <td>6</td>\n      <td>75</td>\n      <td>200</td>\n      <td>2</td>\n    </tr>\n    <tr>\n      <th>21</th>\n      <td>11/26/85</td>\n      <td>6</td>\n      <td>76</td>\n      <td>200</td>\n      <td>0</td>\n    </tr>\n    <tr>\n      <th>22</th>\n      <td>1/12/86</td>\n      <td>6</td>\n      <td>58</td>\n      <td>200</td>\n      <td>1</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
          },
          "metadata": {}
        }
      ],
      "execution_count": 8
    },
    {
      "cell_type": "markdown",
      "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.",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "%matplotlib inline\npd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\nimport matplotlib.pyplot as plt\n\ndata[\"Frequency\"]=data.Malfunction/data.Count\ndata.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\nplt.grid(True)",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 640x480 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwQ0lEQVR4nO3de1yUdf7//+fIUVyhREVMRExNS7cSt8JTWiuu2sHcLdM1PGAbN01UtNJ1Nw9ZdlhZO3moVDLLjx2sj/WjlMrU0G3TxE7+zDxhCkuiBcpHGOH6/uE6t0ZQYK6RGd497rcbt7re1/uaeV0vLsen13XNjMOyLEsAAACGaODrAgAAALyJcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjOLTcLNp0ybdeuutatmypRwOh955551qt9m4caPi4+MVGhqqtm3bavHixRe/UAAAUG/4NNycPHlSV199tZ577rkazd+/f78GDhyoXr16aceOHfrrX/+q1NRUvfXWWxe5UgAAUF84/OWLMx0Oh95++20NHjz4vHMeeughrV27Vrt27XKNpaSkaOfOndq6dWsdVAkAAPxdoK8LqI2tW7cqMTHRbax///5aunSpnE6ngoKCKm1TWlqq0tJS13JFRYWOHTumyMhIORyOi14zAACwz7IsFRcXq2XLlmrQ4MIXnupVuMnPz1dUVJTbWFRUlE6fPq2jR48qOjq60jbz5s3T7Nmz66pEAABwER06dEitWrW64Jx6FW4kVTrbcvaq2vnOwkyfPl1paWmu5Z9//lmtW7fW/v371bhx44tXqA84nU5t2LBBffv2rfIsFqpHD+2hf/bRQ3von33+2sPi4mLFxcXV6O/uehVuWrRoofz8fLexgoICBQYGKjIyssptQkJCFBISUmm8SZMmCg8Pvyh1+orT6VRYWJgiIyP96oCsT+ihPfTPPnpoD/2zz197eLaWmtxSUq8+5yYhIUFZWVluY+vXr1e3bt386hcAAAB8x6fh5sSJE8rJyVFOTo6kM2/1zsnJUW5urqQzl5SSkpJc81NSUnTw4EGlpaVp165dWrZsmZYuXaqpU6f6onwAAOCHfHpZatu2berbt69r+ey9MSNHjlRGRoby8vJcQUeS4uLilJmZqcmTJ+v5559Xy5Yt9cwzz+iPf/xjndcOAAD8k0/DTZ8+fXShj9nJyMioNHbjjTfqiy++uIhVAQCA+qxe3XMDAABQHcINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFF8Hm4WLlyouLg4hYaGKj4+Xps3b77g/FdffVVXX321wsLCFB0drdGjR6uwsLCOqgUAAP7Op+Fm9erVmjRpkmbMmKEdO3aoV69eGjBggHJzc6uc/+mnnyopKUnJycn65ptv9MYbb+jzzz/X2LFj67hyAADgr3wabtLT05WcnKyxY8eqU6dOWrBggWJiYrRo0aIq5//rX/9SmzZtlJqaqri4OPXs2VP33Xeftm3bVseVAwAAfxXoqycuKyvT9u3bNW3aNLfxxMREbdmypcptunfvrhkzZigzM1MDBgxQQUGB3nzzTQ0aNOi8z1NaWqrS0lLXclFRkSTJ6XTK6XR6YU/8x9n9MW2/6hI9tIf+2UcP7aF/9vlrD2tTj8OyLOsi1nJeR44c0WWXXabs7Gx1797dNf7YY4/p5Zdf1u7du6vc7s0339To0aN16tQpnT59WrfddpvefPNNBQUFVTl/1qxZmj17dqXx1157TWFhYd7ZGQAAcFGVlJRo+PDh+vnnnxUeHn7BuT47c3OWw+FwW7Ysq9LYWd9++61SU1P18MMPq3///srLy9MDDzyglJQULV26tMptpk+frrS0NNdyUVGRYmJilJiYWG1z6hun06msrCz169fvvGEPF0YP7aF/9tFDe+ifff7aw7NXXmrCZ+GmadOmCggIUH5+vtt4QUGBoqKiqtxm3rx56tGjhx544AFJ0m9/+1s1atRIvXr10ty5cxUdHV1pm5CQEIWEhFQaDwoK8qtfmjeZvG91hR7aQ//so4f20D/7/K2HtanFZzcUBwcHKz4+XllZWW7jWVlZbpepfqmkpEQNGriXHBAQIOnMGR8AAACfvlsqLS1NL730kpYtW6Zdu3Zp8uTJys3NVUpKiqQzl5SSkpJc82+99VatWbNGixYt0r59+5Sdna3U1FRdd911atmypa92AwAA+BGf3nMzdOhQFRYWas6cOcrLy1Pnzp2VmZmp2NhYSVJeXp7bZ96MGjVKxcXFeu655zRlyhRdcskluummm/TEE0/4ahcAAICf8fkNxePGjdO4ceOqXJeRkVFpbMKECZowYcJFrgoAANRXPv/6BQAAAG8i3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYxefhZuHChYqLi1NoaKji4+O1efPmC84vLS3VjBkzFBsbq5CQEF1++eVatmxZHVULAAD8XaAvn3z16tWaNGmSFi5cqB49emjJkiUaMGCAvv32W7Vu3brKbe666y795z//0dKlS9WuXTsVFBTo9OnTdVw5AADwVz4NN+np6UpOTtbYsWMlSQsWLNC6deu0aNEizZs3r9L8Dz74QBs3btS+ffvUpEkTSVKbNm3qsmQAAODnfBZuysrKtH37dk2bNs1tPDExUVu2bKlym7Vr16pbt2568skn9corr6hRo0a67bbb9Mgjj6hhw4ZVblNaWqrS0lLXclFRkSTJ6XTK6XR6aW/8w9n9MW2/6hI9tIf+2UcP7aF/9vlrD2tTj8/CzdGjR1VeXq6oqCi38aioKOXn51e5zb59+/Tpp58qNDRUb7/9to4ePapx48bp2LFj573vZt68eZo9e3al8fXr1yssLMz+jvihrKwsX5dQ79FDe+ifffTQHvpnn7/1sKSkpMZzfXpZSpIcDofbsmVZlcbOqqiokMPh0KuvvqqIiAhJZy5t/elPf9Lzzz9f5dmb6dOnKy0tzbVcVFSkmJgYJSYmKjw83It74ntOp1NZWVnq16+fgoKCfF1OvUQP7aF/9tFDe+ifff7aw7NXXmrCZ+GmadOmCggIqHSWpqCgoNLZnLOio6N12WWXuYKNJHXq1EmWZemHH35Q+/btK20TEhKikJCQSuNBQUF+9UvzJpP3ra7QQ3von3300B76Z5+/9bA2tfjsreDBwcGKj4+vdNorKytL3bt3r3KbHj166MiRIzpx4oRr7LvvvlODBg3UqlWri1ovAACoHzwKN/v37/fKk6elpemll17SsmXLtGvXLk2ePFm5ublKSUmRdOaSUlJSkmv+8OHDFRkZqdGjR+vbb7/Vpk2b9MADD2jMmDHnvaEYAAD8unh0Wapdu3bq3bu3kpOT9ac//UmhoaEePfnQoUNVWFioOXPmKC8vT507d1ZmZqZiY2MlSXl5ecrNzXXN/81vfqOsrCxNmDBB3bp1U2RkpO666y7NnTvXo+cHAADm8Sjc7Ny5U8uWLdOUKVN0//33a+jQoUpOTtZ1111X68caN26cxo0bV+W6jIyMSmMdO3b0uzu4AQCA//DoslTnzp2Vnp6uw4cPa/ny5crPz1fPnj111VVXKT09XT/++KO36wQAAKgRWzcUBwYG6o477tDrr7+uJ554Qnv37tXUqVPVqlUrJSUlKS8vz1t1AgAA1IitcLNt2zaNGzdO0dHRSk9P19SpU7V37159/PHHOnz4sG6//XZv1QkAAFAjHt1zk56eruXLl2v37t0aOHCgVqxYoYEDB6pBgzNZKS4uTkuWLFHHjh29WiwAAEB1PAo3ixYt0pgxYzR69Gi1aNGiyjmtW7fW0qVLbRUHAABQWx6Fmz179lQ7Jzg4WCNHjvTk4QEAADzm0T03y5cv1xtvvFFp/I033tDLL79suygAAABPeRRuHn/8cTVt2rTSePPmzfXYY4/ZLgoAAMBTHoWbgwcPKi4urtJ4bGys2ycKAwAA1DWPwk3z5s315ZdfVhrfuXOnIiMjbRcFAADgKY/Czd13363U1FRt2LBB5eXlKi8v18cff6yJEyfq7rvv9naNAAAANebRu6Xmzp2rgwcP6uabb1Zg4JmHqKioUFJSEvfcAAAAn/Io3AQHB2v16tV65JFHtHPnTjVs2FBdunRxfZs3AACAr3gUbs7q0KGDOnTo4K1aAAAAbPMo3JSXlysjI0MfffSRCgoKVFFR4bb+448/9kpxAAAAteVRuJk4caIyMjI0aNAgde7cWQ6Hw9t1AQAAeMSjcPM///M/ev311zVw4EBv1wMAAGCLR28FDw4OVrt27bxdCwAAgG0ehZspU6bo6aeflmVZ3q4HAADAFo8uS3366afasGGD3n//fV111VUKCgpyW79mzRqvFAcAAFBbHoWbSy65RHfccYe3awEAALDNo3CzfPlyb9cBAADgFR7dcyNJp0+f1ocffqglS5aouLhYknTkyBGdOHHCa8UBAADUlkdnbg4ePKg//OEPys3NVWlpqfr166fGjRvrySef1KlTp7R48WJv1wkAAFAjHp25mThxorp166bjx4+rYcOGrvE77rhDH330kdeKAwAAqC2P3y2VnZ2t4OBgt/HY2FgdPnzYK4UBAAB4wqMzNxUVFSovL680/sMPP6hx48a2iwIAAPCUR+GmX79+WrBggWvZ4XDoxIkTmjlzJl/JAAAAfMqjy1L//Oc/1bdvX1155ZU6deqUhg8frj179qhp06ZatWqVt2sEAACoMY/CTcuWLZWTk6NVq1bpiy++UEVFhZKTk/XnP//Z7QZjAACAuuZRuJGkhg0basyYMRozZow36wEAALDFo3CzYsWKC65PSkryqBgAAAC7PAo3EydOdFt2Op0qKSlRcHCwwsLCCDcAAMBnPHq31PHjx91+Tpw4od27d6tnz57cUAwAAHzK4++WOlf79u31+OOPVzqrAwAAUJe8Fm4kKSAgQEeOHPHmQwIAANSKR/fcrF271m3Zsizl5eXpueeeU48ePbxSGAAAgCc8CjeDBw92W3Y4HGrWrJluuukmzZ8/3xt1AQAAeMSjcFNRUeHtOgAAALzCq/fcAAAA+JpHZ27S0tJqPDc9Pd2TpwAAAPCIR+Fmx44d+uKLL3T69GldccUVkqTvvvtOAQEB6tq1q2uew+HwTpUAAAA15FG4ufXWW9W4cWO9/PLLuvTSSyWd+WC/0aNHq1evXpoyZYpXiwQAAKgpj+65mT9/vubNm+cKNpJ06aWXau7cubxbCgAA+JRH4aaoqEj/+c9/Ko0XFBSouLjYdlEAAACe8ijc3HHHHRo9erTefPNN/fDDD/rhhx/05ptvKjk5WUOGDPF2jQAAADXm0T03ixcv1tSpUzVixAg5nc4zDxQYqOTkZD311FNeLRAAAKA2PAo3YWFhWrhwoZ566int3btXlmWpXbt2atSokbfrAwAAqBVbH+KXl5envLw8dejQQY0aNZJlWd6qCwAAwCMehZvCwkLdfPPN6tChgwYOHKi8vDxJ0tixY3kbOAAA8CmPws3kyZMVFBSk3NxchYWFucaHDh2qDz74wGvFAQAA1JZH99ysX79e69atU6tWrdzG27dvr4MHD3qlMAAAAE94dObm5MmTbmdszjp69KhCQkJsFwUAAOApj8JN7969tWLFCteyw+FQRUWFnnrqKfXt29drxQEAANSWR5elnnrqKfXp00fbtm1TWVmZHnzwQX3zzTc6duyYsrOzvV0jAABAjXl05ubKK6/Ul19+qeuuu079+vXTyZMnNWTIEO3YsUOXX365t2sEAACosVqfuXE6nUpMTNSSJUs0e/bsi1ETAACAx2p95iYoKEhff/21HA7HxagHAADAFo8uSyUlJWnp0qXergUAAMA2j24oLisr00svvaSsrCx169at0ndKpaene6U4AACA2qpVuNm3b5/atGmjr7/+Wl27dpUkfffdd25zuFwFAAB8qVbhpn379srLy9OGDRsknfm6hWeeeUZRUVEXpTgAAIDaqtU9N+d+6/f777+vkydPerUgAAAAOzy6ofisc8OOJxYuXKi4uDiFhoYqPj5emzdvrtF22dnZCgwM1DXXXGO7BgAAYI5ahRuHw1Hpnho799isXr1akyZN0owZM7Rjxw716tVLAwYMUG5u7gW3+/nnn5WUlKSbb77Z4+cGAABmqtU9N5ZladSoUa4vxzx16pRSUlIqvVtqzZo1NXq89PR0JScna+zYsZKkBQsWaN26dVq0aJHmzZt33u3uu+8+DR8+XAEBAXrnnXdqswsAAMBwtQo3I0eOdFseMWKEx09cVlam7du3a9q0aW7jiYmJ2rJly3m3W758ufbu3auVK1dq7ty51T5PaWmpSktLXctFRUWSznzSstPp9LB6/3R2f0zbr7pED+2hf/bRQ3von33+2sPa1FOrcLN8+fJaF3M+R48eVXl5eaV3WkVFRSk/P7/Kbfbs2aNp06Zp8+bNCgysWenz5s2r8msi1q9fr7CwsNoXXg9kZWX5uoR6jx7aQ//so4f20D/7/K2HJSUlNZ7r0Yf4edO59+xYllXlfTzl5eUaPny4Zs+erQ4dOtT48adPn660tDTXclFRkWJiYpSYmKjw8HDPC/dDTqdTWVlZ6tevn4KCgnxdTr1ED+2hf/bRQ3von33+2sOzV15qwmfhpmnTpgoICKh0lqagoKDKz80pLi7Wtm3btGPHDt1///2SpIqKClmWpcDAQK1fv1433XRTpe1CQkJc9wj9UlBQkF/90rzJ5H2rK/TQHvpnHz20h/7Z5289rE0ttt4KbkdwcLDi4+MrnfbKyspS9+7dK80PDw/XV199pZycHNdPSkqKrrjiCuXk5Oj666+vq9IBAIAf8+llqbS0NN1zzz3q1q2bEhIS9MILLyg3N1cpKSmSzlxSOnz4sFasWKEGDRqoc+fObts3b95coaGhlcYBAMCvl0/DzdChQ1VYWKg5c+YoLy9PnTt3VmZmpmJjYyVJeXl51X7mDQAAwC/5/IbicePGady4cVWuy8jIuOC2s2bN0qxZs7xfFAAAqLd8ds8NAADAxUC4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwo0X7fvxhDbsLtD+oyd9XQoA/Cod+O/r78HCEh9XAl8K9HUBJvippEypq3K0ac+PrrHe7Zvp2WHXKiIsyIeVAcCvw9nX4c/2FejJ66RBz27W9W2b8zr8K8WZGy9IXZWj7O+Puo1lf39UE1bt8FFFAPDrwuswfolwY9O+H09o054fVW5ZbuPllqVNe37kEhUAXGS8DuNchBubDh678HXdA4X8oQKAi4nXYZyLcGNTbJOwC65vE9mojioBgF8nXodxLsKNTW2b/Ua92zdTgMPhNh7gcKh3+2aKa8ofKgC4mHgdxrkIN17w7LBr1aNdU7exHu2a6tlh1/qoIgD4deF1GL/EW8G9ICIsSCuSr9P+oyd1oPCk2kQ24l8KAFCHzr4Of5//s7757BP9fxN6qV2LCF+XBR8h3HhRXFNCDQD4UmxkmL7573/x68VlKQAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRfB5uFi5cqLi4OIWGhio+Pl6bN28+79w1a9aoX79+atasmcLDw5WQkKB169bVYbUAAMDf+TTcrF69WpMmTdKMGTO0Y8cO9erVSwMGDFBubm6V8zdt2qR+/fopMzNT27dvV9++fXXrrbdqx44ddVw5AADwVz4NN+np6UpOTtbYsWPVqVMnLViwQDExMVq0aFGV8xcsWKAHH3xQv/vd79S+fXs99thjat++vd599906rhwAAPirQF89cVlZmbZv365p06a5jScmJmrLli01eoyKigoVFxerSZMm551TWlqq0tJS13JRUZEkyel0yul0elC5/zq7P6btV12ih/bQP/vooT30zz5/7WFt6vFZuDl69KjKy8sVFRXlNh4VFaX8/PwaPcb8+fN18uRJ3XXXXeedM2/ePM2ePbvS+Pr16xUWFla7ouuJrKwsX5dQ79FDe+ifffTQHvpnn7/1sKSkpMZzfRZuznI4HG7LlmVVGqvKqlWrNGvWLP3v//6vmjdvft5506dPV1pammu5qKhIMTExSkxMVHh4uOeF+yGn06msrCz169dPQUFBvi6nXqKH9tA/++ihPfTPPn/t4dkrLzXhs3DTtGlTBQQEVDpLU1BQUOlszrlWr16t5ORkvfHGG/r9739/wbkhISEKCQmpNB4UFORXvzRvMnnf6go9tIf+2UcP7aF/9vlbD2tTi89uKA4ODlZ8fHyl015ZWVnq3r37ebdbtWqVRo0apddee02DBg262GUCAIB6xqeXpdLS0nTPPfeoW7duSkhI0AsvvKDc3FylpKRIOnNJ6fDhw1qxYoWkM8EmKSlJTz/9tG644QbXWZ+GDRsqIiLCZ/sBAAD8h0/DzdChQ1VYWKg5c+YoLy9PnTt3VmZmpmJjYyVJeXl5bp95s2TJEp0+fVrjx4/X+PHjXeMjR45URkZGXZcPAAD8kM9vKB43bpzGjRtX5bpzA8snn3xy8QsCAAD1ms+/fgEAAMCbCDfwmn0/ntCG3QXaf/SkT7YHPTRB9p4fJUlb9x71cSVA/eXzy1Ko/34qKVPqqhxt+u+LsiT1bt9Mzw67VhFh1b91z+72oIcmOFh4UoOfz1ZJaZmevE6695XtCgsJ1trxPRUTaeYHjgIXC2duYFvqqhxlf+/+r8zs749qwqqafaGp3e1BD00w+PlsHS9x/3j54yVO3fb8pz6qCKi/CDewZd+PJ7Rpz48qtyy38XLL0qY9P1Z7ecTu9qCHJti4u6BSsDnreIlTm39xRg5A9Qg3sOXgsQt/18eBwgv/xWp3e9BDE+T88NMF13+Re7xuCgEMQbiBLbFNLnwvQJvIRhd1e9BDE1zT6pILru/a+tK6KQQwBOEGtrRt9hv1bt9MAed82WmAw6He7ZsprumF/2K1uz3ooQluvKK5Lj3Pjd+XhgWpV/tmdVwRUL8RbmDbs8OuVY92Td3GerRrqmeHXVsn24MemmDt+J6VAs6lYUFaO76njyoC6i/eCg7bIsKCtCL5Ou0/elIHCk+qTWSjWp0tsLs96KEJYiLDtOPhRG36//N0fPe/9eI98erdMdrXZQH1EuEGXhPX1N5fqHa3Bz00QcLlTZW5+8x/AXiGy1IAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABjF5+Fm4cKFiouLU2hoqOLj47V58+YLzt+4caPi4+MVGhqqtm3bavHixXVUKYCLbd+PJ7Rhd4H2Hz3p0far/52rSat36I1th3xWg93ts/f8KEnauveoR9vb5ev9t+vAf5/3YGGJR9tv3F2gpz/6Tpv/+3vwhfr+O/D180tSoM+eWdLq1as1adIkLVy4UD169NCSJUs0YMAAffvtt2rdunWl+fv379fAgQN17733auXKlcrOzta4cePUrFkz/fGPf/TBHgDwhp9KypS6KkebfvEXSu/2zfTssGsVERZU7fZf/fCT7li4RacrLEnSOzuOaPqar7R2fA9deVlEndRgd/uDhSc1+PlslZSW6cnrpHtf2a6wkGCtHd9TMZFhNdoHO3y9/3adff7P9hXoyeukQc9u1vVtm9e6/8dLnK6xS8OC6qz/kjm/A189/y/59MxNenq6kpOTNXbsWHXq1EkLFixQTEyMFi1aVOX8xYsXq3Xr1lqwYIE6deqksWPHasyYMfrHP/5Rx5UD8KbUVTnK/t79TEX290c1YdWOGm3/y2Bz1ukKS7c9n11nNdjd/ty/WCXpeIlTtz3/aY22t8vX+29Xfe+/xO/Am3x25qasrEzbt2/XtGnT3MYTExO1ZcuWKrfZunWrEhMT3cb69++vpUuXyul0KiiocjIsLS1VaWmpa/nnn3+WJB07dkxOp7PS/PrM6XSqpKREhYWFVfYC1aOH9njSv9zCk9qy64AcqvyCtGXXCe3YE63WTc7/L+e1OYelshPnfTF7+eOvdMvVLS9qDXa3/9feoyou+kmBkgIrLJWUVCjQ2UDlFQ4VF0nvb/tO18VFXnAf7PD1/tvl9vwN3PtX2/6fqy76L/nX78AXf45rori4WJJkWVY1M89M8onDhw9bkqzs7Gy38UcffdTq0KFDldu0b9/eevTRR93GsrOzLUnWkSNHqtxm5syZliR++OGHH3744ceAn0OHDlWbMXx6z40kORwOt2XLsiqNVTe/qvGzpk+frrS0NNdyRUWFjh07psjIyAs+T31UVFSkmJgYHTp0SOHh4b4up16ih/bQP/vooT30zz5/7aFlWSouLlbLlhc+Eyv58LJU06ZNFRAQoPz8fLfxgoICRUVFVblNixYtqpwfGBioyMiqTxmGhIQoJCTEbeySSy7xvPB6IDw83K8OyPqIHtpD/+yjh/bQP/v8sYcRERE1muezG4qDg4MVHx+vrKwst/GsrCx17969ym0SEhIqzV+/fr26devG/REAAECSj98tlZaWppdeeknLli3Trl27NHnyZOXm5iolJUXSmUtKSUlJrvkpKSk6ePCg0tLStGvXLi1btkxLly7V1KlTfbULAADAz/j0npuhQ4eqsLBQc+bMUV5enjp37qzMzEzFxsZKkvLy8pSbm+uaHxcXp8zMTE2ePFnPP/+8WrZsqWeeeYbPuPmvkJAQzZw5s9JlONQcPbSH/tlHD+2hf/aZ0EOHZdXkPVUAAAD1g8+/fgEAAMCbCDcAAMAohBsAAGAUwg0AADAK4aaemTVrlhwOh9tPixYtXOtHjRpVaf0NN9zgw4r90+HDhzVixAhFRkYqLCxM11xzjbZv3+5ab1mWZs2apZYtW6phw4bq06ePvvnmGx9W7F+q6x/H4YW1adOmUn8cDofGjx8vieOvOtX1j+OveqdPn9bf/vY3xcXFqWHDhmrbtq3mzJmjiooK15z6fBz6/OsXUHtXXXWVPvzwQ9dyQECA2/o//OEPWr58uWs5ODi4zmqrD44fP64ePXqob9++ev/999W8eXPt3bvX7ZOrn3zySaWnpysjI0MdOnTQ3Llz1a9fP+3evVuNGzf2XfF+oCb9kzgOL+Tzzz9XeXm5a/nrr79Wv379dOedd0ri+KtOdf2TOP6q88QTT2jx4sV6+eWXddVVV2nbtm0aPXq0IiIiNHHiREn1/Dis9tun4FdmzpxpXX311eddP3LkSOv222+vs3rqo4ceesjq2bPneddXVFRYLVq0sB5//HHX2KlTp6yIiAhr8eLFdVGiX6uuf5bFcVhbEydOtC6//HKroqKC488Dv+yfZXH81cSgQYOsMWPGuI0NGTLEGjFihGVZ9f91kMtS9dCePXvUsmVLxcXF6e6779a+ffvc1n/yySdq3ry5OnTooHvvvVcFBQU+qtQ/rV27Vt26ddOdd96p5s2b69prr9WLL77oWr9//37l5+crMTHRNRYSEqIbb7xRW7Zs8UXJfqW6/p3FcVgzZWVlWrlypcaMGSOHw8HxV0vn9u8sjr8L69mzpz766CN99913kqSdO3fq008/1cCBAyXV/9dBwk09c/3112vFihVat26dXnzxReXn56t79+4qLCyUJA0YMECvvvqqPv74Y82fP1+ff/65brrpJpWWlvq4cv+xb98+LVq0SO3bt9e6deuUkpKi1NRUrVixQpJcX8567he4RkVFVfri1l+j6voncRzWxjvvvKOffvpJo0aNksTxV1vn9k/i+KuJhx56SMOGDVPHjh0VFBSka6+9VpMmTdKwYcMkGXAc+vrUEew5ceKEFRUVZc2fP7/K9UeOHLGCgoKst956q44r819BQUFWQkKC29iECROsG264wbIsy8rOzrYkWUeOHHGbM3bsWKt///51Vqe/qq5/VeE4PL/ExETrlltucS1z/NXOuf2rCsdfZatWrbJatWplrVq1yvryyy+tFStWWE2aNLEyMjIsy6r/xyFnbuq5Ro0aqUuXLtqzZ0+V66OjoxUbG3ve9b9G0dHRuvLKK93GOnXq5Poes7PvPjv3XycFBQWV/hXza1Rd/863DcdhZQcPHtSHH36osWPHusY4/mquqv5VheOvsgceeEDTpk3T3XffrS5duuiee+7R5MmTNW/ePEn1/zgk3NRzpaWl2rVrl6Kjo6tcX1hYqEOHDp13/a9Rjx49tHv3brex7777zvWFrXFxcWrRooWysrJc68vKyrRx40Z17969Tmv1R9X1ryoch1Vbvny5mjdvrkGDBrnGOP5qrqr+VYXjr7KSkhI1aOAeAQICAlxvBa/3x6GvTx2hdqZMmWJ98skn1r59+6x//etf1i233GI1btzYOnDggFVcXGxNmTLF2rJli7V//35rw4YNVkJCgnXZZZdZRUVFvi7db/z73/+2AgMDrUcffdTas2eP9eqrr1phYWHWypUrXXMef/xxKyIiwlqzZo311VdfWcOGDbOio6Ppo1V9/zgOa6a8vNxq3bq19dBDD1Vax/FXvfP1j+OvZkaOHGlddtll1nvvvWft37/fWrNmjdW0aVPrwQcfdM2pz8ch4aaeGTp0qBUdHW0FBQVZLVu2tIYMGWJ98803lmVZVklJiZWYmGg1a9bMCgoKslq3bm2NHDnSys3N9XHV/ufdd9+1OnfubIWEhFgdO3a0XnjhBbf1FRUV1syZM60WLVpYISEhVu/eva2vvvrKR9X6nwv1j+OwZtatW2dJsnbv3l1pHcdf9c7XP46/mikqKrImTpxotW7d2goNDbXatm1rzZgxwyotLXXNqc/HocOyLMvXZ48AAAC8hXtuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AXJDD4bjgz6hRo3xdotf16dNHkyZN8nUZADwU6OsCAPi3vLw81/+vXr1aDz/8sNsXZzZs2NAXZXnE6XQqKCjI2OcDcAZnbgBcUIsWLVw/ERERcjgcbmObNm1SfHy8QkND1bZtW82ePVunT592be9wOLRkyRLdcsstCgsLU6dOnbR161Z9//336tOnjxo1aqSEhATt3bvXtc2sWbN0zTXXaMmSJYqJiVFYWJjuvPNO/fTTT261LV++XJ06dVJoaKg6duyohQsXutYdOHBADodDr7/+uvr06aPQ0FCtXLlShYWFGjZsmFq1aqWwsDB16dJFq1atcm03atQobdy4UU8//bTr7NSBAweUkZGhSy65xO3533nnHTkcjkp1L1u2TG3btlVISIgsy9LPP/+sv/zlL2revLnCw8N10003aefOnV76DQE4F+EGgMfWrVunESNGKDU1Vd9++62WLFmijIwMPfroo27zHnnkESUlJSknJ0cdO3bU8OHDdd9992n69Onatm2bJOn+++932+b777/X66+/rnfffVcffPCBcnJyNH78eNf6F198UTNmzNCjjz6qXbt26bHHHtPf//53vfzyy26P89BDDyk1NVW7du1S//79derUKcXHx+u9997T119/rb/85S+655579Nlnn0mSnn76aSUkJOjee+9VXl6e8vLyFBMTU+OenK37rbfeUk5OjiRp0KBBys/PV2ZmprZv366uXbvq5ptv1rFjx2r8uABqwcdf3AmgHlm+fLkVERHhWu7Vq5f12GOPuc155ZVXrOjoaNeyJOtvf/uba3nr1q2WJGvp0qWusVWrVlmhoaGu5ZkzZ1oBAQHWoUOHXGPvv/++1aBBAysvL8+yLMuKiYmxXnvtNbfnfuSRR6yEhATLsixr//79liRrwYIF1e7XwIEDrSlTpriWb7zxRmvixIkX3HfLsqy3337b+uXL6MyZM62goCCroKDANfbRRx9Z4eHh1qlTp9y2vfzyy60lS5ZUWxuA2uOeGwAe2759uz7//HO3MzXl5eU6deqUSkpKFBYWJkn67W9/61ofFRUlSerSpYvb2KlTp1RUVKTw8HBJUuvWrdWqVSvXnISEBFVUVGj37t0KCAjQoUOHlJycrHvvvdc15/Tp04qIiHCrsVu3bm7L5eXlevzxx7V69WodPnxYpaWlKi0tVaNGjey2Q5IUGxurZs2auZa3b9+uEydOKDIy0m3e//3f/7ldigPgPYQbAB6rqKjQ7NmzNWTIkErrQkNDXf//y5tqz96jUtVYRUXFeZ/r7ByHw+Ga9+KLL+r66693mxcQEOC2fG5omT9/vv75z39qwYIF6tKlixo1aqRJkyaprKzs/DsqqUGDBrIsy23M6XRWmnfu81VUVCg6OlqffPJJpbnn3sMDwDsINwA81rVrV+3evVvt2rXz+mPn5ubqyJEjatmypSRp69atatCggTp06KCoqChddtll2rdvn/785z/X6nE3b96s22+/XSNGjJB0Jnzs2bNHnTp1cs0JDg5WeXm523bNmjVTcXGxTp486QowZ++puZCuXbsqPz9fgYGBatOmTa1qBeAZwg0Ajz388MO65ZZbFBMTozvvvFMNGjTQl19+qa+++kpz58619dihoaEaOXKk/vGPf6ioqEipqam666671KJFC0ln3pmUmpqq8PBwDRgwQKWlpdq2bZuOHz+utLS08z5uu3bt9NZbb2nLli269NJLlZ6ervz8fLdw06ZNG3322Wc6cOCAfvOb36hJkya6/vrrFRYWpr/+9a+aMGGC/v3vfysjI6Pa/fj973+vhIQEDR48WE888YSuuOIKHTlyRJmZmRo8eHCly2YA7OPdUgA81r9/f7333nvKysrS7373O91www1KT09XbGys7cdu166dhgwZooEDByoxMVGdO3d2e6v32LFj9dJLLykjI0NdunTRjTfeqIyMDMXFxV3wcf/+97+ra9eu6t+/v/r06aMWLVpo8ODBbnOmTp2qgIAAXXnllWrWrJlyc3PVpEkTrVy5UpmZma63j8+aNava/XA4HMrMzFTv3r01ZswYdejQQXfffbcOHDjguv8IgHc5rHMvIgOAj82aNUvvvPNOjS77AMC5OHMDAACMQrgBAABG4bIUAAAwCmduAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBR/h9DigUvvZMUYwAAAABJRU5ErkJggg=="
          },
          "metadata": {}
        }
      ],
      "execution_count": 9
    },
    {
      "cell_type": "markdown",
      "source": "## Logistic regression\n\nLet'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.",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "import statsmodels.api as sm\n\ndata[\"Success\"]=data.Count-data.Malfunction\ndata[\"Intercept\"]=1\n\nlogmodel = sm.GLM(\n    data['Frequency'],\n    data[['Intercept','Temperature']],\n    family=sm.families.Binomial(link=sm.families.links.logit())\n).fit()\n\n\nlogmodel.summary()",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": "/lib/python3.13/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"
        },
        {
          "execution_count": 13,
          "output_type": "execute_result",
          "data": {
            "text/plain": "<class 'statsmodels.iolib.summary.Summary'>\n\"\"\"\n                 Generalized Linear Model Regression Results                  \n==============================================================================\nDep. Variable:              Frequency   No. Observations:                   23\nModel:                            GLM   Df Residuals:                       21\nModel Family:                Binomial   Df Model:                            1\nLink Function:                  logit   Scale:                          1.0000\nMethod:                          IRLS   Log-Likelihood:                -3.9210\nDate:                Fri, 19 Dec 2025   Deviance:                       3.0144\nTime:                        22:15:03   Pearson chi2:                     5.00\nNo. Iterations:                     6   Pseudo R-squ. (CS):            0.04355\nCovariance Type:            nonrobust                                         \n===============================================================================\n                  coef    std err          z      P>|z|      [0.025      0.975]\n-------------------------------------------------------------------------------\nIntercept       5.0850      7.477      0.680      0.496      -9.570      19.740\nTemperature    -0.1156      0.115     -1.004      0.316      -0.341       0.110\n===============================================================================\n\"\"\"",
            "text/html": "<table class=\"simpletable\">\n<caption>Generalized Linear Model Regression Results</caption>\n<tr>\n  <th>Dep. Variable:</th>       <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td> \n</tr>\n<tr>\n  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td> \n</tr>\n<tr>\n  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td> \n</tr>\n<tr>\n  <th>Link Function:</th>         <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n</tr>\n<tr>\n  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -3.9210</td>\n</tr>\n<tr>\n  <th>Date:</th>            <td>Fri, 19 Dec 2025</td> <th>  Deviance:          </th> <td>  3.0144</td>\n</tr>\n<tr>\n  <th>Time:</th>                <td>22:15:03</td>     <th>  Pearson chi2:      </th>  <td>  5.00</td> \n</tr>\n<tr>\n  <th>No. Iterations:</th>          <td>6</td>        <th>  Pseudo R-squ. (CS):</th>  <td>0.04355</td>\n</tr>\n<tr>\n  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n</tr>\n</table>\n<table class=\"simpletable\">\n<tr>\n       <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n</tr>\n<tr>\n  <th>Intercept</th>   <td>    5.0850</td> <td>    7.477</td> <td>    0.680</td> <td> 0.496</td> <td>   -9.570</td> <td>   19.740</td>\n</tr>\n<tr>\n  <th>Temperature</th> <td>   -0.1156</td> <td>    0.115</td> <td>   -1.004</td> <td> 0.316</td> <td>   -0.341</td> <td>    0.110</td>\n</tr>\n</table>",
            "text/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:}            & Fri, 19 Dec 2025 & \\textbf{  Deviance:          } &    3.0144   \\\\\n\\textbf{Time:}            &     22:15:03     & \\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}"
          },
          "metadata": {}
        }
      ],
      "execution_count": 13
    },
    {
      "cell_type": "markdown",
      "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).",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "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\nlogmodel.summary()",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "execution_count": 15,
          "output_type": "execute_result",
          "data": {
            "text/plain": "<class 'statsmodels.iolib.summary.Summary'>\n\"\"\"\n                 Generalized Linear Model Regression Results                  \n==============================================================================\nDep. Variable:              Frequency   No. Observations:                   23\nModel:                            GLM   Df Residuals:                       21\nModel Family:                Binomial   Df Model:                            1\nLink Function:                  Logit   Scale:                          1.0000\nMethod:                          IRLS   Log-Likelihood:                -23.526\nDate:                Fri, 19 Dec 2025   Deviance:                       18.086\nTime:                        22:16:08   Pearson chi2:                     30.0\nNo. Iterations:                     6   Pseudo R-squ. (CS):             0.2344\nCovariance Type:            nonrobust                                         \n===============================================================================\n                  coef    std err          z      P>|z|      [0.025      0.975]\n-------------------------------------------------------------------------------\nIntercept       5.0850      3.052      1.666      0.096      -0.898      11.068\nTemperature    -0.1156      0.047     -2.458      0.014      -0.208      -0.023\n===============================================================================\n\"\"\"",
            "text/html": "<table class=\"simpletable\">\n<caption>Generalized Linear Model Regression Results</caption>\n<tr>\n  <th>Dep. Variable:</th>       <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td> \n</tr>\n<tr>\n  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td> \n</tr>\n<tr>\n  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td> \n</tr>\n<tr>\n  <th>Link Function:</th>         <td>Logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n</tr>\n<tr>\n  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -23.526</td>\n</tr>\n<tr>\n  <th>Date:</th>            <td>Fri, 19 Dec 2025</td> <th>  Deviance:          </th> <td>  18.086</td>\n</tr>\n<tr>\n  <th>Time:</th>                <td>22:16:08</td>     <th>  Pearson chi2:      </th>  <td>  30.0</td> \n</tr>\n<tr>\n  <th>No. Iterations:</th>          <td>6</td>        <th>  Pseudo R-squ. (CS):</th>  <td>0.2344</td> \n</tr>\n<tr>\n  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n</tr>\n</table>\n<table class=\"simpletable\">\n<tr>\n       <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n</tr>\n<tr>\n  <th>Intercept</th>   <td>    5.0850</td> <td>    3.052</td> <td>    1.666</td> <td> 0.096</td> <td>   -0.898</td> <td>   11.068</td>\n</tr>\n<tr>\n  <th>Temperature</th> <td>   -0.1156</td> <td>    0.047</td> <td>   -2.458</td> <td> 0.014</td> <td>   -0.208</td> <td>   -0.023</td>\n</tr>\n</table>",
            "text/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:}            & Fri, 19 Dec 2025 & \\textbf{  Deviance:          } &    18.086   \\\\\n\\textbf{Time:}            &     22:16:08     & \\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}"
          },
          "metadata": {}
        }
      ],
      "execution_count": 15
    },
    {
      "cell_type": "markdown",
      "source": "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\nThe 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**.",
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": "## Predicting failure probability\nThe temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "%matplotlib inline\ndata_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\ndata_pred['Frequency'] = logmodel.predict(data_pred)\ndata_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\nplt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\nplt.grid(True)",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 640x480 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzZUlEQVR4nO3deXBUVf7+8aezE0KCBMiiAcIii+BCGDQoiwtBQNRyRlBkE3DMgCLGDWSURQUdlUFLiYhsEhdUlK84EchYsggoEogDhgKUYFA7k5+gJBqTdNL39weTlqazdGfhQHi/qrqKe+65957+5EIe7mqzLMsSAACAIX6mBwAAAM5thBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABglM9hZPPmzRo2bJhiY2Nls9m0Zs2aGpfZtGmTEhISFBISovbt2+uVV16pzVgBAEAj5HMY+e2333TJJZfopZde8qp/Tk6OhgwZor59+2r37t169NFHNWXKFK1evdrnwQIAgMbHVpcX5dlsNn3wwQe6+eabq+zzyCOP6MMPP9S+fftcbcnJyfrqq6+0ffv22m4aAAA0EgENvYHt27crKSnJrW3QoEFasmSJHA6HAgMDPZYpKSlRSUmJa9rpdOrYsWOKjIyUzWZr6CEDAIB6YFmWCgsLFRsbKz+/qk/GNHgYycvLU1RUlFtbVFSUysrK9NNPPykmJsZjmXnz5mn27NkNPTQAAHAaHDlyRBdccEGV8xs8jEjyOJpRcWaoqqMc06dPV0pKimv6+PHjatOmjXJyctSsWbN6GZNlWSosKtGmzZvVv18/BQSellKc1cocZdTLS9TKe9TKe9TKe9TKexW1GnTtAAUFBdXrugsLCxUfH1/j7+4G/wlFR0crLy/PrS0/P18BAQGKjIysdJng4GAFBwd7tLdo0ULh4eH1NrYIh0PnNQvVBTGtKz1dBHcO6uU1auU9auU9auU9auW9ilq1bNmy3mtVsb6aLrFo8OeMJCYmKiMjw61tw4YN6tWrFzsIAADwPYz8+uuvysrKUlZWlqQTt+5mZWUpNzdX0olTLGPGjHH1T05O1nfffaeUlBTt27dPS5cu1ZIlS/Tggw/WzzcAAABnNZ9P0+zcuVNXX321a7ri2o6xY8dq+fLlstvtrmAiSfHx8UpPT9f999+vl19+WbGxsXrxxRf15z//uR6GDwAAznY+h5EBAwaoukeTLF++3KOtf//+2rVrl6+bAgA0EuXl5XI4HKdtew6HQwEBASouLlZ5eflp2+7ZqC61CgwMlL+/f53HwCXGAIAGY1mW8vLy9Msvv5z27UZHR+vIkSM8n6oGda1V8+bNFR0dXac6E0YAAA2mIoi0bt1aoaGhpy0YOJ1O/frrrwoLC6v2YVuofa0sy1JRUZHy8/MlqdLnhnmLMAIAaBDl5eWuIFLVoxwaitPpVGlpqUJCQggjNahLrZo0aSLpxCM7WrduXetTNvyEAAANouIakdDQUMMjQUOq+PnW5ZogwggAoEFxzUbjVh8/X8IIAAAwijACAACMIowAAHCKcePGyWazeXy++eYb00NrlLibBgCASlx//fVatmyZW1urVq3cpktLS+v9TbfnIo6MAABQieDgYEVHR7t9rr32Wt1zzz1KSUlRy5YtNXDgQElSdna2hgwZorCwMEVFRWn06NH66aefXOv67bffNGbMGIWFhSkmJkbPP/+8BgwYoKlTp7r62Gw2rVmzxm0MzZs3d3uy+Q8//KARI0bovPPOU2RkpG666SYdPnzYNX/cuHG6+eab9dxzzykmJkaRkZGaPHmy250uJSUlevjhhxUXF6fg4GB17txZK1eulGVZ6tixo5577jm3Mezdu1d+fn769ttv617UKhBGAACnjWVZKiotOy2f30vLXX+u7jUmvlqxYoUCAgK0detWLVq0SHa7Xf3799ell16qnTt3at26dfrvf/+r4cOHu5Z56KGH9Omnn+qDDz7Qhg0btHHjRmVmZvq03aKiIl199dUKCwvT5s2b9dlnnyksLEzXX3+9SktLXf0+/fRTffvtt/r000+1YsUKLV++3C3QjBkzRm+//bZefPFF7du3TwsXLlTTpk1ls9k0fvx4j6NBS5cuVd++fdWhQ4faFcwLnKYBAJw2vzvK1e3x9ad9u9lzBik0yLdfeR999JHCwsJc04MHD5YkdezYUf/4xz9c7Y8//rh69uypuXPnutqWLl2quLg4HThwQLGxsVqyZIlef/1115GUFStW6IILLvBpPG+//bb8/Pz02muvuW6nXbZsmZo3b66NGzcqKSlJknTeeefppZdekr+/v7p06aKhQ4fqk08+0V133aUDBw7onXfeUUZGhq677jpJUrt27VRQUCBJuvPOO/X4449rx44d6t27txwOh9LS0vTss8/6NFZfEUYAAKjE1VdfrdTUVNd006ZNdfvtt6tXr15u/TIzM/Xpp5+6BZcK3377rX7//XeVlpYqMTHR1d6iRQt17tzZp/FkZmbqm2++UbNmzdzai4uL3U6hXHTRRW5PQo2JidGePXskSVlZWfL391f//v0r3UZMTIyGDh2qpUuXqnfv3vroo49UXFysW2+91aex+oowAgA4bZoE+it7zqAG347T6VRhQaGahTeTn5+fmgT6/pjypk2bqmPHjpW2n7qtYcOG6ZlnnvHoGxMTo4MHD3q1PZvN5nE66eRrPZxOpxISEvTGG294LHvyhbWBgYEe63U6nZL+eHx7dSZOnKjRo0frn//8p5YtW6YRI0Y0+FN0CSMAgNPGZrP5fLqkNpxOp8qC/BUaFNDg76bp2bOnVq9erXbt2ikgwPO7dezYUYGBgfr888/Vpk0bSdLPP/+sAwcOuB2haNWqlex2u2v64MGDKioqctvOqlWr1Lp1a4WHh9dqrD169JDT6dSmTZtcp2lONWTIEDVt2lSpqan6+OOPtXnz5lptyxdcwAoAQB1MnjxZx44d0+23364dO3bo0KFD2rBhg8aPH6/y8nKFhYVpwoQJeuihh/TJJ59o7969GjdunEdIuuaaa/TSSy9p165d2rlzp5KTk92Octxxxx1q2bKlbrrpJm3ZskU5OTnatGmT7rvvPn3//fdejbVdu3YaO3asxo8frzVr1ignJ0cbN27UBx984Orj7++vcePGafr06erYsaPb6aWGQhgBAKAOYmNjtXXrVpWXl2vQoEHq3r277rvvPkVERLgCx7PPPqt+/frpxhtv1HXXXaerrrpKCQkJbut5/vnnFRcXp379+mnkyJF68MEH3U6PhIaGavPmzWrTpo1uueUWde3aVePHj9fvv//u05GS1NRU/eUvf9GkSZPUpUsX3X333W5HYCRpwoQJKi0t1fjx4+tQGe9xmgYAgFOcfCvsyTZu3Fhpe6dOnfT+++9Xub6wsDCtXLlSK1eudLX961//cusTGxur9evd7zT65Zdf3Kajo6O1YsUKn8a9YMECt+mQkBDNnz9f8+fPl3TilFbF3TQV7Ha7AgICNGbMmCq3VZ8IIwAAQNKJB6IdOXJEjz32mIYPH66oqKjTsl1O0wAAAEnSW2+9pc6dO+v48eNuz1JpaBwZAQDAgKpO+Zg0btw4jRs37rRvlyMjAADAKMIIAKBB1ed7YXDmqY+fL2EEANAgKp6Rcepto2hcKn6+pz751RdcMwIAaBD+/v5q3ry58vPzJZ14TkbFC94amtPpVGlpqYqLixv8Caxnu9rWyrIsFRUVKT8/X82bN3d7H46vCCMAgAYTHR0tSa5AcrpYlqXff/9dTZo0OW0B6GxV11o1b97c9XOuLcIIAKDB2Gw2xcTEqHXr1m4vfWtoDodDmzdvVr9+/ep0+uBcUJdaBQYG1umISAXCCACgwfn7+9fLLy1ftldWVqaQkBDCSA3OhFpxIg0AABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRtQojCxcuVHx8vEJCQpSQkKAtW7ZU2/+NN97QJZdcotDQUMXExOjOO+/U0aNHazVgAADQuPgcRlatWqWpU6dqxowZ2r17t/r27avBgwcrNze30v6fffaZxowZowkTJujrr7/Wu+++qy+//FITJ06s8+ABAMDZz+cwMn/+fE2YMEETJ05U165dtWDBAsXFxSk1NbXS/p9//rnatWunKVOmKD4+XldddZXuvvtu7dy5s86DBwAAZ78AXzqXlpYqMzNT06ZNc2tPSkrStm3bKl2mT58+mjFjhtLT0zV48GDl5+frvffe09ChQ6vcTklJiUpKSlzTBQUFkiSHwyGHw+HLkKtVsa76XGdjRr28R628R628R628R62815C18nadNsuyLG9X+uOPP+r888/X1q1b1adPH1f73LlztWLFCu3fv7/S5d577z3deeedKi4uVllZmW688Ua99957CgwMrLT/rFmzNHv2bI/2N998U6Ghod4OFwAAGFRUVKSRI0fq+PHjCg8Pr7KfT0dGKthsNrdpy7I82ipkZ2drypQpevzxxzVo0CDZ7XY99NBDSk5O1pIlSypdZvr06UpJSXFNFxQUKC4uTklJSdV+GV85HA5lZGRo4MCBVQYj/IF6eY9aeY9aeY9aeY9aea8ha1VxZqMmPoWRli1byt/fX3l5eW7t+fn5ioqKqnSZefPm6corr9RDDz0kSbr44ovVtGlT9e3bV08++aRiYmI8lgkODlZwcLBHe2BgYIPsVA213saKenmPWnmPWnmPWnmPWnmvIWrl7fp8uoA1KChICQkJysjIcGvPyMhwO21zsqKiIvn5uW/G399f0okjKgAA4Nzm8900KSkpeu2117R06VLt27dP999/v3Jzc5WcnCzpxCmWMWPGuPoPGzZM77//vlJTU3Xo0CFt3bpVU6ZMUe/evRUbG1t/3wQAAJyVfL5mZMSIETp69KjmzJkju92u7t27Kz09XW3btpUk2e12t2eOjBs3ToWFhXrppZf0wAMPqHnz5rrmmmv0zDPP1N+3AAAAZ61aXcA6adIkTZo0qdJ5y5cv92i79957de+999ZmUwAAoJHj3TQAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIyqVRhZuHCh4uPjFRISooSEBG3ZsqXa/iUlJZoxY4batm2r4OBgdejQQUuXLq3VgAEAQOMS4OsCq1at0tSpU7Vw4UJdeeWVWrRokQYPHqzs7Gy1adOm0mWGDx+u//73v1qyZIk6duyo/Px8lZWV1XnwAADg7OdzGJk/f74mTJigiRMnSpIWLFig9evXKzU1VfPmzfPov27dOm3atEmHDh1SixYtJEnt2rWr26gBAECj4VMYKS0tVWZmpqZNm+bWnpSUpG3btlW6zIcffqhevXrpH//4h1auXKmmTZvqxhtv1BNPPKEmTZpUukxJSYlKSkpc0wUFBZIkh8Mhh8Phy5CrVbGu+lxnY0a9vEetvEetvEetvEetvNeQtfJ2nT6FkZ9++knl5eWKiopya4+KilJeXl6lyxw6dEifffaZQkJC9MEHH+inn37SpEmTdOzYsSqvG5k3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/Xw+TSNJNpvNbdqyLI+2Ck6nUzabTW+88YYiIiIknTjV85e//EUvv/xypUdHpk+frpSUFNd0QUGB4uLilJSUpPDw8NoMuVIOh0MZGRkaOHCgAgMD6229jRX18h618h618h618h618l5D1qrizEZNfAojLVu2lL+/v8dRkPz8fI+jJRViYmJ0/vnnu4KIJHXt2lWWZen7779Xp06dPJYJDg5WcHCwR3tgYGCD7FQNtd7Ginp5j1p5j1p5j1p5j1p5ryFq5e36fLq1NygoSAkJCR6HcjIyMtSnT59Kl7nyyiv1448/6tdff3W1HThwQH5+frrgggt82TwAAGiEfH7OSEpKil577TUtXbpU+/bt0/3336/c3FwlJydLOnGKZcyYMa7+I0eOVGRkpO68805lZ2dr8+bNeuihhzR+/PgqL2AFAADnDp+vGRkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7Xbm5ua7+YWFhysjI0L333qtevXopMjJSw4cP15NPPll/3wIAAJy1anUB66RJkzRp0qRK5y1fvtyjrUuXLlzRDAAAKsW7aQAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGFWrMLJw4ULFx8crJCRECQkJ2rJli1fLbd26VQEBAbr00ktrs1kAANAI+RxGVq1apalTp2rGjBnavXu3+vbtq8GDBys3N7fa5Y4fP64xY8bo2muvrfVgAQBA4+NzGJk/f74mTJigiRMnqmvXrlqwYIHi4uKUmppa7XJ33323Ro4cqcTExFoPFgAAND4BvnQuLS1VZmampk2b5taelJSkbdu2VbncsmXL9O233yotLU1PPvlkjdspKSlRSUmJa7qgoECS5HA45HA4fBlytSrWVZ/rbMyol/eolfeolfeolfeolfcaslbertOnMPLTTz+pvLxcUVFRbu1RUVHKy8urdJmDBw9q2rRp2rJliwICvNvcvHnzNHv2bI/2DRs2KDQ01JcheyUjI6Pe19mYUS/vUSvvUSvvUSvvUSvvNUStioqKvOrnUxipYLPZ3KYty/Jok6Ty8nKNHDlSs2fP1oUXXuj1+qdPn66UlBTXdEFBgeLi4pSUlKTw8PDaDLlSDodDGRkZGjhwoAIDA+ttvY0V9fIetfIetfIetfIetfJeQ9aq4sxGTXwKIy1btpS/v7/HUZD8/HyPoyWSVFhYqJ07d2r37t265557JElOp1OWZSkgIEAbNmzQNddc47FccHCwgoODPdoDAwMbZKdqqPU2VtTLe9TKe9TKe9TKe9TKew1RK2/X59MFrEFBQUpISPA4lJORkaE+ffp49A8PD9eePXuUlZXl+iQnJ6tz587KysrS5Zdf7svmAQBAI+TzaZqUlBSNHj1avXr1UmJiol599VXl5uYqOTlZ0olTLD/88INef/11+fn5qXv37m7Lt27dWiEhIR7tAADg3ORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7vcZnjgAAAFSo1QWskyZN0qRJkyqdt3z58mqXnTVrlmbNmlWbzQIAgEaId9MAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjavXWXgCnX7nT0o6cY8ovLFbrZiHqHd9C/n4208PCOY79EvWBMAKcBdbttWv22mzZjxe72mIiQjRzWDdd3z3G4MhwLmO/RH3hNA1whlu3166/pe1y+wdfkvKOF+tvabu0bq/d0MhwLmO/RH0ijABnsHKnpdlrs2VVMq+ibfbabJU7K+sBNAz2S9Q3wghwBtuRc8zjf54nsyTZjxdrR86x0zconPPYL1HfCCPAGSy/sOp/8GvTD6gP7Jeob4QR4AzWullIvfYD6gP7JeobYQQ4g/WOb6GYiBBVdaOkTSfuXugd3+J0DgvnOPZL1DfCCHAG8/ezaeawbpLk8Q9/xfTMYd14rgNOK/ZL1DfCCHCGu757jFJH9VR0hPsh7+iIEKWO6snzHGAE+yXqEw89A84C13eP0cBu0TzpEmcU9kvUF8IIcJbw97MpsUOk6WEAbtgvUR84TQMAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKNqFUYWLlyo+Ph4hYSEKCEhQVu2bKmy7/vvv6+BAweqVatWCg8PV2JiotavX1/rAQMAgMbF5zCyatUqTZ06VTNmzNDu3bvVt29fDR48WLm5uZX237x5swYOHKj09HRlZmbq6quv1rBhw7R79+46Dx4AAJz9fA4j8+fP14QJEzRx4kR17dpVCxYsUFxcnFJTUyvtv2DBAj388MP605/+pE6dOmnu3Lnq1KmT1q5dW+fBAwCAs1+AL51LS0uVmZmpadOmubUnJSVp27ZtXq3D6XSqsLBQLVq0qLJPSUmJSkpKXNMFBQWSJIfDIYfD4cuQq1WxrvpcZ2NGvbxHrbxHrbxHrbxHrbzXkLXydp0+hZGffvpJ5eXlioqKcmuPiopSXl6eV+t4/vnn9dtvv2n48OFV9pk3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/XwKIxVsNpvbtGVZHm2VeeuttzRr1iz93//9n1q3bl1lv+nTpyslJcU1XVBQoLi4OCUlJSk8PLw2Q66Uw+FQRkaGBg4cqMDAwHpbb2NFvbxHrbxHrbxHrbxHrbzXkLWqOLNRE5/CSMuWLeXv7+9xFCQ/P9/jaMmpVq1apQkTJujdd9/VddddV23f4OBgBQcHe7QHBgY2yE7VUOttrKiX96iV96iV96iV96iV9xqiVt6uz6cLWIOCgpSQkOBxKCcjI0N9+vSpcrm33npL48aN05tvvqmhQ4f6skkAANDI+XyaJiUlRaNHj1avXr2UmJioV199Vbm5uUpOTpZ04hTLDz/8oNdff13SiSAyZswYvfDCC7riiitcR1WaNGmiiIiIevwqAADgbORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa73e2ZI4sWLVJZWZkmT56syZMnu9rHjh2r5cuX1/0bAACAs1qtLmCdNGmSJk2aVOm8UwPGxo0ba7MJAABwjqhVGAFw7ih3WtqRc0z5hcVq3SxEveNbyN/P5vV8E87EMdVVaZlTadsPK1LSyu2HNapPBwUF8HoxNA6EEQBVWrfXrtlrs2U/Xuxqi4kI0cxh3XR995ga55twJo6prualZ2vxlhwF+ln6R2/pmfX79eTHB3RX33hNH9LN9PCAOiNWA6jUur12/S1tl9svdUnKO16sv6Xt0rz07Grnr9trP53DlVTzmE2Mqa7mpWdr0eYcOS33dqclLdqco3np2WYGBtQjwggAD+VOS7PXZsuqZJ71v8/iLTlVzpek2WuzVX7qb9AGVNOYTYyprkrLnFq8JafaPou35Ki0zHmaRgQ0DMIIAA87co55HF04VXW/0y1J9uPF2pFzrH4HVo2axmxiTHW1cvvhaussnfg5rNx++LSMB2gohBEAHvILqw8ip3s99bmt0zmmuvrumHfv9fC2H3CmIowA8NC6WcgZtZ763NbpHFNdtW3h3YtBve0HnKkIIwA89I5voZiIEFV3M6yfTVXOt+nEHSy941s0wOgqV9OYTYyprkYntlNNdyT72U70A85mhBEAHvz9bJo57MQto6f+LrT973NX3/gq50vSzGHdTuuzPWoas4kx1VVQgJ+rzlW5q288zxvBWY89GEClru8eo9RRPRUd4X5aIzoiRKmjemr6kG7VzjfxTI+axnw2Pmdk+pBuurtfvMcREj+bdHc/njOCxoGHngGo0vXdYzSwW3SVTzOtaf6ZOOaz0fQh3fRAUhelbftW+jlbjwzqzBNY0agQRgBUy9/PpsQOkbWeb8KZOKa6Cgrw0+jEdkpPz9boxHYKJIigEWFvBgAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUQGmBwAAZ5Nyp6UdOceUX1is1s1C1Du+hfz9bJKk30vLNTc9W4ePFqldZKgeHdJNTYL8vVq2unmSVFrmVNr2w4qUtHL7YY3q00FBAd79f7Kmddc0v7brLi1zauX2w/ruWJHatgjV6MR2Z/yYYUatwsjChQv17LPPym6366KLLtKCBQvUt2/fKvtv2rRJKSkp+vrrrxUbG6uHH35YycnJtR40AJiwbq9ds9dmy3682NUWExGimcO6afWu75WRne9q33JQWvl5rgZ2a63FY/5U7bKSqpx3ffcYzUvP1uItOQr0s/SP3tIz6/fryY8P6K6+8Zo+pFutx3x995ga59d23btzf9biLTlyWn/0fyp93xk95pqWRcPxOYysWrVKU6dO1cKFC3XllVdq0aJFGjx4sLKzs9WmTRuP/jk5ORoyZIjuuusupaWlaevWrZo0aZJatWqlP//5z/XyJQCgoa3ba9ff0nbJOqU973ixktN2VblcRna+bnxpi/Z8X+DTsnnHi/W3tF26rltrt5BTwWlJizbnSFKVv9yrG/Pf0nbpr/3i9ermnCrnp47qWeUv6NrU40wec03LomH5fM3I/PnzNWHCBE2cOFFdu3bVggULFBcXp9TU1Er7v/LKK2rTpo0WLFigrl27auLEiRo/fryee+65Og8eAE6Hcqel2WuzPX6JSaq07VT/qSSI1LSs9b9PZUHkZIu35Ki0zOnRXtOYrf8tW924Zq/NVrnTs0dd63GmjrmqZdHwfDoyUlpaqszMTE2bNs2tPSkpSdu2bat0me3btyspKcmtbdCgQVqyZIkcDocCAwM9likpKVFJSYlr+vjx45KkY8eOyeFw+DLkajkcDhUVFeno0aOVjgPuqJf3qJX3zoZaZX73s/7f0aPGL7ILcFoqKnIqwOGncucf1zgszvhKt/V2PzLt7Zir+x/p/zv6mz7J+lYJbc+r1bqr09Bjvjg2zG2/8mbdVX3fxq4h/w4WFhZKkiyrhpBn+eCHH36wJFlbt251a3/qqaesCy+8sNJlOnXqZD311FNubVu3brUkWT/++GOly8ycObMiBPPhw4cPHz58zvLPkSNHqs0XtQq2Npv7VceWZXm01dS/svYK06dPV0pKimva6XTq2LFjioyMrHY7viooKFBcXJyOHDmi8PDweltvY0W9vEetvEetvEetvEetvNeQtbIsS4WFhYqNja22n09hpGXLlvL391deXp5be35+vqKioipdJjo6utL+AQEBioyMrHSZ4OBgBQcHu7U1b97cl6H6JDw8nJ3VB9TLe9TKe9TKe9TKe9TKew1Vq4iIiBr7+HQBa1BQkBISEpSRkeHWnpGRoT59+lS6TGJiokf/DRs2qFevXmfs+WEAAHD6+Hw3TUpKil577TUtXbpU+/bt0/3336/c3FzXc0OmT5+uMWPGuPonJyfru+++U0pKivbt26elS5dqyZIlevDBB+vvWwAAgLOWz9eMjBgxQkePHtWcOXNkt9vVvXt3paenq23btpIku92u3NxcV//4+Hilp6fr/vvv18svv6zY2Fi9+OKLZ8QzRoKDgzVz5kyPU0KoHPXyHrXyHrXyHrXyHrXy3plQK5tl1XS/DQAAQMPhRXkAAMAowggAADCKMAIAAIwijAAAAKPOiTCSmpqqiy++2PVAl8TERH388ceu+ZZladasWYqNjVWTJk00YMAAff311wZHfGaYN2+ebDabpk6d6mqjVn+YNWuWbDab2yc6Oto1n1q5++GHHzRq1ChFRkYqNDRUl156qTIzM13zqdcJ7dq189ivbDabJk+eLIk6naysrEx///vfFR8fryZNmqh9+/aaM2eOnM4/XsJHvf5QWFioqVOnqm3btmrSpIn69OmjL7/80jXfaK1qeB1No/Dhhx9a//rXv6z9+/db+/fvtx599FErMDDQ2rt3r2VZlvX0009bzZo1s1avXm3t2bPHGjFihBUTE2MVFBQYHrk5O3bssNq1a2ddfPHF1n333edqp1Z/mDlzpnXRRRdZdrvd9cnPz3fNp1Z/OHbsmNW2bVtr3Lhx1hdffGHl5ORY//73v61vvvnG1Yd6nZCfn++2T2VkZFiSrE8//dSyLOp0sieffNKKjIy0PvroIysnJ8d69913rbCwMGvBggWuPtTrD8OHD7e6detmbdq0yTp48KA1c+ZMKzw83Pr+++8tyzJbq3MijFTmvPPOs1577TXL6XRa0dHR1tNPP+2aV1xcbEVERFivvPKKwRGaU1hYaHXq1MnKyMiw+vfv7woj1MrdzJkzrUsuuaTSedTK3SOPPGJdddVVVc6nXlW77777rA4dOlhOp5M6nWLo0KHW+PHj3dpuueUWa9SoUZZlsV+drKioyPL397c++ugjt/ZLLrnEmjFjhvFanROnaU5WXl6ut99+W7/99psSExOVk5OjvLw8JSUlufoEBwerf//+2rZtm8GRmjN58mQNHTpU1113nVs7tfJ08OBBxcbGKj4+XrfddpsOHTokiVqd6sMPP1SvXr106623qnXr1rrsssu0ePFi13zqVbnS0lKlpaVp/Pjxstls1OkUV111lT755BMdOHBAkvTVV1/ps88+05AhQySxX52srKxM5eXlCgkJcWtv0qSJPvvsM+O1OmfCyJ49exQWFqbg4GAlJyfrgw8+ULdu3Vwv8Tv1RX9RUVEeL/g7F7z99tvatWuX5s2b5zGPWrm7/PLL9frrr2v9+vVavHix8vLy1KdPHx09epRaneLQoUNKTU1Vp06dtH79eiUnJ2vKlCl6/fXXJbFvVWXNmjX65ZdfNG7cOEnU6VSPPPKIbr/9dnXp0kWBgYG67LLLNHXqVN1+++2SqNfJmjVrpsTERD3xxBP68ccfVV5errS0NH3xxRey2+3Ga+Xz4+DPVp07d1ZWVpZ++eUXrV69WmPHjtWmTZtc8202m1t/y7I82hq7I0eO6L777tOGDRs80vPJqNUJgwcPdv25R48eSkxMVIcOHbRixQpdccUVkqhVBafTqV69emnu3LmSpMsuu0xff/21UlNT3d5lRb3cLVmyRIMHD/Z4/Tp1OmHVqlVKS0vTm2++qYsuukhZWVmaOnWqYmNjNXbsWFc/6nXCypUrNX78eJ1//vny9/dXz549NXLkSO3atcvVx1StzpkjI0FBQerYsaN69eqlefPm6ZJLLtELL7zguvvh1OSXn5/vkRAbu8zMTOXn5yshIUEBAQEKCAjQpk2b9OKLLyogIMBVD2pVuaZNm6pHjx46ePAg+9UpYmJi1K1bN7e2rl27ut5jRb08fffdd/r3v/+tiRMnutqok7uHHnpI06ZN02233aYePXpo9OjRuv/++11HdqmXuw4dOmjTpk369ddfdeTIEe3YsUMOh0Px8fHGa3XOhJFTWZalkpIS1w8hIyPDNa+0tFSbNm1Snz59DI7w9Lv22mu1Z88eZWVluT69evXSHXfcoaysLLVv355aVaOkpET79u1TTEwM+9UprrzySu3fv9+t7cCBA64XbFIvT8uWLVPr1q01dOhQVxt1cldUVCQ/P/dfY/7+/q5be6lX5Zo2baqYmBj9/PPPWr9+vW666SbztWrwS2TPANOnT7c2b95s5eTkWP/5z3+sRx991PLz87M2bNhgWdaJ25kiIiKs999/39qzZ491++23n7O3fp3q5LtpLItaneyBBx6wNm7caB06dMj6/PPPrRtuuMFq1qyZdfjwYcuyqNXJduzYYQUEBFhPPfWUdfDgQeuNN96wQkNDrbS0NFcf6vWH8vJyq02bNtYjjzziMY86/WHs2LHW+eef77q19/3337datmxpPfzww64+1OsP69atsz7++GPr0KFD1oYNG6xLLrnE6t27t1VaWmpZltlanRNhZPz48Vbbtm2toKAgq1WrVta1117rCiKWdeL2r5kzZ1rR0dFWcHCw1a9fP2vPnj0GR3zmODWMUKs/VNyDHxgYaMXGxlq33HKL9fXXX7vmUyt3a9eutbp3724FBwdbXbp0sV599VW3+dTrD+vXr7ckWfv37/eYR53+UFBQYN13331WmzZtrJCQEKt9+/bWjBkzrJKSElcf6vWHVatWWe3bt7eCgoKs6Ohoa/LkydYvv/zimm+yVjbLsqyGP/4CAABQuXP2mhEAAHBmIIwAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAjRCNput2k/FK+kbkwEDBmjq1KmmhwGgFgJMDwBA/bPb7a4/r1q1So8//rjbi+qaNGliYli14nA4FBgY2Gi3B4AjI0CjFB0d7fpERETIZrO5tW3evFkJCQkKCQlR+/btNXv2bJWVlbmWt9lsWrRokW644QaFhoaqa9eu2r59u7755hsNGDBATZs2VWJior799lvXMrNmzdKll16qRYsWKS4uTqGhobr11lv1yy+/uI1t2bJl6tq1q0JCQtSlSxctXLjQNe/w4cOy2Wx65513NGDAAIWEhCgtLU1Hjx7V7bffrgsuuEChoaHq0aOH3nrrLddy48aN06ZNm/TCCy+4jv4cPnxYy5cvV/Pmzd22v2bNGtlsNo9xL126VO3bt1dwcLAsy9Lx48f117/+Va1bt1Z4eLiuueYaffXVV/X0EwJwMsIIcI5Zv369Ro0apSlTpig7O1uLFi3S8uXL9dRTT7n1e+KJJzRmzBhlZWWpS5cuGjlypO6++25Nnz5dO3fulCTdc889bst88803euedd7R27VqtW7dOWVlZmjx5smv+4sWLNWPGDD311FPat2+f5s6dq8cee0wrVqxwW88jjzyiKVOmaN++fRo0aJCKi4uVkJCgjz76SHv37tVf//pXjR49Wl988YUk6YUXXlBiYqLuuusu2e122e12xcXFeV2TinGvXr1aWVlZkqShQ4cqLy9P6enpyszMVM+ePXXttdfq2LFjXq8XgJdOy+v4ABizbNkyKyIiwjXdt29fa+7cuW59Vq5cacXExLimJVl///vfXdPbt2+3JFlLlixxtb311ltWSEiIa3rmzJmWv7+/deTIEVfbxx9/bPn5+Vl2u92yLMuKi4uz3nzzTbdtP/HEE1ZiYqJlWZaVk5NjSbIWLFhQ4/caMmSI9cADD7imT33DdGXf3bIs64MPPrBO/qdv5syZVmBgoJWfn+9q++STT6zw8HCruLjYbdkOHTpYixYtqnFsAHzDNSPAOSYzM1Nffvml25GQ8vJyFRcXq6ioSKGhoZKkiy++2DU/KipKktSjRw+3tuLiYhUUFCg8PFyS1KZNG11wwQWuPomJiXI6ndq/f7/8/f115MgRTZgwQXfddZerT1lZmSIiItzG2KtXL7fp8vJyPf3001q1apV++OEHlZSUqKSkRE2bNq1rOSRJbdu2VatWrVzTmZmZ+vXXXxUZGenW7/fff3c7NQWgfhBGgHOM0+nU7Nmzdcstt3jMCwkJcf355Is4K66xqKzN6XRWua2KPjabzdVv8eLFuvzyy936+fv7u02fGjKef/55/fOf/9SCBQvUo0cPNW3aVFOnTlVpaWnVX1SSn5+fLMtya3M4HB79Tt2e0+lUTEyMNm7c6NH31GtQANQdYQQ4x/Ts2VP79+9Xx44d633dubm5+vHHHxUbGytJ2r59u/z8/HThhRcqKipK559/vg4dOqQ77rjDp/Vu2bJFN910k0aNGiXpRFg4ePCgunbt6uoTFBSk8vJyt+VatWqlwsJC/fbbb67AUXFNSHV69uypvLw8BQQEqF27dj6NFYDvCCPAOebxxx/XDTfcoLi4ON16663y8/PTf/7zH+3Zs0dPPvlkndYdEhKisWPH6rnnnlNBQYGmTJmi4cOHKzo6WtKJO1emTJmi8PBwDR48WCUlJdq5c6d+/vlnpaSkVLnejh07avXq1dq2bZvOO+88zZ8/X3l5eW5hpF27dvriiy90+PBhhYWFqUWLFrr88ssVGhqqRx99VPfee6927Nih5cuX1/g9rrvuOiUmJurmm2/WM888o86dO+vHH39Uenq6br75Zo/TSADqhrtpgHPMoEGD9NFHHykjI0N/+tOfdMUVV2j+/Plq27ZtndfdsWNH3XLLLRoyZIiSkpLUvXt3t1t3J06cqNdee03Lly9Xjx491L9/fy1fvlzx8fHVrvexxx5Tz549NWjQIA0YMEDR0dG6+eab3fo8+OCD8vf3V7du3dSqVSvl5uaqRYsWSktLU3p6uut24FmzZtX4PWw2m9LT09WvXz+NHz9eF154oW677TYdPnzYdf0MgPpjs049oQoAtTBr1iytWbPGq9MgAHAyjowAAACjCCMAAMAoTtMAAACjODICAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjPr/yAWkb36Wd7EAAAAASUVORK5CYII="
          },
          "metadata": {}
        }
      ],
      "execution_count": 16
    },
    {
      "cell_type": "markdown",
      "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.**",
      "metadata": {
        "hideCode": false,
        "hidePrompt": false,
        "scrolled": true
      }
    },
    {
      "cell_type": "markdown",
      "source": "## Computing and plotting uncertainty",
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot.",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "#sns.set(color_codes=True)\nplt.xlim(30,90)\nplt.ylim(0,1)\n#sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\nplt.show()",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 640x480 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGiCAYAAADEJZ3cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAd4UlEQVR4nO3df3DX9X3A8VdIIMEfSSfUkFiIwdoWRakkJyWKXWubHTJvXHcTdRUoerfc2DCkOorsijK3dL1b13YdWCvUOWnLOdGxjgrpbkUstlZKHENOWaEGbdIcOBLULZTw2R+e31uaUPMNsbxLHo+7zx3f9/f9+X7feV/ueN7n+yMFWZZlAQCQsFGnewEAAG9HsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJyztYnnzyybj++uujsrIyCgoK4vHHH3/bc7Zt2xY1NTVRUlISkydPjvvuu28oawUARqi8g+X111+PadOmxVe+8pVBzT9w4EBcd911MWvWrNi1a1fcddddsWTJknj00UfzXiwAMDIVnMofPywoKIjHHnss5s6de9I5y5Yti02bNsXevXtzYw0NDfHcc8/F008/PdSnBgBGkKJ3+gmefvrpqK+v7zP2O7/zO7F27dr4xS9+EaNHj+53Tk9PT/T09ORunzhxIl599dUYN25cFBQUvNNLBgCGQZZlcfTo0aisrIxRo07tbbPveLB0dHREeXl5n7Hy8vI4fvx4HDp0KCoqKvqd09zcHPfcc887vTQA4Nfg4MGD8Z73vOeUHuMdD5aI6HdV5K1XoU52tWT58uXR1NSUu93V1RWTJk2KgwcPRmlp6Tu3UABg2HR3d8fEiRPj3HPPPeXHeseDZcKECdHR0dFnrLOzM4qKimLcuHEDnlNcXBzFxcX9xktLSwULAPyGGY63c7zj38Myc+bMaGlp6TO2devWqK2tHfD9KwAAvyzvYHnttdeitbU1WltbI+LNjy23trZGW1tbRLz5cs78+fNz8xsaGuKll16Kpqam2Lt3b6xbty7Wrl0bd9xxx/D8BADAGS/vl4SeffbZ+MhHPpK7/dZ7TRYsWBAPPvhgtLe35+IlIqK6ujo2b94cS5cujb//+7+PysrK+PKXvxy///u/PwzLBwBGglP6HpZfl+7u7igrK4uuri7vYQGA3xDD+f+3vyUEACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyhhQsq1evjurq6igpKYmamprYvn37r5y/fv36mDZtWpx11llRUVERn/rUp+Lw4cNDWjAAMPLkHSwbNmyIxsbGWLFiRezatStmzZoVs2fPjra2tgHnP/XUUzF//vy49dZbY8+ePfHII4/Ej370o7jttttOefEAwMiQd7B84QtfiFtvvTVuu+22mDJlSnzxi1+MiRMnxpo1awac/4Mf/CAuvPDCWLJkSVRXV8fVV18df/RHfxTPPvvsKS8eABgZ8gqWY8eOxc6dO6O+vr7PeH19fezYsWPAc+rq6uLll1+OzZs3R5Zl8fOf/zz+6Z/+KebMmXPS5+np6Ynu7u4+BwAwcuUVLIcOHYre3t4oLy/vM15eXh4dHR0DnlNXVxfr16+PefPmxZgxY2LChAnxrne9K/7u7/7upM/T3NwcZWVluWPixIn5LBMAOMMM6U23BQUFfW5nWdZv7C3PP/98LFmyJD772c/Gzp0744knnogDBw5EQ0PDSR9/+fLl0dXVlTsOHjw4lGUCAGeIonwmjx8/PgoLC/tdTens7Ox31eUtzc3NcdVVV8Wdd94ZERGXX355nH322TFr1qy49957o6Kiot85xcXFUVxcnM/SAIAzWF5XWMaMGRM1NTXR0tLSZ7ylpSXq6uoGPOeNN96IUaP6Pk1hYWFEvHllBgDg7eT9klBTU1M88MADsW7duti7d28sXbo02traci/xLF++PObPn5+bf/3118fGjRtjzZo1sX///vj+978fS5YsiSuvvDIqKyuH7ycBAM5Yeb0kFBExb968OHz4cKxatSra29tj6tSpsXnz5qiqqoqIiPb29j7fybJw4cI4evRofOUrX4lPf/rT8a53vSs++tGPxl//9V8P308BAJzRCrLfgNdluru7o6ysLLq6uqK0tPR0LwcAGITh/P/b3xICAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5ggUASJ5gAQCSJ1gAgOQJFgAgeYIFAEieYAEAkidYAIDkCRYAIHmCBQBInmABAJInWACA5AkWACB5QwqW1atXR3V1dZSUlERNTU1s3779V87v6emJFStWRFVVVRQXF8dFF10U69atG9KCAYCRpyjfEzZs2BCNjY2xevXquOqqq+KrX/1qzJ49O55//vmYNGnSgOfccMMN8fOf/zzWrl0b733ve6OzszOOHz9+yosHAEaGgizLsnxOmDFjRkyfPj3WrFmTG5syZUrMnTs3mpub+81/4okn4sYbb4z9+/fHeeedN6RFdnd3R1lZWXR1dUVpaemQHgMA+PUazv+/83pJ6NixY7Fz586or6/vM15fXx87duwY8JxNmzZFbW1tfP7zn48LLrgg3ve+98Udd9wR//M//3PS5+np6Ynu7u4+BwAwcuX1ktChQ4eit7c3ysvL+4yXl5dHR0fHgOfs378/nnrqqSgpKYnHHnssDh06FH/8x38cr7766knfx9Lc3Bz33HNPPksDAM5gQ3rTbUFBQZ/bWZb1G3vLiRMnoqCgINavXx9XXnllXHfddfGFL3whHnzwwZNeZVm+fHl0dXXljoMHDw5lmQDAGSKvKyzjx4+PwsLCfldTOjs7+111eUtFRUVccMEFUVZWlhubMmVKZFkWL7/8clx88cX9zikuLo7i4uJ8lgYAnMHyusIyZsyYqKmpiZaWlj7jLS0tUVdXN+A5V111VfzsZz+L1157LTf24osvxqhRo+I973nPEJYMAIw0eb8k1NTUFA888ECsW7cu9u7dG0uXLo22trZoaGiIiDdfzpk/f35u/s033xzjxo2LT33qU/H888/Hk08+GXfeeWcsWrQoxo4dO3w/CQBwxsr7e1jmzZsXhw8fjlWrVkV7e3tMnTo1Nm/eHFVVVRER0d7eHm1tbbn555xzTrS0tMSf/umfRm1tbYwbNy5uuOGGuPfee4fvpwAAzmh5fw/L6eB7WADgN89p+x4WAIDTQbAAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8oYULKtXr47q6uooKSmJmpqa2L59+6DO+/73vx9FRUXxwQ9+cChPCwCMUHkHy4YNG6KxsTFWrFgRu3btilmzZsXs2bOjra3tV57X1dUV8+fPj2uvvXbIiwUARqaCLMuyfE6YMWNGTJ8+PdasWZMbmzJlSsydOzeam5tPet6NN94YF198cRQWFsbjjz8era2tJ53b09MTPT09udvd3d0xceLE6OrqitLS0nyWCwCcJt3d3VFWVjYs/3/ndYXl2LFjsXPnzqivr+8zXl9fHzt27DjpeV//+tfjJz/5SaxcuXJQz9Pc3BxlZWW5Y+LEifksEwA4w+QVLIcOHYre3t4oLy/vM15eXh4dHR0DnrNv3774zGc+E+vXr4+ioqJBPc/y5cujq6srdxw8eDCfZQIAZ5jBFcQvKSgo6HM7y7J+YxERvb29cfPNN8c999wT73vf+wb9+MXFxVFcXDyUpQEAZ6C8gmX8+PFRWFjY72pKZ2dnv6suERFHjx6NZ599Nnbt2hV/8id/EhERJ06ciCzLoqioKLZu3Rof/ehHT2H5AMBIkNdLQmPGjImamppoaWnpM97S0hJ1dXX95peWlsbu3bujtbU1dzQ0NMT73//+aG1tjRkzZpza6gGAESHvl4Samprilltuidra2pg5c2bcf//90dbWFg0NDRHx5vtPXnnllXjooYdi1KhRMXXq1D7nn3/++VFSUtJvHADgZPIOlnnz5sXhw4dj1apV0d7eHlOnTo3NmzdHVVVVRES0t7e/7XeyAADkI+/vYTkdhvNz3ADAr8dp+x4WAIDTQbAAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8oYULKtXr47q6uooKSmJmpqa2L59+0nnbty4MT7+8Y/Hu9/97igtLY2ZM2fGli1bhrxgAGDkyTtYNmzYEI2NjbFixYrYtWtXzJo1K2bPnh1tbW0Dzn/yySfj4x//eGzevDl27twZH/nIR+L666+PXbt2nfLiAYCRoSDLsiyfE2bMmBHTp0+PNWvW5MamTJkSc+fOjebm5kE9xqWXXhrz5s2Lz372swPe39PTEz09Pbnb3d3dMXHixOjq6orS0tJ8lgsAnCbd3d1RVlY2LP9/53WF5dixY7Fz586or6/vM15fXx87duwY1GOcOHEijh49Guedd95J5zQ3N0dZWVnumDhxYj7LBADOMHkFy6FDh6K3tzfKy8v7jJeXl0dHR8egHuNv/uZv4vXXX48bbrjhpHOWL18eXV1duePgwYP5LBMAOMMUDeWkgoKCPrezLOs3NpBvfvObcffdd8c///M/x/nnn3/SecXFxVFcXDyUpQEAZ6C8gmX8+PFRWFjY72pKZ2dnv6suv2zDhg1x6623xiOPPBIf+9jH8l8pADBi5fWS0JgxY6KmpiZaWlr6jLe0tERdXd1Jz/vmN78ZCxcujG984xsxZ86coa0UABix8n5JqKmpKW655Zaora2NmTNnxv333x9tbW3R0NAQEW++/+SVV16Jhx56KCLejJX58+fHl770pfjQhz6UuzozduzYKCsrG8YfBQA4U+UdLPPmzYvDhw/HqlWror29PaZOnRqbN2+OqqqqiIhob2/v850sX/3qV+P48eOxePHiWLx4cW58wYIF8eCDD576TwAAnPHy/h6W02E4P8cNAPx6nLbvYQEAOB0ECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAkT7AAAMkTLABA8gQLAJA8wQIAJE+wAADJEywAQPIECwCQPMECACRvSMGyevXqqK6ujpKSkqipqYnt27f/yvnbtm2LmpqaKCkpicmTJ8d99903pMUCACNT3sGyYcOGaGxsjBUrVsSuXbti1qxZMXv27Ghraxtw/oEDB+K6666LWbNmxa5du+Kuu+6KJUuWxKOPPnrKiwcARoaCLMuyfE6YMWNGTJ8+PdasWZMbmzJlSsydOzeam5v7zV+2bFls2rQp9u7dmxtraGiI5557Lp5++ukBn6Onpyd6enpyt7u6umLSpElx8ODBKC0tzWe5AMBp0t3dHRMnTowjR45EWVnZqT1Yloeenp6ssLAw27hxY5/xJUuWZNdcc82A58yaNStbsmRJn7GNGzdmRUVF2bFjxwY8Z+XKlVlEOBwOh8PhOAOOn/zkJ/nkxoCKIg+HDh2K3t7eKC8v7zNeXl4eHR0dA57T0dEx4Pzjx4/HoUOHoqKiot85y5cvj6amptztI0eORFVVVbS1tZ16oY1wb9Wuq1Wnxj4OH3s5fOzl8LCPw+etV0jOO++8U36svILlLQUFBX1uZ1nWb+zt5g80/pbi4uIoLi7uN15WVuaXZ5iUlpbay2FgH4ePvRw+9nJ42MfhM2rUqX8oOa9HGD9+fBQWFva7mtLZ2dnvKspbJkyYMOD8oqKiGDduXJ7LBQBGoryCZcyYMVFTUxMtLS19xltaWqKurm7Ac2bOnNlv/tatW6O2tjZGjx6d53IBgJEo72s0TU1N8cADD8S6deti7969sXTp0mhra4uGhoaIePP9J/Pnz8/Nb2hoiJdeeimamppi7969sW7duli7dm3ccccdg37O4uLiWLly5YAvE5Efezk87OPwsZfDx14OD/s4fIZzL/P+WHPEm18c9/nPfz7a29tj6tSp8bd/+7dxzTXXRETEwoUL46c//Wl873vfy83ftm1bLF26NPbs2ROVlZWxbNmyXOAAALydIQULAMCvk78lBAAkT7AAAMkTLABA8gQLAJC8ZIJlzZo1cfnll+e+WXDmzJnxne98J3d/lmVx9913R2VlZYwdOzZ++7d/O/bs2XMaV/ybobm5OQoKCqKxsTE3Zi8H5+67746CgoI+x4QJE3L328f8vPLKK/HJT34yxo0bF2eddVZ88IMfjJ07d+but5+Dc+GFF/b7vSwoKIjFixdHhH0crOPHj8ef//mfR3V1dYwdOzYmT54cq1atihMnTuTm2MvBO3r0aDQ2NkZVVVWMHTs26urq4kc/+lHu/mHZy1P+a0TDZNOmTdm//uu/Zi+88EL2wgsvZHfddVc2evTo7D//8z+zLMuyz33uc9m5556bPfroo9nu3buzefPmZRUVFVl3d/dpXnm6nnnmmezCCy/MLr/88uz222/PjdvLwVm5cmV26aWXZu3t7bmjs7Mzd799HLxXX301q6qqyhYuXJj98Ic/zA4cOJB997vfzf7rv/4rN8d+Dk5nZ2ef38mWlpYsIrJ///d/z7LMPg7Wvffem40bNy779re/nR04cCB75JFHsnPOOSf74he/mJtjLwfvhhtuyC655JJs27Zt2b59+7KVK1dmpaWl2csvv5xl2fDsZTLBMpDf+q3fyh544IHsxIkT2YQJE7LPfe5zufv+93//NysrK8vuu+++07jCdB09ejS7+OKLs5aWluzDH/5wLljs5eCtXLkymzZt2oD32cf8LFu2LLv66qtPer/9HLrbb789u+iii7ITJ07YxzzMmTMnW7RoUZ+xT3ziE9knP/nJLMv8TubjjTfeyAoLC7Nvf/vbfcanTZuWrVixYtj2MpmXhP6/3t7e+Na3vhWvv/56zJw5Mw4cOBAdHR1RX1+fm1NcXBwf/vCHY8eOHadxpelavHhxzJkzJz72sY/1GbeX+dm3b19UVlZGdXV13HjjjbF///6IsI/52rRpU9TW1sYf/MEfxPnnnx9XXHFFfO1rX8vdbz+H5tixY/Hwww/HokWLoqCgwD7m4eqrr45/+7d/ixdffDEiIp577rl46qmn4rrrrosIv5P5OH78ePT29kZJSUmf8bFjx8ZTTz01bHuZVLDs3r07zjnnnCguLo6GhoZ47LHH4pJLLsn98cRf/gOL5eXl/f6wIhHf+ta34sc//nE0Nzf3u89eDt6MGTPioYceii1btsTXvva16OjoiLq6ujh8+LB9zNP+/ftjzZo1cfHFF8eWLVuioaEhlixZEg899FBE+L0cqscffzyOHDkSCxcujAj7mI9ly5bFTTfdFB/4wAdi9OjRccUVV0RjY2PcdNNNEWEv83HuuefGzJkz4y/+4i/iZz/7WfT29sbDDz8cP/zhD6O9vX3Y9rJoWFd9it7//vdHa2trHDlyJB599NFYsGBBbNu2LXd/QUFBn/lZlvUbG+kOHjwYt99+e2zdurVf7f5/9vLtzZ49O/fvyy67LGbOnBkXXXRR/MM//EN86EMfigj7OFgnTpyI2tra+Ku/+quIiLjiiitiz549sWbNmj5/e8x+5mft2rUxe/bsqKys7DNuH9/ehg0b4uGHH45vfOMbcemll0Zra2s0NjZGZWVlLFiwIDfPXg7OP/7jP8aiRYviggsuiMLCwpg+fXrcfPPN8eMf/zg351T3MqkrLGPGjIn3vve9UVtbG83NzTFt2rT40pe+lPtkxi+XWGdnZ79iG+l27twZnZ2dUVNTE0VFRVFUVBTbtm2LL3/5y1FUVJTbL3uZv7PPPjsuu+yy2Ldvn9/JPFVUVMQll1zSZ2zKlCnR1tYWEWE/h+Cll16K7373u3Hbbbflxuzj4N15553xmc98Jm688ca47LLL4pZbbomlS5fmrkzby/xcdNFFsW3btnjttdfi4MGD8cwzz8QvfvGLqK6uHra9TCpYflmWZdHT05P7gVtaWnL3HTt2LLZt2xZ1dXWncYXpufbaa2P37t3R2tqaO2pra+MP//APo7W1NSZPnmwvh6inpyf27t0bFRUVfifzdNVVV8ULL7zQZ+zFF1+MqqqqiAj7OQRf//rX4/zzz485c+bkxuzj4L3xxhsxalTf/wILCwtzH2u2l0Nz9tlnR0VFRfz3f/93bNmyJX7v935v+PbylN8ePEyWL1+ePfnkk9mBAwey//iP/8juuuuubNSoUdnWrVuzLHvzI1FlZWXZxo0bs927d2c33XSTj5cN0v//lFCW2cvB+vSnP51973vfy/bv35/94Ac/yH73d383O/fcc7Of/vSnWZbZx3w888wzWVFRUfaXf/mX2b59+7L169dnZ511Vvbwww/n5tjPwevt7c0mTZqULVu2rN999nFwFixYkF1wwQW5jzVv3LgxGz9+fPZnf/ZnuTn2cvCeeOKJ7Dvf+U62f//+bOvWrdm0adOyK6+8Mjt27FiWZcOzl8kEy6JFi7KqqqpszJgx2bvf/e7s2muvzcVKlr35EbOVK1dmEyZMyIqLi7Nrrrkm271792lc8W+OXw4Wezk4b31PwOjRo7PKysrsE5/4RLZnz57c/fYxP//yL/+STZ06NSsuLs4+8IEPZPfff3+f++3n4G3ZsiWLiOyFF17od599HJzu7u7s9ttvzyZNmpSVlJRkkydPzlasWJH19PTk5tjLwduwYUM2efLkbMyYMdmECROyxYsXZ0eOHMndPxx7WZBlWTb8F4UAAIZP0u9hAQCIECwAwG8AwQIAJE+wAADJEywAQPIECwCQPMECACRPsAAAyRMsAEDyBAsAkDzBAgAk7/8AnfbuPYXm5gwAAAAASUVORK5CYII="
          },
          "metadata": {}
        }
      ],
      "execution_count": 20
    },
    {
      "cell_type": "markdown",
      "source": "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\".",
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": "import numpy as np\n\n# Scatter plot of observed frequencies\nplt.scatter(data['Temperature'], data['Frequency'], label=\"Observed\")\n\n# Logistic curve from the fitted model\nx = np.linspace(30, 90, 200)\nX = np.column_stack([np.ones_like(x), x])\ny = logmodel.predict(X)\n\nplt.plot(x, y, color='red', label=\"Logistic fit\")\n\nplt.xlabel(\"Temperature\")\nplt.ylabel(\"Failure frequency\")\nplt.xlim(30, 90)\nplt.ylim(0, 1)\nplt.legend()\nplt.show()\n",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<Figure size 640x480 with 1 Axes>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAG2CAYAAACXuTmvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABY1UlEQVR4nO3de3yP9f/H8cdnm52wOY5hZnOW80QjUckxJUJRiEiUUwfpNJR0UlKREDkk5RQljL7Opww5Jjk0h605bo5j2/X74/r5sLax2We7tn2e99vtutnn+lyH1+eyfJ5d1/tgMwzDQERERMSJuFhdgIiIiEh2UwASERERp6MAJCIiIk5HAUhEREScjgKQiIiIOB0FIBEREXE6CkAiIiLidBSARERExOkoAImIiIjTUQASERERp2NpAFqzZg1t27alVKlS2Gw2Fi5ceNt9Vq9eTUhICJ6engQHB/PVV19lfaEiIiKSp1gagC5evEitWrX44osv0rX94cOHad26NY0bN2b79u28/vrrDBgwgHnz5mVxpSIiIpKX2HLKZKg2m40FCxbQrl27NLcZOnQoixYtYt++ffZ1ffv25Y8//mDjxo3ZUKWIiIjkBW5WF5ARGzdupHnz5snWtWjRgilTpnDt2jXy5cuXYp/4+Hji4+Ptr5OSkjhz5gxFixbFZrNlec0iIiKSeYZhcP78eUqVKoWLS+YfYOWqABQdHU2JEiWSrStRogQJCQmcOnUKf3//FPuMHj2aESNGZFeJIiIikoWOHj1KmTJlMn2cXBWAgBR3ba4/wUvrbs6wYcMYMmSI/XVsbCxly5bl6NGj+Pj4ZF2hIiIi4jBxcXEEBARQsGBBhxwvVwWgkiVLEh0dnWxdTEwMbm5uFC1aNNV9PDw88PDwSLHex8dHAUhERCSXcVTzlVw1DlBoaCjh4eHJ1i1fvpx69eql2v5HREREJDWWBqALFy6wY8cOduzYAZjd3Hfs2EFkZCRgPr7q1q2bffu+ffvyzz//MGTIEPbt28c333zDlClTePnll60oX0RERHIpSx+Bbd26lfvvv9/++npbne7duzNt2jSioqLsYQggKCiIJUuWMHjwYL788ktKlSrFuHHj6NChQ7bXLiIiIrlXjhkHKLvExcXh6+tLbGys2gCJiORSiYmJXLt2zeoyxMHc3d3T7OLu6O/vXNUIWkREnJthGERHR3Pu3DmrS5Es4OLiQlBQEO7u7ll+LgUgERHJNa6HHz8/P7y9vTWgbR6SlJTEiRMniIqKomzZsln+d6sAJCIiuUJiYqI9/KQ19InkbsWLF+fEiRMkJCRkee/uXNUNXkREnNf1Nj/e3t4WVyJZ5fqjr8TExCw/lwKQiIjkKnrslXdl59+tApCIiIg4HQUgERGRHKJcuXKMHTvW6jIcJid/HgUgERGRbHD06FF69epFqVKlcHd3JzAwkIEDB3L69GmrS3NKCkAiIuJUEpMMNh48zU87jrPx4GkSk7J+POBDhw5Rr149/vrrL2bPns3ff//NV199xcqVKwkNDeXMmTNZXkNqEhMTSUpKsuTcVlMAEhERp7F0dxT3fvAbT07axMDvd/DkpE3c+8FvLN0dlaXn7d+/P+7u7ixfvpwmTZpQtmxZWrVqxYoVKzh+/DhvvPGGfdvz58/TpUsXChQoQKlSpfj888+THWv48OGULVsWDw8PSpUqxYABA+zvXb16lVdffZXSpUuTP39+GjRowKpVq+zvT5s2jUKFCvHzzz9TrVo1PDw8mDRpEp6enikGlxwwYABNmjSxv96wYQP33XcfXl5eBAQEMGDAAC5evGh/PyYmhrZt2+Ll5UVQUBCzZs1y0NXLGgpAIiLiFJbujuL5mduIir2SbH107BWen7kty0LQmTNnWLZsGf369cPLyyvZeyVLlqRr167MmTOH6zNTffTRR9SsWZNt27YxbNgwBg8eTHh4OABz587l008/ZeLEiRw4cICFCxdSo0YN+/GeeeYZ1q9fz/fff8/OnTvp2LEjLVu25MCBA/ZtLl26xOjRo5k8eTJ79uzhqaeeolChQsybN8++TWJiIj/88ANdu3YFYNeuXbRo0YL27duzc+dO5syZw7p163jhhRfs+/To0YMjR47w22+/MXfuXMaPH09MTIzjL6iDaCBEERHJ8xKTDEYs3ktqD7sMwAaMWLyXh6qVxNXFsV2xDxw4gGEYVK1aNdX3q1atytmzZzl58iQAjRo14rXXXgOgUqVKrF+/nk8//ZSHHnqIyMhISpYsSbNmzciXLx9ly5alfv36ABw8eJDZs2dz7NgxSpUqBcDLL7/M0qVLmTp1Ku+99x5gjqc0fvx4atWqZa+hc+fOfPfdd/Tq1QuAlStXcvbsWTp27AiYoaxLly4MGjQIgIoVKzJu3DiaNGnChAkTiIyM5Ndff2XTpk00aNAAgClTpqT5mXMC3QESEZE8b8vhMynu/NzMAKJir7DlcPa3xbl+5+f6GDihoaHJ3g8NDWXfvn0AdOzYkcuXLxMcHEzv3r1ZsGABCQkJAGzbtg3DMKhUqRIFChSwL6tXr+bgwYP247m7u1OzZs1k5+jatSurVq3ixIkTAMyaNYvWrVtTuHBhACIiIpg2bVqy47Zo0YKkpCQOHz7Mvn37cHNzo169evZjVqlShUKFCjnwSjmW7gCJiEieF3M+7fBzJ9tlRIUKFbDZbOzdu5d27dqleP/PP/+kcOHCFCtWLM1jXA9HAQEB7N+/n/DwcFasWEG/fv346KOPWL16NUlJSbi6uhIREYGrq2uy/QsUKGD/2cvLK8WAg/Xr16d8+fJ8//33PP/88yxYsICpU6fa309KSuK5555L1t7ourJly7J///5kdeYGCkAiIpLn+RX0dOh2GVG0aFEeeughxo8fz+DBg5O1A4qOjmbWrFl069bNHh42bdqUbP9NmzZRpUoV+2svLy8eeeQRHnnkEfr370+VKlXYtWsXderUITExkZiYGBo3bpzhOrt06cKsWbMoU6YMLi4utGnTxv5e3bp12bNnDxUqVEh136pVq5KQkMDWrVvtj+T279+fomF1TqJHYCIikufVDyqCv68nad2fsAH+vp7UDyqSJef/4osviI+Pp0WLFqxZs4ajR4+ydOlSHnroIUqXLs2oUaPs265fv54PP/yQv/76iy+//JIff/yRgQMHAmYvrilTprB7924OHTrEjBkz8PLyIjAwkEqVKtG1a1e6devG/PnzOXz4ML///jsffPABS5YsuW2NXbt2Zdu2bYwaNYrHH38cT88bYXDo0KFs3LiR/v37s2PHDg4cOMCiRYt48cUXAahcuTItW7akd+/ebN68mYiICJ599tkUjb5zEgUgERHJ81xdbIS1rQaQIgRdfx3WtprDG0BfV7FiRbZu3Ur58uXp3Lkz5cuXp0+fPtx///1s3LiRIkVuBK+XXnqJiIgI6tSpwzvvvMOYMWNo0aIFAIUKFWLSpEk0atSImjVrsnLlShYvXkzRokUBmDp1Kt26deOll16icuXKPPLII2zevJmAgIB01Xj33Xezc+dOe++v62rWrMnq1as5cOAAjRs3pk6dOrz11lv4+/vbt5k6dSoBAQE0adKE9u3b06dPH/z8/Bxx+bKEzbje+spJxMXF4evrS2xsLD4+PlaXIyIi6XTlyhUOHz5MUFBQsrsTGbF0dxQjFu9N1iDa39eTsLbVaFnd/xZ7Sna41d+xo7+/1QZIREScRsvq/jxUrSRbDp8h5vwV/Aqaj72y6s6P5FwKQCIi4lRcXWyEli9qdRliMbUBEhEREaejACQiIiJORwFIREREnI4CkIiIiDgdBSARERFxOgpAIiIi4nQUgERERMTpKACJiIjkcuXKlWPs2LF3vP+0adMoVKhQpmr4+uuvCQgIwMXFhbFjxzJ8+HBq166dqWNmJU2FISIiuYIjpsKwQo8ePTh37hwLFy7MsnOcPHmS/Pnz4+3tfdtty5Urx6BBgxg0aJB93eXLlzl//vwdz90VFxdHsWLF+OSTT+jQoQO+vr4kJSURHx9vn6csPddBU2GIiIhIuhUvXjxT+3t5eWVq5vbIyEiuXbtGmzZtkk2QWqBAgUzVlZX0CExERMRCq1evpn79+nh4eODv789rr71GQkKC/f3z58/TtWtX8ufPj7+/P59++ilNmzZNdgfnv4/Ahg8fTtmyZfHw8KBUqVIMGDAAgKZNm/LPP/8wePBgbDYbNps5B1pqj8AWLVpEvXr18PT0pFixYrRv3z7V+qdNm0aNGjUACA4OxmazceTIkWSPwIYPH863337LTz/9ZD/vqlWrMnfhMkl3gEREJHcyDLh0yZpze3uDLfMTqB4/fpzWrVvTo0cPpk+fzp9//knv3r3x9PRk+PDhAAwZMoT169ezaNEiSpQowdtvv822bdvSbF8zd+5cPv30U77//nvuuusuoqOj+eOPPwCYP38+tWrVok+fPvTu3TvNun755Rfat2/PG2+8wYwZM7h69Sq//PJLqtt27tyZgIAAmjVrxpYtWwgICEhxR+rll19m3759xMXFMXXqVACKFCmSwavlWApAIiKSO126BFY9YrlwAfLnz/Rhxo8fT0BAAF988QU2m40qVapw4sQJhg4dyttvv83Fixf59ttv+e6773jwwQcBmDp1KqVKlUrzmJGRkZQsWZJmzZqRL18+ypYtS/369QEzdLi6ulKwYEFKliyZ5jFGjRrFE088wYgRI+zratWqleq2Xl5e9nY+xYsXT/W4BQoUwMvLi/j4+FueNzvpEZiIiIhF9u3bR2hoqP1RFECjRo24cOECx44d49ChQ1y7ds0eYAB8fX2pXLlymsfs2LEjly9fJjg4mN69e7NgwYJkj9TSY8eOHfbAlVfpDpCIiORO3t7mnRirzu0AhmEkCz/X1wHYbLZkP6e2TWoCAgLYv38/4eHhrFixgn79+vHRRx+xevVq8uXLl666MtMgOrfQHSAREcmdbDbzMZQViwPa/wBUq1aNDRs2JAs0GzZsoGDBgpQuXZry5cuTL18+tmzZYn8/Li6OAwcO3PK4Xl5ePPLII4wbN45Vq1axceNGdu3aBYC7uzuJiYm33L9mzZqsXLkyE58spfScNzvpDpCIiEgWi42NZceOHcnWFSlShH79+jF27FhefPFFXnjhBfbv309YWBhDhgzBxcWFggUL0r17d1555RWKFCmCn58fYWFhuLi4pLgrdN20adNITEykQYMGeHt7M2PGDLy8vAgMDATMHmNr1qzhiSeewMPDg2LFiqU4RlhYGA8++CDly5fniSeeICEhgV9//ZVXX331jq9BuXLlWLZsGfv376do0aL4+vqm+45UVtAdIBERkSy2atUq6tSpk2x5++23KV26NEuWLGHLli3UqlWLvn370qtXL9588037vp988gmhoaE8/PDDNGvWjEaNGlG1atU0B4MsVKgQkyZNolGjRvY7OYsXL7Y3VB45ciRHjhyhfPnyaY4f1LRpU3788UcWLVpE7dq1eeCBB9i8eXOmrkHv3r2pXLky9erVo3jx4qxfvz5Tx8ssjQQtIiK5Qm4dCdrRLl68SOnSpRkzZgy9evWyuhyH0kjQIiIiAsD27dv5888/qV+/PrGxsYwcORKARx991OLKcjcFIBERkRzu448/Zv/+/bi7uxMSEsLatWtTbbsj6acAJCIikoPVqVOHiIgIq8vIc9QIWkRERJyOApCIiOQqTtZ3x6lk59+tApCIiOQK18eMuWTVBKiS5a5evQqAq6trlp9LbYBERCRXcHV1pVChQsTExADg7e2d5mCAkvskJSVx8uRJvL29cXPL+niiACQiIrnG9ZnEr4cgyVtcXFwoW7ZstgRbBSAREck1bDYb/v7++Pn5ce3aNavLEQdzd3fHxSV7WucoAImISK7j6uqaLe1EJO9SI2gRERFxOgpAIiIi4nQUgERERMTpKACJiIiI03HeALR7t9UViIiIiEWcNwDdfz989RVoSHURERGn47wB6OpVeP556NQJzp2zuhoRERHJRs4bgEaNAjc3mDsX6tSBLVusrkhERESyifMGoBdegPXrISgIjhyBRo1gzBhISrK6MhEREclizhuAAOrXh+3boWNHSEiAl1+Gtm3h1CmrKxMREZEs5NwBCMDXF+bMMRtEe3jAkiV6JCYiIpLHKQAB2Gzw3HNm6KlUCY4dg8aNYfJkqysTERGRLKAAdLOaNeH33+HRR81eYr17Q58+EB9vdWUiIiLiQApA/+XjA/Pnm73EbDaYNAnuu8+8KyQiIiJ5guUBaPz48QQFBeHp6UlISAhr16695fazZs2iVq1aeHt74+/vzzPPPMPp06cdW5SLC7z+Ovz6KxQubD4aq1sXVq1y7HlERETEEpYGoDlz5jBo0CDeeOMNtm/fTuPGjWnVqhWRkZGpbr9u3Tq6detGr1692LNnDz/++CO///47zz77bNYU2KIFbN0KtWrByZPQrBl8+qlGjxYREcnlLA1An3zyCb169eLZZ5+latWqjB07loCAACZMmJDq9ps2baJcuXIMGDCAoKAg7r33Xp577jm2bt2adUUGB8OGDdC1KyQmwpAh8MwzahckIiKSi1kWgK5evUpERATNmzdPtr558+Zs2LAh1X0aNmzIsWPHWLJkCYZh8O+//zJ37lzatGmT5nni4+OJi4tLtmSYtzfMmAFjx5qPx779Fh56SOMFiYiI5FKWBaBTp06RmJhIiRIlkq0vUaIE0dHRqe7TsGFDZs2aRefOnXF3d6dkyZIUKlSIzz//PM3zjB49Gl9fX/sSEBBwZwXbbDBwIPzyi9lQeu1aaNAA9u27s+OJiIiIZSxvBG2z2ZK9Ngwjxbrr9u7dy4ABA3j77beJiIhg6dKlHD58mL59+6Z5/GHDhhEbG2tfjh49mrmCW7aEjRvNKTQOHYLQUFi+PHPHFBERkWzlZtWJixUrhqura4q7PTExMSnuCl03evRoGjVqxCuvvAJAzZo1yZ8/P40bN+bdd9/F398/xT4eHh54eHg4tvhq1WDzZnjsMXM+sdatYdw46NfPsecRERGRLGHZHSB3d3dCQkIIDw9Ptj48PJyGDRumus+lS5dwcUlesqurK2DeOcpWxYvDypXQrZvZOLp/fxgwwJxTTERERHI0Sx+BDRkyhMmTJ/PNN9+wb98+Bg8eTGRkpP2R1rBhw+jWrZt9+7Zt2zJ//nwmTJjAoUOHWL9+PQMGDKB+/fqUKlUq+z+AhwdMmwbvvWe+/vxzeOQRuHAh+2sRERGRdLPsERhA586dOX36NCNHjiQqKorq1auzZMkSAgMDAYiKiko2JlCPHj04f/48X3zxBS+99BKFChXigQce4IMPPrDqI5iNo4cNM+cQe/ppc/DE++83G0v7+VlXl4iIiKTJZmT7syNrxcXF4evrS2xsLD4+Po49+KZN8PDDcPo0VKgAS5dC+fKOPYeIiIgTcvT3t+W9wPKUe+4xG0WXKwd//w0NG0JEhNVViYiIyH8oADla5crmyNG1a0NMDDRpAsuWWV2ViIiI3EQBKCv4+8Pq1ebcYRcvmo/Fpk+3uioRERH5fwpAWcXHx2wI3aWL2TW+e3d4/31NpCoiIpIDKABlJXd3cw6x/x+4kWHDzJ8VgkRERCylAJTVXFzgww/hk0/M12PGQN++5uCJIiIiYgkFoOwyeDBMmWIGoq+/NscMunbN6qpERESckgJQdurZE2bPBjc3888OHeDKFaurEhERcToKQNmtUyf46Sfw9ITFi6FNG02dISIiks0UgKzQurU5ZUaBAvDbb/DQQ3D2rNVViYiIOA0FIKs0bQorVkDhwuYUGvffbw6cKCIiIllOAchKDRrAqlVQogT88Yc5anR0tNVViYiI5HkKQFarWRPWroUyZeDPP807QVFRVlclIiKSpykA5QQVK5p3ggICFIJERESygQJQTlG+/I0QtH+/QpCIiEgWUgDKSYKDzRBUtqwZgpo2hRMnrK5KREQkz1EAymluDkF//WXeCVIIEhERcSgFoJwoKMgMQYGBZghq2hSOH7e6KhERkTxDASinujkEHTigO0EiIiIOpACUk5UrZ4agcuXMENSsGZw8aXFRIiIiuZ8CUE5Xrhz873/mOEH79kHz5po2Q0REJJMUgHKDcuVg5UpzxOgdO6BVKzh/3uqqREREci0FoNyiUiVz7rAiRWDzZmjbFi5dsroqERGRXEkBKDepXh2WLwcfH1i9Gtq3h/h4q6sSERHJdRSAcpuQEFiyBLy9YdkyeOIJuHbN6qpERERyFQWg3KhRI1i0CDw8YOFC6N4dEhOtrkpERCTXUADKrR58EObNAzc3mD0b+vYFw7C6KhERkVxBASg3a9PGDD8uLjB5MrzxhtUViYiI5AoKQLnd44/DxInmz6NHw6efWluPiIhILqAAlBc8+6wZfgCGDIEZM6ytR0REJIdTAMorhg6FwYPNn595Bn75xdp6REREcjAFoLzCZoOPP4annzZ7hHXsCOvXW12ViIhIjqQAlJe4uMCUKdC6NVy+DA8/DLt2WV2ViIhIjqMAlNfkywc//ggNG8K5c9CiBRw+bHVVIiIiOYoCUF7k7Q0//2xOnREVZc4gHxNjdVUiIiI5hgJQXlW4MCxdCoGB8PffmjxVRETkJgpAeVnp0uZ8YUWKwJYt5rxhCQlWVyUiImI5BaC8rnLlG/OGLV4MAwZoygwREXF6CkDOoFEjmDXL7Co/YQJ8+KHVFYmIiFhKAchZdOhwY5qM116D776zth4RERELKQA5k4EDb4wW3aMH/O9/lpYjIiJiFQUgZ/Pxx+Yo0deuQbt2GihRRESckgKQs3FxgenT4d57IS7OHDX62DGrqxIREclWCkDOyNMTfvoJqlQxw8/DD8P581ZXJSIikm0UgJxVkSLw66/g5wd//AFPPmlOoioiIuIEMhyApk2bxiWNKJw3lCtnjhHk6Qm//AIvvWR1RSIiItkiwwFo2LBhlCxZkl69erFhw4asqEmyU4MGZpsggM8+gy+/tLYeERGRbJDhAHTs2DFmzpzJ2bNnuf/++6lSpQoffPAB0dHRWVGfZIeOHeG998yfBwww5xATERHJwzIcgFxdXXnkkUeYP38+R48epU+fPsyaNYuyZcvyyCOP8NNPP5GUlJQVtUpWeu01c2ygpCTo1End40VEJE/LVCNoPz8/GjVqRGhoKC4uLuzatYsePXpQvnx5Vq1a5aASJVvYbDBxIjRpYvYIe/hh0F09ERHJo+4oAP377798/PHH3HXXXTRt2pS4uDh+/vlnDh8+zIkTJ2jfvj3du3d3dK2S1dzdYf58qFQJIiPh0UdBDd5FRCQPshlGxqYGb9u2LcuWLaNSpUo8++yzdOvWjSJFiiTb5sSJE5QpUyZHPgqLi4vD19eX2NhYfHx8rC4nZ/r7b7Nx9Jkz8PjjMGeOOYCiiIiIRRz9/Z3hbzU/Pz9Wr17N7t27GTRoUIrwA+Dv78/hw4czXZxYpEIFWLgQ8uWDuXNh+HCrKxIREXGoDN8Byu10BygDvv3WbBgN5l2gTp0sLUdERJyX5XeABgwYwLhx41Ks/+KLLxg0aFCmC5IcpHv3G4Mj9ugB27ZZWo6IiIijZDgAzZs3j0aNGqVY37BhQ+bOneuQoiQH+eADaNUKLl82G0WrZ5iIiOQBGQ5Ap0+fxtfXN8V6Hx8fTp065ZCiJAdxdYXZs6FyZXPi1McegytXrK5KREQkUzIcgCpUqMDSVEYK/vXXXwkODnZIUZLD+PrC4sVQqBBs2gTPPQfO1XRMRETyGLeM7jBkyBBeeOEFTp48yQMPPADAypUrGTNmDGPHjnV0fZJTVKwIP/xgPg6bPh1q1tTkqSIikmvdUS+wCRMmMGrUKE6cOAFAuXLlGD58ON26dXN4gY6mXmCZNG4cDBxojgv0889mIBIREclijv7+zlQ3+JMnT+Ll5UWBAgUyXUh2UQDKJMOAPn1g8mTw8TEfiVWtanVVIiKSx1neDf5mxYsXz1XhRxzAZoMvv4TGjSEuDh55xBwxWkREJBfJcAD6999/efrppylVqhRubm64uromW8QJuLvDvHkQGGhOm9G5MyQkWF2ViIhIumU4APXo0YNt27bx1ltvMXfuXObPn59syajx48cTFBSEp6cnISEhrF279pbbx8fH88YbbxAYGIiHhwfly5fnm2++yfB5JZOKF4dFiyB/flixAoYMsboiERGRdMtwL7B169axdu1aateunemTz5kzh0GDBjF+/HgaNWrExIkTadWqFXv37qVs2bKp7tOpUyf+/fdfpkyZQoUKFYiJiSFBdx+sUbMmzJgB7dvD559D7drQs6fVVYmIiNxWhhtBV6tWjVmzZlGnTp1Mn7xBgwbUrVuXCRMm2NdVrVqVdu3aMXr06BTbL126lCeeeIJDhw6lOglreqgRdBYYORLCwsxHY2vXQv36VlckIiJ5jOWNoMeOHctrr73GkSNHMnXiq1evEhERQfPmzZOtb968ORs2bEh1n0WLFlGvXj0+/PBDSpcuTaVKlXj55Ze5fPlymueJj48nLi4u2SIO9uab0K4dXL1q3g3SdBkiIpLDZfgRWOfOnbl06RLly5fH29ubfPnyJXv/TDp7BJ06dYrExERKlCiRbH2JEiWITuML9NChQ6xbtw5PT08WLFjAqVOn6NevH2fOnEmzHdDo0aMZMWJEumqSO+TiYs4c36AB/PkndOwIK1ead4RERERyoAwHIEeP9myz2ZK9NgwjxbrrkpKSsNlszJo1yz4f2SeffMLjjz/Ol19+iZeXV4p9hg0bxpCbGujGxcUREBDgwE8ggDkm0MKF5uOvdetg8GCzu7yIiEgOlOEA1L17d4ecuFixYri6uqa42xMTE5PirtB1/v7+lC5dOtlkrFWrVsUwDI4dO0bFihVT7OPh4YGHh4dDapbbqFwZZs40xwYaPx5CQtQoWkREcqQ7Ggjx4MGDvPnmmzz55JPExMQAZgPlPXv2pPsY7u7uhISEEB4enmx9eHg4DRs2THWfRo0aceLECS5cuGBf99dff+Hi4kKZMmXu4JOIw7VtC9cfOT7/PGzZYm09IiIiqchwAFq9ejU1atRg8+bNzJ8/3x5Gdu7cSVhYWIaONWTIECZPnsw333zDvn37GDx4MJGRkfTt2xcwH1/dPL9Yly5dKFq0KM888wx79+5lzZo1vPLKK/Ts2TPVx19iETWKFhGRHC7DAei1117j3XffJTw8HPebGrnef//9bNy4MUPH6ty5M2PHjmXkyJHUrl2bNWvWsGTJEgIDAwGIiooiMjLSvn2BAgUIDw/n3Llz1KtXj65du9K2bVvGjRuX0Y8hWel6o+gqVeD4cbNR9NWrVlclIiJil+FxgAoUKMCuXbsICgqiYMGC/PHHHwQHB3PkyBGqVKnClStXsqpWh9A4QNlo/36zUXRcHPTrp0bRIiJyxywfB6hQoUJERUWlWL99+3ZKly6d6YIkD6lcGWbNMn8ePx40ZYmIiOQQGQ5AXbp0YejQoURHR2Oz2UhKSmL9+vW8/PLLydrriADw8MNqFC0iIjlOhh+BXbt2jR49evD9999jGAZubm4kJibSpUsXpk2bluNnhNcjMAskJUGHDuY4QaVLQ0QEpDHUgYiISGoc/f2d4QB03cGDB9m+fTtJSUnUqVMn1TF4ciIFIIvExd0YKbppUwgPB7cMD0MlIiJOKscEoNxKAchCf/4Jd98NFy7AK6/Ahx9aXZGIiOQSjv7+zvD/gve8zci+ac3JJUKVKjBtGjz+OHz0kXlHqEMHq6sSEREnlOFG0GfPnk22xMTE8NtvvzF//nzOnTuXBSVKntKhg3n3B6BHD/OukIiISDbL8B2gBQsWpFiXlJREv379CA4OdkhRkse99x78/jusWmWOFL15MxQsaHVVIiLiRO5oLrAUB3FxYfDgwXz66aeOOJzkdW5u8P33Zo+wffugVy9wrqZoIiJiMYcEIDB7hSUkJDjqcJLXlSgBP/4I+fKZfyo8i4hINsrwI7AhQ4Yke20YBlFRUfzyyy90797dYYWJEwgNhbFjoX9/ePVVCAmBJk2srkpERJxAhgPQ9u3bk712cXGhePHijBkz5rY9xERSeP552LQJZsyATp1g2zbz0ZiIiEgW0jhAYr1Ll8y7QTt3QsOG8L//gbu71VWJiEgOYvlkqCIO5+0N8+eDry9s2AAvv2x1RSIiksdl+BFYnTp1sNls6dp227ZtGS5InFT58jBzJrRtC59/bg6S2LWr1VWJiEgeleE7QC1btuTgwYN4eHjQtGlTmjZtiqenJwcPHqR58+Y8+uij9kUkQx5+GN56y/y5d2/zkZiIiEgWyPAdoJMnTzJgwADeeeedZOvDwsI4evSopsKQzAkLgy1bYNkyc9To33+HQoWsrkpERPKYDDeC9vX1ZevWrSlmfz9w4AD16tUjNjbWoQU6mhpB5wKnT5td4v/5x3wktnAhuKi5moiIM7O8EbSXlxfr1q1LsX7dunV4enpmuiARihaFefPAwwMWL4b337e6IhERyWMy/Ahs0KBBPP/880RERHDPPfcAsGnTJr755hvefvtthxcoTiokBMaPN6fJePNNqFcPmje3uioREckj7mgcoB9++IHPPvuMffv2AVC1alUGDhxIp06dHF6go+kRWC7Tpw9MmmTeFYqIgMBAqysSERELOPr7WwMhSs525Qo0bgxbt5p3hdatAz1qFRFxOpa3AQI4d+4ckydP5vXXX+fMmTOAOebP8ePHM12QSDKenjB37o07QAMGWF2RiIjkARkOQDt37qRSpUp88MEHfPTRR5w7dw6ABQsWMGzYMEfXJ2I+9vruO7DZzMdhGmpBREQyKcMBaMiQIfTo0YMDBw4k6/XVqlUr1qxZ49DiROyaN4frY0/162dOmioiInKHMhyAfv/9d5577rkU60uXLk10dLRDihJJ1bBh5rhA8fHmIIn///hVREQkozIcgDw9PYmLi0uxfv/+/RQvXtwhRYmkysUFpk+H4GA4cgSeegqSkqyuSkREcqEMB6BHH32UkSNHcu3aNQBsNhuRkZG89tprdOjQweEFiiRTqJA5c7ynJ/z6643HYiIiIhmQ4QD08ccfc/LkSfz8/Lh8+TJNmjShQoUKFCxYkFGjRmVFjSLJ1aoFEyeaP48YAUuXWluPiIjkOnc8DtBvv/3Gtm3bSEpKom7dujRr1szRtWUJjQOUhzz/PHz1FRQubHaRDwqyuiIREckilg6EmJCQgKenJzt27KB69eqZPrkVFIDykPh4uO8+c/b4unVh/XoNkigikkdZOhCim5sbgYGBJCYmZvrEIpnm4WEOklismNkt/oUXrK5IRERyiQy3AXrzzTcZNmyYfQRoEUsFBMDs2WYPsSlTzEVEROQ2MtwGqE6dOvz9999cu3aNwMBA8ufPn+z9bTl8gDo9AsujRo+G11837wqtX2/OGyYiInmGo7+/3TK6Q7t27TJ9UhGHGzoUNm2CRYvMQRIjIsz5w0RERFKRrjtA48aNo0+fPnh6ehIZGUmZMmVwcbmjeVQtpztAedi5c1CvHhw8CC1bws8/g6ur1VWJiIgDWNIIesiQIfbRn4OCgjh16lSmTyzicNcHSfTyMscG0iCJIiKShnQ9AitVqhTz5s2jdevWGIbBsWPHuHLlSqrbli1b1qEFimRIzZrw9dfw9NPmIIn160Pr1lZXJSIiOUy6HoF9/fXXvPjiiyQkJKS5jWEY2Gy2HN9FXo/AnET//jB+vAZJFBHJIywbCPH8+fP8888/1KxZkxUrVlA0jQamtWrVynRRWUkByElcvQpNmpgNo+vUMXuGeXlZXZWIiNwhy3qBFSxYkOrVqzN16lQaNWqEh4dHpk8ukmXc3eHHH80RordvN+8ITZkCNpvVlYmISA6Q4a5c3bt3V/iR3KFMGfj+e3OQxKlTYfJkqysSEZEcInf2ZRdJrwcegFGjzJ9feAF+/93aekREJEdQAJK8b+hQaNfObBf0+OOgYRxERJyeApDkfTYbTJsGFStCZCR07Qo5vLeiiIhkrTsOQFevXmX//v237BovkmP4+sK8eWZPsOXLzTGCRETEaWU4AF26dIlevXrh7e3NXXfdRWRkJAADBgzg/fffd3iBIg5TowZMmmT+/M475lQZIiLilDIcgIYNG8Yff/zBqlWr8PT0tK9v1qwZc+bMcWhxIg7XtavZGBrM0aIPHbK2HhERsUSGA9DChQv54osvuPfee7HdNKZKtWrVOHjwoEOLE8kSY8ZAaKg5eWqHDnD5stUViYhINstwADp58iR+fn4p1l+8eDFZIBLJsdzd4YcfwM8PduyA55+H9A2ILiIieUSGA9Ddd9/NL7/8Yn99PfRMmjSJ0NBQx1UmkpVuHiTx22/NCVRFRMRppHsqjOtGjx5Ny5Yt2bt3LwkJCXz22Wfs2bOHjRs3snr16qyoUSRr3H8/jB5tjhM0YIA5bcbdd1tdlYiIZIMM3wFq2LAhGzZs4NKlS5QvX57ly5dTokQJNm7cSEhISFbUKJJ1XnkFHnvMHCSxQwcNkigi4iTSPRs8wLVr1+jTpw9vvfUWwcHBWVlXltFs8JJCXJx55+evv6BZM1i6FFxdra5KRERu4ujv7wzdAcqXLx8LFizI9ElFchQfH5g/H7y9YcUKCAuzuiIREcliGX4E9thjj7Fw4cIsKEXEQnfddWO2+FGjYPFia+sREZEsleFG0BUqVOCdd95hw4YNhISEkD9//mTvDxgwwGHFiWSrJ5+ETZtg3DhzkMStW6FCBaurEhGRLJChNkAAQUFBaR/MZuNQDh9ZV22A5JauXjV7h23YANWrw8aNUKCA1VWJiDg9R39/Z/gO0OHDhzN9UpEcy90dfvwRQkJg92545hlz0EQN8ikikqfc8WzwInlWqVLmzPH58sHcufDBB1ZXJCIiDpbhO0A9e/a85fvffPPNHRcjkmM0bAhffAHPPQevvw61a0PLllZXJSIiDpLhAHT27Nlkr69du8bu3bs5d+4cDzzwgMMKE7Fcnz4QEWFOk/Hkk/D772oULSKSR2Q4AKU2DlBSUhL9+vXLtYMjiqRp3DjYtctsDN2unflnwYJWVyUiIpnkkDZALi4uDB48mE8//TTD+44fP56goCA8PT0JCQlh7dq16dpv/fr1uLm5Ubt27QyfUyTdPDzM9kD+/rBnj9koWjPHi4jkeg5rBH3w4EESEhIytM+cOXMYNGgQb7zxBtu3b6dx48a0atWKyMjIW+4XGxtLt27dePDBBzNTskj6+PvfaBQ9b545gaqIiORqGR4HaMiQIcleG4ZBVFQUv/zyC927d+eLL75I97EaNGhA3bp1mTBhgn1d1apVadeuHaNv8SXzxBNPULFiRVxdXVm4cCE7duxI9zk1DpDcsUmTzHZBNhv8/DO0bm11RSIiTsPycYC2b9+e7LWLiwvFixdnzJgxt+0hdrOrV68SERHBa6+9lmx98+bN2bBhQ5r7TZ06lYMHDzJz5kzefffd254nPj6e+Ph4++u4uLh01yiSTO/eZqPoiROhSxezUXTFilZXJSIidyDDAeh///ufQ0586tQpEhMTKVGiRLL1JUqUIDo6OtV9Dhw4wGuvvcbatWtxc0tf6aNHj2bEiBGZrlcEuNEoesMGs1H0pk1qFC0ikgtZPhCi7T8j7BqGkWIdQGJiIl26dGHEiBFUqlQp3ccfNmwYsbGx9uXo0aOZrlmcmLu7OTiivz/s3Qs9eqhRtIhILpSu2yh16tRJNZSkZtu2benarlixYri6uqa42xMTE5PirhDA+fPn2bp1K9u3b+eFF14AzO73hmHg5ubG8uXLUx2HyMPDAw8Pj3TVJJIu1xtFN2kC8+fDu+/CW29ZXZWIiGRAugJQu3btHH5id3d3QkJCCA8P57HHHrOvDw8P59FHH02xvY+PD7t27Uq2bvz48fz222/MnTv3lpO0ijhcaCiMH2+2C3r7bbjrLmjf3uqqREQkndIVgMLCwrLk5EOGDOHpp5+mXr16hIaG8vXXXxMZGUnfvn0B8/HV8ePHmT59Oi4uLlSvXj3Z/n5+fnh6eqZYL5Itnn3WbA80bhw8/TQEB5tTZoiISI6X4UbQjtS5c2dOnz7NyJEjiYqKonr16ixZsoTAwEAAoqKibjsmkIilxoyBffsgPBwefdTsGebnZ3VVIiJyG+kaB6hIkSL89ddfFCtWjMKFC9+yPdCZM2ccWqCjaRwgcbizZ6FBAzhwABo1gpUrzRGkRUTEYSwZB+jTTz+l4P939R07dmymTyqSpxQuDIsWwT33wPr18PzzMGWKOWCiiIjkSBkeCTq30x0gyTLLlpmjQyclwSefwODBVlckIpJnOPr7O1PjAF2+fJm4uLhki4jTatECPv7Y/Pnll81AJCIiOVKGA9DFixd54YUX8PPzo0CBAhQuXDjZIuLUBg2Cnj3Nu0CdO8Off1pdkYiIpCLDAejVV1/lt99+Y/z48Xh4eDB58mRGjBhBqVKlmD59elbUKJJ72Gzm+ECNGkFsLDzyiNlIWkREcpQMtwEqW7Ys06dPp2nTpvj4+LBt2zYqVKjAjBkzmD17NkuWLMmqWh1CbYAkW8TEwN13Q2QkNGsGv/4K6Zy/TkREUrK8DdCZM2fsoy77+PjYu73fe++9rFmzJtMFieQJfn5mzzBvb1ixAl56yeqKRETkJhkOQMHBwRw5cgSAatWq8cMPPwCwePFiChUq5MjaRHK3WrVg5kzz53HjYOJEa+sRERG7DAegZ555hj/++AMwp6q43hZo8ODBvPLKKw4vUCRXe+wxeOcd8+f+/WH5cmvrERERIANtgA4dOkRQUFCKUaAjIyPZunUr5cuXp1atWllSpCOpDZBkO8OAbt3Mu0E+PrBhgzl5qoiIpJtlbYAqVqzIyZMn7a87d+7Mv//+S9myZWnfvn2uCD8ilrDZYPJkaNwY4uKgTRv491+rqxIRcWrpDkD/vVG0ZMkSLl686PCCRPIkDw9YsAAqVIB//jEnTr182eqqREScVqZGghaRDChaFH75xZw7bPNm6N7dHDBRRESyXboDkM1mS9H+51azwotIKipVMu8E5csHP/4Ib72VLadNTDLYePA0P+04zsaDp0lMcqopACUP0u+0ZFa6R2YzDIMePXrg4eEBwJUrV+jbty/58+dPtt38+fMdW6FIXtOkCUyaBD16wHvvmY/Fnnkmy063dHcUIxbvJSr2in2dv68nYW2r0bK6f5adVySr6HdaHCHdvcCeSec/0FOnTs1UQVlNvcAkx3jzTRg1yhwhevlyuP9+h59i6e4onp+5jf/+R3793u2Ep+rqC0NyFf1OOy9Hf39neCqM3E4BSHKMpCTo0gXmzIFChWD9eqhWzWGHT0wyuPeD35L9X/LNbEBJX0/WDX0AVxc9zpacT7/Tzs3yqTBExEFcXGDqVAgNhXPnoFUrOHHCYYffcvhMml8UAAYQFXuFLYfPOOycIllJv9PiSApAIlby8jLnDKtUyZw4tXVrc6wgB4g5n/YXxZ1sJ2I1/U6LIykAiVitWDFztng/P/jjD3j8cbh2LdOH9Svo6dDtRKym32lxJAUgkZwgONgcIyh/fggPh2efNafQyIT6QUXw9/UkrZYQNsyeM/WDimTqPCLZRb/T4kgKQCI5Rb168MMP4OoK06fD229n6nCuLjbC2pqNqv/7hXH9dVjbamosKrmGfqfFkRSARHKS1q3hq6/Mn999F77+OlOHa1ndnwlP1aWkb/JHAiV9PdVdWHIl/U6Lo6gbvEhOFBYGI0eaPcV++gkefjhTh0tMMthy+Awx56/gV9B8RKD/S5bcTL/TzkfjAGWSApDkCoYBvXqZ3eS9veG336BBA6urEhGxjMYBEnEGNhtMnAgtW8KlS+ajsX37rK5KRCTPUAASyamuT5havz6cOQPNm5tjBYmISKYpAInkZAUKmN3jq1aFY8fMEHTypNVViYjkegpAIjldsWKwbBkEBMD+/ebjsPPnra5KRCRXUwASyQ0CAswZ44sWha1b4bHHID7e6qpERHItBSCR3KJKFXPKjAIFYOVK6NoVEhOtrkpEJFdSABLJTe6+GxYuBHd3mDcP+vXL9JQZIiLOSAFIJLd58EGYNcvsKv/11/Dmm1ZXJCKS6ygAieRGjz9+Y8qM996D0aOtrUdEJJdRABLJrfr0gQ8/NH9+/XX47DNr6xERyUUUgERys1deMecNAxg0KNOTp4qIOAsFIJHcLizMDEIAffvCjBnW1iMikgsoAInkdjYbfPAB9O9v9gjr0cOcQkNERNKkACSSF9hsMG4c9OwJSUnQpQv8/LPVVYmI5FgKQCJ5hYuL2QboySchIQE6dIAVK6yuSkQkR1IAEslLXF3h22+hXTu4ehUeeQRWr7a6KhGRHEcBSCSvyZcPvv8eWraEy5fNyVNXrbK6KhGRHEUBSCQv8vCABQvMEHTpkhmCfvvN6qpERHIMBSCRvMrT0wxBrVqZd4IefticRFVERBSARPK06yGodesbIUgNo0VEFIBE8jwPD5g/H9q0gStXoG1bWL7c6qpERCylACTiDDw8YN48M/xcuWL2Dlu2zOqqREQsowAk4iw8PGDuXHj0UYiPN//89VerqxIRsYQCkIgzcXeHH34wxwm6HoLmzbO6KhGRbKcAJOJsroegjh3h2jXo1MkcPFFExIkoAIk4o3z5YPZs6NXLnDusRw/44gurqxIRyTYKQCLOytUVJk2CwYPN1y++CKNGmTPKi4jkcQpAIs7MZoMxY2D4cPP1m2/C0KEKQSKS5ykAiTg7mw3CwuCTT8zXH30Ezz8PiYnW1iUikoUUgETENHgwTJ5sBqKJE+Hpp81G0iIieZACkIjc0KuXOZO8m5vZSLptWzh/3uqqREQcTgFIRJLr1AkWLQJvb3O06KZNITra6qpERBxKAUhEUmrVClatguLFYds2aNgQ/vrL6qpERBxGAUhEUnf33bBhA5QvD4cPmyFo0yarqxIRcQgFIBFJW4UKZgiqVw9On4YHHoDFi62uSkQk0xSAROTW/Pzgf/+D1q3h8mVzHrGvv7a6KhGRTFEAEpHbK1AAfvoJevY0p8547jlz0MSkJKsrExG5IwpAIpI+bm7mOEFvv22+HjUKOneGS5esrUtE5A5YHoDGjx9PUFAQnp6ehISEsHbt2jS3nT9/Pg899BDFixfHx8eH0NBQli1blo3Vijg5mw1GjICpU80JVefOhSZN4MQJqysTEckQSwPQnDlzGDRoEG+88Qbbt2+ncePGtGrVisjIyFS3X7NmDQ899BBLliwhIiKC+++/n7Zt27J9+/ZsrlzEyfXoAStXQtGisHWr2WNs2zarqxIRSTebYVg362GDBg2oW7cuEyZMsK+rWrUq7dq1Y/To0ek6xl133UXnzp15+/pt+duIi4vD19eX2NhYfHx87qhuEfl/hw6Zo0Xv3QteXjBjBnToYHVVIpIHOfr727I7QFevXiUiIoLmzZsnW9+8eXM2bNiQrmMkJSVx/vx5ihQpkuY28fHxxMXFJVtExEGCg81u8i1bmj3EHn/cbBuk2eRFJIezLACdOnWKxMRESpQokWx9iRIliE7nsPtjxozh4sWLdOrUKc1tRo8eja+vr30JCAjIVN0i8h++vubYQAMHmq/ffBO6dTMDkYhIDmV5I2ibzZbstWEYKdalZvbs2QwfPpw5c+bg5+eX5nbDhg0jNjbWvhw9ejTTNYvIf7i5wdix8NVX4OoKM2fCvffCkSNWVyYikirLAlCxYsVwdXVNcbcnJiYmxV2h/5ozZw69evXihx9+oFmzZrfc1sPDAx8fn2SLiGSR556D5cuhWDGzUXRICISHW12ViEgKlgUgd3d3QkJCCP/PP47h4eE0bNgwzf1mz55Njx49+O6772jTpk1WlykiGfXAAxARYU6fceaM2T5o9Gi1CxKRHMXSR2BDhgxh8uTJfPPNN+zbt4/BgwcTGRlJ3759AfPxVbdu3ezbz549m27dujFmzBjuueceoqOjiY6OJjY21qqPICKpKVsW1q6FXr3M0aJffx3atwd1QhCRHMLSANS5c2fGjh3LyJEjqV27NmvWrGHJkiUEBgYCEBUVlWxMoIkTJ5KQkED//v3x9/e3LwOvN74UkZzD09McOfrrr8HdHRYuhPr1zS7zIiIWs3QcICtoHCARC2zZYo4PdOyYOa/YlClwi96bIiL/5ejvbzcH1CQicmv165P4+1YuPPY4vpvWQefOJK1YgcvYseDtDUBiksGWw2eIOX8Fv4Ke1A8qgqvL7XuEZqfcUKMjXE1IYsbGI/xz5hKBRbx5OrQc7m6WdxoWcSjdARKRLLd0dxQjFu8l5uxFhqydyfOb5uKCwfkKlSm4cB5LjSKMWLyXqNgr9n38fT0Ja1uNltX9Laz8huufISfX6Aijl+xl0trDJN30zeBig96NgxjWupp1hYnTc/T3twKQiGSppbujeH7mNm7+h+bew9v59JcxFL94jmsenrx1f2++r9ncnGz1/13/acJTdS0PGKl9BshZNTrC6CV7mbjmcJrvP3efQpBYJ89MhSEieV9iksGIxXtTBId1QXVo9cznrClXh3zxV3h/6ed8vuhDCsZftG9zfZ8Ri/eSmGTd/6el9Rkg59ToCFcTkpi0Nu3wAzBp7WGuJiRlU0UiWUsBSESyzJbDZ5I9MrrZqfyF6d5pBO836UGCzYW2f67l52kDqRn1l30bA4iKvcKWw2eyqeKUbvUZIGfU6AgzNh7hdhkuyTC3E8kLFIBEJMvEnE87OAAYNhe+uudxOnX9gGM+fgSei2bezFfov2EOrkmJ6T5OVkrvua2s0RH+OXPJoduJ5HQKQCKSZfwKeqZru22lq9L6mXH8UrkR+ZISeWXtDObOfJXg08cydJyskN5zW1mjIwQW8XbodiI5nQKQiGSZ+kFF8Pf15FYdxV1sZmPiOM8C9H/0NQY9/BJxHvmpE7WfX6YN5MW9S6kfWCibKk7pdp/BhtkbrH5Qkewsy+GeDi3H7Xr0u9jM7UTyAgUgEckyri42wtqavYb++91q+/+ld+OgG+/bbCy8635a9PyCdYG18UqI56XFX+DasgUcPZqNld9wu88AENa2Wq4fD8jdzcX+d5GW3o2DNB6Q5Bn6TRaRLNWyuj8TnqpLSd/kj4hK+noy4am6DGtdLcX7UT7FebX3h+x9fRR4ecHKlVCjBkyfbsmkqrf7DHmhCzzAsNbVeO6+oBR3glxs6gIveY/GARKRbHG7UZTTfP+vv6B7d9i0ydywbVsYPx7KlMlxnyGv0EjQkhNpIMRMUgASyYUSEuDDD2H4cLh2DQoWhNGj4fnnwUVfzCLOQAMhiojzcXOD11+H7dvhnnvg/Hl44QW4917Ys8fq6kQkF1IAEpHc4667YN06+Pxzc1b5jRuhTh0IC4P4eKurE5FcRAFIRHIXV1fz7s/evWZ7oGvXYORIqF0b1q61ujoRySUUgEQkdwoIgJ9+gh9+gBIl4M8/4b77zAbT0dFWVyciOZwCkIjkXjYbdOwI+/ZBr17muunToVIl+OQT8+6QiEgqFIBEJPcrXBgmTza7yterZzaSfuklqFULVqywujoRyYEUgEQk72jQADZvhkmToFgx887QQw/B44/DP/9YXZ2I5CAKQCKSt7i4wLPPmgMovvii+XrePKhaFd5+27w7JCJOTwFIRPKmwoVh3Dhz7KD77oPLl+Gdd6BCBZgwQe2DRJycApCI5G01a8KqVfDjj2b4iYmBfv2genVYsMCSucVExHoKQCKS99lsZjugvXvhiy+geHHzEVn79uZo0uvXW12hiGQzBSARcR758kH//vD33/Dmm+ZM8xs2mCHoscdg506rKxSRbKIAJCLOx8fHbA904IDZYNrFBRYuNLvNd+wIu3dbXaGIZDEFIBFxXqVLm13md+40gw/A3Llmu6HOnc1HZiKSJykAiYjcdZc5pcbOndChg9kw+ocfzIbSXbqY02yISJ6iACQicl2NGuYdoB07zDZBhgGzZ5sBqXNn2LbN6gpFxEEUgERE/qtWLZg/3ww8jz4KSUnmHaGQEHNk6RUr1H1eJJdTABIRSUudOmbj6B07oGtXcHU1w89DD8Hdd5uhKDHR6ipF5A4oAImI3E6tWjBzptl9/sUXze7zERHmY7HKlc2RpS9etLpKEckABSARkfQqV86cXiMyEsLCoEgROHjQHFm6dGkYMsR8LSI5ngKQiEhGFSsGw4ebQWjcOHOKjdhY+PRTqFgR2raF5cvVTkgkB1MAEhG5U/nzm4/E9u+HJUugZUsz9Pz8M7RoYc5A/8UXEBdndaUi8h8KQCIimeXiAq1awa+/mmFowAAoWND8+cUXwd8fnnnGnHNMd4VEcgQFIBERR6pUCT77DI4fN+/+VKsGly7BtGnmnGPVqsGYMXDypNWVijg1BSARkaxQsKA58eru3eaEqz17gre3Oar0yy+bjaY7dYKlSyEhwepqRZyOzTCc635sXFwcvr6+xMbG4uPjY3U5IuJM4uLg++9h8mT4/fcb60uUgCeegKeeMgdbtNmsq1Ekh3L097cCkIiIFf74A6ZMMafaOHXqxvrKlc1BF7t2heBg6+oTyWEUgDJJAUhEcpRr18wu8zNnmqNOX7ly473QUHOwxfbtISDAshJFcgIFoExSABKRHCsuDhYsgFmzYOVKcw6y6+65Bx5/3Jytvlw5y0oUsYoCUCYpAIlIrnDihDnX2Lx5KbvPh4SYYah9e7PXmYgTUADKJAUgEcl1Tpww7wzNmwerVye/M1Spkjny9MMPm93s3dysq1MkCykAZZICkIjkajExZluhuXNh1SqzDdF1hQqZAzK2bWuOSl24sEVFijieAlAmKQCJSJ4RFwfLlplTb/zyC5w+feM9V1do0AAeegiaN4f69XV3SHI1BaBMUgASkTwpMRE2bTLD0OLFsGdP8vd9fODBB80w9NBDUL68NXWK3CEFoExSABIRp/DPPxAebnaxX7ECzp5N/n65ctC0qbk0aaKeZZLjKQBlkgKQiDidxETYts0MQ8uXm1Nz/Hf6jcDAG2GoSRMICtKI1JKjKABlkgKQiDi98+fNrvWrV5sNqX//3QxJNytZ0hyI8Z57zD/r1QMvL0vKFQEFoExTABIR+Y8LF8y7QqtWmaFoy5aUd4jc3KB2bTMMXV8CA3WXSLKNAlAmKQCJiNzG5csQEQEbN95YoqNTbleypDkoY926N/4sU0ahSLKEAlAmKQCJiGSQYZiNqjdtuhGItm9PeZcIoFgxMwhdX0JC1J5IHEIBKJMUgEREHODyZTMEbdt2Y9mzJ/VQ5OsL1aunXIoVy/66JddSAMokBSARkSxy5Qrs3p08FO3cCfHxqW9fokTyQHTXXVC5MhQpkr11S66gAJRJCkAiItno2jX4808zGN28HDqU9j7FiplBqHJlc66z6z8HB4OHR/bVLjmKAlAmKQCJiOQAFy7Avn2wa9eNULRnjznxa1pcXMz2RJUqmWEoKOjGn0FB5qM2ybMUgDJJAUhEJAe7cAH++stc9u83l+s/X7hw632LFLkRhm4ORsHBEBAAnp7Z8xkkSygAZZICkIhILmQYZlf8/fvhwAE4fNh8jHb9z1Onbn+MYsXMIFSmjLnc/PP1RYM95lgKQJmkACQikgedPw9HjtwIRf8NSJcvp+84RYveCEMlS6a9FCiQpR9HUlIAyiQFIBERJ2MY5mSwR4/CsWPmcvPP119fupT+Y3p7pwxFxYubd5n+uxQtqjtLDuDo7283B9QkIiKSc9lsZvugIkWgVq3UtzEMOHfuRhg6fhz+/dd87BYdnfznCxfMsHTo0K17s93M2zv1cHQ9IBUubC6FCiVf1G4pyygAiYiI2Gw3QkiNGrfe9sKF5IHo+nLq1I3l9OkbP1+7ZgamyEhzyQhPz5ShqFCh5GGpYMFbLwUKmD3oJBkFIBGRHCQxyWDL4TPEnL+CX0FP6gcVwdXlxjQSl68m8t6SvRw5fYlyRb15vXU1vNxdM3SM271/NSGJGRuP8M+ZSwQW8ebp0HK4u2XsCzSzNWT2+I74HGmeo0ABKFCAxKDg238GwzDbJ90cjm5akk6e5GxkFElnz+J18Tz5L1/AdvYsxMaa+165ciNgZUb+/LcPSt7et1/y50/+2tMz105zYnkboPHjx/PRRx8RFRXFXXfdxdixY2ncuHGa269evZohQ4awZ88eSpUqxauvvkrfvn3TfT61ARKRnGrp7ihGLN5LVOwV+zp/X0/C2lajZXV/ek//nfC9MSn2e6iaH5O63Z2uY9zu/dFL9jJp7WGSbvpmcLFB78ZBDGtdzSGf43bvZ/b4QKY/R1Z/htueo1oJMzidO2e2Xzp37sZy8+uzZ83t0lqSktJVS6akFpS8vMxBKz09HfZnXEICviEheaMR9Jw5c3j66acZP348jRo1YuLEiUyePJm9e/dStmzZFNsfPnyY6tWr07t3b5577jnWr19Pv379mD17Nh06dEjXORWARCQnWro7iudnbuO//yBf/3/rGmV82HksLs39H6rmR4e6ZW55jD73BfH1msNpvt+sml+qAeu65+67fXi43ee4XQ0Tnqp7ywBxu+NPeKou2yPPMnHN4Tv+HFn9GdL7OdIbpNJkGGbvt1sFpOvLxYvmY7r0LBcvwtWrmavtDsQBvpA3AlCDBg2oW7cuEyZMsK+rWrUq7dq1Y/To0Sm2Hzp0KIsWLWLfvn32dX379uWPP/5g48aN6TqnApCI5DSJSQb3fvBbsjsBd8KvQD5iLlxL9T0b5pOKpEz8i+9igz/faZXmY6T0fA6XW9RgA0r6erJu6AOpPg673fFtgF9Bd05euHrLz3mrz5HVnyG9n+N2x7BcQoIZrtIKSFeumEt8vMP+jLt8Gd+4uNzfC+zq1atERETw2muvJVvfvHlzNmzYkOo+GzdupHnz5snWtWjRgilTpnDt2jXy5cuXYp/4+Hjib5qILzY2FjCDkIhITrDl0BmOx5zJ9HGi05hz1FGSgK9X7KJbw6BU30/P57jdA5njMZf4385/qB+cckLU9Bw/Kv72Xdlv9Tmy+jOk9xy3O0aO4eVlLkWLZvmp4uLiICAAh923MSxy/PhxAzDWr1+fbP2oUaOMSpUqpbpPxYoVjVGjRiVbt379egMwTpw4keo+YWFhBqBFixYtWrRoyQPLwYMHHZJDLO8FZvtP63HDMFKsu932qa2/btiwYQwZMsT++ty5cwQGBhIZGYmvJs7LlLi4OAICAjh69KgeJ2aSrqVj6Do6jq6l4+haOkZsbCxly5alSBHH3BWzLAAVK1YMV1dXov/TtS8mJoYSJUqkuk/JkiVT3d7NzY2iadx+8/DwwMPDI8V6X19f/SI6iI+Pj66lg+haOoauo+PoWjqOrqVjuDhoTCPLRkZyd3cnJCSE8PDwZOvDw8Np2LBhqvuEhoam2H758uXUq1cv1fY/IiIiIqmxdGjIIUOGMHnyZL755hv27dvH4MGDiYyMtI/rM2zYMLp162bfvm/fvvzzzz8MGTKEffv28c033zBlyhRefvllqz6CiIiI5EKWtgHq3Lkzp0+fZuTIkURFRVG9enWWLFlCYGAgAFFRUUTeNGx4UFAQS5YsYfDgwXz55ZeUKlWKcePGpXsMIDAfiYWFhaX6WEwyRtfScXQtHUPX0XF0LR1H19IxHH0dLR8JWkRERCS7aXY0ERERcToKQCIiIuJ0FIBERETE6SgAiYiIiNPJswFowoQJ1KxZ0z7wVGhoKL/++qv9fcMwGD58OKVKlcLLy4umTZuyZ88eCyvOHUaPHo3NZmPQoEH2dbqW6TN8+HBsNluypWTJkvb3dR0z5vjx4zz11FMULVoUb29vateuTUREhP19Xc/bK1euXIrfSZvNRv/+/QFdw4xISEjgzTffJCgoCC8vL4KDgxk5ciRJSTdmDtP1TJ/z588zaNAgAgMD8fLyomHDhvz+++/29x12HR0yoUYOtGjRIuOXX34x9u/fb+zfv994/fXXjXz58hm7d+82DMMw3n//faNgwYLGvHnzjF27dhmdO3c2/P39jbi4OIsrz7m2bNlilCtXzqhZs6YxcOBA+3pdy/QJCwsz7rrrLiMqKsq+xMTE2N/XdUy/M2fOGIGBgUaPHj2MzZs3G4cPHzZWrFhh/P333/ZtdD1vLyYmJtnvY3h4uAEY//vf/wzD0DXMiHfffdcoWrSo8fPPPxuHDx82fvzxR6NAgQLG2LFj7dvoeqZPp06djGrVqhmrV682Dhw4YISFhRk+Pj7GsWPHDMNw3HXMswEoNYULFzYmT55sJCUlGSVLljTef/99+3tXrlwxfH19ja+++srCCnOu8+fPGxUrVjTCw8ONJk2a2AOQrmX6hYWFGbVq1Ur1PV3HjBk6dKhx7733pvm+ruedGThwoFG+fHkjKSlJ1zCD2rRpY/Ts2TPZuvbt2xtPPfWUYRj6nUyvS5cuGa6ursbPP/+cbH2tWrWMN954w6HXMc8+ArtZYmIi33//PRcvXiQ0NJTDhw8THR1N8+bN7dt4eHjQpEkTNmzYYGGlOVf//v1p06YNzZo1S7Ze1zJjDhw4QKlSpQgKCuKJJ57g0KFDgK5jRi1atIh69erRsWNH/Pz8qFOnDpMmTbK/r+uZcVevXmXmzJn07NkTm82ma5hB9957LytXruSvv/4C4I8//mDdunW0bt0a0O9keiUkJJCYmIinp2ey9V5eXqxbt86h1zFPB6Bdu3ZRoEABPDw86Nu3LwsWLKBatWr2CVX/O+lqiRIlUky2KvD999+zbds2Ro8eneI9Xcv0a9CgAdOnT2fZsmVMmjSJ6OhoGjZsyOnTp3UdM+jQoUNMmDCBihUrsmzZMvr27cuAAQOYPn06oN/LO7Fw4ULOnTtHjx49AF3DjBo6dChPPvkkVapUIV++fNSpU4dBgwbx5JNPArqe6VWwYEFCQ0N55513OHHiBImJicycOZPNmzcTFRXl0Oto6VQYWa1y5crs2LGDc+fOMW/ePLp3787q1avt79tstmTbG4aRYp2zO3r0KAMHDmT58uUpEvnNdC1vr1WrVvafa9SoQWhoKOXLl+fbb7/lnnvuAXQd0yspKYl69erx3nvvAVCnTh327NnDhAkTks0fqOuZflOmTKFVq1aUKlUq2Xpdw/SZM2cOM2fO5LvvvuOuu+5ix44dDBo0iFKlStG9e3f7drqetzdjxgx69uxJ6dKlcXV1pW7dunTp0oVt27bZt3HEdczTd4Dc3d2pUKEC9erVY/To0dSqVYvPPvvM3vPmv2kxJiYmRap0dhEREcTExBASEoKbmxtubm6sXr2acePG4ebmZr9eupYZlz9/fmrUqMGBAwf0O5lB/v7+VKtWLdm6qlWr2ucO1PXMmH/++YcVK1bw7LPP2tfpGmbMK6+8wmuvvcYTTzxBjRo1ePrppxk8eLD9zrmuZ/qVL1+e1atXc+HCBY4ePcqWLVu4du0aQUFBDr2OeToA/ZdhGMTHx9svYnh4uP29q1evsnr1aho2bGhhhTnPgw8+yK5du9ixY4d9qVevHl27dmXHjh0EBwfrWt6h+Ph49u3bh7+/v34nM6hRo0bs378/2bq//vrLPpGyrmfGTJ06FT8/P9q0aWNfp2uYMZcuXcLFJflXqqurq70bvK5nxuXPnx9/f3/Onj3LsmXLePTRRx17HTPVXDsHGzZsmLFmzRrj8OHDxs6dO43XX3/dcHFxMZYvX24YhtmNztfX15g/f76xa9cu48knn1R3xHS6uReYYehaptdLL71krFq1yjh06JCxadMm4+GHHzYKFixoHDlyxDAMXceM2LJli+Hm5maMGjXKOHDggDFr1izD29vbmDlzpn0bXc/0SUxMNMqWLWsMHTo0xXu6hunXvXt3o3Tp0vZu8PPnzzeKFStmvPrqq/ZtdD3TZ+nSpcavv/5qHDp0yFi+fLlRq1Yto379+sbVq1cNw3DcdcyzAahnz55GYGCg4e7ubhQvXtx48MEH7eHHMMwuiWFhYUbJkiUNDw8P47777jN27dplYcW5x38DkK5l+lwfqyJfvnxGqVKljPbt2xt79uyxv6/rmDGLFy82qlevbnh4eBhVqlQxvv7662Tv63qmz7JlywzA2L9/f4r3dA3TLy4uzhg4cKBRtmxZw9PT0wgODjbeeOMNIz4+3r6Nrmf6zJkzxwgODjbc3d2NkiVLGv379zfOnTtnf99R19FmGIbh2JtWIiIiIjmbU7UBEhEREQEFIBEREXFCCkAiIiLidBSARERExOkoAImIiIjTUQASERERp6MAJCIiIk5HAUhEREScjgKQiNySzWa75dKjRw+rS3S4pk2bMmjQIKvLEJEs5GZ1ASKSs0VFRdl/njNnDm+//XayiUi9vLysKOuOXLt2jXz58uXZ84lI+ukOkIjcUsmSJe2Lr68vNpst2bo1a9YQEhKCp6cnwcHBjBgxgoSEBPv+NpuNiRMn8vDDD+Pt7U3VqlXZuHEjf//9N02bNiV//vyEhoZy8OBB+z7Dhw+ndu3aTJw4kYCAALy9venYsSPnzp1LVtvUqVOpWrUqnp6eVKlShfHjx9vfO3LkCDabjR9++IGmTZvi6enJzJkzOX36NE8++SRlypTB29ubGjVqMHv2bPt+PXr0YPXq1Xz22Wf2u1xHjhxh2rRpFCpUKNn5Fy5ciM1mS1H3N998Q3BwMB4eHhiGQWxsLH369MHPzw8fHx8eeOAB/vjjDwf9DYnInVAAEpE7tmzZMp566ikGDBjA3r17mThxItOmTWPUqFHJtnvnnXfo1q0bO3bsoEqVKnTp0oXnnnuOYcOGsXXrVgBeeOGFZPv8/fff/PDDDyxevJilS5eyY8cO+vfvb39/0qRJvPHGG4waNYp9+/bx3nvv8dZbb/Htt98mO87QoUMZMGAA+/bto0WLFly5coWQkBB+/vlndu/eTZ8+fXj66afZvHkzAJ999hmhoaH07t2bqKgooqKiCAgISPc1uV73vHnz2LFjBwBt2rQhOjqaJUuWEBERQd26dXnwwQc5c+ZMuo8rIg7msOlbRSTPmzp1quHr62t/3bhxY+O9995Lts2MGTMMf39/+2vAePPNN+2vN27caADGlClT7Otmz55teHp62l+HhYUZrq6uxtGjR+3rfv31V8PFxcWIiooyDMMwAgICjO+++y7Zud955x0jNDTUMAzDOHz4sAEYY8eOve3nat26tfHSSy/ZXzdp0sQYOHDgLT+7YRjGggULjJv/GQ0LCzPy5ctnxMTE2NetXLnS8PHxMa5cuZJs3/LlyxsTJ068bW0ikjXUBkhE7lhERAS///57sjs+iYmJXLlyhUuXLuHt7Q1AzZo17e+XKFECgBo1aiRbd+XKFeLi4vDx8QGgbNmylClTxr5NaGgoSUlJ7N+/H1dXV44ePUqvXr3o3bu3fZuEhAR8fX2T1VivXr1krxMTE3n//feZM2cOx48fJz4+nvj4ePLnz5/ZywFAYGAgxYsXt7+OiIjgwoULFC1aNNl2ly9fTvbYT0SylwKQiNyxpKQkRowYQfv27VO85+npaf/55obA19vMpLYuKSkpzXNd38Zms9m3mzRpEg0aNEi2naura7LX/w02Y8aM4dNPP2Xs2LHUqFGD/PnzM2jQIK5evZr2BwVcXFwwDCPZumvXrqXY7r/nS0pKwt/fn1WrVqXY9r9tikQk+ygAicgdq1u3Lvv376dChQoOP3ZkZCQnTpygVKlSAGzcuBEXFxcqVapEiRIlKF26NIcOHaJr164ZOu7atWt59NFHeeqppwAzoBw4cICqVavat3F3dycxMTHZfsWLF+f8+fNcvHjRHnKut/G5lbp16xIdHY2bmxvlypXLUK0iknUUgETkjr399ts8/PDDBAQE0LFjR1xcXNi5cye7du3i3XffzdSxPT096d69Ox9//DFxcXEMGDCATp06UbJkScDscTVgwAB8fHxo1aoV8fHxbN26lbNnzzJkyJA0j1uhQgXmzZvHhg0bKFy4MJ988gnR0dHJAlC5cuXYvHkzR44coUCBAhQpUoQGDRrg7e3N66+/zosvvsiWLVuYNm3abT9Hs2bNCA0NpV27dnzwwQdUrlyZEydOsGTJEtq1a5fiEZ2IZA/1AhORO9aiRQt+/vlnwsPDufvuu7nnnnv45JNPCAwMzPSxK1SoQPv27WndujXNmzenevXqybq5P/vss0yePJlp06ZRo0YNmjRpwrRp0wgKCrrlcd966y3q1q1LixYtaNq0KSVLlqRdu3bJtnn55ZdxdXWlWrVqFC9enMjISIoUKcLMmTNZsmSJvev88OHDb/s5bDYbS5Ys4b777qNnz55UqlSJJ554giNHjtjbQ4lI9rMZ/32oLSJiseHDh7Nw4cJ0PWISEbkTugMkIiIiTkcBSERERJyOHoGJiIiI09EdIBEREXE6CkAiIiLidBSARERExOkoAImIiIjTUQASERERp6MAJCIiIk5HAUhEREScjgKQiIiIOB0FIBEREXE6/wfirNpum/D5pwAAAABJRU5ErkJggg=="
          },
          "metadata": {}
        }
      ],
      "execution_count": 21
    },
    {
      "cell_type": "code",
      "source": "import sys, platform\nimport numpy as np\nimport pandas as pd\nimport matplotlib\nimport statsmodels\n\nprint(\"=== SYSTEM ===\")\nprint(\"OS:\", platform.platform())\nprint(\"Python:\", sys.version.split()[0])\n\nprint(\"\\n=== LIB VERSIONS ===\")\nprint(\"numpy:\", np.__version__)\nprint(\"pandas:\", pd.__version__)\nprint(\"matplotlib:\", matplotlib.__version__)\nprint(\"statsmodels:\", statsmodels.__version__)\n\nprint(\"\\n=== MODEL RESULTS (GLM Binomial Logit) ===\")\nprint(\"params:\\n\", logmodel.params)\nprint(\"\\nstd errors:\\n\", logmodel.bse)\nprint(\"\\nAIC:\", logmodel.aic)\nprint(\"Deviance:\", logmodel.deviance)\nprint(\"Null deviance:\", logmodel.null_deviance)\n",
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "=== SYSTEM ===\nOS: Emscripten-4.0.9-wasm32-32bit\nPython: 3.13.2\n\n=== LIB VERSIONS ===\nnumpy: 2.2.5\npandas: 2.3.2\nmatplotlib: 3.8.4\nstatsmodels: 0.14.4\n\n=== MODEL RESULTS (GLM Binomial Logit) ===\nparams:\n Intercept      5.084977\nTemperature   -0.115601\ndtype: float64\n\nstd errors:\n Intercept      3.052484\nTemperature    0.047024\ndtype: float64\n\nAIC: 51.051946343967785\nDeviance: 18.086326742497448\nNull deviance: 24.230361814971374\n"
        }
      ],
      "execution_count": 22
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "source": "",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "execution_count": null
    }
  ]
}