{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## DS-GA 1003, Machine Learning Spring 2021\n", "### Lab 3 : 17-Feb-2021, Wednesday\n", "### Prostate Cancer Analysis with LASSO" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this demo, we illustrate the classic technique of **LASSO regularization**. You will learn to:\n", "* Fit a LASSO model using the `sklearn` package\n", "* Determine the regularization level with cross-validation\n", "* Draw the coefficient path as a function of the regularization level" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Background\n", "\n", "We use a classic prostate cancer dataset from the paper:\n", "\n", "> Stamey, Thomas A., et al. \"[Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate. II. Radical prostatectomy treated patients](http://www.sciencedirect.com/science/article/pii/S002253471741175X).\" The Journal of urology 141.5 (1989): 1076-1083.\n", "\n", "In the study, **the level of [prostate specific antigen](https://en.wikipedia.org/wiki/Prostate-specific_antigen)** was measured in 102 men before they had a prostatectomy. Elevated values of the PSA are believed to be associated with the presence of prostate cancer and other disorders. To study this hypothesis, various features of the prostate were measured after the prostatectomy. Data analysis is then used to understand the relation between the PSA level and prostate features. The study is old and much more is known about PSA today. But, the analysis is typical for medical problems and illustrates the basic tools well." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The specific analysis presented in this demo taken from the class text: \n", "\n", "> Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. [Elements of statistical learning](https://www.amazon.com/exec/obidos/ASIN/0387952845/trevorhastie-20), New York: Springer series in statistics, 2001.\n", "\n", "The text provides an excellent discussion of LASSO and other methods on this dataset. \n", "\n", "Special thanks to [Phil Schniter](http://www2.ece.ohio-state.edu/~schniter/) at Ohio State for pointing on error in an earlier version of this demo.\n", "\n", "First, we load the regular 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", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading the Data\n", "\n", "Our analysis begins by getting the data from Tibshirani's website. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Get data\n", "url = 'https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data'\n", "df = pd.read_csv(url, sep='\\t', header=0)\n", "df = df.drop('Unnamed: 0', axis=1) # skip the column of indices" ] }, { "cell_type": "code", "execution_count": 4, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lcavollweightagelbphsvilcpgleasonpgg45lpsatrain
0-0.5798182.76945950-1.3862940-1.38629460-0.430783T
1-0.9942523.31962658-1.3862940-1.38629460-0.162519T
2-0.5108262.69124374-1.3862940-1.386294720-0.162519T
3-1.2039733.28278958-1.3862940-1.38629460-0.162519T
40.7514163.43237362-1.3862940-1.386294600.371564T
\n", "
" ], "text/plain": [ " lcavol lweight age lbph svi lcp gleason pgg45 lpsa \\\n", "0 -0.579818 2.769459 50 -1.386294 0 -1.386294 6 0 -0.430783 \n", "1 -0.994252 3.319626 58 -1.386294 0 -1.386294 6 0 -0.162519 \n", "2 -0.510826 2.691243 74 -1.386294 0 -1.386294 7 20 -0.162519 \n", "3 -1.203973 3.282789 58 -1.386294 0 -1.386294 6 0 -0.162519 \n", "4 0.751416 3.432373 62 -1.386294 0 -1.386294 6 0 0.371564 \n", "\n", " train \n", "0 T \n", "1 T \n", "2 T \n", "3 T \n", "4 T " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this data set, the target variable is `lpsa`, the log of the PSA. The goal is to try to predict the `lpsa` from various prostate features." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Print the names of the target and predictor\n", "names = df.columns.tolist()\n", "names_x = names[0:7]\n", "name_y = names[8]\n", "\n", "\n", "# Convert the dataframe values to data matrices\n", "X0 = np.array(df[names_x])\n", "y0 = np.array(df[name_y])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Data Cleaning\n", "\n", "- When performing any regularized estimate, it is critical to standardize the values. \n", "\n", "- For this purpose, we use `sklearn` built-in `scale` command, which makes each feature / label to have zero mean and unit standard deviation." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import sklearn.preprocessing\n", "\n", "X = sklearn.preprocessing.scale(X0)\n", "y = sklearn.preprocessing.scale(y0)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Adding \"bad\" features\n", "\n", "- noise: a vector randomly sampled from the standard normal distribution\n", "- redundant: exactly same as lcp\n", "- dependent: lcp + gleason" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "random_noise = np.random.randn(X.shape[0])\n", "redundant_feature = X[:,5]\n", "dependent_feature = X[:,5]+X[:,6]\n", "X = np.concatenate([X, random_noise.reshape(-1,1), \n", " redundant_feature.reshape(-1,1),\n", " dependent_feature.reshape(-1,1)], axis=1)\n", "names_x = names_x + [\"noise\", \"redundant\", \"dependent\"]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Target variable: lpsa\n", "Predictors: ['lcavol', 'lweight', 'age', 'lbph', 'svi', 'lcp', 'gleason', 'noise', 'redundant', 'dependent']\n", "\n", "num samples = 97, num features = 10\n" ] } ], "source": [ "# Print the number of samples and features\n", "nsamp = X0.shape[0]\n", "nfeatures = X.shape[1]\n", "print(\"Target variable: %s\" % name_y)\n", "print(\"Predictors: \"+str(names_x))\n", "print(\"\")\n", "print(\"num samples = %d, num features = %d\" % (nsamp, nfeatures))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit a Linear Model with No Regularization\n", "\n", "First, we try to fit a multiple linear model with no regularization. We begin by importing the appropriate package." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from sklearn import linear_model\n", "from sklearn.model_selection import train_test_split" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We next split the data into training and test -- will use roughly half the samples for each." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "num samples train = 48, test = 49\n" ] } ], "source": [ "X_tr, X_ts, y_tr, y_ts = train_test_split(X,y,test_size=0.5,shuffle=True)\n", "ntr = X_tr.shape[0]\n", "nts = X_ts.shape[0]\n", "print(\"num samples train = %d, test = %d\" % (ntr, nts))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the model on the training data." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearRegression()" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "regr = linear_model.LinearRegression()\n", "regr.fit(X_tr,y_tr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can measure the normalized RSS on the training data. " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R^2 training = 0.676953\n" ] } ], "source": [ "# Measure normalized RSS\n", "y_tr_pred = regr.predict(X_tr)\n", "rsq_tr = 1-np.mean((y_tr_pred-y_tr)**2)/(np.std(y_tr)**2)\n", "print(\"R^2 training = %f\" % rsq_tr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ths `R^2` value is about `0.68`. However, we need to evaluate the model on the test data. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Normalized test RSS = 0.387242\n", "Normalized test R^2 = 0.612758\n" ] } ], "source": [ "y_ts_pred = regr.predict(X_ts)\n", "rss_ts = np.mean((y_ts_pred-y_ts)**2)/(np.std(y_ts)**2)\n", "rsq_ts = 1-rss_ts\n", "print(\"Normalized test RSS = %f\" % rss_ts)\n", "print(\"Normalized test R^2 = %f\" % rsq_ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also plot the actual vs. predicted values. We see a clear fit." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhV1dXH8e9KSCAIiAoyicapCIiIYNXSWlAoYh3Q1oq0lcokMihqERDnCQRFZFABUVuKQlEEFRW1EKfiAKIgUgR5VQyIIIZBwpTs949zk5sLmXPvPXf4fZ6H5032Hc56T+Ne5+yz99rmnENERJJPit8BiIiIP5QARESSlBKAiEiSUgIQEUlSSgAiIkmqmt8BVES9evVcZmam32EA8PPPP3PYYYf5HYbvdB48Og8enYfYPAfLli3b6pyrf3B7XCWAzMxMli5d6ncYAGRlZdGhQwe/w/CdzoNH58Gj8xCb58DMvimuXUNAIiJJSglARCRJKQGIiCQpJQARkSSlBCAikqSUAEREkpQSgIhILMvNhREjIAJT4ONqHYCISFJ5913o0we+/BIWLPCSQHp62L5edwAiIrFm504YOBDOPdfr/AFWroSZM8N6GN0BiIjEktdeg2uvhQ0bgm21a8NDD0HPnmE9lBKAiEgs+PFHuPFGmDEjtP2ii+Dxx+GYY8J+SCUAERE/OQfPPw+DBsEPPwTb69WDCROge3cwi8ihlQBERPyyaRMMGADz5oW29+gB48dD/UMKeIaVEoCISLQ5B08/DTfdBNu3B9ubNPGGey6+OCphaBaQiEg0rV8PnTtD796hnf+118KqVVHr/EF3ACIi0ZGXBxMnwsiRsHt3sP3EE2HaNOjYMeohKQGIiETaF194V/wffBBsS0nxhoDuvhtq1vQlLCUAEZFI2bcPHnwQ7rvP+7lAq1YwfTqceaZ/saEEICISGUuXelf9K1YE29LS4LbbYPjwsJZ0qCzfHgKbWVMzW2xmq81slZnd4FcsIiLhkrJnD9xyC5x1Vmjnf9ZZsHw53HFHTHT+4O8dwAHgZufcJ2ZWG1hmZm86577wMSYRkcp7+23a9ekD2dnBtowMuP9+uP56SE31L7Zi+HYH4Jzb5Jz7JPDzTmA10MSveEREKm3HDrjuOujQgZpFO//zzoPPP/dKPMRY5w9gzjm/Y8DMMoF3gFOdczsOeq0f0A+gQYMGbWfNmhX1+Iqza9cuatWq5XcYvtN58Og8eJLxPBy5ZAm/eOQRamzZUth24LDDWHfddXx/4YURK+NQER07dlzmnGt3yAvOOV//AbWAZcDlZb23bdu2LlYsXrzY7xBigs6DR+fBk1TnYcsW5/78Z+e8db2F/7a0b+9cdrbf0YUAlrpi+lRfZwGZWRrwAjDTOTfXz1hERMrFOZg9GwYPhq1bg+3168OkSXxevz4dGjf2L74K8HMWkAHTgdXOuXF+xSEiUm7Z2dCtG1x1VWjn/5e/wOrV8Kc/xcSQT3n5WQuoPfBX4Dwz+zTw70If4xERKZ5zXrmGFi3gpZeC7ccc423VOGMGHHWUf/FVkm9DQM6594D4SZUikpy++gr69oXFi0Pbr7sORo+GOnX8iSsMtBJYRKQ4eXnw6KPeyt3c3GD7ySfDk096+/XGOSUAEZGDff65V8bho4+Cbamp8Pe/w513eou7EoASgIhIgX37YNQob+Xu/v3B9tatveJtbdv6F1sEKAGIiIB3td+7t3f1XyA93avdc8stXiG3BKMEICLJbfduuP12bw/e/Pxg+znneFf9zZv7F1uEKQGISPJavBj69PG2aSxQs6Y3DDRwYEzW7wknJQARST7bt8PQod7c/qI6d4YpU+D44/2JK8qUAETCbN7ybMYuXMPGnFwa181gaJdmdGujQrcx4+WXoX9/2Lgx2Fa3LjzyCPTsGVcreatKCUAkjOYtz2bE3JXk7s8DIDsnlxFzVwIoCfhtyxavJv/BFYUvuwwmT4ZGjfyJy0d+loIQSThjF64p7PwL5O7PY+zCNT5FJDgHM2d6D3OLdv4NGsDzz8PcuUnZ+YPuAETCamNOboXaJcI2bPBKNixYENresyeMGwdHHulPXDFCdwAiYdS4bvErREtqTxTzlmfTfvQiVmZvp/3oRcxbnl32hyIpPx+eeAJatgzt/I89Fl5/HZ55Juk7f1ACEAmroV2akZEWOnUwIy2VoV2a+RRR5BU898gO3OUUPPfwLQmsXettxXjddbBzp9dmBoMGeYu8unTxJ64YpAQgEkbd2jRh1OWtaFI3AwOa1M1g1OWtEvoBcMw89zhwAMaOhdNOg7ffDrY3awbvvAMTJ0Lt2tGNKcbpGYBImHVr0yShO/yDxcRzjxUrvDIOS5cG21JTvRIOd9wBNWpEL5Y4ojsAEakSX5977N3rdfBt24Z2/qefDh9/DA88oM6/FEoAIlIlvj33WLIE2rSBe+/1hn8Aqlf3Ov2PPvJek1JpCEgkyVV15XLBe70x/500ifTq559/hpEjYcIEb45/gfbtvY1aTjklMsdNQEoAIkksXCuXC557ZGVlMfjPHSIRquett7ztGb/+urBpd3oGa28cSesHRkCKBjUqQmdLJInFzAyesuTkeA95O3cO6fzfPv4MOveaTPe0tsz7bJN/8cUp3QGIJLGYmMFTlnnzYMAA2BTs4HNq1OKe8/syt+V53hz/QNJKptlX4aAEIJLEGtfNKFzAdXC77zZvhsGDYc6ckOZXm7Xnjs792XrYESHtMZW04oQSgEiMiWY56aFdmoU8A4AYWLnsHMyYAUOGwE8/BdsbNoTJk7n/y7psjdWkFWf0DEAkhhQtq+CIfFmFmFu5/O23cOGFXrG2op3/NdfAF1/A5ZcnZbmNSNEdgEgMKe2hbKQ65ZhYuZyfD48/DsOHw65dwfbMTJg61Xv4G1B02qk23akaJQCRGBIXD2XDbc0ab1/e994Ltpl5m7fcdx/UqnXIR2IiaSUADQGJxJCkKie9fz+MHg2tW4d2/s2bw/vvw/jxxXb+Ej5KACIxJGnGt5cvh7POghEjvHo+ANWqwW23ea+dc46/8SUJDQGJxJCEH9/es8er3fPgg5BX5FlH27Ywfbp3NyBRowQgEmMSdnz7/fe91bxriqwyrlED7r4bbrrJuwOQqNIZF5HI2rULbr0VJk0KLd527rkwbRr84hf+xZbklABEJHLeeAP69YNvvgm21a4NY8Z47SreVqZILgxUAhCR8Nu2DW6+2dt8vaiuXWHKFGja1Jew4k24qrWWROlXRMLrhRegRYvQzv+oo+Bf/4IFC9T5V0Ckq7X6egdgZk8BFwE/OOdO9TMWEamiTZtoeccd8O67oe1XXult3nL00f7EFccivTDQ7zuAZ4ALfI5BRKrCOe9qv0UL6hft/Bs18ko5z5qlzr+SIr0w0NcE4Jx7B9jmZwwiUgVffw1dunjF2nJygu19+3rF2y691LfQEkGkFwaaKzotywdmlgm8UtIQkJn1A/oBNGjQoO2sWbOiF1wpdu3aRS0tU9d5CEi685CfT5N58zhh2jRS9+wpbP65YUPWDh1Kzhln+Bicv8L9t5CTu5/N2/ewLy+f9NQUGhxeg7oZaRX6jo4dOy5zzrU7uD3mZwE556YCUwHatWvnOnTo4G9AAVlZWcRKLH7SefAk1XlYvdor3vbf/wbbUlJgyBCWderEuV27+hdbDIinvwW/nwGISLzYvx/uvx9OPz2082/Rwvv94YfJz0jAonUJLObvAEQkBnzyCfTqBZ99FmyrVg1GjvQKulWv7l9sUmm+3gGY2XPAEqCZmX1nZr39jEdEQr20ZB0zOlzFgXZnhnb+7dp5SeGuu9T5xzFf7wCcc1f5eXwRKdm70+fS6pbrOX5bcDvK3GrV+WrwLZw65g4Vb0sA+l9QRELt3AnDh/Obxx4LaV5ybCuGXzCYA0efyNCVmxO3ZHUSUQIQkaDXXoNrr4UNGwqbdqTXZFTHXsxq/TucpUCgHk1x9Wnq+hK0VJYSgEiciGRVSH78EW68EWbMCGl+68Qzue13A/m+Tr3CtlSzEuvT3H+2JhbGEyUAkTgQsaqQzsGcOTBoEGzZUtj8Y0Yd7u50LS81P9fboL2IvBIWj3r1aQ6rfCwSdUoAInGgtKqQlU4AGzfCwIFevZ6ievTg6uP/wKoDFZvdk5Ab1yc43a+JxIGwVoV0ztt/t0WL0M6/SROWjH+G9q16s+pAdazkbzhEQm5cnwSUAERi3Lzl2aRY8d1xha+616+Hzp29Ug7btwfbr72WV559k15bG5AdSCoOypUEmtTNYNTlrTQLKA5pCEgkhhWM/Rc37l6hq+68PJg40Vu5u3t3sP3EE+HJJ6FDB0aNXnTIMJPD6+CBwsRQVJO6Gbw//Lxy//8jsUV3ACIxrLixf/Bm4pT7qnvVKmjf3pvlU9D5p6TA3/8OK1ZAoHBZacNMkS5LLP7QHYBIDCupU853ruzOf98+ePBBuPder5BbgVatvGcAZ54Z8vbGdTOKvcpvXDej8Fh+LP6K6PTXJKcEIBLDSuuUS/Xxx9C7N6xcGWxLS4PbboPhwyE9/ZCPDO3SLGSqKYRe5Xdr0yTqHW+kN0VPdhoCEolhFR562b0bbrkFzj47tPM/6yxYvhzuuKPYzh+8DnXU5a1oUjcDIzYe7kZ6U/RkpzsAkRhWoaGXt9/2ZvesWxdsq1nTq+E/eDCkph76mWKOF0tX1pHeFD3ZKQGIxLgyO+Xt22HYMJgyJbT9/PNh6lQ44YTIBhhBlR4Ck3LREJBIlM1bnk370Ys4fvgC2o9exLzl2WV/qCQLFkDLlqGd/+GHe1M733wzrjt/iPym6MlOdwAiURS2h5pbtsCQIfDss6Htl14Kjz0GjRuHK2Rf+Tn7KBmUmQDM7GxgItAcSAdSgZ+dc3UiHJtIwqlyTR/nYPZsb0x/69Zge/36MGkSXHHFIcXb4l2sPZdIJOUZApoEXAWsBTKAPngJQUQqqEoPNbOzvSv8q64K7fz/+ldYvRr+9KeE6/wlssr1DMA5tw5Idc7lOeeeBjpGNiyRxFTSw8tSH2o6B9OmecXbXn452N60Kbz6Kvzzn3DUUWGOVJJBeZ4B7DazdOBTMxsDbEJFv0UqpazFVhBc+Zqdk8sJOZu477WJ/OrbFaFfNGAAjBoFdTQSK5VXngTwV7w7hUHAjUBT4A+RDEokUZX1ULPgIfHevfvos3Q+N787k4wDews/v+vY46k14xk491w/wpcEU2YCcM59E/hxj5lNAJoGhoREpBJKe6g5duEamm78ijGvPcrpm9YWth+wFKb98nL+/fteLFbnL2FSnllAWcAlgfd+Cmwxs7edczdFODaR5LJvH1e88iQDlswhPf9AYfMXRx/PLV1v4POGJ2E/5/sYoCSa8gwBHe6c22FmfYCnnXN3mtmKMj8lIuX30UfQqxdDVq0qbNqbWo0Jv7qKKWf9gQOp3n+q8bACVtU740d5EkA1M2sE/AkYGeF4RJLL7t1w++0wfjzkB6/ulzU+hVu63sBX9ZoWtsXDCtic3P2M+I+qd8aL8iSAe4CFwPvOuY/N7AS8NQEiUhWLFkHfvt42jQUOO4wVA4cx5PCz2bBjH6lm5DlHkzi5kt68fQ+5+0Nnl1d583qJmPI8BJ4DzCny+3o0C0giJCmGD3JyYOhQr15PUZ07w9SpnJaZybv+RFZl+/LyKW55kap3xqYyF4KZ2Qlm9rKZbTGzH8xsvpkdH43gJLkUTIHMzsnFERw+qFKxtFjz0kte8bainX/duvD007BwIWRm+hZaOKSnFt+lxMOzi2RUnpXAzwL/BhoBjfHuBmZFMihJTomw+UeJlT5/+AG6d/dKOWzcGPzAZZfBF1/A3/6WEGUcGhxeQ9U740h5EoA552Y45w4E/v0LcJEOTJJPvG/+UewdzAsrWHr/RK+Mw+zZwTc3aADPPw9z50KjRr7FHG51M9JiblcxKVl5HgIvNrPheFf9DrgSWGBmRwI457ZFMD5JIvG++cfBdzCNdmzhvjceo91XH4e+sWdPGDcOjjwyyhFGh6p3xo/yJIArA//32oPae+ElhPjecUJiRnnq5MSygjsVc/n0+PR1hmc9Te19RRLascd6O3R16eJThCKhyjMLSA98JSriffOPxnUzSFu/jgdfn8hZGz4vbM83I2XgQHjgAahd28cIRUKVmADM7PLSPuicm1vVg5vZBcCjeJvMPOmcG13V75T4FqvDB2VNT7W8PKZsWsRJTz9EjQP7CtvXH3UM342ZwLm9LvMjbJFSlXYHcHEprzmgSgnAzFKByUBn4DvgYzN7yTn3RVW+VyTcytzG8bPPOGPAAGp/+WXhZw5YCs/+9krqjrqXS84+0Ze4RcpSYgJwzl0T4WP/ElgXWFiGmc0CLgWUACSmlDQ9dfyClXSb+wSMHk3tA8HibbRpQ7Xp07m6TZsoRypSMX5uCt8E2FDk9++As3yKRaRExU1DPSN7NQ++NgF+LPInXL063HUX3HwzpKVFL0CRSvIzARS36uWQ9QVm1g/oB9CgQQOysrIiHFb57Nq1K2Zi8VMynIfhp+cHShxA2p5czp4zk9ZvLMBc8M/1xxYtWDdsGLnHHgvvv+9XqL5Lhr+HssTTOfAzAXyHt7tYgWOAjQe/yTk3FZgK0K5dO9ehQ4eoBFeWrKwsYiUWPyXDecgJPANo++VSRi2cRNPtm4Mv1qoFo0ezsnlzOpx3nn9Bxohk+HsoSzydAz9nAX0MnByoK5QNdAd6VPE7RcKuW2ZN2qx8huPmzw59oUsXmDIFjjsOqnjFlxRF8CTmlGcW0NHAr4BFgd87AllUcRaQc+6AmQ3CKzWdCjzlnFtVxsdEouvFF2HAAI77/vtg2xFHePX7//rXsNTvKXOWkUiElDkLyMxeAVo45zYFfm+EN32zypxzrwKvhuO7RMJq82YYPBjmzAltv+IKmDjRq+UTJqUVwVMCkEgqzzOAzILOP2Az8IsIxSPiL+dgxgwYMgR++inY3rAhPPaYV70zzOK9CJ7Er/IkgCwzWwg8hzdLpzuwOKJRifjhm2+gf394/fXQ9l694KGHvKGfCKhKETw9O5CqKE8toEFmdhlwbqBpqnPuxciGJRI5h3SanU+m2wcvwfDhsGtX8I2ZmTBtGnTqFNF4KlsET88OpKrKsx8AwCfAAufcjcBCM1NFK4lLB9fsr/7VWo65rCsMGhTs/M3ghhtg5cqId/7gddZFa+gfUTON6tVSuHH2p6GbyhwkETbQEX+VeQdgZn3xFmIdCZyIt4L3CeD8yIYmEn4FnWa1vAP0/fhFhrz3LNXz9gff0Lw5TJ8O55wT1uOWNVRTUASvIlf1enYgVVWeZwAD8er2fAjgnFtrZkdHNCpJCLE4Pr0xJ5eWm7/iwdcmcOrmrwrb96ekkjbyVhg50ivpEEYV6dQrMiMo3jfQEf+VZwhor3OusL6tmVVDW0JKGWJyg/c9e7jrw2eZ/48bQzr/FQ1Pos+gx+Gee8Le+UPFhmoqclU/tEsz7b8rVVKeO4C3zexWIMPMOgMDgJcjG5ZURixdcZf3SjZqMb//PvTuTc81wU53T7V0xv36zzz7qz9w3x9PD/8xAyrSqVfkqj7eN9AR/5UnAQwD+gAr8baFfBV4MpJBScXF2oyQ8nR6UYl550649VaYPNmb4x+w/PjTuLnTQPaecBL3RbjTLKtTL5oED89IIy3V2J8XjLW0q/pY3UBH4kOpCcDMUoAVzrlTgWnRCUkqI9ZWk5bnSjbiMS9cCP36wbffBttq14YxY2jTrx+LUso7Ca5qSpvmeXASzMndT1qKcUTNNHJ279dVvURUqQnAOZdvZp+Z2bHOuW9Le6/4K9ZmhJRnbnvEYt62DW66Cf7xj9D2rl294m1Nmxb/uQgpbaim/ehFhyTB/fmOmunVWH7H76IapySf8gwBNQJWmdlHwM8Fjc65SyIWlVRYrM0IKc/4dERifuEFGDjQq+VT4Kij4NFHoUePsBRvq4yShmpKSnbFnReRcCtPArg74lFIlVV2NWkklTU+HdaYN23yFnPNPahIbffuXud/dGzOXC4pCRreswEN/UgklTgIamY1zGwIcAVwCvC+c+7tgn9Ri1DK5eDVpE3qZjDq8lYR7UDmLc9mzfc7OX74glJXrJYkLDE7B888Ay1ahHb+jRvD/Pnw3HMx2/mDlwRL2hpPK3ol0kq7A/gHsB94F+gKtABuiEZQUjnRnBFS8PBywCn5OFIqPYOnSjF//bX3kPfNN0Pb+/aFMWOgbt3KfW8UdWvThCGzPy32Na3olUgrbRpEC+fcX5xzU4A/Ar+JUkwSB3ytQ5OXBxMmwKmnhnb+J5wA//kPTJ0aF51/gSYlPPPQil6JtNISQGGBFOfcgSjEInEk2rOO5i3Ppv3oRXTq8wQrTmztFWv7OTAnISXFm/WzYgXE4b68WtErfiltCKi1me0I/Gx4K4F3BH52zrk6EY9OYlY0Zx3NW57N7XOWc/V7/+b6/z5H9bwi1yMtW3rF2846K+zHjRat6BW/lLYlZGpJr4kUzOCBYGccqavW+U+9zOx/j6HFD/9X2LYvpRozOvag96vTID097MeMNq3oFT+UZxqoyCEKOqvNaz7BIDJXrbm5cPfdTJs8lmouv7D500YnM6zrDXxZP5PeCdD5i/hFCUAqrVubJmRtX8v/je4Q/i9/5x3o0wfWri38I82tVp2Hf/Nnnmp3KfkpqSU+PBWR8lECkLCqcnXPHTtgxAhvA/YiPjzuNG7pMohvjmgM6CGpSDgoAURQLJVnjoYqV/d87TW49lrYsCHYVqcOy68fyXWpp7Ftj/e9dTPSuOuSlgl9LkWiITrlEJNQTG6IEmGVXhvw449w9dVw4YWhnf/FF/P6v/9DD4KdP8DeA/nFfImIVJQSQIQk44bd5V0bUDCn//hhr3B7jzvYe3IzmDEj+IZ69bwSDvPnc+/yHUl3HkWiRUNAERJr5ZmjoTxrAwrujGpv+4Epbz7O79Z+EPrmHj284m316gHJeR5FokV3ABFS0oKoRF7eX54VrWNf/x8XL32Nt6YPCOn8f6hTD15+GWbOLOz8ITnPo0i0KAFESDIu7y+zuuf69YyZchNjXp9Anb2FW0sw8/QL6NRrMlx00SHfmYznUSRaNAQUIcm6vL/YFa15eTBxIowcSfvduwub/++IRoy4YDAfHHtaiXP6k/U8ikSDEkAEJfPy/oIpsIet/R+PvDmJlhtWF76WZylMO7Mb43/dgz1pNcq8ok/m8ygSSUoAEnbzlmdzx5xP6PnubAb/dzbp+UWKt7VqxbvDRjNjQw325uTSRFf0Ir5RApCwe/nJ+cyeM5bmW74ubNuXUo1/nPcX+i6YQof0dN73LzwRCVACkPDZvRvuvJOpj48jtUjxtuWNmnFL1+tZV/84+qp4m0jMUAKQ8MjK8rZiXLeOgjk7u9Oq89BvruaZthepeJtIDPJlGqiZXWFmq8ws38za+RGDhEfqrl3Qvz907Ajr1hW2L8k8nS69JvPUmV7lTk3dFIk9ft0BfA5cDkzx6fgSDgsW8MtrroGtW4Nthx8O48ax+fTfkf/Gl1gMT91MtmJ9IgfzJQE451YDmJkfh5eq2rIFhgyBZ5+letH2Sy/1yjg3bkw3oNsZx/gUYNmqXLlUJAGYc86/g5tlAX93zi0t5T39gH4ADRo0aDtr1qwoRVe6Xbt2UatWLb/DiJqc3P1szsnluPfe5rczniRj547C1/YdcQRrr7+eLb/9LcRJUl/z/U725R1aVTQ9NYVmDWtX+PuS7e+hJDoPsXkOOnbsuMw5d8hwe8TuAMzsLaBhMS+NdM7NL+/3OOemAlMB2rVr5zp06BCeAKsoKyuLWIkl0uYtz2b8i4u4bcFEOn31cchr33fuTMPnnqPlUUf5FF3lXDN8Aa6YR2AGldrhLJn+Hkqj8xBf5yBiCcA51ylS3y1RlJ/Pmnsf5qUFU6izL1jGIbt2fcb94UYu7nkODeOs84fyVS4VSXQqBiclW7cOzj+fYS8+EtL5/7PN7+nSezJzG5zmY3BVoyJzIj49BDazy4CJQH1ggZl96pzr4kcsUoy8PBg/Hm6/HXKDV8nrj2jMsK7X83HTUwHiel6/isyJ+DcL6EXgRT+OLWX4/HPo1Qs+Do7156em8uRZf+Dhs69kb5o376fwann7Wr8irTIVmZNkpyEg8ezdC3fdBWecEdL507o1KR9+yNGTxlGvft3i6/yHSeFWkcMX0H70ooTeP1kkFqgURJIqugiq046vefiNidT5qsg+u+npcOedMHQopKV58/ojeLWsefki0acEkIQKOlt+/pmR786g19KXSKHIepBzzoHp06F586jFNHbhmhI3f69sAtBKX5HSKQEkobEL13D6uk8Y/fpEjsv5vrB9d3oNaj40BgYMgNTUUr4h/MK9+bvuKETKpgSQbHJyGPzcg3Rf8UZI8zuZbbj1gkG8N7iXL2GFe15+JO4o/KY7Ggk3PQROJi+9BC1bhnT+26sfxt8vHMLVf7oHd1ymb6GFe15+uO8o/FZwR5Odk4sjeEejB+VSFUoAyeCHH6B7d69Y28aNhc2v/eJXdOrzBM+36kRGejVfF0F1a9OEUZe3okndjLDMNCrpziFeV/qWdkcjUlkaAkpkzsHMmXDDDbBtW7C9QQM++vs93HfgJLbG0L684ZyXP7RLs5BnABDfK30T7Y5GYoMSQKLasMHbqOXVV0Pbe/aEceP45ZFHJvS+vIm20le1iyQSlAASTX4+TJkCw4bBzp3B9uOO89q7JE/FjURa6ZtodzQSG5QAEsnatdCnD7zzTrDNDAYNggcegBirUS7ll2h3NBIblAASwYEDMG6ct3J3z55ge7Nm3oKu9u39i60ITWOsmkS6o5HYoAQQ7z77DHr3hmXLgm2pqd4Q0O23Q40axX4s2p2xFmaJxB5NA41Xe/d6HXy7dqGdf5s2sHQp3H9/qZ1/tOeUaxqjSOxRAohHS5Z4Hf1993nDPwDVq8OoUfDhh3D66aV+3I/OWNMYRWKPEkA82bULhgzxxvRXrw62//rX3lDQ8OGQllbm1/jRGSfawiyRRKAEEC/efBNatYJHH/UWeIE3q2fyZMBaENsAAApGSURBVHj7be+Bbzn50RlrC0aR2KMEEOt++sl7yPu738HXXwfbL7jA271rwABIqdj/jH50xuEu9SAiVadZQLHsxRe9Dv77YMlmjjzS26/3L3/x5vhXgl9zyjWNUSS2KAHEou+/h8GD4fnnQ9uvuAImToQGDap8CHXGIqIEEEucgxkzvAe9P/0UbG/YEB57DC67zL/YRCTh6BlArPjmG+ja1SvWVrTz790bvvhCnb+IhJ3uAPyWnw+PP+5N4dy1K9iemQnTpkGnTr6FJiKJTQnAT2vWeMXb3nsv2Gbm1e+/7z447DD/YhORhJfwCSAmC5Dt3w8PPwx33eWVdCjQvLlXvO2cc3wLTUSSR0IngJgsQLZ8uTeuv3x5sK1aNRgxAkaO9Eo6iIhEQUI/BI6pAmR79sCtt8KZZ4Z2/m3besXb7rlHnb+IRFVC3wHETAGy99/3rvrXFEk8NWp4nf6NN3p3ACIiUZbQdwC+FyDbudNb0PWb34R2/ueeCytWwNCh6vxFxDcJnQA6nlK/Qu1htXAhnHoqTJoULN5Wu7Y35XPxYjj55MjHICJSioROAK98tqlC7WGxbRv87W9esbZvvw22X3ghrFoF/ftXuHibiEgkJPT4Q07u/gq1V9nzz8PAgfDDD8G2o47ySjj36FHp4m0iIpGQ0AkgajZtgkGDYO7c0Pbu3b3O/+ij/YlLRKQUCT0WcUTN4nfHKqm9wpyDp5+GFi1CO//GjWH+fHjuOXX+IhKzfEkAZjbWzP5nZivM7EUzqxuJ49x5cUvSUkOHXdJSjTsvblnl767x/ffQpQv06gU5OcEX+vb1irddckmVjyEiEkl+3QG8CZzqnDsN+BIYEYmDdGvThLF/bB2yC9XYP7au2irgvDyYMIEzr7nG26axwAknwH/+A1OnwuGHVzl2EZFI8+UZgHPujSK/fgD8MVLHCuvGJ6tXewu6liyhcEPFlBSvfv+990LNmuE5johIFJgrmKPuVwBmLwOznXP/KuH1fkA/gAYNGrSdNWtWNMPzYjhwgKazZpH5z3+Ssj84g+jnzEz+N3QoO1u0iHpMsWLXrl3UqlXL7zB8p/Pg0XmIzXPQsWPHZc65dge3RywBmNlbQMNiXhrpnJsfeM9IoB1wuStHIO3atXNLly4Nb6BlWbbMG+dfsSLYlpbG1z16kDl1KqSnRzeeGJOVlUWHDh38DsN3Og8enYfYPAdmVmwCiNgQkHOu1J1MzKwncBFwfnk6/6jLzYW774aHHvLG/QuceSZMn87XP/5IZpJ3/iIS3/yaBXQBMAy4xDm3248YSvXOO9C6NTz4YLDzz8jwavgvWQKtWvkbn4hIGPi1EGwSUB1407zVsR845/r7FEvQjh1eXf7HHgtt79jR257xxBP9iUtEJAL8mgV0kh/HLdWrr3p1ejZsCLbVqeMNAfXpozIOIpJwVApi61avJv+/DpqEdPHFXuXOJj5vHykiEiHJmwCcgzlzvBo+W7YE2+vVg4kT4corddUvIgktORPAxo0wYIBXr6eoP/8Zxo/3koCISIJLrgTgHDz1FNx8M2zfHmw/5hh44gn4/e/9i01EJMqSJwGsX+8Valu0KLS9f39vumedOv7EJSLik4QuBw148/gfecSbu1+08z/pJMjK8h70qvMXkSSU+HcAV18Nzz4b/D0lxRsCuusuFW8TkaSW+HcAffsGf27VCj78EMaMUecvIkkv8e8AOnSAwYOhfn0YNizpi7eJiBRI/AQAMGGC3xGIiMScxB8CEhGRYikBiIgkKSUAEZEkpQQgIpKklABERJKUEoCISJJSAhARSVIWi/uxl8TMtgDf+B1HQD1gq99BxACdB4/Og0fnITbPwXHOufoHN8ZVAoglZrbUOdfO7zj8pvPg0Xnw6DzE1znQEJCISJJSAhARSVJKAJU31e8AYoTOg0fnwaPzEEfnQM8ARESSlO4ARESSlBKAiEiSUgKoAjMba2b/M7MVZvaimdX1OyY/mNkVZrbKzPLNLC6mv4WLmV1gZmvMbJ2ZDfc7Hj+Y2VNm9oOZfe53LH4ys6ZmttjMVgf+e7jB75jKogRQNW8CpzrnTgO+BEb4HI9fPgcuB97xO5BoMrNUYDLQFWgBXGVmLfyNyhfPABf4HUQMOADc7JxrDpwNDIz1vwclgCpwzr3hnDsQ+PUD4Bg/4/GLc261c26N33H44JfAOufceufcPmAWcKnPMUWdc+4dYJvfcfjNObfJOfdJ4OedwGqgib9RlU4JIHx6Aa/5HYREVRNgQ5HfvyPG/4OX6DCzTKAN8KG/kZQuOfYErgIzewtoWMxLI51z8wPvGYl3+zczmrFFU3nOQxKyYto0rzrJmVkt4AVgiHNuh9/xlEYJoAzOuU6lvW5mPYGLgPNdAi+qKOs8JKnvgKZFfj8G2OhTLBIDzCwNr/Of6Zyb63c8ZdEQUBWY2QXAMOAS59xuv+ORqPsYONnMjjezdKA78JLPMYlPzMyA6cBq59w4v+MpDyWAqpkE1AbeNLNPzewJvwPyg5ldZmbfAecAC8xsod8xRUNgAsAgYCHeA79/O+dW+RtV9JnZc8ASoJmZfWdmvf2OySftgb8C5wX6g0/N7EK/gyqNSkGIiCQp3QGIiCQpJQARkSSlBCAikqSUAEREkpQSgIhIklICkIQXmKbqzOyUcrx3iJnVrMKx/mZmk8rbLuInJQBJBlcB7+Et1CrLEKDSCUAknigBSEIL1GVpD/SmSAIws1Qze8jMVgb2cxhsZtcDjYHFZrY48L5dRT7zRzN7JvDzxWb2oZktN7O3zKxBBWJ6xsyeMLN3zexLM7so0N7SzD4KLCBaYWYnB9rnmdmyQI35flU/KyIe1QKSRNcNeN0596WZbTOzMwIle/sBxwNtnHMHzOxI59w2M7sJ6Oic21rG974HnO2cc2bWB7gFuLkCcWUCvwVOxEs4JwH9gUedczMDpSVSA+/tFYgtA/jYzF5wzv1YgWOJFEsJQBLdVcD4wM+zAr9/AnQCnijYz8E5V9F69scAs82sEZAO/F8FP/9v51w+sNbM1gOn4JVTGGlmxwBznXNrA++93swuC/zcFDgZUAKQKlMCkIRlZkcB5wGnmpnDu6J2ZnYLXinn8tRBKfqeGkV+ngiMc869ZGYdgLsqGN7Bx3bOuWfN7EPg98DCwJ1FPl6yOsc5t9vMsg6KQ6TS9AxAEtkfgX86545zzmU655riXan/GngD6G9m1QDM7MjAZ3biFfgrsNnMmptZCnBZkfbDgezAzz0rEdsVZpZiZicCJwBrzOwEYL1zbgJeVdHTAsf5KdD5n4K31aBIWCgBSCK7CnjxoLYXgB7Ak8C3wAoz+yzQBjAVeK3gITAwHHgFWARsKvI9dwFzzOxdoKznBcVZA7yNt4tcf+fcHuBK4HMz+xRvSOifwOtANTNbAdyLt/WoSFioGqhIlAVmEr3inHve71gkuekOQEQkSekOQEQkSekOQEQkSSkBiIgkKSUAEZEkpQQgIpKklABERJLU/wP36ekN5rybbgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(y_ts,y_ts_pred)\n", "plt.xlabel('Actual lpsa')\n", "plt.ylabel('Pred lpsa')\n", "ymin = np.min(y_ts)\n", "ymax = np.max(y_ts)\n", "plt.plot([ymin,ymax], [ymin,ymax], 'r-', linewidth=3)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also plot the coefficients in the regression model. Remember that all the parameters are normalized so that the coefficients can be compared. We see that `lcavol` has the highest weight, but there are non-zero weights on all the predictors. This makes it hard to see if other factors are significant or not. " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " lcavol 0.5892\n", " lweight 0.0979\n", " age -0.0897\n", " lbph 0.2182\n", " svi 0.3113\n", " lcp -0.0394\n", " gleason 0.0736\n", " noise 0.1214\n", " redundant -0.0394\n", " dependent 0.0341\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":5: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the \"use_line_collection\" keyword argument to True.\n", " plt.stem(w)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASxElEQVR4nO3df5Bd533X8feHtVTWhoyaWiXRSooFCBVDMAqLE5NSColHcuggNXQGu9BftKOYqfuDYWRsOtN/+CPJKMM0ULdCY9wfQ7Ez46qqpqhZIKUEpvkhJQp2bHdbjdNYKzlYcVBC053RD3/5Y1fO1WZl7+qe3XO1z/s1o/E9z3n2PN85V/r47nPOfU6qCknS2ven+i5AkrQ6DHxJaoSBL0mNMPAlqREGviQ14qa+C3gtt956a9122219lyFJN4zPfOYzX66qjYvtG+nAv+222zhx4kTfZUjSDSPJF6+1zykdSWqEgS9Jjegk8JPsTjKd5FSSh67R57uTfC7JM0n+RxfjSpKWbug5/CRjwCPA3cAMcDzJ0ap6dqDPBuAXgN1V9UKSbx92XEnS8nTxCf9O4FRVPV9VF4AngD0L+nw/cLiqXgCoqpc6GFeStAxd3KUzAZwe2J4B3r6gz18C1iX5XeDPAh+uql9d7GBJ9gH7ALZu3brsYo6cPMOBqWnOnp9l04Zx9u/awd6dE8s+jiStNV0EfhZpW7gE503A3wDeBYwDn0jyyar6g2/6wapDwCGAycnJZS3leeTkGR4+/DSzFy8DcOb8LA8ffhrA0JfUvC6mdGaALQPbm4Gzi/T5aFV9vaq+DHwcuKODsa9yYGr61bC/YvbiZQ5MTXc9lCTdcLoI/OPA9iTbkqwH7gWOLujzm8DfTnJTkpuZm/J5roOxr3L2/Oyy2iWpJUNP6VTVpSQPAFPAGPBYVT2T5P75/Qer6rkkHwWeAl4BHq2qzw879kKbNoxzZpFw37RhvOuhJOmG08nSClV1DDi2oO3ggu0DwIEuxruW/bt2XDWHDzC+boz9u3as5LCSdEMY6bV0luvKhdkHn3yKC5dfYcK7dCTpVWsq8GEu9B//9AsAfOR9d/VcjSSNDtfSkaRGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUiE4CP8nuJNNJTiV5aJH9353kq0k+N//nZ7sYV5K0dEM/0zbJGPAIcDcwAxxPcrSqnl3Q9X9W1fcMO54k6fp08Qn/TuBUVT1fVReAJ4A9HRxXktShLgJ/Ajg9sD0z37bQXUn+d5LfTvJXrnWwJPuSnEhy4ty5cx2UJ0mCbgI/i7TVgu3PAm+pqjuAfwccudbBqupQVU1W1eTGjRs7KE+SBN0E/gywZWB7M3B2sENVfa2q/nj+9TFgXZJbOxhbkrREXQT+cWB7km1J1gP3AkcHOyR5U5LMv75zftyXOxhbkrREQ9+lU1WXkjwATAFjwGNV9UyS++f3HwS+D/hnSS4Bs8C9VbVw2keStIKGDnx4dZrm2IK2gwOvfx74+S7GkiRdH79pK0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRnQS+El2J5lOcirJQ6/R728muZzk+7oYV5K0dEMHfpIx4BHgHuB24L4kt1+j3weBqWHHlCQtXxef8O8ETlXV81V1AXgC2LNIv58Afh14qYMxJUnL1EXgTwCnB7Zn5ttelWQC+F7g4OsdLMm+JCeSnDh37lwH5UmSoJvAzyJttWD754B/WVWXX+9gVXWoqiaranLjxo0dlCdJAripg2PMAFsGtjcDZxf0mQSeSAJwK/CeJJeq6kgH40uSlqCLwD8ObE+yDTgD3At8/2CHqtp25XWSXwZ+y7CXpNU1dOBX1aUkDzB3980Y8FhVPZPk/vn9rztvL0laeV18wqeqjgHHFrQtGvRV9cNdjClJWh6/aStJjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mN6OSJV9KN4MjJMxyYmubs+Vk2bRhn/64d7N050XdZ0qox8NWEIyfP8PDhp5m9eBmAM+dnefjw0wCGvprhlI6acGBq+tWwv2L24mUOTE33VJG0+joJ/CS7k0wnOZXkoUX270nyVJLPJTmR5Du7GFdaqrPnZ5fVLq1FQwd+kjHgEeAe4HbgviS3L+j2MeCOqvrrwD8FHh12XGk5Nm0YX1a7tBZ18Qn/TuBUVT1fVReAJ4A9gx2q6o+rquY3bwEKaRXt37WD8XVjV7WNrxtj/64dPVUkrb4uAn8COD2wPTPfdpUk35vk94H/zNyn/EUl2Tc/7XPi3LlzHZQnzV2Yff9738r6sbm/8hMbxnn/e9/qBVs1pYu7dLJI2zd9gq+q3wB+I8l3Af8aePdiB6uqQ8AhgMnJSX8TUGf27pzg8U+/AMBH3ndXz9VIq6+LT/gzwJaB7c3A2Wt1rqqPA38hya0djC1JWqIuAv84sD3JtiTrgXuBo4MdkvzFJJl//TZgPfByB2NLkpZo6CmdqrqU5AFgChgDHquqZ5LcP7//IPAPgR9MchGYBf7RwEVcSdIq6OSbtlV1DDi2oO3gwOsPAh/sYixJ0vXxm7aS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1opP18KXXcuTkGQ5MTXP2/CybNoyzf9cOHx4u9cDA14o6cvIMDx9+mtmLlwE4c36Whw8/DWDoS6vMKR2tqANT06+G/RWzFy9zYGq6p4qkdhn4WlFnz88uq13SyjHwtaI2bRhfVrukldNJ4CfZnWQ6yakkDy2y/x8neWr+z+8luaOLcTX69u/awfi6savaxteNsX/Xjp4qkto19EXbJGPAI8DdwAxwPMnRqnp2oNsXgL9TVf83yT3AIeDtw46t0XflwuyDTz7FhcuvMOFdOlJvurhL507gVFU9D5DkCWAP8GrgV9XvDfT/JLC5g3F1g9i7c4LHP/0CAB953109VyO1q4spnQng9MD2zHzbtfwo8NsdjCtJWoYuPuFnkbZatGPyd5kL/O+85sGSfcA+gK1bt3ZQniQJuvmEPwNsGdjeDJxd2CnJXwMeBfZU1cvXOlhVHaqqyaqa3LhxYwflSZKgm8A/DmxPsi3JeuBe4OhghyRbgcPAD1TVH3QwpiRpmYae0qmqS0keAKaAMeCxqnomyf3z+w8CPwt8G/ALSQAuVdXksGNLkpauk7V0quoYcGxB28GB1z8G/FgXY0mSro/ftJWkRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjehkeWRJul5HTp7hwNQ0Z8/PsmnDOPt37WDvztd6LLaul4EvqTdHTp7h4cNPM3vxMgBnzs/y8OGnAQz9FeCUjqTeHJiafjXsr5i9eJkDU9M9VbS2GfiSenP2/Oyy2jUcA19SbzZtGF9Wu4Zj4Evqzf5dOxhfN3ZV2/i6Mfbv2tFTRWubF20l9ebKhdkHn3yKC5dfYcK7dFZUJ4GfZDfwYWAMeLSqPrBg/3cAvwS8DfiZqvpQF+OOMm81k5Zm784JHv/0CwB85H139VzN2jZ04CcZAx4B7gZmgONJjlbVswPdvgL8JLB32PFuBN5qJmkUdTGHfydwqqqer6oLwBPAnsEOVfVSVR0HLnYw3sjzVjNJo6iLwJ8ATg9sz8y3XZck+5KcSHLi3LlzQxfXB281kzSKugj8LNJW13uwqjpUVZNVNblx48YhyuqPt5pJGkVdBP4MsGVgezNwtoPj3rC81UzSKOoi8I8D25NsS7IeuBc42sFxb1h7d07w/ve+lfVjc6d3YsM473/vW71gK6lXQ9+lU1WXkjwATDF3W+ZjVfVMkvvn9x9M8ibgBPAG4JUkPw3cXlVfG3b8UeWtZhp13jrcnk7uw6+qY8CxBW0HB15/ibmpHkkjwFuH2+TSClKDvHW4TQa+1CBvHW6TgS81yFuH22TgSw3y1uE2uVqm1CBXqWyTgS81yluH2+OUjiQ1wsCXpEYY+JLUCANfkhrhRVtJGhErvb6RgS9JI2A11jdySkeSRsBqrG9k4EvSCFiN9Y0MfEkaAauxvpGBL0kjYDXWN/KirSSNgNVY38jAl6QRsdLrGzmlI0mNMPAlqRGdBH6S3Ummk5xK8tAi+5Pk387vfyrJ27oYV5K0dEMHfpIx4BHgHuB24L4kty/odg+wff7PPuAXhx1XkrQ8XVy0vRM4VVXPAyR5AtgDPDvQZw/wq1VVwCeTbEjy5qp6sYPxv8nu3/1PvOncab74v96wEodfsh9+8WsAvdbxRy9/HYDbvu2W3mqA0TgXo1SH78vo1TFK78mXNm6BFbho20XgTwCnB7ZngLcvoc8E8E2Bn2Qfc78FsHXr1usq6I23fAs3f3Xs9TuusJvX91/Dn1y4/PqdVsEonAsYnTp8X642CnWM0nvyxlu+ZUWO3UXgZ5G2uo4+c41Vh4BDAJOTk4v2eT17Hv3Q9fxY597SdwHAg//+E0D/j7AbhXMBo1OH78vVRqGOFt6TLi7azgBbBrY3A2evo48kaQV1EfjHge1JtiVZD9wLHF3Q5yjwg/N367wD+OpKzd9LkhY39JROVV1K8gAwBYwBj1XVM0nun99/EDgGvAc4BfwJ8CPDjitJWp5OllaoqmPMhfpg28GB1wX8eBdjSZKuj9+0laRGGPiS1AgDX1plR06e4eQL5/nUF77COz/wOxw5eabvkprXynti4Eur6MqDqi9cfgX4xoOq12rA3Ahaek8MfGkVrcaDqrU8Lb0nBr60ilbjQdVanpbeEwNfWkWr8aBqLU9L74mBL62i1XhQtZanpffEZ9pKq+jKA6kPTE1z9vwsm1bgQdVanpbeEwNfWmV7d06syTC5kbXynjilI0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwz8NayVJV8lLY2Bv0a1tOSrpKUx8NeolpZ8lbQ0QwV+kjcm+a9J/nD+v996jX6PJXkpyeeHGU9L19KSr5KWZthP+A8BH6uq7cDH5rcX88vA7iHH0jK0tOSrpKUZNvD3AL8y//pXgL2LdaqqjwNfGXIsLUNLS75KWpphV8v8c1X1IkBVvZjk24ctKMk+YB/A1q1bhz1cs1pa8lXS0rxu4Cf5b8CbFtn1M92XA1V1CDgEMDk5WSsxRitaWfJV0tK8buBX1buvtS/J/0ny5vlP928GXuq0OklSZ4adwz8K/ND86x8CfnPI40mSVsiwgf8B4O4kfwjcPb9Nkk1Jjl3plORx4BPAjiQzSX50yHElScs01EXbqnoZeNci7WeB9wxs3zfMOJKk4flNW0lqRKpG90aYJOeAL17nj98KfLnDcm5knoureT6u5vn4hrVwLt5SVRsX2zHSgT+MJCeqarLvOkaB5+Jqno+reT6+Ya2fC6d0JKkRBr4kNWItB/6hvgsYIZ6Lq3k+rub5+IY1fS7W7By+JOlqa/kTviRpgIEvSY1Yc4GfZHeS6SSnklzrgSxNSLIlyX9P8lySZ5L8VN819S3JWJKTSX6r71r6lmRDkieT/P7835G7+q6pT0n++fy/k88neTzJn+67pq6tqcBPMgY8AtwD3A7cl+T2fqvq1SXgX1TVXwbeAfx44+cD4KeA5/ouYkR8GPhoVX0HcAcNn5ckE8BPApNV9VeBMeDefqvq3poKfOBO4FRVPV9VF4AnmHsqV5Oq6sWq+uz86//H3D/oZhfIT7IZ+PvAo33X0rckbwC+C/gPAFV1oarO91tV724CxpPcBNwMnO25ns6ttcCfAE4PbM/QcMANSnIbsBP4VL+V9OrngAeBV/ouZAT8eeAc8EvzU1yPJrml76L6UlVngA8BLwAvAl+tqv/Sb1XdW2uBn0Xamr/vNMmfAX4d+Omq+lrf9fQhyfcAL1XVZ/quZUTcBLwN+MWq2gl8HWj2mleSb2VuNmAbsAm4Jck/6beq7q21wJ8Btgxsb2YN/lq2HEnWMRf2v1ZVh/uup0fvBP5Bkj9ibqrv7yX5j/2W1KsZYKaqrvzG9yRz/wNo1buBL1TVuaq6CBwG/lbPNXVurQX+cWB7km1J1jN30eVozzX1JkmYm6N9rqr+Td/19KmqHq6qzVV1G3N/L36nqtbcJ7ilqqovAaeT7JhvehfwbI8l9e0F4B1Jbp7/d/Mu1uBF7KEegDJqqupSkgeAKeausj9WVc/0XFaf3gn8APB0ks/Nt/2rqjr2Gj+jdvwE8GvzH46eB36k53p6U1WfSvIk8Fnm7m47yRpcZsGlFSSpEWttSkeSdA0GviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrE/wecCQNRr7oR6gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "w = regr.coef_\n", "for name, wi in zip(names_x, w):\n", " print('%10s %9.4f' % (name, wi))\n", " \n", "plt.stem(w) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LASSO: Linear Model with L1 Regularization\n", "\n", "The failure of the linear model motivates us to use regularization to try to select only the \"useful\" features. We will demonstrate how to use the Lasso technique. The `sklearn` package has several excellent routines for this. We first import the `model_selection` sub-package for the k-fold cross validation." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "import sklearn.model_selection " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When considering the LASSO, we consider a set of models with different levels of regularization `alpha`. Higher values of `alpha` imply greater regularization. Similar to the [polynomial example](./polyfit.ipynb), we use k-fold cross validation to determine the appropriate `alpha`. That is, for each `alpha` value, we evaluate the test error on different training / tests spilt. For larger data sets, this exhaustive search is time-consuming. But, it should finish very fast for this small set. " ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# Create a k-fold cross validation object\n", "nfold = 10\n", "kf = sklearn.model_selection.KFold(n_splits=nfold,shuffle=True)\n", "\n", "# Regularization values to test\n", "nalpha = 100\n", "alphas = np.logspace(-3,1,nalpha)\n", "\n", "# MSE for each alpha and fold value\n", "mse = np.zeros((nalpha,nfold))\n", "for ifold, ind in enumerate(kf.split(X)):\n", " \n", " # Get the training data in the split\n", " Itr,Its = ind\n", " X_tr = X[Itr,:]\n", " y_tr = y[Itr]\n", " X_ts = X[Its,:]\n", " y_ts = y[Its]\n", " \n", " # Compute the lasso path for the split\n", " for ia, a in enumerate(alphas):\n", " \n", " # Create a LASSO model object\n", " model = linear_model.Lasso(alpha=a)\n", " \n", " # Fit the model on the training data\n", " model.fit(X_tr,y_tr)\n", " \n", " # Compute the prediction error on the test data\n", " y_ts_pred = model.predict(X_ts)\n", " mse[ia,ifold] = np.mean((y_ts_pred-y_ts)**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now plot the test MSE as a function of the regularization parameter." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZRcdZ338fe39zXpTjrprJCAISyyZoEAQisqiwoujIAMIKMggziLjiPI86iMj8M4HueMCwIRGWQGQRmDIAYFkQaVEBIihCQQCCEJWcjW6X2tqu/zR1U3TdNdqe6u21Vd9/M6p09+t+5Sn7rdud+62++auyMiIuGVl+kAIiKSWSoEIiIhp0IgIhJyKgQiIiGnQiAiEnIFmQ4wXDU1NT5nzpwRzdvW1kZ5eXl6A6VBtuaC7M2mXMOjXMOTi7mee+65fe4+ZdCR7j6ufhYsWOAj9cQTT4x43iBlay737M2mXMOjXMOTi7mA1T7EdlWHhkREQk6FQEQk5FQIRERCToVARCTkVAhEREJOhUBEJORUCEREQk6FQEQk5AIrBGZ2p5ntMbN1Q4y/1MzWJn6eNrPjg8oiItnlottXcNHtKwJrj3b+m1d2BJpvpFlvXtkxzDWdmiD3CO4Czkky/nXgTHc/DvgmsDTALCIiMoTA+hpy96fMbE6S8U/3G3wGmBVUFhHJvN5vtj//3JIMJ5GBsqXTuc8Ajww10syuBq4GqK2tpb6+fkRv0traOuJ5g5StuSB7synX8GRDrsbG+GGN+vr6vnZra3TQ10fbHur9Um1Ho1EaGxvTmikdWaPRaCC/x4wXAjN7L/FCcPpQ07j7UhKHjhYuXOh1dXUjeq/6+npGOm+QsjUXZG825RqebMh168b4HkFd3ZK+dkVFF1VVxe94fbTtod4v1fbNKx+hqqoqrZnSkbWxsTGQ32NGC4GZHQfcAZzr7vszmUVEJKwydvmomR0CLAMuc/dXMpVDRCTsAtsjMLN7gTqgxsy2A18HCgHc/Tbga8Bk4EdmBhBx94VB5RERkcEFedXQJQcZ/1ngs0G9v4iIpEZ3FouIhJwKgYhIyKkQiEhgBnalINlJhUBEJORUCEREQk6FQEQk5FQIRERCToVARCTkVAhEREJOhUBEJORUCEREQk6FQEQk5FQIRCStdDfx+KNCICIScioEIiIhp0IgIhJyKgQiIiGnQiAiEnIqBCIiIadCICIScioEIiIhp0IgIqOmm8jGNxUCEZGQUyEQEQk5FQIRkZBTIRARCTkVAhGRkFMhEBEJuYJMBxCR8HL3+A/QE40Rc4c0tIFRze/uxNzTmikdWT0xnG4W1IKDsnDhQl+9evWI5q2vr6euri69gdIgW3NB9mZTruEJOlfvPQQ//9ySpG135+ozDucL9/6Fjp5oYHly1aRiWHPTh0Y0r5k95+4LBxunPQIRGZH+G/lUtHdH2Lq/nc/evZqSwjxmTCxhUU2U5/YXgMHFi2Zz36o3gNG3gVHNf9dTr1JSWpLWTGnJ2tOZ0roeLp0jEJHAdfZEeWlXC+3dUW46/xiOnTmR2ZPK+Ni8ImZWlzKzqpTr3jePmVXpaY92/smllvZM6chaXmiB/H5UCEQkcL9/aTeRmHP4lAquOHUOeRbMBk1GRoVARAL3v89tpyg/j4mlOhqdjfRbEZFAdUdirNrSwLQJJZj2BLKS9ghEJFD7WruIOdRUFmc6igwhsEJgZnea2R4zWzfEeDOz75vZJjNba2YnBZVFRDLD3dnb2sWCQ6spLczPdBwZQpB7BHcB5yQZfy4wL/FzNXBrgFlEJAPauqN09sS4cMGsTEeRJAIrBO7+FNCQZJILgLs97hmgysymB5VHRMbe3pYuzOBDx+m/djbL5MnimcAb/Ya3J17bNXBCM7ua+F4DtbW11NfXj+gNW1tbRzxvkLI1F2RvNuUaniByNTZ2APG7lgdrP/aHJ9jX2kVFAax55s+DTtPaGj3ockbSTiVfsnY0GqWxsTGtmdKRNRqNBvL3lclCMNjlA4P2d+HuS4GlEO9iYqS3yof19v/RyNZsyjU8QeS6dWP8zuK6uiWDtm3aYcR8NTMmVVBXd+ag01RUdFFVVZx0OSNpp5IvWfvmlY9QVVWV1kzpyNrY2BjI31cmrxraDszuNzwL2JmhLCKSZg+v3UlBnjGhtDDTUeQgMlkIHgIuT1w9dArQ5O7vOCwkIuNPNOY8tmE31WVFuot4HAjs0JCZ3QvUATVmth34OlAI4O63AcuB84BNQDtwZVBZRGT0htPJXGNHD23dUWZPKgs6lqRBYIXA3S85yHgHPh/U+4tI5jS0dlFTUcSEEnVeMB7ozmIRSatozDnQ0cN5x05XlxLjhAqBiKTVgfZu3OEjx8/IdBRJkQqBiKTV/rZuCvONBYdUZzqKpEiFQETSpqGtm6b2HiaXF5GXp8NC44UKgYikzf2r38CBKeppdFxRIRCRtHB3/mflVipLCigr0tVC44kKgYikRVNHD280dFCrvYFxR4VARNJid3MXNRXFVJcXZTqKDJMKgYiMWldPlMaOHi5ZPFtdSoxDKgQiMmq7W7oAuGTxIRlOIiOhQiAio9IVibK3pYvqskJmVJVmOo6MgAqBiIzKb9e9SSTmTK0syXQUGSEVAhEZlXue2UZxQR4TS3XJ6Hil35yIDOlgXU+3d0d4cUczs6tL1cHcOKY9AhEZsT0tXRTl5+lO4nFOhUBERiQac/a1dnPusdMozNemZDzTb09ERmR/WzfRmHPpyYdmOoqMkgqBiIzInuZOSgvzWTRH3U2PdyoEIjJsr+xuoa07ytTKYp0kzgEqBCIybH98dR8A1eWFGU4i6aBCICLDtuK1fRQX5FFckJ/pKJIGKgQiMizRmLPy9QYmlGpvIFeoEIjIsKzf2URLZ4QJJbofNVfoNykiw/L0a/sBmFAy8j2C3juV6+vr33bXcrrao53/hpNLqatbktZM6chaX19PELRHICLDsuK1/bxragVFBdp85Ioh9wjMbIK7Nw8x7hB33xZcLBHJRjF3Vm1p4BMnzeKV3S3Dmneo/ook85KV9Prehpk9PmDcrwJJIyIZd9HtK/o6mxuorStCe3eUJYdPHuNUEqRkhaD/XSKTkowTkZBo7owAcMphKgS5JFkh8CHagw2LSAg0d/Rw5LRKJukB9Tkl2VVDU83si8S//fe2SQxPCTyZiGSVmDstXREdFspByQrBj4HKQdoAdwSWSESyUmtnBHdYosNCOWfIQuDuN41lEBHJbi1d8fMDi+cOPGUo492Q5wjM7Cozm5dom5ndaWZNZrbWzE4cu4gikg1aOiOUFuZTVabzA7km2aGhvwfuSrQvAY4HDgNOBL4PvCfQZCKSNaIxp6Wzh5qK4T2SUvcOjA/JrhqKuHtPov1h4G533+/uvwfKg48mItnipV3NxBwq1b9QTkpWCGJmNt3MSoCzgN/3G1eaysLN7Bwz22hmm8zs+kHGTzSzX5vZC2a23syuHF58ERkLq7Y0ACoEuSpZIfgasBrYAjzk7usBzOxMYPPBFmxm+cAtwLnA0cAlZnb0gMk+D2xw9+OBOuC7ZqYDkCJZZvWWAxTl6/kDuSrZVUMPm9mhQKW7H+g3ajVwUQrLXgxscvfNAGZ2H3ABsKH/2wCVFn/WXQXQAESG9xFEZLR6u5QY7Ji+u/PslgbtDeSwZJ3Ofbxfe7BJlh1k2TOBN/oNbwdOHjDND4GHgJ3E71O4yN1jg2S5GrgaoLa2dsRdsba2tgbWjetoZGsuyN5syjU8B8vV2NgBxLuFHtjujjp7W5zaMqOxsbFvOcnm6W2PNlemhC1XshL/v8DziR94e/9CzsELwWDVY2DXFGcnlv8+4HDgMTP748BeT919KbAUYOHChV5XV3eQtx5cfX09I503SNmaC7I3m3INz8Fy3boxvkdQV7fkHe29LV1AG1OrKykrKujrpz/ZPL3t0ebKlLDlSlYIPkH8ENBxwIPAve6+aRjL3g7M7jc8i/g3//6uBP7N3R3YZGavA0cCzw7jfUQkQC2dPUwsLaS0MLXzA7pkdPwZ8mSxuz/g7hcDZwKvET+R+6fEyeJUrALmmdncxAngi4kfBupvG/ErkjCzWmA+KZyIFpGx09IVYdGc6qEOEUsOSOURQ51AE9BM/P6BklQW7O4R4Drgd8BLwC/cfb2ZXWNm1yQm+yZwqpm9CDwOfMXd9w3zM4hIQHqiMTp7Yiyco24lclmyk8XvJX5H8WLi9xB8z91XD2fh7r4cWD7gtdv6tXcCHxzOMkVk7PQ+f2DRnEk88fKeDKeRoCQ7R/A4sBb4E1AMXG5ml/eOdPe/CzibiGRYU3sP+WYcN2tipqNIgJIVAt3lKxJi7k5jRzcTSgsozNeD6nNZshvKfjqWQURkbCW7iQxg4+4WeqJOVVnhQZelK4XGN5V5ERnUkxv3AjCxVL2+5DoVAhEZVP3GvZQW5lNcoM1Erjvob9jMTkvlNRHJHdGYs3prQ0qHhWT8S6XU/yDF10QkRzR19MTPD5SqEIRBsvsIlgCnAlPM7Iv9Rk0A1BetSA5r6uihvCifCvU4GgrJ9giKiHcNXUC8Z9Den2bgwuCjiUgmuDuN7T2c9q4a8tStRCgku3z0SeBJM7vL3bcCmFkeUDGwd1ARGR8OdskoQEdPjO5ojLr5U3nw+R1DTqdLRnNHKucIbjazCWZWTvyhMhvN7MsB5xKRDGls7wbgzPlTMpxExkoqheDoxB7AR4n3G3QIcFmgqUQkY/a3dVNenM/MqpQeTS45IJVCUGhmhcQLwYPu3sM7HzAjIjlg895W2rujTC4vznQUGUOpFILbiT/Avhx4KvEcY50jEMlBD6/dBcDkct1NHCYHvTbM3b8PfL/fS1sTXVSLSI55eO1OKosLKNLdxKGSyp3FtWb2EzN7JDF8NHBF4MlEJC0uun0FN6/sOOh07d0RXtndyuQK7Q2ETSpl/y7iTxmbkRh+BfiHoAKJSGbsb+0mz2CSDguFzpCFwMx6DxvVuPsvgBj0PYIyOgbZRGSMuDv727o59fAaPXsghJL9xp9N/NtmZpNJXClkZqcQf4axiOSI9u4oXZEYHzl+eqajSAYkO1nce2/5F4GHgMPN7M/AFNTFhEhO2dfahQFnHzONZWsGv5tYdxLnrmSFoH9ncw8Qv5nMgC7g/cSfZywiWSiVriR6tXdH2NvazaTyIqrKdH4gjJIVgnzinc4N7HWqLLg4IjLWHnx+J9GYUztBN5GFVbJCsMvd/2XMkojImHN3/nvFVsqK8qkoVpfTYZXsZLH6nxXJca1dETbsamZqZTGmLqdDK1khOGvMUohIRuxu7qKyuICaCh0WCrMhC4G7N4xlEBEZnYtuX9F3kjgVPdEYDW3dfGLBLPLztDcQZrpzRCSk9jR34cBfn3JopqNIhqkQiIRQY3s3u5o7qSot5F1TKzIdRzJMhUAkhG55YhPRmDN7kh4+IyoEIqHT1RPlp09vpaaiiLIiXTIqKTyPQESy13DuIO61vbEDDGZVH3xvQN1KhIP2CERCpL07wr7Wbj596hyKC/IzHUeyhAqBSEi4O1v3t5OfZ1xbd3im40gWUSEQCYlfrH6D5s4Is6tL1bmcvI0Kgcg4M9wbxwAiMef//eYlKksKmFqpu4jl7VQIRHKcu7O73emOxDisplx9Csk7BFoIzOwcM9toZpvM7Pohpqkzs+fNbL2ZPRlkHpEwamjvobUHvviBIygp1AlieafACoGZ5QO3AOcCRwOXmNnRA6apAn4EnO/uxwB/FVQekTBqaOtmy742ivPhM6fPzXQcyVJB7hEsBja5+2Z37wbuAy4YMM2ngGXuvg3A3fcEmEdk3BrJeQGAm369nmjMmV5uFOih9DKEIG8omwm80W94O3DygGmOAArNrB6oBL7n7ncPXJCZXQ1cDVBbW0t9ff2IArW2to543iBlay7I3mxhy9XY2AFAfX19yu3Wbmfl687kEqOAWMrz/u18+tpBC9vvcbSCyhVkIRjsjJQP8v4LiD/7oBRYYWbPuPsrb5vJfSmwFGDhwoVeV1c3okD19fWMdN4gZWsuyN5sYct168b43kBd3ZKU2j/Y8Gc272jiqOkTqCjOp7mpibq6upTmHUth+z2OVlC5gtxX3A7M7jc8C9g5yDS/dfc2d98HPAUcH2AmkXFjpIeDALY2tBOJOt+58DjydJWQHESQhWAVMM/M5ppZEXAx8NCAaR4E3mNmBWZWRvzQ0UsBZhLJeb9dt4t9rd3MqCrh3TMnZjqOjAOBHRpy94iZXQf8DsgH7nT39WZ2TWL8be7+kpn9FlgLxIA73H1dUJlEcl13JMYNy16kvCifGVXqYlpSE2jvo+6+HFg+4LXbBgx/B/hOkDlExouR9Cbay93ZvK+NrkiUI6ZW6pCQpEzXk4nkiD0tXTR19PDV846itEg3jknq9DwCkRywYWczWxvamVhayGWnHMpv1u4a1vx67kC4aY9AJINGc2VQr0gsxrX3PEdBnnH4FPUlJMOnPQKRcczd2by3jebOCPNrKyjU3cMyAvqrERlj6dgL6PVmcxcH2nv4yjnzqSwpTMsyJXxUCETGQDo3/r1WbWngjYZ2qssKueo9h6V12RIuKgQiAQli49+rOxLj2nvWUFyQx2E6LyCjpHMEIml088oObt24ItCrcGLuvLqnlWjMmVdbQUHeyL7P6Uoh6aU9ApFRCvKb/2C2NbTT2hXh3y88jrIifZeT0VMhEElR/w3+WG/8e/2ofhO7m7uYNqGYjxw/Y8zfX3KTCoEIQ2/kM7XBH8ybTZ38+283Mrm8iEMmlWU6juQQ7VdKqPTvy2c0/fqMtT0tXWxtaOeDR9dyoL17xCeHx8NnlbGnPQLJCePhG/1IuDtLn3qN1/e1MbG0kB986kR1Jidppz0CGVd6N+p/O390PXWOBzF3tu5v51+Xv8yk8iIOrymnuECdyUn6haYQXHT7ChobO6ire+cGZKjDBWPVHrhRy6Z8g62zTOYLi+bOHl7Z3UJTR4Rr6w5n9ZYG3SsggdGhIZEs0xWJ8cnbVtDcEWFuTTn/fM6RoyoCP//cEm44WQ+pkaGFZo9AZDxo64qwcXcLxQX5zJ9WycRS9R8kwdMegUiW+O26XWzY1QzA/dcsURGQMaM9ApEMi7nzjYfWc9fTWygvymdebSVHTZ+Q6VgSIioEIhnU3h1h8942Vm05wJWnzWHdjqa0XB4aphPrMno6NCSSAZFYjH/59QZe3NFMZyTGbX99El//yDG6R0AyQnsEImOoubOHXU0d7GzsZM22RqZWFjOrupRz3j0909EkxLRHIDIGOnuibN3fzqk3/4FtDR2UFeXz6+tOZ25NuR4vKRmnPQKRgESiMRraurnw1qd5YXsTABecMIPX9rRSXlzAu2dOTOv76byAjJS+ioikUU/M2d3cyWU/WcmabY28vr+dxo4eZleXcuLsKr538YmUF+v7l2QXFQKRUWju7KGxvZttDe2c859PsbnJ2bK/ne0HOpg2sYRjZkzgsX88gxlVpRQV6L+bZCd9NRFJUcydju4o9z67jdf3tdHaFeH4mx7FHQyYW1POlFJj2qRKHrrudC5e+gxAoH0E6XCQpIO+ooj04+5EojFauyI8+PwOth9o59U9rXzgP55k9ZYDrNvZzA3LXmR/WzcFecbfnzWPI6dVsuDQan521SlMKjHKigrUQZyMK9ojkFBq6uihpbOHzp4YNy9/iY27W+jqiXHsNx6ltSsCwN/f9zwAxQV5nHRINW1dEcqK8vnJpxfx5ftfwMz4h/cfwYrX9mfyo4iMmgqB5LRILEZHd5R7Vm7llTdb2LCrmc6eKMff9GjfNP/15y3k5xnFhXl89ISZ/OHl3RQX5HPLpSdx47IXycsz7rhiYV932IdOLs/oN34dDpJ0UyGQnNATjdHRE+W+Z7exdX87nT1RTr35cXY2dQJw4wPrKC/KB4Oq0kI++57D+OWa7ZQW5rPs2tP41I/jx/O/cf4xvJTo+O2I2kry8nSIR3KfCoGMK+5OR0+Up3dG2Lo/Snt3hMXf+j17WroAuH7Zi5hBaWE+Z86dwpqtBygtyufOTy9iZlVp3wncz515OH94eQ8A+drYS8ipEMioxGJOLOY40NoVIRKLgcePwUeiMTBo6ewhGnOM+Dd3dx/00Ip7fDkH2rrpikTpiToPPr+DHY0ddHRHOfd7f+TlXc048OIOMIOywnzeM28Kq7bsp7QwnzuuWMSXfvE8Zsb3Lj6x73DOrOqysVwtaafDQRKk0BSCSDRGZ8RZv7OJtsTJwPU743d79g6v29E04vaL2+NtB9Zub+w74fjCG2+1/7LtAK2d8WnWbDtAS2cEcF49EKWlsweAZ19voLmzBxye2byfmDtNHfFxf3p1H43tPYDzxMY9NLZ348BjG3bT0NYNwG/W7mJ/axcO/PK57ext6cJx7lm5lTebO8Fh6VOvsbOxA3f47qMb2dbQjrvzf3+1js1723Ccz9+zho27W+jujnHhrU+zYVczMXdO//Yf2N3cScxh/v95hK5IrG8dv/vrv+tr9z8Gf+w33mrPu/GRvva7vrqcaMwBOOyG35BocuI3H+ubpveEbWG+8e6ZE2lsL6G0MJ8r58d4aEcJeWZ895PH923wZ08q0xU7IsMUmkLQ1Blha4vzoe//qe+1/m2AD//gTyNuf+SHb7XP/+Gf+9oX3PJW+2M/erqv/fF+7Q27AOLHsj+Z2KABfYcxev31T1b2ta/8r1V97avuXt3X/vzP1vS1v3T/C33tGx9Y19f+1+Uv97V/+MQmIH4d/MNrd9LaFSHPjJffbKY7EiPmUFSQR1F+HnkGi+dM4o+v7iPP4KMnzqSkMJ9la7YDcPmSOfz3M1sA44pT5/DTp7cAcNkph3L3inj7kwtnc9+qbTjwsRNnsmzNDgA+cdIsHvjLdsyMz5w+l7tXbKUwz/jhpSfx1WUvkp9n/PRvFr+1wa/syvmeOrUXIGMlNPcRVBYXMKPcuP2yBcybWsG8qRXcftmCtw0vPUj7iKkV/PjyhRwxsF1bwR2XL+SI2nj7J1e81b7z02+1/+vTi5hfW8H82gruuvKt9pcWFDO/tpIjp1XyP585mSOnxds/++zJ3HvVKRw1rZKjplVy/zVLOHp6JUdPr2TZtadyzPQJHDNjAg9/4XTePWMCx86cwKP/eAbHzZzIcbMm8uSX6zhh1kROmD2RZ796FicdUsWCQ6pYd9PZLDq0msVzqnn95g+xeM4kFs2ZxF++9kFOOqSaE2ZX8fiX6jh25kQOnZDHz646hfnTKplXW8l/XHQCh00pZ05NOTecdxT/+IEjmFFVyoyqUq464zCmTyxl+sQSPnP6XKZPLGH6xBKuOuOwvmm+cNY8ZlWXMbu6jC+ffSSHTCrjkEll/NPZ85lVXcbMqlKuPG0uUyuLqS4v4ojaSh3DFwlYoIXAzM4xs41mtsnMrk8y3SIzi5rZhUFlKSrIo7LIOPuYaUwqL2JSeRFnHzPtbcMfPEi7uryIDxxdS/XAdlkR7z+6luqyePuso95qv+/It9rvPXIqVWVFVJUVUTf/rfaxUwqoKitkYmkhp8+rYWJpvH3qu2pYcvhkJpQWMqG0kEVzJlFZUkhlSSEnHVJNRUkBFYnOy8qLCygrKuCI2kpKi/IpLczn0MnlFBfmU1yQz9QJJRTm51GQn0dFcQF5eaZDKCICBHhoyMzygVuADwDbgVVm9pC7bxhkum8Dv3vnUkTCRYeDJBOC3CNYDGxy983u3g3cB1wwyHRfAH4J7Akwi4iIDCHIk8UzgTf6DW8HTu4/gZnNBD4GvA9YNNSCzOxq4GqA2tpa6uvrhx2msbGDaDRKfX09jY0dAH3L6T+ciXZra3TIaTKdb7B1lg35kq2zTGaNRqM0NjamPP3fzudt+YLS2toa+HuMhHINT1C5giwEgx2A9gHD/wl8xd2jyY5Xu/tSYCnAwoULva6ubthhbt24gsbGRurq6rh1Y/zKk7q6JX3jeocz0a6o6KKqqnjQaTKdb7B1lg35kq2zTGa9eeUjVFVVJZ1mBH++o1ZfX89I/t8ETbmGJ6hcQRaC7cDsfsOzgJ0DplkI3JcoAjXAeWYWcfdfBZhLRET6CbIQrALmmdlcYAdwMfCp/hO4+9zetpndBTysIiC5SCeBJZsFVgjcPWJm1xG/GigfuNPd15vZNYnxtwX13iLZQBt/GS8CvbPY3ZcDywe8NmgBcPdPB5lFZCzccHJp3/kIkfEiNF1MiKRT/2/7+uYv450KgUiKtMGXXKVCIDKAvu1L2KgQSGhpgy8Sp0IgoaINvsg7qRBIztE3fZHhUSGQcUsbfJH0UCGQcaV3g5+NHYKJjFcqBJJWQ31LT1dbRNJPhSDHjXYD3PvNO9mGWRtqkXHO3cfVz4IFC3yknnjiiRHPG6RszeWevdmUa3iUa3hyMRew2ofYrobm4fUiIjI4FQIRkZBTIRARCTkVAhGRkFMhEBEJORUCEZGQUyEQEQk5FQIRkZBTIRARCTmL33A2fpjZXmDrCGevAfalMU66ZGsuyN5syjU8yjU8uZjrUHefMtiIcVcIRsPMVrv7wkznGChbc0H2ZlOu4VGu4QlbLh0aEhEJORUCEZGQC1shWJrpAEPI1lyQvdmUa3iUa3hClStU5whEROSdwrZHICIiA6gQiIiEXE4XAjP7ppmtNbPnzexRM5sxxHTnmNlGM9tkZtePQa7vmNnLiWwPmFnVENNtMbMXE/lXZ1GusV5ff2Vm680sZmZDXjo31utrmNnGep1NMrPHzOzVxL/VQ0wX+Do72Ge3uO8nxq81s5OCyDGCXHVm1pRYN8+b2dfGKNedZrbHzNYNMT7962uoR5flwg8woV/774DbBpkmH3gNOAwoAl4Ajg441weBgkT728C3h5huC1AzhuvroLkytL6OAuYD9cDCJNON6fpKNVuG1tm/A9cn2tdn6m8slc8OnAc8AhhwCrByDH5vqeSqAx4ey7+nxPueAZwErBtifNrXV07vEbh7c7/BcmCwM+OLgU3uvtndu4H7gAsCzvWou0cSg88As4J8v1SlmCsT6+sld98Y5HuMVIrZxnydJZb/00T7p8BHAxU2cJUAAAWHSURBVH6/oaTy2S8A7va4Z4AqM5ueBbkywt2fAhqSTJL29ZXThQDAzL5lZm8AlwKD7drNBN7oN7w98dpY+Rvi1X0wDjxqZs+Z2dVjmAmGzpXp9ZVMJtdXMplYZ7Xuvgsg8e/UIaYLep2l8tkzsX5Sfc8lZvaCmT1iZscEnClVaV9fBaOKkwXM7PfAtEFG3ejuD7r7jcCNZnYDcB3w9YGLGGTeUV9Te7BciWluBCLAPUMs5jR332lmU4HHzOzlxLeFTObK2PpKQdrXV5qyjfk6G8ZiAlln/aTy2QNZPweRynuuId4/T6uZnQf8CpgXcK5UpH19jftC4O7vT3HSnwG/4Z2FYDswu9/wLGBn0LnM7Argw8BZnjjwN8gydib+3WNmDxDfnR3Vf9I05MrI+kpxGWlfX2nKNubrzMx2m9l0d9+VOGywZ4hlBLLO+knlsweyfkabq/+hZXdfbmY/MrMad890Z3RpX185fWjIzPpX7/OBlweZbBUwz8zmmlkRcDHwUMC5zgG+Apzv7u1DTFNuZpW9beIncge9imAsc5GB9ZWKTKyvYcjEOnsIuCLRvgJ4x57LGK2zVD77Q8DliathTgGaeg9rBeigucxsmplZor2Y+PZyf8C5UpH+9TXWZ8TH8gf4JfE/7LXAr4GZiddnAMv7TXce8ArxqwhuHINcm4gf43s+8XPbwFzEr2Z4IfGzPltyZWh9fYz4t6AuYDfwu2xYX6lmy9A6mww8Drya+HdSptbZYJ8duAa4JtE24JbE+BdJcmXYGOe6LrFeXiB+8cSpY5TrXmAX0JP42/pM0OtLXUyIiIRcTh8aEhGRg1MhEBEJORUCEZGQUyEQEQk5FQIRkZBTIZCcY2ato5j3ukSvjm5mNf1eH7LHRzMrNbMnzSx/tLmGm93MPmxmNw1nHpGBVAhE3u7PwPuBrQNeP5d49wLzgKuBW/uN+xtgmbtHxyTh2/0GON/MyjLw3pIjVAgkZyW+xX/HzNZZvM/9ixKv5yW6C1hvZg+b2XIzuxDA3f/i7lsGWVyyHh8vJXHnrplVmNnjZrYm8Z7v6NHS4v3cP2XxZz5sMLPbzCyv3/hvJTo6e8bMahOvfcTMVprZX8zs972ve/xGoHri3YKIjIgKgeSyjwMnAMcT/5b/ncTG++PAHOBY4LPAkhSWNWiPj4nuCQ7rVzw6gY+5+0nAe4Hv9nZTMMBi4EuJDIcnMkG8u/Rn3P144n3+XJV4/U/AKe5+IvEuk/+537JWA+9J4TOIDGrcdzonksTpwL2JQza7zexJYFHi9fvdPQa8aWZPpLCsoXp8rAEaB0z3r2Z2BhAjXkBqgTcHzPusu28GMLN7E5n+F+gGHk5M8xzwgUR7FvDzRCErAl7vt6w9xLuOEBkR7RFILhts453s9WSG6vGxAyjp9/qlwBRggbufQLz/of7jew3s26V3uMff6vclyltf1n4A/NDdjwU+N2CZJYkcIiOiQiC57CngIjPLN7MpxB8B+CzxwyyfSJwrqCX+SMKDGbTHR3c/AOSbWe+GeSKwx917zOy9wKFDLG9xoufLPOCiRKZkJgI7Eu0rBow7guzpaVXGIRUCyWUPEO959gXgD8A/u/ubxHul3U5843k7sBJoAjCzvzOz7cS/8a81szsSy1oObCbeQ+uPgWv7vc+jxA/tQPxhPgst/iD4Sxm863OAFcC/JTK8nsiazDeA+83sj8DA/vDfS/zqIZERUe+jEkpmVuHxJ09NJr6XcFqiSIxkWScCX3T3y1Kcvg74J3cf9ZU+iT2an7n7WaNdloSXThZLWD1sZlXET7x+c6RFAOKXnJrZE2aWn4F7CQ4hfvWRyIhpj0BEJOR0jkBEJORUCEREQk6FQEQk5FQIRERCToVARCTk/j9wfIApWYxJagAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Compute the mean and standard deviation over the different folds.\n", "mse_mean = np.mean(mse,axis=1)\n", "mse_se = np.std(mse,axis=1) / np.sqrt(nfold-1)\n", "\n", "# Plot the mean MSE and the mean MSE with 1 SE error bars\n", "plt.errorbar(np.log10(alphas), mse_mean, yerr=mse_se)\n", "plt.xlabel('log10(alpha)')\n", "plt.ylabel('Test MSE')\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We find the optimal `alpha`, by the following steps:\n", "* Find the `alpha` with the minimum test MSE\n", "* Set `mse_tgt = ` minimum MSE + 1 std dev MSE\n", "* Find the least complex model (highest `alpha`) such that `MSE < mse_tgt`" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimal alpha = 0.183074\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU5fn/8fedhSwQSAAJskNVVlkkoAhicAWqVqstRb8uRQti/drWVqu1YKuldWvrioDLT+3XIlrXWlCLEhEFBGQREBABIYLsAbKRZOb+/XEmMctkMknmZCYz9+u6cuWZOc8585lDmHvO9hxRVYwxxsSuuHAHMMYYE15WCIwxJsZZITDGmBhnhcAYY2KcFQJjjIlxVgiMMSbGuVYIRORZEdknIusD9MkWkTUiskFEPnQrizHGmNqJW9cRiMhoIB94QVUH+JmeDnwCjFXVnSLSQVX3uRLGGGNMrVzbIlDVxcChAF2uBF5T1Z2+/lYEjDEmDBLC+NqnAIkikgOkAY+o6gv+OorIZGAyQEpKytCuXbs26AW9Xi9xcZF3WCRSc0HkZrNc9WO56icac23ZsuWAqp7gd6KquvYD9ADW1zLtcWAZ0BJoD3wJnFLXMocOHaoNtWjRogbP66ZIzaUaudksV/1YrvqJxlzASq3lczWcWwS5wAFVLQAKRGQxMAjYEsZMxhgTc8K57fMmcJaIJIhIKnA68EUY8xhjTExybYtAROYC2UB7EckF7gYSAVR1lqp+ISLvAOsAL/C0qtZ6qqkxxhh3uFYIVHViEH0eBB5s7GuVlpaSm5tLcXFxwH5t2rThiy8ib6MjUnNB1WzJycl06dKFxMTEMKcyxoRSOI8RhExubi5paWn06NEDEam137Fjx0hLS2vCZMGJ1FzwXTZV5eDBg+Tm5tKzZ89wxzLGhFDknR/VAMXFxbRr1y5gETCNIyK0a9euzq0uY0zzExWFALAi0ARsHRsTnaKmEBhjjGkYKwQhIiJcffXVFY/Lyso44YQTuOiiiwDYu3cvF110EYMGDaJfv36MHz8egB07dtChQwcGDx5c8fPCC34vsDbGGFdExcHieunYEfburfl8ZiZ8+22DF9uyZUvWr19PUVERKSkp/Pe//6Vz584V06dPn87555/PL37xCwDWrVtXMa1nz56sWbOmwa9tjDGNEXtbBP6KQKDn62HcuHH85z//AWDu3LlMnPjdGbR79uyhS5cuFY8HDhzY6NczxphQiM5CkJ1d82fmzODmPXCg5rxB+slPfsJLL71EcXEx69at4/TTT6+Y9vOf/5zrr7+eMWPGMGPGDHbv3l0xbfv27VV2DX300UdBv6YxxjRW7O0actHAgQPZsWMHc+fOrTgGUO7CCy9k27ZtvPPOOyxYsIAhQ4awfr1zIbXtGjLGhFN0FoKcHP/PHztW97zt29c+fxAuueQSfvOb35CTk8PBgwerTGvbti1XXnklV155JRdddBGLFy9m6NChDX4tY4wJhejcNRRGkyZNYvr06Zx66qlVnv/ggw8oLCwEnKt1v/rqK7p16xaOiMYYU0XsFYLMzPo9X09dunSpODOoslWrVpGVlcXAgQMZMWIEN9xwA8OGDQNqHiN49NFHQ5LFGGOCEZ27hgJpxCmigeTn59d4Ljs7m2zfwebbbruN2267rUafHj16sG/fvogda8gYE/1ib4vAGGNMFVYIjDEmxlkhMMaYGGeFwBhjYpwVAmOMiXFWCIwxJsZZIQgRt4ehfvjhhysuSHPTmjVrmD9/vuuvY4yJHLFZCF58EXr0gLg45/eLLzZ6kZWHoQZqHYZ67dq1bNy4kfvuu69iWvlYQ+U/11xzTY3lN6QQeDyeer8PKwTGxJ7YKwQvvgiTJ8PXX4Oq83vy5JAUA7eGoX700UfZvXs3Y8aMYcyYMQBMnTqVrKws+vfvz913313Rt0ePHtxzzz2MGjWKV155hfnz59OnTx9GjRrFLbfcUrGFUlBQwKRJkxg2bBhDhgzhzTffpKSkhOnTpzNv3jwGDx7MvHnzGrU+jDHNQ3QWgkDDUN95J1T/Zl1YCOXDQkTgMNS33HILnTp1YtGiRSxatAiAGTNmsHLlStatW8eHH35Y5UY3ycnJLFmyhEsvvZQpU6awYMEClixZwv79+yv6zJgxg3POOYcVK1awaNEibrvtNkpLS7nnnnuYMGECa9asYcKECUG/d2NM8xV7Q0zk5vp/vtpIoQ3RlMNQv/zyy8yZM4eysjL27NnDxo0bK7Yyyj/AN23aRK9evejZsycAEydOZM6cOQC89957vPXWWzz00EMAFBcXs3Pnzoa/eWNMsxWdhSDQMNTdujm7g6rr3t353QyGod6+fTsPPfQQK1asICMjg+uuu47i4uKK6S1btgRAVWtdhqry6quv0rt37yrPL1++vN55jDHNW3TuGgpkxgxITa36XGqq83wIuDUMdVpaGsd891M4evQoLVu2pE2bNuzdu5cFCxb4nadPnz5s27aNHTt2AFTZ53/hhRfy2GOPVRSL1atX13gdY0xscK0QiMizIrJPRNbX0W+YiHhE5Aq3slRx1VUwZ46zBSDi/J4zx3k+BNwahnry5MmMGzeOMWPGMGjQIIYMGUL//v2ZNGkSI0eO9JslJSWFmTNnMnbsWEaNGkVmZiZt2rQBYNq0aZSWljJw4EAGDBjAtGnTABgzZgwbN260g8XGxBJVdeUHGA2cBqwP0Cce+ACYD1wRzHKHDh2q1W3cuLHGc/4cPXo0qH5Nzc1cx44dU1VVr9erU6dO1b/97W/1mr96tmDXtdsWLVoU7gh+Wa76sVz105hcwEqt5XPVtS0CVV0MHKqj2/8CrwL73MoR65566ikGDx5M//79OXLkCFOmTAl3JGNMhBENcECx0QsX6QG8raoD/EzrDPwTOAd4xtfvX7UsZzIwGSAzM3PoSy+9VGV6mzZtOOmkk+rM4/F4iI+Pr9+baAKRmgtqZtu6dStHjhwJYyJHfn4+rVq1CneMGixX/Viu+mlMrjFjxqxS1Sy/E2vbVAjFD9CDWnYNAa8AZ/jaz2G7hiKS7RqqH8tVP5arftzaNRTO00ezgJdEBKA9MF5EylT1jTBmMsaYmBO2QqCqPcvbIvIczq4hKwLGGNPEXCsEIjIXyAbai0gucDeQCKCqs9x63WBNmL0UgHlTRoQ5iTHGhJebZw1NVNUTVTVRVbuo6jOqOstfEVDV67SWA8XNXY8ePThw4EDAPgMGDKizTyj9+c9/brLXMsZEvti7sthYITDGVGGFIEQuvfRShg4dSv/+/SsGdqtsx44d9OnTh2uvvZaBAwdyxRVXVLm/wGOPPcZpp53GqaeeyqZNmwD49NNPOfPMMxkyZAhnnnkmmzdvrrFcVeW2225jwIABnHrqqRVXA+fk5DB69Gguu+wy+vXrx4033ojX6+WOO+6gqKiIwYMHc1WIrqY2xjRvUTfo3B//vYGNu4/6nVb5nPiNe5w+5ccKAunXqTV3X9w/YJ9nn32Wtm3bUlRUxLBhw7j88stp165dlT6bN2/mmWeeYeTIkUyaNImZM2fym9/8BoD27dvz2WefMXPmTB566CGefvpp+vTpw+LFi0lISGDhwoX87ne/49VXX62yzNdee401a9awdu1aDhw4wLBhwxg9ejTgFJKNGzfSvXt3xo4dy2uvvcZ9993H448/Xu/RTo0x0cu2CELk0UcfZdCgQZxxxhns2rWLL7/8skafrl27VowL9D//8z8sWbKkYtoPf/hDAIYOHVoxSNyRI0f40Y9+xIABA/jVr37Fhg0baixzyZIlTJw4kfj4eDIzMzn77LNZsWIFAMOHD6dXr17Ex8czceLEKq9njDHlom6LINA392PHjpGWlgaE9qyhnJwcFi5cyNKlS0lNTSU7O7vKsNDlfNdM+H2clJQEQHx8PGVlZYAzMNyYMWN4/fXX2bFjB9l+bpKjAa4MD/R6xhhTzrYIQuDIkSNkZGSQmprKpk2bWLZsmd9+O3fuZOlSpwDNnTuXUaNG1bnc8vseP/fcc377jB49mnnz5uHxeNi/fz+LFy9m+PDhgLNraPv27Xi9XubNm1fxeomJiZSWljbkrRpjopAVghAYO3YsZWVlDBw4kGnTpnHGGWf47de3b1+ef/55Bg4cyKFDh5g6dWrA5d5+++3ceeedjBw5stYb0V922WUMHDiQQYMGcc455/DAAw/QsWNHAEaMGMEdd9zBgAED6NmzJ5dddhngDGk9cOBAO1hsjAGicNdQOCQlJdV6c5jy/f35+fnExcUxa1bNa+nWr19fscsqKyuLHN8d0kaMGMGWLVsq+t1777015hURHnzwQR588MEa01JTU/3eU+D+++/n/vvvr/N9GWNiQ8wWArui2BhjHDFbCJpajx49Km5W3xSys7P9Hlw2xpjqouYYQaCzZ0xo2Do2JjpFRSFITk7m4MGD9kHlIlXl4MGDJCcnhzuKMSbEomLXUJcuXcjNzWX//v0B+xUXF0fkB1mk5oKq2ZKTk+nSpUuYExljQi0qCkFiYiI9e/ass19OTg5DhgxpgkT1E6m5ILKzGWNCIyp2DRljjGk4KwTGGBPjrBAYY0yMs0JgjDExzgqBMcbEOCsExhgT46wQGGNMjLNCYIwxMc4KgTHGxDgrBMYYE+OsEBhjTIyzQmCMMTHOtUIgIs+KyD4R8Xs3FhG5SkTW+X4+EZFBbmUxxhhTOze3CJ4DxgaYvh04W1UHAvcCc1zMYoyJIBNmL2XC7KWutRs7/1+WF7mar6FZ/7K8qJ5rOjiuDUOtqotFpEeA6Z9UergMsIHujYli5R9odr/wyBMp9yO4HlhQ20QRmQxMBsjMzCQnJ6dBL5Kfn9/ged0UqbkgcrNZrvqJhFx5ec632ZycnIp2fr7H7/ONbdf2esG2PR4PeXl5Ic0Uiqwej8eVf8ewFwIRGYNTCEbV1kdV5+DbdZSVlaUNvSl7Tk5ORN7QPVJzQeRms1z1Ewm5ntzsbBFkZ4+oaLdqdZz09KQazze2XdvrBdv+y/IFpKenhzRTKLLm5eW58u8Y1kIgIgOBp4FxqnownFmMMSZWhe30URHpBrwGXK2qW8KVwxhjYp1rWwQiMhfIBtqLSC5wN5AIoKqzgOlAO2CmiACUqWqWW3mMMcb45+ZZQxPrmH4DcINbr2+MMSY4dmWxMcY11c+XN5HJCoExxsQ4KwTGGBPjrBAYY0yMs0JgjDExzgqBMcbEOCsExhgT46wQGGNMjLNCYIwJKbt2oPmxQmCMMTHOCoExxsQ4KwTGGBPjrBAYY0yMq7UQiEjrANO6uRPHGGNMUwu0RZBT3hCR96tNe8OVNMYYY5pcoEIgldptA0wzxsQ4O2W0eQtUCLSWtr/HxhhjmqlAdyjrICK34nz7L2/je3yC68mMMcY0iUCF4CkgzU8b4GnXEhljjGlStRYCVf1jUwYxxhgTHrUWAhH5GZCjql+KiADPAJcDXwPXqurqJspojIkyqooqHC9TyrzOIcdjxaUhawONmt+jTq5QZgpFVq+6c3hWtJYFi8h6YIiqlorIlcCvgQuAIcDdqnqWK4nqkJWVpStXrmzQvDk5OWRnZ4c2UAhEai6I3GyWq37czlV+xtC8KSMCtks9Xs7tm8nDC7dQ6rFzTuqrbRJ89sfvN2heEVmlqln+pgU6RlCmqqW+9kXAC6p6EFgoIg80KIkxJmpU/pAPxq7Dhew5UsxnO/NonZxAx9aJjOhQxrL9zsfQNSN68MLSHSFpA42af/YHm0hJSQ5pplBk1ZJi3BDo9FGviJwoIsnAucDCStNSXEljjIlK63Lz2J1XTHpKIu/9ajR9T2xNp/QUxvdswYltUjixTQo3nNUrZO3Gzt82WUKeKRRZUxPduYQr0BbBdGAlEA+8paobAETkbGCbK2mMMVHpiUVbiY8Tep3QilMy0+qewTSpQGcNvS0i3YE0VT1cadJKYILryYwxUaGwpIzl2w/ROT2ZhDgblCASBTpr6IeV2v66vBZowSLyLM6xhX2qOsDPdAEeAcYDhcB1qvpZcLGNMc3F7rxiUlvEk9k6OdxRTC0C7Rr6F7DG9wNVxxdS6igEwHPA48ALtUwfB5zs+zkdeNL32xgTJYpLPRwsKGHy6F6s3ZUX7jimFoEOFl8ObAEGAtuBGar6U9/PpLoWrKqLgUMBuvwA50wkVdVlQLqInFiP7MaYCLc7rxgRuOGsnuGOYgIIdIzgdeB1EWmJ86H9VxFpB9ylqh+G4LU7A7sqPc71PbenekcRmQxMBsjMzCQnJ6dBL5ifn9/ged0UqbkgcrNZrvpxI1deXhHgXKPgr/3GOx+wP/846UmwcdUyv33y8z11Lqch7WDyBWp7PB7y8vJCmikUWT0ejyt/X4F2DZUrBo4AR4FuQKh29Pk78OD3ChNVnQPMAeeCsoZeGBOrF/s0RqRms1z140auJzc71xFkZ4/w294W1xbYSvcT2pCdPcpvn1atjpOenhRwOQ1pB5MvUPsvyxeQnp4e0kyhyJqXl+fK31egg8VjgInAcJxrCB5R1YZd0utfLtC10uMuwO4QLt8YEyZer/J/y3eSkZpIcmJ8uOOYOgTaIngfWAcsAZKAa0TkmvKJqnpLI1/7LeBmEXkJ5yDxEVWtsVvIGBMZ6nMl8YGC4xwqKKFvR7tmoDkIVAh+2pgFi8hcIBtoLyK5wN1AIoCqzgLm45w6uhXn9NFGvZ4xJjKoKt8eOU7fE1uTlmRbA81BoIPFzzdmwao6sY7pCvy8Ma9hjIk8R4vLKCr1cP2onryyclfdM5iwC3T6qDHG1Nu3R4pJjBcuHmRngzcXVgiMMSGz/psj5BWV0iEtmaQE2y3UXNRZCERkZDDPGWPMg+9uJj5O6Ng6KdxRTD0Es0XwWJDPGWNi2NGiUj7cst8ZXC7edjY0J4GuIxgBnAmcICK3VprUGmdoamOMAZwzhXYdLqRj62Qy02xwueYmUNluAbTCKRZplX6OAle4H80Y01wcLiwl/7iHX553MnE21HSzE+j00Q+BD0XkOVX9GkBE4oBWqnq0qQIaYyKbx6vkHi4iOTGOK4Z24fXV34Q7kqmnYHbk/UVEWvsGn9sIbBaR21zOZYxpJt5Y/Q1FpR66ZKTasYFmKph/tX6+LYBLca4G7gZc7WoqY0yz4FXl7wu3kNoinrapieGOYxoomEKQKCKJOIXgTVUtpZZRQo0x0WXC7KUVYwz5s//YcXIPF9E1I6W2OxmaZiCYQjAb2AG0BBb77mNsxwiMiXEer/JNXhHDemTQJsW2BpqzOguBqj6qqp1VdbzvbmJfA2OaIJsxJoLtO1ZMqUe57cI+tjXQzAVzZXGmiDwjIgt8j/sB17qezBgTsfKPl7E7r5g2KYkM79k23HFMIwWza+g54F2gk+/xFuCXbgUyxkS+d9d/S5lX6ZxuF49Fg1oLgYiUX2PQXlVfBrwAqloGeJogmzEmQi1Yv4cW8XG0Sgrmbrcm0gXaIvjU97vAd9N6BRCRM3DuYWyMiUHHiktZvOUAbVsm2rGBKBGonJf/C9+Kc1vJ74nIx8AJ2BATxsSsDzbto8TjpW3LFuGOYkJEnBuF+Zng3F7yb76HcTj3LRbgOOBR1b/5ndFlWVlZunLlygbNm5OTQ3Z2dmgDhUCk5oLIzWa56qehuSrfp7i8nZHags92HqZHu1REpMq0YNuNzeW2aMwlIqtUNcvftEBbBPE4g85V3/ZLbVAKY0yz5/EqOVv2MSGrK5u+PVaveYO56b0Jj0CFYI+q3tNkSYwxEcHfN/dyeUWlFJd6GTvgxHoXAhO5Ah0stqNAxpgqDhWU0K5lC7t2IMoEKgTnNlkKY0zE83qVvMISLhzQkXi750BUqbUQqOqhpgxijIlsR4pK8SqMG9Ax3FFMiNng4caYoBwuLCU+TjijV7twRzEhZpcFGmPq5PUqhwtLSE9JJLEeN5+xM4WaB1e3CERkrIhsFpGtInKHn+ltROTfIrJWRDaIyE/dzGOMaZg1uXmUeZV0u/lMVHKtEIhIPPAEMA7oB0z0jVxa2c+Bjao6CMgG/ioidrmiMRHm/S/2ApBu9x2ISm5uEQwHtqrqNlUtAV4CflCtjwJp4gxY0go4BJS5mMkY40dddyJbuHEfrZMT7J7EUcrNf9XOwK5Kj3N9z1X2ONAX2A18DvxCVb0uZjLG1FNxqYfNe4+Rnmob69HKzYPF/k40rj6w0YXAGuAc4HvAf0XkI1WtcitMEZkMTAbIzMwkJyenQYHy8/MbPK+bIjUXRG42y1U/deXKyysCnLFsqrcPFzv/bePLisjLK65YTqB5ytuNzRUusZbLzUKQC3St9LgLzjf/yn4K3KfOyHdbRWQ70IfvhsAGQFXnAHPAGXSuoYMuReNAUm6L1GyWq37qyvXkZme3UHb2iBrtPXuOckpmMhm+LYLs7BEB56nP22+u6ytc3Mrl5q6hFcDJItLTdwD4JzjDWVe2E98VzCKSCfQGtrmYyRhTD2UeL0eLyzivb2a4oxgXuVYIfHcyuxnnNpdfAC+r6gYRuVFEbvR1uxc4U0Q+B94HfquqB9zKZIypn7yiUgDOtUIQ1Vy9oExV5wPzqz03q1J7N3CBmxmMMQ13IL+EFglxDOmaHu4oxkV2LpgxMaquU0b3HS3mSFEp7Vu1IM4GmYtqNsSEMcavN9c453a0b5VUZ18bSqJ5sy0CY4xfr36WS8ukeFIS48MdxbjMCoExpobCkjI2fXssqK0B0/xZITDG1HAgv4SEOKFdS7uaOBZYITDGVKGqHMg/zpg+Heo15LRpvuxgsTExJNCN6csdKSql1KP8cEhnnvtkR6397ABx9LByb4ypYu+x4yTECef07RDuKKaJWCEwxlTYfqCAvMJSMlsnkZRgZwvFCisExpgKz328HQE6pCWHO4ppQnaMwJgoN2H2UvLyiuocFbTM4+WVVbm0a9WCFgn2HTGW2L+2MQaAfceOU1jioWNr2xqINbZFYIxBVdl79Dhn9GqLVr99lIl6tkVgjOFQQQklHi+TRvYMdxQTBrZFYEyMK/N4yc0rIjkxjnP7ZvLMku1++9l1A9HLtgiMiUJ1DTFd2WuffUNxqZeuGanE23DTMckKgTExzOtV/r5wCy2T4slITQx3HBMmVgiMiWF7jxWz50gx3TJSEbGtgVhlxwiMiRLBjCNUWZnHyzd5xZx9ygkUl3rcjGYinG0RGBOjvskrxuNVbh/bO9xRTJhZITAmBi3fdpBvjxbTIS2J/p3ahDuOCTMrBMbEmDKvcuvLa0lKiKNb29RwxzERwI4RGNOM1fe4AMDOgwUcLCihT8fWdZ4uatcOxAbbIjAmhhwqKGF/fgk3ZZ9EWrJ9DzQOKwTGxIgte4/x1f58WraI55ZzTw53HBNBrBAY08zU56rhcmVe5frnVxAfJ5yc2cqGmTZVuPrXICJjRWSziGwVkTtq6ZMtImtEZIOIfOhmHmNikVeV3QXO6KKndEizO4+ZGlwrBCISDzwBjAP6ARNFpF+1PunATOASVe0P/MitPMbEIlVlx4FCisrgwSsG0sqOCxg/3NwiGA5sVdVtqloCvAT8oFqfK4HXVHUngKruczGPMc1WQ3YHAcxZvI39+cdplww/GNzZhWQmGrj59aAzsKvS41zg9Gp9TgESRSQHSAMeUdUXqi9IRCYDkwEyMzPJyclpUKD8/PwGz+umSM0FkZst1nLl5RUBkJOTE3T7WIly3/ZDpCVCeqIGPe9U34XGTbF+Y+3fsbHcyuVmIfB3gnL1ex8lAEOBc4EUYKmILFPVLVVmUp0DzAHIysrS7LpuvlqLnJwcGjqvmyI1F0RutljL9eRmZ2sgO3tEUO2H1n3Et3lHGdQ1ncQ44ejRI2RnZwc1b1OKtX/HxnIrl5u7hnKBrpUedwF2++nzjqoWqOoBYDEwyMVMxjQbDd0dlHu4kM17j5EQF8dT12QRZ/cYMHVwsxCsAE4WkZ4i0gL4CfBWtT5vAmeJSIKIpOLsOvrCxUzGRLVSj5drnv0Ur0Lvjq04IS0p3JFMM+DariFVLRORm4F3gXjgWVXdICI3+qbPUtUvROQdYB3gBZ5W1fVuZTIm0jVkyIhyHq+yZW8+JR4vp2S2IrWFnSFkguPqX4qqzgfmV3tuVrXHDwIPupnDmGhXUuZl67588o+XMet/TuP/fbwj3JFMM2JfGYwJo8ZsAZTzqnLL3NXkFZXSo10qYwecWO9CYIPLxTYrBMY0Y6rKV/sLWLHjMN3bppLZOjnckUwzZAOOGNPEGno2UHVlHi9f7S/gUEEJd43vS8c2VgRMw1ghMKYJhOrDv5yqc3OZgwUldM1I4Weje4Vs2Sb22K4hY1wSiv3//qgqW/cX8OmOw3TNSKFTekqDlmPHBUw52yIwppEqf9v/y/KikH7zr6641MPWffkcKijhjnF9GlwEjKnMCoExAVT+kK+t3VRKyrxMmLOMQ4WldGubwo1nf69JX99Er+jfNdSxI+zdC0B25eczM+Hbb8ORyIRR5d01wbTDYfbtF8GNh5hX/sSN8NtOvZnyw99TkNGekzu0om3LFmHJZqJT9BcCXxEI+nnTLEXqh3pDpB89VNEukzhmnXEFj4ycSMdjB/nHTWdy95sbGrzs5rQeTNOJ/kIQSHY2/PjHcNNNUFgI48fX7HPddc7PgQNwxRU1p0+dChMmwK5dcPXVNaf/+tdw8cWweTNMmVJz+u9/DwkJsGYN/PKXNaf/+c9w5pnwySfwu9/VnP7wwzB4MCxcCH/6U83ps2dD797w73/DX/9ac/o//gFdu8K8efDkkzUmJ/7qV07jueecn+rmz4fUVJg5E15+ueb08iFzH3oI3n676rSUFFiwwGnfey+8/37V6e3awauvOu0774SlS9mw+ygAg1O9fJSUyeOT/sC8KSO49uWH6b7rS5jbmum+PqzKgqE/ddqTJ8OWLd9Nm9uaa+Mzef7Hzjq/+dk/0Pbw/qrz77gQelzitC+/nOmfb6+Yd/ruo6zvkwW+D9Y7HrsV5ibxwM5DJCQkwNzWXNTuVN6+4EoApv/151XmBcBzA8QPoUVJsfO3CN9NAza3785vxv+Sz088me9/sZg/vfckGbOur7mOjWkkO/7fm4UAABKLSURBVEaA823y6meWA7Bh99GKD5sNu48yc9FWAG54fkWV58vbDy/cUvGts/q85e0Js5fyy5dW++3zp7c3AnD7v9b6nbd8/mlvrPc7/+3/WgvAn97e6Hf+X760Oqh8Dy/c4rfPw6uccepnLtoaMN8zH23zO3/5a/9j6de1vnZt+ZZvP1TR543V3/j7p4tKZRLHzNOv4OJrH2Z36xOY+cZfeOKtB8goPhbuaCZKxfYWQfm31dlLKWmRDDk53FNpV0J5+ybgWKt07vn1E1WenzdlBEvLDxh27co9v36ixrzzLh4Bs5eyp2P3Wpc/iuN83fUUv8sH4POlbPneqfCA//kBPu87jM/7Dqsx/57yPhdfzD2729fM17UrkMvSrPP45VPTaiz/aF4eAB+e+X0+PPP7/vMB72VfzvVTHqo139sXXMnbF1xZc/7yPtOmcU+H82rm880/97KpXFrp+am9j/Pk5u9G1iz/Zl/r8ufMAagy7flK+R6f9IfA87/6qt/3NsE3/33/+zfmTRnB7fcvID09nXlTRvB2peX7/dvwLb/8bw/gzseXEL98Oas792H8piXc+95M2hV9V3QbYt6UERF5kxUTOWK7EBgTIYpLPTz90TbWf3OEtLadeeStB7jki8V+7+5kTKhFfyHIzPR7YDivdVvSwxDHmMpUlcOFpZz/9w/ZdaiIjNREXnvqNnoezK3RtyC5JS3DkNFEv+gvBL5TRCfMXkpeXh7v/nbcd2eShDOXiWmqyptrvmH97qMUlng4JbMV/7zhdB55/0vumPEy81pt4+DNv6Jt3gGkezcePec6Pj79wqD/Zu3sIFMf0V8IjIkg2w8UkHu4iP35x/n0pTUkJ8bRs31L5t9yFgnxcTzy/pdOx6uu4qZ8Z/ygeVNG8HETX7xmYosVAmNcpKp8secou/OKOFRQwpiHcgBIS07giStP46nFXyEiJMTXPIFvxMqFTsO+3RuXWSEwJoRUleJSD6+s3MVX+/M5UlTKuEc+AiC1RTy/G9+H+Z/vISkhnvP7ZfK077Rbf87/8HVfa1pQr227g0xDWSEwphHyCks4WlRKQUkZU/9vFduOKGV5R1j7r3UkxAlpyQncOa4vcz/dSYuEOCaP/h7vf7Ev3LGNqcIKgTF1UFV2HSrkSFEpxaUe/vjvDWz69hiFJWUMvue/3/UDUhKhbVoqc67OYtobnyMi/HhYV179rOZZQKFgWwEmFKwQmJimquQeLuRoUSmlHi/PLtnOzkOFlJR5+fHspazdlUdJmZezHlhUMc9Ln+5CBFonJzJ5dC9e+yyXlkkJvHbTSC68fwHprZPp3TENEbsKwDQPVghMVFNVSj3Kih2H2HmwkNzDhRwv8zJh9lJW+z7kR93/3Yf8PW9vRIDEhDhUlVZJCSS1jOOWc0/muU92kJwYx2tTRzLxqWUATDn7e3ywyXb1mObNCoFp9o4UlVJwvIziMi9PLd7GjoMFHC/zcsHfP2Trvny8Cj+a9d3ply3i4/B4lTTfh/z/nnsyz328gxYJwguTTmfKP1YiIlVGL/3J8G687hvvKC6uab7p/33KDACerva87Q4yoWaFwEQ8VaXMq6zeeZgD+ccpLvUye20J6w8VU1zqZdAf36voO2P+F8QJJCXE071dS44WlZGUGMcfL+lP93Ytuf1fa4mr9iE/cXi3ikHtMlq2iJhdOsda2bXvpmlYIYhwqoqqAlDm8eL1tY+XefB6FQUKjpfh8SqKcqTQ2dcNsP/YcUrKvCjOwc7iUg8KbN3nHOhUhXW5eRwrLkNRPtl6gLzCEhRY8PkeDuQfp+C48uLyr/n2SDFeVR7/4EtyDxfiVbjn3xsp8XjYtj8fBW59eQ1f7c8H4NZ537V//fJavtqfjwB3vraO7QcKEOAPb23g64MFKHD3m+vZfqAAVfjFS6vZsvcYZV4l+8FF7DzkvN5lMz+pWC/FxUK8xNOuZQt+dlYvXl65i6SEOF684Qx+9sIKRISnrsmq+LDP7t0BgLgI+ZAPxtmf/MdpTBlhWwHGVTFTCA4VlLD1sJd+09+huNQDQL/p7wBUPO477R2OlwXf7jNtASVl3iptBXr/fgElvg/jU36/wPlgVjjlrgUVH9In/W6+78MbVuwArxYA0PPO/+D7rKfHHf+p8h5OumtBRbv379+paPe/+92K9qB7vvt2PGzGwop25YOd5/1tcUX7ksc/rmhf+fTyivbUFz+raN/1+vqK9kPvOcNVC/DKyl20SIjjaHEpgvDp9kMcKy4DYMXX37WXbTvI0SKn2Cz8Yh+HC0oAeH31N+QXl4HAm2t3k19chgis2ZVHcamXhDhhQOc2FJd6aZEQx/SL+vHI+1+SnBDHz/uWVIw++rPRvVj4hTOeVJvUxIj5Rt9YZy+d72v5uc+EMSHk6v0IRGSsiGwWka0ickeAfsNExCMifu78EhpJCXGkJ8FVp3ejQ1oyHdKSuer0bhWPM1snc/WI7mS2Dr597YgeNdodWydz3cgedPS1J43s6bTbJHP9WT3p2MZpTx7dixPbJNOpTTLf75lIJ1/75jEn0Tk9mc7pydxy7sn84tyT6ZyeQuf0FG49/xS6pKfQJSOF2y7sTdeMFLpmpHDnuD50a5tCt7apTL+oH93bptK9XSr3XjqAHu1S6dk+lQeuGEiv9i353gkteWziEE7q0IqTO7TimWuz6J3Zij4d05g3+Qz6nZhG/06tWfCLsxjYuQ09WwvLf3cuQ7ulk9U9gy9njGN4jwyG92zL53+8kFXTzue0bhkM6ZbOkt+ew+Cu6Qzums5Ht3/X/viOcxjSLZ3TumWw4q7zGNo9g6HdM1h79wVk9cggq3sGa6ZfwNDuGZzWLYMPbxvDwC5t6NepNY9feRrd26VyYptkzuuXSWqL+CbbRx9u/Tu1pn+n1uGOYWKAa1sEIhIPPAGcD+QCK0TkLVXd6Kff/cC7NZcSOi2TEuiQGsdd3+/HutwjANz1/X4AFY9/N74va3flBd2+c3xf1vhrj+vLmp1O+45xfVi98zAAvx3bh8++dtq3j+3DKl/78lOOs8/37fbXF/TmU98NWW49/xTA+UYNcMu5J/Px1gMA/HzMSSzesh+oeubKpFE9eXeDM9De1Wd05+21uwH4cVZXXl3lnMt+8aBO/N+yrwE4t28mcxY7V7ee3qsdacmJAPQ9sTUpLeI5Hi9ktk6uGAIhMT4uar5xG2Mcbu4aGg5sVdVtACLyEvADYGO1fv8LvAoMczGLMc2CHQsw4eBmIegM7Kr0OBc4vXIHEekMXAacQ4BCICKTgckAmZmZDbrbUl5eER6Ph5ycHPLynNsvli+n8uNwtPPzPbX2CXc+f+ssEvIFWmfhzOrxeMjLywu6/9TeVMlX2WDf3eHWhODuYvn5+RF5lzLLVT9u5XKzEPjbf6DVHj8M/FZVPYF2N6jqHGAOQFZWlmb7bvRdH09udu5HkJ2dzZObfWeSZI+omFb+OBztVq2Ok56e5LdPuPP5W2eRkC/QOgtn1r8sd25VGahP0H++nzhnSWWnpgY5Q+1ycnJoyP8bt1mu+nErl5uFIBfoWulxF2B3tT5ZwEu+ItAeGC8iZar6hou5jGlyDdrlE4ICYEww3CwEK4CTRaQn8A3wE+DKyh1UtWd5W0SeA962ImCiRaP398+c6fy+6abGhzEmANcKgaqWicjNOGcDxQPPquoGEbnRN32WW69tTLjceXpKxW6oRnv5Zee3FQLjMlcvKFPV+cD8as/5LQCqep2bWYwJpcrf9u1MH9PcxcyVxcY0RNg+8F98EZYtg+PHoUcPmDEDrrqq6V7fxBQrBMYQYd/wX3wRJk92igDA1187j8GKgXGFFQIT1ap/qEfUB35t7roLCgurPldY6DxvhcC4QMpHtmwuRGQ/8HUDZ28PHAhhnFCJ1FwQudmiNtdQGFrbtFWwqoGLjdr15ZJozNVdVU/wN6HZFYLGEJGVqpoV7hzVRWouiNxslqt+LFf9xFouV0cfNcYYE/msEBhjTIyLtUIwJ9wBahGpuSBys1mu+rFc9RNTuWLqGIExxpiaYm2LwBhjTDVWCIwxJsZFdSEQkXtFZJ2IrBGR90SkUy39grq3cghzPSgim3zZXheR9Fr67RCRz335V0ZQrqZeXz8SkQ0i4hWRWk+da+r1Vc9sTb3O2orIf0XkS9/vjFr6ub7O6nrv4njUN32diJzmRo4G5MoWkSO+dbNGRKY3Ua5nRWSfiKyvZXro15eqRu0P0LpS+xZglp8+8cBXQC+gBbAW6OdyrguABF/7fuD+WvrtANo34fqqM1eY1ldfoDeQA2QF6Nek6yvYbGFaZw8Ad/jad4TrbyyY9w6MBxbg3MzqDGB5E/y7BZMrG2do/Cb7e/K97mjgNGB9LdNDvr6ieotAVY9WetiSmndIg0r3VlbVEqD83spu5npPVct8D5fh3LQn7ILMFY719YWqbnbzNRoqyGxNvs58y3/e134euNTl16tNMO/9B8AL6lgGpIvIiRGQKyxUdTFwKECXkK+vqC4EACIyQ0R2AVcB/jbt/N1buXNTZPOZhFPd/VHgPRFZ5btvc1OqLVe411cg4VxfgYRjnWWq6h4A3+8OtfRze50F897DsX6Cfc0RIrJWRBaISH+XMwUr5Our2Q86JyILgY5+Jt2lqm+q6l3AXSJyJ3AzcHf1RfiZt9Hn1NaVy9fnLqAMeLGWxYxU1d0i0gH4r4hs8n1bCGeusK2vIIR8fYUoW5Ovs3osxpV1Vkkw792V9VOHYF7zM5zxefJFZDzwBnCyy7mCEfL11ewLgaqeF2TXfwL/oWYhCObeyiHPJSLXAhcB56pvx5+fZez2/d4nIq/jbM426j9pCHKFZX0FuYyQr68QZWvydSYie0XkRFXd49ttsK+WZbiyzioJ5r27sn4am6vyrmVVnS8iM0WkvaqGezC6kK+vqN41JCKVq/clwCY/3SrurSwiLXDurfyWy7nGAr8FLlHVwlr6tBSRtPI2zoFcv2cRNGUuwrC+ghGO9VUP4VhnbwHX+trXAjW2XJponQXz3t8CrvGdDXMGcKR8t5aL6swlIh1FRHzt4TiflwddzhWM0K+vpj4i3pQ/wKs4f9jrgH8DnX3PdwLmV+o3HtiCcxbBXU2QayvOPr41vp9Z1XPhnM2w1vezIVJyhWl9XYbzLeg4sBd4NxLWV7DZwrTO2gHvA1/6frcN1zrz996BG4EbfW0BnvBN/5wAZ4Y1ca6bfetlLc7JE2c2Ua65wB6g1Pe3db3b68uGmDDGmBgX1buGjDHG1M0KgTHGxDgrBMYYE+OsEBhjTIyzQmCMMTHOCoGJOiKS34h5b/aN6qgi0r7S87WO+CgiKSLyoYjENzZXfbOLyEUi8sf6zGNMdVYIjKnqY+A84Otqz4/DGV7gZGAy8GSlaZOA11TV0yQJq/oPcImIpIbhtU2UsEJgopbvW/yDIrJenDH3J/iej/MNF7BBRN4WkfkicgWAqq5W1R1+FhdoxMer8F25KyKtROR9EfnM95o1RrQUZ5z7xeLc82GjiMwSkbhK02f4BjpbJiKZvucuFpHlIrJaRBaWP6/OhUA5OMOCGNMgVghMNPshMBgYhPMt/0Hfh/cPgR7AqcANwIggluV3xEff8AS9KhWPYuAyVT0NGAP8tXyYgmqGA7/2ZfieLxM4w6UvU9VBOGP+/Mz3/BLgDFUdgjNk8u2VlrUSOCuI92CMX81+0DljAhgFzPXtstkrIh8Cw3zPv6KqXuBbEVkUxLJqG/GxPZBXrd+fRWQ04MUpIJnAt9Xm/VRVtwGIyFxfpn8BJcDbvj6rgPN97S7APF8hawFsr7SsfThDRxjTILZFYKKZvw/vQM8HUtuIj0VAcqXnrwJOAIaq6mCc8YcqTy9XfWyX8sel+t24Lx6++7L2GPC4qp4KTKm2zGRfDmMaxAqBiWaLgQkiEi8iJ+DcAvBTnN0sl/uOFWTi3JKwLn5HfFTVw0C8iJR/MLcB9qlqqYiMAbrXsrzhvpEv44AJvkyBtAG+8bWvrTbtFCJnpFXTDFkhMNHsdZyRZ9cCHwC3q+q3OKPS5uJ8eM4GlgNHAETkFhHJxfnGv05EnvYtaz6wDWeE1qeAmyq9zns4u3bAuZlPljg3gr8K/0OfAywF7vNl2O7LGsgfgFdE5COg+nj4Y3DOHjKmQWz0UROTRKSVOneeaoezlTDSVyQasqwhwK2qenWQ/bOB36hqo8/08W3R/FNVz23sskzssoPFJla9LSLpOAde721oEQDnlFMRWSQi8WG4lqAbztlHxjSYbREYY0yMs2MExhgT46wQGGNMjLNCYIwxMc4KgTHGxDgrBMYYE+P+P3AC5d5QkVp2AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Find the minimum MSE and MSE target\n", "imin = np.argmin(mse_mean)\n", "mse_tgt = mse_mean[imin] + mse_se[imin]\n", "alpha_min = alphas[imin]\n", "\n", "# Find the least complex model with mse_mean < mse_tgt\n", "I = np.where(mse_mean < mse_tgt)[0]\n", "iopt = I[-1]\n", "alpha_opt = alphas[iopt]\n", "print(\"Optimal alpha = %f\" % alpha_opt)\n", "\n", "# Plot the mean MSE and the mean MSE + 1 SE errorbars \n", "plt.errorbar(np.log10(alphas), mse_mean, yerr=mse_se)\n", "\n", "# Plot the MSE target\n", "amin_log = np.log10(alpha_min)\n", "aopt_log = np.log10(alpha_opt)\n", "plt.plot([amin_log,aopt_log], [mse_tgt,mse_tgt], 'rs--')\n", "\n", "# Plot the optimal alpha line\n", "plt.plot([aopt_log,aopt_log], [0.35,mse_mean[iopt]], 'ro--')\n", "\n", "plt.legend(['MSE', 'MSE target','alpha opt'],loc='upper left')\n", "plt.xlabel('log10(alpha)')\n", "plt.ylabel('Test MSE')\n", "plt.ylim([0.35,1.6])\n", "plt.grid()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we recompute the coefficients using all the training data at the correct alpha. We see that the model selects three non-zero cofficients: `lcavol`, `lweight` and `svi` (description of the features can be found in https://rafalab.github.io/pages/649/prostate.html). These features are presumably the most relevant in determining the PSA level. Interestingly, the first feature -- `lcavol` -- is the log of the cancer volume suggesting that the cancer volume does indeed influence the PSA level." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " lcavol 0.412157\n", " lweight 0.108583\n", " age 0.000000\n", " lbph 0.000000\n", " svi 0.090242\n", " lcp 0.000000\n", " gleason 0.000000\n", " noise 0.000000\n", "redundant 0.000000\n", "dependent 0.054286\n" ] } ], "source": [ "model.alpha = alpha_opt\n", "model.fit(X,y)\n", "\n", "# Print the coefficients\n", "for i, c in enumerate(model.coef_):\n", " print(\"%8s %f\" % (names_x[i], c))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the three chosen features, we now use linear regression method directly with cross validation to evaluate the test error and to determine the mean regression coefficients." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mse_mean=0.669781 mse SE=0.154929\n", "[0.51885655 0.24720759 0.22870017 0.02052603]\n" ] } ], "source": [ "# Find coefficients greater than a small threshold\n", "Isel = np.where(np.abs(model.coef_) > 1e-3)[0]\n", "\n", "\n", "# Select those features\n", "X1=X[:,Isel]\n", "nfea1 = len(Isel)\n", "\n", "# Run 10-fold validation with reduced feature matrix\n", "nfold = 10\n", "kf = sklearn.model_selection.KFold(n_splits=nfold,shuffle=True)\n", "\n", "# MSE for each alpha and fold value\n", "RSS_ts = np.zeros((nfold,1))\n", "coef=np.zeros((nfold,nfea1))\n", "for ifold, ind in enumerate(kf.split(X1)):\n", " \n", " \n", " # Get the training data in the split\n", " Itr,Its = ind\n", " X_tr = X1[Itr,:]\n", " y_tr = y[Itr]\n", " X_ts = X1[Its,:]\n", " y_ts = y[Its]\n", " \n", " regr.fit(X_tr,y_tr)\n", " y_ts_pred = regr.predict(X_ts)\n", " RSS_ts[ifold] = np.mean((y_ts_pred-y_ts)**2)/(np.std(y_ts)**2)\n", " coef[ifold]=regr.coef_\n", " \n", "mse_mean = np.mean(RSS_ts,axis=0)\n", "mse_se = np.std(RSS_ts,axis=0) / np.sqrt(nfold-1)\n", "coef_mean=np.mean(coef,axis=0)\n", "print(\"mse_mean=%f mse SE=%f\" % (mse_mean, mse_se))\n", "print(coef_mean)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## LASSO path\n", "\n", "To further illustrate the effect of regularization, we conclude by drawing the *LASSO path*. This is simply a plot of the coefficients as a function of the regularization `alpha`. We do not need to do this for the analysis, but the path demonstrates the effect of regularization well. \n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEOCAYAAACXX1DeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3yURd7Av7M9m91seggBEloCoQRIAEFKAAlFAbGAeJ5yniKeYjvvDu71PBFf9Tw9fS0nYoNDPBCxgIIFkaIgGCBID52EFkJI2SSbbfP+sSEEpISQZFPmy2c+8zwz88z8dng2v532+wkpJQqFQqFQXAyNvwVQKBQKRf1GKQqFQqFQXBKlKBQKhUJxSZSiUCgUCsUlUYpCoVAoFJdEKQqFQqFQXBK/KgohxHAhxG4hxF4hxNSLlEkVQmQIIbYLIVbVtYwKhULR1BH+OkchhNACmcBQIBv4GZggpdxRqUwwsBYYLqU8LISIlFLm+EVghUKhaKL4c0TRC9grpdwvpXQC84Ex55W5HfhESnkYQCkJhUKhqHv8qShigKxK99nlaZWJB0KEECuFEBuFEHfWmXQKhUKhAEDnx7bFBdLOnwfTAcnAECAAWCeE+ElKmfmryoSYBEwCCAgISG7ZsmW1hPJ6vWg0ao2/qqj+ujJqsr+OF3sBaBbYePtfvV9XxtX0V2ZmZq6UMuJCef5UFNlA5b/mLYCjFyiTK6UsBoqFEKuBJHxrG+cgpZwFzAJISUmR6enp1RJq5cqVpKamVuvZpojqryujJvtr/FvrAFhwX58aqa8+ot6vK+Nq+ksIcehief5UFD8D7YUQrYEjwG341iQq8znwuhBCBxiA3sDLdSqlQlFPmTK4vb9FUDQR/KYopJRuIcSDwNeAFnhPSrldCDG5PH+mlHKnEOIr4BfAC7wjpdzmL5kVivpEv/bh/hZB0UTw54gCKeVSYOl5aTPPu/8n8M+6lEuhaAhsP1oAQKfmNj9Lomjs+FVRKBSK6vP0Et+Ro8a0RuFyucjOzsbhcABgs9nYuXOnn6VqOFSlv0wmEy1atECv11e5XqUoFApFvSE7Oxur1UpcXBxCCIqKirBarf4Wq8Fwuf6SUnLq1Cmys7Np3bp1letV+84UCkW9weFwEBYWhhAX2j2vuFqEEISFhVWM2KqKGlFUYuaqfezd72Svdj86jUCv06DXajCWxwatBr3Od28oj406rS/WazDptZh0WvRaoV50haKaqO9O7VKd/lWKohIvf5tJmdsLmVc3J6oRYNJrCdBrfbFBi9nguzYbtAQadL7YqCPQ6IstRh1Wkw6LUY/FqCMoQEeQSU9QgB6rUYdGo748CkVdYLFYsNvttd7OxIkTueGGG7jllltqva2rRSmKSmyfPowVK1fR59p+uDwSl8eL0+3F6fGevS4PZR4vZS4vZW6P797txeHyVMSlTg8Ot4dSp5dSl5tSp4cSp4dTdidZzhJKnB6Ky9wUOz14vJc2zCgEBJn02AL0BJv1BJsNhJj1hJgNhAUaCAk0EG4xEG4xEm4xEhlkxGxQ/7WNnT8PT/C3CIomgvprUgmdVoNBK7Caqr4b4GqRUlLm9lLkcFNc5sZe5qbQ4aLI4aaw1EVBqasizi91kV/iIr/EycHcYk4XOykqc1+w3kCDlsggE5FWI9E2E81sAUTbTDQPDqB5sIkWwWZs5rr7nIqaJzk21N8iNHpeeOEF5s6di0ajYcSIETz//PO8/fbbzJo1C6fTSbt27Zg7dy4ul4ukpCT279+PRqOhpKSEhIQE9u/fz/bt25k8eTIlJSW0bduW9957j5CQEH9/tCtCKQo/I4TwrW3otURYjVf8fJnbw+liF7n2svLg5GRRGTlFDnKKysgpdLDx8GlOFBzH6fGe86zVpKNliJlWoWZiw83EhQUSFxZI28hAIixGNVdcz9l4KA9ovApj+pLtbM06jVarrbE6E5sH8fdRnapUdtmyZXz22WesX78es9lMXp6vv2+66SbuvfdeAJ544gneffddpkyZQlJSEqtWrWLQoEEsWbKEYcOGodfrufPOO3nttdcYOHAgTz75JNOnT+eVV16psc9UFyhF0cAx6rQ0s2lpZjNdspyUkly7k6P5pRzNLyX7dClZp0vIyishM6eIFbtyzlEkVpOOthEW4qMsxEdZSWhmpWN0EOGWK1dmitrhha92A43rHEV9Yvny5fzud7/DbDYDEBrqU8jbtm3jiSeeID8/H7vdzrBhwwAYP348CxYsYNCgQcyfP58//OEPFBQUkJ+fz8CBAwG46667uPXWW/3zga4CpSiaCEIIIqxGIqxGkloG/yrf45UczS/lQG4x+0/a2Z9bzJ4TdlbsyuGj9OyKcpFWI52aB9ElxkbXFsHYy7y/qkuhqAn+PqqTX89RSCkvOKqeOHEin332GUlJScyePZuVK1cCMHr0aKZNm0ZeXh4bN25k8ODBdbIoXhcoRaEAQKsRtAw10zLUzID4cy0Nn7KXsft4ETuOFbLjaCHbjxayKvMkZ9bgX8xYQY9WISTHhtCrdSgJUVa1S0vR4ElLS+Ppp5/m9ttvr5h6Cg0NpaioiOjoaFwuF/PmzSMmxudGx2Kx0KtXLx5++GFuuOEGtFotNpuNkJAQ1qxZQ//+/Zk7d27F6KIhoRSF4rKEWYz0bWekb7uzRuhKnG62Hy1k0fcbKTTYWH/gFIu3+KzEB5v19IoL5dp24fRrH06b8EC13qFocAwfPpyMjAxSUlIwGAyMHDmSZ599lhkzZtC7d29iY2Pp0qULRUVFFc+MHz+eW2+9tWKUATBnzpyKxew2bdrw/vvv++HTXB1KUSiqhdmgo2dcKMWt9aSmJiOlJPt0KRsO5LH+wCnW7T/FNztOABATHMCA+Aiu6xjJte3CMelrbnFSoahpKk8XTZ06lalTp56Tf//993P//fdf8NlbbrkFKc/d7t6tWzd++umnX5WdPXv21QtbRyhFoagRhDg7dXVzcgsADp0qZs2eXFZnnmRxxhH+u+EwJr2G/u0jGNmlGUM6RhFUh1uRGxtPjkr0twiKJoJSFIpaIzYskNiwQO64JpYyt4f1+/P4bucJvt5+gm93nECvFQyMj+DG7jFc1zFKjTSuEGVeXFFXKEWhqBOMOi0D4iMYEB/B30d1IiM7n6W/HGPJL0dZvjMHq1FH7zZhhAUaCDbrsZn1WE16gkxnTJnosAXosQUYsAXoMeiUPcsf9uQCyoGRovZRikJR52g0gh6tQujRKoRpIzvy0/5TfLr5CFuzC9h6JJ/TJS6c7ktvu7UYfYojzGIgxGwgNNBnziTM4tsC3DHaSkKUFZ228SqU11bsAZSiUNQ+SlEo/IpWI7i2XTjXVtpRJaXE4fJS5HBR6PCZNDljxqSgwoyJi9MlzvLgYn+unVN2JyVOT0U9Jr2GLjE2burRgrHdY9TUlkJRTfyqKIQQw4H/w+cz+x0p5fPn5acCnwMHypM+kVI+XadCKuocIQQBBp/V3cigK3u2xOnmeIGDbUcLyTicz9p9uUz7ZCsvfZPJ766N447escrGlUJxhfhtXC6E0AJvACOARGCCEOJC2zjWSCm7lQelJBSXxGzQ0SbCwuik5jw5KpFlD/fnw3t6k9g8iH9+vZs+z3/HjC92cCS/1N+iKuopFovlqusYOXIk+fn5lyyTmppKenr6r9IzMjJYunTpVctQk/hzRNEL2Cul3A8ghJgPjAF2+FEmRSNDCEHfduH0bRfOjqOFvL1mP7PXHmTO2oP0ax/O4A6RDEqIpGWo2d+iKhoRV/OHPiMjg/T0dEaOHFmDEl0d/lzpiwGyKt1nl6edTx8hxBYhxDIhRNXMPioUFyCxeRAvj+/Gqj+l8vt+rTmYW8yTn2+n/wvfM/yV1fzf8j1knij61YGp+sqzN3Xh2Zu6+FuMRssf/vAHFi9eDMDYsWO5++67AXj33Xd54oknAPjggw/o1asX3bp147777sPj8a2RxcXFkZvr25U2Y8YMOnTowNChQ5kwYQIvvvhiRRsLFy6kV69exMfHs2bNGpxOJ08++SQLFiygW7duLFiwoC4/8kXx54jiQjYdzv+GbgJipZR2IcRI4DOg/QUrE2ISMAkgKirqnCP0V4Ldbq/2s02RhtpffczQp6fgeHEAm3M8bDpRzCvLM3l5eSbNAwV9muu4JlpHhLlmf0vVRn9lXb5Ig8Fms1WYxDB+/3cCTmzHXYPWX7yRnSgbNP2y5YqKiujZsyffffcdgwYN4vDhw2RnZ1NUVMT333/PzTffTHp6OvPmzeOrr75Cr9fz6KOP8s4773D77bcjpcRut7N9+3YWLlzI6tWrcbvd9O/fn86dO1NUVITH46GkpITvvvuOr7/+mieffJLFixfz17/+lU2bNvHSSy9VyFJVPB5Plco7HI4reg/9qSiygZaV7lsARysXkFIWVrpeKoT4txAiXEqZe35lUspZwCyAlJQUmZqaWi2hVq5cSXWfbYo0hv66rTzOKXTw9fbjLN5ylEV7TrNoj4uU2BDGdGvOyC7RhNWAifWa7K/l5SZSrkuMqpH66gM7d+48ay1Wb8AtQKetwT9TegOGKlijtVqtpKWl8dZbb5GVlUWXLl04ffo0drud9PR03nzzTebMmcOWLVsYPHgwAKWlpbRo0QKr1YoQAovFwubNmxk7diyRkZEAjBkzBqPRiNVqRavVctttt2G1Wunfvz9Tp07FarViMpkwGAzVsppbVWu7JpOJ7t27V7lefyqKn4H2QojWwBF839fbKxcQQjQDTkgppRCiF76pslN1LqmiSRAZZOK3feL4bZ84svJKWLzlKJ9nHOFvn2/nqSU76NcunDHdmpPWqRkWo/93lr+9Zj/QuBTFOYx4nlI/mhmPiYnh9OnTfPXVVwwYMIC8vDw++ugjLBYLVqsVKSV33XUXzz333EXruNw0ptHo+/Gh1Wpxuy/srbI+4Lc1CimlG3gQ+BrYCXwkpdwuhJgshJhcXuwWYJsQYgvwKnCbbCgTyIoGTctQMw8Masc3jw7kq0f6M2lAG/bm2Hnsoy0kz/iWBz7cxKrMk5f1d65o2PTp04dXXnmFAQMG0L9/f1588UX69+8PwJAhQ/j444/JyckBIC8vj0OHDp3zfL9+/ViyZAkOhwO73c6XX3552TatVusVTTfVBX79WSSlXAosPS9tZqXr14HX61ouhaIyHZoF0WF4EH8elsCmw6f5POMoi7cc5ctfjhFtM3Frcgsm9G5FtC3A36Iqapj+/fvzzTff0K5dO2JjY8nLy6tQFImJiTzzzDOkpaXh9XrR6/W88cYbxMbGVjzfs2dPRo8eTVJSErGxsaSkpGCzXdpG16BBg3j++efp1q0b06ZNY/z48bX6GauCaIw/0FNSUuSF9idXhcYw516XNNX+KnN7WL4jhwXpWazZcxKNEKQlRnFnnziuaRN6Uf8bNdlf499aBzQuV6g7d+6kY8eOFff+9HBXU9jtdiwWCyUlJQwYMIBZs2bRo0ePWmmrqv11fj8DCCE2SilTLlTe/xOtCkUDxKjTcn3XaK7vGs3hUyXMW3+IBelZLNt2nK4tbNw3oC3DOzdDqzz9NXkmTZrEjh07cDgc3HXXXbWmJGoTpSgUiqukVZiZaSM78ujQeD7ZdIS31+zngQ83ERtmZsrg9tzYrXmtGCd8eXy3Gq9TUfN8+OGH/hbhqmm8pjUVijrGpNdye+9WLH9sIDPv6IHFqOPxhVsY+vJqPt2cXeML382DA2gerNZFFLWPUhQKRQ2j1QiGd47miyn9eOu3yRh1Gh5dsIXrX13DLyfdNXbye8mWoyzZcvTyBRWKq0QpCoWilhBCMKxTM5Y+1J/XJnSnxOnhXxvL+M0768nIurTBuKrwwU+H+OCnQ5cvqFBcJUpRKBS1jEYjGJXUnOWPDeQ3HQ3sOl7EjW/8yD1z0tl5rPDyFSgUfkYpCoWijjDoNAyN1bP6z4P449B41h84xchX1/D4wi3kFDr8LZ5CcVGUolAo6hiLUceUIe354c+DmdS/DZ9nHCH1xZW88f1eHC7P5StQKOoYpSgUCj9hM+uZNrIj3z46kGvbhfPPr3dz/atr+Plgnr9Fa/LceOONJCcn06lTJ2bNmgX4zIvHx8eTmprKvffey4MPPgjAyZMnufnmm+nZsyc9e/bkxx9/9KfotYI6R6FQ+Jm48EDevjOFlbtz+J9Pt3HrzHXccU0r/jK8A1bTxd22vnlHch1KWff8Y8M/2H5yO1ptzfk67xDagb/0+stly7333nuEhoZSWlpKz549uf7665kxYwabNm3CarUyePBgkpKSAHj44Yd59NFH6devH4cPH2bYsGHs3LmzxmSuDyhFoVDUE1ITIvnm0QG89E0m7689wKrMk7w8rhspcaEXLB8aaKhjCZsOr776Kp9++ikAWVlZzJ07l4EDBxIa6vu/uPXWW8nMzARg+fLl7Nhx1jFnYWFhozA9UhmlKBSKekSgUceToxK5vmszHlmQwbi31vHAoHY8NKQ9+vNOdy9M97ksujWl5YWqavD8pddf/PIHd+XKlSxfvpx169ZhNptJTU0lISHhoqMEr9fLunXrCAhovIcf1RqFQlEPSY4NZelD/bmpRwteW7GX22b9xPGCc3dGfbwxm483ZvtJwsZLQUEBISEhmM1mdu3axU8//URJSQmrVq3i9OnTuN1uFi1aVFE+LS2N118/a+Q6IyPDH2LXKkpRKBT1FKtJz4u3JvHqhO7sPFbIDa+tYe3eXzl3VNQww4cPx+1207VrV/72t79xzTXXEBMTw1//+ld69+7NddddR2JiYoW58FdffZX09HS6du1KYmIiM2fOvEwLDQ819aRQ1HNGJzUnMdrKfXM3cse765k6ogP39m/jb7EaLUajkWXLlv0qPSUlhUmTJuF2uxk7dixpaWkAhIeHs2DBgroWs05RikKhaAC0i7Ty+YP9+PPHW3h26S725RQjJVzE7YWiFnjqqadYvnw5DoeDtLQ0brzxRn+LVGcoRaFQNBAsRh2vT+jByxGZvLZiL0EmHfFRjWdnTX3nxRdf9LcIfkOtUSgUDQiNRvDHtAT+NS6JUqeHvBLnrxa5FYqaxq+KQggxXAixWwixVwgx9RLlegohPEKIW+pSPoWivnJTjxbMvac3OYVl3PzmWg7mFvtbJEUjxm+KQgihBd4ARgCJwAQhROJFyv0D+LpuJVQo6jd7ThRxxzWxlDjd3DJzHTuOKku0itrBnyOKXsBeKeV+KaUTmA+MuUC5KcAiIKcuhVMo6jtf/HKMzYdPs3ByH3Qawe3v/KTMlitqBX8uZscAWZXus4HelQsIIWKAscBgoOelKhNCTAImAURFRbFy5cpqCWW326v9bFNE9deVUZP9lZ9fCkD2jo081k3w3Ho34/69hqm9AoixNszlR5vNRlFRUcW9x+M5574uiI6O5tixY6xZs4ZXX32VhQsXVvnZkSNH8swzz9CjR49alPDiVLW/HA7HFb2H/lQUF9rYd76PyFeAv0gpPeIy+wCllLOAWQApKSkyNTW1WkKtXLmS6j7bFFH9dWXUZH+9uXsdAKmpfQDo2auY8W+t4+UtXuZP6kW7SEuNtFOX7Ny58xyTHf6ymWS1WjGbzeh0uitqX6vVEhgY6Dc7T1XtL5PJRPfu3atcrz9/dmQDlY3UtADOdwCcAswXQhwEbgH+LYRoOpuXFYoroHV4IB/eew0gueOd9RwpH3Eoqk9hYSFjx44lMTGRyZMn4/V6AbBYLPzxj3+kR48eDBkyhJMnT1Y8s3DhQnr16kV8fDxr1qzxl+g1ij9HFD8D7YUQrYEjwG3A7ZULSClbn7kWQswGvpBSflaXQioUDYl2kRbm/r4342auY+J7G/h4cl9s5oubKq/PHH/2WYq3bSevBs2MGzt2oNlf/1rl8hs2bGDHjh3ExsYyfPhwPvnkE2655RaKi4vp0aMHL730Ek8//TTTp0+vsPfkdrvZsGEDS5cuZfr06SxfvrzG5PcXfhtRSCndwIP4djPtBD6SUm4XQkwWQkz2l1wKRUNhwX19WHBfn1+ld4wO4q07kzl0qoR7/5OuvOZdBb169aJNmzZotVomTJjADz/8AIBGo2H8+PEA3HHHHRXpADfddBMAycnJHDx4sM5lrg38ejJbSrkUWHpe2gUtakkpJ9aFTApFY6Bv23BeGpfElP9u5tEFGbxxew80moZl76PZX//qd78O56+NXmyttHK60WgEfOsVbre79oSrQxrm1giFQsGs1fuYtXrfRfNHJTXnf0Z2ZNm247y6Yk8dStZ42LBhAwcOHMDr9bJgwQL69esH+HxQfPzxxwB8+OGHFemNFWXrSaFooHy303e0aNKAthctc0//1uw6XsQry/fQoVkQwzs3qyvxGgV9+vRh6tSpbN26lQEDBjB27FgAAgMD2b59O8nJydhsNmU9VqFQNFyEEPzv2M7sPWnnsY8yiAvvS4dmQf4Wq15jt9sBSE1NveRW5hkzZjBjxoxz0iqfTQgPD280axRq6kmhaOSY9Fpm/TYZi1HHpP9spNDh8rdIigaGUhQKRRMgKsjEm3ckczS/lKmLfkHK88+2Kq6EM6OOpoJSFApFA8Wk12LSV/2MQXJsCI8PS2Dp1uPMW3+4FiVTNDbUGoVC0UCZc3evK35mUv82rNt3iqe/2EFybAgdo9V6heLyqBGFQtGE0GgE/xqXRHCAngc+3KQO4ymqhFIUCkUD5dXv9vDqd1d+PiLMYuTFW5PYf7KY2WsP1rxgikaHUhQKRQPlx725/Lg3t1rPDoiPYEiHSN5YsZe8YmcNS9Z0GDlyJPn5+f4Wo9ZRikKhaKJMG9mBEpeH/1ue6W9RGixLly4lODjY32LUOkpRKBRNlHaRVib0asm89YfZd7Jpbfe8FMXFxVx//fUkJSXRuXNn5syZw7hx4yryV65cyahRowCIi4sjN7d6o7qGhNr1pFA0YR65Lp7PNh/l+WW7ePvOFH+Lcw5rPsrkxMECtDVoZjy8pYX+4+IvWearr76iefPmfPnllwAUFBTwt7/9jeLiYgIDA1mwYEGF5dimghpRKBQNlBCzgRCz4arqCLcYuT+1Ld/uOMH2owU1JFnDpkuXLixfvpy//OUvrFmzBpvNxvDhw1myZAlut5svv/ySMWPG+FvMOkWNKBSKBsrM3ybXSD2392rFy99m8tnmI3RqbquROmuC/uPi/WJmPD4+no0bN7J06VKmTZtGWloa48eP54033iA0NJSePXv61fS5P1AjCoWiiRMSaCA1IYLFW47i8SrTHkePHsVsNnPHHXfw+OOPs2nTJlJTU9m0aRNvv/12k5t2AjWiUCgaLP/4ahcAfxne4arrGtMthuU7c1i//xR924VfdX0Nma1bt/KnP/0JjUaDXq/nzTffRKvVcsMNNzB79mzmzJnjbxHrHKUoFIoGyqZDp2usrus6RhFo0PJZxpEmryiGDRvGsGHDfpX++uuvV/jFPkNjMSN+Ofw69SSEGC6E2C2E2CuEmHqB/DFCiF+EEBlCiHQhRON2I6VQ+IkAg5ZhnZuxbOtxZdZD8Sv8piiEEFrgDWAEkAhMEEIknlfsOyBJStkNuBt4p26lVCiaDjd2i6GozM33u3L8LYqinuHPEUUvYK+Ucr+U0gnMB87ZcyaltMuzhvMDAbXSplDUEn3bhhFuMfJZxhF/i6KoZ/hzjSIGyKp0nw30Pr+QEGIs8BwQCVx/scqEEJOASQBRUVHnuCS8Eux2e7WfbYqo/royarK/NA4HQI32f/cwD9/tOMGX335PoF7UWL1VxWazUVRUVHHv8XjOuVdcmqr2l8PhuKL3xp+K4kJv4a9GDFLKT4FPhRADgBnAdReqTEo5C5gFkJKSIi/l6/ZSrFy58pJ+chXnovrryqjJ/qqNbre1Oc23/15LWVh7ru/RouYbuAw7d+4854yCP85RNGSq2l8mk4nu3btXuV5/Tj1lAy0r3bcAjl6ssJRyNdBWCNG0t2QoFLVIUotgIq1Gvtl+wt+iKOoR/hxR/Ay0F0K0Bo4AtwG3Vy4ghGgH7JNSSiFED8AAnKotgbJ3n6Y4R3J8fwFanQatXoNWp0GnL7/Wa9DpNAhN3Q/JFYrzmb5kOwB/H9WpxurUaARDE6P4dPMRHC7PFblabSxYLJYm5xP7cvhNUUgp3UKIB4GvAS3wnpRyuxBicnn+TOBm4E4hhAsoBcbLWvQK/+XrW3C7JAdXbLxkOY1OoNNr0ek16AwadAYtOoMWvdF3rTdo0RsrBZMWvVGHIUCL4UwcoMNg0mE06zAE6NDq1CF5xZWx42hhrdQ7NDGKeesPs3ZfLoM7RNVKG4qGhV8P3EkplwJLz0ubWen6H8A/6kqe0Y90Z2P6Jrp06orH7fUFlxe367xrlwe303ftKvP40p0eXE4PxflluJ1eXA43rvK4KqpNp9dgNOswBuoxBeoxmnWYLHoCLHpMgQYCrHpMFj3mIAPmIAMBVoNSLopaoU/bMCxGHd9sP9HkFcULL7zA3Llz0Wg0jBgxgueff57U1FS6devGhg0bKCws5L333qNXryv3X96QUCezKxHd1oYlSxDbOazG6pRS4nF5cTo8OB1uXA4PzlI3TocbZ6mbslIPzlIXZSXuiuAodlFwspQTBwtx2F14PRfWNKZAPWabgcBgI4HBRizBRiwhRiwhJiyhRqyhJgwm9V+suDKMOi0DEyJYvjMHr1ei8dNU6/ezZ3Fs3x60upqb/oqMbcOgiZOqVHbZsmV89tlnrF+/HrPZTF5eXkVecXExa9euZfXq1dx9991s27atxmSsj6i/IrWMEKJiasocdOUmoaWUuBweSoqcOOwuSgqdlBQ6KS1yUlzgpKSgjOL8MvKOFlNSUPar0YvRrCMoPICgcBNB4QEER5qxRQYQHGXGHGRACLXeovg1aYlRfPnLMTZn5ZMcG+JvcfzC8uXL+d3vfofZbAYgNDS0Im/ChAkADBgwgMLCQvLz8xu1pzulKOo5QgjfekaAzneS5BJ4PV6KC5zYT5dhz3NQlOeg6JSDwlOlnDpSzIEtueeMTgwBOkKamQlpZia0uYWw5oGExVgw25QCaQi0iQistbpTEyLRaQTf7IRdZCkAACAASURBVDjuN0UxaOIkv26PlVJe9Htwfnpj/74oRdGI0Gg1WENNWENN0PbXfgW8Xok9z0F+Tgn5J0o4fdwXDm3PY9e64xXlAqx6wltaiWhpITI2iMi4ICwhxkb/ZWhoPHdT11qr2xagp0/bML7dcYJpIzrWWjv1mbS0NJ5++mluv/32iqmnM6OKBQsWMGjQIH744QdsNhs2W/3x41EbKEXRhNBoRPk0VACtEs9dhym1O8k7UkzuETu52XZys4rIWJ5VMQIJCDLQrHUQzdraiG4bTGSsVS2mN3KGJkbx5Ofb2Ztjp12kxd/i1DnDhw8nIyODlJQUDAYDI0eO5NlnnwUgJCSEvn37VixmN3aUolAAEGAxEJNgICbh7DSDx+UlN9vOiYOFnDhYwPH9hRzY4nMkr9NraNbWhlMvORFbSESs1W+Lnk2VaZ/8AtTeyGJIR5+iWLk7p0kpispnKKZOncrUqb8ybM3NN9/Mc889V5di+RWlKBQXRavXENU6iKjWQfgOzkNJoZPj+wo4suc0R3bnc+qI5OOt6RjNOlp0CCWuSxixXcIIsFydL2fF5dl/srhW648JDqB9pIVVmSe5p3+bWm1LUb9RikJxRZiDDLTpHkGb7hEALP/qe2LDEjm84xSHd+Sxb1MOQkCzNjZad4ugbY8IgsIC/Cy1orqkJkQwZ+0hSpxuzAb15wJq1ghjQ0H9zyuuCp1J0L5nFO17RiGl5OThIg7+ksuBX3JZu2gvaxftJTIuiPjyMtXZIqzwHwPjI3l7zQF+2n+qyR++a8ooRaGoMYQQvl1SsUH0GtWG/JwS9m3KYe/GHH5YuIcfF+2lZcdQOvRpRptuEWoxvAGQEhdCgF7Lqt0nlaJowlxSUQghjFLKsroSRtG4CI40kzw8juThceQdLWb3+uNkbjjON+9sJ8Cqp2PfaBL7xWCLUFNT1SGxeVCtt2HSa+nTNoxVmSdrvS1F/eVyI4p1QA8hxFwp5W/rQiBF4yS0eSB9xral95g2ZO3IY/uaI2z+NotN3xwmrks4SYNbEJMQos5qXAE1aTX2UgyMj2DFrhwO5hYTF157h/wU9ZfLKQqDEOIuoK8Q4qbzM6WUn9SOWIrGikbjs6UV2zkM++kytq85wrbVRzj4Sy5mmwFbeADWMN+hQUuoyWe/qtxuldGs97f4TZLUBN/GhVWZJ5u0opg4cSI33HADt9xyi79FqXMupygmA78BgoFR5+VJQCkKRbWxhBjpPboNySNi2fNzDkczT1N4ysGxvQXsyc9Bes81XGUwabGGldutigjAdsZ2VVQA1hBTk/MT8sj8zQC8clvVPZVVh9iwQOLCzKzKPMldfeNqtS1F/eRyiiJaSnm/EGJzuatRhaLG0em1dOwbTce+0RVpXq+kpMCJPd+BPa+MolM+21WFp0rJzynl8I48PC5vRXmtXkNwpJmQaDOh0YG07BhKVFxQo1YexwocddbWwPgIPkrPbjLOjGbMmMG8efNo2bIl4eHhJCcnn5O/ceNGHnvsMex2O+Hh4cyePZvo6GjefvttZs2ahdPppF27dsydOxez2czChQuZPn06Wq0Wm83G6tWrcTgc3H///aSnp6PT6fjXv/7FoEGDmD17NosXL6akpIR9+/YxduxYXnjhBT/1hI/LKYppwEJ8IwulKBR1hkYjyk2mG6H1r/OlV1Jc4KQgp+Qc21U5BwvZuzGHDUsOYA4yEJcUTrchLQlp1nSnTGqCgQkRzFl3iJ8P5tG/fUSdtJm/ZB+lWQWUamtuc6aheSDBo9peskx6ejqLFi1i8+bNuN1uevTocY6icLlcTJkyhc8//5yIiAgWLFjA//zP//Dee+9x0003ce+99wLwxBNP8O677zJlyhSefvppvv76a2JiYsjPzwfgjTfeAGDr1q3s2rWLtLQ0MjMzAcjIyGDz5s0YjUYSEhKYMmUKLVu2xF9c7n/glBDie6C1EGLx+ZlSytG1I5afOPgDwae3wkEdCA0gfHFFEOVBc17Qns3XaM9N02h915ryNI0WNLqz12rxtlqISoqkstkRAEexi0PbTnFgSy6Z64+za+0xkga3JOX6OOWfo5pc0yYMvVbw495TdaYo/MUPP/zAmDFjCAjw7cYbNercWffdu3ezbds2hg4dCoDH4yE62jca3rZtG0888QT5+fnY7XaGDRsGwLXXXsvEiRMZN24cN910U0U7U6ZMAaBDhw7ExsZWKIohQ4ZUGBpMTEzk0KFD9VpRXA/0AOYCL9W+OH7mg1vo5i6FLXXYptBUUhw60Op88flBq68U632xVg9aQ6W4POiMlWKjL64IAb5YHwA6E+jNoD8Tm33perOvTANVYqZAPQm9m5HQuxklhU7WfbaPzd8eZveG4/Qd25b4Xs0a9ZRUbWA26OjWMph1+3LrrM3gUW3R+sHM+OW8LUsp6dSpE+vWrftV3sSJE/nss89ISkpi9uzZFae4Z86cyfr16/nyyy/p1q0bGRkZl2zHaDRWXGu1Wtxud/U+TA1xSUUhpXQCPwkh+kopTwohAqWUNWZgRggxHPg/fD6z35FSPn9e/m+Av5Tf2oH7pZS192f8t5+QsWkj3ZK6gvQCEqT0XVfE3kp5XvB6KuV7KqV5zs0/k+Z1XyStPN3rBo/Ll+5xg9d1Nq0idvliV6nv2u0Ej9OX5ikDd9nZa281XzChBUNgpWABo/W8EESrIyfh530QEAwmG5hCfNcBIWAK9o2k/Ig5yMCQOzvSuX8Mqxdksnz2TrauOkK/ce1p1rphm4buUcd+Ivq0Def1FXsoKHVhC2i8O9D69evHfffdx7Rp03C73Xz55ZcV00kACQkJnDx5knXr1tGnTx9cLheZmZl06tSJoqIioqOjcblczJs3j5iYGAD27dtH79696d27N0uWLCErK4sBAwYwb948Bg8eTGZmJocPHyYhIYFNmzb566NflKqOw9sJIVYBFqCVECIJuE9K+YfqNiyE0AJvAEOBbOBnIcRiKeWOSsUOAAOllKeFECPwrZP0rm6blyW2L/kHnNBmYK01Ued4PT7F4XZUCmU+JeN2+OKKUOILzuKzsbMYnHYos/vi/CwoK4SyIigrpI3XDQfmXrhtofEpC3OYLwSGl4cIsET5Ymuz8hDtG8XUElGtg7jlz8ns3nCcdZ/uY9E/NtK8fTCtk8JpnRTRIA/9/WV4hzptr2/bMF79bg8bDuQxNLHxntLu2bMno0ePJikpidjYWFJSUs7xN2EwGPj444956KGHKCgowO1288gjj9CpUydmzJhB7969iY2NpUuXLhQVFQHwpz/9iT179iClZMiQISQlJdGhQwcmT55Mly5d0Ol0zJ49+5yRRH1CXG6YBSCEWA/cAiyWUnYvT9smpexc7YaF6AM8JaUcVn4/DUBKeUHbvUKIEGCblDLmcnWnpKTI9PT0asm1cuVKUlNTq/Vsk0NKVq/4mgE9k8BRAI58X1x6GkryyuNTZ0NxLhSf9F1zgffOHAa2FmBrCcGtIDgWQmIhJM4X9DXzx9zpcPPLiiz2bszh1BHfADmshYV2yZG0S44kONJcI+1ciIb8fpW5PXR96htu792q1g777dy5k44dzzpK8peHO7vdjsVioaSkhAEDBjBr1ix69OhR53JcKVXtr/P7GUAIsVFKmXKh8lVe2ZNSZp13atZT1WcvQgyQVek+m0uPFn4PLLtYphBiEjAJICoqqtoWHu12e5O0Dlld7KVuVm7aXSnFCDTzBQ2+Meh5rgyE14PeVYjBeRqDMx+DMw9jWR7GslyMZScxHd6CKfNbtN6z1mMkgjJjGKUBzSkxx1BibkGJuQXFgbE4DcFXvp5ihmb9IdQuKDwChYftrP/czvrP9xMQCrbWAlsr0Blrdi2jJt+v1zb7tsdO6W6qkfqqQjsbfLvlEAOttWPSw2azVfwKB99CceX7uuLuu+9m9+7dOBwObr/9dtq3b+8XOa6UqvaXw+G4ovewqooiSwjRF5BCCAPwELCzyq1cmAt9Ay84vBFCDMKnKPpdrLLycx6zwDeiqO6vtob8i88f1Fp/SekbgeQfgrwDiLz9mPL2YTq1l5DcH+Fo4dmyAaEQmQjNupwNkR19i/xXQFGeg70bc9i9/jjHN9rJyRC0Tgqn66AYotsF14h5kZrsrzd3+xZTU1P71Eh9VWG73Ms/v95Nl5Q+hFlqfppk586d5/wi9teIYuHChXXeZk1Q1f4ymUx07171g5pVVRST8S06xwBHgK+BB6rcyoXJBirv92oBHD2/kBCiK/AOMEJKeeoq21Q0FIQAS4QvtDhvNCwl2E/AyV2QswtydsCJ7bBpjm9tBXy7vZp1gebdfc+36AmhbS458rCGmug+tBXdh7YiN7uIXeuOs2vdMfZtyiEsJpCug1qS0LsZWn3TtXrbt63Phe5P+/O4vmv0ZUorGgtVUhRSylx8pjxqkp+B9kKI1viUz23A7ZULCCFa4TMT8lspZWYNt69oqAhxdhG8TerZdK8H8vbDsS1wdDMczYAt/4Wf3/blB4RCy17Q6hpoeQ1EJ4HhwusR4S2s9LvVSu8xbdjz8wm2rszm+w92sX7JfroNaUWnAc2b5JmMLjE2LEYda/flKkXRhKjSmy6EaAG8BlyLb3roB+BhKWV2dRuWUrqFEA/iG51ogfeklNuFEJPL82cCTwJhwL/Lh/3uiy22KBRotBDe3he6lBtu83p8I4/snyHrZ8j6CTK/8uUJLUR08I06Ynr4Rh2Rib6zLOXoDVoSr21Ox77RZO8+zaavDrH2k71s/PogPdJi6ZLaAr2x8Zu0OINOq6FX61DW7VOD+6ZEVX8SvQ98CNxafn9HedrQq2lcSrkUWHpe2sxK1/cA91xNG4omjkYLUZ18IXmiL604F7LWw5FNcCwDMpdBxge+PH0gtEiGVn0htq9PgRitCCFo2SGUlh1COXGgkA1fHGDdp/vI+C6LlBFxdBrQHK22bqekrm0XXqftnaFv2zBW7MrhWEEp0baGt61YceVUVVFESCnfr3Q/WwjxSG0IpFDUOoHh0OF6XwDfmkf+IchOh6wNcHgdrPoHFXsrrM0hIt63RTcwkqjACEYNjuBo9wjW/yhYsyCTrSuzufaWdsR2DqsznxoPDWlfJ+2cT5/ydYp1+05xU48WfpGhPvHkk08yYMAArrvuOn+LUmtUVVHkCiHuAP5bfj8BUGNPReNAiLNnNc5MWTkK4PB6OLEVcvfAyd2w60vfGRDps1rbHLhRwsHgFNaeupsv3yihVfAh+nfbQ3DzUN95kJBY35kQSzO/n1KvKTo2CyLIpOPng3lKUQBPP/20v0WodaqqKO4GXgdexvczay3wu9oSSqHwOyYbxKf5QmW8Ht9hQvsJsB9HFJ2gddExWhVsZtuug2w40J35q5qTEvgR3QOfQyvKTajoTBCZSLyMAOsh6HLrRRfSq8pd720AYM7dva6qnitFoxGkxIXy88HTddpuXXHw4EFGjBhBv379WLt2LTExMXz++efs3r2byZMnU1JSQtu2bXnvvfcICQk5x6HR1KlTWbx4MTqdjrS0NF588UVOnjzJ5MmTOXz4MACvvPIK1157rZ8/5ZVRVUUxA7hLSnkaQAgRCryIT4EoFE0Hjfbstl3OGibQAkmjoF1+GT8s3MP6jb8hM+AuhqQ5iTIdhFP74fgvRGT/AEu+htUvwoh/QIeR1RbF4braM6/VJyUuhBW7csgrdhIaaKiVNpYtW8aRI0fQamtus0CzZs0YMWLEZcvt2bOH//73v7z99tuMGzeORYsW8cILL/Daa68xcOBAnnzySaZPn84rr7xS8UxeXh6ffvopu3btQghRYU784Ycf5tFHH6Vfv34cPnyYYcOGsXPn1R5Dq1uqOhbuekZJAEgp84DadaulUDRAAoONDLu3M9c/0BWXW8ui/2pZf2wwnqHPwMQv+PHaeTDxSzBaYP4E+PA2KDjib7GvmJ5xoQBsPNQ4RxWtW7emW7duACQnJ7Nv3z7y8/MZONBnB+6uu+5i9erV5zwTFBSEyWTinnvu4ZNPPsFs9o0Yly9fzoMPPki3bt0YPXo0hYWFDeKUd2WqOqLQCCFCzhtRNL1N5ApFFYnrEk7032ys+WgP6UsPcmjbKYbenehbD4nrB/ethvUz4fvn4M0+cP2/zq6PNAC6xNgwaDWkH6w9A4EjRozw28ns8818nxkdXAqdTseGDRv47rvvmD9/Pq+//jorVqzA6/Wybt26Cv8WDZGqjiheAtYKIWYIIZ7Gt0bhX998CkU9x2jWc93EREbc14WiUw4+evZnTu+XPj8EWj30nQKT10B4PCz6PSy6x7eI3gAw6bV0bWHj54N5/halTrDZbISEhLBmzRoA5s6dWzG6OIPdbqegoICRI0fyyiuvkJGRAUBaWhqvv/56Rbkz6Q2Jqp7M/o8QIh0YjM9G003nmQNXKBQXoU33CCLjglj+/naObMjnW7mDAbfFYwrUQ1hb+N1XsOYl35bcY1tgwnxf+mUY0jGyDqS/OClxobz7w/4m40d7zpw5FYvZbdq04f333z8nv6ioiDFjxuBwOJBS8vLLLwPw6quv8sADD9C1a1fcbjcDBgxg5syZF2qi3lIlM+MNDWVmvO5Q/VV1vF7JojdXcnI7BFj0DJyQQJvuldyKHlgDH93p2347bs655knqId/tPMHv56SzYNI19G4TViN11hcz4w2V2jIz3jg2disUDQCNRhDRSXDr1BTMNgPL3trKV7O2UVrk9BVo3R/uXeFz4jT3Jtg426/yXo7kcg976Y10QVtxFqUoKuE+dQpht+MpLMRbXIy3rAzpdl/Wh65CcSVEtLJyy9QUrrmxDQd+Ocn8ZzaQtbN8rj+0Nfz+G2g7GJY8DKte8J0cvwDj31rH+Ld+7be5rgg2G4iPsjSZdYqmjNq5VIm91w0lsrSUC5qp1WoRWi1CpwOdDlE56PUIgx50et915WAwVORrjEaE3oAwGhFGg+/eYDx7bzIhjEY0AQG+6wAzmgATmoAAREAAmsBAX111ZCJCUXtotRqSh8cR2zmMb97ZzuJXM+g+tBW9x7RBawqCCf+FxVPg+//1He4b8YLvDEc9IyUulCVbjuLxSrQa9V42VpSiqETU1Klk7thOu9atweNBerzgcSNdbqTXA2430u3xjTLcLqTbDS6XL9/tRrpcZ2OnE29xMdLp9IXytDPB63SC233lQup0aMxmNIGBaALNaAMtvmuLBY3VgtZiRRNkRWsN8sVBNrTBNrQ2G9rgYLQ2m0/ZKeoF4S2s3PrXnvywcA+bvznMiQOFDLu3M+YgA9z4ps+v+NpXoTQfxr51jmXb+kDPuBA+XH+Y3ceLSGwe5G9xFLVE/Xrr/EzI+HGUrlxJWB0tzkq3G1lWhtfpRDoceB0OX1zqQJY58JaW4i0pRTp8sbekxJdWXOwLJSV47Xa8djuunBN4i+x4i4rwlpRcsl1NUBDakGB0oWFoQ0PRhYaiiwhHGxaGLjwCfVQkushIdOHhCEPtnLpVnEVv0DLoNx2IaR/Mirm7WPj8z4y8vysRLa2QNgMCguG7pwEJY2fVK2WREus7eJd+KE8pikZM/XnjmiBnpq40gYE1Wq90u/EUFeEtKsJTUICnoBBPQT6e0/l48stDXh7u03m4srIozcjAk5f367lwIdCGh6GPaoY+uhm66GgMMTHomjfH0LIl+hYtLyyAolrE92pGcJSZpW9u5ZMXNjL07k6+XVH9/whCA8uf8hWsR8qiRUgAkVYjmw6d5s4+cf4WR1FL1I+3TVGjCJ0OXUgIhIRU+RnpduM5fRr3yZO4T57ElZOD+0QO7hPHcR07TtmBA9h/XIs8b7QSYbVwsF17DHFxvtCmNca2bTG0bInQX5nPagVExgZx67QUls3cyrJZWxkwPp4uqS2g36M+ZfHtk7547CxuqAce5oQQdGsZzJbshnFQsLZ46qmnsFgsPP7441dd1+zZs0lPTz/nkF5Vyc/PZ+7cufzhD3+4ajkqoxSFAihXLhER6CIiLlpGSoknPx/XkaO4srNxZh3m0PoNWJxOin/8kYJPPz1bWK/HGBeHMSEBY3w8po4dMHXsiC7cP852GhKBNiNjHu3ON29vY/X8TIrzy+g9pg3i2ofB6/ZNQ5nD+O3w5y/pA7yuSGoZzDc7TpBf4iTY3LimKqX0naTXNBAT8QUFBfz73/9WikLhP4QQ6EJC0IWEENC5EwDb2rcnuXxNx1tcTNn+Azj376Ns7z7K9uyhZNNGCr/4oqIOXWQkps6dCejalYCuXTB17YrWYvHHx6nX6A1aRkzuwqoPd7Pxq0OU2l2k3p6A6PcYFJ+Cn96g1BQB1z5KgMG/u6G6twwGICMrn9QE/54WrwnOmBkfNGgQ69at48Ybb+SLL76grKyMsWPHMn36dAD+93//l//85z+0bNmSiIgIkpOTAUhNTeXFF18kJSWF3NxcUlJSOHjwILNnz2bx4sWUlJSwb98+xo4dywsv+Cwhvf/++zz33HNER0cTHx9fYWtqyZIlPPPMMzidTsLCwpg3bx5RUVE89dRTHD58mP3793P48GEeeeQRHnroIf7+97+zb98+unXrxtChQ/nnP/9ZI33iV0UhhBgO/B8+K83vSCmfPy+/Az6Xqz2A/5FSvlj3UiqqiiYwkIAunQno0vmcdE9hIY5du3Ds2OELW7dhX7Gi/CENxoQEzD16ENCjO+bu3dE3b+4H6esfGq2G1Ds6YLIa2PTVIbweL4N+2xFN2jNQcoqJ30jYtJgFfxzrVzm7tLAhBGzJKqhRRZGZOYP8gq1oa3A9xmrpSHz83y5bbvfu3bz//vvceOONfPzxx2zYsAEpJaNHj2b16tUEBgYyf/58Nm/ejNvtpkePHhWK4lJkZGSwefNmjEYjCQkJTJkyBZ1Ox9///nc2btyIzWZj0KBBdO/uM87dr18/fvrpJ4QQvPPOO7zwwgu89NJLAOzatYvvv/+eoqIiEhISuP/++5k+fTq7d++ucXtSflMUQggt8AY+v9vZwM9CiMXn2ZDKAx4CbvSDiIoaQhsURGCvXgT2Outgx1NQQOnWbZRu3kzJpo3kf/opp+fNA0DXrBnmHt0x9+qFuVcvDK1bN9mzI0IIrhnTBq1Ow89fHEB6YfBdHdGMeR22/wdO7YN930PbQX6T0WrS0z7SQkZW4zmhHRsbyzXXXMPjjz/ON998U/GH2263s2fPHoqKihg7dmyFKfHRo0dXqd4hQ4Zgs9kASExM5NChQ+Tm5pKamkpE+bTv+PHjycz0nebKzs5m/PjxHDt2DKfTSevWrSvquv766zEajRiNRiIjIzlx4kSNff7z8eeIohewV0q5H0AIMR8YA1QoCillDpAjhLi+LgRy7M4jIBcc+/IRWoHQakCn8V3rNAidAK2m/NqXrqgeWpsNS79rsfTzefqSbjeO3bsp3bSZ0s2bKEnfSOHSZb6yEeFY+vYl8NprCezbt8mtcwgh6HVDazQawfrF+wEYcldHCE+A47/AwrvgnhUQ3s5vMia1CGb5zhNIKWtMqcfH/81vtp4Cy3ciSimZNm0a99133zn5r7zyykU/p06nw+v1uct1OBzn5J1vvtxdfpbqYnVNmTKFxx57jNGjR7Ny5Uqeeuqpy9ZVG/hTUcQAWZXus4He1a1MCDEJmAQQFRXFypUrr7iOvSu24pWSrVt+AQQCgabi6uw/DQIhff+5AuGLKweNAM3Z6zMBjcZ3rT2jcMrL6QToNKA/GwudBo3m3HCm/vqE3W6vVl9flJYtfGHUKLQnT6LfswfDrl24vltBweeLAXDFxlLWtQtlXbvibtGiXizoVpWr6i8zRHYR7F5/nJOnj5PvKkNjaoHTI3G/O4pNPf6JW++f9R5zqYvTJS4WLvueSHP1F35tNts5Tn08Hk+dO/mx2+14vV6Kioro378/zzzzDKNHj8ZisXD06FH0ej3Jycncf//9PPDAA7jdbj7//HPuvvtuioqKiImJ4ccff6Rjx4588MEHSCkpKirC4XDgdDorPo/b7aakpIROnTrx0EMPcfDgQYKCgpg/fz6dO3emqKiI06dPExwcTFFREe+8805Ff5SVlaHX6yvq8nq92O12zGZzlRwjORyOK3oP/akoLvTtrrZRJSnlLGAW+KzHVsei6cHM/eTn52M2m/F6vUivxCu9SK+3PPbde70SKb145dnYW3Ev8Xq8UEteKnVCi06r8wWdDr1ej96gx2A0YDAZ0Rv06PV6DAbDr8KZYWrlYDKZMBgM1VZAdWU9Vnq9OHbupHjNGuzfr6T0iy+xLPkCfcuWWNOGEjRsGKYuXeqdIj2fq+0vOVDy46K9bFmehTbWhDU0GMPwjzDMGUW/Y+/AHYv8Yuoj4mgBc3b8gLF5Aq0TwlmRV8TKvEKusVm4v1XV1y127tx5zgjCHyMKi8WCRqPBarVy4403cujQIdLS0iryPvjgA/r378+ECRPo378/sbGxDBw4EKPRiNVqZdq0aYwbN46FCxcyePBghBBYrdaK79qZz6PT6TCbzbRv357p06eTlpZGdHQ0PXv2xOPxYLVaefrpp5k4cSIxMTFcc801ZGdnY7VaK76/Z+rSaDRYLBbCwsLo168fffr0YcSIERddzDaZTBXTaVXBb2bGhRB9gKeklMPK76cBSCmfu0DZpwB7VRez/W1mXErpUzTlscfjqYg9bg/uMhduhy94HC5cDicehwt3mRt3mROXw4XH6cJV5sLt9AWX04XL6cbjduF2e3DjwSO8uPHixvP/7d15fJTVvfjxz5l9y8xkD1kgAUICYRMQRFDBfUHRCqJevWJV3L3aFbtYr71tub1eb3+i1VLrdWvVgkvVeq1FjCggsglhT4CwZd8zmSWznN8fM8YASUhCksly3q/X83pmnufMPN8chvnOeZ5zzkNABMOPNUECIkRABPHLICEZOm28QghMJhMmkwmzOwvrOwAAIABJREFU2XzCYrFYWhar1dqyWCwWNBpN1KYZD1RX4/r0Uxr+8TFNGzZAIIB+xHAcV1+D4+p5GEaM6POYOqOnPl9rXt3Lqs1HGTM9hUdumRieafb9f4MLlsLcR3sk1s6q9Qf4vKaRBz7ejSXVSl2kQeHUaakLBHl+3AiuTe7cmB41zfiZ6a1pxqPZotgEZAshsoDjwI3AzVGMp8cIIXr0hvAnk8EQIXeAkNtP0OUn1OQn5PITdDUTavQTbGwOLw3N+F1e/DKInwB+EVkbJUGLIGABv0niN4Twa0P4RABfsBmP10NtbS1ut/uUc6yt/0abzYaUkpKSEmJiYoiJicFut+NwOFrWhl6aAkQXH49zwQKcCxYQrK+ncfVq6t97n6pnn6XqmWewzJhB7I2LiLnookE3DYkQgrn/koOvyU/xuiqKJ1WROeU2OLwhfPOj4TPCs8/2Aiklx3x+ttQ3sbG+iS/rXOxt8iIBTaoF3AGemDyci+LtZJgMLPz6AI/sPcJIi5GJMZZeiUnpfVFLFFLKgBDiAeAfhLvHviil3CWEuCey/3khRAqwGbADISHEw8A4KWVDtOLuD4RWgzbGgDbGgP40tyuWQUnQ1Uyw3kewLrwEar0EayPrci/S36rVoRXoEszokyzocs1ok8wEHAKfSeL2umlqaqKpqQmXy4XL5eLw4cM0NDRw7Ngx3G3MMWW1WnE6nTidTuLi4oiNjSU+Pp74+HisVmuPnCrSOhw4r78e5/XX4y8tpf5v71H3179y/JHvoU1IIPbmm4i7+Wa0TucZH6u/0Gg1TF2UTVmVm4//tIvrfziV+HlPhe+Q99Zd4Vus2s+8m3Fls5/tjR52NLrZ3uhma4ObyubwRVOLVsPZditXJzmZ5bTxj3WHeW3jERZfMQWDLtyseGF8Jpdv3s/igkP8Y9oYEg1qtP5AFNVxFFLKD4EPT9r2fKvHZUB6X8c1mAitQOcwonMYYfip+6WUhBr9BKo8BKo9+Cs9BCrd+Eub8Oys+vaqkVZgSbHiGGbFkBqHPtuGYZiNz9atbTmVEggEaGxspL6+vmWpq6ujtraWkpISdu/efcK9PUwmEwkJCSQmJpKYmEh8fDxxcXE4nU703Zz+Qz9sGAn33E38XXfStG4dNa+9RtXTy6l+4U/E3nADcbffjj554A8KA3ho5dcE40LMc2n5+7M7WLB0GpYbXoYVc2HVHbD4g05dr5BSUtbs54DbxyGPjyK3j70uL7ubPC1JQQCjLEYuiI1hqsPKFLuFPKsZXaupxasy4njx82L2ljUwMT2clBMNel6akMU1Wwt5cPcR3ph8+lu8Kv2PGpk9xAkh0NoNaO0GjCMdJ+yT/hD+Sjf+sqbwUtKEd3c17s2R/tpaQbpVQ13DAQwjYjCMCN+APradOaaCwSD19fVUV1dTXV1NVVUVVVVV7N+/n23btp1Q1m63t7Q84uPjSUxMJCEhAYfD0alWiNBqsZ1/Prbzz8e7bx/VL/yJmldfpfb114m96Sbi77oTXXzP3L4zmrQ6DVfel8c7T27l/54v4NrvnYX2qv+Gd+8JT08++5E2X1fibeaDyjq+qm9iU30T5c3fdq00aQRjrCYuirMzzmZiQoyFCTYzNl3HSWdSRvjz8/XRupZEATAhxsLSkcP4RVEJX9Q2MjtWXXMYaFSiUNol9BoMqTYMqd92uZRSEqz34T/movloI+6CozRtKsO1vgQAbZwJ40gHxpEOTKOdaO0n9vWOi4sjLi6O7OzsE47ldruprq6mpqaG2tralsc7d+484TqJwWAgKSmJpKQkkpOTGTZsGMnJySf0KT+ZKSeHtP/6LYkPPUjVs7+n5pVXqH3zTeIW30b8HXeitfXs7L19LWmEnQtvG8vHL+xi/dtFnLfwRtj3Iaz5FYy6CIZNBML/dpvqm3jheBV/r6wjKCHdpGdWbAxT7BbGWExkWYykGfVounFKMM1pJsFm5OujdfzrzBP33ZaawB+OVvKbg6V8MMXW73unKSdSiULpEiEEOqcJndOEeXwC28yHueC88/CXuvAVN+A7VI9n17etDn2KBeOYWEw5cRgz7eFBjG34pmdVRsaJU5dLKXG73VRWVlJVVUVFRQUVFRXs2bOHrVu3tpRLSEggNTWVtLQ00tPTSUlJOaVDgSEjg9RlvyH+7iVULV9O9XPPU7dyFYkPPYjz+usRvdgBobdlT0um7GA9O9YcI2Wkg+x5v4OjG+Gdu6lZvJq3qt38pbSaPU1eHDotS9ITuS0tgUxz+wm2q4QQTEizs7vk1EuIJq2G72Wm8IN9R/lndQOXJjjaeAelv1KJQjljQiswpMdgSI8hZnYaMiTxlzbhK6rFu78W17oSXGuPI0w6TLmxmMfFY8qNQ9OJyeyEEC1dcjMzM1u2fzOIqbS0tGU5ePAgO3bsAMItj+HDh5OZmUl2djZJSUktv2KNWVmkPfUUcbfdRvmy/6TssV9Q+/obDHv8F5gnTeqVOuoL535nNBXFDXz66l7sP57C1kv/yFt7t/Hxl3tpFlomx1h4MieD65KdWHspKealOlhbWIXXH8SkP/EYi1LieOZIOcsOlnJxvL1brZZo6MkpxLsiPz+fJ598kg9aTarZFb/+9a/5yU9+0iOxqESh9DihERjSbBjSbMRckEHIF8RXVItndw3evTV4vq5E6DWYcuMwT0zEPDYOoevaaF4hBHa7HbvdTk5ODhBOHg0NDRw9epTDhw9z6NAhVq9ezerVq3E6neTk5JCbm8uIESPQaDSYJ01ixF/+TONHH1H+m2UU33gTzkU3kPTII2gd/f8X7y3nnDhWJKAB46JM/vbZQZbtLMKtdxAfP4Nbj7/LzedeS96oMb0e0/g0O8GQZG9ZI5MzTuxlptcIfpQ1jPt2H+a9irpOj61QukclCmVA0Ri1mPMSMOclIEMyfHqqoArPzio8BVUIkw7LxATMefHoU21oY7o37kEIgcPhwOFwMH58eAbbhoYGCgsL2bdvH1u2bGHjxo1YLBZyc3PJzc0lKysL+xVXYD3vPKqWL6fm1ddwfbKGlF8+QUwUBhJ2xdWTUqn1B3i7vJZ/VtWzurqBxmAI63AjIw95uSbWzr2XZqF/7rvw8RpYkg/a3u2empcaTrA7j9efkigArk1ysvxwOU8Wl3F1khNtP21VtDWF+IEDB7j//vuprKzEYrHwxz/+kdzcXBYvXozJZGLXrl2Ul5fz1FNPMW/ePILBIEuXLiU/Px+fz8f999/P3Xff3TJnU0JCAjt37mTq1Km89tprCCH46KOPePjhh0lISGDKlCkt8TQ1NfHggw9SUFBAIBDg8ccfZ/78+adMXX7VVVfxu9/9jqVLl+LxeJg8eTJ5eXn8OTLhZnepRKH0KaERmEY5MY1y4rxmFL6iOtzbKnBvq6DpqzIANDF69MNsGFKt6IfZ0Kda0cWbw/NldZHdbmfq1KlMnToVn89HUVERe/bsYefOnWzduhWdTsfIkSPDrY2HHsJ+zTWULn2UY/fci+M73yH50aVo+9HI4EBIsq3RzWc1jaw+XsOOBjchs454vY6rEp1clejg/LgYNpYVUfDRMcpyksm48rfw5i2w4VmY/XCvxpcea8Zh1rOrjesUABoheDgzmbt3HebDynquTmp/bMvPC4+xvc6F9jS9rbpivM3ML7M77nG/ZcuWNqcQX7JkCc8//zzZ2dls3LiR++67jzWR6fKLi4v57LPPOHDgAHPnzqWoqIhXXnkFh8PBpk2b8Pl8zJo1q2UqkG3btrFr1y5SU1OZNWsW69atY9q0adx1112sWbOG0aNHs2jRopaYfvWrX3HhhRfy4osvUldXx/Tp07n44ouBE6cuHzNmDN///vdZtmwZzzzzTI9NN64ShRI1QiMwjYnFNCaW0LWjaT7WiL+kCX+pC39JE41FdRAKj7sQRi2GNBv6jBiMw2MwDLd3ueVhNBrJy8sjLy+PQCBAcXEx+/fvb1k++OADRowYQd5Pf0LyunXU/+lFmj7/nKQf/RD7vHlR6anjCYbY0ehmY30TG+pcfFXfRFMwhACcW2oYodfy7OJpTLZbTjjnP/M7ozi2t4ZPXtrNjY9djil3HuQvg3HXQNzIXotXCMH4NDu7Stq/Neq8RCcjzWU8fbiceYmd6+7clz7//PNTphD3er2sX7+ehQsXtpTz+Xwtj2+44QY0Gg3Z2dmMHDmSvXv38vHHH7Njxw5WrVoFhO8+V1hYiMFgYPr06aSnhxPW5MmTKS4uxmazkZWV1dIj8JZbbmHFihUAfPzxx7z33ns8+WR4FiOv18uRI0eAE6cuz8nJ4fDhw6d0CjlTKlEo/YLGqG1paXxDBkL4K9z4j7toPu6i+Vgjri+O4wqGk4c23oQxM9wV1zjKgc5p6vTxdDodo0ePZvTo0VxxxRWUlZWxe/dudu/ezd8/+giNRsPIf3uItI1f4V36KDFvvknKz36GKTe3x//2b4Sk5IDbx/ZGN183utlc72aXy4M/Mkgx12piYUocs5w2ZsfauGf/JgCmOE7t3qs3aLnku3msWraZz/6yj8tu+C08OwM++B7c+k6vzribl+rgpXXF+IMh9G30ctMKwQMjkvje3qN8WtPIhfH2Nt/nl9npUZvr6eTkFQqFcDqd7f5CP7m8EAIpJcuXL+eyyy47YV9+fn6XpxuXUvLWW2+1XI/7xsaNG/tkunGVKJR+S+i+HcdhPTu8TfpDNB9vpPlII77iBjy7q3FvCXfF1SWYMWY7MWXHYsqN6/SpKiEEw4YNY9iwYVx44YWUlZWxY8cOCgoKKBqegSFzBGlHjzL8nnsZPf1skh96CEP6mU0Y0BQIsq/Jy54mLztdHnZFlqZgeDoVs0YwKcbC3RmJTLNbOdthJd7Qtf+uicNjOPvqLDb+7SDZ05IZedHP4f9+BHs/gLFXn1H8HclLtdMcDFFY7mJcattJYEFyLE8eCrcq2ksU0XL++eezePFili5dSiAQ4P333+fuu+8mKyuLlStXsnDhQqSU7Nixg0mRXnIrV67ktttu49ChQxw8eJCcnBwuu+wynnvuOS688EL0ej379+8nLS2t3ePm5uZy6NAhDhw4wKhRo3j99ddb9l122WUsX76c5cuXI4Rg27Ztp539Va/X4/f7uz3LQWsqUSgDitBrwq2ITAcx5xPuilvuxnegLny9Y0s5TRtKMWQ5iFs0pkutDDgxaVxyySUUFxdTUFDAbp2OQ2lpbPR6yfjxUibm5jL2gfvRtTMKvbXmUIhdLi9bGpr4UFr48YbdHPY2t+y3aTXk2czckBLHpBgzk2IsZFtMJ0yP0V1nXTqcoi0VfPbGPtJ+fhvGzf8LH/8Msi8FXc+NoWhtfFrkgnZJfbuJwqDRcN/wJH5WeJyNdS5mOPvPfdOnTJnCokWLmDx5MiNGjOC8884D4M9//jP33nsv//Ef/4Hf7+fGG29sSRQ5OTlccMEFlJeX8/zzz2MymbjzzjspLi5mypQpSClJTEzk3Xffbfe4JpOJFStWcNVVV5GQkMDs2bPZuXMnAD//+c95+OGHmThxIlJKMjMzT9ttdsmSJUycOJEpU6ac8cXsqE0z3puiPc34UNLf6ksGQri/rqDu/YMgIPba0Vgmn/ncTn6/n8LCQrZv2kzhwYOEBDgbGpg8ejTT//VfsdjCX3SBkKTQ7aXAFZ5Ib1uDm50uD77ItZY4QsxMjGW8zcxYq5lcm4nhJkO3xhQs+sMGAN68e2aH5SoON7Bq2WbGnZfGnBnH4dXr4OLH253e40yFQpIJj/+DBVPT+ff549st5w6GmLZhF2fFWPnzpPB1k4E4zfjixYuZN28eCxYsiHYog3KacUXpcUKnwTotBeNIJzVv7qPmjX2EvEFs5ww7o/fV6/WMGzeOcePG4fF42PrJJ3y6ZRuv+mH5X94mNCyN+oQUDjYH8EaSglmjYVKMme+mJXCW3cpUu4XCL9czZ/yU0xytc+46r3MXpZNG2Jl4UQbbVx9lzNlTSM25EtY+CZNuhpjTTD/cDRqNYFyqnZ3t9Hz6hkWr4c70RH57qIzCJi/Z1q61/pS+oxKFMijp4kwkLplI1Us7qf/wEKbc2E6dhpJS4glJ3MEQnlCIWn+AiuYAFT4/R73NHPT4OOT2ccA+nKbzv71OYfZ5iSs9zjQhmTtyOBeNyiTbajplnEBhD/6NF4/r/Jf8jKtHcnBbJZ++tpcb7/sl2j+cA2uegPnP9mBE38pLdfDXzUcJhiTaDk6h3Zoaz++Ky3nhWCX/mdOzPXX6yksvvRTtEHqdShTKoCW0gtjrsin/3Rbq3ikifnEeQgiqmwP8eP9Ritw+fKEQ3pDEG0kM37QG2qIBMkwGRlqMTHdaybaYGGM1MSoUoOnp37Ht6+0Ujc2lcruR/JQUArNnM27cODSa7t9DuiMHKl0AjEo8/fl9vVHLBTfl8MEz29m+Xc+Uc+6B9c/AjHshpf3TQ92Vl2rH3RzkUFUTo5Pajy/RoOe65Fj+WlbLoyPPrNWn9B6VKJRBTRdnwn5ZJvXvH8T9dSUVOQ5u3nGAEp+fuXExmDUaTFpNeK3RYNIKLBoNZq0Gi1aDQ6cl2aAn0aAj2ajH2N6X/uOPk7xlC0d/8lOKtFoKZ0xn1apVJCUlMXfuXHJzc3t8vMBP3i4ATn+N4hsjxseTOSGezR8Wk/PoQ1i3vAL5v4Ebz+xCZ1u+uaC9q6S+w0QBcFd6Am+W1fCX0hrmEm7V9bexFYNJd65Lq0ShDHq2mal4tlfyRf5BHnFZQcCqyaM5u43xB2fCMnUq2W+/hX3ZMjLfeJOyc2awy27nzTffJDU1lUsuuaRHj9cdsxZk8/oTG/nyoyouOvcB+PRXcHwLpE3t0eOMTrJh0GkoOFbP/MntdwkFGB9j4RyHlRePV3J5ioXq6mri4+NVsugFUkqqq6sxmbp2PUglCmXQExrBzsvTWHLwGImeIG/MymV0L1041VitDPvlL7HNmYP+Jz8ldecu6v7tIb6srOTll18mPj6evLw8EhMTe+X4p+NMtjDpogy2fXyE8efcQrL5ufB9K259u0ePo9dqmJ4Zx6qtx7jr/JEk2zuu77syErljZzH7zDHkNjZSWVkJhEcgd/VLbSjrTH2ZTKaWUeGdFdVEIYS4HPh/hO+Z/YKUctlJ+0Vk/5WAG1gspdx6yhspSgdWVzdwx9ESMrU6nl5bT4qtBs498/tJdyTmooswvZ3LsYcfwfkfv+KG22/n4Ny55K9dy+9//3umTZvGnDlzsFr7/qZJ067MZN+XZax9p4Trzv0eujU/g8PrYcS5PXqcf5+fx1VPf84PVm7nle9O77CFcHmCg3STnt8fr+Xts0ajj1wAz8/PP+3AMuVbvVVfvXOVrROEEFrgWeAKYBxwkxBi3EnFrgCyI8sS4Lk+DVIZ8D6qrOf2gkPkWEy8PSuXtJGx1H14kObSpl4/tj4tjRF/fo3Ym2+i4X//l4znnmfmmDFMnTqVzZs38/TTT/PFF1+cMGdQXzCYdMxemE1FcQPvfDGdRmMufPJLCPp79DijEm389MqxfF5YxSsbDndYVisEj4xIYVNDEwu/LqKyuWdjUc5M1AbcCSFmAo9LKS+LPH8UQEr5m1Zl/gDkSylfjzzfB8yRUpZ29N7dHXBX6mtm/foNzJzZuYuDCmzY0DP1FZIh3HV1SBkgGPIgQ93/8nSF4DNXiH82Sgq8kokmwR8ytNi1gpAnSNPKIwijBssVqaDt3HlwfzCAy9/xuIAOX7+tAO8bbyORmG9aiDsxnoJDxZTV1KDVaEhPTCQzJQlzF06zbK8O/86bFB/qVkz1ZX4ObfeiIcgI7YdYDB5IyoHEMaDtmdM9EvhsXxWVjR7m5CZjPc3Nqr4wxbPCMRJ7yM8DdQcQFSUkDIJ7m/eVmppqltzbvRssdTTgLpqJYgFwuZTyzsjzW4EZUsoHWpX5AFgmpfwi8vwT4MdSyg6zQHcTRdZnO/CEuvefTul/smQRZ/Mll/J/mPGe/gVKv3CILP6HH1MtonMdZyCzyzr2XzinW6/tryOz2/opd3LW6kyZcEEhlhA+PUVycjL5+fldDug2acDt82IyqotnneU9TX1JJNAI1AO1IGuR1AANfPNPGfRp8LtNaPQJIPUgdUjZ/bOiWhliTHMJCcFGAELM4uQTTSa/EX2o/cnSJCFqArU0BMOtCIvGjEVj6XZMLe8bDGH0+hEn/UALIXFroSs/U477wzPtpunrzjAqAVIXWUu6FkXndfY3qR34gXYlOx0jCUbv7PiApA8G6cZX32lFM1EcA1oPxUwHSrpRBgAp5QpgBYRbFN2Zf2gO38xd1LMX9Qaz/Px8LrhgJs3NVXg8h/F4DuN2F+P2FON2H8TtPkSo1WkkkykdmzUHT42Z/WsLqT8iGX/Bd5h9w63o+0nvFk/Aww8/+yGfHdvKwjELuT3vdjLsPTNquCfnxursXE8DWX+bS6y/6636imai2ARkCyGygOPAjcDNJ5V5D3hACPEGMAOoP931CaXnBYNufL5yvL4yfN4yfL4yvN7jeH0lBEOF5H9WQyj07akdIbSYTGlYLCOJi52F1To6smRTeaiET19aQWnRPoZlj2XhT+4hZVR2FP+6E9V563hgzQPsqNzBz2b8jEW5i07/IkUZ5KKWKKSUASHEA8A/CHePfVFKuUsIcU9k//PAh4S7xhYR7h57e7TiHUxCoWb8/nr8/lr8/jr8/lqa/dX4m2to9lfT3FwVWSrx+SoJBl2nvIdeH4vJlAokk552BSZzBmZTOhZLJiZTGhrNiXefqyk5xod/Wk7hxvVYnbFcft8jjDtvLqKXprfojuL6Yh5c8yAlrhKemvMUF4+4ONohKUq/ENVxFFLKDwkng9bbnm/1WAL393Vc/ZGUIUIhH8Ggu9XiIRhsIhh0Ewg2EQy4IutGAgEXgWBkHWggEGjA768nEKgnGHS3exydLgaDIQGDIRGbbSxxcedjNCZjNCRiNKZgMg3DaExGqw2fr8/Pzyc7e06771dXXsbGd/7Krs9WozMYmbngZqbNuxaD+czP9/ek9cfX84O1P0AndKy4dAVTk3t2pLKiDGRqZHYrvuYqpKzF6y1BygBSBk9avt0Wkv7w45C/ZV9I+pGhQKvH/m/XoebI4+bI48g65Dtp7Y0kBC+hkLdl3frUzukIoUWrjUGns6HTxaDT2TGbMoixjUOnd6LXOdDrnZElFr0+FoMhHr0+9pSWQHdVHz/KV++uZM8X+Wi0WiZfdhXnXLcIi8N5+hf3ISklr+x+hae2PMVo52ievvBp0mwdTzmhKEONShStrF9/ASHpZd363nl/IXRoNAaEMKDR6NFoDGg0xvBaGNBoTeh0djQaI1qNCY3WhEZjQqs1odWY0WjNaLVmtBpLeK21oNVa0eos6LQ2tDobOq0VjcYUlXlypJQc3bWDLX9/l4NbN6EzGplyxTVMm3cdtrj+1xe+3lfPY+seY83RNVw8/GJ+NftXWPT9q6XTkR9dnnP6QorSA1SiaGXMmMfYt28vuTnjEEKDEDoQGoTQohE6EFqE0Ia/8IUO8c2i0bXapkej0UeSgh4h9K229Z/z8T2p2eNmzxf5bP/4QyqPFGO2O5i54CYmX3pVv2tBfKOgsoAfrv0h5U3l/HDaD7l13K0DbhK6qSPioh2CMkSoRNFKWuoiCvfnk5o6J9qh9HtSSkoL93I4/2N2vPgMfp+XxBFZXLLkQcadNxedoWdOYfWGdwrf4YkvnyDJnMTLV7zMxMSJ0Q6pW7YcrgFUwlB6n0oUSpfUlBxn34a17F67hrqyUjQ6HWNnz2HSxVeQMnpMv/5VHggF+O/N/81re17jnGHn8OQFT+IwOqIdVrf99qN9wOAeR6H0DypRKB2SUlJ97AhFX21g/5dfUHmkGICMvInMuG4R5c0hLrr00ugG2Qm13lp+tPZHfFn6JbeMvYXvT/s+Oo36+CtKZ6j/KcopAs3NHNuzk0Nfb+HAlo3Ul5cBkDpmLHNvu4vsGbOIiU8AoKo35gvoYQWVBXzvs+9R46nhiXOf4Lrs66IdkqIMKCpRKIRCQSqLD3Fk1w6O7NzOsd07CTT70Or1DB8/ibOvvp6RU88mJi4h2qF2iZSSlftXsuyrZSRZknjlylfIi8+LdliKMuCoRDEE+X1eyg8UcXzfbkr27+H4vt34msLT5sWmpjPhokvJnDSFjLET+s38S11V1lTGExue4PPjnzM7bTbLzls2oK9HKEo0qUQxyAWam6k6UkxF8UHKDxZRemA/VUeKkZHp1ONS0xkzYxYZeRPJGDehX4536AopJe8UvcN/bfovAqEAS6cv5abcm9AMwq7Jj1198n2+FKV3qEQxSAQDAeoryqg+fpTqo0eoOlJM1dHD1JQca0kKRquVlFFjmHHtQlJGjyF1zFjMMfYoR95zCioLWLZpGTsqdzAteRpPnPtEj8362h/lpaoWktI3VKIYQIIBPw1VldSXl1FfUUZtaQm1ZSXUlZZQV15GKBhoKetISiY+YwSjz55JUtZIkjJH4UhK7tfdV7ur1FXK8m3Lef/g+ySYE3ji3CeYP3r+oGxFtPZFYRUAs7MH1rUjZeBRiaKfCAb8NNXV4qqpoam2hsaaKhqrq3DVVNNQWUFDVQWu2poT7v6i0xtwpgwjLi2DUWefQ3xaBnFp6cSnZfS7Sfd6Q7WnmhcKXuDNfW8CcMf4O7hr4l1Y9dYoR9Y3lq8pBFSiUHqfShS9QIZC+DxuvC4XviYXnsYGPK5GPA0N4ccN9Xga6nE31uOuq8NdX4e36dSpvLV6Pba4eOzxiYyYMJmYhCQcSck4k1JwJKdgi43rV9N095UqTxWv7H6FN/a+gS/oY/6o+dw76V6G2YZFOzRFGZRUomjFVVuDr76OyiPFBHw+/D5vePF6afZ+s3bT7PHg93po9njwedw0u9343E343G58bhc+t7vD+z6abDEmj7LtAAAMpklEQVSY7Q4sdjvxGcPJGD8Jq8OJNTYOW1wcVmccMfEJmGPsg/JUUXcdbTzKy7te5p3CdwjIAJdlXsa9k+4ly5EV7dAUZVBTiaKVPz14JwF/Mzv/0nE5rV6PwWzBaLagN5sxWizYE5MwmC2YrDaMVitGixWTLSay2DDb7JjtdkxWGxqttm/+oEFASsmmsk28tuc18o/mo9VomT9qPt8d/12G24dHOzxFGRJUomjlojvuZd/+/UyYPBm9wYjeaEJvMqEzGjGYzOhNJvRGE1qdqrbeVuet470D7/FW4VscrD9IrDGWuybexaKcRSRZkqIdnqIMKeobr5Xxcy+hSugZM2NWtEMZkgKhAOtL1vP+gff55Mgn+EN+JiZM5Ilzn+CKrCsw6Qbm4L/e8uvvTIh2CMoQEZVEIYSIA94EMoFi4AYpZW0b5V4E5gEVUsrxfRmj0ncO1R/irf1v8f7B96nx1uA0OlkwZgHXZ19PTpy6OU97RiXaoh2CMkREq0WxFPhESrlMCLE08vzHbZR7CXgGeKUPY1P6gNvv5pMjn/B24dtsLt+MTuiYkzGHa0Zdw+y02ei1+miH2O+t3l0OwMXjkqMciTLYRStRzAfmRB6/DOTTRqKQUq4VQmT2VVBK7/KH/Gwq3cTfD/2dfx7+J56Ah3RbOg9PeZj5o+eTYFbjAbrij58fBFSiUHpftBJFspSyFEBKWSqEUFcnB6nmYDNflX3F6sOr+eTIJ9T56ojRx3Bl1pVcPepqzko6a9CPoFaUgU7IDvr7n9EbC7EaSGlj10+Bl6WUzlZla6WUse28TybwwemuUQghlgBLAJKTk6e+8cYb3Yrb5XJhs6lzv53VVn3VB+rZ493DLs8u9nj24JM+jMLIBPMEzrKexVjzWPRiaJ5a6snP1282egB4dIa5R96vP1L/H7vmTOpr7ty5W6SU09ra12stCinlxe3tE0KUCyGGRVoTw4CKHjjeCmAFwLRp0+ScOXO69T75+fl097VDUX5+PtPOncbWiq1sKtvEupJ1FNaGp5ZIMidxTfY1zMmYw4xhMzBqjVGONvp68vP13L4NAMyZM3hvhar+P3ZNb9VXtE49vQfcBiyLrP8WpTiUbqhwV/B1xddsq9jG2tK1HH/jOEEZRK/Rc1bSWTw85WFmpc1iTOwYdVpJUQaBaCWKZcBfhRB3AEeAhQBCiFTgBSnllZHnrxO+6J0ghDgG/EJK+afohDw0VXmq2Fezjz01e9hZtZOCqgIq3OEGoFFrJEOXwZ0T7uTslLOZlDhJjXXoQ/+zaHK0Q1CGiKgkCillNXBRG9tLgCtbPb+pL+MaylzNLg7VH+JA/QGKaosoqitif+1+Kj2VLWWGxwxnavJUxsePZ3LSZMbGjWXd5+uYc9ac6AU+hKU6B++1CaV/USOzhxBXs4tjrmMcbzzOkcYjHG44HF7XH6bC8+1lIoPGwCjnKGamziQnNoex8WPJicvBbhg8NzkaDN7fXgLA1ZNSoxyJMtipRDEISClx+V1UeiqpcldR7i4PL03llLnLKHWVUtpUSkNzwwmvizXGMtw+nHNSzyHLkcVIx0iyHFkMjxmOVqMmLuzvXvvyMKAShdL7VKLoh/whP43NjdT76mlobqDOW0edL7zUemup9dVS462hxlNDtbeaak813qD3lPeJMcSQYk0h1ZrK5KTJpNpSSbelkx4TXlQLQVGUzlCJ4gyFZAhf0Icv4MMb9OINePEGvXgCHjx+D56AB3fAHV773TQFmmjyf7u4/C6amsPrhuYGGpsb8QQ87R5PJ3TEmmKJNcUSZ4pjuH048aZ4EswJJFoSSTQnkmhJJNmSjEU/+O9ypyhK71OJopXH1j3G4crDvPvpuzQHm/GH/CesfUEf/pA/nBgiyaE51Nzl41h0Fqx6K1a9FZvehtVgJcGcQIwhhhhDDDaDDYfBgcPowG6w4zQ6cRqdOEwOYvQx6mZGiqL0KZUoWimoKgj/qm9sRCd0GLQGDFoDZp255bFRa0Sv0WPSmVqem7QmTDoTRq0Rs86MSWfCrDVj1psx68KLVW/ForNg0pnU2AJFUQYUlShaeWf+O2okqDJgPHfL1GiHoAwRKlEoygAVZzVEOwRliFDnQBRlgFq5+SgrNx+NdhjKEKAShaIMUKu2HGPVlmPRDkMZAlSiUBRFUTqkEoWiKIrSIZUoFEVRlA6pRKEoiqJ0SHWPVZQB6qXbp0c7BGWIUIlCUQYos0HN8Kv0DXXqSVEGqFc3FPPqhuIoR6EMBSpRKMoA9cGOUj7YURrtMJQhQCUKRVEUpUNRSRRCiDghxD+FEIWRdWwbZTKEEJ8KIfYIIXYJIf4tGrEqiqIMddFqUSwFPpFSZgOfRJ6fLAB8X0o5FjgHuF8IMa4PY1QURVGIXqKYD7wcefwycO3JBaSUpVLKrZHHjcAeIK3PIlQURVEAEFLKvj+oEHVSSmer57VSylNOP7XanwmsBcZLKRvaKbMEWBJ5mgPsizx2APWtip7ueQJQ1ak/pOtOPlZPvqajcu3ta2v76bap+uraNlVfXd/W+rmqr76rrxFSysQ290gpe2UBVgM721jmA3Unla3t4H1swBbgO92MY0UXn2/uxTpZ0Vuv6ahce/va2n66baq+VH31Zn21UX+qvvpBffXagDsp5cXt7RNClAshhkkpS4UQw4CKdsrpgbeAP0sp3+5mKO938Xlv6s6xOvuajsq1t6+t7afbpuqra9tUfXV9W1/VmaqvTorWqaf/AqqllMuEEEuBOCnlj04qIwhfv6iRUj7ch7FtllJO66vjDXSqvrpG1VfXqPrqmt6qr2hdzF4GXCKEKAQuiTxHCJEqhPgwUmYWcCtwoRDi68hyZR/EtqIPjjGYqPrqGlVfXaPqq2t6pb6i0qJQFEVRBg41MltRFEXpkEoUiqIoSodUolAURVE6pBJFJwkhxgohnhdCrBJC3BvteAYCIcS1Qog/CiH+JoS4NNrx9HdCiJFCiD8JIVZFO5b+SghhFUK8HPlc/Uu04+nveuozNSQShRDiRSFEhRBi50nbLxdC7BNCFEW66bZLSrlHSnkPcAMw6Lvr9VCdvSulvAtYDCzqxXCjrofq66CU8o7ejbT/6WLdfQdYFflcXdPnwfYDXamvnvpMDYlEAbwEXN56gxBCCzwLXAGMA24SQowTQkwQQnxw0pIUec01wBeEJzIc7F6iB+os4meR1w1mL9Fz9TXUvEQn6w5IB45GigX7MMb+5CU6X189YkjcClVKuTYyX1Rr04EiKeVBACHEG8B8KeVvgHntvM97wHtCiL8Df+m9iKOvJ+osMmhyGfB/MjLB42DVU5+xoagrdQccI5wsvmbo/NA9QRfra3dPHHNIVnREGt/+MoHwB7Dd2WmFEHOEEE8LIf4AfNheuUGuS3UGPAhcDCwQQtzTm4H1U139jMULIZ4HzhJCPNrbwfVz7dXd28D1Qojn6NvpUfq7Nuurpz5TQ6JF0Q7RxrZ2Rx9KKfOB/N4KZoDoap09DTzde+H0e12tr2pgKCbUtrRZd1LKJuD2vg5mAGivvnrkMzWUWxTHgIxWz9OBkijFMlCoOusaVV/dp+qua3q1voZyotgEZAshsoQQBuBG4L0ox9TfqTrrGlVf3afqrmt6tb6GRKIQQrwObAByhBDHhBB3SCkDwAPAPwjfPe+vUspd0YyzP1F11jWqvrpP1V3XRKO+1KSAiqIoSoeGRItCURRF6T6VKBRFUZQOqUShKIqidEglCkVRFKVDKlEoiqIoHVKJQlEURemQShSK0sOEEMVCiIQzLaMo/YVKFIqiKEqHVKJQlDMghHhXCLFFCLFLCLHkpH2ZQoi9InxHth0ifHdES6siDwohtgohCoQQuZHXTBdCrBdCbIusc/r0D1KUNqhEoShn5rtSyqmE73r4kBAi/qT9OcAKKeVEoAG4r9W+KinlFOA54AeRbXuB86WUZwGPAb/u1egVpRNUolCUM/OQEGI78CXh2TuzT9p/VEq5LvL4NWB2q31vR9ZbgMzIYwewMnKby/8B8nojaEXpCpUoFKWbhBBzCN+YaaaUchKwDTCdVOzkydRaP/dF1kG+vTfML4FPpZTjgavbeD9F6XMqUShK9zmAWimlO3KN4Zw2ygwXQsyMPL6J8D3XT/eexyOPF/dIlIpyhlSiUJTu+wjQCSF2EG4JfNlGmT3AbZEycYSvR3Tkt8BvhBDrAG1PBqso3aWmGVeUXiKEyAQ+iJxGUpQBS7UoFEVRlA6pFoWiKIrSIdWiUBRFUTqkEoWiKIrSIZUoFEVRlA6pRKEoiqJ0SCUKRVEUpUMqUSiKoigd+v8K6nDpqI8legAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "alphas1, coeffs, _ = sklearn.linear_model.lasso_path(X, y, alphas=alphas)\n", "\n", "# Plot the paths of the coefficients\n", "plt.semilogx(alphas1,coeffs.T)\n", "plt.grid()\n", "plt.legend(names_x, loc='upper right')\n", "\n", "\n", "# Plot a line on the optimal alpha\n", "plt.semilogx([alpha_opt,alpha_opt], [-0.2,0.6], '--')\n", "plt.ylim([-0.2,0.6])\n", "plt.xlabel('alpha')\n", "plt.ylabel('coeff')\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reference:\n", "\n", "- DS-GA 1003 Machine Learning Spring 2020" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 2 }