{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Chapter 10: Model Selection\n", "\n", "Load the packages:" ] }, { "cell_type": "code", "execution_count": 1, "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": [ "## Testing based procedures\n", "\n", "Load in the data:" ] }, { "cell_type": "code", "execution_count": 2, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PopulationIncomeIlliteracyLifeExpMurderHSGradFrostArea
AL361536242.169.0515.141.32050708
AK36563151.569.3111.366.7152566432
AZ221245301.870.557.858.115113417
AR211033781.970.6610.139.96551945
CA2119851141.171.7110.362.620156361
\n", "
" ], "text/plain": [ " Population Income Illiteracy LifeExp Murder HSGrad Frost Area\n", "AL 3615 3624 2.1 69.05 15.1 41.3 20 50708\n", "AK 365 6315 1.5 69.31 11.3 66.7 152 566432\n", "AZ 2212 4530 1.8 70.55 7.8 58.1 15 113417\n", "AR 2110 3378 1.9 70.66 10.1 39.9 65 51945\n", "CA 21198 5114 1.1 71.71 10.3 62.6 20 156361" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "statedata = pd.read_csv(\"data/statedata.csv\",index_col=0)\n", "statedata.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the base 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: LifeExp R-squared: 0.736
Model: OLS Adj. R-squared: 0.692
Method: Least Squares F-statistic: 16.74
Date: Wed, 26 Sep 2018 Prob (F-statistic): 2.53e-10
Time: 10:37:14 Log-Likelihood: -51.855
No. Observations: 50 AIC: 119.7
Df Residuals: 42 BIC: 135.0
Df Model: 7
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", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 70.9432 1.748 40.586 0.000 67.416 74.471
Population 5.18e-05 2.92e-05 1.775 0.083 -7.1e-06 0.000
Income -2.18e-05 0.000 -0.089 0.929 -0.001 0.000
Illiteracy 0.0338 0.366 0.092 0.927 -0.705 0.773
Murder -0.3011 0.047 -6.459 0.000 -0.395 -0.207
HSGrad 0.0489 0.023 2.098 0.042 0.002 0.096
Frost -0.0057 0.003 -1.825 0.075 -0.012 0.001
Area -7.383e-08 1.67e-06 -0.044 0.965 -3.44e-06 3.29e-06
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 2.385 Durbin-Watson: 1.929
Prob(Omnibus): 0.303 Jarque-Bera (JB): 1.420
Skew: -0.081 Prob(JB): 0.492
Kurtosis: 2.190 Cond. No. 1.85e+06


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.85e+06. This might indicate that there are
strong multicollinearity or other numerical problems." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: LifeExp R-squared: 0.736\n", "Model: OLS Adj. R-squared: 0.692\n", "Method: Least Squares F-statistic: 16.74\n", "Date: Wed, 26 Sep 2018 Prob (F-statistic): 2.53e-10\n", "Time: 10:37:14 Log-Likelihood: -51.855\n", "No. Observations: 50 AIC: 119.7\n", "Df Residuals: 42 BIC: 135.0\n", "Df Model: 7 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 70.9432 1.748 40.586 0.000 67.416 74.471\n", "Population 5.18e-05 2.92e-05 1.775 0.083 -7.1e-06 0.000\n", "Income -2.18e-05 0.000 -0.089 0.929 -0.001 0.000\n", "Illiteracy 0.0338 0.366 0.092 0.927 -0.705 0.773\n", "Murder -0.3011 0.047 -6.459 0.000 -0.395 -0.207\n", "HSGrad 0.0489 0.023 2.098 0.042 0.002 0.096\n", "Frost -0.0057 0.003 -1.825 0.075 -0.012 0.001\n", "Area -7.383e-08 1.67e-06 -0.044 0.965 -3.44e-06 3.29e-06\n", "==============================================================================\n", "Omnibus: 2.385 Durbin-Watson: 1.929\n", "Prob(Omnibus): 0.303 Jarque-Bera (JB): 1.420\n", "Skew: -0.081 Prob(JB): 0.492\n", "Kurtosis: 2.190 Cond. No. 1.85e+06\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.85e+06. 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('LifeExp ~ Population + Income + Illiteracy + Murder + HSGrad + Frost + Area', data=statedata).fit()\n", "lmod.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Backward elimination using p-values:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 0.00\n", "Population 0.08\n", "Income 0.93\n", "Illiteracy 0.93\n", "Murder 0.00\n", "HSGrad 0.04\n", "Frost 0.08\n", "Area 0.96\n", "dtype: float64" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod.pvalues.round(2)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 0.00\n", "Population 0.08\n", "Income 0.92\n", "Illiteracy 0.93\n", "Murder 0.00\n", "HSGrad 0.02\n", "Frost 0.06\n", "dtype: float64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Population + Income + Illiteracy + Murder + HSGrad + Frost', data=statedata).fit()\n", "lmod.pvalues.round(2)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 0.00\n", "Population 0.07\n", "Income 0.92\n", "Murder 0.00\n", "HSGrad 0.01\n", "Frost 0.02\n", "dtype: float64" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Population + Income + Murder + HSGrad + Frost', data=statedata).fit()\n", "lmod.pvalues.round(2)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 0.00\n", "Population 0.05\n", "Murder 0.00\n", "HSGrad 0.00\n", "Frost 0.02\n", "dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Population + Murder + HSGrad + Frost', data=statedata).fit()\n", "lmod.pvalues.round(2)" ] }, { "cell_type": "code", "execution_count": 8, "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: LifeExp R-squared: 0.713
Model: OLS Adj. R-squared: 0.694
Method: Least Squares F-statistic: 38.03
Date: Wed, 26 Sep 2018 Prob (F-statistic): 1.63e-12
Time: 10:37:14 Log-Likelihood: -53.987
No. Observations: 50 AIC: 116.0
Df Residuals: 46 BIC: 123.6
Df Model: 3
Covariance Type: nonrobust
\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 71.0364 0.983 72.246 0.000 69.057 73.016
Murder -0.2831 0.037 -7.706 0.000 -0.357 -0.209
HSGrad 0.0499 0.015 3.286 0.002 0.019 0.081
Frost -0.0069 0.002 -2.824 0.007 -0.012 -0.002
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 5.102 Durbin-Watson: 1.794
Prob(Omnibus): 0.078 Jarque-Bera (JB): 2.569
Skew: -0.290 Prob(JB): 0.277
Kurtosis: 2.053 Cond. No. 1.19e+03


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.19e+03. This might indicate that there are
strong multicollinearity or other numerical problems." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: LifeExp R-squared: 0.713\n", "Model: OLS Adj. R-squared: 0.694\n", "Method: Least Squares F-statistic: 38.03\n", "Date: Wed, 26 Sep 2018 Prob (F-statistic): 1.63e-12\n", "Time: 10:37:14 Log-Likelihood: -53.987\n", "No. Observations: 50 AIC: 116.0\n", "Df Residuals: 46 BIC: 123.6\n", "Df Model: 3 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 71.0364 0.983 72.246 0.000 69.057 73.016\n", "Murder -0.2831 0.037 -7.706 0.000 -0.357 -0.209\n", "HSGrad 0.0499 0.015 3.286 0.002 0.019 0.081\n", "Frost -0.0069 0.002 -2.824 0.007 -0.012 -0.002\n", "==============================================================================\n", "Omnibus: 5.102 Durbin-Watson: 1.794\n", "Prob(Omnibus): 0.078 Jarque-Bera (JB): 2.569\n", "Skew: -0.290 Prob(JB): 0.277\n", "Kurtosis: 2.053 Cond. No. 1.19e+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.19e+03. This might indicate that there are\n", "strong multicollinearity or other numerical problems.\n", "\"\"\"" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Murder + HSGrad + Frost', data=statedata).fit()\n", "lmod.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Variables that are eliminated may have some significance." ] }, { "cell_type": "code", "execution_count": 9, "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: LifeExp R-squared: 0.674
Model: OLS Adj. R-squared: 0.653
Method: Least Squares F-statistic: 31.69
Date: Wed, 26 Sep 2018 Prob (F-statistic): 2.91e-11
Time: 10:37:14 Log-Likelihood: -57.147
No. Observations: 50 AIC: 122.3
Df Residuals: 46 BIC: 129.9
Df Model: 3
Covariance Type: nonrobust
\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 74.5567 0.584 127.611 0.000 73.381 75.733
Illiteracy -0.6018 0.299 -2.013 0.050 -1.203 -5.18e-05
Murder -0.2800 0.043 -6.454 0.000 -0.367 -0.193
Frost -0.0087 0.003 -2.937 0.005 -0.015 -0.003
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
Omnibus: 0.019 Durbin-Watson: 1.638
Prob(Omnibus): 0.990 Jarque-Bera (JB): 0.071
Skew: 0.026 Prob(JB): 0.965
Kurtosis: 2.823 Cond. No. 638.


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: LifeExp R-squared: 0.674\n", "Model: OLS Adj. R-squared: 0.653\n", "Method: Least Squares F-statistic: 31.69\n", "Date: Wed, 26 Sep 2018 Prob (F-statistic): 2.91e-11\n", "Time: 10:37:14 Log-Likelihood: -57.147\n", "No. Observations: 50 AIC: 122.3\n", "Df Residuals: 46 BIC: 129.9\n", "Df Model: 3 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 74.5567 0.584 127.611 0.000 73.381 75.733\n", "Illiteracy -0.6018 0.299 -2.013 0.050 -1.203 -5.18e-05\n", "Murder -0.2800 0.043 -6.454 0.000 -0.367 -0.193\n", "Frost -0.0087 0.003 -2.937 0.005 -0.015 -0.003\n", "==============================================================================\n", "Omnibus: 0.019 Durbin-Watson: 1.638\n", "Prob(Omnibus): 0.990 Jarque-Bera (JB): 0.071\n", "Skew: 0.026 Prob(JB): 0.965\n", "Kurtosis: 2.823 Cond. No. 638.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Illiteracy + Murder + Frost', data=statedata).fit()\n", "lmod.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Criterion-based procedures.\n", "\n", "Can use `sklearn` but only makes sense if data are scaled. Note this method is not in \"Linear Models with R\"." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/anaconda/lib/python3.7/site-packages/sklearn/utils/__init__.py:4: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working\n", " from collections import Sequence\n" ] }, { "data": { "text/plain": [ "Population 0.172276\n", "Income -0.009981\n", "Illiteracy 0.015357\n", "Murder -0.828079\n", "HSGrad 0.294402\n", "Frost -0.222074\n", "Area -0.004693\n", "dtype: float64" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.preprocessing import scale\n", "scalstat = pd.DataFrame(scale(statedata), index=statedata.index, columns=statedata.columns)\n", "from sklearn import linear_model\n", "reg = linear_model.LinearRegression(fit_intercept=False)\n", "X = scalstat.drop('LifeExp',axis=1)\n", "smlmod = sm.OLS(scalstat.LifeExp, X).fit()\n", "smlmod.params" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.17227607, -0.00998072, 0.0153566 , -0.82807917, 0.29440196,\n", " -0.22207364, -0.004693 ])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "reg.fit(X, scalstat.LifeExp)\n", "reg.coef_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Expect the intercept to be zero since was not included." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "reg.intercept_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Chooses importance of predictors based on the coefficient size:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4, 6, 5, 1, 2, 3, 7])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.feature_selection import RFE\n", "selector = RFE(reg, 1, step=1)\n", "selector = selector.fit(X, scalstat.LifeExp)\n", "selector.ranking_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order of importance:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['Murder', 'Frost', 'HSGrad', 'Population', 'Income', 'Illiteracy',\n", " 'Area'],\n", " dtype='object')" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X.columns[selector.ranking_-1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use 10-fold crossvalidation" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2, 4, 3, 1, 1, 1, 5])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.feature_selection import RFECV\n", "selector = RFECV(reg, step=1, cv=10)\n", "selector = selector.fit(X, scalstat.LifeExp)\n", "selector.ranking_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Picks only 3 predictors." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['Murder', 'HSGrad', 'Frost'], dtype='object')" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X.columns[selector.support_]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### All subsets regression\n", "\n", "Compute the RSS for every subset of predictors." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "import itertools" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([88.299, 34.461, 29.77 , 25.372, 23.308, 23.302, 23.298, 23.297])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pcols = list(statedata.columns)\n", "pcols.remove('LifeExp')\n", "rss = np.empty(len(pcols) + 1)\n", "rss[0] = np.sum((statedata.LifeExp - np.mean(statedata.LifeExp))**2)\n", "selvar = ['Null']\n", "for k in range(1,len(pcols)+1):\n", " RSS = {}\n", " for variables in itertools.combinations(pcols, k):\n", " predictors = statedata.loc[:,list(variables)]\n", " predictors['Intercept'] = 1\n", " res = sm.OLS(statedata.LifeExp, predictors).fit()\n", " RSS[variables] = res.ssr\n", " rss[k] = min(RSS.values())\n", " selvar.append(min(RSS, key=RSS.get))\n", "rss.round(3)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['Null',\n", " ('Murder',),\n", " ('Murder', 'HSGrad'),\n", " ('Murder', 'HSGrad', 'Frost'),\n", " ('Population', 'Murder', 'HSGrad', 'Frost'),\n", " ('Population', 'Income', 'Murder', 'HSGrad', 'Frost'),\n", " ('Population', 'Income', 'Illiteracy', 'Murder', 'HSGrad', 'Frost'),\n", " ('Population', 'Income', 'Illiteracy', 'Murder', 'HSGrad', 'Frost', 'Area')]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selvar" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8leX9//HXJwlb9t57BYgKAVQQUFBEhgr4q6PuVm2ttVpFqYqIoyLuOuqoOFvrospegiiKElRW2BtlBJAVVsbn90cO36YxHCDJyX2SvJ+Px3mQ3PN9UM7nXNd139dt7o6IiMixxAQdQEREopsKhYiIhKVCISIiYalQiIhIWCoUIiISlgqFiIiEpUIhIiJhBVIozOxSM1tqZplmlpjL+kZmtt/M7gwin4iI/FdQLYolwGBgzjHWPw1MLrw4IiJyLHFBnNTdlwGY2S/WmdnFwFog9USPV6NGDW/SpElBxRMRKREWLFiww91rHm+7QArFsZhZBeBu4DwgbLeTmd0I3AjQqFEjkpKSIh9QRKQYMbMNJ7JdxLqezGyGmS3J5XVRmN0eBJ529/3HO767v+Luie6eWLPmcQuiiIjkUcRaFO7eJw+7dQWGmtnjQBUg08wOufvzBZtOREROVFR1Pbn72Ud/NrORwH4VCRGRYAV1eewlZrYZOBOYaGZTg8ghIiLHF9RVT+OAccfZZmThpBERkXB0Z7aIiISlQiEiImGV6EKx+8ARHhy/lD0H04KOIiIStUp0odi46wBvfrWev05aFnQUEZGoVaILRUKDKvy2RzPem7+JL1ftCDqOiEhUKtGFAuD2Pq1oWqMC93y8iNTD6UHHERGJOiW+UJQtFcvjQxP4cfdBxkxdEXQcEZGoU+ILBUDnJtW4+ozGvPn1euav3xV0HBGRqKJCETLsgjbUq1yOuz9cxKG0jKDjiIhEDRWKkApl4nhsSAfW7kjlmRmrgo4jIhI1VCiyObtlTX6V2JBX5qxh0ebdQccREYkKKhQ5/KV/W2pWLMOwDxdxJD0z6DgiIoFTocihcrlSPHxxB5Zv3ceLs1cHHUdEJHAqFLk4L742g06txwuzVrNi676g44iIBEqF4hhGDmpHpbKlGPbhQtIz1AUlIiWXCsUxVKtQmpGD2rFw8x7+8eW6oOOIiARGhSKMAQl1OT++Nk9NX8nalP1BxxERCYQKRRhmxsMXt6dMXAx3f7SIzEwPOpKISKFToTiOWpXKcv+AeOav/5l3vtkQdBwRkUKnQnEChnZqQI9WNXls8nI27ToQdBwRkUIVSKEws0vNbKmZZZpZYo51CWb2dWj9YjMrG0TGHJl49JL2GPCXcYtxVxeUiJQcQbUolgCDgTnZF5pZHPAOcLO7twN6AVHxnNIGVctzT782fLFqBx8kbQ46johIoQmkULj7MnfP7eEP5wOL3H1haLud7h41U7le2bUxXZpW46GJyWzbeyjoOCIihSLaxihaAW5mU83sOzMbFnSg7GJijNFDEjiSnsm945aoC0pESoSIFQozm2FmS3J5XRRmtzigO3Bl6M9LzKz3MY5/o5klmVlSSkpKBN5B7prWqMCfz2/FjGXbGL9oS6GdV0QkKHGROrC798nDbpuBz919B4CZTQI6AjNzOf4rwCsAiYmJhfrV/obuzZi4eCsjP11Kt+bVqX5KmcI8vYhIoYq2rqepQIKZlQ8NbPcEkgPO9AuxMcaYoQnsO5TGyPFRF09EpEAFdXnsJWa2GTgTmGhmUwHc/WfgKWA+8APwnbtPDCLj8bSqXZFbz23J+IU/MW3p1qDjiIhEjBWHAdnExERPSkoq9POmZWQy8G9fsiv1CNPv6EnlcqUKPYOISF6Z2QJ3TzzedtHW9VSklIqNYczQU9mZeoRHJqoLSkSKJxWKfOrQoDI39mjG+0mbmbOy8K6+EhEpLCoUBeC23i1pVrMCwz9ezP7D6UHHEREpUCoUBaBsqVjGDE3gpz0HeXzK8qDjiIgUKBWKAtKpcTWuPasJb329gW/X7Qo6johIgVGhKEB39W1Nw2rluPujRRxKi5opqkRE8kWFogCVLx3HY4MTWLcjlaenrww6johIgVChKGDdWtTg8i4NefWLtSzctDvoOCIi+aZCEQHDL2xLrYpluevDhRxOVxeUiBRtKhQRUKlsKR4d3J6V2/bzwqw1QccREckXFYoIObdNbS45vT4vzlrNsi17g44jIpJnKhQRNGJAPFXKl2LYh4tIz8gMOo6ISJ6oUERQ1QqlGXVRexb/uIdXv1gXdBwRkTxRoYiwCzvU5YJ2dXh6xkrWpOwPOo6IyElToSgEoy5uR7lSsQz7cBEZmUV/WncRKVlUKApBrYplGTEgngUbfuatr9cHHUdE5KSoUBSSwR3r06t1TR6fsoJNuw4EHUdE5ISpUBQSM+PRSzoQG2Pc8/EiisOTBUWkZFChKET1qpRj+IVtmLt6J/+evynoOCIiJ0SFopBd3rkRZzSrxiMTl7F1z6Gg44iIHFcghcLMLjWzpWaWaWaJ2ZaXMrM3zWyxmS0zs+FB5IukmBhj9JAE0jIzuXfcYnVBiUjUC6pFsQQYDMzJsfxSoIy7dwA6ATeZWZPCjRZ5jatX4M7zWzNz+XY+XfhT0HFERMIKpFC4+zJ3X5HbKqCCmcUB5YAjQLGcKOm6bk05vVEVRn66lB37DwcdR0TkmKJtjOJDIBXYAmwEnnD3Yvlc0dgY4/EhCaQezuCBT5cGHUdE5JgiVijMbIaZLcnldVGY3boAGUA9oCnwZzNrdozj32hmSWaWlJKSEoF3EHkta1fkj71bMHHRFqYs2Rp0HBGRXMVF6sDu3icPu10BTHH3NGC7mc0FEoG1uRz/FeAVgMTExCI7InxTz+ZMWryV+z9ZwhnNqlGlfOmgI4mI/I9o63raCJxrWSoAZwDLA84UUaViY3h8aAK7Uo/w0IRlQccREfmFoC6PvcTMNgNnAhPNbGpo1QvAKWRdFTUfGOvui4LIWJja16/M73o256PvNjN7xfag44iI/A8rDtfxJyYmelJSUtAx8uVwegb9n/uSA4fTmXp7DyqWLRV0JBEp5sxsgbsnHm+7aOt6KrHKxMXy+NAEtuw9xOgpxbq3TUSKGBWKKNKxUVWu79aUd+Zt5Os1O4OOIyICqFBEnTvPb03j6uW55+NFHDySEXQcEREVimhTrnQsjw1OYMPOAzw5Lbeb10VECpcKRRQ6s3l1ruzaiNfnruP7jT8HHUdESjgViih1T7821KlUlmEfLuJwurqgRCQ4KhRRqmLZUjwyuAOrtu/n+c9WBx1HREowFYoodk7rWgzuWJ+XZq9h6U97go4jIiWUCkWUGzEgnirlSzPsw0WkZWQGHUdESiAViihXpXxpHr64HUt/2ssrc34xN6KISMSpUBQBF7SvS/8OdXl2xipWb98XdBwRKWFUKIqIkYPaUb5MLMM+XERGZtGfn0tEig4ViiKiZsUyjBzYju827uaNr9YHHUdEShAViiLkotPqcW6bWoyZupwNO1ODjiMiJYQKRRFiZjxySXtKxcRwz0eLKQ5TxItI9FOhKGLqVi7HX/q35eu1O/nXt5uCjiMiJYAKRRF0WeeGdGtRnUcnLeOn3QeDjiMixZwKRRFkZjw2OIGMTOfeceqCEpHIUqEoohpWK89dfVsza0UK477/Meg4IlKMqVAUYdec1YROjasyakIyKfsOBx1HRIqpQAqFmY0xs+VmtsjMxplZlWzrhpvZajNbYWZ9g8hXVMTGGKOHJHDgSAYPfLok6DgiUkwF1aKYDrR39wRgJTAcwMzigcuAdsAFwItmFhtQxiKhRa1T+FOflkxavJXJi7cEHUdEiqFACoW7T3P39NCv84AGoZ8vAt5z98Puvg5YDXQJImNRcuPZzWhfvxL3f7KUvYfSgo4jIsVMNIxRXA9MDv1cH8h+c8Dm0DIJIy42hscGJ7Az9TDPzVgVdBwRKWYiVijMbIaZLcnldVG2be4F0oF3jy7K5VC5XvtpZjeaWZKZJaWkpBT8Gyhi2tevzK8SG/LGV+tZk7I/6DgiUoxErFC4ex93b5/L6xMAM7sGGABc6f+9EWAz0DDbYRoAPx3j+K+4e6K7J9asWTNSb6NIubNva8qViuWhCclBRxGRYiSoq54uAO4GBrn7gWyrPgUuM7MyZtYUaAl8G0TGoqjGKWW4rU9LZq9I4bPl24KOIyLFRFBjFM8DFYHpZvaDmf0dwN2XAu8DycAU4BZ3zwgoY5F09ZlNaFazAg9NWMaRdD06VUTyL6irnlq4e0N3Py30ujnbukfcvbm7t3b3yeGOI79UOi6G+wfEs25HKm98tS7oOCJSDETDVU9SwM5pXYtz29TiuZmrdce2iOSbCkUxdV//thxOz2DM1OVBRxGRIk6FophqVvMUruvWlA8WbGbR5t1BxxGRIkyFohi79dwWVK9QmgfHJ2sqchHJMxWKYqxi2VIM69uGBRt+5tOFud6OIiJyXCoUxdzQTg3oUL8yf520nANH0o+/g4hIDioUxVxMjDFyUDxb9x7ipdlrgo4jIkWQCkUJ0KlxNS4+rR4vz1nLpl0Hjr+DiEg2KhQlxN392hBrxqOTlgUdRUSKmGMWCjMra2a/mG3PzGqZWdnIxpKCVrdyOW45pzmTl2zlqzU7go4jIkVIuBbFc8DZuSw/D3g6MnEkkn5zdjMaVC3HqPHJpGdoHigROTHhCkV3d/8450J3fxfoEblIEillS8VyX/+2LN+6j399uzHoOCJSRIQrFLk9ROhE9pMo1rddHc5sVp0np69k94EjQccRkSIg3Af+djP7xfOqzawzoEfKFVFmxgOD4tl7MI2np68MOo6IFAFxYdbdBbxvZm8AC0LLEoGrgcsinEsiqE2dSlzZtTHvfLORK7o2pnWdikFHEpEodswWhbt/C3Qlqwvq2tDLgK7u/k1hhJPIueO8VpxSJo5RE5ZqHigRCStciwJ33wY8UEhZpBBVrVCaO85rxQOfLmVa8jb6tqsTdCQRiVLHLBRmthjI7aumAZnufmrEUkmhuLJrI979ZgMPT0ymZ6ualC0VG3QkEYlC4QazBwADc7wGAb8Dfox8NIm0uNgYHhjYjk27DvKPL/XYVBHJXbgxig1HX0BV4BZgNvAQMKlw4kmkdWtRg77tavPCrNVs3XMo6DgiEoXCTeHRysxGmNky4HlgE2Dufo67P5+fk5rZGDNbbmaLzGycmVUJLT/PzBaY2eLQn+fm5zxyYu69MJ70TGf0FD02VUR+KVzX03KgNzDQ3bu7+9+AjAI673SgvbsnACuB4aHlO0Ln6wBcA7xdQOeTMBpVL89vz27KuO9/ZMGGn4OOIyJRJlyhGAJsBWaZ2atm1pvwd2ufMHef5u5Hn6IzD2gQWv69ux99FNtSoKyZlSmIc0p4v+/VgtqVyjBq/FIyM3W5rIj8V7gxinHu/iugDVljE7cDtc3sJTM7vwAzXA9MzmX5EOB7dz9cgOeSY6hQJo57+rVh4eY9fPTd5qDjiEgUOe6cTe6e6u7vuvsAsr75/wDcc7z9zGyGmS3J5XVRtm3uBdKBd3Ps2w4YDdwU5vg3mlmSmSWlpGhGkYJw8Wn16dioCqOnrGDfobSg44hIlLCg7so1s2uAm4He7n4g2/IGwGfAde4+90SOlZiY6ElJSZEJWsIs3LSbi16Yy009mzG8X9ug44hIGIfSMjiclknl8qXytL+ZLXD3xONtF8gssGZ2AXA3MChHkagCTASGn2iRkIJ1asMqXNqpAa9/uY51O1KDjiMiuXB3Ji/ewnlPf86D45dG/HxBTRf+PFARmG5mP5jZ30PL/wC0AO4PLf/BzGoFlLHEuuuC1pSJi+WRiclBRxGRHJb+tIfLX53H7979jvKl4hjSqUHEzxl2rqdIcfcWx1j+MPBwIceRHGpVLMut57bgr5OX8/nKFHq2+sUTcUWkkO3Yf5gnp63gvfmbqFKuFA9d3J7LOzckLjby3/cDKRQS/a7t1oR/fbuRUeOXMuVPPShVCP8zisgvHUnP5M2v1vPczFUcTMvgurOaclvvlnkel8gL/euXXJWJi+X+AfGsSUnlra83BB1HpMRxd2Ykb6PvM3N4ZNIyEptUZcqfejBiYHyhFglQi0LCOLdNLXq0qskzM1Zy8Wn1qH6K7n0UKQwrt+3joQnJfLFqB81rVuCN6zrTq3Vww7VqUcgxmRkjBrTl4JEMnpimx6aKRNrPqUcY8ckS+j37BQs37eaBgfFM+VOPQIsEqEUhx9GiVkWuPrMJY79ax6/PaES7epWDjiRS7KRlZPLOvA08M2MV+w6l8eszGnN7n1ZUrVA66GiAWhRyAm7r05Kq5Uvz4KfJemyqSAGbvWI7/Z79ggfHJ9OhfmUm39aDURe1j5oiASoUcgIqlyvFnee35tv1u5i4eEvQcUSKhTUp+7lu7LdcO3Y+6RmZvHZ1Im/f0IXWdSoGHe0X1PUkJ+RXnRvyzrwNPDpxGb3b1KZcaT02VSQv9hxI49mZq3jr6/WUKxXLvRe25ZqzmlA6Lnq/t0dvMokqsTHGyEHt+GnPIV6esyboOCJFTnpGJm/P20CvJ2Yx9qt1XJrYgFl39eK3PZpFdZEAtSjkJHRpWo0BCXX5++druDSxIfWrlAs6kkiRMHf1DkaNT2bFtn10bVqNEQPji9SFIdFdxiTqDL8wa0bZv05aFnASkei3fkcqv30riStf+4bUI+m8dGVH3rvxjCJVJEAtCjlJ9auU4+aezXlmxiquOmMnXZtVDzqSSNTZdyiN5z9bzdi564mLNe7q25obujelbKmiObanQiEn7aYezXl//iZGjk9mwq3diY0pkCfkihR5GZnOB0mbeGLaCnbsP8LQTg0Y1rc1tSqVDTpavqjrSU5audKx/KV/W5Zt2cu/528KOo5IVPhm7U4GPf8l93y8mCbVK/DpH7rxxKWnFvkiAWpRSB7171CXt5pu4IlpK+jfoW6hT1ImEi027TrAY5OXM3HxFupVLstzl5/OwIS6mBWflrZaFJInZsYDA+P5+cARnp25Kug4IoUu9XA6T0xdQe+nPmfm8m3c3qcVM//ci0Gn1itWRQLUopB8aFevMpd1bsRbX6/niq4NaVEr+u4oFSlomZnOuO9/ZPSU5Wzfd5iLT6vH3f3aULdy8b1cXC0KyZc7z29FudKxPDhe80BJ8bdgw89c8tJX/PmDhdStUo6Pf38Wz1x2erEuEqAWheRT9VPKcHufVoyakMzMZdvpE1876EgiBe6n3QcZPWU5n/zwE7UqluHJS0/lktPrE1NCrvhToZB8u+rMxvzz2408PDGZs1vVoExc0bxWXCSng0cyeHnOGv7++Rrc4dZzW3Bzz+ZUKFOyPjrV9ST5Vio2hvsHxLN+5wHGzl0fdByRfHN3PvnhR3o/OZtnZqyid5vazLijJ38+v3WJKxIQUKEwszFmttzMFpnZODOrkmN9IzPbb2Z3BpFPTl7PVjXp07YWf5u5iu37DgUdRyTPFm7azdC/f81t7/1A1Qql+feNZ/DClR1pWK180NECE1SLYjrQ3t0TgJXA8BzrnwYmF3oqyZf7+sdzJCOTx6esCDqKyEnbvvcQf35/IRe9MJcNO1MZPaQDn/6hu6apIaAxCneflu3XecDQo7+Y2cXAWiC1sHNJ/jSpUYHruzfl5c/XctUZjTm1YZXj7yQSsENpGfzjy3W8MGs16RnOTT2b8YdzWlCxrG4iPSoaxiiuJ9R6MLMKwN3Ag8fbycxuNLMkM0tKSUmJcEQ5Ubee25KaFcswcvxSMjN1uaxEL3dn0uIt9Hnqc8ZMXUH3FjWYfkcPhvdrqyKRQ8QKhZnNMLMlubwuyrbNvUA68G5o0YPA0+6+/3jHd/dX3D3R3RNr1qwZmTchJ+2UMnEM69ua7zfu5pOFPwYdRyRXS3/aw2WvzOP3737HKWXi+OdvuvLK1Yk0rl4h6GhRKWJdT+7eJ9x6M7sGGAD09v/eqdUVGGpmjwNVgEwzO+Tuz0cqpxS8IR0b8M68DTw2eTnnx9cpkVeJSHRK2XeYp6av4L35m6hSrhQPX9yeyzo3JC42GjpXolcg/4LN7AKyuph6uvuBo8vd/exs24wE9qtIFD0xMcYDg9ox+MWveHH2au7q2yboSFLCHU7P4I256/nbZ6s5lJbB9d2a8sfeLalcTl1MJyKor3rPA2WA6aHJs+a5+80BZZEI6NioKoNPr8+rX6zjV4mNaFS95F5aKMFxd6Ylb+PRScvYsPMAvdvU4t7+bWlW85SgoxUpQV311OIEthlZCFEkgu7u14YpS7fyyKRkXr4qMeg4UsIs27KXhyYk89WanbSsdQpvXd+FHq00npkX6jyWiKldqSy3nNOCMVNXMHf1Drq1qBF0JCkBdu4/zJPTV/LetxupVK4Uoy5qxxVdGmkcIh9UKCSibujelH/P38SD45cy6Y9n6x+rRMyR9Eze/Go9z81cxYG0DK4+swl/6tOSKuVLBx2tyFOhkIgqWyqWe/u35aa3F/DuNxu55qwmQUeSYsbdmbFsO49MTGb9zgP0al2T+/q31fNRCpAKhUTc+fG16daiOk9NX8mgU+tRtYK+4UnBWLF1Hw9NSObL1TtoXrMCY6/rzDmtawUdq9hRP4BEnJkxYkA79h9O56npK4OOI8XArtQj3PefxfR7dg6Lf9zDAwPjmfKnHioSEaIWhRSK1nUq8uuujXh73gau6NqItnUrBR1JiqAj6Zm8PW8Dz85YSeqRDK46ozF/6tNKrdQIU6GQQnP7ea34ZOFPjBqfzD9/27XYPYBeIsfdmbViOw9PWMbaHamc3bIG9w+Ip1VtjUMUBnU9SaGpUr40fz6vFV+v3cmUJVuDjiNFxKpt+7j69W+5/o0kAF6/NpG3ru+iIlGI1KKQQnV5l0a8+81GHpm0jHPa1KJsKT02VXL3c+oRnpmxkne+2UiF0rHcPyCeq85oTOk4fb8tbPobl0IVFxvDiIHxbP75IK/OWRt0HIlCaRmZjJ27jl5PzObteRu4vEtDZt91Djd0b6oiERC1KKTQndW8Bv3a1+HF2WsYmtiAupXLBR1JokTWOEQya1JS6d4iaxyidR11MQVN5VkC8ZcL25LhzmOTlwcdRaLA6u37uXbst1w3dj4Zmc6rVyfy9g1dVCSihFoUEoiG1cpzU49m/O2z1Vx1RmMSm1QLOpIEYPeBIzwzYxXvzNtAuVKx3HthW645q4m6mKKMCoUE5ne9mvNB0mYeHJ/MJ7d0IyZGl8uWFOkZmfzz2408NX0lew+mcVmXRtxxXitqnFIm6GiSCxUKCUz50nEMv7ANt733Ax8u2Mz/69ww6EhSCOasTOGhCcms2r6fM5tVZ8TAeN2AGeVUKCRQg06tx9tfb+Dxqcu5oEMdKumh9sXWmpT9PDpxGTOXb6dx9fK8fFUnzo+vrRsviwB1BEqgzIwHBrZjZ+oRnv9sddBxJAL2HEjjoQnJ9H16Dt+s28Xwfm2YdnsP+raroyJRRKhFIYHr0KAy/69TQ8bOXcdlnRvqMZXFRHpGJv+av4mnpq1g98E0LuvckDvOa03NihqHKGrUopCocGff1pSNi+XhicuCjiIF4MtVO+j/3Jfc/58ltKpdkQm3duevgxNUJIootSgkKtSsWIY/9m7JI5OWMWvFdk0XXUSt25HKIxOXMWPZNhpULcdLV3bkgvbqYirqAikUZjYGGAgcAdYA17n77tC6BOBloBKQCXR290NB5JTCdc1ZTfjXtxt5aEIy3ZrX0LX0RcjeQ2n8beYq3vhqPaVjYxh2QWuu79ZUc3kVE0H9S5wOtHf3BGAlMBzAzOKAd4Cb3b0d0AtICyijFLLScTHcPyCetSmpvDR7De4edCQ5joxM591vNnDOmNm89uU6Ljm9PrPu7MXve7VQkShGAmlRuPu0bL/OA4aGfj4fWOTuC0Pb7SzsbBKsc9rUol/7Ojw9YyU/bPqZRy7pQL0qmgsqGn21egejJiSzfOs+OjepyhsDutChQeWgY0kEREPb/npgcujnVoCb2VQz+87Mhh1rJzO70cySzCwpJSWlUIJK4Xj+io7cPyCeeWt3cd5Tn/PW1+vJzFTrIlps2JnKjW8lccVr37DvUDovXNGR9286U0WiGLNINe/NbAZQJ5dV97r7J6Ft7gUSgcHu7mZ2J3AL0Bk4AMwE7nP3meHOlZiY6ElJSQWaX4K3adcB/jJuMV+s2kFi46o8NiSBFrV06WxQ9h1K4/nPVjN27nriYo1bzmnBDd01DlGUmdkCd0883nYR63py9z7h1pvZNcAAoLf/t1ptBj539x2hbSYBHckqGFLCNKxWnreu78LH3/3IQxOTufDZL/hj7xbc1LM5pWKjoTFcMmRkOh8kbeKJaSvYsf8IQzo2YNgFraldqWzQ0aSQBHXV0wXA3UBPdz+QbdVUYJiZlSfriqiewNMBRJQoYWYM6dSAHq1q8uD4pTwxbSUTFm1h9JAETm1YJeh4xd7Xa3YyakIyy7bspVPjqvzjms76ey+BItb1FPakZquBMsDRwep57n5zaN2vyboKyoFJ7n7McYqj1PVUckxP3sb9/1nC9n2HuL5bU+44vxXlS+t2oIK2YWcqj05axtSl26hfpRz39GvDgIS6uh+imDnRrqdACkVBU6EoWfYeSmP05OW8+81GGlYrx2ODE+jWokbQsYqFfYfSeH7WasZ+uZ7YGOP3vZrz2x7NNA5RTKlQSLE3b+1Ohn+8mHU7Urm0UwPu6x9P5fKafTYvMjKd95M28aTGIUoUFQopEQ6lZfDszFW8MmctVcuX5qGL2tGvQ92gYxUpOcchRgyI1zhECaFCISXKkh/3cM/Hi1jy4176tqvNqIva69vwcWgcQlQopMRJz8jktS/X8fT0lZSOi+EvF7blss4N9cGXQ/ZxiLjYrHGI35ytcYiSSIVCSqz1O1K55+NFzFu7izOaVeOxwQk0qVEh6FiB0ziE5KRCISWau/Pe/E08OmkZR9Izuf28Vvyme1PiSuiNetnHIRIbV2XEwHgSGmgcoqRToRABtu09xP3/WcK05G20r1+J0UMSaFev5MxJpHEICUeFQiTE3Zm8ZCsjPlnKzweOcGOPZtzWu2WV7lxnAAALzElEQVSx7pPXOISciMDnehKJFmbGhR3qclbz6jwycRkvzV7DlCVbeWxwB7o2qx50vAKlcQiJBLUopMT5ctUOho9bxKZdB7miayPu6deGSmWL/o16GoeQk6WuJ5EwDhxJ56lpK3l97jpqVSzLQxe357z42kHHyhONQ0heqVCInICFm3Zz90eLWL51H/0T6jJyYDtqViwTdKwTonEIyS+NUYicgFMbVuHTP3Tn5c/X8LfPVjN39Q7u6x/PkI71o/YbucYhpLCpRSESsnr7Pu75aDFJG37m7JY1ePSSDjSsVj7oWP9D4xBSkNT1JJIHmZnOO99sYPTk5WQ63Nm3Ndee1YTYmGBbFznHIYZf2Ib+HTQOIfmjQiGSDz/uPsh94xYza0UKpzWswughCbSuU7HQc2gcQiJJhUIkn9ydTxf+xIPjk9l3KI3f9WrBLec0p0xc5D+kc45DDO3UgLv6ahxCCpYGs0Xyycy46LT6nN2yJqPGL+W5mauYvHgLjw1JoFPjqhE7b85xiNev7axxCAmUWhQiJ2jWiu3c+/Fituw9xDVnNuGuvq2pUKbgvmtpHEIKm7qeRCJg/+F0xkxZzlvzNlCvcjkeuaQ9vVrXytcxNQ4hQYnqQmFmY4CBwBFgDXCdu+82s1LAa0BHsrrF3nL3vx7veCoUUtgWbNjFsA8XsSYllUtOr8/9A+KpVqH0SR1D4xAStGgvFOcDn7l7upmNBnD3u83sCmCQu19mZuWBZKCXu68PdzwVCgnC4fQMXvhsNS/OXkPlcqUYMTCeQafWO6GuIt0PIdHgRAtFIE9xcfdp7p4e+nUe0ODoKqCCmcUB5chqcewNIKLIcZWJi+WO81sz4Y/daVCtPLe99wO/eTOJn3YfPOY+G3amctPbSVz+6jz2Hkzj+StO54Obz1SRkKgW+BiFmY0H/u3u74S6nt4GegPlgdvd/ZXjHUMtCglaRqYzdu46npy2ktgY4+4LWnNl18bEhG7UyzkOccs5Lbihe1ONQ0igAr881sxmAHVyWXWvu38S2uZeIB14N7SuC5AB1AOqAl+Y2Qx3X5vL8W8EbgRo1KhRwb8BkZMQG2P85uxm9G1Xh+EfL+b+T5by6cKfePSSDiRt+Pl/xiGG9W1NLY1DSBESWIvCzK4BbgZ6u/uB0LIXgHnu/nbo99eBKe7+frhjqUUh0cTd+XDBZh6euIw9B9MA6NykKiMGtKNDg5LzGFaJfoG3KMIxswuAu4GeR4tEyEbgXDN7h6yupzOAZwKIKJJnZsaliQ3p2bomL81eQ6fGVXU/hBRpQV31tBooA+wMLZrn7jeb2SnAWCAeMGCsu4853vHUohAROXlR3aJw9xbHWL4fuLSQ44iISBiBXB4rIiJFhwqFiIiEpUIhIiJhqVCIiEhYKhQiIhKWCoWIiISlQiEiImEFPilgQTCzFGBDPg5RA9hRQHGCVFzeB+i9RKPi8j5A7+Woxu5e83gbFYtCkV9mlnQidydGu+LyPkDvJRoVl/cBei8nS11PIiISlgqFiIiEpUKR5bgPRyoiisv7AL2XaFRc3gfovZwUjVGIiEhYalGIiEhYJbZQmNnrZrbdzJYEnSW/zKyhmc0ys2VmttTMbgs6U16ZWVkz+9bMFobey4NBZ8oPM4s1s+/NbELQWfLDzNab2WIz+8HMivTDX8ysipl9aGbLQ/9mzgw6U16YWevQf4+jr71m9qeInKukdj2ZWQ9gP/CWu7cPOk9+mFldoK67f2dmFYEFwMXunhxwtJNmWY+Bq+Du+82sFPAlcJu7zws4Wp6Y2R1AIlDJ3QcEnSevzGw9kOjuRf7eAzN7E/jC3V8zs9JAeXffHXSu/DCzWOBHoKu75+eeslyV2BaFu88BdgWdoyC4+xZ3/y708z5gGVA/2FR541n2h34tFXoVyW8zZtYA6A+8FnQWyWJmlYAewD8A3P1IUS8SIb2BNZEoElCCC0VxZWZNgNOBb4JNkneh7pofgO3AdHcvqu/lGWAYkBl0kALgwDQzW2BmNwYdJh+aASnA2FCX4GtmViHoUAXgMuBfkTq4CkUxEnrm+EfAn9x9b9B58srdM9z9NKAB0MXMilzXoJkNALa7+4KgsxSQbu7eEegH3BLqui2K4oCOwEvufjqQCtwTbKT8CXWfDQI+iNQ5VCiKiVB//kfAu+7+cdB5CkKoS2A2cEHAUfKiGzAo1Lf/HnCumb0TbKS8c/efQn9uB8YBXYJNlGebgc3ZWqkfklU4irJ+wHfuvi1SJ1ChKAZCA8D/AJa5+1NB58kPM6tpZlVCP5cD+gDLg0118tx9uLs3cPcmZHULfObuvw44Vp6YWYXQRRKEumnOB4rk1YLuvhXYZGatQ4t6A0Xuoo8cLieC3U6Q1QwrkczsX0AvoIaZbQYecPd/BJsqz7oBVwGLQ337AH9x90kBZsqrusCboas4YoD33b1IX1paDNQGxmV9HyEO+Ke7Twk2Ur7cCrwb6rJZC1wXcJ48M7PywHnATRE9T0m9PFZERE6Mup5ERCQsFQoREQlLhUJERMJSoRARkbBUKEREJCwVCok6ZuZm9mS23+80s5EFdOw3zGxoQRzrOOe5NDQz6awcy5uY2cHQbJ/JZvZ3M8vzv0Mzu9bMng/9fLOZXR1m2yZmdkVezyUllwqFRKPDwGAzqxF0kOxC93acqBuA37v7ObmsWxOaoiQBiAcuzsd5/o+7/93d3wqzSRPgpAqFmZXYe63kv1QoJBqlk/V4x9tzrsjZIjCz/aE/e5nZ52b2vpmtNLPHzOzK0LMtFptZ82yH6WNmX4S2GxDaP9bMxpjZfDNbZGY3ZTvuLDP7J7A4lzyXh46/xMxGh5aNALoDfzezMcd6k+6eDnwFtMjtPGb261D+H8zs5aMFxMyuC2X/nKybLY9mGWlmd4Z+bmFmMyzruR7fhd7/Y8DZoePdblnP/hgbyv+9mZ0T2vdaM/vAzMaTNRFgXTObE9pviZmdfZz/flLM6NuCRKsXgEVm9vhJ7HMq0Jas6ePXAq+5exfLepDTrcDRh7o0AXoCzYFZZtYCuBrY4+6dzawMMNfMpoW27wK0d/d12U9mZvWA0UAn4GeyPlQvdvdRZnYucKe7H/MhP6G7ansDI3Kex8zaAr8iazK+NDN7EbjSzKYDD4bOuQeYBXyfy+HfBR5z93FmVpasL4X3hDIdLY5/BnD3DmbWJpS/VWj/M4EEd98V2m6quz8SKlblj/WepHhSoZCo5O57zewt4I/AwRPcbb67bwEwszXA0Q/6xUD2LqD33T0TWGVma4E2ZM1flJCttVIZaAkcAb7NWSRCOgOz3T0ldM53yXrWwX+Ok7N5aKoVBz5x98lm1ivHeXqTVQzmh6bOKEfWtOtdc5zz30Cr7AcPzctU393HAbj7odDynDm6A38LbbPczDZkO9Z0dz/6vJb5wOuWNfHkf9z9h5wHkuJNhUKi2TPAd8DYbMvSCXWZWtYnX+ls6w5n+zkz2++Z/O//6znnrXHAgFvdfWr2FaEP8NRj5PvFJ+8JOjpGkVP28xjwprsPz5HnYo7/IKcTzRVuu//L4u5zLGta8f7A22Y25jhjIVLMaIxColboG+37ZA0MH7WerG/aABeR9QS8k3WpmcWE+u2bASuAqcDvQt+aMbNWdvwH2nwD9DSzGqEumcuBz/OQJzczgaFmViuUp5qZNQ6ds5eZVQ9lvTTnjqFnkWwOFRXMrEyom2sfUDHbpnOAK0PbtAIakfV38T9C593u7q+SNUtxUZ+WW06SWhQS7Z4E/pDt91eBT8zsW7I+TI/1bT+cFWR9oNcGbnb3Q2b2GlljF9+FWiop5LgaKSd332Jmw8kaJzBgkrt/koc8uR072czuI2vcIAZIA25x93mWdanw18AWslpcuV0ldRXwspmNCu17KbAISDezhcAbwItkDbgvJquldq27H86li6oXcJeZpZH1nPljXoIrxZNmjxURkbDU9SQiImGpUIiISFgqFCIiEpYKhYiIhKVCISIiYalQiIhIWCoUIiISlgqFiIiE9f8B5TeUUu4zbpIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "aic = 50 * np.log(rss/50) + np.arange(1,9)*2\n", "plt.plot(np.arange(1,8), aic[1:])\n", "plt.xlabel('Number of Predictors')\n", "plt.ylabel('AIC')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VdW5//HPQwJhDCTMEiABwiCDCgG0OCCiUmettlrbqm21vb16W3ulFVvUom2tVvuz99paa7X1WmvVKmC14FCnWlEmhYSZMAcIJswQyPD8/tgbPMaQE0JOdobv+/XKKzn77H3294g5T/Zae61l7o6IiEh1WkQdQEREGj4VCxERiUvFQkRE4lKxEBGRuFQsREQkLhULERGJS8VCRETiUrEQEZG4VCxERCSu5KgD1JUuXbp4ZmZm1DFERBqV+fPnf+zuXePt12SKRWZmJvPmzYs6hohIo2Jm62qyn5qhREQkLhULERGJS8VCRETiUrEQEZG4VCxERCQuFQsREYlLxUJEROJqMuMsRBqiigrnxUUFHCirYFTfNPp1aYeZRR1L5KipWIgkyOpte/jhc4uYt2774W2d2rZkVJ80RvZNY1TfNE7I6ESbVkkRphSpGRULkTpWVl7B799Zw69eW0Hr5Bb88ooTOLF3R+av23746/VlhQAktzCOPy6VkX2C4jGqbxrHdWoT8TsQ+Sxz96gz1ImcnBzXdB8StWVbdjH52UUs3rSTc4d2566Lh9EttfVn9tu+9yALN3xSPD7csIOS0goAenZsHVx5hAXk+ONSaZmk7kVJDDOb7+458fbTlYVIHThYVsFDb6ziN2+uIrV1Sx768kjOG97jiP0Tae1aMWFwdyYM7g5AaXkFyzbvZv66Yuav38GCddt5adFmAFq3bMGIjE7BlUfYhJXerlW9vTcR0JWFyDFbtHEHP3huEcu27ObiE4/jjguH1smH+ead+1mwbkdw9bF+O3mbdlJWEfy+9uvS7nC/R07fNPp3bU+LFuo4l6NX0ysLFQuRWiopLedXr63g92/n07VDCj+9ZDgTj++esPPtP1jO4k07DzddLVi/neK9BwFIbZ38qaarE3p3ol2KGg4kPjVDiSTQ3LXF/PC5ReR/vJcv5fTmtvOH0LFNy4Ses02rJMZkpTMmKx0Ad2dt0b6YjvNi3ly+DYAWBkN6ph7uNB/ZJ42MtDa6bVdqTVcWIkdh74Ey7pu9nD+9t5bjOrbhni8M57TsuOvG1Jud+0tZuH47C8Kmqw/X72DvwXIAunVI+aR49E1j6HGppCTrtt3mTlcWInXsXys/5tbnF7Fpx36uOSWTyecOanBNPR3btGT8oG6MH9QNCG7jXb51d1A8wgLyj9wtALRKbsGIXh0PF4+RfdLo2iElyvjSgOnKQiSOXSWl/OylpTw9dwP9urTjF5ePYHRmetSxaq1wVwkL1n9y227upl0cLA9u2+3bue2nBg0O7N6BJHWcN2nq4BapA68v3cptLyxm2+4DXH96P26eOJDWLZtW001JaTl5BTtj+j528PGeAwC0T0nmpD6dDg8aPLFPJ1JbJ7ZvRuqXmqFEjkHx3oNMezGP6R8WMKh7Bx75ag4n9O4UdayEaN0yiVF90xnV95OO8w3F+5m/vvhw8fiff66kwsEMBvdIZUxmGqOz0hmTmV7loENpenRlIRLD3Xl58RZun5HLzv2l3DhhAN8ZP4BWyc17BPXuklI+2rCTeeuKmbu2mAXrdrC/NOg479u5LaMzg8IxOiudzM5tdddVI6IrC5GjVLi7hKnTc5mdt5XhvTry5DfHMqRnatSxGoQOrVtyanYXTs3uAgQjzvMKdjF3TTEfrC3m9aVbeW7+RgC6dkhhdGYaozPTGZ2ZzpCeqer3aAJ0ZSHNnrvz/IJNTPv7EvaXlvP9swfyzVOzSNZ8TDVWUeGs3raHD9YWM3dNMXPXbmfTjv0AdEgJBgyOyQqKx4iMjk2u36cxUwe3SA1s2rGf255fzFsrtpHTN41fXD6C/l3bRx2rSdi0Y//hK4+5a4pZWbgHCG7ZPSGjY3DlkZXOqL5p6jSPkIqFSDUqKpy/zF3Pz19eRnmF88NJg/jaKZmaXymBivceZN7aoM/jg7WfzHXV4lCneXjlMTorjW4d1GleX1QsRI5gXdFefvi3RczJL2bcgM7cc9kIeqe3jTpWs7PvYBkL1+/ggzVBAVm4/pNO88yw0/zQHVd91WmeMOrgFqmkvMJ5/N01/PKV5bRs0YJ7LhvOl0b31odQRNq2SmbcgC6MG/BJp3nupp3Blcea7by6dCvPxnSaj8lMDzrOs9IZ3EOd5vUtoVcWZjYJeBBIAh5193sqPf8r4MzwYVugm7t3Cp+7Bvhx+Nzd7v6n6s6lKwupzqrC3Ux+bhEL1+/grMHduPvSYfTsqBXpGrKKCmfVtj2HrzzmrimmYGcJAB1aJzOqb3DH1ZisoNNc81zVTuTNUGaWBKwAzgY2AnOBq9x9yRH2vwk4yd2/bmbpwDwgB3BgPjDK3bdXdSyoWEjVSssreOTtfB58bSVtU5K488KhXHzicbqaaKQ2bt93+Mpj7tpiVsV0mp+Y0YnRWUEBGdU3jQ7qNK+RhtAMNQZY5e75YaCngYuBKosFcBVwR/jzucCr7l4cHvsqMAn4SwLzShOTV7CTyc8uYsnmXZw/vCd3XjRUE+U1chlpbclIa8ulJ2UAQaf53MO36xbz8Fv5PPTG6sNTtB+68hidma5/+2OUyGLRC9gQ83gjMLaqHc2sL5AF/LOaY3slIKM0QQfKyvmf11fx8Fur6dS2FQ9/ZSSThvWMOpYkQHq7Vpw7tAfnDu0BBFPIL1y/4/Dtuk/PXc8f/70WgKwu7Q4PFhyTlU6fdHWaH41EFouq/hWO1OZ1JfCcu5cfzbFmdgNwA0CfPn1qk1GamIXrt/OD5xaxsnAPl43sxe0XHE+ntlqvurlol5L8qZHmB8sqyC3YefjKY3beVp6ZF3Sad09N4eR+nQ9/aZqS6iWyWGwEesc8zgAKjrDvlcB/Vjp2fKVj36x8kLs/AjwCQZ9F7aNKY7f/YDn3v7Kcx95dQ/fU1jx+3WjODNd0kOarVXILRvYJ1ur41hn9qahwVhYGI80/WFPMv1cXMePD4GOpR2prTu6Xfrh46HbdT0tkB3cyQQf3WcAmgg7uL7t7XqX9BgGzgSwPw4Qd3POBkeFuCwg6uIuPdD51cDdfc/KL+OHfFrGuaB9fHtuHKZ8frM5NqRF3J//jvczJL+K91UXMyS8+PD37oeJxSv+geDTVZqvIO7jdvczMbiQoBEnAY+6eZ2bTgHnuPjPc9SrgaY+pWu5ebGZ3ERQYgGnVFQppnvYcKOOefyzlyTnr6ZPelqeuH8vn+neJOpY0ImZG/67t6d+1PVeP7Yu7s3pbWDzyi/jXqo+ZHl559OzYmpP7deaU8Mqjd3rzWtNcI7ilUXprxTZue34xBTv38/VxWfz3OQNp20pjTKVuBcVjD+/lFzMnv4j384v4eM9BAI4Li8fJ/YMCkpHWOItH5OMs6puKRfOwc18pd720hOfmb6R/13bce/kJjOqbFnUsaSbcnVWFe5iTHzRZzckvomhvUDx6dWrD2H7pMVcejWMKGRULaXJm523hx9NzKd57kG+f0Y+bJmRrqmuJlHvQYR4Uj6CAFMcUj5P7dQ77PNLJSGuYxUPFQpqMoj0HuGNmHn9ftJkhPVO57/IRDOvVMepYIp9xaIqSoLM8+Nq+rxSAjLQ2h++0OqV/Z3p1ahjTzahYSKPn7sz8qIA7Z+ax50AZ/zUhm2+P709LLUokjcShW3XfW/0xc/KLeX/NJ8Wjd3obTs76pHgcF1HxULGQRq1wdwm3PZ/La0u3ckLvTtx3+QgGdu8QdSyRY1JR4awo3M2c1cHdVu+vKWZHWDz6pLf91DiP+ioeKhbSaJVXOF/47b9ZunkXk88dxHXjsjQdtTRJFRXO8q27D4/zeH9NMTv3B8Wjb+e2wZVH/6CAJGqW5MjHWYjU1p/fX8eHG3bw4JUncvGJmhJMmq4WLYwhPVMZ0jOV68ZlUVHhLNuy+3B/x6y8Lfx1XjBNXmbntp+anqRHx/pdTVBXFtKgbNlZwsQH3uKkPp144utjGuV96yJ1paLCWbpl1+HbdN/PL2JXSRkQTIwY22zVPbV2xUNXFtIo/eTFPErLK7j7kmEqFNLstWhhDD2uI0OP68g3Ts2ivMJZunnX4dt0/75oM3/5YAODe3Rg1vdOT2gWFQtpMF5bspV/5G5h8rmD6Nu5XdRxRBqcpBbGsF4dGdarI988rd/h4rGrpDTh51axkAZh74Eybp+Ry6DuHbjh9H5RxxFpFA4Vj/qgYiENwgOvrqBgZwl/+/JJGkch0gDpt1Iil7tpJ4+/u4arx/ZhVN/0qOOISBVULCRSZeUVTHl+MZ3bp/CDSYOjjiMiR6BiIZF64r11LN60kzsuPJ6ObbRgkUhDpWIhkSnYsZ/7X1nO+EFdOX94z6jjiEg1VCwkMnfMzKPcnbsu1pgKkYZOxUIiMSt3C68u2crNEwc2mkViRJozFQupd7tLSrlzZh5Deqby9VOzoo4jIjWgYiH17v5XVrB1dwk/v2y4xlSINBL6TZV69dGGHfzpvbV87eS+nNi7U9RxRKSGVCyk3hwaU9GtQwq3nDso6jgichQ03YfUm8ffXcuSzbt4+Csj6dBaYypEGhNdWUi92FC8jwdeXcHEId04d2iPqOOIyFFSsZCEc3dun5GLGfxEYypEGiUVC0m4lxdv4Y3l2/j+2QPpVU+L0ItI3VKxkITaVVLKnS/mMaxXKtd+LjPqOCJSSwktFmY2ycyWm9kqM7v1CPt80cyWmFmemT0Vs/3ecNtSM/u1qe2iUbp31jKK9hzg55eOIFljKkQarYTdDWVmScBDwNnARmCumc109yUx+2QDU4Bx7r7dzLqF2z8HjANGhLv+CzgDeDNReaXuzV+3nT+/v57rPpfF8Iz6Wc1LRBIjkX/qjQFWuXu+ux8EngYurrTP9cBD7r4dwN0Lw+0OtAZaASlAS2BrArNKHSstr+C25xfTI7U13z9nYNRxROQYJbJY9AI2xDzeGG6LNRAYaGbvmtkcM5sE4O7vAW8Am8Ov2e6+tPIJzOwGM5tnZvO2bduWkDchtfPoO2tYvnU30y4eRvsUDecRaewSWSyq6mPwSo+TgWxgPHAV8KiZdTKzAcAQIIOgwEwws9M/82Luj7h7jrvndO3atU7DS+2tL9rHg6+v4Nyh3Tn7+O5RxxGROpDIYrER6B3zOAMoqGKfGe5e6u5rgOUExeNSYI6773H3PcA/gJMTmFXqiLvz4xm5JLdowU8uGhZ1HBGpI4ksFnOBbDPLMrNWwJXAzEr7TAfOBDCzLgTNUvnAeuAMM0s2s5YEndufaYaShmfmRwW8vWIbt5wzkB4dW0cdR0TqSMKKhbuXATcCswk+6J9x9zwzm2ZmF4W7zQaKzGwJQR/FZHcvAp4DVgOLgY+Aj9z9xURllbqxc18pd/19CSdkdOSrp2RGHUdE6lBCex7d/WXg5Urbbo/52YHvh1+x+5QD30pkNql798xayvZ9pfzp62NIaqFhMSJNiUZJSZ2Yu7aYv3ywgW+cmsXQ4zSmQqSpUbGQY3awLBhT0atTG743MTvqOCKSALoBXo7ZI2+vZmXhHh6/djRtW+l/KZGm6IhXFmbW28yeNrN3zOy28K6kQ89Nr5940tCt+Xgvv/7nKs4f3pMzB3eLOo6IJEh1zVCPEczFdBPQE3jLzDqHz/VNcC5pBNydH72wmJSkFtxx4fFRxxGRBKquzaCruz8c/nyTmX0FeDu87bXySGxphl5YuIl/ry7irkuG0S1VYypEmrLqikVLM2vt7iUA7v6kmW0hGBvRrl7SSYNVvPcgd7+0lJP6dOLqMX2ijiMiCVZdM9SjwNjYDe7+GnAFkJvIUNLw/fzlpezaX8rPLxtOC42pEGnyjnhl4e6/OsL2hQRrVEgz9d7qIp6dv5H/GN+fwT1So44jIvUg7jgLM6s8rbg0YwfKyvnRC4vpnd6G/5qgMRUizUW1xcLMhhPM0yQCwG/eWE3+x3u5+5LhtGmVFHUcEakn1Y2zOJNgdbuv1l8cachWFe7ht2+u5qITjuOMgVo/RKQ5qe5uqJnAWHdfVV9hpOE6NKaidcsWTL1AYypEmpvqmqGeAm43M80fJTw7fyPvrylmynlD6NohJeo4IlLPjlgI3P1bBLfIPll/caQhKtpzgJ+9vJTRmWl8Kad3/ANEpMmp9qrB3e8GZtVTFmmgfvrSUvYeKONnl2pMhUhzFbeJyd2fqI8g0jD9a+XHPL9wE98+oz/Z3TtEHUdEInLU/RFmNsjMfp+IMNKwlJSW8+Ppi8ns3Jb/PHNA1HFEJELV3To7wsxeMbNcM7vbzLqb2d+A14El9RdRovK//1zF2qJ9/PTS4bRuqTEVIs1ZdVcWvye4I+oLwDZgAZAPDDjSVCDSdKzYupvfvb2ay07qxbgBXaKOIyIRq26cRYq7/zH8ebmZ3QLc6u7liY8lUaqocG57fjHtUpL50flDoo4jIg1AdcWitZmdBBy6/WUPMMLMDMDdFyQ6nETjr/M2MG/ddu69fASd22tMhYhUXyw2Aw/EPN4S89iBCYkKJdEp3F3Cz19eytisdK4YlRF1HBFpIKqbovzM+gwiDcPdf19KSWkFP7tsOOFFpIjI0d86K03Xm8sLmflRAd85sz/9u7aPOo6INCAqFgLA/oPlTJ2RS7+u7fiP8f2jjiMiDUxCi4WZTTKz5Wa2ysxuPcI+XzSzJWaWZ2ZPxWzvE47zWBo+n5nIrM3dg6+vZEPxfn526XBSkjWmQkQ+7Yh9FmY2sroD490NZWZJwEMES7BuBOaa2Ux3XxKzTzYwBRjn7tvNrFvMSzwB/NTdXzWz9kBF3HcjtbJsyy4efSefK0ZlcHK/zlHHEZEGqLq7oe4Pv7cGcoCPCG6jHQG8D5wa57XHAKvcPR/AzJ4GLubTo7+vBx5y9+0A7l4Y7ns8kOzur4bb9xzFe5KjUFHhTHl+MaltWnLbeRpTISJVq26K8jPDO6LWASPdPcfdRwEnATVZEKkXsCHm8cZwW6yBwEAze9fM5pjZpJjtO8zseTNbaGb3hVcqUsf+/MF6Fq7fwY/PH0Jau1ZRxxGRBqomfRaD3X3xoQfungucWIPjqrrv0is9TgaygfHAVcCjZtYp3H4acAswGugHXPuZE5jdYGbzzGzetm3bahBJYm3dVcK9/1jGuAGdufSkynVcROQTNSkWS83sUTMbb2ZnhDPOLq3BcRuB2JVyMoCCKvaZ4e6l7r4GWE5QPDYCC909393LgOnAZ/pQ3P2R8Ionp2tXrQl9tKa9uIQD5RXcfYnGVIhI9WpSLK4D8oDvAt8j6HO4rgbHzQWyzSzLzFoBVxKs6x1rOnAmgJl1IWh+yg+PTTOzQxVgAprptk79c9lWXlq8mf+aMICsLu2ijiMiDVx1HdwAuHuJmT0MvOzuy2v6wu5eZmY3ArOBJOAxd88zs2nAPHefGT53jpktAcqBye5eBBBOXPh6OBfVfIJZcKUO7DtYxtTpeWR3a88Np2tMhYjEF7dYmNlFwH1AKyDLzE4Eprn7RfGOdfeXgZcrbbs95mcHvh9+VT72VYI7r6SO/erVFWzasZ9nv30KrZI1LlNE4qvJJ8UdBLfB7gBw9w+BzARmkgTK3bSTx95dy1VjejM6Mz3qOCLSSNSkWJS5+86EJ5GEK69wbnthMWltW3LrJI2pEJGaq0mxyDWzLwNJZpZtZv8D/DvBuSQB/u+9tSzauJOpFxxPx7Yto44jIo1ITYrFTcBQ4ADBMqs7Ce6MkkZk88793Dd7OacP7MpFJxwXdRwRaWTidnAD57v7j4AfHdpgZlcAzyYsldS5O2fmUe7O3RcP05gKETlqNbmymFLDbdJAvZK3hdl5W/nuWQPp07lt1HFEpBGqbtbZzwPnAb3M7NcxT6UCZYkOJnVjz4Ey7piZx+AeHfjmaVlRxxGRRqq6ZqgCYB5wEcGguEN2AzcnMpTUnftfWc6WXSX875dH0jJJYypEpHaqW4P7I+AjM3vK3UsBzCwN6H1oSnFp2BZt3MGf/r2Wq8f2YVTftKjjiEgjVpM/NV81s1QzSydY0+JxM3sgwbnkGJWVVzDl+cV0bp/CDyYNjjqOiDRyNSkWHd19F3AZ8Hi4psXExMaSY/XHf68lr2AXd144lNTWGlMhIsemJsUi2cx6Al8E/p7gPFIHNu3YzwOvrmDC4G6cN7xH1HFEpAmoSbGYRjA77Cp3n2tm/YCViY0lteXu3D49F3f4yUVDNaZCROpETaYof5aYAXjhmtpfSGQoqb1ZuVt4fVkhPzpvCL3TNaZCROpGTaYof5zPLoeKu389IYmk1naVlHLHzDyO75nKdeMyo44jIk1ITab7iO2naA1cymeXR5UG4IFXVrBtzwF+/7UckjWmQkTqUE2aof4W+9jM/gK8lrBEUiu5m3byxHtr+crYvpzQu1PUcUSkianNn5/ZQJ+6DiK1V1HhTJ2RS1rbVtxyzqCo44hIE1STPovdBH0WFn7fAvwwwbnkKDwzbwML1+/gl1ecoHUqRCQhatIM1aE+gkjtFO89yD2zljEmM50vjOwVdRwRaaKqm3V2sLsvM7ORVTztQLG7r0tcNKmJe2ctY3dJGdMu0ZgKEUmc6q4s/hu4Hrj/CM93NrOP3P2rdR9LamLB+u08PXcD15+WxeAeqVHHEZEmrLpZZ68Pv595pH3M7JVEhJL4ysormDo9l+6pKXx34sCo44hIE1ddM9Rl1R3o7s+7+zl1H0lq4sk568gr2MVDXx5J+5SaDJcREam96j5lLgy/dwM+B/wzfHwm8CbwfOJiSXUKd5dw/ysrOC27iyYKFJF6UV0z1HUAZvZ34Hh33xw+7gk8VD/xpCo/e2kpB8oqNFGgiNSbmgzKyzxUKEJbATWSR+S91UVM/7CAb53Rj35d20cdR0SaiZoUizfNbLaZXWtm1wAvA2/U5MXNbJKZLTezVWZ26xH2+aKZLTGzPDN7qtJzqWa2ycz+tybna+oOllUwdUYuGWlt+M74AVHHEZFmpCaD8m40s0uB08NNv3P3F+IdZ2ZJBM1VZwMbgblmNtPdl8Tskw1MAca5+3Yz61bpZe4C3qrZW2n6Hnt3DasK9/CHa3Jo0yop6jgi0ozUaG4od3/B3W9295uBbWZWkz6LMQQLJuW7+0HgaeDiSvtcDzzk7tvD8xQeesLMRgHdAd2eCxTs2M+Dr61k4pDunDWke9RxRKSZqVGxMLMTzewXZraW4K/9ZTU4rBewIebxxnBbrIHAQDN718zmmNmk8HwtCAYDTq5JvuZg2otLcJw7Ljw+6igi0gxVN85iIHAlcBVQBPwVsOoG6VV+iSq2VV5EKZlgFtvxQAbwjpkNA74CvOzuG6q728fMbgBuAOjTp+lOhPvG8kJm5W1h8rmDtPqdiESiuj6LZcA7wIXuvgrAzG4+itfeCPSOeZzBZxdN2gjMcfdSYI2ZLScoHqcAp5nZd4D2QCsz2+Pun+okd/dHgEcAcnJyPrOaX1NQUlrOnTPz6NelHd88LSvqOCLSTFXXDPUFgunI3zCz35vZWVR9tXAkc4FsM8sys1YEVykzK+0znWCQH2bWhaBZKt/dr3b3Pu6eCdwCPFG5UDQXD7+1mnVF+5h28TBSktWpLSLROGKxCDu1vwQMJhixfTPQ3cx+a2Zxp/lw9zLgRmA2sBR4xt3zzGyamV0U7jYbKDKzJQS3405296JjekdNyLqivfzmzdVcMKInp2Z3iTqOiDRj5l7z1hszSweuAL7k7hMSlqoWcnJyfN68eVHHqDPuznV/nMu8tdt5/b/PoHtq66gjiUgTZGbz3T0n3n5Htayquxe7++8aWqFoimbnbeXN5dv43sRsFQoRiVxt1uCWBNt3sIxpL+YxuEcHrv1cZtRxRERULBqiX7++ioKdJdx9yTCSk/RPJCLR0ydRA7Ny624efSefy0dlkJOZHnUcERFAxaJBcXemzsilXUoyUz4/OOo4IiKHqVg0IDM+LGBOfjGTzx1E5/YpUccRETlMxaKB2FVSyt0vLeWEjI5cNabpTl0iIo2TFm9uIB54ZQVFew/w2LU5JLXQ6nci0rDoyqIByN20kyfeW8tXxvZlREanqOOIiHyGikXEKiqCTu20tq245ZxBUccREamSikXEnpm3gYXrdzDlvCF0bNsy6jgiIlVSsYhQ8d6D3DNrGWMy0/nCyMrrQomINBwqFhG6d9YydpeUMe2SoVS3yJOISNRULCKyYP12np67ga+Py2Rwj9So44iIVEvFIgJl5RVMnZ5L99QUvjtxYNRxRETiUrGIwJNz1pFXsIvbLxhK+xQNdRGRhk/Fop4V7i7h/ldWcFp2F84b3iPqOCIiNaJiUc9+9tJSDpRV8JOL1KktIo2HikU9em91EdM/LOBbZ/SjX9f2UccREakxFYt6crCsgqkzcslIa8N3xg+IOo6IyFFR72o9eezdNawq3MMfrsmhTaukqOOIiBwVXVnUg4Id+3nwtZVMHNKds4Z0jzqOiMhRU7GoB9NeXILj3HHh8VFHERGpFRWLBHtjeSGz8rZw04Rseqe3jTqOiEitqFgkUElpOXfOzKNfl3Z887SsqOOIiNSaOrgT6OG3VrOuaB9PfmMsKcnq1BaRxiuhVxZmNsnMlpvZKjO79Qj7fNHMlphZnpk9FW470czeC7ctMrMvJTJnIqwr2stv3lzNBSN6cmp2l6jjiIgck4RdWZhZEvAQcDawEZhrZjPdfUnMPtnAFGCcu283s27hU/uAr7n7SjM7DphvZrPdfUei8tYld+eOmXm0SmrB1AvUqS0ijV8iryzGAKvcPd/dDwJPAxdX2ud64CF33w7g7oXh9xXuvjL8uQAoBLomMGudmp23lTeXb+N7E7Ppnto66jgiIscskcWiF7Ah5vHGcFusgcBAM3vXzOaY2aTKL2JmY4BWwOqEJa1D+w6WMe3FPAb36MC1n8uMOo6ISJ1IZAd3VbPkeRXnzwbGAxnAO2Y27FBzk5n1BP4PuMbdKz7WQO5OAAAMwElEQVRzArMbgBsA+vTpU3fJj8GvX19Fwc4SHrzqJJKTdLOZiDQNifw02wj0jnmcARRUsc8Mdy919zXAcoLigZmlAi8BP3b3OVWdwN0fcfccd8/p2jX6VqqVW3fz6Dv5XD4qg9GZ6VHHERGpM4ksFnOBbDPLMrNWwJXAzEr7TAfOBDCzLgTNUvnh/i8AT7j7swnMWGfcnakzcmmXksyUzw+OOo6ISJ1KWLFw9zLgRmA2sBR4xt3zzGyamV0U7jYbKDKzJcAbwGR3LwK+CJwOXGtmH4ZfJyYqa12Y8WEBc/KLmXzuIDq3T4k6johInTL3yt0IjVNOTo7PmzcvknPvKillwi/folen1jz/nXEktdCiRiLSOJjZfHfPibefRnDXgQdeWUHR3gM8dm2OCoWINEm6XecY5W7ayRPvreXqsX0YkdEp6jgiIgmhYnEMKiqCTu20tq2YfI46tUWk6VKxOAbPzNvAwvU7mHLeEDq2bRl1HBGRhFGxqKXivQe5Z9YyxmSm84WRlQemi4g0LSoWtXTvrGXsLilj2iVDMVOntog0bSoWtbBg/XaenruBr4/LZHCP1KjjiIgknIrFUSorr2Dq9Fy6p6bw3YkDo44jIlIvVCyO0pNz1pFXsIupFxxP+xQNUxGR5kHF4igU7i7h/ldWcFp2F84f3jPqOCIi9UbF4ij87KWlHCir4CcXqVNbRJoXFYsaem91EdM/LOBbZ/SjX9f2UccREalXKhY1cLCsgqkzcslIa8N3xg+IOo6ISL1TD20NPPbuGlYV7uEP1+TQplVS1HFEROqdriziKNixnwdfW8nEId05a0j3qOOIiERCxSKOaS8uwXHuuPD4qKOIiERGxaIabywvZFbeFm6akE3v9LZRxxERiYyKxRGUlJZz58w8+nVpxzdPy4o6johIpNTBfQQPv7WadUX7ePIbY0lJVqe2iDRvurKowrqivfzmzdVcMKInp2Z3iTqOiEjkVCwqcXfumJlHq6QWTL1AndoiIqBi8Rmz87by5vJtfG9iNt1TW0cdR0SkQVCxiLHvYBnTXsxjcI8OXPu5zKjjiIg0GCoWMX79+ioKdpZw1yXDSE7SfxoRkUP0iRhauXU3j76Tz+WjMhidmR51HBGRBkXFgqBTe+qMXNqlJDPl84OjjiMi0uCoWAAzPixgTn4xk88dROf2KVHHERFpcBJaLMxskpktN7NVZnbrEfb5opktMbM8M3sqZvs1ZrYy/LomURl3lZRy90tLOSGjI1eN6ZOo04iINGoJG8FtZknAQ8DZwEZgrpnNdPclMftkA1OAce6+3cy6hdvTgTuAHMCB+eGx2+s6Z0lpOSP7dOLGCQNIaqHV70REqpLIK4sxwCp3z3f3g8DTwMWV9rkeeOhQEXD3wnD7ucCr7l4cPvcqMCkRIbt1aM0jX8thREanRLy8iEiTkMhi0QvYEPN4Y7gt1kBgoJm9a2ZzzGzSURwrIiL1JJETCVbVpuNVnD8bGA9kAO+Y2bAaHouZ3QDcANCnj/obREQSJZFXFhuB3jGPM4CCKvaZ4e6l7r4GWE5QPGpyLO7+iLvnuHtO165d6zS8iIh8IpHFYi6QbWZZZtYKuBKYWWmf6cCZAGbWhaBZKh+YDZxjZmlmlgacE24TEZEIJKwZyt3LzOxGgg/5JOAxd88zs2nAPHefySdFYQlQDkx29yIAM7uLoOAATHP34kRlFRGR6pn7Z7oCGqWcnByfN29e1DFERBoVM5vv7jnx9tMIbhERiUvFQkRE4moyzVBmtg1Ydwwv0QX4uI7iRKmpvA/Qe2momsp7aSrvA47tvfR197i3kzaZYnGszGxeTdrtGrqm8j5A76Whairvpam8D6if96JmKBERiUvFQkRE4lKx+MQjUQeoI03lfYDeS0PVVN5LU3kfUA/vRX0WIiISl64sREQkrmZfLMzsMTMrNLPcqLMcCzPrbWZvmNnScNXB70adqbbMrLWZfWBmH4Xv5SdRZzoWZpZkZgvN7O9RZzkWZrbWzBab2Ydm1qinSzCzTmb2nJktC39nTok6U22Y2aDw3+PQ1y4z+15CztXcm6HM7HRgD/CEuw+LOk9tmVlPoKe7LzCzDsB84JLYlQkbCzMzoJ277zGzlsC/gO+6+5yIo9WKmX2fYNXHVHe/IOo8tWVma4Ecd2/0YxPM7E/AO+7+aDjRaVt33xF1rmMRrk66CRjr7scy5qxKzf7Kwt3fBhr9JIXuvtndF4Q/7waW0kgXjPLAnvBhy/CrUf5VY2YZwPnAo1FnkYCZpQKnA38AcPeDjb1QhM4CVieiUICKRZNkZpnAScD70SapvbDp5kOgkGCJ3cb6Xv4f8AOgIuogdcCBV8xsfrjwWGPVD9gGPB42Dz5qZu2iDlUHrgT+kqgXV7FoYsysPfA34HvuvivqPLXl7uXufiLBwldjwhUUGxUzuwAodPf5UWepI+PcfSTweeA/wybcxigZGAn81t1PAvYCt0Yb6diETWkXAc8m6hwqFk1I2L7/N+DP7v581HnqQtg88CYwKc6uDdE44KKwrf9pYIKZPRltpNpz94LweyHwAjAm2kS1thHYGHO1+hxB8WjMPg8scPetiTqBikUTEXYK/wFY6u4PRJ3nWJhZVzPrFP7cBpgILIs21dFz9ynunuHumQRNBP90969EHKtWzKxdeOMEYZPNOUCjvIPQ3bcAG8xsULjpLKDR3QhSyVUksAkKErhSXmNhZn8BxgNdzGwjcIe7/yHaVLUyDvgqsDhs6we4zd1fjjBTbfUE/hTe3dECeMbdG/Vtp01Ad+CF4G8SkoGn3H1WtJGOyU3An8Pmm3zguojz1JqZtQXOBr6V0PM091tnRUQkPjVDiYhIXCoWIiISl4qFiIjEpWIhIiJxqViIiEhcKhbSIJmZm9n9MY9vMbM76+i1/2hml9fFa8U5zxXhjKZvVNqeaWb7w1lCl5jZw2ZW699FM7vWzP43/PnbZva1avbNNLMv1/Zc0nypWEhDdQC4zMy6RB0kVjj2o6a+AXzH3c+s4rnV4XQmI4DjgUuO4TyHufvD7v5ENbtkAkdVLMys2Y/HEhULabjKCJaKvLnyE5WvDMxsT/h9vJm9ZWbPmNkKM7vHzK4O18ZYbGb9Y15mopm9E+53QXh8kpndZ2ZzzWyRmX0r5nXfMLOngMVV5LkqfP1cM/tFuO124FTgYTO770hv0t3LgH8DA6o6j5l9Jcz/oZn97lARMbPrwuxvEQzIPJTlTjO7Jfx5gJm9ZsG6IAvC938PcFr4ejdbsHbI42H+hWZ2ZnjstWb2rJm9SDB5YE8zezs8LtfMTovz7ydNjP5ikIbsIWCRmd17FMecAAwhmHY+H3jU3cdYsBjUTcChhWEygTOA/sAbZjYA+Bqw091Hm1kK8K6ZvRLuPwYY5u5rYk9mZscBvwBGAdsJPlgvcfdpZjYBuMXdj7hQUDj69izg9srnMbMhwJcIJvArNbPfAFeb2avAT8Jz7gTeABZW8fJ/Bu5x9xfMrDXBH4e3hpkOFcj/BnD34WY2OMw/MDz+FGCEuxeH+81295+GBavtkd6TNE0qFtJgufsuM3sC+C9gfw0Pm+vumwHMbDVw6MN+MRDbHPSMu1cAK80sHxhMMN/RiJirlo5ANnAQ+KByoQiNBt50923hOf9MsFbC9Dg5+4fTsjgww93/YWbjK53nLIKCMDecZqMNwZTtYyud86/AwNgXD+dx6uXuLwC4e0m4vXKOU4H/CfdZZmbrYl7rVXc/tNbLXOAxCyarnO7uH1Z+IWnaVCykoft/wALg8ZhtZYRNqBZ8+rWKee5AzM8VMY8r+PT/75XnuXHAgJvcfXbsE+GH+N4j5PvMp28NHeqzqCz2PAb8yd2nVMpzCfEXg6pprur2O5zF3d+2YEry84H/M7P74vSNSBOjPgtp0MK/bJ8h6Cw+ZC3BX9wAFxOspHe0rjCzFmE7fj9gOTAb+I/wr2fMbKDFXxTnfeAMM+sSNs9cBbxVizxVeR243My6hXnSzaxveM7xZtY5zHpF5QPDtUw2hoUFM0sJm7x2Ax1idn0buDrcZyDQh+C/xaeE5y10998TzG7c2Kf0lqOkKwtpDO4Hbox5/Htghpl9QPCBeqS/+quznOBDvTvwbXcvMbNHCfoyFoRXLNuodJdSZe6+2cymEPQbGPCyu8+oRZ6qXnuJmf2YoB+hBVAK/Ke7z7HgNuL3gM0EV15V3T31VeB3ZjYtPPYKYBFQZmYfAX8EfkPQCb+Y4IrtWnc/UEVz1XhgspmVEqxZf8Tbc6Vp0qyzIiISl5qhREQkLhULERGJS8VCRETiUrEQEZG4VCxERCQuFQsREYlLxUJEROJSsRARkbj+PwNwyuv/ZflvAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "adjr2 = 1 - (rss/(50 - np.arange(1,9)))/(rss[0]/49)\n", "plt.plot(np.arange(1,8),adjr2[1:])\n", "plt.xlabel('Number of Predictors')\n", "plt.ylabel('Adjusted R^2')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VHXWwPHvSSEQILSEUEIISUgBRIGgoPQqoKKuu2t9d9fCqvvau+vqFgvq2tYVXSu6+tp1XaVLFxENiFISAqH30FsSUs77xwwQEZJhkpk75XyeJ88kNzf3ngvJnFvO73dEVTHGGBO+IpwOwBhjjLMsERhjTJizRGCMMWHOEoExxoQ5SwTGGBPmLBEYY0yYs0RgjDFhzhKBMcaEOUsExhgT5qKcDsAT8fHxmpKS4nQYxhgTVBYuXLhDVRNqWi8oEkFKSgq5ublOh2GMMUFFRNZ5sp7dGjLGmDBnicAYY8KcJQJjjAlzlgiMMSbMWSIwxpgwZ4nAGGPCnM8SgYi8LiLbRWTpcctvEpEVIrJMRJ7w1f6NMcZ4xpdXBOOBc6suEJGBwGigq6p2Bv7uw/2zYPVOXpxV6MtdGGNM0PNZIlDVOcCu4xbfAIxV1VL3Ott9tX+Aacu38cSUfJZv3ufL3RhjTFDz9zOCDKCviCwQkdki0tOXO7tpUEeaNIjmkYnLUVVf7soYY4KWvxNBFNAM6AXcBXwgInKiFUVkjIjkikhuUVGRVztrEhvNrYM7Mm/VTmau8OnFhzHGBC1/J4KNwCfq8i1QCcSfaEVVfVlVc1Q1JyGhxjmTTuqKXu1JjW/IIxPyKKuo9Ho7xhgTqvydCP4DDAIQkQygHrDDlzuMjozgvpHZFBYd5L1v1/tyV8YYE5R8WT76LjAfyBSRjSJyDfA6kOouKX0P+I364eb9kOyW9E5twTNfrmRvcZmvd2eMMUHFl1VDl6lqa1WNVtUkVX1NVQ+r6pWq2kVVu6vqDF/tvyoR4Y+jstl96DDjZq7yxy6NMSZohM3I4i5tm3BJ9yTemLeWDbsOOR2OMcYEjLBJBAB3Ds8kMkIYOznf6VCMMSZghFUiSIyrz+/7pzLhxy0sXHf8WDdjjAlPYZUIAMb0SyUxLoa/fZFng8yMMYYwTASx9aK4c1gmizfs4fMftzgdjjHGOC7sEgHAL7on0blNHI9PyqekrMLpcIwxxlFhmQgiIlzlpJv2FPP6vDVOh2OMMY4Ky0QAcHZaPEM7JTJuZiE7DpQ6HY4xxjgmbBMBwH0jsigpq+CZaQVOh2KMMY4J60SQmtCIK3u1591v11Owbb/T4RhjjCPCOhEA3DK4I41ionhkQp7ToRhjjCPCPhE0a1iPmwd3ZHZBEbMLvOt7YIwxwSzsEwHA//ROoX2LWB6ZsJxy61lgjAkzlgiAelER3Dcii4JtB/ggd6PT4RhjjF9ZInAb3rkVZ6Y05+lpK9hfYj0LjDHhwxKBm4jwwHnZ7DhwmJdmFzodjjHG+I0lgiq6JjXlom5teWXuGjbutp4FxpjwYIngOHcNz0SAJ6escDoUY4zxC0sEx2nTtAFj+qXy2eLNLN6wx+lwjDHG53zZvP51EdnublR//PfuFBEVkXhf7b82ft8/jfhGMTz8xXLrWWCMCXm+vCIYD5x7/EIRaQcMBdb7cN+10igmijuHZZC7bjeTlm51OhxjjPEpnyUCVZ0DnKgf5DPA3UBAn2r/MqcdWa0a89ikPErLrWeBMSZ0+fUZgYhcAGxS1R/8uV9vREYID4zqxIZdxbz59VqnwzHGGJ/xWyIQkVjgj8CDHq4/RkRyRSS3qMiZOYD6dIxnYGYCz89Yxa6Dhx2JwRhjfM2fVwRpQAfgBxFZCyQBi0Sk1YlWVtWXVTVHVXMSEhL8GOZP3T8ym0OHK3juS+tZYIwJTX5LBKq6RFVbqmqKqqYAG4HuqhrQT2M7Jjbm8jOTeXvBelZtP+B0OMYYU+d8WT76LjAfyBSRjSJyja/25Wu3DulIbHQkj020ngXGmNDjy6qhy1S1tapGq2qSqr523PdTVHWHr/Zfl1o0iuEPg9KZnr+deauCImRjjPGYjSz20G/PTiGpWQMenpBHRWVAV74aY8wpsUTgofrRkdw7Iou8Lfv4eKH1LDDGhA5LBKdg1Gmt6Z7clCenruBgabnT4RhjTJ2wRHAKXD0LOlG0v5R/zVntdDjGGFMnLBGcou7JzTj/9Da8PKeQLXuLnQ7HGGNqzRKBF+4enkmlWs8CY0xosETghXbNY7n6nA58smgTSzbudTocY4ypFUsEXrpxYBotGtbj4QnWs8AYE9wsEXgprn40tw3NYMGaXUxdvs3pcIwxxmuWCGrh0p7t6NiyEY9NzONweaXT4RhjjFcsEdRCVGQE94/KZu3OQ7z9zTqnwzHGGK9YIqilARkJ9O0Yz3PTV7LnkPUsMMYEH0sEtSQi/HFUNvtLyvjH9FVOh2OMMafMEkEdyGoVx697JvPW/LWs2XHQ6XCMMeaUWCKoI7cPzSAmKoKxk6xngTEmuFgiqCMJjWO4cWA6U5Zt45vVO50OxxhjPGaJoA5d06cDbZrU5+EJy6m0ngXGmCBhiaAO1Y+O5J4RWSzdtI9Pv9/kdDjGGOMRSwR17PyubTg9qQlPTllB8eEKp8Mxxpga+bJ5/esisl1EllZZ9qSI5IvIjyLyqYg09dX+nRIR4epZsHVfCa/MtZ4FxpjA58srgvHAucctmwZ0UdWuQAFwnw/375ieKc0ZeVorXpxVyLZ9JU6HY4wx1fJZIlDVOcCu45ZNVdUjPR6/AZJ8tX+n3XNuFuWVlTw11XoWGGMCm5PPCK4GJjm4f59q36Ihvz07hQ8XbmT55n1Oh2OMMSflSCIQkT8C5cA71awzRkRyRSS3qKjIf8HVof8d1JGmDaKtZ4ExJqD5PRGIyG+A84ArtJp3R1V9WVVzVDUnISHBfwHWoSYNorl1SAZfF+5kRv52p8MxxpgT8msiEJFzgXuAC1T1kD/37ZTLz0omNaEhj07Mo6zCehYYYwKPL8tH3wXmA5kislFErgH+CTQGponIYhF5yVf7DxTRkRHcPyKbwqKDvPvteqfDMcaYn4ny1YZV9bITLH7NV/sLZIOzW3J2WguemVbA6DPa0qRBtNMhGWPMUTay2A+O9CzYU1zGCzOtZ4ExJrBYIvCTzm2acEn3JMbPW8v6nWHxeMQYEyQsEfjRncMziYwQHp+c73QoxhhzlCUCP0qMq8/1/dOYsGQLuWt31fwDxhjjB5YI/Oy6fh1oFVefv03Is54FxpiAYInAz2LrRXHn8Ex+2LCHz3/c7HQ4xhhjicAJF3drS5e2cTwxeQUlZdazwBjjLEsEDoiIEP44shOb9hTz2ldrnA7HGBPmLBE4pHdaC4Z2SmTczFUU7S91OhxjTBirMRGISEMRiajydYSIxPo2rPBw34gsSssreebLAqdDMcaEMU+uCKYDVd/4Y4EvfRNOeElNaMRVvdvz3rfrWbF1v9PhGGPClCeJoL6qHjjyhftzuyKoI7cM7kjj+tE8MjHP6VCMMWHKk0RwUES6H/lCRHoAxb4LKbw0ja3HTYPSmVNQxKwV1rPAGON/niSCW4EPRWSuiMwF3gf+17dhhZf/6Z1CSotYHp2YR7n1LDDG+FmNiUBVvwOygBuAG4FsVV3o68DCSb2oCO4dkU3BtgO8n7vB6XCMMWHmpIlARAa5Xy8GzgcygI7A+e5lpg4N75zImR2a8/TUAvaXlDkdjjEmjFR3RdDf/Xr+CT7O83FcYUdEeGBUNjsPHmbcrEKnwzHGhJGTdihT1Yfcn/5VVX8y/FVEOvg0qjDVNakpF3dry2tfreGKs5JJambFWcYY3/PkYfHHJ1j2UV0HYlzuHJ5JhMATk1c4HYoxJkxU94wgS0R+ATQRkYurfPwWqF/ThkXkdRHZLiJLqyxrLiLTRGSl+7VZnRxFCGnTtAFj+qby3x828/363U6HY4wJA9VdEWTiehbQlJ8+H+gOXOfBtscD5x637F5guqp2xDVi+d5TjDcs/L5/GgmNY3h4Qh6q1rPAGONb1T0j+Az4TER6q+r8U92wqs4RkZTjFo8GBrg/fxOYBdxzqtsOdQ1jorhzWAb3fLyEiUu2Mqpra6dDMsaEME+eEVwkInEiEi0i00Vkh4hc6eX+ElV1C4D7taWX2wl5l/RoR1arxoydnGc9C4wxPuVJIhimqvtw3SbaiGs8wV0+jQoQkTEikisiuUVFRb7eXcCJjBAeGNWJDbuKefPrtU6HY4wJYZ4kgmj360jgXVWtTdf1bSLSGsD9etLJdVT1ZVXNUdWchISEWuwyePXpGM+grJb8c8Yqdh6wngXGGN/wJBF8LiL5QA4wXUQSgBIv9/df4Dfuz38DfObldsLG/SOzOFRWwXPTVzodijEmRHky19C9QG8gR1XLgIO4HvpWS0TeBeYDmSKyUUSuAcYCQ0VkJTDU/bWpRnrLxlxxVjLvLFjPqu3Ws8AYU/dOWjUkIoNUdUbVeYVEpOoqn1S3YVW97CTfGnxKERpuGdyRT7/fxF+/yOPN3/U8/v/BGGNqxeYaCgItGsVw+9AM5hQUMWnpVqfDMcaEGJtrKEhc1as9H+Zu5K+fL6dfRgKNYk76X2eMMafE5hoKElGRETxyURe27S/h2WnW7N4YU3eqe0aQBXTGPddQlW/F4cFcQ6budUtuxqU9k3nj67X8okcS2a3jnA7JGBMCfDnXkPGBe87NpEmDaB74z1IqK20eImNM7flsriHjG01j63HfiCzu+uhHPlq4kV/1bOd0SMaYIOfJE8fvReQPuG4THb0lpKpX+ywqU61LeiTxYe5GHpuUx9BOiTRrWM/pkIwxQcyTh8X/BloBw4HZQBJgI5scJCL87cIu7C8p5/HJ+U6HY4wJcp4kgnRV/RNwUFXfBEYBp/k2LFOTzFaNuaZPB977bgML19Vm+idjTLjzJBGUuV/3iEgXoAmQ4rOIjMduHtyRNk3q88dPl1JeUel0OMaYIOVJInjZ3VLyAVyTxi0HHvdpVMYjDWOiePD8zuRv3c94m6raGOMlTxLBdFXdrapzVDVVVVsCU30dmPHM8M6JDMxM4JlpBWzd6+2ksMaYcGYji4OciPCXC7pQXqn87YvlTodjjAlCNrI4BCS3iOV/B6bz1LQCflVQRP+M8GzkY4zxjo0sDhFj+qeSGt+Qhz5baj2OjTGnxEYWh4iYqEj+dmEXrnh1AS/NLuTWIRlOh2SMCRKePCO4SETiRCRaRKaLyA4RudLnkZlTdk56PBec3oZxswpZu+Og0+EYY4KEJ4lgmKruw3WbaCOQAdzl06iM1x4YlU1MZAQP/ncZqjYpnTGmZp4kgmj360jgXVW1YawBrGVcfe4Y5upmNnGJdTMzxtTMk0TwuYjkAznAdBFJAGpVsC4it4nIMhFZKiLviohVIdWhK3u1p3ObOP76xTIOlJY7HY4xJsDVmAhU9V6gN5CjqmXAIWC0tzsUkbbAze7tdQEigUu93Z75uajICB6+sAvb95fyjHUzM8bUwJMrAtwjiyvcnx9U1drec4gCGohIFBALbK7l9sxxuiU347Izkxn/9VqWb97ndDjGGC/4q/mUR4mgLqnqJuDvwHpgC7BXVW3KCh+4Z3gWTRtE88B/llg3M2OCRXkpunoWhf93B6sf7sb65d/6fJd+TwTuCexGAx2ANkDDE5WjisgYEckVkdyioiJ/hxkSmsRGc9/IbBat38OHCzc4HY4x5kRUoagAvnkR3vklFWPbI2+Npt2KNzgojSgt9n0puCcdynBPMdEHUOArVf20FvscAqxR1SL3tj8BzgberrqSqr4MvAyQk5Njp7Ne+kX3tnzw3QYem5TP0E6taG7dzIxxXvFuWD0bCqdD4UzY6zpR2xqdxOSSvvxYvwd9Bl/E6LMyiIwQn4dTYyIQkXFAOvCue9HvRWSIqv7By32uB3qJSCxQDAwGcr3clqmBiPDwRV0Y+dxcHp+Uz+OXdHU6JGPCT0U5bFoIhTNcb/6bFoJWQkwTSpP7MDXuUp4qTKKoshXXD0zj4b4diK3n0Xl6nfBkT/2BLuoenSQibwJLvN2hqi4QkY+ARUA58D3uM3/jGxmJjbmmbwf+NXs1v8xJIieludMhGRP69qyHVdNdb/5rZkPJXpAIaNMd+t1Fafv+vLa2BeNmr6O4rIJLe7bj1iEZJDSO8XuoniSCFUAysM79dTvgx9rsVFUfAh6qzTbMqbl5UEc+X7yZB/6zlC9u6kNUpN8fDxkT2koPwNqvjp3171zlWh7XFrIvgPTB0KE/FfWb8en3m3jq/RVs2VvIkOxE7h2RRXrLRo6F7kkiaAHkiciRR9c9gfki8l8AVb3AV8GZutMwJoqHLujM7/+9kPFfr+XavqlOh2RMcKushG1Ljp31r/8GKssgqgGk9IGca1xv/vEZIK77/HNXFvHoxK/I27KP05Oa8Oyvz+Cs1BYOH4hnieBBn0dh/GJYp0QGZbXkmWkFjOramtZNGjgdkjHBZf82WD3T9ea/eiYcdFc0JnaBXje43vjb9YLon06WkLdlH49NymdOQRFJzRrwj8u6cd5prYnww4NgT1TXmCYdSFTV2cct7wtsVtVCXwdn6parm1lnhjw9m799sZxxV/RwOiRjAlt5Kayf7zrjXzXDdQUAEBsPaQMhbbDrtXGrE/741r0lPDV1BR8t2khc/WgeGJXNVb3bExMV6ceDqFl1VwTPAvefYHmx+3vn+yQi41Ptmsdy06B0/j61gFkrtjMgs6XTIRkTOFRhx0p3WecM1z3/skMQEQ3JvWDwQ5A2CFp1hYiTP2c7UFrOv2YX8src1VRWwrV9OvCHgek0jQ3M8u3qEkGKqv7sobCq5opIis8iMj53Xb9UPvl+Ew/9dxlTbm1B/ejAOjsxxq9OUtNPi3TodqXrrD+lD8TU/DC3rKKS977bwHNfFrDjwGHOP70Ndw/PpF3zWB8fRO1UlwiqmxHUbi4HsZioSB4e3YXLX13Ai7MKuW2odTMzYeRoTb/7rL9KTT+p/aDvHa6z/mbtPd6kqjJt+TbGTs5nddFBzuzQnNd+k83p7Zr68EDqTnWJ4DsRuU5VX6m6UESuARb6Nizja2enxzP6jDa8OKuQC7u1pUN8Q6dDMsZ3jtb0T4fVc6D0pzX9pA2Gtj0g8tQHcS3esIdHJ+Tx7dpdpCU05JX/yWFIdktEAuNBsCfkZF2sRCQR+BQ4zLE3/hygHnBRHcxA6rGcnBzNzbXBx3Vt+/4SBv99NmckN+Wtq88Mql9cY6p10pr+JEgf5Drj79AfYr0fXLl+5yGemJLPFz9uIb5RPW4dksGlPdsF1BgdEVmoqjk1rVdd8/ptwNkiMhDo4l48QVVn1FGMxmEtG9fnzuGZPPTfZUxYsoXzurZxOiRjvFNTTX/Pa11v/lVq+r2159Bhnp+xirfmryUyQrh5UDpj+qfRKMZ/U0LUtRojV9WZwEw/xGIccGWv9ny4cAN//Xw5/TMSaFw/uuYfMiYQnLSm/7Rqa/q9VVJWwb/nr+P5GSs5UFrOL3u04/ZhGSTGBX+DxeBNYaZOREYID194GheNm8cz01by4PmdnA7JmBOrtqbffbunmpp+b1VWKp//uJknp6xg4+5i+mckcN/ILLJaxdXpfpxkicBwRrumXH5mMuO/XsMverSlc5smTodkjGc1/emDXVcA1dT018Y3q3fy6MQ8fty4l06t43j7mq706Rjvk305yRKBAeDu4VlMXrqVP/1nKR9df3bADH03YaYOa/prY9X2/YydlM+Xedtp3aQ+T/3ydC7q1jZk/y4sERjA1c3s/pHZ3PHhD3yQu4FLz0x2OiQTDnxQ018b2/eX8OyXK3n/uw3ERkdy97mZXH1Oh5AfdGmJwBx1cfe2vJ+7gbGT8xnW2bqZGR/Zve5YWWcd1/R769Dhcl6Zs4Z/zSnkcHklV/Vqz02D0mnRyP+9AZxgicAcJSI8fKGrm9nYSXk8ccnpTodkQkF1Nf2dR9dJTb+3KiqVD3M38PS0ArbvL+Xczq24Z0RW2A2wtERgfiIjsTHX9k3lpdmF/CqnnXUzM6fOjzX93lJVZhUUMXZiPiu27ad7clPGXdE9bH/fLRGYn7l5cDqf/7CZP366lC9u7kN0AI2UNAHKzzX9tbF0014em5THvFU7ad8ilnFXdGdEl1ZhPbLeEoH5mdh6UTx0fifG/Hsh4+et5bp+1s3MHOdITf8qd3WPn2r6a2PTnmKemrKCTxdvommDaB46vxNXnNWeelF2ouNIIhCRpsCruKauUOBqVZ3vRCzmxIZ2SmRwVkue+bKA8063bmZhTxV2FLgHc0133fMvL/ZrTb+39pWUMW5mIa/PWwPA7/ulccOANJo0sFH0Rzh1RfAcMFlVLxGRekBgT9YdhkSEP1/QmaHPzOavny/nxSutm1nYKd4Nq2cdG8m7b6NreYt06P4/rrN+P9T0e+tweSXvLFjHP6avZPehMi7u1pY7hmfStqmd1BzP74lAROKAfsBvAVT1MK4ZTk2AcXUz68iTU1Ywc8V2Blo3s9BWU01/vzv9WtPvLVVl0tKtPDE5n7U7D3F2WgvuH5lNl7Y2Yv5knLgiSAWKgDdE5HRcU1zfoqoHHYjF1OC6vql8smgjD322jN63WTezkHOymv62PRyr6a+Nhet28ejEfBau201GYiPe+G1PBmQmhPWDYE+ctB+Bz3YokgN8A5yjqgtE5Dlgn6r+6bj1xgBjAJKTk3usW7fOr3GaY75etYPLX13AzYM7crt1Mwtufpin3wlrdhzk8Un5TF62lZaNY7h9aAaX9EgKqN4ATqh1PwIf2ghsVNUF7q8/Au49fiVVfRl4GVyNafwXnjne2enxXHhGG16aVciFZ7QhNSEw7wmbE/Copn8wxHd0rKa/NnYeKOUf01fyzoL11IuK4PahGVzbtwOx9YLjCiZQ+P1fS1W3isgGEclU1RXAYGC5v+Mwp+b+UdlMz9/Og58t49/XWDezgFZdTX/vG11n/cm9ISp4p08oPlzB6/PW8OKsQorLKri0ZztuHZJBQuPgPSYnOZU2bwLecVcMrQZ+51AcxkMtG9fnruGZPPjZMr74cQvnn27dzAJGTTX96YMhdSA0TnQ2zjpQUal8smgjT00tYOu+EoZ2SuSec7NIb2lXqbXhSCJQ1cW4+h+bIHLFWe35MHcjf/tiOQMyrZuZY4K4pr82ZhcU8djEPPK37uf0dk157tIzOCu1hdNhhQS7kWY8FhkhPHJRF0a/MI+npxXw0PmdnQ4pfAR5TX9tLN+8j8cm5TF35Q7aNW/A85d147yure32ZB2yRGBOSdekplx5Vnve/Hotl/RIsm5mvhIiNf21sXlPMX+fuoJPv99EkwbR/Om8TlzZK5mYKCthrmuWCMwpu3N4JpOWbuGB/yzlY+tmVndCrKbfW/tKynhxViGvf7UGBcb0TeXGgek2JYQPhfZvlPGJJg2i+eOobG57/wfez93AZdbNzDtHa/rdZ/0BNE+/E46fEuKibm25Y1gGSc1sBhpfs0RgvHLhGW15/7sNjJ2Uz7BOiWHTyalWKith64/us/7Qq+n3lk0J4TxLBMYrR7qZnfvsXMZOyufJX1o3sxPav+3YG3+I1vTXRu7aXTw6MY9F6/eQmdiYN37XkwEZNiWEv1kiMF5Lb9mY6/ql8uKsQn7Vsx09w7S700+UlcCGb46N5N221LU8BGv6a2N10QEen5zPlGXbaNk4hsd/cRqX9GhHpD1vcoTf5xryRk5Ojubm5jodhjmB4sMVDHl6No1iosKzm1lNNf1H3vxDrKbfWzsOlPLclyv5v2/XUz8qguv7p3GNTQnhM4E815AJIQ3qRfLnCzpz3Vu5vDFvDWP6pTkdku8dqek/MpI3jGr6vVV8uILXvlrNS7NXU1xWwWVntuOWwTYlRKCwRGBqbWinRIZkJ/Lslys5r2sb2oRa44+qNf2rpsPmRWFX0++tikrl40Ubedo9JcSwToncMyKLNJu4MKBYIjB14qHzOx3tZvbSVSHQzWz3umNlnT+r6b/b9cYfBjX93lJVZhcUMXZSPvlb93NGu6b847JunNnBniMFIvstNnWiXfNYbh7ckScmr2Bm/nYGZgVZN7Maa/oHQ4d+YVPTXxvLNu/lsYn5fLVqB8nNY/nn5d0YdZpNCRHILBGYOnNtn1Q+WbSJh/67jN5pAd7N7GQ1/dGxYV3TXxub9hTz1JQVfLrYNSXEg+d14spe7akXZQ/JA50lAlNn6kVF8LfRXbjslW8YN3MVtw/LdDqkn7Kafp/YV1LGuJmFvD5vDQBj+qVy4wCbEiKYWCIwdap3Wgsu6taWl2av5sJubZ3tZmY1/T51uLySt79Zx/MzVrKnuIyLzmjLHcMzaRtqxQJhwBKBqXP3j8zmy7xt/u9mFqbz9PubqjJxyVaemJLPup2HOCe9BfeNsCkhgpklAlPnEhrHcPfwTP702TI+/3ELF/iym1lNNf3pg6H9OVbTX0e+W7uLRybksXiDa0qI8b/rSX+bEiLoWSIwPnH5We35cOGxbmZxddXNzGr6HVFYdIDHJ+Uzdfk2EuNieOIXXflFjySbEiJEWCIwPhEZ4ZqUbvQL83h6agF/vqAW3cyspt8xx08JceewDK7pk0qDegFcEWZOmWN/OSISCeQCm1T1PKfiML7TNakpV/Vqz1vzXd3MPL6HbDX9jis+XMGrc1fz0uxCSsorufzMZG4Z0pF4m248JDl5CnULkAfEORiD8bE7hmUycclWHvjPUj654STdzI7W9Lvv81tNv2MqKpWPF27kqWkr2LavlOGdE7n7XJsSItQ5kghEJAkYBTwC3O5EDMY/mjSI5oFR2dz6/mLe+24Dl5/l7mZ2tKbf/eZ/aIdrudX0O0JVmVVQxNiJ+azYtp9uyU355+XdbWrxMOHUFcGzwN1AY4f2b/xo9Blt+OTbQmZN+oALd+widv1sq+kPIEs37eWxSXnMW7WT9i1iGXdFd0Z0aWWVQGHE74lARM4DtqvqQhEZUM16Y4AxAMnJ1hM36Byp6V81HSmcwfiiuURQQnluFLTvbTX9AWDj7kM8NbWAT7/fRLNM1nocAAAQYklEQVTYaB46vxNXnGVTQoQjJ64IzgEuEJGRQH0gTkTeVtUrq66kqi8DL4OrMY3/wzSn7NAuWDP7hDX9Ed1/w0d7M3jwh2aM7z/QZqF00N5DZYybtYo3vl6LADcMSOOGAWl1V+Jrgo6jHcrcVwR31lQ1ZB3KAlRFOWzKPTaS9yc1/f1dt3yq1PQf6WbWMCaSCTf3Db9uZg4rLa/g3/PX8c+Zq9hbXMbF3ZK4Y1hG6PWPMEdZhzLjGzXV9KcPhjbdT1jT36BeJH+5oDPXvpXLiOfmcn3/NEaf0cYSgo9VVipfLNnCk1Py2bCrmH4ZCdx7bhad2ljBnnGxnsWmelVr+ldNh12FruVxSZA+yFXWmdofGjTzeJMTl2zhH9NXkr91P22bNmBMv1R+3bNdYE9bHaS+Wb2Txybm8cPGvWS3juP+kVn07ZjgdFjGTzy9IrBEYH6qppr+tMGu2z21rOlXVWau2M4LMwtZuG438Y3q8btzOnBV7/Z2r7oOrNy2n7GT8pmev502Tepzx7BMLurW9sTjOEzIskRgPFddTf+Rs/7kXj6p6VdVvl2zixdmFTKnoIjGMVFc1bs9V/fpYKNYvbB9XwnPfFnA+99toGG9KG4cmM7vzkmxq60wZYnAnFxZCayff6xJy5Ga/oYJxx7wOlDTv3TTXl6cVcjEpVuoFxnBpT3bcV2/VJKaxfo1jmB0oLScl+es5pU5qymvrOTKXu25aVBHmjes53RoxkGWCMwxVWr6KZzx83n60923ewKkpr+w6AD/ml3IJ4s2ATD6jLbcMCCV9JY2/vB45RWVvPfdBp79ciU7DpQyqmtr7h6eSfsWDZ0OzQQASwTh7qQ1/R2PjeQN8Hn6N+8p5pW5q3n32/WUllcyvFMrbhyYRtekpk6H5jhVZdrybYydnM/qooOcmdKc+0Zm0S3Z84f2JvRZIgg3p1jTH0x2Hihl/NdrGf/1WvaXlNO3Yzw3DEijd2qLsJwG4fv1u3lsYj7frt1FWkJD7h2RzZDslmH5b2GqZ4kgHFRX0582uNqa/mC0v6SMdxas59W5a9hxoJRuyU25cUA6g7NahkU1zLqdB3li8gomLNlCfKMYbhvakV/ntCPKxmGYk7BEEIpKD8DaucfO+uugpj8YlZRV8OHCjfxrdiEbdxeTmdiYGwemMeq01iH5prjr4GH+MX0l7yxYR1REBGP6pTKmXyoNY0IjwRvfsUQQCvxU0x+syioq+eLHzYybWcjK7QdIbh7LmH6pXNIjKSTKJUvKKnh93hpenFnIwcPl/LpnMrcN6UjLuPpOh2aChCWCYLV/q+tN34Ga/mBVWal8mbeNF2YV8sOGPSQ0juHaPh24old7GgXhWXNFpfLp95t4auoKtuwtYUh2S+4dkWVVU+aUWSIIFgFa0x+MVJX5hTt5YdYq5q3aSZMG0fymd3t+e06HoKmnn1NQxGOT8snbso/Tk5pw38hseqW2cDosE6QsEQSqn9T0T4e18wK6pj9Y/bBhD+NmrWLKsm00iI7ksjOTua5fB1o3CcyZNpdv3sdjk/KYu3IH7Zo34O7hWZzXtbVVAplasUQQSEKgpj9Yrdy2nxdnF/LZ4s1ECFzcLYnrB6TRIT4wBlxt3lPM36eu4NPvN9GkQTQ3DerIlb2SiYkK/mccxnmWCJx0pKb/yEje42v6j7RlDMKa/mC1YdchXpm7mve+20B5RSUjTmvNjQPS6NymiSPx7CspY9zMQt6YtwYFfndOCjcOSKdJA5twz9QdSwT+tnvtsbLONXOgdJ+7pj/n2Fl/CNX0B6ui/aW8Pm8N/56/jgOl5QzITODGAel+65h2uLySt79Zx/MzVrKnuIyLzmjLHcMzaWvNYYwPWCLwtZPV9DdpV+Uhb+jX9AervcVlvP3NOl77ag27Dh6mZ0ozbhyYzoCMBJ/cl1dVJizZwhOTV7B+1yH6pMdz74gsurR15orEhAdLBHXNavpDUvHhCt7/bj0vz1nN5r0ldGodxw0D0hh5Wmsi62i08rdrdvHIxDx+2LCHrFaNuW9kNv06xtuDYONzlgjqwslq+lud5j7rt5r+UHG4vJLPFm/ixdmFrC46SIf4hlzfP5WLuiVRL8q76q1V2w8wdlI+X+Zto1Vcfe4YlsHF3ZPqLMEYUxNLBN6wmv6wV1GpTF22lRdmrWLppn20iqvPtX07cPlZycTW8+z5zvb9JTz35Ure+24DDaIjuWFAGlef04EG9awSyPhXwCYCEWkHvAW0AiqBl1X1uep+xmeJwKOa/sGQ2MVq+sOMqjJ35Q5emLmKBWt20Sw2mt+e3YHfnN2eprEnHpx2sLScV+au5uU5qzlcfqQ5TDotrNOacUggJ4LWQGtVXSQijYGFwIWquvxkP1OniaC6mv4jg7lS+kC9wKgzN85buG4X42YWMj1/Ow3rRXJFr/Zc26fD0Tl/yisq+SB3I898WUDR/lJGntaKu4dnkRIgYxVM+ArYRPCzAEQ+A/6pqtNOtk6tEoEnNf1pg6BpspdHYMJF/tZ9vDirkM9/2ExURASX5CRxVofmPD9jFau2HyCnfTPuH5VNd2sOYwJEUCQCEUkB5gBdVHXfydbzOhHMfgK+ft5q+k2dWrfzIP+as5qPcjdyuKKS1PiG3DMii2GdEq0SyASUgE8EItIImA08oqqfnOD7Y4AxAMnJyT3WrVt36jv5/m3Y8K3V9Buf2LavhOVb9tEnPZ7oEOyDYIJfQCcCEYkGvgCmqOrTNa0fEOMIjDEmyHiaCPx+GiOua+fXgDxPkoAxxhjfcuJ69hzgKmCQiCx2f4x0IA5jjDGA35+WqupXgD1RM8aYAGFPuIwxJsxZIjDGmDBnicAYY8KcJQJjjAlzlgiMMSbMOT7XkCdEpAjwYmgxAPHAjjoMx0l2LIEnVI4D7FgCVW2Opb2qJtS0UlAkgtoQkVxPRtYFAzuWwBMqxwF2LIHKH8dit4aMMSbMWSIwxpgwFw6J4GWnA6hDdiyBJ1SOA+xYApXPjyXknxEYY4ypXjhcERhjjKlGyCYCEWknIjNFJE9ElonILU7H5C0RqS8i34rID+5j+YvTMdWGiESKyPci8oXTsdSGiKwVkSXuGXSDumGGiDQVkY9EJN/9N9Pb6ZhOlYhkVpnReLGI7BORW52Oy1sicpv7732piLwrIvV9tq9QvTUkIq2B1qq6SEQaAwuBC1V1ucOhnTJ3D4eGqnrA3dTnK+AWVf3G4dC8IiK3AzlAnKqe53Q83hKRtUCOqgZ9vbqIvAnMVdVXRaQeEKuqe5yOy1siEglsAs5SVW/HIDlGRNri+jvvpKrFIvIBMFFVx/tifyF7RaCqW1R1kfvz/UAe0NbZqLyjLgfcX0a7P4Iyg4tIEjAKeNXpWIyLiMQB/XA1jEJVDwdzEnAbDBQGYxKoIgpoICJRQCyw2Vc7CtlEUJWIpADdgAXORuI99+2UxcB2YJqqBuuxPAvcDVQ6HUgdUGCqiCx099gOVqlAEfCG+5bdqyLS0OmgaulS4F2ng/CWqm4C/g6sB7YAe1V1qq/2F/KJQEQaAR8Dt6rqPqfj8ZaqVqjqGUAScKaIdHE6plMlIucB21V1odOx1JFzVLU7MAL4g4j0czogL0UB3YEXVbUbcBC419mQvOe+tXUB8KHTsXhLRJoBo4EOQBugoYhc6av9hXQicN9P/xh4R1U/cTqeuuC+ZJ8FnOtwKN44B7jAfW/9PVztSt92NiTvqepm9+t24FPgTGcj8tpGYGOVq8yPcCWGYDUCWKSq25wOpBaGAGtUtUhVy4BPgLN9tbOQTQTuB6yvAXmq+rTT8dSGiCSISFP35w1w/ZLkOxvVqVPV+1Q1SVVTcF26z1BVn53l+JKINHQXIeC+jTIMWOpsVN5R1a3ABhHJdC8aDARdUUUVlxHEt4Xc1gO9RCTW/V42GNdzTp/we89iPzoHuApY4r63DnC/qk50MCZvtQbedFdCRAAfqGpQl16GgETgU9ffKFHA/6nqZGdDqpWbgHfct1VWA79zOB6viEgsMBT4vdOx1IaqLhCRj4BFQDnwPT4cYRyy5aPGGGM8E7K3howxxnjGEoExxoQ5SwTGGBPmLBEYY0yYs0RgjDFhzhKB8TkRURF5qsrXd4rIn+to2+NF5JK62FYN+/mle1bOmcctTxGRYvdsl8tF5CURcfzvyj2b6I1Ox2GCg+O/sCYslAIXi0i804FU5R6X4alrgBtVdeAJvlfonv6jK9AJuNDD/YsPk0ZT4JQSgY/jMQHM/tONP5TjGgxz2/HfOP6MXkQOuF8HiMhsEflARApEZKyIXOHuy7BERNKqbGaIiMx1r3ee++cjReRJEflORH4Ukd9X2e5MEfk/YMkJ4rnMvf2lIvK4e9mDQB/gJRF58mQHqarlwNdAuog0EpHpIrLIvb3R7m2luK8sxuEaLNRORF4UkVw5rteEuPodPCoi893f7y4iU0SkUESur7LeXVWO88jPjwXS3FcqT55svZPEM959/EtE5Gf/ZyYEqap92IdPP4ADQBywFmgC3An82f298cAlVdd1vw4A9uAaVR2Da275v7i/dwvwbJWfn4zrpKYjrnlz6gNjgAfc68QAubgm8BqAa1K1DieIsw2uof0JuEYLz8DVwwJc8zvlnOBnUoCl7s9jge9wzXUThavfAkA8sAoQ9/qVQK8q22jufo1076er++u1wA3uz58BfgQau+Pb7l4+DFeSFfe/wRe4ppQ+GpcH6x2NB+iBa3bbIz/X1OnfH/vw/UcoTzFhAoiq7hORt4CbgWIPf+w7Vd0CICKFwJFpeJcAVW/RfKCqlcBKEVkNZOF64+ta5WqjCa5EcRj4VlXXnGB/PYFZqlrk3uc7uN4s/1NDnGnuaUwU+ExVJ7knPHxUXDOSVuLqhZHoXn+d/rSp0K/ENY11FK7E1wnXmz7Af6sccyN19dbYLyIl7vmnhrk/vnev18h9nOuPi7G69arGsxpIFZHngQkc+zc3IcwSgfGnZ3HdfnijyrJy3Lco3ZNr1avyvdIqn1dW+bqSn/7uHj9PiuI6871JVadU/YaIDMB1RXAiUuMRnNiRZwRVXYHrzL2HqpaJa8bVI60Gj+5fRDrgukLqqaq7RWR8lfXgp8d8/L9HlDvmx1T1Xz85EFcPjp8sqma9o/G4YzgdGA78AfgVcPWJD9uECntGYPxGVXcBH+B68HrEWly3I8A1/3q0F5v+pYhEuJ8bpAIrgCnADe4zc0QkQ2putrIA6C8i8e4HyZcBs72IB1xXINvdSWAg0P4k68XheiPeKyKJuG4rnYopwNXi6ruBiLQVkZbAfly3kWpa7yfcD/QjVPVj4E8E93TUxkN2RWD87Sngf6t8/QrwmYh8C0zn5Gfr1VmB6w07EbheVUtE5FVc978Xua80iqihmkdVt4jIfcBMXGfQE1X1My/iAXgH+FxcTe0Xc5Jpw1X1BxH5HliG67bMvFPZiapOFZFsYL7rMDkAXKmqhSIyT0SWApNU9a4TrQdUHLfJtrg6lR05SbzvVOIxwclmHzXGmDBnt4aMMSbMWSIwxpgwZ4nAGGPCnCUCY4wJc5YIjDEmzFkiMMaYMGeJwBhjwpwlAmOMCXP/D0stIaWj38bhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "lmod = smf.ols('LifeExp ~ Population + Income + Illiteracy + Murder + HSGrad + Frost + Area', data=statedata).fit()\n", "cpstat = rss/lmod.scale + 2 * np.arange(1,9) - 50\n", "plt.plot(np.arange(2,9),cpstat[1:])\n", "plt.xlabel('Number of Parameters')\n", "plt.ylabel('Cp statistic')\n", "plt.plot([2,8],[2,8])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the hat values. Although the parameter estimates are the same, the hat values are different from R." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "NY 0.225227\n", "HI 0.239792\n", "AK 0.247279\n", "NV 0.288609\n", "CA 0.384759\n", "dtype: float64" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ Population + Murder + HSGrad + Frost', data=statedata).fit()\n", "diagv = lmod.get_influence()\n", "hatv = pd.Series(diagv.hat_matrix_diag, statedata.index)\n", "hatv.sort_values().iloc[-5:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Leave out AK as in LMR book." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0. , 0.59232605, 0.66032808, 0.69488547, 0.70867029,\n", " 0.71044047, 0.70730267, 0.70088995])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "state49 = statedata.drop(['AK'])\n", "rss = np.empty(len(pcols) + 1)\n", "rss[0] = np.sum((state49.LifeExp - np.mean(state49.LifeExp))**2)\n", "selvar = ['Null']\n", "for k in range(1,len(pcols)+1):\n", " RSS = {}\n", " for variables in itertools.combinations(pcols, k):\n", " predictors = state49.loc[:,list(variables)]\n", " predictors['Intercept'] = 1\n", " res = sm.OLS(state49.LifeExp, predictors).fit()\n", " RSS[variables] = res.ssr\n", " rss[k] = min(RSS.values())\n", " selvar.append(min(RSS, key=RSS.get))\n", "adjr2 = 1 - (rss/(49 - np.arange(1,9)))/(rss[0]/48)\n", "adjr2" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('Population', 'Murder', 'HSGrad', 'Frost', 'Area')" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selvar[adjr2.argmax()]" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adjr2.argmax()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Scale data and strip plot:" ] }, { "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", " \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", "
PopulationIncomeIlliteracyLifeExpMurderHSGradFrostArea
AL-0.142867-1.3345521.541248-1.3760232.113047-1.476772-1.641325-0.237101
AK-0.8782253.0892950.546895-1.1803731.0732161.6998880.9238535.868329
AZ-0.4603150.1548591.044071-0.2472720.1154760.624326-1.7384910.505283
AR-0.483394-1.7389611.209797-0.1644970.744848-1.651863-0.766833-0.222457
CA3.8355281.114921-0.1160080.6256290.7995761.187120-1.6413251.013678
\n", "
" ], "text/plain": [ " Population Income Illiteracy LifeExp Murder HSGrad Frost \\\n", "AL -0.142867 -1.334552 1.541248 -1.376023 2.113047 -1.476772 -1.641325 \n", "AK -0.878225 3.089295 0.546895 -1.180373 1.073216 1.699888 0.923853 \n", "AZ -0.460315 0.154859 1.044071 -0.247272 0.115476 0.624326 -1.738491 \n", "AR -0.483394 -1.738961 1.209797 -0.164497 0.744848 -1.651863 -0.766833 \n", "CA 3.835528 1.114921 -0.116008 0.625629 0.799576 1.187120 -1.641325 \n", "\n", " Area \n", "AL -0.237101 \n", "AK 5.868329 \n", "AZ 0.505283 \n", "AR -0.222457 \n", "CA 1.013678 " ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.preprocessing import scale\n", "scalstat = pd.DataFrame(scale(statedata), index=statedata.index, columns=statedata.columns)\n", "scalstat.head()" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD7CAYAAACYLnSTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd8VFX6+PHPlMxMkkkvJISEFkqoofcmCiwsYgFEEfuq6KqrfnWx/FxXd13dXXRXXcCyWFARRUW6UqT3GgghpEBCQnov02d+fwwpk5kUIpk7yZ7368Vr9965k3mEzDNnzj3neWQ2m82GIAiC0K7JpQ5AEARB+PVEMhcEQegARDIXBEHoAEQyFwRB6ABEMhcEQegARDIXBEHoAEQyFwRB6ABEMhcEQegARDIXBEHoAJTueqFTp06hVqvd9XKCIAgdgsFgID4+vtnr3JbM1Wo1cXFx7no5QRCEDiEpKalF14lpFkEQhA5AJHNBEIQOQCRzQRCEDkAkc0EQhA6g1TdAP/jgA3bu3InJZOLOO+9k3rx51zMuQRAE4Rq0KpkfPnyYkydPsnr1anQ6HStXrrzecQmCIAjXoFXJfN++ffTu3ZvHH3+cyspKnn/++esdlyAIgnANWpXMS0pKuHLlCitWrCArK4vFixezdetWZDJZo88xGAwtXi8pCILg6axWKzk5OVitVjp37oxCoZA0nlYl88DAQHr06IFKpaJHjx6o1WqKi4sJCQlp9Dli05AgCB2FyWRi5cqV5OTkABASEsJDDz2Et7f3dX+tNt00NGzYMPbu3YvNZiMvLw+dTkdgYGBrfpQgCEK7k5SUVJvIAYqKikhISJAwolaOzKdMmcLRo0eZO3cuNpuNV155RfKvGIIgCO5iMBhadM6dWr00Udz0FAThf1W/fv345ZdfqK6uBuzTyAMHDpQ0JrcV2hIEQegofH19efjhhzl+/DhWq5WhQ4cSFBQkaUwimQuCILRCYGAgU6dOlTqMWmI7vyAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIHIJK5IAhCByCSuSAIQgcgkrkgCEIH8KuSeVFREZMmTSItLe16xSMIgiC0QquTuclk4pVXXkGj0VzPeARBEIRWaHUyf+utt1iwYAHh4eHXMx5BEAShFZStedL3339PcHAwEyZM4MMPP2zRcwwGA0lJSa15OUEQBKEZMpvNZrvWJy1cuBCZTIZMJiMpKYlu3bqxfPlywsLCGn1OUlIScXFxrQ70RGYJ+1MK6R/lzw19O7X65wiCILQnLc2drRqZf/nll7X/f9GiRbz66qtNJvJf65tjl3l+bULt8e8mdOelWf3a7PUEQRDam3axNPHDPekOx58dzEBvskgUjSAIgudp1ci8vlWrVl2POK6JzO2vKAiC4Nnaxcj80Uk9HY7vG9sNjZdComgEQRA8z68embvD3GFdiA3Xsj+1kH6d/ZnSRyyHFARBqK9dJHOA+OhA4qMDpQ5DEATBI7WLaRZBEAShaSKZC4IgdAAimQuCIHQAIpkLgiB0ACKZC4IgdAAimQuCIHQAIpkLgiB0ACKZC4IgdAAimQuCIHQAIpkLgiB0AO1iO/++lEKOZRQzNCaIib3brm66IAhCe+XxyXzF7jTe3HK+9vi56X14fEqshBG1UlURnFwFJh0MXgDB3aWOSBCEDsTjp1kaNqZoeNwuGCrho8mw/U+w+034YCIUt8P/DkEQPJbHJ3O5rOnjdiF5C5Rm1h0byuHkl41fLwiCcI08Ppkvnuw4pfLY5HY4xaJUuTindn8cgiB0WB4/Z/7g+O4M7hLA0UslDOsaxMjuwVKHdO16z4DIwZBz2n7sFwlD75E2JkEQOhSPT+YAw7sFM7xbO0ziNZRqeOBnSN4ExmqImw3eotGGIAjXT7tI5h2ClwYG3C51FIIgdFAeP2cuCIIgNE8kc0EQhA5AJHPBSZWpijMFZ9CZdVKHIghCC4k5c8HBnqw9PL/neapMVQSoA3hn8juMiBghdViCIDRDjMwFB28cfoMqUxUAZYYy3jrylsQRCYLQEiKZC7UsVgu5VbkO565UXpEoGkEQroVI5kIthVzBDTE3OJyb1m2aRNEIgnAtxJy54OD1ca/TRduFhMIEhnUaxsODHpY6JEEQWkAkc8GBr5cvzwx/RuowBEG4RmKaRRAEoQMQyVwQBKEDEMlcEAShAxDJXBAEoQMQN0AFQRB+pWPHjnHu3DmCg4OZOHEi/v7+bo9BJHNBEIRf4fDhw2zZsgWA9PR0MjIyeOyxx5DJ3NvjUkyzCIIg/AqJiYkOxwUFBeTn57s9DpHMBeFXMhstmE0WqcMQJBIQEOBwrFAo8PPzc3scYppFEFrJZrOx79sUzu7JRi6TEX9TDKNu7iF1WIKbTZ48mczMTMrKypDL5UydOhUfHx+3xyGS+f+4Un0pScVJxAXHEajxnL6kNquViu3bMaSkoJ04Ee+BA6UOyUn6qQISdmYBYMXGsc2XiOoTRJc+QRJHJrhTSEgITz75JDk5OQQEBEgyKodWJnOTycSLL75IdnY2RqORxYsXM3Xq1Osdm9DGfsn8hef2PIfBYkCj0PD3iX9nSswUqcMCIPfVP1P6zTcAFL7/H6KW/hP/mTMljspR4eVKF+cqRDL/H6RQKOjSpYukMbRqznz9+vUEBgby1Vdf8dFHH/H6669f77gEN/j70b9jsBgA0Fv0/OPYPySOyM5SVkbpd9/VnbDZKFr5iXQBNaJL3wZJWwZRIpELEmnVyHzGjBlMnz699lihUFy3gAT3KdQVOhwXVBdIFEkDMpn9T31yz7tXH9U7iEl39eH0jsvI5DKGzehKWLQ0X7EFoVXJ3NfXF4DKykqefPJJ/vCHPzT7HIPBQFJSUmteTmgj44LHsaNgh8PxgdMHyKjOIFYbi1aplS64m26Cq2t3kcvRz5jukb8/ijAYusC+msFKCUlJJRJHJPyvktlsNltrnpiTk8Pjjz/OXXfdxdy5c5u9Pikpibi4uNa8lNBGTBYTn537jNP5pxkcPphQTSh/PvRnzFYz3kpv/jXlX4ztPFaS2Gw2G5W7d9tvgE6YgKZvX0niEASptTR3tiqZFxYWsmjRIl555RXGjBlzXQNqzpVSHf/5JZXsUh2/HdSZucOkvenQUVhtVqZ8M4VifXHtub7Bffl29rcSRtW+mY0Wygv1BEb4IJe7dzeg0HG0NHe2applxYoVlJeXs2zZMpYtWwbARx99hEajac2PazGL1cbdHx8mvdDecHhXcgE2m415w6Pb9HX/F5itZsoMZQ7nPGYOvR26dKaQ7Z+cw1BtRhukZuZjgzxyPv3kTxs5tuEHZHIZI+fMZdDUGVKHJLRSq5L5yy+/zMsvv3y9Y2nW2eyy2kReY/3pKyKZX6PLFZe5WHaRrv5duVR2iX4h/QjzCWNat2lsubil9rqbe94sYZTtl9VqY9cX5zFUmwGoLDGw/9sUbnlmqKRx6SsrsVot+Pjb5/gvnzvDzpUrah/f9uH7hHXtTmRsH6lCFH6FdrVpaHtSrtO5zgHeEkTSfn2W+BlLjy3FRt3smlKu5K/j/srr416nV2AvEosSGRExggV9FkgYaftlNlioKjM6nCvNq5YoGrtdn3/Eya0bsVqt9B07kRmPPc3lxASn6y4nnhHJvJ1qV8n850Tn4jULRolReUtVmap4/+T7Dokc7FMsS48tZWaPmfxu0O8kiq7jUHkrieoTRHZy3cqW7vFhksWTeTaB45t+rD0+v3833QYPJaJnb6drI3r2cmdozTJmV2LVm1F3C0CmEPcdmtKuknmon4rkvLpjby85fTp53jykp6o2VaO36F0+VmwoxmqzIpd53nru9iQ3vYy0kwXE9AvGL0hN0ZUqouOCGDGru2QxFWVnOp/LymTiwvsZecs8TmxZj0wmZ8Ts24gZMFiCCB3ZTBb0F0qo2H8FY7r9Po4yzJuwRwah0Kokjs5ztatk/uy0Phy9eAijxQpAfHQQPqp29Z8gqTCfMMZEjuFgzkGnx2b3mI1cJsdsNXMk9wg+Sh/iw+MliLKOKS+P3D+/hu7UKXyGDSPilf+HMky6EW5zLiYUsnl5AjVffLr0DWL+iyOkDQroNngocoUCq6WusmOPIfa4Jtx5L+Pm3w2A3AM2/1kqjOQvP42l2HHQYS7QUXkwh4CbukoUmedrV5mwTGeqTeQAB9OL2Hk+jxv6dpIwqnpsNsg8CFYzdB0HcunfHA29PfltVp1bRUpJCmqlGoPFwKDQQSyMW0iZoYx7t9xLWlkaAOOjxvOfqf+RbLR+ZckSqg8eAqBi2zZ0CQl0+/YbvMLDJYmnOWd3Z1F/BivrfAnFOVUER/pKFxQQFNGZW577fxxe9w0Wk4n4GbMpL8znpxX/JqJnLwbeMN0jEjlA1ZFcp0Rew1ppdHlesGtXyfzclXKnc18dymRsz1A0XhL/MpqNsOoWyNhvP46Mh/s3g0raN3JDWpWWxfGLXT629tza2kQOsC97HwevHGRc1Dh3heeg+tBhh2NzXh7ZzzxDty++kCSe5ihd/A4qvTxj2qr7kOF0HzIcgP3ffMmh71YDcPaXbeRfSuem3/1eyvBqWXVm1w/IwGeIZ36IewrP+E1roXGxoU7ntp/P55FVxyWIpoHzG+oSOUDOKUhYI108rVCid96KXn8Tkbtp+vVzOqc7dpzs559Hd/q0BBE1bci0GJSqurdU3zER+Id63mqrMzu2Ohwn7tqOxdxIEnUzn6HhoKy70SnzkqMZGErogwNRdwto4plCuxqZx0cH8o+5g3jlx0R09Tq77L5QwMXCKrqHSjgKripyPpe2E4Y/4P5YrsHFsouklaYxvNNwZvWYxZfnv8Rstb+xg9RBTIqeJFlskW/8lUsL7sSm0zmcL1+/gfItW+n+zRo0EpeIuHyumH1rU6gqNdBrRCfufGUUWedL8AvVeGwpXLWPL1WldR/cXt4+yD2kkJmqs5bwxfFUHctFrlbiOyYSZYBa6rDaBc/4F7wG84ZHO43QZTJQKSX+T4n7LdBg6VTWUUlCaamVZ1dy87qbeXrX00z7bhpVpio+nfEpt8beyl197+LLmV/ir3J/l/Eamj596P7NGrxc1Yk2mShbv8H9QdVj0JnZ8sEZiq9UYag2c3Z3NinH8ug3vjPRfYPd3tC3pfpPnlpXlVImY8KCe5B5SDIHUEVp8b8hBqveTPGa85T8mIoxu0LqsDxeuxqZ13h0Ug/2pRagN9lvhs4b1oWoQIm/zvp3BqUGzPVGkSYdHP0vHHgXZHIY/zQMvUe6GOupNlWz/NTy2mOdWcey08tYOX0lg8OkX55WQ92rFz1//omyH34g5yXHXceW8nJK167Fd8JEvDq5fz61ILMCk8Gx9+eVlDKGXd0RfzGhkP3fplBVZqD3yAgmLuiNQsJBR2luDqe3bebElvX2m/WAQqEkesAgyWJqTOFniZiy7M0/jOnlVB3MwSvSh/Anh3rsh6TU2mUyH94tmJ3PTuaX5Hy6BvsyLjZE6pDsRjwIB9+vO+41DTY9U3e8/gkI7w9dhrk/tgZ0Zp3TmvNSQ6lE0TRNJpejjIhAGR6OuabruUxG2XffUfbdd8g0Grp++gne8e5dShnaRYvSS47ZVLfCqlN3+zcZfZWJnz86W/vYuX1X8A/VMGxGN7fGWCP/UjqrX3kOs8HgcN5iNpFy+AAj5zRf+dRdzMX62kRenymnmrL1aQTOiZUgKs/nOd+trlHnQG8WjurK+F6hnvNJfdPrcMsK+zz5bR9BcE/nay7tcX9cLoR4hzCpi+N8+K2xt0oUjWs2i4XKvXspWPEBlx98qC6RQ+3IEsCm11P48cduj0/j68VND/RHG6RGLpfRe2Qnhk6LAaDgcoVDkgfITStz9WPc4vS2zU6JvEbmmVOUF3pOUTW5rxcylevUpLsg3Q15T9cuR+YeSy6H+DvtfwCStzpfY6wGiwkUXu6NzYWGI/NKo/NoSCqm7GzS583DWtyyZg+2al3zF7WBHkPC6DEkDKvV5lDmNqyLHwovOZZ6CT2ip3SrMWRN7BXIOHOKlU89zPDf3sKYeXehUEr7uylXKwiY1YPSdak0qDyBpchA9qsH8JvUBf8pMdIE6KHa7cjcY1mtYLh6s6bPDPs8ubLe3fg9f4dVtzqMLKWQWZ7J4RzHddwfn/kYq83ayDPcq+C991qcyAGCFt7VhtE0r2G98lPbM7FZ7P/GMrmMuHGRxE+VLvnET5+FyrvuvpLK28fhcYvZxOF137Lv61XuDs0l7ahIIl8aBS7W6dv0Fsp/ykCfKro61deuRub5FXqW/ZJGZnE1MwZEMN/TSt+mbIcNT0J5NmgCoXM8THweqovhxGd1113aC5cPQ8xoyUItNzpvwDJZTVhsFo+oz2LMvNzoYz5jx2IpK8N48SKq6GjCX1iCdrR0f5f1mU0WLp4u5PjWjNpzNquN7oNDUUi4gSg0uiv3LV1OyuH9ePv5k3U+kYTtzt8cU48cZNLd0i2ntdlsVB3NRZdQiFdnXzS9g9Anulj2CxgvV6CJ9czln1JoN8ncZrNxz3+PcD7XPurdeT4fo9nK3aM9pFaD2QBrFoL56tSFvhTSd0HGQeg7y/l6i8mt4TUU5h2G1ktLpaluaqV3UG+O5h5ldORoyRN64G23ojtxwuFcwNzbCV64UPK15a5YrTb2rrnAuX1XcNW8K/lgLt0HSVtXxi8klKEz5wAQFdefpH27Mekdp6cCIztLEVqtsp8uUbkrCwBDainyIDXqHgEY0p3vN4hNRI7aTTK/kFdZm8hrrD91xXOS+ckv6hJ5fRYDBMZcXbZ49fHIeHvtFokcuHKAJ3Y8gdFqr3Xhp/LDZrORXJLMI9seYWTESD686UMUEtaWCZw7F0tVFcWffgYWC4Hz5xH6+OOec7O7gYSdlzm7O7vRx7XBnrXxxT803OWHjpSjcoCqw449C6wlBgIWxuEV4UvloRwqdmeBDPwmdkHd3bOSeVlZGUeOHMFgMDBkyBCioqLc+vrtJpmHalV4KWSYLHW/gBEBbdum7pqUZDT+WPeJEL8QEr8H31AYtMB+s1Qiy08tr03kYF+mWLPrE+BI7hH2X9nPxC4TpQivVsi99xJy772SxtASFouVoxsvNnlNeFfpNl81RhsURGluTu2xf2g4odESD45cfFabK4youvjhNz4Kv/HuTZAtZTAY+Pjjj6mosA84T548yUMPPURkZKTbYpB+crSFfNVK7h7dlZr7TBH+Gp660YMK6XdtpLH1kLuh5w0Q1hsmL4ERD4Fa697YGqgwOn7DqZ/IG7tGCrqEBDLuXkTqDVPJ+8c/sHlI/ZCG8i+WY9Rbmrwm5bhzYxWpTb7nIZRe9vrgCi8vJt/zkMQRge+ICKdz1jLPr5aYkpJSm8gBLBYLp06dcmsM7WJkvjelgCdWn6S02kSQj4pnburNHSOipd/CX1+f38CUl+DQMlCoYfj99kQeUG8renWxvYqiUtqv3HN7z+Wto2/VHo+KGMXJ/JO1o/Uw7zCnNejuZjUYuPzoYizF9nXFxf9diSIgkOC7F3Llj0uo2LEDua8vob//PSH3SrurVhusQSZreoFSRaGeyhID2iDPmW7pOWwUv1v2CfkX0wjv3rO2N6iUFH7OzScMaaVoR7tvhHutDAYDly5dcjrv7e3eXentIpm/9MNZSqvtNwxLqo18c+wyi8Z4yFx5fZOet/9pSF8Oa++H1O2gCYBpf4Whi9wf31U3x97MJ4mfkF9tHy1WGCv4dManbEzfiI+XD3f0uQOtStpvD4akpNpEXqNq/34spSVUbNsGgLWigvy//Q11zx5ox4+XIkwA/II1DJ/ZjWObLzWa0IuyK/nqz4e4/blhhERJ+3cLcDkxgRNbN5CblkJ1aSmRvfow/dEnCYqUdhpDn+68C9mUU4nNYkWm8KDB21U2m43PP/+c7GzH+yVBQUEMHz7crbF43t9OAxarjawSx2a4l4qqJIqmlfa9Y0/kAPoy2Pg0VOQ1/Zw2tCFtQ20iBzhXfI7LFZd5YdQLPDX0KSJ8nb/qupuqe3dkGsd7Ipq+fanat8/p2uIvvnRXWI0aObsH97wxltueG+pQBrc+k95CwtWVGlJKP3GEb15/idQjB6ksKsRqMZN9PpEt/3lb6tCQu6gJby7UU3XEuZm7J8jKynJK5DExMTz22GNote790Pb4kblCLmNqXCe2natLftP7S59snOjLYO/bkJdonyMf9Uhdp6G8RMdrrSYoTAY/aTokudrpWX+Joicw5eWh7tULfVISmM34jh9P6GOLsZSWYriQ4nCtylVVRQnYbHB2dzZmY+Mbr6wWaTaL5V9K5/jmH8lOSqQs33VizElJxmq1IJdwFZN2bGd0CQVOOz9NVzxzAOfl5bxbNjg42OX5tubxyRzgn/MG8/bPyZzKKmNU92D+cGMv9CYLKoXcaeedZNY+UDf6Tt0GyZsg/zzIFPbNQ/VpAiFKumJbv+n+Gz4+83Htdv5gTTA3dr1RsngaspSXc2nuPGzGuhtfqj69kWk0RLz6J6qOHsV85QoA8sBAQh64X6pQATAbLWz96CwZZ1xvbqkhV8gYOMn90xg5qcl8/crzDj1AXenUs5ekiRxA3dWf0IcGUvjxGYeEro4NlC6oJkRERNC3b1/Onz8PgFqtZsyYRhZDtLF2kcwDvL3485wBAOhNFv7v29NsOZtLkI8XL8/qxy1DJF6uVF1cl8hrXKo3HZDyk32tOXIIiIIb/yxpO7m/H/17bSLXKDQsnbSUYE2wZPE0VPTpZw6JHKDkvyvRHztO19Vf0WvnDqpPnMBSXIzvuHHI3XyjqaFz+3OaTeQAY+fGSrJE8cyOn5pN5AB9xkxwQzTN0/QMJGRRP8q3ZdjbyNmg5IcUKo/kEDQnFq9wn+Z/iBvNnz+ftLQ0Kioq6N27t9unV2p4/Jx5Qx/tSWdjQg4Wq43CSiPPrT1NfoXrBrBuo9LaR9tNKc20L028fzNES9ex/VT+KXZn7a491lv0bM/c3sQzJGB2vTtWd/o0Oa+8QtmGDXhFROB3442SJ3KA0rzqZq8JivRlgFRrpFu40aosz3Pmpb37hRD2yCAs5QYsZQZsegvGtDLyP0rAZpW2rlFDcrmcXr16MXToUMkSObTDZJ6Q7bit12SxcT5H4jXRShVMf6P5SojFaU0/7gY6s3N1wSqTZ81HBi1cCI10iy9b+x1Xnnue1Kk3UvzVV26OzLWaGuYNKZRyBk+NZtJdfZj7x2GS1Wbx9W/ZFEVkrz5tHMm1KfkhBRrcfrBVmDAXNP/h2ZYsFguHDh1izZo1HDhwALOH7H9od8l8TA/HRhQ+KgWDoz1gPm3IQhjcROU+mRz6zHRfPI0YETGCHgE9ao+VciW397odAJPFxPaM7WxO30y1Sbo3jFenTnR5/31kTY1ybDYK3n4HayM1ut1l15fn2f7pOQCUDZK1xWzFarYyYGIUKo10M5rqJv4eFSoVSrWaYb+9lX4TprgxquYZ0p2LwQEoAqXd+f3TTz+xdetWkpKS+Pnnn9m8ebOk8dRoF3Pm9d07tht55Xq+P5lNuJ+aF34TR4C39LXBAcg85HxOJofOQ2HCsxA11P0xNaCUK/lsxmesTVlLsb6YWT1m0T+kPwaLgUWbF5FUnARAlDaKr2d9TWBz00dtpHzLZmyVTa+wsVZW2ps9q6XZiJOTWkri3iu1xw2bUQDIPWBjW9z4yRzb8L1DE+caFqORaY88ycAbpkkQWdMUfl5YKxzvnfgM74RcLe1N2tOnTzsdz549W/K6QdL/pl0jhVzGCzPjOPrSjWx6cgLje4U2/yR3SPjGvtywIZsVitOhm3SFtRoK1ATy0MCHeH7E8/QP6Q/AzsydtYkcILsym3Wp66QKEX3CmWav0U6ejCJQum9l5UVN36tR+yoZIMHqlYZ8A4NY9Na7DJt1i8vHf/7gXb557UXMRs/aNu832bHEtaqbP8Fze0sUTb04VI67VBUKBTKZjOLiYk6cOOG07txd2t3I3GM1XM1Sn64YflgMpmroNh7GPmmfZ/cgrubSXZ1zF58RIzDW2yKt7NKFkAcfoPrQIax6A96DBhJy332SxQcQ0z8YpUrutK48INyb+Btj6BEfho+/Z/w7+wYGMfmeh8g+n0huWorT45cTEzh/YA8DJnvGElWb1Ybcx4vAW2KxlBsw5lRhuFjGldcOookLIWB6VxT+0nwja5jMTSYTZ8+e5fvvv8dqtf8uTJkyhUmT3FsSQyTz6yW8X9OPJ2+y/2/6L1BVAL95q+nr21CJvoTV51dTpCvi9l63o7foqTBWEKwOpthg30Kv9dIyu+dsyWIMf/459Enn0J9NBJkMTd8+BM2dS/Cdd0oWU0PeWhWznxjMundOOaywiI4LZsBE6Ufkrsx66o+s+uOTGHXO90QqiwoliMiZ1WCh4IPTdRuFZNSuObcC1cfzqD6RR/BdffEZ6P4a8YoGN+dtNhs//fRTbSIH2Lt3L2PHjnXr5qF2m8zNFivLd6Wx60IBvTv58cxNvQnzk7CI0ciHIeOAfU15c86slSyZXyq7xPyN82tH3d9c+Kb2MYVMwawes+js25lbY2+li590OyvNBYXoE+03FrHZqNy+g9K1awnyoGQO0LlXENMf6s/u1cnoKkxE9Qli1OwezT9RAid/2sjR9d+h8vGh+9DhpB4+iOXqMlCFUkmvUZ4xFVh9Is9xx6erlYg2KNtySZJkPnLkSDZu3Ohwrn7FRLCveLFYLCKZt8TSbRdYvsu+1O94RgkX8ir4bvFY6QJS+UB4XMuSeYB0SfLLpC8bnT6x2CxkV2Tz5oQ837OyAAAgAElEQVQ33RyVM/25c05lCHWJiXhik7CeQ8PpPjgUk8GC2sf+5jXqzCTsyqKsQEeP+DC6D5L23k5W0ll2rlxRe5y8fw83PfwElxMTsFmtDJkxm5AuntGG0VrVsi5c1kppunUNHz4cLy8v1q1b57LBB0D37t3RaNy76qbd3QCtsfWs4waH4xkl0m8eyk9q/DHl1c0t6gCY/lf3xOOCwdL0Ur7mHncXn2FDocGoxlJahrmgQKKImiZXyGsTOcDG/5zm8I/pnD+Qw+ZlCSQflnZDTta5s07nkg/soVP3nkx9cDFRfZuZJnQj78FhoGx+ZYjvcGlqGwEcOHDAKZHXn35JT09n7969bo2pXSbz5bvSyCh03ujSsLqi2/W6qfHHfvcL3L8Fnk2y3wSVyLze85rs7zk5erL7gmmCV2QkXd55G6/outFi5fbtXLrzLsnXljenNK+anFTHzW3n9l1p5Gr3iIh1XgWSefY0u79YyZcvPoPJIPFAqB6vMB80/UIafVwZ6UvgzT0JmCXNdFZxcTF5ec5VTxvOpe/fv99dIQHtMJmfvlzKW1vPN9wYBsCOJIm7uXSbYO/12ZBCBctHw7rFTY/e3WBg2EC+nf0tN8TcwPBOw3lt3Gv4KuvqxCw/vZztGZ6xvd/vxhvxHe84j2vKyqLKzW+Sa+WlUSBrUABO7SPtjGa3wUMZM/dOvDTeyJWOsZTl5/Le/Xdw8qeNjTzbfWwWG0Wrz6NPcH0zVh6gInzxYLRjOyNTSLOu29fX12lFS69evZzmx9297rzdJPPU/ErmrTjA7csPNHqN5D1B1z7QSFPnq+t3Sy7BD4+6NSRXegf15t9T/s0nMz5BJVdRZXb8lvPO8XckisyZQuvndE7uK31zh4ZsNhtWi32I4RugJv7Gum8UKm8lw2d2kyiyOmPnLeT3n3zNyJvnOj1ms1jYuXIFBZlN9LJ1A93ZQnSnG5lKU8oIf2QQcpXElR3Van7zm9+gvPqhqNFosNlshIU53ozt1auXwwqXttZuboD+Yc1Jzma73t5bQ91IPQ+3yDwE+YnNX1eUAiY9eHlGM2pX2/Zd9QSVStDCuyhbvx7z1a+1vhMm4DNSukJlrpw/mMOB71MxVJvpMyqCSQv7MPa2WHqN6ER5gY6oPkFofD1jl7JcrmDQTTM4s/MnlztCD69bw2+fdNEty03MhU3sbTDbMGZWoAyWvrjakCFD6N27Nx988AHl5eWkpqY6XZOQkIBCoWDOnDluialdjMx1RkuziRygoFLCudTMgy27LrBrXdMKD3Bbr9vwVjq+Oe4fIG198Pq8IiLouWUzUf/6F1H/eodOL74g+bbp+iqK9excdR5dhQmrxUbSgRzO7rLvAAyL9qPn0HCPSeQ1/IJDuW/pclQ+zmWYfQOlLYWsiQtuMitJtVHIlby8PMrLm85LJ0+epLKZshTXS6uTudVq5ZVXXuGOO+5g0aJFZGS03dczb5WCuMim60DLgBkDJOxA1GVky64rzYC3+0HqjraNp4UUcgWbb93MtK7TiA+L5+1Jb7Og7wKpw3Ig9/FBf+4c2c88S/pvZnJh4kQMly9LHRYAuellTiVZ8y41P/CQWm5qMr4BjqUQlCoVY+ctlCgiO1VnLSH39EfmovWed3wY6h7SN52u4ePTsrrqV6645+Z3q5P59u3bMRqNrFmzhmeffZY332y7tckGs4UBnf1RXy1a5Gpc5q1SkJBVilWqWsfdxsFNr4G8BTNXVfnw4+/B2nzDAHcI9Qll6eSlrJq5ipu6NbEiRyL65GSKPvwQrs4/WvILuHjzHHTnkrBUSFf+2Kgzc3h9utP5zr08oIpnE7LPn+O7N1+lJKeuhohPQCCPLP8ctQfUh1d10WJrUCJBplFgLtRR+Hki5jLpVzMlJyezfft21C0o8paTk+OGiH7FnPnx48eZMMHemSQ+Pp6zZ53XsV4v72xL4dvjdY1wQ7QqyqpNmOol7mqjhafXnOZUZmltVyK3G/cUpP1i37LfnIor9r6hPp7T4cdTGV1867PpdFy67TZk3t5E/OkVAm9xXUSqLV04mkdZvuMcb0z/YPqP7+z2WJqTkXCKoxu/Iz89DV2F8zeH6rJSygvz0UjYXKGG3McLRYAaS72kbdNbMGVVYgLyUo7R+ZXRyFw0f3aH5ORkVq9e3eLrO3Vyz3r4VifzyspKh64aCoUCs9lce4e3IYPBQFJS65blbU1w7GheWNl4dbevDmdye085XhItW/KJvpWYS3uRNXMTsTpkABkZeYDzelV3KDYWszFnI/mGfIJVwcRqYxkaOJRDxYdIqkgixieGqWFT0Sqlf3MT1Pi+T5tOR85rr5PTrZvbS+Fmpjt/K/AKNHE+2d4P0lBlQeUtd1qm6G6lWZkc+O/7TjtqG0pLSaFIJ/2oF0AxWoNmtwG5i8VhNpOV1I2nMfeV5lvEgQONr6hrKCgoCKvV2urcdy1ancy1Wi1VVXVL2qxWa6OJHOzLeeLi4lr1WgNO6rhU0rKvKhqVgv794lBI9QaKi4NB4+H0ajj8IRjK7G27Rj0GNguk7YSIgfjc9DpxAdIUY0orTeOpDU9hsNZ74+bb+4HW9AalCLYWbGXtzWsJ9Xb/VnTj5cvk/3Mp+vPn0fTvj2zePMq//db1xdXV9AwNQ9XFvX+flekXuYjjza3Izp2IDA5lywdnKL5ShUbrReywcAbfEE1gJ2l6V+46uq/ZRA7gK6fV79HrLg7yzpxwrNFST+duXfCNC3dzUHaXLl3icgvv2UyaNIl+/X7d7tqWfhC0OpkPHTqUX375hZkzZ3Lq1Cl69267OsNLftOXtIIqknKav7H05A29pEvkNUJ6Qu5ZeyIH+xtJJoMZ0lVKrG/1+dWOifyq2kR+VZG+iPVp63lgwAPuCg2wr9m+vHgxxlR77R1TRgYoFKhiYzHn5mJtsDpA3S/O7YkccNkG7mJCIUc2XsSos98P0VeaOLs7m8S92cz+fTzR/dw/raYNbtmHcWmutLtUG7LqXd9Tkvt54TPI/QW2aowdO5bk5GRKS0ubvXbHjh3Exsa6pTdoq2+A3nTTTahUKhYsWMDf/vY3XnjhhesZl4MuQT5seWoCQ2Oav7H0u4kSVqyz2eyNm3UlcGGr42Nnv5MmJheuZR25RYKbtKbLl2sTeV0gFoypqc6JvE8fopcvd2N0dfqMjHBYdqhQyslOLq1N5PXZrHDiZ2k25Ay8YTqRsc339+zvIbXMa8i9ncea2huiiVwyUrLdnwABAQE88cQTzJ8/n5iYmCYrI1ZUVPDZZ5+5Ja5Wj8zlcjmvvfba9YylUZvP5PCXjefIKWu6fkTXYGm+xgJQmApf3wmFF8AnFNT+dSNzkLRSYkML+i7gx9QfMdsck7qX3AuTta4SnZ/ST5Ka5srwcGTe3vaWcM1Q9eyBl5tuMDXkG6hm/ksjSDqQg1wORzddavL66nJpOvmofXy48y//ZO9Xn3J0vfOgQq5QMG7+IkKiPKNqYg1zqfP73X9yNDKF9NtjFAoF/fr149y5c2RmZjZ5bUFBAeXl5fj7N728+tfy+B2ghZUGnvr6JCZL03N+gd5e/GPeYDdF5cLGp+yJHKC6EFRaUKjBYgDvIJj2F+lia6BvcF9+vOVHPkz4kJyqHEK9Q4kPjydEE8Lze57HYrOPLOM7xRPh6/61+3KNBp8RI6jas6f5i6VainqVX7CGETO7kbg3G7lChrWJ39OS3GpK86olmTuXyWR07tMPcEzmY+ctZNisOai8JRwINcIr3AfjxbqpVUWgGpmLqS2ppKSktGgVn1wud0tdc49P5olXyp0SuUYpR292XIdqtdkY1EXCDQVZxxyPjZXw+2NQXQQRg+z1zj1IjH8Mfxnv+AHzxM4nahM5wN7svWSWZxLjH+Pu8FAEtOzfUq6Wvi3b1o/Okn7SsZ6IUiWn2+AQUo/WnbdZbWQkFkl2I7TH0OH0HD6atGP2xuNBkZ25cHg/5/b+Qt+xExh9+wIUSs/ZrRrwm+4UfXEOa7kJua8XQbf38pjdv2azmR9++KFF144ZMwZvN6zf95yPuUb4uFhL2iXIm7hIxwJM5Xoz3xyTcFegq9UC/lEQM7oukaf9Agf/A3nn3BtbC7kqtG+1ua9QUH0Bc+bYbxo3RSZDO3mKewJqhL7S5JTIAcxGq1NJVECyRA72uiy3PPcy9/z9PeKnzaIk5wqFmZcozb3Coe/XsO/rVZLF1pAuqYiiTxOxlptQBKoJfXAAml6e05rk4sWLVFc3X3I7JiaGm25yz0Y8j0/mvmrnLw9pBVWkFzjXO7hcLGE9897THI+1nRwrKP78Mqy6BX56EVaMg3M/ujc+Fy5XXOaTs5+wIW0DiYWJmC2Oc+iToyfTLaCbJLFpx48j+qOP8B7huqiWMjISZDKyn3mGrKefxiZVZ/kmPm90lUZ6DrUvn5PJoN+EzsRIsJqlobCu3Unav9vpfOrRFtYXamNVJ/Mp+iIJa7X999FSaqDk+xQKPz9H4efn0Kc1v4qkLWVnZ/PNN980fyEQGRnZxtHU8fhpln6dnW8a2ACD2XkUqVVL+BVx9rtQXQIZ++zHlXnw6Sx4dD+YquDwB3XX2qyw923o555qaq4kFiVy35b7nJYj1lDIFDw++HE3R+XIdDkT3dGjLh8z19siXbFlK2XjxxN4++3uCq2Wxte+jjz1uHMt/eBILeNuj6WyRI9MLsM3wHOKRMnlzuO4oAjpd67qkoooWZPsdN6++9M+gNMnF9PpqaF4hUvzLWffvn2YTC1rWRcY6L7SDh4/MgdatCQRYNWhS5gs0kwL4BMMvg3W8+afg8uH7cm74RI/ieuyrE5a3WgiB3s/0LUpa90YkbOSr9e4fsDF9IshzblGiruMvqUHIV20yBWy2pF6ZGwAw6Z3BUAbpPGoRA4wco5jTXO1r5bYEWMovpLdyDPcQ9dIUwoHFhu6xKK2D6YRxmv4FuiqI1FbaRfJfOn8eIJ8mh91F1YaqdBLWIu73MUbQe0HiT9AcIP172Mec09M7ZisQTcXFAp8Ro92eX9CO3mSm6Jy9tNHiRRlVdpXstggdng4t/3fMGw2G8e3XuLgD6kUN7KTUSrDZ9/GHa++ydj5Cxl/170olEq2ffQ+nzz9CPu+/lyyuBRBLfvQU7bwurYwfPhwl+e7dHFeftyW1WQbahfJvHuoL6seHEVzGzvjuwQS7Cvh6gajizfsgXdhw1P2phQA0aPh/q0Qf5d7Y2tgbNRYFDReqMhL7sWCPtKUwrVZLFz54xL0Z844nJdrtdgn2Rxpb7wR35EtLEF8nV06U0hBpmONlty0MsxGC2v/fpxD69I58VMm3/7tKEXZ7qlr3VJd4gYw5vY7KbqcSXVZ3Tz0kXVrqShuwQi5DWjHRaHs7FxnXe5bNyOsiQvGe6D7S0zUiIuLY+DAgU7ns7KynMriRkS4b2lvu0jmAAOiAlhx97AmC2hlFFdRKGWDipCejsdKDSSuczxXeAG6jnFfTC7szdrLC3tewELdVI9aUTfSkSHjtbGvERsUK0V4VOzYQdmPzjeIrWVlVB867HRe1dX9SycBTEYL21Y6d5cK7aIlM7GY8oK6TU9mk5Vz+z1ruzxATmoyxdmOhexsNiu6ZpoutBWFrxednhhCwOweyDT2wYY8QEXI3f3o9OwwOj0zjNB7+0u+cWj69OkulxvGxcXVzpNHR0czY8YMt8XUbpJ5dqmOKqOZVQ+M5NYhUQR4O0+7lFSb2HhawjfMlJdAe/WTWKaAqX+yT7PU56rhs5utOrcKa4OW2AZLvXKj2Nh8cbO7w6plvHjpmq73HjiobQJpRkFmhdPWfS+1gvHze6NzUdlTLnXNoKv0VZWYDHq+ff0lvnrpWfLSUxweD+vanbCu3SWKzr7ByW9cFJFLRqLuE4S13EjBRwlUHcyR7KZnQ1qtlltclF0eOHAgTz31FC+88AIPPvggAS3cL3E9ePxqFoCfE3NZ/MVxmtkECoBKKWFLtvA4ePIkHPsv+IbDgNvs7eSS1tddY6yQvAeoXNb8Z3hLrmkr2smTKHj3XbC4vknsM2oUuhMnsNlsBM6bi980aRpqBEf6ovCSYzHVfTAOnBxFQJg3Z/c4lyLwC5b2g1xXUc7Gf71F5tnTqLx9MOocl/J26hFLzIDBDJ99m0dsztElFmFIvtqn1AaVB66g6ReMJlba9eYWi4X9+/eTmJiIWq3GYLAPhPz9/fHz80Mmk7WoacX11i6S+Yvfn2lRIo8O8mb2YPet63Ri0sNnsyH76m7QA++CpsEns6ECcs9AtDRNiXVmncuStgqZonb3p5fci3v73+vu0Gpp+vShy/vvUfTxfzFkZGAtrJu/VYSGErPyv9j0emw2UGid51fdFqevF1PvjWPP18noK8107hVIaIwfxVeqULsoEuUXKm0XnwPffkXm2dMATokcIHbEGEbfdoe7w3LJmF1J+U7nmifmfB1InMy3b9/OwYPOa/LLy8tZtmwZY8eOJSYmhl69erk1Lo9P5j+eyqawqmVLgaqNlmZruLSpcz/WJXKAvLPQfaLjNUqN89y6m9hsNu7beh/nipx3oNYkcrVczYqbVjA8wvUde3fxmzIFvyn23Z26s4mUfv89XpERBN9zDzKFApmvdEm8Pv8Q79pVpldSSrmSYr+ROHByF/xCNFQU2Zd/RsYG0LW/tBuGCjIaX76p9FLRe/R4N0bTOKveTMFHZ7C5WJnmqjeouzVVj8VisbB3714AQkNDGT9+PPHx8W6Jy+OT+e5k563SjSmqMrLuVDYPjJNovk/vYmdal5H2NeUZ++2VFGe9LVmruFMFp1wm8voMVgPH8o5Jnszr8x7QH+8B/aUOw6UjG9Mx6pyTzpldWcgVMnwD1Qyc3IUh02IknzPvOnAI2ecd//0VSiURsb2ZvOghgjtL0yylPpvJStn2TJeJHMCYVYmvxL+agYGBVLSg92xhYSHr1q2jurqasWPHtnlc0n/MNaN3hF/zF9Wz4ZSEN0D7zQFNvQ1OKi1EDIArp+zHhnL7RiKJKGQtu5+gbElTagHAqQdofVaLjapSA6e2ZWKVajNbPSNvmcuw397qUCHRYjaTm5ZCQCf3V8dsyGa12W907mt845JS4qkqgGnTpjXZVa2h06dPt2E0dTw+md87phtaF/VZGnMmuwyzVG8cvwj43U4Y/TiMfAQe2gHHP7Nv56+x7234dDbo3F9fYlDYIIZ1GtbsdV8nfc2yk8tcFt4SHFlaMK2nrzJRmidh3aCrFEovJi96kMhejo0qLCYTRdkSFqm7ypBehjGz8RGvuncQviOl/9CJjo7mmWeeoWtX+w5fuVzusqhaDV83TQl6/BDMW6Vg1YMjWfzFCXLLm25OARCqVaOUcg1qSE+Y8Ubdsa7E+ZpLe2DPP2D6X90X11WLBy/m4W0PN1kNMU+Xx/KE5SQVJ/He1PfcGF37I3ex7yEsxs9hI5FG6yVptcSGYgYMJiPhZO2xRutHp+7S3Mdx4GLw4B0fRuBve2AzWlFKvBqoPh8fH+6//37MZjNyuZyEhATWrVvndJ1arWbKFPdU9vT4ZA4wJCaI/UtuIKdMR1axjhe+P8PFIufdliqlnFdm/7rmqb/a3qVw5CPw8obJL8DQe2DTKefrchPcHxvw7YVvW1zWdnfWbgwWg8OGIsFRv3GRHFpXd2Mxdng4U+7uy64vk0k/VUBAmDcTF/RG6aKUs1SG//ZW9JUVnD+wB7+QMCbdfT9eaukTpbpnIF5RWkxXd8rKVHL8J0ej0Epfs74xNdMt8fHx+Pv7c+HCBcLCwoiJiaGkpISYmBg0Gvf83cpsbvounZSUdN06fxvNVv75czJrj2dhsdoY3SOYWwZHMapniLTb+ZM2wJq7645lcnjssH0p4rrF9q5DNaa+AhOedXuIS/YuYVP6phZdK0PG8buP46WQphqlKS+P3D+9SvXJk3jHDyby1VfxcmNJ0ZZKPpzL5XPFhEZrGTApyqMSd3tjNVqoPpWPtdqMXKVAd7YQmZccvynRqLtJ2HxGQi3Nne1iZN6QSinnxZlxvDjz+nw4XDeX9jse26z2VSzD74egbrDt/9kbPve/FcY+JUmIi/otYmfmTnTm5vtrDgkfIlkiB8h58SWq9tv/Tqt27+HKiy/S9ZNPJIunMWExfpgMFsK7+jkkcl2FkYsJhfj4q+jaPwSZh+wA9WRylQLtyEj0aaUUflRXm8eQXkbEc8NR+ItviY1pl8ncY0UNbfxcl2Fwv3Rb5Gv0D+nP+lvWsyNjB7lVufyQ+gNlxjKX157IP8G7J97lyaFPujlKu+ojRxocu65tLqXkw7ns+PRc7XTv2NtjGXJTDPkZ5fz4r1O1yxa7DQpl1mPSlB1oj/TnHEvc2kxW9Mkl+I6Q/gaop/L41SztyoC5MGqxfWOQJsBeqyVpA6x/0nnULqEI3wgW9lvI1K5TG03kNb5I+sJNUTnTDBjgcCxTKrn86GIM6dLVLm/o6KaLDvftjm2+xNFN6ax985jD+vNLCYXkZ0hTvKo9crUEsWRdKjl/P4runHS1zD2ZSObXm1oL3sEQEAPHPrGvWjnxGXz2W0jfJXV0DloyheIll26aJfKvf0HTr+6Gts1goHLXLi4/8ig2qxVrC3owtjVDlWPHGaPezJENl1y2hDWbpF9r3l74Du+Epk+DbfsWG5ZiPcWrz2Otblmnn/8lIplfTyc+tSfviiuQd8b+vzVsVjgp3SjXlf4h/ZnYpa7cgI/SBx+l4xK638f/3t1h1VL36EH3779D3c9xhZLp8mXSZ80ieegw0ufcgiElpZGf0PaU6gY3OxtZThAS5Utkj//NG3itIfNSEHr/ACKeG4461rHTmM1kxehhzT48gZgzv17MRvsGoaZ4e0538RrvTnmXPVl7KNAVMCV6Cmqlms8TPye3Kpf5veczKFz6eV5Nnz4YztXbOSuX15bJNSQnc+Xll+m+ppEWc23MW6uisrj5Gvr9J0SJG6CtoAzxRt0zEENq3SY7mZcclYsGFv/rRDK/Xn58DHJcrCevL3EdGCrtm4o8JLEr5AqmxDhuavj9EOlG466E/eEPGFJT0Z85gyIgAEuDuhj6ROlKJAy+oQvbP01q9rqQKK0bomkZq8VC2vHDVJWUEDtiNNrgEKlDapLfhCgMKSUYLpYhU8rxn94VeQvaSP6vEcn8etCVwNnvHM/J5Paplfqq8uH0V2DWwbxP3RZee1X8xZeUb9qEMqITnd96E7mvlpI1ayhavtzhOt9RoySKEPqMjsQvxJsLR3JJ3Ou6LpBPgAqjlL1pAavVQl5aKtrgELZ9+B4XTx0HYO/qT1nw579L2oyiOfrUUgzp9hv1NpOV8u2Z+A7tJBJ6AyKZXw8Klf2PuV65AbkSLI2U7k3d6Z642rGSb74h7y9/qT3WnTxF9AcrKFq2zOE6VY8eRL7xRsOnu1XnXoF07hVIeaGOy0nO5Ruqy4xsWX6GBa+MJCjC/dMD5QX5fPv6S5Tm5SCTyRxq7hh1On5a8W9Co7vRc/hIeo1s++p+16p8u2NTZJvegj6tFJ+BYRJF5JnEDdDrQeUL45+pO5Z7OY/K6+skccmBdqBi23aHY3NuLpW7dztdp+nfH69O4e4Kq0m9mlgDbbXayDgrzZK6Qz+soTQvB8Bl8bS89FQSd29n/dI3SNjxk7vDa5LNYsWU63yzUxksffVETyOS+fUy+Y/w8C6Ys8zeOm7gPNfXyeQwXdqRZHugulqRrpZCgd/Uqci1jnPPflNvcGNUTQuNbnpeXIpROUBFYct7Apzd+XMbRtIK1qt/6pGpFag86B6EpxDJ/HrqPASGLITAaPjtOzBpiXNDZ5vVI5o6e7rQRx5GfbUehczLi/Bnn0XdsycxK/+L76SJaAYNIuLVV/F3Y/fzxhh1Zn5ZlcSW5Wfs1RGvLlrR+CpBBjK5jP4To4iRqNNQrIupE5nc9Vtfo/WsJCnzkuM7vJPDOf8boiWKxrOJOfO24uUNU14AYyUcfL/uvLYThLq3N2B7pAgJQTt5EubCApRh4bWbh7wHDSLmgw8kjs7R3m9TOH8gp/Y4LMaPmYsHoQ1SU11uRCYDbz/pCsANmjqdPV9+grG6brpCLldgsToOeeVKJaNvv9Pd4TUrcE4sqmg/jNmVqHsG4jPQuYetIEbmbW/KSzDkbnsHoqhhsOArkLB4VXtR8vXXFC1fgaWgEMO5c2Q99hiWcs/cDn/5XLHDcUFmBeveOcH5gzn4+KskTeQAMpmM6Y8+iVJlL1Kl9vElMMJ5ft83IJDOvfu6OzyXrEYLZVsvkr/8NGVbLiLTKDAX66k+kYfhYtMlKP5XiZF5W1P5wJz/2P8ILVbdoPu5tboa3enTaCdMkCiixoVEaakqddw4VJavY8dnSQR39iW8q79EkdXpPWoc0f0HUZSVSXlBPlveX+p0jbe/5+xQLV2fRvWxPACMGeWwr+4xfUoJEc8M96hmFZ5AjMzbQuYhe4OKggtSR9JuaRps4UepRN27tzTBNGPCHb0a3RSUdd5FpymJeGv96NK3P6d/dq7eKVcqGXfH3S6e5X7mEj3Vp/KbuMCGPkkU22pIJPPr7Zc3YOV02Px/sGyUfdencM2C77sPv2nTQCYDuRyFVkvF1q1Sh+VSYLgPC/7fSKbe57zktLkVLlKQO/WrlHH33/5FjyEjJImnPkNGOblLj4O56Z45ihCxNLEhkcyvJ5MO9r9bd2yzwp5/ShdPOyb39ib08cftfSGtViylpeT97U0q9+5r/skS6TOqE/E3RqNQylF4yRk6vSsx/Txnq3zh5QzSTx5lyMybHRJ679FjCYvpJl1g9VTsyQKziz0aXnWpyivSFy+Jlnl6MjFnfj1ZLWBtUJrT0nwRJsG1hs0pAKoPH0I7YbwE0TRPJpMxbm4vRt3cA2R4VPu4Xav+y/GNPwDgExDI1AcfY/eq/2LUVXPh0K3EmYkAABIvSURBVH52rFzO1AcWSxwlYGlkRG6y2pd82sCUU0XBslN0emYYco1IYTXEyPx6UmvtK1fqG/WoNLF0AJoB/Z3P9Xc+52mUKoVHJfLygnyOb6qb7qsuK+Xo+rUYdXX14E/9tImSnGwpwnOgHdu58axUL89byo3oEsW8eX3iY+16m/U2dB1nb+Lccwr09Jwdiu2Nz5AhhP3hDxR9+CE2s5nABXfg5wGbhNobXUU5DbtlGFw09qgqKyMoMspdYbmk6R1E+GPx6BKL0CcVYcptvAGJzEuMRetrVTKvqKjgueeeo7KyEpPJxJIlSxgyZMj1jq19kitg0Hz7H+FXC330EUIefACbzYZcJe167fYqvHtPwrp2pyDjYu25sK7dyTzjWLLZpG++ybc7qLr4oerih3ZMZ0p+SMFwsRxVZ19MhTqs5fbidV5RWrw96H6EJ2hVMv/kk08YPXo09913H+np6Tz77LP88MMP1zs2QQDs2/lFW4fWk8lkzH35L5zY/CNl+Xn0GTuR7POJTsn8Wmq4uIPCX0Xovf2x2WwYUktrC24pgzRo+gYjU4qReX2tSub33XcfqqujJIvFglqtvq5BCYJwffn4BzB+wT0Ox8c3rsN2tbqnUqWm+5DhUoXXpNJ1qVQdzrUfKGSE3j9AJHIXmk3m3377LZ995tgO7Y033mDQoEEUFBTw3HPP8eKLLzb7QgaDgaSk5juyCEKzSkpg9deQmQnDhsJtt4HT2mmhOUMX3EvGkf0ovLzoMW4yWfkFkO9Zo3OZzorvkeK6b2YWG7mbzqOb5jm7VT2FzOaqwHELJCcn88wzz/D8888zadKkZq9PSkoi7moVPEH4NS7ecQf60wm1xyGPPkL4H/4gYURCWzGXGsh903GJqqqbP+GPDpYoIvdrae5s1XeV1NRUnnrqKZYuXdqiRP4/zWqF9F2QtBFM+mYvFxpXdegwBSs+cEjkABU/b5MoIqGtKQPVaPrXu9EpA+24ztIF5MFaNWe+dOlSjEYjf/3rXwHQarUsb9CXUcC+iWjVrXDxaoecwK7w0A7QinZX1yrvb29S3GC6r4YqJsbN0QjuFHJnX6pP5mMu1KHpF4LaAwqXeaJWJXORuFso7Ze6RA5QmgHHP4FJz0sXUztkKS2l+MsvHU/K5WC1ouwcSfj/PStNYIJbyJRyfJtoySfYiU1DbUlf6uKcqMV8rWwWC1gsDufUffvS+Y2/oo6NRaYUv8aCINb3tKXe08Evsu5YoYbBC6SLp51ShoTgP2uWw7mQ++9D07evSOSCcJV4J7QltZ99jvzox/b2cUPuhoiBUkfVLnX+2xv4jhuHMS0V7aRJ+IyQvlyrIHgSkczbWkAU3PgnqaNo92ReXgTeeovUYQiCxxLTLIIgCB2ASOaCIAgdgEjmgiAIHYBI5oIgCB2ASOaCIAgdgEjmgiAIHYBI5oIgCB2ASOaCIAgdgNs2DYnmFIIgCNfOYDC06LpWN6cQBEEQPIeYZhEEQegARDIXBEHoAEQyFwRB6ABEMhcEQegARDIXBEHoANo8mR8+fJgxY8awaNEiFi1axPz581m1atV1+dnjxo1r8vE1a9ZgMplISkri/fffb/XrZGVlMX/+/FY/v63VxLdo0SLS0tJ47733WL16tcN/97Zt28jLy5M4Uvvvw9NPP+1w7umnn8ZoNHL58mXmzJnDH//4R5fPzcrKYujQobW/SzV/LA1ayl3PWPv06cPmzZsdzs+ePZslS5Zc8897+umnOXz48PUKz+Xf5T//+U++//57fvjhB+655x7uv/9+7rvvPvbt21d7zYEDB7j33nu58847WbRoEUuWLKGioqJVMdS8Xmu4+vf8Ne/Tmvd7W/rwww8ZP358i5cLupNb1pmPHj2ad955BwCj0ciMGTOYM2cO/v5t22X7gw8+4JZbbiEuLo64uLg2fS1PVP+/+/PPP+fVV1+lU6dOEkflrOZ348SJE4wZM6bJRBkbG3vdBgMt0aNHDzZu3MjMmTMBSE5ORqfTue31W6OiooIvvviCTZs2oVKpyMvLY968eezatYsLFy7wj3/8gxUrVtT+Lnz66ad8/PHHTh8M7nA9/z1r3u9tacOGDcycOZNNmzZx2223telrXSu3dxqqrKxELpdz4cIFli5dikKhQK1W8/rrr2O1WnnqqacICwsjLy+PiRMn8vTTT7NkyRJmzpzJxIkT2bNnD5s3b+bNN9+s/ZlHjhyp/UTX6/W89dZbHDt2jIKCAp5++mnuvfdevv76a9555x3Wr1/PZ599hkqlolu3brz22mts2LCB3bt3o9fryczM5He/+53Lf6hFixbRt29fUlJSqKys5N///jdRUVEsW7aM7du3Y7FYuPPOO1mwYAErV65k06ZNKJVKhg8fznPPPcd7771HRkYGJSUllJWVcdddd/Hzzz9z8eJF3nrrLeL/fzvnHhRl9cbxzyL3lqtcgoQAqVlqBkVgyNHshg1awESBsAlhIRM6XHTkMkqE47JI3DJnWqLAgTVQV2vCQJtoGh1oINMpnMwCgQSjiLgoIbCw+/uD4f2JQinRZZz389fumXPe833Pc57nPOe87+7y5ajVaj7++GMkEgnr168nJiZm3mPd3NzMoUOHCA0N5bvvviM9PZ2qqioOHz58Sx8ZGRkMDg4yODiISqWioKCAn3/+mYGBAdasWUNKSgqdnZ1kZmai1WoxNTWlsLCQqKgoNBoN1tbWVFVVMTIyQlxc3B3pfPLJJzl48CAqlYrR0VFcXV3x9fVFoVAAYG1tjVKp/MNrJCUlsWrVKkJCQpDL5eTk5FBZWYler6enp4eRkRHy8vJYunTpHWmTyWR0dnZy9epVLC0tqampITg4mJ6eHlatWkVjYyMwlXVHRkZy5coVjh07hk6nIykpifb2djQaDfb29vz2228AaLVaXn/9dX788Ud0Oh0pKSkEBATw7LPP4ubmhrGxMUVFRXek80bMzc2ZnJykurqaJ554AldXV+rr6zEwMKC6upqEhIQZi3psbKzw+UYNaWlpZGdnMzY2xuDgIFu3biUwMJBPPvkElUqFra0tWq0WDw+PeWu9mebmZgoKCjAyMiIiIgJ7e3vefPNNTExMhHkwMTFBSkoKer0erVbL7t27aWlpEfz97bffXjA9N2tzdXUlMjKS1NRUwsLCiI6OxsbGhqtXr1JaWkp2dvYtdj158iTvv/++cJ19+/Zha2u74Pr+kWDe1NREdHQ0EokEIyMjXnvtNZRKJTk5OXh5eVFfX8/evXtJS0vjypUrlJWVYWFhgVwu59tvv/3T67e2tpKfn4+joyMlJSWcPHmShIQEVCoVxcXFfP311wAMDAywf/9+PvzwQ6RSKUqlksOHD2Nubs7w8DBlZWV0dnby6quvzrnqent7s2vXLoqLi6mtrWX16tWcPn0ajUbD+Pg4hYWFfP/995w4cYJDhw5haGhIYmIin3/+OQCmpqaUlZVRWlrKqVOnKCkp4dixY9TW1iKVSqmrq6OqqgqJREJsbCyrV6/+y87y+OOP4+XlRXZ2NpcvX561D5jaQcXGxtLd3c3y5csJDw9nbGxMCOZ5eXnEx8ezZs0a6urquHjxIsHBwdTW1vLiiy9SU1Mz723y4sWLiY+Pp729HblcTkREBEqlEk9PTzQaDe+99x7h4eG0tbURHR0ttHv44YfJyMhAoVAgl8tpaGhgw4YNPPTQQwC4uLiQl5fHqVOnhIz0Tlm7di2ffvopYWFhtLS0sHnzZnp6euasb2lpiUql4tq1a2RnZ3P8+HEkEokwpzQaDTY2NiiVSgYGBti4cSO1tbWMjIywZcsWQfvtMO1b03R1dZGUlMSBAweoqKggLi4OrVbL5s2bkcvldHd34+rqKtTduXMner1eCP43avjiiy/YtGkTAQEBnDt3jv379xMYGEh+fr6wgMfHx9/xeN7IzfacnnMajQa9Xs9TTz1FdXU1jo6OVFRUoFKpCAgIwMLCgsLCQtra2hgeHiY8PFzw978LjUZDeHg4Hh4eGBsb88033wBTx25r166lqqpqVrt2dnZSWlqKmZkZWVlZNDQ0EBISsuD6/vFjlml27dolHAH4+/tTWFgITGVC1tbWwFTg7OjomNFuth+sOjo6kpOTg7m5Ob/88gsrVqyYVUdXVxeenp5IpVKh34aGBpYtW4ZMJgPAycmJ8fHxOe9l2tHuvfde+vr66OjowNvbm0WLFmFmZkZmZiYnTpxg2bJlGBkZAeDn50dra+uM9hYWFnh6egJgZWXF2NgYP/zwAz/99JOQKQ0NDXH58uUFzXzm6gPA3d0dmMqEz58/T1NTE1KpVBiPjo4OfHx8AIRjBw8PD7Zt24a/vz92dnbY2dktiM5Lly6xe/duYCqTndY217bc0tKSkJAQDhw4QEFBgVD+yCOPAODj4/On2f1cBAcHk52djYuLC35+frPWuXFeTmttb2/H09MTY2NjYGo+w5QNzp49S0tLCwATExMMDAzMaHu73OxbBQUF9Pb2Mjo6SlZWFjBlt7i4OHx9fXFycqK7uxuZTIaLiwtqtZqxsTHWrVt3i357e3tUKhVHjx5FIpEwMTFBX18fUqkUGxsbAGE+zJeb7dnc3Cz0PzAwgFQqFXYR/v7+FBUVkZqaSmdnJ1u2bMHQ0JCEhIS/pOF2GBoa4vTp0/T396NWqxkeHubgwYPA/8drLrsuXryY9PR07rnnHtrb21m+fPnfovFfe5vFwcGBixcvAnDmzBnc3NyAKSe+fv06k5OTtLS0CM7w66+/AnDhwoVbrpWZmYlSqWTv3r04ODgIjiWRSNDpdEK9JUuWcOnSJUZGRoCp45lpQ0gkknndh4eHBxcuXECn06HVatm0aRPu7u60tLQwMTGBXq/nzJkzt9WPh4cHnp6eVFZWolarCQsL48EHH5yXrpuRSCTo9fo/7GNa2wcffCBkPi+//DKjo6Po9XqWLl3K+fPnAaipqUGtVuPs7IyFhQUlJSW88MILC6IVphwkLy8PtVpNamoqjz322B/W7+rqora2lujoaPLy8oTy6Z3duXPneOCBB+alxcXFhZGREdRq9YyMamJigt9//53x8XHa2tqEcgMDA6FdW1sbo6OjTE5OCv9N5OHhwTPPPINarebdd98lKCgIKyurGW3/CiYmJuzYsYOhoSEA7rvvPmxsbDAyMiIyMhKVSkVvb69Qv6mpaUb7aQ379u0jNDSU/Px8AgIC0Ov1WFtbc+3aNfr7+wGE+bCQTPdvY2PD8PCwoPXLL7/Ezc2N5uZmHBwcKC8vJyEhQTiSutnfF5Kamhqef/55ysvLKSsr48iRIzQ2NtLf3y/4zWx2NTQ05K233qK4uBiFQoGJicmsCelC8I+fmU+jUCjYs2cPer2eRYsWCVmTkZERycnJ9PX1ERQUhEwmIzw8nJ07d3L8+HEh6N9IaGgoERERWFpaYmdnJxjfz8+P+Ph4tm7dCoCtrS2JiYnExMRgYGCAq6srO3bsoLa2dt734eXlxaOPPkpUVBQ6nY6oqChkMhnr1q0Tynx9fQkMDBQWr7mQyWSsXLmSqKgoxsfH8fb2XrAHlj4+PqSlpVFeXv6nfaxcuZLt27dz9uxZzMzMuP/+++nt7SUtLY2srCxUKhWmpqbk5+cDEBERgUKhEL7/GY2NjTOOsWbbCWVnZ5Oeni68qZKTkwPcui0H2LNnD+np6WRmZuLn50dsbCz19fUAnD59ms8++wydTkdubu5tjtatrF+/no8++gh3d3e6uroAiImJYcOGDSxZsgRnZ+db2tja2pKcnExkZCS2traYmZkBEBkZSWZmJhs3bmR4eBi5XL4gQXwaKysrYmJieOmllzA1NWVyclI4HgBIS0sjIyMDrVbL9evXcXZ2prS09JbrBAUFkZOTwzvvvIOTkxMDAwMYGhqSm5vLK6+8gpWVFYaGf18IkUgkKBQKEhMTkUgkWFlZkZubi0QiYdu2bVRUVGBgYCD497S/V1ZWzjs5mwuNRsMbb7whfDczM+Ppp5/m6NGjQtlsdpVKpaxYsYLnnnsOc3NzLC0tZyykC8l/6o+2uru72b59O0eOHPm3pYjcAXV1dbS2tpKcnPxvS5nBjQ/ORUTudv61zFzk7qCoqIivvvrqb3uDQERE5Pb4T2XmIiIiIiLzQ/w5v4iIiMhdgBjMRURERO4CxGAuIiIichcgBnMRERGRuwAxmIuIiIjcBYjBXEREROQu4H+6oq99q60xRAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import seaborn as sns\n", "sns.set_style(\"whitegrid\")\n", "sns.stripplot(data=scalstat, jitter=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do variable selection on the transformed model." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0. , 0.60158926, 0.64849909, 0.69392297, 0.71256902,\n", " 0.7061129 , 0.69932684, 0.69218231])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lmod = smf.ols('LifeExp ~ np.log(Population) + Income + Illiteracy + Murder + HSGrad + Frost + np.log(Area)', data=statedata).fit()\n", "rss = np.empty(len(pcols) + 1)\n", "rss[0] = np.sum((statedata.LifeExp - np.mean(statedata.LifeExp))**2)\n", "selvar = ['Null']\n", "for k in range(1,len(pcols)+1):\n", " RSS = {}\n", " for variables in itertools.combinations(pcols, k):\n", " predictors = statedata.loc[:,list(variables)]\n", " predictors['Intercept'] = 1\n", " res = sm.OLS(statedata.LifeExp, predictors).fit()\n", " RSS[variables] = res.ssr\n", " rss[k] = min(RSS.values())\n", " selvar.append(min(RSS, key=RSS.get))\n", "adjr2 = 1 - (rss/(50 - np.arange(1,9)))/(rss[0]/49)\n", "adjr2\n" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('Population', 'Murder', 'HSGrad', 'Frost')" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selvar[adjr2.argmax()]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Would be worth creating a function for the all subsets regression." ] }, { "cell_type": "code", "execution_count": 31, "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
Wed Sep 26 10:37:15 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|}{Wed Sep 26 10:37:15 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", "Wed Sep 26 10:37:15 2018 BST" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%load_ext version_information\n", "%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels" ] } ], "metadata": { "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": 2 }