"
+ ]
+ },
+ "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:
Thu, 21 Sep 2023
Deviance:
3.0144
\n",
+ "
\n",
+ "
\n",
+ "
Time:
13:47:36
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: Thu, 21 Sep 2023 Deviance: 3.0144\n",
+ "Time: 13:47:36 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:
Thu, 21 Sep 2023
Deviance:
18.086
\n",
+ "
\n",
+ "
\n",
+ "
Time:
13:47:36
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: Thu, 21 Sep 2023 Deviance: 18.086\n",
+ "Time: 13:47:36 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/XmS2ZTPaQBUgIW9gRLYiiIoqyyK5IFahaFau22m9rvVe9tV61V63312tvbXtb0VrX0gouLFG0goIbimskEPZAWDIJ2TOZ9Zzz++MkIyEBkpDJLHk/H488MufMmZPPh5B5z2d7fxRd13WEEEKIE5jCXQAhhBCRSQKEEEKIdkmAEEII0S4JEEIIIdolAUIIIUS7JEAIIYRoV8gCxH333cekSZOYM2dOu8/rus5//dd/MW3aNObOnUtxcXGoiiKEEKILQhYgrrrqKp555pmTPr9582ZKS0t55513+PWvf82DDz4YqqIIIYTogpAFiHPPPZeUlJSTPr9hwwYWLFiAoiicffbZ1NfXU1FREariCCGE6KSwjUE4nU5ycnKCxzk5OTidztO+zudX8QdU/AENf0AjoGqoqoaq6WiaLAoXQojuYgnXD24vw4eiKKd9XW2DF2dlw2mvUxRQmh+YABQwNd9fURRMivFdUVqOjccmk/HYZKL53OnL1J0yM5Oo7ED9opXUL3rFct2gd9Svs8IWIHJycigvLw8el5eXk5WV1W3313XQmx9o353t9H2U5sBiNilG8DAZj40vE2azEgw8QggRS8IWIKZOncpLL73E7Nmz+eabb0hKSurWANFddB1UXUc9RfeVSQGL2YTZpGCxmLA0Bw6LWWYRCyGiV8gCxF133cVnn31GTU0NF198MXfeeSeBQACAxYsXM2XKFDZt2sS0adOw2+08+uijoSpKyGk6+ALN7RSfGjyvYAQOi8WEzWLCajFJ0BBCRA0l2tJ9V1Q3dWgMIlKZFLBZzdgsZmzWtgGjN/SDSv2iUyzXDXpH/TorbF1MvZWmg8en4mluaVhMCnE2M/E2M1aLOcylE0KI70iACLOAphPwBHB5AphNCja7DX9Aw2qRrighRHhJgIggqqbT6PZTXe/BajaREG8h3mbu8am2QggBEiAill/VqHP5aGiChHgrCXEWTCYJFEKIniMBIsJpOjS6/bg8fhLiLDjirRIohBA9QgJElNB1cHkCNHkDJNqNFoV0PQkhQklGQqOMrkNDk5+qOg/e49ZcCCFEd5MAEaUCmk5No5eaBi+qpp3+BUII0UkSIKKc169yrM5DkycQ7qIIIWKMBIgYoOtQ3+Sjut4jrQkhRLeRABFDfAGNqjoPbq+0JoQQZ04CRIzRdKhz+ahr9La754YQQnSUBIgY5fapVNfLALYQouskQMQwv2p0Ofn8Mh1WCNF5EiBinKZDTYNXxiWEEJ0mAaIX0DHGJZo8/nAXRQgRRSRA9CL1TX4a3RIkhBAdIwGil2l0+6lv8oW7GEKIKCABohdq8gSkJSGEOC0JEL1Uo1u6m4QQpyYBohdrdPtl4FoIcVISIHq5+iY/Hp9MgRVCtCUBQlDn8uEPyIprIURrEiAEug41jZKWQwjRmgQIAYCm6dQ2+CTBnxAiSAKECPKrGvUuWSMhhDBIgBCtuH2q5G0SQgASIEQ76pt8BFQZjxCit5MAIdrQdahtkA2HhOjtoi5APPTMFr7efUzevEIsoOnUN8kiOiF6M0u4C9BZB50NHHQ28NXuShZMHkRaUny4ixSz3N4A8VYzcTZzuIsihAiDqGtBWMwKALsP1fG/K4v46NujaNKaCJm6Jp/8+wrRS0VdgHhw2SQG5iQB4A9oFH5ygKfXbudYnTvMJYtNmqbTIFNfheiVoi5A9O3jYNncUSyYPIg4q9H1caC8gT+s+lZaEyHi9ql4fbKvtRC9TdQFCACTojBxZDb/b9FZFOSmAMYir8JPDvBs4Q5qG71hLmHsqWvyoWkSfIXoTUIaIDZv3syMGTOYNm0ay5cvb/N8Q0MDt912G/PmzWP27Nm8+uqrnbp/amIcP7xiBFddPDjYmth3pJ4nVxXx9R6Z6dSdNE2XneiE6GVCFiBUVeXhhx/mmWeeobCwkHXr1rFnz55W17z88ssMGTKENWvW8OKLL/L444/j83XuTUhRFCaMyOKnV5/FwL7G2ITHp/LKxj288t4eSWXdjTyyylqIXiVkAaKoqIj8/Hzy8vKw2WzMnj2bDRs2tLpGURRcLhe6ruNyuUhJScFi6drM27SkOJbNHsUV5w3AbDJmOn2zp4o/vvotZRWNZ1wfYWho8knWVyF6iZCtg3A6neTk5ASPs7OzKSoqanXN0qVLuf3225k8eTIul4vf/e53mEynj1np6Y6TPjf/0gLGj8rhmTXbKK9qorrBy1NrilkwZQiXTxyASVG6Xqkecqr6RQKrzUxGir3Lr8/MTOrG0kSeWK5fLNcNYr9+nRWyANFe/79ywpvzhx9+yMiRI3nhhRc4ePAgN954IxMmTCAxMfGU966udp3y+QSridvmjabwkwNsLalA03Ree28P2/ce4+pLhpIQH7nrA9PTHaetXyRwNXiwx3X+3zEzM4nKyoYQlCgyxHL9Yrlu0Dvq11kh62LKycmhvLw8eOx0OsnKymp1zWuvvcb06dNRFIX8/Hxyc3PZt29ft/x8m9XMlRcPZvHlBcEB7JKDtfzxtSLKKmL3P0FPaXRLGg4hYl3IAsTYsWMpLS2lrKwMn89HYWEhU6dObXVN3759+eSTTwA4duwY+/fvJzc3t3vLMTiDOxaOpV9GAgC1jT6Wr9nOluJymeV0BlRNlwFrIWJcyPpaLBYLDzzwAMuWLUNVVRYuXEhBQQErVqwAYPHixfz4xz/mvvvuY+7cuei6zt133016enq3lyUjOZ5b54/hzS0H+HS7E1XTWfNRKWUVjSyYPBirJSqXg4Sdy+3vUjeTECI6KHqUfYyuqG7CeQb9hF/vPsbrm/fhb97voF9GAkunD4uYpH/RMgbRIsVh61SQ6A39vLFav1iuG/SO+nVWr/vofHZBH25bMJr0pDgAjlQ18afXtrHvSF2YSxadXDIWIUTM6nUBAqBvhoOfXDWWYXmpADR5AzxbWMKW7eWneaU4UUDGIoSIWb0yQADY4yxcP2M4F4/rB4Cm66z5sJQ3PtgnC8E6SVoRQsSmXhsgAEwmhZnnDeD7U4cG95n4bEcFz71VIp+KO0FaEULEpl4dIFqcPbQPP5o3muQEKwB7D9fz5ze2UVXvCXPJooe0IoSIPRIgmuVmJnL7ld+tlzhW5+HPr2+jtLw+zCWLDtKKECL2SIA4TorDxo/mjWbUwDSgZfB6B9/uqwpzyaKDtCKEiC0SIE5gs5pZMm0Yk8/qC0BA1fnHu7v5sOiorLw+DWlFCBFbJEC0w6QoXHF+PnMvHIiigA68ueUAhZ8ckC1NT0NaEULEDgkQpzBpdA5Lpw3Dajb+mT7eVs4rG/cQUGUa7MkENF02aRIiRkiAOI1RA9NZNnckCc3pJIr2VvHC+p14fWqYSxa5mjwSIISIBRIgOiAvK4lb548mNdEGwJ7DdTy9brukvD4JX0CTVpYQMUACRAdlptq5df4YstOMndSOHHOxfE0xtY3eMJcsMrmkFSFE1JMA0Qkt02Dzc4ysiMfqPDy1uphjte4wlyzyeHwBGdAXIspJgOgke5yFG2eNCCb6q3P5eGpNMUeORU+K7p6g68iUVyGinASILrBZzPxg+jDGDs4AjO6UZ9Ztl61MTyCD1UJENwkQXWQxm7hm6lDOHWHss+3xqfy1cAf7jkhqjhaqpstsLyGimASIM2AyKSyYPIgLx+QA4PNrPP9WCbvKasNcssjRJN1MQkQtCRBnSFEUZk3K55Jz+gPgVzVefHsnJQdrwlyyyODzq2iaDFYLEY0kQHQDRVGYfm4e08/NA4yulZff2cWO0uowlyz8dMAtK6uFiEoSILrRJef054rzBgDNQeJfuyneL0HCLYPVQkQlCRDdbPK4fsyelA8Y25iueHc323p5uvCApuMPyGC1ENFGAkQIXDi2L3MvGAgYQeIfG/b0+pZEk1cChBDRRgJEiEwak8O8CwcC37UktvfiMQmPLyD7aQgRZSRAhND5o3NatSRWvLu71w5c67qxVkQIET0kQITYpDE5zLnAGJNQNZ2/v7u7166TkNQbQkQXCRA94IIxfYMD16qm89I7O9lzuC7Mpep5voCGKmnAhYgaEiB6yIVj+zJjorFOIqDqvPj2TvYf7X1pOdzSzSRE1JAA0YOmnN2fy8bnAuAPaDy/vqTXJfjzSDeTEFFDAkQPm/q9/kw5ux9g5G7625slHK3qPanCvZJ6Q4ioIQGih7Wk5WhJ8OfxqTxbuIOKXrTpkMxmEiI6SIAIg5YEfxOaU4W7PAGeLdxBdb0nzCXrGV6/BAghooEEiDBRFIUFFw3irCHGpkP1Lh/PFu6grhfsce3zq7IdqRBRQAJEGJlMCosuHcLI/DQAqhu8/P6fX8X8Tmw6yEZCQkSBDgWIm2++mffee6/TqRI2b97MjBkzmDZtGsuXL2/3mk8//ZT58+cze/ZsfvCDH3Tq/rHAbDJx7WUFDO6XDMCRShfPvbUj5t9AZRxCiMjXoQBxzTXX8Pzzz3P55ZezfPlyampOvxmOqqo8/PDDPPPMMxQWFrJu3Tr27NnT6pr6+noeeugh/vznP1NYWMjvf//7rtUiylktJq6bPpy8rEQADlW6ePGdnfgDsbuoTLqZhIh8HQoQ06dP57nnnuPpp5+moqKCOXPm8O///u9s27btpK8pKioiPz+fvLw8bDYbs2fPZsOGDa2uWbt2LdOmTaNfP2PaZ0ZGxhlUJbrF2czcMHME/TIdAOw7Us8rG/fE7JRQ6WYSIvJZuvIiq9VKXFwc99xzD5MnT+bee+9tc43T6SQnJyd4nJ2dTVFRUatrSktLCQQCXHfddbhcLq6//noWLFhw2p+fnu7oSrEjXjrw/645h//vpS84VuumuLSat7aW8YOZI1AUJdzF6zYtv794m5mMFHuYS9P9MjOTwl2EkInlukHs16+zOhQg3nnnHV566SWqqqpYsmQJhYWFOBwOAoEA06dPbzdAtDdeceKbnKqqFBcX89xzz+HxeLj22msZN24cgwYNOmV5qqtjd2FZerqDG2YO56nVxTS6/Xz0zRHMwMzmneqiXXq6I/j7U4CA148phoJfZmYSlZWxuTo+lusGvaN+ndWhALFq1SpuueUWJk+e3PrFFgv3339/u6/JycmhvLw8eOx0OsnKympzTVpaGgkJCSQkJDBhwgRKSkpOGyBiXUZyPDfOGsHTa7fj8als/uYIDruFyWf1C3fRulVLN5M9rksNWSFEiHVoDOKpp55qExxaTJ06td3zY8eOpbS0lLKyMnw+H4WFhW2uveyyy/j8888JBAK43W6KiooYMmRIJ6sQm/pmOLh+5nAsZuPT9VtbDvLlrsowl6r7yWwmISJXhwLEkiVLqKv7Lj11bW0tS5cuPeVrLBYLDzzwAMuWLWPWrFlcccUVFBQUsGLFClasWAHAkCFDmDx5MvPmzWPRokVcffXVDBs27AyqE1sG5iSz+PJhmJp7YF7btJedB08/gyyayGwmISKXondgccP8+fNZvXr1ac/1hIrqJpwx3E94fB99iy92VvDqpn0AWM0mbp4zkgHZ0TmY1l79Uhy2mOlmiuV+7FiuG/SO+nVWh1oQmqbR1NQUPHa5XKiqdA30lPHDs5g50Rik9qsaz6/fSUVN7CT3k24mISJThwLEnDlzuOmmm1i9ejWrV6/m5ptvZt68eaEumzjO5HF9uWhsX8DYuvNvb8ZO3iafpAAXIiJ1qF1/6623kpWVxcaNG9F1nWuvvbZD6xVE91EUhZnnD6DR7efrPceoc/n421sl3DpvdNR3z+gYrYiE+OiuhxCxpsN/kVdeeSVXXnllKMsiTsOkKFw1ZTAuj5/dh+qoqHHzwts7uWnWSKyW6M676PEFJEAIEWE69BdZVVXFiy++SFlZGYHAd5lGe2vupHCymE0suXwYz6zbzuFjLg6UN/DPjbtZcvkwTKboXXDmC2gEVA2LOboDnRCxpEMB4s4772TIkCFMmjQJs9kc6jKJ04izmbnhihH8ZfU2quu9bC+tYe3Hpcy7cGBUp+Ro8gZITrCFuxhCiGYdChD19fX8+te/DnVZRCck2q3cOGskf3ljGy5PgE+3O0lx2LjknP7hLlqXebwBkuzWqA5yQsSSDrXnCwoKcDqdoS6L6KSM5HhuuGIEtubxh3e2lvHFzoowl6rrNF2mvAoRSTrcgpg3bx7nnHMOcXFxwfMyBhF+uZmJLJk2jBfW70TTdV7fvI9Eu5XhA9LCXbQucXsDUT8rS4hY0aG/xDlz5jBnzpxQl0V00bC8VK6aMphV7+9F02HFu7tZNncUuZmJ4S5ap8lgtRCRo0MBQqa3Rr7vDcuk3uXjna1l+ALGauvb5o8mIzk+3EXrNLc3QJIMVgsRdh36mFZaWsrixYuD2ViLi4v5wx/+ENKCic6bcnY/zh+VDYDL7ee5N0todPvDXKrOc3sDnd7/XAjR/ToUIB588EFuv/12kpKMZE8jR45k/fr1IS2Y6DxFUZhzwUBGDTTGH6rqPbywvgSfP7oGfmWwWojI0KEA0dDQwMUXXxycfmgymbBarSEtmOgak0nhmqkF5Ddnez1U6WLFht2oUZbryOWJvpaPELGmQwHCbDbj9/uDAcLpdGIyySBipLJaTFw3Yxh9Uozxh50Ha1nz4f6o6rYJqHrUtXyEiDUd3jDojjvuoKamhj/84Q8sWbKEm266KdRlE2cgId7KjbNGkGQ3WnpbSyrY+OXhMJeqc1yewOkvEkKETIdmMS1YsIDc3Fzee+893G43jz/+OBMmTAh12cQZSksyFtItX1uMz6+x4YtDpDhsTBiRdfoXRwCvX5Upr0KEUYdXJE2YMEGCQhTq18fB0mnDeP4tYyHdGx/sIykhehbSSX4mIcKnQwFi4cKF7ebHWbVqVbcXSHS/gtxUFk4ZzMooXEjn9gZItFsxSX4mIXpchwLEPffcE3zs9XopLCwkKys6uimE4ZxhmdRF4UI6XTeS+CXEy6w5IXpahwLExIkTWx1fdNFFMkgdhaac3Y86l49PtzuDC+lunT+aRHtkv/k2SYAQIiy6NPrX2NhIWVlZd5dFhJiiKMyNwoV0AVXHH4jsMgoRizo9BqFpGocOHeLGG28MacFEaLQspPtr4XYOOhuDC+l+MH045gjeka7Jq5Jikc2qhOhJnR6DMJvN5Obmkp2dHbJCidCyWkxcP2M4T60pprLWw86Dtaz+YB9XXjw4Yjfr8fgCJCXIYLUQPalLYxAi+iXEW/nhFSP5y+ptNDT5+XxnJckOG5dPyAt30dplDFarJMTLXhFC9JQO/bWdf/757X6y1HUdRVH45JNPur1gIvTSkuL44RUjWL5mO16/ysYvD5OUYOO8UZHZOnR7AxIghOhBHfprW7x4MbW1tVxzzTXous6rr75KdnY2s2bNCnX5RIj1zXDwg+nDeO6tElRNZ81H+0lKsDJqYHq4i9aGX9XwBzSsFllZLURP6NBf2tatW/nP//xPRowYwciRI7n//vvZtGkT/fv3p3///qEuowixIf1TuPqSIYDRlfOPDbs5UN4Q5lK1z+2V/ExC9JQOBYiKigqqq6uDx9XV1VRWVoasUKLnjRvah9mT8gFjWunz60twVjeFuVRtuX0BtChLXS5EtOpQF9MNN9zA/PnzufTSSwHYtGkTt956a0gLJnrehWP7Uu/y8UHRUTw+lefeMhbSpSbGhbtoQbpu7BUhW5IKEXodChBLly5l/PjxbN26FV3XWbp0KcOHDw912UQYzDhvAI1uP1/tPkady8dzb5Xwo7mjI2pwuMkbwBFvxRTB6zaEiAUd/qvPzc1FVVVGjx4dyvKIMDMpCldNGUyj28/uQ3VU1Lh54e0Sbpo9EluELFTTdSNIRHqKECGiXYfGIDZt2sTs2bO58847Afj222+57bbbQlowET5mk4kl04aRm+kA4KCzkRXv7kbVtDCX7Dsujx8tinbIEyIadShAPPnkk6xatYrk5GQAxo4dy8GDB0NaMBFecVYzN1wxotW2pa9v3hcx25bqOjTJjnNChFSHJ5RnZma2OrbZZJAw1jnirdw4ayTJDuN3/eWuY7z16cGICRJN0ooQIqQ6FCAcDgfHjh0Lrqb+9NNPSUpKOu3rNm/ezIwZM5g2bRrLly8/6XVFRUWMHDmS9evXd7DYoqe0rLa2xxnjDx8WHWXzN0fCXCqDpsu6CCFCqUMB4he/+AW33HILhw4d4rrrruPuu+9ulcCvPaqq8vDDD/PMM89QWFjIunXr2LNnT7vX/fa3v+Wiiy7qWg1EyOWkJ3D9jBFYm/eGfvuzMraWVIS5VAaXJxAxLRohYk2HZjGNGzeOF154gS+//BKAc845JzgecTJFRUXk5+eTl2ckf5s9ezYbNmxg6NChra578cUXmTFjBt9++21Xyi96SH5OEkumFfDi27uCe1vbbWbGDM4Ia7k0TcctSfyECInT/lWpqsr3v/99Xn31VaZMmdLhGzudTnJycoLH2dnZFBUVtbnm3Xff5fnnn+9UgEhPd3T42mgUqfWblO7AYrPw7JpidB1eeW8PfTIcjBrUuSDR3fUzmxQyMyLn3ywz8/Tdr9EqlusGsV+/zjptgDCbzaSlpeH1eomL6/iK2vaa/SdmhH3kkUe4++67MZs7N7++utrVqeujSXq6I6LrNyQniTkXDmTtR6UEVJ0/v1rEzbNHMiC7Y39Yoaqfz+3DHhf+VkRmZhKVlZGZx+pMxXLdoHfUr7M69Bc1cOBAli5dyowZM0hISAieX7p06Ulfk5OTQ3l5efDY6XSSlZXV6ppt27Zx1113AVBTU8OmTZuwWCxcfvnlnaqE6FmTRufg9gZ49/ND+AMaz71Vwi1zR9E3jJ/iXR5/RAQIIWJJh/6iXC4XBQUF7Nu3r8M3Hjt2LKWlpZSVlZGdnU1hYSH/8z//0+qajRs3Bh/fe++9XHLJJRIcosSl5/TH41X58Fsjb9Pf3izhR3NH0SfVHpbyBFQdr08lzhYZq72FiAWnDBC/+c1vuPfee3nsscf46KOPuPDCCzt+Y4uFBx54gGXLlqGqKgsXLqSgoIAVK1YAxh4TInopisIV5w/A7Qvwxc5KGt1+/lq4gx/NG01aUniS+7k8fgkQQnQjRT/FHMErr7yS119/vc3jcKqobsIZw/2EkT4GcSJN0/nnxt18u89IB5+RHM8t80aRfJJsq6GuX1pSHHHW8AWJWO7HjuW6Qe+oX2edch3E8bFD5pqL9phMCosuHcrwAakAVNV7eLZwBy6PPyzlcbnD83OFiEWnDBA+n4+9e/eyZ8+eVo9bvoQAsJhNLLl8GIP7GWtjKmrc/K1wR1hWOfsCGl6/2uM/V4hYdMoupqlTp578hYrChg0bQlKoU5Eupsjl9av87c0dHHQ2ApCXlchNs0a2GhfoifrZLCbSk+ND+jNOJpa7KWK5btA76tdZpxykPn6WkRCnE2c188MrRvDXdTs4fMxFWUUjz68v4YdXjMDWg+MCLa2IcI5FCBELOpzNVYiOiLdZuHHWCHLSjfUypeUNvPD2TnyBnu32kbEIIc6crCwS3S4h3spNs0fy9NrtVNa62Xeknpfe3sV1M85sm9rdh2r5vKSCmgYvaUlxTBiRRUFuarvX+gIaHl+AeJv8F4822/ZX8WHRUSpr3WSm2rnorL6M6WQ6F9E9pAUhQiLRbuXmOSODGw7tOVzHS+/sxN/FlsTuQ7W8/VkZVfVeNB2q6r28/VkZuw/VnvQ1jU3Siog22/ZX8eqmfThr3Gg6OGvcvLppH9v2V4W7aL2SBAgRMskJNpbNGUVGc5DYfaiOv7z2Lf5A57cu/fwk6cVPdh4goOmy61yU+bDoaKfOi9CSACFCKtlhBIn0ZGN1dfG+quaWROeCRE2Dt1PnWzTKrnNRpbLWfZLznh4uiQAJEKIHpJwQJHYfquPFTg5cnyx9x+nSemjSiogqmSfJ5ZWZGp5py72dBAjRI1IT47hlziiy0ow3gD2H63hh/U58HVzUNmFEVqfOH8/l8aNqne/WEj3vorP6duq8CC0JEKLHpCTGcdeS8cGB631H6nnurRK8vtMHiYLcVGZMzCMjOQ6TAhnJccyYmHfSWUzH03VwuaUVEQ3GDMpg4ZTBZKfZMSkK2Wl2Fk4ZLLOYwuSUK6kjkaykjm7p6Q5KD9Xw13Xbg/3KeVmJ/PCKESHdz0EBMlLisZhD+5kollfjxnLdoHfUr7OkBSF6XHKCjVvmjg4upiuraOSZddtpDOHiNh1ZPCdEZ0mAEGGRaLeybM4o+mcau9AdrWri6bXbqXP5QvYz3T6VgCpjEUJ0lAQIETYJ8Zbm/awTAWOK4/I1xVTVh25KYyhbKULEGgkQIqyM3E0jGdo/BTDWNSxfXUx5dVNIfp7Hp3ZpoZ4QvZEECBF2cVYz188czqiBaQA0uP08vbaYg87QDBhKK0KIjpEAISKCxWxi8eXD+N6wPgC4vSp/XbeDnQdruv1nef1qh6bWCtHbSYAQEcNsUrhqyhAuHJsDgF/VePHtnXy1q7Lbf1ZDk0+20RXiNCRAiIhiUhRmnZ/PjIl5AGg6rHx/L5u/PtKtb+gBTacpDFuiChFNJECIiKMoClPO7s/CKYMxKca59Z8dZO1HpWha9wWJRre/W+8nRKyRACEi1vjhWfxg+nCsFuO/6ZbtTl7+165u251O142uJiFE+yRAiIg2Ij+NZXNG4Yg30nDsOFDDX9ft6LY3drdP7XDCQCF6GwkQIuLlZSVy24IxwY2Hyioa+cvqYpzdtFaizuWTPSOEaIcECBEVMpLjuW3+aPKzjYRjNQ1e/rK6+JRbjnaUquk0hDDFhxDRSgKEiBqOeGOf67OHGmslvH6V598q4ZPi8jOe4eT2qbhlVpMQrUiAEFHFYjax6NIhXDY+FzCmwa79qJTVH+4/40R8DU0+2VhIiONIgBBRR1EULhufy7WXDcViNubBfrajgmff3HHzuuRsAAAdzklEQVRGaTQ0HeoaZQGdEC0kQIioddaQPtw6bzTJDhsApUcb+L/Xv+VwZWOX7+kLaNQ3Sa4mIUAChIhy/TMT+cmVY8jLMlKG1zb6eGpNMV+eQXoOtzeAyyNBQggJECLqJSXYuGXuKM4dkQVAQNVZ9f7eMxqXaGjyS0I/0etJgBAxwWI2ceXFg7ly8iDMzfk5Pt3uZPmaYmoavF26Z63Li1cW0YleTAKEiCnnjszmR/NGkdI8LnGo0sUfXyuipAtpw3Udahu80pIQvVZIA8TmzZuZMWMG06ZNY/ny5W2eX7NmDXPnzmXu3Llce+21lJSUhLI4opfIy0rijoVjKcg1dqlze1VeWL+Tt7Yc6HSXkw7UNnrx+GSNhOh9QhYgVFXl4Ycf5plnnqGwsJB169axZ8+eVtfk5uby0ksvsXbtWm6//XZ+9atfhao4opdxxFu5YeYILhufS3NCWD4oOsryNcVUd3LPax1j+qsspBO9TcgCRFFREfn5+eTl5WGz2Zg9ezYbNmxodc33vvc9UlKMT3lnn3025eXloSqO6IVMJmO9xI2zR5JktwJGl9MfXv2Wr/cc69S9dIycTbJdqehNLKG6sdPpJCcnJ3icnZ1NUVHRSa9ftWoVF198cYfunZ7uOOPyRTKpX/eamO5gxOA+PF+4neJ9VXj9Kq9s3MP+8gYWTx9OQry1U/ezxFtIS4o/6fOZmUlnWuSIFct1g9ivX2eFLEC0txpVUZR2roQtW7awatUq/v73v3fo3tXVrjMqWyRLT3dI/UJk8WVD+Tgrkbc/O4iq6Wzd7mTXgRoWXTqEwf1SOnyfaqDS2khqoq3N/+nMzCQqKxu6ueSRIZbrBr2jfp0Vsi6mnJycVl1GTqeTrKysNteVlJRw//3383//93+kpaWFqjhCYFIULjqrLz++cgxZaXbA6Db667odrPu4tFMbEXn9KtX1XtmRTsS0kAWIsWPHUlpaSllZGT6fj8LCQqZOndrqmiNHjnDnnXfy3//93wwaNChURRGilb4ZDn5y5VguGGN0gerAx9vK+cOqbzlQ3vFPkH5Vo6rec8ZJAoWIVCHrYrJYLDzwwAMsW7YMVVVZuHAhBQUFrFixAoDFixfzpz/9idraWh566CEAzGYzr732WqiKJESQ1WJizgUDGZGfxmub9lLb6KOq3sPyNcVMGpPDtHPziLOaT3sfVdOprveQlGDDHheyPychwkLRoyx1ZUV1E84Y7ieUMYieV7y/ire2HKT6uBXXiXYLfVLiCag6aUlxTBiRRUFu6invYzWbGJyfTl1t53e627a/ig+LjlJZ6yYz1c5FZ/VlzKCMTt8nFNZ9Usr7Xx3G5QngiLdwyTn9mTNpYLiL1e1kDKIt+cgjerXdh2p576sjxMdZSFeM9Q6qptPoDtDobsQeZ8av6rz9WRnAKYOEX9WorHXjbvSSlGDDZGp/UsaJtu2v4tVN+4LHzhp38DjcQWLdJ6Ws+6gUMCaZNDb5g8exGCREa5JqQ/Rqn5dUBB/H2yxkptmDuZzAWIVdUdOEy+1n6w5nh+7p9qlU1rlxefwd2lviw6KjnTrfk97/6nCnzovYIgFC9GonJvIzKQomE5jNRgJAMHIy1bl87Cyro6yiY10Qum5khK2q89Do9p9yp7rKWvdJznduxXconGxhoEsWDPYKEiBEr5aWFNfmnNlkwmo2k5kaT7LDRstSB39A489vFLPq/b00NPk6dP+AptPo9lNZ66G63oPbG2jTqshMtbf72szUky/G6ymJ9vYXETpOcl7EFgkQolebMKLt2pyEeAuOeAuKopBot5KVaife9t2Mpi93VfLEP79h09eH8Qc6PsXVF9Coc/moqHVT7/IFX3vRWX3bvf5k53vSJef079R5EVvMDz744IPhLkRnuNx+XB389BaN7HYb7hhuvkda/TKS40lLigum9U5PjmPq93IZPiAteK5PSjwzJg5g/PBMjlS6cHkCqJrO3sP1fL27EofdSlaaHUVROlw/v6rh9gbw+lSyUhPISbdTXe/F7VXJSrMz87wBYR+gBhiWlwoKHKlyEVA1HHYr0ycOiMkBaocjjqYYfm9xONq2lk9HprlGmEicBtqdor1+qqbz2XYn735xqFV21/6ZDmacO4CJZ/XrUv0UBew2CwnxluDYR6TpDdNAY71+nSXTXIXoBLNJYdKYHM4u6MP7Xx3m423lqJrO4UoXz765g4+Ky5l6Tv/gHtkdpevQ5A3Q5A0QZzWTEG/BZjGdNH+ZED1BAoQQXWCPs3DF+fmcNyqbf31exjd7qgDYeaCGnQdqGDEgjcvG96d/ZucCBRh5nrx+FZMCcVYzcTYzcVazBAvR4yRACHEG0pPjuWZqAReP68c7W8vYebAWgJKDNZQcrGFkfhqXfq8/uV0IFJpurKlw+1QUxVinEd8cLIToCRIghOgGfTMc3DBzBNVNfl7buJt9R+oB2HGghh0HaijITWHK2f0Z1DepSy0BXQe3N4DbG0BRwGYxAoXNaorYMQsR/SRACNGNhuamsmzOKPYfrWfDF4eCgWL3oTp2H6pjQHYik8/qx8j8tA6n4jiRrn/XDQXGuIjNasZmMRFnNXf5vkKcSAKEECEwqG9yMFBs+voIu8qMrqeDzkZe/tcuMpLjuXBsDt8blontDLuMVE1vbl0Yx2aTgs1iwmqRFoY4MxIghAihQX2TGdQ3mSPHXGz6+jDb9lej61BV72HNR6W8s7WMCSOyOH9UNunJ3bNyWtX04NgFGHtzt7QupIUhOkMChBA9oF8fB4svH0Z1vYePtpXzRUkFvoCGx6fyYdFRPio6yvABaUwclcWw3NRufRPXNB2PT8XTHDCM1oXRsmj5LkR7JEAI0YPSk+OZe8FALh+fy+clFXxSXE5tow+d72Y+pSbamDAii/HDMklJ7Pzq19PxBTR8x6UIUTC6pcxmE2aTgsmkGMcmBbNZwWySANJbyUrqCBPtK41PR+rXmqbplBysYUuxkz2H61o9pwBDc1MYPzyLkflpWC3heaNWFLCYTGRnJVFb24TZpGBpDhyx1F0lK6nbkhaEEGFkMimMGpjOqIHpVNV5+GyHky92VdLkCaDz3eyneJuZMYMzOKegD/k5SZh6cNGcrhu5o5q8gTbpvxWlufVhMmE2K1hM3wUOs1np0XKK7icBQogIkZESzxXn5zPt3Dx2HKjhi52V7D5Ui66Dx6fyeUkFn5dUkOKwMXZIBmcNyaB/H0dYV1jrOgRUnYCqQjs5Ck0KrQKG2WQEje/23TBJEIlgEiCEiDAWs4mxgzMYOziDOpePr3dX8tXuY1TUGBsL1bl8fFh0lA+LjpKeFMfoQemMGZxO/8zEiHuz1XTQVA1U2g0gYLSiLCYFi9kUHP+wmFvGQmT8I5wkQAgRwVIcNqac3Z+Lx/XjaFUTX+8+xrf7qqhzGWmpqxu8fFB0lA+KjpLssDEyP42R+WkM7pccNbOTNE3Hp+mtBs5bKIBiUoyNa5SWcwqKYuyRrSjHXdMcHFvO0XwdGDsFohgtGoKvb74HSJ6rk5AAIUQUUBSFfn0c9OvjYOb5AyhzNlK0t4pt+6toaDI+mte7fHy63cmn253YLCaG5qYwPC+VYXmpIZkN1RN0QNd0tDZnu59fUaipaYLjgoURPIwAc3xAalW+E4rTOvAcH8iM75quo2o6uma8UGluNSmcKkjprX6Wruut/hVMx/0Mk8n4mS3nWp7vCgkQQkQZk6KQn5NEfk4Ssy/I51BFI8X7qyneX0118x7bvoDG9tIatpfWAJCdZqcgN5WhuSkMzEk649XbsUjXjS6xNu/4xrM9XZxul52d3OnXSIAQIoqZFIUB2UkMyE5i5nkDqKh1U9KcILCsojH4XuesceOscfPht0cxmxTyshMZ3DeZwf2SyctKCtsUWhHZJEAIESMURSE7LYHstASmnN2fJk+APYdr2VVWy+6yOhqap6iqmk7p0QZKjzaw8cvDWMwKuZmJ5OckMTDHCDb2OHlrEBIghIhZCfEWzhrSh7OG9EHXdZw1bvYermPPoTpKyxuC2WADqk5peQOl5Q1san5tZqqdAdmJDMhKJDcrkay0hPBVRIRN1AWIxAQr7gQrrWYinDArAU4YlDlxVOkU9A70NZ7p2vNTvT4zzY4SUOlIn2dwwKrVyZZvers/67vj72p6qvK0LLQ/7mXHPW79XNufpZ9wTHPuH+WEsrbWPNGk1WwVaJ6BcvzMlVaDhu39f1CCA4wc95qTlbG9f8e2Bfvuv5Mx/19D1XR8fhV/O7NwIoWiKOSkJ5CTnsCFY/uiajpHjrnYd6SO0qMNHHA2BHM1AVTWuqmsdfPFzkoArGYTA/omkZVqp3/zYHlmqh1zDK2kFm1FXaoNIOaXw0v9opOu6ySnJnC0vN6YZaIb5zTdmMrZMnslEmm6jrO6iQPlDRx0NnKwooHqeu8pX2MxK2SlJdA3PYGcjASy0xPITrOTaLdG5bTRWE8DM3Z4dqdfE3UtCCEilaIoxNssJNqtJ71G1TSaPAHcPhUtgoKFSVHom+Ggb4aD80cb5xrdfg5XNlJW0cjhShdHqlzBKbVgtKCOHHNx5FjrN9WEeAtZaXayUu1kBr/iSUmMi7iFfOLUJEAI0YPMJhNJCTYS7To+v4Y3oBIIaPgDWsRNpEy0Wxk+II3hA9IASEtLoLSshiPHXBw+5qK8uomjVU3UNLRuaTR5AsFB8ONZzAp9UuxkJMeTkRJHRnI86cnxpCfHkeyIk+6qCCQBQogwUBSFOJuZOJuxHkHXjZXEXr+Kz6cSiKDWRQtFUUhJjCMlMY6RA9OD5z2+AM5qN+XVTThrmnBWG+MXJyb2C6g65dVNlFc3tbm3SVFITbKRmhhHWlJc8HtKonEuxWGLmpXhsUQChBARQFGU4I5vJEBA1fD5NfwBFW9Ai6juqBPF2yzBhXvHa/IEgoPdx+rcVNZ6OFbnobre02YsRtN1quu9pxz3cMRbSHHYSHbYSEowvicnWI0WWct3u0XyN3UjCRBCRCCLuWWnN+NPVNOMGVMBVcOv6vgDKgE1coMGGGMR7QUOTdOpb/JxrM5DTYOX6nojaNQ2+qhp8LZpebRweQK4PAGOVLVtgRzPHmeMAyXarTjsFhzxxuOEeAuOeAsJ8VYS4iwkxBtfVrMpKgfVe4IECCGigMmkYDOZW6XI0HSdQEAzpttqGupx024jmcmkkJpodCO1xxdQqW30Udfopa7RR22jlzqXj3qXjzqXj7pGX3ANR3vc3gBur9F66QiLWTGCSoINm8WE3WbBHmcm3mYhPs5MvK35sc0c/IqzWoizmoizmbFZYnefbwkQQkQpk6Jgs5qxnTBpSm+eTquqxtTalqm2/oCGX43s7ioAm8VMVqoxC+pkfH6VhiY/dS4fjW4fDU1+Gpr8NLpbf7nc/tMGzICqB1/fVVaLiTirGZu15bsZm8XU5rvVYsJmMWNtfnz8V8se4daWvcKbH1vMprBtvhTSALF582YeeeQRNE1j0aJF/OhHP2r1vK7rPPLII2zatIn4+Hh+85vfMHr06FAWSYiYtW1/FR8WHaWy1k1mqp2LzurLmEEZ/LVwO1t3VOBXNaxmExNGZDJ+eBYffXuUY3Vu0pPiOXdkFkP7p1JysIbPSyqoafCSlhTHhBFZHKps5LPtTpp8Kgk2MxNHZXPpObknLcfuQ7Vt7gG0OVeQm9rh1xfkpvLeV4eMcngDJMRZTlmOlntU13uMQfX8NDJT7ew5XEfJgRoamvzYrCbSk+KxWky4fQH8AZ2GJh9ubwCvT+3UrDJ/80w0OtZo6RJz854ZFrMSDBoWs8nYxa/luHlnv5Zd/iwtj82mLq2DCNlCOVVVmTFjBn/729/Izs7m6quv5oknnmDo0KHBazZt2sSLL77I008/zTfffMMjjzzCypUrT3vvWF1oBbG9kAykfqGybX8Vr27a1+a8I97CjuaMri10jPGBzBM+oY8fnskXOytbrS6vbfTicvuD6a5bnrv83Dwu+17ud4sBMb7vKqvl7c/KWt3X4wugAHG21p9HZ0zMaxMkdh9q+3qAnHQ73+6tanP+ku/1bxMkTnaP0YPSKN5f0+Z8SzmOXyi3q6yG9Z+WtWqB6ZrOmMHppCbG4fGpeP3NX82P/c2z0Lw+FV9AM74iaIX92v+Z3+nXhKwFUVRURH5+Pnl5eQDMnj2bDRs2tAoQGzZsYMGCBSiKwtlnn019fT0VFRVkZWWFqlhCxKQPi462e77kQNs3RAC3J9Dm3PtfHSYpwdZqwNblDqDrxp4Fx5//dLuTRZcMbXOP1zfv+y6VSsvP8ho/y2G3tkrZUrT3GBOGG3/rLTHpmz3H2unP19sNDgBbd1Qw6/yBrfK8fLW7kvZ6Y7buqMDRziLGr3ZVMnpgutE11JzV9uvdx4xP4ifs0VDn8jHvwkHtlqU9RgpxY6zI5zcChz+g4Qt897jNl2oElZZULi1dg+pxxwFVI6DpqKpxndp8HGi+TtWaz59hd2LIAoTT6SQnJyd4nJ2dTVFR0SmvycnJwel0njZAZGYmnfL5aCf1i27hqF9No6/dlN2aTts3y+Z8Wide7/IESE+OP+H1zZvaBHdrM743eQLt1rOm0ddmrwlNAxTanG9wBxiQm9bmXLyt7V4VAU3HZml73uNTGTow44R6qO1moz1W56FPO+MaTT6V4UMyAYLPu7zb272H+7hre4OQBYj2eq5OlSjtZNcIIU7viZ9NCXcRgDMvR3fUI1LuEQtCtqIkJyeH8vLy4HF7LYMTrykvL5fuJSGEiBAhCxBjx46ltLSUsrIyfD4fhYWFTJ06tdU1U6dO5Y033kDXdb7++muSkpIkQAghRIQIWReTxWLhgQceYNmyZaiqysKFCykoKGDFihUALF68mClTprBp0yamTZuG3W7n0UcfDVVxhBBCdFJU7gchhBAi9CSrlRBCiHZJgBBCCNGuiM7F5PV6Wbp0KT6fL7gy+6c//Sm1tbX8/Oc/5/Dhw/Tv35///d//JSUlJdzF7ZKW8Zns7GyeeuqpmKrb1KlTcTgcmEwmzGYzr732WkzVr76+nvvvv59du3ahKAqPPvoogwYNion67du3j5///OfB47KyMn7605+yYMGCmKjfc889x8qVK1EUhWHDhvHYY4/hdrtjom4Azz//PCtXrkTXdRYtWsQPf/jDLv3tRXQLwmaz8fzzz7NmzRreeOMNPvjgA77++muWL1/OpEmTeOedd5g0aRLLly8Pd1G77IUXXmDIkCHB41iqGxj/UVevXs1rr70GxFb9HnnkESZPnsz69etZvXo1Q4YMiZn6DR48mNWrVwd/d3a7nWnTpsVE/ZxOJy+88AKvvvoq69atQ1VVCgsLY6JuALt27WLlypWsXLmS1atX8/7771NaWtql+kV0gFAUBYfDAUAgECAQCKAoSjBFB8CCBQt49913w1nMLisvL+f999/n6quvDp6LlbqdTKzUr7Gxka1btwZ/dzabjeTk5Jip3/E++eQT8vLy6N+/f8zUT1VVPB4PgUAAj8dDVlZWzNRt7969jBs3DrvdjsVi4dxzz+Vf//pXl+oX0QECjF/k/PnzueCCC7jgggsYN24cVVVVwfUSWVlZVFdXh7mUXfPoo4/yb//2b5iO2wErVurW4uabb+aqq67in//8JxA79SsrKyM9PZ377ruPBQsW8Mtf/pKmpqaYqd/xCgsLmTNnDhAbv7/s7GxuuukmLr30Ui666CISExO56KKLYqJuAMOGDePzzz+npqYGt9vN5s2bKS8v71L9Ij5AmM1mVq9ezaZNmygqKmLXrl3hLlK3eO+990hPT2fMmDHhLkrIrFixgtdff52nn36al19+ma1bt4a7SN0mEAiwfft2Fi9ezBtvvIHdbo/aLolT8fl8bNy4kZkzZ4a7KN2mrq6ODRs2sGHDBj744APcbjerV68Od7G6zZAhQ1i2bBk33XQTy5YtY/jw4ZjNbfNYdUTEB4gWycnJnHfeeXzwwQdkZGRQUVEBQEVFBenp6ad5deT58ssv2bhxI1OnTuWuu+5iy5Yt3H333TFRtxbZ2Ub++YyMDKZNm0ZRUVHM1C8nJ4ecnBzGjRsHwMyZM9m+fXvM1K/F5s2bGT16NH369AGIifp9/PHH5Obmkp6ejtVqZfr06Xz11VcxUbcWixYt4vXXX+fll18mNTWV/Pz8LtUvogNEdXU19fX1AHg8Hj7++GMGDx4cTNEB8MYbb3DZZZeFs5hd8otf/ILNmzezceNGnnjiCc4//3x++9vfxkTdAJqammhsbAw+/uijjygoKIiZ+mVmZpKTk8O+fcYeDJ988glDhgyJmfq1KCwsZPbs2cHjWKhfv379+Oabb3C73ei6HpO/u6oqIz36kSNHeOedd5gzZ06X6hfRK6lLSkq49957UVUVXdeZOXMmd9xxBzU1NfzsZz/j6NGj9O3bl9///vekpra/O1U0+PTTT3n22Wd56qmnYqZuZWVl/OQnPwGMcaQ5c+Zw++23x0z9AHbs2MEvf/lL/H4/eXl5PPbYY2iaFjP1c7vdXHLJJbz77rskJRmpvWPl9/fkk0/y5ptvYrFYGDlyJI888ggulysm6gawZMkSamtrsVgs3HfffUyaNKlLv7uIDhBCCCHCJ6K7mIQQQoSPBAghhBDtkgAhhBCiXRIghBBCtEsChBBCiHZFdDZXIU5l0aJF+Hw+/H4/paWlFBQUADBq1Cgee+yxMJeuY4qLiykrK4uplcoidsg0VxH1Dh06xMKFC/n000/DXZQ2AoEAFsvJP4etXLmSjz/+mN/97nfdfm8hzpT87xIxadWqVfzjH/9AVVWSk5N56KGHGDhwICtXrmT9+vU4HA527dpF3759+Y//+A8ef/xxysrKGDduHI8//jiKonD33Xdjt9s5ePAg5eXlnHfeefzqV7/CarXS0NDAo48+yu7du/F6vVxwwQXcc889mEwmFi9ezMSJE/nqq69ISEjgySefDC4S9Hq9jBs3joceeoj6+nr+9Kc/4XK5mD9/Pueddx5Lly5lyZIlfPTRRwAcOHAgeHzgwAEWL17MNddcw5YtW7jqqquYP38+TzzxBJ9//jk+n4+RI0fy4IMPYrfbw/wbEDFBFyLKlZWV6RMnTgweb9myRb/11lt1r9er67qub9iwQV+6dKmu67r+yiuv6BMnTtTLy8t1Xdf1m266SV+wYIHe0NCg+3w+fdasWfqWLVt0Xdf1X/ziF/r8+fN1l8ul+3w+/frrr9f//ve/67qu6/fcc4++du1aXdd1XVVV/ac//am+atUqXdd1/dprr9V//OMf64FAIPh8bW1t8PFdd92lv/LKK8Hy/OxnPwuWvbS0VL/gggvaPS4tLdWHDRumr1+/Pvj8k08+qT/11FPB48cee0z//e9/f2b/oEI0kxaEiDkbN25k+/btLFq0CABd13G5XMHnx48fH0wkOGrUKDweD4mJiQAMHz6cgwcPct555wEwa9YsEhISACOH/vvvv8/ixYt57733KC4u5umnnwaMXGEDBgwI/oy5c+cGM2hqmsby5cv58MMP0TSN2traLu9UlpCQwIwZM1rV1e12U1hYCBjZV0ePHt2lewtxIgkQIubous73v/997rjjjnafj4uLCz42mUxtjgOBwEnvqygKYLzpP/XUU/Tr16/da1uCCsDq1aspKiri73//Ow6Hgz/+8Y8cPXq03deZzWY0TQsee73ek963pUy//vWvOffcc9u9nxBnQqa5ipjTkrXS6XQCRrLAbdu2deleb731Fm63G7/fz9q1a4Mti6lTp7J8+XJUVQWMzMNlZWXt3qOhoYG0tDQcDgd1dXXBT/sADoeDhoaG4HFWVhYejyd4r3Xr1p22rs8++2wwkDQ2NrJ3794u1VWIE0mAEDHn/PPP54477uDWW29l3rx5zJ07l/fff79L9xo/fjy33347c+bMIS8vL7jF6K9+9Ss0TWP+/PnMnTuXW265hcrKynbvceWVV1JbW8ucOXO46667Wn3av/DCC2loaGDevHk8+uij2Gw27r33Xm644Qauu+46rFbrKct32223MWTIEK6++mrmzp3L0qVL2b9/f5fqKsSJZJqrECdx9913M378eBYvXhzuoggRFtKCEEII0S5pQQghhGiXtCCEEEK0SwKEEEKIdkmAEEII0S4JEEIIIdolAUIIIUS7/n9/w7e4SFWOPQAAAABJRU5ErkJggg==\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
+}
--
2.18.1