{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" ] }, { "cell_type": "markdown", "metadata": { "hidePrompt": false }, "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": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "System Info\n", "\t 3.6.4 |Anaconda, Inc.| (default, Mar 13 2018, 01:15:57) \n", "[GCC 7.2.0]\n", "\t uname_result(system='Linux', node='c0de44a43350', release='4.4.0-164-generic', version='#192-Ubuntu SMP Fri Sep 13 12:02:50 UTC 2019', machine='x86_64', processor='x86_64')\n", "Imported modules\n", "\t IPython 7.12.0\n", "\t IPython.core.release 7.12.0\n", "\t PIL 7.0.0\n", "\t PIL.Image 7.0.0\n", "\t PIL._version 7.0.0\n", "\t _csv 1.0\n", "\t _ctypes 1.1.0\n", "\t _curses b'2.2'\n", "\t decimal 1.70\n", "\t argparse 1.1\n", "\t backcall 0.1.0\n", "\t cffi 1.13.2\n", "\t csv 1.0\n", "\t ctypes 1.1.0\n", "\t cycler 0.10.0\n", "\t dateutil 2.8.1\n", "\t decimal 1.70\n", "\t decorator 4.4.1\n", "\t distutils 3.6.4\n", "\t ipaddress 1.0\n", "\t ipykernel 5.1.4\n", "\t ipykernel._version 5.1.4\n", "\t ipython_genutils 0.2.0\n", "\t ipython_genutils._version 0.2.0\n", "\t ipywidgets 7.2.1\n", "\t ipywidgets._version 7.2.1\n", "\t jedi 0.16.0\n", "\t json 2.0.9\n", "\t jupyter_client 6.0.0\n", "\t jupyter_client._version 6.0.0\n", "\t jupyter_core 4.6.3\n", "\t jupyter_core.version 4.6.3\n", "\t kiwisolver 1.1.0\n", "\t logging 0.5.1.2\n", "\t matplotlib 2.2.3\n", "\t matplotlib.backends.backend_agg 2.2.3\n", "\t matplotlib.pylab 1.15.2\n", "\t numexpr 2.6.9\n", "\t numpy 1.15.2\n", "\t numpy.core 1.15.2\n", "\t numpy.core.multiarray 3.1\n", "\t numpy.lib 1.15.2\n", "\t numpy.linalg._umath_linalg b'0.1.5'\n", "\t numpy.matlib 1.15.2\n", "\t optparse 1.5.3\n", "\t pandas 0.22.0\n", "\t _libjson 1.33\n", "\t parso 0.6.0\n", "\t patsy 0.5.1\n", "\t patsy.version 0.5.1\n", "\t pexpect 4.8.0\n", "\t pickleshare 0.7.5\n", "\t platform 1.0.8\n", "\t prompt_toolkit 3.0.3\n", "\t ptyprocess 0.6.0\n", "\t pygments 2.5.2\n", "\t pyparsing 2.4.6\n", "\t pytz 2019.3\n", "\t re 2.2.1\n", "\t scipy 1.1.0\n", "\t scipy._lib.decorator 4.0.5\n", "\t scipy._lib.six 1.2.0\n", "\t scipy.fftpack._fftpack b'$Revision: $'\n", "\t scipy.fftpack.convolve b'$Revision: $'\n", "\t scipy.integrate._dop b'$Revision: $'\n", "\t scipy.integrate._ode $Id$\n", "\t scipy.integrate._odepack 1.9 \n", "\t scipy.integrate._quadpack 1.13 \n", "\t scipy.integrate.lsoda b'$Revision: $'\n", "\t scipy.integrate.vode b'$Revision: $'\n", "\t scipy.interpolate._fitpack 1.7 \n", "\t scipy.interpolate.dfitpack b'$Revision: $'\n", "\t scipy.linalg 0.4.9\n", "\t scipy.linalg._fblas b'$Revision: $'\n", "\t scipy.linalg._flapack b'$Revision: $'\n", "\t scipy.linalg._flinalg b'$Revision: $'\n", "\t scipy.ndimage 2.0\n", "\t scipy.optimize._cobyla b'$Revision: $'\n", "\t scipy.optimize._lbfgsb b'$Revision: $'\n", "\t scipy.optimize._minpack 1.10 \n", "\t scipy.optimize._nnls b'$Revision: $'\n", "\t scipy.optimize._slsqp b'$Revision: $'\n", "\t scipy.optimize.minpack2 b'$Revision: $'\n", "\t scipy.signal.spline 0.2\n", "\t scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", "\t scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", "\t scipy.special.specfun b'$Revision: $'\n", "\t scipy.stats.mvn b'$Revision: $'\n", "\t scipy.stats.statlib b'$Revision: $'\n", "\t seaborn 0.8.1\n", "\t seaborn.external.husl 2.1.0\n", "\t seaborn.external.six 1.10.0\n", "\t six 1.14.0\n", "\t statsmodels 0.9.0\n", "\t statsmodels.__init__ 0.9.0\n", "\t traitlets 4.3.3\n", "\t traitlets._version 4.3.3\n", "\t urllib.request 3.6\n", "\t zlib 1.0\n", "\t zmq 17.1.2\n", "\t zmq.sugar 17.1.2\n", "\t zmq.sugar.version 17.1.2\n" ] } ], "source": [ "def print_imported_modules():\n", " import sys\n", " print(\"Imported modules\")\n", " for name, val in sorted(sys.modules.items()):\n", " if(hasattr(val, '__version__')): \n", " print(\"\\t\",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(\"System Info\")\n", " print(\"\\t\",sys.version)\n", " print(\"\\t\",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": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
56/18/82672500
68/30/836731000
711/28/836701000
82/03/846572001
94/06/846632001
108/30/846702001
1110/05/846782000
1211/08/846672000
131/24/856532002
144/12/856672000
154/29/856752000
166/17/856702000
177/2903/856812000
188/27/856762000
1910/03/856792000
2010/30/856752002
2111/26/856762000
221/12/866582001
\n", "
" ], "text/plain": [ " Date Count Temperature Pressure Malfunction\n", "0 4/12/81 6 66 50 0\n", "1 11/12/81 6 70 50 1\n", "2 3/22/82 6 69 50 0\n", "3 11/11/82 6 68 50 0\n", "4 4/04/83 6 67 50 0\n", "5 6/18/82 6 72 50 0\n", "6 8/30/83 6 73 100 0\n", "7 11/28/83 6 70 100 0\n", "8 2/03/84 6 57 200 1\n", "9 4/06/84 6 63 200 1\n", "10 8/30/84 6 70 200 1\n", "11 10/05/84 6 78 200 0\n", "12 11/08/84 6 67 200 0\n", "13 1/24/85 6 53 200 2\n", "14 4/12/85 6 67 200 0\n", "15 4/29/85 6 75 200 0\n", "16 6/17/85 6 70 200 0\n", "17 7/2903/85 6 81 200 0\n", "18 8/27/85 6 76 200 0\n", "19 10/03/85 6 79 200 0\n", "20 10/30/85 6 75 200 2\n", "21 11/26/85 6 76 200 0\n", "22 1/12/86 6 58 200 1" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = pd.read_csv(\"./shuttle.csv\")\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAF9JJREFUeJzt3X2UXXV97/H3d5IACYlAg02VQAFJsVyBCOFJtDfx6Qa7JPUCBbyCl940ZUlul9y2htvVa6m1a1V8qHpFY+SiQldNVRBom14e1Ii0IASM4UHBuYBhEhogBshASGYy3/vH2bN7Mkxmzhlmz5lzeL/WmpWz9/mdne939pz5zN5nn9+JzESSJICuVhcgSZo8DAVJUslQkCSVDAVJUslQkCSVDAVJUqmyUIiIqyPiqYh4YC/3R0R8PiK6I2JDRJxQVS2SpMZUeaTwNWDxCPefAcwrvpYBX6qwFklSAyoLhcy8HfjlCEOWANdkzV3AgRHxuqrqkSSNbmoL/+9DgCfqlnuKdU8OHRgRy6gdTTB9+vQTDz300AkpsFEDAwN0dXXmyzOd2pt9tZ9O7W2i+nrkkUeeyczXjjaulaEQw6wbds6NzFwFrAJYsGBBrlu3rsq6mrZ27VoWLlzY6jIq0am92Vf76dTeJqqviPhFI+NaGbs9QP2f/HOBzS2qRZJEa0PhJuDC4iqkU4HnMvNlp44kSROnstNHEfENYCFwcET0AH8OTAPIzJXAGuA9QDfwInBRVbVIkhpTWShk5vmj3J/AJVX9/5Kk5nXeS/mSpDEzFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklSqNBQiYnFEPBwR3RFx2TD3HxAR/xARP4mIByPioirrkSSNrLJQiIgpwJXAGcAxwPkRccyQYZcAD2Xm8cBC4NMRsU9VNUmSRlblkcLJQHdmPpqZu4DVwJIhYxKYFREBzAR+CfRXWJMkaQSRmdVsOOJsYHFmLi2WLwBOyczldWNmATcBbwRmAedm5j8Ns61lwDKAOXPmnLh69epKah6r3t5eZs6c2eoyKtGpvdlX++nU3iaqr0WLFt2bmQtGGze1whpimHVDE+g/AeuBtwNvAG6NiB9m5vN7PChzFbAKYMGCBblw4cLxr/YVWLt2LZOtpvHSqb3ZV/vp1N4mW19Vnj7qAQ6tW54LbB4y5iLg+qzpBh6jdtQgSWqBKkPhHmBeRBxRvHh8HrVTRfU2Au8AiIg5wNHAoxXWJEkaQWWnjzKzPyKWAzcDU4CrM/PBiLi4uH8l8JfA1yLifmqnm1Zk5jNV1SRJGlmVrymQmWuANUPWray7vRl4d5U1SJIa5zuaJUklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVDIUJEklQ0GSVKo0FCJicUQ8HBHdEXHZXsYsjIj1EfFgRPygynokSSOb2sigiHhTZj7QzIYjYgpwJfAuoAe4JyJuysyH6sYcCHwRWJyZGyPiV5v5PyRJ46vRI4WVEXF3RHyo+EXeiJOB7sx8NDN3AauBJUPGvB+4PjM3AmTmUw1uW5JUgcjMxgZGzAN+DzgHuBv4ambeOsL4s6kdASwtli8ATsnM5XVjPgtMA/4DMAv4XGZeM8y2lgHLAObMmXPi6tWrG+tugvT29jJz5sxWl1GJTu3NvtpPp/Y2UX0tWrTo3sxcMOrAzGz4C5gCnAVsAn4K/Az4z3sZew5wVd3yBcD/HjLmC8BdwP7AwcDPgd8YqYYTTzwxJ5vvf//7rS6hMp3am321n07tbaL6AtZlA7/nG31N4TjgIuC3gVuB92bmfRHxeuBO4PphHtYDHFq3PBfYPMyYZzLzBeCFiLgdOB54pJG6JEnjq9HXFL4A3Accn5mXZOZ9AJm5GfizvTzmHmBeRBwREfsA5wE3DRlzI/C2iJgaETOAU6gdgUiSWqChIwXgPcCOzNwNEBFdwH6Z+WJmXjvcAzKzPyKWAzdTO+10dWY+GBEXF/evzMyfRsT/BTYAA9RONzV1lZMkafw0Ggq3Ae8EeovlGcAtwFtGelBmrgHWDFm3csjyJ4FPNliHJKlCjZ4+2i8zBwOB4vaMakqSJLVKo6HwQkScMLgQEScCO6opSZLUKo2ePvow8K2IGLx66HXAudWUJElqlYZCITPviYg3AkcDAfwsM/sqrUySNOEaPVIAOAk4vHjMmyOCHObdx5Kk9tXom9euBd4ArAd2F6sTMBQkqYM0eqSwADimeKu0JKlDNXr10QPAr1VZiCSp9Ro9UjgYeCgi7gZ2Dq7MzDMrqUqS1BKNhsLlVRYhSZocGr0k9QcR8evAvMy8rZi8bkq1pUmSJlpDrylExO8D3wa+XKw6BLihqqIkSa3R6AvNlwCnA88DZObPAT9PWZI6TKOhsDNrn7MMQERMpfY+BUlSB2k0FH4QEX8KTI+IdwHfAv6hurIkSa3QaChcBjwN3A/8AbXPSNjbJ65JktpUo1cfDQBfKb4kSR2q0bmPHmOY1xAy88hxr0iS1DLNzH00aD/gHOBXxr8cSVIrNfSaQmZurfvalJmfBd5ecW2SpAnW6OmjE+oWu6gdOcyqpCJJUss0evro03W3+4HHgd8d92okSS3V6NVHi6ouRJLUeo2ePvofI92fmZ8Zn3IkSa3UzNVHJwE3FcvvBW4HnqiiKElSazTzITsnZOZ2gIi4HPhWZi6tqjBJ0sRrdJqLw4Bddcu7gMPHvRpJUks1eqRwLXB3RHyH2jub3wdcU1lVkqSWaPTqo7+KiH8G3lasuigzf1xdWZKkVmj09BHADOD5zPwc0BMRR1RUkySpRRr9OM4/B1YA/7NYNQ3426qKkiS1RqNHCu8DzgReAMjMzTjNhSR1nEZDYVdmJsX02RGxf3UlSZJapdFQ+GZEfBk4MCJ+H7gNP3BHkjpOo1cffar4bObngaOBj2bmrZVWJkmacKMeKUTElIi4LTNvzcw/ycw/bjQQImJxRDwcEd0RcdkI406KiN0RcXYzxUuSxteooZCZu4EXI+KAZjYcEVOAK4EzgGOA8yPimL2M+wRwczPblySNv0bf0fwScH9E3EpxBRJAZv7hCI85GejOzEcBImI1sAR4aMi4/w5cR23CPUlSCzUaCv9UfDXjEPacRbUHOKV+QEQcQu1y17czQihExDJgGcCcOXNYu3Ztk6VUq7e3d9LVNF46tTf7aj+d2ttk62vEUIiIwzJzY2Z+fQzbjmHW5ZDlzwIrMnN3xHDDiwdlrgJWASxYsCAXLlw4hnKqs3btWiZbTeOlU3uzr/bTqb1Ntr5Ge03hhsEbEXFdk9vuAQ6tW54LbB4yZgGwOiIeB84GvhgRv9Pk/yNJGiejnT6q//P9yCa3fQ8wr5gjaRNwHvD++gGZWc6fFBFfA/4xM29AktQSo4VC7uX2qDKzPyKWU7uqaApwdWY+GBEXF/evbKpSSVLlRguF4yPieWpHDNOL2xTLmZmvGenBmbkGWDNk3bBhkJn/taGKJUmVGTEUMnPKRBUiSWq9Zj5PQZLU4QwFSVLJUJAklQwFSVLpVRMKW3t38pMnnmVr785WlyKpSVt7d7Kjb7fP3wnwqgiFG9dv4vRPfI8PXPUjTv/E97hp/aZWlySpQYPP38eefsHn7wTo+FDY2ruTFddt4KW+Abbv7OelvgE+ct0G/+KQ2kD983d3ps/fCdDxodCzbQfTuvZsc1pXFz3bdrSoIkmN8vk78To+FOYeNJ2+gYE91vUNDDD3oOktqkhSo3z+TryOD4XZM/flirOOY79pXczadyr7TeviirOOY/bMfVtdmqRR1D9/p0T4/J0AjX7ITls7c/4hnH7UwfRs28Hcg6b7AyW1kcHn79133sG/nPlWn78Ve1WEAtT+4vCHSWpPs2fuy/RpU3wOT4COP30kSWqcoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqSSoSBJKhkKkqRSpaEQEYsj4uGI6I6Iy4a5/79ExIbi618j4vgq65EkjayyUIiIKcCVwBnAMcD5EXHMkGGPAf8xM48D/hJYVVU9kqTRVXmkcDLQnZmPZuYuYDWwpH5AZv5rZm4rFu8C5lZYjyRpFJGZ1Ww44mxgcWYuLZYvAE7JzOV7Gf/HwBsHxw+5bxmwDGDOnDknrl69upKax6q3t5eZM2e2uoxKdGpv9tV+OrW3iepr0aJF92bmgtHGTa2whhhm3bAJFBGLgP8GvHW4+zNzFcWppQULFuTChQvHqcTxsXbtWiZbTeOlU3uzr/bTqb1Ntr6qDIUe4NC65bnA5qGDIuI44CrgjMzcWmE9kqRRVPmawj3AvIg4IiL2Ac4DbqofEBGHAdcDF2TmIxXWIklqQGVHCpnZHxHLgZuBKcDVmflgRFxc3L8S+CgwG/hiRAD0N3LOS5JUjSpPH5GZa4A1Q9atrLu9FHjZC8uCrb076dm2g7kHTWf2zH3HbWw76dS+qtK9ZTvbXuyje8t2jpozq9XlqE1VGgoamxvXb2LFdRuY1tVF38AAV5x1HGfOP+QVj20nndpXVT56w/1cc9dG/ujYfi79m9u58LTD+NiSY1tdltqQ01xMMlt7d7Liug281DfA9p39vNQ3wEeu28DW3p2vaGw76dS+qtK9ZTvX3LVxj3XX3LmR7i3bW1SR2pmhMMn0bNvBtK49d8u0ri56tu14RWPbSaf2VZX1Tzzb1HppJIbCJDP3oOn0DQzssa5vYIC5B01/RWPbSaf2VZX5hx7Y1HppJIbCJDN75r5ccdZx7Deti1n7TmW/aV1ccdZxw77Q2szYdtKpfVXlqDmzuPC0w/ZYd+Fph/lis8bEF5onoTPnH8LpRx3c0JU3zYxtJ53aV1U+tuRYLjz1cO6/9y5uu/RUA0FjZihMUrNn7tvwL8JmxraTTu2rKkfNmUXPjGkGgl4RTx9JkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpZChIkkqGgiSpVGkoRMTiiHg4Iroj4rJh7o+I+Hxx/4aIOKHKeqRmbe3dyU+eeJatvTtHHbvusa185paHWffY1nHbZjNju7dsZ9uLfXRv2T7q2GZUVW+zNezo2z3qdru3bOfb657o2O9BFdsdampVG46IKcCVwLuAHuCeiLgpMx+qG3YGMK/4OgX4UvGv1HI3rt/Eius2MK2ri76BAa446zjOnH/IsGM/cNVd3NFdC4PPf6+btx01m2uXnvqKttnM2I/ecD/X3LWRPzq2n0v/5nYuPO0wPrbk2DF2Xn29Y6nhD3+zj0s/8b29bnfwezCoE78H473d4VR5pHAy0J2Zj2bmLmA1sGTImCXANVlzF3BgRLyuwpqkhmzt3cmK6zbwUt8A23f281LfAB+5bsOwf6Wte2xrGQiDfti99WVHDM1ss5mx3Vu27/HLEOCaOze+4r+Wq6p3rDXsztzrdl8t34Px3O7eRGZWs+GIs4HFmbm0WL4AOCUzl9eN+UfgrzPzjmL5u8CKzFw3ZFvLgGXF4tHAw5UUPXYHA8+0uoiKdGpvI/YV06bPmHrQ634jurqmDK7LgYHd/duefCT7drxYP3bKrINfP2X/A1/2x8zuF559cvf2ZzaPZZvNjO2accDsqa957eEAu198jikzDgCg//mnHx948bmRz2WNoKp6x1rDYG/Dbbf+e1CvTb4H4/azOIpfz8zXjjaostNHQAyzbmgCNTKGzFwFrBqPoqoQEesyc0Gr66hCp/bWyX31P/dUx/UFndvbZPtZrPL0UQ9waN3yXGDzGMZIkiZIlaFwDzAvIo6IiH2A84Cbhoy5CbiwuArpVOC5zHyywpokSSOo7PRRZvZHxHLgZmAKcHVmPhgRFxf3rwTWAO8BuoEXgYuqqqdik/bU1jjo1N7sq/10am+Tqq/KXmiWJLUf39EsSSoZCpKkkqEwBhHxeETcHxHrI2Jdse7yiNhUrFsfEe9pdZ3NiogDI+LbEfGziPhpRJwWEb8SEbdGxM+Lfw9qdZ3N2ktfnbC/jq6rf31EPB8RH273fTZCX52wzy6NiAcj4oGI+EZE7DfZ9pevKYxBRDwOLMjMZ+rWXQ70ZuanWlXXKxURXwd+mJlXFVeMzQD+FPhlZv51MX/VQZm5oqWFNmkvfX2YNt9f9YppZTZRmybmEtp8nw0a0tdFtPE+i4hDgDuAYzJzR0R8k9rFNscwifaXRwoCICJeA/wW8H8AMnNXZj5LbSqSrxfDvg78TmsqHJsR+uo07wD+X2b+gjbfZ0PU99UJpgLTI2IqtT9ONjPJ9pehMDYJ3BIR9xZTcAxaXsz2enWrDwHH4EjgaeCrEfHjiLgqIvYH5gy+d6T491dbWeQY7K0vaO/9NdR5wDeK2+2+z+rV9wVtvM8ycxPwKWAj8CS192XdwiTbX4bC2JyemSdQm+X1koj4LWozvL4BmE9th3+6hfWNxVTgBOBLmflm4AXgZdOdt6G99dXu+6tUnBI7E/hWq2sZT8P01db7rAixJcARwOuB/SPiA62t6uUMhTHIzM3Fv08B3wFOzswtmbk7MweAr1CbJbad9AA9mfmjYvnb1H6Zbhmcubb496kW1TdWw/bVAfur3hnAfZm5pVhu9302aI++OmCfvRN4LDOfzsw+4HrgLUyy/WUoNCki9o+IWYO3gXcDD8SeU36/D3igFfWNVWb+G/BERBxdrHoH8BC1qUg+WKz7IHBjC8obs7311e77a4jz2fMUS1vvszp79NUB+2wjcGpEzIiIoPaz+FMm2f7y6qMmRcSR1I4OoHZq4u8y868i4lpqh7UJPA78QbvN4xQR84GrgH2AR6ld7dEFfBM4jNoP9TmZ+cuWFTkGe+nr87T5/gKIiBnAE8CRmflcsW427b/PhuurE55jfwGcC/QDPwaWAjOZRPvLUJAklTx9JEkqGQqSpJKhIEkqGQqSpJKhIEkqVfbJa9JEKy7F/G6x+GvAbmpTXEDtDYa7WlLYCCLi94A1xfsppJbzklR1pMk0a21ETMnM3Xu57w5geWaub2J7UzOzf9wKlOp4+kivChHxwYi4u5iH/4sR0RURUyPi2Yj4ZETcFxE3R8QpEfGDiHh0cL7+iFgaEd8p7n84Iv6swe1+PCLuBk6OiL+IiHuKefRXRs251N6M9ffF4/eJiJ6IOLDY9qkRcVtx++MR8eWIuJXa5H5TI+Izxf+9ISKWTvx3VZ3IUFDHi4g3UZsW4S2ZOZ/aadPzirsPAG4pJjjcBVxObfqBc4CP1W3m5OIxJwDvj4j5DWz3vsw8OTPvBD6XmScBxxb3Lc7MvwfWA+dm5vwGTm+9GXhvZl4ALAOeysyTgZOoTcx42Fi+P1I9X1PQq8E7qf3iXFebcobp1KZQANiRmbcWt++nNp1xf0TcDxxet42bM3MbQETcALyV2vNnb9vdxb9PhwLwjoj4E2A/4GDgXuCfm+zjxsx8qbj9buA3I6I+hOZRmyZBGjNDQa8GAVydmf9rj5W1Dzqp/+t8ANhZd7v++TH0xbccZbs7snjBrpjH5wvUZmfdFBEfpxYOw+nn34/gh455YUhPH8rM7yKNI08f6dXgNuB3I+JgqF2lNIZTLe+O2mc9z6A2J/6/NLHd6dRC5pliht2z6u7bDsyqW34cOLG4XT9uqJuBDxUBNPi5xtOb7El6GY8U1PEy8/5idsrbIqIL6AMupvZRiI26A/g7ah/ycu3g1UKNbDczt0btc6IfAH4B/Kju7q8CV0XEDmqvW1wOfCUi/g24e4R6vkxtVs31xamrp6iFlfSKeEmqNIriyp43ZeaHW12LVDVPH0mSSh4pSJJKHilIkkqGgiSpZChIkkqGgiSpZChIkkr/HzHofwgP0tIHAAAAAElFTkSuQmCC\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": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
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: Mon, 06 Apr 2020 Deviance: 3.0144
Time: 11:51:34 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: Mon, 06 Apr 2020 Deviance: 3.0144\n", "Time: 11:51:34 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": 4, "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": 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: -23.526
Date: Mon, 06 Apr 2020 Deviance: 18.086
Time: 11:51:34 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: Mon, 06 Apr 2020 Deviance: 18.086\n", "Time: 11:51:34 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": 5, "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": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3dyb7QmLYISA7yA5hEXEBrYK2KiriinVBpHWp7SNVn199tE+16oNt1VZxQ3GpgisupYJa44JbQBBkX8UEkJ0kkD33748ZMGAgQzLJLPm8rivXzDlzn3O+dwY+c3LmnPuYcw4REYkunlAXICIiwadwFxGJQgp3EZEopHAXEYlCCncRkSikcBcRiUI1hruZPW1mW83s28O8bmb2sJmtMbPFZjYw+GWKiMjRCGTPfTow+givnwl09f9MBKbWvSwREamLGsPdOfcxsPMITc4FnnM+XwDpZtY6WAWKiMjRiwnCOtoC31eZzvXP23xoQzObiG/vnsTExKx27drVaoOVlZV4PNHxdYH6Ep6ipS/R0g9QX/ZbtWrVdudc85raBSPcrZp51Y5p4Jx7AngCYNCgQW7+/Pm12mB2djYjRoyo1bLhRn0JT9HSl2jpB6gv+5nZd4G0C8bHYC5QdRc8E9gUhPWKiEgtBSPc3wKu8J81czywxzn3k0MyIiLScGo8LGNmLwEjgGZmlgvcCcQCOOceA2YDZwFrgH3AVfVVrIiIBKbGcHfOXVLD6w64PmgViUhEKCsrIzc3l+Li4gbZXlpaGsuXL2+QbdW3QPqSkJBAZmYmsbGxtdpGML5QFZFGKDc3l9TUVDp06IBZdedVBFdBQQGpqan1vp2GUFNfnHPs2LGD3NxcOnbsWKttRMd5RSLS4IqLi2natGmDBHtjY2Y0bdq0Tn8VKdxFpNYU7PWnrr9bhbuISBTSMXcRiVher5c+ffocmJ41axYdOnQIXUFhROEuIhErMTGRRYsWHfb18vJyYmIaZ8zpsIyIRJXp06dz4YUXcvbZZ3PGGWcAMGXKFAYPHkzfvn258847D7S955576N69Oz/72c+45JJLeOCBBwAYMWIE+4dH2b59+4G/BioqKpg8efKBdT3++OPAj8MJjB07lh49enDZZZfhO0sccnJyOOGEE+jXrx9DhgyhoKCAUaNGHfShNHz4cBYvXhzU30Pj/EgTkaD649tLWbYpP6jr7NmmCXee3euIbYqKiujfvz8AHTt25I033gDg888/Z/HixWRkZDB37lxWr17NV199hXOOc845h48//pjk5GRmzJjBwoULKS8vZ+DAgWRlZR1xe9OmTSMtLY2cnBxKSkoYPnz4gQ+QhQsXsnTpUtq0acPw4cOZN28eQ4YM4aKLLmLmzJkMHjyY/Px8EhMTueKKK5g+fToPPvggq1atoqSkhL59+wbht/YjhbuIRKzDHZY5/fTTycjIAGDu3LnMnTuXAQMGAFBYWMjq1aspKCjgvPPOIykpCYBzzjmnxu3NnTuXxYsX8+qrrwKwZ88eVq9eTVxcHEOGDCEzMxOA/v37s2HDBtLS0mjdujWDBw8GoEmTJgCcd955DB8+nClTpvD0009z5ZVX1u0XUQ2Fu4jUWU172A0tOTn5wHPnHLfffjvXXXfdQW0efPDBw55uGBMTQ2VlJcBB55o75/j73//OqFGjDmqfnZ1NfHz8gWmv10t5eTnOuWq3kZSUxOmnn86bb77Jyy+/TG1HyD0SHXMXkag2atQonn76aQoLCwHIy8tj69atnHzyybzxxhsUFRVRUFDA22+/fWCZDh06sGDBAoADe+n71zV16lTKysoAWLVqFXv37j3stnv06MGmTZvIyckBfFemlpeXAzBhwgRuuukmBg8efOCvjGDSnruIRLUzzjiD5cuXM2zYMABSUlJ44YUXGDhwIBdddBH9+/fn2GOP5aSTTjqwzC233MK4ceN4/vnnOfXUUw/MnzBhAhs2bGDgwIE452jevDmzZs067Lbj4uKYOXMmN954I0VFRSQmJvL+++8DkJWVRZMmTbjqqnoaa9E5F5KfrKwsV1sffvhhrZcNN+pLeIqWvtRnP5YtW1Zv665Ofn5+va7/zjvvdFOmTKnXbeyXn5/v8vLyXNeuXV1FRcVh21X3OwbmuwAyVodlREQa2IsvvsjQoUO555576u3WgTosIyIC3HXXXQ22rUsvvfQnX/AGm/bcRaTWnKv2dskSBHX93SrcRaRWEhIS2LFjhwK+Hjj/eO4JCQm1XocOy4hIrWRmZpKbm8u2bdsaZHvFxcV1CrtwEkhf9t+JqbYU7iJSK7GxsbW+S1BtZGdnH7jKNNI1RF90WEZEJAop3EVEopDCXUQkCincRUSikMJdRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSiUMSF+77Sct7bUEZ5RWWoSxERCVsRF+7vLN7MP1eUMu7xz/lux+HvXSgi0phFXLiPG9SOSX3jWbO1kDMf+oSZORs15KiIyCEiLtwBjm8Tw7s3n0z/dunc+toSbnhxIXuKykJdlohI2IjIcAdok57IC9cM5dbRPZizdAtnPfQJX2/cFeqyRETCQsSGO4DHY/xqRGde/dUJeDww7rHPefLjdTpMIyKNXkDhbmajzWylma0xs9uqeT3NzN42s2/MbKmZXRX8Ug+vf7t03rnxJE47rgX3zF7OxOcX6DCNiDRqNYa7mXmBR4AzgZ7AJWbW85Bm1wPLnHP9gBHAX8wsLsi1HlFaYiyPXZ7F//yiJx+u2Mq5//iUFVvyG7IEEZGwEcie+xBgjXNunXOuFJgBnHtIGwekmpkBKcBOoDyolQbAzLj6xI7MmHg8+0orGPPIPN76ZlNDlyEiEnJW0/FpMxsLjHbOTfBPjweGOuduqNImFXgL6AGkAhc55/5VzbomAhMBWrZsmTVjxoxaFV1YWEhKSsoR2+wuqeTRRSWs2lXJWR1jGdstFo9ZrbZXnwLpS6RQX8JPtPQD1Jf9Ro4cucA5N6jGhs65I/4AFwJPVZkeD/z9kDZjgb8BBnQB1gNNjrTerKwsV1sffvhhQO1Kyirc/3tjsTv21nfcFdO+dHuKSmu9zfoSaF8igfoSfqKlH86pL/sB810Nue2cC+iwTC7Qrsp0JnDosY6rgNf9217jD/ceAay7XsXFeLh7TB/uPb8P89Zs5/xHP9NVrSLSKAQS7jlAVzPr6P+S9GJ8h2Cq2gicBmBmLYHuwLpgFloXlwxpz/PXDGV7YQnnPjKPL9btCHVJIiL1qsZwd86VAzcAc4DlwMvOuaVmNsnMJvmb/Qk4wcyWAB8AtzrnttdX0bUxrHNT3rx+OE2T4xg/7UveWJgb6pJEROpNTCCNnHOzgdmHzHusyvNNwBnBLS34jm2azOu/Gs6kFxbw25nfsHFHETed1gULwy9aRUTqIqKvUK2NtKRYnr16COcPbMvf3l/Fba8t0fDBIhJ1AtpzjzZxMR7+cmE/MtMTefg/a9hWWMI/Lh1AUlyj/HWISBRqdHvu+5kZvzujO/ec15vslVu59Mkv2bW3NNRliYgERaMN9/0uG3osUy/PYtnmfC58/HM27S4KdUkiInXW6MMdYFSvVjx71RC27Clm7NTPWLetMNQliYjUicLdb1jnpsyYeDzF5ZWMe/xzlm/WoGMiErkU7lX0bpvGy9cNI8bj4aLHP2ehbv4hIhFK4X6ILi1SeGXSMNKT4hg/7StyNuwMdUkiIkdN4V6NdhlJvHzdMFo0ieeKaV/x2ZqwuthWRKRGCvfDaJWWwMyJw2ifkcRV03P4eNW2UJckIhIwhfsRNE+N56WJx9OxWTITnpuvgBeRiKFwr0FGchwvXns8nZunMOG5+XykgBeRCKBwD0BGchwvThhK5+YpTHxuPvN0DF5EwpzCPUDHJMfxzwlD6dA0mWuezdGY8CIS1hTuRyEjOY5/XjuUzGOSuHp6Dgu+02mSIhKeFO5HqVlKPC9eO5SWTRK48ukcvs3bE+qSRER+QuFeCy1SE/jnhKE0SYxl/LQvWbmlINQliYgcROFeS23SE3nx2qHExXi4fNqXuvG2iIQVhXsdHNs0mReuGUp5RSWXPfUlW/YUh7okERFA4V5nXVum8uzVQ9i9r4zLp+mGHyISHhTuQdA3M50nrxjExp37uHJ6DntLykNdkog0cgr3IBnWuSn/uGQAS3J3M+mFBZSW66bbIhI6CvcgOqNXK+47vy+frN7Of73yDZWVLtQliUgjFRPqAqLNuMHt2LG3lPvfXUHzlHju+MVxmFmoyxKRRkbhXg8mndKJrQXFPD1vPS2bxHPdKZ1DXZKINDIK93pgZtzx855sKyjh3n+voEWTeM4bkBnqskSkEVG41xOPx/jLuH7sKCzl968upkVqAsO7NAt1WSLSSOgL1XoUH+PlsfFZdGqWwqTnF7B8c36oSxKRRkLhXs/SEmN55qrBJMfHcNUzOWzeUxTqkkSkEVC4N4A26Yk8c9VgCkvKueqZHAqKy0JdkohEOYV7AzmudRMevWwgq7cWcv2LCymr0EVOIlJ/FO4N6ORuzfnzeb35eNU2/ufNpTini5xEpH7obJkGdtHg9ny3Yx+PZq+lY7MkuoW6IBGJStpzD4FbzujOz/u05t5/r2D+Fg0yJiLBF1C4m9loM1tpZmvM7LbDtBlhZovMbKmZfRTcMqPL/nPg+7dL54nFJSzO3R3qkkQkytQY7mbmBR4BzgR6ApeYWc9D2qQDjwLnOOd6ARfWQ61RJSHWyxPjB5EaZ0x4dr5OkRSRoApkz30IsMY5t845VwrMAM49pM2lwOvOuY0AzrmtwS0zOjVPjee3WQnsK63gmunzNQ68iASN1XTGhpmNBUY75yb4p8cDQ51zN1Rp8yAQC/QCUoGHnHPPVbOuicBEgJYtW2bNmDGjVkUXFhaSkpJSq2XDTWFhIeuKEvjbghL6t/By44B4PBE6imS0vS/R0Jdo6QeoL/uNHDlygXNuUE3tAjlbprqkOfQTIQbIAk4DEoHPzewL59yqgxZy7gngCYBBgwa5ESNGBLD5n8rOzqa2y4ab7OxsbvrFCJq0Wc9dby8jp6Q1t47uEeqyaiXa3pdo6Eu09APUl6MVSLjnAu2qTGcCm6pps905txfYa2YfA/2AVUhAfnlCB1ZvLWRq9lq6NE/hgiyNIikitRfIMfccoKuZdTSzOOBi4K1D2rwJnGRmMWaWBAwFlge31OhmZtx1Ti9O6NyU219fwoLvdoa6JBGJYDWGu3OuHLgBmIMvsF92zi01s0lmNsnfZjnwLrAY+Ap4yjn3bf2VHZ1ivR4evWwgrdMTuO75BeTt1hk0IlI7AZ3n7pyb7Zzr5pzr7Jy7xz/vMefcY1XaTHHO9XTO9XbOPVhfBUe79KQ4pv1yECVllUx4VmfQiEjt6ArVMNSlRSoPXzqAlVvy+a+XdaNtETl6CvcwNbJ7C/77rON4d+kWHv7P6lCXIyIRRgOHhbFrTuzIii0FPPj+arq3TOXMPq1DXZKIRAjtuYcxM+Oe83ozsH06v3v5G5Zt0m36RCQwCvcwt/8+rGmJsVz73Hx2FJaEuiQRiQAK9wjQIjWBJ67IYnthCb/+59e6i5OI1EjhHiH6ZqZz/wV9+XL9Tv749tJQlyMiYU5fqEaQMQPasnxzPo9/vI6erdO4dGj7UJckImFKe+4R5veje3BKt+bc+da3zN+gIQpEpHoK9wjj9RgPXzyAzGOSmPTCAjZpiAIRqYbCPQKlJcXy5BVZFJdVct3zCyguqwh1SSISZhTuEapLi1T+dlF/luTt4fbXl1DTTVdEpHFRuEew03u25Hend+ONhXlM+3R9qMsRkTCicI9wN4zswuherfjz7OV8snpbqMsRkTChcI9wHo/xwLh+dGmRwg0vLmTjjn2hLklEwoDCPQqkxMfw5BW+++Ve+5zGgBcRhXvUOLZpMv+4dACrtxZwyyvf6AtWkUZO4R5FTuranNvO7MG/v93CIx+uCXU5IhJCCvcoc+1JnRjTvw1/eW8V7y/7IdTliEiIKNyjjJlx3wV96dWmCTfPXMSarYWhLklEQkDhHoUSYr08Pn4Q8TEeJj43nz1FZaEuSUQamMI9SrVNT2Tq5Vls3LmP38xYSIVusi3SqCjco9iQjhnceU4vsldu44G5K0Ndjog0II3nHuUuH9qeZZvymZq9lp6tm3B2vzahLklEGoD23KOcmfHHc3ox6NhjmPzqN3ybtyfUJYlIA1C4NwJxMR6mXp7FMUlxTHxuPtt1k22RqKdwbySap8bz5BWD2LG3lF+9sIDSct1kWySaKdwbkd5t0/i/sX3J2bCLO99aqiEKRKKYvlBtZM7t35blmwt47KO19GydyvhhHUJdkojUA+25N0KTR3Xn1B4t+OPby/hs7fZQlyMi9UDh3gh5PcZDF/enQ7Nkrv/n1xoDXiQKKdwbqdSEWJ66YhCVDiY8l0NBsYYoEIkmCvdGrEOzZB69bCBrt+3l5hmLNESBSBRRuDdyw7s0486ze/LBiq1MmaMhCkSihc6WEcYffywrt/jOoOnWMoXzB2aGuiQRqaOA9tzNbLSZrTSzNWZ22xHaDTazCjMbG7wSpb6ZGXed04thnZpy22tLWPDdzlCXJCJ1VGO4m5kXeAQ4E+gJXGJmPQ/T7n5gTrCLlPoX6/Uw9fKBtElPYOJzC8jdpTNoRCJZIHvuQ4A1zrl1zrlSYAZwbjXtbgReA7YGsT5pQOlJcUy7cjClFZVMeHa+zqARiWBW0yXo/kMso51zE/zT44GhzrkbqrRpC7wInApMA95xzr1azbomAhMBWrZsmTVjxoxaFV1YWEhKSkqtlg034diXpdsr+MuCYno383LzwHg8ZgEtF459qa1o6Uu09APUl/1Gjhy5wDk3qKZ2gXyhWt3/7EM/ER4EbnXOVdgRgsA59wTwBMCgQYPciBEjAtj8T2VnZ1PbZcNNOPZlBJCW+R1/mPUtnxS24M6zewW0XDj2pbaipS/R0g9QX45WIOGeC7SrMp0JbDqkzSBghj/YmwFnmVm5c25WUKqUBnf58ceybttenp63nk7NkjUGjUiECSTcc4CuZtYRyAMuBi6t2sA513H/czObju+wjII9wv2/nx/Hdzv2cudbS8nMSGJk9xahLklEAlTjF6rOuXLgBnxnwSwHXnbOLTWzSWY2qb4LlNDxeoyHLxnAca2bcMM/v2bZpvxQlyQiAQroPHfn3GznXDfnXGfn3D3+eY855x6rpu2V1X2ZKpEpOT6Gab8cTGpCLNc8m8OWPcWhLklEAqDhB6RGrdISmHblIPKLyrh6eg6FJeWhLklEaqBwl4D0apPGI5cNZOUPBdzw4teUV+g2fSLhTOEuARvRvQV/Orc32Su3cceb3+o2fSJhTAOHyVG5dGh7cnft49HstbRNT+SGU7uGuiQRqYbCXY7a5FHd2bKnmAfmrqJVWiJjszSKpEi4UbjLUTMz7rugL1sLSrjttcU0T43nlG7Na7WuWQvzmDJnJZt2F9EmPZHJo7ozZkDbIFcs9UXvX/jSMXeplbgY3yiS3Vqm8qsXFvDN97uPeh2zFuZx++tLyNtdhAPydhdx++tLmLUwL/gFS9Dp/QtvCneptdSEWKZfPZimKXFcPT2HLXuP7gyaKXNWUlRWcdC8orIK3REqQuj9C28Kd6mTFqkJPHvVEBzwwPxifsgP/CKnTbuLjmq+hBe9f+FN4S511ql5CtOvGkxhqeOKaV+xZ19g48C3SU88qvkSXvT+hTeFuwRF38x0bhqYwPrte7n62Rz2ldZ8FevkUd1JjPUeNC8x1svkUd3rq0wJIr1/4U3hLkHTs6mXBy/uz8KNu5j0wteUlh/5GPyYAW259/w+tE1PxIC26Ynce34fnW0RIfT+hTedCilBdVaf1tx3fl9+/9pifjtzEQ9fMgCv5/A3cBkzoK3CIILp/QtfCncJunGD25FfXMbd/1pOcryX+87vi+cIAS8iwadwl3ox4aROFBSX89AHq0mM9XLXOb040i0YRSS4FO5Sb27+WVf2lZbz5CfrSYjzctvoHgp4kQaicJd6Y2b891nHsa+0gsc/Wkd8jJffnd4t1GWJNAoKd6lXZsafzu1NaXklD3+wmliPceNpGklSpL4p3KXeeTy+gcYqKh1/eW8VMV4PvxrROdRliUQ1hbs0CK/HmHJhP8orHfe/uwKH49cjuoS6LJGopXCXBuP1GH8d1w8z+L93V+IcXD9SAS9SHxTu0qBivB7+Oq4/hm9UwfIKx02nddFZNCJBpnCXBuf1GH8Z1x+vx8Pf3l9FSXkFk0d1V8CLBJHCXULC6zGmjO1LXIyHR7PXUlxWyR2/OE4BLxIkCncJGY/H+PN5vYmP8fD0vPXsKy3nnvP6HHEsGhEJjMJdQsrMuPPsnqQmxPD3/6yhsKScv47rT1yMBiwVqQuFu4ScmfFfZ3QnJT6Ge/+9goLicqZePpCkOP3zFKkt7R5J2LjulM7cd34fPlm9jUuf/JJde0tDXZJIxFK4S1i5eEh7Hr0si2Wb87nw8c/J0/04RWpF4S5hZ3TvVjx71RB+yC/mvEfmsWxTfqhLEok4CncJS8M6N+WVScPwmDHu8c/5ZPW2UJckElEU7hK2erRqwhvXn0DmMYlc9UwOM3M2hrokkYihcJew1jotkVcmDWNY56bc+toS7n93BZWVLtRliYQ9hbuEvdSEWJ6+cjCXDGnP1Oy1/OqfC9hXWh7qskTCWkDhbmajzWylma0xs9uqef0yM1vs//nMzPoFv1RpzGK9Hv58Xm/u+EVP3lv2A2Onfs4mnUkjclg1hruZeYFHgDOBnsAlZtbzkGbrgVOcc32BPwFPBLtQETPjmhM7Mu3KwXy/cx9n//1Tvlq/M9RliYSlQPbchwBrnHPrnHOlwAzg3KoNnHOfOed2+Se/ADKDW6bIj0Z2b8Eb1w8nLTGWS5/8guc/34BzOg4vUpXV9J/CzMYCo51zE/zT44GhzrkbDtP+FqDH/vaHvDYRmAjQsmXLrBkzZtSq6MLCQlJSUmq1bLhRX2pvX5nj8cUlfLOtguFtYriiVxzx3uAMOhYt70u09APUl/1Gjhy5wDk3qMaGzrkj/gAXAk9VmR4P/P0wbUcCy4GmNa03KyvL1daHH35Y62XDjfpSNxUVle5v7610HW57x43620du/bbCoKw3Wt6XaOmHc+rLfsB8V0O+OucCOiyTC7SrMp0JbDq0kZn1BZ4CznXO7QhgvSJ15vEYN/+sG89cOZjNe4o5+++f8q/Fm0NdlkjIBRLuOUBXM+toZnHAxcBbVRuYWXvgdWC8c25V8MsUObIR3Vvwr5tOpHOLFK5/8WvumPUtxWUVoS5LJGRqDHfnXDlwAzAH3yGXl51zS81skplN8jf7H6Ap8KiZLTKz+fVWschhZB6TxMvXDePakzry/BffMeaReaz+oSDUZYmEREADZjvnZgOzD5n3WJXnE4CffIEq0tBmL9nM7CVbAFj1QwFnPfwJY/q3Zd6a7WzeU0yb9EQmj+rOmAFtg77tWQvzmDJnJZt2F9XrdgLxh1lLeOnL77m5dxnX3D6bS4a24+4xfUJSi4SG7oYgUWPWwjxuf30JRf7DMZUOXIXjlQW5B9rk7S7i9teXAAQ1eA/ddn1tJxB/mLWEF774cRyeCucOTCvgGw8NPyBRY8qclQfCdb/qTvQtKqtgypyV9b7t+thOIF768vujmi/RSeEuUeNohiMI9k1ADrftUAyRUHGYa1cON1+ik8Jdokab9MSA23oM3lyUF7QrWw+37aOpKVi8Vv2FXIebL9FJ4S5RY/Ko7iTGeg+aF+sxYg+5ajU+xkPmMUn8ZsYiJjw7n9xd++pl24mxXiaP6l7ndR+tS4a2O6r5Ep0U7hI1xgxoy73n96FteiIGtE1PZMqF/Zgytt9B8+6/oC8f3jKCP/z8OD5bu4PT//oxT3y8lrKKyqBu+97z+4TkbJm7x/Th8uPbH9hT95px+fHt9WVqI6OzZSSqjBnQttpArW7ehJM6Mbp3K+56ayl/nr2C1xbk8cdze3F8p6ZB3XYo3D2mD3eP6UN2djZrLxsR6nIkBLTnLo1a5jFJPPXLwTwxPovCknIufuILbnppITuLa78XLxIOtOcuApzRqxUndW3O1I/W8thHa3m3spLvvKu47pROJMXpv4lEHu25i/glxnn53end+OB3p9C/hZeHPljNyAeymZmzkQrdt1UijMJd5BDtMpL4df8EXpk0jNZpidz62hLOfOhj3lv2g24KIhFD4S5yGIM7ZPDGr0/g0csGUlbhuPa5+Yx59DM+Xb1dIS9hT+EucgRmxll9WvPeb0/m/gv6sC2/mMunfcm4xz9n3hqFvIQvhbtIAGK8Hi4a3J7/3DKC/z23F9/vLOKyp77kgqmf8Z8VOlwj4UfhLnIUEmK9XDGsA9mTR/Cnc3vxQ34JV0+fz1kPf8obC3PrdCGUSDAp3EVqISHWy3h/yD9wYT/KKir57cxvOPn/PuTxj9ayZ19ZqEuURk4n8IrUQazXw9isTM4f0JbsVVt54uN13PvvFTz4/mouyGrLFcM60K1laqjLlEZI4S4SBB6PcWqPlpzaoyVLN+1h+rwNvDw/lxe+2MjQjhlcfvyxnNGrJfEx3ppXJhIECneRIOvVJo0pF/bj9rOO4+X53/PCF99x40sLyUiO44KBbRk3qB1dtTcv9UzhLlJPMpLjmHRKZyae1IlP12znpa828sy8DTz5yXr6ZaYxNiuTn/dtQ0ZyXKhLlSikcBepZx6PcXK35pzcrTnbC0uYtTCPVxfkcsebS/nj28sY0b05Z/drw8+Oa0lyvP5LSnDoX5JIA2qWEs+Ekzox4aROLN+cz6yFeby5aBPvL99KQqyHU3u04MzerRnZowUpCnqpA/3rEQmR41o34bjWTbh1dA/mf7eLt7/ZxL+/3cLsJVuIi/FwUpdmnN6zJacd15LmqfGhLlcijMJdJMQ8HmNIxwyGdMzgrnN6seC7Xcxespn3lv3AByu2YraEvpnpnNajBSO6N6d3mzQ8Ht0PVY5M4S4SRrxVgv7Os3uyfHMBHyz3hfzf3l/FX99bRUZyHCd2acaJXZtxYpc9/m12AAANA0lEQVRmIbkJt4Q/hbtImDIzerZpQs82TbjxtK5sLyzh09XbyV65lU/XbOetbzYB0LFZMsd3yuD4Tk0Z0jGD1mkKe1G4i0SMZinxB+7T6pxj5Q8FfLp6O5+v3cE732zmpa++B6BdRiKDj81g4LHH4Aoqqah0eHUYp9FRuItEIDOjR6sm9GjVhAkndaK8opLlmwv4asNOvlq/g49Xb+P1hXkA3Jczh76Z6fRrl06/zDT6tkunTVoCZgr8aKZwF4kCMV4PfTLT6JOZxjUndsQ5x8ad+3jh3c8oSWnNwo27mfbpOsoqfEMTH5MUS++2afRs7Tvsc1zrJnRslkysV2MJRguFu0gUMjOObZrM8LaxjBjRG4CS8gpWbC5gce5ulm7K59tNe3hm3gZK/cMUx3k9dGqeTPdWqXRtkULXlql0aZFC+4wkhX4EUriLNBLxMV7foZl26QfmlVVUsm7bXpZvzmfFlgJWbMln/oZdvLlo04E2sV6jfUYSHZul0Kl5Msc2TaJj02TaN02idVqijueHKYW7SCMW6/XQvVUq3VsdPJBZYUk5a7YWsnZrIWu2FbJ+217Wb9/Lx6u3UVpeWWV5o216Iu0yksg8JonMYxJpm55I22MSaZOeSMvUeGK01x8SCncR+YmU+Bj6t0unf5W9fIDKSseW/GI2bN/Lhh37+H7XPjbu2Efurn3M3bSFHXtLD2rvMWiRmkCrtARaNfE9tmgST4vUBFqkxtOiSTzNUuLJSIrThVlBpnAXkYB5PEabdN9e+Qldfvr6vtJyNu0uIndXEZv3FLN5dxGb9hTzQ34xa7YVMm/tdgqKy3+ynNdjZCTH0TQ5jmYp8TRNiSMjOY6MpDgyUuI4JimO73ZU0HJzPsckxZGWGEtCrEdn/ByBwl1EgiYpLoYuLVLp0uLw49XvKy1na34JWwtK2FpQzPaCErYXlrKtoIQde0vZXljCxp372LW3lIKSgz8I7s/55MDzOK+HtKRYmiTEkJYYS5PEWJokxJKaEEPqgccYUuJ//En2//iee0mKi4na7wwCCnczGw08BHiBp5xz9x3yuvlfPwvYB1zpnPs6yLWKRK1ZC/OYMmclm3YX0SY9kcmjuvPK/I3MW7vzQJvhnTO4cFD7n7QDfjJv/nc7eenL77m5dxnX3D6bS4a24+4xfQLa7pgBbQ87P5Dl92+7wjm8Zj/ZdlJcDB2axbDo+9019uXOs3tycrfm7NxXyofzcujQrSe79pXx2drtZK/cxraCEgqKy/B6jLIKx4bte8kvLqeguOzAaZ81iY/xkBTnC/qkOC+JcV4SY70kxXlJiPU9T4jzkhDjJSHWQ0Lsj4/xMR7iY/yPsT8+j9v/4/3xebzXS2yMEev14FxgtdVFjeFuZl7gEeB0IBfIMbO3nHPLqjQ7E+jq/xkKTPU/ikgNZi3M4/bXl1BUVgFA3u4ibp656Cft5q3deVDY5+0uYvKr34CDskp3YN7vZi6isspyFc7xwhcbAQ4K2eq2e/vrS5j/3U5eW5D3k/nAQQFf3fJ12fbkV74B40Ao5+0u4o43l3Lv+X0YM6AtW5p6GdGnNbMW5vHB8q0Hli0uq+T7nUUH2u1XXFZBYUk5hcXlFBSXU1hSzt6ScvaWlrO3pIK9JeXsK61gb2k5+0p9z4tKKw487thbemC6pLyC4rJKisoqqKisezCf2TGWkSPrvJojCmTPfQiwxjm3DsDMZgDnAlXD/VzgOef7OPrCzNLNrLVzbnPQKxaJMlPmrDwQVEerur3TymraAbz05fcHBWx12y0qqziw133o/ClzVh4UntUtX5dtl1UTmoFut7p2vj1sL81SgjtccllFJcVlFZSUV1JS7nte6n/ue6ygpKyS0grf9IHH8krKKnw/nl0bg1pTdaymPw/MbCww2jk3wT89HhjqnLuhSpt3gPucc5/6pz8AbnXOzT9kXROBif7J7sDKWtbdDNhey2XDjfoSnhqsL3GtumTV17or9u3Bm5R2YLp0y5oFddluXZYPwrLNgO1HWrbqNsJcXf59Heuca15To0D23Kv7tuHQT4RA2uCcewJ4IoBtHrkgs/nOuUF1XU84UF/CU7T0xczml+/ZGvH9gOh5T6Bh+hLI1QW5QLsq05nAplq0ERGRBhJIuOcAXc2so5nFARcDbx3S5i3gCvM5Htij4+0iIqFT42EZ51y5md0AzMF3KuTTzrmlZjbJ//pjwGx8p0GuwXcq5FX1VzIQhEM7YUR9CU/R0pdo6QeoL0elxi9URUQk8mhEHxGRKKRwFxGJQmEf7maWYGZfmdk3ZrbUzP7on59hZu+Z2Wr/4zGhrjUQZuY1s4X+awMiuR8bzGyJmS0ys/n+eZHal3Qze9XMVpjZcjMbFol9MbPu/vdj/0++md0coX35rf//+7dm9pI/ByKuHwBm9ht/P5aa2c3+efXel7APd6AEONU51w/oD4z2n5FzG/CBc64r8IF/OhL8BlheZTpS+wEw0jnXv8r5upHal4eAd51zPYB++N6fiOuLc26l//3oD2ThO7nhDSKsL2bWFrgJGOSc643vRI6LibB+AJhZb+BafFf69wN+YWZdaYi+OOci5gdIAr7GN27NSqC1f35rYGWo6wug/kz/G3kq8I5/XsT1w1/rBqDZIfMiri9AE2A9/pMLIrkvh9R/BjAvEvsCtAW+BzLwndH3jr8/EdUPf50X4htscf/0HcDvG6IvkbDnvv9QxiJgK/Cec+5LoKXzn0vvf2wRyhoD9CC+N7bqEByR2A/wXYE818wW+IeVgMjsSydgG/CM/3DZU2aWTGT2paqLgZf8zyOqL865POABYCOwGd91M3OJsH74fQucbGZNzSwJ3ynj7WiAvkREuDvnKpzvT81MYIj/T52IYma/ALY65yJl7IuaDHfODcQ3Iuj1ZnZyqAuqpRhgIDDVOTcA2EsE/Ll/JP6LDc8BXgl1LbXhP/58LtARaAMkm9nloa2qdpxzy4H7gfeAd4FvgJ/eraQeRES47+ec2w1kA6OBH8ysNYD/cWsISwvEcOAcM9sAzABONbMXiLx+AOCc2+R/3IrvuO4QIrMvuUCu/69BgFfxhX0k9mW/M4GvnXM/+KcjrS8/A9Y757Y558qA14ETiLx+AOCcm+acG+icOxnYCaymAfoS9uFuZs3NLN3/PBHfG78C35AHv/Q3+yXwZmgqDIxz7nbnXKZzrgO+P5n/45y7nAjrB4CZJZtZ6v7n+I6HfksE9sU5twX43sy6+2edhm8464jrSxWX8OMhGYi8vmwEjjezJDMzfO/JciKvHwCYWQv/Y3vgfHzvTb33JeyvUDWzvsCz+L4x9wAvO+f+18yaAi8D7fH9Y7jQObfz8GsKH2Y2ArjFOfeLSOyHmXXCt7cOvsMaLzrn7onEvgCYWX/gKSAOWIdv+AwPkdmXJHxfRnZyzu3xz4u498V/yvNF+A5hLAQmAClEWD8AzOwToClQBvzOOfdBQ7wnYR/uIiJy9ML+sIyIiBw9hbuISBRSuIuIRCGFu4hIFFK4i4hEoUBukC3SoPyniX3gn2wFVOAbIgBgiHOuNCSFHYGZXQ3M9p83LxJyOhVSwpqZ3QUUOuceCINavM65isO89ilwg3Nu0VGsL8Y51yCXokvjo8MyElHM7JfmG99/kZk9amYeM4sxs91mNsXMvjazOWY21Mw+MrN1ZnaWf9kJZvaG//WVZvaHANd7t5l9hW9coz+aWY5/fO7HzOcifMNRz/QvH2dmuVWurD7ezN73P7/bzB43s/fwDVYWY2Z/9W97sZlNaPjfqkQjhbtEDP+AcecBJ/gHkovBN5QDQBow1z+YWSlwF77L1i8E/rfKaob4lxkIXGpm/QNY79fOuSHOuc+Bh5xzg4E+/tdGO+dmAouAi5xvPPWaDhsNAM52zo0HJuIbUG4IMBjfIGzta/P7EalKx9wlkvwMXwDO9w05QiK+S+0Bipxz7/mfL8E3TGy5mS0BOlRZxxzn3C4AM5sFnIjv/8Hh1lvKj0MtAJxmZpOBBKAZsAD491H2403nXLH/+RnAcWZW9cOkK75L0kVqTeEukcSAp51zdxw00ywGXwjvV4nvDl77n1f9d37ol0yuhvUWOf8XU/5xW/4BDHTO5ZnZ3fhCvjrl/PiX8aFt9h7Sp1875z5AJIh0WEYiyfvAODNrBr6zampxCOMM890zNQnfmOHzjmK9ifg+LLb7R8W8oMprBUBqlekN+G51xyHtDjUH+LX/g2T/fVATj7JPIj+hPXeJGM65Jf7RAt83Mw++UfYmAZuOYjWfAi8CnYHn95/dEsh6nXM7zOxZfMMbfwd8WeXlZ4CnzKwI33H9u4AnzWwL8NUR6nkc38iAi/yHhLbi+9ARqROdCimNhv9MlN7OuZtDXYtIfdNhGRGRKKQ9dxGRKKQ9dxGRKKRwFxGJQgp3EZEopHAXEYlCCncRkSj0/wHRUJwHFwSFegAAAABJRU5ErkJggg==\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": { "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": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/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": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8lOW9///XPWsmkz1kARLCYtgRLQiiIoqyyO5CFag7Vm2139Z6jvbUetQetZ5fjz217WlFa11LK7iwpKIVFERFsS6RQNgDYckkZJ1MZr3v+/fHJAMhAZKQySz5PB+PmMyde+5clyHznvtaFV3XdYQQQoiTGCJdACGEENFJAkIIIUS7JCCEEEK0SwJCCCFEuyQghBBCtEsCQgghRLvCFhA/+9nPmDRpEnPmzGn3+7qu81//9V9MmzaNuXPnUlJSEq6iCCGE6IKwBcQ111zD888/f8rvb9q0ibKyMt577z1++ctf8sgjj4SrKEIIIbogbAFxwQUXkJqaesrvr1+/ngULFqAoCueddx4NDQ1UVlaGqzhCCCE6KWJ9EA6Hg9zc3NDj3NxcHA7HGZ/n8wfw+VV8fhV/QMUf0AioGqqqoWo6MjFcCCG6hylSP7i9F3JFUc74vDqnD0eV84znKc3/URQFBVAUUFCCnxUFQ8v3FDAYlND3DAYFg6JgMIBBUTpUpu6UlZVMVQfqF6ukfrErnusGvaN+nRWxgMjNzaWioiL0uKKiguzs7G67vt78n9ZB1Pm7C0UJBkUwNIIfxuYQMRoUjMbg554OEiGECLeIBcTUqVN59dVXmT17Nt988w3JycndGhDdRddB1XVUdFBPfZ5BAaPBEAoMk9EQ+mwwSHgIIWJP2ALivvvu4/PPP6e2tpZLL72Ue++9l0AgAMCiRYuYMmUKGzduZNq0adhsNp544olwFaVHaDpoqoa/nRAxGBTMRgNmU/DDYjLIHYcQIuopsbbcd2VNU4f6IKKZAsGgMBuxmg2YTcbQ93pDO6jULzbFc92gd9SvsyLWxNSb6YAvoOELaDS6g81TVouJBIvxjM8VQoieIgERBTQd3N4Abm8A4zEXTS4fNqsJs0lWQhFCRI4ERJTRdJ0mb4AmbwCz0UBiQvDOQvoshBA9TQIiivlVjXqXD6dbwZ5gwmY1YZCgEEL0EGnDiAGapuNs8nOszk2j2y+zxYUQPUICIoZoOjS6/VTVe3B7A5EujhAizklAxCBN06l3+aiu9+APaJEujhAiTklAxDC/qlHT4JFmJyFEWEhAxDidYLNTTYNX7iaEEN1KAiJOyN2EEKK7SUDEkZa7ieoG6ZsQQpw9CYg4FFB1aho8NHn8kS6KECKGSUDEKR1oaPJT6/SiadLkJIToPAmIOOf1q1Q3eAio0uQkhOgcCYheQNV0apxeCQkhRKdIQPQSmqZT6/SiahISQoiOkYDoRVRNp7ZB+iSEEB0jAdHLBDSdukavzJUQQpyRBEQv5AtoNDTJEFghxOlJQPRSbm9A5kkIIU5LAqIXczb58fnVSBdDCBGlJCB6MR2oa5ROayFE+yQgejlNh3qXL9LFEEJEIQkIgdev0uSRHeqEEK1JQAgAnE0+mWkthGhFAkIAwf6I+kafzI8QQoRIQIgQv6rhkqYmIUQzCQjRisvtl6YmIQQgASFOohOcHyGEEDEXEI8+v4Wvdx+TtvIw8vpV3F5pahKitzNFugCdddDh5KDDyVe7q1gweRDpyQmRLlJccrr9WC1GDIoS6aIIISIk5u4gTMbgC9buQ/X874piPv72KJrcTXQ7TdNplKYmIXq1mAuIR5ZOYmBuMgD+gEbRpwd4bs12jtW7I1yy+NPkDeAPyFpNQvRWMRcQffvYWTp3JAsmD8JqNgJwoMLJ71Z+K3cTYVDvkrkRQvRWMRcQAAZFYcKIHP7fwnMpzEsFgmP4iz49wAtFO6hr9Ea4hPEjoOoyN0KIXiqsAbFp0yZmzJjBtGnTWLZsWZvvO51O7rrrLubNm8fs2bN54403OnX9tCQrt1w1nGsuHRy6m9h3pIFnVhbz9R4Z6dRdZG6EEL1T2AJCVVUee+wxnn/+eYqKili7di179uxpdc5rr73GkCFDWL16Na+88gpPPfUUPl/nVhZVFIXxw7P50XXnMrBvsG/C41N5fcMeXv9gDx6fvPs9WzI3QojeKWwBUVxcTEFBAfn5+VgsFmbPns369etbnaMoCi6XC13XcblcpKamYjJ1beRterKVpbNHctXEARgNwZFO3+yp5vdvfEt5ZeNZ16e38/pV2VxIiF4mbPMgHA4Hubm5occ5OTkUFxe3OmfJkiXcfffdTJ48GZfLxW9+8xsMhjNnVkaG/ZTfm395IeNG5vL86m1UVDdR4/Ty7OoSFkwZwpUTBsTEuP7T1S+SrGYjfdJsZ32drKzkbihN9Irn+sVz3SD+69dZYQuI9tr/lZNenDdv3syIESN4+eWXOXjwILfeeivjx48nKSnptNeuqXGd9vuJZgN3zRtF0acH2FpaiabpvPnBHrbvPcZ1l51DYkL0zg/MyLCfsX6R5HN7MZuMXX5+VlYyVVXObixRdInn+sVz3aB31K+zwtbElJubS0VFReixw+EgOzu71Tlvvvkm06dPR1EUCgoKyMvLY9++fd3y8y1mI1dfOphFVxaGOrBLD9bx+zeLKa+M338E4dbolj4dIXqLsAXEmDFjKCsro7y8HJ/PR1FREVOnTm11Tt++ffn0008BOHbsGPv37ycvL697yzE4k3uuHUO/zEQA6hp9LFu9nS0lFTLKqQu8fhV/QEY0CdEbhC0gTCYTDz/8MEuXLmXWrFlcddVVFBYWsnz5cpYvXw7AD37wA7766ivmzp3LLbfcwv33309GRka3lyUzJYE7549m4sgcAFRNZ/XHZaz8cK+82HVBo1tGNAnRGyh6jL2NrqxpwnEW7YRf7z7GW5v24W8e198vM5El04dGzaJ/0d4H0SIzJQGzqfPvL3pDO2+81i+e6wa9o36dFZMzqc/GeYV9uGvBKDKSrQAcqW7iD29uY9+R+giXLLY0eeQuQoh41+sCAqBvpp0fXjOGoflpQHBRuheKStmyveIMzxQtPD4VVZPmOSHiWa8MCACb1cRNM4Zx6dh+AGi6zurNZbz90T554esAHWSNJiHiXK8NCACDQWHmxAF8d+o5oX0mPt9RyYvvlMqOah3g9gZk9Vwh4livDogW553Th+/PG0VKohmAvYcb+OPb26hu8ES4ZNFN15EgFSKOSUA0y8tK4u6rj8+XOFbv4Y9vbaOsoiHCJYtuLk9A5pMIEackIE6Qarfw/XmjGDkwHWjpvN7Bt/uqI1yy6KVpOh6fLOInRDySgDiJxWxk8bShTD63LxDcMOdv7+9mc/FRead8Ci4Z8ipEXJKAaIdBUbjqwgLmXjwQRQmO2PnHlgMUfXpAOmXbEVB1WQpciDgkAXEak0blsmTaUMzG4P+mT7ZV8PqGPbK7WjtkyKsQ8UcC4gxGDsxg6dwRJFqDS4QX763m5XU78Uq7eytevyrBKUSckYDogPzsZO6cP4q0JAsAew7X89za7bJo3Uma5C5CiLgiAdFBWWk27pw/mpz04I5qR465WLa6hLpGb4RLFj3cPpk4J0Q8kYDohJZhsAW5wVURj9V7eHZVCcfq3BEuWXSQiXNCxBcJiE6yWU3cOmt4aKG/epePZ1eXcORY9C/R3RNk4pwQ8UMCogssJiPfmz6UMYMzgeCL4vNrt8tWpgQnznllyKsQcUECootMRgPXTz2HC4YH99n2+FT+XLSDfUdkaQ7prBYiPkhAnAWDQWHB5EFcPDoXAJ9f46V3StlVXhfhkkWWL6DJkFch4oAExFlSFIVZkwq47Pz+APhVjVfe3UnpwdoIlyyy5C5CiNgnAdENFEVh+gX5TL8gHwBV03ntvV3sKKuJcMkiR4a8ChH7JCC60WXn9+eqiQOA5pD4525K9vfOkNB18Hils1qIWCYB0c0mj+3H7EkFQHAb0+Xv72ZbL10uvMkrM82FiGUSEGFw8Zi+zL1oIBAMib+t39Mr7yRklVchYpsERJhMGp3LvIsHAsfvJLb3wj4JmVktROySgAijC0fltrqTWP7+7l7Xce3xqdJZLUSMkoAIs0mjc5lzUbBPQtV0/vr+7l41T0JHOquFiFUSED3gotF9Qx3Xqqbz6ns72XO4PsKl6jnSzCREbJKA6CEXj+nLjAnBeRIBVeeVd3ey/2jvWJbDr8rMaiFikQRED5pyXn+uGJcHgD+g8dK60l6zwF+T3EUIEXMkIHrY1O/0Z8p5/YDg2k1/+UcpR6vjf6lwj1eWARci1khA9LCWZTlaFvjz+FReKNpBZZxvOqTpyDLgQsQYCYgIaFngb3zzUuEuT4AXinZQ0+CJcMnCyy2jmYSIKRIQEaIoCgsuGcS5Q4KbDjW4fLxQtIP6ON7j2utXUaWzWoiYIQERQQaDwsLLhzCiIB2AGqeX3/79q7heKrvRLeszCRErOhQQt99+Ox988EGnOxk3bdrEjBkzmDZtGsuWLWv3nM8++4z58+cze/Zsvve973Xq+vHAaDBwwxWFDO6XAsCRKhcvvrMDry8+m2Pc0lktRMzoUEBcf/31vPTSS1x55ZUsW7aM2tozb4ajqiqPPfYYzz//PEVFRaxdu5Y9e/a0OqehoYFHH32UP/7xjxQVFfHb3/62a7WIcWaTgRunDyM/OwmAQ1UuXnlvJ/5A/DXHqJqOJ07DT4h406GAmD59Oi+++CLPPfcclZWVzJkzh3//939n27Ztp3xOcXExBQUF5OfnY7FYmD17NuvXr291zpo1a5g2bRr9+gWHfWZmZp5FVWKb1WLk5pnD6ZdlB2DfkQZe37AHTYu/d9sys1qI2GDqypPMZjNWq5UHHniAyZMn8+CDD7Y5x+FwkJubG3qck5NDcXFxq3PKysoIBALceOONuFwubrrpJhYsWHDGn5+RYe9KsaNeBvD/rj+f/+/Vf3Gszk1JWQ3vbC3nezOHoyhKpIvXbZJSbKSlJ2I2xWcXWFZWcqSLEDbxXDeI//p1VocC4r333uPVV1+lurqaxYsXU1RUhN1uJxAIMH369HYDor125pNf5FRVpaSkhBdffBGPx8MNN9zA2LFjGTRo0GnLU1MTvxPLMjLs3DxzGM+uKqHR7efjb45gBGY271QX6zIy7NTUuPC4vKTYLZEuTrfLykqmqio+Z8fHc92gd9SvszoUECtXruSOO+5g8uTJrZ9sMvHQQw+1+5zc3FwqKipCjx0OB9nZ2W3OSU9PJzExkcTERMaPH09paekZAyLeZaYkcOus4Ty3Zjsen8qmb45gt5mYfG6/SBet27h9AZISzRji6M5IiHjToXv8Z599tk04tJg6dWq7x8eMGUNZWRnl5eX4fD6KioranHvFFVfwxRdfEAgEcLvdFBcXM2TIkE5WIT71zbRz08xhmIzBF9B3thzky11VES5V99F14naklhDxokMBsXjxYurrjy9PXVdXx5IlS077HJPJxMMPP8zSpUuZNWsWV111FYWFhSxfvpzly5cDMGTIECZPnsy8efNYuHAh1113HUOHDj2L6sSXgbkpLLpyKIbmN9lvbtzLzoNnHkEWK6SzWojopugdGJQ+f/58Vq1adcZjPaGypglHHLcTtrTRn+hfOyt5Y+M+AMxGA7fPGcGAnNjsTDu5fn1SEzAZ46ezOp7bseO5btA76tdZHfrL1DSNpqam0GOXy4WqSvNATxk3LJuZE4Kd1H5V46V1O6msjY/F/eQuQojo1aGAmDNnDrfddhurVq1i1apV3H777cybNy/cZRMnmDy2L5eM6QsEX1T/8o/4WLfJLf0QQkStDo1iuvPOO8nOzmbDhg3ous4NN9zQofkKovsoisLMCwfQ6Pbz9Z5j1Lt8/OWdUu6cNwqbtUvTWaKCpul4/SpWszHSRRFCnKTDryxXX301V199dTjLIs7AoChcM2UwLo+f3Yfqqax18/K7O7lt1oiYnnTm9gYkIISIQh0KiOrqal555RXKy8sJBI63GffWtZMiyWQ0sPjKoTy/djuHj7k4UOHk7xt2s/jKoRgMsTmnwOtT0XRd5kQIEWU6FBD33nsvQ4YMYdKkSRiN8k4v0qwWIzdfNZw/rdpGTYOX7WW1rPmkjHkXD4zJJTl0wONVSUyI3aYyIeJRh/4iGxoa+OUvfxnusohOSLKZuXXWCP709jZcngCfbXeQardw2fn9I120LvH4AhIQQkSZDjVcFxYW4nA4wl0W0UmZKQncfNVwLM39D+9tLedfOysjXKqu8QU0VC3+ljcXIpZ1+A5i3rx5nH/++Vit1tBx6YOIvLysJBZPG8rL63ai6TpvbdpHks3MsAHpkS5ap7m9Kkm22O1sFyLedCgg5syZw5w5c8JdFtFFQ/PTuGbKYFZ+uBdNh+Xv72bp3JHkZSVFumid4vEFSLKZI10MIUSzDgWEDG+Nft8ZmkWDy8d7W8vxBYKzre+aP4rMlIRIF63DAqpOQNXiaukNIWJZh/4Sy8rKWLRoUWg11pKSEn73u9+FtWCi86ac148LR+YA4HL7efEfpTS6/REuVefIdqRCRI8OBcQjjzzC3XffTXJycLGnESNGsG7durAWTHSeoijMuWggIwcG+x+qGzy8vK4Unz92XnQ9sjaTEFGjQwHhdDq59NJLQ2PsDQYDZrO0FUcjg0Hh+qmFFDSv9nqoysXy9btRY2Rv64Cm4w/ETqAJEc86FBBGoxG/3x8KCIfDgcEg7cTRymwycOOMofRJDfY/7DxYx+rN+9vdBjYayQJ+QkSHDm8YdM8991BbW8vvfvc7Fi9ezG233RbusomzkJhg5tZZw0luHhW0tbSSDV8ejnCpOkb6IYSIDh0axbRgwQLy8vL44IMPcLvdPPXUU4wfPz7cZRNnKT05OJFu2ZoSfH6N9f86RKrdwvjh2Wd+cgRpmo7Xp2K1yLIuQkRSh9c2GD9+vIRCDOrXx86SaUN56Z3gRLq3P9pHcmL0T6Rr8gYkIISIsA4FxLXXXtvuInArV67s9gKJ7leYl8a1UwazIoYm0vn8Kpqmx+wKtULEgw4FxAMPPBD62uv1UlRURHZ2dDdTiNbOH5pFfQxNpNMBty+APUFGywkRKR0KiAkTJrR6fMkll0gndQyacl4/6l0+PtvuCE2ku3P+qKhd3sLtkYAQIpK6NFa1sbGR8vLy7i6LCDNFUZgbQxPpApoetWUTojfodB+EpmkcOnSIW2+9NawFE+HRMpHuz0XbOehoDE2k+970YRijsL3f7Q1gke1IhYiITvdBGI1G8vLyyMnJCVuhRHiZTQZumjGMZ1eXUFXnYefBOlZ9tI+rLx0cdTvSefwqybIdqRAR0aU+CBH7EhPM3HLVCP60ahvOJj9f7KwixW7hyvH5kS5aK7ou25EKESkd+qu78MIL231nqes6iqLw6aefdnvBRPilJ1u55arhLFu9Ha9fZcOXh0lOtDBxZHTdHbq9sh2pEJHQob+6RYsWUVdXx/XXX4+u67zxxhvk5OQwa9ascJdPhFnfTDvfmz6UF98pRdV0Vn+8n+REMyMHZkS6aCF+VcMf0DCbZP0vIXpSh/7itm7dyn/+538yfPhwRowYwUMPPcTGjRvp378//fv3D3cZRZgN6Z/KdZcNAYJNOn9bv5sDFc4Il6o1tywDLkSP61BAVFZWUlNTE3pcU1NDVVVV2Aolet7Yc/owe1IBENzZ7aV1pThqmiJcquM8vkDMrEYrRLzoUBPTzTffzPz587n88ssB2LhxI3feeWdYCyZ63sVj+tLg8vFR8VE8PpUX3wlOpEtLska6aGh6cJVXm1X6IoToKR36a1uyZAnjxo1j69at6LrOkiVLGDZsWLjLJiJgxsQBNLr9fLX7GPUuHy++U8r3546Kik5itzcgASFED+rwX1teXh6qqjJq1KhwlkdEmEFRuGbKYBrdfnYfqqey1s3L75Zy2+wRWEyRnbDmC2gEVA2TUTqrhegJHfpL27hxI7Nnz+bee+8F4Ntvv+Wuu+4Ka8FE5BgNBhZPG0pelh2Ag45Glr+/G1XTIlwy6awWoid1KCCeeeYZVq5cSUpKCgBjxozh4MGDYS2YiCyr2cjNVw1vtW3pW5v2Rbyj2O0NoElntRA9osP36llZWa0eWyyWbi+MiC72BDO3zhpBij34u/5y1zHe+exgRENC06HJI3cRQvSEDgWE3W7n2LFjodnUn332GcnJyWd83qZNm5gxYwbTpk1j2bJlpzyvuLiYESNGsG7dug4WW/SUltnWNmuw/2Fz8VE2fXMkomVyefxomtxFCBFuHQqIn/70p9xxxx0cOnSIG2+8kfvvv7/VAn7tUVWVxx57jOeff56ioiLWrl3Lnj172j3v17/+NZdccknXaiDCLjcjkZtmDMfc3Dn87uflbC2tjFh5dD0YEkKI8OrQKKaxY8fy8ssv8+WXXwJw/vnnh/ojTqW4uJiCggLy84OLv82ePZv169dzzjnntDrvlVdeYcaMGXz77bddKb/oIQW5ySyeVsgr7+4K7W1tsxgZPTgzIuVpat5MSLYkFSJ8zhgQqqry3e9+lzfeeIMpU6Z0+MIOh4Pc3NzQ45ycHIqLi9uc8/777/PSSy91KiAyMuwdPjcWRWv9JmXYMVlMvLC6BF2H1z/YQ59MOyMHdS4kuqt+1gQzacmRn8R3sqysMze/xqp4rhvEf/0664wBYTQaSU9Px+v1YrV2/I+xvY7Mk1eEffzxx7n//vsxGjs3vr6mxtWp82NJRoY9qus3JDeZORcPZM3HZQRUnT++Uczts0cwIKdjf1jdWb9awJOWgNEQPfMisrKSqaqKrnWsuks81w16R/06q0NNTAMHDmTJkiXMmDGDxMTE0PElS5ac8jm5ublUVFSEHjscDrKzs1uds23bNu677z4Aamtr2bhxIyaTiSuvvLJTlRA9a9KoXNzeAO9/cQh/QOPFd0q5Y+5I+mb27J2PDjQ2+UmNgqVAhIhHHQoIl8tFYWEh+/bt6/CFx4wZQ1lZGeXl5eTk5FBUVMT//M//tDpnw4YNoa8ffPBBLrvsMgmHGHH5+f3xeFU2fxtct+kv/yjl+3NH0ifN1qPlcPtUEmUpcCHC4rQB8atf/YoHH3yQJ598ko8//piLL7644xc2mXj44YdZunQpqqpy7bXXUlhYyPLly4HgHhMidimKwlUXDsDtC/CvnVU0uv38uWgH3583ivQe7hdodPt7/GcK0Rso+mlmPV199dW89dZbbb6OpMqaJhxx3E4Y7X0QJ9M0nb9v2M23+4LLwWemJHDHvJGkJLY/kTJc9ctItmIxR3atKIjvdux4rhv0jvp11mnvy0/MjkgvsSCik8GgsPDycxg2IA2A6gYPLxTt6PF5Cs4mmRchRHc7bUD4fD727t3Lnj17Wn3d8iEEgMloYPGVQxncLzg3prLWzV+KdvTownp+VZOF/IToZqdtYpo6deqpn6gorF+/PiyFOh1pYopeXr/KX/6xg4OORgDys5O4bdYIrJbjTT/hrJ/BoNAnNQGDErnJc/HcTBHPdYPeUb/OOm0n9YmjjIQ4E6vZyC1XDefPa3dw+JiL8spGXlpXyi1XDe+R/gFN03G6fDLsVYhuImMDRbdKsJi4ddZwcjOC82XKKpy8/O5OfAG1R36+26fi9fXMzxIi3sn+jaLbJSaYuW32CJ5bs52qOjf7jjTw6ru7uHHG2W1Tu/tQHV+UVlLr9JKebGX88GwK89LanFff5KOPObJNTaLrtu2vZnPxUarq3GSl2bjk3L6M7uRyLqJ7yB2ECIskm5nb54wIbTi053A9r763E38X7yR2H6rj3c/LqW7woulQ3eDl3c/L2X2ors25LU1NIvZs21/NGxv34ah1o+ngqHXzxsZ9bNtfHemi9UoSECJsUhItLJ0zkszmkNh9qJ4/vfkt/kDnty794hTLi5/quNun4vNLU1Os2Vx8tFPHRXhJQIiwSrEHQyIjJdhxXLKvuvlOonMhUev0duo4QEOT3EXEmqo69ymOe3q4JAIkIEQPSD0pJHYfqueVTnZcn2opjdMtsRFQddmeNMZknWItr6y0hB4uiQAJCNFD0pKs3DFnJNnpwReAPYfreXndzg43A40fnt2p4y0a3T40WQUgZlxybt9OHRfhJQEhekxqkpX7Fo8LdVzvO9LAi++UdmhYamFeGjMm5JOZYsWgQGaKlRkT8tsdxXQiTQeXW5bhiBWjB2Vy7ZTB5KTbMCgKOek2rp0yWEYxRchpZ1JHI5lJHdsyMuyUHarlz2u3h9qV87OTuOWq4dis4Rl1rQCZqQmYjOF/PxTPs3HjuW7QO+rXWXIHIXpcSqKFO+aOCk2mK69s5Pm122kM0zt9HWiQYa9CdJoEhIiIJJuZpXNG0j8ruAvd0eomnluznfowvZD7AlqPrzArRKyTgBARk5hgat7POgkIDnFctrqE6obwDGlsbPITUDs/B0OI3koCQkRUcO2mEZzTPxUIzmtYtqqEipqmbv9ZOlDfKE1NQnSUBISIOKvZyE0zhzFyYDoATref59aUcNDR/R2GflULW1+HEPFGAkJEBZPRwKIrh/KdoX0AcHtV/rx2BzsP1nb7z2p0+/HKMhxCnJEEhIgaRoPCNVOGcPGYXCD4bv+Vd3fy1a6qbv9Z9Y1e6Y8Q4gwkIERUMSgKsy4sYMaEfCA40W3Fh3vZ9PWRbt0XXdOhzumVWdZCnIYEhIg6iqIw5bz+XDtlMIbmLR3WfX6QNR+XoWnd94Ie0HTptBbiNCQgRNQaNyyb700fhtkU/Ge6ZbuD1/65q1t3p/P6VZkfIcQpSECIqDa8IJ2lc0ZiTwguw7HjQC1/XrsDZzcu5S3zI4RonwSEiHr52UnctWB0aOOh8spG/rSqBEc3zZWQpTiEaJ8EhIgJmSkJ3DV/FAU5wQXHap1e/rSqpN0tR7vCF9BokqYmIVqRgBAxw54Q3Of6vHOCcyW8fpWX3inl05KKbhnh5HRLU5MQJ5KAEDHFZDSw8PIhXDEuDwgOV13zcRmrNu8/6xd3XQ8uxRFjK+ALETYSECLmKIrCFePyuOGKczAZg+NgP99RyQv/2HHWy2j4VY2GJmn8sDMpAAAeDElEQVRqEgIkIEQMO3dIH+6cN4oUuwWAsqNO/u+tbzlc1XhW13V7A9IfIQQSECLG9c9K4odXjyY/O7hkeF2jj2dXl/DlWS7P4Wzyd3i/bCHilQSEiHnJiRbumDuSC4ZnAxBQdVZ+uPes+iV0oE7WaxK9nASEiAsmo4GrLx3M1ZMHYWxen+Oz7Q6WrS6h1unt0jU1HWqcEhKi95KAEHHlghE5fH/eSFKb+yUOVbn4/ZvFlHZx2XBN0yUkRK8V1oDYtGkTM2bMYNq0aSxbtqzN91evXs3cuXOZO3cuN9xwA6WlpeEsjugl8rOTuefaMRTmBXepc3tVXl63k3e2HOjSC72m6dRKSIheKGwBoaoqjz32GM8//zxFRUWsXbuWPXv2tDonLy+PV199lTVr1nD33Xfzi1/8IlzFEb2MPcHMzTOHc8W4PJoXhOWj4qMsW11CTRf2vFY1nep6D84mmScheo+wBURxcTEFBQXk5+djsViYPXs269evb3XOd77zHVJTg+/yzjvvPCoqKsJVHNELGQzB+RK3zh5Bss0MBJucfvfGt3y951inr6cDLk+AY/UevD4Z4STinylcF3Y4HOTm5oYe5+TkUFxcfMrzV65cyaWXXtqha2dk2M+6fNFM6te9JmTYGT64Dy8VbadkXzVev8rrG/awv8LJounDSEwwd+m69mRru8/Nyko+2yJHrXiuG8R//TorbAHR3m24oijtnAlbtmxh5cqV/PWvf+3QtWtqXGdVtmiWkWGX+oXJoivO4ZPsJN79/CCqprN1u4NdB2pZePkQBvdL7fT1ampcpCRaSEw4/meUlZVMVZWzO4sdNeK5btA76tdZYWtiys3NbdVk5HA4yM7ObnNeaWkpDz30EP/3f/9Henp6uIojBAZF4ZJz+/KDq0eTnW4DoN7l489rd7D2k7IubUTU0OSTWdciboUtIMaMGUNZWRnl5eX4fD6KioqYOnVqq3OOHDnCvffey3//938zaNCgcBVFiFb6Ztr54dVjuGh0sAlUBz7ZVsHvVn7LgYrOv4NsaPJT1+jt1u1QhYgGYWtiMplMPPzwwyxduhRVVbn22mspLCxk+fLlACxatIg//OEP1NXV8eijjwJgNBp58803w1UkIULMJgNzLhrI8IJ03ty4l7pGH9UNHpatLmHS6FymXZCP1Wzs8PU8PhWf301yqi2MpRaiZyl6jI3Zq6xpwhHH7YTSB9HzSvZX886Wg9ScMOM6yWaiT2oCAVUnPdnK+OHZFOalnfFaGRl2GhvcJNnMWDoRMNv2V7O5+ChVdW6y0mxccm5fRg/K7FJ9utvaT8v48KvDuDwB7AkmLju/P3MmDYx0sbqd9EG0FbY7CCFiwe5DdXzw1RESrCYylOB+EKqm0+gO0OhuxGY14ld13v28HKBDIeELaNQ4vVjNRpJsZsym07fkbttfzRsb94UeO2rdoceRDom1n5ax9uMyIDjIpLHJH3ocjyEhWpOlNkSv9kVpZejrBIuJrHRbaC0nCM7CrqxtwuX2s3WHo1PX9vpVqhs81Dq9+E/TAb65+GinjvekD7863KnjIr7IHYTo1U5eyM+gKBgMgAIKBgKqFtxpzuWjqTxAeaWT/OzO3ap7/Spev4rFZMBuM7fp26iqc7f7vKq6zs/47m6n2oDJdZYbM4nYIHcQoldLT7a2OWY0GDAbjWSlJZBit9Ayfccf0Pjj2yWs/HAvziZfp3+WL6BR6/RS0+DBe8JeE1lp7XdsZ6UldPpndLckW/uTCO2nOC7iiwSE6NXGD287NycxwYQ9wYSiKCTZzGSn2UiwHH/X/+WuKp7++zds/Pow/kDnF/A7OSguObdvu+ed6nhPuuz8/p06LuKL8ZFHHnkk0oXoDJfbj6sL795ihc1mwR3Ht+/RVr/MlATSk63UOb14fSoZKVamfiePYQPSQ8f6pCYwY8IAxg3L4kiVC5cngKrp7D3cwNe7q7DbzGSn21AUpVP1UzUdj08lLclKbkYidY1e3F6V7HQbMycOiHgHNcDQ/DRQ4Ei1i4CqYbeZmT5hQFx2UNvtVpri+LXFbm97t3wmMsw1ykTjMNDuFOv1UzWdz7c7eP9fh3B7A6Hj/bPszLhgABPO7dfl+pmNBhITTNis0dk12BuGgcZ7/TpLAiLKxPoL6JnES/3c3gAffnWYT7ZVoJ4wg3pYQTpTz+8f2iO7K4wGBZvVhNVsPOMQ2Z7UG15A471+nRWdb1WEiHI2q4mrLixg4sgc/vlFOd/sqQZg54Fadh6oZfiAdK4Y15/+WZ0PiuA8DD+Nbj8Gg4LVbMRmMXZq4p0Q3UECQoizkJGSwPVTC7l0bD/e21rOzoN1AJQerKX0YC0jCtK5/Dv9yetCUEBwNzu3N4DbG8BkULAlmEiwGDEaoufOQsQvCQghukHfTDs3zxxOTZOfNzfsZt+RBgB2HKhlx4FaCvNSmXJefwb1TT7lsvdnEtB0nE1+nE1+LCYDVosRq9mIyShhIcJDAkKIbnROXhpL54xk/9EG1v/rUCgodh+qZ/ehegbkJDH53H6MKEjHYOhaUEBwqKwvoOHEj0EBk9GAyWTAbDRgNhkkNES3kIAQIgwG9U0JBcXGr4+wqzzY9HTQ0chr/9xFZkoCF4/J5TtDs866b0HTjwdGi1BoNH+YTQomo6HLdy+id5KAECKMBvVNYVDfFI4cc7Hx68Ns21+DrkN1g4fVH5fx3tZyxg/P5sKROWSkdN/M6fZCQwGMRgWzyYjFJHca4swkIIToAf362Fl05VBqGjx8vK2Cf5VW4gtoeHwqm4uP8nHxUYYNSGfCyGyG5qWdVfPTqehAQNUJqAHczUtQGQzKCc1SSuiOQwiQgBCiR2WkJDD3ooFcOS6PL0or+bSkgrpGHzrHRz6lJVkYPzybcUOzSE3q/OzXztA0Ha+mtlobquVOI9Q8ZTRgMknTVG8kE+WiTLxMJDsVqV9rmqZTerCWLSUO9hyub/U9BTgnL5Vxw7IZUZAe8UlzfTKTaKhvatWvYTQaMMRJv4ZMlGtL7iCEiCCDQWHkwAxGDsygut7D5zsc/GtXFU2eADrHRz8lWIyMHpzJ+YV9KMhNjsiLsqbrbfo1WupgMigYDcHAMDZ/bTAo0lwV4yQghIgSmakJXHVhAdMuyGfHgVr+tbOK3Yfq0PXgntdflFbyRWklqXYLY4Zkcu6QTPr3sUd8ZJKm6fhCy4203hhJARSDglEJBkYwRJo/GwwYjUrc3IHEIwkIIaKMyWhgzOBMxgzOpN7l4+vdVXy1+xiVtcGNhepdPjYXH2Vz8VEykq2MGpTB6MEZ9M9KiroXWx3QNR0N/eTsCFEUQoFhUIJ3JAaDEgyXk+oTbBDXg9dtziSDEjyv5a7FaFAiHprxQgJCiCiWarcw5bz+XDq2H0erm/h69zG+3VdNvSu4LHWN08tHxUf5qPgoKXYLIwrSGVGQzuB+KTHTvKPrLaOrTr0ta2cZlODugIqioDQHSEuQKC0hdEKonGnUmK7r6HqwmU3XdVpumIzN1wnHqLNoIAEhRAxQFIV+fez062Nn5oUDKHc0Ury3mm37q3E2BfefaHD5+Gy7g8+2O7CYDJyTl8qw/DSG5qeFfTRUtNGaX8yD9zAdoxoM1Ne5g+1iBJ+q0xIKp3/uyU1pwQA6HkQtdzZnCpNg+Bz/eYoCCsdDrqdJQAgRYwyKQkFuMgW5ycy+qIBDlY2U7K+hZH8NNc17bPsCGtvLatleVgtATrqNwrw0zslLZWBusqwM2w5dD6531aXncuamtBYtYdKSE7p+/A7ldD9daf5PS1AoJxxruVs6OXs0PdhHpOu6jGISorcxKAoDcpIZkJPMzIkDqKxzU9q8QGB5ZWPonaij1o2j1s3mb49iNCjk5yQxuG8Kg/ulkJ+dHPEhtL3J8TDp/PNoDpNTfLfbSUAIEScURSEnPZGc9ESmnNefJk+APYfr2FVex+7yepzNW6Gqmk7ZUSdlR51s+PIwJqNCXlYSBbnJDMwNhk207monepb8KxAiTiUmmDh3SB/OHdIHXddx1LrZe7iePYfqKatwhmZPB1SdsgonZRVONjY/NyvNxoCcJAZkJ5GXnUR2emLkKiIiJuYCwmoxkmAxhvqRUJTjbXO07shRlJb+pmAnz/FjoW80P6ftOc1H0DneYdQyxK7l69AtHzoKJzX+Ka0+ha55vE0w2Pml6XrwdrP5Z1jNwYXU2nNy3ULXbxmpccI57T9WmtsrW9f3RK3bQ5s/nzC08HTdZCcOPQxeSw8db/kiMcGE22oKnaiHzj25rid8fdJB5aRvHP9/3H7pdF0/4XcV/GEnljVUTr2llsfPPfH80PNjkKIo5GYkkpuRyMVj+qJqOkeOudh3pJ6yo04OOJx4fMcbz6vq3FTVufnXzioguF/2gL7JZKfZ6N/cWZ6VZsMYp6N3RFDMLbUBxP10eKlfdNNDAaOjac1B35zwGZlJVFU5UTUdVdUIaMe/F800XcdR08SBCicHHY0crHRS0+A97XNMRoXs9ET6ZiSSm5lITkYiOek2kmzmmJyHEO/LwIwZltPp58TcHYQQkaaE7loVTp5qkGQz47aZWx3zBzQ8vgAen4oapWFhUBT6Ztrpm2nnwlHBY41uP4erGimvbORwlYsj1a7QkFoINk0dOebiyLHWL6qJCSay021kp9nICn0kkJpkjbqJfOL0JCCECDOzyYDZZCE5EQKq1ubuouWjq0MswyXJZmbYgHSGDUgHID09kbLyWo4cc3H4mIuKmiaOVjdR62x9p9HkCYQ6wU9kMir0SbWRmZJAZqqVzJQEMlISyEixkmK3SnNVFJKAEKIHBVdCBdqZh6BpOl6/GlwQzx99dxuKopCaZCU1ycqIgRmh4x5fAEeNm4qaJhy1TThqgv0XjW5/q+cHVJ2KmiYqapraXNugKKQlW0hLspKebA19Tk0KHku1W2JmZng8kYAQIkoYDAo2qwlb86TngKoFA8OvoWoamh4cPx9dsQEJFlNo4t6JmjyBUGf3sXo3VXUejtV7qGnwtAk/TdepafCett/DnmAi1W4hxW4hOTH4OSXRTHKihaSWzzYTRoMESXeRgBAiSrXsu2A/aSfSgKrhb77L8AW0qLvTaJGY0H5waJpOQ5OPY/Ueap1eahqCoVHX6KPW6W1z59HC5Qng8gQ4Ut32DuRENquJJJuZJJsZu82EPSH4dWKCCXuCicQEM4lWE4kJwQ+z7NV9ShIQQsSYluBomcymaTr+5tAIqMEPVY2+O40WBoNCWlKwGak9voBKXaOP+kYv9Y0+6hq91Lt8NLh81Lt81Df6Wu2AdzK3N4DbG7x76QiTMXjnlpRowWIyYLOYsFmNJFhMJFiDw+oTLKbmz8EPq9mE1WzAajFiMRllsT4hRHQyGBSsBiPWE/o1dF1Hbe78VrXgfBtV1YMBEuVDby0mI9lpwVFQp+Lzqzib/NS7fDS6fTib/Dib/DS6W3+43P4z3mEFVD30/K4ymwzBOUzmls/B+UwnfzabDFhMxuaBC60/TM17g4f2CG/+2mSM3L4ZYQ2ITZs28fjjj6NpGgsXLuT73/9+q+/rus7jjz/Oxo0bSUhI4Fe/+hWjRo0KZ5GEiFvb9lezufgoVXVustJsXHJuX0YPyuTPRdvZuqMSv6phNhq4YEQ2E0Zk81HzuX1SEpg4KpfhA9LZXlbDlpIKqhs8ZKQkMG5YFocqXXy2vYImn0qixciEkTlcfn7eKcux+1AdX5RWUuv0kp5sZfzwbIA2xwrz0jr8/MK8ND746hCfb3fQ5A2QaDWdthwt16hp8AQ71QvSyUqzsedwPaUHanE2+bGYDWQkJ2A2GXD7AvgDOs4mH25vAK9P7dQdmD8QvIOjYzctXWI0tOwTroRCw2Q0BHfza3ncvAlTy/4aJuPxnf66Mg8ibBPlVFVlxowZ/OUvfyEnJ4frrruOp59+mnPOOSd0zsaNG3nllVd47rnn+Oabb3j88cdZsWLFGa8d6xOtTiceJpKdjtQvPLbtr+aNjfvaHLcnmNjRvKJrC51g/0DWSe/Qxw3LCs2cblHX6MXl9odm5be8XMy8sIAZFwxARw+tDKDrOjsO1lL06YHQD9IJjnJSAKul9fvRGRPy24TE7kN1vPt5eZt65GbY+HZvdZvjl32nf5uQONU1Rg1Kp2R/bZvjLeU4caLcrvJa1n1W3mr/B13TGT04g7QkKx6fitff/NH8tT+ghR63bM3qaz4eDdb8z/xOPydsdxDFxcUUFBSQn58PwOzZs1m/fn2rgFi/fj0LFixAURTOO+88GhoaqKysJDs7O1zFEiIubS4+2u7x0gNtXxAB3J5Am2MffnWY5ERLq2ONTX50Xcd0Ukfux98e5erJg9tc4+vdx9oMR61rDPYXJCVaWq1E+u2+ai4cmQscX+q6eG91u/Mh2gsHgK07KpkzaSBwfBmUr3cfa7Psdcu5Sc2TGE98V/zVripGD8ps7lswhq5hNCoYT1pcpsHlY97Fg9oty6lout5qUIH/hODwqxp+f/CzL9Dytdrcn6SHBiT4m/uVAif2NbXMp1GPz6tp6X9StebjZ9mUGLaAcDgc5Obmhh7n5ORQXFx82nNyc3NxOBxnDIiurGseS6R+sS0S9att9LW7ZLem03bNreYX45PPd3kCZKS0HjKlNb+gH1/TK/i5yRNot57tlUNVdVDa/ryGJj/9+qa2OWa1tJ0jEtB0LKa2xz0+lcEFma2ONXoCJLSzGm1VvYc+7fRrNPlUhg7uA0BmavD7Lq/a7oq2TT6VYUOy2hyPV2ELiPZartruL3vmc4QQZ/b0j6dEugjA2ZejO+oRLdeIB2GbUZKbm0tFRUXocXt3BiefU1FRIc1LQggRJcIWEGPGjKGsrIzy8nJ8Ph9FRUVMnTq11TlTp07l7bffRtd1vv76a5KTkyUghBAiSoSticlkMvHwww+zdOlSVFXl2muvpbCwkOXLlwOwaNEipkyZwsaNG5k2bRo2m40nnngiXMURQgjRSTG5H4QQQojwk1WthBBCtEsCQgghRLuiei0mr9fLkiVL8Pl8oZnZP/rRj6irq+MnP/kJhw8fpn///vzv//4vqampZ75gFGrpn8nJyeHZZ5+Nq7pNnToVu92OwWDAaDTy5ptvxlX9GhoaeOihh9i1axeKovDEE08waNCguKjfvn37+MlPfhJ6XF5ezo9+9CMWLFgQF/V78cUXWbFiBYqiMHToUJ588kncbndc1A3gpZdeYsWKFei6zsKFC7nlllu69LcX1XcQFouFl156idWrV/P222/z0Ucf8fXXX7Ns2TImTZrEe++9x6RJk1i2bFmki9plL7/8MkOGDAk9jqe6QfAf6qpVq3jzzTeB+Krf448/zuTJk1m3bh2rVq1iyJAhcVO/wYMHs2rVqtDvzmazMW3atLion8Ph4OWXX+aNN95g7dq1qKpKUVFRXNQNYNeuXaxYsYIVK1awatUqPvzwQ8rKyrpUv6gOCEVRsNvtAAQCAQKBAIqihJboAFiwYAHvv/9+JIvZZRUVFXz44Ydcd911oWPxUrdTiZf6NTY2snXr1tDvzmKxkJKSEjf1O9Gnn35Kfn4+/fv3j5v6qaqKx+MhEAjg8XjIzs6Om7rt3buXsWPHYrPZMJlMXHDBBfzzn//sUv2iOiAg+IucP38+F110ERdddBFjx46luro6NF8iOzubmpqaCJeya5544gn+7d/+DcMJO2DFS91a3H777VxzzTX8/e9/B+KnfuXl5WRkZPCzn/2MBQsW8POf/5ympqa4qd+JioqKmDNnDhAfv7+cnBxuu+02Lr/8ci655BKSkpK45JJL4qJuAEOHDuWLL76gtrYWt9vNpk2bqKio6FL9oj4gjEYjq1atYuPGjRQXF7Nr165IF6lbfPDBB2RkZDB69OhIFyVsli9fzltvvcVzzz3Ha6+9xtatWyNdpG4TCATYvn07ixYt4u2338Zms8Vsk8Tp+Hw+NmzYwMyZMyNdlG5TX1/P+vXrWb9+PR999BFut5tVq1ZFuljdZsiQISxdupTbbruNpUuXMmzYMIzGtutYdUTUB0SLlJQUJk6cyEcffURmZiaVlZUAVFZWkpGRcYZnR58vv/ySDRs2MHXqVO677z62bNnC/fffHxd1a5GTE1x/PjMzk2nTplFcXBw39cvNzSU3N5exY8cCMHPmTLZv3x439WuxadMmRo0aRZ8+zYvZxUH9PvnkE/Ly8sjIyMBsNjN9+nS++uqruKhbi4ULF/LWW2/x2muvkZaWRkFBQZfqF9UBUVNTQ0NDAwAej4dPPvmEwYMHh5boAHj77be54oorIlnMLvnpT3/Kpk2b2LBhA08//TQXXnghv/71r+OibgBNTU00NjaGvv74448pLCyMm/plZWWRm5vLvn3BPRg+/fRThgwZEjf1a1FUVMTs2bNDj+Ohfv369eObb77B7Xaj63pc/u6qq4PLox85coT33nuPOXPmdKl+UT2TurS0lAcffBBVVdF1nZkzZ3LPPfdQW1vLj3/8Y44ePUrfvn357W9/S1pa+7tTxYLPPvuMF154gWeffTZu6lZeXs4Pf/hDINiPNGfOHO6+++64qR/Ajh07+PnPf47f7yc/P58nn3wSTdPipn5ut5vLLruM999/n+Tk4NLe8fL7e+aZZ/jHP/6ByWRixIgRPP7447hcrrioG8DixYupq6vDZDLxs5/9jEmTJnXpdxfVASGEECJyorqJSQghRORIQAghhGiXBIQQQoh2SUAIIYRolwSEEEKIdkX1aq5CnM7ChQvx+Xz4/X7KysooLCwEYOTIkTz55JMRLl3HlJSUUF5eHlczlUX8kGGuIuYdOnSIa6+9ls8++yzSRWkjEAhgMp36fdiKFSv45JNP+M1vftPt1xbibMm/LhGXVq5cyd/+9jdUVSUlJYVHH32UgQMHsmLFCtatW4fdbmfXrl307duX//iP/+Cpp56ivLycsWPH8tRTT6EoCvfffz82m42DBw9SUVHBxIkT+cUvfoHZbMbpdPLEE0+we/duvF4vF110EQ888AAGg4FFixYxYcIEvvrqKxITE3nmmWdCkwS9Xi9jx47l0UcfpaGhgT/84Q+4XC7mz5/PxIkTWbJkCYsXL+bjjz8G4MCBA6HHBw4cYNGiRVx//fVs2bKFa665hvnz5/P000/zxRdf4PP5GDFiBI888gg2my3CvwERF3QhYlx5ebk+YcKE0OMtW7bod955p+71enVd1/X169frS5Ys0XVd119//XV9woQJekVFha7run7bbbfpCxYs0J1Op+7z+fRZs2bpW7Zs0XVd13/605/q8+fP110ul+7z+fSbbrpJ/+tf/6rruq4/8MAD+po1a3Rd13VVVfUf/ehH+sqVK3Vd1/UbbrhB/8EPfqAHAoHQ9+vq6kJf33ffffrrr78eKs+Pf/zjUNnLysr0iy66qN3HZWVl+tChQ/V169aFvv/MM8/ozz77bOjxk08+qf/2t789u/+hQjSTOwgRdzZs2MD27dtZuHAhALqu43K5Qt8fN25caCHBkSNH4vF4SEpKAmDYsGEcPHiQiRMnAjBr1iwSExOB4Br6H374IYsWLeKDDz6gpKSE5557DgiuFTZgwIDQz5g7d25oBU1N01i2bBmbN29G0zTq6uq6vFNZYmIiM2bMaFVXt9tNUVEREFx9ddSoUV26thAnk4AQcUfXdb773e9yzz33tPt9q9Ua+tpgMLR5HAgETnldRVGA4Iv+s88+S79+/do9tyVUAFatWkVxcTF//etfsdvt/P73v+fo0aPtPs9oNKJpWuix1+s95XVbyvTLX/6SCy64oN3rCXE2ZJiriDstq1Y6HA4guFjgtm3bunStd955B7fbjd/vZ82aNaE7i6lTp7Js2TJUVQWCKw+Xl5e3ew2n00l6ejp2u536+vrQu30Au92O0+kMPc7Ozsbj8YSutXbt2jPW9YUXXggFSWNjI3v37u1SXYU4mQSEiDsXXngh99xzD3feeSfz5s1j7ty5fPjhh1261rhx47j77ruZM2cO+fn5oS1Gf/GLX6BpGvPnz2fu3LnccccdVFVVtXuNq6++mrq6OubMmcN9993X6t3+xRdfjNPpZN68eTzxxBNYLBYefPBBbr75Zm688UbMZvNpy3fXXXcxZMgQrrvuOubOncuSJUvYv39/l+oqxMlkmKsQp3D//fczbtw4Fi1aFOmiCBERcgchhBCiXXIHIYQQol1yByGEEKJdEhBCCCHaJQEhhBCiXRIQQggh2iUBIYQQol3/Pxb9I/eeQAYgAAAAAElFTkSuQmCC\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/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/5c9dbef11b4d7638b7ddf2ea71026e7bf00fcfb0/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] } ], "metadata": { "celltoolbar": "Hide code", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }