{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this document we reperform some of the analysis provided in \n", "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n", "\n", "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Technical information on the computer on which the analysis is run" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#!pip install statsmodels \n", "#!pip install seaborn" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.6.4 |Anaconda, Inc.| (default, Mar 13 2018, 01:15:57) \n", "[GCC 7.2.0]\n", "uname_result(system='Linux', node='H2-SCHMIDER', release='5.3.0-46-generic', version='#38~18.04.1-Ubuntu SMP Tue Mar 31 04:17:56 UTC 2020', machine='x86_64', processor='x86_64')\n", "IPython 7.13.0\n", "IPython.core.release 7.13.0\n", "_csv 1.0\n", "_ctypes 1.1.0\n", "_curses b'2.2'\n", "decimal 1.70\n", "argparse 1.1\n", "backcall 0.1.0\n", "csv 1.0\n", "ctypes 1.1.0\n", "cycler 0.10.0\n", "dateutil 2.8.1\n", "decimal 1.70\n", "decorator 4.4.2\n", "distutils 3.6.4\n", "ipaddress 1.0\n", "ipykernel 5.1.4\n", "ipykernel._version 5.1.4\n", "ipython_genutils 0.2.0\n", "ipython_genutils._version 0.2.0\n", "ipywidgets 7.5.1\n", "ipywidgets._version 7.5.1\n", "jedi 0.16.0\n", "json 2.0.9\n", "jupyter_client 6.1.2\n", "jupyter_client._version 6.1.2\n", "jupyter_core 4.6.3\n", "jupyter_core.version 4.6.3\n", "kiwisolver 1.1.0\n", "logging 0.5.1.2\n", "matplotlib 3.1.3\n", "matplotlib.backends.backend_agg 3.1.3\n", "mkl 2.3.0\n", "numpy 1.18.1\n", "numpy.core 1.18.1\n", "numpy.core._multiarray_umath 3.1\n", "numpy.lib 1.18.1\n", "numpy.linalg._umath_linalg b'0.1.5'\n", "numpy.matlib 1.18.1\n", "optparse 1.5.3\n", "pandas 0.22.0\n", "_libjson 1.33\n", "parso 0.6.2\n", "patsy 0.5.1\n", "patsy.version 0.5.1\n", "pexpect 4.8.0\n", "pickleshare 0.7.5\n", "platform 1.0.8\n", "prompt_toolkit 3.0.4\n", "ptyprocess 0.6.0\n", "pygments 2.6.1\n", "pyparsing 2.4.6\n", "pytz 2019.3\n", "re 2.2.1\n", "scipy 1.1.0\n", "scipy._lib.decorator 4.0.5\n", "scipy._lib.six 1.2.0\n", "scipy.fftpack._fftpack b'$Revision: $'\n", "scipy.fftpack.convolve b'$Revision: $'\n", "scipy.integrate._dop b'$Revision: $'\n", "scipy.integrate._ode $Id$\n", "scipy.integrate._odepack 1.9 \n", "scipy.integrate._quadpack 1.13 \n", "scipy.integrate.lsoda b'$Revision: $'\n", "scipy.integrate.vode b'$Revision: $'\n", "scipy.interpolate._fitpack 1.7 \n", "scipy.interpolate.dfitpack b'$Revision: $'\n", "scipy.linalg 0.4.9\n", "scipy.linalg._fblas b'$Revision: $'\n", "scipy.linalg._flapack b'$Revision: $'\n", "scipy.linalg._flinalg b'$Revision: $'\n", "scipy.ndimage 2.0\n", "scipy.optimize._cobyla b'$Revision: $'\n", "scipy.optimize._lbfgsb b'$Revision: $'\n", "scipy.optimize._minpack 1.10 \n", "scipy.optimize._nnls b'$Revision: $'\n", "scipy.optimize._slsqp b'$Revision: $'\n", "scipy.optimize.minpack2 b'$Revision: $'\n", "scipy.signal.spline 0.2\n", "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", "scipy.special.specfun b'$Revision: $'\n", "scipy.stats.mvn b'$Revision: $'\n", "scipy.stats.statlib b'$Revision: $'\n", "seaborn 0.10.0\n", "seaborn.external.husl 2.1.0\n", "six 1.14.0\n", "statsmodels 0.9.0\n", "statsmodels.__init__ 0.9.0\n", "traitlets 4.3.3\n", "traitlets._version 4.3.3\n", "urllib.request 3.6\n", "zlib 1.0\n", "zmq 17.1.2\n", "zmq.sugar 17.1.2\n", "zmq.sugar.version 17.1.2\n" ] } ], "source": [ "def print_imported_modules():\n", " import sys\n", " for name, val in sorted(sys.modules.items()):\n", " if(hasattr(val, '__version__')): \n", " print(val.__name__, val.__version__)\n", "# else:\n", "# print(val.__name__, \"(unknown version)\")\n", "def print_sys_info():\n", " import sys\n", " import platform\n", " print(sys.version)\n", " print(platform.uname())\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import seaborn as sns\n", "\n", "print_sys_info()\n", "print_imported_modules()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading and inspecting data\n", "Let's start by reading data." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "file downloaded on Tue Apr 14 10:06:35 2020\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
56/18/82672500
68/30/836731000
711/28/836701000
82/03/846572001
94/06/846632001
108/30/846702001
1110/05/846782000
1211/08/846672000
131/24/856532002
144/12/856672000
154/29/856752000
166/17/856702000
177/2903/856812000
188/27/856762000
1910/03/856792000
2010/30/856752002
2111/26/856762000
221/12/866582001
\n", "
" ], "text/plain": [ " Date Count Temperature Pressure Malfunction\n", "0 4/12/81 6 66 50 0\n", "1 11/12/81 6 70 50 1\n", "2 3/22/82 6 69 50 0\n", "3 11/11/82 6 68 50 0\n", "4 4/04/83 6 67 50 0\n", "5 6/18/82 6 72 50 0\n", "6 8/30/83 6 73 100 0\n", "7 11/28/83 6 70 100 0\n", "8 2/03/84 6 57 200 1\n", "9 4/06/84 6 63 200 1\n", "10 8/30/84 6 70 200 1\n", "11 10/05/84 6 78 200 0\n", "12 11/08/84 6 67 200 0\n", "13 1/24/85 6 53 200 2\n", "14 4/12/85 6 67 200 0\n", "15 4/29/85 6 75 200 0\n", "16 6/17/85 6 70 200 0\n", "17 7/2903/85 6 81 200 0\n", "18 8/27/85 6 76 200 0\n", "19 10/03/85 6 79 200 0\n", "20 10/30/85 6 75 200 2\n", "21 11/26/85 6 76 200 0\n", "22 1/12/86 6 58 200 1" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")\n", "data_url ='https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/7e01583b99527ad27cbeae0f9d2085fe8f2f1d15/data/shuttle.csv?inline=false'\n", "data_file='/home/aschmide/Documents/formation_MOOC_RR/mooc-rr/module4/shuttle.csv'\n", "import os\n", "import urllib.request\n", "\n", "if not os.path.exists(data_file):\n", " urllib.request.urlretrieve(data_url, data_file)\n", " \n", "data = pd.read_csv(data_file)\n", "\n", "import time\n", "modtime=os.path.getmtime(data_file)\n", "modtime=time.ctime(modtime)\n", "print('file downloaded on',modtime)\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYMklEQVR4nO3df5TU9X3v8edrAWERIgRSagWiVuqNRy3KBjX2B0STg54TaK6a4D3RNC2h90ROTkyaaHtzreWm59zYJLa5sYnEaxrtSYhKo9xeev0RJak9/gCVgGKwWzW4YFA3qKwg7LLv+8d8txmG2eU7y35ndubzepyzZ+f7/X7mO+/PfHf2Nd8f8xlFBGZmlq62RhdgZmaN5SAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcYUEg6VZJr0h6epDlkvR1SZ2SNks6u6hazMxscEXuEfw9sGiI5RcBc7Kf5cA3C6zFzMwGUVgQRMRPgF8O0WQJcFuUPApMkXR8UfWYmVl1Yxv42CcAL5VNd2XzXq5sKGk5pb0G2tvb582aNasuBebV399PW1vrnW5p1X5B6/bN/Wo+9erbc88991pEvKvaskYGgarMqzreRUSsAlYBdHR0xMaNG4usq2br169nwYIFjS5jxLVqv6B1++Z+NZ969U3Szwdb1siI7QLK39rPBHY2qBYzs2Q1MgjWAldmVw+dC7wREYcdFjIzs2IVdmhI0veBBcB0SV3AXwDjACLiW8A64GKgE9gLfKKoWszMbHCFBUFEXH6E5QFcVdTjm5lZPq15Gt7MzHJzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWOAeBmVniHARmZolzEJiZJc5BYGaWuEKDQNIiSdskdUq6tsry2ZIekvSUpM2SLi6yHjMzO1xhQSBpDHATcBFwGnC5pNMqmn0RuCMizgKWAn9XVD1mZlZdkXsE84HOiHg+Ig4Aq4ElFW0CeEd2+zhgZ4H1mJlZFYqIYlYsXQosiohl2fQVwDkRsaKszfHAfcBU4Fjgwoh4osq6lgPLAWbMmDFv9erVhdQ8XD09PUyaNKnRZYy4Vu0XtG7f3K/mU6++LVy48ImI6Ki2bGyBj6sq8ypT53Lg7yPiq5LOA26XdHpE9B9yp4hVwCqAjo6OWLBgQRH1Dtv69esZbTWNhFbtF7Ru39yv5jMa+lbkoaEuYFbZ9EwOP/Tzx8AdABHxCDABmF5gTWZmVqHIINgAzJF0kqRjKJ0MXlvRZjtwAYCk91AKglcLrMnMzCoUFgQR0QesAO4FnqV0ddAzklZKWpw1+xzwSUk/Bb4P/GEUddLCzMyqKvIcARGxDlhXMe+6sttbgfOLrMHMzIbmTxabmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJcxCYmSXOQWBmlrhCg0DSIknbJHVKunaQNh+RtFXSM5K+V2Q9ZmZ2uLF5Gkk6PSKermXFksYANwEfALqADZLWRsTWsjZzgD8Dzo+I3ZJ+rZbHMDOzo5d3j+Bbkh6X9ClJU3LeZz7QGRHPR8QBYDWwpKLNJ4GbImI3QES8knPdZmY2QhQR+RqW3r3/EXAZ8DjwnYi4f4j2lwKLImJZNn0FcE5ErChrczfwHHA+MAa4PiL+X5V1LQeWA8yYMWPe6tWr8/WuTnp6epg0aVKjyxhxrdovaN2+uV/Np159W7hw4RMR0VF1YUTk/qH0z/oSYAfwLPAz4D8P0vYy4Jay6SuA/1XR5p+AHwLjgJMoHUKaMlQN8+bNi9HmoYceanQJhWjVfkW0bt/cr+ZTr74BG2OQ/6u5Dg1JOlPSjdk///cDH4qI92S3bxzkbl3ArLLpmcDOKm3uiYjeiHgB2AbMyVOTmZmNjLznCL4BPAn8dkRcFRFPAkTETuCLg9xnAzBH0kmSjgGWAmsr2twNLASQNB34LeD52rpgZmZHI9dVQ8DFwL6IOAggqQ2YEBF7I+L2aneIiD5JK4B7KR1SujUinpG0ktIuytps2QclbQUOAp+PiO6j7JOZmdUgbxA8AFwI9GTTE4H7gPcNdaeIWAesq5h3XdntAD6b/ZiZWQPkPTQ0ISIGQoDs9sRiSjIzs3rKGwRvSTp7YELSPGBfMSWZmVk95T009BngTkkDV/0cD3y0mJLMzKyecgVBRGyQ9J+AUwEBP4uI3kIrMzOzusi7RwDwXuDE7D5nSSIibiukKjMzq5u8g87dDvwmsInSZZ4AATgIzMyaXN49gg7gtOxyTzMzayF5rxp6Gvj1IgsxM7PGyLtHMB3YKulxYP/AzIhYXEhVZmZWN3mD4PoiizAzs8bJe/nojyW9G5gTEQ9Imkhp/CAzM2tyeYeh/iRwF3BzNusESiOHmplZk8t7svgqSt8i9iZARPwb4O8XNjNrAXmDYH+UvncYAEljKX2OwMzMmlzeIPixpD8H2iV9ALgT+D/FlWVmZvWSNwiuBV4FtgB/Quk7Bgb7ZjIzM2siea8a6ge+nf2YmVkLyTvW0AtUOScQESePeEVmZlZXtYw1NGACcBnwzpEvx8zM6i3XOYKI6C772RERfwO8v+DazMysDvIeGjq7bLKN0h7C5EIqMjOzusp7aOirZbf7gBeBj4x4NWZmVnd5rxpaWHQhZmbWGHkPDX12qOUR8bWRKcfMzOqtlquG3guszaY/BPwEeKmIoszMrH5q+WKasyNiD4Ck64E7I2JZUYWZmVl95B1iYjZwoGz6AHDiiFdjZmZ1l3eP4HbgcUk/pPQJ4w8DtxVWlZmZ1U3eq4b+StI/A7+bzfpERDxVXFlmZlYveQ8NAUwE3oyIvwW6JJ1UUE1mZlZHeb+q8i+Aa4A/y2aNA/6hqKLMzKx+8u4RfBhYDLwFEBE78RATZmYtIW8QHIiIIBuKWtKxxZVkZmb1lDcI7pB0MzBF0ieBB/CX1JiZtYS8Vw19Jfuu4jeBU4HrIuL+QiszM7O6OOIegaQxkh6IiPsj4vMR8ad5Q0DSIknbJHVKunaIdpdKCkkdg7UxM7NiHDEIIuIgsFfScbWsWNIY4CbgIuA04HJJp1VpNxn4NPBYLes3M7ORkfeTxW8DWyTdT3blEEBEfHqI+8wHOiPieQBJq4ElwNaKdv8DuAH407xFm5nZyMkbBP83+6nFCRw6OmkXcE55A0lnAbMi4p8kDRoEkpYDywFmzJjB+vXrayylWD09PaOuppHQqv2C1u2b+9V8RkPfhgwCSbMjYntEfHcY61aVeVG27jbgRuAPj7SiiFgFrALo6OiIBQsWDKOc4qxfv57RVtNIaNV+Qev2zf1qPqOhb0c6R3D3wA1Ja2pcdxcwq2x6JrCzbHoycDqwXtKLwLnAWp8wNjOrryMFQfm7+pNrXPcGYI6kkyQdAyzlV19sQ0S8ERHTI+LEiDgReBRYHBEba3wcMzM7CkcKghjk9hFFRB+wArgXeBa4IyKekbRS0uLayjQzs6Ic6WTxb0t6k9KeQXt2m2w6IuIdQ905ItYB6yrmXTdI2wW5KjYzsxE1ZBBExJh6FWJmZo1Ry/cRmJlZC3IQmJklzkFgZpY4B4GZWeKSCYLunv389KXX6e7Z3+hSzKwG3T372dd70K/dAiURBPds2sH5X36Qj93yGOd/+UHWbtrR6JLMLIeB1+4Lr77l126BWj4Iunv2c82azbzd28+e/X283dvPF9Zs9rsLs1Gu/LV7MMKv3QK1fBB07d7HuLZDuzmurY2u3fsaVJGZ5eHXbv20fBDMnNpOb3//IfN6+/uZObW9QRWZWR5+7dZPywfBtEnjueGSM5kwro3J48cyYVwbN1xyJtMmjW90aWY2hPLX7hjJr90C5f1imqa2eO4JnH/KdLp272Pm1Hb/IZk1iYHX7uOPPMy/Lv4dv3YLkkQQQOndhf+IzJrPtEnjaR83xq/fArX8oSEzMxuag8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tcoUEgaZGkbZI6JV1bZflnJW2VtFnSjyS9u8h6zMzscIUFgaQxwE3ARcBpwOWSTqto9hTQERFnAncBNxRVj5mZVVfkHsF8oDMino+IA8BqYEl5g4h4KCL2ZpOPAjMLrMfMzKpQRBSzYulSYFFELMumrwDOiYgVg7T/BvCLiPhSlWXLgeUAM2bMmLd69epCah6unp4eJk2a1OgyRlyr9gtat2/uV/OpV98WLlz4RER0VFs2tsDHVZV5VVNH0seADuD3qy2PiFXAKoCOjo5YsGDBCJU4MtavX89oq2kktGq/oHX75n41n9HQtyKDoAuYVTY9E9hZ2UjShcB/A34/IvYXWI+ZmVVR5DmCDcAcSSdJOgZYCqwtbyDpLOBmYHFEvFJgLWZmNojCgiAi+oAVwL3As8AdEfGMpJWSFmfN/hqYBNwpaZOktYOszszMClLkoSEiYh2wrmLedWW3Lyzy8ZtZd89+unbvY+bUdqZNGj9ibZtJq/arKJ279rB7by+du/ZwyozJjS7HmkihQWDDc8+mHVyzZjPj2tro7e/nhkvOZPHcE466bTNp1X4V5bq7t3Dbo9v53Bl9XH3jT7jyvNmsXHJGo8uyJuEhJkaZ7p79XLNmM2/39rNnfx9v9/bzhTWb6e45/Dx6LW2bSav2qyidu/Zw26PbD5l32yPb6dy1p0EVWbNxEIwyXbv3Ma7t0M0yrq2Nrt37jqptM2nVfhVl00uv1zTfrJKDYJSZObWd3v7+Q+b19vczc2r7UbVtJq3ar6LMnTWlpvlmlRwEo8y0SeO54ZIzmTCujcnjxzJhXBs3XHJm1ZOltbRtJq3ar6KcMmMyV543+5B5V5432yeMLTefLB6FFs89gfNPmZ7ripla2jaTVu1XUVYuOYMrzz2RLU88ygNXn+sQsJo4CEapaZPG5/7nV0vbZtKq/SrKKTMm0zVxnEPAauZDQ2ZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIKDQJJiyRtk9Qp6doqy8dL+kG2/DFJJxZZj1mtunv289OXXqe7Z/+Q7Ta+0M3X7tvGxhe6R2ydtbbt3LWH3Xt76dy154hta1FUvbU8/r7eg7mfg7s2vtRyz0GR6wUYO+JrzEgaA9wEfADoAjZIWhsRW8ua/TGwOyJOkbQU+DLw0aJqMqvFPZt2cM2azYxra6O3v58bLjmTxXNPOKzdx255lIc7SwHw9Qc7+d1TpnH7snOPap21tr3u7i3c9uh2PndGH1ff+BOuPG82K5ecMcyeF19vrY//6ff0cvWXH8z1HAxoleegyPUOKHKPYD7QGRHPR8QBYDWwpKLNEuC72e27gAskqcCazHLp7tnPNWs283ZvP3v29/F2bz9fWLP5sHdjG1/o/o8QGPAvnd1V9wzyrrPWtp279hzyDxDgtke2H/W74qLqHc7jH4xI8jkocr3lFBEjtrJDVixdCiyKiGXZ9BXAORGxoqzN01mbrmz637M2r1WsazmwPJs8FdhWSNHDNx147Yitmk+r9guO0DeNa584durxv6W2tjED86K//2Df7pefi959ewfmjZk8/TfGHDvl+Mr7H3zr9ZcP7nlt53DWWWvbtonHTRv7jnedCHBw7xuMmXgcAH1vvvpi/943jnys6iifg1rbDufxB/qV5zko1yTPwYj8Lebw7oh4V7UFhR0aAqq9s69MnTxtiIhVwKqRKKoIkjZGREej6xhprdovaN2+SdrY98Yr7lcTGQ1/i0UeGuoCZpVNzwR2DtZG0ljgOOCXBdZkZmYVigyCDcAcSSdJOgZYCqytaLMW+Hh2+1LgwSjqWJWZmVVV2KGhiOiTtAK4FxgD3BoRz0haCWyMiLXA/wZul9RJaU9gaVH1FGzUHrY6Sq3aL2jdvrlfzafhfSvsZLGZmTUHf7LYzCxxDgIzs8Q5CIZB0ouStkjaJGljNu96STuyeZskXdzoOmslaYqkuyT9TNKzks6T9E5J90v6t+z31EbXWatB+tUK2+vUsvo3SXpT0meafZsN0a9W2GZXS3pG0tOSvi9pQnZBzWPZ9vpBdnFNfevyOYLaSXoR6Cj/4Juk64GeiPhKo+o6WpK+C/xLRNyS/TFOBP4c+GVE/M9svKipEXFNQwut0SD9+gxNvr3KZUO67ADOAa6iybfZgIp+fYIm3maSTgAeBk6LiH2S7gDWARcD/xgRqyV9C/hpRHyznrV5j8AAkPQO4PcoXclFRByIiNc5dBiQ7wJ/0JgKh2eIfrWaC4B/j4if0+TbrEJ5v1rBWKA9+9zUROBl4P2UhtiBBm0vB8HwBHCfpCey4S8GrJC0WdKtzbY7DpwMvAp8R9JTkm6RdCwwIyJeBsh+/1ojixyGwfoFzb29Ki0Fvp/dbvZtVq68X9DE2ywidgBfAbZTCoA3gCeA1yOiL2vWBYzcaHI5OQiG5/yIOBu4CLhK0u8B3wR+E5hLaSN/tYH1DcdY4GzgmxFxFvAWcNjQ4U1osH41+/b6D9nhrsXAnY2uZSRV6VdTb7MsuJYAJwG/ARxL6X9Ipbofr3cQDENE7Mx+vwL8EJgfEbsi4mBE9APfpjT6ajPpAroi4rFs+i5K/0B3SToeIPv9SoPqG66q/WqB7VXuIuDJiNiVTTf7NhtwSL9aYJtdCLwQEa9GRC/wj8D7gCnZoSKoPhRP4RwENZJ0rKTJA7eBDwJPD7zwMh8Gnm5EfcMVEb8AXpJ0ajbrAmArhw4D8nHgngaUN2yD9avZt1eFyzn08ElTb7Myh/SrBbbZduBcSRMliV+9xh6iNMQONGh7+aqhGkk6mdJeAJQOO3wvIv5K0u2UdlkDeBH4k4HjtM1C0lzgFuAY4HlKV2m0AXcAsyn9IV8WEU01MOAg/fo6Tb69ACRNBF4CTo6IN7J502j+bVatX63wGvtLSl++1Qc8BSyjdE5gNfDObN7HImLkv4ZsqLocBGZmafOhITOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxBX55fVmdZVdNvmjbPLXgYOUhpeA0of+DjSksCFI+iNgXfZ5B7OG8OWj1pJG02iwksZExMFBlj0MrIiITTWsb2zZ2DRmR82HhiwJkj4u6fFsHPu/k9Qmaayk1yX9taQnJd0r6RxJP5b0/MB495KWSfphtnybpC/mXO+XJD0OzJf0l5I2ZOPQf0slH6X0AakfZPc/RlKXpCnZus+V9EB2+0uSbpZ0P6UB9MZK+lr22JslLav/s2qtwkFgLU/S6ZSGJHhfRMyldEh0abb4OOC+bBDBA8D1lD76fxmwsmw187P7nA38F0lzc6z3yYiYHxGPAH8bEe8FzsiWLYqIHwCbgI9GxNwch67OAj4UEVcAy4FXImI+8F5Kgx/OHs7zY+ZzBJaCCyn9s9xYGuKFdkrDFwDsi4j7s9tbgDciok/SFuDEsnXcGxG7ASTdDfwOpdfPYOs9wK+GIgG4QNLngQnAdErDD/9zjf24JyLezm5/EHiPpPLgmUNpSAmzmjgILAUCbo2I/37IzNKIj+XvwvuB/WW3y18flSfT4gjr3RfZCbhs3JxvUBr1dIekL1EKhGr6+NWeemWbtyr69KmI+BFmR8mHhiwFDwAfkTQdSlcXDeMwygdV+u7jiZTGlP/XGtbbTilYXstGrr2kbNkeYHLZ9IvAvOx2ebtK9wKfGhi+WKXv+W2vsU9mgPcILAERsSUb9fEBSW1AL/BfqW3c94eB71H6YpTbB67yybPeiOhW6XuTnwZ+DjxWtvg7wC2S9lE6D3E98G1JvwAeH6KemymNLropOyz1CqWAMquZLx81O4LsipzTI+Izja7FrAg+NGRmljjvEZiZJc57BGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmifv/4gtsoMkmdNcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n", "import matplotlib.pyplot as plt\n", "\n", "data[\"Frequency\"]=data.Malfunction/data.Count\n", "data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression\n", "\n", "Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Wed, 15 Apr 2020 Deviance: 3.0144
Time: 15:20:47 Pearson chi2: 5.00
No. Iterations: 6 Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740
Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110
" ], "text/plain": [ "\n", "\"\"\"\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: Frequency No. Observations: 23\n", "Model: GLM Df Residuals: 21\n", "Model Family: Binomial Df Model: 1\n", "Link Function: logit Scale: 1.0000\n", "Method: IRLS Log-Likelihood: -3.9210\n", "Date: Wed, 15 Apr 2020 Deviance: 3.0144\n", "Time: 15:20:47 Pearson chi2: 5.00\n", "No. Iterations: 6 Covariance Type: nonrobust\n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import statsmodels.api as sm\n", "\n", "data[\"Success\"]=data.Count-data.Malfunction\n", "data[\"Intercept\"]=1\n", "\n", "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", " family=sm.families.Binomial(sm.families.links.logit)).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Wed, 15 Apr 2020 Deviance: 18.086
Time: 15:21:05 Pearson chi2: 30.0
No. Iterations: 6 Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err z P>|z| [0.025 0.975]
Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068
Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" ], "text/plain": [ "\n", "\"\"\"\n", " Generalized Linear Model Regression Results \n", "==============================================================================\n", "Dep. Variable: Frequency No. Observations: 23\n", "Model: GLM Df Residuals: 21\n", "Model Family: Binomial Df Model: 1\n", "Link Function: logit Scale: 1.0000\n", "Method: IRLS Log-Likelihood: -23.526\n", "Date: Wed, 15 Apr 2020 Deviance: 18.086\n", "Time: 15:21:05 Pearson chi2: 30.0\n", "No. Iterations: 6 Covariance Type: nonrobust\n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", " family=sm.families.Binomial(sm.families.links.logit),\n", " var_weights=data['Count']).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\n", "The Goodness of fit (Deviance) indicated for this model is $G^2=18.086$ with 21 degrees of freedom (Df Residuals).\n", "\n", "**I have therefore managed to fully replicate the results of the Dalal *et al.* article**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predicting failure probability\n", "The temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU5dn/8c81k8lGFiBAWMKmBhDZswBiLVgFtIobCoi4FMQ+VVtrpZU+Vq3VLg99fu5VKOBaRWoFqfURBMUFEQKCrLIjJOxLQkL25Pr9MQPGGMgkmcksud6vV16Zc+Y+51x3JvnOyZlz7iOqijHGmNDnCHQBxhhjfMMC3RhjwoQFujHGhAkLdGOMCRMW6MYYEyYs0I0xJkzUGugiMltEDonIhjM8LyLytIhsF5F1IjLA92UaY4ypjTd76C8BI8/y/OVAqudrMvB8w8syxhhTV7UGuqp+Ahw7S5OrgVfU7QuguYi081WBxhhjvBPhg3V0APZWmc72zNtfvaGITMa9F09MTExax44d67XByspKHI7wOPxvfQk+4dIPsL4Eq4b0ZevWrUdUtXVNz/ki0KWGeTWOJ6CqM4AZAOnp6bpq1ap6bXDp0qUMHTq0XssGG+tL8AmXfoD1JVg1pC8i8s2ZnvPF2102UHVXOwXY54P1GmOMqQNfBPoC4BbP2S6DgDxV/d7hFmOMMf5V6yEXEXkDGAq0EpFs4GHABaCqLwDvAVcA24FC4HZ/FWuMMebMag10VR1Xy/MK3OWziowxIaGsrIzs7GyKi4sbZXuJiYls3ry5Ubblb970JTo6mpSUFFwul9fr9cWHosaYJig7O5v4+Hi6dOmCSE3nRvhWfn4+8fHxft9OY6itL6rK0aNHyc7OpmvXrl6vNzzOATLGNLri4mKSkpIaJcybGhEhKSmpzv/9WKAbY+rNwtx/6vOztUA3xpgwYcfQjTEhy+l00rt379PT8+fPp0uXLoErKMAs0I0xISsmJoa1a9ee8fny8nIiIppOzNkhF2NMWHnppZe44YYbuOqqqxg+fDgA06ZNIyMjgz59+vDwww+fbvv444/TvXt3Lr30UsaNG8df//pXAIYOHcqpoUmOHDlyeq+/oqKCKVOmnF7X9OnTgW8v5R89ejQ9evRg/PjxuM/ohqysLC688EL69u1LZmYm+fn5jBgx4jtvREOGDGHdunUN7nvTeesyxvjN7/+9kU37Tvh0nT3bJ/DwVRectU1RURH9+vUDoGvXrsybNw+A5cuXs27dOlq2bMmiRYvYtm0bK1euRFUZNWoUn3zyCc2aNWPOnDmsWbOG8vJyBgwYQFpa2lm3N2vWLBITE8nKyqKkpIQhQ4acftNYs2YNGzdupH379gwZMoRly5aRmZnJmDFjePPNN8nIyODEiRPExMRwyy238NJLL/Hkk0+ydetWSkpK6NOnT4N/ZhboxpiQdaZDLpdddhktW7YEYNGiRSxatIj+/fsDUFBQwLZt28jPz+faa68lNjYWgFGjRtW6vUWLFrFu3TreeustAPLy8ti2bRuRkZFkZmaSkpICQL9+/di9ezeJiYm0a9eOjIwMABISEgC49tprGTJkCNOmTWP27NncdtttDftBeFigG2MarLY96cbWrFmz049VlalTp3LnnXd+p82TTz55xlMDIyIiqKysBPjOueCqyjPPPMOIESO+037p0qVERUWdnnY6nZSXl6OqNW4jNjaWyy67jHfeeYe5c+dS35Fnq7Nj6MaYsDZixAhmz55NQUEBADk5ORw6dIiLL76YefPmUVRURH5+Pv/+979PL9OlSxdWr14NcHpv/NS6nn/+ecrKygDYunUrJ0+ePOO2e/Towb59+8jKygLcV4iWl5cDMGnSJH7+85+TkZFx+r+JhrI9dGNMWBs+fDibN29m8ODBAMTFxfHaa68xYMAAxowZQ79+/ejcuTM/+MEPTi9z//33c+ONN/Lqq69yySWXnJ4/adIkdu/ezYABA1BVWrduzfz588+47cjISN58803uueceioqKiImJYfHixQCkpaWRkJDA7bf7cDxDVQ3IV1pamtbXRx99VO9lg431JfiESz9U/duXTZs2+W3dNTlx4oRf1//www/rtGnT/LqNU06cOKE5OTmampqqFRUVZ2xX088YWKVnyFU75GKMMY3s9ddfZ+DAgTz++OM+va2eHXIxxhjgkUceabRt3XTTTd/7kNYXbA/dGFNvqjXePtj4QH1+thboxph6iY6O5ujRoxbqfqCe8dCjo6PrtJwdcjHG1EtKSgrZ2dkcPny4UbZXXFxc54ALVt705dQdi+rCAt0YUy8ul6tOd9NpqKVLl56+2jPU+asvdsjFGGPChAW6McaECQt0Y4wJExboxhgTJizQjTEmTFigG2NMmLBAN8aYMGGBbowxYcIC3RhjwoQFujHGhImQC/RDJ4r5OLvMBgQyxphqQi7QX1uxhxc3lDLp5VUczi8JdDnGGBM0Qi7Q7/1RKjf1iOTT7UcY+eQnLN50MNAlGWNMUAi5QHc4hOFdXPznnotITohm0iurePidDRSXVQS6NGOMCaiQC/RTUpPjmXfXhUy8qCsvL/+Ga55bxvZDBYEuyxhjAiZkAx0gKsLJ767syYu3Z3Aov4RRz37GO2tzAl2WMcYEhFeBLiIjRWSLiGwXkQdqeL6TiHwkImtEZJ2IXOH7Us9sWPc2/OfnF3FB+wR+MWctv523npJyOwRjjGlaag10EXECzwGXAz2BcSLSs1qzB4G5qtofGAv8zdeF1qZdYgxv3DGIn/7wXF5fsYcbX1hOTm5RY5dhjDEB480eeiawXVV3qmopMAe4ulobBRI8jxOBfb4r0XsRTgcPXN6D6RPS2Hn4JFc+/Smfbz8SiFKMMabRSW0X6IjIaGCkqk7yTE8ABqrq3VXatAMWAS2AZsClqrq6hnVNBiYDJCcnp82ZM6deRRcUFBAXF3fWNgdOVvL0mmIOnFTGdI9keOcIRKRe2/Mnb/oSKsKlL+HSD7C+BKuG9GXYsGGrVTW9xidV9axfwA3AzCrTE4BnqrW5D/iV5/FgYBPgONt609LStL4++ugjr9rlF5fpHS9naeffvKv3vblWi8vK671Nf/G2L6EgXPoSLv1Qtb4Eq4b0BVilZ8hVbw65ZAMdq0yn8P1DKhOBuZ43iOVANNDKi3X7VVxUBC/cnMa9l6byry+zGf/3FRwpsKtLjTHhyZtAzwJSRaSriETi/tBzQbU2e4AfAYjI+bgD/bAvC60vh0O499JuPHfTADbsy+PqZ5ex5UB+oMsyxhifqzXQVbUcuBtYCGzGfTbLRhF5VERGeZr9CrhDRL4C3gBu8/xrEDR+3Kcdc+8cTFlFJaOf/5xPtwXF+40xxviMV+ehq+p7qtpNVc9V1cc98x5S1QWex5tUdYiq9lXVfqq6yJ9F11eflObMv2sIHVrEcPuLWbyZtSfQJRljjM+E9JWi9dG+eQz//OlgLjyvFb/513qeXLzVhuI1xoSFJhfoAPHRLmbdms7otBSeXLyNqW+vp7yiMtBlGWNMg0QEuoBAcTkdTBvdh/aJ0Tz94XaOFJTw7E0DiHY5A12aMcbUS5PcQz9FRLhveHf+cPUFLPn6ELfMWkleUVmgyzLGmHpp0oF+yoTBXXh6bH/W7D3OmOnL7U5IxpiQZIHucVXf9sy6NYPdR08yZroN7GWMCT0W6FVc3K01r04cyOH8Em58YTm7j5wMdEnGGOM1C/RqMrq05I3Jgygqq+DG6cvtLkjGmJBhgV6DXh0SmTN5EJUKY2cs5+sDJwJdkjHG1MoC/Qy6Jccz985BRDgcjJ3xBRv35QW6JGOMOSsL9LM4p3Ucb945iFiXk/EzV7Ahx0LdGBO8LNBr0TmpGXMmD7ZQN8YEPQt0L3RKimXO5ME0i3Ry86wVbN5vx9SNMcHHAt1LnZJieWPyIKIjnNw8cwXbDtqY6saY4GKBXgedk5rx+h0DcTqEcX9fwc7DdkqjMSZ4WKDX0Tmt43j9joGoKuNnrmDvscJAl2SMMYAFer2c1yaeVycOpLC0gvEzV3AgrzjQJRljjAV6ffVsn8DLP8nk2MlSbp61gmMnSwNdkjGmibNAb4B+HZsz89Z09h4r5NbZK8kvtqF3jTGBY4HeQIPOSeL5mwewef8JJr28iuKyikCXZIxpoizQfeCSHsn87419Wbn7GHe/vsZuZ2eMCQgLdB+5ul8HHrnqAhZvPsjUt9fbjaeNMY2uyd5T1B9uvbALR0+W8vSSbSTFRfHA5T0CXZIxpgmxQPexX16ayrGTJbzw8Q7axEfxk4u6BrokY0wTYYHuYyLC70f14kh+KX/4zyZax0dxVd/2gS7LGNME2DF0P3A6hCfH9iOjc0vum7uWz7cfCXRJxpgmwALdT6JdTv5+SzpdWzXjzldX212PjDF+Z4HuR4mxLl68PZPYKCe3v5jF/ryiQJdkjAljFuh+1qF5DC/elkl+cTm3v5hlV5MaY/zGAr0R9GyfwPM3D2D7oQJ+9o8vKbMLj4wxfmCB3kh+kNqaP17bm0+3HeHBeRvswiNjjM/ZaYuN6MaMjuw9XsgzH26nU1Isdw07L9AlGWPCiAV6I7vvsm7sOVbItIVb6JwUS1ygCzLGhA075NLIRIS/XN+H9M4tuG/uV2w/bqMzGmN8w6tAF5GRIrJFRLaLyANnaHOjiGwSkY0i8rpvywwv0S4nM25Jp11iNE+tKbbb2BljfKLWQBcRJ/AccDnQExgnIj2rtUkFpgJDVPUC4F4/1BpWWjaLZPZtGVRUwsSXszhhpzMaYxrImz30TGC7qu5U1VJgDnB1tTZ3AM+p6nEAVT3k2zLD07mt47i7fzQ7D5+0cdSNMQ0mtZ0+JyKjgZGqOskzPQEYqKp3V2kzH9gKDAGcwCOq+n4N65oMTAZITk5OmzNnTr2KLigoIC4uPD5OLCgoYPXxKF7cWMqlnSK4uWdUoEuqt3B5XcKlH2B9CVYN6cuwYcNWq2p6Tc95c5aL1DCv+rtABJAKDAVSgE9FpJeq5n5nIdUZwAyA9PR0HTp0qBeb/76lS5dS32WDzdKlS3n4yqE4393EzM92MXRAD24e1DnQZdVLuLwu4dIPsL4EK3/1xZtDLtlAxyrTKcC+Gtq8o6plqroL2II74I2Xpl5xPsO6t+bhBRttdEZjTL14E+hZQKqIdBWRSGAssKBam/nAMAARaQV0A3b6stBw53QIT4/rz7mtm/HT11az68jJQJdkjAkxtQa6qpYDdwMLgc3AXFXdKCKPisgoT7OFwFER2QR8BExR1aP+KjpcxUe7mHlLBk6HMPHlLPKK7MwXY4z3vDoPXVXfU9Vuqnquqj7umfeQqi7wPFZVvU9Ve6pqb1Wt36edhk5JsTx/cxp7jhZyzxt25osxxnt2pWgQGnROEn+4phefbD3Mn/7v60CXY4wJETaWS5Aal9mJLQfymfXZLnq0jeeG9I61L2SMadJsDz2IPfjj8xlyXhL/PW8Dq785HuhyjDFBzgI9iEU4HTx30wDaNY/mzldX2y3sjDFnZYEe5JrHRjLzlnSKyyqY/MpqistsdEZjTM0s0ENAanI8T47px4Z9efzmX+vsbkfGmBpZoIeIS3smc//w7ryzdh/TP7Frtowx32eBHkJ+NvRcftynHX95/2uWbrEBLY0x32WBHkJEhGmj+9CjbQL3vLGGnYcLAl2SMSaIWKCHmNjICGZMSMPldHDHK6vsxhjGmNMs0ENQx5axPHfTAHYfLeS+N9dSWWkfkhpjLNBD1uBzk3joyp4s3nyIJxZvDXQ5xpggYJf+h7BbBndm4748nvlwO+e3S+CK3u0CXZIxJoBsDz2EiQh/uKYX/Ts15/5/fsXXB04EuiRjTABZoIe4qAgnL9ycRlxUBHe8sorcwtJAl2SMCRAL9DCQnBDNCxPSOJhXwt2v2xjqxjRVFuhhYkCnFjx2TS8+236EP9sY6sY0SfahaBi5MaMjG/blMfOzXVzQIYFr+6cEuiRjTCOyPfQw87srezKwa0se+Nd61mXnBrocY0wjskAPMy6ng7+NH0CruCjufHU1h/NLAl2SMaaRWKCHoaS4KKZPSON4YSk/+8dqSsvtQ1JjmgIL9DDVq0Mif7m+D1m7j/PIvzcGuhxjTCOwD0XD2NX9OrBp/wmmf7yTC9onMH5g50CXZIzxI9tDD3O/HtGDH3ZrzcPvbGTlrmOBLscY40cW6GHO6RCeHtefji1j+a/XVpOTazeaNiZcWaA3AYkxLv5+Szql5ZVMfmUVRaV2o2ljwpEFehNxXps4nhrXj037TzDlra/sRtPGhCEL9Cbkkh7JTBnRnXfX7edvS3cEuhxjjI9ZoDcx//XDc7mqb3v+umgLSzYfDHQ5xhgfskBvYkSE/7m+Dxe0T+AXc9ay7WB+oEsyxviIBXoTFBPpZMaEdKJdTibZGOrGhA0L9CaqffMYpk9IY39uMT/7x5eU2RjqxoQ8C/QmLK1zC/54XW8+33GUP7y7KdDlGGMayC79b+JGp6Ww9WA+Mz7ZSWpyPBMG2fAAxoQq20M3/GZkD4Z1b80jCzby+fYjgS7HGFNPXgW6iIwUkS0isl1EHjhLu9EioiKS7rsSjb+dGh7gnFbN+K9/fMmuIycDXZIxph5qDXQRcQLPAZcDPYFxItKzhnbxwM+BFb4u0vhffLSLWbdm4HQIE1/KIq+wLNAlGWPqyJs99Exgu6ruVNVSYA5wdQ3t/gD8D1Dsw/pMI+qUFMsLN6ex93ghd71uZ74YE2qktjE9RGQ0MFJVJ3mmJwADVfXuKm36Aw+q6vUishS4X1VX1bCuycBkgOTk5LQ5c+bUq+iCggLi4uLqtWywCca+fJpdxqwNpQztGMGtPSMREa+WC8a+1Ee49AOsL8GqIX0ZNmzYalWt8bC2N2e51PTXfPpdQEQcwBPAbbWtSFVnADMA0tPTdejQoV5s/vuWLl1KfZcNNsHYl6FA5Ptf8/zSHfygbzcmXtTVq+WCsS/1ES79AOtLsPJXX7w55JINdKwynQLsqzIdD/QClorIbmAQsMA+GA1tU4Z3Z+QFbXnsP5tYvMnGfDEmFHgT6FlAqoh0FZFIYCyw4NSTqpqnqq1UtYuqdgG+AEbVdMjFhA6HQ3hiTD96tU/k53PWsCEnL9AlGWNqUWugq2o5cDewENgMzFXVjSLyqIiM8neBJnBiIp3MujWd5jEuJr6cxf48u9uRMcHMq/PQVfU9Ve2mqueq6uOeeQ+p6oIa2g61vfPw0SYhmtm3Z3CypILbX8wiv9hOZzQmWNmVoqZWPdom8Nz4AWw7VMBdr6+x0xmNCVIW6MYrP+zWmsev6cUnWw/z0Dsb7BZ2xgQhG5zLeG1sZif2Hi/kuY92kNIilruGnRfokowxVVigmzr51WXdyT5exLSFW2iXGM11A1ICXZIxxsMC3dSJwyFMG92Xw/kl/PqtdbSJj+ai1FaBLssYgx1DN/UQGeHghQlpnNcmjp++tpqN++p/jvr8NTkM+fOHdH3gPwz584fMX5Pjw0qNv9nrF1ws0E29JES7eOn2TBKiI7jtxSz2Hius8zrmr8lh6tvrycktQoGc3CKmvr3eQiFE2OsXfCzQTb21TYzmlYmZlJZXcsvslZworduZL9MWbqGorOI784rKKpi2cIsvyzR+Yq9f8LFANw1yXpt4Zt+Wzv68Ip5YVUxBSbnXy+7LrfnK0zPNN8HFXr/gY4FuGiytc0v+Nn4A3+RXMvmVVZSUV9S+ENC+eUyd5pvgYq9f8LFANz5xSY9kJvaK5PMdR7l3zloqKms//DJlRHdiXM7vzItxOZkyoru/yjQ+ZK9f8LFANz4zpIOL313Zk//bcICpb6+r9WrSa/p34E/X9aZD8xgE6NA8hj9d15tr+ndonIJNg9jrF3zsPHTjUxMv6kpeURlPL9lGQrSL//7x+We949E1/TtYAIQwe/2CiwW68blfXprKiaIyZn62i/hoF7+4NDXQJRnTJFigG58TER66sicFJeU8sXgrsZFO7rj4nECXZUzYs0A3fuFwCH+5vg9FZRU8/t5mol0OJgzuEuiyjAlrFujGb5wO4Ykb+1FcWsHv3tlIZISDMRmdAl2WMWHLznIxfhUZ4eC58QO4uFtrHnh7Pf9anR3okowJWxboxu+iXU5mTEjjwnOTmPLWV7yz1sb6MMYfLNBNo4h2OZl5SwaZXVvyyzfX2gBOxviBBbppNDGRTmbflsHArkncN3ct89bY4RdjfMkC3TSq2MgIZt+WwaBzkrhv7lfMXbU30CUZEzYs0E2ji4l0MuvWDC46rxW/fmsdr33xTaBLMiYsWKCbgIiJdPL3W9L5UY82PDh/A7M+2xXokowJeRboJmCiXU6evzmNy3u15Q/vbuKpxdtqHdDLGHNmFugmoCIjHDwzrj/XD0jhicVbefw/my3Ujaknu1LUBFyE08G00X2Ij45g5me7yCsq40/X9SbCafsbxtSFBboJCg6H8PBVPUmMcfHUkm0cLyzj2Zv6E13tBgrGmDOzXSATNESEX17Wjd+PuoAlXx9kwqwV5BaWBrosY0KGBboJOrde2IWnx/bnq715jH5hOdnHCwNdkjEhwQLdBKWr+rbn5Z9kcvBEMdf97XM25OQFuiRjgp4Fuglag89N4q2fXojTIdw4fTkffn0w0CUZE9Qs0E1Q6942nvl3DaFrq2ZMenkVL3++O9AlGRO0LNBN0EtOiGbunYO5pEcbHl6wkYfe2UB5RWWgyzIm6HgV6CIyUkS2iMh2EXmghufvE5FNIrJORJaISGffl2qasmZREUyfkM7ki8/hleXfcNuLWeQVlgW6LGOCSq2BLiJO4DngcqAnME5EelZrtgZIV9U+wFvA//i6UGOcDuG3V5zPtNF9WLHrKKOe+4ytB/MDXZYxQcObPfRMYLuq7lTVUmAOcHXVBqr6kaqeOrfsCyDFt2Ua860b0jsyZ/IgCksruPa5Zby/YX+gSzImKEht42aIyGhgpKpO8kxPAAaq6t1naP8scEBVH6vhucnAZIDk5OS0OXPm1KvogoIC4uLi6rVssLG+1N/x4kqeWVPCzrxKrujq4vpUF06HNHi99poEJ+uL27Bhw1aranqNT6rqWb+AG4CZVaYnAM+coe3NuPfQo2pbb1pamtbXRx99VO9lg431pWGKy8p16tvrtPNv3tWx05fr4fziBq/TXpPgZH1xA1bpGXLVm0Mu2UDHKtMpwL7qjUTkUuC/gVGqWuLtu40xDREV4eSP1/bmrzf05cs9x7niqU9ZvuNooMsyJiC8CfQsIFVEuopIJDAWWFC1gYj0B6bjDvNDvi/TmLMbnZbC/LuGEBcVwfiZX/D0km1UVNowvKZpqTXQVbUcuBtYCGwG5qrqRhF5VERGeZpNA+KAf4rIWhFZcIbVGeM357dLYME9FzGqb3v+3wdbGT/zC/bnFQW6LGMajVfD56rqe8B71eY9VOXxpT6uy5h6WbzpICt3HQNgxc5j/Oh/P2ZsRkcWbjzIvtwi2jePYcqI7lzTv4PPtz1/TQ7TFm7x+3a88eD89byxYi/39ipj4tT3GDewI49d0zsgtZjGY+Ohm7Axf00OU99eT1FZBQAKFJVWMHvZ7tNtcnKLmPr2egCfhm31bftrO954cP56Xvtiz+npCtXT0xbq4c0u/TdhY9rCLacD9ZSajqIXlVUwbeEWv2/bH9vxxhsr9tZpvgkfFugmbOzL9f54eU4d2jZk23WpyVcqznBtyZnmm/BhgW7CRvvmMV63dTqEz7Yd8fu261KTrzil5ourzjTfhA8LdBM2pozoTky1e5C6HILL+d0gi3Q6aNkskptnreD+f37FsZMNv81dTduOcTmZMqJ7g9ddV+MGdqzTfBM+7ENREzZOffhY/UyTmuaN7NWWp5dsY8YnO1my+SC/veJ8RqelIPXciz3TtgNxlsupDz5PHTN3ithZLk2EBboJK9f071BjiNY079cjezCqX3senLeBKW+tY+6qvfx+VC+fbzsQHrumN49d05ulS5eyY/zQQJdjGokdcjFNWo+2Ccy9czB/ub43Ow6f5MpnPuXVTSXkFjb8MIwxjc0C3TR5DocwJqMTH/1qKDcP6syHe8r54bSlvLhsF2V2ZyQTQizQjfFIjHXx6NW9eHRIDL06JPD7f29ixBOf8P6GA6dGEzUmqFmgG1NNx3gHr00cyMxb0nE4hJ++tprRLyxnxU4bxdEENwt0Y2ogIlzaM5n3f/ED/nRdb/YeK2TMjC+4ZfZK1mXnBro8Y2pkgW7MWUQ4HYzL7MTHU4bx2yt6sC47l1HPLmPiS1kW7CboWKAb44WYSCeTLz6XT389jPuHd2PVN8cZ9ewybntxJVm7jwW6PGMAC3Rj6iQ+2sXdl6Ty2W+GMWVEd9Zn53HDC8u58YXlLN50kEq7qYYJIAt0Y+ohPtrFXcPO47PfXMJDV/YkJ7eISa+s4rInPub1FXsoKq2ofSXG+JgFujENEBPp5CcXdWXplKE8NbYf0S4nv523nsF/XsKf/+9r9h4rDHSJpgmxS/+N8QGX08HV/Towqm97snYf58Vlu5jxyQ6mf7KDS7q3YfygTvywWxucDhvx0PiPBboxPiQiZHZtSWbXluzLLeKNlXt4Y+Velry0ivaJ0dyQ3pHRaSl0bBkb6FJNGLJAN8ZP2jeP4VfDu3PPJaks2XyQ11fu4ekPt/HUkm0MPieJ69NSGNmrLXFR9mdofMN+k4zxs8gIB5f3bsflvduRfbyQt7/M4a3V2dz/z694cP56LuvZllF923Nxt1ZERThrX6ExZ2CBbkwjSmkRy89/lMo9l5zHl3uOM29NDu+u28+/v9pHfHQEw3u25fJebbkotRXRLgt3UzcW6MYEgIiQ1rklaZ1b8vBVF7Bs+xH+/dV+Pth0gH99mU1cVAQ/7N6a4T2TGdq9DYkxrkCXbEKABboxAeZyOhjavQ1Du7ehtLw3n+84wvsbDrB48yH+s24/ToeQ0aUFl/Rwt0ltE1fvOyuZ8GaBbkwQiYz4NtwrK5U1e3P58OuDLNl8iD++9zV/fO9r2iVG84PUVlyU2poLz02iVVxUoMs2QcIC3Zgg5XAIaZ1bkNa5BcegSTwAAA0gSURBVFNG9GBfbhGfbD3Mx1sP8/6GA8xdlQ1Aj7bxDDoniUHnJJHRpQVJFvBNlgW6MSGiffMYxmZ2YmxmJyoqlfU5eSzbfoTlO44yJ2sPL32+G4Dz2sSR7nkjqDhZiaraIZomwgLdmBDkdAj9OjanX8fm3DXsPErKK1ifncfK3cfI2nWM99bvZ07WXgD+vPoD+qY0p2/H5vRNSaR3SiJt4qMD3APjDxboxoSBqAgn6V1akt6lJQyFykplx+ECXl/0BUWxyazdm8uzH27j1GCQbeKj6NUhkZ7tEujZPoHz2yXQqWWsDU0Q4izQjQlDDoeQmhzPDzu6GDq0DwCFpeVs3HeCddl5bNyXx4acPD7eepgKT8pHuxyktomnW3I83ZLjSE2O47zW8XRoEWNBHyIs0I1pImIjI8jo0pKMLi1Pzysuq2D7oQI27T/BlgP5bDmQz6fbDvOvL7NPt4mKcNC1VTO6tmpGl1bN6JrUjM5JsXRp1YzWcVE4LOyDhgW6MU1YtMtJrw6J9OqQ+J35uYWlbD9UwI7DBWw/VMCuIyfZcjCfDzYdpLzKTTyiIhx0bBlLSosYOraIpUOLGDo0jzn9vVVclO3dNyILdGPM9zSPjfz2mHwV5RWV7MstZtfRk+w5VsjeY4V8c/Qk2ceLWLMnl7yisu+0j3AIyQnRtE2Mpm1CNMkJ0SQnRNEmIYo28dG0iY+iVVwUzWNddiaOD1igG2O8FuF00Ckplk5JNQ//m19cRk5uETnHi9iXV8z+3CIO5BVz4EQxm/ef4KMthyis4W5OLqeQ1CyKpLhIkuKiSGoWScsqXy1iI/nmWAXtDuTTItZFQozLxrqpgQW6McZn4qNd9GjrokfbhDO2KSgp5+CJYg6dKOFQfjFHCko5UlDC4fwSjp10P955uIBjJ0u/F/5/WvnJ6cfRLgeJMS4Sol3u7zEuEqIjiI92Ee/5HhcdQXxUBM2iIojzfDWLchIXFUFsVASxLmdYfQbgVaCLyEjgKcAJzFTVP1d7Pgp4BUgDjgJjVHW3b0s1JnzNX5PDtIVb2JdbRPvmMUwZ0Z1/rtrDsh3HTrcZcm5Lbkjv9L12wPfmrfrmGG+s2Mu9vcqYOPU9xg3syGPX9PZquzWt75r+Hbyu+9S2K1Rxinxv23FREcS1jmN9dl6t237kqlR+0K0Vx06W8vHyVXRKPZ/jhWV8seMoH289zMETJZ7DPLEUlVWw/VA5J4rLyC8uP332Tm1iXE5iI53ERJ76HkGMy0FsZAQxLidRLgcxLifRLifRLgfREd8+jopwPx8V4Xkc4SDK5SDS6SQywvHtl/Pb7y6noOqfm4nXGugi4gSeAy4DsoEsEVmgqpuqNJsIHFfV80RkLPAXYIw/CjYm3Mxfk8PUt9dTVObeG83JLeLeN9d+r92yHce+E/A5uUVMeesrUCjzhFdObhH3vbmWyirLVajy2hd7AL4TrDVtd8o/vwKBsopv1zf17fUA3wv1mpb39bYfXrCRP13Xm2v6d+BwkpOhfdozf00OH3596PSyxWWVZB8vOt0OQFUpLqskv6SMguJyCko8X8XlFJZWcLK0nJMl5ZwsqeBkSTmFZRUUlVZQWFpOUVklRaXlHM4vocgzv7isgqIy93cv3yfOakLPSIY1fDXf480eeiawXVV3AojIHOBqoGqgXw084nn8FvCsiIj6623ImDAybeGW0+FUV6fCr6rKGtoBvLFi73dCtabtltWQVkVlFUxbuOV7gV7T8o2x7ZqWrd5ORIjx7HW3iT9DUfWgqpRVKMXlFZSUVVJcVkFpRSUlZZWUlFdQWl5JSXklpeWV7vnlFZSVKyUV7nllFZWUlVcSX7DHd0VVIbVlroiMBkaq6iTP9ARgoKreXaXNBk+bbM/0Dk+bI9XWNRmY7JnsDmypZ92tgCO1tgoN1pfg06j9iGx7Xpq/1l1RmIcz9ttTEksPbF9d3+1WXbahy9dz2VbAkbMtW73GINaQ37HOqtq6pie82UOv6ROD6u8C3rRBVWcAM7zY5tkLElmlqukNXU8wsL4En3DpB7j7Up53KGz6Ek6viz/64vCiTTbQscp0CrDvTG1EJAJIBI5hjDGm0XgT6FlAqoh0FZFIYCywoFqbBcCtnsejgQ/t+LkxxjSuWg+5qGq5iNwNLMR92uJsVd0oIo8Cq1R1ATALeFVEtuPeMx/rz6LxwWGbIGJ9CT7h0g+wvgQrv/Sl1g9FjTHGhAZvDrkYY4wJARboxhgTJoI+0EUkWkRWishXIrJRRH7vmd9VRFaIyDYRedPzgW3QExGniKwRkXc906Haj90isl5E1orIKs+8liLygacvH4hIi0DX6Q0RaS4ib4nI1yKyWUQGh2JfRKS75/U49XVCRO4N0b780vP3vkFE3vDkQKj+rfzC04+NInKvZ55fXpOgD3SgBLhEVfsC/YCRIjII9/ACT6hqKnAc9/ADoeAXwOYq06HaD4Bhqtqvyvm0DwBLPH1Z4pkOBU8B76tqD6Av7tcn5Pqiqls8r0c/3OMqFQLzCLG+iEgH4OdAuqr2wn0yxqkhRULqb0VEegF34L7ivi9wpYik4q/XRFVD5guIBb4EBuK+yirCM38wsDDQ9XlRf4rnxbsEeBf3BVkh1w9PrbuBVtXmbQHaeR63A7YEuk4v+pEA7MJzgkAo96Va/cOBZaHYF6ADsBdoiftMvHeBEaH4twLcgHtAw1PTvwN+7a/XJBT20E8dplgLHAI+AHYAuapa7mmSjfuXINg9ifvFPDXkRRKh2Q9wXwm8SERWe4Z0AEhW1f0Anu9tAlad984BDgMveg6FzRSRZoRmX6oaC7zheRxSfVHVHOCvwB5gP5AHrCY0/1Y2ABeLSJKIxAJX4L4I0y+vSUgEuqpWqPvfyBTc/7qcX1Ozxq2qbkTkSuCQqlYda8KrIROC1BBVHQBcDtwlIhcHuqB6igAGAM+ran/gJEF+SKI2nmPLo4B/BrqW+vAcT74a6Aq0B5rh/j2rLuj/VlR1M+5DRR8A7wNfAeVnXagBQiLQT1HVXGApMAho7hlmAGoejiDYDAFGichuYA7uwy5PEnr9AEBV93m+H8J9nDYTOCgi7QA83w8FrkKvZQPZqrrCM/0W7oAPxb6ccjnwpaoe9EyHWl8uBXap6mFVLQPeBi4kdP9WZqnqAFW9GPeFl9vw02sS9IEuIq1FpLnncQzuF3sz8BHuYQbAPezAO4Gp0DuqOlVVU1S1C+5/hz9U1fGEWD8ARKSZiMSfeoz7eO0GvjsEREj0RVUPAHtFpLtn1o9wDw0dcn2pYhzfHm6B0OvLHmCQiMSKiPDtaxJyfysAItLG870TcB3u18Yvr0nQXykqIn2Al3F/0u0A5qrqoyJyDu493ZbAGuBmVS0JXKXeE5GhwP2qemUo9sNT8zzPZATwuqo+LiJJwFygE+4/yhtUNegHaRORfsBMIBLYCdyO53eN0OtLLO4PFM9R1TzPvJB7XTynJ4/BfXhiDTAJ9zHzkPpbARCRT3F/XlYG3KeqS/z1mgR9oBtjjPFO0B9yMcYY4x0LdGOMCRMW6MYYEyYs0I0xJkxYoBtjTJjw5ibRxjQqzyldSzyTbYEK3JfnA2SqamlACjsLEfkJ8J7nvHZjAsJOWzRBTUQeAQpU9a9BUItTVSvO8NxnwN2qurYO64uoMjaJMQ1mh1xMSBGRW8U9Pv5aEfmbiDhEJEJEckVkmoh8KSILRWSgiHwsIjtF5ArPspNEZJ7n+S0i8qCX631MRFYCmSLyexHJ8oxv/YK4jcE9tPObnuUjRSS7yhXOg0RksefxYyIyXUQ+wD0gWISI/D/PtteJyKTG/6macGGBbkKGZ2zpa4ELPYO1RfDtDckTgUWeAcNKgUdwXzJ+A/BoldVkepYZANwkIv28WO+XqpqpqsuBp1Q1A+jteW6kqr4JrAXGqHs88toOCfUHrlLVCcBk3IO2ZQIZuAc661Sfn48xdgzdhJJLcYfeKvcQH8TgvswdoEhVP/A8Xg/kqWq5iKwHulRZx0JVPQ4gIvOBi3D/HZxpvaV8O8wBwI9EZAoQDbTCPazr/9WxH++oarHn8XDgfBGp+gaSivtycGPqxALdhBIBZqvq774z0z0CX9W94krcd7o69bjq73n1D420lvUWqeeDJs84Kc8CA1Q1R0Qewx3sNSnn2/+Aq7c5Wa1PP1PVJRjTQHbIxYSSxcCNItIK3GfD1OPwxHBx30M0FveY28vqsN4Y3G8QRzyjTV5f5bl8IL7K9G7ct4GjWrvqFgI/OzUsrLjvCxpTxz4ZA9geugkhqrreMwrfYhFx4B697qfUbVzsz4DXgXOBV0+dleLNelX1qIi8jHuo4G+AFVWefhGYKSJFuI/TPwL8XUQOACvPUs903CPurfUc7jmE+43GmDqz0xZNk+E5g6SXqt4b6FqM8Qc75GKMMWHC9tCNMSZM2B66McaECQt0Y4wJExboxhgTJizQjTEmTFigG2NMmPj/K5iDbLoKp3YAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n", "data_pred['Frequency'] = logmodel.predict(data_pred)\n", "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The predictiction is not correct. Let's see what went wrong" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhV5bn38e+dnZkwDxEIGFREkTkQRNSCE9Qqjgg41KpI+1Zrraee1reD1uprW3p6HOqAA9rqUbQeRVQqqAW1KAgIMsqMGkCZh0Dm3O8fewcDJmRn3AO/z3XtK3sNe637yUp+WXn22s8yd0dERGJfQqQLEBGRhqFAFxGJEwp0EZE4oUAXEYkTCnQRkTihQBcRiRM1BrqZTTazrWa2rJrlV5nZktDjQzPr2/BliohITcI5Q38GGHmE5RuA77h7H+D3wOMNUJeIiNRSYk0ruPv7ZpZ9hOUfVpqcC2TVvywREamtGgO9lm4A/lndQjObAEwASEtLy+nSpUuddlJeXk5CQnx0/6st0Sle2hIv7QC1pcLq1au3u3v7Khe6e40PIBtYVsM6w4GVQNtwtpmTk+N1NWvWrDq/NtqoLdEpXtoSL+1wV1sqAAu8mlxtkDN0M+sDPAl81913NMQ2RUSkdur9/4uZdQVeAa5x99X1L0lEROqixjN0M3sBGAa0M7M84E4gCcDdHwN+C7QFHjEzgFJ3H9hYBYuISNXCucplXA3LxwPjG6wiEYkJJSUl5OXlUVhY2CT7a9myJStXrmySfTW2cNqSmppKVlYWSUlJYW+3oa9yEZGjRF5eHs2bNyc7O5vQf+eNat++fTRv3rzR99MUamqLu7Njxw7y8vLo1q1b2NuNj2uARKTJFRYW0rZt2yYJ86ONmdG2bdta//ejQBeROlOYN566fG8V6CIicUJ96CISswKBAL179z44PXXqVLKzsyNXUIQp0EUkZqWlpbF48eJql5eWlpKYePTEnLpcRCSuPPPMM4wePZoLL7yQ8847D4CJEycyaNAg+vTpw5133nlw3XvvvZcePXpwzjnnMG7cOP785z8DMGzYMBYsWADA9u3bD571l5WVcfvttx/c1qRJkwCYPXs2w4YN4/LLL+ekk07iqquuqhgShfnz53PaaafRt29fcnNz2bdvHyNGjDjkD9HQoUNZsmRJvdt+9PzpEpFG87vXl7Ni894G3WbPTi2488JTjrhOQUEB/fr1A6Bbt268+uqrAHz00UcsWbKENm3aMHPmTNasWcPHH3+MuzNq1Cjef/99mjVrxpQpU1i0aBGlpaUMGDCAnJycI+7vqaeeomXLlsyfP5+ioiKGDh168I/GokWLWL58OZ06dWLo0KHMmTOH3NxcxowZw4svvsigQYPYu3cvaWlpfP/73+eZZ57h/vvvZ/Xq1RQVFdGnT596f88U6CISs6rrcjn33HNp06YNADNnzmTmzJn0798fgPz8fNasWcO+ffu45JJLSE9PB2DUqFE17m/mzJksWbKEl19+GYA9e/awZs0akpOTyc3NJSsrOHp4v3792LhxIy1btqRjx44MGjQIgBYtWgBwySWXMHToUCZOnMjkyZP5wQ9+UL9vRIgCXUTqraYz6abWrFmzg8/dnTvuuIMf/vCHh6xz//33V3tpYGJiIuXl5QCHXAvu7jz00EOMGDHikPVnz55NSkrKwelAIEBpaSnuXuU+0tPTOffcc3nttdd46aWXDnbv1Jf60EUkro0YMYLJkyeTn58PwKZNm9i6dStnnnkmr776KgUFBezbt4/XX3/94Guys7NZuHAhwMGz8YptPfroo5SUlACwevVq9u/fX+2+TzrpJDZv3sz8+fOB4CdES0tLARg/fjy33HILgwYNOvjfRH3pDF1E4tp5553HypUrGTJkCAAZGRk899xzDBgwgDFjxtCvXz+OPfZYzjjjjIOv+fnPf84VV1zBs88+y1lnnXVw/vjx49m4cSMDBgzA3Wnfvj1Tp06tdt/Jycm8+OKL/OQnP6GgoIC0tDTeeecdAHJycmjRogXXXXddwzW2uoHSG/uhG1wEqS3RKV7a0pjtWLFiRaNtuyp79+5t1O3feeedPnHixEbdR4W9e/f6pk2bvHv37l5WVlbtelV9jznCDS7U5SIi0sSef/55Bg8ezL333tugt9VTl4uICHDXXXc12b6uvPLKb71J2xB0hi4ideahD89Iw6vL91aBLiJ1kpqayo4dOxTqjcBD46GnpqbW6nXqchGROsnKyiIvL49t27Y1yf4KCwtrHXDRKpy2VNyxqDYU6CJSJ0lJSbW6m059zZ49++CnPWNdY7VFXS4iInFCgS4iEicU6CIicUKBLiISJxToIiJxQoEuIhInFOgiInFCgS4iEidqDHQzm2xmW81sWTXLzcweNLO1ZrbEzAY0fJkiIlKTcM7QnwFGHmH5d4HuoccE4NH6lyUiIrVVY6C7+/vAziOschHw99DY63OBVmbWsaEKFBGR8DREH3pn4MtK03mheSIi0oQsnKEvzSwbeMPde1Wx7E3gPnf/d2j6XeA/3X1hFetOINgtQ2ZmZs6UKVPqVHR+fj4ZGRl1em20UVuiU7y0JV7aAWpLheHDhy9094FVLqzu3nSVH0A2sKyaZZOAcZWmVwEda9qm7ikapLZEp3hpS7y0w11tqUAj31N0GvD90NUupwJ73H1LA2xXRERqocbx0M3sBWAY0M7M8oA7gSQAd38MmA6cD6wFDgDXNVaxIiJSvRoD3d3H1bDcgZsarCIREakTfVJURCROKNBFROKEAl1EJE4o0EVE4oQCXUQkTijQRUTihAJdRCROKNBFROKEAl1EJE4o0EVE4oQCXUQkTijQRUTihAJdRCROKNBFROKEAl1EJE4o0EVE4oQCXUQkTijQRUTihAJdRCROKNBFROKEAl1EJE4o0EVE4oQCXUQkTsRcoJeXO5/vLYt0GSIiUSfmAv2VRZu488NCbntpMVv3Fka6HBGRqBFzgT6y1zF8r1sSb3y6heF/ns1j762juLQ80mWJiERczAV6Rkoio3skM/NnZzLk+Hb84Z+fcf6DH/Dhuu2RLk1EJKJiLtArZLdrxpPXDmTyDwZSVFrGlU/M47YXF7MjvyjSpYmIRERYgW5mI81slZmtNbNfVrG8q5nNMrNFZrbEzM5v+FKrdtZJmbz9s+/wk7NO4PUlmzn7L+/x0oIvcfemKkFEJCrUGOhmFgAeBr4L9ATGmVnPw1b7NfCSu/cHxgKPNHShR5KaFOA/zuvB9FvOoHuHDP7z5SV8f/LHfLnzQFOWISISUeGcoecCa919vbsXA1OAiw5bx4EWoectgc0NV2L4umc258UJQ/j9xb345PNdnPff7/P3jzZSXq6zdRGJf1ZT14SZXQ6MdPfxoelrgMHufnOldToCM4HWQDPgHHdfWMW2JgATADIzM3OmTJlSp6Lz8/PJyMg44jo7Csp5enkxy7aXcXKbBK7vlUL79Oh7yyCctsQKtSX6xEs7QG2pMHz48IXuPrDKhe5+xAcwGniy0vQ1wEOHrXMb8B+h50OAFUDCkbabk5PjdTVr1qyw1isvL/cX5n3up/z2Le/5m3/6ix9/4eXl5XXeb2MIty2xQG2JPvHSDne1pQKwwKvJ1XBOWfOALpWms/h2l8oNwEuhPxAfAalAuzC23ajMjLG5XXnr1jPondWS//zfJUx4dqGuhBGRuBROoM8HuptZNzNLJvim57TD1vkCOBvAzE4mGOjbGrLQ+shqnc7z40/lV+efzHurtjHygQ/4YE3UlCci0iBqDHR3LwVuBmYAKwlezbLczO42s1Gh1f4DuNHMPgVeAH4Q+tcgaiQkGDeeeRxTbxpKy7QkrnnqY+6bvlKfMhWRuJEYzkruPh2Yfti831Z6vgIY2rClNY6enVrw+s2nc8+bK5j0/nrmbtjJX8f1p0ub9EiXJiJSL9F32UcTSEsOcO8lvXnkqgGs35rP+Q9+wFvLtkS6LBGRejkqA73C+b07Mv2nZ3Bcu2b86LlPuOeNFZSUqQtGRGLTUR3oAF3apPPSj4Zw7ZBjefLfGxj3+Fy+1rC8IhKDjvpAB0hJDPC7i3rx0Lj+rNiyl+89+G/mrd8R6bJERGpFgV7JhX07MfWmoTRPTeTKJ+fx9JwNGuRLRGKGAv0wJ2Y257WbhzK8Rwd+9/oKbn95CYUluuWdiEQ/BXoVWqQm8fg1Ofz07O68vDCPMY/P1e3uRCTqKdCrkZBg/OzcE3ns6hzWfL2PUX+dw9K8PZEuS0SkWgr0GozsdQwv/+g0AgnG6Ekf8uYSXa8uItFJgR6Gnp1aMPWmoZzSqSU3Pf8JD89aqzdLRSTqKNDD1L55Cv8zfjAX9evExBmruP3lJRoHRkSiSlhjuUhQalKA+8f0o1u7Ztz/zho27y7g0atzaJmWFOnSRER0hl5bZsat55zIX67oy/yNOxn92Ifk7dK9S0Uk8hTodXTpgCz+dl0uW/YUcukjH7Ji895IlyQiRzkFej2cdkK7g1fAXDHpI+as3R7pkkTkKKZAr6cexzTnlR+fRudWafzg6Y+Z9unhd+cTEWkaCvQG0LFlGi/9aAj9u7bmlhcW8fScDZEuSUSOQgr0BtIyLYm/X5/LeT0z+d3rK5g44zNdqy4iTUqB3oBSkwI8ctUAxg7qwsOz1vGrqcsoK1eoi0jT0HXoDSwxkMB9l/amdbNkHp29jr0FJfzlin4kJ+pvp4g0LgV6IzAzfjHyJFqlJXHfPz8jv6iUR6/KIS05EOnSRCSO6bSxEf3wO8dz36W9eW/1Nq59+mP2FZZEuiQRiWMK9EY2LrcrD4ztzyef7+KqJ+exa39xpEsSkTilQG8Co/p24rGrc/hsyz7GPTGX7flFkS5JROKQAr2JnNMzk8k/GMTnOw5wxaSP+GqP7oAkIg1Lgd6ETu/ejr/fkMvWvUVcMekjDeolIg1Kgd7EBmW34dkbctl1oJgxk+ay9YDGVBeRhhFWoJvZSDNbZWZrzeyX1axzhZmtMLPlZvZ8w5YZX/p3bc0LN57K/uJS7ptXyIbt+yNdkojEgRoD3cwCwMPAd4GewDgz63nYOt2BO4Ch7n4KcGsj1BpXenVuyQs3nkppuTNm0kes25Yf6ZJEJMaFc4aeC6x19/XuXgxMAS46bJ0bgYfdfReAu29t2DLj08kdW/CL3DTKyp2xj89l7dZ9kS5JRGKY1TSAlJldDox09/Gh6WuAwe5+c6V1pgKrgaFAALjL3d+qYlsTgAkAmZmZOVOmTKlT0fn5+WRkZNTptdEmPz+fPaTzx4+DV738IjeVzhmx+dZGvB2XeGhLvLQD1JYKw4cPX+juA6tc6O5HfACjgScrTV8DPHTYOm8ArwJJQDcgD2h1pO3m5OR4Xc2aNavOr402FW1Z8/Vez/n9257z+7d99Vd7I1tUHcXjcYl18dIOd7WlArDAq8nVcE4F84AulaazgMPv4pAHvObuJe6+AVgFdA/rz40AcEKH5kyZcCpmMO6Juaz5Wt0vIlI74QT6fKC7mXUzs2RgLDDtsHWmAsMBzKwdcCKwviELPRqc0CEjFOrGuCfmsXar3igVkfDVGOjuXgrcDMwAVgIvuftyM7vbzEaFVpsB7DCzFcAs4HZ339FYRcez49tn8MKNg4HgmbqufhGRcIX17pu7T3f3E939eHe/NzTvt+4+LfTc3f02d+/p7r3dvW7vdgoQ7H554cbBuDvjHp+r69RFJCyxeTnFUaB7ZnOeD12nPu7xuXy+Q6EuIkemQI9iJ2Y253/GD6aotIxxj8/ly50a+0VEqqdAj3Ind2zBc+MHs7+4jHFPzGXT7oJIlyQiUUqBHgNO6dSS524YzJ6CEq58Yi5b9ijUReTbFOgxondWS/5+fS478ou56ol5bN2r8dRF5FAK9BjSv2trnrluEF/tLeSqJ+fpzkcicggFeowZmN2Gp64dxJe7DnC17lEqIpUo0GPQkOPb8sT3B7J++36umTyPPQUlkS5JRKKAAj1GndG9PZOuzmHVV/u4dvLH5BeVRrokEYkwBXoMG35SB/565QCWbdrD9U/P50CxQl3kaKZAj3EjTjmG+8f2Y8HnOxn/twUUlpRFuiQRiRAFehy4oE8n/uuKvny0fgc/fHYhRaUKdZGjkQI9TlzSP4v7LunNe6u3cfPziygpK490SSLSxBTocWRsblfuvugU3l7xNbdOWUypQl3kqJIY6QKkYX1/SDbFpeXc8+ZKkhMT+PPovgQSLNJliUgTUKDHofFnHEdRaTkTZ6wiOZDAfZf2JkGhLhL3FOhx6qbhJ1BUUsaD/1pLUqLx+4t6YaZQF4lnCvQ49rNzT6SorJxJ760nORDgNxecrFAXiWMK9DhmZvxy5EkUl5Yzec4GkhMT+MXIHgp1kTilQI9zZsZvL+hJSVk5j723juSAcdt5PSJdlog0AgX6UcDMuHtUL0pKPdinHkjgJ2d3j3RZItLAFOhHiYQE4/9d2puSsnL+6+3VJAYS+D/Djo90WSLSgBToR5FAgjFxdF9Ky50/vvUZiQnGjWceF+myRKSBKNCPMoEE4y9X9KWs3Ll3+koCCcb1p3eLdFki0gAU6EehxEAC94/tR1m5c/cbKwgkGNeelh3pskSknjSWy1EqKZDAg+P6c27PTO6ctpxnP9oY6ZJEpJ4U6Eex5MQEHr5yAOecnMlvXlvOc3M/j3RJIlIPYQW6mY00s1VmttbMfnmE9S43MzezgQ1XojSm5MQEHrlqAOec3IFfT13G/8xTqIvEqhoD3cwCwMPAd4GewDgz61nFes2BW4B5DV2kNK7kxAQevmoAZ5/UgV+9ukxn6iIxKpwz9Fxgrbuvd/diYApwURXr/R74E1DYgPVJE0lJDPDI1d+cqatPXST2mLsfeQWzy4GR7j4+NH0NMNjdb660Tn/g1+5+mZnNBn7u7guq2NYEYAJAZmZmzpQpU+pUdH5+PhkZGXV6bbSJtraUlDsPLypi8bYyrj45mXOOTQr7tdHWlvqIl7bESztAbakwfPjwhe5edbe2ux/xAYwGnqw0fQ3wUKXpBGA2kB2ang0MrGm7OTk5XlezZs2q82ujTTS2paikzMf/bb4f+4s3/KkP1of9umhsS13FS1vipR3uaksFYIFXk6vhdLnkAV0qTWcBmytNNwd6AbPNbCNwKjBNb4zGroo3Skeecgx3v7GCx99fF+mSRCQM4QT6fKC7mXUzs2RgLDCtYqG773H3du6e7e7ZwFxglFfR5SKxIymQwENX9ud7fTry/6Z/xl//tSbSJYlIDWr8pKi7l5rZzcAMIABMdvflZnY3wVP/aUfegsSqpEACD4zpR3IggT/PXE1RaTm3nXuixlMXiVJhffTf3acD0w+b99tq1h1W/7IkWiQGgjeaTgoYD/1rLYUlZfzf83XnI5FopLFcpEaBBOMPl/YhNSnAEx9soKCkjLtH9dKNp0WijAJdwpKQYPxu1CmkJQeY9N56DhSX8afL+pAY0OgRItFCgS5hq7hHabPkRP7y9moOFJXxwLh+pCQGIl2aiKBAl1oyM245uzsZKYnc/cYKxv9tAZOuySE9uW4/SlMXbWLijFVs3l1Ap1Zp3D6iBxf379zAVUtj0fGLLvp/Werk+tO78afL+jBn7XaufnIeuw8U13obUxdt4o5XlrJpdwEObNpdwB2vLGXqok0NX7A0OB2/6KNAlzq7YlAXHrlqAMs27WXMpLnsLiyv1esnzlhFQUnZIfMKSsqYOGNVQ5YpjUTHL/oo0KVeRvbqyNPXDSJv1wHumVfIhu37w37t5t0FtZov0UXHL/oo0KXehp7QjhcmnEpRqXP5ox+yJG93WK/r1CqtVvMluuj4RR8FujSIPlmt+NWpaaQmBRj3+FzeX72txtfcPqIHaUmHXiGTlhTg9hE9GqtMaUA6ftFHgS4N5phmCbzy49Po0iad65+Zz/8uzDvi+hf378x9l/amc6s0DOjcKo37Lu2tqyRihI5f9NFli9KgMluk8tKPhvCjZxfyH//4lC17Crhp+AnVDhVwcf/OCoAYpuMXXXSGLg2uRWoSz1yXy8X9OvHnmau545WllJTV7goYEak9naFLo0hOTOC/x/Qjq3U6f521ls17Cnn4yv40Tw3/DkgiUjs6Q5dGY2b8fEQP/nBpb+as3c7lj35E3q4DkS5LJG4p0KXRjc3tyt+uy2XzngIufvhDFn2xK9IlicQlBbo0idO7t+PVH59GenKAMY/P5bXF+ni4SENToEuTOaFDc6beNJR+XVrx0ymL+dNbn1Fe7pEuSyRuKNClSbVplsxzNwxm7KAuPDJ7HTf+fQF7C0siXZZIXFCgS5NLTkzgvkt7c/dFp/De6m1c/PAc1m3Lj3RZIjFPgS4RYWZ8f0g2z40fzO4DJVz01znMWP5VpMsSiWkKdImoU49ry+s/OZ3j2jfjh88u5E9vfUaZ+tVF6kSBLhHXuVUaL/1wyMF+9e9Pnsf2/KJIlyUScxToEhVSkwL84bI+/OmyPizYuIvvPfgB8zfujHRZIjFFgS5R5YpBXXj1x0NJSwow9vG5PDxrrS5tFAmTAl2iTs9OLXj9J6dzfu+OTJyximuf/pit+wojXZZI1FOgS1RqnprEg2P78YdLe/Pxhp2c/8AHzPpsa6TLEolqCnSJWmbG2NyuvP6T02mXkcJ1z8znrmnLKTzsxsQiEhRWoJvZSDNbZWZrzeyXVSy/zcxWmNkSM3vXzI5t+FLlaHViZnDIgOuHduOZDzdywUP/ZmnenkiXJRJ1agx0MwsADwPfBXoC48ys52GrLQIGunsf4GXgTw1dqBzdUpMC/PbCnjx3w2DyC0u55JE5PPjuGt04Q6SScM7Qc4G17r7e3YuBKcBFlVdw91nuXjHQ9Vwgq2HLFAk6vXs7Ztx6Jt/r05G/vL2aSx6Zw6qv9kW6LJGoEE6gdwa+rDSdF5pXnRuAf9anKJEjaZmexANj+/PY1QPYsruQCx76gAfeWUNxqc7W5ehm7ke+xtfMRgMj3H18aPoaINfdf1LFulcDNwPfcfdvfdTPzCYAEwAyMzNzpkyZUqei8/PzycjIqNNro43aUj97i53nVxYxd0sZWRnGdb1SOL5VoN7bjZfjEi/tALWlwvDhwxe6+8AqF7r7ER/AEGBGpek7gDuqWO8cYCXQoaZtujs5OTleV7Nmzarza6ON2tIw3lnxlQ++9x3P/uUb/pupS31PQXG9thcvxyVe2uGutlQAFng1uRpOl8t8oLuZdTOzZGAsMK3yCmbWH5gEjHJ3XSwsTe7skzN5+7YzuXZINs/O/Zxz/us9pn26ueJkQ+SokFjTCu5eamY3AzOAADDZ3Zeb2d0E/1JMAyYCGcA/zAzgC3cf1Yh1i3xL89Qk+nVpxZtLUti6r4hbXljEQ++u4bIBWTw793M27y6gU6s0bh/Rg4v7H+ltoLqZumgTE2esavT9hOPXU5fywrwvubVXCTfcMZ1xg7twz8W9I1KLNJ0aAx3A3acD0w+b99tKz89p4LpEam3qok3c8cpSCip98GjN1nz+8NZnB6c37S7gjleWAjRo2B6+78baTzh+PXUpz8394uB0mfvBaYV6fNMnRSVuTJyx6pAwr05BSRkTZ6xq9H03xn7C8cK8L2s1X+KHAl3ixubdBWGvu2l3QYP2r1e379rU1FDKqmlXdfMlfijQJW50apVWq/VHP/YRCxpozPXq9l3bmhpCIPg+VtjzJX4o0CVu3D6iB2lJh16DnpRgJAUODbLUxARG52Tx+c4DXP7YR1z/zHyWbarf2DBV7TstKcDtI3rUa7t1MW5wl1rNl/gR1puiIrGg4s3Hw680qWrexf07c6C4lKfnbOTx99dzwUP/ZsQpmdxydndO6dSywfYdiatcKt74rOgzD5jpKpejhAJd4srF/TtXGaJVzUtPTuSm4SdwzZBjeeqDDUz+9wZmLP+ac3tmclrLMoY10L4j4Z6Le3PPxb2ZPXs2664aFulypIko0OWo1yI1iZ+deyLXn96Np+cEg/3twlLe3TaPHw87niHHt8XU/ywxQH3oIiEt05K49ZwTmfPLs7jixCQ++2ofVz45j1F/ncO0TzdTqqF6Jcop0EUO0zw1ifOPS+bfvxjOfZf2Zn9xKbe8sIgz/jSLR2evY/eB4kiXKFIldbmIVCM1KcC43K6MGdiFf322lac/3MAf3/qM+99ZzUX9OnHNqdn0zqr9G6gijUWBLlKDhATjnJ6ZnNMzk8++2svfPvycqYs28dKCPPpktWRcblcu7NuJjBT9OklkqctFpBZOOqYF913am3m/Ops7L+xJUUk5d7yylNx73+H2f3zK/I07NcKjRIxOKUTqoEVqEtcN7cYPTsvmky928+L8L3hzyRb+sTCPY9umc3G/zlzSvzPZ7ZpFulQ5iijQRerBzMg5tjU5x7bmzgtP4Z/LvuKVT/J48F9reODdNfTt0ooL+3Tkgj6dOKZlaqTLlTinQBdpIM1SErk8J4vLc7LYsqeAaYs38/qSzdzz5krueXMlOce25ru9jmHEKcfQpU16pMuVOKRAF2kEHVum8cPvHM8Pv3M867flM33pFt5c+tXBcO/ZsQXn9szknJMz6dW5hT64JA1CgS7SyI5rn8HNZ3Xn5rO6s3H7fmau+IoZy78+2C2T2SKFYSd2YFiP9gzt3o4WqUmRLllilAJdpAllt2vGhDOPZ8KZx7Mjv4jZq7bx7mdfM33ZFl5c8CWBBKNfl1ac0b0dQ09oR9+sViQn6mI0CY8CXSRC2makcFlOFpflZFFSVs6iL3bz/uptfLBmGw+8u4b731lDWlKAgdmtOfW4tgzu1obeWS1JSQzUvHE5KinQRaJAUiCB3G5tyO3Whp+P6MHuA8XM27CTD9duZ+76nQdvZZecmEC/rFbkZLdmQNfW9O/ainYZKRGuXqKFAl0kCrVKT2bEKcErYgB27i/m4w07WLBxF/M/38UT76+ntDz4AaYubdLom9WKfl1a0btzS07p3FKfWj1K6aiLxIA2zZIZ2asjI3t1BKCguIylm/bwyRe7WJK3m0Vf7OaNJVsAMINubZtxcqcWpBUW48ds5eSOLchskaKraeKcAl0kBqUlBw520VTYtlCiT2kAAAxJSURBVK+IZZv2sHTTHpZv3sOSvN18ubOEl1fPB4LDA/fIbE73zAy6d8jghA7NOaFDhoI+jijQReJE++YpDD+pA8NP6nBw3vS3Z9G+e19WbtnLqq/2sfrrfbz+6Wb2FpYeXKdZcoDj2meQ3a4Z3dqmk92uGce2bcaxbdNp2yxZYR9DFOgicSw9yRiU3YZB2d+cybs72/KLWPt1Puu25bNu237Wbctn8Ze7eGPJZiqPLZaeHKBrm3SyWqeT1TqNrNZpdG6VRqfQo22zZBISFPjRQoEucpQxMzo0T6VD81ROO6HdIcuKSsv4cmcBX+zcz8btB/hy1wG+3HmAvF0H+GjddvYXlx2yfnIggcyWKRzTIpXMFqkc0yKVDi1SyGyRSvvmKXRonkL7jFRapCXqTL8JKNBF5KCUxAAndMjghA4Z31rm7uwpKGHT7gI27y5ky55vvn61p5Blm/bwzsqvKSz59q36kgMJtM1Ipm1GMm2apdC2WTJtKj1apyfTOj2J1s2SaZWWRMv0JF1vXwcKdBEJi5nRKj2ZVunJnNKp6js1uTv7ikrZureIbfuK2LqvkG37itieX8z2/CJ27i9mR34R67bms+tAMQcOO+OvLDUpgZZpSbRMS6JFahIlBwp59atFtEhNonlqIhmpiTRPTaJ5SiLNUhJplhIgI/S84mt6UuCo6hIKK9DNbCTwABAAnnT3Pxy2PAX4O5AD7ADGuPvGhi1VJH5NXbSJiTNWsXl3AZ1apXH7iB78Y8EXzFm38+A6Q49vw+iBXb+1HvCteQs+38kL877k1l4l3HDHdMYN7sI9F/cOa78X9+9c7fxwXl+x7zJ3AmZH3Pef3vqMzXsKaZ+RwuiBWRwoLmPqok3sLighMSGB7LbNaJ2ezN7CErYVOYu+2M2+whL2FJRQHuZ9RFKTEkhPTiQtKUB6cvCRlhw4OC8lKYG0pACpSQFSkxJITfzmeUpicHlKYuh5YgLJoUdKYuDg86SAkRyoeJ5AYoJFpIupxkA3swDwMHAukAfMN7Np7r6i0mo3ALvc/QQzGwv8ERjTGAWLxJupizZxxytLKSgJnq1u2l3ArS8u/tZ6c9btPCTgN+0u4PaXPwWHklC6bdpdwG0vLqZyp0eZO8/N/QLgkGCtar93vLKUBZ/v5H8XbvrWfOCQUK/q9XXd97b8Ip54fz0YlJQF25JfVMqSvD3cd2lvLu7fmdmzZzNs2DCmLtrEL/93CYWl3+wpJZDA+DO60a9ra/YXlbK/uJT9RaUcKC7jQHEZ+4tKKQg9P1BSRmFxGVv3FVJYUk5BcRmFJaFHaTll4f6lqEFyIBj0SYkJJCYkkBwwEkPzBrUpYdiwBtnNIcI5Q88F1rr7egAzmwJcBFQO9IuAu0LPXwb+ambmuheXSI0mzlh1MNhqqyL8Kvt2D3bQC/O+PCRUq9pvQUnZwbPrw+dPnLHqkECv6vX12XdJFUFa3X4rhzlAUVk5Uxdv5vaRJ1VTQfhKysopLCmjqLQ8FPTlFJeWU1QanBd8Xk5JWfB5cWk5xaHnJZW+lpQ7xaXllJaVU1zmlJQFn5eUOy3YWXMhdWA1Za6ZXQ6MdPfxoelrgMHufnOldZaF1skLTa8LrbP9sG1NACaEJnsAq+pYdztge41rxQa1JTo1WVuSjzkhp7G2XXZgD4H0b/q7i79au7A++63P6xvgte2A7Ud6beV9RLn6/Hwd6+7tq1oQzhl6VR1Bh/8VCGcd3P1x4PEw9nnkgswWuPvA+m4nGqgt0Sle2mJmC0r3bI35dkD8HBNovLaEM9ByHtCl0nQWsLm6dcwsEWgJjfQ/hYiIVCmcQJ8PdDezbmaWDIwFph22zjTg2tDzy4F/qf9cRKRp1djl4u6lZnYzMIPgZYuT3X25md0NLHD3acBTwLNmtpbgmfnYxiyaBui2iSJqS3SKl7bESztAbalRjW+KiohIbNDNCkVE4oQCXUQkTkR9oJtZqpl9bGafmtlyM/tdaH43M5tnZmvM7MXQG7ZRz8wCZrbIzN4ITcdqOzaa2VIzW2xmC0Lz2pjZ26G2vG1mrSNdZzjMrJWZvWxmn5nZSjMbEottMbMeoeNR8dhrZrfGaFt+Fvp9X2ZmL4RyIFZ/V34aasdyM7s1NK9RjknUBzpQBJzl7n2BfsBIMzuV4PAC/+3u3YFdBIcfiAU/BVZWmo7VdgAMd/d+la6n/SXwbqgt74amY8EDwFvufhLQl+Dxibm2uPuq0PHoR3BcpQPAq8RYW8ysM3ALMNDdexG8GKNiSJGY+l0xs17AjQQ/cd8XuMDMutNYx8TdY+YBpAOfAIMJfsoqMTR/CDAj0vWFUX9W6OCdBbxB8ANZMdeOUK0bgXaHzVsFdAw97wisinSdYbSjBbCB0AUCsdyWw+o/D5gTi20BOgNfAm0IXon3BjAiFn9XgNEEBzSsmP4N8J+NdUxi4Qy9optiMbAVeBtYB+x294r7aOUR/CGIdvcTPJgVA1G0JTbbAcFPAs80s4WhIR0AMt19C0Doa4dqXx09jgO2AU+HusKeNLNmxGZbKhsLvBB6HlNtcfdNwJ+BL4AtwB5gIbH5u7IMONPM2ppZOnA+wQ9hNsoxiYlAd/cyD/4bmUXwX5eTq1qtaauqHTO7ANjq7pXHmghryIQoNdTdBwDfBW4yszMjXVAdJQIDgEfdvT+wnyjvkqhJqG95FPCPSNdSF6H+5IuAbkAnoBnBn7PDRf3viruvJNhV9DbwFvApUHrEF9VDTAR6BXffDcwGTgVahYYZgKqHI4g2Q4FRZrYRmEKw2+V+Yq8dALj75tDXrQT7aXOBr82sI0Do69bIVRi2PCDP3eeFpl8mGPCx2JYK3wU+cfevQ9Ox1pZzgA3uvs3dS4BXgNOI3d+Vp9x9gLufSfCDl2topGMS9YFuZu3NrFXoeRrBg70SmEVwmAEIDjvwWmQqDI+73+HuWe6eTfDf4X+5+1XEWDsAzKyZmTWveE6wv3YZhw4BERNtcfevgC/NrEdo1tkEh4aOubZUMo5vulsg9tryBXCqmaWbmfHNMYm53xUAM+sQ+toVuJTgsWmUYxL1nxQ1sz7A3wi+050AvOTud5vZcQTPdNsAi4Cr3b0ocpWGz8yGAT939wtisR2hml8NTSYCz7v7vWbWFngJ6Erwl3K0u0f9IG1m1g94EkgG1gPXEfpZI/bakk7wDcXj3H1PaF7MHZfQ5cljCHZPLALGE+wzj6nfFQAz+4Dg+2UlwG3u/m5jHZOoD3QREQlP1He5iIhIeBToIiJxQoEuIhInFOgiInFCgS4iEifCuUm0SJMKXdL1bmjyGKCM4MfzAXLdvTgihR2BmV0PTA9d1y4SEbpsUaKamd0F5Lv7n6OgloC7l1Wz7N/Aze6+uBbbS6w0NolIvanLRWKKmV1rwfHxF5vZI2aWYGaJZrbbzCaa2SdmNsPMBpvZe2a23szOD712vJm9Glq+ysx+HeZ27zGzj4FcM/udmc0PjW/9mAWNITi084uh1yebWV6lTzifambvhJ7fY2aTzOxtggOCJZrZX0L7XmJm45v+uyrxQoEuMSM0tvQlwGmhwdoS+eaG5C2BmaEBw4qBuwh+ZHw0cHelzeSGXjMAuNLM+oWx3U/cPdfdPwIecPdBQO/QspHu/iKwGBjjwfHIa+oS6g9c6O7XABMIDtqWCwwiONBZ17p8f0TUhy6x5ByCobcgOMQHaQQ/5g5Q4O5vh54vBfa4e6mZLQWyK21jhrvvAjCzqcDpBH8PqttuMd8McwBwtpndDqQC7QgO6/rPWrbjNXcvDD0/DzjZzCr/AelO8OPgIrWiQJdYYsBkd//NITODI/BVPisuJ3inq4rnlX/OD3/TyGvYboGH3mgKjZPyV2CAu28ys3sIBntVSvnmP+DD19l/WJt+7O7vIlJP6nKRWPIOcIWZtYPg1TB16J44z4L3EE0nOOb2nFpsN43gH4jtodEmL6u0bB/QvNL0RoK3geOw9Q43A/hxxbCwFrwvaFot2yQC6AxdYoi7Lw2NwveOmSUQHL3uR9RuXOx/A88DxwPPVlyVEs523X2Hmf2N4FDBnwPzKi1+GnjSzAoI9tPfBTxhZl8BHx+hnkkER9xbHOru2UrwD41IremyRTlqhK4g6eXut0a6FpHGoC4XEZE4oTN0EZE4oTN0EZE4oUAXEYkTCnQRkTihQBcRiRMKdBGROPH/AQog38hyUT00AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1.2])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing and plotting uncertainty" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/aschmide/miniconda3/envs/envStats9/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXgUVb4//ndV9ZZ96exhExQMskoEFRUHEFARUK/iD6/LKKBXrjpcmStXHBZBeXDujDKiMs4gI4PL/HCFyIACM1c2AyiyCKIDgbB0ts7e6bXqfP/opKVZO0WW7s779Tx5SHeqqj+HXt5d51SdkoQQAkRERM0kt3cBREQUmRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLq0SYAsWrQIw4cPR69evfDjjz+ecxlVVTFv3jyMHDkSt9xyC1atWtUWpRERkU5tEiAjRozAu+++i9zc3PMus2bNGhQXF+OLL77A3/72N7z22ms4ceJEW5RHREQ6tEmA5OfnIzs7+4LLrF27Fvfccw9kWUZqaipGjhyJdevWtUV5RESkQ9iMgdhsNuTk5ARuZ2dno6SkpB0rIiKiCwmbACEioshiaO8CmmRnZ+PUqVPo168fgLP3SEJVVeWApkXn9F5Wazzs9vr2LqPVRHP7orltANsXyWRZQkpKnK51wyZAxowZg1WrVmHUqFGorq7Ghg0b8O677zZ7O5omojZAAER124Dobl80tw1g+zqiNunCWrBgAW666SaUlJTgl7/8JW6//XYAwJQpU7Bv3z4AwPjx49GpUyeMGjUK9957L6ZNm4bOnTu3RXlERKSDFG3Tudvt9VH7TSE9PQHl5XXtXUarieb2RXPbALYvksmyBKs1Xt+6LVwLERF1EAwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6GNrqgYqKijBz5kxUV1cjOTkZixYtQrdu3YKWsdvt+J//+R/YbDZ4vV5ce+21eP7552EwtFmZREQUojbbA5kzZw4mTZqE9evXY9KkSZg9e/ZZyyxduhQ9evTAmjVrsGbNGnz//ff44osv2qpEIiJqhjYJELvdjgMHDmDs2LEAgLFjx+LAgQOorKwMWk6SJDgcDmiaBo/HA6/Xi8zMzLYokYiImqlN+oZsNhsyMzOhKAoAQFEUZGRkwGazITU1NbDcE088gSeffBI33HADnE4n7r//fgwaNKhZj2W1xrdo7eEmPT2hvUtoVdHcvmhuG8D2dURhNbiwbt069OrVC++88w4cDgemTJmCdevWYcyYMSFvw26vh6aJVqyy/aSnJ6C8vK69y2g10dy+aG4bwPZFMlmWdH/xbpMurOzsbJSWlkJVVQCAqqooKytDdnZ20HIrV67EuHHjIMsyEhISMHz4cBQWFrZFiURE1ExtEiBWqxV5eXkoKCgAABQUFCAvLy+o+woAOnXqhK+++goA4PF4sH37dlxxxRVtUSIRETVTmx2FNXfuXKxcuRKjR4/GypUrMW/ePADAlClTsG/fPgDAc889h2+++QZ33HEHJkyYgG7duuHee+9tqxKJiKgZJCFEVA0YcAwkckVz+6K5bQDbF8nCfgyEiIiiDwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6RJygKxYsQKVlZWtWQsREUWQkANk27ZtGDFiBB577DGsXbsWHo+nNesiIqIwF3KALF26FJs2bcJNN92Ed955B0OHDsWsWbOwc+fO1qyPiIjCVLPGQFJSUnD//ffjb3/7G/76179i3759ePDBBzF8+HC8+eabcDgc5123qKgIEydOxOjRozFx4kQcPXr0nMutXbsWd9xxB8aOHYs77rgDFRUVzWoQERG1DUNzV9i+fTtWr16NjRs3ok+fPpg8eTJycnKwYsUKTJkyBe+9994515szZw4mTZqE8ePH47PPPsPs2bOxYsWKoGX27duHJUuW4J133kF6ejrq6upgMpn0tYyIiFpVyAGyaNEifP7550hISMD48eOxZs0aZGZmBv7ev39/DB48+Jzr2u12HDhwAMuXLwcAjB07FvPnz0dlZSVSU1MDy/3lL3/BI488gvT0dABAQkKCrkYREVHrCzlA3G43lixZgn79+p3z70ajER9++OE5/2az2ZCZmQlFUQAAiqIgIyMDNpstKEAOHz6MTp064f7770dDQwNuueUW/Md//AckSWpOm4iIqA2EHCCPPfYYLBZL0H01NTVwuVyBPZEePXpcUjGqquLQoUNYvnw5PB5PoHtswoQJIW/Dao2/pBrCXXp6dO+VRXP7orltANvXEYUcIE888QReeuklJCUlBe4rKSnB888/j1WrVl1w3ezsbJSWlkJVVSiKAlVVUVZWhuzs7KDlcnJyMGbMGJhMJphMJowYMQJ79+5tVoDY7fXQNBHy8pEkPT0B5eV17V1Gq4nm9kVz2wC2L5LJsqT7i3fIR2EVFRWhV69eQff16tULR44cuei6VqsVeXl5KCgoAAAUFBQgLy8vqPsK8I+NbNmyBUIIeL1efP3117jyyitDLZGIiNpQyAFitVpx7NixoPuOHTuG5OTkkNafO3cuVq5cidGjR2PlypWYN28eAGDKlCnYt28fAOD222+H1WrFbbfdhgkTJuDyyy/Hv/3bv4VaIhERtSFJCBFSf8/SpUuxdu1aTJ8+HZ07d0ZxcTEWL16MW2+9FY8//nhr1xkydmFFrmhuXzS3DWD7ItmldGGFPAYydepUGAwGLFq0CCUlJcjKysI999yDX/7yl7oemIiIIlvIASLLMiZPnozJkye3Zj1ERBQhmnUm+pEjR/DDDz+goaEh6H6OUxARdTwhB8jSpUvx+uuv48orrww6H0SSJAYIEVEHFHKAvPPOO1i1ahUPqyUiIgDNOIzXYrGge/furVkLERFFkJAD5Omnn8aCBQtQVlYGTdOCfoiIqOMJuQtr5syZABA0bYkQApIk4eDBgy1fGRERhbWQA2Tjxo2tWQcREUWYkAMkNzcXAKBpGioqKpCRkdFqRRERUfgLeQyktrYWzzzzDPr164dRo0YB8O+VvPLKK61WHBERha+QA2TOnDmIj4/Hpk2bYDQaAQADBw7E3//+91YrjoiIwlfIXVjbt2/H5s2bYTQaA1cITE1Nhd1ub7XiiIgofIW8B5KQkICqqqqg+06dOhW4fjkREXUsIQfIPffcg6eeegpff/01NE3D7t278eyzz+K+++5rzfqIiChMhdyFNWXKFJhMJrzwwgvw+Xx47rnnMHHiRDz00EOtWR8REYWpkANEkiQ8/PDDePjhh1uxHCIiihTNGkQ/n+uuu65FiiEiosgRcoDMmjUr6HZVVRW8Xi8yMzN5ljoRUQcUcoBs2rQp6LaqqnjzzTcRFxfX4kUREVH4C/korDMpioLHH38cf/7zn1uyHiIiihC6AwQAtm7dGjipkIiIOpaQu7CGDRsWFBZOpxMejwdz5sxplcKIiCi8hRwgv/3tb4Nux8TE4LLLLkN8fHyLF0VEROEv5AAZPHhwa9ZBREQRJuQA+fWvfx3SeMfLL798SQUREVFkCHkQPTExERs2bICqqsjKyoKmadi4cSMSExPRpUuXwA8REXUMIe+BHD16FG+99Rby8/MD9+3atQtvvvkmli1b1irFERFR+Ap5D+S7775D//79g+7r378/du/e3eJFERFR+As5QHr37o3f//73cLlcAACXy4VXXnkFeXl5rVYcERGFr5C7sBYuXIgZM2YgPz8fiYmJqK2tRZ8+fc46vJeIiDqGkAOkU6dO+OCDD2Cz2VBWVob09HTk5OS0Zm1ERBTGmjWVSVVVFQoLC7Fjxw7k5OSgtLQUJSUlrVUbERGFsZADZMeOHRgzZgzWrFmDN954AwBw7NgxzJ07t7VqIyKiMBZygLz00kt49dVXsWzZMhgM/p6v/v37Y+/eva1WHBERha+QA+TkyZOBKw82nZFuNBqhqmpI6xcVFWHixIkYPXo0Jk6ciKNHj5532SNHjqB///5YtGhRqOUREVEbCzlAevTogc2bNwfdt23bNvTs2TOk9efMmYNJkyZh/fr1mDRpEmbPnn3O5VRVxZw5czBy5MhQSyMCAAgBCIj2LoOowwj5KKyZM2fisccew8033wyXy4XZs2dj06ZNgfGQC7Hb7Thw4ACWL18OABg7dizmz5+PyspKpKamBi371ltv4eabb0ZDQwMaGhqa2Rzq6Dw+DWaD0t5lEHUIIQfIgAEDsHr1aqxevRp33303srOz8eGHHyIrK+ui69psNmRmZkJR/G9sRVGQkZEBm80WFCA//PADtmzZghUrVoQUTOditUb39PLp6QntXUKrupT2qZpAdZ0L1qSYFqyo5fC5i2zR3j49QgoQVVXx8MMPY9myZZgyZUqrFOL1evGb3/wGCxcuDASNHnZ7PTQtOrsx0tMTUF5e195ltJpLbZ8QQGWdCz63F3KYXSmTz11ki+b2ybKk+4t3SAGiKApOnDgBTdN0PUh2djZKS0uhqioURYGqqigrK0N2dnZgmfLychQXF2Pq1KkAgNraWgghUF9fj/nz5+t6XOp4NE3A49NgMbIbi6i1hdyFNW3aNMydOxdPPvkksrKygq4NIssXHou3Wq3Iy8tDQUEBxo8fj4KCAuTl5QV1X+Xk5KCwsDBw+7XXXkNDQwOeffbZ5rSHOjgBwOH0wGyMQXjtgxBFn5AD5PnnnwcAfPrpp4HwEEJAkiQcPHjwouvPnTsXM2fOxBtvvIHExMTAIbpTpkzBU089hb59++qpn+gsXp+A26tyL4SolUlCiAsOGJSXlyM9PR0nT5487zK5ubktXpheHAOJXC0xBlJe44SmCRgMEqwJMQiXoRA+d5Etmtt3KWMgFz0PZPTo0QD8IZGbm4uFCxcGfm/6IQo3Pp+A0+tr7zKIotpFA+TMHZQdO3a0WjFELamhwcsTC4la0UUDRAqXPgCiZvJpAg3u0KbaIaLmu+gguqqq+PrrrwN7Ij6fL+g2gMAcWUThxuH0IsakhN15IUTR4KIBYrVa8dxzzwVuJycnB92WJAkbN25sneqILpGmCThcPiTEGNu7FKKoc9EA2bRpU1vUQdRqnC4fYs0KlIucr0REzRPyeSBE4Wrv4QqsKyyGy6MizmJA/pUZ6NUlJfB3TQjUOb1IjjMBPL0wojU91xU1LqQlWTBmSBf065HW3mV1WPxKRhFt7+EKvPvlj6h2eBBjUVDr9GL11iIcKq4KWs7tVuH26puKh8LD6c91rMWAaocH7375I/Yermjv0josBghFtHWFxVAUGWajAkmSYDIoUBQZm/ecClpOAKhr8EC78HmzFMbOfK7NRv9zva6wuL1L67AYIBTRKmpcMBmCX8ZGRUZVnfusZX2qQG2DB+C5IRHpXM+1ySCjosbVThURA4QiWlqSBR5fcNeUV9WQkmA+5/Iutwqnh+eGRKJzPdcen4a0JEs7VUQMEIpoY4Z0gapqcHtVCCHg8alQVQ039s857zoOJ89Qj0RnPtdur/+5HjOkS3uX1mHxKCyKaE1H4KwrLIbTpSIxxojhA3ODjsI6k0/1n6EeZ+bLP5Kc/lzzKKzwwHcQRbx+PdLQr0da0Gy8F9Pg9CLGKF/0WjYUXpqeawoPfPdQh6RqArUN3vYugyiiMUCow3J5VDS4OeU7kV4MEOrQ6p1eaBpPMCTSgwFCHZqmCdQ6vWFz5UKiSMIAoQ7P7VHh5rkhRM3GAKEOTwig1ulp7zKIIg4DhAj+a6g7OKBO1CwMEKJGDqcXKgfUiULGACFqpGkCdQ1ecLJFotAwQCgq2Gtc2LDrOOqdl3ZyoNvDyRaJQsWpTCgqrN5ahM17bUhJMOOhMVciIyVG13YEgDqHF4osnzV1OBEF4zuEosLQvtkwGf3XAVn62X7860SN7m1pQqDW4eaMvUQXwQChqNCzczKenXQ1EmKNcHlU/OXvB1F4oFT39nyqQIOLR2URXQgDhKJGt6xETLurL7KtsdAE8NmWIqzZdhRqCLPznovD5eM0J0QXwAChqJIcb8bUcVchr6v/eiDb95fgnb//AKeOczw0TaDawYtPEZ0PA4Sijtmo4P5beuKmxqsS/utkDd74ZD9KKxuavS2PV0V1Hc9SJzoXBghFJVmWMGZIF9z7i8thUCTYa11487P92F9U2extub0qahs8nHCR6AwMEIpqA65Iw9RxVyEpzgSPV8N7X/6I9TuKQ7pq4emcLh+nOiE6AwOEol6n9HhMu6svLstOBAD833ensPzvB5t10qH//BAPXF6eZEjUhAFCHUJ8jBGP3J6HG/pmAwAOn6zFko/34VhJXcjbEAKoqXfD6fGB050QtWGAFBUVYeLEiRg9ejQmTpyIo0ePnrXM66+/jttvvx3jxo3DXXfdhc2bN7dVedQBKLKE267riv9v5BUwGWXUOjz405oD+GrPKWgitEDwh4gH1fWekNchilZtFiBz5szBpEmTsH79ekyaNAmzZ88+a5l+/frhww8/xOrVq/HSSy9h+vTpcLlcbVUidRB9u1sx7c6+yEyJgSYE1hUW46/rDjWrS8vlURki1OG1SYDY7XYcOHAAY8eOBQCMHTsWBw4cQGVl8BExN954I2Ji/HMY9erVC0IIVFdXt0WJ1MGkJ8fgP+7sg/xe6QCAQ8erseSjvTh8KvQpUDxeFVV1LmiCJxtSx9QmAWKz2ZCZmQlFUQAAiqIgIyMDNpvtvOt8+umn6NKlC7KystqiROqATAYFdw3rgXuHX+7v0mrw4u2Cg1i/oxg+NbRQ8PoEKmvd8HhVHuZLHU5Yzsa7Y8cOLF68GG+//Xaz17Va41uhovCRnp7Q3iW0qktpn6oJ+GQJze1VGj44Dn2uSMey1d/jmK0W//fdKRTZ6vDIuKuQZY0LeTuS0YDkeBMMyrm/l/G5i2zR3j492iRAsrOzUVpaClVVoSgKVFVFWVkZsrOzz1p29+7d+PWvf4033ngD3bt3b/Zj2e31zT7GP1KkpyegvDz0o4YizaW2Twigqsap6/k3AJh8+5XYsOsEvvruFIpL67Dg7UKMGdIF116VBTnE3YsyRUJinPmsqeD53EW2aG6fLEu6v3i3SReW1WpFXl4eCgoKAAAFBQXIy8tDampq0HJ79+7F9OnT8Yc//AFXXXVVW5RGFKDIMkYP7oLJd/RGcrwJPlWgYNsxvP35QVTVhXYwh08VqK5zN54vEp1fZIiaSEK0zWEkhw8fxsyZM1FbW4vExEQsWrQI3bt3x5QpU/DUU0+hb9++uPvuu3Hy5ElkZmYG1nv55ZfRq1evkB+HeyCRqyX2QMp17oGcyeXx4fNtx/DNj+UAAJNRxpghXTA4LzOkvREJgNmkID7WCIMs87mLcNHcvkvZA2mzAGkrDJDIFU4B0mTDruP4v+9OBaaEz0yJwfV9srDnXxWoqnMjJcGMG/vnoFeXlHOuL0sSLGYFXXJTUF3l0FXD3sMVWFdYjIoaF9KSLBgzpAv69UjT3aaWtHrLEXyx8wRcXhUWo4JR13TCuBua3/Uc7qL5vRf2XVhEkehQcRV2/1SOpHgTYsz+IwhLq5z4ZHMRSqqcMJsU1Dq9WL21CIeKq865DU34L0xVUd0Ap8fX7E6tvYcr8O6XP6La4UGsxYBqhwfvfvkj9h6uuMTWXbrVW45g9bajcHtVGGT/pJOrtx3F6i1H2rs0aiMMEKLz2LznFBRFhsVkQEqCBamJ5sDfGlw+VNS4AAEoiozNe05dcFta4xnslbVOuJsxPrKusBiKIsNsVCBJEsxGBYoiY11h8aU0rUV8sfMEJEhQZAmSJPv/hYQvdp5o79KojYTlYbxE4aCqzg2L+ee3iMVkgAR34KPfpwpU1LgQY1bgDXGSRa/PP8huMiqIsxhgMioXXL6ixoVYS/Db1GSQ/eHVzlweHxQ5eDxIlvz3U8fAPRCi80hJMMN7xgmFBkWCUZGQlmSBsfF8D6dbRV2DF18fKAlp/EXA391TVeeGvdZ5wa6ttCQLPL7gGjw+DWlJFj1NalEWkwFnNlcT/vupY2CAEJ3Hjf1zoKoaPD4VQgh4fCpMRgUmkwGQAGuSGXExBkjwf3Cu3nIUr3+yD0dO1Ya0fQH/HklT15bLq0ITIuiM9jFDukBVNbi9/hrcXhWqqmHMkC6t0ubmGHVNJwgIqJqAEJr/XwiMuqZTe5dGbUSZO3fu3PYuoiU5nZ5mn4kcKeLizGhoiN7Lq7ZE+xrcvhZ7/tOSYpCWZEGpvQH1DV4kx5lw67Vd0btbiv8+pxdpSRbcck1nWEwG2Brv+/bHcpTYG5CbFhfofoqJMcF5gckaNc0/QaPTo8Lr1QBZgkGRkZkai8yUGJwoq0eNw4PUBDPuuql7WByF1atLCiAEjpXUw6MKWIwKbru2S1QehRXN7z1JkhAba9K3Lg/jjRzRfCghEJ6H8TbHibJ6FGw/iuLSegD+6eOH9M7EL67OReecZFRWNu8wXpNRQXyMEWajHPZfivjajFw8jJcoDHTKiMdj467CxOGXIzneBFUT2La/BP/7/ndYu7Wo8eir0DXN9ltV74bHq/G8dgo7HO0iakGSJKH/5Wno3S0V278vwT93n4TLo2L15iPYuOs4bh6Qg8F5mTAaQvvuJgTgcqtwu1UoBgmxZiPMBhkGQ/jvlVD0Y4AQtQKjQcZN/XOQ3ysDX+05ie3fl8Lh9OLz7cewea8NNw/IQf6VGeedufdMAoDPJ1Dr80CSAKNBgcWs+MNEYZhQ++AYSASJ5n5YIPLHQC5EMir4ZNNP+OZQeWBalMQ4E27qn438KzNgMlz4fJDzkSUJJoMMi8UAs1GGhPa5KAlfm5HrUsZAuAdC1AZSEiyYcGN3DBuQg3/sPoVvD5Wj1uFBwbZj+MfuUxjaJwtDemcixty8t6QmBFxeFS6vClmWYDEpMBkVGBQJBplDnNS6GCBEbSglwYK7buqOXwzMwf99dwrfHCqHw+nFFzv9kzZec2UGru+bheR488U3dgZN88+71eDyQZL8R4GZTQb/9Ceyf7p6opbEACFqB017JMOv7oSt+2woPFgKt1fFln02bNtvQ5/uVgztm4XOGfqugieEf6oVn9MLh9PrDxRFgqUxUIwhjr0QXQgDhKgdJTaenHjzwFwUHijF9v0lqHN6sfewHXsP29EpPQ7X9clC3+7WkAfcz0UI/yB8vc8Lh+SFIkv+s+obw8SgNP9SwEQcRI8g0TyQB0T3IHpqalxIJxL6VA17D9uxdZ8NNntD4P5YiwH5vdJxTV4mrIktOw+WLEkwGGQYG38MigRFlps1HM/XZuTiIDpRlDAoMq7umY6BV6ThaEkdtn9fggNFlWhw+fDVHhu+2mNDj9xE5PfKQO9uqSGfT3IhmhDweFV4Gk90lCRAkSQYDQoMRhlGWYIs+0NFlsE9FQpggBCFIUmScFl2Ii7LTkStw4OdP5Rh1w9lqHF4cPhkLQ6frIXFpKD/5Wm4umcaOqXHQwrhUruhEALwCQGfxwd4murx16TIkr/LyyDD0BQsigTROAkkw6VjYYAQhbnEOBNGDOqEXwzMxU8nqrHzhzL8cKwaLo+KwgOlKDxQirQkC/pfnob+l1uRlhTT4jUIAQghoGkCXp8GuP33NwWLKsmoqXHBIMvBey2KBFni+Eq0YoAQRQhZltCrSwp6dUlBXYMHe/5l98/8W9mAihoXNn5zAhu/OYHctDj07WFF3+6pSElo3euGBIJFAF6fBi+0c+y1AEZFgaz492AU2T/OIktSYOp6BkxkYoAQRaCEWBNu6JeNG/plw2Z34LufKrDnsB21Dg9OVjhwssKBdYXFyE2Lw1WXpaJ3t1RkpLT8nsmF/LzXAnh9P1+lsClYJACQ/FcxNMhN4y0ylMagofDHACGKcNnWOGRb4zB6SBccK6nD3sN2fF9UiXqnNxAmX+w8jrQkC/K6puDKrinokpnQbh/STcHSRAXgxc/jLbIkQZb9BxQoigxZliA3nhjp/1vw3suZ2z4f7u20PAYIUZSQTxt4v+P6bjhaUof9RXYcPFqFGocHFTUubN5rw+a9NsSYFVzRKRk9Oyfj8k5JSNR5QaHWoAkBTQV8qgp/vPxMkgAJUuNeTNN9EhT/H/y3ERwuQgACAk3z4UuSBFnxh5J0WhhJ8N93oYCiYAwQoigkyxK65ySie44/TE5WOHDwaBV+KK6Czd4Ap1sNnKwIAFmpsbi8UxJ65CSiW3YizEZ9kzu2tjPDoPFenP9aj81zevfa6SElFAV1Dd7AOE7TsiZF6dBBwwChqGI2Kj93jwhAQ9OvTR864ufPnp/vwml3neX0z4fTPyy007tizvxMCyOSJKFTejw6pcfjlms6o7rejUPF1fjpRDX+dbIGHq+GksoGlFQ2YMteG2RJQqeMOHTPScJl2QnokpEAsyk8A6Wlndm91sTj0+BwBceULEtIS7K02wzI4YABQlFDkoCkOFPQ7VCdlgPBv/zcM3LW9oQAVE2DEP4wgRDwqBo8Xg0+VQvLM+IBIDnejCG9MzGkdyZ8qobjZfX46UQNDp+swYnyemhCoLi0HsWl9fjnbv8gd3ZaHLpmJqBrVgK6ZCYE/T93WOH59LYpBghFLT2DpdJZv1x4e2fOcGsyKpBi/P34bq8Gj0eFpomw7eYwKHJg3ATXdIbL40ORrQ5HTtWgyFYHm90BTQAnyx04We7Atv0lAPznpnROj0fnjHjkZsQhJq75swdT5GOAEPJL8gMAABIESURBVLUwIfwDshajAovR30eekBQDp8MNV2OghCuLyYC8rinI65oCAHB5fDhWUuf/Ka3DiTIHvKqGWocH3zsq8f3RysY1DyItyYKctDjkWOOQnRaLbGsc4mOM7dcYanUMEKJWJoR/bCYx1oTEWMDtVeF0+eATGjTN3+ceroeWWkyGwMmLgL/LrsTegOKyepwoq8fxsnrYa1wQACpqXKiocQUG5gEgIcaILGssMlNjkZkSg8yUWKSnxITtID01DwOEqI2ZjQosjYPSQgBeVYPQAA0CHo8Kr6pC1RCWeyqKLCM3PR656fHAVf77XB4f6twafiiqgK2iAScrHKiocUIIoM7pRd2JGvx0oiZoO8nxJqQnxyA9OQZpyRakJcUgLcmCxDgT5HDt76OzMECI2sHpexxGRQYav5A3dXlpQsDlUdHg9sHn08J2DwXw76XkZMUhPeHngXWvT0NpVQNK7P6ju0qrGlBS6YTD6T+Sqbreg+p6z1nBYlRkpCaakZpogTXRgpREM1ITzEhJtCAl3twisw9Ty2GAEIWZpjGUGJMBsWYDfD7Nv0fSeNiP16fB7fFB1cK368tokAOHDp/O4fKirMqJsionKqqdKK9xorzaheo6NwT8e2OlVU6UVjnPud34GCOS401ITjAjOc6MpHgTkuLNSIozISnOhPgYI2ROg9JmGCBEYUwIQFFkKKcNGViMChJijPBpGjw+/5FePlULnJcSrqECAHEWIy7LNvqP+jqN16ehss4Fe40L9lr/v1V1blTWulFd74ba2J1X7/Si3unFifJzX5xLkvzjLglxJiTEmJAQa2z88YdLQqwR8TFGxMUYYTLILTYFfkfFACGKUAZZhsEkI87sfxtrQkBVBbyqBo9Hg9en+kMF4R0qgH+PJTMlFpkpsWf9TdMEahs8qKpzo6rOHyjV9R7U1LtR4/Cgpt4Dd+PFsIQAahu8qG3wArjwFSCNioxYiwFxMUbEWQyIsxgRY/Hv9cWe9m+M2QBVluFx+2A2KRyjOQ0DhCjCNYWDBAkGRYJBkRFj8p/n1nQND00TUBunXfd5fw4XLdyTBf4zvpPjzUiON+Oy7HMv4/L4UOvwotbhQW2DB3UNHtQ6vKhzelDX4EV9g/93j1cLrONVNX8AOTwh1yIBMJv8B0HEWYy4c1h3XH1F+iW2MHK1WYAUFRVh5syZqK6uRnJyMhYtWoRu3boFLaOqKhYsWIDNmzdDkiRMnToV99xzT1uVSBRV/PM5+ScO3Hu0AusKi1FZ60aWNQajrumMXl1S8P9v+gnfF1XDq2nQVIFeXZLQu1sqtu2zwV7rRmKcEdf3yUavLik4VFyFzXtOoarOjZQEM27sn4OT5fXYsrcEbp8Ks0HBDf2yMHxQ53PWc671z7fdpsOGQ93Gtn02fx1eFWajv447hnY7a32PT8W+f9mx/fsS1Dg8iDUb0CUzHi6vimMldXB7VMiyBLNRgU8TcLl9OP1gOAHA5VHh8qiorvdg+ecHsSHjOMYM6YJ+PdIu/UmLMJI418QvreDBBx/E3XffjfHjx+Ozzz7DRx99hBUrVgQt8+mnn2LNmjX405/+hOrqakyYMAHvvfceOnXqFPLj2O31YXn4Y0tIT09AeXlde5fRaqK5fe3Ztr2HK/Dulz9CUWSYDDI8Pg2qqsGaYMYPx2v8U6UrMgyNU6ebDRKsyf5zNVR/Hxj6X27Ftz9WBK7f4fZqqK13wenRgMbL2Xp9ApoQuHlAzlkhcqi4Cqu3FkFRZBgVGV7VX8Ognun45sfys+4fN/Sys0LkfNvomhGPPUcqIcFfW2PJGHF1bkh1OF0+QJIQY1bOqqFn52S4vSpMFhNKyurgdPtw+GQNdvxQCpNBQUZKDGobvFBVDfff0jMiQ0SWJVit8Rdf8BzaZA/EbrfjwIEDWL58OQBg7NixmD9/PiorK5GamhpYbu3atbjnnnsgyzJSU1MxcuRIrFu3DpMnTw75saL9CAy2L3K1V9sKD5QiIzUWJsPPI/Een4ryKicyUmKCZm1p+uoVH+M/JNfYuOzX35ciMd4Es8ngn/JckiBJElIkCUaD//wQqXH9ynoPrMkx/u6xxoH9E+UO9OySAqMiQ9P8E1p6fBpO2hvQLScJBkUOFOD2+vDT8Wr065F2+tSXOFhchey0OBgCRxQIeLwqSqqcyEyJCfr/1TSBQ8drMHpI16D/i++LKs/6v7DX+I/4sp52KWCPT8X3RZXo090Kk1FBcnIszI0l7vlXBTpnJiA+xoQYs4K4GB88Pv/lhQdEYHfWpbwu2yRAbDYbMjMzoTQ+8YqiICMjAzabLShAbDYbcnJyArezs7NRUlLSrMdKSYlrmaLDlN5vCpEimtvXXm177pFr2+VxT/erSYMueRu/6X7pH85zeujfRpY17pK3EW14Vg4REenSJgGSnZ2N0tJSqKr/UDtVVVFWVobs7Oyzljt16lTgts1mQ1ZWVluUSEREzdQmAWK1WpGXl4eCggIAQEFBAfLy8oK6rwBgzJgxWLVqFTRNQ2VlJTZs2IDRo0e3RYlERNRMbXYU1uHDhzFz5kzU1tYiMTERixYtQvfu3TFlyhQ89dRT6Nu3L1RVxQsvvICtW7cCAKZMmYKJEye2RXlERNRMbRYgREQUXTiITkREujBAiIhIFwYIERHpwgAhIiJdInY23ieeeAInTpyALMuIjY3Fb37zG+Tl5YU0aWOkWLJkCV577TWsWbMGPXv2xHfffYfZs2fD7XYjNzcXv/3tb2G1Wtu7TF2GDx8Ok8kEs9kMAJgxYwZuvPHGqGij2+3GSy+9hO3bt8NsNmPAgAGYP39+VLw2T5w4gWnTpgVu19XVob6+Hjt27IiK9gHAP/7xDyxevLhxJmMNTz75JEaNGhU17fvnP/+JxYsXw+fzISkpCQsXLkTnzp31tU9EqNra2sDvX375pZgwYYIQQogHHnhAfPrpp0IIIT799FPxwAMPtEt9l2r//v3i0UcfFTfffLM4dOiQ0DRNjBw5UuzcuVMIIcTrr78uZs6c2c5V6veLX/xCHDp0KOi+aGnj/PnzxYsvvig0TRNCCFFeXi6EiJ7X5ukWLFgg5s2bJ4SIjvZpmiby8/MDr82DBw+KAQMGCFVVo6J91dXVYvDgweLIkSNCCH87HnnkESGEvucvYgPkdJ988om48847RUVFhRg0aJDw+XxCCCF8Pp8YNGiQsNvt7Vxh87jdbnHvvfeK4uLiwAftnj17xO233x5Yxm63iwEDBrRjlZfmXAESDW2sr68XgwYNEvX19UH3R8tr83Rut1sMGTJE7N+/P2rap2maGDx4sNi1a5cQQogdO3aIUaNGRU379uzZI2677bbA7aqqKtGzZ0/d7YvYLiwAmDVrFrZu3QohBP785z+HPGljuFu8eDHGjRuHzp1/nor6zIkmU1NToWlaYHczEs2YMQNCCAwaNAj/9V//FRVtPH78OJKTk7FkyRIUFhYiLi4OTz/9NCwWS1S8Nk+3adMmZGZm4qqrrsL+/fujon2SJOHVV1/FE088gdjYWDgcDvzxj3+Mms+Wyy67DBUVFdi7dy/69euHNWvWAAh9wtszRfQg+osvvoh//vOfmD59Ol5++eX2LqdF7N69G/v27cOkSZPau5RW9e6772L16tX46KOPIITACy+80N4ltQifz4fjx4+jd+/e+PjjjzFjxgw8+eSTaGhoaO/SWtxHH32Eu+++u73LaFE+nw9//OMf8cYbb+Af//gH3nzzTUyfPj1qnr+EhAS88sorWLhwIe666y7Y7XYkJibqbl9EB0iTCRMmoLCwEFlZWSFN2hjOdu7ciSNHjmDEiBEYPnw4SkpK8Oijj+LYsWNBE01WVlZCkqSI+WZ+pqbnxGQyYdKkSfj222/PmkwzEtuYk5MDg8GAsWPHAgD69++PlJQUWCyWiH9tnq60tBQ7d+7EHXfcASD0CVPD3cGDB1FWVoZBg/zTzw8aNAgxMTEwm81R0T4AuP766/H+++/j448/xr//+7/D5XIhNzdXV/siMkAcDgdsNlvg9qZNm5CUlBTypI3hbOrUqdiyZQs2bdqETZs2ISsrC8uWLcPkyZPhcrmwa9cuAMAHH3yAW2+9tZ2r1aehoQF1df6r8wkhsHbtWuTl5aFPnz4R38bU1FQMGTIkMJ9bUVER7HY7unXrFvGvzdN98sknGDZsGFJS/FcNjIb3HgBkZWWhpKQER44cAeCfw6+iogJdu3aNivYBQHl5OQBA0zT8/ve/x3333Yfc3Fxd7YvIubAqKirwxBNPwOl0QpZlJCUl4dlnn8VVV1113kkbI9Xw4cOxdOlS9OzZE99++y3mzJkTdIhrWlrkXULz+PHjePLJJ6GqKjRNQ48ePfD8888jIyMjKtp4/PhxPPfcc6iurobBYMCvfvUrDBs2LKpem6NHj8asWbNw0003Be6LlvatXr0af/rTnyBJ/iv1PfXUUxg5cmTUtG/WrFn49ttv4fV6MXToUDz33HMwm8262heRAUJERO0vIruwiIio/TFAiIhIFwYIERHpwgAhIiJdGCBERKQLA4SIiHSJ6LmwiE43cODAwO9OpxMmkykwt8+8efMwbty49irtkg0dOhSLFy9Gfn5+e5dCFMAAoaixe/fuwO/Dhw/HggULcP3117djRaHx+XwwGFr3rdgWj0EdD7uwqMNQVRWvv/46RowYgSFDhuCZZ55BbW0tAP9Z1L1798aqVatw4403YsiQIfjwww+xe/dujB07Fvn5+Vi4cGFgW++//z4efPBBzJ49G1dffTVuu+027Ny5M/D36upq/Pd//zeGDh2KYcOGYcmSJdA0LWjdefPm4ZprrsFbb72Fw4cP44EHHsDgwYNx7bXX4tlnn0V9fT0A/5nQdrsdjz76KAYOHIgVK1bgq6++wi233BLUvqFDhwamgfnf//1fPPPMM/jVr36FgQMH4vPPP79g+4n0YIBQh7Fs2TJs3boV7733Hr766isYjcagUFBVFYcOHcLGjRvx0ksvYcGCBXj77bfx17/+FatXr8bHH3+MPXv2BJbftWsXevXqhcLCQkydOhXTpk0LfOjPmDEDCQkJ2LBhA1atWoUNGzbgs88+C1o3Ly8PX3/9NR555BEA/qtsbtmyBWvWrEFRURGWLl0KAPjDH/4Aq9WKZcuWYffu3XjwwQdDau/69etx55134ptvvsHo0aMv2n6i5mKAUIfxwQcf4JlnnkFmZibMZjOmTZuGtWvX4vTZfKZNmwaTyYQRI0YAAMaPH4+UlBTk5ORg4MCBOHDgQGDZrKws3H///TAajZgwYQIyMzOxefNmnDx5Ert27cLMmTMRExODjIwMPPDAA/j8888D63bu3Bn33nsvFEWBxWJBjx49cN1118FkMiE9PR0PPfRQ0B6NHoMHD8awYcMgyzIsFktI7SdqDnaKUocghEBJSQmmTp0amCQP8M9IWlVVBcB/EZ2m2WUBwGw2B03kaLFYgq6bkJWVFfQYubm5KCsrw6lTp+B2u3HdddcFPU7Xrl3Pu25paSlefPFF7N69Gw6HA0IIpKenX1KbT3+Mi7U/EmeVpfbHAKEOQZIkZGZm4rXXXkOfPn3O+ntTiDRHSUlJ0O1Tp04hIyMDWVlZiI2Nxc6dO4M+rM+s53Qvv/wyYmNjUVBQgKSkJHz++ed49dVXz7t8bGwsnE5n4LbX60VNTc15H+Ni7SfSg11Y1GHcd999+N3vfhe4lozdbsemTZt0b6+kpATvv/8+fD4fPvvsM9hsNtxwww3o3LkzBgwYgJdffhn19fXQNA1Hjx4NDHCfi8PhQGxsLOLj43Hq1CksX7486O9WqxUnTpwI3O7evTtqamqwfft2eL1evPbaa4FB+rZqPxEDhDqMyZMn47rrrsNDDz2EgQMH4r777gsa02iu/Px8HDx4EIMHD8bSpUuxZMkSJCQkAAB+97vfoa6uDrfeeisGDx6M6dOnw263n3dbTz/9NL755hvk5+fjP//zPzFq1Kigvz/++ON45ZVXkJ+fj5UrVyI1NRWzZs3CjBkzMGzYMKSlpQV1v7VF+4l4PRAiHd5//32sX78ef/nLX9q7FKJ2wz0QIiLShQFCRES6sAuLiIh04R4IERHpwgAhIiJdGCBERKQLA4SIiHRhgBARkS4MECIi0uX/AQOnm/nkzh6fAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.set(color_codes=True)\n", "plt.xlim(30,90)\n", "plt.ylim(0,1)\n", "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
InterceptTemperatureFrequency
0130.00.834373
1130.50.826230
2131.00.817774
3131.50.809002
4132.00.799911
5132.50.790500
6133.00.780766
7133.50.770712
8134.00.760339
9134.50.749648
10135.00.738645
11135.50.727334
12136.00.715721
13136.50.703816
14137.00.691626
15137.50.679164
16138.00.666441
17138.50.653471
18139.00.640269
19139.50.626851
20140.00.613235
21140.50.599439
22141.00.585485
23141.50.571391
24142.00.557181
25142.50.542876
26143.00.528501
27143.50.514078
28144.00.499631
29144.50.485186
............
91175.50.025508
92176.00.024110
93176.50.022787
94177.00.021535
95177.50.020350
96178.00.019229
97178.50.018169
98179.00.017166
99179.50.016217
100180.00.015321
101180.50.014473
102181.00.013671
103181.50.012913
104182.00.012197
105182.50.011520
106183.00.010880
107183.50.010275
108184.00.009703
109184.50.009163
110185.00.008653
111185.50.008171
112186.00.007716
113186.50.007286
114187.00.006879
115187.50.006496
116188.00.006133
117188.50.005791
118189.00.005467
119189.50.005162
120190.00.004873
\n", "

121 rows × 3 columns

\n", "
" ], "text/plain": [ " Intercept Temperature Frequency\n", "0 1 30.0 0.834373\n", "1 1 30.5 0.826230\n", "2 1 31.0 0.817774\n", "3 1 31.5 0.809002\n", "4 1 32.0 0.799911\n", "5 1 32.5 0.790500\n", "6 1 33.0 0.780766\n", "7 1 33.5 0.770712\n", "8 1 34.0 0.760339\n", "9 1 34.5 0.749648\n", "10 1 35.0 0.738645\n", "11 1 35.5 0.727334\n", "12 1 36.0 0.715721\n", "13 1 36.5 0.703816\n", "14 1 37.0 0.691626\n", "15 1 37.5 0.679164\n", "16 1 38.0 0.666441\n", "17 1 38.5 0.653471\n", "18 1 39.0 0.640269\n", "19 1 39.5 0.626851\n", "20 1 40.0 0.613235\n", "21 1 40.5 0.599439\n", "22 1 41.0 0.585485\n", "23 1 41.5 0.571391\n", "24 1 42.0 0.557181\n", "25 1 42.5 0.542876\n", "26 1 43.0 0.528501\n", "27 1 43.5 0.514078\n", "28 1 44.0 0.499631\n", "29 1 44.5 0.485186\n", ".. ... ... ...\n", "91 1 75.5 0.025508\n", "92 1 76.0 0.024110\n", "93 1 76.5 0.022787\n", "94 1 77.0 0.021535\n", "95 1 77.5 0.020350\n", "96 1 78.0 0.019229\n", "97 1 78.5 0.018169\n", "98 1 79.0 0.017166\n", "99 1 79.5 0.016217\n", "100 1 80.0 0.015321\n", "101 1 80.5 0.014473\n", "102 1 81.0 0.013671\n", "103 1 81.5 0.012913\n", "104 1 82.0 0.012197\n", "105 1 82.5 0.011520\n", "106 1 83.0 0.010880\n", "107 1 83.5 0.010275\n", "108 1 84.0 0.009703\n", "109 1 84.5 0.009163\n", "110 1 85.0 0.008653\n", "111 1 85.5 0.008171\n", "112 1 86.0 0.007716\n", "113 1 86.5 0.007286\n", "114 1 87.0 0.006879\n", "115 1 87.5 0.006496\n", "116 1 88.0 0.006133\n", "117 1 88.5 0.005791\n", "118 1 89.0 0.005467\n", "119 1 89.5 0.005162\n", "120 1 90.0 0.004873\n", "\n", "[121 rows x 3 columns]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_pred" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 1.0\n", "1 1.0\n", "2 1.0\n", "3 1.0\n", "4 1.0\n", " ... \n", "116 1.0\n", "117 1.0\n", "118 1.0\n", "119 1.0\n", "120 1.0\n", "Length: 121, dtype: float64" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n", "logmodel.predict(data_pred)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Hide code", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }