"
]
},
"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": 4,
"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:
Mon, 06 Apr 2020
Deviance:
3.0144
\n",
"
\n",
"
\n",
"
Time:
07:33:03
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: Mon, 06 Apr 2020 Deviance: 3.0144\n",
"Time: 07:33:03 Pearson chi2: 5.00\n",
"No. Iterations: 6 Covariance Type: nonrobust\n",
"===============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"-------------------------------------------------------------------------------\n",
"Intercept 5.0850 7.477 0.680 0.496 -9.570 19.740\n",
"Temperature -0.1156 0.115 -1.004 0.316 -0.341 0.110\n",
"===============================================================================\n",
"\"\"\""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import statsmodels.api as sm\n",
"\n",
"data[\"Success\"]=data.Count-data.Malfunction\n",
"data[\"Intercept\"]=1\n",
"\n",
"logmodel=sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \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": 5,
"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:
Mon, 06 Apr 2020
Deviance:
18.086
\n",
"
\n",
"
\n",
"
Time:
07:33: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: Mon, 06 Apr 2020 Deviance: 18.086\n",
"Time: 07:33: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": 5,
"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": 6,
"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": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/conda/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9+P/XmS2ZTBKykAVIWA07ogVRVERRFtkVqQJVq2LVVvttrfeqt9ar9qr1/nrtrW1vK1rrWlrBhSWKVlBwQ3GNBMIeCEsmIXsms53l98ckAyEBMiGTWfJ+Ph4hc07OnHw+JJP3fLb3RzEMw0AIIYQ4gSnSBRBCCBGdJEAIIYRolwQIIYQQ7ZIAIYQQol0SIIQQQrRLAoQQQoh2hS1A3H///UycOJHZs2e3+3XDMPiv//ovpk6dypw5cyguLg5XUYQQQnRC2ALE1VdfzbPPPnvSr2/atInS0lLeffddfv3rX/PQQw+FqyhCCCE6IWwB4rzzzqNXr14n/fr69euZP38+iqJwzjnnUF9fT0VFRbiKI4QQIkSWSH1jp9NJbm5u8Dg3Nxen00l2dvYpn+fxqQAoihL4DAQeKihK4BhFCZ5vuU4IIURoIhYg2svw0ZE/5vWNPpyVDSF9r5bAoQQDhxIMHiblxOPAY5Mp8Nhkovlc9wSarKwUKkOsXyyR+sWueK4b9Iz6hSpiASI3N5fy8vLgcXl5+WlbD51lGGC0PAicCfkeigJmRQkEDpOCOfhhwmRSsJi7L4gIIUR3iFiAmDJlCi+//DKzZs3i22+/JSUlJWwBoisYBqiGAfrJg4vJpGAxKVjMpuYPBYvFhEkChxAiBoUtQNx99918/vnn1NTUcMkll3DXXXehqoHxg0WLFjF58mQ2btzI1KlTsdvtPPbYY+EqSrfRdQOfbuBT9VbnLSYFq8WE1WLGZg0EDyGEiHZKrKX7rqhuCnkMItqYTAo2i4kEq5kEqxmT6VgLoyf0g0r9YlM81w16Rv1CFbEupp5M1w08Pg2PTwPAajaRYAsECyGEiBYSIKKAX9Pxu3Ua3X5MNguNLh9JCWasFgkYQojIkQARZTTdwO1VcXtVLGaFpAQLiQkWGegWQnQ7CRBRTNUM6pv8NLj9JCVYcCRaW41XCCFEOEmAiAGGAS6PSpNHxZ5gIdkugUIIEX4SIGKIATR5Vdw+FUeiFUeiRRbnCSHCRibkxyDDgEa3n8o6TzA3lRBCdDUJEDFM1w1qG33UNHjRdP30TxBCiBBIgIgDXr/G0ToPLo8/0kURQsQRCRBxwjCgocnP0To3flWLdHGEEHFAAkScUTWDqnov9S4femxlURFCRBkJEHGqyatSXeeR1oQQotMkQMQxVTeorvfS6JaxCSFE6CRAxDmDwJTY6noP+in2shBCiBNJgOghfKpOVb0HvyrTYYUQHSMBogfRdIPqeg9uryyuE0KcngSIHsYA6lw+CRJCiNOSANFDSZAQQpyOBIgeTIKEEOJUJED0cPUun6yVEEK0SwJED2cAdY0+DFl1LYQ4gQQIgaoHdq4TQojjSYAQALi9Kl6fdDUJIY6RACGC6lxeWW0thAiSACGCdAMamnyRLoYQIkpIgBCtuH0aPr90NQkhJECIdtQ3yawmIYQECNEOVTNokgV0QvR4EiBEuxrdfjRdMr8K0ZPFXIB4+NnNfLPrqHSBhFnLHtdCiJ7LEukChOqAs4EDzga+3lXJ/EmDSE9JjHSR4pbHp+FXNawWc6SLIoSIgJhrQVjMCgC7DtbxvyuK+Pi7I+jSmgibepe0IoToqWIuQDy0dCIDc1MA8Ks6hZ/u55k12zha545wyeKTX9Ml46sQPVTMBYg+vR0snTOS+ZMGkWANdH3sL2/gDyu/k9ZEmDS4/TLmI0QPFHMBAsCkKEwYkcP/W3g2BXm9gMA73cJP9/Nc4XZqG70RLmF80XUDl0daEUL0NGENEJs2bWL69OlMnTqVZcuWtfl6Q0MDt99+O3PnzmXWrFm89tprId0/LTmBH145nKsvGRxsTew9XM9TK4v4ZrfMdOpKLo9f8jQJ0cOELUBomsYjjzzCs88+S2FhIWvXrmX37t2trnnllVcYMmQIq1ev5qWXXuKJJ57A5wstF5CiKIwfns1PrzmbgX0CYxMen8arG3bz6vu78fjknW9XMIxAkBBC9BxhCxBFRUUMGDCA/Px8bDYbs2bNYv369a2uURQFl8uFYRi4XC569eqFxdK5mbfpKQksnTWSK8/vj9kUmOn07e4q/vjad5RVNJ5xfQQ0eVVpRQjRg4RtHYTT6SQ3Nzd4nJOTQ1FRUatrlixZwh133MGkSZNwuVz87ne/w2Q6fczKyHCc9GvzLitg3Mhcnl29lfKqJqobvDy9upj5k4dwxYT+mBSl85XqJqeqX6TZk2ykOmxndI+srJQuKk10iuf6xXPdIP7rF6qwBYj2+v+VE/44f/TRR4wYMYIXX3yRAwcOcNNNNzF+/HiSk5NPee/qatcpv55kNXH73FEUfrqfLSUV6LrB6+/vZtueo1xz6VkkJUbv+sCMDMdp6xdJtTUueqfZOx1os7JSqKxs6OJSRY94rl881w16Rv1CFbYuptzcXMrLy4PHTqeT7OzsVte8/vrrTJs2DUVRGDBgAHl5eezdu7dLvr/NauaqSwaz6IqC4AB2yYFa/vh6EWUV8ftLEG66AU0yo0mIHiFsAWLMmDGUlpZSVlaGz+ejsLCQKVOmtLqmT58+fPrppwAcPXqUffv2kZeX17XlGJzJnQvG0DczCYDaRh/LVm9jc3G5zHLqpCaPrIsQoicIW1+LxWLhwQcfZOnSpWiaxoIFCygoKGD58uUALFq0iB//+Mfcf//9zJkzB8MwuOeee8jIyOjysmSmJnLbvNG8tXk/n21zoukGqz8upayikfmTBmO1xORykIjRjcAe1kmJ1kgXRQgRRooRY28FK6qbcJ5BP+E3u47yxqa9+LVAKuu+mUksmTY0apL+RfsYRAuTSSGrV2KbcaXT6Qn9vPFav3iuG/SM+oWqx711PqegN7fPH0VGSgIAh6ua+NPrW9l7uC7CJYstum7g8cnWpELEsx4XIAD6ZDr4ydVjGJqfBgTm9z9XWMLmbeWneaY4nsstC+eEiGc9MkAA2BMs3DB9GJeM7QuAbhis/qiUNz/cKzupdZCqG7JSXYg41mMDBAT60Wec35/vTzkruM/E59sreP7tEklx3UEut/w/CRGvenSAaHHOWb350dxRpCYFZuXsOVTPn9/cSlW9J8Ili35+Tcfrl7EIIeKRBIhmeVnJ3HHVsfUSR+s8/PmNrZSW10e4ZNFPFs4JEZ8kQBynl8PGj+aOYuTAdKBl8Ho73+2tinDJopvXr6FqMm4jRLyRAHECm9XM4qlDmXR2HwBUzeAf7+3io6Ijsnr4FKQVIUT8kQDRDpOicOUFA5hz0UAUBQzgrc37Kfx0v2xpehJunyr/N0LEGQkQpzBxVC5Lpg7Fag78N32ytZxXN+yW7pR2GAZ4ZOaXEHFFAsRpjByYwdI5I0hKCKStKtpTxYvrduCVVcRtSDeTEPFFAkQH5GencNu8UaQlBzbK2X2ojmfWbqNRVhK3ouqGBE4h4ogEiA7KSrNz27zR5KTbATh81MWy1cXUNnojXLLo0iTdTELEDQkQIWiZBjsgN5AV8Widh6dXFXO01h3hkkUPmfIqRPyQABEie4KFm2YODyb6q3P5eHp1MYePRn+K7u4iWV6FiA8SIDrBZjHzg2lDGTM4EwCXR+XZtdtkK9NmMptJiPggAaKTLGYT1045i/OGB/bZ9vg0/lq4nb2HJTWHqhv4VWlFCBHrJECcAZNJYf6kQVw0OhcAn1/nhbdL2FlWG+GSRZ7bKwFCiFgnAeIMKYrCzIkDuPTcfkAgu+lL7+yg5EBNhEsWWR6fKqlJhIhxEiC6gKIoTDsvn2nn5QOg6QavvLuT7aXVES5Z5OgGkgZciBgnAaILXXpuP648vz/QHCT+tYvifT03SEg3kxCxTQJEF5s0ti+zJg4AAtuYLn9vF1t7aLpwn19D16WbSYhYJQEiDC4a04c5Fw4EAkHiH+t398iWhAGyZ7UQMUwCRJhMHJ3L3IsGAsdaEtt64JiEpN4QInZJgAijC0bltmpJLH9vV48buFY1WRMhRKySABFmE0fnMvvCwJiEphv8/b1dPW6dRJMMVgsRkyRAdIMLR/cJDlxrusHL7+5g96G6CJeq+3hktzkhYpIEiG5y0Zg+TJ8QWCehagYvvbODfUd6RlqOwG5z0ooQItZIgOhGk8/px+Xj8gDwqzovrCvpMQn+mryyuZIQsUYCRDeb8r1+TD6nLxDI3fS3t0o4UhX/qcJlsFqI2CMBopu1pOVoSfDn8Wk8V7idih6w6ZAMVgsRWyRAREBLgr/xzanCXR6V5wq3U13viXDJwksS+AkRWyRARIiiKMy/eBBnDwlsOlTv8vFc4Xbq4niPa8MAr+w2J0TMkAARQSaTwsLLhjBiQDoA1Q1efv/Pr2nyxO/qY7esrBYiZnQoQNxyyy28//77IXcPbNq0ienTpzN16lSWLVvW7jWfffYZ8+bNY9asWfzgBz8I6f7xwGwycd3lBQzumwrA4UoXz7+9PW7faXt8mnQzCREjOhQgrr32Wl544QWuuOIKli1bRk3N6TfD0TSNRx55hGeffZbCwkLWrl3L7t27W11TX1/Pww8/zJ///GcKCwv5/e9/37laxDirxcT104aRn50MwMFKFy+9uwO/qke4ZF1PNwx8/virlxDxqEMBYtq0aTz//PM888wzVFRUMHv2bP793/+drVu3nvQ5RUVFDBgwgPz8fGw2G7NmzWL9+vWtrlmzZg1Tp06lb9/AtM/MzMwzqEpsS7CZuXHGcPpmOQDYe7ieVzfsjst02ZLhVYjYYOnMk6xWKwkJCdx7771MmjSJ++67r801TqeT3Nzc4HFOTg5FRUWtriktLUVVVa6//npcLhc33HAD8+fPP+33z8hwdKbYUS8D+H/Xnsv/9/KXHK11U1xazdtbyvjBjOEoihLp4nUZR4qd3plJcVWn42VlpUS6CGETz3WD+K9fqDoUIN59911efvllqqqqWLx4MYWFhTgcDlRVZdq0ae0GiPb6mU/8g6BpGsXFxTz//PN4PB6uu+46xo4dy6BBg05Znurq+F1YlpHh4MYZw3h6VTGNbj8ff3sYMzCjeae6WJeR4eBoVSOa10+CzRzp4nS5rKwUKivjc3V8PNcNekb9QtWhALFy5UpuvfVWJk2a1PrJFgsPPPBAu8/Jzc2lvLw8eOx0OsnOzm5zTXp6OklJSSQlJTF+/HhKSkpOGyDiXWZqIjfNHM4za7bh8Wls+vYwDruFSWf3jXTRuozHr8VlgBAinnRoDOLpp59uExxaTJkypd3zY8aMobS0lLKyMnw+H4WFhW2uvfzyy/niiy9QVRW3201RURFDhgwJsQrxqU+mgxtmDMNiDrS63t58gK92Vka4VF3HK4vmhIh6HQoQixcvpq7uWHrq2tpalixZcsrnWCwWHnzwQZYuXcrMmTO58sorKSgoYPny5SxfvhyAIUOGMGnSJObOncvChQu55pprGDp06BlUJ74MzE1l0RVDMTX3zL2+cQ87Dpx+Blks0A1kNpMQUU4xOvA2bt68eaxateq057pDRXUTzjjuJ8zIcLQZY/lyRwWvbdwLgNVs4pbZI+ifE5uDacfXL9FmJi05IcIl6lrx3I8dz3WDnlG/UHWoBaHrOk1NTcFjl8uFpsXnQq5oNG5YNjMmBAap/ZrOC+t2UFET+8n9vD5NNhISIop1KEDMnj2bm2++mVWrVrFq1SpuueUW5s6dG+6yieNMGtuHi8f0AQLpKv72VuznbTKQjYSEiGYdmsV02223kZ2dzYYNGzAMg+uuu65D6xVE11EUhRkX9KfR7eeb3Uepc/n429sl3DZ3FPaETi1niQoen0pSYuyWX4h41uFX5lVXXcVVV10VzrKI0zApCldPHozL42fXwToqaty8+M4Obp45AqslNvMu+lQdVdOxmGOz/ELEsw4FiKqqKl566SXKyspQ1WNpEnpq7qRIsphNLL5iKM+u3cahoy72lzfwzw27WHzFUEym2FyZ7PFpJNslQAgRbToUIO666y6GDBnCxIkTMZtlcVOkJdjM3HjlcP6yaivV9V62ldaw5pNS5l40MCbTV7i9Ksl2a6SLIYQ4QYcCRH19Pb/+9a/DXRYRgmS7lZtmjuAvb27F5VH5bJuTXg4bl57bL9JFC5mmG3j9GglWefMhRDTpULu+oKAAp9MZ7rKIEGWmJnLjlcOxNY8/vLuljC93VES4VJ0jGwkJEX063IKYO3cu5557LgkJxxY2yRhE5OVlJbN46lBeXLcD3TB4Y9Neku1WhvVPj3TRQuL1aTJYLUSU6VCAmD17NrNnzw53WUQnDc1P4+rJg1n5wR50A5a/t4ulc0aSl5Uc6aJ1mAE0eVRSHbZIF0UI0axDAUKmt0a/7w3Not7l490tZfjUwGrr2+eNIjM1MdJF67CWwepYnY0lRLzpUHu+tLSURYsWBbOxFhcX84c//CGsBROhm3xOXy4YmQOAy+3n+bdKaHT7I1yqjjMAlyd2yitEvOtQgHjooYe44447SEkJJHsaMWIE69atC2vBROgURWH2hQMZOTAw/lBV7+HFdSX4/LGTzqLJq0p+JiGiRIcCRENDA5dccklwjr3JZMJqlXnr0chkUrh2SgEDmrO9Hqx0sXz9LrQY2dvaMAJjEUKIyOtQgDCbzfj9/mCAcDqdmEwy2yRaWS0mrp8+lN69AuMPOw7UsvqjfTGzQU+Txx8zZRUinnV4w6A777yTmpoa/vCHP7B48WJuvvnmcJdNnIGkRCs3zRxOSvMK5S0lFWz46lCES9UxuhHoahJCRFaHZjHNnz+fvLw83n//fdxuN0888QTjx48Pd9nEGUpPCSykW7amGJ9fZ/2XB+nlsDF+ePbpnxxhLo9KUoIlJlOHCBEvOpzNdfz48RIUYlDf3g6WTB3KC28HFtK9+eFeUpKifyGdrhs0eVUciTLWJUSkdChALFiwoN13citXruzyAomuV5CXxoLJg1kRYwvppBUhRGR1KEDce++9wcder5fCwkKys6O/m0Icc+7QLOpibCGdrhu4vSpJ0ooQIiI6FCAmTJjQ6vjiiy+WQeoYNPmcvtS5fHy2zRlcSHfbvFFRnWq70aNil1aEEBHRqbmqjY2NlJWVdXVZRJgpisKcGFtI19KKEEJ0v5DHIHRd5+DBg9x0001hLZgIj5aFdH8t3MYBZ2NwId0Ppg3DHKU5kKQVIURkhDwGYTabycvLIycnJ2yFEuFltZi4Yfownl5dTGWthx0Haln14V6uumRwVP4RlrEIISKjU2MQIvYlJVr54ZUj+MuqrTQ0+fliRyWpDhtXjM+PdNHaJa0IIbpfhwLEBRdc0O4L0zAMFEXh008/7fKCifBLT0ngh1cOZ9nqbXj9Ghu+OkRKko3zR0Zf61DWRQjR/ToUIBYtWkRtbS3XXnsthmHw2muvkZOTw8yZM8NdPhFmfTId/GDaUJ5/uwRNN1j98T5SkqyMHJgR6aK1IesihOheHZrFtGXLFv7zP/+T4cOHM2LECB544AE2btxIv3796NevX7jLKMJsSL9eXHPpECCQTfUf63exv7whwqVqq6UVIYToHh0KEBUVFVRXVwePq6urqaysDFuhRPcbe1ZvZk0cAICqGbywrgRndVOES9WWy6NKplchukmHuphuvPFG5s2bx2WXXQbAxo0bue2228JaMNH9LhrTh3qXjw+LjuDxaTz/dmAhXVpyQqSLFqTrBl6/RqKtw2nEhBCd1KFX2ZIlSxg3bhxbtmzBMAyWLFnCsGHDwl02EQHTz+9Po9vP17uOUufy8fzbJfxoziiSEqPnD3KTR5UAIUQ36PCrLC8vD03TGDVqVDjLIyLMpChcPXkwjW4/uw7WUVHj5sV3Srh51ghsFnOkiweAT9VRNR2LWTatEiKcOvQK27hxI7NmzeKuu+4C4LvvvuP2228Pa8FE5JhNJhZPHUpelgOAA85Glr+3C03XI1yyY2RbUiHCr0MB4qmnnmLlypWkpqYCMGbMGA4cOBDWgonISrCaufHK4a22LX1j096oGSB2+1T0KCmLEPGqw230rKysVsc2m63LCyOiiyPRyk0zR5DqCPysv9p5lLc/OxAVQcIwwCNTXoUIqw4FCIfDwdGjR4MLlD777DNSUlJO+7xNmzYxffp0pk6dyrJly056XVFRESNGjGDdunUdLLboLi2rre0JgfGHj4qOsOnbwxEuVYB0MwkRXh0KEL/4xS+49dZbOXjwINdffz333HNPqwR+7dE0jUceeYRnn32WwsJC1q5dy+7du9u97re//S0XX3xx52ogwi43I4kbpg/H2jwo/M7nZWwpqYhwqUDVjahOVS5ErOvQLKaxY8fy4osv8tVXXwFw7rnnBscjTqaoqIgBAwaQnx9I/jZr1izWr1/PWWed1eq6l156ienTp/Pdd991pvyimwzITWHx1AJeemdncG9ru83M6MGZES2X26dhs0bH7Coh4s1pA4SmaXz/+9/ntddeY/LkyR2+sdPpJDc3N3ick5NDUVFRm2vee+89XnjhhZACREaGo8PXxqJord/EDAcWm4XnVhdjGPDq+7vpnelg5KDQgkRX1k9RoHemI6ryM2Vlnb77NVbFc90g/usXqtMGCLPZTHp6Ol6vl4SEjq+obW8g88QX8aOPPso999yD2RzaO8DqaldI18eSjAxHVNdvSG4Ksy8ayJqPS1E1gz+/VsQts0bQP6djL6xw1E/1+LEnRMfCuaysFCoroy+PVVeI57pBz6hfqDr0qho4cCBLlixh+vTpJCUlBc8vWbLkpM/Jzc2lvLw8eOx0OsnOzm51zdatW7n77rsBqKmpYePGjVgsFq644oqQKiG618RRubi9Ku99cRC/qvP82yXcOmckfTIj0/Lx+LSoCRBCxJMOvapcLhcFBQXs3bu3wzceM2YMpaWllJWVkZOTQ2FhIf/zP//T6poNGzYEH993331ceumlEhxixGXn9sPj1fjou0Depr+9VcKP5oykd5q928vi9Wtouo7ZJCurhehKpwwQv/nNb7jvvvt4/PHH+fjjj7nooos6fmOLhQcffJClS5eiaRoLFiygoKCA5cuXA4E9JkTsUhSFKy/oj9un8uWOShrdfv5auJ0fzR1Fekr3J/dzezWS7RIghOhKinGKVU9XXXUVb7zxRpvHkVRR3YQzjvsJo30M4kS6bvDPDbv4bm8gHXxmaiK3zh1JalL7CynDVT+LSYlI6+VE8dyPHc91g55Rv1Cd8i3X8bEjGlbPiuhjMiksvOwshvVPA6Cq3sNzhdtxefzdWg5VN/CrsiZCiK50ygDh8/nYs2cPu3fvbvW45UMIAIvZxOIrhjK4b2BtTEWNm78VbsfdzakwZGW1EF3rlF1MU6ZMOfkTFYX169eHpVCnIl1M0cvr1/jbW9s54GwEID87mZtnjiDBdmwac7jrl5maiNUSubGIeO6miOe6Qc+oX6hOOUh9/CwjIU4nwWrmh1cO569rt3PoqIuyikZeWFfCD68c3m2rnRuafGSkJnbL9xIi3sm0D9GlEm0Wbpo5nNyMwHqZ0vIGXnxnB75uGh/wqTpen4xFCNEVZHWR6HJJiVZunjWCZ9Zso7LWzd7D9bz8zk6un35m29TuOljLFyUV1DR4SU9JYPzwbAry0tpc1+D2kWCL/Iwm0Tlb91XxUdERKmvdZKXZufjsPowOMZ2L6BrSghBhkWy3csvsEcENh3YfquPld3d0eqbRroO1vPN5GVX1XnQDquq9vPN5GbsO1ra5VtWMbh8gF11j674qXtu4F2eNG90AZ42b1zbuZeu+qkgXrUeSACHCJjXJxtLZI8lsDhK7Dtbxl9e/w6+GvnXpFydJL36y8w1uv+w4F4M+KjoS0nkRXhIgRFilOgJBIiM1sLq6eG9Vc0sitCBR0+AN6byuGzS6u3cthjhzlbXuk5z3dHNJBEiAEN2g1wlBYtfBOl4KceD6ZOk7TpXWw+1RUbXQWysicrJOsho+K01mpkWCBAjRLdKSE7h19kiy0wN/AHYfquPFdTs6vCPc+OHZIZ0HMIB6ly/ksorIufjsPiGdF+ElAUJ0m17JCdy9eFxw4Hrv4Xqef7ukQ9NSC/LSmD4hn8zUBEwKZKYmMH1CfruzmI7nU3UZsI4howdlsmDyYHLS7ZgUhZx0OwsmD5ZZTBFyypXU0UhWUse2jAwHpQdr+OvabcF+5fzsZH545fCw7elgUqB3WuAPTrjF82rceK4b9Iz6hUpaEKLbpSbZuHXOqOBiurKKRp5duy1sg8q6IXmahOgMCRAiIpLtVpbOHkm/rMAudEeqmnhmzTbqwjRm0OSRaa9ChEoChIiYpERL837WyUBgiuOy1cVU1Xf9lEbdQMYihAiRBAgRUYHcTSM4q18vILCuYdmqYsqrm7r8e7k8quxrIkQIJECIiEuwmrlhxjBGDkwHAqugn1lTzAFn1w4Y6rqB2yuJ/IToKAkQIipYzCYWXTGU7w3tDQT2mP7r2u3sOFDTpd+nu3e6EyKWSYAQUcNsUrh68hAuGpMLgF/TeemdHXy9s7LLvoemSyI/ITpKAoSIKiZFYeYFA5g+IR8IDC6v+GAPm7453GXjB41uv4xFCNEBEiBE1FEUhcnn9GPB5MGYmte2rfv8AGs+LkXXz/wPuyaJ/IToEAkQImqNG5bND6YNC+4xvXmbk1f+tbNLdqdrkkR+QpyWBAgR1YYPSGfp7JE4EgNpOLbvr+Gva7fT0HRmC+okkZ8QpycBQkS9/Oxkbp8/OrjxUFlFI39ZVYzzDNdKSCI/IU5NAoSICZmpidw+bxQDcgIJx2oavPxlVXG7W46GoqHJh6ZLV5MQ7ZEAIWKGIzGwz/U5ZwXWSnj9Gi+8XcKnxeWdnpWkG1BT7+2SwW8h4o0ECBFTLGYTCy8bwuXj8oDAH/g1H5ey6qN9nR50VnWDmgavJPMT4gQSIETMURSFy8flcd3lZ2ExB+bBfr69gufe2t7p6at+TafI5GQKAAAeZElEQVS2wSvrI4Q4jgQIEbPOHtKb2+aOItVhA6D0SAP/98Z3HKps7NT9fKpObaNPgoQQzSRAiJjWLyuZn1w1mvzsQMrw2kYfT68u5qtOpufw+jWZ/ipEMwkQIualJNm4dc5IzhueDYCqGaz8YE+nxyXcPi1sGxcJEUskQIi4YDGbuOqSwVw1aRDm5vwcn21zsmx1MTUN3pDv5/aq1J/hYjwhYp0ECBFXzhuRw4/mjqRX87jEwUoXf3y9iJJOpA1v8qhU1Xnwq7JOQvRMYQ0QmzZtYvr06UydOpVly5a1+frq1auZM2cOc+bM4brrrqOkpCScxRE9RH52CncuGENBXmCXOrdX48V1O3h78/6Qu5z8mk51vYeGJhm8Fj1P2AKEpmk88sgjPPvssxQWFrJ27Vp2797d6pq8vDxefvll1qxZwx133MGvfvWrcBVH9DCORCs3zhjO5ePyaE4Iy4dFR1i2upjqEPe8NghsVyprJURPE7YAUVRUxIABA8jPz8dmszFr1izWr1/f6prvfe979OoVeJd3zjnnUF5eHq7iiB7IZAqsl7hp1ghS7FYg0OX0h9e+45vdR0O+n08NtCZk1bXoKSzhurHT6SQ3Nzd4nJOTQ1FR0UmvX7lyJZdcckmH7p2R4Tjj8kUzqV/XmpDhYPjg3rxQuI3ivVV4/RqvbtjNvvIGFk0bRlKiNaT7GSaF9DQ7FnP776+yslK6othRKZ7rBvFfv1CFLUC011+rKEo7V8LmzZtZuXIlf//73zt07+pq1xmVLZplZDikfmGy6PKz+CQ7mXc+P4CmG2zZ5mTn/hoWXjaEwX17hXSvqqpGUpJs2BNav4SyslKorGzoymJHjXiuG/SM+oUqbF1Mubm5rbqMnE4n2dnZba4rKSnhgQce4P/+7/9IT08PV3GEwKQoXHx2H3581Wiy0+0A1Ll8/HXtdtZ+UhrSRkS6EXhuTYNXssGKuBW2ADFmzBhKS0spKyvD5/NRWFjIlClTWl1z+PBh7rrrLv77v/+bQYMGhasoQrTSJ9PBT64aw4WjA12gBvDJ1nL+sPI79peH9g7S69eoqvPQ5JF9rkX8CVsXk8Vi4cEHH2Tp0qVomsaCBQsoKChg+fLlACxatIg//elP1NbW8vDDDwNgNpt5/fXXw1UkIYKsFhOzLxzI8AHpvL5xD7WNPqrqPSxbXczE0blMPS+fBKu5Q/fSDahv8tPkUUnplRTmkgvRfRQjxt72VFQ34YzjfkIZg+h+xfuqeHvzAaqPW3GdbLfQu1ciqmaQnpLA+OHZFOSlnfZeGRkOGuvdOOzWDgcYgK37qvio6AiVtW6y0uxcfHYfRg/K7FR9utraT0v54OtDuDwqjkQLl57bj9kTB0a6WF1OxiDaClsLQohYsOtgLe9/fZjEBAsZCtQ1+tB0g0a3SqO7EXuCGb9m8M7nZQAdChI+VcfX4MVmMZFst2I7TaDYuq+K1zbuDR47a9zB40gHibWflrL241IgMMmksckfPI7HICFak1Qbokf7oqQi+DjRZiEr3R7M5QSBVdgVNU243H62bHeGdG+fqlPd4KW63oPXf/IB8I+KjoR0vjt98PWhkM6L+CItCNGjnZjIz6QomEyAAgomVE3HaJ6x1FSmUlbRQH52aE3141sU9gQLCTYzpuOmfFfWutt9XmVtaCu+w+FkGzC5Orkxk4gt0oIQPVp6SkKbc2aTCavZTFZaIqkOGy1/y/2qzp/fLGblB3to6ESmV5+qU+fyUVnjpqbBi9cXaFVkpdnbvT4rLTHk79HVku3tLyJ0nOS8iC8SIESPNn5427U5SYkWHIkWFEUh2W4lO81Oou3YOMJXOyt58p/fsvGbQ53K9GoQmB5b0+ilzuXjojG57V538dl9Qr53V7v03H4hnRfxxfzQQw89FOlChMLl9uOK4zz9drsNdxw336OtfpmpiaSnJFDb/I4+IzWBKd/LY1j/9OC53r0SmT6hP+OGZXG40oXLo6LpBnsO1fPNrkocdivZ6XYURQm5fqqmk5pko2/vJGobvLi9Gtnpdmac3z/iA9QAQ/PTQIHDVS5UTcdhtzJtQv+4HKB2OBJoiuO/LQ5H29by6cg01ygTjdNAu1Ks10/TDT7f5uS9Lw/i9qrB8/2yHEw/rz8Tzu7b6folWM0kJVpCmh7bnXrCNNB4r1+oZJBaiBCYTQoTR+dyTkFvPvj6EJ9sLUfTDQ5Vunjure18XFzOlHP7BffIDoXXr+H1a1hMCvZEC3abBZOp/fxlQnQHCRBCdII9wcKVFwzg/JE5/OuLMr7dXQXAjv017Nhfw/D+6Vw+rh/9skIPFKpu0NDkp6HJf9KZT0J0BwkQQpyBjNRErp1SwCVj+/LuljJ2HKgFoORADSUHahgxIJ3LvtePvE4ECmieIqv6UFxgs5pJtJklWIhuIwFCiC7QJ9PBjTOGU93k5/UNu9h7uB6A7ftr2L6/hoK8Xkw+px+D+qScNO39qbTMfPL6NQkWottIgBCiC52Vl8bS2SPZd6Se9V8eDAaKXQfr2HWwjv45yUw6uy8jBqR3enyhvWCRYDWTYDNhNsnMddF1JEAIEQaD+qQGA8XGbw6zsyzQ9XTA2cgr/9pJZmoiF43J5XtDs06bq+lUjg8WNIHFrAQDhs1i6lRrRYgWEiCECKNBfVIZ1CeVw0ddbPzmEFv3VWMYUFXvYfXHpby7pYzxw7O5YGQOGalnvnJa1QxUTaXJo6IQSGueYDNjs5ixWqR1IUIjAUKIbtC3t4NFVwylut7Dx1vL+bKkAp+q4/FpfFR0hI+LjjCsfzoTRmYzNC+tS6a3GrQMcuuAH5NJIcFiwmY1Y7NKd5Q4PQkQQnSjjNRE5lw4kCvG5fFFSQWfFpdT2+jD4NjMp7RkG+OHZzNuaBa9kkNf/Xoyum7g9mm4m3NAmUwKNosJq8WExRz4LAPe4ngxt5La59dwVtRjGGAYBroBumEEjwFOVaP2qhu8F2DoBpH8D4n1lcanI/VrTdcNSg7UsLnYye5Dda2+pgBn5fVi3LBsRgxI75YuIotJwWoxYW3ukjr+e/aElcbxXr9QxVwLIjC9L7zFNpoDTiDwNAch/VjwMTCCQcgwAgHFaP6n5XHgeccCV+B5QrRmMimMHJjByIEZVNV5+Hy7ky93VtLkUTE4Nvsp0WZm9OBMzi3ozYDclLC901d1A/W4VoYCmM0KVrOJREcCXp+G2axgMUv3VE8Qcy0IIKajfEvQOb4FxHGtoMxMB5VHG1E1HVXVUfWY+/GckrQgTk/VdLbvr+HLHZXsOljbpkXcy2FjzJBMzh6SSb/ejm6bqXR83RQFLCYTFrOC2WzCajbFfOCQFkRbMdeCiHUmRYFTvKCTk2y4j8u1rxsGmmYca83ogRaMrrcOMie2eqTFErssZhNjBmcyZnAmdS4f3+yq5OtdR6moCWwsVOfy8VHRET4qOkJGSgKjBmUwenAG/bKSu20MwTDAr+kENso7tlteS4vD3Bw8LGYTZlPgs+SVij3SgogyXfkupnXQaHvcXmum1fFx9+gq0oLoHMMwOFLVxDe7jvLd3irqXG3TUqc6bIwYkM6IAekM7pva5e/mz7RuigJmRcFkUlAUBZMS6GIzm5Tmz6bg40iQFkRb0oKIY4qioChg4sxfcKcNJu0cQ9uxGEUJvMuMuXclEaYoCn17O+jb28GMC/pT5mykaE8VW/dV0dAU2H+i3uXjs21OPtvmxGYxcVZeL4blpzE0P61LZ0N1lmGAGnhncsrrFI4FjlbBw6wEurRkem63kQAhOuR0XWMdldU7GWtwtlnz5+A/ge6zwNcIfm4JMrpO82fj2OeWxz0o4pgUhQG5KQzITWHWhQM4WNFI8b5qivdVU928x7ZP1dlWWsO20hoActLtFOSlcVZeLwbmppzR6u1wMwjsu6Gd5IfaKoCYTZgUmlskSvBXVGk+p3DsjVLgzUnzO5TjvxmB3ztV01G1YzsEtvwOnvhrH2yNNz+/5Xe25d5K8zXHT3I5Xkt5jwU/JWpXvEuAEBHT8qJQgv+0ehCSdrvToE3gCVzb3vMDV7XMWGsJQFpzEIrW+GNSFPrnpNA/J4UZ5/enotZNSXOCwLKKxmBdnTVunDVuPvruCGaTQn5OMoP7pDK4byr52Skxtcq6VQDpxJavJ6OZzFTXebrsfqFQlMDP8viXQcuREgyAoLS0qpoDosl0fABUWr16uiLoSIAQcaEru9Paowff0R7rTlM1A03T8Ws6qhb5EKIoCjnpSeSkJzH5nH40eVR2H6plZ1ktu8rqaGjeClXTDUqPNFB6pIENXx3CYlbIy0pmQG4KA3MDwcaeIH8aupNhgNbmnUvX/k7JGIQQYWI6zeCppgfSZlijaJpnUqKFs4f05uwhvTEMA2eNmz2H6th9sI7S8oZAgj8C+ZtKyxsoLW9gY/Nzs9Ls9M9Jpn92MnnZyWSnJ0WuIiJiJEAI0QXMJhOORBNZGUloPj8+v45f1fCrOloUdFEpikJuRhK5GUlcNKYPmm5w+KiLvYfrKD3SwH5nAx7fsemqlbVuKmvdfLmjEgCr2UT/Pilkp9np1zxYnpVmxyxTV+OaBAghupjFbGqeYnrs5dUyruFvTtDn82sRDRpmk0J+djL52clMPicwduOsbmJ/eQMHnI0cqGigut4bvN6v6ew5WMeeg8fSgVjMCtnpSfTJSCI3M4mcjCRy0u0k261RO+gqQiMBQohuYDIpmAgsGLMnWNB1A49Pxe3V8GtdN9Da6fIpCn0yHfTJdHDBqMC5RrefQ5WNlFU0cqjSxeEqV3BKLQS6pg4fdXH4aOu1EUmJFrLT7WSn2ckKfiTSKzlBkgHGGAkQQkSAyaSQlGglKdGKquk0eVW8Pu2kUzsjIdluZVj/dIb1TwcgPT2J0rIaDh91ceioi/LqJo5UNVHT4G31vCaPGhwEP57FrNC7l53M1EQyeyWQmZpIRmoiGakJpDoSpLsqCkmAECLCLGYTqUk2SArkYfL5Nbx+PeLdUCdSFIVeyQn0Sk5gxMCM4HmPT8VZ7aa8uglnTRPO6sD4RaPb3+r5qmZQXt1EeXVTm3ubFIW0FBtpyQmkpyQEP/dKDpzr5bDFdJ6nWCUBQogo0jJ+kZQYGBfwNo9X+DUdTYv8YHd7Em2W4MK94zV51OBg99E6N5W1Ho7Weaiu97RpKemGQXW9t9W4x4kciRZ6OWykOmykJAU+pyZZSUmykdzy2W6RldZdSAKEEFHKpCjYEyyt1iS0rPb1q8c+ojFoQGAsor3AoesG9U0+jtZ5qGnwUl0fCBq1jT5qGrxtWh4tXB4Vl0flcFXbFsjx7AkWku1Wku1WHHYLjsTA46REC45ES6BrL8FCUmLgw2qWvbtPRgKEEDGkpYWRaAsctyzY86t6cysjEECiaCijDZNJIS050I3UHp+qUdvoo67RS12jj9pGL3UuH/UuH3UuH3WNvuAajva4vSpub6D10hEWcyAQJyfZsFlM2G0W7AmBfWcSE8wk2pof28zBjwSrhQTrsf2+4zVTrQQIIWKYoihYLUqbVBmargdXeqtac56h5rQh0c5mMZOdFpgFdTI+v0ZDk586l49Gt4+GJj8NTX4a3a0/XG7/aQf+Vc0IPr+zrBYTCc17fQc+m7G17P993GerxYTtuN36jv9o2fbV2vzZ0vzY0rzXRiRmgIU1QGzatIlHH30UXddZuHAhP/rRj1p93TAMHn30UTZu3EhiYiK/+c1vGDVqVDiLJETc2rqvio+KjlBZ6yYrzc7FZ/dh9KBM/lq4jS3bK/BrOhaTwvgR2Ywfls0n3x2hss5DZmoi543IpqBfGiUHqtlSUkF1vZf0lATGD8/mYGUjn29z0uTTSLKZmTAyh8vOzTtpOXYdrOWLkgpqGo7dA2hzriAvrcPPL8hL4/2vDwbK4VVJSrCcshwt96iu9wQG1Qekk5VmZ/ehOkr219DQ5MdmNZGRkojVYsLtU/GrBg1NPtzNM8pCCaUt3X10rNHSKS37agT32Wj+bGlOWmg2K1ias96aTcf25GhJajhmWE7I3zNs+0Fomsb06dP529/+Rk5ODtdccw1PPvkkZ511VvCajRs38tJLL/HMM8/w7bff8uijj7JixYrT3jvec7ZL/WJXpOq3dV8Vr23c2+a8I9HC9uaMri0MAuMDWSe8Qx83LCu4crrlz0JtoxeX29+cKVUJnp82oT9XjM9v3rjqWILEkrIa1n1W1uq+Hp+KAiScsFXw9An5bYLEroO1vPN56+cD5GbY+W5PVZvzl36vX5sgcbJ7jBqUTvG+mjbnW8px/H4XO5vrcXxGVkM3GD04g7TkBDw+Da+/+aP5sV/Vg8c+VQ98NJ+PBmv+Z17IzwlbC6KoqIgBAwaQn58PwKxZs1i/fn2rALF+/Xrmz5+Poiicc8451NfXU1FRQXZ2driKJURc+qjoSLvnS/a3/YMI4Paobc598PUhUpICgxstg7Yut4phBLKIHj+Q+2lxOQsmD2lzj60b9wS7u1qCidsb+F6O5p0SW96RFu2pYvywY691w4Bvdx9tpz/faDc4AGzZXsHMCwa2StH79c7KdjPTb9leESzD8b7eWcnoQZkk2gLdQADf7DoaeCd+QvLHOpePuRcNarcsJ6MbBv7macu+5paGTz32uM2HFggqLV2Dx8aXjh23dBm2dCFqzcdq83XBLsYz7FIMW4BwOp3k5uYGj3NycigqKjrlNbm5uTidztMGiM5kJYwlUr/YFon61TT62k3ZrRvtbOPRvB3tide7PCoZqYknPL95r4OW1OzNn5s8arv1bK8cug4otNmDosHtp39e+gnnVBJtbfeqUHUDm6XteY9P46yBma3r4dXazUZ7tM5D73bGNZp8GkMH9wYgs5e9+R7b2r2H26cxbEhWm/PxKmwBor2eqxOnknXkGiHE6T35s8mRLgJw5uXoinpEyz3iQdhWlOTm5lJeXh48bq9lcOI15eXl0r0khBBRImwBYsyYMZSWllJWVobP56OwsJApU6a0umbKlCm8+eabGIbBN998Q0pKigQIIYSIEmHrYrJYLDz44IMsXboUTdNYsGABBQUFLF++HIBFixYxefJkNm7cyNSpU7Hb7Tz22GPhKo4QQogQhW2aqxBCiNgmWa2EEEK0SwKEEEKIdkV1Liav18uSJUvw+XzBldk//elPqa2t5ec//zmHDh2iX79+/O///i+9evWKdHE7pWV8Jicnh6effjqu6jZlyhQcDgcmkwmz2czrr78eV/Wrr6/ngQceYOfOnSiKwmOPPcagQYPion579+7l5z//efC4rKyMn/70p8yfPz8u6vf888+zYsUKFEVh6NChPP7447jd7rioG8ALL7zAihUrMAyDhQsX8sMf/rBTr72obkHYbDZeeOEFVq9ezZtvvsmHH37IN998w7Jly5g4cSLvvvsuEydOZNmyZZEuaqe9+OKLDBlybEVqPNUNAr+oq1at4vXXXwfiq36PPvookyZNYt26daxatYohQ4bETf0GDx7MqlWrgj87u93O1KlT46J+TqeTF198kddee421a9eiaRqFhYVxUTeAnTt3smLFClasWMGqVav44IMPKC0t7VT9ojpAKIqCw+EAQFVVVFVFUZRgig6A+fPn895770WymJ1WXl7OBx98wDXXXBM8Fy91O5l4qV9jYyNbtmwJ/uxsNhupqalxU7/jffrpp+Tn59OvX7+4qZ+maXg8HlRVxePxkJ2dHTd127NnD2PHjsVut2OxWDjvvPP417/+1an6RXWAgMAPct68eVx44YVceOGFjB07lqqqquB6iezsbKqrqyNcys557LHH+Ld/+zdMx+2AFS91a3HLLbdw9dVX889//hOIn/qVlZWRkZHB/fffz/z58/nlL39JU1NT3NTveIWFhcyePRuIj59fTk4ON998M5dddhkXX3wxycnJXHzxxXFRN4ChQ4fyxRdfUFNTg9vtZtOmTZSXl3eqflEfIMxmM6tWrWLjxo0UFRWxc+fOSBepS7z//vtkZGQwevToSBclbJYvX84bb7zBM888wyuvvMKWLVsiXaQuo6oq27ZtY9GiRbz55pvY7faY7ZI4FZ/Px4YNG5gxY0aki9Jl6urqWL9+PevXr+fDDz/E7XazatWqSBerywwZMoSlS5dy8803s3TpUoYNG4bZ3DaPVUdEfYBokZqayvnnn8+HH35IZmYmFRUVAFRUVJCRkXGaZ0efr776ig0bNjBlyhTuvvtuNm/ezD333BMXdWuRkxPIP5+ZmcnUqVMpKiqKm/rl5uaSm5vL2LFjAZgxYwbbtm2Lm/q12LRpE6NGjaJ37+ZkdnFQv08++YS8vDwyMjKwWq1MmzaNr7/+Oi7q1mLhwoW88cYbvPLKK6SlpTFgwIBO1S+qA0R1dTX19fUAeDwePvnkEwYPHhxM0QHw5ptvcvnll0eymJ3yi1/8gk2bNrFhwwaefPJJLrjgAn7729/GRd0AmpqaaGxsDD7++OOPKSgoiJv6ZWVlkZuby969gT0YPv30U4YMGRI39WtRWFjIrFmzgsfxUL++ffvy7bff4na7MQwjLn92VVWB9OiHDx/m3XffZfbs2Z2qX1SvpC4pKeG+++5D0zQMw2DGjBnceeed1NTU8LOf/YwjR47Qp08ffv/735OW1v7uVLHgs88+47nnnuPpp5+Om7qVlZXxk5/8BAiMI82ePZs77rgjbuoHsH37dn75y1/i9/vJz8/n8ccfR9f1uKmf2+3m0ksv5b333iMlJZDaO15+fk899RRvvfUWFouFESNG8Oijj+JyueKibgCLFy+mtrYWi8XC/fffz8SJEzv1s4vqACGEECJyorqLSQghRORIgBBCCNEuCRBCCCHaJQFCCCFEuyRACCGEaFdUZ3MV4lQWLlyIz+fD7/dTWlpKQUEBACNHjuTxxx+PcOk6pri4mLKysrhaqSzih0xzFTHv4MGDLFiwgM8++yzSRWlDVVUslpO/D1uxYgWffPIJv/vd77r83kKcKfntEnFp5cqV/OMf/0DTNFJTU3n44YcZOHAgK1asYN26dTgcDnbu3EmfPn34j//4D5544gnKysoYO3YsTzzxBIqicM8992C32zlw4ADl5eWcf/75/OpXv8JqtdLQ0MBjjz3Grl278Hq9XHjhhdx7772YTCYWLVrEhAkT+Prrr0lKSuKpp54KLhL0er2MHTuWhx9+mPr6ev70pz/hcrmYN28e559/PkuWLGHx4sV8/PHHAOzfvz94vH//fhYtWsS1117L5s2bufrqq5k3bx5PPvkkX3zxBT6fjxEjRvDQQw9ht9sj/BMQccEQIsaVlZUZEyZMCB5v3rzZuO222wyv12sYhmGsX7/eWLJkiWEYhvHqq68aEyZMMMrLyw3DMIybb77ZmD9/vtHQ0GD4fD5j5syZxubNmw3DMIxf/OIXxrx58wyXy2X4fD7jhhtuMP7+978bhmEY9957r7FmzRrDMAxD0zTjpz/9qbFy5UrDMAzjuuuuM3784x8bqqoGv15bWxt8fPfddxuvvvpqsDw/+9nPgmUvLS01LrzwwnaPS0tLjaFDhxrr1q0Lfv2pp54ynn766eDx448/bvz+978/s/9QIZpJC0LEnQ0bNrBt2zYWLlwIgGEYuFyu4NfHjRsXTCQ4cuRIPB4PycnJAAwbNowDBw5w/vnnAzBz5kySkpKAQA79Dz74gEWLFvH+++9TXFzMM888AwRyhfXv3z/4PebMmRPMoKnrOsuWLeOjjz5C13Vqa2s7vVNZUlIS06dPb1VXt9tNYWEhEMi+OmrUqE7dW4gTSYAQcccwDL7//e9z5513tvv1hISE4GOTydTmWFXVk95XURQg8Ef/6aefpm/fvu1e2xJUAFatWkVRURF///vfcTgc/PGPf+TIkSPtPs9sNqPrevDY6/We9L4tZfr1r3/Neeed1+79hDgTMs1VxJ2WrJVOpxMIJAvcunVrp+719ttv43a78fv9rFmzJtiymDJlCsuWLUPTNCCQebisrKzdezQ0NJCeno7D4aCuri74bh/A4XDQ0NAQPM7Ozsbj8QTvtXbt2tPW9bnnngsGksbGRvbs2dOpugpxIgkQIu5ccMEF3Hnnndx2223MnTuXOXPm8MEHH3TqXuPGjeOOO+5g9uzZ5OfnB7cY/dWvfoWu68ybN485c+Zw6623UllZ2e49rrrqKmpra5k9ezZ33313q3f7F110EQ0NDcydO5fHHnsMm83Gfffdx4033sj111+P1Wo9Zfluv/12hgwZwjXXXMOcOXNYsmQJ+/bt61RdhTiRTHMV4iTuuecexo0bx6JFiyJdFCEiQloQQggh2iUtCCGEEO2SFoQQQoh2SYAQQgjRLgkQQggh2iUBQgghRLskQAghhGjX/w/VTvdCkDkrbgAAAABJRU5ErkJggg==\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/gitlab/moocrr-session1/moocrr-reproducibility-study/raw/5c9dbef11b4d7638b7ddf2ea71026e7bf00fcfb0/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
}