{ "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, and numpy library." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.6.5rc1 (default, Mar 14 2018, 06:54:23) \n", "[GCC 7.3.0]\n", "uname_result(system='Linux', node='icarus', release='4.15.0-2-amd64', version='#1 SMP Debian 4.15.11-1 (2018-03-20)', machine='x86_64', processor='')\n", "\n", "IPython 5.5.0\n", "IPython.core.release 5.5.0\n", "PIL 4.3.0\n", "PIL.version 4.3.0\n", "_csv 1.0\n", "_ctypes 1.1.0\n", "_curses b'2.2'\n", "decimal 1.70\n", "argparse 1.1\n", "csv 1.0\n", "ctypes 1.1.0\n", "cvxopt 1.1.9\n", "cycler 0.10.0\n", "dateutil 2.7.3\n", "decimal 1.70\n", "decorator 4.3.0\n", "distutils 3.6.5rc1\n", "ipaddress 1.0\n", "ipykernel 4.8.2\n", "ipykernel._version 4.8.2\n", "ipython_genutils 0.2.0\n", "ipython_genutils._version 0.2.0\n", "ipywidgets 6.0.0\n", "ipywidgets._version 6.0.0\n", "joblib 0.11\n", "json 2.0.9\n", "jupyter_client 5.2.3\n", "jupyter_client._version 5.2.3\n", "jupyter_core 4.4.0\n", "jupyter_core.version 4.4.0\n", "logging 0.5.1.2\n", "matplotlib 2.1.1\n", "matplotlib.backends.backend_agg 2.1.1\n", "matplotlib.pylab 1.14.5\n", "numexpr 2.6.5\n", "numpy 1.14.5\n", "numpy.core 1.14.5\n", "numpy.core.multiarray 3.1\n", "numpy.core.umath b'0.4.0'\n", "numpy.lib 1.14.5\n", "numpy.linalg._umath_linalg b'0.1.5'\n", "numpy.matlib 1.14.5\n", "optparse 1.5.3\n", "pandas 0.22.0\n", "_libjson 1.33\n", "patsy 0.5.0\n", "patsy.version 0.5.0\n", "pexpect 4.2.1\n", "pickleshare 0.7.4\n", "pkg_resources._vendor.packaging.__about__ 16.8\n", "pkg_resources._vendor.six 1.10.0\n", "pkg_resources._vendor.appdirs 1.4.0\n", "pkg_resources._vendor.packaging 16.8\n", "pkg_resources._vendor.pyparsing 2.1.10\n", "pkg_resources._vendor.six 1.10.0\n", "platform 1.0.8\n", "prompt_toolkit 1.0.15\n", "ptyprocess 0.5.2\n", "py 1.5.3\n", "py._vendored_packages.apipkg 1.4\n", "pytest 3.3.2\n", "pygments 2.2.0\n", "pyparsing 2.2.0\n", "pytz 2018.5\n", "re 2.2.1\n", "scipy 0.19.1\n", "scipy._lib.decorator 4.3.0\n", "scipy._lib.six 1.2.0\n", "scipy.fftpack 0.4.3\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", "simplejson 3.15.0\n", "six 1.11.0\n", "statsmodels 0.9.0\n", "statsmodels.__init__ 0.9.0\n", "traitlets 4.3.2\n", "traitlets._version 4.3.2\n", "urllib.request 3.6\n", "zlib 1.0\n", "zmq 17.0.0\n", "zmq.sugar 17.0.0\n", "zmq.sugar.version 17.0.0\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", "\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(\"https://app-learninglab.inria.fr/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/master/data/shuttle.csv\")\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAGBNJREFUeJzt3XuQnXWd5/H3t5MACYmAxMkwBAQGhpUCJkDLReaSCFqBKpN1AYUtwXEGM1uSskRHYWZdhmGdqpVRmXFlFGRxhC2NXEbIzmaWixAdprgFiOEmTA8gdEDAGCANIemkv/vHefrxpOnLOZ1++vQ5vF9VqZznOb9++vvtp09/+rn070RmIkkSQFerC5AkTR2GgiSpZChIkkqGgiSpZChIkkqGgiSpVFkoRMTVEfFSRDwywvMREV+PiJ6IWBcRR1dViySpMVUeKfwDsHiU508BDin+LQO+WWEtkqQGVBYKmfkT4FejDFkKXJM19wB7RsQ+VdUjSRrb9BZ+7n2B5+qWe4t1LwwdGBHLqB1NMHPmzGP222+/SSmwUQMDA3R1dd7lmU7tCzq3N/tqP5PV25NPPvnLzHzXWONaGQoxzLph59zIzCuBKwG6u7tzzZo1VdbVtNWrV7Nw4cJWlzHhOrUv6Nze7Kv9TFZvEfHzRsa1Mnp7gfpf+ecDz7eoFkkSrQ2FlcA5xV1IxwOvZuZbTh1JkiZPZaePIuL7wEJgbkT0An8JzADIzG8Bq4BTgR7gDeATVdUiSWpMZaGQmWeN8XwC51X1+SVJzevMy/mSpHExFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklQyFCRJJUNBklSqNBQiYnFEPBERPRFx4TDP7x8Rd0bEQxGxLiJOrbIeSdLoKguFiJgGXA6cAhwGnBURhw0Z9kXgusw8CjgT+Puq6pEkja3KI4VjgZ7MfCoztwIrgKVDxiTwjuLxHsDzFdYjSRpDZGY1G444HVicmecWy2cDx2Xm8rox+wC3AnsBuwMnZ+YDw2xrGbAMYN68ecesWLGikprHq6+vj9mzZ7e6jAnXqX1B5/ZmX+1nsnpbtGjRA5nZPda46RXWEMOsG5pAZwH/kJlfjYgTgGsj4vDMHNjhgzKvBK4E6O7uzoULF1ZR77itXr2aqVbTROjUvqBze7Ov9jPVeqvy9FEvsF/d8nzeenroT4DrADLzbmA3YG6FNUmSRlFlKNwPHBIRB0bELtQuJK8cMuZZ4CSAiHgPtVB4ucKaJEmjqCwUMnMbsBy4BXic2l1Gj0bEJRGxpBj2OeCTEfFT4PvAH2VVFzkkSWOq8poCmbkKWDVk3UV1jx8DTqyyBklS4/yLZklSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUMBUlSyVCQJJUqDYWIWBwRT0RET0RcOMKYj0TEYxHxaER8r8p6JEmjm97IoIg4PDMfaWbDETENuBz4ANAL3B8RKzPzsboxhwB/DpyYmRsj4jea+RySpInV6JHCtyLivoj4VETs2eDHHAv0ZOZTmbkVWAEsHTLmk8DlmbkRIDNfanDbkqQKRGY2NrD2W/0fA2cA9wHfyczbRhl/OrA4M88tls8GjsvM5XVjbgKeBE4EpgEXZ+b/G2Zby4BlAPPmzTtmxYoVjXU3Sfr6+pg9e3ary5hwndoXdG5v9tV+Jqu3RYsWPZCZ3WMOzMyG/1H7wX0asB54HPgZ8J9GGHsGcFXd8tnA/xwy5p+AHwIzgAOpnWbac7QajjnmmJxq7rzzzlaXUIlO7Suzc3uzr/YzWb0Ba7KBn/MNnT6KiCMj4rIiCN4PfCgz31M8vmyED+sF9qtbng88P8yYmzOzPzOfBp4ADmmkJknSxGv0msI3gAeB383M8zLzQYDMfB744ggfcz9wSEQcGBG7AGcCK4eMuQlYBBARc4HfAZ5qrgVJ0kRp6O4j4FRgc2ZuB4iILmC3zHwjM68d7gMyc1tELAduoXba6erMfDQiLqF2GLOyeO6DEfEYsB34fGZu2MmeJEnj1Ggo3A6cDPQVy7OAW4H3jfZBmbkKWDVk3UV1jxP4bPFPktRijZ4+2i0zBwOB4vGsakqSJLVKo6HwekQcPbgQEccAm6spSZLUKo2ePvoMcH1EDN49tA/w0WpKkiS1SkOhkJn3R8R/AA4FAvhZZvZXWpkkadI1eqQA8F7ggOJjjooIMvOaSqqSJLVEoxPiXQv8NrCW2q2jAAkYCpLUQRo9UugGDituIZUkdahG7z56BPjNKguRJLVeo0cKc4HHIuI+YMvgysxcUklVkqSWaDQULq6yCEnS1NDoLak/joh3A4dk5u0RMYvafEaSpA7S6NTZnwRuAK4oVu1LbYZTSVIHafRC83nU3h3tNYDM/DfA91OWpA7TaChsydr7LAMQEdOp/Z2CJKmDNBoKP46IvwBmRsQHgOuB/1NdWZKkVmg0FC4EXgYeBv6U2nskjPSOa5KkNtXo3UcDwLeLf5KkDtXo3EdPM8w1hMw8aMIrkiS1TDNzHw3aDTgDeOfElyNJaqWGrilk5oa6f+sz82+B91dcmyRpkjV6+ujousUuakcOcyqpSJLUMo2ePvpq3eNtwDPARya8GklSSzV699GiqguRJLVeo6ePPjva85n5tYkpR5LUSs3cffReYGWx/CHgJ8BzVRQlSWqNZt5k5+jM3AQQERcD12fmuVUVJkmafI1Oc7E/sLVueStwwIRXI0lqqUaPFK4F7ouIH1L7y+YPA9dUVpUkqSUavfvoryPin4HfL1Z9IjMfqq4sSVIrNHr6CGAW8Fpm/h3QGxEHVlSTJKlFGn07zr8ELgD+vFg1A/jfVRUlSWqNRo8UPgwsAV4HyMzncZoLSeo4jYbC1sxMiumzI2L36kqSJLVKo6FwXURcAewZEZ8Ebsc33JGkjtPo3UdfKd6b+TXgUOCizLyt0sokSZNuzCOFiJgWEbdn5m2Z+fnM/LNGAyEiFkfEExHRExEXjjLu9IjIiOgeaYwkqXpjhkJmbgfeiIg9mtlwREwDLgdOAQ4DzoqIw4YZNwf4NHBvM9uXJE28Rv+i+U3g4Yi4jeIOJIDM/PQoH3Ms0JOZTwFExApgKfDYkHH/HbgU+LNGi5YkVaPRUPi/xb9m7MuOs6j2AsfVD4iIo4D9MvOfImLEUIiIZcAygHnz5rF69eomS6lWX1/flKtpInRqX9C5vdlX+5lqvY0aChGxf2Y+m5nfHce2Y5h1WbftLuAy4I/G2lBmXglcCdDd3Z0LFy4cRznVWb16NVOtponQqX1B5/ZmX+1nqvU21jWFmwYfRMSNTW67F9ivbnk+8Hzd8hzgcGB1RDwDHA+s9GKzJLXOWKFQ/9v+QU1u+37gkIg4MCJ2Ac7k12/SQ2a+mplzM/OAzDwAuAdYkplrmvw8kqQJMlYo5AiPx5SZ24DlwC3A48B1mfloRFwSEUuaK1OSNBnGutD8uxHxGrUjhpnFY4rlzMx3jPbBmbkKWDVk3UUjjF3YUMWSpMqMGgqZOW2yCpEktV4z76cgSepwhoIkqWQoSJJKhoIkqfS2CYUNfVv46XOvsKFvS6tLkdSEDX1b2Ny/3dfuJHlbhMLNa9dz4pfv4GNX3cuJX76DlWvXt7okSQ0YfO0+/fLrvnYnSceHwoa+LVxw4zre7B9g05ZtvNk/wBduXOdvHdIUV//a3Z7pa3eSdHwo9G7czIyuHduc0dVF78bNLapIUiN87bZGx4fC/L1m0j8wsMO6/oEB5u81s0UVSWqEr93W6PhQ2Hv2rlx62pHsNqOLObtOZ7cZXVx62pHsPXvXVpcmaRT1r91pEb52J0mjb7LT1pYs2JcTD55L78bNzN9rpt9UUpsYfO3ed/dd/OuS3/O1OwneFqEAtd86/IaS2s/es3dl5oxpvn4nScefPpIkNc5QkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUslQkCSVDAVJUqnSUIiIxRHxRET0RMSFwzz/2Yh4LCLWRcSPIuLdVdYjSRpdZaEQEdOAy4FTgMOAsyLisCHDHgK6M/NI4Abg0qrqkSSNrcojhWOBnsx8KjO3AiuApfUDMvPOzHyjWLwHmF9hPZKkMURmVrPhiNOBxZl5brF8NnBcZi4fYfw3gF9k5peGeW4ZsAxg3rx5x6xYsaKSmserr6+P2bNnt7qMCdepfUHn9mZf7Weyelu0aNEDmdk91rjpFdYQw6wbNoEi4mNAN/CHwz2fmVcCVwJ0d3fnwoULJ6jEibF69WqmWk0ToVP7gs7tzb7az1TrrcpQ6AX2q1ueDzw/dFBEnAz8V+APM3NLhfVIksZQ5TWF+4FDIuLAiNgFOBNYWT8gIo4CrgCWZOZLFdYiSWpAZaGQmduA5cAtwOPAdZn5aERcEhFLimF/A8wGro+ItRGxcoTNSZImQZWnj8jMVcCqIesuqnt8cpWfv51t6NtC78bNzN9rJnvP3nXCxraTTu2rKj0vbmLjG/30vLiJg+fNaXU5alOVhoLG5+a167ngxnXM6Oqif2CAS087kiUL9t3pse2kU/uqykU3Pcw19zzL547YxvmX/YRzTtifS5Ye0eqy1Iac5mKK2dC3hQtuXMeb/QNs2rKNN/sH+MKN69jQ99Zr8M2MbSed2ldVel7cxDX3PLvDumvufpaeFze1qCK1M0NhiunduJkZXTvulhldXfRu3LxTY9tJp/ZVlbXPvdLUemk0hsIUM3+vmfQPDOywrn9ggPl7zdypse2kU/uqyoL99mxqvTQaQ2GK2Xv2rlx62pHsNqOLObtOZ7cZXVx62pHDXmhtZmw76dS+qnLwvDmcc8L+O6w754T9vdiscfFC8xS0ZMG+nHjw3IbuvGlmbDvp1L6qcsnSIzjn+AN4+IF7uP384w0EjZuhMEXtPXvXhn8QNjO2nXRqX1U5eN4cemfNMBC0Uzx9JEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpJKhIEkqGQqSpFKloRARiyPiiYjoiYgLh3l+14j4QfH8vRFxQJX1SM3a0LeFnz73Chv6tow6bs3TG/jarU+w5ukNE7bNZsf2vLiJjW/00/PipjHHNqOqepv5/Jv7tzf8NbhhzXMd9zWocrtDTa9qwxExDbgc+ADQC9wfESsz87G6YX8CbMzMgyPiTODLwEerqklqxs1r13PBjeuY0dVF/8AAl552JEsW7PuWcR+76h7u6qmFwdfv6OH3D96ba889fqe22ezYi256mGvueZbPHbGN8y/7CeecsD+XLD1inJ1XX2+zn//T7+nn/C/f0dDXYFCnfA2q3O5wqjxSOBboycynMnMrsAJYOmTMUuC7xeMbgJMiIiqsSWrIhr4tXHDjOt7sH2DTlm282T/AF25c95bf0tY8vaEMhEH/0rNh2COGRrfZ7NieFzft8MMQ4Jq7n93p35arqnc8n3975tvya1DldkcSmVnNhiNOBxZn5rnF8tnAcZm5vG7MI8WY3mL534sxvxyyrWXAsmLxUOCJSooev7nAL8cc1X46tS8Yo7eYMXPW9L32+Z3o6po2uC4HBrZv2/jCk9m/+Y3BddPmzP2tabvvuc/Qj9/++isvbN/0y+fHs81mx3bN2mPv6e941wEA2994lWmz9gBg22svPzPwxqtjn8/aya9Bs2PH8/kH+2rka1CvTb4GE/K92IB3Z+a7xhpU2ekjYLjf+IcmUCNjyMwrgSsnoqgqRMSazOxudR0TrVP7gs7tLSLWbHv1JftqI1Pte7HK00e9wH51y/OB50caExHTgT2AX1VYkyRpFFWGwv3AIRFxYETsApwJrBwyZiXw8eLx6cAdWdX5LEnSmCo7fZSZ2yJiOXALMA24OjMfjYhLgDWZuRL4X8C1EdFD7QjhzKrqqdiUPbW1kzq1L+jc3uyr/Uyp3iq70CxJaj/+RbMkqWQoSJJKhsI4RMQzEfFwRKyNiDXFuosjYn2xbm1EnNrqOpsVEXtGxA0R8bOIeDwiToiId0bEbRHxb8X/e7W6zmaN0Fcn7K9D6+pfGxGvRcRn2n2fjdJXJ+yz8yPi0Yh4JCK+HxG7FTfj3Fvsrx8UN+a0rkavKTQvIp4Buuv/yC4iLgb6MvMrraprZ0XEd4F/ycyrim/MWcBfAL/KzP9RzF+1V2Ze0NJCmzRCX5+hzfdXvWJamfXAccB5tPk+GzSkr0/QxvssIvYF7gIOy8zNEXEdsAo4FfjHzFwREd8CfpqZ32xVnR4pCICIeAfwB9TuCCMzt2bmK+w4Fcl3gf/YmgrHZ5S+Os1JwL9n5s9p8302RH1fnWA6MLP4u6xZwAvA+6lN8wNTYH8ZCuOTwK0R8UAxBceg5RGxLiKubrdDduAg4GXgOxHxUERcFRG7A/My8wWA4v/faGWR4zBSX9De+2uoM4HvF4/bfZ/Vq+8L2nifZeZ64CvAs9TC4FXgAeCVzNxWDOsFqpnprkGGwvicmJlHA6cA50XEHwDfBH4bWEBth3+1hfWNx3TgaOCbmXkU8DrwlunO29BIfbX7/ioVp8SWANe3upaJNExfbb3PihBbChwI/BawO7WfIUO19Jy+oTAOmfl88f9LwA+BYzPzxczcnpkDwLepzRLbTnqB3sy8t1i+gdoP0xcjYh+A4v+XWlTfeA3bVwfsr3qnAA9m5ovFcrvvs0E79NUB++xk4OnMfDkz+4F/BN4H7FmcToLhpwOaVIZCkyJi94iYM/gY+CDwyOCLsPBh4JFW1DdemfkL4LmIOLRYdRLwGDtORfJx4OYWlDduI/XV7vtriLPY8RRLW++zOjv01QH77Fng+IiYFRHBr19jd1Kb5gemwP7y7qMmRcRB1I4OoHZq4nuZ+dcRcS21w9oEngH+dPC8bruIiAXAVcAuwFPU7vboAq4D9qf2TX1GZrbVpIUj9PV12nx/AUTELOA54KDMfLVYtzftv8+G66sTXmN/Re2NxLYBDwHnUruGsAJ4Z7HuY5lZ7durjVajoSBJGuTpI0lSyVCQJJUMBUlSyVCQJJUMBUlSqbJ3XpMmW3Er5o+Kxd8EtlOb4gJqf2C4tSWFjSIi/hhYVfw9hdRy3pKqjjSVZq2NiGmZuX2E5+4Clmfm2ia2N71urhxpQnn6SG8LEfHxiLivmIf/7yOiKyKmR8QrEfE3EfFgRNwSEcdFxI8j4qnB+foj4tyI+GHx/BMR8cUGt/uliLgPODYi/ioi7i/m0f9W1HyU2h9j/aD4+F0iojci9iy2fXxE3F48/lJEXBERt1Gb3G96RHyt+NzrIuLcyf+qqhMZCup4EXE4tWkR3peZC6idNj2zeHoP4NZigsOtwMXUph84A7ikbjPHFh9zNPCfI2JBA9t9MDOPzcy7gb/LzPcCRxTPLc7MHwBrgY9m5oIGTm8dBXwoM88GlgEvZeaxwHupTcy4/3i+PlI9ryno7eBkaj8419SmnGEmtSkUADZn5m3F44eBVzNzW0Q8DBxQt41bMnMjQETcBPwetdfPSNvdyq+nQwE4KSI+D+wGzKU2ZfI/N9nHzZn5ZvH4g8B7IqI+hA6hNq2FNG6Ggt4OArg6M//bDitrM1PW/3Y+AGype1z/+hh68S3H2O7mLC7YFfP4fIPa7KzrI+JL1MJhONv49RH80DGvD+npU5n5I6QJ5OkjvR3cDnwkIuZC7S6lcZxq+WDU3ut5FrU58f+1ie3OpBYyvyxm2D2t7rlNwJy65WeAY4rH9eOGugX41OCUy1F7X+OZTfYkvYVHCup4mflwMTvl7RHRBfQD/4Xm5q2/C/getTd5uXbwbqFGtpuZG6L2PtGPAD8H7q17+jvAVRGxmdp1i4uBb0fEL4D7RqnnCmqzoK4tTl29RC2spJ3iLanSGIo7ew7PzM+0uhapap4+kiSVPFKQJJU8UpAklQwFSVLJUJAklQwFSVLJUJAklf4/GKF1l7kqzyEAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n", "import matplotlib.pyplot as plt\n", "\n", "data[\"Frequency\"]=data.Malfunction/data.Count\n", "data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression\n", "\n", "Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature." ] }, { "cell_type": "code", "execution_count": 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: Sun, 23 Sep 2018 Deviance: 3.0144
Time: 22:50:48 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: Sun, 23 Sep 2018 Deviance: 3.0144\n", "Time: 22:50:48 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']], 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.*\n", "**I have therefore managed to partially 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": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VOXZ//HPNUs2sgABwhI2NYDIngUQa8EqoFXcUEDEpSD2qUutlVb6WLVWuzz0+blXoYBrFakVROsjCIoLIgQEWWVHSNiXhITsyfX7YwYMMZAhmWSWXO/XK6/MOXOfc647J/nOyZkz9xFVxRhjTHhxBLoAY4wx/mfhbowxYcjC3RhjwpCFuzHGhCELd2OMCUMW7sYYE4ZqDHcRmSkiB0Rk3WmeFxF5RkS2isgaEenn/zKNMcacDV+O3F8Ghp/h+cuBFO/XROCFupdljDGmLmoMd1X9DDhyhiZXA6+qx1dAUxFp468CjTHGnD2XH9bRDthdaTrLO29v1YYiMhHP0T3R0dGp7du3r9UGKyoqcDjC4+0C60vwCZd+gPUlWNWlL5s3bz6kqi1rauePcJdq5lU7poGqTgOmAaSlpemKFStqtcHFixczePDgWi0bbKwvwSdc+gHWl2BVl76IyHe+tPPHy2AWUPkQPBnY44f1GmOMqSV/hPs84BbvVTMDgFxV/cEpGWOMMQ2nxtMyIvImMBhoISJZwCOAG0BVXwQ+AK4AtgIFwO31Vawxxhjf1BjuqjqmhucVuMtvFRljQkJpaSlZWVkUFRU1yPYSEhLYuHFjg2yrvvnSl6ioKJKTk3G73bXahj/eUDXGNEJZWVnExcXRqVMnRKq7rsK/8vLyiIuLq/ftNISa+qKqHD58mKysLDp37lyrbYTHdUXGmAZXVFREYmJigwR7YyMiJCYm1um/Igt3Y0ytWbDXn7r+bC3cjTEmDNk5d2NMyHI6nfTs2fPk9Ny5c+nUqVPgCgoiFu7GmJAVHR3N6tWrT/t8WVkZLlfjjDk7LWOMCSsvv/wyN9xwA1dddRVDhw4FYMqUKaSnp9OrVy8eeeSRk22feOIJunbtyqWXXsqYMWP429/+BsDgwYM5MTzKoUOHTv43UF5ezqRJk06ua+rUqcD3wwmMHDmSbt26MXbsWDxXiUNmZiYXXnghvXv3JiMjg7y8PIYNG3bKi9KgQYNYs2aNX38OjfMlzRjjV394bz0b9hzz6zq7t43nkasuOGObwsJC+vTpA0Dnzp2ZM2cOAEuXLmXNmjU0b96cBQsWsGXLFpYvX46qMmLECD777DOaNGnCrFmzWLVqFWVlZfTr14/U1NQzbm/GjBkkJCSQmZlJcXExgwYNOvkCsmrVKtavX0/btm0ZNGgQS5YsISMjg1GjRvHWW2+Rnp7OsWPHiI6O5pZbbuHll1/mqaeeYvPmzRQXF9OrVy8//NS+Z+FujAlZpzstc9lll9G8eXMAFixYwIIFC+jbty8A+fn5bNmyhby8PK699lpiYmIAGDFiRI3bW7BgAWvWrOHtt98GIDc3ly1bthAREUFGRgbJyckA9OnTh507d5KQkECbNm1IT08HID4+HoBrr72WQYMGMWXKFGbOnMltt91Wtx9ENSzcjTF1VtMRdkNr0qTJyceqyuTJk7nzzjtPafPUU0+d9nJDl8tFRUUFwCnXmqsqzz77LMOGDTul/eLFi4mMjDw57XQ6KSsrQ1Wr3UZMTAyXXXYZ7777LrNnz6a2I+SeiZ1zN8aEtWHDhjFz5kzy8/MByM7O5sCBA1x88cXMmTOHwsJC8vLyeO+9904u06lTJ1auXAlw8ij9xLpeeOEFSktLAdi8eTPHjx8/7ba7devGnj17yMzMBDyfTC0rKwNgwoQJ3HvvvaSnp5/8L8Of7MjdGBPWhg4dysaNGxk4cCAAsbGxvP766/Tr149Ro0bRp08fOnbsyI9+9KOTyzzwwAPceOONvPbaa1xyySUn50+YMIGdO3fSr18/VJWWLVsyd+7c0247IiKCt956i3vuuYfCwkKio6NZuHAhAKmpqcTHx3P77fU01qKqBuQrNTVVa+uTTz6p9bLBxvoSfMKlH6r125cNGzbU27qrc+zYsXpd/yOPPKJTpkyp122ccOzYMc3OztaUlBQtLy8/bbvqfsbACvUhY+20jDHGNLA33niD/v3788QTT9TbrQPttIwxxgCPPvpog23rpptu+sEbvP5mR+7GmFpTrfZ2ycYP6vqztXA3xtRKVFQUhw8ftoCvB+odzz0qKqrW67DTMsaYWklOTiYrK4uDBw82yPaKiorqFHbBxJe+nLgTU21ZuBtjasXtdtf6LkG1sXjx4pOfMg11DdEXOy1jjDFhyMLdGGPCkIW7McaEIQt3Y4wJQxbuxhgThizcjTEmDFm4G2NMGLJwN8aYMGThbowxYcjC3RhjwlDIhfuBY0V8mlVqgxUZY8wZhFy4v75sFy+tK2HCKys4mFcc6HKMMSYohVy43/eTFG7qFsHnWw8x/KnPWLhhf6BLMsaYoBNy4e5wCEM7ufnPPReRFB/FhFdX8Mi76ygqLQ90acYYEzRCLtxPSEmKY85dFzL+os68svQ7rnl+CVsP5Ae6LGOMCQohG+4AkS4nv7+yOy/dns6BvGJGPPcF767ODnRZxhgTcD6Fu4gMF5FNIrJVRB6s5vkOIvKJiKwSkTUicoX/Sz29IV1b8Z97L+KCtvH8ctZqfjdnLcVldprGGNN41RjuIuIEngcuB7oDY0Ske5VmDwGzVbUvMBr4u78LrUmbhGjevGMAP//xubyxbBc3vriU7JzChi7DGGOCgi9H7hnAVlXdrqolwCzg6iptFIj3Pk4A9vivRN+5nA4evLwbU8elsv3gca585nO+3HooEKUYY0xASU0fBhKRkcBwVZ3gnR4H9FfVuyu1aQMsAJoBTYBLVXVlNeuaCEwESEpKSp01a1atis7Pzyc2NvaMbfYdr+CZVUXsO66M6hrB0I4uRKRW26tPvvQlVIRLX8KlH2B9CVZ16cuQIUNWqmpajQ1V9YxfwA3A9ErT44Bnq7S5H/i19/FAYAPgONN6U1NTtbY++eQTn9rlFZXqHa9kasffvq/3v7Vai0rLar3N+uJrX0JBuPQlXPqhan0JVnXpC7BCa8htVfXptEwW0L7SdDI/PO0yHpjtfbFYCkQBLXxYd72KjXTx4s2p3HdpCv/+Ooux/1jGoXz7VKsxJvz5Eu6ZQIqIdBaRCDxvmM6r0mYX8BMAETkfT7gf9GehteVwCPdd2oXnb+rHuj25XP3cEjbtywt0WcYYU69qDHdVLQPuBuYDG/FcFbNeRB4TkRHeZr8G7hCRb4A3gdu8/z4EjZ/2asPsOwdSWl7ByBe+5PMtQfHaY4wx9cKn69xV9QNV7aKq56rqE955D6vqPO/jDao6SFV7q2ofVV1Qn0XXVq/kpsy9axDtmkVz+0uZvJW5K9AlGWNMvQjpT6jWRtum0fzr5wO58LwW/Pbfa3lq4WYbPtgYE3YaXbgDxEW5mXFrGiNTk3lq4RYmv7OWsvKKQJdljDF+4wp0AYHidjqYMrIXbROieObjrRzKL+a5m/oR5XYGujRjjKmzRnnkfoKIcP/Qrvzx6gtY9O0BbpmxnNzC0kCXZYwxddaow/2EcQM78czovqzafZRRU5faHZ6MMSHPwt3rqt5tmXFrOjsPH2fUVBt0zBgT2izcK7m4S0teG9+fg3nF3PjiUnYeOh7okowxplYs3KtI79ScNycOoLC0nBunLrW7OxljQpKFezV6tEtg1sQBVCiMnraUb/cdC3RJxhhzVizcT6NLUhyz7xyAy+Fg9LSvWL8nN9AlGWOMzyzcz+CclrG8decAYtxOxk5fxrpsC3hjTGiwcK9Bx8QmzJo40ALeGBNSLNx90CExhlkTB9IkwsnNM5axca+dgzfGBDcLdx91SIzhzYkDiHI5uXn6MrbstzHhjTHBy8L9LHRMbMIbd/TH6RDG/GMZ2w/aZZLGmOBk4X6WzmkZyxt39EdVGTt9GbuPFAS6JGOM+QEL91o4r1Ucr43vT0FJOWOnL2NfblGgSzLGmFNYuNdS97bxvPKzDI4cL+HmGcs4crwk0CUZY8xJFu510Kd9U6bfmsbuIwXcOnM5eUU2XLAxJjhYuNfRgHMSeeHmfmzce4wJr6ygqLQ80CUZY4yFuz9c0i2J/72xN8t3HuHuN1bZLfuMMQFn4e4nV/dpx6NXXcDCjfuZ/M5au+m2MSagGu09VOvDrRd24vDxEp5ZtIXE2EgevLxboEsyxjRSFu5+9qtLUzhyvJgXP91Gq7hIfnZR50CXZIxphCzc/UxE+MOIHhzKK+GP/9lAy7hIrurdNtBlGWMaGTvnXg+cDuGp0X1I79ic+2ev5suthwJdkjGmkbFwrydRbif/uCWNzi2acOdrK+1uTsaYBmXhXo8SYty8dHsGMZFObn8pk725hYEuyRjTSFi417N2TaN56bYM8orKuP2lTPsUqzGmQVi4N4DubeN54eZ+bD2Qzy/++TWl9iEnY0w9s3BvID9Kacmfru3J51sO8dCcdfYhJ2NMvbJLIRvQjent2X20gGc/3kqHxBjuGnJeoEsyxoQpC/cGdv9lXdh1pIAp8zfRMTGG2EAXZIwJS3ZapoGJCH+9vhdpHZtx/+xv2HrURpE0xvifT+EuIsNFZJOIbBWRB0/T5kYR2SAi60XkDf+WGV6i3E6m3ZJGm4Qonl5VZLfqM8b4XY3hLiJO4HngcqA7MEZEuldpkwJMBgap6gXAffVQa1hp3iSCmbelU14B41/J5JhdImmM8SNfjtwzgK2qul1VS4BZwNVV2twBPK+qRwFU9YB/ywxP57aM5e6+UWw/eNzGgTfG+JXUdEmeiIwEhqvqBO/0OKC/qt5dqc1cYDMwCHACj6rqh9WsayIwESApKSl11qxZtSo6Pz+f2NjweCsyPz+flUcjeWl9CZd2cHFz98hAl1Rr4bJfwqUfYH0JVnXpy5AhQ1aqalpN7Xy5WkaqmVf1FcEFpACDgWTgcxHpoao5pyykOg2YBpCWlqaDBw/2YfM/tHjxYmq7bLBZvHgxj1w5GOf7G5j+xQ4G9+vGzQM6BrqsWgmX/RIu/QDrS7BqiL74clomC2hfaToZ2FNNm3dVtVRVdwCb8IS98dHkK85nSNeWPDJvvY0iaYypM1/CPRNIEZHOIhIBjAbmVWkzFxgCICItgC7Adn8WGu6cDuGZMX05t2UTfv76SnYcOh7okowxIazGcFfVMuBuYD6wEZitqutF5DERGeFtNh84LCIbgE+ASap6uL6KDldxUW6m35KO0yGMfyWT3EK7gsYYUzs+Xeeuqh+oahdVPVdVn/DOe1hV53kfq6rer6rdVbWnqtbunVJDh8QYXrg5lV2HC7jnTbuCxhhTO/YJ1SA04JxE/nhNDz7bfJA//9+3gS7HGBOCbGyZIDUmowOb9uUx44sddGsdxw1p7WteyBhjvOzIPYg99NPzGXReIv89Zx0rvzsa6HKMMSHEwj2IuZwOnr+pH22aRnHnayvtNn3GGJ9ZuAe5pjERTL8ljaLScia+upKiUhtF0hhTMwv3EJCSFMdTo/qwbk8uv/33GruLkzGmRhbuIeLS7kk8MLQr767ew9TP7PNhxpgzs3APIb8YfC4/7dWGv374LYs32cCbxpjTs3APISLClJG96NY6nnveXMX2g/mBLskYE6Qs3ENMTISLaeNScTsd3PHqCrvJhzGmWhbuIah98xiev6kfOw8XcP9bq6mosDdYjTGnsnAPUQPPTeThK7uzcOMBnly4OdDlGGOCjA0/EMJuGdiR9XtyefbjrZzfJp4rerYJdEnGmCBhR+4hTET44zU96NuhKQ/86xu+3Xcs0CUZY4KEhXuIi3Q5efHmVGIjXdzx6gpyCkoCXZIxJghYuIeBpPgoXhyXyv7cYu5+w8aAN8ZYuIeNfh2a8fg1Pfhi6yH+YmPAG9Po2RuqYeTG9Pas25PL9C92cEG7eK7tmxzokowxAWJH7mHm91d2p3/n5jz477WsycoJdDnGmACxcA8zbqeDv4/tR4vYSO58bSUH84oDXZIxJgAs3MNQYmwkU8elcrSghF/8cyUlZfYGqzGNjYV7mOrRLoG/Xt+LzJ1HefS99YEuxxjTwOwN1TB2dZ92bNh7jKmfbueCtvGM7d8x0CUZYxqIHbmHud8M68aPu7TkkXfXs3zHkUCXY4xpIBbuYc7pEJ4Z05f2zWP4r9dXkp1jN9k2pjGwcG8EEqLd/OOWNErKKpj46goKS+wm28aEOwv3RuK8VrE8PaYPG/YeY9Lb39hNto0Jcxbujcgl3ZKYNKwr76/Zy98Xbwt0OcaYemTh3sj814/P5arebfnbgk0s2rg/0OUYY+qJhXsjIyL8z/W9uKBtPL+ctZot+/MCXZIxph5YuDdC0RFOpo1LI8rtZIKNAW9MWLJwb6TaNo1m6rhU9uYU8Yt/fk2pjQFvTFixcG/EUjs240/X9eTLbYf54/sbAl2OMcaPbPiBRm5kajKb9+cx7bPtpCTFMW6ADVFgTDiwI3fDb4d3Y0jXljw6bz1fbj0U6HKMMX7gU7iLyHAR2SQiW0XkwTO0GykiKiJp/ivR1LcTQxSc06IJ//XPr9lx6HigSzLG1FGN4S4iTuB54HKgOzBGRLpX0y4OuBdY5u8iTf2Li3Iz49Z0nA5h/MuZ5BaUBrokY0wd+HLkngFsVdXtqloCzAKurqbdH4H/AYr8WJ9pQB0SY3jx5lR2Hy3grjfsChpjQpnUNMaIiIwEhqvqBO/0OKC/qt5dqU1f4CFVvV5EFgMPqOqKatY1EZgIkJSUlDpr1qxaFZ2fn09sbGytlg02wdiXz7NKmbGuhMHtXdzaPQIR8Wm5YOxLbYRLP8D6Eqzq0pchQ4asVNUaT337crVMdX/ZJ18RRMQBPAncVtOKVHUaMA0gLS1NBw8e7MPmf2jx4sXUdtlgE4x9GQxEfPgtLyzexo96d2H8RZ19Wi4Y+1Ib4dIPsL4Eq4boiy+nZbKA9pWmk4E9labjgB7AYhHZCQwA5tmbqqFt0tCuDL+gNY//ZwMLN9gYNMaEGl/CPRNIEZHOIhIBjAbmnXhSVXNVtYWqdlLVTsBXwIjqTsuY0OFwCE+O6kOPtgncO2sV67JzA12SMeYs1BjuqloG3A3MBzYCs1V1vYg8JiIj6rtAEzjREU5m3JpG02g341/JZG+u3cXJmFDh03XuqvqBqnZR1XNV9QnvvIdVdV41bQfbUXv4aBUfxczb0zleXM7tL2WSV2SXSBoTCuwTqqZG3VrH8/zYfmw5kM9db6yySySNCQEW7sYnP+7Skieu6cFnmw/y8Lvr7DZ9xgQ5GzjM+Gx0Rgd2Hy3g+U+2kdwshruGnBfokowxp2Hhbs7Kry/rStbRQqbM30SbhCiu65cc6JKMMdWwcDdnxeEQpozszcG8Yn7z9hpaxUVxUUqLQJdljKnCzrmbsxbhcvDiuFTOaxXLz19fyfo9tb8Gfu6qbAb95WM6P/gfBv3lY+auyvZjpaa+2f4LXhbuplbio9y8fHsG8VEubnspk91HCs56HXNXZTP5nbVk5xSiQHZOIZPfWWsBESJs/wU3C3dTa60Tonh1fAYlZRXcMnM5x0rO7gqaKfM3UVhafsq8wtJypszf5M8yTT2x/RfcLNxNnZzXKo6Zt6WxN7eQJ1cUkV9c5vOye3Kq/8Tr6eab4GL7L7hZuJs6S+3YnL+P7cd3eRVMfHUFxWXlNS8EtG0afVbzTXCx/RfcLNyNX1zSLYnxPSL4ctth7pu1mvKKmk/RTBrWlWi385R50W4nk4Z1ra8yjR/Z/gtuFu7Gbwa1c/P7K7vzf+v2MfmdNTV+ivWavu3483U9adc0GgHaNY3mz9f15Jq+7RqmYFMntv+Cm13nbvxq/EWdyS0s5ZlFW4iPcvPfPz3/jHdyuqZvOwuDEGb7L3hZuBu/+9WlKRwrLGX6FzuIi3Lzy0tTAl2SMY2OhbvxOxHh4Su7k19cxpMLNxMT4eSOi88JdFnGNCoW7qZeOBzCX6/vRWFpOU98sJEot4NxAzsFuixjGg0Ld1NvnA7hyRv7UFRSzu/fXU+Ey8Go9A6BLsuYRsGuljH1KsLl4Pmx/bi4S0sefGct/16ZFeiSjGkULNxNvYtyO5k2LpULz01k0tvf8O5qG3vEmPpm4W4aRJTbyfRb0sno3JxfvbXaBpcypp5ZuJsGEx3hZOZt6fTvnMj9s1czZ5WdojGmvli4mwYVE+Fi5m3pDDgnkftnf8PsFbsDXZIxYcnC3TS46AgnM25N56LzWvCbt9fw+lffBbokY8KOhbsJiOgIJ/+4JY2fdGvFQ3PXMeOLHYEuyZiwYuFuAibK7eSFm1O5vEdr/vj+Bp5euKXGwcaMMb6xcDcBFeFy8OyYvlzfL5knF27mif9stIA3xg/sE6om4FxOB1NG9iIuysX0L3aQW1jKn6/rictpxx7G1JaFuwkKDofwyFXdSYh28/SiLRwtKOW5m/oSVeVmEMYY39ihkQkaIsKvLuvCH0ZcwKJv9zNuxjJyCkoCXZYxIcnC3QSdWy/sxDOj+/LN7lxGvriUrKMFgS7JmJBj4W6C0lW92/LKzzLYf6yI6/7+JeuycwNdkjEhxcLdBK2B5yby9s8vxOkQbpy6lI+/3R/okowJGRbuJqh1bR3H3LsG0blFEya8soJXvtwZ6JKMCQkW7iboJcVHMfvOgVzSrRWPzFvPw++uo6y8ItBlGRPUfAp3ERkuIptEZKuIPFjN8/eLyAYRWSMii0Sko/9LNY1Zk0gXU8elMfHic3h16Xfc9lImuQWlgS7LmKBVY7iLiBN4Hrgc6A6MEZHuVZqtAtJUtRfwNvA//i7UGKdD+N0V5zNlZC+W7TjMiOe/YPP+vECXZUxQ8uXIPQPYqqrbVbUEmAVcXbmBqn6iqieuV/sKSPZvmcZ874a09syaOICCknKufX4JH67bG+iSjAk6UtM4HiIyEhiuqhO80+OA/qp692naPwfsU9XHq3luIjARICkpKXXWrFm1Kjo/P5/Y2NhaLRtsrC+1d7SogmdXFbM9t4IrOru5PsWN0yF1Xq/tk+BkffEYMmTISlVNq7Ghqp7xC7gBmF5pehzw7Gna3oznyD2ypvWmpqZqbX3yySe1XjbYWF/qpqi0TCe/s0Y7/vZ9HT11qR7MK6rzOm2fBCfriwewQmvIV1X16bRMFtC+0nQysKdqIxG5FPhvYISqFvuwXmPqLNLl5E/X9uRvN/Tm611HueLpz1m67XCgyzIm4HwJ90wgRUQ6i0gEMBqYV7mBiPQFpuIJ9gP+L9OYMxuZmszcuwYRG+li7PSveGbRFsorbOhg03jVGO6qWgbcDcwHNgKzVXW9iDwmIiO8zaYAscC/RGS1iMw7zeqMqTfnt4ln3j0XMaJ3W/7fR5sZO/0r9uYWBrosYwLCpyF/VfUD4IMq8x6u9PhSP9dlTK0s3LCf5TuOALBs+xF+8r+fMjq9PfPX72dPTiFtm0YzaVhXrunbzu/bnrsqmynzN9X7dnzx0Ny1vLlsN/f1KGX85A8Y0789j1/TMyC1mMCw8dxN2Ji7KpvJ76ylsLQcAAUKS8qZuWTnyTbZOYVMfmctgF+Dt+q262s7vnho7lpe/2rXyely1ZPTFvCNhw0/YMLGlPmbTobrCdWddS8sLWfK/E31vu362I4v3ly2+6zmm/Bk4W7Cxp4c38+vZ59F27ps+2xq8pfy03x25XTzTXiycDdho23TaJ/bOh3CF1sO1fu2z6Ymf3FK9R/kOt18E54s3E3YmDSsK9FV7rnqdghu56mhFuF00LxJBDfPWMYD//qGI8frfiu/6rYd7XYyaVjXOq/7bI3p3/6s5pvwZG+omrBx4o3LqlesVDdveI/WPLNoC9M+286ijfv53RXnMzI1Ganl0e3pth2Iq2VOvGl64hy7U8SulmmELNxNWLmmb7tqA7W6eb8Z3o0Rfdry0Jx1THp7DbNX7OYPI3r4fduB8Pg1PXn8mp4sXryYbWMHB7ocEwB2WsY0at1axzP7zoH89fqebDt4nCuf/ZzXNhSTU1D3UzXGBJKFu2n0HA5hVHoHPvn1YG4e0JGPd5Xx4ymLeWnJDkrtjk8mRFm4G+OVEOPmsat78NigaHq0i+cP721g2JOf8eG6fSdGPTUmZFi4G1NF+zgHr4/vz/Rb0nA4hJ+/vpKRLy5l2XYbbdKEDgt3Y6ohIlzaPYkPf/kj/nxdT3YfKWDUtK+4ZeZy1mTlBLo8Y2pk4W7MGbicDsZkdODTSUP43RXdWJOVw4jnljD+5UwLeRPULNyN8UF0hJOJF5/L578ZwgNDu7Diu6OMeG4Jt720nMydRwJdnjE/YOFuzFmIi3Jz9yUpfPHbIUwa1pW1Wbnc8OJSbnxxKQs37KfCbhBigoSFuzG1EBfl5q4h5/HFby/h4Su7k51TyIRXV3DZk5/yxrJdFJaU17wSY+qRhbsxdRAd4eRnF3Vm8aTBPD26D1FuJ7+bs5aBf1nEX/7vW3YfKQh0iaaRsuEHjPEDt9PB1X3aMaJ3WzJ3HuWlJTuY9tk2pn62jUu6tmLsgA78uEsrnA4bmdE0DAt3Y/xIRMjo3JyMzs3Zk1PIm8t38eby3Sx6eQVtE6K4Ia09I1OTad88JtClmjBn4W5MPWnbNJpfD+3KPZeksGjjft5YvotnPt7C04u2MPCcRK5PTWZ4j9bERtqfofE/+60ypp5FuBxc3rMNl/dsQ9bRAt75Opu3V2bxwL++4aG5a7mse2tG9G7LxV1aEOly1rxCY3xg4W5MA0puFsO9P0nhnkvO4+tdR5mzKpv31+zlvW/2EBflYmj31lzeozUXpbQgym1Bb2rPwt2YABARUjs2J7Vjcx656gKWbD3Ee9/s5aMN+/j311nERrr4cdeWDO2exOCurUiIdge6ZBNiLNyNCTC308Hgrq0Y3LUVJWU9+XLbIT5ct4+FGw/wnzV7cTqE9E7NuKSbp01Kq9ha3zHKNB4W7sYEkQjX90FfUaGs2p3Dx9/uZ9HGA/zpg2/50wff0iYhih+ltOCilJZceG4iLWIjA122CUIW7sbIwl0aAAANK0lEQVQEKYdDSO3YjNSOzZg0rBt7cgr5bPNBPt18kA/X7WP2iiwAurWOY8A5iQw4J5H0Ts1ItLA3WLgbEzLaNo1mdEYHRmd0oLxCWZudy5Kth1i67TCzMnfx8pc7ATivVSxp3heF8uMVqKqdxmmELNyNCUFOh9CnfVP6tG/KXUPOo7isnLVZuSzfeYTMHUf4YO1eZmXuBuAvKz+id3JTerdvSu/kBHomJ9AqLirAPTD1zcLdmDAQ6XKS1qk5aZ2aw2CoqFC2HcznjQVfURiTxOrdOTz38RZODFrZKi6SHu0S6N4mnu5t4zm/TTwdmsfY8AhhxMLdmDDkcAgpSXH8uL2bwYN7AVBQUsb6PcdYk5XL+j25rMvO5dPNByn3Jn6U20FKqzi6JMXRJSmWlKRYzmsZR7tm0Rb6IcjC3ZhGIibCRXqn5qR3an5yXlFpOVsP5LNh7zE27ctj0748Pt9ykH9/nXWyTaTLQecWTejcogmdWjShc2ITOibG0KlFE1rGRuKw4A9KFu7GNGJRbic92iXQo13CKfNzCkrYeiCfbQfz2Xognx2HjrNpfx4fbdhPWaUbkkS6HLRvHkNys2jaN4uhXbNo2jWNPvm9RWykHfUHiIW7MeYHmsZEfH8Ov5Ky8gr25BSx4/Bxdh0pYPeRAr47fJyso4Ws2pVDbmHpKe1dDiEpPorWCVG0jo8iKT6KpPhIWsVH0iouilZxkbSIjaRpjNuu6PEzC3djjM9cTgcdEmPokFj9kMV5RaVk5xSSfbSQPblF7M0pZF9uEfuOFbFx7zE+2XSAgmruUuV2ColNIkmMjSAxNpLEJhE0r/TVLCaC746U02ZfHs1i3MRHu23snRpYuBtj/CYuyk231m66tY4/bZv84jL2HyviwLFiDuQVcSi/hEP5xRzMK+bIcc/j7QfzOXK85AcvBH9e/tnJx1FuBwnRbuKj3J7v0W7io1zERbmJ836PjXIRF+miSaSLWO9Xk0gnsZEuYiJdxLidYfuegU/hLiLDgacBJzBdVf9S5flI4FUgFTgMjFLVnf4t1ZjwNXdVNlPmb2JPTiFtm0YzaVhX/rViF0u2HTnZZtC5zbkhrcMP2gE/mLfiuyO8uWw39/UoZfzkDxjTvz2PX9PTp+1Wt75r+rbzue4T2y5XxSnyg23HRrqIbRnL2qzcGrf96FUp/KhLC44cL+HTpSvokHI+RwtK+WrbYT7dfJD9x4q9p4JiKCwtZ+uBMo4VlZJXVHbyKqCaRLudxEQ4iY448d1FtNtBTISLaLeTSLeDaLeTKLeTKLeDKNf3jyNdnucjXd7HLgeRbgcRTicRLsf3X87vv7udgmr930i9xnAXESfwPHAZkAVkisg8Vd1Qqdl44Kiqnicio4G/AqPqo2Bjws3cVdlMfmcthaWeo9TsnELue2v1D9ot2XbklLDPzilk0tvfgEKpN8iycwq5/63VVFRarlyV17/aBXBKyFa33Un/+gYESsu/X9/kd9YC/CDgq1ve39t+ZN56/nxdT67p246DiU4G92rL3FXZfPztgZPLFpVWkHW08GQ7AFWlqLSCvOJS8ovKyC/2fhWVUVBSzvGSMo4Xl3G8uJzjxWUUlJZTWFJOQUkZhaUVFJaUcTCvmELv/KLScgpLPd99fM04o3HdIxhS99WckS9H7hnAVlXdDiAis4CrgcrhfjXwqPfx28BzIiLaEC9PxoS4KfM3nQyqs3UiCCurqKYdwJvLdp8SsNVtt7Sa5CosLWfK/E0/CPfqlm+IbVe3bNV2IkK092i8VdxpiqoFVaW0XCkqK6e4tIKi0nJKyisoLq2guKyckrIKissqKCmr8MwvK6e0TCku98wrLa+gtKyCuPxd/ivqNKSm/BWRkcBwVZ3gnR4H9FfVuyu1Wedtk+Wd3uZtc6jKuiYCE72TXYFNtay7BXCoxlahwfoSfBq0HxGtz0utr3WXF+TijPn+MseSfVtX1na7lZet6/K1XLYFcOhMy1atMYjV5Xeso6q2rKmRL0fu1b3bUPUVwZc2qOo0YJoP2zxzQSIrVDWtrusJBtaX4BMu/QBPX8pyD4RNX8Jpv9R3Xxw+tMkC2leaTgb2nK6NiLiABOAIxhhjAsKXcM8EUkSks4hEAKOBeVXazANu9T4eCXxs59uNMSZwajwto6plInI3MB/PpZAzVXW9iDwGrFDVecAM4DUR2YrniH10fRaNH07tBBHrS/AJl36A9SVY1XtfanxD1RhjTOjx5bSMMcaYEGPhbowxYSjow11EokRkuYh8IyLrReQP3vmdRWSZiGwRkbe8b/YGPRFxisgqEXnfOx2q/dgpImtFZLWIrPDOay4iH3n78pGINAt0nb4QkaYi8raIfCsiG0VkYCj2RUS6evfHia9jInJfiPblV96/93Ui8qY3B0L1b+WX3n6sF5H7vPPqfZ8EfbgDxcAlqtob6AMMF5EBeIY4eFJVU4CjeIZACAW/BDZWmg7VfgAMUdU+la7XfRBY5O3LIu90KHga+FBVuwG98eyfkOuLqm7y7o8+eMZ5KgDmEGJ9EZF2wL1Amqr2wHMhx4lhTULqb0VEegB34Pmkf2/gShFJoSH2iaqGzBcQA3wN9Mfz6S6Xd/5AYH6g6/Oh/mTvjrwEeB/Ph79Crh/eWncCLarM2wS08T5uA2wKdJ0+9CMe2IH34oJQ7kuV+ocCS0KxL0A7YDfQHM8Vfe8Dw0LxbwW4Ac9giyemfw/8piH2SSgcuZ84lbEaOAB8BGwDclS1zNskC88vRLB7Cs+OPTEERyKh2Q/wfAJ5gYis9A4rAZCkqnsBvN9bBaw6350DHARe8p4umy4iTQjNvlQ2GnjT+zik+qKq2cDfgF3AXiAXWElo/q2sAy4WkUQRiQGuwPOBz3rfJyER7qparp5/NZPx/HtzfnXNGraqsyMiVwIHVLXy2Bc+DdsQpAapaj/gcuAuEbk40AXVkgvoB7ygqn2B4wT5aYuaeM9FjwD+FehaasN7/vlqoDPQFmiC5/esqqD/W1HVjXhOJ30EfAh8A5SdcSE/CYlwP0FVc4DFwACgqXeoA6h+SIRgMwgYISI7gVl4Ts08Rej1AwBV3eP9fgDPed0MYL+ItAHwfj8QuAp9lgVkqeoy7/TbeMI+FPtywuXA16q63zsdan25FNihqgdVtRR4B7iQ0P1bmaGq/VT1Yjwf8txCA+yToA93EWkpIk29j6Px7PiNwCd4hjoAz9AH7wamQt+o6mRVTVbVTnj+Zf5YVccSYv0AEJEmIhJ34jGe87vrOHUYipDoi6ruA3aLSFfvrJ/gGc465PpSyRi+PyUDodeXXcAAEYkREeH7fRJyfysAItLK+70DcB2efVPv+yToP6EqIr2AV/C8Y+4AZqvqYyJyDp4j4ObAKuBmVS0OXKW+E5HBwAOqemUo9sNb8xzvpAt4Q1WfEJFEYDbQAc8f6A2qGvQDyIlIH2A6EAFsB27H+7tG6PUlBs+bkeeoaq53XsjtF+8lz6PwnMJYBUzAc449pP5WAETkczzvr5UC96vqoobYJ0Ef7sYYY85e0J+WMcYYc/Ys3I0xJgxZuBtjTBiycDfGmDBk4W6MMWHIlxtkG9OgvJeJLfJOtgbK8QwRAJChqiUBKewMRORnwAfe6+aNCTi7FNIENRF5FMhX1b8FQS1OVS0/zXNfAHer6uqzWJ+r0lgpxviVnZYxIUVEbhXP+P6rReTvIuIQEZeI5IjIFBH5WkTmi0h/EflURLaLyBXeZSeIyBzv85tE5CEf1/u4iCwHMkTkDyKS6R2f+0XxGIVnOOq3vMtHiEhWpU9WDxCRhd7Hj4vIVBH5CM9gZS4R+X/eba8RkQkN/1M14cjC3YQM79jY1wIXegeSc/H9zdgTgAXewcxKgEfxfGz9BuCxSqvJ8C7TD7hJRPr4sN6vVTVDVZcCT6tqOtDT+9xwVX0LWA2MUs946jWdNuoLXKWq44CJeAaUywDS8QzC1qE2Px9jKrNz7iaUXIonAFd4hhwhGs9H7QEKVfUj7+O1QK6qlonIWqBTpXXMV9WjACIyF7gIz9/B6dZbwvdDLQD8REQmAVFACzxD0f7fWfbjXVUt8j4eCpwvIpVfTFLwfCTdmFqzcDehRICZqvr7U2Z6RgqsfLRcgecOXiceV/49r/omk9aw3kL1vjHlHbflOaCfqmaLyON4Qr46ZXz/n3HVNser9OkXqroIY/zITsuYULIQuFFEWoDnqppanMIYKp57psbgGTN8yVmsNxrPi8Uh76iY11d6Lg+IqzS9E8+t7qjSrqr5wC9ODGUrnvugRp9ln4z5ATtyNyFDVdd6RwtcKCIOPKPs/ZyzG9f7C+AN4FzgtRNXt/iyXlU9LCKv4Bne+DtgWaWnXwKmi0ghnvP6jwL/EJF9wPIz1DMVz8iAq72nhA7gedExpk7sUkjTaHivROmhqvcFuhZj6pudljHGmDBkR+7GGBOG7MjdGGPCkIW7McaEIQt3Y4wJQxbuxhgThizcjTEmDP1/pGenSMj5bcYAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n", "data_pred['Frequency'] = logmodel.predict(data_pred)\n", "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "No confidence region was given in the original article. I have tried to compute and draw the confidence region in python but I haven't found how to do so. **I have failed so far to obtain the confidence region**. Here are my attempts" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " low up OR\n", "Intercept -9.569730 19.739685 5.084977\n", "Temperature -0.341358 0.110156 -0.115601\n" ] } ], "source": [ "# Inspiring from http://blog.yhat.com/posts/logistic-regression-and-python.html\n", "# odds ratios and 95% CI\n", "params = logmodel.params\n", "conf = logmodel.conf_int()\n", "conf['OR'] = params\n", "conf.columns = ['low', 'up', 'OR']\n", "\n", "#conf.low.Temperature = conf.OR.Temperature-2*0.047 ## I know my previous estimates of error are wrong. What if I fixed them ?\n", "#conf.up.Temperature = conf.OR.Temperature+2*0.047\n", "#conf.low.Intercept = conf.OR.Intercept-2*3.052\n", "#conf.up.Intercept = conf.OR.Intercept+2*3.052\n", "\n", "print(conf)\n", "def logit_inv(x):\n", " return(np.exp(x)/(np.exp(x)+1))\n", "\n", "data_pred['Prob']=logit_inv(data_pred['Temperature'] * conf.OR.Temperature + conf.OR.Intercept)\n", "\n", "# mean_temp = np.mean(data.Temperature)\n", "# mean_prob_logit = mean_temp * conf.OR.Temperature + conf.OR.Intercept\n", "# # (np.power((data_pred.Temperature-mean_temp),2) / \n", "# # ((np.sum(np.power(data_pred.Temperature,2))) - n*(np.power(mean_temp,2))))\n", "\n", "data_pred['Prob_low']=logit_inv(np.minimum(\n", " data_pred['Temperature'] * conf.low.Temperature + conf.low.Intercept/2,\n", " data_pred['Temperature'] * conf.up.Temperature + conf.low.Intercept/2))\n", "data_pred['Prob_up']=logit_inv(np.maximum(\n", " data_pred['Temperature'] * conf.low.Temperature + conf.up.Intercept/2,\n", " data_pred['Temperature'] * conf.up.Temperature + conf.up.Intercept/2))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VPW9//HXZ9bsCVmEsMkiLoioEFARFLVubbX99bYVW7rc1vKz1ra2Wutu3brYzS5ay/XX5dZbl3pvW8q1FUsFxRYFNwQRREQJBLOH7MnMfH9/zIAxBjIJk8yS9/Px4JGZM+ec+XwZ5p3D93zP95hzDhERySyeZBcgIiKJp3AXEclACncRkQykcBcRyUAKdxGRDKRwFxHJQP2Gu5n9ysyqzWzjAV43M/upmW0zsw1mNivxZYqIyEDEc+T+G+C8g7x+PjAt9mcJ8ItDL0tERA5Fv+HunHsSqD/IKh8C/tNFrQWKzKw8UQWKiMjA+RKwj3HAzh7PK2PLqnqvaGZLiB7dk52dPXvChAmDesNwOIyZDWrbVOOcU1tSTKa0A9SWVOWcw+v1DmrbrVu31jrnyvpbLxHh3tffdp9zGjjnlgJLASoqKtz69esH9YbLly9n/Pjxg9o21VRVVVFenhn/0cmUtmRKO0BtSVWVlZV88IMfHNS2ZvZmPOslYrRMJdDzEHw8sDsB+xURkUFKRLgvAz4dGzVzMtDknHtPl4yIiAyffrtlzOwBYCFQamaVwM2AH8A5dy/wKPB+YBvQBvz7UBUrIiLx6TfcnXMX9/O6A76UsIpERPrhnMPjSd9rMMeOHcvmzZsPuk5WVhbjx4/H7/cP6j0ScUJVRGRYeTweSktLKSoqSssRNF1dXRQVFR3wdeccdXV1VFZWMnny5EG9R/r+6hORES1dgz0eZkZJSQkdHR2D3ofCXUTSUqYG+z6H2j6Fu4hIBlK4i4gMUH5+PieffDIVFRUsXryYtra2AW0/HBdhKtxFRAYoOzubtWvXsn79evx+P/fdd9+7XnfOEYlEklRdlMJdROQQnHrqqWzfvp0333yTWbNmccUVVzBv3jwqKyt5+OGHmTNnDhUVFdxwww3v2u7KK69k1qxZnHXWWdTU1CS8Lg2FFJG09p3HXuPVPS0J3efRY/K49txp/a4XCoVYsWIFZ599NgBbt27l3nvv5a677qKqqoobb7yRNWvWMGrUKC644AL+8pe/cMEFF9Da2sqsWbP44Q9/yK233sott9zCz3/+84S2QUfuIiID1N7ezsknn8z8+fOZMGECn/nMZwCYOHEic+fOBeC5555jwYIFlJWV4fP5WLRoEWvWrAGi4/QvuugiABYvXrx/eSLpyF1E0lo8R9iJtq/PvbecnJz9j6MX78dnKIZ16shdRGQIVFRUsGbNGmprawmHwzz88MMsWLAAgEgkwiOPPALA73//e+bPn5/w99eRu4jIECgvL+eWW27h/PPPxznHueeeu38O99zcXDZt2sTs2bMpLCzkoYceSvj7K9xFRAaourr6PcsOP/xwet+A6KKLLtrft95TZWUlRUVF3HbbbUNWo7plREQykMJdRCQDKdxFJC0NZDRKOjrU9incRSQtNTY2ZmzA75vPPSsra9D70AlVEUk7kUiE2tpaamtrk13KoITDYbKzsw+6zr47MQ2Wwl1E0o6ZpfVR++7du/cPixwq6pYREclACncRkQykcBcRyUAKdxGRDKRwFxHJQAp3EZEMpHAXEclACncRkQykcBcRyUAKdxGRDJR24V69t4N/7sn8GeFERA5F2oX7/c+8xe+3GXc8VU9jRzjZ5YiIpKS0C/crzprGv012vLSng6/+tZpnd7UnuyQRkZSTduHu8RhnjIMfnXsYxdlevv1UPUufa6QrrG4aEZF90i7c95lQ6OfOs8u48KhcHn2tlW+sqKZyb3eyyxIRSQlpG+4Afq/xuROLuPG0Eho6Ily1oobVO9qSXZaISNLFFe5mdp6ZbTGzbWZ2TR+vTzSzJ8zsBTPbYGbvT3ypBzZ7bBY/OvcwJhf5+fHaBn6xroFuddOIyAjWb7ibmRe4GzgfmA5cbGbTe612A/Cwc+5EYBFwT6IL7U9pjpfbzyzlI8fk8djrbVy3soaa1tBwlyEikhLiOXKfC2xzzm13znUBDwIf6rWOAwpijwuB3YkrMX5ej/Hp4wu5Zn4xu5pDfP2xGja83ZmMUkREksr6uxjIzD4KnOecuyT2/FPASc65y3usUw6sAEYBucD7nHPP9bGvJcASgNGjR89+8MEHB1V0U1MTgUDgoOtUtznu3eR4uw0+MtU4c1z0vouppru7G7/fn+wyEiJT2pIp7QC1JVV1dXVRWFg4qG3POOOM55xzFf2tF88NsvtKxN6/ES4GfuOc+6GZnQL8zsxmOOci79rIuaXAUoCKigq3cOHCON7+vZYvX055eflB1ykHfjQxwl1rG3jk9Q7qwjlcNqcIvze1Ar6qqqrftqSLTGlLprQD1JZUVVlZyWDzL17xdMtUAhN6PB/Pe7tdPg88DOCc+xeQBZQmosBDke338M35xSyakc8TO9q46YlaXdUqIiNCPOG+DphmZpPNLED0hOmyXuu8BZwFYGbHEA33mkQWOlgeMxbNKOAb84p5vaGbb6yo4c1GjYcXkczWb7g750LA5cBjwGaio2I2mdmtZnZhbLUrgS+Y2UvAA8BnXYrN7HXqxGy+fVYpoYjj2pU1vLinI9kliYgMmXj63HHOPQo82mvZTT0evwKcmtjSEu+I4gB3nl3G7U/WcdvqOi6tKOLsqbnJLktEJOHS+grVwSjL9fGd95Uxc3SQu9c18uDGvZo+WEQyzogLd4Acv4frTyvhzMk5PLixmXvWNRKOKOBFJHPE1S2TiXwe48tziyjN8fLwpmYaOyJcNa+YoC+1hkqKiAzGiDxy38fM+MRxBSyZXcj63R3csrqWlq5I/xuKiKS4ER3u+7x/Wh5XzhvF1roubvhHjcbCi0jaU7jHzJ+Yw/ULStjdHOa6lbWadExE0prCvYcTy7P41sISGjuiAV/VrIAXkfSkcO9lelmQ284opTPsuG5lje7uJCJpSeHeh6nFAW4/sxQHXL+ylh2arkBE0ozC/QAmFvq548xSfB648R81bG/oSnZJIiJxU7gfxLgCP7efWUbQ5+HmJ2oV8CKSNhTu/SjP93H7maUEfR5uUsCLSJpQuMdhTF404LNiR/DqgxeRVKdwj9OYPB+3nVFKwGvc/EQtO5sU8CKSuhTuA1Ce7+PWM0rxGNz4RC27NExSRFKUwn2AxhX4ufWMUpyDm56o4+0WXegkIqlH4T4IEwr9fGthKZ3hCDc9UUtdm+aiEZHUonAfpMmj/Nx0eil7OyPcvKqWvZ0KeBFJHQr3Q3BkSYDrTyuhujXELavqaOvWdMEikhoU7odoxmFBrj61hB2N3dzxZB2dId3RSUSST+GeABVjs/jqyaN4paaLH/6rXrfsE5GkU7gnyGmH53DJrEKe3dXBPesaddNtEUmqEXsP1aHwgSPz2NsZ4aFNzRRmefj08YXJLklERiiFe4ItmpFPU2eE/9ncwqgsLxcclZfskkRkBFK4J5iZ8YVZhTR2hPnVC02MyvYwf2JOsssSkRFGfe5DwOsxvn5KMceUBbhrbQMb3u5MdkkiMsIo3IdIwGtct6CEsfk+vrumTjNJisiwUrgPobyAhxtPKyHLZ9y2uo5aTVMgIsNE4T7EynJ93HhaKW3dEW5bXaurWEVkWCjch8HkUX6+Ob+Yyr0h7ny6npAuchKRIaZwHyYnjMnii3OKeHFPJ/eu10VOIjK0NBRyGL1vSi5vt4T5wyvNjMnz8dHp+ckuSUQylMJ9mH3iuHzebg1x/4a9jMnzMtWf7IpEJBOpW2aYmRmXzx3FMaUBfrK2ge1N6p4RkcSLK9zN7Dwz22Jm28zsmgOs83Eze8XMNpnZ7xNbZmYJeI1rFxRTkuPlF5ucbtUnIgnXb7ibmRe4GzgfmA5cbGbTe60zDbgWONU5dyxwxRDUmlEKgl5uPK2EsIM7nqqjtUtDJEUkceI5cp8LbHPObXfOdQEPAh/qtc4XgLudcw0AzrnqxJaZmcYV+Fky3di1N8QP/ql54EUkceI5oToO2NnjeSVwUq91jgQws6cBL/At59zfeu/IzJYASwBGjx7NqlWrBlEyhMNhqqqqBrVtqpmaF2LRNB//tbWTnz29m4uOSN/TIN3d3RnxuWRKO0BtSVXhcHjQ+ReveMLd+ljW+xDTB0wDFgLjgafMbIZzrvFdGzm3FFgKUFFR4RYuXDjQegFYvnw55eXlg9o21VRVVfGxWeU0WxPLtrRwdHkB5x2Rm+yyBqWqqiojPpdMaQeoLamqsrKSweZfvOI5TKwEJvR4Ph7Y3cc6f3bOdTvn3gC2EA17idNnji9gdnmQpc81ahZJETlk8YT7OmCamU02swCwCFjWa50/AWcAmFkp0W6a7YksNNN5PcaV84oZX+Dje2vq2N2sETQiMnj9hrtzLgRcDjwGbAYeds5tMrNbzezC2GqPAXVm9grwBPAN51zdUBWdqXL8Hq5bUILHjDuerKNFI2hEZJDiOnvnnHvUOXekc26qc+6O2LKbnHPLYo+dc+7rzrnpzrnjnHMPDmXRmWxMno9vzi9mT0uIH2oEjYgMUvoOzchgMw4L8n8rinhhTye/fakp2eWISBrS3DIp6pypubzZ2M2yLa0cXujnrCnpOYJGRJJDR+4p7HMnFjJzdJBfrG/k1VqNoBGR+CncU5jXY3xjXjGlOV6+u6Zet+kTkbgp3FNcftDD9QtK6Aw5vvNUHZ0hnWAVkf4p3NPAhEI/XztlFNsburl7XYPu4iQi/VK4p4m547L5xHEFPPlmO398tSXZ5YhIilO4p5GPTs/j1AnZ/O6lvTxf1ZHsckQkhSnc04iZ8eWTiji8yM8P/lnPrr3dyS5JRFKUwj3NZPk8XDu/GJ/H+M6aet3kQ0T6pHBPQ6PzfHxjXjG7m0PctbaBiE6wikgvCvc0ddzoIJ8/sZB1uzt44OXmZJcjIilG4Z7G3j8tl7Mm5/CHV5r55872ZJcjIilE4Z7GzIxLK4o4qsTPT59pYEejTrCKSJTCPc35vcY355eQ7TO+81QdzZ06wSoiCveMUJzt5Zr5JdS1h/mB5oAXERTuGeOo0gCXVhTx0tud/PalvckuR0SSTPO5Z5D3Tclle0M3y7a0MGWUn4WTcpJdkogkiY7cM8znTizk2LIA96xrYFt9V7LLEZEkUbhnGJ/HuPrUYgqDXr7zVD2NHZoDXmQkUrhnoMIsL9cuKKa5K8L31tTTHdYJVpGRRuGeoaaMCnD53CI213Zx3/ONyS5HRIaZTqhmsNMOz+GNhm7++GoLU0YFOPcI3WRbZKTQkXuGWzyzgFnlQZY+18imat1kW2SkULhnOK/H+PopxYzO83Ln0/XUtIaSXZKIDAOF+wiQF/Bw3YISuiOO76yppzOkKQpEMp3CfYQYX+Dn66cU80ZDNz97tlE32RbJcAr3EaRibBaLZxaw5q12/nuzbrItkskU7iPMR47JY8HEbP5rw17W7dIc8CKZSuE+wpgZl88tYvIoPz/6VwM7mzQHvEgmUriPQMHYTbaDPuMOzQEvkpEU7iNUWa6Pa+YXU9sW5s6n6wlpDniRjKJwH8GOLg1y2ZwiXq7u5FcvNCW7HBFJIE0/MMKdOTmXt5pC/OnVFiYU+Dh/Wl6ySxKRBNCRu/CpmQXMLg/yH883seFtTVEgkgniCnczO8/MtpjZNjO75iDrfdTMnJlVJK5EGWpej3HlvGLG5fu48+k6djdrigKRdNdvuJuZF7gbOB+YDlxsZtP7WC8f+ArwTKKLlKGX4/dw/WkleMy4/ck6Wro0gkYkncVz5D4X2Oac2+6c6wIeBD7Ux3q3AXcCHQmsT4bRmDwf3zy1mOrWEN/XCBqRtBbPCdVxwM4ezyuBk3quYGYnAhOcc8vN7KoD7cjMlgBLAEaPHs2qVasGXDBAOBymqqpqUNummu7u7pRqSzHwiWnGf27p5CdP7ebiaYaZxbVtqrVlsDKlHaC2pKpwODzo/ItXPOHe1zd7/yGdmXmAHwOf7W9HzrmlwFKAiooKt3DhwriK7G358uWUl5cPattUU1VVlXJt+Ug5tHqa+O/NLUwbU8CFR8U3giYV2zIYmdIOUFtSVWVlJYPNv3jF0y1TCUzo8Xw8sLvH83xgBrDKzHYAJwPLdFI1vX1yZgEnj8/i1y808azmoBFJO/GE+zpgmplNNrMAsAhYtu9F51yTc67UOTfJOTcJWAtc6JxbPyQVy7DwmPG1k0cxJTYHzfaGrmSXJCID0G+4O+dCwOXAY8Bm4GHn3CYzu9XMLhzqAiV5gr7oCJq8gIfbn6yjti2c7JJEJE5xjXN3zj3qnDvSOTfVOXdHbNlNzrllfay7UEftmaM428sNp5XQ3u24/cla2ro1RFIkHegKVenXpCI/V59azFtNGiIpki4U7hKXE8uz+GJFES/s6eSX63WbPpFUp4nDJG5nT83l7dYQj7zSwmG5Pj52bH6ySxKRA1C4y4B84rgCqlvD/NfLeynN8XLG5JxklyQifVC4y4B4zPjy3FE0tEf4+bMNFGd7OH5MVrLLEpFeFO4yYH6vcc38Yq5bWcN319Rzx1mlTBkVGNS+Vu9o5f4NzdS2hSnN8bJ4Zj6nT8pNcMUyVPT5pS6dUJVByQ14uOn0UnIDHm5bXcfbLQOfJnj1jlbuWddETVsYB9S0hblnXROrd7QmvmBJOH1+qU3hLoNWkuPl5tNL6I44blldS3PXwEbQ3L+hmc7wu7fpDDvu39CcyDJliOjzS20KdzkkEwr93HBaCbVtEe7e6GgfwEVOB7riVVfCpgd9fqlN4S6H7OjSIFefWszOZvj2U/V0h+M7gi/N8Q5ouaQWfX6pTeEuCVExNotPH228XN3Jj/5VTziOq1gXz8wn6H33jNJBr7F4psbPpwN9fqlN4S4Jc9Jo43MnFvKvyg7uWdf/VaynT8rlsjmFlOV4MaAsx8tlcwo12iJN6PNLbRoKKQl14VF5tHZFeGhTM7kBD/9+QsFB7+R0+qRchUEa0+eXuhTuknCLZuTT2h1h2ZYWcv3GRTMKkl2SyIijcJeEM4t2z7R1Ox7Y2EzQZ3z4aPXDigwnhbsMCY8ZX5pTRGfI8ZsX9xLwGu+fFt+9WEXk0CncZch4PcYVJ4+iM+xY+lwTfo9x9lT1z4oMB42WkSHl9xpXn1rMiWOC3LOukSfeaEt2SSIjgsJdhlzAa1wzv4TjRgf52bMNrN6hgBcZagp3GRZBn3H9gmKmlwX4yTMKeJGhpnCXYRP0ebjhtBKOjQX8KgW8yJBRuMuwytof8EF+sraBlds1PazIUFC4y7CLHsEXc/yYID97tpG/bVPAiySawl2SIujzcN2CEirGZnHv+kaWbWlJdkkiGUXhLkkT8BrfPLWYU8Zn8asXmnho495+JxsTkfgo3CWp/F7jqnnFnDEphwc2NvPrFxXwIomgK1Ql6bwe48snFZHjN5ZtaaG1K8Jlc4rweg48m6SIHJzCXVKCx4xLZhWSF/Dw0KZmmrsiXHlKMUGfAl5kMNQtIynDzLj4uAK+MKuQdbs6+NaqWpo7478nq4i8Q+EuKecDR+Zx5bxRvFbfxbUra6huDSW7JJG0o3CXlDR/Yg43n15KfXuYbz5ew/aGrmSXJJJWFO6Sso4bHeQ7Z5XhMeO6lbWs392R7JJE0obCXVLa4UV+vn9OGWPzfXz7qTr+d6sudhKJh8JdUl5xtpc7zixldnkW//F8E0ufayQc0Vh4kYOJK9zN7Dwz22Jm28zsmj5e/7qZvWJmG8xspZkdnvhSZSTL9nu4Zn4xHz46j0dfa+XW1XW0dGkkjciB9BvuZuYF7gbOB6YDF5vZ9F6rvQBUOOdmAo8Adya6UBGvx/jsCYV8eW4Rm2o6uWpFNW81dSe7LJGUFM+R+1xgm3Nuu3OuC3gQ+FDPFZxzTzjn9k3OvRYYn9gyRd5x1pRcbj+zlM6Q4+rHa/jXzvZklySScuK5QnUcsLPH80rgpIOs/3ngr329YGZLgCUAo0ePZtWqVfFV2Us4HKaqqmpQ26aa7u5utWUQCoGrT3D8cpPje0/Xc84EuHCy4bVDv6JVn0lqyqS2hMPhQedfvOIJ976+LX2ezTKzxUAFcHpfrzvnlgJLASoqKtzChQvjq7KX5cuXU15ePqhtU01VVZXaMkjlwPcnOO57vpHHXm+jqtPPlfOKKcryHtJ+9ZmkpkxqS2VlJYPNv3jF0y1TCUzo8Xw8sLv3Smb2PuB64ELnXGdiyhM5OL/X+OKcUXzlpCK21HXxtb9V8/Lb+ucnEk+4rwOmmdlkMwsAi4BlPVcwsxOBXxIN9urElylycGdOzuXOsw8j2+/h5lW1PLRxr4ZLyojWb7g750LA5cBjwGbgYefcJjO71cwujK32fSAP+IOZvWhmyw6wO5EhM6nIzw/PKWPBxGwe2NjMzatqqW0LJ7sskaSIa8pf59yjwKO9lt3U4/H7ElyXyKA8u6udTdXRbpmN1V1c/ujbnD0lm7WVndS2hSnN8bJ4Zj6nT8pN+Huv3tHK/Ruah/x94nHvunpWbG/nazNCfOnJXZwzJZtL5xQnpRZJDs3nLhlj9Y5W7lnXRGf4ne6YjpDjL1vb9j+vaQtzz7omgIQGb+/3Hqr3ice96+r52+vvDA+NOGLP6xXwI4imH5CMcf+G5ncF+4F0hh33b2ge8vceiveJx4rtfY/7P9ByyUwKd8kYA+lfr0lwX/yB3jsZff4HOo+s88sji8JdMkZpTvzj2z0GL+1J3BTCB3rvgdSUKAe69axuSTuyKNwlYyyemU/Q++4E8xr4ev0r93ugMOjh5lV1/PSZBvZ2HvrRdV/vHfQai2fmH/K+B+qcKdkDWi6ZSSdUJWPsO3HZe8RKX8tOmZDDQxv38qdXW1i3q4PPnlDAmZNzsEFOX3Cg907GaJnoSdP6/X3sHkOjZUYghbtklNMn5fYZqH0t+9TxhSw4PIdfrm/kZ8828vftbSyZXURWgt87GS6dU8ylc6KX7P/PRZlxyb4MjLplZESbVOTnjrNK+dKcInY1h7hyRTUPvRahuVNzxUt6U7jLiOcx4+ypudzzgdGcd0Quq3fDF/93D8u3thDSEBNJUwp3kZi8gIcls4u4frYxZVSA+55v4it/rWZtZTvOKeQlvSjcRXoZl2fcsrCE6xYU4zH47pp6rl1Zu39aA5F0oBOqIn0wM+aOy2Z2eRYr32jjgZf3cv0/ajlxTJBPzizgiOJAsksUOSiFu8hBeD3GOVNzOf3wbP66rZX/fqWZq1bUUDE2i0Uz8hXykrIU7iJxCPo8fPjofM6Zmsv/bm3lz1uauWpFB7PKg3x0ej7Ty4LJLlHkXRTuIgOQ4/fwsWPz+cCRuTz6Wit/2dLCdStrmV4W4MNH51ExNgtPAu7jKnKoFO4ig5Dj9/DR6flccGQuK15vY9mWFr79VD3jC3xccGQeCydlE+w974HIMFK4ixyCoM/DBUflcf60XP65s50/vdrCL9Y3cv+GJs6emsu5U3MZnaevmQw//asTSQCfxzjt8BwWTMzmlZoulm9t4U+vtvDHzS3MHpvFeUfkcuKYIF5NzSjDROEukkBmxrGHBTn2sCA1rSFWvN7G49tbuf3JDkpzvJw1OYczJ+foaF6GnP6FiQyRslwfn5xZwMePzWfd7g5WvN7Kw5uaeWhTM8cdFuCMyTmcMj6bbL/65iXxFO4iQ8zvNeZNyGbehGyqW0M88UYb/3ijjZ8+08i965s4aVwWCw7P5sQxWfi96raRxFC4iwyjw3J9XDQjejS/pa6LVTvaefqtNp56q50cv3HSuGxOmZDFCWOyCCjo5RAo3EWSwMw4ujTI0aVBLplVyIa3O3nqzXae3dXOEzvayPYZs8qzmDsui9ljs8gLqOtGBkbhLpJkPk80yGeVZ9EdLuLl6k7+tbOddbs7eHpnOx6D6WUBKsZG15lQ4Bv0HaNk5FC4i6QQv/edoI84x9a6Ltbv7mD97g5+8+JefvPiXkqyvZwwJsgJY4IcNzpIUdbw34RbUp/CXSRFeXp03SyeWUhNa4gX93TyfFUHayvbWflGGwCHF/o4bnSQGYcFOaY0QKHCXlC4i6SNslwfZ0/1cfbUXMIRx+sN3by0p5ON1Z2seL2N5VtbARhf4OOY0gBHlwYocY4xzqkbZwRSuIukIa/HOLIkwJElAT52bD7dYce2+i5eqenilZpO/rmznce3R4/s81+qYlpJgGnFAY4o9nNEcYBR2Tq6z3QKd5EM4Pcax5QFOaYsyL+RT8Q5du0Nsfb1GqpD2Wyt6+LFPc3suyXsqCwPU4v9TC4KMHmUn0lFfkbnejU9QgZRuItkII8ZEwr9+MqN8vJRAHSEImxv6GZbfTdvNHTxekM3z1e9E/gBrzGhwMfEQj8TC31MKPQzvsBHWY5CPx0p3EVGiCyfh+llwXfdWKQr7NjZ1M2Oxm7ebOrmzcYQL+7p4Ikdkf3rBLxQnudjbH70T3m+j/K86M+iLI/mr09RCneRESzgNaYWB5ja63aBzZ0RKvd2U7k3ROXeELubQ7zVFOLZXR2EXc/to1fdHpbrZXSuj7JcL4fleinLiT4uDHp01J8kCncReY/8oGd/H35P4Yijpi1MVXOIPS0h3m4Ns6clRHVrmK11bbR0uXet7zUozvZSkuON/sz2UJztZVS2l1FZ0ceFWR7yAx6N6EkwhbuIxM3rMcbk+RhzgCmL27ojVLeGqWkNUdsWprYtTF17mLq2CDsau3m+KkxHyL1nO58HCoMeCrOiR/sFQQ+FWR4Kgl4KYs+7mh1d2d3kBz3kBTyae6cfCncRSZgcv4dJRR4mFfkPuE57d4T69jD17REaOsI0dkRo6gjT0BFhb2f0+a7mEHs7I+/9RfBS9f6HAa+RFzBy/dGwz409zvF7yPEbuQEP2T4jx+8h229k+/b9NLL90deCPsvYcwZxhbuZnQf8BPAC9znnvtvr9SDwn8BsoA64yDm3I7GlimSu1TtauX9DM7VtYUpzvCyemc/K7a1sqO6Ss4KhAAAKDklEQVTev87Mw/ycNSX3PesB71m2uaaTFdvb+dqMEF96chfnTMnm0jnFcb1vX/s7fVJu3HXve++IA4/xnvfO9nsY5/ewrb7/9/7CrDyOH5PN3s4wb1bV4s8rorkzwsbqTp6v6qS+PUJLVwTw0Rk2djaFaOuO0Nrt9o8C6k/Qa2TFgj7LZwS97zwOeD0EvdHbKQa80V8ovf/4PbGf+5dF5wvye6LL/J7oUNXosuhrLs7aDkW/4W5mXuBu4GygElhnZsucc6/0WO3zQINz7ggzWwR8D7hoKAoWyTSrd7Ryz7omOmNnKmvawvx4beN71ttQ3c2G6neW17SF+dmzjTjH/pOcNW1h7lrbSM/siDj42+vtQP27Qrav9/3pM42YQSjyzv7uWdcE8J6A72v7RL/30uf2ctkc4/RJuWR3GOXlOaze0cr63Z37t+0KQ3VrmMvmFO6v0TlHV9jR1u1o647QHnK0xx53hBwdIUd7KEJ7t6M95OgMRZd3hh2dsdcb2iN0hsN0hqL76gxHf8b7S+NgPj4VLjj03RxUPEfuc4FtzrntAGb2IPAhoGe4fwj4VuzxI8DPzcycG47fTyLp7f4NzfuDaqBCkfcuO9CeVmxv59I5B3/fsHvvDjrDjvs3NL8n3Pvafjjeu69te69nFj36DvpI6NW4zjlCkegQ0u5I9BdBKOLoCkN3xNEd+wXQHXlnvVBseXcEQpHo8/HevQmr6UCsv/w1s48C5znnLok9/xRwknPu8h7rbIytUxl7/npsndpe+1oCLIk9PQrYMpiig8HglKysrI7BbJtqurq6cgKBQFuy60iETGnLcLejPbusfKj2HW5rwptTuP95dntN1WDft+e2h7r9YLbd97kcbNveNaaqjo6OrM7Ozu2D3Pxw51xZfyvFc+Te19mG3r8R4lkH59xSYGkc73nwgszWd3R0VBzqflKBma1va2tTW1JIprQDom0JNVVnTFsy6XNxzg1pW+K5vUslMKHH8/HA7gOtY2Y+oBCoT0SBIiIycPGE+zpgmplNNrMAsAhY1mudZcBnYo8/CvxD/e0iIsnTb7eMcy5kZpcDjxEdCvkr59wmM7sVWO+cWwb8P+B3ZraN6BH7oqEsmgR07aQQtSX1ZEo7QG1JVUPeln5PqIqISPrRLdVFRDKQwl1EJAOlfLibWZaZPWtmL5nZJjO7JbZ8spk9Y2avmdlDsZO9Kc/MvGb2gpktjz1P13bsMLOXzexFM1sfW1ZsZo/H2vK4mY1Kdp3xMLMiM3vEzF41s81mdko6tsXMjop9Hvv+7DWzK9K0LV+Lfd83mtkDsRxI1+/KV2Pt2GRmV8SWDflnkvLhDnQCZzrnjgdOAM4zs5OJTnHwY+fcNKCB6BQI6eCrwOYez9O1HQBnOOdO6DFe9xpgZawtK2PP08FPgL85544Gjif6+aRdW5xzW2KfxwlE53lqA/5ImrXFzMYBXwEqnHMziA7k2DetSVp9V8xsBvAFolf6Hw980MymMRyfiXMubf4AOcDzwElALeCLLT8FeCzZ9cVR//jYB3kmsJzoxV9p145YrTuA0l7LtgDlscflwJZk1xlHOwqAN4gNLkjntvSq/xzg6XRsCzAO2AkUEx3Rtxw4Nx2/K8DHiE62uO/5jcDVw/GZpMOR+76ujBeBauBx4HWg0TkXiq1SSfQfRKq7i+gHu29GkBLSsx0QvQJ5hZk9F5tWAmC0c64KIPbzsKRVF78pQA3w61h32X1mlkt6tqWnRcADscdp1Rbn3C7gB8BbQBXQBDxHen5XNgKnmVmJmeUA7yd6weeQfyZpEe7OubCL/ldzPNH/3hzT12rDW9XAmNkHgWrn3HM9F/exakq3o4dTnXOzgPOBL5nZackuaJB8wCzgF865E4FWUrzboj+xvugLgT8ku5bBiPU/fwiYDIwFcon+O+st5b8rzrnNRLuTHgf+BrwEhA66UYKkRbjv45xrBFYBJwNFsakOoO8pEVLNqcCFZrYDeJBo18xdpF87AHDO7Y79rCbarzsXeNvMygFiP6sPvIeUUQlUOueeiT1/hGjYp2Nb9jkfeN4593bsebq15X3AG865GudcN/A/wDzS97vy/5xzs5xzpxG9yPM1huEzSflwN7MyMyuKPc4m+sFvBp4gOtUBRKc++HNyKoyPc+5a59x459wkov9l/odz7pOkWTsAzCzXzPL3PSbav7uRd09DkRZtcc7tAXaa2VGxRWcRnc467drSw8W80yUD6deWt4CTzSzHzIx3PpO0+64AmNlhsZ8TgY8Q/WyG/DNJ+StUzWwm8FuiZ8w9wMPOuVvNbArRI+Bi4AVgsXOuM3mVxs/MFgJXOec+mI7tiNX8x9hTH/B759wdZlYCPAxMJPoF/ZhzLuUnkDOzE4D7gACwHfh3Yv/WSL+25BA9GTnFOdcUW5Z2n0tsyPNFRLswXgAuIdrHnlbfFQAze4ro+bVu4OvOuZXD8ZmkfLiLiMjApXy3jIiIDJzCXUQkAyncRUQykMJdRCQDKdxFRDJQPDfIFhlWsWFiK2NPxwBholMEAMx1znUlpbCDMLPPAY/Gxs2LJJ2GQkpKM7NvAS3OuR+kQC1e51z4AK+tAS53zr04gP35esyVIpJQ6paRtGJmn7Ho/P4vmtk9ZuYxM5+ZNZrZ983seTN7zMxOMrPVZrbdzN4f2/YSM/tj7PUtZnZDnPu93cyeBeaa2S1mti42P/e9FnUR0emoH4ptHzCzyh5XVp9sZn+PPb7dzH5pZo8TnazMZ2Y/ir33BjO7ZPj/ViUTKdwlbcTmxv4/wLzYRHI+3rkZeyGwIjaZWRfwLaKXrX8MuLXHbubGtpkFfMLMTohjv8875+Y65/4F/MQ5Nwc4Lvbaec65h4AXgYtcdD71/rqNTgQucM59ClhCdEK5ucAcopOwTRzM349IT+pzl3TyPqIBuD465QjZRC+1B2h3zj0ee/wy0OScC5nZy8CkHvt4zDnXAGBmfwLmE/0eHGi/Xbwz1QLAWWb2DSALKCU6Fe1fB9iOPzvnOmKPzwGOMbOev0ymEb0kXWTQFO6STgz4lXPuxnctjM4U2PNoOUL0Dl77Hvf8d977JJPrZ7/tLnZiKjZvy8+BWc65XWZ2O9GQ70uId/5n3Hud1l5tusw5txKRBFK3jKSTvwMfN7NSiI6qGUQXxjkWvWdqDtE5w58ewH6zif6yqI3NivlvPV5rBvJ7PN9B9FZ39Fqvt8eAy/ZNZWvR+6BmD7BNIu+hI3dJG865l2OzBf7dzDxEZ9m7lIHN670G+D0wFfjdvtEt8ezXOVdnZr8lOr3xm8AzPV7+NXCfmbUT7df/FvAfZrYHePYg9fyS6MyAL8a6hKqJ/tIROSQaCikjRmwkygzn3BXJrkVkqKlbRkQkA+nIXUQkA+nIXUQkAyncRUQykMJdRCQDKdxFRDKQwl1EJAP9fwhNYzzKNBjLAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "## http://markthegraph.blogspot.com/2015/05/using-python-statsmodels-for-ols-linear.html\n", "data_pred.plot(x=\"Temperature\",y=\"Prob\",kind=\"line\",ylim=[0,1])\n", "plt.fill_between(data_pred.Temperature,data_pred.Prob_low,data_pred.Prob_up,color='#888888', alpha=0.4)\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "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.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }