Analysis.py.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.ensemble import AdaBoostRegressor, GradientBoostingRegressor, ExtraTreesRegressor, RandomForestRegressor\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn import model_selection\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "from sklearn import preprocessing\n",
    "from sklearn.cross_validation import cross_val_score\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import r2_score, mean_squared_error\n",
    "from sklearn.linear_model import ElasticNet\n",
    "from sklearn.linear_model import Ridge \n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.linear_model import RidgeCV\n",
    "import xgboost as xgb\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import KFold, GridSearchCV\n",
    "%matplotlib inline\n",
    "\n",
    "#Stop warnings\n",
    "import warnings\n",
    "def ignore_warn(*args, **kwargs):\n",
    "    pass\n",
    "warnings.warn = ignore_warn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv('ATM_Training.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>Shops</th>\n",
       "      <th>ATMs</th>\n",
       "      <th>Downtown</th>\n",
       "      <th>Weekday</th>\n",
       "      <th>Center</th>\n",
       "      <th>High</th>\n",
       "      <th>Withdraw</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10.18</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>72.750556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>9.74</td>\n",
       "      <td>10</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>66.720482</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.96</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>19.189516</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>9.58</td>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>67.388669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.03</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>15.813127</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Shops  ATMs  Downtown  Weekday  Center  High   Withdraw\n",
       "0  10.18    10         1        0       0     0  72.750556\n",
       "1   9.74    10         1        1       0     0  66.720482\n",
       "2   0.96     2         0        0       0     1  19.189516\n",
       "3   9.58     9         1        1       0     1  67.388669\n",
       "4   1.03     4         0        1       0     1  15.813127"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EDA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#I'll leave this to you if that's okay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The variable type of  Shops  is float64\n",
      "The variable type of  ATMs  is int64\n",
      "The variable type of  Downtown  is int64\n",
      "The variable type of  Weekday  is int64\n",
      "The variable type of  Center  is int64\n",
      "The variable type of  High  is int64\n",
      "The variable type of  Withdraw  is float64\n"
     ]
    }
   ],
   "source": [
    "for x in list(data):\n",
    "    print('The variable type of ', x, ' is',data[x].dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Downtown, Weekday, Center, and high are dummy variables.\n",
    "\n",
    "Shops, ATMs, and withdraw are continuous variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Shops      -0.877958\n",
       "ATMs       -0.318768\n",
       "Downtown   -0.883353\n",
       "Weekday    -0.947692\n",
       "Center      2.622120\n",
       "High        0.864683\n",
       "Withdraw   -0.772358\n",
       "dtype: float64"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Only worthwhile to look at shops, ATMs, and withdraw\n",
    "data.skew()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Shops        7.316373\n",
       "ATMs         7.937455\n",
       "Downtown     0.702000\n",
       "Weekday      0.714091\n",
       "Center       0.102455\n",
       "High         0.301591\n",
       "Withdraw    54.652818\n",
       "dtype: float64"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Shops        9.890000\n",
       "ATMs         9.000000\n",
       "Downtown     1.000000\n",
       "Weekday      1.000000\n",
       "Center       0.000000\n",
       "High         0.000000\n",
       "Withdraw    68.240749\n",
       "dtype: float64"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Construction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#By test, I refer to the validation set in which we will test our models on.\n",
    "\n",
    "final_train = data.sample(frac=0.6, random_state=450411920)\n",
    "final_test = data[data.index.isin(final_train.index)==False]\n",
    "#Now we have final train/test which only has predictors whilst y_train/test are the response\n",
    "y_train = final_train.pop('Withdraw')\n",
    "y_test = final_test.pop('Withdraw')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see more information about the models I use in this: http://scikit-learn.org/stable/modules/ensemble.html#gradient-boosting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hyperparameter Tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we try to tune the hyperparameter in our models, we use this package called randomizedsearchcv.\n",
    "\n",
    "This allows us to try out different hyperparameters.\n",
    "\n",
    "http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html#sklearn.model_selection.RandomizedSearchCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#In order to store the results from our models\n",
    "modelname = []\n",
    "models = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Standardising the data\n",
    "scaler = StandardScaler()\n",
    "scaler.fit(final_train)\n",
    "train_scale = scaler.transform(final_train)\n",
    "scaler.fit(final_test)\n",
    "test_scale = scaler.transform(final_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Try setting from 1 - 25 neighbors\n",
    "k_range = range(1, 25)\n",
    "\n",
    "# we create a list. This allows us to see whether we should weigh all neighbours equally or weigh closer ones more\n",
    "weight_options = ['uniform', 'distance']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Creating the objects necessary for us to try cross-validation in order to locate best hyperparameters\n",
    "param_dist = dict(n_neighbors=k_range, weights=weight_options)\n",
    "knn = KNeighborsRegressor()\n",
    "rand = RandomizedSearchCV(knn, param_dist, cv=10, scoring='neg_mean_squared_error', n_iter=10, random_state=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Fitting the model and testing out random hyperparameters. This saves on computation time\n",
    "knnopt = rand.fit(final_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of this model is 0.479\n",
      "Best parameters from this is  {'weights': 'distance', 'n_neighbors': 9}\n",
      "This is the best KNN estimator KNeighborsRegressor(algorithm='auto', leaf_size=30, metric='minkowski',\n",
      "          metric_params=None, n_jobs=1, n_neighbors=9, p=2,\n",
      "          weights='distance')\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print('MSE of this model is %.3f'%np.negative(rand.best_score_))\n",
    "print('Best parameters from this is ',rand.best_params_)\n",
    "print('This is the best KNN estimator', rand.best_estimator_)\n",
    "models.append(knnopt)\n",
    "modelname.append('KNN')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.62975447467504286, -0.49265331541510449, -0.47201821877027417, -0.48343268344064799, -0.51228980857997686, -0.5482003247759184, -0.58442889195058823, -0.64126882365103521, -0.69626566178245475, -0.75309255703045574, -0.8439416629723302, -0.92346914690549986, -0.99552476713283689, -1.0576872512738387, -1.1253266265406536, -1.2075369939491591, -1.2937798170150558, -1.3819374370122393, -1.4721009238087881, -1.5475567132612265, -1.6308032455387997, -1.694247371502867, -1.7801438833549124, -1.859940424471676]\n"
     ]
    }
   ],
   "source": [
    "#Now constructing a graph for this\n",
    "k_range = range(1, 25)\n",
    "\n",
    "# list of scores from k_range\n",
    "k_scores = []\n",
    "\n",
    "# 1. we will loop through reasonable values of k\n",
    "for k in k_range:\n",
    "    # 2. run KNeighborsRegressor with k neighbours\n",
    "    knn = KNeighborsRegressor(n_neighbors=k, weights = 'uniform')\n",
    "    # 3. obtain cross_val_score for KNeighborsRegressor with k neighbours\n",
    "    scores = cross_val_score(knn, final_train, y_train, cv=10, scoring='neg_mean_squared_error')\n",
    "    # 4. append mean of scores for k neighbors to k_scores list\n",
    "    k_scores.append(scores.mean())\n",
    "print(k_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11455eb38>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecFPX9x/HXh96LUqSDVOnlqKJiwx7sgAVBBDXGaIxG\nkx+xxRRLjEaTEEREUaoNuyL2RMrRe5HeezvKcbef3x87xL3zbm+B29sr7+fjcY/bnfnOzOeGZT77\nLfMdc3dERESOKZboAEREJH9RYhARkQyUGEREJAMlBhERyUCJQUREMlBiEBGRDJQYpMAws+Fm9vsY\ny442syeirHcza5J70YkUHkoMkmfMbI2ZbTOz8hHLbjOzr2LZ3t3vcPc/xC3AQiI4z6lmVi3T8jlB\nQmwYvK9rZm+Z2Q4z22tmC81sYLCuYVD2QKafvnn+B0meU2KQvFYcuCfRQeQVMyuRoEOvBvpHxNEG\nKJepzBhgPdAAOBW4GdiaqUwVd68Q8TMhjjFLPqHEIHntaeB+M6uS1Uoza2FmU8xsl5ktM7PrI9Zl\naB4ys9+Y2WYz2xTUPDI3D1U1sw/NbL+ZTTezxpkOd6mZrQq+MT9tZsWC/RYzs2Fmtjao4bxmZpWD\ndb3MbEOmmNeY2QXB60fN7E0ze93M9gEDzayLmSWb2T4z22pmz2bzty8xs8sj3pcws+1m1tHMygT7\n3Glme8xsppnVjHKexwADIt7fAryWqUxnYLS7p7h7mrvPcfePo+xTigglBslrycBXwP2ZVwRNTFOA\nsUANoB/wTzNrmUXZi4H7gAuAJkCvLI7VD3gMqAqsBP6Yaf1VQBLQEegD3BosHxj8nAucDlQAXozx\n7yPY15tAFeAN4HngeXevBDQGJmaz3TgivuUDFwE73H024Qt7ZaAe4W/3dwCHosQwDahkZmeYWXHC\n5+L1LMr8w8z6mVn94/j7pJBTYpBEeBi428yqZ1p+ObDG3V859g0WeAu4Lot9XA+84u6L3P0g8GgW\nZd5x9xnunkb4At0+0/on3X2Xu68DnuPHi/KNwLPuvsrdDwC/BfodR7PQ9+7+rruH3P0QcBRoYmbV\n3P2Au0/LZruxwM/M7FiTzw2EkwXBPk4Fmrh7urvPcvd9OcRxrNZwIbAE2Jhp/XXAt8DvgdVmNtfM\nOmcqsyOooRz7OSOHY0ohoMQgec7dFwIfAA9lWtUA6Bp5ISJ8kT4ti93UJtw+fsz6LMpsiXh9kPA3\n/0iR26wN9nls32szrSsBRGu6yW6/AIOBZsDSoAno8iy2wd1XEr6AXxEkh58RThYQvsh/CowPms6e\nMrOSOcQxhnByGchPm5Fw993u/pC7tyL8t80F3jUziyhWzd2rRPwsyeGYUggoMUiiPAIMAepELFsP\nfJ3pQlTB3e/MYvvNQN2I9/VOIIbIbeoDm4LXmwgnqch1aYQ7ZlOI6MQNmmky13wyTFns7ivcvT/h\n5rEngTcjR2Zlcqw5qQ+wOEgWuPtRd3/M3VsCPQjXrgZks49jx11LuBP6UuDtHMruAJ4hnBRPiVZW\nCj8lBkmI4II3AfhlxOIPgGZmdrOZlQx+OmfTfDERGBS0oZcj3BxyvB4ws6pmVo/wSKljI27GAb8y\ns0ZmVgH4EzAhaJJaDpQxs8uCb+zDgNLRDmJmN5lZdXcPAXuCxaFsio8HegN38mNtATM718zaBIlo\nH+Gmpez2EWkwcJ67p2QR15Nm1jro5K4YHHOlu++MYb9SiCkxSCI9Dvzvm7O77yd8UexH+Fv7FsLf\nsH9y4Q1Gz/wd+JJwx/Kxdvsjx3H8ycAswk0oHwIvB8tHEW6G+YbwN+7DwN3BcfcCPwdGEm6zTwEy\njFLKwsXAIjM7QLgjul/Q9/AT7r4Z+J5wrSByaOhphDu09xFubvo6iDEqd//B3ZOzWV0OeIdwslpF\nuJb0s0xl9mS6j+G+nI4pBZ/pQT1SGAS1ioVA6eCbvYicINUYpMAys6vMrLSZVSVcs3hfSUHk5Ckx\nSEF2O7AN+AFIJ9xGLiInSU1JIiKSgWoMIiKSQaIm+Dph1apV84YNGyY6DBGRAmXWrFk73D3zPTdZ\nKnCJoWHDhiQnZzf6TkREsmJma3MuFaamJBERyUCJQUREMlBiEBGRDJQYREQkAyUGERHJQIlBREQy\nUGIQEZEMlBhERPI5d+eFqStYvCmnp7nmjgJ3g5uISFFyJC2d37w5n8lzN5GSmk7L2pXifkwlBhGR\nfGpXSiq3j0lm5prdPHBRc37eq3GeHFeJQUQkH1q9I4VBr8xg097DvNC/A1e0q51nx1ZiEBHJZ2as\n3sXQMckUM2PckK50anBKnh5fiUFEJB95Z84GHnxzAXVPKcsrAzvT4NTyOW+Uy5QYRETyAXfn+akr\neO7zFXQ7/RT+fVMSlcuVTEgscRuuamajzGybmS3MZn1lM3vfzOaZ2SIzGxSvWERE8rMjaencN3Ee\nz32+gms71eW1W7smLClAfO9jGA1cHGX9XcBid28H9AL+amal4hiPiEi+s+dgKje/PIN35mzk/t7N\nePratpQqkdhbzOLWlOTu35hZw2hFgIpmZkAFYBeQFq94RETymzU7Urh19Ew27D7E8/3a06d9nUSH\nBCS2j+FF4D1gE1AR6OvuoawKmtlQYChA/fr18yxAEZF4SV6ziyGvhZ9G+caQrnRumLcjj6JJZH3l\nImAuUBtoD7xoZlne0ufuI9w9yd2TqleP6ZGlIiL51uS5G7nhpelUKVeKd35+Zr5KCpDYxDAIeNvD\nVgKrgRYJjEdEJO5GfruKe8bPpX39Krx9Zw8aVsv74ag5SWRiWAecD2BmNYHmwKoExiMiEjfuzlOf\nLOWJD5dwWZtajBncharl8+d4m7j1MZjZOMKjjaqZ2QbgEaAkgLsPB/4AjDazBYABD7r7jnjFIyKS\nKOkhZ9i7Cxk3Yx03dK3PH/q0pngxS3RY2YrnqKT+OazfBPSO1/FFRPKDI2np3DdhHh8u2Mxd5zbm\n/t7NCQ/GzL9057OISJykHEnjjtdn8e2KHQy77AxuO+v0RIcUEyUGEZE42J2SyqDRM5m/YQ9PXduW\n65PqJTqkmCkxiIjksi17D3Pzy9NZu+sg/7qpExe1Oi3RIR0XJQYRkVy0ekcKN42czt5DRxk9qDM9\nGldLdEjHTYlBRCSXLNy4l4GvzCDkMG5IN9rUrZzokE6IEoOISC6YsXoXg0fPpGKZEoy5rSuNq1dI\ndEgnTIlBROQkTV2ylZ+/MZu6VcsyZnBXalcpm+iQTooSg4jISXhnzgbunzSfVrUrMXpQF07Jp3cz\nHw8lBhGRE5CaFuLvU1fw4pcr6X76qbx0SxIVSheOS2rh+CtERPLQ0i37uG/CPBZv3sd1neryhytb\nU6Zk8USHlWuUGEREYpQeckZ8s4q/TVlOpbIleGlAEhe2rJnosHKdEoOISAzW7Ejh15PmMWvtbi5p\nfRpPXNmaUyuUTnRYcaHEICISRSjkvD59LX/+aCklixvP9W1Pn/a18/1EeCdDiUFEJBub9hziN2/O\n57uVOzi7WXWeuqYtp1Uuk+iw4k6JQUQkE3fn7dkbefT9RaSHnD9e1ZobutQv1LWESEoMIiIRdhw4\nwu/eXsBni7fSuWFVnrmuHQ1OzX+P34wnJQYRkcAnCzfzu3cWcuBIGv936Rnc2rNRvn7SWrwoMYhI\nkZeWHuLh9xYxdvo6WtepxLPXt6dZzYqJDithlBhEpEg7lJrOL8bOZurSbdx+zunc37s5JYsXS3RY\nCRW3v97MRpnZNjNbGKVMLzOba2aLzOzreMUiIpKVXSmp3DByGl8u28YTV7bmt5ecUeSTAsS3xjAa\neBF4LauVZlYF+CdwsbuvM7MacYxFRCSD9bsOcsuoGWzcc6hAPmUtnuKWGNz9GzNrGKXIDcDb7r4u\nKL8tXrGIiERatGkvA1+ZSWpaiDdu60pSw1MSHVK+ksg6UzOgqpl9ZWazzGxAdgXNbKiZJZtZ8vbt\n2/MwRBEpbP6zcgd9/z2NksWMN+/orqSQhUR2PpcAOgHnA2WB781smrsvz1zQ3UcAIwCSkpI8T6MU\nkULjvXmb+PXEuZxerQKjb+1MrcoF+4E68ZLIxLAB2OnuKUCKmX0DtAN+khhERE7WyG9X8cSHS+jS\n6BReGpBE5bIlEx1SvpXIpqTJQE8zK2Fm5YCuwJIExiMihVAo5Pzxw8U88eESLm1zGq/d2kVJIQdx\nqzGY2TigF1DNzDYAjwAlAdx9uLsvMbNPgPlACBjp7tkObRUROV6paSEeeHMek+du4pbuDXj4ilZF\n8k7m4xXPUUn9YyjzNPB0vGIQkaJr/+Gj3Pn6bL5buYPfXNycO89pXGQmwTtZuvNZRAqdbfsPM3DU\nTJZv3c9fr2vHNZ3qJjqkAkWJQUQKlcWb9jHktWR2H0xl5C1J9Gque2ePlxKDiBQaH87fzP2T5lG5\nbEnGD+1G27pVEh1SgaTEICIFXijkPDtlOS9+uZJODaryr5s6UqNi4X/SWrwoMYhIgbb/8FF+NWEu\nny/ZRt+kejx+ZStKlyie6LAKNCUGESmwVu9IYchryazekcLjfVpxc7cGGnmUC7K9wc3Mzot43SjT\nuqvjGZSISE6+WraNPi9+x66UVF4f3JUB3RsqKeSSaHc+PxPx+q1M64bFIRYRkRy5O//++gduHT2T\nOlXLMfmuM+ne+NREh1WoRGtKsmxeZ/VeRCTuDh9N58G35jN57iYua1OLp69rS7lSahHPbdHOqGfz\nOqv3IiJxtWnPIW4fM4uFm/bywEXN+Xkv3ckcL9ESw+lm9h7h2sGx1wTvG2W/mYhI7pq5Zhd3vj6L\nw0dDvHRzEhe0rJnokAq1aImhT8TrZzKty/xeRCQuxs1Yx8OTF1K3ajnGD+1EkxoVEx1SoZdtYnD3\nryPfm1lJoDWwUY/hFJG88PepK3h2ynLObladF/p1oHI5TZedF6INVx1uZq2C15WBecBrwBwzy3Hm\nVBGRE+XuPPPpMp6dspyrO9Rh1C1JSgp5KNpw1bPcfVHwehCw3N3bEH4c52/iHpmIFEnuzp8+WsKL\nX66kX+d6PHNdO0oUT+QzxYqeaH0MqRGvLwQmAbj7Fo0EEJF4CIWcx95fxKvfr+WW7g145IpWFNOD\ndfJctMSwx8wuBzYCZwKDAcysBKAnaItIrgqFnP97dwHjZqxnyFmN+N2lZ2g4aoJESwy3A38HTgPu\ndfctwfLzgQ/jHZiIFB3pIeeBN+fx9uyN3HVuY+7v3VxJIYGijUpaDlycxfJPgU/jGZSIFB1H00Pc\nN3Ee78/bxH0XNuOX5zdNdEhFXraJwcz+Hm1Dd/9l7ocjIkVJalqIu8fN5tNFW3nokhbccU7jRIck\nRB+VdAfQE9gEJAOzMv1EZWajzGybmS3MoVxnM0szs2tjD1tECrrDR9O54/VZfLpoK49c0VJJIR+J\n1sdQC7gO6AukAROAN919T4z7Hg28SPjehyyZWXHgSeCzGPcpIoXAodR0ho5J5tsVO/jjVa25sWuD\nRIckEbKtMbj7Tncf7u7nEr6PoQqw2MxujmXH7v4NsCuHYncTntJbd1KLFBEpR9IYNHoG363cwVPX\ntlVSyIdynK/WzDoC/Qnfy/AxMTQjxcLM6gBXAecCnXMoOxQYClC/fv3cOLyIJMD+w0cZ9MpM5qzf\nw3N929OnfZ1EhyRZiNb5/DhwGbAEGA/81t3TcvHYzwEPunsop2Fp7j4CGAGQlJSkKb9FCqC9B48y\n4JUZLNq4lxf6d+DSNrUSHZJkI1qNYRiwGmgX/PwpuIAb4O7e9iSPnQSMD/ZZDbjUzNLc/d2T3K+I\n5DM7Dhzh5pdn8MO2A/zrpk5cqGmz87VoiSGuz1xw9//t38xGAx8oKYgUPpv3HuLGkdPZtOcQI29J\n4uxm1RMdkuQg2g1ua09mx2Y2DugFVDOzDcAjQMlg38NPZt8iUjCs3ZnCjSOns/fgUcYM7krnhqck\nOiSJQdweluruMU/N7e4D4xWHiCTGiq37uXHkdI6mhxg7pBtt6lZOdEgSIz1FW0Ry3cKNe7n55emU\nKF6MCbd3p1lNPXWtIFFiEJFclbxmF4NemUmlsiV547auNKxWPtEhyXGKNlx1AZDt0NBcGJUkIoXM\ndyt2MOS1ZGpVLsPrt3WldhXN0F8QRasxXB78viv4PSb4fWP8whGRguqzRVv4xdg5nF69PGMGd6V6\nxdKJDklOUI6jkszsQnfvELHqITObDTwU7+BEpGCYPHcj902cR+s6lXl1UGeqlCuV6JDkJMTyIFUz\nszMj3vSIcTsRKQLGzVjHvRPmktSgKm/c1lVJoRCIpfN5MDDKzI6NNdsD3Bq/kESkoBj57Sqe+HAJ\nvZpX5183dqJsqeKJDklyQY6Jwd1nAe2OJQZ33xv3qEQkX3N3/j51JX/7fDmXtD6N5/t1oFQJNSQU\nFjn+S5pZTTN7GRjv7nvNrKWZDc6D2EQkH3J3/vTREv72+XKu6ViXF/orKRQ2sfxrjib8jOfawfvl\nwL3xCkhE8q+09BAPvjWfl75dzYDuDXj62raUKK6kUNjE8i9azd0nAiGAYOrt9LhGJSL5zpG0dH4x\ndg4Tkzfwy/Ob8tjPWlGsWPQp86VgiqXzOcXMTiW42c3MugHqZxApQg4cSeP2Mcn8Z+VOHr68Jbf2\njOvky5JgsSSG+4D3gMZm9h+gOuFnQYtIEbA7JZWBo2eycONe/npdO67pVDfRIUmcxZIYFgHnAM0J\nP6RnGbqPQaRI2LL3MDe/PJ21uw4yXA/YKTJiSQzfu3tHwgkCgODO545xi0pEEm71jhRuGjmdvYeO\n8uqgLnRvfGqiQ5I8Em0SvdOAOkBZM+tAuLYAUAkolwexiUiCLN60jwGjZhByZ5yepVDkRKsxXAQM\nBOoCz0Ys3w/8Lo4xiUgCJa/ZxaDRM6lQugRjBnejSY0KiQ5J8li0SfReBV41s2vc/a08jElEEuTL\nZdu48/VZ1K5cljG3daWOps0ukmKZEuMtM7sMaAWUiVj+eDwDE5G8NXnuRn49cR4talVk9KAuVKug\nabOLqlimxBgO9AXuJtzPcB3QIIbtRpnZNjNbmM36G81svpktMLP/mlm744xdRHLJmGlruXfCXDo2\nqMq4Id2UFIq4WIad9nD3AcBud38M6A40i2G70cDFUdavBs5x9zbAH4ARMexTRHKRu/PC1BX8/t2F\nnNe8Bq/d2oWKZUomOixJsFiGqx4Kfh80s9rATqBWThu5+zdm1jDK+v9GvJ1GuJNbRPJIWnqIR95b\nxBvT13FVhzo8dW1bSmreIyG2xPCBmVUBngZmE54aY2QuxzEY+Di7lWY2FBgKUL9+/Vw+tEjRcyg1\nnbvHzebzJdu4s1djHujdXPMeyf+Yu8de2Kw0UCbWZzIENYYP3L11lDLnAv8Eerr7zpz2mZSU5MnJ\nybEFLCI/sfPAEQa/msz8DXt47GetuLl7w0SHJHnAzGa5e1IsZaPd4HZ1lHW4+9snElym/bQlXPu4\nJJakICInZ+3OFG4ZNYPNew8z/KZO9G51WqJDknwoWlPSFcHvGkAP4Ivg/bnAf4GTSgxmVj/Yx83u\nvvxk9iUiOZu7fg+DR88k5M7YId3o1KBqokOSfCraDW6DAMzsM6Clu28O3tciPOIoKjMbB/QCqpnZ\nBuARoGSw7+HAw8CpwD/NDCAt1mqOiByfqUu2ctfY2VSvWJpXB3Xh9Oq6m1myF0vnc71jSSGwFcix\nB9jd++ew/jbgthiOLyInYez0dQx7dwGt61Tm5Vs6U72i7lGQ6GJJDFPN7FNgXPC+L/B5/EISkdzg\n7jw7ZTkvfLGSc5tX58UbOlK+dCz/5aWoi2VKjF8EHdFnBYtGuPs78Q1LRE7G0fQQD721gLdmb6Bf\n53o8cWVrPZtZYhbT14dgBNJJj0ISkfjbf/goP39jNt+u2MGvLmjGL89vQtCPJxKTaMNVv3P3nma2\nn+B5z8dWAe7uleIenYgcl637DjPolZks27qfp65ty/VJ9RIdkhRA0UYl9Qx+V8y7cETkRM1et5u7\nx85h98FURg3szDnNqic6JCmgotUYTom2obvvyv1wROR4pYec4V//wLNTllOrchkm3t6d1nX0xDU5\ncdH6GGYRbkLKqnHSgdPjEpGIxGzrvsP8asJc/vvDTi5rW4s/XdWGymU1O6qcnGhNSY3yMhAROT5f\nLN3K/ZPmcyg1nSevacP1SfXUySy5IqZRSWZWFWhKxie4fROvoEQke0fS0vnLx0t55T9rOKNWJV7o\n354mNdQVKLknx8RgZrcB9xB+XsJcoBvwPXBefEMTkcx+2H6Au8fOYfHmfQzs0ZCHLmlBmZLFEx2W\nFDKx1BjuAToD09z9XDNrAfwpvmGJSCR3Z9KsDTwyeRFlShZj5IAkLmhZM9FhSSEVS2I47O6HzQwz\nK+3uS82sedwjExEA9h0+yrB3FvLevE10O/0UnuvbgdMql8l5Q5ETFEti2BA8we1dYIqZ7QbWxjcs\nEQGYs243vxw/h017DnN/72bc2asJxfWkNYmzWOZKuip4+aiZfQlUBj6Ja1QiRVwo5Az/5gee/Ww5\nNSuVYeLt3ejUIOqtRSK5JtoNbh8BY4F33f0AgLt/nVeBiRRVa3ak8MCb85i5ZjeXtanFn67WvQmS\nt6LVGP4N9AP+FtQUxgEfuntqnkQmUsSEQs7o/67hqU+XUrJ4MZ65rh3XdKyjexMkz0W7wW0yMNnM\nyhF+zOcA4F9m9jEw1t2n5FGMIoXe2p0pPPDmfGas3sW5zavz56vbqoNZEiaWPoaDwARggpm1BV4l\nnCQ0eFrkJIVCzmvfr+HJT5ZRopjx9LVtubZTXdUSJKFiucGtJnA94WalWsBEYGB8wxIp/NbtPMgD\nb85j+updnNOsOn+5pg21KpdNdFgiUTufhwD9gebAW8AD7v7fvApMpLAKhZzXp6/lLx8vpbgZT13T\nluuSVEuQ/CNajaE78GdgqruHjnfHZjYKuBzY5u6ts1hvwPPApcBBYKC7zz7e44gUJOt3hWsJ01bt\n4qym1XjymrbUrqJaguQv0Tqfb828zMwedfdHY9z3aOBF4LVs1l9CeGK+pkBX4F/Bb5FCJxRy3pix\njj9/tIRiZvzl6jb07azZUCV/iml21Qg/Ax6NpaC7f2NmDaMU6QO85u4OTDOzKmZWy903H2dMIvna\nxj2HeGDSPP77w056NqnGk9e2pY5qCZKPHW9iyM2vN3WA9RHvNwTLfpIYzGwoMBSgfv36uRiCSHx9\nvXw794yfw9G0EH+6qg39u6iWIPnf8SaGTnGJIgfuPgIYAZCUlOSJiEHkeIRCzgtfrOS5qctpVqMi\nw2/uRKNq5RMdlkhMiuVUwMyeMrNKZlaS8CR6283splw49kagXsT7usEykQJtz8FUbn11Jn/7fDlX\ntq/DO3f1UFKQAiXHxAD0dvd9hEcYrQGaAA/kwrHfAwZYWDdgr/oXpKBbuHEvl7/wHf9ZuYM/XNma\nZ69vR7lSx1sxF0msWD6xx8pcBkxy972xtJGa2TigF1DNzDYAjwAlAdx9OPAR4aGqKwkPVx10vMGL\n5CcTZq7j95MXUa18KSbd0YP29aokOiSRExJLYvjAzJYCh4A7zaw6cDinjdy9fw7rHbgrpihF8rHD\nR9N5ePJCJiZv4Kym1Xi+XwdOKV8q0WGJnLBY5kp6yMyeItzUk25mKYSHmooUeet2HuTON2axaNM+\n7j6vCfde0EwP0pECL5a5kq4DPgmSwjCgI/AEsCXewYnkZ1OXbOVXE+YCMGpgEue10DOYpXCIpfP5\n9+6+38x6AhcALxO+S1mkSEoPOX/9bBmDX02mbtVyfHD3WUoKUqjE0seQHvy+DBjh7h+a2RNxjEkk\n39qVkso94+fw7YodXJ9Ul8f7tKZMSc1AL4VLLIlho5n9G7gQeNLMShNbTUOkUFmxdT+DRs9k2/4j\nPHlNG/p21l34UjjFkhiuBy4GnnH3PWZWi9y5j0GkwPhuxQ7ufH0WZUoVZ9Lt3WmnoahSiMX0BDcz\n+wG4yMwuAr5198/iH5pI/jB+xjqGvbuQJjUq8PLAzpoATwq9WKbEuAd4A6gR/LxuZnfHOzCRRAuF\nnD9/vISH3l7AmU2qMemO7koKUiTE0pQ0GOjq7ikAZvYk8D3wQjwDE0mkQ6np3DdxLh8v3MJN3erz\n6BWtKFFcXWtSNMSSGIwfRyYRvNYdPFJobdt/mCGvJjN/416GXXYGg3s20lTZUqTEkhheAaab2TvB\n+ysJ38sgUugs27KfW0fPZFdKKv++qRO9W52W6JBE8lwsnc/PmtlXQM9g0SB3nxPXqEQS4Jvl27nr\njdmULVWcibd3p03dyokOSSQhoiYGMysOLHL3FsDsvAlJJO+9MX0tD09eRNMaFRg1sDO11cksRVjU\nxBDMj7TMzOq7+7q8CkokrxwbefTSt6vp1bw6L97QkQql9fwEKdpi+R9QFVhkZjOAlGML3f1ncYtK\nJA8cTE3j3vFz+WzxVgZ0b8DDl7fUyCMRYksMv497FCJ5bHdKKgNfmcH8jXt5+PKWDDqzoUYeiQSy\nTQxm1gSo6e5fZ1reE9AjOKXA2nngCDeOnM6qHSkaeSSShWj15ueAfVks3xusEylwtu8/Qv+XprF6\nRwojByQpKYhkIVpTUk13X5B5obsvMLOGcYtIJE627jvMDS9NY9Oew7wysDM9mlRLdEgi+VK0GkO0\n6SNjGstnZhcHo5pWmtlDWayvbGbvm9k8M1tkZoNi2a/I8dq05xB9//09W/Ye5tVbuygpiEQRLTEk\nm9mQzAvN7DZgVk47Du6B+AdwCdAS6G9mLTMVuwtY7O7tgF7AX81MT1GXXLV+10H6jvienQdSeW1w\nF7o0OiXRIYnka9Gaku4F3jGzG/kxESQBpYCrYth3F2Clu68CMLPxQB9gcUQZBypaeDhIBWAXkHZc\nf4FIFOt2HqT/S9PYf/goY27rSns9R0EkR9kmBnffCvQws3OB1sHiD939ixj3XQdYH/F+A9A1U5kX\ngfeATUBFoK+7hzLvyMyGAkMB6tfXU7MkNqt3pHDDS9M4dDSdsUO60bqOprgQiUUscyV9CXwZp+Nf\nBMwFzgMpV/cgAAARK0lEQVQaA1PM7Ft3zzAayt1HACMAkpKSPE6xSCGyctsBbnhpGmkhZ+xt3WhZ\nu1KiQxIpMOJ5m+dGoF7E+7rBskiDgLc9bCWwGmgRx5ikCFi2ZT/9RnxPyGH8UCUFkeMVz8QwE2hq\nZo2CDuV+hJuNIq0Dzgcws5pAc2BVHGOSQm7xpn30f2kaxcwYP7QbzWpWTHRIIgVO3GYLc/c0M/sF\n8ClQHBjl7ovM7I5g/XDgD8BoM1tA+OE/D7r7jnjFJIXbwo17uenl6ZQtWZyxQ7rRqFr5RIckUiDF\ndRpJd/8I+CjTsuERrzcBveMZgxQNc9fvYcDL06lYpiTjhnSj/qnlEh2SSIGl+YWlQHN33py1gcff\nX0yV8uGkULeqkoLIyVBikAJr2Zb9DHt3ATPX7KZTg6q80L+DHrAjkguUGKTASTmSxvNTV/Dyd6up\nVKYET13Tlms71aVYMU2bLZIblBikwHB3Pl20lcffX8SmvYfpm1SPBy9pwSnlNYuKSG5SYpACYd3O\ngzzy3kK+XLadFqdV5IUbOtCpgeY8EokHJQbJ146kpTPi61W8+OVKShQzhl12BgN7NNQjOEXiSIlB\n8q3/rNzB7ycvZNX2FC5rU4thl59BrcrqXBaJNyUGyXe27TvMEx8u4b15m2hwajlGD+pMr+Y1Eh2W\nSJGhxCD5Rlp6iNe+X8vfpiznSFqIe85vyp29GlOmZPFEhyZSpCgxSL7w/Q87efS9RSzbup+zmlbj\n8T6tNaWFSIIoMUhCbd57iD9+uIQP5m+mbtWy/PvmTvRuWZPws5tEJBGUGCQhjqSlM/Lb1bz4xUpC\n7tx7QVPuOEfNRiL5gRKD5Lkvl27jsfcXsWbnQS5qVZNhl7Wk3ima30gkv1BikDyzdmcKj7+/mKlL\nt3F69fK8dmsXzm5WPdFhiUgmSgwSdwdT0/jnlz8w4ptVlCxu/O7SFgzs0YhSJXSTmkh+pMQgcePu\nfLRgC3/8cDGb9h7myva1+e2lZ1CzUplEhyYiUSgxSFz8sP0AD09eyH9W7uSMWpV4rl8HujTS3EYi\nBUGRSQypaSE+WrCZPu1rayhkHB0+ms4/vlzJv79eRemSxXi8Tytu6FJfcxuJFCBFJjG8PXsDD729\nAIArO9RJcDSF05dLt/HwewtZv+sQV3Wow+8uPYPqFUsnOiwROU5FJjFcl1SPSbM28PDkhXRvfKra\nuXPRpj2HeOz9RXy6aCuNq5dn7JCu9GhcLdFhicgJimv93swuNrNlZrbSzB7KpkwvM5trZovM7Ot4\nxVK8mPHMde1ITQ/x4Fvzcfd4HarIOJoeYsQ3P3DBs1/z9fLtPHBRcz6+52wlBZECLm41BjMrDvwD\nuBDYAMw0s/fcfXFEmSrAP4GL3X2dmcV1Cs1G1crz0MUtePT9xUxMXk/fzvXjebhCbeaaXQx7ZyHL\ntu7ngjNq8MgVrXSTmkghEc+mpC7ASndfBWBm44E+wOKIMjcAb7v7OgB33xbHeAAY0L1h8HjIxfRo\nXE0Xs+O088AR/vLxUibN2kCdKmUZcXMnerc6LdFhiUguimdTUh1gfcT7DcGySM2Aqmb2lZnNMrMB\nWe3IzIaaWbKZJW/fvv2kgipWzHjq2rYA/ObN+YRCalKKRSjkjJuxjvP++jXvzNnIHec0Zsp9Zysp\niBRCie58LgF0As4HygLfm9k0d18eWcjdRwAjAJKSkk76Sl7vlHL8/vKWPPT2AsZMW8stPRqe7C4L\nrQNH0vhg3iZen76WhRv30bXRKTxxZWua1qyY6NBEJE7imRg2AvUi3tcNlkXaAOx09xQgxcy+AdoB\ny4mzvp3r8cmiLfz54yWc3ay65v6P4O7MWb+HCTPW8/78TRxMTadpjQr89bp2XN2xju4DESnk4pkY\nZgJNzawR4YTQj3CfQqTJwItmVgIoBXQF/hbHmP7HzPjL1W3p/bevuX/SPCbe3p3ixYr2BW9XSirv\nzNnIhJnrWL71AGVLFueKdrXo27k+HetXUUIQKSLilhjcPc3MfgF8ChQHRrn7IjO7I1g/3N2XmNkn\nwHwgBIx094Xxiimz0yqX4bE+rfjVhHmM/HYVt5/TOK8OnW+EQs5/ftjB+JnrmbJoK6npIdrXq8Kf\nr27D5W1rUbFMyUSHKCJ5zAraeP6kpCRPTk7Otf25O3e8Posvl27ng1/2pFkRaTvfvPcQk5I3MDF5\nPRt2H6JKuZJc1aEOfTvXo8VplRIdnojkMjOb5e5JMZUt6okBYMeBI/T+2zfUqVKWt3/eg5KFdF6f\n1LQQU5dsZULyer5Zvp2QQ88m1bi+cz16t6ypp6eJFGLHkxgSPSopX6hWoTRPXNman78xm3999QO/\nPL9pokPKVSu27mfCzPW8M2cjO1NSOa1SGe46twnXdapH/VN1H4eIZKTEELi0TS1+1q42f5+6gvNa\n1KB1ncqJDumkHBtmOiF5PXPW7aFEMeOCM2rSt3M9zm5Wvch3tItI9pQYIjzepxXfr9rJ/ZPmMfkX\nZ1K6RMFqWnF3Zq3dzYSZ6/lwwWYOpqbTpEYF/u/SM7iqYx2qVdBMpyKSMyWGCFXKleLJa9pw6+hk\nnv98Bb+5uEWiQ4rJ9v1HeHt2uCP5h+0plC9VnCva1ub6zvU0zFREjpsSQybntajJ9Ul1Gf71D1zY\nsiYd6ldNdEhZ2nf4KF8t284H8zbxxdJtpIWcpAZVeeqaxlzWthblS+ufVkROjK4eWRh2eUu+W7GD\nX0+ax0e/PCvfjNbZsvcwU5Zs5bNFW5i2aidH051qFUpza89GXJ9UjyY1KiQ6RBEpBJQYslCpTEme\nurYdN708nac/XcbvL2+ZkDjcnRXbDvDZoi1MWbyVeRv2AuHpw289s9H/ajTqSBaR3KTEkI2eTatx\nc7cGjPrPavYfPkqnBlXpWL8qjatXoFgcL8TpoXAH8pTFW/hs8VbW7jwIQPt6VXjgouZc1KomjatX\nUL+BiMSNEkMUv720BbsOpvLJwi1MTN4AQKUyJWhfvyod61ehY/2qtK9fhUonOG3E4aPpbN57mA27\nD7Jx9yFmrd3NF0u3sTMllVLFi9G98akMPft0Ljijph5FKiJ5RokhinKlSvCPGzoSCjmrdqQwe91u\n5qzbzey1e3h+6grcwQya1qhAx/rhGkWH+lX+V6s4mJrGxt2H2LDnEBt2Hwq/3n2QjcH77fuPZDhe\nxTIlOK9FDS5sWZNzmlXXPEUikhCaEuME7Tt8lPnr9zJ73e4gYexh76GjQLhWUaJ4MXalpGbYpmRx\no3aVstStWpY6VcpSt2o56lQpS52q4WWnVSpDiUI6HYeIJJamxMgDlcqUpGfTavRsGn7w/bFaxZx1\nu5mzfg9AcPE/lgjKUaNi6bj2T4iI5AYlhlxSrJjRpEYFmtSowHVJ9XLeQEQkn1K7hYiIZKDEICIi\nGSgxiIhIBkoMIiKSgRKDiIhkENfEYGYXm9kyM1tpZg9FKdfZzNLM7Np4xiMiIjmLW2Iws+LAP4BL\ngJZAfzP7yWx0Qbkngc/iFYuIiMQunjWGLsBKd1/l7qnAeKBPFuXuBt4CtsUxFhERiVE8b3CrA6yP\neL8B6BpZwMzqAFcB5wKds9uRmQ0FhgZvD5jZMqAasCM3Ay6gdB5+pHMRpvMQpvMQduw8NIh1g0Tf\n+fwc8KC7h6JNI+3uI4ARkcvMLDnWeT8KM52HH+lchOk8hOk8hJ3IeYhnYtgIRM4NUTdYFikJGB8k\nhWrApWaW5u7vxjEuERGJIp6JYSbQ1MwaEU4I/YAbIgu4e6Njr81sNPCBkoKISGLFLTG4e5qZ/QL4\nFCgOjHL3RWZ2R7B++EkeYkTORYoEnYcf6VyE6TyE6TyEHfd5KHDPYxARkfjSnc8iIpKBEoOIiGRQ\nIBNDrFNtFHZmtsbMFpjZXDNL/PNO84iZjTKzbWa2MGLZKWY2xcxWBL+rJjLGvJDNeXjUzDYGn4m5\nZnZpImPMC2ZWz8y+NLPFZrbIzO4Jlhepz0SU83Dcn4kC18cQTKGxHLiQ8E1zM4H+7r44oYElgJmt\nAZLcvUjdxGNmZwMHgNfcvXWw7Clgl7v/JfiyUNXdH0xknPGWzXl4FDjg7s8kMra8ZGa1gFruPtvM\nKgKzgCuBgRShz0SU83A9x/mZKIg1hlin2pBCyt2/AXZlWtwHeDV4/Srh/xCFWjbnochx983uPjt4\nvR9YQnjmhSL1mYhyHo5bQUwMWU21cUJ/fCHgwOdmNiuYNqQoq+num4PXW4CaiQwmwe42s/lBU1Oh\nbj7JzMwaAh2A6RThz0Sm8wDH+ZkoiIlBftTT3dsTnsH2rqBpocjzcPtowWojzT3/Ak4H2gObgb8m\nNpy8Y2YVCE/Iea+774tcV5Q+E1mch+P+TBTExBDLVBtFgrtvDH5vA94h3MxWVG0N2liPtbUWydl6\n3X2ru6e7ewh4iSLymTCzkoQvhm+4+9vB4iL3mcjqPJzIZ6IgJob/TbVhZqUIT7XxXoJjynNmVj7o\nYMLMygO9gYXRtyrU3gNuCV7fAkxOYCwJc+xCGLiKIvCZsPBkay8DS9z92YhVReozkd15OJHPRIEb\nlQQQDLd6jh+n2vhjgkPKc2Z2OuFaAoSnNhlbVM6DmY0DehGeeHEr8AjwLjARqA+sBa5390LdMZvN\neehFuMnAgTXA7RHt7IWSmfUEvgUWAKFg8e8It68Xmc9ElPPQn+P8TBTIxCAiIvFTEJuSREQkjpQY\nREQkAyUGERHJQIlBREQyUGIQEZEMlBgk3wlmiLwo07J7zexfOWx3IM5xVTez6WY2x8zOyrTuKzNL\nCl43Cmb0vCiLfTwdzHz59AnG0MvMPoh4/4SZfWJmpYMYkiPWJZnZVxHbuZldEbH+AzPrdSJxSOGm\nxCD50TjCNy5G6hcsT6TzgQXu3sHdv82qgJnVBT4Bfu3un2ZRZCjQ1t0fiOWAZpbt43fNbBhwJnCV\nux8JFtcws0uy2WQD8H+xHFeKNiUGyY/eBC4L7mw/NiFYbeBbM6tgZlPNbLaFn0Xxk5l1s/hW/aKZ\nDQxedzKzr4OJBz/NdFfosfINzeyLYNKxqWZW38zaA08BfYI57ctmEXct4DPg/9z9J3fjm9l7QAVg\nlpn1zeo4QbnRZjbczKYHx/wJM/s14TmyrnD3QxGrnib7i/88YK+ZXZjNehFAiUHyoeDu1BmEL3wQ\nri1MDCZCO0z4G3JH4Fzgr8FUADkK5pF5AbjW3TsBo4Cs7hZ/AXjV3dsCbwB/d/e5wMPABHdvn+li\nfMyrwIvu/mY2f9fPgEPB9hOyOk5E8bpAD3e/L4tdnQncAVzi7pmbz74HUs3s3KxiCP7eYdmsEwGU\nGCT/imxOimxGMuBPZjYf+JzwlOuxTqfcHGgNTDGzuYQvkHWzKNcdGBu8HgP0jHH/nwM3mVm5GMtH\nO84kd0/PZruVhM9Ddt/8nyCbi3/wDIdj0yeIZEmJQfKrycD5ZtYRKOfus4LlNwLVgU7BlONbgTKZ\ntk0j42f72HoDFgXf2Nu7ext3752LMT9FeJLHSdH6BmKUEmXdVuBS4Lmsagbu/gVQFuiWzfaqNUhU\nSgySLwVNJF8Sbu6J7HSuDGxz96PBRbFBFpuvBVoGI3WqEO40BlgGVDez7hBuWjKzVlls/19+rK3c\nSHhisljdC+wDXo6hieuEj+Puy4GrgdeD/o/MngB+k822nwFVgbaxHk+KFiUGyc/GAe3ImBjeAJLM\nbAEwAFiaeSN3X094Vs2Fwe85wfJU4FrgSTObB8wFemRx3LuBQUFz1c3APbEGHPSD3EK4IzrLjuPc\nOE5wrJnAIOA9M2ucad1HwPYom/+RjM81Efkfza4qIiIZqMYgIiIZKDGIiEgGSgwiIpKBEoOIiGSg\nxCAiIhkoMYiISAZKDCIiksH/A6lzeRkN44sXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1142c9da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the value of K for KNN (x-axis) versus the cross-validated accuracy (y-axis)\n",
    "#remember that MSE is negative due to implementation in python, so we take negative\n",
    "k_scores = np.negative(k_scores)\n",
    "plt.plot(k_range, k_scores)\n",
    "plt.xlabel('Value of K for KNN')\n",
    "plt.ylabel('Cross-Validated MSE')\n",
    "plt.title('Neighbours vs MSE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decision Tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the parameters and distributions to sample from: param_dist\n",
    "param_dist = {\"max_depth\": [3, None],\n",
    "              \"min_samples_leaf\": range(1, 9),\n",
    "              \"criterion\": [\"mse\", \"mae\"]}\n",
    "\n",
    "# Instantiate a Decision Tree regressor: tree\n",
    "tree = DecisionTreeRegressor(random_state=450411920)\n",
    "\n",
    "# Instantiate the RandomizedSearchCV object: tree_cv\n",
    "tree_cv = RandomizedSearchCV(tree,param_dist, cv=10, scoring='neg_mean_squared_error', n_iter=10, random_state=5)\n",
    "\n",
    "#Fitting the model and testing out random hyperparameters. This saves on computation time\n",
    "treeopt = tree_cv.fit(final_train, y_train)\n",
    "treepred = treeopt.predict(final_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of this model is 0.375\n",
      "Best parameters from this is  {'min_samples_leaf': 3, 'max_depth': None, 'criterion': 'mse'}\n",
      "This is the best Decision Tree estimator DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,\n",
      "           max_leaf_nodes=None, min_impurity_decrease=0.0,\n",
      "           min_impurity_split=None, min_samples_leaf=3,\n",
      "           min_samples_split=2, min_weight_fraction_leaf=0.0,\n",
      "           presort=False, random_state=450411920, splitter='best')\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print('MSE of this model is %.3f'%np.negative(tree_cv.best_score_))\n",
    "print('Best parameters from this is ',tree_cv.best_params_)\n",
    "print('This is the best Decision Tree estimator', tree_cv.best_estimator_)\n",
    "models.append(tree_cv)\n",
    "modelname.append('Decision Tree')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the parameters and distributions to sample from: param_dist\n",
    "param_dist = {\"max_depth\": [3, None],\n",
    "              \"min_samples_leaf\": range(4, 9),\n",
    "              \"criterion\": [\"mse\", \"mae\"]}\n",
    "\n",
    "# Instantiate a Extremely Random Forest regressor: randext\n",
    "randFor = RandomForestRegressor(random_state=450411920)\n",
    "\n",
    "# Instantiate the RandomizedSearchCV object: randFor_cv\n",
    "randFor_cv = RandomizedSearchCV(randFor,param_dist, cv=10, scoring='neg_mean_squared_error', n_iter=10, random_state=5)\n",
    "\n",
    "#Fitting the model and testing out random hyperparameters. This saves on computation time\n",
    "optrand = randFor_cv.fit(final_train, y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rfpred = optrand.predict(final_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of this model is 0.347\n",
      "Best parameters from this is  {'min_samples_leaf': 4, 'max_depth': None, 'criterion': 'mse'}\n",
      "This is the best Random Forest estimator RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,\n",
      "           max_features='auto', max_leaf_nodes=None,\n",
      "           min_impurity_decrease=0.0, min_impurity_split=None,\n",
      "           min_samples_leaf=4, min_samples_split=2,\n",
      "           min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,\n",
      "           oob_score=False, random_state=450411920, verbose=0,\n",
      "           warm_start=False)\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print('MSE of this model is %.3f'%np.negative(randFor_cv.best_score_))\n",
    "print('Best parameters from this is ',randFor_cv.best_params_)\n",
    "print('This is the best Random Forest estimator', randFor_cv.best_estimator_)\n",
    "models.append(randFor_cv)\n",
    "modelname.append('Random Forest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient Boosting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the parameters and distributions to sample from: param_dist\n",
    "param_dist = {\"loss\": [\"huber\", \"ls\"],\n",
    "              \"n_estimators\": range(300,550,50),\n",
    "              \"learning_rate\": [0.1],\n",
    "              \"min_samples_split\" : [2,4],\n",
    "              \"min_samples_leaf\" : range(1,3)\n",
    "              }\n",
    "\n",
    "# Instantiate a GB boosting classifier: randGB\n",
    "randGB = GradientBoostingRegressor()\n",
    "\n",
    "# Instantiate the RandomizedSearchCV object: randGB_cv\n",
    "randGB_cv = RandomizedSearchCV(randGB,param_dist, cv=10, scoring='neg_mean_squared_error', n_iter=10, random_state=5)\n",
    "\n",
    "#Fitting the model and testing out random hyperparameters. This saves on computation time\n",
    "optgb = randGB_cv.fit(final_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gbpred = treeopt.predict(final_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of this model is 0.266\n",
      "Best parameters from this is  {'n_estimators': 500, 'min_samples_split': 4, 'min_samples_leaf': 2, 'loss': 'ls', 'learning_rate': 0.1}\n",
      "This is the best GradientBoost estimator GradientBoostingRegressor(alpha=0.9, criterion='friedman_mse', init=None,\n",
      "             learning_rate=0.1, loss='ls', max_depth=3, max_features=None,\n",
      "             max_leaf_nodes=None, min_impurity_decrease=0.0,\n",
      "             min_impurity_split=None, min_samples_leaf=2,\n",
      "             min_samples_split=4, min_weight_fraction_leaf=0.0,\n",
      "             n_estimators=500, presort='auto', random_state=None,\n",
      "             subsample=1.0, verbose=0, warm_start=False)\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print('MSE of this model is %.3f'%np.negative(randGB_cv.best_score_))\n",
    "print('Best parameters from this is ',randGB_cv.best_params_)\n",
    "print('This is the best GradientBoost estimator', randGB_cv.best_estimator_)\n",
    "models.append(randGB_cv)\n",
    "modelname.append('Gradient Boost')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extreme Gradient Boosting\n",
    "\n",
    "https://cambridgespark.com/content/tutorials/hyperparameter-tuning-in-xgboost/index.html\n",
    "\n",
    "http://xgboost.readthedocs.io/en/latest/parameter.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['base_score', 'colsample_bylevel', 'colsample_bytree', 'gamma', 'learning_rate', 'max_delta_step', 'max_depth', 'min_child_weight', 'missing', 'n_estimators', 'nthread', 'objective', 'reg_alpha', 'reg_lambda', 'scale_pos_weight', 'seed', 'silent', 'subsample'])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xgb_model.get_params().keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 3 folds for each of 180 candidates, totalling 540 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done 540 out of 540 | elapsed:  8.9min finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.999575695547\n",
      "{'learning_rate': 0.1, 'max_depth': 4, 'min_child_weight': 4, 'n_estimators': 400}\n"
     ]
    }
   ],
   "source": [
    "xgb_model = xgb.XGBRegressor()\n",
    "clf = GridSearchCV(xgb_model,\n",
    "                   {'max_depth': range(3,6),\n",
    "                     'min_child_weight': range(1,4),\n",
    "                     \"n_estimators\": range(300,550,50),\n",
    "                     \"learning_rate\": [0.1,0.3,0.5],\n",
    "                      \"min_child_weight\" : [1,2,3,4],},\n",
    "                     verbose=1)\n",
    "clf.fit(final_train,y_train)\n",
    "print(clf.best_score_)\n",
    "print(clf.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xgpred = clf.predict(final_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LASSO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#LASSO\n",
    "lasso = Lasso(alpha = 1)\n",
    "optlasso = lasso.fit(final_train, np.ravel(y_train)) \n",
    "modelname.append('LASSO')\n",
    "models.append(optlasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rg = Ridge(alpha = 10)\n",
    "rigopt = rg.fit(final_train, y_train)\n",
    "modelname.append('Ridge')\n",
    "models.append(rigopt)\n",
    "\n",
    "#Need to standardise features later on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Generates us a table for results with hyperparameter optimised models\n",
    "def getResultTable(rows, modelsUsed):\n",
    "    columns=['R-Square', 'MSE',]\n",
    "    results=pd.DataFrame(0.0, columns=columns, index=rows)\n",
    "    predictions = []\n",
    "    for clf in modelsUsed:\n",
    "        pred = clf.predict(final_test)\n",
    "        predictions.append(pred)\n",
    "        \n",
    "    for row,pred in zip(range(0,len(rows)),predictions):\n",
    "        \n",
    "        results.iloc[row,0] = r2_score(y_test,pred)\n",
    "        results.iloc[row,1] = mean_squared_error(y_test,pred)\n",
    "        \n",
    "    return results.round(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "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>R-Square</th>\n",
       "      <th>MSE</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>KNN</th>\n",
       "      <td>0.999</td>\n",
       "      <td>0.575</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Decision Tree</th>\n",
       "      <td>0.999</td>\n",
       "      <td>0.371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>0.999</td>\n",
       "      <td>0.344</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Random Forest</th>\n",
       "      <td>0.999</td>\n",
       "      <td>0.344</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Gradient Boost</th>\n",
       "      <td>1.000</td>\n",
       "      <td>0.263</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LASSO</th>\n",
       "      <td>0.978</td>\n",
       "      <td>14.019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Ridge</th>\n",
       "      <td>0.990</td>\n",
       "      <td>6.129</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                R-Square     MSE\n",
       "KNN                0.999   0.575\n",
       "Decision Tree      0.999   0.371\n",
       "Random Forest      0.999   0.344\n",
       "Random Forest      0.999   0.344\n",
       "Gradient Boost     1.000   0.263\n",
       "LASSO              0.978  14.019\n",
       "Ridge              0.990   6.129"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getResultTable(modelname,models)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-square of XGboost:  1.0\n",
      "MSE of XGBoost  0.266\n"
     ]
    }
   ],
   "source": [
    "print(\"R-square of XGboost: \",r2_score(y_test,xgpred).round(3))\n",
    "print(\"MSE of XGBoost \", mean_squared_error(y_test,xgpred).round(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "final = (gbpred+rfpred+xgpred)/3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-square of XGboost:  1.0\n",
      "MSE of XGBoost  0.3\n"
     ]
    }
   ],
   "source": [
    "print(\"R-square of XGboost: \",r2_score(y_test,final).round(3))\n",
    "print(\"MSE of XGBoost \", mean_squared_error(y_test,final).round(3))"
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}