Air-Quality-Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Air Quality Prediction\n",
    "____\n",
    "\n",
    "In this example we would look at the task of predicting air quality. We would use the following dataset: \n",
    "\n",
    "https://archive.ics.uci.edu/ml/datasets/Air+Quality\n",
    "\n",
    "This dataset contains the responses of a gas multisensor device deployed on the field in an Italian city. Hourly responses averages are recorded along with gas concentrations references from a certified analyzer. The dataset contains 9358 instances of hourly averaged responses from an array of 5 metal oxide chemical sensors embedded in an Air Quality Chemical Multisensor Device. The device was located on the field in a significantly polluted area, at road level, within an Italian city. Data was recorded from March 2004 to February 2005 (one year) representing the longest freely available recordings of on field deployed air quality chemical sensor devices responses. Ground Truth hourly averaged concentrations for CO, Non Metanic Hydrocarbons, Benzene, Total Nitrogen Oxides (NOx) and Nitrogen Dioxide (NO2) are provided by a co-located reference certified analyzer.\n",
    "\n",
    "\n",
    "Attribute Information:\n",
    "\n",
    "0. Date\t(DD/MM/YYYY) \n",
    "1. Time\t(HH.MM.SS) \n",
    "2. True hourly averaged concentration CO in mg/m^3 (reference analyzer) \n",
    "3. PT08.S1 (tin oxide) hourly averaged sensor response (nominally CO targeted)\t\n",
    "4. True hourly averaged overall Non Metanic HydroCarbons concentration in microg/m^3 (reference analyzer) \n",
    "5. True hourly averaged Benzene concentration in microg/m^3 (reference analyzer) \n",
    "6. PT08.S2 (titania) hourly averaged sensor response (nominally NMHC targeted)\t\n",
    "7. True hourly averaged NOx concentration in ppb (reference analyzer) \n",
    "8. PT08.S3 (tungsten oxide) hourly averaged sensor response (nominally NOx targeted) \n",
    "9. True hourly averaged NO2 concentration in microg/m^3 (reference analyzer)\t\n",
    "10. PT08.S4 (tungsten oxide) hourly averaged sensor response (nominally NO2 targeted)\t\n",
    "11. PT08.S5 (indium oxide) hourly averaged sensor response (nominally O3 targeted) \n",
    "12. Temperature in Â°C\t\n",
    "13. Relative Humidity (%) \n",
    "14. AH Absolute Humidity \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the dataset from the link and save it in the same directory as your code. Next we import all the required modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from sklearn.tree import DecisionTreeClassifier, export_graphviz\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from time import time\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.linear_model import SGDRegressor\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.cross_validation import train_test_split, cross_val_score\n",
    "from sklearn import cross_validation, metrics\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.metrics import roc_auc_score , classification_report, mean_squared_error, r2_score\n",
    "from sklearn.metrics import precision_score, recall_score, accuracy_score, classification_report\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.feature_selection import *\n",
    "from sklearn import metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read .csv from provided dataset\n",
    "csv_filename=\"AirQualityUCI.csv\"\n",
    "\n",
    "# df=pd.read_csv(csv_filename,index_col=0)\n",
    "df=pd.read_csv(csv_filename, sep=\";\" , parse_dates= ['Date','Time'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `AirQualityUCI.csv` file as our dataset. It is a ';' seperated file so we'll specify it as a parameter for the read_csv function. We'll also use parse_dates parameter so that pandas recognizes the 'Date' and 'Time' columns and format them accordingly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>Time</th>\n",
       "      <th>CO(GT)</th>\n",
       "      <th>PT08.S1(CO)</th>\n",
       "      <th>NMHC(GT)</th>\n",
       "      <th>C6H6(GT)</th>\n",
       "      <th>PT08.S2(NMHC)</th>\n",
       "      <th>NOx(GT)</th>\n",
       "      <th>PT08.S3(NOx)</th>\n",
       "      <th>NO2(GT)</th>\n",
       "      <th>PT08.S4(NO2)</th>\n",
       "      <th>PT08.S5(O3)</th>\n",
       "      <th>T</th>\n",
       "      <th>RH</th>\n",
       "      <th>AH</th>\n",
       "      <th>Unnamed: 15</th>\n",
       "      <th>Unnamed: 16</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>18.00.00</td>\n",
       "      <td>2,6</td>\n",
       "      <td>1360.0</td>\n",
       "      <td>150.0</td>\n",
       "      <td>11,9</td>\n",
       "      <td>1046.0</td>\n",
       "      <td>166.0</td>\n",
       "      <td>1056.0</td>\n",
       "      <td>113.0</td>\n",
       "      <td>1692.0</td>\n",
       "      <td>1268.0</td>\n",
       "      <td>13,6</td>\n",
       "      <td>48,9</td>\n",
       "      <td>0,7578</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>19.00.00</td>\n",
       "      <td>2</td>\n",
       "      <td>1292.0</td>\n",
       "      <td>112.0</td>\n",
       "      <td>9,4</td>\n",
       "      <td>955.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>1174.0</td>\n",
       "      <td>92.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>972.0</td>\n",
       "      <td>13,3</td>\n",
       "      <td>47,7</td>\n",
       "      <td>0,7255</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>20.00.00</td>\n",
       "      <td>2,2</td>\n",
       "      <td>1402.0</td>\n",
       "      <td>88.0</td>\n",
       "      <td>9,0</td>\n",
       "      <td>939.0</td>\n",
       "      <td>131.0</td>\n",
       "      <td>1140.0</td>\n",
       "      <td>114.0</td>\n",
       "      <td>1555.0</td>\n",
       "      <td>1074.0</td>\n",
       "      <td>11,9</td>\n",
       "      <td>54,0</td>\n",
       "      <td>0,7502</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>21.00.00</td>\n",
       "      <td>2,2</td>\n",
       "      <td>1376.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>9,2</td>\n",
       "      <td>948.0</td>\n",
       "      <td>172.0</td>\n",
       "      <td>1092.0</td>\n",
       "      <td>122.0</td>\n",
       "      <td>1584.0</td>\n",
       "      <td>1203.0</td>\n",
       "      <td>11,0</td>\n",
       "      <td>60,0</td>\n",
       "      <td>0,7867</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>22.00.00</td>\n",
       "      <td>1,6</td>\n",
       "      <td>1272.0</td>\n",
       "      <td>51.0</td>\n",
       "      <td>6,5</td>\n",
       "      <td>836.0</td>\n",
       "      <td>131.0</td>\n",
       "      <td>1205.0</td>\n",
       "      <td>116.0</td>\n",
       "      <td>1490.0</td>\n",
       "      <td>1110.0</td>\n",
       "      <td>11,2</td>\n",
       "      <td>59,6</td>\n",
       "      <td>0,7888</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Date      Time CO(GT)  PT08.S1(CO)  NMHC(GT) C6H6(GT)  PT08.S2(NMHC)  \\\n",
       "0 2004-10-03  18.00.00    2,6       1360.0     150.0     11,9         1046.0   \n",
       "1 2004-10-03  19.00.00      2       1292.0     112.0      9,4          955.0   \n",
       "2 2004-10-03  20.00.00    2,2       1402.0      88.0      9,0          939.0   \n",
       "3 2004-10-03  21.00.00    2,2       1376.0      80.0      9,2          948.0   \n",
       "4 2004-10-03  22.00.00    1,6       1272.0      51.0      6,5          836.0   \n",
       "\n",
       "   NOx(GT)  PT08.S3(NOx)  NO2(GT)  PT08.S4(NO2)  PT08.S5(O3)     T    RH  \\\n",
       "0    166.0        1056.0    113.0        1692.0       1268.0  13,6  48,9   \n",
       "1    103.0        1174.0     92.0        1559.0        972.0  13,3  47,7   \n",
       "2    131.0        1140.0    114.0        1555.0       1074.0  11,9  54,0   \n",
       "3    172.0        1092.0    122.0        1584.0       1203.0  11,0  60,0   \n",
       "4    131.0        1205.0    116.0        1490.0       1110.0  11,2  59,6   \n",
       "\n",
       "       AH  Unnamed: 15  Unnamed: 16  \n",
       "0  0,7578          NaN          NaN  \n",
       "1  0,7255          NaN          NaN  \n",
       "2  0,7502          NaN          NaN  \n",
       "3  0,7867          NaN          NaN  \n",
       "4  0,7888          NaN          NaN  "
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df.dropna(how=\"all\",axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df.dropna(how=\"all\",axis=0,inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data contains null values. So we drop those rows and columns containing nulls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9471, 15)"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last few lines(specifically 9357 to 9471) of the dataset are empty and of no use. So we'll ignore them too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = df[:9357]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>Time</th>\n",
       "      <th>CO(GT)</th>\n",
       "      <th>PT08.S1(CO)</th>\n",
       "      <th>NMHC(GT)</th>\n",
       "      <th>C6H6(GT)</th>\n",
       "      <th>PT08.S2(NMHC)</th>\n",
       "      <th>NOx(GT)</th>\n",
       "      <th>PT08.S3(NOx)</th>\n",
       "      <th>NO2(GT)</th>\n",
       "      <th>PT08.S4(NO2)</th>\n",
       "      <th>PT08.S5(O3)</th>\n",
       "      <th>T</th>\n",
       "      <th>RH</th>\n",
       "      <th>AH</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9352</th>\n",
       "      <td>2005-04-04</td>\n",
       "      <td>10.00.00</td>\n",
       "      <td>3,1</td>\n",
       "      <td>1314.0</td>\n",
       "      <td>-200.0</td>\n",
       "      <td>13,5</td>\n",
       "      <td>1101.0</td>\n",
       "      <td>472.0</td>\n",
       "      <td>539.0</td>\n",
       "      <td>190.0</td>\n",
       "      <td>1374.0</td>\n",
       "      <td>1729.0</td>\n",
       "      <td>21,9</td>\n",
       "      <td>29,3</td>\n",
       "      <td>0,7568</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9353</th>\n",
       "      <td>2005-04-04</td>\n",
       "      <td>11.00.00</td>\n",
       "      <td>2,4</td>\n",
       "      <td>1163.0</td>\n",
       "      <td>-200.0</td>\n",
       "      <td>11,4</td>\n",
       "      <td>1027.0</td>\n",
       "      <td>353.0</td>\n",
       "      <td>604.0</td>\n",
       "      <td>179.0</td>\n",
       "      <td>1264.0</td>\n",
       "      <td>1269.0</td>\n",
       "      <td>24,3</td>\n",
       "      <td>23,7</td>\n",
       "      <td>0,7119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9354</th>\n",
       "      <td>2005-04-04</td>\n",
       "      <td>12.00.00</td>\n",
       "      <td>2,4</td>\n",
       "      <td>1142.0</td>\n",
       "      <td>-200.0</td>\n",
       "      <td>12,4</td>\n",
       "      <td>1063.0</td>\n",
       "      <td>293.0</td>\n",
       "      <td>603.0</td>\n",
       "      <td>175.0</td>\n",
       "      <td>1241.0</td>\n",
       "      <td>1092.0</td>\n",
       "      <td>26,9</td>\n",
       "      <td>18,3</td>\n",
       "      <td>0,6406</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9355</th>\n",
       "      <td>2005-04-04</td>\n",
       "      <td>13.00.00</td>\n",
       "      <td>2,1</td>\n",
       "      <td>1003.0</td>\n",
       "      <td>-200.0</td>\n",
       "      <td>9,5</td>\n",
       "      <td>961.0</td>\n",
       "      <td>235.0</td>\n",
       "      <td>702.0</td>\n",
       "      <td>156.0</td>\n",
       "      <td>1041.0</td>\n",
       "      <td>770.0</td>\n",
       "      <td>28,3</td>\n",
       "      <td>13,5</td>\n",
       "      <td>0,5139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9356</th>\n",
       "      <td>2005-04-04</td>\n",
       "      <td>14.00.00</td>\n",
       "      <td>2,2</td>\n",
       "      <td>1071.0</td>\n",
       "      <td>-200.0</td>\n",
       "      <td>11,9</td>\n",
       "      <td>1047.0</td>\n",
       "      <td>265.0</td>\n",
       "      <td>654.0</td>\n",
       "      <td>168.0</td>\n",
       "      <td>1129.0</td>\n",
       "      <td>816.0</td>\n",
       "      <td>28,5</td>\n",
       "      <td>13,1</td>\n",
       "      <td>0,5028</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           Date      Time CO(GT)  PT08.S1(CO)  NMHC(GT) C6H6(GT)  \\\n",
       "9352 2005-04-04  10.00.00    3,1       1314.0    -200.0     13,5   \n",
       "9353 2005-04-04  11.00.00    2,4       1163.0    -200.0     11,4   \n",
       "9354 2005-04-04  12.00.00    2,4       1142.0    -200.0     12,4   \n",
       "9355 2005-04-04  13.00.00    2,1       1003.0    -200.0      9,5   \n",
       "9356 2005-04-04  14.00.00    2,2       1071.0    -200.0     11,9   \n",
       "\n",
       "      PT08.S2(NMHC)  NOx(GT)  PT08.S3(NOx)  NO2(GT)  PT08.S4(NO2)  \\\n",
       "9352         1101.0    472.0         539.0    190.0        1374.0   \n",
       "9353         1027.0    353.0         604.0    179.0        1264.0   \n",
       "9354         1063.0    293.0         603.0    175.0        1241.0   \n",
       "9355          961.0    235.0         702.0    156.0        1041.0   \n",
       "9356         1047.0    265.0         654.0    168.0        1129.0   \n",
       "\n",
       "      PT08.S5(O3)     T    RH      AH  \n",
       "9352       1729.0  21,9  29,3  0,7568  \n",
       "9353       1269.0  24,3  23,7  0,7119  \n",
       "9354       1092.0  26,9  18,3  0,6406  \n",
       "9355        770.0  28,3  13,5  0,5139  \n",
       "9356        816.0  28,5  13,1  0,5028  "
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cols = list(df.columns[2:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you might have noticed, the values in our data don't contain decimal places but have weird commas in place of them. For example *9.4* is written as *9,4*. We'll correct it using the following piece of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Date</th>\n",
       "      <th>Time</th>\n",
       "      <th>CO(GT)</th>\n",
       "      <th>PT08.S1(CO)</th>\n",
       "      <th>NMHC(GT)</th>\n",
       "      <th>C6H6(GT)</th>\n",
       "      <th>PT08.S2(NMHC)</th>\n",
       "      <th>NOx(GT)</th>\n",
       "      <th>PT08.S3(NOx)</th>\n",
       "      <th>NO2(GT)</th>\n",
       "      <th>PT08.S4(NO2)</th>\n",
       "      <th>PT08.S5(O3)</th>\n",
       "      <th>T</th>\n",
       "      <th>RH</th>\n",
       "      <th>AH</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>18.00.00</td>\n",
       "      <td>2.6</td>\n",
       "      <td>1360.0</td>\n",
       "      <td>150.0</td>\n",
       "      <td>11.9</td>\n",
       "      <td>1046.0</td>\n",
       "      <td>166.0</td>\n",
       "      <td>1056.0</td>\n",
       "      <td>113.0</td>\n",
       "      <td>1692.0</td>\n",
       "      <td>1268.0</td>\n",
       "      <td>13.6</td>\n",
       "      <td>48.9</td>\n",
       "      <td>0.7578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>19.00.00</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1292.0</td>\n",
       "      <td>112.0</td>\n",
       "      <td>9.4</td>\n",
       "      <td>955.0</td>\n",
       "      <td>103.0</td>\n",
       "      <td>1174.0</td>\n",
       "      <td>92.0</td>\n",
       "      <td>1559.0</td>\n",
       "      <td>972.0</td>\n",
       "      <td>13.3</td>\n",
       "      <td>47.7</td>\n",
       "      <td>0.7255</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>20.00.00</td>\n",
       "      <td>2.2</td>\n",
       "      <td>1402.0</td>\n",
       "      <td>88.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>939.0</td>\n",
       "      <td>131.0</td>\n",
       "      <td>1140.0</td>\n",
       "      <td>114.0</td>\n",
       "      <td>1555.0</td>\n",
       "      <td>1074.0</td>\n",
       "      <td>11.9</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0.7502</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>21.00.00</td>\n",
       "      <td>2.2</td>\n",
       "      <td>1376.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>9.2</td>\n",
       "      <td>948.0</td>\n",
       "      <td>172.0</td>\n",
       "      <td>1092.0</td>\n",
       "      <td>122.0</td>\n",
       "      <td>1584.0</td>\n",
       "      <td>1203.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>0.7867</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2004-10-03</td>\n",
       "      <td>22.00.00</td>\n",
       "      <td>1.6</td>\n",
       "      <td>1272.0</td>\n",
       "      <td>51.0</td>\n",
       "      <td>6.5</td>\n",
       "      <td>836.0</td>\n",
       "      <td>131.0</td>\n",
       "      <td>1205.0</td>\n",
       "      <td>116.0</td>\n",
       "      <td>1490.0</td>\n",
       "      <td>1110.0</td>\n",
       "      <td>11.2</td>\n",
       "      <td>59.6</td>\n",
       "      <td>0.7888</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        Date      Time  CO(GT)  PT08.S1(CO)  NMHC(GT)  C6H6(GT)  \\\n",
       "0 2004-10-03  18.00.00     2.6       1360.0     150.0      11.9   \n",
       "1 2004-10-03  19.00.00     2.0       1292.0     112.0       9.4   \n",
       "2 2004-10-03  20.00.00     2.2       1402.0      88.0       9.0   \n",
       "3 2004-10-03  21.00.00     2.2       1376.0      80.0       9.2   \n",
       "4 2004-10-03  22.00.00     1.6       1272.0      51.0       6.5   \n",
       "\n",
       "   PT08.S2(NMHC)  NOx(GT)  PT08.S3(NOx)  NO2(GT)  PT08.S4(NO2)  PT08.S5(O3)  \\\n",
       "0         1046.0    166.0        1056.0    113.0        1692.0       1268.0   \n",
       "1          955.0    103.0        1174.0     92.0        1559.0        972.0   \n",
       "2          939.0    131.0        1140.0    114.0        1555.0       1074.0   \n",
       "3          948.0    172.0        1092.0    122.0        1584.0       1203.0   \n",
       "4          836.0    131.0        1205.0    116.0        1490.0       1110.0   \n",
       "\n",
       "      T    RH      AH  \n",
       "0  13.6  48.9  0.7578  \n",
       "1  13.3  47.7  0.7255  \n",
       "2  11.9  54.0  0.7502  \n",
       "3  11.0  60.0  0.7867  \n",
       "4  11.2  59.6  0.7888  "
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for col in cols:\n",
    "    if df[col].dtype != 'float64':\n",
    "        str_x = pd.Series(df[col]).str.replace(',','.')\n",
    "        float_X = []\n",
    "        for value in str_x.values:\n",
    "            fv = float(value)\n",
    "            float_X.append(fv)\n",
    "\n",
    "            df[col] = pd.DataFrame(float_X)\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "features=list(df.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will define our features and ignore those that might not be of help in our prediction. For example, date is not a very useful feature that can assist in predicting the future values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "features.remove('Date')\n",
    "features.remove('Time')\n",
    "features.remove('PT08.S4(NO2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = df[features]\n",
    "y = df['C6H6(GT)']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will try to predict the C6H6(GT) values. Hence we set it as our target variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the dataset to 60% training and 40% testing sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# split dataset to 60% training and 40% testing\n",
    "X_train, X_test, y_train, y_test = cross_validation.train_test_split(X,y, test_size=0.4, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5614, 12) (5614,)\n"
     ]
    }
   ],
   "source": [
    "print(X_train.shape, y_train.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "# Regression\n",
    "\n",
    "Please see the previous examples for better explanations. We have already implemented Decision Tree Regression and Random Forest Regression to predict the Electrical Energy Output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decision tree regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE train: 5.931, test: 5.361\n",
      "R^2 train: 0.996, test: 0.997\n"
     ]
    }
   ],
   "source": [
    "from sklearn.tree import DecisionTreeRegressor\n",
    "\n",
    "tree = DecisionTreeRegressor(max_depth=3)\n",
    "tree.fit(X_train, y_train)\n",
    "\n",
    "y_train_pred = tree.predict(X_train)\n",
    "y_test_pred = tree.predict(X_test)\n",
    "\n",
    "print('MSE train: %.3f, test: %.3f' % (\n",
    "        mean_squared_error(y_train, y_train_pred),\n",
    "        mean_squared_error(y_test, y_test_pred)))\n",
    "print('R^2 train: %.3f, test: %.3f' % (\n",
    "        r2_score(y_train, y_train_pred),\n",
    "        r2_score(y_test, y_test_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random forest regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE train: 0.006, test: 0.009\n",
      "R^2 train: 1.000, test: 1.000\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "forest = RandomForestRegressor(n_estimators=1000, \n",
    "                               criterion='mse', \n",
    "                               random_state=1, \n",
    "                               n_jobs=-1)\n",
    "forest.fit(X_train, y_train)\n",
    "y_train_pred = forest.predict(X_train)\n",
    "y_test_pred = forest.predict(X_test)\n",
    "\n",
    "print('MSE train: %.3f, test: %.3f' % (\n",
    "        mean_squared_error(y_train, y_train_pred),\n",
    "        mean_squared_error(y_test, y_test_pred)))\n",
    "print('R^2 train: %.3f, test: %.3f' % (\n",
    "        r2_score(y_train, y_train_pred),\n",
    "        r2_score(y_test, y_test_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-squared: 1.0\n"
     ]
    }
   ],
   "source": [
    "regressor = LinearRegression()\n",
    "regressor.fit(X_train, y_train)\n",
    "y_predictions = regressor.predict(X_test)\n",
    "print('R-squared:', regressor.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The R-squared score of 1 indicates that 100 percent of the variance in the test set is explained by the model. The performance can change if a different set of 75 percent of the data is partitioned to the training set. Hence Cross-validation can be used to produce a better estimate of the estimator's performance. Each cross-validation round trains and tests different partitions of the data to reduce variability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average of scores:  1.0\n",
      "Cross validation scores:  [ 1.  1.  1.  1.  1.]\n"
     ]
    }
   ],
   "source": [
    "scores = cross_val_score(regressor, X, y, cv=5)\n",
    "print (\"Average of scores: \", scores.mean())\n",
    "print (\"Cross validation scores: \", scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect some of the model's predictions and plot the true quality scores against the predicted scores:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting models with gradient descent\n",
    "Gradient descent is an optimization algorithm that can be used to estimate the local minimum of a function.\n",
    "\n",
    "We can use gradient descent to find the values of the model's parameters that minimize the value of the cost function. Gradient descent iteratively updates the values of the model's parameters by calculating the partial derivative of the cost function at each step. Although the calculus behind the cost function is not entirely required to implement it with scikit-learn, having an intuition for how gradient descent will always help to you use it effectively.\n",
    "\n",
    "There are two varieties of gradient descent that are distinguished by the number of training instances that are used to update the model parameters in each training iteration.\n",
    "\n",
    "* Batch gradient descent, which is sometimes called only gradient descent, uses all of the training instances to update the model parameters in each iteration.\n",
    "* Stochastic Gradient Descent (SGD), in contrast, updates the parameters using only a single training instance in each iteration. The training instance is usually selected randomly, hence the name stochastic. Stochastic gradient descent is often preferred to optimize cost functions when there are large number of training instances, as it will converge more quickly than batch gradient descent.\n",
    "\n",
    "Batch gradient descent is a deterministic algorithm, and will produce the exact same parameter values given the same training set. As a random algorithm, SGD can produce different parameter estimates each time it is run. SGD may not minimize the cost function as well as gradient descent because it uses only single training instances to update the weights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## SGDRegressor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use Stochastic Gradient Descent to estimate the parameters of a model with\n",
    "Scikit-Learn. SGDRegressor is an implementation of SGD that can be used even for\n",
    "regression problems with large number of features. It can be used\n",
    "to optimize different cost functions to fit different linear models; by default, it will\n",
    "optimize the residual sum of squares. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Scaling the features using StandardScaler:\n",
    "X_scaler = StandardScaler()\n",
    "y_scaler = StandardScaler()\n",
    "X_train = X_scaler.fit_transform(X_train)\n",
    "y_train = y_scaler.fit_transform(y_train)\n",
    "X_test = X_scaler.transform(X_test)\n",
    "y_test = y_scaler.transform(y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross validation r-squared scores: [ 0.99746801  0.9982387   0.99870596  0.99788837  0.99826564]\n",
      "Average cross validation r-squared score: 0.998113337934\n",
      "Test set r-squared score 0.998330577369\n"
     ]
    }
   ],
   "source": [
    "regressor = SGDRegressor(loss='squared_loss')\n",
    "scores = cross_val_score(regressor, X_train, y_train, cv=5)\n",
    "print ('Cross validation r-squared scores:', scores)\n",
    "print ('Average cross validation r-squared score:', np.mean(scores))\n",
    "regressor.fit_transform(X_train, y_train)\n",
    "print ('Test set r-squared score', regressor.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## Selecting the best features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size=0.25, random_state=33)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes there are a lot of features in the dataset, so before learning, we should try to see which features are more relevant for our learning task, i.e. which of them are better prize predictors. We will use the `SelectKBest` method from the `feature_selection` package, and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Date', 'Time', 'CO(GT)', 'PT08.S1(CO)', 'NMHC(GT)', 'C6H6(GT)',\n",
       "       'PT08.S2(NMHC)', 'NOx(GT)', 'PT08.S3(NOx)', 'NO2(GT)', 'PT08.S4(NO2)',\n",
       "       'PT08.S5(O3)', 'T', 'RH', 'AH'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "feature_names = list(df.columns[2:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "feature_names.remove('PT08.S4(NO2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([False,  True, False,  True, False, False, False, False, False,\n",
      "        True,  True,  True], dtype=bool), ['CO(GT)', 'PT08.S1(CO)', 'NMHC(GT)', 'C6H6(GT)', 'PT08.S2(NMHC)', 'NOx(GT)', 'PT08.S3(NOx)', 'NO2(GT)', 'PT08.S5(O3)', 'T', 'RH', 'AH'])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0xf9f29e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAcAAABkCAYAAAD661mlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJztnXm4XEWZh9/vZk8kgQQkIKsGE0AhIQhhEdAgySCLjDNs\nsg46gDgBIqBxCwoOssgimyAQEQQEBR1CTECRRXYua4gkIAlBtkASkpAFQvLNH1Wde+5J973dfbtv\nd9/ze5+nnu5Tp06dr8/5fae6q6u+MndHCCGEEEIIIYQQ2aWp1gYIIYQQQgghhBCitqhzQAghhBBC\nCCGEyDjqHBBCCCGEEEIIITKOOgeEEEIIIYQQQoiMo84BIYQQQgghhBAi46hzQAghhBBCCCGEyDjq\nHBBCCCGEEEIIITKOOgeEEEIIIYQQQoiMo84BIYQQQgghhBAi46hzQAghhBBCCCGEyDjqHOhimNnR\nZra6QPrfKp7338xsYrXqrwRmdpyZzTCz5WY2y8y+VWubRGWR/vNjZiea2a1m9mq8FtfV2iZRHeQD\na2Nmm5jZRDN7zMwWmNk7ZvY3Mxtda9tEZZH+18bMepvZtWb2vJm9Z2ZLzOwZMxtnZt1rbZ+oHNJ/\n+5jZ7vF6rDKzgbW2px7RQ6Fr4sAPgTmp/OlVPOe+wDeBH1fxHGVjZscDVwK3AT8HPg/8wsz6uPv5\nNTVOVBrpf23OAD4GPA4MrrEtovrIB1pzIHA68Efg14TvPkcB95jZse5+fQ1tE5VH+m9NH2Br4C7C\nNVkN7ApcBOwEHFEzy0Q1kP4LYGYGXAq8D/SrsTl1izoHui5T3f2pTjyfVaVSs77uvqyDdfQGzgbu\ndPdDYva1ZtYN+KGZXe3uizpqq6grpP/W7OHur8U6l1SgPlH/yAdauBfYzN0XJOq9CngG+AmgzoGu\nh/QfcfeFhM6AJFeb2WLgJDMb7+7zOnIOUXdI//k5HvgEcA1wcgXr7VJoWkGGMbMjzOxJM1tmZvPN\n7GYz2yRVZvfEcOQVZjbXzC6MP7hzZSYRegxJDF9aFbf3itt7pOrdPOYflcj7dRzu9kkzmxIbrhsT\n+3c2s6lxWNxSM7vPzNINXj6+AAwErkjlX074N/XLxVwv0bXIkP7JdQwIkSQrPuDu/0h2DMS8D4Ep\nwCZmpn+QMkhW9N8Gr8bXdTtQh2hQsqZ/M1sPOIswqkJ/CLaBRg50XQaY2aBkhrvPz703s+8T/jG5\nBfgVsAEwDrjfzEa4++JY9D8JQ9KuAOYThqD9D6HnLfcv/C+BjYG9ga/RugfRYyoGJ2hyGvAg8G1g\nWbT3i4Qvck8CZxKGxR0L3Gtmu7v7k23UOyK+Nqfym2M9I4CbirRRNAbSv8g68oH22SjWX8l/pkR9\nIP2nMLMeQP/4eT4X658DvFykfaJxkP7X5mzgTeBq4EdF2pRN3F2pCyXgaILTpNOqRJnNgJXAd1LH\nbgN8CHw3kdcrzzm+A3wEbJLIuzR5jkT+nsAqwrDmZP7m0a6jEnmTYtmz89QzE7grldcL+Cdh+FRb\n1+RS4MMC+94Gflvr+6ZUmST9F3WNlgDX1fpeKVUnyQeKvk5DCF88J9X6nilVLkn/bV6bQ1LX5DFg\n21rfM6XKJem/4HXZLn7m0XF7YjzXwFrfs3pMGjnQNXHCEJ+XCuz/KqFn77ZUz+K8eMwXgJ8BuPsH\nuZ1m1pfQg/gIYUrKCOBflTae0Au5BjMbDmwFnJWy14C/0n4wnT6EB14+VsT9ousg/YusIx9oAzPr\nQwhOuwyY0DFTRR0i/efnXsK/u+sCo4HtCVMrRddC+l+bXxA6F/5aMSu7MCV3DpjZwcAfPczXS+b3\nAA5y91srZZzoEE944WAkQwiOnW8omZP4IW1mmxLm6OwPrJcqN6AyprbiI3dPP2y2iq+/KXDMajMb\n4IWDCi4HehbY1zvuLwrpv2GQ/quEfKBhkA/kwcyagN8Bw4Cx7v5WKcZJ/w2D9J/C3d8hdBAA3G5m\nEwgrdgzxEgISygcaAuk/YmaHAKOAbStkY5ennJEDNxPm6aUfJP3jPj0U6p8mwnCesfE1zfuw5kvU\nXwi9zOcQhvUsJcw1up7iAloWmmvUrUD+B3nycuf5NvBsgePeb8OGN4FuZra+u7+by4wN2SDgjTaO\nTSP9Nz5Z03+lkQ80Pln2gWsIy24d7u73F3lMEum/8cmy/pP8HvgpYanPX5VwnHygscma/s8jjBT7\nyMw2j3m5jo7NzKyXu7/ZxvGZo5zOASP/zd4IWJwnX9Qf/yTcxznu3lYgms8SeuyOdPff5jLNbO88\nZQs9ABbGc6Wj4W5RtLXBXoAl7n5vmyXz80y0YUdgaiL/c4SHzjMl1CX9Nz5Z03+lkQ80Ppn0ATM7\nnzAn9+QO/Lsp/Tc+mdR/HnJTKkv9B1g+0NhkTf+bAocTgiWmeYrwG2CHMurtshS9lKGZPWJmDxME\nMMXMHk6kx4AHgL9Vy1BRUW4n9BZOzLfTzAbGt6via1onp7D2g2BpPLZ/Kv/VWM8eqfxv5qmjEM2E\nh8NplmfJKTNbv53j7wUWACem8k8k2H1XewZI/12KrOm/IsgHuhSZ8wEzO53wz9NP3f2yIs+bPF76\n7zpkSv/pqPUJvhFtKGqlD/lAlyFT+ge+AhwUX3Ppd/H8RwCnFmlHZihl5MB98XUUIRjF0sS+DwkB\nJH5XGbNEB7G2drr7K2b2A+B/zWxL4I+ECOafJDjNVcCFwIsEh/y5hbVPFxMCmeRbE7c5nvdSM5tG\niFr6O3dfbGa3AePMjFjffoRlU4rC3d3Mvk5YxuQFC2uqvk4Y2vQFwnqlB7Zx/Aoz+yFwmZndSlgm\nZQ9CT+L33P29Isy4L75K//WP9J/CzPYjBJ8yoAewvYWljAD+5O7TizDlvvgqH6h/5AMJzOwg4Fxg\nFjDTzNL/IN0d52O3xX3xVfqvf6T/1hxhZifEz/kKsA4whhCc8P/c/b4iTcmVkw/UN9J/6+P/L51n\nZrklzqe6+4JibckMpS5vABwP9K7V8gpK7d6fowm9dDsUUfYrwP0Eh18MvABcAgxJlBlK+DG9iLDs\n35XAZ+I5kkuQNAEXA28RljhJLpsyiDAHbQnwLnA5sHWeOiYBi9qwdzvCvKF5hCjTrxDmt+1V5LU5\nDphBCEA4C/ifMq6v9F/HSfoveGxuiaB86aj2jk/VJR+o4yQfyHtcbtmqQmmP9q5Voi7pv46T9J/3\nuJGE9exnx+MWA08Q1rVvKuMaywfqNEn/RV8nLWXYRrJ4kUrCzD4WRfUp4BfuvtDMPgO84+5vl1yh\nEA2E9C+yjnxAZBnpX2Qd+YAQXZdyljLchhC98iNC8JEbCAEnjgA2BI6tpIFC1BPSv8g68gGRZaR/\nkXXkA0J0bYoOSJjgYsKwjs2BFYn8ycBeFbCpw5jZRDNbnUozEvt7mdnlZvaumS0xs9+b2cdrabNo\nGKR/kXXq2gekf1Fl6lr/IB8QVaeufUD6F6JjlNM5sBNwqa89H+FfwOCOm1QxphN6MAfHtHti38XA\nlwmBNfYANgb+0NkGioZE+hdZpxF8QPoX1aIR9A/yAVE9GsEHpH8hyqTkaQXASmCtpSQI847qKeLj\nR54n+nBcZuO/gEPd/f6YdyzwDzPbyd0f72Q7RWMh/Yus0wg+IP2LatEI+gf5gKgejeAD0r8QZVJO\n58BdwPfN7LC47Wa2EXAOcEfFLOs4W5nZ64QhT48AE9z9NULU1u7AX3MF3X2mmc0FdgHWejBYWCN2\nDDCH1kOoRPZ4BDjfzL5HGHmzjZkNISz78ncz26FK5+0NbAFMc/f5RZSvmP5BPiBaUQsfkP5FvaA2\nQD6QddQGCNH5lOoD5VPq8gbAQMLSF/MIwUheIiwN9wiwTq2XX4g2jiEMF/oM8CXgIcISLv2Aw4Dl\nsdxJMX85YYmNXxeo73DAlZTqIB3e2fqXDyjVUaqY/lM+sJqwbvLnpH+lOk9qA5SynNQGKGU9tesD\nHU0ljxxw9wXAnma2N2HNyY8BTwFT3H11qfVVA3eflticbmaPA68CBxN7/MzsEODnwH8TegrvBw4x\ns9Pc/d1UlXMARowYwTrrrNNqx5gxYxg7dmyb9px66qlcdNFFZX+eSlOsPQ8//DC33XYbb775JvPn\nz2fhwoW5h2QHMIK2C20Xym8iPLuLqSv52gfYBJjVRl3p7SHAx4FL4/YU4IeEALw3EkbU5diYMJXt\nLkLQ3rcIy6+uBs4C9k2UzdVzB3A7IcDvXbSeovcWYSrcx4Bc5/vbwMxcgTm0QxX0v+a8pfjAyJEj\nyX9/DdiA0L8YuOOOO9hss83a+2gVo1F9srMoxp7HHnuMl156iWXLljFs2DB23313mprKCWPTmqlT\npzJt2rRWeUuWLOHpp5+GCukf1vKBbxOcdpqZfbqjbcAFF1zAzTffTP5nVO51NWA0Nz/Z3keqOI2o\nt86mPZuqpX/I7wPPPPMMixcvhoZrA5Lk9J//O8DQoUO56aab2vt4FaHeNNeI9qgNCOTTf4v2C3/f\nbG5ubu+jVI1G1FtnU0ubOtoGdJjO+Ce/HhKh8fsp8AVgFfAEcEli/xzgPeCMPMfuAHhzc7OXw/77\n71/WcdWiPXtefvll799/oENTqrfKEnmW2k6XbcpzTDeHAQ43OsyNrwNSx3ZL1fNFh6tiuf5xO1d/\nuq6BDsMd1nPoFdO6MT99zuFxe3ieetZ1+LiDxzQ3nnPdWHf6nPs63BDLnJ+wfW6ijmQ9Uxzui+9v\nTJW5IXH8rNRxOLCDd7L+vQwfaLmXvRx65Ln+re9zZ9NoPtnZ1Js9zc3NldZ/f+DRnA9E/Z9MCKjV\n4TZg7WfU3vE192zqFp+NeL9+/Sp6rYqh3u5vvdnjXn827bHHHg3VBgwbNizhB7m2oFeeNrR1e9BZ\n1Nv9lT1t05htQPI7795r6b2W1Nv9rTd73OvPpo62AaWkkkcOmNl/F9jlhB65l4FH3YM3dTZmNgdI\n/gXpwERCoJTrgWbCg2EksJ2ZfQX4XTxmCmHOUabZeefdYu/UOsDlhGCuWxEuWz9gS+AfsXRfwmo2\nc4FfxLIPAP9DmNa1OJbZEHgx1ve1eOzXCLfnSGBXwj/usxPnfAAYR5hmcxlwNPA0cD5weht15fYD\nfED4tz9fuX8Az7Sx/6X4ue+P+94jxLFZCvw55h0ITAK2iduzE1fygUS9JOoZAjxG6EE+KZ5vz7j/\nZOCLwL0EV0qevzgK+MAKWvT/ETAC6GNmywkBhDai4vpfFdMI8l9fAGPTTTep3ClF1bn66qvz5psZ\nvXv3ZsiQIYwaNQozy1uu2hSp/28QIm6PMLP/IOj/EcLa3RXygdW0PKN+x9raDyxdurQypxOdQr3r\nP9oyhxq3AS+++GJ8t5rwHWBZ3L6Wwu0B3HPPPXzpS1+qhAmiStS7D9RPG7CKYtoAIeqNcgIS/hgY\nQPjFtjzm9SE43oeE3rgXzWxvd3+jIlaWhgN/J/xCfJ3g8OMJD4Nb4v6PgB7Azwgt1v8Sfom9QPhV\nmlmmTZvG/Plvx63cQ+1awo9sgB/Q8sMb4Oy4XegHdrIMrH1594yv7wPPtVHP4YQvGZcSwl60VVd6\nudpC5R5tZ/8jhB/x3wK2JnQm3AncTJB7H1qGR14cX6+Ir8MJHRvJH/7jCD/8H0vsW0LrhmJf4ABC\n50A/wvX4FiWyLnANcCVh3sPp8QPc4u6Lzew2wsVdAhxB6E0fDHQDBpV6sjRrfyFIV7lnq625c+d2\n9JSiE5k4cSKLFi1ixYoV9OnTB4Dly5fTu3dvevbsyeLFixk2bBh/+ctf2HjjjWthYnv6/w1wHmE8\n80TgGIL+tyPM4RlaOVMKPVtEo9IA+ocatwFrsyzxvm2fuPzyy9U5UOc0gA+oDRCiA5TTOfAtwt+b\nJ7r7CwBmti3hl+QltPzCuxA4tEJ2lko/wkNhEPAOobPg6+4+38xOJHRq9CF8jm6ECentTpQ69dRT\nGTBgQKu8ww47jMMOO6zAEY3HY489ltjKPdTuS+Slf3jnttt6ACaPKfRveu4HZaF6bk/sz/VJFapr\nHq0pVG5UO/uPjq8bAocAZxJi2TxNkNe2BKkfSesYB8sIrjWc1j/8mwg/+u8FesU8I8QXOAHYnzAC\nIzc4J3z27t170KtXv1L+YexOmFt3JC36P9ZbopvOiK9bA78GpgLT4kV4qa2Ky/OBdFDVlpEQ3/jG\n19s6nahDLrvsMi655BKuvPJKtt12WwBeeOEFTjrpJE4++WS22247jjjiCMaPH88tt9xScv0333xz\nnK/fwqJFi0qpoj39P0foyO5NCAIyhaD/8YTev4KUrv9CzxbRqFRb/5DfB6ZPn15KFXXWBiRp2yeW\nLFlSZD2iVqgNUBsgqksF2oCOUeo8BELDsdZ8B2BH4KX4fnfgjWrPiShg32zgDeBdQqDE04Buif3X\nEyLCrQQOiHl7Ecb//Ba4I0+dHYo5cNNNN5V1XLVoy56pU6cm5oLn5sNfk8hLzqlPbrc1dz5XZriH\nOfo3xPmGN8T5V3s6/Lydevqm9u+bp658MQdyc77S50zHHEjvt3i+UQ6/TMwTa3I4K1W2ycO8+nyx\nF5LxGvLFZMDNeqy1b6uthvp3v/tdv/vuu93d/ac//Wluf7tzjYrwgd8QhmEckMjbK+bdWaDOon2g\n9Wfs5dAvz/UNn79WNJJP1oK27BkyZEheHTzxxBM+ZMgQd3d/8MEHfaONNqqYPaXMN62HNqBF/7ln\ny6XxNfds6t/K5zubRtJbrShkUy30795YbYB7sh3ItX39vSXmQP72APBJkyZV9LoVot4010j2qA1o\nH1rFHEi2AfURc6CR9FYr6s2mUtqAjqZyfnwvz2cYYQ7/svh+C+D9ihsL3yMsSbIUWFCgzJmE8eDL\ngEXR3gsS+6cBfyKMY/+IMGrgu/Gh8CZwep46O9Q50GgMGrShQ/fUj+Zeibzh3jq4UL4f/es5bJAq\ns258Tf4Q7u7hx3dbP9RzAY2SX7afy1NXMjhiLmBivh/s6aCH6e18P/AL/ehPBl5sva979+6+4YYb\n+uDBg3348OF+0EEH+ejRo/2UU07xu+++26dMmeKzZs1yd/e7777bx40b52efffaavCSJhvGyCvnA\nAsLwuRUJH3DgvAJ1ltgw5gKu9cxz7WrbMSA6Ru/evfPq4Mknn/Q+ffq4u/vs2bMrGmivCvqvehsQ\ndN4zj/Z7Oqyz5pkmGota6N+98doAd/eWgMTdo+bztQc9Y7lu3q1br4peM1Ed1AYUB2u+46a/8+r7\njyidjgblLCWV8wN9GmGy9taJvK2jI06N218GXqiYkXAOoUc7l3IXKLe9Cvg0Ydz289HGzxLWOl0c\nnb9HrOt+Qm/hn+ID4/q4vRpYCGyQ5/yZ6hx45ZVXqrRaQVv1FftDPV8d+X60VzY1NTV5nz59vEeP\nHt6zZ08fNGiQb7311n7ggQf6uHHjfPjw4T5mzBifNWuWn3nmmWt6zyvFMccck7SnEj7wUdT9GcBP\n4vEOjPKKdA60fT9F47LPPvv4zjvv7DNmzFiTN2PGDB81apSPGTPG3d0nT57s22yzTcXOWQX9V70N\naPt51ZIvGota6N+98doAd/f11luvwHeH/G36008/XdFrJqqD2oDiyP8dqPajBkRj0pmdA+XEHPgv\nQmC/F8wsNwm6L6En77i4/SHwnTLqLsQFhJDwOb5CGEWwUyLvFcJDYBjwBQ9rlD5vZhfSslrBi4Ro\ndsvc/UAz+yahYWwiTP7+iru/U8iILMQcANhyyy1ZtGg+99xzD1dccQWzZ8/m7bffZt68eaxenVuv\n1VOvq1O1JMsVKpPcl2NVgXryb/fr15fddtuNYcOGsd9++7HFFlvw8ssvM2TIELbaaqs2PmV5vP76\n6xx66KE89NBDfPDBB7z22mssW7aM3XbbjWuvvZZPfOIT7LrrrnnWdy6efHON3n13zZK7BxHmi3bU\nB5YCEwj635DQg96flmUo8lKsD3hoTNcKTuievr+i0bjuuus49NBD2XbbbenXrx9AKx8A6NmzJ+ee\ne25Z9XeS/qveBuT3gZZnYW6/aCyqrX/oGm0AwIIFC4CcD6S/M+RYzaRJkzjmmGPaOq2oI9QGlN8G\nuH/U5mcXAioSd6NDWLlfUMxse1oies5092crZlXhc36PMCphR0IPYFNq/4+BH7F262PA7u7+kJnN\nIqy9N50QUW4uoZHczt27FTjvDkBzc3MzO+ywQyU/kmhQnn32WWbOnAnA0KFD2X777at6vqeeeirX\n4XA5YQkq+YCoKZ3pA9K/qDfUBoisozZAiM4j4QMj3f2pap6rpJEDZtaDsDD8V2NnQNU7BFL0AG4l\nDBEak2f/dvH1B4QhRTsRhsutTwgJD+Ez9ySEi/934KuE0RBuZr3c/QOEyMPKlSsZPnw4f/jDH9h+\n++2r/mWwAN2RD4gaUQc+IP2LmlEH+gf5gKghdeAD0r8QVaakzgF3X2lmA9svWTxmdg5tT0FwQnyD\nWe7+43jMtgXK5sYsTwC+T4hYehFwdqLMAEIE0/nA7+P7V4BPtmdrVqYViPz06NFjzTDJSjJhwoS1\nht+lR/TsuOOOubdXu/tT8gFRC6rhA9K/aBTUBsgHso7aAOlfVJ9aTysoJ+bAVcB4MzvB3dMTwssh\nHU8gH68UWdcLhB7A9wjrly4lxD+AEIEUoBswx933zB1kZtcAW7bXW3jRRRdpOFHGOf7447nwwgv5\n5S9/SVNTU/sHFMFpp53Gscce22aZ9957j5133rmY6uQDoqpU2gekf9FIqA2QD2QdtQHSv6gu+Tqc\nEtMKqk+pEQyBmwnDeeYQIn3elEzVjJ5ImCN0TTy/Ay8RlizJRSAdS+g1PAFopiX66IeJMvMJ65/+\nJyHwznJCIB4HehU4b4dWK6i3tTLrzR73+rOpkD2HHnqor7POOr755pv7AQcc4IcddlirVC0S65ve\nQegsWxm31+jfW/vAc4QlqhYSovQuS/nAwoT+c1OEquID9XZv3evPpkaypxY+kIjS23D6d2+s+1sL\n6s0ed7UB6aQ2oLo0kj1qA0qnke5vLag3e9zrz6ZEG1D11QrK7fK7i5Y1RC2VSsLMzjGz1W2kVWb2\n6Vh8d2AQLbEOfgF8Czgvbt9NiGJ6MfA6cCKxYwA4OpZ5jzD36Bbgj8C5QJ+4b1ip9hdDemhIrak3\ne6D+bGrLni9/+cvssssu9O3bN18DUjITJkygqampYOrWrRt33HFHrnhfwtKef43b3wKuNrN+cfth\nQofYQOAQ4ApCLznuvjKWWQ6sS+gkO4jgK9ux9tIQFaHe7i3Un02NZk8lfaAY/b/66qu54g2nf2i8\n+9vZ1Js9oDagkjTa/a0FjWaP2oDSaLT729nUmz1QfzZNmzat085V8rQCd6/0xJpSphXsDRxASyfE\nxfH9fwKnuvtqM7uZ0JP4RWBn4E7CEKPxhFEH/yJ87gHAKXH7WkIwkuOBbxYyYuTIkey///6t8jTX\nKFtU42GRHlI3efJkJk+e3KrMM888k3u7AvglLT4wEDiG4EMPAF8DlhA60H5LGFL3IjDMzHrExrEn\noVFcl9AL/y9gHjDA2xlSJx8QlfaBYvR/0UUX5d5K/6KmqA2QD2QdtQHSv6gu+WIOTJ8+vdPOX07M\ngYri7vMJQ3yKKXsscKyZHQ1c5O4DzexsYJ9EsaHAZHf/dwAz+37cP9TMBhBGPOwCnObuv4hlbiIE\nLdmlPRvuvHMKsEprVIuKMWjQIAYNGrRme/z48YwfP75VmT333JMHHngA4Mce1uY9mhBk5wpgH3d/\nIBYdBdyX0z+smUs3jNDjvogQoXelu49MlHmeMD+vXeQDopIUo//EXDvpX3Q51AaILKM2QIjW5Otw\nSrQBVaeszgEz2w84GNiM0AO3BnfftQJ2FTrvpoRews2BbtGOccAZCbtGAtPN7FOEToEJwHXAScBg\nQo/jd4DRZjYNGA38B3Bl/ExtcBZhoMP7Ff5kopGYPHkyt956K3PnzuXDDz9ste/hhx+u9uk3NLPt\nCT7Qg6D/n5tZP3dfCnwWWBQj+a4g+MBhhHlKgwkNY2+gp5mdS/CN0cDWcV87yAdETX1A+hc1R22A\nfCDrqA2Q/kXXpeTOATM7keAZNxH+ab8JGEJwyF+VY0SxyxkSfugflci/M+6bEbdXEjostgL2BV4m\nTB14iDAvCXefY2Yro73PEIYTHUfL+qf5iL2JWwKnAT/EzGhubi7q8y1atIinnnqqqLKdQb3ZA/Vn\nUyF7brvtNi6++GLGjh3Lww8/zNixY3nttdd4+eWXOeigg8r6DJdeeinXX399wf1mxtChQ3ObJwH/\nFt83EfT/I+BvhCF1q4HPEObdGcEHziGs85vkZ4QpOuMIPnADIZBPIcr2gXq7t1B/NjWSPZX2gWL0\nn1jmquH0D411f2tBvdkDagPyoDagijSSPWoD1AZUmnqzB+rPpiVLluTeFjXCpUOUGsGQEN3zyPh+\nCfDJ+P5nwMXlREUkBBn8dDupe6L8xsBMYFKeuq4Hbk/l7UWIXjogbr8KjEuVORN4uoB9hxMeQEpK\ntU6H08n6lw8o1VGS/pWynuQDSllO0r9S1tPh5fzWLiWVM61gc0LvHIQhO+vE99cSeupOKbVCLyHu\ngJl9ArgXeIIQRDDNI8DZZtbN3VfFvH2Ame6+KFFmNGG1gxxfivn5mEYIcjKH8JlFdnmIEODyLUK0\n3BMIS2puCvyaoKtq0BvYghBk5290rv5BPiBaqIUPSP+iXlAbIB/IOmoDhOh8cj5Q/WULSu1NIDjG\n8Pi+GTguvv8isLCaPRmE3sKXCEsWbgxsmEuJMv0JUUivB7YhLGXyfs7OWGYX4APCCgZDCT2GK4Bt\nqt0bo9TYSfpXynqqlQ9I/0r1kNQGKGU9qQ1QUuraqfQDwnIhP4jvT4lOdyfwLvCbqhoLRxOGBiXT\namBVqtxngPuBZcBcwsoE6bq+SljeZDnwHDCm1jdDqf6T9K+U9VQrH5D+leohqQ1QynpSG6Ck1LWT\nuTulYGY9CfP/l8XtY4BdCb15l7q7htuILov0L7KOfEBkGelfZB35gBBdm6I7B8zsR8AFuYeBEFlC\n+hdZRz6ZeZ/2AAAL4UlEQVQgsoz0L7KOfECIbFBK58AqYCN3n1ddk4SoP6R/kXXkAyLLSP8i68gH\nhMgGTSWUtapZ0YUxs5PMbLaZLTezR83sc1U4x0QzW51KMxL7e5nZ5Wb2rpktMbPfm9nHU3VsamZ3\nmdlSM3vLzM4zs6L1YWafN7P/M7PX4/kPyFPmJ2b2hpktM7N7zGxIav96ZvZbM1tkZgvN7Boz65cq\ns52ZPRCv56tmdno59pjZpDzXbEohewi+cmEH7JlgZo+b2WIze9vM7jCzT6fKVOQ+mdleZtZsZivM\nbJaZHZ3Pps6gM/Qfz1NTH+jq+jezhQQf6FumPdK/9N9p+i/GplJ9ALUBZSEfaFVGbYD0L/03aBvQ\nEZsaTv/FBicgBP3YoNZBEhopEaKkrgCOAoYBVwELgPUrfJ6JhIAqGwAfj2lgYv+VhOiyewIjCEtO\nPpjY3wQ8T1ge47PAGGAecHYJNowFfgIcSAgSc0Bq/3fiZ9+PECzmj8A/gZ6JMn8GngJ2JMxfmwXc\nmNi/DvAmIQrt1sDBwFLg62XYMwm4K3XNBqTKJO1ZHe0t154pwJGx3GeByfGe9KnkfSIsc/I+cB4h\nCu9JwErgS11V//XgAxnQ/66E9XVvk/6l/3rXf5V8QG2AfKBhfKAK+lcbIP03jP6r6ANl2dRo+i9F\neKuBhfHmFkyd7Xz1nIBHgUsS2wb8CzijwueZCDxVYF9/wpItByXyhsb7uVPc/rconvUTZY6P97t7\nGfaszuOEbwCnpuxaDhwct7eOx41IlBkDfAQMjtsnEqLhdk+UOQeYUYY9k4Db2zhmWNKe+H4JoXHM\n+cHSmJ/0gXbtifWtH4/dvZL3CTgXeC51rpuBKV1V//XmA11R/4l62tP/culf+q8n/VfKB1AbIB9o\nUB+ohP4T9agNkP4bSv8V9oFK+WRd67+UaQUQxHdqO0kAZtYDGAn8NZfn4S79hbDGaqXZKg6d+aeZ\n3Whmm8b8kUD3lB0zCcu75OwYBTzv7u8m6psGDAC27ahhZrYlMDhlw2LgsZQNC9396cShfyE0RDsn\nyjzg7h+l7BxqZgPKMG2vOLznRTO7wswGJvbtkseeHxB6H68naP1pQo9iUv/F2rNu/GwL4nal7tMo\nwnUjVaYamitIDfQPdeoDXUj/EBqqtvR/cZH2SP+VR/ovHbUBVUI+0EId+4DagCoh/bdQx/qH0n2g\nUjbVtf5L7Ry4xd2vbyuVWF9XZn2gG/B2Kv9tgpNUkkeBYwg9WicAWwIPxHkxg4EPoyMWsmNwATuh\nMrYOJjhBW9diMGF4zBrcfRXBcaph558JQ72+CJxBGMYzxcxysTXWsofQ+zYf+EfU+lLgiZT+27Un\nnuNi4O/uPiNRvhL3qVCZ/mbWq5BNVaAz9Q/17QNdRf/Qvv5/05490r/0X+D8nal/UBtQbeQDLdSj\nD6gNqC7Sfwv1qH8owwcqYVMj6L97sQUJN1bUIe4+LbE53cweB14lzH1ZURur6ht3vzWx+YKZPU+Y\n/7QX8Ld8h1Tw9FcA2wC7V7DOTCMfKI0y9A+V8wHpv8JI/6WjNqBrIR8oDbUBXQvpv3TK9IFKUPf6\n12oF1eNdwvDDDVP5GwJvVfPE7r6IEDRjSDxXTzPr34YdbxWwEypj61sE/bR1Ld4iBANZg5l1AwYS\ngn1U1U53n024Z7noqWl7jOAvHbLHzC4D9gX2cvc3Ers6ep/as2mxu3+Qz6YqUTP9Q935QFfQP4TP\nsF5H7JH+pf82zl8z/YPagCogH2ih7n1AbUDFkf5bqHv9Q3E+0FGbGkX/RXcOuHuTa23TonH3lUAz\nMDqXF4eSjCZEoKwaZvYx4FOEACDNhOAZSTuGApsl7HgE+KyZrZ+oZh9gETCDDhId7q2UDf0Jc3aS\nNqxrZiMSh44mPFAeT5TZIzpn0s6Z8UFYNma2CTCIFgdrZY+7NwHbdcSe+FA4EPiCu89N7e7offpH\nosxoWrNPzO80aqn/eK668YGuoP/I2I7YI/1L/4nz15X+o01qAyqIfKCFRvABtQGVRfpvoRH0H20q\nxgfKtqmh9O8Vjpip1CpC5MHAMlovYzKfCi8JCZwP7AFsTlhq4x7CHJNBcf8VwGzCUJmRwEOsvTzG\ns4T5N9sR5iy9DZxVgg39gO2B4YSANafE7U3j/jPiZ9+fsATHH4GXaL2MyRTgSeBzwG7ATOCGxP7+\nhAfd9YQhOYcQluw4rhR74r7zCA+mzQmO9CTBuXpUyZ4rCBFFP0/oxcul3qkyHbpPhGVMlhAilg4F\nvgl8COzdVfVfDz4g/Uv/0n/96F8+IB/Iug9I/9J/lvVfbz7QaPrvVEfJYoo3Zg5hyY5HgB2rcI6b\nCcujLCdEtrwJ2DKxvxdwKWG4zBLgNuDjqTo2Jay7+X4U27lAUwk27Bmdb1UqXZcoc2Z0omWE6JlD\nUnWsC9xI6AVbCPwK6Jsq8xng/ljHXOC0Uu0BegNTCT2ZK4BXCOuLblBFe/LZsgo4qtL3idBANEc9\nvAQc2ZX1Xw8+IP1L/9J//ehfPiAfyLoPSP/Sf5b1X28+0Gj6t1iREEIIIYQQQgghMkqpSxkKIYQQ\nQgghhBCii6HOASGEEEIIIYQQIuOoc0AIIYQQQgghhMg46hwQQgghhBBCCCEyjjoHhBBCCCGEEEKI\njKPOASGEEEIIIYQQIuOoc0AIIYQQQgghhMg46hwQQgghhBBCCCEyjjoHhBBCCCGEEEKIjKPOASGE\nEEIIIYQQIuOoc0BgZpPMbLWZrYqvufefrGD9t1eiLiGqgXxAZBnpX2QZ6V9kHfmASNK91gaIuuHP\nwDGAJfLeqY0p+TGzJsDd3Wtti+iSyAdElpH+RZaR/kXWkQ8IQCMHRAsfuPs77j4vkdwCE8zsFTNb\nZmZPm9lXcweZWZOZXZPY/6KZjUvsnwgcDRyY6Incw8z2jNv9E2W3j3mbxe2jzWyhme1vZi8AK4BN\n476vm9kMM1seX0/srAsluizyAZFlpH+RZaR/kXXkAwLQyAHRPt8DDgf+G3gZ2AO4wczmufuDhA6m\n14CvAguAXYGrzewNd/89cAGwNbAOLT2SC4DdgHw9f+m8vsAZwHHAfGCemX0NOBM4CXgGGAH8ysze\nd/cbKvOxhViDfEBkGelfZBnpX2Qd+UDGUOeAyLG/mS1JbE8BjgQmAKPd/bGYP8fMPg8cDzzo7h8B\nP04c96qZ7QocDPze3Zea2XKgp7uvGZ5klhy11CbdgRPdfXri2DOBb7v7nxLn3BY4AdBDQZSLfEBk\nGelfZBnpX2Qd+YAA1DkgWriX4FQ5b10KDCH02N1jrb24B/B0bsPMTgKOBTYD+gA9k/s7yIepB0Jf\n4FPAtWZ2TaJcN+C9Cp1TZBP5gMgy0r/IMtK/yDryAQGoc0C0sNTdZyczzGzz+HZf4I1U+Q9imUOB\n84FTgUeBJYThPzu1c77VudMk8nrkKbc8tf2x+Pp14PHUvlXtnFOItpAPiCwj/YssI/2LrCMfEIA6\nB0TbzCA4/+bu/vcCZXYFHnL3q3IZZvapVJkPCT16Sd4hPBA2AhbFvBHtGeTu88zsDeBT7n5L+x9B\niA4hHxBZRvoXWUb6F1lHPpBB1DkgCuLu75vZBcBFZtYN+DswgBBEZFEM+vEScKSZ7QPMJsxP+hzw\nSqKqOcA+ZvZpQjCRRYSgJq8BZ5rZD4ChwPgiTZsIXGJmi4GpQC9gR2Bdd7+4Ax9ZiFbIB0SWkf5F\nlpH+RdaRD2QTLWUo2sTdfwicBXyX0IP4Z8LwotzQo6uA24FbCMOJBgKXp6r5FTATeBKYB+waA5gc\nCgwDngVOB75fpE3XEoYTHQs8B9xHWCZldhuHCVEW8gGRZaR/kWWkf5F15APZw9zzrSIhhBBCCCGE\nEEKIrKCRA0IIIYQQQgghRMZR54AQQgghhBBCCJFx1DkghBBCCCGEEEJkHHUOCCGEEEIIIYQQGUed\nA0IIIYQQQgghRMZR54AQQgghhBBCCJFx1DkghBBCCCGEEEJkHHUOCCGEEEIIIYQQGUedA0IIIYQQ\nQgghRMZR54AQQgghhBBCCJFx1DkghBBCCCGEEEJknP8HlFN5ckaT7FQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xfb6e780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.feature_selection import *\n",
    "\n",
    "fs=SelectKBest(score_func=f_regression,k=5)\n",
    "X_new=fs.fit_transform(X_train,y_train)\n",
    "print((fs.get_support(),feature_names))\n",
    "\n",
    "x_min, x_max = X_new[:,0].min(), X_new[:, 0].max()\n",
    "y_min, y_max = y_train.min(), y_train.max()\n",
    "\n",
    "fig=plt.figure()\n",
    "#fig.subplots_adjust(left=0, right=1, bottom=0, top=1, hspace=0.05, wspace=0.05)\n",
    "\n",
    "# Two subplots, unpack the axes array immediately\n",
    "fig, axes = plt.subplots(1,5)\n",
    "\n",
    "fig.set_size_inches(12,12)\n",
    "\n",
    "for i in range(5):\n",
    "    axes[i].set_aspect('equal')\n",
    "    axes[i].set_title('Feature ' + str(i))\n",
    "    axes[i].set_xlabel('Feature')\n",
    "    axes[i].set_ylabel('Target')\n",
    "    axes[i].set_xlim(x_min, x_max)\n",
    "    axes[i].set_ylim(y_min, y_max)\n",
    "    plt.sca(axes[i])\n",
    "    plt.scatter(X_new[:,i],y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In regression tasks, is very important to normalize data (to avoid that large-valued features weight too much in the final result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.26478671796 -4.95114422825 -5.82246078184e-18 1.49107035114 -4.87256850449 -8.10081500083e-18\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "scalerX = StandardScaler().fit(X_train)\n",
    "scalery = StandardScaler().fit(y_train)\n",
    "\n",
    "X_train = scalerX.transform(X_train)\n",
    "y_train = scalery.transform(y_train)\n",
    "X_test = scalerX.transform(X_test)\n",
    "y_test = scalery.transform(y_test)\n",
    "\n",
    "print(np.max(X_train), np.min(X_train), np.mean(X_train), np.max(y_train), np.min(y_train), np.mean(y_train))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## A Linear Model\n",
    "\n",
    "Let's try a lineal model, SGDRegressor, that tries to find the hyperplane that minimizes a certain loss function (typically, the sum of squared distances from each instance to the hyperplane). It uses Stochastic Gradient Descent to find the minimum. \n",
    "\n",
    "Regression poses an additional problem: how should we evaluate our results? Accuracy is not a good idea, since\n",
    "we are predicting real values, as it is almost impossible for us to predict exactly the final value. There are several measures that can be used. The most common is the R2 score, or coefficient of determination that measures the proportion of the outcomes variation explained by the model, and is the default score function for regression methods in scikit-learn. This score reaches its maximum value of **1** when the model perfectly predicts all the test target values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import *\n",
    "def train_and_evaluate(clf, X_train, y_train):\n",
    "    \n",
    "    clf.fit(X_train, y_train)\n",
    "    \n",
    "    print (\"Coefficient of determination on training set:\",clf.score(X_train, y_train))\n",
    "    \n",
    "    # create a k-fold croos validation iterator of k=5 folds\n",
    "    cv = KFold(X_train.shape[0], 5, shuffle=True, random_state=33)\n",
    "    scores = cross_val_score(clf, X_train, y_train, cv=cv)\n",
    "    print (\"Average coefficient of determination using 5-fold crossvalidation:\",np.mean(scores))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.998521032316\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.998347675549\n",
      "[-0.00203934  0.037831   -0.00643619  0.28580185  0.15658046  0.00879088\n",
      "  0.05961933  0.00366377  0.0145986   0.17708452  0.08886078  0.26907647]\n"
     ]
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "clf_sgd = linear_model.SGDRegressor(loss='squared_loss', penalty=None,  random_state=42)\n",
    "train_and_evaluate(clf_sgd,X_train,y_train)\n",
    "print( clf_sgd.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You probably noted the `penalty=None` parameter when we called the method. The penalization parameter for linear regression methods is introduced to avoid overfitting. It does this by penalizing those hyperplanes having some of their coefficients too large, seeking hyperplanes where each feature contributes more or less the same to the predicted value. This parameter is generally the L2 norm (the squared sums of the coefficients) or the L1 norm (that is the sum of the absolute value of the coefficients). Let's see how our model works if we introduce an L2 or L1 penalty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.998519126977\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.998346089899\n"
     ]
    }
   ],
   "source": [
    "clf_sgd1 = linear_model.SGDRegressor(loss='squared_loss', penalty='l2',  random_state=42)\n",
    "train_and_evaluate(clf_sgd1,X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.998523707937\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.998350925502\n"
     ]
    }
   ],
   "source": [
    "clf_sgd2 = linear_model.SGDRegressor(loss='squared_loss', penalty='l1',  random_state=42)\n",
    "train_and_evaluate(clf_sgd2,X_train,y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "_____\n",
    "## Support Vector Machines for regression\n",
    "The regression version of SVM can be used instead to find the hyperplane (note how easy is to change the classification method in scikit-learn!). We will try a linear kernel, a polynomial kernel, and finally, a rbf kernel. For more information on kernels, see http://scikit-learn.org/stable/modules/svm.html#svm-kernels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.997200877829\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.997130094148\n"
     ]
    }
   ],
   "source": [
    "from sklearn import svm\n",
    "clf_svr= svm.SVR(kernel='linear')\n",
    "train_and_evaluate(clf_svr,X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.996715736848\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.996467997134\n"
     ]
    }
   ],
   "source": [
    "clf_svr_poly= svm.SVR(kernel='poly')\n",
    "train_and_evaluate(clf_svr_poly,X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.998175246372\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.997722562258\n"
     ]
    }
   ],
   "source": [
    "clf_svr_rbf= svm.SVR(kernel='rbf')\n",
    "train_and_evaluate(clf_svr_rbf,X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 0.997507601395\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.997418409054\n"
     ]
    }
   ],
   "source": [
    "clf_svr_poly2= svm.SVR(kernel='poly',degree=2)\n",
    "train_and_evaluate(clf_svr_poly2,X_train,y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "## Random Forests for Regression Analysis\n",
    "Finally, let's try again Random Forests, in their Extra Trees, and Regression version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination on training set: 1.0\n",
      "Average coefficient of determination using 5-fold crossvalidation: 0.999983293438\n"
     ]
    }
   ],
   "source": [
    "from sklearn import ensemble\n",
    "clf_et=ensemble.ExtraTreesRegressor(n_estimators=10,random_state=42)\n",
    "train_and_evaluate(clf_et,X_train,y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An interesting side effect of random forest classification, is that you can measure how 'important' each feature is when predicting the final result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imp_features = (np.sort((clf_et.feature_importances_,features),axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.000 <-> CO(GT)\n",
      "0.001 <-> PT08.S1(CO)\n",
      "0.000 <-> NMHC(GT)\n",
      "0.208 <-> C6H6(GT)\n",
      "0.011 <-> PT08.S2(NMHC)\n",
      "0.001 <-> NOx(GT)\n",
      "0.099 <-> PT08.S3(NOx)\n",
      "0.000 <-> NO2(GT)\n",
      "0.003 <-> PT08.S5(O3)\n",
      "0.194 <-> T\n",
      "0.194 <-> RH\n",
      "0.291 <-> AH\n"
     ]
    }
   ],
   "source": [
    "for rank,f in zip(imp_features[0],imp_features[1]):\n",
    "    print(\"{0:.3f} <-> {1}\".format(float(rank), f))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, evaluate our classifiers on the testing set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficient of determination:1.000 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "def measure_performance(X,y,clf, show_accuracy=True,\n",
    "                        show_classification_report=True,\n",
    "                        show_confusion_matrix=True,\n",
    "                        show_r2_score=False):\n",
    "    y_pred=clf.predict(X)   \n",
    "    if show_accuracy:\n",
    "        print (\"Accuracy:{0:.3f}\".format(metrics.accuracy_score(y,y_pred)),\"\\n\")\n",
    "\n",
    "    if show_classification_report:\n",
    "        print (\"Classification report\")\n",
    "        print (metrics.classification_report(y,y_pred),\"\\n\")\n",
    "        \n",
    "    if show_confusion_matrix:\n",
    "        print (\"Confusion matrix\")\n",
    "        print (metrics.confusion_matrix(y,y_pred),\"\\n\")\n",
    "        \n",
    "    if show_r2_score:\n",
    "        print (\"Coefficient of determination:{0:.3f}\".format(metrics.r2_score(y,y_pred)),\"\\n\")\n",
    "\n",
    "        \n",
    "measure_performance(X_test,y_test,clf_et,\n",
    "                    show_accuracy=False,\n",
    "                    show_classification_report=False,\n",
    "                    show_confusion_matrix=False,\n",
    "                    show_r2_score=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}