notebook/Chapter_6_Linear_Model_Selection_and_Regularization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Model Selection and Regularization\n",
    "\n",
    "Topics covered in this chapter of the book-\n",
    "\n",
    "* Subset Selection Models\n",
    "  * Best Subset Selection\n",
    "  * Forward Stepwise Selection \n",
    "  * Backward Stepwise Selection \n",
    "* Shrinkage Methods\n",
    "  * Ridge Regression\n",
    "  * The Lasso \n",
    "* Dimension Reduction Methods- PCR and PLS Regression \n",
    "  * Principal Components Regression\n",
    "  * Partial Least Squares\n",
    "\n",
    "**Following is the summary of concepts along with data and python code-**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import pandas as pd \n",
    "import math\n",
    "import time\n",
    "import random\n",
    "import itertools\n",
    "\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "from sklearn import datasets, linear_model\n",
    "from sklearn.decomposition import PCA \n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.cross_decomposition import PLSRegression, PLSSVD\n",
    "from sklearn import linear_model \n",
    "from sklearn.model_selection import cross_val_predict \n",
    "from sklearn.metrics import mean_squared_error, r2_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter, we discuss some ways in which the simple linear model can be improved, by replacing plain least squares fitting with some alternative fitting procedures.\n",
    "\n",
    "There are many alternatives, both classical and modern, to using least squares to fit. In this chapter, we discuss three important classes of methods- \n",
    "* **Subset Selection**- This approach involves identifying a subset of the p predictors that we believe to be related to the response. We then fit a model using least squares on the reduced set of variables.\n",
    "\n",
    "* **Shrinkage**- This approach involves fitting a model involving all p pre- dictors. However, the estimated coefficients are shrunken towards zero relative to the least squares estimates. This shrinkage (also known as regularization) has the effect of reducing variance. Depending on what type of shrinkage is performed, some of the coefficients may be esti- mated to be exactly zero. Hence, shrinkage methods can also perform variable selection.\n",
    "\n",
    "* **Dimension Reduction**- This approach involves projecting the p predic- tors into a M-dimensional subspace, where M < p. This is achieved by computing M different linear combinations, or projections, of the variables. Then these M projections are used as predictors to fit a linear regression model by least squares."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Subset Selection Methods\n",
    "\n",
    "### Best Subset Selection\n",
    "To perform best subset selection, we fit a separate least squares regression for each possible combination of the    predictors. The process for selecting best model looks like this-\n",
    "\n",
    "1. Let M0 denote the null model, which contains no predictors. This model simply predicts the sample mean for each observation.\n",
    "\n",
    "\n",
    "2. For k=1,2,...p:\n",
    "* (a) Fit all models that contain exactly k predictors.\n",
    "* (b) Pick the best among these k models, and call it Mk. Here *best* is defined as having the smallest RSS, or equivalently largest R .\n",
    "\n",
    "3. Select a single best model from among M0, . . . , Mp using cross- validated prediction error, Cp (AIC), BIC, or adjusted R ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(322, 20)\n",
      "   AtBat  Hits  HmRun  Runs  RBI  Walks  Years  CAtBat  CHits  CHmRun  CRuns  \\\n",
      "0    293    66      1    30   29     14      1     293     66       1     30   \n",
      "1    315    81      7    24   38     39     14    3449    835      69    321   \n",
      "2    479   130     18    66   72     76      3    1624    457      63    224   \n",
      "3    496   141     20    65   78     37     11    5628   1575     225    828   \n",
      "4    321    87     10    39   42     30      2     396    101      12     48   \n",
      "\n",
      "   CRBI  CWalks League Division  PutOuts  Assists  Errors  Salary NewLeague  \n",
      "0    29      14      A        E      446       33      20     NaN         A  \n",
      "1   414     375      N        W      632       43      10   475.0         N  \n",
      "2   266     263      A        W      880       82      14   480.0         A  \n",
      "3   838     354      N        E      200       11       3   500.0         N  \n",
      "4    46      33      N        E      805       40       4    91.5         N  \n"
     ]
    }
   ],
   "source": [
    "Hitters = pd.read_csv('../data/Hitters.csv', header=0, na_values='NA')\n",
    "print(Hitters.shape) # get the dimension of this \n",
    "print(Hitters.head()) # pull a sample of this dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(263, 20)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "Hitters.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = Hitters.Salary  # the response variable \n",
    "\n",
    "# change categorical variables into dummies\n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']]) \n",
    "\n",
    "#keep n-1 dummy variables for each categorical variable with n unique value/category\n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Functions to run on a subset of feature and extract RSS\n",
    "def getRSS(y, X, feature_list):\n",
    "    model = sm.OLS(y, X[list(feature_list)]).fit()\n",
    "    RSS = ((model.predict(X[list(feature_list)]) - y) ** 2).sum()\n",
    "    return {'Model':model, \"RSS\":RSS}\n",
    "\n",
    "## Functions to select the best model on basis of RSS for given number of features\n",
    "def bestModel(y, X, K):\n",
    "    results = []\n",
    "    for c in itertools.combinations(X.columns, K):\n",
    "        results.append(getRSS(y, X, c))     \n",
    "    model_all =  pd.DataFrame(results)\n",
    "    \n",
    "    best_model = model_all.loc[model_all[\"RSS\"].argmin()]\n",
    "    return best_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_feature = 3\n",
    "models = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "for i in range(1,(max_feature+1)):   \n",
    "    models.loc[i] = bestModel(y, X, i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            RSS                                              Model\n",
      "1  4.321393e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "2  3.073305e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "3  2.943854e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "                                 OLS Regression Results                                \n",
      "=======================================================================================\n",
      "Dep. Variable:                 Salary   R-squared (uncentered):                   0.772\n",
      "Model:                            OLS   Adj. R-squared (uncentered):              0.769\n",
      "Method:                 Least Squares   F-statistic:                              292.7\n",
      "Date:                Fri, 23 Oct 2020   Prob (F-statistic):                    5.01e-83\n",
      "Time:                        23:48:21   Log-Likelihood:                         -1902.0\n",
      "No. Observations:                 263   AIC:                                      3810.\n",
      "Df Residuals:                     260   BIC:                                      3821.\n",
      "Df Model:                           3                                                  \n",
      "Covariance Type:            nonrobust                                                  \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Hits           2.3164      0.318      7.295      0.000       1.691       2.942\n",
      "CRBI           0.6665      0.065     10.293      0.000       0.539       0.794\n",
      "PutOuts        0.2614      0.077      3.381      0.001       0.109       0.414\n",
      "==============================================================================\n",
      "Omnibus:                      114.300   Durbin-Watson:                   2.008\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              666.384\n",
      "Skew:                           1.652   Prob(JB):                    1.98e-145\n",
      "Kurtosis:                      10.063   Cond. No.                         8.59\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] R² is computed without centering (uncentered) since the model does not contain a constant.\n",
      "[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "# The summary() function also returns R2, RSS, adjusted R2, Cp, and BIC.\n",
    "# We can examine these on each model to try to select the best overall model.\n",
    "print(models)\n",
    "print(models.loc[3, 'Model'].summary()) # model with 3 variables- Hits, CRBI, PutOuts is the best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAERCAYAAAB2CKBkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAneUlEQVR4nO3deZxVdf3H8ddn9g1mZx1g2BSRZYCLmVqmZvIzEzMXVAySxMwsq1+LLaZWvzJLzcyFhB+I5poVaZmm+MtcmZFFQURkUXBhhp1hh8/vj3sGh2GGmYG598yd+34+HvfRufd8770fTsf7me/5fs/3Y+6OiIgkr5SwAxARkXApEYiIJDklAhGRJKdEICKS5JQIRESSnBKBiEiSS8hEYGbTzGyNmb3egrY3m9m84LHEzDbEIUQRkYRhiXgfgZl9EtgC3OPuQ1rxviuBEe5+ScyCExFJMAnZI3D3fwPr6r9mZv3N7AkzqzKz58xsUCNvvQC4Py5BiogkiLSwA2hDU4CvuPtbZvYx4Hbg5LqdZtYH6As8E1J8IiLtUodIBGaWBxwHPGxmdS9nNmg2DnjE3ffEMzYRkfauQyQCope4Nrh7xUHajAOuiE84IiKJIyHHCBpy903AcjM7F8CihtftD8YLCoEXQwpRRKTdSshEYGb3E/1RP9LMVpnZJOAiYJKZzQcWAmPrvWUc8IAn4hQpEZEYS8jpoyIi0nYSskcgIiJtJ+EGi0tKSry8vDzsMEREEkpVVVWNu5c2ti/hEkF5eTmVlZVhhyEiklDMbGVT+3RpSEQkySkRiIgkOSUCEZEkp0QgIpLklAhERJKcEoGISJJTIhARSXJJkwiW19RywxOL0ZIaIiL7S5pE8NSiD7jj2be59emlYYciItKuJNydxYfq0k/0480PtnDzv5ZQXpLD2IqeYYckItIuJE2PwMz4xdlDOaZvEd95ZAFVK9c1/yYRkSSQNIkAICMthbvGj6JHfhaT76ni3XVbww5JRCR0SZUIAApzM5g2cTS79zpfmj6Hjdt2hR2SiEioki4RAPQrzePO8aNYUVPL1/74Krv27A07JBGR0CRlIgD4eP9i/ufsoTz3Vg0/mbVQ00pFJGklzayhxpwX6cXymlruePZt+pXk8uVP9As7JBGRuEvqRADwnc8cyYqaWn7+9zcoL87l04O7hh2SiEhcJe2loTopKcZN51UwtGc+X39gLgvf2xh2SCIicZX0iQAgOyOVu78YoSA7nUnTK/lw0/awQxIRiRslgkCXzlncPWE0m7fvYtKMOWzduTvskERE4kKJoJ7BPTrzuwtHsOi9TVz1wDz27tVMIhHp+JQIGjh5UFd+fMZgnlz0ITc8sTjscEREYi7pZw01ZuJx5SyrruWufy+jb0ku447pHXZIIiIxE/MegZmlmtlcM3uskX3fMrNFZrbAzJ42sz6xjqclzIyffG4wJx5Ryo/+8jrPL60JOyQRkZiJx6WhbwBvNLFvLhBx92HAI8Cv4hBPi6SlpnDbhSPoX5rHV+6tYumaLWGHJCISEzFNBGZWBnwWuLux/e4+293rlgB9CSiLZTyt1SkrnakTI2SmpXDJ9Dmsq90ZdkgiIm0u1j2CW4DvAi1Z1W0S8I/GdpjZZDOrNLPK6urqNgyveWWFOUz5YoQPN23nspmV7Ni9J67fLyISazFLBGZ2BrDG3ata0HY8EAFubGy/u09x94i7R0pLS9s40uaN7F3Ib84bzpwV6/n+n17TAnUi0qHEctbQ8cCZZnY6kAV0NrN73X18/UZm9mngh8CJ7r4jhvEcljOG9WBFTS2/fnIJfUty+fopA8MOSUSkTcSsR+DuV7t7mbuXA+OAZxpJAiOAu4Az3X1NrGJpK1ecNICzR/bkpqeWMGv+e2GHIyLSJuJ+Q5mZXW9mZwZPbwTygIfNbJ6ZzYp3PK2xr+5xeRH//fB8qlauDzskEZHDZol2vTsSiXhlZWWoMayr3cnnb3+eLdt385crjqdXUU6o8YiINMfMqtw90tg+LTFxCIqCuse79uzlkulz2LRddY9FJHEpERyi/qV53HnxKJbX1HLFfa+yW3WPRSRBKREchuP6l/Dzzw/hubdquPZvqnssIolJi84dpvNH92ZZTS13/d8y+pXkcckJfcMOSUSkVZQI2sD3ThvEypqt/PTxRfQpzuGUo1T3WEQShy4NtYGUFOPm8ysY0iOfK++fy6L3NoUdkohIiykRtJHsjFTunhAhPzudSTPmqO6xiCQMJYI21LVzFlMnjGbjtl18eUal6h6LSEJQImhjg3t05ncXjGDhexv55oOqeywi7Z8SQQycclRXfvjZwfxz4Yfc8E/VPRaR9k2zhmLkkuPLWV6zJZhWmsv5o1X3WETaJyWCGDEzrv3c0axcu5Uf/vl1ehXmcNyAkrDDEhE5gC4NxVBaagq/v2gkfUty+cq9VbxdrbrHItL+KBHEWOesdKZNHE2G6h6LSDulRBAHvYpyuOviCO9v3M5XZlap7rGItCtKBHEyqk8hvzl3OK+sWMfVqnssIu2IBovj6HPDe7C8ppabnlpCv9Jcvnay6h6LSPiUCOLsypMHsLymll8/uYTyklzOGNYj7JBEJMnp0lCcmRm//MJQRpcX8q2H5vPqO6p7LCLhinkiMLNUM5trZo81si/TzB40s6Vm9rKZlcc6nvYgMy2Vuy6O0K1zFpPvqeTddVvDDklEklg8egTfAN5oYt8kYL27DwBuBm6IQzztQl3d45279zJphuoei0h4YpoIzKwM+CxwdxNNxgIzgu1HgFPMzGIZU3syoEsed4wfxbLqWr72x7mqeywioYh1j+AW4LtAU79wPYF3Adx9N7ARKG7YyMwmm1mlmVVWV1fHKNRwHD+ghJ+dNYR/L6nmur8t0rRSEYm7mCUCMzsDWOPuVYf7We4+xd0j7h4pLS1tg+jal3HH9OayT/Zj5ksrmf7CirDDEZEkE8sewfHAmWa2AngAONnM7m3QZjXQC8DM0oB8YG0MY2q3vjdmEJ8Z3JWfPraIZxZ/GHY4IpJEYpYI3P1qdy9z93JgHPCMu49v0GwWMCHYPidok5TXRlJSjFvGVTC4R2eu/KPqHotI/MT9PgIzu97MzgyeTgWKzWwp8C3g+/GOpz3JyUhj6oTRdMqK1j1eo7rHIhIHlmh/gEciEa+srAw7jJh6ffVGzrvrRQZ0yePByR8nOyM17JBEJMGZWZW7RxrbpzuL26EhPfO5ddwIXlutusciEntKBO3Upwd35YenH8UTCz/gxiffDDscEenAtOhcOzbphL4sr6nljmffpm9JLudFeoUdkoh0QEoE7ZiZce2ZR/POuq384NHXKCvM5rj+qnssIm1Ll4baufTUFG67MFr3+PJ7X2WZ6h6LSBtTIkgA+dnRusdpKcYl0+ewXnWPRaQNKREkiF5FOUz54ije27idy+5V3WMRaTtKBAlkVJ8ibjxnGK8sX8cPHn1dC9SJSJvQYHGCGVvRkxU1W7n5X9G6x1ecNCDskEQkwSkRJKCvnzKA5TVbuPGfb1JenMtnh3UPOyQRSWC6NJSAonWPhxHpU8i3HprHXNU9FpHDoESQoLLSU7nr4lF07ZzFpfdUsWq96h6LyKFRIkhgxXmZTJsYYcfuPUyaXslm1T0WkUOgRJDgBnTpxB0XjWJp9RauvF91j0Wk9ZQIOoATBkbrHj/7ZjU/fWxR2OGISILRrKEO4oJjerOsegt/eG45fUtymXh837BDEpEEoUTQgXz/v45ixdqtXP/YIvoU53LSoC5hhyQiCUCXhjqQ1BTjt+MqOKp7Z772x1d5433VPRaR5ikRdDB1dY/zstKYNH0Oazar7rGIHJwSQQfULT+LqRNGs37rLi6dUcm2nVqgTkSaFrNEYGZZZvaKmc03s4Vmdl0jbXqb2Wwzm2tmC8zs9FjFk2yG9Mzn1gtGsGD1Rr79sOoei0jTYtkj2AGc7O7DgQpgjJkd26DNj4CH3H0EMA64PYbxJJ1Tg7rHf3/tA36tusci0oSYzRry6BrJdeW00oNHwz9LHegcbOcD78UqnmQ16YS+vF1dy+1B3eNzVfdYRBqI6RiBmaWa2TxgDfCUu7/coMm1wHgzWwX8Hbiyic+ZbGaVZlZZXV0dy5A7HDPj+rFHc8KAEn7w59d4adnasEMSkXYmponA3fe4ewVQBhxjZkMaNLkAmO7uZcDpwEwzOyAmd5/i7hF3j5SWlsYy5A4pPTWF3180kj7FuVw2s0p1j0VkP3GZNeTuG4DZwJgGuyYBDwVtXgSygJJ4xJRs8rPTmTZhNKkpxqQZlWzYqrrHIhIVy1lDpWZWEGxnA6cCixs0ewc4JWhzFNFEoGs/MdK7OIcpF49i9fptXDazip27tUCdiMS2R9AdmG1mC4A5RMcIHjOz683szKDNt4FLzWw+cD8w0VWIN6Yi5UXceO4wXl6+jh/8+TXVPRaRmM4aWgCMaOT1a+ptLwKOj1UM0rixFT1ZVl3Lb59+i36luXz1U6p7LJLMtOhckrrq0wNZXlPLr56I1j0+fajqHoskKy0xkaTMjF+dM4xRfQr55oPzmPfuhrBDEpGQKBEksaz0VKZcPIounTP58oxKVm/YFnZIIhICJYIkV5yXybQJo9mxaw+Tps9R3WORJKREIAzs2onbx4/krTVb+LrqHoskHSUCAeATA0v56dghzH6zmp89/kbY4YhIHGnWkOxz4ceidY/v/k+07vGE48rDDklE4kCJQPZz9enRusfX/W0hvYtzOOlI1T0W6eh0aUj2U1f3eFC3zlz5x7ks/kB1j0U6OiUCOUBuZhpTJ0bIzUxl0vRK1T0W6eCUCKRR3fOzmTphNOtqd3LpPVVs36W6xyIdlRKBNGlIz3xuGVfBglUb+PZD81X3WKSDOmgiMLPPmVmfes+vCYrRzzKzvrEPT8J22tHduPq/BvH4a+9z01NLwg5HRGKguR7BzwnqA5jZGcB44BJgFnBnbEOT9uLST/TjgmN6cdvspTxStSrscESkjTWXCNzdtwbbZwNT3b3K3e8GVDMySUTrHg/h+AHFXP3oAl5W3WORDqW5RGBmlhfUET4FeLrevqzYhSXtTXpqCrdfOIreRTlcdm8Vy2tqww5JRNpIc4ngFmAeUAm84e6VAGY2Ang/ppFJu5Ofk860iaMxYNL0Oap7LNJBHDQRuPs04ESiReZPr7frfeBLMYxL2qk+xblM+WKEVeu38ZV7VfdYpCNobtZQH2CLu891971mdpKZ/Ra4EPggLhFKuzO6vIhfnTOMl5at40d/Ud1jkUTX3KWhh4BcADOrAB4G3gGGA7fHNDJp184a0ZOvnzKQhypXcef/LQs7HBE5DM0tOpft7u8F2+OBae7+m2DweN7B3mhmWcC/gczgex5x95800u484FrAgfnufmGr/gUSmm8GdY9veGIxfUtyGDNEdY9FElGzs4bqbZ9MMGvI3VtyYXgHcLK7DwcqgDFmdux+H242ELgaON7djwaualnY0h6YGTeeM4yRvQu46sF5LFi1IeyQROQQNJcInjGzh4JxgULgGQAz6w4cdMqIR20JnqYHj4YXky8Ffu/u64P3rGll/BKyrPRUpnwxQkleJpNmVPKe6h6LJJzmEsFVwKPACuAEd68raNsN+GFzH25mqWY2D1gDPOXuLzdocgRwhJk9b2YvmdmYJj5nsplVmllldXV1c18rcVaSl8n/ThzN9p17uGT6HLbs2B12SCLSCs1NH3V3f8Ddb3b31fV2zQdKmvtwd9/j7hVAGXCMmQ1p0CQNGAh8CrgA+IOZFTTyOVPcPeLukdJS3dDcHg3s2onfX/RR3eM9WqBOJGE0N320s5ldbWa3mdlnLOpKYBlwXku/xN03ALOBhn/xrwJmufsud18OLCGaGCQBffKIUq4782ieWbyGnz2+KOxwRKSFmrs0NBM4EngN+DLRH/NzgLPcfezB3mhmpXV/3ZtZNnAqsLhBs78Q7Q1gZiVELxVpLmICG39sHyad0Jf/fX4FM19cEXY4ItICzU0f7efuQwHM7G6idxT3dveWlKzqDswws1SiCechd3/MzK4HKt19FvBP4DNmtgjYA3zH3bWiWYL7welHsXJtLdf+bRG9inL4lOoei7RrdrC7Qs3sVXcf2dTzMEQiEa+srAwzBGmB2h27OffOF3ln3Vb+dPlxHNmtU9ghiSQ1M6ty90hj+5q7NDTczDYFj83AsLptM1NVc2lSXd3jnIxULpk+h+rNO8IOSUSa0NysoVR37xw8Orl7Wr3tzvEKUhLT/nWPK1X3WKSdUs1iiamhZfncfH4F81dt4L8fVt1jkfZIiUBibsyQbnx/zCAeW/A+N/9LdY9F2pvmZg2JtInJn+zHsupafvfMUsqLc/nCqLKwQxKRgHoEEhdmxk/PGsJx/Yv5/qMLeGX5urBDEpGAEoHETUZaCndcNIpeRTlcNrOSFap7LNIuKBFIXOXnpDNtwmgALpk+h41bdzXzDhGJNSUCibvyklzuujjCu+u3qu6xSDugRCChOKZvETd8YRgvLlvLj//yuuoei4RIs4YkNGePLGNFTS23PrOUfqW5XHZi/7BDEklKSgQSqm+eegTLamr55ROL6VOcy5gh3cIOSSTp6NKQhMrM+PW5w6noVcBVD87ltVUbww5JJOkoEUjostJTmXJxhOLcTCbNmMP7G1X3WCSelAikXSjtlMn/fmk023bu4ZLpldSq7rFI3CgRSLtxRNdO3HbRSJZ8uFl1j0XiSIlA2pUTjyjl2jOP5unFa/j542+EHY5IUtCsIWl3Lj62D8uqtzDt+eX0Lc3l4mP7hB2SSIemRCDt0o8+O5iVa7dy7ayF9C7K4cQjSsMOSaTDitmlITPLMrNXzGy+mS00s+sO0vYLZuZm1mg9TUk+qSnGrReMYGCXPL5236ss+XBz2CGJdFixHCPYAZzs7sOBCmCMmR3bsJGZdQK+Abwcw1gkAeVlpjFt4miygrrHNVtU91gkFmKWCDxqS/A0PXg0Ng3kp8ANwPZYxSKJq0dBNlMnRKjZskN1j0ViJKazhsws1czmAWuAp9z95Qb7RwK93P3xZj5nsplVmllldXV17AKWdmlYWQG3nF/B3Hc28J1HFmiBOpE2FtNE4O573L0CKAOOMbMhdfvMLAW4Cfh2Cz5nirtH3D1SWqpBw2Q0Zkh3vjdmEH+b/x43/+utsMMR6VDich+Bu28AZgNj6r3cCRgCPGtmK4BjgVkaMJamfOXEfpwXKePWp9/iz3NXhR2OSIcRy1lDpWZWEGxnA6cCi+v2u/tGdy9x93J3LwdeAs5098pYxSSJzcz42VlD+Xi/Yr73yGvMWaG6xyJtIZY9gu7AbDNbAMwhOkbwmJldb2ZnxvB7pQPLSEvhjvEjKSvMZvI9laxcq7rHIofLEm3gLRKJeGWlOg3JbkVNLWfd/jzFuRk8evnx5Oekhx2SSLtmZlXu3uild601JAmpvCSXu8aP4p11W7n8vip27VHdY5FDpUQgCetj/Yr55dnDeOFt1T0WORxaa0gS2hdGlbG8ppbbZkfrHk/+pOoei7SWEoEkvG+degTL19byi39E6x6fdrTqHou0hi4NScJLSTF+c+5whpcVcNUD83h9teoei7SGEoF0CFnpqfzhixGKcjNU91iklZQIpMMo7ZTJtImjqd2xh0mqeyzSYkoE0qEc2a0Tt104gsUfbOIbD8xT3WORFlAikA7nU0d24dozj+Zfb3zIL/6uuscizdGsIemQvvjxcpZV13L3f6J1jy/6mOoeizRFiUA6rB+fMZiVa2u55q/RusefGKglzEUao0tD0mGlphi/u3AkA7vk8dV7X+Ut1T0WaZQSgXRoeZlpTJ04msz0VC6ZobrHIo1RIpAOr2dQ97h68w4mq+6xyAGUCCQpDO9VwE3nVfDqOxv4ruoei+xHiUCSxulDu/PdMUcya/573KK6xyL7aNaQJJXLT+zP8upafvv0W/QtyeWsET3DDkkkdOoRSFIxM37++aEc26+I7z6ygErVPRZRIpDkk5GWwp3jR9GzMJvJM6t4Z+3WsEMSCZUSgSSlgpwMpk0czV53vjT9FTZu2xV2SCKhiVkiMLMsM3vFzOab2UIzu66RNt8ys0VmtsDMnjYzrQMgcdO3JJc7g7rHV9z3quoeS9KKZY9gB3Cyuw8HKoAxZnZsgzZzgYi7DwMeAX4Vw3hEDnBsv2J+cfYw/rO0hmv+ulDTSiUpxSwReNSW4Gl68PAGbWa7e90F2peAsljFI9KUc0aVccVJ/bn/lXeY+p/lYYcjEncxHSMws1QzmwesAZ5y95cP0nwS8I8mPmeymVWaWWV1dXUMIpVk9+1Tj+T0od34+d/f4MmFH4QdjkhcxTQRuPsed68g+pf+MWY2pLF2ZjYeiAA3NvE5U9w94u6R0lKtICltLyXFuOm8CoaVFfAN1T2WJBOXWUPuvgGYDYxpuM/MPg38EDjT3bUimIQmWvd41L66xx9s3B52SCJxEctZQ6VmVhBsZwOnAosbtBkB3EU0CayJVSwiLdWlUxZTJ0aidY9nzFHdY0kKsewRdAdmm9kCYA7RMYLHzOx6MzszaHMjkAc8bGbzzGxWDOMRaZFB3TrzuwtH8Mb7qnssycESbbpcJBLxysrKsMOQJDDjhRX8ZNZCLv1EX3742cFhhyNyWMysyt0jje3TonMiTZhwXDnLqrfwh+eW07ckjws/1jvskERiQolA5CB+fMZgVq7byo//+jq9i3I4YWBJ2CGJtDmtNSRyEGmpKfzughEMKM3j8vuqWLpGdY+l41EiEGlGp6x0pk6MkJmWypemz2Gt6h5LB6NEINICZYU53D0hwppNO7hsZpXqHkuHokQg0kIVvQq4+fwKKleu53t/Ut1j6TiUCERa4fSh3fnOaUfy13nvcevTS8MOR6RNaNaQSCt99VP9WVZdy83/WkJ5SQ5jK1T3WBKbegQirWRm/OLsoXysbxHfeWQBVStV91gSmxKByCHYV/e4IJvJ91Tx7jrVPZbEpUQgcogKczOYOiHC7r3Ol6bPUd1jSVhKBCKHoV9pHneOH8XKtbVcfm8Vzyz+kCUfbtaqpZJQNFgscpg+3r+Y//n8UL73pwW88Pbafa8X5KTTsyCbssJsehbk0LMwe9/zssJs8rPTMbMQIxeJUiIQaQPnRnpx0qAurFy7ldUbtrF6/TZWrY9uL6uu5bm3ati6c/+b0HIzUulZmE1ZYQ49C7L3JYqeQaIozctUopC4UCIQaSMleZmU5GUyqk/hAfvcnfVbd7F6/TZWb9jKqvXbWLV+276kUbVy/QFjDBlpKdHEsK9X8VGyKCvKoWunTNJSdXVXDp8SgUgcmBlFuRkU5WYwtCy/0Tabt+/alxg+6lVsY9WGbfzrjTXUNFjjKDXF6NY5a18PomxfosihrDCb7gVZZKalxuOfJwlOiUCkneiUlc6gbukM6ta50f3bd+1pJFFELz+99PZaPti0nYbF1Lp0yjzg8lNZXQ+jMJucDP0EiBKBSMLISk+lf2ke/UvzGt2/a89ePti4fb9LTnWJYsGqDTzx+vvs2rN/pijMSQ+SQ84BYxRlBTl0zk7TOEUSUCIQ6SDSU1PoVZRDr6KcRvfv3eus2bxj3xjF6g3BOMX6bSyt3sL/LalmW4NVVfMy0/ZLDvuNUxTmUJKXoUTRASgRiCSJlBSjW34W3fKzGNXnwP3uzrranftdfqo/qF25Yh2btu9/f0Rm3YD2AYkiOk7RtXMWqSlKFO1dzBKBmWUB/wYyg+95xN1/0qBNJnAPMApYC5zv7itiFZOINM3MKM7LpDgvk2FlBY222bQ9mPlUd/lpQ3D5af02nnp/EzVbdu7XPi1IPnU9iJ71BrXLCrPpnp9NRppmPoUtlj2CHcDJ7r7FzNKB/5jZP9z9pXptJgHr3X2AmY0DbgDOj2FMInIYOmel07l7Okd1b3xAe9vOPfsSRP2psqvXb+OFt2v4YNN26pdxMAsGtOslirpeRa+gZ5GdoZlPsRazRODRqh1bgqfpwaNhJY+xwLXB9iPAbWZmroofIgkpOyOVAV3yGNCl8QHtnbuDAe0NW/e7/LR6/TbmvbuBfzQyoF2Um9H4vRRB4sjPTo/HP61Di+kYgZmlAlXAAOD37v5ygyY9gXcB3H23mW0EioGaBp8zGZgM0Lt371iGLCIxlJGWQu/iHHoXNz6gvWevs2bz9kbHKJZ8uJnZb65h+669+72nU2Zak2MUPQuzKc7VgHZzYpoI3H0PUGFmBcCfzWyIu79+CJ8zBZgCEIlE1FsQ6aBSU4zu+dGxg0gj+92dtbU7G72XYtX6bby8fB2bGwxoZ6Wn0KPgo15EWYNpsl06aUA7LrOG3H2Dmc0GxgD1E8FqoBewyszSgHyig8YiIgcws31LeQzvVdBom43bdtVLFB8lidUbtrHovQ9YW3vggHb3gnoD2vvdeJdDt/ysDj+gHctZQ6XAriAJZAOnEh0Mrm8WMAF4ETgHeEbjAyJyOPKz08nPTmdwj8YHtLfu3M17G7Y1uPEuuv2ft2r4cPOBA9pdO2U1eS9Fz4LshB/QjmWPoDswIxgnSAEecvfHzOx6oNLdZwFTgZlmthRYB4yLYTwiIuRkpDGgSycGdOnU6P6du/fy/sYgQeyXKLby6jvreXzB++xusJZHcW7G/omiIJuehR+NU3TOat8D2pZof4BHIhGvrKwMOwwRSVJ79jofbtre4Ma7rfv1MHbsbjCgnZXW5BhFz4JsiuIwoG1mVe7e2NCL7iwWEWmN1BSjR0E2PQqyGV1+4H53p2bLzkbvpVi1fisvL1vL5gYV7LLTU+lRkPVRL6LedNmywhy6dMokJYYD2koEIiJtyMwo7ZRJaadMKg4yoF13R3bDcYrXV29kXYMB7fTU6Gyqb3/mCMZW9GzzmJUIRETiLDqgnc/RPRqvTbF15+79xijqZj6V5GXGJB4lAhGRdiYnI42BXTsxsGvjA9ptrWNPjhURkWYpEYiIJDklAhGRJKdEICKS5JQIRESSnBKBiEiSUyIQEUlySgQiIkku4RadM7NqYOUhvr2EBtXP2gnF1TqKq/Xaa2yKq3UOJ64+7l7a2I6ESwSHw8wqm1p9L0yKq3UUV+u119gUV+vEKi5dGhIRSXJKBCIiSS7ZEsGUsANoguJqHcXVeu01NsXVOjGJK6nGCERE5EDJ1iMQEZEGlAhERJJch0gEZjbNzNaY2etN7Dczu9XMlprZAjMbWW/fBDN7K3hMiHNcFwXxvGZmL5jZ8Hr7VgSvzzOzyjjH9Skz2xh89zwzu6bevjFm9mZwLL8f57i+Uy+m181sj5kVBftiebx6mdlsM1tkZgvN7BuNtIn7OdbCuOJ+jrUwrrifYy2MK+7nmJllmdkrZjY/iOu6RtpkmtmDwTF52czK6+27Onj9TTM77ZCCcPeEfwCfBEYCrzex/3TgH4ABxwIvB68XAcuC/y0MtgvjGNdxdd8H/FddXMHzFUBJSMfrU8BjjbyeCrwN9AMygPnA4HjF1aDt54Bn4nS8ugMjg+1OwJKG/+4wzrEWxhX3c6yFccX9HGtJXGGcY8E5kxdspwMvA8c2aPNV4M5gexzwYLA9ODhGmUDf4NiltjaGDtEjcPd/A+sO0mQscI9HvQQUmFl34DTgKXdf5+7rgaeAMfGKy91fCL4X4CWgrK2++3DiOohjgKXuvszddwIPED22YcR1AXB/W333wbj7++7+arC9GXgDaFhBPO7nWEviCuMca+HxakrMzrFDiCsu51hwzmwJnqYHj4azeMYCM4LtR4BTzMyC1x9w9x3uvhxYSvQYtkqHSAQt0BN4t97zVcFrTb0ehklE/6Ks48CTZlZlZpNDiOfjQVf1H2Z2dPBauzheZpZD9Mf0T/VejsvxCrrkI4j+1VZfqOfYQeKqL+7nWDNxhXaONXe84n2OmVmqmc0D1hD9w6HJ88vddwMbgWLa6HipeH07YGYnEf2P9IR6L5/g7qvNrAvwlJktDv5ijodXia5LssXMTgf+AgyM03e3xOeA5929fu8h5sfLzPKI/jBc5e6b2vKzD0dL4grjHGsmrtDOsRb+/xjXc8zd9wAVZlYA/NnMhrh7o2NlsZAsPYLVQK96z8uC15p6PW7MbBhwNzDW3dfWve7uq4P/XQP8mUPo7h0qd99U11V1978D6WZWQjs4XoFxNOiyx/p4mVk60R+P+9z90UaahHKOtSCuUM6x5uIK6xxryfEKxP0cCz57AzCbAy8f7jsuZpYG5ANraavj1dYDH2E9gHKaHvz8LPsP5L0SvF4ELCc6iFcYbBfFMa7eRK/pHdfg9VygU73tF4AxcYyrGx/dbHgM8E5w7NKIDnb25aOBvKPjFVewP5/oOEJuvI5X8G+/B7jlIG3ifo61MK64n2MtjCvu51hL4grjHANKgYJgOxt4DjijQZsr2H+w+KFg+2j2HyxexiEMFneIS0Nmdj/RWQglZrYK+AnRARfc/U7g70RndSwFtgJfCvatM7OfAnOCj7re9+8Kxjqua4he57s9Ou7Dbo+uLNiVaPcQov9h/NHdn4hjXOcAl5vZbmAbMM6jZ91uM/sa8E+iszumufvCOMYF8HngSXevrffWmB4v4HjgYuC14DouwA+I/siGeY61JK4wzrGWxBXGOdaSuCD+51h3YIaZpRK9SvOQuz9mZtcDle4+C5gKzDSzpUST1Lgg5oVm9hCwCNgNXOHRy0ytoiUmRESSXLKMEYiISBOUCEREkpwSgYhIklMiEBFJckoEIiJJTolAkoqZ/cLMTjKzs8zs6ibalAYrPM41s08cwndMNLMehx+tSHwoEUiy+RjRxddOBJpaHuAU4DV3H+Huzx3Cd0wEWpUIgrtFRUKhRCBJwcxuNLMFwGjgReDLwB1Wbx38oF0F8CtgbLDufLaZfcbMXjSzV83s4WCtGszsGjObY9F166dY1DlABLiv3vtXBMsnYGYRM3s22L7WzGaa2fNEbxYqNbM/BZ85x8yOD9qdaB+tkT/XzDrF5aBJ8miLW6T10CMRHkSTwO+I3q38/EHaTQRuC7ZLiPYccoPn3wOuCbaL6r1nJvC5YPtZIFJv3wqCdeyJJolng+1rgSogO3j+R6ILm0H0btc3gu2/AccH23lAWtjHUo+O9VB3VJLJSKLrsgwiuhZ9SxxLtPjH88HyAhlEexQAJ5nZd4EcomsKLST6o90as9x9W7D9aWBw8D0AnYPex/PATWZ2H/Cou69q5XeIHJQSgXR4weWe6URXZqwh+sNtwXozH6/3Q9zo24muD39Bg8/MAm4n+pf/u2Z2LZDVxGfs5qPLsA3b1F/PJoVoZartDdr80sweJ7qW0fNmdpq7Lz5IzCKtojEC6fDcfZ67VxCUJgSeAU5z94pmkgBEB5aPN7MBAGaWa2ZH8NEPek3wV/s59d6zmWgpxDorgFHB9hcO8l1PAlfWPQkSGGbW391fc/cbiC5eN6iZmEVaRYlAkoKZlQLr3X0vMMjdF7Xkfe5eTXTM4P5gsPnF4P0bgD8ArxNdKXNOvbdNB+6sGywGrgN+a9GC5wdbGfLrQMSixeYXAV8JXr8qGJBeAOxi/ypjIodNq4+KiCQ59QhERJKcEoGISJJTIhARSXJKBCIiSU6JQEQkySkRiIgkOSUCEZEk9/8K+0V/KPM8fQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot to show an example to plot the RSS of best models with different number of parameters\n",
    "plt.figure()\n",
    "plt.plot(models[\"RSS\"])\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('RSS')\n",
    "plt.show()\n",
    "# RSS is falling with features- 3 variable model is the best one on the basis of RSS- we can keep adding features further to find best model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1    0.663357\n",
       "2    0.759667\n",
       "3    0.768905\n",
       "dtype: float64"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rsquared_adj = models.apply(lambda a: \n",
    "                            a[1].rsquared_adj, axis=1)\n",
    "rsquared_adj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnyElEQVR4nO3deXxU9b3/8deHhCRAIGGVzQAKimBFYFgUrVur/NqK19/FCi4FhYtoq8XfT6+t3l+1tt57axetrZUqWgEFVNSW1talbrelLkkA2SlIUSJbMBI2Q0jy+f0xhzKEISSQmTPJvJ+PxzyYOcvkneHkfOZ8v+d8j7k7IiIitbUIO4CIiKQmFQgREYlLBUJEROJSgRARkbhUIEREJK7MsAM0lk6dOnnv3r3DjiEi0qQUFxdvd/fO8eY1mwLRu3dvioqKwo4hItKkmNlHR5qnJiYREYlLBUJEROJSgRARkbhUIEREJC4VCBERiUsFQkRE4lKBEBGRuJrNdRAiIumgqrqGbbv2sbn8czbtqGBz+efkZrfkqhEFjf6zVCBERFJETY3z6Z7KQ3b+m8sr+GTH52zeEX2+bdc+qmsOvY/P4IJ8FQgRkabK3dn5eRWbyj8/pABs2lHBpmDnv6W8gsrqmkPWy85sQff8VnTLy+HskzvRPT+Hbnmt6JafQ/fg33Y5LROSWQVCRKQR7K2sOvitf0fwrT84AjhQAPZWVh+yTmYL44R2OXTPz+HME/Pp9oVgp5+X88+i0KFNFmYWyu+kAiEichT7qqrZWr7vkG//B3b6B/4t/3z/IeuYQefcbLrlt+KUE9py3ildDvn23yO/FZ1ys8loEc7Ovz5UIEQkrVXXONt2VcTs9A9t/9+0o4Ltu/cdtl771i3plteKnu1bMax3h382+Rz45n9CuxyyMpv2iaIqECLSbLlHO3037ai904/+u3nH52yN0+mbm51Jt7wcuuW3YkC3dnTLa0X3/IPNPt3yWtEqKyOk3yp5VCBEpElyd3ZWVB3+rX9HRdAUVMHm8goqqw7t9M3KbEH3YCc/8uSO9MhvlbRO36ZGBUJEUlLtTt9NtXf+Oz5nT61O34wWRtd2OXTLy2FQz3xGn55anb5NjQqEiCRdZVUNW8orDjvlM1oAos937D2807dTbjbd81vRr0suX+zX+bBTPju3Te1O36ZGBUJEGlVsp2/tb/+byz9nU3kFpbsO7/TNb92S7nmt6JGfQ6RX+4NNPsG3/+bQ6dvUqECISL0d6PQ9uNMPOn2DJp/N5RVs2VlxWKdvm6yMaBNPfitOCzp9Y9v8u6dJp29TowIhIsDBTt+4bf4xbf91dfqOOKnDITv9bkETULucTLX7N0EqECJp4vPK6sN2+geafDYFY/3U1el7Rs98Rg/M+efpnwcKQEd1+jZbKhAizUBlVQ1bd8Zc3XugEOz4/IidvgCd22bTPS+Hvp1zObdfp0O+9ffIV6dvuktogTCz0cDPgQxghrv/d635DwAXBC9bA13cPT+YVwDMAE4EHPiKu29IZF6RVFRd45Tu2nfot/7Yq37Lo1f6+qHN/uQHV/p2z8thaK/8gxd7BVf7qtNXjiZhBcLMMoCHgS8DJUChmS1w95UHlnH3W2OWvxkYHPMWs4D73P01M8sFDm34FGkG3J2yPZWHXN37z2/9weutOyuoitPp2y04r79/13bRNv+YZp9ueTm0zlIDgRyfRG5Bw4F17r4ewMzmAZcBK4+w/Hjg7mDZAUCmu78G4O67E5hTJKFKd+1j+abymCafg0cCm8sr2Ben0zc6nEMOI/p0CM7+UaevJF8iC0QPYGPM6xJgRLwFzawX0Ad4I5h0CrDDzF4Ipv8Z+I67V9dabwowBaCgoPFvliFyvNZt28Xlv/obuyqqgGin7wltoyN8nt4jj0sGdlWnr6SsVDkGHQfMjykAmcC5RJucPgaeASYCj8eu5O6PAo8CRCKRWi2wIuEq21PJ9U8WkZ2Zwa//bSh9OrWhS9scdfpKk5HIHqpPiHYwH9AzmBbPOGBuzOsSYIm7r3f3KuC3wJBEhBRJhH1V1UydXczWnRU89o2hnH1yJ7rltVJxkCYlkQWiEOhnZn3MLItoEVhQeyEz6w+0B96ptW6+mXUOXl/IkfsuRFKKu3PnC8t5f0MZP/36IAYXtA87ksgxSViBCL75fwt4BVgFPOvuK8zsXjMbE7PoOGCe+8GT9IKmptuA181sGWDAY4nKKtKYHnn7Q55fVML/+fIpfO2M7mHHETlm5rVPnm6iIpGIFxUVhR1D0tzLyzcz9alFXHZmdx688kx1NkvKM7Nid4/Em6erZEQaybKScqY9s4QhBfn86F/PUHGQJk8FQqQRbCmvYPKsQjrlZvPoNyLktNTIpNL0pcppriJN1t7KKibNLGTPvmqev3EEnXKzw44k0ih0BCFyHGpqnGnzlrBq805+cdVgTu3aNuxIIo1GBULkONz/yhpeXbmV731tABec2iXsOCKNSgVC5Bg9W7SR6W9/yLUjezHh7N5hxxFpdCoQIsfg3fWfcteLyzi3XyfuvnSAzliSZkkFQqSBNmzfw9Sniino0JpfXjWEzAz9GUnzpC1bpAHK9+7n+pmFGPDExGHktWoZdiSRhNFpriL1tL+6hpvmFLOxbC9PTx5Jr45two4kklAqECL14O5873crWLjuU35yxSCG9+kQdiSRhFMTk0g9PLFwA3Pf/5ibzj+ZsUN7hh1HJClUIESO4vVVW/nhSysZPbArt118athxRJJGBUKkDqs27+SWuYs5vXseP7tyEC10wx9JIyoQIkewbVcFk2cW0TanJTMmRGidpS47SS/a4kXiqNhfzZRZxZTtqeS5qWdxQrucsCOJJJ0KhEgt7s5tz33AByU7mH7NUE7vkRd2JJFQqIlJpJYH/7yWPyzdzB2j+3PJwK5hxxEJjQqESIzfLfmEn7++liuG9uSGL54UdhyRUKlAiASKP/qM2+cvZXifDtx3+Rc0AJ+kPRUIEWBj2V5umF1Et7wcfn3NULIy9achok5qSXu7KvYzeWYRlVU1zJsyjPZtssKOJJISVCAkrVVV13Dz3MWsK93NrOuH07dLbtiRRFKGjqMlrd33x1W8taaUH1x2OqP6dgo7jkhKUYGQtDX73Y/4zcINTDqnD1eNKAg7jkjKUYGQtPSXtaXcs2AFF/bvwp1fOS3sOCIpSQVC0s66bbu46elF9OuSy0PjB5OhAfhE4lKBkLRStqeS658sIjszgxkTIuRm6zwNkSNJaIEws9FmtsbM1pnZd+LMf8DMlgSPv5vZjlrz25lZiZn9MpE5JT3sq6pm6uxituys4LFvDKVn+9ZhRxJJaQn7+mRmGcDDwJeBEqDQzBa4+8oDy7j7rTHL3wwMrvU2PwD+J1EZJX24O3e+sJz3N5Txi/GDGVzQPuxIIikvkUcQw4F17r7e3SuBecBldSw/Hph74IWZDQVOAF5NYEZJE4+8/SHPLyrh1i+dwqWDuocdR6RJSGSB6AFsjHldEkw7jJn1AvoAbwSvWwA/BW6r6weY2RQzKzKzotLS0kYJLc3Py8s3c//LaxgzqDu3XNQ37DgiTUaqdFKPA+a7e3Xw+ibgj+5eUtdK7v6ou0fcPdK5c+eEh5SmZ1lJOdOeWcKQgnzuH3uGBuATaYBEnsLxCXBizOuewbR4xgHfjHl9FnCumd0E5AJZZrbb3Q/r6BY5ki3lFUyeVUjHNtn8+toIOS0zwo4k0qQkskAUAv3MrA/RwjAOuKr2QmbWH2gPvHNgmrtfHTN/IhBRcZCG2FtZxaSZhezZV83zN46gc9vssCOJNDkJa2Jy9yrgW8ArwCrgWXdfYWb3mtmYmEXHAfPc3ROVRdJLTY0zbd4SVm3eyS/GD+bUrm3DjiTSJFlz2S9HIhEvKioKO4akgP/+02qmv/0hd186gOtG9Qk7jkhKM7Nid4/Em5cqndQijeLZoo1Mf/tDrhlZwMSze4cdR6RJU4GQZuPd9Z9y14vLOLdfJ+6+dKDOWBI5TioQ0ixs2L6HqU8VU9ChNb+8aggtM7Rpixwv/RVJk1e+dz/XzyzEgCcmDiOvVcuwI4k0CxrKUpq0/dU13DSnmI1le3l68kh6dWwTdiSRZkMFQposd+d7v1vBwnWf8pMrBjG8T4ewI4k0K2pikibriYUbmPv+x9x0/smMHdoz7DgizY4KhDRJr6/ayg9fWsnogV257eJTw44j0iypQEiTs2rzTm6Zu5jTu+fxsysH0UK3DBVJCBUIaVK27apg8swi2ua0ZMaECK2z1I0mkij665Imo2J/NVNmFVO2p5Lnpp7FCe1ywo4k0qypQEiT4O7c9twHfFCyg+nXDOX0HnlhRxJp9tTEJE3Cg39eyx+WbuaO0f25ZGDXsOOIpAUVCEl5v1vyCT9/fS1XDO3JDV88Kew4ImlDBUJSWvFHn3H7/KUM79OB+y7/ggbgE0kiFQhJWRvL9nLD7CK65eXw62uGkpWpzVUkmdRJLSlpV8V+Js8sorKqhnlThtG+TVbYkUTSjgqEpJyq6hpunruYdaW7mXX9cPp2yQ07kkha0jG7pJz7/riKt9aU8oPLTmdU305hxxFJWyoQklJmv/sRv1m4gUnn9OGqEQVhxxFJayoQkjL+sraUexas4ML+XbjzK6eFHUck7R2xQJhZOzP7LzObbWZX1Zr3q8RHk3Sybtsubnp6Ef265PLQ+MFkaAA+kdDVdQTxG8CA54FxZva8mWUH80YmPJmkjbI9lVz/ZBHZmRnMmBAhN1vnToikgroKxMnu/h13/627jwEWAW+YWcckZZM0sK+qmqmzi9mys4LHvjGUnu1bhx1JRAJ1fVXLNrMW7l4D4O73mdknwP8AOu9Qjpu7c+cLy3l/Qxm/GD+YwQXtw44kIjHqOoL4PXBh7AR3fxL4v0BlAjNJmnjk7Q95flEJt37pFC4d1D3sOCJSyxGPINz9348w/WWgX8ISSVp4eflm7n95DWMGdeeWi/qGHUdE4jjqaa5m1jYZQSR9LCspZ9ozSxhSkM/9Y8/QAHwiKarOAmFmPYA/JCmLpIEt5RVMnlVIxzbZ/PraCDktM8KOJCJHUNd1EAOBl4E7jvXNzWy0ma0xs3Vm9p048x8wsyXB4+9mtiOYfqaZvWNmK8xsqZldeawZJHXsraxi0sxC9uyr5omJw+jcNvvoK4lIaOo6i+lN4F/c/d1jeWMzywAeBr4MlACFZrbA3VceWMbdb41Z/mZgcPByL/ANd19rZt2BYjN7xd13HEsWCV9NjTNt3hJWbd7J4xOGcWpXtVyKpLq6mpgKgcuP472HA+vcfb27VwLzgMvqWH48MBfA3f/u7muD55uAbUDn48giIbv/lTW8unIr/+9rA7igf5ew44hIPdRVIMYAeWZ2/zG+dw9gY8zrkmDaYcysF9AHeCPOvOFAFvBhnHlTzKzIzIpKS0uPMaYk2rNFG5n+9odcM7KAiWf3DjuOiNTTEQuEu1e7+xRgdxJyjAPmu3t17EQz6wbMBq47cMFerYyPunvE3SOdO+sAIxW9u/5T7npxGef268Tdlw7UGUsiTchRT3N193uP8b0/AU6Med0zmBbPOILmpQPMrB3wEnDXsfaDSLg2bN/D1KeKKejQml9eNYSWGRo8WKQpafBfrJnlm9ld9Vi0EOhnZn3MLItoEVgQ5/36A+2Bd2KmZQEvArPcfX5DM0r4yvfu5/qZhRjwxMRh5LVqGXYkEWmguk5zPdHMHjWzP5jZZDNrY2Y/Bf4OHLWX0d2rgG8BrwCrgGfdfYWZ3WtmY2IWHQfMc3ePmfZ14IvAxJjTYM9s+K8nYdhfXcM35yxiY9lefn1thF4d24QdSUSOgR26X46ZYfYm8DbRb/ajg8cS4FZ335KsgPUViUS8qKgo7Bhpz92567fLmfPex/zkikGMHdoz7EgiUgczK3b3SLx5dV0H0cHd7wmev2JmVwBXx+ssFjngNws3MOe9j7np/JNVHESauDrvzGJm7YneNAjgU6KnvRqAu5clOJs0MW+s3soPX1rJ6IFdue3iU8OOIyLHqa4CkQcUc7BAQPSmQQAOnJSoUNL0rN6yk5vnLGZg9zx+duUgWuiWoSJNXl3DffdOYg5pwkp37WPSk0W0zWnJjAkRWmfplqEizYH+kuW4VOyvZsrsIsr2VPLc1LM4oV1O2JFEpJGoQMgxc3dun7+UJRt3MP2aoZzeIy/sSCLSiHRpqxyzn7++lt9/sIk7RvfnkoFdw44jIo3siEcQZtbO3XeaWYc4sx3YWXvsJEkfv1vyCQ/+eS1XDO3JDV/U+QoizVFdTUxzgK8RPZPJOfRsJoBcM3vM3e9MVDhJTYs+/ozb5y9leJ8O3Hf5FzQAn0gzVddZTF8L/u0Tb35wQ6DlgApEGin5bC9TZhXRLS+HX18zlKxMtVKKNFd1NTENqWtFd18EnNboiSRl7arYz6Qni6isqmHelGG0b5MVdiQRSaC6mph+GvybA0SAD4g2M50BFAFnJTaapJLqGueWuYtZV7qbWdcPp2+X3LAjiUiC1XXDoAvc/QJgMzAkuDHPUKL3jT7SfR2kmbrvpVW8uaaUH1x2OqP6dgo7jogkQX0akE9192UHXrj7ctS0lFaeevcjnlj4D64f1YerRhSEHUdEkqQ+F8otNbMZwFPB66uBpYmLJKnkr2u3c/eCFVzYvwt3fVXfC0TSSX0KxHXAjcC3g9dvA48kLJGkjHXbdnPj08X065LLQ+MHk6EB+ETSSn3uSV3h7g+4++XufjnR6yJ+lvhoEqayPZVMmllIdmYLZkyIkJutUVlE0k29/urNbDAwnuitQP8BvJDIUBKufVXVTJ1dzObyCuZNGUnP9q3DjiQiIajrOohTiBaF8cB24Bmityi9IEnZJATuzp0vLOf9DWU8NH4wQwrahx1JREJS1xHEauAvwNfcfR2Amd2alFQSmulvr+f5RSVM+1I/xgzqHnYcEQlRXX0Q/5voNRBvmtljZnYRh4/HJM3Iy8s386OXVzNmUHe+fVG/sOOISMjqulDut+4+DugPvAlMA7qY2SNmdnGS8kmSLCspZ9ozSxhckM/9Y8/QAHwiUq+zmPa4+xx3vxToCSwG7kh4MkmaLeUVTJ5VSMc22Tx6bYSclhlhRxKRFNCgoTjd/TN3f9TdL0pUIEmuvZVVTJ5VyO6KKh6fGKFz2+ywI4lIitDJ7Wmspsa59ZklrNy0kxkTIvTv2i7sSCKSQjSYfxr78atreGXFVv7jqwO4sP8JYccRkRSjApGmnivayCNvfcjVIwq4blTvsOOISApSgUhD763/lDtfXMY5fTtxz5iBOmNJROJKaIEws9FmtsbM1pnZd+LMf8DMlgSPv5vZjph5E8xsbfCYkMic6WTD9j3c8FQxBR1a8/DVQ2iZoe8IIhJfwjqpg3tWPwx8GSgBCs1sgbuvPLCMu98as/zNRG9GhJl1AO4meic7B4qDdT9LVN50UL53P9fPLMSAJyYOI69Vy7AjiUgKS+TXx+HAOndf7+6VwDzgsjqWHw/MDZ5fArzm7mVBUXgNGJ3ArM3e/uoavjlnERvL9jL9mqH06tgm7EgikuISWSB6ABtjXpcE0w5jZr2APsAbDVnXzKaYWZGZFZWWljZK6ObI3bl7wQr+um47/3n5FxhxUsewI4lIE5AqDdDjgPnuXt2QlYKL9iLuHuncuXOCojV9v1m4gTnvfcyN55/MFZETw44jIk1EIgvEJ0Ds3qhnMC2ecRxsXmroulKHN1Zv5YcvreSSgSdw+8Wnhh1HRJqQRBaIQqCfmfUxsyyiRWBB7YXMrD/QHngnZvIrwMVm1t7M2gMXB9OkAVZv2cnNcxYzoHs7HrjyTFrolqEi0gAJO4vJ3avM7FtEd+wZwBPuvsLM7gWK3P1AsRgHzHN3j1m3zMx+QLTIANzr7mWJytocle7ax6Qni8jNyWTGN4bROkujqohIw1jMfrlJi0QiXlRUFHaMlFCxv5rxj73L6s27eG7qWZzeIy/sSCKSosys2N0j8ebpa2Uz4+7cPn8piz/ewfRrhqo4iMgxS5WzmKSR/Pz1tfz+g03cMbo/o0/vGnYcEWnCVCCakd8t+YQH/7yWsUN7MvW8k8KOIyJNnApEM7Ho48+4ff5ShvfpwH9e/gUNwCcix00Fohko+WwvU2YV0S0vh+nXDCUrU/+tInL81EndxO2q2M+kJ4vYV1XDvCnD6NAmK+xIItJMqEA0YdU1zi1zF7OudDczrxtO3y65YUcSkWZEbRFN2H0vreLNNaV8f8xAzunXKew4ItLMqEA0UU+9+xFPLPwH14/qwzUje4UdR0SaIRWIJuiva7dz94IVXNi/C3d99bSw44hIM6UC0cSs27abG58upl+XXB4aP5gMDcAnIgmiAtGElO2pZNLMQrIzWzBjQoTcbJ1jICKJoz1ME7Gvqpqps4vZXF7BvCkj6dm+ddiRRKSZ0xFEE+Du3PnCct7fUMZPrhjEkIL2YUcSkTSgAtEETH97Pc8vKmHal/oxZlD3sOOISJpQgUhxLy/fzI9eXs2YQd359kX9wo4jImlEBSKFLSspZ9ozSxhckM/9Y8/QAHwiklQqEClqS3kFk2cV0rFNNo9eGyGnZUbYkUQkzahApKC9lVVMnlXI7ooqHp8YoXPb7LAjiUga0mmuKaamxrn1mSWs3LSTGRMi9O/aLuxIIpKmdASRYn786hpeWbGV//jqAC7sf0LYcUQkjalApJDnijbyyFsfcvWIAq4b1TvsOCKS5lQgUsR76z/lzheXcU7fTtwzZqDOWBKR0KlApIAN2/dww1PFFHRozcNXD6Flhv5bRCR82hOFrHzvfq6fWYgBT0wcRl6rlmFHEhEBdBZTqPZX1/DNOYvYWLaXpyaNoFfHNmFHEhH5JxWIkLg7dy9YwV/XbefHY89gxEkdw44kInIINTGF5DcLNzDnvY+58fyTuSJyYthxREQOowIRgjdWb+WHL63kkoEncPvFp4YdR0QkroQWCDMbbWZrzGydmX3nCMt83cxWmtkKM5sTM/3+YNoqM3vImsl5n6u37OTmOYsZ0L0dD1x5Ji10y1ARSVEJ64MwswzgYeDLQAlQaGYL3H1lzDL9gO8Co9z9MzPrEkw/GxgFnBEs+lfgPOCtROVNhtJd+5j0ZBG5OZnM+MYwWmepC0hEUlcijyCGA+vcfb27VwLzgMtqLfNvwMPu/hmAu28LpjuQA2QB2UBLYGsCsyZcxf5qpswuomxPJY9PGEbXvJywI4mI1CmRBaIHsDHmdUkwLdYpwClmttDM3jWz0QDu/g7wJrA5eLzi7qtq/wAzm2JmRWZWVFpampBfojG4O7fPX8rij3fwwJVncnqPvLAjiYgcVdid1JlAP+B8YDzwmJnlm1lf4DSgJ9GicqGZnVt7ZXd/1N0j7h7p3LlzEmM3zM9fX8vvP9jEHaP7M/r0rmHHERGpl0QWiE+A2PM3ewbTYpUAC9x9v7v/A/g70YJxOfCuu+92993An4CzEpg1YX635BMe/PNaxg7tydTzTgo7johIvSWyQBQC/cysj5llAeOABbWW+S3RowfMrBPRJqf1wMfAeWaWaWYtiXZQH9bElOoWffwZt89fyvA+HfjPy7+gAfhEpElJWIFw9yrgW8ArRHfuz7r7CjO718zGBIu9AnxqZiuJ9jnc7u6fAvOBD4FlwAfAB+7++0RlTYSSz/YyZVYR3fJymH7NULIyw27NExFpGHP3sDM0ikgk4kVFRWHHAGBXxX7GPvIOm8o/58WbRtG3S27YkURE4jKzYnePxJunE/EbWXWNc8vcxawr3c3M64arOIhIk6V2j0Z230ureHNNKd8fM5Bz+nUKO46IyDFTgWhET737EU8s/AfXj+rDNSN7hR1HROS4qEA0kr+u3c7dC1ZwYf8u3PXV08KOIyJy3FQgGsG6bbu58eli+nXJ5aHxg8nQAHwi0gyoQBynsj2VTJpZSHZmC2ZMiJCbrX5/EWketDc7Dvuqqpk6u5jN5RXMmzKSnu1bhx1JRKTR6AjiGLk7d76wnPc3lPGTKwYxpKB92JFERBqVCsQxmv72ep5fVMK0L/VjzKDuYccREWl0KhDH4OXlm/nRy6sZM6g7376oX9hxREQSQgWigZaVlDPtmSUMLsjn/rFnaAA+EWm2VCAaYEt5BZNnFdKxTTaPXhshp2VG2JFERBJGBaKe9lZWMXlWIbsrqnh8YoTObbPDjiQiklA6zbUeamqcW59ZwspNO5kxIUL/ru3CjiQiknA6gqiHH7+6hldWbOU/vjqAC/ufEHYcEZGkUIE4iueKNvLIWx9y9YgCrhvVO+w4IiJJowJRh/fWf8qdLy7jnL6duGfMQJ2xJCJpRQXiCDZs38MNTxVT0KE1D189hJYZ+qhEJL1orxdH+d79XD+zEIDHJwwjr1XLkBOJiCSfCkQt+6tr+OacRWws28v0a4bSu1ObsCOJiIRCp7nGcHfuXrCCv67bzv1jz2DkSR3DjiQiEhodQcT4zcINzHnvY6aedzJfj5wYdhwRkVCpQATeWL2VH760kosHnMC/X3Jq2HFEREKnAgGs3rKTm+cs5rRu7Xhw3Jm00C1DRURUIEp37WPSk0Xk5mTy+IRhtM5St4yICKiTmpYZRv+ubZn2pVPompcTdhwRkZSR9gUiv3UWj08cFnYMEZGUk/ZNTCIiEl9CC4SZjTazNWa2zsy+c4Rlvm5mK81shZnNiZleYGavmtmqYH7vRGYVEZFDJayJycwygIeBLwMlQKGZLXD3lTHL9AO+C4xy98/MrEvMW8wC7nP318wsF6hJVFYRETlcIo8ghgPr3H29u1cC84DLai3zb8DD7v4ZgLtvAzCzAUCmu78WTN/t7nsTmFVERGpJZIHoAWyMeV0STIt1CnCKmS00s3fNbHTM9B1m9oKZLTazHwdHJIcwsylmVmRmRaWlpQn5JURE0lXYndSZQD/gfGA88JiZ5QfTzwVuA4YBJwETa6/s7o+6e8TdI507d05SZBGR9JDIAvEJEDugUc9gWqwSYIG773f3fwB/J1owSoAlQfNUFfBbYEgCs4qISC2JLBCFQD8z62NmWcA4YEGtZX5L9OgBM+tEtGlpfbBuvpkdOCy4EFiJiIgkTcLOYnL3KjP7FvAKkAE84e4rzOxeoMjdFwTzLjazlUA1cLu7fwpgZrcBr1v0Pp/FwGN1/bzi4uLtZvbRcUTuBGw/jvUTRbkaRrkaRrkapjnm6nWkGebux/iezYuZFbl7JOwctSlXwyhXwyhXw6RbrrA7qUVEJEWpQIiISFwqEAc9GnaAI1CuhlGuhlGuhkmrXOqDEBGRuHQEISIicalAiIhIXM2+QJjZE2a2zcyWH2G+mdlDwZDkS81sSMy8CWa2NnhMSHKuq4M8y8zsb2Y2KGbehmD6EjMrSnKu882sPPjZS8zsezHzjjq8ewJz3R6TabmZVZtZh2BeIj+vE83szZgh678dZ5mkbmP1zBTW9lWfbEnfxuqZK+nbmJnlmNn7ZvZBkOv7cZbJNrNngs/kPYu5NYKZfTeYvsbMLmlwAHdv1g/gi0SH6Vh+hPlfAf4EGDASeC+Y3oHoVd0dgPbB8/ZJzHX2gZ8H/K8DuYLXG4BOIX1e5wN/iDM9A/iQ6LhZWcAHwIBk5aq17KXAG0n6vLoBQ4LnbYkOFzOg1jJJ3cbqmSms7as+2ZK+jdUnVxjbWLDN5AbPWwLvASNrLXMTMD14Pg54Jng+IPiMsoE+wWeX0ZCf3+yPINz9f4CyOha5DJjlUe8SHeKjG3AJ8Jq7l3l0OPLXgNF1vE+j5nL3vwU/F+BdomNZJVw9Pq8jqc/w7snKNR6Y21g/uy7uvtndFwXPdwGrOHzU4qRuY/XJFOL2VZ/P60gSto0dQ66kbGPBNrM7eNkyeNQ+s+gyYGbwfD5wkZlZMH2eu+/z6Fh364h+hvXW7AtEPRxpWPL6DFeeLJOIfgM9wIFXzazYzKaEkOes4JD3T2Y2MJiWEp+XmbUmupN9PmZyUj6v4NB+MNFvebFC28bqyBQrlO3rKNlC28aO9pklexszswwzWwJsI/qF4ojbl0cHNy0HOtIIn1fCxmKSxmFmFxD9Az4nZvI57v6JRe/A95qZrQ6+YSfDIqCXu+82s68QHXCxX5J+dn1cCix099ijjYR/Xha96+HzwDR339mY732s6pMprO3rKNlC28bq+f+Y1G3M3auBMy16K4QXzex0d4/bF9fYdARx5GHJ6zNceUKZ2RnADOAyDwYxBHD3T4J/twEv0sDDxuPh7jsPHPK6+x+BlhYdiTf0zyswjlqH/on+vMysJdGdytPu/kKcRZK+jdUjU2jb19GyhbWN1eczCyR9GwveewfwJoc3Q/7zczGzTCAP+JTG+Lwau1MlFR9Ab47c6fpVDu1AfD+Y3gH4B9HOw/bB8w5JzFVAtM3w7FrT2wBtY57/DRidxFxdOXiB5XDg4+CzyyTaydqHgx2IA5OVK5ifR7Sfok2yPq/gd58FPFjHMkndxuqZKZTtq57Zkr6N1SdXGNsY0BnID563Av4CfK3WMt/k0E7qZ4PnAzm0k3o9DeykbvZNTGY2l+hZEZ3MrAS4m2hHD+4+Hfgj0bNM1gF7geuCeWVm9gOi96YAuNcPPaRMdK7vEW1H/FW0v4kqj47WeALRw0yI/sHMcfeXk5hrLHCjmVUBnwPjPLo1xh3ePYm5AC4HXnX3PTGrJvTzAkYB1wLLgnZigDuJ7oDD2sbqkymU7aue2cLYxuqTC5K/jXUDZlr0lsstiO78/2CH3jbhcWC2ma0jWrzGBZlXmNmzRO+lUwV806PNVfWmoTZERCQu9UGIiEhcKhAiIhKXCoSIiMSlAiEiInGpQIiISFwqECKAmf2XmV1gZv9iZt89wjKdg9EyF5vZucfwMyaaWffjTyuSHCoQIlEjiA5adx5wpCESLgKWuftgd//LMfyMiUCDCkRwZaxIKFQgJK2Z2Y/NbCkwDHgHmAw8YjH3IAiWOxO4H7gsGPO/lZldbGbvmNkiM3suGMcHM/uemRVa9J4Bj1rUWCACPB2z/oZgCAnMLGJmbwXP7zGz2Wa2kOgFUJ3N7PngPQvNbFSw3Hl28P4Ei82sbVI+NEkfjXUJvR56NNUH0eLwC6JXZi+sY7mJwC+D552IHmm0CV7fAXwveN4hZp3ZwKXB87eASMy8DQT3ECBaPN4Knt8DFAOtgtdziA4GB9Ere1cFz38PjAqe5wKZYX+WejSvhw5fRaI3IvoA6E/0PgD1MZLoDVkWBkMsZBE9AgG4wMz+HWhNdLylFUR35g2xwN0/D55/CRgQ/ByAdsHRykLgZ2b2NPCCu5c08GeI1EkFQtJW0Gz0JNFRLrcT3aFbMBbPWTE76LirEx2bf3yt98wBfkX0SGGjmd0D5BzhPao42Mxbe5nYsX5aEL2LWEWtZf7bzF4iOs7TQjO7xN1X15FZpEHUByFpy92XuPuZBLeXBN4ALnH3M49SHCDaoT3KzPoCmFkbMzuFgzv67cG3/LEx6+wiejvLAzYAQ4Pn/1rHz3oVuPnAi6CwYWYnu/syd/8R0QH/+h8ls0iDqEBIWjOzzsBn7l4D9Hf3lfVZz91LifZJzA06ud8J1t8BPAYsJzrqaGHMak8C0w90UgPfB35u0Zvc1zXK5i1AxMyWmtlKYGowfVrQEb4U2M+hd4UTOW4azVVEROLSEYSIiMSlAiEiInGpQIiISFwqECIiEpcKhIiIxKUCISIicalAiIhIXP8fwM2fQSYwU7oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot to show an example to plot the RSS of best models with different number of parameters\n",
    "plt.figure()\n",
    "plt.plot(rsquared_adj)\n",
    "plt.xlabel('# features')\n",
    "plt.ylabel('Adj R^2')\n",
    "plt.show()\n",
    "# Adj R2 is increasing with features - 3 variable model is the best one on the basis of Adj R2 - we can keep adding features further to find best model \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward Stepwise Selection\n",
    "\n",
    "Forward stepwise selection is a computationally efficient alternative to best subset selection. While the best subset selection procedure considers all 2p possible models containing subsets of the p predictors, forward step- wise considers a much smaller set of models. Forward stepwise selection begins with a model containing no predictors, and then adds predictors to the model, one-at-a-time, until all of the predictors are in the model. In particular, at each step the variable that gives the greatest additional improvement to the fit is added to the model. \n",
    "\n",
    "1. Let M0 denote the null model, which contains no predictors.\n",
    "\n",
    "\n",
    "2. For k=0,...,p−1:\n",
    "* (a) Consider all p − k models that augment the predictors in Mk with one additional predictor.\n",
    "* (b) Choose the best among these p − k models, and call it Mk+1. Here best is defined as having smallest RSS or highest R .\n",
    "\n",
    "3. Select a single best model from among M0, . . . , Mp using cross-validated prediction error, Cp (AIC), BIC, or adjusted R .\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_select(y, X, feature_list):\n",
    "    remaining_predictors = [p for p in X.columns if p not in feature_list]\n",
    "    results = []\n",
    "    for p in remaining_predictors:\n",
    "        results.append(getRSS(y, X, feature_list+[p]))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_feature = 15\n",
    "models2 = pd.DataFrame(columns=[\"RSS\", \"Model\"])\n",
    "feature_list = []\n",
    "for i in range(1,max_feature+1):\n",
    "    models2.loc[i] = forward_select(y, X, feature_list)\n",
    "    feature_list = models2.loc[i][\"Model\"].model.exog_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             RSS                                              Model\n",
      "1   4.321393e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "2   3.073305e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "3   2.943854e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "4   2.826074e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "5   2.717538e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "6   2.621956e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "7   2.596156e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "8   2.516030e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "9   2.483520e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "10  2.455470e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "11  2.446856e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "12  2.441054e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "13  2.438377e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "14  2.435901e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "15  2.433831e+07  <statsmodels.regression.linear_model.Regressio...\n",
      "                                 OLS Regression Results                                \n",
      "=======================================================================================\n",
      "Dep. Variable:                 Salary   R-squared (uncentered):                   0.772\n",
      "Model:                            OLS   Adj. R-squared (uncentered):              0.769\n",
      "Method:                 Least Squares   F-statistic:                              292.7\n",
      "Date:                Sat, 24 Oct 2020   Prob (F-statistic):                    5.01e-83\n",
      "Time:                        00:02:48   Log-Likelihood:                         -1902.0\n",
      "No. Observations:                 263   AIC:                                      3810.\n",
      "Df Residuals:                     260   BIC:                                      3821.\n",
      "Df Model:                           3                                                  \n",
      "Covariance Type:            nonrobust                                                  \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Hits           2.3164      0.318      7.295      0.000       1.691       2.942\n",
      "CRBI           0.6665      0.065     10.293      0.000       0.539       0.794\n",
      "PutOuts        0.2614      0.077      3.381      0.001       0.109       0.414\n",
      "==============================================================================\n",
      "Omnibus:                      114.300   Durbin-Watson:                   2.008\n",
      "Prob(Omnibus):                  0.000   Jarque-Bera (JB):              666.384\n",
      "Skew:                           1.652   Prob(JB):                    1.98e-145\n",
      "Kurtosis:                      10.063   Cond. No.                         8.59\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] R² is computed without centering (uncentered) since the model does not contain a constant.\n",
      "[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(models2)\n",
    "print(models2.loc[3, 'Model'].summary())\n",
    "# RSS is falling with features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1     0.663357\n",
      "2     0.759667\n",
      "3     0.768905\n",
      "4     0.777294\n",
      "5     0.785017\n",
      "6     0.791772\n",
      "7     0.793015\n",
      "8     0.798617\n",
      "9     0.800436\n",
      "10    0.801910\n",
      "11    0.801822\n",
      "12    0.801504\n",
      "13    0.800929\n",
      "14    0.800332\n",
      "15    0.799698\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "rsquared_adj = models2.apply(lambda a: \n",
    "                            a[1].rsquared_adj, axis=1)\n",
    "print(rsquared_adj)\n",
    "# Adj R^2 is increasing with features\n",
    "# R^2 increasing & RSS lowering but adj R^2 is falling post 10 p hence, we dont need more than 10 p, others are insignificant.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backward Stepwise Selection\n",
    "\n",
    "Like forward stepwise selection, backward stepwise selection provides an efficient alternative to best subset selection. However, unlike forward stepwise selection, it begins with the full least squares model containing all p predictors, and then iteratively removes the least useful predictor, one-at-a-time. \n",
    "\n",
    "1. Let Mp denote the full model, which contains all p predictors.\n",
    "\n",
    "\n",
    "2. For k=p,p−1,...,1:\n",
    "\n",
    "* (a) Consider all k models that contain all but one of the predictors in Mk, for a total of k − 1 predictors.\n",
    "* (b) Choose the best among these k models, and call it Mk−1. Here best is defined as having smallest RSS or highest R.\n",
    "\n",
    "3. Select a single best model from among M0, . . . , Mp using cross-validated prediction error, Cp (AIC), BIC, or adjusted R ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backward_select(y, X, feature_list):\n",
    "    results = []\n",
    "    for combo in itertools.combinations(feature_list, len(feature_list)-1):\n",
    "        results.append(getRSS(y, X, combo))\n",
    "\n",
    "    models = pd.DataFrame(results)\n",
    "    best_model = models.loc[models['RSS'].idxmin()]\n",
    "    return best_model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "models3 = pd.DataFrame(columns=[\"RSS\", \"Model\"], index = range(1,len(X.columns)))\n",
    "feature_list = X.columns\n",
    "\n",
    "while(len(feature_list) > 1):\n",
    "    models3.loc[len(feature_list)-1] = backward_select(y, X, feature_list)\n",
    "    feature_list = models3.loc[len(feature_list)-1][\"Model\"].model.exog_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1     4.32139e+07\n",
      "2     3.13002e+07\n",
      "3     2.96395e+07\n",
      "4     2.84592e+07\n",
      "5     2.75168e+07\n",
      "6     2.67055e+07\n",
      "7     2.59362e+07\n",
      "8     2.51603e+07\n",
      "9     2.48352e+07\n",
      "10    2.45547e+07\n",
      "11    2.44686e+07\n",
      "12    2.44105e+07\n",
      "13    2.43838e+07\n",
      "14     2.4359e+07\n",
      "15    2.43383e+07\n",
      "16    2.43192e+07\n",
      "17    2.43037e+07\n",
      "18    2.42965e+07\n",
      "Name: RSS, dtype: object\n",
      "1     0.663357\n",
      "2     0.755232\n",
      "3     0.767327\n",
      "4     0.775731\n",
      "5     0.782317\n",
      "6     0.787913\n",
      "7     0.793218\n",
      "8     0.798617\n",
      "9     0.800436\n",
      "10    0.801910\n",
      "11    0.801822\n",
      "12    0.801504\n",
      "13    0.800929\n",
      "14    0.800332\n",
      "15    0.799698\n",
      "16    0.799045\n",
      "17    0.798357\n",
      "18    0.797593\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "rsquared_adj = models3.apply(lambda a: \n",
    "                            a[1].rsquared_adj, axis=1)\n",
    "print(models3.iloc[:,0])\n",
    "print(rsquared_adj)\n",
    "# R^2 increasing & RSS lowering but adj R^2 is falling post 10 p hence, we dont need more than 10 p, others are insignificant"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAATQAAAA1CAYAAADPqheOAAASCUlEQVR4Ae2dhbdXxRbH3/9g67MDno0FBqI+C7ERsUXFRr02BiYGFhYGmBiYGNggtii2KLYitoLdtd/6jG//1nDuiTn1u/O7d/Za956aM/E9e/bs2bP3/P4lgQICAYGAQCdB4F+dpB2hGQGBgEBAQIJAC0xQGoFHH31Ubhx3oxxx5BHy/vvvl84vZBAQKIpAEGhFkQvvNRDYfdDu5vzee++Vtra2xv1wEhBoNgLeC7QHHnhA3nvvvWbj0uXL+/vvv2X0mNFOOPz2228m3TnnniNXXXWV0zshUUCgDgS8FmhvvPGGbPDfDeSPP/6oo+0hzwwEjj32WLnzrjszUv3zePLkyXL++efLn3/+6ZQ+JAoI1IGAtwINIbb+BuvLQw89VEe7Q54OCMyePVt69uopX3/9dWpqvtFRRx8lCLVRl4xKTRseBgTqRMBbgTZmzBjZZttt6mx7yNsBgYsuukgOHHJgakqE2IgRI8zflGempKYNDwMCdSLgpUD79ddfpVv3bjJt2rQ62x7ydkAA+1iPHj0kCCoHsEKSDkfAS4F2++23y04779Th4DSjAhjfr732WmGF8KabbpJnnn2mGcXmKuPqq6+WvfbaK9c7zU784cwP5ZJLLzEmCur7ww8/NLsKnb48ML344otl4sSJMnbsWJnx4Qzv2lypQHvllVdk3LhxMunhScKUUVe/8ra6/3b95YYbb8j7Wkumf/vtt+XWW2+Vgw8+2NR/wIAB3rXj888/lwUWXEC++uor7+qmFWI1/JBDDpHnn39e7rvvvrDaqsBUeHx48sNy6WWXypVXXimffPKJwbvC7CvJqlKBBlNtuNGG8s477xgJfv/99+eu5McffyzzLzC/zJo9K/e7rfqCGfUmTRQExw477uBlMzbeZGO54IILvKybVmrrbbY2p2efc7bRePV+OFaHwN777G0GNrS0008/vbqMK8qpUoH23XffyW6772aqRsNffvll+f3333NV9dzzzpXNt9g81zutnnjnXXYWjO+4SaCx+UjnjTxP1lxrTR+rZur02WefSa81e5kZwrBhw4SpfKDqEQDj666/To486kihv/tGlQo0VFJ1xqST4peUl/qs18fYQvK+16rp8dtiNbfo9LxZ7cYncK6555KPPvqoWUXmKodpO3YdFpQC1YMAM6/DDj/Ma4wrFWj4junISAfVc1d40eaYbr7wwguur7R8uieffFKYIn3zzTfet2XpZZYWFmx8o7/++ktGjRrlZd18w6pMfRg0wDlvvy5TZt53KxVoeQuPpn/t9deMFtAKnTta965wjR3tuOOO6wpNDW1sUQQKCTQkNMvkLNvm+WNlJI1uueUWQQvojFQXZs3Eat/99pWttt6qmUW2K+unn37KxXPKn2GQbAdl4o1W5tVCAu3HH3+UJZZcwmhT+IudcOIJ7f4Yydl5YZdddzGOmdhfmE5ivE2iE086UTbZdJOkxy19vy7MmgkK0QAr91i5mUW2K+uee+4xfDff/PMZbTHKe8NOGCZHDz1ahgwZIqx6Lrb4Yib9lltt2S6vcCMegVbm1UICDRgIWkZI9enTx8mg/dJLL8l2A7ZLXerFB2v/A/aPR7kT3K0Ds2bCgg1l3vnm7fDNAvbcc0/Dey5uA9h1qfeKK60omDRc6a6775IrrrhChg4dKt9//73ra50mXTN4tY599AoLNL7coYcdahiLwGQXwniL82PSit5qq68mw4cPd8mq1jRffPGFfPvtt7WUUTVmrJIyrWrGjiSEPzGIUV6dBJ+kEe4ChGPNM+888sQTT6QlbTzD3HHSySc1rqMnTGXxgYTwiEfDgxCaRB40i2jbzz//bIrTIxd8Z3Vstu8XrRfTSv7SqGpetTGm3Dr20Ssl0H755RdZa+21DJPjne1CjHZJAo1pLP5OHUms4hCCxPI0I0jVVCVms2bNkpNPOdloy7jJ1B3u8+qrr5pvjX9hnUS0Cdp8mmB77rnnhGnncssv5+yErQIhWnfcjXBuJjqF6Sqk/pN0OsoqSwhGput77LGHYFphYCeaxm7jqcNPNS5L+CRir7zs8stMsXdPuFt4hmBFY1RBUKZO7DG43vrrGSf4pHyq5NU4jFUOVLmPXimBBhDTp0+XBRdaUJZaeikTDpEETtZ9RgtGXOLxOoqYmmhA/COPPCI4aNZBVWGGMFPNrO3QNhP2U0d9NU/8kNDQnn76ab2V+3jzzTc7vcPWUVna18iRI019ykRXMChoBAQ8uO666zbqR4gPU6+q9nhbe521TVggBfDduIbnIEwyKky5Hn/HeBNmRJ369evXEHzMHLbffnvzTtI/3lVhkZSG+/iJsm16GlXBq2kYV72PXmmBBhiMHDB6v837Ff74qKPkcc011yTi++mnnzaeZanLjYQ5TpjuMmLCwPvtv1+tgeJlMXvrrbfMPv50DBZatthyi1QmZjoFZmhxRTU5nGr5RhMnTcyB6pxJszqQpmZg0/hWvRc98q227b+tqRMxhkUIzYwOB4533HmH0Z7Ih/sIzAcffNBo7EXyjr5jCzSesVCBxgWNHz/eLIipvQ6BNGHCBEFLWniRhQXtWCnLF5CBTvPRd+KO8MQy3ZZpaKNxabhXlleTMK5jH71KBBqNRg2G2c8888wkXFLvs6zO+0kj+IUXXmhsdqjfEPN7nFKrIoQloyW2k4EDB8pBBx1UVdaJ+ZTBDIM1YVKbbbaZbLTxRo2RP1oYHRVNkx8xweWCVcHlV1i+odlF06ddE1/LN8JgXpRcBRp2TLT+LM9/0tEpCZ5nc4S8BF/Bc2CDVoh9ilW+EWf9s78bK7szZrS3GSJgMLPE/aHdxxECjWiG16e/bnZYIcRPXZmwnRFWtORSS5q+xGqu0jHHHGOm1ygMmGQQcmnkKtDIg4HQJea6DK/GYUzZdeyjV5lAQyCttPJKZto4derUNLxjnzFK0llUYNmJsGfQIRFiOjqtssoqDcM9NocowZSo8XF/tqan7yHMbJcSmAfGmTJlimFCRusvv/xSkzeOrJzFlQHTZlEZzOw2ozlpDC1lMrpfPvpyUzxCf9pr08yusyrI2A1FCXuVbVtCUCYZnVWLThp0NE/7yDs2PoMGDZrjOqlzvvnmm6ZzuwhPNMa555lbVl9j9VQt1a4X59SNaaVSnmgDBBp2obg/dvyIIwQaZfBNsHnS0W1Cq8KOx0o/ZhxtO5o1UzMGI/rY4MGD7dfMuc2HzC6eeuqpBs42X9svwidouKwaZ1FRXi2DcVad4p5XJtDIHHtD9/90LxTvxzSof//+8uKLL8bV00wFWXoHIFbZGE0hPgpaSJSw98AccX9xI7ktIMhrnd7rGIHJB4cQUNHtpSkbQRtXRpxgjtaR6yKYMdWiYyjR+TGiKyHMo8ROKFGXGPKxBUoSlpoXOONa89hjj+mtzCM/a2fjwzTLvo5bTSZudPDeg2XsdWONH2NmISKmU+LzmMcUgZBA2CuhlVGmEtofO6DEEXzIs7i/pJ1i7Ckngyo2Yx0kEYw2H6PJIfz5lggnJQafRRZdZA7BrQO+4soMg4FJr9UurHlw5FuzgcSzU5+VxZdY3GmKWoRXszCmHvBvVVSZQEMgMfVJGp2yKgxjo6GpBhZNT2QCxlEIpmMahZp+1tlnZW4RHc0r7nqfffdp3GZVC6aAwQiWh2A8VqiqpKKYwaDYJZTAAuaF0Li41tEd4c2UCdcabCEwkK6eoXWqQZ0BAK0urY1M//hGasjW8vMcs6acGKHRjtEI6MysfGf9pgHuJJv23dQMdnnqwsrhzJkzzSsIdra+srUZVnPjBC4vYG9jdS7uD5zjyBZofAewRJBRJt/E/v0MBg1WRhGqBxx4QCM77Lvgk0ZZU06+44477dhYxWc6qfyTlG9RXs3CGDveu+++2yiWPh2ncDQSZJxUItAAeeAOA422kVFe4mNGGT6wPULaiXmOzYHnq662qqBx8CMerOhEVXf7PZdz7Bio8Wg9dGqW1Bll0WB01YsRtcrNF8tgxuh7yqmnGCxwAbDbT1u4VsM9hnUMzwg0jrwHAyGs0YbRhCDeQ4jA3EkEU/ONwLwopQk0BBhuQPaPFWPLTFv5/uCDD8ygYwsi17qBHVs20Zkpx95mHP5KEkyu+dvpsMUR1rfrbrs2Nl8gygb3DL7J6NGjzTN4EHxZyWQQR6DhvE4EBAsU1JNZQRplCTR+I4JdkpXY26xv37562e5YhlfTMMY0xcwI4a5EG7Ej5tG09V2OlQg0OgujfVkiNEptP9G82KKalT06HqOxgsDW0BhQ4+xi0TySrtE4mDowckWN0ExjsCnhQlClj1wZzHR6zGgW9+HpKAgfJdJBaDowpxLL9vjc6eIKBm7V3jSNfURLRaDR2YtS1pZSWlfNn/YlaWgIYH7mMG5Kpe8nHZky8hui5B/VwhAkzDQYwGy8kvIqc5+6I8gZnCGODDTUD2JQpf3wO5qLi+DmG0b52GT2/3/R9nJby7fT6XlRXk3DGAWFVVzsgtHvB08WpdICDfDUq9qlEmgPcYDyLnP5JIbHnvD444+b+D3bD4rRxV4RcqlDNI09fYs+g7kZEZnKqRCNpsl7XQYzYz9L+ak4GIUtzF0I+xFaiNpuWE1LCw9iVRWBZtt0XMqpIw3tZB85RnQXooNjnlDCtpO2cMO3tu2S+l5XO5bh1SyMwRIbtT1osDCieyoWwbqUQIOZYCqYK4voiKj2aUzCKlxSfB6jFHYFl7Ky6hJ9nmTEjaar4rosZuCIZhJHTJVPO+20OaZOcemi93AzOP7441NDg3iHKQICDTtXRxNTLzqbCzGAYgDH0K4EP8Vpt/qc6ecZZ5yhl13yWJZXszBmFsFqr03MhnSAte+7nhcWaKiJK6y4gtGaOI/+oR7zC0ZoTxhO1+i5hukMt912W2Ld2GkDLaGzUh2Y2Vhh2E6bathpo+dJQtJOhzaMQINRO5IwbzCQRnlOr5myMY3GvYStohddbFFZ6N8LOa3koZkx3cT2k6atdmT7m1F23bzKLIvtwvirkgoJNLSk3r17G+MdBjzXP9wubDeBaEMYRXXVLfqs1a/rwqyZuFx/w/Wy7HLLNrPIdmXhZoDDrSvPaTq0VxfCH5LftbB99Vze60xp6uZVNGMGpTz+jK74FhJorpnnTUfoESuYgfxEgN0qqlzp9bOVoVatjIBXAo3VDfbbqsr43sofxse6Y+9gKhYoIOArAl4JNNwysNFU6TnsK/CtWK+evXpWbvNoRRxCnf1FoLRAU5+otCayOoZTp/o7JaVlbo3rBv4pnZnwCEd4pxErczjQ+rJIwooU+4/5ZCi3HYqTsHTlvaT3u/r9tNhexcYnjEsJNKaGrkusbGaHc2wW4Ul9+BGHZyVr2ee4XSQFf0cbhXE2a++r6Dt1XbMjQ7fu3Zy/d1310Hxx2kxzu9B0HF15z34nnCfHScdh4wvGhQUaHt12HCUrYKwO2X9ob3RgyLXBLOMS4K7vxYHXqvfohEQbECeKMyE7Nth46blqQT4JNNwfiBH1gXDYJESIkCGoKt7zoW2+1CEa29sqGBcWaHnjKF0FGn5ULMuzbU9nIxx40Wpdf4HIJ4HG6nOaZ30zvxXhbzhZ49/oQq6855JXV0njEttrY+ELxoUFGo2x4yjjJDg2DtW0iNF03VIHr3VfbEf2R6viHEbRbXyyNDSEe1pkRRX1cckDJ0hiJn0iIko0mLxK3vOpjR1dFzu2t1UwLiXQXOMo6bjsJcY2ycQDZhFpCIPqjO4buKZkbdUCPtjZiLNkv/uysapZeGc9J44W47BPhKB34Y+8vOdTGzu6LlmxvVo/nzAuJdC0QXUc2QMqKVC9jvLqzpPdQGgTu+66dMS66+OaP/uFEeLGjhA+EFo/ttu0357woZ6tXAfX2F4f2+itQCNekN/p1B09fQQvb51c4iXz5ll3esLRXDTKuuuh+bOqbm+NpPfDsVoEWpFXQcBbgUbl2BW0ra2t2i8VcnNGgN8l9cGG51zhkLDLI+C1QOPrsCW0i92ty3/JigHAx4tfQupMGnLFEIXsPETAe4HmIWZdpkqujqtdBpDQUO8RCALN+08UKhgQCAi4IhAEmitSIV1AICDgPQJBoHn/iUIFAwIBAVcEgkBzRSqkCwgEBLxHIAg07z9RqGBAICDgikAQaK5IhXQBgYCA9wj8D+LUuAMNCl1DAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shrinkage Methods\n",
    "\n",
    "### Ridge Regression\n",
    "\n",
    "Ridge regression is very similar to least squares, except that the coefficients are estimated by minimizing a slightly different quantity. In particular, the ridge regression coefficient estimates βˆR are the values that minimize\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "where λ ≥ 0 is a tuning parameter, to be determined separately. Equation trades off two different criteria. As with least squares, ridge regression seeks coefficient estimates that fit the data well, by making the RSS small. However, the second term, called a shrinkage penalty, is small when β1, . . . , βp are close to zero, and so it has the effect of shrinking the estimates of βj towards zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import scale \n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import Ridge, RidgeCV, Lasso, LassoCV\n",
    "from sklearn.metrics import mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hitters = pd.read_csv('../data/Hitters.csv', header=0, na_values='NA')\n",
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']])\n",
    "\n",
    "y = Hitters.Salary  # the response variable \n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test , y_train, y_test = train_test_split(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AtBat           0.094214\n",
      "Hits            0.425688\n",
      "HmRun           1.393435\n",
      "Runs            0.667159\n",
      "RBI             0.696568\n",
      "Walks           1.053478\n",
      "Years           2.479727\n",
      "CAtBat          0.008734\n",
      "CHits           0.035571\n",
      "CHmRun          0.290251\n",
      "CRuns           0.078347\n",
      "CRBI            0.071455\n",
      "CWalks          0.080231\n",
      "PutOuts         0.073132\n",
      "Assists         0.009901\n",
      "Errors         -0.120890\n",
      "League_A       -8.035691\n",
      "Division_E     38.755126\n",
      "NewLeague_A    -5.881910\n",
      "dtype: float64\n",
      "75728.85761302235\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Ridge(alpha=4, copy_X=True, fit_intercept=True, max_iter=None, normalize=True,\n",
       "      random_state=None, solver='auto', tol=0.001)"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridge = Ridge(fit_intercept=True, normalize=True, alpha=4)\n",
    "ridgemodel = ridge.fit(X_train, y_train)             # Fit a ridge regression on the training data\n",
    "pred = ridge.predict(X_test)           # Use this model to predict the test data\n",
    "print(pd.Series(ridge.coef_, index=X.columns)) # Print coefficients\n",
    "print(mean_squared_error(y_test, pred))        # Calculate the test MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "RidgeCV?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridgecv =  RidgeCV(alpha, normalize = True)\n",
    "ridgecv.fit(X_train, y_train)\n",
    "ridgecv.alpha_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80032.16337972828"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridge_best = Ridge(alpha=ridgecv.alpha_, normalize=True)\n",
    "ridge_best.fit(X_train, y_train)\n",
    "mean_squared_error(y_test, ridge_best.predict(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AtBat           -1.937886\n",
       "Hits             6.125915\n",
       "HmRun           -0.915264\n",
       "Runs            -1.018186\n",
       "RBI              0.791576\n",
       "Walks            6.247442\n",
       "Years          -15.849746\n",
       "CAtBat          -0.085709\n",
       "CHits            0.251573\n",
       "CHmRun           1.654119\n",
       "CRuns            0.985110\n",
       "CRBI             0.041817\n",
       "CWalks          -0.539329\n",
       "PutOuts          0.256748\n",
       "Assists          0.316338\n",
       "Errors          -1.251228\n",
       "League_A       -98.739059\n",
       "Division_E     144.956790\n",
       "NewLeague_A     56.504161\n",
       "dtype: float64"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(ridge_best.coef_, index=X.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lasso Regression\n",
    "\n",
    "As with ridge regression, the lasso shrinks the coefficient estimates towards zero. However, in the case of the lasso, the l1 penalty has the effect of forcing some of the coefficient estimates to be exactly equal to zero when the tuning parameter λ is sufficiently large. Hence, much like best subset se- lection, the lasso performs variable selection. As a result, models generated from the lasso are generally much easier to interpret than those produced by ridge regression. We say that the lasso yields sparse models—that is, models that involve only a subset of the variables. As in ridge regression, selecting a good value of λ for the lasso is critical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "76680.11673280169"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## We saw that ridge regression with a wise choice of λ can outperform least squares as well as the null model on the Hitters data set. \n",
    "## We now ask whether the lasso can yield either a more accurate or a more interpretable model than ridge regression.\n",
    "\n",
    "lassocv = LassoCV(alphas=None, cv=10, max_iter=1e5, normalize=True)\n",
    "lassocv.fit(X_train, y_train)\n",
    "\n",
    "lasso.set_params(alpha=lassocv.alpha_)\n",
    "lasso.fit(X_train, y_train)\n",
    "mean_squared_error(y_test, lasso.predict(X_test)) ## Its lower than the ridge model- more accurate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AtBat           -1.792745\n",
       "Hits             5.933814\n",
       "HmRun           -0.000000\n",
       "Runs            -0.000000\n",
       "RBI              0.000000\n",
       "Walks            5.254110\n",
       "Years          -13.832284\n",
       "CAtBat          -0.000000\n",
       "CHits            0.000000\n",
       "CHmRun           1.540912\n",
       "CRuns            0.820241\n",
       "CRBI             0.000000\n",
       "CWalks          -0.400841\n",
       "PutOuts          0.257400\n",
       "Assists          0.201540\n",
       "Errors          -0.000000\n",
       "League_A       -45.697486\n",
       "Division_E     143.328128\n",
       "NewLeague_A      0.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Some of the coefficients reduced to exact zero like RBI, CRBI- more interpretable\n",
    "pd.Series(lasso.coef_, index=X.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimension Reduction Methods- PCR and PLS Regression\n",
    "\n",
    "The methods that we have discussed so far in this chapter have controlled variance in two different ways, either by using a subset of the original vari- ables, or by shrinking their coefficients toward zero. All of these methods are defined using the original predictors, X1,X2,...,Xp. We now explore a class of approaches that transform the predictors and then fit a least squares model using the transformed variables. We will refer to these tech- niques as dimension reduction methods.\n",
    "\n",
    "The term dimension reduction comes from the fact that this approach reduces the problem of estimating the p+1 coefficients β0,β1,...,βp to the simpler problem of estimating the M + 1 coefficients θ0, θ1, . . . , θM , where M < p. In other words, the dimension of the problem has been reduced from p+1 to M +1.\n",
    "\n",
    "### Principal Components Regression (PCR)\n",
    "Principal components analysis (PCA) is a popular approach for deriving a low-dimensional set of features from a large set of variables. The first principal component/feature direction of the data is that along which the observations vary the most. \n",
    "\n",
    "The principal components regression (PCR) approach is PCA + regression. It involves constructing the first M principal components,Z1,...,ZM, and then using these components as the predictors in a linear regression model that is fit using least squares. The key idea is that often a small number of principal components suffice to explain most of the variability in the data, as well as the relationship with the response."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hitters = pd.read_csv('../data/Hitters.csv', header=0, na_values='NA')\n",
    "Hitters = Hitters.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "y = Hitters.Salary  # the response varddiable \n",
    "dummies = pd.get_dummies(Hitters[['League', 'Division', 'NewLeague']])\n",
    "X_prep = Hitters.drop (['Salary', 'League', 'Division', 'NewLeague'], axis = 1).astype('float64')\n",
    "X = pd.concat([X_prep,  dummies[['League_A', 'Division_E', 'NewLeague_A']]], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 38.31423957,  60.15531541,  70.8416746 ,  79.03419435,\n",
       "        84.29027516,  88.63477938,  92.26288766,  94.96304322,\n",
       "        96.28269124,  97.25541298,  97.97775424,  98.64870031,\n",
       "        99.15178693,  99.4730334 ,  99.74659124,  99.89398793,\n",
       "        99.96815949,  99.99375108, 100.        ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Scale the variables as unit of the variables are different and generate principal components\n",
    "pca = PCA()\n",
    "X_pca = pca.fit_transform(scale(X))\n",
    "regr = linear_model.LinearRegression()\n",
    "\n",
    "## Variation explained by each component\n",
    "np.cumsum(pca.explained_variance_ratio_) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pcr(X,y,pc):\n",
    "    \n",
    "    # Define PCA object\n",
    "    pca = PCA()\n",
    "\n",
    "    # Preprocessing first derivative-Scale the variables as unit of the variables are different and generate principal components\n",
    "\n",
    "    X_pca = pca.fit_transform(scale(X))[:,:pc]\n",
    "\n",
    "    # Define linear regression object\n",
    "    regr = linear_model.LinearRegression()\n",
    "    \n",
    "    # Fit\n",
    "    regr.fit(X_pca, y)\n",
    "\n",
    "    # Predict\n",
    "    y_train = regr.predict(X_pca)\n",
    "\n",
    "    # Cross-validation\n",
    "    y_cv = cross_val_predict(regr, X_pca, y, cv=20)\n",
    "\n",
    "    # Calculate scores for training and cross-validation\n",
    "    score_train = r2_score(y, y_train)\n",
    "    score_cv = r2_score(y, y_cv)\n",
    "\n",
    "    # Calculate mean square error for training and cross validation\n",
    "    mse_train = mean_squared_error(y, y_train)\n",
    "    mse_cv = mean_squared_error(y, y_cv)\n",
    "\n",
    "    return(y_cv, score_train, score_cv, mse_train, mse_cv)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABI/klEQVR4nO3dd3hUVfrA8e+bSUhooYRAgNCL0lsoKgoKCgIqKgrYQLDrYll/iKur6Lq74lpW1F3WgsCuUtRFWREQEQQLJSi9NyVIM/ROkvf3x70Dk2TSZzIp7+d55pk757Z3JpO8Oeeee46oKsYYY0yghYU6AGOMMSWTJRhjjDFBYQnGGGNMUFiCMcYYExSWYIwxxgSFJRhjjDFBYQnGmCJIRHaISM9Qx2FMQViCMSaIRKSriHwvIodF5ICIfCciHUMdlzGFITzUARhTUolINPA5cD8wDSgDXAqcDuI5w1U1JVjHNyYvrAZjTPA0BVDVyaqaqqonVfVLVV0lIo1E5GsRSRaR30TkAxGp7O8gItJJRH4QkUMisltE3hSRMj7rVUQeFJHNwGYReUtEXslwjBki8mgw36wxGVmCMSZ4NgGpIjJRRK4WkSo+6wT4K1ALaAbUAUZncZxU4FGgGnAR0AN4IMM2/YHOQHNgIjBYRMIARKQa0BP4sOBvyZjcswRjTJCo6hGgK6DAO8B+tyZRQ1W3qOpcVT2tqvuBV4FuWRxnuaouVtUUVd0B/MvPtn9V1QNuLWkpcBgnEQEMAhao6t7Av0tjsmYJxpggUtX1qjpUVeOBljg1lr+LSA0RmSIiu0TkCPAfnBpKJiLSVEQ+F5E97rZ/8bPtzgyvJwK3ucu3Af8O1HsyJrcswRhTSFR1AzABJ9H8Badm00pVo3GSgGSx6z+BDUATd9s/+Nk247Do/wGuE5E2OE1wnwbgLRiTJ5ZgjAkSEblQRH4vIvHu6zrAYGAxUBE4BhwWkdrA/2VzqIrAEeCYiFyI0ystW6qaBCzDqbl8oqonC/RmjMkHSzDGBM9RnAvvS0TkOE5iWQP8HngOaI9zrWQm8N9sjvM4cIt7vHeAqbk8/0SgFdY8ZkJEbMIxY0omEbkMp6msntovugkBq8EYUwKJSATwMPCuJRcTKpZgjClhRKQZcAioCfw9pMGYUs2ayIwxxgSF1WCMMcYEhQ126apWrZrWr18/1GEYY0yxsnz58t9UNdbfOkswrvr165OYmBjqMIwxplgRkZ+zWmdNZMYYY4LCEowxxpigsARjjDEmKCzBGGOMCQpLMMYYY4LCepEZY0Kiz+uLWLf7SKby5jWj+eLhS0MQkQk0q8EYY0Kifd3KRHjST2sT4RHa16uSxR6muLEEY4wJiRE9mhAm6ROMR4QRPRqHKCITaJZgjDGFLjVNWf7zQaKjzrfSewQGJNShesWoEEZmAsmuwRhjCs3RU2eZlpjEhO+3s/PASWpVjiI8TEhJU1IVBibEhzpEE0BWgzHGBN3OAyf40+fruPivX/Onz9cRFx3FuNvas2jkFQzqWAfBqcG8MncTNsJ7yWE1GGNMUKgqP/5ykPe+3c7sNXsIE6Fv65oM79qA1vGVz203okcTNu07xqWNq/HK3E3898dd3NjBajIlgSUYY0y+ZNXNuFnNitzfvTHvfbudlTsPER0Vzj2XNWLIxfWoWalspu2rR0cx7d6LSEtTvtm0n+c/X8elTavZtZgSwJrIjDH54q+bsUfg5+QTjJj8E0dOnuVP17Vg8R96MOrqC/0mF19hYcKYAa05eTaVZz5dG8zQTSGxBGOMyRd/3YxTFZrFRfPuHQnMe6wbt19Un3Jlct9Q0ii2Ao/2bMrstXv4YvXuQIdsCpklGGNMnv2cfJypy3YSGX7+T0iYQN9WcXzywMX0bF6DsDDJ5ghZu/vSBrSqXYlnPlvDweNnAhWyCQFLMMaYXNl75BTvfbud6978lm5/W8ArczfRsFoFwt1EUsYTxrPXtijwecI9YYy5sTWHTpzl+c/XFfh4JnSClmBEZLyI7BORNT5lfxORDSKySkSmi0hlt7y+iJwUkRXuY5zPPh1EZLWIbBGRsSJOnVxEqorIXBHZ7D5XccvF3W6Le572wXqPxpR0h06cYfLSXxj89mK6/HUef/p8HSlpyh/6XMj3o67g04cucboZB/gmyea1onng8sZM/2kXX2/YG5BjmsInwepzLiKXAceASara0i27CvhaVVNEZAyAqj4hIvWBz73bZTjOUmAEsAT4AhirqrNE5CXggKq+KCKjgCrusfoAvwP6AJ2B11W1c07xJiQkaGmbMtkGGyzdsvr5XxhXkfu7N2LGil9ZuHk/Z1OVhtXKc02bWlzbthaNYiuk237fkVM8NPkn3rylXUB7fp1JSaPfG4s4cjKFLx+7jOioiIAd2wSOiCxX1QR/64LWTVlVF7qJw7fsS5+Xi4EB2R1DRGoC0aq62H09CegPzAKuA7q7m04EFgBPuOWT1Mmci0WksojUVNUid8Uw1H/g29etzOZ9Rzmbev6fDBtssPTw9/MXYPPeozw8ZQVx0VEMvbg+17apTcva0Yj4v6bi7WYcaGXCw3hpQBtu+Md3/PWLDfz1hlYBP4cJrlDeBzMMmOrzuoGI/AQcAZ5W1UVAbSDJZ5sktwyghk/S2APUcJdrAzv97JMpwYjIPcA9AHXr1i3Qm8mPgv6Bz0uCOnEmhaSDJ0k6eIKdB5znpIMnSUlNX4O1wQZLjxE9mjAtMQk4/x1Q4Nq2tRnUsQ4d61fN94X6QGlbpzJ3XdqQtxdu45rWNbm4cbWQxmPyJiQJRkSeAlKAD9yi3UBdVU0WkQ7ApyKS66uFqqoikue2PlV9G3gbnCayvO5fUCN6NOGj5el/wQE61a/CD1uTKRMuRHjCKBMe5jx7nOcIjxARHkbr+EqZEpQnTKgQGc5fZ613EsoBJ5EkZ+iNExkeRnyVssRVimLvkVOkuYeoF1OOipHWFFGSpaUp32zez8Tvd3AmNe1cuSdMGNSxDn++vmjVFB7t2ZQv1+7hif+uYs4jl+Wp27MJrUL/SYnIUKAf0MNtxkJVTwOn3eXlIrIVaArsAnzHjIh3ywD2epu+3Ka0fW75LqBOFvsUKdWjo7ipQzwfLPnlXIo5m6qMmLIi38dMTVOW7jjAip2HqF2lLPFVynJVrUrEu8t1qpYjvkpZYitEIiLsO3KKS1+az+mUNDxhwsa9x+gzdhEv39SaDvWqBuR9mqLh6KmzfLw8iUk//Mz2344TWzGSuy9twKQffuZ0ShoRYcLDPZuEOsxMypbxMObG1gx8ezEvz9nEM9c0D3VIJpcKNcGISG9gJNBNVU/4lMfiXLBPFZGGQBNgm6oeEJEjItIF5yL/HcAb7m4zgCHAi+7zZz7lD4nIFJyL/IeL4vUXL28t5nRKGmU8Ybx9ewcqRIVzJjWNs6nK2ZQ0zqamnX+d6r5Occpmr97Dml8Pk6bgCYOezWrw3LUtqV4xMlfNG+eS3NJfGNypLn1axvF/H6/ipnE/cM9ljXj0yiZEhnsK4ZMwwbJ1/zEmfb+Dj5cncfxMKu3qVub1QW25umVNyoSHcfJMKh8s/aVID5XfuWEMt3epx/vfb6dv6zj756eYCGYvssk4F+GrAXuBZ4EngUgg2d1ssareJyI3As8DZ4E04FlV/Z97nARgAlAW5+L+79wmsRhgGlAX+Bm42U1IArwJ9AZOAHeqao7dw0LZi+zp6av5YOkv3Nq5Hi/0z9SRLlu+NZCo8DAWPnF5nv9IZOwFdPTUWf48cz1Tlu2kaY0KvHpzW1rWrpSnY5rQSktTFmzax/vf7WDR5t8o4wmjX+uaDLm4Pm3qVE63bbB6gQXasdMp9HptIVERYcwccSlREfaPT1GQXS+yoCWY4iaUCaagv+AFSVDZmb9hH098sooDx8/w0BWNefDyxkR47N7coiKrTh41oiOJivDwc/IJakRHclvnegzqVJfYipEhiDKwFm7azx3jl/Lg5Y34v14XhjocQ/YJxv5aFAHebp75/e9xRI8mdKxfNeC9vy6/sDpfPnoZfVvX5O9fbeb6f3zHpr1HA3oOk3/+BpsE2HvkNLEVInljcDu+feIKftejSYlILgCXNY3lpg7xjPtmG2t2HQ51OCYHVoNxlcYbLfNi1urdPPXpGo6dSuH3VzXlrksb4glxF9bSbt+RU3R9aT5nUnx6gglMuLMTlzaNDWFkwXX4xFl6vvYN1SpEMuOhS6xWHWIhudHSlCxXt6pJxwZVeWr6av46awNfrtvLyze1oUG18qEOrVQ6cuos7327nRSf5BLhEQZ2rFuikwtApXIRvNC/Jff+eznjFmzldz3y1/Mt1Dc6lwaWYEyuVasQybjbOvDZil955rM1XPHyAvzVf+0XNHhSUtOYsmwnr83dRPLxM/RpVZOv1u/lTEpaqbpJtleLOPq1rskbX2+hV8s4mtaomOdj2EgWwWcJxuSJiNC/XW26NIzh+n98x+7Dp9Ktt1/Q4Fm4aT8vzFzHpr3H6FS/KhPubE6r+ErnOnkU5W7GwbBp71HOpKZx1WsL05X7+wdHVUk+foZfD53k10Mn2XXoFKmqpKbZSBbBZAnG5EtcpSg+feBiLhkznxSfX1L7BQ28LfuO8eeZ65i/cT91q5bjn7e2p3fLuHNjg3nntC9tn3un+lXZuv8YPoMREB4mVCobwWtzNznJ5PBJfj10il2HTqa7VgVQNsJDhchwjp5KQXHmsyltSTrYLMGYfKtRqSyDOtZhauLOc80MbepUtl/QADl4/Ax//2oT/1nyC+UiPDx59YUMvaR+phtfgzXYZFHnvUk5Ne184khJU37Ylszi7cnUqBhFrcpRtKgVzZXNa1CrUhS1KpelVuWy1K5clsrlIth/9PS5+8jSFAZ3rJPNGU1eWYIxBeI7nlqYwJLtB3jh83U82aeZ9TLLpzMpaUz6YQdj523m2OkUBneqy6NXNqVahZLR1ThQvKNQeP/B8YRBjwtr8Md+zakRHUWZ8Jx7l/kO1wTw4dJfitxYbMWZJRhTIBmHmgkPE979djs7ko/z+qB2lI+0r1hWsurFVMYjnElVLm1Sjaf7NueCuLxfwC4tfP/BiQgL44XrW+a5Bu1tYqxTpSxTl+3kvm6NqFO1XHACLmWsA7kpMO+Nng/3bMJz17XkuWtb8PWGfQwY9wO/HjoZ6vCKrKxulIyK8PD+0I5MGtbJkksOvP/gFGRGTW8T48jeF+IJE/7+1eYgRFo6WYIxBZZxJIIhF9dn/NCOJB04wXVvfceKnYdCG2ARNaJHE4T0CSY8TJj9yGVcfmH1LCf4MukFaiSLGtFR3N6lHtN/SmLLvmMBiq50swRjgqL7BdX55IGLiQwPY+C/fuCL1UV2QOuQ+Dn5OC/N2ZhuPpYIjzCoU11qVS4bwsiKn4IOteTr/u6NiIrw8NpXmwIQmbEEY4KmaY2KfPrgJbSsXYkHPviRt+ZvobQPTbTzwAlGfrySK175hv+t/JVBHesQ6V6Mti7eoRdTIZJhlzRg5qrdrPs18/UxkzeWYExQVasQyQd3daZ/21r8bc5Gfj9tJadTUkMdVqFLOniCJ/+7istfXsCnK37l9i71WDTycl68sXWBryGYwLr7soZER4Xz6tyNoQ6l2LMuPibooiI8vDawLQ1jK/Dq3E3sPHiCf92eQNXyZUIdWtDtOnSSt+Zv4aPEnQjCLZ3r8kD3xsRVOp9ISuuNkkVVpbIR3HNZQ17+chM//nKQ9nVtZIr8stGUXTaacuGYsfJXHv9oJXHRUYwfmkDj6iWzl9Tuwyf5x/ytTFnm3F8xsGMdHuje2K6vFBPHT6dw2UvzaVYzmv/c1TnU4RRpNpqyKTKubVOL+CpluWdSIle9tpA0P//fFJfBMrO6j6VKuQiOn04lTZWbO9bhwcsbU9sSS7FSPjKc+7s34oWZ6/lhazIXNYoJdUjFkiWYUBrXFfaszlwe1wru+7bw4ykk7etW4dMHL6Hf2G85dPJsunXFabBMf6PxAhw6cZaBbmKxG/aKr9u61OOdRdt45cuNfHTfRdZtPB/sIn8oxXcCT4brEJ4yTnkJF1+lHP994GIyjiZTnHpSPXR55jg9Av+9/2JevLG1JZdiLirCw++uaELizwdZsGl/qMMploKWYERkvIjsE5E1PmV/E5ENIrJKRKaLSGWfdU+KyBYR2SgivXzKe7tlW0RklE95AxFZ4pZPFZEybnmk+3qLu75+sN5jgXUbCZLhRyBh0O2J0MRTyBrGVmBwp7rpkkyTGhUyDeZY1Gzae5QXZ23g+n9+n672Eu4RBneuR7tiUgMzObs5oQ7xVcryypcbS30X+/wIZg1mAtA7Q9lcoKWqtgY2AU8CiEhzYBDQwt3nHyLiEREP8BZwNdAcGOxuCzAGeE1VGwMHgeFu+XDgoFv+mrtd0VQxDtreCuLzB7XuxVCxRuhiKmQP92hybsrbMIHVu47QdczXjJ23maOnzuawd+HZd/QU7y7aRt+xi7jqtYW8s2gbzWpG86f+Lc/dxxJejGpfJnfKhIfxSM+mrNl1hDlr94Q6nGInaAlGVRcCBzKUfamqKe7LxUC8u3wdMEVVT6vqdmAL0Ml9bFHVbap6BpgCXCdOY+gVwMfu/hOB/j7Hmugufwz0kKLceNptJHgi3BcC276Gj4fBiQPZ7lZS+I4ldUvneswc0ZXODWJ4de4mLn1pPv9csJUTZ1JyPlAQnDiTwqc/7eKO8Uvp8pd5vDBzPZ4w4dlrmrPkDz0YP7Qjt3epZ/exlHD929aiYWx5Xp27KdMEZSZ7obzIPwyY6i7Xxkk4XkluGcDODOWdgRjgkE+y8t2+tncfVU0RkcPu9r9lDEBE7gHuAahbt24B304+eWsxy9+H9kMhuiZ8MwZ2fAfXvgFNrwpNXIXI9z6Q6hWjeHdIAit3HuK1rzYxZvYG3l20jfu7N+K2LvWIighc81lWvcDqVS1Hh/pVmLNmD8fPpFK7clnu796I69vF07h6hWzjNyVPuCeMx65sykMf/sT/Vv5K/3a1c97JACFKMCLyFJACfBCK83up6tvA2+DcBxOyQLqNhP3rofsop3msyVUw/T748CZofwf0+gtElsz7RcD/hFlt6lRmwp2dWP7zAV6bu5kXZq7nXwu38WD3RgzqVDcgiSarXmA/HzjBgRNnuKZNLa5vV5uO9asSls3cNqV1wq/SpE/LmjSruZW/f7WJvq1rnmvWNdkr9AQjIkOBfkAPPX/VbBfgO5VcvFtGFuXJQGURCXdrMb7be4+VJCLhQCV3+6KrYhzcOev861pt4d5vYP6f4buxsG0BXPcPaFD07w0JtA71qvKfuzqzeFsyr87dxOj/rXMSzeWN+WDJz6zffTTTPlndR6Oq7Dlyii37jrF57zGOnU5NN90zONeBXujfihva1w5obckUb2Fhwu+vbMpdkxL5ZHkSgzqFqMWjmCnUBCMivYGRQDdVPeGzagbwoYi8CtQCmgBLAQGaiEgDnMQxCLhFVVVE5gMDcK7LDAE+8znWEOAHd/3XWhy7f4RHwpXPwwV9Yfq9MLEfdHkAejwDEaXvpr0uDWOYek8Xvt+azCtfbuTpT9dQvowHTxjp5mSP8Ajt61bml+QTbNl/lM17j7F53zG2uI9jp89fz6lcLoJq5cuQfPwMaersO7BjXW7pbH88TGY9mlWnbZ3KjJ23mevb1y7yvR2LgqANFSMik4HuQDVgL/AsTq+xSM7XKBar6n3u9k/hXJdJAR5R1VlueR/g74AHGK+qf3bLG+Ikl6rAT8BtqnpaRKKAfwPtcDoZDFLVbTnFW6SHijlzHOY+C8vegWpN4fpxULtDqKMKGVXlm037GTN7Q6YajABlwoXTKee/17EVI2lSvQKNq1dwnyvSpEYFYsqXSTcne1R4GAufuNwu1Jssfbv5N257bwmjr2nO0EsahDqcAsnqGmReR9LIbqgYG4vMVaQTjNfW+fDZg3B0D1z6GFw2EsJL/oCRWVFV7nx/Wbqb4GpWiqJPq5rnEkrj6hWoXC77z+jp6av5YOkv3Nq5Hi/0bxnssE0xpqoMfmcxW/YdZ+HI7pQrU3wHQ3l6+mqmJu5Mdw3SW4vPy++BJZhcKBYJBuDkIZj9JKz8EMKjIOVU5m1K+FAzvvYdOVXgGsi+I6d4aPJPvHlLO6u9mBwl7jjAgHE/MOrqC7mvW6NQh5Nvvr87Xvn5HbLBLkuSspXh+n/ChX3hk+GZ15eSoWa8vPfRfLD0lwLPyW5MbiTUr0r3C2IZ981Wbulcl+ioiJx38iNQTVT5FV02gnox5di015keOsIjAb+Xy/raFVfN+sHd80v1UDNegZqT3Zjc+v2VF3DoxFnGf7s938doX7cyEZ703d8La7DXLfuO0v+t79i09xgetwt+MMYBtARTnNVoDh2GQphPbxZVWPo2HE4KWViFLZBzshuTG63iK9G7RRzvLtrOweNn8nWMe7o1zFQW7MFeVZVpiTu55o3v2Hf0NO8P7cjgjnWCNhKFNZEVd92egBUfQlqq0zxWvyssegW+fc1pRut0j1NWhEfLMaY4euyqpsxZt4d/LdzGqKsvzHH75GOnWf7zQZb/fJBlOw6wZteRTDf5xlWK4uSZ4EwpfvTUWZ6avoYZK3/looYx/H1QW2pER9GiVnTQRqKwi/yuYnOR35/PH3OGmulwJ/R7FQ7ugMTx8OMkOHkQYptBp7uh9UCIzDzUiTEm77K7hjJzRFe2/Xac5TsOkvjzARJ3HGTbb8cBKOMJo1V8JRLqVaFx9Qo8/ekaTqek4QkTwt37um7pXJffXdGE2IqRAYl15c5D/G7yT+w6dJJHezbh/u6NzzWNFZT1IsuFYp1gju6Bj++EARPSj8R89iSs+QSW/Av2rILIStDuVuh4F8QU394vxhQFT09fzZTEnaT41EI8AjUrl+XEmVQOuE1nlctFkFCvCh3qVSWhfhVa1a6UbpQI327yI65ozOvzNjNl2U4iw8MY3rUB91zWkIr57EiQlqa8++02Xpq9kRrRUbw+qC0J9asW7I1nYAkmF4p1gsmJKiQtc67NrP0U0s5C456QvMWp7WRUiro5G5Nf/rr5AsRXKUvnBjF0rF+FhPpVaFitQrZj2fnrJr/9t+O88uVGPl+1myrlInjw8sZ5Huz1t2OneWzaShZu2k/vFnGMubE1lcrlL1FlxxJMLpToBOPr6F74caLThHZ0d+b1njLQ7nanqc0Yky3fmxU9YcL17Wrx8k1tA3b81UmHeWnOBhZt/o3alcvySM8m3NA+PsfmrW83/8aj01Zw+ORZnunXnFs71w3alM+WYHKh1CQYr9Sz8NN/YOZjoD7/gYVHwcOrStWkZ8bkVyBu9M2N77b8xpjZG1iVdJimNSrwf70upGez6pmSxtnUNF6du4lx32ylUWwF3rylHRfGRQc8Hl92o6XJzBMBCXfCntVOZ4A0d/bImCZQLrBttMaUVIG40Tc3Lmlcjc8evIRZa/bw8pyN3D0pkbIRHk6e9d/jbFDHOjxzTfOQD2Vj98GUdt1Gnr+PRjywdzW836dU3UdjTEEU1o2+IkKfVjWZ8+hl/OX6Vij+W5+6No7hxRtbhzy5gCUY451RU8KcmzYHjId962FcV9g0J9TRGVPkFfaNvhGeMG7pXJc5j1xGeIZrMZHhYbw6sG2hxJEblmCMU4up28W5abPljc5kZ5Xi4cOb4cs/OtdrjDFFSr2Y8gzqWOfccDMRHuGmIDbT5YclGHN+Rk3vhf2YRjD8K0gYDt+PtSYzY4qoET2aECbBG0usoCzBGP8iopyuytZkZkyR5e1kEKyxxArKEozJnjWZGVOkFeXRxC3BmJxZk5kxRVZRHk3cEozJHWsyM8bkUdA6SovIeKAfsE9VW7plNwGjgWZAJ1VNdMvrA+uBje7ui1X1PnddB2ACUBb4AnhYVVVEqgJTgfrADuBmVT0ozq2trwN9gBPAUFX9MVjvs9RpeSPUbAsfDXGazPyxscyMMQS3BjMB6J2hbA1wA7DQz/ZbVbWt+7jPp/yfwN1AE/fhPeYoYJ6qNgHmua8BrvbZ9h53fxNI3iazWD9zYJSyKZuNMVkLWoJR1YXAgQxl61V1Yxa7ZCIiNYFoVV2szqBpk4D+7urrgInu8sQM5ZPUsRio7B7HBFJEFNzxGYRlGJ21FE7ZbIzxryhdg2kgIj+JyDcicqlbVhvwvZqc5JYB1FBV73DAe4AaPvvszGKfdETkHhFJFJHE/fv3B+RNlCoV46D9HemTTOuBNlCmMQYoOglmN1BXVdsBjwEfikiuhwB1azd5HhZaVd9W1QRVTYiNjc3r7gbSj2UGcGQXpKVlvb0xptQoEglGVU+rarK7vBzYCjQFdgHxPpvGu2UAe71NX+7zPrd8F1Ani31MoPmOZVanC2z5Cua/EOqojDFFQJFIMCISKyIed7khzgX6bW4T2BER6eL2DrsD+MzdbQYwxF0ekqH8DnF0AQ77NKWZYPCOZXbzRGg/BBa9AqumhToqY0yIBbOb8mSgO1BNRJKAZ3Eu+r8BxAIzRWSFqvYCLgOeF5GzQBpwn6p6Owg8wPluyrPcB8CLwDQRGQ78DHj7zH6B00V5C0435TuD9R6NyzuWGUCflyF5K3z2EFRpAHU6hjY2Y0zI2IyWrlI3o2UwnTgA71wBZ47D3V9D5To572OMKZaym9GySDSRmRKmXFW4ZSqknIIpg51EY4wpdSzBmOCIvQAGvA9718J/77GeZcaUQpZgTPA06QlX/Rk2fA7z/xzqaIwxhSz0kzabkq3L/bB/Ayx62anVtM5i/DJjTIljNRgTXCJOz7J6XZ2eZUnWkcKY0sISjAm+8DIw8N8QXRMmD7a5ZIwpJSzBmMJRrioMdnuWTR5kPcuMKQUswZjCU/1C61lmTCliCcYULutZZkypYb3ITOHrcj/sX+/2LLsQWt8U6oiMMUFgQ8W4bKiYQpZyBsbUh7N+rsXYlMvGFBs2VIwpesLLQIvrM5fblMvGlBiWYEzo9PgjeCLTl9mUy8aUGJZgTOhUjIN2t0GY91KgQJvBNuWyMSWEJRgTWt1G+iQYhZMHwa4LGlMiWIIxoeU75XKtDrDuU/j6T6GOyhgTAJZgTOh5p1we/CF0GOpMubz0nVBHZYwpoGwTjIjc5rN8SYZ1DwUrKFPKeKdcrhgHfV6BplfDF/8H6z8PdWTGmALIqQbzmM/yGxnWDQtwLMaAJxwGjIfaHeCT4fDLklBHZIzJp5wSjGSx7O91+pUi40Vkn4is8Sm7SUTWikiaiCRk2P5JEdkiIhtFpJdPeW+3bIuIjPIpbyAiS9zyqSJSxi2PdF9vcdfXz+E9mqKmTDlnyuXo2jB5IPy2OdQRGWPyIacEo1ks+3ud0QSgd4ayNcANwELfQhFpDgwCWrj7/ENEPCLiAd4CrgaaA4PdbQHGAK+pamPgIDDcLR8OHHTLX3O3M8VN+Wpw2ydOD7P/3ABH94Q6ImNMHuWUYC4UkVUistpn2fv6gux2VNWFwIEMZetVdaOfza8DpqjqaVXdDmwBOrmPLaq6TVXPAFOA60REgCuAj939JwL9fY410V3+GOjhbm+Km6oN4JZpcDwZPrgJTh8NdUTGmDzIabDLZoUSBdQGFvu8TnLLAHZmKO8MxACHVDXFz/a1vfuoaoqIHHa3/y3jSUXkHuAegLp16wbkjZgAq90ebp4IHw6EaXc4c8qElwl1VMaYXMi2BqOqP/s+gGNAe6Ca+7pYU9W3VTVBVRNiY2NDHY7JSpMr4dqxsPVrmPE7uxHTmGIip27Kn4tIS3e5Js41lGHAv0XkkQDGsQuo4/M63i3LqjwZqCwi4RnK0x3LXV/J3d4UZ+1ug8ufhlVTYN7zoY7GGJMLOV2DaaCq3l5gdwJzVfUanGaqQHZTngEMcnuANQCaAEuBZUATt8dYGZyOADPUmWNgPjDA3X8I8JnPsYa4ywOAr9XmJCgZLnvcuRHz21ftRkxjioGcrsGc9VnuAbwDoKpHRSTb+W5FZDLQHagmIknAszgX/d8AYoGZIrJCVXup6loRmQasA1KAB1U11T3OQ8AcwAOMV9W17imeAKaIyAvAT8B7bvl7ODWsLe75BuXwHk1xIeLciHl0r3MjZsWa0KxfqKMyxmQh2wnHROR/wJc4F9HH49RoDolIWSBRVVsUTpjBZxOOFSNnTsDEa2DvGrhjBtTtHOqIjCm1sptwLKcazHDgeaAnMFBVD7nlXYD3AxahMXnhvRHz1WYw/qrM621GTGOKhGwTjKruA+7zUz4f5xqIMaFRvho0vw5Wf5S+3GbENKbIyDbBiMiM7Nar6rWBDceYPLjqBVj3GaSeOV+maRDT2BnDLKYxlKvqXLvxZ1xX2LM6c7nVgIwJiJyayC7CuWlxMrCEHMYfM6ZQVYyDdrfDjxMhzb3nNi0V5jx5fpuoyhDTyEk2MY2d5aqNnOf4TrB/Y/oElZcakCUoY7KVU4KJA64EBgO3ADOByT49uYwJrW4jYcUHToIJj4LfLYeU05C8BZK3us9b4OfvYdXU9PuWiz2fmLxUoUYL2DjbudZTpjyUqQAR3uXyThISsQRlTA5yugaTCswGZotIJE6iWSAiz6nqm4URoDHZ8s6Iufx957lSvFMe0yjztmdOwMHtPslnK2z+Eo7vO79N2lmY+VjmfX2FhUNEeSehpZ5Nvy4t1UkyP/zDaZ4rF3P+uWxViKx4vsmuoAnKmCIupxoMbmLpi5Nc6gNjgenBDcuYPOg2Evavh25PZL9dmXJO7aSGT+/6o3vg9TaQcspJGMNmO7WVM8echHTmuLN81rvs8zh7HHZ8Bwd3cG5w8bBwWPov51qQP2ERbtKJcWpGGWtQEpbz+zCmmMjpIv8koCXwBfCcz139xhQd3hkx87uvbw2oVru87Z8xQT28CsrHwunDcOIAnEh2Hz7LJw+cf12morOtV+OeULFG/t6LMUVMTjWY24DjwMPACJ9R7wVQVY0OYmzGFI7c1oD8yZigvMmhbBXn4a+pzpdvgkJgw+fOKAU9nnGa04wpxnK6BpPTWGXGFH8FqQFB4BJUu9ud5rkl42DDF3DN352RpI0ppiyBGFNQ3gSV36atbiOhbhe4/Cm4+kUY/iVEVoAPBsAndzsTrhlTDFmCMSbUMiaoOp3g3oXQbRSsnQ5vdYTVH9s8OKbYsQRjTFEUHgmXP+kkmir14ZPhMHkQHN6V467GFBWWYIwpymo0h+FzoddfYPtCeKszLHsP0rKdLcOYIsESjDFFXZgHLnoQ7v8eard3bgSd2A9+2xLqyIzJVo43WhpjioiqDeCOz5yhceb8Ad5M4NwNnr5yO9SMDVVjgswSjDHFiQi0u825IfO9XnBoR/r1ngio3cHpEJDVKNJeNpaaCTJLMMYURxXjYPgceK2lM36aV+pZWD4BfpzkjJdWplz6gTq9yxHlnNul01LTH1fVGc9t+UQnWYVFgCfcfY5whsLxlHGWKzeAfRvSn9/GUjM+LMEYU1xVjIP2d8BP/3ZqIeJxujg3ucodK+1E+mfvuGrH9jnjqJ05TqYZONLOwrzn8h9T6lk4fQQWj3M6KFRvAeVj/G9rNaASzxKMMcWZd7oCcGoVN03M2w2fGcdSu2chlK3sJJpU9+FdTkvJ/HrJv2D7AqcmJGFQvgZsmZd+ptEKcW6yae4MNFq9OcReaKNJlwJBSzAiMh7oB+xT1ZZuWVVgKs6ozDuAm1X1oIh0Bz4Dtru7/1dVn3f36Q28DniAd1X1Rbe8ATAFiAGWA7er6hl39OdJQAcgGRioqjuC9T6NCamsxkLL7/7VL8jb/jXbOAnKO03Bvd9AhepwbC/sXQv71sHedbB3DSx9B1JPO/uJB6rUs9GkS7hg1mAmAG/i/LH3GgXMU9UXRWSU+9r7bVqkqv18DyAiHuAtnEnPkoBlIjJDVdcBY4DXVHWKiIwDhgP/dJ8PqmpjERnkbjcwWG/SmJAryFhoBd0/qwRXMc55NO5xftvUFDiwDfatdZPOWqe57swxZ31YeP6SpCmyRIM4/ISI1Ac+96nBbAS6q+puEakJLFDVC9wazON+EsxFwGhV7eW+9s6F+yKwH4hT1RTf7URkjrv8g4iEA3uAWM3hjSYkJGhiYmKA3rkxpcjRPfDxnTBgQt6Tw9E98HprZxZSBB5ZBZXrBiNKEyQislxVE/ytK+wbLWuo6m53eQ/g+228SERWisgsEfHOCFUb2OmzTZJbFgMcUtWUDOXp9nHXH3a3z0RE7hGRRBFJ3L9/fwHfmjGlVEEG+6wYB21vw50BJPO01qZYC9md/G6Nwlur+BGop6ptgDeATwsphrdVNUFVE2JjYwvjlMaYjLqNhHoXQdNesPBld4ZQUxIUdoLZ6zaN4T7vA1DVI6p6zF3+AogQkWrALqCOz/7xblkyUNltAvMtx3cfd30ld3tjTFHkrQH1fc25+D9rVKgjMgFS2AlmBjDEXR6C03MMEYkTd7pMEenkxpUMLAOaiEgDESkDDAJmuLWf+cCAjMfKcI4BwNc5XX8xxhQBlWpD91GwaZYz4Zop9oKWYERkMvADcIGIJInIcJyL81eKyGagp/sanESwRkRWAmOBQepIAR4C5gDrgWmqutbd5wngMRHZgnON5T23/D0gxi1/DKenmjGmOOhyP8Q2g1lPwJkToY7GFFBQe5EVJ9aLzJgiYsd3MKEPXPp76PFMqKMxOShKvciMMSZ79S+BNoPhu7Hw2+ZQR2MKwBKMMaboufJ5Z0DOmb+3qaKLMUswxpiip0J16PFH2P4NrP1vqKMx+WQJxhhTNCUMg5ptYfYf4NSRUEdj8sESjDGmaArzQN9XnYEzF7yY8/amyLEEY4wpuuI7QIehsGQc7FkT6mhMHlmCMcYUbT2eceaomfkYpKWFOhqTB5ZgjDFFW7mqTq+ynUtg5YehjsbkgSUYY0zR1+YWqNMF5j4DJw6EOhqTS5ZgjDFFX1gY9H0FTh6Cec+HOhqTS5ZgjDHFQ1xL6HwfLJ8ASctDHY3JBUswxpjio/soZ3j/mY9CWmqoozE5sARjjCk+oqKh159h90pIHB/qaEwOLMEYY4qXFjdAw+4w709wbF+oozHZsARjjCleRKDPy3D2BHz5x1BHY7JhCcYYU/xUawKXPAyrpsCOb0MdjclCeM6bGGNMEbRplvM8oW/68rhWcJ8lnaLAEowxpniq0wX2b4S0lPNlnjIQ3yl3+4/rCntWZy63BBUwlmCMMcVTt5Gw4oP0CSYtBY7uhTlPOd2ZK9Z0nivEOc+RFc5vG9/JSVCpZ86X5SVBmRwFNcGIyHigH7BPVVu6ZVWBqUB9YAdws6oeFBEBXgf6ACeAoar6o7vPEOBp97AvqOpEt7wDMAEoC3wBPKyqmtU5gvlejTGFrGIctL0VfpoEqWeBMCgfC/vXwdZ5kHIy8z5lKrqJJw6iKmW+l0bCoNsThRJ+aRDsGswE4E1gkk/ZKGCeqr4oIqPc108AVwNN3Edn4J9AZzdZPAskAAosF5EZbsL4J3A3sAQnwfQGZmVzDmNMSeKtxXAWwsvAvYugYg1nmuXTR+DoHji626nVHN19/vWxvbA3w/D/YeFOwqpYIyRvpSQKaoJR1YUiUj9D8XVAd3d5IrAA54//dcAkVVVgsYhUFpGa7rZzVfUAgIjMBXqLyAIgWlUXu+WTgP44CSarcxhjShJvLWb5++mTg4hTQ4mqBLEXZL3/kd0wtg2knHZqM10fLZy4S4lQdFOuoaq73eU9gPffhdrATp/tktyy7MqT/JRnd450ROQeEUkUkcT9+/fn8+0YY0Kq20io2yV/TVvRNaHtbYAAChs+D3R0pVpI74NxaysaqnOo6tuqmqCqCbGxscEMwxgTLBXj4M5Z+W/a6jYS6l0EdS92pma26QACJhQJZq/b9IX77B3rYRdQx2e7eLcsu/J4P+XZncMYY9LzJqi+rzjXbb4ZE+qISoxQJJgZwBB3eQjwmU/5HeLoAhx2m7nmAFeJSBURqQJcBcxx1x0RkS5uD7Q7MhzL3zmMMca/Gs2hw1BY9i7s3xTqaEqEoCYYEZkM/ABcICJJIjIceBG4UkQ2Az3d1+D0AtsGbAHeAR4AcC/u/wlY5j6e917wd7d5191nK84FfrI5hzHGZO3ypyCiHHz5dM7bmhyJc4nCJCQkaGJiYqjDMMaE2nevO1Mz3z4dGl0R6miKPBFZrqoJ/tbZYJfGGOOr831Qpb4zGkBqSo6bm6xZgjHGGF/hkXDl87BvnTNKgMk3SzDGGJNRs2uh3iXw9Z/h1OFQR1NsWYIxxpiMRJypmU8kw8KXQx1NsWUJxhhj/KnVDtreAkvGwYHtoY6mWLIEY4wxWbnijxAW4fQqM3lmCcYYY7ISXRO6PgLrZ8CO70IdTbFjCcYYY7Jz0UMQXRvmPAlpaaGOplixBGOMMdkpUw56Pge7V8LKyaGOplixKZOzcfbsWZKSkjh16lSoQykxoqKiiI+PJyIiItShGJN7rQY4F/vnPQ/Nr0s/9bLJkiWYbCQlJVGxYkXq16+PM56mKQhVJTk5maSkJBo0aBDqcIzJPRHo/Vd470pnKJkrngp1RMWCNZFl49SpU8TExFhyCRARISYmxmqEpniq0wla3gjfj4VDO3Pe3liCyYkll8Cyz9MUaz1HO8/zngtpGMWFJRhjjMmtynWdXmWrP4IkG309J3YNJkD6vL6IdbuPZCpvXjOaLx6+NAQRGWOCouuj8NO/YfaTMPxL5/qM8ctqMAHSvm5lIjzpv2gRHqF9vSoFPvaePXsYNGgQjRo1okOHDvTp04ewsDA2btyYbrtHHnmEMWP8T/e6YMECRIR33333XNmKFSsQEV5+2RlrafHixXTu3Jm2bdvSrFkzRo8eDcCECROIjY2lbdu25x7r1q0r8PsypliKrODc4Z+0FNZ8EupoijSrweTSc/9by7pfM9dQvM6kpJGSln7ytpQ0Ze2uwwz81w9+92leK5pnr2mR7XlVleuvv54hQ4YwZcoUAFauXMnJkyeZMmUKzz77LABpaWl8/PHHfPdd1ncbt2zZkmnTpnHXXXcBMHnyZNq0aXNu/ZAhQ5g2bRpt2rQhNTU1XQIbOHAgb775ZraxGlNqtL0Flv4LvhoNF/aFiLKhjqhIshpMgJQJDyO2QiTeOowAsRUiKRNesI94/vz5REREcN99950ra9OmDWPHjmXq1KnnyhYuXEi9evWoV69elseqV68ep06dYu/evagqs2fP5uqrrz63ft++fdSsWRMAj8dD8+bNCxS7MSVWmAd6/RUO74Qf7B+vrFgNJpdyqmkA7Dtyiktfms/plDQiw8P4fERXqleMKtB516xZQ4cOHTKVt2rVirCwMFauXEmbNm2YMmUKgwcPzvF4AwYM4KOPPqJdu3a0b9+eyMjIc+seffRRLrjgArp3707v3r0ZMmQIUVFO/FOnTuXbb789t+0PP/xA2bL2X5spxeY86Tx//YLz8IprBfd963+fUiYkNRgReVhE1ojIWhF5xC0bLSK7RGSF++jjs/2TIrJFRDaKSC+f8t5u2RYRGeVT3kBElrjlU0WkTGG8r+rRUdzUIR4RGJBQp8DJJSeDBw9mypQppKSk8Omnn3LTTTfluM/NN9/MRx99xOTJkzMlpGeeeYbExESuuuoqPvzwQ3r37n1u3cCBA1mxYsW5hyUXU+rFd3JGWvblKeOUGyAECUZEWgJ3A52ANkA/EWnsrn5NVdu6jy/c7ZsDg4AWQG/gHyLiEREP8BZwNdAcGOxuCzDGPVZj4CAwvJDeHiN6NKFj/aqM6NE4541zoUWLFixfvtzvukGDBjFt2jS++uorWrduTY0aNXI8XlxcHBEREcydO5cePXpkWt+oUSPuv/9+5s2bx8qVK0lOTi7wezCmROo20mkq8yVh0O2J0MRTBIWiBtMMWKKqJ1Q1BfgGuCGb7a8DpqjqaVXdDmzBSU6dgC2quk1VzwBTgOvEuZPvCuBjd/+JQP/gvJXMqkdHMe3eiwJWe7niiis4ffo0b7/99rmyVatWsWjRIho1akS1atUYNWpUrprHvJ5//nnGjBmDx5P+l2PmzJmoOh0VNm/ejMfjoXLlygF5H8aUOBXjoO2tTq3FSxU2f2mjLrtCkWDWAJeKSIyIlAP6AHXcdQ+JyCoRGS8i3v69tQHfcRmS3LKsymOAQ27y8i3PRETuEZFEEUncv39/IN5bwIkI06dP56uvvqJRo0a0aNGCJ598kri4OMBpJtuwYQM33JBdjk7v4osvpn///pnK//3vf3PBBRfQtm1bbr/9dj744INzSWjq1Knpuil///33AXl/xhRr3UY6tRYATyTUaAEzHoLxV8GvK0IaWlEg3v9YC/WkIsOBB4DjwFrgNPBX4DdAgT8BNVV1mIi8CSxW1f+4+74HzHIP1VtV73LLbwc6A6Pd7Ru75XWAWaraMruYEhISNDEx/Z2569evp1mzZgV/wyYd+1xNifL5Y7D8fehwJ/R9BVZOgbl/hOO/QcIwuOJpKFc11FEGjYgsV9UEf+tCcpFfVd9T1Q6qehnONZJNqrpXVVNVNQ14B6cJDGAX52s4APFuWVblyUBlEQnPUG6MMYHXbSTU7eJcexGBtoPhoUTofK+TeN7oAMsnlspms1D1IqvuPtfFuf7yoYjU9NnkepymNIAZwCARiRSRBkATYCmwDGji9hgrg9MRYIY6VbL5wAB3/yHAZ8F+T0XFnDlz0jVltW3bluuvvz7UYRlTclWMgztnQUWfTjZlK8PVY+DeRRB7AfxvBLzXE3b9GLIwQyFU98F8IiIxwFngQVU9JCJviEhbnCayHcC9AKq6VkSmAeuAFHf7VAAReQiYA3iA8aq61j3+E8AUEXkB+Al4r9DeWYj16tWLXr165byhMSb44lo6yWfVNKfZ7J0roMMQ6PFsiW428wrJNZiiyK7BFB77XE2pdOoILHjRmRkzKhoiysERP633xexGzSJ3DcYYY0qdqGjo/RcneVRv7j+5lLAbNS3BGGNMYarRHIbOhD6vZF5Xwm7UtLHIAmVcV9izOnN5MavuGmMKgQh0ugt2r4AVH4C6PcxiLyxRIzNbDSZQ4julv6MXClzdTU5OPtcTLC4ujtq1a597febMmWz3TUxMZMSIEfk6r4hw2223nXudkpJCbGws/fr1A2Dv3r3069ePNm3a0Lx5c/r0cYaN27FjB2XLlk3Xg23SpEn5isGYUuGKp8//3ZAwJ+GMbQuLx0HK6VBGFhBWg8mtWaP811C8Us5AWkr6srQUZ5/3+/rfJ64VXP1iloeMiYlhxYoVAIwePZoKFSrw+OOPnz9lSgrh4f5/hAkJCSQk+L3ulqPy5cuzZs0aTp48SdmyZZk7dy61a58fDOGZZ57hyiuv5OGHHwacoWu8GjVqdC5mY0wOvMPNeG/UbHcbfPUszH4CFv/Dmdis5Y0QVjzrAsUz6qIovAyUrw6+M8KUr565VlNAQ4cO5b777qNz586MHDmSpUuXctFFF9GuXTsuvvjic5OELViw4FyNY/To0QwbNozu3bvTsGFDxo4dm+N5+vTpw8yZMwEyjby8e/du4uPjz71u3bp1IN+iMaWL742atdvDHTPgtv86nQL+exe83Q22fh3qKPPFajC5lU1N45yje+D1NpByCsIj4d6F6W++CpCkpCS+//57PB4PR44cYdGiRYSHh/PVV1/xhz/8gU8+yTyN64YNG5g/fz5Hjx7lggsu4P777yciIsLP0R2DBg3i+eefp1+/fqxatYphw4axaNEiAB588MFzM1z27NmTO++8k1q1agGwdetW2rZte+44b7zxBpdeemlgPwBjShLvjZpeItC4BzS8HNZ8DF//Cf59PTTsDj1HQ612oYo0zyzBBJJvdbftrUFJLgA33XTTuUEoDx8+zJAhQ9i8eTMiwtmzZ/3u07dvXyIjI4mMjKR69ers3bs3XS0ko9atW7Njxw4mT5587hqLV69evdi2bRuzZ89m1qxZtGvXjjVrnIEXrInMmAAJC4PWN0Pz6yBxPHzzErzdHVoOcK7dVG0Q6ghzZAkm0LqNhP3rg9rVsHz58ueW//jHP3L55Zczffp0duzYQffu3f3u4ztzpcfjISUlxe92vq699loef/xxFixYkGlemKpVq3LLLbdwyy230K9fPxYuXOh35k1jTAGFR0KX+6HtLfDdWPjhLVj3mdOEdsLPfE257blaCD1f7RpMoPkblyiIDh8+fO4C/IQJEwJ67GHDhvHss8/SqlWrdOVff/01J06cAODo0aNs3bqVunXrBvTcxpgMoipBjz/CiJ+czgD+kkteeq4GoedrRlaDKeZGjhzJkCFDeOGFF+jbN4veavkUHx/vt6vz8uXLeeihhwgPDyctLY277rqLjh07smPHjkzXYIYNG5bv7tLGGD+ia8I1f4dWN8PEvufvoQFIPQsbZ8OWuZzrcCTiLItPByQRp5draoYm9QDf6GljkblsLLLCY5+rMQHy+WPw40T3FokwiGkE8QnOzJoAqLvsvs64/OtPcHCHUxYWAe3vgH6v5imE7MYisxqMMcYUV91GOiMBpKU4t0oMnZm35nnfnq9hnoBfO7YEU0olJyfTo0ePTOXz5s0jJiYmBBEZY/KsoD1Xg9zz1RJMDlQVOdd2WXL4jhJQmKxJ1pgAK2jP1SD2fLUEk42oqCiSk5OJiYkpkUmmsKkqycnJREVFhToUY0qOjDdqFvb+2bAEk434+HiSkpLYv39/qEMpMaKiorK9wdMYU3JYgslGREQEDRoU/btljTGmKLIbLY0xxgSFJRhjjDFBYQnGGGNMUNid/C4R2Q/8XIBDVAN+C1A4wWDxFYzFVzAWX8EU5fjqqWqsvxWWYAJERBKzGi6hKLD4CsbiKxiLr2CKenxZsSYyY4wxQWEJxhhjTFBYggmct0MdQA4svoKx+ArG4iuYoh6fX3YNxhhjTFBYDcYYY0xQWIIxxhgTFJZg8khEeovIRhHZIiKj/KyPFJGp7volIlK/EGOrIyLzRWSdiKwVkYf9bNNdRA6LyAr38Uxhxeeef4eIrHbPnehnvYjIWPfzWyUi7Qsxtgt8PpcVInJERB7JsE2hfn4iMl5E9onIGp+yqiIyV0Q2u89Vsth3iLvNZhEZUojx/U1ENrg/v+kiUjmLfbP9LgQxvtEissvnZ9gni32z/V0PYnxTfWLbISIrstg36J9fgamqPXL5ADzAVqAhUAZYCTTPsM0DwDh3eRAwtRDjqwm0d5crApv8xNcd+DyEn+EOoFo26/sAs3AmFO8CLAnhz3oPzk1kIfv8gMuA9sAan7KXgFHu8ihgjJ/9qgLb3Ocq7nKVQorvKiDcXR7jL77cfBeCGN9o4PFc/Pyz/V0PVnwZ1r8CPBOqz6+gD6vB5E0nYIuqblPVM8AU4LoM21wHTHSXPwZ6SCFNJqOqu1X1R3f5KLAeqF0Y5w6g64BJ6lgMVBaRmiGIowewVVULMrpDganqQuBAhmLf79hEoL+fXXsBc1X1gKoeBOYCvQsjPlX9UlVT3JeLgZDNz5DF55cbufldL7Ds4nP/btwMTA70eQuLJZi8qQ3s9HmdROY/4Oe2cX/JDgOFPgex2zTXDljiZ/VFIrJSRGaJSIvCjQwFvhSR5SJyj5/1ufmMC8Mgsv7FDuXnB1BDVXe7y3sAf/PcFpXPcRhOjdSfnL4LwfSQ24Q3PosmxqLw+V0K7FXVzVmsD+XnlyuWYEogEakAfAI8oqpHMqz+EafZpw3wBvBpIYfXVVXbA1cDD4rIZYV8/hyJSBngWuAjP6tD/fmlo05bSZG810BEngJSgA+y2CRU34V/Ao2AtsBunGaoomgw2ddeivzvkiWYvNkF1PF5He+W+d1GRMKBSkByoUTnnDMCJ7l8oKr/zbheVY+o6jF3+QsgQkSqFVZ8qrrLfd4HTMdpivCVm8842K4GflTVvRlXhPrzc+31Nhu6z/v8bBPSz1FEhgL9gFvdJJhJLr4LQaGqe1U1VVXTgHeyOG+oP79w4AZgalbbhOrzywtLMHmzDGgiIg3c/3IHATMybDMD8PbYGQB8ndUvWKC5bbbvAetV9dUstonzXhMSkU4434FCSYAiUl5EKnqXcS4Gr8mw2QzgDrc3WRfgsE9zUGHJ8j/HUH5+Pny/Y0OAz/xsMwe4SkSquE1AV7llQScivYGRwLWqeiKLbXLzXQhWfL7X9K7P4ry5+V0Ppp7ABlVN8rcylJ9fnoS6l0Fxe+D0ctqE08PkKbfseZxfJoAonKaVLcBSoGEhxtYVp7lkFbDCffQB7gPuc7d5CFiL0ytmMXBxIcbX0D3vSjcG7+fnG58Ab7mf72ogoZB/vuVxEkYln7KQfX44iW43cBbnOsBwnGt684DNwFdAVXfbBOBdn32Hud/DLcCdhRjfFpzrF97voLdXZS3gi+y+C4UU37/d79YqnKRRM2N87utMv+uFEZ9bPsH7nfPZttA/v4I+bKgYY4wxQWFNZMYYY4LCEowxxpigsARjjDEmKCzBGGOMCQpLMMYYY4LCEkwJJiIqIq/4vH5cREYH6NgTRGRAII6Vw3luEpH1IjK/AMd4V0Sa53Pf7wtw3gUikpDf/YsLEemfl89XRBJEZGwwYyoMBflulBaWYEq208ANIbjTPFvuXcq5NRy4W1Uvz+e5PKp6l6quy8/+qnpxfvYrZfoDuU4wqpqoqiPyezIR8eR3X3f/vHz/smTfjZxZginZUnDm8n4044qMNRAROeY+dxeRb0TkMxHZJiIvisitIrLUnXuikc9heopIoohsEpF+7v4eceYDWeYOJnivz3EXicgMINMfexEZ7B5/jYiMccuewbl59D0R+VuG7buLyEIRmSnOnB3jRCTM+15E5BURWYkzMOW5moS77s/uYJWLRaSGW15DnLlLVrqPi/18Llmd75/u57BWRJ7L6YciIh1F5Hv3PEtFpKKIRInI++5n8JOIXO5uO1REPhVn3pcdIvKQiDzmbrNYRKq62y0QkdfFmRtkjTijDHjnjvnU/VksFpHWbvlocQZ6XOD+nEf4xHebG9cKEfmX9w+6v8/O/ZyuBf7mbt9IREaIMyfRKhGZ4uf9dxeRz3OKI8M+GX+mWcU43P0+LhWRd0TkTbd8gvszWwK85MY5W5yBIheJyIXudje5n99KEVnolrXwOdcqEWmS4bsh4nzn17g/v4E+73OBiHwszvw4H4gUzsjqRUao7/S0R/AewDEgGmfeiErA48Bod90EYIDvtu5zd+AQztwykTjjLz3nrnsY+LvP/rNx/klpgnMXchRwD/C0u00kkAg0cI97HGjgJ85awC9ALBAOfA30d9ctwM/d/O7xTuHc0ezBGY5+gLtOgZt9tj13DHfdNe7ySz6xTsUZHBT3eJX8fC5Zna+qz34LgNZZxY4zt8g2oKP7Otp9z78HxrtlF7qfRxQwFOfO+Iru53OY86MKvOYT8wLgHXf5Mtz5RXAG5HzWXb4CWOEujwa+d39G1XBGL4gAmgH/AyLc7f4B3JHDZzeB9N+lX4FId7lyFj+7z7OLw88+536mWcWI8z3agTMHTgSwCHjTJ8bPAY/7eh7QxF3ujDOkEzh3+Nf2jd39DG/1+fmVzfDduBHn++DBGdn6F5zfn+7uzyse5/fkB5wBKkP+t6GwHlaDKeHUGU15EpCXJoll6swtcxpnmIwv3fLVQH2f7aapapo6w4lvw/nDeBXOWGIrcKYKiMFJQABLVXW7n/N1BBao6n51pjj4AOePZE6WqjNfRyrOkBtd3fJUnAE//TmD84cGYLnP+7kCZ5Rd1BkI8XAezneziPwI/AS0IPvmoguA3aq6zD3XEfc9dwX+45ZtAH4Gmrr7zFfVo6q6H+cP1v/c8ow/j8nu/guBaHFmkuyKMzQKqvo1ECMi0e72M1X1tKr+hjNgZg2ceXA6AMvcn2EPnKSa3WeX0SrgAxG5DacWnRN/cWTk+zPNKsZOwDfqzIFzlsyjYX+kqqnijDZ+MfCRu/+/cBICwHfABBG5GydhgJMY/iAiT+CMpH0yw3G7ApPd781e4Buc7zQ435kkdQbWXEHWn1mJFJC2SFPk/R1nmPn3fcpScJtI3aaeMj7rTvssp/m8TiP9dybjOEOKM5bY71Q13cCKItIdpwYTSP7OD3DKTQL+nFX3306cP1p5+R3IdD4RaYBTM+yoqgdFZAJOzSOQCvLzyO1xvZ+FABNV9Uk/2+f2s+uL8w/CNcBTItJKz09Alts4MvL9mfqNUUT6Z3MOOP/9CwMOqWrbjBuo6n0i0tl9D8tFpIOqfug2rfUFvhCRe91knRu5eW8lltVgSgFVPQBMw7lg7rUD579AcNrQI/Jx6JtEJEyc6zINgY04I/beL860AYhIU3FGe83OUqCbiFRz29IH4/wXmJNO4ox2GwYMBL7Nx3vwmgfc78bsEZFKuTxfNM4frsPiXM+5OofzbARqikhH91wVxbnovAi41S1rCtR1t80Lb9t/V5xRqA9nOG534DfNPEeQr3nAABGp7u5TVUTq5XDeozhNeN5/Vuqo6nzgCZym2Qp5fB85ySrGZTjfoyruZ3qjv53d979dRG5y9xcRaeMuN1LVJar6DLAfqCMiDYFtqjoWZ+Tq1hkOuQgY6H5vYnGS69IAv+diyRJM6fEKThu31zs4v4wrgYvIX+3iF5xfpFk41wVOAe/iXMT/UUTW4DQ/ZPtfmzrD8Y8C5uOMDrtcVf0NQZ/RMuBNnKmht+PMiZFfDwOXi8hqnOYff81cmc6nqitxmsY2AB/iNLFkSZ3pdwcCb7if/VycGs8/gDD3/FOBoW4TZV6cEpGfgHGc/2diNNBBRFYBL3J+mP+s4lsHPI0zU+IqN76cpqyeAvyfe+4mwH/c9/ETMFZVD+XxfWQrqxjVmR/lLzjfye9w/ony19QJTtId7v4M1nJ+OuS/uRfq1+BcG1qJM23xGrc5rSVOk7Ov6TjNgitxrh+OVNU9AXirxZ6NpmyKJfe/8cdVtV9JPF9eicgCnPgSQx1LKIlIBVU95tZgpuN0nCjIPx6mAKwGY4wpSUa7NY01OLXMT0MaTSlnNRhjjDFBYTUYY4wxQWEJxhhjTFBYgjHGGBMUlmCMMcYEhSUYY4wxQfH/OudHUEETeLYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mse_train = []\n",
    "mse_cv = []\n",
    "\n",
    "# Calculate MSE using CV for the 19 principle components, adding one component at the time.\n",
    "for i in np.arange(1, 20):\n",
    "    results =  pcr(X,y,i)\n",
    "    mse_train.append(results[3])\n",
    "    mse_cv.append(results[4])\n",
    "    \n",
    "# Plot results    \n",
    "plt.plot(np.arange(1, 20), mse_cv, '-v', label = 'CV_MSE')\n",
    "plt.plot(np.arange(1, 20), mse_train, '-v', label = 'Train_MSE')\n",
    "plt.xlabel('Number of principal components in regression')\n",
    "plt.ylabel('MSE')\n",
    "plt.title('Salary')\n",
    "plt.xlim(xmin=-1)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the CV MSE starts rising post 6 components and falls again at multiple points later, while train MSE keep falling. 6 components explain 92% of the variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "133288.1611562523\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.66)\n",
    "\n",
    "pca2 = PCA()\n",
    "\n",
    "# Scale the data\n",
    "X_reduced_train = pca2.fit_transform(scale(X_train))\n",
    "X_reduced_test = pca2.transform(scale(X_test))[:,:6]\n",
    "\n",
    "# Train regression model on training data \n",
    "regr = linear_model.LinearRegression()\n",
    "a = regr.fit(X_reduced_train[:,:6], y_train)\n",
    "\n",
    "# Prediction with test data\n",
    "pred = regr.predict(X_reduced_test)\n",
    "print(mean_squared_error(y_test, pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial Least Squares\n",
    "\n",
    "Like PCR, PLS is a dimension reduction method, which first identifies a new set of features Z1,...,ZM that are linear combinations of the original features, and then fits a linear model via least squares using these M new features. But unlike PCR, PLS identifies these new features in a supervised way—that is, it makes use of the response Y in order to identify new features that not only approximate the old features well, but also that are related to the response. Roughly speaking, the PLS approach attempts to find directions that help explain both the response and the predictors.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "119839.95920752938"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train, X_test , y_train, y_test = train_test_split(X, y, test_size=0.66)\n",
    "pls = PLSRegression(n_components=6)\n",
    "pls.fit(scale(X_train), y_train)\n",
    "\n",
    "mean_squared_error(y_test, pls.predict(scale(X_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}