05.2_Mean_and_median_imputation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean and median imputation\n",
    "\n",
    "Imputation is the act of replacing missing data with statistical estimates of the missing values. The goal of any imputation technique is to produce a **complete dataset** that can then be then used for machine learning.\n",
    "\n",
    "Mean/median imputation consists of replacing all occurrences of missing values (NA) within a variable by the mean (if the variable has a Gaussian distribution) or median (if the variable has a skewed distribution).\n",
    "\n",
    "\n",
    "### Assumptions\n",
    "\n",
    "Mean/median imputation has the assumption that the data are missing completely at random (MCAR). If this is the case, we can think of replacing the NA with the  most frequent occurrence of the variable, which is the mean if the variable has a Gaussian distribution, or the median otherwise.\n",
    "\n",
    "The rationale is to replace the population of missing values with the most frequent value, since this is the most likely occurrence.\n",
    "\n",
    "### Advantages\n",
    "\n",
    "- Easy to implement\n",
    "- Fast way of obtaining complete datasets\n",
    "\n",
    "### Limitations\n",
    "\n",
    "- Distortion of original variance\n",
    "- Distortion of covariance with remaining variables within the dataset\n",
    "\n",
    "When replacing NA with the mean or median, the variance of the variable will be distorted if the number of NA is big respect to the total number of observations (since the imputed values do not differ from the mean or from each other). Therefore leading to underestimation of the variance.\n",
    "\n",
    "In addition, estimates of covariance and correlations with other variables in the dataset may also be affected.  This is because we may be destroying intrinsic correlations since the mean/median that now replace NA will not preserve the relation with the remaining variables.\n",
    "\n",
    "\n",
    "### Final note\n",
    "Replacement of NA with mean/median is widely used in the data science community and in various data science competitions. 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",
    "\n",
    "Typically, mean/median imputation is done together with adding a variable to capture those observations where the data was missing (see lecture \"Creating a new variable with the missing data\"), thus covering 2 angles: if the data was missing completely at random, this would be contemplated by the mean imputation, and if it wasn't this would be captured by the additional variable.\n",
    "\n",
    "In addition, both methods are extremely straight forward to implement, and therefore are a top choice in data science competitions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "===============================================================================\n",
    "\n",
    "## Real Life example: \n",
    "\n",
    "### Predicting Survival on the Titanic: understanding society behaviour and beliefs\n",
    "\n",
    "Perhaps one of the most infamous shipwrecks in history, the Titanic sank after colliding with an iceberg, killing 1502 out of 2224 people on board. Interestingly, by analysing the probability of survival based on few attributes like gender, age, and social status, we can make very accurate predictions on which passengers would survive. Some groups of people were more likely to survive than others, such as women, children, and the upper-class. Therefore, we can learn about the society priorities and privileges at the time.\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 mean/median imputation using the Titanic and 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\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "# for classification\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.svm import SVC\n",
    "\n",
    "# to split and standarize the datasets\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# to evaluate regression models\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "# to evaluate classification models\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Titanic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>Survived</th>\n",
       "      <th>Age</th>\n",
       "      <th>Fare</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>7.2500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>38.0</td>\n",
       "      <td>71.2833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>26.0</td>\n",
       "      <td>7.9250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>35.0</td>\n",
       "      <td>53.1000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>8.0500</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Survived   Age     Fare\n",
       "0         0  22.0   7.2500\n",
       "1         1  38.0  71.2833\n",
       "2         1  26.0   7.9250\n",
       "3         1  35.0  53.1000\n",
       "4         0  35.0   8.0500"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load the Titanic Dataset with a few variables for demonstration\n",
    "\n",
    "data = pd.read_csv('titanic.csv', usecols = ['Age', 'Fare','Survived'])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Survived    0.000000\n",
       "Age         0.198653\n",
       "Fare        0.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's look at the percentage of NA\n",
    "\n",
    "data.isnull().mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imputation important\n",
    "\n",
    "Imputation should be done over the training set, and then propagated to the test set. This means that the mean/median to be used to fill missing values both in train and test set, should be extracted from the train set only. And this is to avoid overfitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((623, 3), (268, 3))"
      ]
     },
     "execution_count": 4,
     "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, data.Survived, test_size=0.3,\n",
    "                                                    random_state=0)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's make a function to create 2 variables from Age:\n",
    "# one filling NA with median, and another one filling NA with zeroes\n",
    "\n",
    "def impute_na(df, variable, median):\n",
    "    df[variable+'_median'] = df[variable].fillna(median)\n",
    "    df[variable+'_zero'] = df[variable].fillna(0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "median = X_train.Age.median()\n",
    "median"
   ]
  },
  {
   "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>Survived</th>\n",
       "      <th>Age</th>\n",
       "      <th>Fare</th>\n",
       "      <th>Age_median</th>\n",
       "      <th>Age_zero</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>857</th>\n",
       "      <td>1</td>\n",
       "      <td>51.0</td>\n",
       "      <td>26.5500</td>\n",
       "      <td>51.0</td>\n",
       "      <td>51.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>1</td>\n",
       "      <td>49.0</td>\n",
       "      <td>76.7292</td>\n",
       "      <td>49.0</td>\n",
       "      <td>49.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>386</th>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>46.9000</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>124</th>\n",
       "      <td>0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>77.2875</td>\n",
       "      <td>54.0</td>\n",
       "      <td>54.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>578</th>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14.4583</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>549</th>\n",
       "      <td>1</td>\n",
       "      <td>8.0</td>\n",
       "      <td>36.7500</td>\n",
       "      <td>8.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>118</th>\n",
       "      <td>0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>247.5208</td>\n",
       "      <td>24.0</td>\n",
       "      <td>24.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>8.0500</td>\n",
       "      <td>20.0</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>157</th>\n",
       "      <td>0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>8.0500</td>\n",
       "      <td>30.0</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>127</th>\n",
       "      <td>1</td>\n",
       "      <td>24.0</td>\n",
       "      <td>7.1417</td>\n",
       "      <td>24.0</td>\n",
       "      <td>24.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>653</th>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7.8292</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>235</th>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7.5500</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>785</th>\n",
       "      <td>0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>7.2500</td>\n",
       "      <td>25.0</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>241</th>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>15.5000</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>351</th>\n",
       "      <td>0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>35.0000</td>\n",
       "      <td>29.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     Survived   Age      Fare  Age_median  Age_zero\n",
       "857         1  51.0   26.5500        51.0      51.0\n",
       "52          1  49.0   76.7292        49.0      49.0\n",
       "386         0   1.0   46.9000         1.0       1.0\n",
       "124         0  54.0   77.2875        54.0      54.0\n",
       "578         0   NaN   14.4583        29.0       0.0\n",
       "549         1   8.0   36.7500         8.0       8.0\n",
       "118         0  24.0  247.5208        24.0      24.0\n",
       "12          0  20.0    8.0500        20.0      20.0\n",
       "157         0  30.0    8.0500        30.0      30.0\n",
       "127         1  24.0    7.1417        24.0      24.0\n",
       "653         1   NaN    7.8292        29.0       0.0\n",
       "235         0   NaN    7.5500        29.0       0.0\n",
       "785         0  25.0    7.2500        25.0      25.0\n",
       "241         1   NaN   15.5000        29.0       0.0\n",
       "351         0   NaN   35.0000        29.0       0.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "impute_na(X_train, 'Age', median)\n",
    "X_train.head(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "impute_na(X_test, 'Age', median)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean/median imputation alters the variance of the original distribution of the variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original Variance:  14.506993804588443\n",
      "Variance after median imputation:  13.024746215318798\n"
     ]
    }
   ],
   "source": [
    "# we can see a change in the variance after imputation\n",
    "\n",
    "print('Original Variance: ', X_train['Age'].std())\n",
    "print('Variance after median imputation: ', X_train['Age_median'].std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xc46fbd5f98>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOXZ//HPlZ0QQsgCJCQQdglC2BfRilIUEcUqVfRx\nb7XUtbU+Vmtt9Vef2tpVW5e6IkpBQQTcqlZFBdkS9oSdBBKSsCQEAmTP/fvjnsQQQphAZs4kc71f\nr7ySuc+ZmYsl+eacexNjDEoppdTpBDhdgFJKqdZBA0MppZRbNDCUUkq5RQNDKaWUWzQwlFJKuUUD\nQymllFs0MJRSSrlFA0MppZRbNDCUUkq5JcjpAlpSbGysSU5OdroMpZRqNdLT0w8aY+LcObdNBUZy\ncjJpaWlOl6GUUq2GiOx291y9JaWUUsotGhhKKaXcooGhlFLKLW2qD0Mp5bsqKyvJzc2lrKzM6VL8\nUlhYGImJiQQHB5/xa2hgKKW8Ijc3lw4dOpCcnIyIOF2OXzHGUFhYSG5uLj179jzj19FbUkoprygr\nKyMmJkbDwgEiQkxMzFlf3WlgKKW8RsPCOS3xd6+BoRRAWRn885+wZInTlSjlszQwlAK491778f3v\nw5o1TlejPGjhwoWICFu2bHG6lFZHA0OpvXth5ky45RaIiICnnnK6IuVBc+bM4fzzz2fOnDlOl9Lq\naGAoNXs2VFXBb34Dt90GixfD0aNOV6U84OjRoyxdupRXX32VuXPnAlBTU8Ndd93FOeecw8SJE5k8\neTLz588HID09nQsvvJDhw4dz6aWXkp+f72T5jtNhtUp9/DEMHgy9esEVV8Df/w6ffw5TpzpdWZv1\nxPsZZOYdadHXTEmI5LdXDGzynEWLFjFp0iT69etHTEwM6enpZGVlkZ2dTWZmJvv372fAgAHcfvvt\nVFZWcu+997Jo0SLi4uJ4++23efTRR3nttddatO7WRAND+bfSUli6FH7+c/v4/POhXTvb+a2B0ebM\nmTOH+++/H4Dp06czZ84cqqqq+OEPf0hAQABdu3bloosuAmDr1q1s2rSJiRMnAlBdXU18fLxjtfsC\nDQzl39autbejzjvPPg4JgeHDYcUKZ+tq4053JeAJRUVFfPHFF2zcuBERobq6GhHhBz/4QaPnG2MY\nOHAgy5cv93Klvkv7MJR/W73afh458ru2MWPsSKnycmdqUh4xf/58brrpJnbv3k12djY5OTn07NmT\n6Oho3n33XWpqati3bx9LXEOr+/fvz4EDB+oCo7KykoyMDAf/BM7TwFD+bfVqSEiAbt2+axs2DCoq\nYNs25+pSLW7OnDknXU1cc801FBQUkJiYSEpKCjfeeCPDhg2jY8eOhISEMH/+fH75y1+SmprKkCFD\n+Pbbbx2q3jfoLSnl39LSTry6AEhJsZ8zMmDQIO/XpDziyy+/PKntvvvuA+zoqYiICAoLCxk1ahSD\nXP/uQ4YM4euvv/Zqnb5MA0P5r7Iy2L4drrvuxPb+/SEgwAaG8gtTpkyhuLiYiooKHnvsMbp27ep0\nST5JA0P5r+3boaYGBgw4sT0sDPr0gcxMZ+pSXrdEl4Rxi/ZhKP9VuzTEOeecfGzgQL3CUKoBDQzl\nv7ZsARHo1+/kYwMGwI4dUFnp/bqU8lEaGMp/bd4MPXpAePjJx3r3hupqyMnxfl1K+SgNDOW/tmxp\n/HYU2GVCAHbt8l49Svk4DQzln2pqYOvWUwdG7TaWWVneq0kpH6eBofzT3r1w/PipAyMxEYKC9Aqj\nDWoN+2HMnDmTe+65B4AXX3yRWbNmOVyRpYGh/FPtlUPtraeGAgNt/4YGRpvT2vbDmDFjBjfffLPT\nZQA6D0P5q+xs+zk5+dTn9Oqlt6Q85Wc/g3XrWvY1hwyxS9M3oXY/jC+//JIrrriCJ554gpqaGu65\n5x6++OILkpKSCA4O5vbbb2fatGmkp6fzwAMPcPToUWJjY5k5c+YpV6wdP348Q4cO5ZtvvuHYsWPM\nmjWLp556io0bN3Ldddfx5JNPAvDWW2/x7LPPUlFRwejRo3n++ecJDAzk9ddf56mnniIqKorU1FRC\nQ0MBePzxx4mIiODBBx/k5Zdf5qWXXqKiooI+ffrw5ptvEh4ezq233kpkZCRpaWkUFBTw9NNPM23a\ntJb9+0WvMJS/qg2MpKRTn9Ozp15htDGN7YexYMGCuv0w3nzzzRMWG7z33nuZP38+6enp3H777Tz6\n6KNNvn5ISAhpaWnMmDGDqVOn8txzz7Fp0yZmzpxJYWEhmzdv5u2332bZsmWsW7eOwMBAZs+eTX5+\nPr/97W9ZtmwZS5cuJfMUk0avvvpqVq9ezfr16xkwYACvvvpq3bH8/HyWLl3KBx98wMMPP9xyf2n1\nePQKQ0QmAc8AgcArxpg/NDguruOTgePArcaYNfWOBwJpwF5jzBRP1qr8zO7dEB9vZ3WfSq9eUFgI\nR45AZKT3avMHp7kS8BRP74dx5ZVXAjBo0CAGDhxYd36vXr3Iyclh6dKlpKenM9K1fllpaSmdO3dm\n5cqVjB8/nri4OACuu+46tjWy+OWmTZv49a9/TXFxMUePHuXSSy+tO3bVVVcREBBASkoK+/btO5u/\nplPyWGC4ftg/B0wEcoHVIrLYGFM/Oi8D+ro+RgMvuD7Xuh/YDOh3q2pZ2dlN346CE0dKpaZ6uiLl\nYd7YD6P2NlJAQEDd17WPq6qqMMZwyy238FSDfeMXLlzo1uvfeuutLFy4kNTUVGbOnHnCkib1388Y\n43bNzeHJW1KjgB3GmF3GmApgLtBwC7OpwCxjrQCiRCQeQEQSgcuBVzxYo/JXu3fbTu2mdO9uP+vk\nvTbBF/bDmDBhAvPnz2f//v2ADbHdu3czevRovvrqKwoLC6msrGTevHmNPr+kpIT4+HgqKyuZPXv2\nWdVyJjwZGN2A+t9pua42d8/5O/AQUNPUm4jInSKSJiJpBw4cOLuKlX+oroY9e05/hVHbv5Gb6/GS\nlOf5wn4YKSkpPPnkk1xyySUMHjyYiRMnkp+fT3x8PI8//jhjx45l3LhxDGi4IKbL7373O0aPHs24\nceM451RDwj3JGOORD2Aatt+i9vFNwD8bnPMBcH69x58DI4ApwPOutvHAB+685/Dhw41Sp5WTYwwY\n88ILTZ9XVWVMYKAxjz7qnbrauMzMTKdLOKWSkhJjjDEHDx40vXr1Mvn5+Q5X5BmN/RsAacbNn+ue\n7PTeC9QfgpLoanPnnGuAK0VkMhAGRIrIW8aYGz1Yr/IXu3fbz6e7wggMtB3jeoXR5ul+GO7xZGCs\nBvqKSE9sCEwHbmhwzmLgHhGZi+3sPmyMyQcecX0gIuOBBzUsVItxZw5GrcREDQw/0Jz9MO6++26W\nLVt2Qtv999/Pbbfd1sJV+R6PBYYxpkpE7gE+wQ6rfc0YkyEiM1zHXwQ+wg6p3YEdVtv2/8aV82qv\nMGo7tZuSmAgbN3q2Hj9ijMGOpm+9nnvuOadLOCOmBUZOeXQehjHmI2wo1G97sd7XBrj7NK+xBFji\ngfKUv8rOhs6dG1/WvKHERPj4YzDG7p2hzlhYWBiFhYXExMS0+tBobYwxFBYWEtbUvCM36NIgyv+4\nMwejVmIiHDsGhw9DVJQnq2rzEhMTyc3NRUczOiMsLIzExMSzeg0NDOV/du92fyJe7TdYbq4GxlkK\nDg6mZ+1kSNUq6VpSyr/U1NjAaM4VBmjHt1JoYCh/s28flJeffpZ3LZ28p1QdDQzlX9ydg1ErPt52\ndmtgKKWBofxMc+ZgAAQHQ9euGhhKoYGh/E3tFYa7t6TA9mPoAoRKaWAoP5OdDTExEBHh/nO6dbN7\ngCvl5zQwlH/Jzm7e1QVoYCjlooGh/EtzhtTWSkiA4mI4ftwjJSnVWmhgKP9hTPNmedfq5tqiJS+v\npStSqlXRwFD+4+BBKC09s1tSoLellN/TwFD+o7lDamtpYCgFaGAof1IbGHqFodQZ0cBQ/uNM5mAA\nREbaYbgaGMrPaWAo/5GdbVecPZNVZ3VorVIaGMqPnMkcjFoJCRoYyu9pYCj/cSZzMGp166bDapXf\n08BQ/qF2DsaZXmHUBkZNTYuWpVRrooGh/MOhQ3D06NkFRmWlncuhlJ/SwFD+ISvLfj7TLUJ1aK1S\nGhjKT2hgKHXWNDCUf9DAUOqsaWAo/5CdDZ06QceOZ/b8rl3tVq06Ukr5MQ0M5R+yss58SC1AUBB0\n6aJXGMqvaWAo/5CVdea3o2rpbG/l5zQwVNtXOwdDA0Ops6KBodq+ffugrEwDQ6mzpIGh2r7aEVJn\n04cBNjCKiuwmTEr5IQ0M1fad7ZDaWrpVq/JzGhiq7TvTnfYa0sBQfk4DQ7V9WVl2SGx4+Nm9jk7e\nU35OA0O1fWc7B6NWQoL9rIGh/JQGhmr7WmIOBthZ4uHhGhjKb3k0MERkkohsFZEdIvJwI8dFRJ51\nHd8gIsNc7WEiskpE1otIhog84ck6VRtWWWk3TurV6+xfS0SH1iq/5rHAEJFA4DngMiAFuF5EUhqc\ndhnQ1/VxJ/CCq70cuNgYkwoMASaJyBhP1arasOxsqK6Gfv1a5vU0MJQf8+QVxihghzFmlzGmApgL\nTG1wzlRglrFWAFEiEu96fNR1TrDrw3iwVtVWbd9uP/ft2zKvp4Gh/JgnA6MbkFPvca6rza1zRCRQ\nRNYB+4HPjDErPViraqtOERjGnOHvH7VbtZ7p85VqxYKcLuBUjDHVwBARiQLeE5FzjTGbGp4nIndi\nb2fRvXt3L1epfN727bazOjaWvOJSXliyk88y97GvpIzo8BDO6xPLj8/vSWpSlHuv160bVFRAYSHE\nxnq2dqV8jCcDYy+QVO9xoqutWecYY4pF5EtgEnBSYBhjXgJeAhgxYoT+2qdOtG0b9O3LfzIKeOCd\n9VTVGCYO6EKvuPbkHy7j04wC3l+fx81je/CryQMICw5s+vXqD63VwFB+xpOBsRroKyI9sSEwHbih\nwTmLgXtEZC4wGjhsjMkXkTig0hUW7YCJwB89WKtqq7ZvJ2/AUH46ew1Dk6J4ZvpQkqK/m8BXckUK\nf/tsO68ty2JLfgmv3DqCyLDgU79e/cl7qakeLl4p3+KxPgxjTBVwD/AJsBl4xxiTISIzRGSG67SP\ngF3ADuBl4C5XezzwpYhswAbPZ8aYDzxVq2qjyssxe/awoKQdQ5KimP3jMSeEBUCHsGB+c0UK/7h+\nKGtzDnHra6soq6w+9WvqbG/lxzzah2GM+QgbCvXbXqz3tQHubuR5G4ChnqxNtX0V23YQUlPDvi5J\nvHLzCNqFnPp20xWpCQSIcPe/1/DgvPX84/qhiMjJJ8bH2/kYGhjKD+lMb9VmffH+UgCmXPM9YiJC\nT3v+5YPjeWhSfz7YkM8b32Y3flJwMHTurIGh/JIGhmqTjpRVkvlVGgCjL3V/zudPL+zNRf3j+P3H\nW9i2r6Txk2qH1irlZzQwVJv06jdZJOZnUxnXGaKj3X6eiPD0tFQiw4J4cN56qmsaGXink/eUn9LA\nUG1OWWU1b67YzYhj+QSfO7DZz4/rEMpjU1LYkHuYuav3nHyCBobyUxoYqs35aGM+RUfL6b4vGwY2\nPzAArkxNYEyvaJ7+z1aKjlWceDAhAQ4ehPLysy9WqVZEA0O1ObOW72ZM8DGCjh0748AQEX439VyO\nlVfxt8+2nXhQd95TfkoDQ7UpWQePsS6nmFsij9mGlIYLJLuvb5cOTB+VxJxVe9hTePy7AzoXQ/kp\nDQzVpry/3v7WP65iv204wyuMWvde3JegQOFv/613laGBofyUBoZqM4wxLF6fx6jkaCJ3bbP7eMfE\nnNVrdokM45bzklm4bi9bC1zDbPWWlPJTGhiqzdhSUMKO/Ue5YkgCZGSc1e2o+n56YW/CgwN5YckO\n2xAVBe3a6RWG8jtuBYaILBCRy0VEA0b5rE8yChCBy1I6Q2bmWd+OqhUVHsINo7vz/oZ8coqO26VB\nkpJgTyNDbpVqw9wNgOexK81uF5E/iEh/D9ak1Bn5cusBhiRFEZu3G44ehWHDWuy1f3R+LwIEXv5m\nl23o2ROyslrs9ZVqDdwKDGPMf40x/wMMA7KB/4rItyJym4g0sRa0Ut5ReLScDbnFjO/XGVatso0j\nR7bY63ftGMbVQxN5e3UOB0rKNTCUX3L7FpOIxAC3Aj8G1gLPYAPkM49UplQzfLP9IMbA+P5xsHo1\ntG8PAwa06Hv85MJeVFTX8ObybEhOtrvulZxivSml2iB3+zDeA74BwoErjDFXGmPeNsbcC0R4skCl\n3LFk635i2ocwqFtHGxjDh0PgaXbPa6ZecRFc3L8z/161h8ruPWyjXmUoP+LuFcbLxpgUY8xTxph8\nABEJBTDGjPBYdUq5oabG8M32g3yvXxwBVZWwdm2L3o6q7+bzkjl4tIJl1R1sgwaG8iPuBsaTjbQt\nb8lClDpT2/cfpfBYBef1joGNG6GiAkaN8sh7XdAnll6x7XmtdgqGBobyI00Ghoh0FZHhQDsRGSoi\nw1wf47G3p5Ry3MqsQgDG9IrxSId3fQEBwk1je/B1oaG6fYQGhvIrp9ui9VJsR3ci8Nd67SXArzxU\nk1LNsnJXEQkdw0js1M72X8TG2k5pD7lmeCJ/+mQr+2LiSdDAUH6kycAwxrwBvCEi1xhj3vVSTUq5\nzRjDyqxCLugbZ/fgXr3aXl00th93C4kMC+bqYd3InBtDl527aNmudaV8V5OBISI3GmPeApJF5IGG\nx40xf23kaUp5zc4Dxzh4tIJRPaPtZL3MTLj6ao+/7/SR3VkV2ZmazA0EGuPRgFLKV5yu07u963ME\n0KGRD6UctSqrCIDRPaNhzRqoqfFYh3d9AxMiqejeg+DS43YzJaX8wOluSf3L9fkJ75SjVPOs2XOI\nmPYh9Ixtb29Hgcc6vOsTEXqPTYX3IGv5enpe+X2Pv6dSTnN34t7TIhIpIsEi8rmIHBCRGz1dnFKn\nsy6nmCFJUbb/YtUq6NEDOnf2ynuPumQMAOu/WOmV91PKae7Ow7jEGHMEmIJdS6oP8L+eKkopdxwu\nrWTH/qMMSYqyDbUd3l7ScWB/qgKDKErfSHlVtdfeVymnuBsYtbeuLgfmGWMOe6gepdy2IbcYgKHd\nO8GBA3ZOhBcDg6Agynv0JGH/Hj7fvN9776uUQ9wNjA9EZAswHPhcROKAMs+VpdTprdtTjAgMTuoI\naWm20Qsd3vWFn5tCv+I85qXlePV9lXKCu8ubPwycB4wwxlQCx4CpnixMqdNZl1NM77gIIsOC7e0o\nEbvooBdJ/350L8pj6dZ9FB4t9+p7K+VtzdlB7xzgOhG5GZgGXOKZkpQ6PWMMa10d3oDt8B4wADp4\nebR3//4EVVXStXg/H23M9+57K+Vl7o6SehP4M3A+MNL1oavUKsfkFJVSdKyCod2jwBivd3jX6W83\nn/xeTSEL1+Wd5mSlWrfTrSVVawSQYowxnixGKXetd3V4pyZGQU4O7N/vTGD06wfAZSGHmb37EDlF\nx0mK1nU5Vdvk7i2pTUBXTxaiVHNk5B0hOFDo16XDdyvUernDG4C4OIiKYshxO0pq0bq93q9BKS9x\nNzBigUwR+UREFtd+eLIwpZqSkXeYfl06EBIUYG9HBQfD4MHeL0QEUlKI2LGVkcmdWLguD70QV22V\nu7ekHvdkEUo1hzGGzLwjTBjgmtG9ejUMGQKhoc4UlJoK//43U/+YwK8XZZCZf4SBCR2dqUUpD3J3\nWO1X2Bnewa6vVwNrPFiXUqdUcKSMwmMV9odyTY2dg+FE/0WtwYPh8GGuiKokKEBYpJ3fqo1yd5TU\nHcB84F+upm7AQjeeN0lEtorIDhF5uJHjIiLPuo5vEJFhrvYkEflSRDJFJENE7nf/j6Tauoy9RwC7\nYixbt0JJifOBAXTcsYXx/eNYvC6P6hq9LaXaHnf7MO4GxgFHAIwx24EmV3gTkUDgOeAyIAW4XkRS\nGpx2GdDX9XEn8IKrvQr4hTEmBRgD3N3Ic5Wfysg7gggMiI90tsO71qBB9vP69Uwd0o2CI2V1y64r\n1Za4GxjlxpiK2gciEgSc7leoUcAOY8wu13PncvLs8KnALGOtAKJEJN4Yk2+MWQNgjCkBNmOvapQi\nI+8wPWPa0z40yPZfRETUzYdwRIcO0KsXbNjA9wd0ITwkkMXrdbSUanvcDYyvRORXQDsRmQjMA94/\nzXO6AfUX2Mnl5B/6pz1HRJKBoYCuIa0Ae4UxsJurU3n1ahgxAgId3ig1NRXWr6ddSCCXDuzKRxsL\ndAVb1ea4GxgPAweAjcBPgI+AX3uqqFoiEgG8C/zMtbx6Y+fcKSJpIpJ24MABT5ekHFZ8vIK9xaW2\n/6KiAtatc7b/otbgwbB9Oxw/zpVDEjhcWsnX23QnPtW2uDtKqgbbyX2XMWaaMeZlN2Z97wWS6j1O\ndLW5dY6IBGPDYrYxZkETtb1kjBlhjBkRFxfnzh9HtWKZefU6vDdssKHhC4GRmmqXKNmwgfP7xBLd\nPkQn8ak2p8nAcI1ielxEDgJbga2u3fZ+48Zrrwb6ikhPEQkBpgMNJ/stBm52vc8Y4LAxJl9EBHgV\n2GyM+Wuz/1SqzcqoC4yO323J6mSHd63Ro+3nlSsJDgxgyuB4/rt5H0fLq5ytS6kWdLorjJ9jR0eN\nNMZEG2OigdHAOBH5eVNPNMZUAfcAn2A7rd8xxmSIyAwRmeE67SNgF7ADeBm4y9U+DrgJuFhE1rk+\nJp/Bn0+1MZvyDhPfMYzo9iF2hFRcHHTv7nRZkJAAiYmwYgUAU4ckUFZZw6cZBQ4XplTLOd1M75uA\nicaYupuxxphdrv28PwX+1tSTjTEfYUOhftuL9b422CG7DZ+3FJDTVq/8TkbeEXs7Cr6bsCc+8l9l\nzJi6wBjWvROJndqxaF0eVw9LdLgwpVrG6a4wguuHRS1jzAEg2DMlKdW40opqdh04SkpCRzh2DDIz\n7QgpXzF6NGRnw759iAhXpiawdMdBDurGSqqNOF1gVJzhMaVa3OaCI9QYV4f3+vV2WRAv77DXpDFj\n7OeVdgT41CHdqK4xurGSajNOFxipInKkkY8SYJA3ClSqVm2H97nd6u3h7UtXGMOGQVBQXWD079qB\nc7p20LWlVJvRZGAYYwKNMZGNfHQwxugtKeVVmXmHiQoPJqFjmA2M+Hjb2ewrwsPt8Nply+qapg7p\nRrprYyWlWrvm7OmtlKNqO7xFxAaGL11d1Bo/HpYvh9JSAK5IjQdg8Xq9ylCtnwaGahUqq2vYUlBi\n51+UlMCWLb4ZGBddZCcTLl8OQGKncEYmd2Kx3pZSbYAGhmoVdh44SkVVje3wXrvWzqr2pQ7vWhdc\nYNe1+vLLuqYrh3Rj674SNuc3urqNUq2GBoZqFU7YAyM93Tb6YmBERtorny++qGu6fFC8bqyk2gQN\nDNUqZOQdISw4gJ6xEbb/IjERunZ1uqzGXXSRnYV+9CgA0e1D+F6/OBau3UtVdY3DxSl15jQwVKuQ\nkXeYc7pGEhjgwx3etSZMgKoq+OqruqZrRyRRcKSMr7bpisqq9dLAUD7PGENmvmtJkMOHYds23w6M\nCy6A9u3hww/rmiYM6ExsRChzVuU08USlfJsGhvJ5OUWllJRV2RFSa9bYRl/sv6gVGgoTJ8IHH9jO\neSA4MIAfjkjky637KThc5nCBSp0ZDQzl8zLyDgOuDu/aGd6+HBgAU6ZATg5s3FjXNH1kEtU1hnlp\nepWhWicNDOXzMvOPEBgg9O/awY6Q6tHDLmvuyya7VuOvd1uqR0x7zusdw9tpOdTUnG7/MaV8jwaG\n8nkZeUfoHdeesOBA3+/wrhUfb6+CFi48oXn6qO7kHirVzm/VKmlgKJ+XkXfY9l8cOgQ7d/r+7aha\n06bZ4bXZ2XVNkwZ2pUtkKK8uzXKuLqXOkAaG8mkHj5az70i57b9oDR3e9V17rf38zjt1TSFBAdxy\nXjJLdxxkS4HO/FatiwaG8mm1S5qn+PoM78b06mV3BHz77ROabxjVnXbBgbz6jV5lqNZFA0P5tLoR\nUvGuPTCSkyEmxtmimuO66+yV0fbtdU1R4SFMG57IonV5HCjR3fhU66GBoXxaRt4RukW1o2N4sL3C\naC1XF7Vqb0vNnXtC8+3n96Sypkb7MlSrooGhfNqG3GJSk1wd3rt2tb7ASEqye2S8/rrdUtalZ2x7\nrkxNYNbybAp1z2/VSmhgKJ9VeLScnKJSUhOjWl+Hd3133AFZWSesYAtw78V9Ka2s5mXty1CthAaG\n8lkb9tr+i8GJUa2vw7u+q6+G6Gh4+eUTmvt0jtCrDNWqaGAon7U+pxgRGJTY8bsZ3q2pw7tWWBjc\ndBO89x4cOHHCXu1VxgtLdjpUnFLu08BQPmt9TjF9O0cQERrUOju867vjDqishFdfPaG5T+cIfjg8\nkTeWZ7O78JgztSnlJg0M5ZOMMWzIPWxvRxUXt64Z3o0ZONCuYPvMM1B+4u2nX1zSn+DAAP7w8RaH\nilPKPRoYyiflHiql8FgFqUmtvMO7vocegoICeOutE5q7RIYx48LefLypgFVZRQ4Vp9TpaWAon7Q+\ntxiAIa29w7u+CRNgyBD4859PGGILcMcFvegaGcbvPsikWleyVT5KA0P5pLV7igkJCvhuSfPu3SE2\n1umyzo6IvcrYsgXeffeEQ+1CAnlk8jls3HuYf6/c7VCBSjVNA0P5pNXZRQxJiiIkKMAGRmtY0twd\n115r+zMefdR2gtdzZWoC4/rE8PR/trK/RHflU75HA0P5nGPlVWTkHWFUcrTt8N6xo/XfjqoVGAhP\nPWXXlnrttRMOiQi/m3ou5VU1/O6DzQ4VqNSpaWAon7NmzyGqawyjekbb/SQARo92tqiWNGUKjBsH\njz8OJSUnHOoVF8FdF/Xm/fV5fK2bLCkfo4GhfM7qrCICBIb16AQrVth7/yNHOl1WyxGBv/4V9u2D\nX//6pMMzLuxNz9j2PLZoE2WV1Q4UqFTjNDCUz1mZVcTAhI52wt6KFZCSApGRTpfVskaNgrvugn/8\nA1avPuFQWHAgT151LrsLj/Ps59tP8QJKeZ8GhvIp5VXVrMspZmRyNBgDK1fCmDFOl+UZv/+93fv7\nxz+GshOvy4oPAAAVwElEQVQ7ucf1iWXa8ET+9fWuuj1BlHKaRwNDRCaJyFYR2SEiDzdyXETkWdfx\nDSIyrN6x10Rkv4hs8mSNyrekZR+ivKqGcX1ibGd3UVHbDYzISHjpJdiwAX75y5MO//ryAXQKD+GX\n726gqrqmkRdQyrs8FhgiEgg8B1wGpADXi0hKg9MuA/q6Pu4EXqh3bCYwyVP1Kd/09bYDBAcKY3rF\n2NtR0HYDA+Dyy+FnP4Nnn4XFi084FBUewhNXDmTT3iO8ohstKR/gySuMUcAOY8wuY0wFMBeY2uCc\nqcAsY60AokQkHsAY8zWg6yT4ma+3H2REj2ja1/ZfdOgAAwY4XZZn/eEPMGyYXdF2y4nrSU0e1JVL\nUrrwt8+2kXVQFydUzvJkYHQDcuo9znW1NfecJonInSKSJiJpBw7oMMTWbH9JGZvzj3BBP9eM7hUr\nbOdwYKCzhXlaaKhd+jwsDK64wt6GcxERfnfVuYQEBfDIgg0Yo8uGKOe0+k5vY8xLxpgRxpgRcXFx\nTpejzsLnm/cDML5fZzhyBNatg/POc7gqL+ne3YbGnj12Nni9WeBdIsP41eQBrNhVxNzVOU28iFKe\n5cnA2Ask1Xuc6Gpr7jnKT3y0MZ/kmHAGxHeApUvtAn3jxztdlvecd57tBP/8c7j//hMOTR+ZxJhe\n0fz+w80UHNZlQ5QzPBkYq4G+ItJTREKA6cDiBucsBm52jZYaAxw2xuR7sCblow4dq+DbnYVcNige\nEYElSyAkpG13eDfmllvsAoUvvAD//Gdds4jwh6sHU1Fdw2OLNumtKeUIjwWGMaYKuAf4BNgMvGOM\nyRCRGSIyw3XaR8AuYAfwMnBX7fNFZA6wHOgvIrki8iNP1aqc92lmAdU1hsnnxtuGJUvsciDh4Y7W\n5Yjf/x6uvNJeZXz6aV1zcmx7HpjYj88y9/HRxgIHC1T+StrSbyojRowwaWlpTpehzsAPnl/G4dJK\nPn/gQqSkBDp1siu6/r//53RpzigpgfPPh927YfnyupFiVdU1/OD5b8k/XMp/H7iQqPAQhwtVrZ2I\npBtj3FoOutV3eqvWLzPvCGv3FPM/o3vY21H+2H/RUIcOdl5GaKgdOVVYCEBQYAB/uGYQh45X8uSH\nuqKt8i4NDOW4N1fsJjQogGuGuUZUf/KJHWLqb/0XDfXoAQsXQk4OTJtWN3JqYEJHfvK9XsxPz+Wb\n7TqUXHmPBoZy1N7iUuan53DN8ER7e8UY+PBDuPhi/+y/aGjsWHjlFdun89hjdc33TehLr9j2PLJg\nI0fLq5yrT/kVDQzlqL98shVBuPuiPrZh2zbYudMumaGsm26CO++EP/4RPvgAsCva/nHaYPKKS3l8\ncYbDBSp/oYGhWoQxhk17D/Pe2lz+symf3YWnX8biyy37WbB2L3d+rxfdotrZRtcPRA2MBp55BoYM\ngZtvth3hwMjkaO65qA/z03NZvD7P4QKVPwhyugDV+q3OLuKxhZvYUtBw97j2TE3txlVDE+gR0/6E\nYxtyi7lvzlrO6dqBey7u892BBQtg8GB7/159JywM5s2za05ddx18/TWEhHDfhL4s3XGQR9/byNCk\nKJKi9Tae8hwdVqvOyry0HB5esJH4jmHce3EfhvfoRFllDem7D/HxpnxW7LLrIg3tHsXU1AQSotqR\ntvsQb3ybTWxEKPNmjCWh9uoiKwt69bJ7Xj980mr4CmxoXHutnaPx978DkFN0nMue+YZ+XSKYe+dY\nQoL0xoFyX3OG1WpgqDP2WeY+fvJmGuP6xPL8/wyjQ1jwSefkFZeyeH0eC9furbsCCRCYMjiB316R\nQkxE6Hcn/9//2S1Ls7P1CqMp991nd+qbN8+OngI+3JDP3f9eww2ju/P7HwxyuEDVmmhgKI/bXXiM\ny575hj6dI5h75xjCQ05/d3N34TGKj1eSHNOejuENwqWmxk5O69LF3m5Rp1ZRAd/7HmRmQno69O0L\nwB//s4UXluzkyavO5cYxGrjKPTpxT3lUdY3hgXfWExggvHjDUMLffMN2yIaGQlKSvV2y9+Q1JHvE\ntCc1KerksAC7BMa2bXY0kGpaSAi88w4EB9srjOPHAXjwkv5c1D+O3y7O4Ist+xwuUrVFGhiq2eau\n3kP67kP8blIfEn5yq92TOjjYBsXo0fD883DOOfDmm+6/6J//bPe3vvZaj9XdpnTvDrNnw8aNcM89\nAAQGCM9eP5SU+Eh++tYaVmXp/mOqZWlgqGY5XFrJXz7dxqgenZj6rydh/nz4059g1Sp4+mn7eOtW\nO5rn5pvhttugtLTpF/3sM7uk9//+r/3tWbln0iTb5/P663Ur23YIC2bmbSPp1qkdP5q5mvTdGhqq\n5WhgqGZ57ssdHDpewd9K0pCZM+3s4wcfBJHvTurVC774wh6bOdMuoped3fgLHjtmf0NOToa77mr8\nHHVqv/2tXdn2vvvsMiJATEQob/1oNDERIdz06iqW7TjocJGqrdDAUG47UFLOG99mc0uPYLo9+Wu7\nfMfjjzd+cmCgXWn2/fftzO3hw+0aUfVVVdkrkO3b4bXXbB+Iap7AQJgzB0aOhOuvtws3AglR7Xhn\nxli6R4dz2+ur+c8m3WZGnT0NDOW2V77ZRWV1DQ9++rJdCO+llyDgNP+FpkyBtDRISLC3UK6+Gt5+\nG+bOtavRzptnb2lddJFX/gxtUni4Debu3e3f8ZIlAHTuEMbcO8cwsFskP529hheW7NSNl9RZ0cBQ\nbik6VsGbK3ZzR6fjRLw3H37xC+jd270n9+kDK1bY2ydffAHTp9vfhrdtg1mz7Gups9O5sw2KHj3g\nsstsXxIQFR7CnDvGcPmgeP74ny387/wNVFTVOFurarU0MJRbXluaxfGKau5e/jZERMADDzTvBdq3\nt7ev9u+3cwfS0yEvzy6sp1pGfLwNjSFD4Ic/tB3iVVWEBQfyj+uHcv+EvsxPz+XGV1ZSdKzC6WpV\nK6SBoU7rcGklb3ybzW0xZUQuWgB33w3R0Wf2YiEhdgTVsGEQpEuZtbi4OBsat91mZ85fcAFs346I\n8POJ/Xj2+qGszy1m6nNL2bav5LQvp1R9GhjqtGYuy6akvIr7V86Hdu2af3WhvCs0FF591c7T2LIF\nUlPtPuHl5VyZmsDbPxlLeWUNVz//rU7wU82igaGaVFJWyWvLsrghuoKo996BGTPs/XLl20Tghhsg\nI8P2aTz6KAwcCB9+yJCkKBbdM47k2HB+9EYaL3+9SzvDlVs0MFST3lyxm8OllfxizQJ7O+nBB50u\nSTVHQgK8+65deiU42I5amzKF+KIC3vnJWC47tyv/99FmHpq/gfKqaqerVT5OA0Od0rHyKl75Jotp\nMZXEvDsX7rjDdqyq1mfiRFi/3i7B8tVXMHAg4X//K/+cNoj7JvRlnqszvPBoudOVKh+mgaFO6a0V\nuyk6VsEv1y+y8y0eesjpktTZCAmxQ5gzM+GSS+CXvyRg1Ege6HSEf1w/lA25h5n63DK2FmhnuGqc\nBoZqVGlFNS99vYupMdXEvTMbbr8dEhOdLku1hKQku4zIggVw8CCMHcsVH8/inR+PoqKqhqufX8bn\nm7UzXJ1MA0M1avbK3RQeq+BXm94HY3QHvLboBz+wq91edRU88gipt03jgx8k0ysugh/PSuNfX+nM\ncHUiDQx1kpKySl5YspMroirpMncW3Hqr7oDXVkVH2701XnsNVq+m8/mjmN91H5PPjeepj7fw4Dzt\nDFff0cBQJ3l+yU4Kj1XwxEa7+imPPeZsQcqzROxEv7VroXdvQq+dxj+/eYlfjEvk3TW5XPevFeQU\nHXe6SuUDNDDUCXKKjvPq0izuSKgh+p3Zdge87t2dLkt5Q9++sGwZPPQQ8q9/ce/DN/DmyHbs3H+U\nyc98w6J1J++iqPyLBoaqY4zhyQ8zCRB44IvX7bj9X/3K6bKUN4WEwB//aJeiLyzkgv+ZzNfhGfTr\nEsH9c9dx97/XUHC4zOkqlUM0MFSdDzbk80nGPv7UcT/tFr0Hjzyi8y781SWXwIYNMGECnR7+BfMW\nPsGTA4L4b+Y+JvxlCS8s2cnxiiqnq1ReJm1pFMSIESNMWlqa02W0SvuOlHHZM9/Qu30A7/zrLqS6\nGjZtgrAwp0tTTjIGnnvOLi1y/DhH7ryL3/SfzMK8KmLah/DjC3pxw6judAwPdrpSdYZEJN0YM8Kt\nczUwVHlVNdNfWsHWghK+zV1A1Osv26UkJk50ujTlK/bts7cnX38dQkPZf810nus1njfKYwgNDmTS\nuV25dkQSo3tGExSoNy5aEw0M5bbqGsOD89bz3tq9LIzKZsgj99jVaP/yF6dLU75oyxb7f2PWLKio\noLxnb9LPHcvsiL581WUAwZ06cvE5XbhkYBfG9YklIlSXsPd1GhjKLRVVNTy8YAML1uzl2ej9XPnY\nDBg71nZ46v7aqilFRXam+NtvwzffQHk5JiCAvMTeLI/tw6rOfdmY2J/2qecypm8XzusTw7DunQgL\nDnS6ctWAzwSGiEwCngECgVeMMX9ocFxcxycDx4FbjTFr3HluYzQw3Jd18BgPzltPenYRL5sMJv79\nMRgwwG6+06mT0+Wp1qSsDL791v7fWbkSs2oVUlwMQGloO9Z26cu6hH5kdDuH8qHDSDy3L6ndoxic\nGEXPmPYEBIiz9fs5nwgMEQkEtgETgVxgNXC9MSaz3jmTgXuxgTEaeMYYM9qd5zZGA+P08opLmflt\nNjO/zWbQwWye27qQrks+hYsusr8xRkU5XaJq7WpqYPt2WLkSVq2ievkKZMN6AqrsqKrDYRFsi+nO\n9tjuZMclYbp1I7RHEpG9k4nu24OusR3p2jGM+I5htNdbWh7XnMDw5L/GKGCHMWaXq6i5wFSg/g/9\nqcAsY1NrhYhEiUg8kOzGc1UTKqtrOHS8gqJjFezcd5RdWflsWr+TsswtDCrYwaf560jevhE6dICn\nn7b9FoF6u0C1gIAA6N/fftx8M4Fgr0LWroX0dDps2sTAdRsZvHk5oev/c9LTS0LacahdJDvaRVIS\n0ZHyyCgqo2Oo6RSNiYkhMC6W4C6dCevamfbdutChSyxhkRGEh4XQLiSQ0KAA7M0L1dI8GRjdgJx6\nj3OxVxGnO6ebm89tMTuS+hNcUQ7Yqy1xXXWJq6Xucd3V2Mnn2eaG7abu7IbnijEN2mufW/dijb5W\n7aeTn9fgStEYwoAEILm6ksurKk48PnQoPPUU/OQnegtKeV5YmO0fGzuWACAc7PdLURHs3Qt791K1\new8lu3MpKziA7NtPp8JCYg8VEZKTSXjmYdqXHWvyLSoCgigPCuZoUAiVQcGYgED73SGCEQFcn+se\nU6/do3966v2U8IjjHTqSsn2dR98DPBsYXiEidwJ3AnQ/wyUsinv0JqCyEgRM7T+snPjZiNgAqf+b\ni+s/mz1Oo8/jpPOxz5ETX/uk4029Vt1pASfUVv/r4KBAwoIDCA0JpENEOIHJ3Qju2gV69rR7POut\nJ+U0EYiJsR+DBxMENPmrS0UFFBVRXrCfo3sLOJa3j7KC/VQVHqK6tJSa0lJMaZm9mikrw1RVY0yN\n/eXMGBtQNfZz7WPB1P2i15pVdejolffxZGDsBZLqPU50tblzTrAbzwXAGPMS8BLYPowzKXTE0o/O\n5GlKKW8KCYGuXQnt2pXQIYOJcboeP+TJGTargb4i0lNEQoDpwOIG5ywGbhZrDHDYGJPv5nOVUkp5\nkceuMIwxVSJyD/AJdmjsa8aYDBGZ4Tr+IvARdoTUDuyw2tuaeq6nalVKKXV6OnFPKaX8WHOG1eqi\nL0oppdyigaGUUsotGhhKKaXcooGhlFLKLRoYSiml3NKmRkmJyAFgdzOfFgsc9EA5Z8tX6wLfrU3r\naj5frU3rar4zra2HMSbOnRPbVGCcCRFJc3dImTf5al3gu7VpXc3nq7VpXc3njdr0lpRSSim3aGAo\npZRyiwaGa+FCH+SrdYHv1qZ1NZ+v1qZ1NZ/Ha/P7PgyllFLu0SsMpZRSbvHrwBCRX4iIEZHYem2P\niMgOEdkqIpc6UNOfRGSLiGwQkfdEJKreMadrm+R67x0i8rC3379eHUki8qWIZIpIhojc72qPFpHP\nRGS767MjWwmKSKCIrBWRD3ysrigRme/6/7VZRMb6Qm0i8nPXv+MmEZkjImFO1SUir4nIfhHZVK/t\nlLV463vyFHV5/WeF3waGiCQBlwB76rWlYPfeGAhMAp4XEW9vdP0ZcK4xZjCwDXjEF2pzvddzwGVA\nCnC9qyYnVAG/MMakAGOAu121PAx8bozpC3zueuyE+4HN9R77Sl3PAP8xxpwDpGJrdLQ2EekG3AeM\nMMaci93OYLqDdc3Efn/V12gtXv6ebKwur/+s8NvAAP4GPETdLtkATAXmGmPKjTFZ2H06RnmzKGPM\np8aYKtfDFdjdBn2htlHADmPMLmNMBTDXVZPXGWPyjTFrXF+XYH/wdXPV84brtDeAq7xdm4gkApcD\nr9Rr9oW6OgLfA14FMMZUGGOKfaE27L487UQkCLvdd55TdRljvgaKGjSfqhavfU82VpcTPyv8MjBE\nZCqw1xizvsGhbkBOvce5rjan3A587Pra6dqcfv9GiUgyMBRYCXRx7dgIUAB0caCkv2N/Eamp1+YL\ndfUEDgCvu26XvSIi7Z2uzRizF/gz9ko/H7vr5qdO19XAqWrxpe8Jr/ys8OSe3o4Skf8CXRs59Cjw\nK+ztKEc0VZsxZpHrnEext15me7O21kREIoB3gZ8ZY46ISN0xY4wREa8OARSRKcB+Y0y6iIxv7Bwn\n6nIJAoYB9xpjVorIMzS4zePQ31kn7G/EPYFiYJ6I3Oh0XafiS7XU8ubPijYbGMaY7zfWLiKDsP85\n17t+wCQCa0RkFLAXSKp3eqKrzSu11avxVmAKMMF8N+7ZK7U1wen3P4GIBGPDYrYxZoGreZ+IxBtj\n8kUkHtjv5bLGAVeKyGQgDIgUkbd8oC6wv2XmGmNWuh7PxwaG07V9H8gyxhwAEJEFwHk+UFd9p6rF\n8e8Jb/+s8LtbUsaYjcaYzsaYZGNMMvYbaZgxpgBYDEwXkVAR6Qn0BVZ5sz4RmYS9pXGlMeZ4vUNO\n17Ya6CsiPUUkBNupttiL719HbNK/Cmw2xvy13qHFwC2ur28BFnmzLmPMI8aYRNf/q+nAF8aYG52u\ny1VbAZAjIv1dTROATB+obQ8wRkTCXf+uE7B9Uk7XVd+panH0e9KRnxXGGL/+ALKB2HqPHwV2AluB\nyxyoZwf2/uM618eLPlTbZOxojJ3Y22dO/Zudjx2ssKHe39NkIAY7imU78F8g2sEaxwMfuL72ibqA\nIUCa6+9tIdDJF2oDngC2AJuAN4FQp+oC5mD7Uiqxv0z+qKlavPU9eYq6vP6zQmd6K6WUcovf3ZJS\nSil1ZjQwlFJKuUUDQymllFs0MJRSSrlFA0MppZRbNDCUUkq5RQNDKaWUWzQwlFJKueX/Aw+WAG0C\nILgCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xc46f7874a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we can see that the distribution has changed slightly with now more values accumulating towards the median\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "X_train['Age'].plot(kind='kde', ax=ax)\n",
    "X_train.Age_median.plot(kind='kde', ax=ax, color='red')\n",
    "lines, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(lines, labels, loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, the median imputation distorts the original distribution of the variable Age. The transformed variable shows more values around the median value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xc46fcd3f98>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAD8CAYAAABdCyJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8lfXZ+PHPlc0IhEACIQkk7C0jIggKDqYDqQutq9pa\nt7X2qbY+tvXXxz5WW+solWrd9QGtC2xRVMSBzIARCDOMQCCBEEYIEEKS6/fHfQIhOUlOkrOSXO/X\n67zOOff9vc99HSW58t2iqhhjjDENFRLoAIwxxjRtlkiMMcY0iiUSY4wxjWKJxBhjTKNYIjHGGNMo\nlkiMMcY0iiUSY4wxjWKJxBhjTKNYIjHGGNMoYYEOwB86deqkKSkpgQ7DGGOalFWrVu1X1bi6yrWI\nRJKSkkJ6enqgwzDGmCZFRLI9KWdNW8YYYxrFEokxxphGsURijDGmUVpEH4kxpvk6efIkOTk5FBcX\nBzqUJisqKoqkpCTCw8MbdL1PE4mITAaeBUKBf6jqE1XOi+v8VOAYcIuqrhaRKOBrINIV47uq+lvX\nNbHA20AKsAO4RlUP+vJ7GGOCV05ODtHR0aSkpOD8SjH1oaoUFBSQk5NDampqgz7DZ01bIhIKzASm\nAAOA60RkQJViU4DersftwAuu4yeAC1X1LGAoMFlERrnOPQwsVNXewELXe2NMC1VcXEzHjh0tiTSQ\niNCxY8dG1eh82UcyEshS1W2qWgLMAaZVKTMNeEMdy4AYEUlwvS9ylQl3PbTSNa+7Xr8OXOHD72CM\naQIsiTROY//7+bJpKxHYVel9DnCOB2USgVxXjWYV0AuYqarLXWU6q2qu63Ue0NndzUXkdpxaDt26\ndWvE1zAthaqyfPsBvtt5iLaRoVzYvzOJMa0CHZYxQS9oR22papmqDgWSgJEiMshNGeV0TaXquRdV\nNU1V0+Li6pyYaVq4gqIT3PTKCma8uIw/frKRR+dmMu7JRfz5002Ulbv9J2bMGT788ENEhI0bNwY6\nFL/zZSLZDSRXep/kOlavMqp6CFgETHYd2isiCQCu531ejNm0QEUnSrnplRWs2H6A3102gLW/m8iX\nvxjP5UO78vwXWTzwdoYlE1On2bNnM3bsWGbPnh3oUPzOl4lkJdBbRFJFJAKYAcyrUmYecJM4RgGH\nVTVXROJEJAZARFoBE4CNla652fX6ZmCuD7+DaQF+M3cdG/OOMOvGEdwyJpXoqHBSOrXh6WuG8svJ\nfZn3/R7+8tnmQIdpglhRURGLFy/m5ZdfZs6cOQCUl5dz11130a9fPyZMmMDUqVN59913AVi1ahXj\nxo1jxIgRTJo0idzc3No+Puj5rI9EVUtF5B5gAc7w31dUNVNE7nCdnwXMxxn6m4Uz/PdHrssTgNdd\n/SQhwDuq+m/XuSeAd0TkNiAbuMZX38E0f4u37Of91bu554JeXNA3vtr5u8b3Ysf+o8z8MotRPToy\ntnenAERpPPXYR5ms31Po1c8c0LUdv71sYK1l5s6dy+TJk+nTpw8dO3Zk1apVbN++nR07drB+/Xr2\n7dtH//79ufXWWzl58iT33nsvc+fOJS4ujrfffptHHnmEV155xatx+5NP55Go6nycZFH52KxKrxW4\n2811a4BhNXxmAXCRdyM1LZGq8uSCjSTGtOLei3rVWO6xyweRnn2QX3+wlk8fOJ+o8FA/Rmmagtmz\nZ3P//fcDMGPGDGbPnk1paSlXX301ISEhdOnShQsuuACATZs2sW7dOiZMmABAWVkZCQkJAYvdG2xm\nu2mxvt6ynzU5h3nyyiFEhtWcHFpFhPL7aYP44T+W8/evtnH/xb39GKWpj7pqDr5w4MABvvjiC9au\nXYuIUFZWhogwffp0t+VVlYEDB7J06VI/R+o7QTtqyxhf++eybDq1jeCKYYl1lh3TqxNTBnXhpW+2\ncehYiR+iM03Fu+++y4033kh2djY7duxg165dpKamEhsby3vvvUd5eTl79+7lyy+/BKBv377k5+ef\nSiQnT54kMzMzgN+g8SyRmBYp9/BxFm7YyzVpyUSEefZj8LOL+1B0opSXF2/3cXSmKZk9e3a12seV\nV15JXl4eSUlJDBgwgBtuuIHhw4fTvn17IiIiePfdd3nooYc466yzGDp0KEuWLAlQ9N5hTVumRXp/\n9W7KFa4b6flk1b5dopk6uAuvfbuDH4/tQfvWDVvgzjQvixYtqnbsvvvuA5zRXG3btqWgoICRI0cy\nePBgAIYOHcrXX3/t1zh9yWokpkWavzaX4d1iSI5tXa/r7rmgN0dOlDJn5U4fRWaak0svvZShQ4dy\n3nnn8eijj9KlS5dAh+QTViMxLc7OgmNk7inkkan9633tgK7tGJkayz+XZ/Pj83oQGmJrPJmaVfSL\nNHdWIzEtzoLMPAAmD2rYX4c3j05h14HjfLnJFlUwBiyRmBZo4ca99OsSXe9mrQoTB3amc7tI3lia\n7eXIjGmaLJGYFuV4SRmrsw9xfp+GL+QZHhrC1SOS+WZLPvsKbVc+YyyRmBZl5Y4DlJSVc27Pjo36\nnOnDEylXmJuxx0uRGdN0WSIxLcq3W/cTHiqMTI1t1Of0jGvLWckxvLc6x0uRGdN0WSIxLcrSrQUM\n69aB1hGNH7B45fBENuYd8foigaZpsv1IjGkBjpWUkrmnkHMaWRupcOmQroSHCh9mVN1mx7REwbYf\nSWlpqd/uZfNITIuxJucwZeXK8G4dvPJ5sW0iGNOrE/PX5vKrKf1s3/Bg8LOfQUaGdz9z6FB45pla\ni1TsR7Jo0SIuu+wyHnvsMcrLy7nnnnv44osvSE5OJjw8nFtvvZWrrrqKVatW8fOf/5yioiI6derE\na6+95nYF4D179jB16tRT79euXcu2bdto3bo1d9xxBzt3OhNjn3nmGcaMGcPvfvc7tm7dyrZt2+jW\nrRuvvvoqd955J+np6YSFhfH000+fWoXYmyyRmBZj9c6DAAxNjvHaZ04dlMAv31tD5p5CBiW299rn\nmqbFV/uRdO3alQxXYpw5cyZfffUV3bt35/rrr+eBBx5g7Nix7Ny5k0mTJrFhwwYA1q9fz+LFi2nV\nqhV//vOfERHWrl3Lxo0bmThxIps3byYqKsqr398SiWkxVmcfokenNnRoE+G1z5wwoDOhHwgfr8u1\nRBIM6qg5+Iqv9yP59ttveemll1i8eDEAn3/+OevXrz91vrCwkKKiIgAuv/xyWrVqBcDixYu59957\nAejXrx/du3dn8+bNDBkyxIvf3hKJaSFUlYxdBxnXp/ouiI3RoU0Eo3rE8vG6PH4xsa81b7VAvt6P\nJDc3l9tuu4158+bRtm1bwNnGd9myZW5rFm3atGn4l2kg62w3LcLuQ8fZX1TC0G7ea9aqMGVQAtvy\nj7JlX5HXP9sEP1/uR3Ly5Emuvvpq/vjHP9KnT59TxydOnMjzzz9/6n1GDf1C5513Hm+99RYAmzdv\nZufOnfTt29cbX/sMlkhMi1AxRHdg13Ze/+yJAzsD8KlrDS/TsvhyP5IlS5aQnp7Ob3/7W4YOHcrQ\noUPZs2cPzz33HOnp6QwZMoQBAwYwa9Yst9ffddddlJeXM3jwYK699lpee+01IiMjvf7fQJxt05u3\ntLQ0TU9PD3QYJoCe+Xwzzy7cQuZjk7wyh6SqaX9dTGiI8P5dY7z+2aZ2GzZsoH//+q/k7A9V9yP5\n9ttvg3YpeXf/HUVklaqm1XWt9ZGYFmH9nkJSO7XxSRIBuKBfPM8u3MKBoyXEerEz3zRtl156KYcO\nHaKkpMT2IzGmqduQV8hZSd7vH6lwYb94nvl8C19t3sf0YUk+u49pWuqzH8ndd9/Nt99+e8ax+++/\nnx/96Edejsr7LJGYZu/w8ZPsOnC8Xtvq1tegru3p1DaSLzbmWyIJAFVt8iPmZs6cGbB7N7aLwzrb\nTbO3MdfpaB+Q4P2O9gohIcIFfeP4atM+SsvKfXYfU11UVBQFBQWN/mXYUqkqBQUFjZqk6NMaiYhM\nBp4FQoF/qOoTVc6L6/xU4Bhwi6quFpFk4A2gM6DAi6r6rOua3wE/AfJdH/NrVZ3vy+9hmrb1FYnE\nByO2KruwXzz/WpXDquyDnNOjccvUG88lJSWRk5NDfn5+3YWNW1FRUSQlNbwm7bNEIiKhwExgApAD\nrBSReaq6vlKxKUBv1+Mc4AXXcynwoCupRAOrROSzStf+RVX/5KvYTfOyMfcIHdtEEB/t3WUhqhrb\nuxNhIcKXm/MtkfhReHg4qampgQ6jRfNl09ZIIEtVt6lqCTAHmFalzDTgDXUsA2JEJEFVc1V1NYCq\nHgE2AIk+jNU0Y1v2HaF357Y+v090VDjDu3Xg26z9Pr+XMcHEl4kkEdhV6X0O1ZNBnWVEJAUYBiyv\ndPheEVkjIq+IiNulXEXkdhFJF5F0q/K2XKpK1r4iesX7PpEAjOnVibW7D3PoWIlf7mdMMAjqznYR\naQu8B/xMVSt2D3oB6AEMBXKBP7u7VlVfVNU0VU2Li2v4/tymacsvOkFhcSm94vyTSMb27oiqs4GW\nMS2FLxPJbiC50vsk1zGPyohIOE4SeUtV368ooKp7VbVMVcuBl3Ca0IxxK8u1/lWv+Gi/3G9IUgxt\nI8NYbM1bpgXxZSJZCfQWkVQRiQBmAPOqlJkH3CSOUcBhVc11jeZ6Gdigqk9XvkBEKq+3PB1Y57uv\nYJq6racSiX9qJOGhIYzqEWuJxLQoPkskqloK3AMswOksf0dVM0XkDhG5w1VsPrANyMKpXdzlOj4G\nuBG4UEQyXI+KbcKeFJG1IrIGuAB4wFffwTR9WfuKaBsZRud23l+oriZjenUiu+AYuw4c89s9jQkk\nn84jcc3vmF/l2KxKrxW42811iwG301RV9UYvh2masaz8InrGt/XrrOfzencC4Nus/czw4Wx6Y4JF\nUHe2G9NYWfuK/NbRXqFnXFs6t4u05i3TYlgiMc3WkeKT7C084bf+kQoiwphenViy1ZbtMC2DJRLT\nbO3Y7/RR9Ijz/9ajo3p05MDRklOjxoxpziyRmGZrR8FRAFI6ephIDhyABx+Erl1BBFJT4b//Gw4f\nrve9z0mNBWD59gP1vtaYpsYSiWm2sl2JpFts67oLb9gAI0bAM8/AmDHw6KMwcCA8/jicdRbUsCd2\nTbrFtqZzu0hLJKZFsP1ITLO1ff8xurSLolVEaO0F9+yBiRPh5ElYuhRGVprjumwZXH01XHABfPYZ\npNW56yjg9JOck9qR5dsLmsVeGcbUxmokptnKLjhK94511EZU4dZbnWatTz45M4kAjBoF33wDMTFw\n2WVO0vHQyNRY9haeYKfNJzHNnCUS02ztKDhWd//Ia6/BggXw5JMwdKj7Mikp8NFHcOSIUzspK/Po\n/qf6SbZZ85Zp3iyRmGap6EQp+4tO0L1TLTWS48fhkUdg9Gi4887aP3DQIJg1C5Ysgeee8yiGXvFt\niW0TYf0kptmzRGKapWxPRmzNmgW5ufCHP0CIBz8KP/whXHKJk3y2bq2zuIgwMiWWFTtsJWDTvFki\nMc1SdoHTL1FjH8nJk/DUU04n+vjxnn2oiJN8QkLgoYc8umRkaiy7Dhxnz6Hjnt3DmCbIEolplrbv\nr6NGMneuUxt58MH6fXBSEvzXf8F77zkjvOow0tVPssKat0wzZonENEvZBUeJi46kTWQNI9z/9jen\nE33y5Pp/+IMPQpcuTkKpYwmU/gntiI4KY/l2a94yzZclEtMsOSO2amjW2roVFi2C22+H0DrmmLjT\ntq0z4/3bb+Grr2otGhoijOjegVXZB+t/H2OaCEskplnKLjhKt9gamrXmzHGef/jDht/g1lshPh7+\n93/rLDqiWwc27y3i8PGTDb+fMUHMEolpdk6UlrG38ETNS6PMmQNjx0K3RuwV0qoVPPAAfPoprFpV\na9ER3TsAsHqn1UpM82SJxDQ7ew4VA5DUoVX1k5mZsG4dzJjR+BvdeSe0a+dMZqzFWckxhIYIq615\nyzRTlkhMs5Nz0Bn66zaRzJ3rPP/gB42/Ufv28JOfwPvv17p0SpvIMPonRFs/iWm2LJGYZifnoDNn\nI8ld09b8+c4qvwkJ3rnZnXdCaSm8+GKtxUZ060DGrkOUlpV7577GBBFLJKbZyTl4jLAQoXN05Jkn\nCgqcuR9Tp3rvZj17OkOIX3zRmeRYgxEpsRwrKWNj3hHv3duYIGGJxDQ7OQePkxATRVholX/en34K\n5eXOMifedPfdzuTGDz+ssUhFh7s1b5nmyBKJaXZyDh4nKaaGZq1OnTzeU8RjU6Y4I8BefrnGIokx\nrUhoH0W6JRLTDFkiMc1OzsFj1TvaVZ2NqSZObNgkxNqEhsLNNzufv3t3jcWGd+9gI7dMs2SJxDQr\nFXNIkjpUqZFs2gR79zqLNPrCzTc7zWZvvlljkRHdOrD70HFyD9sCjqZ58WkiEZHJIrJJRLJE5GE3\n50VEnnOdXyMiw13Hk0VkkYisF5FMEbm/0jWxIvKZiGxxPXfw5Xcw9ZSZCT/9qTMkNgBqnEPy5ZfO\ns6cr/dZXz57OJMfXXqtx/a20FOsnMc2TzxKJiIQCM4EpwADgOhEZUKXYFKC363E78ILreCnwoKoO\nAEYBd1e69mFgoar2Bha63ptgsGcPnH++M4LpyivhnXf8HkKNc0i+/BISE51f+L5yyy1OzWfFCren\n+ye0Iyo8xBKJaXZ8WSMZCWSp6jZVLQHmANOqlJkGvKGOZUCMiCSoaq6qrgZQ1SPABiCx0jWvu16/\nDlzhw+9g6uOxx5ztaNeudTq0H3gASkr8GoLbOSSqTiIZP97ZU8RXrr7aWTrltdfcng4PDeGspBjr\nJzHNji8TSSKwq9L7HE4nA4/LiEgKMAxY7jrUWVVzXa/zgM7ubi4it4tIuoik5+fnNyR+Ux+FhfDP\nfzp9BYMGwe9/79RQ3n3Xr2G4nUNS0T/iq2atCu3awbRpzncuLXVbJC2lA5l7Cjle4tm+78Y0BUHd\n2S4ibYH3gJ+pamHV86qqgNsGaVV9UVXTVDUtLi7Ox5Ea5s6FY8ecVXEBJk1yNoHyc/OW2zkk33zj\nPJ9/vu8DuPZa2L8fvvjC7ekR3TtQWq58n3PI97EY4ye+TCS7geRK75NcxzwqIyLhOEnkLVWt3HO7\nV0QSXGUSgH1ejts0xCefOMuqn3OO814Epk+HBQvg6FG/heF2DsmyZRAbC717+z6AyZMhOhreftvt\n6eHdrMPdND++TCQrgd4ikioiEcAMYF6VMvOAm1yjt0YBh1U1V0QEeBnYoKpPu7nmZtfrm4G5vvsK\nxiPl5c6s8YkTnf3MK0yfDsXFzjk/cTuHZPlyGDXKt/0jFaKinOat99932z8U0zqCXvFtLZGYZsVn\niURVS4F7gAU4neXvqGqmiNwhIne4is0HtgFZwEvAXa7jY4AbgQtFJMP1qFgg6QlggohsAS52vTeB\nlJHhNOdMmnTm8TFjnM7niqG3PuZ2Dsnhw7B+vZNI/OXaa+HQIWeCohtp3TuQvuMA5eW1b9NrTFNR\nw4bW3qGq83GSReVjsyq9VuBuN9ctBtz++aiqBcBF3o3UNMrSpc7zuHFnHo+IgNGj4euv/RKG2zkk\nK1c6o7b8mUgmToSYGKd5y826XiO6d2DOyl1k5RfRp3O0/+IyxkeCurPdNBHLl0OXLk7nelXjxsH3\n3zt/ofuY2zkky5Y5TVojR/r8/qdERDjNeh9+6DTtVZGWEgtA+g5r3jLNgyUS03grVji/qN31QZx3\nnlMjqKi1+JDbOSTLl0O/fs4mVP509dXOnJqFC6udSunYmo5tIkjPPuDfmIzxEUskpnEOHXLmaVSM\n1qpq+HDnuY59zb2h2hwSVadG4s9mrQoXXuiM3nKztLyIMKJ7B+twN82GJRLTON9/7zxXJIyq2rd3\nht2uXu3zUKrNIdm+3RkEUFOS86XISGcDrXnzoKz65MOzU2LJLjhG/pET/o/NGC+zRGIaJzPTeR40\nqOYyI0b4qUZSZQ7Jd9+dvn8gTJ8O+/a5bdYbcWoBR2veMk2fJRLTOOvWObWOxKqr31QyfDjs3OnU\nDnyo2hyS77939goZONCn963RlCkQHu62eWtQ1/ZEhoWw0jrcTTNgicQ0Tmam84u6tsl+FTUCHzZv\nuZ1DkpEBffs6c1kCoV07uOgi+OCDakvLR4Q5CzjajommObBEYhpO1amR1NasBTBsmPPsw0Tidg5J\nRgYMHeqze3pk+nTYts3571TFiJQOZO4+bAs4mibPEolpuL174cCBupuOOnRw5pi4+WXqLdXmkBw4\nALt2BT6RXH65U1tz07yVZgs4mmbCEolpuIrEUFeNBJxkU9Ex7wPV5pBUjCY76yyf3dMjXbo4s/s/\n+KDaqRHdbQFH0zxYIjENt3Gj89y/f91lBw50yrsZCusN1eaQBEsiAbjiCmcEWXb2GYdjWkfQO74t\n6Tts5JZp2jxKJCLyvohcIiKWeMxpW7ZAmzbOX911GTjQWS5k+3afhFJtDklGhhNXZ7f7nvnXNNfG\noPOqLn7tbHS1KvugLeBomjRPE8PfgOuBLSLyhIj09WFMpqnIyoJevTxbnr2iH8VHzVvV5pB8/33g\n+0cq9OnjLNMyt/qOByO6x1JYXMqWfUUBCMwY7/Aokajq56r6Q2A4sAP4XESWiMiPXBtQmZYoK8vz\nzaIGDHCefZZIKs0hKSlx7hMMzVoVpk2Dr76qtnhlmqufxNbdMk2Zx01VItIRuAX4MfAd8CxOYnG/\n6YJp3kpLnWaqXr08Kx8dDd26+SSRVJtDsmEDnDwZPDUScBJJaSnMP2NXBbp3bE18dCTLt1kiMU2X\np30kHwDfAK2By1T1clV9W1XvBdr6MkATpHbtcn5Ze5pIwKmV+CCRVMwhSayokVR0tAdTIjnnHKe/\npkrzlogwqkdHlm4rQNX6SUzT5GmN5CVVHaCq/6uquQAiEgmgqmk+i84Er6ws57k+icRHI7cq5pAk\nVySSjAxnNrs/9mj3VEgIXHYZfPwxnDhzocbRPTuSf+QEW/P9t7e9Md7kaSL5HzfHfL/BhAleDUkk\nAwY4v0S9PHLL7RySwYOddbaCybRpzh4lVbYeHt2jIwBLt/p2LTJjfKXWRCIiXURkBNBKRIaJyHDX\nYzxOM5dpqbKynL/6ExI8v6ZivsmGDV4N5Yw5JKpOjSSYOtorXHQRtG5drXmre8fWJLSPYum2ggAF\nZkzj1FUjmQT8CUgCngb+7Hr8HPi1b0MzQa1i6G9IPaYWVSSSiomMXrLrQKU5JLt3O8ujBFP/SIVW\nrWDSJGc+SaX+EBFhdM+OLNt2wOaTmCap1t8Cqvq6ql4A3KKqF1R6XK6q7/spRhOMKhJJfcTEOJME\nfVAjSa4YsZWR4TwHY40EnOat3bur7c8yukdHDhwtYfO+IwEKzJiGq6tp6wbXyxQR+XnVhx/iM8FI\n1ennSE2t/7X9+/sgkRw/PYekYsTWkCFevYfXXHKJU4urMst9dM+KfhJr3jJNT13tEm1cz22BaDcP\n0xLt3w/Hj0P37vW/tl8/J5F4aahr8cky9h2pNIckIwN69nTmrQSjTp1gzJhq/SRJHVqTHNvKEolp\nksJqO6mqf3c9P+afcEyTULH4YEMSSf/+cPgw5OXVr6O+BrsPuUZsVR76G4z9I5VNmwa/+EW1Wt3o\nHh35ZF0epWXlp9cMM6YJ8HRC4pMi0k5EwkVkoYjkV2r2qu26ySKySUSyRORhN+dFRJ5znV8jIsMr\nnXtFRPaJyLoq1/xORHaLSIbrMdWT72C8qLGJBLzWvFUx9Dc5trUztHbr1qaRSKBa89b5feIoLC4l\nY5ftT2KaFk//7JmoqoXApThrbfUC/qu2C0QkFJgJTAEGANeJyIAqxaYAvV2P24EXKp17DZhcw8f/\nRVWHuh7zayhjfCWoEkmlDa3WrnWazIK1o71Cr17OnJoqzVvn9YojRODLTfkBCsyYhvE0kVQ0gV0C\n/EtVD3twzUggS1W3qWoJMAeYVqXMNOANdSwDYkQkAUBVvwZsAaJglJ3t9EHExNT/2q5dnWu9WCMJ\nDxXio6OCc2mUmkybBl9/7QxVdmnfOpzh3Trw5eZ9AQzMmPrzNJH8W0Q2AiOAhSISBxTXcU0isKvS\n+xzXsfqWcedeV1PYKyLSwYPyxpuys53aiCfLx1cl4tWRWzkHj5MY04rQEHH6Ryq29Q1206Y5S8VU\nWcRxfN841u0uZN+Run68jAkeni4j/zBwLpCmqieBo1SvXfjLC0APYCiQizNBshoRuV1E0kUkPT/f\nmgq8qiKRNFT//l6blLjrwLHTI7Yq9iBpSILzt7PPdgYbVGneGt83HoCvN9tyKabpqM/QkH7AtSJy\nE3AVMLGO8ruB5Ervk1zH6lvmDKq6V1XLVLUceAmnCc1duRdVNU1V0+Li4uoI1dSLNxLJnj3O6K1G\nOjWHpKwM1qwJ/v6RChWLOH7yyRmLOA5IaEentpF8ucmat0zT4emorTdxlkoZC5ztetS16u9KoLeI\npIpIBDADqLrX6DzgJtforVHA4YrVhWuJpfKY0enAuprKGh84cgQOHmx8IoFG10qKT5axv+iEk0i2\nbHHmtjSF/pEK06ZBURF88cWpQyEhwrg+cXyzZT+lZeUBDM4Yz9U6j6SSNGCA1mPDBFUtFZF7gAVA\nKPCKqmaKyB2u87OA+cBUIAs4Bvyo4noRmQ2MBzqJSA7wW1V9GXhSRIYCijOC7KeexmS8oGLEVrdu\nDf+MyiO3zjmnwR9zxtDf7xc7B5tSIrnwQmfP+7lzYcqUU4cv7h/Pe6tzWLHjAOf27BTAAI3xjKeJ\nZB3QBadPwmOuobnzqxybVem1AnfXcO11NRy/sT4xGC9rzNDfCqmpEBHR6A73XZWH/mZkQHj46STV\nFERFweTJ8NFH8Le/nVoAc1zfOCLDQliwLs8SiWkSPO0j6QSsF5EFIjKv4uHLwEyQ8kYiCQtzNp1q\nZCI5tQ9Jh9ZOR/uAAU6Cakouv9zpL6q0iGPriDDG9YljQeZeWw3YNAme1kh+58sgTBOSne38su7S\npXGf07//6ZV6Gyjn4DEiQkOIaxvpfNaECY2LKRAuucTZgGvuXGckl8uUwV34dP1evs85xLBuNsLd\nBDdPh/+pfYMqAAAgAElEQVR+hdMfEe56vRJY7cO4TLDKzobk5PrtQ+JO//6wbRsUN3y+xM6CYyTH\ntiJkfz7k5jat/pEKHTvC2LHVhgFf2K8zYSHCJ5l5AQrMGM95OmrrJ8C7wN9dhxKBD30VlAlijR36\nW6F/fygvd0ZbNdCOgmOkdGxzekZ7Uxn6W9W0abBunZNYXdq3CufcXp1YsC6PeoxxMSYgPP2z8m5g\nDFAIoKpbgHhfBWWCmDcTCTR4CLCqkl1wlO4d28B33zkHm2KNBE4v4lilVjJ5YBd2FBwjc09hAIIy\nxnOeJpITrvWyABCRMJzht6YlOXHCaULyRiLp29eZgd7ADvf8ohMcKymje8fWTiLp3h1iYxsfVyD0\n6AGDBlVLJFMGdSE8VJibUescXWMCztNE8pWI/BpoJSITgH8BH/kuLBOUdrmWRfNGImnVClJSGpxI\nsgucob+nEklTrY1UmDYNvvkGCk5vbNWhTQTj+8YzN2MPZTZ6ywQxTxPJw0A+sBZnAuB84L99FZQJ\nUt4Y+ltZIxZvrEgkqZEKmzfDsGHeiSlQrrjC6TOqUiuZPiyRfUdOsGSrrb1lgpeno7bKcTrX71LV\nq1T1pfrMcjfNhDdmtVfWvz9s2uSsk1XfUAqOEhoiJO7a4uxB0tQTyYgRzhbBs2efcfjCfvFER4bx\nwXfWvGWCV62JxLUG1u9EZD+wCdjk2h3xN/4JzwSVnTudfo3k5LrLeqJ/f2f4b0WCqocdBcdIjGlF\n2No1zoGmnkhE4PrrnXW3ck8vIBEVHsrUwQksWJfHsZLSAAZoTM3qqpE8gDNa62xVjVXVWOAcYIyI\nPODz6Exwyc52lj731uzxRuyW6IzYcvWPxMY2jT1I6nLddU7z1ttvn3F4+vBEjpaU8ck6m1NiglNd\nieRG4DpV3V5xQFW3ATcAN/kyMBOEvDX0t0K/fs5zAxLJjv1HnTkk333n1Eaawh4kdenf3/ku//d/\nZxw+JzWWlI6tmbNiVw0XGhNYdSWScFWt1sunqvlAuG9CMkHL24kkNhbi4+s9l+TQsRIKi0tJaR/u\n7NPe1Ju1Krv+eli58oyJmiLCtWd3Y8WOA2TtOxLA4Ixxr65EUtLAc6a5KS93hv96M5FAg0Zu7XCN\n2BpwaA+UlDSvRDJjhlO7qtLpftWIJMJCxGolJijVlUjOEpFCN48jwGB/BGiCRG4unDzpu0RSj0GA\nO/YfBSB11ybnQFOfQ1JZUhKcfz689dYZ/03ioiOZMKAz763O4URp/Ue5GeNLtSYSVQ1V1XZuHtGq\nak1bLYm355BU6N/f2XFxn+dby27NLyJEIG7DGmjb1pkl35zccIMzN2b58jMOzxjZjYPHTvJp5t4A\nBWaMe41cwtW0GL5MJFCv5q2sfUV079iG0PSVztLroaHejSnQrrkGWreGl18+4/B5vTqRGNOKOSt3\nBigwY9yzRGI84+3JiBUamEj6tg93Vv0dOdK78QSDdu3g2mthzhxnT3eXkBBhxtnJfJtVwLb8olo+\nwBj/skRiPLNzJ3ToANHR3v3cxETnMz1MJCfLytlRcJQxRTudPptG7Pke1H78YyeJvPPOGYevHZlM\neKjwz2VWKzHBwxKJ8Yy3h/5WEHHmk3iYSLILjnGyTBmye7NzoDnWSABGj3b+u/zjH2ccjo+OYvKg\nBP61apfNdDdBwxKJ8YyvEgnAkCHOxEIPRm5l7XOadLplrXVqM4mJvokp0EScWsnSpbB+/Rmnbhrd\nnSPFpczN2BOg4Iw5kyUSUzdV3yaSs892lk/fsaPOoltdfQPt13zXfGsjFW68EcLD4e9/P+NwWvcO\n9OsSzZtLs233RBMULJGYuh086LTX+zKRAKxYUWfRrH1FDJYiQrZvgzFjfBNPsIiPd0ZwvfoqHDk9\no11EuHF0d9bnFrJ658EABmiMwxKJqdtOV8eurxLJ4MEQGeksDVKHrH1FTD7gWj5k3DjfxBNM7rvP\nSSKvv37G4SuGJhIdGcYbS+u/crIx3ubTRCIik0Vkk4hkicjDbs6LiDznOr9GRIZXOveKiOwTkXVV\nrokVkc9EZIvruYMvv4PBd0N/K4SHO8uc1JFISsvK2bz3COfsWueM9GpOM9prMnKkMzLt+eedZWpc\n2kSGceWIJOavzSX/yIkABmiMDxOJiIQCM4EpwADgOhEZUKXYFKC363E78EKlc68Bk9189MPAQlXt\nDSx0vTe+5KvJiJWdfTasWgWlNY9E2lFwlBOl5fTZuNpp1goL8108weS++5yZ7p9+esbhG0d352SZ\n8k66rb9lAsuXNZKRQJaqblPVEmAOMK1KmWnAG+pYBsSISAKAqn4NHHDzudOAinr+68AVPonenJad\n7eyxHhfnu3uMHQtHjzrJpAaZewrpePQQ7bZvaRnNWhWuugq6dIHnnjvjcM+4tozp1ZG3lmVTWlZe\nw8XG+J4vE0kiUPlPpRzXsfqWqaqzqlZsIZcHdG5MkMYDO3c6zVq+3PNj/HjnedGiGouszy3kouzv\nnDcTJvgulmATEQF33gkff1xtyf0bR6Ww53AxCzd6vlaZMd7WpDvbXfvGux3/KCK3i0i6iKTn5+f7\nObJmxpdDfyvEx8OgQc5WszVYv6eQS3NWO3+dN6el4z1xxx0QFQVPPXXG4Yv7x5PQPop/LrNOdxM4\nvkwku4HKm3snuY7Vt0xVeyuav1zPbv8UU9UXVTVNVdPifNkk0xL4I5EAXHghLF4MJ6p3Hqsqm3cd\n4OxNK2HqVAhp0n8D1V98PNx2G7z5JuTknDocFhrC9SO78c2W/bb+lgkYX/40rgR6i0iqiEQAM4B5\nVcrMA25yjd4aBRyu1GxVk3nAza7XNwNzvRm0qeL4cWeJd38lkuPH4dtvq53aW3iC1M0ZtDpWBJdc\n4vtYgtEvfuGM3Hr66TMOzxjZjfBQ4U2rlZgA8VkiUdVS4B5gAbABeEdVM0XkDhG5w1VsPrANyAJe\nAu6quF5EZgNLgb4ikiMit7lOPQFMEJEtwMWu98ZXdrm6sHw19Leyiy92OvXfe6/aqdU7D3L5+q8p\na9UaJk70fSzBKCUFrrsOXnzRWQnAJS46kimDEnh3VY6tv2UCwqftA6o6X1X7qGpPVX3cdWyWqs5y\nvVZVvdt1frCqple69jpVTVDVcFVNUtWXXccLVPUiVe2tqherqruRXcZb/DH0t0KbNk6z1XvvQdmZ\nuwBmbMnjkk2LYfp0ZzOrluqhh5zRbX/96xmHbf0tE0gtrKHZ1Js/Ewk4e5bv3Qvz559xuPWH79G+\nuIjQW252f11LMWgQXHaZMxS40l4lI1zrb71h62+ZALBEYmqXne3sQOivVXanTXP2La/UD1B84iSX\nzX+dval9neavlu5Xv4IDB+CF0/N3RYSbRqewIbeQVdm2/pbxL0skpnbZ2U4S8dcs8vBweOAB+PJL\n+M9/ANj7x7/QsyCHvLse8O1clqZi9Ginn+jJJ8+olVwxrCvRUWG8tmRH4GIzLZIlElO7HTucTl5/\nuucepwnn+uvh9ttJ+n+P8GWPEXS74xb/xhHMHnsM9u8/o6+kdUQYM85O5uN1eeQePh7A4ExLY4nE\n1G7bNujZ07/3jIhw+kgGD4aXXyZ90LnMuuNxOrSN9G8cwWzUKGdgwlNPQWHhqcM3jU5BVW1VYONX\nlkhMzY4fh927oUcP/987ORkWL6aw6DjXT32ItCEp/o8h2D32mNNX8vzzpw4lx7Zm0sAu/N/ynRwv\nKavlYmO8xxKJqVnFjoWBSCQuX27eT1m5Mq6vrU5QTVqaM4LrT3+Cw4dPHb51bCqHj5/k/e9yarnY\nGO+xRGJqtm2b8xzARDJ/TS7x0ZGM6Gbbzrj12GNw6BA888ypQ2ndOzA4sT2vLN5OebkNBTa+Z4nE\n1GzrVufZ330kLkUnSlm0aR9TBycQEmKjtdwaNsyZpPn006dmu4sIt45NYWv+Ub7J2h/gAE1LYInE\n1GzbNmcWeadOAbn9+6tzOFFazhXD/DSHpan6/e+dYcCPP37q0CWDuxIXHcnLi7cHMDDTUlgiMTXb\nts1p1vLS3I3dh46zcscB8g4X11m2vFx5fckOzkpqz9DkGK/cv9kaOBBuuQVmzjzVrxURFsIt56bw\n9eZ81u0+XOvlxjSWJRJTs61bvdI/sq+wmFtfW8mYJ77g6llLGfW/C7nuxWW1zsB+/7vdbM0/ym3n\nBa5/pkl57DFnaf1HHz116MbR3YmOCuOvX2QFMDDTElgiMe6pemUOyb7CYq6ctYRl2wp4cEIfXr91\nJP81qS9b9hVx5QtL+Pk7GeQfOXP/kbzDxTzx8UbOSo7h0sEJjbp/i5GUBPffD2+9BRkZALSLCudH\n56bwSWYem/KOBDhA05xZIjHu5eVBcXGjaiRl5co9s79j/5ES/u8no7j3ot6M6xPH3Rf04qv/Gs9d\n43vy0fd7uPDPX/Ly4u3sPnScjF2HuPmVFRwvKeXJK4dYJ3t9PPwwxMQ4zy4/GpNK64hQZi6yWonx\nHUskxj0vDP19e+UuVmw/wO+vGFStn6NNZBi/nNyPj+8/nyFJ7fn9v9cz5okvuGLmt+QVFvPiTWn0\n7RLdmG/Q8sTEwCOPwIIFsHAhAB3aRHDjqO78e80ettoOisZHpCUsOZ2Wlqbp6el1FzSnvfEG3Hwz\nbNoEffrU+/LC4pNc8NSX9Ixry9s/HYXU0mGvqqzbXcj3OYdoFR7KxQM6075VeGOib7mKi6FfP2jf\nHlatgrAw8o+cYNxTizi/dxyzbhwR6AhNEyIiq1Q1ra5yViMx7m3d6ozWauA+JG8s2UHB0RIevXRA\nrUkEnHkPg5Pac8Oo7lw5IsmSSGNERTkz3descXZSxNlB8afn9+STzDxWZds+cMb7LJEY9zZvdlb9\njaz/QoknSst4bUk24/rEMTipvfdjM7W78kq44AL47/8+NUnxJ+enEhcdyeP/2WAbXxmvs0Ri3Nu4\n0WkiaYC5GXvYX3SCn9jQ3cAQgWefdVYF/s1vAGeJ+Z9P6MPqnYf4z9rcAAdomhtLJKa68nKnRtK3\nb4Mun7NiJ73j2zKmV0cvB2Y8Nngw3HUXzJoF338PwNUjkhjYtR3/76P1FBafDHCApjmxRGKq270b\njh1rUI1ka34Rq3ce4uq0pDr7RoyPPfYYxMbCnXdCWRlhoSH87w8Gs7/oBE99sinQ0ZlmxBKJqW7j\nRue5ATWS91blEBoiXDHU1scKuA4dnMUcly6Fv/0NgCFJMdw0OoV/Ls+2vd2N11giMdVtcv21Ws8a\nSVm58v7q3YzrE0d8uygfBGbq7YYbYNIk+NWvINvZNfEXk/qS0C6KB9/J4OiJ0gAHaJoDSySmuo0b\noV076Ny5Xpet3nmQvMJiW603mIjA3//uvL7jDlClbWQYT187lOwDx3jso8zAxmeaBZ8mEhGZLCKb\nRCRLRB52c15E5DnX+TUiMryua0XkdyKyW0QyXI+pvvwOLdKmTU6zVj37OD5Zl0dEaAgX9ov3UWCm\nQbp3hz/8AT75BF59FYBRPTpy1/ievJOew3wbxWUayWeJRERCgZnAFGAAcJ2IDKhSbArQ2/W4HXjB\nw2v/oqpDXY/5vvoOLdamTfVu1lJVPlmXx3m9O9E2MsxHgZkGu/tuGDcO7rsPtmwB4GcX9+GspPY8\n9N4asguOBjhA05T5skYyEshS1W2qWgLMAaZVKTMNeEMdy4AYEUnw8FrjC0VFsGtXvTvaM/cUsvvQ\ncSYN7OKjwEyjhIbCm29CRARcfz2UlBAeGsJfrx9OiAg/fXMVx0qsv8Q0jC8TSSKwq9L7HNcxT8rU\nde29rqawV0TENvP2pkxXm/ngwfW67JN1eYQIXDygfv0qxo+Sk+GllyA9/dS+JcmxrXnuumFs2nuE\nX72/1ma9mwZpip3tLwA9gKFALvBnd4VE5HYRSReR9Pz8fH/G17StWeM81zORLMjM45zUjsS2ifBB\nUMZrrrwSfvpTePJJeP99AMb1ieMXE/syN2MPr367I7DxmSbJl4lkN5Bc6X2S65gnZWq8VlX3qmqZ\nqpYDL+E0g1Wjqi+qapqqpsXFxTXqi7Qoa9ZAdHS9FmvM2lfEln1FTB5kzVpNwrPPwjnnwE03naqB\n3jmuJxMHdObx+RtYtq0gwAGapsaXiWQl0FtEUkUkApgBzKtSZh5wk2v01ijgsKrm1natqw+lwnRg\nnQ+/Q8uzZg0MGuRs2+qhBZl5AEwcaM1aTUJkJLz3nvMHw/TpcOgQISHCn685i+4dW3PP/60m73Bx\noKM0TYjPEomqlgL3AAuADcA7qpopIneIyB2uYvOBbUAWTu3irtqudV3zpIisFZE1wAXAA776Di2O\nKqxdC0OG1OuyBZl5DE2OIaF9Kx8FZrwuMRHefRe2b4err4aSEqKjwvn7DSM4VlLGXW+toqS0PNBR\nmibCp30kqjpfVfuoak9Vfdx1bJaqznK9VlW923V+sKqm13at6/iNrrJDVPVyVw3GeMPu3XDwYL0S\nye5Dx1mTc9iatZqiMWPgH/+Azz+HH/8YVOndOZqnrjqL1TsP8T//WR/oCE0TYQP+zWkZGc5zPRLJ\ngnVOs5YN+22ibr7ZGe796KPOqK7HH+eSIQlk7ErlpW+2MzQ5hh8MTwp0lCbIWSIxp61c6fSNDBvm\n8SULMvPo2zma1E5tfBiY8alHHnGSyR/+AElJcOedPDS5H2tyDvPrD9bSr0s7BnRtF+goTRBrisN/\nja+sXAkDB0Ibz5JCQdEJVu44wCRr1mraRGDmTLj0UmcG/OzZhLkmK8a0iuCOf67i8DHbv8TUzBKJ\ncajCihVw9tkeX/L5hr2UK0yy0VpNX1gYvPOOs4zKjTfC3LnERUcy84fDyT18nJ+9/R3l5TZZ0bhn\nicQ4duxw9vce6XZajlufrMsjObYVAxKs2aNZaNUK5s2DtDS45hr49FNGdO/Aby4dwKJN+Tz/RVag\nIzRByhKJcaxY4Tx7mEiOFJ/k26wCJg/sYjshNifR0fDxx9C/P1xxBXzzDTeM6s4PhiXyzMLNfLlp\nX6AjNEHIEolxLFni/EU6aJBHxRdtyqekrNxGazVHHTrAp59Ct25wySVIejqPTx9M387R3D8ng10H\njgU6QhNkLJEYx1dfOfMKwsM9Kr5gXR5x0ZEM72ZrZjZL8fGwcCF06gQTJ9JqzXf8/cYRlKty51ur\nKD5ZFugITRCxRGLgwAFnaZTx4z0qXnSilIUb9zJ5YBdCQqxZq9lKTIRFiyAmBi6+mO5bM/nLNUNZ\nt7uQ3861nRXNaZZIDHz9tTNqy8NEsmBdHsUny7liWFffxmUCr3t3p7YaGwsTJnDxkR3cc0Ev3k7f\nxZwVOwMdnQkSlkgMfPml0z/i4dDfDzN2kxzbypq1Wopu3Zxk0qkTTJjAA20LOK93J34zN5OVOw4E\nOjoTBCyRGGcv77Fjnd3z6rCvsJhvs/YzfWiijdZqSZKTnWTSuTOhkycxs9sxkjq04idvpJO1ryjQ\n0ZkAs0TS0m3a5Dwuv9yj4vO+30O5wrRhVTe7NM1eUpJTe+3alXbTLuHthP2EhQi3vLqCfUds2fmW\nzBJJS/fRR87zZZfVWVRVmb1iJ0OTY+gZ19bHgZmglJgI33wDAwcSd/MMPojcSEFRCT98aTn5R04E\nOjoTIJZIWrp58+CsszzaEXHp1gK25h/lxlGe755omqH4eGc018UXk/zL+/j00EJ2HzjKjBeXsrfQ\naiYtkSWSlmznTli82JnB7IE3lmbToXU4lwxJqLuwad7atnVqsz/6Ecl//RPfLnueor0F/OBvS8jc\nczjQ0Rk/s0TSkr3xhjPs9+ab6yy668AxPtuwl2vOTiYqPNQPwZmgFx4OL78Mzz9Ph68+5+t3H6Ln\n7i1c9cJS5n2/J9DRGT+yRNJSqcJrr8EFF0Bqap3FX/hqK6Ei3HJuis9DM02ICNxzDyxaROTxo7z+\n93t59Pv3eeCtdO6b/R2HjpUEOkLjB5ZIWqpPPoGtW+HWW+ssuvvQcf6Vvotrzk6yfdmNe2PHwrp1\nyJVXcv1HL7H0g1+xd/7nXPTnr3hz6Q5Oltn+782ZJZKWSBX+53+ciWbXXltn8T9+vJEQEe4c38sP\nwZkmq2NHmD0b/vUv4o8f5u1//pLn5/2RV19ZwORnvuaD73IsoTRTttVuS/T5585qvzNn1rlI4/Jt\nBcz7fg/3XdSbxBirjRgPXHUVTJkCTz7J6CefZOGqL/jqrPE8te4K/tR3ELeNTeWqtCTaRXm2QKgJ\nfqLa/Hc9S0tL0/T09ECHERyKi53hvqWlkJkJUVE1Fj1SfJKpz32DKnz2wDhaRVgnu6mnvXvhmWfQ\nmTORI0fYnDKQl/tdyOdDLuCis3tw3chuDE2OsVUSgpSIrFLVtDrLWSJpYX75S3jqKViwACZOrLFY\nebly75zv+HhtLv+6YzQjusf6MUjT7Bw6BK++Ci+9BBs2UBIZxaIeafy712h2n3sBF4/qzaWDu9Kt\nY+tAR2oqsURSiSUSlzffhJtugjvugBdeqLGYqvLYR+t5bckOfjWlHz8d19OPQZpmTRWWLoU336T8\ngw8I2buXk2HhLEkezJLuQ9iXNobek87j3L6dGZzYnlDbpiCggiKRiMhk4FkgFPiHqj5R5by4zk8F\njgG3qOrq2q4VkVjgbSAF2AFco6oHa4ujxScSVSdx3HsvnH++s/tdDX0jR0+U8qv31zLv+z3cOiaV\nRy/tb80OxjfKypyk8t57nPx4AeGbNgBwOLINqxL7szmpLyeHD6ftmFH0GNKbwYntiW1T98KixnsC\nnkhEJBTYDEwAcoCVwHWqur5SmanAvTiJ5BzgWVU9p7ZrReRJ4ICqPiEiDwMdVPWh2mJp0YlkzRr4\n9a/hP/+BSy6Bt9+GNm2qFSspLWfe93v4y2eb2XP4OL+Y2Je7xve0JGL8JzcXFi2ieMFnlHy7lLbb\nthCiziiv/DYxbItNIq9LN46m9kb79qHtgL50HtiLHt07Excdaf9WfcDTROLLUVsjgSxV3eYKaA4w\nDVhfqcw04A11stkyEYkRkQSc2kZN104Dxruufx34Eqg1kbQIZWVw8CDk58P69ZCRAR9/DKtWOctZ\n/OUvTo0kNJTik2UcPFbCnkPFbMo7wuqdB1m0cR8FR0vo1yWaf80YTVqK9YkYP0tIgOuvJ+r664kC\nOHoUMjI4vnQ55cvS6bZpMwM3LqHtyvlnXHYoqi1b2nXiYMcuHI9PQOPjCY/vRKvOcUR0iScqvhOt\nE+KJ7hhDm5hoQtu0hshIZzKl8QpfJpJEYFel9zk4tY66yiTWcW1nVc11vc4DOnsr4Gp+/3uYPZv8\nI8UcP1EKgIDTVAQIrmdVnJdKxT9NUYVK5+VUxa/SMddrtMpnVXo+4xha6d6nj4doOW1PHDsj9HIR\n1iX25fMJP+b94ZMpPNQWfv85JWXlFJ88cyx/u6gwxveNZ/rwRMb3ibO/7ExwaNMGxoyh1ZgxnDHw\n/MABdONGDmVu4vCmbRRvz4acHLrs3UNM+gaijxUSUkdLSznCibAISsIjKAsNRSWEcpHqzyEhKIKK\n86jMmz8lVT/bW9pEhNLprdfgvPN88vkVmvQ8ElVVEXH7L0ZEbgduB+jWrVvDbtC1KwwaxJH8oxw8\nfvL0/+xK/9NVXHM65fRxFUFwpZVK11S7Xpz0IZWvkerXnDpecaNK9wEgRChu3Zbj0R041i6GgsQU\n8pN7UBrp/PhdXCne8FAhpnUEsW0iiI+OpE/naJI6tLLkYZqO2Fjk3HPpcO65uN2js6yMo3v3s39n\nLsdy91Kcl8/J/AJKCo9Qduw4cvw4UnyckOJiQoqL0bIypLwctBwpL0e0HCnXM54p9+FESh/2U5e3\niaBTdLTPPr+CLxPJbiC50vsk1zFPyoTXcu1eEUlQ1VxXM9g+dzdX1ReBF8HpI2nQN7jtNrjtNno0\n6GJjTECEhtKma2fadPVdY4U5ky+XSFkJ9BaRVBGJAGYA86qUmQfcJI5RwGFXs1Vt184DKparvRmY\n68PvYIwxpg4+q5GoaqmI3AMswBnC+4qqZorIHa7zs4D5OCO2snCG//6otmtdH/0E8I6I3AZkA9f4\n6jsYY4ypm01INMYY45anw39t9V9jjDGNYonEGGNMo1giMcYY0yiWSIwxxjSKJRJjjDGN0iJGbYlI\nPs5Q4YboBOz3YjjeYnHVj8VVPxZX/QRrXNC42LqralxdhVpEImkMEUn3ZPibv1lc9WNx1Y/FVT/B\nGhf4JzZr2jLGGNMolkiMMcY0iiWSur0Y6ABqYHHVj8VVPxZX/QRrXOCH2KyPxBhjTKNYjcQYY0yj\nWCKpg4g8KCIqIp0qHfuViGSJyCYRmeTneJ4SkY0iskZEPhCRmGCIy3X/ya57Z4nIw/6+f6U4kkVk\nkYisF5FMEbnfdTxWRD4TkS2uZ7f7Ivk4tlAR+U5E/h0sMbniiBGRd13/tjaIyOhgiE1EHnD9P1wn\nIrNFJCoQcYnIKyKyT0TWVTpWYxz++lmsIS6//46wRFILEUkGJgI7Kx0bgLM/ykBgMvA3EQn1Y1if\nAYNUdQiwGfhVMMTlutdMYAowALjOFVMglAIPquoAYBRwtyuWh4GFqtobWOh672/3AxsqvQ+GmACe\nBT5R1X7AWTgxBjQ2EUkE7gPSVHUQzpYSMwIU12s4P1eVuY3Dzz+L7uLy++8ISyS1+wvwSyo2WndM\nA+ao6glV3Y6zl8pIfwWkqp+qaqnr7TKc3SMDHpfrXlmquk1VS4A5rpj8TlVzVXW16/URnF+Kia54\nXncVex24wp9xiUgScAnwj0qHAxqTK672wPnAywCqWqKqh4IhNpw9k1qJSBjQGtgTiLhU9WvgQJXD\nNcXht59Fd3EF4neEJZIaiMg0YLeqfl/lVCKwq9L7HNexQLgV+Nj1OtBxBfr+bolICjAMWA50du3A\nCZAH+Hsv1mdw/jCpvAF4oGMCSAXygVddzW7/EJE2gY5NVXcDf8JpEcjF2UH100DHVUlNcQTTz4Jf\nfuNjJ8AAAAJGSURBVEf4cs/2oCcinwNd3Jx6BPg1TrOW39UWl6rOdZV5BKcJ5y1/xtaUiEhb4D3g\nZ6paKCKnzqmqiojfhiyKyKXAPlVdJSLj3ZXxd0yVhAHDgXtVdbmIPEuV5qJAxObqc5iGk+gOAf8S\nkRsCHZc7wRJHZf78HdGiE4mqXuzuuIgMxvnH+73rl08SsFpERgK7geRKxZNcx3weV6X4bgEuBS7S\n0+O3fR5XHQJ9/zOISDhOEnlLVd93Hd4rIgmqmisiCcA+P4Y0BrhcRKYCUUA7EflngGOqkAPkqOpy\n1/t3cRJJoGO7GNiuqvkAIvI+cG4QxFWhpjgC/rPg798R1rTlhqquVdV4VU1R1RScH7ThqpoHzANm\niEikiKQCvYEV/opNRCbjNI9crqrHKp0KaFzASqC3iKSKSAROp948P97/FHGy/8vABlV9utKpecDN\nrtc3A3P9FZOq/kpVk1z/nmYAX6jqDYGMqVJsecAuEenrOnQRsD4IYtsJjBKR1q7/pxfh9HcFOq4K\nNcXR8n5HqKo96ngAO4BOld4/AmwFNgFT/BxLFk47Z4brMSsY4nLdfyrOKJGtOM1wgfr/NRZngMSa\nSv+dpgIdcUbXbAE+B2IDFN944N+u18ES01Ag3fXf7EOgQzDEBjwGbATWAW8CkYGIC5iN009zEucP\ny9tqi8NfP4s1xOX33xE2s90YY0yjWNOWMcaYRrFEYowxplEskRhjjGkUSyTGGGMaxRKJMcaYRrFE\nYowxplEskRhjjGkUSyTGGGMa5f8DGkd3aOS81wwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xc46f787358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# filling NA with zeroes creates a peak of population around 0, as expected\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "X_train['Age'].plot(kind='kde', ax=ax)\n",
    "X_train.Age_zero.plot(kind='kde', ax=ax, color='red')\n",
    "lines, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(lines, labels, loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filling NA with 0s also distorts the distribution of the original variable, generating an accumulation of values around 0. We will see in the next lecture a method of NA imputation that preserves variable distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Machine learning model performance on different imputation methods\n",
    "\n",
    "#### Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train set zero imputation\n",
      "Logistic Regression roc-auc: 0.6863462831608859\n",
      "Test set zero imputation\n",
      "Logistic Regression roc-auc: 0.7137499999999999\n",
      "\n",
      "Train set median imputation\n",
      "Logistic Regression roc-auc: 0.6795188824537428\n",
      "Test set median imputation\n",
      "Logistic Regression roc-auc: 0.7244940476190476\n"
     ]
    }
   ],
   "source": [
    "# Let's compare the performance of Logistic Regression using Age filled with zeros or alternatively the median\n",
    "\n",
    "# model on NA imputed with zeroes\n",
    "logit = LogisticRegression(random_state=44, C=1000) # c big to avoid regularization\n",
    "logit.fit(X_train[['Age_zero', 'Fare']], y_train)\n",
    "print('Train set zero imputation')\n",
    "pred = logit.predict_proba(X_train[['Age_zero', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set zero imputation')\n",
    "pred = logit.predict_proba(X_test[['Age_zero', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))\n",
    "print()\n",
    "\n",
    "# model on NA imputed with median\n",
    "logit = LogisticRegression(random_state=44, C=1000) # c big to avoid regularization\n",
    "logit.fit(X_train[['Age_median', 'Fare']], y_train)\n",
    "print('Train set median imputation')\n",
    "pred = logit.predict_proba(X_train[['Age_median', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set median imputation')\n",
    "pred = logit.predict_proba(X_test[['Age_median', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that median imputation leads to better performance of the logistic regression. Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average total survival: 0.3884430176565008\n"
     ]
    }
   ],
   "source": [
    "print('Average total survival:', X_train.Survived.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average real survival of children:  0.5740740740740741\n",
      "Average survival of children when using Age imputed with zeroes:  0.38857142857142857\n",
      "Average survival of children when using Age imputed with median:  0.5740740740740741\n"
     ]
    }
   ],
   "source": [
    "print('Average real survival of children: ', X_train[X_train.Age<15].Survived.mean())\n",
    "print('Average survival of children when using Age imputed with zeroes: ', X_train[X_train.Age_zero<15].Survived.mean())\n",
    "print('Average survival of children when using Age imputed with median: ', X_train[X_train.Age_median<15].Survived.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Children were more likely to survive the catastrophe (0.57 for children vs 0.38 for the entire Titanic). Thus, smaller values of Age are a good indicator of survival.\n",
    "\n",
    "When we replace NA with zeroes, we are masking the predictive power of Age. After zero imputation it looks like children did not have a greater chance of survival, and therefore the model loses predictive power.\n",
    "\n",
    "On the other hand, replacing NA with the median, preserves the predictive power of the variable Age, as smaller Age values will favour survival."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Support Vector Machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train set zero imputation\n",
      "Logistic Regression roc-auc: 0.687935185787727\n",
      "Test set zero imputation\n",
      "Logistic Regression roc-auc: 0.7243452380952381\n",
      "\n",
      "Train set median imputation\n",
      "Logistic Regression roc-auc: 0.692989306088805\n",
      "Test set median imputation\n",
      "Logistic Regression roc-auc: 0.7417559523809523\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Let's compare the performance of SVM using Age filled with zeros or alternatively the median\n",
    "\n",
    "SVM_model = SVC(random_state=44, probability=True, max_iter=-1, kernel='linear',)\n",
    "SVM_model.fit(X_train[['Age_zero', 'Fare']], y_train)\n",
    "print('Train set zero imputation')\n",
    "pred = SVM_model.predict_proba(X_train[['Age_zero', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set zero imputation')\n",
    "pred = SVM_model.predict_proba(X_test[['Age_zero', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))\n",
    "print()\n",
    "SVM_model = SVC(random_state=44, probability=True,  max_iter=-1, kernel='linear')\n",
    "SVM_model.fit(X_train[['Age_median', 'Fare']], y_train)\n",
    "print('Train set median imputation')\n",
    "pred = SVM_model.predict_proba(X_train[['Age_median', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set median imputation')\n",
    "pred = SVM_model.predict_proba(X_test[['Age_median', 'Fare']])\n",
    "print('Logistic Regression roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For SVM as well, median imputation leads to better performance of the model, compared to replacing NA with zeroes.\n",
    "\n",
    "#### Random Forests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train set zero imputation\n",
      "Random Forests roc-auc: 0.7555855621353116\n",
      "Test set zero imputation\n",
      "Random Forests zero imputation roc-auc: 0.7490476190476191\n",
      "\n",
      "Train set median imputation\n",
      "Random Forests roc-auc: 0.7515021366130887\n",
      "Test set median imputation\n",
      "Random Forests roc-auc: 0.7697619047619048\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Let's compare the performance of Random Forests using Age filled with zeros or alternatively the median\n",
    "\n",
    "rf = RandomForestClassifier(n_estimators=100, random_state=39, max_depth=3)\n",
    "rf.fit(X_train[['Age_zero', 'Fare']], y_train)\n",
    "print('Train set zero imputation')\n",
    "pred = rf.predict_proba(X_train[['Age_zero', 'Fare']])\n",
    "print('Random Forests roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set zero imputation')\n",
    "pred = rf.predict_proba(X_test[['Age_zero', 'Fare']])\n",
    "print('Random Forests zero imputation roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))\n",
    "print()\n",
    "rf = RandomForestClassifier(n_estimators=100, random_state=39, max_depth=3)\n",
    "rf.fit(X_train[['Age_median', 'Fare']], y_train)\n",
    "print('Train set median imputation')\n",
    "pred = rf.predict_proba(X_train[['Age_median', 'Fare']])\n",
    "print('Random Forests roc-auc: {}'.format(roc_auc_score(y_train, pred[:,1])))\n",
    "print('Test set median imputation')\n",
    "pred = rf.predict_proba(X_test[['Age_median', 'Fare']])\n",
    "print('Random Forests roc-auc: {}'.format(roc_auc_score(y_test, pred[:,1])))\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random Forests, as well as SVMs, perform better with median imputation, compared with replacing with zeroes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## House Sale Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# we are going to train a model on the following variables,\n",
    "# they contain the 3 numerical variables with NA\n",
    "\n",
    "cols_to_use = ['OverallQual', 'TotalBsmtSF', '1stFlrSF', 'GrLivArea','WoodDeckSF', 'BsmtUnfSF',\n",
    "               'LotFrontage', 'MasVnrArea', 'GarageYrBlt', 'SalePrice']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1460, 10)\n"
     ]
    },
    {
     "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>LotFrontage</th>\n",
       "      <th>OverallQual</th>\n",
       "      <th>MasVnrArea</th>\n",
       "      <th>BsmtUnfSF</th>\n",
       "      <th>TotalBsmtSF</th>\n",
       "      <th>1stFlrSF</th>\n",
       "      <th>GrLivArea</th>\n",
       "      <th>GarageYrBlt</th>\n",
       "      <th>WoodDeckSF</th>\n",
       "      <th>SalePrice</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>65.0</td>\n",
       "      <td>7</td>\n",
       "      <td>196.0</td>\n",
       "      <td>150</td>\n",
       "      <td>856</td>\n",
       "      <td>856</td>\n",
       "      <td>1710</td>\n",
       "      <td>2003.0</td>\n",
       "      <td>0</td>\n",
       "      <td>208500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>80.0</td>\n",
       "      <td>6</td>\n",
       "      <td>0.0</td>\n",
       "      <td>284</td>\n",
       "      <td>1262</td>\n",
       "      <td>1262</td>\n",
       "      <td>1262</td>\n",
       "      <td>1976.0</td>\n",
       "      <td>298</td>\n",
       "      <td>181500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>68.0</td>\n",
       "      <td>7</td>\n",
       "      <td>162.0</td>\n",
       "      <td>434</td>\n",
       "      <td>920</td>\n",
       "      <td>920</td>\n",
       "      <td>1786</td>\n",
       "      <td>2001.0</td>\n",
       "      <td>0</td>\n",
       "      <td>223500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>60.0</td>\n",
       "      <td>7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>540</td>\n",
       "      <td>756</td>\n",
       "      <td>961</td>\n",
       "      <td>1717</td>\n",
       "      <td>1998.0</td>\n",
       "      <td>0</td>\n",
       "      <td>140000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>84.0</td>\n",
       "      <td>8</td>\n",
       "      <td>350.0</td>\n",
       "      <td>490</td>\n",
       "      <td>1145</td>\n",
       "      <td>1145</td>\n",
       "      <td>2198</td>\n",
       "      <td>2000.0</td>\n",
       "      <td>192</td>\n",
       "      <td>250000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   LotFrontage  OverallQual  MasVnrArea  BsmtUnfSF  TotalBsmtSF  1stFlrSF  \\\n",
       "0         65.0            7       196.0        150          856       856   \n",
       "1         80.0            6         0.0        284         1262      1262   \n",
       "2         68.0            7       162.0        434          920       920   \n",
       "3         60.0            7         0.0        540          756       961   \n",
       "4         84.0            8       350.0        490         1145      1145   \n",
       "\n",
       "   GrLivArea  GarageYrBlt  WoodDeckSF  SalePrice  \n",
       "0       1710       2003.0           0     208500  \n",
       "1       1262       1976.0         298     181500  \n",
       "2       1786       2001.0           0     223500  \n",
       "3       1717       1998.0           0     140000  \n",
       "4       2198       2000.0         192     250000  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's load the House Sale Price dataset\n",
    "\n",
    "data = pd.read_csv('houseprice.csv', usecols=cols_to_use)\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OverallQual    0.000000\n",
       "TotalBsmtSF    0.000000\n",
       "1stFlrSF       0.000000\n",
       "GrLivArea      0.000000\n",
       "WoodDeckSF     0.000000\n",
       "BsmtUnfSF      0.000000\n",
       "LotFrontage    0.177397\n",
       "MasVnrArea     0.005479\n",
       "GarageYrBlt    0.055479\n",
       "SalePrice      0.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# let's plot the histograms to have an impression of the distribution of the numerical variables\n",
    "data[cols_to_use].isnull().mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 4 variables with missing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1022, 10), (438, 10))"
      ]
     },
     "execution_count": 20,
     "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, data.SalePrice, test_size=0.3,\n",
    "                                                    random_state=0)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's impute the NA with the median or comparatively, zeros\n",
    "# remember that we need to impute with the median for the train set, and then propagate to test set\n",
    "# we use the function that we made a few cells ago.\n",
    "\n",
    "median = X_train.LotFrontage.median()\n",
    "impute_na(X_train, 'LotFrontage', median)\n",
    "impute_na(X_test, 'LotFrontage', median)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "median = X_train.MasVnrArea.median()\n",
    "impute_na(X_train, 'MasVnrArea', median)\n",
    "impute_na(X_test, 'MasVnrArea', median)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "median = X_train.GarageYrBlt.median()\n",
    "impute_na(X_train, 'GarageYrBlt', median)\n",
    "impute_na(X_test, 'GarageYrBlt', median)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['OverallQual',\n",
       " 'TotalBsmtSF',\n",
       " '1stFlrSF',\n",
       " 'GrLivArea',\n",
       " 'WoodDeckSF',\n",
       " 'BsmtUnfSF',\n",
       " 'LotFrontage',\n",
       " 'MasVnrArea',\n",
       " 'GarageYrBlt']"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a list with the untransformed columns\n",
    "cols_to_use.remove('SalePrice')\n",
    "cols_to_use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['OverallQual',\n",
       " 'TotalBsmtSF',\n",
       " '1stFlrSF',\n",
       " 'GrLivArea',\n",
       " 'WoodDeckSF',\n",
       " 'BsmtUnfSF',\n",
       " 'LotFrontage_median',\n",
       " 'MasVnrArea_median',\n",
       " 'GarageYrBlt_median']"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create a list with the transformed columns\n",
    "cols_zero = [col+'_zero' if col in ['LotFrontage', 'MasVnrArea', 'GarageYrBlt'] else col for col in cols_to_use]\n",
    "cols_median = [col+'_median' if col in ['LotFrontage', 'MasVnrArea', 'GarageYrBlt'] else col for col in cols_to_use]\n",
    "cols_median"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Machine learning model performance\n",
    "#### Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train set zero imputation\n",
      "Linear Regression mse: 1232243282.8714073\n",
      "Test set zero imputation\n",
      "Linear Regression mse: 2162158044.601015\n",
      "\n",
      "Train set median imputation\n",
      "Linear Regression mse: 1161895545.4832032\n",
      "Test set median imputation\n",
      "Linear Regression mse: 2212393764.74631\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Let's compare the performance of Linear Regression on zero vs median imputation of variables\n",
    "\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(X_train[cols_zero], y_train)\n",
    "print('Train set zero imputation')\n",
    "pred = linreg.predict(X_train[cols_zero])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_train, pred)))\n",
    "print('Test set zero imputation')\n",
    "pred = linreg.predict(X_test[cols_zero])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_test, pred)))\n",
    "print()\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(X_train[cols_median], y_train)\n",
    "print('Train set median imputation')\n",
    "pred = linreg.predict(X_train[cols_median])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_train, pred)))\n",
    "print('Test set median imputation')\n",
    "pred = linreg.predict(X_test[cols_median])\n",
    "print('Linear Regression mse: {}'.format(mean_squared_error(y_test, pred)))\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, models trained on the dataset in which missing data was replaced by the median show a worse predictive performance. If we look at the mean squared error of the test sets, we see that when we use median imputation, the difference between the estimated and real values is bigger (bigger mse).\n",
    "\n",
    "This suggests that the missingness of the data is a good indicator of house price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "House average price when data is missing:  160000.0\n",
      "House average price on complete cases:  164900.0\n"
     ]
    }
   ],
   "source": [
    "print('House average price when data is missing: ', data[(data.LotFrontage.isnull())|(data.GarageYrBlt.isnull())|(data.MasVnrArea.isnull())]['SalePrice'].median())\n",
    "print('House average price on complete cases: ', data.dropna()['SalePrice'].median())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "House average price when LotFrontage is missing:  172400.0\n",
      "House average price when GarageYrBlt is missing:  100000.0\n",
      "House average price when MasVnrArea is missing:  203287.0\n"
     ]
    }
   ],
   "source": [
    "print('House average price when LotFrontage is missing: ', data[(data.LotFrontage.isnull())]['SalePrice'].median())\n",
    "print('House average price when GarageYrBlt is missing: ', data[(data.GarageYrBlt.isnull())]['SalePrice'].median())\n",
    "print('House average price when MasVnrArea is missing: ', data[(data.MasVnrArea.isnull())]['SalePrice'].median())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that Houses for which the LotFrontage is missing are more expensive than the mean house. Similarly, houses where MasVnrArea is missing are also more expensive than the mean house. On the other hand, houses where GarageYrBlt is missing are substantially cheaper than the mean house. \n",
    "\n",
    "Therefore, missingness in these variables is a strong indicator of price, missingnesst would take the price up or down, and thus, replacing by zero (and therefore flagging the missingness), leads to better performance."
   ]
  },
  {
   "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": {
    "height": "556px",
    "left": "0px",
    "right": "1156.4px",
    "top": "107px",
    "width": "278px"
   },
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}