Lab04-2_107062558.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>No</th>\n",
       "      <th>year</th>\n",
       "      <th>month</th>\n",
       "      <th>day</th>\n",
       "      <th>hour</th>\n",
       "      <th>pm2.5</th>\n",
       "      <th>DEWP</th>\n",
       "      <th>TEMP</th>\n",
       "      <th>PRES</th>\n",
       "      <th>cbwd</th>\n",
       "      <th>Iws</th>\n",
       "      <th>Is</th>\n",
       "      <th>Ir</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-21</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>1021.0</td>\n",
       "      <td>NW</td>\n",
       "      <td>1.79</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-21</td>\n",
       "      <td>-12.0</td>\n",
       "      <td>1020.0</td>\n",
       "      <td>NW</td>\n",
       "      <td>4.92</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-21</td>\n",
       "      <td>-11.0</td>\n",
       "      <td>1019.0</td>\n",
       "      <td>NW</td>\n",
       "      <td>6.71</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-21</td>\n",
       "      <td>-14.0</td>\n",
       "      <td>1019.0</td>\n",
       "      <td>NW</td>\n",
       "      <td>9.84</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-20</td>\n",
       "      <td>-12.0</td>\n",
       "      <td>1018.0</td>\n",
       "      <td>NW</td>\n",
       "      <td>12.97</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   No  year  month  day  hour  pm2.5  DEWP  TEMP    PRES cbwd    Iws  Is  Ir\n",
       "0   1  2010      1    1     0    NaN   -21 -11.0  1021.0   NW   1.79   0   0\n",
       "1   2  2010      1    1     1    NaN   -21 -12.0  1020.0   NW   4.92   0   0\n",
       "2   3  2010      1    1     2    NaN   -21 -11.0  1019.0   NW   6.71   0   0\n",
       "3   4  2010      1    1     3    NaN   -21 -14.0  1019.0   NW   9.84   0   0\n",
       "4   5  2010      1    1     4    NaN   -20 -12.0  1018.0   NW  12.97   0   0"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Do not modify\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import Image\n",
    "\n",
    "# inline plotting instead of popping out\n",
    "%matplotlib inline\n",
    "\n",
    "df = pd.read_csv(\n",
    "    'http://archive.ics.uci.edu/ml/machine-learning-databases/00381/PRSA_data_2010.1.1-2014.12.31.csv',\n",
    "    sep=',')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>No</th>\n",
       "      <th>year</th>\n",
       "      <th>month</th>\n",
       "      <th>day</th>\n",
       "      <th>hour</th>\n",
       "      <th>pm2.5</th>\n",
       "      <th>DEWP</th>\n",
       "      <th>TEMP</th>\n",
       "      <th>PRES</th>\n",
       "      <th>Iws</th>\n",
       "      <th>Is</th>\n",
       "      <th>Ir</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>25</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>129.0</td>\n",
       "      <td>-16</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>1020.0</td>\n",
       "      <td>1.79</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>26</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>148.0</td>\n",
       "      <td>-15</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>1020.0</td>\n",
       "      <td>2.68</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>27</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>159.0</td>\n",
       "      <td>-11</td>\n",
       "      <td>-5.0</td>\n",
       "      <td>1021.0</td>\n",
       "      <td>3.57</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>28</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>181.0</td>\n",
       "      <td>-7</td>\n",
       "      <td>-5.0</td>\n",
       "      <td>1022.0</td>\n",
       "      <td>5.36</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>29</td>\n",
       "      <td>2010</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>138.0</td>\n",
       "      <td>-7</td>\n",
       "      <td>-5.0</td>\n",
       "      <td>1022.0</td>\n",
       "      <td>6.25</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    No  year  month  day  hour  pm2.5  DEWP  TEMP    PRES   Iws  Is  Ir\n",
       "24  25  2010      1    2     0  129.0   -16  -4.0  1020.0  1.79   0   0\n",
       "25  26  2010      1    2     1  148.0   -15  -4.0  1020.0  2.68   0   0\n",
       "26  27  2010      1    2     2  159.0   -11  -5.0  1021.0  3.57   0   0\n",
       "27  28  2010      1    2     3  181.0    -7  -5.0  1022.0  5.36   1   0\n",
       "28  29  2010      1    2     4  138.0    -7  -5.0  1022.0  6.25   2   0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Do not modify\n",
    "df = df.drop(['cbwd'], axis=1)  # drop non-scalar feature\n",
    "df = df.dropna(axis=0, how='any')  # drop samples who has nan feature\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>No</th>\n",
       "      <th>day</th>\n",
       "      <th>hour</th>\n",
       "      <th>DEWP</th>\n",
       "      <th>TEMP</th>\n",
       "      <th>PRES</th>\n",
       "      <th>Iws</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>34320</th>\n",
       "      <td>34321</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>-13</td>\n",
       "      <td>-2.0</td>\n",
       "      <td>1026.0</td>\n",
       "      <td>0.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34321</th>\n",
       "      <td>34322</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-14</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>1025.0</td>\n",
       "      <td>3.13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34322</th>\n",
       "      <td>34323</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>-13</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>1025.0</td>\n",
       "      <td>0.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34323</th>\n",
       "      <td>34324</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>-13</td>\n",
       "      <td>-4.0</td>\n",
       "      <td>1025.0</td>\n",
       "      <td>0.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34324</th>\n",
       "      <td>34325</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>-13</td>\n",
       "      <td>-3.0</td>\n",
       "      <td>1025.0</td>\n",
       "      <td>1.78</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          No  day  hour  DEWP  TEMP    PRES   Iws\n",
       "34320  34321    1     0   -13  -2.0  1026.0  0.89\n",
       "34321  34322    1     1   -14  -3.0  1025.0  3.13\n",
       "34322  34323    1     2   -13  -3.0  1025.0  0.89\n",
       "34323  34324    1     3   -13  -4.0  1025.0  0.89\n",
       "34324  34325    1     4   -13  -3.0  1025.0  1.78"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Do not modify\n",
    "idx = np.logical_or(\n",
    "    np.logical_and(df['year'].values == 2014, df['month'].values < 3),\n",
    "    np.logical_and(df['year'].values == 2013, df['month'].values == 12))\n",
    "X = df.loc[idx].drop(['pm2.5','Ir','Is','year','month'], axis=1)\n",
    "y = df.loc[idx]['pm2.5'].values\n",
    "# print(X)\n",
    "# print(X.shape)\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.preprocessing import PolynomialFeatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Training data points: 1491\n",
      "#Testing data points: 639\n"
     ]
    }
   ],
   "source": [
    "# X = df.iloc[idx].values\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.3, random_state=0)\n",
    "\n",
    "print('#Training data points: %d' % X_train.shape[0])\n",
    "print('#Testing data points: %d' % X_test.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Standardization\n",
    "sc = StandardScaler()\n",
    "X_train_std = sc.fit_transform(X_train)\n",
    "sc_x_test = StandardScaler()\n",
    "X_test_std = sc_x_test.fit_transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class LinearRegressionGD(object):\n",
    "\n",
    "    def __init__(self, eta=3.5e-10, n_iter=1000000, random_state=1,shuffle = True,lamda=1000):\n",
    "        self.eta = eta\n",
    "        self.n_iter = n_iter\n",
    "        self.random_state = random_state\n",
    "        self.shuffle = shuffle\n",
    "        self.lamda = lamda\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        rgen = np.random.RandomState(self.random_state)\n",
    "        self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1 + X.shape[1])\n",
    "        self.cost_ = []\n",
    "        if self.shuffle:\n",
    "            X, y = self._shuffle(X,y)\n",
    "        for i in range(self.n_iter):\n",
    "            output = self.net_input(X)\n",
    "            errors = (y - output)\n",
    "            self.w_[1:] += self.eta * X.T.dot(errors)\n",
    "            self.w_[0] += self.eta * errors.sum()\n",
    "            cost = (errors**2).sum() / 2.0\n",
    "            self.cost_.append(cost)\n",
    "        return self\n",
    "    \n",
    "    def net_input(self, X):\n",
    "        return np.dot(X, self.w_[1:]) + self.w_[0]\n",
    "\n",
    "    def predict(self, X):\n",
    "#         print(self.w_)\n",
    "        return self.net_input(X)\n",
    "    def sub_gradient(self, W):\n",
    "        sub_W = []\n",
    "        for w in W:\n",
    "            if w < 10**-5:\n",
    "                sub_W.append(self.lamda * np.sign(w))\n",
    "            else:\n",
    "                sub_W.append(self.lamda * np.sign(w))\n",
    "        return sub_W\n",
    "    def _shuffle(self,X,y):\n",
    "        r = np.random.permutation(len(y))\n",
    "        return X[r],y[r]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "[degree4]\n",
      "#Features: 330\n",
      "MSE train: 3080.24\n",
      "test: 3911.65454\n",
      "\n",
      "R^2 train: 0.81\n",
      "test: 0.74898\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYVNWZ7/HvS9PYIHgBnEi4H8cYGyMNdjx6MCaTOIJk\nJjoZk5hplWgmmLRm0DNxxKAmMQeHZ5wk4gWVURQn/WiYJF6eCY6jxnNiEm+IoIgaMIKCqIBGuYh0\n43v+2Lu6d1fXZVd1XXZX/z7P009XrX1bu2rXfvdae621zd0RERFJmgHVzoCIiEgmClAiIpJIClAi\nIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIA6udgXIbOXKkT5gwodrZEBHp\n155++ult7n5IIcvUfICaMGECK1asqHY2RET6NTPbWOgyquITEZFEUoASEZFEUoASEZFEqvl7UJm0\nt7ezadMm9uzZU+2sSKihoYExY8ZQX19f7ayISEL0ywC1adMmhg0bxoQJEzCzamen33N3tm/fzqZN\nm5g4cWK1syMiCdEvq/j27NnDiBEjFJwSwswYMWKESrQi0k2/DFCAglPC6PsQkXT9NkCJiEiyKUBV\nwfbt22lqaqKpqYlDDz2U0aNHd77fu3dvrHWcc845vPTSSznnueGGG2hraytFlrt56KGHOO2003LO\ns3LlSv7rv/6r5NsWkf6jXzaSKFRbG8ybB6++CuPGwfz50NJS/PpGjBjBqlWrAPj+97/P0KFD+c53\nvtNtHnfH3RkwIPM1xG233ZZ3O+eff37xmeyllStXsmbNGmbMmFG1PIhI36YSVB5tbTB7NmzcCO7B\n/9mzg/RSW79+PY2NjbS0tDBp0iS2bNnC7NmzaW5uZtKkSVx55ZWd855wwgmsWrWKjo4ODjroIObO\nncvkyZM5/vjjeeuttwC47LLLuOaaazrnnzt3LsceeyxHHHEEv//97wHYtWsXf/u3f0tjYyOnn346\nzc3NncEz6le/+hVHHHEEU6dO5d577+1Mf/zxxzn++OOZMmUK06ZNY926dbz//vtceeWVtLW10dTU\nxM9//vOM84mI5KIAlce8ebB7d/e03buD9HJ48cUXueiii1i7di2jR49mwYIFrFixgtWrV/Pggw+y\ndu3aHsu8++67fPrTn2b16tUcf/zxLFmyJOO63Z0nn3ySq6++ujPYXXfddRx66KGsXbuWyy+/nGee\neabHcrt37+a8885j+fLlPP3007z++uud04488kgeffRRnnnmGS6//HIuu+wyBg8ezBVXXEFLSwur\nVq3i9NNPzzifiEguquLL49VXC0vvrcMOO4zm5ubO93feeSe33norHR0dvP7666xdu5bGxsZuywwe\nPJhTTjkFgGOOOYZHH30047q/+MUvds6zYcMGAH77299yySWXADB58mQmTZrUY7m1a9fysY99jMMO\nOwyAlpYW7rjjDgD+9Kc/cfbZZ/Pyyy/n3K+484mIpKgElce4cYWl99b+++/f+XrdunUsXLiQX//6\n1zz77LPMmDEjY1+hQYMGdb6uq6ujo6Mj47r322+/vPMUat68eUyfPp01a9Zwzz33ZO3LFHc+EZEU\nBag85s+HIUO6pw0ZEqSX23vvvcewYcM44IAD2LJlCw888EDJtzFt2jSWLVsGwHPPPZexCrGxsZF1\n69bxyiuv4O7ceeedndPeffddRo8eDcDtt9/emT5s2DB27NiRdz4RkWwUoPJoaYHFi2H8eDAL/i9e\n3LtWfHFNnTqVxsZGPv7xj3P22Wczbdq0km/j29/+Nps3b6axsZEf/OAHNDY2cuCBB3abZ8iQIdx0\n002ccsopNDc3M2rUqM5pl1xyCRdffDFTp07F3TvTP/vZz7J69WqmTJnCz3/+86zziYhkY7V+smhu\nbvb0Bxa+8MILHHnkkVXKUbJ0dHTQ0dFBQ0MD69at4+STT2bdunUMHFj525P6XkRql5k97e7N+efs\nokYS/dzOnTv53Oc+R0dHB+7OzTffXJXgJCKSTmeifu6ggw7i6aefrnY2RER6qNo9KDMba2aPmNla\nM3vezOaE6cPN7EEzWxf+PziyzKVmtt7MXjKz6dXKu4iIlF81G0l0AP/o7o3AccD5ZtYIzAUedvfD\ngYfD94TTzgAmATOARWZWV5Wci4hI2VUtQLn7FndfGb7eAbwAjAZOBZaGsy0FUqOSngrc5e4fuPsr\nwHrg2MrmWkREKiURzczNbAIwBXgC+Ii7bwknvQF8JHw9GngtstimMC3T+mab2QozW7F169ay5FlE\nRMqr6gHKzIYCvwAudPf3otM8aANfcDt4d1/s7s3u3nzIIYeUKKelU4rHbQAsWbKEN954I+9869ev\np6mpKec8f/zjH7nrrrtib1tEpNyq2orPzOoJglObu/8yTH7TzEa5+xYzGwW8FaZvBsZGFh8TppXV\noncW0U57j/R66mk9uLWodcZ53EYcS5YsYerUqRx66KFF5SMqFaDOOOOMXq9LRKQUqtmKz4BbgRfc\n/ceRSfcBs8LXs4B7I+lnmNl+ZjYROBx4stz5zBSccqX31tKlSzn22GNpamqitbWVDz/8kI6ODs46\n6yw+8YlPcNRRR3Httdfys5/9jFWrVvGVr3wlY8nrqaee4uijj6apqYmbbrqpM/3ll1/mU5/6FFOm\nTOGYY47hiSeeAGDu3Lk88sgjNDU1ce2112adT0SkUqpZgpoGnAU8Z2apBxB9F1gALDOzrwMbgS8D\nuPvzZrYMWEvQAvB8d99X+WyXz5o1a7j77rv5/e9/z8CBA5k9ezZ33XUXhx12GNu2beO5554DgpHB\nDzroIK677jquv/76jNV3X/va11i8eDHTpk3joosu6kwfNWoUDz74IA0NDbz44ovMmjWLJ554ggUL\nFnD99ddzzz33AMEjNjLNJyJSKVULUO7+W8CyTP5clmXmAxUYprU6HnroIZ566qnOx228//77jB07\nlunTp/PSSy/xD//wD3z+85/n5JNPzrmebdu28f7773eO3XfWWWfxyCOPAPDBBx9wwQUXsHr1agYO\nHJj18Rdx5xMRKReNJJEg7s65557LD3/4wx7Tnn32We6//35uuOEGfvGLX7B48eKitvGjH/2IsWPH\n8tOf/pT29naGDh3aq/lERMql6q34pMtJJ53EsmXL2LZtGxC09nv11VfZunUr7s6XvvQlrrzySlau\nXAn0fKRFysiRIxk8eDCPPfYYAG2R59O/++67jBo1CjNj6dKlnSOLZ3o8Rqb5REQqRQEqj3rqC0rv\njU984hN873vf46STTuLoo4/m5JNP5s033+S1117jxBNPpKmpiXPOOYerrroKgHPOOYe///u/z9hI\n4rbbbuO8886jqamJAQO6vuYLLriAW265hcmTJ/PKK690PsRwypQp7Nu3j8mTJ3PttddmnU9EpFL0\nuA1JDH0vIrWrmMdtqAQlIiKJpAAlIiKJ1G8DVK1XbfY1+j5EJF2/DFANDQ1s375dJ8WEcHe2b99O\nQ0NDtbMiIgnSL/tBjRkzhk2bNqGRzpOjoaGBMWPGVDsbIpIg/TJA1dfXM3HixGpnQ0REcuiXVXwi\nIpJ8ClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJI\nClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIVQ1QZrbEzN4y\nszWRtOFm9qCZrQv/HxyZdqmZrTezl8xsenVyLSIilVDtEtTtwIy0tLnAw+5+OPBw+B4zawTOACaF\nyywys7rKZVVERCqpqgHK3X8DvJ2WfCqwNHy9FDgtkn6Xu3/g7q8A64FjK5JRERGpuGqXoDL5iLtv\nCV+/AXwkfD0aeC0y36YwrQczm21mK8xsxdatW8uXUxERKZskBqhO7u6AF7HcYndvdvfmQw45pAw5\nExGRcktigHrTzEYBhP/fCtM3A2Mj840J00REpAYlMUDdB8wKX88C7o2kn2Fm+5nZROBw4Mkq5E9E\nRCpgYDU3bmZ3Ap8BRprZJuB7wAJgmZl9HdgIfBnA3Z83s2XAWqADON/d91Ul4yIiUnZVDVDu/tUs\nkz6XZf75wPzy5UhERJIiiVV8IiIiClAiIpJMClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJI\nClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIlKwtjaYMAEGDAj+t7WVfhtVHc1cRESS\np60N5s2DV1+FceNg/nxoaek+ffZs2L07eL9xY/B+28mLYGB7j/XVU19UPlSCEhGpUcWUclLBZ+NG\ncA/+n3MOjBzZtZ45c7qCU8ru3WQMTgDttAMjhxeafwUoEZEalCnQzJ7dPUhlCmDz5vUMPu3tsH17\n13q2by8mR+PGF7qEqvhERCooV/VZpmm/+x3cfDN8+GEwz/77B++jVW6ZzJsHV7ywiIZh3Us124BF\n79Rz4PLWjNV06cGpdKzgApEClIhIaNE7i8LqqO7qqaf14NbO921t8PqnFzFo//zzRmW6d3PmmXDe\neWAGO3d2zbtxI8yaBfv2dV/Hrl1BOvS8LxQNbhs30iM4pbTT3q2absHGnoEMYM+OeuaOb409T6kp\nQIlIzclUEoHcN/6BjMEpPb21FW66CX6yPf+86TJVnxVz4t+3L1hXKv8Lty6Cme1cNDPrpnu4eOWi\nzvVnC2TR9DjzlJoClIgkUr6WZLmWSy+lnHNOUELZu7crbdOJi1j4TmEn17a2IDi5559vzpyuezUj\nRkBTU7BdyB6UovJN37gxuHc0bhxc9EzhQaIcgWXPjvqswbYYClAikjjZmjFD9yC1cGuGZs0zYf5r\nQUBKFy2VNAwt/AQ9b17+4ARBQGyPrH77dnj44a73cYPDNW8v7JEW3YdUo4WkyF3Vt6Tg9SlAiUjR\noqWc4WEj4rffzn3zP32+mTNh4nfT7ufMhKs2dT8Z794dlEq63fvJcgbLFJygd6WGCRPiB4P28tV6\n0TCsPWPgKkap1pPJkCHB97BrV/HrUIASqTHZbvTTUc9PPtlacJVZNumlnGjT41SJ53e/g2XLuk9L\nn+/GG+Ga+fHub2zfTsaGCZWQpJJKJV3z9sLYVXQDBgSlumgLxBtvLH7bClAiZZDr/knclmL5ZA1E\n2Qxs7zzJbtwIZ50VnEAWLcqf50zTd+7M3SR59+7sJ6c492CSpNh7KLUi7nd13nndj6elS3u3XQUo\nkRJra4Nzz+1+Q/7cc4PXLS25W4qZQV0dHHEEvPhiz74v784sMCilWbBxUbf7FzfeGPyNGAE7dsCV\n67oCxzZg4Tvhgh317Dkx3s34OM2Okx6cLhw+J+882RoEQHmrzgqRCqxxPu8Lh8/pdb6XL+96nanF\nYqEUoKRfKKR/S6HNk3useyb8yxvdT9R79wb9Xc48E655O3de9+2DtWu73qdKG9sK3uuesp2oUtVu\nWU9kA9tpGBp/G+knut70lYmuq5Qlmd62OJs7vrXkJcHU5zR+PFz0TOHBIl9gLXfgfPXVzK+LpQAl\niRe3uXGu+XKVWha+E/nRzoSLZnadKM48s/v8qaqxM8+E8eODbbTPjHf/JM7JrBJX3te8vRD37A0J\nyqFUJ/HericafOIGzFz9lEoZnKLBpRT3u6pRjTpuXPfXvd0PBSgpq97eb8nU3Piss+C222D67T1b\nfkWDS6ZmyXGllwLSm/am8jJ7dtDaLO46k6KSwanUig0M6aWLoUO7j9yQTTU6qBYjveRXSP5KUTId\nMqSrxgGC170dOilWgDKzacAqd99lZmcCU4GF7l7xdi1mNgNYCNQBt7j7gkrnQeJpbYUjsrTO6lFy\nSemoZ+8HkZZaM+MHgJRocOl2H6UXsv3YyzdumWST+i727Kxn7rjuFzmFlEB37YofpPqCYqpQ04N2\nvmrPbNNGjICFC7tfDKZep2o13Nv3Fpq/uCWoG4HJZjYZ+EfgFuAO4NOFbrA3zKwOuAH4S2AT8JSZ\n3efua3MvKcVKDesS7ZyY6WBM9+Mti7IGp5wGtjNI5fqaNmJEaarHGoa2M3588dVI48bBRU9lfn5R\nOceXi2PBxkVV2W5qn4cODQJ49Hc/ZAgsXhy8Tq9K35blBmlLS9d5wuzZ5wrNj3mMbtFmttLdp5rZ\nFcBmd781lVboBnvDzI4Hvu/u08P3lwK4+z9nW6bZzFdUKH8iIpKZwdPu3lzIMnGHP98RBoMzgV+Z\n2QAo8hGJvTMaeC3yflOY1o2ZzTazFWam2CQikgiFP7AwbmXKV4C/A77u7m+Y2Tjg6kI3VinuvhhY\nDNDc3OysKCxOZbuxv3dXPZeMay2oF37G+yw5ROuE998/KGJH73OkitnRIWS2zcy+jZ9MmdNZFN++\nPXd9e676+2zNV5PS3yOO1D5UMs/Rz62Srar27Kjn6qmtfPe5RdQ19NxmMa34Ri6f0znGXL7PME4j\nmFSry2gVXaHfzZyDuz7ffL+F6LwpuX6fqe+uN4+ZGD8+df+l5zSzrn5uKW1twaM0frS18N8iBH3o\nli4Nzg9xOnJfOHxOj3NKajSIOOrqoKOje1r2IaEm9yhM5BMrQLn7G8CPI+9fJbgHVWmbgbGR92PC\ntNjSmyLPnBl0LovWp2ZrNjxo//bOwRnTmxr3ZsiYqGhHykxjWO3eDZs+FRmFOc/w+tGRA6S65o5v\njXUCTgUXCC4qci2T6aQLwMFwyTaA1h4ja+eSa1stLfHWM3L5nJL9HuJKtfYstEFNHMXei0q1aksP\nwinDM5QnUp9bMX3e0gNN68GtOQPwnh31Gc9fw4fHf2Ju+rOqIFf/p/pB8dbaJWeAMrMdQKZYaoC7\n+wGFbrCXngION7OJBIHpDIKSXSyZmixHh2KJ02w40w84Y0uxjnrmHFL4gR3nCruQq/BSdphMKWdJ\noDdNiPOd/ONenZeqpJXedHdAzAr1hmHtvP02/Pu/d41A0RupG9W5Ls4ynSwzeTtPJ+PU9lJyjQv4\n3dmtPVpBFvL910fuMqRGLci2fH2F70ikn/jTRzeHYOSOtraeF7ctLfFano4YETRmyNU/sJ76rN08\n5oxr5ZIN8fcpk7q6rtep7/onacGt65xT4lZ87j6s0BWWk7t3mNkFwAMEzcyXuPvzcZePM/RGyZoN\nDwyGrck3akCllSKwlOs5MqnAWc1qw2L7g+TrwZ+6uo17ZTxuXGFX03E6M0dbVGWT78RYaOfLrFVM\nA9sz/tZSx0C0xVh6X5r0kgLAt36TvRruknHZL8iyncDpqGfIkMLPB5nylq3kuXdv94cOxslX6vgc\nMiR/S1qgoLEdU+JchKSk+hpC9u+6YVh7+Fm+XlBtFxTYUdfM/gxoSL0Pq/oqyt2XA8vzzhiKftjZ\nip5JHrjSLH59cG/FGfpl/PjSb/efDp3TOW5dtZWyaXFdXXCPIRow4vbJmjkzqMt/9VV6XJGmi/vs\npDhyXXFDV+fLcoqWPiZM6Bkkdu/ueWIvtjNtrhP4yMU9S5y5HlaYq7o/20k/2zkpPV/RC5BS31ZI\nl+0iZOhQeP/9oFqvrg4+85mgBN710MTs65w1C268cVvBl+txO+p+AfgR8FHgLWA88AIwqdANVtrG\njV3F6GwffJKD0ze/GRwEhVy1FjvoY5yT8/z5ua/oi2lMUYrgVF/f+741velNn361nelKGnJcsae5\n5ZauKqFc+7XwnYVZn52U7eo8l3xX3Kn1vb6rPuNjL3pblTZ+PGzY0PU+2wm8FOO85ZOtxJkepDJ9\n1+kl2mz3daJDAxWTl3LINALEkCHBfkcbZ6VfFOWyPHaRoru4JagfAscBD7n7FDP7C4Im54n34Ydd\nP9RSDL1RKangdNT8Iju8lsGIEYWVAqLLFaKYQHPbbTDvxO7POsrVoivOaNVx1VPP4sXxxguMBoBc\nN7Cj9ytSQWfEiKBDZLbl0j+zcp3Eg/3qfUkzU1CPDpUD2S8q457YeytT1em0afkfS5J+8h40KLiI\nin6vmfa3VHnsTTBLHwEi0zoLHam82GMxboBqd/ftZjbAzAa4+yNmdk1xm6y81IeT+oDTBwDNp9LN\nk6NF+M7WelWWqvPOJ/1HWF8fLFdIq6TUCTm9ejNbVWyxV+29aT7co/VcEVe4+e4zpIvbsiqlUifx\nYsUJ6tmu5ktxYs8nW9XprFm5l8t08t67N6gi++hHSxdIcuUReh+kci1faMApduDYuAHqT2Y2FPgN\n0GZmbwG9eJBvZQ0f3lWfP25ccCVa6I+9vh4+2FnPfkMLCxjFlAZ6c+CmTm69fZQAdAWI9DrvbJ/D\nBzvrg5JMhpPOwq31GYeU2bcn883oESPgy18O+nSkpmUKGmZBa7dMP9KrNmTeZuoGeL77Flm/u47S\ntAjLVp1mWVrULdi4KPYFS7aTeDFX2725Qs91Tyt1Eky1/kpvDVtPPa0twWdUyhJCXJkCze7d3av4\nMgWEbCfvnTu7V5NB9s827meeLY/FVO8WIlPAyXrO2VnP/PmFFwwg/lBH+wN7CJqXtwAHAm3uXuBp\nvvIGDGj2+voV3e5zDBoUdC5LdZLLVzL64eFz2LGj615JnEYV0U69O3f2DIi5rt5vPLG1sx4+VzXQ\nd8fMyVvMTtXpZ2vum15aGDEiuKGb60cR7ciXq1ST60Zv3B9jps6c6fsHmaen389Iz0uuKsCoPTvq\nuXRCYR20e2PkyMwXUHFK8BeNmJM1n+lX25D9XllvlilUrmM8az+vNLmO7xtPbC3qeyukw2r0WMve\nUbX7fNk+21mzul+YpdIzfebZ8pipE3Ap5TouIPPv2cwKHuooVoDqywYObPZ9+3qOJDFgQNcXmC/g\n/GTKnIJPgOnbKuRjjh5cuX68I5fP6TwQhg+nWxCF7gd1nB7zcfcn+gMsqhNpDrmCWLYfxFlnFfcj\nLWSUj2L2pVjpT+SF4KLqX97In99c+cx24sz1vRezTKFKEaCiShVUcwWadNFjra0te2khOl+29dfV\nZe4Am+kzr8T3k02hJetiAlSsroNmtsPM3gv/9pjZPjN7r5ANVcuhjVu55u2FPf6ueqVrtOB89xt6\n25Ko0HsBhbTs2bAhOOC3bYMlS4ID0yz4X8gPspB6/fnzg/lLLXVi2biRzhE7Zs/uaoW5eHHm/cv2\neSX5HkxbW3ByGTAg+N/W1jWtpaXnd7lkSf517t2Vu+qxmOO4mq3oipWr2qsQmY7zbMNDRY+1lpbs\nDYOi82X7DDMFp2zzZ8pjvt/yoncWsfCdhT3+Fr1T2AjqqYZn48YFeZs3r/txXApxhzrq7LBrZgac\nStCqL/EG1GUuuqSXmHLds+ltS6JC6l/TD658/VKiim2KWmi/imgrn1J6/dOLuGpTz319fVc90Jp1\n/6p5I70YcW5sZ9rXXK0nLxw+BzO4OEeJsZjjuNqt6IpRqqCaqTXbzJmZq9/Sj7WFC3Mfk21twcVJ\npmCUrQSV6TOP0+IuXa6nSxeiXA00ogp+8o4HdYL3mNn3gLmlyUZ1ZKqeSr8nkzp59+YEmGsMs3zD\nlRTTE7xQxVQFpE6gpXgYYEqmfjW50qN5gcJvpMftk1Rq5bqxnS9oFHMc97XgD6UNqpkuFPI1M08t\nB7mrqzMFoVz3oLJ95pXsIxVViQYacTvqfjHydgDQTNBoouZES1Gpg6LYE2BUtiuqOMOVSH7F/EjT\ng3+hI88Xq9gr/FzN0uMEjWKO41Ic+/kUUksQR7mDatxjLdt82foQ1dV1VVvHCYLVVonq37glqL+O\nvO4ANhBU89Uks54HRW+vUirxQ8+l1CeBTh3Zm3L3NWX7jNIUe4UfDajpN6jj3m8s5jgu9xV6qWsJ\nCvmtlbqTaxzZTuAffli6800lVKL6N+49qHNKt8nkK1fzzPSDLrhZmfmEWOofbbmqCn/yydasrYjm\nbCjLJsvmwOWtFTlZleIKvy+cwKopzudTiXsomfTF+3qZVKL6N2crPjO7zsyuzfZXumyUz4f7Cnwq\nWwWV6mZlNfXFVl6Z5GpBWGq5WiRKeWRqNVmq1n6FKqblXSllqxEotKagEsdxzn5QZpYa1GMa0Aj8\nLHz/JWCtu3+zdFkpD7Nmh579oErdf6cYpe7/UUpxqz5K1Q+jrQ02nbiIhgwjVJSjRJmumv1JpLyy\n9YvK1sm93J1cU3lK+j2mUiumH1S+50EtDVf8LeAEd+8I398EPFpsRitp0KDMo2XvLdNozLWgkKqP\nUhXzg9EiegahSgWIWikJSk/ZSkqFNOcuNVXRxhO3kcTBwAFA6nkeQ8O0xBs9Gt58s+cJ9KP/r1UH\nSBaFNB8tVeOPageIWrkvID3l6hAbZ0R1qZ6YD6FmAfCMmd1uZkuBlcBV5ctW6Qwfrvr+QhUaLKIj\nWmzYUNxnW+3RIKp9X0DKJ9sxlDoX6NyQXLEClLvfBvxP4G7gl8Dxqeq/vqAUJ9ByKNXNylKrRrCo\ndoBQw4XalevYSuq5QQL5Gkl83N1fNLOpmaa7+8qy5axEmpubfcWKno0kJLtKjGCdbbv97caxVIaO\nreor+WjmZrbY3Web2SMZJru7f7bQTFaaAlRxauUHXSv7IdLX6XEbGShA9V/VKgmKSE/lfNzGl8xs\nWPj6MjP7pZlNKSaTIpVSrY6Y0vflehxKf1PNzyJuK77L3X2HmZ0AnATcCtxUvmyJ9F61m65L31TJ\nUUWSrtqfRdwAlerO9nlgsbv/ChhUniyJlEa1m65L36SSd5dqfxZxA9RmM7sZ+Aqw3Mz2K2BZkaqo\ndtN16ZtU8u5S7c8ibpD5MvAAMN3d/wQMBy4uW65ESkB9m6QYKnl3qfZnEbej7m7gLeCEMKkDWFeu\nTImUijpiSqFU8u5S7c8ibiu+7wGXAJeGSfXAT4vdaNgq8Hkz+9DMmtOmXWpm683sJTObHkk/xsye\nC6dda2bJfY6GiPRZlSh595VWgtWuhYjVD8rMVgFTgJXuPiVMe9bdjy5qo2ZHAh8CNwPfcfcVYXoj\ncCdwLPBR4CHgY+6+z8yeBP4BeAJYDlzr7vfn25b6QYlIkvTX/nll6wcF7PUgknm4of0LzVyUu7/g\n7i9lmHTH+zSCAAANlklEQVQqcJe7f+DurwDrgWPNbBRwgLs/HubjDuC03uRBRKQaqt0yri+JG6CW\nha34DjKzbxCUbG4pQ35GA69F3m8K00aHr9PTRUT6lGq3jOtLYj0Pyt3/1cz+EngPOAK4wt0fzLWM\nmT0EHJph0jx3v7fgnBbAzGYDswHG9cemNyKSWHr2WHxxH1hIGJAeBDCzAWbW4u5Zb+25+0lF5Gcz\nMDbyfkyYtjl8nZ6ebduLgcUQ3IMqIh8iImVRqqdQ9wc5q/jM7ICwVd31ZnayBS4A/kjQN6rU7gPO\nMLP9zGwicDjwpLtvAd4zs+PC1ntnA2UthYmIlEO1W8b1Jfket3Ev8A7wGPA54M8AA+a4+6qiN2r2\nN8B1wCHAn4BV7j49nDYPOJegr9WFqZZ6YXP024HBwP3Atz1GE0S14hMRqb5yPA/qOXf/RPi6DtgC\njHP3Pb3KaQUpQImIVF85mpm3p164+z5gU18KTiIi0nflayQx2czeC18bMDh8bwRP1D2grLkTEZF+\nK2eAcve6SmVEREQkSo/MEBGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKA\nEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGR\nRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRFKAEhGRRKpKgDKzq83s\nRTN71szuNrODItMuNbP1ZvaSmU2PpB9jZs+F0641M6tG3kVEpDKqVYJ6EDjK3Y8G/gBcCmBmjcAZ\nwCRgBrDIzOrCZW4EvgEcHv7NqHSmRUSkcqoSoNz9v929I3z7ODAmfH0qcJe7f+DurwDrgWPNbBRw\ngLs/7u4O3AGcVvGMi4hIxSThHtS5wP3h69HAa5Fpm8K00eHr9PSMzGy2ma0wsxVbt24tcXZFRKQS\nBpZrxWb2EHBohknz3P3ecJ55QAfQVsptu/tiYDFAc3Ozl3LdIiJSGWULUO5+Uq7pZvY14K+Az4XV\ndgCbgbGR2caEaZvpqgaMpouISI2qViu+GcA/AV9w992RSfcBZ5jZfmY2kaAxxJPuvgV4z8yOC1vv\nnQ3cW/GMi4hIxZStBJXH9cB+wINha/HH3f2b7v68mS0D1hJU/Z3v7vvCZVqB24HBBPes7u+xVhER\nqRlVCVDu/uc5ps0H5mdIXwEcVc58iYhIciShFZ+IiEgPClAiIpJIClAiIpJIClAiIpJIClAiIpJI\nClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAi\nIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJIClAiIpJI\nClAiIpJIClAiIpJIClAiIpJIClAiIpJIVQlQZvZDM3vWzFaZ2X+b2Ucj0y41s/Vm9pKZTY+kH2Nm\nz4XTrjUzq0beRUSkMqpVgrra3Y929ybgP4ErAMysETgDmATMABaZWV24zI3AN4DDw78ZFc+1iIhU\nTFUClLu/F3m7P+Dh61OBu9z9A3d/BVgPHGtmo4AD3P1xd3fgDuC0imZaREQqamC1Nmxm84GzgXeB\nvwiTRwOPR2bbFKa1h6/T07OtezYwG2DcuHGly7SIiFRM2UpQZvaQma3J8HcqgLvPc/exQBtwQSm3\n7e6L3b3Z3ZsPOeSQUq5aREQqpGwlKHc/KeasbcBy4HvAZmBsZNqYMG1z+Do9XUREalS1WvEdHnl7\nKvBi+Po+4Awz28/MJhI0hnjS3bcA75nZcWHrvbOBeyuaaRERqahq3YNaYGZHAB8CG4FvArj782a2\nDFgLdADnu/u+cJlW4HZgMHB/+CciIjXKgkZxtau5udlXrFhR7WyIiPRrZva0uzcXsoxGkhARkURS\ngBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIR\nkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkURSgBIRkUSq+Ue+m9kO4KVq56NC\nRgLbqp2JCtG+1q7+tL/9aV+PcPdhhSwwsFw5SZCX3L252pmoBDNboX2tPf1pX6F/7W9/29dCl1EV\nn4iIJJIClIiIJFJ/CFCLq52BCtK+1qb+tK/Qv/ZX+5pDzTeSEBGRvqk/lKBERKQPUoASEZFEqskA\nZWZfMrPnzexDM2tOm3apma03s5fMbHq18lhKZjYj3J/1Zja32vkpNTNbYmZvmdmaSNpwM3vQzNaF\n/w+uZh5LxczGmtkjZrY2PIbnhOk1t79m1mBmT5rZ6nBffxCm19y+pphZnZk9Y2b/Gb6v5X3dYGbP\nmdmqVBPzQve3JgMUsAb4IvCbaKKZNQJnAJOAGcAiM6urfPZKJ8z/DcApQCPw1XA/a8ntBN9X1Fzg\nYXc/HHg4fF8LOoB/dPdG4Djg/PD7rMX9/QD4rLtPBpqAGWZ2HLW5rylzgBci72t5XwH+wt2bIn29\nCtrfmgxQ7v6Cu2caPeJU4C53/8DdXwHWA8dWNncldyyw3t3/6O57gbsI9rNmuPtvgLfTkk8Floav\nlwKnVTRTZeLuW9x9Zfh6B8HJbDQ1uL8e2Bm+rQ//nBrcVwAzGwN8HrglklyT+5pDQftbkwEqh9HA\na5H3m8K0vqwW9ymOj7j7lvD1G8BHqpmZcjCzCcAU4AlqdH/DKq9VwFvAg+5es/sKXAP8E/BhJK1W\n9xWCi42HzOxpM5sdphW0v312qCMzewg4NMOkee5+b6XzI9Xj7m5mNdVfwsyGAr8ALnT398ysc1ot\n7a+77wOazOwg4G4zOyptek3sq5n9FfCWuz9tZp/JNE+t7GvECe6+2cz+DHjQzF6MToyzv302QLn7\nSUUsthkYG3k/Jkzry2pxn+J408xGufsWMxtFcAVeE8ysniA4tbn7L8Pkmt1fAHf/k5k9QnCvsRb3\ndRrwBTObCTQAB5jZT6nNfQXA3TeH/98ys7sJbkcUtL/9rYrvPuAMM9vPzCYChwNPVjlPvfUUcLiZ\nTTSzQQSNQO6rcp4q4T5gVvh6FlATpWYLikq3Ai+4+48jk2puf83skLDkhJkNBv4SeJEa3Fd3v9Td\nx7j7BILf6K/d/UxqcF8BzGx/MxuWeg2cTNB4rbD9dfea+wP+huBezAfAm8ADkWnzgJcJHsFxSrXz\nWqL9nQn8IdyvedXOTxn2705gC9Aefq9fB0YQtAJaBzwEDK92Pku0rycQ1N0/C6wK/2bW4v4CRwPP\nhPu6BrgiTK+5fU3b788A/1nL+wr8D2B1+Pd86rxU6P5qqCMREUmk/lbFJyIifYQClIiIJJIClIiI\nJJIClIiIJJIClIiIJJIClPR7ZrYvHHF5jZn9h5kN6cW6PhMZqfoLuUaXN7ODzKy1iG1838y+U2we\nS70ekXJRgBKB9z0YcfkoYC/wzehECxT8W3H3+9x9QY5ZDgIKDlAi/YUClEh3jwJ/bmYTwmds3UHQ\niXSsmZ1sZo+Z2cqwpDUUOp/H9aKZrSR4zAth+tfM7Prw9UfM7O7w2Uerzex/AQuAw8LS29XhfBeb\n2VNm9mzq+Uhh+jwz+4OZ/RY4Ij3TZnagmW1MBdKwJ/9rZlZvZt8I17nazH6RqYRoZv/XwmenmdlI\nM9sQvq4zs6sjeTovTB9lZr+JlDw/VYoPXyRKAUokZGYDCZ6r9VyYdDiwyN0nAbuAy4CT3H0qsAL4\n32bWAPwb8NfAMWQewBjgWuD/efDso6kEvevnAi+HpbeLzezkcJvHEjwf6RgzO9HMjiEYHqeJYFSJ\nT6av3N3fJRh14tNh0l8RjKDSDvzS3T8ZbvsFgpE44vo68K67fzLc7jfCYcL+Llx/EzA53LZISfXZ\nwWJFSmhw+MgHCEpQtwIfBTa6++Nh+nEED4T8XTiy+CDgMeDjwCvuvg4gHAB0Nj19FjgbOkfwftd6\nPk305PDvmfD9UIKANQy42913h9vINtbiz4CvAI8QBLRFYfpRZvZ/CKoUhwIP5PowMuTpaDM7PXx/\nYJinp4Al4cC297i7ApSUnAKUSHgPKpoQBqFd0SSC5xV9NW2+bsv1kgH/7O43p23jwpjL3wdcZWbD\nCUpzvw7TbwdOc/fVZvY1grHg0nXQVaPSkJanb7t7j6BmZicSPIDvdjP7sbvfETOfIrGoik8knseB\naWb259B5j+djBKNvTzCzw8L5vppl+YeBb4XL1pnZgcAOgtJRygPAuZF7W6MteJbOb4DTzGxwOEL0\nX2fagAdPp30KWEgwGOm+cNIwYEtY2mnJkr8NBEEN4PRI+gPAt8JlMbOPhfs+HnjT3f+N4AmxU7Os\nV6RoKkGJxODuW8PSx51mtl+YfJm7/8GCp4X+ysx2E1QRDsuwijnAYjP7OrAP+Ja7P2ZmvzOzNcD9\n4X2oI4HHwhLcTuBMd19pZj8jGBn6LYIglM3PgP+geynpcoKn8m4N/2fK378Cy1L7Ekm/BZgArLQg\nU1sJHtP9GeBiM2sP83l2jjyJFEWjmYuISCKpik9ERBJJAUpERBJJAUpERBJJAUpERBJJAUpERBJJ\nAUpERBJJAUpERBLp/wNt1ff3IJsuGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff4775e7160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "regr = LinearRegressionGD(9.45e-9,2000000)\n",
    "X_lin = X_train_std\n",
    "# Create quadratic features\n",
    "# quadratic = PolynomialFeatures(degree=3)\n",
    "cubic = PolynomialFeatures(degree=4)\n",
    "# X_quad = quadratic.fit_transform(X_lin)\n",
    "X_cubic = cubic.fit_transform(X_lin)\n",
    "\n",
    "print('\\n[degree4]')\n",
    "X_cubic_train = cubic.fit_transform(X_train_std)\n",
    "X_cubic_test = cubic.fit_transform(X_test_std)\n",
    "print('#Features: %d' % X_cubic_train.shape[1])\n",
    "regr = regr.fit(X_cubic_train, y_train)\n",
    "y_train_pred = regr.predict(X_cubic_train)\n",
    "y_test_pred = regr.predict(X_cubic_test)\n",
    "print('MSE train: %.2f\\ntest: %.5f\\n' %\n",
    "      (mean_squared_error(y_train, y_train_pred),\n",
    "       mean_squared_error(y_test, y_test_pred)))\n",
    "print('R^2 train: %.2f\\ntest: %.5f\\n' % (r2_score(y_train, y_train_pred),\n",
    "                                       r2_score(y_test, y_test_pred)))\n",
    "plt.scatter(\n",
    "    y_train_pred,\n",
    "    y_train_pred - y_train,\n",
    "    c='blue',\n",
    "    marker='o',\n",
    "    label='Training data')\n",
    "plt.scatter(\n",
    "    y_test_pred,\n",
    "    y_test_pred - y_test,\n",
    "    c='lightgreen',\n",
    "    marker='s',\n",
    "    label='Test data')\n",
    "plt.xlabel('Predicted values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.legend(loc='upper left')\n",
    "plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='red')\n",
    "plt.xlim([-10, 50])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE train: 201.89, test: 1489.11\n",
      "R^2 train: 0.99, test: 0.90\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUVOWZ7/HvjzuERhLMKOHqGGNsVBrs8ejBmJuDaGai\nk9GJCV6iJjhRk0xmYoRB4yVDcC0nWfEaw3g/9lGZJF7WRONB4xxN4g0QFFECRkEUFcmJTUSBxuf8\nsXc3RXd1dVV3V9Wm+vdZi9VV7749L9W9n3r3fvf7KiIwMzPLmn7VDsDMzCwfJygzM8skJygzM8sk\nJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8ukAdUOoNz23HPPmDhxYrXDMDPr05Ys\nWfJWRHy4lG1qPkFNnDiRxYsXVzsMM7M+TdLaUrfxJT4zM8skJygzM8skJygzM8ukmr8Hlc/27dtZ\nv3497733XrVDsdSQIUMYO3YsAwcOrHYoZpYRfTJBrV+/nrq6OiZOnIikaofT50UEmzZtYv369eyz\nzz7VDsfMMqJPXuJ77733GDVqlJNTRkhi1KhRbtGa2S76ZIICnJwyxp+HmbXXZxOUmZllmxNUFWza\ntImGhgYaGhrYe++9GTNmTNv7bdu2FbWP008/nVWrVhVc55prrqGpqak3Qt7Fgw8+yPHHH19wnaVL\nl/KrX/2q149tZn1Hn+wkUaqmJpg7F9atg/HjYd48mDmz+/sbNWoUy5YtA+Diiy9m+PDhfOc739ll\nnYggIujXL/93iJtuuqnL45xzzjndD7KHli5dyooVK5gxY0bVYjCz3ZtbUF1oaoJZs2DtWohIfs6a\nlZT3tjVr1lBfX8/MmTOZNGkSGzZsYNasWTQ2NjJp0iQuvfTStnWPOOIIli1bRktLCyNHjmT27NlM\nnjyZww8/nDfffBOACy64gB//+Mdt68+ePZtDDz2U/fffn9/97ncAvPPOO/z93/899fX1nHDCCTQ2\nNrYlz1y//OUv2X///Zk6dSr33HNPW/njjz/O4YcfzpQpU5g2bRqrV6/m3Xff5dJLL6WpqYmGhgZ+\n9rOf5V3PzKwQJ6guzJ0LW7bsWrZlS1JeDi+88ALf/va3WblyJWPGjOGyyy5j8eLFLF++nEWLFrFy\n5coO27z99tt88pOfZPny5Rx++OHceOONefcdETz55JNcfvnlbcnuqquuYu+992blypVceOGFPP30\n0x2227JlC2eddRb33XcfS5Ys4bXXXmtbdsABB/Doo4/y9NNPc+GFF3LBBRcwdOhQvve97zFz5kyW\nLVvGCSeckHc9M7NCfImvC+vWlVbeU/vuuy+NjY1t72+//XZuuOEGWlpaeO2111i5ciX19fW7bDN0\n6FCOOeYYAA455BAeffTRvPv+whe+0LbOyy+/DMBvfvMbzj//fAAmT57MpEmTOmy3cuVKPvaxj7Hv\nvvsCMHPmTG699VYA/vSnP3Hqqafy4osvFqxXseuZmbVyC6oL48eXVt5TH/jAB9per169miuuuIJf\n//rXPPPMM8yYMSPvs0KDBg1qe92/f39aWlry7nvw4MFdrlOquXPncvTRR7NixQruvvvuTp9lKnY9\nM7NWVUtQksZJeljSSknPSfpWWv4hSYskrU5/fjBnmzmS1khaJenoSsQ5bx4MG7Zr2bBhSXm5NTc3\nU1dXx4gRI9iwYQMPPPBArx9j2rRpLFy4EIBnn3027yXE+vp6Vq9ezUsvvUREcPvtt7cte/vttxkz\nZgwAN998c1t5XV0dmzdv7nI9M7POVLMF1QL8S0TUA4cB50iqB2YDD0XEfsBD6XvSZScBk4AZwLWS\n+pc7yJkzYcECmDABpOTnggU968VXrKlTp1JfX8/HP/5xTj31VKZNm9brx/jGN77Bq6++Sn19PZdc\ncgn19fXsscceu6wzbNgwrrvuOo455hgaGxsZPXp027Lzzz+f8847j6lTpxIRbeWf+cxnWL58OVOm\nTOFnP/tZp+uZmXVGWTlZSLoHuDr996mI2CBpNPDfEbG/pDkAETE/Xf8B4OKIeKzQfhsbG6P9hIXP\nP/88BxxwQDmqsdtpaWmhpaWFIUOGsHr1aqZPn87q1asZMKDytyf9uZjVLklLIqKx6zV3ykQnCUkT\ngSnAE8BeEbEhXfQ6sFf6egzweM5m69OyfPubBcwCGF+um0U14s9//jOf/exnaWlpISL46U9/WpXk\nZGbWXtXPRJKGAz8H/ikimnPHZIuIkFRyEy8iFgALIGlB9VastWjkyJEsWbKk2mGYmXVQ1V58kgaS\nJKemiPhFWvxGemmP9OebafmrwLiczcemZWZmVoOq2YtPwA3A8xHxo5xF9wKnpa9PA+7JKT9J0mBJ\n+wD7AU9WKl4zM6usal7imwacAjwrqXVsnX8FLgMWSjoTWAv8A0BEPCdpIbCSpAfgORGxo/Jhm5lZ\nJVQtQUXEb4DOJgH6bCfbzAMq8ASSmZlVm0eSqILemG4D4MYbb+T111/vcr01a9bQ0NBQcJ0//OEP\n3HHHHUUf28ys3Kreiy/rRswfweZtmzuU1w2qo3lOc7f2Wcx0G8W48cYbmTp1KnvvvXe34sjVmqBO\nOumkHu/LzKw3uAXVhXzJqVB5T91yyy0ceuihNDQ0cPbZZ/P+++/T0tLCKaecwkEHHcSBBx7IlVde\nyZ133smyZcv44he/mLfl9dRTT3HwwQfT0NDAdddd11b+4osv8olPfIIpU6ZwyCGH8MQTTwAwe/Zs\nHn74YRoaGrjyyis7Xc/MrFLcgsqQFStWcNddd/G73/2OAQMGMGvWLO644w723Xdf3nrrLZ599lkg\nGRl85MiRXHXVVVx99dV5L9995StfYcGCBUybNo1vf/vbbeWjR49m0aJFDBkyhBdeeIHTTjuNJ554\ngssuu4yrr76au+++G0im2Mi3nplZpThBZciDDz7IU0891Tbdxrvvvsu4ceM4+uijWbVqFd/85jf5\n3Oc+x/Tp0wvu56233uLdd99tG7vvlFNO4eGHHwZg69atnHvuuSxfvpwBAwZ0Ov1FseuZmZWLE1SG\nRARnnHEG3//+9zsse+aZZ7j//vu55ppr+PnPf86CBQu6dYwf/vCHjBs3jttuu43t27czfPjwHq1n\nZlYuvgeVIUcddRQLFy7krbfeApLefuvWrWPjxo1EBCeeeCKXXnopS5cuBTpOadFqzz33ZOjQoTz2\nWDKOblPO/PRvv/02o0ePRhK33HJL28ji+abHyLeemVmlOEF1oW5QXUnlPXHQQQdx0UUXcdRRR3Hw\nwQczffp03njjDV555RWOPPJIGhoaOP300/nBD34AwOmnn85Xv/rVvJ0kbrrpJs466ywaGhro12/n\nx3zuuedy/fXXM3nyZF566aW2SQynTJnCjh07mDx5MldeeWWn65mZVUpmptsoF0+3sfvw52JWu7oz\n3YZbUGZmlklOUGZmlkl9NkHV+qXN3Y0/DzNrr08mqCFDhrBp0yafFDMiIti0aRNDhgypdihmliF9\n8jmosWPHsn79ejZu3FjtUCw1ZMgQxo4dW+0wzCxD+mSCGjhwIPvss0+1wzAzswL65CU+MzPLPico\nMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPL\nJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCcoMzPLJCco\nMzPLJCcoMzPLJCcoMzPLpKomKEk3SnpT0oqcsg9JWiRpdfrzgznL5khaI2mVpKOrE7WZmVVCtVtQ\nNwMz2pXNBh6KiP2Ah9L3SKoHTgImpdtcK6l/5UI1M7NKqmqCiohHgD+2Kz4OuCV9fQtwfE75HRGx\nNSJeAtYAh1YkUDMzq7hqt6Dy2SsiNqSvXwf2Sl+PAV7JWW99WtaBpFmSFktavHHjxvJFamZmZZPF\nBNUmIgKIbmy3ICIaI6Lxwx/+cBkiMzOzcstignpD0miA9OebafmrwLic9camZWZmVoOymKDuBU5L\nX58G3JNTfpKkwZL2AfYDnqxCfGZmVgEDqnlwSbcDnwL2lLQeuAi4DFgo6UxgLfAPABHxnKSFwEqg\nBTgnInZUJXAzMyu7qiaoiPhSJ4s+28n684B55YvIzMyyIouX+MzMzJygzMwsm5ygzMwsk5ygzMws\nk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5yg\nzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwsk5ygzMwskwZUOwAzM9s9jZg/gs3bNncorxtUR/Oc\n5h7v3y0oMzPrlnzJqVB5qdyCMjOzXqdL1Pa6blBdt/bhBGVmlnHlvJRW7st00P0WVVGX+CRNk/SB\n9PXJkn4kaUK3jmhmVuOammDiROjXL/nZ1NSz/RW6lFbOfVdbsS2onwCTJU0G/gW4HrgV+GS5AjMz\n2x01NcGsWbBlS/J+7drkPcDMmYW37aw1U0j7fZfSIuppciu3YhNUS0SEpOOAqyPiBklnljMwM+vb\nevPSU7En/t64rDV37s7k1GrLlqT86+sK16k7rZbWfbcmqFJaRHPnAqeXfMg2dYPqytrSKjZBbZY0\nBzgZOFJSP2Bg2aIysz6vlBNtV8ms2JNovvVyb/a3FxdFh7J16/Kvu24dRJkup3V2zHJt16p5TnPB\n/5+eKjZBfRH4MnBmRLwuaTxwedmiMjMrQW/eR5k4EV47bQTb+3W9bfuTc92gOsaPb2bt2o7rjh8P\neYp77mIRQL+5dcQPmuHi4jftKqbW+vW0ZVk3qI7NlP5ZFJWgIuJ14Ec579eR3IMyM+u2Qi2fQpqa\nur6f011rvzwCikhO+WzetplBf+5YPmwYzJsHJ6/pYXAFxKCuYx4xf8QuiWbePDj5uToYXHjb7rbw\ncluY+tfSW1oFE5SkzUDHNiwIiIgYUfIRzSyTmpqSexLr1iXfrOfNS5JAObshd7flU2yng27p4mTd\nlU2bdn0/ahRccUUS68mX9GjXPdb+/3XmTPjt2c1cdx1EABd3nkTKeSmvMwUTVER07+kqM8ucQonm\nJ+ObO+15VkoS6SzJFTp+d7TvGJBlw4cnPydOpGCHhKYmYGvXrZkulbiP++5Lk1MZ7PKZj+aQUrdX\nlBCZpL8AhrS+Ty/1ZVpjY2MsXry42mGYlU2hpJCr0Dfg/t8PduzoWD5hAqw9vchvztvqGPTvzWzb\nllM2Z0TPT7i7m4s7nlOHDUuTf6EWyra6oi7TFXXsAsdp37GjX7+cBFVgu1J16OH3U4jXoqQDFHUP\nStLngR8CHwHeBCYAzwOTSjmYmXUuX6IBOOsseOeddKXOTvhfroP5zaz98ghOXrO546WkbXUwqPNj\n77gw57wRJBfxKfGm/qDNbMu9z9AbrYFedttHY2dLsRdPxl1p3+08nx4nJ+hWncaPJ2+njp5oTYI9\nvSxYbC++7wOHAQ9GxBRJnybpcm5mvSDfw50nP5cmo/OK2MHgzYVbK6Wc/HrrvF3h5FTMSfHkNYLv\nViqiPDpJ2l22nnJbZb3cKp03L+d3L2NfKoodzXx7RGwC+knqFxEPA41ljMtst9LV0Dbtl599dvJT\nggED4OST83zLLvVEkaETSzX0xpBCZTe/OUk2F0eSDFJFtZ7mjEhaSO0/563FdxVo7R05Yv4IdInQ\nJeLkNWLLd7Wz9dUaX55LlcUaMb93+s8V24L6k6ThwCNAk6Q3gXe62MYsM7rqidZ6eW3tWujfH3bs\nSO6/HHssLFy4s2dWbo+sVvlaP2ecAd/6VrKdtOtN6LVr4Sc/2fk+370fK93atUmiL+U5oF5XQrIo\n6QtFoUt3Re5nwk3Rdtm4084q+ZJfN7749FZnmKI6SaQDxb5H0vifCewBNKWtqkxzJ4ndT6FLNBNu\nig73aHraQaBNzr2XopS6/tbkPlHRKniPpCYU0UGg4gp95pWO8+LoncuDW+uY8L+bu+w8c9tHY9eh\nnbrRSaKkXny7IyeobOiyi3H6S588zFjgdziAS/L8zl6k/Mmi1CRSbp3Fn0+WTrRZl7XPOVdW7utc\nHL32O6VLguHzuv6b1mXNSSeM01XWXny5D+wOIhmH751qPKgraQZwBdAfuD4iLqt0DNZ5j7PWy1q7\nKOZb2+DNOR0DCqwnSvsjy9pJKzf+9t+u+2KX7GIU0zLK2uecqxKfaYWT4LBhRYzDN3gzET3rIVjs\nUEdtF1YlCTiOpFdfRUnqD1wD/DWwHnhK0r0RsbLSsfRlnd1zaWmB99/Ps0Epfzh96QQ9eHNy0m1N\nVH2p7q1KuSzX2Um4pyfn1vtGu/P/f1exl3JvrAjvvFNih5StddCNsfi6fYlP0tMRMaVbG3eTpMOB\niyPi6PT9HICImN/ZNo1S+AKfmVl1iUOIWFyWS3xfyHnbj6SL+XulHKiXjAFeyXm/Hvgf7VeSNAuY\nBZQ+toaZmWVCsd3M/zbndQvwMsllvkyKiAXAAkg6SeBOEr1q4sQSrysXe8+oF2/iWsa1f8amO5/7\ntjr4QXP3t28Xy8CBsH1uEfvJ93xQZ/cPe3D58baPRtujD7s8qlBKXdvH2ku9+JjfzIQJ8PLLRfSQ\nbYuh9Edni70H1YM5F3vVq8C4nPdj0zKroF2ePE8NGlTgHpRVX2/3cuvtm/Ld2d+gzTvHuOthPK3P\nt3V7OoxCjw90M3nNnLnzkYncTkk96nc9v7n4+33zm9ueCdylDum907XAiPl1ZZ1Vt6vpNq6iwP9H\nRHyz1yMq7ClgP0n7kCSmk0gmUrQKav2jKboXX6V6GNXCze7eku9bfi+2TkctaO5eb83O5DvBFxHv\nggVw2mmwo9gTbx65t+F7Mh3GhAm7/j389rdJfDvSE/2sWTBtGpxySumtodxkpW7GOGpUMrL62kKd\nTdLPYdgwWHBbEivQ6ee6edtm4qIo21QcXbWgWq+NTQPqgTvT9ycCFe85FxEtks4FHiDpZn5jRDxX\n6Ths1z+Y9uUd7Tz5FJw+O5JvZKV+G8s37XarpiY4ZeWIooaSaZ12Yu7cEkbwbguCvC2U3DmTunoW\nrDe/iY4aBX/8464PLxd7Ymv9/yw0+sYVV3RsRTO/meHDkx5ebc++dELaGdvJPRjVs/X3bdYsKGI8\n1g7aT4zY3c9g1Kjkclf72K69tuO6PalvSXJ67g0bljsCSjN77pnni2RqwoSdvzOtlxirpav5oG4B\nkPR14IiIaEnfXwc8Wv7w8sZ0H3BfNY5tPdfZCaD1RJE7AV4x8wd1NfPqzJkwk9Im1SvmZN6dyfqa\n5zQXNflfT7+N1g2q4623erQLgKLqV2gUj86+bNQNqqP5/V330ZOTYOsxu7o8V+iLTKtiPoOBA2H7\n9l3fX3FFl7tuM2FCTn07ac0U+r3uNInm3Bs69thknqd1yv/Z5PuCMWxY0uLLXa/tcn4XdSoYUw8U\nO9TRKuDwiPhj+v6DwOMRsX+Pjl4BHknCuqOcs8h2pZQE1Tq+WjET9xW732JO5L2p/XN1kJws3z9/\nBO9F4S8ou0wpXqB+pX5uXU3uWMzwWp3prL7tk0O5FTuPWFNT4dFdurqC0XbvLCZvi1he6DH8DopN\nUKeTDMH4MMmFjCNJnke6pZSDVYMTlO1uSklQpSSTYme0rXSCgu7NxNs+6VTzS0Wpik0OWVHw0nyR\nvy+SlkRESV35in5QV9Le7Hzm6ImIeL2UA1WLE5TtbopNUD058e5OJ3Orvt74felOguqqF9/HI+IF\nSVPTotaHZD8i6SMRsbSUg5lZz/RW68ZJyEpRrd+Xrnrx/TPJiAw/zLMsgM/0ekRmfVxXHUnM+oqu\nevHNSn9+ujLhmJlbN2aJoqZ8l3SipLr09QWSfiGpogPFmplZ31JUggIujIjNko4AjgJuAK4rX1hm\nZtbXFTtY7I705+eABRHxS0n/VqaYzMxqkntPlqbYFtSrkn4KfBG4T9LgErY1MzPo9Dm0cg22ursr\nNsn8A8n4d0dHxJ+ADwHnlS0qMzPr84pKUBGxBXgTOCItagFWlysoMzOzYnvxXQScD8xJiwYCt5Ur\nKDMzs2Iv8f0d8HngHYCIeA3wU4NmZlY2xSaobZEM2hcAkj5QvpDMzGpTZ6OBeJSQ/IrtZr4w7cU3\nUtLXgDOA68sXlplZ7XFX8tIUlaAi4t8l/TXJ1Kj7A9+LiEVljczMzPq0YltQpAlpEYCkfpJmRkRT\n2SIzM7M+reA9KEkjJM2RdLWk6UqcC/yB5NkoMzOzsuiqBfW/gP8HPAZ8FfhXkhl1j4+IZWWOzczM\n+rCuEtRfRsRBAJKuBzYA4yPivbJHZmZmfVpX3cy3t76IiB3AeicnMzOrhK5aUJMltfaLFDA0fS8g\nImJEWaMzM7M+q6sZdftXKhAzM7NcnjLDzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwy\nyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyyQnKzMwyqSoJStKJkp6T\n9L6kxnbL5khaI2mVpKNzyg+R9Gy67EpJqnzkZmZWKdVqQa0AvgA8klsoqR44CZgEzACuldQ6J9VP\ngK8B+6X/ZlQsWjMzq7iqJKiIeD4iVuVZdBxwR0RsjYiXgDXAoZJGAyMi4vGICOBW4PgKhmxmZhWW\ntXtQY4BXct6vT8vGpK/bl+claZakxZIWb9y4sSyBmplZeRWc8r0nJD0I7J1n0dyIuKdcxwWIiAXA\nAoDGxsYo57HMzKw8ypagIuKobmz2KjAu5/3YtOzV9HX7cjMzq1FZu8R3L3CSpMGS9iHpDPFkRGwA\nmiUdlvbeOxUoayvMzMyqq1rdzP9O0nrgcOCXkh4AiIjngIXASuBXwDkRsSPd7GzgepKOEy8C91c8\ncDMzqxglneJqV2NjYyxevLjaYZiZ9WmSlkREY9dr7pS1S3xmZmaAE5SZmWWUE5SZmWWSE5SZmWWS\nE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZ\nmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWS\nE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZ\nmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWWSE5SZmWVSVRKUpMslvSDpGUl3SRqZs2yOpDWSVkk6\nOqf8EEnPpsuulKRqxG5mZpVRrRbUIuDAiDgY+D0wB0BSPXASMAmYAVwrqX+6zU+ArwH7pf9mVDpo\nMzOrnKokqIj4PxHRkr59HBibvj4OuCMitkbES8Aa4FBJo4EREfF4RARwK3B8xQM3M7OKycI9qDOA\n+9PXY4BXcpatT8vGpK/bl+claZakxZIWb9y4sZfDNTOzShhQrh1LehDYO8+iuRFxT7rOXKAFaOrN\nY0fEAmABQGNjY/Tmvs3MrDLKlqAi4qhCyyV9Bfgb4LPpZTuAV4FxOauNTcteZedlwNxyMzOrUdXq\nxTcD+C7w+YjYkrPoXuAkSYMl7UPSGeLJiNgANEs6LO29dypwT8UDNzOziilbC6oLVwODgUVpb/HH\nI+IfI+I5SQuBlSSX/s6JiB3pNmcDNwNDSe5Z3d9hr2ZmVjOqkqAi4qMFls0D5uUpXwwcWM64zMws\nO7LQi8/MzKwDJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8sk\nJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygzM8skJygz\nM8skRUS1YygrSZuBVdWOo0L2BN6qdhAV4rrWrr5U375U1/0joq6UDQaUK5IMWRURjdUOohIkLXZd\na09fqiv0rfr2tbqWuo0v8ZmZWSY5QZmZWSb1hQS1oNoBVJDrWpv6Ul2hb9XXdS2g5jtJmJnZ7qkv\ntKDMzGw35ARlZmaZVJMJStKJkp6T9L6kxnbL5khaI2mVpKOrFWNvkjQjrc8aSbOrHU9vk3SjpDcl\nrcgp+5CkRZJWpz8/WM0Ye4ukcZIelrQy/R3+Vlpec/WVNETSk5KWp3W9JC2vubq2ktRf0tOS/it9\nX8t1fVnSs5KWtXYxL7W+NZmggBXAF4BHcgsl1QMnAZOAGcC1kvpXPrzek8Z/DXAMUA98Ka1nLbmZ\n5PPKNRt4KCL2Ax5K39eCFuBfIqIeOAw4J/08a7G+W4HPRMRkoAGYIekwarOurb4FPJ/zvpbrCvDp\niGjIedarpPrWZIKKiOcjIt/oEccBd0TE1oh4CVgDHFrZ6HrdocCaiPhDRGwD7iCpZ82IiEeAP7Yr\nPg64JX19C3B8RYMqk4jYEBFL09ebSU5mY6jB+kbiz+nbgem/oAbrCiBpLPA54Pqc4pqsawEl1bcm\nE1QBY4BXct6vT8t2Z7VYp2LsFREb0tevA3tVM5hykDQRmAI8QY3WN73ktQx4E1gUETVbV+DHwHeB\n93PKarWukHzZeFDSEkmz0rKS6rvbDnUk6UFg7zyL5kbEPZWOx6onIkJSTT0vIWk48HPgnyKiWVLb\nslqqb0TsABokjQTuknRgu+U1UVdJfwO8GRFLJH0q3zq1UtccR0TEq5L+Algk6YXchcXUd7dNUBFx\nVDc2exUYl/N+bFq2O6vFOhXjDUmjI2KDpNEk38BrgqSBJMmpKSJ+kRbXbH0BIuJPkh4muddYi3Wd\nBnxe0rHAEGCEpNuozboCEBGvpj/flHQXye2Ikurb1y7x3QucJGmwpH2A/YAnqxxTTz0F7CdpH0mD\nSDqB3FvlmCrhXuC09PVpQE20mpU0lW4Ano+IH+Usqrn6Svpw2nJC0lDgr4EXqMG6RsSciBgbERNJ\n/kZ/HREnU4N1BZD0AUl1ra+B6SSd10qrb0TU3D/g70juxWwF3gAeyFk2F3iRZAqOY6oday/V91jg\n92m95lY3wM+FAAAD+UlEQVQ7njLU73ZgA7A9/VzPBEaR9AJaDTwIfKjacfZSXY8guXb/DLAs/Xds\nLdYXOBh4Oq3rCuB7aXnN1bVdvT8F/Fct1xX4S2B5+u+51vNSqfX1UEdmZpZJfe0Sn5mZ7SacoMzM\nLJOcoMzMLJOcoMzMLJOcoMzMLJOcoKzPk7QjHXF5haT/lDSsB/v6VM5I1Z8vNLq8pJGSzu7GMS6W\n9J3uxtjb+zErFycoM3g3khGXDwS2Af+Yu1CJkv9WIuLeiLiswCojgZITlFlf4QRltqtHgY9KmpjO\nsXUryUOk4yRNl/SYpKVpS2s4tM3H9YKkpSTTvJCWf0XS1enrvSTdlc59tFzS/wQuA/ZNW2+Xp+ud\nJ+kpSc+0zo+Uls+V9HtJvwH2bx+0pD0krW1NpOmT/K9IGijpa+k+l0v6eb4WoqT/Vjp3mqQ9Jb2c\nvu4v6fKcmM5Ky0dLeiSn5fmJ3vjPN8vlBGWWkjSAZF6tZ9Oi/YBrI2IS8A5wAXBUREwFFgP/LGkI\n8B/A3wKHkH8AY4Argf8bydxHU0merp8NvJi23s6TND095qEk8yMdIulISYeQDI/TQDKqxF+133lE\nvE0y6sQn06K/IRlBZTvwi4j4q/TYz5OMxFGsM4G3I+Kv0uN+LR0m7Mvp/huAyemxzXrVbjtYrFkv\nGppO+QBJC+oG4CPA2oh4PC0/jGRCyN+mI4sPAh4DPg68FBGrAdIBQGfR0WeAU6FtBO+31XE20enp\nv6fT98NJElYdcFdEbEmP0dlYi3cCXwQeJklo16blB0r6N5JLisOBBwr9Z+SJ6WBJJ6Tv90hjegq4\nMR3Y9u6IcIKyXucEZZbeg8otSJPQO7lFJPMVfanderts10MC5kfET9sd45+K3P5e4AeSPkTSmvt1\nWn4zcHxELJf0FZKx4NprYecVlSHtYvpGRHRIapKOJJmA72ZJP4qIW4uM06wovsRnVpzHgWmSPgpt\n93g+RjL69kRJ+6brfamT7R8Cvp5u21/SHsBmktZRqweAM3LubY1RMpfOI8DxkoamI0T/bb4DRDI7\n7VPAFSSDke5IF9UBG9LWzsxO4nuZJKkBnJBT/gDw9XRbJH0srfsE4I2I+A+SGWKndrJfs25zC8qs\nCBGxMW193C5pcFp8QUT8Xslsob+UtIXkEmFdnl18C1gg6UxgB/D1iHhM0m8lrQDuT+9DHQA8lrbg\n/gycHBFLJd1JMjL0myRJqDN3Av/Jrq2kC0lm5d2Y/swX378DC1vrklN+PTARWKokqI0k03R/CjhP\n0vY0zlMLxGTWLR7N3MzMMsmX+MzMLJOcoMzMLJOcoMzMLJOcoMzMLJOcoMzMLJOcoMzMLJOcoMzM\nLJP+P/wZI7CB6SGNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff476ccd080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "forest = RandomForestRegressor(\n",
    "    n_estimators=1000, criterion='mse', random_state=1, n_jobs=-1)\n",
    "forest.fit(X_train, y_train)\n",
    "y_train_pred = forest.predict(X_train)\n",
    "y_test_pred = forest.predict(X_test)\n",
    "\n",
    "print('MSE train: %.2f, test: %.2f' %\n",
    "      (mean_squared_error(y_train, y_train_pred),\n",
    "       mean_squared_error(y_test, y_test_pred)))\n",
    "print('R^2 train: %.2f, test: %.2f' % (r2_score(y_train, y_train_pred),\n",
    "                                       r2_score(y_test, y_test_pred)))\n",
    "\n",
    "# Residual plot\n",
    "plt.scatter(\n",
    "    y_train_pred,\n",
    "    y_train_pred - y_train,\n",
    "    c='blue',\n",
    "    marker='o',\n",
    "    label='Training data')\n",
    "plt.scatter(\n",
    "    y_test_pred, y_test_pred - y_test, c='green', marker='s', label='Test data')\n",
    "\n",
    "plt.xlabel('Predicted values')\n",
    "plt.ylabel('Residuals')\n",
    "plt.legend(loc='upper left')\n",
    "plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='red')\n",
    "plt.xlim([-10, 50])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}