{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 3. Inference\n", "\n", "Read in the packages, data and exclude an unused variable." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy as sp\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.2 Testing examples\n", "\n", "Read in the Galapagos data" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SpeciesAreaElevationNearestScruzAdjacent
Baltra5825.093460.60.61.84
Bartolome311.241090.626.3572.33
Caldwell30.211142.858.70.78
Champion250.10461.947.40.18
Coamano20.05771.91.9903.82
\n", "
" ], "text/plain": [ " Species Area Elevation Nearest Scruz Adjacent\n", "Baltra 58 25.09 346 0.6 0.6 1.84\n", "Bartolome 31 1.24 109 0.6 26.3 572.33\n", "Caldwell 3 0.21 114 2.8 58.7 0.78\n", "Champion 25 0.10 46 1.9 47.4 0.18\n", "Coamano 2 0.05 77 1.9 1.9 903.82" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gala = pd.read_csv(\"data/gala.csv\", index_col=0)\n", "gala.drop('Endemics', axis=1, inplace=True)\n", "gala.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the model:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: Species R-squared: 0.766
Model: OLS Adj. R-squared: 0.717
Method: Least Squares F-statistic: 15.70
Date: Fri, 07 Sep 2018 Prob (F-statistic): 6.84e-07
Time: 15:12:52 Log-Likelihood: -162.54
No. Observations: 30 AIC: 337.1
Df Residuals: 24 BIC: 345.5
Df Model: 5
Covariance Type: nonrobust
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 7.0682 19.154 0.369 0.715 -32.464 46.601
Area -0.0239 0.022 -1.068 0.296 -0.070 0.022
Elevation 0.3195 0.054 5.953 0.000 0.209 0.430
Nearest 0.0091 1.054 0.009 0.993 -2.166 2.185
Scruz -0.2405 0.215 -1.117 0.275 -0.685 0.204
Adjacent -0.0748 0.018 -4.226 0.000 -0.111 -0.038
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 12.683 Durbin-Watson: 2.476
Prob(Omnibus): 0.002 Jarque-Bera (JB): 13.498
Skew: 1.136 Prob(JB): 0.00117
Kurtosis: 5.374 Cond. No. 1.90e+03


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.9e+03. This might indicate that there are
strong multicollinearity or other numerical problems." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: Species R-squared: 0.766\n", "Model: OLS Adj. R-squared: 0.717\n", "Method: Least Squares F-statistic: 15.70\n", "Date: Fri, 07 Sep 2018 Prob (F-statistic): 6.84e-07\n", "Time: 15:12:52 Log-Likelihood: -162.54\n", "No. Observations: 30 AIC: 337.1\n", "Df Residuals: 24 BIC: 345.5\n", "Df Model: 5 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 7.0682 19.154 0.369 0.715 -32.464 46.601\n", "Area -0.0239 0.022 -1.068 0.296 -0.070 0.022\n", "Elevation 0.3195 0.054 5.953 0.000 0.209 0.430\n", "Nearest 0.0091 1.054 0.009 0.993 -2.166 2.185\n", "Scruz -0.2405 0.215 -1.117 0.275 -0.685 0.204\n", "Adjacent -0.0748 0.018 -4.226 0.000 -0.111 -0.038\n", "==============================================================================\n", "Omnibus: 12.683 Durbin-Watson: 2.476\n", "Prob(Omnibus): 0.002 Jarque-Bera (JB): 13.498\n", "Skew: 1.136 Prob(JB): 0.00117\n", "Kurtosis: 5.374 Cond. No. 1.90e+03\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "[2] The condition number is large, 1.9e+03. This might indicate that there are\n", "strong multicollinearity or other numerical problems.\n", "\"\"\"" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols(formula='Species ~ Area + Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", "lmod.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Total sum of squares and residual sum of squares" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(381081.36666666664, 89231.36633005117)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.centered_tss, lmod.ssr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Degrees of freedom for the numerator and denominator of the F-test" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(5.0, 24.0)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.df_model, lmod.df_resid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mean square of numerator and denominator" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(58370.0000673231, 3717.9735970854654)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.mse_model, lmod.mse_resid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "F-statistic reproduced" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "15.699412204831035" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.mse_model/ lmod.mse_resid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute p-value" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.837892995159578e-07" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1-sp.stats.f.cdf(lmod.fvalue, lmod.df_model, lmod.df_resid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "F-test for comparing two models. (Turn off the warnings about NaNs in the output - we don't care)." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
df_residssrdf_diffss_diffFPr(>F)
025.093469.0839900.0NaNNaNNaN
124.0254156.6807241.0-160687.596734-15.173721.0
\n", "
" ], "text/plain": [ " df_resid ssr df_diff ss_diff F Pr(>F)\n", "0 25.0 93469.083990 0.0 NaN NaN NaN\n", "1 24.0 254156.680724 1.0 -160687.596734 -15.17372 1.0" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%capture --no-display\n", "lmods = smf.ols(formula='Species ~ Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", "sm.stats.anova_lm(lmods,lmod)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
df_residssrdf_diffss_diffFPr(>F)
026.0158291.6285680.0NaNNaNNaN
124.089231.3663302.069060.2622389.2873520.00103
\n", "
" ], "text/plain": [ " df_resid ssr df_diff ss_diff F Pr(>F)\n", "0 26.0 158291.628568 0.0 NaN NaN NaN\n", "1 24.0 89231.366330 2.0 69060.262238 9.287352 0.00103" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%capture --no-display\n", "lmods = smf.ols(formula='Species ~ Elevation + Nearest + Scruz', data=gala).fit()\n", "sm.stats.anova_lm(lmods,lmod)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Use I() notation to test replacing by the sum of two variables" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
df_residssrdf_diffss_diffFPr(>F)
025.0109591.1208010.0NaNNaNNaN
124.089231.3663301.020359.7544715.4760350.027926
\n", "
" ], "text/plain": [ " df_resid ssr df_diff ss_diff F Pr(>F)\n", "0 25.0 109591.120801 0.0 NaN NaN NaN\n", "1 24.0 89231.366330 1.0 20359.754471 5.476035 0.027926" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%capture --no-display\n", "lmods = smf.ols(formula='Species ~ I(Area+Adjacent) + Elevation + Nearest + Scruz', data=gala).fit()\n", "sm.stats.anova_lm(lmods,lmod)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Need to use glm function to get offset functionality" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "11.318196837955563" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.glm(formula='Species ~ Area + Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", "lmods = smf.glm(formula='Species ~ Area + Nearest + Scruz + Adjacent', offset=(0.5*gala['Elevation']), data=gala).fit()\n", "fstat = (lmods.deviance-lmod.deviance)/(lmod.deviance/lmod.df_resid)\n", "fstat" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.002573836486092773" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1-sp.stats.f.cdf(fstat, 1, lmod.df_resid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute t-statistic and corresponding p-value" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-3.3642527904358723, 0.0025738364860927276)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols(formula='Species ~ Area + Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", "tstat=(lmod.params['Elevation']-0.5)/lmod.bse['Elevation']\n", "tstat, 2*sp.stats.t.cdf(tstat, lmod.df_resid)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "11.318196837955554" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tstat**2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.3 Permutation tests\n", "\n", "Permutation of the response" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 8, 104, 25, 44, 51, 16, 280, 31, 285, 70, 58, 58, 18,\n", " 5, 2, 2, 108, 24, 3, 347, 12, 40, 97, 93, 2, 10,\n", " 62, 21, 444, 237])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols(formula='Species ~ Nearest + Scruz', data=gala).fit()\n", "np.random.permutation(gala.Species)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the f-statistic for a sample of permutations" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "fstats = np.zeros(4000)\n", "for i in range(0,4000):\n", " gala['ysamp'] = np.random.permutation(gala.Species)\n", " lmodi = smf.ols(formula='ysamp ~ Nearest + Scruz', data=gala).fit()\n", " fstats[i] = lmodi.fvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Proportion of permuted f-statistics that exceed the observed value" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.568" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(fstats > lmod.fvalue)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Which is close to the observed p-value" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5549254563908441" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.f_pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do t-test with permutation" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "tstats = np.zeros(4000)\n", "for i in range(0, 4000):\n", " gala['ssamp'] = np.random.permutation(gala.Scruz)\n", " lmodi = smf.ols(formula='Species ~ Nearest + ssamp', data=gala).fit()\n", " tstats[i] = lmodi.tvalues[2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Proportion of permuted t-statistcs which exceed the observed t-statistic in absolute value" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.28025" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.mean(np.fabs(tstats) > np.fabs(lmod.tvalues[2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Very close to observed p-value:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.2833295186486556" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.pvalues[2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.5 Confidence intervals for beta" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.07713396, 0.07945568])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols(formula='Species ~ Area + Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", "qt = np.array(sp.stats.t.interval(0.95,24))\n", "lmod.params[1] + lmod.bse[1]*qt" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01
Intercept-22.800754127.659321
Area-0.0771340.079456
Elevation-0.2241030.156443
Nearest1.0878956.825321
Scruz-0.7863850.434932
Adjacent-0.0052520.121375
\n", "
" ], "text/plain": [ " 0 1\n", "Intercept -22.800754 127.659321\n", "Area -0.077134 0.079456\n", "Elevation -0.224103 0.156443\n", "Nearest 1.087895 6.825321\n", "Scruz -0.786385 0.434932\n", "Adjacent -0.005252 0.121375" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.conf_int()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Drawing the confidence ellipses appears to require some effort in coding" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.6 Bootstrapped Confidence Intervals\n", "\n", "Bootstrap the residuals" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "breps = 4000\n", "coefmat = np.empty((breps,6))\n", "resids = lmod.resid\n", "preds = lmod.predict()\n", "for i in range(0,breps):\n", " gala['ysamp'] = preds + np.random.choice(resids,30)\n", " lmodi = smf.ols(formula='ysamp ~ Area + Elevation + Nearest + Scruz + Adjacent', data=gala).fit()\n", " coefmat[i,:] = lmodi.params" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Turn into pandas dataframe and add columns" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "coefmat = pd.DataFrame(coefmat, columns=(\"intercept\", \"area\",\"elevation\",\"nearest\",\"Scruz\",\"adjacent\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the quantiles by column for bootstrap confidence intervals" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptareaelevationnearestScruzadjacent
0.025-7.527822-0.064221-0.1922201.710021-0.6338710.005495
0.975120.5659760.0882430.1405166.6839280.4502770.124617
\n", "
" ], "text/plain": [ " intercept area elevation nearest Scruz adjacent\n", "0.025 -7.527822 -0.064221 -0.192220 1.710021 -0.633871 0.005495\n", "0.975 120.565976 0.088243 0.140516 6.683928 0.450277 0.124617" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coefmat.quantile((0.025,0.975))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the kernel density estimate" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl4nXWZ8PHvnX1fm6RJl4QudKMsJSCLigpoUcER9cUFB9QZZMZxFHVU1BlmfOd1m/dS0ddhrIriMqgso6LICHWQRbYUSmkbaNrSLE3aJM2+L+d+/zjnQBrS5uTkPNs59+e6ejV5zpPndz+/puc+z28VVcUYY0zqSvM6AGOMMd6yRGCMMSnOEoExxqQ4SwTGGJPiLBEYY0yKs0RgjDEpzhKBMcakOEsExhiT4iwRGGNMisvwOoBYLFmyROvq6rwOw/jQwa5hAFZV5HscSXKw+kwuO3bs6FbVivnOC0QiqKuro6GhweswjA9d9d3HAPjFh8/3OJLkYPWZXESkOZbzrGnIGGNSXCCeCIw5kY++Ya3XISQVq8/UZInABNqr1y7xOoSkYvWZmqxpyATanvZ+9rT3ex1G0rD6TE2WCEygffGevXzxnr1eh5E0rD5TkyUCY4xJcY4lAhG5VUQ6RWT3HK99SkRURKxB0hhjPObkE8GPgK2zD4rICuBSoMXBso3xxPbGo3z/4YMc6R/zOhRjYuZYIlDVh4CeOV76BvBpwDZLNknl5gea+NBtDfzr7xp567cf5sXuYa9DMiYmrg4fFZErgMOq+qyIuFm0SVKf3rrO6xAA2Hd0kJu37+NtZ9bw169ZxdU/eIJP/HInd//NBQTpd90v9Wnc5VoiEJE84PPAG2M8/zrgOoCVK1c6GJkJsrNry7wOAYDvPXSQrIw0/vnyTZTmZ3HjZev5zF3P8eALXbx+faXX4cXML/Vp3OXmqKHVwCnAsyJyCFgOPC0iS+c6WVW3qWq9qtZXVMy7ZpJJUTuae9jRPFcLpHvGp6a597kO3nbGMkrzswC4cstylpXkcsuDBzyNbaH8UJ/Gfa4lAlV9TlUrVbVOVeuANmCLqh5xKwaTfL523wt87b4XPI3hzweOMTwxzdbNL3+myUxP4+rzannyUA8tx0Y8jG5h/FCfxn1ODh+9HXgMWCcibSLyIafKMsZLf9hzlPysdC5YXX7c8SvOrAHgnl3tXoRlTMycHDX0HlWtVtVMVV2uqj+Y9XqdqnY7Vb4xbnn84DHOX72E7Iz0444vK8mlvraU3+3q8CgyY2JjM4uNWYRjQ+O82D1MfV3pnK+/fn0lezsG6BocdzkyY2JnicCYRXi6pQ+As2vnTgSvXRse6PDI/i7XYjJmoWwZahNo/3T5Rk/L39HcS2a6sHlZ8Zyvb6opoiw/i4f3dfP2s5a7HN3CeV2fxhuWCEygbaqZ+w3YLXva+1m3tJCczPQ5X09LE151ShlPBWRIptf1abxhTUMm0B5p6uaRJu/GHLxwZJB1VUUnPefs2lJae0bpHPD/+kNe16fxhj0RmED79h+bAG921uodnqBzcJz1SwtPet6WSP/B0y29bD2t2o3Q4uZlfRrv2BOBMXF64eggAKfOkwg21RSRlZHGjuZeN8IyZsEsERgTp32RRDDfE0F2RjobqovYfXjAjbCMWTBLBMbE6fkjgxTnZlJZmD3vuRurC9nbMYCqrb5u/McSgTFx2ndkkHVVhTEtM72xuoj+0UnabcMa40PWWWwC7UtXbvas7IPdw7xp05yL577CxprwyKLG9gGWleQ6GdaieFmfxjuWCEygra4o8KTcofEpeoYnWFmWF9P565YWIQJ7Owa4ZGOVw9HFz6v6NN6ypiETaA/sPcoDe4+6Xm5rT3hp6RVlsX26L8jOoLYsj73t/u4w9qo+jbfsicAE2vcePgjg+qfsaCKI9YkAYEN1EY0d/k4EXtWn8ZY9ERgTh9beUQBWlMaeCNZUFtDaO8r41LRTYRkTF0sExsShtWeEguwMSvIyY/6Z1RUFTIc0UDuWmdRgicCYOLT2jLC8NDemoaNRqyryATjQNexUWMbExRKBMXFo7R1ZUP8AwKrIiJwDXUNOhGRM3Kyz2ATaN6460/UyVZXWnlFeE9l0JlYF2RlUFWVz0MdPBF7Up/Gek5vX3yoinSKye8axfxOR50Vkl4j8l4iUOFW+SQ01JbnUuDxBq3togtHJaVaULrzc1RUFvn4i8KI+jfecbBr6EbB11rH7gdNU9XRgH3Cjg+WbFHDPs+3c82y7q2W29kaGjpYvrGkIXk4Efl1zyIv6NN5zLBGo6kNAz6xjf1DVqci3jwP+37vP+NpPH2/mp483u1rmS5PJFjB0NOqUJfkMjk1xbHgi0WElhBf1abznZWfxB4Hfe1i+MXGJJoLlcSSC2shTRPQaxviBJ4lARD4PTAE/O8k514lIg4g0dHV1uRecMfNo7RllSUE2uVlz71N8MtGRRi2WCIyPuJ4IROQa4K3A+/QkDaWquk1V61W1vqJiYaMzjHFSeOhofB2q0acIeyIwfuJqIhCRrcBngCtU1f4nmEBq6RlhxQLnEETlZqVTWZhtTwTGVxybRyAitwOvA5aISBtwE+FRQtnA/ZEZmY+r6vVOxWCS3y1Xn+1qeVPTITr6x+LqKI5aWZbn20Tgdn0af3AsEajqe+Y4/AOnyjOpqSw/y9XyOvrHmA7pgmcVz7SyLI8nXuyZ/0QPuF2fxh9siQkTaHc0tHJHQ6tr5UU/yS+Ps48AYEVZHu39o0xMhRIVVsK4XZ/GHywRmEC7c0cbd+5oc628xcwhiFpZlocqHO4bTVRYCeN2fRp/sERgzAK09o6QniZUF+fEfY0VNoTU+IwlAmMWoKVnlGUluWSkx/9fp6YknEQ6fPhEYFKTJQJjFqC1ZyTmfYpPpKoohzSBdksExicsERizAG29I4vqHwDITE+jsjCH9v6xBEVlzOLYfgQm0H70gXNdK2t4fIruoYm4J5PNVFOS48snAjfr0/iHJQITaPGs9xOvtuiG9QlIBNUluextH1j0dRLNzfo0/mFNQybQfvLYIX7y2CFXynp56OjiN25ZVpJLe9+o7/YlcLM+jX9YIjCB9ttdHfx2V4crZUWHey5mVnFUdXEO41Mheny2L4Gb9Wn8wxKBMTFq7R0hLys9IcswRLeDbO+zDmPjPUsExsSotWeUFaV5RBZMXJSa4nAi8OPsYpN6LBEYE6NEzCGIemlSWb8lAuM9SwTGxEBVae2Nfx+C2crys8jOSPPlEFKTemz4qAm0X3z4fFfK6RmeYGRietGTyaJEhJqSXN9NKnOrPo2/2BOBMTFoTeAcgqjqYn9OKjOpxxKBCbRtDx1g20MHHC8nkUNHo5YW5dA5MJ6w6yWCW/Vp/MUSgQm07Y2dbG/sdLyc6GSy5QmYTBZVVZxD5+AYoZB/JpW5VZ/GXywRGBODtt4RyvOzyM9OXLdaVWE2k9NK74i/JpWZ1ONYIhCRW0WkU0R2zzhWJiL3i0hT5O9Sp8o3JpFaekZYnsBmIQgvRw1wZMBfHcYm9Tj5RPAjYOusY58FtqvqWmB75HtjfK+lZ4TaRCeCyC5nfusnMKnHsUSgqg8BPbMOvw24LfL1bcBfOFW+SQ05menkZDq7YubkdIj2vjFqy515IjjqoycCN+rT+I/b8wiqVLUDQFU7RKTS5fJNkrntg86vn3+4d5TpkCZ06ChARUE24K+mITfq0/iPbzuLReQ6EWkQkYauri6vwzEprDkyYijRTUNZGWmU52dx1JqGjMfcTgRHRaQaIPL3Ccepqeo2Va1X1fqKigrXAjTB8q3tTXxre5OjZUTnENSW5yf82lVFOXT66InAjfo0/uN2IvgNcE3k62uAX7tcvkkyj+7v5tH93Y6W0XJsmOyMNCoLsxN+7aqibI4O+icRuFGfxn+cHD56O/AYsE5E2kTkQ8BXgEtFpAm4NPK9Mb7WfCy82Fxa2uKXn56tqiiHI/3WNGS85Vhnsaq+5wQvXexUmcY4wYmho1FVRTkcGx5ncjpEZrpvu+xMkrPfPGNOQlVp6RlhZYKHjkZVFeWgCt1D9lRgvGPLUJtAK81b/LaRJ9M9FF5+OpGLzc1UVRTudzg6ME51ceLWMYqX0/Vp/MkSgQm0/3j/2Y5e/+URQ849EQAc6R+DFY4UsSBO16fxJ2saMuYkWnqGAVhZlviho/ByIuj00cghk3osEZhA++p9z/PV+5537PrNx0YQSezy0zOV52eRnia+WWbC6fo0/mRNQybQnm7udfT6h7qHqSnOdWz9nbQ0obIw2zdDSJ2uT+NP9kRgzEkc6BpmVYUzzUJRlUU51jRkPGWJwJgTUFUOdA2xuqLA0XKWFmX7pmnIpCZLBMacwJGBMUYmplld6WwiCM8utkRgvGN9BCbQqiObuzjhYFd4xNBqh5uGqopyGBibYmxy2vO9AJysT+NflghMoH3z3Wc5du0DXUMAjjcNRRez6xwYd2wGc6ycrE/jX9Y0ZMwJHOgcoiA7w5FVR2d6aacy6zA2HrFEYALtX+7Zw7/cs8eRax/oGmZ1RT4iiV91dCY/bVnpZH0a/7KmIRNoe9sHHLv2ga4hzl9V7tj1o2auN+Q1J+vT+Jc9ERgzh+HxKTr6xxyfQwBQnJtJVkaar3YqM6nFEoExc3jh6CAAp1YVOl6WiIR3KrNEYDxiicCYOTR2hJtINlQXuVJeVWGOL5qGTGqyPgITaE413TR2DFCYneHYYnOzVRXl0HjE+/Z5N5rCjP9YIjCB9uUrT3fkus93DLK+utDxEUNRlUXZ/Gmf908ETtWn8beYmoZE5C4ReYuIJKQpSURuEJE9IrJbRG4XEZvOaHwjFFKePzLoWrMQhJ8IhsanGBqfcq1MY6JifWO/BXgv0CQiXxGR9fEWKCLLgL8H6lX1NCAdeHe81zOp7ca7d3Hj3bsSes223lGGxqdYv9S9RPDy7GJvO4ydqE/jfzE1DanqA8ADIlIMvAe4X0Rage8BP1XVyTjKzRWRSSAPaF/gzxsDvLweUCLtfamj2PkRQ1EvTyobZ5XDS1qcjBP1afwv5qYeESkHrgX+CngGuBnYAty/kAJV9TDwf4EWoAPoV9U/LOQaxjjp+SMDiMC6pW4mgsgTgS0zYTwQax/B3cDDhD+9X66qV6jqL1T1o8CCPr6ISCnwNuAUoAbIF5Gr5zjvOhFpEJGGrq6uhRRhzKLsPtzPKUvyyctybyxFZXTvYhtCajwQ6xPB91V1o6p+WVU7AEQkG0BV6xdY5iXAi6raFWlSuhu4YPZJqrpNVetVtb6iomKBRRgTH1VlZ2sfZ64ocbXcwuwMcjPTbVKZ8USsH3n+Fbh31rHHCDcNLVQLcJ6I5AGjwMVAQxzXMYaNNYnt0G3rHaV7aIKzVpYm9LrzeWl28aC3TwSJrk8TDCdNBCKyFFhGuGP3LCA6qLqIcDPRgqnqEyJyJ/A0MEW4v2FbPNcy5qbLNyX0es+09gFwlstPBBBuHvL6iSDR9WmCYb4ngjcR7iBeDnx9xvFB4HPxFqqqNwE3xfvzxjhlZ0sf2RlprnYUR1UV5fBcW5/r5Rpz0kSgqrcBt4nIO1T1LpdiMiZmH//5M0DidtZ6prWX05cXk5nu/jJcVYXZPDAwjqq6NqN5tkTXpwmG+ZqGrlbVnwJ1IvKJ2a+r6tfn+DFjXNORwE3fJ6ZC7Gkf4JrzaxN2zYWoKsphdHKawfEpinIyPYkhkfVpgmO+pqHoClTezXAxxiWNHQNMTIU4c4W7HcVRlUUvzy72KhGY1DRf09B3I3//izvhGOOdZ1p6AThrpfsdxXD87OI1le73UZjUFeuEsq+JSJGIZIrIdhHpnmsSmDFBtrO1j8rCbKqLvVkD0U97F5vUEmuP2BtVdQB4K9AGnAr8g2NRGROjLbWlbKlNTFNOdCKZVx210YXnvNygJpH1aYIj1gll0QbLNwO3q2qPV/9ZjJnpM1vjXgj3OL3DExw6NsL/OmdFQq4Xj/zsDAqzMzx9IkhUfZpgiTUR3CMizxOeCfy3IlIB2POrSRo7I+P33V5aYrbKomxbeM64LqamIVX9LHA+4T0EJoFhwgvHGeOp63+yg+t/smPR19nZ0ocInL7c20RQVeTt3sWJqk8TLAtZXnED4fkEM3/mxwmOx5gF6R2ZSMh1nm3r49TKQgqyvd29taooh6cO9XhWfqLq0wRLTL/1IvITYDWwE5iOHFYsEZgkoKo829rHpRurvA6FysJsOj2eXWxST6wff+qBjaqqTgZjjBdaekboHZn0bCLZTJVFOUxMh+gbmaQ0P8vrcEyKiHX46G5gqZOBGOOVna3+6CiGl3cqO2odxsZFsT4RLAH2isiTwEs9Wap6hSNRGROjC9csWfQ1drb2kZuZzqlV3q+kMnN28XoPPnoloj5N8MSaCP7ZySCMidffX7x20dfY2drH5mXFZHiw4uhsVYXezi5ORH2a4Il1+OifgENAZuTrpwhvLGNMoE2HlMaOATYvL/Y6FODlhee6PN6pzKSWWNca+mvgTuC7kUPLgF85FZQxsbrm1ie55tYn4/75F7uHGZsMsbHaH1s05mSmU5ybSUf/qCflL7Y+TTDF+iz8EeBCYABAVZuASqeCMiZWY5PTjE1Oz3/iCeztGABgg08SAUBNSS4dfd40DS22Pk0wxZoIxlX1pZkmkUllNpTUBF5jxwCZ6cKaSu87iqOWleRwuM+bJwKTmmJNBH8Skc8R3sT+UuAO4J54CxWREhG5U0SeF5FGETk/3msZsxh72wdYU1lIVob3HcVRNSW5tFsiMC6K9bf/s0AX8BzwYeBe4AuLKPdm4D5VXQ+cATQu4lrGxK2xY8A3/QNRNSW5DIxNMTg26XUoJkXENHxUVUMi8ivgV6ratZgCRaQIeC1wbeTaE4AtcGLicvGG+LuquofG6RwcZ0O1v3YDqynJBcL7Bxe6vGXlYurTBNd8m9cLcBPwd4BEDk0D31bVL8ZZ5irCTxc/FJEzgB3Ax1R1OM7rmRR23WtXx/2zjZGO4o01/noiWFYSnktwuG+UU6vcTVKLqU8TXPM1DX2c8Gihc1S1XFXLgFcBF4rIDXGWmQFsAW5R1bMIL2n92dknich1ItIgIg1dXYt6CDFmTi8cGQRg/VJ/JYLoE4H1Exi3zJcI/hJ4j6q+GD2gqgeBqyOvxaMNaFPVJyLf30k4MRxHVbepar2q1ldUVMRZlEl2V333Ma767mNx/ez+ziHK87Mo89nibpWFOaSniSeJYDH1aYJrvkSQqardsw9G+gniarxU1SNAq4isixy6GNgbz7WMWYz9nUOs9tGw0aj0NGFpUQ7tHs0lMKlnvs7ik3XiLqaD96PAz0QkCzgIfGAR1zJmwVSV/V1DXHZatdehzGlZSa7NJTCumS8RnCEiA3McFyAn3kJVdSfhPQ6M8UTP8AR9I5O+mkg2U01JDg3NvV6HYVLESROBqqa7FYgxbtrfOQTg40SQy5FdHUyHlPQ026nMOMvbDVqNWaS3nh5f087+Lv8ngqmQ0jU4ztLiuB++Fyze+jTBZonABNr7z6+L6+f2dw6Rl5VOdZF7b7ILsSw6hLR/1NVEEG99mmDzzwIrxsRhdGKa0YmFr5a5v3OIVRX5pPm02cWruQTx1qcJNksEJtCu/eGTXPvDha+ff7BrmDUV/mwWAlhWGk4ErT3uJoJ469MEmyUCk3KGx6c43Dfq2/4BgILsDMrzs2jpsZVXjPMsEZiUc7Ar/Obq50QAsLI8j0PdI16HYVKAJQKTcvZ3hdcY8nsiqC3Lo6XHEoFxniUCk3L2dw6RniasLMv3OpSTWlmeT3v/KONT1nlrnGXDR02gvfPs5Qv+mQOdw9SW5/lqV7K51JXnoQptvaOsdqljO576NMFnicAE2rvqVyz4Z/Z3Dfl6xFBUbXkeAC3HRlxLBPHUpwk+f38kMmYePcMT9AzHvv7h5HSIQ93Dvu8fAF5qujp0zL2RQwutT5Mc7InABNrf/HQHAL/48Pkxnd98bISpkAYiESwpyCIvK53mY+51GC+0Pk1ysCcCk1Kii8251dSyGCLCShs5ZFxgicCklAORxeb8uCHNXOrK82l2sWnIpCZLBCal7O8coro4h4LsYLSK1i4JPxFMTYe8DsUkMUsEJqXs7xwKRP9A1NrKQianlWZrHjIOCsbHImNO4OrzamM+NxRSDnQNcdU5wRkiuTaStJqODrnSr7GQ+jTJwxKBCbTLz6iJ+dz2/lFGJqYD9UQQ7cvY3zkILHW8vIXUp0kenjUNiUi6iDwjIr/1KgYTfO19ozGv2R8dMbS2stDJkBKqIDuDZSW5NEVid9pC6tMkDy/7CD4GNHpYvkkCN/xiJzf8YmdM5/p9n+ITWVNZQNNRdxLBQurTJA9PEoGILAfeAnzfi/JNatrfOURZfhZl+Vleh7IgaysLONA1xHRIvQ7FJCmvngi+CXwasDFxxjVBGzEUdWpVIeNTIdp6beSQcYbriUBE3gp0quqOec67TkQaRKShq6vLpehMslJVmgKaCNZUvTxyyBgnePFEcCFwhYgcAn4OvEFEfjr7JFXdpqr1qlpfUVHhdowmyXQPTdA/OhmIVUdniyavF44OehyJSVauDx9V1RuBGwFE5HXAp1T1arfjMMnhr1+zKqbzgtpRDFCUk8nKsjz2tg84Xlas9WmSi80jMIF2ycaqmM5r7Ai/ia6vDs7Q0ZlOW1bE7vZ+x8uJtT5NcvF0iQlVfVBV3+plDCbYDnQNvbSQ3Mns7RhgSUE2lYU5LkSVeJtqimk+NkL/6KSj5cRanya52FpDJtA+d/dzfO7u5+Y9r7FjgI01RS5E5IzTlhUDsMfhp4JY69MkF0sEJulNTIVoOjrEhoA2CwFsiiSxPYed7ycwqccSgUl6B7qGmJgOsbE6uE8ESwqyqS7OcaWfwKQeSwQm6UU7ijcFuGkIwv0Euw9bIjCJZ4nAJL1dbf3kZqZTV57vdSiLcvryYg52DzveYWxSjw0fNYH20Tesnfecp1t6OWNFMRnpwf7cc3ZtKarwTEsvr1tX6UgZsdSnST6WCEygvXrtkpO+Pjoxzd72AT58UfAnSp2xooQ0gR3NziWC+erTJKdgf0QyKW9Pe/9Jh1TuautjKqScXVvqYlTOKMjOYEN1ETuaex0rY776NMnJEoEJtC/es5cv3rP3hK/vaAm/aZ61IviJAKC+tpSdrX2ObWY/X32a5GSJwCS1hkO9rKrIpzRgexCcyNl1ZYxMTNPYYQvQmcSxRGCS1vjUNI8fPMar1yRPu3d9pInrqUM9HkdikoklApO0dhzqZWRimotOTZ5lzGtKcllRlstjB495HYpJIpYITNJ6cF8XWelpnLeq3OtQEurVa5bw+IFjjvUTmNRjw0dNoH1667o5j6sq9z7Xwfmry8nPTq5f8wvXLOH2J1vZdbifLSsT2wl+ovo0yS25/oeYlHN2bdmcx59u6aOtd5QbLjnV5Yicd8HqcJ/Ho03dCU8EJ6pPk9ysacgE2o7mHnY0v7Lj9K6n28jKSOONm5Jvo5Wy/Cw21RTxyP7uhF/7RPVpkpslAhNoX7vvBb523wvHHesbmeDup9t4+5nLKMzJ9CgyZ7167RKebullaHwqodedqz5N8rNEYJLO7U+2MjYZ4gOvrvM6FMdcvL6KyWnlwRc6vQ7FJAFLBCapDI1P8f2HD/KatUtYvzTYy06fzNm1pZTlZ/GHPUe9DsUkAdcTgYisEJH/EZFGEdkjIh9zOwaTvLb96QDHhif41BuTe/RLeppwyYZK/uf5TiambBipWRwvngimgE+q6gbgPOAjIrLRgzhMkjk6MMb3Hn6Ry8+o4YwVJV6H47g3bVrK4PgUj9vkMrNIrg8fVdUOoCPy9aCINALLAFvpyizYP13+8meIb9y/j6lQiH9I8qeBqAvXLCEvK53f7ergtQmaPT2zPk3q8LSPQETqgLOAJ7yMwwTXpppiNtUUs+/oIL9saOX959WxsjzP67BckZOZzls2V/PbXe2MTCRm9FC0Pk1q8SwRiEgBcBfwcVUdmOP160SkQUQaurq63A/QBMIjTd080tTNl+9tJD87g4++YY3XIbnqXfUrGJ6Y5t7njiTketH6NKnFk5nFIpJJOAn8TFXvnuscVd0GbAOor69XF8MzAfLtPzYxMDpJ45FBbrxsfdIsNx2rc+pKOWVJPv/5RDPv2LIMEVnU9b79xybAdipLNV6MGhLgB0Cjqn7d7fJN8mnpHWFZSS7XXFDndSiuExGuvaCOp1v6eOJFmxFs4uNF09CFwPuBN4jIzsifN3sQh0kCA6OTDI9P87evX01OZrrX4XjiqnNWsKQgm//3x/1eh2ICyvVEoKqPqKqo6umqembkz71ux2GSw+G+UTLThXdsWe51KJ7JyUzn+otW8cj+bu7faxPMzMLZzGITWM+29jEwNkV1cU7KPg1EXXNBHeuqCvnCr56ja3Dc63BMwFgiMIH17w/uJz8rne+872yvQ/FcZnoaX7/qDPpHJ/mrHzfQNzIR13W+dOVmvnTl5gRHZ/zOEoEJpP2dg/z3nqN88NWnsHmZjXuH8ByAb79nC40dA7zjlj+z+3D/gq+xuqKA1RUFDkRn/MwSgQmkWx48SE5mGquW5POAtYu/5NKNVfz4g+cyODbF277zKJ/85bO82D0c888/sPeo1WcKsh3KTOC09Y7w652Hufq8Wn7+VCsAl2xMvg1o4nXeqnLuv+Eibt7exM+eaOaup9vYsrKE81aVU1mYzVRI6R+dpGd44qU/wxNTVBfn0tg+QFl+ltVnirFEYALnlgcPkCbChy9axcd/vtPrcHypOC+Tf7p8I9e/bhV3NLTx+90dbHvoIFOh8NzMNIHSvCxK87Moy89iSUE2h7qHaesbpa1vlPf/4AluunwjayoLPb4T4wZLBCZQjvSPcUdDG++sX051ca7X4fheZWEOH3n9Gj7y+jVMTIUYHJskIy2NgpwM0tNeOQv5yn9/lK7BcXa19fPmmx/hpis28r5X1XoQuXGT9RGYQPnuQwcIqfI3F632OpTAycpIo7wgm+K8zDmTAIRHH9WU5PLAJy7igjXlfP6/dvPPv9lDKGSrvCQzSwTWn970AAAJxklEQVQmMFp7RvjZEy28/axlrChLjRVGvVJRmM0PrjmHD154Cj/68yFuvPs5SwZJzJqGTGB86d5G0kX45Iz9Br5x1ZkeRpR8ZtZneprwj2/dQEF2Ot/6436mVfnqO04/4dOECS5LBCYQ/rDnCL/ffYRPXHoqS4tzXjpeU2L9BIk0uz5FhE+8cR1pacI3H2gipMq/vfMMSwZJxhKB8b0j/WN85q5dbKop4vpZfQP3PNsOwOVn1HgRWtI5UX1+/JJTEYRvPLCPNBF7MkgylgiMrw2MTXLtD59kfCrEze8+k6yM47u1fvp4M2CJIFFOVp8fu2QtIVVu3t6EAF99x+mkWTJICpYIjG8d6R/jQ7c9xYGuIW699hwb0+4DN1x6Kgp8a3sTIvCVKy0ZJANLBMaXHmnq5pN37GRobIptf1nPa9YmZnN2s3g3XLIWVPnWH/ejCl++cjMZ6TYAMcgsERhfGZmY4iu/f54fP9bMqop8fvSBc9lQXeR1WGYGEeGGS08FEb61vYmjg+N8571nUZiT6XVoJk6WCIxvPLSvi3/89W6aj43wgQvr+MzW9Sm/z4BfiQifuPRUqotz+MKvdvOOW/7Md967hbVV1nwXRKLq/0ki9fX12tDQ4HUYxiGdA2P87981cs+z7ZyyJJ//8/bTuGB1bJun9wyH190vS7FN650ST30+0tTN3//8GYbHp/jCWzbwvlfVWr+BT4jIDlWtn/c8SwTGK6MT09z66Iv8x4MHGJ8K8bevX831F6Xu3sNB1jk4xqfu2MVD+7o4c0UJN12+kbNWlnodVsrzdSIQka3AzUA68H1V/crJzrdEkFyGxqe4o6GVWx48QOfgOBevr+Tzb9nAqjg2RLmjIbwM9bvqVyQ6zJS0mPoMhZS7nm7ja//9Al2D41ywupy/es0pvHZthXUmeyTWROB6H4GIpAPfAS4F2oCnROQ3qrrX7ViMe1SV3YcH+NXOw/zyqVYGx6c4p66U77xvC+fUlcV93Tt3tAGWCBJlMfWZlia8q34Fl22u5mePN/PDRw/xwR81sKQgi7dsrubNm6vZUltKpiUF3/Gis/hcYL+qHgQQkZ8DbwMsEXhMVZkKKdMhJU2E9DQhTcIdgwsRCindw+O09oxwoGuYHYd6+fPBblp7RslIEy7bXM0HL6yzpoMkVZCdwYcvWs0HLjyF7Y1H+c2z7dz+VCu3PdZMQXYG560qp76ulHVLC1lbWcCSgmxrDvSYF4lgGdA64/s24FVOFPSt7U38JjJlPtoEdlxDmPKKY3OdN7P1TCOvHHdsjta1mU1uOsd5c15nzjL1FcfmiveE5ZzkvgGmZ7z5T59gdcn0tHBSSBchI01IS5v1twgi4Tb/kYlpRienj/v5krxM6mtL+bvXr+FNm5ZSkmcdu6kgKyONyzZXc9nmagbHJnl0fzcPN4X/PNB4/HaYBdkZ5GSmk5kuZKQLmWlpLPDzxyss9APMK35+ccUnzJeu3Lyop+ZYeJEI5qrfV7wDich1wHUAK1eujKugysJs1s0cziavDCD6y3L8sROfd9zxGSdI5Bs57thc15RXHjuuRuZ4PcbrzDx3rv8EMke8aQIZ6WlkRN7sM9PDb/CqMDWtTKsyHQoxHYLpUIipkBIKhZNHSDV8Tiic1nKz0snPSicvK4Oy/CxWluVRW55HXXm+jSJJcYU5mWw9rZqtp1UD0D8yyfNHBjjYPUzP8ATdQ+OMTYaYmg7/jk1Oh175prAQi+z61MVeIIFyXXha8iIRtAEzGyCXA+2zT1LVbcA2CHcWx1PQu89dybvPjS+JGGOcU5yXyatWlfOqVeVeh2LwYNSQiGQA+4CLgcPAU8B7VXXPiX7GRg2ZExmdCDdD5WZZG3MiWH0mF9+OGlLVKRH5O+C/CQ8fvfVkScCYk7E3rMSy+kxNniwxoar3Avd6UbZJLj957BAA7z+/zsswkobVZ2qyAb0m0H67q4Pf7urwOoykYfWZmiwRGGNMirNEYIwxKc4SgTHGpDhLBMYYk+ICsQy1iHQBzXO8tATodjkcNyXz/dm9BZPdW7DUquq8+7wGIhGciIg0xDJZIqiS+f7s3oLJ7i05WdOQMcakOEsExhiT4oKeCLZ5HYDDkvn+7N6Cye4tCQW6j8AYY8ziBf2JwBhjzCIFLhGISJmI3C8iTZG/X7HfoYjUisgOEdkpIntE5HovYl2oGO/tTBF5LHJfu0TkKi9iXahY7i1y3n0i0iciv3U7xoUSka0i8oKI7BeRz87xeraI/CLy+hMiUud+lPGJ4d5eKyJPi8iUiLzTixjjFcO9fUJE9kb+f20XkVov4nRT4BIB8Flgu6quBbZHvp+tA7hAVc8kvA3mZ0WkxsUY4xXLvY0Af6mqm4CtwDdFpMTFGOMVy70B/BvwfteiipOIpAPfAS4DNgLvEZGNs077ENCrqmuAbwBfdTfK+MR4by3AtcB/uhvd4sR4b88A9ap6OnAn8DV3o3RfEBPB24DbIl/fBvzF7BNUdUJVxyPfZhOc+4zl3vapalPk63agE5h3wogPzHtvAKq6HRh0K6hFOBfYr6oHVXUC+Dnhe5xp5j3fCVwsi91I1x3z3puqHlLVXUDIiwAXIZZ7+x9VHYl8+zjhXRSTWlDeIGeqUtUOgMjflXOdJCIrRGQX0Ap8NfKm6Xcx3VuUiJwLZAEHXIhtsRZ0bwGwjPDvVlRb5Nic56jqFNAPBGFvxljuLagWem8fAn7vaEQ+4MnGNPMRkQeApXO89PlYr6GqrcDpkSahX4nInap6NFExxisR9xa5TjXwE+AaVfXFp7JE3VtAzPXJfvYQvFjO8aOgxh2LmO9NRK4G6oGLHI3IB3yZCFT1khO9JiJHRaRaVTsib4ad81yrXUT2AK8h/HjuqUTcm4gUAb8DvqCqjzsU6oIl8t8tANqAFTO+Xw7MfuqMntMW2au7GOhxJ7xFieXegiqmexORSwh/gLloRjNz0gpi09BvgGsiX18D/Hr2CSKyXERyI1+XAhcCL7gWYfxiubcs4L+AH6vqHS7Gtljz3lvAPAWsFZFTIv8m7yZ8jzPNvOd3An/UYEzcieXegmreexORs4DvAleoatA/sMRGVQP1h3Ab63agKfJ3WeR4PfD9yNeXAruAZyN/X+d13Am8t6uBSWDnjD9neh17Iu4t8v3DQBcwSvjT25u8jv0k9/RmYB/hPprPR459kfAbCEAOcAewH3gSWOV1zAm8t3Mi/z7DwDFgj9cxJ/DeHgCOzvj/9RuvY3b6j80sNsaYFBfEpiFjjDEJZInAGGNSnCUCY4xJcZYIjDEmxVkiMMaYFGeJwBhjUpwlAmOMSXGWCIwxJsX9fwUM0/0+fWiZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "coefmat.area.plot.density()\n", "xci = coefmat.area.quantile((0.025,0.975)).ravel()\n", "plt.axvline(x=xci[0], linestyle='--')\n", "plt.axvline(x=xci[1], linestyle='--') \n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "application/json": { "Software versions": [ { "module": "Python", "version": "3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]" }, { "module": "IPython", "version": "6.5.0" }, { "module": "OS", "version": "Darwin 17.7.0 x86_64 i386 64bit" }, { "module": "pandas", "version": "0.23.4" }, { "module": "numpy", "version": "1.15.1" }, { "module": "matplotlib", "version": "2.2.3" }, { "module": "seaborn", "version": "0.9.0" }, { "module": "scipy", "version": "1.1.0" }, { "module": "patsy", "version": "0.5.0" }, { "module": "statsmodels", "version": "0.9.0" } ] }, "text/html": [ "
SoftwareVersion
Python3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]
IPython6.5.0
OSDarwin 17.7.0 x86_64 i386 64bit
pandas0.23.4
numpy1.15.1
matplotlib2.2.3
seaborn0.9.0
scipy1.1.0
patsy0.5.0
statsmodels0.9.0
Fri Sep 07 15:13:43 2018 BST
" ], "text/latex": [ "\\begin{tabular}{|l|l|}\\hline\n", "{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n", "Python & 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE\\_401/final)] \\\\ \\hline\n", "IPython & 6.5.0 \\\\ \\hline\n", "OS & Darwin 17.7.0 x86\\_64 i386 64bit \\\\ \\hline\n", "pandas & 0.23.4 \\\\ \\hline\n", "numpy & 1.15.1 \\\\ \\hline\n", "matplotlib & 2.2.3 \\\\ \\hline\n", "seaborn & 0.9.0 \\\\ \\hline\n", "scipy & 1.1.0 \\\\ \\hline\n", "patsy & 0.5.0 \\\\ \\hline\n", "statsmodels & 0.9.0 \\\\ \\hline\n", "\\hline \\multicolumn{2}{|l|}{Fri Sep 07 15:13:43 2018 BST} \\\\ \\hline\n", "\\end{tabular}\n" ], "text/plain": [ "Software versions\n", "Python 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)]\n", "IPython 6.5.0\n", "OS Darwin 17.7.0 x86_64 i386 64bit\n", "pandas 0.23.4\n", "numpy 1.15.1\n", "matplotlib 2.2.3\n", "seaborn 0.9.0\n", "scipy 1.1.0\n", "patsy 0.5.0\n", "statsmodels 0.9.0\n", "Fri Sep 07 15:13:43 2018 BST" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%load_ext version_information\n", "%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels" ] } ], "metadata": { "anaconda-cloud": {}, "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.7.0" } }, "nbformat": 4, "nbformat_minor": 1 }