Principal components regression

PCR-Examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demo: Principal Components Regression\n",
    "                                                By Arin Ghosh, Karanjit Singh Tiwana, Manoj Karthick Selva Kumar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What is PCR ?**\n",
    "\n",
    "The idea behind principal components regression is to calculate an optimal number of principal components using a training dataset and use those components as features in a linear regression model to make predictions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Implementation of PCR**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "code_folding": [
     13
    ],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.base import RegressorMixin\n",
    "from sklearn.utils.validation import check_is_fitted\n",
    "from sklearn.utils.validation import check_random_state\n",
    "from sklearn.utils.validation import check_X_y\n",
    "import copy\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "\n",
    "class PCR(BaseEstimator, RegressorMixin):\n",
    "    \"\"\"\n",
    "    Principal Components Regression algorithm as presented in\n",
    "    Introduction to Statistical Learning by Tibshirani et al (2013).\n",
    "    The underlying model is a linear regression model operating on\n",
    "    reduced number of components general by PCA.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_components : int, float, None or string\n",
    "        Number of components to keep.\n",
    "        if n_components is not set all components are kept::\n",
    "            n_components == min(n_samples, n_features)\n",
    "\n",
    "    copy : bool (default True)\n",
    "        If False, data passed to fit are overwritten and running\n",
    "        fit(X).transform(X) will not yield the expected results,\n",
    "        use fit_transform(X) instead.\n",
    "\n",
    "    whiten : bool, optional (default False)\n",
    "        When True (False by default) the `components_` vectors are multiplied\n",
    "        by the square root of n_samples and then divided by the singular values\n",
    "        to ensure uncorrelated outputs with unit component-wise variances.\n",
    "        Whitening will remove some information from the transformed signal\n",
    "        (the relative variance scales of the components) but can sometime\n",
    "        improve the predictive accuracy of the downstream estimators by\n",
    "        making their data respect some hard-wired assumptions.\n",
    "\n",
    "    svd_solver : string {'auto', 'full', 'arpack', 'randomized'}\n",
    "        auto :\n",
    "            the solver is selected by a default policy based on `X.shape` and\n",
    "            `n_components`: if the input data is larger than 500x500 and the\n",
    "            number of components to extract is lower than 80% of the smallest\n",
    "            dimension of the data, then the more efficient 'randomized'\n",
    "            method is enabled. Otherwise the exact full SVD is computed and\n",
    "            optionally truncated afterwards.\n",
    "        full :\n",
    "            run exact full SVD calling the standard LAPACK solver via\n",
    "            `scipy.linalg.svd` and select the components by postprocessing\n",
    "        arpack :\n",
    "            run SVD truncated to n_components calling ARPACK solver via\n",
    "            `scipy.sparse.linalg.svds`. It requires strictly\n",
    "            0 < n_components < min(X.shape)\n",
    "        randomized :\n",
    "            run randomized SVD by the method of Halko et al.\n",
    "\n",
    "    tol : float >= 0, optional (default .0)\n",
    "        Tolerance for singular values computed by svd_solver == 'arpack'.\n",
    "\n",
    "    iterated_power : int >= 0, or 'auto', (default 'auto')\n",
    "        Number of iterations for the power method computed by\n",
    "        svd_solver == 'randomized'.\n",
    "\n",
    "    random_state : int, RandomState instance or None, optional (default None)\n",
    "        If int, random_state is the seed used by the random number generator;\n",
    "        If RandomState instance, random_state is the random number generator;\n",
    "        If None, the random number generator is the RandomState instance used\n",
    "        by `np.random`. Used when ``svd_solver`` == 'arpack' or 'randomized'.\n",
    "\n",
    "    fit_intercept : boolean, optional\n",
    "        whether to calculate the intercept for this model. If set\n",
    "        to false, no intercept will be used in calculations\n",
    "        (e.g. data is expected to be already centered).\n",
    "\n",
    "    copy_X : boolean, optional, default True\n",
    "        If True, X will be copied; else, it may be overwritten.\n",
    "\n",
    "    normalize : boolean, optional, default False\n",
    "        If True, the regressors X will be normalized before regression.\n",
    "\n",
    "    n_jobs : int, optional, default 1\n",
    "        The number of jobs to use for the computation.\n",
    "        If -1 all CPUs are used. This will only provide speedup for\n",
    "        n_targets > 1 and sufficient large problems.\n",
    "\n",
    "    Attributes\n",
    "    ----------\n",
    "    components_ : array, shape (n_components, n_features)\n",
    "        Principal axes in feature space, representing the directions of\n",
    "        maximum variance in the data. The components are sorted by\n",
    "        ``explained_variance_``.\n",
    "    explained_variance_ : array, shape (n_components,)\n",
    "        The amount of variance explained by each of the selected components.\n",
    "        Equal to n_components largest eigenvalues\n",
    "        of the covariance matrix of X.\n",
    "        .. versionadded:: 0.18\n",
    "    explained_variance_ratio_ : array, shape (n_components,)\n",
    "        Percentage of variance explained by each of the selected components.\n",
    "        If ``n_components`` is not set then all components are stored and the\n",
    "        sum of the ratios is equal to 1.0.\n",
    "    singular_values_ : array, shape (n_components,)\n",
    "        The singular values corresponding to each of the selected components.\n",
    "        The singular values are equal to the 2-norms of the ``n_components``\n",
    "        variables in the lower-dimensional space.\n",
    "    mean_ : array, shape (n_features,)\n",
    "        Per-feature empirical mean, estimated from the training set.\n",
    "        Equal to `X.mean(axis=0)`.\n",
    "    n_components_ : int\n",
    "        The estimated number of components. When n_components is set\n",
    "        to 'mle' or a number between 0 and 1 \n",
    "        (with svd_solver == 'full') this number is estimated from \n",
    "        input data. Otherwise it equals the parameter n_components, \n",
    "        or the lesser value of n_features and n_samples if n_components\n",
    "        is None.\n",
    "    noise_variance_ : float\n",
    "        The estimated noise covariance following the Probabilistic PCA model\n",
    "        from Tipping and Bishop 1999. See \"Pattern Recognition and\n",
    "        Machine Learning\" by C. Bishop, 12.2.1 p. 574 or\n",
    "        http://www.miketipping.com/papers/met-mppca.pdf. It is required to\n",
    "        computed the estimated data covariance and score samples.\n",
    "        Equal to the average of (min(n_features, n_samples) - n_components)\n",
    "        smallest eigenvalues of the covariance matrix of X.\n",
    "    coef_ : array, shape (n_features, ) or (n_targets, n_features)\n",
    "        Estimated coefficients for the linear regression problem.\n",
    "        If multiple targets are passed during the fit (y 2D), this\n",
    "        is a 2D array of shape (n_targets, n_features), while if only\n",
    "        one target is passed, this is a 1D array of length n_features.\n",
    "    intercept_ : array\n",
    "        Independent term in the linear model.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, fit_intercept=True, copy_X=True, normalize=False,\n",
    "                 n_jobs=1, n_components=None, copy=True, whiten=False,\n",
    "                 svd_solver='auto', tol=0.0, iterated_power='auto',\n",
    "                 random_state=None):\n",
    "\n",
    "        self.fit_intercept = fit_intercept\n",
    "        self.copy_X = copy_X\n",
    "        self.normalize = normalize\n",
    "        self.n_jobs = n_jobs\n",
    "        self.n_components = n_components\n",
    "        self.copy = copy\n",
    "        self.whiten = whiten\n",
    "        self.svd_solver = svd_solver\n",
    "        self.tol = tol\n",
    "        self.iterated_power = iterated_power\n",
    "        self.random_state = random_state\n",
    "\n",
    "    def _get_dict(self, obj):\n",
    "        \"\"\"\n",
    "         Get the attributes associated with the PCA or Linear Regression model.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        obj : The PCA or regression model whose attributes needs to looked up.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        for attr, value in obj.__dict__.items():\n",
    "            self.__dict__[attr] = value\n",
    "\n",
    "    def _set_attributes(self, model):\n",
    "        \"\"\"\n",
    "         Set the attributes associated with the PCA or Linear Regression\n",
    "         model to passed estimator.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        model : Specifies the PCA or linear regression model whose attributes\n",
    "        to be set\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        if model:\n",
    "            self._get_dict(model)\n",
    "\n",
    "    def get_regression_model(self):\n",
    "        \"\"\"\n",
    "        The associated linear regression model for PCR.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        _lr: Returns the unfitted linear regression model of PCR.\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        return self._lr\n",
    "\n",
    "    def get_transformed_data(self):\n",
    "        \"\"\"\n",
    "        Generates the reduced principal components from the training data\n",
    "        using PCA.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        transformed: Returns a transformed numpy array or sparse matrix. The\n",
    "        data has been reduced to the specified number of components using the\n",
    "        `n_components` parameter of the PCR model.\n",
    "\n",
    "        Notes\n",
    "        -------\n",
    "        The algorithm should have first been executed on a training dataset.\n",
    "\n",
    "        \"\"\"\n",
    "\n",
    "        transformed = self._X_reduced\n",
    "        return transformed\n",
    "\n",
    "    def _reduce(self, X):\n",
    "        \"\"\"\n",
    "        Fit the Principal Components Analysis model.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : numpy array or sparse matrix of shape [n_samples,n_features]\n",
    "            Training data\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        \"\"\"\n",
    "        if not self.random_state:\n",
    "            random_state = 0\n",
    "        else:\n",
    "            random_state = self.random_state\n",
    "\n",
    "        self._pca = PCA(n_components=self.n_components, copy=self.copy,\n",
    "                        whiten=self.whiten, svd_solver=self.svd_solver,\n",
    "                        tol=self.tol, iterated_power=self.iterated_power,\n",
    "                        random_state=random_state)\n",
    "\n",
    "        self._pca.fit(X)\n",
    "        self._X_reduced = self._pca.transform(X)\n",
    "        self._set_attributes(self._pca)\n",
    "\n",
    "    def _regress(self, X, y):\n",
    "        \"\"\"\n",
    "        Fit the Linear Regression model.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : numpy array or sparse matrix of shape [n_samples,n_features]\n",
    "            Training data\n",
    "        y : numpy array of shape [n_samples, n_targets]\n",
    "            Target values\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        \"\"\"\n",
    "        self._model.fit(X, np.ravel(y))\n",
    "        self._set_attributes(self._model)\n",
    "\n",
    "    def fit(self, X, y):\n",
    "        \"\"\"\n",
    "        Fit the Principal Components Regression model.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : numpy array or sparse matrix of shape [n_samples,n_features]\n",
    "            Training data\n",
    "        y : numpy array of shape [n_samples, n_targets]\n",
    "            Target values\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        self : returns an instance of self.\n",
    "        \"\"\"\n",
    "        X, y = check_X_y(X, y)\n",
    "\n",
    "        self._model = LinearRegression(fit_intercept=self.fit_intercept,\n",
    "                                       copy_X=self.copy_X,\n",
    "                                       normalize=self.normalize,\n",
    "                                       n_jobs=self.n_jobs)\n",
    "\n",
    "        self._lr = copy.deepcopy(self._model)\n",
    "\n",
    "        self._reduce(X)\n",
    "        self._regress(self._X_reduced, y)\n",
    "\n",
    "        return self\n",
    "\n",
    "    def predict(self, X):\n",
    "        \"\"\"\n",
    "        Predict using the PCR model\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix}, shape = (n_samples, n_features)\n",
    "            Samples.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        predictions : array, shape = (n_samples,)\n",
    "            Returns predicted values.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, \"coef_\")\n",
    "        X_reduced_test = self._pca.transform(X)\n",
    "        predictions = self._model.predict(X_reduced_test)\n",
    "        return predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above class inherits the BaseEstimator and RegressorMixin components of the scikit-learn library to create a PCR Estimator. It provides the `fit()` and `predict()` methods for learning from data and predicting based on unseen data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example using the inbuilt diabetes dataset in sklearn** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean squared error on test data using just OLS regression: 2817.80157017\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn import datasets\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.cross_validation import train_test_split\n",
    "diabetes = datasets.load_diabetes()\n",
    "\n",
    "X = diabetes.data\n",
    "y = diabetes.target\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.33, random_state=42)\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train, y_train)\n",
    "p = lr.predict(X_test)\n",
    "e = mean_squared_error(p, y_test)\n",
    "print(\"Mean squared error on test data using just OLS regression:\", e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean squared error on the test data using PCR: 2855.17572618\n"
     ]
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "from sklearn.cross_validation import train_test_split\n",
    "\n",
    "diabetes = datasets.load_diabetes()\n",
    "\n",
    "X = diabetes.data\n",
    "y = diabetes.target\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.33, random_state=42)\n",
    "\n",
    "pcr = PCR(n_components=7)\n",
    "pcr.fit(X_train, y_train)\n",
    "preds = pcr.predict(X_test)\n",
    "error = mean_squared_error(preds, y_test)\n",
    "print(\"Mean squared error on the test data using PCR:\", error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Example using Hitters data from the book [\"Introduction to Statistical learning\"](http://www-bcf.usc.edu/~gareth/ISL/getbook.html) by Hastie, et al**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example uses k-Fold Cross validation with Principal Components Regression to show variation of MSE with different number of components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold, cross_val_score\n",
    "from sklearn.preprocessing import scale\n",
    "from sklearn import datasets\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.cross_validation import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "df = pd.read_csv(\"Hitters.csv\").dropna().drop(\"Player\", axis=1)\n",
    "dummies = pd.get_dummies(df[['League', 'Division', \"NewLeague\"]])\n",
    "\n",
    "y = df.Salary\n",
    "X_ = df.drop(['Salary', 'League', 'Division', 'NewLeague'],\n",
    "             axis=1).astype('float64')\n",
    "X = pd.concat([X_, dummies[['League_N', 'Division_W', 'NewLeague_N']]], axis=1)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.5, random_state=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The possible number of components are 0 to 19. The following program iterates by adding each of the principal components and calculates the MSE. The scikit-learn `KFold` and `cross_val_score` objects are used with our implementation of PCR to do the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEWCAYAAACwtjr+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XecVNX5+PHPsx220HaHjiBFpSgC\nKraIYhRLsHw1tiga8zMaFY0aS0zsJmoSW4wajcaOorGQxIYFNSrSkaYCitJZOguyy+4+vz/OmfWy\nzO7OzM7sbHner9e8dubMveeee/fOPHPvPfc5oqoYY4wxyZSW6gYYY4xp/izYGGOMSToLNsYYY5LO\ngo0xxpiks2BjjDEm6SzYGGOMSToLNiYhRERFpI9//rCI/D6aaeNYzlki8na87WzsRKSjiHwoIltE\n5C+pbo8xiWLBJoVE5EwRmSYiJSKyUkTeEJFDUt2u+lLVC1X11vrWIyI9fWDKCNT9rKoeVd+6G7EL\ngLVAgape2dALF+dOEVnnH3eJiNQw7QgRWVaPZZ3r/7+/qVa+TERGxFtvcyUiT4jIbaluR7ws2KSI\niFwB3Av8AegI9AAeBE6oYfqMSOWm4UX6X8Tz/xGR9AjFuwHztYa7rRtgP7gAOBHYB9gbOB74ZRKX\ntx64RkQKkrgM0xioqj0a+AG0AUqAU2uZ5ibgJeAZYDPwCyAbF6BW+Me9QLafvhD4D7AR9wH+CEjz\n710DLAe2AF8CIyMsbziwCkgPlJ0EfO6f7w986utfCTwAZAWmVaCPf/4EcFvgvd/4eVYAP6827XHA\nTL+OS4GbAvN956ct8Y8DgXOB/wWmOQiYCmzyfw8KvDcJuBX42K/720BhLdv8eGCWX8dPgL0D7y3x\n2/FzoBTIqKFsL7/cjcA8YHSgjieAh4DXga3AkdWW/wSwAyjz63tkHPvBCGAZcDWwxm/3E4Fjga/8\nvvHbWrbBJ8AFgdfnA5MjTJcLfA9UBv4/XWprW4Q6zgX+B/wbuDFQvgwYUcM8rYC/AN/6//n/gFb+\nvdF+m2/0/4O9qv3/fuP/V1uBx3A/8t7w+8Y7QDs/bU/cfneBX4eVwJWBuqLZ/lcGtv951eb9M27f\nXg08HGh/jfP6tgT3jX9H+9luLI+UN6AlPoBRQDmQUcs0N/md60TcEWgr4BZgMhACivwXw61++j/6\nHTfTPw4FBNgD9yXexU/XE+hdwzIXAz8OvH4RuNY/H4oLSBm+jgXA5YFpIwYbv66rgYG4L6jnqk07\nAhjk13FvP+2JgbZqcDsRCDZAe2ADcLZv1xn+dQf//iS/Tv389psE3FHDug/xH/ADgHRgDO4LKvwl\nsgQXiLrzw5fDTmV+uy8CfgtkAUf4L4E9AttlE3CwX9+cCO2o2nZx7gcjcPvWDb49/w8o9ts9HxgA\nbAd2r2E7bAIOCLweBmypYdoRwLJqZTW2LcL85+KCxWBcgGjvy2sLNn/z/8eu/v90EO4LvB8uiPzY\nr/fV/n+RFfhfTcYFmK7+fz0D2NfP/x4+4PHDfjcOt88O8tvwyLrWMbD9b/HtOBbYxg+B7F5gAm7f\nzccF2j9GOW/1fSPqz3ZjeKS8AS3xAZwFrKpjmpuAD6uVLQaODbw+Gljin98CvIb/Eg9M08d/sI4E\nMutY5m3A4/55vv/w7lbDtJcDrwRe1xRsHifwBe+/FLR6OwPv3wvc45+HP/Q1BZuzgSnV5v8UONc/\nnwT8LvDer4A3a1juQ1T7UsT9UjzMP18C/Lza+zuV4QL8KvwRpS8bhz9a89vlqTr+B9W/UGLdD0bg\njjjSA/9HZecAMh0f0CMsvwLYM/C6r59fIkw7gl2DTY1tizB/8H85HrjTP48YbHDB9ntgnwjv/R4Y\nX23a5eF6/P/qrMD7/wIeCry+FHi12n4X3A53AY/FsP2D++wa3A81wX2megfeOxD4pq55a9g3ov5s\nN4aHXbNJjXVAYRTn35dWe90Fd/og7FtfBvAn3C+5t0XkaxG5FkBVF+ECw03AGhF5XkS6ENlzwMki\nkg2cDMxQ1W8BRKSfiPxHRFaJyGbctabCKNa1S7X1CLYfETlARN4XkWIR2QRcGGW94bq/rVb2Le6X\na9iqwPNtQF4Nde0GXCkiG8MP3BFLcFtV/39UL+sCLFXVylraE6mOusSyHwCsU9UK//x7/3d14P3v\nqXk7lADB6ycFQIn6b7co1NW2mtwAXCQinWqZphDIwX3Z17pc/z9Yys7bvvo2qGubVN9vw+sRzfYv\nD7wO73dFQGtgemAfe9OX1zXvLmL8bKecBZvU+BR3KuPEOqar/gFfgftSDOvhy1DVLap6paruDvwE\nuEJERvr3nlPVQ/y8CtwZcWGq83EfnGOAM3HBJ+wh4Augr6oW4E4VReylVM1K3Jd2sM1Bz+FOK3RX\n1Ta4U4Hheuv6gqu+PcL1L4+iXdUtBW5X1baBR2tVHReYJlJ7gmUrgO4iEvxcVW9PtF/aNS0jvJyI\n+0ECzMN1Dgjbx5dF0y6Is22q+gXwMm6/qsla3Oemd13L9T3ouhPfvhBWfb8Nr0e8238tLqgNCOxj\nbVS1psBf3S7bO9rPdmNgwSYFVHUT7pfc30TkRBFpLSKZInKMiNxVy6zjgN+JSJGIFPo6ngEQkeNF\npI//kG3GnQ6pEJE9ROQIf7SyHbezV9RQP7gv/7HAj3DXbMLyfb0lIrIncFGUqzseOFdE+otIa+DG\nau/nA+tVdbuI7I8LcmHFuAvQu9dQ9+tAP9+FPENETgP64zpKxOpR4EJ/pCUikisix4lIfgx1fIY7\nTXK1/3+OwAX+5+NoT21q3A8S4CncD5Wu/lfylbjTN5GsBjqISJsEte1m4DygbaQ3/dHK48DdItJF\nRNJF5EC/b48HjhORkSKS6dtdirueEq/f+8/mAN+uF3x5XOvo2/8ocI+IhAD8dj46yvasJvBZiOOz\nnVIWbFJEVe8GrgB+h/tSXQpcArxay2y3AdNwPWrm4C5whvvd98X1qCnBHTk9qKqTcBc/78D9qlqF\nu6hZ26/Hcbhzx++p6tpA+VW4QLAF94F5YddZI67nG7jrMO/hTvO9V22SXwG3iMgW3Id2fGDebcDt\nwMf+tMPwanWvw/UguxJ3avJq4Phq7Y6Kqk7DXUx/ANfJYBHumkIsdZThekQdg9veDwLn+F/tiVTb\nflBff8ddtJ4DzAX+68t24ddrHPC1//90qU/bVPUb4GncRfmaXOXrnYrrWXcn7hrZl8DPgL/itv1P\ngJ/4/0m8PsDtB+8Cf1bV8M3E9dn+1/g6J/vT0e/gLvRH4zGgv9/WrxL7ZzulJPpTscYY0/yJSE/g\nG9xF9/LapzbRsiMbY4wxSWfBxhhjTNLZaTRjjDFJZ0c2xhhjks6SO3qFhYXas2fPVDfDGGOalOnT\np69V1aK6pktasBGR7rg++51w90o8oqr3Bd6/CnfXe5GqrvX3h9zHD/mAzlXVGX7aMbguwuDSNTzp\ny4fi7gFohbvn4jJVVRFpj+ua2xOXpuKnqrqhtvb27NmTadOmJWDNjTGm5RCR6lk8IkrmabRyXKbU\nvXB5gS4Wkf6+cd1xCfO+C0x/DO5ekb64DKcP+Wnb424EPACXefhGEWnn53nITxueb5QvvxZ4V1X7\n4vrIX5ukdTTGGBOFpAUbVV0ZPjJR1S24LMHhPEX34G7AC/ZOOAGXpFBVdTLQVkQ645LcTVTV9f7o\nZCIwyr9XoKqf+rxNT/FD+pcTgCf98yepOy2MMcaYJGqQDgL+Jql9gc9EZDSwXFVnV5usKzsnvlvm\ny2orXxahHKCjqq4EF/Rwd9ZGatcF4kbKnFZcXBzHmhljjIlG0oONiOTh0nlfjju1dj0uLckuk0Yo\n0zjKo6aqj6jqMFUdVlRU5/UtY4wxcUpqsPEJ8f4FPKuqL+OytfYCZovIEqAbMMOnFV/GzllWu+Ey\nqdZW3i1COcBqf5oN/3dNYtfMGGNMLJLZG01wieMW+KSTqOocAqe0fMAZ5nujTQAuEZHncZ0BNqnq\nShF5C/hDoFPAUcB1qrpeRLb45IyfAefgkvCBS1k/BpekbgxuULGEOva+j5i/cvMu5f07F/D6ZYcm\nenHGGNOkJfPI5mDcSIpHiMgs/zi2lulfB77GZUR9FJcNGFVdjxtHfqp/3OLLwKW5/4efZzFuPHFw\nQebHIrIQ1+vtjkSuGMCQHm3JTN/5TF5mujBkt3Y1zGGMMS2Xpavxhg0bprHcZ7Nm83YOvet9Sst/\nGJQxJyOND685nFB+TjKaaIwxjY6ITFfVYXVNZ+lq4hQqyOHUod1I8wc3menCKcO6W6AxxpgILNjU\nw9iRfUkTF23SRBg7sk+KW2SMMY2TBZt6CBXksF/P9gCMGtDJjmqMMaYGFmzq6Yz9Xa/sYwd1SnFL\njDGm8bJgU0/9OuUDsKPSOloYY0xNLNjUU/jU2ZrNpSluiTHGNF4WbOqpXetMMtOFNVss2BhjTE0s\n2NSTiFCUl82aLdtT3RRjjGm0LNgkQFFBDsV2ZGOMMTWyYJMAofxsu2ZjjDG1sGCTAKF8O41mjDG1\nsWCTAEX52WzYtoOyQJ40Y4wxP7BgkwDh7s/FJXYqzRhjIrFgkwCh/GzAZYI2xhizKws2CRAq8MHG\neqQZY0xEFmwSoCqLgAUbY4yJyIJNAhTmZSECxXYazRhjIrJgkwAZ6Wl0yM2yDgLGGFMDCzYJUpSf\nYzd2GmNMDSzYJIi7sdOCjTHGRGLBJkEsi4AxxtTMgk2ChAqyWVtSRoUNomaMMbuwYJMgofwcKiqV\n9VvLUt0UY4xpdCzYJEhVFgE7lWaMMbuwYJMglkXAGGNqZsEmQaqScVr3Z2OM2YUFmwQpstNoxhhT\nIws2CZKTmU5+ToadRjPGmAgs2CSQDQ9tjDGRWbBJoFB+jp1GM8aYCCzYJFCoINuScRpjTARJCzYi\n0l1E3heRBSIyT0Qu8+V/EpEvRORzEXlFRNoG5rlORBaJyJcicnSgfJQvWyQi1wbKe4nIZyKyUERe\nEJEsX57tXy/y7/dM1noGhU+jqVoWAWOMCUrmkU05cKWq7gUMBy4Wkf7ARGCgqu4NfAVcB+DfOx0Y\nAIwCHhSRdBFJB/4GHAP0B87w0wLcCdyjqn2BDcD5vvx8YIOq9gHu8dMlXSg/h9LySjZvL2+IxRlj\nTJORtGCjqitVdYZ/vgVYAHRV1bdVNfxtPBno5p+fADyvqqWq+g2wCNjfPxap6teqWgY8D5wgIgIc\nAbzk538SODFQ15P++UvASD99UoVv7Cy26zbGGLOTBrlm409j7Qt8Vu2tnwNv+OddgaWB95b5sprK\nOwAbA4ErXL5TXf79TX766u26QESmici04uLieFZtJ1X32liPNGOM2UnSg42I5AH/Ai5X1c2B8utx\np9qeDRdFmF3jKK+trp0LVB9R1WGqOqyoqKjmlYhSOIuA3WtjjDE7y0hm5SKSiQs0z6rqy4HyMcDx\nwEj94Wr6MqB7YPZuwAr/PFL5WqCtiGT4o5fg9OG6lolIBtAGWJ/IdYvkh/xodhrNGGOCktkbTYDH\ngAWqenegfBRwDTBaVbcFZpkAnO57kvUC+gJTgKlAX9/zLAvXiWCCD1LvA6f4+ccArwXqGuOfnwK8\npw3QRSw/O4OczDQ7jWaMMdUk88jmYOBsYI6IzPJlvwXuB7KBif6a/WRVvVBV54nIeGA+7vTaxapa\nASAilwBvAenA46o6z9d3DfC8iNwGzMQFN/zfp0VkEe6I5vQkrmcVEfE3dlqwMcaYoKQFG1X9H5Gv\nnbxeyzy3A7dHKH890nyq+jWut1r18u3AqbG0N1FseGhjjNmVZRBIsKL8bDuyMcaYaizYJFgoP9vG\ntDHGmGos2CRYqCCHLaXlfF9WkeqmGGNMo2HBJsHCN3YW26k0Y4ypYsEmwUI2YqcxxuzCgk2CWRYB\nY4zZlQWbBKvKIrDZjmyMMSbMgk2CtW+dRUaa2JGNMcYEWLBJsLQ0oTDP7rUxxpggCzZJECqwYGOM\nMUEWbJLADQ9t12yMMSas1mDjh2V+pqEa01wU5efYfTbGGBNQa7DxWZeLfGp/E6VQfjbrtpaxo6Iy\n1U0xxphGIZqsz0uAj0VkArA1XBgco8bsLNz9eW1JKZ3btEpxa4wxJvWiCTYr/CMNyE9uc5qHqhs7\nN1uwMcYYiCLYqOrNACKS715qSdJb1cQVVaWsses2xhgDUfRGE5GBIjITmAvME5HpIjIg+U1ruiw/\nmjHG7Cyars+PAFeo6m6quhtwJfBocpvVtBXmWeZnY4wJiibY5Krq++EXqjoJyE1ai5qBrIw02udm\n2Wk0Y4zxoukg8LWI/B542r/+GfBN8prUPLgbOy3YGGMMRHdk83OgCHjZPwqB85LZqOagKD+bYrtm\nY4wxQB1HNiKSDvxWVcc2UHuajVB+DovWWMc9Y4yB6DIIDG2gtjQroYJsireUUlmpqW6KMcakXDTX\nbGb67AEvsnMGgZeT1qpmIJSfTXmlsmFbGR187zRjjGmpogk27YF1wBGBMsVdvzE1CA4PbcHGGNPS\nRXPN5nNVvaeB2tNsVA0PvaWUvTqnuDHGGJNi0VyzGd1AbWlWqrII2Lg2xhgT1Wm0T0TkAeAFdr5m\nMyNprWoGgqfRjDGmpYsm2Bzk/94SKFN2voZjqmmVlU5+doalrDHGGKLL+nx4QzSkOSoqyLZknMYY\nQ3RZnzuKyGMi8oZ/3V9Ezk9+05q+orxsO7IxxhiiS1fzBPAW0MW//gq4vK6ZRKS7iLwvIgtEZJ6I\nXObL24vIRBFZ6P+28+UiIveLyCIR+VxEhgTqGuOnXygiYwLlQ0Vkjp/nfhGR2pbR0EIFOXbNxhhj\niC7YFKrqeKASQFXLgYoo5isHrlTVvYDhwMUi0h+4FnhXVfsC7/rXAMcAff3jAuAhcIEDuBE4ANgf\nuDEQPB7y04bnG+XLa1pGgwon41S1LALGmJYtmmCzVUQ64DoFICLDgU11zaSqK8M91lR1C7AA6Aqc\nADzpJ3sSONE/PwF4Sp3JQFsR6QwcDUxU1fWqugGYCIzy7xWo6qfqvs2fqlZXpGU0qFB+Nt/vqKCk\ntDwVizfGmEYjmt5oVwATgN4i8jEuA/QpsSxERHoC+wKfAR1VdSW4gCQiIT9ZV2BpYLZlvqy28mUR\nyqllGdXbdQHuyIgePXrEskpRCd7YmZ+TmfD6jTGmqYimN9oMETkM2AMQ4EtV3RHtAkQkD/gXcLmq\nbvaXVSJOGmnxcZRHTVUfwY1EyrBhwxJ+rqvqXpvNpfQuykt09cYY02REcxoNVS1X1XmqOjfGQJOJ\nCzTPBhJ3rvanwPB/1/jyZUD3wOzdgBV1lHeLUF7bMhpUVRYB6/5sjGnhogo28fA9wx4DFqjq3YG3\nJgDhHmVjgNcC5ef4XmnDgU3+VNhbwFEi0s53DDgKeMu/t0VEhvtlnVOtrkjLaFDhIxvr/myMaemi\nuWYTr4OBs4E5IjLLl/0WuAMY7+/V+Q441b/3OnAssAjYhh8NVFXXi8itwFQ/3S2qut4/vwjXNbsV\n8IZ/UMsyGlRBqwyyMtKs+7MxpsWrMdgE73OJpK7caKr6PyJfVwEYGWF6BS6uoa7HgccjlE8DBkYo\nXxdpGQ1NRHz3ZzuNZoxp2Wo7svmL/5sDDANm44LH3rheZYckt2nNQyg/245sjDEtXo3XbFT1cJ8X\n7VtgiKoOU9WhuC7MixqqgU1dKN+yCBhjTDQdBPZU1TnhF6o6FxicvCY1L6ECO41mjDHRdBBYICL/\nAJ7B3cfyM1w2ABOFUH42m7eXs31HBTmZ6alujjHGpEQ0RzbnAfOAy3AJOOf7MhOFIn+vjXV/Nsa0\nZNFkENguIg8Dr6vqlw3QpmYlOGJn9/atU9waY4xJjWjGsxkNzALe9K8Hi8iEZDesufjhyMau2xhj\nWq5oTqPdiEvtvxFAVWcBPZPYpmYlmIzTGGNaqmiCTbmq1jmkgImsQ242aeKScRpjTEsVTW+0uSJy\nJpAuIn2BscAnyW1W85GeJhTmZVsyTmNMixbNkc2lwACgFHgON3BancNCmx+ECiyLgDGmZav1yEZE\n0oGbVfU3wPUN06TmJ5Sfw6pNdmRjjGm5aj2yUdUKYGgDtaXZsvxoxpiWLpprNjN9V+cXga3hwsBg\naKYOofxs1m0tpbyikoz0pA0hZIwxjVY0waY9sA44IlCmgAWbKBUV5KAK67aW0bEgJ9XNMcaYBhdN\nBgFLTVNPVcNDby61YGOMaZHqDDYikgOcj+uRVvVNqao/T2K7mpWqYLNlO9AmtY0xxpgUiOYCwtNA\nJ+Bo4AOgG7AlmY1qbkL+aMaScRpjWqpogk0fVf09sFVVnwSOAwYlt1nNS2FeFmApa4wxLVc0wWaH\n/7tRRAbizgP1TFqLmqHsjHTats60LALGmBYrmt5oj4hIO+D3wAQgD7ghqa1qhkL52ZYfzRjTYkXT\nG+0f/ukHwO7JbU7zFcrPsdNoxpgWK5reaBGPYlT1lsQ3p/kK5WfzzdqtdU9ojDHNUDSn0YLfkDnA\n8cCC5DSn+SoqyKZ4SymqioikujnGGNOgojmN9pfgaxH5M+7ajYlBKD+HsopKNm7bQbvcrFQ3xxhj\nGlQ8ibpaY9duYvbDjZ123cYY0/JEc81mDi4XGkA6UATY9ZoYBbMI7NEpP8WtMcaYhhXNNZvjA8/L\ngdWqWp6k9jRb4SwC1v3ZGNMSRRNsqqemKQhe4FbV9QltUTNlp9GMMS1ZNMFmBtAd2AAI0Bb4zr+n\n2PWbqORmZ5CblW5ZBIwxLVI0HQTeBH6iqoWq2gF3Wu1lVe2lqhZoYhAqyLFknMaYFimaYLOfqr4e\nfqGqbwCH1TWTiDwuImtEZG6gbLCITBaRWSIyTUT29+UiIveLyCIR+VxEhgTmGSMiC/1jTKB8qIjM\n8fPcL/7cnoi0F5GJfvqJPtVOo1Bkw0MbY1qoaILNWhH5nYj0FJHdROR63MiddXkCGFWt7C7gZlUd\njMuvdpcvPwbo6x8XAA+BCxzAjcABwP7AjYHg8ZCfNjxfeFnXAu+qal/gXf+6UQjlZ9uRjTGmRYom\n2JyB6+78CvCqf35GXTOp6odA9c4DChT4522AFf75CcBT6kwG2opIZ9wYOhNVdb2qbgAmAqP8ewWq\n+qmqKvAUcGKgrif98ycD5SkXys9hzWa7ZmOMaXmiySCwHrgMQETSgVxV3Rzn8i4H3vJZCNKAg3x5\nV2BpYLplvqy28mURygE6qupK3/aVIhKqqTEicgHu6IgePXrEuUrRK8rPZmtZBVtLy8nNjqZvhjHG\nNA91HtmIyHMiUiAiucA84EsR+U2cy7sI+LWqdgd+DTwWXkyEaTWO8pio6iOqOkxVhxUVFcU6e8ys\n+7MxpqWK5jRaf38kcyLwOtADODvO5Y0BXvbPX8RdhwF3ZNI9MF033Cm22sq7RSgHWO1Ps+H/romz\nrQkXKvDBxk6lGWNamGiCTaaIZOKCzWuquoM4jiK8FfzQk+0IYKF/PgE4x/dKGw5s8qfC3gKOEpF2\nvmPAUcBb/r0tIjLc90I7B3gtUFe419qYQHnKhfJ9FgE7sjHGtDDRXDj4O7AEmA18KCK7AXVesxGR\nccAIoFBEluF6lf0/4D4RyQC246+X4I6YjgUWAduA88BdLxKRW4GpfrpbAhkLLsL1eGsFvOEfAHcA\n40XkfNzNp6dGsY4Nwk6jGWNaqmg6CNwP3B9+LSLfAYdHMV9NPdaGRphWgYtrqOdx4PEI5dOAgRHK\n1wEj62pfKrRtnUlWepplETDGtDgxd4nygcESccZBRCjKz6bYknEaY1qYeMazMfVgWQSMMS2RBZsG\nFsrPttNoxpgWJ6rTaCJyENAzOL2qPpWkNjVroYJspi6xURmMMS1LNCN1Pg30BmYBFb44nCLGxCiU\nn8OGbTsoK68kK8MOLI0xLUM0RzbDcDd2xntvjQkId38uLimla9tWKW6NMcY0jGh+Ws8FOiW7IS2F\nZREwxrRE0RzZFALzRWQKUNWNSlVHJ61VzZhlETDGtETRBJubkt2IlsSyCBhjWqJoMgh80BANaSna\n52YhAsV2Gs0Y04JEM8TAcBGZKiIlIlImIhUiEu94Ni1eRnoaHXLtxk5jTMsSTQeBB3Ajcy7EJb38\nhS8zcQpZFgFjTAsT1U2dqrpIRNJVtQL4p4h8kuR2NWuhAssiYIxpWaIJNttEJAuYJSJ3ASuB3OQ2\nq3kL5Wczf4WdiTSmIR1730fMX7nr565/5wJev+zQFLSoZYnmNNrZfrpLgK24kTP/L5mNau5C+Tms\nLSmlotLukzWmoQzp0ZbM9J1HlM9MF4bs1i5FLWpZ6gw2qvotIEBnVb1ZVa9Q1UXJb1rzFSrIplJh\n3Va7bmNMQznvkF67lKWLMHZknxS0puWJJjfaT4A/A1lALxEZjBsx027qjFNVypotpVU3eRpjahbL\nKTBVZfXmUuav3MS85ZuZv3Iz81Zs5rv123aaLiNdOGVYd/sMNpBob+rcH5gEoKqzRKRn0lrUAhQF\nsggMSHFbjGkKhvRoy8I1W9hR8cOp58x0Yd8ebVm0poR5KzYxf+Vm5q9wj3Vby6qm61WYy6CubTht\nv+50aZvDNS/NoayikooK5byDeqZgbVqmaIJNuapuEpG6pzRRqTqysRE7jYnK2JF9eXH6MlzCeae8\nUnlp+lKe/ew7ALLS0+jXKY8j9+pI/y4FDOhSwJ6dC8jL3vlrbvqSDVXz3DBhLk+ctz+Z6ZaBPdmi\nCTZzReRMIF1E+gJjAev6XA9FVSlrrPuzMdFYt7WMTgU5fBs4FRbKy+a4vbswoEsB/bsU0CeUF1XQ\nGDuyL1+tKWHUwE7c8u/53DRhHredOBD7QZ1c0QSbS4HrcUk4xwFvAbcms1HNXU5mOm1aZdqNncbU\nYfbSjTzw/iImzl9N66x00tOEikolJyONf489JK7rLaGCHMb/8kAA1mwu5eEPFtMnlMd5B+/agcAk\nTjS50bbhgs31yW9OyxHKz2aNnUYzJqLPvl7HA+8v4qOFa2nTKpNfH9mPcw/qyZ/e+oJnp3yXsAv7\nVx+9B4uLS7j1P/PpVZjLiD1CCWi9iaTGYCMiE2qb0Xqj1Y9lETBmZ6rKhwvX8rf3FjFlyXoK87K4\n7pg9OWv4blXXXcKnwBLVXTn5NbzoAAAgAElEQVQtTbj3tMGc8vCnXPrcTF7+1UH07ZifkLrNzmo7\nsjkQWIo7dfYZ7l4bkyBFedlM+3ZDqpthTMpVVirvLFjNA+8v4vNlm+jcJoebRw/gtP26k5OZvtO0\nwVNgiZKbncFjY4Yx+oGP+fmTU3n1VwfTIS87ocswtQebTsCPcUk4zwT+C4xT1XkN0bDmLlSQw5ot\npaiqXZg0LVJFpfKfz1fw4PuL+XL1Fnbr0Jo7Th7EyUO6kZXRsL3DurRtxaPnDOW0RyZz4TPTeeYX\nB5CdkV73jCZqNQYbn3TzTeBNEcnGBZ1JInKLqv61oRrYXIXysykrr2Tz9+W0aZ2Z6uYYkzQ13ZCZ\nlS6UVSh9Q3nce9pgjt+7Mxkp7IK8b492/PnUfRg7bibXvzKXP52yt/0QTKBaOwj4IHMcLtD0BO4H\nXk5+s5q/YPdnCzamOYt0QyZAXk4GfzhpEEf170RaWuP4Uh+9TxcWrynhvncX0ieUx4WH9U51k5qN\n2joIPAkMBN4AblbVuQ3WqhYgFMgiYBckTXNVWakc0reQcVO+26k8M11487IfESpofKliLhvZl0XF\nJdz55hfsXpjLUQM6pbpJzUJtRzZn47I89wPGBg4nBVBVLUhy25q1UIHd2Gmar8XFJbwyYzmvzFzO\n8o3fk57msv5Wqgs0p+3Xo1EGGnA91P5y6j4sW7+Ny1+YxYsXHsiALm1S3awmr7ZrNpa/IYnCKWvs\nXhvTXKwtKeU/s1fwyszlzF62iTSBQ/oWcdXR/RjSox1H3fMhpeWVTSLTck5mOo+e43qo/b8np/Hq\nJQdbws56imqkTpN4edkZtMpMp9iyCJgmbPuOCibOX80rM5fzwVfFVFQq/TsX8Lvj9mL0Pl12Ono5\ndWi3hN6QmWyhghz+MWYYpz78KRc8NZ3nLxi+S1dsE72kBRsReRw4HlijqgMD5ZfiBmIrB/6rqlf7\n8uuA84EKYKyqvuXLRwH3AenAP1T1Dl/eC3geaA/MAM5W1TLfqeEpYCiwDjhNVZckaz3jJSL+xk4L\nNqbxqqkn2W4dWnNAr/a8MWcVW0rL6VSQwy8O7cXJ+3Zjj06Rr0Em+obMhjCwaxvuOW0wFz4znatf\n+pz7Th9sPdTilMwjmyeAB3Bf/ACIyOHACcDeqloqIiFf3h84HRgAdAHeEZF+fra/4e73WQZMFZEJ\nqjofuBO4R1WfF5GHcYHqIf93g6r2EZHT/XSnJXE94xbKtywCpnGrqSfZt+u2sXZLKccM6sxJ+3Zl\n+O4dSK+jR1kybshsCKMGduLqUXtw15tf0rsoj8uO7JvqJjVJSQs2qvphhHFvLgLuUNVSP80aX34C\n8Lwv/0ZEFuHG0AFYpKpfA4jI88AJIrIAOAJ3synAk7hxdx7ydd3ky18CHhARUdVGNwZzKD+HBat2\n/dVoTCLEMuBYmKqyZkspi9aUsGhNCdvLKyivNnx5msDNowdwytDutMpqGaeVLjqsN4vWlHDPO19x\nzztf7fJ+bdvUOA19zaYfcKiI3A5sB65S1alAV2ByYLplvgxcypxg+QFAB2CjqpZHmL5reB5VLReR\nTX76tdUbIyIXABcA9OjRo94rF6ui/Gw+/MpOo5nkqGnAsSG7taO8opJv129j8ZoSFhWXsHjNVhYV\nl/D1mhK2lJZXTZ+XnUH71lls2Fa2U0+ysw/smYI1Sh0R4Y8nD+Kjr4opLinb6b3wNjW1a+hgkwG0\nA4YD+wHjRWR3IuddU1xvyUjlNU1PHe/tXKj6CPAIwLBhwxr8yCdUkM2W0nK+L6toMb8QTcOJNOBY\nRaXy8cJi9rrhzZ2CUMeCbHoX5XHSkK70CeXRuyiPPqE8QvnZFG8p5dC73m8yPcmSJTsjnafPP4BR\n9320U3lL3iaxaOhgswx42Z/SmiIilUChL+8emK4bsMI/j1S+FmgrIhn+6CY4fbiuZSKSAbQB1idp\nferlhxs7t7Nbh9wUt8Y0N6GCHI4e0JEJs1dWleVmZdCnYz5HD+xMn5ALKLsX5VKQU3MWi1BBTpPr\nSZYse3Yu4PhBnfnPHLdNM9OlxW+TaDV0sHkVd61lku8AkIULHBOA50TkblwHgb7AFNxRSl/f82w5\nrhPBmaqqIvI+cAquR9oY4DW/jAn+9af+/fca4/UaCNxrs6XUgo1JuI8WFvPeF2uqXmdnpPHuVYfF\n9cXYFHuSJcsNP+nPW/NXsaNCqVRsm0QpaTduisg43Bf+HiKyTETOBx4HdheRufggoc48YDwwH5f8\n82JVrfBHLZfgRgddAIwPZJ2+BrjCdyboADzmyx8DOvjyK4Brk7WO9VVkN3aaJFBVnvj4G87951S6\ntm3NiYO7IAKn1uMXeLgnmf2Cd9vitGHuhEtFpTJ3+aYUt6hpSGZvtDNqeOtnNUx/O3B7hPLXgdcj\nlH/NDz3WguXbgVNjamyKhPItZY1JrLLySm6cMI9xU77jyL06cu/pg9lWWs6KTdvtF3gCjR3Zly9W\nbWHT9zu4YvxsXh97KF3atkp1sxo1S0mTQu1aZ5GRJnZjp0mI9VvLOPuxzxg35TsuGtGbR84eSl52\nhh2VJEGoIIeXLjqIR84Zxo7ySsaOm8mOispUN6tRs2CTQmlpQlF+tp1GM/X21eotnPC3/zFz6Ubu\nOW0frhm1Z6NJ29+c9SrM5Q8nD2Latxu4e+Ku99+YH1hutBSzLAKmvt5dsJrLnp9Fq6x0XrhgOPv2\nsHs+GtIJg7sy+et1PDRpMQf0as+IPUKpblKjZEc2KVaUn2PJOE1cVJW/f7CYXzw1jZ6FrZlwycEW\naFLkxp8MYM9O+VwxfjarNtmPx0gs2KRYqCDbgo2J2fYdFVz54mz++MYXHDuwMy/+8iA6t7EL1KmS\nk5nOA2cOYfuOCsY+P5Nyu36zCws2KRbKz2bd1jK7uGiitmbLds58dDIvz1jOr4/sxwNn7msZKBqB\nPqE8bjtxIFO+Wc997y5MdXMaHbtmk2LhHkJrS0rtl6mp09zlm7jgqWms31bGg2cN4dhBnVPdJBNw\n8pBufLp4HQ+8v4gDenXgkL6FqW5So2HBJsWCI3ZasDFhNWVsFqBTmxxeuvAgBna1oYobo5tPGMCs\npRu5/IWZvD720EY7/HVDs9NoKRYq+CFljTFhQ3q0JTN9167L7fOyeO2Sgy3QNGKtszL421lDKCkt\n57LnZ1FR2SizZTU4CzYpFkzGaUzY2JF9Sas2ImSawKu/OthuzmwC+nXM55YTBvLp1+v463t2/QYs\n2KRcYV4WIpYfzewsVJDD6H26VI2XkS5w5v496N6+dUrbZaJ36tBunLxvV+57dyGfLN5lOK0WRxpp\nQuQGN2zYMJ02bVqDLjOekRRNy7B84/ec+cinfLv+ewByMtL48JrD7aimidlaWs7oB/7H5u3lvD72\n0Krku82JiExX1WF1TWdHNikU6by8jfpnvli1mZMf/Jj123bw471CiGBjpjRRudnu+s3m73dwxfhZ\nVLbg6zfWGy2FIo2kCHD83p2pqFTSGyi3VSKOsOwoLTE++3odv3hqGq2z0nnxwgNp3zqLTeNmWsbm\nJmzPTgXcNHoA1708hwcnLeKSI/rWu86m+HmzYJNC4REQx01dWtVjZUeFcvojk8nNSmdA1zYM7t6W\nvbu1YZ9ubenWrhVS7aJxIna62saqj1Yi6mjp3py7krHPz6J7u1Y8+fP96dbOXZ8Z/8sDU9wyU1+n\n79edTxev4+6JX7Ffz/YcsHuHetXXFD9vFmxSLHx0U1Gp5GSk8eTP92fphu/5fNlGZi/bxBMfL6HM\nZxdon5vFoK5t2KdbG/bp3pa9u7WNeqcrK69k3dZS1m4pY+3WUtZuKWVtSRnrSkpZW1K6S/fM8grl\nk0XFHH3Ph4iAiCBAWhoIQpoA4v6Kn768Wh02Nnv0np78LTe8Npd9u7flsTH70S43K9VNMgkkIvzh\n5EHMWb6Jsc+7+2865MV2/aa8opLFxVuZt2IT5ZVN7/NmwSbFqo/vfsDuHTgAOGVoN8AFiS9XbWHW\nso18vnQjny/bxEcLiwnvZx3zs3cJFJUKKzdu46d//5S1JS6wbN5eHnH5OZlpFOZl0651Fuu3lqG4\n4NG9fSv6dSxAcUPfun4k4eeK+uWoKqqgKKH8bFb7XnVp4u6mtusMtVNV7p74FX99bxEj9wzxwJlD\nLPVMM5WXncEDZ+7Lcff/j6G3vbPL+8GzEaXlFXy1qoR5KzYxd8Um5i7fzIKVmyktdz88czLTaJ+b\nxfqSsqqT8B3ysvhu3bZG+5mzYNMI1Da+e1ZGGoO6tWFQtzYwfDfA9XCZt2Jz1dHP+1+spqS0omqe\ndIFv1m2jMC+bvToV0KFPFoV52RTmZdMhL/zc/c3NdrvAms3bOfSu9yktryQ7I42XLjoo5p02WEel\nwvyVm1lbUkphjL/gWoryikp+9+pcnp+6lNOGdef2kwaSkW59dpqzAV3asF/P9kxdsn6n8ow0ITc7\nnatfms3c5Zv5avWWqiOX/OwM+ncp4GfDd2Ng1wIGdmlDr8Jc1m8tq/q8ZaQJJaXlnPLwpwzbrR0X\nHtabI/YMNaoxjazrs5eKrs+JUj1QfBRnF9nfvTKHZ6d8x1kH7MZtJw6Mqy3hOg7q3YFpSzZQmJfN\nI+cMZUAXu+M96PuyCi4dN4N3Fqzh0iP6cMWP++1yPc40T6s3fc+Bd7xHpI5p7XOzGNi1DQO6uKAy\nsGsB3du1rjFoBD+zvz12T8ZPXcqjH33D8o3f0zeUxy8P683ofbqQlZG8HzHRdn22YOM15WADiQkU\nazZv55JxM3ngzH3jPhQP1rF6UykXPD2NDdvK+POp+3D83l3iqrO52bC1jPOfnMrMpRu5ZfQAzj6w\nZ6qbZBrY1S/N5sVpy1DcKecf9SvijycPolNBTkw/OiJ9ZndUVPLfz1fy8AeL+WLVFjq3yeH8Q3px\nxv49qs5kJJIFmxg19WCTiECRDMVbSrnomelM+3YDFx/emyt/vEejOrRvaMs3fs+Yx6fw3fpt3H/6\nYEYNtKzNLVHwbESybthVVSZ9VczDkxbz2TfradMqk7OH78a5B/dM6KltCzYxaurBpjErK6/kxglz\nGTdlKSP3DHHv6YPJz8lMdbMa3BerNjPm8SlsK6vg0XOGMbye3V9N05aIsxHRmvndBh7+YDFvz19N\nVnoapw7rxuTF61lUXLLLtLHeqxNtsLEOAibpsjLS+MNJg+jfuYCb/z2fkx78hEfPGUavwtxUN63B\nVL9Zc89OBalukkmx2joGJdq+Pdrx97OHsbi4hEc++JoXpi5lR4Ui7HxLeTLv1bEjG8+ObBrGp4vX\ncfFzMyivqOSvZw7hsH5FqW5SwtV0o21WuvDeVSOqbtY0JlVWb97OX99dyDOffbdTeTyn9Cw3mmmU\nDuzdgdcuPpgubVtx3j+n8OiHX9PcfvBEynknwImDu1qgMY1Cx4IcbjtpEKcO60Z4V81Ml6Tm4LNg\nYxpc9/ateflXBzFqYCduf30BV4yfzfYdFXXP2ESMHdl3lx5FWRlpXDVqjxS1yJjIfnPUHlX3diU7\nA4EFG5MSrbMy+NuZQ7jqqH68MnM5P/37p6za1PQHkNuyfQePffwN5f5Ob3C/GE+1rM2mEQpnMGmI\nzOLWQcCkjIhwyRF92aNTAZc/P5OD73wv4hC6jTmTbVhFpTJ+2lL+8vaXrC0p49hBnXl3wWpKyysb\nfc4q07I1VEcFO7IxKffj/h155eKDaZW56+7Y2DPZAnyyaC3H3f8R1708h16FuUy45GAePGtIg/1i\nNKY+QgU5jP/lgUnfR+3IxjQK/Trm88qvDuboez/cKY1HYz4q+GbtVv7w+gImzl9Nt3atePCsIRwz\nsFPV9ZqG7NpqTGOXtCMbEXlcRNaIyNwI710lIioihf61iMj9IrJIRD4XkSGBaceIyEL/GBMoHyoi\nc/w894v/hItIexGZ6KefKCKN+2exqdK3Yz6n79edYEeutq0z+bp4a6Pqsbbp+x3c/t/5HHXPB3yy\naC1Xj9qDd644jGMHdd6pY0BD/WI0pilI5mm0J4BR1QtFpDvwYyDYwfsYoK9/XAA85KdtD9wIHADs\nD9wYCB4P+WnD84WXdS3wrqr2Bd71r00TcfmR/ap6x2SkCWXlbjC5Ux/+lPe/XJPSoFNeUcnTk7/l\n8D9P4h//+4aT9+3G+78Zwa9G9CEn04YFMKY2STuNpqofikjPCG/dA1wNvBYoOwF4St03yWQRaSsi\nnYERwERVXQ8gIhOBUSIyCShQ1U99+VPAicAbvq4Rvt4ngUnANQlcNZNEwfF9Tt+/B787bi/GT1vK\nw5MWc94/pzKwawGXHN6Ho/p3SlqOtZpuyszOSKO0vJIDerXn98f3Z2BXy2RtTLQa9JqNiIwGlqvq\n7Gr3IXQFlgZeL/NltZUvi1AO0FFVVwKo6koRCSV0JUzSBa915GSmc86BPTl9vx68OnM5D05axIXP\nzKBvKI+LD+/D8Xt3TvgYMJFGPwXXWeG+04dy9ICONhyAMTFqsN5oItIauB64IdLbEco0jvJY23SB\niEwTkWnFxcWxzm6SJNK1jqyMNH66X3feueIw7jt9MCJw+QuzGHn3Bzw/5TvKAve11EdZeSU/2WfX\nTMwZacKbl/+IUYEOAMaY6DXkkU1voBcQPqrpBswQkf1xRybdA9N2A1b48hHVyif58m4RpgdYLSKd\n/VFNZ2BNTQ1S1UeAR8DlRot3xUzDyUhP44TBXfnJ3l2YuGA1D7y3iGtfnsP97y7kl4f1ZtyU7/hi\n1ZZd5qt+r87GbWUsLi5h8ZqtLF7r/n5dXMK367ftcq9PZrpw2n49LNWMMfXQYMFGVecAVae0RGQJ\nMExV14rIBOASEXke1xlgkw8WbwF/CHQKOAq4TlXXi8gWERkOfAacA/zVTzMBGAPc4f8Grw2ZZiIt\nTTh6QCeO6t+RDxeu5YH3FnLjhHlkZ6SRLhA8A5aeJrTKSuO6lz93waW4hHVby6rez0pPo1dhLnt0\nyufYQZ3pHcqlXessfvn0dLsp05gESVqwEZFxuKOSQhFZBtyoqo/VMPnrwLHAImAbcB6ADyq3AlP9\ndLeEOwsAF+F6vLXCdQx4w5ffAYwXkfNxPd5OTeBqmUZGRDisXxGH9Svis6/XcffEr/jsm53Hd6+o\nVKZ/u5Ela7fRuyiPowZ0ZPfCPHqHculdlEe3dq1Jj9DZINxRwW7KNKb+bIgBz4YYaD4ufGY6b89b\nRaVCehocuVdH7jh5b9rlZsVUT2Md/dSYxsSGGDAt1i2jB5Dpe6hlpqVx64kDYw40YDdlGpNIFmxM\ns9OQmWyNMdGx3GimWbK8ZMY0LhZsTLMUPgVmjGkc7DSaMcaYpLNgY4wxJuks2BhjjEk6CzbGGGOS\nzoKNMcaYpLMMAp6IFAPfprAJhcBaq6PZ1dEY2mB1WB3JrGM3VS2qayILNo2EiEyLJuWD1dG06mgM\nbbA6rI6GqKMudhrNGGNM0lmwMcYYk3QWbBqPR6yOZllHY2iD1WF1NEQdtbJrNsYYY5LOjmyMMcYk\nnQUbY4wxSWfBJsVE5HERWSMic+tRR3cReV9EFojIPBG5LI46ckRkiojM9nXcHGdb0kVkpoj8J875\nl4jIHBGZJSJxDZ0qIm1F5CUR+cJvk5jSP4vIHn754cdmEbk8jnb82m/LuSIyTkRiHlhHRC7z88+L\ntg2R9ikRaS8iE0Vkof/bLo46TvXtqBSROrvJ1lDHn/z/5XMReUVE2sZRx61+/lki8raIdIm1jsB7\nV4mIikhhHO24SUSWB/aTY2Ntg4hcKiJf+u16VxxteCGw/CUiMiuOOgaLyOTwZ05E9q+tjripqj1S\n+AB+BAwB5tajjs7AEP88H/gK6B9jHQLk+eeZwGfA8DjacgXwHPCfONdlCVBYz236JPAL/zwLaFuP\nutKBVbgb12KZryvwDdDKvx4PnBtjHQOBuUBr3HAg7wB949mngLuAa/3za4E746hjL2APYBIwLM52\nHAVk+Od3xtmOgsDzscDDsdbhy7sDb+Fu5q51n6uhHTcBV0X5v4w0/+H+f5rtX4fiWY/A+38Bboij\nHW8Dx/jnxwKTYtlPo33YkU2KqeqHwPp61rFSVWf451uABbgvu1jqUFUt8S8z/SOm3iMi0g04DvhH\nLPMlkogU4D5QjwGoapmqbqxHlSOBxaoaT3aJDKCViGTgAsaKGOffC5isqttUtRz4ADiprplq2KdO\nwAVh/N8TY61DVReo6pdRtr2mOt726wIwGegWRx2bAy9zqWM/reUzdg9wdV3z11FHVGqY/yLgDlUt\n9dOsibcNIiLAT4FxcdShQIF/3obY99OoWLBpZkSkJ7Av7sgk1nnT/WH4GmCiqsZax724D29lrMsO\nUOBtEZkuIhfEMf/uQDHwT3867x8ikluP9pxOHR/gSFR1OfBn4DtgJbBJVd+OsZq5wI9EpIOItMb9\n6uwea1u8jqq60rdtJRCKs55E+jnwRjwzisjtIrIUOAu4IY75RwPLVXV2PMsPuMSf0nu8rlOTEfQD\nDhWRz0TkAxHZrx7tOBRYraoL45j3cuBPfnv+GbiuHu2okQWbZkRE8oB/AZdX+/UXFVWtUNXBuF+b\n+4vIwBiWfTywRlWnx7rcag5W1SHAMcDFIvKjGOfPwJ0meEhV9wW24k4bxUxEsoDRwItxzNsOdzTR\nC+gC5IrIz2KpQ1UX4E41TQTeBGYD5bXO1ESIyPW4dXk2nvlV9XpV7e7nvyTGZbcGrieOIFXNQ0Bv\nYDDuB8VfYpw/A2gHDAd+A4z3RyjxOIM4fhR5FwG/9tvz1/izAolmwaaZEJFMXKB5VlVfrk9d/rTT\nJGBUDLMdDIwWkSXA88ARIvJMHMte4f+uAV4BYr1YuQxYFjgqewkXfOJxDDBDVVfHMe+RwDeqWqyq\nO4CXgYNirURVH1PVIar6I9zpj3h+uQKsFpHOAP5vradskklExgDHA2epv1BQD88B/xfjPL1xPwJm\n+/21GzBDRDrFUomqrvY/0CqBR4lvX33Zn8KegjsjUGtHhUj8adqTgRdindcbg9s/wf2wSkoHAQs2\nzYD/NfQYsEBV746zjqJwzyARaYX7svwi2vlV9TpV7aaqPXGnnt5T1Zh+yYtIrojkh5/jLibH1EtP\nVVcBS0VkD180EpgfSx0B9fm1+B0wXERa+//PSNy1tJiISMj/7YH7Qom3PRNwXyr4v6/FWU+9iMgo\n4BpgtKpui7OOvoGXo4lhPwVQ1TmqGlLVnn5/XYbrYLMqxnZ0Drw8iRj3VeBV4AhfVz9cZ5Z4Mi8f\nCXyhqsvimBfcNZrD/PMjiP8HTe2S0evAHtE/cF8eK4EduJ3+/DjqOAR3reNzYJZ/HBtjHXsDM30d\nc6mjV0sddY0gjt5ouOsts/1jHnB9nMsfDEzz6/Iq0C6OOloD64A29dgON+O+COcCT+N7HcVYx0e4\nYDkbGBnvPgV0AN7FfZG8C7SPo46T/PNSYDXwVhx1LAKWBvbTunqSRarjX36bfg78G+gaax3V3l9C\n3b3RIrXjaWCOb8cEoHOM82cBz/h1mQEcEc96AE8AF9Zj3zgEmO73sc+AofHu87U9LF2NMcaYpLPT\naMYYY5LOgo0xxpiks2BjjDEm6SzYGGOMSToLNsYYY5LOgo1p0nzG3r8EXl8lIjclqO4nROSURNRV\nx3JOFZed+v1kLyvVROS3qW6DSQ0LNqapKwVOritFfEMTkfQYJj8f+JWqHp6s9jQiFmxaKAs2pqkr\nx42f/uvqb1Q/MhGREv93hE98OF5EvhKRO0TkLHHj+cwRkd6Bao4UkY/8dMf7+dPFjcsy1Sdh/GWg\n3vdF5DnczX7V23OGr3+uiNzpy27A3VT3sIj8KcI8V/t5ZovIHb4sPP5IeEyYdr58kojcIyIf+iOl\n/UTkZXFj2Nzmp+kpbjyZJ/38L/lcYYjISJ+8dI5PLJnty5eIyM0iMsO/t6cvz/XTTfXzneDLz/XL\nfdMv+y5ffgcuC/YsEXnWz/9fv25zReS0GP7vpqlJxp2i9rBHQz2AElx69CW49OhXATf5954ATglO\n6/+OADbixgHKBpYDN/v3LgPuDcz/Ju5HWV/cHdc5wAXA7/w02bhsBb18vVuBXhHa2QWXwqYIl4Dx\nPeBE/94kIowPg8vN9gnQ2r9u7/9+Dhzmn98SaO8k/Pgwfj1WBNZxGS6LQE9ctomD/XSP+22Wg7uz\nv58vfwqX0BW/bS/1z38F/MM//wPwM/+8LW4cpVzgXOBr///IwY0X0z34P/DP/w94NPA67mwN9mj8\nDzuyMU2eugzXT+EG0orWVHXjAJUCi3EDSIE7IukZmG68qlaqS93+NbAnLmfbOeKGY/gM9yUeztc1\nRVW/ibC8/XCDUhWrG8/lWdy4O7U5Evin+hxiqrpeRNrgBoP7wE/zZLV6JgTWY15gHb/mh+EJlqrq\nx/75M7gjqz1wiUO/qqHecKLG6fywfY4CrvXbYRIusPTw772rqptUdTsu3c5uEdZvDu7I8U4ROVRV\nN9WxPUwTlpHqBhiTIPfi8kv9M1BWjj9V7JNhZgXeKw08rwy8rmTnz0X1fE6KG9X0UlV9K/iGiIzA\nHdlEEk/qeImw/LoE16P6OobXq6Z1iqbeikA9AvyfVhtQTUQOqLbs4Dw/LFT1KxEZihun548i8raq\n3lJHO0wTZUc2pllQ1fW4oZfPDxQvAYb65yfgRh+N1akikuav4+wOfIkbSvgiccM6ICL9pO4B2j4D\nDhORQt954AzcyJu1eRv4eeCaSnv/63+DiBzqpzk7inqq6yEiB/rnZwD/wyUM7SkifWKo9y3gUh/I\nEZF9o1j2jsB26wJsU9VncIN2xTsUhGkC7MjGNCd/YeeBtB4FXhORKbhMxzUdddTmS9yXbkdcZt3t\nIvIP3KmkGf6Ltpi6h1leKSLXAe/jjgheV9Va0/yr6psiMhiYJiJlwOu43lxjcB0KWuNOj50X4zot\nAMaIyN9xWaAf8ut1HvdJ5xAAAAB1SURBVPCiuPFRpgIP11HPrbgjys/9dliCG6emNo/46WfgTn3+\nSUQqcVmIL4pxPUwTYlmfjWlBxA0b/h9VjXoUVmMSwU6jGWOMSTo7sjHGGJN0dmRjjDEm6SzYGGOM\nSToLNsYYY5LOgo0xxpiks2BjjDEm6f4/i6TfwD8v//YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a18fa9358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pcr = PCR()\n",
    "kf_10 = KFold(n_splits=10, shuffle=True, random_state=1)\n",
    "\n",
    "pcr.fit(scale(X_train), y_train)\n",
    "\n",
    "X_reduced_train = pcr.get_transformed_data()  # Get the reduced dataset\n",
    "regr = pcr.get_regression_model()  # Get the undelying regression model\n",
    "n = len(X_reduced_train)\n",
    "\n",
    "mse = list()\n",
    "\n",
    "score = -1*cross_val_score(regr, np.ones((n, 1)), y_train.ravel(),\n",
    "                           cv=kf_10, scoring='neg_mean_squared_error').mean()\n",
    "mse.append(score)\n",
    "\n",
    "for i in np.arange(1, X.shape[1]):\n",
    "    score = -1*cross_val_score(regr, X_reduced_train[:, :i], y_train.ravel(\n",
    "    ), cv=kf_10, scoring='neg_mean_squared_error').mean()\n",
    "    mse.append(score)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(mse, '-v')\n",
    "ax.xaxis.set_ticks(np.arange(1, X.shape[1], 1))\n",
    "ax.set_title('Cross validation error from 0 to N components')\n",
    "ax.set_xlabel(\"Number of components\")\n",
    "ax.set_ylabel(\"Mean squared error\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Explained variance values: [ 39.92  61.42  73.4   80.8   86.33  90.14  93.38  95.78  96.82  97.76\n",
      "  98.4   98.81  99.18  99.5   99.75  99.89  99.95  99.98  99.98]\n"
     ]
    }
   ],
   "source": [
    "variances = np.cumsum(np.round(pcr.explained_variance_ratio_, decimals=4)*100)\n",
    "print(\"Explained variance values:\", variances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VNX9//HXOxskQEhAQCBAXHBX\nFOO+i120rVorWruI1tZ+rYq17be1td+q1fanVVvbb/vVWq3iUhV3a63autXWhU0QBFcI+xIIeyAh\nyef3xz2BIUySmUmGSTKf5+Mxj7n3zr3nnjuZ3M+955x7jswM55xzrrmcTGfAOedc5+QBwjnnXFwe\nIJxzzsXlAcI551xcHiCcc87F5QHCOedcXB4gughJ90q6IcF1/y5pXBryUC7JJOV1dNpx9nWcpA/S\nvZ+OJmm4pA2ScltZxyTtuTPz5VwqPEB0MEmVkjaFk0TT6/c7Mw9mdqqZTdiZ++xoZva6me2d6Xwk\ny8wWmFlvM2sAkPSqpG+2J01JV0paJmmtpD9L6tHCeu0K4JJODNv/odnyf0u6IJU0uzNJ10p6INP5\nSCcPEOnxhXCSaHpdlukMdSU74w6lq5D0GeAqYAxQDuwOXJfGXW4EzpdUnsZ9uC7CA8ROJOl2SY/F\nzN8k6SVFTpS0SNJPJK0MdyJfbSGdUknPSqqStDpMl8V8vvWqVdIF4QrwlrDuPEmnxqzbV9LdkpZK\nWizphqbiEUm5YbuVkuYCn2vl2K6KPbaw7LeSfhemL5Q0R9J6SXMlfTtmvaZj/5GkZcA9Tcuapf9J\n2H62pC/GfNbWMfaTdI+kJeHzp2I++7yk6ZLWSHpD0kEtHN91kv43TOdL2ijpV2G+UNLm8HfZehUv\n6RfAccDv49xJniLpo5CfP0hSC1/tOOBuM3vPzFYD1wMXtLDuv8L7mrC/oyTlSPqppPmSVki6T1Lf\nFrYHWAPcC1zTyjqx30tu+M02/W2mShoWPjta0uRw5zNZ0tEx270afmtvhLz+VVJ/SQ9KWhfWL49Z\n3ySND7+dlZJulpQTPmvxGGP+HuMkLQjbXh2Tbk7Mb2uVpImS+rW1raTPAj8Bzg35nxGWXxDyuD78\nDuP+D3cZZuavDnwBlcApLXxWBHxI9A9+HLASKAufnQjUA78GegAnEF3N7R0+vxe4IUz3B74U0usD\nPAo8FbOfV4FvhukLgC3At4Bc4BJgCaDw+VPAH4FewEBgEvDt8Nl/Ae8Dw4B+wCuAAXlxjm0EUAMU\nh/lcYClwZJj/HLAHoHBsNcDoZsd+Uzj2wrBsUUz6Y4EhRBc154bvZnCCx/g34BGgFMgHTgjLRwMr\ngCPCduPC369HnOM7GZgZpo8GPgHejvlsRpguj/2OYv8WMWkZ8CxQAgwHqoDPtvCbmQGcGzO/S9i+\nf5x1t9t3WPYN4GOiO4/ewBPA/S3s60RgEbArsI5tv71/Axe0sM1/AzOBvcPfdhTR77MfsBr4OpAH\nnBfm+8d8Lx+H30RfYDbR/8YpYf37gHuafWevhHSHh3W/2dYxxnwnfyL6XY0CaoF9w+ffBd4Cyoh+\ne38EHkpw22uBB2Ly2KvZ9zYY2D/T56R2nc8ynYHu9iI6wWwguhJren0r5vPDgWpgPnBezPITiU6S\nvWKWTQT+J0zfSwgQcfZ5MLA6Zv5Vtg8QH8d8VhR+9LsCg8IPvjDm8/OAV8L0y8B/xXz2aVoIEOHz\nfwPnh+lPAZ+08j09BVwRc+x1QM9m38eiVrafDpyRwDEOBhqB0jhp3A5c32zZB4QA0mx5IbCZ6OR3\nFdHV4yKiE9J1wO/Cek0nlbYCxLHN/s5XtXCcnxATPIgCnAHlcdbdbt9h2UvAd2Lm9yYKpvGC/Nbv\nHPgV8EjM3/WCFvL3QdPfodnyrwOTmi17symd8L1cHfPZrcDfY+a/AExv9p3Ffg/fAV5q6xhjvpOy\nmM8nAV8O03OAMTGfDU5i22vZMUCsIbp4K2z+nXTFlxcxpceZZlYS8/pT0wdmNgmYS3S1NbHZdqvN\nbGPM/Hyiq+btSCqS9MdwS72OqGihRC23nFkWs/+aMNmb6Ko/H1gailjWEF1BDQzrDAEWNstPa/5C\nFGAAvhLmm/J8qqS3JFWH/ZxGdDXcpMrMNreUsKTzY4qC1gAHNNu+pWMcBlRbVDzT3Ajg+01phnSH\nEec7N7NNwBSiu5/jgdeAN4BjwrLXWsp7C5bFTNeEvMazASiOmW+aXp/gfoaw/d9tPtHJb1Ab290E\nfEbSqDbWG0YUxNrab9O+h8bML4+Z3hRnvvl30vy32PR3SuQYW/q+RwBPxvz95wANCW67nfC/ey7R\nnfdSSX+TtE+8dbsKDxA7maRLiW5llwA/bPZxqaReMfPDw3rNfZ/oKukIMysmOmFBFHSSsZDoDmKX\nmGBWbGb7h8+XEp0AYvPTmkeBExXVh3yRECAUtbp5HLgFGGRmJcBzzfLbYrfCkkYQ3eZfRlREUQLM\nIrHjXQj0k1TSwme/aBbMi8zsoRbSeo2oOOkQYHKY/wzRXeG/Wtimvd0lv0dUtNFkFLDczFYluK8l\nRCfBJsOJ7lSXx1l3W0JR+rcR1Xm0ZiFRMVFb+23a9+I20mtN899i0/9GSscYLARObfYb6GlmieRz\nh+/bzF4ws08R3Ym8T/S77bI8QOxEkvYCbgC+RnQL/kNJBzdb7TpJBZKOAz5PdNJtrg/RFdaaUKF2\nTSr5MbOlwIvArZKKQ4XdHpJOCKtMBMZLKpNUSlS00lp6VURFB/cA88xsTviogCgoVgH1iiqQP51E\nVnsR/TNWQVThTXQHkegx/h34v1CJnC+pKaD+CfgvSUco0kvS5yT1aSG514DzgdlmVheO9ZvhWKta\n2GY5Udl4qu4DLpK0X/gb/JSouDGeKqLitNj9PQRcKWk3Sb2BXxIVHdUnsO9fE9W37NvKOncB10sa\nGb7DgyT1J7oA2EvSVxRV2J8L7EdU95Kq/w5/w2HAFUT1StC+Y7wD+EW4CEHSAElnJJif5UB5TGX5\nIEmnh4u8WqK7v4aEj64T8gCRHn/V9s9BPKmo6eYDwE1mNsPMPiIqx75f29q1LyOqyFsCPEhU/v9+\nnPRvIyoTX0lUwfZ8O/J6PtEJfHbY92NEVz8QnUBfIKoonUZU+deWvxBVNG4tXjKz9cB4ooCzmqj4\n6ZlEM2hms4nKqN8k+qc8EPhPotsTBeMtRFd0K4gqJjGzKUQV278P+fqYllsIQVSkVMi2u4XZRPUS\nLd09APwWOFtRa6XfJZFnQh6fJ6oPeIWo6GQ+LVwQhKK1XwD/CUUmRwJ/Bu4PeZwX8nt5gvteF/bd\nr5XVfk30d32RqIL2bqLy91VEFzjfB1YR3S1/3sxWJrLvFjwNTCWqf/pb2Be04xiJ/j7PAC9KWk/0\n/3REgts2XbytkjSN6Hz6faL/32qiosfvJJhWp9TUysNlmKQTiSq8ytpa17lsI8mAkWb2cabzkk38\nDsI551xcHiCcc87F5UVMzjnn4vI7COecc3F16U7RdtllFysvL890NpxzrkuZOnXqSjMb0NZ6XTpA\nlJeXM2XKlExnwznnuhRJbfWKAHgRk3POuRZ4gHDOOReXBwjnnHNxeYBwzjkXlwcI55xzcaWtFZOk\nPxN11rXCzA4Iy/oR9cBYTjSwzjlmtlqSiDrNOo2ov/ULzGxauvLmnHPJOO23rzN76bodlu83uJjn\nrjiuS6WRjHQ2c72XqJfM+2KWXUU0CtSNkq4K8z8CTgVGhtcRRCN9JdqjonOum+osJ9XRw0v4aMV6\ntjRs63kiP1eMHlGa0PadKY1kpC1AmNm/FDPoeHAG0bCGABOI+tP/UVh+n0X9frwlqUTS4NCXv3Ou\nC8rEibmh0airb6SuvpHahgbq6hvZc2AvPly+nvrGbWnk5YghJT15/aMqGg0aG42GRqPBLJq2aL7R\njIZGKO/fi+a9EpnBoD49uPvf88IQndBohhHetw7pDI0GhQW5NDZLo9GgR6749YsfJPR99MjL2SGN\nXInxY/ZMaPtk7ewH5QY1nfTNbKmkpqEth7L9cIKLwrIdAoSki4GLAYYPb2uAM+dcKtJ5ct93cB/m\nrdzIxtp6NtbWU1PXwIbaemrq6tlQ20BNbT0b6uqpqW1g1cY6GpqdEesbjH99uIJjbnyZuobGrQGh\nrqFxh3VbUt9o/HPOCv45Z0VC67eUxq3/+DDl7SEKaH9+ozKpbWIDVX6uOLtiGAP79GxXPlrSWZ6k\njjd0ZNy/tJndCdwJUFFR4T0NOtdMOk/u+w+JTu7rNm1h3eYtrNtUH953nF+1oY76hu3/Rbc0GI9P\nW8zj01of0bMgN4dePXIpKsijT8981m3aghGdKIaWFnJQWQkFeTn0yMuhIDeHgrzwys3dNp2XQ4/w\n2ePTFvHGJ6toaDTycsSJew/gm8ftTo5Ebg7hXVvft5uWyMmB1TV1nH37m9TWN9IjL4e/Xn4sA/v0\nQAiFNER4F9ELkSOQoveq9bUc96tXqK1vpGdeDv/60UlJn9xXrNu8NY103j3Azg8Qy5uKjiQNJhrd\nC6I7htjxZsuIPxazc93aziiW2bylgTU1W1hdU8eami2sqaljzaYt26ZrtrBs7abtimQgOrk/OnUx\nj06Nf3LPzRF9euZR3DOf4sLovay0kMVrNtFokCM4eFgJXxxdRu9w8u/dI4+igtzovUcevQqi5QV5\n2xpYxp4Qe+Tl8MR3jk76pHr0Hv057levbA0QvzzrwKTTKCstYuyhZTw4aQFjK4ax16CWRqZt2cDi\nnlvTSPXKvyPSSNTODhDPAOOAG8P70zHLL5P0MFHl9Fqvf3DZKNky9011DVTX1LF6Yx2rNkbvu/Tp\nsUM5dX2D8cKsZTw2dSGbtzS2uP+CvBxKi/IpKSxgQO8eVK2vxYhO7gcO7ctZo8u2nvyLC/O3CwZF\nBblEDRK3iT25F+TmcMfXD036hNaZTqrjx4zkwxUb2nXV3lnSSETaxoOQ9BBRhfQuROMIXwM8RTR+\n7XBgATDWzKpDM9ffA58lauZ6YRgvuFUVFRXmnfW57mT52k0cf/Or1NZvO4nn5YizDy2jtr5xaxCo\nDq9NWxraTLOpWOao3ftT2quAvoX5lBYVUFKUH70KCyjtFb0XFuRu3S725J5qcQjAT5+cyYOTFvDV\nI0Zww5kHJL19U14ue+gdfv+VQ1I+uXdEGt2FpKlmVtHmel15wCAPEK4zSaR4yMyo3ljHotWbwquG\nRas3sXjNtumauh1P+r0KcunXu4B+RQX061VAaa8C+of3pmWxr811DZxwy6vd5uTuOlaiAaKzVFI7\n1+XFKx7KFdQ3NnLBPZNYHIJC86v+4p55lJUWUd6/F8fuOYC+hXn8/pWP2dJg9MjL4Z/fO4Fh/YqS\ny0wRnaZIZWBxTyZ++6iUt3eZ4wHCuRSZGYtWb+L9ZeuZs3Qdi9ds2qHVToPB0rWbyc/NYfcBvTh+\nrwGUlRZSVlpEWWkhQ0sLKe6Zv0PaVetrt1aGJh0cAj+5u/byAOFc0FoR0eOXHM2Hy6NAEL3WM2fZ\nOtZvrt+6Xnn/6KS/ZM0mGiyqO/jS6KHcdPaopPPiJ3fXGXiAcC6IV0QkYOHqGva/5vmtLYN6FeSy\nz+BiTh81hH0HF7Pv4GL22bUPvXrkba3YbahvJC9HfP8ze6eUFz+5u87AA4TLaivWbebdRWt5d/Fa\n5q7csF1waDJ6eAmjhpWy3+A+7Du4mGGlReTkxHu2c+e2UXcu3TxAuG4hkRZEKzfUMnPRWt5dtJaZ\ni9cwc/Falq+rBaJ2/iMH9mH3XXoxf9VGGix6/uDcw4Yn3XpnZ7VRdy7dPEC4biFe8VBejsjPFd++\nfwozF61lydrNQNQFwh4DenPMHrtwwNC+HFTWl/2GFFNUsH0RUardGHjxkOsuPEC4bmH8mJFMnLJw\nu2X1jcaMRWtZt7meivJ+HFTWlwOH9mX/oX3p3SP+T9+LiJzbxgOE69LMjFc+WMEdr86lrtnzByft\nM5BbzzmYvoU7NiNtjRcRORfxAOG6pPqGRp59dyl3vPYJ7y9bz5C+PbnylJH836ufUFvfSH5uDr88\n68CkgwN4EZFzTTxAuC5lU10DE6cs5M5/zWXxmk2MHNibW8aO4vRRQyjIy9n6gJkXDznXfh4gXJew\npqaO+96cz71vVFK9sY5DR5Ry7en7M2afgds1OfXiIec6jgcI16ktXbuJu16fx0OTFlBT18DJ+wzk\nkhP34LDyfnHX9+Ih5zqOBwiXcS09w9C3MJ+aunoaDU4fNYRvn7A7++xanIEcOpedPEC4jIv3DAPA\nhs1b+PpR5Vx07G4pd1jnnEudBwiXcdEzDIuIHYY8N0c8N/5Y9vY7BucyxgOEy6jNWxp44K35bGnY\nNoJaUxcXHhycyywPEC5j3vxkFVc/OZO5Kzdy6gG78vL7K6htRxcXzrmOlZPpDLjss6amjh8+NoPz\n/vQW9Y3G/Rcdzu1fO5Sxh5Yh4c8wONdJ+B2E22nMjGdmLOHnf53Nmk1buOTEPRh/8kgKC3IBf4bB\nuc7GA4TbKRasquHqp2by+kcrGTWshAfOOpB9B29fx+DPMDjXuXiAcGm1paGRu/89j9v++SG5Eted\nvj9fO3IEuS0MuOOc6zw8QLi0mb5wDVc9/i7vL1vPp/cbxHVn7M/gvoWZzpZzLkEeIFyH21Bbzy0v\nfMCENysZ2KcHd3ztUD57wK6ZzpZzLkkZCRCSrgC+RTQm/J/M7DZJ/YBHgHKgEjjHzFZnIn8ucS11\nkwEw7qgR/OAze9OnZ/JdbjvnMm+nN3OVdABRcDgcGAV8XtJI4CrgJTMbCbwU5l0nN3p4Cfm5O9Yn\nfHb/QVx3xgEeHJzrwjLxHMS+wFtmVmNm9cBrwBeBM4AJYZ0JwJkZyJtL0vgxI3dY1jMvh5+feUAG\ncuOc60iZCBCzgOMl9ZdUBJwGDAMGmdlSgPA+MN7Gki6WNEXSlKqqqp2WabejzVsa+MMrH2/XyV5+\nrvxBN+e6iZ0eIMxsDnAT8A/geWAGUJ/E9neaWYWZVQwYMCBNuXRtmbdyI1+6/Q0mvDmf8w4fRo+8\n6Kfk3WQ4131kpJLazO4G7gaQ9EtgEbBc0mAzWyppMLAiE3lzbXt6+mJ+8sRM8vNyuOv8Ck7ZbxC5\nkg/16Vw3k6lWTAPNbIWk4cBZwFHAbsA44Mbw/nQm8uZatqmugev++h4PT15IxYhSfnfeIQwpiZ5r\n8G4ynOt+MvUcxOOS+gNbgEvNbLWkG4GJki4CFgBjM5Q3F8eHy9dz2V+m8dGKDVx60h5cecpe5OVu\nK6H0bjKc634yVcR0XJxlq4AxGciOa4WZ8eiURfzsmVn07pHHhAsP5/i9vO7HuWzgT1K7Fm2oreen\nT87kqelLOHqP/tx27sEMLPb6BeeyhQcIF9d7S9Zy2V/eYf6qjXzvU3tx6Ul7egd7zmUZDxBuO2bG\n/W/N54a/zaG0KJ+HvnUkR+zeP9PZcs5lgAeILNZaP0on7j2AW8eOon/vHjs5V865zsKHHM1iLfWj\nNHp4CX8ed5gHB+eynAeILDZ+zEhytH2AKMjN4Y6vH0qO1zc4l/U8QGSxgcU9GVVWsnU+P1ecc5g/\nCe2ci3iAyGLPzFjCpMpqmm4WvB8l51wsDxBZ6pX3V/C9R6Zz+G79OKdiGBLej5JzbjveiikLTa6s\n5pIHp7L3rn24a1wFm+samLtyo989OOe24wEiy8xeso5v3DuZIX0LmfCNwynumU9xz3zvR8k5twMv\nYsoi81Zu5Pw/T6J3jzzu/+YR7OLNWJ1zrfAAkSWWrd3M1+56m0Yz7r/oCIaGbrqdc64lHiCyQPXG\nOr5299us3bSFCRcezp4De2c6S865LsDrILq5DbX1XHjPJBZU13DfNw7nwLK+mc6Sc66LaPMOQtIg\nSXdL+nuY3y8M6uM6uc1bGrj4vinMWrKO//vKaI70Tvecc0lIpIjpXuAFYEiY/xD4broy5DpGfUMj\n4x96hzc+WcUtYw/ilP0GZTpLzrkuJpEAsYuZTQQaAcysHmhIa65cuzQ2Glc9MZMXZy/nmi/sxxcP\nKct0lpxzXVAiAWJjGD/aACQdCaxNa65cysyMXzw3h8emLuK7p4zkwmN2y3SWnHNdVCKV1N8DngH2\nkPQfYABwdlpz5VL2h1c+5u5/z+OCo8u5YszITGfHOdeFtRkgzGyapBOAvQEBH5jZlrTnzCXt/jcr\nueXFDznrkKH87PP7IXmX3c651LUZICRdCjxoZu+F+VJJ55nZ/6U9d65FLY0G17tHHjedfZCP5+Cc\na7dE6iC+ZWZrmmbMbDXwrfRlySUi3mhwAr4wajD5uf78o3Ou/RI5k+QopqxCUi5QkL4suUTEHQ0u\nL4crP7VXhnLknOtuEgkQLwATJY2RdDLwEPB8e3Yq6UpJ70maJekhST0l7SbpbUkfSXpEkgehVgws\n7skXRg3eOp+fK8b6eA7OuQ6USID4EfAycAlwKfAS8MNUdyhpKDAeqDCzA4Bc4MvATcBvzGwksBrw\np7VbYWasXF+7dd5Hg3POdbQ2A4SZNZrZ7WZ2tpl9ycz+aGbtfVAuDyiUlAcUAUuBk4HHwucTgDPb\nuY9u7fFpi3n1w5WMHl7io8E559Iikb6YjpH0D0kfSporaZ6kuanu0MwWA7cAC4gCw1pgKrAmPKUN\nsAgY2kJ+LpY0RdKUqqqqVLPRpS2sruHaZ97jiN368YevjOaw8n5+9+Cc63CJPCh3N3Al0Um83V1s\nSCoFzgB2A9YAjwKnxlnV4m1vZncCdwJUVFTEXac7a2g0vj9xBgC3njOKwSWFPhqccy4tEgkQa83s\n7x24z1OAeWZWBSDpCeBooERSXriLKAOWdOA+u427Xp/LpMpqbh07irLSokxnxznXjSVSSf2KpJsl\nHSVpdNOrHftcABwpqSg0nx0DzAZeYVsXHuOAp9uxj25pztJ13Prih3x2/105a3TcEjjnnOswidxB\nHBHeK2KWGVGlctLM7G1JjwHTgHrgHaIio78BD0u6ISy7O5X0u6va+gaufGQ6xYX5/PKsA70bDedc\n2iXSF9NJHb1TM7sGuKbZ4rnA4R29r+7i1y9+yPvL1nPPBYfRr5c/IuKcS7+EhhyV9Dlgf2BrO0oz\n+3m6MuW299bcVdz5+ly+esRwTtpnYKaz45zLEok0c70DOBe4nKi7n7HAiDTnywXrN2/h+xNnMKJf\nEVd/bt9MZ8c5l0USqaQ+2szOB1ab2XXAUcCw9GbLNbnur7NZunYTvz73YIoKErrhc865DpFIgNgU\n3mskDQG2ED3D4NLs+VlLeWzqIi47aU9GDy/NdHacc1kmkUvSZyWVADcTtTwy4K605sqxYv1mfvzE\nTA4c2pfLfWQ451wGJNKK6fow+bikZ4GeZuZjUqeRmXHV4zOpqWvgN+eO8vEdnHMZ0WKAkHSymb0s\n6aw4n2FmT6Q3a9nroUkLefn9FVz7hf3Yc2CfTGfHOZelWruDOIGom+8vxPnMAA8QaVC5ciPXPzub\n40buwvlHlWc6O865LNZigDCzayTlAH83s4k7MU9Zq76hkSsnTic/V9x89igfV9o5l1GtFm6bWSNw\n2U7KS9a747VPeGfBGm744oHs2tfHdnDOZVYitZ//kPQDScMk9Wt6pT1nWWbmorXc9s+POH3UEE4f\nNSTT2XHOuYSauX4jvF8as8yA3Ts+O9lp85YGvvvIO+zSuwfXn3FAprPjnHNAYs1c/aG4NLvx7+/z\nSdVGHrjoCPoW5Wc6O845ByTeWd8BwH5s31nffenKVHd32m9fZ/bSdTss/+Vzc3juiuMykCPnnNtR\nIp31XQP8b3idBPwKOD3N+erWRg8vIT93+xZK+bli9AjvTsM513kkUkl9NtGob8vM7EJgFNAjrbnq\n5saPGUlOswF/ciXGj9kzQzlyzrkdJdRZX2juWi+pGFiBV1C3y8Dintu1VMrPFWdXDGNgH2/a6pzr\nPBKpg5gSOuv7EzAV2ABMSmuussCAPttuwvzuwTnXGSXSiuk7YfIOSc8DxWb2bnqz1b3VNzTy1DuL\n2bW4B8vX1/rdg3OuU0qkkvppSV+R1MvMKj04tN8/Zi9nydrNfP/Te3NYeT+/e3DOdUqJ1EH8GjgW\nmC3pUUlnS/LL3Xa4541KykoLOWt0GRO/fZTfPTjnOqU2A4SZvRaKmXYH7gTOIaqodil4b8laJs2r\nZtxR5eR6Z3zOuU4s0QflCom6/T4XGA1MSGemurMJb1RSmJ/LORU+rLdzrnNrM0BIegQ4Ange+APw\namj26pJUvbGOp6YvYeyhZd6lhnOu00vkDuIe4Ctm1tARO5S0N/BIzKLdgZ8B94Xl5UAlcI6Zre6I\nfXYWD09eQF19IxccXZ7prDjnXJsSqYN4vqOCQ0jvAzM72MwOBg4FaoAngauAl8xsJPBSmO826hsa\nuf/N+RyzZ39GDvJhRJ1znV8irZjSaQzwiZnNB85gW93GBODMjOUqDV6cvZylazdzwdHeOa5zrmvI\ndID4MvBQmB5kZksBwvvAeBtIuljSFElTqqqqdlI22+/e/1QyrF8hJ+8T97Ccc67TabEOQtLo1jY0\ns2nt2bGkAqJeYX+czHZmdidRc1sqKiqsPXnYWd5bspZJldX89HP7etNW51yX0Vol9a3hvSdQAcwA\nBBwEvE308Fx7nApMM7PlYX65pMFmtlTSYLrRsxZNTVvHetNW51wX0mIRk5mdZGYnAfOB0WZWYWaH\nAocAH3fAvs9jW/ESwDPAuDA9Dni6A/aRcU1NW88aPZS+hd601TnXdSRSB7GPmc1smjGzWcDB7dmp\npCLgU8ATMYtvBD4l6aPw2Y3t2Udn8dAkb9rqnOuaEnkOYo6ku4AHAAO+Bsxpz07NrAbo32zZKqJW\nTd3GloZGHnhrPsfuuYs3bXXOdTmJ3EFcCLwHXAF8F5gdlrk2vPhe1LR1nN89OOe6oETGg9gs6Q7g\nOTP7YCfkqduY8IY3bXXOdV2JjAdxOjCdqC8mJB0s6Zl0Z6yrm7U4atrqvbY657qqRIqYrgEOB9YA\nmNl0ov6SXCu8aatzrqtLJEDUm9natOekG1m1oZanZyzhS4d601bnXNeVSCumWZK+AuRKGgmMB95I\nb7a6tocnL6SuvpFxR5VnOivypukKAAAVm0lEQVTOOZeyRO4gLgf2B2qJHmxbR9SaycXhTVudc91F\nIq2YaoCrw8u1oalp6/VnHJDprDjnXLskMqLcXsAPiCqmt65vZienL1td171vzGN4vyJO8qatzrku\nLpE6iEeBO4C7gA4bOKg7mrV4LZMrV3uvrc65biGRAFFvZrenPSfdgDdtdc51J4lUUv9V0nckDZbU\nr+mV9px1Md601TnX3SRyB9HUBfd/xywzYPeOz07X5U1bnXPdTSKtmHwQ5TZsaWjk/je9aatzrntp\nbcjRk83sZUlnxfvczJ6ItzwbvfDeMpat28wNZ3rTVudc99HaHcQJwMvAF+J8Zmw/2E9Wm/BGpTdt\ndc51Oy0GCDO7Jrz72A+t8KatzrnuKpFKaiR9jqi7jZ5Ny8zs5+nKVFdy7xuVFBV401bnXPeTyHgQ\ndwDnEvXJJGAsMCLN+eoSVm6o5ZnpSzhrtDdtdc51P4ncQRxtZgdJetfMrpN0K1lc/3Dab19n9tJ1\n2y174K0FTJu/hueuOC5DuXLOuY6XyINym8J7jaQhwBYga5u+jh5eQn7u9nUN+bli9IjSDOXIOefS\nI5EA8aykEuBmYBpQCTyczkx1ZuPHjCRH2weIXInxY/bMUI6ccy492gwQZna9ma0xs8eJ6h72MbP/\nSX/WOqeBxT0Ze2gZTSEiP1ecXTGMgX16trqdc851Na09KBf3AbnwWVY/KHf5yXvywNsLAL97cM51\nX61VUsd7QK5Jux6UC0VWdwEHhLS+AXwAPEI07kQlcI6ZrU51H+lUW29A1KTL7x6cc91Vaw/KpfMB\nud8Cz5vZ2ZIKgCLgJ8BLZnajpKuAq4AfpTEPKZtUWQ3A/kOL/e7BOddtJfIcRH9Jv5M0TdJUSb+V\n1D/VHUoqBo4H7gYwszozWwOcAUwIq00Azkx1H+k2eV41fQvzeebSY/3uwTnXbSXSiulhoAr4EnB2\nmH6kHfvcPaRxj6R3JN0lqRcwyMyWAoT3uB0bSbpY0hRJU6qqqtqRjdRNqqzmsPJScrxrDedcN5ZI\ngOgXWjLNC68bgJJ27DMPGA3cbmaHABuJipMSYmZ3mlmFmVUMGDCgHdlIzYr1m5m3ciOHlfuYSc65\n7i2RAPGKpC9Lygmvc4C/tWOfi4BFZvZ2mH+MKGAslzQYILyvaMc+0mZKZVRvfthuHiCcc91bIgHi\n28BfgNrwehj4nqT1kta1umUcZrYMWChp77BoDDAbeIZto9eNA55ONu2dYdK8anrm53DAkL6Zzopz\nzqVVIiPKpWOItMuBB0MLprnAhUTBaqKki4AFRJ0CdjqTK6s5ZFgpBXmJxFbnnOu6EmnFdFGz+VxJ\n17Rnp2Y2PdQjHGRmZ5rZajNbZWZjzGxkeK9uzz7SYd3mLcxZuo7DvXjJOZcFErkMHiPpOUmDJR0I\nvAVk5cDLU+evptHwAOGcywqJFDF9RdK5wEygBjjPzP6T9px1QpPnVZOXIw4Z3p5GXM451zUkUsQ0\nErgCeJyoC4yvSypKc746pcmV1ew/tC9FBQkNxOecc11aIkVMfwX+x8y+DZwAfARMTmuuOqHNWxqY\nsXAth5f7uA/OueyQyKXw4Wa2DsDMDLhV0jPpzVbn8+6itdQ1NPoDcs65rNHiHYSkHwKY2TpJzZuc\nprMjv05pcuigzwOEcy5btFbE9OWY6R83++yzachLp/b2vGr2GtSb0l4Fmc6Kc87tFK0FCLUwHW++\nW2toNKbNX+13D865rNJagLAWpuPNd2tzlq5jQ229P//gnMsqrVVSjwp9LQkojOl3SUBWDYIwaZ7X\nPzjnsk9rI8rl7syMdGaTK6spKy1kSElhprPinHM7jfc41wYzY3JlNYf73YNzLst4gGjD3JUbWbmh\nzsd/cM5lHQ8QbZjs9Q/OuSzlAaINkyqr6d+rgD0G9Mp0VpxzbqfyANGGyZXVVJSXImXVox/OOecB\nojXL1m5mYfUmDt+tf6az4pxzO50HiFZMCv0veQsm51w28gDRiknzVtGrIJd9B2flAHrOuSznAaIV\nk+etZvSIUvJy/WtyzmUfP/O1YE1NHR8sX+/FS865rOUBogVTKlcDeAd9zrms5QGiBZMrqynIzWHU\nsJJMZ8U55zLCA0QLJlVWc1BZX3rme5+FzrnslJEAIalS0kxJ0yVNCcv6SfqHpI/Ce2km8gZQU1fP\nzEVrvf8l51xWy+QdxElmdrCZVYT5q4CXzGwk8FKYz4jpC9ZQ32heQe2cy2qdqYjpDGBCmJ4AnJmp\njEyqrEaC0SMydhPjnHMZl6kAYcCLkqZKujgsG2RmSwHC+8B4G0q6WNIUSVOqqqrSkrnJldXsu2sx\nfQvz05K+c851BZkKEMeY2WjgVOBSSccnuqGZ3WlmFWZWMWDAgA7P2JaGRqbNX+PNW51zWS8jAcLM\nloT3FcCTwOHAckmDAcL7ikzk7b0l69i0pcHHf3DOZb2dHiAk9ZLUp2ka+DQwC3gGGBdWGwc8vbPz\nBlH/SwCH7eb1D8657JaXgX0OAp4M4yvkAX8xs+clTQYmSroIWACMzUDemDRvNeX9ixjYp2cmdu+c\nc53GTg8QZjYXGBVn+SpgzM7OT6zGRmPK/Go+te+gTGbDOec6hc7UzDXjPq7awJqaLV5B7ZxzeIDY\nzqR5YYAgDxDOOecBItbkymoG9unB8H5Fmc6Kc85lnAeIwMyYNK+aw3brR6hAd865rOYBIli0ehNL\n1272/peccy7wABFMrvT6B+eci+UBIphcWU1xzzz2HtQn01lxzrlOwQNEMGleNRXl/cjJ8foH55wD\nDxAArNxQyydVG73/Jeeci+EBApiytf7B+19yzrkmHiCI+l/qkZfDgUNLMp0V55zrNDxAEFVQHzK8\nhII8/zqcc65J1p8RN9TW896Stf78g3PONZP1AWLa/NU0Ghzmzz8459x2sj5ATK6sJjdHjB7uFdTO\nORcr6wPE2/Oq2X9IMb16ZGLsJOec67yyOkDU1jcwfeEar39wzrk4sjpAzFy0lrr6Rq9/cM65OLI6\nQEwKD8j5E9TOObejrA4Qk+dVs+fA3vTrVZDprDjnXKeTtQGiodGYMn+13z0451wLsjZAvL9sHes3\n13v/S84514KsDRCT5zV10Nc/wzlxzrnOKWMBQlKupHckPRvmd5P0tqSPJD0iKa0VA5MrVzO0pJCh\nJYXp3I1zznVZmXw67ApgDlAc5m8CfmNmD0u6A7gIuL0jd3jab19n9tJ12y0rv+pv7De4mOeuOK4j\nd+Wcc11eRu4gJJUBnwPuCvMCTgYeC6tMAM7s6P2OHl5Cfu72I8bl54rRI7wewjnnmstUEdNtwA+B\nxjDfH1hjZvVhfhEwtKN3On7MSHK0fYDIlRg/Zs+O3pVzznV5Oz1ASPo8sMLMpsYujrOqtbD9xZKm\nSJpSVVWV1L4HFvdk7KFl5IZxp/NzxdkVwxjYp2dS6TjnXDbIxB3EMcDpkiqBh4mKlm4DSiQ11YmU\nAUvibWxmd5pZhZlVDBgwIOmdjx8zkrwQIPzuwTnnWrbTA4SZ/djMysysHPgy8LKZfRV4BTg7rDYO\neDod+2+6i5DwuwfnnGtFZ3oO4kfA9yR9TFQncXe6djR+zEgOK+/ndw/OOdcKmcUt6u8SKioqbMqU\nKZnOhnPOdSmSpppZRVvrdaY7COecc52IBwjnnHNxeYBwzjkXlwcI55xzcXmAcM45F1eXbsUkqQqY\nn8Es7AKs9DQ8jU6aB0/D02jJCDNr80njLh0gMk3SlESainka2ZdGZ8iDp+FptJcXMTnnnIvLA4Rz\nzrm4PEC0z52ehqfRifPgaXga7eJ1EM455+LyOwjnnHNxeYBwzjkXlweIFEj6s6QVkma1I41hkl6R\nNEfSe5KuSCGNnpImSZoR0rguxbzkSnpH0rMpbl8paaak6ZJS6l5XUomkxyS9H76To5Lcfu+w/6bX\nOknfTSEfV4bvcpakhyQlPWCIpCvC9u8lmod4vylJ/ST9Q9JH4b3VwdNbSGNsyEejpDabRLaQxs3h\n7/KupCcllaSQxvVh++mSXpQ0JNk0Yj77gSSTtEsK+bhW0uKY38lpqeRD0uWSPgjf7a+SzMMjMfuv\nlDQ9heM4WNJbTf9zkg5vLY2UmZm/knwBxwOjgVntSGMwMDpM9wE+BPZLMg0BvcN0PvA2cGQKefke\n8Bfg2RSPpRLYpZ3f6QTgm2G6AChpR1q5wDKih4GS2W4oMA8oDPMTgQuSTOMAYBZQBOQB/wRGpvKb\nAn4FXBWmrwJuSiGNfYG9gVeBihTz8WkgL0zflGI+imOmxwN3JJtGWD4MeIHoAdlWf3Mt5ONa4AdJ\n/D3jpXFS+Lv2CPMDkz2OmM9vBX6WQh5eBE4N06cBrybzO0305XcQKTCzfwHV7UxjqZlNC9PrgTlE\nJ6hk0jAz2xBm88MrqVYHksqAzwF3JbNdR5JUTPRPcDeAmdWZ2Zp2JDkG+MTMUnnKPg8oVDT8bREt\nDH3bin2Bt8ysxszqgdeAL7a1UQu/qTOIAifh/cxk0zCzOWb2QYJ5bymNF8OxALxFNCRwsmmsi5nt\nRRu/01b+x34D/LCt7dtII2EtpHEJcKOZ1YZ1VqSSB0kCzgEeSiEPBhSH6b4k/ztNiAeITkBSOXAI\n0R1AstvmhlvUFcA/zCzZNG4j+odrTHbfMQx4UdJUSRensP3uQBVwTyjquktSr3bk58u08U8Xj5kt\nBm4BFgBLgbVm9mKSycwCjpfUX1IR0dXdsGTzEgwys6Uhb0uBgSmm05G+Afw9lQ0l/ULSQuCrwM9S\n2P50YLGZzUhl/zEuC8Vdf26r2K4FewHHSXpb0muSDksxH8cBy83soxS2/S5wc/g+bwF+nGIeWuUB\nIsMk9QYeB77b7CorIWbWYGYHE13VHS7pgCT2/XlghZlNTXa/zRxjZqOBU4FLJR2f5PZ5RLfQt5vZ\nIcBGoiKVpEkqAE4HHk1h21Kiq/bdgCFAL0lfSyYNM5tDVAzzD+B5YAZQ3+pGXYSkq4mO5cFUtjez\nq81sWNj+siT3XQRcTQqBpZnbgT2Ag4kuAm5NIY08oBQ4EvhvYGK4G0jWeaRwIRNcAlwZvs8rSdMQ\nzR4gMkhSPlFweNDMnmhPWqFI5lXgs0lsdgxwuqRK4GHgZEkPpLDvJeF9BfAkkGyF2SJgUczdz2NE\nASMVpwLTzGx5CtueAswzsyoz2wI8ARydbCJmdreZjTaz44mKBlK5QgRYLmkwQHhvsSgj3SSNAz4P\nfNVCwXc7/AX4UpLb7EEUuGeE32sZME3SrskkYmbLw0VVI/Ankv+tQvR7fSIU8U4iuvtutcK8uVCE\neRbwSAr7BxhH9PuE6GIoLZXUHiAyJFxx3A3MMbNfp5jGgKYWJZIKiU5w7ye6vZn92MzKzKycqFjm\nZTNL6opZUi9JfZqmiSo0k2rdZWbLgIWS9g6LxgCzk0kjRnuuyhYAR0oqCn+fMUR1Q0mRNDC8Dyc6\nCaSan2eITgSE96dTTKddJH0W+BFwupnVpJjGyJjZ00nidwpgZjPNbKCZlYff6yKiRh7LkszH4JjZ\nL5LkbzV4Cjg5pLcXUaOKZHtVPQV438wWpbB/iOocTgjTJ5P6RUjr0lHz3d1fRP/wS4EtRD/Ui1JI\n41iisvt3genhdVqSaRwEvBPSmEUbrSHaSOtEUmjFRFR/MCO83gOuTnH/BwNTwrE8BZSmkEYRsAro\n247v4Tqik9cs4H5CS5Uk03idKMDNAMak+psC+gMvEf3zvwT0SyGNL4bpWmA58EIKaXwMLIz5nbbV\nAileGo+H7/Rd4K/A0GTTaPZ5JW23YoqXj/uBmSEfzwCDU0ijAHggHM804ORkjwO4F/ivdvw2jgWm\nht/Y28Chqf7mW3t5VxvOOefi8iIm55xzcXmAcM45F5cHCOecc3F5gHDOOReXBwjnnHNxeYBwO13o\nifPWmPkfSLq2g9K+V9LZHZFWG/sZq6jX2VfSva9Mk/STTOfBZYYHCJcJtcBZbXXXvLNJyk1i9YuA\n75jZSenKTyfiASJLeYBwmVBPNJ7ulc0/aH4HIGlDeD8xdIw2UdKHkm6U9FVF42HMlLRHTDKnSHo9\nrPf5sH2uonENJoeO2r4dk+4rkv5C9ABV8/ycF9KfJemmsOxnRA8q3SHp5jjb/DBsM0PSjWFZU//9\nTWMqlIblr0r6jaR/hTuSwyQ9oWgMiBvCOuWKxmOYELZ/LPRNhKQxoYPDmaHzuR5heaWk6yRNC5/t\nE5b3CutNDtudEZZfEPb7fNj3r8LyG4l6t50u6cGw/d/Csc2SdG4Sf3fX1aTj6Tt/+au1F7CBqKvi\nSqKuin8AXBs+uxc4O3bd8H4isIZoHI0ewGLguvDZFcBtMds/T3TxM5LoydOewMXAT8M6PYie2t4t\npLsR2C1OPocQdb8xgKiDtpeBM8NnrxJnfAWivqDeAIrCfL/w/i5wQpj+eUx+XyWMrxCOY0nMMS4i\nepq6nOip+2PCen8O31lPoiec9wrL7yPq9JHw3V4epr8D3BWmfwl8LUyXEI1D0gu4AJgb/h49icZb\nGBb7NwjTXwL+FDOf8lPr/ur8L7+DcBlhUc+19xENHpOoyRaNo1ELfEI0aApEV/7lMetNNLNGi7pR\nngvsQ9RH1PmKukZ/m+jE29Q/0CQzmxdnf4cRDcRSZdF4CA8SjVvRmlOAeyz0WWRm1ZL6Eg2A9FpY\nZ0KzdJ6JOY73Yo5xLtu6Cl9oZv8J0w8Q3cHsTdS54IctpNvUmdtUtn0/nwauCt/Dq0TBYHj47CUz\nW2tmm4m6ChkR5/hmEt2h3STpODNb28b34bqwvExnwGW124j6srknZlk9oegzdJhXEPNZbcx0Y8x8\nI9v/lpv3H2NEo+9dbmYvxH4g6USiO4h4UunCWXH235bY42h+jE3H1dIxJZJuQ0w6Ar5kzQYRknRE\ns33HbrNtp2YfSjqUaJyL/yfpRTP7eRv5cF2U30G4jDGzaqJhPS+KWVwJHBqmzyAaJS9ZYyXlhHqJ\n3YEPiIapvERRF+tI2kttD0r0NnCCpF1CBfZ5RCPEteZF4BsxdQT9wlX2aknHhXW+nkA6zQ3XtnG6\nzwP+TdSpYLmkPZNI9wXg8hB8kXRIAvveEvO9DQFqzOwBooFqUu2W3XUBfgfhMu1Wth885k/A05Im\nEfVg2tLVfWs+IDpRDiLqMXOzpLuIilmmhZNjFW0P4blU0o+BV4iuvJ8zs1a73Daz5yUdDEyRVAc8\nR9QKaBxRpXYRUdHRhUke0xxgnKQ/EvXuens4rguBRxWNLzAZuKONdK4nunN7N3wPlUTjPLTmzrD+\nNKJiwZslNRL1LnpJksfhuhDvzdW5Tk7RkLTPmlnCowU61xG8iMk551xcfgfhnHMuLr+DcM45F5cH\nCOecc3F5gHDOOReXBwjnnHNxeYBwzjkX1/8HLn1RExrSO5kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a19648cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(variances, '-v')\n",
    "ax.xaxis.set_ticks(np.arange(1, X.shape[1], 1))\n",
    "ax.set_title('Explained variance with 0 to N components')\n",
    "ax.set_xlabel(\"Number of components\")\n",
    "ax.set_ylabel(\"Explained variance\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Inferences**:\n",
    "\n",
    "* In the above plot, we can see that the lowest cross-validation error occurs at N=6 components. \n",
    "* Also, just 6 components can explain 90% of the variance in our dataset.\n",
    "\n",
    "This shows us that we can just use 6 components in our PCR model to get a good enough result. This provides the following advantages:\n",
    "\n",
    "1. Lesser execution time\n",
    "2. Avoiding multi-collinearity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**\n",
    "\n",
    "1. \"Introduction to Statistical Learning by Tibshirani et al (2013).\" [http://www-bcf.usc.edu/~gareth/ISL/getbook.html](http://www-bcf.usc.edu/~gareth/ISL/getbook.html)\n",
    "2. R-Bloggers - [\"https://www.r-bloggers.com/performing-principal-components-regression-pcr-in-r/\"](https://www.r-bloggers.com/performing-principal-components-regression-pcr-in-r/)\n",
    "3. PLS on CRAN- [https://cran.r-project.org/web/packages/pls/index.html](https://cran.r-project.org/web/packages/pls/index.html)\n",
    "4. [Hitters data set](https://github.com/selva86/datasets/blob/master/Hitters.csv)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}