{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from sklearn import tree, model_selection, datasets, metrics, ensemble, linear_model, neighbors\n",
"import graphviz as gv\n",
"import seaborn as sns"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"7. In the lab, we applied random forests to the Boston data using mtry=6 and using ntree=25 and ntree=500. Create a plot displaying the test error resulting from random forests on this data set for a more comprehensive range of values for mtry and ntree. You can model your plot after Figure 8.10. Describe the results obtained."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"df = pd.DataFrame({'sqrt_rmses' : sqrt_rmses,\n",
" 'log2_rmses' : log2_rmses,\n",
" 'all_fx_rmses': all_fx_rmses,\n",
" 'trees' : rnge})\n",
"sns.lineplot(x='trees', y='sqrt_rmses', data=df, color='tab:blue')\n",
"sns.lineplot(x='trees', y='log2_rmses', data=df, color='tab:green')\n",
"sns.lineplot(x='trees', y='all_fx_rmses', data=df, color='tab:red')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"8. In the lab, a classification tree was applied to the Carseats data set after converting Sales into a qualitative response variable. Now we will seek to predict Sales using regression trees and related approaches, treating the response as a quantitative variable."
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Sales
\n",
"
CompPrice
\n",
"
Income
\n",
"
Advertising
\n",
"
Population
\n",
"
Price
\n",
"
ShelveLoc
\n",
"
Age
\n",
"
Education
\n",
"
Urban
\n",
"
US
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
9.50
\n",
"
138
\n",
"
73
\n",
"
11
\n",
"
276
\n",
"
120
\n",
"
0
\n",
"
42
\n",
"
17
\n",
"
True
\n",
"
True
\n",
"
\n",
"
\n",
"
1
\n",
"
11.22
\n",
"
111
\n",
"
48
\n",
"
16
\n",
"
260
\n",
"
83
\n",
"
2
\n",
"
65
\n",
"
10
\n",
"
True
\n",
"
True
\n",
"
\n",
"
\n",
"
2
\n",
"
10.06
\n",
"
113
\n",
"
35
\n",
"
10
\n",
"
269
\n",
"
80
\n",
"
1
\n",
"
59
\n",
"
12
\n",
"
True
\n",
"
True
\n",
"
\n",
"
\n",
"
3
\n",
"
7.40
\n",
"
117
\n",
"
100
\n",
"
4
\n",
"
466
\n",
"
97
\n",
"
1
\n",
"
55
\n",
"
14
\n",
"
True
\n",
"
True
\n",
"
\n",
"
\n",
"
4
\n",
"
4.15
\n",
"
141
\n",
"
64
\n",
"
3
\n",
"
340
\n",
"
128
\n",
"
0
\n",
"
38
\n",
"
13
\n",
"
True
\n",
"
False
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Sales CompPrice Income Advertising Population Price ShelveLoc Age \\\n",
"0 9.50 138 73 11 276 120 0 42 \n",
"1 11.22 111 48 16 260 83 2 65 \n",
"2 10.06 113 35 10 269 80 1 59 \n",
"3 7.40 117 100 4 466 97 1 55 \n",
"4 4.15 141 64 3 340 128 0 38 \n",
"\n",
" Education Urban US \n",
"0 17 True True \n",
"1 10 True True \n",
"2 12 True True \n",
"3 14 True True \n",
"4 13 True False "
]
},
"execution_count": 116,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"car_df = pd.read_csv('carseats.csv')\n",
"car_df = car_df.drop(car_df.columns[0], axis=1)\n",
"car_df['Urban'] = car_df['Urban'] == 'Yes'\n",
"car_df['US'] = car_df['US'] == 'Yes'\n",
"car_df['ShelveLoc'] = car_df['ShelveLoc'].map({'Bad' : 0, 'Medium': 1, 'Good' : 2})\n",
"car_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) Split the data set into a training set and a test set."
]
},
{
"cell_type": "code",
"execution_count": 118,
"metadata": {},
"outputs": [],
"source": [
"train, test = model_selection.train_test_split(car_df, test_size=0.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) Fit a regression tree to the training set. Plot the tree, and interpret the results. What test MSE do you obtain?"
]
},
{
"cell_type": "code",
"execution_count": 119,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"test rmse: 2.2041731241021694\n"
]
}
],
"source": [
"reg = tree.DecisionTreeRegressor(max_depth=3) \n",
"train_x = train.drop('Sales', axis=1)\n",
"train_y = train.Sales\n",
"reg.fit(train_x,train_y)\n",
"test_x = test.drop('Sales', axis=1)\n",
"test_y = test.Sales\n",
"preds = reg.predict(test_x)\n",
"rmse = np.sqrt(metrics.mean_squared_error(test_y, preds))\n",
"print('test rmse: {}'.format(rmse))"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dot_dat = tree.export_graphviz(reg, out_file=None, rotate=True)\n",
"graph = gv.Source(dot_dat)\n",
"graph"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) Use cross-validation in order to determine the optimal level of tree complexity. Does pruning the tree improve the test MSE?\n",
"\n",
".. omitting this as tree pruning isn't easily available in python world, the python community prefer to control variance with boosting, bagging and random forest."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(d) Use the bagging approach in order to analyze this data. What test MSE do you obtain? Use the importance() function to de- termine which variables are most important."
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"training rmse: 1.5726550643736208\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4IAAAHICAYAAAAbVVRCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8bXVdL/zPl014SY9p7MoLuNHQQk3ULVaamoeMnlKsNEEzPI8nsiN28fSc43nsUaTj69Gup5OampG3CgXTtkYaeSG0TDaKICiJeIOsUMxLKgR+zx9zbJku1t6sDWvstff+vd+v13qtOX7jMr/zt8YaY37muMzq7gAAADCOAza6AAAAAPYsQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwmAM3uoD1cvDBB/eWLVs2ugwAAIANcd55532muzevZdr9Jghu2bIl27dv3+gyAAAANkRVfWKt0zo1FAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYA7c6AL2pAf8P6/a6BL2Guf9xs9sdAkAAMAGcUQQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYzKxBsKqOqapLqurSqnrmKuOfWlUXVtX5VfWuqjpiadz/mOa7pKp+eM46AQAARjJbEKyqTUlelORHkhyR5PjloDf5k+6+T3cfmeTXk/z2NO8RSY5Lcq8kxyR58bQ8AAAAbqY5jwgeleTS7r6su69JclqSY5cn6O4vLA1+c5KeHh+b5LTuvrq7P5bk0ml5AAAA3EwHzrjsOyf51NLw5UketHKiqnpakmckOSjJI5bmfc+Kee88T5kAAABj2fCbxXT3i7r77kn+e5Jf3Z15q+rEqtpeVduvvPLKeQoEAADYz8wZBK9IcsjS8F2mtp05Lcljdmfe7n5Zd2/t7q2bN2++meUCAACMYc4geG6Sw6vqsKo6KIubv2xbnqCqDl8a/NEkH5keb0tyXFXdoqoOS3J4kvfOWCsAAMAwZrtGsLuvraqTkrw1yaYkp3b3RVV1SpLt3b0tyUlVdXSSf0/yuSQnTPNeVFWvS3JxkmuTPK27r5urVgAAgJHMebOYdPeZSc5c0fbspce/uIt5n5fkefNVBwAAMKYNv1kMAAAAe5YgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAM5sCNLoB91ydPuc9Gl7BXOPTZF250CQAAsFscEQQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwmFmDYFUdU1WXVNWlVfXMVcY/o6ourqoLquptVXXXpXHXVdX508+2OesEAAAYyYFzLbiqNiV5UZIfSnJ5knOralt3X7w02fuTbO3uL1fVzyf59SSPn8Z9pbuPnKs+AACAUc15RPCoJJd292XdfU2S05IcuzxBd7+ju788Db4nyV1mrAcAAIDMGwTvnORTS8OXT20785Qkf7k0fMuq2l5V76mqx8xRIAAAwIhmOzV0d1TVTyfZmuRhS8137e4rqupuSd5eVRd290dXzHdikhOT5NBDD91j9QIAAOzL5jwieEWSQ5aG7zK1fYOqOjrJs5I8uruv3tHe3VdMvy9L8s4k91s5b3e/rLu3dvfWzZs3r2/1AAAA+6k5g+C5SQ6vqsOq6qAkxyX5hrt/VtX9krw0ixD4L0vtt6+qW0yPD07y4CTLN5kBAADgJprt1NDuvraqTkry1iSbkpza3RdV1SlJtnf3tiS/keQ2SU6vqiT5ZHc/Osl3J3lpVX0ti7D6/BV3GwUAAOAmmvUawe4+M8mZK9qevfT46J3M97dJ7jNnbQAAAKOa9QvlAQAA2PsIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYzaxCsqmOq6pKqurSqnrnK+GdU1cVVdUFVva2q7ro07oSq+sj0c8KcdQIAAIxktiBYVZuSvCjJjyQ5IsnxVXXEisnen2Rrd39PkjOS/Po07x2SPCfJg5IcleQ5VXX7uWoFAAAYyZxHBI9Kcml3X9bd1yQ5LcmxyxN09zu6+8vT4HuS3GV6/MNJzuruq7r7c0nOSnLMjLUCAAAMY84geOckn1oavnxq25mnJPnLmzgvAAAAa3TgRheQJFX100m2JnnYbs53YpITk+TQQw+doTIAAID9z5xHBK9IcsjS8F2mtm9QVUcneVaSR3f31bszb3e/rLu3dvfWzZs3r1vhAAAA+7M5g+C5SQ6vqsOq6qAkxyXZtjxBVd0vyUuzCIH/sjTqrUkeWVW3n24S88ipDQAAgJtptlNDu/vaqjopiwC3Kcmp3X1RVZ2SZHt3b0vyG0luk+T0qkqST3b3o7v7qqr6tSzCZJKc0t1XzVUrAADASGa9RrC7z0xy5oq2Zy89PnoX856a5NT5qgMAABjTrF8oDwAAwN5HEAQAABiMIAgAADAYQRAAAGAwgiAAAMBg1hwEq+qu05e/p6puVVW3na8sAAAA5rKmIFhVP5vkjCy+/D1J7pLkjXMVBQAAwHzWekTwaUkenOQLSdLdH0nybXMVBQAAwHzWGgSv7u5rdgxU1YFJep6SAAAAmNNag+DZVfX/JrlVVf1QktOTvGm+sgAAAJjLWoPgM5NcmeTCJD+X5MwkvzpXUQAAAMznwDVOd6skp3b3HyRJVW2a2r48V2EAAADMY61HBN+WRfDb4VZJ/nr9ywEAAGBuaw2Ct+zuL+0YmB7fep6SAAAAmNNag+C/VdX9dwxU1QOSfGWekgAAAJjTWq8R/KUkp1fVPyapJN+R5PGzVQUAAMBs1hQEu/vcqvquJPecmi7p7n+frywAAADmstYjgknywCRbpnnuX1Xp7lfNUhUAAACzWVMQrKpXJ7l7kvOTXDc1dxJBEAAAYB+z1iOCW5Mc0d09ZzEAAADMb613Df1gFjeIAQAAYB+31iOCBye5uKrem+TqHY3d/ehZqgIAAGA2aw2CJ89ZBAAAAHvOWr8+4uy5CwEAAGDPWNM1glX1vVV1blV9qaquqarrquoLcxcHAADA+lvrzWJemOT4JB9Jcqsk/znJi+YqCgAAgPmsNQimuy9Nsqm7r+vuP0pyzHxlAQAAMJe13izmy1V1UJLzq+rXk3w6uxEiAQAA2HusNcw9aZr2pCT/luSQJD8xV1EAAADMZ61B8DHd/dXu/kJ3P7e7n5Hkx+YsDAAAgHmsNQiesErbk9exDgAAAPaQXV4jWFXHJ3lCkrtV1balUbdNctWchQEAADCPG7tZzN9mcWOYg5P81lL7F5NcMFdRAAAAzGeXQbC7P1FVlyf5anefvYdqAgAAYEY3eo1gd1+X5GtVdbs9UA8AAAAzW+v3CH4pyYVVdVYWXx+RJOnuX5ilKgAAAGaz1iD4Z9MPAAAA+7g1BcHufmVVHZTkHlPTJd397/OVBQAAwFzWFASr6uFJXpnk40kqySFVdUJ3/818pQEAADCHtZ4a+ltJHtndlyRJVd0jyZ8mecBchQEAADCPG71r6OSbdoTAJOnuf0jyTfOUBAAAwJzWekRwe1W9PMlrpuEnJtk+T0kAAADMaa1B8OeTPC3Jjq+LOCfJi2epCAAAgFmt9a6hV1fVC5O8LcnXsrhr6DWzVgYAAMAs1nrX0B9N8pIkH83irqGHVdXPdfdfzlkcAAAA62937hr6g919aZJU1d2T/EUSQRAAAGAfs9a7hn5xRwicXJbkizPUAwAAwMx2566hZyZ5XZJO8rgk51bVTyRJd//ZTPUBAACwztYaBG+Z5J+TPGwavjLJrZI8KotgKAgCAADsI9Z619D/NHchAAAA7BlrvWvoYUmenmTL8jzd/eh5ygIAAGAuaz019I1J/jDJm7L4HsE1qapjkvxukk1JXt7dz18x/qFJ/leS70lyXHefsTTuuiQXToOfFDoBAADWx1qD4Fe7+3/vzoKralOSFyX5oSSXZ3FzmW3dffHSZJ9M8uQkv7LKIr7S3UfuznMCAABw49YaBH+3qp6T5K+SXL2jsbvft4t5jkpyaXdfliRVdVqSY5N8PQh298encWs+yggAAMDNs9YgeJ8kT0ryiFx/amhPwztz5ySfWhq+PMmDdqO2W1bV9iTXJnl+d79xN+YFAABgJ9YaBB+X5G7dfc2cxaxw1+6+oqruluTtVXVhd390eYKqOjHJiUly6KGH7sHSAAAA9l0HrHG6Dyb5lt1c9hVJDlkavsvUtibdfcX0+7Ik70xyv1WmeVl3b+3urZs3b97N8gAAAMa01iOC35Lkw1V1br7xGsFd3cnz3CSHT189cUWS45I8YS1PVlW3T/Ll7r66qg5O8uAkv77GWgEAANiFtQbB5+zugrv72qo6Kclbs/j6iFO7+6KqOiXJ9u7eVlUPTPKGJLdP8qiqem533yvJdyd56XQTmQOyuEbw4p08FQAAALthTUGwu8++KQvv7jOTnLmi7dlLj8/N4pTRlfP9bRY3qAEAAGCd7TIIVtW7uvshVfXFLO4S+vVRSbq7/8Os1QEAALDudhkEu/sh0+/b7plyAAAAmNta7xoKAADAfkIQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADGbWIFhVx1TVJVV1aVU9c5XxD62q91XVtVX12BXjTqiqj0w/J8xZJwAAwEhmC4JVtSnJi5L8SJIjkhxfVUesmOyTSZ6c5E9WzHuHJM9J8qAkRyV5TlXdfq5aAQAARjLnEcGjklza3Zd19zVJTkty7PIE3f3x7r4gyddWzPvDSc7q7qu6+3NJzkpyzIy1AgAADGPOIHjnJJ9aGr58apt7XgAAAHZhn75ZTFWdWFXbq2r7lVdeudHlAAAA7BPmDIJXJDlkafguU9u6zdvdL+vurd29dfPmzTe5UAAAgJHMGQTPTXJ4VR1WVQclOS7JtjXO+9Ykj6yq2083iXnk1AYAAMDNdOBcC+7ua6vqpCwC3KYkp3b3RVV1SpLt3b2tqh6Y5A1Jbp/kUVX13O6+V3dfVVW/lkWYTJJTuvuquWqFjfbg33vwRpew13j309+90SUAAOz3ZguCSdLdZyY5c0Xbs5cen5vFaZ+rzXtqklPnrA8AAGBE+/TNYgAAANh9giAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwRy40QUArLezH/qwjS5hr/Gwvzl7o0sAAPZCjggCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYNw1FIBdeuF/fdNGl7DXOOm3HrXRJQDAunBEEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADCYAze6AAAYxfN++rEbXcJe41mvOWOjSwAY2qxHBKvqmKq6pKourapnrjL+FlX12mn831fVlql9S1V9parOn35eMmedAAAAI5ntiGBVbUryoiQ/lOTyJOdW1bbuvnhpsqck+Vx3f2dVHZfkBUkeP437aHcfOVd9AAAAo5rziOBRSS7t7su6+5okpyU5dsU0xyZ55fT4jCT/sapqxpoAAACGN2cQvHOSTy0NXz61rTpNd1+b5PNJvnUad1hVvb+qzq6qH5ixTgAAgKHsrTeL+XSSQ7v7s1X1gCRvrKp7dfcXlieqqhOTnJgkhx566AaUCQAAsO+Z84jgFUkOWRq+y9S26jRVdWCS2yX5bHdf3d2fTZLuPi/JR5PcY+UTdPfLuntrd2/dvHnzDC8BAABg/zNnEDw3yeFVdVhVHZTkuCTbVkyzLckJ0+PHJnl7d3dVbZ5uNpOquluSw5NcNmOtAAAAw5jt1NDuvraqTkry1iSbkpza3RdV1SlJtnf3tiR/mOTVVXVpkquyCItJ8tAkp1TVvyf5WpKndvdVc9UKAAAwklmvEezuM5OcuaLt2UuPv5rkcavM9/okr5+zNgAAgFHN+oXyAAAA7H0EQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBHLjRBQAA3BQfet7bN7qEvcZ3P+sRG10CsI9xRBAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYA7c6AIAANh4J5988kaXsNfQF4zAEUEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGc+BGFwAAAPuT151+1EaXsNf4qce9d6NLYCcEQQAAYK913zPeutEl7DU+8NgfXrdlzXpqaFUdU1WXVNWlVfXMVcbfoqpeO43/+6rasjTuf0ztl1TV+r1iAACAwc0WBKtqU5IXJfmRJEckOb6qjlgx2VOSfK67vzPJ7yR5wTTvEUmOS3KvJMckefG0PAAAAG6mOY8IHpXk0u6+rLuvSXJakmNXTHNskldOj89I8h+rqqb207r76u7+WJJLp+UBAABwM80ZBO+c5FNLw5dPbatO093XJvl8km9d47wAAADcBNXd8yy46rFJjunu/zwNPynJg7r7pKVpPjhNc/k0/NEkD0pycpL3dPdrpvY/TPKX3X3Giuc4McmJ0+A9k1wyy4tZXwcn+cxGF7Ef0Z/rS3+uH325vvTn+tKf60t/rh99ub705/raF/rzrt29eS0TznnX0CuSHLI0fJepbbVpLq+qA5PcLsln1zhvuvtlSV62jjXPrqq2d/fWja5jf6E/15f+XD/6cn3pz/WlP9eX/lw/+nJ96c/1tb/155ynhp6b5PCqOqyqDsri5i/bVkyzLckJ0+PHJnl7Lw5Rbkty3HRX0cOSHJ7El5AAAACsg9mOCHb3tVV1UpK3JtmU5NTuvqiqTkmyvbu3JfnDJK+uqkuTXJVFWMw03euSXJzk2iRP6+7r5qoVAABgJLN+oXx3n5nkzBVtz156/NUkj9vJvM9L8rw569sg+9SprPsA/bm+9Of60ZfrS3+uL/25vvTn+tGX60t/rq/9qj9nu1kMAAAAe6c5rxEEAABgLyQIrqKqnlVVF1XVBVV1flU9qKo+XlUH78YyHl5Vb76Jz//kqnrhTZl3f1BV1039/sGqOr2qbr2T6c6sqm/Z0/VtlKr6jqo6rao+WlXnTa//HjM918Or6vPT3+FDVfWcnUx3p6o6Y7VxI6iqx1RVV9V3bXQtG+XG+qCqXjF9ndB6PNeTq+pOS8Mvr6ojdjH9KVV19Ho8996kqr600TXsK9a6P7kZy7/R/fW0Pf3+peGnVtXPrGcde9JSn+74eeYq09zk90C7eN79qh/XW1Vtmb6Wbbnt5Kr6lar63qr6+6V9+skbVOZeYVd9tcq067YP2xvNeo3gvqiqvi/JjyW5f3dfPYW/gza4rNF8pbuPTJKq+uMkT03y2ztGVlVlcVrz/7VB9e1x02t+Q5JXdvdxU9t9k3x7kn+Y6WnP6e4fq6pvTnJ+Vb2pu9+3VNOB3f2PWdzxd1THJ3nX9HvVsDyAPdIHVbUpyZOTfDDJPybJju+p3Znla9IZ1i73J3vIw5N8KcnfJkl3v2QPP/96+3qf7mEPz/7Vj3vSK5P8VHd/YNqW3nOjC9oXTF9tt19zRPCG7pjkM919dZJ092emN7tJ8vSqel9VXbjj0++q+uaqOrWq3ltV76+qY5cXVlUHTEcTv2Wp7SNV9e1VtbmqXl9V504/D95VYVV1/PTcH6yqFyy1HzPV9YGqett6dcRe4pwk3zl9enNJVb0qizeChywfpa2qn5mO4H6gql49te1W/+7lfjDJvy/v+Lr7A0neVVW/Ma0TF1bV45Ovf3J6dlX9eVVdVlXPr6onTuvphVV192m6V1TVS6pqe1X9Q1X92Mon7u5/S3JeFn+HJ1fVtqp6e5K3LX+qVlWbquo3p1ouqKqnT+0PmGo5r6reWlV3nL239oCquk2ShyR5SqY7Hk//7y+uqg9X1Vm1OGr72GncftcPO+mDqqoXTv+vf53k26b2Y6rq9KV5v37EoKoeWVV/N23HTp+Wm+l//AVV9b4sgubWJH9ci0+1b1VV76yqrdO694ql/4Nfnub/+ie507Keu8o2fPP0t7qoFkcYP1G7cfbHRpr68J1Vdca0zv1xVdU07oFV9bfTNvG9VXXbqrplVf3R9PrfX1U/OE375Kp649QPH6+qk6rqGdM076mqO0zT3b2q3jKtw+fUvnck/Jwk35kk0+v74PTzS1PblqV+/NDUr7eexi3vb7ZW1TtXLryqHlWLoy7vr6q/rsV+fksW4fOXp/X2B2rpyENVHTn18QVV9Yaquv3U/s5p3X9vLbbNPzB/99w80//4h6f/159Yav+GIy1Tn2+ZHq+27x66H2fybUk+nSTdfV13X7zB9ey1pnXmf1XV9iS/ODUfXSveJ03bi3Omfcr7ajpavavt8t5IELyhv8oiZPxDLd7QPWxp3Ge6+/5Jfj/Jjo3as7L4/sOjsniz/hu1OIKSJOnuryX58yQ/niRV9aAkn+juf07yu0l+p7sfmOQnk7x8Z0XV4nSoFyR5RJIjkzywFqdkbU7yB0l+srvvm53chXVfVItPYn4kyYVT0+FJXtzd9+ruTyxNd68kv5rkEVMf7PjHXXP/7gPunUUYW+knslgf7pvk6CzWvx0B475Z7Di/O8mTktxjWk9fnuTpS8vYkuSoJD+a5CVVdcvlJ6iqb03yvUkumprun+Sx3b38v5EkJ07LOrK7vyeLN+zflOT3pukfkOTU7D93Az42yVu6+x+SfLaqHpDF32NLkiOy6PPvS5L9uB9W64Mfz+LT5iOS/EySHady/XWSBy1tHx+f5LTpzfWvJjl62r5uT/KMpef4bHffv7tfM417Yncf2d1fWZrmyCR37u57d/d9kvzRTupdbRv+nCy24fdKckaSQ29aV2yY+yX5pSz6+25JHlyL7+59bZJfnLaJRyf5SpKnJempj45P8sql//d7Z7H+PjCLdfPL3X2/JH+Xxd8xWdwt7+nTOvwrSV68B17fuljen0zr6X9K8qAstm0/W1X3mya9Zxb7me9O8oUk/2U3nuZdSb536rfTkvy37v54kpdksS86srvPWTHPq5L892mbeWG+8aj6gdM2+5eyd51xcKv6xlNDHz+tR3+Q5FFJHpDkO25sIbvYd4/Sj3vS7yS5ZArJP7dyP88NHNTdW7v7t6bhLbnh+6R/SfJD0z7l8Un+99L8N9gu76nCd9d+f8hzd3X3l6adxA9kEexeW9ef//5n0+/zcv2nXY9M8uilT7tumRu+kXhtkmdn8ebkuGk4Weycj1j6oOA/1PRJ+CoemOSd3X1l8vVTXB6a5Lokf9PdH5vqv2r3XvFe6VZVdf70+Jwsvm/yTlkE6PesMv0jkpze3Z9JvqEPVu3f7t6frq95SJI/nb5n85+r6uws1pUvJDm3uz+dJFX10Sw+5EgWO8kfXFrG66YPLD5SVZcl2fF78K6hAAAJbklEQVQp/w9U1fuTfC3J86fv93xgkrN2sp4dneQl3X1tsvg7VNW9s3iDedb0d9iU6VPJ/cDxWXzYkCzerByfxTb19Kk//6mq3jGNv2f2z37YWR/sWCf/sRZHj3d8t+xbkjyqFteV/miS/5bkYVnsLN899c1BWYSPHV6bG3dZkrtV1e8l+Ytcv66vtNo2/CGZPqjr7rdU1efW8Hx7k/d29+VJMm03tyT5fJJPd/e5SdLdX5jGPySLDyTS3R+uqk8k2XGd8Tu6+4tJvlhVn0/ypqn9wiTfM+2bvj/J6Uvb1FvM/NrWw2r7k59P8obpbIdU1Z9lsc/fluRT3f3uafrXJPmFJL+5xue6SxbvGe6YxXr8sV1NXFW3S/It3X321PTKJKcvTbK8vm5ZYw17wg1ODa2qI5N8rLs/Mg2/JosPB3dlZ/vuUfpxve3sawC6u0+Z3jc+MskTsthWP3xPFbYX2mlfTb9X7ndWe5/0sSQvnNb963L9tjRZfbv8rnWqfV0JgquY3sC8M8k7q+rCJCdMo66efl+X6/uusjgad8nyMqrq25cG/y6L0+o2J3lMkv85tR+QxadeX10x7zq9kn3WajuZJPm33VzOqv27j7oou38t3tVLj7+2NPy1fOP//soN4o7hc7r7BqeKZvf+DpXkou7+vt2YZ69Xi1PlHpHkPlXVWQS7zuI6zlVnyX7WDzehD5JFWDwpyVVJtnf3F6dTZs7q7uN3Ms+Nrm/d/blaXDP7w1kcBf+pJP/3KpOutg3f1y3/n9+c13Vj24sDkvzrBl0bdnPsbH+yMzvbHl6b68+i2tnRlN9L8tvdva2qHp7k5N2q9Ib2l/V1ue+SnfffDvrxpvlsktuvaLtDpiDd3R9N8vtV9QdJrqyqb+3uz+7hGvcWu+yr3HC/s9p24ZeT/HMWZ18dkGT5veZ6bZdn59TQFarqnlV1+FLTkUk+sbPpk7w1i2sHd1yXcb+VE/TiyxrfkMUF6h9a+sf7qyydojd9qrAz703ysKo6uBYX+h6f5Owk70ny0Ko6bFrGHW7kJe6P3p7kcdMpjMt9sDv9u7d7e5JbVNXXP2Gtqu9J8q9JHl+La6Q2Z3GU+L27uezH1eLatrtncQrDJTc2w06cleTnplOwdvwdLkmyuRY3YUpVfdN0OtC+7rFJXt3dd+3uLd19SBY7kKuS/OTUn9+e6z9x3R/7YWd98Nlcv07eMd949PnsLE4t/tksQmGy2IY9uKp2XLv1zbXzu+F+McltVzZOp5ce0N2vz+JUs/vvxut4dxbBMVX1yNzwzcG+6JIkd5yO4KcW1wcemMURsSdObffI4uyVNf2/T0cVP1ZVj5vmryl874vOSfKYqrr1dKryj09tSXLojv/TLI6c7PgU/+NZnPKYLC41WM3tklwxPT5hqX3V9ba7P5/kc3X9dWtPyuJ/ZF/04SRbpv1IsniPssPHM/1PVtX9kxw2te9s3z1yP95k09lOn66qRyRf789jsriXwI/ueJ+axWU212Xx/mFIu+qrncyy2vuk22Vx5sXXsljnNs1f+foTBG/oNllcN3FxVV2QxSlLJ+9i+l9L8k1JLqiqi6bh1bw2yU/nGw83/0KSrbW4uPniLD7J3uHJVXX5jp8sVrBnJnlHkg8kOa+7/3w6VfTEJH9WVR/I2k6j2q9090VZXNNy9tQHO+4It6v+3adMHyb8eBYXLH90Wtf+/yR/kuSCLNaJt2dxLcU/7ebiP5lFePzLJE+9GUdQXz4t64Lp7/CE7r4mi8Dwgqnt/Fx/zdi+7Pjc8MjX67O4LubyJBdncVrZ+5J8fj/th531wR2TfCSLPnhVlk7znM62eHMW12q9eWq7Mou7gf7ptM39u1x/evJKr8ji+ozzq+pWS+13zuIMjvOz6Pf/sRuv47lJHlmLmx49Lsk/ZfGGc581rW+PT/J70/p2VhZHYV6c5IDpTJfXJnlyTzdGW6MnJnnKtMyLsrhGdJ/Ti7sfvyKL7d7fJ3l5d79/Gn1JkqdV1Yey+FDg96f25yb53VrcQOK6nSz65CxOnT0vyWeW2t+U5Men9XblzUpOyOLa7guy+OD5lJvz2vaQldcIPn/ab5yY5C9qcbOYf1ma/vVJ7jDtt07KdKfrXey7T84Y/TiHn0ny/03bwrcnee50JPBJWVwjeH6SV2dxrfXO1uNR7KyvVrPa+6QXJzlhWne/K7t/1tpeoRbvL4ERVdUrkry5u4f9LsD1VtN1qNOn3O9N8uCbEM7ZQ6rqFkmum65h/L4kv78Pnv7IOqjFnSnf3N333uBSAPaIvfacVYB91Jtr8XUxByX5NSFwr3doktdV1QFJrsnitFUA2O85IggAADAY1wgCAAAMRhAEAAAYjCAIAAAwGEEQAG6mqvrS9PtOVbXLu/BW1S9V1a33TGUAsDo3iwGAVVTVprV+11ZVfam7b7PGaT+eZGt3f+bGpr0ptQDAWjgiCMBwqmpLVX24qv64qj5UVWdU1a2r6uNV9YLpS7EfV1V3r6q3VNV5VXVOVX3XNP9hVfV3VXVhVf3PFcv94PR4U1X9ZlV9sKouqKqnV9UvJLlTkndU1Tum6Y6flvPBqnrB0rK+VFW/NX1h8fftyf4BYP/newQBGNU9kzylu99dVacm+S9T+2e7+/5JUlVvS/LU7v5IVT0oyYuTPCLJ72bx5fOvqqqn7WT5JybZkuTI6Qvr79DdV1XVM5L8YHd/pqrulOQFSR6Q5HNJ/qqqHtPdb0zyzUn+vrv/6yyvHoChOSIIwKg+1d3vnh6/JslDpsevTZKquk2S709yelWdn+SlSe44TfPgJH86PX71TpZ/dJKXdve1SdLdV60yzQOTvLO7r5ym++MkD53GXZfk9TflhQHAjXFEEIBRrbxIfsfwv02/D0jyr9195BrnX29fdV0gAHNxRBCAUR1aVTuuvXtCknctj+zuLyT5WFU9Lklq4b7T6HcnOW56/MSdLP+sJD9XVQdO899hav9ikttOj9+b5GFVdXBVbUpyfJKzb97LAoAbJwgCMKpLkjytqj6U5PZJfn+VaZ6Y5CnTDVsuSnLs1P6L07wXJrnzTpb/8iSfTHLBNP8TpvaXJXlLVb2juz+d5JlJ3pHkA0nO6+4/v/kvDQB2zddHADCcqtqS5M3dfe8NLgUANoQjggAAAINxRBAAAGAwjggCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwfwfmXi8fCM53dAAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"reg = ensemble.RandomForestRegressor(n_estimators=100, max_features=None)\n",
"train_x = train.drop('Sales', axis=1)\n",
"train_y = train.Sales\n",
"reg.fit(train_x,train_y)\n",
"test_x = test.drop('Sales', axis=1)\n",
"test_y = test.Sales\n",
"preds = reg.predict(test_x)\n",
"rmse = np.sqrt(metrics.mean_squared_error(test_y, preds))\n",
"print('training rmse: {}'.format(rmse))\n",
"\n",
"_,_ = plt.subplots(figsize=(15, 7.5))\n",
"bar_df = pd.DataFrame({'predictor': train_x.columns, 'importance' : reg.feature_importances_})\n",
"sns.barplot(x='predictor', y='importance', data=bar_df.sort_values('importance', ascending=False))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) Use random forests to analyze this data. What test MSE do you obtain? Use the importance() function to determine which variables are most important. Describe the effect of m, the number of variables considered at each split, on the error rate obtained."
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"training rmse: 1.7675543885549883\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4IAAAHICAYAAAAbVVRCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm4ZXdZJ/rvm8QwCCKYqMwVEJAIEqAIKghIR4ytEFQiiYihmzZiExxob1/64oUQ2+eCQ3ttZsQogwokOBQYwcgQA4pJBUImiCRhSkQNBJkJJnn7j7WK7JycUzmVnF2nqn6fz/Ocp/b6rWG/+1drr7W/ew27ujsAAACMY7/NLgAAAIDdSxAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADOaAzS5goxx00EG9ZcuWzS4DAABgU5xzzjmf7u6D1zPtPhMEt2zZku3bt292GQAAAJuiqj6+3mmdGgoAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDAHbHYBu9ND/6/XbnYJe4xzfvNnNrsEAABgkzgiCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMJilBsGqOrKqLq6qS6rqOauMf3ZVXVRV51XVO6rqngvjrq2qc+e/bcusEwAAYCQHLGvBVbV/kpcm+cEklyc5u6q2dfdFC5N9IMnW7v5yVf18kt9I8uR53Fe6+7Bl1QcAADCqZR4RPDzJJd19WXd/Lckbkhy1OEF3v6u7vzwPvi/J3ZZYDwAAAFluELxrkk8uDF8+t63l6Un+amH41lW1vareV1VPXG2Gqjp+nmb7lVdeecsrBgAAGMDSTg3dFVX100m2Jnn0QvM9u/uKqrpXkndW1fndfenifN39qiSvSpKtW7f2bisYAABgL7bMI4JXJLn7wvDd5rYbqKojkjw3yRO6++od7d19xfzvZUneneTBS6wVAABgGMsMgmcnuU9VHVJVByY5JskN7v5ZVQ9O8spMIfBfF9rvWFW3mh8flOQRSRZvMgMAAMDNtLRTQ7v7mqo6Icnbk+yf5OTuvrCqTkqyvbu3JfnNJLdLckpVJcknuvsJSe6f5JVVdV2msPrCFXcbBQAA4GZa6jWC3X1aktNWtD1v4fERa8z3d0keuMzaAAAARrXUH5QHAABgzyMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYA7Y7ALYe33ipAdudgl7hHs87/zNLgEAAHaJI4IAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEsNQhW1ZFVdXFVXVJVz1ll/LOr6qKqOq+q3lFV91wYd1xVfWT+O26ZdQIAAIxkaUGwqvZP8tIkP5zk0CTHVtWhKyb7QJKt3f3dSU5N8hvzvHdK8vwkD09yeJLnV9Udl1UrAADASJZ5RPDwJJd092Xd/bUkb0hy1OIE3f2u7v7yPPi+JHebH/9QktO7+6ru/myS05McucRaAQAAhrHMIHjXJJ9cGL58blvL05P81a7MW1XHV9X2qtp+5ZVX3sJyAQAAxrBH3Cymqn46ydYkv7kr83X3q7p7a3dvPfjgg5dTHAAAwD5mmUHwiiR3Xxi+29x2A1V1RJLnJnlCd1+9K/MCAACw65YZBM9Ocp+qOqSqDkxyTJJtixNU1YOTvDJTCPzXhVFvT/K4qrrjfJOYx81tAAAA3EIHLGvB3X1NVZ2QKcDtn+Tk7r6wqk5Ksr27t2U6FfR2SU6pqiT5RHc/obuvqqpfyxQmk+Sk7r5qWbUCAACMZGlBMEm6+7Qkp61oe97C4yN2Mu/JSU5eXnUAAABj2iNuFgMAAMDuIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYzLqDYFXds6qOmB/fpqpuv7yyAAAAWJZ1BcGq+tkkpyZ55dx0tyR/vqyiAAAAWJ71HhF8ZpJHJPl8knT3R5J867KKAgAAYHnWGwSv7u6v7RioqgOS9HJKAgAAYJnWGwTPqKr/J8ltquoHk5yS5C3LKwsAAIBlWW8QfE6SK5Ocn+TnkpyW5FeXVRQAAADLc8A6p7tNkpO7+/eSpKr2n9u+vKzCAAAAWI71HhF8R6bgt8NtkvzNxpcDAADAsq03CN66u7+4Y2B+fNvllAQAAMAyrTcIfqmqHrJjoKoemuQryykJAACAZVrvNYK/lOSUqvqnJJXk25M8eWlVAQAAsDTrOiLY3Wcn+c4kP5/kGUnu393n3NR8VXVkVV1cVZdU1XNWGf+oqnp/VV1TVU9aMe7aqjp3/tu2vpcDAADATVnvEcEkeViSLfM8D6mqdPdr15p4vrPoS5P8YJLLk5xdVdu6+6KFyT6R5GlJfmWVRXyluw/bhfoAAABYh3UFwap6XZJ7Jzk3ybVzcydZMwgmOTzJJd192byMNyQ5KsnXg2B3f2wed92uFg4AAMDNs94jgluTHNrdvQvLvmuSTy4MX57k4bsw/62ranuSa5K8sLv/fBfmBQAAYA3rDYIXZLpBzKeWWMtK9+zuK6rqXkneWVXnd/elixNU1fFJjk+Se9zjHruxNAAAgL3XeoPgQUkuqqqzkly9o7G7n7CTea5IcveF4bvNbevS3VfM/15WVe9O8uAkl66Y5lVJXpUkW7du3ZWjlQAAAMNabxA88WYs++wk96mqQzIFwGOS/NR6ZqyqOyb5cndfXVUHJXlEkt+4GTUAAACwwrqCYHefsasL7u5rquqEJG9Psn+Sk7v7wqo6Kcn27t5WVQ9L8mdJ7pjk8VX1gu7+riT3T/LK+SYy+2W6RvCiNZ4KAACAXbDeu4Z+T5IXZwpoB2YKdl/q7m/a2XzdfVqS01a0PW/h8dmZThldOd/fJXngemoDAABg16zrB+WTvCTJsUk+kuQ2Sf5Lpt8IBAAAYC+z3iCY7r4kyf7dfW13/0GSI5dXFgAAAMuy3pvFfLmqDkxyblX9RqafkVh3iAQAAGDPsd4w99R52hOSfCnTz0L8+LKKAgAAYHnWGwSf2N1f7e7Pd/cLuvvZSX50mYUBAACwHOsNgset0va0DawDAACA3WSn1whW1bGZfgT+XlW1bWHU7ZNctczCAAAAWI6bulnM32W6McxBSX57of0LSc5bVlEAAAAsz06DYHd/vKouT/LV7j5jN9UEAADAEt3kNYLdfW2S66rqDruhHgAAAJZsvb8j+MUk51fV6Zl+PiJJ0t2/sJSqAAAAWJr1BsE/nf8AAADYy60rCHb3a6rqwCT3nZsu7u5/X15ZAAAALMu6gmBVPSbJa5J8LEkluXtVHdfdf7u80gAAAFiG9Z4a+ttJHtfdFydJVd03yZ8keeiyCgMAAGA5bvKuobNv2BECk6S7/zHJNyynJAAAAJZpvUcEt1fVq5O8fh5+SpLtyykJAACAZVpvEPz5JM9MsuPnIs5M8rKlVAQAAMBSrfeuoVdX1UuSvCPJdZnuGvq1pVYGAADAUqz3rqE/kuQVSS7NdNfQQ6rq57r7r5ZZHAAAABtvV+4a+gPdfUmSVNW9k/xlEkEQAABgL7Peu4Z+YUcInF2W5AtLqAcAAIAl25W7hp6W5E1JOsnRSc6uqh9Pku7+0yXVBwAAwAZbbxC8dZJ/SfLoefjKJLdJ8vhMwVAQBAAA2Eus966h/2nZhQAAALB7rPeuoYckeVaSLYvzdPcTllMWAAAAy7LeU0P/PMnvJ3lLpt8RBAAAYC+13iD41e7+30utBAAAgN1ivUHwd6vq+Un+OsnVOxq7+/1LqQoAAIClWW8QfGCSpyZ5bK4/NbTnYQAAAPYi6w2CRye5V3d/bZnFAAAAsHz7rXO6C5J88zILAQAAYPdY7xHBb07y4ao6Oze8RtDPRwAAAOxl1hsEn7/UKgAAANht1hUEu/uMZRcCAADA7rHTIFhV7+nuR1bVFzLdJfTro5J0d3/TUqsDAABgw+00CHb3I+d/b797ygEAAGDZ1nvXUAAAAPYRgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwB2x2AUDyiBc/YrNL2GO891nv3ewSAAD2eUs9IlhVR1bVxVV1SVU9Z5Xxj6qq91fVNVX1pBXjjquqj8x/xy2zTgAAgJEsLQhW1f5JXprkh5McmuTYqjp0xWSfSPK0JH+8Yt47JXl+kocnOTzJ86vqjsuqFQAAYCTLPCJ4eJJLuvuy7v5akjckOWpxgu7+WHefl+S6FfP+UJLTu/uq7v5sktOTHLnEWgEAAIaxzCB41ySfXBi+fG7bsHmr6viq2l5V26+88sqbXSgAAMBI9uq7hnb3q7p7a3dvPfjggze7HAAAgL3CMoPgFUnuvjB8t7lt2fMCAACwE8sMgmcnuU9VHVJVByY5Jsm2dc779iSPq6o7zjeJedzcBgAAwC20tCDY3dckOSFTgPtQkjd194VVdVJVPSFJquphVXV5kqOTvLKqLpznvSrJr2UKk2cnOWluAwAA4BZa6g/Kd/dpSU5b0fa8hcdnZzrtc7V5T05y8jLrAwAAGNFefbMYAAAAdp0gCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwB2x2AQAb7YxHPXqzS9hjPPpvz9jsEgCAPZAjggAAAIMRBAEAAAbj1FAAduol/+0tm13CHuOE3378ZpcAABvCEUEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABnPAZhcAAKP49Z9+0maXsMd47utP3ewSAIbmiCAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGL8jCADslT706+/c7BL2GPd/7mM3uwRgL+OIIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAaz1CBYVUdW1cVVdUlVPWeV8beqqjfO4/+hqrbM7Vuq6itVde7894pl1gkAADCSA5a14KraP8lLk/xgksuTnF1V27r7ooXJnp7ks939HVV1TJIXJXnyPO7S7j5sWfUBAACMaplHBA9Pckl3X9bdX0vyhiRHrZjmqCSvmR+fmuQ/VFUtsSYAAIDhLTMI3jXJJxeGL5/bVp2mu69J8rkk3zKPO6SqPlBVZ1TV96/2BFV1fFVtr6rtV1555cZWDwAAsI/aU28W86kk9+juByd5dpI/rqpvWjlRd7+qu7d299aDDz54txcJAACwN1pmELwiyd0Xhu82t606TVUdkOQOST7T3Vd392eSpLvPSXJpkvsusVYAAIBhLDMInp3kPlV1SFUdmOSYJNtWTLMtyXHz4ycleWd3d1UdPN9sJlV1ryT3SXLZEmsFAAAYxtLuGtrd11TVCUnenmT/JCd394VVdVKS7d29LcnvJ3ldVV2S5KpMYTFJHpXkpKr69yTXJXlGd1+1rFoBAABGsrQgmCTdfVqS01a0PW/h8VeTHL3KfG9O8uZl1gYAADCqPfVmMQAAACyJIAgAADCYpZ4aCgDA3uHEE0/c7BL2GPqCETgiCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYARBAACAwQiCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDCCIAAAwGAEQQAAgMEIggAAAIMRBAEAAAYjCAIAAAxGEAQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwmAM2uwAAANiXvOmUwze7hD3GTx591maXwBocEQQAABiMIAgAADAYQRAAAGAwgiAAAMBgBEEAAIDBCIIAAACD8fMRAADAHutBp759s0vYY3zwST+0YctyRBAAAGAwgiAAAMBgBEEAAIDBCIIAAACDEQQBAAAGIwgCAAAMRhAEAAAYjCAIAAAwGEEQAABgMIIgAADAYJYaBKvqyKq6uKouqarnrDL+VlX1xnn8P1TVloVx/2Nuv7iqfmiZdQIAAIxkaUGwqvZP8tIkP5zk0CTHVtWhKyZ7epLPdvd3JPmdJC+a5z00yTFJvivJkUleNi8PAACAW2iZRwQPT3JJd1/W3V9L8oYkR62Y5qgkr5kfn5rkP1RVze1v6O6ru/ujSS6ZlwcAAMAttMwgeNckn1wYvnxuW3Wa7r4myeeSfMs65wUAAOBmqO5ezoKrnpTkyO7+L/PwU5M8vLtPWJjmgnmay+fhS5M8PMmJSd7X3a+f238/yV9196krnuP4JMfPg/dLcvFSXszGOijJpze7iH2I/txY+nPj6MuNpT83lv7cWPpz4+jLjaU/N9be0J/37O6D1zPhAUss4ookd18Yvtvctto0l1fVAUnukOQz65w33f2qJK/awJqXrqq2d/fWza5jX6E/N5b+3Dj6cmPpz42lPzeW/tw4+nJj6c+Nta/15zJPDT07yX2q6pCqOjDTzV+2rZhmW5Lj5sdPSvLOng5RbktyzHxX0UOS3CfJWUusFQAAYBhLOyLY3ddU1QlJ3p5k/yQnd/eFVXVSku3dvS3J7yd5XVVdkuSqTGEx83RvSnJRkmuSPLO7r11WrQAAACNZ5qmh6e7Tkpy2ou15C4+/muToNeb99SS/vsz6NsledSrrXkB/biz9uXH05cbSnxtLf24s/blx9OXG0p8ba5/qz6XdLAYAAIA90zKvEQQAAGAPJAiuoqqeW1UXVtV5VXVuVT28qj5WVQftwjIeU1VvvZnP/7SqesnNmXdfUFXXzv1+QVWdUlW3XWO606rqm3d3fXuzqnpiVXVVfedm17JZbqoPquoP55+/2YjnelpV3WVh+NVVdehOpj+pqo7YiOfenarq26vqDVV1aVWdM78377uk53pMVX1u3kZ8qKqev8Z0d6mqU1cbt7erqi9udg17i/XuT27B8m9yfz2vs9+3MPyMqvqZjaxjd1ro0x1/z1llmpv9GWgnz7tP9eNGq6ot88+yLbadWFW/UlXfU1X/sLDdPHGTytwj7KyvVpl2wz4T7ImWeo3g3qiqvjfJjyZ5SHdfPYe/Aze5rNF8pbsPS5Kq+qMkz0jyv3aMrKrKdFrzf9yk+vZmxyZ5z/zvqh+gB7Bb+qCq9k/ytCQXJPmnJNnxu6prWbyGem8xvx//LMlruvuYue1BSb4tyT8u6WnP7O4frapvTHJuVb2lu9+/UNMB3f1Pme5Gzdh2uj/ZTR6T5ItJ/i5JuvsVu/n5N9rX+3Q3e0z2rX7cnV6T5Ce7+4Pzvul+m13Q3mD+abt9miOCN3bnJJ/u7quTpLs/PX+gSJJnVdX7q+r8HUcTquobq+rkqjqrqj5QVUctLqyq9puPJn7zQttHqurbqurgqnpzVZ09/z1iZ4VV1bHzc19QVS9aaD9yruuDVfWOjeqIPcSZSb5j/vbm4qp6baYP1ndfPEpbVT8zH8H9YFW9bm7bpf7d11XV7ZI8MsnTM9+hd14/X1ZVH66q0+cjOU+axz20qs6Yj/C8varuvInlb4g1+qCq6iXz+vU3Sb51bj+yqk5ZmPfr33BX1eOq6u/n990p83Izr5Mvqqr3ZwqaW5P80fwt7G2q6t1VtbWq9p+/Zbxgfk//8jz/1795nJf1glW2OQfP/1cX1nSE8eO1C2crLMEPJPn3xQ9l3f3BJO+pqt9ceI1PTr7ej2dU1V9U1WVV9cKqesq8DT2/qu49T/eHVfWKqtpeVf9YVT+68om7+0tJzsm0jXhaVW2rqncmeUctfOM79/dvzbWcV1XPmtv36nV87st3V9Wp83v4j6qq5nEPq6q/m7eJZ1XV7avq1lX1B3M/f6CqfmCe9mlV9efzevWxqjqhqp49T/O+qrrTPN29q+ptc3+dWXvfmQVnJvmOJJlf3wXz3y/NbVsW+vFDc7/edh63uL/ZWlXvXrnwqnp8TUddPlBVf1PTfn5LpvD5y/N24Ptr4chDVR029/F5VfVnVXXHuf3d87bkrHn9//7ld88tU9M288Pz9u/HF9pvcKRl7vMt8+PV9t1D9+OSfGuSTyVJd1/b3Rdtcj17rHmd+f+ranuSX5ybj1i5L5q3F2fWtI9+f81Hq3e2Xd4TCYI39teZQsY/1vQB+dEL4z7d3Q9J8vIkOzZqz830+4eHZ/pA9Js1fUudJOnu65L8RZIfS5KqeniSj3f3vyT53SS/090PS/ITSV69VlE1nV72oiSPTXJYkofVdIrbwUl+L8lPdPeDssZdWPdGNX0T88NJzp+b7pPkZd39Xd398YXpvivJryZ57NwHO9646+7fQRyV5G3d/Y9JPlNVD820s96S5NAkT03yvUlSVd+Q5MVJntTdD01ycvaNu/iu1gc/lunb0UOT/EySHace/U2Shy+8n5+c5A3zh8FfTXLEvD3YnuTZC8/xme5+SHe/fh73lO4+rLu/sjDNYUnu2t0P6O4HJvmDNepdbZvz/EzbnO9KcmqSe9y8rtgwD8gUxlb68Uyv80FJjsi0bdwRtB6U6UPd/TOtd/edt6GvTvKshWVsSXJ4kh9J8oqquvXiE1TVtyT5niQXzk0PybTOLm63k+T4eVmHdfd3Zwrn+8o6/uAkv5Rp/b1XkkfU9Nu9b0zyi/M28YgkX0nyzCQ9r3PHJnnNQp8+INP/2cMy9cOXu/vBSf4+0/sime6W96y5v34lyct2w+vbEIv7k/l9/5+SPDzT+vOzVfXgedL7ZdrP3D/J55P81114mvck+Z65396Q5L9398eSvCLTvuiw7j5zxTyvTfJ/z+vl+bnhWQoHzO+LX8qedQbHbeqGp4Y+eV6Pfi/J45M8NMm339RCdrLvHqUfd6ffSXLxHJJ/buW2lBs5sLu3dvdvz8NbcuN90b8m+cF5H/3kJP97Yf4bbZd3V+G7ap8/5LmruvuL807i+zMFuzfW9ee//+n87zm5/tuuxyV5wsK3XbfOjT+YvTHJ8zJ92DtmHk6mnfOhC18UfFPNRxZW8bAk7+7uK5Ovn+LyqCTXJvnb7v7oXP9Vu/aK90i3qapz58dnZvq9ybtkCtDvW2X6xyY5pbs/ndygD1bt3+4e9fqaYzOF42TauR6baRtwyvyFxT9X1bvm8ffL9MHw9Ln/9s/8beJebq0++JP5t0r/qaYjSjt+C/VtSR5f07VmP5Lkvyd5dKaN+3vnvjkw04flHd6Ym3ZZkntV1YuT/GWmL6BWs9o255GZv1jq7rdV1WfX8Xyb4ZG5vl//parOyLQd+3ySs7v7U0lSVZfm+td/fqbt7g5vmtfNj1TVZUl2HIH6/qr6QJLrkrxw/u3ZhyU5fY1t4BFJXtHd1yTTNqKqHpB9Yx0/q7svT5J5u7klyeeSfKq7z06S7v78PP6RmcJvuvvDVfXxJDuu5XxXd38hyReq6nNJ3jK3n5/ku+d90/clOWVhm3qrJb+2jbDa/uTnk/zZfEQ5VfWnmfb525J8srvfO0//+iS/kOS31vlcd8v0meHOmbYLH93ZxFV1hyTf3N1nzE2vSXLKwiSL7/8t66xhd7jRqaFVdViSj3b3R+bh12f6AmZn1tp3j9KPG22tnwHo7j5p/tz4uCQ/lWnf95jdVdhAPExOAAAIDUlEQVQeaM2+mv9duR9fbV/00SQvmdf9a3P9tjRZfbv8ng2qfUMJgquYP7i8O8m7q+r8JMfNo66e/7021/ddZToad/HiMqrq2xYG/z7TqUsHJ3likv85t++X6Vuvr66Yd4NeyV5rtZ1MknxpF5ezav+OqKZTux6b5IFV1Zk+9Hama7tWnSXJhd39vbupxKW7GX2QTGHxhCRXJdne3V+YT/E4vbuPXWOem1xPu/uzNV1H90OZjoz9ZJL/vMqkq21z9jQXZtevxbt64fF1C8PX5Yavc+XOesfwmd19o1NFs2vbiH1lHV/sy1uyntzU/8l+Sf5tk64NuyXW2p+sZa117ppcfxbVWkdTXpzkf3X3tqp6TJITd6nSG9sb3v/rsdh3ydr9t4N+vHk+k+SOK9rulDlId/elSV5eVb+X5Mqq+pbu/sxurnFPsdO+yo33JattF345yb9kOsNlvySLnzU3aru8dE4NXaGq7ldV91loOizJx9eaPsnbM107uOO6jAevnKCnH2v8s0wXqH9o4Y3311k4DWr+VmEtZyV5dFUdVNOFvscmOSPJ+5I8qqoOmZdxp5t4ifuidyY5ej5NbLEPdqV/93VPSvK67r5nd2/p7rtn2uBdleQnarpW8Nty/TeEFyc5uKabJ6WqvmE+jWdvtlYffCbJk2u6juzOueERqTMynW74s5lCYTK95x5RVTuuNfrGWvsOmV9IcvuVjfPppft195sznRr1kF14He/NFBxTVY/LjXdmu9s7k9yqqr7+7X9VfXeSf8v1/XpwpjMYztrFZR89r5v3znR6zcU3NcMaTk/yc/PpgTu2EfviOr7DxUnuPB8lTU3XBx6Q6YjYU+a2+2Y6e2VdfTofVfxoVR09z1/zlxl7ozOTPLGqbjuf+v1jc1uS3GPHOpHpyMmOb/E/lumUx2S61GA1d0hyxfz4uIX2VbcD3f25JJ+t669be2qmbc7e6MNJtszv1WT6jLLDxzJv46rqIUkOmdvX2neP3I8323y206eq6rHJ1/vzyEzXa//Ijs+pmS6zuTbTNnpIO+urNWZZbV90h0xnXlyXaZ3bf/mVbzxB8MZul+m6iYuq6rxMp4CduJPpfy3JNyQ5r6ounIdX88YkP50bHm7+hSRba7q4+aJMRwZ2eFpVXb7jL9MK9pwk70rywSTndPdfzKeKHp/kT6vqg1nfaWn7lO6+MNM1LWfMfbDjjnA769/RHJsbH/l6c6brOC5PclGm06Den+Rz3f21TMHpRXOfnpvrr53bW63VB3dO8pFMffDaLJzmOZ8d8NZM1xa9dW67MtPdQP9k3kb8fa4/ZXGlP8x0PcG5VXWbhfa7Zjrj4NxM/f4/duF1vCDJ42q6EcrRSf450wekTTF/0fVjmS6mv3TeDv5/Sf44yXmZtlfvzHSdzz/v4uI/kSk8/lWSZ9yCo/uvnpd13rw+/9Q+uo4nSebX9uQkL55f2+mZjsK8LMl+85kub0zytJ5vjLZOT0ny9HmZF2a65nav09MdZv8w07r1D0le3d0fmEdfnOSZVfWhTF+yvHxuf0GS363pBhLXrrHoEzOdOntOkk8vtL8lyY/N24GVNys5LtP1s+dl+uL5pFvy2naTldcIvnB+bx6f5C9rulnMvy5M/+Ykd5q3DSdkvpvwTvbdJ2aMflyGn0ny/877lncmecF8JPCpma4RPDfJ6zJdu77WejyKtfpqNavti16W5Lh53f3O7PpZa3uEmvbhwKhqvm5y/lb2rCSPuBkf2NlNqupWSa6dr2H83iQv3wtP17tJVfWHSd7a3fvkbwGy56npzpRv7e4HbHIpALvFHnvOKrDbvLWmnzc5MMmvCYF7vHskeVNV7Zfka5lOWwUA2CWOCAIAAAzGNYIAAACDEQQBAAAGIwgCAAAMRhAEgFuoqr44/3uXqtrpnU6r6peq6ra7pzIAWJ2bxQDAKqpq//X+1lZVfbG7b7fOaT+WZGt3f/qmpr05tQDAejgiCMBwqmpLVX24qv6oqj5UVadW1W2r6mNV9aL5R7GPrqp7V9Xbquqcqjqzqr5znv+Qqvr7qjq/qv7niuVeMD/ev6p+q6ouqKrzqupZVfULSe6S5F1V9a55umPn5VxQVS9aWNYXq+q35x8s/t7d2T8A7Pv8jiAAo7pfkqd393ur6uQk/3Vu/0x3PyRJquodSZ7R3R+pqocneVmSxyb53SQv7+7XVtUz11j+8Um2JDmsu6+pqjt191VV9ewkP9Ddn66quyR5UZKHJvlskr+uqid2958n+cYk/9Dd/20prx6AoTkiCMCoPtnd750fvz7JI+fHb0ySqrpdku9LckpVnZvklUnuPE/ziCR/Mj9+3RrLPyLJK7v7miTp7qtWmeZhSd7d3VfO0/1RkkfN465N8uab88IA4KY4IgjAqFZeJL9j+Evzv/sl+bfuPmyd82+0r7ouEIBlcUQQgFHdo6p2XHv3U0nesziyuz+f5KNVdXSS1ORB8+j3JjlmfvyUNZZ/epKfq6oD5vnvNLd/Icnt58dnJXl0VR1UVfsnOTbJGbfsZQHATRMEARjVxUmeWVUfSnLHJC9fZZqnJHn6fMOWC5McNbf/4jzv+UnuusbyX53kE0nOm+f/qbn9VUneVlXv6u5PJXlOkncl+WCSc7r7L275SwOAnfPzEQAMp6q2JHlrdz9gk0sBgE3hiCAAAMBgHBEEAAAYjCOCAAAAgxEEAQAABiMIAgAADEYQBAAAGIwgCAAAMBhBEAAAYDD/B1iKMYvXtjnGAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"line_df = pd.DataFrame({'m' : rnge * car_df.shape[1] - 1, 'rmse' : sqrt_rmses})\n",
"sns.lineplot(x='m', y='rmse', data=line_df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"9. This problem involves the OJ data set which is part of the ISLR package."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Purchase
\n",
"
WeekofPurchase
\n",
"
StoreID
\n",
"
PriceCH
\n",
"
PriceMM
\n",
"
DiscCH
\n",
"
DiscMM
\n",
"
SpecialCH
\n",
"
SpecialMM
\n",
"
LoyalCH
\n",
"
SalePriceMM
\n",
"
SalePriceCH
\n",
"
PriceDiff
\n",
"
Store7
\n",
"
PctDiscMM
\n",
"
PctDiscCH
\n",
"
ListPriceDiff
\n",
"
STORE
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
1
\n",
"
237
\n",
"
1
\n",
"
1.75
\n",
"
1.99
\n",
"
0.00
\n",
"
0.0
\n",
"
0
\n",
"
0
\n",
"
0.500000
\n",
"
1.99
\n",
"
1.75
\n",
"
0.24
\n",
"
0
\n",
"
0.000000
\n",
"
0.000000
\n",
"
0.24
\n",
"
1
\n",
"
\n",
"
\n",
"
1
\n",
"
1
\n",
"
239
\n",
"
1
\n",
"
1.75
\n",
"
1.99
\n",
"
0.00
\n",
"
0.3
\n",
"
0
\n",
"
1
\n",
"
0.600000
\n",
"
1.69
\n",
"
1.75
\n",
"
-0.06
\n",
"
0
\n",
"
0.150754
\n",
"
0.000000
\n",
"
0.24
\n",
"
1
\n",
"
\n",
"
\n",
"
2
\n",
"
1
\n",
"
245
\n",
"
1
\n",
"
1.86
\n",
"
2.09
\n",
"
0.17
\n",
"
0.0
\n",
"
0
\n",
"
0
\n",
"
0.680000
\n",
"
2.09
\n",
"
1.69
\n",
"
0.40
\n",
"
0
\n",
"
0.000000
\n",
"
0.091398
\n",
"
0.23
\n",
"
1
\n",
"
\n",
"
\n",
"
3
\n",
"
0
\n",
"
227
\n",
"
1
\n",
"
1.69
\n",
"
1.69
\n",
"
0.00
\n",
"
0.0
\n",
"
0
\n",
"
0
\n",
"
0.400000
\n",
"
1.69
\n",
"
1.69
\n",
"
0.00
\n",
"
0
\n",
"
0.000000
\n",
"
0.000000
\n",
"
0.00
\n",
"
1
\n",
"
\n",
"
\n",
"
4
\n",
"
1
\n",
"
228
\n",
"
7
\n",
"
1.69
\n",
"
1.69
\n",
"
0.00
\n",
"
0.0
\n",
"
0
\n",
"
0
\n",
"
0.956535
\n",
"
1.69
\n",
"
1.69
\n",
"
0.00
\n",
"
1
\n",
"
0.000000
\n",
"
0.000000
\n",
"
0.00
\n",
"
0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Purchase WeekofPurchase StoreID PriceCH PriceMM DiscCH DiscMM \\\n",
"0 1 237 1 1.75 1.99 0.00 0.0 \n",
"1 1 239 1 1.75 1.99 0.00 0.3 \n",
"2 1 245 1 1.86 2.09 0.17 0.0 \n",
"3 0 227 1 1.69 1.69 0.00 0.0 \n",
"4 1 228 7 1.69 1.69 0.00 0.0 \n",
"\n",
" SpecialCH SpecialMM LoyalCH SalePriceMM SalePriceCH PriceDiff \\\n",
"0 0 0 0.500000 1.99 1.75 0.24 \n",
"1 0 1 0.600000 1.69 1.75 -0.06 \n",
"2 0 0 0.680000 2.09 1.69 0.40 \n",
"3 0 0 0.400000 1.69 1.69 0.00 \n",
"4 0 0 0.956535 1.69 1.69 0.00 \n",
"\n",
" Store7 PctDiscMM PctDiscCH ListPriceDiff STORE \n",
"0 0 0.000000 0.000000 0.24 1 \n",
"1 0 0.150754 0.000000 0.24 1 \n",
"2 0 0.000000 0.091398 0.23 1 \n",
"3 0 0.000000 0.000000 0.00 1 \n",
"4 1 0.000000 0.000000 0.00 0 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"oj_df = pd.read_csv('oj.csv')\n",
"oj_df = oj_df.drop(oj_df.columns[0], axis=1)\n",
"oj_df.Purchase = oj_df.Purchase.map({'CH' : 1, 'MM': 0})\n",
"oj_df.Store7 = oj_df.Store7.map({'Yes' : 1, 'No': 0})\n",
"oj_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations."
]
},
{
"cell_type": "code",
"execution_count": 133,
"metadata": {},
"outputs": [],
"source": [
"train, test = model_selection.train_test_split(oj_df, test_size=oj_df.shape[0] - 800)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) Fit a tree to the training data, with Purchase as the response and the other variables as predictors. Use the summary() function to produce summary statistics about the tree, and describe the results obtained. What is the training error rate? How many terminal nodes does the tree have?"
]
},
{
"cell_type": "code",
"execution_count": 162,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"train error: 0.16249999999999998\n"
]
}
],
"source": [
"reg = tree.DecisionTreeClassifier(max_depth=3) \n",
"train_x = train.drop('Purchase', axis=1)\n",
"train_y = train.Purchase\n",
"reg.fit(train_x,train_y)\n",
"\n",
"# train error\n",
"preds = reg.predict(train_x)\n",
"train_error = 1 - ((train_y == preds).sum() / train_y.shape[0])\n",
"print('train error: {}'.format(train_error))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) Type in the name of the tree object in order to get a detailed text output. Pick one of the terminal nodes, and interpret the information displayed. (d) Create a plot of the tree, and interpret the results."
]
},
{
"cell_type": "code",
"execution_count": 138,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 138,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dot_dat = tree.export_graphviz(reg, out_file=None, rotate=True)\n",
"graph = gv.Source(dot_dat)\n",
"graph"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) Predict the response on the test data, and produce a confusion matrix comparing the test labels to the predicted test labels. What is the test error rate?"
]
},
{
"cell_type": "code",
"execution_count": 163,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
0
\n",
"
1
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
75
\n",
"
22
\n",
"
\n",
"
\n",
"
1
\n",
"
27
\n",
"
146
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1\n",
"0 75 22\n",
"1 27 146"
]
},
"execution_count": 163,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"test_x = test.drop('Purchase', axis=1)\n",
"test_y = test.Purchase\n",
"preds = reg.predict(test_x)\n",
"conf = pd.DataFrame(metrics.confusion_matrix(test_y, preds))\n",
"conf"
]
},
{
"cell_type": "code",
"execution_count": 164,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"test error rate: 0.18148148148148147\n"
]
}
],
"source": [
"print('test error rate: {}'.format( 1 - ((conf[0][0] + conf[1][1]) / test_y.shape[0])))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(f) Apply the cv.tree() function to the training set in order to determine the optimal tree size.\n",
"\n",
".. not easily available in python"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(g) Produce a plot with tree size on the x-axis and cross-validated classification error rate on the y-axis. (h) Which tree size corresponds to the lowest cross-validated classification error rate?"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"rnge = range(1, 50)\n",
"\n",
"errors = [run_cv(oj_df,\n",
" 5,\n",
" lambda df: df.drop('Purchase', axis=1),\n",
" lambda df: df.Purchase,\n",
" lambda x, y: tree.DecisionTreeClassifier(max_depth=i) .fit(x,y),\n",
" lambda preds, true: 1 - ((preds == true).sum() / preds.shape[0])).mean()[0] for i in rnge]"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8m9WV8PHfkbzvdux4iZeskJjESUiAsBQKTWloEygtbaGkU/oyL3Shw5TpdOjbKZ3SaWdKWzplyszQUroBpcB0oZQ9hDUkJEDiJM5mJ7Hjfd9XSff9Q5Kj2LIl23osyz7fz8cfrEfPY18FWee599x7rhhjUEoppcZjC3cDlFJKzXwaLJRSSgWkwUIppVRAGiyUUkoFpMFCKaVUQBoslFJKBaTBQimlVEAaLJRSSgWkwUIppVRAUeFuQKhkZmaahQsXhrsZSikVUd55551mY0xWoPNmTbBYuHAhe/bsCXczlFIqoohIZTDn6TCUUkqpgDRYKKWUCkiDhVJKqYA0WCillApIg4VSSqmANFgopZQKSIOFUkqpgDRYqIjSP+Rk+5FGdDtgpaaXBgsVMXoHHXzul7v53C9383RpXbibo9ScosFCRYTuAQc3PbSbXSdaSEuI5uGdQS06VUqFiAYLNeN19A3xmV/s4p2qNu67YS2fv2wJu060cqyhK9xNUxaobuulqWsg3M1QI2iwUDNae+8gWx/cxYGaDu7/9LlsLsnjE+vyibHbtHcxS93623e4688Hwt0MNYIGCzVjtXQPcMPPd3Gkvov/2bqOTStzAJiXFMuHV+Xwh3dr6BlwhLmVKpQcThdHG7qoau0Nd1PUCLOm6qwKv67+IQ7X+x8aSoqNYkVuStA/q6lrgBsf3EllSy8PfnY9l551ZgXlrRuK+NPeWp7aV8sN5xdOqd1q5qhq7WXIaWjojJxhKIfTRX1nP/npCeFuiqU0WKiQuevPB/njezVjPv/c37+P5TnBBYwfPH+YypZefnnTeVy0NHPU8+uK0lmek8zDOyu5/rwCRGTS7VYzR3ljNwAtPQMMOV1E22f+4Mejb1fxnafLeOmOyyialxju5lhGg4UKmd0nW7l46Ty+cNnSM4639Axw+2N72XeqPehgse9UB5cszfQbKABEhBs3FPHNPx1g76l21hamT7n9KvwqmnoAMMbdu8xLiw9ziwJ75UgTQ07Do7uq+PqHV4S7OZaZ+WFbhYUxhq7+oaDPb+0ZpLqtj0uXZXHJsswzvraU5JEYY6estjOon9U/5KS8qTvgsNW1axeQGGPn4Z1VQbdTzWzengVAQ2d/GFsSHIfTxdsnWgF4fM8p+oecYW6RdTRYqDMYY3iprIGP/tcOzv3Oi1S1BJdoLK1uB6AkP23UczabsDw3hbK64ILFsYZunC5Dcd74wSIpNoqPrl3A06W1tPUMBvWz1cxW0dRNWkI0QETkLfbXdNA94ODTFxTS1jvEM/tn72JRDRYKAJfL8Oz+Oj5y3xv87W/2UN/Rx5DTsKOiOajr91d3IAIrF/j/gC/OTeFQXRcuV+AyHWV1HcPXBLJ1QxEDDhdPvlMdVDvVzGWMoaKpmwsXzwMio2exo6IFgK9sPIvFmYmzejq3Bos5zuky/HlvDR/6j9f4wiPv0j/k5IefWM0b/3QFGYkxvFPZFtTP2VfdweLMRJLjov0+X5yXQveAg+q2voA/61BdF4kxdgozAs8uWZGbwrqidB7ZVRlUIJrNnBH++pu6Bujqd3D+ogzsNomIYLHzeAvLc5LJSo7l0xcU8m5Ve9DDrZFGg8UcNeR08cSeU2y891Vuf2wvInDfDWt58Y7LuG5dPtF2G+cWpgcdLEqr2/0OQXl58w/eXsN4ymo7WZGbgs0W3AynrRsKOdnSy5tB9oKm4icvHePKH786owKTw+nijt/v5QM/eiWiCyyWN7nzFcvmJzM/OXbGD0MNOJzsPtnKhUvcPaHr1uUTG2Xj4V2zs3ehwWKOGXA4eWRXJZf/8BX+8clSEmLs/M/Wc3nu9ku5enUedp8P6HVF6Rxv7qE1QD6gobOfxq4BSvJTxzzn7OxkbELAuy6Xy1BW1zmhNRlXrcwlfZrqRf3xvWqONnSzJ8ggarVBh4sv/+49/vBeDSdbemnvDX5SwkxT4UluL52fxPyUOBq7ZnbPYm9VO/1DruFhs7SEGLaszuNP79VMaHJIpNBgMUf0DTr55ZsnuOyeV/jGHw+QmRTLQzet5+kvX8Kmlbl+7+LXL3RPR303wAfjvlNjJ7e94mPsLM5KCpjkrm7ro3vAETC57Ssu2s4n1xfw0qFG6jus+4A53tTNSU/C/+nSWst+T7AGHE6++Mg7PHugniuWzwcIaphvpqpo6iEpNorslFhyUmIt/X8ZCjsqWrAJXOAJFuDOofUOOvnTOOuNIpUGiznA6TJc/dM3+PZfyijMSOC3N5/PH794EVcszx53MduqBalE2yXgXfT+mg7sNgmYkPYmucczkeS2r09fUIjLGH73tnXTaF8+3AjAmoI0ntlfh8Ppsux3BdI/5OSW37zDS4ca+c5HV/IPV54FwKm24MtkfOnRd7n3xaNWNXHCKpq6WZKViIiQnRI343MWbx1vYeWCVFLjT+fpVuensmpBKr/dWRnRQ4L+WBosRGSTiBwRkXIRudPP83eISJmIlIrINhEpGvF8iohUi8hPrWznbFfb3sexxm7+8UNn8/jnL+R9y7KCWvEcF23nnLzUwD2L6g7Oyk4mPsY+7nnFeSnUtPfR3jv2sFZZXRc2gbNzkgO2z1fRvEQuXZbFY7urGLLoQ3z7kUaWzU/i1ksX09w9yM7jrZb8nkB6Bx3c/OvdvHasie9/fBWf2VBEgWcywKkgayoZY9h2qIFfvH58xgyZlDd2syQrCYDslDg6+x30Dc7MdQt9g07eq2obHoLyEhG2bijkaEM3u0/OjKHKULEsWIiIHbgfuAooBm4QkeIRp70HrDfGlABPAveMeP47wGtWtXGu8BZlW1sw9jDRWNYXpbOvup1Bh/8PYGMM+6vbKVkwdr7C63SSe+yhqLLaThZnJREXPX7g8WfrhiIaOgd4qaxhwtcG0j3g4O0TrVyxfD6XL59PYow9LENR3QMObvrlbt6qaOFHn1jNp85z18VKiYsmNT466J5FY9cA/UMuembIkEn3gIO6jn6WzD8dLIAZm7fYU9nKkNMMJ7d9bVmdR3Jc1KybRmtluY/zgXJjzHEAEXkMuAYo855gjNnuc/5OYKv3gYisA7KB54D1FrZz1qv0jLMXzpt4obN1Rek8+MYJDtZ2+C2pUd3WR1vvECUFgYOFd2iprLaTi5b4L+NxqK6TdUWTK91xxfL55KXG8fCuSq5alTupnzGWN465Szpcvnw+cdF2PliczbMH6rn7mpXERIX2nqu+o5/P/GIXvX7uqnsGHXT1O/jJ9WvZsjrvjOcKMuI51RpczsL7noiJsvHwziq2bigKa32t456ZUKd7FrGA+99iKvWW/uWpg2QkxvB3H1g29Ub6eKuihSibcN7CjFHPJcRE8fFz83lkVyXN3cVkJsWe8fzuk6088GoFZ2Un87VNy0PaLitZOQy1ADjl87jac2wsNwPPAoiIDfgR8NXxfoGI3CIie0RkT1NT0xSbO3tVtvYQbRdyUydeZ8f7wT3WFNp93pXbCwL3WrKSY8lKjh0zb9HeO0hNe9+Ektu+7DbhhvMLebO8ZfjDJ1RePtxIclzU8L/HltV5dPQN8WZ56KfrvnSogWON3awrSufCJfPO+Nq4IpuHbjpvVKAAKEhPCLpnUdnirsH0fy5exJGGrrDP7qpo8s6EcgcGb8+iYQqbIO091c6vdpzk568dD3kZjh0VLawpSCMx1v/99tYNhQw5DY/vcX8EGmPYUd7M9T97i0/8z1tsO9zIL944MWOGAIMxIwoJishW3L2HyzyHvgg8Y4ypHu9uxxjzM+BnAOvXr59d2aQQqmrppSA94YxpscGanxJHQUY871S28bfvG/38/uoOYuy2oHMMxeOU/fAen2hy29enzi/gJ9uO8ciuKr65eeSo5+S4XIbtR5q49Kys4Sqo71uWRUpcFH/ZV8vlnplIofJWRQu5qXH85Po1E7rbL8hIYNvhRlwuE3CNSlVrLzaBL7x/CY/squThnZV+75KnS3ljN1E2Ge5FZCd7hqGmkOS+98WjRNmErgEHrx5t4kPn5ISkrV39Q+yv6eCL718y5jlL5yezYXEGj+6qYkVOCv/58jHerWpnfnIs39xczNnZyWz9xS5eOtTAtWvzQ9Iuq1nZs6gBCnwe53uOnUFENgLfAK42xnhvIy4EbhORk8APgb8RkX+3sK2zWmVL76SGoLzWFaazp7LN7+yOfdXtrMhNDnooZkVuCuWNXX5zIN4ex0TWWIw0PzmOD52Tw5PvVIfsbvJgbSdNXQNccfbpoBATZWPTyhxeKGsI6V2ry2V463gLFy6ZN+Fhofz0eAYdLpq7A9+NV7b0kpcWT2p8NB8/N59n99fTEsR1Vqlo7KFwXsJwME6JjyIu2jbpGVG7T7by2tEm7rjyLDISY/jLvtDll3afbMXp8p+v8LV1QxHVbX187le7aegc4DsfXclrX7ucmy9ZxEVL5pGXGsfT+yKnlpSVwWI3sExEFolIDHA98JTvCSKyFngAd6Bo9B43xtxojCk0xizEPRT1G2PMqNlUKjBjDFWtvRQFUTpjLOsWZtDUNTBqDr/LZThQ0znu+oqRivNSGHIajjWOHooqq+0cHqqaihs3FNLRNxSyD4hthxsQgfeffeYGTFtW59E94OCVI6EbAj3a2EVrz+CoWTbBKPBsvhPMUFRlay9FnhuIGy8oZNDp4vE94auvVdHUzVJPvgIYnj5bP8lV3D964QhZybF87qJFXLUyh22HGukdDM2uijvKW4iJclc4GM+HzsnhsxcW8f2Pr2L7V9/PZzYUDU/csNmEzavzeO1Y07izA2cSy4KFMcYB3AY8DxwCHjfGHBSRu0Xkas9pPwCSgCdEZK+IPDXGj1OT1NozSPeAg8IpJAnXFfrPWxxv7qF7wDHuyu2RvENM/vIWZXWdUxqC8rpw8TyWZCXy8K7QrLnYfriRNQVpzBuRqLxw8TzmJcbwlxDOitpR7i5MF+iu1Z+CDHdOKpgkd1VLD4UZ7vfEsuxkLliUwaNvh6e+lsPp4mRLz/BMKK/s5MmttdhR3szO46188f1LiI+xs7kkj74hJ9sONQa+OJifX9HCusL0gDP2ou02vn3NSj51XqHfnvfmklyGnIbnD9aHpF1Ws3SdhTHmGWPMWcaYJcaY73qO3WWMecrz/UZjTLYxZo3n62o/P+NXxpjbrGznbFbpmTY7lZ7F2TnJJMVGjQoW45UlH8uizETiom2jyn4MOlyUN3ZNOrntyz3XvYh9p9o5UBO4FtV4mroG2FfdccYQlFeU3cZVq3J4OZR3rRUtFM1LmNQWnd5rAq216Owfoq13aLhnAe4hk1Otfbx6bPonini3Ul2SNSJYpMZNOGdhjOFHLx4lNzVueLvd8xdlMD85NiRTndt6BjlU38lFkwjmI61akErRvASeLo2MoShdwT3LefejKJpCzsJuE9YWpo2aMVNa3UF8tJ2lI+4IA/2s5TkpowoKHmvsYshpQtKzAPjYufnER9unPNf9lSPuu9GxkthbPHetL4XgrtXpMuw60TLpD6K4aDtZybEBh6GG3xM+NxAfOieHzKRYHgnD2oByn5pQvrKTY6nv7J/QSuhXjjbxTmUbt12xdPjO324TPlKSy/YjTXROcfbRrhMtGAMXLZ16sBARtpTk8WZ5c1B5pnDTYDHLeefTF0yhZwFwbmE6R+o7z5jqV1rdzsoFKROeZbUiN4Wy2s4zPgRCkdz2lRofzdWr8/jz3topfUBsP9JIdkos54zR4zlvYQbZKbEhyY8crO2gq9/BhknkK7wK0gOvtfC37iYmysb15xXw8uFGatqnt76UdyvVxVlnDpVmp8TRP+Sisz+4Xpsxhh+/eJT89Hg+sa7gjOc2l+Qx6HDx4sGpLdjcUdFCQox9Qr3p8WxenYvLwLMHZv5QlAaLWa6ytYeclLhJrYj2tX5hOi7j3hsb3OPMB2snltz2Ks5LobPfQa1Pobiy2k7iom0sygzdhvdbNxTRN+TkD5PcGGnI6eL1o81cfvb8MWcm2WzCR1bl8WoI7lq9G+lMJl/hVZCRQHX7+D2Lylb3h/PIxW43XOAetvndBHM9xhgee7uK1yc5hFXR1E12SiwpI/ZCme9ZmBfsUNSLZQ2UVndw+weWjcoRnFuYxoK0+CkPRe2oaOG8hRnDs7am6uzsZJbNTwrpbC2raLCY5aqmOG3Wa01BGiLuMgcARxu6GXC4JpTc9vJdye1VVtfB8pyJ91LGsyo/ldX5qTy8q2pSRd12n2yla8ARcB3FltW5DDpdvBCCu9Zl85OY71ljMBn56fHUtvePW+SwqqWXeYkxJI1YULYgLZ4rls/nsd2nxizvMpIxhn9/7jB3/mE/n//tO0HXpvLlWxPKV453YV4QM6JcLsO9Lx5lcWYi164dvfZXRNhcksvrx5onvQVvY1c/5Y3dIclXnNmuPHafbJ3xVXY1WMxylVOcNuuVHBfN2dnJw0nuySS3vZbnJCM+e1sYYyir7QxJcnukGzcUUd7Yza4TEy/6t/1wIzF2G5cs9V+axGtNQRr56fFTujscdLjY47ORzmQVpCfgdBnqxvngqWzpHTOHdeOGIpq7B3ihLPCwiDGGu58u44FXj3Pt2gXYRPiHJ/ZNaEaVdytVf8HCu4q7PoiexTMH6jhc38XtG5cRNcZd/5bVeThchucmOfvoLU/Pb6xSNZO1eXUuxsBfZ/j+3RosItD+6o7hN+54egcdNHUNTCm57Wv9wnTeq2rH6TKU1nSQHBfFwkn87MTYKBbOSxxOctd29NPZ7whZvsLXlpI8UuKi+O0kErcvH27kgsUZY5Z08PLeHb5Z3hxwo6ixlFa30zvonPJd63D12XGS3FWtvWPWW7psWRYFGfH84o0T4w6ruVyGb/75AL988ySfu3gh935yNXdtKebtE6089OaJoNvr3UrV3ySJ4ZIfAYKF0+XOVZyVncSWktFlULzOyUthUWbipIei3qpoISUuKuQ3NUuykjgnL2XGD0VpsIhA9zx/mK8+sS/ged5qs1NZY+FrXVE63QMOjjZ0ebZRTZ108TnfvS28PYxQzYTyFR9j52Pn5vPiwQY6+oLPKVS29FDR1MPlfqbM+vOxcxfgMobvP3t4Uu3cUdGCCFywaOo9C4DqMZLcAw4ntR19Y+5vbrMJX7hsKe9VtXPxv7/MvS8cGTVs43QZvv6H/Ty8s4pbL1vMXZuLERGuW5fPB4uzuef5IxxtGH/fEq/yEQUEfcXH2EmJiwqYszhU10lFUw+3XLpk3DIn3qGotypaaJpgzamO3iH+ur+OS8/KCulQqdfmkjz2nmqf1DDedNFgEYFq2/uoae8LWJ7BO+tlMnf//qwvctcO2lHRwpH6rinNCCnOS6GqtZfO/iHKajsRcQ9PWeGaNXmenELwww/ejY6uCLLu01nZydxy6RJ+v+cULx+eeO5iR0UzK3JSSE+MmfC1vnLT4rDJ2D2L6rY+jBl/KvWnLyjk6S9fwsVLMrnv5XIu+f7L/Nuzh2juHsDhdPGPT+zj93tO8XdXLOXOTcuHbxhEhH/72CqSYqO44/G9Qe0r4p0JNdb0a/cmSOO/z701xYKpVrxldZ5n9tHEhnx+/vpxuvodfOnypRO6LlibS9xVkmfymgsNFhHImwgrDbDg7PR8+tD0LPLT48lKjuXRXZUMOU1Qe1iMxduLOFzXRVldB4vmJQYc7pksb05hIn+ILx9uZHFmIgsnMDvrKx9cxvKcZP7pf/dPKInaP+Tk3ar2kCROo+02clPjx7xDDXbdzcoFqfzPZ9bx/N9fygdWZPPz145zyfdf5uP/vYM/vFfDP3zwLO648uxRPcvMpFi+d+1KDtR08tOXywO2t6Kxe3grVX/cJT/G71mU1XaSEGMPKjd3VnYyZ2VPbPZRa88gv3zzBB8pybVkqBTcw4drCtJmxHa9Y9FgEWG6+ofo8exzUHpq/GBR2dpDanw0qQnR454XLBFhXWH68N1gySQ2U/LyjvuW1XZQVtfJCguS217enMIbQeYUmroG2Hm8hY3F2RP6PbFRdn70ydW09w7yz38+EPR171a2MehwhWShF7jLfoy1F7e3NHlhkDcQZ+ckc98Na3npjsvYXJLHobouvn7Vcr48zv4Qm1bmcu3aBfx0e/nwRIix+G6l6s/8lNiAw1BldZ0sz0kOWGnXa0tJHrtPtlEb5HqSB16toG/IyVc2hnZPjFHtWp3HwdrOkJfXDxUNFhHGd3rd/prx/xDHm/UyWesXurv68xJjyEud/BTP+cmxZCTG8PbJVk619lmSr/C1ZXUuTpfhuSAWPz2+5xRDTsOnzisIeO5I5+Sl8vcbz+KvpXU8FeTd646KFuxjbKQzGfnj7GtxsqWXhBg7mUkTG+5anJXEDz+xmrK7P8Stl41dmtvrX64+h6ykWO54fN+4VXnHmjbrlZMSR2PXwJgzrIwxHKqb2Ey6zZ69QJ4JYvZRY1c/v37rJB9ds4Cl860ZJvX6yKpcRGbuUJQGiwjjnRJZkBHPvuqOcdcPVLb0jpnInKxzPePCU0lug/tuvzg3hZfK3LkBq4NFcW4KizMTAw4/OF2GR3dVcdGSeeN+iI3n1ksXs7YwjW/+6UBQhfDeOt7CqgWpJMeFpgdYkJ5AQ+eA3w/pqlb3e2Ky/+/GmpY6Ump8NPdcV0J5Yzc/eP6I33NGbqXqT3ZKHA6XoXWMyqzVbX109Tsozg1+SHRRZiIrF6QEFcz/a3sFQ04T8p32/MlJjeO8hRk8ta92QuuCDtR0sOt44NmRU6XBIsJ4x2+vLM6hqWtgzOTfkNNFTXtfyHsW5+SlMC8xhosDrD0IRnFeCoOeJKgVayx8ibhLQu880TLusMYrR9zlLrZuKJr074qy2/jRJ1Yz4HDyT/9bOu4ffveAg32nQpOv8PJWn/VXtqOypSfk74mxXHpWFp/ZUMQv3jjBg68fH/X8yK1U/fHdXtWf4Q2zJvj+uXZtPqXVHfzntmNjnlPb3seju6r4xLr8CeWupuJjaxcMB9hgAsZ7VW3c8POdfPPPB3BaXDFYg0WE8f7RfNAznr5vjDHh2vY+nC4TsuS2V2yUnde+djmfu3jRlH+WtzcxLzGG+VPcwyIYW0rci5/GG354eGcl85Njh/99J2txVhJ3blrOK0eaeGz3qTHP232yFUcQG+lMxPBaixFJbpfLcKqtb0p7Wk/UNzcXc9XKHP71r4f471cqznhu5Faq/njXWjR2jREsajuxibtsxkTcdNFCPrZ2AT968Sj3vuD/g/mn28sxGG67wpoZUP58cn0BN15QyH+9UsF3/3po3ICx+2Qrn/nF26QnxPDQTedZMqXXlwaLCFPX0c+8xBjWFKQRZZMxE4j+isWFSmJsVEjemN6ZJStyU6Y0pBWsZdnJLM9J5i9jjAmfau3llaNNXH9eQUhq//zNhQu5aMk8vvN0GSebe/ye81ZFC9F2GZ6WHAqnN0E6s2dR39nPoMMV8qHJ8cRE2fjPG9Zy9eo8vv/cYe7zuZOvaOw5YytVf7IDlPwoq+tkUWYi8TETq31mtwk/+MRqPrW+gPteLuffnzt8xgdzVUsvj+8+xfXnFU6qXPxk2WzCv350JTddtJAH3zjBvzx10G++5q2KFj770NvMT47l8VsvnJY2arCIMA2d/eSkugsDnpWdTGm1/xlRw/tYTNOQw2QszkokNT6acwtDU8EzGFtW5/FOZZvfIZpH365CgOs9+yBMlc3zgWS3CR++73W+98yhUXfIb1W0sLYwfcIfduOZnxxLTJSN6hE9i8oQlKufjCi7jR9/ag0fO3cB9754lB96hljKG7vP2ErVH++uiWMOQ9V2Upw3uSncdpt7XcjWDYU88Opx7n66bDhg3PfyMew2mdZehZeI8K0txdxy6WJ+/VYl3/jT/jMCxhvHmvncr95mQVo8j926gZwpTDSZCGsmtivL1HX0D89CKslP5bmD9RhjRt2ZV7X0EBNlG974fiaKttt44SuXkhofmsRuMDaX5PKD54/w19Jabrn09KyeAYeTx3efYuOKbPLS4kP2+xakxfOnL13Mf247xoOvH+fXO05yw/mF3HrZYhKiozhQ666SGko2m5CfNnr6bJW32myIhyaDYbcJP7xuNTF2Gz/dXs6Q00X5iK1U/Ym228hMivE7DNXRNzTl/JLNJnznmpVE22388s2TDDld3HTRIv7wbjWfu3jRcM9muokIX79qOdF24f7tFQw6DPdcV8JrR5u49eF3WJyZyCN/e8Go3RutpMEiwtR39A3fiZfkp/HY7lOcau0bNdzknQkV7NzzcJnuP8aieYmU5Kfyl311ZwSL5w7U09IzOKUPnrEsyUriP65fy+0bz+K/Xynn4Z2VPLKrkvVFGRjDpPbbDmRBevyo6bOVLb1E2YS8tPB8ANpswveuXUW03cYDr7kT3sHkhsZaxX3Ik9xekTu1Ka0iwl2bi4nxtOuZ/fXERtn5wvsDTxG2kojw1SvPJsZu58cvHaWuo4/dJ1s5OyeZ3/6fC6a82n+idBgqgvQPOWnrHSLXp2cB/pPcVSGqNjsbbSnJY39Nxxl5hId3VlI0LyFghdmpWJSZyD3XrWb7V9/PJ9cX8E5lG0mxUayxYBiuICNhVIK7srWXBenxQU9/tYLNJtx9zTncfIl7gkQwJV7cwWJ0z2K4pliItuK986rlfPmKpbT2DHLTxQvJnMa79vHadfvGZXxt09nsqGihOC+VR/52w7QHCtCeRUTx/sF478bPyk4mJsrG/poOtqw+XW3TGENVa29IZ9jMJh8pyeW7zxzi6dJabrtiGYfrO9l9so3/9+Hl09ITK8hI4LvXruLvPrCMrv4hYqNCl68Y/h3pCbT1DtE94Bjet6LKgnU3kyEi/PNHVvDhVTmsWhA4UGanxPqdyFFW10lmUuyU9v8Y2a5/uPJsPlicbfm6n4n64vuXcvGSTJZlJ5FHB0leAAAdVklEQVQQE56Pbe1ZRBDvgrzcVPeYekyUjRW5Kew7deYfUlP3AL2DTu1ZjCEvLZ71Ren8ZZ97VtQjO6uIibKN2orTatkpcZatCvautfDtXUznGotARIR1RRmjdrTzJzsljubuwVGFCa3aA6UkPy2sva+xrC5IC1ugAA0WEcXbs/Cd/bA6P5UDNR1nzJY4XSxu+hOZkWLL6jyONHTxXlUbf3yvhs2rcsPStbfK8PRZT7Bo7x2ks98RluT2VHl70r5lxQcdLsobu2dcD2A202ARQbw9C99gUZKfRs+gk+PNp4uPWbnGYra4alUONoGv/H4v3QMObrQgsR1O3oV53hlRkfye8K7i9s1bVDR1M+h0TTm5rYKnwSKC1Hf0kxwbdcbeyd4kt+96i8rWXkTcJcWVf/OT49iweB4nW3pZkZsyrWs9pkN6QjSJMfbhGVGRsO5mLN6chG+w8Ca3z7G4TIw6zdJgISKbROSIiJSLyJ1+nr9DRMpEpFREtolIked4kYi8KyJ7ReSgiHzeynZGivqO/lELcJZkJZEQYz8jWFS19JCXGm9J4nQ28U4K2LqhcFpWkE8nEfHMiHL3LKqGS5NHXrDwvud9p8+W1XUSF21jUebkij2qibMsWyIiduB+4INANbBbRJ4yxpT5nPYesN4Y0ysiXwDuAT4F1AEXGmMGRCQJOOC5dubuDDIN6jpHBwu7TViZl3rG9NnK1tCXJp+Nrl3r3gp1uhPb0yU/PZ5qb8+ipZf5ybFhTZBOVkZCDFE2GdWzODsnxfJ6SOo0K3sW5wPlxpjjxphB4DHgGt8TjDHbjTHe6Ro7gXzP8UFjjPc2ItbidkaM+o4+cvwsYivJT6WstnN4tkiVBftYzEZx0XZuvKAoqBk5kSg/3b3WwhgT0TcQNpswPzl2uGdhjOFQfSfFmq+YVlb+lSwAfMttVnuOjeVm4FnvAxEpEJFSz8/4/lzvVTicLpq6BoYX5PlalZ/KgMPF0YYuugcctPQMBr0Tmpq9CjIS6Bl0L+R0r7GI3PdEdurphXl1Hf209w7pTKhpNiP6pCKyFVgPXOY9Zow5BZSISB7wJxF50hjTMOK6W4BbAAoLQ1P8baZq6h7AZdx/NCOtzncnZ/f75C0i9S5ShU6BZ4JDeWM39Z39Ef2eyE6OGy5pHsqV2yp4VvYsagDfweB8z7EziMhG4BvA1T5DT8M8PYoDwPv8PPczY8x6Y8z6rKyskDV8Jjq9IG90sCial0BKXBT7qjuG11hEYiJThZZ3+uyOimYgsm8gslNih3sWZXWdiMDZORosppOVwWI3sExEFolIDHA98JTvCSKyFngAd6Bo9DmeLyLxnu/TgUsA/3szzhEN3jUWKaOnw4oIJflp7K9pj+gpkiq0hoNFuXvLzUi+gZifEkdnv4O+QSdltZ0UZSScMYVcWc+yYGGMcQC3Ac8Dh4DHjTEHReRuEbnac9oPgCTgCc80WW8wWQHsEpF9wKvAD40x+61qayTwtyDP16r8VA7XdXGsoZuMxJiQ7eesIldSbBTpCdG8W9UGRPaK/pyU02stDtVbU+ZDjc/S0GyMeQZ4ZsSxu3y+3zjGdS8CJVa2LdLUd/YTE2UjPcF/EFidn4rDZdh2uIGFEfyhoEIrPz2B/TUdJHsCR6TylvyoaOqmsqWXT6zLD3OL5p7ZOWdwFqrv6Cc3NW7MxWMlniR3e++QDkGpYd6CgoXzEiJ64aG35MerR5sATW6HgwaLCFHf0T/uRkG5qXFkJrkL4Wm1WeXlLSgY6TcQ3lmALx92pzaLcye3laqaPA0WEaKus8/vTCgvb5IboFCHoZRHvufGIZLXWAAkx0YRH22nuq2P9ITo4Z6Gmj4aLCKAMYaGjgG/q7d9rVrgvtuK9LtIFTretRaR/p4QkeEAUZyXEtFDapFKg0UEaO0ZZNDpGnMmlNdHSnK5ZGmmrmxVw9YUpHHRknmWbhc7XeZ7bpb0/R0eOlE5Aoy3IM/XWdnJPPy3F0xHk1SESEuI4dH/uyHczQgJb89ak9vhoT2LCDBy722l5iLvMNQK7VmEhfYsIsDIvbeVmosuXDKPd6vaWZKle1iEgwaLCFDf0Y/dJmQl6wwQNXddsTybK5Znh7sZc5YOQ0WA+s5+spJidaMXpVTYaLCIAP62U1VKqemkwSIC1HWMvyBPKaWspsEiAjR0DuhMKKVUWGmwmOG6+ofoHnBoz0IpFVYaLGa4+gD7WCil1HTQYDHD1Xd6d8jTYKGUCh8NFjOcLshTSs0EGixmOO8w1HwtyayUCiMNFjNcfWc/GYkxxEXbw90UpdQcpsFihqvv6Nd8hVIq7DRYzHB1nr23lVIqnDRYzHANnf3D+w8rpVS4aLCYwfqHnLT2DJKrw1BKqTDTYDGDeTc90gV5SqlwszRYiMgmETkiIuUicqef5+8QkTIRKRWRbSJS5Dm+RkTeEpGDnuc+ZWU7Zypdva2UmiksCxYiYgfuB64CioEbRKR4xGnvAeuNMSXAk8A9nuO9wN8YY84BNgH/ISJpVrV1pvKu3tYEt1Iq3KzsWZwPlBtjjhtjBoHHgGt8TzDGbDfG9Hoe7gTyPcePGmOOeb6vBRqBLAvbOiOd7lno6m2lVHhZGSwWAKd8Hld7jo3lZuDZkQdF5HwgBqgIaesiQF1HP0mxUSTF6u63SqnwmhGfQiKyFVgPXDbieC7wW+CzxhiXn+tuAW4BKCwsnIaWTi/dIU8pNVNY2bOoAQp8Hud7jp1BRDYC3wCuNsYM+BxPAf4KfMMYs9PfLzDG/MwYs94Ysz4ra/aNUtV36oI8pdTMYGWw2A0sE5FFIhIDXA885XuCiKwFHsAdKBp9jscAfwR+Y4x50sI2zmj1Hf26Q55SakYIGCxExC4iX5noDzbGOIDbgOeBQ8DjxpiDInK3iFztOe0HQBLwhIjsFRFvMPkkcClwk+f4XhFZM9E2RLLuAQeNXVoXSik1MwTMWRhjnCJyA/Djif5wY8wzwDMjjt3l8/3GMa57GHh4or9vNvnuX8sA2FicHeaWKKVU8AnuN0Xkp8DvgR7vQWPMu5a0ao57+XADv3v7FJ+/bAlrCubc8hKl1AwUbLDwDgHd7XPMAFeEtjmqrWeQf/rf/SzPSeYrH1wW7uYopRQQZLAwxlxudUOU2zf/fID23kF+9bnziI3SDY+UUjNDULOhRCRVRO4VkT2erx+JSKrVjZtr/rKvlqdL67j9A8s4J0//eZVSM0ewU2cfArpwz1L6JNAJ/NKqRs1FDZ39fPPPB1hTkMbnL1sS7uYopdQZgs1ZLDHGfNzn8bdFZK8VDZqLjDH80/+W0j/k5N5PribKrpXjlVIzS7CfSn0icon3gYhcDPRZ06S557Hdp3jlSBN3blrO4qykcDdHKaVGCbZn8XngNz55ijbgs9Y0aW6pbuvlX58u46Il8/ibCxeGuzlKKeVXwGAhIjbgbGPMak+9JowxnZa3bI54/mADPYNOvnftKmw2CXdzlFLKr4DDUJ5qr1/zfN+pgSK0TjR3kxIXRdG8hHA3RSmlxhRszuIlEfmqiBSISIb3y9KWzREnm3tZlJmIiPYqlFIzV7A5C+8e2F/yOWaAxaFtztxzormH8xamh7sZSik1rmBzFluNMW9OQ3vmlP4hJ7UdfSzMzA93U5RSalzB5ix+Og1tmXMqW3oxBhZlJoa7KUopNa5gcxbbROTjogPrIXWi2V3AV4OFUmqmCzZY3Ao8DgyISKeIdImIzoqaIm+wWKjBQik1wwWb4E4FbgQWGWPuFpFCINe6Zs0NJ5t7yEyKISUuOtxNUUqpcQXbs7gf2ADc4HncheYxpuxEc48OQSmlIkKwweICY8yXgH4AY0wbEGNZq+aIEy09LJynwUIpNfMFGyyGRMSOe20FIpIFuCxr1RzQ1T9EU9cAi7I0WCilZr5gg8V9wB+B+SLyXeAN4HuWtWoOqGzpBWCR9iyUUhEg2G1VHxGRd4APAAJ81BhzyNKWzXLHvdNmtWehlIoAwc6GwhhzGDhsYVvmlJOeYFGUocFCKTXz6ZZsYXKiuYe81DjiY+zhbopSSgVkabAQkU0ickREykXkTj/P3yEiZSJSKiLbRKTI57nnRKRdRJ62so3hcqK5R4eglFIRw7Jg4Zk9dT9wFVAM3CAixSNOew9Yb4wpAZ4E7vF57gfAZ6xqX7idaNZps0qpyGFlz+J8oNwYc9wYMwg8Blzje4IxZrsxptfzcCeQ7/PcNtyL/2adtp5BOvqGdEGeUipiWBksFgCnfB5Xe46N5WbgWQvbM2Mc1wKCSqkIE/RsKCuJyFZgPXDZBK+7BbgFoLCw0IKWWUOrzSqlIo2VPYsaoMDncb7n2BlEZCPwDeBqY8zARH6BMeZnxpj1xpj1WVlZU2rsdDrZ3IPdJhRk6L7bSqnIYGWw2A0sE5FFIhIDXA885XuCiKwFHsAdKBotbMuMcqK5h4L0eKLtOnNZKRUZLPu0MsY4gNuA54FDwOPGmIMicreIXO057QdAEvCEiOwVkeFgIiKvA08AHxCRahH5kFVtnW4nmnt0DwulVESxNGdhjHkGeGbEsbt8vt84zrXvs7BpYWOM4WRLDxcszgh3U5RSKmg6DjLNGrsG6B10anJbKRVRNFhMs+NNOhNKKRV5NFhMs5Mtnn23dfW2UiqCaLCYZieae4iJspGXFh/upiilVNA0WEyzE809FGUkYLdJuJuilFJB02AxCadae/nuX8twusyErz3R3KP5CqVUxNFgMQlP7DnFz18/QVVrb+CTfThdhqqWXg0WSqmIo8FiEkprOgBo6ppQdRJq2/sYdLo0WCilIo4GiwkyxlBa7Q4Wzd0TCxZaQFApFak0WExQTXsfrT2DwMR7FhoslFKRSoPFBHl7FTC5YJEYYycrOTbUzVJKKUvNiP0sIsm+6nai7UJSbNSkgsXCzEREdNqsUiqyaLCYoP3VHazITcHpMjRNMGdxsqWHVQtSLWqZUkpZR4ehJsDlMuyv7qAkP5Ws5NgJ9SwGHS5Oteq0WaVUZNJgMQEnW3roGnBQsiCNrKSJBYuq1l5cRpPbSqnIpMFiArzJ7ZKCVDKTY2nuHsAV5Cruk56ZULrpkVIqEmmwmIDS6g7io+0szUoiKykWh8vQ3jcU1LXeabOLNVgopSKQBosJKK1u55y8FKLstuHpr8EuzDvR0kNaQjRpCTFWNlEppSyhwSJIDqeLA7UdlOSnAQwHi2DzFpUtPRTpHhZKqQilwSJI5U3d9A+5KMl3T32daLCoa+8nX/ewUEpFKA0WQSo95UluTyJYGGOoae8jNzXOugYqpZSFNFgEqbSmneTYqOHtUJNjo4iNsgW1MK+td4gBh4tc7VkopSKUBosglVZ3sCo/FZtnhzsRCXphXm17HwAL0rRnoZSKTBosgjDgcHKorpNV+WeW6gg2WNR19AOQm6o9C6VUZLI0WIjIJhE5IiLlInKnn+fvEJEyESkVkW0iUuTz3GdF5Jjn67NWtjOQI/VdDDkNqz0zobyCXcVd1+HuWeRqz0IpFaEsCxYiYgfuB64CioEbRKR4xGnvAeuNMSXAk8A9nmszgG8BFwDnA98SkXSr2hrIvuozk9temcmxQeUsatr7iLYLmYlamlwpFZms7FmcD5QbY44bYwaBx4BrfE8wxmw3xng3st4J5Hu+/xDwojGm1RjTBrwIbLKwrePaX91ORmIMC0YkqLOSYmnrHWTI6Rr3+rr2fnJS44bzHUopFWmsDBYLgFM+j6s9x8ZyM/DsJK+1VKmn0uzIfSiykmMxhuGd88ZS19FHnuYrlFIRbEYkuEVkK7Ae+MEEr7tFRPaIyJ6mpiZL2tY36ORoQxclfvahCHatRW17P3k6bVYpFcGsDBY1QIHP43zPsTOIyEbgG8DVxpiBiVxrjPmZMWa9MWZ9VlZWyBru62BtBy7DcJkPX8EEC6fL0NDZrwvylFIRzcpgsRtYJiKLRCQGuB54yvcEEVkLPIA7UDT6PPU8cKWIpHsS21d6jk270jGS2+DOWcD4waKpawCHy+iCPKVURLNsW1VjjENEbsP9IW8HHjLGHBSRu4E9xpincA87JQFPePIBVcaYq40xrSLyHdwBB+BuY0yrVW0dT2l1OzkpccxPGd0zGO5ZjDMjqtYzbTZPexZKqQhm6R7cxphngGdGHLvL5/uN41z7EPCQda0Ljje57U9ctJ3kuKhxexZ17e4FeZqzUEpFshmR4J6pOvuHON7cM2awgMCruOuGexYaLJRSkUuDxTgODOcrRie3vQKt4q5t7ychxk5KvKWdOKWUspQGi3GU1riDxSo/02a9vHtxj6XWU5p85BoNpZSKJHP+drejd4hPP7jT73O17X0UZiSQnjj2VqhZSbG8FmAYSvMVSqlIN+eDhdgYcw1EbmocV56TM+71WcmxdA046Bt0Eh9jH/V8bUc/y3NSQtJWpZQKlzkfLFLionnws+dN+nrv9Nnm7gEKMhLOeG7Q4aK5e0CrzSqlIp7mLKbIGywa/QxFNXT2Y4zOhFJKRT4NFlM03irumnbdx0IpNTtosJii+eOs4h7e9Eh7FkqpCKfBYooyEmMQ8d+zqB1eva09C6VUZNNgMUVRdhvzEmP8rrWo6+gjLSGahJg5P49AKRXhNFiEQOYYq7hr2/t1CEopNStosAiBsepD1bb3abVZpdSsoMEiBMaqD1XX0a8zoZRSs4IGixDISo6lqXsAY8zwsd5BBx19Q1rqQyk1K2iwCIGs5FgGHS46+x3Dx4ZnQmnOQik1C2iwCAF/e3HXehfkac5CKTULaLAIAX+ruIc3PdJhKKXULKDBIgT87cVd296PCGT72btbKaUijQaLEBiuPDuiZ5GVFEtMlP4TK6Uin36ShUBqfDTRdjmjZ+GeNqtDUEqp2UGDRQiIyKi1FjW6IE8pNYtosAiRTJ9V3MYY6rTUh1JqFtFgESK+PYuOviH6hpxabVYpNWtYGixEZJOIHBGRchG508/zl4rIuyLiEJHrRjz3fRE54Pn6lJXtDAXvKm7wLU2uPQul1OxgWbAQETtwP3AVUAzcICLFI06rAm4CHh1x7UeAc4E1wAXAV0Ukxaq2hkJWciwt3QM4XcZn0yPtWSilZgcrexbnA+XGmOPGmEHgMeAa3xOMMSeNMaWAa8S1xcBrxhiHMaYHKAU2WdjWKctKjsVloLVncHj1tvYslFKzhZXBYgFwyudxtedYMPYBm0QkQUQygcuBghC3L6R8V3HXdvQTZRMyPceUUirSzcgt3IwxL4jIecAOoAl4C3COPE9EbgFuASgsLJzWNo40vDCve4C69j6yU+Kw2ySsbVJKqVCxsmdRw5m9gXzPsaAYY75rjFljjPkgIMBRP+f8zBiz3hizPisra8oNngrfYoK1Hf0s0CEopdQsYmWw2A0sE5FFIhIDXA88FcyFImIXkXme70uAEuAFy1oaAt4hp6buAeo6+nTTI6XUrGLZMJQxxiEitwHPA3bgIWPMQRG5G9hjjHnKM9T0RyAd2CIi3zbGnANEA6+LCEAnsNUY4/D/m2aGxNgoEmPsNHT2U9+hC/KUUrOLpTkLY8wzwDMjjt3l8/1u3MNTI6/rxz0jKqJkJsdyqK6TIafRBXlKqVlFV3CHUFZSLAdqOgG0Z6GUmlU0WIRQVnIs3QPu0TLtWSilZhMNFiHknREFuve2Ump20WARQt6FeXHRNtISosPcGqWUCh0NFiHk7VnkpcbjmcmllFKzggaLEPIGC11joZSabTRYhJBvz0IppWYTDRYhdLpnocFCKTW7aLAIoezkOLZuKOSqlTnhbopSSoXUjKw6G6lsNuFfP7oq3M1QSqmQ056FUkqpgDRYKKWUCkiDhVJKqYA0WCillApIg4VSSqmANFgopZQKSIOFUkqpgDRYKKWUCkiMMeFuQ0iISBNQGeC0TKB5GpozU83l1z+XXzvM7devr318RcaYrEA/aNYEi2CIyB5jzPpwtyNc5vLrn8uvHeb269fXHprXrsNQSimlAtJgoZRSKqC5Fix+Fu4GhNlcfv1z+bXD3H79+tpDYE7lLJRSSk3OXOtZKKWUmoQ5EyxEZJOIHBGRchG5M9ztsZqIPCQijSJywOdYhoi8KCLHPP9ND2cbrSIiBSKyXUTKROSgiNzuOT7rX7+IxInI2yKyz/Pav+05vkhEdnne/78XkZhwt9UqImIXkfdE5GnP47n02k+KyH4R2SsiezzHQvK+nxPBQkTswP3AVUAxcIOIFIe3VZb7FbBpxLE7gW3GmGXANs/j2cgB/IMxphjYAHzJ8/97Lrz+AeAKY8xqYA2wSUQ2AN8HfmyMWQq0ATeHsY1Wux045PN4Lr12gMuNMWt8psyG5H0/J4IFcD5Qbow5bowZBB4DrglzmyxljHkNaB1x+Brg157vfw18dFobNU2MMXXGmHc933fh/uBYwBx4/cat2/Mw2vNlgCuAJz3HZ+VrBxCRfOAjwIOex8Icee3jCMn7fq4EiwXAKZ/H1Z5jc022MabO8309kB3OxkwHEVkIrAV2MUdev2cYZi/QCLwIVADtxhiH55TZ/P7/D+BrgMvzeB5z57WD+8bgBRF5R0Ru8RwLyfte9+Ceo4wxRkRm9VQ4EUkC/hf4e2NMp/sm0202v35jjBNYIyJpwB+B5WFu0rQQkc1AozHmHRF5f7jbEyaXGGNqRGQ+8KKIHPZ9cirv+7nSs6gBCnwe53uOzTUNIpIL4PlvY5jbYxkRicYdKB4xxvzBc3jOvH4AY0w7sB24EEgTEe/N4Wx9/18MXC0iJ3EPNV8B/IS58doBMMbUeP7biPtG4XxC9L6fK8FiN7DMMysiBrgeeCrMbQqHp4DPer7/LPDnMLbFMp5x6l8Ah4wx9/o8Netfv4hkeXoUiEg88EHcOZvtwHWe02blazfGfN0Yk2+MWYj7b/xlY8yNzIHXDiAiiSKS7P0euBI4QIje93NmUZ6IfBj3eKYdeMgY890wN8lSIvI74P24q042AN8C/gQ8DhTirtD7SWPMyCR4xBORS4DXgf2cHrv+f7jzFrP69YtICe4kph33zeDjxpi7RWQx7rvtDOA9YKsxZiB8LbWWZxjqq8aYzXPltXte5x89D6OAR40x3xWReYTgfT9ngoVSSqnJmyvDUEoppaZAg4VSSqmANFgopZQKSIOFUkqpgDRYKKWUCkiDhVJKqYA0WCg1zTxlpDMnee1NIpIXip+l1ERosFAqstwE5AU6SalQ02Ch5iwRWSgih0XkVyJyVEQeEZGNIvKmZ6OY8z1fb3k209khImd7rv2KiDzk+X6ViBwQkYQxfs88EXnBsxnRg4D4PLfVs1nRXhF5wLP3CiLSLSI/9lyzzVPG4zpgPfCI5/x4z4/5soi869n0Zk4UDVTTT4OFmuuWAj/CXZl1OfBp4BLgq7hLhBwG3meMWQvcBXzPc91PgKUici3wS+BWY0zvGL/jW8AbxphzcJdjKAQQkRXAp4CLjTFrACdwo+eaRGCP55pXgW8ZY54E9gA3eja36fOc22yMORf4b0+7lQo5LVGu5roTxpj9ACJyEPeOYkZE9gMLgVTg1yKyDPdeAdEAxhiXiNwElAIPGGPeHOd3XAp8zHPdX0WkzXP8A8A6YLenfHo8pyuCuoDfe75/GPgDY/M+94739ygVahos1FznW1DO5fPYhfvv4zvAdmPMtZ6NlF7xOX8Z0M3kcwgC/NoY8/Ugzh2viJu3zU70b1pZRIehlBpfKqf3P7jJe1BEUoH7cPca5nnyCWN5DffwFiJyFZDuOb4NuM6zUQ0ikiEiRZ7nbJwuq/1p4A3P911A8hRej1KTosFCqfHdA/ybiLzHmXftPwbuN8YcBW4G/t37oe/Ht4FLPcNcHwOqAIwxZcA/494GsxT3Fqi5nmt6gPNF5ADuTXzu9hz/FfA/IxLcSllOS5QrNQOJSLcxJinc7VDKS3sWSimlAtKehVIhIiKfA24fcfhNY8yXwtEepUJJg4VSSqmAdBhKKaVUQBoslFJKBaTBQimlVEAaLJRSSgWkwUIppVRA/x/LM3kv95ak7gAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"line_df = pd.DataFrame({'max_depth' : rnge, 'error' : errors})\n",
"sns.lineplot(x='max_depth', y='error', data=line_df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(i) Produce a pruned tree corresponding to the optimal tree size obtained using cross-validation. If cross-validation does not lead to selection of a pruned tree, then create a pruned tree with five terminal nodes. (j) Compare the training error rates between the pruned and un- pruned trees. Which is higher? \n",
"(k) Compare the test error rates between the pruned and unpruned trees. Which is higher?\n",
"\n",
".. pruning not easily available in python world"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"10. We now use boosting to predict Salary in the Hitters data set."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
AtBat
\n",
"
Hits
\n",
"
HmRun
\n",
"
Runs
\n",
"
RBI
\n",
"
Walks
\n",
"
Years
\n",
"
CAtBat
\n",
"
CHits
\n",
"
CHmRun
\n",
"
CRuns
\n",
"
CRBI
\n",
"
CWalks
\n",
"
League
\n",
"
Division
\n",
"
PutOuts
\n",
"
Assists
\n",
"
Errors
\n",
"
Salary
\n",
"
NewLeague
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1
\n",
"
315
\n",
"
81
\n",
"
7
\n",
"
24
\n",
"
38
\n",
"
39
\n",
"
14
\n",
"
3449
\n",
"
835
\n",
"
69
\n",
"
321
\n",
"
414
\n",
"
375
\n",
"
1
\n",
"
1
\n",
"
632
\n",
"
43
\n",
"
10
\n",
"
6.163315
\n",
"
1
\n",
"
\n",
"
\n",
"
2
\n",
"
479
\n",
"
130
\n",
"
18
\n",
"
66
\n",
"
72
\n",
"
76
\n",
"
3
\n",
"
1624
\n",
"
457
\n",
"
63
\n",
"
224
\n",
"
266
\n",
"
263
\n",
"
0
\n",
"
1
\n",
"
880
\n",
"
82
\n",
"
14
\n",
"
6.173786
\n",
"
0
\n",
"
\n",
"
\n",
"
3
\n",
"
496
\n",
"
141
\n",
"
20
\n",
"
65
\n",
"
78
\n",
"
37
\n",
"
11
\n",
"
5628
\n",
"
1575
\n",
"
225
\n",
"
828
\n",
"
838
\n",
"
354
\n",
"
1
\n",
"
0
\n",
"
200
\n",
"
11
\n",
"
3
\n",
"
6.214608
\n",
"
1
\n",
"
\n",
"
\n",
"
4
\n",
"
321
\n",
"
87
\n",
"
10
\n",
"
39
\n",
"
42
\n",
"
30
\n",
"
2
\n",
"
396
\n",
"
101
\n",
"
12
\n",
"
48
\n",
"
46
\n",
"
33
\n",
"
1
\n",
"
0
\n",
"
805
\n",
"
40
\n",
"
4
\n",
"
4.516339
\n",
"
1
\n",
"
\n",
"
\n",
"
5
\n",
"
594
\n",
"
169
\n",
"
4
\n",
"
74
\n",
"
51
\n",
"
35
\n",
"
11
\n",
"
4408
\n",
"
1133
\n",
"
19
\n",
"
501
\n",
"
336
\n",
"
194
\n",
"
0
\n",
"
1
\n",
"
282
\n",
"
421
\n",
"
25
\n",
"
6.620073
\n",
"
0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns \\\n",
"1 315 81 7 24 38 39 14 3449 835 69 321 \n",
"2 479 130 18 66 72 76 3 1624 457 63 224 \n",
"3 496 141 20 65 78 37 11 5628 1575 225 828 \n",
"4 321 87 10 39 42 30 2 396 101 12 48 \n",
"5 594 169 4 74 51 35 11 4408 1133 19 501 \n",
"\n",
" CRBI CWalks League Division PutOuts Assists Errors Salary \\\n",
"1 414 375 1 1 632 43 10 6.163315 \n",
"2 266 263 0 1 880 82 14 6.173786 \n",
"3 838 354 1 0 200 11 3 6.214608 \n",
"4 46 33 1 0 805 40 4 4.516339 \n",
"5 336 194 0 1 282 421 25 6.620073 \n",
"\n",
" NewLeague \n",
"1 1 \n",
"2 0 \n",
"3 1 \n",
"4 1 \n",
"5 0 "
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"hit_df = pd.read_csv('hitters.csv')\n",
"hit_df = hit_df.drop(hit_df.columns[0], axis = 1)\n",
"hit_df = hit_df.dropna()\n",
"hit_df.Salary = np.log(hit_df.Salary)\n",
"hit_df.League = hit_df.League.map({'N':1, 'A':0})\n",
"hit_df.Division = hit_df.Division.map({'W':1, 'E':0})\n",
"hit_df.NewLeague = hit_df.NewLeague.map({'N':1, 'A':0})\n",
"hit_df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) Create a training set consisting of the first 200 observations, and a test set consisting of the remaining observations."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"train, test = model_selection.train_test_split(hit_df, test_size=hit_df.shape[0] - 200)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) Perform boosting on the training set with 1,000 trees for a range of values of the shrinkage parameter λ. Produce a plot with different shrinkage values on the x-axis and the corresponding training set MSE on the y-axis."
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"def hitters_boosting_test_error(lr):\n",
" train_x = train.drop('Salary', axis=1)\n",
" train_y = train.Salary\n",
" reg = ensemble.GradientBoostingRegressor(n_estimators=1000, learning_rate=lr).fit(train_x, train_y)\n",
" test_x = test.drop('Salary', axis=1)\n",
" test_y = test.Salary\n",
" preds = reg.predict(test_x)\n",
" return np.sqrt(metrics.mean_squared_error(test_y, preds))\n",
"\n",
"lambdas = np.linspace(0.01, 1, 100)\n",
"test_errors = [hitters_boosting_test_error(lr) for lr in lambdas]"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [],
"source": [
"def hitters_boosting_training_error(lr):\n",
" train_x = train.drop('Salary', axis=1)\n",
" train_y = train.Salary\n",
" reg = ensemble.GradientBoostingRegressor(n_estimators=1000, learning_rate=lr).fit(train_x, train_y)\n",
" preds = reg.predict(train_x)\n",
" return np.sqrt(metrics.mean_squared_error(train_y, preds))\n",
"\n",
"lambdas = np.linspace(0.01, 1, 100)\n",
"train_errors = [hitters_boosting_training_error(lr) for lr in lambdas]"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lOW5//HPlZ1ACJCFPYQlbIICBhRQ3C0uR2pdCtaqrS2ntrZW6zm1P1vr0r3Hc7qotaK2tVVxqVVUlLohm+x7WJMASVgDWSFkv39/zCTOZCFDyGQC+b5fL17MPHPP81xPCHPNvZtzDhERkTphoQ5AREQ6FiUGERHxo8QgIiJ+lBhERMSPEoOIiPhRYhARET9KDCIi4keJQURE/CgxiIiIn4hQB9AaiYmJLjU1NdRhiIicVtasWXPYOZfUUrnTMjGkpqayevXqUIchInJaMbM9gZRTU5KIiPhRYhARET9KDCIi4keJQURE/CgxiIiIHyUGERHxo8QgIiJ+lBhERNpJzpEyPtxyMNRhtEiJQUSkncxZnM3sv6/mUGl5qEM5ISUGEZF2sr+4nFoHb2/YH+pQTkiJQUSkneR7awpvrd8b4khOTIlBRKSdHCypICo8jI15xWTlHw11OM1SYhARaQc1tY78oxXMGNcPM3hrXcetNSgxiIi0g4JjldTUOsYOiGfq0ETeXL8P51yow2qSEoOISDs4WOLpX0iOi2HGuH7kFJSxNqcoxFE1TYlBRKQd5JdWAJDcPZrpY/oQHRHWYTuhlRhERNpBXY2hd/cY4mIiuXx0b97esI/qmtoQR9aYEoOISDs4WOKpMSR1iwbgC2f1obCsis37SvzKFZVV8uQnmdTUhq7/QYlBRKQdHCotp1fXKKIiPB+7k4ckALAs67BfuRdX5PDbBdvZmBe6/gclBhGRdnCwpILkuOj650lx0YzoHcdnWUf8yi3cfgiAnIKydo3PlxKDiEg7yC8tJ7l7jN+xyUMTWLW7gIrqGsDTjLRmTyEAe46cwYnBzKab2XYzyzSzB5opc7OZbTGzDDN7KdgxiYi0t4MlFfT2qTEATB2WSHlVLeu8w1YX7TxMrYMwC21iiAjmyc0sHHgSuALIA1aZ2Tzn3BafMmnAj4CpzrlCM0sOZkwiIu2tbtZzcnf/xDBpcC/CDJZlHub8IQl8su0QvbpGMSSxK7lncFPSJCDTOZftnKsE5gIzGpT5JvCkc64QwDl3KMgxiYi0q7pZz70bNCXFd4lk7IAeLMs6Qk2tY+H2Q1w8PInUxK7sKTgWomiDnxj6A7k+z/O8x3wNB4ab2VIzW25m05s6kZnNNrPVZrY6Pz8/SOGKiLQ931nPDU0ZmsD63CI+yzpCYVkVF49MZlCvWA6WVFBeVdPeoQIdo/M5AkgDLgZmAXPMrEfDQs65Z5xz6c659KSkpHYOUUSk9XxnPTc0ZWgC1bWO3yzYRniYcVFaEikJsUDoRiYFOzHsBQb6PB/gPeYrD5jnnKtyzu0CduBJFCIiZwTfWc8NpQ/qVb8U97kpPYmPjSSllzcxhKgDOtiJYRWQZmaDzSwKmAnMa1DmTTy1BcwsEU/TUnaQ4xKRDuqRtzM67BpCgXrwX5v43svr6p83nPXsq0tUOONTPI0kF4/0tIYMSugKwJ4zscbgnKsG7gYWAFuBV51zGWb2qJld5y22ADhiZluAT4D/cs4dafqMInImKy6r4q/LdvP4v3e0y5LUtbWOdza27XpFB0vKmbsql/c27+doRTXQeNZzQxcMSwTg0pGeQZk9YyOJi44g50hoOqCD3sfgnJvvnBvunBvqnPu599hDzrl53sfOOXefc260c26sc25usGMSkY5pbU4hznna1lftLgz69T7edoi7X1rH/M0H2uycL63IoabWUVXjWJHt+Y7bcNZzQ3dMTWXObemM7NMdADNjYK/YM7PGICJyMlbtLiAizIiNCuefa/ICek/BsUoyD5W26npLMj3rFC3LPNxCycBU1dTy8socpg5LICYyjMU7Pedtatazr7iYSK4Y3dvv2KCE2DO281lEJGCrdxdyVv94rh7bl3c37aessvqE5fMKy5jx5BJmzVnRquvVrVO0NKttEsMHWw5yqLSCOy8YzPlDEli00zO0vqlZzy1JSYglr+B4SFZZVWIQkQ6horqG9XlFTBzUkxvPHcDRimoWZDTfxJNXWMbMZ5aTW3Cc/NIKjhytOKnr5ZdWsP1gKQN7dSG34HibzDR+4bPdDOjZhYuGJ3NhWhLZ+cfILShrctZzSwb16kplTS0HvCOa2pMSg4h0CJv3FlNZXUt6ai8mpfZiYK8u/HNN06OT6pJCyfEq/usLIwDYeejoSV3vM2/7//cvGw7A0lNsTtp5sJTl2QV85bxBhIcZ09I8Hcpvrtvb5KznltQNWd0Tgg5oJQaRDsQ51yF39GoPdZ3N6ak9CQszbpgwgKVZh9lbdLxR2R+9sYni41X84xvncf14z2IKJ5sYlmUeJi4mghnj+pEcF82yrFMbDPmP5XuICg/j5vQBAAxL7kaf7jH8c62nr6SpWc8nMsg7yS0UayYpMYh0IK+vyWPSLz7iWMWJ29bPRKt3FzAksSuJ3rH+N0wYgHPwRoNO6D1HjrF452FmXziEswf0oG98DF2jwsk8eHId0MuyjnD+kAQiwsOYMjSBZVlHWhwiW1ZZzcPzMliy07924ZxjQcZBLh+dTII3fjPjwrREdnsnqZ1sU1Lf+Bgiwiwkq6wqMYh0IAt35FNwrJLVe4I/VLMjqa11rN5TSHpqz/pjA3vFcsGwRF5YvsdvzaBXVuUSZnBTumdRBTNjWO+4k6ox5BaUkVNQxtShnl3UpgxL5PDRCnYcbP4cBccqmTVnBX9dtps/L8pqcL7jHCgpZ/LQRL/jFw7/fPmek21KiggPo3/PLiEZsqrEINIGfrtgGx9vO3jK51nvXZe/4XaPZ7qs/KMUlVWRntrL7/h3Lx1GfmkFL63IATzDQV9bk8elI5PpE//5B21acreTSgx1o5GmeieWTfEmiLp+hgPF5cx+YTX3v7aBt9bvZfPeYm58ehnb9pcwPqWH3+Y6AMt3ec533mD/+C8YloiZ53FTs55bktIrNiTLYigxiJyirftLePKTLP6xPOeUznOotLy+PX35KbZ3n27q+hcmNkgM5w1JYMrQBP70aRblVTV8vO0Q+aUVfHliil+5tORu5JdWUFxWFdD1lmYdJikummHJ3QAY0DOWQQmxLMs6wqHScm55djlLMg/zwZaD3DN3Pdf+cQmHSyv4+53n8e2Lh/ltrgOwclcBPWMjGZbUze86vbpGMbZ//AlnPZ9IqOYyBHWjHpHO4O/L9wCwZV9JwO/5dEc+mYeOcucFg+uP1dUWLhiWyLKsw5SUV9E9JrJVMTnn+MGrG/jShAFckJbY8htCbPXuAhK7RZHq7XD1dc9laXz5meW8uCKHJTvz6d09mktG+K+wnNbb84GcmV/KuYN6NTqHL+ccy7KOMGVoAlb3dR6YMjSRdzbs49ZnV7C/qJwX7pzEhJSebNpbzNo9hVw0IomhSd0oPl7l2VzH20cBnsQwMbUXYWHW6HrfvTSN3YdbN7JoUK+uFB+vorisivjY1v0utIZqDCKnoKS8ijfX7SU2KpwDJeUBjaU/WlHND15dzy/mb6XwWGX98fW5RUSEGd+4cDC1DlZmF7Q6rrzC47yxbi8PvbU55KOc9jUxqshXdU0tK3YVkD6ol98HdZ26WsMTH+/k0x353HTuQCLC/T+60pLjANh5gj6COhvziskvrahvPqozdVgCpRXV7DlSxnO3pzMxtRfhYca4gT34+gWDGeqtDdRvruNtdtpffJycgjLOG5LQ6FoAV4zuzTenDWkxrqbULb+9o5Uzu1tLiUHkFLyxJo+yyhruucyzUvyW/S3XGp5emMXho54dvT7e9vmGhetyihjVtzvnD0kgOiLslIZPZnhrL9mHj/GvdaFbqfTllTlM+dXHvLNxX7NlfjF/G3uLjnPduH7Nlvn+5cMpLKui1sGXJw5s9Hr/Hl2IiQxrsZ9h4fZDfPW5FfSIjeSSkf67CE8bnsTlo3rz7O3pTBl24lpW3eY6xyqqWbnLk8Ab9i+0hfMHJxAXHcGcRe274LQSg0grOef4+/I9nDOwBzd7R8i01Jy0r+g4cxZnM2NcP/rGx9TP7K2pdWzMK2LcwB7ERIaTntrzlDqgt+wvIcxgZJ84fv/RTiqr26bWsHlvMf/wNp21ZFNeMT99KwOAuStzmyzz2upcnl+6i69NTeXqsX2bPdekwb24fFRvrhrTh4G9Gjc3hYUZw07QAe2c46mFmXztr6vo16ML875zQaN5Bd1jInn29nQuTGt5I7CpQxOprnWs3F3Ail0FdIuOYFTf7i2+72TFx0byjQuH8O8tB9mYV9TyG9qIEoNIK32WdYSs/GPcdv4genaNol98TP039eb8dsF2AP57+kiuHN2bRTvzOV5ZQ+ahoxyrrKlfl3/ykAS2HSilwKep6WRs2VfC0KRu/PCqkeQVHufV1U1/MJ+MNXsKmfnMcn785uZGi9btKzrOj97YyNLMwzjnKCqr5K4X15DYLYrbJw9iadbhRk1K63IKefBfm5kyNIEHrx7V4vWfvT2dP916brOvpyXHNTuX4c+LsvnN+9u5Zmxf3vj2lPommtZKT+1JVHgYyzIPs3JXAempPQlvon+hLXz9glR6xEby+L93BOX8TVFiEGmlFz7bQ8/YSK452/NNd3S/+BM2JW3ILeJf6/byjQsH079HF648qw/lVbUs3pnPuhzPqJxxA72JwTsefnl265qTtu4vYVTf7lw8PIlzB/XkiY8zT2n/4DV7Crj9+ZX07OrpAH1vk/8aRs8t2cXLK3P5yrMrmPHkUma/sIaDJeU8deu53HnBEJzDr0mruKyKb/1jDb3jo3nylgmN+gxaY1hyN/YVl1Na7j8yaefBUv733zuYflYf/jhrPLFRpz7mJiYynAmDevDe5gNkHjrKeYOb7l9oC3ExkXzroqF8uiOfVbtb3+90MpQYRE5SVU0tj7ydwfsZB7jlvBRiIsMBGN2vO9n5Rzle2fQH8P99uIPEblHcdfEwwNM80j0mgn9vOcj63CLiu0QyONGzc9fZA+KJjQqvH29/MorKKtlbdJzR/bpjZvzgyuEcKCnn5ZWtG067LqeQ255bSVJcNK/95xTOHdTTb/+C6ppa5m3Yx6Ujk/nF9WMpOV7Fyt0F/OTa0Ywb2IOUhFjOG9yL19fk1c8s/uV7Wzl8tJI/feVcenaNalVcDaV5h55m5X8+Aqi6ppb7X9tA1+hwfnb9mCY7t1tr6tBE8go9taBJQehf8HX75FQSu0XzPwu2t8sGRkoMIifhUGk5X5mzgr8s3c3Xpw7m+5cPr3/trH7dqXWw7UDjWsPeouN8uiOfWyal0C3a8401MjyMy0b15qOtB1mzp5BzBvao/+CKDA9j0uBerepnqKu1jPa2eU8Zmkj6oJ78ZeluaptZwrmyupaF2w+xKa/Y77hzjp/OyyC+SyRzZ59Pn/gYrhrTh637S+qHYC7NOkJ+aQU3pw/glvNS+OgHF/PhfdO4bXJq/XluOHcAuw4fY21OISuyjzB3VS7fuGAwY/rHn/T9NSetd93IpM+bk55ZnM2GvGIenTGmfqmNtlLXQR0TGcbYNryPpnSJCufuS4ayYlfBKa/pFAglBpEAlVfVcP2Ty9i4t4jfzxzHQ/8xmkifJpC6D+KmmpNeX52Hc58v41DnytG9KSyrYueho4z3NiPVmTI0gaz8Y/UbyQdq637PB6NvZ+jtU1LJKShj4Y5DfmU37y3m/tc2kP6zD7jjL6u49bkVfkNol2YeYWNeMd+9LK1+SYfpY/oA8J631vDmur10j4moH+UTHmYM8w4frXP12L50iQznpRW5/OhfmxjQswv3XJ52UvfVkoE9uxAVHkZmvqcDel1OIb/7YCdXjenDtWc337HdWmcPiKdrVDgTUnq2avLayZp1Xgq3np9C/x5dgn4tJQaRAK3PLWJv0XEev2kcM8b1b/T6gJ5d6B4T0agDurbW8dqaXKYOS2g0omba8KT6D5VxKf6JoW65hkU78k8qzi37SkiOiybJZ2OY6WP6kBwXzV+XfT6i6FBpObfMWc6CjANcMboPv/zSWErLq/jdh593cj75SSa9u0fzpQmf3++AnrGcMyCe9zfv51hFNe9vPsA1Z/cjOiK82Zi6RUdw1dg+/HNtHtn5x/jZF8e0SVu/r4jwMIYkdWXLvhL+Z8F2bnz6M3p1jeKxL7ZtE1KdyPAwfjdzPD+6quWO87YQHRHOz744llRvc2MwKTGIBGh59hHMaHYmsZkxul/3RkNWl2UdIa/weP2QVl9doyO40JsAxg3wTwyj+3YnKS6aRTtPrjlpi7fj2VdkeBhfOW8Qi3bkk+39Rv3IvC2UV9fy1nem8vjN5zBrUgq3nJfCP1bksONgKetyCvks+wjfvHBIow/96WP6siGvmL8u283xqpr6pa9P5MZzPctR/8c5/bh4RHILpVtnWHI3Fu88zBOfZPLFcf15//sXtnkTkq8rRvdm7IDgNiOFghKDSIBWZBcwum934rs0vzTB6L7xbDtQ4rcd4yurc4nvEskXzurT5HvuuTyNB68e1agTtm7Z5iU78xtt79jcCKPK6loyD5Uyul/jMfWzzhtIZLjxwmd7+HDLQd7dtJ/vXTqMIT7r+9x7+XBio8J57J0tPLUwix6xkcyalNLoXFd5m5N+9+EOBvTsQvqgno3KNDR5SAK/nzmOn31xTItlW+vyUb0Z0TuOv3xtIo/ffA49YtumY7uzCXpiMLPpZrbdzDLN7IEmXr/DzPLNbL33zzeCHZM01h4jHTqaHQdLeeGz3VQFsGRERXUNa3MKWxyWeFa/7pRX1bLrsOdbeVFZJQsyDvDFcf3qRy81dPaAHs0umXDR8CQKy6rYvPfzTuE1ewoY+/ACVjQxlHXnoVKqalx9f4ev5LgYrhnbl9fX5PGTtzYzonccs6cN9SuT0C2aey5LY/FOzwJyd0xJpWt04yaf1MSujOrbnaoax/Xj+ze5RlBDZsaMcf1PmFhP1RfH92fBvdO4JEg1ks4iqInBzMKBJ4GrgNHALDMb3UTRV5xz47x/ng1mTNJYVU0ts+Ys5+F5Ge1yvbLKan7z/jbyCtt/1UjwjPH/9otr+MLvFvHQWxk89NbmFhPjhtxiKqprOX/IiYcl1n1Tz9hXgnOOV1blUlldy81NLOMQiLplmz/16Wf486fZVNU4nluyq1H5umaspmoMALdNSeVoRTUHSsr51Q1jm+w0vW1yKkMSuxIbFc4dU1Kbje3as/ti5vkwljNLsFdXnQRkOueyAcxsLjAD2BLk655RXlyxh5+/u5WYyHB6dIlkaHI3Hr/5nFavvNnQX5buYnm2Z2r/TekDOKtfcNtMn1u8i6cWZrE8+wiv/ufkgCY3Zecf5VhFDcP7dGu2k/OJj3eSlX+MX90wtsky1TW1/Pr9bcxZvIu46Ai+c/EwjlfV8NySXQxJ7HbChc5WePsXWhqvPiy5G1HhYfzfBzt47J2tHD5awbiBPVr9M03oFs3Y/vEs2pHP9y5LI7egjA+2HqRX1yg+3HqQfUXH6eczSmXr/lK6RIaTmtB0B+X4gT24akwf0pK7MT6l6eafqIgwnr09ncKyyhM2xXzjwsFcNDypfnE5OXMEOzH0B3zn4ucB5zVR7gYzmwbsAO51zp36/P0zxNqcQn76VgbjU3owsk93Co5V8u6m/Ty/ZJffGPrWyiss4/8+2MmFaYlszCvm1+9v54WvT2qDyJtWXFbFM4uzSekVy9qcIp5amMX3LjvxsMXisiq++ORSSsqriQw3RvSJ44YJA7hjSmr9aJOXV+bwP94lAyqra/nDrPF+SxQUl1Vx98trWbzzMF89fxD3XzmC+NhIamsd+4uP84v3ttKvRxfCDP65No+lmUeYc1t6fUfz8l1HGNmne4tt1pHhYVw6MplNe4u5MC2Riam96od3tta0tCT+9GkWJeVV/G3ZbsLMeOar53LTnz/j5ZU5/ODKEfVlt+wvZmTfuGaXZzCzEy4rUWdIAB/20RHhbToPQTqOjrAfw9vAy865CjP7T+BvwKUNC5nZbGA2QEpK486wUCivquGz7CN8ut2z3k1khBEZHsbE1F5cMbq33xj3Okcrqvnnmjx6xEY2OeTRV+GxSu5+cS19e8Tw7G0T69djr3phNc8t2cXXpg6ub691zrF1fyldo8MZ0DM2oHVbnHP89K0MzOBXN5zN/I37+fn8rSzNPFw/VPJkHT5awfxN++keE0mf+BgGJcTSN/7zb7R/XpTF0YpqXvvWZJ5emMXvP9rJtOFJ9UtBNOXpRVmUVlTz2Iyz2FtUzopdR3jk7S2syyni1zeczbqcQn7y5mYuGp7ElKEJ/PK9bfSIjeRnXxxDZU0tyzKP8Og7W8grLONXXxrLTJ/O1LAw4/GbxrG38DO+89JaAJLjoukWE8Gj72Qw/3sXUuu86wRNDOz37umvtvzBezKmDU/iiU8y+SDjIK+szuXqsX1JT+3FZSOTeXllLt+9NI2oiDCcc2zZV8K15zS/SqlIIIKdGPYCvo2rA7zH6jnnfHvQngV+09SJnHPPAM8ApKenh7SntLyqhh+9sYkFGQcoq6yhS2Q48V0iqaqppayyhr8s3U1SXDQzJw5kwiDPYlthZny87SBzV+ZS6t3o/fDRSr+NWnzV1jrue3U9h49W8vpdk/026bjn8jT+/YeDPL9kF/de4ak1PLUwq36BtuiIMIb3juPBa0bVbyTSkHOOdzft56Nth/jxNaPo36MLX508iL8u282v3tvGW9+ZGlCHoq/cgjJufW5Fo83Lv3fpMO65fDgFxyr5y9LdXHt2P0b26c4jM8awanch35+7jne/d2GTnZyHSsv5y9JdXHdOP77qnUnrnONPn3rud+eho+wrOs7gxK788ZbxdI+JpLCsiqc/zWLHwVK27i/laEU1id2ieemb5zfaIQw8s0rn3J7Os4t3MWVoAhcMS+SDLQe568W1vLo6jxF9ulFeVdvszzLYxqf0oFt0BD97dwul5dV8bWoqALeeP4gPt67i/YwDXD2mDw+/nUFJeTXnNtNEJBKoYCeGVUCamQ3GkxBmArf4FjCzvs65/d6n1wFbgxXMsYpqSsur/faKbaiuqSa/pJzDxypxDr5/eZrfRt6/eX87/1q3l1mTUvjCWb05f0hC/YiTmlrHpzsO8eLyHJ74JBPfPs3wMOOasX25fUoqzy7O5rF3tlBTW8vsaUM5WlHNkp2HWZtTSHb+UTIPHWX3kTIem3EWZzcY335Wv3i+cFZvnl+6i69fMJglOw/z2wWelSOnDU9k58GjfLD1IF/7yyqev2Mik70bkmw7UMJDb2aQmX+U4uNV1NR6Rq/UdTDGRIZz3xXD+cFrG3hxxR6+PDEl4BmdmYdKufXZlZRVVvPyN88nKS6agyXlvLF2L3/4OJP1ecX0joumorqGe70zXuO7RHrGz89Zzh8+3tnkRKEnPs6kusZxr0+zmZnx7YuHMapPd743dx2R4WE8f8fE+j6XH04fQWl5FQsyDnLN2L5MH9OHKcMSTjgBKzkuhv/ns8Ln9DF9SB/Uk//9YAc3p3vG3wd7PZzmRIaHMWVoAv/ecpBxA3swwfvBPy0tiUEJsTy3ZBevrc5l8c7DfOuioQHNKRA5EQv2MEUzuxr4HRAOPO+c+7mZPQqsds7NM7Nf4kkI1UABcJdzbtuJzpmenu5Wr1590rE89NZm3lq/j59fP4Zrz/avbpdVVvP8kl08/Wk2RyuqMYNesVGUVlQzNKkbr31rMt2iI1iefYRZc5bz1fMH8eiME4/HPlBczv7i41TVOCqraxmW3K0+KVXV1HLvK+t5Z+N+zhnYgy37iqmqcURFhDE4oStDkrpy/pAEbps8qMlZm1v2lXD1HxZz9dg+fLT1EGP6x/PiN86rT1D5pRXcMmc5eYXHee6OdLbsK+E372+ne5dIvnBWb3rGRtEjNpL/OKefX9KrqXXMeHIJm/eW0MW7L0ByXAz7i4+zv7ic+C6R3HNZGhePSMLMqK6pZUHGQX785iYiwsP4+52TGNnn8xExzjleXpnLw/MyqKyp5aZzB/Dbm87xu5f7XvX8HD7+wUUM6Pn5zODcgjIufXwhN6UP5BfXj23yZ3yopJwa5/yaq9rKupxCrn9qGeFhRlpyN97//rQ2v0agXlyxhwf/tZnfz/SfdT1nUTY/n7+ViDDj59ePabQXsogvM1vjnEtvsdzpOH69tYkhO/8o9766gQ25RVw/vj8/nD6SrftLWLzzMG9v3Ed+aQVXju7NfVcOJy3Z04G3cPsh7vzbaqYOS+SPM8dz7ROLCTPjvXsuPOUp/dU1tTzy9hZW7S5g2vAkLh2ZTPqgngEvQfytv6/h/YwDDOzVhTe/PZWEBjM865JD3eYlV4zuza++NLZRuYbqai/Ls4+wPPsIRWVV9O0RQ7/4LmzaW0xOQRmTBvdi6tBEXlmVw77icoYkdeX52yc2O11/Q24Rzy7Zxf+7emSjD/F9Rce55H8WctWYPvxu5vj64/e+sp53N+1n0X9dcsJaXjB99+V1vL1hH7dPHsQjLXwRCKbyqhre33yA/zinn3+n+vEqfvLmZmZOHNjirmMiSgzNqK6p5YlPMvnjx5n1s0mjIjxV9e9eOqzJjcTnrszhgTc20bt7NIdKK3jtPyeT3kRbdXvLyj/KY+9s4cGrR9WvLNlQfmkFP3lzM5eMTOLm9IGnvGZMZXUtc1fl8IePdnL4aCVThiZwx5RULhvV+5Q2KvnN+9t4amEW8+6eyph+8fx8/laeW7KLuy4eyg+njzylmE9FbkEZs+Ys57c3nlPfJCdyulJiaMGG3CI+3ZHPhJSepKf2bHZWap3H/72dP36cyexpQ/zaojursspqisqq/MbQn4rS8iou/u1ChiZ3I6lbNO9u2s8dU1L5ybWjg7Yzlkhno8TQxpxzrM0p5JwBPdpktylp7IXPdvOQd4/gH18zijsvGByUVTFFOqtAE0NHmMdwWjCzJpuZpO3MmpRIkjPKAAASyklEQVTC1v2lTEtL5KoTbAwvIsGlxCAdRmR4GL/8UtOjj0Sk/ahNRERE/CgxiIiIHyUGERHxo8QgIiJ+lBhERMSPEoOIiPhRYhARET9KDCIi4keJQURE/CgxiIiIHyUGERHxo8QgIiJ+lBhERMSPEoOIiPhRYhARET9KDCIi4ifoicHMppvZdjPLNLMHTlDuBjNzZtbitnMiIhI8QU0MZhYOPAlcBYwGZpnZ6CbKxQH3ACuCGY+IiLSsxcRgZuFm9kkrzz8JyHTOZTvnKoG5wIwmyj0G/Boob+V1RESkjbSYGJxzNUCtmcW34vz9gVyf53neY/XMbAIw0Dn3bivOLyIibSwiwHJHgU1m9gFwrO6gc+57p3JxMwsD/he4I4Cys4HZACkpKadyWREROYFAE8Mb3j8nay8w0Of5AO+xOnHAGGChmQH0AeaZ2XXOudW+J3LOPQM8A5Cenu5aEYuIiAQgoMTgnPubmUUBw72HtjvnqgJ46yogzcwG40kIM4FbfM5bDCTWPTezhcD9DZOCiIi0n4ASg5ldDPwN2A0YMNDMbnfOLTrR+5xz1WZ2N7AACAeed85lmNmjwGrn3LxTCV5ERNpeoE1JjwNXOue2A5jZcOBl4NyW3uicmw/Mb3DsoWbKXhxgPCIiEiSBzmOIrEsKAM65HUBkcEISEZFQCrTGsNrMngX+4X3+FUD9ACIiZ6BAE8NdwHeAuuGpi4GnghKRiIiEVIuJwbusxfPOua/gmXMgIiJnsEBnPg/yDlcVEZEzXKBNSdnAUjObh//MZ9UgRETOMIEmhizvnzA8s5VFROQMFWgfQ5xz7v52iEdEREIs0D6Gqe0Qi4iIdACBNiWt9/YvvIZ/H0NrFtYTEZEOLNDEEAMcAS71OeZo3YqrIiLSgQW6uurXgh2IiIh0DAGtlWRmw83sIzPb7H1+tpn9OLihiYhIKAS6iN4c4EdAFYBzbiOevRVEROQME2hiiHXOrWxwrLqtgxERkdALNDEcNrOheDqcMbMbgf1Bi0pEREIm0FFJ38Gz3/JIM9sL7MKz9LaIiJxhAh2VlA1cbmZdgTDnXKnv695tPv8WjABFRKR9BdqUBIBz7ljDpOB1TxvFIyIiIXZSieEErI3OIyIiIdZWicG10XlERCTEVGMQERE/bZUYljb3gplNN7PtZpZpZg808fq3zGyTma03syVmNrqNYhIRkVYIaFSSmUUDNwCpvu9xzj3q/fvuZt4XDjwJXAHkAavMbJ5zbotPsZecc097y1+HZ1/p6Sd9JyIi0iYCncfwFlAMrAEqTuL8k4BM73BXzGwuMAOoTwzOuRKf8l1Rf4WISEgFmhgGOOda8y2+P5Dr8zwPOK9hITP7DnAfEIX/0t6+ZWYDswFSUlJaEYqIiAQi0D6GZWY2NlhBOOeedM4NBX4INLlqq3PuGedcunMuPSkpKVihiIh0eoHWGC4A7jCzXXiakgxwzrmzW3jfXmCgz/MB3mPNmQv8KcCYREQkCAJNDFe18vyrgDQzG4wnIcwEbvEtYGZpzrmd3qfXADsREZGQOWFiMLPu3s7hppbBaJFzrtrM7gYWAOHA8865DDN7FFjtnJsH3G1ml+PZ66EQuL011xIRkbbRUo3hJeBaPKORHP4T2RwwpKULOOfmA/MbHHvI57HWWRIR6UBOmBicc9d6/x7cPuGIiEioBdrHgJn1BNKAmLpjzrlFwQhKRERCJ9CZz9/As7T2AGA9cD7wGc3MORARkdNXoPMY7gEmAnucc5cA44GioEUlIiIhE2hiKHfOlYNn3STn3DZgRPDCEhGRUAm0jyHPzHoAbwIfmFkhsCd4YYmISKgEuufz9d6HD5vZJ0A88H7QohIRkZBpMTF4l87OcM6NBHDOfRr0qEREJGRa7GNwztUA281MS5qKiHQCgfYx9AQyzGwlcKzuoHPuuqBEJSIiIRNoYojBszRGHQN+3fbhiIhIqAWaGCIa9i2YWZcgxCMiIiHW0uqqdwHfBoaY2Uafl+KApcEMTEREQiOQ1VXfA34JPOBzvNQ5VxC0qEREJGRaWl21GCgGZrVPOCIiEmqBLokhIiKdhBKDiIj4UWIQERE/SgwiIuJHiUFERPwoMYiIiJ+gJwYzm25m280s08weaOL1+8xsi5ltNLOPzGxQsGMSEZHmBTUxeJfsfhK4ChgNzDKz0Q2KrQPSnXNnA68DvwlmTCIicmLBrjFMAjKdc9nOuUpgLjDDt4Bz7hPnXJn36XJgQJBjEhGREwh2YugP5Po8z/Mea86deJbgEBGREAl0ddWgM7NbgXTgomZenw3MBkhJ0Z5BIiLBEuwaw15goM/zAd5jfszscuBB4DrnXEVTJ3LOPeOcS3fOpSclJQUlWBERCX5iWAWkmdlgM4sCZgLzfAuY2Xjgz3iSwqEgxyMiIi0IamJwzlUDdwMLgK3Aq865DDN71MzqtgX9LdANeM3M1pvZvGZOJyIi7SDofQzOufnA/AbHHvJ5fHmwYxARkcBp5rOIiPhRYhARET9KDCIi4keJQURE/CgxiIiIHyUGERHxo8QgIiJ+lBhERMSPEoOIiPhRYhARET9KDCIi4keJQURE/CgxiIiIHyUGERHxo8QgIiJ+lBhERMSPEoOIiPhRYhARET9KDCIi4keJQURE/CgxiIiIn6AnBjObbmbbzSzTzB5o4vVpZrbWzKrN7MZgxyMiIicW1MRgZuHAk8BVwGhglpmNblAsB7gDeCmYsYiISGAignz+SUCmcy4bwMzmAjOALXUFnHO7va/VBjkWEREJQLCbkvoDuT7P87zHRESkgzptOp/NbLaZrTaz1fn5+aEOR0TkjBXsxLAXGOjzfID32Elzzj3jnEt3zqUnJSW1KpjjmzZR9PrrrXqviEhnEezEsApIM7PBZhYFzATmBfmazSp5/30OPPIorqoqVCGIiHR4QU0Mzrlq4G5gAbAVeNU5l2Fmj5rZdQBmNtHM8oCbgD+bWUaw4okZOQpXVUVF9q5gXUJE5LQX7FFJOOfmA/MbHHvI5/EqPE1MQRczehQAFdu2EjNieHtcUkTktHPadD63hajUVCwmhvKt20IdiohIh9WpEoOFhxM9fDjl25QYRESa06kSA0DMyJFUbN2Kcy7UoYiIdEidLzGMGklNcTHV+/eHOhQRkQ6p0yWG6JEjAdScJCLSjE6XGGKGDwczyrduDXUoIiIdUqdLDGFduxI1aBAVqjGIiDSp0yUGgOhRIzVkVUSkGZ0yMcSMHEVVXh41JSWhDkVEpMPpnImhbgb09u0hjkREpOPpnImhbmSSmpNERBrplIkhIimJ8MREDVkVEWlCp0wM4Kk1aMiqiEhjnTcxjBpJRWYmNaWloQ5FRKRD6bSJofvVV0NVFUWvvBLqUEREOpROmxhiRo2i65QpFPztBWorK0MdjohIh9FpEwNArzu/TnV+PiVvvx3qUEREOoxOnRi6TplC9KhRHHnueVxtbajDERHpEDp1YjAzEu68k8rsbI4uXBjqcEREOoROnRgAuk//ApH9+nHk2ee0eY+ICEoMWEQEve78OsfXruXAw4/gqqtDHZKISEgFPTGY2XQz225mmWb2QBOvR5vZK97XV5hZarBjaqjnrFkkfPObFL3yCnl3f5fasrL2DkFEpMMIamIws3DgSeAqYDQwy8xGNyh2J1DonBsG/B/w62DG1BQLCyP5B/fR56cPcXTRIvbc+lWK33lXk99EpFOKCPL5JwGZzrlsADObC8wAtviUmQE87H38OvCEmZkLQYN/z1mziOjdm/0//Sn77r8fIiPpOmkSMWPHEJ2WRvSwYUQkJBDevTsWFdXe4YmItItgJ4b+QK7P8zzgvObKOOeqzawYSAAOBzm2JsVdeindLrqI4xs2UPrBhxxbspgjc5ZDTY1fOYuJISw6GouO9iSJ8DAMg7AwMAtF6B1XUz+Pts77ddc40wcQnMzvVjB/Fi3F0ZbXbu5aJ7rGyb7nZH9HT+X/eKA/m2au0efHD9J18uTWXz8AwU4MbcbMZgOzAVJSUoJ7rfBwYidMIHbCBPjhf1NbWUnlrl1UZmVRXVREbUkpNSUluIoKXGUFrrISV1MLtbU4dxrNh3BAW+Ww5s51ov8DzZU/2ZgaXuNkzhvI9dry53Qq523NZ+2p/Ixb82/a0rXb+lon8x7nPB+2Dd/Tmt/RUxXI71wzwrp2bYMATizYiWEvMNDn+QDvsabK5JlZBBAPHGl4IufcM8AzAOnp6e36tTAsKoqYESOIGTGiPS8rIhISwR6VtApIM7PBZhYFzATmNSgzD7jd+/hG4ONQ9C+IiIhHUGsM3j6Du4EFQDjwvHMuw8weBVY75+YBzwF/N7NMoABP8hARkRAJeh+Dc24+ML/BsYd8HpcDNwU7DhERCUynn/ksIiL+lBhERMSPEoOIiPhRYhARET9KDCIi4sdOxykDZpYP7DmJtyQSoiU2Qqwz3ndnvGfonPfdGe8ZTu2+BznnkloqdFomhpNlZqudc+mhjqO9dcb77oz3DJ3zvjvjPUP73LeakkRExI8Sg4iI+OksieGZUAcQIp3xvjvjPUPnvO/OeM/QDvfdKfoYREQkcJ2lxiAiIgE6oxKDmU03s+1mlmlmDzTxerSZveJ9fYWZpbZ/lG0rgHu+z8y2mNlGM/vIzAaFIs621tJ9+5S7wcycmZ32o1cCuWczu9n7751hZi+1d4zBEMDveIqZfWJm67y/51eHIs62ZGbPm9khM9vczOtmZn/w/kw2mtmENg3AOXdG/MGzrHcWMASIAjYAoxuU+TbwtPfxTOCVUMfdDvd8CRDrfXzX6X7Pgd63t1wcsAhYDqSHOu52+LdOA9YBPb3Pk0Mddzvd9zPAXd7Ho4HdoY67De57GjAB2NzM61cD7+HZC+58YEVbXv9MqjFMAjKdc9nOuUpgLjCjQZkZwN+8j18HLjM7rTdobvGenXOfOOfKvE+X49lF73QXyL81wGPAr4Hy9gwuSAK5528CTzrnCgGcc4faOcZgCOS+HdDd+zge2NeO8QWFc24Rnv1pmjMDeMF5LAd6mFnftrr+mZQY+gO5Ps/zvMeaLOOcqwaKgYR2iS44ArlnX3fi+ZZxumvxvr1V64HOuXfbM7AgCuTfejgw3MyWmtlyM5vebtEFTyD3/TBwq5nl4dn75bvtE1pInez//ZMS9I16pGMws1uBdOCiUMcSbGYWBvwvcEeIQ2lvEXiaky7GUzNcZGZjnXNFIY0q+GYBf3XOPW5mk/HsCDnGOVcb6sBOV2dSjWEvMNDn+QDvsSbLmFkEnmrnkXaJLjgCuWfM7HLgQeA651xFO8UWTC3ddxwwBlhoZrvxtMHOO807oAP5t84D5jnnqpxzu4AdeBLF6SyQ+74TeBXAOfcZEINnPaEzWUD/91vrTEoMq4A0MxtsZlF4OpfnNSgzD7jd+/hG4GPn7ck5TbV4z2Y2HvgznqRwJrQ5Qwv37Zwrds4lOudSnXOpePpWrnPOrQ5NuG0ikN/vN/HUFjCzRDxNS9ntGWQQBHLfOcBlAGY2Ck9iyG/XKNvfPOA27+ik84Fi59z+tjr5GdOU5JyrNrO7gQV4RjI875zLMLNHgdXOuXnAc3iqmZl4OnZmhi7iUxfgPf8W6Aa85u1nz3HOXReyoNtAgPd9RgnwnhcAV5rZFqAG+C/n3OlcIw70vn8AzDGze/F0RN9xmn/hw8xexpPkE719Jz8FIgGcc0/j6Uu5GsgEyoCvten1T/Ofn4iItLEzqSlJRETagBKDiIj4UWIQERE/SgwiIuJHiUFERPwoMYh4mdnRNjrPw2Z2fwDl/mpmN7bFNUXakhKDiIj4UWIQacDMunn3rlhrZpvMbIb3eKqZbfN+099hZi+a2eXeRet2mtkkn9OcY2afeY9/0/t+M7MnvHsLfAgk+1zzITNbZWabzeyZ03zVXznNKTGINFYOXO+cm4BnP4vHfT6ohwGPAyO9f24BLgDuB/6fzznOBi4FJgMPmVk/4HpgBJ49A24DpviUf8I5N9E5NwboAlwbpHsTadEZsySGSBsy4BdmNg2oxbOccW/va7ucc5sAzCwD+Mg558xsE5Dqc463nHPHgeNm9gmefQWmAS8752qAfWb2sU/5S8zsv4FYoBeQAbwdtDsUOQElBpHGvgIkAec656q8K7TGeF/zXZ221ud5Lf7/nxquNdPs2jNmFgM8hWeXuVwze9jneiLtTk1JIo3FA4e8SeESoDX7ZM8wsxgzS8CzGNoqPNuMftnMwr27bV3iLVuXBA6bWTc8K/+KhIxqDCKNvQi87W0eWg1sa8U5NgKf4NkX4DHn3D4z+xeefocteJaK/gzAOVdkZnOAzcABPElEJGS0uqqIiPhRU5KIiPhRYhARET9KDCIi4keJQURE/CgxiIiIHyUGERHxo8QgIiJ+lBhERMTP/wc4Ac9AtlyamQAAAABJRU5ErkJggg==\n",
"text/plain": [
"