{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 07 - Ensemble Methods\n", "\n", "by [Alejandro Correa Bahnsen](http://www.albahnsen.com/) & [Iván Torroledo](http://www.ivantorroledo.com/)\n", "\n", "version 1.3, June 2018\n", "\n", "## Part of the class [Applied Deep Learning](https://github.com/albahnsen/AppliedDeepLearningClass)\n", "\n", "\n", "This notebook is licensed under a [Creative Commons Attribution-ShareAlike 3.0 Unported License](http://creativecommons.org/licenses/by-sa/3.0/deed.en_US). Special thanks goes to [Kevin Markham](https://github.com/justmarkham)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why are we learning about ensembling?\n", "\n", "- Very popular method for improving the predictive performance of machine learning models\n", "- Provides a foundation for understanding more sophisticated models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lesson objectives\n", "\n", "Students will be able to:\n", "\n", "- Define ensembling and its requirements\n", "- Identify the two basic methods of ensembling\n", "- Decide whether manual ensembling is a useful approach for a given problem\n", "- Explain bagging and how it can be applied to decision trees\n", "- Explain how out-of-bag error and feature importances are calculated from bagged trees\n", "- Explain the difference between bagged trees and Random Forests\n", "- Build and tune a Random Forest model in scikit-learn\n", "- Decide whether a decision tree or a Random Forest is a better model for a given problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 1: Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ensemble learning is a widely studied topic in the machine learning community. The main idea behind \n", "the ensemble methodology is to combine several individual base classifiers in order to have a \n", "classifier that outperforms each of them.\n", "\n", "Nowadays, ensemble methods are one \n", "of the most popular and well studied machine learning techniques, and it can be \n", "noted that since 2009 all the first-place and second-place winners of the KDD-Cup https://www.sigkdd.org/kddcup/ used ensemble methods. The core \n", "principle in ensemble learning, is to induce random perturbations into the learning procedure in \n", "order to produce several different base classifiers from a single training set, then combining the \n", "base classifiers in order to make the final prediction. In order to induce the random permutations \n", "and therefore create the different base classifiers, several methods have been proposed, in \n", "particular: \n", "* bagging\n", "* pasting\n", "* random forests \n", "* random patches \n", "\n", "Finally, after the base classifiers \n", "are trained, they are typically combined using either:\n", "* majority voting\n", "* weighted voting \n", "* stacking\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are three main reasons regarding why ensemble \n", "methods perform better than single models: statistical, computational and representational . First, from a statistical point of view, when the learning set is too \n", "small, an algorithm can find several good models within the search space, that arise to the same \n", "performance on the training set $\\mathcal{S}$. Nevertheless, without a validation set, there is \n", "a risk of choosing the wrong model. The second reason is computational; in general, algorithms \n", "rely on some local search optimization and may get stuck in a local optima. Then, an ensemble may \n", "solve this by focusing different algorithms to different spaces across the training set. The last \n", "reason is representational. In most cases, for a learning set of finite size, the true function \n", "$f$ cannot be represented by any of the candidate models. By combining several models in an \n", "ensemble, it may be possible to obtain a model with a larger coverage across the space of \n", "representable functions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![s](images/ch9_fig1.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example\n", "\n", "Let's pretend that instead of building a single model to solve a binary classification problem, you created **five independent models**, and each model was correct about 70% of the time. If you combined these models into an \"ensemble\" and used their majority vote as a prediction, how often would the ensemble be correct?" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 0 1 1]\n", "[1 1 1 1 1 1 1 0 1 0 0 0 1 1 1 0 1 0 0 0]\n", "[1 1 1 1 0 1 1 0 0 1 1 1 1 1 1 1 1 0 1 1]\n", "[1 1 0 0 0 0 1 1 0 1 1 1 1 1 1 0 1 1 1 0]\n", "[0 0 1 0 0 0 1 0 1 0 0 0 1 1 1 1 1 1 1 1]\n" ] } ], "source": [ "import numpy as np\n", "\n", "# set a seed for reproducibility\n", "np.random.seed(1234)\n", "\n", "# generate 1000 random numbers (between 0 and 1) for each model, representing 1000 observations\n", "mod1 = np.random.rand(1000)\n", "mod2 = np.random.rand(1000)\n", "mod3 = np.random.rand(1000)\n", "mod4 = np.random.rand(1000)\n", "mod5 = np.random.rand(1000)\n", "\n", "# each model independently predicts 1 (the \"correct response\") if random number was at least 0.3\n", "preds1 = np.where(mod1 > 0.3, 1, 0)\n", "preds2 = np.where(mod2 > 0.3, 1, 0)\n", "preds3 = np.where(mod3 > 0.3, 1, 0)\n", "preds4 = np.where(mod4 > 0.3, 1, 0)\n", "preds5 = np.where(mod5 > 0.3, 1, 0)\n", "\n", "# print the first 20 predictions from each model\n", "print(preds1[:20])\n", "print(preds2[:20])\n", "print(preds3[:20])\n", "print(preds4[:20])\n", "print(preds5[:20])" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1 1 1 1 0 0 1 0 1 1 1 1 1 1 1 1 1 0 1 1]\n" ] } ], "source": [ "# average the predictions and then round to 0 or 1\n", "ensemble_preds = np.round((preds1 + preds2 + preds3 + preds4 + preds5)/5.0).astype(int)\n", "\n", "# print the ensemble's first 20 predictions\n", "print(ensemble_preds[:20])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.713\n", "0.665\n", "0.717\n", "0.712\n", "0.687\n" ] } ], "source": [ "# how accurate was each individual model?\n", "print(preds1.mean())\n", "print(preds2.mean())\n", "print(preds3.mean())\n", "print(preds4.mean())\n", "print(preds5.mean())" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.841\n" ] } ], "source": [ "# how accurate was the ensemble?\n", "print(ensemble_preds.mean())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** As you add more models to the voting process, the probability of error decreases, which is known as [Condorcet's Jury Theorem](http://en.wikipedia.org/wiki/Condorcet%27s_jury_theorem)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is ensembling?\n", "\n", "**Ensemble learning (or \"ensembling\")** is the process of combining several predictive models in order to produce a combined model that is more accurate than any individual model.\n", "\n", "- **Regression:** take the average of the predictions\n", "- **Classification:** take a vote and use the most common prediction, or take the average of the predicted probabilities\n", "\n", "For ensembling to work well, the models must have the following characteristics:\n", "\n", "- **Accurate:** they outperform the null model\n", "- **Independent:** their predictions are generated using different processes\n", "\n", "**The big idea:** If you have a collection of individually imperfect (and independent) models, the \"one-off\" mistakes made by each model are probably not going to be made by the rest of the models, and thus the mistakes will be discarded when averaging the models.\n", "\n", "There are two basic **methods for ensembling:**\n", "\n", "- Manually ensemble your individual models\n", "- Use a model that ensembles for you" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Theoretical performance of an ensemble\n", " If we assume that each one of the $T$ base classifiers has a probability $\\rho$ of \n", " being correct, the probability of an ensemble making the correct decision, assuming independence, \n", " denoted by $P_c$, can be calculated using the binomial distribution\n", "\n", "$$P_c = \\sum_{j>T/2}^{T} {{T}\\choose{j}} \\rho^j(1-\\rho)^{T-j}.$$\n", "\n", " Furthermore, as shown, if $T\\ge3$ then:\n", "\n", "$$\n", " \\lim_{T \\to \\infty} P_c= \\begin{cases} \n", " 1 &\\mbox{if } \\rho>0.5 \\\\ \n", " 0 &\\mbox{if } \\rho<0.5 \\\\ \n", " 0.5 &\\mbox{if } \\rho=0.5 ,\n", " \\end{cases}\n", "$$\n", "\tleading to the conclusion that \n", "$$\n", " \\rho \\ge 0.5 \\quad \\text{and} \\quad T\\ge3 \\quad \\Rightarrow \\quad P_c\\ge \\rho.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 2: Manual ensembling\n", "\n", "What makes a good manual ensemble?\n", "\n", "- Different types of **models**\n", "- Different combinations of **features**\n", "- Different **tuning parameters**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![Machine learning flowchart](https://raw.githubusercontent.com/justmarkham/DAT8/master/notebooks/images/crowdflower_ensembling.jpg)\n", "\n", "*Machine learning flowchart created by the [winner](https://github.com/ChenglongChen/Kaggle_CrowdFlower) of Kaggle's [CrowdFlower competition](https://www.kaggle.com/c/crowdflower-search-relevance)*" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# read in and prepare the vehicle training data\n", "import zipfile\n", "import pandas as pd\n", "with zipfile.ZipFile('../datasets/vehicles_train.csv.zip', 'r') as z:\n", " f = z.open('vehicles_train.csv')\n", " train = pd.io.parsers.read_table(f, index_col=False, sep=',')\n", "with zipfile.ZipFile('../datasets/vehicles_test.csv.zip', 'r') as z:\n", " f = z.open('vehicles_test.csv')\n", " test = pd.io.parsers.read_table(f, index_col=False, sep=',')\n", "\n", "train['vtype'] = train.vtype.map({'car':0, 'truck':1})\n", "# read in and prepare the vehicle testing data\n", "test['vtype'] = test.vtype.map({'car':0, 'truck':1})" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
priceyearmilesdoorsvtype
02200020121300020
11400020103000020
21300020107350040
3950020097800040
4900020074700040
\n", "
" ], "text/plain": [ " price year miles doors vtype\n", "0 22000 2012 13000 2 0\n", "1 14000 2010 30000 2 0\n", "2 13000 2010 73500 4 0\n", "3 9500 2009 78000 4 0\n", "4 9000 2007 47000 4 0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "train.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Train different models" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from sklearn.linear_model import LinearRegression\n", "from sklearn.tree import DecisionTreeRegressor\n", "from sklearn.naive_bayes import GaussianNB\n", "from sklearn.neighbors import KNeighborsRegressor\n", "\n", "models = {'lr': LinearRegression(),\n", " 'dt': DecisionTreeRegressor(),\n", " 'nb': GaussianNB(),\n", " 'nn': KNeighborsRegressor()}" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Train all the models\n", "X_train = train.iloc[:, 1:]\n", "X_test = test.iloc[:, 1:]\n", "y_train = train.price\n", "y_test = test.price\n", "\n", "for model in models.keys():\n", " models[model].fit(X_train, y_train)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# predict test for each model\n", "y_pred = pd.DataFrame(index=test.index, columns=models.keys())\n", "for model in models.keys():\n", " y_pred[model] = models[model].predict(X_test)\n", " " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "lr 2138.3579028745116\n", "dt 1414.213562373095\n", "nb 5477.2255750516615\n", "nn 1671.3268182295567\n" ] } ], "source": [ "# Evaluate each model\n", "from sklearn.metrics import mean_squared_error\n", "\n", "for model in models.keys():\n", " print(model,np.sqrt(mean_squared_error(y_pred[model], y_test)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Evaluate the error of the mean of the predictions" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1193.164765760328" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(mean_squared_error(y_pred.mean(axis=1), y_test))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing manual ensembling with a single model approach\n", "\n", "**Advantages of manual ensembling:**\n", "\n", "- Increases predictive accuracy\n", "- Easy to get started\n", "\n", "**Disadvantages of manual ensembling:**\n", "\n", "- Decreases interpretability\n", "- Takes longer to train\n", "- Takes longer to predict\n", "- More complex to automate and maintain\n", "- Small gains in accuracy may not be worth the added complexity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 3: Bagging\n", "\n", "The primary weakness of **decision trees** is that they don't tend to have the best predictive accuracy. This is partially due to **high variance**, meaning that different splits in the training data can lead to very different trees.\n", "\n", "**Bagging** is a general purpose procedure for reducing the variance of a machine learning method, but is particularly useful for decision trees. Bagging is short for **bootstrap aggregation**, meaning the aggregation of bootstrap samples.\n", "\n", "What is a **bootstrap sample**? A random sample with replacement:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]\n", "[ 6 12 13 9 10 12 6 16 1 17 2 13 8 14 7 19 6 19 12 11]\n" ] } ], "source": [ "# set a seed for reproducibility\n", "np.random.seed(1)\n", "\n", "# create an array of 1 through 20\n", "nums = np.arange(1, 21)\n", "print(nums)\n", "\n", "# sample that array 20 times with replacement\n", "print(np.random.choice(a=nums, size=20, replace=True))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**How does bagging work (for decision trees)?**\n", "\n", "1. Grow B trees using B bootstrap samples from the training data.\n", "2. Train each tree on its bootstrap sample and make predictions.\n", "3. Combine the predictions:\n", " - Average the predictions for **regression trees**\n", " - Take a vote for **classification trees**\n", "\n", "Notes:\n", "\n", "- **Each bootstrap sample** should be the same size as the original training set.\n", "- **B** should be a large enough value that the error seems to have \"stabilized\".\n", "- The trees are **grown deep** so that they have low bias/high variance.\n", "\n", "Bagging increases predictive accuracy by **reducing the variance**, similar to how cross-validation reduces the variance associated with train/test split (for estimating out-of-sample error) by splitting many times an averaging the results.\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array([13, 2, 12, 2, 6, 1, 3, 10, 11, 9, 6, 1, 0, 1]),\n", " array([ 9, 0, 0, 9, 3, 13, 4, 0, 0, 4, 1, 7, 3, 2]),\n", " array([ 4, 7, 2, 4, 8, 13, 0, 7, 9, 3, 12, 12, 4, 6]),\n", " array([ 1, 5, 6, 11, 2, 1, 12, 8, 3, 10, 5, 0, 11, 2]),\n", " array([10, 10, 6, 13, 2, 4, 11, 11, 13, 12, 4, 6, 13, 3]),\n", " array([10, 0, 6, 4, 7, 11, 6, 7, 1, 11, 10, 5, 7, 9]),\n", " array([ 2, 4, 8, 1, 12, 2, 1, 1, 3, 12, 5, 9, 0, 8]),\n", " array([11, 1, 6, 3, 3, 11, 5, 9, 7, 9, 2, 3, 11, 3]),\n", " array([ 3, 8, 6, 9, 7, 6, 3, 9, 6, 12, 6, 11, 6, 1]),\n", " array([13, 10, 3, 4, 3, 1, 13, 0, 5, 8, 13, 6, 11, 8])]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# set a seed for reproducibility\n", "np.random.seed(123)\n", "\n", "n_samples = train.shape[0]\n", "n_B = 10\n", "\n", "# create ten bootstrap samples (will be used to select rows from the DataFrame)\n", "samples = [np.random.choice(a=n_samples, size=n_samples, replace=True) for _ in range(1, n_B +1 )]\n", "samples" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
priceyearmilesdoorsvtype
131300199713800040
21300020107350040
121800199916300021
21300020107350040
63000200417700040
11400020103000020
3950020097800040
102500200319000021
11500020016200040
91900200316000040
63000200417700040
11400020103000020
02200020121300020
11400020103000020
\n", "
" ], "text/plain": [ " price year miles doors vtype\n", "13 1300 1997 138000 4 0\n", "2 13000 2010 73500 4 0\n", "12 1800 1999 163000 2 1\n", "2 13000 2010 73500 4 0\n", "6 3000 2004 177000 4 0\n", "1 14000 2010 30000 2 0\n", "3 9500 2009 78000 4 0\n", "10 2500 2003 190000 2 1\n", "11 5000 2001 62000 4 0\n", "9 1900 2003 160000 4 0\n", "6 3000 2004 177000 4 0\n", "1 14000 2010 30000 2 0\n", "0 22000 2012 13000 2 0\n", "1 14000 2010 30000 2 0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# show the rows for the first decision tree\n", "train.iloc[samples[0], :]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Build one tree for each sample" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from sklearn.tree import DecisionTreeRegressor\n", "\n", "# grow each tree deep\n", "treereg = DecisionTreeRegressor(max_depth=None, random_state=123)\n", "\n", "# DataFrame for storing predicted price from each tree\n", "y_pred = pd.DataFrame(index=test.index, columns=[list(range(n_B))])\n", "\n", "# grow one tree for each bootstrap sample and make predictions on testing data\n", "for i, sample in enumerate(samples):\n", " X_train = train.iloc[sample, 1:]\n", " y_train = train.iloc[sample, 0]\n", " treereg.fit(X_train, y_train)\n", " y_pred[i] = treereg.predict(X_test)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0123456789
01300.01300.03000.04000.01300.04000.04000.04000.03000.04000.0
15000.01300.03000.05000.05000.05000.04000.05000.05000.05000.0
214000.013000.013000.013000.013000.014000.013000.013000.09500.09000.0
\n", "
" ], "text/plain": [ " 0 1 2 3 4 5 6 7 \\\n", "0 1300.0 1300.0 3000.0 4000.0 1300.0 4000.0 4000.0 4000.0 \n", "1 5000.0 1300.0 3000.0 5000.0 5000.0 5000.0 4000.0 5000.0 \n", "2 14000.0 13000.0 13000.0 13000.0 13000.0 14000.0 13000.0 13000.0 \n", "\n", " 8 9 \n", "0 3000.0 4000.0 \n", "1 5000.0 5000.0 \n", "2 9500.0 9000.0 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_pred" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results of each tree" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 1621.7274740226856\n", "1 2942.7877939124323\n", "2 1825.7418583505537\n", "3 1000.0\n", "4 1276.7145334803704\n", "5 1414.213562373095\n", "6 1414.213562373095\n", "7 1000.0\n", "8 1554.5631755148024\n", "9 1914.854215512676\n" ] } ], "source": [ "for i in range(n_B):\n", " print(i, np.sqrt(mean_squared_error(y_pred[i], y_test)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results of the ensemble" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 2990.0\n", "1 4330.0\n", "2 12450.0\n", "dtype: float64" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_pred.mean(axis=1)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "998.5823284370031" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.sqrt(mean_squared_error(y_test, y_pred.mean(axis=1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bagged decision trees in scikit-learn (with B=500)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# define the training and testing sets\n", "X_train = train.iloc[:, 1:]\n", "y_train = train.iloc[:, 0]\n", "X_test = test.iloc[:, 1:]\n", "y_test = test.iloc[:, 0]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# instruct BaggingRegressor to use DecisionTreeRegressor as the \"base estimator\"\n", "from sklearn.ensemble import BaggingRegressor\n", "bagreg = BaggingRegressor(DecisionTreeRegressor(), n_estimators=500, \n", " bootstrap=True, oob_score=True, random_state=1)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 3344.2, 5395. , 12902. ])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# fit and predict\n", "bagreg.fit(X_train, y_train)\n", "y_pred = bagreg.predict(X_test)\n", "y_pred" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "657.8000304043775" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# calculate RMSE\n", "np.sqrt(mean_squared_error(y_test, y_pred))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating out-of-sample error\n", "\n", "For bagged models, out-of-sample error can be estimated without using **train/test split** or **cross-validation**!\n", "\n", "On average, each bagged tree uses about **two-thirds** of the observations. For each tree, the **remaining observations** are called \"out-of-bag\" observations." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([13, 2, 12, 2, 6, 1, 3, 10, 11, 9, 6, 1, 0, 1])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# show the first bootstrap sample\n", "samples[0]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{0, 1, 2, 3, 6, 9, 10, 11, 12, 13}\n", "{0, 1, 2, 3, 4, 7, 9, 13}\n", "{0, 2, 3, 4, 6, 7, 8, 9, 12, 13}\n", "{0, 1, 2, 3, 5, 6, 8, 10, 11, 12}\n", "{2, 3, 4, 6, 10, 11, 12, 13}\n", "{0, 1, 4, 5, 6, 7, 9, 10, 11}\n", "{0, 1, 2, 3, 4, 5, 8, 9, 12}\n", "{1, 2, 3, 5, 6, 7, 9, 11}\n", "{1, 3, 6, 7, 8, 9, 11, 12}\n", "{0, 1, 3, 4, 5, 6, 8, 10, 11, 13}\n" ] } ], "source": [ "# show the \"in-bag\" observations for each sample\n", "for sample in samples:\n", " print(set(sample))" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[4, 5, 7, 8]\n", "[5, 6, 8, 10, 11, 12]\n", "[1, 5, 10, 11]\n", "[4, 7, 9, 13]\n", "[0, 1, 5, 7, 8, 9]\n", "[2, 3, 8, 12, 13]\n", "[6, 7, 10, 11, 13]\n", "[0, 4, 8, 10, 12, 13]\n", "[0, 2, 4, 5, 10, 13]\n", "[2, 7, 9, 12]\n" ] } ], "source": [ "# show the \"out-of-bag\" observations for each sample\n", "for sample in samples:\n", " print(sorted(set(range(n_samples)) - set(sample)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How to calculate **\"out-of-bag error\":**\n", "\n", "1. For every observation in the training data, predict its response value using **only** the trees in which that observation was out-of-bag. Average those predictions (for regression) or take a vote (for classification).\n", "2. Compare all predictions to the actual response values in order to compute the out-of-bag error.\n", "\n", "When B is sufficiently large, the **out-of-bag error** is an accurate estimate of **out-of-sample error**." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7986955133989982" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute the out-of-bag R-squared score (not MSE, unfortunately!) for B=500\n", "bagreg.oob_score_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating feature importance\n", "\n", "Bagging increases **predictive accuracy**, but decreases **model interpretability** because it's no longer possible to visualize the tree to understand the importance of each feature.\n", "\n", "However, we can still obtain an overall summary of **feature importance** from bagged models:\n", "\n", "- **Bagged regression trees:** calculate the total amount that **MSE** is decreased due to splits over a given feature, averaged over all trees\n", "- **Bagged classification trees:** calculate the total amount that **Gini index** is decreased due to splits over a given feature, averaged over all trees" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 4: Random Forests\n", "\n", "Random Forests is a **slight variation of bagged trees** that has even better performance:\n", "\n", "- Exactly like bagging, we create an ensemble of decision trees using bootstrapped samples of the training set.\n", "- However, when building each tree, each time a split is considered, a **random sample of m features** is chosen as split candidates from the **full set of p features**. The split is only allowed to use **one of those m features**.\n", " - A new random sample of features is chosen for **every single tree at every single split**.\n", " - For **classification**, m is typically chosen to be the square root of p.\n", " - For **regression**, m is typically chosen to be somewhere between p/3 and p.\n", "\n", "What's the point?\n", "\n", "- Suppose there is **one very strong feature** in the data set. When using bagged trees, most of the trees will use that feature as the top split, resulting in an ensemble of similar trees that are **highly correlated**.\n", "- Averaging highly correlated quantities does not significantly reduce variance (which is the entire goal of bagging).\n", "- By randomly leaving out candidate features from each split, **Random Forests \"decorrelates\" the trees**, such that the averaging process can reduce the variance of the resulting model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 5: Building and tuning decision trees and Random Forests\n", "\n", "- Major League Baseball player data from 1986-87: [data](https://github.com/justmarkham/DAT8/blob/master/data/hitters.csv), [data dictionary](https://cran.r-project.org/web/packages/ISLR/ISLR.pdf) (page 7)\n", "- Each observation represents a player\n", "- **Goal:** Predict player salary" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
AtBatHitsHmRunRunsRBIWalksYearsCAtBatCHitsCHmRunCRunsCRBICWalksLeagueDivisionPutOutsAssistsErrorsSalaryNewLeague
131581724383914344983569321414375NW6324310475.0N
2479130186672763162445763224266263AW8808214480.0A
3496141206578371156281575225828838354NE200113500.0N
43218710394230239610112484633NE80540491.5N
55941694745135114408113319501336194AW28242125750.0A
\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 NewLeague \n", "1 414 375 N W 632 43 10 475.0 N \n", "2 266 263 A W 880 82 14 480.0 A \n", "3 838 354 N E 200 11 3 500.0 N \n", "4 46 33 N E 805 40 4 91.5 N \n", "5 336 194 A W 282 421 25 750.0 A " ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# read in the data\n", "with zipfile.ZipFile('../datasets/hitters.csv.zip', 'r') as z:\n", " f = z.open('hitters.csv')\n", " hitters = pd.read_csv(f, sep=',', index_col=False)\n", "\n", "# remove rows with missing values\n", "hitters.dropna(inplace=True)\n", "hitters.head()" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
AtBatHitsHmRunRunsRBIWalksYearsCAtBatCHitsCHmRunCRunsCRBICWalksLeagueDivisionPutOutsAssistsErrorsSalaryNewLeague
131581724383914344983569321414375006324310475.00
2479130186672763162445763224266263108808214480.01
349614120657837115628157522582883835401200113500.00
432187103942302396101124846330180540491.50
559416947451351144081133195013361941028242125750.01
\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 NewLeague \n", "1 414 375 0 0 632 43 10 475.0 0 \n", "2 266 263 1 0 880 82 14 480.0 1 \n", "3 838 354 0 1 200 11 3 500.0 0 \n", "4 46 33 0 1 805 40 4 91.5 0 \n", "5 336 194 1 0 282 421 25 750.0 1 " ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# encode categorical variables as integers\n", "hitters['League'] = pd.factorize(hitters.League)[0]\n", "hitters['Division'] = pd.factorize(hitters.Division)[0]\n", "hitters['NewLeague'] = pd.factorize(hitters.NewLeague)[0]\n", "hitters.head()" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# allow plots to appear in the notebook\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.style.use('fivethirtyeight')" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAADxCAYAAAB4begPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsnXl4VNX5xz/n3tmyJwQICWuAIIuALIoLiktFxQVtxX1vtXVfal1aWxeKWn/Wpa2t1q1WcV/R4oKKIgqySZA9gbAkgQDZM5PZ7j2/P2ZIAuZOQphkMsn5PM88mbk5995z5t55v/c95z3vEVJKFAqFQqGIB7RYV0ChUCgUitaiREuhUCgUcYMSLYVCoVDEDUq0FAqFQhE3KNFSKBQKRdygREuhUCgUcUOHiJYQor8QYr4QYp0QYo0Q4ubw9vuEECVCiJXh17Qm+9wthCgUQmwQQpzSEfVUKBQKRedGdMQ8LSFENpAtpVwhhEgBlgNnA+cBdVLKR/crPxJ4DTgCyAE+B4ZJKY12r6xCoVAoOi22jjiJlHIHsCP8vlYIsQ7oG2GX6cDrUkofUCSEKCQkYIv2FqiurlazohUKRdyRlpYmDmb/ttq+gz1vZ6HDx7SEEIOAccD34U03CCFWCSFeEEJkhLf1BbY32a2YyCKnUCgUim5Ah4qWECIZeAe4RUpZA/wLGAIcRsgT++veos3srjwrhUKh6OZ0SPcggBDCTkiwZksp3wWQUpY1+f+zwEfhj8VA/ya79wNKO6iqnZaCggLy8vJiXY0OQ7W3a9Pd2tse/CM9PeL/b6iq6qCadBwdFT0ogOeBdVLKx5psz25S7Bxgdfj9HOACIYRTCJEL5AFLOqKuCoVCES8ktPDqinSUp3UMcCnwoxBiZXjb74ELhRCHEer62wL8GkBKuUYI8SawFggC16vIQYVCodgXe6wrEAM6KnpwIc2PU82NsM8sYFa7VUqhUCjinA4b3+lEdMc2KxQKRZdAeVoKhUKhiBu6owHvjm2OO+wr38a+8i0GewOInk8gMwbEukoKhaIToDwtRafDtup9XHPuQPNU0AMwXjiXuus+g4TIoa4KhaLr01UjBCOhsrx3chwr30LzVDR81nZvxLZlcQxrpFAoOgv2Fl5dEeVpdXKkM3nfDfYEZFLP2FRGoVB0KrqjAVeeVifHe/osjOxDkZqNoD2RwKgzMfpPiHW1FApFJ0B5WopOh0zuSd21n6EX/8DWXVX0mzQNRJdI1qxQKA6S7mjAu2Ob4w9HIsbgY6g3CpRgKRSKBrpjIIYSLYVCoYhTumoXYCSUaCkUCkWc0h0NeHdss0KhUHQJlKelUCgUirihOxrw7thmhUKh6BIoT0uhUCgUcYOKHlS0iffft/H22w5SU01mzfKRkSFjXSWFQtENsHdDC94Nmxxd3njDzl13uaisDCUXWbNG55NP3CR0x0cghULRodhasOBd8fFZpXE6SN5+294gWADr1un8+KMewxopFIrugl2P/OqKKE/rIHE69/3sckmSk7vi841CoehstORp+TumGh2KEq2DZObMejZs0Cgo0ElMlEybFmTkSDPW1VIoFN0Au7PlMl0NJVoHSW6u5LPP3Hz7rU6vXpIjjjBiXSWFQtFd6IYWvBs2OfpkZEjOOCMY62ooFIruRje04N2wyQqFQtFF6IYWvBs2WaFQKLoIXTRCMBIq5L2DMQy47joXRx2VzHHHJfHhh+q5QaFQtBFbC68uSBdtVudl1iwnb77pIBgMLeZ4992CiRPdZGdbh8kHTdhUp7HLbyevoyqqUCg6Pyp6UNHerFmjNwgWQHGxxtq1GtnZzUcdVgfgFwuT2FijYZMjODcgeeQwb8Rz1NT4WbCghJQUB8cem4OmqdWOFYouSTe04N2wybFl8GATISRShoSkTx9JXp71vK57VrlYVrH3Mjl4bavJJYP8jElvfp9duzycffZHrF1bicOhcdJJ/Zk9+xQlXApFV6QbWnA1ptXB3H+/l2nTAgwcaDBsmMHdd3sZMMC6a7DCt+8lqg1qFHusL9u9937P2rWVAPj9Jl98sZ2vvy6OTuUVCkXnQm/h1QXphjodWxwOmD27HtMErRWPDKdmB+iTOptTBr6H33Dw7Oo/Mykz27K8xxPY57Pfb1JV1RWTuSgUiu5owbthkzsHrREsgIuHvc4vHPfgslUBcFLOBlz1c0BmNFv+yitHsnjxTsrK6gEYPjyDE07oF5U6KxSKTkY3tODdsMlxhv2jBsECcNoLkf4liOApzRY//vh+/OtfJ/DSS+twuWzcf/8k0tO7YYiRQtEd6IY/bSVanR0zdb/PSWD2ibjLiSf258QT+7djpRQKRaegG1pwFYjRyRG++yA4BqSLYCAFAtPRzLGxrpZCoegMtGFysRCivxBivhBinRBijRDi5vD2HkKIeUKIgvDfjPB2IYT4mxCiUAixSggxvsmxLg+XLxBCXN6ubQ3TDXU6vhAyA9xzQFtP6ZZyBg04MdZVUigUnYW2RQgGgd9KKVcIIVKA5UKIecAVwBdSyoeFEHcBdwF3AqcBeeHXJOBfwCQhRA/gXmAioUWSlwsh5kgpKw+uUZHpEE8rmsreFfjf/+q59NIKrr22kvLylpcyETgR5lgCPtXlFwuqCwtZeOWVfHP55exZsSLW1VEoGmmDpyWl3CGlXBF+XwusA/oC04GXwsVeAs4Ov58O/FeGWAykCyGygVOAeVLKirBQzQNOjXYT96ejPK2oKHsH1bVdmTOnnltvraK8PDQ3a/XqAJ9+2pPERNVT2xlxl5by9YwZ1BUVAVC+bBnHvf46PUaPjnHNFAoOOhBDCDEIGAd8D2RJKXdASNiEEL3DxfoC25vsVhzeZrW9XekQSxlFZY97Xn3V0yBYAGvWBFm2zHoelUTytusH/pL0GW8M38CPeklHVFMRpui11xoEC8BTUkLBs8/GsEYKRRMOImGuECIZeAe4RUpZE6loM9tkhO3tSoePaR2ksu9o7pgFBQXtVd2oEwg4afq1O52SqqoSCgqaT8u0rE8Z3yfvwNAlJMNrviUENg0jJeDooBrHllhf2+pgMDSpzmy8PnVStlu9Yt3ejqa7tDcvr51SXbfRggsh7IQEa7aU8t3w5jIhRHbYFmcDu8Lbi4GmYxP9gNLw9uP32/5V22rUejpUtPZXdiEs8+EdkIK32w3RDjz+eJBzzy1n40YDlwtOPTWBM8/Mweq7mJ9YFhKsMLXOAPrQNPKCXX/CcEFBQcyv7eBbb6X+m2/YvWgR0jTpMW4cxz34IPaUlKifqzO0tyPpbu1tF9pgwUXI2DwPrJNSPtbkX3OAy4GHw38/aLL9BiHE64SGaarDwvYp8ODeWARgKnB3W5pxIHSYaEVJ2eOeAQNsfPZZLxYs8NKzp85RRzksBQsg00gKXaVwkWTTQZ/9524p2g3d4eDE999n59dfYwYC9JkyBVtCQqyrpVCEaFv04DHApcCPQoiV4W2/JyRWbwohfglsA2aE/zcXmAYUAh7gSgApZYUQYiawNFzuASllRZtqdAB0iGhFS9k7oq4dQXq6xllnJbaq7Fm+MezS69ihVWEEDI6Xw8lSotWhaDYbOSedFOtqKBQ/pQ0WXEq5kOZ7swB+cqNLKSVwvcWxXgBeOPBatJ2O8rSiouzdETs613qOJYDB5oJNHJI3LNZVUigUnQWVxql9iKayd1fs6GiWX6FCoeiWdMP0EN2wyQqFQtFF6IYWvBs2WaFQKLoI3dCCd8MmKxQKRRehi65OHAklWl2QWlHFOscKnNLFKP/h2LDHukoKhaI96IYWvBs2uWtTqe1hbuJsavRKkLDZvp6z3Jehq0utUHQ9VPSgojPyyg923ltjJ+AdwuM9NIZkNp/yCWCZ86uQYAEI2KlvY5utkNzg8A6qrUKh6DC6oQXvhk2OL95YZee3nyTg8wsgg7NmGyy82k1GQvNZrWSz2a7aPYelQqGIBd3Qgqv1MDo5Tyx3hgUrREmlxsebrEdfx/uOJdlMa/jc2+hL/+DQdq2jQqGIEXoLry5IN9Tp+KJmfydJE+yR1s8aPc0+nOG+hNWOpTili3G+ySoQQ6HoqnRDC94NmxxfnDghyMt7NKgUoIO9n8nU3GDEfTLMXhzrndZBNVQoFDHDFesKdDyqe7CTc2O2j4QKoAqokPSvlQxLtA7EACis1bhzpYuH1jhxR9Y3hUIRz6juQUVn4+E5Luqrw2NaUlC8TWNxoc7Rw4xmy6+p0jhvkZ0STygWdt4uwdzjvLi66A2sUHRruqEFV55WJ8cf3DdJbiAIdX7rxLmzNvop8TT2GayosPNJWfMCp1Ao4hxbC68uiBKtTs4vp/jISm3sDjy0n8kxedZ9fnXa7n0+CyTVenG71U/RPBKJqaYaKNob1T2o6GycMMLguas8/GehA8NXzeOXO0mKMAv+4lHLWbxtIH6ZAFKSmVDOCVnejquwgpecu5lnr8YEDgsm8jtvNkItK6NoD7qhBe8WTa6uDpKfX0NWlpNDDkmKdXUOmGMPMTj2kHoKCraSkZQXseyGH87Av90OyYAh8BnJuPccAj07pq7dnXzdzVuOCuq0kHe8R6tmsOniXH+PGNdM0SVRaZy6HgUFbi66aBWbNnlITbVx8cXZzJrVdVf/LdidCEFbKNoQqCWR/F1uRvRUYYRtxvBg3/kGwvTj73M+2NMti67V6xsEC8AvYL1e3xG1VHRHurwF/yldfkzr978voKDAg2lCVVWQt94qo6Sk63aX5WYFEXqj0XQkGozvEzlEPt7ZKYLMdFXxZ1c1ZSLK4mx4SPrhTBI23EpCwZ0kr5iGCFRYFj8smESa2fizcknBYcHE6NZJodhLNwzE6KLNasTv39dg19cbVFcH6ds3RhVqZ3ZP2MFwbx27tvdB00yGHrGOTaljGRZhVPYLzeA1WxCbhLuDNgbKyCO4X+brPPuxE5su+eOFPob1i50o7hIGl7kqKHOEIiQXaD5m12fSq4U2tBZH6cvoNcsbRqR091ocRY/gG/Zws+VHmAlc4e3FR84qTCRHBJI5PWDtmSkUB0WXt+A/pcs3+eSTM1m+vIa6upBRGzo0kaFD4+vJ9/GvHHy2wU7AN4xHXILx/a1FIiV5GdOmb2z4LCXsrk0HObjZ8gs1gzvsfvaEnYP1WoB3fBqZFoEDi9bpXPv3RMqqQjus3qLz8Z/d9MmwjpRbYt/IcscmBDDFeygjjP4ttLj1/B1Pg2ABlNkN/lHv4X5SonMCs/4n34QwfRF3OSfQg3MC8TuG9UONxr2FLgKm4MTMIL/LjdxeRQzpohGCkejyonXDDQNxOjU+/7yctDQ7Dz88DIcjdr2ihiG5+eYSVqyox24X3HFHb04/PdWy/POL7fx1vpM6vwak8svXDT671k2v5OZFYpiAHU0+SwSHaQIspmq9rgcbBAugSJN8rhmcbzZ/a7z4maNBsACKynQ+WGTn19P8zZZfrxcz17UcjxYyfHsSaviVZyp9zAzLNh8Ia6t02O/rW1v9021tJZB9KY4dr6J7Qg8ChmsgvgE3RefgnZByv+Cq1YkU1YesYX6tTqrN5Nf9AzGumaJZurwF/yldfkwL4Oqr+/PGG4fx73+PokeP2CaPffDBXbz+ehVr1/rIz/dy55072LHD2iB8sdEeFqwQRRU6S7ZaP179wnckScFQlncpoX8wh0ONgZblU/bzI+wycqBhxn5iadclvdKsPb+VjqIGwQKo0etZY98W4QwHxhHlqfgrGkOofBVOjiyPkmIB0pGJe9yH+HKuxJd9Ce7D3kMm5kbt+J2NFTU6RfWN95vHFMyvUAmXOy3OFl5dkG6o07Fl9ep6gk1iBYqLA6xb5yU7u3nD0DtlX0FIc5kM7GEtEhkyjWs8v2CVfT0J0sWYwHC0CM8mdwXs/ChM1mgSB3CcoXGCaS2Kf7jQy9KNOvlFOnYdjhsd5OyjrIMfepkpCAkyrI12qdPLSLMsf6DcMjDI/OV9WdejFhCMrEjmxonRjdaTziy8wx+P6jE7K/1dJuk2SVWTTCy97F07kCeu6YYWvBs2OfoUFHh5440KcnLsXHZZT2w264mkubkOyHLC5N5QHSBrbTlDh1o/Ev15mpcNZTob92hoMsiFE0wOzY5sRKo8yXxVfDQZdsmI/n70CP3eaQje9jtZoZkkSsFYKSJOhE1NhLkz3SzdoONywoShBloEf/143xi26Lsp1ssRwLBAX0YHrT2/A8Wpw0ejd7Bn2YsA9Jp4JboW3cAHNyZvOtz4BZznTyIjwtIw8c7wZJMr+vp5Y4cdnynISzKZNazrRtvGPd3Qggsp4zPVTHV1daeo+NKldVx55RaKiwPoOhx/fApvvTUETWve8C/bA2d84cKb7ARDMixYz6JzAugREiaYJhRXC8qKN3H46OYDKvZS6NaYsSw0JiGQTEo3mHOEmxgO4yGRVAsPOoIU2fogmIKCAvLyIk+mFt5Kkt4/Hb1iLQBG5ijqpv8PXNERLg8mVyWVs84W8iYHGzrPuTPJjFJ0YlNa096OoiIgqAtCP5fE4lY+aDpTezuStLS0g/pGm9q+NH/k+7zaURW183YWuu4jYwfx6KNlFBeHxqQMAxYtcrNqlXX31GNbEkOCBaALilwJLK2MbAA1DQZkSNJdLSe+fajA2TCILhEsqdL5YndsH8cEgnSZdECC1VocK59qECwAvXwNzlX/jNrx33V4GgQLYLNu8G9nXdSO31npYZcMSGg/wVJEB6lHfnVFuqFzGV32d1RNU2Ka1k7g9qr9srYbUBvFiOL9Ow5NLAMH24zD93fsgTmAhtd5I4b9jCif4QCQzYynmdFrcXNHUiM8is6Cv6VFILtg0KcSrYPk+ut7sWqVh507gwgBEycmMWaMtUeRsFs2RvZIoA6Eu+XzVHnBZ7T82HvrIB+fvmjDs1UDDfofa3LSydHLEmH3v4PT939o1ACQ6N2KW8vD1A9ptrxE8pSrjCU2NwI4xZ/GRf7IiRD/+JKTT5fZCfhHcuFJGnecZ63q/rHXYS/6H3pVAQBG+jD8Y37TtsY1wzn+ROY46inUQ9/hAEPnKl/081d6DKiLEACjUDRHUG+hs0yJlmJ/pkxJ5dVXB/PKKxVkZ9u56abeEQMxBuiSJRuA3kAAUsol2SdYe2beIFz4fiLr9ugI81CuqJTceZS1Ef/0fTuBVQLCOuWdLyj5hcbQAdHxD2zBjxsEC0CTO7EFP8dvIVr/s1fxkaMKrwi18TVnOaOMBMYazRv+txbYeeFjB+5qDUjkqQ9MJg4LcuJhzXtPMrE37rM/wrHiSQD8429BJvY6iBbuSyoaz7szecFRh19IrvAl0UdG92fzxwIn75U58AVGMdELL4+ux6Y67hWtwLB1PxPe/VrcDowfn8T48a17+r7jRB8frbfj3QIC6JtlMipCbsA/fe1i/lYbodJOnl5hMn1YgOGZze/zY6FOoEm48u5KjWXroidapjYUiY4Id5yZJGNoIyzLr7J5GgQLoFYzybfVW4rW/JU67nUaeAAB1Ts1FqyyWYoWgEzMwjf5wbY1qBWkS43bfNGb+9WU7yp1/lPioNbQAAef7ZH8pcjkD0NUFgpFyxiRQoPpmkELXbFNnZpHFzrxCgF2gbQLSus11u6yvgzFtRo0CUGv9GpsqrQuP2qIgd3WKBK9MkwmjIjeKIzPeTsB21RMkYUhcgjYL8Kwn2hZfnQwEZdsrH+KqTEmmGBZvm67ALcMdZ2aQK3Eu6tTBIq2C2vqtLBghTAQbPKon6WidRjoEV9dEeVpdTDVXhHSoLAdr/MJdtUJRvZuvvzR/YLM32qjPuw99U81Gd/H2uv43eU+Fq/TWbLahq5Lrr/QR16UvCwAhI36xNcQshKJDUTkHH9nBNLZrPtYFh7TOjmQymEWXhZAD6eEfeaJCZLjLIRtq0cwc4OLgAk3DPZxeIb1939CD4Nsh8EOf8jApOgmP8tUy8goWkewBWHqirlMlGh1MNMOCbBou06NL/Q0PTTTZFyOtQhdP8FPeb1gwTYbAV89M08SZFvkHQT4ZqPOGqHjHhhSxjcKHFxZ7yfN2rlpE1K0LnegQHCzt0+rj3vxuX4+/tJOWdj77NvH5Nwz42c0ebdP8IslSRS6Q8ZkaZXO7AkexqU3L1xDk0wePaSeJ7e68Hi9nN3PxkU58dNeRWzxt5CrKco/+06BEq39MAzJrbduJz/fg9MpmDmzL5MmJVuWN03J7beXsnx5KAHufff1YfJka0/i8nEBAgbM3WjHqcODU+tJayFs1V8LRg0YfoGvhaGO5xfYqTq5noQRfjAFG+cmMH+tjbMnxMfT++GHmfx9loen/+uk3uPmnt/aGJ4XP0Hmb5fYGwQLoNSr8+8tDv51mHVWidN7G5ze2x2abJvb/SbbKtpOV+0CjESHdJ4LIV4QQuwSQqxusu0+IUSJEGJl+DWtyf/uFkIUCiE2CCFO6Yg67uW++0p59dVy8vPrWbLEw3XXbaWqytrgz5q1i5dfriQ/38uyZfXceGMJFRWRBeJXEwO8e5GH1873kBthSQ+Avy108PxSJytLbazZk8xtcxIprbHuLisb58d+Qj36IAN9cBDnRW78adGeqdV+SCmxz72NM/OHc87aMciP/xTrKh0QGQ4THYnmDKK7goAkRT0aKtqJto5pRcsmCyFODW8rFELc1W4NbUJHjfj+Bzi1me2PSykPC7/mAgghRgIXAKPC+/xTCHFQjxPV1QEWLtzNxo21LZb98cd9E9qWlATYvNnavcnPryfQpDdn2zY/GzZEL/Lr+202vE2iAYurNVYUW38dWRP8iCaem9bDRB8aWUS9Jiyq1llZp/1ksnRHs/bVV1kzezbVmzfj2b6dVc8/T9Gnn8a2UgfAz3MCDD+2lOzTisg+rYghJ23njkOim8BXodhLED3iKwL/4SBtctguPwWcBowELgyXbVc6RLSklAsA6zXK92U68LqU0ielLAIKgSPaeu7Cwjp+9rOvOOushUyd+jV/+MOPEcv36bPvY3Fmpo2+fR2W5XNy9h3q7NXLRv/+0Rv+HJhhImhUksxEk2G9rLvLDrEJtCbCky5hRIQEuNVBOG1VEmetTmLaj0lctj4hpsJVungxQY+n4bO/pobS77+PXYWA+iC8tN7Os2vtVDe/bFgDi+xegn1rsSUa6AkGZm8PbyV3/bRPithgYIv4siJKNvkIoFBKuVlK6QdeD5dtV2LdcXGDEOIyYBnwWyllJdAXWNykTHF4myUFBQWW/7v55i0UFISMRlVVgNdeK+K003SyspoXlt/8RlJQoLN9u4nDAZdcIqip2UJNTbPF+dWvJOvWCbZulTgccMEFkvr6rUSo0gFxyWBB/tYhbKpKQNck5+TtRlSVUVDVfPmfA8sGplOYYEdHMq3ci313HVbVebhmAIW6ZPSh6wkGbcwrHMV/1pQx2WnR4HbGMWIEta6hfOW9Cg2Tk5KeRxs6NOI1bk+8huA3qw5hTW0CIHhmlYd/jdlAur35LtdlvQXe5MYxTVPA2toyCrbubNX5YtXOWNFd2tteiYH9WD9Qt5EDtcnb99s+KdoV2p9Yita/gJmEZuTMBP4KXAXNugURn/0j3RB2exnQ+KQbCECPHv3Iy7OeLPrFF1BXZ5CQoKFHSr8eZt48cLtNXC7RqvIHytzh4PZ72VZUwIhD8oi0LO9urYTTE17Cp7vRMMlw9iE3/ZfYLG5uuc3HMUPmk5QY8m56Ze6CXZPJ650V9Xa0BpdrML955Ax2eEPRiWUJZ/Grqf3o2TM2c5eeWeNgTa2LvbdloTuRd6pH8LBFVpIZJe/wWd1gdiaHIibTvZWcV/oVeXnXtXiu7pb1vLu1tz1oKeT9ADlQm9zcj7Ld+2liNotRSlkmpTSklCbwLI1dgMVA/yZF+wGlbT3PySdnkZzcqM1DhqQwdKh1NKCUkpkzNzJjxjLOPXcZGze2rmsnKal1AtdWkhy0KrXPp9q3mLZq7CKILkyqbaWsE4WW5Qflrm0QLID0tEpysoqjUeU28fe/B9hR1RhOv21PBi++GDkE/N2v7Zx+RxKn35HE+wui+xxWH4T9f7PeCDkgh7jLeXDBvRxRupSJO5Zz29K/c3LZasvyCsXB0NbuweZog02Oqq1uLTHztIQQ2VLKHeGP5wB7f9lzgFeFEI8BOUAesKSt57nhhjwcDo3PPy8jPd3Oww+PwRFhcamHHirkiSc2Y4R7f37+82UsXjx5H+HrzGz12ejVJOLelILVHgejLSZs5NgFW5tukJChN+Yu3B+J5C8OP9/YgmjARX47FwYjd1E4v56JvfATJBr+sZcSmHiNZVlHM4dyRpiKsniNxl1Pu9hVFbqmhcUafXu7OXx4dMLkLz0kwKsFDjZWh55oB6UY3DzGOtAmMOg8pmx6gRM/vAwAI3EAnilvRqUuCsX+RDPkvQ02WQB5QohcoIRQsMZFUauQBR1iiYUQrwHHAz2FEMXAvcDxQojDCLmTW4BfA0gp1wgh3gTWEjKd10spDypm+5prhnDNNUNaVXb27JIGwQIoLvayfHk1U6ZkHkwV9iEYNCkqqiU11UFWVnSn/+3eORVc2+mVuAfDFGytHsLIwFDLWYbn+Iax2VZJme4GCUONDMYHrScDv2ELMNsRwBN2Nh53+hlj6oyyyFBuX/06jqVPowVCqey1hX/ByB6H2ffwZsvfdpuTL780WLs2JDpjx2pcdZW1KL63wNEgWABllRrvfW3n8OHRieDMdEk+nObmoRVOAibcNtZHbmqEHhBHOu4TPsT540NgePGPuBkzrflkwp0Zv6jAEB5cZjaiG84FihfaKlrRsslCiBuATwEdeEFKueZg2tMaOkS0pJQXNrP5+QjlZwGz2q9G1vj9PzVIHk/05jlVVfk599zPKSysxum0cf75uTzwwMSoHf+PWSlcuPYGEtOXEzAS6ekZx3G51t1rmTKRO91H8519OwnSxjGB/tgi9BovthkNggVQrsFi3bAULVvR/AbBAtDqy7EXfYXPQrR69BDMnZvI7Nl+Kir2cOutOSQnW3fH5fU3sOuSQLjLzmGT5PWP7mTkrETJE5Nbv+S8dPXCe/hjrS5vmvDtGp3ColRy+kNSS2sktTNbXc8Ebw2oAAAgAElEQVRRbv8aiR+X2Y9D3A9gI/rLsSgOnraKVrRscjgsfm6bKtFGVGbO/ZgyZd+1nlJTHUycGJ2l2wHuvnsJy5ZVU1WVQVlZAi+9VMiGDRahgG0g0yZ5b7DOLc6jeCBlLP/JDSBaGGpLk05O8w/l+MAg7C38CAYH9X2GWnUTRgetb6Ng9gSk3ti/ZzpTCfaNLNLp6YLrr3dywQWeiIIFcOFJAZITJaFKSZISJOefGD9pkAwDzp+VyC8eSOK2F/M47fdJVMUwQt6jbWGP/XOCWiWG5sZt28A217Oxq5AiIj6cEV+dESHETUKIyIvqRUCJ1n48+eRIRo9Ow+XSSUtzMHPmMHr1il5Y6ZYtQSAXyAJyqK7OYfv2yKtALlqic/WNCTz8xACqWqFvSTocn2owIcloUbAOlBK/nUDAhmkITENQ77NTbFg77IEJV+MfNQMjbSBGei7+8b/EyD0havV5eo6Dytq9GYgFlbWCf8+JehhwuzF3qY35K234gwKJYFWRjftfjq6rVY3kDzY/N9l9LBORew38WjmG2Fc1g1rLk/IVsSFOs7z/DNgihPhICHG+EOKA1DU+ogs6kM8+28WOHR68XgOv1+DVV4u5+OIcbBFC96666ge+/roch0PjlVfGMWGCdTLZmpoMGhdsF0ASQliPay34TufqGxPDCWR7U7g5yMfvukmMUSbMWinwehJodLcElSLCGI8QeE//B5gGCAEius9J5dVNUuaH61NRGz9Z4StqBUFz3/rW1kev/h4kFzp8rNFD12iR5ucffgeTZPMGLcnIw2nm4NNDQWCamUBaYELU6qOILp1YmCyRUp4lhMgkFLhxC/C0EOId4L/hSc8RUZ7Wfrz8cjF79jR2L+Xn17B2rXV/za9+9QPvvltGeXmQHTv8TJu2hN27rYMABgzYN9zeZtNISLDOoPHvF50NGc8BVq3R+XZx7G7UKwVksXf5EMFgJGe1xsZqetQFC+CqaQH6924cw+rf2+DyU1tIW9HOfBTUOKXewdR6B0/5I1+rMyYFyc5pDNVMyTC4Kor1/1YzWNMkRcouDV60W6f1sstUhnruJjUwjpTgaHJ8F5AVmGZZXhFbDiKNU0yRUpZLKZ+SUh4FTAEOB+YLIbYIIf4ghLCcl9QtPC0pJVVVQZKTdez2yIazqGjfAXe/XxIMWnsS8+eX7/PZ5zP58MOdXHXVwGbL33RTFvn59ZSWhoRx4sQkDj/cet6Y3b7vue02SEqM2IR25XANnjHhWSmxAb8X0KsVolUbBF1AYit+R/fmO/mk1I7fP5KL/Bq/G2X9EDC0n8nxfwzwXo0dJJyUHmRI38jzGx+pc/KuN/SgcJozwL0p0csVuc4Q3O63sVOG7rMNpqCvJjnb1nxwSE2qgf5HN9orCWCAPs1L6UiDaD1PJhD6kTeVKXsL0z+TzFyGe2ZG5fyK9uVA52J1JoQQJwGXEEr9tAx4BNgG3Ax8DBzb3H7x2+JWUlkZ4PzzV7Ntmw+XS+OWW/pzxRXZluV79EigqKjpGJOGlNZW+adzvgT9+ln33R11VDKvvz6YF18sp3dvnVtu6YPdbn38e+/0smq1zqYiHV03mTLZ4MjDY5u1/WgNjm5lWUPCL1cksKRSRxNwalaQRw+1jsR7a6ud5woduIMakMhTG0wmZgY5wWLhyw/rbLyv2alLC12HdzU7U90BpiU1X36eV+cpj4PqsKgUewRjbCbnJEQneGOeoTUIFkA1grlB3VK0vgbKs02cvwvdcz7gcwnnRqU2cLSpc4xp8K1mEhSQawruDHbFpQG7J/HYPSiEeJRQ12A18F/gHillSZP/LwYqrfbv8qJ1220FLFnSOJD86KPbOPPMnmRmNv/Dzc1NYvnyOsAABL17J5CTYz1O+PrrE5g6dTE+n4kQghEjEpk61WIZ4jBjxiTy+OOtc5cGDZR8+p6bOR/bCPhK+dUVvdDiqFP3qc0O/ldmJxAW/leLBT/rFeDUrOZFZUGZHhasEFUBjW922SxF65t6G7VNRKLG1FhQb7MUrW8CtgbBAqhDY4Ffj5pojdAkiUg84XE2G5KhmnUI/jAgSYI7/NyiSWjeR28bGoLnfQ7+UqFRFhTc0sOkn9KsLoMv+rkHOwIXcI6Ucmlz/5RSBoQQliHGcWT+2kbT8SkIeV6lpdbdQY8+OpQjjkijd+9EBgxI4qab+pOdbS1aY8em8fbbRzJlSj8uvngIX301OWp130vPTMlVlwQ48bgq9FY8WG0TAf7qrOJZZw31xHYBxbW1eoNgAXgMjR9rrJ+Vju5lkKg39l+l2U2O6mk9BnO0K0iyaGxjijA5xmXtiR5tD5LapHwSJkc7oue5nmwzuchmkCNMsjCZqpvcZpFcF+BoNC4CciRkBoOcANweISv/gWJKuCQ/iSfzU3h5dQoXLkumxBs/gSqKyEQzjVNHEF7O5DRgVaRyUsr1Vv/rfK2KMqNHJ7NoUXXDGlk5OU4GD7buvktPt/Ppp4exe3eAlBSdhITIKvHtt25+/esdlJYKNM3Htm1bee+9Qe2ahzAShZqfKwo9FH2ajnCafHFJBS85MnFGMITfVek8U+LApcG9uV5yXNHLeXlaVoCPy2xUh72n3g6Tqb2tvZoLcwOsqtL5rNRGMODnvCGCk3Osjf7ZKUFW+PzM9YTchzOSApyZbC1yp7oMrgn6+cBrRwKnOAPMiJKXtZebNIM9lXb8wG1pPmwt3AoPoHEHkvVbtjNhyGBEFEXr20qd+RW2hgeHjR6dewtdPHeoWuOrKxBv3YNSSkMIYRDytto0mNxq0RJCXAislFKuE0IcQiihYhC4LpIqxpoHHhiM222walUdCQkaDz88hKSkyBdaCEHv3q1zu598cjelpSEjaZqwZImH/Px6xo+PTbTEAxsDrPtbDkZlyIh/sSmBefdWc4aj+bk/i6s1frkugR3hKLcf6nQ+PcxNRkuj9a1kenaQbfU+3i+1IwRcM9DP2LTI3t9D47w8NK71WcAf6OnjgQO4/+9J9nFPcvSCL5qyOyj4+Y4kCgOh73Ol18bsbDfjnJHbnIwgzTSjKlgANUHB/klefKbytLoK8SZaYZ4A3hRCPEgo6W7DHSql3NzSzgfiaf2ZxvH3RwklTKwD/gmceADH6VBsNsGTTw5rt+Nv3blvD6vPL/D4YmcUVn+a2iBYAIFtLvJX1XOGRQ/xsyXOBsGC0JP4x+U2LuoTPe/jxsF+bhzc+jDuf2x28E6pHZ93OL92alw+IPYZLvYujNnSZO236+xhwQrtUGpo/LvKwb+yWp8GKppM6RFkVJLJGnfoGmc5TH7Vt30EW9HxxKlo/SP89+T9tktouUEHIlq9pJRlQggXMJlQgFMA2HMAx4gLTFOyY0eA1FSdlJTI32FKbl9Y54VA2CgnpSASkiCKY0lSwk6PoDbY8g060q7T9FFF2EwOd9oJBZb8lKbjRxAKHEizxW7p4rk7dR4tdFIV0IBkZq43OSTZ5MgesYuYnPU/J+8utyMlnDA8yKPneS3FK12T/CrrX0zr8SGaMFlWewS7qn/fsRVuQrIN5ox3c1+hkzpDcFU/P8dmxDb6VBE9OmuqpkhIKQ8qluJARGu3EGIoMBpYKqX0CSESaX6BsLilsjLIjBmb2LrVj8sluPrqXtx0k/WCiAOGJLNsyCioKAO7g/QhWfTpGb3xAncAzv8siQ1VGpo5kovqJPcebv2k/NfzAqwtsLN5rQ3NYXLCMQY/G2VtpO4Z7GGeu56dtelowuCQHlWcmhm7oc5PdtnDghViT0Bj3i5bRNEyJKysDAn6YRkG0RxOnL9e55mvHdTUh+q043uN0f0Mrjimee/vjPSlHJrwEkm20MrPWfadpCf2g+D50avUAZLpkPx9ZGw8PUX7Eqee1kFxINZpJrCc0CP73l/gSUB+tCsVS3772+0sW9a4KOI//7mLc8/NICen+TGuv9zopbA4kc3FA0l0wQVT/QzpFz0v649LXCzcufcyOXlxvcmMIQFG9mj+HA4JictAWytxOiEtU0bs0vo+OZ+TJm+ieFc/7PYAuZllbKk/jiFG5LD99mJsqoFTkw3jLomaZEyatWD5TfjFwkSWloe+o0mZQd6a7CHCkmkHxPebbQ2CBVAfECwp0i1FK2hb3SBYAE7dj2ZfZ7k+mUJxMMSjaAkhbMB1hDJh9KSJ4yOlPK6l/Vv905ZS/gfIBvpJKeeFN39Po4B1CSor9zWQVVUGO3daj6n0ypDM+4ebT//hZv4zdTzwm+iOF+z07HuJqvwaRTXWl+1P9zlZ/WMooW29R2Pux3a+/c76xq7Q3Lh3JFH6Xn9KP+xPXb2gTKuxLN/eXDUwwDnZfvq5DPrYfFzYz8/0bGuL/68CBwt32/CaAq8p+Ga3jWcLozd35dhcNz30xizFyVodxw+qtiyfEByHZjauCiCki4TgYVGrj0LRlDhN4/Q4obW6FgATgHeA3sCXrdn5QKIHf5BSjmu6TUq5SwixDIjeglAxZuLERL77rhZfWHv693eQlxc567bTAaMGt967+tG2iVWOQhzSxqneo0iR1pGGU3ICfFVqwxMMPYwMTDaY2Nva86ip2W9peJ+grMza1UrcPJjP/68fdbtTANi9qh93/8FLNLvKt+vFLHMuR0jBZN8x9DStF9QUAp4+zEu9AZsKCzn0kKERj73TG8qOvhcTwY4ozkM6wfs8d6fu5L91l2Oic0bCR1wa3IaPh5ot7zJHke69kjrnHCQmiYEjSQlMj1p9FIqmdMa5WK3g58BRUsptQoj7pZRPCiE+BZ4B7mtp5wNp8U+shxBCAIMP4Bidnrvvzqa+XrJ0aR1Op8bDD/drMRjjQMi3FXLP4z42zTsK3WGw8pHv+NPgKThpPk3BNSMDVPk0viy1EfB6ePBYQVaidaDE+ef5WfidjfLykDc2dIjBCcdbi9y894ZQt7vRM9ldlMl3S90MmByd/qwd2g4+TPgfdXoo6XCZvosL3eeRKlMj7pegg1NrOSDk0kEB5pQ4KAl34fVNMLl4YBSjDYP13J7+OLenP96wyRe4IuIu6YHzSA+cF706dDD5RRr3vurCHxT8bGyA286ObQJihTXx2D0IJALbw+/rhRCJUsr1QohxkXbaS4uiJYT4b/ito8n7vQwC2n155Y5E0wR//nPfA97P7zex2wWihZjovz0jWPzwRAJuGyCYfX4qv/h8F+PTms+HKATcOd7HneN9oXlLWZHnLU07zcDvr+e1NxzY7TDzvnp69LA2/vp+wiCQOFrxKBMwQBOgt9DBvNKxqkGwAKr1atba13Ok/4iWT9IKRqaZPHeEhyc3OpESfjvcy4gW5oEdCIFDL8Wx5lX0io0AGGkD8U28KWrH72yU1wiufDKRzWUhY5hfpJGcILnmlNhPO1D8FH98pnFaRyir+xJCiXLvE0LUACUR9wrTGk9rk8V7CXwLvNW6enZNPB6DSy5Zw8aN9TidGjfc0Jcrr8yxLL/+s94E3KuAWkCjqmAAm1emMn5K9Op09vQgZ09vnaf0h3N9LN1kY9NOHYHkiDyD0ydY72uacM3sBL4v0tE1OHdcgHtOtx7HS5D7dq0KKUiK0B3aFo7qaXBUT0/LBcNICcVVAlPCgIzIgSoyIRP3uXNwLnoIjAD+I25DZnSpzoV9WLFJZ3NZ45OI26fx5Sq7Eq1OSicet4rEzTTOwbkN+BeQAlzTmp1bFC0p5f0Qyrwrpfy0jZXsstx99ya+/LJxoP7//m8bU6dm0rdv84NC3s3lwO7GDWIzmebodq6lNT0zJZnXS0oWmdhcMPRYA2eEhKpPfungf2sl/YZuJBiw8ey3uZwwPMgxQ5rvgjzcexRfmKUkunYAAm/9QEYGRsZsooRpwhUvJfDtJhtSwhGDDGb/0hPRY5TJffCe/GTHVTKG9Otpkp4kqXI3XqCeqbHNX6mwJh7HtJomypVSFhBaybjVRGyxEOK4JitJBoQQzWa+kFK2KuojVqxbV8W7724lNzeZCy4YjKZFz2IWF+/rZZSV+dm61WspWlm96thU0GSD9COlh2hGPpSU1PHKKwWkpzu44orhOJ3WT2Mza50sSdDhxNB38p7h4HxfgOOczYvQunI/V/3uObIHlGCaGkXrh/LD5ks5Zkjzx/97ZQpPV/2SnKRtGFKj3D2Aw3p7OSclNjHgs5fY+WSNHb8Rau/n6wVPL3Bw/fFq3AZgRH+Ty07089ZCB74ADM02efBSNcersxIvY1pW2rE/rdGSlmT6n8Ch4ffPW52HThyM8dVXO7juusWUltbjcGjMnVvCyy8f2+LYU2s59NAkvv66siEhb//+LoYNs07IW17ugYaVfwFMpIxehoLCwirOOvszSotDy7G8/W4Rcz+aZrn45Q5To6nb40FQZGgcZ5FB4/BTPiE4sBgAXTcYMqKAkRn5wKhmy28NCgLSzta6RlXbHNSJ1cSlTbu1BsECCJqCTbu7/GIHB8QDF/m4+Qw/tV7o31O2OG6piB3xIlpY60dTWqUlEW9HKeWhTd7nWrw6rWABPPnkOkpLQxkq/H6Tb7/dRVFRXQt7tZ4//SmX8eN7k5ycSHp6MvfcM4SePa0HRwcN0ghlvzIIGe56bC2lAT8Abr9nZYNgASz9voz/fVpsWf4sV4CMJkt19NdMTnJaC0pu9r5jR7rNpFcP6+/z3OQAmU3Wk8rRDc5Mit34yDnjAvRp0t3VK9lkxgQ1XrM/mamSQb2VYHV24mWeVgT9OGAtaU304Dc0ycJrUaEWZzHHCrlfzU1TEgxGr4/+v/+VrF+fR13Ybj/2GJx2miQpqXkhOv/8viz4Zh31nlCXy+AhKUycmBa1+hQ0MydreRGcbVH+HFeQStPLu14HupDcneSln259uQ/1TaTEVoRHCwljmtGDIcHmvSyAkxINZmV6mV3rACS3ZfgY5ojdGMnYfiaPz/Dwz69D3bFXHe3nqMEqF58iPvHHYe7Bg6U1o3jPNXkvgKcIpeCICy65JJfvvivDH16fYdCgJIYOjTxHSEpJdXWApCSbZbfaXubOldQ0SSBRUACrV8OkSc2X3+oZSH2vRKjYBrqN6qwxeANBkiyeCwwDbrjHxfJVNkxzJLf/Bi6Ybu0ZjDx5AiUby8Ad9rZ6ZjFpcv+IbbgqMcAMRwC7AFcLD2d9zAH8zPNzfnR+jyZ1JnlPJFEmR9zngtQAF6R2Hm/mR83GnjwNKWCVXecclWNJEafEUfdgA0KIVEKTiJtL4zSgpf1bEz340n4nfGz/bZ2ZBQu24ffXAXbApLq6itpaP2lpzT+hVFcHuOCCRWzZ4sbp1Lj++jyuvtraa03Yb/gqKQkyMqzr89InDkgdEHoB5bWSL5Z7OO+k5g3nw085eesjB8GgABK5/zGToyYGGdi3eZF7/NoUppVPZ9vCVQink2N+Pp7Txlob5YAJly1PIL9GR0fw8xw/94+InIqqr5FLX09uxDKdlc936jxV6KA6nJS3xCMYm25wdj8lXIr4Ix5Fi1CsRD/gAeAV4BLgd4TSObVI/MVLHiBFRTWElgkJGeLduyXbttUxenTzonXHHfksWlTe8Pnxxzdw5pk59OnTfCqnv/xFo7DQpKAAkpPhF78QDBtmPUb1E39Di7x855oNeliwQuzYpfHjWp2BfZs3sv0yJPPvd/HBqqNJS5BMHxOMOC7x8EYnn+2yY4Qfdl7Y6uD0rCBHRHEpELcPPllpw2mDU8YGscfwrluw29YgWAC1QY2vd9mUaCniks40bnUATAVGSCnLhRCGlPKDcDrADwnlJYxIlxetAQP27brq1SuB/v2tu7PKyvYN762o8FNaWm8pWtnZgnnzNNasCXlYeXmRgyqmnRzg2Vc08AsQkNxbcvxYa4E4ZKjBZwtsDcKV1dNk1CGRx4QykyVXHd267riVlXqDYAHUGhqrqrWoiVaVB6b/JYn8bTY0ITlqmMF7t7tblXWjPTi6Z5D/bDapCYaEK0k3ObqnGtNSxCfxOE+LUADg3qzTdUKIdGAHzaQKbI7WBGLsH19vE0KcwL79kJ12ntbJJ+fy1ltF+P1BQJCenkxamnV039ix6SxcuIdgMNT91rdvAkOGRB6z+etfN/P661UkJgreeWcEAwdah7w/eLGPnVUa+Zt1XA645Swf/XtaBz784UYfW4s1Vq7WkdLPTb+U5A6IXiDDqt3avmuFCsmP23TIjc4Y1F8+cJK/LXSbmVLw3Qad1xbaufz42IxxnZptcPUQPx+U2JHAqX0CzOgEKyMrFG0hTtM45RMaz/oCWEgoTqIO2NianVsj0/vH15cDLzT53Knnaf3nP5vw+xu9pC1bfGzaVGsZjPGnP42itjbIDz9U4nRq/PnPo0lLs04RcffdBTz9tBspQ92NRx65lnXrxpCe3vw+dhu8fHPrF4m02eCFv4bKFxQUkJcXOffggeKUQALgJ/QYkiAwo9hTVlu/r+cpEVTVxygdRpg/Hurjj4eqJecV8U+cdg9eTaPTcyPwEJAGXNaanVsTiBGfI+6WRDaYui547LHWr3/05ps1DYIFUF9v54MPdnH55QeedDcW3J7r5bZdiZgZoe/FUSG5/ajoZUC47mQ/X62xU1wR6o4bkmVwYSu7LuOVN2vsPFPjwAROTwxyew8lkIr2IZ66B4UQEwCflHJ1+HMv4AlgNLAI2Naa43T5qYPXXTeCrKxQd52uCyZN6smQISlRO77tJ/eMSUZGhOR9nYwrxgR5pKeHPmUGA3YG+Xx8HQPTW96vtYzsbzL7Jjc/P8LPeUf6ee92N73TWl5yJF5Z7dO4p9zFcp+NH3w2nqxy8FZt/NwPivjCQI/46mQ8AfRp8vk5YBihdbRGAY+05iDxI9NtZOrUHF5//TjefnsrgwYlc9VVQ1tM4bRjh5uPPtrGwIHJTJ0aeY7TP/85kPPP30Yw6AAM+vUzOOOMnlFswYHj9UJ+vkZiouTQQyNnMfeb8LHPgZkgCAp4dY+DMVmRPS2/P3R8u10ydmzk4wOMHWjywrWt7xKNZ77w2NhlNolOlBrzPDZmpFh7l1LCqnUahZuSGDAwtKioQtEaOqEwRWIE8A1AOPjiNOBQKeVGIcQc4DtaMQe4y4sWwLhxmYwbZ71ablMWLdrJ9Onz8PtNhBCMGpXOwoVnWZY/6aSeLFqUwDPPlJCb6+Laa/uhadF1YHft8vD++wX4/VUMHjwEPUIMe1UVnHGGk7XrNBx2mDbN4Pnn/ZbCMmuDky/32BpW/31lu+CsPgGOzmw+oq6uDs45x8mqVRo2G0yZYvDKK36i3OS4ZYzTIFmY1MnQF2JDMtxhHZ0oJVx2cwJffmvD5x/G2Jcl773gJjVy7I9CAcTdmJaN0Og5wJHATinlRgAp5fawkLVItzA1Cxfu4IYbFvDIIz/g90cOb77iiq/D2TMEUsKaNVUsXLgj4j5VVRKfz05VlYY/ysnCi4qqOfXUd7jjjm+4994fmTHjQwzDOnrw1t86WL1axzQEXq/g/fdtLFhgfZm31Wv7LFdfZ2hsrLMuP2uWnaVLdXw+gdstmDdP58MP4+s22uUR3LHAxe1fuyipjW5QyAmJBpenBuivm/TVDU5PCnBzuvVN8fF8G59+Zcft0QgGdZb/aGPm481Pr2grdSbcW+vk5moXqwOxv1ZePzz0opOZLw3k+9Wxr08848cZ8dXJWAPMCL+/APh87z+EEH1pDIOPSId4WkKIF4AzgF17k/AKIXoAbxBa/XgLcJ6UslKE+u6eBKYBHuAKKeWKtp77gw+KuP3279i924umhTypd9451XJ5Er9/X0GQUlJSYr3A4JdfVnPddUXs3Bnq/vn221rmzBketSS4DzywiM2bQ9fSMOCbb0pYtKiUyZP7NVt+ybJ9n7xME75coDNlSvNCd0rvIF/stjXMW+rrMjihp3X4YHn5vu0KBAQ7dmiEJnB3fvbUC854P4mNlaHvaf52G3Omu+mbEr1xtlk9vfwp00tQQlILNnnnLoE/sO93WlEVPSH1SphemcTyYOin/rnfzstpbsbHKP+jYcC5dybxbb6OlC6+X2fy9995mHqkmivXFuKse/BO4EMhxNOEMoZPbvK/8wktKtwiHfWY8x/g1P223QV8IaXMIxSvf1d4+2lAXvh1DaFVLdvMiy+uZ/fu0BiNacKKFXsoLLQW9MmTs/b57HLpnHJK8wIB8PTTZQ2CBbBihZv8fLdleSklf/3rCqZNm8M553zE6tV7ItY/ENjXuASDJvX11j/wvoPYN0BSg2HDrQ3yBf0C3D7Ux6SMIEdlBPnbmHoGJlmXv/TSIL17N9Zp0CCTs86KH4Pz3I+OBsEC2FSt81R+9AeRnKJlwQI4/aQguQMav79emSYX/zx67vrXPhsrgo3tLTE1nvDE7gl89WaNFet1pAzdpGUVGv9+r9N5BHFDPAViSCkXAgOAk4HBUsoNTf79P+DW1hynQ0QrvJBkxX6bpwN7cxi+RGMi8unAf2WIxUC6ECK7refefyxH08Bms272yy+fyDnnDCAz00Hfvgl8+eXppKdb/6h+enyBrls/KT/99Goee2wl3323g/nzS7j88s+prLQOfLj22rFkZzcuTz92bC8mT7YOp591v4+MvoBdgEMwbDxMnxZ54tVNQ/x8erSbj492c1KvyAJ07LEmf/ubn1NOCXLaaUFefdVHTk78RAM6tJ/WtYWcyO1KVi/Ja095mHZSgGMmVvHE/R5OPCZ6DwF2IX/yI4+lKbPpoO13DaK0tF23pK1LkwghXhBC7BJCrG6yrYcQYp4QoiD8NyO8XQgh/iaEKBRCrBJCjG+yz+Xh8gVCiMtbqq+UslZKuVxKWbvf9g1SytLWtFnI/dfuaCeEEIOAj5p0D1ZJKdOb/L9SSpkhhPgIeDisygghvgDulFIua3q86urqhooXFDRdCnhfli6t5L77NrBrlx+7XXDssT14+OGRUVsEcs2aIHfd5WbnTomuw4tyPmYAACAASURBVKRJNh5/PMmy+/Hmm/P57rtyQqZDAib/+MdhTJrUw/IcK1dW8sor5aSlwa23DiQ5OXII9dqNibwxpzeJCSbXXlZCakr8eELtTV1Q49dLD2FjbRIAQ5I9PHP4BtLssf2Odgfs+ExBjsNPFBfWJgjckDiMH/QUTCHoa3h5wlPAIBmbuWNSwm+fGsp3a1IxDI0+PXw8eM0mRg+27oLvCjRNCpCWlnZQV7ip7TsqbXPEsouqG/M+ND2vEOI4Qlko/tvEJj8CVEgpHxZC3AVkSCnvFEJMIzQJeBowCXhSSjkpPMSzDJhIyJgtByZIKSsPpn0t0RmjB5u7oBGVNVKWiLw8mDAhj/ff38KgQSnMmDHEUlDaQl4ejB3r5e23K8jOtnPhhT0jelqlpRsIZeO3EWqWj5ycXPLyejdbvrZW8txzJhs3gq4H6dfPzoMPRnYN8vJg+ul7P3XaZCWW+P0G3367k507S5kxY0JEz7gtfDEkyH/X1mOYcNmoAOnO2H1HUsLNK1x8ssOON2hwWA948xhPi0vEHAifSMnL9fVUmBoXuvz0tbW4+kO78v5j8PpnXlatr+Sac1MY0j8+JuJ3RtraBSilXBB2JJoyHTg+/P4l4CtC41ANvV/AYiHE3t6v44F5UsoKACHEPELDQK/9f3tnHh5FlfX/z63qJenOxpawJIQAYVFZRkRFQNxQERERZEQdURhfmZ/zOuq8Ksq4izLuMioyOirquLCpjGyig8rgIBAhgGxJIIQsJCGE7Omt6vdHZSGY6gRo6KRzP8/TD1Rzqb5dXXW/95x77jkn1almEkzRyhdCdNF1Pa/mAhTUvJ8NHLs5Kh5oltloRt++7XjoIT/1Qk6RpKQwHniga7PaHjrUnoaeUgubN2tcdnyGxxr+8hedn36qb/vRRzpTpugMGBA44c33CT6qshEhdKY63IQF2F1TVCr4YK0VmwVuv9yN009wXGWll+uvX8nmzQUoCnz8cSFLllyFzRa4UTzSBncPbv66UbkLFvzXhtsLU4d5aO9nze9EWVugsuigjSqfABR+KNSZ/YudpwcGzhKyCqNmWktBVeGWMR7O751HLz/JqyVN4wps7sE4XdfzAGrG5dqZdDfg4DHtsmveM3v/tBJM0VoGTAXm1Pz55THv/1EI8SmGKVpSeyHPFN98U8Bnn2UTGxvGrFl9cTgCN2AK8esBr7ra3DVVWNiwfWkpZGUFTrSyvIIJR51k+Ax35VKXlWXtKrAHSLgKjgrGPeNkT45xDRevt/LVYxWmwvXKK6ls3GjMXzQN1q/PY8GCPdx551mB6dAJUu6Ca193sjXbeFQ+22zjX3dX0ClA0Yb7y5UawapFkFslw8AlzeMMpXEy836dsFcsEJyRp0MI8QlGbqm+QohsIcR0DLEaLYRIw4gmmVPTfAWwD0gH3uYUqyTrus6bb+5g0qRVTJu2lsOH/Wdm+OKLXGbM2MqiRbm88cY+brhhA16v//Dg997L5MYbf2Lq1BTy8vxnk7jiiliOveyqGs7vfmeeQSMpqaFg2myCgQP9K8qSJRqTJ/u45RYfGRn++/5shb1GsAAEGz0qy6sD9yA8v9ReJ1gAW/Zb+Ph78zW5oqKG10/TaPI3O518sMFWJ1gAu/NVXlwTuGi3Kzp7iQ+v/43aWTWu69ZyrCJJyybA0YP5tUFvzfR+Bdwr1hzOiEzruj7F5J8ub6StDtwdqM+eO3c7zz+/hYoKI4IuI6OENWvGmbqbPvggi8OH611HW7eWsHNnGQMHRjfa/r33Mnn88d2UlhrnT0srZ82a4TidjV/aDz9MZsIELykph1FVhVdfPYukJPNB8MABB8ZyujGQud3h7N7tIiGhcets+XKNBx/UKaqpY7lnj87q1TodOjQudD7911nYA7k/2tNIN91ec9GdPr0/q1cfJCfH2DaQmBjBzTf3CWCPToxizQM9wuGcmjfSoVAL3BVKdOr8fWglf91lp7yyilv6WBgvC1JKmkmAw9pPyPslhFgNPFsbZYhR3PHhQHaoMVpiIEZA+fbb7DrBAsjIKCU9vYSzzmo8Wm///oZ7rFwuza+ltWJFfp1ggSFaO3aU+o0G/Pzz/s3tPlZrbe0QI+mvzaYTHm6+3vHZZ/WCBZCeDt99pzNxYuNCcb/TxU8elSzNuPkHWrxcaw/coHn/dS6+32Ehs8A4f/94H7dcYj7on312e0ZeN4mv1lrQNI0xE3QSE4MXpK2fXwBRERBWYx230yDhMGD++54oF3Xy8WWnSqP0TM/Alp6RhDY+7eSejRrv1yVARyFENvA4hlgtrPGEZVGfvWIFRuRgOkbChzsAdF0/IoR4GthU0+6p2qCM00nIi1bYcWFYDofqd99V+/aC/fs1DBeejhAaum4uWuHhDc/vdFpo1y5wi6OPPVbNtm0qGRkqqqoxcqSPiy4yXwMzolrr3cphYRAXZ9qc/laNJTGVvFlpBGL8X4SLiAA6jRPjdL58pIK5y+3YLDoPTHAT4zRvv3illa/WRVLuMjrx6SqNqy6t5NIgZUzIroqpFywAu0Kut/FabBLJmcZVfXKu6kB5v3Rdf5eG9RVPOyEvWnPmXMiBA2vIyCglKsrGTTcl07Wr+ajZo4eDlJQj1IpWp062Bpt7j+e5584mPb2ctLQKnE6ViRO70qdP4CKikpJ0Vq6s4MsvLbjdecyY0dFvctqnnxZs366zcyfY7XDVVTB8uP81sGSrxivRgauhdTyJcTovTWve+b/fpFJeUf8Fj5YqrNtkCZpoXW0NY7Hiw1szoxVojLaEU+uulUiCic8bzK3iwSHkRatnz2heemk477yzl+TkSGbO9F/g8c47k/jqqwO4XMYg2blzGF27hpu279YtnDVrRrBjRynt2tmaJVi7d3v48MNK4uIU7rorAnsToXqxsTp33ukhLa0YVfVf9iQmRrBqlcL27eBwwFlnEbCN1GeCiwb7WLJKp7KmunFUhMaFvwneGs/4KB9ryr18Xa6jIbgg3Me9HaVgSVoGTYqW2nqy1TSXkBetr7/O5p57/suhQ1VYLIJt24r59NPLTAfyl1/+GZfrKEY0p05WVhVpaUdJTjbPmu90WvyuYR3L5s1upk49Qk6OhhDw9dcuvviiQ8AS7AKEhQmGDm1++w0bNObP17HbBU88Iejc2X9fNm9WeOMNOzYbPPpoNfHxgXswbhrnIXWPypr/WPB43EweK7hyRPCyVQgBb3Sr5qgPvLqgoyX0BgFJ68XraUq0Qi+oJ+RF6803d3HokBEy7fXq/PRTIRkZZfTu3fi6RH3QhTE4uVw+KisD98O//HIZOTnGZ+g6bNzoZtMmN8OGmfum3/3QyrKVVjzu3rz0rEK/PoHL0P3TTxp33KGTZ2wpZMsWndWrFWJiGheuLVsUbrvNSW6u4cJLSVFZvbqCDh0CM5gLARf/3sPusQqVlW4uTm4Z7o8YFc7AFhSJ5ITQfE0N4aEnWnIX43FMnpxMu3b1AnLWWe3p3//0ZdNoio8XWnn4yXC+W2dl/U8xTLjFwZHiwFllf/97rWAZ7NkDK1aYD87z59vrBAsgPV1l6VL/D85r/7Zx2ctOLn/Fyaeb/OdN/KlM4b59Dr4rs7LRF82MfeHsqmxdt+nilVYu/52Ty2518vI/ZBliyWnEq/p/hSCtazQ4CWbM6EfnzsaalMUiOP/8TvTqFWnafvz4nnTv3oHwcCfR0ZFMmzYgoCmE7r8/kvh447ILAUOH2hg61Hxge22+HRcCHEA45B1WWbEmcP1xOGoE0CpAFagqRPkJjnM4GgqaqupERZmL3BdbLby8xs7PBy2kZFl47F9hpGab33afFNrI99b/e65bZUmRf6FrSfyyV2HWS2Gk7LDw8y8WXn3fzuKVwe9/WRkUFBjWvSSEqLb4f4UgofmtjuHqqxP4+ONwFi3aT1JSJNOn9/EbmHDffRtITTWSFFdVaTzyyCauvDKejh0DU032vPNsLF7coS4QY8aMCL/rWblHhPEr1TYRkJGlEiiz/+FZgqUZTipUC+g6Xaxuxowx3wf26KMuNm60sGOHiqLoDB/uZeJE876s+sVKSXW9CBWUKazZaWVQfOOfkWDXUDCCHgCs6HS3t44CkwDf/NdCflH99y0tV1j1g4VJY4IXvPHkkxYWLrTg8UDfvhoLF7oJN48tkrQmQs/71yQhL1oA557bkXPP9R91V8u//32wwfHRo242by7g6qsDlxm7Xz8rs2c3nmHjeGLa65QdmxhFgc5xgRvE/7YyjAq7lVpVzLOG8dNuHxed3XjwQ7t2OitXlrNmjYXwcBg92ovFz110dlcfVlXHU5Nfz2nTGNDN/En73y5u1pVaSClXQdMYGaNxa6fWE63n9uo0TMumU1gSvOjN7dsF775roaTEENKCAsFjj1l54YXWc00lfpCiFZqsX5/Np5/uIjExinvuOc+vu89iOX7Q0bAHsbDqjdd6ePnvSl1/IpwaV18WuDt11Z56wQLweQRfbLeYihZAZCTccEPz+nD3KDdbslR+ylQRCK4Z4OEqP+e2KbC0XyV7qxQOZh3gij4JAS8SWHBY8PI8G5om+NP/uOjWJXA+szCHMFytHozLqgg6xQXPJ7dvn6gTLANBXl7r2QIhaQIpWqHHl1+m8cADaykoqERRYP36HJYsmWBaU2vy5ERee20Xum5kPe/SReW88zqf2U4fw1VXeXhriY3KEgECOvXU6dzJ/yC4er/KBzvshFl0nhxeTbyfNad2CRr7M9T6mz8MunQP3CCrKPDu1CoqXKAqENaM5R1FQD+HhqpWB1ywDhcJRk92cqAmrdTyf1tY81kFXTsH5juPvtDLm9008goNoYiJ0hg3KrBWzcEyweObwnD5BLf3czHaJA8lwPnn63TvrpGVZfTH6dQZNSq47tYjR7zMmpVNbm4VN99cxG9/2yGo/WnVtEGDOeRF6733tlFQYFRF1TT4+edDpKcX06dP4/uqHnvsItxuL5s2HSIsTOW550YRGWkeKKHrOi+8kMF33xURFqbw7LP96NfPPNAD4OWXS/j22yqsVsHs2e04+2zz87++JIzKSAVqTplZrvCf7SqXD2l8oFqTqTLtcRcVqYdBFWy8vgPrHvQRY7Ik9/B1LqYeUqnMUUBA3EAft/zG/5OwYIGVRYtsCKHzwAMuLr646X1UziBaq8fy4jw7B4pUqBHPnKMqT75sZ/7zgckI0i9J46X7K/nbp2FoGky43M34SwI3HS52CSaudrK3xBDdTQUqb42q5LL4xn+DLl105s1zM3u2Fa8XLr/cx513Bm967nJpTJiQTmqqsQ0lNTUHtxt+9zspXCdFcApQB5WQF63jgy4URfithKsogmefHdXs87/66j7++tdd+HzGoHHDDSX897+jiI5u3KR4661SXnqphIoKY2Z/222FfPNNZ9q1a9xlmZV1XBZ2DSorzM2PZ9+qpuLrIqg2ZtMH3/Xw+RWx3HFR45bE6AQf86dX8eEeGxYFHj+vmthwc6tjxQqVJ54Io7jYuIb79ql89VU5SUmBdYH5dOMVaHZkKA1jZhXYlRnY0OBrRvq4ZmRF0w1PgjUH1TrBAiioVvhgr43L4s3LtwwfrrFiRcsY3XburGbPnvoJwtGjPr74oliK1snSBt2DIR/yfu+9Q+na1UitZLMpjBgRT1JS84IgmsP77+/H5/NirIPp5OaW8eOP5omO166tqhMsgH37vGzdap71PKpch9pnXAfKBHqpeX8Kt1TUCRYARz1kpfqvRzU20cvzF1Xx8vAq+rbz7zr6/HNbnWAB5OQorAxgSLeuw59Swjh3VQQTd5/DrG2Bidqs5fIRtb9V3Sdy6UWt58nvEKZjUxqqubMVZemIiFAIC2s47DSVxkziB28TrxAk5EVr1KgEli6dwMyZF/Daa1ewYMG1Ac3FV1X1a1eavwwakZENP9tm0+ja1dzg7dlZg91uSMuGvQVE5PpI6mEuLJee1fD8apjgmvPMv2+1D8Z/7eDSryK4+F8R3Ptf/yIRH681qL4cHq7Tq1fg0ix9kGnlsywbBypVcjxhLNhvZXlu4Cyhu29zc+FgH4rQEUJnyNk+Zs5oGVZIc7i0m48r4j2EqToCnXPa+3j6/NbT/+TkMK67LpqImlICycl2nnnmtFdoD13aoGiFvHsQoF+/DsycOey0nDs8XOP4EOdu3fxZHulABEZ9LA04DLQDGs9teM8fj7Lws/9QdbgYUOmcHM85Zw+l8UrXMOeZSPbscrNzlxerVTDx+jCGDjVfUHr6ZzvrDtVvBFu4z8oNPTxc3KVxIZo508WWLSqpqSoWC4we7eHKKwMnWpuPqFRr9d+t3KuQcsTC2K6B+QybFb56u4Iff1bxaTB8iA9b8Pf+NhtFwD+vqGJDvpsKLwyL8+FsRf0HmDs3kWnTOvHzz5nccEMfYmLaxDB0eghRYfKHvFtOkY4dvRw44IG6CqJVdetbjVFWVgUUAomAG5crk9zcJPr2bVy05sxJpaqyuObIR072QTZt6sn55ze+78zpVFi5vCNZWT7CwwWdO/u3UnIr6sPpASq9CvtKFVPRstth6dJKsrIENht07RpY19ToOC9fZlsprcmK0d6mcXlcYJ/MQ4fcLPqwCJ9Pp0dcB7+Vo88EaRkKt04Lp7zsHGY/5eP6a/1/3zLFy089cqgSGgmuDvTVzEvnnAylpTB7tpXiYsG0aR4uvDDw7sfBgx04nVYpWKeKFC3JiVJW5gHK646F8J8qp1u3BMCJkZdJx2ptT//+5hniq6p8xx1rFBf7L/euqoKkpOb9tOMSPazNs3DUbYhEglPjiibKvSsK9OhxetZRrov3sqfMzbIcK25XNbf3EQzvFDhLLj/fy/jxB8nIMNy6P/5YxbJlCXTvHpwcgQdzYMSoCFzlArAwbZqO751KJl7X+G9QiY//c2aSoRoLnVstFTxZ0Z1+ARKu6moYP97Oli3GZOf77xXmz3dzySWtJytJm+L0lcFrsYT8mpbPp/HQQ99x5ZWfMX78EnbvLmr6P50APXo0DG+3WhXCw80FY/PmjhjuQQtgxePpyLZt5gIwZUoiHTvWWwL9+kUxbFgnv32aMCGX+PidJCTs4vXX/URtADckeRme6SXyM52ohRp3hVUT7wzuwv4D/V2su6Kcj/rs4g/Jgd2I8o9/FNcJFkBmpoc33yz28z9OL3NeDqsRLAPNLXjyafN1xY2WsjrBAjiseFloPxyw/vz0k8LW7SqEG6/8Ygvz58u5bYtFrmmFHk8+uZ53392Gx2PMFKdOXc7atVNwOAKzEDBjxtmkphaRn2/sBRsypBODBpmH7x46pHJ8zPX+/ebnHzu2G5qm88knB/B6K3nttRFERZn3/a67Cli7Ng9jvQwee2w/w4b1Z8iQxi2JT7+0sm6phbJyAQjmvWRn3GAviQGskdWSCA//9VpgY++dKi6Xjs+n43D4nxc6HcdnYAGL1fza21FQdDhm2Q+ryfrmyVBVLdDDLNTt6lZ1DuSFZrbwkCBEhckfIW9pbd2aXydYADk5ZezbdzRg5z///Di6dLHjcAiio1WuuSbB7z6wQYMazhNUFYYN8y+g48bF8/HHw5k9uwdduvjPdPrDD6XUChaAprn55z/Nra3V31soLT8mhP2QytofAzuX8Wmw+7BC+hEl6FnG/+d/2nPuufWWzMCBdu69N7B7hGbNKuG88/IZOrSA228/gs/PhrPHZrqI6lgrXGBz6nzwD/MtCud7IxnkdVIb9Z7gszG9Oi5gfVcsomGMjxB0iJUh6S0WaWmFHjHHpYKIjrYTF+cM2PlnzlzP1q317pk33tjG9df3JCGh8awYb7wRzrhxlaSl6djtMHGilYEDA/czdOxobVAfC1SGDjV3NyV0NULYdb0+t2GfnoFbv3B54cZFDrbmq6gCRnb3suD6qoCnZ2ouTqfCV191Z/HiEnw+mDQpisjIwFkSa9e6WLCgkvJyQ1Xy86uZO7ec++5r/H6IcMLO1DJmPWUnJ7eCF56xkpRoLnIqgjmVPfjOUkK58HGpN5poPXD3T2J3jXYxOsVH63+gnklyPavFEqLC5I+QF63p0weyatU+3G7jwevaNYJOnQIXbZWX1zDzQWFhFdnZ5aaiFRur8O9/O9m2zUdMjKB//8C6Xv71rwQGD66ipKQKIQSDBrVnypQI0/az7nGxY7fKjj0qFlVn3GgvF50XuMCHF3+0sy7Lgl4zfV+VbmXJLg+TzjJ/2goLBW+9ZaO0tAuzZkFM44GVJ43DoXDbbaensOe2be46wQLwemHXLv/rchFOeO2vLtLS9pOUmNzkZ1gQXOEN8EWpoW+yxu9vc/PJYituj6B3T43Zj7XB1f7WghSt0GPu3JQ6wQLYt6+EAwdKSEwMTFaMc8+NZd263LrP6N49kr59/Q+ITqdg2LDTc+ljYhT27evLzp1eoqMFCQn+P8dugyVvV1JwWGC368T4KQB5MuSUiTrBAnBrggMl5vXACgoE117rZO9eFejGjz/6WL68PODCdbq4/HI7b75ZQX6+cT9ERgpGjw5sVo/TzawHXPzvDBeVlYK4WD1oVrGkGfhPdhOShLxoeb0NXRtut4+qqsBNTx5++DyKi6vZvLkAm03lyScvoH374A5SiiI455zmB5oIAXFNZI4/WX57joc1+6wUVhrrZt0iNcb3Mbc8XnnFViNYBr/8ovLeezbuu89/mH9L4ZxzbDz9dBTz51fg88G114Zx442B3Ud1ouze7eLRR4vwemH8eCe33970hC0qEqIiQzMYJ6QInFOk1RDyojVhQh9SUws4etRIddO/fwd69za3hHw+jQcfXM/PPxdityvMnj2MIUP8L3QrShRCeBFCQYjg7PepRdfhmRQ7a3MsqAr83+BqrvJTuuJ0MyrRxwujq3h/qw1FwAMXuejdwXyNpLFAjWAHb5wokyc7mDw5uEJVS1GRj1tuOURGhjFR27KlGodDMHlygE1qSXCQ7sHQ4/bbB+B0Wlm0aDddukQwe/bFfqP7nnpqEwsW7MLrNUbKu+5ay7ffTiA6uvGsCXPm7OL99/fjdhvt//CHFL79dhTt25tnWdB1Y93G4dCJMF9uOine3mVl3i92Kr2GT+f+9eGsHFtB9yDOmq/v5+X6fs17uu67z82331pJSzOsrbPO8jFtWuuwsloi69dX1QkWwNGjOsuWVUjRChXaoGiFfMh7UVEVb721iy1bSlm7Np9//nOv3/apqYV1ggWQnV1ORkaJafuUlOI6wQLIyqpgz54y0/bl5TB2rJMRIyK48MJInn++6RRC1dU6mza52bev6cWF/+RZ6wQLIKdCZUN+69lnExens3x5BffeW820abmsWNF61rNaIl27qjidDe8bszI4klaIDHkPPe67bx0pKQV1x3PnbuWGG3rRqVPj+51iYxu+366dna5dzUPkj6+bZbOpxMebu4YeeSSMH3/UgSLAwltvRTFxoodevRp3mR05onHttQXs2VOK1SqYMOEI8+aZp31KjPAhqI/Wi7Fp9ItpXSHLsbE6TzzhIi0tl5iYpqPpTpSqKp3PP69C0+D668PqMo6HIuedF87EiRF89VUF1dU6/frZeOYZWbsqZGiDgZ0hL1q1a1nHHufnV5qK1m239eOLL+pD5Dt1CicuzlyEXK5fC8Lx+QKPZd8+N5BFbcnRI0eiyMyMo1evxtvfc08hO3fmAj58Pli4sIpp0xyme68ePc/F7qMqvxxRUBXBlN5uBnZsXaJ1Oqmo0Lj2ugK2pPgAwVtvKyxfFkt0dOgK19y5sdx/v4eyMo2+fW3YbDIcMGQIUWvKH6H7pNYwaFBHwI0RG1pNbGw4PXua+/Nfey21QYh8VlYZe/eaZ9A4ti0Y0Yr+EtoWFR2iYY3sMoQwj1tNSSni2BAhn8/FunXlpu1tKiy6qpINE8vZPKmMR4a0nlpLZ4J5/zjKlhSN2rQPO1I1Xpob2HyULZEePawMGGCXghVqtEH3YMiL1oEDxw76Onl5R2hYubYhx6fccbk0XC5zy2nMmM4NKq927x7OOeeYhxT36OGr6Y+n5uXDZjM/f3j4r3+iLl2a/tmibGBvIUsXn/7byrhHnIyb5WDtluB2an9V9q/ey6zKDUJPJJIA4GniFYKEvGj98MPBBsfV1V7Wrj1o0hqio1WOFTVd10hMbDy7BUB5+RE07QhQCVTg8+Xj9ZqLUJcugmNzA9rtGr16me+pio2No6EXNxyn09xS1HWdhx+uZOTIEi69tJTPPw+upfXtzyqz/hHGuu0W1m2zcvdr4ezJCt5tN356GY6udnD0AEcPwuOcTPjfI0Hrj6T1sbNU4aofnYz4IYLbUsKpDOZeKV8TrxAk6KIlhMgUQmwXQmwVQmyuea+9EGKNECKt5s+TzrmjKA3dIUJgup4FcPRoJYY70YsxVXH5TbD7ww8H8XhKMQo7HiY7+yi7dpm7m3JyGk5/3G6d3bvNhSUhIQroiVHduCN2ezKxseY/2+uvV/P++y62b9fYssXHI49UkZkZvLt38fdWikrr+5tbpPLVhuCV2lX3XoiudAE1HNRwdDUWZcfFQeuPpHXh1WD6Fgc/FVvYUaay7JCNe7f5T2J9ejvUxCsECbpo1XCpruuDdV0/r+Z4JvCtruvJwLc1xydFv34NI+1UVZCcbB59FxVlx7CEPICX6Gib31yFxyfkjYmx+w3caN++oXssOlqha1fzQfzppxX6948AkrDb4xkzxsIFF5ivS2zc6KPqmCWyvDydDRv8372VlTpr12ps3KihN2Mnb3U1fP+9yoYNKloTMR6JnXUUUX/OMKtOz67BE9Hvv3NQVVI/yFSXhfHdt4FLoCwJbQ65BIXuhs/f/sogDqNVTbxCkJYaPTgeuKTm7wuA74CHTuZEqtrwhrLbLWRnl9GuXePRdy++eCmZmSVkZpbgcFi5/fZziI83dw8+//wlZGQUs39/CeHh6HrR9QAADrdJREFUFm6++SySksw3Fj33XBx797rJyHATFiaYPDmavn3N92p17iz4+muFH3/UKSvLY+LEBISfZHC9e6uoqgdfjS60by8YMMB8HenoUZ3x4zW2bwebDS67TOejj5RfWai1lJbC9ROcbN2qYrXCxRd7+fSTSlSTj/jzjS427lLZkq6iKnDJIC/XDw/eFHDIYB9Oh0ZFzUATHqYz5Dch6keRBJxONp1oi87hY2KtOtiCmLKlDd66ojkz69PaASH2A8UYC0nzdV3/uxDiqK7rMce0KdZ1vYGLsKSkpK7jaWlppuefM2c7S5dm1aUCSkhw8OGHI3E6zfXa69XIzq4kOtpKu3ZNb/71ejVyciqJjLT6zYRR314nJwciIqBDh8BGc3m98Mc/Oti5ExQFbrrJx4wZ5u7HZ5/twOef1weOWK0af/1rPiNHNj5Ne/HFBD5bWJ/WSlU1nn5qH6NHm7tQdR1yjtiwKNC5XeCzW+QV2nhncRcA7pyUS+dO/legX/xbAuv+G4MOXDCklEfuPyCTwkqazX/Ko5l7OJ5KTSHO4ubFrhm0s/ifiCUn1+83jI6OPqW77dixL+Z3/vNIHv2wPjHCqX5uS6ElWFrDdV3PFULEAmuEELtP9ATH3hDHM29eT4T4mp07DxMWZuHpp0cyeHB8k+fs3//E+nAi7f/5z3y+/LKIsDCV557rTrdu/hPsLlqUz+LFhXg8FbzyykASE8196Js2ucjLK6CqSqs5tvLUU/E4nY27MIRoOFXzeBRstq4kJzfeXqfhZ/t8ChZLV5KTO/n9Dn38/mvjpKWl+f1tAXLyBff+2cm+g4ap90tGe/71dgVdY80nY/Pngtdr7Mq0WGxA4DcwnwzN+b6hRGv9vsnA7bobt1YboZsUvM6cpNNCCJEJlGHYal5d188TQrQHPgN6AJnAZF3Xi4Xh2nkNuAYj4ux2Xdd/PtWunyxBFy1d13Nr/iwQQnwOnA/kCyG66LqeJ4ToAhT4PYkfbDaVd94ZE6DenjqfflrAX/6SSXGxIRZpaZV8/fVA00KES5cW8OCDGRQXG3fn5Mk7WL36N8TENP7TvfFGKbm59QtN27Z5+OabKsaPb3zdZupUwfr1Ovn5xnHv3nD11eYTsml3uPjhBwuHDhmilpTkY/z44Ln75n9sqxMsgIwslb9/auOJe/xHTVqCfudLWjNCtJAtJacW1n6pruuHjzmujSWYI4SYWXP8EDAGQ6uTgQuAeTV/BoWgBmIIIZxCiMjavwNXAjuAZcDUmmZTgS+D08PmU17u+VUZlMb44ouiOsEC2LOnipQU81yFS5YU1glWbfsNG8xzIVqPi+lQVQgPNxehUaMUXn9dcNVVOtddp7N4saBjR/P2F16oMf+tSsZe4+G669wsWlhJXFzwXMzhjRipztMQzOVy6XXWa7D56CM3I0ZUcOGF5dx/f3WzgmckIUpgQ97HY8QQUPPn9ce8/4FusAGIqTEmgkKw55txwOc1gQUW4GNd11cJITYBC4UQ0zFyHt0YxD76paLCy5Qp35OeXorNpjJjRl9mzOhr2v54N53TqdKhg3n04PEWmMOh0LGjeftHH41hyxY36eleVBVGjLBz2WXmo7jXq/P++1527NBQVfjwQ5VHH/U/lxk1yseoUZV+25wKmg7puiBLtdFbx+96092/c7HqBwupu41beVB/LzNuDuzetIceOsLy5cb3HTbMzvz5HU0DVU43aWk+nnrKTUGBIVT793tITBT86U9Nr6VKQpCTzz2oA18LIepiCYA4XdfzAGq8XLE1bbsBx25uza55L++kP/0UCKpo6bq+DxjUyPtFwOVnvkcGe/Yc4Z13ttOnTzumTz8HRTEfxGfOTOGHH/Lrjl99dSdjx8aTkNC4O+7ZZ5PYvbuKPXsqcThUxo/vwIAB5iHXs2f34pdfKti1qwK7XTBuXEeGDDGPZuze3crKlZ1ZurSC6GiFSZOcWCzmA+ycOV5WrdLqog3fecfH1VcrDB0aHN+HR4ffVltJ0RS0Tr0Y6RJ8ZPdgphFREbD8HxUsWmFFADeO9QTU0lq5spKPPiqnosIQiWXLKjn33DL+8IfglPbYuFGjoECjdhrtcqn8/HMbDCGTGJy8e/BEYgkae/qCZt4H29JqcSxfvo+pU9fg9QLovPHGNrZuvdW0fV5eQ4ujsLCa7OwKU9GKi7Px9dcDSE0tJybGwlln+d8j1KGDldWrB5OaWkZxcS7XXNOvye/QqZPKXXc1b1BNT9frBAugpAR++UVn6NBm/feA8zePyneagoYAxcIan87HXoVbreauuQgH3DHp9OSsSU111wkWgMsF27cHr77XoEFgtbrxeGqvh5fExOBt1pYEmZOcr5xgLEE2kHDMf48Hgpb7rKVsLm4x/OlP39cIFoDgwIFyvv02y7R9x442jp10WCzQu7e5JQSGS/Cii6KbFKxaHA6VYcNiSE4O/OA0apSC85hudOkCI0cG9rbYtq2M6dN3cdddu8nK8r/j8YAuDMGqwYMgUw9epO7o0eF06lR/PaKjBWPG+DflsrN9zJhRxfTpVWzdGtggldxcDZ/vWAHXKSyURTLbLCeREeMkYgmWAbcJgwuBklo3YjCQltZxaFpDq1fX4cgRc8fxzp3ZGFvPrYCOy1XJvn0lflNFtSTuuMNCbq7ON99oKArcc4+FXr0CJ1o7d5Zz8807yc421plSUspYsWIQsbG2RttPsfhY6VUoqJlPdRUaE9XgBUAMGWLn8cdjePfdcjQNJkxwMG6c+WTj8GGNCROqSEsz7qMNG3x88kkYAwcG5lErLdV+lYWkug3WVJLUcHJzohONJViBEe6ejhHyfscp9vqUkKJ1HNdc04OPPtpLrRs3LExhzJgepu1zckqBigbvpaUd4YILOp++TgaYWbOszJp1es795pu5dYIFkJ5exSef5POnPyU02v5CVecVu5e3PSqVlZU8Em2jvxrc6Lhbb43k1lv9W8+1LF7sqRMsgJwcnXnzPMybF5hH7Yor7PTrV8nu3YZfKC5OYfr01jFBkpwGTmLCcqKxBLoRnnr3SfTutCBF6zg6dLBhJMw1AhFiY8N/Va7kWJKToykqqhctVRUMHNjxNPey9RAR0TCgQwiIifEf5DHWojHWopGWnUlyh9a1+TQqSkFRaGANHV/u/lSIiVH44osYnniigqoqjd//3sGIEY1brZI2QIgmxfWHXNM6jtTUAuoT5nooKKggPb3YtP3cuVfQrZsxC7dYBNde24sBA/xnh2hLzJyZWBcdqSgwbFgUU6a0Hiv0RJk0ycKwYUpdLsYBAxRmzQpsOHrnzipvvRXFggUxjBwpBatN0wazvEtL6zhiYxtmaG/fPrxOlBqjT5/2fPfdFNasyaRzZyeXXtrdb0LbtkZMjIWVKwezalURNptgzJgOWK2hO1ey2QRffOFg1Sov1dVw9dUWIiPl/SA5TYRooUd/SNE6jhdeuJSsrFIyM0sIC1P5f//vXDp39h/l16mTg5tvPusM9TD4ZGaW8/rr6YSHq/z5z32JifE/24+IUJk0KdZvm1DCahWMG9d2wtCr3fDqQjuHjgimXu3mN31aRuaQNkEb3KInRes4oqPtrFx5IyUlLpxOK1ZrS0gw1nLIyCjnhhvWc+CAsT/tu+8KWbFiJJGRbWeQltTj9cENs5z8uEMFBKs3Wpj/QBUXD2qDo2kwCFEXoD9C109zCgghiIkJk4LVCC+9tKdOsAC2by9hyZLsIPZIEky2pqmk7DEECyCvSOWNpTKl1BlDFoGUSPyjqr9en/GXJkoS2iiCX6XYkku6Z5A2aNBKSysI+Hw6mZkuiopa3yrqzJn96N07ou54yJB2TJzY+J4rSegzONnHsHO8qIqxLSQhVuPBKXK38xlDb+IVgkhL6wxTUuLlxhv3kp5ejd2uMGVKRx57rOmilC2Fbt0crFgxkvffz8ThUJk2LYnwcOlGbasoCix8qpIPVlk5dETht5e56dUtREdLSYtAitYZ5pFHDrJxY+1mZB/vvVfATTd1oE+f1pPVIDY2jAcfbDpxr6RtYFFh2tjW5zWQtE6ke/AMc/hww4f76FEfOTky4alEIpE0BylaZ5hRo6JwOOove1KSncGDm5ftXSKRSBrS9sIHpXvwDPOHP8RRXu5j7dpS7HbBU08l0K6d/BkkEsnJ0PbcsnK0PMMIIXjwwW48+GC3YHdFIpG0epraXRx6QVJStCQSiaTV0pSlJUVLIpFIJC2GpkQr7Iz04kwiRUsikUhaLU0FWzSveGlrQoqWRCKRtFraXsZcKVoSiUTSapHRgxKJRCJpNUhLSyKRSCStBmlpSSQSiaTVIC0tiUQikbQaQjNVkz+kaEkkEkmrRboHJRKJRNJqkO5BiUQikbQapKUlkUgkklaDFK2QxKtDNtAOiBbB7o1EIpEECukeDDkKdLhZg4OAA/i9gLtl6UuJRBIStL3owZAfvh/UYBuCYgQ5CN7WDSGTSCSS1o+niVfoEfKWVmkjx4VArJ//o6Z8j23lJ2gdO+O64yGwh5++DkokEslJI92DLQohxNXAaxiVzN7RdX3OiZ7jXAEbdR03xmJWN6Cnn/bqD1/hePE+lOJCdMCyYyMVry4DS4u+VBKJpE0SmtaUP1rsSCyEUIE3gNEYcRSbhBDLdF3feSLneUSAD9is64QDcwSE+wnGsH/5HkpxodEHQNmbirJ/J1rywJP7IhKJRHLaaHuWltD1lrnAI4QYBjyh6/pVNccPA+i6/hxASUlJy+y4RCKR+CE6+tRimE927DvVz20ptORAjG4YQX+1ZNe8J5FIJJI2SksWrcZmBdK6kkgkkjZMi13TwrCsEo45jgdyaw9CxdSVSCSSE6Gtj30t2dLaBCQLIZKEEDbgJmBZkPskkUgkkiDSYi0tXde9Qog/AqsxQt7f1XX9lyB3SyKRSCRBpMVGD0okEolEcjwt2T0okUgkEkkDpGhJJBKJpNUgRUsikUgkrQYpWhKJRCJpNUjRkkgkEkmrQYqWRCKRSFoNUrQkEolE0mr4/3U/26phLD2CAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# scatter plot of Years versus Hits colored by Salary\n", "hitters.plot(kind='scatter', x='Years', y='Hits', c='Salary', colormap='jet', xlim=(0, 25), ylim=(0, 250))" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Index(['AtBat', 'Hits', 'HmRun', 'Runs', 'RBI', 'Walks', 'Years', 'League',\n", " 'Division', 'PutOuts', 'Assists', 'Errors', 'NewLeague'],\n", " dtype='object')" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# define features: exclude career statistics (which start with \"C\") and the response (Salary)\n", "feature_cols = hitters.columns[hitters.columns.str.startswith('C') == False].drop('Salary')\n", "feature_cols" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# define X and y\n", "X = hitters[feature_cols]\n", "y = hitters.Salary" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predicting salary with a decision tree\n", "\n", "Find the best **max_depth** for a decision tree using cross-validation:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# list of values to try for max_depth\n", "max_depth_range = range(1, 21)\n", "\n", "# list to store the average RMSE for each value of max_depth\n", "RMSE_scores = []\n", "\n", "# use 10-fold cross-validation with each value of max_depth\n", "from sklearn.model_selection import cross_val_score\n", "for depth in max_depth_range:\n", " treereg = DecisionTreeRegressor(max_depth=depth, random_state=1)\n", " MSE_scores = cross_val_score(treereg, X, y, cv=10, scoring='neg_mean_squared_error')\n", " RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'RMSE (lower is better)')" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEGCAYAAADWjcoaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XecpFWV//HPqeqcJzM5NkHCIAgMioqA5CwIBkBBXd1VYZU17JpXl8UsGNYV+AkYAJEsCCgCq8yAgMOQBnryDJNTh+lcdX5/PE9PV1Wn6p6u6qru7/v1qlfXE6rqdE1Nnb73ufdcc3dERETyQWSkAxAREUmXkpaIiOQNJS0REckbSloiIpI3lLRERCRvFIx0AENVX1+vYY8iIqNYdXW1pe5TS0tERPKGkpaIiOQNJa0sqqurG+kQBk0xZ0e+xZxv8YJizpZMx6ykJSIieUNJS0RE8oaSloiI5A0lLRERyRtKWiJD0BZzlm5vZ3dbfKRDERlT8nZyschI2doS46QHtrGuKUZhBE6YVsy5c8s4fVYJ1UX6O1Akk5S0RAbpO0sbWdcUA6AjDg9vaOPhDW0UReDE6SWcN7eUU2eWUKUEJjLssvq/ysyiZvYPM3sgZf/1ZtaUsF1sZreb2Qoze9rM5mQzTpG+bG2JcWvdnl6PtcfhofWtfOzJXdTetokP/HkHd65qprFDXYgiwyXbLa0rgVeBqq4dZvYWoCblvCuAXe6+wMwuBq4FLspalCJ9+PkrTbTGBj6vLQZ/WNfKH9a1UhKFk2cELbCTZ5RQXqgWmMhQZe1/j5nNAM4AbkjYFwW+A3wu5fRzgJvD+3cCJ5pZj8KJItnU0B7nF8uTW1mfOqSCzx5WwbzKaJ+Pa43BfWtb+fDju1jw28186C87uXdNC82daoGJDFY2W1o/JEhOlQn7Pgnc5+6bUnLSdGA9gLt3mlk9MAHYnqVYRXq4afkeGtq7FxcYXxzhC4dXUl4Y4UtHVLFsZwf3rG7h7jUtrGnsvTnWEnPuWdPCPWtaKC8wTp1ZwrlzS5mTRutNRMDcM7/Ch5mdCZzu7v9sZscDVwMfA+4Ajg8TU5O7V4Tnvwyc4u4bwu2VwNHuvqPrOROXJsnH+lySX9ricPbfS9nZ0f3H1T/Naucjszp7nOsOy/cYf9pWwJ+2R9nYNnCHRnHEOaYmxtvHxzhufIyJRcMavkjeqK2t3Xu/t6VJspW0rgEuATqBEoJrWm3hrTU8bRawKryO9TDwNXdfbGYFwGZgkicEm4/radXV1SX9g+QDxRy4afkePrN4997tigLjxffux7ji/hOSu/P89g7uXh20rjbsSa9JdeTEQk6dWcJps0o5eFwBudY7rs9Fdoz1mHtLWlnpHnT3LwJfBOhqabn7mYnnhC2tBeHmfcBlwGLgAuAxz0Z2FelFZ9y57qXGpH0fOqB8wIQFYGYcOamIIycV8Y2jqnhuWwd3r2nmntUtbGzu+5rWc9s7eG57B9/6RyMzyqOcNrOE02aV8Lb9iimO5lYCE8mmXJ2ndSNwq5mtAHYCF49wPDKG3ZNyjaowAv98cMWgnydixlGTizhqchHfPKqaZ7a2c/fqFv6wrrXfFtiGPTF+sXwPv1i+h4oC44TpxZw2q5STZxQzoaTvASAio1HWk5a7Pw483sv+ioT7rcCF2YtKpHfuzg+WJbey3regjGnl+5YsImYsmlLMoinF/PcxzoNLV/KyTeaP61t5fntHn49r6nTuW9vKfWtbiRgcM7mIU2eWcOrMEvavzr1uRJHhlqstLZGc8OiGNl7e1T3YwoArD6ns+wFDYGbsX+GcUVvF5w6vYnNzjIfXt/LQ+lYe39ja57ywuMPiLe0s3tLOV59tYF5llCMnFTGvqoD5VQUsqCpgXlUBNWl0Y4rkCyUtkX784MXkVtY5c0qZX53Z/zb7lUW57IByLjugnObOOE9sbOOP61v54/pWtrT0fR1sVWOMVY0tPfZPKI4wv6qA+dVBMptfFWV+mNAqNNFZ8oySlkgfFm9pY/GW9qR9Vx06+GtZ+6KsIMJps0o5bVYpcXeWbu/gobAV9tLOvrsRE+1oi7NjWzvPbGvvcWy/0gjzqgpYECa0rlbatLIo1UWm7kbJOUpaIn34Ycq1rBOmFXP4CE6giphxxKQijphUxH8cUcW6pk4eDltg/7epjfYhFNjY3BJnc0s7T23pmdCKIjC5NMrk0kjSzymlEWK7o2yvamNKaZRJpREqCpTgJDuUtER68dLODh7e0Ja0718PG95rWftqVkUBHz2ogo8eVEFjR5yl2ztY1dDJioZOVjZ0siq8DSWZQVAAeMOeWDiyMbVVVwzLuwvUlBVYkNRKupPb9PIo52ahO1XGFn2aRHrxo5RrWW+ZVMhx++VumYrKwghvn1rM26cWJ+2PxZ0Ne2I9ktmK+k7WNsWIDdPsx+ZOZ01jrEf5qmv+0cA3jqrmE28qV0tMhoWSlkiKNY2d/H518oCGfz20Mi+/dKMRY3ZlAbMrC3jX9ORjHXFnXWMsKZmtbOhkTWMnW1vi7Onc94zW6fDvz9Tz963tXHdcDZUa+CH7SElLJMX1LzURT/i+PrCmgNNmlYxcQBlSGLFgRGEf3XdNHXG2tcTZ0hJja0ucrS0xtrTE2dYSY9X2BvZEStnaGuxvG6A61d1rWnhlVwe3nDCeA2oKM/Db5A9356kt7Ty7rZ2BCv3v2FHAhObG/k/KMYdiZLLwlJKWSIItzTF+lbLI45WHVhLJw1bWvqoojFBRGGFuVc+vibq67dTWzgaCL+H6dmdba5DUtjbHWNcU4/vLGmno6M7+r9V3cuL927j+uBrOm1uWtd8jV7g7j2xo49tLG3iunwnkyYpgbUNG4xpu1xyY2da0kpZIgv95pSmp1TCjPMoF80pHLqA8YGbUFBs1xRFqq7v3nzW7lEv+soNXEiZnN3U6H358F89sbecbR1VTGBn9fwy4Ow+ua+XbLzTywo50k5X0RR3MIqH69jg39rLI41j4Ys2E+dUF/OnMSVw0v2fS/9krezjroe1sah69C4nF3bl3TQtvv28bH3hspxLWMFFLSyR00/I9Sd1ZE4ojXLL/2OvGGk5lBRH+5+3jOHpyEV94up6OhGs4S7a28877tnLT8eM5br/ivp8kz8TiQbL6zguNvLq753prXU6ZWcLB4/r/Ct65cxfjx48b7hAzak60Z1WW4aSkJQK0dDo/fbkpad/H31ROWYE6I/aVmXHFgRUcPqGIy/6yM6mi/daWOOf8cTtfO7KKTx5SkZcjNLt0xp27Vrfw3Rcaeb2+72R11uwSrl5YycIJA0+hqKvbSm1in2seqKvbmtHnV9ISAX5dt4dtrd3NgIoC46MHZbdk02h35KQinjh7Eh95Yhd/2dg9cTvm8OVnG3hmWzs/OW4cVUX59YdCZ9y5Y2Uz31vWyMqG3rs7DTh3TilXL6zk4PFje/TkvlLSkjEvWOQxuZV1+YHlqo6eARNKotz57glcs7SR776QPJT7/rWtvLJrK7eeMIE3jcv9L/aOuHPbiiBZpU6q7hIxeM/cUj67sJIDx/hQ/+GipCVj3l2rW1jX1P2lUxSBTwxhkUdJTzRifOmIKo6aVMTHntxJfXv3dcSVDTFOemAbP3prDRfOz83rie0x5zdhslrf1Heyeu+8IFnVVitZDSclLRnT3L1HYdz3LyhjaplWBM60U2aW8MTZk7nksZ28mFCxvrnT+eiTu3hmWzvfOqqaoujgr3O5Ozva4ryxJ8am5hgb98RZ/kYBNU0NwcRxhzhO3Om+JWw74HuPOXG6z3tyU1ufK01HDS5eUMZnD6tkXi/z22Tf6V2VMe3hDa28kjDCK2Lw6UNzqzDuaDansoBHzpjEvy3Zza/qmpOO/eLVPSzd3s4v3zWB6QkrRcfizpaWOBubYwlJKZa0vam5tyodRbAuM9UlCgw+UFvGvx5WyZxKfa1mkt5dGdN+uCz5Wta5c0r1F3KWlRYYPz4uGBb/b0t2JyWbv2/r4B33buVt+xWxsTnGpj1xNrcMX6HffVUUgUv2L+fKQyuYVaHPTTboXZYx66nNbSzZmryO1JVZXuRRul26fzmHjS/k0r/sTLrGuKMtzn1rW0cwsp6Ko3DZ/uVceWhlUitQMk9JS/LC/WtbuH9tCzPLo7x3ftmwFF39Qcq1rJOmF6c1d0Yy5/CJRTxx9mT+6cmdPJKyntlgVRYa08qiTCsPboUt9UyZOJ6IBaWAImZELBiOHjEwC/eF293HLGm7vMA4bmoxk0uVrEaCkpbkvLtXN/Phx3ft3f7esiaOmFjI+xaU8Z65pYwvGfyXx4s7O3j0jdxe5HGsGlcc4baTJvDdFxq55h+N9NYTOL44wrTyKNPLgp9Tw+Q0Pfw5tSzaY75XUOS3Kju/hGSMkpbktO2tMa5eXN9j//PbO3h+ez3//kw9p8wo4X0Lyjh5ZknadQJTRwwePamIt05RKytXRMz43OFVnDqzhMVb2hm3N0lF2a8sSmlB/lbOkH2jpCU57QtP17Ojre9Fhzri8MC6Vh5Y18qE4ggXzCvluGJjgXufJYFWN3Ry95qURR4Py+8SQqPVYROKOExdtpJAU/4lZz20roU7VyUnlwX9jOzb0Rbn56/u4ZKlpbztnq1c/2Ijm3upIn7dS41JizweVFPAKTNH3yKPIqNRvy0tMysAzgbOABYCNcBu4AXgIeAed++7MqTIEO1ui/OZxbuT9h06vpDHzprEpuYYd6xs4bcr9vRZ6+2V3Z18+dkGvvpcAydOK+biBWWcPquU+vY4v06ZD3TVYWNzkUeRfNRn0jKzfwL+A3gVeAJ4AGgEKoGDgI8C3zez/3L3/8lCrDKGfOXZejY1d3cLRg1+fFwNhRFjVkUBVy+s5LOHVfD3be38dkUzd61uSSoH1CXu8OgbbTz6RhtVRbuZU1FAe0Jv46yKKO+Zq0UeRfJFfy2t/YGj3X1zL8fuBv7LzKYCn81IZDJmPbGxlVteT24NXXloRY/h6GbG0ZOLOXpyMdccXcND61u4bUUzf9rQSoyeLaeGdmfZzuSF+D51SAUFWuRRJG/0mbTc/bMAZhYBjgf+6u7tKedsAq7OZIAytuzpiPPpvyV3C+5fXcDnFvY/VLmkwDhvbhnnzS1j8ct1POf78dsVzby8q+/e64klET5YWz4scYtIdgw4EMPd48C9qQlLJBP+8/kG1iZUQzDg+rfVUDKIIc4Ti+CTh1Tyt3On8OTZk/jEm8qZWNLzo/4vB1do6LRInkl39OCTZrZoX1/MzKJm9g8zeyDc/rWZvWZmL5nZTWZWGO43M7vOzFaY2TIzO2JfX1ty39Nb2vj5K3uS9n3soHKOmTL0pdgPm1DENcfU8OpF+/HbE8dzzpwS5lZG+fABZXzyEJVsEsk36c7TWgs8ZGb3Auuhe5K6u39lEK93JcHAjq6+nl8DHwzv/wb4CPAz4DSgNrwdE+47ZhCvI3mmtdP51N92J1U/mFUR5ctHDk8Fg8KIcdqsUk6bpUEXIvks3ZZWKXAPQbKaAcxMuKXFzGYQDJ2/oWufuz/oIeCZ8LkBzgFuCQ8tAWrCQR8ySn3nhQZer0++/nTd22qoKNRUQhHpZkG+yMILmd0JXEMwZP5qdz8z4Vgh8DRwpbv/X9h9+N/u/tfw+J+Bz7v7s12Pqa+v3xt4XV1dVn4HyYzXmozLlpYkjfg7Z0onX6rVZVSRsaa2tnbv/erq6h4XndMu42RmBwEXAFPc/ZNmdgBQ7O7L0njsmcBWd3/OzI7v5ZSfAk+6+/91PaSXc/rMrom/ZC6rq6vLm1i7ZDrmjrhz+f3biNE9FH1qWYQfnTiTmuKhtbL0PmdevsULijlbMh1zWt8KZnYh8CQwHbg03F0JfD/N13kbcLaZrQFuA04ws1+Fz/1VYBLwmYTzN5Dc9TgD2Jjma0keue7FpqSl1gG+d2zNkBOWiIxu6X4zfAN4t7t/HOgaj/wCQWmnAbn7F919hrvPAS4GHnP3D5rZR4BTgPeFQ+u73AdcGo4iXATUh3PCZBR5bXcH1y5tSNr3nrmlnK7BEiLSh3S7BycTJCno7qZz+umyS9P/EIxMXBxW2L7L3b8BPAicDqwAmoEP7+PrSI6JxZ1P/XV3UkmlCcURrl1UPXJBiUjOSzdpPQdcAtySsO9ighF/g+LujwOPh/d7ff1wNOG/DPa5JX/876t7eGZb8kCLaxdVM3EICzqKyNiRbtL6NPCImV0BlJvZwwS1CU/OWGQyaq1p7OQ/n0/uFjx1ZokK14rIgNJKWu6+3MwOBM4kqPa+HnjA3ZsyGZyMPu7OlX/bTXNnd89yVaHx/WNrtAijiAworaRlZte5+6eBO1L2/9Ddr8pIZDIq3VrXzBOb2pL2ffPoaqaVq1tQRAaW7ujBD/Wx/5JhikNyiLvzzecaWPi7zXzixWL+95Um3tjT+2KLg7FxT4wvPVOftO+dU4u5pLZsn59bRMaGgVYuvrzrvIT7XeYB2zMSlYyoh9a38t1ljQCsJcqzT9fzuafrecukQs6aXcqZs0qZX532vHQgSISfWbybho7ubsGyAuNHb1O3oIikb6Bvnq6WVBHJrSoHtgCXZSIoGVmpy9F3eXZbB89u6+CrzzbwpnEFnDm7lLNml3LIuIIBE8/vV7fwx/WtSfu+fEQVcyoHl/xEZGzr9xvD3d8FYGbfdPcvZSckGUm72+I8uqF1wPNe2dXJK7sa+fbSRuZURjlrdilnzS7hLZOKiKQksO2tMT6/JLlb8OhJRXzsIC3AKCKDk+41rTN622lmz/a2X/LX/Wtbkib8TiyKc9Skwn4fs6YxxvUvNXHyH7bzpts3c/Xi3TyxsZWOeNAV+Pkl9exo637Soghcf1wNUS1zLyKDlG7fzPzUHRb0B80b3nBkpN25qiVp+9wpMb570mQ27onxh3Ut3L+2lb9tbiPWRy2UzS1xbli+hxuW72FcsbFocjEPpXQLfv7wKg6o6T8Rioj0ZqCBGF0VMIoT7neZA7yciaBkZGxujvFkynD0UyYFa1xNK4/y0YMq+OhBFexojfHQ+lbuX9vKX95oTWqZJdrV5j0S1qHjC/n0oVoxWESGZqCW1so+7jvwN+B3wx6RjJi7VrckFZNcOKGQOWU9B2VMKInywdpyPlhbTkN7nD9tCBLYIxta2dPZdznKqMGPj6uhUN2CIjJEAw3E+DqAmS1x94ezE5KMlDtXJSeoC+aVAvW9nxyqKopw/rwyzp9XRkun8/jGIIE9uK6F3e3JCeyqQytYOKFouMMWkTEk3TJOD5vZuwmK5E5297PM7C1Albs/ltEIJStWNXTy/Pbuda0MOH9uGc2DWMWstMA4bVYpp80qpSNew1Ob27hvbStrGjtZNLmIzy6sHP7ARWRMSbeM06eAK4EbCFYvBmgBrgPempnQJJtSW1lv3a+I6eVR6ob4fIUR453TSnjntJJ9D05EJJTukPergJPc/b+Brsvuy4EDMhKVZJW787uUUYMXzlNpJRHJPekmrUqCyu7QvfBjIdDe++mST5bt7KCuvnPvdmEEzp6tFpKI5J50k9aTwBdS9n0a+MvwhiMjIXVu1gnTSxivxRhFJAelO7n4U8D9ZvZRoNLMXgMagLMyFplkRdydu3p0DWoxRhHJTemOHtxkZkcBRwGzCboKn3H3PqaVSr5YvKWdN5q7lx0pKzBOm6muQRHJTel2D3ad21V7J0owKlryXOqowTNmlVBeOJiPhYhI9qQ75P0w4B6gGHgDmAG0mtl57v5CBuOTDGqPOfesSe4avECjBkUkh6X7J/VNwE+AGe5+NDAd+HG4X/LUYxtb2dXWXbViXLHxrmnFIxiRiEj/0k1a+wM/dHcHCH/+CKjNVGCSeT0qus8ppSiqXl8RyV3pJq0HgbNT9p0F/GF4w5Fs2dMR58F1yRXY1TUoIrmuz2taZnYr3ROJo8BtZvYcwcjBmcCRwL0Zj1Ay4qH1rTQnVGSfXhbl2CkqZisiua2/gRgrUrZfSrj/CqCq73kstWzTe+aVEjF1DYpIbuszaXUtSyKjz87WGH/ekNw1+B5NKBaRPKAJOWPQvWtaSVyrcf/qAg4bX9j3A0REcoSS1hh05+qeiz2augZFJA8oaY0xG5o6eWpzcnF+jRoUkXyR1aRlZlEz+4eZPRBuzzWzp82szsxuN7OicH9xuL0iPD4nm3GOZnevbiGhZ5AjJhYyryrduskiIiMrraRlZu8zs4PC+weY2ZNm9piZHTjI17sSeDVh+1rgB+5eC+wCrgj3XwHscvcFwA/C82QYpI4aVCtLRPJJui2tbwI7w/vfBZ4hWGPrp+m+kJnNAM4Abgi3DTgBuDM85Wbg3PD+OeE24fETTRdd9tnruztYtrNj77YB58/VqEERyR8WVmbq/ySzBnevMrMSYBOwH9ABbHf38Wm9kNmdwDUEqyBfDXwIWBK2pjCzmcBD7n6Imb0EnOruG8JjK4Fj3H171/PV19fvDbyuri6dEMa8n68t5Ib13aMEj6qO8dND20YwIhGRZLW13dUBq6urezRW0r2Ysc3MFgCHAn939zYzKyPN5UnM7Exgq7s/Z2bHd+3u5VRP41gPib9kLqurqxuxWN2dx17YAnSvnXXpIROorS3v93EjGfNQKebMy7d4QTFnS6ZjTjdp/SfwHME33kXhvhOBdJcleRtwtpmdDpQAVcAPgRozK3D3ToLlTjaG528gKBW1wcwKgGq6uydlCP6xvYNVjd0JqygCZ81W16CI5Je0rmm5+y+BqQRLkzwa7n4auDjNx3/R3We4+5zwMY+5+weAvwAXhKddRnctw/vCbcLjj3k6/ZjSp9+lLPb47hkl1BRrxoOI5Jc+v7USBz6YWQRoJVj4MRJubwe27uPrfx74jJmtACYAN4b7bwQmhPs/A3xhH19nTIvFnbtWJ48avFCjBkUkD/XXPVhP0I0H0EnPa0oW7osO5gXd/XHg8fD+KuDoXs5pBS4czPNK3/66uZ0tLfG92xUFxikzS0YwIhGRoekvaR2ccH9upgORzLkzpWvwjNkllBZoBoGI5J/+qryvT7i/NjvhyHBrizn3rlXXoIiMDroSP8o9uqGVhvbunt2JJRHeOa14BCMSERk6Ja1R7s6Usk3nzSmlMKKuQRHJT0pao1hDe5w/rk+tNai5WSKSvwZMWmFl9pVmpj6lPPPgulZau+cTM7MiytGTi0YuIBGRfTRg0nL3GEElDI2RzjOpowYvmKvFHkUkv6VbxumHwB1m9l8EJZb2XtkP51pJjtnWEuMvG5OL4WoZEhHJd+kmrR+HP9+dsn/Qk4slO+5Z00IsYTr4m2oKOHh8Yd8PEBHJA2klLXfXgI08kzpq8IL5amWJSP4bVDIys5lmtihTwcjwWNvYydNb25P2abFHERkN0kpaZjbLzP4GLAf+FO67wMxuyGRwMjSpxXGPnlTEnMp0e4JFRHJXui2tnwN/IFh1uGu99kfpeY1LckDqMiSamyUio0W6f34fDZzh7nEzcwB3rzez6syFJkPx8s4OXtnVuXc7anCuugZFZJRIt6W1BViQuMPM3gSsG/aIZJ/8fnVyK+udU4uZXKoBniIyOqSbtL4LPGBmHwYKzOx9wO3AtRmLTAbN3XuOGlTXoIiMIukOeb/JzHYCHwPWA5cCX3b3ezIZnAzO37e1s66pu25TcRTOnK2kJSKjR9pDysIEpSSVw25fmdzKOnVmCVVFmmInIqNHWknLzP4BPA48ATzh7rsyGZQMXnvMuSvlepYWexSR0SbdP8OvBhqAq4A3zGyZmV1vZhdkLjQZjEc2tLKrrbtu07hi4+QZqnEsIqNLute0/gz8GcDMJgCfAT4J/DOqPZgTbl+Z3Mo6f24ZRVFVdBeR0SXd7sFTgXeGt5nAYuCLBN2FMsJ2t8V5eH1r0r6L5msAhoiMPukOxHgQWAlcA9zi7p0DnC9ZdPfqFtrj3dvzKqMcNUmLPYrI6JPuNa13ADcBFwLrzOwRM/sPM3t75kKTdKV2Db53fpkWexSRUSmtpOXuf3X3a9z9NOBw4O/A5whGFMoIWtPYyZKUiu4XaRkSERml0r2mdR5wPME1rf2B5wgWhtQ1rRGW2so6ZnIRc6tU0V1ERqd0v92uJEhQnwEWu3vLAOdLFrg7t69ITlpqZYnIaJbukPfjMxyHDMGz2zpY1dhdtqkoAueporuIjGLpLgJZaGZfN7NVZtYa/vy6mWmI2ghK7Ro8eUYJ44pVtklERq90v+G+DZwEfBxYGP48gTSrvJtZiZk9Y2YvmNnLZvb1cP+JZva8mS01s7+a2YJwf7GZ3W5mK8zsaTObM8jfa9Rrj3mPZUguWqCuQREZ3dK9pnUhsNDdd4Tbr5nZ88ALwL+m8fg24AR3bzKzQuCvZvYQ8DPgHHd/1cz+GfgS8CHgCmCXuy8ws4sJkuNFaf9WY8CjKWWbaopUtklERr90W1p9TfpJazKQB5rCzcLw5uGtKtxfDWwM758D3BzevxM40TTxKElvZZuKVbZJREa5dFtavwPuD7v11gGzCVpFd6T7QmYWJRgqvwD4ibs/bWYfAR40sxaCgryLwtOnE6zbhbt3mlk9MAHYnu7rjWa72+L8UWWbRGQMMncf+KRgwMWXgPcD0whaRL8FvunubYN6QbMa4G7gU8A3gGvDBPZvwAHu/hEzexk4xd03hI9ZCRyd0D1JfX393sDr6uoGE0KS5hhsbzdmlQ78PuSKuzZHuWZF8d7tGSVx7jqyFbVFRSTf1dbW7r1fXV3d41st3SHv7cBXwts+cffdZvY4cBrBdbKnw0O3A38M728gKMy7wcwKCLoOd/b1nIm/ZDo2Nce47sVGFm9p58WdHSycUMhjZ00e5G8yeHV1dYOOtTeP120DuqtgfODAavbff+Y+P29vhivmbFLMmZdv8YJizpZMx9xn0jKzE9J5And/bKBzzGwS0BEmrFKCkYjXAtVmtr+7vw68G3g1fMh9wGUE1eQvAB7zdJqEaYoa/OyVPXu3X9jRQVNHnIrC3B8uvqaxk8VbVLZJRMam/lpaN6bxeAfmpXHeVODm8LpWBLjD3R8ws48CvzezOLALuDzhtW81sxUELayL03iNtE0ujTK/Ksrme7QEAAAWFUlEQVTKhmBibszhuW3tvHNa7o++uyNlAMbRk4qYp7JNIjJG9Plt5+5zh+tF3H0Z8OZe9t9NcH0rdX8rwTD7jDl2SjErG7oTwOItuZ+03L3HqMGLFmgAhoiMHbnfH5Yhi6YkF/NI7XLLRc9t79jbOgQojMB5c5S0RGTs6DNpmdnfzezCvko1mVmRmb3XzJ7u7XiuO3ZycdL2s9va6Yjn9gjC1OK4J88oYXxJdISiERHJvv4uhlxGMCT9Z2H1i9eARqCSYHmSI4DHCCpY5J15VVEml0bY2hIs+bun03lxRwdH5OiKv0HZpuTi+hqAISJjTZ8tLXd/xd0vAA4BbgVagIlAM3ALcLC7X+Tur/b1HLnMzFg0OaWLcGvudhH+6Y1WdrbF927XFBmnzMzta3AiIsNtwGFn7r6ZIGmNOsdOKea+td2VJZZsaeNfDq4YwYj6ljoA47y5pSrbJCJjzpgdiAFwbC+DMYZxOtiw6b1sk7oGRWTsGdNJ65DxhVQUdLdWtrfGWdnQOYIR9e7eNS20dQ8aZE5llGMm5+a1NxGRTBrTSasgYhyV8uX/VA4Ofb8tpWvwvfPLUNF7ERmLxnTSgp7ztZbkWNJa21vZpnnqGhSRsanfpGVm16VsX5Gy/ftMBJVNx05Jnq+1ZMugitZnXGrZpqMmFTK/WmWbRGRsGqil9aGU7e+kbL97+EIZGUdOLCThsharGmNsaY71/YAsCso2aW6WiEiXgZJW6oWTUXchpbwwwsIJhUn7luTIfK3nt3ewImFgSGEEzp+rsk0iMnYNlLRSx3/n3njwYZDaRfjU5tzoIkwdgPFulW0SkTFuoIsjBWb2LrpbWKnbo+IbdNGUIn78cvd2LrS0OuLOXavUNSgikmigpLUVuClhe0fK9tZhj2gEpI4gfHFnB40dcSpHcFHIP7/Ryo6Esk3VRcapKtskImNcv0nL3edkKY4RNbEkyv7VBbxeH1w/ijv8fWs7J0wfuSRx+4rkVtZ5c1S2SURk0E0JMzvAzM4zs9mZCGik5NL6WvXtcR5cn9I1uEBdgyIiA83T+p6ZfTBh+1LgZeB/geVmdlqG48uaXJqvlVq2aXZFtEdFehGRsWiglta5wJMJ2/8FfNrdJwEfB76aqcCyLbV47rPbOmiPjcxgydSK7irbJCISGChpTXL3dQBmdggwAbgxPPYrgsUgR4XZFVH2K+1+O1pizrKdHVmPY11TJ3/bnNw1ebFGDYqIAAMnrXozmxLefzvwrLt39ZsVMoomG5tZjy7CxSMwX+t3KRUw3qKyTSIiew2UtO4AbjOzTwNfAH6TcOwYYGWmAhsJPQZjZHm+VlC2KblrUHOzRES6DZS0vgA8TlBj8H+BnyccOzzcN2r0VvE9m4tCLt3RsXfYPUCBqWyTiEiigeZpdQBf7+PYjzIS0Qg6ZFwhlYVGY0eQqHa2xXm9vpMDagoHeOTwuG1Fz7JNE1S2SURkr36TVjjEvV/ufsvwhTOyohHj6MlF/PmN7mtZS7a0ZyVpdcSd369W2SYRkf4MdIX/l8AKYDO9D7pwYNQkLQjmayUmrae2tHHZAeUZf93H3mhje2t32aYqlW0SEelhoKR1HXAB0EiQnO5JGD04Ko3USsapAzDOnVNKScGoGZwpIjIs+h2I4e5XAbOBnwLnA2vM7Bdmdlw2ghsJR04sIrFO7tqmGBv3ZHZRyPr2OA+uU9egiMhABqw96O4xd/+Du18EHADsAh4PlygZdUoLjDdPSG1tZbZxedeqFloT8uLMimiPCh0iIpJmwVwzqzazfwL+CJwH/CewNJOBjaRsztdyd25Y3pS076L5ZURUtklEpIeBCuaeaWa/A14F3gz8m7vXuvvX3X1Xui9iZiVm9oyZvWBmL5vZ18P9ZmbfMrPXzezVcBJz1/7rzGyFmS0zsyP24XcctNRWTiavaz2ztZ2Xd3XPzYoYXLa/ugZFRHoz0ECM+4DXgF8DLcApZnZK4gnu/pU0XqcNOMHdm8ysEPirmT0EHATMBA5097iZTQ7PPw2oDW/HAD8Lf2bFMSkV1V/a2UF9e5zqouFfFPLG5XuStk+ZUcLMCpVtEhHpzUDfjrcQDGufuC8v4kFZia4+sMLw5sAngPe7ezw8r2sl5HOAW8LHLTGzGjOb6u6b9iWOdI0viXJgTQHLdwctICdYFPKkGcM7BH17a4x71iQPwPjIQZkfXi8ikq8GqojxoeF6ITOLAs8BC4CfuPvTZjYfuMjMzgO2ESx7UgdMB9YnPHxDuC8rSQuCLsKupAWweEvbsCetX73eTHv31CzmVkZ517Tivh8gIjLGDbkfyswOA77s7hemc767x4DDzawGuDtc6qQYaHX3t5jZ+cBNBNXk+5rI3Ku6urpBxz+QOR4Nwws8tqaei6u29v2ANHXFGnP435dKSLyseNaEFlauWLHPrzHcMvH+Zppizrx8ixcUc7bsS8y1tbX9Hh+ojFMZ8EWC4rh1wNcIugq/R1BE9+bBBuTuu83sceBUghbU78NDdwP/L7y/geBaV5cZwMa+nnOgX3Iozt2vk6++vmXv9qt7osyat4Di6NBH9dXV1e2N9ZH1rWxs27H3WHEUrjp2FuNzrNZgYsz5QjFnXr7FC4o5WzId80AjC34CnAW8ApxEkGCeAF4G5rj7v6TzImY2KWxhYWal4XMtB+4BTghPeyfwenj/PuDScBThIqA+W9ezusyqiDK9rDuBtMZg6fbhG0V4Y8ow9/PnluVcwhIRyTUDdQ+eAhzu7lvN7HpgHfBOd/+/Qb7OVODm8LpWBLjD3R8ws78CvzazfyUYqPGR8PwHgdMJ6h42Ax8e5OvtMzNj0ZSipCK2S7a2c8yUfb/mtLaxk0c2JE9YvuJADcAQERnIQEmromtEn7tvMLOmISQs3H0ZwTyv1P27gTN62e9AWq24TDo2JWk9taWdKw/d9+f95Wt7ki7QLZxQyJETs7P8iYhIPhsoaRWE5Zr2XshJ3Xb3xzIU24hblNKqenpLG3H3fapW0RZzbnk9uTjuFQeWY6qAISIyoIGS1laCEX1ddqRsOzBvuIPKFQfVFFBVZDS0B+2i3e3Oa7s7OWjc0FtF961pYUdb8hIk79HqxCIiaRlontacLMWRk6IRY9HkoqTrT4u3tO9T0kqtgPH+BWWUFw5/pQ0RkdFI35YDSO0i3JeK73V7jCUpxXc1AENEJH1KWgNILZ67LxXff78puWH7jqnF1FZrAIaISLqUtAbw5glFJNbJXd8UY0NTZ98P6ENDe5wHtyYnLbWyREQGR0lrACUFxhETU5YqGUJr6/aVzbTEu0cITi2LcPqs4a1lKCIy2ilppaFHF+Eg19dyd25KGYBx6f7lFEY0zF1EZDCUtNKQOhhj8SAHYzy1pZ1XEyrGRw0u219dgyIig6WklYZjJhcllZ1/dVcnuxPmWg0ktZV1xqwSppWrzqCIyGApaaWhpjjCQeO6B1E48HSa17W2NMe4b23yQo8agCEiMjRKWmk6dojztW6ta6YjoVFWW13AO6ZqoUcRkaFQ0krTUAZjxOLOL19L7hq8/ADVGRQRGSolrTQtmpyctJ7f3k5rZ5+LKQPw8IZWNuyJ7d0ujjjvW1CWkfhERMYCJa00zagoYEbC4In2OPxjR/+trdQ6g6dOilFTrLdcRGSo9A06CG8dRBfh6oZO/vxG8nWvC6Z2ZCQuEZGxQklrEAZTPPemlGtZR04s5MCK/rsTRUSkf0pag5A6GGPJ1nbi3jMRtXQ6v6pLTloa5i4isu+UtAbhgJoCaoq6R/41tDuv7OpZPPeeNS3sautOZuOKjfPmagCGiMi+UtIahIgZx6TRRXjj8qak7Q8sKKe0QMPcRUT2lZLWIA00GGPp9nae3ZY84OJydQ2KiAwLJa1BSp2vtXhLG55wXSt1AMaJ04uZV5W8jpaIiAyNktYgHT6xiOKEWrcbm+OsDycQ726L87uVyXUGLz9ArSwRkeGipDVIxVHjyIm9dxHetrKZllh3q2tGeZRTZmqhRxGR4aKkNQQ9hr6HXYSpFTA+dEA5BVroUURk2ChpDUFqxffFW9p5clM7dfXdw98LDC6p1TB3EZHhpKQ1BEelLAq5fHcn31vWmHTO2XNKmVKmhR5FRIaTktYQVBdFOHh8YdK+Jzclz9dSBQwRkeGnpDVEqde1Eh1UU9BjPpeIiOw7Ja0hOnZy30np8gO10KOISCZkJWmZWYmZPWNmL5jZy2b29ZTj15tZU8J2sZndbmYrzOxpM5uTjTgHI7Xie5fyAuOi+RqAISKSCdlqabUBJ7j7QuBw4FQzWwRgZm8BalLOvwLY5e4LgB8A12YpzrRNK48yu6LnQIv3zi+lqkgNWBGRTMjKt6sHulpSheHNzSwKfAf4XMpDzgFuDu/fCZxoOdjftqiX61aXH1gxApGIiIwN5r2sB5WRFwoS1HPAAuAn7v55M7sSiLj7D8ysyd0rwnNfAk519w3h9krgGHff3vV89fX1ewOvq6vLyu+Q6p7NUb61orubcGFVjBsO63thSBER6V9tbe3e+9XV1T0aK1mr5OruMeBwM6sB7jazdwAXAsf3cnpvrao+s2viL5lNn5gT5+ZNW9mwJ0ZhBP77uCnU9nGtC4LkOlKxDpVizo58iznf4gXFnC2Zjjnr5cfdfbeZPQ68i6DVtSLs+SszsxXhdawNwExgg5kVANXAzmzHOpCKwghPnTuZx95o45DxBSyoLhz4QSIiMmTZGj04KWxhYWalwEnAc+6+n7vPcfc5QHOYsADuAy4L718APObZ6sccpKqiCOfOLVXCEhHJgmy1tKYCN4fXtSLAHe7+QD/n3wjcamYrCFpYF2chRhERyXFZSVruvgx48wDnVCTcbyW43iUiIrKXJhSJiEjeUNISEZG8oaQlIiJ5I2uTi4db4uRiEREZfXqbXKyWloiI5A0lLRERyRt52z0oIiJjj1paIiKSN5S0hpGZzTSzv5jZq+Fil1f2cs7xZlZvZkvD21dGItaUmNaY2YthPM/2ctzM7LpwUc5lZnbESMSZEM8BCe/fUjNrMLOrUs7JiffZzG4ys63hygVd+8ab2aNmVhf+HNfHYy8Lz6kzs8t6OydL8X7HzJaH//Z3d5Vk6+Wx/X6Oshzz18zsjYR//9P7eOypZvZa+Nn+wgjHfHtCvGvMbGkfjx2p97nX77esf57dXbdhuhGUqzoivF8JvA68KeWc44EHRjrWlJjWABP7OX468BBB9f1FwNMjHXNCbFFgMzA7F99n4B3AEcBLCfu+DXwhvP8F4NpeHjceWBX+HBfeHzdC8Z4MFIT3r+0t3nQ+R1mO+WvA1Wl8dlYC84Ai4IXU/6/ZjDnl+PeAr+TY+9zr91u2P89qaQ0jd9/k7s+H9xuBV4HpIxvVsDgHuMUDS4AaM5s60kGFTgRWuvvakQ6kN+7+JD1XKEhc5PRm4NxeHnoK8Ki773T3XcCjwKkZCzTUW7zu/oi7d4abS4AZmY5jMPp4j9NxNLDC3Ve5eztwG8G/Tcb1F3O44O17gd9mI5Z09fP9ltXPs5JWhpjZHIJ6i0/3cvhYM3vBzB4ys4OzGljvHHjEzJ4zs4/1cnw6sD5hewO5k4wvpu//3Ln2PneZ4u6bIPgiACb3ck6uvueXE7S6ezPQ5yjbPhl2ad7UR5dVrr7Hbwe2uHtfq9uO+Puc8v2W1c+zklYGmFkF8HvgKndvSDn8PEFX1kLgeuCebMfXi7e5+xHAacC/hAt0JhrUopzZYmZFwNnA73o5nIvv82Dk3HtuZv8BdAK/7uOUgT5H2fQzYD5wOLCJoLstVc69x6H30X8ra0Tf5wG+3/p8WC/7hvReK2kNMzMrJPgH/bW735V63N0b3L0pvP8gUGhmE7McZmpMG8OfW4G7CbpNEnUtytllBrAxO9H16zTgeXffknogF9/nBFu6ulfDn1t7OSen3vPwwvmZwAc8vEiRKo3PUda4+xZ3j7l7HPhFH7Hk1HsMYMGit+cDt/d1zki+z318v2X186ykNYzCvugbgVfd/ft9nLNfeB5mdjTBv8GO7EXZI55yM6vsuk9w0f2llNPuAy4NRxEuAuq7ugNGWJ9/keba+5wicZHTy4B7eznnYeBkMxsXdm2dHO7LOjM7Ffg8cLa7N/dxTjqfo6xJueZ6Xh+x/B2oNbO5Yav9YoJ/m5F0ErDc3Tf0dnAk3+d+vt+y+3nO9giU0XwDjiNo8i4Dloa304GPAx8Pz/kk8DLBSKUlwFtHOOZ5YSwvhHH9R7g/MWYDfkIw0upF4C058F6XESSh6oR9Ofc+EyTVTUAHwV+bVwATgD8DdeHP8eG5bwFuSHjs5cCK8PbhEYx3BcH1iK7P9P+E504DHuzvczSCMd8aflaXEXypTk2NOdw+nWAU3MqRjjnc/8uuz3DCubnyPvf1/ZbVz7MqYoiISN5Q96CIiOQNJS0REckbSloiIpI3lLRERCRvKGmJiEjeUNISyTMWVLDvdR7PMDz3HDPzcJKrSM5R0hIZw8JlLk4a6ThE0qWkJSIieUNJSyQNYYvk38Kq4XvM7EYzmxJWkG80sz91VRI3s9+Z2WYLFqF8sqvCvJkVhQv3fSrcjprZ32yABSrNrNTMfmlmu8zsFeColOPTzOz3ZrbNzFab2acTjn3NzO60YIHBRjN73swWhsduBWYB95tZk5l9LuFpP2Bm68xse1goVyQnKGmJpO89wLuB/YGzCJbo+HdgIsH/pa5k8RBQS7BEw/OEVdE9WLPpg8A3zOwgggXzosC3BnjdrxJULJ9PsC7R3lVfzSwC3E9Q1mc6wfpiV5nZKQmPP4egEv544DfAPWZW6O6XAOuAs9y9wt2/nfCY44ADwuf7ShivyIhT0hJJ3/UeVA9/A/g/ghWc/+HubQTVtt8M4O43uXtjuP9rwEIzqw6PvQR8Mzz/auASd48N8LrvBb7lwQJ664HrEo4dBUxy92+4e7u7ryKoan5xwjnPufud7t4BfB8oIViBuj9fd/cWd++qc7dwgPNFskJJSyR9iUugtPSyXRF2+f23ma00swaCpdEhaI11uRmYQ1AEta+F/hJNI3kBvcRVmmcD08xsd9eNoPU3JeGcvY/1YKmODeFz9mdzwv1moCKNOEUyTklLZHi9n6A77iSgmiA5QfIieD8FHgBOMbPj0njOTSSvRTQr4f56YLW71yTcKt399IRz9j427E5MXMtIFbMlryhpiQyvSqCNYNmUMuC/Eg+a2SXAkcCHCK6B3RyuBNufO4AvhmsRzQA+lXDsGaDBzD4fDtiImtkhZpY4WONIMzs/nHt1VRjfkvDYFoLlLkTygpKWyPC6haD77g3gFbqTA2Y2C/ghcKm7N7n7b4BngR8M8JxfD59zNfAIwVpRAITXw84iWFZ+NbAduIGgldflXuAiYBdwCXB+eH0L4BrgS2HX4tVD+YVFsknraYmMYmb2NWCBu39wpGMRGQ5qaYmISN5Q0hLJAeEk5aZebv8+0rGJ5BJ1D4qISN5QS0tERPKGkpaIiOQNJS0REckbSloiIpI3lLRERCRvKGmJiEje+P9TKq9f6l8hkQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot max_depth (x-axis) versus RMSE (y-axis)\n", "plt.plot(max_depth_range, RMSE_scores)\n", "plt.xlabel('max_depth')\n", "plt.ylabel('RMSE (lower is better)')" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(340.034168704752, 2)" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# show the best RMSE and the corresponding max_depth\n", "sorted(zip(RMSE_scores, max_depth_range))[0]" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "DecisionTreeRegressor(criterion='mse', max_depth=2, max_features=None,\n", " max_leaf_nodes=None, min_impurity_decrease=0.0,\n", " min_impurity_split=None, min_samples_leaf=1,\n", " min_samples_split=2, min_weight_fraction_leaf=0.0,\n", " presort=False, random_state=1, splitter='best')" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# max_depth=2 was best, so fit a tree using that parameter\n", "treereg = DecisionTreeRegressor(max_depth=2, random_state=1)\n", "treereg.fit(X, y)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
featureimportance
0AtBat0.000000
2HmRun0.000000
3Runs0.000000
4RBI0.000000
5Walks0.000000
7League0.000000
8Division0.000000
9PutOuts0.000000
10Assists0.000000
11Errors0.000000
12NewLeague0.000000
6Years0.488391
1Hits0.511609
\n", "
" ], "text/plain": [ " feature importance\n", "0 AtBat 0.000000\n", "2 HmRun 0.000000\n", "3 Runs 0.000000\n", "4 RBI 0.000000\n", "5 Walks 0.000000\n", "7 League 0.000000\n", "8 Division 0.000000\n", "9 PutOuts 0.000000\n", "10 Assists 0.000000\n", "11 Errors 0.000000\n", "12 NewLeague 0.000000\n", "6 Years 0.488391\n", "1 Hits 0.511609" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute feature importances\n", "pd.DataFrame({'feature':feature_cols, 'importance':treereg.feature_importances_}).sort_values('importance')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Predicting salary with a Random Forest" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "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=1, min_samples_split=2,\n", " min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,\n", " oob_score=False, random_state=None, verbose=0, warm_start=False)" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.ensemble import RandomForestRegressor\n", "rfreg = RandomForestRegressor()\n", "rfreg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tuning n_estimators\n", "\n", "One important tuning parameter is **n_estimators**, which is the number of trees that should be grown. It should be a large enough value that the error seems to have \"stabilized\"." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# list of values to try for n_estimators\n", "estimator_range = range(10, 310, 10)\n", "\n", "# list to store the average RMSE for each value of n_estimators\n", "RMSE_scores = []\n", "\n", "# use 5-fold cross-validation with each value of n_estimators (WARNING: SLOW!)\n", "for estimator in estimator_range:\n", " rfreg = RandomForestRegressor(n_estimators=estimator, random_state=1, n_jobs=-1)\n", " MSE_scores = cross_val_score(rfreg, X, y, cv=5, scoring='neg_mean_squared_error')\n", " RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'RMSE (lower is better)')" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEGCAYAAADWjcoaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4HNXZ9/HvvdKqWpbcO9jGphcDARxCCBgIHQKhpQF5SCGNJIQ3pBIglTypBAJ5kpAEEkLvndAhhOJgTEduYONeJFlW373fP2YkrWSVlayt+n2ua6/dOTuzc47W1q05c59zzN0RERHJBZFMV0BERCRZCloiIpIzFLRERCRnKGiJiEjOUNASEZGcUZjpCgxWbW2t0h5FRPJcZWWlJW7rSktERHKGgpaIiOSMvAta1dXVma5C2qit+UltzU9q69DIu6AlIiL5S0FLRERyhoKWiIjkDAUtERHJGQpaIiKSMxS0REQkZ+TsjBjbIu7Oyi0xVjXEWdkQo7Ylzpk7lme6WiIi0o9hGbRaYrD7zWs6tiMGH59VRmHE+jhKREQybVh2D5YUGqOKOwNU3GFdUzyDNRIRkWSkJWiZWYmZPW9mL5vZa2Z2SVj+DzN7y8xeNbNrzCwalh9iZrVmtiB8XDTUdZpUVtBle9WW2FCfQkREhli6rrSagXnuvhcwBzjKzOYC/wB2BvYASoHPJBzzlLvPCR+XDnWFJncPWg0KWiIi2S4t97Tc3YH6cDMaPtzd72vfx8yeB6amoz4AExW0RERyjgXxJA0nMisA5gOzgCvd/cKE96LAc8BX3f0pMzsEuBVYAawELnD31xI/L3E9rcFMznj1O1H+vDzasf3pqa18cXrrgD9HRESG1uzZszted19PK23Zg+4eA+aYWRVwu5nt7u6vhm//HnjS3Z8Kt/8LbO/u9WZ2DHAHMHvrTw0kNrC6urrLdm92j22B5TUd280llcyePWqArcqsZNuaD9TW/KS25qdUtjXt2YPuXgM8DhwFYGY/AMYB5yfsU+fu9eHr+4ComY0dynpMLOvadHUPiohkv3RlD44Lr7Aws1LgcOBNM/sMcCTwMXePJ+w/0cwsfL1/WM8NQ1mnrbIHFbRERLJeuroHJwF/C+9rRYCb3P0eM2sD3gGeDWPUbWGm4CnAF8L3G4EzfIhvvk0uV9ASEck16coeXAjs3UN5j+d39yuAK1JZp7ElEQoN2sJQWNviNLTFKSscluOtRURywrD9DR0x2yrtfXWDZsUQEclmwzZoAUzqloyxUl2EIiJZbVgHra0GGGsqJxGRrDasg1b3DMLVutISEclqCloJ1D0oIpLdFLQSrFIihohIVlPQSqDuQRGR7DbMg5ayB0VEcsnwDlrlW19ppWvWexERGbhhHbQqohEqop2z3rfEYWOz7muJiGSrYR20YOuxWiuVjCEikrWGfdBSMoaISO5Q0NK6WiIiOUNBq3v3oKZyEhHJWgpa6h4UEckZClpawVhEJGcoaCl7UEQkZyhodUvEUPegiEj2GvZBa0JZAZawva4pTktMs2KIiGSjYR+0ohFjXGnXH8OaRl1tiYhko2EftEDJGCIiuUJBi62nctK6WiIi2UlBC5isWTFERHJCYV9vmlkhcAJwLLAXUAXUAC8D9wN3uHtbqiuZalt1D2pWDBGRrNTrlZaZfR5YAnweWAz8GDg3fF4MfBZYYmbn9ncSMysxs+fN7GUze83MLgnL/2Fmb5nZq2Z2jZlFw3Izs8vNbJGZLTSzfba5pX3YuntQQUtEJBv1daW1I7C/u6/u4b3bgZ+Y2STgG0mcpxmY5+71YWB62szuB/4BfDLc53rgM8BVwNHA7PBxQFh2QBLnGZTJCloiIjmh1ystd/+Gu682s4iZzTOzoh72WeXuF/R3Eg/Uh5vR8OHufl/4ngPPA1PDfU4Erg3f+g9QFQbIlNg6e1CJGCIi2ajfRAx3jwN3unvLtpzIzArMbAGwFnjY3Z9LeC8KfAp4ICyaAixPOHxFWJYSPS1PEsRRERHJJn0mYiR40szmhlc9g+LuMWCOmVUBt5vZ7u7+avj274En3f2pcNt6+ojePru6urrP7f7rBkVWSosHp93S5ix4cxEjkv3pZNBA25rL1Nb8pLbmp21p6+zZs3t9L9lfy+8A95vZnQRXQB0BxN0vGkhl3L3GzB4HjgJeNbMfAOMIEj7arQCmJWxPBVb29pmJDayuru6zwb2ZvHA1yzZ33ssqmzSd2VXRAX9OOg22rblIbc1Pamt+SmVbkx2nVQrcQRCsphIElPZHv8xsXHiFhZmVAocDb5rZZ4AjgY+F3ZDt7gLODLMI5wK17r4qyboOitbVEhHJfkldabn7p7fxPJOAv5lZAUGgvMnd7zGzNoKruGfNDOA2d78UuA84BlgENADbev7+K6glSkREsl7Sd23MbBfgFGCCu3/ZzHYCit19YX/Hhvvs3UN5j+cPswm/lGzdhsJEzYohIpL1kuoeNLNTgScJMvjODIsrgF+lqF5pp0lzRUSyX7L3tC4FjnD3c4H23+YvE0ztlBe2GmCsqZxERLJOskFrPEGQgs7MQaePNPRco6mcRESyX7JBaz7B4N9EZxDMYpEXul9prVYihohI1kk2EeM84CEzOwcoN7MHCeYm/HDKapZm3a+0VjfGiMWdgkhP45xFRCQTkk15f9PMdgaOA+4hGGB8T8J8gjmvtNCoKjJqWoIez5jDuqb4VsFMREQyJ9nswcvdvcHdb3L3/3X3G8IZ23+T6gqm09ZdhLqvJSKSTZK9p3V2L+Xd73PltEnl3QcYK2iJiGST/lYu/p/2/RJet5sJrE9JrTJkq/taSsYQEckq/d3Tar+SKqLrVZUDa4CzUlGpTNl6KiddaYmIZJM+g5a7HwpgZj9y9++lp0qZoxWMRUSyW7L3tI7tqdDMXhzCumRc9/kHlYghIpJdkg1aO3QvsGBa9plDW53M0lROIiLZrb9EjGvDl8UJr9tNB15LRaUyRfe0RESyW3+JGIt7ee3AM8DNQ16jDBpbEqHAgoHFADUtTmObU1qoWTFERLJBf4kYlwCY2X/c/cH0VClzCiLGxNIC3ku4wlrdEGPGyKSXHRMRkRRK6p6Wuz9oZkeY2Z/N7G4AM3ufmc1LbfXST4tBiohkr2SncfoKcBVQDRwcFjcCP0pRvTJGi0GKiGSvZLMHvwYc7u4/A9qniXgT2CkltcogTeUkIpK9kg1aFQQzu0Pnwo9RoGXIa5Rh3a+0NJWTiEj2SDZoPQl8q1vZecBjQ1udzFP3oIhI9ko2Le4rwN1m9lmgwszeAuqA41NWswyZrEQMEZGslewikKvMbD9gP2B7gq7C59097/rOus/0rqAlIpI9ku0ebN83Gr4uAPJyxG1P3YPu3sveIiKSTkldaZnZnsAdQDHwHjAVaDKzk9z95RTWL+1GFkUYUWjUtwWBqjkGm5rjjC4p6OdIERFJtWSvtK4BrgSmuvv+wBTgirC8X2ZWYmbPm9nLZvaambXPtPFlM1tkZm5mYxP2P8TMas1sQfi4aGDN2jZbdxHmXS+oiEhOSjZo7Qj8xsN+svD5t8DsJI9vBua5+17AHOAoM5tLMH/h4cA7PRzzlLvPCR+XJnmeITFJyRgiIlkp2aB1H3BCt7LjgXuTOdgD9eFmNHy4u7/k7suSrEPadB9grKAlIpIder2nZWbX0TmQuAC4wczmE2QOTgP2Be5M9kRmVgDMB2YBV7r7c/0c8n4zexlYCVzg7mlbBmVSqYKWiEg2st4y48zsB8l8QPtM8Emf0KwKuB34iru/GpYtA97n7uvD7ZFA3N3rzewY4Lfu3qUrsra2tqPi1dXVA6lCv25YWcgvlxR1bJ88sZVvz2od0nOIiEjPZs/u/HVfWVnZJVO91yutgQajZLl7jZk9DhwFvNrLPnUJr+8zs9+b2dj2oNZdYgOrq6u7bA/GntFGWLKxY7shWsHs2WO26TNTYSjamivU1vyktuanVLZ1IOO0Bs3MxoVXWJhZKUHyxZt97D/RzCx8vX9Yzw3pqCsoEUNEJFulJWgBk4DHzGwh8ALwsLvfY2bnmdkKgnFfC83sT+H+pwCvhve0LgfO8N76MVNRWc2KISKSldKyJK+7LwT27qH8coKg1L38CoJxYBkxoVsixrrGOK1xJxrJy0lARERyRrqutHJKUYExrqTzR+PAGl1tiYhkXLIrF3/MzHYJX+9kZk+a2aNmtnNqq5c5W62r1ahZMUREMi3ZK60fAe3pdL8AnidYY+v3qahUNuiejLFyi660REQyLdl7WuPcfY2ZlQAHESRKtAI9pqDnAyVjiIhkn2SD1jozmwXsAbzg7s1mVkaeLk8CW0/ltFpBS0Qk45INWj8kmIIpBpwelh0G5NWyJIm6X2mtVNASEcm4ZFcu/quZ3RS+bgiLnwPOSFXFMm3r7kElYoiIZFpfE+Za+4BeM4sATQmvIY/vZ0EP2YO60hIRybi+rrRqgZHh6zY6Z3xvZ2FZXi7pq6mcRESyT19Ba7eE1zNSXZFsM7o4QlEEWsJewc2tzubWOBVRjccWEcmUvmZ5X57wuqeVhfOamTGxrIB36zuvsFY3xKioVNASEckU/Qbuw+TuGYRblIwhIpJJClp92HoqJ93XEhHJJAWtPkzsnoyhqZxERDKq36BlZgVmttjMitNRoWyyVfegMghFRDKq36Dl7jGCmTBKUl+d7KKpnEREskuy0zj9BrjJzH4CrCBhzJa7L0lFxbLBRE2aKyKSVZINWu2rCB/RrTxvBxfD1t2DmspJRCSzkp17cFgmbHRPxFjdECPuTsTydnJ7EZGsNqBgZGbTzGxuqiqTbcoKI1QWdQaoNof1TbraEhHJlKSClpltZ2bPAG8C/wrLTjGzP6Wyctlg6y5C3dcSEcmUZK+0/gDcC1QQrFgM8DBb3+PKO1rBWEQkeySbiLE/cKy7x83MAdy91swqU1e17LBVBqGmchIRyZhkr7TWALMSC8xsV+DdIa9Rltmqe1BTOYmIZEyyQesXwD1m9mmg0Mw+BtwIXJaymmUJTeUkIpI9kgpa7n4N8E3gVGA5cCbwfXf/RzLHm1mJmT1vZi+b2WtmdklY/mUzW2RmbmZjE/Y3M7s8fG+hme0z4JYNEd3TEhHJHsne08Ld7wDuGOR5moF57l5vZlHgaTO7H3gGuAd4vNv+RwOzw8cBwFXhc9pNLlfQEhHJFkkFLTN7iSCwPAE84e6bBnISd3egPtyMhg9395fCz+9+yInAteFx/zGzKjOb5O6rBnLeobD1VE5KxBARyRQL4kI/O5kdBhwMfIggk3ARnQHslqROZFYAzCdI6LjS3S9MeG8Z8D53Xx9u3wP8zN2fDrcfAS509xfbj6mtre2oeHV1dTJVGJSYwweeKSVGZ2B9+sAGioflHCEiIqk3e/bsjteVlZVdrmqSncbpEeARADMbA5wPfBn4IknOPRjOFj/HzKqA281sd3d/tZfde5onqdfomtjA6urqLttDYcJLq1iZcIU1YvIMplck3bOaMqloa7ZSW/OT2pqfUtnWZGfEOMrMfmpm/wYWENxr+jaw50BP6O41BF2NR/Wx2wpgWsL2VGDlQM81VDTbu4hIdkj2cuE+YDHwU4J7TW0DOYmZjQNa3b3GzEqBw+k7Xf4u4MtmdgNBAkZtJu5ntQsyCFs7trWulohIZiR7Z+Zg4BqClPd3zewhM/uumX0wyeMnAY+Z2ULgBeBhd7/HzM4zsxUEV1ILE+YyvA9YQnDv7I8E3ZAZs/UKxkrGEBHJhGTvaT0NPA381MzGA18lGLd1KUnc03L3hcDePZRfDlzeQ7kDX0qmbumw9VROutISEcmEZFPeTwIOIcge3JEgC/AKggzCvDep26wYL29oyVBNRESGt2TvaX2VIECdDzzr7o2pq1L22XVUtMv2U6tbeGJlEx+aXJKhGomIDE/JTuN0iLv/wN0fHW4BC2CvMVEOnFDUpex7L9QRT2KMm4iIDJ1kU96jZnaJmS0xs6bw+RIzK+r/6NxnZvxov66rsLyysZUbFw+7+C0iklHJZg/+nCBN/Vxgr/B5HsNglvd2+4wr4qMzSruU/Wh+HY1tutoSEUmXZIPWqcAJ7v6Qu7/l7g8BJwGnpa5q2ef7+46kKOEn9l5DjKter+/9ABERGVLJBq2eplXqqzwvTa8o5HO7jOhS9uuFm1nfpBR4EZF0SDZo3QzcbWZHmtkuZnYUwTIlN6Wuatnpgr0qqCrqjNWbW53LXtqcwRqJiAwfyQatbwL/Aq6kc4zWY8CFfR2Uj6qKI/y/OSO7lP3lrS1U17b2coSIiAyVZFPeW9z9Inef5e5l4fP33b051RXMRp/ZuZzpFZ2zZLQ5XPxiXQZrJCIyPPQ6uNjM5iXzAe7+6NBVJzcUFxg/2Hckn368cy3Me99t4t+rmzlwYnEGayYikt/6mhHjz0kc78DMIapLTvnI9FKuHFfPi+s6uwW//0It/zpuXE8rMYuIyBDoNWi5+4x0ViTXmBk/3K+So+9b31E2f30rty9t5OSZZRmsmYhI/tKi8dvg/ROKOW67rvMPXjK/juaYBhyLiKRCr0HLzF4ws1N7m6rJzIrM7DQzey511ct+F79vJIUJvYHv1Mf44xsacCwikgp9XWmdBZwOrAwXffydmf0kfH4QeA/4KHB2GuqZtWZVRvn0zuVdyn7x8mY2NWuhSBGRodZr0HL31939FGB34DqgERgLNADXAru5++nu/kZaaprFLpxTwcho5+VWTYvzi5c14FhEZKj1u56Wu68mCFrSi7ElBXx9zwoumd85VuuPb9Tz2V3KmV6R7JJlIiLSHyViDJFzdx3B1PLOAcctcbh0vgYci4gMJQWtIVJaaHxvn67TO922tJEX17VkqEYiIvlHQWsInbZDKXuOjnYp+/4LtbhWOBYRGRIKWkMoEg44TvTsmhb+9d6wnKJRRGTI9Rm0zOzybtvndNu+NRWVymUfmlzMh6d2nX/whkUNGaqNiEh+6e9K6+xu2//bbfuIoatK/vjaHhVdtu9f3sSWVo3bEhHZVv0Fre4zvw5qJlgzKzGz583sZTN7zcwuCctnmNlzZlZtZje2z75hZmeb2TozWxA+PjOY82bK3AlFTC7r/NE2tDkPrWjKYI1ERPJDf0GrewbBYDMKmoF57r4XMAc4yszmApcBv3b32cAmILH78UZ3nxM+/jTI82ZExIyTZnSdNPe2pY0Zqo2ISP7oL2gVmtmhZjYvXF+r+3ZBP8cD4IH2Cfmi4cOBecAtYfnfgI8MvAnZ6eQZpV22H1rRRF2LughFRLZFf9M1rAWuSdje0G17bbInMrMCYD4wC7gSWAzUuHtbuMsKYErCIR81s4OBt4Gvu/vyZM+VDfYZG2X7EQW8Ux8DoDkW3Ns6fQctWyIiMliW7jFEZlYF3A5cBPzF3WeF5dOA+9x9DzMbA9S7e7OZnQuc5u5dVlKura3tqHh1dXX6GjAAVy6L8tcVneO2DhoV49e7Kf1dRKQvs2fP7nhdWVnZJZdiwBPjmdlOwK7Af939nYEe7+41ZvY4MBeoMrPC8GprKrAy3GdDwiF/JLj31avEBlZXV3fZzqRzxrTy1xWdF6P/qSlg7HY7MKp4aIbHZVNbU01tzU9qa35KZVv7G6f1SzP7ZML2mcBrwP8Bb5rZ0cmcxMzGhVdYmFkpcDjwBvAYcEq421nAneE+kxIOPyHcN+fsPqqQHSs7/y5oc7j7HSVkiIgMVn9/8n8EeDJh+yfAee4+DjgX+EGS55kEPGZmC4EXgIfd/R7gQuB8M1sEjAH+HO5/Xpga/zJwHjm6ZpeZcVK3hAxlEYqIDF5/3YPj3P1dADPbna6B5e/Ar5M5ibsvBPbuoXwJsH8P5d8Gvp3MZ2e7k2eUctmCzrW1nlzVzLrGGONKk0q8FBGRBP1dadWa2YTw9QeBF929PZMgyiAHGw8nO1VF2W1U598GcYc7l+lqS0RkMPoLWjcBN5jZecC3gOsT3juAIG1d+vHRmRpoLCIyFPoLWt8CHieYY/D/gD8kvDcnLJN+dB9o/OyaFt7bEstQbUREclef97TcvRW4pJf3fpuSGuWh6RWF7DM2yn/XtwLBVCB3LGvkS7uNyGzFRERyTJ9BK0xx75O7Xzt01clfJ88o7QhaALcvbVDQEhEZoP6yB/8KLAJW03PShQMKWkn4yPRSvvdCXcf2i+taWba5jekVAx7fLSIybPV3T+tyoAzYDFwBHO7uH0x4HJzyGuaJqSMKef+Eoi5ldyghQ0RkQPoMWu7+NWB74PfAycAyM/ujmR2Ujsrlm5Omd03IuFVBS0RkQPqdBM/dY+5+r7ufDuxEsO7V42Z2aMprl2dOnF5KJKGT9ZWNrVTXtvZ+gIiIdJHUzK1mVmlmnwceAE4CfggsSGXF8tGEsgIOmljcpUxjtkREktffhLnHmdnNBBPW7g38P3ef7e6XuPumtNQwz3y0+1yESxpJ9/IwIiK5qr/UtbuAt4B/AI3AkWZ2ZOIO7n5RiuqWl47fvoRvPBvM+A7wVm0br29qY7fR0b4PFBGRfoPWtQRp7WPTUJdhYXRJAYdOLubh9zoXg7xtaQO7ja7MYK1ERHJDfzNinJ2megwrJ88s6xa0GvnePiMx0/zDIiJ9GfQSuma2Z3i/SwbomO1KKEr4yS/dHGPBBmURioj0p79EjDIz+6GZ3W1mvzKzkWY208xuB/4NrO3reOlZZVGEI6aWdCm7dYmyCEVE+tPfldaVwPHA68DhwK3AE8BrwHR3/1Jqq5e/us/8fseyRuLKIhQR6VN/iRhHAnPcfa2Z/Q54F/iQuz+V+qrlt6OmlVBWaDSEaYQrtsR4fm0LcycU93OkiMjw1V/QGuHuawHcfYWZ1StgDY3yaIQjp5Zwe8IqxrctbVTQEpGsVd8aZ21jnDWNMdY0BM9rG2OsaYyzpiF4/tyu5RyQwjr0F7QKw+maOtLaum+7+6MpqlveO3lmaZegdeeyRn66fyUFEWURigxH65tiPL6ymbdr2ygvNEYVRxjd/igJnkcVR4j28TvC3alt8SCwdASTWEewWdsYpzHs4fEuxyW8Tngn7rCpOQhW9W3938JYWtfGASMH3PSk9Re01gLXJGxv6LbtwMyhrtRwccSUEiqixubW4B/CmsY4z6xp4eBJutoSGQ6a2pzn1jbz2MpmHn2vmYUbk8siHhkNA1oYyKIRY10YpNY2xmiJp7jifVjdGIdMBS13n566U0tJoXHMdiXcuDihi3BJg4KWSJ5yd96oaePR95p4bGUz/17dQmNs4AlYda1OXWuMd+pjKajltlnbkNo6aQXCDPvojLIuQevOdxr5+dwqigrURSiSD1Y3xHhiVTN3vV3E/PmrgyuRHFUUgfGlBUwojTChLHgeX1rAxNICxodlU8oLqH8vdVPTKmhl2CGTi6kqMmpagr+2NjU7l79azwV7VWS4ZiIyGDXNcZ5Z3cwTq5p5clUzb9a0he8UAv0HrL3GRDtWg9jYHGdjc5xNTcHzhuYYNc1Of9dm5YXWEVjGh4FlQntgKS2gImokTsBj3Z6BLu+PLAqOqyqypGbuqe53j8FT0MqwogLj1Jll/PHNLR1lP19Qx3Hbl7BzlSbRFUmH5pjTHPPwl/nAejka25zn1wZB6omVzby0oZX4AHr8ppQVcOiUYg6dXMyHJhcztqSgz/1jcae2Jd4R0DY2x2mO0RGQxpdGGBEd9GRHWS8tQcvMSoAngeLwnLe4+w/MbAZwAzAa+C/wKXdvMbNigsl69yVI/jjd3Zelo66ZcOHeFdy2tJENzcFfYS1x+MrTm3jgmHE5kUm4YH0LCza0csTUEqaU9/0fTiRTalviLK1rY9nmGEs3t7F0cxtLwu33tsRwoMCgqijI0BtVHCQ7VBV3Zu2NCt+rKDJe29jGE6uaeW5tM80DuI1TXmgcNLGIQ6eUcOjkYnasLBxQoCyIGKNLChjdT3DLV+m60moG5rl7vZlFgafN7H7gfODX7n6DmV0NnANcFT5vcvdZZnYGcBlweprqmnZjSwr4+dxKznmisx/4hXWt/OGNLXxxtxEZrFn/rnt7C1/9dw1xh7JC46oPjuLE6aX9HyiSIq1x5/GVzbywroWldUFwWloX6/ijsC8xhw3N8aT2TVbEYO8xUfYsaeCje0xi/3FFume9DdIStDxY5bA+3IyGDwfmAR8Py/8GXEwQtE4MXwPcAlxhZuZ5vFriyTNKuWVJI/cvb+oo++H8Oo6eVsKMkdnZi/vg8ia+FgYsgIY256zHNvKtORV8c04FEc1aL2kSd+fZNS3cuqSRO5Y1snEIg85g7FJVyMGTijl4UjEfmFhMVXGE6upqZk9UZvC2SttvQzMrAOYDswjmNFwM1Lh7+13KFcCU8PUUYDmAu7eZWS0wBlifrvqmm5nxqwOreOb2NdSFSRmNMee8ZzZx11Fjs27ZkhfXtXD2YxvpKVv3Zws280ZNK78/aBTledy3Lpnl7ryysZVbljRy65JG3tuGVGsDShOmVRuoaSMK+NCkYj4UBqoJZcOz6y4dLN0XL2ZWBdwOXAT8xd1nheXTgPvcfQ8zew040t1XhO8tBvZ39w3tn1NbW9tR8erqVOaqpNedqwv40aKuf419e1YzJ0/MnvEY7zQa57xcQm1b34F0x/I4v9ylmYkleXuBnPea47ClDbbEjPoYbGmz4HVYtiUG9W3hc8woK3BmlTk7jYgzqyxOeQr+LF7RaDy4roAH1xWytDH5P4oKzZlc4kwrcaaWxJlS4kwtjTO1JCgvjgT3k+vaoK7NqGs1atusY7u2NXhd22ZsbjOqos4+lTH2rww+K8v+rsxps2fP7nhdWVnZ5Seb9n4nd68xs8eBuUCVmRWGV1tTgZXhbiuAacAKMysEKoGNvX1mYgOrq6u7bOeab8xynm7YwOMrOxeJvOKdEj65z4Stkhwy0dY1DTHOv3cdtW1dg+i5u5Zz8+LGLvcC3t4S4dOvlnPdvNG8fxvnVMz173UgMt3WmuY411Vv4a9vbWFx3bb9sTSjooA9RkeDx5goe4wuYnJZpKPnoK+2NrTF2dgUZ1OLs7EpzuubWrl1aQMvrut/1ohRxcYJ25eyz9gipldNlJW5AAAVTElEQVQUMmNkAVPKCjKa2JTp7zWdUtnWdGUPjgNaw4BVSrDMyWXAY8ApBBmEZwF3hofcFW4/G77/aD7fz0pkZvz2wCoOvGMtW8Kuis2tzvn/3sQNh4/JaDdhXUucUx/ewLvdRuFfsGcF39t3JF/YdQQfe2QDr29q63hvfVOcEx5Yzy/fX8WZO5anu8oyAItr27j6jXqur27o+Le3rZZujrF0c4y73um8Vzu6OMIeo6PsPjpKfW0UW7upS/p2TfjcNMB4WRbOMHPKzFLmTS5RskOeSteV1iTgb+F9rQhwk7vfY2avAzeY2Y+Al4A/h/v/GbjOzBYRXGGdkaZ6ZoXtKwq5aN+RXPhcbUfZgyuauXlJI6ftUJaROrXEnDMf27jV3GifmF3Gd/cJBkJvX1HIQ8eO49wnN3HPu52/pFrjcN4zNby+qZUf7VdJYQ6k8Q8X7s6Tq5r5/etbeGh5U7+DVofCxuZ4MKZpVTNBTlbDoD+r0OCwqSWcOrOUo6eV6B7qMJCu7MGFwN49lC8B9u+hvAk4NQ1Vy1qf3aWc25c28p+1LR1lFz5XwyGTixlfmt6bvHF3vvz0pi5dlgBHTCnmNwdWdbn6GxGNcO280fxswWZ+vmBzl/2vfn0Lb9a08ZdDRjOquOdfLq1x5+2aNl7Z2MqrG1t5dVMrb9e0Uk4Jv61o5kBlXw2Jxjbn5iUNXP1aPa/XtPW6XyQctzSyyKiIRqiIGhVFESrD54qoMTJ8HhGNsKYhxisbW3llYyuL6toGNMh2IA6cUMSpM8s4cXrJsB2vNFxlZy61EDHjdwdVcdCdazsGLm5qdr75n1r+eujotNbl4hfruGlJY5eyfcZG+euho3tcIiFixnf2HsmuVVG+8NSmLhOCPr6ymcPuXss/Dx/D+NKCjuDU/vxWTWsvM1RHOOXhDdxx5Bj2H6/ABcFV0ootMWIO0YhRXBA+R4yiAnoccrC6Icaf3tzCX97c0udYpImlET67ywjO3qmMMYMMCg1tcd7Y1NYRxF7ZEPwRkmyGXjRCx7Ico4ojjCmJsN+4Ik6aUcq0EfrVNVzpm89isyujfHvOSC6eX9dRdseyRu5+p5Hjt0/PAN6rXqvn8lfru5TNrCjgpiPG9NsV85EZpcwYWcAnHtnIii2dNyiWbI5x4B1re0yX70tDm3Pqwxu49+hx7D56+E5xtaSujRsWN3DjooY+Z/kutGCasKJI+7OxpjFGax9DmPYeG+WLu47gxOml23xPqKwwwr7jith3XFFHWSzuLN0cBLI3a9qo2biB2ZPHdc44kbDcRnnhwKdUkvynoJXlvrz7CO5Y1siCDZ33ki54toYPpqGb7PalDXzn+douZeNKItx25Nh+50drt9eYIh49fhyfenQjzyV0dQ5iNQYAaluckx9azwPHjGNmlg66ToVNzXFuX9rIDYsaeH5dS/8HAG0ObW0e3jHq/QceMThh+1K+sGs5+48vSmmgKIgYsyqjzKoM/uiorl7D7NnZPeuLZJfh878+RxVGjCsOGsUhd62lvVdlTWOc7zxfy9cnpu68T61q5vNPburyq25EoXHzEWOYXjGwfzbjSwu466ixfOPZGv5e3f9N90llEXYfFWSX7TE6ysKNrfzmlc6rvbWNcT7yYBC4JufxXIctMefhFU3csLiBB5c3DfnCfpVFxlk7lvOZXcrZTt1tkiP0LzUH7D46yvl7VXRJbLh+UQNziyOkYiTEaxtb+cSjG7r8kiw0uHbeaOaMLer9wD4UFxi/+0AVu42KcvH8WppjwT2LHSsLO4JTexp093soJ81wVq3fxI2rOrsE362PcfJD67nv6LFZeyN+dUOM6xc1sLSujZHhRKtVRUZVcaRjUtaqoghVxUZlUYTCiOHuvLY5wv89W8OtS/ufjqi4ACaWFtASd1piBM9x73MC11kjCzl313LOmFWW17OBS35S0MoRF+xZwd3LGnkjIdPrJ4uKOGnvOBVD+ItnSV0bpzy8vmMqqXZXHDSKeVNKtumzzYwv7DaC03YoZUNTnOkVhUndNzEzzp/ZCmUjuyyY+WZNG6c+vIE7jho7pD+DbdXY5vz+tXp+tXDzgMY7jYwaRQXG+qYSYEuf+x44oYgzZpVxwvalVPWQienutHmw5EZrPHhuiTuFZkxKGNwrkmsUtHJEUUHQTXjEves60ohXN0f4xrM1XHXQqCEZ6b+4to3jHljHqoauf91fvO9Izpg1dOPDxpQUDDgjLWJB4Kxr8S6TCs9f38onHtnITYePoaQws7+I3Z07ljVy0Yt1LB/EMuh1rQ6tvQe5HUYWcPoOZZy2Q1m/XbRmRtToMbtTJJdlz5+n0q99xxXxpW5Lldy0uJFPP76Rpm2cwWBRbWuPAetzu5Tz1T2y40Z5NGL85ZDRfHBi1y7KJ1c1c84TG2lL1aCgJLy0voWj71vPpx/fNKiA1ZuqIuOcnct5+NhxvHjyBL45Z+SA7ymK5BP9688x39l7JPe929hlTri73mliw8Pr+ce8MT12FfWnuraV4+9fz+rGrgHr47PK+On+lVnVlVRSaFx/+BhOeGA9L63vzKi8990mznumhisOqkrrkiirGmJcOr+Ofy7qOcFkVLHxhV1HUFxgbAqnKKpp8eB1S/t2nNqE7thCc46cVsoZs8r48NQSijUdkUgHBa0cU1poXH/YGE5+cEOXpRieWd3CMfev49YPj2XSAJZFeLumlRMe2DpgfWJ2Gb/7QHoDQLIqohFuOWIMx9y3nrdqO+/xXb+ogcoi4ydpCLSNbc4Vr27m16/U9zhYttCCWU0unDMyqT8kYnGnrjVYRr3uvaXsufPUVFRbJOepezAH7VQV5cFjxzKjtGugeX1TGx++dx3Vtf3Pgg1BwDq+h4D1qSwOWO3GlBRw25FjmTaia4C+6vUt/O/Lm3s5atu5O7cuaWC/29bw45c29xiwjpxWwrMnjeenB1QlfeVbEAmWdp9eUUiaZ+kSySkKWjlq6ohC/rhnE/uP63p/Z3l9jCPvXc+L/QxAfaumleMeWM+abgHrzB3L+G2WB6x2U8oLuPPIsYwv7frP+Ccvbebq1+sZyoUBNjTF+POb9cy7Zx3nPLGpywwf7XapKuT2D4/hxsPHMLty+M7YIZJKClo5rDIKdxw1hqOmdU1F39gcLAfyUEKWXaI3a1o57v71rO0WsM7asYzfHJgbAavdzJGF3PrhsVQWda3zt56rZe9b1/Dt52p4clUzrYNI0mhoi3PrkgZO/9cGdrphNd94trbLfbR2Y4oj/Or9VTx14ngO3cZhASLSN93TynFlhRH+Pm80X/93DdclzDbR0OZ87JEN/O4DVXx8duc6Vm9sCu5hrWvqGrA+vVMZv3x/bgWsdnuMjnLT4WP4yIMbukzOu2xzjKte38JVr2+hssj48NQSjppWwmFTSnrttmuLB0t13Li4gXvfaaK+j6zMaAQ+v8sILtirYlAJMCIycApaeaAwYlz+gSomlBXwi4T7OTGHLz5dw9rGOF/dYwRv1LRxwgPrWd8tYJ2zczn/O7cyJwNWuwMmFPP3w0Zzxr829DghbG2Lc/OSRm5e0kihwQcmFnP0dkEQ235EAS+tb+WmJQ3ctrRxqyvQnhyzXQk/fF8lO1Tqv5BIOul/XJ4wM763z0gmlEb45n9qu8wZePH8Ot6ubeOhFU1bBazPhAErm9LaB+uwKSU8cMw4LltQxxOrmnudyqjN6ViE8FvP1TKmONLnMh3tdqws5NSZpZyaxOBeEUkN/c/LM5/dZQTjSwv47BMbu8wdeH0P44g+u3M5P8+TgNVu33FF3HTEWOpb4zy2spn7323iweVNfQal/taV+ujMMk6dWcpeY6J59bMSyUUKWnnoxOmljC4eyyce2RBMDdSDz+1SzmUH5FfASjQiGuH47Us5fvtSYnHnhXUt3P9uE/cvb+Lt2t5X6gWoiBonTC/ltJmlHDSxeEimyBKRoaGglac+OKmYe48ZxykPbZ3W/vldyvlZHges7goixtwJxcydUMwl+1WyqLaV+5c3cf+7TfxnbQtxD5IqjphawmkzyzhyWgmlGZ7HUER6pqCVx/YYHeWhY8dxysMbqA6vLr64Wzk/3m/4BKyezKqM8pXKKF/ZvYKNTTFWN8aZXFagDECRHKCglee2ryjkqRPG8+CKJiaWRjhgQupXPM4lo0sKsnY9LhHZmoLWMFBSaJw4vTTT1RAR2WbqDxERkZyhoCUiIjkjLUHLzKaZ2WNm9oaZvWZmXw3L9zKzZ83sFTO728xGhuXTzazRzBaEj6vTUU8REclu6bqn1QZ8w93/a2YVwHwzexj4E3CBuz9hZv8D/D/g++Exi919TprqJyIiOSAtV1ruvsrd/xu+3gy8AUwBdgKeDHd7GPhoOuojIiK5yYZyzaGkTmg2nSBQ7Q48AFzm7nea2fnAJe5eEe7zGvA2UAd8z92fSvyc2tra9FZcRETSrrKyssug0rQmYpjZCOBW4GvuXgf8D/AlM5sPVADtKxeuArZz972B84Hr2+93iYjI8JW2cVpmFiUIWP9w99sA3P1N4MPh+zsCx4blzUBz+Hq+mS0GdgReTFd9RUQk+6QlaFkwZ9CfgTfc/VcJ5ePdfa2ZRYDvAVeH5eOAje4eM7OZwGxgSeJndr9kFBGR/JeuK60PAJ8CXjGzBWHZd4DZZvalcPs24C/h64OBS82sDYgB57r7xjTVVUREspW7580DOAp4C1gEfCvT9UlB+5YBrwALgBfDstEEmZfV4fOoTNdzkG27BlgLvJpQ1mPbAAMuD7/nhcA+ma7/NrbzYuC98HtdAByT8N63w3a+BRyZ6foPsK3TgMcIsoVfA76ax99rb23Nu+8WKAGeB14O23pJWD4DeC78Xm8EisLy4nB7Ufj+9G06f6Z/AEP4gywAFgMzgaLwB7prpus1xG1cBoztVvbz9gANfIsgGzPjdR1E2w4G9un2y7zHtgHHAPeHv+TmAs9luv7b2M6LCcYrdt931/DfcXH4C2ExUJDpNgygrZPaAw9BotXbYZvy8Xvtra15992G38+I8HU0DERzgZuAM8Lyq4EvhK+/CFwdvj4DuHFbzp9P0zjtDyxy9yXu3gLcAJyY4Tqlw4nA38LXfwM+ksG6DJq7Pwl07wLurW0nAtd64D9AlZlNSk9Nt00v7ezNicAN7t7s7ksJ/lLdP2WVG2Le+/jMfPxee2trb3L2uw2/n/pwMxo+HJgH3BKWd/9e27/vW4DDbBvWRsqnoDUFWJ6wvYK+/9HkIgceMrP5Zva5sGyCu6+C4D8OMD5jtRt6vbUtH7/rL5vZQjO7xsxGhWV5085w7OXeBH+V5/X32q2tkIffrZkVhPkJawm6eBcDNe7evix4Yns62hq+XwuMGey58ylo9RS5820A8gfcfR/gaILxbQdnukIZkm/f9VXADsAcgjGKvwzL86KdPYzP7HXXHspyqr09tDUvv1t3j3kwzd5UgivEXXraLXwe0rbmU9BaQXAztN1UYGWG6pIS7r4yfF4L3E7wj2VNexdK+Lw2czUccr21La++a3dfE/4SiAN/pLObKOfb2dP4TPL0e+1lLGrefrcA7l4DPE5wT6vKzNoz0hPb09HW8P1Kku8i30o+Ba0XCFLoZ5hZEcENv7syXKchY2bl4WTDmFk5waDsVwnaeFa421nAnZmpYUr01ra7gDMtMBeobe9uykXd7tucRPC9QtDOM8ys2MxmEIxXfD7d9Rus3sZnkoffax9jUfPuuzWzcWZWFb4uBQ4nuIf3GHBKuFv377X9+z4FeNTDrIxByXQmylA+CLKP3iboX/1upuszxG2bSZBt1J5m+t2wfAzwCEGa6SPA6EzXdZDt+ydB90krwV9m5/TWNoLuhivD7/kV4H2Zrv82tvO6sB0Lw//gkxL2/27YzreAozNd/wG29SCCbqCFJKR85+n32ltb8+67BfYEXgrb9CpwUVg+kyDwLgJuBorD8pJwe1H4/sxtOX/aJ8wVEREZrHzqHhQRkTynoCUiIjlDQUtERHKGgpaIiOQMBS0REckZCloiIpIzFLREsoiZ3W9mZ/W/p8jwpHFaIhliZhcDs9z9k2k413RgKRD1zklNRXKOrrREpF8Jc8qJZJSClkg3ZrbMzC4Il5OoNbMbzaykn2OOM7MFZlZjZv82sz0T3rvQzN4zs81m9paZHWZmRwHfAU43s3ozeznc93Ez+0z4+mwze8bMfh1+7hIzOzAsX25maxO7Es3sWDN7yczqwvcvTqjik+FzTXi+95tZxMy+Z2bvhJ91rZlVhp813czczM4xs3eBR82sxMz+bmYbwvq8YGYThuJnLpIsBS2Rnp0GHEWwquyewNm97Whm+wDXAJ8nmFfvD8Bd4WSoOwFfBvZz9wrgSGCZuz8A/IRgFdcR7r5XLx9/AMEcb2OA6wkWN90PmAV8ErgiXA4DYAtwJlAFHAt8wczaF+JrX8amKjzfs2GbzgYOJZg3bgRwRbfzf4hg2YkjCSY9rSSYsXsMcC7Q2NvPRSQVFLREena5u690943A3QTrIfXms8Af3P05D5ah+BvQTLBcQ4xgSfVdzSzq7svcffEA6rHU3f/i7jHgRoKAcakHK94+BLQQBDDc/XF3f8Xd4+6+kGBy3g/18dmfAH7lwWrf9cC3CWYeT+wKvNjdt7h7I8Ekv2MI7sPF3H2+970+lsiQU9AS6dnqhNcNBFchvdke+EbYZVZjZjUEwWWyuy8CvgZcDKw1sxvMbPIA6rEm4XUjBGs0dSsbAWBmB5jZY2a2zsxqCa6Exvbx2ZOBdxK23wEKgcQuv8TVda8DHgRuMLOVZvbzcA0pkbRR0BLZdsuBH7t7VcKjzN3/CeDu17v7QQTBzYHLwuOGOnX3eoLlL6a5eyVwNZ2rxvZ0rpVhndptB7TRNVB2HOfure5+ibvvChwIHEfQHSmSNgpaItvuj8C54ZWOhQt2HmtmFWa2k5nNM7NioIngyigWHrcGmG5mQ/X/sALY6O5NZrY/8PGE99YBcYJ7V+3+CXw9XDh1BJ332HpMiTezQ81sDzMrAOoIugtjPe0rkioKWiLbyN1fJLivdQWwiWCxu7PDt4uBnwHrCbocxxNkDUKwMB7ABjP77xBU5YvApWa2GbgIuCmhjg3Aj4Fnwi7MuQTJI9cRZBYuJQiqX+nj8ycCtxAErDeAJ4C/D0G9RZKmwcUiIpIzdKUlIiI5Q0FLJAlm9p1wUG73x/2ZrpvIcKLuQRERyRm60hIRkZyhoCUiIjlDQUtERHKGgpaIiOQMBS0REckZ/x/xL7seUuNtYQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot n_estimators (x-axis) versus RMSE (y-axis)\n", "plt.plot(estimator_range, RMSE_scores)\n", "plt.xlabel('n_estimators')\n", "plt.ylabel('RMSE (lower is better)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Tuning max_features\n", "\n", "The other important tuning parameter is **max_features**, which is the number of features that should be considered at each split." ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# list of values to try for max_features\n", "feature_range = range(1, len(feature_cols)+1)\n", "\n", "# list to store the average RMSE for each value of max_features\n", "RMSE_scores = []\n", "\n", "# use 10-fold cross-validation with each value of max_features (WARNING: SLOW!)\n", "for feature in feature_range:\n", " rfreg = RandomForestRegressor(n_estimators=150, max_features=feature, random_state=1, n_jobs=-1)\n", " MSE_scores = cross_val_score(rfreg, X, y, cv=10, scoring='neg_mean_squared_error')\n", " RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'RMSE (lower is better)')" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEGCAYAAADWjcoaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XecXGXZ//HPtb1vCElIhSRkk5CEBAhVipIiWFEEFAsgNh4LKvp7LI/lwe6jIiogIAGCYgEUBTsdVFpoIX0TCJDet9eZ6/fHOZudnW2zyU7b/b5fr3nNnPucmXOd3WSuvcu5b3N3REREskFOugMQERFJlJKWiIhkDSUtERHJGkpaIiKSNZS0REQka+SlO4ADVVNTo2GPIiJDXGVlpcVuq6YlIiJZQ0lLRESyhpJWilVXV6c7hKTTNQ4NusahYahdo5KWiIhkDSUtERHJGkpaIiKSNZS0REQkayhpiYhI1lDSEhGRrJG1M2IcLHdn5d52HtjczCt1Ea563Yh0hyQiIv0YlkmrsT3KCb/fwebGyP6yLx1bzuji3DRGJSIi/UlJ86CZFZnZU2b2gpmtNLMrw/IlYdlyM7vLzMrC8kIz+52ZrTezJ81s8mDGU5KXw6jirpf+0JaWwTyFiIgkQar6tFqABe4+DzgGONvMTgY+6+7z3H0u8CrwyfD4DwF73X0a8GPg+4Md0KIJhV2279/UPNinEBGRQZaSpOWB+nAzP3y4u9cCmJkBxUDHzO3nAEvD13cBC8NjBs2iiUVdth/Y3ELUNXG8iEgmM0/RF7WZ5QLPANOAa939C2H5LcCbgVXAW9y90cxWAGe7+6bwmA3ASe6+q+PzYpcmOZC5tdodFj9RTH2kMxfeOq+Z2eXRA7k8EREZJFVVVftfxy9NkrKBGO4eAY4xsxHA3WY2x91XuPsHw4T2M+DdwC1AT7WqXrNr7AUOxMLNu/nTxs5mwXU5o3lHVcUBfVaiqqurDzjebKFrHBp0jUPDULvGlN+n5e77gIeBs2PKIsDvgHeFRZuASQBmlgdUAnsGO5aFE7o2Ed6/SYMxREQyWapGD44Oa1iYWTGwCFhrZtPCMgPeBqwJ33IPcHH4+jzgQU9CO+aiuKS1bFcre1vUPCgikqlS1Tw4DlgaNgPmAHcAfwEeM7MKgubAF4D/Co9fAvzSzNYT1LDek4ygxpfmMuuQPFbtbQcg6vDQ5mbOnVqSjNOJiMhBSknScvflwLE97Dq1l+ObgfOTGlRo8YQiVu2t3799/+YWJS0RkQw17OceXNht6Huzhr6LiGSoYZ+0Th5TQFle52DF7U1RVuxpS2NEIiLSm2GftApyjTPGx82OsVmjCEVEMtGwT1oQ9GvF0pROIiKZSUkLWBA3D+GTO1qpadXQdxGRTKOkBRxRnsf0ys6BlBGHRzTru4hIxlHSCi2aGN+vpSZCEZFMo6QVip8d44FNLaRqMmEREUmMklbodYcVUpzbOfR9c2OENfva0xiRiIjEU9IKFeUZp48r6FKmUYQiIplFSStGfBOh7tcSEcksSlox4lczfnx7C/VtGvouIpIplLRiTK3IY2p57v7t1ig8tlW1LRGRTKGkFSd+Al01EYqIZA4lrTjxUzrdt6lZQ99FRDKEklac08YVUNjZQsir9RHW12rou4hIJlDSilOSl8Oph8XNjrFJTYQiIplASasH3fu1dL+WiEgmUNLqweK4Wd//va2Fpnb1a4mIpJuSVg+qKvOYVNbZsdUcgX9tUxOhiEi6KWn1wMxYNCG+X0tNhCIi6aak1YvuUzopaYmIpJuSVi/OGF9IfsxPZ0NthJc19F1EJK2UtHpRnp/DyWPiZn1XbUtEJK36TFpmlmdm55rZEjNbZmbrw+clZnaemeX19f5sFz+BrqZ0EhFJr16Tlpl9DHgJ+BiwAfg2cFn4vAH4CPCSmV2WgjjTIr5f67GtLTRr6LuISNr0VVOaDpzo7tt62Hc38B0zGwd8LimRZYBZh+QxviSHLY3B8iSN7c4TO1p4w/iift4pIiLJ0GtNy90/5+7bzCzHzBaYWUEPx2x1988nN8T0MTMWxo8i1JROIiJp0+9ADHePAn9y99YUxJNxuvdraTCGiEi6JDp68FEzOzmpkWSo148rJNc6t9fsa+e1eg19FxFJh0RH/70C/M3M/gS8BuwfjeDuX0tGYJliRGEOJ44p4PHtnRXNBza3cMmMIT1wUkQkIyVa0yoG/kiQrCYCk2IeQ1632TE0pZOISFokVF1w9w8mO5BMtnBCId98tnP7ka0ttEacgth2QxERSbqEZ8Qws6PM7Ktmdk24PcPM5iYvtMwx99B8Rhd1/qjq2pyndg7LcSkiImmVUNIys/OBR4EJwEVhcTlwVZLiyig5ZizUrO8iImmXaE3rG8Bid78MiIRlLwDzkhJVBtKUTiIi6Zdo0hpDkKSgc+Sgx7we8haMLyS2B2vFnja2NkZ6PV5ERAZfoknrGeADcWXvAZ4a3HAy18iiXOaPzu9S9oBuNBYRSalEk9blwLfM7BGg1Mz+AXwT+GzSIstA3Ye+q4lQRCSVEkpa7r4GmAlcC3wFuAU42t2rE3m/mRWZ2VNm9oKZrTSzK8Py281srZmtMLObzSw/LDcz+2m4FMpyMzvugK5ukMX3az20pZn26LBpIRURSbtERw/+1N0b3f0Od/+Bu//W3evN7OoEz9MCLHD3ecAxwNnhtFC3EyTDowluYP5wePybgKrw8VHg54lfUvIce2g+Iws7f2Q1rc4yDX0XEUmZRJsHL+mlPL6fq0ceqA8388OHu/tfw31O0D82MTzmHOC2cNcTwIhwGZS0ys0xFsQPfdcoQhGRlOlzRgwzu7TjuJjXHaYCuxI9kZnlEgzomAZc6+5PxuzLJ0iAnw6LJhDMcdhhU1i2NdHzJcuiCUXc9VLT/u37NzXzleMq0hiRiMjw0d80Th01qQK61qoc2A5cnOiJ3D0CHGNmI4C7zWyOu68Id18HPOruj4XbPc2P1GvnUXV1Ql1rg2JyK0DJ/u3nd7fx5MpqRnZbbax3qYw3XXSNQ4OucWjItmusqqrqdV+fScvdzwQws2+5+1cGIxh332dmDwNnAyvM7OvAaOBjMYdtoutkvBOBLb19Zl8XONiqgGM27OD53W37yzYWTeCkI0t6f1OM6urqlMabDrrGoUHXODQMtWtMtE/rLT0VmtmyRN5sZqPDGhZmVgwsAtaY2YeBs4ALw8UmO9wDXBSOIjwZqHH3tDcNdtCs7yIi6ZFo0joyvsDMjKBfKxHjgIfMbDnwNHCfu/8ZuB44DHjczJ43s461uf4KvASsB34BfDzB86TEwoldB2M8sLmFiIa+i4gkXX8DMW4LXxbGvO4wGViZyEncfTlwbA/lPZ4/HE34iUQ+Ox1OGF1ARYFR2xokqj0tUZ7f3cb80QPo2BIRkQHrr6a1IXzEvt5AUAO6nWBo+rCTl2OcOT5+6LuaCEVEkq2/gRgdM1c84e7/SE1I2WHRhCL+tLEzUd2/qZkvHKOh7yIiyZToNE7/MLPFZrbEzO4FMLPjzWxBcsPLXAvjBmM8s6uNPc2a9V1EJJkSncbpUwRTKVUDZ4TFTcC3khRXxhtfmsvsQzorqlGHh7ZodgwRkWRKdPTgZ4BF7v49oGNo+hpgRlKiyhLdhr5rSicRkaRKNGmV0zmtUsfY7nxgWM8WGz/r+wObm4m6hr6LiCRLoknrUeCLcWWXAw8NbjjZ5aQxBZTldc44taMpyot72vp4h4iIHIxEk9angHea2Uag3MzWAucDVyQrsGxQkGu8Pn7ouxaGFBFJmkRHD24FTgAuAN5LMFHuSe6+LYmxZYXu/Vq6X0tEJFkSrWl1HJsfvs6l55nYh534KZ2e2tHKvpZoL0eLiMjBSHTI+1yC4e53AP8PuBOoNrN5SYwtKxxelseMys6h7xGHR7aqiVBEJBkSrWndDFwLTHT3EwkWZLwmLB/24mtbmvVdRCQ5Ek1a04Grw4lsOya0/QnB8lLD3uIJ3Ye+u4a+i4gMukST1l+Bt8eVvQ34y+CGk51OOayQkpih71sao6ze157GiEREhqZeJ8w1s1/SeSNxLvBbM3uG4CbjScB84E9JjzALFOUZp48t4B8xw93v39TMrEPy+3iXiIgMVF+zvK+P214R83oVoFnfYyycUNQ1aW1u4fKjy9MYkYjI0NNr0upYlkQSs3hiETxZs3/78e0t1LVFKc8fyF0FIiLSF32jDpIpFXlMLc/dv90Whcc09F1EZFApaQ2i+Al0NaWTiMjgUtIaRD1N6aSh7yIig0dJaxCdNq6Aws4WQl6tj1Bdo6HvIiKDJdFpnC40s6PC1zPM7FEze9DMZiY3vOxSkpfDqYfFzY6hhSFFRAZNojWtbwF7wtc/BJ4iWGPrumQElc2692tpSicRkcHS131asUa7+3YzKwJOA84D2oBdSYssSy2aUMiXY7b/vb2FxvYoJXlqiRUROViJfpPuNLNpwJuAp929BShCy5N0U1WZx6Syzo6tlgj8a2trGiMSERk6Ek1a3wSeAZYAPwjLFgIvJCOobGZm3SbQ1cKQIiKDI9GVi28FxhEsTXJfWPwk8J4kxZXVFk7QUiUiIsnQa9IyM4t5nQM0A81mlhNu7wJ2JD/E7HPG+EJiZ296qS7CS7Ua+i4icrD6qmnVxLxuJxh4EfvoKJM45fk5nDymoEuZalsiIgevr6Q1O+b1FGBq3KOjTHqweGL3hSFFROTg9DXL+2sxr19JTThDx8IJRXxtWe3+7ce2tdLcrimdREQOhm4eSpJZh+QxvqTzx9vY7jy+XbNjiIgcDCWtJDEzFnYb+q6kJSJyMJS0kkhTOomIDK5+k5aZ5ZrZBjMr7O9Y6er14wrJjZkzZG1NO1ubNYmIiMiB6jdpuXsEiBBM2yQDMKIwhxPjhr4/vje3l6NFRKQ/iTYPXg3cYWavN7MjzWxqxyOZwQ0F8QtDPr5XLbIiIgcq0VnerwmfF8eVO6CqQx8WTijkm892bj9Vk0trxCnIVTOhiMhAJTr3YE4vDyWsfsw9NJ8xxTFD3yPGkzs067uIyIEYUFuVmU0ys5MHehIzKzKzp8zsBTNbaWZXhuWfNLP1ZuZmNirmeDOzn4b7lpvZcQM9Z6bIMWPBeE2gKyIyGBJKWmZ2uJn9G1gD3B+WnWdmNyV4nhZggbvPA44Bzg6T37+BRUD8jBtvAqrCx0eBnyd4nowUP6WTlioRETkwida0bgD+ApTTOUnufXTv4+qRB+rDzfzw4e7+nLtv7OEt5wC3he97AhhhZuMSjDXjnDm+kJyYLqyVe9t5rV6zvouIDFSiSetE4HvuHiUYfIG71wCViZ4ovN/reYLlTO5z9yf7OHwC8FrM9qawLCuNLMpl/qj8LmXXr2pIUzQiItkr0dGD24FpwLqOAjObBbya6InC+72OMbMRwN1mNsfdV/RyeE9D63qdbba6ujrRMNJmYUUeT+/svGdryeo63lG2gxH5fbwpi2XD7+Rg6RqHBl1j5qmqqup1X6JJ64fAn83su0CemV0IfBn43kCDcfd9ZvYwcDbQW9LaBEyK2Z4IbOntM/u6wEzxmSnOLVu2sb0pCkBz1Liv5TC+PKsizZENvurq6qz4nRwMXePQoGvMPokOeb8Z+G/gfIJmu4uAr7r77Ym838xGhzUszKyYYPDFmj7ecg9wUTiK8GSgxt23JnKuTFWUZ3xydlmXshtW1VPXFk1TRCIi2SfhIe/u/kd3f7O7z3b3N7n7HwdwnnHAQ2a2HHiaoE/rz2Z2uZltIqhJLY8ZjfhX4CVgPfAL4OMDOFfGumRmKRV5na2cNa3OrWvUtyUikqiEmgfN7DngYeAR4BF33zuQk7j7cuDYHsp/Cvy0h3IHPjGQc2SD8vwcLhjXzk2vdXZkXbOyno8cVUZRnmbIEBHpT6I1rc8DtcBngM3hDb8/M7Pzkhfa0PTu8W2UxiSo7U1Rfr2+MY0RiYhkj0T7tB5w96+7+xsIBkjcS9Cv9bskxjYkjciHS2aUdim7+sU62qO9Do4UEZFQojNinG1m3zWz/wDPE8xU8SVgbjKDG6o+MbuMgpif/Kv1EX7/clP6AhIRyRKJNg/+FTgPuAmY4u4XuPt17r4yeaENXeNLc3nvtJIuZVcvryPqqm2JiPQl0aR1BnAzwZD3V83sn2b2P2Z2evJCG9o+fXR5l6mdVu9r52+vak5CEZG+JNqn9S93/667v4lgwtunCe7bejiJsQ1pUyryOHdKcZeyq5bX4aptiYj0KtE+rXea2U/CuQM3EtS8riGYjV0O0GeOLu+y/cyuNh7d2pKmaEREMl+izYOfBvYBVwCHuvvp7v4/7v7P5IU29M0Zmc/Zk7ouW/Kj5fW9HC0iIgndXBwOdZck+Nzccv7+Wmdf1qNbW1i2s5XjRxf08S4RkeEp0ebBfDO70sxeMrPm8PlKM9M360E6YUwBp43t+mO8anldmqIREclsiTYP/h/BJLeXAfPC5wXA95MU17Dyubld+7b++mozq/a29XK0iMjwlWjSOh94u7v/093Xhn1Z7wQuSF5ow8cbxhdybNwikVertiUi0k2iSau32Vw1y+sgMDOuiKtt3fVyExvr2tMUkYhIZko0ad0J3GtmZ5nZUWZ2NvBH4I7khTa8vOXwImZUdo6LiTr89EWNJBQRiZVo0vpv4H7gWuAZgnu0HgK+kKS4hp0cMz4bV9v6VXUD2xojaYpIRCTzJDojRqu7f83dp7l7Sfj8VXfXnbCD6F1Ti5lUlrt/uzUK165UbUtEpEOv92mZ2YJEPsDdHxy8cIa3/Bzj03PK+PwTNfvLbl7TwBVzyzmkMOFFpkVEhqy+bi5eksD7HZg6SLEI8L6qUv7vhTp2NEUBaGh3blhVzxePrUhzZCIi6dfrn+/uPiWBhxLWICvOMz4xu6xL2fWr6qlvi6YpIhGRzKE2pwz0wRmlVBZ03k2wr9W5dW1DGiMSEckMvSYtM3vazM7vbaomMyswswvM7MnkhTc8VRTk8JGjuta2rllRT0tEy5aIyPDWV03rYuDdwJZw0cefmdl3wud/AJuBdwGXpCDOYee/ZpVSktdZ29rWFOU36xvTGJGISPr11ae1yt3PA+YAvwSagFFAI3AbMNvd3+3uq1MS6TBzaFEuF08v6VJ29Yt1tEdV2xKR4avfpUncfRtB0pIU++Sccm5a00DHGIyNdRHufrmJ848s6fuNIiJDlAZiZLAJpblcOK1rgvrx8jqirtqWiAxPSloZ7tNzysmJmZZ41b52/hGzaKSIyHCipJXhjqzM4x2Ti7uU/Wh5Ha7alogMQ0paWSB+It1lO9t4bFtrmqIREUmfPpOWmf00bvtDcdu/T0ZQ0tXRI/M5a2Jhl7KrtEikiAxD/dW0Lonb/kHc9uLBC0X6Er9I5MNbWnh2p2pbIjK89Je04lcm1krFaXLSYYWcOrbr5CSqbYnIcNNf0orv7VfvfxrF17b+/Goza/a1pSkaEZHU6+/m4jwzO5POGlb8dm7Pb5NkWDC+kHmH5vPC7s5E9ePlddxwxsg0RiUikjr9Ja0dwM0x27vjtncMekTSKzPjirnlXPzQnv1ld73UxJeObWdyeb+Tm4iIZL0+v+ncfXKK4pAEve2IIqZX5rGuph2AiAczwP/wlBFpjkxEJPkGfJ+Wmc0ws3ea2RHJCEj6lmPGZ47uumzJL6sb2N4YSVNEIiKp0999Wj8ys/fHbF8ErARuBNaY2ZuSHJ/04PwjS5hY2tmd2BKB61bWpzEiEZHU6K+m9Q7g0Zjt7wCXu/to4DLg68kKTHqXn2NcPqdrbWvJmgb2tUTTFJGISGr0l7RGu/urAGY2BzgUWBLu+xUwPYmxSR8+ML2U0UWdv776dufG1apticjQ1l/SqjGzw8LXpwPL3L0l3M5HNxunTXGe8fHZXWtb169qoKFNtS0RGbr6S1p3AL81s8uBLwK/jtl3ErAhkZOYWZGZPWVmL5jZSjO7MiyfYmZPmlm1mf3OzArC8sJwe324f/JAL2w4uHRmKRUFnX837GmJsnRdYxojEhFJrv6S1heBhwnmGLwRuCFm3zFhWSJagAXuPi9839lmdjLwfeDH7l4F7AU6JuT9ELDX3acBPw6PkziVBTl8dGbX2tY1K+poiWjiEhEZmvpMWu7e5u5Xuvvb3P3b7h6N2fcTd08oaXmgo8MlP3w4sAC4KyxfSjDwA+CccJtw/0IzU1NkDy6bXUpxbuePZktjlN9tUG1LRIYm62sxwXCIe5/c/baETmSWCzwDTAOuJZgx/omwNoWZTQL+5u5zzGwFcLa7bwr3bQBOcvddHZ9XU1OzP/Dq6upEQhiyfrghn99tzd+/Pakoyp3zm8lVmheRLFRVVbX/dWVlZZdvsv7m/rkVWA9so+dBFw4klLTcPQIcY2YjgLuBo3r5PPo4V49iLzDTVVdXD3q8Xx3Xzu/v2k57+BN6rTmHVXkTOHdqyaCeJ1HJuMZMo2scGnSN2ae/Pq2fAiVAHXANsMjdT495nDHQE7r7PoJ+spOBEWbWkTgnAlvC15uASQDh/kpgD9KjiWV5vGda1wT1o+V19FWLFhHJRv31aX0GOAK4DjgX2GhmvzCz0wZyEjMbHdawMLNiYBGwGngIOC887GLgT+Hre8Jtwv0Pur6B+/SZo8u6VE9X7m3nn5taej1eRCQb9Tv3oLtH3P0v7v5uYAbBKL+HwyVKEjUOeMjMlgNPA/e5+5+BLwBXmNl6ut64vAQ4NCy/gmAUo/RhWmU+50wu7lL2oxdU2xKRoSWh9SzMrBJ4D0HtZzTwTeD5RE/i7suBY3sofwk4sYfyZuD8RD9fAp+dW8YfNzbt335qZyv/3t7KaWML0xiViMjg6W/C3Lea2Z0ETXnHAv/P3avCYfB7UxKhJGzeoQUsntA1QV31Ql2aohERGXz91bTuAdYCtwNNwFlmdlbsAe7+tSTFJgfginnl3Le5sy/rwS0tPL+rlWNGFaQxKhGRwdFfn9ZtwBPAKILRfD09JIOcclghpxzWNUFdtVy1LREZGvpbufiSFMUhg+hzc8s5777d+7fvfaWZBzY3s2B8IZpYRESy2YBXLu5gZnPD/i7JMAsnFDJ3ZOcMGQ6865+7OfPendyxoZFWzU0oIlmqv4EYJWb2TTO718yuMrMKM5tqZncD/wF2pCZMGQgz44q55d3Kn9/dxkcf3cu8u7Zx1fI69jRH0hCdiMiB66+mdS3wNmAVwQ3BvwceAVYCk939E8kNTw7U2ycX8fYjinrct7UxyjeeqWX2Hdv57H/2sm5fW4qjExE5MP2NHjwLOMbdd5jZz4BXgde7+2PJD00ORo4ZS88cyb+2tXLdynr+/lpzt8kbmyLOLWsbuWVtI4snFPLx2WW8Qf1eIpLB+ktaZe6+A8DdN5lZvRJW9jAzTh9XyOnjCtlQ0871q+v5dXUjDe3d+7Tu29zCfZtbmDUij8tml3HB1BKK8pS8RCSz9Nc8mGdmZ5rZAjNbABC73VEmme/Iyjx+cPIIVl4wlm8cX8HE0twej1u1r53L/72POXdu4zvP1bKjSf1eIpI5+qtp7QBujtneHbftwNTBDkqSZ0RhDpcfXc5/zS7j3o1NXLeqnmU7u/dp7WqO8n/P13H18jrOm1rCx2eXMSdmRKKISDr0d5/W5BTFISmWn2OcO7WEc6eW8NSOFq5b2cA9rzQRjWs5bI3Cr9c38uv1jZwxrpCPzy7ljROLyFG/l4ikQUIT5srQduKYQk4cU8ir9e3cuKqB29Y1UNvWvd/r0a0tPLq1hWkVeVw2q5QLp5VQmn/At/qJiAyYvnFkv8PL8vjWiZWsfPdYvndSJZPLe+73Wl/bzuefqGH2Hdv432U1bG5Qv5eIpIaSlnRTnp/DZbPKeObcw/jVgpGcOrbnyXb3tTpXv1jPvDu38eFH9vDsztYURyoiw42aB6VXuTnGW48o5q1HFPP8rlZ+vqqeP7zcRFu063HtDne91MRdLzVx0pgC3jkyl7GTo5Sr6VBEBpmSliTkmFEF3HDGSP73+AhLVjdw89oG9rREux335I5WntxRyJfXbmXmiDyOH13A8aMLmD+qgJkj8sjN0QAOETlwSloyIONKcvnK/AqumFfGHRua+PnKetbWtHc7Luqwam87q/a2c9u6RgDK8oxjR+UHSSxMZmNLeu43ExHpiZKWHJCSvBwumVHKRdNLeHBzC9etrOfBLS19vqe+3XlsWyuPbevs+5pYmhsmsSCZzTs0n5I8NSuKSM+UtOSg5JixaGIRiyYWsXpvGzetaeD+V+p4pSmxxLOpIcKmhib+uLEJgFyD2Yfkc8KYAuaHtbJplXm6L0xEACUtGURHHZLPj04ZQfWonYw+/Eie3dXKsp0dj7Ye+8DiRRyW72lj+Z42loRlFQXG/FFBk+IJows4fnQ+hxapWVFkOFLSkqQYUZjDgglFLJgQLI/i7mysi8QksVZe3NNGa/95jNpW56EtLTwU0/w4uTx3/wCP40cXcER5LpUFORTmqkYmMpQpaUlKmBlTKvKYUpHH+UeWANAScV7c09YlkW2sS+xG5Y11ETbWBcPsYxXmQkV+DhUFRkVBTtxro7IgJ3hdYFTk51DZw3FKfCKZS0lL0qYw1/YPie+wqznCMzs7E9kzu1qpbe0+pVRvWiKwMxJlZzPAgc3UUZRLjwkvNtnl1edyyaQIo9RMKcOcu9PY7uxrdfa1RBlXksPIJP6/UNKSjDKqKJezJuVy1qSgWTHqzvqa9jCBBclsxZ42IonnsQFrjkBzU5QdTX21XRbyg5e3cf7UEj56VClzD+151hCRbODu1LU5+1qj7GuJ7k9A+1qj1LREqWkN9+3fH2VfS1BW0xrtMuHAdaeN4L1VpUmLVUlLMlqOGdNH5DN9RD7vrQrKGtujvLA7SGDP7GzjxT2t7G1xaluj9LC+ZdK0ROBX1Y38qrqR1x1WwMdmlfGWw4vI0w3UA9IWdf7+WjPLdrRy3OgC3np4kW5CPwiRqLOzOcqWhgibGyOs2pZLUXNdt2QTm4BqWr3bCg8Hat8AWkYOhJKWZJ2SvBxOOayQUw4r7FLu7jRFnJrWIIHVtjq1bdHO161RatrC5y77w+NiuhvBAAARAUlEQVTD1wdSi/vP9lb+s30PE0tz+chRpVw0vZRDCnW/WV821rVz27oGflXd2KVWO2tEHl+ZX8GbJhVhutWhi9aIs60pwpaGyP6ktLUxwpaGIEltCbe7/hsuBGpTFuO+REZXHQQlLRkyzIySPKMkL5i540B0tM/XtnmY2GKTWvC8uznKndW1bGnpnpQ2NUT4+rJavvdcHRccWczHZpUx6xAtntmhNeL87bVmbl3b0GU0aKxV+9p57wN7OGF0Pl+dX8kZ4wp7PG6oaWyPsrUhGpOIOhNTR0La2RQlhY0JCSvMhREFOYwoyKGyILl/rClpicQwM0rzjdL8vhPfhZU7WF84kRtW1XeZ4aNDU8RZuq6RpeuCxTM/dlQpZ08avs1eL9W2s3RtA79e38jO5sT+En96Zxtv//suzhxfyNfmV3DsqOzuN6xtjfLMzlY2hQmoIyltaQwee1vSm45K8owRBRYknsIgAY0oDEbYdrwOni1mX1BWnJe6f9dKWiIHINfYPwP+yj1t3Li6nt9taKS5hwGLHYtnHl4WNB1+oKqUEcOg6bAl4vzllSZuXdfIo1v7nuKrPN+YP7qAR7a0dKtJBPfo7eTtRxTxleMqmD4iu2quL+xu5eY1Ddz5UhONKep0HVmYw7iSHCaU5lLU1sDho0b0mYAqC3IoyJJbPZS0RA7S7JH5/OTUQ/j6/ApuW9fITWsa2NTDwpiv1kf46tO1fPe5Oi6cFow6nJFlX8CJWF/TxtJ1jfy6upHd/cyCcvzofC6eXsq5U4opzc9h5Z42vvVsLX97rbnbsfe80syfX23mwmklfPGYciaVZe7XV1O784eXG7llbQPLdrYN2ucaMKY4h/GluYwvye32PKE0l3EluV1qPtXVe6mqqhy0GNItc3/rIllmZFEun5lbzifnlPGXV5u5flU9j2/v3nTY2O4sWdPAkjUNnDm+kMtmlbF4YmFWz6/Y3O7c+0oTS9c18K8emktjVeQb7z6yhItmlHL0yK5Je/bIfH6z6FCe2tHCN56p7fZZUYfbqxu5c0Mjl84s5XNzyxldnDn3yq2vaePmtQ38urpxwKPo8gzGxiSe8aU5+xNRR2IaW5JL/jBtYu6gpCUyyPJyjHMmF3PO5GJe2N3KjasbuOulRlp6aDrsmJ5qSnkuHz2qjPdVlVCR5I7swbR2XxtL1zXwm/WN/fbJnDSmgIunl/COKcX9zuR/4phC7j17FA9tCZLX87u71lZao3D9qgZ+ua6R/5pdxqfmlCV9AEBv2qLOX19t5uY1DTzSTzPo+JIcTh1b2CURdTyPLsoZtn2eA6GkJZJE8w4t4NrTCrjy+AqWrm3kpjX1bG3s3mT2cl2ELz1Vw7efreW9VUHT4bTKzGw6bGp37nmliVvXNvRYk4xVWWC858gSLp5ROuBRlGbGgglFnDm+kHteaeZbz9ZSHbd2W0O788MX6liypp7PHl3OR44qS9mggM0NEZaua+CX6xp6/J3GWjC+kEtnBoNxdB/fwVHSEkmBUUW5fG5eOZcfXca9G5u4YXUDT+7o/oVf3+7cuLqBG1c3sHhCIR+bVcaCCZnRdLh6b1Cr+u36/pu+TjmsgIunl3LO5OKDTiJmQc31LYcX8dsNjXzvubpufYZ7W5yvLavl56vq+e95Fbx/eklSmtGi7jyypYWb1jTw99ea+7yn75BC4/1VpXxwRilTK/RVO1j0kxRJofwc49ypJZw7tYTndrVyw6p6/vByU4+z3d+3uYX7NrdQVZnHR48q5fRxhRTlWpdHYS5JvQG3sT3KH19uYum6xh6TbKxDCo0Lp5Vw0fRSZiZhgEleTpAEzp9aws1rGvjR8jp2xQ2f39oY5bOP7+NnK+r48nEVnDuleFAS/p7mCLevb+SWNQ281M+kzieOLuDSmaW8Y3IxRSkcCj5cmHsm3qrWv5qamqwMvLq6mqqqqnSHkVS6xoHZ0RThlrUN3Lymge19znfYs8JcYpKYUdzxnBc8F8XtL8rrnvg6ju14b14O/GHldv6+u6DfCYtPHVvAJdNLedsRqf2SrmuL8vOV9Vyzop7atp5jnDMyn68eV8EbJxb2mNz7+j26O8t2trFkTT13b2zqsU+yQ2meccGRxVw6s6zb4JJ0y/b/j5WVlV1+cappiaTZmOJcvnBMBZ89upw/bWzi+lX1PLMr8WHSLZHgnqiaQZ8rIR96+cyRhTm8d1oJF00vSdt9U+X5Ofz3MRV8eGYpP36xnl+sru92n9yKPW28+/7dnDymgK/Nr+B1Y/ufXaO+LcqdG5q4eW0DL+7p+/cwa0Qel84s5YIjs2sATTZT0hLJEAW5xvlHlnD+kSUs2xk0Hd79clNKJwHuz+ljC7hkRilvPaI4Y9YdG1mUyzdPqOSyWWX84Plaflnd2K2v6Ykdrbz5b7tYPKGQr8yvYF4Ps/Kv2tvGzWsa+N2GRup6qbkBFOTAOZOLuXRmKSePKdD8iCmWkqRlZpOA24CxQBS40d1/YmbzgOuBMmAj8D53rw3f8yXgQwSLIl3u7v9IRawimeD40QUc//qRfOOEoOnwoc3N1LQ6zRGnJRJMDNwS8T6brAbLqKLOWlWmjmgEmFCay9WnHsKn5pTznedq+f3LTd2OCfoJd3LulGK+fGw5rVG4c0MjNycwEvKIslw+OKOU908v0TpqaZSSPi0zGweMc/dnzawceAZ4B7AU+Ly7P2JmlwJT3P2rZjYL+A1wIjAeuB+Y7u77/4uqTytz6RpTJ+q+v3mwI5E1tQfPzbGPdrolvO7H0aU8t7WB980Zw5sPL8qaKX5iLd/dyreereWfm3q+dyrXoCzXqWnv/dpyDN44sYgPzSxlYYaM4hyoTPm3eqDS0qfl7luBreHrOjNbDUwAZgCPhofdB/wD+CpwDvBbd28BXjaz9QQJ7PFUxCuSLXLMKM4LBlKMGOTPrq7eS9WU4kH+1NSZe2gBdywexX+2tfDNZ2u71aQiTq8Ja0xxDhdVlXLxjJKMni5qOEr56EEzm0yQqOYAfwe+7+5/MrMrgCvdvdzMrgGecPdfhe9ZAvzN3e/q+JzYmlZ1dXUKr0BEso07/GdvDte+UkB1Q+8DJuZXRnjX2HbecGiEfI2rSJvYmmFaRw+aWRnwe+Az7l4bNgn+1My+BtwDdPwp1NOfP71m12yq+mZ7VT0RusahYahd43TgohOdu19u4tvP1u6/36qiwLjwyBIunTk0JzAear/HlCUtM8snSFi3u/sfANx9DfDGcP904C3h4ZuASTFvnwhsSVWsIjI05ZjxrqklvH1yMQ9ubuHVzVt47/zJlKpalTVS8puyYEzoEmC1u18VUz4mfM4BvkIwkhCCWtd7zKzQzKYAVcBTqYhVRIa+/BzjrElFvOHQiBJWlklVTetU4APAi2b2fFj2ZaDKzD4Rbv8BuAXA3Vea2R3AKqAd+ETsyEERERmeUjV68F/03E8F8JNe3vNt4NtJC0pERLKO6sUiIpI1lLRERCRrKGmJiEjW0NIkIiKSseJvLlZNS0REsoaSloiIZI2sbR4UEZHhRzUtERHJGkpaKWJmk8zsITNbbWYrzezT6Y4pWcws18yeM7M/pzuWZDCzEWZ2l5mtCX+fp6Q7psFmZp8N/52uMLPfmFlRumM6WGZ2s5ntMLMVMWUjzew+M6sOnw9JZ4wHq5dr/EH4b3W5md1tZoO9ik1KKWmlTjvwOXc/CjgZ+ES42OVQ9GlgdbqDSKKfAH9395nAPIbYtZrZBOBy4Hh3nwPkAu9Jb1SD4lbg7LiyLwIPuHsV8EC4nc1upfs13gfMcfe5wDrgS6kOajApaaWIu29192fD13UEX3QT0hvV4DOziQSz9d+U7liSwcwqgDMIJoDG3VvdfV96o0qKPKDYzPKAEobAKgvu/iiwJ674HIIV1Amf35HSoAZZT9fo7v909/Zw8wmCVTOylpJWGoQLYR4LPJneSJLiauC/gWi6A0mSqcBO4JawCfQmMytNd1CDyd03Az8EXiVYcbzG3f+Z3qiS5rBwZfWOFdbHpDmeZLsU+Fu6gzgYSlopFr8QZrrjGUxm9lZgh7s/k+5YkigPOA74ubsfCzSQ/U1KXYT9OucAU4DxQKmZvT+9UcnBMrP/IeimuD3dsRwMJa0U6mkhzCHmVODtZrYR+C2wwMx+ld6QBt0mYJO7d9SS7yJIYkPJIuBld9/p7m0Eywa9Ls0xJct2MxsHED7vSHM8SWFmFwNvBd7nWX6fk5JWivS2EOZQ4u5fcveJ7j6ZoOP+QXcfUn+hu/s24DUzmxEWLSRY920oeRU42cxKwn+3Cxlig01i3ANcHL6+GPhTGmNJCjM7G/gC8HZ3b0x3PAdLSSt1OhbCXGBmz4ePN6c7KDkgnwJuN7PlwDHAd9Icz6AKa5F3Ac8CLxJ8T9yY1qAGgZn9BngcmGFmm8zsQ8D3gMVmVg0sDrezVi/XeA1QDtwXfu9c3+eHZDjNiCEiIllDNS0REckaSloiIpI1lLRERCRrKGmJiEjWUNISEZGsoaQlIiJZQ0lLJM0scIuZ7TWzp9Idj0gmU9ISSb/TCG5snejuJx7MB5nZJWb2r8EJSyTzKGmJpN8RwEZ3b0h3IOFSJCIZS0lLBDCzjWb2/8LVXRvMbImZHWZmfzOzOjO7v2NVWzO708y2mVmNmT1qZrPD8oJwmpxPhdu5ZvZvM/taH+f9EMHaY6eYWb2ZXRmWvzX8rH1m9h8zmxvzni+a2YYwrlVm9s6w/Cjg+pjP2heWP2xmH455f5famJm5mX0inMqoOiybGa7ku8fM1prZBTHHvzk8b52ZbTazzx/0L0AkUe6uhx7D/gFsJFgg7zCCxTl3EMy9dyxQCDwIfD089lKCudwKCdYPez7mc+YAe4GjgP8JPzO3n3NfAvwrZvu48PwnEawafHEYX2G4/3yCJUNygHcTLI8yrqfPCsseBj7cx/mcYHXbkUAxUAq8BnyQzqVYdgGzw+O3AqeHrw8Bjkv370+P4fNQTUuk08/cfbsHiyA+Bjzp7s+5ewtwN0ECw91vdve6sPx/gXlmVhnuWwF8Kzz+88AH3D0ywDg+Atzg7k+6e8TdlwItwMnhOe509y3uHnX33xHUjg6qLwz4rrvvcfcmgiUsNrr7Le7e7sGK278HzguPbQNmmVmFu+8N94ukhJKWSKftMa+betguC5v8vhc2z9US1IAARsUcuxSYDPzV3asPII4jgM+FTYP7wma+SQS1K8zsopimw30EtbtRfXxeIl6LO/9Jced/HzA23P8u4M3AK2b2iJmdcpDnFkmYOl1FBua9BKv6LiJIWJUEzYEWc8x1wJ+Bs8zsNHcf6Gi+14Bvu/u343eY2RHALwjWuHrc3SNm9nzM+XtatqEBKInZHtvDMbHvew14xN0X9xScuz8NnBMuavpJ4A6CpCqSdKppiQxMOUFT3W6CRNBlLS0z+wAwn6Df6HJgqZmVDfAcvwAuM7OTwnu4Ss3sLWZWTtDf5MDO8HwfJKhpddgOTDSzgpiy54Fzw0UdpwEf6uf8fwamm9kHzCw/fJxgZkeFg03eZ2aVHqxqXAsMtPlT5IApaYkMzG3AK8BmghWLn+jYYWaHEwzMuMjd693918Ay4McDOYG7LyPo17qGoBa3niAJ4u6rgB8RLPS3HTga+HfM2x8EVgLbzGxXWPZjoDU8filwez/nrwPeSLD69BZgG/B9goEnECxmujFsHr0MGFKrU0tm0yKQIiKSNVTTEhGRrKGkJZIC4U3K9T08vpzu2ESyiZoHRUQka6imJSIiWUNJS0REsoaSloiIZA0lLRERyRpKWiIikjX+Pwlhrl/gUgWrAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot max_features (x-axis) versus RMSE (y-axis)\n", "plt.plot(feature_range, RMSE_scores)\n", "plt.xlabel('max_features')\n", "plt.ylabel('RMSE (lower is better)')" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(290.0078511328435, 10)" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# show the best RMSE and the corresponding max_features\n", "sorted(zip(RMSE_scores, feature_range))[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fitting a Random Forest with the best parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# max_features=8 is best and n_estimators=150 is sufficiently large\n", "rfreg = RandomForestRegressor(n_estimators=150, max_features=8, max_depth=3, oob_score=True, random_state=1)\n", "rfreg.fit(X, y)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
featureimportance
7League0.003603
12NewLeague0.004290
8Division0.005477
10Assists0.023842
11Errors0.028618
2HmRun0.044607
9PutOuts0.060063
3Runs0.071800
0AtBat0.094592
4RBI0.130965
5Walks0.139899
1Hits0.145264
6Years0.246980
\n", "
" ], "text/plain": [ " feature importance\n", "7 League 0.003603\n", "12 NewLeague 0.004290\n", "8 Division 0.005477\n", "10 Assists 0.023842\n", "11 Errors 0.028618\n", "2 HmRun 0.044607\n", "9 PutOuts 0.060063\n", "3 Runs 0.071800\n", "0 AtBat 0.094592\n", "4 RBI 0.130965\n", "5 Walks 0.139899\n", "1 Hits 0.145264\n", "6 Years 0.246980" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute feature importances\n", "pd.DataFrame({'feature':feature_cols, 'importance':rfreg.feature_importances_}).sort_values('importance')" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5274187002769267" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# compute the out-of-bag R-squared score\n", "rfreg.oob_score_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reducing X to its most important features\n" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(263, 13)" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check the shape of X\n", "X.shape" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,\n", " max_features=8, max_leaf_nodes=None, min_impurity_decrease=0.0,\n", " min_impurity_split=None, min_samples_leaf=1,\n", " min_samples_split=2, min_weight_fraction_leaf=0.0,\n", " n_estimators=150, n_jobs=1, oob_score=True, random_state=1,\n", " verbose=0, warm_start=False)" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rfreg" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(263, 4)\n", "(263, 5)\n", "(263, 7)\n" ] } ], "source": [ "# set a threshold for which features to include\n", "from sklearn.feature_selection import SelectFromModel\n", "print(SelectFromModel(rfreg, threshold=0.1, prefit=True).transform(X).shape)\n", "print(SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X).shape)\n", "print(SelectFromModel(rfreg, threshold='median', prefit=True).transform(X).shape)" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# create a new feature matrix that only includes important features\n", "X_important = SelectFromModel(rfreg, threshold='mean', prefit=True).transform(X)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "284.35550515135395" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check the RMSE for a Random Forest that only includes important features\n", "rfreg = RandomForestRegressor(n_estimators=150, max_features=3, random_state=1)\n", "scores = cross_val_score(rfreg, X_important, y, cv=10, scoring='neg_mean_squared_error')\n", "np.mean(np.sqrt(-scores))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparing Random Forests with decision trees\n", "\n", "**Advantages of Random Forests:**\n", "\n", "- Performance is competitive with the best supervised learning methods\n", "- Provides a more reliable estimate of feature importance\n", "- Allows you to estimate out-of-sample error without using train/test split or cross-validation\n", "\n", "**Disadvantages of Random Forests:**\n", "\n", "- Less interpretable\n", "- Slower to train\n", "- Slower to predict" ] } ], "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }