06.4_Adding_a_category_to_capture_NA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a category to capture NA\n",
    "\n",
    "This is perhaps the most widely used method of missing data imputation for categorical variables. This method consists in treating missing data as if they were an additional label or category of the variable. All the missing observations are grouped in the newly created label 'Missing'. \n",
    "\n",
    "The beauty of this technique resides on the fact that it does not assume anything on the missingness of the values. It is very well suited when the number of missing data is high.\n",
    "\n",
    "\n",
    "### Advantages\n",
    "\n",
    "- Easy to implement\n",
    "- Captures the importance of missingess if there is one\n",
    "\n",
    "### Disadvantages\n",
    "\n",
    "- If the number of NA is small, creating an additional category may cause trees to over-fit\n",
    "\n",
    "I would say that for categorical variables this is the method of choice, as it treats missing values as a separate category, without making any assumption on their missingness. It is used widely in data science competitions and business settings. See for example the winning solution of the KDD 2009 cup: \"Winning the KDD Cup Orange Challenge with Ensemble Selection\" (http://www.mtome.com/Publications/CiML/CiML-v3-book.pdf).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "===============================================================================\n",
    "\n",
    "## Real Life example: \n",
    "\n",
    "### Predicting Sale Price of Houses\n",
    "\n",
    "The problem at hand aims to predict the final sale price of homes based on different explanatory variables describing aspects of residential homes. Predicting house prices is useful to identify fruitful investments, or to determine whether the price advertised for a house is over or underestimated, before making a buying judgment.\n",
    "\n",
    "=============================================================================\n",
    "\n",
    "In the following cells, I will demonstrate NA imputation by adding an additional label using the House Price datasets from Kaggle.\n",
    "\n",
    "If you haven't downloaded the datasets yet, in the lecture \"Guide to setting up your computer\" in section 1, you can find the details on how to do so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "% matplotlib inline\n",
    "\n",
    "# for regression problems\n",
    "from sklearn.linear_model import LinearRegression, Ridge\n",
    "\n",
    "# to split and standarize the datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# to evaluate regression models\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### House Price dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SalePrice      0.000000\n",
       "BsmtQual       0.025342\n",
       "GarageType     0.055479\n",
       "FireplaceQu    0.472603\n",
       "dtype: float64"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's load the dataset with a few columns for the demonstration\n",
    "cols_to_use = ['BsmtQual', 'FireplaceQu', 'GarageType', 'SalePrice']\n",
    "\n",
    "data = pd.read_csv('houseprice.csv', usecols=cols_to_use)\n",
    "\n",
    "# let's inspect the percentage of missing values in each variable\n",
    "data.isnull().mean().sort_values(ascending=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1022, 3), (438, 3))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's separate into training and testing set\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(data[['BsmtQual', 'FireplaceQu', 'GarageType']],\n",
    "                                                    data.SalePrice, test_size=0.3,\n",
    "                                                    random_state=0)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's create a variable to replace NA with the most frequent label or a random sample\n",
    "\n",
    "def impute_na(df_train, df_test, variable):\n",
    "    df_train[variable+'_NA'] = np.where(df_train[variable].isnull(), 'Missing', df_train[variable])\n",
    "    df_test[variable+'_NA'] = np.where(df_test[variable].isnull(), 'Missing', df_test[variable])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# and let's replace the NA\n",
    "for variable in ['BsmtQual', 'FireplaceQu', 'GarageType',]:\n",
    "    impute_na(X_train, X_test, variable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BsmtQual           24\n",
       "FireplaceQu       478\n",
       "GarageType         54\n",
       "BsmtQual_NA         0\n",
       "FireplaceQu_NA      0\n",
       "GarageType_NA       0\n",
       "dtype: int64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's check that data have been completed\n",
    "X_train.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BsmtQual</th>\n",
       "      <th>FireplaceQu</th>\n",
       "      <th>GarageType</th>\n",
       "      <th>BsmtQual_NA</th>\n",
       "      <th>FireplaceQu_NA</th>\n",
       "      <th>GarageType_NA</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>Gd</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Attchd</td>\n",
       "      <td>Gd</td>\n",
       "      <td>Missing</td>\n",
       "      <td>Attchd</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>682</th>\n",
       "      <td>Gd</td>\n",
       "      <td>Gd</td>\n",
       "      <td>Attchd</td>\n",
       "      <td>Gd</td>\n",
       "      <td>Gd</td>\n",
       "      <td>Attchd</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>960</th>\n",
       "      <td>TA</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>TA</td>\n",
       "      <td>Missing</td>\n",
       "      <td>Missing</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1384</th>\n",
       "      <td>TA</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Detchd</td>\n",
       "      <td>TA</td>\n",
       "      <td>Missing</td>\n",
       "      <td>Detchd</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1100</th>\n",
       "      <td>TA</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Detchd</td>\n",
       "      <td>TA</td>\n",
       "      <td>Missing</td>\n",
       "      <td>Detchd</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     BsmtQual FireplaceQu GarageType BsmtQual_NA FireplaceQu_NA GarageType_NA\n",
       "64         Gd         NaN     Attchd          Gd        Missing        Attchd\n",
       "682        Gd          Gd     Attchd          Gd             Gd        Attchd\n",
       "960        TA         NaN        NaN          TA        Missing       Missing\n",
       "1384       TA         NaN     Detchd          TA        Missing        Detchd\n",
       "1100       TA         NaN     Detchd          TA        Missing        Detchd"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's see how the new variable looks like, where data was missing we have\n",
    "# not the label 'Missing'\n",
    "\n",
    "X_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEmCAYAAABs7FscAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEv1JREFUeJzt3X+wZ3Vdx/Hnq+WXk6Ewe6MddvNiLdniCOKGP7IiyYGk\nXKyExaKtoSjDSe2HQTX9sNYoralMKvy59kPcwnTFRqMVIpsSFwFhUWQTCJgVVq0kq1WWd398z3a/\nXvfuvd+79+659/N9PmZ2vud8zjnf876Hy+ue8zm/UlVIktr1VX0XIElaXAa9JDXOoJekxhn0ktQ4\ng16SGmfQS1LjDHpJapxBL0mNM+glqXFH9F0AwMqVK2tycrLvMiRpWbn55ps/U1UTs823JIJ+cnKS\nHTt29F2GJC0rSe6by3x23UhS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIatyRumFoI\nk5e9r+8SALj3inP7LkGSvox79JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS\n1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcc08plhTfGSzpGHu0UtS4wx6SWqcQS9JjTPoJalxBr0k\nNc6gl6TGGfSS1DiDXpIaZ9BLUuPmHPRJViS5Jcm13fjxSa5Lcnf3edzQvJcn2ZXkriRnL0bhkqS5\nGWWP/uXAx4fGLwO2V9VaYHs3TpJ1wEbgFOAc4MokKxamXEnSqOYU9ElWA+cCbxpq3gBs6Ya3AOcN\ntV9dVXur6h5gF3DGwpQrSRrVXPfofx94FfDYUNsJVbW7G/40cEI3fCJw/9B8D3RtXybJJUl2JNmx\nZ8+e0aqWJM3ZrEGf5HuAh6vq5pnmqaoCapQVV9VVVbW+qtZPTEyMsqgkaQRzeUzxtwIvTPIC4Bjg\n2CR/DjyUZFVV7U6yCni4m/9BYM3Q8qu7NklSD2bdo6+qy6tqdVVNMjjJ+sGq+iFgG7Cpm20T8J5u\neBuwMcnRSU4C1gI3LXjlkqQ5OZQXj1wBbE1yMXAfcD5AVe1MshW4E3gUuLSq9h1ypZKkeRkp6Kvq\nBuCGbvizwFkzzLcZ2HyItUmSFoB3xkpS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa\nZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEG\nvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BL\nUuMMeklqnEEvSY0z6CWpcbMGfZJjktyU5LYkO5P8etd+fJLrktzdfR43tMzlSXYluSvJ2Yv5A0iS\nDm4ue/R7gedV1anAacA5SZ4FXAZsr6q1wPZunCTrgI3AKcA5wJVJVixG8ZKk2c0a9DXwX93okd2/\nAjYAW7r2LcB53fAG4Oqq2ltV9wC7gDMWtGpJ0pzNqY8+yYoktwIPA9dV1YeBE6pqdzfLp4ETuuET\ngfuHFn+ga5v+nZck2ZFkx549e+b9A0iSDm5OQV9V+6rqNGA1cEaSp06bXgz28uesqq6qqvVVtX5i\nYmKURSVJIxjpqpuq+g/gegZ97w8lWQXQfT7czfYgsGZosdVdmySpB3O56mYiyRO74ccBzwc+AWwD\nNnWzbQLe0w1vAzYmOTrJScBa4KaFLlySNDdHzGGeVcCW7sqZrwK2VtW1Sf4Z2JrkYuA+4HyAqtqZ\nZCtwJ/AocGlV7Vuc8iVJs5k16KvqY8DTD9D+WeCsGZbZDGw+5OokSYfMO2MlqXEGvSQ1zqCXpMYZ\n9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEv\nSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLU\nOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LhZgz7JmiTXJ7kzyc4kL+/aj09y\nXZK7u8/jhpa5PMmuJHclOXsxfwBJ0sHNZY/+UeBnq2od8Czg0iTrgMuA7VW1FtjejdNN2wicApwD\nXJlkxWIUL0ma3axBX1W7q+qj3fAjwMeBE4ENwJZuti3Aed3wBuDqqtpbVfcAu4AzFrpwSdLcjNRH\nn2QSeDrwYeCEqtrdTfo0cEI3fCJw/9BiD3Rt07/rkiQ7kuzYs2fPiGVLkuZqzkGf5PHANcArqurz\nw9OqqoAaZcVVdVVVra+q9RMTE6MsKkkawZyCPsmRDEL+L6rqXV3zQ0lWddNXAQ937Q8Ca4YWX921\nSZJ6MJerbgK8Gfh4Vf3e0KRtwKZueBPwnqH2jUmOTnISsBa4aeFKliSN4og5zPOtwEXA7Ulu7dp+\nEbgC2JrkYuA+4HyAqtqZZCtwJ4Mrdi6tqn0LXrkkaU5mDfqq+hCQGSafNcMym4HNh1CXJGmBeGes\nJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS\n4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXO\noJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuFmDPslbkjyc5I6h\ntuOTXJfk7u7zuKFplyfZleSuJGcvVuGSpLmZyx7924BzprVdBmyvqrXA9m6cJOuAjcAp3TJXJlmx\nYNVKkkY2a9BX1Y3A56Y1bwC2dMNbgPOG2q+uqr1VdQ+wCzhjgWqVJM3DfPvoT6iq3d3wp4ETuuET\ngfuH5nuga5Mk9eSQT8ZWVQE16nJJLkmyI8mOPXv2HGoZkqQZzDfoH0qyCqD7fLhrfxBYMzTf6q7t\nK1TVVVW1vqrWT0xMzLMMSdJs5hv024BN3fAm4D1D7RuTHJ3kJGAtcNOhlShJOhRHzDZDkncAZwIr\nkzwA/CpwBbA1ycXAfcD5AFW1M8lW4E7gUeDSqtq3SLVLkuZg1qCvqgtnmHTWDPNvBjYfSlGSpIXj\nnbGS1DiDXpIaZ9BLUuNm7aOXlrPJy97Xdwnce8W5fZegMecevSQ1zqCXpMYZ9JLUOINekhpn0EtS\n4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXO\noJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ17oi+C5B0eExe9r6+S+DeK87tuwRg\n/LaFe/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcYsW9EnOSXJXkl1JLlus9UiSDm5Rgj7JCuAN\nwHcD64ALk6xbjHVJkg5usfbozwB2VdWnquqLwNXAhkValyTpIFJVC/+lyQ8A51TVj3XjFwHPrKqX\nDc1zCXBJN/pNwF0LXsjoVgKf6buIJcJtMcVtMcVtMWUpbIsnVdXEbDP19giEqroKuKqv9R9Ikh1V\ntb7vOpYCt8UUt8UUt8WU5bQtFqvr5kFgzdD46q5NknSYLVbQfwRYm+SkJEcBG4Fti7QuSdJBLErX\nTVU9muRlwAeAFcBbqmrnYqxrgS2prqSeuS2muC2muC2mLJttsSgnYyVJS4d3xkpS4wx6SWqcQS9J\njTPohyR5bpI39F2HJC2ksX9nbJKnAy8BXgzcA7yr34oOrySnH2x6VX30cNWylCS5uKrePDS+Avjl\nqvr1HsuS5mUsgz7JycCF3b/PAO9kcAXSd/ZaWD9+t/s8BlgP3AYEeBqwA3h2T3X17awk3w9cDBwP\nvA34h14r6kmS9wLTL8/7Twa/H39aVf97+KvqR5JnAa8Hvhk4isHl41+oqmN7LWwWYxn0wCeAfwS+\np6p2ASR5Zb8l9WP/H7ck7wJOr6rbu/GnAr/WY2m9qqqXJLkAuB34AvCSqvqnnsvqy6eACeAd3fgF\nwCPAycAbgYt6qqsPf8TgBtC/YrBj9MMMtsOSNq5B/30M/mNdn+T9DJ6umX5L6t037Q95gKq6I8k3\n91lQn5KsBV4OXMNg7+2iJLdU1X/3W1kvnlNV3zI0/t4kH6mqb0myHG6EXFBVtSvJiqraB7w1yS3A\n5X3XdTDjGvTXVtW7k3w1g8cnvwL42iR/DPxNVf1dv+X14mNJ3gT8eTf+g8DHeqynb+8FLq2q7UkC\n/AyDR3uc0m9ZvXh8kq+vqn8DSPL1wOO7aV/sr6xe/Hf3WJdbk/wOsJtlcFHLWN4Zm+SjVXX6tLbj\nGJyQvaCqzuqnsv4kOQZ4KfBtXdONwB9X1d7+qupPkmOr6vPT2k6uqk/2VVNfkrwA+BPgXxkc+Z4E\n/BRwA/DjVfX7/VV3eCV5EvAQg/75VwJPAK7c3wW8VI1r0N9SVU/vu46lIMkGYHVVvaEbv4lBf2wB\nr6qqv+6zvsMtyauq6ne64RdX1V8NTXtNVf1if9X1J8nRwFO60bvG6QQsDI5i9h/RLEfjGvQPAL83\n0/SqmnFaa5L8E7Cxqu7vxm8Fnsfg0Pyt43Z0M3y0N/3I70BHguMiyXOASYa6e6vq7b0VdJhN+724\npqq+v++aRjGuffQrGATZuJ+ABThqf8h3PlRVnwM+153DGDeZYfhA42MhyZ8B3wDcCuzrmgsYm6Dn\ny//bP7m3KuZpXIN+d1W9uu8ilojjhkeGX/fIoAtn3NQMwwcaHxfrgXU1jof/Uw72e7HkjWvQj+We\n2Qw+nOTHq+qNw41JfgK4qaea+nRqks8z+B15XDdMN35Mf2X16g7g6xhcYTKuDvZ7UUv9hqlx7aM/\nvuueGHtJvhZ4N7AX2P+4g2cARwPnVdVDfdWmpSHJ9cBpDP7w//9VWFX1wt6K0kjGMuj1lZI8j6lr\nxHdW1Qf7rEdLR5LvOFB7VY3lIyGWI4Nekho3rn30kmaR5ENV9dwkj/DlJyCXRb+0prhHL0mNW/LP\naJDUryTf0N0ZS5Izk/x0kif2XZfmzqCXNJtrgH1JvhG4ClgD/GW/JWkUBr2k2TxWVY8CLwJeX1U/\nD6zquSaNwKCXNJsvJbkQ2ARc27Ud2WM9GpFBL2k2P8rglZKbq+qeJCcBf9ZzTRqBV91ImrPuvQ1r\nqmqcX0qz7LhHL+mgktyQ5NgkxzN4TMYbk4zNo7xbYNBLms0TurdtfR/w9qp6JvBdPdekERj0kmZz\nRJJVwPlMnYzVMmLQS5rNq4EPALuq6iNJngzc3XNNGoEnYyWpcT7UTNIB7X9RepLXc4C3KlXVT/dQ\nlubBoJc0k493nzt6rUKHzK4bSWqce/SSDijJtoNN91WCy4dBL2kmzwbuB94BfJjBC0e0DNl1I+mA\nkqwAng9cCDwNeB/wjqra2WthGpnX0Us6oKraV1Xvr6pNwLOAXcANSV7Wc2kakV03kmbUvVnqXAZ7\n9ZPAHwJ/02dNGp1dN5IOKMnbgacCfwtcXVV39FyS5smgl3RASR4DvtCNDgdFgKqqYw9/VZoPg16S\nGufJWElqnEEvSY0z6CWpcQa9lrQk+5LcmuS2JB9N8pwF+M7TkrxgWtt5ST6W5BNJ7kjyA4fw/ZNJ\nZrxCJcmZSSrJ9w61XZvkzKHxlUm+lOQn51uHtJ9Br6Xuf6rqtKo6Fbgc+K0F+M7TgP8P+iSnAq8D\nNlTVU4DvBX47yTMWYF0zeQD4pYNMfzHwLwyuX5cOiUGv5eRY4N8BkqxKcmO3t39Hkm/r2v8ryWuT\n7Ezy90nO6F5u/akkL0xyFIM3Jl3QLXsB8HPAa6rqHoDu8zXAz3bfeUOS9d3wyiT3dsOTSf6xO9IY\n9WjjNuA/kzx/hukXdus/McnqkbaSNI1Br6XucV0gfwJ4E/AbXftLgA9U1WnAqcCtXftXAx+sqlOA\nR4DfZPC8lhcBr66qLwK/AryzO1J4J3AKcPO09e4A1s1S28PA86vqdOACBneNjmIz8MvTG5OsAVZV\n1U3A1u67pXnzEQha6v6nC3OSPBt4e5KnAh8B3pLkSODdVbU/6L8IvL8bvh3YW1VfSnI7g1v4F9KR\nwB8lOQ3YB5w8ysJVdWMSkjx32qQLGAQ8wNXAW4DfPdRiNb7co9eyUVX/DKwEJqrqRuDbgQeBtyX5\n4W62L9XUXYCPAXu7ZR9j5h2bO4Hp/fHPYOrNSo8y9f/KMUPzvBJ4iMERxXrgqHn8WAfaq78Q+JGu\ni2gb8LQka+fx3RJg0GsZSfIUYAXw2SRPAh6qqjcy6NI5fYSvegT4mqHx1wGXJ5ns1jMJvAJ4bTf9\nXqb+EAxfjfMEYHf3R+SirraRVNXfAccxeAwwSU4GHl9VJ1bVZFVNMjgB7UlZzZtBr6Vufx/9rcA7\ngU1VtQ84E7gtyS0Mujr+YITvvB5Yt/9kbNft8wvAe5N8Evgk8NKququb/3XAS7t1rRz6niuBTUlu\nA57C1HNhRrUZWNMNX8hXPh3yGgx6HQKfdSNNk+QK4JnA2d3JW2lZM+glqXFedSMtkiRnA789rfme\nqnpRH/VofLlHL0mN82SsJDXOoJekxhn0ktQ4g16SGvd/wffx/8ADWecAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8cb93f8710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEmCAYAAABs7FscAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFWlJREFUeJzt3X+0XWV95/H3pwGBpYNCiWkWCU3axrZBBWpEbWlrpRZa\nrMGqGJw66ZQxjqLizLQMdGZWf610UWdNV6uV1ihqxBFMi0iEGSiNUH+MJYSfEjAlQ6DABBKlU1t0\nUOJ3/jg7crjk5p6be+49uc99v9bKOns/+9lnf/cSP+e5z95nn1QVkqR2fd+oC5AkTS+DXpIaZ9BL\nUuMMeklqnEEvSY0z6CWpcQa9JDVuoKBPcn+SryS5PcmWru3oJNcnubd7Paqv/4VJtifZluS06Spe\nkjSxyYzof66qTqyqFd36BcCmqloGbOrWSbIcWAUcD5wOXJxk3hBrliRNwiFT2Hcl8MpueT1wI/Af\nu/bLq+oJYEeS7cDJwJfHe6NjjjmmlixZMoVSJGnuueWWW75WVfMn6jdo0Bfw10n2AB+sqnXAgqra\n2W1/BFjQLR8L/G3fvg91bU+TZA2wBuC4445jy5YtA5YiSQJI8sAg/QYN+lOq6uEkzweuT/LV/o1V\nVUkm9dCc7sNiHcCKFSt84I4kTZOB5uir6uHudRdwJb2pmEeTLAToXnd13R8GFvftvqhrkySNwIRB\nn+TZSf7F3mXgF4C7gI3A6q7bauCqbnkjsCrJYUmWAsuAzcMuXJI0mEGmbhYAVybZ2/+TVXVtkpuB\nDUnOAR4AzgKoqq1JNgB3A08C51bVnmmpXpI0oQmDvqruA07YR/vXgVPH2WctsHbK1UmSpsxvxkpS\n4wx6SWqcQS9JjZvKN2NHbskF18zo8e6/6IwZPZ4kDYMjeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0\nktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9J\njTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4\ngYM+ybwktyW5uls/Osn1Se7tXo/q63thku1JtiU5bToKlyQNZjIj+vOAe/rWLwA2VdUyYFO3TpLl\nwCrgeOB04OIk84ZTriRpsgYK+iSLgDOAD/c1rwTWd8vrgTP72i+vqieqagewHTh5OOVKkiZr0BH9\nHwPnA9/ta1tQVTu75UeABd3yscCDff0e6tqeJsmaJFuSbNm9e/fkqpYkDWzCoE/yGmBXVd0yXp+q\nKqAmc+CqWldVK6pqxfz58yezqyRpEg4ZoM9PAa9N8kvA4cCRST4BPJpkYVXtTLIQ2NX1fxhY3Lf/\noq5NkjQCE47oq+rCqlpUVUvoXWT9XFX9KrARWN11Ww1c1S1vBFYlOSzJUmAZsHnolUuSBjLIiH48\nFwEbkpwDPACcBVBVW5NsAO4GngTOrao9U65UknRAJhX0VXUjcGO3/HXg1HH6rQXWTrE2SdIQ+M1Y\nSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJek\nxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqc\nQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY2bMOiTHJ5kc5I7kmxN8rtd\n+9FJrk9yb/d6VN8+FybZnmRbktOm8wQkSfs3yIj+CeBVVXUCcCJwepKXAxcAm6pqGbCpWyfJcmAV\ncDxwOnBxknnTUbwkaWITBn31/HO3emj3r4CVwPqufT1wZre8Eri8qp6oqh3AduDkoVYtSRrYQHP0\nSeYluR3YBVxfVTcBC6pqZ9flEWBBt3ws8GDf7g91bWPfc02SLUm27N69+4BPQJK0fwMFfVXtqaoT\ngUXAyUleOGZ70RvlD6yq1lXViqpaMX/+/MnsKkmahEnddVNV/xe4gd7c+6NJFgJ0r7u6bg8Di/t2\nW9S1SZJGYJC7buYneV63fATwauCrwEZgdddtNXBVt7wRWJXksCRLgWXA5mEXLkkazCED9FkIrO/u\nnPk+YENVXZ3ky8CGJOcADwBnAVTV1iQbgLuBJ4Fzq2rP9JQvSZrIhEFfVXcCJ+2j/evAqePssxZY\nO+XqJElT5jdjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJek\nxhn0ktQ4g16SGmfQS1LjBnkevUZkyQXXzOjx7r/ojBk9nqSZ4Yhekhpn0EtS4wx6SWqcQS9JjTPo\nJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuN8BIJGxkc8SDPDEb0kNc6gl6TGGfSS1DiDXpIaZ9BL\nUuMMeklqnEEvSY0z6CWpcRMGfZLFSW5IcneSrUnO69qPTnJ9knu716P69rkwyfYk25KcNp0nIEna\nv0FG9E8C/6GqlgMvB85Nshy4ANhUVcuATd063bZVwPHA6cDFSeZNR/GSpIlNGPRVtbOqbu2W/wm4\nBzgWWAms77qtB87sllcCl1fVE1W1A9gOnDzswiVJg5nUHH2SJcBJwE3Agqra2W16BFjQLR8LPNi3\n20NdmyRpBAYO+iTPAa4A3lNV3+jfVlUF1GQOnGRNki1JtuzevXsyu0qSJmGgoE9yKL2Q/+9V9emu\n+dEkC7vtC4FdXfvDwOK+3Rd1bU9TVeuqakVVrZg/f/6B1i9JmsAgd90EuAS4p6r+qG/TRmB1t7wa\nuKqvfVWSw5IsBZYBm4dXsiRpMgZ5Hv1PAW8BvpLk9q7tt4CLgA1JzgEeAM4CqKqtSTYAd9O7Y+fc\nqtoz9MolSQOZMOir6otAxtl86jj7rAXWTqEuSdKQ+M1YSWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS\n1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN\nM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiD\nXpIaZ9BLUuMMeklqnEEvSY0z6CWpcRMGfZKPJNmV5K6+tqOTXJ/k3u71qL5tFybZnmRbktOmq3BJ\n0mAGGdF/DDh9TNsFwKaqWgZs6tZJshxYBRzf7XNxknlDq1aSNGkTBn1VfR54bEzzSmB9t7weOLOv\n/fKqeqKqdgDbgZOHVKsk6QAc6Bz9gqra2S0/Aizolo8FHuzr91DX9gxJ1iTZkmTL7t27D7AMSdJE\npnwxtqoKqAPYb11VraiqFfPnz59qGZKkcRxo0D+aZCFA97qra38YWNzXb1HXJkkakQMN+o3A6m55\nNXBVX/uqJIclWQosAzZPrURJ0lQcMlGHJJcBrwSOSfIQ8NvARcCGJOcADwBnAVTV1iQbgLuBJ4Fz\nq2rPNNUuSRrAhEFfVWePs+nUcfqvBdZOpShJ0vD4zVhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLU\nOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z\n6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNO2TUBUitWnLBNTN6vPsvOmNGj6fZ\nwxG9JDXOoJekxjl1I2nSnJaaXQx6SRqjtQ8yp24kqXEGvSQ1zqCXpMZNW9AnOT3JtiTbk1wwXceR\nJO3ftAR9knnAB4BfBJYDZydZPh3HkiTt33SN6E8GtlfVfVX1beByYOU0HUuStB+pquG/afIG4PSq\n+jfd+luAl1XVO/v6rAHWdKs/CmwbeiHjOwb42gweb6Z5frNby+fX8rnBzJ/fD1bV/Ik6jew++qpa\nB6wbxbGTbKmqFaM49kzw/Ga3ls+v5XODg/f8pmvq5mFgcd/6oq5NkjTDpivobwaWJVma5FnAKmDj\nNB1LkrQf0zJ1U1VPJnkncB0wD/hIVW2djmMdoJFMGc0gz292a/n8Wj43OEjPb1ouxkqSDh5+M1aS\nGmfQS1LjDHpJapxBr1knySlJPjDqOqTZwh8emeWS/MT+tlfVrTNVy3RKchLwZuCNwA7g06OtSINK\nck5VXdK3Pg/4z1X1uyMsa6iSnAD8dLf6haq6Y5T1jDVngj7JZ4Gxtxj9I7AF+GBV/b+Zr2oo/lv3\nejiwArgDCPBieuf2ihHVNWVJXgCc3f37GvApeneK/dxIC5sGSV4OvB/4ceBZ9G5LfryqjhxpYcNx\napLXA+cARwMfA/5mpBUNUZLzgLfy1ODjE0nWVdX7R1jW08yZ2yuT/AkwH7isa3oT8A164X9kVb1l\nVLUNQ5JPA79dVV/p1l8I/E5VvWG0lR24JN8FvgCcU1Xbu7b7quqHRlvZ8CXZQu+LhX9B7wP7XwEv\nqKoLR1rYkCR5E70n2j4OvLmqvjTikoYmyZ3AK6rq8W792cCXq+rFo63sKXNmRA/8ZFW9tG/9s0lu\nrqqXJjmYvsx1oH50b8gDVNVdSX58lAUNwa/QC78bklxL7ymoGW1J06eqtieZV1V7gI8muQ2Y9UGf\nZBlwHnAFvb9Y3pLktqr65mgrG5oAe/rW93CQ/Xc6l4L+OUmOq6q/B0hyHPCcbtu3R1fW0NyZ5MPA\nJ7r1fwncOcJ6huHqqvpMN0JaCbwHeH6SPwOurKq/Gm15Q/XN7nEhtyd5L7CTdm6W+CxwblVtShLg\n39N7TMrxoy1raD4K3JTkym79TOCS/fSfcXNp6uaXgD8H/je9T9ulwDuAG4G3VtUfj666qUtyOPB2\nnrog9Hngz6rqidFVNTVJbq2qnxjTdhS9C7JvqqpTR1PZ8CX5QeBRevPz/w54LnDx3imr2SzJkVX1\njTFtL6iqvxtVTcPW3RRxSrf6haq6bZT1jDVngh4gyWHAj3Wr22bxBdjvSbISWFRVH+jWN9O7FlHA\n+VX1l6Osbyq6P+9PGnUd06n/r8zWJDm/qt7bLb+xqv6ib9sfVNVvja66qesGV/8W+BHgK8AlVfXk\naKvat7kW9D8JLKFvyqqqPj6ygoYgyZeAVVX1YLd+O/AqetNSH53No94kDwF/NN72qhp322zR/1dL\nkiuq6vWjrmlYxpzb0/4629dfa7NNkk8B36F3w8AvAvdX1XtGW9W+zZk5+iSXAj8M3M5TF04KmNVB\nDzxrb8h3vlhVjwGPdXPbs9k8eh9YB9WFrSHrP7fW7ibKOMv7Wp+NllfViwCSXAJsHnE945ozQU/v\nlrXl1d6fMEf1r/T/XCO9KZzZbGdV/d6oi5hmNc5yC/Z3bi2c63f2LnSPZh9lLfs1l4L+LuAH6N3N\n0JKbkry1qj7U35jkbRzEI4wBHbz/zxmeE5J8g965HtEt063XLP/C1P7O7fDRlTU0J4w5pyP6zveg\n+t9uzszRJ7kBOJFe+H3vTpSqeu3IihqCJM8HPkPvnPY+7uAlwGHAmVX16Khqm6okR3fTUJKmYC4F\n/c/uq72qmvgqdpJX8dR9yVur6nOjrEfSwWPOBL0kzVXNz9En+WJVnZLkn3j6BaCDbh5NkqaDI3pJ\nalwrz9KYUJIf7r4ZS5JXJnl3kueNui5Jmm5zJujpPTlvT5IfAdYBi4FPjrYkSZp+cynov9s9h+J1\nwPur6jeBhSOuSZKm3VwK+u8kORtYDVzdtR06wnokaUbMpaD/1/R+Vm9tVe1IshS4dMQ1SdK0m5N3\n3XTPNF9cVbP9hzkkaUJzZkSf5MYkRyY5mt6jAj6UZNY/5laSJjJngh54bvcrN78CfLyqXgb8/Ihr\nkqRpN5eC/pAkC4GzeOpirCQ1by4F/e8B1wHbq+rmJD8E3DvimiRp2s3Ji7GSNJfMhYeanV9V703y\nfvbxqzZV9e4RlCVJM6b5oAfu6V63jLQKSRoRp24kqXHNj+iTbNzf9tn+U4KSNJHmg57eYw8eBC4D\nbmJu/OC0JH1P81M3SeYBrwbOBl4MXANcVlVbR1qYJM2Q5u+jr6o9VXVtVa0GXg5sB25M8s4RlyZJ\nM2IuTN3Q/bLUGfRG9UuA9wFXjrImSZopc2Hq5uPAC4H/AVxeVXeNuCRJmlFzIei/CzzerfafbICq\nqiNnvipJmjnNB70kzXXNX4yVpLnOoJekxhn0ktQ4g14HhSR7ktze929JkhVJ3jfEY9yf5JhhvV/3\nnqck2Zzkq0m2JXnHAb7PkiSV5F19bX+a5Nf61g9JsjvJRUMoXXPInLiPXrPCt6rqxDFt97OPp44m\nOaSqnpyRqvYjyQ8AnwTOrKpbuw+R65LsrKoD+Z7GLuC8JB+sqm/vY/urgb8D3pjkwvJOCg3IEb0O\nWklemeTqbvl3klya5EvApUnmJfmvSW5OcmeSt/Xt8/kk13Qj7D9P8oz/zpN8JsktSbYmWdPXfnqS\nW5PckWRT1/bsJB/pRu63JVnZdT8X+FhV3QpQVV8Dzgd+s9vvY0ne0Pfe/zzBKe8GNgGrx9l+NvAn\nwN/Te4aTNBBH9DpYHJHk9m55R1W9bh99lgOnVNW3unD+x6p6affN5y8l+auu38ld3weAa+n9IPxf\njnmvX6+qx5IcAdyc5Ap6A58PAT9TVTuSHN31/U/A56rq15M8D9ic5K+B44H1Y953S3fsA/WHwP9M\n8pH+xiSH0/sx+7cBz6MX+v9rCsfRHGLQ62Cxr6mbsTZW1be65V8AXtw3Yn4usAz4NrC5qu4DSHIZ\ncArPDPp3J9n7YbK423c+8Pmq2gFQVY/1Heu1SX6jWz8cOG6yJziIqrovyU3Am8dseg1wQ/chdwXw\nX5K8p6r2TEcdaotBr9nk8b7lAO+qquv6OyR5Jc/8ycjaR5+fB15RVd9MciO98B5PgNdX1bYx73M3\n8BLgqr7ml/DUdYUn6aZHu+mjZ+3nGP3+gN4H09/0tZ0NnJLk/m79+4FXAdcP+J6aw5yj12x1HfD2\nJIcCJHlBkmd3205OsrQL1zcBXxyz73OBf+hC/sfoPdUU4G+Bn0mytHvPvVM31wHvSpKu/aSu/QPA\nryU5sWv/fmAt8Pvd9vvpBT/Aa4FDBzmxqvoqcDfwy937Hgn8NHBcVS2pqiX0rg+cPcj7SQa9ZqsP\n0wvDW5PcBXyQp/5CvRn4U3q/F7yDZz6p9FrgkCT3ABfRC3iqajewBvh0kjuAT3X9f59eSN+ZZGu3\nTlXtBH4VWJdkG/B/gPdV1d6R+IeAn+3e6xU8/S+SiawFFnXLr6N3jeCJvu1XAb/cXZ+Q9stn3agp\n3bTMb1TVa0Zw7HcAb6d3MfcfZvr40ngc0UtDUlUXV9WLDHkdbBzRSzMoyYuAS8c0P1FVLxtFPZob\nDHpJapxTN5LUOINekhpn0EtS4wx6SWrc/weJruW0ibHWyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8cb949a470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEtCAYAAAAGK6vfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHsRJREFUeJzt3Xu4XFV9//H3B6KIFxAeTiMSMGiDNngBG1OsPlZBBasV\nvGHwFhWNrXipl18btP3Zqqm0/VVrtaioaFSEBhWIaBEIUqtVQ0AUAiKpgYYIJN5RFAp8fn/sNTA5\nnJyZOZfsmcXn9Tx5zszae2a+k5x8Zs/aa68l20RERL12aruAiIiYXQn6iIjKJegjIiqXoI+IqFyC\nPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicnPaLgBgr7328vz589suIyJipFx88cU/tj3Wa7+hCPr5\n8+ezbt26tsuIiBgpkq7tZ7903UREVC5BHxFRuQR9RETlEvQREZVL0EdEVC5BHxFRuQR9RETlEvQR\nEZXrK+glPVDS5yR9X9KVkh4vaU9J50m6uvzco2v/4yVtkHSVpMNnr/yIiOil3ytj3w+cY/v5ku4N\n3Bd4G7DG9gmSlgPLgb+UtBBYAhwIPBg4X9IBtm+fqaLnL//STD3VhK454Zmz+vwRETtSzyN6SbsD\nTwI+DmD7Vts/B44EVpbdVgJHldtHAqfZvsX2RmADsHimC4+IiP7003WzP7AV+ISk70j6mKT7AXNt\nX1/2uQGYW27vA2zqevx1pS0iIlrQT9DPAR4LfMj2wcCvabpp7mTbgAd5YUnLJK2TtG7r1q2DPDQi\nIgbQT9BfB1xn+9vl/udogv9GSXsDlJ9byvbNwL5dj59X2rZh+yTbi2wvGhvrOctmRERMUc+gt30D\nsEnSw0vTYcAVwGpgaWlbCpxVbq8GlkjaRdL+wAJg7YxWHRERfet31M3rgVPKiJsfAq+g+ZBYJelY\n4FrgaADb6yWtovkwuA04biZH3ERExGD6CnrblwKLJth02Hb2XwGsmEZdERExQ3JlbERE5RL0ERGV\nS9BHRFQuQR8RUbkEfURE5RL0ERGVS9BHRFQuQR8RUbkEfURE5RL0ERGVS9BHRFQuQR8RUbkEfURE\n5RL0ERGVS9BHRFQuQR8RUbkEfURE5RL0ERGVS9BHRFQuQR8RUbkEfURE5RL0ERGVS9BHRFQuQR8R\nUbkEfURE5foKeknXSLpM0qWS1pW2PSWdJ+nq8nOPrv2Pl7RB0lWSDp+t4iMiordBjuifYvsg24vK\n/eXAGtsLgDXlPpIWAkuAA4EjgBMl7TyDNUdExACm03VzJLCy3F4JHNXVfprtW2xvBDYAi6fxOhER\nMQ39Br2B8yVdLGlZaZtr+/py+wZgbrm9D7Cp67HXlbaIiGjBnD73e6LtzZJ+BzhP0ve7N9q2JA/y\nwuUDYxnAfvvtN8hDIyJiAH0d0dveXH5uAc6g6Yq5UdLeAOXnlrL7ZmDfrofPK23jn/Mk24tsLxob\nG5v6O4iIiEn1DHpJ95P0gM5t4OnA5cBqYGnZbSlwVrm9GlgiaRdJ+wMLgLUzXXhERPSnn66bucAZ\nkjr7f9b2OZIuAlZJOha4FjgawPZ6SauAK4DbgONs3z4r1UdERE89g972D4HHTND+E+Cw7TxmBbBi\n2tVFRMS05crYiIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegj\nIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6\niIjKJegjIiqXoI+IqFyCPiKicgn6iIjK9R30knaW9B1JZ5f7e0o6T9LV5eceXfseL2mDpKskHT4b\nhUdERH8GOaJ/I3Bl1/3lwBrbC4A15T6SFgJLgAOBI4ATJe08M+VGRMSg+gp6SfOAZwIf62o+ElhZ\nbq8EjupqP832LbY3AhuAxTNTbkREDKrfI/p/Bv4CuKOrba7t68vtG4C55fY+wKau/a4rbRER0YKe\nQS/pWcAW2xdvbx/bBjzIC0taJmmdpHVbt24d5KERETGAfo7onwA8W9I1wGnAoZI+A9woaW+A8nNL\n2X8zsG/X4+eVtm3YPsn2ItuLxsbGpvEWIiJiMj2D3vbxtufZnk9zkvUC2y8BVgNLy25LgbPK7dXA\nEkm7SNofWACsnfHKIyKiL3Om8dgTgFWSjgWuBY4GsL1e0irgCuA24Djbt0+70oiImJKBgt72hcCF\n5fZPgMO2s98KYMU0a4uIiBmQK2MjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6\niIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyC\nPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIirXM+gl3UfSWknflbRe0t+W9j0lnSfp\n6vJzj67HHC9pg6SrJB0+m28gIiIm188R/S3AobYfAxwEHCHpEGA5sMb2AmBNuY+khcAS4EDgCOBE\nSTvPRvEREdFbz6B341fl7r3KHwNHAitL+0rgqHL7SOA027fY3ghsABbPaNUREdG3vvroJe0s6VJg\nC3Ce7W8Dc21fX3a5AZhbbu8DbOp6+HWlLSIiWtBX0Nu+3fZBwDxgsaRHjttumqP8vklaJmmdpHVb\nt24d5KERETGAgUbd2P458FWavvcbJe0NUH5uKbttBvbteti80jb+uU6yvcj2orGxsanUHhERfehn\n1M2YpAeW27sCTwO+D6wGlpbdlgJnldurgSWSdpG0P7AAWDvThUdERH/m9LHP3sDKMnJmJ2CV7bMl\nfRNYJelY4FrgaADb6yWtAq4AbgOOs3377JQfERG99Ax6298DDp6g/SfAYdt5zApgxbSri4iIacuV\nsRERlUvQR0RULkEfEVG5BH1EROUS9BERletneGXMsPnLvzSrz3/NCc+c1eePiNGSI/qIiMol6CMi\nKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXII+IqJyCfqI\niMol6CMiKpegj4ioXII+IqJyCfqIiMol6CMiKpegj4ioXM+gl7SvpK9KukLSeklvLO17SjpP0tXl\n5x5djzle0gZJV0k6fDbfQERETK6fI/rbgLfYXggcAhwnaSGwHFhjewGwptynbFsCHAgcAZwoaefZ\nKD4iInrrGfS2r7d9Sbl9E3AlsA9wJLCy7LYSOKrcPhI4zfYttjcCG4DFM114RET0Z6A+eknzgYOB\nbwNzbV9fNt0AzC239wE2dT3sutIWEREt6DvoJd0f+Dzw57Z/2b3NtgEP8sKSlklaJ2nd1q1bB3lo\nREQMoK+gl3QvmpA/xfYXSvONkvYu2/cGtpT2zcC+XQ+fV9q2Yfsk24tsLxobG5tq/RER0UM/o24E\nfBy40vZ7uzatBpaW20uBs7ral0jaRdL+wAJg7cyVHBERg5jTxz5PAF4KXCbp0tL2NuAEYJWkY4Fr\ngaMBbK+XtAq4gmbEznG2b5/xyiMioi89g9721wFtZ/Nh23nMCmDFNOqKiIgZkitjIyIql6CPiKhc\ngj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIq\nl6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionIJ+oiI\nyiXoIyIq1zPoJZ0saYuky7va9pR0nqSry889urYdL2mDpKskHT5bhUdERH/6OaL/JHDEuLblwBrb\nC4A15T6SFgJLgAPLY06UtPOMVRsREQPrGfS2vwb8dFzzkcDKcnslcFRX+2m2b7G9EdgALJ6hWiMi\nYgqm2kc/1/b15fYNwNxyex9gU9d+15W2iIhoybRPxto24EEfJ2mZpHWS1m3dunW6ZURExHZMNehv\nlLQ3QPm5pbRvBvbt2m9eabsb2yfZXmR70djY2BTLiIiIXqYa9KuBpeX2UuCsrvYlknaRtD+wAFg7\nvRIjImI65vTaQdKpwJOBvSRdB7wDOAFYJelY4FrgaADb6yWtAq4AbgOOs337LNUeERF96Bn0to/Z\nzqbDtrP/CmDFdIqKiIiZkytjIyIql6CPiKhcgj4ionIJ+oiIyiXoIyIql6CPiKhcgj4ionI9x9FH\njDd/+Zdm9fmvOeGZs/r8Efc0OaKPiKhcgj4ionLpuol7nHQ9xT1NjugjIiqXoI+IqFyCPiKicgn6\niIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqXoI+IqFyCPiKicgn6iIjKJegjIiqX2SsjRkxm\n34xBzdoRvaQjJF0laYOk5bP1OhERMblZCXpJOwP/CjwDWAgcI2nhbLxWRERMbra6bhYDG2z/EEDS\nacCRwBWz9HoRMQJGvdtpVOufra6bfYBNXfevK20REbGDyfbMP6n0fOAI268q918K/IHt13XtswxY\nVu4+HLhqxgu5y17Aj2fx+Wdb6m9X6m/PKNcOs1//Q2yP9dpptrpuNgP7dt2fV9ruZPsk4KRZev1t\nSFpne9GOeK3ZkPrblfrbM8q1w/DUP1tdNxcBCyTtL+newBJg9Sy9VkRETGJWjuht3ybpdcBXgJ2B\nk22vn43XioiIyc3aBVO2vwx8ebaef0A7pItoFqX+dqX+9oxy7TAk9c/KydiIiBgemesmIqJyCfqI\niMol6CMqJGmXftriniGzVw4ZSXtOtt32T3dULdMlaR/gIXT9ntn+WnsV3aN8E3hsH21DR9KuwG9t\nW9LDaC6oPNf2bS2X1hdJa2wf1qttR6oq6CXdBGz37LLt3XZgOVN1Mc17ELAf8LNy+4HA/wD7t1da\n/yT9PfBCmvmNbi/NBkYi6CV9kbv/Lv0CWAd8xPZvd3xVvUl6EM10I7tKOpjmdwdgN+C+rRU2mP8E\nniRpd+AC4BKaa3Fe1mpVPUi6D83f8V6S9mDbv/tWp4CpKuhtPwBA0ruA64FP0/xlvxjYu8XS+mZ7\nfwBJHwXOKMNUkfQM4Kg2axvQUcDDbd/SdiFT9ENgDDi13H8hcBNwAPBR4KUt1dXL4cDLaa5Gf29X\n+03A29ooaAp2sn2zpFcCH7J9gqRL2y6qD68B/hx4MM0BWyfofwl8sK2ioNLhlZK+a/sxvdqGmaTL\nbD+qV9uwkvTvwAts/6rtWqZC0kW2HzdRm6T1tg9sq7Z+SHqe7c+3XcdUlFB/NfAvwKttXz4qv/tl\niva32X5X27V0q+qIvsuvJb0YOI3m6/cxwK/bLWlgP5L0V8Bnyv0XAz9qsZ5B3QxcKmkNcOdRve03\ntFfSQO4vaT/b/wMgaT/g/mXbre2V1bezJb0ImM+250je2VpF/Xsz8LfA2SXkH0rTnTP0bN8u6blA\ngn4HeBHw/vLHwDdK2yg5BngHcEa5/7XSNipWM9rzG70F+Lqk/6b5Cr4/8FpJ9wNWtlpZf86iOadw\nMV0ftKPA9gXABZ1RQmVdi9e2W9VA1kh6HvAFD0mXSZVdNxEzoQTNI8rdq4b1BOxEJF1u+5Ft1zEV\nkhYDHwd2t72fpMcAr7L9+pZL60sZFHI/mkEIv6E5UHCbg0GqPKKXNEbTxzefbb+2vrKtmgYl6QDg\nrdz9PRzaVk39kHQZk498evQOLGe6fp+7/v4fIwnbn2q3pL79l6RH2b6s7UKm4F+AZwFnAtj+rqSn\ntFtS/zqDQoZJlUFP87X1P4HzuWto36g5Hfgw8DFG6z08q+0CZoKkTwMPAy5l2+GhoxL0TwReLmkj\nTddN56hyFD5od7J9raTutlH6P4CkZwNPKncvtH12m/XUGvT3tf2XbRcxTbfZ/lDbRQzK9rXQjKMf\n/29QxtaPyr/LImDhsPSxTsEz2i5gGjaV7huXUSyvB37Qck19k3QC8DjglNL0RklPsH18WzXVOgXC\n2ZL+uO0ipkLSnuXq2C9Keq2kvTttva6aHTJPm6BtlMLncuBBbRcxVeUDd1/g0HL7Zkbn//uf0Yy8\n2Q/YAhxS2kbFHwNPs32y7ZOBI4DZXbW8h6pOxnZdGSuakyG3AP/LEJwM6Vf5qt15D+PZ9kN3cEkD\nkfRnNCMkHgr8d9emBwDfsP2SVgobkKSvAgcBa9l2eOizWytqAJLeQfOt5OG2D5D0YOB0209oubTq\nSfoe8OTOdCXlAO3CNrvNquq6GcaTIIPqXBk7wj4L/DvwHmB5V/tNozRPD/A3bRcwTc8BDqaZPgDb\nP5I0Ev8/JM0H3gc8vjR9A3iL7WtaKmlQ7wG+Uw4WRNNXv3zyh8yuqo7oOyQ9B7jA9i/K/QfSfMKe\n2W5l/ZN0HHCK7Z+X+3sAx9g+sd3KJlfTpGyjTNJa24slXWL7sWX8/zdH4WSspG/SrMzU6eN+EfAa\n24/f/qOGi6S9afrpAdbavqHVeioN+kttHzSu7Tu2D26rpkGN6nvo6nqCu3c/jULX09dtP3GCCfJG\npvsPQNJbgQU050reA7wS+KztD7RaWB8kfW/8B9IoTGEi6Xdo5hP6XeAy4D22f9luVY1ag36iX5SR\nmCujo4xHf3Rn1EcZffC9YZ9jJYaHpKcBT6f5kPqK7fNaLqkvZdTKj7lrCpMXAnsBJwAMS3iOJ+kc\nmiuRv0YzzPgBtl/ealFFrUF/MvBz4F9L03HAnsPyl94PSf+PZtTBR0rTa4BNtt/SXlW9SXqE7e9L\nmnDec9uX7OiapqLMg36d7VskPRl4NPCpTlfaqJC0G9tecDf0XWeSNk2y2bb322HFDGD8t45Ot1mb\nNXXUGvT3A/4aeCrNEcF5wLtt39xqYQOQtBOwjOY9QPMePmr7jvaq6k3SSbaXlRNR43nYr+ztKDMo\nLqK5MvbLNBfhHWh7JIbtSnoNzcRgvwXu4K6up6HuOhtlkr4LPJm7uiy/2n2/zQ/ZWoP+BbZP79U2\nzCS90fb7e7XF7Og6ifl/aFY7+sAonCPpkHQ18HjbP267lkFJ+hZwMnCq7Zvarqdfkq7hrg/V8Vr9\nkK016O/2lWmYvkb1YzvvYZSCZsLVgEZlrhhJ3wb+GXg78Ce2N47SRGGlv/i5o/QttkPSI4BXAC8A\n/gv4hO017VbVHzXzNuzrMr31sKgq6NWswvTHwNHAv3Vt2o3mcvbFrRQ2AEnH0AwneyLbzsH9AOAO\nt7ju5CAkdY/uuA9wGHCJ7ee3VNJAJC0E/pRmSOKpkvYHjrb99y2X1hc1ywh+Avg2o7keQGcAwrNp\nVme6leYo/wPDfp5kGAd+VHXBFM3CHOtofjku7mq/CXhTKxUN7r9olkHcC/inrvabgO+1UtEUeNyU\nsuVahtNaKmdgtq8A3gB3XsPwgFEJ+eIjNOutXkbTnTBSygftK4A/oTk/cgrNwc8FDP8C55dIepzt\ni9oupKOqI/oOSX9h+x/GtY1c/7akhwALbJ8vaVdgzij1WXaTdC/gctsPb7uWfki6kOaAYQ7NQcMW\nmikc3txmXf0apW6+8SStpZmb52SaaRt+07Vt9bBPQyHp+zRj6a+lWdmu9ZlDaw36ke7fBpD0appR\nN3vafpikBcCHR6jr5ovcdcHRTsBCYJXtVi8F71fn90XSq2j6XN8x0fUZw0rS3wHXAF9k266boR1e\nKem5tr8g6QDbIzNb5XjlAO1uXGZ2bUNVQT9J//ZuNNP+PnXCBw6hMrxvMfDtzgfUMPb9bY+kP+q6\nextwre3r2qpnUOWCtafTLBv4dtsXjVjQb5ygeaiHV47agIleypWy9+ncb/MEbW199Nvr3+5cXTdK\nbrF9q8riC5LmMMnKTcPG9n90bkvaC/hJi+VMxTuBrwBfLyH/UODqlmvqm0d/cryRpWbRkX8CHkzT\n5fcQ4Eqgtavaqzqi71ZGHbyIZojWRuDztj/YblX9k/QPNFf3voxm4YXXAlfYfnurhfUg6RCaS9V/\nCrwL+DTNB+9OwMtsn9NiefcYkl4AnGP7Jkl/RXMC8122v9Nyadsl6WZgw0SbGJ3VsToXTh0KnF+6\n/54CvMT2sa3VVFPQq1ln9Zjy58c0QyzfanvCPrNhVq6MPZauuUqAj3nI/8EkraOZ2Gl3mhkIn2H7\nW2Vs9KnDfp6kcyK/DA+929/1qAxP7HQzSXoi8G7gH4H/a/sPWi5tuyStpxkePaE2+7gHIWmd7UUl\n8A+2fUfbk7LV1nXzfZq++WfZ3gAgaVSGVW6j/HKcCZxpe2vb9Qxgju1zASS90/a3AMr8N+1W1p8r\ny891rVYxfZ01Vp8JnGT7S5Le3WZBfbh1VMK8h59Luj/N5GanSNpCM/qmNbUF/XOBJcBXy5WBpzHx\n5chDq1xZ9w7gdZSl3yTdTnOhyDvbrK1P3WO2fzNu21B/GwGw/cXyc2XbtUzTZkkfoZmm+O8l7cLw\nLyX4jfLNbx+aQQi/6myQdMSwd/tJ+l1gLnAkze/+m4AX0/TRv36Sh866qrpuOsqkZkfSdOEcCnwK\nOKNzpDnMJL2ZZm3VZbY3lraHAh+i6XN9X5v19VI+lDpjh3elGQ9NuX8f2/dqq7Z+SFo92fZhH8Pd\nIem+NGuVXmb7ajULYTxqmP8PSHoDzUyzV9Is4/hG22eVbUM/IkfS2cDxti8b1/4o4O9s/0k7lVUa\n9N3KVY0vAF44CmPQJX2HZmHhH49rHwPOHfY+7lEnaSuwCTiVZvqAbb4Rdo8mGgXDNMSvlzKk9fG2\nf6VmOcHPAZ+2/f5RuA5G0kW2H7edba0Oja6t6+ZubP+M5qTgSW3X0qd7jQ95ANtby9WlMbseRNPd\n0bkm40s0J5HXt1rVgCYY4rcfzTmsYV64ZqdOd43ta9SsA/C5cgHSKHTBPnCSbbvusComMOx9dvdE\nt05xW8wA27fbPsf2UuAQmuF+F0p6XculDepdNPX/oIypfyrwrXZL6ulGSXcun1lC/1k0w3NH4ULB\ndeWK9m2Uq6svnmD/Hab6rptR09XHfbdNjEAfdw3Kictn0hzVzwdWAyfb3txmXYMYxiF+vUiaR3MF\n+90W0pb0BNvfaKGsvkmaC5xBc0DWCfZFwL2B50z0vnZYbQn6iLtI+hTwSJpVpU6zfXnLJU2JpPOB\no2gWBt+Lpvvmcbb/sNXC7gHKBVKddQvW276gzXogQR+xDUl3cNc3qu7/HJ2rM3fb8VUNrow8+w1N\n9+yLaS5gO8X2qE1FETMgQR9Ruc5cQ8N+VXXMnpyMjaiIpEMkXSjpC5IOlnQ5cDnNic4j2q4v2pEj\n+oiKjPpcQzE7ckQfUZc5ts+1fTpwQ/dcQy3XFS1K0EfUZaTnGorZka6biIqM+lxDMTsS9BERlUvX\nTURE5RL0ERGVS9BHRFQuQR9DQdJcSZ+V9ENJF0v6pqTntFDHKyRdWv7cKumycvuEWX7d35VkSX/W\n1fZhSS/pun9vST8dgSUBY8gk6KN1ZfnEM4Gv2X6o7d+nWRJyXp+Pn7F1FWx/wvZBtg8CfgQ8pdxf\nPlOvMYkbgTdN8n4OB64AXrgDaomKJOhjGBxKszD0hzsNtq+1/QFJ8yX9p6RLyp8/BJD05NK+mib8\nkHRm+TawXtKyznNJOlbSDyStlfRRSR8s7WOSPi/povLnCdsrUNJOkjZI2rPc37l8+9hT0mckfai8\n9g8kPaPsM0fSe8vrfq/MSz6ZG2gWt3/pdrYfA7wXuEHS4h7PFXGn6leYipFwIHDJdrZtoVla8beS\nFtAs8beobHss8MjO2rrAK23/VNKuwEWSPg/sAvx12fcm4ALgu2X/9wPvs/11SfsBXwF+b6Iiynzu\np9KsOvVBmqPri8rrAewLPA5YAJxfFoo+Fthie3GZ4/5bks7tsZzfCcBZkrZZnLysAftk4JU0q2Ad\nA6yd5Hki7pSgj6Ej6V+BJ9Is4PBU4INl5aHbgQO6dl3bFfIAb+jq19+XJnQfBPyH7Z+W5z696zme\nCiwsQQ2wm6T7d5azm8DHgdNpgv6VwMe6tq2yfQdwlaRN5bWfDvyepCVln91L+3aDvizkfSl37555\nNnBe+cA7HbhY0lvKa0ZMKkEfw2A98LzOHdvHlal11wFvoum7fgxNV+Nvux5350pcZX3Rp9IsLn2z\npAvpWhR7O3YCDrH92x77deq6RtLPysISBwPndm8evzvN1aivtb2mn+fvsgL4LNsu/XcMcIika8r9\nMeCPgK8O+NxxD5Q++hgGFwD36R5xAty3/NwduL4cub4U2Hk7z7E78LMS8o+gWS8V4CLgjyTtUU5y\nPq/rMecCr+/c6V6vdBIfB06hWX2q+2j6BWocQPNt4mqarqDXdk6uSnp46VaaVFmI/L+BTl//A8v7\nmWd7vu35wBtowj+ipwR9tK4siHEUTSBvlLQWWAn8JXAisFTN2qePYOL1dAHOAeZIupKmn7sza+Nm\n4O9o+rO/AVwD/KI85g3AonKi9ArgT/so9wyaD5VPjmvfTPMN5IvAMtu3Ah+hCfxLy7zwH6L/b9Hv\npvnAgObD6Tzb/9u1/UzgKEmZuyZ6ylw3Ub1Ov3s5sj6DZqHvM6b4XIcA77H9lK62zwCfs33mzFQc\nMbNyRB/3BH9TTnBeDmykORoemKS3A/9Gs7BHxMjIEX3EDlTOA3xyXPPNtv+whXLiHiJBHxFRuXTd\nRERULkEfEVG5BH1EROUS9BERlUvQR0RU7v8Dz2thssqBmaAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8cba8b32e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# let's look at the number of observations per label on each of the variables\n",
    "for col in ['BsmtQual_NA', 'FireplaceQu_NA', 'GarageType_NA']:\n",
    "    X_train.groupby([col])[col].count().sort_values(ascending=False).plot.bar()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe in the plots the label indicating the number of observations with missing values, within each of the 3 categorical variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's transform the categories into numbers quick and dirty so we can use them in scikit-learn\n",
    "\n",
    "# the below function numbers the labels from 0 to n, n being the number of different labels \n",
    "#  within the variable\n",
    "\n",
    "for col in ['BsmtQual_NA', 'FireplaceQu_NA', 'GarageType_NA',]:\n",
    "    labels_dict = {k:i for i, k in enumerate(X_train[col].unique(), 0)}\n",
    "    X_train.loc[:, col] = X_train.loc[:, col].map(labels_dict )\n",
    "    X_test.loc[:, col] = X_test.loc[:, col].map(labels_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train set random imputation\n",
      "Linear Regression mse: 4810016310.466396\n",
      "Test set random imputation\n",
      "Linear Regression mse: 5562566516.826057\n"
     ]
    }
   ],
   "source": [
    "# Let's evaluate the performance of Linear Regression\n",
    "\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(X_train[['BsmtQual_NA', 'FireplaceQu_NA', 'GarageType_NA',]], y_train)\n",
    "print('Train set random imputation')\n",
    "pred = linreg.predict(X_train[['BsmtQual_NA', 'FireplaceQu_NA', 'GarageType_NA',]])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_train, pred)))\n",
    "print('Test set random imputation')\n",
    "pred = linreg.predict(X_test[['BsmtQual_NA', 'FireplaceQu_NA', 'GarageType_NA',]])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_test, pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous lectures we trained linear regressions on data where missing observations were replaced by i) random sampling, or ii) random sampling plus a variable to indicate missingness. We obtained the following mse for the testing sets:\n",
    "\n",
    "- frequent label imputation mse: 6456070592\n",
    "- random sampling + additional category: 4911877327\n",
    "- adding 'missing' label: 5562566516\n",
    "\n",
    "Therefore, adding an additional 'Missing' category lies between the 2 other approaches in terms of performance.\n",
    "\n",
    "A next step could be to investigate which approach works better for each variable individually, to try and optimise the performance of logistic regression even more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**That is all for this demonstration. I hope you enjoyed the notebook, and see you in the next one.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}