{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this document we reperform some of the analysis provided in \n", "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n", "\n", "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Technical information on the computer on which the analysis is run" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#!pip install statsmodels \n", "#!pip install seaborn" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.7.6 (default, Jan 8 2020, 19:59:22) \n", "[GCC 7.3.0]\n", "uname_result(system='Linux', node='H2-SCHMIDER', release='5.3.0-46-generic', version='#38~18.04.1-Ubuntu SMP Tue Mar 31 04:17:56 UTC 2020', machine='x86_64', processor='x86_64')\n", "IPython 7.13.0\n", "IPython.core.release 7.13.0\n", "html 6.0.3\n", "_csv 1.0\n", "_ctypes 1.1.0\n", "_curses b'2.2'\n", "decimal 1.70\n", "argparse 1.1\n", "backcall 0.1.0\n", "csv 1.0\n", "ctypes 1.1.0\n", "cycler 0.10.0\n", "dateutil 2.8.1\n", "decimal 1.70\n", "decorator 4.4.2\n", "distutils 3.7.6\n", "ipykernel 5.1.4\n", "ipykernel._version 5.1.4\n", "ipython_genutils 0.2.0\n", "ipython_genutils._version 0.2.0\n", "jedi 0.16.0\n", "json 2.0.9\n", "jupyter_client 6.1.2\n", "jupyter_client._version 6.1.2\n", "jupyter_core 4.6.3\n", "jupyter_core.version 4.6.3\n", "kiwisolver 1.2.0\n", "logging 0.5.1.2\n", "matplotlib 3.2.1\n", "matplotlib.backends.backend_agg 3.2.1\n", "numpy 1.18.2\n", "numpy.core 1.18.2\n", "numpy.core._multiarray_umath 3.1\n", "numpy.lib 1.18.2\n", "numpy.linalg._umath_linalg b'0.1.5'\n", "pandas 1.0.3\n", "parso 0.6.2\n", "patsy 0.5.1\n", "patsy.version 0.5.1\n", "pexpect 4.8.0\n", "pickleshare 0.7.5\n", "platform 1.0.8\n", "prompt_toolkit 3.0.4\n", "ptyprocess 0.6.0\n", "pygments 2.6.1\n", "pyparsing 2.4.6\n", "pytz 2019.3\n", "re 2.2.1\n", "scipy 1.4.1\n", "scipy._lib._uarray 0.5.1+5.ga864a57.scipy\n", "scipy._lib.decorator 4.0.5\n", "scipy._lib.six 1.2.0\n", "scipy.integrate._dop b'$Revision: $'\n", "scipy.integrate._ode $Id$\n", "scipy.integrate._odepack 1.9 \n", "scipy.integrate._quadpack 1.13 \n", "scipy.integrate.lsoda b'$Revision: $'\n", "scipy.integrate.vode b'$Revision: $'\n", "scipy.interpolate._fitpack 1.7 \n", "scipy.interpolate.dfitpack b'$Revision: $'\n", "scipy.linalg 0.4.9\n", "scipy.linalg._fblas b'$Revision: $'\n", "scipy.linalg._flapack b'$Revision: $'\n", "scipy.linalg._flinalg b'$Revision: $'\n", "scipy.ndimage 2.0\n", "scipy.optimize._cobyla b'$Revision: $'\n", "scipy.optimize._lbfgsb b'$Revision: $'\n", "scipy.optimize._minpack 1.10 \n", "scipy.optimize._nnls b'$Revision: $'\n", "scipy.optimize._slsqp b'$Revision: $'\n", "scipy.optimize.minpack2 b'$Revision: $'\n", "scipy.signal.spline 0.2\n", "scipy.sparse.linalg.eigen.arpack._arpack b'$Revision: $'\n", "scipy.sparse.linalg.isolve._iterative b'$Revision: $'\n", "scipy.special.specfun b'$Revision: $'\n", "scipy.stats.mvn b'$Revision: $'\n", "scipy.stats.statlib b'$Revision: $'\n", "seaborn 0.10.0\n", "seaborn.external.husl 2.1.0\n", "six 1.14.0\n", "statsmodels 0.11.1\n", "statsmodels.__init__ 0.11.1\n", "statsmodels.api 0.11.1\n", "statsmodels.tools.web 0.11.1\n", "traitlets 4.3.3\n", "traitlets._version 4.3.3\n", "urllib.request 3.7\n", "zlib 1.0\n", "zmq 18.1.1\n", "zmq.sugar 18.1.1\n", "zmq.sugar.version 18.1.1\n" ] } ], "source": [ "def print_imported_modules():\n", " import sys\n", " for name, val in sorted(sys.modules.items()):\n", " if(hasattr(val, '__version__')): \n", " print(val.__name__, val.__version__)\n", "# else:\n", "# print(val.__name__, \"(unknown version)\")\n", "def print_sys_info():\n", " import sys\n", " import platform\n", " print(sys.version)\n", " print(platform.uname())\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "import seaborn as sns\n", "\n", "print_sys_info()\n", "print_imported_modules()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading and inspecting data\n", "Let's start by reading data." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "file downloaded on Tue Apr 14 10:06:35 2020\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
56/18/82672500
68/30/836731000
711/28/836701000
82/03/846572001
94/06/846632001
108/30/846702001
1110/05/846782000
1211/08/846672000
131/24/856532002
144/12/856672000
154/29/856752000
166/17/856702000
177/2903/856812000
188/27/856762000
1910/03/856792000
2010/30/856752002
2111/26/856762000
221/12/866582001
\n", "
" ], "text/plain": [ " Date Count Temperature Pressure Malfunction\n", "0 4/12/81 6 66 50 0\n", "1 11/12/81 6 70 50 1\n", "2 3/22/82 6 69 50 0\n", "3 11/11/82 6 68 50 0\n", "4 4/04/83 6 67 50 0\n", "5 6/18/82 6 72 50 0\n", "6 8/30/83 6 73 100 0\n", "7 11/28/83 6 70 100 0\n", "8 2/03/84 6 57 200 1\n", "9 4/06/84 6 63 200 1\n", "10 8/30/84 6 70 200 1\n", "11 10/05/84 6 78 200 0\n", "12 11/08/84 6 67 200 0\n", "13 1/24/85 6 53 200 2\n", "14 4/12/85 6 67 200 0\n", "15 4/29/85 6 75 200 0\n", "16 6/17/85 6 70 200 0\n", "17 7/2903/85 6 81 200 0\n", "18 8/27/85 6 76 200 0\n", "19 10/03/85 6 79 200 0\n", "20 10/30/85 6 75 200 2\n", "21 11/26/85 6 76 200 0\n", "22 1/12/86 6 58 200 1" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")\n", "data_url ='https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/raw/7e01583b99527ad27cbeae0f9d2085fe8f2f1d15/data/shuttle.csv?inline=false'\n", "data_file='/home/aschmide/Documents/formation_MOOC_RR/mooc-rr/module4/shuttle.csv'\n", "import os\n", "import urllib.request\n", "\n", "if not os.path.exists(data_file):\n", " urllib.request.urlretrieve(data_url, data_file)\n", " \n", "data = pd.read_csv(data_file)\n", "\n", "import time\n", "modtime=os.path.getmtime(data_file)\n", "modtime=time.ctime(modtime)\n", "print('file downloaded on',modtime)\n", "data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAYAElEQVR4nO3dfZRddX3v8fdnkiEPJEJMvCMlQbBErlRyEYaHXGqd1IcGu0xq4Sp0IeoV4r2X3Cq6KtTrolyqaxVapVcvVVNKFbQOD7E4rWlRiqPYK5BEQoDQ4BQwTMAAYwgZCMlM5nv/OHv0zHBmZp/J7HPmnN/ntVZWzt77d/b5/s6eM5/ZD+e3FRGYmVm6WupdgJmZ1ZeDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscYUFgaQbJD0j6aExlkvSFyT1SNoq6ZSiajEzs7EVuUfwVWDlOMvPBpZm/9YAXyqwFjMzG0NhQRARPwR+MU6T1cCNUXIPcKSko4qqx8zMKptZx9c+GniybLo3m/f06IaS1lDaa2DOnDmnLlmypCYF5jU0NERLS/OdbmnWfkHz9s39ajy16tujjz76XES8ptKyegZBbhGxDlgH0N7eHps2bapzRSN1d3fT0dFR7zKmXLP2C5q3b+5X46lV3yT9bKxl9YzYnUD5n/aLs3lmZlZD9QyCLuDC7OqhM4E9EfGKw0JmZlaswg4NSfom0AEsktQL/AnQChARXwY2AO8CeoCXgA8VVYuZmY2tsCCIiPMnWB7AJUW9vpmZ5dOcp+HNzCw3B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiSs0CCStlLRdUo+kyyssP0bS9yXdL2mrpHcVWY+Zmb1SYUEgaQZwHXA2cCJwvqQTRzX7NHBLRLwZOA/4q6LqMTOzyorcIzgd6ImIxyLiANAJrB7VJoBXZY+PAJ4qsB4zM6tAEVHMiqVzgZURcVE2/X7gjIhYW9bmKOC7wALgcODtEbG5wrrWAGsA2traTu3s7Cyk5snq7+9n3rx59S5jyjVrv6B5++Z+NZ5a9W3FihWbI6K90rKZhb/6+M4HvhoRn5O0HLhJ0psiYqi8UUSsA9YBtLe3R0dHR+0rHUd3dzfTraap0Kz9gubtm/vVeKZD34o8NLQTWFI2vTibV+7DwC0AEfFjYDawqMCazMxslCKDYCOwVNJxkg6jdDK4a1SbHcDbACS9kVIQPFtgTWZmNkphQRARg8Ba4A7gEUpXBz0s6SpJq7JmnwAulvQA8E3gg1HUSQszM6uo0HMEEbEB2DBq3hVlj7cBZxVZg5mZjc/fLDYzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscYUGgaSVkrZL6pF0+Rht3itpm6SHJf1dkfWYmdkrzSxqxZJmANcB7wB6gY2SuiJiW1mbpcAfA2dFxG5J/6GoeszMrLJcewSSTprEuk8HeiLisYg4AHQCq0e1uRi4LiJ2A0TEM5N4HTMzOwSKiIkbSXcDs4CvAt+IiD05nnMusDIiLsqm3w+cERFry9rcDjwKnAXMAK6MiH+usK41wBqAtra2Uzs7OyfuWQ319/czb968epcx5Zq1X9C8fXO/Gk+t+rZixYrNEdFeaVmuQ0MR8ZbsMM5/BTZLug/424j43iHWNhNYCnQAi4EfSjopIp4f9frrgHUA7e3t0dHRcYgvO7W6u7uZbjVNhWbtFzRv39yvxjMd+pb7ZHFE/BT4NHAZ8FbgC5L+TdLvj/GUncCSsunF2bxyvUBXRAxExOOU9g6W5q3JzMwOXd5zBMskXQs8Avw28O6IeGP2+NoxnrYRWCrpOEmHAecBXaPa3E5pbwBJi4A3AI9V2wkzM5u8vFcNfRG4HvhUROwbnhkRT0n6dKUnRMSgpLXAHZSO/98QEQ9LugrYFBFd2bJ3StoGHAT+KCL6DqE/ZmZWpbxB8LvAvog4CCCpBZgdES9FxE1jPSkiNgAbRs27ouxxAB/P/pmZWR3kPUdwJzCnbHpuNs/MzBpc3iCYHRH9wxPZ47nFlGRmZrWUNwhelHTK8ISkU4F947Q3M7MGkfccwceAWyU9BQh4LfC+wqoyM7OayfuFso2S/iNwQjZre0QMFFeWmZnVSjWDzp0GHJs95xRJRMSNhVRlZmY1kysIJN0E/DqwhdL1/gABOAjMzBpc3j2CduDEyDNCnZmZNZS8Vw09ROkEsZmZNZm8ewSLgG3ZqKP7h2dGxKpCqjIzs5rJGwRXFlmEmZnVT97LR38g6XXA0oi4U9JcSgPJmZlZg8s7DPXFwG3AV7JZR1MaQtrMzBpc3pPFl1C6neQL8Mub1PhG82ZmTSBvEOzPbkAPgKSZlL5HYGZmDS5vEPxA0qeAOZLeAdwK/ENxZZmZWa3kDYLLgWeBB4GPULrZTMU7k5mZWWPJe9XQEPDX2T8zM2sieccaepwK5wQi4vVTXpGZmdVUNWMNDZsN/Bfg1VNfjpmZ1VqucwQR0Vf2b2dE/CWlG9qbmVmDy3to6JSyyRZKewjV3MvAzMymqby/zD9X9ngQeAJ475RXY2ZmNZf3qqEVRRdiZmb1kffQ0MfHWx4Rn5+acszMrNaquWroNKArm343cB/w0yKKMjOz2skbBIuBUyJiL4CkK4HvRMQFRRVmZma1kXeIiTbgQNn0gWyemZk1uLx7BDcC90n6+2z694CvFVOSmZnVUt6rhj4r6Z+At2SzPhQR9xdXlpmZ1UreQ0MAc4EXIuL/AL2SjiuoJjMzq6G8t6r8E+Ay4I+zWa3A14sqyszMaifvHsF7gFXAiwAR8RQwv6iizMysdvIGwYGICLKhqCUdXlxJZmZWS3mD4BZJXwGOlHQxcCe+SY2ZWVOYMAgkCbgZuA1YD5wAXBERX8zx3JWStkvqkXT5OO3OkRSS2sdqY2ZmxZjw8tGICEkbIuIk4Ht5VyxpBnAd8A6gF9goqSsito1qNx/4KHBvVZWbmdmUyHto6CeSTqty3acDPRHxWEQcADqB1RXa/SlwNfByles3M7MpkPebxWcAF0h6gtKVQ6K0s7BsnOccDTxZNt2breeXshveLImI70j6o7FWJGkNsAagra2N7u7unGXXRn9//7SraSo0a7+gefvmfjWe6dC3cYNA0jERsQP4nal+YUktwOeBD07UNiLWAesA2tvbo6OjY6rLOSTd3d1Mt5qmQrP2C5q3b+5X45kOfZtoj+B2SqOO/kzS+og4p4p17wSWlE0vzuYNmw+8CegunY/mtUCXpFURsamK1zEzs0Mw0TkClT1+fZXr3ggslXScpMOA8/jV/QyIiD0RsSgijo2IY4F7AIeAmVmNTRQEMcbjCUXEILAWuAN4BLglIh6WdJWkVdWVaWZmRZno0NB/kvQCpT2DOdlj+NXJ4leN9+SI2ABsGDXvijHaduSq2MzMptS4QRARM2pViJmZ1Uc1w1CbmVkTchCYmSXOQWBmljgHgZlZ4pIJgr7+/Tzw5PP09e+vdylmVoW+/v3sGzjoz26BkgiCb2/ZyVlX38UF19/LWVffRdeWnRM/yczqbviz+/izL/qzW6CmD4K+/v1ctn4rLw8MsXf/IC8PDPHJ9Vv914XZNFf+2T0Y4c9ugZo+CHp376O1ZWQ3W1ta6N29r04VmVke/uzWTtMHweIFcxgYGhoxb2BoiMUL5tSpIjPLw5/d2mn6IFg4bxbXnLOM2a0tzJ81k9mtLVxzzjIWzptV79LMbBzln90Zkj+7Bcp7Y5qGturkoznr+EX07t7H4gVz/INk1iCGP7v3/fhH/Ouq3/RntyBJBAGU/rrwD5FZ41k4bxZzWmf481ugpj80ZGZm43MQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiXMQmJklzkFgZpY4B4GZWeIcBGZmiSs0CCStlLRdUo+kyyss/7ikbZK2SvoXSa8rsh4zM3ulwoJA0gzgOuBs4ETgfEknjmp2P9AeEcuA24BriqrHzMwqK3KP4HSgJyIei4gDQCewurxBRHw/Il7KJu8BFhdYj5mZVaCIKGbF0rnAyoi4KJt+P3BGRKwdo/3/BX4eEZ+psGwNsAagra3t1M7OzkJqnqz+/n7mzZtX7zKmXLP2C5q3b+5X46lV31asWLE5ItorLZtZ+KvnIOkCoB14a6XlEbEOWAfQ3t4eHR0dtSsuh+7ubqZbTVOhWfsFzds396vxTIe+FRkEO4ElZdOLs3kjSHo78L+At0bE/gLrMTOzCoo8R7ARWCrpOEmHAecBXeUNJL0Z+AqwKiKeKbAWMzMbQ2FBEBGDwFrgDuAR4JaIeFjSVZJWZc3+HJgH3Cppi6SuMVZnZmYFKfQcQURsADaMmndF2eO3F/n6jayvfz+9u/exeMEcFs6bNWVtG0mz9qsoPbv2svulAXp27eX4tvn1LscayLQ4WWwjfXvLTi5bv5XWlhYGhoa45pxlrDr56ENu20iatV9FueL2B7nxnh184qRBLr32h1y4/BiuWn1SvcuyBuEhJqaZvv79XLZ+Ky8PDLF3/yAvDwzxyfVb6et/5Xn0ato2kmbtV1F6du3lxnt2jJh344930LNrb50qskbjIJhmenfvo7Vl5GZpbWmhd/e+Q2rbSJq1X0XZ8uTzVc03G81BMM0sXjCHgaGhEfMGhoZYvGDOIbVtJM3ar6KcvOTIquabjeYgmGYWzpvFNecsY3ZrC/NnzWR2awvXnLOs4snSato2kmbtV1GOb5vPhcuPGTHvwuXH+ISx5eaTxdPQqpOP5qzjF+W6Yqaato2kWftVlKtWn8SFZx7Lg5vv4c5Lz3QIWFUcBNPUwnmzcv/yq6ZtI2nWfhXl+Lb59M5tdQhY1XxoyMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS1yhQSBppaTtknokXV5h+SxJN2fL75V0bJH1mFWrr38/Dzz5PH39+8dtt+nxPj7/3e1serxvytZZbdueXXvZ/dIAPbv2Tti2GkXVW83r7xs4mPs9uG3Tk033HhS5XoCZU77GjKQZwHXAO4BeYKOkrojYVtbsw8DuiDhe0nnA1cD7iqrJrBrf3rKTy9ZvpbWlhYGhIa45ZxmrTj76Fe0uuP4eftRTCoAv3NXDW45fyE0XnXlI66y27RW3P8iN9+zgEycNcum1P+TC5cdw1eqTJtnz4uut9vX/8I0DXHr1Xbneg2HN8h4Uud5hRe4RnA70RMRjEXEA6ARWj2qzGvha9vg24G2SVGBNZrn09e/nsvVbeXlgiL37B3l5YIhPrt/6ir/GNj3e98sQGHZ3T1/FPYO866y2bc+uvSN+AQLc+OMdh/xXcVH1Tub1D0Yk+R4Uud5yiogpW9mIFUvnAisj4qJs+v3AGRGxtqzNQ1mb3mz637M2z41a1xpgTTZ5ArC9kKInbxHw3IStGk+z9gsm6Jta58ydueCoN6ilZcbwvBgaOji4++lHY2DfS8PzZsxf9GszDj/yqNHPP/ji808f3PvcU5NZZ7VtW+YesXDmq15zLMDBl/YwY+4RAAy+8OwTQy/tmfhY1SG+B9W2nczrD/crz3tQrkHegyn5WczhdRHxmkoLCjs0NJUiYh2wrt51jEXSpohor3cdU61Z+wXN2zdJmwb3PON+NZDp8LNY5KGhncCSsunF2byKbSTNBI4AJp3eZmZWvSKDYCOwVNJxkg4DzgO6RrXpAj6QPT4XuCuKOlZlZmYVFXZoKCIGJa0F7gBmADdExMOSrgI2RUQX8DfATZJ6gF9QCotGNG0PWx2iZu0XNG/f3K/GU/e+FXay2MzMGoO/WWxmljgHgZlZ4hwEVZL0hKQHJW2RtCmbd6Wkndm8LZLeVe86J0PSkZJuk/Rvkh6RtFzSqyV9T9JPs/8X1LvOao3Rr4bfZpJOKKt/i6QXJH2s0bfZOP1qhm12qaSHJT0k6ZuSZmcX1NybDbVzc3ZxTW3r8jmC6kh6Amgv/9KbpCuB/oj4i3rVNRUkfQ24OyKuz34Y5wKfAn4REX+WjRe1ICIuq2uhVRqjXx+jCbbZsGxIl53AGcAlNPg2GzaqXx+igbeZpKOBHwEnRsQ+SbcAG4B3Ad+KiE5JXwYeiIgv1bI27xEYAJKOAH6L0pVcRMSBiHiekcOAfA34vfpUODnj9KvZvA3494j4GQ2+zUYp71czmAnMyb43NRd4GvhtSkPsQJ22l4OgegF8V9LmbOiLYWslbZV0Q6PtimeOA54F/lbS/ZKul3Q40BYRT2dtfg601a3CyRmrX9D426zcecA3s8eNvs3KlfcLGnibRcRO4C+AHZQCYA+wGXg+IgazZr3A1I0ml5ODoHq/GRGnAGcDl0j6LeBLwK8DJ1PawJ+rY32TNRM4BfhSRLwZeBEYMXR49mW/RjuWOFa/mmGbAZAd7loF3Dp6WYNuM6Bivxp6m2XBtZrSHye/BhwOrKxrURkHQZWyVCcingH+Hjg9InZFxMGIGAL+mtLIq42mF+iNiHuz6dso/QLdJekogOz/Z+pU32RV7FeTbLNhZwM/iYhd2XSjb7NhI/rVBNvs7cDjEfFsRAwA3wLOAo7MDhVB5aF4CucgqIKkwyXNH34MvBN4aPhDl3kP8FA96jsUEfFz4ElJJ2Sz3gZsY+QwIB8Avl2H8iZtrH41wzYrcz4jD5809DYrM6JfTbDNdgBnSporSfzqM/Z9SkPsQJ22l68aqoKk11PaC4DSIYe/i4jPSrqJ0u5qAE8AHyk7RtswJJ0MXA8cBjxG6SqNFuAW4BjgZ8B7I+IXdStyEsbo1xdojm12OKVfMK+PiD3ZvIU0/jar1K+G/5xJ+t+Ubr41CNwPXETpnEAn8Ops3gURMfW3IRuvLgeBmVnafGjIzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDXHzerO8sksn/yWbfC1wkNIQE1D68t+BuhRWgaQO4EBE/L9612JpcxBYU4mIPkrXmk+LUWElzSwbR2a0DqAfyB0EE6zPbFJ8aMianqRTJf0gGyjwjrLhF7olXStpU3afgtMkfSsbx/8zWZtjs/sYfCNrc5ukuTnW+5cq3a/io5LenY03f7+kOyW1SToW+G/ApdnY+m+R9FVJ55bV3Z/93yHpbkldlL4VPUPSn0vamA3A9pFavp/WfBwE1uwEfBE4NyJOBW4APlu2/EBEtANfpvTV/kuANwEfzA4zAZwA/FVEvBF4AfgfklonWO9hEdEeEZ+jNAb9mdmgd53AJyPiiew1r42IkyPi7gn6cQrw0Yh4A/BhYE9EnAacBlws6bjq3xqzEh8asmY3i9Iv9u+VhndhBqWRK4d1Zf8/CDw8PGSBpMeAJcDzwJMR8a9Zu68Dfwj88wTrvbns8WLg5myP4TDg8Un0476IGH7eO4FlZXsPRwBLJ7leMweBNT1R+gW/fIzlw2O6DJU9Hp4e/nyMHoclcqz3xbLHXwQ+HxFd2QniK8d4ziDZXrqkFkqhUWl9Av5nRNwxxnrMquJDQ9bs9gOvkbQcQFKrpN+och3HDD8f+ANKh3q2V7HeI/jV0MIfKJu/F5hfNv0EcGr2eBXQOsb67gD+e3Z4CklvKLvZjlnVHATW7IYoDfF7taQHgC3Af65yHdsp3YToEWABpZvcHKhivVcCt0raDDxXNv8fgPcMnyymNMb+W7P1LWfkXkC56ykNX/wTSQ8BX8F793YIPPqo2Tiyq3v+MSLeVOdSzArjPQIzs8R5j8DMLHHeIzAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS9z/Bxqb4z828mHUAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n", "import matplotlib.pyplot as plt\n", "\n", "data[\"Frequency\"]=data.Malfunction/data.Count\n", "data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression\n", "\n", "Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/aschmide/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", "Use an instance of a link class instead.\n", " import sys\n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Wed, 15 Apr 2020 Deviance: 3.0144
Time: 15:47: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: Wed, 15 Apr 2020 Deviance: 3.0144\n", "Time: 15:47:34 Pearson chi2: 5.00\n", "No. Iterations: 6 \n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n", "Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import statsmodels.api as sm\n", "\n", "data[\"Success\"]=data.Count-data.Malfunction\n", "data[\"Intercept\"]=1\n", "\n", "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", " family=sm.families.Binomial(sm.families.links.logit)).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/aschmide/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: DeprecationWarning: Calling Family(..) with a link class as argument is deprecated.\n", "Use an instance of a link class instead.\n", " \n" ] }, { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Generalized Linear Model Regression Results
Dep. Variable: Frequency No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Wed, 15 Apr 2020 Deviance: 18.086
Time: 15:47: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: Wed, 15 Apr 2020 Deviance: 18.086\n", "Time: 15:47:34 Pearson chi2: 30.0\n", "No. Iterations: 6 \n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n", "Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", "===============================================================================\n", "\"\"\"" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n", " family=sm.families.Binomial(sm.families.links.logit),\n", " var_weights=data['Count']).fit()\n", "\n", "logmodel.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\n", "The Goodness of fit (Deviance) indicated for this model is $G^2=18.086$ with 21 degrees of freedom (Df Residuals).\n", "\n", "**I have therefore managed to fully replicate the results of the Dalal *et al.* article**." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predicting failure probability\n", "The temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAZqUlEQVR4nO3df3TV9Z3n8efbECT8ECwqIwaF2SKOg/Ij4VdxbGirYNtB3GFE6tDaU0p3t7Q6jsyRM93qOHrOduOOzjqOI6usnXU1IMdJsYfT0Dpk7bqrBgqIwAbQppJoB8XyIzZICO/94/u96SUkuTc39+be++H1OCcn9/u9n+/3+3nnS15887nf+7nm7oiISPE7L98dEBGR7FCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEImWgm9laMztkZm/18LyZ2X81swNm9qaZTc9+N0VEJJV0rtCfARb08vxNwMT4awXwRP+7JSIifZUy0N39FeCjXprcDPyTR14DRpnZpdnqoIiIpGdQFvZxGXAwabk5Xvd+14ZmtoLoKp6ysrKKcePG9flgH51wPulwLLO+FhwH1VJgQqkDVEuhKj0PRpdl9hLmvn37PnT3i7t7LhuBnjZ3XwOsAaisrPStW7dmtJ/6+nqqqqqy2LP8US2FJ5Q6QLUUqv7UYma/6um5bNzl0gIkX2qXx+tERGQAZSPQNwJfje92mQ0cdfezhltERCS3Ug65mNnzQBVwkZk1A/cBpQDu/o/AJuCLwAHgt8DXc9VZERHpWcpAd/elKZ534NtZ65GIFI329naam5s5ceJEzo81cuRI9u7dm/PjDIR0ahkyZAjl5eWUlpamvd8BfVFURMLS3NzMiBEjGD9+PGa5vQfl+PHjjBgxIqfHGCipanF3Dh8+THNzMxMmTEh7v3rrv4hk7MSJE4wePTrnYX6uMTNGjx7d5798FOgi0i8K89zI5OeqQBcRCYTG0EWkqJWUlHDNNdd0LtfW1jJ+/Pj8dSiPFOgiUtTKysrYsWNHt8+5O+7OeeedG4MR50aVInLOaGpqYtKkSXz1q19l8uTJHDx4kOrqambMmMG1117Lfffd19n2oYce4sorr+S6665j6dKlPPzwwwBUVVWRmJrkww8/7Lzi7+joYNWqVZ37evLJJ4HfvZV/8eLFXHXVVdx+++1Ed3RDQ0MDn/nMZ5gyZQozZ87k+PHjLFiw4Iz/hK677jp27tzZ79p1hS4iWfHXL+1mz3vHsrrPq8dewH1//Ie9tmlra2Pq1KkATJgwgUceeYT9+/fzwx/+kNmzZ7N582b279/PG2+8gbuzcOFCXnnlFYYNG0ZNTQ07duzg1KlTTJ8+nYqKil6P9fTTTzNy5EgaGhr45JNPmDt3LjfeeCMA27dvZ/fu3YwdO5a5c+fy6quvMnPmTJYsWcK6deuYMWMGx44do6ysjGXLlvHMM8/w6KOPsm/fPk6cOMGUKVP6/fNSoItIUes65NLU1MQVV1zB7NmzAdi8eTObN29m2rRpALS2trJ//36OHz/OLbfcwtChQwFYuHBhymNt3ryZN998kw0bNgBw9OhR9u/fz+DBg5k5cybl5eUATJ06laamJkaOHMmll17KjBkzALjgggsAuOWWW5g7dy7V1dWsXbuWO+64Iys/CwW6iGRFqivpgTRs2LDOx+7O6tWr+da3vnVGm0cffbTH7QcNGsTp06cBzrgX3N157LHHmD9//hnt6+vrOf/88zuXS0pKOHXqVI/7Hzp0KDfccAM/+tGPWL9+Pdu2bUuvsBQ0hi4iQZs/fz5r166ltbUVgJaWFg4dOsT1119PbW0tbW1tHD9+nJdeeqlzm/Hjx3eGbOJqPLGvJ554gvb2dgD27dvHxx9/3OOxJ02axPvvv09DQwMQvUM0EfTLly/nu9/9LjNmzODCCy/MSq26QheRoN14443s3buXOXPmADB8+HCeffZZpk+fzpIlS5gyZQqXXHJJ57AIwD333MOtt97KmjVr+NKXvtS5fvny5TQ1NTF9+nTcnYsvvpja2toejz148GDWrVvHd77zHdra2igrK+NnP/sZABUVFVxwwQV8/etZnM8wcVvPQH9VVFR4prZs2ZLxtoVGtRSeUOpwz30te/bsyen+kx07diyn+7/vvvu8uro6p8dIOHbsmLe0tPjEiRO9o6Ojx3bd/XyBrd5DrmrIRURkgD333HPMmjWLhx56KKv3yGvIRUQEuP/++wfsWF/5ylfOepE2G3SFLiL94vEbaCS7Mvm5KtBFJGNDhgzh8OHDCvUs83g+9CFDhvRpOw25iEjGysvLaW5u5oMPPsj5sU6cONHngCtU6dSS+MSivlCgi0jGSktL+/SJOv1RX1/f+W7PYperWjTkIiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBCKtQDezBWbWaGYHzOzebp6/3My2mNl2M3vTzL6Y/a6KiEhvUga6mZUAjwM3AVcDS83s6i7Nvgesd/dpwG3AP2S7oyIi0rt0rtBnAgfc/R13PwnUADd3aePABfHjkcB72euiiIikw1J9WreZLQYWuPvyeHkZMMvdVya1uRTYDFwIDAO+4O7butnXCmAFwJgxYypqamoy6nRrayvDhw/PaNtCo1oKTyh1gGopVP2pZd68edvcvbLbJ9291y9gMfBU0vIy4O+7tLkb+Iv48RxgD3Beb/utqKjwTG3ZsiXjbQuNaik8odThrloKVX9qAbZ6D7mazpBLCzAuabk8XpfsG8D6+D+I/wsMAS5KY98iIpIl6QR6AzDRzCaY2WCiFz03dmnzLvB5ADP7A6JA/yCbHRURkd6lDHR3PwWsBOqAvUR3s+w2swfMbGHc7C+Ab5rZTuB54I74TwMRERkgg9Jp5O6bgE1d1n0/6fEeYG52uyYiIn2hd4qKiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEoi0At3MFphZo5kdMLN7e2hzq5ntMbPdZvZcdrspIiKpDErVwMxKgMeBG4BmoMHMNrr7nqQ2E4HVwFx3/42ZXZKrDouISPfSuUKfCRxw93fc/SRQA9zcpc03gcfd/TcA7n4ou90UEZFUzN17b2C2GFjg7svj5WXALHdfmdSmFtgHzAVKgPvd/Sfd7GsFsAJgzJgxFTU1NRl1urW1leHDh2e0baFRLYUnlDpAtRSq/tQyb968be5e2d1zKYdc0jQImAhUAeXAK2Z2jbsfSW7k7muANQCVlZVeVVWV0cHq6+vJdNtCo1oKTyh1gGopVLmqJZ0hlxZgXNJyebwuWTOw0d3b3f2XRFfrE7PTRRERSUc6gd4ATDSzCWY2GLgN2NilTS3R1TlmdhFwJfBOFvspIiIppAx0dz8FrATqgL3AenffbWYPmNnCuFkdcNjM9gBbgFXufjhXnRYRkbOlNYbu7puATV3WfT/psQN3x18iIpIHeqeoiEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBCKtQDezBWbWaGYHzOzeXtr9iZm5mVVmr4siIpKOlIFuZiXA48BNwNXAUjO7upt2I4A7gdez3UkREUktnSv0mcABd3/H3U8CNcDN3bT7G+AHwIks9k9ERNJk7t57A7PFwAJ3Xx4vLwNmufvKpDbTgb9y9z8xs3rgHnff2s2+VgArAMaMGVNRU1OTUadbW1sZPnx4RtsWGtVSeEKpA1RLoepPLfPmzdvm7t0Oaw/qV68AMzsP+FvgjlRt3X0NsAagsrLSq6qqMjpmfX09mW5baFRL4QmlDlAthSpXtaQz5NICjEtaLo/XJYwAJgP1ZtYEzAY26oVREZGBlU6gNwATzWyCmQ0GbgM2Jp5096PufpG7j3f38cBrwMLuhlxERCR3Uga6u58CVgJ1wF5gvbvvNrMHzGxhrjsoIiLpSWsM3d03AZu6rPt+D22r+t8tERHpK71TVEQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEP1+p6hIf9Rub6G6rpH3jrQxdlQZq+ZPYtG0y/LdLUmTzl9hUaBL3tRub2H1i7toa+8AoOVIG6tf3AWgUCgCOn+FR0MukjfVdY2dYZDQ1t5BdV1jnnokfaHzV3gU6JI37x1p69N6KSw6f4VHgS55M3ZUWZ/WS2HR+Ss8CnTJm1XzJ1FWWnLGurLSElbNn5SnHklf6PwVHr0oKnmTeOFMd0kUJ52/wqNAl7xaNO0yBUAR0/krLBpyEREJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJRFqBbmYLzKzRzA6Y2b3dPH+3me0xszfN7GUzuyL7XRURkd6kDHQzKwEeB24CrgaWmtnVXZptByrd/VpgA/Cfs91RERHpXTpX6DOBA+7+jrufBGqAm5MbuPsWd/9tvPgaUJ7dboqISCrm7r03MFsMLHD35fHyMmCWu6/sof3fA7929we7eW4FsAJgzJgxFTU1NRl1urW1leHDh2e0baFRLYUnlDpAtRSq/tQyb968be5e2d1zg/rVqy7M7M+ASuCz3T3v7muANQCVlZVeVVWV0XHq6+vJdNtCo1oKTyh1gGopVLmqJZ1AbwHGJS2Xx+vOYGZfAP4K+Ky7f5Kd7omISLrSGUNvACaa2QQzGwzcBmxMbmBm04AngYXufij73RQRkVRSBrq7nwJWAnXAXmC9u+82swfMbGHcrBoYDrxgZjvMbGMPuxMRkRxJawzd3TcBm7qs+37S4y9kuV8iGand3kJ1XSPvHWlj7KgyVs2fBHDWukXTLhuQY+fiOOn4Xu0unn/9IHdNbucbqzexdNY4Hlx0TV76IgMnqy+KiuRT7fYWVr+4i7b2DgBajrSx6oWdYNDe4Z3rVr+4CyCrYdvdsXNxnHR8r3YXz772budyh3vnskI9bHrrvwSjuq6xM1AT2k97Z5gntLV3UF3XmPNj5+I46Xj+9YN9Wi/hUKBLMN470paTtv3ZX7aPk46OHt5b0tN6CYcCXYIxdlRZTtr2Z3/ZPk46Ssz6tF7CoUCXYKyaP4my0pIz1pWeZ5SWnBlkZaUlnS+W5vLYuThOOpbOGten9RIOvSgqwUi8+JiPu1x6OnY+7nJJvPCZGDMvMdNdLucIBboEZdG0y7oN0YEI1p6OnQ8PLrqGBxddQ319PW/fXpXv7sgA0ZCLiEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiEHpNDKzBcDfASXAU+7+n7o8fz7wT0AFcBhY4u5N2e2qSLhqt7dQXdfIe0faGDuqjFXzJ/HC1nd59e2POtvM/Tef4k8rLz+rHXDWuq2/+ojnXz/IXZPb+cbqTSydNY4HF12T1nG729+iaZel3e/EsTvcKTHLybG727anPp5LUga6mZUAjwM3AM1Ag5ltdPc9Sc2+AfzG3T9tZrcBPwCW5KLDIqGp3d7C6hd30dbeAUDLkTbuWrfjrHavvv3RGQHfcqSNVRt2gkP7ae9cd/e6HZxO2q7DnWdfexfgjGDt7rirXtgJBu0dv9vf6hd3AZwVmN1tPxDH7m7bnvp4rklnyGUmcMDd33H3k0ANcHOXNjcDP4wfbwA+b2aWvW6KhKu6rrEznPqqvcM7wzzhdA9tn3/9YMrjtp/2zkBNaGvvoLqu8az9dbf9QBy7u2176uO5xty99wZmi4EF7r48Xl4GzHL3lUlt3orbNMfLb8dtPuyyrxXAinhxEpDpGbgI+DBlq+KgWgrPgNYx+Pc+XZGrfXf89iglQ0d2Lp/89YFtmR43edv+bp/hthcBH/a2bdc+FrD+/Bu7wt0v7u6JtMbQs8Xd1wBr+rsfM9vq7pVZ6FLeqZbCE0odENVy6uihYGoJ6bzkopZ0hlxagHFJy+Xxum7bmNkgYCTRi6MiIjJA0gn0BmCimU0ws8HAbcDGLm02Al+LHy8G/sVTjeWIiEhWpRxycfdTZrYSqCO6bXGtu+82sweAre6+EXga+B9mdgD4iCj0c6nfwzYFRLUUnlDqANVSqHJSS8oXRUVEpDjonaIiIoFQoIuIBKLgA93MhpjZG2a208x2m9lfx+snmNnrZnbAzNbFL9gWPDMrMbPtZvbjeLlY62gys11mtsPMtsbrPmVmPzWz/fH3C/Pdz3SY2Sgz22Bm/8/M9prZnGKsxcwmxecj8XXMzO4q0lr+PP59f8vMno9zoFh/V+6M69htZnfF63JyTgo+0IFPgM+5+xRgKrDAzGYTTS/wiLt/GvgN0fQDxeBOYG/ScrHWATDP3acm3U97L/Cyu08EXo6Xi8HfAT9x96uAKUTnp+hqcffG+HxMJZpX6bfAP1NktZjZZcB3gUp3n0x0M0ZiSpGi+l0xs8nAN4necT8F+LKZfZpcnRN3L5ovYCjwC2AW0busBsXr5wB1+e5fGv0vj0/e54AfA1aMdcR9bQIu6rKuEbg0fnwp0JjvfqZRx0jgl8Q3CBRzLV36fyPwajHWAlwGHAQ+RXQn3o+B+cX4uwL8KfB00vJ/BP4yV+ekGK7QE8MUO4BDwE+Bt4Ej7n4qbtJM9I+g0D1KdDITU16MpjjrAHBgs5lti6d0ABjj7u/Hj38NjMlP1/pkAvAB8N/jobCnzGwYxVlLstuA5+PHRVWLu7cADwPvAu8DR4FtFOfvylvAH5nZaDMbCnyR6E2YOTknRRHo7t7h0Z+R5UR/ulyV5y71mZl9GTjk7sUy10Qq17n7dOAm4Ntmdn3ykx5dehTDPbGDgOnAE+4+DfiYLn/+FlEtAMRjywuBF7o+Vwy1xOPJNxP9ZzsWGAYsyGunMuTue4mGijYDPwF2AB1d2mTtnBRFoCe4+xFgC9GfW6PiaQag++kICs1cYKGZNRHNWPk5orHbYqsD6LyKwt0PEY3TzgT+1cwuBYi/H8pfD9PWDDS7++vx8gaigC/GWhJuAn7h7v8aLxdbLV8AfunuH7h7O/Ai0e9Psf6uPO3uFe5+PdHY/z5ydE4KPtDN7GIzGxU/LiOal30vUbAvjpt9DfhRfnqYHndf7e7l7j6e6M/hf3H32ymyOgDMbJiZjUg8JhqvfYszp4Aoilrc/dfAQTObFK/6PLCHIqwlyVJ+N9wCxVfLu8BsMxsaT8OdOCdF97sCYGaXxN8vB/4t8Bw5OicF/05RM7uWaK71EqL/gNa7+wNm9vtEV7qfArYDf+bun+Svp+kzsyrgHnf/cjHWEff5n+PFQcBz7v6QmY0G1gOXA78CbnX3j3rYTcEws6nAU8Bg4B3g68T/1ii+WoYRBeLvu/vReF3RnZf49uQlwCmi34vlRGPmRfW7AmBmPyd6vawduNvdX87VOSn4QBcRkfQU/JCLiIikR4EuIhIIBbqISCAU6CIigVCgi4gEYkA/JFokXfFtXS/Hi79H9O66D+Llme5+Mi8d60Z8G+pJd/8/+e6LnNsU6FKQ3P0w0eyamNn9QKu7P5yv/pjZoKR5RLqqAlqBtAM9xf5EMqIhFykaZlZhZv8rnhCsLumt0/Vm9oiZbY3nM59hZi/Gc00/GLcZH893/j/jNhviyZJS7fdRi+Z7v9PM/jiej3u7mf3MzMaY2Xjg3wF/Hs9B/kdm9oyZLU7qd2v8vcrMfm5mG4E98aRz1WbWYGZvmtm3BvLnKeFRoEuxMOAxYLG7VwBrgYeSnj/p0bzs/0j0NupvA5OBO+LhG4BJwD+4+x8Ax4D/YGalKfY72N0r3f2/AP8bmB1P4lUD/KW7N8XHfMSjuch/nqKO6cCd7n4l0XzeR919BjAD+KaZTej7j0YkoiEXKRbnEwX0T6PpPSghmlo1YWP8fRewOzE1qZm9QzRd6RHgoLu/Grd7luhDFH6SYr/rkh6XA+viK/jBRPOo99Ub7p7Y7kbg2qSr+ZHAxAz3K6JAl6JhREE9p4fnE3N6nE56nFhO/DvvOs+Fp7Hfj5MePwb8rbtvjF8Ivb+HbU4R//VrZucRhX93+zPgO+5e18N+RPpEQy5SLD4BLjazOQBmVmpmf9jHfVye2B74CtEQSmMf9juS303Z+rWk9ceBEUnLTUQfAQfRvOSlPeyvDvj38bAPZnZlPLmWSEYU6FIsThNNnfoDM9tJ9EEBn+njPhqJPoxjL3Ah0YdanOzDfu8HXjCzbUQfh5bwEnBL4kVR4L8Bn433N4czr8qTPUU0LewvzOwt4En0V7P0g2ZblHNCfDfKjz360GGRIOkKXUQkELpCFxEJhK7QRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQC8f8Bc2QgqTmZRacAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n", "data_pred['Frequency'] = logmodel.predict(data_pred)\n", "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The predictiction is not correct. Let's see what went wrong" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAapUlEQVR4nO3df3xV9Z3n8dfHABJ+W9QMGDTMFnEckB/hZ3FsaKtg20HcYUTq0OpDSne3tDqOzEMe2xXr6OMx3bijs661ssraWVcCskyKfTANrUPWrh0VKCC/NoA2lUQ7KDZAbIAQPvvHOTe9xCT3JrmX3Pvl/Xw88uCec773nO8nh7xz8r3nfq+5OyIikv8u6u0OiIhIZijQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCkTLQzWy1mR0xsz0dbL/DzN4ys91m9gszm5D5boqISCrpXKE/D8ztZPuvgM+6+3jgb4BVGeiXiIh0UZ9UDdz9VTMr6WT7L5IWXweKe94tERHpqpSB3kV3A//U0UYzWwosBSgsLCwdNWpUtw5y9uxZLroojOF/1ZKbQqkllDpAtSQcOHDgQ3e/rN2N7p7yCygB9qRoMxvYDwxPZ5+lpaXeXVu2bOn2c3ONaslNodQSSh3uqiUB2OYd5GpGrtDN7DrgWeBmdz+aiX2KiEjX9PjvFzO7EtgALHb3Az3vkoiIdEfKK3QzWwOUAZeaWR2wEugL4O4/AB4EhgPfNzOAM+4+JVsdFhGR9qVzl8uiFNuXAEsy1iMRyRvNzc3U1dVx8uTJrB9r6NCh7N+/P+vHOR/SqaV///4UFxfTt2/ftPeb6btcROQCUldXx+DBgykpKSH+Cz1rTpw4weDBg7N6jPMlVS3uztGjR6mrq2P06NFp7zeMe4BEpFecPHmS4cOHZz3MLzRmxvDhw7v8l48CXUR6RGGeHd35virQRUQCoTF0EclrBQUFjB8/vnW5srKSkpKS3utQL1Kgi0heKywsZOfOne1uS7yDMpQpA1K5MKoUkQtGbW0tY8eO5atf/Srjxo3j8OHDlJeXM3XqVK677jpWrlzZ2vbRRx/l6quv5vrrr2fRokU89thjAJSVlbFt2zYAPvzww9Yr/paWFpYvX966r2eeeQaA6upqysrKWLBgAddccw133HFHYkoUtm7dymc+8xkmTJjAtGnTOHHiBHPnzj3nl9D111/Prl27ely7rtBFJCO++/Je9r13PKP7vHbkEFb+6R932qapqYmJEycCMHr0aB5//HEOHjzID3/4Q2bMmMHmzZs5ePAgb775Ju7OvHnzePXVVxk4cCAVFRXs3LmTM2fOMHnyZEpLSzs91nPPPcfQoUPZunUrp06dYtasWdx0000A7Nixg7179zJy5EhmzZrFa6+9xrRp01i4cCFr165l6tSpHD9+nMLCQhYvXszzzz/PE088wYEDBzh58iQTJvT8oyQU6CKS19oOudTW1nLVVVcxY8YMADZv3szmzZuZNGkSAI2NjRw8eJATJ05w6623MmDAAADmzZuX8libN2/mrbfeYv369QAcO3aMgwcP0q9fP6ZNm0ZxcTR7+MSJE6mtrWXo0KGMGDGCqVOnAjBkyBAAbr31VmbNmkV5eTmrV6/mzjvvzMj3QoEuIhmR6kr6fBo4cGDrY3dnxYoVfOMb3zinzRNPPNHh8/v06cPZs2cBzrkX3N158sknmTNnzjntq6urufjii1uXCwoKOHPmTIf7HzBgADfeeCM/+tGPWLduHdu3b0+vsBQ0hi4iQZszZw6rV6+msbERgPr6eo4cOcINN9xAZWUlTU1NnDhxgpdffrn1OSUlJa0hm7gaT+zr6aefprm5GYADBw7w8ccfd3jssWPH8v7777N161YgeodoIuiXLFnCt7/9baZOncoll1ySkVp1hS4iQbvpppvYv38/M2fOBGDQoEG88MILTJ48mYULFzJhwgQuv/zy1mERgPvvv5/bbruNVatW8aUvfal1/ZIlS6itrWXy5Mm4O5dddhmVlZUdHrtfv36sXbuWb33rWzQ1NVFYWMjPfvYzAEpLSxkyZAh33XVX5ortaKL0bH/pAy4iqiU3hVJLtuvYt29fVvef7Pjx41nd/8qVK728vDyrx0g4fvy419fX+5gxY7ylpaXDdu19f+nkAy405CIicp69+OKLTJ8+nUcffTSj98hryEVEBHjooYfO27G+8pWvfOJF2kzQFbqI9IjHb6CRzOrO91WBLiLd1r9/f44ePapQzzCP50Pv379/l56nIRcR6bbi4mLq6ur44IMPsn6skydPdjngclU6tSQ+sagrFOgi0m19+/bt0ifq9ER1dXXruz3zXbZq0ZCLiEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiJSBbmarzeyIme3pYLuZ2X81s0Nm9paZTc58N0VEJJV0rtCfB+Z2sv1mYEz8tRR4uufdEhGRrkoZ6O7+KvBRJ01uAf4hnnv9dWCYmY3IVAdFRCQ9mZjL5QrgcNJyXbzu/Qzs+xO++/JefrGviadr/iUbuz/vGhpUSy4KpZZQ6oCwahly9hRlZZnf73mdnMvMlhINy1BUVER1dXWX91FXd4qWlhYaGhoy3LveoVpyUyi1hFIHhFVLYWFLt/IvpY4+my75CygB9nSw7RlgUdJyDTAi1T71maIR1ZKbQqkllDrcVUsCWf5M0Y3AV+O7XWYAx9w9K8MtIiLSsZRDLma2BigDLjWzOmAl0BfA3X8AbAK+CBwCfgfcla3OiohIx1IGursvSrHdgW9mrEciItIteqeoiEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiLQC3czmmlmNmR0yswfa2X6lmW0xsx1m9paZfTHzXRURkc6kDHQzKwCeAm4GrgUWmdm1bZp9B1jn7pOA24HvZ7qjIiLSuXSu0KcBh9z9HXc/DVQAt7Rp48CQ+PFQ4L3MdVFERNJh7t55A7MFwFx3XxIvLwamu/uypDYjgM3AJcBA4Avuvr2dfS0FlgIUFRWVVlRUdKvTjY2NDBo0qFvPzTWqJTeFUksodYBqSZg9e/Z2d5/S7kZ37/QLWAA8m7S8GPhvbdrcB/xV/HgmsA+4qLP9lpaWendt2bKl28/NNaolN4VSSyh1uKuWBGCbd5Cr6Qy51AOjkpaL43XJ7gbWxb8g/gXoD1yaxr5FRCRD0gn0rcAYMxttZv2IXvTc2KbNu8DnAczsj4gC/YNMdlRERDqXMtDd/QywDKgC9hPdzbLXzB42s3lxs78Cvm5mu4A1wJ3xnwYiInKe9EmnkbtvAja1Wfdg0uN9wKzMdk1ERLpC7xQVEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQKQV6GY218xqzOyQmT3QQZvbzGyfme01sxcz200REUmlT6oGZlYAPAXcCNQBW81so7vvS2ozBlgBzHL335rZ5dnqsIiItC+dK/RpwCF3f8fdTwMVwC1t2nwdeMrdfwvg7kcy200REUnF3L3zBmYLgLnuviReXgxMd/dlSW0qgQPALKAAeMjdf9LOvpYCSwGKiopKKyoqutXpxsZGBg0a1K3n5hrVkptCqSWUOkC1JMyePXu7u09pb1vKIZc09QHGAGVAMfCqmY1394bkRu6+ClgFMGXKFC8rK+vWwaqrq+nuc3ONaslNodQSSh2gWtKRzpBLPTAqabk4XpesDtjo7s3u/iuiq/UxmemiiIikI51A3wqMMbPRZtYPuB3Y2KZNJdHVOWZ2KXA18E4G+ykiIimkDHR3PwMsA6qA/cA6d99rZg+b2by4WRVw1Mz2AVuA5e5+NFudFhGRT0prDN3dNwGb2qx7MOmxA/fFXyIi0gv0TlERkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJRFqBbmZzzazGzA6Z2QOdtPszM3Mzm5K5LoqISDpSBrqZFQBPATcD1wKLzOzadtoNBu4B3sh0J0VEJLV0rtCnAYfc/R13Pw1UALe00+5vgO8BJzPYPxERSZO5e+cNzBYAc919Sby8GJju7suS2kwG/qO7/5mZVQP3u/u2dva1FFgKUFRUVFpRUdGtTjc2NjJo0KBuPTfXqJbcFEotodQBqiVh9uzZ29293WHtPj3qFWBmFwF/B9yZqq27rwJWAUyZMsXLysq6dczq6mq6+9xco1pyUyi1hFIHqJZ0pDPkUg+MSloujtclDAbGAdVmVgvMADbqhVERkfMrnUDfCowxs9Fm1g+4HdiY2Ojux9z9UncvcfcS4HVgXntDLiIikj0pA93dzwDLgCpgP7DO3fea2cNmNi/bHRQRkfSkNYbu7puATW3WPdhB27Ked0tERLpK7xQVEQmEAl1EJBAKdBGRQPT4PnSRnqjcUU95VQ3vNTQxclghy+eMZf6kK3q7W5Imnb/cokCXXlO5o54VG3bT1NwCQH1DEys27AZQKOQBnb/coyEX6TXlVTWtYZDQ1NxCeVVNL/VIukLnL/co0KXXvNfQ1KX1klt0/nKPAl16zchhhV1aL7lF5y/3KNCl1yyfM5bCvgXnrCvsW8DyOWN7qUfSFTp/uUcvikqvSbxwprsk8pPOX+5RoEuvmj/pCgVAHtP5yy0achERCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKRVqCb2VwzqzGzQ2b2QDvb7zOzfWb2lpm9YmZXZb6rIiLSmZSBbmYFwFPAzcC1wCIzu7ZNsx3AFHe/DlgP/OdMd1RERDqXzhX6NOCQu7/j7qeBCuCW5AbuvsXdfxcvvg4UZ7abIiKSSjqBfgVwOGm5Ll7XkbuBf+pJp0REpOvM3TtvYLYAmOvuS+LlxcB0d1/WTtu/AJYBn3X3U+1sXwosBSgqKiqtqKjoVqcbGxsZNGhQt56ba1RLbgqlllDqANWSMHv27O3uPqXdje7e6RcwE6hKWl4BrGin3ReA/cDlqfbp7pSWlnp3bdmypdvPzTWqJTeFUksodbirlgRgm3eQq+kMuWwFxpjZaDPrB9wObExuYGaTgGeAee5+pDu/dUREpGf6pGrg7mfMbBlQBRQAq919r5k9TPSbYiNQDgwCXjIzgHfdfV4W+y3Srsod9ZRX1fBeQxMjhxWyfM5YgE+smz+ps5eBMnfsbBwnHd+p3M2aNw5z77hm7l6xiUXTR/HI/PG90hc5f1IGOoC7bwI2tVn3YNLjL2S4XyJdVrmjnhUbdtPU3AJAfUMTy1/aBQbNLd66bsWG3QAZDdv2jp2N46TjO5W7eeH1d1uXW9xblxXqYdM7RSUY5VU1rYGa0HzWW8M8oam5hfKqmqwfOxvHSceaNw53ab2EQ4EuwXivoSkrbXuyv0wfJx0tHdy51tF6CYcCXYIxclhhVtr2ZH+ZPk46CqLXsdJeL+FQoEswls8ZS2HfgnPW9b3I6FtwbpAV9i1ofbE0m8fOxnHSsWj6qC6tl3Ck9aKoSD5IvPjYG3e5dHTs3rjLJfHCZ2LMvMBMd7lcIBToEpT5k65oN0TPR7B2dOze8Mj88TwyfzzV1dW8fUdZb3dHzhMNuYiIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISCAW6iEggFOgiIoFQoIuIBEKBLiISiD7pNDKzucDfAwXAs+7+t222Xwz8A1AKHAUWunttZrsqEq7KHfWUV9XwXkMTI4cVsnzOWF7a9i6vvf1Ra5tZ/+ZT/PmUKz/RDvjEum2//og1bxzm3nHN3L1iE4umj+KR+ePTOu78SVd0uD6d5yeO3eJOgVmXjt1eLeket712F5qUgW5mBcBTwI1AHbDVzDa6+76kZncDv3X3T5vZ7cD3gIXZ6LBIaCp31LNiw26amlsAqG9o4t61Oz/R7rW3Pzon4Osbmli+fhc4NJ/11nX3rd3J2aTntbjzwuvvApwTrO0dd8WG3Wz79Uf87+31n1gPnBOa7T2/J8de/tIuMGhu+X0t6R63vXYXonSGXKYBh9z9HXc/DVQAt7Rpcwvww/jxeuDzZmaZ66ZIuMqralrDqauaW7w1zBPOdtB2zRuHUx63qbmFNW8cbnd9eVVNyuf35NjNZ701zLt63PbaXYjM3TtvYLYAmOvuS+LlxcB0d1+W1GZP3KYuXn47bvNhm30tBZbGi2OB7p6BS4EPU7bKD6olN523Wvr9wadLs7Xvlt8do2DA0Nbl0785tL0nx+3J8zPw3EuBDzt7bvIxclxP/n9d5e6XtbchrTH0THH3VcCqnu7HzLa5+5QMdKnXqZbcFEotZrbtzLEjeV8HhHNOIHu1pDPkUg+MSloujte128bM+gBDiV4cFRGR8ySdQN8KjDGz0WbWD7gd2NimzUbga/HjBcA/e6qxHBERyaiUQy7ufsbMlgFVRLctrnb3vWb2MLDN3TcCzwH/08wOAR8RhX429XjYJoeoltwUSi2h1AGqJaWUL4qKiEh+0DtFRUQCoUAXEQlEzge6mfU3szfNbJeZ7TWz78brR5vZG2Z2yMzWxi/Y5jwzKzCzHWb243g5X+uoNbPdZrbTzLbF6z5lZj81s4Pxv5f0dj/TYWbDzGy9mf0/M9tvZjPzsRYzGxufj8TXcTO7N09r+cv4532Pma2JcyBff1buievYa2b3xuuyck5yPtCBU8Dn3H0CMBGYa2YziKYXeNzdPw38lmj6gXxwD7A/aTlf6wCY7e4Tk+6nfQB4xd3HAK/Ey/ng74GfuPs1wASi85N3tbh7TXw+JhLNq/Q74B/Js1rM7Arg28AUdx9HdDNGYkqRvPpZMbNxwNeJ3nE/AfiymX2abJ0Td8+bL2AA8EtgOtG7rPrE62cCVb3dvzT6XxyfvM8BPwYsH+uI+1oLXNpmXQ0wIn48Aqjp7X6mUcdQ4FfENwjkcy1t+n8T8Fo+1gJcARwGPkV0J96PgTn5+LMC/DnwXNLyfwL+OlvnJB+u0BPDFDuBI8BPgbeBBnc/EzepI/pPkOueIDqZiSkvhpOfdQA4sNnMtsdTOgAUufv78ePfAEW907UuGQ18APyPeCjsWTMbSH7Wkux2YE38OK9qcfd64DHgXeB94Biwnfz8WdkD/ImZDTezAcAXid6EmZVzkheB7u4tHv0ZWUz0p8s1vdylLjOzLwNH3D1f5ppI5Xp3nwzcDHzTzG5I3ujRpUc+3BPbB5gMPO3uk4CPafPnbx7VAkA8tjwPeKnttnyoJR5PvoXol+1IYCAwt1c71U3uvp9oqGgz8BNgJ9DSpk3GzkleBHqCuzcAW4j+3BoWTzMA7U9HkGtmAfPMrJZoxsrPEY3d5lsdQOtVFO5+hGicdhrwr2Y2AiD+90jv9TBtdUCdu78RL68nCvh8rCXhZuCX7v6v8XK+1fIF4Ffu/oG7NwMbiH5+8vVn5Tl3L3X3G4jG/g+QpXOS84FuZpeZ2bD4cSHRvOz7iYJ9Qdzsa8CPeqeH6XH3Fe5e7O4lRH8O/7O730Ge1QFgZgPNbHDiMdF47R7OnQIiL2px998Ah81sbLzq88A+8rCWJIv4/XAL5F8t7wIzzGxAPA134pzk3c8KgJldHv97JfBvgRfJ0jnJ+XeKmtl1RHOtFxD9Alrn7g+b2R8SXel+CtgB/IW7n+q9nqbPzMqA+939y/lYR9znf4wX+wAvuvujZjYcWAdcCfwauM3dP+pgNznDzCYCzwL9gHeAu4j/r5F/tQwkCsQ/dPdj8bq8Oy/x7ckLgTNEPxdLiMbM8+pnBcDMfk70elkzcJ+7v5Ktc5LzgS4iIunJ+SEXERFJjwJdRCQQCnQRkUAo0EVEAqFAFxEJxHn9kGiRdMW3db0SL/4B0bvrPoiXp7n76V7pWDvi21BPu/sversvcmFToEtOcvejRLNrYmYPAY3u/lhv9cfM+iTNI9JWGdAIpB3oKfYn0i0acpG8YWalZvZ/4gnBqpLeOl1tZo+b2bZ4PvOpZrYhnmv6kbhNSTzf+f+K26yPJ0tKtd8nLJrv/R4z+9N4Pu4dZvYzMysysxLg3wF/Gc9B/idm9ryZLUjqd2P8b5mZ/dzMNgL74knnys1sq5m9ZWbfOJ/fTwmPAl3yhQFPAgvcvRRYDTyatP20R/Oy/4DobdTfBMYBd8bDNwBjge+7+x8Bx4H/YGZ9U+y3n7tPcff/AvxfYEY8iVcF8NfuXhsf83GP5iL/eYo6JgP3uPvVRPN5H3P3qcBU4OtmNrrr3xqRiIZcJF9cTBTQP42m96CAaGrVhI3xv7uBvYmpSc3sHaLpShuAw+7+WtzuBaIPUfhJiv2uTXpcDKyNr+D7Ec2j3lVvunvieTcB1yVdzQ8FxnRzvyIKdMkbRhTUMzvYnpjT42zS48Ry4v9523kuPI39fpz0+Eng79x9Y/xC6EMdPOcM8V+/ZnYRUfi3tz8DvuXuVR3sR6RLNOQi+eIUcJmZzQQws75m9sdd3MeViecDXyEaQqnpwn6H8vspW7+WtP4EMDhpuZboI+Agmpe8bwf7qwL+fTzsg5ldHU+uJdItCnTJF2eJpk79npntIvqggM90cR81RB/GsR+4hOhDLU53Yb8PAS+Z2Xaij0NLeBm4NfGiKPDfgc/G+5vJuVflyZ4lmhb2l2a2B3gG/dUsPaDZFuWCEN+N8mOPPnRYJEi6QhcRCYSu0EVEAqErdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQPx/tRIgu/jEeqIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1.2])\n", "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n", "plt.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "hideCode": false, "hidePrompt": false, "scrolled": true }, "source": [ "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing and plotting uncertainty" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXgUZZ4H8G9V9Zn7vrhBwcipgBkVZQxIMhAQdRGXWd31wFlFUWecR0RHDhUH3VURr1FnxnFhdZf1AlRAVBQQOTwIEA4NgXB0rs6d9Fn17h+dtAlnp5J00p3v53l4THeqqn+vneTbVe9b7ysJIQSIiIjaSO7qAoiIKDQxQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0CUqALF26FNnZ2RgyZAgOHTp0xm1UVcWiRYswceJEXHvttVi1alUwSiMiIp2CEiATJkzAypUr0atXr7Nus2bNGhQXF2PDhg34n//5HyxfvhzHjx8PRnlERKRDUAJkzJgxSE9PP+c2n3zyCWbMmAFZlpGQkICJEydi3bp1wSiPiIh06DZ9IDabDRkZGf7H6enpKCkp6cKKiIjoXLpNgBARUWgxdHUBzdLT03Hy5EmMGDECwOlnJIGqqmqApoXn9F6JiVGw2+u7uoxOE87tC+e2AWxfKJNlCfHxkbr27TYBkpubi1WrVmHSpEmorq7Gxo0bsXLlyjYfR9NE2AYIgLBuGxDe7QvntgFsX08UlEtYTz75JK6++mqUlJTgtttuw5QpUwAAs2fPxp49ewAA1113HXr37o1Jkybhpptuwpw5c9CnT59glEdERDpI4Tadu91eH7afFJKTo1FeXtfVZXSacG5fOLcNYPtCmSxLSEyM0rdvB9dCREQ9BAOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuDBAiItKFAUJERLowQIiISBcGCBER6cIAISIiXRggRESkCwOEiIh0YYAQEZEuhmC9UFFREebNm4fq6mrExcVh6dKl6N+/f6tt7HY7HnnkEdhsNni9XmRlZeGxxx6DwRC0MomIKEBBOwNZsGABZs2ahfXr12PWrFl4/PHHT9vmtddew6BBg7BmzRqsXr0a+/btw4YNG4JVIhERtUFQAsRut6OgoAB5eXkAgLy8PBQUFKCysrLVdpIkoaGhAZqmwe12w+PxIDU1NRglEhFRGwXl2pDNZkNqaioURQEAKIqClJQU2Gw2JCQk+Le75557cN9992HcuHFwOBz47W9/i9GjR7fptRITozq09u4mOTm6q0voVOHcvnBuG8D29UTdqnNh3bp1GDJkCP7xj3+goaEBs2fPxrp165CbmxvwMez2emia6MQqu05ycjTKy+u6uoxOE87tC+e2AWxfKJNlSfcH76BcwkpPT0dpaSlUVQUAqKqKsrIypKent9puxYoVmDZtGmRZRnR0NLKzs7F9+/ZglEhERG0UlABJTExEZmYm1q5dCwBYu3YtMjMzW12+AoDevXvj66+/BgC43W5s27YNF154YTBKJCKiNgraKKyFCxdixYoVyMnJwYoVK7Bo0SIAwOzZs7Fnzx4AwPz58/Hdd99h6tSpmD59Ovr374+bbropWCUSEVEbSEKIsOowYB9I6Arn9oVz2wC2L5R1+z4QIiIKPwwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEgXBggREenCACEiIl0YIEREpAsDhIiIdAk4QDZu3Aiv19uZtRARUQgJOEBefPFFjBs3DosXL8bu3bs7syYiIgoBAQfI6tWr8dZbb8FsNuO+++5DTk4OXnnlFRw/fjyg/YuKijBz5kzk5ORg5syZOHLkyBm3++STTzB16lTk5eVh6tSpqKioCLREIiIKIkkIIdq6kxAC27Ztw5///Gf89NNPuPTSSzFz5kzk5eVBls+cSbfeeituvPFGXHfddfjoo4/w3nvv4e233261zZ49e/Dwww/jH//4B5KTk1FXVweTyQSz2RxwbXZ7PTStzU0KCcnJ0Sgvr+vqMjpNOLcvnNsGsH2hTJYlJCZG6du3rTsUFxfj5ZdfxsKFC+FyuTB37lzMmDEDK1euxNy5c8+4j91uR0FBAfLy8gAAeXl5KCgoQGVlZavt3nrrLdx+++1ITk4GAERHR7cpPIiIKHgMgW64cuVKfPTRRzh69Ch+85vf4JlnnsGoUaP838/JycEVV1xxxn1tNhtSU1OhKAoAQFEUpKSkwGazISEhwb9dYWEhevfujd/+9rdobGzEtddei7vvvhuSJOltHxERdZKAA+Trr7/GbbfdhgkTJsBkMp32favViuXLl7erGFVVcfDgQfz973+H2+3GnXfeiYyMDEyfPj3gY+g9FQsVycnRXV1Cpwrn9oVz2wC2rycKOEBefPFFyLIMo9Hof87j8UAI4Q+UcePGnXHf9PR0lJaWQlVVKIoCVVVRVlaG9PT0VttlZGQgNzcXJpMJJpMJEyZMQH5+fpsChH0goSuc2xfObQPYvlAWlD6Q22+/Hfv27Wv13L59+3DHHXecd9/ExERkZmZi7dq1AIC1a9ciMzOz1eUrwNc3smXLFggh4PF48O233+Kiiy4KtEQiIgqigAPk4MGDGDlyZKvnRowYgQMHDgS0/8KFC7FixQrk5ORgxYoVWLRoEQBg9uzZ2LNnDwBgypQpSExMxOTJkzF9+nRccMEF+Kd/+qdASyQioiAK+BJWTEwMKioq/COkAKCiogJWqzWg/QcNGoRVq1ad9vwbb7zh/1qWZTzyyCN45JFHAi2LiIi6SMBnIJMmTcIf/vAHHDp0CA6HAwcPHsTDDz+M3/zmN51ZHxERdVMBB8iDDz6IQYMGYcaMGf4bBwcMGIDf//73nVkfERF1U22+E10IgaqqKsTHx3fL+zM4Cit0hXP7wrltANsXytozCivgPhAAqKurQ1FRERoaGlo9f/nll+t6cSIiCl0BB8j777+PxYsXIyIiAhaLxf+8JEn4/PPPO6U4IiLqvgIOkOeffx7Lli3D+PHjO7MeIiIKEQF3oquqetY7zYmIqOcJOEBmz56NV199FZqmdWY9REQUIgK+hPXWW2+hoqICb775JuLi4lp9b9OmTR1dFxERdXMBB8izzz7bmXUQEVGICThALrvsss6sg4iIQkzAfSButxvPP/88JkyYgNGjRwMAtmzZghUrVnRacURE1H0FHCBLlizBoUOH8B//8R/+O9AvvPBCvPPOO51WHBERdV8BX8LauHEjNmzYgIiICMiyL3dSU1NRWlraacUREVH3FfAZiNFohKqqrZ6rrKw8bUQWERH1DAEHSG5uLh5++GEcO3YMAFBWVobFixdjypQpnVYcERF1X22azr13796YNm0aamtrkZOTg5SUFMyZM6cz6yMiom6qzdO5A75LV5zOPfjCeUppILzbF85tA9i+UBaU6dybL101azmle58+fXS9OFFHkiSg7R+HiEivgAPk2muvhSRJaHnC0nwGsn///o6vjKiNNA1QNQGD0v3OjInCUcABcuDAgVaPy8vL8dJLL2HMmDEdXhSRXo1OD2KjTDwTIQqCgDvRT5WcnIxHH30Uzz33XEfWQ9Qubq8Gl5szRhMFg+4AAYDDhw/D4XB0VC1E7aYJgUanB91wfAdR2An4EtasWbNajbpyOBz4+eefOYyXuh2XV4XHy74Qos4WcIDMmDGj1WOr1YqLLroI/fv37+iaiNpFCMDh9iLaauzqUojCWsABcv3113dmHUQdyuHyItJigMxrWUSdJuAAWbZsWUDb3X///bqLIeoomibgdKuIMAf8I05EbRTwb9fRo0exYcMGDBs2DL169cLJkyexZ88eTJo0CWazuTNrJNKl0emB1WwAz0GIOkfAASKEwH/+538iJyfH/9yGDRuwbt06PP30051SHFF7eFUBt0eF2ah0dSlEYSngYbxff/01Jk6c2Oq57OxsfPXVVx1eFFFHaXR6OaSXqJMEHCD9+vXDypUrWz33zjvvoG/fvh1eFFFHcXtVuL28sZCoMwR8CevJJ5/EvffeizfffNO/EqHBYMDy5cs7sz6idhHCNyLLGGHq6lKIwk7AAXLxxRdj/fr12L17N8rKypCcnIxRo0bBaORYe+renC4VkRYBRea1LKKOpHsqk7Fjx8Lj8aCxsbEj6yHqcJoQcLq9XV0GUdgJ+Azk4MGDuPvuu2EymVBaWorJkydj586d+OCDD/DCCy90Zo1E7eZwehFp4dkyUUcK+Axk4cKFmDt3LtatWweDwZc7Y8eOxXfffRfQ/kVFRZg5cyZycnIwc+ZMHDly5KzbHj58GCNHjsTSpUsDLY/onLyagNOjdnUZRGEl4AD5+eefcd111wH4ZSGpiIgIuFyugPZfsGABZs2ahfXr12PWrFl4/PHHz7idqqpYsGDBaUOGidqLQ3qJOlbAAdKrVy/s3bu31XP5+fkBDeO12+0oKChAXl4eACAvLw8FBQWorKw8bdvXX38dv/71rzlJI3U4j1eFh0N6iTpMwH0g999/P373u9/h5ptvhsfjwV/+8he8++67eOKJJ867r81mQ2pqKhTFd0ewoihISUmBzWZDQkKCf7sDBw5gy5YtePvtt/HKK6/oaA50Lw4fKpKTo7u6hE7VnvapmoBXls65GqHVYkRcdNdMvcP3LrSFe/v0CDhArrnmGrz55pv43//9X4wdOxYnTpzA8uXLMWzYsA4pxOPx4E9/+hOefvppf9DoYbfXQ9PCcz3T5ORolJfXdXUZnaa97RMCqKpxnPP9r5UluB3WoF/K4nsX2sK5fbIs6f7gHVCAqKqKnJwcfPLJJ1i4cGGbXyQ9PR2lpaVQVRWKokBVVZSVlSE9Pd2/TXl5OYqLi3HXXXcBAGprayGEQH19fUBnOUSB8GoCHq8KE+fHImq3gAJEURQoigKXywWTqe139CYmJiIzMxNr167Fddddh7Vr1yIzM7PV5auMjAxs377d/3j58uVobGzEww8/3ObXIzoXp4cBQtQRAu5Ev/XWW/HAAw9gx44dKC4uxrFjx/z/ArFw4UKsWLECOTk5WLFiBRYtWgQAmD17Nvbs2aOveiIdnG4VAuF5mZMomCQhztXl6Lu0lJycjIsuusi3gySh5S6SJGH//v2dW2UbsA8kdHVEH0j5efpAmiVEm4N6FsL3LrSFc/va0wdy3jOQ5vU/Dhw4gAMHDiA7O9v/9YEDB7pVeBAFijcVErXfeQPk1BOUnTt3dloxRMHCy1hE7XfeAJFOGe94niteRCFB0wTXCSFqp/OOwlJVFd9++60/OE59DACXX35551VI1EkcTi8s0co5bzwkorM7b4AkJiZi/vz5/sdxcXGtHkuShM8//7xzqiPqRG6PBq/KdUKI9DpvgHzxxRfBqIMo6DQh4PKqiDAFPCEDEbXA3xwKefmFFVi3vRhOt4pIiwFjLkrBkL7xAe3rcHoYICGk+b2uqHEiKdaC3Ky+GDEoqavL6rF0r0hI1B3kF1Zg5WeHUN3ghtWioNbhweqtRThYXBXQ/l6VnemhouV7HWExoLrBjZWfHUJ+YUVXl9ZjMUAopK3bXgxFkWE2KpAkCSaDAkWRsXn3yYD2FwJc7jZEnPpem42+93rd9uKuLq3HYoBQSKuoccJkaP1jbFRkVNUFttAZ4LsnRONQrG7vTO+1ySCjosbZRRURA4RCWlKs5bRLUB5VQ3wb1vzgPSGh4UzvtdurISnW0kUVEQOEQlpuVl+oqgaXR4UQAm6vClXVcNXIjDYdx8Hlbru9U99rl8f3XudmnX9VVOocHH5CIa15BM667cVwOFXEWI3IvqRXwKOwmnm8GryqBkXmZ6ruquV7zVFY3QMDhELeiEFJGDEoqU2z8Z5KEwIuj4YIMwOkO2t+r6l74G8LUROHy9PVJRCFFAYIURPeE0LUNgwQoiZCAA4XO9OJAsUAIWrB5VbhVXkWQhQIBghRC5oQcHoYIESBYIAQncLh9HCNEKIAMECITuHrTOea6UTnwwAhOoNGpxfgmulE58QAIToDt0dlXwjReTBAKCzYa5zYuOsY6h0dczOgAFDf6IbgWQjRWXEqEwoLq7cWYXO+DfHRZvxr7kVIibe2+5heVaDe4UW01dgBFRKFH56BUFi4cng6TEbfOiCvfbQXPx+v6ZDjNjo9vDud6CwYIBQWBveJw8OzLkV0hBFOt4q3Pt2P7QWl7T6uEEBto4sLThGdAQOEwkb/tBjMuWE40hMjoAngoy1FWPPNEag6ZudtyesVqGv0gKOyiFpjgFBYiYsy465pQ5HZz7ceyLa9JfjHpwfgcLVv3XOHy4tGN+8NIWqJAUJhx2xU8NtrB+PqplUJfz5Rg1c+2IvSysZ2Hbe+wQMP58ki8mOAUFiSZQm5WX1x0zUXwKBIsNc68epHe7G3qFL3MTUhUNvA/hCiZgwQCmujLkzCXdOGIjbSBLdHw39/dgjrdxTrWrUQADxegcpaJ0dmEYEBQj1A7+QozLlhOAakxwAAvvrxJP7+6X7dNx16VYHqOhfcHvaJUM/GAKEeIcpqxO1TMjFueDoAoPBELV56fw+OltTpOp4mBKrr3XBzuhPqwYIWIEVFRZg5cyZycnIwc+ZMHDly5LRtXn75ZUyZMgVTp07FDTfcgM2bNwerPOoBFFnC5Mv74Z8nXgiTUUZtgxtvrCnA17tP6urX0IRATYML3nYOEyYKVUELkAULFmDWrFlYv349Zs2ahccff/y0bUaMGIH/+7//w5o1a7BkyRI8+OCDcDqdwSqReojhAxMx5/rhSI23QhMC67YX47/WHdR1SUvVBKrrnAwR6pGCEiB2ux0FBQXIy8sDAOTl5aGgoACVla1HxFx11VWwWn1zGA0ZMgRCCFRXVwejROphkuOsuPv6YRgzJBkAcPBYNV56Lx+FJ9s+BYpXFaiqdbJPhHqcoASIzWZDamoqFEUBACiKgpSUFNhstrPu8+GHH6Jv375IS0sLRonUA5kMCm4YPwg3ZV/gu6TV6MHf1u7H+h3FbV4XXdUEqupdaHBxHRHqObrlbLw7duzAsmXL8Le//a3N+yYmRnVCRd1HcnJ0V5fQqdrTPlUT8MpSm5ejzb4sEsMuTMZfV+/DUVstvvrxJIpsdbh92lCkJUa2uQ7JqCAu2gyD0vrzGd+70Bbu7dMjKAGSnp6O0tJSqKoKRVGgqirKysqQnp5+2rY//PAD/vjHP+KVV17BwIED2/xadnu97jH+3V1ycjTKy/WNGgoF7W2fEEBVjUPX+28AcOeUi7Bx13F8/eNJFJfW4cm/bUduVl/8amgaZElq0/Eq7DLio02QJV+I8L0LbeHcPlmWdH/wDsolrMTERGRmZmLt2rUAgLVr1yIzMxMJCQmttsvPz8eDDz6IF198EUOHDg1GaUR+iiwj57K+uHPqxYiLMsGrCqz95ij+9vF+VNW1bTCHx6uhsoYjtCi8SUIEZ16GwsJCzJs3D7W1tYiJicHSpUsxcOBAzJ49G3PnzsXw4cNx44034sSJE0hNTfXv98wzz2DIkCEBvw7PQEJXR5yBlOs8AzmV0+3Fx98cxXeHygEAJqOM3Ky+uCwztU1nIwZFQkK0GampsXzvQlg4t689ZyBBC5BgYYCEru4UIM027jqGr3486Z8SPjXeiiuGpWH3zxWoqnMhPtqMq0ZmYEjf+LMew6BI6NMrDo31LrTtQphPfmEF1m0vRkWNE0mxFuRm9cWIQUk6W9SxVm85jA07j8PpUWExKpg0tjemjWv7pefuLpx/97r9JSyiUHSwuAo//FSO2CgTrGbfCMLSKgc+2FyEkioHzCYFtQ4PVm8twsHiqrMex6sK1NS7Ya91oMHlbdP6JPmFFVj52SFUN7gRYTGgusGNlZ8dQn5hRbvb116rtxzG6m+OwOVRYZABl0fF6m+OYPWWw11dGgUJA4ToLDbvPglFkWExGRAfbUFCjNn/vUanFxU1TkAAiiJj8+6T5z2e1ytQ1+CGvdaJOkfz1PDnDpN124uhKDLMRgWSJMFsVKAoMtZtL25v89ptw87jkCBBkSVIkuz7LyRs2Hm8q0ujIOmWw3iJuoOqOhcs5l9+RSwmAyS4/H/yvapARY0TVrMCTxtuItQ0gQaHB41ODwyKBKvZiAjzmX8VK2qciLC0/p7JIPvCq4s53V4ocuuLcrLke556Bp6BEJ1FfLT5tAWkDIoEoyIhKdYCY9N9Hg6XirpGD74tKGlT/4sQvunhaxvcqGlwQ5zhbCQp1nLa1PFur4akWIuOFnUsi8mAU5urCd/z1DMwQIjO4qqRGVBVDW6vCiEE3F4VJqMCk8kASEBirBmRVgMk+P5wrt5yBC9/sAeHT9a2+bUcLi/s1U40nvLpPTerL1RVg8vjq8HlUaGqGnKz+nZQK/WbNLY3BARUTUAIzfdfCEwa27urS6MgURYuXLiwq4voSA6Hu813IoeKyEgzGhvdXV1Gp+mI9jW6vB32/ifFWpEUa0GpvRH1jR7ERZrwm1/1w8X9433POTxIirXg2rF9YDEZYGt67vtD5SixN6JXUqT/8pPVaoLjPJM1agJwuVXfJSDZ17eQlhCJ1HgrjpfVo6bBjYRoM264emC3GIU1pG88IASOltTDrQpYjAom/6pvWI7CCuffPUmSEBFh0rcvh/GGjnAeSgh0z2G8bXG8rB5rtx1BcWk9AN/08VkXp+KaS3uhT0YcKisb2nQ8WZZgNRsQYVagyN37YgF/NkMXh/ESdQO9U6Lwu2lDMTP7AsRFmaBqAt/sLcF/vPMjPtlaBFcbZ+tt7my317hQ2+gOaNQWUTCxt4uoA0mShJEXJOHi/gnYtq8Em344AadbxerNh/H5rmP49agMXJaZCqMh8M9umhBodHrhcHlhUGRYzAaYDQoMip7bEok6DgOEqBMYDTKuHpmBMUNS8PXuE9i2rxQNDg8+3nYUm/Nt+PWoDIy5KOW0GXvPxTdqS4PH60a9BCiKBKvJ4L83hHFCwcY+kBASztdhgdDvAzkXyajggy9+wncHy/13osdEmnD1yHSMuSgFJoOi/9gSIEsSrBYDIswGKDqmtG8v/myGrvb0gfAMhCgI4qMtmH7VQIwflYEvfziJ7w+Wo7bBjbXfHMWXP5zElcPSkHVxKqxnuaHwXIQAVCFQ3+hBo9N3mcto8N0ZbjLKMCpy2I5MpK7FACEKovhoC264eiCuuSQDX/14Et8dLEeDw4MNO32TNo69KAVXDE9DXJT5/Ac7A00TcGuqf3ldSfKtvGi1GGCQJV7qog7FACHqAs1nJNmX9sbWPTZs318Kl0fFlj02fLPXhmEDE3Hl8DT0SWnfKnhC+CY5dHlU/6Uuo0GGyajAoMgwtaEzn+hUDBCiLhTTdHPiry/phe0Fpdi2twR1Dg/yC+3IL7Sjd3IkLh+WhuEDE9vU4X4mzZe6VLcKp1uFBECSJZgMMkwmBQZZhkGR/PNb8bIXnQ870UNIOHfkAeHdiZ6QEBnQjYReVUN+oR1b99hgszf6n4+wGDBmSDLGZqYiMaZz5sGS4BuGLMuAQZabvvY9lmUJiiRDUQBZliFLrQOGP5uhi53oRGHCoMi4dHAyLrkwCUdK6rBtXwkKiirR6PTi6902fL3bhkG9YjBmSAou7p/QpvtJzkcAEEJAUwGveuabHiXJFzLNHfW+MxYZqqpBknjW0tMwQIi6IUmSMCA9BgPSY1Db4MbOA2XYdaAMNQ1uFJ6oReGJWlhMCkZekIRLByehd3IUpDYstauXEL6QObWjXjIqqK91NY36UmAwSF0ynJiCiwFC1M3FRJowYXRvXHNJL/x0vBo7D5ThwNFqON0qtheUYntBKZJiLRh5QRJGXpCIpFhrUOvzhcovnfWAx3fJSwaMigJZkZrOWiQoku+MRZZ5thIOGCBEIUKWJQzpG48hfeNR1+jG7p/tvpl/KxtRUePE598dx+ffHUevpEgMH5SI4QMTEB/dNeuGaJqApgEeb+vp6c92CawtnfeSBGgaICAgSxKaT7wYSMHHACEKQdERJowbkY5xI9Jhszfgx58qsLvQjtoGN05UNOBERQPWbS9Gr6RIDB2QgIv7JyAlPrhnJmdy1ktgkgRZ8nXey0pT531TJ76m+dYcaf7a7VWhqr+khSz5+o4U5ZdAkpumepElXkbrTAwQohCXnhiJ9MRI5GT1xdGSOuQX2rGvqBL1Do8/TDbsPIakWAsy+8Xjon7x6JsafdpytF2lOVQ0NHXen3vZlNP498MvHf/+UGoaUdYcLr7v+YJFkZuDimcvejFAiMKE3KLjfeoV/XGkpA57i+zYf6QKNQ1uVNQ4sTnfhs35NljNCi7sHYfBfeJwQe9YxOhcUKi78oeSBnjROlyAX4YsSxIgKxKMsgylafoXpensxxdCgITm4Gk6LoQ/cHwZ3D2CuCswQCisyJIEIZ3942RP+aQpyxIGZsRgYIYvTE5UNGD/kSocKK6Czd4Ih0v136wIAGkJEbigdywGZcSgf3oMzEb9kzuGguYhyxCAqgl4oAEu3/f8odGUCxIAVZJQVe1oChDfASQZSIyxoJucyHUJBgiFDUkCEmLMZ72RUAAQGqAKzXddXRXwagKaprX6w9DyE2Y4kCQJvZOj0Ds5CteO7YPqehcOFlfjp+PV+PlEDdweDSWVjSipbMSWfBtkSULvlEgMzIjFgPRo9E2JhtkU3oHSUvNZRsu1u1QBeNXWPxSy1oOTowkDhMKKLEmQz7XQkgK0XIjz1BE8qiagCd8/VW3+p0GF8G/U/GdEaC2+FqETOnFRZmRdnIqsi1PhVTUcK6vHT8drUHiiBsfL61Ku2nEAABJZSURBVKEJgeLSehSX1mPTD77LNOlJkeiXGo1+adHomxqN2MjwuuRF+jBAqEc79Y++IktQmq9dNP12nO3+PP8nVfiCx6sKf5AIIeBweuHthtOqtGRQZH+/Ccb2gdPtRZGtDodP1qDIVgebvQGaAE6UN+BEeQO+2VsCwHdvSp/kKPRJiUKvlEhYI/XNHkyhjQFCdB7nOrNo7mA1yBJOnVUkwmKAy6NB1QSsJgOMBhleVevWZyoWkwGZ/eKR2S8eAOB0e3G0pM73r7QOx8sa4FE11Da4sa+hEvuOVDbtuR9JsRZkJEUiIzES6UkRSE+MRJTV2HWNoU7HACHqJBIkWJo6oxNiLdA8HnhVDS6vBpdLhUfVuv2lL4vJ4L95EQBUTUOJvRHFZfU4XlaPY2X1sNc4IQBU1DhRUeP0d8wDQLTViLTECKQmRCA13orU+Agkx1vDvpO+p2CAEAWJEIAiy4gwyYg0G3657AXhu3Nb+O7gVr0aPKoGrRuGiyLL6JUchV7JUcBQ33NOtxd1Lg0Hiipgq2jEiYoGVNT4RizVOTyoO16Dn47XtDpOXJQJyXFWJMdZkRRnQVKsFUmxFsREmiAHYU4v6hgMEKIuIISvw99kOP2PZfP9BqqmQdV8/3W5VLi9vlDpbiwmAzLSIpEc/UvHuserobSqESV23+iu0qpGlFQ60ODw3SVYXe9Gdb37tGAxKjISYsxIiLEgMcaC+BgzEqLNiI+xID7K3KGzD1P7MUCIupnmjFBkGb41pGREmA3wqgIur+rrnO/mfSlGg+wfOtxSg9ODsioHyqocqKh2oLzGgfJqJ6rrXBAAPKqG0ioHSqscZzxulNWIuCgT4qLNiIs0IzbKhNgoM2IjTYiNNCHKaoTck2/MCDIGCFEI8F3+khBhMiDCpMCrCnhUDR6PBrdXhaahW56dnCrSYsSAdKNv1FcLHq+Gyjon7DVO2Gt9/62qc6Gy1oXqehfUptFs9Q4P6h0eHC8/8+JckuTrd4mONCHaakJ0hLHpny9coiOMiLIaEWk1wmSQgzIFfjhjgBCFHN/U6AZFhrXpqpGq+e5Xcasa3B6t6Qyl+/WhnI3RICM1PgKp8RGnfU/TBGob3aiqc6Gqzhco1fVu1NS7UNPgRk29u2kaeV/Q1jZ6UNvoAXDuFSCNiowIiwGRViMiLQZEWoywWgyIMBsQ0eK/VrMBqizD7fLCbFLYR9MCA4QoDPjmcFJgMiqQrL6zEbVpFluP13em4u2mHfPnI8sS4qLMiIsyY0D6mbdxur2obfCgtsGN2kY36hrdqG3woM7hRl2jB/WNvq/dHs2/j0fVfAHU4A64FgmA2aTAYlIQaTHi+vEDcemFye1sYegKWoAUFRVh3rx5qK6uRlxcHJYuXYr+/fu32kZVVTz55JPYvHkzJEnCXXfdhRkzZgSrRKKwIIRvCHHzvSlmg4J9JXZ8+f0J1Da4kRhnxeVD0zAgPQb/9+VP+PEnOzyqBlXVMLR/PEZemIzNu0+iqs6F+GgzrhqZgSF943GwuOq050+U12NLfglcXhVmg4JxI9KQPbrPGes60/5nO27zsOFAj/HNHpuvDo8Ks9FXx9Qr+5+2v9urYs/PdmzbV4KaBjcizAb0TY2C06PiaEkdXG4VsizBbFTg1QScLi9a3gsqADjdKpxuFdX1bvz94/3YmHIMuVl9MWJQUge8e6FFEiI4n0duvfVW3Hjjjbjuuuvw0Ucf4b333sPbb7/dapsPP/wQa9aswRtvvIHq6mpMnz4d//3f/43evXsH/Dp2e/1Z50IKdcnJ0Sgvr+vqMjpNOLevK9uWX1iBlZ8dgqLIMBlk32guVUNKnBWFtrqmhZ2a1tJQZCiSQFTTcFqnW4XHq2LEgER8/3MFJMm3qqDbq6Km3gmnW/jW3pABtWlqlwmX9jotRA4WV2H11iIoigyjIvsDa/TgZHx3qPy056ddOeC0EDnbMfqlRGH34UpI8E27oom21eFwegFJgtWsnFbD4D5xcHlUmCwmlJTVweHyovBEDXYcKIXJoCAl3oraRg9UVcNvrx0ckiEiyxISE6POv+EZBOUMxG63o6CgAH//+98BAHl5eXjiiSdQWVmJhIQE/3affPIJZsyYAVmWkZCQgIkTJ2LdunW48847A36tcB+BwfaFrq5q2/aCUqQkRMBk+OXmPbdXRVmVA9ERxlaTkauqBiEB0VYTJFlCTATg1QQOl9ShT2o0LEbFPw16RY0TAPyz0WqagFfVUFLlQITVCAjfPS5CAwpP1qBvWoxvGG7TH3i3V8VhWx16pUTBICtA0ySWHlVFwdEqDBuY2Kod+4oqz9gOW6UDqfHWVv9/NU3g4LEa5GT1O+8x7DW+EV+JLZYCdntV7CuqxLCBiTAZFcTFRcDcNIJ4988V6JMajSirCVazgkirF26vb3nhUSF4Oas9P5dBCRCbzYbU1FQoiu9NUxQFKSkpsNlsrQLEZrMhIyPD/zg9PR0lJSVteq34+MiOKbqb0vtJIVSEc/u6qm3zb/9Vl7xuSw8PaP8n8wWD2v/HuT3HSEuM7LA6wgXvyiEiIl2CEiDp6ekoLS2FqvqG2qmqirKyMqSnp5+23cmTJ/2PbTYb0tLSglEiERG1UVACJDExEZmZmVi7di0AYO3atcjMzGx1+QoAcnNzsWrVKmiahsrKSmzcuBE5OTnBKJGIiNooaKOwCgsLMW/ePNTW1iImJgZLly7FwIEDMXv2bMydOxfDhw+HqqpYvHgxtm7dCgCYPXs2Zs6cGYzyiIiojYIWIEREFF7YiU5ERLowQIiISBcGCBER6cIAISIiXUJ2Nt577rkHx48fhyzLiIiIwJ/+9CdkZmYGNGljqHjppZewfPlyrFmzBoMHD8aPP/6Ixx9/HC6XC7169cKzzz6LxMTE8x+oG8rOzobJZILZbAYAPPTQQ7jqqqvCoo0ulwtLlizBtm3bYDabMWrUKDzxxBNh8bN5/PhxzJkzx/+4rq4O9fX12LFjR1i0DwC+/PJLLFu2rGk6fIF7770XkyZNCpv2bdq0CcuWLYPX60VsbCyefvpp9OnTR1/7RIiqra31f/3ZZ5+J6dOnCyGEuOWWW8SHH34ohBDiww8/FLfcckuX1Ndee/fuFXfccYe45pprxMGDB4WqqmLixIli586dQgghXn75ZTFv3rwurlK/5na1FC5tfOKJJ8RTTz0lNE0TQghRXl4uhAifn82WnnzySbFo0SIhRHi0T9M0MWbMGP/P5v79+8WoUaOEqqph0b7q6mpx2WWXicOHDwshfO24/fbbhRD63r+QDZCWPvjgA3H99deLiooKMXr0aOH1eoUQQni9XjF69Ghht9u7uMK2cblc4qabbhLHjh3z/6HdvXu3mDJlin8bu90uRo0a1YVVts+ZAiQc2lhfXy9Gjx4t6uvrWz0fLj+bLblcLpGVlSX27t0bNu3TNE1cdtllYteuXUIIIXbs2CEmTZoUNu3bvXu3mDx5sv9xVVWVGDx4sO72hewlLAB49NFHsXXrVggh8OabbwY8aWN3t2zZMkybNq3VNPanTjSZkJAATdP8p5uh6KGHHoIQAqNHj8bvf//7sGjjsWPHEBcXh5deegnbt29HZGQk7r//flgslrD42Wzpiy++QGpqKoYOHYq9e/eGRfskScILL7yAe+65BxEREWhoaMDrr78eNn9bBgwYgIqKCuTn52PEiBFYs2YNgMAnvD1VSHeiP/XUU9i0aRMefPBBPPPMM11dTof44YcfsHfvXsyaNaurS+lUK1euxOrVq/Hee+9BCIHFixd3dUkdQlVVHDt2DBdffDHef/99PPTQQ7jvvvvQ2NjY1aV1uPfeew833nhjV5fRobxeL/7yl7/glVdewZdffolXX30VDzzwQNi8f9HR0Xj++efx9NNP44YbboDdbkdMTIzu9oV0gDSbPn06tm/fjrS0tIAmbezOdu7cicLCQkyYMAHZ2dkoKSnBHXfcgaNHj7aaaLKyshKyLIfMJ/NTNb8nJpMJs2bNwvfff3/aZJqh2Mb09HQYDAbk5eUBAEaOHIn4+HhYLJaQ/9lsqbS0FDt37sTUqVMBBD5hane3f/9+lJWVYfTo0QCA0aNHw2q1wmw2h0X7AOCKK67AO++8g/fffx//8i//AqfTiV69eulqX0gGSENDA2w2m//xF198gdjY2IAnbezO7rrrLmzZsgVffPEFvvjiC6SlpeGvf/0r7rzzTjidTuzatQsA8O677yI3N7eLq9WnsbERdXW+1fmEEPjkk0+QmZmJYcOGhXwbExISkJWV5Z/PraioCHa7Hf379w/5n82WPvjgA4wfPx7x8b5VA8Phdw8A0tLSUFJSgsOHDwPwzeFnt9vRr1+/sGgfAJSXlwMANE3Dc889h5tvvhm9evXS1b6QnAuroqIC99xzDxwOB2RZRmxsLB5++GEMHTr0rJM2hqrs7Gy89tprGDx4ML7//nssWLCg1RDXpKTQW0Lz2LFjuO+++6CqKjRNw6BBg/DYY48hJSUlLNp47NgxzJ8/H9XV1TAYDHjggQcwfvz4sPrZzMnJwaOPPoqrr77a/1y4tG/16tV44403IEm+lfrmzp2LiRMnhk37Hn30UXz//ffweDy48sorMX/+fJjNZl3tC8kAISKirheSl7CIiKjrMUCIiEgXBggREenCACEiIl0YIEREpAsDhIiIdAnpubCITnXJJZf4v3Y4HDCZTP75fRYtWoRp06Z1VWm6ZWdn48knn8QVV1zR1aUQtcIAobDyww8/+L8OhT+8Xq8XBkPn/hoG4zWoZ+IlLOoRNE3D66+/jokTJyIrKwv3338/qqurAfgWSRoyZAjee+89jB8/HmPHjsU777yD/Px8TJ06FWPGjGk12eP777+Pm2++GYsXL8bo0aORm5uLbdu2+b9fV1eH+fPnY9y4cbjqqqvw/PPP++cYat53yZIlyMrKwvLly1FcXIxbb70VWVlZyMrKwh/+8AfU1tYCAP74xz/i5MmT+Pd//3dccskleOONN7B9+/ZWd4ADvrD85ptvAADLly/H3Llz8dBDD+HSSy/FBx98cM6aiPRigFCP8F//9V/YuHEjVqxYgc2bNyM2Nva0GYB3796NDRs24Pnnn8eSJUvw2muv4a233sLHH3+MTz/9FDt27PBvm5+fj759++Lbb7/F3Llzce+99/oDad68eTAYDNiwYQM+/PBDbN26FatWrWq1b58+fbB161bcfffdEELgd7/7HTZv3oxPP/0UJSUlWL58OQDg2WefRUZGBl577TX88MMPmD17dkDt/fzzz5Gbm4tdu3Zh6tSp562JSA8GCPUI7777Lh588EGkpaXBZDLh3nvvxfr16+H1ev3bzJkzB2azGePGjUNERATy8vKQmJiI1NRUjBkzBgUFBf5tExIS8K//+q8wGo2YPHkyBgwYgE2bNqGiogJfffUV5s+fj4iICCQmJuLf/u3f8PHHH/v3TUlJwS233AKDwQCLxYJ+/frhyiuvhMlkQkJCAm677Tbs3LmzXe0dNWoUJk6cCFmWUV9ff96aiPTghVHqEU6ePIk5c+ZAln/5zCTLMux2u/9xy7XXzWbzaY9brpmQmprqn2wPADIyMlBWVoaTJ0/C6/Vi3Lhx/u9pmtZqWuy0tLRWtVVUVOCpp57Crl270NDQACEEYmJi2tXelq8RSE1EejBAqEdIS0vDkiVL/Os8tHT8+PE2H6+0tBRCCH+I2Gw2ZGdn+89wvv3227N2XLcMHgB47rnnIEkS1qxZg7i4OGzcuPGcC2xZrVY4nU7/Y1VVUVlZedbXCKQmIj14CYt6hH/+53/GCy+8gBMnTgDwLVa1ceNG3cerrKzE22+/DY/Hg08//RSFhYUYP348UlJScOWVV+LPf/4z6uvroWkaiouLW/WfnKqhoQERERGIjo5GaWkp3nzzzVbfT0pKwrFjx/yPBwwYAJfLhU2bNsHj8eDVV1+F2+0+6/H11EQUCAYI9Qi33norsrOzcfvtt+OSSy7BTTfdhPz8fN3HGzFiBI4ePYpf/epXeOGFF/Diiy/6F1d65pln4PF4MHnyZIwdOxZz5871L+JzJvfeey8KCgowZswY3HXXXZg0aVKr799111149dVXMWbMGPz1r39FdHQ0FixYgMceewxXX301rFbraZfFTtXWmogCwfVAiNro/fffx6pVq/DOO+90dSlEXYpnIEREpAsDhIiIdOElLCIi0oVnIEREpAsDhIiIdGGAEBGRLgwQIiLShQFCRES6MECIiEiX/wfZJXJyJQSzcgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sns.set(color_codes=True)\n", "plt.xlim(30,90)\n", "plt.ylim(0,1)\n", "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"." ] } ], "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.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }