"
+ ]
+ },
+ "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": 24,
+ "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:
Tue, 15 Nov 2022
Deviance:
3.0144
\n",
+ "
\n",
+ "
\n",
+ "
Time:
14:16:51
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: Tue, 15 Nov 2022 Deviance: 3.0144\n",
+ "Time: 14:16:51 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": 24,
+ "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": 25,
+ "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:
Tue, 15 Nov 2022
Deviance:
18.086
\n",
+ "
\n",
+ "
\n",
+ "
Time:
14:16:54
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: Tue, 15 Nov 2022 Deviance: 18.086\n",
+ "Time: 14:16:54 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": 25,
+ "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": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz2UlEQVR4nO3deXRV5b3/8c/JnBBCmDKAgSBSIAoBQpNfaFFbAkG5iO1tpYBMIveqsArkqhArhEgxjhStCFcU9IIottfSKhhMUwMoqRFoqBQKAkGQZsBSCBCSHHL27w9vjhwynZPpgfB+rXXWYj/72Xs/55sN+bBHm2VZlgAAAAzxMj0AAABwfSOMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKM8DiPbt2/X2LFj1a1bN9lsNm3atKnBZXJycjRkyBD5+/vrpptu0uuvv96IoQIAgLbI4zBy4cIFxcbGasWKFW71Lygo0JgxY/SDH/xA+fn5mjt3ru6//35t3brV48ECAIC2x9aUF+XZbDb97ne/0913311nn/nz52vz5s3at2+fs+1nP/uZzpw5o8zMzMZuGgAAtBE+Lb2B3NxcJSUlubQlJydr7ty5dS5TUVGhiooK57TD4dDp06fVuXNn2Wy2lhoqAABoRpZl6dy5c+rWrZu8vOo+GdPiYaSoqEjh4eEubeHh4SotLdXFixcVGBhYY5mMjAylp6e39NAAAEArOHHihG644YY657d4GGmM1NRUpaSkOKfPnj2rHj16qKCgQO3bt2+WbViWpXNlFdq2fbtuu/VW+fhelaW4alyyX6JWHqBe7qNW7qNW7qNW7quuVfKI2+Xn59es6z537px69erV4O/uFv8JRUREqLi42KWtuLhYISEhtR4VkSR/f3/5+/vXaO/UqZNCQkKabWwd7HZ1bB+kGyLD5Ovr22zrbYvs1Moj1Mt91Mp91Mp91Mp91bXq0qVLs9eqen0NXWLR4s8ZSUxMVHZ2tktbVlaWEhMTW3rTAADgGuBxGDl//rzy8/OVn58v6Ztbd/Pz83X8+HFJ35ximTJlirP/Aw88oKNHj+rRRx/V3//+d7388st65513NG/evOb5BgAA4JrmcRjZtWuXBg8erMGDB0uSUlJSNHjwYC1atEiSVFhY6AwmktSrVy9t3rxZWVlZio2N1fPPP69XX31VycnJzfQVAADAtczja0Zuv/121fdoktqernr77bfrL3/5i6ebAgC0EVVVVbLb7a22PbvdLh8fH5WXl6uqqqrVtnstakqtfH195e3t3eQxcIkxAKDFWJaloqIinTlzptW3GxERoRMnTvB8qgY0tVahoaGKiIhoUp0JIwCAFlMdRMLCwhQUFNRqwcDhcOj8+fMKDg6u92FbaHytLMtSWVmZSkpKJEmRkZGNHgNhBADQIqqqqpxBpHPnzq26bYfDocrKSgUEBBBGGtCUWlU/oqOkpERhYWGNPmXDTwgA0CKqrxEJCgoyPBK0pOqfb1OuCSKMAABaFNdstG3N8fMljAAAAKMIIwAAwCjCCAAAV5g2bZpsNluNz+HDh00PrU3ibhoAAGoxevRorV271qWta9euLtOVlZXN/qbb6xFHRgAAqIW/v78iIiJcPiNGjNDs2bM1d+5cdenSxflqk3379umOO+5QcHCwwsPDNXnyZH399dfOdV24cEFTpkxRcHCwIiMj9fzzz+v222/X3LlznX1sNps2bdrkMobQ0FCXJ5ufOHFC99xzj0JDQ9WpUyeNGzdOx44dc86fNm2a7r77bj333HOKjIxU586dNWvWLJc7XSoqKjR//nxFRUXJ399f3/nOd7Ru3TpZlqWbbrpJzz33nMsY8vPzW/yoEGEEANBqLMtSWeWlVvlcrKxy/rm+15h46o033pCfn58++eQTrVq1SmfOnNEPf/hDDR48WLt27VJmZqaKi4t1zz33OJd55JFHtG3bNv3+97/Xhx9+qJycHO3Zs8ej7drtdiUnJ6t9+/basWOHPvnkEwUHB2v06NGqrKx09vvoo4905MgRffTRR3rjjTf0+uuvuwSaKVOm6K233tKLL76oAwcOaOXKlWrXrp1sNpvuu+++GkeD1q5dq1tvvVU33XRT4wrmBk7TAABazUV7lWIWbW317e5/IllBfp79ynv//fcVHBzsnL7jjjskSX369NEzzzzjbP/lL3+pwYMH68knn3S2rVmzRlFRUTp06JC6deum1157TevXr9eIESMkfRNobrjhBo/Gs3HjRjkcDr366qvO22nXrl2r0NBQ5eTkaNSoUZKkjh076qWXXpK3t7f69eunMWPGKDs7WzNnztShQ4f0zjvvKCsrS0lJSZKk6OholZaWSvrmyMqiRYuUl5en+Ph42e12bdiwocbRkuZGGAEAoBY/+MEPtHLlSud0u3btNGHCBMXFxbn027t3rz766COX4FLtyJEjunjxoiorK5WQkOBs79Spk/r27evRePbu3avDhw+rffv2Lu3l5eU6cuSIc/rmm292eRJqZGSkPv/8c0nfnHLx9vbWbbfdVus2unXrpjFjxmjNmjWKj4/Xe++9p4qKCv30pz/1aKyeIowAAFpNoK+39j+R3OLbcTgcOld6Tu1D2svLy0uBvp4/prxdu3a1nppo166dy/T58+c1duxYPf300zX6RkZGun2thc1mq3E66fJrPc6fP6+4uDi9+eabNZa9/MJaX1/fGut1OBySvn18e33uv/9+TZ48Wb/61a+0du1ajR8/vsWfoksYAQC0GpvN5vHpksZwOBy65OetID+fFn83zZAhQ/S///u/io6Olo9Pze/Wu3dv+fr66tNPP1WPHj0kSf/617906NAhlyMUXbt2VWFhoXP6iy++UFlZmct2Nm7cqLCwMIWEhDRqrAMGDJDD4dC2bducp2mudOedd6pdu3ZauXKlMjMztX379kZtyxNcwAoAQBPMmjVLp0+f1oQJE/TZZ5/pyJEj2rp1q6ZPn66qqioFBwdrxowZeuSRR/SnP/1J+/bt07Rp02qEpB/+8Id66aWX9Je//EW7du3SAw884HKUY9KkSerSpYvGjRunHTt2qKCgQDk5Ofr5z3+ur776yq2xRkdHa+rUqbrvvvu0adMm5zp+97vfOft4e3tr2rRpSk1NVZ8+fZSYmNg8haoHYQQAgCbo1q2bPvnkE1VVVWnUqFEaMGCA5s6dq9DQUGfgePbZZzV8+HCNHTtWSUlJ+v73v1/j2pPnn39eUVFRGj58uCZOnKiHH37Y5fRIUFCQtm/frh49eujHP/6x+vfvrxkzZqi8vNyjIyUrV67UT37yEz300EPq16+f/vM//9PlCIwkzZgxQ5WVlZo+fXoTKuM+TtMAAHCFy2+FvVxOTk6t7X369NG7775b5/qCg4O1bt06rVu3ztm2efNmlz7dunXT1q2udxqdOXPGZToiIkJvvPGGR+Nevny5y3RAQICWLVumZcuWSfrmlFb13TTVTp48KV9fX02ZMqXObTUnwggAAJD0zQPRTp06pcWLF+unP/2pwsPDW2W7nKYBAACSpLfeeks9e/bUmTNnXJ6l0tI4MgIAgAF1nfIxadq0aZo2bVqrb5cjIwAAwCjCCACgRTXne2Fw9WmOny9hBADQIqqfkXHlbaNoW6p/vlc++dUTXDMCAGgR3t7eCg0NVUlJiaRvnpNR/YK3luZwOFRZWany8vIWfwLrta6xtbIsS2VlZSopKVFoaKjL+3A8RRgBALSYiIgISXIGktZiWZYuXryowMDAVgtA16qm1io0NNT5c24swggAoMXYbDZFRkYqLCzM5aVvLc1ut2v79u269dZbm3T64HrQlFr5+vo26YhINcIIAKDFeXt7N8svLU+2d+nSJQUEBBBGGnA11IoTaQAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIxqVBhZsWKFoqOjFRAQoISEBOXl5dXbf/ny5erbt68CAwMVFRWlefPmqby8vFEDBgAAbYvHYWTjxo1KSUlRWlqa9uzZo9jYWCUnJ6ukpKTW/hs2bNCCBQuUlpamAwcO6LXXXtPGjRv12GOPNXnwAADg2udxGFm2bJlmzpyp6dOnKyYmRqtWrVJQUJDWrFlTa/+dO3fqe9/7niZOnKjo6GiNGjVKEyZMaPBoCgAAuD74eNK5srJSu3fvVmpqqrPNy8tLSUlJys3NrXWZYcOGaf369crLy1N8fLyOHj2qLVu2aPLkyXVup6KiQhUVFc7p0tJSSZLdbpfdbvdkyPWqXldzrrOtolaeoV7uo1buo1buo1bua8laubtOm2VZlrsr/cc//qHu3btr586dSkxMdLY/+uij2rZtmz799NNal3vxxRf18MMPy7IsXbp0SQ888IBWrlxZ53YWL16s9PT0Gu0bNmxQUFCQu8MFAAAGlZWVaeLEiTp79qxCQkLq7OfRkZHGyMnJ0ZNPPqmXX35ZCQkJOnz4sObMmaMlS5Zo4cKFtS6TmpqqlJQU53RpaamioqI0atSoer+Mp+x2u7KysjRy5Ej5+vo223rbImrlGerlPmrlPmrlPmrlvpasVfWZjYZ4FEa6dOkib29vFRcXu7QXFxcrIiKi1mUWLlyoyZMn6/7775ckDRgwQBcuXNB//Md/6Be/+IW8vGpetuLv7y9/f/8a7b6+vi2yU7XUetsiauUZ6uU+auU+auU+auW+lqiVu+vz6AJWPz8/xcXFKTs729nmcDiUnZ3tctrmcmVlZTUCh7e3tyTJgzNEAACgjfL4NE1KSoqmTp2qoUOHKj4+XsuXL9eFCxc0ffp0SdKUKVPUvXt3ZWRkSJLGjh2rZcuWafDgwc7TNAsXLtTYsWOdoQQAAFy/PA4j48eP16lTp7Ro0SIVFRVp0KBByszMVHh4uCTp+PHjLkdCHn/8cdlsNj3++OM6efKkunbtqrFjx2rp0qXN9y0AAMA1q1EXsM6ePVuzZ8+udV5OTo7rBnx8lJaWprS0tMZsCgAAtHG8mwYAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGNCiMrVqxQdHS0AgIClJCQoLy8vHr7nzlzRrNmzVJkZKT8/f31ne98R1u2bGnUgAEAQNvi4+kCGzduVEpKilatWqWEhAQtX75cycnJOnjwoMLCwmr0r6ys1MiRIxUWFqbf/va36t69u7788kuFhoY2x/gBAMA1zuMwsmzZMs2cOVPTp0+XJK1atUqbN2/WmjVrtGDBghr916xZo9OnT2vnzp3y9fWVJEVHRzdt1AAAoM3wKIxUVlZq9+7dSk1NdbZ5eXkpKSlJubm5tS7zhz/8QYmJiZo1a5Z+//vfq2vXrpo4caLmz58vb2/vWpepqKhQRUWFc7q0tFSSZLfbZbfbPRlyvarX1ZzrbKuolWeol/uolfuolfuolftaslburtOjMPL111+rqqpK4eHhLu3h4eH6+9//XusyR48e1Z/+9CdNmjRJW7Zs0eHDh/XQQw/JbrcrLS2t1mUyMjKUnp5eo/3DDz9UUFCQJ0N2S1ZWVrOvs62iVp6hXu6jVu6jVu6jVu5riVqVlZW51c/j0zSecjgcCgsL0yuvvCJvb2/FxcXp5MmTevbZZ+sMI6mpqUpJSXFOl5aWKioqSqNGjVJISEizjc1utysrK0sjR450nkJC7aiVZ6iX+6iV+6iV+6iV+1qyVtVnNhriURjp0qWLvL29VVxc7NJeXFysiIiIWpeJjIyUr6+vyymZ/v37q6ioSJWVlfLz86uxjL+/v/z9/Wu0+/r6tshO1VLrbYuolWeol/uolfuolfuolftaolburs+jW3v9/PwUFxen7OxsZ5vD4VB2drYSExNrXeZ73/ueDh8+LIfD4Ww7dOiQIiMjaw0iAADg+uLxc0ZSUlK0evVqvfHGGzpw4IAefPBBXbhwwXl3zZQpU1wucH3wwQd1+vRpzZkzR4cOHdLmzZv15JNPatasWc33LQAAwDXL42tGxo8fr1OnTmnRokUqKirSoEGDlJmZ6byo9fjx4/Ly+jbjREVFaevWrZo3b54GDhyo7t27a86cOZo/f37zfQsAAHDNatQFrLNnz9bs2bNrnZeTk1OjLTExUX/+858bsykAANDG8W4aAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGNSqMrFixQtHR0QoICFBCQoLy8vLcWu7tt9+WzWbT3Xff3ZjNAgCANsjjMLJx40alpKQoLS1Ne/bsUWxsrJKTk1VSUlLvcseOHdPDDz+s4cOHN3qwAACg7fE4jCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qypc5mqqipNmjRJ6enpuvHGG5s0YAAA0Lb4eNK5srJSu3fvVmpqqrPNy8tLSUlJys3NrXO5J554QmFhYZoxY4Z27NjR4HYqKipUUVHhnC4tLZUk2e122e12T4Zcr+p1Nec62ypq5Rnq5T5q5T5q5T5q5b6WrJW76/QojHz99deqqqpSeHi4S3t4eLj+/ve/17rMxx9/rNdee035+flubycjI0Pp6ek12j/88EMFBQV5MmS3ZGVlNfs62ypq5Rnq5T5q5T5q5T5q5b6WqFVZWZlb/TwKI546d+6cJk+erNWrV6tLly5uL5eamqqUlBTndGlpqaKiojRq1CiFhIQ02/jsdruysrI0cuRI+fr6Ntt62yJq5Rnq5T5q5T5q5T5q5b6WrFX1mY2GeBRGunTpIm9vbxUXF7u0FxcXKyIiokb/I0eO6NixYxo7dqyzzeFwfLNhHx8dPHhQvXv3rrGcv7+//P39a7T7+vq2yE7VUutti6iVZ6iX+6iV+6iV+6iV+1qiVu6uz6MLWP38/BQXF6fs7Gxnm8PhUHZ2thITE2v079evnz7//HPl5+c7P3fddZd+8IMfKD8/X1FRUZ5sHgAAtEEen6ZJSUnR1KlTNXToUMXHx2v58uW6cOGCpk+fLkmaMmWKunfvroyMDAUEBOiWW25xWT40NFSSarQDAIDrk8dhZPz48Tp16pQWLVqkoqIiDRo0SJmZmc6LWo8fPy4vLx7sCgAA3NOoC1hnz56t2bNn1zovJyen3mVff/31xmwSAAC0URzCAAAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARjXqrb0AWl+Vw1JewWmVnCtXWPsAxffqJG8vm+lh4TrGPonmQhgBrgGZ+wqV/t5+FZ4td7ZFdghQ2tgYjb4l0uDIcL1in0Rz4jQNcJXL3FeoB9fvcflHX5KKzpbrwfV7lLmv0NDIcL1in0RzI4wAV7Eqh6X09/bLqmVedVv6e/tV5aitB9D82CfREggjwFUsr+B0jf99Xs6SVHi2XHkFp1tvULiusU+iJRBGgKtYybm6/9FvTD+gqdgn0RIII8BVLKx9QLP2A5qKfRItgTACXMXie3VSZIcA1XWzpE3f3MEQ36tTaw4L1zH2SbQEwghwFfP2siltbIwk1fjHv3o6bWwMz3ZAq2GfREsgjABXudG3RGrlvUMU0cH1sHdEhwCtvHcIz3RAq2OfRHPjoWfANWD0LZEaGRPB0y5x1WCfRHMijADXCG8vmxJ7dzY9DMCJfRLNhdM0AADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwqlFhZMWKFYqOjlZAQIASEhKUl5dXZ9/Vq1dr+PDh6tixozp27KikpKR6+wMAgOuLx2Fk48aNSklJUVpamvbs2aPY2FglJyerpKSk1v45OTmaMGGCPvroI+Xm5ioqKkqjRo3SyZMnmzx4AABw7fM4jCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qyptf+bb76phx56SIMGDVK/fv306quvyuFwKDs7u8mDBwAA1z4fTzpXVlZq9+7dSk1NdbZ5eXkpKSlJubm5bq2jrKxMdrtdnTp1qrNPRUWFKioqnNOlpaWSJLvdLrvd7smQ61W9ruZcZ1tFrTxDvdxHrdxHrdxHrdzXkrVyd502y7Isd1f6j3/8Q927d9fOnTuVmJjobH/00Ue1bds2ffrppw2u46GHHtLWrVv1t7/9TQEBAbX2Wbx4sdLT02u0b9iwQUFBQe4OFwAAGFRWVqaJEyfq7NmzCgkJqbOfR0dGmuqpp57S22+/rZycnDqDiCSlpqYqJSXFOV1aWuq81qS+L+Mpu92urKwsjRw5Ur6+vs223raIWnmGermPWrmPWrmPWrmvJWtVfWajIR6FkS5dusjb21vFxcUu7cXFxYqIiKh32eeee05PPfWU/vjHP2rgwIH19vX395e/v3+Ndl9f3xbZqVpqvW0RtfIM9XIftXIftXIftXJfS9TK3fV5dAGrn5+f4uLiXC4+rb4Y9fLTNld65plntGTJEmVmZmro0KGebBIAALRxHp+mSUlJ0dSpUzV06FDFx8dr+fLlunDhgqZPny5JmjJlirp3766MjAxJ0tNPP61FixZpw4YNio6OVlFRkSQpODhYwcHBzfhVAADAtcjjMDJ+/HidOnVKixYtUlFRkQYNGqTMzEyFh4dLko4fPy4vr28PuKxcuVKVlZX6yU9+4rKetLQ0LV68uGmjBwAA17xGXcA6e/ZszZ49u9Z5OTk5LtPHjh1rzCYAAMB1olXvpgFw7alyWMorOK2Sc+UKax+g+F6d5O1lc3u+CVfjmJqq8pJD63OPqbOkdbnHdO+w3vLz4fViaBsIIwDqlLmvUOnv7Vfh2XJnW2SHAKWNjdHoWyIbnG/C1TimpsrYsl+rdxTI18vSM/HS01sP6pcfHNLM4b2UemeM6eEBTUasBlCrzH2FenD9Hpdf6pJUdLZcD67fo4wt++udn7mvsDWHK6nhMZsYU1NlbNmv/95eIMcVj6d0WNJ/by9Qxpb9ZgYGNCPCCIAaqhyW0t/br9oez2z932f1joI650tS+nv7VXXlb9AW1NCYTYypqSovObR6R0G9fVbvKFDlJUcrjQhoGYQRADXkFZyucXThSvX9TrckFZ4tV17B6eYdWD0aGrOJMTXVutxj9dZZ+ubnsC73WKuMB2gphBEANZScqz+ItPZ6mnNbrTmmpvrydFmz9gOuVoQRADWEta/73VEm1tOc22rNMTVVz07uvRjU3X7A1YowAqCG+F6dFNkhQPXdDOtlU53zbfrmDpb4Xp1aYHS1a2jMJsbUVJMTo9XQHcletm/6AdcywgiAGry9bEob+80to1f+LrT932fm8F51zpektLExrfpsj4bGbGJMTeXn4+Wsc11mDu/F80ZwzWMPBlCr0bdEauW9QxTRwfW0RkSHAK28d4hS74ypd76JZ3o0NOZr8TkjqXfG6D9v7VXjCImXTfrPW3nOCNoGHnoGoE6jb4nUyJiIOp9m2tD8q3HM16LUO2P0X6P6af3OI9K/9mt+cl+ewIo2hTACoF7eXjYl9u7c6PkmXI1jaio/Hy9NTozWli37NTkxWr4EEbQh7M0AAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDKx/QAAOBaUuWwlFdwWiXnyhXWPkDxvTrJ28smSbpYWaUnt+zXsX+WKbpzkB67M0aBft5uLVvfPEmqvOTQ+txj6ixpXe4x3Tust/x83Pv/ZFO229j1Vo95Xe4xfXm6TD07BWlyYnSzjLmlx43W16gwsmLFCj377LMqKipSbGysfv3rXys+Pr7O/r/5zW+0cOFCHTt2TH369NHTTz+tO++8s9GDBgATMvcVKv29/So8W+5si+wQoLSxMfrfPV8pa3+Js33HF9K6Px/XyJgwrZ7y3XqXlVTnvNG3RCpjy36t3lEgXy9Lz8RLT289qF9+cEgzh/dS6p0xjR5zQ9tt7HovH7PD+naZpVsONHnMo2+JbHB+U9YNMzwOIxs3blRKSopWrVqlhIQELV++XMnJyTp48KDCwsJq9N+5c6cmTJigjIwM/du//Zs2bNigu+++W3v27NEtt9zSLF8CAFpa5r5CPbh+j6wr2ovOluuB9XvqXC5rf4nuemmHPv+q1KNli86W68H1e5QUE+YScqo5LOm/txdIUp2/3Bsz5urtrrx3SJ2/nOtbb0uO+cH1e/Qft/bSK9sL6pzflHHXtyxalsfXjCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qyptf8LL7yg0aNH65FHHlH//v21ZMkSDRkyRC+99FKTBw8AraHKYSn9vf01folJqrXtSn+tJYg0tKz1f5/afqlfbvWOAlVectRob+yYq+elv7dfVY6aPRtab0uO2fq/Zev7To0dd33LouV5dGSksrJSu3fvVmpqqrPNy8tLSUlJys3NrXWZ3NxcpaSkuLQlJydr06ZNdW6noqJCFRUVzumzZ89Kkk6fPi273e7JkOtlt9tVVlamf/7zn/L19W229bZF1Moz1Mt910Ktdn/5L5365z+NX2Tn47BUVuaQj91LVY5vr3FYnbVXP4vv4dK3qWM+9c8Lys4/orieHZt1vdWaMub6/hddPe6B3YJd9it31l3Xd27rWvLv4Llz5yRJltVAyLM8cPLkSUuStXPnTpf2Rx55xIqPj691GV9fX2vDhg0ubStWrLDCwsLq3E5aWlp1CObDhw8fPnz4XOOfEydO1JsvTAf9WqWmprocTXE4HDp9+rQ6d+4sm635rnguLS1VVFSUTpw4oZCQkGZbb1tErTxDvdxHrdxHrdxHrdzXkrWyLEvnzp1Tt27d6u3nURjp0qWLvL29VVxc7NJeXFysiIiIWpeJiIjwqL8k+fv7y9/f36UtNDTUk6F6JCQkhJ3VTdTKM9TLfdTKfdTKfdTKfS1Vqw4dOjTYx6MLWP38/BQXF6fs7Gxnm8PhUHZ2thITE2tdJjEx0aW/JGVlZdXZHwAAXF88Pk2TkpKiqVOnaujQoYqPj9fy5ct14cIFTZ8+XZI0ZcoUde/eXRkZGZKkOXPm6LbbbtPzzz+vMWPG6O2339auXbv0yiuvNO83AQAA1ySPw8j48eN16tQpLVq0SEVFRRo0aJAyMzMVHh4uSTp+/Li8vL494DJs2DBt2LBBjz/+uB577DH16dNHmzZtuiqeMeLv76+0tLQap4RQE7XyDPVyH7VyH7VyH7Vy39VQK5tlNXS/DQAAQMvhRXkAAMAowggAADCKMAIAAIwijAAAAKOuizCycuVKDRw40PlAl8TERH3wwQfO+eXl5Zo1a5Y6d+6s4OBg/fu//3uNB7Vdj5566inZbDbNnTvX2UatvrV48WLZbDaXT79+/ZzzqZWrkydP6t5771Xnzp0VGBioAQMGaNeuXc75lmVp0aJFioyMVGBgoJKSkvTFF18YHLEZ0dHRNfYrm82mWbNmSWK/ulxVVZUWLlyoXr16KTAwUL1799aSJUtc3oPCfvWtc+fOae7cuerZs6cCAwM1bNgwffbZZ875RmvV0Pto2oI//OEP1ubNm61Dhw5ZBw8etB577DHL19fX2rdvn2VZlvXAAw9YUVFRVnZ2trVr1y7r//2//2cNGzbM8KjNysvLs6Kjo62BAwdac+bMcbZTq2+lpaVZN998s1VYWOj8nDp1yjmfWn3r9OnTVs+ePa1p06ZZn376qXX06FFr69at1uHDh519nnrqKatDhw7Wpk2brL1791p33XWX1atXL+vixYsGR976SkpKXPaprKwsS5L10UcfWZbFfnW5pUuXWp07d7bef/99q6CgwPrNb35jBQcHWy+88IKzD/vVt+655x4rJibG2rZtm/XFF19YaWlpVkhIiPXVV19ZlmW2VtdFGKlNx44drVdffdU6c+aM5evra/3mN79xzjtw4IAlycrNzTU4QnPOnTtn9enTx8rKyrJuu+02ZxihVq7S0tKs2NjYWudRK1fz58+3vv/979c53+FwWBEREdazzz7rbDtz5ozl7+9vvfXWW60xxKvWnDlzrN69e1sOh4P96gpjxoyx7rvvPpe2H//4x9akSZMsy2K/ulxZWZnl7e1tvf/++y7tQ4YMsX7xi18Yr9V1cZrmclVVVXr77bd14cIFJSYmavfu3bLb7UpKSnL26devn3r06KHc3FyDIzVn1qxZGjNmjEtNJFGrWnzxxRfq1q2bbrzxRk2aNEnHjx+XRK2u9Ic//EFDhw7VT3/6U4WFhWnw4MFavXq1c35BQYGKiopc6tWhQwclJCRcl/WqVllZqfXr1+u+++6TzWZjv7rCsGHDlJ2drUOHDkmS9u7dq48//lh33HGHJPary126dElVVVUKCAhwaQ8MDNTHH39svFZX5Vt7W8Lnn3+uxMRElZeXKzg4WL/73e8UExOj/Px8+fn51XgRX3h4uIqKiswM1qC3335be/bscTmPWK2oqIhaXSYhIUGvv/66+vbtq8LCQqWnp2v48OHat28ftbrC0aNHtXLlSqWkpOixxx7TZ599pp///Ofy8/PT1KlTnTWpfpJzteu1XtU2bdqkM2fOaNq0aZL4O3ilBQsWqLS0VP369ZO3t7eqqqq0dOlSTZo0SZLYry7Tvn17JSYmasmSJerfv7/Cw8P11ltvKTc3VzfddJPxWl03YaRv377Kz8/X2bNn9dvf/lZTp07Vtm3bTA/rqnLixAnNmTNHWVlZNdIzaqr+35ckDRw4UAkJCerZs6feeecdBQYGGhzZ1cfhcGjo0KF68sknJUmDBw/Wvn37tGrVKk2dOtXw6K5er732mu64444GX79+vXrnnXf05ptvasOGDbr55puVn5+vuXPnqlu3buxXtVi3bp3uu+8+de/eXd7e3hoyZIgmTJig3bt3mx7a9XE3jfTNG4dvuukmxcXFKSMjQ7GxsXrhhRcUERGhyspKnTlzxqV/cXGxIiIizAzWkN27d6ukpERDhgyRj4+PfHx8tG3bNr344ovy8fFReHg4tapHaGiovvOd7+jw4cPsV1eIjIxUTEyMS1v//v2dp7Wqa3LlXSHXa70k6csvv9Qf//hH3X///c429itXjzzyiBYsWKCf/exnGjBggCZPnqx58+Y5X9TKfuWqd+/e2rZtm86fP68TJ04oLy9PdrtdN954o/FaXTdh5EoOh0MVFRWKi4uTr6+vsrOznfMOHjyo48ePKzEx0eAIW9+IESP0+eefKz8/3/kZOnSoJk2a5Pwztarb+fPndeTIEUVGRrJfXeF73/ueDh486NJ26NAh9ezZU5LUq1cvRUREuNSrtLRUn3766XVZL0lau3atwsLCNGbMGGcb+5WrsrIylxezSpK3t7ccDock9qu6tGvXTpGRkfrXv/6lrVu3aty4ceZr1eKXyF4FFixYYG3bts0qKCiw/vrXv1oLFiywbDab9eGHH1qW9c2tcj169LD+9Kc/Wbt27bISExOtxMREw6O+Olx+N41lUavL/dd//ZeVk5NjFRQUWJ988omVlJRkdenSxSopKbEsi1pdLi8vz/Lx8bGWLl1qffHFF9abb75pBQUFWevXr3f2eeqpp6zQ0FDr97//vfXXv/7VGjdu3HV7C2ZVVZXVo0cPa/78+TXmsV99a+rUqVb37t2dt/a+++67VpcuXaxHH33U2Yf96luZmZnWBx98YB09etT68MMPrdjYWCshIcGqrKy0LMtsra6LMHLfffdZPXv2tPz8/KyuXbtaI0aMcAYRy7KsixcvWg899JDVsWNHKygoyPrRj35kFRYWGhzx1ePKMEKtvjV+/HgrMjLS8vPzs7p3726NHz/e5bkZ1MrVe++9Z91yyy2Wv7+/1a9fP+uVV15xme9wOKyFCxda4eHhlr+/vzVixAjr4MGDhkZr1tatWy1JtX5/9qtvlZaWWnPmzLF69OhhBQQEWDfeeKP1i1/8wqqoqHD2Yb/61saNG60bb7zR8vPzsyIiIqxZs2ZZZ86ccc43WSubZV32qDoAAIBWdt1eMwIAAK4OhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGgDbIZrPV+1m8eLHpITa76OhoLV++3PQwADSCj+kBAGh+hYWFzj9v3LhRixYtcnlRXXBwsIlhecyyLFVVVcnHp/X+qaqsrJSfn1+rbQ8AR0aANikiIsL56dChg2w2m0vb22+/rf79+ysgIED9+vXTyy+/7Fz22LFjstlseueddzR8+HAFBgbqu9/9rg4dOqTPPvtMQ4cOVXBwsO644w6dOnXKudy0adN09913Kz09XV27dlVISIgeeOABVVZWOvs4HA5lZGSoV69eCgwMVGxsrH7729865+fk5Mhms+mDDz5QXFyc/P399fHHH+vIkSMaN26cwsPDFRwcrO9+97v64x//6Fzu9ttv15dffql58+Y5j/5I0uLFizVo0CCX2ixfvlzR0dE1xr106VJ169ZNffv2lSSdOHFC99xzj0JDQ9WpUyeNGzdOx44da44fD4ArEEaA68ybb76pRYsWaenSpTpw4ICefPJJLVy4UG+88YZLv7S0ND3++OPas2ePfHx8NHHiRD366KN64YUXtGPHDh0+fFiLFi1yWSY7O1sHDhxQTk6O3nrrLb377rtKT093zs/IyND//M//aNWqVfrb3/6mefPm6d5779W2bdtc1rNgwQI99dRTOnDggAYOHKjz58/rzjvvVHZ2tv7yl79o9OjRGjt2rI4fPy5Jevfdd3XDDTfoiSeeUGFhocuRIXdkZ2fr4MGDysrK0vvvvy+73a7k5GS1b99eO3bs0CeffKLg4GCNHj3aJVwBaCat8jo+AMasXbvW6tChg3O6d+/e1oYNG1z6LFmyxPka+oKCAkuS9eqrrzrnv/XWW5YkKzs729mWkZFh9e3b1zk9depUq1OnTtaFCxecbStXrrSCg4Otqqoqq7y83AoKCrJ27tzpsu0ZM2ZYEyZMsCzLsj766CNLkrVp06YGv9fNN99s/frXv3ZO9+zZ0/rVr37l0ictLc2KjY11afvVr35l9ezZ02Xc4eHhLm96XbdundW3b1/L4XA42yoqKqzAwEBr69atDY4NgGe4ZgS4jly4cEFHjhzRjBkzNHPmTGf7pUuX1KFDB5e+AwcOdP45PDxckjRgwACXtpKSEpdlYmNjFRQU5JxOTEzU+fPndeLECZ0/f15lZWUaOXKkyzKVlZUaPHiwS9vQoUNdps+fP6/Fixdr8+bNKiws1KVLl3Tx4kXnkZGmGjBggMt1Inv37tXhw4fVvn17l37l5eU6cuRIs2wTwLcII8B15Pz585Kk1atXKyEhwWWet7e3y7Svr6/zz9XXYFzZ5nA4PN725s2b1b17d5d5/v7+LtPt2rVzmX744YeVlZWl5557TjfddJMCAwP1k5/8pMFTJl5eXrIsy6XNbrfX6Hfl9s6fP6+4uDi9+eabNfp27dq13m0C8BxhBLiOhIeHq1u3bjp69KgmTZrU7Ovfu3evLl68qMDAQEnSn//8ZwUHBysqKkqdOnWSv7+/jh8/rttuu82j9X7yySeaNm2afvSjH0n6JixceTGpn5+fqqqqXNq6du2qoqIiWZblDFT5+fkNbm/IkCHauHGjwsLCFBIS4tFYAXiOC1iB60x6eroyMjL04osv6tChQ/r888+1du1aLVu2rMnrrqys1IwZM7R//35t2bJFaWlpmj17try8vNS+fXs9/PDDmjdvnt544w0dOXJEe/bs0a9//esaF89eqU+fPnr33XeVn5+vvXv3auLEiTWOykRHR2v79u06efKkvv76a0nf3GVz6tQpPfPMMzpy5IhWrFihDz74oMHvMWnSJHXp0kXjxo3Tjh07VFBQoJycHP385z/XV1991fgCAagVYQS4ztx///169dVXtXbtWg0YMEC33XabXn/9dfXq1avJ6x4xYoT69OmjW2+9VePHj9ddd93l8oC1JUuWaOHChcrIyFD//v01evRobd68ucFtL1u2TB07dtSwYcM0duxYJScna8iQIS59nnjiCR07dky9e/d2nkrp37+/Xn75Za1YsUKxsbHKy8vTww8/3OD3CAoK0vbt29WjRw/9+Mc/Vv/+/TVjxgyVl5dzpARoATbryhOqANAI06ZN05kzZ7Rp0ybTQwFwjeHICAAAMIowAgAAjOI0DQAAMIojIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMCo/w84A34Aq5eZXgAAAABJRU5ErkJggg==\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": 27,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG6CAYAAAALTELXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVPklEQVR4nO3dd3wUdf4/8NfMbE12UyEJHQQhhBJAiCAgRUW/gJ69nHhYED31OEUFPRue51kRK3J62OXnnYKoByL2ShUNKh1pAdJ7ts98fn9ssiQkQLI7ySaZ1/PxyEMyM7vz2beb7Cufz2c+IwkhBIiIiIgMSo52A4iIiIiiiWGIiIiIDI1hiIiIiAyNYYiIiIgMjWGIiIiIDI1hiIiIiAyNYYiIiIgMjWGIiIiIDI1hiIiIiAytVYWhf/3rX7jqqquOe0xJSQluv/12jBgxAllZWXjwwQfhdrtbqIVERETU3pii3YAab7/9Np5++mkMHz78uMfNmjULbrcbr732GsrLy3HPPffA5XLhsccea6GWEhERUXsS9TCUl5eHBx54AOvWrUPPnj2Pe+xPP/2E9evXY+XKlejduzcA4O9//ztmzJiB2bNnIzU1tQVaTERERO1J1IfJfvvtN5jNZnz44YfIzMw87rEbN25Ex44dQ0EIALKysiBJEn788cfmbioRERG1Q1HvGZo4cSImTpzYqGPz8vLQqVOnOtssFgsSEhJw+PDh5mgeERERtXNR7xlqCrfbDYvFUm+71WqF1+sN+3mFEJE0i4iIiNqwqPcMNYXNZoPP56u33ev1IiYmJuznlSQJ5eVuqKoWSfMMT1FkxMXZWcsIsY76YS31w1rqg3XUT3y8HbKsT59OmwpDaWlp+Oyzz+ps8/l8KC0tRUpKSkTPraoaAgG+MfXAWuqDddQPa6kf1lIfrGPk9BzUaVPDZCNGjEBubi727dsX2rZ+/XoAwCmnnBKtZhEREVEb1qrDkKqqKCgogMfjAQBkZmZi2LBhuO2227B582asXbsW999/P84//3xeVk9ERERhadVh6PDhwxgzZgxWrlwJIDi35/nnn0fXrl0xffp03HrrrTj99NMxb9686DaUiIiI2ixJ8FIqAEBJSRXHbyNkMslITIxlLSPEOuqHtdQPa6kP1lE/SUmxUBR9+nRadc8QERERUXNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDi3oY0jQNzz77LMaOHYshQ4bg+uuvx4EDB455fFFREW6//XaMHDkSp556Km677Tbk5eW1YIuJiIioPYl6GFq4cCGWLFmChx56CO+88w40TcOMGTPg8/kaPP7WW2/FoUOH8Oqrr+LVV1/FoUOHcPPNN7dwq4mIiKi9iGoY8vl8eOWVVzBr1iyMHz8e6enpWLBgAXJzc7F69ep6x5eXl2P9+vW4/vrr0b9/f2RkZGDmzJn45ZdfUFpa2vIvgIiIiNo8UzRPvm3bNlRVVWHUqFGhbXFxccjIyMCGDRswderUOsfbbDbExsZi+fLlyMrKAgB88MEH6NWrF+Li4iJqi6JEvZOszaupIWsZGdZRP6ylflhLfbCO+pEk/Z4rqmEoNzcXANCpU6c621NSUkL7arNYLHj00Udx//33Y/jw4ZAkCSkpKXjrrbcgy5G9seLi7BE9no5gLfXBOuqHtdQPa6kP1rF1iWoYcrvdAIIhpzar1YqysrJ6xwshsHXrVgwdOhQzZsyAqqpYsGABbrrpJvy///f/4HA4wm5LebkbqqqF/XgK/qUTF2dnLSPEOuqHtdQPa6kP1lE/8fH2iDtCakQ1DNlsNgDBuUM1/wYAr9cLu71+av7444/x1ltv4csvvwwFn0WLFmHChAl47733cPXVV4fdFlXVEAjwjakH1lIfrKN+WEv9sJb6YB0jJ4R+zxXVQcua4bH8/Pw62/Pz85Gamlrv+I0bN6JXr151eoDi4+PRq1cv7Nu3r3kbS0RERO1SVMNQeno6HA4H1q1bF9pWXl6OLVu2YMSIEfWOT0tLw759++D1ekPbXC4XcnJy0LNnz5ZoMhEREbUzUQ1DFosF06ZNw5NPPonPP/8c27Ztw2233Ya0tDRMmjQJqqqioKAAHo8HAHD++ecDCK41tG3bNmzbtg2zZ8+G1WrFhRdeGMVXQkRERG1V1K/tmzVrFi6++GLce++9uOKKK6AoChYvXgyz2YzDhw9jzJgxWLlyJYDgVWZLliyBEALTp0/HNddcA7PZjCVLlsDpdEb5lRAREVFbJAmh5xSktqukpIqT2SJkMslITIxlLSPEOuqHtdQPa6kP1lE/SUmxuq3XFPWeISIiIqJoYhgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ4t6GNI0Dc8++yzGjh2LIUOG4Prrr8eBAweOebzf78f8+fNDx0+bNg1bt25twRYTERFRexL1MLRw4UIsWbIEDz30EN555x1omoYZM2bA5/M1ePy8efOwbNky/POf/8TSpUuRlJSE66+/HhUVFS3cciIiImoPohqGfD4fXnnlFcyaNQvjx49Heno6FixYgNzcXKxevbre8QcOHMDSpUvx8MMPY+zYsejduzf+8Y9/wGKx4Ndff43CKyAiIqK2zhTNk2/btg1VVVUYNWpUaFtcXBwyMjKwYcMGTJ06tc7x33//PZxOJ04//fQ6x3/xxRcRt0VRot5J1ubV1JC1jAzrqB/WUj+spT5YR/1Ikn7PFdUwlJubCwDo1KlTne0pKSmhfbXt2bMH3bp1w+rVq/HSSy8hLy8PGRkZuOuuu9C7d++I2hIXZ4/o8XQEa6kP1lE/rKV+WEt9sI6tS1TDkNvtBgBYLJY6261WK8rKyuodX1lZiX379mHhwoWYM2cO4uLi8OKLL+KPf/wjVq5cieTk5LDbUl7uhqpqYT+egn/pxMXZWcsIsY76YS31w1rqg3XUT3y8HbKsTw9bWGHof//7HyZNmlQvxDSVzWYDEJw7VPNvAPB6vbDb66dmk8mEyspKLFiwINQTtGDBAowbNw7vv/8+ZsyYEXZbVFVDIMA3ph5YS32wjvphLfXDWuqDdYycEPo9V1iRas6cORg9ejTmzZuHzZs3h33ymuGx/Pz8Otvz8/ORmppa7/i0tDSYTKY6Q2I2mw3dunVDTk5O2O0gIiIi4worDH3xxRe49tprsXbtWlx22WWYPHkyFi9ejIKCgiY9T3p6OhwOB9atWxfaVl5eji1btmDEiBH1jh8xYgQCgQB++eWX0DaPx4MDBw6gR48e4bwUIiIiMriwwlBaWhr+/Oc/Y9WqVXj77bcxfPhwvPzyy5gwYQJuvPFGrF69GoFA4ITPY7FYMG3aNDz55JP4/PPPsW3bNtx2221IS0vDpEmToKoqCgoK4PF4AADDhw/Haaedhrlz52Ljxo3YtWsX5syZA0VR8Ic//CGcl0JEREQGJwmhz6jbL7/8gscffxwbNmwAAHTo0AHTp0/HtddeC0VRjvk4VVXx1FNPYdmyZfB4PBgxYgTuv/9+dO3aFTk5OTjjjDPwyCOP4MILLwQQnET95JNPYtWqVfB4PBg2bBj+9re/oU+fPhG1v6SkiuO3ETKZZCQmxrKWEWId9cNa6oe11AfrqJ+kpFjdliiIKAwdPHgQH3zwAT744APs378f3bt3x0UXXYTx48fjq6++wsKFC3H22Wfjscce06WxzYlvzMjxh1wfrKN+WEv9sJb6YB31o2cYCutqsnfffRcffPABNm3aBKvVinPOOQcPP/wwhg8fHjqmb9++KCkpwTvvvNMmwhAREREZU1hh6L777kNmZibmzZuHyZMnw+FwNHhcv379cNlll0XUQCIiIqLmFPY6Q3369IGqqqH5QB6PB36/H06nM3Tc+eefr0sjiYiIiJpLWINtPXv2xAMPPIBLL700tG3Tpk0YNWoUHnvsMWgax0GJiIiobQgrDD377LP48MMP69xINSMjA3fccQf++9//4t///rduDSQiIiJqTmENk3300UeYO3cuLr/88tC2hIQEXH311TCZTHjjjTcwc+ZM3RpJRERE1FzC6hkqKSlBt27dGtx30kknNXjHeSIiIqLWKKwwdNJJJ+GTTz5pcN8XX3zBW2MQERFRmxHWMNmf/vQn3HXXXSgtLcWZZ56J5ORkFBcX48svv8THH3+MRx55RO92EhERETWLsMLQ+eefj6qqKixcuBCrV68ObU9MTMR9993HS+qJiIiozQgrDAHAlVdeiT/+8Y/Ys2cPSktLERcXh5NOOgmyrM/S2EREREQtIewwBACSJOGkk07Sqy1ERERELS6sMFRcXIyHH34YX331FdxuN46+16skSdiyZYsuDSQiIiJqTmGFob///e/48ssvMWXKFKSlpXFojIiIiNqssMLQN998g7/97W+8CSsRERG1eWF16ZjN5mMuukhERETUloQVhs466yz873//07stRERERC0urGGyjIwMPP300zhw4AAyMzNhs9nq7JckCTfffLMuDSQiIiJqTpI4+lKwRkhPTz/+k0oStm7dGnajoqGkpAqBgBbtZrRpJpOMxMRY1jJCrKN+WEv9sJb6YB31k5QUC0XR5wKusHqGtm3bpsvJiYiIiKIt4khVUVGB3bt3w+fzQVVVPdpERERE1GLCDkPr1q3DJZdcgqysLJx77rnYuXMnbr/9djz66KN6to+IiIioWYUVhtasWYPrrrsONpsNd9xxR2gF6vT0dLzxxht49dVXdW0kERERUXMJKww9/fTTOOOMM/Dmm29i+vTpoTB04403YsaMGXj33Xd1bSQRERFRcwkrDG3duhUXXXQRgOCVY7WNHj0aBw8ejLxlRERERC0grDDkdDpRUFDQ4L7Dhw/D6XRG1CgiIiKilhJWGDrjjDOwYMEC/PLLL6FtkiQhNzcXixYtwvjx4/VqHxEREVGzCmudodtvvx3Z2dm49NJL0aFDBwDA7NmzkZubi06dOmH27Nm6NpKIiIiouYQVhuLj4/Huu+9i+fLlWLt2LUpLS+F0OnHVVVfhwgsvhN1u17udRERERM0irDAEABaLBZdeeikuvfRSPdtDRERE1KLCCkPLly8/4THnn39+OE9NRERE1KLCCkN33XVXg9slSYKiKFAUhWGIiIiI2oSwwtDnn39eb5vL5cLGjRvx8ssv44UXXoi4YUREREQtIaww1KVLlwa3n3zyyfD7/XjooYewZMmSiBpGRERE1BIivmv90fr164fffvtN76clIiIiaha6hiGfz4f33nsPycnJej4tERERUbMJa5hs4sSJ9e5JpmkaSkpK4PV6MXfuXF0aR0RERNTcwgpDWVlZ9cIQADgcDkyYMAGnnXZaxA0jIiIiaglhhaFHH31U73YQERERRUVYYejQoUNNOr5z587hnIaIiIio2ek2Z+h4tm7dGs5piIiIiJpdWGHo6aefxgMPPIABAwbgvPPOQ2pqKkpKSvDFF1/g448/xp///OdjrkVERERE1JqEFYY++OADTJgwod7cocmTJyM5ORmbNm3CLbfcoksDiYiIiJpTWOsMrVmzBlOnTm1w3+mnn44ff/wxokYRERERtZSwwlBiYiKys7Mb3LdmzRqkpqZG1CgiIiKilhLWMNnFF1+MF198EW63GxMnTkRSUhIKCwuxatUq/L//9/9w33336d1OIiIiomYRVhi66aabUFFRgddeew2LFy8GAAghYLfbcdttt+Hyyy/XtZFEREREzSWsMCRJEu666y7cdNNN+Pnnn1FWVobExEQMGTIEDodD7zYSERERNZuwwlANh8OBlJQUAMCQIUMQCAR0aRQRERFRSwk7DH3wwQeYP38+CgoKIEkS3n33XTz33HMwm82YP38+LBaLnu0kIiIiahZhXU22cuVKzJ07FyNHjsRTTz0FTdMAAGeddRa+/vprLFy4UNdGEhERETWXsHqGFi1ahMsvvxzz5s2Dqqqh7RdddBGKi4vx3//+F7feeqtebSQiIiJqNmH1DO3ZswdnnXVWg/syMzORl5cXUaOIiIiIWkpYYSg5ORm7d+9ucN/u3buRnJwcUaOIiIiIWkpYYWjy5Ml49tlnsWrVKvh8PgDBy+1//fVXLFy4EOecc46ujSQiIiJqLmHNGbr11luxY8cO3HrrrZDlYJ666qqr4HK5MHz4cPz1r3/VtZFEREREzSWsMGSxWPDvf/8b33//PdauXYvS0lI4nU5kZWVh3LhxkCRJ73YSERERNYuwwtB1112HGTNmYPTo0Rg9erTebSIiIiJqMWHNGdq0aRN7f4iIiKhdCCsMjR07Fh9++CH8fr/e7SEiIiJqUWENk1mtVnz44Yf4+OOP0bt3b8TExNTZL0kSXn/9dV0aSETNo6ZzV4jotoOIKNrCCkO5ubkYOnRo6Htx1G/To78notZH1QRkmcPdRESNDkOrV6/GyJEjERcXhzfffFO3Bmiahueffx7vvvsuKioqMGLECNx///3o1q3bCR/74Ycf4s4778Tnn3+Orl276tYmIiMQAvAHNJiVsEbLiYjajUb/FvzrX/+KvXv31tn28ssvo6ioKKIGLFy4EEuWLMFDDz2Ed955B5qmYcaMGaHFHI/l4MGD+Pvf/x7RuYmMzutXwWshiMjoGh2Gjh76UlUVTz31FHJzc8M+uc/nwyuvvIJZs2Zh/PjxSE9Px4IFC5Cbm4vVq1cf83GapuHOO+/EgAEDwj43EQF+v8Y5Q0RkeGHNGaoR6dygbdu2oaqqCqNGjQpti4uLQ0ZGBjZs2ICpU6c2+LhFixbB7/fjlltuwdq1ayNqQw2FQwURq6khaxmZlqqjpGmADGgCsJjb5/8zvif1w1rqg3XUj5692hGFoUjV9Cp16tSpzvaUlJRj9jht3rwZr7zyCt577z3k5eXp1pa4OLtuz2V0rKU+mruOXr8Kv5BgtpiQGGdr1nNFG9+T+mEt9cE6ti5RDUNutxtA8PYetVmtVpSVldU73uVy4Y477sAdd9yBnj176hqGysvdUFVNt+czIkWRERdnZy0j1FJ1VDUNZeUeyJIESVWb7TzRxPekflhLfbCO+omPt4fujxqpiMNQJCtR22zBv0Z9Pl/o3wDg9Xpht9dPzf/4xz/Qq1cvXH755WGf81hUVUMgwDemHlhLfTR3HVVNQFUFAkJDlccPq0lptnNFG9+T+mEt9cE6Rk7P+Y5NCkM333xzvV6cG2+8EWazuc42SZLw2WefnfD5aobH8vPz0b1799D2/Px89OvXr97xS5cuhcViCa1xpFb/NTt16lTceOONuPHGG5vycogIwV8oHk8ANqeJa4QRkSE1OgxdcMEFup88PT0dDocD69atC4Wh8vJybNmyBdOmTat3/NFXmGVnZ+POO+/ESy+9hL59++rePiKj8AU0BFQNChdhJCIDanQYeuSRR3Q/ucViwbRp0/Dkk08iKSkJXbp0wRNPPIG0tDRMmjQJqqqiuLgYTqcTNpsNPXr0qPP4mknWnTt3RkJCgu7tIzIKVRPw+gKItZt5qT0RGU7Ur+2bNWsWLr74Ytx777244ooroCgKFi9eDLPZjMOHD2PMmDFYuXJltJtJ1O65fSqDEBEZkiQ4SQAAUFJSxclsETKZZCQmxrKWEWqpOqqaQFG5B5oW/BUgSxLinZZ2NZGa70n9sJb6YB31k5QUq9t6TVHvGSKi1kETAm5PIKIrRImI2iKGISIK8QU0+APtc80hIqJjYRgiohBNE3B72TtERMbCMEREdXj8KlSNcxmIyDgYhoioDlWt6R2KdkuIiFoGwxAR1ePxqlyNmogMg2GIiOoJaBo8fg6VEZExMAwRUT1CgENlRGQYDENE1CB/QIOPi8IRkQEwDBFRgzRNwMPeISIyAIYhIjomj1+FxonURNTOMQwR0TEFe4e4IjURtW+maDeAiFqWJgT251XA41Ph9anomGiHfIyxMCEAty+AGJuJd7SnVq3mfV3p8sMRY0b3VOcx39dER2MYIjKQrXuLsWLtPuQWu5AcZ0Ol24cYmxnjMzujd9eEBh8TUAW8fg0WEzuSqXWq/b5WVQFFkZCWFIMpI3ugf8+kaDeP2gD+diMyiK17i/H6J9uRU1AJq1lBbIwZFrOC3GI33v9uD3bnlDb4uCMTqflXNrU+R7+v4xwWWM0Kcgqq8Pon27F1b3G0m0htAMMQkQFoQmDF2n3w+AJIcFhhMSuQJQlmk4K4GDO8fg1fZR865mRpb0CFqnGcjFqXY72vLWYFCQ4LPD4VK9bu40UAdEIMQ0QGsD+vArnFLsTazPV6eCRJQoxVQWGpG4cLqxp8vKYKeLx+XmZPrcqJ3texNhNyi13Yn1cRpRZSW8EwRGQAlS4/VFXAdIx5P4oiQ9UAlyfQ4H4BwM37lVErc6L3tckkQ1UFKl3+Fm4ZtTUMQ0QG4IgxQ1EkBI6xorSqalBkIMZ27GsqApoGj4+X2VPrcaL3dSCgQVEkOGLMLdwyamsYhogMoHuqE2lJMajyBOr17ggh4PKq6JBgR6cOscd8DiEAlzcAcKiMWokTva+rPAGkJcWge6ozSi2ktoJhiMgAZEnClJE9YLMoKK30wVe9srQ/oKLc5YfVLGN8ZucTrssSCAh4/ewdotbhWO9rn19FaaUPNouCKSN7cL0hOiGGISKD6N8zCdPP7oeuHWPh9auocvnh86tIS7LjgjG9jrnOUG2aEMF5RfxsoVbi6Pd1eaUPXr+Krh1jMf3sflxniBqFiy4SGUj/nkno1yOx0StQN8TnV+H2BhBj5arU1DrUfl9zBWoKB8MQkcHIkoSeaXFQNYGicg+0Jq4fJATgcgdgNZsg87OGWoma9zVRODhMRkRN5lc1VLl9XJWaiNoFhiEiCovbp8IX4GRqImr7GIaIKCyaJlDp9nMyNRG1eQxDRBS2msnUHC0joraMYYiIwlYzmZo3wiSitoxhiIgi4lc1VLoDnExNRG0WwxARRcztDXBlaiJqsxiGiChimiZQ5fYheH97IqK2hWGIiHThC2io8nAyNRG1PQxDRKSLmrvaB1T2DhFR28IwRES6UVUBl8fPydRE1KYwDBGRrjw+FX6uTE1EbQjDEBHpStUEqtx+yLyLKxG1EQxDRKQ7b0CDiytTE1EbwTBERLrTNIFKl5+TqYmoTWAYIqJmEVA1lLt80W4GEdEJMQwRUbPx+VRUuv0cLiOiVo1hiIiajQDg9gUQCHC4jIhaL1O0G0BELUsIgRVr9mFvbgVOG5iGzh1im/V8qipQ6fYhwWkFb25PRK0RwxCRwezNrcCyb34HAPz6exGm/186enWKa9Zzev0aPH4VVpPSrOchIgoHh8mIDCYpzgarORhKfAENr328DbsPlTXrOTURXHtI8EauRNQKMQwRGUx8rAW3XDQIZlPwx98f0PDGx9uxK6d5A5HPr6HSHeCtOoio1WEYIjKgAT2T8JeLBh8JRKqGNz7Zhu37S5r1vG5vAF5/oFnPQUTUVAxDRAbVv0cirv6/dFiqA1FAFXhr9Q5s2VvcbOfUNIEKF4fLiKh1YRgiMrDeXeJxzeT+oTlEqiaw5NMdyN5V2GznDAQ0VLl5qw4iaj0YhogMrkeaE9dOSYfNEgxEmgD++8UubNyW3yznC609xFt1EFErwTBEROiW4sSMqRmIsQVX2xAAln3zO77bfLhZzhdce8jPydRE1CowDBERAKBzh1hcf24GnDHm0LaVa/fhs40HIJphtUSvX+VkaiJqFRiGiCgkNTEGN5w3AIlOa2jbF5sO4qMf9kLTORBpmkAF1x4iolaAYYiI6kiKs2HmeQPQMcEe2rb2tzy89+VuqJqm67n8fg0VVbyRKxFFF8MQEdUTH2vBzPMy0LXjkfuW/byrEG9+sgM+v6rruTw+FW4vry4jouhhGCKiBsXazLhuSgZ6dzly37IdB0qxeMVWVHn8up1HEwKVbt7Znoiih2GIiI7JalEw/Zx0DDwpKbTtQH4l/vXBbyip8Oh2noCqodzl5fwhIooKhiEiOi6TIuPyiSfj1IzU0LbCMg8WLf8NhwqrdDuPl/OHiChKGIaI6IRkWcJ5o3vizOFdQ9sq3H689NFv2JlTqtt53L4AKlwMRETUshiGiKhRJEnCxGFdceHpJ0GuDis+v4bXP96OH7frs1q1EIDLE+CCjETUohiGiKhJhqen4Kqz+4Vu8KoJgaVf/45PN+izOKMmBKo8AVS5fewhIqIWwTBERE3Wr3sirj83Aw77kdWqv/zpIP775S4E1MjXItK04BVmHDIjopYQ9TCkaRqeffZZjB07FkOGDMH111+PAwcOHPP4nTt3YubMmTj11FMxatQozJo1C4cOHWrBFhMRAHTp6MCfz6+7OGP2riIs/t9WVLojv/Q+2EPkR3mVL+LnIiI6nqiHoYULF2LJkiV46KGH8M4770DTNMyYMQM+X/1fgCUlJbjmmmtgs9nw5ptv4uWXX0ZxcTFmzJgBr9cbhdYTGVui04Yb/zAAJ3U+shbRvrwKvLj8V+SVuCJ+fiEAlzeAsiofwMvuiaiZRDUM+Xw+vPLKK5g1axbGjx+P9PR0LFiwALm5uVi9enW94z/77DO4XC48/vjj6Nu3LwYOHIgnnngCu3fvxqZNm6LwCojIbjXhmsnpGN6vY2hbSYUXi5b/hu37SyJ+fiEAjzeA0kof1yEiomZhiubJt23bhqqqKowaNSq0LS4uDhkZGdiwYQOmTp1a5/hRo0Zh4cKFsNlsoW2yHMxz5eXlEbVFUaLeSdbm1dSQtYxMS9VR0jQoiqTLnBxFUXDxhN5ISYrBx2v2QSB4V/o3PtmOySN7YGxmp4ivDvOrGspdfsTFWmBuZG34ntQPa6kP1lE/es4njGoYys3NBQB06tSpzvaUlJTQvtq6du2Krl271tn20ksvwWazYcSIERG1JS7OfuKDqFFYS300dx29fhU+IUHPm9GfN64PuneKx6sf/QavX4UQwIo1+1BU4cWV56TDbFIiPocKCXaLGY4Yc6MDFt+T+mEt9cE6ti5RDUNutxsAYLFY6my3Wq0oKys74ePffPNNvPXWW7j33nuRlJR0wuOPp7zcDVWHq2CMTFFkxMXZWcsItVQdVU1DebkHmqbv0FOv1Fj8+YIBeP3j7SipCM7lW/trLg7mV+Kqs/si3mGN+BylZRJsFgXOGDMU+dh/YfM9qR/WUh+so37i4+2h0aFIRTUM1Qx3+Xy+OkNfXq8XdvuxU7MQAs888wxefPFF/PnPf8ZVV10VcVtUVUMgwDemHlhLfTR3HVVNQFWF7mEIAFISYvDn8wdiyac7sDe3AkDwnmbPvvcLrjjzZPTqFHeCZzgRgcqABo9XhTPWDJvZdNw1jvie1A9rqQ/WMXJ69mpHddCyZngsP7/u6rX5+flITU1t6CHw+/248847sWjRItx999249dZbm7uZRBQGh92Ma6f0r3NPs0q3H4v/txU//JqrywKNAVVDWaWPCzQSUUSiGobS09PhcDiwbt260Lby8nJs2bLlmHOA5syZg1WrVmH+/Pm4+uqrW6ilRBQOkyLjD2N64fyxvaBU38NDEwL/+2Ev3v1yN3x+NeJzaJpAhZvrERFR+KI6TGaxWDBt2jQ8+eSTSEpKQpcuXfDEE08gLS0NkyZNgqqqKC4uhtPphM1mw7Jly7By5UrMmTMHWVlZKCgoCD1XzTFE1Ppk9U9FWlIMlny6A+Wu4IKMP+8qRG6xC38882R0SIhsMmnNekSaJhDnsEACu4mIqPGifm3frFmzcPHFF+Pee+/FFVdcAUVRsHjxYpjNZhw+fBhjxozBypUrAQD/+9//AACPP/44xowZU+er5hgiap26pzpx84WD0LOTM7Qtt9iFF97/Fb/8XhTx8wsBuH0qSiu8UJthHhQRtV+S0GPgvh0oKaniZLYImUwyEhNjWcsItVQdVU2gqBmuJjvxeTV8sv4Avtt8uM72UQPS8H8ju8Okw/orZkWGM9aCWLsJCQl8T+qBP9/6YB31k5QUq9t6TVEdJiMi41FkGeec2h0OmwmfbzoIf/UHwprfcrE/rwKXn3kykuNs0ITA4cIquDwBxNhM6NQhFnIjZ0n7VQ1llV6omgans/nWc9GEwP68ClS6/HDEmNE91dnoNrYlAU3Dul/y4PJqiLHKOKVfR5h0uqSZqDVgGCKiFrU7pxRfZR9CYakbZkWGqmnQqv9APlhYheeX/oLTBqRif0ElCkvdUDVAkYEOCXaMz+yM3l0TGnUeVROodPtRUOaG8KswK/ouMLl1bzFWrN2H3GIXVFVAUSSkJcVgysge6N8zsnXPWpNV6/ZhxZp9cHsDEAAkAG9ZTZgyqgfOObVHtJtHpAtGeyJqMbtzSvH+d3uQW+yCxawgzmFBktNaZ2jM61fx5c+HsC+3AmaTDEeMGRazgtxiN97/bg9255Q26ZyqKlBW6UVJhRcBVYv4tiBAMAi9/sl25BRUwlr9OqxmBTkFVXj9k+3Yurc44nO0BqvW7cPSr39HlScAWZZgUiTIsoQqTwBLv/4dq9bti3YTiXTBMERELUITAl9lH4LXryIuxgKzSYEsSbCYTeiYYIPVrNRZKyigCpRW+hBQNZhNCuJizPD6NXyVfQhaE7t4NCHg8akoqfCirMoLVRVhhyJNCKxYuw8eXwAJDiss5prXoSDBYYHHp2LF2n1NbmNrE9A0rFizD6omYK4OQbIkQ5YlmBUJqiawYs0+BDTOe6G2j2GIiFrE4cIqFJa6EWM11QsikiQhLtYMm0WpM+cmoAoUlnpQ4QquIRRjVVBY6sbhwqqw2qBqAi5PAMUVHpRXBa86a2om2p9XgdxiF2Jt9e+NJkkSYm0m5Ba7sD+vIqw2thbrt+TB7Q3AJEsNvk6TLMHtDWD9lrwotZBIPwxDRNQiXJ5AcP7PMa7+UBQZQgPMpmCgqK3C5UdRmQcAoGrB54qEqglUVYeiCpcfmmh8KKp0+aGqAiZTw6/DZJKhqgKV1esptVXF5R4IAMdcskkCRPVxRG0dwxARtYgYmwmKjGPenFJVNSgmCSaTghibCclxNsjykU9iX0BDYakHqqrBblV0aZOqBidZF5UHQ1Gwp+j4qcgRY4aiSMe8LDoQ0KAoEhwxZl3aGC1JcbZgDjrWaJ8I5qSkOC52S20fwxARtYhOHWLRIcEOl1etd18yIQRcXhWpiXakJQWPsZhlpCTYYbMcCT4CwYUVv/zpUGjoTA81oai43IPSCg/8qlYniNXWPdWJtKQYVHkCDb6OKk8AaUkx6J7qbPDxbUVWRirsVhMCmmjwdQY0AbvVhKyMhu8jSdSWMAwRUYuQJQnjMzvDapZR7vLDH1ChCQF/QEW5yw+rWcaEIV3qHKNqGuIdFjiP6mXZuq8Ez7y7Gdm7CnW54WsNVRNwV0+0LikPhqKje4pkScKUkT1gsygorfTB5w++Dp9fRWmlDzaLgikje7T59YZMsowpo3pAkSX4VQFNE9CEBk0T8KsCiixhyqgeXG+I2gWuQF2Nq4FGjiur6qO9r0Bde52hY60h1NAxCQ4r/AENh4pcdZ4vo2ci/jCmF5wxlnrnUhQJ8fExKCsLrgXUVLIswWZREGszw6TIdYKXkdcZsnOdobDx96R+9FyBmmGoGt+YkeMPuT7aexgC0KjVpRs6RgKwbmseVq3dD1+t2titCqaM6omhJ3eo05MTaRiq/Tw2swK7zQxzrVBkpBWof9xewBWodcDfk/phGGoGfGNGjj/k+jBCGIpUcbkHS7/+HXsOl9fZ3qdLPM4f2ys0qVevMFSjJhTFNNBT1N7x51sfrKN+9AxDjPZE1OYkxdlw3dT++MOYXrCaj0yw3nWwDM+8uxnfZB+C2gyLAapq7UvyfRDHvNSKiNoShiEiapNkScKpGan46yWD0a97Qmi7X9Wwat1+vLDsV+zLbZ6FD49cfeaFx69CIPwVrYko+hiGiKhNS3BY8aez++HyM/og1n7kqrPcYhcWvv8r3vx4K6o8zbMAoj+goazSi8Ky4CX5Hr8KAE1e1ZqIoot3rSeiNk+SJAzu3QF9uiTgk/X7sWFbfmjf99mH8NO2fJw1ohtGpKccc/2gcAkR7Clyqyo8fhVVsgy7VYHNaoIiG2teEVFbxZ4hImo3YmwmXHD6SZh5XgZSE+2h7S5vAB98twcvLm++oTMgGIz8qoZyV3BV6/Iqb/VaRc12SiLSAcMQEbU7PdPicMtFgzBlVA9Ya61gfbCwCv/68Df894tdKKvSbwXrhtRMti6p8KK43AOPL8C5RUStFIfJiKhdUmQZpw/pjLHDuuGd1dvw887C0L6fdxXit73FOD2zM8ZmdoLFpM+9zhqiaQJeTcAX8EFxy7CYZVgtCqwmGZIkgaNoRNHHniEiatcSnFZccebJuP7cDHRKjglt9wc0fP5jDp76TzY27Sho9vWWhAACqgaXJ4DSCi8Ky7wor/LBF+Cka6JoYxgiIkPo1SkON18wCOeP7YVY25FO8fIqH977ajdeeP8X7DhQ2iITnmuCUZUngNIKH4rKPahy+6EJDqMRRQOHyYjIMGRZQlb/VAzunYyvfjqI73/JhVrdI3S4yIXXPt6G3l3icHZWd3Tt6GiRNmlCQAsI+AMaqryB4ArXVjNMpmAo4jAaUfNjzxARGY7NYsI5p/bA7MsyMbh3cp19uw+WY+H7v+LtT3cgv8Tdou2qmXRdVOFBUbkHlW4/ArwajajZsWeIiAwr0WnD5WecjDGDOmHV+v34/dCRe539tqcYW/YWY0ifDph4SlckV9/vrCVomoBPE/D5NVTJAZgVCTarCVazwrWLiJoBwxARGV7XFAeum9IfO3PKsHr9fhwqcgEIDlH9tLMQ2bsKMbRvR0wY2iV0E9iWUvtqNFmWYFZkWM0KzGYZJkWGBA6lEUWKYYiICMFVrPt2S0CfrvH4bU8xPt1wAIVlHgCAJoAftxfgpx0FGHJyR4wf2hkd4u0neEZ91ax0raoqPD4VsixBkSVYzAqs1cFIkTnPiCgcDENERLXIkoRBJyUjo2cSsncV4osfc1Bc4QUQDEWbdhTgp50FGNw7GeOGdEFaUswJnrF5aJqApgUnXrs8wXbLsgSbRYHVrMBkkiEF+40YjohOgGGIiKgBiixhWN+OyOyTjJ93FuLLTQdDoUgIIHtXEbJ3FSG9eyLGDemMHmnOqLVVCEAVAmp1OKqSA5AlCWZFgskkQ6nuNTIpUmjOEQMS0REMQ0REx6HIMk7pl4IhJ3dA9q4ifPnTQRRVD58BwLb9Jdi2vwQ90pw4fXAn9OuRCDnKl39pmoAGgYAKwKdCkoLDgJIUDHkWU3DOkTk0tMYeJDI2hiEiokZQZBnD+nbEkD4d8MvvRfj650PILXaF9u/LrcCbuRXoEG/D6EGdMLRvh2a9zUdTCIHQFWiqGrxKTfIE112SEOxBMptkmEyce0TGxDBERNQEsiwhs08HDO6djO0HSvHNz4ewN7citL+wzIMPvtuD1RsOIKt/Ck7NSEWCwxrFFjesZkI2qnuQ3NU9SDVzj8yKHLpirSYcEbVXDENERGGQJAnp3ROR3j0R+/Mq8O3mw9iypxg1nSlubwBf/3wI32YfQkbPJIwckIZenZyt+nYbR889kryAJAeH10yKjAAkeNx+AIBJkSBLEofZqF1gGCIiilD3VCeuPMuJojIPvv/1MDZtL4AvoAEIXoH2655i/LqnGKmJdpyakYohJ3eAzdL6f/0KACJ0A1sNvoCGCpcPmiZCc5AkBCdmm01ycKhNUaBU39uA4YjaCklwKVMAQElJFQLVv7woPCaTjMTEWNYyQi1VR1UTKCr3NPvd2qNJUSTEx8egrMxVPSTUMtzeADZuz8fa3/JQUn0FWm0Wk4zBfTogKz0FXTrGtureohqNqaUkBYcRg1ey1cxBCl7BJsuos3q2EMHjjfYJxN+T+klKioWi6HNXsdb/pwkRURtjt5owdnBnjB7YCdv3l2DNb3nYdbAstN8X0LBxWz42bstHp+QYDO+Xgsw+HRBja9u/kkMLQyI4zAYvQleyBeciofoS/+DK2QCA6jlKcnVPE4A6V+MFH1/z7+BxRwIUh+ZIH237J4+IwibLEmJsJqgBDQFNQBMCQgveRZ30IcsS+vdMQv+eSSgsdWPdljxs2lkAt1cNHXO4yIWPftiLj9ftQ/8eSTilX0f06RIPuZ1MWq65kk2DAEIvW61zTCjsQAJqvWyp1j8khA6CLAGKJEExyaFwFZzDJIMBicLBMERkUBIAp92MmsmvqgYEVC34FdDq3gy0+sNIkgCE5opUf1jVfJIJAQ2AGtDgVwUEguFKCAF+NgEdEuyYclpPTMrqjl9+L8KGbfnYV+sqtIAq8MvvRfjl9yLExZiR2acDhpzcAZ2SY6PY6pZR81YTEDj2m6WBHT41+D6smeQtSzCbgqtv1ywyGXrfVr/Pa5+PqAbDEJGBBT8Ugp8MshScy2IxyZBsQJ0/0Y8TZ2o+WKTaf8YjeEWSqgkEVA2qKuBXNWha8K/2UDwSCH1vlA8osym4XtGwvh2RX+LGj9vzsWlnIaqqr9ICgHKXH99uPoxvNx9GWlIMMvskY3DvDkh0tr5L9KOt9iRvVRXw+rU6i0zW9Cgp1XOWTIoEudZyAZoQkCBBkSRIMkLLCzA8GQsnUFfjZLbIcWKgPtpjHWv+Mpek4IdPze0ghAh+r6oa/AENPr8GVQjdJnVHawJ1U6mahu37S7FpRwG27Ss95lBl91QHBvdOxsBeyYiLtbRoG9tKLRujJiyhOojXDM/VCU9KMCApSvWcJvlIL1PN/Kbgit5yk+Ywtcef72jhBGoialNqeqBCvUihITdAQfDKI7s1+He4P6DB71fhDQSH7Gp6k9ozRZaR0TMJGT2TUOXxY/PuIvy0owA5BVV1jtufV4n9eZVY8cM+9EhzYuBJSRjQMwnxrXBRx9as9orcwNHDc8F/BGpNa6oZZas9p+noXidZlqAocmjRSrnWlXUSqlf7rjVJXJIkCAhoWjAMA8EJ5rJ0ZDJ5e3/ftybsGarGlB45/sWjD9bxSE+Sqmnw+VX4/MGeI63W/KPG3Gy0rfdmFJa5q28IW4jCWvdDO1rXjrEY0CsYpjom2JulLW29li2pduCpHaQkCTCbZTgcdpSXuxGoeU9Xv5GlWsEp0WmN+j3uWjs9e4YYhqoZ+YNHL/wQ1wfrWF/1iEaol0irXiVZ0wS8vgD8asNDa+3lA1wIgcNFLmzeHZxg3dDaRTU6xNuQ0TMR6T0S0S3FqdutNNpLLaOtMXVUZAlJcTbeBuUEOExGRIZS8yebXH0ZW3BoLRiSYm0m+AIaXJ4AfH4VajtcRFKSJHTuEIvOHWJxdlY3HCyswq+/F+O3PcUoKq/bY1RY5sE32YfxTfZh2K0m9OuWgL7dE9C3azxibOYovQKi1o1hiIjarJqQZFZkJDisCAQ0uLx+uH1qu11ZW5IkdO3oQNeODpyd1Q25xS5s2VuCLXuLcbjIVedYtzeAn3cV4uddhZAkoGtHB/p2S8DJXePRpaODPQ9E1RiGiKhdEEJAUSTExVphtQRQ4fajvc8CkCQJnZJj0Sk5Fmec0hUlFR5s3VeKbftKsOdweZ1eMiGAA/mVOJBfic9/zIHNoqB353j07hqHPl3ikRxnaxO3BSFqDgxDRNSuCCFgMSlIcsoIaAIWswJFlgwxzyXRacNpA9Nw2sA0eHwB7Mwpw44DpdixvxQVtdYxAgCPT8Vve4vx295iAEB8rAUndY6r/ornmkZkKAxDRNQuSZBgtyhITLADagClFT54/e13+OxoNosJg05KxqCTkqEJgdwiF3bmlGLHgVLsz6usN7eqrMqHn3YW4qedhQCABIcFvTrFoVenOPRMcyIlqXmuUiNqDRiGiKjdU+TgnCKvP4AKlz94E1EDkWtNwB43pAu8fhV7DpVj98Ey7DxYhvwSd73HlFbWDUexNhP6dEtElw4x6NbRgc4dYmE26XMlD1G0MQwRkSGEhs/iZFS4/PC040nWJ2I1K0jvEbz8HgAqXD78fqgcvx8qx57D5Q2uaVTlCSB7ZwGydwa/V+RgwOqW4gh9JTqtnHdEbRLDEBEZigQJ8bFW2CwqqjzBXqLGLODYnjljLMjs0wGZfToAAMqrfNhzOBiM9uZWNNhzpGoiNCG7RozVhC4dY9GlowNdOwZ7ouJjLQxI1OoxDBGR4QR7iWRYnVb4AlrohrI+vwpVFdWrAke7ldETF1s3HLk8AeQUVCK31I3te4uRk18Fv1p/qNHlDU7a3plTFtoWazOhc4fgFW+dO8QgLTkWHeJskHlZP7UiDENEZFhCBNcoMiuAZFEAuxkBVcCvBu+LpgY0BLRgODLCPdKOJcZmQv+eiRgZ3wVlQ1zw+VXkFrmwP78SOfmV2J9fiaJj3C6kylM/IJkVGalJdqQlxyItyY7UpBikJsbAYeeikBQdDENERDiygKMiS1BkJRiOIAEIhiFfQIPPr4VWuTbqfCMgOCG9S0cHunR0AAOC29zeYO/RocIq5ORX4WBhJUorfQ0+3q9qyCmoqncj2libCSmJMUhNtCMl0Y6OiXakJNjhsJs51EbNimGIiKgBwXBUfQNNSLCaFNjMCgTM8AeCocjjVRHQNMP2GNVmt5pwctcEnNw1IbStyuPHocIqHC504VBRFQ4XVaGwzHPMelV5AqG5SrXZLAo6JtjRMcGGDvF2dIi3oUOCHclxNl7RRrpgGCIiaqTat/+wmGTE2kzw+DV4vAH4Vc3QQ2kNibWZ6wUkX0BFfrEbucUuHC5yIa/EhdxiF1yewDGfx+NT603WrhEfa0FSnA3J8TYkx1mRFGcLfjmtsFv5EUeNw3cKEVEYgqFHgs2swG4xQdU0eP0qfD4VflVAZY9RgywmBV1THOia4ghtE0Kg0u1HfokbeSVu5Je4kF/qRkGJG1XHCUlAcLHIsuqr345msyhIclqR6LQh0WlFgtOKRIcFCU4rEhxW2CwKh98IAMMQEVHEhBCQJQl2iwkxVjM0ocHr1+D1qfAHNAajE5AkCc4YC5wxFvTuEl9nX5XHj8JSDwpK3SgodaOwzIPCMjeKy731VtE+msen4lCRC4eOuoFtDatZQbzDgvhYC+IdVsTHWpDgsCAuNvgVH2uB1czAZAQMQ0REOhJCQKrVY6QJDb6ABq363mhCCPgDGvyqgBC8jP9EYm1mxKaZ0SPNWWe7pgmUVnpRVO5BUbkHxWXBfxeXe1Bc4W3UKuNev4r8EneD6yjVsJhkxMVa4IwxVwe2o/5rN8MZY4bNaoLM0NRmMQwRETWTmmBkNSmh37bBz0sJAgJq9XBaIBAMTAH1yGX8dHyyLIXmB5181L6aYbeSCm+9r9LK4FegkTfu9QW06t6ohpcOCLVHkuCwm+CwmxFrNx/5r82MWLspGOrsJjhjLLBYLRBMwK0KwxARUQuqfZVazWX8VrOCWEjQRDAQBQJasPcooIV6joQQ4Mdn49Qeduue6qy3XwiBKk+gOhj5UFbpRVmVD+VVPpRV+lDuCv77RMNwtWlCoNzlR7nL36jjZVlCjNWEGJupzn/t1mCgOnVAGrrXmldFzYthiIgoymoCkgQJZkWCWZERYwuucRRcHbu6B0kNrpR9ZJVsDrGFQ5IkOKp7b7p2bPgYTQi4PAGUVwXDUYXLj/IqHyqq/13h8qHS7UeFy9+k0BR6fi3Ye1Xpbjg8ff5jDv45cySS4mxNfm5qOoYhIqJWqPYwyqHCSrg8AThjzOiWEuzpCGgCmqbB41fx/ebDKCpzw2G3IGtACkyyEgpKmhA4XFgFlyeAGJsJnTrE1pvbcqJjVE1D9u5ieHwabBYZA3slQZEbv76PHm2I9ByqpmHzrkKUVvqQ4LBgcJ8Ox30Ncq3A1BmxxzyHhOBE7Qp3MCDl5AcXm1Q1DRIkuDwBVHn81V8BuD2BRvXw+QIadh8qQ4LTyrlILUASUR641DQNzz//PN59911UVFRgxIgRuP/++9GtW7cGjy8pKcE//vEPfPPNN5AkCVOmTMGcOXNgt9sjakdJSRUCjZhwR8dmMslITIxlLSPEOuqnrddy695irFi7D7nFLqiqgKJISEuKwZSRPdC/ZxLe/GQb1m7JgyzLMCsyTCYJJkXG8PSOmDyyJ34/VIbvf8lFfqkbPp8KQCDBacVpGWno2TkOALA7pxRfZx9GfokLqgYoMtAhwY7xmZ3Ru2sCvs0+iK9+OgSPT4VAcE1um0XB+KGdMTazywlfw+6cUnyVfQiFpe4Gn7+xx0RyjkhfQ2PO0djXIEmA2WrG4fwKVLqCIcnlCcDtDeBgQSX25lbAr2pwxlgQULU6/7+prqSkWCiKPotuRj0MPf/883jrrbfw6KOPIi0tDU888QRycnLw0UcfwWKx1Dv+qquugtvtxoMPPojy8nLcc889GDFiBB577LGI2tFWf1m2Jm39g6e1YB3105ZruXVvMV7/ZDs8vgBibWaYTDICAQ1VngBsFgVpSXZs3l3c4GMlCRjSJxmFZV74VRVxMRZYLAo0VcDr12CzKpg6sgc0Aaxcuw/egAq7WQFkCV6fiip3ACYF6J7qxKbtBQhoAjJEaNK3XwVkCTg7q9txw8TunFK8/90eeP0qYqwmKIoMVdXg8qqwmmVcMKYXAJzwmOMFohOdo3fnOPy4vQCaABQJNXdYgSoa9xoac45T+6dg3db8Rr0GRZEQHx+DsrJgwG3oHB3j7YiPsyK/yBX6/z397H4MREfRMwxFdR1zn8+HV155BbNmzcL48eORnp6OBQsWIDc3F6tXr653/E8//YT169fjsccew4ABAzBq1Cj8/e9/xwcffIC8vLwovAIiIv1pQmDF2n3w+AJIcFhhMSuQJQkWs4IEhwVVbl+dICTV+gKCc5B+2lmE4nIXLCYFXr+GiqrgMI0/oOJAXiVWrN2HVev2Iq+kCpIIDsuoqoDFrCApzgIJwO6D5Yh3WJGaYEPHpBikJjnQqaMD3VNikZoUi9/2lEBRJJhNMkxK8EtRJCiKBEkC1m3LhyxL6JhgR6zdDKtZgd1qRpLTgoAq8HX2IXy9+RC8/mBgM5uCr9NsUhAXY4bXr+Gr7EPQjvE3uyYEvso+9uM93gA2VgchkxyctCxLEmRZgkkGNAF89dMhqNqxg/KJzuH1qfjqp/BfQ0PnMJnkOv+/PT4VK9buO+5zUGSi2jO0efNmXHLJJVi1ahV69eoV2n7FFVegb9++ePDBB+sc//LLL+P111/Hd999F9rm8/mQmZmJ+fPnY/LkyWG3ReOiaBGTJECWZdYyQqyjftpqLQOqhrIqXzDgNDBfRNVEoy7NliQJilz/8bUnXktSw+fQNBH68K1ZDqChmSuOmGDIaeg1VLr9kKSGHydwZF6UdOQkdQ4QCLYz1m5qcH6Pqmqo8viDbWvgJMHXgCPnqNeG4M5YmxmWBl5Dzeuoqp7k3FCdhAhOcldk6Zj7BQCH3QxTdS+GLEt1lk8I1ar6HDX/T2qOqXmO+FhL6DkoWEe9FsSM6gTq3NxcAECnTp3qbE9JSQntqy0vL6/esRaLBQkJCTh8+HBEbZGbMBmQjo+11AfrqJ+2VsuaNXBqPhgj0fDjpVDQkRtzjuqr3cTRmxAMHA0NVQSvgBOQpWO1AaGrsBRZglTzhNKR/0qQIISAIssN3pBViGCN6kwwlmo1TgagiVCYa/ClVT/HsUKGECLUo3Sk6632fgmABkWWj1trRa57Dlk5cnDN/qMnSh/5VqoOjpJuw0JUV1TDkNsdXPXz6LlBVqsVZWVlDR7f0Dwiq9UKr9fbPI0kImphVouCFEtMtJsRkZZ4DVaLAqslsotnTnwOE6yW5v2otFlMsDXzOej4ohoxbbbg+gk+n6/Odq/X2+DVYTabrd6xNcfHxLTtXxxEREQUHVENQzVDXvn5+XW25+fnIzU1td7xaWlp9Y71+XwoLS1FSkpK8zWUiIiI2q2ohqH09HQ4HA6sW7cutK28vBxbtmzBiBEj6h0/YsQI5ObmYt++faFt69evBwCccsopzd9gIiIianeiOkhpsVgwbdo0PPnkk0hKSkKXLl3wxBNPIC0tDZMmTYKqqiguLobT6YTNZkNmZiaGDRuG2267DfPmzYPL5cL999+P888/v8GeJCIiIqITifqii6qq4qmnnsKyZcvg8XhCK1B37doVOTk5OOOMM/DII4/gwgsvBAAUFRXhwQcfxLfffgur1YpzzjkHd999N6xWazRfBhEREbVRUQ9DRERERNHEBQuIiIjI0BiGiIiIyNAYhoiIiMjQGIaIiIjI0BiGiIiIyNAYhoiIiMjQDBOGioqKcOedd2LkyJEYOnQoZs6cid27d4f2b926FdOmTcOQIUMwceJEvPHGG1FsbduxZ88eDB06FMuWLQttYy0bLy8vD/369av3VVNP1rLxli9fjsmTJ2PQoEGYMmUKPv7449C+nJwc3HDDDRg2bBjGjBmDp59+GqqqRrG1rdO6desafD/269cPZ5xxBgDWsikCgQCeeeYZTJgwAUOHDsWVV16Jn3/+ObSfP9+NU1lZiQceeABjxoxBVlYW7rjjDhQVFYX2r1mzBhdeeCEyMzNxzjnnYMWKFU0/iTCIyy67TFxyySUiOztb7Nq1S/zlL38RY8aMES6XSxQXF4tTTz1V3H333WLXrl3ivffeE4MGDRLvvfdetJvdqvl8PnHhhReKvn37iqVLlwohBGvZRF999ZUYNGiQyMvLE/n5+aEvt9vNWjbB8uXLRUZGhnjrrbfEvn37xMKFC0V6errYtGmT8Pl8YtKkSWLmzJli+/bt4tNPPxVZWVnimWeeiXazWx2v11vnfZifny9Wr14t+vXrJ9577z3WsomeffZZMXr0aPHtt9+KvXv3invuuUeccsopIi8vjz/fTXDttdeKcePGia+++krs2LFD3HTTTWLy5MnC6/WKXbt2iUGDBomnnnpK7Nq1S/z73/8WGRkZ4ocffmjSOQwRhkpLS8Xs2bPF9u3bQ9u2bt0q+vbtK7Kzs8WiRYvEmDFjhN/vD+2fP3++mDRpUjSa22bMnz9f/OlPf6oThljLpnnppZfEueee2+A+1rJxNE0TEyZMEI8++mid7ddee61YtGiR+Oijj8TAgQNFaWlpaN8777wjhg0bJrxeb0s3t02pqqoSEyZMEHfddZcQQrCWTXTeeeeJRx55JPR9RUWF6Nu3r/jkk0/4891IW7ZsEX379hVff/11aFtlZaUYPny4WLZsmbjvvvvExRdfXOcxs2fPFtdee22TzmOIYbL4+HjMnz8fffv2BQAUFxfjtddeQ1paGvr06YONGzciKysLJtORW7WNHDkSe/fuRWFhYbSa3apt2LAB//nPf/Doo4/W2c5aNs327dvRu3fvBvexlo2zZ88eHDx4EOeee26d7YsXL8YNN9yAjRs3YsCAAYiPjw/tGzlyJCorK7F169aWbm6bsmjRIrjdbsydOxcAWMsmSk5OxpdffomcnByoqor//Oc/sFgsSE9P5893I+3duxcAMHz48NC22NhY9OjRA+vXr8fGjRsxatSoOo8ZOXIkfvzxR4gm3GDDEGGotvvuuw+jRo3CihUr8PDDDyMmJga5ublIS0urc1xKSgoA4PDhw9FoZqtWXl6OOXPm4N5770WnTp3q7GMtm2bHjh0oLi7GlVdeidNOOw1XXHEFvvnmGwCsZWPt2bMHAOByuXDddddh1KhRuOSSS/DFF18AYB3DVfNH44033oiEhAQArGVT3XPPPTCbzTjjjDMwaNAgLFiwAM8++yy6d+/OWjZSQzVRVRW5ubkoLi4+Zh3dbjdKSkoafR7DhaHp06dj6dKlmDp1Km6++Wb89ttv8Hg8sFgsdY6rufGr1+uNRjNbtXnz5mHo0KH1/hIHwFo2QSAQwO+//46ysjL85S9/wUsvvYQhQ4Zg5syZWLNmDWvZSJWVlQCAuXPnYurUqXjllVcwevRo3HTTTaxjBJYsWQKn04nLLrsstI21bJpdu3bB6XTihRdewH/+8x9ceOGFuOOOO7B161bWspEGDRqEk046CQ888ADy8vLg8Xgwf/58lJSUwO/3N1jHmu99Pl+jz2M68SHtS58+fQAADz/8MLKzs/HWW2/BZrPVK1rNmzEmJqbF29iaLV++HBs3bsRHH33U4H7WsvFMJhPWrVsHRVFgs9kAAAMHDsTOnTuxePFi1rKRzGYzAOC6667DBRdcAADo378/tmzZgldffZV1DNPy5ctx/vnnh96bAH++m+Lw4cO4/fbb8dprr4WGeAYNGoRdu3bhueeeYy0byWKx4Pnnn8ecOXNw+umnw2w249xzz8WECRMgyzKsVmu9OtZ8b7fbG30eQ4Sh4uJirFmzBmeffXZofFaWZfTp0wf5+flIS0tDfn5+ncfUfJ+amtri7W3Nli5diqKiIowfP77O9gceeAArV65kLZsoNja23raTTz4Z3333HWvZSDW1qJkTWKNPnz746quvkJWVhR07dtTZxzoe37Zt23DgwIF6vb9paWmsZSNlZ2fD7/dj0KBBdbZnZmbim2++QefOnfnz3Ui9e/fG0qVLUVpaCpPJBIfDgYsvvhgjR45Ep06dGqxjTEwMnE5no89hiGGywsJCzJ49G2vWrAlt8/v92LJlC3r37o0RI0bgxx9/rLNWxtq1a9GrVy8kJydHo8mt1pNPPomVK1di+fLloS8AmDVrFh5++GHWsgl27tyJYcOGYd26dXW2//rrr+jTpw9r2UgDBgxAbGwssrOz62zfsWMHunfvjhEjRmDLli2h4TQgWMfY2Fikp6e3dHPbhI0bNyI5OblefVjLxquZx7J9+/Y623fs2IGePXvy57uRKisrMW3aNGzbtg0JCQlwOBzIycnBli1bMHr0aAwfPhzr16+v85i1a9di2LBhkOUmRBw9Ln1rC2bMmCEmTZok1q9fL7Zv3y5mz54tRowYIQ4ePCgKCwvFiBEjxNy5c8XOnTvF0qVLxaBBg8SyZcui3ew2ofal9axl46mqKi666CIxefJksWHDBrFr1y7xz3/+UwwcOFBs376dtWyCF154QQwdOlR89NFHddYZWrt2rfB4POLMM88U1113ndi6dWtobZznnnsu2s1ute6++25x9dVX19vOWjaeqqriiiuuEOecc45Ys2aN2LNnj1iwYIHo37+/+Pnnn/nz3QR//OMfxbRp08SOHTvE5s2bxdSpU8U111wjhBBix44dYsCAAeKJJ54Qu3btEosXL+Y6Q8dTXl4uHnjgATF69GgxePBgce2114odO3aE9mdnZ4tLL71UDBw4UEyYMEG8+eabUWxt21I7DAnBWjZFQUGBuOuuu8To0aPFoEGDxGWXXSY2bNgQ2s9aNt4rr7wiJk6cKAYMGCDOO+888emnn4b27d27V1xzzTVi0KBBYsyYMeLpp58WqqpGsbWt24wZM8Stt97a4D7WsvFKS0vFvHnzxPjx48XQoUPFZZddJtatWxfaz5/vxsnNzRU333yzOOWUU8SoUaPEAw88ICorK0P7v/76azF16lQxcOBAcc4554gVK1Y0+RySEE24EJ+IiIionTHEnCEiIiKiY2EYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiJqAi7NRtT+GOJGrUQUubvuugvvv//+cY/JysrCm2++2UItannvvvsudu/ejbvuuivaTSEiHXEFaiJqlP3796O4uDj0/cKFC7FlyxY8//zzoW0OhwN9+vSJRvNaxMSJE5GVlYVHH3002k0hIh2xZ4iIGqV79+7o3r176PukpCRYLBYMGTIkeo0iItIB5wwRkW42btyIadOmITMzE1lZWZg7d26d3qRly5Zh0KBB2LhxIy666CIMGjQIZ599Nr744gv8/vvvmD59OjIzM3HWWWdhxYoVdR7Xr18/ZGdn44ILLsDgwYNx7rnnYtWqVXXO7/V68fjjj2PcuHEYOHAgzj33XKxcubLOMRMnTsQ///lPTJ8+HYMHD8Y999wDANi2bRtuueUWjBw5EgMGDMDYsWPxj3/8Ax6PJ/S4gwcP4v3330e/fv2Qk5OD5557Dv369atXh379+uG5554DAOTk5KBfv3549dVXcc455yAzMxNLly4FAOzYsQM33HADhg0bhmHDhuHmm2/GgQMHdPg/QURNwTBERLrYsGEDrr76athsNjz99NP429/+hvXr1+NPf/pTKFAAQCAQwO23347LL78cL774Iux2O+644w7ceOONGD9+PBYtWoSUlBTMnTsXubm5dc5xww034IwzzsDzzz+PXr164dZbb8XXX38NIDix+eabb8Y777yDa665Bi+++CKGDh2K2267DcuXL6/zPG+//TYGDRqEhQsX4uKLL0Z+fj6uvPJKuN1uPProo3j55ZcxZcoUvPnmm3jjjTcAAM8//zw6duyIcePG4T//+Q9SUlKaVJ/nnnsO119/PR5//HGMHj0ae/bsweWXX46ioiI89thjePjhh3HgwAFcccUVKCoqCuP/ABGFi8NkRKSL+fPno1evXvjXv/4FRVEAAJmZmZgyZQqWLl2KK6+8EgCgaRpuvPFGXHLJJQCA8vJy3HbbbZg+fTquueYaAIDT6cRFF12EX3/9FWlpaaFzXHXVVbj55psBAGPHjsUFF1yAF154AePGjcMPP/yAb7/9FgsWLMDkyZNDx7jdbjz55JOYOnUqTKbgr7zOnTvjjjvuCD3vd999h/79++OZZ56Bw+EAAJx22mn4/vvvsW7dOsycORMZGRmwWCxISkoKa2jw//7v/3DRRReFvr/99ttht9vx2muvhc45atQonHnmmfj3v/+NuXPnNvkcRBQe9gwRUcTcbjeys7Mxbtw4CCEQCAQQCATQrVs39O7dG99//32d44cOHRr6d3JyMoBgcKqRkJAAIBiUarvgggtC/5YkCWeddRY2b94Mj8eDNWvWQJIkjBs3LnT+QCCAiRMnoqCgADt37gw9tn///nWed8yYMXjrrbdgtVqxa9cufP7553jxxRdRXFwMn88XWXGOcc61a9ciKysLNpst1FaHw4Hhw4fjhx9+0OWcRNQ47BkiooiVl5dD0zS8/PLLePnll+vtt1qtdb6v6QmpzW63n/A8Rw9NJScnQwiB8vJylJaWQgiBYcOGNfjY/Pz8UCCJiYmps0/TNDz11FN4++234XK50KlTJwwePLheuyNx9DlLS0uxcuXKenOagODkdCJqOQxDRBSx2NhYSJKEq6++GlOmTKm3vzFBpzFKS0vRoUOH0PeFhYVQFAUJCQlwOp2IiYkJzfE5Wo8ePY75vC+99BJee+01PPjgg5g0aRKcTicA4OKLLz5ueyRJAgCoqhoaGqyqqmrUa3E6nTjttNNCQ4O11QznEVHL4DAZEUXM4XAgIyMDv//+OwYNGhT6Ovnkk/Hcc89h3bp1upzns88+C/1bCIHVq1fjlFNOgcViQVZWFlwuF4QQddqwY8cOvPDCCwgEAsd83h9//BF9+vTBRRddFApCeXl52LFjBzRNCx0ny3V/Zdb0cNWe6P3jjz826rVkZWVh165d6N+/f6itAwcOxGuvvYZPP/20Uc9BRPpgGCIiXcyePRvfffcdbr/9dnz99df44osvMGPGDKxZswYDBgzQ5RyPP/44Xn/9dXzzzTeYNWsWdu/ejb/+9a8AgHHjxmHEiBG46aabsGTJEqxbtw4vv/wy5s2bB1mWjzv0NHjwYGzfvh0vvfQS1q9fj3fffRdXXnklfD4f3G536Li4uDhs2bIF69evh8fjwbhx4wAA999/P3744QcsXboU8+bNQ2xs7Alfy0033YT9+/fjhhtuwGeffYZvv/0Wf/nLX7BixQqkp6dHWCkiagqGISLSxZgxY7B48WLk5uZi1qxZmDNnDhRFwauvvqrbwozz5s3Df//7X9xyyy0oKCjAK6+8guHDhwMI9tq89NJLmDJlCv71r3/huuuuC11mv2DBguM+7w033IArrrgCb7zxBq6//nosXrwYf/jDH3DLLbdg586doYnc1157LQoLC3Hdddfh119/Ra9evfDYY48hJycHM2fOxBtvvIGHHnqoUZfdp6en4+2334YkSZgzZw5mzZqFgoICvPDCC5g0aVLkxSKiRuPtOIio1Vu2bBnuvvtufP755+jatWu0m0NE7Qx7hoiIiMjQGIaIiIjI0DhMRkRERIbGniEiIiIyNIYhIiIiMjSGISIiIjI0hiEiIiIyNIYhIiIiMjSGISIiIjI0hiEiIiIyNIYhIiIiMrT/D7iec7e/Eo0mAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "sns.set(color_codes=True)\n",
+ "plt.xlim(30,90)\n",
+ "plt.ylim(0,1)\n",
+ "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}