Machine Learning

Lab06 Ridge Regression and Lasso/Lab06 Ridge Regression Lasso.ipynb

{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# 6.6: Ridge Regression and the Lasso" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from sklearn.preprocessing import scale \n", "from sklearn.model_selection import train_test_split\n", "from sklearn.linear_model import Ridge, RidgeCV, Lasso, LassoCV\n", "from sklearn.metrics import mean_squared_error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the sklearn package in order to perform ridge regression and the lasso. The main functions in this package \n", "that we care about are Ridge(), which can be used to fit ridge regression models, and Lasso() which will fit lasso models. \n", "They also have cross-validated counterparts: RidgeCV() and LassoCV(). We'll use these a bit later.\n", "\n", "Before proceeding, let's first ensure that the missing values have been removed from the data, as described in the previous lab." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "<class 'pandas.core.frame.DataFrame'>\n", "Int64Index: 263 entries, 1 to 321\n", "Data columns (total 20 columns):\n", " # Column Non-Null Count Dtype \n", "--- ------ -------------- ----- \n", " 0 AtBat 263 non-null int64 \n", " 1 Hits 263 non-null int64 \n", " 2 HmRun 263 non-null int64 \n", " 3 Runs 263 non-null int64 \n", " 4 RBI 263 non-null int64 \n", " 5 Walks 263 non-null int64 \n", " 6 Years 263 non-null int64 \n", " 7 CAtBat 263 non-null int64 \n", " 8 CHits 263 non-null int64 \n", " 9 CHmRun 263 non-null int64 \n", " 10 CRuns 263 non-null int64 \n", " 11 CRBI 263 non-null int64 \n", " 12 CWalks 263 non-null int64 \n", " 13 League 263 non-null object \n", " 14 Division 263 non-null object \n", " 15 PutOuts 263 non-null int64 \n", " 16 Assists 263 non-null int64 \n", " 17 Errors 263 non-null int64 \n", " 18 Salary 263 non-null float64\n", " 19 NewLeague 263 non-null object \n", "dtypes: float64(1), int64(16), object(3)\n", "memory usage: 43.1+ KB\n" ] } ], "source": [ "# Read Hitter.csv and remove the missing values\n", "df = pd.read_csv('Hitters.csv').dropna().drop('Player', axis = 1)\n", "df.info()\n", "dummies = pd.get_dummies(df[['League', 'Division', 'NewLeague']]) \n", "# Convert categorical variable into dummy/indicator variables." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "<div>\n", "<style scoped>\n", " .dataframe tbody tr th:only-of-type {\n", " vertical-align: middle;\n", " }\n", "\n", " .dataframe tbody tr th {\n", " vertical-align: top;\n", " }\n", "\n", " .dataframe thead th {\n", " text-align: right;\n", " }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", " <thead>\n", " <tr style=\"text-align: right;\">\n", " <th></th>\n", " <th>League_A</th>\n", " <th>League_N</th>\n", " <th>Division_E</th>\n", " <th>Division_W</th>\n", " <th>NewLeague_A</th>\n", " <th>NewLeague_N</th>\n", " </tr>\n", " </thead>\n", " <tbody>\n", " <tr>\n", " <th>1</th>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>2</th>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>1</td>\n", " <td>0</td>\n", " </tr>\n", " <tr>\n", " <th>3</th>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>0</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>4</th>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>0</td>\n", " <td>1</td>\n", " </tr>\n", " <tr>\n", " <th>5</th>\n", " <td>1</td>\n", " <td>0</td>\n", " <td>0</td>\n", " <td>1</td>\n", " <td>1</td>\n", " <td>0</td>\n", " </tr>\n", " </tbody>\n", "</table>\n", "</div>" ], "text/plain": [ " League_A League_N Division_E Division_W NewLeague_A NewLeague_N\n", "1 0 1 0 1 0 1\n", "2 1 0 0 1 1 0\n", "3 0 1 1 0 0 1\n", "4 0 1 1 0 0 1\n", "5 1 0 0 1 1 0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dummies.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now perform ridge regression and the lasso in order to predict Salary on the Hitters data. Let's set up our data:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "<class 'pandas.core.frame.DataFrame'>\n", "Int64Index: 263 entries, 1 to 321\n", "Data columns (total 19 columns):\n", " # Column Non-Null Count Dtype \n", "--- ------ -------------- ----- \n", " 0 AtBat 263 non-null float64\n", " 1 Hits 263 non-null float64\n", " 2 HmRun 263 non-null float64\n", " 3 Runs 263 non-null float64\n", " 4 RBI 263 non-null float64\n", " 5 Walks 263 non-null float64\n", " 6 Years 263 non-null float64\n", " 7 CAtBat 263 non-null float64\n", " 8 CHits 263 non-null float64\n", " 9 CHmRun 263 non-null float64\n", " 10 CRuns 263 non-null float64\n", " 11 CRBI 263 non-null float64\n", " 12 CWalks 263 non-null float64\n", " 13 PutOuts 263 non-null float64\n", " 14 Assists 263 non-null float64\n", " 15 Errors 263 non-null float64\n", " 16 League_N 263 non-null uint8 \n", " 17 Division_W 263 non-null uint8 \n", " 18 NewLeague_N 263 non-null uint8 \n", "dtypes: float64(16), uint8(3)\n", "memory usage: 35.7 KB\n" ] } ], "source": [ "y = df.Salary\n", "\n", "# Drop the column with the independent variable (Salary), and columns for which we created dummy variables\n", "X_ = df.drop(['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n", "\n", "# Define the feature set X.\n", "X = pd.concat([X_, dummies[['League_N', 'Division_W', 'NewLeague_N']]], axis = 1)\n", "\n", "X.info()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(263, 20)\n", "(263, 19)\n", "(263,)\n" ] } ], "source": [ "print(df.shape)\n", "print(X.shape)\n", "print(y.shape)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1 475.0\n", "2 480.0\n", "3 500.0\n", "4 91.5\n", "5 750.0\n", "Name: Salary, dtype: float64" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.head()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# 6.6.1 Ridge Regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Ridge() function has an alpha argument (λ, but with a different name!) that is used to tune the model. We'll generate \n", "an array of alpha values ranging from very big to very small, essentially covering the full range of scenarios from the \n", "null model containing only the intercept, to the least squares fit:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5.00000000e+09, 3.78231664e+09, 2.86118383e+09, 2.16438064e+09,\n", " 1.63727458e+09, 1.23853818e+09, 9.36908711e+08, 7.08737081e+08,\n", " 5.36133611e+08, 4.05565415e+08, 3.06795364e+08, 2.32079442e+08,\n", " 1.75559587e+08, 1.32804389e+08, 1.00461650e+08, 7.59955541e+07,\n", " 5.74878498e+07, 4.34874501e+07, 3.28966612e+07, 2.48851178e+07,\n", " 1.88246790e+07, 1.42401793e+07, 1.07721735e+07, 8.14875417e+06,\n", " 6.16423370e+06, 4.66301673e+06, 3.52740116e+06, 2.66834962e+06,\n", " 2.01850863e+06, 1.52692775e+06, 1.15506485e+06, 8.73764200e+05,\n", " 6.60970574e+05, 5.00000000e+05, 3.78231664e+05, 2.86118383e+05,\n", " 2.16438064e+05, 1.63727458e+05, 1.23853818e+05, 9.36908711e+04,\n", " 7.08737081e+04, 5.36133611e+04, 4.05565415e+04, 3.06795364e+04,\n", " 2.32079442e+04, 1.75559587e+04, 1.32804389e+04, 1.00461650e+04,\n", " 7.59955541e+03, 5.74878498e+03, 4.34874501e+03, 3.28966612e+03,\n", " 2.48851178e+03, 1.88246790e+03, 1.42401793e+03, 1.07721735e+03,\n", " 8.14875417e+02, 6.16423370e+02, 4.66301673e+02, 3.52740116e+02,\n", " 2.66834962e+02, 2.01850863e+02, 1.52692775e+02, 1.15506485e+02,\n", " 8.73764200e+01, 6.60970574e+01, 5.00000000e+01, 3.78231664e+01,\n", " 2.86118383e+01, 2.16438064e+01, 1.63727458e+01, 1.23853818e+01,\n", " 9.36908711e+00, 7.08737081e+00, 5.36133611e+00, 4.05565415e+00,\n", " 3.06795364e+00, 2.32079442e+00, 1.75559587e+00, 1.32804389e+00,\n", " 1.00461650e+00, 7.59955541e-01, 5.74878498e-01, 4.34874501e-01,\n", " 3.28966612e-01, 2.48851178e-01, 1.88246790e-01, 1.42401793e-01,\n", " 1.07721735e-01, 8.14875417e-02, 6.16423370e-02, 4.66301673e-02,\n", " 3.52740116e-02, 2.66834962e-02, 2.01850863e-02, 1.52692775e-02,\n", " 1.15506485e-02, 8.73764200e-03, 6.60970574e-03, 5.00000000e-03])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# alphas represent lambda in Ridge Regression. alpha values range very big to very small.\n", "alphas = 10**np.linspace(10,-2,100)*0.5\n", "# numpy.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None, axis=0)[source]¶\n", "# Return evenly spaced numbers over a specified interval.\n", "alphas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Associated with each alpha value is a vector of ridge regression coefficients, which we'll store in a matrix coefs. \n", "In this case, it is a 19×100 matrix, with 19 rows (one for each predictor) and 100 columns (one for each value of alpha). \n", "Remember that we'll want to standardize the variables so that they are on the same scale. To do this, we can use the \n", "(normalize = True) parameter:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(100, 19)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ridge = Ridge(normalize = True) \n", "# normalize - standardize the variables so that they are on the same scale by squeezing data in [0,1]\n", "coefs = []\n", " \n", "for a in alphas:\n", " ridge.set_params(alpha = a)\n", " ridge.fit(X, y) \n", " coefs.append(ridge.coef_) # coefs matrix - 100 rows x 19 columns\n", " \n", "np.shape(coefs)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array([ 2.41792033e-10, 8.77087820e-10, 3.53418911e-09, 1.48322166e-09,\n", " 1.56674826e-09, 1.84396625e-09, 7.54102361e-09, 2.07601888e-11,\n", " 7.64034224e-11, 5.76187223e-10, 1.53282320e-10, 1.58190714e-10,\n", " 1.67364484e-10, 9.68457013e-11, 1.58184041e-11, -7.37555427e-11,\n", " -2.57640933e-09, -3.46788982e-08, -5.11835167e-10]),\n", " array([ 3.19634838e-10, 1.15945848e-09, 4.67199001e-09, 1.96073175e-09,\n", " 2.07114900e-09, 2.43761486e-09, 9.96878940e-09, 2.74437478e-11,\n", " 1.01000828e-10, 7.61685599e-10, 2.02630206e-10, 2.09118814e-10,\n", " 2.21245998e-10, 1.28024318e-10, 2.09109993e-11, -9.75004868e-11,\n", " -3.40586150e-09, -4.58434625e-08, -6.76615963e-10]),\n", " array([ 4.22538444e-10, 1.53273587e-09, 6.17609583e-09, 2.59197198e-09,\n", " 2.73793708e-09, 3.22238339e-09, 1.31781529e-08, 3.62790194e-11,\n", " 1.33517150e-10, 1.00690354e-09, 2.67865207e-10, 2.76442765e-10,\n", " 2.92474189e-10, 1.69240613e-10, 2.76431104e-11, -1.28889906e-10,\n", " -4.50234845e-09, -6.06023595e-08, -8.94446472e-10]),\n", " array([ 5.58570956e-10, 2.02618662e-09, 8.16443521e-09, 3.42643441e-09,\n", " 3.61939169e-09, 4.25980119e-09, 1.74207426e-08, 4.79587286e-11,\n", " 1.76501815e-10, 1.33106721e-09, 3.54102039e-10, 3.65441066e-10,\n", " 3.86633664e-10, 2.23726131e-10, 3.65425651e-11, -1.70384870e-10,\n", " -5.95183965e-09, -8.01127527e-08, -1.18240558e-09]),\n", " array([ 7.38397930e-10, 2.67849946e-09, 1.07929028e-08, 4.52954464e-09,\n", " 4.78462280e-09, 5.63120648e-09, 2.30291964e-08, 6.33986167e-11,\n", " 2.33325012e-10, 1.75959252e-09, 4.68102056e-10, 4.83091582e-10,\n", " 5.11106950e-10, 2.95752778e-10, 4.83071204e-11, -2.25238771e-10,\n", " -7.86798170e-09, -1.05904344e-07, -1.56307055e-09])]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coefs[0:5]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[array([-1.71961785e+00, 6.06622348e+00, 1.27540705e+00, -6.76054919e-01,\n", " -1.21922719e-01, 5.45628756e+00, -9.49631946e+00, -7.26036515e-02,\n", " 1.96162636e-01, 6.09184001e-01, 8.08032861e-01, 4.25005674e-01,\n", " -6.46655431e-01, 2.80320693e-01, 3.06142702e-01, -3.69563350e+00,\n", " 6.16685846e+01, -1.21878443e+02, -2.77610000e+01]),\n", " array([-1.78822160e+00, 6.35234831e+00, 1.69242754e+00, -9.26987185e-01,\n", " -2.43350576e-01, 5.61390256e+00, -8.80100460e+00, -8.59233028e-02,\n", " 1.96979402e-01, 5.46795916e-01, 8.90564305e-01, 4.59823055e-01,\n", " -6.75778781e-01, 2.80934805e-01, 3.17755891e-01, -3.64766963e+00,\n", " 6.19333364e+01, -1.21176136e+02, -2.76369152e+01])]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coefs[98:100]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We expect the coefficient estimates to be much smaller, in terms of l2 norm, when a large value of alpha is used, as \n", "compared to when a small value of alpha is used. Let's plot and find out:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'weights')" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxcd33v/9dnzuwa7YsXWbLlfYtjxyIOhAQnJAFSQgiFQqDsXAMlt6W3vfdCaW/59RZul0vb2+XShiU3bE24CUtugIQkJEDIZsdbJHmVV1nyIslaRttsn98fM3ZkR7LHtkZnZvR5Ph7Hc7aZeX8leT5zvmcTVcUYY4zJhsftAMYYYwqHFQ1jjDFZs6JhjDEma1Y0jDHGZM2KhjHGmKxZ0TDGGJM1r9sBcq2mpkYXLFjgdgxjjCkoL7/8creq1p4/v+iLxoIFC9iyZYvbMYwxpqCIyOGJ5lv3lDHGmKxZ0TDGGJM1KxrGGGOyZkXDGGNM1qxoGGOMyZoVDWOMMVkr+kNuL9fovtMgglPqwxPx4wl7ERG3YxljjKusaEyi75F2EqdGXp3hFbzlAZyKAE5lEF9dGG9tCN+sEpzKgBUUY8yMYEVjEtUfWklyIEYqGic5GCM5ECPZN0qyb4zR3b0Mbzlxdl1P2ItvXimBxlKCy6rw1UcQjxURY0zxsaIxCV9tGF9teNLlqeE48e4R4p1DxDoGiXcMMvDUaQaePIKnxEdwWSUlzbPwN5XbVogxpmhY0bhMnrCPQKOPQGMZMAeA5FCcsX2nGd3dy0hbD8NbT+KtCVHyutmUXDcbT8B+3MaYwmafYlPIKfERXltHeG0dqViSkZZuhjYfp/9nBxn8VQelNzUQ2TAH8dlBa8aYwmRFI0c8foeSa2ZRcs0sYkcH6X/8EP2PHiD662NU3LWY0PIqtyMaY8wls6+808DfUErtJ66i5hNXIUGHnv/TSu9De0mNJtyOZowxl8SKxjQKLq5g1n9cR+nGeQy/fIIT/7CV2NFBt2MZY0zWrGhMM/F6KH9rE7WfvhqAk/+2k+Gdp1xOZYwx2bGi4ZJAYxl196zFXx+h93u7GXjqCKrqdixjjLmgvC0aInJIRF4Rke0isiUzr0pEnhCRfZnHSrdzXgkn4qf2E1cRXlvLwBOH6Xuk3QqHMSav5W3RyLhJVdeqanNm+nPAU6q6BHgqM13QxOeh8r3LiNxQz9DzXVY4jDF5rdAOub0T2JgZvx94BvivboWZKiJC+e1NAER/fQyAincssjPJjTF5J5+LhgI/FxEF/k1V7wVmqWoXgKp2iUjdRE8UkU3AJoDGxsbpyntFzi8cHr9D+duaXE5ljDHnyueicb2qdmYKwxMisjvbJ2YKzL0Azc3NBdPXc6ZwaDzF4C87cCoCRF4/1+1YxhhzVt7u01DVzszjSeCHwLXACRGZA5B5POlewtwQESresYjgiir6HmlnpLXb7UjGGHNWXhYNESkRkdIz48BtQAvwCPDhzGofBn7sTsLcEo9Qdfdy/PNK6fn3PYwdGXA7kjHGAHlaNIBZwLMisgN4CfiJqj4G/BVwq4jsA27NTBclj9+h+sMrccr89HyrjUT/mNuRjDEGKfbDO5ubm3XLli1ux7hs8RNDnPyXHXjrQtR9cg3ic9yOZIyZAUTk5XGnO5yVr1saJsM3q4Sq9y4j3hHl9MP77BwOY4yrrGgUgNCqaspum8/w9lNnz+Mwxhg3WNEoEKU3NRC6qob+nx1ktL3P7TjGmBnKikaBEBEq370Eb22I3u/tJtFnO8aNMdPPikYB8QS8VP/uSjSeoue7u9BEyu1IxpgZxopGgfHVhal8z1LiRwfp+3/tbscxxswwVjQKUPiqGiJvmsfQi8cZ2nLC7TjGmBnEikaBKr9tAYGF5Zz+0X5inVG34xhjZggrGgVKHKHq/ctxwl56vrOL1HDc7UjGmBnAikYBcyJ+qn53Bcn+MXof3IOm7MQ/Y0xuWdEocIHGMiruWMTontMMPHnY7TjGmCJnRaMIlGyYTbh5FoO/OMpIi11K3RiTO1Y0ioCIUPnOxfgbSun9/l7iJ4bcjmSMKVJWNIqEeD1U/+4KxO+h51tttmPcGJMTVjSKiFMeoPp3V5DoG6Pne7vRpO0YN8ZMLSsaRSawoJzKu5Ywtr+P/p8ccDuOMabIeN0OYKZeSfMs4seHiD57DO+sMJENc9yOZIwpEnm5pSEiDSLytIjsEpFWEfmDzPwvisgxEdmeGW53O2u+Kr+9ieCySvp+3M7ovtNuxzHGFIm8LBpAAvgjVV0BXAd8RkRWZpb9vaquzQw/dS9ifhOPUHX3cnx1IXq+s8uOqDLGTIm8LBqq2qWqWzPjg8AuoN7dVIXHE/RS/ZHViN+h+75WkoMxtyMZYwpcXhaN8URkAbAOeDEz6x4R2Ski3xSRSteCFQhvRYCaD68kNRSn+/5WUmNJtyMZYwpYXhcNEYkADwOfVdUB4KvAImAt0AV8ZZLnbRKRLSKy5dSpU9OWN1/555VSdfdy4sei9H5vF5q0mzcZYy5P3hYNEfGRLhjfVdUfAKjqCVVNqmoK+Bpw7UTPVdV7VbVZVZtra2unL3QeC62spuKuxYzuOc3ph/ehaudwGGMuXV4ecisiAnwD2KWqfzdu/hxV7cpM3gW0uJGvUEWunUNqIMbAk0dwygOUv2WB25GMMQUmL4sGcD3wQeAVEdmemfcnwN0ishZQ4BDwSXfiFa7SNzeSHIgx+PRRPCU+St9oxxcYY7KXl0VDVZ8FZIJFdojtFRIRKu5cTGooTv+jB/CEvJSsn+V2LGNMgcjbfRomd8RJn8MRWFzB6Yf3MtLW43YkY0yBsKIxQ4nXQ/UHV+CbG6Hne7sY3W9njRtjLs6KxgzmCXip+ehqvNUheu5vY+xQv9uRjDF5zorGDOeU+Kj9xFU45QG672sldnTQ7UjGmDxmRcPglPqp+Q9X4SnxceqbLcSORd2OZIzJU1Y0DADe8gC1n7gKT8Dh1NdeIdZhWxzGmNeyomHO8lYFqd20Bk/Q4dTXX7GuKmPMa1jRMOfwVgWp/eQaPGEfp77+iu0cN8acw4qGeQ1vZZDaTVfhlPrp/kYLI7t73Y5kjMkTVjTMhLwVQWo/tQZvXZieb7UytO2k25GMMXkgLy8jYnJHVRka2sfAwHbi8dPE430kU2OEgvWEQg2EwwsJhxchIjgRP7X/4Sp6vt3G6Qf3kDw9SulNDaSvJ2mMmYmsaMwQw8MHOdb5IKdO/ZyRkcNn53s8fkT8JJOvHmYbCMymunojtTVvpqrqBmo+sprTD+9l4OeHiR8fovLdS/H4HTeaYYxxmRWNIheP93Hg4D9y7Nh3AaGy8joaGz9BVeX1BAK1eDwhRIR4vI+RkSNEo3vo7nmGEyf+H52dD+D31zBnzruZ84734JuzgP7HDpHoGaX6AyvwVgXdbp4xZppJsd+Mp7m5Wbds2XLpTzzwS/B4oXQ2RGZBIDL14XKsq+th9u77EonEIPVz30vTws8S8Ndk9dxUaoze3t9wrPNBurt/AaSoqb6ZWbyHxA9KEISKdy4mvLbWuquMKUIi8rKqNp8/37Y0JrH1a39GbLAXvyeJ35MgEPATLK8iUDGLYM1cQrMW4qtbDFVNULUQfCG3I5+VSsXZt/9LdHR8m4qKDSxb+udEIssu6TU8ngA1NTdTU3Mzo6NddHY+SMex79Id/wWRty6n4sCtJL8/yuiuOVTcuRinxJej1hhj8oltaUzi/s9+gu6u4xdcxytJQk6csDdOOOgQjpRQUllFSe1cInMWUtK4itK5iyipqsLxTs+HaizWyyst99DX9yKNDR9n0aL/gsczNd8NkslRjp/4MUeP3sfQ0D58VFHWvpHKUzdT/aY1lFw3B3HsgDxjisFkWxpWNC4gmYgTGxlhbHiYseEhxoaijA0NMRIdYKT3JKM9xxjuPs5Ify/Dg1GGhscYHhNSE9w/Khx0iJSVUFpVQ2RWPaWz5xOpriVSVU2kqprSqmr8ofAVtXUs1s3Wre9ndPQoy5d9mTlz7rqi15uMqtLb+2uOHr2Pnt5fgTpETq6jOnorcza8nfDqWYjHuqyMKWRWNKaJJpOMdO0lenAHQx17GOw6RLT7ONG+fqKjMJgIEE34GU2+dsvD5/cRqaigpKqWcGU1JeUVhMsrCJeXEyorJ1Ralh4ipQRKIjjeV7cgYrFetm77ACMjR1l79TeprLz20rMnUiT7xkhGY6SicZJDcTSWQmNJNJEC1fSNdgEcQXweRp0OTvETuuM/IyH9OGPllJ3eQN2c25l97a14S2xnuTGFqGiKhoi8FfhfgAN8XVX/6kLrX27ReOboM/g9fhrLGplTMgfHc4WHmKrCUDd074HufcRP7CXasY/oqU6ip3sZjDkMJfxE436GEn6GUyGGEj5iycm/sfsCAQLhMMHyALPesBVvSZTRA7fhTSzCFwziC4bwBQL4/AF8wSBefyA97QvgHfLiGRQ8fQp9SVL9CXQgPnl+DzB+h3fy3L+blMSJ1m1jcPaLDNXsRJ04nniYcN8KIsNXUeq7mpKSxXgjJThlfpyIH6fMjyfiw4n4EZ91axmTT4piR7iIOMC/ALcCHcBmEXlEVdum+r3+9NGfc3okQcoD6nEoDZVRHiqnMlhOZbiCylAl1eFKSoMRHEdwPJ70oyM4jgePBzzxODo0RGpoBIZHSA6NkBgZRUerSY01k4xdDeE46k3C2CjEYxCLIfE4JBMEE0kCqSTJVIJUKoGSIkkSNEVKUyQ0QSIWZ/a6l3HCg7Q/tYK+Y32gL6GahFQSSOJRpco3i7mh+dQFG6kJzEU8fhQYTY3RH+smGj9NNNHHUKKf0eQQo8kRxlKjJDROUhMkRVEPpERQUVIeQRwP4nhxvH583iC+jgD+402EA8uoqRkkXHWC0cp2orUvcxwg5RCIziNwoh7fcB2+kTp8I9V4Y+U4Wg7BEpywD0+JD2+JH29JAG/YjyfsxRP0IkEvnqCD+B08gfSj+D2Iz7HuMGOmSUEVDeBaYL+qHgAQkQeAO4EpLxqffq6BEV8dgqa3EtDMngoFHUQYBA6n55/dWlOSKKnM+q8uywyqOJw7LePG01/nA4j6X13nnOWvHcpveppQ3QCnn3gzkYMLiGT6j5RXMwFoXOmJJ+kZOAgcGL8E5Mz6ghIECeAAYUAFVBRIoZoghaKpFCRTpOKpzDul0stJEidKH0pfuwLlwDqqKv1UzUkSro4iVb0MV+0mMfe51/7Qkw4SD0IiAHE/OuzAoBdNOJDyoCkP6SouoIKqkA5/5hHkzP4kPVNEJPPvZEVFQNKtT6GoyNmfnopmfutnps/8RF/dytJxL6vj/h1PX40x8TJjcmTFNZ9nxbpL76q+kEIrGvXA0XHTHcCG81cSkU3AJoDGxsbLeqMB3yGcZG/mP/X4rhPJfJBI5oNAMh8049YTAfGcXX/SQSTzHHn1uWfGZdz4Oc97dX7lqmcILdlP947b6D31Zohw3jrnfnDCmdozguoIqqOoDmce09PoWGY6lhmPARfotspC56n0cEbA00RpcC0VFQFCpYIvnMQbjuMJxBF/HLwx8MVQJ4EGYhCKo54kKknUkwBJoZICSaXbRArkTCHlbBEY/+85H/WiE35Yj/+JXaor6eW1bSSTK3293VP+moVWNCb6//Wa/66qei9wL6T3aVzOG70y92UGvcOoeFBxSHocEh4fCY+PmMfPmCfAqPgZFj/D6kdx0lE8McQzinhGEGcs8zgCnjHEGUXkyr9bCsLrIn7eX3maQ8lZHLi6iuZ5oyw6UIX3xBjDGiVWmyRZpYw5I0T7eon29jB0uoeh/r70lsJ5PI6XcFkZwdIygpEagiURAiURAuESAuEw/mAovY8kGMQXCOD1B/D6/Xj9fhyvD8fnw/H68DgOjteLeDx4HAePx3N2XMSTqbWCeLLfh6GqkCKzI14Zt2HHQHSArq4uOru66Onp5lRPN729vSSTyXNew3EcIpEI4XCYUChEKBQiGAwSCATw+/1nB5/Ph8/nw+v14vV6cRzn7ODxeM4+jh9E5Ozj+HHg7LwzgzGFrtCKRgfQMG56HtCZizcaWfRljp8ewed48Hk9+B0h6HUIeD0EfB6CXoeg3yHkcwj7HUJ+h7DPoSTgJRLwEs48RgJeSgIOpUEfIZ8wlBikf6yf/lg//WP9jCZGGUuOMZIYAdIfkIrid/z4HT++lIN/VPDFwBlN4h1JkRg+ynDon0gN1RL6RTPzuzs4Et/FAY292oDMTyVYEjl7WG9N43wilVWUVFRSUllFuLzi7BFa/lA4bz/URCR92APCwMAA7e3tHDhwgEOHDjE4OHh2ncrKSmpra1myZAmVlZVUVFRQXl5OaWkpwWAwb9tnTCEptKKxGVgiIk3AMeB9wPtz8UYfZBtRevDgxaMOnlRmSDjp6aQHiTvjvk2/+q2azLfKJNCfGTSVQlMpUqkUqWSCZCJBMhGHWAxnbIxgLEZsdITY8DCx0RH6hqKMRaMk4rFzcnm8KZbcdQifJ8nRR+bhGxujurqe8qa5lDXNobSyKl0kqmuIVFbhCxT+Ia8DAwO0trbS0tLCsWPHACgpKaGpqYl58+ZRX1/P7Nmz8fnsrHRjcq2gioaqJkTkHuBx0t89v6mqrbl4r8GeU5zu6iSVTJJKJkglU6RSybPTmkql5yUTpFKpS+rU9jgOHseL4/XiDaS7eXz+AL5QCH84TKS6+mz3ULAkQqgsfX6GZ9DDqf5/Zqh8L/N2/DEb3nErkevri/LCgarKoUOHeP7559m7dy8As2fP5pZbbmHx4sXU1dXhuYQuLmPM1Ci48zQu1XSd3HdmKwIUTSmq5+43OLslIpfYn59SRlq7GfxlB6f0UU6s+j/M5SMsu+5zeMLF981aVdm9eze//OUvOX78OOFwmPXr13P11VdTU5PdxRaNMVeuKM7TyGfi8eBM4Tff1FiCoS0niD7XSbJnlETDSU6u/B6VFdezfN2fkD5lpbh0dXXx+OOPc+jQIaqrq7njjjtYs2aNdTsZk0esaOSZ+Mlhhl7sYmjLCXQsiX9+GZG3zKJl6C/wJStYvfrviq5gxONxnnjiCV566SXC4TC/9Vu/xTXXXIPjFFc7jSkGVjQmEesYxCkP4JT6c/5eqZEEI63dDG0+QezwAHiE0FU1lL6xHt+8CG1tf8zwyCGuWfdt/FneD6NQHD9+nIcffphTp06xYcMGNm7cSCiUP5eZN8acy4rGJHof3EPi1AjemhD+BWUEFpThbyjFWxu+4ktWqCrJ02OM7j/NSEsPY+19kFS8tSHKb28ivK7ubLHq7HqI4yd+RFPTH1BZed1UNC1vbN68mccee4xQKMQHP/hBFi1a5HYkY8xFWNGYRNXvLGPsYD9jB/sZaelheMsJAMTv4Jtbgm9WGG9dGG9NCKcskL74Xth7zrkAmlJSIwlSQ3ESp4aJnxohfnyI2MF+kv3pQ2mdqiCR6+sJra7G31B6zvOj0T3s2fNFKio20LTgM9P7A8ghVeWpp57i2WefZcmSJbzzne+kpKTE7VjGmCxY0ZiEv6EUf0MppTfOQ1NKonuE2NFBYh2DxLuGGN7ZjY4kXvtEJ30hPwR0LPmaxZ4yP4H5ZQQWlhNoKsc7a+KT6hKJKK+03IPXW8LqVf9QNPsxkskkjzzyCDt27KC5uZnbb7/dDp01poBY0ciCeARfXRhfXZiS9bOA9LflVDROomeE5GCM5ECM1FAckoom05e78ITSV2Z1Snx4a0J4a0N4ghf/kacPO/0Cw8OHWLfuWwQCdblu4rRIJpM8+OCD7N27l5tuuokbb7zRztI2psBY0bhMIoJT6s/JjvKOY9/hxMlHWbTwj6mqfP2Uv74bVJVHHnmEvXv3cvvtt3PttVN75U1jzPSwfoE809e3hX37vkR19U3Mn/9Jt+NMmaeeeoodO3awceNGKxjGFDArGnlkZKSDna98mmCwnlUrv5K+KmwReOGFF3j22Wdpbm7mTW96k9txjDFXIKtPJRFZJCKBzPhGEfl9EanIbbSZJZGIsnPnJlQTXL3ma/h85W5HmhJHjhzh8ccfZ/ny5dx+++22D8OYApftV9mHgaSILAa+ATQB38tZqhlGNUlr2x8xNLyf1av/iZKShW5HmhLDw8M89NBDVFRU8M53vtOOkjKmCGT7vzilqgngLuAfVPUPgTm5izVzqKbYvftP6e5+kiWLv0B11RvdjjQlVJUf//jHRKNR3vOe9xAMFt+VeI2ZibItGnERuRv4MPBoZp5dRe4KqSp79/0lnV3fZ8GCe2ho+LDbkabMiy++yJ49e7jtttuYO3eu23GMMVMk26LxUeD1wJdU9WDmJkjfyV2smeHAgb+jo+N+Gho+xsKmz7odZ8r09PTwxBNPsHTpUjZseM0t3I0xBSzb8zRuVdXfPzORKRwjOcpU9FST7Nv/Pzh69D7mzn0fSxb/SdHsIFZVfvrTn+L1enn7299eNO0yxqRlu6UxUb/JR6Ywx4yRTA6z85Xf4+jR+2ho+CjLl/1FUX2wtra20t7ezs0330xZWZnbcYwxU+yCWxqZ/RjvB5pE5JFxi0qBnlwEEpG/Be4AYkA78FFV7RORBcAuYE9m1RdU9VO5yJAro6Od7Hzl9xgcbGXp0j+nYd6H3I40pUZHR3nssceYM2cOr3vd69yOY4zJgYt1Tz0HdAE1wFfGzR8EduYo0xPA5zP3A/9r4PPAf80sa1fVtTl635xRVTq7vs++fV8GlDVr/pXamje7HWvKPf3000SjUe6++247vNaYInXBoqGqh4HDpHeCTwtV/fm4yReAd0/Xe+fC8PBB9uz9/+jt/TWVFdexYsVfEwrNczvWlOvu7uall16iubmZ+vp6t+MYY3Ikqx3hIvIu4K+BOkAyg6pqrjutPwY8OG66SUS2AQPAn6rqr3P8/pctGt3LocP/mxMnfoLHE2Dp0i8yr/4DRXNpkPM988wzeL1eNm7c6HYUY0wOZXv01N8Ad6jqrql4UxF5Epg9waIvqOqPM+t8AUgA380s6wIaVbVHRNYDPxKRVao6MMHrbwI2ATQ2Nk5F5KyMjnZx6tTjnDz1OH19L+E4YeY3foKGxo8TKLLbtI53/PhxWlpaeOMb30gkEnE7jjEmh7ItGiemqmAAqOotF1ouIh8G3g68WVU185wxYCwz/rKItANLgS0TvP69wL0Azc3NejkZBwd3AYrjhHCcMCJeVFNAimRymLFYN7FYNyMjRxgcbCUabWN4+CAAJSVLWdj0h8yb9wF8vsrLefuC8vTTTxMIBLj++uvdjmKMybGLHT31rszoFhF5EPgRmQ9uAFX9wVQHEpG3kt7x/SZVHR43vxboVdWkiCwElgAHpvr9z2hp/SzDw/uzWjcYrKc0spI5s99Nbe1tRXPtqGx0dHSwZ88ebrrpJkKhkNtxjDE5drEtjTvGjQ8Dt42bVmDKiwbwz0AAeCJz/sKZQ2tvBP5CRBJAEviUqvbm4P0BWL78L4nHe0kmR0gmh1FNIngQ8eBxQvj9Nfj9NQQDc4rmirSX4xe/+AXhcJjrrrvO7SjGmGlwsaOnPjpdQca95+JJ5j9M+mq706Kyws4zuJiOjg4OHDjArbfeSiAQcDuOMWYaZHv01D9OMLsf2HJmx7WZeZ5//nkCgQDNzc1uRzHGTJNsj/8MAmuBfZlhDVAFfFxE/iFH2Uwe6+vro62tjfXr19tWhjEzSLZHTy0Gbs7cUwMR+Srwc+BW4JUcZTN57MUXXwSw+30bM8Nku6VRD5SMmy4B5qpqknFHU5mZYXR0lK1bt7Jq1SoqKuyuv8bMJJdyct92EXmG9NngNwJfFpES4MkcZTN5atu2bYyNjfH610/b1WWMMXkiq6Khqt8QkZ8C15IuGn+iqp2Zxf85V+FM/kkmk7zwwgvMnz/frjFlzAx0we4pEVmeebyG9D3BjwJHgNmZeWaG2bdvH/39/XZehjEz1MW2NP4T6Ws4fWWCZQrcPOWJTF7bunUrkUiEpUuXuh3FGOOCi53ctynzeNP0xDH5rL+/n3379nH99dfjOI7bcYwxLsjq6CkRCYvIn4rIvZnpJSLy9txGM/lm+/btqCrXXGM9k8bMVNkecnsf6duvviEz3QH8ZU4SmbyUSqXYunUrTU1NVFVVuR3HGOOSbIvGIlX9GyAOoKojpI+iMjNEe3s7/f39rF+/3u0oxhgXZVs0YiISIr3zGxFZhJ3UN6Ns3bqVUCjE8uXL3Y5ijHFRtif3/TnwGNAgIt8Frgc+kqtQJr9Eo1H27NnDhg0b8Hqz/ZMxxhSjbD8BPgT8BHiI9I2P/kBVu3OWyuSVlpYWUqkU69atczuKMcZl2RaN+4A3kr5A4ULSlxT5lar+r5wlM3lj586dzJ49m7q6OrejGGNcltU+DVX9BfAl4M+ArwPNwKdzmMvkie7ubjo7O1mzZo3bUYwxeSDbmzA9RfrKts8DvwZep6oncxnM5IedO3ciIqxevdrtKMaYPJDt0VM7SZ+nsZr0DZhWZ46mmnIi8kUROSYi2zPD7eOWfV5E9ovIHhF5Sy7e37xKVXnllVdoamqirKzM7TjGmDyQ7VVu/xBARCLAR0nv45gN5OqWbX+vqv9z/AwRWQm8D1gFzAWeFJGlmXt6mBzo6Ojg9OnTvOlNb3I7ijEmT2TbPXUPcAOwHjgMfJN0N9V0uhN4QFXHgIMisp/0pdqfn+YcM8bOnTvxer12boYx5qxsj54KAX8HvHzmlq85do+IfAjYAvyRqp4mfffAF8at05GZ9xoison01XlpbGzMcdTilEwmaWlpYfny5QSDQbfjGGPyRLZHT/2tqr44VQVDRJ4UkZYJhjuBrwKLgLVAF69eln2iy5boJHnvVdVmVW2ura2disgzTnt7OyMjI1x11VVuRzHG5BFXTu9V1VuyWU9EvgY8mpnsABrGLZ4HdL7mSWZKtLa2EgwGWbRokdtRjDF5JNujp6aNiMwZN3kX0JIZfwR4n4gERKQJWAK8NN35ZoJEIsHu3btZsWKFXTbEGHOOfPxE+BsRWUu66+kQ8EkAVW0Vke8DbUAC+IwdOZUb7e3tjI2NsXLlSrejGGPyTN4VDVX94M4xhOkAAA8LSURBVAWWfYn0mekmh850TS1cuNDtKMaYPJN33VPGXfF4/GzXlN3S1RhzPisa5hzt7e3EYjFWrVrldhRjTB6yomHO0draSigUoqmpye0oxpg8ZEXDnBWPx9mzZ491TRljJmVFw5y1f/9+65oyxlyQFQ1zVltbG6FQiAULFrgdxRiTp6xoGCB9Qt/evXtZvny5dU0ZYyZlRcMAcODAATuhzxhzUVY0DJDumgoEAnbUlDHmgqxoGJLJJLt372bZsmV2rSljzAVZ0TAcPHiQ0dFR65oyxlyUFQ3Drl278Pv9dhl0Y8xFWdGY4VKpFLt27WLp0qX4fD634xhj8pwVjRnu8OHDDA8Ps2LFCrejGGMKgBWNGa6trQ2v18uSJUvcjmKMKQBWNGawM11TS5Yswe/3ux3HGFMArGjMYEePHiUajdpRU8aYrOXdQfki8iCwLDNZAfSp6loRWQDsAvZklr2gqp+a/oTFo62tDcdxWLp0qdtRjDEFIu+Khqq+98y4iHwF6B+3uF1V105/quJzpmtq8eLFBAIBt+MYYwpE3nZPiYgAvwP8u9tZitGxY8cYGBiwriljzCXJ26IB3ACcUNV94+Y1icg2EfmliNzgVrBi0NbWhsfjYdmyZRdf2RhjMlzpnhKRJ4HZEyz6gqr+ODN+N+duZXQBjaraIyLrgR+JyCpVHZjg9TcBmwAaGxunNnwRUFXa2tpYtGgRwWDQ7TjGmALiStFQ1VsutFxEvMC7gPXjnjMGjGXGXxaRdmApsGWC178XuBegublZpy55cejs7KS/v5+NGze6HcUYU2DytXvqFmC3qnacmSEitSLiZMYXAkuAAy7lK2jWNWWMuVx5d/RUxvt47Q7wG4G/EJEEkAQ+paq9056swKkqLS0tLFq0iHA47HYcY0yBycuioaofmWDew8DD05+muHR0dNDf38/NN9/sdhRjTAHK1+4pkyMtLS04jmNdU8aYy2JFYwZJpVK0traydOlSO2rKGHNZrGjMIIcPHyYajbJ69Wq3oxhjCpQVjRmkpaUFn89nl0E3xlw2KxozRDKZpK2tjWXLltll0I0xl82Kxgxx4MABRkZGrGvKGHNFrGjMEDt37iQQCLB48WK3oxhjCpgVjRlgdHSUXbt2sXr1arzevDw1xxhTIKxozABtbW0kEgnWrrVbkRhjrowVjRlg+/btVFdXM2/ePLejGGMKnBWNItfb28uRI0dYu3Yt6ftaGWPM5bOiUeR27NgBwJo1a1xOYowpBlY0ilgqlWL79u0sXLiQ8vJyt+MYY4qAFY0idvjwYfr7+20HuDFmyljRKGLbtm3D7/ezfPlyt6MYY4qEFY0iNTQ0RGtrK1dffbVdNsQYM2WsaBSpbdu2kUwmed3rXud2FGNMEbGiUYRSqRSbN29mwYIF1NXVuR3HGFNEXCsaIvIeEWkVkZSINJ+37PMisl9E9ojIW8bNf2tm3n4R+dz0py4Me/fupb+/n2uvvdbtKMaYIuPmlkYL8C7gV+NnishK4H3AKuCtwP8WEUdEHOBfgLcBK4G7M+ua82zevJnS0lK7pasxZsq5dvU6Vd0FTHSW8p3AA6o6BhwUkf3Ama/M+1X1QOZ5D2TWbZuexIWhu7ub9vZ2brrpJhzHcTuOMabI5OM+jXrg6Ljpjsy8yea/hohsEpEtIrLl1KlTOQuajzZv3ozH42H9+vVuRzHGFKGcbmmIyJPA7AkWfUFVfzzZ0yaYp0xc4HSiF1DVe4F7AZqbmydcpxgNDQ2xdetWVq9eTSQScTuOMaYI5bRoqOotl/G0DqBh3PQ8oDMzPtl8Azz//PPE43FuuOEGt6MYY4pUPnZPPQK8T0QCItIELAFeAjYDS0SkSUT8pHeWP+JizrwyPDzMSy+9xKpVq6itrXU7jjGmSLm2I1xE7gL+CagFfiIi21X1LaraKiLfJ72DOwF8RlWTmefcAzwOOMA3VbXVpfh558UXXyQWi3HjjTe6HcUYU8TcPHrqh8APJ1n2JeBLE8z/KfDTHEcrOKOjo7zwwgssX76cWbNmuR3HGFPE8rF7ylyiF198kbGxMdvKMMbknBWNAheNRnnuuedYunQpc+fOdTuOMabIWdEocE899RTxeJxbb73V7SjGmBnAikYB6+joYNu2bVx33XV2xJQxZlpY0ShQqVSKn/3sZ0QiEduXYYyZNlY0CtSOHTs4duwYt9xyC8Fg0O04xpgZwopGARoYGOCJJ55g3rx5rFmzxu04xpgZxIpGgUmlUvzwhz8kHo9z55134vHYr9AYM33sE6fAPPvssxw8eJC3ve1ttvPbGDPtrGgUkCNHjvD000+zatUq1q1b53YcY8wMZEWjQPT39/PQQw9RXl7OHXfcMdHNq4wxJuesaBSAoaEhvvWtbzE2NsZ73/teO1rKGOMaKxp5bnR0lG9/+9v09/fz/ve/nzlz5rgdyRgzg7l2lVtzcUNDQzzwwAOcPHmSu+++m/nz57sdyRgzw1nRyFNdXV088MADRKNRfvu3f5slS5a4HckYY6xoTCaVSrlyDoSqsmPHDh599FHC4TAf+9jHqK+vn/YcxhgzESsak7j//vtJJpPMnz+f+fPn09DQQCgUyul7dnV18fOf/5yDBw+yYMEC3v3udxOJRHL6nsYYcymsaExiwYIFHDhwgOeff57f/OY3AMyaNYv58+fT2NhIQ0MD5eXlV/w+qkpnZyebN29m+/bthEIhbr/9dtavX4/jOFf8+sYYM5VEVaf/TUXeA3wRWAFcq6pbMvNvBf4K8AMx4D+r6i8yy54B5gAjmZe5TVVPXuy9mpubdcuWLZedNRaL0dHRwZEjRzh8+DAdHR3E43EASktLqa+vZ/bs2cyePZu6ujrKy8sv+mE/OjpKV1cXR44cYefOnfT09OA4Dhs2bOCGG27I+RaNMcZcjIi8rKrN5893a0ujBXgX8G/nze8G7lDVThFZDTwOjO/Q/8CZAjNd/H4/CxcuZOHChQAkk0mOHz9OR0cHR48epbOzk927d59dX0QoKyujrKwMv9+P3+/HcRzGxsYYGRkhGo3S19d3dv358+fzhje8gZUrV1qxMMbkPVeKhqruAl5zVrOqbhs32QoERSSgqmPTGO+CHMehvr6e+vp6NmzYAMDY2BgnT57k1KlT9PX10dfXx+DgIKOjowwMDJBMJgkGgwSDQerr67nmmmuYM2cOc+fOpaSkxOUWGWNM9vJ5n8ZvA9vOKxj3iUgSeBj4S52kb01ENgGbABobG3MeNBAI0NDQQENDQ87fyxhj3JSzY0pF5EkRaZlguDOL564C/hr45LjZH1DVq4AbMsMHJ3u+qt6rqs2q2mxXgjXGmKmTsy0NVb3lcp4nIvOAHwIfUtX2ca93LPM4KCLfA64FvjUVWY0xxmQnr649JSIVwE+Az6vqb8bN94pITWbcB7yd9M50Y4wx08iVoiEid4lIB/B64Cci8nhm0T3AYuDPRGR7ZqgDAsDjIrIT2A4cA77mRnZjjJnJXDlPYzpd6XkaxhgzE012nkZedU8ZY4zJb1Y0jDHGZM2KhjHGmKwV/T4NETkFHHY7xxWoIX15lWJQLG0plnaAtSVf5UNb5qvqa050K/qiUehEZMtEO6MKUbG0pVjaAdaWfJXPbbHuKWOMMVmzomGMMSZrVjTy371uB5hCxdKWYmkHWFvyVd62xfZpGGOMyZptaRhjjMmaFQ1jjDFZs6JhjDEma1Y0CpiILBSRb4jIQ25nuVSFnP18IrJCRP5VRB4SkU+7nedKiMhGEfl1pj0b3c5zuUTkhkwbvi4iz7md50qIyEoR+b6IfFVE3u12HisaLhGRb4rISRFpOW/+W0Vkj4jsF5HPXeg1VPWAqn48t0mzdyltyrfs57vEtuxS1U8BvwPk3QlZl/i3pkAUCAId0531Qi7xd/LrzO/kUeB+N/JeyCX+Tt4G/JOqfhr40LSHPZ+q2uDCANwIXAO0jJvnAO3AQsAP7ABWAleR/uMfP9SNe95DbrfnUtuUb9mvtC3AO4DngPe7nf0K/9Y8meWzgO+6nX0K/r6+D5S5nf0Kfyd1wL8Afwv8xu3stqXhElX9FdB73uxrgf2a/hYeAx4A7lTVV1T17ecNJ6c99EVcSpumPdwlutS2qOojqvoG4APTm/TiLvFvLZVZfpr0zc/yxqX+TkSkEehX1YHpTXpxl/g7OamqnwE+h/vXo7KikWfqgaPjpjsy8yYkItUi8q/AOhH5fK7DXaYJ21Qg2c83WVs2isg/isi/AT91J9olm6wt78q049vAP7uS7NJc6P/Mx4H7pj3R5Zvsd7JARO4FvkV6a8NVXrcDmHPIBPMmPftSVXuAT+UuzpSYsE0Fkv18k7XlGeCZ6Y1yxSZryw+AH0x3mCsw6f8ZVf3zac5ypSb7nRwCNk1zlknZlkZ+6QAaxk3PAzpdyjJViqlN1pb8UyztgAJpixWN/LIZWCIiTSLiB94HPOJypitVTG2ytuSfYmkHFEhbrGi4RET+HXgeWCYiHSLycVVNAPcAjwO7gO+raqubOS9FMbXJ2pJ/iqUdUNhtsQsWGmOMyZptaRhjjMmaFQ1jjDFZs6JhjDEma1Y0jDHGZM2KhjHGmKxZ0TDGGJM1KxrG5JCIHBKRmitdx5h8YUXDGGNM1qxoGDNFRORHIvKyiLSKyKbzli0Qkd0icr+I7Mzc5S88bpX/KCJbReQVEVmeec61IvKciGzLPC6b1gYZMwErGsZMnY+p6nrSd+/7fRGpPm/5MuBeVV0DDAC/N25Zt6peA3wV+OPMvN3Ajaq6DvhvwJdzmt6YLFjRMGbq/L6I7ABeIH210iXnLT+qqr/JjH8HeOO4ZWcuR/4ysCAzXg7838wtQf8eWJWL0MZcCisaxkwBEdkI3AK8XlWvBraRvs/2eOdf6G389FjmMcmr97n578DTqroauGOC1zNm2lnRMGZqlAOnVXU4s0/iugnWaRSR12fG7waezeI1j2XGPzIlKY25QlY0jJkajwFeEdlJegvhhQnW2QV8OLNOFen9FxfyN8D/EJHfAM5UhjXmctml0Y2ZBiKyAHg009VkTMGyLQ1jjDFZsy0NY4wxWbMtDWOMMVmzomGMMSZrVjSMMcZkzYqGMcaYrFnRMMYYkzUrGsYYY7L2/wNAr2l4eLf16AAAAABJRU5ErkJggg==\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = plt.gca() # gca - get the current axes\n", "ax.plot(alphas, coefs)\n", "ax.set_xscale('log')\n", "plt.axis('tight')\n", "plt.xlabel('alpha')\n", "plt.ylabel('weights')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For pyplot tutorial including gca(), visit: https://matplotlib.org/3.1.1/tutorials/introductory/pyplot.html#sphx-glr-tutorials-introductory-pyplot-py." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# Split data into training and test sets\n", "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=1)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AtBat 0.098658\n", "Hits 0.446094\n", "HmRun 1.412107\n", "Runs 0.660773\n", "RBI 0.843403\n", "Walks 1.008473\n", "Years 2.779882\n", "CAtBat 0.008244\n", "CHits 0.034149\n", "CHmRun 0.268634\n", "CRuns 0.070407\n", "CRBI 0.070060\n", "CWalks 0.082795\n", "PutOuts 0.104747\n", "Assists -0.003739\n", "Errors 0.268363\n", "League_N 4.241051\n", "Division_W -30.768885\n", "NewLeague_N 4.123474\n", "dtype: float64\n", "106216.52238005563\n" ] } ], "source": [ "# fit a ridge regression model on the training set and get MSE, using lambda=4\n", "ridge2 = Ridge(alpha = 4, normalize = True)\n", "ridge2.fit(X_train, y_train) # Fit a ridge regression on the training data\n", "pred2 = ridge2.predict(X_test) # Use this model to predict the test data\n", "print(pd.Series(ridge2.coef_, index = X.columns)) # Print coefficients\n", "print(mean_squared_error(y_test, pred2)) # Calculate the test MSE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The test MSE when alpha = 4 is 106216. Now let's see what happens if we use a huge value of alpha, say 10^10:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AtBat 1.317464e-10\n", "Hits 4.647486e-10\n", "HmRun 2.079865e-09\n", "Runs 7.726175e-10\n", "RBI 9.390640e-10\n", "Walks 9.769219e-10\n", "Years 3.961442e-09\n", "CAtBat 1.060533e-11\n", "CHits 3.993605e-11\n", "CHmRun 2.959428e-10\n", "CRuns 8.245247e-11\n", "CRBI 7.795451e-11\n", "CWalks 9.894387e-11\n", "PutOuts 7.268991e-11\n", "Assists -2.615885e-12\n", "Errors 2.084514e-10\n", "League_N -2.501281e-09\n", "Division_W -1.549951e-08\n", "NewLeague_N -2.023196e-09\n", "dtype: float64\n", "172862.23580379886\n" ] } ], "source": [ "ridge3 = Ridge(alpha = 10**10, normalize = True)\n", "ridge3.fit(X_train, y_train) # Fit a ridge regression on the training data\n", "pred3 = ridge3.predict(X_test) # Use this model to predict the test data\n", "print(pd.Series(ridge3.coef_, index = X.columns)) # Print coefficients\n", "print(mean_squared_error(y_test, pred3)) # Calculate the test MSE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This big penalty shrinks the coefficients to a very large degree, essentially reducing to a model containing just the \n", "intercept. This over-shrinking makes the model more biased, resulting in a higher MSE." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "532.7090148153795\n" ] } ], "source": [ "print(ridge3.intercept_)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AtBat -1.821115\n", "Hits 4.259156\n", "HmRun -4.773401\n", "Runs -0.038760\n", "RBI 3.984578\n", "Walks 3.470126\n", "Years 9.498236\n", "CAtBat -0.605129\n", "CHits 2.174979\n", "CHmRun 2.979306\n", "CRuns 0.266356\n", "CRBI -0.598456\n", "CWalks 0.171383\n", "PutOuts 0.421063\n", "Assists 0.464379\n", "Errors -6.024576\n", "League_N 133.743163\n", "Division_W -113.743875\n", "NewLeague_N -81.927763\n", "dtype: float64\n", "116690.46856660438\n" ] } ], "source": [ "# least squares is simply ridge regression with alpha = 0.\n", "ridge2 = Ridge(alpha = 0, normalize = True)\n", "ridge2.fit(X_train, y_train) # Fit a ridge regression on the training data\n", "pred = ridge2.predict(X_test) # Use this model to predict the test data\n", "print(pd.Series(ridge2.coef_, index = X.columns)) # Print coefficients\n", "print(mean_squared_error(y_test, pred)) # Calculate the test MSE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like we are indeed improving over regular least-squares!\n", "Instead of arbitrarily choosing alpha=4, it would be better to use cross-validation to choose the tuning parameter alpha. \n", "We can do this using the cross-validated ridge regression function, RidgeCV(). By default, the function performs generalized \n", "cross-validation (an efficient form of LOOCV), though this can be changed using the argument cv." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5748784976988678" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ridgecv = RidgeCV(alphas = alphas, scoring = 'neg_mean_squared_error', normalize = True)\n", "ridgecv.fit(X_train, y_train)\n", "ridgecv.alpha_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, we see that the value of alpha that results in the smallest cross-validation error is 0.57. \n", "What is the test MSE associated with this value of alpha?" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "99825.6489629273" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ridge4 = Ridge(alpha = ridgecv.alpha_, normalize = True)\n", "ridge4.fit(X_train, y_train)\n", "mean_squared_error(y_test, ridge4.predict(X_test))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This represents a further improvement over the test MSE that we got using alpha=4. Finally, we refit our ridge regression \n", "model on the full data set, using the value of alpha chosen by cross-validation, and examine the coefficient estimates." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AtBat 0.055838\n", "Hits 0.934879\n", "HmRun 0.369048\n", "Runs 1.092480\n", "RBI 0.878259\n", "Walks 1.717770\n", "Years 0.783515\n", "CAtBat 0.011318\n", "CHits 0.061101\n", "CHmRun 0.428333\n", "CRuns 0.121418\n", "CRBI 0.129351\n", "CWalks 0.041990\n", "PutOuts 0.179957\n", "Assists 0.035737\n", "Errors -1.597699\n", "League_N 24.774519\n", "Division_W -85.948661\n", "NewLeague_N 8.336918\n", "dtype: float64" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ridge4.fit(X, y)\n", "pd.Series(ridge4.coef_, index = X.columns)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, none of the coefficients are exactly zero - ridge regression does not perform variable selection!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# 6.6.2 The Lasso" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We saw that ridge regression with a wise choice of alpha can outperform least squares as well as the null model on the \n", "Hitters data set. We now ask whether the lasso can yield either a more accurate or a more interpretable model than ridge \n", "regression. In order to fit a lasso model, we'll use the Lasso() function; however, this time we'll need to include the \n", "argument max_iter = 10000. Other than that change, we proceed just as we did in fitting a ridge model:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "lasso = Lasso(max_iter = 10000, normalize = True) \n", "# normalize means all data are squeezed in [0,1]\n", "coefs = []\n", "\n", "for a in alphas:\n", " lasso.set_params(alpha=a)\n", " lasso.fit(scale(X_train), y_train) # scale - Standardize a dataset along any axis\n", " coefs.append(lasso.coef_)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is quite questionable whether we have to scale the data set for fitting a model. If we fit with the scaled data set, we need to roll back to the unscaled prediction later. It might be OK without scaling even though there might be different sensitivities with the data set weighting." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'weights')" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3wc9Z3/8ddntmlVVpJVbVmy3HG3QRCwA4EAoYWjBBLiBDjCndO43F1y+YWUu+TClfQEUkgIR0K5gxBC6KFDSLABF8AVFyzZkmWrd+2utnx/f8zKyLZkq8+u9HnymMfufnd29zNao7fm+535jhhjUEoppUbCcroApZRSqU/DRCml1IhpmCillBoxDROllFIjpmGilFJqxDRMlFJKjZjb6QKckp+fb8rLy50uQymlUsrGjRsbjTEFR7dP2jApLy9nw4YNTpehlFIpRUT29deu3VxKKaVGTMNEKaXUiGmYKKWUGjENE6WUUiOmYaKUUmrENEyUUkqNmIbJEDU3ryUYrHG6DKWUSioaJkMQj0fY8c5XeWP9pTQ0vuB0OUoplTQ0TIbAsjycvOJe/P5SNm9ew5493yUejzhdllJKOU7DZIj8/jJOOfn3lJR8gn377+Dtt28kFgs6XZZSSjlKw2QYXC4fJ83/NgtO+i7NLWt5e/MaDRSl1KSmYTIC06ZdxcIF36elZR2bN3+aWCzkdElKKeUIDZMRmjr1ChYu+B7NLWvZsvXzGBN3uiSllBp3GiajYOrUK5k/799panqZqqpfOF2OUkqNOw2TUVJSspriosvYW3krzS3rnC5HKaXGlYbJKBER5s+/hfT0mWzb9s+EexqdLkkppcaNhskocrszWLL4p0SjHWzf9iWMMU6XpJRS40LDZJRlZs5n7pyv0tzyVxoannG6HKWUGhcaJmNg2rRryMyYz+493yEWCztdjlJKjTkNkzFgWW7mzv06oVA11TW/dbocpZQacxomY2TKlFXk559HVdUvCIcbnC5HKaXGlIbJGJo752bi8TB79/7I6VKUUmpMaZiMofT0mUyf/klqDz6k10BRSk1oGiZjrKz0U4gINTX3OF2KUkqNGQ2TMZaWNo2CggupPfgg0Win0+UopdSY0DAZB2WlNxCNdnDw0MNOl6KUUmNCw2QcZGevIBBYTnX13TqrsFJqQtIwGSel068nGKyiqenPTpeilFKjztEwEZG7RKReRLb2aZsiIs+JyO7EbW6iXUTkNhHZIyKbReTkPq+5PrH+bhG53oltOZHCwovweYvYX32X06UopdSoc3rP5LfAhUe13Qy8YIyZC7yQeAxwETA3sawBbgc7fIBvAu8DTgO+2RtAycSyPJSWXk9Ly1oOHXrM6XKUUmpUORomxphXgOajmi8D7k7cvxu4vE/7Pcb2GpAjIlOBC4DnjDHNxpgW4DmODaikUFr6KbKzK9jxztfo7NzldDlKKTVqnN4z6U+RMeYgQOK2MNFeAlT3Wa8m0TZQe9KxLA9LFt+G253Blq2fIxrtcLokpZQaFckYJgORftrMcdqPfQORNSKyQUQ2NDQ4M1+Wz1fE4kW3EQzuZ/uOm4nHo47UoZRSoykZw6Qu0X1F4rY+0V4DlPZZbzpQe5z2Yxhj7jDGVBhjKgoKCka98MHKzX0fc2Z/hYaGp9n05mqCwQOO1aKUUqMhGcPkMaD3iKzrgUf7tF+XOKrrdKAt0Q32DPAhEclNDLx/KNGW1MrKbmTRwh/T2bmTN9ZfQl39U06XpJRSw+b0ocH3A+uA+SJSIyI3At8BzheR3cD5iccATwF7gT3Ar4HPARhjmoFbgPWJ5duJtqRXXPw3vO+0x0lPn8XWrf9AVdXtTpeklFLDIpP1OuUVFRVmw4YNTpcBQDweYceOr3Co7lHKZ3yWWbO+hEh/Q0FKKeUsEdlojKk4ut3tRDHqSJblYeHC72O50qjadzuxWJC5c7+hgaKUShkaJklCxMVJ8/8Tlyud6urf4PHkMnPmTU6XpZRSg6JhkkREhLlzvk5PTxN7K39CVmAx+XlnO12WUkqdUDIezTWpiQgLTvovMjNPYtu2f6a7e5/TJSml1AlpmCQhl8vP0iW/AIQtWz5LLNbtdElKKXVcGiZJyu8vY/Gin9DZtZO9lbc5XY5SSh2XhkkSy8s7i6lTr6K6+rd0de11uhyllBqQhkmSmz37y1iWj927b2GynhOklEp+GiZJzufNZ9asf6Kp+RUaG19wuhyllOqXhkkKmF7ySTIy5rJ7938Si4WdLkcppY6hYZICLMvDvLn/SjC0n0OH/uh0OUopdQwNkxSRm7sSrzef1rb1TpeilFLH0DBJESJCdmAFbW1vOl2KUkodQ8MkhWRnryAY3EdPT5PTpSil1BE0TFJIIPtkANra33K4EqWUOpKGSQoJZC1BxK1dXUqppKNhkkJcrjQyMxfQ1rbJ6VKUUuoIGiYpJjt7Be3tm4nHo06XopRSh2mYpJjswAri8SBdXTudLkUppQ7TMEkx2b2D8DpuopRKIhomKSYtrQSvt0DDRCmVVDRMUox98uJy2to1TJRSyUPDJAXpyYtKqWSjYZKC9ORFpVSy0TBJQXryolIq2WiYpCD75MX5dLRvdroUpZQCNExSViCwjPaOzRgTd7oUpZTSMElVgaylRKMddHdXOV2KUkolb5iISJWIbBGRt0RkQ6Jtiog8JyK7E7e5iXYRkdtEZI+IbBaRk52tfuwFAksBaG9/2+FKlFIqicMk4RxjzHJjTEXi8c3AC8aYucALiccAFwFzE8sa4PZxr3ScZWTMweVKp71Dx02UUs5L9jA52mXA3Yn7dwOX92m/x9heA3JEZKoTBY4XERdZWUto10F4pVQSSOYwMcCzIrJRRNYk2oqMMQcBEreFifYSoLrPa2sSbRNaILCUjo7txOM9TpeilJrk3E4XcByrjDG1IlIIPCci7xxnXemnzRyzkh1KawDKyspGp0oHBQLLMKaHzs53Do+hKKWUE5J2z8QYU5u4rQf+CJwG1PV2XyVu6xOr1wClfV4+Hajt5z3vMMZUGGMqCgoKxrL8cRHI6h2E164upZSzkjJMRCRDRLJ67wMfArYCjwHXJ1a7Hng0cf8x4LrEUV2nA2293WETWVraNDyePD2iSynluGTt5ioC/igiYNf4f8aYp0VkPfCgiNwI7AeuTqz/FHAxsAfoBm4Y/5LHnz2D8DLadM9EKeWwpAwTY8xeYFk/7U3Auf20G+Dz41Ba0gkEltLY9BLRaAdud5bT5SilJqmk7OZSgxcILAMM7e1bnC5FKTWJaZikuPfOhNeuLqWUczRMUpzHk4PPV0xX926nS1FKTWIaJhOA3z+DYHC/02UopSYxDZMJwO8vIxisPvGKSik1RjRMJgC/v5SengZisW6nS1FKTVIaJhOA329PDaN7J0oppyTleSZqaNL9MwAIBveT2dIOHQchpxSySyE9D6S/qcuUUmr0aJhMAH6/PS1ZcMNt8PorRz6ZWQQzPwCzPgDzLoSMfAcqVEpNdBomqc4Y3NWbcccg2LgJVn4BFn8E2mqgrRpqNsDel2DLg+DPhSt/DXPPd7pqpdQEo2GSqkJtsOX3sPG3yKEt+CsKCc45Dc64xX5+2vL31jUGDr4Fj94E/3s1nPVlOPtmsFzO1K6UmnA0TFKJMXBgE2y8C7Y+DJFuKFoCl/wQv3cDHQOduCgC01bAjc/Bk1+CV74HhzbDR+8Ft3d8t0EpNSFpmDgpFoGeLogE7WAId7y3hFoh2ALdzdBZZy/NldC0GzwZsOQqOPlvoeRkEMG/p5mG5pcxJobIAHsc3nS4/Bd2sPzpy/DYP8AVv9QBeqXUiGmYDNFXd9UQiRsWZflZlJHGoiw/Ga4TdBeF2mDfOtj3Vzj4NnQcgo46CLed+APFgoxCyCqCvNlw+mdhydWQFjhiNb+/DGMihEKH8PuPc8ViEXjfGruml/4DAlPhvG+duA6llDoODZMhqgtHWNvayX0HmwBwCyzLSuf0nEzOyvazkiY8zXuh4R04tMXuTmrcBSYOLi8UL4GCk2DW2ZBRAN4M8KTbiy8TfFngzbQHy/254AuAdeLTgQ4f0RXcd/ww6XXWv0D7AfjrjyFrmh0wSik1TIMKExGZDdQYY8IicjawFLjHGNM6lsUlo7v872LCB6jtamdrGDZG0lgXLOCO1un83HKTE2nnQ00buazhJc6JHsAqXgKLroAZq2B6BXj8Y1KX3z+DSMzN64//geW+3zE1LxvcaXDSJTD1mEvD2Hsol/zQ7j57+mYoPdXu/lJKqWEQ+7pSJ1hJ5C2gAigHnsG+TO58Y8zFY1rdGKqoqDAbNmwY+gt//j57rwPssYvMQsieTnegjFeyl/OEbx7PxrJpjwvlfi+fKsnnmql5BNwuGjvDPL31ELWtQQ62hWgPRijLS2d2QSZzCjNZNC1AVppnWNvz+t4GbrnvPv4Yu4WYuPC5BImFIS0HPv86ZBX3/8JgC/ziDHsvaM3L4PYN6/OVUpODiGw0xlQc0z7IMNlkjDlZRL4MhIwxPxWRN40xKfun7LDDpHGP/Qs3I3/AvYyeeJynGtr4n5pG1rd3kWkJ5wbdrHuthvZQFLclFAXSyEpzs6+pm2AkBtg7C7PyM1hWmsPK2fmcNTefwkDagKUcagvxemUTL75Tz6Nv1fL1tPv4O57i7PCPuOq8s/iHZQK/fD/MPAtW/27ggfZdz8L/XQ1nfgnO/TeIRaF2ExzYaI/vdDVAsNXeXm+63Q2Xlt1nyUl0y+WAf4p9X48SU2pCGihMBjtmEhGRjwPXA5cm2ob3J3Sqy59zwlW8lsXlRblcXpTLH3fX8c3Ht/NMfTeevDS++JGF/O38qQS89o8+Hjccag+xs66DLTVtbK5p5c87G3h40wEATirO4oJFxVy6bBpzCjOpaw/xh001PLzpAHvqOwHI8rn5zKoSrt/6Am3ZUzgpaxm//PO7XHPaORSc9y14+ivw5r1w8nX9FzzvQ7DwCvjLj2D/a/ZYT7jdfs7y2HtfaTkQDdlHn/V0QU/HCX4IWZBbDvlzIX8ezDkXpp+qR44pNUENds9kIfAZYJ0x5n4RmQl8zBjznbEucKwMe89kkDbtb+HnL+7hhXfqyfZ7uOoDM3kpPcr2rjA+S3h/ThZn5mbid1n0/noVAUEwxlDX2M2e/W3srGqh6kAHBsjPTqO5PUTcwGnlU/jQoiLeNzOPhdMCuLY8CH9cw+ZlxWS+fyPn//gVVp9Wxi1/sxDu+RuofRM+uxZyZ7xXZGe9HTJbH4a6rXab5YFl19i//Gessg8S6C8A4jE7cIKt7x3G3Hsoc7AFuhqhpdI++KBlH2Agfz6s+CQs/wRk5I3Zz14pNXZG2s31j8aYW0/UlkrGIkxau3t4fPNBHt5Uw5v7W8lJ93Djqplct7KcbL+HmDGsa+3k2cZ2nm5sY3+oZ3BvHIrhqgtiNYbIzE3jP86ex0dmFx65zq8/SKRzP68sjXPWWZv49pM1/N8b+3n2n89itrsJbl9pH002bTkULoT2WtjxOMQjUHYGzL/Y7r770/+Dy34BKz4xej+YUDtsfwQ23Qs1b4DbDydfC2fcdGS4KaWS3qiMmRzVNinHTILbtmFCIQCMMRxoC/P2/hbeqmllS00bkTiUTUnngwsKOX9RMf5Ed5a923F4HwQEWqMx4on9EoPB2Hcwxr5vgHg8fviz6sNRfrq/jgPhCFcU5XJDaSE+twsadsDDf093xZVssV5k8ZKfE/XM4aO/foOVs/P49mWLMdXrke2PQssupOVdLJ8XOWU1cuqnoGAeiQ+BX51lnzR50wZwjcGR4/U7YO3PYPPv7MOll6+2x2kyC0/8WqWU44YVJolxktXA+4G/9HkqC4gZY84b7ULHy3DDZMv5F+Gurhr9ghwi6elYfv97i9WD1boTa8ZyrKnzkbQ0rDQf4vfjyszEyszEyszClZODKycHd24OroICLO8QB9zbDsDan8L6O+2B/Q98BU5bowP3SiW54Q7ArwUOAvnAD/u0dwCbR6+81PHjpVfRXNpOcVYa5blpzC3MZElJDsW5fixL7L/ugcQuRuKuwfS532cX5Mg3P2JsQo692/u8MWxp7+bOmnqK2iq5ee8dMPtDxJdcyY5tX6Yg/0IK8s7lQGMnt76wm8tXlLBydj6YOCYaxUSimHCIeHeQeHe3vYSCmGCQeFc38bYDRPfsJFbdhgmHMaEQ8VAIEntJ/XHl5uIuKsIzvQRv2Qy8ZWX45s7BN28erqysY1+QXQIXfQdO/Tv7PJdnvw5vPwBX/8YetFdKpZRBdXNNRMPdM/nxv/2VeF0XPuJ2Nw0xjIkBMSxXDI83jstjcHviuD0xXO4YLlcUscJAD0IEl8dlL24XIgJi2beAiAUiSG+bWHYbdjeZEMXEw/QEu2nvCXJW+F5iItw09Wuc+/KTzPvIW3QezKRx02IycqZwu2sVXrfF3y81HMzOZ3ogk/LcHOblTyE7KwuXu5+D8nY8Ab/7BFzxK3swnkTXWzBIrLOTeEcHsdZWYq2tRJubiTY0EK2vJ3qojp6aaiL7qzE9740HeUpK8K9YQcbKlWSccTqeqVOP/cx3nrRnNY6G4dKfwNKPDvm7UUqNvREdGiwiVwLfBQqx/04WwBhjAsd94QSU2XwP7a21hAd4PjTgKwXEB+JJ7JnE4fDISO9jEnsr5qjn+wa+BeLDcvn4UPFWpmbW8d/5X+D1GcsovLyAZe678cx+iywW0dXUzZLW/TwTO4kvu/OJu3zQBXR14nu3geXb3uB9OzeS5xL8mVn4AwH8Wdlk5k4hs6eCjN/firc1A29BOd40P560NDxpaXinFpM2sxxrgCnsTTxO9NAhQrt2EX5nJ6EdO+hat472J54AIG3xYnI+ejWBiy/BlZlhv+ikS2DqcnjoU/Dw39vnuFz4HT2UWKkUMdgB+D3ApcaYHWNf0vgY7p7Jm//2Eu5IYtC8T/eT6Y1YywJL7MVlIS4L8bjAZRFDiGGII8R6X++2ELcFXgsrzYX4PLjTXXj8Hjw+C8tlHR67NyZOqCtGV2sPvurnWN7w/9jUdSVrO67lTxUZbJzt46Z4kDPkk5TP+BLFMz7N59+q5KU/7KRswRRuPbecA83N7Gvr4OVQnFctP24MZ7Qd4vSD71J8cD/BtlY6W5qJ9gwUl+/x+v140zPsIMoK4M8KkJE7xQ6j3ClkFxWTUzwNf5b9N0d41266/vpX2h55hPDu3Vjp6eRcfTV5n/k07txc+01jUXjma/DGr+wwOf2zQ/6OlFJjZ6RHc71qjFk1JpU5ZNgD8P/xGlYwiiR+bAKJ+wYxiTxJtNmdUwYLsIb4F3bEGCIGIkDUZRHzWJg0N1bAizvbkLbrZwRyOsj81F20NUWp3NrIN2JtvJNjsbxnG42ebA4ynZgllG9opqc7yvOfez+5hemHP2Nvd5jbq+t56FAzwbhhTrqPjxVP4eqiXHLjEbpqdtLz6i/p2fEsPXE3EXcmEfz0iI8e0ggbL+G4h5DxEYy46A7H6OoKEQkfecizLyOD/NIZFM2cQ9GsORTPmUfaoQZaHrif9ieexMrIIP/Ta8i99losn88em/ndJ2HX03Dtw/akmEqppDDco7muTNz9AFAMPALv9fAYYx4e5TqHTUQuBG4FXMCdJzqhcrhh8sxt/4tp6SEQyyDbBAhIJi6PG3FbiCexeF2Iz4XlS9z63UiaC/G6wGNhLCEuEEeIxw3xaJxYd5RY8L0l3h3FBCMQjiHhGK5IHE/c9B2WJ2IMbTFDt8siPiUN18wAv3WFeDsjRI53K1O7pzFvTz6HKlt51N/DVZ1e3jc9l7kVhcw9tYiMbHsero5ojMcbWvndwWZeb+vCJXB+XoCri6dQEcigqO1d2PyAfb5ITxf0dNpLuMM+abG9FqLBw3X1ePPoyFlKa8ZJtLqn0dJtUb+/ioaqysN7POnZOZTMX0h54TQCL/+V4Kuv4i0vp+TWn5A2f779WXeeB1319pxhueVD/q6UUqNvuGHym+O8pzHGfGo0ihspsa8GtQs4H6gB1gMfN8ZsH+g1ww2TD99/Mft6qg8/duMmX3LJlykUmXymxYsoiRUyPVzEjOBUvCE3JhTFRAY+Egq3hSvDjZXlxZXlxRXw4spNwz0lsRSkY/lcmLgh9tT3Cb32HF1z/4lQTxmxhm7cHT1Yia+xI2ZojRmCrm7igTaK5q0gLT+NL7+4k4wMDxfF/XTUdCNui8UfKOGUC2fgz3rvcNw93SHuP9jMg4eaaeiJAlDodbMgw0+2x0WGyyLDZZHlcpHtdhHwuJjqdTPNBCnpPkBm/ebEvF6boG4bxhi6XSW0TTmH1vTl1IWm0NjcRltDJV3NlUR7WhFxM50pLNi1FVckTPpnvsS0G1fj7doHvz4H0vPta9dPP2XI35dSanSNqJsr2YnIGcC3jDEXJB5/FcAY898DvWa4YbK7ZTcNwQZaQ620hFto6G6grruOuu46ajtrOdh1kLixg8MSi/JAOQvzFnJy/gpOCSynRKZhglFMMEqsO0K8O0q8O0K8M0KsM0K8vYdYe5h4d/SIz3VNSSM9azPZdV8kNvcarI/dbo+1ACZuiNR20r2zhc53monVd+MKx467HQaIGUMcEK8LV6YHb44Pb04arhwfBLxUeg073IYNEmFbT4TOWIyuWJzOaIyOWP/hWORyMTPuoqTLkFfXTcbeTtKbzOGwgzh+q4M0XwxfhoeI6aK9ZRddrdtxh9pZXt1KfkcbzXkLcJ16FrPOLaJg17eRzkP2GfPnfG3MpvFXSp3YSMdMbuunuQ3YYIx5dBTqGxERuQq40Bjzd4nH1wLvM8bcNNBrhhsmjXf8mlhLC+L1Il4PVlpa4sQ+P5Y/jWi6lzpPkEppZk/sILuD+9nSsp2mkH0xrQxPBm5xYUXDSHTgY7/6HsBlJcZjPCbOgrCbJe0f4+TQMorc+YfXiVthOjLeoTX9bYK+aqK+JshoQ1y9b2QfISCH7/e2ynu3IvT+d2w5I/mjw9gHHoghfszRaUqp8Zbt+T5nnDu8K4iMdNbgNOAk4PeJxx8BtgE3isg5xph/GlZVo6e/0e1jfmOJyBpgDUBZWdmwPqjjhecJ79ptn0cRG/iv/5mJ5XzAeNzUl2ezY7aPg8VdeP11YEWJRDMwcStxNLCxz1vpPcExbj/unVolLkJIoDLTQ23e73mT3zPTuJnmNQQye/DlBpHEBRl7Ot1EOryYej/GWIStGB0SIyr2mItlLMS4seIuLONKPO4NEftH6RYXHvHgETcePLjF3TsLzKDEgbAVJWhF6bESU+wDnrgLlxn8wQjS51vUg4SVGh0zluWfeKUhGmyYzAE+aIyJAojI7cCz2L8rt4x6VUNXA5T2eTwdqD16JWPMHcAdYO+ZDOeDGv/9JrqC7UTjUaKRMLFQiFjQPoM8FuwmFuzChLowwW57CQWxwiHoCTHP3cipoSChTg/NXQF6jBvjAeNJnH7iBfEaxGuwvPHEbQyXL4YrLYY7LYrl7jqinp5ON53tHqp3B+huSiO7yk0ww8c7FVOoL4R97fvojobIimSxtGkpEYmwdcpWivOLObv0bM4pPYelBUux5MhLA//ouV3c9sJuzplfwJcvOImF045/SlE8HqehoYG9e/fy7rvvUlVVRTQaJTs7m6VLlzJ//nymTp2Ky9X/uSlKqdQ22DApATKwu7ZI3J9mjImJyIlPSBh764G5ianxDwDXYM8pNupqt36G9EAYBNxicCfORcQHZJ/49fYPMEYGPWQMsE4sKkQjFpGIRU/UIhQRgm0egk1e2mJCcxwaDey3oMUTx225kUwhkt6Jq8yFJV1E4nXQDB7Lw8Lshcx/ez6Z5Zm8/7z3syh/Edm+gYvtCEX47auVXLCoiF9de+TebDwep7Ozk6amJpqammhsbKS2tpaDBw8SiUQAyM/P55RTTmHBggWUlZVhDeIa9kqp1DbYMPke8JaIvIzd23AW8F8ikgE8P0a1DZoxJioiN2FfUtgF3GWM2TYWnxWKVxBs6QBjn0ViT4NiHV5E3AguEBcWHsCNBNug7SCSlgd58xDxIeIF8QIejPFgjAtjfBjjBqwj5vLyWnHEHcFrRclBKE+MbbgsF26XG5flwuVy0RHtYFfrLowYCvwFFGYWkufPo7GhkV1mFwvzFmKqDDuqd+ByuXC73bjdbizL/rx4PE48HufF7QcpjdRyTiDOk08+SVdXF11dXbS1tdHe3n54JmMAt9tNcXExK1asYNq0aZSXl5OTkzMWP3qlVBIb9NFcIjIVOA07TN4wxhzTjZRKhjsA/8ADD9DS0kLfn1t/903v2Ee4A9PVgHGnQ2bBEXM7GmP6zMl17O3R9+WoEx/7BkAsFju89H52PG6PucSOM7ZzIn6/n/T0dDIyMsjOziYQCJCdnU1eXh55eXkEAgHd81BqEhnWALyInGSMeUdEeq9l0nuCRbGIFBtjNo12ocnummuuGdyK8Ri88G149Sf2Ndg//jv7+unjrKurix/84AesXLmSs88++4jgiUajRKNR4vE4IoJlWTy++SDfe/ZdfnXdaayaX6RjHEqpQTlRN9cXsY9++mE/zxngg6Ne0UQQaoeH18CuP0HFp+Ci74Grn9l5x8H27dsxxrBkyRI8nuPXEI3FuWvDZuaVFnDmgqnH7AkppdRAjhsmxpg1idtzxqecCWDXM/Dkl+wpRi7+AZz2946Ws3XrVvLz8ykqKjrhuo9vrqW6Oci/XrJQg0QpNSSD6uwWkXQR+YaI3JF4PFdEPjy2paWYtgPw4HXwfx8Fbwbc8JTjQdLe3s6+fftYvHjxCcOhKxzle0/vZNG0AOctOHHwKKVUX4M9mus3wEZgZeJxDfYJjE+MRVEppeldePVWePt+QOCD/worv5AUl5/dts0+oG3x4sUnXPfnL+3hYFuIn358hX3FSKWUGoLBhslsY8zHEteExxgTlMncDxINw+5n7cvM7nwKLA+suBZW/SPkznC6usO2b99OcXEx+fnHP9u1srGLO/9SyZUnl1BRPmWcqlNKTSSDDZMeEfGTmKJERGbDgBcbnLiMscdDtj4EoTbIKLT3Qk7/HGQlV9dQKBSipqaGM88885ip0A0AABFbSURBVLjrGWP498e34XVb3HzRSeNUnVJqohlsmHwTeBooFZH/BVYBfztWRSUtEehqgHkXwdKrYebZ4Brsj3B8VVVVYYxh5syZx13v+R31vLyzgW9csoDCrLRxqk4pNdEM9jfhdcCTwEPAXuAfjTGNY1ZVMvvoPSlxXfLKykrcbjelpaUDrtMWjPCvj2xlXlEm168sH7/ilFITzlAG4N+PPbHjLOypVV4xxtw6ZpUlqxQIEoC9e/dSVlaG2z3wV3zLE9tp6Azzq2tPwePSs9iVUsM3qDAxxrwoIn8GTgXOAT4DLMK+TK5KMh0dHTQ0NLBs2bIB13nxnToe2ljD58+ZzbJSnUtLKTUygwoTEXkBe6bgdcBfgFONMfVjWZgavsrKSoABx0vauiN89eEtzC/K4gvnzh3P0pRSE9Rg+zY2Az3AYmApsDhxdJdKQpWVlaSlpTF16tRjnovFDV/5w2YaO3v4wdXL8Ll17i2l1MgNtpvrnwFEJBO4AXsMpRj7Kh4qiRhj2Lt3L+Xl5cfM5tt7GPDT2w7xjUsWsGT6IC7AopRSgzDYbq6bgDOBU4B9wF3Y3V2Tzhd/9xZ1HSH8HjfpXhcZPjeBNDeZPjfZ6R5y0r3kpnvIy/BRGPAxJd07rmeUt7S00NbWxqpVq4557mcv7uGedfv49Fmz+LszZ41bTUqpiW+wR3P5gR8BG3sv3TtZWZYQisRp7goS7InSGY7REYoQjsb7Xd9tCUWBNKbn+imbkk5xdho+t4XHZS9et7343BY+t4s0j4XfY4dUb1hl+z2keQbXHbV3714AZs16LyziccP//LWSHz63iytPLuErF+rJiUqp0TXYbq7vj3UhqeIHV/d/hFRPNE57KEJrdw/NXRGaOsPUd4Spaw9R2xqkuiXIn3c1UN8xvIkDvG6L3HQP84qyWFGaw/KyHOYUZDE1J+2Iw3orKyvJysoiLy8PgB0H2/naH7fw5v5Wzl9YxHc/slTn3lJKjbrkPH07BXndFvmZPvIzjz+MZIwhGjdEYnEiUUM4FqMnGiccjROKxAhF7NvunhjdPVE6w1HaghHaghGaOnvYVtvOz17aQzxxxUZLYGq2n8UlAVbNyedg1T4Kp07ngfXVrK9s5tG3a8n2e/jRR5dxxYoSnVpeKTUmNEzGmYjgcYm9N+EFGPpFs7rCUbbVtlPV2EVNSzf7mrvZUNXCS9sOsDqtk4d2dLJt6xZy0z1cc2opX75gPjnpzs9irJSauDRMUlCGz81pM6dw2sz3Zvg1xvDGll386eE3+eiZi/hAxVJm5KXrnohSalxomEwQIoIVbgfg8jMWkp2d4XBFSqnJRCdkmkDq6+vx+XwEAgGnS1FKTTIaJhNIQ0MDBQUF2rWllBp3GiYTSH19PYWFhU6XoZSahDRMJojOzk66u7s1TJRSjtAwmSDq6+1JnAsKChyuRCk1GWmYTBC9YaJ7JkopJ2iYTBANDQ34/X4yMzOdLkUpNQlpmEwQvYPveiSXUsoJSRcmIvItETkgIm8llov7PPdVEdkjIjtF5II+7Rcm2vaIyM3OVO4cYwz19fU6XqKUckyyngH/Y2PMD/o2iMhC4Brsa89PA54XkXmJp38OnA/UAOtF5DFjzPbxLNhJ7e3thMNhHS9RSjkmWcOkP5cBDxhjwkCliOwBTks8t8cYsxdARB5IrDtpwqShoQHQwXellHOSrpsr4SYR2Swid4lIbqKtBKjus05Nom2g9mOIyBoR2SAiG3p/AU8EeliwUsppjoSJiDwvIlv7WS4DbgdmA8uBg8APe1/Wz1uZ47Qf22jMHcaYCmNMxUT6xVtfX09GRgYZGTq5o1LKGY50cxljzhvMeiLya+CJxMMaoLTP09OB2sT9gdonBZ1GRSnltKTr5hKRqX0eXgFsTdx/DLhGRHwiMhOYC7wBrAfmishMEfFiD9I/Np41OykajVJfX09RUZHTpSilJrFkHID/nogsx+6qqgI+DWCM2SYiD2IPrEeBzxtjYgAichPwDOAC7jLGbHOicCfU1tYSjUYpKytzuhSl1CSWdGFijLn2OM/9J/Cf/bQ/BTw1lnUlq8rKSgDKy8udLUQpNaklXTeXGpqqqiqKiopIT093uhSl1CSmYZLCotEo1dXVzJw50+lSlFKTnIZJCqupqSEajWoXl1LKcRomKayyshIRYcaMGU6XopSa5DRMUlhlZSXFxcX4/X6nS1FKTXIaJimqp6eHmpoaHS9RSiUFDZMUVV1dTTwe1/ESpVRS0DBJUVVVVTpeopRKGhomKaqyspJp06bh8/mcLkUppTRMUlE4HKa2tlbHS5RSSUPDJAXt37+feDyuYaKUShoaJimosrISy7IoLS098cpKKTUONExSUGVlJaWlpXi9XqdLUUopQMMk5QSDQQ4ePKhdXEqppKJhkmKqqqoANEyUUklFwyTFVFZW4na7KSkpcboUpZQ6TMMkxVRWVjJjxgzc7qS7rplSahLTMEkhnZ2dNDQ0aBeXUirpaJikEB0vUUolKw2TFFJZWYnP56O4uNjpUpRS6ggaJimksrKS8vJyXC6X06UopdQRNExSRF1dHc3NzdrFpZRKShomKWLdunV4PB6WLl3qdClKKXUMDZMU0NHRwebNm1mxYgXp6elOl6OUUsfQMEkBr7/+OsYYTj/9dKdLUUqpfmmYJLlwOMyGDRtYsGABU6ZMcbocpZTql4ZJknvzzTcJhUKsXLnS6VKUUmpAGiZJLBaLsW7dOsrKypg+fbrT5Sil1IAcCRMRuVpEtolIXEQqjnruqyKyR0R2isgFfdovTLTtEZGb+7TPFJHXRWS3iPxORCbMRT7Wr19PW1sbq1atcroUpZQ6Lqf2TLYCVwKv9G0UkYXANcAi4ELgFyLiEhEX8HPgImAh8PHEugDfBX5sjJkLtAA3js8mjK2Ojg5efPFFZs+ezbx585wuRymljsuRMDHG7DDG7OznqcuAB4wxYWNMJbAHOC2x7DHG7DXG9AAPAJeJiAAfBB5KvP5u4PKx34Kx9+yzzxKLxbj44ouxN1MppZJXso2ZlADVfR7XJNoGas8DWo0x0aPa+yUia0Rkg4hsaGhoGNXCR1NlZSVbtmxh1apV5OXlOV2OUkqd0JhdFENEngf6m5Hw68aYRwd6WT9thv5Dzxxn/X4ZY+4A7gCoqKgYcD0nRaNRnnzySXJycjjzzDOdLkcppQZlzMLEGHPeMF5WA5T2eTwdqE3c76+9EcgREXdi76Tv+inpmWeeobGxkdWrV+PxeJwuRymlBiXZurkeA64REZ+IzATmAm8A64G5iSO3vNiD9I8ZYwzwEnBV4vXXAwPt9SS9TZs2sX79elauXKmD7kqplOLUocFXiEgNcAbwpIg8A2CM2QY8CGwHngY+b4yJJfY6bgKeAXYADybWBfgK8EUR2YM9hvI/47s1o6OmpoYnn3ySWbNmce655zpdjlJKDYnYf9xPPhUVFWbDhg1OlwFAW1sbd955Jy6XizVr1uhkjkqppCUiG40xFUe3J1s316TT3NzMb37zG8LhMNdcc40GiVIqJY3ZALw6sfr6eu655x5isRjXX3+9Xo5XKZWyNEwccuDAAe677z5cLhc33HADhYWFTpeklFLDpmHigMrKSu6//37S09O57rrrdGp5pVTK0zAZZzt27OChhx5iypQpXHvttQQCAadLUkqpEdMwGUc7duzgwQcfpKSkhNWrV+tgu1JqwtAwGSctLS088sgjTJs2jWuvvRafz+d0SUopNWr00OBxEIvF+MMf/gDAVVddpUGilJpwNEzGwZ///Gdqamr48Ic/TG5urtPlKKXUqNMwGWN79+7llVdeYfny5SxZssTpcpRSakxomIyhAwcO8MADD5Cfn89FF13kdDlKKTVmNEzGSF1dHffddx/p6ek64K6UmvA0TMZAU1MT9957L263m+uuu47s7GynS1JKqTGlhwaPsnfffZeHHnoIEeGGG27Qs9uVUpOChskoMcawdu1ann/+eQoKCvjYxz6m129XSk0aGiajoLW1lT/96U/s3LmThQsXctlll+kYiVJqUtEwGYFoNMratWt55ZVXEBEuuOACTj/9dETE6dKUUmpcaZgMQzgcZtOmTbz22mu0tbWxYMECLrjgAnJycpwuTSmlHKFhMgTxeJyXXnqJ9evXEwqFKCsr49JLL2XOnDlOl6aUUo7SMBkCy7Kora1l5syZrFq1iunTpztdklJKJQUNkyFavXo1LpfL6TKUUiqp6EmLQ6RBopRSx9IwUUopNWIaJkoppUZMw0QppdSIaZgopZQaMQ0TpZRSI6ZhopRSasQ0TJRSSo2YGGOcrsERItIA7HO6juPIBxqdLmKU6LYkn4myHaDbMt5mGGMKjm6ctGGS7ERkgzGmwuk6RoNuS/KZKNsBui3JQru5lFJKjZiGiVJKqRHTMEledzhdwCjSbUk+E2U7QLclKeiYiVJKqRHTPROllFIjpmGilFJqxDRMlFJKjZiGSQoSkVki8j8i8pDTtQxHqtffS0QWiMgvReQhEfms0/WMhIicLSJ/SWzP2U7XMxIicmZiO+4UkbVO1zMSIrJQRB4UkdtF5Cqn6zkeDZNxJiJ3iUi9iGw9qv1CEdkpIntE5ObjvYcxZq8x5saxrXRohrJdyVh/ryFuxw5jzGeAjwJJd6LZEP+tGaATSANqxrvWExni9/KXxPfyBHC3E/UezxC/l4uAnxpjPgtcN+7FDoUxRpdxXICzgJOBrX3aXMC7wCzAC7wNLASWYP8P0Xcp7PO6h5zenuFsVzLWP9ztAP4GWAusdrr2Ef5bsxLPFwH/63Tto/Tv60Eg4HTtI/xeCoGfA98HXnW69uMtumcyzowxrwDNRzWfBuwx9l/sPcADwGXGmC3GmA8ftdSPe9GDMJTtGvfihmCo22GMecwYsxL4xPhWemJD/LcWTzzfAvjGscxBGer3IiJlQJsxpn18Kz2xIX4v9caYzwM3k+RzdmmYJIcSoLrP45pEW79EJE9EfgmsEJGvjnVxI9DvdqVQ/b0G2o6zReQ2EfkV8JQzpQ3ZQNtyZWI77gV+5khlQ3e8/29uBH4z7hUN30DfS7mI3AHcg713krTcThegAJB+2gY8m9QY0wR8ZuzKGTX9blcK1d9roO14GXh5fEsZsYG25WHg4fEuZoQG/P/GGPPNca5lpAb6XqqANeNcy7DonklyqAFK+zyeDtQ6VMtomijbNVG2A3RbklXKb4uGSXJYD8wVkZki4gWuAR5zuKbRMFG2a6JsB+i2JKuU3xYNk3EmIvcD64D5IlIjIjcaY6LATcAzwA7gQWPMNifrHKqJsl0TZTtAtyVZTaRt6UsnelRKKTViumeilFJqxDRMlFJKjZiGiVJKqRHTMFFKKTViGiZKKaVGTMNEKaXUiGmYKOUAEakSkfyRrqNUstAwUUopNWIaJkqNMRF5REQ2isg2EVlz1HPlIvKOiNwtIpsTV21M77PKP4jIJhHZIiInJV5zmoisFZE3E7fzx3WDlOqHholSY+9TxphTsK/G+AURyTvq+fnAHcaYpUA78Lk+zzUaY04Gbgf+JdH2DnCWMWYF8G/Af41p9UoNgoaJUmPvCyLyNvAa9sywc496vtoY82ri/n3A+/s81zst/EagPHE/G/h94rKvPwYWjUXRSg2FholSY0hEzgbOA84wxiwD3sS+znpfR0+Q1/dxOHEb473rD90CvGSMWQxc2s/7KTXuNEyUGlvZQIsxpjsx5nF6P+uUicgZifsfB/46iPc8kLj/t6NSpVIjpGGi1Nh6GnCLyGbsPYrX+llnB3B9Yp0p2OMjx/M94L9F5FXANZrFKjVcOgW9Ug4SkXLgiUSXlVIpS/dMlFJKjZjumSillBox3TNRSik1YhomSimlRkzDRCml1IhpmCillBoxDROllFIjpmGilFJqxP4/v6a2LESHDY0AAAAASUVORK5CYII=\n", "text/plain": [ "<Figure size 432x288 with 1 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = plt.gca()\n", "ax.plot(alphas*2, coefs)\n", "ax.set_xscale('log')\n", "plt.axis('tight')\n", "plt.xlabel('alpha')\n", "plt.ylabel('weights')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that in the coefficient plot that depending on the choice of tuning parameter, some of the coefficients are exactly \n", "equal to zero. We now perform 10-fold cross-validation to choose the best alpha, refit the model, and compute the associated \n", "test error:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "104960.65853895503" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lassocv = LassoCV(alphas = None, cv = 10, max_iter = 100000, normalize = True)\n", "lassocv.fit(X_train, y_train)\n", "\n", "lasso.set_params(alpha=lassocv.alpha_)\n", "lasso.fit(X_train, y_train)\n", "mean_squared_error(y_test, lasso.predict(X_test))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is substantially lower than the test set MSE of the null model and of least squares, and only a little worse than \n", "the test MSE of ridge regression with alpha chosen by cross-validation.\n", "\n", "However, the lasso has a substantial advantage over ridge regression in that the resulting coefficient estimates are sparse. \n", "Here we see that 13 of the 19 coefficient estimates are exactly zero:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AtBat 0.000000\n", "Hits 1.082446\n", "HmRun 0.000000\n", "Runs 0.000000\n", "RBI 0.000000\n", "Walks 2.906388\n", "Years 0.000000\n", "CAtBat 0.000000\n", "CHits 0.000000\n", "CHmRun 0.219367\n", "CRuns 0.000000\n", "CRBI 0.513975\n", "CWalks 0.000000\n", "PutOuts 0.368401\n", "Assists -0.000000\n", "Errors -0.000000\n", "League_N 0.000000\n", "Division_W -89.064338\n", "NewLeague_N 0.000000\n", "dtype: float64" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Some of the coefficients are now reduced to exactly zero.\n", "pd.Series(lasso.coef_, index=X.columns)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above data result has been a motive for the Movie Moneyball.\n", "\n", "**Movie Moneyball:** Billy Beane (Brad Pitt), general manager of the Oakland A's, one day has an epiphany: Baseball's conventional wisdom is all wrong. Faced with a tight budget, Beane must reinvent his team by outsmarting the richer ball clubs. Joining forces with Ivy League graduate Peter Brand (Jonah Hill), Beane prepares to challenge old-school traditions. He recruits bargain-bin players whom the scouts have labeled as flawed, but have game-winning potential. Based on the book by Michael Lewis." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now it's time to test out these approaches (ridge regression and the lasso) and evaluation methods (validation set, cross \n", "validation) on other datasets. You may want to work with a team on this portion of the lab. You may use any of the datasets \n", "included in ISLR, or choose one from the UCI machine learning repository (http://archive.ics.uci.edu/ml/datasets.html). \n", "Download a dataset, and try to determine the optimal set of parameters to use to model it! You are free to use the same \n", "dataset you used in Lab 9, or you can choose a new one." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Lab06 Ridge Lasso Homework in a separate file." ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }