GBM in Azure ML Jupyter Notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GBM in Azure ML Jupyter Notebook\n",
    "Lixun Zhang  \n",
    "Nov 10, 2015\n",
    "\n",
    "## 1 Introduction\n",
    "Before starting this example, it is assumed that you are familiar with Azure Machine Learning (Azure ML) experiments. If not, you can refer to [Data Scientists' Guide][guide link]. \n",
    "\n",
    "The purpose of this example is to demonstrate how to use Azure ML Jupyter notebook to fit a Gradient Boosting Machine (GBM) model.\n",
    "\n",
    "GBM is well-known among data scientists and as a [Kaggle Profile][kaggle link] explains, it has several major advantages compared with traditional statistical models like linear regression:\n",
    "    1. it automatically approximates non-linear transformations and interactions\n",
    "    2. it treats missing values without having to fill in values or remove observations\n",
    "    3. monotonic transformation of features won't influence the model's performance\n",
    "\n",
    "In an Azure ML experiment, we can use the \"Boosted Decision Tree Regression\" module to fit a GBM model when the response variable is continuous. This module, however, does not allow users to specify the types of loss functions (for statisticians, this means that you can't specify the distribution for the response variable). One way to address this challenge is to make use of the [sklearn][sklearn link] package for Python. In what follows we'll illustrate how to do this within an Azure ML Jupyter notebook. \n",
    "\n",
    "[guide link]: https://gallery.cortanaanalytics.com/Experiment/Tutorial-for-Data-Scientists-3\n",
    "[kaggle link]: http://blog.kaggle.com/2015/06/22/profiling-top-kagglers-owen-zhang-currently-1-in-the-world/\n",
    "[sklearn link]: http://scikit-learn.org/stable/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2 Read Data\n",
    "Following the same procedure as [Tutorial for Using Azure Machine Learning Notebook][tutorial link], we'll read data from an experiment with the following code.\n",
    "\n",
    "[tutorial link]: link"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from azureml import Workspace\n",
    "ws = Workspace(\n",
    " workspace_id='b2bbeb56a1d04e1599d2510a06c59d87',\n",
    " authorization_token=<token>,\n",
    " endpoint='https://studioapi.azureml.net'\n",
    ")\n",
    "experiment = ws.experiments['b2bbeb56a1d04e1599d2510a06c59d87.\\\n",
    "f-id.911630d13cbe4407b9fe408b5bb6ddef']\n",
    "ds = experiment.get_intermediate_dataset(\n",
    " node_id='a0a931cf-9fb3-4cb9-83db-f48211be560c-323',\n",
    " port_name='Results dataset',\n",
    " data_type_id='GenericCSV'\n",
    ")\n",
    "frame = ds.to_dataframe()\n",
    "\n",
    "mydata = frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Estimate Hyperparameters\n",
    "In a GBM model, there are several hyperparameters and we need to estimate them first. One way to estimate these parameters is to use cross validation on a parameter-grid. In our example, we'll optimize the following parameters over a grid: number of estimators, maximum tree depth, minimum number of samples on a split, and learning rate. To do this we start by providing several values for each of them and create a set of combinations, each combination consisting of one value for each parameter. Then for each combination we use cross validation to estimate the performance, using mean squared error as performance metric.\n",
    "\n",
    "[gbm link]: https://gradientboostedmodels.googlecode.com/git/gbm/inst/doc/gbm.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Tuning hyper-parameters for mean_squared_error\n",
      "The grid cross validation lasted 865.39 seconds\n",
      "Best parameters set found on development set:\n",
      "{'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "Grid scores on development set:\n",
      "-58.797 (+/-51.215) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-41.242 (+/-40.364) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-18.555 (+/-25.300) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-58.797 (+/-51.215) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-41.242 (+/-40.364) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-18.555 (+/-25.300) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-58.797 (+/-51.215) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-41.242 (+/-40.364) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-18.509 (+/-25.283) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 2}\n",
      "-53.845 (+/-42.896) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-38.126 (+/-31.772) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-21.410 (+/-20.525) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-53.845 (+/-42.896) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-38.126 (+/-31.772) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-21.410 (+/-20.525) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-53.910 (+/-42.933) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-37.759 (+/-32.024) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-21.224 (+/-20.414) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 4}\n",
      "-52.514 (+/-36.711) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-37.161 (+/-28.747) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-29.040 (+/-30.215) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-52.514 (+/-36.711) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-37.161 (+/-28.747) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-29.040 (+/-30.215) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-51.178 (+/-36.634) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-35.353 (+/-28.416) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-27.474 (+/-30.702) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.001, 'max_depth': 8}\n",
      "-20.181 (+/-27.408) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.655 (+/-24.999) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.857 (+/-23.010) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-20.181 (+/-27.408) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.655 (+/-24.999) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.857 (+/-23.010) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-20.181 (+/-27.408) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.619 (+/-25.211) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-18.613 (+/-22.535) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 2}\n",
      "-22.122 (+/-20.959) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-21.494 (+/-20.573) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-22.283 (+/-22.320) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-22.122 (+/-20.959) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-21.494 (+/-20.573) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-22.283 (+/-22.320) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-21.559 (+/-20.461) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-20.998 (+/-20.496) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-21.632 (+/-20.797) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 4}\n",
      "-29.367 (+/-30.096) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-29.215 (+/-30.054) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-29.207 (+/-30.058) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-29.367 (+/-30.096) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-29.215 (+/-30.054) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-29.207 (+/-30.058) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-27.097 (+/-30.452) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-26.983 (+/-30.476) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-26.980 (+/-30.467) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.01, 'max_depth': 8}\n",
      "-18.197 (+/-23.375) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-18.811 (+/-24.261) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-19.255 (+/-25.165) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-18.197 (+/-23.375) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-18.811 (+/-24.261) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-19.255 (+/-25.165) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-18.177 (+/-24.002) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-18.581 (+/-24.695) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-19.166 (+/-25.754) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 2}\n",
      "-21.753 (+/-21.421) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.760 (+/-21.475) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.761 (+/-21.476) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.753 (+/-21.421) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.760 (+/-21.475) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.761 (+/-21.476) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.553 (+/-19.929) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.567 (+/-19.970) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-21.567 (+/-19.972) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 4}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 1, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 1, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 1, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 2, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 2, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-29.011 (+/-29.515) for {'min_samples_split': 2, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-26.989 (+/-28.936) for {'min_samples_split': 4, 'n_estimators': 500, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-26.989 (+/-28.936) for {'min_samples_split': 4, 'n_estimators': 1000, 'learning_rate': 0.1, 'max_depth': 8}\n",
      "-26.989 (+/-28.936) for {'min_samples_split': 4, 'n_estimators': 10000, 'learning_rate': 0.1, 'max_depth': 8}\n"
     ]
    }
   ],
   "source": [
    "# create X and y\n",
    "feature_cols = ['crim', 'zn', 'indus', 'chas', 'nox', 'rm', 'age', \n",
    "                'dis', 'rad', 'tax', 'ptratio', 'black', 'lstat']\n",
    "X = mydata[feature_cols]\n",
    "y = mydata.medv\n",
    "\n",
    "#%% use cross validation to estimate hyperparameters on a grid\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from time import time\n",
    "\n",
    "tuned_params = {'n_estimators': [500, 1000, 10000], \n",
    "                'max_depth': [2, 4, 8], \n",
    "                'min_samples_split': [1, 2, 4],\n",
    "                'learning_rate': [0.001, 0.01, 0.1]}\n",
    "          \n",
    "print(\"# Tuning hyper-parameters for mean_squared_error\")\n",
    "\n",
    "gscv = GridSearchCV(GradientBoostingRegressor(loss = 'ls', random_state=0), \n",
    "                    tuned_params, cv=5, scoring='mean_squared_error')\n",
    "\n",
    "start = time()\n",
    "gscv.fit(X, y)\n",
    "print(\"The grid cross validation lasted {0:0.2f} seconds\"\n",
    "      .format(time() - start))\n",
    "print(\"Best parameters set found on development set:\")\n",
    "print(gscv.best_params_)\n",
    "print(\"Grid scores on development set:\")\n",
    "for params, mean_score, scores in gscv.grid_scores_:\n",
    "    print(\"{0:0.3f} (+/-{1:0.3f}) for {2:s}\"\n",
    "          .format(mean_score, scores.std() * 2, params))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Fit Model with Estimated Parameters\n",
    "With the selected parameter values from above, we can fit a GBM model and check model performance. A plot comparing the importance of variables is also generated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean Absolute Error: 0.844351\n",
      "Root Mean Squared Error: 1.071576\n",
      "Relative Absolute Error: 0.127023\n",
      "Relative Squared Error: 0.013602\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAAGrCAYAAAAxVgQJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8Zfd8//HXOzOISGSEiiCc8JO6/FwSl7o7CD91Ca17\nXRJNi1IxLX4u1Z9JtYpqO0RVEcm4lYQIQTURDql7iCQkEcJBkAnhSIJcTD6/P9Y6yZ6Tc2b2TGbv\nfb4zr+fjcR6zbnvvz9qzztrv893f9V2pKiRJkqRW7TTpAiRJkqRrw0ArSZKkphloJUmS1DQDrSRJ\nkppmoJUkSVLTDLSSJElqmoFWkpaJJC9P8vZJ1yFJrYnj0EraHiSZBW4KbOgXFbBvVZ1/LZ/zT6vq\n09e6wMYkWQPctqqeMelaJGlzVk66AEnaRgp49DYOnwVkax+cZEVVbdj8lstLEj8bJDXFLgeStmtJ\ndk9yRJKfJDkvyauT7NSvu22STyf5eZKfJXlPkt37de8GbgUcn+TiJC9OMp3kRwuefzbJQ/rpNUk+\nmOTdSX4FHLSp11+k1jX965JkKsmVSQ5O8sMkFyZ5bpJ7Jjk9yS+THD7w2IOTfD7J4Unmkpw1X1e/\n/uZJPto/z3eS/NmC1x2s+znAy4En9/t+ar/ds5KcmeSiJOcmefbAc0z3+/fXSdb3+3vwwPrrJ/nn\n/v2aS3Jykp37dfdO8oV+n76R5EFb9Z8taYdloJW0PVmsNfUo4HLgtsB+wMOBPxtY/w/AXsAdgL2B\nNQD9V+0/pGv13a2q3rDEay7st3UgcExV7Q68b4jX39RzAdwL+F/AU4A3Aq8AHgLcCXhSkgcu2Pa7\nwI2BVwHHJlnVr3t/vz97AU8AXpPkwUvUfQTwGuD9/b7v12+zHnhUVd0QeBbwr0n2G3iOPYEbAjcH\nDgH+bf4PBOAN/f7fB9gDeAlwZZJbAB8D/q6qbgS8GPhQkpss8R5J0jUYaCVtLwIc17fy/TLJsUn2\nBP4Q+Kuq+m1V/QxYSxcOqapzq+qkqrqiqn4O/CtwbVsHv1BVH+2nd9/U6y+xDwu9uqour6oTgYuB\n91XVz6vqJ8DJdCFx3gVV9caq2lBVRwPfBh6dZG/gvsBL++c6DXgH8MzF6q6qS/taNqqnqj5RVd/v\npz8HnAA8YGCTK+iC6Yaq+i/gEuD3+xbpZwEvrKqfVtWVVfWlqroceDrwiar6ZP+8nwJOAR65xHsk\nSddgPylJ24sCHjvYhzbJvYDrAD9NrspmO9G1VNIH3jcC9wd269f94lrWcd7A9K039fpDWj8w/dtF\n5m8wMP/jBY/9AV2L7F7AL6rq1wPrfgjcY4m6F5XkD+lafm9Htx+7AKcPbHJhVV05MP8bYFfgJsDO\nwLmLPO2tgScmeczAspXADnchnqStZ6CVtD37EXAZcOMFQWvea+hGRfjfVTWX5HHA4QPrF3YB+DVd\niAO6i76A31uwzeBjNvf6C13bYWdusWD+1sBHgJ8AeyTZtaou6dfdio1D7MLX3qjeJNcDPkTXovqR\nqtqQ5MMMd9Hcz4FL6bpOnL5g3Q+Bd1fVs6/xKEkakl0OJG23quqndF+L/0uS3ZLs1F8INt/vdFe6\nkHpR35fzJQueYj1d39d55wA7J3lkkusArwSudy1ef6GtGVFh8DE3TXJokuskeSJwe7qv888DvgD8\nY5LrJbkL8KfAezbxvOuBqVzdtHzd/ufndH1f/5CuP/Bm9WH+nXTvw15JViS5T5Lr9jU8JsnD++U7\n9xeYLQznkrQkA62k7d0z6YLYmXTdCY4BbtavOwzYH/gVcDxdC+RgS+U/Aq/s++T+dVX9CngeXf/T\n8+j6iA6OelBcs6VzU6+/0MLHD9NiO7jNl+m6A/wMeDXw+Kr6Zb/uqcAUXWvtscD/G+iesVjdx/T/\nXpjklKq6GDgUOLrfj6fStf4uVctCLwbOAL4KXEj33u7Uh+3H0l3sdgFdi+2L8PNJ0hYY2Y0VkrwT\neBTdRQp37pftAXyA7muwWeBJVTXXr3s5XYvBBuDQqjphJIVJ0naoHyLrkKp6wOa2laTtzSj/Aj4S\neMSCZS8DTqyqfYGT+nmS3BF4MnDH/jFvWWqcRkmSJGnQyEJjVZ0M/HLB4gOBdf30OuBx/fRjgf/s\nh86ZpRtH8V6jqk2StkOLdRuQpB3CuFtB96yq+SFn1tMNwg3dINyDV9uexzWv1pUkLaGq1lXVUheb\nSdJ2bWLDdlVVJdlUa8I11m1me0mSJG3HqmrR0WDGHWjXJ7lZVZ2fZC+6K1qhGwx874Htbsk1BwgH\n4FWvetVICpuZmWF6enqLHzc7O8tRRx21zesZ1po1a1izZs3EXn9rWfd4Wff4tFgzWPe4Wfd4Wff2\nYeAGNdcw7i4HHwUO6qcPAo4bWP6UJNdNsg/dsDNfGXNtkiRJatDIWmiT/CfdPdFvkuRHwP8DXgsc\nneQQ+mG7AKrqzCRH043T+DvgeTWq8cQkSZK0XRlZoK2qpy6x6oAltn8N3W0oJ2JqampSL32tbE03\nieXAusfLusenxZrBusfNusfLurd/I7uxwigkqVH1od1ak+5DK0mStCNIsuRFYd68QJIkSU0z0EqS\nJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJTTPQSpIkqWkG\nWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmrZy0gXsCFavXs3c3Nyky9jIqlWrWLt27aTLkCRJutYM\ntGMwNzfH1NTUpMvYyOzs7KRLkCRJ2ibsciBJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqa\ngVaSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS0wy0kiRJapqBVpIk\nSU0z0EqSJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTVk66AC1fq1evZm5ubtJl\nbGTVqlWsXbt20mVIkqRlxECrJc3NzTE1NTXpMjYyOzs76RIkSdIyY5cDSZIkNc1AK0mSpKYZaCVJ\nktQ0A60kSZKaZqCVJElS0wy0kiRJapqBVpIkSU0z0EqSJKlpBlpJkiQ1zUArSZKkphloJUmS1DQD\nrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJTTPQSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmS\nmmaglSRJUtMMtJIkSWqagVaSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKZNJNAmeXmSbyU5I8n7klwv\nyR5JTkxyTpITkqyaRG2SJElqy9gDbZIp4M+B/avqzsAK4CnAy4ATq2pf4KR+XpIkSdqkSbTQXgRc\nAeySZCWwC/AT4EBgXb/NOuBxE6hNkiRJjRl7oK2qXwD/DPyQLsjOVdWJwJ5Vtb7fbD2w57hrkyRJ\nUntWjvsFk9wWWA1MAb8Cjkny9MFtqqqS1GKPn5mZuWp6amqKqampUZUqSZKkCZmZmdko923K2AMt\ncA/gC1V1IUCSY4H7AOcnuVlVnZ9kL+CCxR48PT09tkIlSZI0GdPT0xvlvsMOO2zJbSfRh/Zs4N5J\nrp8kwAHAmcDxwEH9NgcBx02gNkmSJDVm7C20VXVakncBpwBXAl8H3gbsBhyd5BBgFnjSuGuTJElS\neybR5YCqej3w+gWLf0HXWitJkiQNzTuFSZIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmag\nlSRJUtMMtJIkSWqagVaSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS\n0wy0kiRJapqBVpIkSU0z0EqSJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSS\nJElqmoFWkiRJTTPQSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqa\ngVaSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS0wy0kiRJapqBVpIk\nSU0z0EqSJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJTTPQ\nSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqagVaSJElNM9BKkiSp\naQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS0wy0kiRJapqBVpIkSU0z0EqSJKlpEwm0\nSVYl+WCSs5KcmeQPkuyR5MQk5yQ5IcmqSdQmSZKktkyqhfaNwCeq6g7AXYCzgZcBJ1bVvsBJ/bwk\nSZK0SWMPtEl2Bx5QVe8EqKrfVdWvgAOBdf1m64DHjbs2SZIktWcSLbT7AD9LcmSSryd5e5IbAHtW\n1fp+m/XAnhOoTZIkSY1ZOaHX3B/4y6r6apK1LOheUFWVpBZ78MzMzFXTU1NTTE1Nja5SSZIkTcTM\nzMxGuW9TJhFozwPOq6qv9vMfBF4OnJ/kZlV1fpK9gAsWe/D09PR4qpQkSdLETE9Pb5T7DjvssCW3\nHXuXg6o6H/hRkn37RQcA3wKOBw7qlx0EHDfu2iRJktSeSbTQArwAeG+S6wLnAs8CVgBHJzkEmAWe\nNKHaJEmS1JCJBNqqOg245yKrDhh3LZIkSWqbdwqTJElS0wy0kiRJapqBVpIkSU0z0EqSJKlpBlpJ\nkiQ1zUArSZKkpg0VaJNMJTmgn94lyQ1HW5YkSZI0nM0G2iTPBo4B/qNfdEvgw6MsSpIkSRrWMC20\nzwfuD1wEUFXnADcdZVGSJEnSsIYJtJdV1WXzM0lWAjW6kiRJkqThDRNoP5vkb4BdkjyMrvvB8aMt\nS5IkSRrOMIH2ZcDPgDOA5wCfAF45yqIkSZKkYa0cYpudgSOq6m0ASVYA1wd+M8rCJEmSpGEM00L7\naboAO28X4FOjKUeSJEnaMsME2utV1SXzM1V1MV2olSRJkiZumED76yR3n59Jcg/gt6MrSZIkSRre\nMH1oVwNHJ/lpP78X8OTRlSRJkiQNb7OBtqq+muQOwO/TjT/77aq6YuSVSZIkSUMYpoUW4B7APv32\n+yehqt41urIkSZKk4Ww20CZ5D3Ab4BvAhoFVBlpJkiRN3DAttHcH7lhV3u5WkiRJy84woxx8k+5C\nMEmSJGnZGaaF9veAM5N8BbisX1ZVdeDoypIkSZKGM0ygXTPqIiRJkqStNcywXTNjqEOSJEnaKpvt\nQ5vkPkm+muSSJFckuTLJReMoTpIkSdqcYS4KezPwJ8B3gJ2BQ4C3jLIoSZIkaVjDBFqq6jvAiqra\nUFVHAo8YbVmSJEnScIa5KOzXSa4HnJbk9cD5QEZbliRJkjScYVpon9Fv95fAb4BbAo8fZVGSJEnS\nsIZpoX1cVb0R+C39EF5JXgi8cYR1SVtt9erVzM3NTbqMjaxatYq1a9dOugxJkrZLwwTag7lmeH3W\nIsukZWFubo6pqalJl7GR2dnZSZcgSdJ2a8lAm+SpdKMb7JPk+IFVuwEXjrowSZIkaRibaqH9AvBT\n4CbAG7j6QrCLgNNHXJckSZI0lCUDbVX9IMmPgcuq6rNjrEmSJEka2iZHOaiq3wEbkqwaUz2SJEnS\nFhlqHFrgjCQn9tMAVVWHjq4sSZIkaTjDBNpj+5/q5zMwLUmSJE3UZgNtVR3V3yls337R2VV1xWjL\nkiRJkoaz2UCbZBpYB/ygX3SrJAd5oZgkSZKWg2G6HPwL8PCq+jZAkn2B9wP7j7IwaUfjHc4kSdo6\nwwTalfNhFqCqzkkyzOMkbQHvcCZJ0tYZJph+Lck7gPfQXRD2NOCUkVYlSZIkDWmYQPsXwPOB+WG6\nTgbeMrKKJEmSpC0wzCgHlyZ5M/ApuuG6zq6qy0demSRJkjSEYUY5eBTwVuB7/aLbJHlOVX1ipJVJ\nasJyu5jNC9kkaccz7CgHD66q7wIkuS3wif5H0g5uuV3M5oVskrTj2WmIbS6aD7O97wEXjageSZIk\naYsMO8rBJ4Cj+/knAqck+WOAqjp2VMVJkiRJmzNMoN0ZuAB4UD//s37ZY/p5A60kSZImZphRDg4e\nQx2SJEnSVhlmlIPbAC8Apga2r6o6cIR1SZIkSUMZpsvBccA7gOOBK/tlNbKKJEmSpC0wTKC9tKre\nNPJKJEmSpK0wTKA9PMka4L+By+YXVtXXR1WUJEmSNKxhAu2dgGcAD+bqLgf085IkSdJEDRNonwjs\nU1WXj7oYSZIkaUsNc6ewM4AbjboQSZIkaWsM00J7I+DsJF/l6j60DtslSZKkZWGYQPuqkVchSZIk\nbaVh7hQ2M4Y6JEmSpK2yZB/aJJ/v/70kycULfi4aX4mSJEnS0pZsoa2q+/X/7jq+ciRJkqQtM8wo\nB5IkSdKyZaCVJElS0yYWaJOsSHJqkuP7+T2SnJjknCQnJFk1qdokSZLUjkm20L4QOBOofv5lwIlV\ntS9wUj8vSZIkbdJEAm2SWwKPBN4BpF98ILCun14HPG4CpUmSJKkxk2qh/VfgJcCVA8v2rKr1/fR6\nYM+xVyVJkqTmDHOnsG0qyaOBC6rq1CTTi21TVZWkFls3MzNz1fTU1BRTU1MjqFKSJEmTNDMzs1Hu\n25SxB1rgvsCBSR4J7AzcMMm7gfVJblZV5yfZC7hgsQdPT0+Pr1JJkiRNxPT09Ea577DDDlty27F3\nOaiqV1TV3lW1D/AU4NNV9Qzgo8BB/WYHAceNuzZJkiS1ZzmMQzvfteC1wMOSnAM8pJ+XJEmSNmkS\nXQ6uUlWfBT7bT/8COGCS9UiSJKk9y6GFVpIkSdpqBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppm\noJUkSVLTDLSSJElqmoFWkiRJTTPQSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJ\nUtMMtJIkSWraykkXIEmTsHr1aubm5iZdxlVWrVrF2rVrJ12GJDXJQCtphzQ3N8fU1NSky7jK7Ozs\npEuQpGbZ5UCSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKa5rBdktSI5TZ2Ljh+\nrqTlwUArSY1YbmPnguPnSloe7HIgSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJ\nTTPQSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqagVaSJElNM9BK\nkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS0wy0kiRJapqBVpIkSU0z0EqSJKlp\nBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJTTPQSpIkqWkGWkmS\nJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqagVaSJElNM9BKkiSpaQZaSZIkNW3s\ngTbJ3kk+k+RbSb6Z5NB++R5JTkxyTpITkqwad22SJElqzyRaaK8A/qqq7gTcG3h+kjsALwNOrKp9\ngZP6eUmSJGmTxh5oq+r8qvpGP30JcBZwC+BAYF2/2TrgceOuTZIkSe2ZaB/aJFPAfsCXgT2ran2/\naj2w54TKkiRJUkNWTuqFk+wKfAh4YVVdnOSqdVVVSWqxx83MzFw1PTU1xdTU1GgLlSRJ0tjNzMxs\nlPs2ZSKBNsl16MLsu6vquH7x+iQ3q6rzk+wFXLDYY6enp8dUpSRpW1i9ejVzc3OTLmMjq1atYu3a\ntZMuQ9ImTE9Pb5T7DjvssCW3HXugTdcUewRwZlUNnk0+ChwEvK7/97hFHi5Jaszc3Nyy+zZtdnZ2\n0iVI2oYm0UJ7P+DpwOlJTu2XvRx4LXB0kkOAWeBJE6hNkiRJjRl7oK2q/2Hpi9EOGGctkiRJap93\nCpMkSVLTDLSSJElqmoFWkiRJTZvYOLSSJC1nDjcmtcNAK0nSIhxuTGqHXQ4kSZLUNAOtJEmSmmag\nlSRJUtMMtJIkSWqagVaSJElNM9BKkiSpaQZaSZIkNc1AK0mSpKYZaCVJktQ0A60kSZKaZqCVJElS\n0wy0kiRJapqBVpIkSU0z0EqSJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSS\nJElqmoFWkiRJTTPQSpIkqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtNWTroASZK0\n7axevZq5ublJl7GRVatWsXbt2kmXoe2YgVaSpO3I3NwcU1NTky5jI7Ozs5MuQds5uxxIkiSpaQZa\nSZIkNc1AK0mSpKbZh1aSJE2cF7Pp2jDQSpKkifNiNl0bdjmQJElS0wy0kiRJapqBVpIkSU0z0EqS\nJKlpBlpJkiQ1zUArSZKkphloJUmS1DQDrSRJkppmoJUkSVLTDLSSJElqmoFWkiRJTTPQSpIkqWkG\nWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtMMtJIkSWqagVaSJElNM9BKkiSpaSsnXYAk\nSVKrVq9ezdzc3KTLuMqqVatYu3btpMsYOwOtJEnSVpqbm2NqamrSZVxldnZ20iVMhF0OJEmS1DQD\nrSRJkppmoJUkSVLTDLSSJElq2rIKtEkekeTsJN9J8tJxvnarnaite7yse7xarLvFmsG6x826x8u6\nt3/LJtAmWQG8GXgEcEfgqUnuMK7Xb/Wgse7xsu7xarHuFmsG6x436x4v697+LZtAC9wL+G5VzVbV\nFcD7gcdOuCZJkiQtc8sp0N4C+NHA/Hn9MkmSJGlJqapJ1wBAkscDj6iqP+/nnw78QVW9YGCb5VGs\nJEmSxq6qstjy5XSnsB8Dew/M703XSnuVpXZCkiRJO67l1OXgFOB2SaaSXBd4MvDRCdckSZKkZW7Z\ntNBW1e+S/CXw38AK4IiqOmvCZUmSJGmZWzZ9aCVJkqStsZy6HGhISdYkeVGSw5I8dNL1DEpyyWbW\nv2LI5xlqux1R3y3njEWWzyS5+1Y838FJDt821W0bSR4z7purbIkkuyf5i0nXsT1LsjrJ9bficQcl\n2Wtg/u2jHtM8yee3cPvpJMePqp5x2972Z9ySHNVfGK9rYYcOtOm0+B4UQFW9qqpOmnQxC2yuyf/l\nQz7PsNvpasXm3/+lHrdsJFlRVcdX1esmXcsm3Ah43qSL2M69ENhlsRWbOW8fDNx8fqaq/nzU3deq\n6n6jfP5J6T8jvRh79JbVObhVLYa5a6Vv3fp2knXAxcB3kxzZL3tvkocn+XySc5Lcc9L1zkvyN32N\nJwO/3y87cv6vuiSvTfKtJKcl+aeJFtvVs1eSzyU5NckZSe6f5LXA9ftl7+63Oy7JKUm+mWR+yLZr\nbDfGuj+8SD2H9O/9l/vWnsP75b+X5INJvtL/3HeMpa5M8p4kZyY5ZmFLVpK3JPlqvx9rBpbfsz++\nv5HkS0l2BTKw/lFJvpBkj1EWn+SZ/bH6jSTv6o/ltyb5EvD6vpVt/n0+qt+fLyY5t28NWtfv+5Gj\nrHMTXgvctj9G/yXJp5J8LcnpSQ7s675nv4/XS3KD/v/ijhOqd/7cd1aSt/W1/HeSnZPcrT8WTkty\nbJJVSW7dnwNvnGSnJCcnOWCEdZ294Hh+AV0o/UySk/rtLknyhiTfAO6T5G/737szkvxHv80TgHsA\n703y9X7/rvrmIslT+/+jM/rzzLbah0v6f6f71zumf6/fM7DNI/plXwP+aGD5miQvGpj/ZpJb9cfM\nx/vfkTOSPGlb1buZfRn8jDwDOGKJc8mi+7McJHlO/7t5apLvJfl0kouT/H3/fn4xyU0nWN9G579+\n8QP7c/O5ufpzfdclzi0TOTaWvaraoX6AKWAD3Z3Jbg1cAdyJ7kP9FLqL0QAOBD486Xr7Wu4OnA7s\nDOwGfAd4EXAk8MfAjYGzB7a/4QRrvbj/90XAK/rpnYBdB9cPbH+j/t/r0508b7TYdmOsf2E9Nwe+\nD6yiu4jyc8Cb+m3eB9yvn74VcOaYapwCrgTu088f0b/fnwH2X7AfK/rldwauC5wL3L1ft2u//iDg\ncLoPpc8Bu4+4/jsB3wb2mK+1P5Y/ytX9+g8CDu+njwLe108fCFy04Hf2rhM4Tm4NnDHwHu/WT98E\n+M7Adq8G/onutt4vncQxveC4uQK4Sz//AeBpwGnAA/plhwH/2k8fAhwNvAT49wkcz9+fP0b65VcC\nTxiYv9HA9LuAR/fTV/0eDM73v8s/oDtfrgBOAh67jfZh/rw3Dcz1rxXgC8B96c7dPwRuO/Def7Sf\nfhXwooHnOqM/vh4PvG1g+VjO6wx8Rg6+z2x8Lllyf5bTD1efsx/dHz+P6pe/DvibCdW01PnvA/38\nHebPIUudWyZ1bCz3nx2uhbb3g6r6Ct0J5/tV9a3qjopvAZ/qt/km3S/2cvAA4NiqurSqLuaaw5nN\nAZcmOSLJHwG/HXuF1/QV4FlJXgXcuaqW6lv7wr7F5Yt0Yw/fblwFDlnPM4CZqpqrqt8Bx3B1i+YB\nwJuTnAp8BNgtyaJfkY7Aj6rqi/30e4D7L1j/5L7l5Ot0J9A70rXs/7SqvgZQVZdU1Qa6/XkI8H+B\nR1bVr0Zc+0OAo6vqF30dv+yXH9P/Hi5UwHz/vG8C5y/4nZ0acb2LGfwadifgH5OcBpwI3Hyg9efv\ngIfTtRq+frwlLur7VXV6P/014LbAqqo6uV+2DnggQFUdAewOPAd48YjrWng8P6CfHnyfNwAfGph/\nSN+yfDrdMTXY+r3wa/IA96T7Xb6wP+7fS7+v29hXquon/fH5DWAf4PZ07/25/TbvWaTGQUXXiPGw\ndN++3b+qLhpBrUuZ/4yExc8lW7o/k/Im4KSq+hhweVV9vF/+NSb3+b7U+e+4fv4sYM9+2VLnlkke\nG8vWjhpofz0wfdnA9JXA5QPTy2VYs2Lpk0X6k/O9gA/S/SX6yXEVtpT+A/IBdDfMOCrJMxZuk2Qa\neChw76q6G3Aq3V/+E7FEPWez8Xsfru7vFLq72e3X/+xdVb8ZU7mDwW+wJpLsQ9fC9ZCquivwcbr3\ndal+WkUB15GWAAAGG0lEQVTXcrsrfXeWEVvqeN7Uezf4e7nwd3bSv6dPo2s92b+q9gMu4Orj+CbA\nDeje2y2+wGkEBt+7DXTfPAwa7H6yC3BLuv+v3UZc18Lj+cpFll86/wdPkp2BfwMeX1V3Ad7OxueO\npf4wGjSqALbwPV65mdf+HRt/Fu8MUFXfAfaja7H9+yR/u+1LXdKvYYvOJcsuzCY5GNi7qg7rF10x\nsHqS542lzn+XD0zPr1/03DLhY2PZ2lEDbWs+Bzyu7w+2G/CYwZVJbkDXyvJfwF8Dd51AjRtJcivg\nZ1X1DrqvEPfrV12RZP5EckPgl1V1aZLbA/ceeIrB7cZlsXpuADyo71e4ku6rnnknAIfOzyS52xhr\nvVWS+ffrT4D/mS+Dbj9+DVyUZE/gD+lOot8G9kpyj77e3ZKs6B/zA+AJwLsy+n6enwaemL6fbhbv\nr7vsPiAXuJirQ97uwAVVtSHJg+m+Lp73H8Ar6bqnLMeL3H4F/CLJfAv/M4CZfvp1wLvpvhJ/+4jr\nWOx4vpjuWF7MfHi9MF0/8CcOrFvscUX3rdGD0vULXgE8hav3dZSK7g/jqSS36Zc9dWD9LF2XCJLs\nT9eiS7qRGi6tqvcCb5jfZsyWOpdsan8mru8z/SK643m5Geb8N++GLHJuWSbHxrIz6ZaNSaklpje3\nbiKq6tQkH6Dr63YB3Yn5qtV0H6wf6VstAvzV+KvcqB6ABwMvTnIF3QfMM/vlbwNO77/COgR4bpIz\n6cLWFwee56rtqmpcJ6VPLlLPecBr6N7zX9CdyOe/3jkU+Lf+66CVwGcZz5Xv8+H0+UneSfe1+7/T\n/aFTVXVa3w3ibOBH9GG3qq5I8mTg8HQXkf0GeFj/fFVV307yNOCYJI+uqu+PpPiqM5P8A/DZJBvo\nWsIXjtCw2Pxi04vNj1xVXdhfwHEG8FXg9v1X36cAZ9FdIP5M4LKqen+6q/K/kGS6qmbGXe+Axd67\ng4G39i2y59J1FXoQXd/9Q6uqkjw+yUFVtW5EdS12PF8OfDLJj6vqoYO1V9VckrfTd0EBvjzwXEf1\n+/Mbuv6r8485P8nL6PqBBvhYVW2roaY2+blRVZcleTbw8b6uk+n+WIauG8Uzk3yz349v98vvDPxT\nkvlvDsc5TFz1dS91LtnU/iwHz6frm/qZdIM0nMKmzy9jswXnP+i6xRy/4NwCkz02li1vrCBtRpIb\nVNWv+xbaY+kuHPzIpOuStgdJpoDjq+rOEy5FUsPsciBt3pq+leIM4HuGWWmbs2VF0rViC60kSZKa\nZgutJEmSmmaglSRJUtMMtJIkSWqagVaSJElNM9BK0gglOTTJmUnevYWPu3WSZTVgvSQtVwZaSRqt\nvwAO2IobhOxDd9esLdLfxEGSdiie+CRpRJK8FbgN3R2vXpHkiCRfTvL1JAf220wl+VySr/U/9+kf\n/lrgAUlOTbI6yUFJDh947o8leWA/fUmSNyT5BnCfJE/vX+fUJG9NslOSFUmOSnJGktOTrB7z2yFJ\nI2OglaQRqarnAj8BpuluDfrpqvoD4CF0t67cBVgPPKyq7g48BXhT//CXAidX1X5VtXaxpx+Y3gX4\nUlXdje4WzU8C7ltV+wEbgKcBdwVuXlV3rqq7AEdu272VpMlZOekCJGkHEOD/AAcmeXG/7HrA3sD5\nwJuT3JUufN5u4DHD2gB8qJ9+KHB34JT+PvbXpwvNxwO3SfIm4OPACVu9N5K0zBhoJWl8/riqvjO4\nIMka4KdV9YwkK4BLl3js79j4W7WdB6YvrY1v+7iuql6x8AmS3AV4BPBculbcQ7Z8FyRp+bHLgSSN\nx38Dh87PJNmvn7whXSstwDOBFf30xcBuA4+fBe6Wzt7AvZZ4nZOAJyT5vf519khyqyQ3BlZW1bHA\n3wL7X/tdkqTlwRZaSRqt6n9eDaxNcjpdY8L3gAOBtwAfSvJM4JPAJf3jTgM29Bd6HVlVb0zyfeBM\n4Czgawteo5uoOivJK4ET+hEPrgCeR9fye+TAKAgvG8neStIEZONvqSRJkqS22OVAkiRJTTPQSpIk\nqWkGWkmSJDXNQCtJkqSmGWglSZLUNAOtJEmSmmaglSRJUtP+P2va19vYJAo4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x18ab1390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# fit model with the best set of parameter values\n",
    "params = {'n_estimators': 500, 'max_depth': 2, 'min_samples_split': 4,\n",
    "          'learning_rate': 0.1, 'loss': 'ls', 'random_state': 0}\n",
    "          \n",
    "gbm = GradientBoostingRegressor(**params)\n",
    "\n",
    "gbm.fit(X, y)\n",
    "\n",
    "# check additional performance metrics\n",
    "# assign test data\n",
    "import pandas as pd\n",
    "newX = X\n",
    "newY = y\n",
    "\n",
    "# join predictions with original data\n",
    "predicted = gbm.predict(newX)\n",
    "predicted_df = pd.DataFrame({\"predicted\": predicted})\n",
    "mydata_with_pd = newX.join(newY).join(predicted_df)\n",
    "mydata_with_pd.head()\n",
    "\n",
    "# check performance metrics\n",
    "import numpy as np\n",
    "obs = mydata_with_pd.medv\n",
    "pred = mydata_with_pd.predicted\n",
    "\n",
    "mae = np.mean(abs(pred-obs))\n",
    "rmse = np.sqrt(np.mean((pred-obs)**2))\n",
    "rae = np.mean(abs(pred-obs))/np.mean(abs(obs-np.mean(obs)))\n",
    "rse = np.mean((pred-obs)**2)/np.mean((obs-np.mean(obs))**2)\n",
    "\n",
    "# print additional performance metrics\n",
    "print(\"Mean Absolute Error: {0:0.6f}\".format(mae))\n",
    "print(\"Root Mean Squared Error: {0:0.6f}\".format(rmse))\n",
    "print(\"Relative Absolute Error: {0:0.6f}\".format(rae))\n",
    "print(\"Relative Squared Error: {0:0.6f}\".format(rse))\n",
    "\n",
    "# plot variable importance\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "feature_importance = gbm.feature_importances_\n",
    "feature_importance = 100.0 * (feature_importance / feature_importance.max())\n",
    "sorted_idx = np.argsort(feature_importance)[::-1]\n",
    "\n",
    "plt.figure(figsize=(6 * 1.618, 6))\n",
    "index = np.arange(len(feature_cols))\n",
    "bar_width = 0.5\n",
    "plt.bar(index, feature_importance[sorted_idx], color='black', alpha=0.5)\n",
    "plt.xlabel('features')\n",
    "plt.ylabel('importance')\n",
    "plt.title('Feature importance')\n",
    "plt.xticks(index + bar_width, np.array(feature_cols)[sorted_idx])\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Conclusion\n",
    "Using the Boston housing dataset, we started the analysis by estimating the parameters in the GBM model. Then we fitted the model and examined variable importance. \n",
    "\n",
    "In addition to the least squares (ls) loss function, [GradientBoostingRegressor][boostingreg link] from [sklearn][sklearn link] allows for several other loss functions: least absolute deviation (lad), huber loss function which is a combination of ls and lad, and quantile loss function. For example, the least squares loss function is commonly used when the response variable follows a Gaussian distribution and the least absolute deviation loss function is appropriate when the response variable follows a Laplace distribution. \n",
    "\n",
    "\n",
    "[boostingreg link]: http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html\n",
    "[sklearn link]: http://scikit-learn.org/stable/index.html"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Anaconda2]",
   "language": "python",
   "name": "Python [Anaconda2]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}