"
+ ]
+ },
+ "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": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "
Generalized Linear Model Regression Results
\n",
+ "
\n",
+ "
Dep. Variable:
Frequency
No. Observations:
23
\n",
+ "
\n",
+ "
\n",
+ "
Model:
GLM
Df Residuals:
21
\n",
+ "
\n",
+ "
\n",
+ "
Model Family:
Binomial
Df Model:
1
\n",
+ "
\n",
+ "
\n",
+ "
Link Function:
logit
Scale:
1.0000
\n",
+ "
\n",
+ "
\n",
+ "
Method:
IRLS
Log-Likelihood:
-3.9210
\n",
+ "
\n",
+ "
\n",
+ "
Date:
Thu, 22 Aug 2024
Deviance:
3.0144
\n",
+ "
\n",
+ "
\n",
+ "
Time:
09:07:03
Pearson chi2:
5.00
\n",
+ "
\n",
+ "
\n",
+ "
No. Iterations:
6
Pseudo R-squ. (CS):
0.04355
\n",
+ "
\n",
+ "
\n",
+ "
Covariance Type:
nonrobust
\n",
+ "
\n",
+ "
\n",
+ "
\n",
+ "
\n",
+ "
coef
std err
z
P>|z|
[0.025
0.975]
\n",
+ "
\n",
+ "
\n",
+ "
Intercept
5.0850
7.477
0.680
0.496
-9.570
19.740
\n",
+ "
\n",
+ "
\n",
+ "
Temperature
-0.1156
0.115
-1.004
0.316
-0.341
0.110
\n",
+ "
\n",
+ "
"
+ ],
+ "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: Thu, 22 Aug 2024 Deviance: 3.0144\n",
+ "Time: 09:07:03 Pearson chi2: 5.00\n",
+ "No. Iterations: 6 Pseudo R-squ. (CS): 0.04355\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": 13,
+ "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'],\n",
+ " data[['Intercept','Temperature']],\n",
+ " family=sm.families.Binomial(sm.families.links.logit())).fit() # Added \"()\" after \"logit\"\n",
+ "\n",
+ "logmodel.summary()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0850$ 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": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "
Generalized Linear Model Regression Results
\n",
+ "
\n",
+ "
Dep. Variable:
Frequency
No. Observations:
23
\n",
+ "
\n",
+ "
\n",
+ "
Model:
GLM
Df Residuals:
21
\n",
+ "
\n",
+ "
\n",
+ "
Model Family:
Binomial
Df Model:
1
\n",
+ "
\n",
+ "
\n",
+ "
Link Function:
logit
Scale:
1.0000
\n",
+ "
\n",
+ "
\n",
+ "
Method:
IRLS
Log-Likelihood:
-23.526
\n",
+ "
\n",
+ "
\n",
+ "
Date:
Thu, 22 Aug 2024
Deviance:
18.086
\n",
+ "
\n",
+ "
\n",
+ "
Time:
09:09:20
Pearson chi2:
30.0
\n",
+ "
\n",
+ "
\n",
+ "
No. Iterations:
6
Pseudo R-squ. (CS):
0.2344
\n",
+ "
\n",
+ "
\n",
+ "
Covariance Type:
nonrobust
\n",
+ "
\n",
+ "
\n",
+ "
\n",
+ "
\n",
+ "
coef
std err
z
P>|z|
[0.025
0.975]
\n",
+ "
\n",
+ "
\n",
+ "
Intercept
5.0850
3.052
1.666
0.096
-0.898
11.068
\n",
+ "
\n",
+ "
\n",
+ "
Temperature
-0.1156
0.047
-2.458
0.014
-0.208
-0.023
\n",
+ "
\n",
+ "
"
+ ],
+ "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: Thu, 22 Aug 2024 Deviance: 18.086\n",
+ "Time: 09:09:20 Pearson chi2: 30.0\n",
+ "No. Iterations: 6 Pseudo R-squ. (CS): 0.2344\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": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n",
+ " family=sm.families.Binomial(sm.families.links.logit()), # Added \"()\" after \"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": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzZUlEQVR4nO3deXBUVf7+8aezE0KCBMiiAcIii+BCGDQoiwtBQNRyRlBkE3DMgCLGDWSURQUdlUFLiYhsEhdUlK84EchYsggoEogDhgKUYFA7k5+gJBqTdNL39weTlqazdGfhQHi/qrqKe+65957+5EIe7mqzLMsSAACAIX6mBwAAAM5thBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABglM9hZPPmzRo2bJhiY2Nls9m0Zs2aGpfZtGmTEhISFBISovbt2+uVV16pzVgBAEAj5HMY+e2333TJJZfopZde8qp/Tk6OhgwZor59+2r37t169NFHNWXKFK1evdrnwQIAgMbHVpcX5dlsNn3wwQe6+eabq+zzyCOP6MMPP9S+fftcbcnJyfrqq6+0ffv22m4aAAA0EgENvYHt27crKSnJrW3QoEFasmSJHA6HAgMDPZYpKSlRSUmJa9rpdOrYsWOKjIyUzWZr6CEDAIB6YFmWCgsLFRsbKz+/qk/GNHgYycvLU1RUlFtbVFSUysrK9NNPPykmJsZjmXnz5mn27NkNPTQAAHAaHDlyRBdccEGV8xs8jEjyOJpRcWaoqqMc06dPV0pKimv6+PHjatOmjXJyctSsWbN6GZNlWSosKtGmzZvVv18/BQSellKc1cocZdTLS9TKe9TKe9TKe9TKexW1GnTtAAUFBdXrugsLCxUfH1/j7+4G/wlFR0crLy/PrS0/P18BAQGKjIysdJng4GAFBwd7tLdo0ULh4eH1NrYIh0PnNQvVBTGtKz1dBHcO6uU1auU9auU9auU9auW9ilq1bNmy3mtVsb6aLrFo8OeMJCYmKiMjw61tw4YN6tWrFzsIAADwPYz8+uuvysrKUlZWlqQTt+5mZWUpNzdX0olTLGPGjHH1T05O1nfffaeUlBTt27dPS5cu1ZIlS/Tggw/WzzcAAABnNZ9P0+zcuVNXX321a7ri2o6xY8dq+fLlstvtrmAiSfHx8UpPT9f999+vl19+WbGxsXrxxRf15z//uR6GDwAAznY+h5EBAwaoukeTLF++3KOtf//+2rVrl6+bAgA0EuXl5XI4HKdtew6HQwEBASouLlZ5eflp2+7ZqC61CgwMlL+/f53HwCXGAIAGY1mW8vLy9Msvv5z27UZHR+vIkSM8n6oGda1V8+bNFR0dXac6E0YAAA2mIoi0bt1aoaGhpy0YOJ1O/frrrwoLC6v2YVuofa0sy1JRUZHy8/MlqdLnhnmLMAIAaBDl5eWuIFLVoxwaitPpVGlpqUJCQggjNahLrZo0aSLpxCM7WrduXetTNvyEAAANouIakdDQUMMjQUOq+PnW5ZogwggAoEFxzUbjVh8/X8IIAAAwijACAACMIowAAHCKcePGyWazeXy++eYb00NrlLibBgCASlx//fVatmyZW1urVq3cpktLS+v9TbfnIo6MAABQieDgYEVHR7t9rr32Wt1zzz1KSUlRy5YtNXDgQElSdna2hgwZorCwMEVFRWn06NH66aefXOv67bffNGbMGIWFhSkmJkbPP/+8BgwYoKlTp7r62Gw2rVmzxm0MzZs3d3uy+Q8//KARI0bovPPOU2RkpG666SYdPnzYNX/cuHG6+eab9dxzzykmJkaRkZGaPHmy250uJSUlevjhhxUXF6fg4GB17txZK1eulGVZ6tixo5577jm3Mezdu1d+fn769ttv617UKhBGAACnjWVZKiotOy2f30vLXX+u7jUmvlqxYoUCAgK0detWLVq0SHa7Xf3799ell16qnTt3at26dfrvf/+r4cOHu5Z56KGH9Omnn+qDDz7Qhg0btHHjRmVmZvq03aKiIl199dUKCwvT5s2b9dlnnyksLEzXX3+9SktLXf0+/fRTffvtt/r000+1YsUKLV++3C3QjBkzRm+//bZefPFF7du3TwsXLlTTpk1ls9k0fvx4j6NBS5cuVd++fdWhQ4faFcwLnKYBAJw2vzvK1e3x9ad9u9lzBik0yLdfeR999JHCwsJc04MHD5YkdezYUf/4xz9c7Y8//rh69uypuXPnutqWLl2quLg4HThwQLGxsVqyZIlef/1115GUFStW6IILLvBpPG+//bb8/Pz02muvuW6nXbZsmZo3b66NGzcqKSlJknTeeefppZdekr+/v7p06aKhQ4fqk08+0V133aUDBw7onXfeUUZGhq677jpJUrt27VRQUCBJuvPOO/X4449rx44d6t27txwOh9LS0vTss8/6NFZfEUYAAKjE1VdfrdTUVNd006ZNdfvtt6tXr15u/TIzM/Xpp5+6BZcK3377rX7//XeVlpYqMTHR1d6iRQt17tzZp/FkZmbqm2++UbNmzdzai4uL3U6hXHTRRW5PQo2JidGePXskSVlZWfL391f//v0r3UZMTIyGDh2qpUuXqnfv3vroo49UXFysW2+91aex+oowAgA4bZoE+it7zqAG347T6VRhQaGahTeTn5+fmgT6/pjypk2bqmPHjpW2n7qtYcOG6ZlnnvHoGxMTo4MHD3q1PZvN5nE66eRrPZxOpxISEvTGG294LHvyhbWBgYEe63U6nZL+eHx7dSZOnKjRo0frn//8p5YtW6YRI0Y0+FN0CSMAgNPGZrP5fLqkNpxOp8qC/BUaFNDg76bp2bOnVq9erXbt2ikgwPO7dezYUYGBgfr888/Vpk0bSdLPP/+sAwcOuB2haNWqlex2u2v64MGDKioqctvOqlWr1Lp1a4WHh9dqrD169JDT6dSmTZtcp2lONWTIEDVt2lSpqan6+OOPtXnz5lptyxdcwAoAQB1MnjxZx44d0+23364dO3bo0KFD2rBhg8aPH6/y8nKFhYVpwoQJeuihh/TJJ59o7969GjdunEdIuuaaa/TSSy9p165d2rlzp5KTk92Octxxxx1q2bKlbrrpJm3ZskU5OTnatGmT7rvvPn3//fdejbVdu3YaO3asxo8frzVr1ignJ0cbN27UBx984Orj7++vcePGafr06erYsaPb6aWGQhgBAKAOYmNjtXXrVpWXl2vQoEHq3r277rvvPkVERLgCx7PPPqt+/frpxhtv1HXXXaerrrpKCQkJbut5/vnnFRcXp379+mnkyJF68MEH3U6PhIaGavPmzWrTpo1uueUWde3aVePHj9fvv//u05GS1NRU/eUvf9GkSZPUpUsX3X333W5HYCRpwoQJKi0t1fjx4+tQGe9xmgYAgFOcfCvsyTZu3Fhpe6dOnfT+++9Xub6wsDCtXLlSK1eudLX961//cusTGxur9evd7zT65Zdf3Kajo6O1YsUKn8a9YMECt+mQkBDNnz9f8+fPl3TilFbF3TQV7Ha7AgICNGbMmCq3VZ8IIwAAQNKJB6IdOXJEjz32mIYPH66oqKjTsl1O0wAAAEnSW2+9pc6dO+v48eNuz1JpaBwZAQDAgKpO+Zg0btw4jRs37rRvlyMjAADAKMIIAKBB1ed7YXDmqY+fL2EEANAgKp6Rcepto2hcKn6+pz751RdcMwIAaBD+/v5q3ry58vPzJZ14TkbFC94amtPpVGlpqYqLixv8Caxnu9rWyrIsFRUVKT8/X82bN3d7H46vCCMAgAYTHR0tSa5AcrpYlqXff/9dTZo0OW0B6GxV11o1b97c9XOuLcIIAKDB2Gw2xcTEqHXr1m4vfWtoDodDmzdvVr9+/ep0+uBcUJdaBQYG1umISAXCCACgwfn7+9fLLy1ftldWVqaQkBDCSA3OhFpxIg0AABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRtQojCxcuVHx8vEJCQpSQkKAtW7ZU2/+NN97QJZdcotDQUMXExOjOO+/U0aNHazVgAADQuPgcRlatWqWpU6dqxowZ2r17t/r27avBgwcrNze30v6fffaZxowZowkTJujrr7/Wu+++qy+//FITJ06s8+ABAMDZz+cwMn/+fE2YMEETJ05U165dtWDBAsXFxSk1NbXS/p9//rnatWunKVOmKD4+XldddZXuvvtu7dy5s86DBwAAZ78AXzqXlpYqMzNT06ZNc2tPSkrStm3bKl2mT58+mjFjhtLT0zV48GDl5+frvffe09ChQ6vcTklJiUpKSlzTBQUFkiSHwyGHw+HLkKtVsa76XGdjRr28R628R628R628R62815C18nadNsuyLG9X+uOPP+r888/X1q1b1adPH1f73LlztWLFCu3fv7/S5d577z3deeedKi4uVllZmW688Ua99957CgwMrLT/rFmzNHv2bI/2N998U6Ghod4OFwAAGFRUVKSRI0fq+PHjCg8Pr7KfT0dGKthsNrdpy7I82ipkZ2drypQpevzxxzVo0CDZ7XY99NBDSk5O1pIlSypdZvr06UpJSXFNFxQUKC4uTklJSdV+GV85HA5lZGRo4MCBVQYj/IF6eY9aeY9aeY9aeY9aea8ha1VxZqMmPoWRli1byt/fX3l5eW7t+fn5ioqKqnSZefPm6corr9RDDz0kSbr44ovVtGlT9e3bV08++aRiYmI8lgkODlZwcLBHe2BgYIPsVA213saKenmPWnmPWnmPWnmPWnmvIWrl7fp8uoA1KChICQkJysjIcGvPyMhwO21zsqKiIvn5uW/G399f0okjKgAA4Nzm8900KSkpeu2117R06VLt27dP999/v3Jzc5WcnCzpxCmWMWPGuPoPGzZM77//vlJTU3Xo0CFt3bpVU6ZMUe/evRUbG1t/3wQAAJyVfL5mZMSIETp69KjmzJkju92u7t27Kz09XW3btpUk2e12t2eOjBs3ToWFhXrppZf0wAMPqHnz5rrmmmv0zDPP1N+3AAAAZ61aXcA6adIkTZo0qdJ5y5cv92i79957de+999ZmUwAAoJHj3TQAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIyqVRhZuHCh4uPjFRISooSEBG3ZsqXa/iUlJZoxY4batm2r4OBgdejQQUuXLq3VgAEAQOMS4OsCq1at0tSpU7Vw4UJdeeWVWrRokQYPHqzs7Gy1adOm0mWGDx+u//73v1qyZIk6duyo/Px8lZWV1XnwAADg7OdzGJk/f74mTJigiRMnSpIWLFig9evXKzU1VfPmzfPov27dOm3atEmHDh1SixYtJEnt2rWr26gBAECj4VMYKS0tVWZmpqZNm+bWnpSUpG3btlW6zIcffqhevXrpH//4h1auXKmmTZvqxhtv1BNPPKEmTZpUukxJSYlKSkpc0wUFBZIkh8Mhh8Phy5CrVbGu+lxnY0a9vEetvEetvEetvEetvNeQtfJ2nT6FkZ9++knl5eWKiopya4+KilJeXl6lyxw6dEifffaZQkJC9MEHH+inn37SpEmTdOzYsSqvG5k3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/Xw+TSNJNpvNbdqyLI+2Ck6nUzabTW+88YYiIiIknTjV85e//EUvv/xypUdHpk+frpSUFNd0QUGB4uLilJSUpPDw8NoMuVIOh0MZGRkaOHCgAgMD6229jRX18h618h618h618h618l5D1qrizEZNfAojLVu2lL+/v8dRkPz8fI+jJRViYmJ0/vnnu4KIJHXt2lWWZen7779Xp06dPJYJDg5WcHCwR3tgYGCD7FQNtd7Ginp5j1p5j1p5j1p5j1p5ryFq5e36fLq1NygoSAkJCR6HcjIyMtSnT59Kl7nyyiv1448/6tdff3W1HThwQH5+frrgggt82TwAAGiEfH7OSEpKil577TUtXbpU+/bt0/3336/c3FwlJydLOnGKZcyYMa7+I0eOVGRkpO68805lZ2dr8+bNeuihhzR+/PgqL2AFAADnDp+vGRkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7Xbm5ua7+YWFhysjI0L333qtevXopMjJSw4cP15NPPll/3wIAAJy1anUB66RJkzRp0qRK5y1fvtyjrUuXLlzRDAAAKsW7aQAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGFWrMLJw4ULFx8crJCRECQkJ2rJli1fLbd26VQEBAbr00ktrs1kAANAI+RxGVq1apalTp2rGjBnavXu3+vbtq8GDBys3N7fa5Y4fP64xY8bo2muvrfVgAQBA4+NzGJk/f74mTJigiRMnqmvXrlqwYIHi4uKUmppa7XJ33323Ro4cqcTExFoPFgAAND4BvnQuLS1VZmampk2b5taelJSkbdu2VbncsmXL9O233yotLU1PPvlkjdspKSlRSUmJa7qgoECS5HA45HA4fBlytSrWVZ/rbMyol/eolfeolfeolfeolfcaslbertOnMPLTTz+pvLxcUVFRbu1RUVHKy8urdJmDBw9q2rRp2rJliwICvNvcvHnzNHv2bI/2DRs2KDQ01JcheyUjI6Pe19mYUS/vUSvvUSvvUSvvUSvvNUStioqKvOrnUxipYLPZ3KYty/Jok6Ty8nKNHDlSs2fP1oUXXuj1+qdPn66UlBTXdEFBgeLi4pSUlKTw8PDaDLlSDodDGRkZGjhwoAIDA+ttvY0V9fIetfIetfIetfIetfJeQ9aq4sxGTXwKIy1btpS/v7/HUZD8/HyPoyWSVFhYqJ07d2r37t265557JElOp1OWZSkgIEAbNmzQNddc47FccHCwgoODPdoDAwMbZKdqqPU2VtTLe9TKe9TKe9TKe9TKew1RK2/X59MFrEFBQUpISPA4lJORkaE+ffp49A8PD9eePXuUlZXl+iQnJ6tz587KysrS5Zdf7svmAQBAI+TzaZqUlBSNHj1avXr1UmJiol599VXl5uYqOTlZ0olTLD/88INef/11+fn5qXv37m7Lt27dWiEhIR7tAADg3ORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa7vcZnjgAAAFSo1QWskyZN0qRJkyqdt3z58mqXnTVrlmbNmlWbzQIAgEaId9MAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjavXWXgCnX7nT0o6cY8ovLFbrZiHqHd9C/n4208PCOY79EvWBMAKcBdbttWv22mzZjxe72mIiQjRzWDdd3z3G4MhwLmO/RH3hNA1whlu3166/pe1y+wdfkvKOF+tvabu0bq/d0MhwLmO/RH0ijABnsHKnpdlrs2VVMq+ibfbabJU7K+sBNAz2S9Q3wghwBtuRc8zjf54nsyTZjxdrR86x0zconPPYL1HfCCPAGSy/sOp/8GvTD6gP7Jeob4QR4AzWullIvfYD6gP7JeobYQQ4g/WOb6GYiBBVdaOkTSfuXugd3+J0DgvnOPZL1DfCCHAG8/ezaeawbpLk8Q9/xfTMYd14rgNOK/ZL1DfCCHCGu757jFJH9VR0hPsh7+iIEKWO6snzHGAE+yXqEw89A84C13eP0cBu0TzpEmcU9kvUF8IIcJbw97MpsUOk6WEAbtgvUR84TQMAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKNqFUYWLlyo+Ph4hYSEKCEhQVu2bKmy7/vvv6+BAweqVatWCg8PV2JiotavX1/rAQMAgMbF5zCyatUqTZ06VTNmzNDu3bvVt29fDR48WLm5uZX237x5swYOHKj09HRlZmbq6quv1rBhw7R79+46Dx4AAJz9fA4j8+fP14QJEzRx4kR17dpVCxYsUFxcnFJTUyvtv2DBAj388MP605/+pE6dOmnu3Lnq1KmT1q5dW+fBAwCAs1+AL51LS0uVmZmpadOmubUnJSVp27ZtXq3D6XSqsLBQLVq0qLJPSUmJSkpKXNMFBQWSJIfDIYfD4cuQq1WxrvpcZ2NGvbxHrbxHrbxHrbxHrbzXkLXydp0+hZGffvpJ5eXlioqKcmuPiopSXl6eV+t4/vnn9dtvv2n48OFV9pk3b55mz57t0b5hwwaFhob6MmSvZGRk1Ps6GzPq5T1q5T1q5T1q5T1q5b2GqFVRUZFX/XwKIxVsNpvbtGVZHm2VeeuttzRr1iz93//9n1q3bl1lv+nTpyslJcU1XVBQoLi4OCUlJSk8PLw2Q66Uw+FQRkaGBg4cqMDAwHpbb2NFvbxHrbxHrbxHrbxHrbzXkLWqOLNRE5/CSMuWLeXv7+9xFCQ/P9/jaMmpVq1apQkTJujdd9/VddddV23f4OBgBQcHe7QHBgY2yE7VUOttrKiX96iV96iV96iV96iV9xqiVt6uz6cLWIOCgpSQkOBxKCcjI0N9+vSpcrm33npL48aN05tvvqmhQ4f6skkAANDI+XyaJiUlRaNHj1avXr2UmJioV199Vbm5uUpOTpZ04hTLDz/8oNdff13SiSAyZswYvfDCC7riiitcR1WaNGmiiIiIevwqAADgbORzGBkxYoSOHj2qOXPmyG63q3v37kpPT1fbtm0lSXa73e2ZI4sWLVJZWZkmT56syZMnu9rHjh2r5cuX1/0bAACAs1qtLmCdNGmSJk2aVOm8UwPGxo0ba7MJAABwjqhVGAFw7ih3WtqRc0z5hcVq3SxEveNbyN/P5vV8E87EMdVVaZlTadsPK1LSyu2HNapPBwUF8HoxNA6EEQBVWrfXrtlrs2U/Xuxqi4kI0cxh3XR995ga55twJo6prualZ2vxlhwF+ln6R2/pmfX79eTHB3RX33hNH9LN9PCAOiNWA6jUur12/S1tl9svdUnKO16sv6Xt0rz07Grnr9trP53DlVTzmE2Mqa7mpWdr0eYcOS33dqclLdqco3np2WYGBtQjwggAD+VOS7PXZsuqZJ71v8/iLTlVzpek2WuzVX7qb9AGVNOYTYyprkrLnFq8JafaPou35Ki0zHmaRgQ0DMIIAA87co55HF04VXW/0y1J9uPF2pFzrH4HVo2axmxiTHW1cvvhaussnfg5rNx++LSMB2gohBEAHvILqw8ip3s99bmt0zmmuvrumHfv9fC2H3CmIowA8NC6WcgZtZ763NbpHFNdtW3h3YtBve0HnKkIIwA89I5voZiIEFV3M6yfTVXOt+nEHSy941s0wOgqV9OYTYyprkYntlNNdyT72U70A85mhBEAHvz9bJo57MQto6f+LrT973NX3/gq50vSzGHdTuuzPWoas4kx1VVQgJ+rzlW5q288zxvBWY89GEClru8eo9RRPRUd4X5aIzoiRKmjemr6kG7VzjfxTI+axnw2Pmdk+pBuurtfvMcREj+bdHc/njOCxoGHngGo0vXdYzSwW3SVTzOtaf6ZOOaz0fQh3fRAUhelbftW+jlbjwzqzBNY0agQRgBUy9/PpsQOkbWeb8KZOKa6Cgrw0+jEdkpPz9boxHYKJIigEWFvBgAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUQGmBwAAZ5Nyp6UdOceUX1is1s1C1Du+hfz9bJKk30vLNTc9W4ePFqldZKgeHdJNTYL8vVq2unmSVFrmVNr2w4qUtHL7YY3q00FBAd79f7Kmddc0v7brLi1zauX2w/ruWJHatgjV6MR2Z/yYYUatwsjChQv17LPPym6366KLLtKCBQvUt2/fKvtv2rRJKSkp+vrrrxUbG6uHH35YycnJtR40AJiwbq9ds9dmy3682NUWExGimcO6afWu75WRne9q33JQWvl5rgZ2a63FY/5U7bKSqpx3ffcYzUvP1uItOQr0s/SP3tIz6/fryY8P6K6+8Zo+pFutx3x995ga59d23btzf9biLTlyWn/0fyp93xk95pqWRcPxOYysWrVKU6dO1cKFC3XllVdq0aJFGjx4sLKzs9WmTRuP/jk5ORoyZIjuuusupaWlaevWrZo0aZJatWqlP//5z/XyJQCgoa3ba9ff0nbJOqU973ixktN2VblcRna+bnxpi/Z8X+DTsnnHi/W3tF26rltrt5BTwWlJizbnSFKVv9yrG/Pf0nbpr/3i9ermnCrnp47qWeUv6NrU40wec03LomH5fM3I/PnzNWHCBE2cOFFdu3bVggULFBcXp9TU1Er7v/LKK2rTpo0WLFigrl27auLEiRo/fryee+65Og8eAE6Hcqel2WuzPX6JSaq07VT/qSSI1LSs9b9PZUHkZIu35Ki0zOnRXtOYrf8tW924Zq/NVrnTs0dd63GmjrmqZdHwfDoyUlpaqszMTE2bNs2tPSkpSdu2bat0me3btyspKcmtbdCgQVqyZIkcDocCAwM9likpKVFJSYlr+vjx45KkY8eOyeFw+DLkajkcDhUVFeno0aOVjgPuqJf3qJX3zoZaZX73s/7f0aPGL7ILcFoqKnIqwOGncucf1zgszvhKt/V2PzLt7Zir+x/p/zv6mz7J+lYJbc+r1bqr09Bjvjg2zG2/8mbdVX3fxq4h/w4WFhZKkiyrhpBn+eCHH36wJFlbt251a3/qqaesCy+8sNJlOnXqZD311FNubVu3brUkWT/++GOly8ycObMiBPPhw4cPHz58zvLPkSNHqs0XtQq2Npv7VceWZXm01dS/svYK06dPV0pKimva6XTq2LFjioyMrHY7viooKFBcXJyOHDmi8PDweltvY0W9vEetvEetvEetvEetvNeQtbIsS4WFhYqNja22n09hpGXLlvL391deXp5be35+vqKioipdJjo6utL+AQEBioyMrHSZ4OBgBQcHu7U1b97cl6H6JDw8nJ3VB9TLe9TKe9TKe9TKe9TKew1Vq4iIiBr7+HQBa1BQkBISEpSRkeHWnpGRoT59+lS6TGJiokf/DRs2qFevXmfs+WEAAHD6+Hw3TUpKil577TUtXbpU+/bt0/3336/c3FzXc0OmT5+uMWPGuPonJyfru+++U0pKivbt26elS5dqyZIlevDBB+vvWwAAgLOWz9eMjBgxQkePHtWcOXNkt9vVvXt3paenq23btpIku92u3NxcV//4+Hilp6fr/vvv18svv6zY2Fi9+OKLZ8QzRoKDgzVz5kyPU0KoHPXyHrXyHrXyHrXyHrXy3plQK5tl1XS/DQAAQMPhRXkAAMAowggAADCKMAIAAIwijAAAAKPOiTCSmpqqiy++2PVAl8TERH388ceu+ZZladasWYqNjVWTJk00YMAAff311wZHfGaYN2+ebDabpk6d6mqjVn+YNWuWbDab2yc6Oto1n1q5++GHHzRq1ChFRkYqNDRUl156qTIzM13zqdcJ7dq189ivbDabJk+eLIk6naysrEx///vfFR8fryZNmqh9+/aaM2eOnM4/XsJHvf5QWFioqVOnqm3btmrSpIn69OmjL7/80jXfaK1qeB1No/Dhhx9a//rXv6z9+/db+/fvtx599FErMDDQ2rt3r2VZlvX0009bzZo1s1avXm3t2bPHGjFihBUTE2MVFBQYHrk5O3bssNq1a2ddfPHF1n333edqp1Z/mDlzpnXRRRdZdrvd9cnPz3fNp1Z/OHbsmNW2bVtr3Lhx1hdffGHl5ORY//73v61vvvnG1Yd6nZCfn++2T2VkZFiSrE8//dSyLOp0sieffNKKjIy0PvroIysnJ8d69913rbCwMGvBggWuPtTrD8OHD7e6detmbdq0yTp48KA1c+ZMKzw83Pr+++8tyzJbq3MijFTmvPPOs1577TXL6XRa0dHR1tNPP+2aV1xcbEVERFivvPKKwRGaU1hYaHXq1MnKyMiw+vfv7woj1MrdzJkzrUsuuaTSedTK3SOPPGJdddVVVc6nXlW77777rA4dOlhOp5M6nWLo0KHW+PHj3dpuueUWa9SoUZZlsV+drKioyPL397c++ugjt/ZLLrnEmjFjhvFanROnaU5WXl6ut99+W7/99psSExOVk5OjvLw8JSUlufoEBwerf//+2rZtm8GRmjN58mQNHTpU1113nVs7tfJ08OBBxcbGKj4+XrfddpsOHTokiVqd6sMPP1SvXr106623qnXr1rrsssu0ePFi13zqVbnS0lKlpaVp/Pjxstls1OkUV111lT755BMdOHBAkvTVV1/ps88+05AhQySxX52srKxM5eXlCgkJcWtv0qSJPvvsM+O1OmfCyJ49exQWFqbg4GAlJyfrgw8+ULdu3Vwv8Tv1RX9RUVEeL/g7F7z99tvatWuX5s2b5zGPWrm7/PLL9frrr2v9+vVavHix8vLy1KdPHx09epRaneLQoUNKTU1Vp06dtH79eiUnJ2vKlCl6/fXXJbFvVWXNmjX65ZdfNG7cOEnU6VSPPPKIbr/9dnXp0kWBgYG67LLLNHXqVN1+++2SqNfJmjVrpsTERD3xxBP68ccfVV5errS0NH3xxRey2+3Ga+Xz4+DPVp07d1ZWVpZ++eUXrV69WmPHjtWmTZtc8202m1t/y7I82hq7I0eO6L777tOGDRs80vPJqNUJgwcPdv25R48eSkxMVIcOHbRixQpdccUVkqhVBafTqV69emnu3LmSpMsuu0xff/21UlNT3d5lRb3cLVmyRIMHD/Z4/Tp1OmHVqlVKS0vTm2++qYsuukhZWVmaOnWqYmNjNXbsWFc/6nXCypUrNX78eJ1//vny9/dXz549NXLkSO3atcvVx1StzpkjI0FBQerYsaN69eqlefPm6ZJLLtELL7zguvvh1OSXn5/vkRAbu8zMTOXn5yshIUEBAQEKCAjQpk2b9OKLLyogIMBVD2pVuaZNm6pHjx46ePAg+9UpYmJi1K1bN7e2rl27ut5jRb08fffdd/r3v/+tiRMnutqok7uHHnpI06ZN02233aYePXpo9OjRuv/++11HdqmXuw4dOmjTpk369ddfdeTIEe3YsUMOh0Px8fHGa3XOhJFTWZalkpIS1w8hIyPDNa+0tFSbNm1Snz59DI7w9Lv22mu1Z88eZWVluT69evXSHXfcoaysLLVv355aVaOkpET79u1TTEwM+9UprrzySu3fv9+t7cCBA64XbFIvT8uWLVPr1q01dOhQVxt1cldUVCQ/P/dfY/7+/q5be6lX5Zo2baqYmBj9/PPPWr9+vW666SbztWrwS2TPANOnT7c2b95s5eTkWP/5z3+sRx991PLz87M2bNhgWdaJ25kiIiKs999/39qzZ491++23n7O3fp3q5LtpLItaneyBBx6wNm7caB06dMj6/PPPrRtuuMFq1qyZdfjwYcuyqNXJduzYYQUEBFhPPfWUdfDgQeuNN96wQkNDrbS0NFcf6vWH8vJyq02bNtYjjzziMY86/WHs2LHW+eef77q19/3337datmxpPfzww64+1OsP69atsz7++GPr0KFD1oYNG6xLLrnE6t27t1VaWmpZltlanRNhZPz48Vbbtm2toKAgq1WrVta1117rCiKWdeL2r5kzZ1rR0dFWcHCw1a9fP2vPnj0GR3zmODWMUKs/VNyDHxgYaMXGxlq33HKL9fXXX7vmUyt3a9eutbp3724FBwdbXbp0sV599VW3+dTrD+vXr7ckWfv37/eYR53+UFBQYN13331WmzZtrJCQEKt9+/bWjBkzrJKSElcf6vWHVatWWe3bt7eCgoKs6Ohoa/LkydYvv/zimm+yVjbLsqyGP/4CAABQuXP2mhEAAHBmIIwAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAjRCNput2k/FK+kbkwEDBmjq1KmmhwGgFgJMDwBA/bPb7a4/r1q1So8//rjbi+qaNGliYli14nA4FBgY2Gi3B4AjI0CjFB0d7fpERETIZrO5tW3evFkJCQkKCQlR+/btNXv2bJWVlbmWt9lsWrRokW644QaFhoaqa9eu2r59u7755hsNGDBATZs2VWJior799lvXMrNmzdKll16qRYsWKS4uTqGhobr11lv1yy+/uI1t2bJl6tq1q0JCQtSlSxctXLjQNe/w4cOy2Wx65513NGDAAIWEhCgtLU1Hjx7V7bffrgsuuEChoaHq0aOH3nrrLddy48aN06ZNm/TCCy+4jv4cPnxYy5cvV/Pmzd22v2bNGtlsNo9xL126VO3bt1dwcLAsy9Lx48f117/+Va1bt1Z4eLiuueYaffXVV/X0EwJwMsIIcI5Zv369Ro0apSlTpig7O1uLFi3S8uXL9dRTT7n1e+KJJzRmzBhlZWWpS5cuGjlypO6++25Nnz5dO3fulCTdc889bst88803euedd7R27VqtW7dOWVlZmjx5smv+4sWLNWPGDD311FPat2+f5s6dq8cee0wrVqxwW88jjzyiKVOmaN++fRo0aJCKi4uVkJCgjz76SHv37tVf//pXjR49Wl988YUk6YUXXlBiYqLuuusu2e122e12xcXFeV2TinGvXr1aWVlZkqShQ4cqLy9P6enpyszMVM+ePXXttdfq2LFjXq8XgJdOy+v4ABizbNkyKyIiwjXdt29fa+7cuW59Vq5cacXExLimJVl///vfXdPbt2+3JFlLlixxtb311ltWSEiIa3rmzJmWv7+/deTIEVfbxx9/bPn5+Vl2u92yLMuKi4uz3nzzTbdtP/HEE1ZiYqJlWZaVk5NjSbIWLFhQ4/caMmSI9cADD7imT33DdGXf3bIs64MPPrBO/qdv5syZVmBgoJWfn+9q++STT6zw8HCruLjYbdkOHTpYixYtqnFsAHzDNSPAOSYzM1Nffvml25GQ8vJyFRcXq6ioSKGhoZKkiy++2DU/KipKktSjRw+3tuLiYhUUFCg8PFyS1KZNG11wwQWuPomJiXI6ndq/f7/8/f115MgRTZgwQXfddZerT1lZmSIiItzG2KtXL7fp8vJyPf3001q1apV++OEHlZSUqKSkRE2bNq1rOSRJbdu2VatWrVzTmZmZ+vXXXxUZGenW7/fff3c7NQWgfhBGgHOM0+nU7Nmzdcstt3jMCwkJcf355Is4K66xqKzN6XRWua2KPjabzdVv8eLFuvzyy936+fv7u02fGjKef/55/fOf/9SCBQvUo0cPNW3aVFOnTlVpaWnVX1SSn5+fLMtya3M4HB79Tt2e0+lUTEyMNm7c6NH31GtQANQdYQQ4x/Ts2VP79+9Xx44d633dubm5+vHHHxUbGytJ2r59u/z8/HThhRcqKipK559/vg4dOqQ77rjDp/Vu2bJFN910k0aNGiXpRFg4ePCgunbt6uoTFBSk8vJyt+VatWqlwsJC/fbbb67AUXFNSHV69uypvLw8BQQEqF27dj6NFYDvCCPAOebxxx/XDTfcoLi4ON16663y8/PTf/7zH+3Zs0dPPvlkndYdEhKisWPH6rnnnlNBQYGmTJmi4cOHKzo6WtKJO1emTJmi8PBwDR48WCUlJdq5c6d+/vlnpaSkVLnejh07avXq1dq2bZvOO+88zZ8/X3l5eW5hpF27dvriiy90+PBhhYWFqUWLFrr88ssVGhqqRx99VPfee6927Nih5cuX1/g9rrvuOiUmJurmm2/WM888o86dO+vHH39Uenq6br75Zo/TSADqhrtpgHPMoEGD9NFHHykjI0N/+tOfdMUVV2j+/Plq27ZtndfdsWNH3XLLLRoyZIiSkpLUvXt3t1t3J06cqNdee03Lly9Xjx491L9/fy1fvlzx8fHVrvexxx5Tz549NWjQIA0YMEDR0dG6+eab3fo8+OCD8vf3V7du3dSqVSvl5uaqRYsWSktLU3p6uut24FmzZtX4PWw2m9LT09WvXz+NHz9eF154oW677TYdPnzYdf0MgPpjs049oQoAtTBr1iytWbPGq9MgAHAyjowAAACjCCMAAMAoTtMAAACjODICAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjPr/yAWkb36Wd7EAAAAASUVORK5CYII=\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": [
+ "## 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": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG6CAYAAAALTELXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABThUlEQVR4nO3deXxU1f0//te9s2ayJ5CFHdlCIASQVUAQRf2BReWDrViQakVbF4q0gv26QNWKrSAtAuICYl0qKopVEUWtdWMXQYSwCYFA9j2Z/d7z+2OSIUMCJDM3mYT7ej4ekeTeO3PPvDMmr5xz7rmSEEKAiIiISKfkcDeAiIiIKJwYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI11pVGFq5ciVmzJhx3mNKS0vxxz/+EUOHDsXQoUPxyCOPwG63t1ALiYiI6GLTasLQ2rVrsWzZsgseN3v2bJw8edJ//Lfffou//OUvLdBCIiIiuhgZw92A/Px8PPTQQ9i1axe6d+9+3mN3796N7du3Y+PGjejRowcA4LHHHsMdd9yBuXPnIjk5uSWaTERERBeRsPcM/fTTT4iNjcV//vMfZGZmnvfYnTt3on379v4gBADDhg2DJEnYtWtXczeViIiILkJh7xkaP348xo8f36hj8/PzkZqaGrDNbDYjLi4Oubm5zdE8IiIiusiFvWeoKRwOB8xmc73tFosFLpcr6OcVQoTSLCIiImrDwt4z1BRWqxVut7vedpfLBZvNFvTzSpKEigoHFEUNpXm6ZzDIiImJYC1DxDpqh7XUDmupDdZRO7GxEZBlbfp02lQYSklJwWeffRawze12o6ysLOTJ04qiwuvlG1MLrKU2WEftsJbaYS21wTqGTstBnTY1TDZ06FDk5eUhOzvbv23btm0AgMGDB4erWURERNSGteowpCgKCgsL4XQ6AQCZmZkYPHgw7r//fuzduxdbt27FggULcMMNN/CyeiIiIgpKqw5Dubm5GD16NDZu3AjAN7dn+fLl6NSpE2bOnIk5c+bg8ssvx8KFC8PbUCIiImqzJMFLqQAApaXVHL8NkdEoIz4+krUMEeuoHdZSO6ylNlhH7SQkRMJg0KZPp1X3DBERERE1N4YhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0jWGISIiItI1hiEiIiLSNYYhIiIi0rWwhyFVVbFs2TKMGTMGmZmZuP3225GdnX3O4wsLCzF37lwMHz4cw4cPxx/+8Afk5eW1YIuJiIjoYhL2MLRy5Uq8+eabeOKJJ7Bu3TpIkoRZs2bB7XY3ePz999+P3NxcvPzyy3j55ZeRl5eHu+++u4VbTURERBeLsIYht9uNNWvW4L777sPYsWORlpaGpUuXIj8/H5s3b653fEVFBXbs2IFZs2YhPT0d6enpuPPOO/HTTz+htLQ0DK+AiIiI2jpjOE+elZWF6upqjBgxwr8tJiYG6enp2LFjByZNmhRwvMVigc1mw4YNGzBs2DAAwPvvv49u3bohNjY2pLYYDGHvJGvzamvIWoaGddQOa6kd1lIbrKN2JEm75wprGKqd65OamhqwPSkpCbm5ufWOt1gs+Otf/4rHHnsMQ4YMgSRJaN++PV577TXIcmhvrJiYiJAeT2ewltpgHbXDWmqHtdQG69i6hDUMORwOAIDZbA7YbrFYUF5eXu94IQQOHjyIQYMG4Y477oCiKFi6dCnuuece/Pvf/0ZUVFTQbamocEBR1KAfT76/dGJiIljLELGO2mEttcNaaoN11E5sbETIHSG1whqGrFYrAN/codrPAcDlciEion5q/uijj/DGG2/gv//9rz/4rFq1CldccQXWr1+PmTNnBt0WRVHh9fKNqQXWUhuso3ZYS+2wltpgHUMnhHbPFdZBy9rhsYKCgoDtBQUFSElJqXf8rl270L1794AeoNjYWHTv3h3Hjx9v1rYSERHRxSmsYSgtLQ1RUVHYtm2bf1tFRQX279+PIUOG1Ds+NTUV2dnZcLlc/m0OhwM5OTno2rVri7SZiIiILi5hDUNmsxnTp0/H4sWL8fnnnyMrKwv3338/UlJSMGHCBCiKgsLCQjidTgDADTfcAACYM2cOsrKy/MebzWZMmTIljK+EiIiI2qqwX9s3e/ZsTJ06FQ8//DCmTZsGg8GA1atXw2w2Izc3F6NHj8bGjRsB+K4ye+ONNyCEwMyZM3HbbbfBZDLh3//+N2JiYsL8SoiIiKgtkoTQcgpS21VaWs3JbCEyGmXEx0eyliFiHbXDWmqHtdQG66idhIRIzdZrCnvPEBEREVE4MQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRrjEMERERka4xDBEREZGuMQwRERGRroU9DKmqimXLlmHMmDHIzMzE7bffjuzs7HMe7/F4sGTJEowZMwYDBw7E9OnTceDAgRZsMREREV1Mwh6GVq5ciTfffBNPPPEE1q1bB0mSMGvWLLjd7gaPX7hwId555x08/vjjWL9+PeLi4jBr1ixUVla2cMuJiIjoYhDWMOR2u7FmzRrcd999GDt2LNLS0rB06VLk5+dj8+bN9Y4/efIk3nnnHSxatAjjxo1Djx498OSTT8JsNmPfvn1heAVERETU1hnDefKsrCxUV1djxIgR/m0xMTFIT0/Hjh07MGnSpIDjv/nmG8TExODyyy8POP6LL74IuS0GQ9g7ydq82hqylqFhHbXDWmqHtdQG66gdSdLuucIahvLy8gAAqampAduTkpKQm5tb7/jjx4+jc+fO+PTTT/HCCy8gPz8f6enpePDBB9GjR4+Q2hITExHS4+kM1lIbrKN2WEvtsJbaYB1bl7CGIYfDAQAwm80B2y0WC8rLy+sdX1VVhRMnTmDlypWYN28eYmJi8Nxzz+GWW27Bxo0bkZiYGHRbKiocUBQ16MeT7y+dmJgI1jJErKN2WEvtsJbaYB21ExsbAVnWpoctqDC0YsUKTJkypV6PTlNZrVYAvrlDtZ8DgMvlQkRE/dRsMplQWVmJpUuX+nuCli5dirFjx+K9997DHXfcEXRbFEWF18s3phZYS22wjtphLbXDWmqDdQydENo9V1CR6pVXXsGVV16J2267DR988AFcLldQJ68NUwUFBQHbCwoKkJKSUu/4lJQUGI3GgCExq9WKzp07IycnJ6g2EBERkb4FFYa++eYbLF68GCaTCQ8++CBGjRqFRx99FLt3727S86SlpSEqKgrbtm3zb6uoqMD+/fsxZMiQescPGTIEXq8XP/74o3+b0+nEyZMn0bVr12BeChEREelcUMNkZrMZEydOxMSJE1FQUID//Oc/+Pjjj/H222+jW7dumDJlCqZMmXLBOTxmsxnTp0/H4sWLkZCQgI4dO+Lpp59GSkoKJkyYAEVRUFJSgujoaFitVgwZMgSXXXYZ5s+fj8ceewxxcXFYtmwZDAYDrr/++qAKQERERPoW8syjpKQk3Hrrrfjd736HIUOG4NixY3jmmWcwduxYPPLII6iqqjrv42fPno2pU6fi4YcfxrRp02AwGLB69WqYzWbk5uZi9OjR2Lhxo//4Z599FsOGDcO9996LqVOnoqqqCv/617+QkJAQ6kshIiIiHZKECH4K0vbt2/H+++/jk08+gd1ux4gRIzB16lSMHTsW//vf//DYY48hIyMDL774opZtbhalpdWczBYio1FGfHwkaxki1lE7rKV2WEttsI7aSUiI1Gy9pqCGyZYuXYoPPvgAubm5SE1NxW9+8xtMmTIFHTp08B8zceJEHDx4EP/61780aSgRERFRcwgqDL388su46qqr8Pjjj+Oyyy6DdI5lIDMyMjBnzpxQ2kdERETUrIIKQ19//TViY2NRWFjoD0Ll5eXIzc1FWlqa/7irrrpKm1YSERERNZOgBttkWcZtt92GGTNm+Lft2bMHN9xwA+6++27/ytJERERErV1QYejpp5/G4cOHMXfuXP+2ESNGYOXKldi3bx+WLVumWQOJiIiImlNQYeiLL77A/PnzcfXVV/u3mc1mjB8/HnPnzsXHH3+sWQOJiIiImlNQYai6uhoxMTEN7ktMTERpaWlIjSIiIiJqKUGFoX79+mH9+vUN7nv33XfRp0+fkBpFRERE1FKCuprs97//PWbNmoUpU6ZgwoQJSExMRElJCT7//HP89NNPWLVqldbtJCIiImoWQYWhUaNG4bnnnsOyZcuwbNkyCCEgSRL69u2LlStX4vLLL9e6nURERETNIqgwBABjx47F2LFj4XK5UFZWhujoaNhsNi3bRkRERNTsgg5DgG+hRYfDAVVVUVZWhrKyMv++urfmICIiImqtggpDx48fx4MPPog9e/ac85gDBw4E3SgiIiKilhJUGHr88cdx/Phx3HvvvUhJSYEsa3PXWCIiIqKWFlQY2rlzJ/7617/iuuuu07o9RERERC0qqC6dqKgoxMbGat0WIiIiohYXVBi6/vrr8frrr0MIoXV7iIiIiFpUUMNkERER2LVrFyZMmICMjAxYrdaA/ZIk4cknn9SkgURERETNKagw9N577yE6OhqqqjZ4RZkkSSE3jIiIiKglBBWGvvjiC63bQURERBQWIV0Tr6oqsrKy8NVXX6Gqqipg0UUiIiKitiDoFajff/99LFmyBAUFBZAkCe+88w6effZZmEwmLFmyBGazWct2EhERETWLoHqGNm7ciPnz52PEiBFYunSp/6qyq6++Gl999RVWrlypaSOJiIiImktQPUOrVq3CzTffjIULF0JRFP/2KVOmoLi4GG+99RbmzJmjVRuJiIiImk1QPUPHjh3DhAkTGtyXmZmJ/Pz8kBpFRERE1FKCCkOJiYk4evRog/uOHj2KxMTEkBpFRERE1FKCCkMTJ07EsmXLsGnTJrjdbgC+tYX27duHlStX4tprr9W0kURERETNJag5Q3PmzMGhQ4cwZ84c/x3rZ8yYAbvdjiFDhuAPf/iDpo0kIiIiai5BhSGz2YyXXnoJ3377LbZs2YLy8nJER0dj2LBhGDt2LFegJiIiojYj6HWGAGDUqFEYNWqUVm0hIiIianFBhaHly5df8Jh77703mKcmIiIialGah6GoqCgkJSUxDBEREVGbEFQYysrKqrfNbrdj165dWLhwIR555JGQG0ZERETUEkK6UWtdNpsNY8aMwT333IO///3vWj0tERERUbPSLAzVSk1NPeeCjEREREStTUhXk9UlhEBubi5efPFFdOzYUaunJSIiImpWQYWhtLS0c64lJITgMBkRERG1GUGFoXvuuafBMBQVFYVx48ahW7duobaLiIiIqEUEFYbuu+8+rdtBREREFBZBhaHTp0836fgOHToEcxoiIiKiZhdUGBo/fnyT7j924MCBYE5DRERE1OyCCkP/+Mc/sGDBAvTr1w+TJ09GcnIySktL8cUXX+Djjz/G73//e15RRkRERG1CUGFow4YNGD9+PBYtWhSwfeLEiUhMTMT333/P23EQERFRmxDUootbt27Fdddd1+C+yy+/HLt27QqpUUREREQtJagwFB8fjx9++KHBfd9++y2Sk5NDaRMRERFRiwlqmGzq1Kl47rnn4HA4MH78eCQkJKCoqAgbN27Em2++iUcffVTrdhIRERE1i6DC0N13343KykqsXbsWq1evBuBbeToiIgJz587FzTffrGkjiYiIiJpLUGFIkiQ8+OCDuPvuu/HDDz+gvLwc8fHxGDhwIKKiorRuIxEREVGzCelGrVFRUUhKSgIADBw4EF6vV5NGEREREbWUoMPQ+++/jyVLlqCwsBCSJOHtt9/Gs88+C5PJhCVLlsBsNmvZTiIiIqJmEdTVZBs3bsT8+fMxYsQIPPPMM1BVFQBw9dVX46uvvsLKlSs1bSQRERFRcwmqZ2jVqlW4+eabsXDhQiiK4t8+ZcoUFBcX46233sKcOXO0aiMRERFRswmqZ+jYsWOYMGFCg/syMzORn58fUqOIiIiIWkpQYSgxMRFHjx5tcN/Ro0eRmJgYUqOIiIiIWkpQYWjixIlYtmwZNm3aBLfbDcB3uf2+ffuwcuVKXHvttZo2koiIiKi5BDVnaM6cOTh06BDmzJkDWfblqRkzZsBut2PIkCH4wx/+oGkjiYiIiJpLUGHIbDbjpZdewrfffoutW7eirKwM0dHRGDZsGMaOHQtJkrRuJxEREVGzCCoM/e53v8Ott96KUaNGYdSoUVq3iYiIiKjFBDVnaMeOHTAYDFq3hYiIiKjFBRWGRo0ahbfffhsul0vr9hBRC5FlCRzRJiIKcpjMYrHg448/xubNm9GpU6d6l9JLkoRXXnlFkwYSUfNwexWYDEH9PUREdFEJKgzl5eVh0KBB/q+FEAH7z/6aiFofIQC3V2UgIiLda3QY+uCDDzBmzBjExcXh1Vdf1awBqqpi+fLlePvtt1FRUYFLL70UCxYsQNeuXRvVpj/96U/4/PPP0alTJ83aRKQXLo8Cs1EG/34hIj1r9J+E8+bNw4kTJwK2rVq1CkVFRSE1YOXKlXjzzTfxxBNPYN26dZAkCbNmzfIv5ngup06dwl/+8peQzk2kdx6PyiBERLrX6DB09tCXoij45z//GdJ9yNxuN9asWYP77rsPY8eORVpaGpYuXYr8/Hxs3rz5nI9TVRUPPPAA+vXrF/S5iQhQVAG3V7nwgUREF7Gg5gzVCnVuUFZWFqqrqzFixAj/tpiYGKSnp2PHjh2YNGlSg49btWoVPB4P7r33XmzdujWkNtQycN5EyGpryFqGpqXqqAKQJMCrCERGXJxDZXxPaoe11AbrqB0tr4YNKQyFKi8vDwCQmpoasD0pKQm5ubkNPmbv3r1Ys2YN3nnnnZB6pc4WExOh2XPpHWupjeauo8PlgSJJkCAhMioCJuPF+8OZ70ntsJbaYB1bl7CGIYfDAcB3e4+6LBYLysvL6x1vt9vxpz/9CX/605/QrVs3TcNQRYUDiqJq9nx6ZDDIiImJYC1D1FJ1dHtVVFQ4ISAgFC8izGH9cdAs+J7UDmupDdZRO7GxEf77o4Yq5J9+odyHzGq1AvDNHar9HABcLhciIuqn5ieeeALdunXDzTffHPQ5z0VRVHi9fGNqgbXURnPXUVFUKIqAKgSq7B6YDRfnUBnA96SWWEttsI6h0/LnVZPC0D333FOvF+d3v/sdTCZTwDZJkvDZZ59d8Plqh8cKCgrQpUsX//aCggKkpaXVO379+vUwm83+NY4UxTfx87rrrsPkyZPx2GOPNeXlEFENj1eFx6vCyHkMRKRDjQ5DN954o+YnT0tLQ1RUFLZt2+YPQxUVFdi/fz+mT59e7/hPP/004Os9e/bggQcewAsvvIAePXpo3j4ivVBVAadbQbTt4u0dIiI6l0aHoUWLFml+crPZjOnTp2Px4sVISEhAx44d8fTTTyMlJQUTJkyAoigoKSlBdHQ0rFZrvYUYaydgd+jQod4tQYioaZxuBVERRgC8YRkR6UvY+8Rnz56NqVOn4uGHH8a0adNgMBiwevVqmM1m5ObmYvTo0di4cWO4m0l00VNUFU4P5zAQkf5IgjcSAwCUllZzMluIjEYZ8fGRrGWIWqqOHkVFaYULap0fARazAQnR1ovm/oJ8T2qHtdQG66idhIRIzdZrCnvPEBG1Hl4vr3AhIv1hGCIiP0UVcLi9mq7sSkTU2jEMEVEAp1u5aIbJiIgag2GIiAJwIjUR6Q3DEBEFEAJwOD28wp6IdINhiIjqcXtVOFycO0RE+sAwRET1CAHYnd6Ay+6JiC5WDENE1CCvV4Xdyd4hIrr4MQwRUYMEAIdLYe8QEV30GIaI6JwUVYWLV5YR0UWOYYiIzkkIwOXyQuJYGRFdxBiGiOi83IoKReVQGRFdvIzhbgARtSxVCJzIr4TDpcDl8SIp3gb5PD0/qirg8ngRYeaPC2q9at/XVXYPomwmdEmOPu/7mqgu/nQj0pEDx0vw0dZs5JfYkRgbgfIqJ6JsZozL7IAeneIafIwQvlt02Cwm3qaDWqXa93VeiR2KImAwSEhJsGHSiK7o2y0h3M2jNoDDZEQ6ceB4CV755CByCqtgMRkQGWGC2WRAXokD731zDEdzys75WK9XhVfhRGpqfc5+X8dEmWExGZBTWI1XPjmIA8dLwt1EagMYhoh0QBUCH23NhtPtRVyUBWaTAbIEGI0GxNhMcHlUfLnn9Dkvo1dUwRWpqdVp+H0twWwyIC7KDKdbwUdbs7k8BF0QwxCRDpzIr0ReiR2RVlO9K8MkSYLNYkBRmQO5RdXnfA6nm2sOUetyofd1pNWIvBI7TuRXhqmF1FYwDBHpQJXdA0URMBob/l/eYJChqL5bcJyLoqhwcEVqakUu9L42GmUoikCV3dPCLaO2hmGISAeibCYYDBK83obn/SiKCoMM2KznvqaidkVqTqKm1uJC72uvV4XBICHKZmrhllFbwzBEpANdkqORkmBDtdNbL8wIIWB3KWgXF4HUdpHnfR6vqsLhVtg7RK3Chd7X1U4vUhJs6JIcHaYWUlvBMESkA7IkYdKIrrCaDSircsPtUaAKwOtVUGH3wGKSMS6zwwXXZRECsDu84IVl1Bo0/L4WcHsUlFW5YTUbMGlEV643RBfEMESkE327JWDmNX3QqX0kXB4F1Q4P3B4FKQkRuHF093OuM3Q2j6Ki0u5m7xC1Cme/ryuq3HB5FHRqH4mZ1/ThOkPUKFx0kUhH+nZLQJ+u8U1agbohLrcCh1uB1WRoppYSNV7d9zVXoKZgMAwR6YwsSeiWEgOPoqK0whXU5fKqEKh2eGA1yQD4C4fCr/Z9TRQMDpMRUVC8igqnWwl3M4iIQsYwRERBEQJclZqILgoMQ0QUNI9XwOXhpWVE1LYxDBFR0FQh4HB66t0KgYioLWEYIqKQuL3qOVcAJiJqCxiGiCgkiirgcHPuEBG1XQxDRBQy3tGeiNoyhiEiCpmiqHC4eM8yImqbGIaIKGS+O9p7wc4hImqLGIaISBNeRYXTw0UYiajtYRgiIk0IgZrL7MPdEiKipmEYIiLNeBQuwkhEbQ/DEBFpRlW5CCMRtT0MQ0SkKS7CSERtDcMQEWlKUQXsXISRiNoQhiEi0pzT7YVX4XX2RNQ2MAwRkeYURaDC7g53M4iIGoVhiIiahdujoJqX2hNRG2AMdwOIqGV5FRXr/3cUJZUujM3sgIQYa7OcRwjA7vTCbDLAZODfXUTUevEnFJHO/HC4CJ9sP4kdBwqw6v2fkF9qb7ZzKapAlcMDsHeIiFoxhiEinUlOsPmzSZXDg5c+PIC8kuYLRG6PAqfL22zPT0QUKoYhIp3pnBSF2yf19QeiaocHL324H7nF1c1yPiGAaqcXAry6jIhaJ4YhIh0alZGK2yb19U9utju9eOnD/cgprGqW83m9KqodXq5MTUStEsMQkU4NT0/GzeN7Qa7JJw6XgtUfHkB2XqXm5xIA7C4vXB4OlxFR68MwRKRjmT3b4eYre0Gu6bFxeRS8vPEAjp4u1/xcqipQWe2BytEyImplGIaIdK7/JYn49YReMNR0Ebm9Kl75OAsHT5Rqfi6PoqKi2sWry4ioVWEYIiL07ZaAW6/t418PyKsIvPrJIew9Wqz5uVweBZV2LsZIRK0HwxARAQB6dYrDbyamwWIyAABUIbDu88PYkVWg6XmEABxOL6qdnFBNRK0DwxAR+XVPjcFvr+uLCItvcXoB4L2vfsZXe05reh5V+BZjdPDu9kTUCjAMEVGATu2jcOcv0hFtM/m3bdp2Apu2ZUMI7WY/q6pApd0Nl0dlICKisGIYIqJ6khNsuGtyPyTEWPzbvtqTi3e/+hmKhpeDKYpARbUbHq+q2XMSETUVwxARNSghxoq7JvdDaqLNv23XwUK8/ukhuL2KZufxKirKq9xQFF5zT0ThwTBEROcUbTPjjuvS0S012r8t60Qp1nx0AHanR7PzeBQVZdUuTXudiIgai2GIiM4rwmLEbf9fX6R3i/dvO5Ffhef/8xNKK12ancfjVVFWxUBERC2PYYiILshklDHtqt4Ympbk31ZY5sSqDftwqki7G7x6vOwhIqKWxzBERI1ikCXcMKY7rhrSyb+t0uHBi//5SdPVqj0eFeXVLqgaXrlGRHQ+DENE1GiSJGH84E6Ycvkl/hu8ur0qXv3kILbtz9fsPG4Ph8yIqOUwDBFRkw1JS8Kt16bBbPL9CFEF8P43x/Dx1mzNenTcNT1EDERE1NzCHoZUVcWyZcswZswYZGZm4vbbb0d2dvY5jz98+DDuvPNODB8+HCNHjsTs2bNx+rS2q+MS0YX17hyHO3/RDzF1Fmf8em8u3th8CG6PNpfe+3uIeNk9ETWjsIehlStX4s0338QTTzyBdevWQZIkzJo1C263u96xpaWluO222xAZGYnXXnsNL774IkpLS3HHHXfA5dLuqhYiapwO7SLx+xv6B6xFtP94KV74YD/Kq+v/PxyM2qvMPApXqiai5hHWMOR2u7FmzRrcd999GDt2LNLS0rB06VLk5+dj8+bN9Y7/7LPP4HA48NRTT6FXr17o378/nn76aRw9ehTff/99GF4BEcVGWXDnL/qhT+c4/7bTRdVY+d6PyCmo0uQcnpqFGZ0ehYGIiDRnDOfJs7KyUF1djREjRvi3xcTEID09HTt27MCkSZMCjh85ciRWrFgBi8Vy9lOhvLw8pLYYDGHvJGvzamvIWoampeqoAjAYJEgajEDZIoz4zcQ0fLQlG9/szQUAVNo9eOGDn3DTFT0xsFe7kM8h4Lu5qyRJiLQa0ZipSXxPaoe11AbrqB0t/zAKaxjKy8sDAKSmpgZsT0pKQm5ubr3jO3XqhE6dOgVse/7552GxWDB06NCQ2hITExHS4+kM1lIbzV1Hh8sDRZIaFSoaa/rEdHRJjcGbmw9BVQW8isC/PzuMkio3rr+8B2Q59J9ekgQIgwExNnOjf6HwPakd1lIbrGPrEtYw5HA4AABmszlgu8ViaVRPz7/+9S+88cYb+POf/4zExMSQ2lJR4YCi8GaRoTAYZMTERLCWIWqpOrq9KioqnJqv55N5SQKiruuL1z45BLvLCwD4ZGs2sk+X4+areiHCEvqPnTIAxSYDoiPMsJjlcwY6vie1w1pqg3XUTmxsBGRZmx62sIYhq9UKwDd3qPZzAHC5XIiIOHdqFkLgn//8J5577jncdddd+M1vfhNyWxRFhZd3ztYEa6mN5q6joqhQFNEsixt2S4nB72/sj1c/OYiCUt8fPVknyvDs+h8xfUJvJCfYLvAMF+ZQvHC7FURFmBAZYYJ6nkvw+Z7UDmupDdYxdFr+6ArroGXt8FhBQUHA9oKCAqSkpDT4GI/HgwceeACrVq3CvHnzMHfu3GZvJxE1XWKMFb+/vn/APc2Ky514bsM+/PhzsSbnUFSBSrsHFdVugBOriShIYQ1DaWlpiIqKwrZt2/zbKioqsH//fgwZMqTBx8ybNw+bNm3CkiVL8Nvf/ralmkpEQbCYDbhlQm9ceemZuX5ur4p/f3YYH2/N1mRBRVUIVDs9qGQgIqIghXWYzGw2Y/r06Vi8eDESEhLQsWNHPP3000hJScGECROgKApKSkoQHR0Nq9WKd999Fxs3bsS8efMwbNgwFBYW+p+r9hgial1kScKVl3ZCx3aRWPfFEbhqFmT8em8ucgqrcPOVvRBtM1/gWc5PCMDu8kIAiIk0A1yjkYiaIOzX9s2ePRtTp07Fww8/jGnTpsFgMGD16tUwm83Izc3F6NGjsXHjRgDAhx9+CAD4+9//jtGjRwd81B5DRK1TWtd43DOlP5Ljz8wHPJZbieXrf8TPpytCfn4hAIfLi/IqNwTTEBE1gSQEbw0NAKWl1ZzMFiKjUUZ8fCRrGaKWqqNHUVFa0fJ3h3d7FLz39c/Yc+TMvCFJAiYM6YzLB3aArMHiIWaTjGibGTarEXFxfE9qgf9/a4N11E5CQqRm6zWFdZiMiPTHbDJg6rgeiI0045u9uVCFr1fn0x0ncSy3Ajdd0RM2qxG5RdWwO72wWY1IbRfZpJDk9qgoq3TBq6iIim6+9VxUIXAivxJVdg+ibCZ0SY7WJMy1Nl5VxbYf82F3qbBZZFzapz2MGl3STNQaMAwRUYs6mlOGL/ecRlGZAxazAU634r9E9nBOOf7x1h7ERpthd3igqIBBBtrFRWBcZgf06BTX6PMoqm/F6qIyB4RXgcmg7QKTB46X4KOt2cgrsUNRBAwGCSkJNkwa0RV9uyVod6Iw27QtGx9tyYajZk6WBOA1ixGTRnbFtcO7hrt5RJpgtCeiFnM0pwzvfXMMeSV2mE0GxEZZkBBtgaHOytR2lxe5RXZ4vCoiI4wwmwzIK3HgvW+O4WhOWZPPqagC5VUulFW5oQqhyRL+B46X4JVPDiKnsAoWkwExUWZYTAbkFFbjlU8O4sDxktBP0gps2paN9f/7GdVOL2RZgtEgQZYlVDu9WP+/n7FpW3a4m0ikCYYhImoRqhD4cs9puDwKYmxmmIwGyJIEi9mI9nFWmI2BP47sLgUlFS5IkoQYmwkuj4ov95wOao6TKgQcLi9KKlxwuLwhBSJVCHy0NRtOtxdxURaYTb7XYTYZEBdlhtOt4KOt2S0+F0trXlXFR1t8yx+YakKQLMmQZQkmgwRFFfhoSza8Kue9UNvHMERELSK3qBpFZQ7YLEZIZ6URWZYRafWN2te9fZnHq6KwzAGHy4sIs4yiMgdyi6qDboNXUVFR7UFZlRuKKuq1ozFO5Fcir8SOSKup3uNrbyKbV2LHifzKoNvZGmzfnw+HywujLDX4Oo2yBIfLi+3788PUQiLtMAwRUYuwO72+OUDnuvqj5vdtbKQp4P5lQgBlVW5UOrzwKgJ2pzekdvh7iSqdqKh2NTkUVdk9UBQBo7Hh12E0ylAUgSq7J6R2hltJhdO3QMG5SiP5lnMqqXC2XKOImgnDEBG1CJvVCIOMc9+cUtT+3pUQH21BfLQlYDjL5VZgd3lRUObQpD2KIlDtPBOKvI0MRVE2EwwG6ZyXRXu9KgwGCVE2kybtDJeEGKvv+3Gu0b6a71dCDBe7pbaPYYiIWkRqu0i0i4uA3aXg7OXNhBBwe1VYzQa4FRVCCERYjEiKi4DZJNc5DvhoSzbe++pnuNyKJu2qDUWljewp6pIcjZQEG6qd3gZfR7XTi5QEG7okR2vSvnAZlp6MCIsRXlU0+Dq9qu97NCw9OUwtJNIOwxARtQhZkjAuswMsJhkVdg88XgWqEPB4FVTYPbCYDRg3qAMsJoN/vyT7Jk+fPbl6R1YB/vnOHhzJKdesfWf3FClKw6FIliRMGtEVVrMBZVVuuD2+1+H2KCircsNqNmDSiK5tfr0hoyxj0siuMMgSPIqAqgqoQoWqCngUAYMsYdLIrlxviC4KfBcTUYvp0SkON47ujpSECLg9CqrsHrg9ClISInDj6O4Yk9mx3n6PV0Xn5ChMGdMdHdpF+p+rrMqNNRsP4L2vfobTHdo8orrqhqKyKhc8ilrv6rO+3RIw85o+6NQ+Ei6PgooqN1weBZ3aR2LmNX0umnWGrh3eFf839hJEWo1QVQFvTSiKtBrxf2Mv4TpDdNHg7ThqcGn00HGZeW1c7LfjAHyTmM+3wvS59iuqii93n8Z/vz8V0O7YSDNuGNMdfbrEB5zHYJAQG2tDeblvYcRgyLIEk1GG1WyAxWSAQZb9w0Z6WoF618FCrkCtAf6c1I6Wt+NgGKrBN2bo+D+5NvQQhkJ1uqga6/93FLnF9oDtA3ok4rrLuiEqwjd5WYswVEuSfMHIYjQgwmqE2ShruqJ1a8f/v7XBOmpHyzDEaE9EbU6HdpG4+8b+uGpIp4DVq/ceLcbSt/Zg18GCepN+QyWEbwjN7vKitNKFkgonnB7fJO6LsDOISFcYhoioTTLIMsYP7oR7p2Sgc1KUf7vD5btVxIsf7kd+if08zxA8VRVweVSUV7lQXO5EtcNTc6sPpiKitohhiIjatOQEG+6a3A+/uKxbwGX4x3Mr8Y+39+K9L4/A7dHmMvyzCeEbbqywe1Bc4UT5OSZcE1HrxjBERG2eLEsY2T8F99+UifRuZyZRq6rAJ1uzsfjNH/Djz8WaD53VVXcIrbSyNhQxFRG1BQxDRHTRiI2yYPrVfTDj6t6IizL7t5dXufHvzw7j5Y1ZyC9tnqGzWqoq4HQrKK10oazKBbeXPUVErR3DEBFddPp2S8CcX2biikEdAyZYHzlVjmff2YsPvjsOh0u7tYkaoqq+e6CVVdVMtnZ7AQgGI6JWiGGIiC5KZqMB147ogkd+Oxy9Osf6t6sC2LIvD4vf/AFb9uVBUZv38mb/ZOtqN4rLXai0e+CpuT8bh9GIWgeGISK6qKUkRuK3k/pi+tW9kRBt8W93uLz44LvjWPbOXhw4XtKs84mAM5OtqxwelFa6UFTuW+Ha4fZe8H5oRNS8jOFuABFRc5MkCendEtCrUxy+/TEXX/5wCm6Pr3emsMyJVz89hO6p0bh2eNeAy/Sbi6oKqBDwKiqc7ppVrg2+Va5NRgOMBl8w0tOijkThxDBERLphMsoYN6gjLu3THpt35mBXVgFq88ax3Eo8t2Ef+ndPwIShndE+LqJF2lS7mKOiKHB5FEiSBIMswWyUYTIZYDHJkCW52XuuiPSMYYiIdCfaZsaUyy/ByH7J2LTtBA7nlPv37TtWgv3HSzC4TxLGD+6IuCjLeZ5JW0IAQvhuhurxqpBcXsiyBKvJAKtFf7cAIWopDENEpFupiZG4bWJfHM4pwyfbTuB0zb3OVAHszCrAD4cLMbRvMsYN7IBom/kCz6a92l6jasULh1uBySDBajbCbOJQGpGWGIaISPd6dYpDj46x+PFoMTbvPImSChcAwKsIbNmXh50HCjA8PRljMlPDEoqAmqvSVAG3xw1ZlmAw+HqMLGYjTP4eI8FwRBQEhiEiIgCyJCGzZzv0vyQBO7MK8d/dp1BR7Qbguwrsmx9zsW1/Pob1TcKYzA6IiQxPKBIAFFVAUQXcHhUGp28ozSBLMBhkmIyyf54RwxFR4zAMERHVYZBlDE9PxuDe7bH9QD7+98NpVDk8AHyh6Nt9edh2IB+X9knCmAGpSIixhrW9tcHI10IFkuS7Os0oSzCZDDAZZRglGbIBDEhE58AwRETUAJNRxqiMVAztm4Tt+wvw1Z4zocirCGzbn48dB/IxoEc7XD6wA1ISbGFusc+Zq9N8iz1Kkm9pAVnyDa2ZjTLMNSFJAuccEQEMQ0RE52U2GjB6QCqGpSdhx4ECfL3nNCrsvlCkCuCHI0X44UgR+nSJw5gBHdA9NbpVLaDov0INAl4FcLkVyLLXd/m+yQCTQYLRKMNkkP3HE+kNwxARUSOYjQaMykjFsL7J+P5QIb7ecxollS7//oMnynDwRBk6tovEqAGpyLgkAQa5dS7yr6p1Lt+v6TkyyBIsJgPMJhlGg+9DCA6pkT4wDBERNYHJ6JtTNCQtCT/+XIyv95xGbs0l+QBwqqgab31xBJu2mTEiPRlD+yYh0moKY4vPr97aRk74J2QbDb4eI4NBhiT7jvOFJ/Yg0cWFYYiIKAgGWcLAnu2Q2SMRh3PK8c3eXBw5dWbxxopqNz7dcRJffJ+DzJ7tMKJfCjq2iwxjixun7pyj2luWyJIEo1GCFzIcdhdkAAajDKMsQ5YBCRJkWfKvks2gRG0NwxARUQgkSULvznHo3TkOp4uq8e2Pudh7tBiK6ksEXkVg18FC7DpYiM5JURiRnoz+lyTCZGydQ2gNUYWAovquXHO4vFAU4R9ekyRfGJJkwFjTm2QwyL5L/WVfSJJ8U7UZkqjVYhgiItJIh3aRuOmKnrh2eBds25+P7QcK/FegAcDJgiqcLKjCh1uyMbh3Owztm4ykFroHmtZqh9dqvgIUBFzeL0HyByZjzSRtY81aSLVBCTgz0dw39MbAROHBMEREpLFomxlXDemMcYM64sefi7H1p3ycLKjy73e4vPj2xzx8+2MeuqZEY0if9si4JBFmkyGMrdaOEICA8K0QWXMVG9xKQG+SLPnWQoLkm9At4JukbjKdWRfJIJ9ZWbv2eYmaA8MQEVEzMRpkDOrVHoN6tcepomps+ykPe44Ww+NV/cdk51UiO68SH3x3HAMuScTgPu3RNbl1XZ6vlbq9SQpqF4o8w+3xTeA+OzDJBhlyzdeyoWb4rWbtJEi1E8ADn6t2ojeE7/lQ8/jadtS0iAGLADAMERG1iI7tIjFlbA9MHNkVuw8XYWdWQcBVaG6Pip0HC7HzYCESYiwY1Ks9BvZqh8Qwr3Dd0hoOTIp//9lzlQBfL9TZoUby/+fMcWcej5pFKGX/vCaj4UzAkiTfPKmzA1bt4+vOg6ptM7VtDENERC3IajZiZL8UjEhPxqmiauzMKsCeI8Vwec78wi+pcOHzXTn4fFcOOidFIbNnO2RckhC2m8S2JvXmKjXuUefY7qt5YwOW/1hIkGX47wdXOwfKIPmWIJCl2qvrzjymtu0X6vCreyy1HIYhIqIwkCQJndpHoVP7KEwc2RX7j5Xi+0OFOHqqPOBXd+2k64+2HEePDrEY0CMR6d3iYWvFaxe1NU0PWKJOZ5Xi74UKuLrO34MEKDUTw2tX/RYGA+wuL1RV+I9XhfCv96QKIDqC39+WxDBERBRmZqMBA3u1w8Be7VBW5cKeI0XYfbgIBaUO/zFCAEdOlePIqXJs+FpCj44x6H9JIvp2jUcUf3GGVe0c7wsFKi8Ar6LCZDGhvMoFVT3ruJo550aDjKgIY8DwHjUvhiEinZIkwGKqWeumzl+ztTsl/3yLMxdAq6j5K1qt+SsWvrkVomZuhX9oQZx7mIHOLy7KgrEDO+LyzA7ILbZjz5Ei7D1ajPJqt/8YVQgczinH4ZxybJCAbikx6Nc9HundEhAXZQlj66kp+P9H68EwRKRTJoOMuOjaybnn/6l89tyHmq/824QQNd38gFABFQKKGtjtr6oqIHyBCsK3rc5V0/7wxLVmfCRJQod2kejQLhLXDO+CE/mV2Hu0GD/9XILKOmsXCQEcy63AsdwKfPhdNjok2pDWNR59uyWgQ6LtorwqjUhrDENEOlV3/ZamPcb/VcA2f8+SATBAQu2SOQ0FqLrba8NPbQ+TIlQoioBXFVAVFYqoCVZqzbFNavHFQZYkdEuJQbeUGFw3shuy8yux7+cS/HSsGBX2wAvUTxfbcbrYji++P4UYmwm9u8SjT+c49OwYC4v54ljHiEhrDENE1KwaClD1t9esAWMADDAAxtqw5AtYiqrCqwgoiurrZRICqiJ8QUkRAT1TFztZltA9NQbdU2Mw6bKuyCmowv7jJfjpeCmKy50Bx1bYPdiZVYCdWQUwyBK6pkSjd6c49Ooci+QEm3/dHSK9YxgiolaptudKCF+vk8kgwWTwzXE68ztcgoAvENUGJo9Xhduj+O8NdjGTJQldkqPRJTka1wzrgsJyJw4cL0FWdhlOFFQGhENFFfj5dAV+Pl2BTduByAgTenaMQc+OsejZMRaxnGtEOsYwRERtTt0VhIGa9V5kA8xGQLJK8Coq7C4vvIp6wXVdLhaSJCEpLgJJAzti7MCOqHZ6cOhEGQ6eLMPhnDI4XErA8dUOD/YcKcaeI8UAgMRYKy5JjcElHXwfXNOI9IRhiIguKkIIGGQJMTYTDEYZNpsFwqug2u6GRxH1L2e+SEVaTRjUuz0G9W4PRRU4VViFgyfLcCSnHDmFVfWGFIvLnSgud2JHVgEAoF2s1T8c1zUlGomx7DmiixfDEBFdlGov8beYjbBZjLCaDPB4VFQ7PXDpZBitlkE+M5w2YUhn2J1eHD1djqOnynEkpxwlla56jykqd6KoTjiKjTSjV5c4dEi0oXP7KKQk2mCQ5ZZ+KUTNgmGIiHRBVQUMBgmxURZ4vAoqHR64PYouJl2fzWY1IuOSRGRckggAKK104efT5f45RXXXNKpVXu3GzgMF/q9NBhkdkyLRJSkKnZKi0TkpCrGRHFqjtolhiIh0RQgBo0FGfLQF1Q4P7C4vFEWHiaiO+GgLLu2ThEv7JEEIgdJKV83aRZXIzqtEcYWz3mM8iorjuZU4nlsJIBcAEG0zoVP7KHRsH4mONWskce4RtQUMQ0SkTwKIijDBajbC4fLC6Vb86xjpZV5RQyRJQkKMFQkxVlzaJwkAUGF342RBFfJKHTiUXYrTRdUNDjNW2j04kF2KA9ml/m0xkWZ0SLQhtV0kOiRGIjXRhvhoCxeDpFaFYYiIdKv25pnRNhOiIky+9YuEgNutwOFW/Kto1x6rVzE2Mwb0SMSYWBvKy+1wuhScLqquuYlsJU4WVKGsqv7QGgBUVLtRUe1G1oky/zaLyYCUBBtSEm1ISbAhOSECyfE2RFj4K4nCg+88ItK92qAjSxJkSYIpQobNaoJXVaGqvluJuNwK3IrKW4YAMBlldE2JRteUaACpAIAqhwc5hVU4VViNnMIqnC6sDrhtSF0uj4Ls/Epk51cGbI+JNCM5PsK3REB8BNrXfG6z8ka01LwYhoiIzlIbdIyyDMgAICPCYoIqVCiqgKoCHq/C+UZ1REWYkNYlHmld4v3bKuxunC6qxumiauQW2ZFbXN3glWv+42t6kQ7nlAdsj7Qa0S4uAu3jItAu1or2sVYkxkYgIcYCo4FXtFHoGIaIiBpBCAEJEoyyBMiAxSQjwmJEtcMDp0dhKGpAjM2MmC7mgIDkdHuRX+JAbkk18ortyC9xIL/UDqdbOefzVDu9qM7zTeauS5KAuCgLEmOsSIixIDHWisQYK+KjLUiItvJebNRoDENEREEQwjesFhNpQaSqwuHy+tcvUlV9D6Odj9VsrDPE5iOEQEW1G/mlDhSUOlBQakdBme/z84UkIXzLApRWuoBT9ffbrEYkRFsQX/MRF21BfJQFsVG+fxmWqBbDEBFRCIQQkCUJUREmRNvMvvujKQo8HhUuj+K/gaze5xmdjyT51n+KjbKgd+c4/3YhBKocHhSWOVBY5kRRuQNFZb7FIEsrnbjQRX92pxd2pxc5hdUN7reaDYiLsiA2yozYSDPioiyIifR9Xvuv2cTApAcMQ0REGqgNPLIEWIwGWE0GRMN3hZpXEVCUmhvJKqr/KjWGpPOTJAnRNjOibWZc0iE2YJ+iqiitcKG4wlnz4UJJhRMlFU6UVrrgbcSwpdOtIK/EjrwS+zmPsZgMiIk0I9pmQozN92+0/18TomxmREeYYDUbuFxAG8YwRETUDGoDjgQJJoMEk0GuuWmsBMB3Cb+iCni9KjxeFS6P6ru0X8drHDWFQZbRLi4C7eIi6u1ThUCl3YPSSidKK1worXL5h9PKKl0or3Y3+nYsLo9S0zPluEB7fL2DtR+RESZERRgRafV9Hmn1fR4daYIlwuxfsoFaB4YhIqIW4vv95/slWDsZ22iWEWHxbfUqKjweFW6P4ruprBDsOQqCLEmIrRnm6pZSf79aM/xWXuVCWZUb5VVulFf5QlJ5zRVtlXb3BYfh6lJU4X98YxhkCbaagBRhMcJm9d1Dz2Y1ItpmxpgBKWgfZ2t8AygkDENERGFW91J+o8W3xpEQAooi4FFVX++RokJRBASE/ya0/s+pSWRJ8l3pZjOjc1LDx6iqQJXTg8qacFRh96DC7kal3YPKmn+rHB5U2T1Qg/gmKKqoea6G12L6as9p/O13I2HhnKUWwTBERNTK1A6hGAwSJFlGYakddqcXkTYTUhMjARVQalfL9irY8mMuisqdiI4wYWh6Mgyy7O9RUoVAblE17E4vbFYjUttFQq4zt+VC+wHf/Jw9R0vgdKuwmmX0757QpDvWa9GGUM+hqCr2HilCWZUbcVFmDOjZ7ryvQZbPBKaO7c99DgBwuryorAlLpwp9q3ErqoAk+SZxVzk8sDu9qHZ64HCd++q4uqocHhzPrUCvznFNrgU1nSTCPHCpqiqWL1+Ot99+GxUVFbj00kuxYMECdO3atcHjS0tL8cQTT+Crr74CAFx77bX485//DJsttO7E0tJqeL1qSM+hd0ajjPj4SNYyRKyjdtp6LQ8cL8FHW7ORV2KHoggYDBJSEmyYNKIr+nZLwGufZmHLT/mQJAlGgwyjQYbJKGFY3yRMHNkNR0+V49u9ecgvs8PjUSGEQGyUGaP6paB7h1gcySnDl3tOo6jMAUUFDDLQLi4C4zI7oEenOADA13tO4cvdp333boNvxpPVbMC4QR0wJrPjBV/D0Quc40L7G+NCzxHqa9D0dUiA0WxCbkEFquy+kGR3eWF3enC6qBrH8yqhKCqibWZ4FDXg+02BEhIiYdBo0c2wh6Hly5fjjTfewKJFi5CcnIynn34aJ0+exIcffgizuf7djmfMmAGXy4UFCxagoqICDz30EIYOHYq//e1vIbWjrf6wbE3a+i+e1oJ11E5bruWB4yV45ZODcLq9iLSaYDTK8HpVVDu9sJoNSEmIwN6jJed8/KBeiSgqd8GjKIiJtMBskgEVcHlVmI0yBvduh12HiuBye2E2GSDLErxeFXaXAqNBwrVDOyO3pBqf78yBr3S++UuK6rs6TgJwzbDO5w0TR3PK8N43x+DyKLBZjDAYZCiK7xwWk4zhfZOw7UDBOfffOLr7BQPRhc7Ro0MMdh0shCoAg4Ta+etQBCBLF34NWr8Og0FCbM093uou1Fn3HCkJNtisRhSVOf3f75nX9GEgOouWYSis65i73W6sWbMG9913H8aOHYu0tDQsXboU+fn52Lx5c73jd+/eje3bt2PRokXo168fRo4cicceewzvv/8+8vPzw/AKiIi0pwqBj7Zmw+n2Ii7K4gsrkgSzyYC4KDOqHe6AICTV+ai1+3AxSirsMBsNcLkVVFZ7UOnwwOX24nRRFT787jiKyuwwGmT/MIzJKCM20gRVFdh9uAA/HC5GVKQF7WMtaB8fgeSESKS2i0Tn9pFITozEgRNlMJlkWMwGmIyyv3fKN7wHbD9YAINBQlJcBKIizLCYDIiwGBEfZYbHq+J/P5yGy6MgxmaGyeh7jSajATE2E1weFV/uOX3e+TiqEPhyz7mfw+nyYmdNEDLKvqEvWZIgyxKMMqAK4Mvdp6Go5w7KFzqHy63gy93avg6jIfD77XQr+GhrdlBzk6hxwtoztHfvXtx0003YtGkTunfv7t8+bdo09OnTBwsXLgw4/sUXX8Qrr7yCb775xr/N7XYjMzMTS5YswcSJE4Nui6qqnIgYIkkCZFlmLUPEOmqnrdbSq6gor3b7Ak4D80Vq1ym6EEmSYJDrP762h6c2HDS0XxVn5i6dOUTC2UdH2UwNTvL1KiqqHB5IUv3HAPAvIyDLUoOvsXbOU2SEyXcLlNrtdY5RVBXVNedoiKr6Xod/RYOzn6Dmy0ir8ZyLKyr+11H3SQLbecHXAd+922rvoybLUsASCv5awfc9kyTfv7XH1D5HbKSZ92Kr41w1D0ZYJ1Dn5eUBAFJTUwO2JyUlITc3t97x+fn59Y41m82Ii4tr8PimkJswGZDOj7XUBuuonbZWy9oFA2t/MYaiocefWQOp4f2ABIg6vSXizCfirE2qKhocqvAqNYFLOncbVAFIApD8M3nq/iNBwBcyTGcFldqQJrw1R54jcEEGaq+P9//SDDjQF7gkSYLReNZrEGfO5Q+NDb0OVQKgwiDLkGTUC1uABFUIGGQpIMjIhjoBr2b/2cG0bgj1vWZJs2EhChTWMORw+BaxOntukMViQXl5eYPHNzSPyGKxwOU6952QiYjaEovZgCRz215jpjlfQ22wsZgMsDSw6KI2J/H9YzEbYTE376/KljgHnV9YI6bVagXgG+qqy+VyISKi/hvcarXWO7b2+FCvJiMiIiJ9CmsYqh3yKigoCNheUFCAlJT6y4ampKTUO9btdqOsrAzJycnN11AiIiK6aIU1DKWlpSEqKgrbtm3zb6uoqMD+/fsxZMiQescPHToUeXl5yM7O9m+rfezgwYObv8FERER00QnrIKXZbMb06dOxePFiJCQkoGPHjnj66aeRkpKCCRMmQFEUlJSUIDo6GlarFZmZmRg8eDDuv/9+LFy4EHa7HQsWLMANN9zAniEiIiIKStgXXVQUBc888wzeffddOJ1ODB06FI8++ig6deqEnJwcXHnllVi0aBGmTJkCACguLsZf/vIXfP3117BYLP4VqC0WSzhfBhEREbVRYQ9DREREROHEBQuIiIhI1xiGiIiISNcYhoiIiEjXGIaIiIhI1xiGiIiISNcYhoiIiEjXdBOGiouL8cADD2DEiBEYNGgQ7rzzThw5csS//8CBA5g+fToGDhyIcePGYfXq1WFsbdtw7NgxDBo0CO+++65/G+vYNKdOnUKfPn3qfbz99tsAWM+m2LBhAyZOnIiMjAxMmjQJH3/8sX8f69g427Zta/D92KdPH1x55ZUAWMvG8ng8WLp0KcaNG4dBgwbhlltuwffff+/fzzo2XnV1NR5//HGMHTsWl156Ke6++26cOHHCv1+TWgqduOmmm8SvfvUrsXfvXnHkyBFx3333iVGjRgm73S5KSkrE8OHDxUMPPSSOHDki3nnnHZGRkSHeeeedcDe71XK73WLKlCmid+/eYv369UIIwToG4fPPPxcZGRkiPz9fFBQU+D8cDgfr2QQbNmwQffv2FWvXrhXHjx8Xy5cvF2lpaeL7779nHZvA5XIFvA8LCgrEN998I9LT08Vbb73FWjbBP//5TzFq1Cjx9ddfi+PHj4uHHnpIDB48WOTl5bGOTXTHHXeIMWPGiC+++EIcOXJEPPzww+Kyyy4TJSUlmtVSF2GopKRE3H///eLQoUP+bQcOHBC9e/cWe/bsEatWrRJjxowRHo/Hv3/JkiXimmuuCUdz24QlS5aIGTNmBIQh1rHpnnvuOTF58uQG97GejaOqqrjiiivEU089FbD99ttvF6tWrWIdQ+B2u8WkSZPEnDlzhBB8TzbF5MmTxaJFi/xfV1ZWit69e4tNmzaxjk1Q+7v6yy+/9G9TFEVcffXVYvny5ZrVUhfDZPHx8XjmmWfQq1cvAEBRURFWr16NlJQU9OzZEzt37sTQoUNhNJ65VduIESNw7NgxFBcXh6vZrdaOHTuwbt06/O1vfwvYzjo23cGDB9GzZ88G97GejfPzzz/j1KlT+MUvfhGwffXq1bjrrrtYxxC8/vrryM3NxZ///GcAfE82RVxcHP773/8iJycHiqJg3bp1MJvN6Nu3L+vYBMeOHQOAgJu3y7KMtLQ07NixQ7Na6iIM1fXII49g1KhR2LRpE/7617/CZrMhLy8PKSkpAcclJSUBAE6fPh2OZrZaFRUVmDdvHh5++GGkpqYG7GMdm+7QoUMoLi7GLbfcgssuuwzTpk3D119/DYD1bKzjx48DAOx2O377299i5MiRuOmmm/DFF18AYB2D5XK5sGrVKsycOdNfL9ay8R566CEYjUZceeWVyMjIwNKlS/GPf/wDXbp0YR2boH379gB87726Tp06heLiYs1qqbswNHPmTKxfvx6TJ0/GPffcg59++glOpxNmsznguNobv7pcrnA0s9VauHAhBg4cWO+vcACsYxO53W4cP34cVVVVmDNnDl544QVkZGRg1qxZ2LJlC+vZSFVVVQCA+fPn47rrrsOaNWswatQo3H333axjCN5//324XC7MmDHDv421bLyjR48iJiYGK1aswLp16zBlyhTMnz8fWVlZrGMTZGZmokePHliwYAFyc3Phdruxdu1aHDhwAG63W7NaGi98yMWldkji8ccfxw8//IDXXnsNVqsVbrc74LjaItpsthZvY2u1YcMG7Ny5Ex988EGD+1nHpjGbzdixYweMRqP/f+b+/fvj6NGjWL16NevZSCaTCQDw29/+FjfeeCMAoG/fvti/fz9efvll1jFIGzZswNVXX434+Hj/NtaycU6dOoUHHngAa9eu9Q/vZGRk4MiRI3j22WdZxyYwmUxYsWIFHnzwQYwbNw5GoxHjxo3D1KlTsW/fPrjdbk1qqYueoeLiYnz44YdQFMW/TZZl9OjRAwUFBUhJSUFBQUHAY2q/Tk5ObtG2tmbr169HcXGx/1LRQYMGAQAWLFiASZMmsY5BsNls9f6q6d27N/Lz81nPRqrtIu/du3fA9p49eyInJ4d1DEJJSQl2796NiRMnBmxnLRtn79698Hg8yMjICNiemZmJ48ePs45N1L17d6xbtw7bt2/Hli1bsGLFCpSVlaFbt26a1VIXYaigoAB//OMfsX37dv82j8eD/fv3o0ePHhg6dCh27doVEJa2bNmC7t27IzExMRxNbpUWL16MjRs3YsOGDf4PAJg9ezZeeOEF1rGJsrKyMGjQIOzcuTNg+759+9CzZ0/Ws5HS09MRGRmJPXv2BGw/dOgQunTpwjoG4fvvv4ckSRg2bFjAdtaycWrnUx48eDBg+6FDh9C1a1fWsQmqqqowffp07Nu3D7GxsYiJiUFlZSW+++47jBkzRrtaanb9Wyumqqq4/fbbxTXXXCN27NghDh48KO6//34xdOhQcerUKVFUVCSGDh0q5s+fLw4fPizWr18vMjIyxLvvvhvuprd6dS+tZx2bRlEUcdNNN4nrrrtO7NixQxw5ckQ8+eSTon///iIrK4v1bIIVK1aIQYMGiQ8++EBkZ2eLlStXirS0NLF161bWMQjPPvusuPrqq+ttZy0bR1EUccstt4hrr71WbNmyRRw7dkwsXbpU9O3bV+zevZt1bKLp06eLadOmiaysLHHgwAFxyy23iMmTJwuPx6NZLXURhoQQoqKiQixYsECMGjVKDBgwQNx+++0B6w7t2bNH/PKXvxT9+/cXV1xxhXj11VfD2Nq2o24YEoJ1bKri4mLx5z//WYwaNUpkZGSIX/3qV2LHjh3+/axn461Zs0aMHz9e9OvXT0yePFls3rzZv491bJoFCxaIX/7ylw3uYy0bp6ysTCxcuFCMGzdODBo0SPzqV78S27Zt8+9nHRsvPz9f3HfffWLIkCFi2LBhYv78+aK4uNi/X4taSkII0Qw9W0RERERtgi7mDBERERGdC8MQERER6RrDEBEREekawxARERHpGsMQERER6RrDEBEREekawxARERHpGsMQEVEjcVk2oouT7u5aT0RN9+CDD+K999477zEdO3bEF1980UItanmff/45PvnkE/z9738Pd1OISGNcgZqILujEiRMoKSnxf71y5Urs378fy5cv928zm81IT08PR/NaxIwZMwAAr776aphbQkRaY88QEV1Qly5d0KVLF//XCQkJMJvNGDhwYPgaRUSkEc4ZIiJNHDp0CHfddRcGDx6MwYMH45577sHJkyf9+7dt24Y+ffpgy5YtmDFjBgYMGIBx48bh7bffRkFBAe69914MGjQIY8eOxdq1a+s97ptvvsGvf/1rDBgwABMmTMBrr70WcH5VVfHCCy9gwoQJ6N+/P6655pp6vTgzZszAn/70J8yePRuDBw/GnXfeCQDIycnBvHnzMHr0aPTr1w8jR47EvHnzUFpa6n/c9u3bsX37dvTp0wfbtm3zt2vbtm31zlHbiwQA48ePx5NPPomZM2di8ODBePTRRwEAZWVlePTRR3HZZZchIyMDv/zlL7Fly5bQvxFE1GQMQ0QUsmPHjuHmm29GcXExnnrqKfz1r3/FyZMnMW3aNBQXFwccO3fuXIwfPx6rVq1Ct27dsGDBAtx6663o3bs3li1bhn79+mHRokXYu3dvwOPuv/9+pKenY8WKFRg1ahQef/zxgLCzcOFCLFu2DJMnT8aqVatw7bXX4sknn8SKFSsCnufjjz+GyWTCihUrcOutt8LhcODWW2/F0aNHsWDBAqxevRrTp0/Hhx9+iGeeeQYAsGDBAqSnpyM9PR3r1q1Dv379mlSf119/HX369MGzzz6L66+/Hi6XCzNnzsTnn3+O+++/H8uXL0dKSgruuOMOBiKiMOAwGRGFbPny5bBarVi7di2ioqIAACNHjsRVV12Fl156CfPnz/cf+3//93+47bbbAAA2mw2/+tWvMGDAAMyePRsA0L9/f3z++ef4/vvvMWDAAP/jrrrqKjz00EMAgDFjxqCgoADPPfccfv3rXyM7OxtvvfUW5s6d6+/tGT16NCRJwvPPP49bbrkF8fHxAABZlvH444/DZrMBAA4cOICUlBQ89dRT/qHAESNG4Mcff8T27dsBAD179vS/rmCGBpOSkvDggw9Cln1/f7711lvIysrCW2+9hczMTADA5ZdfjhkzZmDx4sVYv359k89BRMFjzxARhWzr1q0YPnw4rFYrvF4vvF4voqKiMGTIEHz33XcBxw4aNMj/ebt27QDAHwgA+ENLZWVlwOOuv/76gK+vvvpqFBcX49ixY9i6dSuEEBg/frz//F6vF+PHj4fL5cKuXbv8j+vUqZM/CAFA37598cYbb6BTp044efIkvv76a6xZswY///wzPB5PiJXx6dGjhz8IAcCWLVvQvn179OvXz99WRVFwxRVXYN++fSgvL9fkvETUOOwZIqKQlZWVYePGjdi4cWO9fQkJCQFf1/aw1BUREXHBcyQlJQV8nZiYCACoqKhAWVkZAGDSpEkNPjY/P9//eW0Aq+vll1/G888/j9LSUrRr1w79+vVDREREvUAWrLPPWVZWhsLCwnMOtxUWFiI2NlaTcxPRhTEMEVHIoqOjcdlll/mHv+oyGrX5MVMbeGrVzkVKTExETEwMAOCVV15BZGRkvcd26NDhnM/7wQcf4KmnnsIf//hHTJ061R/e/vCHP+DHH3885+MkSQLgm7hdV3V1dYNtqCs6OhrdunXD4sWLG9zfqVOn8z6eiLTFYTIiCtmwYcNw5MgR9O3bFxkZGcjIyED//v2xdu1abN68WZNznL2g46ZNm9CxY0d06dIFQ4cOBQCUlpb6z5+RkYGysjL84x//qBek6tq1axeio6Nx5513+oNQdXU1du3aFRB06g5zAWd6uHJzc/3bysvLcfTo0Qu+lmHDhiE3NxeJiYkB7d2yZQteeuklGAyGCz4HEWmHPUNEFLK7774bN998M+666y5MmzYNFosF69atw2effYZly5Zpco61a9fCarVi4MCB+PTTT/Hf//4XS5YsAQD07t0bkydPxiOPPIJTp06hf//+OHbsGJYuXYpOnTqhW7du53zeAQMG4N///jeeeuopXHHFFSgoKMDq1atRVFQUMFQVExOD3bt3Y8uWLUhPT0efPn2QmpqK5cuXIzo6GrIs44UXXmjUkN+UKVPw2muv4bbbbsPvfvc7pKam4rvvvsOLL76I6dOnw2QyhVwvImo8hiEiCllaWhpef/11LF26FPPmzYMQAr1798aKFStw5ZVXanKO//f//h/ee+89PP/887jkkkuwbNkyXHPNNf79ixYtwvPPP48333wTeXl5SExMxMSJEzFnzpzz9rTceOONyMnJwfr16/HGG28gOTkZY8eOxS233IJHHnkER44cQc+ePfHrX/8a+/btw6xZs7Bo0SL84he/wLJly/Dkk09i7ty5aNeuHWbOnImff/4Zx44dO+9rsdlseP3117FkyRI8/fTTqKysRMeOHfHHP/4Rt99+uyb1IqLG4+04iKhV27ZtG2699Vb861//wvDhw8PdHCK6CHHOEBEREekawxARERHpGofJiIiISNfYM0RERES6xjBEREREusYwRERERLrGMERERES6xjBEREREusYwRERERLrGMERERES6xjBEREREuvb/A03Thu9NOyiGAAAAAElFTkSuQmCC\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 (ipykernel)",
+ "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.16"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}