"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n",
"import matplotlib.pyplot as plt\n",
"\n",
"data[\"Frequency\"]=data.Malfunction/data.Count\n",
"data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Logistic regression\n",
"\n",
"Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature."
]
},
{
"cell_type": "code",
"execution_count": 6,
"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:
Wed, 15 Dec 2021
Deviance:
3.0144
\n",
"
\n",
"
\n",
"
Time:
10:11:01
Pearson chi2:
5.00
\n",
"
\n",
"
\n",
"
No. Iterations:
6
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: Wed, 15 Dec 2021 Deviance: 3.0144\n",
"Time: 10:11:01 Pearson chi2: 5.00\n",
"No. Iterations: 6 Covariance Type: nonrobust\n",
"===============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------\n",
"Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n",
"Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n",
"===============================================================================\n",
"\"\"\""
]
},
"execution_count": 6,
"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": 7,
"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:
Wed, 15 Dec 2021
Deviance:
18.086
\n",
"
\n",
"
\n",
"
Time:
10:11:03
Pearson chi2:
30.0
\n",
"
\n",
"
\n",
"
No. Iterations:
6
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: Wed, 15 Dec 2021 Deviance: 18.086\n",
"Time: 10:11:03 Pearson chi2: 30.0\n",
"No. Iterations: 6 Covariance Type: nonrobust\n",
"===============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------\n",
"Intercept 5.0850 3.052 1.666 0.096 -0.898 11.068\n",
"Temperature -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n",
"===============================================================================\n",
"\"\"\""
]
},
"execution_count": 7,
"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": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4VOXd//H3dyb7QmLYISA7yA5hEXEBrYK2KiriinVBpHWp7SNVn199tE+16oNt1VZxQ3GpgisupYJa44JbQBBkX8UEkJ0kkD33748ZMGAgQzLJLPm8rivXzDlzn3O+dwY+c3LmnPuYcw4REYkunlAXICIiwadwFxGJQgp3EZEopHAXEYlCCncRkSikcBcRiUI1hruZPW1mW83s28O8bmb2sJmtMbPFZjYw+GWKiMjRCGTPfTow+givnwl09f9MBKbWvSwREamLGsPdOfcxsPMITc4FnnM+XwDpZtY6WAWKiMjRiwnCOtoC31eZzvXP23xoQzObiG/vnsTExKx27drVaoOVlZV4PNHxdYH6Ep6ipS/R0g9QX/ZbtWrVdudc85raBSPcrZp51Y5p4Jx7AngCYNCgQW7+/Pm12mB2djYjRoyo1bLhRn0JT9HSl2jpB6gv+5nZd4G0C8bHYC5QdRc8E9gUhPWKiEgtBSPc3wKu8J81czywxzn3k0MyIiLScGo8LGNmLwEjgGZmlgvcCcQCOOceA2YDZwFrgH3AVfVVrIiIBKbGcHfOXVLD6w64PmgViUhEKCsrIzc3l+Li4gbZXlpaGsuXL2+QbdW3QPqSkJBAZmYmsbGxtdpGML5QFZFGKDc3l9TUVDp06IBZdedVBFdBQQGpqan1vp2GUFNfnHPs2LGD3NxcOnbsWKttRMd5RSLS4IqLi2natGmDBHtjY2Y0bdq0Tn8VKdxFpNYU7PWnrr9bhbuISBTSMXcRiVher5c+ffocmJ41axYdOnQIXUFhROEuIhErMTGRRYsWHfb18vJyYmIaZ8zpsIyIRJXp06dz4YUXcvbZZ3PGGWcAMGXKFAYPHkzfvn258847D7S955576N69Oz/72c+45JJLeOCBBwAYMWIE+4dH2b59+4G/BioqKpg8efKBdT3++OPAj8MJjB07lh49enDZZZfhO0sccnJyOOGEE+jXrx9DhgyhoKCAUaNGHfShNHz4cBYvXhzU30Pj/EgTkaD649tLWbYpP6jr7NmmCXee3euIbYqKiujfvz8AHTt25I033gDg888/Z/HixWRkZDB37lxWr17NV199hXOOc845h48//pjk5GRmzJjBwoULKS8vZ+DAgWRlZR1xe9OmTSMtLY2cnBxKSkoYPnz4gQ+QhQsXsnTpUtq0acPw4cOZN28eQ4YM4aKLLmLmzJkMHjyY/Px8EhMTueKKK5g+fToPPvggq1atoqSkhL59+wbht/YjhbuIRKzDHZY5/fTTycjIAGDu3LnMnTuXAQMGAFBYWMjq1aspKCjgvPPOIykpCYBzzjmnxu3NnTuXxYsX8+qrrwKwZ88eVq9eTVxcHEOGDCEzMxOA/v37s2HDBtLS0mjdujWDBw8GoEmTJgCcd955DB8+nClTpvD0009z5ZVX1u0XUQ2Fu4jUWU172A0tOTn5wHPnHLfffjvXXXfdQW0efPDBw55uGBMTQ2VlJcBB55o75/j73//OqFGjDmqfnZ1NfHz8gWmv10t5eTnOuWq3kZSUxOmnn86bb77Jyy+/TG1HyD0SHXMXkag2atQonn76aQoLCwHIy8tj69atnHzyybzxxhsUFRVRUFDA22+/fWCZDh06sGDBAoADe+n71zV16lTKysoAWLVqFXv37j3stnv06MGmTZvIyckBfFemlpeXAzBhwgRuuukmBg8efOCvjGDSnruIRLUzzjiD5cuXM2zYMABSUlJ44YUXGDhwIBdddBH9+/fn2GOP5aSTTjqwzC233MK4ceN4/vnnOfXUUw/MnzBhAhs2bGDgwIE452jevDmzZs067Lbj4uKYOXMmN954I0VFRSQmJvL+++8DkJWVRZMmTbjqqnoaa9E5F5KfrKwsV1sffvhhrZcNN+pLeIqWvtRnP5YtW1Zv665Ofn5+va7/zjvvdFOmTKnXbeyXn5/v8vLyXNeuXV1FRcVh21X3OwbmuwAyVodlREQa2IsvvsjQoUO555576u3WgTosIyIC3HXXXQ22rUsvvfQnX/AGm/bcRaTWnKv2dskSBHX93SrcRaRWEhIS2LFjhwK+Hjj/eO4JCQm1XocOy4hIrWRmZpKbm8u2bdsaZHvFxcV1CrtwEkhf9t+JqbYU7iJSK7GxsbW+S1BtZGdnH7jKNNI1RF90WEZEJAop3EVEopDCXUQkCincRUSikMJdRCQKKdxFRKKQwl1EJAop3EVEopDCXUQkCincRUSiUMSF+77Sct7bUEZ5RWWoSxERCVsRF+7vLN7MP1eUMu7xz/lux+HvXSgi0phFXLiPG9SOSX3jWbO1kDMf+oSZORs15KiIyCEiLtwBjm8Tw7s3n0z/dunc+toSbnhxIXuKykJdlohI2IjIcAdok57IC9cM5dbRPZizdAtnPfQJX2/cFeqyRETCQsSGO4DHY/xqRGde/dUJeDww7rHPefLjdTpMIyKNXkDhbmajzWylma0xs9uqeT3NzN42s2/MbKmZXRX8Ug+vf7t03rnxJE47rgX3zF7OxOcX6DCNiDRqNYa7mXmBR4AzgZ7AJWbW85Bm1wPLnHP9gBHAX8wsLsi1HlFaYiyPXZ7F//yiJx+u2Mq5//iUFVvyG7IEEZGwEcie+xBgjXNunXOuFJgBnHtIGwekmpkBKcBOoDyolQbAzLj6xI7MmHg8+0orGPPIPN76ZlNDlyEiEnJW0/FpMxsLjHbOTfBPjweGOuduqNImFXgL6AGkAhc55/5VzbomAhMBWrZsmTVjxoxaFV1YWEhKSsoR2+wuqeTRRSWs2lXJWR1jGdstFo9ZrbZXnwLpS6RQX8JPtPQD1Jf9Ro4cucA5N6jGhs65I/4AFwJPVZkeD/z9kDZjgb8BBnQB1gNNjrTerKwsV1sffvhhQO1Kyirc/3tjsTv21nfcFdO+dHuKSmu9zfoSaF8igfoSfqKlH86pL/sB810Nue2cC+iwTC7Qrsp0JnDosY6rgNf9217jD/ceAay7XsXFeLh7TB/uPb8P89Zs5/xHP9NVrSLSKAQS7jlAVzPr6P+S9GJ8h2Cq2gicBmBmLYHuwLpgFloXlwxpz/PXDGV7YQnnPjKPL9btCHVJIiL1qsZwd86VAzcAc4DlwMvOuaVmNsnMJvmb/Qk4wcyWAB8AtzrnttdX0bUxrHNT3rx+OE2T4xg/7UveWJgb6pJEROpNTCCNnHOzgdmHzHusyvNNwBnBLS34jm2azOu/Gs6kFxbw25nfsHFHETed1gULwy9aRUTqIqKvUK2NtKRYnr16COcPbMvf3l/Fba8t0fDBIhJ1AtpzjzZxMR7+cmE/MtMTefg/a9hWWMI/Lh1AUlyj/HWISBRqdHvu+5kZvzujO/ec15vslVu59Mkv2bW3NNRliYgERaMN9/0uG3osUy/PYtnmfC58/HM27S4KdUkiInXW6MMdYFSvVjx71RC27Clm7NTPWLetMNQliYjUicLdb1jnpsyYeDzF5ZWMe/xzlm/WoGMiErkU7lX0bpvGy9cNI8bj4aLHP2ehbv4hIhFK4X6ILi1SeGXSMNKT4hg/7StyNuwMdUkiIkdN4V6NdhlJvHzdMFo0ieeKaV/x2ZqwuthWRKRGCvfDaJWWwMyJw2ifkcRV03P4eNW2UJckIhIwhfsRNE+N56WJx9OxWTITnpuvgBeRiKFwr0FGchwvXns8nZunMOG5+XykgBeRCKBwD0BGchwvThhK5+YpTHxuPvN0DF5EwpzCPUDHJMfxzwlD6dA0mWuezdGY8CIS1hTuRyEjOY5/XjuUzGOSuHp6Dgu+02mSIhKeFO5HqVlKPC9eO5SWTRK48ukcvs3bE+qSRER+QuFeCy1SE/jnhKE0SYxl/LQvWbmlINQliYgcROFeS23SE3nx2qHExXi4fNqXuvG2iIQVhXsdHNs0mReuGUp5RSWXPfUlW/YUh7okERFA4V5nXVum8uzVQ9i9r4zLp+mGHyISHhTuQdA3M50nrxjExp37uHJ6DntLykNdkog0cgr3IBnWuSn/uGQAS3J3M+mFBZSW66bbIhI6CvcgOqNXK+47vy+frN7Of73yDZWVLtQliUgjFRPqAqLNuMHt2LG3lPvfXUHzlHju+MVxmFmoyxKRRkbhXg8mndKJrQXFPD1vPS2bxHPdKZ1DXZKINDIK93pgZtzx855sKyjh3n+voEWTeM4bkBnqskSkEVG41xOPx/jLuH7sKCzl968upkVqAsO7NAt1WSLSSOgL1XoUH+PlsfFZdGqWwqTnF7B8c36oSxKRRkLhXs/SEmN55qrBJMfHcNUzOWzeUxTqkkSkEVC4N4A26Yk8c9VgCkvKueqZHAqKy0JdkohEOYV7AzmudRMevWwgq7cWcv2LCymr0EVOIlJ/FO4N6ORuzfnzeb35eNU2/ufNpTini5xEpH7obJkGdtHg9ny3Yx+PZq+lY7MkuoW6IBGJStpzD4FbzujOz/u05t5/r2D+Fg0yJiLBF1C4m9loM1tpZmvM7LbDtBlhZovMbKmZfRTcMqPL/nPg+7dL54nFJSzO3R3qkkQkytQY7mbmBR4BzgR6ApeYWc9D2qQDjwLnOOd6ARfWQ61RJSHWyxPjB5EaZ0x4dr5OkRSRoApkz30IsMY5t845VwrMAM49pM2lwOvOuY0AzrmtwS0zOjVPjee3WQnsK63gmunzNQ68iASN1XTGhpmNBUY75yb4p8cDQ51zN1Rp8yAQC/QCUoGHnHPPVbOuicBEgJYtW2bNmDGjVkUXFhaSkpJSq2XDTWFhIeuKEvjbghL6t/By44B4PBE6imS0vS/R0Jdo6QeoL/uNHDlygXNuUE3tAjlbprqkOfQTIQbIAk4DEoHPzewL59yqgxZy7gngCYBBgwa5ESNGBLD5n8rOzqa2y4ab7OxsbvrFCJq0Wc9dby8jp6Q1t47uEeqyaiXa3pdo6Eu09APUl6MVSLjnAu2qTGcCm6pps905txfYa2YfA/2AVUhAfnlCB1ZvLWRq9lq6NE/hgiyNIikitRfIMfccoKuZdTSzOOBi4K1D2rwJnGRmMWaWBAwFlge31OhmZtx1Ti9O6NyU219fwoLvdoa6JBGJYDWGu3OuHLgBmIMvsF92zi01s0lmNsnfZjnwLrAY+Ap4yjn3bf2VHZ1ivR4evWwgrdMTuO75BeTt1hk0IlI7AZ3n7pyb7Zzr5pzr7Jy7xz/vMefcY1XaTHHO9XTO9XbOPVhfBUe79KQ4pv1yECVllUx4VmfQiEjt6ArVMNSlRSoPXzqAlVvy+a+XdaNtETl6CvcwNbJ7C/77rON4d+kWHv7P6lCXIyIRRgOHhbFrTuzIii0FPPj+arq3TOXMPq1DXZKIRAjtuYcxM+Oe83ozsH06v3v5G5Zt0m36RCQwCvcwt/8+rGmJsVz73Hx2FJaEuiQRiQAK9wjQIjWBJ67IYnthCb/+59e6i5OI1EjhHiH6ZqZz/wV9+XL9Tv749tJQlyMiYU5fqEaQMQPasnxzPo9/vI6erdO4dGj7UJckImFKe+4R5veje3BKt+bc+da3zN+gIQpEpHoK9wjj9RgPXzyAzGOSmPTCAjZpiAIRqYbCPQKlJcXy5BVZFJdVct3zCyguqwh1SSISZhTuEapLi1T+dlF/luTt4fbXl1DTTVdEpHFRuEew03u25Hend+ONhXlM+3R9qMsRkTCicI9wN4zswuherfjz7OV8snpbqMsRkTChcI9wHo/xwLh+dGmRwg0vLmTjjn2hLklEwoDCPQqkxMfw5BW+++Ve+5zGgBcRhXvUOLZpMv+4dACrtxZwyyvf6AtWkUZO4R5FTuranNvO7MG/v93CIx+uCXU5IhJCCvcoc+1JnRjTvw1/eW8V7y/7IdTliEiIKNyjjJlx3wV96dWmCTfPXMSarYWhLklEQkDhHoUSYr08Pn4Q8TEeJj43nz1FZaEuSUQamMI9SrVNT2Tq5Vls3LmP38xYSIVusi3SqCjco9iQjhnceU4vsldu44G5K0Ndjog0II3nHuUuH9qeZZvymZq9lp6tm3B2vzahLklEGoD23KOcmfHHc3ox6NhjmPzqN3ybtyfUJYlIA1C4NwJxMR6mXp7FMUlxTHxuPtt1k22RqKdwbySap8bz5BWD2LG3lF+9sIDSct1kWySaKdwbkd5t0/i/sX3J2bCLO99aqiEKRKKYvlBtZM7t35blmwt47KO19GydyvhhHUJdkojUA+25N0KTR3Xn1B4t+OPby/hs7fZQlyMi9UDh3gh5PcZDF/enQ7Nkrv/n1xoDXiQKKdwbqdSEWJ66YhCVDiY8l0NBsYYoEIkmCvdGrEOzZB69bCBrt+3l5hmLNESBSBRRuDdyw7s0486ze/LBiq1MmaMhCkSihc6WEcYffywrt/jOoOnWMoXzB2aGuiQRqaOA9tzNbLSZrTSzNWZ22xHaDTazCjMbG7wSpb6ZGXed04thnZpy22tLWPDdzlCXJCJ1VGO4m5kXeAQ4E+gJXGJmPQ/T7n5gTrCLlPoX6/Uw9fKBtElPYOJzC8jdpTNoRCJZIHvuQ4A1zrl1zrlSYAZwbjXtbgReA7YGsT5pQOlJcUy7cjClFZVMeHa+zqARiWBW0yXo/kMso51zE/zT44GhzrkbqrRpC7wInApMA95xzr1azbomAhMBWrZsmTVjxoxaFV1YWEhKSkqtlg034diXpdsr+MuCYno383LzwHg8ZgEtF459qa1o6Uu09APUl/1Gjhy5wDk3qKZ2gXyhWt3/7EM/ER4EbnXOVdgRgsA59wTwBMCgQYPciBEjAtj8T2VnZ1PbZcNNOPZlBJCW+R1/mPUtnxS24M6zewW0XDj2pbaipS/R0g9QX45WIOGeC7SrMp0JbDqkzSBghj/YmwFnmVm5c25WUKqUBnf58ceybttenp63nk7NkjUGjUiECSTcc4CuZtYRyAMuBi6t2sA513H/czObju+wjII9wv2/nx/Hdzv2cudbS8nMSGJk9xahLklEAlTjF6rOuXLgBnxnwSwHXnbOLTWzSWY2qb4LlNDxeoyHLxnAca2bcMM/v2bZpvxQlyQiAQroPHfn3GznXDfnXGfn3D3+eY855x6rpu2V1X2ZKpEpOT6Gab8cTGpCLNc8m8OWPcWhLklEAqDhB6RGrdISmHblIPKLyrh6eg6FJeWhLklEaqBwl4D0apPGI5cNZOUPBdzw4teUV+g2fSLhTOEuARvRvQV/Orc32Su3cceb3+o2fSJhTAOHyVG5dGh7cnft49HstbRNT+SGU7uGuiQRqYbCXY7a5FHd2bKnmAfmrqJVWiJjszSKpEi4UbjLUTMz7rugL1sLSrjttcU0T43nlG7Na7WuWQvzmDJnJZt2F9EmPZHJo7ozZkDbIFcs9UXvX/jSMXeplbgY3yiS3Vqm8qsXFvDN97uPeh2zFuZx++tLyNtdhAPydhdx++tLmLUwL/gFS9Dp/QtvCneptdSEWKZfPZimKXFcPT2HLXuP7gyaKXNWUlRWcdC8orIK3REqQuj9C28Kd6mTFqkJPHvVEBzwwPxifsgP/CKnTbuLjmq+hBe9f+FN4S511ql5CtOvGkxhqeOKaV+xZ19g48C3SU88qvkSXvT+hTeFuwRF38x0bhqYwPrte7n62Rz2ldZ8FevkUd1JjPUeNC8x1svkUd3rq0wJIr1/4U3hLkHTs6mXBy/uz8KNu5j0wteUlh/5GPyYAW259/w+tE1PxIC26Ynce34fnW0RIfT+hTedCilBdVaf1tx3fl9+/9pifjtzEQ9fMgCv5/A3cBkzoK3CIILp/QtfCncJunGD25FfXMbd/1pOcryX+87vi+cIAS8iwadwl3ox4aROFBSX89AHq0mM9XLXOb040i0YRSS4FO5Sb27+WVf2lZbz5CfrSYjzctvoHgp4kQaicJd6Y2b891nHsa+0gsc/Wkd8jJffnd4t1GWJNAoKd6lXZsafzu1NaXklD3+wmliPceNpGklSpL4p3KXeeTy+gcYqKh1/eW8VMV4PvxrROdRliUQ1hbs0CK/HmHJhP8orHfe/uwKH49cjuoS6LJGopXCXBuP1GH8d1w8z+L93V+IcXD9SAS9SHxTu0qBivB7+Oq4/hm9UwfIKx02nddFZNCJBpnCXBuf1GH8Z1x+vx8Pf3l9FSXkFk0d1V8CLBJHCXULC6zGmjO1LXIyHR7PXUlxWyR2/OE4BLxIkCncJGY/H+PN5vYmP8fD0vPXsKy3nnvP6HHEsGhEJjMJdQsrMuPPsnqQmxPD3/6yhsKScv47rT1yMBiwVqQuFu4ScmfFfZ3QnJT6Ge/+9goLicqZePpCkOP3zFKkt7R5J2LjulM7cd34fPlm9jUuf/JJde0tDXZJIxFK4S1i5eEh7Hr0si2Wb87nw8c/J0/04RWpF4S5hZ3TvVjx71RB+yC/mvEfmsWxTfqhLEok4CncJS8M6N+WVScPwmDHu8c/5ZPW2UJckElEU7hK2erRqwhvXn0DmMYlc9UwOM3M2hrokkYihcJew1jotkVcmDWNY56bc+toS7n93BZWVLtRliYQ9hbuEvdSEWJ6+cjCXDGnP1Oy1/OqfC9hXWh7qskTCWkDhbmajzWylma0xs9uqef0yM1vs//nMzPoFv1RpzGK9Hv58Xm/u+EVP3lv2A2Onfs4mnUkjclg1hruZeYFHgDOBnsAlZtbzkGbrgVOcc32BPwFPBLtQETPjmhM7Mu3KwXy/cx9n//1Tvlq/M9RliYSlQPbchwBrnHPrnHOlwAzg3KoNnHOfOed2+Se/ADKDW6bIj0Z2b8Eb1w8nLTGWS5/8guc/34BzOg4vUpXV9J/CzMYCo51zE/zT44GhzrkbDtP+FqDH/vaHvDYRmAjQsmXLrBkzZtSq6MLCQlJSUmq1bLhRX2pvX5nj8cUlfLOtguFtYriiVxzx3uAMOhYt70u09APUl/1Gjhy5wDk3qMaGzrkj/gAXAk9VmR4P/P0wbUcCy4GmNa03KyvL1daHH35Y62XDjfpSNxUVle5v7610HW57x43620du/bbCoKw3Wt6XaOmHc+rLfsB8V0O+OucCOiyTC7SrMp0JbDq0kZn1BZ4CznXO7QhgvSJ15vEYN/+sG89cOZjNe4o5+++f8q/Fm0NdlkjIBRLuOUBXM+toZnHAxcBbVRuYWXvgdWC8c25V8MsUObIR3Vvwr5tOpHOLFK5/8WvumPUtxWUVoS5LJGRqDHfnXDlwAzAH3yGXl51zS81skplN8jf7H6Ap8KiZLTKz+fVWschhZB6TxMvXDePakzry/BffMeaReaz+oSDUZYmEREADZjvnZgOzD5n3WJXnE4CffIEq0tBmL9nM7CVbAFj1QwFnPfwJY/q3Zd6a7WzeU0yb9EQmj+rOmAFtg77tWQvzmDJnJZt2F9XrdgLxh1lLeOnL77m5dxnX3D6bS4a24+4xfUJSi4SG7oYgUWPWwjxuf30JRf7DMZUOXIXjlQW5B9rk7S7i9teXAAQ1eA/ddn1tJxB/mLWEF774cRyeCucOTCvgGw8NPyBRY8qclQfCdb/qTvQtKqtgypyV9b7t+thOIF768vujmi/RSeEuUeNohiMI9k1ADrftUAyRUHGYa1cON1+ik8Jdokab9MSA23oM3lyUF7QrWw+37aOpKVi8Vv2FXIebL9FJ4S5RY/Ko7iTGeg+aF+sxYg+5ajU+xkPmMUn8ZsYiJjw7n9xd++pl24mxXiaP6l7ndR+tS4a2O6r5Ep0U7hI1xgxoy73n96FteiIGtE1PZMqF/Zgytt9B8+6/oC8f3jKCP/z8OD5bu4PT//oxT3y8lrKKyqBu+97z+4TkbJm7x/Th8uPbH9hT95px+fHt9WVqI6OzZSSqjBnQttpArW7ehJM6Mbp3K+56ayl/nr2C1xbk8cdze3F8p6ZB3XYo3D2mD3eP6UN2djZrLxsR6nIkBLTnLo1a5jFJPPXLwTwxPovCknIufuILbnppITuLa78XLxIOtOcuApzRqxUndW3O1I/W8thHa3m3spLvvKu47pROJMXpv4lEHu25i/glxnn53end+OB3p9C/hZeHPljNyAeymZmzkQrdt1UijMJd5BDtMpL4df8EXpk0jNZpidz62hLOfOhj3lv2g24KIhFD4S5yGIM7ZPDGr0/g0csGUlbhuPa5+Yx59DM+Xb1dIS9hT+EucgRmxll9WvPeb0/m/gv6sC2/mMunfcm4xz9n3hqFvIQvhbtIAGK8Hi4a3J7/3DKC/z23F9/vLOKyp77kgqmf8Z8VOlwj4UfhLnIUEmK9XDGsA9mTR/Cnc3vxQ34JV0+fz1kPf8obC3PrdCGUSDAp3EVqISHWy3h/yD9wYT/KKir57cxvOPn/PuTxj9ayZ19ZqEuURk4n8IrUQazXw9isTM4f0JbsVVt54uN13PvvFTz4/mouyGrLFcM60K1laqjLlEZI4S4SBB6PcWqPlpzaoyVLN+1h+rwNvDw/lxe+2MjQjhlcfvyxnNGrJfEx3ppXJhIECneRIOvVJo0pF/bj9rOO4+X53/PCF99x40sLyUiO44KBbRk3qB1dtTcv9UzhLlJPMpLjmHRKZyae1IlP12znpa828sy8DTz5yXr6ZaYxNiuTn/dtQ0ZyXKhLlSikcBepZx6PcXK35pzcrTnbC0uYtTCPVxfkcsebS/nj28sY0b05Z/drw8+Oa0lyvP5LSnDoX5JIA2qWEs+Ekzox4aROLN+cz6yFeby5aBPvL99KQqyHU3u04MzerRnZowUpCnqpA/3rEQmR41o34bjWTbh1dA/mf7eLt7/ZxL+/3cLsJVuIi/FwUpdmnN6zJacd15LmqfGhLlcijMJdJMQ8HmNIxwyGdMzgrnN6seC7Xcxespn3lv3AByu2YraEvpnpnNajBSO6N6d3mzQ8Ht0PVY5M4S4SRrxVgv7Os3uyfHMBHyz3hfzf3l/FX99bRUZyHCd2acaJXZtxYpc9/m12AAANA0lEQVRmIbkJt4Q/hbtImDIzerZpQs82TbjxtK5sLyzh09XbyV65lU/XbOetbzYB0LFZMsd3yuD4Tk0Z0jGD1mkKe1G4i0SMZinxB+7T6pxj5Q8FfLp6O5+v3cE732zmpa++B6BdRiKDj81g4LHH4Aoqqah0eHUYp9FRuItEIDOjR6sm9GjVhAkndaK8opLlmwv4asNOvlq/g49Xb+P1hXkA3Jczh76Z6fRrl06/zDT6tkunTVoCZgr8aKZwF4kCMV4PfTLT6JOZxjUndsQ5x8ad+3jh3c8oSWnNwo27mfbpOsoqfEMTH5MUS++2afRs7Tvsc1zrJnRslkysV2MJRguFu0gUMjOObZrM8LaxjBjRG4CS8gpWbC5gce5ulm7K59tNe3hm3gZK/cMUx3k9dGqeTPdWqXRtkULXlql0aZFC+4wkhX4EUriLNBLxMV7foZl26QfmlVVUsm7bXpZvzmfFlgJWbMln/oZdvLlo04E2sV6jfUYSHZul0Kl5Msc2TaJj02TaN02idVqijueHKYW7SCMW6/XQvVUq3VsdPJBZYUk5a7YWsnZrIWu2FbJ+217Wb9/Lx6u3UVpeWWV5o216Iu0yksg8JonMYxJpm55I22MSaZOeSMvUeGK01x8SCncR+YmU+Bj6t0unf5W9fIDKSseW/GI2bN/Lhh37+H7XPjbu2Efurn3M3bSFHXtLD2rvMWiRmkCrtARaNfE9tmgST4vUBFqkxtOiSTzNUuLJSIrThVlBpnAXkYB5PEabdN9e+Qldfvr6vtJyNu0uIndXEZv3FLN5dxGb9hTzQ34xa7YVMm/tdgqKy3+ynNdjZCTH0TQ5jmYp8TRNiSMjOY6MpDgyUuI4JimO73ZU0HJzPsckxZGWGEtCrEdn/ByBwl1EgiYpLoYuLVLp0uLw49XvKy1na34JWwtK2FpQzPaCErYXlrKtoIQde0vZXljCxp372LW3lIKSgz8I7s/55MDzOK+HtKRYmiTEkJYYS5PEWJokxJKaEEPqgccYUuJ//En2//iee0mKi4na7wwCCnczGw08BHiBp5xz9x3yuvlfPwvYB1zpnPs6yLWKRK1ZC/OYMmclm3YX0SY9kcmjuvPK/I3MW7vzQJvhnTO4cFD7n7QDfjJv/nc7eenL77m5dxnX3D6bS4a24+4xfQLa7pgBbQ87P5Dl92+7wjm8Zj/ZdlJcDB2axbDo+9019uXOs3tycrfm7NxXyofzcujQrSe79pXx2drtZK/cxraCEgqKy/B6jLIKx4bte8kvLqeguOzAaZ81iY/xkBTnC/qkOC+JcV4SY70kxXlJiPU9T4jzkhDjJSHWQ0Lsj4/xMR7iY/yPsT8+j9v/4/3xebzXS2yMEev14FxgtdVFjeFuZl7gEeB0IBfIMbO3nHPLqjQ7E+jq/xkKTPU/ikgNZi3M4/bXl1BUVgFA3u4ibp656Cft5q3deVDY5+0uYvKr34CDskp3YN7vZi6isspyFc7xwhcbAQ4K2eq2e/vrS5j/3U5eW5D3k/nAQQFf3fJ12fbkV74B40Ao5+0u4o43l3Lv+X0YM6AtW5p6GdGnNbMW5vHB8q0Hli0uq+T7nUUH2u1XXFZBYUk5hcXlFBSXU1hSzt6ScvaWlrO3pIK9JeXsK61gb2k5+0p9z4tKKw487thbemC6pLyC4rJKisoqqKisezCf2TGWkSPrvJojCmTPfQiwxjm3DsDMZgDnAlXD/VzgOef7OPrCzNLNrLVzbnPQKxaJMlPmrDwQVEerur3TymraAbz05fcHBWx12y0qqziw133o/ClzVh4UntUtX5dtl1UTmoFut7p2vj1sL81SgjtccllFJcVlFZSUV1JS7nte6n/ue6ygpKyS0grf9IHH8krKKnw/nl0bg1pTdaymPw/MbCww2jk3wT89HhjqnLuhSpt3gPucc5/6pz8AbnXOzT9kXROBif7J7sDKWtbdDNhey2XDjfoSnhqsL3GtumTV17or9u3Bm5R2YLp0y5oFddluXZYPwrLNgO1HWrbqNsJcXf59Heuca15To0D23Kv7tuHQT4RA2uCcewJ4IoBtHrkgs/nOuUF1XU84UF/CU7T0xczml+/ZGvH9gOh5T6Bh+hLI1QW5QLsq05nAplq0ERGRBhJIuOcAXc2so5nFARcDbx3S5i3gCvM5Htij4+0iIqFT42EZ51y5md0AzMF3KuTTzrmlZjbJ//pjwGx8p0GuwXcq5FX1VzIQhEM7YUR9CU/R0pdo6QeoL0elxi9URUQk8mhEHxGRKKRwFxGJQmEf7maWYGZfmdk3ZrbUzP7on59hZu+Z2Wr/4zGhrjUQZuY1s4X+awMiuR8bzGyJmS0ys/n+eZHal3Qze9XMVpjZcjMbFol9MbPu/vdj/0++md0coX35rf//+7dm9pI/ByKuHwBm9ht/P5aa2c3+efXel7APd6AEONU51w/oD4z2n5FzG/CBc64r8IF/OhL8BlheZTpS+wEw0jnXv8r5upHal4eAd51zPYB++N6fiOuLc26l//3oD2ThO7nhDSKsL2bWFrgJGOSc643vRI6LibB+AJhZb+BafFf69wN+YWZdaYi+OOci5gdIAr7GN27NSqC1f35rYGWo6wug/kz/G3kq8I5/XsT1w1/rBqDZIfMiri9AE2A9/pMLIrkvh9R/BjAvEvsCtAW+BzLwndH3jr8/EdUPf50X4htscf/0HcDvG6IvkbDnvv9QxiJgK/Cec+5LoKXzn0vvf2wRyhoD9CC+N7bqEByR2A/wXYE818wW+IeVgMjsSydgG/CM/3DZU2aWTGT2paqLgZf8zyOqL865POABYCOwGd91M3OJsH74fQucbGZNzSwJ3ynj7WiAvkREuDvnKpzvT81MYIj/T52IYma/ALY65yJl7IuaDHfODcQ3Iuj1ZnZyqAuqpRhgIDDVOTcA2EsE/Ll/JP6LDc8BXgl1LbXhP/58LtARaAMkm9nloa2qdpxzy4H7gfeAd4FvgJ/eraQeRES47+ec2w1kA6OBH8ysNYD/cWsISwvEcOAcM9sAzABONbMXiLx+AOCc2+R/3IrvuO4QIrMvuUCu/69BgFfxhX0k9mW/M4GvnXM/+KcjrS8/A9Y757Y558qA14ETiLx+AOCcm+acG+icOxnYCaymAfoS9uFuZs3NLN3/PBHfG78C35AHv/Q3+yXwZmgqDIxz7nbnXKZzrgO+P5n/45y7nAjrB4CZJZtZ6v7n+I6HfksE9sU5twX43sy6+2edhm8464jrSxWX8OMhGYi8vmwEjjezJDMzfO/JciKvHwCYWQv/Y3vgfHzvTb33JeyvUDWzvsCz+L4x9wAvO+f+18yaAi8D7fH9Y7jQObfz8GsKH2Y2ArjFOfeLSOyHmXXCt7cOvsMaLzrn7onEvgCYWX/gKSAOWIdv+AwPkdmXJHxfRnZyzu3xz4u498V/yvNF+A5hLAQmAClEWD8AzOwToClQBvzOOfdBQ7wnYR/uIiJy9ML+sIyIiBw9hbuISBRSuIuIRCGFu4hIFFK4i4hEoUBukC3SoPyniX3gn2wFVOAbIgBgiHOuNCSFHYGZXQ3M9p83LxJyOhVSwpqZ3QUUOuceCINavM65isO89ilwg3Nu0VGsL8Y51yCXokvjo8MyElHM7JfmG99/kZk9amYeM4sxs91mNsXMvjazOWY21Mw+MrN1ZnaWf9kJZvaG//WVZvaHANd7t5l9hW9coz+aWY5/fO7HzOcifMNRz/QvH2dmuVWurD7ezN73P7/bzB43s/fwDVYWY2Z/9W97sZlNaPjfqkQjhbtEDP+AcecBJ/gHkovBN5QDQBow1z+YWSlwF77L1i8E/rfKaob4lxkIXGpm/QNY79fOuSHOuc+Bh5xzg4E+/tdGO+dmAouAi5xvPPWaDhsNAM52zo0HJuIbUG4IMBjfIGzta/P7EalKx9wlkvwMXwDO9w05QiK+S+0Bipxz7/mfL8E3TGy5mS0BOlRZxxzn3C4AM5sFnIjv/8Hh1lvKj0MtAJxmZpOBBKAZsAD491H2403nXLH/+RnAcWZW9cOkK75L0kVqTeEukcSAp51zdxw00ywGXwjvV4nvDl77n1f9d37ol0yuhvUWOf8XU/5xW/4BDHTO5ZnZ3fhCvjrl/PiX8aFt9h7Sp1875z5AJIh0WEYiyfvAODNrBr6zampxCOMM890zNQnfmOHzjmK9ifg+LLb7R8W8oMprBUBqlekN+G51xyHtDjUH+LX/g2T/fVATj7JPIj+hPXeJGM65Jf7RAt83Mw++UfYmAZuOYjWfAi8CnYHn95/dEsh6nXM7zOxZfMMbfwd8WeXlZ4CnzKwI33H9u4AnzWwL8NUR6nkc38iAi/yHhLbi+9ARqROdCimNhv9MlN7OuZtDXYtIfdNhGRGRKKQ9dxGRKKQ9dxGRKKRwFxGJQgp3EZEopHAXEYlCCncRkSj0/wHRUJwHFwSFegAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n",
"data_pred['Frequency'] = logmodel.predict(data_pred)\n",
"data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n",
"plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"scrolled": true
},
"source": [
"This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing and plotting uncertainty"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9+P/XmS3LZA9ZgISwGHZEC6KoiKKIsitSBapWxaqt9tta71Vvrbfaq9b767W3tr2taK1raQUXlihaQUFUFOsSCYQ9EJZMQtbJ7HPO+f0xyUBIApOQySx5Px+PPDLnzJmTz4eQec9ne38UXdd1hBBCiJMYIl0AIYQQ0UkChBBCiA5JgBBCCNEhCRBCCCE6JAFCCCFEhyRACCGE6FDYAsSDDz7I5MmTmT17dofP67rOf/3XfzF9+nTmzJlDWVlZuIoihBCiG8IWIK699lqee+65Tp/ftGkTFRUVvPfee/zqV7/il7/8ZbiKIoQQohvCFiDOO+880tPTO31+/fr1zJ8/H0VROOecc2hqaqK6ujpcxRFCCNFFERuDsNls5OfnB4/z8/Ox2WynfZ0s/BZCiN5hitQP7uiNXlGU076uodlDTU0zigJKy2sUJfDdEPx+wjkDLcenv3c0yMlJpabGHulihI3UL3bFc92gb9SvqyIWIPLz86mqqgoeV1VVkZube9rX+XwaHp/a5Z+nAIohEESMBgMGBQwGJfClKBhbHhsNsRNMhBAinCIWIKZNm8Yrr7zCrFmz+Oabb0hNTQ0pQHSXDuiajgb41VMHGINBwagoGI2BgGE0GDAaFUzGwGMhhOgLwhYg7r33Xj7//HPq6+u55JJLuOeee/D7/QAsWrSIqVOnsnHjRqZPn05SUhKPP/54uIrSZZqmo6HTUUNFUcBoUDAbDZhMBkwGAyaTBA4hRPxRYi3dd3WdE1sU9hMaWoKG2XT8y9CNrqq+0A8q9YtN8Vw36Bv166qIdTHFG03T8WhqcHxEAUxGA2azgQSTEYvZIGMbQoiYIgEiTHTAp2r4VA0nfhTAbDKQaDFiMRsxGaVLSggR3SRA9BId8Po1vH4N8GEyKiRaTCRaJFgIIaKTBIgI8as6zS4fzS4fFpMhECwSjJEulhBCBEmAiAKBloUXuwssSRZ8fg2zSVoVQojIknehKKLr4HT7qW1yU2/34PNrkS6SEKIPkxZElPL4AjOiEsxGUpPNMk4hhOh18q4T5Tw+ldpGN80unyQqFEL0KgkQMUAHml0+ahvd3cpDJYQQ3SEBIob4NZ16u4cmp1daE0KIsJMAEYNaB7L9qgxiCyHCRwJEjPKrOrWNblwef6SLIoSIUxIgYpgONDq82J3eSBdFCBGHJEDEAYfbT73dgybjEkKIHiQBIk54fCr1TR40TYKEEKJnSICIIz5Vo67JjarJ4LUQ4sxJgIgzfk2nrskjM5yEEGdMAkQcUjWdOrtHWhJCiDMiASJOaZouYxJCiDMiASKOta68ltlNQojukAAR53yqRoPdI6k5hBBdJgGiD/D6NZocsphOCNE1EiD6CJdXxeH2RboYQogYIgGiD7E7fXi8ki5cCBEaCRB9TIND1kgIIUIjAaKP0XWolzUSQogQSIDog1SZ/iqECIEEiD7Kr+o0NsvOdEKIzkmA6MM8PlWmvwohOhVzAeKR57bw9e5j8sm3h8j0VyFEZ0yRLkBXHbTZOWiz89XuGuZPGUJmamKkixTzmp0+TEYDCWZjpIsihIgiMdeCMBkVAHYfauR/V5Ty8bdHZbD1DOlAY7NMfxVCtBVzAeKXSyczOD8VAJ9fo+TTAzy7ZjvHGl0RLlls03RoaJacTUKI42IuQPTvZ2XpnNHMnzIk2CVyoMrO71d+K62JM+RXdZqcMh4hhAiIuQABYFAUJo3K4/8tPJvignQgkLW05NMDPF+yg4ZmT4RLGLtcHr+k4xBCAGEOEJs2bWLGjBlMnz6dZcuWtXvebrdz5513MnfuXGbNmsXrr7/epftnpCTw/atHcu0lQ4OtiX1Hmnh6ZSlf75GZTt3V6JCV1kKIMAYIVVV59NFHee655ygpKWHt2rXs2bOnzTWvvvoqw4YNY/Xq1bz88ss8+eSTeL1dm5evKAoTR+by4+vOZnD/wNiE26vy2oY9vPbBHtxef4/Vqa/QdGhySFeTEH1d2AJEaWkpRUVFFBYWYrFYmDVrFuvXr29zjaIoOBwOdF3H4XCQnp6OydS9mbeZqQksnTWaq88fhNEQmOn0zZ5a/vD6t1RWN59xffoaj0/FKesjhOjTwrYOwmazkZ+fHzzOy8ujtLS0zTVLlizhrrvuYsqUKTgcDn77299iMJw+ZmVlWTt9bt5lxUwYnc9zq7dRVeukzu7hmdVlzJ86jCsmDcKgKN2vVC85Vf16k0FRyM5KxmDo2X+znJzUHr1ftInn+sVz3SD+69dVYQsQHfX/Kye9OW/evJlRo0bx0ksvcfDgQW655RYmTpxISkrKKe9dV+c45fPJZgN3zh1DyacH2FpejabpvPHBHrbvPcZ1l55FcmL0rg/MyrKetn69yeVwk5ps6bH75eSkUlNj77H7RZt4rl881w36Rv26KmxdTPn5+VRVVQWPbTYbubm5ba554403uPLKK1EUhaKiIgoKCti3b1+P/HyL2cg1lwxl0RXFwQHs8oMN/OGNUiqr4/c/QU9zevxomgz2C9EXhS1AjBs3joqKCiorK/F6vZSUlDBt2rQ21/Tv359PP/0UgGPHjrF//34KCgp6thxDs7l7wTgGZCcD0NDsZdnq7Wwpq5JZTiHQdWh2yViEEH1R2PpaTCYTDz/8MEuXLkVVVRYsWEBxcTHLly8HYNGiRfzwhz/kwQcfZM6cOei6zn333UdWVlaPlyU7LZE75o3l7S0H+Gy7DVXTWf1xBZXVzcyfMhSzKSaXg/Qal8dPcqIJk1H+nYToSxQ9xj5GV9c5sZ1BP+HXu4/x5qZ9+FryDg3ITmbJlcOjJulftI1BtEqyGElPSTjj+/SFft54rV881w36Rv26qs99JDynuB93zh9DVmrgze5IrZM/vrGNfUcaI1yy6Obyqvj8snhOiL6kzwUIgP7ZVn507TiGF2YAgYHY50vK2bK96jSv7NvsTtlcSIi+pE8GCICkBBM3zRjBJeMHAKDpOqs3V/DWR/skzUQnvH5NVqYL0Yf02QABYDAoXHX+IL477azgPhOf76jmhXfKcXnkjbAjdqdPZn8J0Uf06QDR6pyz+vGDuWNISzYDsPdwE396axu1Te4Ilyz6qJqOwy3BU4i+QAJEi4KcFO665vh6iWONbv705jYqqpoiXLLo43D7ZPGcEH2ABIgTpFst/GDuGEYPzgRaB6938O2+2giXLLroOjRLIj8h4p4EiJNYzEYWTx/OlLP7A4Fd1v7+/m42lx6VvvcTuCQFhxBxTwJEBwyKwtUXFDHnosEoCujA21sOUPLpAdnStIWuB1pYQoj4JQHiFCaPyWfJ9OGYW1JMfLKtitc27MGvyjRYAKfbJwFTiDgmAeI0Rg/OYumcUSQnBNJWle6t5aV1O2XfZgI7z8l0YCHilwSIEBTmpnLHvDFkpAT2RdhzuJFn126XLKeAw+2XsRkh4pQEiBDlZCRxx7yx5GUmAXDkmINlq8toaPZEuGSRpWk6Lo+0poSIRxIguqB1GmxRfiAr4rFGN8+sKuNYgyvCJYssh1tWVwsRjyRAdFFSgolbZo4MJvprdHh5ZnUZR45FX4ru3qJKK0KIuCQBohssJiPfu3I444ZmA4F++OfWbu/TW5k2SytCiLgjAaKbTEYD1087i/NGBvbZdntV/lKyg31H+mZqjsBYhMxoEiKeSIA4AwaDwvwpQ7hobD4AXp/Gi++Us6uyIcIli4xmmdEkRFyRAHGGFEVh5uQiLj13IAA+VePld3dSfrA+wiXrfZqmy+pqIeKIBIgeoCgKV55XyJXnFQKBQdtX39vFjoq6CJes9zlcMhYhRLyQANGDLj13IFefPwhoCRL/3E3Z/r4VJAKrq2VGkxDxQAJED5syfgCzJhcBgW1Ml7+/m219LF240yMrzIWIBxIgwuCicf2Zc+FgIBAk/r5+T59qSfhVHZ9fWhFCxDoJEGEyeWw+cy8aDBxvSWzvQ2MSTtmWVIiYJwEijC4Yk9+mJbH8/d19ZuDa7VVlQyEhYpwEiDCbPDaf2RcGxiRUTedv7+/uE+skdMDllVaEELFMAkQvuHBs/+DAtarpvPLeTvYcboxwqcLPJd1MQsQ0CRC95KJx/ZkxKbBOwq/qvPzuTvYfje+0HH5Nx+OTwWohYpUEiF409ZyBXD6hAACfX+PFdeVxn+BP8jMJEbskQPSyad8ZyNRzBgCB3E1/fbuco7XxmyrcI4PVQsQsCRC9rDUtR2uCP7dX5fmSHVTH6aZDMlgtROySABEBrQn+JrakCne4/TxfsoO6JneESxYe0s0kRGySABEhiqIw/+IhnD0ssOlQk8PL8yU7aIzDPa5lZbUQsUkCRAQZDAoLLxvGqKJMAOrsHn73j6/ichWyUxL4CRFzQgoQt912Gx988EGX0zhv2rSJGTNmMH36dJYtW9bhNZ999hnz5s1j1qxZfO973+vS/eOB0WDghsuLGTogDYAjNQ5eeGcHHm98vaG6vbKZkBCxJqQAcf311/Piiy9yxRVXsGzZMurrT78ZjqqqPProozz33HOUlJSwdu1a9uzZ0+aapqYmHnnkEf70pz9RUlLC7373u+7VIsaZTQZuvHIEhbkpAByqcfDyezvx+bUIl6zn6HpgQF4IETtCChBXXnklL7zwAs8++yzV1dXMnj2bf//3f2fbtm2dvqa0tJSioiIKCwuxWCzMmjWL9evXt7lmzZo1TJ8+nQEDAtM+s7Ozz6AqsS3BYuTmq0YyIMcKwL4jTby2YU9cTRGVwWohYoupOy8ym80kJCRw//33M2XKFB544IF219hsNvLz84PHeXl5lJaWtrmmoqICv9/PjTfeiMPh4KabbmL+/Pmn/flZWdbuFDvqZQH/7/pz+f9e+RfHGlyUVdTxztZKvnfVSBRFiXTxeoTPr5GTkxrpYoRVPNcvnusG8V+/rgopQLz33nu88sor1NbWsnjxYkpKSrBarfj9fq688soOA0RH/c0nv8mpqkpZWRkvvPACbrebG264gfHjxzNkyJBTlqeuLn4XlmVlWbn5qhE8s6qMZpePj785ghG4qmWnulhnTTTjc3sjXYywyclJpaYmPlfHx3PdoG/Ur6tCChArV67k9ttvZ8qUKW1fbDLx0EMPdfia/Px8qqqqgsc2m43c3Nx212RmZpKcnExycjITJ06kvLz8tAEi3mWnJXLLzJE8u2Y7bq/Kpm+OYE0yMeXsAZEu2hlzenwYdR1DnLSIhIhnIY1BPPPMM+2CQ6tp06Z1eH7cuHFUVFRQWVmJ1+ulpKSk3bWXX345X3zxBX6/H5fLRWlpKcOGDetiFeJT/2wrN101ApMx8Eb6zpaDfLmrJsKlOnO6Dm6Z8ipETAgpQCxevJjGxuPpqRsaGliyZMkpX2MymXj44YdZunQpM2fO5Oqrr6a4uJjly5ezfPlyAIYNG8aUKVOYO3cuCxcu5LrrrmP48OFnUJ34Mjg/jUVXDMfQ8mH7jY172Xnw9DPIop3sWS1EbFD0ECanz5s3j1WrVp32XG+ornNii+N+wqwsa7sxln/trOb1jfsAMBsN3DZ7FIPyYnMwrbV+2WkJmE3GSBenx8VzP3Y81w36Rv26KqQWhKZpOJ3O4LHD4UBVpZugt0wYkctVkwKD1D5V48V1O6muj+3kfvG4WlyIeBNSgJg9eza33norq1atYtWqVdx2223MnTs33GUTJ5gyvj8Xj+sPBNYT/PXt2M7bJHtWCxH9QprFdMcdd5Cbm8uGDRvQdZ0bbrghpPUKoucoisJVFwyi2eXj6z3HaHR4+es75dwxdwxJCd1azhJRrWnArYnmSBdFCNGJkN9ZrrnmGq655ppwlkWchkFRuHbqUBxuH7sPNVJd7+Kld3dy68xRmE2xl3fR4ZYAIUQ0CylA1NbW8vLLL1NZWYnff7zvuK/mTookk9HA4iuG89za7Rw+5uBAlZ1/bNjN4iuGYzDE1toCTdNxefwx2QISoi8I6S/znnvuYdiwYUyePBmjMf5mnsSaBIuRm68eyZ9XbaOuycP2inrWfFLB3IsGx1xKDqdbAoQQ0Sqkv8ympiZ+9atfhbssogtSkszcMnMUf35rGw63n8+220i3Wrj03IGRLlqX+FQNn1+NyymvQsS6kDqui4uLsdls4S6L6KLstERuvnoklpbxh/e2VvKvndURLlXXOWTKqxBRKeQWxNy5czn33HNJSEgInpcxiMgryElh8fThvLRuJ5qu8+amfaQkmRkxKDPSRQuZx6uiahpGQ+wNtAsRz0IKELNnz2b27NnhLovopuGFGVw7dSgrP9yLpsPy93ezdM5oCnJSIl20kOgExiJSky2RLooQ4gQhBQiZ3hr9vjM8hyaHl/e2VuL1B1Zb3zlvDNlpiZEuWkhcHj8pSeaYG2QXIp6F1KavqKhg0aJFwWysZWVl/P73vw9rwUTXTT1nABeMzgPA4fLxwtvlNLtiIzGeJluSChF1QgoQv/zlL7nrrrtITQ0kexo1ahTr1q0La8FE1ymKwuwLBzN6cGD8obbJzUvryvH6YuONV7YkFSK6hBQg7HY7l1xySbD5bzAYMJtlBWw0MhgUrp9WTFFLttdDNQ6Wr9+NGgN5j7x+Db+qRboYQogWIQUIo9GIz+cLBgibzYZBZpxELbPJwI0zhtMvPTD+sPNgA6s37+9wG9ho45RWhBBRI+QNg+6++27q6+v5/e9/z+LFi7n11lvDXTZxBpITzdwycySpSYGW3tbyajZ8eTjCpTo9t8cfE4FMiL4gpFlM8+fPp6CggA8++ACXy8WTTz7JxIkTw102cYYyUwML6ZatKcPr01j/r0OkWy1MHJl7+hdHSOtgtaTfECLyQv4rnDhxogSFGDSgn5Ul04fz4juBhXRvfbSP1OToXkgn+ZmEiA4h/RUuWLCgw/npK1eu7PECiZ5XXJDBgqlDWREjC+kC+Zm0mExhLkQ8CSlA3H///cHHHo+HkpIScnOjt5tCtHfu8BwaY2ghndPjJ90kK6uFiKSQAsSkSZPaHF988cUySB2Dpp4zgEaHl8+224IL6e6YN4aUpOibsuz2+ElNNmOQldVCREy32vDNzc1UVlb2dFlEmCmKwpwYWUinIwvnhIi0Lo9BaJrGoUOHuOWWW8JaMBEerQvp/lKynYO25uBCuu9dOQJjlO1I55ItSYWIqC6PQRiNRgoKCsjLywtboUR4mU0GbpoxgmdWl1HT4GbnwQZWfbSPay4ZGlXJ8vyajsenkmCWzYSEiIRujUGI2JecaOb7V4/iz6u2YXf6+GJnDWlWC1dMLIx00dpwuv0SIISIkJACxAUXXNDhJ0td11EUhU8//bTHCybCLzM1ge9fPZJlq7fj8als+PIwqckWzh8dPa1Dj082ExIiUkIKEIsWLaKhoYHrr78eXdd5/fXXycvLY+bMmeEunwiz/tlWvnflcF54pxxV01n98X5Sk82MHpwV6aIFOVx+0qwy5VWI3hbSx7KtW7fyn//5n4wcOZJRo0bx0EMPsXHjRgYOHMjAgQPDXUYRZsMGpnPdpcMA0HX4+/rdHKiyR7hUx7k8fsnyKkQEhBQgqqurqaurCx7X1dVRU1MTtkKJ3jf+rH7MmlwEgF/VeXFdObY6Z4RLFaADdmdsbHwkRDwJqYvp5ptvZt68eVx22WUAbNy4kTvuuCOsBRO976Jx/WlyePmo9Chur8oL7wQW0mWkJES6aHh8qsxoEqKXhRQglixZwoQJE9i6dSu6rrNkyRJGjBgR7rKJCJhx/iCaXT6+2n2MRoeXF94p5wdzxpCcGPnkeXanl4T0pEgXQ4g+I+S/+oKCAlRVZcyYMeEsj4gwg6Jw7dShNLt87D7USHW9i5feLefWWaOwmCL76d2v6jjd/qgIVkL0BSGNQWzcuJFZs2Zxzz33APDtt99y5513hrVgInKMBgOLpw+nIMcKwEFbM8vf342qRX6g2OGWsQghektIAeLpp59m5cqVpKWlATBu3DgOHjwY1oKJyEowG7n56pFtti19c9O+iO/2pmo6Hm/05Y4SIh6FvPooJyenzbHFIvPS45010cwtM0cF1yB8uesY73x2MOJBQvatFqJ3hBQgrFYrx44dC66m/uyzz0hNTT3t6zZt2sSMGTOYPn06y5Yt6/S60tJSRo0axbp160IstugtrautkxIC4w+bS4+y6ZsjES2Tx6fKugghekFIAeJnP/sZt99+O4cOHeLGG2/kvvvua5PAryOqqvLoo4/y3HPPUVJSwtq1a9mzZ0+H1/3mN7/h4osv7l4NRNjlZyVz04yRmI2B/y7vfl7J1vLqiJZJWhFChF9I00HGjx/PSy+9xJdffgnAueeeGxyP6ExpaSlFRUUUFgaSv82aNYv169dz1llntbnu5ZdfZsaMGXz77bfdKb/oJUX5qSyeXszL7+4K7m2dZDEydmh2RMrj8vhJTTJHVfZZIeLNaQOEqqp897vf5fXXX2fq1Kkh39hms5Gfnx88zsvLo7S0tN0177//Pi+++GKXAkRWljXka2NRtNZvcpYVk8XE86vL0HV47YM99Mu2MnpI14JET9XPmpKANQp3w8vJOX33a6yK57pB/Nevq04bIIxGI5mZmXg8HhISQl9R29FA5smf9h577DHuu+8+jMauza+vq3N06fpYkpVljer6DctPZfZFg1nzcQV+VedPr5dy26xRDMoL7Q+rJ+vX1OikX5QtnMvJSaWmJnryWPWkeK4b9I36dVVIXUyDBw9myZIlzJgxg+Tk5OD5JUuWdPqa/Px8qqqqgsc2m43c3Nw212zbto17770XgPr6ejZu3IjJZOKKK67oUiVE75o8Jh+Xx8/7XxzC59d44Z1ybp8zmv7Zvdvy8auBKa8JFkm/IUQ4hBQgHA4HxcXF7Nu3L+Qbjxs3joqKCiorK8nLy6OkpIT/+Z//aXPNhg0bgo8feOABLr30UgkOMeKycwfi9qhs/jaQt+mvb5fzgzmj6ZfRu5/om5xe+pkTZSxCiDA4ZYD49a9/zQMPPMATTzzBxx9/zEUXXRT6jU0mHn74YZYuXYqqqixYsIDi4mKWL18OBPaYELFLURSuvmAQLq+ff+2sodnl4y8lO/jB3DFkpvZecj9V02l2+UhNlnU5QvQ0RT/FqqdrrrmGN998s93jSKquc2KL437CaB+DOJmm6fxjw26+3RdIB5+dlsjtc0eT1skbdjjqpwDZ6YmYjJHfdS6e+7HjuW7QN+rXVaf8izoxdkR69ayITgaDwsLLzmLEoAwAapvcPF+yo1dzJulAk8Pbaz9PiL7ilAHC6/Wyd+9e9uzZ0+Zx65cQACajgcVXDGfogMDamOp6F38t2YGrFxezef1ar/48IfqCU3YxTZs2rfMXKgrr168PS6FORbqYopfHp/LXt3dw0NYMQGFuCrfOHNVmllE462c0KOT08iD5yeK5myKe6wZ9o35ddcpB6hNnGQlxOglmI9+/eiR/WbuDw8ccVFY38+K6cr5/9UgsvbATnKrpuDx+khJkvwghekLkR/VEXEm0mLhl5kjyswLrZSqq7Lz07k68/t5J0d3skv0ihOgp8lFL9LjkRDO3zhrFs2u2U9PgYt+RJl55dxc3zjizbWp3H2rgi/Jq6u0eMlMTmDgyl+KCjDbXSCsi9m3bX8vm0qPUNLjIyUji4rP7M7aL6VxEz5AWhAiLlCQzt80eFdxwaM/hRl55bye+brYkdh9q4N3PK6lt8qDpUNvk4d3PK9l9qKHdtdKKiF3b9tfy+sZ92OpdaDrY6l28vnEf2/bXRrpofZIECBE2ackWls4eTXZLkNh9qJE/v/EtPn/X93L4opP04h2db21FiNizufRol86L8JIAIcIqzRoIEllpgdXVZftqW1oSXQsS9XZPl85LKyI21TS4Ojnv7uWSCJAAIXpB+klBYvehRl7u4sB1Z+k7OjsvrYjY1Nk05ZyMxF4uiQAJEKKXZKQkcPvs0eRmBt4A9hxu5KV1O/H6QgsSE0fmduk8gENaETHn4rP7d+m8CC8JEKLXpKckcO/iCcGB631HmnjhnXI83tMHieKCDGZMKiQ7LQGDAtlpCcyYVNhuFtOJ/JqO2yutiFgydkg2C6YOJS8zCYOikJeZxIKpQ2UWU4ScciV1NJKV1LEtK8tKxaF6/rJ2e7BfuTA3he9fPTIsU1NNRqVXNxWK59W48Vw36Bv16yppQYhel5Zs4fY5Y4KL6Sqrm3lu7fawDCy3biokhOg6CRAiIlKSzCydPZqBOYFd6I7WOnl2zXYaw5CVVWY0CdE9EiBExCQnmlr2s04BAlMcl60uo7apZ6c0+lRNWhFCdIMECBFRgdxNozhrYDoQWNewbFUZVXXOHv05dqdX9jQRooskQIiISzAbuemqEYwenAmA3eXj2TVlHLT13IChX9NxyroIIbpEAoSICiajgUVXDOc7w/sB4PKo/GXtDnYerO+xn9Hs8qFp0ooQIlQSIETUMBoUrp06jIvG5QOBsYOX393JV7tqeuT+uh5onQghQiMBQkQVg6Iw84IiZkwqBEDTYcWHe9n09ZEeGUNwefzdzigrRF8jAUJEHUVRmHrOQBZMHYpBCZxb9/lB1nxc0SNdRE0OaUUIEQoJECJqTRiRy/euHIHZFPhvumW7jVf/ueuMd6fzqRoOtwQJIU5HAoSIaiOLMlk6ezTWxEAajh0H6vnL2h3YnWe2oK7Z6cOvdn1fCiH6EgkQIuoV5qZw5/yxwY2HKqub+fOqMmxnsFZCB5rCsGpbiHgiAULEhOy0RO6cN4aivEDCsXq7hz+vKutwy9FQef0aTulqEqJTEiBEzLAmBva5PueswFoJj0/lxXfK+bSsqtsznOwu6WoSojMSIERMMRkNLLxsGJdPKAAC02DXfFzBqs37u/VGr+vQ0OyRNBxCdEAChIg5iqJw+YQCbrj8LEzGwDzYz3dU8/zbO7qVudWv6jIeIUQHJEB2kQ/iAAAd/UlEQVSImHX2sH7cMXcMaVYLABVH7fzfm99yuKa5y/dyeVWcbsnVJMSJJECImDYwJ4UfXTOWwtxAyvCGZi/PrC7jy26k57A7vbLKWogTSIAQMS812cLtc0Zz3shcINBltPLDvV0el9AJzI7y+WXQWgiQACHihMlo4JpLhnLNlCEYW/JzfLbdxrLVZdTbPSHfR9Oh3u6WmU1CIAFCxJnzRuXxg7mjSW8ZlzhU4+APb5RS3oW04ZoOdXaPBAnR54U1QGzatIkZM2Ywffp0li1b1u751atXM2fOHObMmcMNN9xAeXl5OIsj+ojC3FTuXjCO4oLALnUuj8pL63byzpYDIb/pa5pOnd2D1ydjEqLvCluAUFWVRx99lOeee46SkhLWrl3Lnj172lxTUFDAK6+8wpo1a7jrrrv4xS9+Ea7iiD7Gmmjm5qtGcvmEAloSwvJR6VGWrS6jLsQ9rzVNp97ukdXWos8KW4AoLS2lqKiIwsJCLBYLs2bNYv369W2u+c53vkN6euBT3jnnnENVVVW4iiP6IIMhsF7illmjSE0yA4Eup9+//i1f7zkW0j10oMnpo6HZgyaL6UQfYwrXjW02G/n5+cHjvLw8SktLO71+5cqVXHLJJSHdOyvLesbli2ZSv541KcvKyKH9eLFkO2X7avH4VF7bsIf9VXYWXTmC5ERzaDcyGshIS8BsMp7yspyc1B4odXSK57pB/Nevq8IWIDpKXaAoSgdXwpYtW1i5ciV/+9vfQrp3XZ3jjMoWzbKyrFK/MFl0+Vl8kpvCu58fRNV0tm63setAPQsvG8bQAekh3aOmxk5qsoXkxI7/dHJyUqmpsfdksaNGPNcN+kb9uipsXUz5+fltuoxsNhu5ubntrisvL+ehhx7i//7v/8jMzAxXcYTAoChcfHZ/fnjNWHIzkwBodHj5y9odrP2kIqSNiAJdTl4amj2omsxyEvEtbAFi3LhxVFRUUFlZidfrpaSkhGnTprW55siRI9xzzz3893//N0OGDAlXUYRoo3+2lR9dM44Lxwa6QHXgk21V/H7ltxyoCu0TpNurcqzRLQPYIq6FrYvJZDLx8MMPs3TpUlRVZcGCBRQXF7N8+XIAFi1axB//+EcaGhp45JFHADAajbzxxhvhKpIQQWaTgdkXDmZkUSZvbNxLQ7OX2iY3y1aXMXlsPtPPKyTBfOqxBl0PDGC7PCrpKRZMRllWJOKLosdYnuPqOie2OO4nlDGI3le2v5Z3thyk7oQV1ylJJvqlJ+JXdTJTE5g4MpfigoxO76EAKclmBhdmdasfe9v+WjaXHqWmwUVORhIXn92fsUOyu1OdHrf20wo+/OowDrcfa6KJS88dyOzJgyNdrB4nYxDtha0FIUQs2H2ogQ++OkJigoksBRqbvaiaTrPLT7OrmaQEIz5V593PKwE6DRI6YHf6qKl34Ve1LrUmtu2v5fWN+4LHtnpX8DjSQWLtpxWs/bgCCEwyaXb6gsfxGCREW9ImFn3aF+XVwceJFhM5mUnBXE4QWIVdXe/E4fKxdYfttPfz+lVqG900NoeeqmNz6dEune9NH351uEvnRXyRACH6tJMT+RkUBYMBjEaCrQBdD8x22lnZSGX16bsgdAL7SxxrdNPQfPrssDUNrk7Oh7biO5w624DJ0Y2NmUTskQAh+rTM1IR254wGA2ajkZyMRNKsFlqX7/j8Gn96q4yVH+7F7gxtBzq3V6W2yU293YOnk7xOORlJnZxPDK0SYZSS1PEiQmsn50V8kQAh+rSJI9uvzUlONGFNNKEoCilJZnIzkki0HJ/R9OWuGp76xzds/PpwyHtHeHwq9XYPxxpdON2+Nmk7Lj67f4ev6ex8b7r03IFdOi/iiwxSiz6tddD5i/Jq6u2elhlLhW3OZWcmMXNyEQaDwtqPK7DVu/D4VN79vJLPttu4ctIgzh6WjaGTTAEn8qs6TU4fdpePRIuJ5ARTcCA6MIvJTU5GYtTMYmodiP7wq8M43X6sSea4ncUk2pNprlEmGqeB9qRYr5+q6Xy+3cb7/zqEy3N8D+uBOVZmnDeISWcP6HL9zEYDSQkmEi1GDIbTB5lI6QvTQOO9fl0lASLKxPob6OnES/1cHj8ffnWYT7ZVoWrH/4RGFGUy7dyBwT2yu0IhsIAv0RKdwaIvvIHGe/26SrqYhOiGpAQTV19QxPmj8/jnF5V8s6cWgJ0H6tl5oJ6RgzK5fMJABuaEHih0wOvX8Pq9NDkhwWwk0WIkwWIMqftKiJ4mAUKIM5CVlsj104q5ZPwA3ttayc6DDQCUH6yn/GA9o4oyuew7AynoQqBo5fGpeHwqigMsrcHCHH0tCxG/JEAI0QP6Z1u5+aqR1Dl9vLFhN/uONAGw40A9Ow7UU1yQztRzBjKkf2qnae87o3M8WEBgzCLBYsRiMmA2Gbp8PyFCJQFCiB50VkEGS2ePZv/RJtb/61AwUOw+1MjuQ40MykthytkDGFWU2e2WgE/V8LkC02sVBSwmIxazgQSzURIGih4lAUKIMBjSPy0YKDZ+fYRdlYGup4O2Zl795y6y0xK5aFw+3xmeg+U0WWNPRdePty7s+DAoYG4JGBaTEbNJAoboPgkQQoTRkP5pDOmfxpFjDjZ+fZht++vQdahtcrP64wre21rJxJG5XDA6j6y0M185rekndkcFAobFbAy2MqSFIbpCAoQQvWBAPyuLrhhOXZObj7dV8a/yarx+DbdXZXPpUT4uPcqIQZlMGp3L8IKMHhuI1vRAug+3NzB+0baFYTjt/tqib5MAIUQvykpLZM6Fg7liQgFflFfzaVkVDc1edI7PfMpIsTBxZC4ThueQntI+V9SZaNvCCIxhmI0GLGYjZqMBs9kgU2pFkAQIISIgKcHElPEDuGhcf8oP1rOlzMaew40ANDR7ef+LQ6z/4hBnFaQzYUQuo4oywzKeoOutay+O55QyGZXg+IXZJN1SfZkECCEiyGBQGD04i9GDs6htdPP5Dhv/2lWD0+1H5/jsp0SLkbFDszm3uB9F+alh/ZTvV3X8qh9aMqEblEDqc0uSBafbh9FgwGhUMBgUaW3EOUm1EWXiJRVFZ6R+p+dXNXYcqOdfO2vYfaiBk/9C060Wxg3L5uxh2QzsZ+21dRAd1U0BFIOCyaBgNBowGwPfTUYFoyG2Wh6SaqM9aUEIEWVMRgPjhmYzbmg2jQ4vX++u4avdx6iuD2ws1Ojwsrn0KJtLj5KVmsCYIVmMHZrFwJyUXv9ErwO6puPVdPBrnLj1kQIYDYGWhqIoGJRAMDGe8BW4qu0ddb3lvnrLYz0wVmJoeU1gUydpufQGCRBCRLF0q4Wp5wzkkvEDOFrr5Ovdx/h2Xy2NjsCGRXV2Dx+VHuWj0qOkWS2MKspkVFEmQwekRXzsQAf8mh4YGQ8Dg0IwWBhOCBxGg4LJqARbVoEgo7ds/KSgKC0BR1FkFfppxFwXk9Pto6bGTmuhOyr9iVXST3qgt72w3XWnu18HL0U/6eYaxz/9dJV0wcS23qifputU2pop3VvLtv212J3tt/+0mAycVZDOiMIMhhdm9MhsqHj83SkcDxbZ/VKor3O0ae0oHA8wHcU5Y5vg1D5gRZM+ke4biKl+Qk3XT2gqB/6TnXysnXCclWXl2LHmNudam9zxIB7fZE7U2/XTdJ1D1c2U7a+jbH8ddSftsd0qLzOJ4oIMzipIZ3B+ardWb8vvrmtax2eMrUHjhC62YCBRlGD3mUJgvxFV1fFrWkvwOum6lu+tdAg5EMkYRBQytLZnQ5SdnoTm9bc73xootGDLREfTCG5d2VnwaXMcuDCuAk5fZ1AUBuWlMigvlavOH0R1g4vylgSBldXNwVasrd6Frd7F5m+PYjQoFOalMLR/GkMHpFGYmyopOcKgdXxGQ4eOtyPvMYYTusxa325ObAkZFIWcbtxXAkSMaP3FG9oN6p251kZk2+63Ex/qJ13f+b0U5fjzrYHsxOuz0xPRff52P+PEn3Oq+7cOWHZ0/sRydnaPdnXtqEuxo5t3cI2m6fg1HS1MfexdpSgKeZnJ5GUmM/WcgTjdfvYcbmBXZQO7KxuxuwJdUaqmU3HUTsVROxu+PIzJqFCQk0JRfiqD8wPBJilB3hpiSWtPRE9/9JP/BSI4UNfmfbfzgzOSaDGRcAbJ6aKRrustawc0rIlmmk0GNE1H1U4Orb0rOdHE2cP6cfawfui6jq3exd7Djew51EhFlT24mtqv6lRU2amosrOx5bU5GUkMykthUG4KBbkp5GYmR64iImIkQAhxhhRFwWxSMJsMZKQm4HMfT7qnaTo+v4ZP1fD6VHx+LSJBQ1EU8rOSyc9K5qJx/VE1nSPHHOw70kjFUTsHbPZgviaAmgYXNQ0u/rWzBgik4xjUP5XcjCQG9rMyoJ+VnIyklqmqIl5JgBAijAwGhQSLkQSMkGRG1/VALiRvIB9SpHqnjAaFwtwUCnNTmHpOoHvCVufkQJWdg7ZmDlbbqWs6PuDtUzX2Hmpk76HG4DmTUSE3M5n+WcnkZyeTl5VMXmYSKUlmmT4aJyRACNGLFEUh0WIi0RL409N0Pdgd5fWpuLxqRMY0DIpC/2wr/bOtXDAmcK7Z5eNwTTOV1c0crnFwpNbRZkqtXw20Qo4cazvzJznRRG5mErkZSeQEvxJJT0mIuqmf4tQkQAgRQQZFwWBUMBkhwWwkNRm8vkB6bo9PRY3gAHhKkpkRgzIZMSgTgMzMZCoq6zlyzMHhYw6q6pwcrXVSf9LUWqfbHxwEP5HJqNAvPYnstESy0xPITkskKy2RrLQE0qwJ0l0VhSRACBFlLGZjcJ2Cz6/haRm78Pkj1yUFgdZPekoC6SkJjBqcFTzv9vqx1bmoqnNiq3diqwuMXzS72i7g86s6VXVOquqc7e5tUBQyUi1kpCSQmZoQ/J6eEjiXbrVEfGV4XyQBQogo1ppyu5Vf1fD6AsHC69ci2sJolWgxUZSfSlF+24VYTrc/ONh9rNFFTYObY41u6prc7cqt6Tp1TZ424x4nsyaaSLdaSLNaSE0OfE9LNpOabCGl9XuSKeaSBEYzCRBCxBCTsXV/hsCfrqpp+P06PlWLilbGiZITOw4cmqbT5PRyrNFNvd1DXVMgaDQ0e6m3e9q1PFo53H4cbj9Hatu3QE6UlGAiJclMSpIZa5IJa2LgcXKiCWuiieREM8kJJpITA19mo0EG1TshAUKIGGY0GDBaCMySaqFqgWDRujbD39LSiJK4gcGgkJES6EbqiNev0tDspbHZQ2Ozl4ZmD40OL00OL40OL43N3uAajo64PH5cnkDrJRQmoxIIKskWLCYDSRYTSQnGwGSCBCOJlpbHFmPwK8FsIsFsIMES2O87XrPLSoAQIs4EgkbbbhZdD8yU8quBwKG2rABX1egKHgAWk5HcjMAsqM54fSp2p49Gh5dmlxe704fd6aPZ1fbL4fKdthvOr+rB13eX2WQgwRzY6zuhZQzJYjK0+242Gdrs1nfil8nY8rjlu6nlscnYskFTBFo5YQ0QmzZt4rHHHkPTNBYuXMgPfvCDNs/rus5jjz3Gxo0bSUxM5Ne//jVjxowJZ5GEiFvb9teyufQoNQ0ucjKSuPjs/owdks1fSrazdUc1PlXDbDRw3qhczh+d1+baC8fmM6ook2/31fLJtiqONbrISkvkvJF5VFY382lZFU6Pn2SLkUmj87js3IJOy7H7UANflFdTb/eQmZrAxJG5AO3OFRdkhPz64oIMPvjqEJ9vtwXKkWA6ZTla71HX5A4MqhdlkpORxJ7DjZQfqMfu9GExG8hKTcRsMuDy+vH5dexOLy6PH49X7VLQDHTvaRBao6VbAmnMA5sxtQYNk9EQ3KzJaFQwtez2F9hvo3XjpsDz40bkdflnhi2bq6qqzJgxg7/+9a/k5eVx3XXX8dRTT3HWWWcFr9m4cSMvv/wyzz77LN988w2PPfYYK1asOO29Yymba1f1hV2tpH49b9v+Wl7fuK/deWuiiR0V9W3O6QTGB3JO+oQ+YUROcOV0q4ZmDw6XL5gITtMCe1fPvHAwV00aFEgCqQXyYOk6bK+oY+0nFW3eXN0ePwpgsbT9PDpjUmG7ILH7UAPvfl7Zrh75WUl8u7e23flLvzOwXZDo7B5jhmRStr++3fnWcpyYzXVXZT3rPqsMZmPW9EDivbFDs8hISQhOQz5x0WPrjDOPVw3u8926ej4arPmfeV1+TdhaEKWlpRQVFVFYWAjArFmzWL9+fZsAsX79eubPn4+iKJxzzjk0NTVRXV1Nbm5uuIolRFzaXHq0w/PlB9q/IQK43O0zBn/41WFSky1tzjU7fei6jsHYkq+rpZvjo2+OMO+iIe3u8eWuGownTUd1taTwsJ507237arlobP9gcAF4Y2MtJqPSLpliR8EBYGt5NXMvGtImn+LXu4/RbiKTDlt3VJOSZG53j69213D20GySEgLjDK33MBoVjCflIWtyeJl38dAOywId7x2j6YF0K15fIHD4TggcPlXD52tJxeI/PtHA59dQ1cDkA7/adkwp8FjD39JF6Fdbvrd0IaqqHpi80NKVeCbCFiBsNhv5+fnB47y8PEpLS095TX5+Pjab7bQBojt5zWOJ1C+2RaJ+9c3eDlN2ax1lv21J937y9Q63n6y0xDbnWtPJBxM6tnx3uv0d1rOjcqiqDkr7n9fo9JGXl9buXEd7Vfg1HYup/Xm3R2VwYVabc81uf3Cl+olqGt3062Bcw+lROWtIP4Bg/R0etcOMtk6vyvCh/dqdj1dhCxAdRdKTp5KFco0Q4vSe+snUSBcBOPNy9EQ9ouUe8SBsK0ry8/OpqqoKHnfUMjj5mqqqKuleEkKIKBG2ADFu3DgqKiqorKzE6/VSUlLCtGnT2lwzbdo03nrrLXRd5+uvvyY1NVUChBBCRImwdTGZTCYefvhhli5diqqqLFiwgOLiYpYvXw7AokWLmDp1Khs3bmT69OkkJSXx+OOPh6s4Qgghuihs01yFEELENslqJYQQokMSIIQQQnQoqnMxeTwelixZgtfrDa7M/vGPf0xDQwM//elPOXz4MAMHDuR///d/SU9Pj3Rxu6V1fCYvL49nnnkmruo2bdo0rFYrBoMBo9HIG2+8EVf1a2pq4qGHHmLXrl0oisLjjz/OkCFD4qJ++/bt46c//WnwuLKykh//+MfMnz8/Lur3wgsvsGLFChRFYfjw4TzxxBO4XK64qBvAiy++yIoVK9B1nYULF/L973+/W397Ud2CsFgsvPjii6xevZq33nqLjz76iK+//pply5YxefJk3nvvPSZPnsyyZcsiXdRue+mllxg2bFjwOJ7qBoH/qKtWreKNN94A4qt+jz32GFOmTGHdunWsWrWKYcOGxU39hg4dyqpVq4K/u6SkJKZPnx4X9bPZbLz00ku8/vrrrF27FlVVKSkpiYu6AezatYsVK1awYsUKVq1axYcffkhFRUW36hfVAUJRFKxWKwB+vx+/34+iKMEUHQDz58/n/fffj2Qxu62qqooPP/yQ6667LnguXurWmXipX3NzM1u3bg3+7iwWC2lpaXFTvxN9+umnFBYWMnDgwLipn6qquN1u/H4/breb3NzcuKnb3r17GT9+PElJSZhMJs477zz++c9/dqt+UR0gIPCLnDdvHhdeeCEXXngh48ePp7a2NrheIjc3l7q6ugiXsnsef/xx/u3f/g3DCYlj4qVurW677TauvfZa/vGPfwDxU7/KykqysrJ48MEHmT9/Pj//+c9xOp1xU78TlZSUMHv2bCA+fn95eXnceuutXHbZZVx88cWkpKRw8cUXx0XdAIYPH84XX3xBfX09LpeLTZs2UVVV1a36RX2AMBqNrFq1io0bN1JaWsquXbsiXaQe8cEHH5CVlcXYsWMjXZSwWb58OW+++SbPPvssr776Klu3bo10kXqM3+9n+/btLFq0iLfeeoukpKSY7ZI4Fa/Xy4YNG7jqqqsiXZQe09jYyPr161m/fj0fffQRLpeLVatWRbpYPWbYsGEsXbqUW2+9laVLlzJixAiMxvZ5rEIR9QGiVVpaGueffz4fffQR2dnZVFdXA1BdXU1WVtZpXh19vvzySzZs2MC0adO499572bJlC/fdd19c1K1VXl4g/3x2djbTp0+ntLQ0buqXn59Pfn4+48ePB+Cqq65i+/btcVO/Vps2bWLMmDH06xdIUBcP9fvkk08oKCggKysLs9nMlVdeyVdffRUXdWu1cOFC3nzzTV599VUyMjIoKirqVv2iOkDU1dXR1NQEgNvt5pNPPmHo0KHBFB0Ab731Fpdffnkki9ktP/vZz9i0aRMbNmzgqaee4oILLuA3v/lNXNQNwOl00tzcHHz88ccfU1xcHDf1y8nJIT8/n337AnswfPrppwwbNixu6teqpKSEWbNmBY/joX4DBgzgm2++weVyoet6XP7uamsD6dGPHDnCe++9x+zZs7tVv6heSV1eXs4DDzyAqqrous5VV13F3XffTX19PT/5yU84evQo/fv353e/+x0ZGR3vThULPvvsM55//nmeeeaZuKlbZWUlP/rRj4DAONLs2bO566674qZ+ADt27ODnP/85Pp+PwsJCnnjiCTRNi5v6uVwuLr30Ut5//31SUwOpvePl9/f000/z9ttvYzKZGDVqFI899hgOhyMu6gawePFiGhoaMJlMPPjgg0yePLlbv7uoDhBCCCEiJ6q7mIQQQkSOBAghhBAdkgAhhBCiQxIghBBCdEgChBBCiA5FdTZXIU5l4cKFeL1efD4fFRUVFBcXAzB69GieeOKJCJcuNGVlZVRWVsbVSmURP2Saq4h5hw4dYsGCBXz22WeRLko7fr8fk6nzz2ErVqzgk08+4be//W2P31uIMyX/u0RcWrlyJX//+99RVZW0tDQeeeQRBg8ezIoVK1i3bh1Wq5Vdu3bRv39//uM//oMnn3ySyspKxo8fz5NPPomiKNx3330kJSVx8OBBqqqqOP/88/nFL36B2WzGbrfz+OOPs3v3bjweDxdeeCH3338/BoOBRYsWMWnSJL766iuSk5N5+umng4sEPR4P48eP55FHHqGpqYk//vGPOBwO5s2bx/nnn8+SJUtYvHgxH3/8MQAHDhwIHh84cIBFixZx/fXXs2XLFq699lrmzZvHU089xRdffIHX62XUqFH88pe/JCkpKcK/AREXdCFiXGVlpT5p0qTg8ZYtW/Q77rhD93g8uq7r+vr16/UlS5bouq7rr732mj5p0iS9qqpK13Vdv/XWW/X58+frdrtd93q9+syZM/UtW7bouq7rP/vZz/R58+bpDodD93q9+k033aT/7W9/03Vd1++//359zZo1uq7ruqqq+o9//GN95cqVuq7r+g033KD/8Ic/1P1+f/D5hoaG4ON7771Xf+2114Ll+clPfhIse0VFhX7hhRd2eFxRUaEPHz5cX7duXfD5p59+Wn/mmWeCx0888YT+u9/97sz+QYVoIS0IEXc2bNjA9u3bWbhwIQC6ruNwOILPT5gwIZhIcPTo0bjdblJSUgAYMWIEBw8e5Pzzzwdg5syZJCcnA4Ec+h9++CGLFi3igw8+oKysjGeffRYI5AobNGhQ8GfMmTMnmEFT0zSWLVvG5s2b0TSNhoaGbu9UlpyczIwZM9rU1eVyUVJSAgSyr44ZM6Zb9xbiZBIgRNzRdZ3vfve73H333R0+n5CQEHxsMBjaHfv9/k7vqygKEHjTf+aZZxgwYECH17YGFYBVq1ZRWlrK3/72N6xWK3/4wx84evRoh68zGo1omhY89ng8nd63tUy/+tWvOO+88zq8nxBnQqa5irjTmrXSZrMBgWSB27Zt69a93nnnHVwuFz6fjzVr1gRbFtOmTWPZsmWoqgoEMg9XVlZ2eA+73U5mZiZWq5XGxsbgp30Aq9WK3W4PHufm5uJ2u4P3Wrt27Wnr+vzzzwcDSXNzM3v37u1WXYU4mQQIEXcuuOAC7r77bu644w7mzp3LnDlz+PDDD7t1rwkTJnDXXXcxe/ZsCgsLg1uM/uIXv0DTNObNm8ecOXO4/fbbqamp6fAe11xzDQ0NDcyePZt77723zaf9iy66CLvdzty5c3n88cexWCw88MAD3Hzzzdx4442YzeZTlu/OO+9k2LBhXHfddcyZM4clS5awf//+btVViJPJNFchOnHfffcxYcIEFi1aFOmiCBER0oIQQgjRIWlBCCGE6JC0IIQQQnRIAoQQQogOSYAQQgjRIQkQQgghOiQBQgghRIf+fzrLaxgk8IZIAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"sns.set(color_codes=True)\n",
"plt.xlim(30,90)\n",
"plt.ylim(0,1)\n",
"sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
]
}
],
"metadata": {
"celltoolbar": "Hide code",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}