{
"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",
" Species | \n",
" Area | \n",
" Elevation | \n",
" Nearest | \n",
" Scruz | \n",
" Adjacent | \n",
"
\n",
" \n",
" \n",
" \n",
" Baltra | \n",
" 58 | \n",
" 25.09 | \n",
" 346 | \n",
" 0.6 | \n",
" 0.6 | \n",
" 1.84 | \n",
"
\n",
" \n",
" Bartolome | \n",
" 31 | \n",
" 1.24 | \n",
" 109 | \n",
" 0.6 | \n",
" 26.3 | \n",
" 572.33 | \n",
"
\n",
" \n",
" Caldwell | \n",
" 3 | \n",
" 0.21 | \n",
" 114 | \n",
" 2.8 | \n",
" 58.7 | \n",
" 0.78 | \n",
"
\n",
" \n",
" Champion | \n",
" 25 | \n",
" 0.10 | \n",
" 46 | \n",
" 1.9 | \n",
" 47.4 | \n",
" 0.18 | \n",
"
\n",
" \n",
" Coamano | \n",
" 2 | \n",
" 0.05 | \n",
" 77 | \n",
" 1.9 | \n",
" 1.9 | \n",
" 903.82 | \n",
"
\n",
" \n",
"
\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",
"OLS Regression Results\n",
"\n",
" Dep. Variable: | Species | R-squared: | 0.766 | \n",
"
\n",
"\n",
" Model: | OLS | Adj. R-squared: | 0.717 | \n",
"
\n",
"\n",
" Method: | Least Squares | F-statistic: | 15.70 | \n",
"
\n",
"\n",
" Date: | Fri, 07 Sep 2018 | Prob (F-statistic): | 6.84e-07 | \n",
"
\n",
"\n",
" Time: | 15:12:52 | Log-Likelihood: | -162.54 | \n",
"
\n",
"\n",
" No. Observations: | 30 | AIC: | 337.1 | \n",
"
\n",
"\n",
" Df Residuals: | 24 | BIC: | 345.5 | \n",
"
\n",
"\n",
" Df Model: | 5 | | | \n",
"
\n",
"\n",
" Covariance Type: | nonrobust | | | \n",
"
\n",
"
\n",
"\n",
"\n",
" | coef | std err | t | P>|t| | [0.025 | 0.975] | \n",
"
\n",
"\n",
" Intercept | 7.0682 | 19.154 | 0.369 | 0.715 | -32.464 | 46.601 | \n",
"
\n",
"\n",
" Area | -0.0239 | 0.022 | -1.068 | 0.296 | -0.070 | 0.022 | \n",
"
\n",
"\n",
" Elevation | 0.3195 | 0.054 | 5.953 | 0.000 | 0.209 | 0.430 | \n",
"
\n",
"\n",
" Nearest | 0.0091 | 1.054 | 0.009 | 0.993 | -2.166 | 2.185 | \n",
"
\n",
"\n",
" Scruz | -0.2405 | 0.215 | -1.117 | 0.275 | -0.685 | 0.204 | \n",
"
\n",
"\n",
" Adjacent | -0.0748 | 0.018 | -4.226 | 0.000 | -0.111 | -0.038 | \n",
"
\n",
"
\n",
"\n",
"\n",
" Omnibus: | 12.683 | Durbin-Watson: | 2.476 | \n",
"
\n",
"\n",
" Prob(Omnibus): | 0.002 | Jarque-Bera (JB): | 13.498 | \n",
"
\n",
"\n",
" Skew: | 1.136 | Prob(JB): | 0.00117 | \n",
"
\n",
"\n",
" Kurtosis: | 5.374 | Cond. No. | 1.90e+03 | \n",
"
\n",
"
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",
" df_resid | \n",
" ssr | \n",
" df_diff | \n",
" ss_diff | \n",
" F | \n",
" Pr(>F) | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 25.0 | \n",
" 93469.083990 | \n",
" 0.0 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" 1 | \n",
" 24.0 | \n",
" 254156.680724 | \n",
" 1.0 | \n",
" -160687.596734 | \n",
" -15.17372 | \n",
" 1.0 | \n",
"
\n",
" \n",
"
\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",
" df_resid | \n",
" ssr | \n",
" df_diff | \n",
" ss_diff | \n",
" F | \n",
" Pr(>F) | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 26.0 | \n",
" 158291.628568 | \n",
" 0.0 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" 1 | \n",
" 24.0 | \n",
" 89231.366330 | \n",
" 2.0 | \n",
" 69060.262238 | \n",
" 9.287352 | \n",
" 0.00103 | \n",
"
\n",
" \n",
"
\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",
" df_resid | \n",
" ssr | \n",
" df_diff | \n",
" ss_diff | \n",
" F | \n",
" Pr(>F) | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 25.0 | \n",
" 109591.120801 | \n",
" 0.0 | \n",
" NaN | \n",
" NaN | \n",
" NaN | \n",
"
\n",
" \n",
" 1 | \n",
" 24.0 | \n",
" 89231.366330 | \n",
" 1.0 | \n",
" 20359.754471 | \n",
" 5.476035 | \n",
" 0.027926 | \n",
"
\n",
" \n",
"
\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",
" 0 | \n",
" 1 | \n",
"
\n",
" \n",
" \n",
" \n",
" Intercept | \n",
" -22.800754 | \n",
" 127.659321 | \n",
"
\n",
" \n",
" Area | \n",
" -0.077134 | \n",
" 0.079456 | \n",
"
\n",
" \n",
" Elevation | \n",
" -0.224103 | \n",
" 0.156443 | \n",
"
\n",
" \n",
" Nearest | \n",
" 1.087895 | \n",
" 6.825321 | \n",
"
\n",
" \n",
" Scruz | \n",
" -0.786385 | \n",
" 0.434932 | \n",
"
\n",
" \n",
" Adjacent | \n",
" -0.005252 | \n",
" 0.121375 | \n",
"
\n",
" \n",
"
\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",
" intercept | \n",
" area | \n",
" elevation | \n",
" nearest | \n",
" Scruz | \n",
" adjacent | \n",
"
\n",
" \n",
" \n",
" \n",
" 0.025 | \n",
" -7.527822 | \n",
" -0.064221 | \n",
" -0.192220 | \n",
" 1.710021 | \n",
" -0.633871 | \n",
" 0.005495 | \n",
"
\n",
" \n",
" 0.975 | \n",
" 120.565976 | \n",
" 0.088243 | \n",
" 0.140516 | \n",
" 6.683928 | \n",
" 0.450277 | \n",
" 0.124617 | \n",
"
\n",
" \n",
"
\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": [
"Software | Version |
---|
Python | 3.7.0 64bit [Clang 4.0.1 (tags/RELEASE_401/final)] |
IPython | 6.5.0 |
OS | Darwin 17.7.0 x86_64 i386 64bit |
pandas | 0.23.4 |
numpy | 1.15.1 |
matplotlib | 2.2.3 |
seaborn | 0.9.0 |
scipy | 1.1.0 |
patsy | 0.5.0 |
statsmodels | 0.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
}