{ "cells": [ { "cell_type": "markdown", "metadata": { "button": false, "hide": true, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "# CS-109A Introduction to Data Science\n", "\n", "\n", "## Lab 7: $k$-NN Classification and Imputation\n", "\n", "**Harvard University**
\n", "**Fall 2019**
\n", "**Instructors:** Pavlos Protopapas, Kevin Rader, Chris Tanner
\n", "**Lab Instructors:** Chris Tanner and Eleni Kaxiras.
\n", "**Contributors:** Kevin Rader\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## RUN THIS CELL TO PROPERLY HIGHLIGHT THE EXERCISES\n", "import requests\n", "from IPython.core.display import HTML\n", "styles = requests.get(\"https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css\").text\n", "HTML(styles)" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "## Learning Goals \n", "In this lab, we'll explore classification models to predict the health status of survey respondents and be able to build a classification decision boundary to predict the resultsing unbalanced classes.\n", "\n", "By the end of this lab, you should:\n", "- Be familiar with the `sklearn` implementations of\n", " - $k$-NN Regression\n", " - ROC curves and classification metrics\n", "- Be able to optimize some loss function based on misclassification rates\n", "- Be able to impute for missing values\n", "- Be comfortable in the different approaches in handling missingness" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "button": false, "hide": true, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/_collections_abc.py:841: MatplotlibDeprecationWarning: \n", "The examples.directory rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2. In the future, examples will be found relative to the 'datapath' directory.\n", " self[key] = other[key]\n", "//anaconda3/lib/python3.7/_collections_abc.py:841: MatplotlibDeprecationWarning: \n", "The savefig.frameon rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n", " self[key] = other[key]\n", "//anaconda3/lib/python3.7/_collections_abc.py:841: MatplotlibDeprecationWarning: \n", "The text.latex.unicode rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2.\n", " self[key] = other[key]\n", "//anaconda3/lib/python3.7/_collections_abc.py:841: MatplotlibDeprecationWarning: \n", "The verbose.fileo rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n", " self[key] = other[key]\n", "//anaconda3/lib/python3.7/_collections_abc.py:841: MatplotlibDeprecationWarning: \n", "The verbose.level rcparam was deprecated in Matplotlib 3.1 and will be removed in 3.3.\n", " self[key] = other[key]\n", "//anaconda3/lib/python3.7/site-packages/seaborn/apionly.py:9: UserWarning: As seaborn no longer sets a default style on import, the seaborn.apionly module is deprecated. It will be removed in a future version.\n", " warnings.warn(msg, UserWarning)\n" ] } ], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import statsmodels.api as sm\n", "from statsmodels.api import OLS\n", "import sklearn as sk\n", "from sklearn.decomposition import PCA\n", "from sklearn.linear_model import LogisticRegression\n", "from sklearn.linear_model import LogisticRegressionCV\n", "from sklearn.utils import resample\n", "from sklearn.model_selection import cross_val_score\n", "from sklearn.metrics import accuracy_score\n", "# %matplotlib inline\n", "import seaborn.apionly as sns" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "## Part 1: General Social Survey Data + EDA\n", "\n", "The dataset contains a subset of data from the General Social Survey (GSS) that is a bi-annual survey of roughly 2000 Americans. We will be using a small subset of the approx 4000 questions they ask. Specifically we'll use: \n", "\n", "- **id:** respondant's unique ID\n", "- **health:** self-reported health level with 4 categories: poor, fair, good, excellent\n", "- **partyid:** political party affiliation with categories dem, rep, or other\n", "- **age:** age in years\n", "- **sex:** male or female\n", "- **sexornt:** sexual orientation with categories hetero, gay, or bisexual/other\n", "- **educ:** number of years of formal education (capped at 20 years)\n", "- **marital:** marital status with categories married, never married, and no longer married\n", "- **race:** with categories black, white, and other\n", "- **income:** in thousands of dollars\n", "\n", "Our goal is to predict whether or not someone is in **poor health** based on the other measures.\n", "\n", "For this task, we will exercise our normal data science pipeline -- from EDA to modeling and visualization. In particular, we will show the performance of 2 classifiers:\n", "\n", "- Logistic Regression\n", "- $k$-NN Regression\n", "\n", "So without further ado..." ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "### EDA" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "Do the following basic EDA (always good ideas):\n", "1. Determine the dimensions of the data set.\n", "2. Get a glimpse of the data set.\n", "3. Calculate basic summary/descriptive statistics of the variables.\n", "\n", "We also ask that you do the following:\n", "4. Create a binary called `poorhealth`. \n", "5. Explore the distribution of the responses, `health` and `poorhealth`, \n", "6. Explore what variables may be related to whether or not some is of poor health. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The dimensions of the data set are: 1569 observations and 10 variables.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
idhealthpartyidagesexsexornteducmaritalraceincome
01goodrep43.0malebisexual/other14.0never marriedwhiteNaN
12excellentdem74.0femalehetero10.0no longer marriedwhiteNaN
25excellentrep71.0malehetero18.0no longer marriedblackNaN
36gooddem67.0femalebisexual/other16.0no longer marriedwhiteNaN
47gooddem59.0femalebisexual/other13.0no longer marriedblack18.75
\n", "
" ], "text/plain": [ " id health partyid age sex sexornt educ \\\n", "0 1 good rep 43.0 male bisexual/other 14.0 \n", "1 2 excellent dem 74.0 female hetero 10.0 \n", "2 5 excellent rep 71.0 male hetero 18.0 \n", "3 6 good dem 67.0 female bisexual/other 16.0 \n", "4 7 good dem 59.0 female bisexual/other 13.0 \n", "\n", " marital race income \n", "0 never married white NaN \n", "1 no longer married white NaN \n", "2 no longer married black NaN \n", "3 no longer married white NaN \n", "4 no longer married black 18.75 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gssdata = pd.read_csv(\"data/gsshealth18.csv\")\n", "\n", "#####\n", "# You code here: EDA\n", "# 1. Determine the dimensions of the data set.\n", "# 2. Get a glimpse of the data set.\n", "####\n", "\n", "print(\"The dimensions of the data set are:\", gssdata.shape[0], \"observations and\", gssdata.shape[1], \"variables.\" )\n", "gssdata.head()\n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "good 771\n", "excellent 359\n", "fair 355\n", "poor 84\n", "Name: health, dtype: int64\n", "dem 708\n", "rep 514\n", "other 347\n", "Name: partyid, dtype: int64\n", "female 872\n", "male 697\n", "Name: sex, dtype: int64\n", "bisexual/other 907\n", "hetero 640\n", "gay 22\n", "Name: sexornt, dtype: int64\n", "married 655\n", "never married 458\n", "no longer married 454\n", "Name: marital, dtype: int64\n", "white 1137\n", "black 259\n", "other 173\n", "Name: race, dtype: int64\n" ] } ], "source": [ "# 3. Calculate basic summary/descriptive statistics of the variables.\n", "gssdata.describe()\n", "\n", "\n", "print(gssdata['health'].value_counts())\n", "print(gssdata['partyid'].value_counts())\n", "print(gssdata['sex'].value_counts())\n", "print(gssdata['sexornt'].value_counts())\n", "print(gssdata['marital'].value_counts())\n", "print(gssdata['race'].value_counts())" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "text/plain": [ "0.05353728489483748" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#####\n", "# You code here: EDA\n", "# 4. Create a binary called `poorhealth`. \n", "# 5. Explore the distribution of the responses, `health` and `poorhealth`, \n", "# 6. Explore what variables may be related to whether or not some is of poor health.\n", "####\n", "\n", "gssdata['poorhealth']=1*(gssdata['health']=='poor')\n", "gssdata['poorhealth'].mean()" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Question**: What classification accuracy could you achieve if you simply predicted `poorhealth` without a model? What classification accuracy would you get if you were to predict the multi-class `health` variable? Is accuracy the correct metric?" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Solution**: Poor health is a quite rare health status: only 5.35\\% of respondents said they were in poor health. If we predicted all persons to be in better than poo health, our naive classifier would have $1-0.0535 = 0.9465 = 94.65\\%$ accuracy. Acuracy is almost certainly not the ideal metric to use here: we'd be better off looking at false positive and false negative rate instead (it is more important to correctly classify those in poor health than those in better than poor health)." ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "### Data Cleaning - Basic Handling of Missingness" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "Let's begin by fitting an unregularized logistic regression model to predict poor health based on all the other predictors in the model and three $k$-NN models with $k=5,10,20$.\n", "\n", "First we need to do a small amount of data clean-up. \n", "1. Determine the amount of missingness in each variable. If there is *a lot*, we will drop the variable from the predictor set (not quite yet). If there is a little, we will impute.\n", "2. Drop any variables with lots of missingnes (in a new data set).\n", "3. Do simple imputations for variables with a little bit of missingness.\n", "4. Create dummies for categorical predictors.\n" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "text/plain": [ "id 0\n", "health 0\n", "partyid 0\n", "age 2\n", "sex 0\n", "sexornt 0\n", "educ 2\n", "marital 2\n", "race 0\n", "income 661\n", "poorhealth 0\n", "dtype: int64" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#########\n", "# 1. Determine the amount of missingness in each variable. \n", "# Use isna() in combination with .sum()\n", "########\n", "\n", "# Your code here\n", "\n", "gssdata.isna().sum()" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:14: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame.\n", "Try using .loc[row_indexer,col_indexer] = value instead\n", "\n", "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", " \n" ] }, { "data": { "text/plain": [ "age 2\n", "educ 2\n", "female 0\n", "marital_never married 0\n", "marital_no longer married 0\n", "race_other 0\n", "race_white 0\n", "sexornt_gay 0\n", "sexornt_hetero 0\n", "partyid_other 0\n", "partyid_rep 0\n", "dtype: int64" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#######\n", "# And then build your predictor set\n", "# 2. Drop any variables with lots of missingnes (in a new data set).\n", "# 3. Do simple imputations for variables with a little bit of missingness.\n", "# 4. Create dummies for categorical predictors.\n", "#########\n", "\n", "# get the predictors without a ton of missingness \n", "# (income was not included since it had so much missingness)\n", "X = gssdata[['partyid','age','sex','sexornt','educ','marital','race']]\n", "\n", "#create dummies (lots of ways to do it, two ways will be in the solutions\n", "# create dummies 2 different ways\n", "X['female'] = 1*(gssdata['sex']==\"female\")\n", "dummies = pd.get_dummies(X[['marital','race','sexornt','partyid']],drop_first=True)\n", "\n", "# add the dummies in via the join command.\n", "X = X.join(dummies)\n", "\n", "# let's drop the redundat variables no longer needed since we created the dummies\n", "X = X.drop(['partyid','sex','sexornt','marital','race'],axis=1)\n", "\n", "# now check the 'nulls'\n", "X.isna().sum()" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "# handle missingness in age and education\n", "# missingness in marital was handled with get_dummies\n", "\n", "# impute the median age\n", "X['age']=X['age'].fillna(X['age'].median())\n", "\n", "# impute the most common education: having a HS degree (13 years)\n", "# see histogram for justification\n", "X['educ']=X['educ'].fillna(13)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAQbElEQVR4nO3df6zddX3H8edrFHRTx8/CWNtYnc2i+0MkDevGZpw4xg9j2SILxowGmzRmmGjcMruZOLfsD9gy2VgWlk6MxTiFqYxGcNoAxiwZ6AX5KWoLqdK1o1WgaIib6Ht/nM91l9tze8+9vedc+vH5SE6+3+/n8zn3++73fO+r3/s5v1JVSJL68jPLXYAkaekZ7pLUIcNdkjpkuEtShwx3SerQiuUuAOC0006rtWvXLncZknRMueeee75TVSuH9b0gwn3t2rVMTU0tdxmSdExJ8q25+pyWkaQOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDr0g3qEq6XBrt966LPvdc9XFy7JfLS2v3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDI4V7kj1JHkxyX5Kp1nZKkp1JdrXlya09Sa5NsjvJA0nOHuc/QJJ0uIVcuf9WVZ1VVevb9lbg9qpaB9zetgEuBNa12xbguqUqVpI0mqOZltkIbG/r24FLZrTfUAN3ASclOfMo9iNJWqBRw72ALyS5J8mW1nZGVe0HaMvTW/sq4PEZ993b2p4nyZYkU0mmDh48uLjqJUlDrRhx3LlVtS/J6cDOJF8/wtgMaavDGqq2AdsA1q9ff1i/JGnxRrpyr6p9bXkAuBk4B3hierqlLQ+04XuBNTPuvhrYt1QFS5LmN2+4J3lJkpdNrwPnAw8BO4BNbdgm4Ja2vgO4vL1qZgNwaHr6RpI0GaNMy5wB3Jxkevy/VNW/J/kKcFOSzcC3gUvb+NuAi4DdwLPAFUtetSTpiOYN96p6DHjtkPbvAucNaS/gyiWpTpK0KL5DVZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0aOdyTHJfkq0k+27ZfkeTuJLuS3JjkhNb+ora9u/WvHU/pkqS5LOTK/d3AIzO2rwauqap1wFPA5ta+GXiqql4FXNPGSZImaKRwT7IauBj4cNsO8EbgU23IduCStr6xbdP6z2vjJUkTMuqV+98BfwL8uG2fCjxdVc+17b3Aqra+CngcoPUfauOfJ8mWJFNJpg4ePLjI8iVJw8wb7kneDByoqntmNg8ZWiP0/X9D1baqWl9V61euXDlSsZKk0awYYcy5wFuSXAS8GPh5BlfyJyVZ0a7OVwP72vi9wBpgb5IVwInAk0teuSRpTvNeuVfVn1bV6qpaC1wG3FFVbwfuBN7ahm0CbmnrO9o2rf+Oqjrsyl2SND5H8zr39wHvTbKbwZz69a39euDU1v5eYOvRlShJWqhRpmV+oqq+CHyxrT8GnDNkzA+AS5egNknSIvkOVUnqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUML+mwZSRqntVtvXZb97rnq4mXZ7zh55S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6tC84Z7kxUm+nOT+JA8n+YvW/ookdyfZleTGJCe09he17d2tf+14/wmSpNlGuXL/H+CNVfVa4CzggiQbgKuBa6pqHfAUsLmN3ww8VVWvAq5p4yRJEzRvuNfA99vm8e1WwBuBT7X27cAlbX1j26b1n5ckS1axJGleI825JzkuyX3AAWAn8CjwdFU914bsBVa19VXA4wCt/xBw6pCfuSXJVJKpgwcPHt2/QpL0PCOFe1X9qKrOAlYD5wCvHjasLYddpddhDVXbqmp9Va1fuXLlqPVKkkawoFfLVNXTwBeBDcBJSVa0rtXAvra+F1gD0PpPBJ5cimIlSaNZMd+AJCuBH1bV00l+FngTgydJ7wTeCnwS2ATc0u6yo23/Z+u/o6oOu3KX9MK0duuty12ClsC84Q6cCWxPchyDK/2bquqzSb4GfDLJXwFfBa5v468HPpZkN4Mr9svGULck6QjmDfeqegB43ZD2xxjMv89u/wFw6ZJUJ0laFN+hKkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA7NG+5J1iS5M8kjSR5O8u7WfkqSnUl2teXJrT1Jrk2yO8kDSc4e9z9CkvR8o1y5Pwf8UVW9GtgAXJnkNcBW4PaqWgfc3rYBLgTWtdsW4Lolr1qSdETzhntV7a+qe9v694BHgFXARmB7G7YduKStbwRuqIG7gJOSnLnklUuS5rSgOfcka4HXAXcDZ1TVfhj8BwCc3oatAh6fcbe9rU2SNCEjh3uSlwKfBt5TVc8caeiQthry87YkmUoydfDgwVHLkCSNYKRwT3I8g2D/eFV9pjU/MT3d0pYHWvteYM2Mu68G9s3+mVW1rarWV9X6lStXLrZ+SdIQo7xaJsD1wCNV9aEZXTuATW19E3DLjPbL26tmNgCHpqdvJEmTsWKEMecCfwA8mOS+1vZnwFXATUk2A98GLm19twEXAbuBZ4ErlrRiSdK85g33qvoPhs+jA5w3ZHwBVx5lXZKko+A7VCWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdGuUdqpLUtbVbb122fe+56uKx/Fyv3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QO+R2q0hEs53drSkdj3iv3JB9JciDJQzPaTkmyM8mutjy5tSfJtUl2J3kgydnjLF6SNNwo0zIfBS6Y1bYVuL2q1gG3t22AC4F17bYFuG5pypQkLcS84V5VXwKenNW8Edje1rcDl8xov6EG7gJOSnLmUhUrSRrNYp9QPaOq9gO05emtfRXw+Ixxe1vbYZJsSTKVZOrgwYOLLEOSNMxSv1omQ9pq2MCq2lZV66tq/cqVK5e4DEn66bbYcH9ierqlLQ+09r3AmhnjVgP7Fl+eJGkxFhvuO4BNbX0TcMuM9svbq2Y2AIemp28kSZMz7+vck3wCeANwWpK9wJ8DVwE3JdkMfBu4tA2/DbgI2A08C1wxhpolSfOYN9yr6m1zdJ03ZGwBVx5tUZKko+PHD0hShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUof8mj0dE/y6O2lhvHKXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIL+vQgvilGdKxwSt3SerQWK7ck1wA/D1wHPDhqrpqHPtZbst1FbvnqouXZb+Sjh1LHu5JjgP+EfhtYC/wlSQ7quprS70v+OmcJvhp/DdLWphxTMucA+yuqseq6n+BTwIbx7AfSdIcxjEtswp4fMb2XuBXZw9KsgXY0ja/n+Qbi9zfacB3FnnfcbKuhbGuhXuh1mZdC5Crj6qul8/VMY5wz5C2Oqyhahuw7ah3lkxV1fqj/TlLzboWxroW7oVam3UtzLjqGse0zF5gzYzt1cC+MexHkjSHcYT7V4B1SV6R5ATgMmDHGPYjSZrDkk/LVNVzSd4FfJ7BSyE/UlUPL/V+ZjjqqZ0xsa6Fsa6Fe6HWZl0LM5a6UnXYdLgk6RjnO1QlqUOGuyR16JgJ9yQXJPlGkt1Jtg7pf1GSG1v/3UnWTqCmNUnuTPJIkoeTvHvImDckOZTkvnb7wLjravvdk+TBts+pIf1Jcm07Xg8kOXsCNf3yjONwX5Jnkrxn1piJHa8kH0lyIMlDM9pOSbIzya62PHmO+25qY3Yl2TTmmv4mydfb43RzkpPmuO8RH/Mx1fbBJP814/G6aI77HvH3dwx13Tijpj1J7pvjvmM5ZnNlw0TPr6p6wd8YPDH7KPBK4ATgfuA1s8b8IfBPbf0y4MYJ1HUmcHZbfxnwzSF1vQH47DIcsz3AaUfovwj4HIP3JWwA7l6Gx/S/gZcv1/ECXg+cDTw0o+2vga1tfStw9ZD7nQI81pYnt/WTx1jT+cCKtn71sJpGeczHVNsHgT8e4bE+4u/vUtc1q/9vgQ9M8pjNlQ2TPL+OlSv3UT7SYCOwva1/CjgvybA3VC2ZqtpfVfe29e8BjzB4h+6xYCNwQw3cBZyU5MwJ7v884NGq+tYE9/k8VfUl4MlZzTPPo+3AJUPu+jvAzqp6sqqeAnYCF4yrpqr6QlU91zbvYvDekYmb43iNYqwfSXKkuloG/D7wiaXa34g1zZUNEzu/jpVwH/aRBrND9Cdj2i/CIeDUiVQHtGmg1wF3D+n+tST3J/lckl+ZUEkFfCHJPRl81MNsoxzTcbqMuX/hluN4TTujqvbD4BcUOH3ImOU8du9g8BfXMPM95uPyrjZl9JE5phmW83j9JvBEVe2ao3/sx2xWNkzs/DpWwn2UjzQY6WMPxiHJS4FPA++pqmdmdd/LYOrhtcA/AP82iZqAc6vqbOBC4Mokr5/Vv5zH6wTgLcC/DuleruO1EMty7JK8H3gO+PgcQ+Z7zMfhOuCXgLOA/QymQGZbtnMNeBtHvmof6zGbJxvmvNuQtgUfr2Ml3Ef5SIOfjEmyAjiRxf0JuSBJjmfw4H28qj4zu7+qnqmq77f124Djk5w27rqqal9bHgBuZvCn8UzL+TERFwL3VtUTszuW63jN8MT09FRbHhgyZuLHrj2p9mbg7dUmZmcb4TFfclX1RFX9qKp+DPzzHPtclnOt5cDvATfONWacx2yObJjY+XWshPsoH2mwA5h+VvmtwB1z/RIslTafdz3wSFV9aI4xvzA995/kHAbH/LtjruslSV42vc7gCbmHZg3bAVyegQ3Aoek/Fydgzqup5Thes8w8jzYBtwwZ83ng/CQnt2mI81vbWGTw5TfvA95SVc/OMWaUx3wctc18nuZ359jncn0kyZuAr1fV3mGd4zxmR8iGyZ1fS/0s8bhuDF7d8U0Gz7q/v7X9JYMTHuDFDP7M3w18GXjlBGr6DQZ/Lj0A3NduFwHvBN7ZxrwLeJjBKwTuAn59AnW9su3v/rbv6eM1s64w+FKVR4EHgfUTehx/jkFYnzijbVmOF4P/YPYDP2RwtbSZwfM0twO72vKUNnY9g28Vm77vO9q5thu4Ysw17WYwBzt9jk2/KuwXgduO9JhP4Hh9rJ0/DzAIrjNn19a2D/v9HWddrf2j0+fVjLETOWZHyIaJnV9+/IAkdehYmZaRJC2A4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI69H/HTRiwVWTnOAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(X['educ']);" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "25 50.0\n", "1478 65.0\n", "Name: age, dtype: float64\n", "529 24.0\n", "1546 75.0\n", "Name: age, dtype: float64\n" ] } ], "source": [ "# we checked these to see if there were any patterns in the missingness. \n", "# Nothing really showed up.\n", "print(gssdata['age'][pd.isna(gssdata['marital'])])\n", "print(gssdata['age'][pd.isna(gssdata['educ'])])" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "text/plain": [ "age 0\n", "educ 0\n", "female 0\n", "marital_never married 0\n", "marital_no longer married 0\n", "race_other 0\n", "race_white 0\n", "sexornt_gay 0\n", "sexornt_hetero 0\n", "partyid_other 0\n", "partyid_rep 0\n", "dtype: int64" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Just to make sure missingness is gone\n", "X.isna().sum()" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "## Part 2: Fit Basic Models\n", "\n", "In this section we ask you to:\n", "\n", "1. Split the data into 70-30 train-test splits (use the code provided...should have been done before EDA :( )\n", "2. Fit an unregularize logistic regression model to predict `poorhealth` from all predictors except income.\n", " \n", " 2b. If you have time: use 'LogisticRegressionCV' to find a well-tuned L2 regularized model.\n", " \n", " \n", "3. Fit $k$-NN classification models with $k=1,15,25$ to predict `poorhealth` from all predictors except income.\n", "4. Report classification accuracy on both train and test set for all models." ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [], "source": [ "#######\n", "# Use the following train_test_split code to: \n", "# 1. Split the data into 70-30 train-test splits\n", "#######\n", "from sklearn.model_selection import train_test_split\n", "itrain, itest = train_test_split(range(gssdata.shape[0]), train_size=0.70)\n", "\n", "# Note: the train-test split above is for the INDICES for splitting in case we \n", "# want to use them again in the future...we can have an identical split\n", "X_train = X.loc[itrain]\n", "X_test = X.loc[itest]\n", "\n", "y_train = gssdata['poorhealth'][itrain]\n", "y_test = gssdata['poorhealth'][itest]" ] }, { "cell_type": "code", "execution_count": 108, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAELCAYAAAAhuwopAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXgUVdr///cJIeyBBJCAhEVxZBVEnQFlTBCRURSVZUBAQQFHIejvERUQlSgqIhEHUfD3yBKUAUQUxA0fhiGByIgioAhRAVkVVBYBSUgIOd8/OokJWehKOqnu5PO6rr6Krjp16u6m03efOnVOGWstIiIiTgS5HYCIiAQeJQ8REXFMyUNERBxT8hAREceUPERExLFgtwMoC8aYzUBz4Hdgp8vhiIgEihZATWC3tfby3BtMRbhU1xjzG1Db7ThERALUcWttndwrKkTLA0+Lo3bt2rXp0KGD27GIiASELVu2cPz4cfB8h+bhdfIwxlwK/A24CrgS+BNggH7W2qXFDc4YMxC4H7gMqAR8C8wDZllrM4tb7zl2Ahd26NCBhIQEH1UpIlK+RUdHk5iYCAWc7nfS8rgfeNBXQQEYY14FRgKngdXAGaAb8ArQzRjTz1p71pfHFBGRknNytdU3wFSgP55OlMSSHNgY0wdP4jgEXGatvdlaeztwCZAM3A7ElOQYIiJSOrxueVhrZ+d+bowp6bHHZy3HWmt35DrOz8aY+4EEYJwxZoYPT1+JiIgPuDLOwxjTGLgCSAfePne7tTYR+BGIADqVbXQiInI+bg0SzL5eeJu1NrWQMl+cU1bEa8aYfA8R8R23kkfzrOXeIsrsO6esiFdyJ4qoqKgC14tIybg1zqNm1vJUEWWyryuuVdBGY8xQYKiXx9Pgjgoo9wBYJQ4R33IreWT/JZdkeHszIOp8haRiyt3iyH6edb16HpmZmRw9epSTJ0+SlpZGRZhxQSoWYwzBwcFUr16d0NBQatasef6dvOBW8jiZtSzqVWRvO1nI9j14f7lwBzQ9SYVybqIoLHHs37+flJSUsgpLpMxZazlz5gzHjx/n+PHjhIeHc8EFF5S4Ne5W8tiTtWxaRJnIc8rmYa2NB+K9OZgxJgG1UiocY0yhLQ6Ao0ePkpKSQnBwMBEREdSoUYOgIE00LeVLZmYm6enpnDx5kiNHjnD06FGqVq1K7dol+z3t1l/K5qxlG2NMtULKXHVOWRGv5D71lDtxnHtK6uRJT6M2IiKCWrVqKXFIuRQUFETVqlWpX78+DRo0AODYsWMlr7fENRSDtXY/sAkIAfqdu90YEwU0xjP6/L9lG52UB9bafI9zpaWlAVCjRo2yDk/EFaGhoQCcPn26xHWVavIwxkw2xnxrjJlcwObsdVOMMS1y7XMBMDPr6fMaXS6lJTuhqMUhFUX2Z90XF4Y4mVW3I398qQO0zlo+Z4x5OHultTb3iPCGwKVZyzystUuNMbPwTLi41Rjzb/6YGDEUWI5ngkQREfEBX16y7qTDPBT4SwHrLynuwa21I40xScAoPB3a2VOyz8W3U7KLiIgPOZkYMYE/xmd4u89QzjOQz1q7EFjopF4REXGXTvaKiIhjSh4iUqjvvvuO6dOnM3jwYFq2bElQUBDGGJYuLfbNQ70WHx+PMYahQ4eW+rGKMnToUIwxxMfHl+pxEhISMMYQHR1dqsfxlYpyD3MRKYZZs2Yxffp0t8MoF5o1a8bevXvZvXs3zZo1czucElPLQ0QK1bZtWx555BHeeustdu7cmW/OsIpg8uTJJCcnc/vtt5fqcf785z+TnJzMG2+8UarH8RW1PESkUMOHD3c7BNc1bNiQhg3zjTbwuerVq9OyZctSP46vqOUhIgFp/fr19OnTh4iICEJCQoiIiKBv37589tlnhe7zyy+/MHLkSBo3bkzVqlVp0aIFEyZMIDU1lejoaIwxJCQk5NmnsD6Ps2fP8tprr3H11VdTu3ZtQkJCaNCgAR07dmTMmDH8+uuvwB99N3v3em5f1Lx58zw3KduzZw9w/j6PI0eO8OSTT3L55ZcTGhpKjRo1uOSSSxg6dCjr168v1ntYEmp5iEjAmTVrFjExMWRmZnLVVVdx3XXXsXPnTt555x2WLVvGa6+9xogRI/Ls89NPP3HNNdewZ88eLrjgAm655RbS0tJ4+eWX8yUMbwwbNoz58+dTrVo1unTpQr169Th8+DC7du1i2rRp9OvXj/r169OiRQuGDBnC0qVLOXXqFH369MkzLbo3U6Rv3ryZnj17cvDgQcLDw4mOjqZq1ars3buXRYsWAXD11Vc7fg0loeQhIj41dOhQ5s+f73g/bzuSv/rqKx544AEAlixZQr9+f0yPt3jxYgYNGsSoUaPo3Lkzbdu2zdk2cuRI9uzZw4033sjbb7+dM6fZoUOH6NatG9u3b/c61r179zJ//nwiIyP54osvciYczLZlyxYaNWoEQJcuXejSpQsJCQmcOnWKuLg4Rx3mJ0+epFevXhw8eJD77ruPadOmUa3aH/PJ/vrrr3z33Xde1+crSh4i4lNdunQp1n7e3qTo5ZdfJiMjg4EDB+ZJHAADBgxg2bJlLFmyhOnTp/P6668Dni/7FStWEBwczMyZM/NMhhkREUFcXBw33XST17H+8ssvAHTs2DFf4gDo0MF3Ny+dM2cOBw4coFOnTsycOTPfFCP169enfv36Pjuet5Q8RIoQiLevdftuiMOHDy/VjvbsafYLG/9xzz33sGTJkjynotauXYu1ls6dOxf4q//GG28kLCzM66nKW7ZsSa1atfjwww957rnnGDRoEE2bFnV7ouJbuXIl4DlN5k+fR3WYixShoKnd/f1R3v3444+Ap+O5IBdffHGecrn/XdQXfJMmTbyOoVatWsydO5dq1aoxYcIEmjVrRuPGjenXrx/x8fE+mfI8W3ZHu79diaWWh4j41OzZs0lKSnK8X1xcHPXq1fO6fGG/wotKoEX9cnc6NX/fvn25/vrree+991i7di2ffvopS5cuZenSpcTGxrJu3ToiIyPPX1GAUvIQEZ9KSkoqVod5bGysV8njwgsvZNeuXfzwww85rYzcdu/enVMuW3bndfav+IIUta0wderUYciQIQwZMgSAXbt2MWLECNasWcPYsWNZuLDkc742bdqUb7/9lu+++67Y/UmlQaetRMSn4uPji3W6zdsrkLJHuRc2EnvevHkAecZL/PWvf8UYw/r16wtMEp988glHjx519kILcPHFFzNhwgTAc1VYbiEhIQBkZGQ4qrNHjx4AzJ07169OSyp5iEhAeeCBBwgODmbRokUsW7Ysz7a3336bJUuWULly5ZzLecHTP9KzZ08yMjIYNWoUKSkpOdt+/vlnHn74YZzYvHkzb731Fqmpqfm2vf/++0D+/pXsllBycrKjYw0fPpxGjRqxfv16Ro8ena8/5ddffy3WacKS0mkrESnUpk2bGDlyZM7z7LEQjz32GHFxcTnrixrV7Wvt27dn+vTpxMTE0Lt3b/7yl79w8cUXs3PnTj7//HOCgoJ45ZVXaNeuXZ79Zs2axddff82HH37IRRddxLXXXktaWhpr1qyhTZs2dOrUic8++yynhVCUvXv3MmDAAKpXr07Hjh2JjIwkPT2dzZs388MPP1CrVi2efvrpPPvcfvvtJCQkMGjQIG644Qbq1KkDwJQpU6hbt26hx6pVqxbvvfcePXv25NVXX2Xx4sVcc801OYMEN2/ezB133FHmp7SUPESkUCdOnGDDhg351u/YscOFaP4wcuRI2rdvz4svvsinn37Kl19+SXh4OL179+bhhx+mc+fO+fZp3Lgxn3/+ORMnTmTFihW89957XHjhhdx///1MnDiRyy67DMCrfpdOnToxefJkEhMT+fbbb/nyyy8JCQkhMjKSMWPGMHr06Hwtj5iYGE6cOMG//vUvPvjgA9LS0gB4/PHHi0weAFdeeSVbt27lpZde4v3332fVqlUEBQXRqFEjBg4cyD/+8Q9v3zqfMf50Dq20GGMSgKioqKhiTUMg5VP26YNWrVq5HIm4bc+ePbRo0YIaNWpw7Ngxx1deBRInn/vo6OjscTWJ1tro3NvK7zskIpKLtZaNGzfmW79//37uvPNOzp49y1133VWuE4cv6bSViFQIZ8+e5aqrrqJJkya0bNmSsLAw9u/fz6ZNmzh9+jRt27blmWeecTvMgKHkISIVQqVKlZgwYQL//ve/2bx5M7/99htVqlShdevW9O7dmwcffNDr+bVEyUNEKghjDM8884xaFz6ik3siIuKYkoeIiDim5CEiIo4peYiIiGNKHiIi4piSh4iIOKbkISIijil5iIiIY0oeIiLimJKHiIg4puQhIoU6c+YMq1evZsyYMXTq1ImGDRsSEhLChRdeSN++fUv9FgexsbEYY4iNjS3V45TE4cOHmTt3Lvfffz9XXXUVVapUwRhDTEyM26GVKs1tJSKFSkxMpHv37gBERERwxRVXUKNGDbZv384777zDO++8wxNPPJHvrnkVSVJSEsOGDXM7jDKnloeIFCooKIg+ffqwdu1aDh48yAcffMBbb73F1q1bWbx4MZUqVWLSpEmsWbPG7VBd06BBA+6//35mz57N5s2bmTBhgtshlQm1PESkUNdddx3XXXddgdv69+/PqlWrmDNnDgsWLKBr165lHJ1/6Ny5c57b3i5fvtzFaMqOWh4iUmyXX345AAcOHCjzYycmJhIeHk5ISAjz58/PWR8dHY0xhoSEBL788kt69epF3bp1qVatGu3bt2fOnDkF1lfc/SoqJQ8RKbYdO3YA0LBhwzI97uLFi+nRowcZGRl8+OGHDBkyJF+ZlStX0rlzZ3bv3s0NN9xAx44d+frrrxk+fDgvvvhioXUXd7+KRslDRIrl0KFDxMfHA9CnT58824YOHYoxxvFjz5495z3u1KlTGThwIHXr1mXdunU5HfrnmjJlCq+99hpbt25l0aJFfPrpp7z55psAPP3006SkpPh0v4pGfR4i4lhGRgaDBw/m+PHjdOvWjVtuuSXP9i5duhSr3qJuA3v27FkeeOABZs6cSZs2bfj444+JjIwstHyfPn2455578qwbPHgwzz33HMnJyWzcuJFrr73WZ/tVNEoeIkUwxrgdgmPW2lI/xn333cfq1auJjIxkwYIF+bYPHz6c4cOH++x4KSkp9O7dmxUrVtC1a1eWLVtG7dq1i9zn5ptvLnB9y5YtSU5O5qeffvLpfhWNkodIEcriizjQPPjgg8yZM4eIiAhWr15NREREqR/zpZdeIiMjgw4dOrBy5UpCQkLOu0+TJk0KXB8aGgrA6dOnfbpfRaPkISJeGzNmDC+//DL169dn9erVXHLJJQWWmz17NklJSY7rj4uLo169evnW9+zZk6SkJLZs2cK0adMYN27ceesKCipel25x96tolDxExCuPPvoo06ZNo27duqxatYrWrVsXWjYpKSnP5bPeio2NLTB5dOjQgUmTJtG9e3fGjx9PSkpKhR7V7g+UYkXkvMaNG8fUqVMJCwtj1apVtG/fvsjy8fHxWGsdP5o1a1Zone3atWPt2rU0btyYSZMm8cgjj/j4VYoTSh4iUqQnnniCKVOmUKdOHVatWpUzMNANf/rTn1i7di3NmzcnLi6OmJgY9Uu5RKetRKRQK1as4JlnngGgRYsWzJgxo8ByLVu29KofwheaN2/OunXr6NatG6+++iqpqam8/vrrrvZVdOrUKeff2aPtly5dysaNG3PWz5w5k44dO5Z5bKVFyUNECnX06NGcf2/cuDHPl2FuUVFRZZY8AC688MKcGX/nzp1Lamoqb7zxBsHB7nylbdiwId+6n3/+mZ9//jnn+YkTJ8oypFJnKkKTzxiTAERFRUWV+v0HJHAkJycD0KpVK5cjESk7Tj730dHRJCYmAiRaa6Nzb1Ofh4iIOKbkISIijil5iIiIY0oeIiLimJKHiIg4puQhIiKOKXmIiIhjSh4iIuKYkoeIiDim5CEiIo4peYiIiGNKHiIi4piSh4iIOKbkISIijil5iEiRZsyYwd///ndatWpF3bp1qVy5MvXr1+f6669nwYIFpXonv6FDh2KMIT4+vtSOUVL79+9n1qxZDBs2jMsuu4zg4GCMMcTFxRW5X2xsLMaYQh9Vq1Yto1dQPLoZlIgUacqUKfzyyy+0bduWq6++mho1arB3717+85//sHr1apYuXcq7777r6p383PTOO+/wP//zP8Xev3379nTo0CHf+sqVK5ckrFKn5CEiRVq8eDGXX345NWrUyLN+27ZtdOvWjffee4/58+dz9913uxShu5o3b86DDz7IFVdcwZVXXsnkyZN58803vd7/tttuIzY2tvQCLCUV86eCiHitS5cu+RIHQJs2bRg1ahQAq1atKuuw/Matt97KP//5T+68805atWpVYVpgFeNVikipyL5nuBvn599++22qVq1KrVq1+OSTT3LWN2vWDGMMe/bsYdWqVXTr1o3atWtTvXp1OnXqxIoVKwqsr7j7VVRKHiJSLLt37+a1114D4JZbbinTY0+bNo3+/fsTHh7O2rVr6dGjR74yc+bMoUePHvz+++/cdNNNtGzZkg0bNnDbbbexdOnSQusu7n7FtWnTJsaOHcu9997LuHHjWLZsGenp6T4/js9Za8v9A0gAbFRUlBXJtn37drt9+3a3wwgYc+fOtUOGDLEDBw601157rQ0ODrZBQUF2/Pjx+cpGRUVZwPHjXEOGDLGAnTdvnrXW2rNnz9rRo0dbwLZu3dru2bMn3z5Nmza1gA0JCbEff/xxnm2TJk2ygG3RooXP9iss5qlTpxZZbuLEiYW+D40bN7YJCQnnPVZxOPnc5/p/TLDnfK+qw1xEvPLpp58yf/78nOfBwcFMmjSJhx56KF/Zv/3tbzRr1synx09NTWXgwIEsX76ca6+9luXLlxMWFlZo+dGjR/O3v/0tz7pHH32UuLg4du7cyb59+2jSpInP9nPq4osvZvLkydx44400b96c9PR0tm7dylNPPUViYiI33XQT69evp3379iU+VmlQ8hApgjHG7RAcs6U07mL27NnMnj2b1NRUdu/ezbx584iNjWXJkiV89NFHNGrUKKfsuHHjfHrsw4cPc9111/HZZ58xYMAA4uPjqVKlSpH73HzzzfnWhYSEcNFFF7F582Z++umnApNAcfdz6s4778y3rmvXrnTt2pW+ffvyzjvvMGHCBD744IMSH6s0qM9DpAjnNtUD4VHaqlWrRuvWrZk6dSqTJ0/mq6++IiYmplSPOX78eD777DN69uzJwoULz5s4gEK/4ENDQwE4ffq0T/fzpSeffBLwXMV25syZUj9ecajlISLFdvfdd/Pwww/z/vvvc+bMmZyBbc8//zzffvut4/oKG0ner18/li1bxsqVK1m4cCGDBg06b13FvWTWHy61bdmyJQDp6ekcPnyYhg0buhxRfkoeIlJsderUITg4mIyMDI4ePUqDBg0AWLlyJYmJiY7rKyx53HDDDYwYMYJevXpx1113kZqayvDhw0sSul87cuRIzr9r1qzpYiSFcz/FikjAWrt2LRkZGdSpU4d69erlrE9ISPD5KbeuXbvyySefUKtWLe69915mzJhR2i/PNUuWLAHg0ksvpVatWi5HUzAlDxEp1Lp16/jXv/5FWlpavm2ffvopw4YNA2DYsGFUqlSp1OO5+uqrWb16NeHh4TzwwAO88MILpX7M0rBv3z4WLlyY73211vLmm28yfvx4gBLNmVXadNpKRAq1a9cu7r77bmJiYujYsSMRERGcPHmSXbt2sX37dgB69uzJpEmTyiymK664goSEBK6//nrGjh1LSkqKq3NDHTx4kNtvvz3n+a5duwDPbMS5BxUuW7Ysp+/i6NGjDBo0iPvuu49LL72UJk2akJ6ezrZt29i9ezcAMTEx/OMf/yjDV+KMkoeIFCoqKoonnniCdevW8f3337N+/XqstURERNCnTx8GDx7MbbfdVuZxtW3blrVr19KtWzeeeuopUlNTmTJlSpnHAZCWlsaGDRvyrd+3bx/79u3LUy5bZGQkjzzyCF988QU7d+5k27ZtZGZmEhERQf/+/bn33nu57rrryiT+4jJlcWmf24wxCUBUVFQUCQkJLkcj/iI5ORmAVq1auRyJSNlx8rmPjo7OvvAh0VobnXub+jxERMQxJQ8REXHMcfIwxgw0xqwzxhw3xvxujNlojBlljHFUlzEm1hhji3iU/jBOEREpFkcd5saYV4GRwGlgNXAG6Aa8AnQzxvSz1p51GMNXwJYC1vvnmHwREfE+eRhj+uBJHIeAa621O7LWNwDWALcDMcB0hzEst9bGOtxHRERc5ORU0/is5djsxAFgrf0ZuD/r6Tinp69ERCTwePVFb4xpDFwBpANvn7vdWpsI/AhEAJ18GaCIiPgfb09bXZ613GatTS2kzBfAhVll1zuIoaMxZgoQBhwFNgAfWmsD4D6MIiIVk7fJo3nWcm8RZbKHUjYvokxBbsl65HbAGDM4q0UjIiJ+xtvkkT0n8KkiyvyetfR2CshdePpRPgZ2AyFAO2AiEAV8ZIy52lr7VUE7G2OGAkO9PFYHL8uJiIgXvE0e2ffi9NlcJtbaNwtYvQZYY4xZCvQBngXy3xPSoxmeJCMiImXM2+RxMmtZ1F1JsredLKKMt57Gkzy6G2MqW2sLGvOxB/D2tFYHoLYP4hIREbxPHnuylk2LKBN5TtmSyL5/ZQhQDzh4bgFrbTwQ701l2RMj+iAuERHB+3Eem7OWbYwx1Qopc9U5ZUuibq5//15oKRERcYVXycNaux/YhKcl0O/c7caYKKAxntHn//VBXH/PWn5nrfXFaTAR8ZHHHnsMYwzGGOLi4krtOLGxsRhjXL3R0/ls3ryZ5557jm7dutGsWTOqVKlCeHg4Xbt2Zd68eWRmZha4X3x8fM57WNjj0KFDZfxqnHEyt9VkPAMEpxhj1ltrdwIYYy4AZmaVed5am/NuGWNi8ExZ8rm19q5c65sAXYB3rLVpudYbYHDWsQBecv6SRKS0fPHFF7zwwgsYY857z/HyLiMjg44dOwJQs2ZNrrrqKjp37syBAwdYt24dCQkJLF68mPfee4+qVasWWMfFF19Mly5dCtxWrVphJ3n8g9fJw1q71BgzC89UJFuNMf/mj4kRQ4HleCZIzK0ecCmeFklu4cC/gNeMMd/hGSMSArThj3Eir1hr/39nL0dESktaWhpDhw6lQYMG/PnPf2b58uVuh+S6K664grFjx9KrVy+qVKmSs37r1q306NGD//u//2Py5Mk89dRTBe7fpUsX4uPjyyha33I0D5W1diQwCM8prCigB7ATT+uij4MZdfcDU4Ev8UxpciPQPSuet4Bu1trRTmITkdL15JNPsn37dl577TVq19bFi8HBwWzcuJF+/frlSRwA7dq144UXXgBgwYIFboRX6hxPYmitXWitvcZaG2qtrWGtvcJa+2ru01W5ysZaa825ty+01h6x1j5qre1qrY201la31la11jaz1g6w1v6nBK9JRHxsw4YNvPjiiwwcOJBbbjl3Qoiyl5iYSHh4OCEhIcyfPz9nfXR0NMYYEhIS+PLLL+nVqxd169alWrVqtG/fnjlz5hRYX3H3K8rll3tmdTpw4EDxXqSfc3Q/DxGpeE6fPs2QIUMIDw9n+nSnd1zwvcWLFzN06FBCQkL48MMP6d69e74yK1euZNq0aVx66aXccMMN7Nu3j/Xr1zN8+HB+++03xowZU2Ddxd2vIDt2eCYfb9iwYaFldu7cyeOPP84vv/xCaGgoHTt2pFevXtSsWdSQOj9hrS33DyABsFFRUVYk2/bt2+327dvdDsPvPfTQQxawixcvzlk3ZMgQC9ipU6fmK5+9zelj9+7deeqZOHGiBezEiRNz1r3wwgvWGGMbNWpkt2zZku/YUVFROfXNmTMnz7Y333zTAjY0NNSeOnXKJ/sVJjMz03bu3NkC9qGHHsq3fd68eYW+D2FhYfbtt9/26jjF4eRzn+t9SbDnfK+q5SEihVq/fj3//Oc/ue222+jfv79X+xR29dD5FPVr++zZszzwwAPMnDmTNm3a8PHHHxMZGVlo+T59+nDPPffkWTd48GCee+45kpOT2bhxI9dee63P9jvXU089xX//+18aNGjA+PHj821v2LAhjz/+OL169eKiiy4iODiY5ORkXnjhBZYtW0b//v356KOP6NGjx3mP5RYlD5EieK4eDyzWR5fQpqamcvfddxMaGsrMmTPPv0OW4cOHM3z4cJ/EAJCSkkLv3r1ZsWIFXbt2ZdmyZeftsL/55oKnxGvZsiXJycn89NNPPt0vtzfeeIOnn36akJAQFi1aRL169fKV6dGjR77E0KlTJ959913GjBnDtGnTGDNmjJKHSKDy1RdxIHrsscf4/vvvmTt3bpHn7UvbSy+9REZGBh06dGDlypWEhIScd58mTZoUuD40NBTw9OP4cr9sb7/9Nvfccw+VKlVi8eLFdO3a9byxnuvxxx9n+vTpbNu2jX379hUak9uUPESkQMuWLSMoKIj58+fnuaIJ4NtvPdPPzZo1iw8++IAWLVowe/ZsAGbPnk1SUpLj48XFxRX4K71nz54kJSWxZcsWpk2bxrhx485bV1BQ8e6GXdz9AN59910GDhyItZYFCxZw++23F6uesLAwLrjgAg4ePMiPP/6o5CEigSczM5PExMInr/7hhx/44Ycf+O2333LWJSUl5Us23oiNjS0weXTo0IFJkybRvXt3xo8fT0pKCk8//bTj+kvT8uXLGTBgAJmZmcybN48BAwYUu66zZ89y/PhxoOh+ILcVP82KSLm2Z88erLUcPnyYb775hi+++IJvvvmGw4cPM2TIEACmTp2KtZYtW7bk7BcfH1+sqyKbNWtWaCzt2rVj7dq1NG7cmEmTJvHII4+U9sv32vvvv8/f//53MjIymD17Nnfdddf5dyrCBx98QEpKCrVq1aJly5Y+itL3lDxEpFBHjhzhp59+IjIyko4dOxIZGclPP/1EWlra+Xf2sT/96U+sXbuW5s2bExcXR0xMjOt9Uh999BF9+/YlIyOD//3f/+Xuu+8+7z4pKSnMmjWL33/PPyGy5G4AABAKSURBVGH4Rx99xIgRIwAYNWoUlStX9nnMvqLTViJSqEOHDtG0adOcDuPQ0FCaNm1KamqqK/E0b96cdevW0a1bN1599VVSU1N5/fXXS9RXUVy//PILvXv3Jj09ncaNG5OUlFRoX0/u+avS09MZOXIkDz30EK1atSIyMpJKlSqRnJyc05fUu3dvvzs1dy4lDxEpVGpqar7z7jVr1iQjI8OliODCCy8kMTGR7t27M3fuXFJTU3njjTcIDi7br7OUlJScFtiBAweK7OfJnTyqV6/O448/zueff853333H999/T3p6OvXr16dXr14MGTKE3r17l3b4JWbcbvaVhew7CUZFRZGQkOByNOIvkpOTAWjVqpXLkfivbdu2ERkZmdPyADhx4gT79++nTZs2LkYmxeXkcx8dHZ19wUSiPWeOQvV5iEihIiIi2Lt3LydOnCAzM5MTJ06wd+9eIiIi3A5NXKbk4ScWLVpE27ZtqVSpEm3btmXRokVuhyRC3bp1adSoEfv372fTpk3s37+fRo0aUbdu3fPvLOWa+jz8wKJFi5gwYQJz5syhS5cuJCUlMWzYMADuuOMOl6OTiq5u3bpKFpKPWh5+4Nlnn2XOnDl07dqVypUr07VrV+bMmcOzzz7rdmgiIgVS8vADycnJ+WYi7dKlS07HloiIv1Hy8AOtWrXKd314UlKSrgISEb+l5OEHJkyYwLBhw1izZg1nzpxhzZo1DBs2jAkTJrgdmoiUI74cmqEOcz+Q3Sk+evRokpOTadWqFc8++6w6y0uZMQZrLZmZma6MUBYpa5mZmYBv7lOj5OEn7rjjDiWLMlalShVOnz7NqVOnqFWrltvhiJS6EydOAFC1atUS16WfW35C4zzKXnbCOHToECdPniQzM9P1ifb80bZt29i4cWPOY9u2bW6HJF7KblmfPn2aX3/9lUOHDgGee4aUlFoefkDjPNwRHh7OqVOnSElJ4cCBA26H45dSU1PJzMwkODiYkJAQ0tPTSU1NZdOmTVSrVs3t8KQYwsPD80w3U1xqefgBjfNwR1BQEJGRkdSvX5+qVasG5P3KS9u2bds4ePAgVapUwRhDlSpVOHjwoFofAcQYQ+XKlalduzaRkZE0aNBAfR7lhcZ5uCcoKIh69eoVeAc7gdatW/Prr7/meX/q169P/fr1dYqvglPLww9onIf4s+xTqIU9l4pJycMPaJyH+Kt27dqxYsUKbr31Vg4fPsytt97KihUraNeunduhict02soPaJyH+Kuvv/6ayy67jBUrVlC/fn3Ak1C+/vprlyMTtyl5+AmN8xB/pUQhBdFpKxERcUzJQ0REHFPyEBEpQ+VlNgklDz9RXj5QIlK47NkkZsyYwenTp5kxYwYTJkwIyL93JQ8/UJ4+UCJSuPI0m4SShx949tlnGThwIKNHj6Zq1aqMHj2agQMHBuQHSsofY0y+hxRPeZpNQsnDD2zfvp2FCxfmaXksXLiQ7du3ux2aVHC5E8XixYsLXC/eK0+zSSh5+IGQkBBiYmLyNGVjYmIICQlxOzQRwDO1d//+/TWfVQmVp9kkNEjQD6SnpxMbG8u4ceM4c+YMlStXpmrVqqSnp7sdmkieFkf28wEDBrgUTWArT7NJmIrwS8IYkwBERUVFkZCQ4HI0+dWtW5djx45xwQUX8Msvv+Qsw8LCOHLkiNvhSQWWfXoq9/dEQeukfIqOjiYxMREg0VobnXubTlv5gRMnThAWFsaiRYtIS0tj0aJFhIWF5dwyUsRtxhjeeust9XX4QHm5AEHJww9kZGQQFxeX52qruLg4MjIy3A5NKrjcrYvcp6rU6iie3Ini5ZdfLnB9oFDy8ANVqlTh2LFjfPPNN5w9e5ZvvvmGY8eOUaVKFbdDE8Fam+8hJWOtZfTo0QH9Xip5+IERI0YwduxYpk2bRkpKCtOmTWPs2LGMGDHC7dBExMdytzgKeh4o1GHuJ3r06MGqVauw1mKMoXv37nzyySduhyUiPhRoFyCow9zPLVq0iB07drB69WrS09NZvXo1O3bs0PQkIuWUMYYZM2YEZF9HNiUPP6DpSUQqhtytiwceeKDA9YFCgwT9wPbt29m5cydpaWkAbNu2jZ07d2qQoEg5FIiJoiBqefiJtLQ0wsLCMMYQFhaWk0hE3FZexiWIbyl5+IHsXyInTpzAWpszOLC8/EKRwJU7Udxyyy0FrpeKSaet/MjZs2fzLEX8RUFXB0nFppaHHwkLC8uzFPEHuVscBT2XvAo6zVeSh79Sy8OPHDt2LM9SxB+8//77RT6XvLw93WyMCehT02p5iMh5GWPo1auXX/8SlrKl5OFHsv8w9Qcq/iL3L+PcLY5A/sUsvqHTVn4k+w9Sf5jiT/R5lIKo5SEiIo4peYiIiGNKHiIi4piSh4iIOKYO8zJS3CuoCttPnZgi4iYljzJS1Jd9UFBQgduNMWRmZpZmWCIixaLTVn5g1KhRjtaLiLhNLQ8/MGPGDABef/110tLSqFKlCiNGjMhZLyLib5Q8/MSMGTNybkt5+vRpt8ORCsbXsxqoT678U/IQEa++7AN9Ij/xLfV5iIiIY0oeIiLimJKHiIg4puQhIiKOKXmIiIhjSh4iIuKYkoeIiDim5CEiIo4peYiIiGMaYS4i4qXw8HCOHTvms/p8NS1MWFgYR48e9Uld3lLy8AF9oEQqhmPHjvnlFC2+npvMG0oePqAPlPit2No+q8pODPVpfcQe911dUuaUPETKMfPUCb/9YWNj3Y5CSkId5iIi4piSh4iIOKbkISIijqnPwwd83pHoI3ZiqNshiEg55Th5GGMGAvcDlwGVgG+BecAsa21mMer7G/AQcCVQFfgBWATEWWvTnNbnBnVKilQM+qH4B0fJwxjzKjASOA2sBs4A3YBXgG7GmH7W2rMO6nsUmAKcBRKAY0AU8AxwszGmm7U2xUmMIiKlRT8U/+B1n4cxpg+exHEIuMxae7O19nbgEiAZuB2IcVDflcDzQApwjbX2emttP+AiYC3QCXjW2/pERKTsOOkwH5+1HGut3ZG90lr7M57TWADjjDHe1jkOMMAUa+2GXPX9DtwNZAIjjTF1HMToGmOM3z3CwsLcfltEpJzy6oveGNMYuAJIB94+d7u1NhH4EYjA02I4X30hwI1ZT/9VQH0/AP8FQoCbvInRTdZanz18WZ+mJhGR0uJtK+HyrOU2a21qIWW+OKdsUS4FqgNHrbW7fFCfiIiUIW87zJtnLfcWUWbfOWW9qW9fEWWKrM8YMxQY6sWxADp4WU5ERLzgbfKombU8VUSZ37OWtcqovmZ4rswKCE4mKfSmrD9e8eEzfngpJBCwE/n54wSZgdwfp/fTw9vkkf1u+eobyxf17QESvSzbAXD1G6lcf9n7WoB+Sfsjfe58S+/nH7xNHiezljWLKJO97WQRZXxWn7U2Hoj34lgYYxIIoFaKiIi/87bDfE/WsmkRZSLPKetNfU18VJ+IiJQhb5PH5qxlG2NMtULKXHVO2aJ8C6QC4caYiwsp82cH9YmISBnyKnlYa/cDm/CMu+h37nZjTBTQGM/o8/96UV868HHW00EF1HcR0BnPuJIPvYlRRETKjpMR5pOzllOMMS2yVxpjLgBmZj19PvfkiMaYGGPMt8aYNwqo73k8HeZjjTF/zrVPTWBuVmwzrbW/OYhRRETKgNfJw1q7FJiFZxT5VmPM+8aYd4EdQGtgOZ4JEnOrh2dAYL6+DWvtF3imKKkOrDfG/J8xZgmwC0/n9gZgguNXJCIipc7RrLrW2pHGmCRgFJ4v+Owp2edSjCnZrbUvGGO+Bsbg6TPJnpL9ZQJoSnYRkYrGVITrlo0xB4ALa9euTYcOGmwuIuKNLVu2cPz4cYAfrbWNc2+rKMnjN1weJCgiEsCOW2vzzHBeUW5DuxvPHFm/AztdjqUo2SPhjwNbXI6lPND76Tt6L30rUN7PFngGbO8+d0OFSB7W2oCYmTfXSPgt1tpod6MJfHo/fUfvpW+Vh/fTyaW6IiIigJKHiIgUg5KHiIg4puQhIiKOKXmIiIhjSh4iIuKYkoeIiDim5CEiIo4peYiIiGMVYoR5AIkHEtCtd30lHr2fvhKP3ktfiifA388KMTGiiIj4lk5biYiIY0oeIiLimJKHi4wxlxpjHjTGLMi613umMcYaY/q6HVugMcZUNsZ0M8a8aIz5zBhz0BiTboz50Riz1BgT7XaMgcYYM9oYs8QYk2yMOWKMOWOM+dUY829jzGBjjHE7xkBmjHku6+/dGmMedjsep9Rh7q77gQfdDqKciAJWZf37EPAlcApoDfQB+hhjJllrn3QpvkA0FrgA+AZYj+f9bApcB3QD+hpjeju9/bSAMeYq4FHAAgGZhNXycNc3wFSgP56briS6G05AywTeAa611ja01t5sre1vrW0HDADOAk8YY7q6GmVgGQCEWWs7WmtvsdYOsNZ2BtoBPwO3AkNcjTAAGWOq4Lna6mfgPXejKT4lDxdZa2dbax+11i6x1u5yO55AZq39j7W2r7V2XQHb3sLzxwowuEwDC2DW2iRr7akC1m8DXs162r1soyoXnsbTIr4Pz50EA5KSh1QUm7OWjV2NovzIyFqedjWKAGOM+QswBlhorX3f7XhKQslDKopLspYHXY2iHDDGNMfzqxkgoL8Ay5IxpiowHzhKOejrVIe5lHvGmAhgaNbTd1wMJSAZY+7Gc0FCZTwtt6vx/PCcbK1d5mZsAeZZ4FJggLX2sNvBlJSSh5RrxphgYAFQG1gd6KcKXHINeTvGM4AngGnuhBN4jDFXA/8fsDyrDy7g6bSVlHev4bmsdD/qLC8Wa+1wa60BqgNtgH8CscBnxphGbsYWCIwx1YB5wAlgpMvh+IySh5RbxpjpwDA84z66WWsPuRxSQLPWplprt1trHwHGA+2BV1wOKxA8B/wJeMhaW2763DQxoh8xxiTgObfcz1q71OVwApox5kXgIeBXINpau93lkMoVY0w4cATPKazq1tozLofkt4wxe4BIIN9l5EBLoAHwA57W8U5r7fCyi6741Och5Y4x5gU8ieMI0F2Jo1T8hidxBAPheAa8SeGC8PwwLMxFWY86ZRNOyem0lZQrxpjngUeAY3gSx1cuh1ReXYsncfwGBPyVQ6XJWtvMWmsKeuC5dBfgkax1HdyM1QklDyk3jDGT8MzH9BuexLH5PLtIIYwxfzXGDMqaSuPcbdcAc7KezrHWni3b6MQf6LSVi4wxHYGZuVa1zlo+l3uWTWttpzINLAAZY3oBj2c93QmMLmTS12+ttc+XWWCB62I8Vwi9YozZhOeig1pZ67M/px/iuWRXKiAlD3eFAn8pYP0lBayTooXn+veVWY+CJAJKHueXCEwC/ornSqGr8cz+egjPQMsF1trl7oUnbtPVViIi4pj6PERExDElDxERcUzJQ0REHFPyEBERx5Q8RETEMSUPERFxTMlDREQcU/IQERHHlDxERMQxJQ8REXHs/wF87ezQG4/D0AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "######\n", "# 2. Fit an unregularize logistic regression model to predict `poorhealth` \n", "# from all predictors except income.\n", "# 2b. If you have time: use 'LogisticRegressionCV' to find a well-tuned L2 regularized model.\n", "# 3. Fit $k$-NN classification models with k=1,15,25 to predict `poorhealth` \n", "# from all predictors except income.\n", "######\n", "\n", "from sklearn.neighbors import KNeighborsClassifier\n", "\n", "# unregularized Logistic Regression\n", "logit = sk.linear_model.LogisticRegression(C=100000)\n", "logit.fit(X_train,y_train)\n", "\n", "# k-NN for k=1, 15, and 25\n", "knn1 = KNeighborsClassifier(1)\n", "knn1.fit(X_train,y_train)\n", "\n", "knn15 = KNeighborsClassifier(15)\n", "knn15.fit(X_train,y_train)\n", "\n", "knn25 = KNeighborsClassifier(25)\n", "knn25.fit(X_train,y_train)\n", "\n", "logit.predict_proba(X_train)[:,1],\n", "\n", "#visualize the predictions via boxplots\n", "plt.boxplot([logit.predict_proba(X_train)[:,1],knn1.predict_proba(X_train)[:,1],\n", " knn15.predict_proba(X_train)[:,1],knn25.predict_proba(X_train)[:,1]])\n", "plt.legend([\"1=logistic\",\"2=knn1\",\"3=knn15\",\"4=knn25\"]);" ] }, { "cell_type": "code", "execution_count": 109, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Classification accuracy for logistic were: \n", " Train = 0.9462659380692168 , Test = 0.9469214437367304\n", "Classification accuracy for knn1 were: \n", " Train = 1.0 , Test = 0.9193205944798302\n", "Classification accuracy for knn15 were: \n", " Train = 0.9462659380692168 , Test = 0.9469214437367304\n", "Classification accuracy for knn25 were: \n", " Train = 0.9462659380692168 , Test = 0.9469214437367304\n" ] } ], "source": [ "######\n", "# 4. Report classification accuracy on both train and test set for all models.\n", "######\n", "\n", "print(\"Classification accuracy for logistic were: \\n Train =\",\n", " logit.score(X_train,y_train),\", Test =\", logit.score(X_test,y_test))\n", "print(\"Classification accuracy for knn1 were: \\n Train =\",\n", " knn1.score(X_train,y_train),\", Test =\", knn1.score(X_test,y_test))\n", "print(\"Classification accuracy for knn15 were: \\n Train =\",\n", " knn15.score(X_train,y_train),\", Test =\", knn15.score(X_test,y_test))\n", "print(\"Classification accuracy for knn25 were: \\n Train =\",\n", " knn25.score(X_train,y_train),\", Test =\", knn25.score(X_test,y_test))\n", "\n", "# Note the severe overfitting of knn1, while the others are identical!" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "## Part 3: Evaluate Models via Confusion matrices and ROC Curves\n", "\n", "In this part we ask that you:\n", "1. Plot the histograms of predicted probabilities for your favorite model from above\n", "2. Create the confusion matrices for (a) the default threshold for classification and (b) a well-chosen threshold for classification to balance errors more equally.\n", "3. Make ROC curves to evaluate a model's overall useability.\n", "4. Use the ROC curves to select a threshold to balance the two types of errors." ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "As a reminder of Confustion Matrices:\n", "- the samples that are +ive and the classifier predicts as +ive are called True Positives (TP)\n", "- the samples that are -ive and the classifier predicts (wrongly) as +ive are called False Positives (FP)\n", "- the samples that are -ive and the classifier predicts as -ive are called True Negatives (TN)\n", "- the samples that are +ive and the classifier predicts as -ive are called False Negatives (FN)\n", "\n", "A classifier produces a confusion matrix which looks like this:\n", "\n", "![confusionmatrix](confusionmatrix_360.png)\n", "\n", "\n", "IMPORTANT NOTE: In sklearn, to obtain the confusion matrix in the form above, always have the observed `y` first, i.e.: use as `confusion_matrix(y_true, y_pred)`\n", "\n" ] }, { "cell_type": "code", "execution_count": 110, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAELCAYAAAD6AKALAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAT0UlEQVR4nO3de7BdVWHH8e/PQFAMhmqH4hjkIWir1YbyEOu0SYud8UGtDFI7UFsc7YwkqB0RiYP2YbUG0Wm1Cp0OpRkqaVV8TC1KFW0iDIhQjSBKB8UAZYRakZQAEiSrf+x9k8Ph3MM5N+s87r3fz8yZdc/ea6299uJyf9lnP05KKUiSVMsTJj0ASdLCYrBIkqoyWCRJVRkskqSqDBZJUlV7TXoAk5bkm8ChwHbgexMejiTNF4cDy4AflFKO7FyRxX65cZJ7geWTHockzVPbSin7dy5Y9EcsNEcqy5cvX87KlSsnPRZJmhe2bNnCtm3boPkb+igGS/Px1zNWrlzJpk2bJj0WSZoXVq9ezebNm6HHKQRP3kuSqjJYJElVGSySpKoMFklSVQaLJKkqg0WSVJXBIkmqyvtY9tAh6y6b2La3rn/FxLYtSbPxiEWSVJXBIkmqymCRJFVlsEiSqjJYJElVGSySpKoMFklSVQaLJKkqg0WSVJXBIkmqymCRJFVlsEiSqjJYJElVGSySpKoMFklSVQMFS5K9kxyf5INJvpbkh0l2JLkzyaVJVj9O+1OSXJlkW5LtSa5PsjZJ3+0neWmSLya5J8kDSb6d5Jwk+wyxj5KkMRr0iGUVcAXwVuBg4D+BzwD3ACcB/5Hk3b0aJvkocAlwNHAl8CXg2cBHgEuTLJml3duBLwC/BXwDuAw4AHgPsCnJvgOOXZI0RoMGy07gU8BvlFKeXko5oZTymlLK84HfBx4B3pXkNzsbJTkJWAPcBbygbXcicATwXeBE4IzujSU5GlgPPAC8uJTyklLKycBhwFeB44D3Dr+7kqRRGyhYSilfKaW8upRyZY91Hwc2tG//oGv1O9ry7FLKLR1t7gZOb9+u6/GR2DogwLmllGs72m0HXkcTdGuS7D/I+CVJ41Pr5P0323LFzIIkK4CjgB3AJ7sblFI2A3cCB9Icgcy0Wwq8rH17SY92twLXAEuBl9cZviSpllrBckRb/rBj2ZFteVMp5cFZ2l3XVRfgOcC+wD2llO8P0U6SNAX2OFiSHAic1r79VMeqQ9vytj7Nb++q2/nz7cyuVztJ0hTYa08aJ9kL+BiwHPhyKeVzHauXteX9fbrY3pb7VWjXOa7T2B12j2flgPUkSQPYo2AB/g44HriDx564T1uWIfuca7tOh9BcIi1JGrM5B0uSDwGvp7mU+PhSyl1dVe5ry2XMbmbdfR3L5tqu01Zgc5/2nVbSHHFJkiqYU7Ak+SDwZuBHNKFyS49qW9vy4D5dHdRVt/PnZw7ZbpdSygZ2XwLdV5JNeHQjSdUMffI+yftp7sD/MfDbpZTvzFJ15hLk5yV50ix1jumqC3Az8CDw1CTPmqXdsT3aSZKmwFDBkmQ9cBbwE5pQ+dZsdUspd9A8imUpcHKPvlbR3PdyF819KTPtdtA8ygXg1B7tDgNeRHN/zGXDjF+SNHoDfxSW5C+Bs4F7aUJlkKOF99HcHHlukqtLKd9r+zoAOL+ts76UsrOr3Xqax72cneTyUsrX23bLgItoAvH8Usq9g45/lC7c+7yxbu8ND5811u1J0jAGCpYkrwTe2b79HvCmJL2q3lxKWT/zppRyaZILaB7fcmOSK4CHaa4kewrwWZqHUT5KKeW6JOuAc4Grk3yFJtBW0TyI8lrgnIH2UJI0VoMesTy14+ej21cvm2mONnYppaxJchWwliYYltCcR7kIuKDH0cpMu/cnuQE4k+ZczBOBW4EPAx8opTw04NglSWM0ULAMc5XVLO03Ahvn0O5y4PK5bleSNH5+g6QkqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVXtNekBaP45ZN1lE9nu1vWvmMh2JQ3HIxZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqDBZJUlV7TXoAGt6Fe5/X/LDx4vFs8JSPj2c7khYEj1gkSVUZLJKkqgYOliTPSfKWJB9LcnOSnUlKklcP0PaUJFcm2ZZke5Lrk6xN0nf7SV6a5ItJ7knyQJJvJzknyT6DjluSNF7DnGM5HXjLsBtI8lFgDfBT4MvAw8DxwEeA45OcXEp5pEe7twPnAo8Am4CfAKuA9wAnJDm+lPLAsOORJI3WMB+FfRs4D3gNcDiw+fEaJDmJJlTuAl5QSjmhlHIicATwXeBE4Iwe7Y4G1gMPAC8upbyklHIycBjwVeA44L1DjF2SNCYDB0sp5cJSyttLKZ8opXx/wGbvaMuzSym3dPR1N80REMC6Hh+JrQMCnFtKubaj3XbgdcBOYE2S/QcdvyRpPEZ28j7JCuAoYAfwye71pZTNwJ3AgTRHIDPtlgIva99e0qPdrcA1wFLg5dUHLknaI6O8KuzItryplPLgLHWu66oL8BxgX+CePkdGvdpJkqbAKIPl0La8rU+d27vqdv58O7Pr1U6SNAVGeef9sra8v0+d7W25X4V2uyQ5DTit//B2WTlgPUnSAEYZLGnLMqZ2nQ6huTRZkjRmowyW+9pyWZ86M+vu61g213adtjLA5dCtlcDyAetKkh7HKINla1se3KfOQV11O39+5pDtdimlbAA29Gm/S5JNeHQjSdWM8uT9N9vyeUmeNEudY7rqAtwMPAg8NcmzZml3bI92kqQpMLJgKaXcAXyD5n6Tk7vXJ1kFrKC5K/+ajnY7gC+0b0/t0e4w4EU098dcVn3gkqQ9MuqnG7+vLc9NcvjMwiQHAOe3b9eXUnZ2tVtPc/L+7CTHdrRbBlxEM+7zSyn3jmzkkqQ5GfgcS5JfZXcYADy3Lf8qydtmFpZSjuv4+dIkF9A8vuXGJFew+yGUTwE+S/MwykcppVyXZB3NQyivTvIV4F6acyEHANcC5ww6dknS+Axz8v4pwAt7LD+iX6NSypokVwFraYJhCc15lIuAC3ocrcy0e3+SG4Azac7FPBG4Ffgw8IFSykNDjF2SNCYDB0spZRO77zEZSillI7BxDu0uBy6fyzYlSZPhN0hKkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqRvkNklooNr7mUW8v3PvukW3qDQ+fNbK+JY2HRyySpKoMFklSVQaLJKkqg0WSVJXBIkmqymCRJFVlsEiSqjJYJElVGSySpKoMFklSVQaLJKkqg0WSVJUPoZQGcMi6yyay3a3rXzGR7Up7wiMWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSarKYJEkVWWwSJKqMlgkSVUZLJKkqgwWSVJVBoskqSqDRZJUlcEiSapqr0kPQOp04d7nzb5y48X1N3jKx+v3KS1yHrFIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRV5Rd9SerpkHWXTWS7W9e/YiLbVT0esUiSqjJYJElVGSySpKo8x6LFbeNrBqp24d53j3ggs9h48dzbnvLxeuOQhuARiySpKoNFklSVwSJJqspgkSRVZbBIkqoyWCRJVRkskqSqvI9Fklo+H62OqT9iSXJKkiuTbEuyPcn1SdYmmfqxS9JiNNV/nJN8FLgEOBq4EvgS8GzgI8ClSZZMcHiSpB6mNliSnASsAe4CXlBKOaGUciJwBPBd4ETgjAkOUZLUwzSfY3lHW55dSrllZmEp5e4kpwObgHVJ/raUsnMSA5Sm2oDPQZvNsM9He8PDZ+3R9rRwTOURS5IVwFHADuCT3etLKZuBO4EDgePGOzpJUj/TesRyZFveVEp5cJY61wHPaOtePZZRSVqULtz7vNFuoPsp1vP8ydTTGiyHtuVtferc3lV3lySnAacNuK0XAWzZsoXVq1cP2GS3u279MWfmjqHb1fBzT146ke3+5P4dE9nupPYX3OdBLOfUKttd/a+T2+fls+zzmWMex3+te/7YtnXcYU+bU7stW7bM/Hh497ppDZZlbXl/nzrb23K/HusOAVYNs8Ft27axefPmYZrscsOcWknSbLaNbUub9/zfxcu6F0xrsKQtyxzbbwUGTYmjgCXAPcD3htzOSmA5zW/Blsepqz3jXI+X8z0+83WuD6cJlR90r5jWYLmvLR+ThB1m1t3XvaKUsgHYUHdIj5VkE82R0ZZSyupRb28xc67Hy/ken4U411N5VRjNEQfAwX3qHNRVV5I0BaY1WL7Zls9L8qRZ6hzTVVeSNAWmMlhKKXcA3wCWAid3r0+yClhBc1f+NeMdnSSpn6kMltb72vLcJLsuZ0tyAHB++3a9d91L0nSZ1pP3lFIuTXIBcDpwY5IrgIeB44GnAJ+leRilJGmKTG2wAJRS1iS5ClhLc9XEEuBm4CLgAo9WJGn6THWwAJRSNgIbJz0OSdJgpvkciyRpHjJYJElVTf1HYVNuA833wmyd6CgWhw041+O0Aed7XDawwOY6pcz1cVySJD2WH4VJkqoyWCRJVRksHZKckuTKJNuSbE9yfZK1SeY0T0lemuSLSe5J8kCSbyc5J8k+tcc+39Sa6yQHJTk9yT8kuSHJz5KUJG8b1djnoxrzneQJSX4tyXvavv47yY4kdyf5fJJXjXIf5ouKv9unJvmnJDcm+VGSh5P8JMlVSc5Isveo9mGPlVJ8NeeZPkrz/S8PAv8GfAb4v3bZp4ElQ/b39rbtz4ArgE8C/9MuuwbYd9L7vBDmGviTtl33622T3s9pedWab5rv35iZ3x8D/w78C/D1juX/SHvudjG+Kv9uXwU8AtwIfB74Z5rvmdrR8XfkyZPe555jn/QApuEFnNT+h/ohcETH8l8AvtOue8sQ/R0N7KT5BswXdixf1v5iFOCvJ73fC2Sufxf4G+C1wC8BFxsso5lv4FnAl4GXdv+BpHkyxva2v9dNer/n+1y37Y4F9u+xfAXw3ba/v5j0fvcc+6QHMA0v4Pr2P9If9li3quOX5QkD9ndp2+ZPe6w7rP1XyEO9fmkW+qv2XPfoY4PBMr757urvnW1/X570fi+CuX5t29/Vk97vXq9Ff44lyQqaryfeQfNx1aOUUjYDdwIHAscN0N9S4GXt20t69HcrzSHsUuDlcx74PFR7rtXfBOZ75ruRVlToa16ZwFz/rC1/WqGv6hZ9sABHtuVNpZQHZ6lzXVfdfp4D7AvcU0r5foX+FpLac63+xj3fR7TlDyv0Nd+Mba6T/DxwVvv2c3vS16h45z0c2pa39alze1fdQfq7vU+dYfpbSGrPtfob23wn2Rd4c/v2U3vS1zw1srlO8js052+WAE8HXgw8keZj36n86hCDpTmhDs2J9tlsb8v9JtDfQuLcjNc45/t8mj+Y3wH+fg/7mo9GOde/AvxR17IPAX9WSnl4yL7Gwo/CIG1Z69k2tftbSJyb8RrLfCd5F80fvm3A75VSHhrl9qbUyOa6lPKeUkqAfYBn01wk8XrgW0meW3t7NRgscF9bLutTZ2bdfX3qjKq/hcS5Ga+Rz3eStwLvpvnX+MtKKTfNpZ8FYORzXUrZUUq5pZTyXuA04GDg4iTp33L8DJbdTxQ9uE+dg7rqDtLfMyv1t5Bsbctac63+trblSOY7yZuAD9LcDHhCKeWaYftYQLa25bh+tz9Nc+PlUcAhFfqrymDZfYnk85I8aZY6x3TV7edmmv/RnprkWbPUOXaI/haS2nOt/kY230nWAh+mudz1le3ltIvZWH+3S3Mzy4/btwfsaX+1LfpgKaXcAXyD5r6Sk7vXJ1lFc13+XTT3nzxefzuAL7RvT+3R32HAi2iud79szgOfh2rPtfob1XwneSPN1UgPAa8qpVxRZcDz2Lh/t5McSnOkshO4dU/7q27Sd2hOwwt4Nbvvij28Y/kBwE30eBQDcAbN0cnFPfo7ht2PdDm2Y/kymi/0WcyPdKk61z3634B33o9svoE/bn+3fwq8fNL7N02vmnMNPBd4I7Bfj+38Mrvv8r900vvdcy4mPYBpedFcLjnz8LjP0XyGua1d9hke+2ykP2/XbZqlv86HUH4R+ARwd7vsayzuh1BWm2ua6/q/1vH6UVv3tq7lT5/0fs/3+QZWtqFSaJ5VtWGW1wcmvc8LYK5Xt8vvB75K8wDKT9McFc38N7gWeNqk97nXy/tYWqWUNUmuAtbSPNdnCc2/JC4CLiil7Byyv/cnuQE4k+YI5ok0h6wfpvkfbzFekglUn+t9gBf2WP5MHn0BxaL9qoKK870/uy+r/cX21cttwKL82oKKc30TzWXFv04zz0fR3Hf4vzQftX8C+Fgp5ZG6e1CHX00sSapq0Z+8lyTVZbBIkqoyWCRJVRkskqSqDBZJUlUGiySpKoNFklSVwSJJqspgkSRVZbBIkqr6fxqBCqeMfQK1AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#####\n", "# 1. Plot the histograms of predicted probabilities on test for your favorite \n", "# model from above\n", "#####\n", "\n", "# We plot them for the logistic and knn15\n", "plt.hist(knn15.predict_proba(X_test)[:,1])\n", "plt.hist(logit.predict_proba(X_test)[:,1],alpha=0.7);\n", "\n", "# Note this illustrates the fact that neither model predicted proabilities above 0.5 \n", "# and thus all prediced classification were 0 by default.. It also shows that knn15 \n", "# predictions are in increments of 1/15, while the logistic has predicted probabilties\n", "# in a more 'continuous' like range of values." ] }, { "cell_type": "code", "execution_count": 111, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 23)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m23\u001b[0m\n\u001b[0;31m print(confusion_matrix(y_test,t_repredict(knn15),0.06,X_test)))\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "#####\n", "# 2. Create the confusion matrices for (a) the default threshold for classification and \n", "# (b) a well-chosen threshold for classification to balance errors more equally.\n", "#####\n", "\n", "from sklearn.metrics import confusion_matrix\n", "\n", "# this function may help to manually make confusion table from a different threshold\n", "def t_repredict(est, t, xtest):\n", " probs = est.predict_proba(xtest)\n", " p0 = probs[:,0]\n", " p1 = probs[:,1]\n", " ypred = (p1 > t)*1\n", " return ypred\n", "\n", "# Using the logistic model throughout:\n", "\n", "# Re-calculating the default confusion matrix\n", "print(confusion_matrix(y_test,t_repredict(knn15,0.5,X_test)))\n", "\n", "#And then looking at smaller threshold values: 0.32 and 0.06\n", "print(confusion_matrix(y_test,t_repredict(knn15,0.32,X_test)))\n", "print(confusion_matrix(y_test,t_repredict(knn15),0.06,X_test)))\n" ] }, { "cell_type": "code", "execution_count": 118, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAFACAYAAADd6lTCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXhU1fnA8e+ZNclkJYEQDCKggCCLiLgrShG1tVLBulQqSn8FFVyrohWXWsEFQRZZ1FatW1VcsFq1uCBWUaQsgmyyBQMhEJLJOvuc3x93MmRnAslMlvfzPHkmc+6Zc98Qnsw75577HqW1RgghhBCiLqZYByCEEEKIlksSBSGEEELUSxIFIYQQQtRLEgUhhBBC1EsSBSGEEELUSxIFIYQQQtRLEgUhhBBC1EsSBSFEm6CUelEppev4KlVK/aiUmq+UOjHCsYaG+m9UShUrpVxKqV1KqTeVUlcopVQj4uqtlJqulFqplMpXSnmVUkVKqdVKqblKqdOP/KcWovkpKbgkhGgLlFIvAtcBPqCwshnI4NCHIi9wrdb6rXrGiAMWAb+v0uwOvS65StsqYIzWOqeBeKzAU8BNgDnUHASKgUTAWqX7UmC01rq0wR9SiBiQGQUhRFvzjda6c+grE4gDLgZ2ATbgBaVUx5ovCr2xf4yRJASB+UBfrXW81joFyARux3ijHwKsUEodV1cASikL8AEwGSNJeAM4B4jTWncA7MAJwN3APmAEkN4UP7wQTU0SBSFEm6a19mmtPwZ+F2pyAKPr6DoNOA8jSbhGa32z1npTlXH2a62fBs4EDgBZwGtKqbr+jj4CXAhoYJzW+iqt9X+11r7QWFprvU1r/STQE3ixKX5WIZqDJApCiPZiBVAW+r5v1QNKqS7AraGnC7TWb9Q3iNZ6I3Bz6OkZwG9qjJUF3BF6+ozW+qWGgtJaV2itrwd2R/JDCBFtkigIIdqTykWI5hrt12OsGQgAjx1ukNAah62hpxPqGMsG+IHpkQamtQ5G2leIaJJEQQjRXpyJcdkBYEeNY8NCj//TWudGON6S0ONZoTUJlc6vMtbeRkcpRAsjiYIQok1TSlmVUiOBV0JNPozFhVVVXopY14ihfwg9JgDdqrRX3oLZmLGEaLEsh+8ihBCtyplKqX2h72veHhkEJtYxa9Ah9HiwEecpqPJ9OrC9yvdw6BZNIVo1SRSEEG2NFeNWxpoKgYu01t830XkiLrokRGsmlx6EEG3Nl1prpbVWGDUUBgGLMWYNnldKpdXxmspP/42pZVC1b9XZg8pZiQ4I0QZIoiCEaLO01h6t9Trgt8AnwACMyos1VdZLGNiI4QeEHiuAqhUaj2QsIVosSRSEEG2eNmrV34Jx++MVSqnzanT5IvR4ilIqO8JhLws9flNZSKmOsbocUcBCtCCSKAgh2gWt9VYO3e3waI3DL2LcDWEGphxuLKXUFUCv0NOaMxQvYuwNYYlkrCpjypoH0SJJoiCEaE+eDD2epZQaVtmotd4DzA09vVEpdWV9A4R2oHwm9HQl8G7V46HaCU+Hnk5SSl3XUEBKqQSl1AtUv8VSiBZDEgUhRLuhtV4LfBp6en+Nw/cC/8X4u/iaUmqeUqpP5UGlVEel1K3AN0BHIB+4WmsdqONUfwY+w7gz4gWl1GtKqWqFmZRSxyul/oRxW+W4JvkBhWgGss20EKJNqLLN9Jda62EN9BsB/Cf09Eyt9Yoqx+KB54Frqrykrm2mV2NsM72zgfPYMGYW/kj1baadQBLVt5n+ACPpKEOIFkZmFIQQ7YrWeimwJvT0/hrHXFrr32Fs9rQI2IKxdsGGsWnT28BVwJCGkoTQWF6t9U3AScATwCqM2yiTMe6UWAPMCY11qSQJoqWSGQUhhBBC1EtmFIQQQghRL0kUhBBCCFEvSRSEEEIIUS9JFIQQQghRL9k9UrQISqk1QHegDNgW43CEaCtMoa+2XPVRh77qqmfRHhwPJAI7tdYnN8cJ5K4H0SIopZxASqzjEKJNUWaUxQom8+H7tlY6iA76we87fN+2rVhrndocA8uMQpQppXoDvwC6Agla61uqHDNjvFlqrXVRjEKMlTIgJSUlhUGDBsU6FiFapWAwiMvlwq+sBM12tMWK2WxFmdruhILWmmAggPb7MAW9mHwuEuLjsFjax9vb2rVrKS4uBuNvaLNoH/+SLYBSKhV4Dri8xqFbqnwfh7FFbbpSqpfWeke04msBtgHHDBo0iGXLlsU6FiFaHb/fz5YtW9hTUEy530RmVhcSk5JpL3tNuV0uDhbsx11WzDFpCfTudQLJycmHf2ErN2zYML788ktoxku2spgxCpRSdoySsZcDB4E3AVfNflrrcmAhxu9lzBGcp7dS6lal1CtKqc1KqaBSSiulGj1WjXGvUUp9pZQqVkqVKaVWKaVuVkrJ/x8hWogDBw6wv6iUcr+J43qeQFJySrtJEgDi4uM5pms3EpLT2FfsIjc3N9YhtRnyhz46bgKGAGuBflrrq4HSevouDj1ecgTnuRGjtvzvgN40wQImpdQzwKsY8X8FLMXYXncesDh0uUQIEWNFRUWUuH1kZnVpN9PudencJRu3X1NaVo7b7Y51OG2CJArRcRXGqtzJWusDh+n7I0Zt+ROP4DwbMLbRvRJjJeyXRzBGmFJqNEaSsw8YoLX+ldb6N8AJGJdIfgNMOppzCCGOntvtprSsHG8AEpPa/nR7Q5RSJCWnUObxU1TU3pZ6NY/2m3ZGVx+MN/9vD9dRax1USpUAjV69qrV+vurzJph2vDf0eI/W+qcq58lXSt0ILAOmKKXmaq2DR3syIcSRcblcuH1BHI6kdnW5oT6JickUHyjF5ap1hVccAUkUosMK+CJ5Mw1d90+ijjUM0aSUygZOwdhe962ax7XWXyql9gDHAKcD30Q3QiFEpUAgQEBrTObWMUlcUOZma34ZJS4fyfFWemUmkpEY12Tjmy1mAloTCLTt0gp5TheF5d5mP0/r+F/V+uUCCUqprAj6no2xpW2s73ioLNzxo9a6vqTl+xp9hRAxoLVGa43J1PL/pBeUufluRyEeX4DUBCseX4DvdhRSUNZ06wmUUuF/k7Yqz+li6cZ8glH4GWVGITo+A3oCfwQerq+TUsoKPI6xnuGj6IRWr+6hx5wG+uyu0bcapdQ4YFyE55PiCUJEUcH+fJ5/ZibLP/2E/fl5JCYlc9KgwVw7/iZOP/u8Ix63rLSEFxbM4dOP3icv92fscXH07tef344dz4W/vIyt+WUk2Mx4/EH2FrsJak2Fy8Nn776Ka+8Wfv5pE8WFBygrLsJqs9PpmG70G3oOF4y+jtSMThHF4K4ooyivmCKvmd3+tlnH7dste9i0egUub/PPmkiiEB0zgOuBe5VSPwMv1uyglDo71G8oRuGMudEMsA6JocfyBvpUFvhIquf4ccCR/8URQjSLrZs28Icrf42zqBAwFkA6Cw+y/NNP+Oqz/3DLPQ8w/ubbGz3uvrw9XD/mEvbsNj5fJDgSKS8rZeXXy1n59XKuuPYGTv3dXaQmWMNJAkDQXcqnzz8aHsdkNhOXkIirrITdP/3I7p9+ZNmSV7nprwvoM/iMJvgXaN2++3IpMx68h5KD++jSo0+zn08ShSjQWm9XSo0HXsIouvQkoTdipdQyjDscMjBuZ/QDv9da74tNtGGVK6KOZl5rF5HfeTEIKeEsRLNzu1zccsPVOIsK6XPSAKY9vYjje59IWWkJC59+gn88O4/Zjz3MiScN5MzzLoh4XK01d064jj27c+jS9Vgem/Mcg4achsft5rUXn2X29Id465W/Y+vcg6Ejr6g2ZR5UVgZf8jtOO/MsevQdRHKHjphMJvw+L5v+9w1vzpvGvt3bWfTgZB597XMS2umdHQf27WX+9Kn899MPD7Xt2dXs55VEIUq01q+GZhNmAwOrHDq3yvfrgUla66+iGlzdKus8JDbQp/JYnTUhtNYvUsfsSV1CCZPMPgjRzN569QX25v5MgiORuX//J5lZXQBjVuFPU/9Kbs5OPv/kQ2Y//nCjEoUvPvmQ9WtWYTKZePq5V+jTbwAA9rg4rp94Cwf25fHK3xbw0Ytz6HP2r/D4gtgsJrz+IMd16ciVTz1V54LGs3qNZuQZg/nlOSdTVlxE2U/fcf7oqxqMpbzMRr65mBOPzaBXr8xG/Ou0TH6/n3nz5jF16lTKyg5Vak5MSSM5NY29Oc27pK3lr3xpQ7TWy0O7ew3CqE/wCPAoRhnnU7XWA1tIkgDGbABAtwb6dK3RVwjRwv37XeMmpktGjQknCVWNm2hUld+0fh07t22NeNwP3zPGPf3sYeEkoarrJkxGKUVhwX6s+zdiNZso9/ixmk2c1qNDg3c9dD2uO8kpxh3jB/JjPdkaXStXrmTo0KHcfvvt1ZKEG264gRX/W0d6x+ZPhGRGIQa01j8AP8Q6jsNYE3rsp5SKr+fOh1Nr9BVCtGDlZaVsXL8WoN7ZggGDTyUpOZnSkhK++3o53Y/vFdHYq1b8t8FxM7O60LPXiWzbspHNq7/lrGtOCx873K2RO7f/REmxE4Bjujb02aXtcDqd/PnPf2bBggXV7t7o27cvCxYs4NxzjcnoDg5bs8ciMwpRoJS6Qyk1oRH9xyul7mjOmA5Ha/0zsBrjVs0rah5XSp0HZGNUbVwR3eiEEJEqKHPzzfYCPt6Qx3vLVobfdHr2qrv4q8lkoluPEwDY8dOWiM5xsOAARYUHGxwXoOcJvQHYvGkjm/NKWZ1TxOa80jpvjQwGgxzI38fH77/DpOt+C0DWMdmc94uLIoqptbvhhhuYP39++PcVFxfHtGnTWLNmTThJiBZJFKJjBvBAI/rfj7HgsdkppaaHNpCaXsfhyrbHlVLHV3lNJ2B+6OljUpVRiJapZs2Cg/sPTdt3yuxc7+sqjx3YH9k0f0GVfh0bGLfy2O7cvfgCQRx2C75AsFodhQfvmsyArqkM6taB4UP6cPfNN/Bzzk769OvPs68vIS4+PqKYWruHH34Ys9nYSufiiy/mxx9/5N5778Vma/4ZhJrk0kP0NHtdVaXUYA69gQP0DT1OU0r9qbJRa316lT5ZGBtI1SoGpbVerJRagLHZ1Hql1KcYpaiHA8nAexibQwkhWqCaNQvyDhaHj63Lq8BsqbuqX1nAeIPKK3Dy3c6Dhz3Ptm154e+3FLgpquc1haHTVVSUY7can1PtVhMJNjNb88vISIwjKSmZ9I6d8Hm94csNffr1Z8pfnqBb956H/6FbIY/Hg9lsrraZV//+/Zk2bRo9evRg9OjRMS3NLTMKLVMHjqyEczJwWpWvyvoGJ9Roj5jW+iaM3ShXY9yVMBJj3/NJwGitdduukSpEK1bi8hFvM3Ow3Bu6HbG5qvhVGTeCN7SqPUxKEW8zU+LyAfCnBx7li9Vb+e+GXXyzcTePz/sbJcVOxo2+mBmP3N/Eccfe559/zoABA5g7t3bpnLvvvpsxY8bEfP8OSRRaEKWURSn1B4w3+Ebf76K1Xqa1Vof7qvGacaH2cQ2M+5rW+iytdbLW2qG1PkVr/YxcchCiZUuOt+LyBsI1C6xxCeFjPm/9JZO9HuNzSlx8Qr19qrLHOw69toGtnSuPWezG5QOTUqQ7bLi8AZLjrbX6JyYlc/Flo3npnU9ITErmH8/O47OP/hVRTC3d/v37GTt2LMOHD2fr1q1MnTqVn3/+OdZh1UkuPTQDpdS9HNp5sVKn0K6Q9b4MiA89auCfzRSeEKKd6JWZGFqjYNQssCdlhI8dG+ele/f0Ol83u9So2Nir+7GcVk+fqk5I7sNfQt93trjqfc1il3EpoUPHTLJT44m3mXF5A1R4A/TPrr/eWmZWFy4Y+UveX/w67775CsMvvvSwMbVUwWCQ559/nnvuuQen0xluN5vNbNiwga5duzbw6tiQGYXmYcMoRlT5pTH+rRMb+HKE+pQAszD2fBBCiCOWkRjHaT06hGsWdOraIzyNvX3rpjpfEwwGydlh7CrfI3SXwuF0SM8grUN6g+MCbA/dRTGofz/sVjPOCh92q/mwdRQAOnU2aj7k5uyKKKaW6IcffuDss89mwoQJ1ZKEK6+8kk2bNnHxxRfHMLr6yYxC85gPLA59rzBqJhQA5zfwmiBQorXe08yxCSHakYzEOPpkHdqOpd+Ak9mwbjUrli/jFxf/ulb/9WtWUVpiTH6e1ojNoU498xz+88F7rPhqGb//46Rax/Pz9oaTiGHnD+fMnhm1+jRkz8+h/SMSHIfp2fKUlZXx8MMPM2vWrGpbX/fo0YP58+czcuTIGEZ3eDKj0Ay01ge01j+GvjZgLARcVaWtrq9NkiQIIZrbxaPGAPDv996qs8rhi4uMRXV9+w+ie88TIh73ksuMcVcs/5wtG9fXOv7yc8+gtaZjp84MPfOcasf8fn+DY+fs3M4Xnxj7Gwwe2ro2hfrxxx/p27cvM2bMCCcJVquV+++/nw0bNrT4JAEkUYgKrfUQrfUlsY5DCNGy5TldfLwhj9e+y+HjDXnkOY/k5qfqCsrc1YobXTDqarpkd6W8rJRJ11/J9q2bAaNq48xHHwgvFrzlntqlXwZ0TWVA11Tmz6xdduX8kb+k/8lDCAaD3PZ/17Ju9fcAeD0eXlo0l1f+tgCAG++8F2uNWgCPPXAPjz1wN2tXfYenymLIkmInS958lRvG/BK324UjMYlr/++mo/43iabu3btXu+1x2LBh/PDDDzzyyCPEt5KaEHLpQQghWoA8p4ulG/NJirOQkWin3ONn6cZ8RvTNJCv1yN5QKgsuVRY38vqDrN3r4uG5L3LXDVewaf06fjP8dBKTkqkoLyMYDKKU4pZ7HmjUhlAASimeWvRSeJvpsZeNIMGRiNfjDs8YXHHtDYy55rpar3W7Knh/8eu89sKzmEwmEpOS0WhKiw/VfcjolMmMBS/SOeuYI/q3iJWEhATmzZvHddddx1NPPcXYsWNjfrtjY0miEGVKqY4YtQy6YCxgrPd/jNZ6ZrTiEkLE1rpcJ0lxFly+ALsOVhAIalzeAC+vyGFwt7QGX+ssLGDX7iIC9gC7vYduadycV4ovEKxV3Cjg6MY7n67g+WdmsvzTT9ifn0dKWgf6DzqFa/9wE6c3Ym1CVZ2zjuGtj7/i7/Nn89nH/2Lvz7tJcCTSu19/rhw7ngt/NarO142/+XZ6nNCblV8vZ/euHRwsOIDf7yO9YyeO730i515wIaOuvJak5Ja9E/2KFSt45513ePLJ6oV1L7nkEnbu3EliYkOb8bZcqupmE6L5KKUygLnAaMB8uO6A1lofrl+bUbnN9HnnnceyZctiHI0Q0ffadzlkJNr5cW8JgaDxd1lrTYnbx3m9OjX4WmdhAbtycgjYk0nPPLQj5OqcIhx2S7gGkkkpenZ04KzwcdFJtYqxthnlZaXk5+7kxGMz6dUrsk2tjkZRURFTpkzh2WefBeC9997jsssua/bzgnEp48svvwT4Ums9rDnOIWsUokAplQh8CfwWqAA+x0gGPMAS4FvAG2orBN4G3olJsEKImOjgsFHu8YeTBAC3L0iSvXYhokgl2IzLDXD44kai8bTWvPLKK/Tu3TucJABMmTKFYLDt1KOTSw/RcStwIrAWGKm1PqCUCgJOrfXlAEqpVOAvGKWRc7XWt8csWiFE1A3MTmXpxnxc3gBxVhNuX5Bu6QkRrVE4cMBEmi4hYEsi65hDxY56dnTw3Y5CEmzmiIsbichs2bKFm266ic8//7xa+6WXXsrcuXMxmdrO5/C285O0bKMwii79SWt9oK4OWmun1voW4BngFqXUmGgGKISIrazUeEb0zcRuMVHi9mG3mI5qISMcKrjU2OJGon5ut5sHH3yQAQMGVEsSsrOzeffdd1myZAndunWLYYRNT2YUouMEjERheY32uvYL/StwE8aOjYvrOC6EaKOyUuOrLVxsbJJQ15qzjMQ4SQyayNKlS7npppvYtm1buM1sNnPrrbfy0EMPkZSU1MCrWy9JFKLDBhRrratWFang0O6OYVrrfKWUExgQreCEENGX53SxLtdJYbmXDg4bA7NTAWMBYqnHR5LdSr8uyRElCyaTCbNJ4W1D18WPRiAQwKxUk0////3vf6+WJAwdOpRFixYxaNCgJj1PSyOXHqJjD5CqlLLWaLMopaqVPlNK2YEUjP0fhBBtUGXNBJc3QEaiHZc3wOJVP7N4dS4ef5DkOCsef5ClG/MjKrpksVgwmxQ+rzcK0bd8Pq8Xs1lVK3TUFGbOnElycjIpKSksWLCAb775ps0nCSCJQrRUpqDHVWn7NvQ4vkbfiRi/l13NG5IQIlaq1kz4cW8J2w+Us/bnYtblOIm3mVFKEW8zkxRnYV2u87DjJSUl4bBbcbsrJFkASkqKcdgspKQc+aLNtWvXUlJSfcPfrKwsFi9ezObNm5k4cSJmc/u4g10Shej4N8atj1WrjSwKtd2llHpDKXWvUuo14CmM9QyvRD9MIUQ0FJZ7cdgt5Jd4wrdD+oNBfPrQhkFmk8Jht1BYfvg3fpPJREpKCg6bBWdRYbPF3Rp43G7crnIS46xHlCiUlpZyxx13cMopp/DAA7XLWI8YMYLOnTs3RaithiQK0fEu8AkQ/t+ltf4GeBgjWbgCYxHjVRi/k0+QbaaFaLPqqplgMZmwKuMTqtmkyEw2yjh3cNS15rm29PR00hxWCgv2tdtkwevxkLNzOxmJdtLS0hq1RkFrzTvvvMOJJ57IrFmzCAaDzJ07l9WrVzdjxK2DLGaMAq11LlBro3Gt9cNKqQ8xEoVsoBhYCrynpWSmEG1WXTUTBnVNAZMiOzUeh91CucdPqdvP6T3SDz8gkJqaSvdjuwI/s3fvboqdRSSnpJKYmITZYml1+wtEQmtNIBDA43ZRUuykpNhJhwQLx3RM49hjj414nF27djF58mQ++OCDau3Dhg0jOTm5qcNudSRRiDGt9SpgVazjEEJET2XNhJdX5FDiNu5wGDOkK2CsXygo89DBYeP0HumNukUyKysLpRQ28x7KPT5KD+SxP+9nAm34ZgiTApvFhMNmoWtqHB3SUujZs2dE6wd8Ph+zZs3i4YcfpqKiItzeqVMnZs6cyTXXXNMmE6zGkkShhVJKWWrcTimEaEPqq5lwNAWWADp37kxGRgZOp5OioiJcLheBQKDOGgttgdlsxmq1kpqaSlpaWsRbN3/99ddMnDiRDRs2VGufMGEC06dPJy2t4Y242hNJFFoYpZQN+CPwJ6rfJSGEaMHqqotQ802/ss+OA2XsKapg2/5y4qxmuqUnRFwzIRIWi4WMjAwyMjKaZLy2Zvny5Zx3XvUdMgcMGMDChQs544wzYhRVyyWLGVsIpVS8UuoOYCcwG+ga45CEEBGqqy5CzRoIlX32Ol1szitla34ZRS4vgWCQn/aXsXh1bkQ1E8TRO/vss8MJgcPhYMaMGfzvf/+TJKEeMqPQjJRSZwO/B/pibC29A3hZa/1xlT5WjI2g7gXSObSr5MtRD1gIcUSq1kXYdbCCQFDj8gZ4eUVO+PLC6pwiPP4ge4sr2FvkBgUJVisuX5AODjuFZR7W5TqbbFZBHBIIBKqtWTCZTCxatIiHHnqIWbNmNWrhY3skMwrNRCn1MMbW0uOBM4HTMG5//FAp9XSoT19gDTADyMC46+Ex4Dit9R9jEbcQovHqqosQZzVR6vGF+5R6fMRZTbi8QQJBjcWkMJvAEwhgt5rwBYIR1UwQkXO5XEydOpWzzz4bv7/6kq/+/fvz9ttvS5IQAZlRaAZKqXOBqaGnJRjJgAIGAcnAZKXUtxiXGDoCezGShee01uXRj1gIcTTqqovg9gVJsh+q2p5kt+L2BYm3Gfsy+IMahSLOYibeasJqNkVcM0Ec3ieffMLNN9/M9u3bAXjmmWe49dZbYxxV6yQzCs1jYujxa6Cn1vp8rfUwjF0kv8FIGv6BMYswHThea/20JAlCtE4Ds1MpdftxeY27C1zeAN3SExh7Rjd+0TeTX/TNZOwZ3eiWnsCpx3WgW3oCNrMi3maiZ8dEAkHokGgPbwwljlxeXh5XXXUVF110UThJAPj3v//dZu/8aG6SKDSP0zHKME/WWh+sbNRaHwAmh56agSe01n/WWrtjEKMQoolU1kWwW0yUuH3YLSZG9M2stt6gsk+X1Hj6ZCXRKzORzKQ4zCbFqcelMWZwtqxPOAqBQIB58+bRp08f3njjjXB7amoqixYt4qOPPpKaCEdILj00j0zAq7VeW8extYAXsGLs9yCEaAPqq4tQs48kA01v9erVTJgwgVWrqteuu/baa5kxYwaZmZkxiqxtkBmF5hEPFNV1IFSaufJYTtQiEkKINuiuu+7i1FNPrZYk9OrVi88++4yXX35ZkoQmIDMKMST7OQjRNuQ5XSzbks9HP+ajNE1eQEnUz2KxEAwaNartdjv33Xcf99xzD3a7PcaRtR0yoyCEEEchz+li8aqfWZXjxGYy7l6QAkrRM3XqVI477jh+8YtfsH79eh544AFJEpqYzCg0n05KqZJ6jiUANHAcjAmHxm+mLoSIqnW5TgorvJgU2K1GUR+llBRQamJer5dZs2Zx9dVXV6t9kJCQwDfffEPnzp1lsWIzkUSh+ZiAxMP0aei4XJYQohUoLPfiCxi3RBp3PiMFlJrYV199xcSJE9m4cSMrVqzgvffeq3Y8KysrRpG1D5IoNI/Jh+8ihGgLOjhsWM0Kj19jsxiJghRQahoFBQXcfffdvPDCC+G2JUuW8J///IcLL7wwhpG1L5IoNAOt9TOxjkEIER0Ds1P5aV8pm/PKcNiNSw92i0kKKB0FrTUvvvgid911FwcPhkvRkJiYyCOPPMIFF1wQw+jaH0kUhBDiKGSlxjNmSFf2l3rIKapAaTj1uDSG9e4k6xOOwMaNG7nxxhtZvnx5tfbLL7+c2bNnk52dHaPI2i9JFIQQ4ihlpcZzXu9O4ee/6Cv37jdWRUUFf/3rX3nyySerbeDUrVs35s2bx69+9asYRte+SaIghBANyHO6jDsbyr10cNgYmJ1aa6Ygz+lidU4RpR4fSXar1FA4AitWrGD69Onh5xaLhTvvvJOpU6ficDhiGJmQOgpCCFGPPKeLpRvzcXkDZCTacXkDLPln+6wAACAASURBVN2YX60+QmUfjz9IcpwVjz9Yq484vOHDh3PllVcCcNZZZ7F69Woee+wxSRJaAJlREEKIeqzLdZIUZ8HlC7DrYAWBoHEb5MsrcsL7OqzOKcLjDxJvMxYyxtvMJMVZpIZCAwKBADt37uT444+v1j5r1ixGjBjB9ddfj8kkn2NbCvlNCCFEPQrLvTjsFvJLPASCRmmTOKuJUo8v3KfU4yPOeuhPqdmkcNgtUkOhHqtWrWLo0KEMGzaM0tLSaseysrIYP368JAktjPw2hBCiHh0cNso9/nCSAOD2BUmyW8PPk+xW3D5jrwGzSZGZbKfc45caCjUUFxczefJkhg4dyurVq9mzZw8PPvhgrMMSEZBLD0IIUY+B2anhNQpxVhNuX5Bu6QmM6JsZvqzQr0sySzfmkxRnwWG3UO7xU+r2c3qP9BhH3zJorXnzzTe57bbb2LdvX7g9Li6OTp06NfBK0VLIjIIQQtQjKzWeEX0zsVtMlLh92C2maklC1T7xNjMFZR7ibeZafdqr7du3c/HFF3PVVVdVSxJGjhzJhg0bmDJlSgyjE5GSGYUYUEr1A7oCCVrrd2IdjxCiflmp8eGFi5XP6+ojicEhHo+HJ598kkcffRS32x1uz8rK4umnn+aKK66QDZxaEUkUokgpdTNwL1C5g4mmyu9AKZUKLAWswPla66KoBymEEEdBa83555/PihUrwm1KKSZNmsQjjzxCSopsitvayKWHKFFKPQfMAboA+4FDW82FaK2dwBqgPzAq2jEKIcTRUkrx+9//Pvx88ODBrFy5kjlz5kiS0EpJohAFSqnfAOOBA8AFWuss4GA93V/BSCAkURBCtHha61ptf/zjH7nwwguZM2cOK1euZMiQITGITDQVSRSiYwLGZYbbtdbLDtN3JRAEBjV3UEIIcTTWr1/Peeedx6pVq6q1m0wmPv74YyZPnozZbI5RdKKpSKIQHadgJArvHq6j1toNFANy35AQokUqLy/nnnvuYfDgwXz11VdMmDCBQCBQrY8sVmw7JFGIjiSgNJQERMKGsYZBCCFalA8++IB+/frxxBNPhHd5XL9+PStXroxxZKK5SKIQHQVAslIq8XAdlVK9AQewt9mjEkKICOXm5nL55Zdz6aWXkpOTE24/99xzWbt2LWeccUYMoxPNSRKF6Pg29Dgmgr73YVym+LL5whFCiMj4/X5mzZrFiSeeyLvvHrp6mp6ezgsvvMCyZcvo27dvDCMUzU3qKETHfOBy4K9Kqe+11j/W7KCUcgCPAmMxEoVnohuiEO1bntPFulwnheVeOjhsDMxOJSs1njyni9U5RZR6fCTZrfTrktxuiiutW7eOcePGsXbt2mrtN9xwA0888QTp6VKmuj2QGYUo0Fp/jvHG3wVYqZRagrFuAaXUQ0qpN4FcYHLoJdO11mvrHEwI0eTynK7wng4ZiXZc3gBLN+azbncRSzfm4/EHSY6z4vEHWboxnzynK9YhR4XZbGbDhg3h53379mX58uX87W9/kyShHZFEIXpuAaZi/JtfCiRg1EuYinFJIgXwAvdpre+PVZBCtEfrcp0kxVlw+QL8uLeE7QfKyTlYwfxl28k5WEG8zYxSinibmaQ4C+tynbEOOSpOOukk7rzzTuLj45k+fTpr1qzhnHPOiXVYIsrk0kOUaKMqyaOhCo1XAWdilHI2AfnACuA1rXVe7KIUon0qLPeSkWhn18GK8JbScVYTTpeX4zsdWoNsNikcdgsFZZ5YhdpsfvrpJ9avX8/ll19erX3q1KlMmDCB7t27xygyEWuSKESZ1no/RinnObGORQhh6OCwUe7xh5MEALcvSGq8DbcvSLzNjNmkyEy2U+7x08Fhi2G0Tcvj8fDYY48xffp0LBYLp556Kl27dg0fdzgckiS0c3LpQQjR7g3MTqXU7cflDaC1xuUN0C09gZuG9aRbegI9Ozro1yWZeKuZUrefgdmpsQ65SXz++ecMGDCAhx56CI/HQ3l5ObfddluswxItjCQKUaCU2qmUelQpdVKsYxFC1JaVGs+IvpnYLSZK3D7sFhMj+mYy8Ng0RvTNJN5mpqDMQ7zNzIi+ma3+rof8/HzGjh3L8OHD2bp1a7h9yJAh3HfffTGMTLREcukhOroBU4ApSqkfgVeBf2qtcxp+mRAiWrJS4xncLa3a88rH1p4YVAoGgzz33HNMmTIFp/PQgszk5GSmTZvGxIkTZW8GUYvMKETHBGAZRn2Ek4BpwA6l1FdKqRuVUhmxDE4I0fatW7eOs846i4kTJ1ZLEq666io2b97MzTffLEmCqJMkClGgtX5Oaz0cyAbuAL7HuDXyLGAesFcp9aFS6nehwktCCNFkgsEgv/3tb/n222/DbT179uTjjz/m9ddfJysrK4bRiZZOEoUo0lrv01o/rbU+HTgeeADYhHEJ6GLgH0C+UuqfSqlfxzBUIUQbYjKZmD17NgBWq5WpU6eyfv16Ro4cGePIRGsgiUKMaK13aK3/qrU+CRgIPA7kYBRi+i3wTizjE0K0Xvv378co3XLIRRddxKOPPsoPP/zAX/7yF+Lj28a6C9H8JFFoAbTW67XW9wLXA6tDzbKZuxCiUXw+HzNmzKB79+4sWbKk1vH77ruPPn36xCAy0ZpJohBjSqmTlVJPKKV2A58DJ4cOFcUwLCFEK7NixQqGDBnCXXfdRUVFBZMnT6a0tDTWYYk2QG6PjAGl1PHANcDVQK/KZsAFfAC8Bvw7NtEJIVqToqIipkyZwrPPPlutPTU1lX379pGUlBSjyERbIYlClCilugBXYiQIgyubAT/wGUZy8K7Wuiw2EQohWhOtNa+++ip33HEHBw4cCLcnJCTw0EMPcdttt2G1WmMYoWgrJFGIAqXU58C5GIlB5dqDFRjJwRta64JYxSaEaH22bNnCjTfeyBdffFGt/dJLL2Xu3Ll069YtRpGJtkgShegYFnrcALyOsUukVGUU7VKe08W6XCeF5V46OGwMzE5tEZUP85wuVucUUerxkWS30q9LcouIq6a3336ba665Bq/XG27Lzs5m7ty5XHbZZSgl66BF05LFjNHxODBAaz1Aaz1dkgTRXuU5XSzdmI/LGyAj0Y7LG2DpxnzynK4WEZfHHyQ5zorHH2wRcdXlrLPOCt/aaDabueOOO9i0aROjRo2SJEE0C5lRiILQrY9CtHvrcp0kxVlw+QLsOlhBIGjs1Pjyipxq+yxE2+qcIjx+YztpgHibmaQ4C+tynS1uVqFz585Mnz6dl156iYULFzJo0KBYhyTaOJlREEJETWG5F4fdQn6Jh0DQKAgUZzVR6vHFNK5Sj48466E/h2aTwmG3UFjubeBVzSsQCLBgwQKmTp1a69iECRP45ptvJEkQUSEzCk1MKXV56NsSrfWnNdoaRWst1RlFm9LBYaPc4w8nCQBuX5Ake2xX5yfZrbh9xoyC2aTITLZT7vHTwWGLSTxr1qxh4sSJrFy5EpPJxKhRozjllFPCx00m+YwnokcShaa3GGOXyC1A3xptjSVbuYk2ZWB2aniNQpzVhNsXpFt6AiP6ZsZ0ir9fl2SWbswnKc6Cw26h3OOn1O3n9B7pUY2jtLSUBx98kNmzZxMMBgFjQ6cnnniCN954I6qxCFFJEoWmtxojKdhVR5sQ7VpWajwj+mby8oocStzG3QWxThKqxrUu10lBmYcODhun90iPWlxaa9577z1uueUWcnNzw+02m417772XKVOmRCUOIeoiiUIT01oPiaRNiPYqKzW+2sLFWCcJlbJS42MSS05ODpMmTeKDDz6o1j58+HDmz59Pr1696nmlENEhiYIQotGq7kzY2LoIeU4X/9tVGK5X0DcrqcUkC9G8vdDv9zNz5kwefvhhKioqwu2dOnVi5syZXHPNNXK7o2gRJFGIgtBiRrfWOqL9G5RSIwGHLGYULUEgEMDpdFJUVER5eTl+vz+cKBSUevhux0ES7GbirRb2+vysXBngtB7pZCTZa41V2T+v2I3NYqLAH+S5vJ/q7R9tJpMJq9VKSkoKaWlpJCYmNtubtclkYsmSJdWShAkTJjB9+nTS0mJ3q6gQNUmiEB2LgTzgmAj7LwS6Ir8fEWMHDx5k585dlHv8lHn8VHj9VLlhgR9ynXj9QQpdPopdRgLh9QfZdqCC4zISao23q6ACfzCIzWLC5TcW65V5A3y1rYAB2anR+rEaZDUrHIWlJO7dR0piAieccAJ2e9MnMSaTiYULF3LyySfTr18/Fi5cyBlnnNHk5xHiaMkbUfQ09mOJzDmKmDp48CDbtu9gj9OFNc5BcnoGWUnJWKzW8Kfs3SqP1AQr2w+UYw/NMmgN5R4/3eoooHQwrgiH3ULlh3STUvTs6MBZ4aPvSVlR+9nqEwwG8XjclDid7HMWUeIuRuut9OrV66iSBa017777Lr/+9a+xWA792e3fvz+ff/45Z5xxhmzgJFosuRm3ZUoBPLEOQrRfFRUV4SQhLaMzx/U4ng7pGVhttmpT8cnxVlzeAMEqaxa8/iAJtro/gyTYLHhDMwkmpUh32HB5AyTHt4w3SZPJRHx8AplZXejZqw8+bOw5WML27duPeMxNmzZx/vnnM3r0aObNm1fr+LnnnitJgmjRZEahhVFKXQikAltjHYtov4qKiihx+0hISiWjU2a9/XplJvLdjkI8PuNygtcfJCsljtN6dCAjMa5W/54dHXy3o5AEm5l4mxmXN0CFN0D/7JTm/HGOiNls5tjuPdm2ZSMlpeW4XK7wHguRcLlcPProozzxxBP4fEblyalTpzJmzBiys7ObK2whmpwkCs1AKXUjcGON5nSl1A8NvQwjQeiCUXPh/WYKT4jDKioqotzjJyuz4YJDGYlGUvDhD/so9/hJsFnqTRKq9t+aX4azwkdyvJX+2Sn19o81s9lMUnIKZZ4yioqKIk4UPvnkE2666SZ27NhRbawbb7yR1NSWsRZDiEhJotA8OgEnVXmuAWuNtvoEgfeAh5o+LCEOz+PxUFbuwq9NxCc4Dts/IzGOPllJ1Z4frn9LTQzqkpScwoG9TkpKSujSpUuDffPy8rj99ttrVVE844wzWLhwIQMGDGjOUIVoFpIoNI/XgFWh7xXG7EAR8PsGXhMESoDNWuuDR3NypdQ1GDMaAzDKQG8GXgAWaK2DjRjnIeDBBrp4tNat5y++iIjf7yegg1itNrmPH7DZ7ASCGr/fX2+fyg2c/vznP1NSUhJuT0tL4/HHH2f8+PGyP4NotSRRaAZa65+AnyqfK6UKgTyt9YfNfW6l1DPATYAb+AzwAcOBecBwpdQVWutAI4ddB6ytoz22W/6JZhEMBglqMJllqxEwFjgGtQ7vvVCXZ555hltvvbVa29ixY5kxYwadOnVq7hCFaFaSKESB1jojGudRSo3GSBL2AeeGEhaUUpnAF8BvgEnA7EYO/Z7W+qEmDFWI1iOCWZU//OEPPP300+zcuZNevXqxYMECLrjggigEJ0Tzk0Shbbk39HhPZZIAoLXODy2wXAZMUUrNbcwlCCEiUXzwAB+9soCHv/+S/fl5JCYlc9KgwVw7/iZOP/u8Ix63rLSEFxbM4dOP3icv92fscXH07tef344dz4W/vKzO18yfOZ2Fsx6PaPxTzzibv735weE7hmitcblcJCQcKiiVkJDA/Pnz+f7777n77rubpUCTELEiiUITU0oNDn3r0lpvqtHWKFrr1Y04bzZwCuAF3qpjrC+VUnswqkOeDnxzJDEJUZfc7Zt56rZrKSsuAiAxKRln4UGWf/oJX332H2655wHG33x7o8fdl7eH68dcwp7dOQAkOBIpLytl5dfLWfn1cq649gamTp9Z63UJCYmkd6x/yj8YDFJ0sACAE08aGHE8O3bsYNKkSVitVpYsWVLt2EUXXcRFF10U8VhCtBaSKDS9VRh3OWwG+tVoawxN434/J4cef9Rau+rp8z1GonAyjUsUBiulHgfSgELgO+BDrbW3EWOINsrtcjHv3j9SVlzEsSf04+kFz3N87xMpKy1h4dNP8I9n5zH7sYc58aSBnHle5NPxWmvunHAde3bn0KXrsTw25zkGDTkNj9vNay8+y+zpD/HWK3/nxP4DGXPNddVeO27iZMZNnFzv2J999C9u/+NYAH59xTWHjcXv9/P888/z3HPP4Xa7AViyZAmXXVb3jIYQbYksw216haGv4jraGvNV1Mjzdg895jTQZ3eNvpG6FLgb+D/gHuAdYLtS6sjnk0Wb8darL3Bw3x7s8Q4mPfYsx/c+ETBmFf409a9cMPKXAMx+/OFGjfvFJx+yfs0qTCYTTz/3CoOGnAaAPS6O6yfewjXXTwBg/lPT8Hkbl7O+v/h1APqcNIBeJ/ZrsO+atWuYMmUKc+fODScJSinWrVvXqHMK0VrJjEITq2vhYpQWMyaGHssb6FMWekxqoE9V2zHWPXwE7ARsQH+MWybPA/6tlDpTa13nX0yl1DhgXITnGhRhP9HC/Ptd40rXab+4lLSOnWsdHzfxFj7/5EM2rV/Hzm1b6X58r4jG/fA9Y9zTzx5Gn3616w9cN2Eyr/59IQX78/nu6+Wcff4vIhq3qPAgX32xFIDLGphNcDqdzJ4zm/eXLEF7D03SDRo0iEWLFjF06NCIzidEayeJQttRuTS7sZc46qW1frmO5i+AL5RSi4HRwKPAr+oZ4jiMhEK0UeVlpWxcb9w522/ouXX2GTD4VJKSkyktKeG7r5dHnCisWvFfgHovVxj7MZzIti0bWflN5InCv997C7/Ph8Vq5eLLxtQ6HtRB/vWvfzF79mxKig9NDMbHxzNt2jQmTZpUbWMnIdo6+d/edpSGHhMb6FN5rLSBPpH6C0aiMEIpZdVa11VTYRfwZYTjDcLYDEu0YAVlbrbml1HiMsovB/K3okMbQhXZO7E5r5SeHR3VKi+aTCa69TiBDWv/x46ftkR0noMFBygqNOqO9ex1Yr39ep7Qm21bNkY8Lhy67HDO+SPokF59ss/jcTNp0iTWrFlTrX3oaafx0J03MWLEiIjPI0RbIYlCC6CUGovxydsOfKy1fvUIhtkVeuzWQJ+uNfoejc2hRxuQAeTV7KC1fhF4MZLBlFLLkNmHFq2gzB3e0Ck1wdg18qvVh/Yu65SZhS8Q5LsdhbX2e+iUaVySOLB/X2TnqtKvY2btyxk1j0U67k+bN7JpvXGlrK7LDnZ7HB2rFEjK6tKFP91xJ507ppKZWXvbbCHaA1nMGAVKqeuVUl6lVK0EQCn1Jsab6fXA74B/KKX+eQSnqfwI1E8pVd/ONafW6Hs0qu4WVFZvL9FmbM0vI8FmxuMPsv1AOblOF7v3H1pza7XbsVtNJNjMbM2v/l8iLt6oOeAqb2gJzSGuiopDr42rfyOmynErIhx3yVuvAZCa1oFzho+ss88dt99Bamoq48aN46033+LMs86KaGwh2ipJFKLjlxh7LlRLFJRSI4ExGOsLPsSofxAErlBKXdGYE2itfwZWY3zCr/Xa0B0K2RhVG1c0/keo5behxy1a66a4lCFauBKXj3ibmYPlXoKhyw3BYPUlMSaliLeZKXEdXXXvyssZQJPtNxEIBPh3aIHkJaPGUFhUyKOP/pXy8upJTUZGBkvef59JkyYTFydbmQghiUJ09A891nyDvg5j8eFTWutfa62vBO7ASBzGHcF5poceH1dKHV/ZqJTqBMwPPX2salVGpdQkpdRmpdQ/qg6klDpWKXWNUspeo12FLpVUnmvWEcQpWqHkeONyQ7DKm7gt7lB1woDPS7rDhssbIDneWu21bpcxQxDvOPxulAAJVfq5XBX19qscNyGCcb/58jMK9ucDYE9OZ8yYMbz77rssXLiwVl9HBLtmCtFeSKIQHZ2ACq11zdoIlcu0F1Rp+3vosdHVHLXWi0NjdQbWK6X+pZR6B2ODqr4Y21fPq/GyDKA3cGyN9g4YMyAHlFLfK6XeVkr9C+OWyX8A8cA8rfWixsYpWqdemYlUeAN4fEG0Bo8vSO8eXcPHO1CG3WKiwhugV2b1NbX78401BB071b/eoKqOmVnh7w/k17/+4EAjxn3/LWMRoz0hiVf/+Wb48sYbb77Jvn2RrXEQoj2SRCE6Emo2hD7xZwA5Wusdle1a63KMYksdjuREWuubMNY6rMZYHDgS2IaxGdToRuwc+TPwJPA/jMTjYmAExv+ZN4DhWuv6S9+JNicjMY7TenTAajZR7vFjNZsYPfy08KWBn7Zswm4111rIGAwGydlhbD3S44TeEZ2rQ3oGaR2MZTDbt26qt9/20N0Ohxs3b+8ePv3oXwD4TbZwe/fu3VmwYAGdO0eWwAjRHkmiEB0FQEJoF8dKlUXh/1tHfztQUkd7RLTWr2mtz9JaJ2utHVrrU7TWz9S1EZTW+iGttdJaD6vRflBrfbfW+nytdVetdYLWOk5rfZzW+iqt9edHGp9ovTIS4+iTlcTgbmn0yUri2MwM+g0wqoe7d63jzJ4Z1ZIEgPVrVlFaYvx3Pq0Rm0OdeuY5AKz4almdx/Pz9oaTiNPOqntcjeaT/3zCFZf9kkDAD4ApLhGb3c7Nkybx2uuvc8rgUyKOSYj2SBKF6Pg+9HgfgFIqBeMTvgaWVu0Y2twpAdgbzQBF25PndPHxhjxe+y6HjzfkkeesbwuQyBWUudmcV8rqnCI255VSUObm4lFG0aJ/v/dWnZcJXlw0F4C+/QfRvecJEZ/rklAxpBXLP2fLxvW1jr/83DNorenYqTNDQ0lFVT/n/szkSZP48333UXrQuHtX2RI485xzeevNN7l+3PVYLdZarxNCVCeJQnQ8g7FAcZJSKh9jz4VeQD6wuEbfynULa6MXnmhr8pwulm7Mx+UNkJFox+UNsHRj/lElC5V1FHyBIA67JVwz4YJRV9MluyvlZaVMuv5Ktm81SmyUl5Uy89EH+Cw05X/LPQ/UGnNA11QGdE1l/szptY6dP/KX9D95CMFgkNv+71rWrTbyba/Hw0uL5vLK34ylPTfeeS9Wm63W67///nu+/fZbtN+L9nkAGPuHicyZM4djjsk+4n8HIdobKbgUBVrrpUqpO4FpQMdQ88/ANXXs9Hh96PHTaMUn2p51uU6S4iy4fAF2HawgENS4vAFeXpHD4G4NFw6qKC9jx64iyrSVfHUw3L45rxRfIIjdany+qKyZkFMcYPbfXuP/rrqMTevX8Zvhp5OYlExFeRnBYBClFLfc80Cjdo4E47bIpxa9FN5meuxlI0hwJOL1uPH7jcsIV1x7Q62dIyuNGjWKf73/PmtXGMVBk5JTuOWOe1A0ze2WQrQXMqMQJVrrWUAWMByj8FEvrfXXVfsopazAImAs8H7UgxRtRmG5F4fdQn6Jh0Co1kGc1USp58jrG1R4/dgsh/5kVK2Z0Ltvf975dAXX3DCB7GOPw+v1kJLWgXOHj2TRa+8x/ubbj+icnbOO4a2Pv+IPk+6k+/G9CPj9JDgSOfXMc5ix4EWmTp8JQJGziJycXdVea1Im/nz//aQmGLMNF/16NDa7veYphBCHoaoWNhEiVipLOJ933nksW7YsxtG0fh9vyMPlDbD9wKGKhS5vALvFFNmMwrafKNNWOnc9tCN51RkFk1KkO2zYLSbsVjNn9ozGBqm1BXWQ95csYfacORzTpQsvvfQSZnPTTpT6fD52bP2RXllpDBhQexdLIWJp2LBhfPnllwBf1lyU3lTk0oMQbdDA7NTwGoU4qwm3L0i39ARG9M0kK7X+ksgApaUJpPkLcPosHNf9UKXunh0d4b0e4m1mXN4AFd4A/bNjs5fXtm3bmDZ9Gj+sM/Zu2FxSwltvvcVVV10dk3iEaKvk0kMUKaVMSqkxSqk3lVLblFLFoa9tobbRSin5nYijlpUaz4i+mdgtJkrcPuwWU0RJQkMq6yjYrWacFb46ayZEg8vtYs7cOfzud78LJwkAXY45huOO697AK4UQR0JmFKJEKdUV4w6HIZVNVQ4nAT0wtm3+Xin1W6317iiHKNqYrNT4apcZjiZJqJSRGBf1xKCq5V8t54nHH69WSdFisTD2979n/A3jZW8GIZqBJApRoJRyAJ8Bx2Ns+vQh8DmQG+qSDVyAsXnUUGCpUupkrXX9Re6FaCZKGfcFtKT1S/n5+5gxYwZffPFFtfaTBw/m3nvvpUf3Hs12bh0MolBNtjmVEK2NJArRcRtGkrAXGKW1XlVHn9lKqSHAklDf2zBupxQiqsxmMyaTwu87uh0gm0p5eRlXXX11uLojQEpKKrfffhu//NWvmv12R3/Aj9lkzFwI0R7J9fDoGI1RhXF8PUkCAKFj4zEuS4yJUmxCVGO324mzWgj4vfi83liHg8ORyOjLLw8/v+yyUbzzztv86leXRqUmQnlpKXFWM/HxR3/pRojWSFLk6DgecGmtPzlcR631x0qpCiDyWrdCNCGTyURqaiqOYhclxU7SO3aK6vkDwQBmk7la2x/+8H9s27aN68aN4+RBJ0c1npJiJxnxFtLSGr6tVIi2SmYUosOEsTYhUhqkfJyInbS0NFLiLRwsyMftOvo9IiKh0Xz88UeMvvzyWts+x8XF8fTTs6OeJBTszyfo95AYbyc5OTmq5xaipZBEITp2Ag6l1FmH66iUOgdwhF4jREykpqaSmZFORoKV3bu24SwqJBCIdIfyxtv9825uvvlm7r//fnJzc5kx48lmO1ckPG43+/bmUlSwj2NS4+natassZhTtllx6iI73gX7A35VSF2mt60wClFInAH/DmFFYEsX4hKhGKUXPnj0BMB04SPH+PeTt2U18QiIWi6XJ3jR9fh+LFy/mrbfeqrZ4cu2a//Hj+nWkpqY2yXkiFQwE8XjcBPweEu0WjkmN5/iePUhPTz/8i4VooyRRiI4ZwDiMtQoblFKvA8uAPYAd6AacD1yO8TvZG3qNEDFTmSykpKRQVFSEs7gEt89PUPsINuZCWj3WrF7N7DlzyM3N/YJezgAAIABJREFUDbeZlGLUqFFcf/31xCfYIBC9O4RV6PwpDoUjLoXU1FTS09NJSkqKWgxCtESSKESB1rpIKTUCeA8jWbieQ7tEVqWAnzBuoXRGMUQh6qSUomPHjnTs2BG/3095eTmBQIDgUWQK+/fv54EHHuDtt9+u1j5w4ECeeuopBg4ceLRhHzGz2YzVasXhcMilBiFCJFGIEq31RqXUAIwEYTQwGKicV3UCqzEqN76otXbHJkrR2uQ5XazLdVJY7qWDw8bA7NRwBcY8p4vVOUWUenwk2a3065J8VNUZLRYLKSlHt6/DG2+8wcSJE3E6D+XBycnJTJs2jYkTJ2I2mxt4tRAiFiRRiKJQArAg9IVSyhxqb75VYqLNynO6WLoxn6Q4CxmJdso9fpZuzGdE30wAlm7Mx+MPkhxnxe0Lho81RSnnI5WamlotSbjqqquYOXMmWVlZMYtJCNEwSRRiSBIEcTTW5TpJirPg8gXYdbCCQFDj8gZ4eUUOAB5/kHib8Qk93mYmKc7CulxnTBOFkSNHcuWVV7Jq1SqeeeYZRo4cGbNYhBCRkUShGSml0oBrgFOAZIxLDN8Br2uty2IZm2j9Csu9ZCTaw0kCQJzV2C0SIDnOGu5rNikcdgsFZZ6oxbdkyRJsNhsX/z975x0eVbH+8e9s3/RGII1ACIQqSECKoCAqQYIYQKkKghfFC4r+BNR7UUS9oNdrQxTFgoWmiIgiCtJUNFQBhYSahFQgm75997y/P87uspvshgR205jP85znZOdMeWdyds97Zt553xEjXNLfffddqNVq7umQw2kmcEXBRzDGRgBYA1FBcOZBAC8xxsYQ0d6Gl4zTUgjzV0BrtDiUBAAwmAUEKuWOv9UKKaQShtZB4tJEmL/C53Ll5OTgsccew+bNmxETE4OMjAyXnQNhYWE+l4HD4XgP7nDJBzDGEgB8BSAY4k6GCgAZALS2z60AfMMYi2g0ITnNnp6xIag0WKA3WUEkLjvEh/vh/gHxuH9APOLD/dChlT+6RQdBLZei0mBBz1jf+SUwm8147bXX0LVrV2zevBkAkJ+fjyVLlvisTQ6H43u4ouAbHgfgB6AIwF1EFEpE3YkoCOJSRAWAcAD/aEQZOc2cqBA17ujaGkqZuNyglEkcxor2a2qFFMVVRqgVUp8aMv7+++9ITk7GvHnzoNNd9n3wj3/8A0899ZRP2uRwOA0DX3rwDcMgeld8lIh+dL5AROtstgvLbfn46xbnqokKUaN3fKjLZ+e/fW24WFJSgmeeeQYffPCBS3qPHj2wYsUKDBw40Kftczgc38MVBd/QFmIQqB88XN8MUVGIbzCJOC2C6n4TooJUXvWVUFeICF988QX+7//+D5cuXXKk+/n5YdGiRZg7dy7kcnktNXA4nOYCX3rwDQEALhGRyd1FIsq3/enXcCJxmjt2vwl6kxURAUoUlOnxwS9ZKNEaEaSSw2gRfSUUlvk+2qNGo8GcOXNclIRRo0bhxIkTmDdvHlcSOJwWBFcUGhfuI5ZTZ5z9JhwvqMCB7BJUGS2o0FvAGHPxleBrIiIisHTpUgBAXFwcNm3ahM2bNyM+nk+ScTgtDb70wOE0E6r7TdCbBPgpJNBZLAB86yvh1KlT6NSpk0vazJkzodfr8Y9//AMBAQFeb5PD4TQNuKLgO8IZY8euIQ8RUeNFx+E0Oar7TVArJNAZBfgrZD7zlVBYWIgnn3wSX375Jfbt24c+ffo4rkkkEjzxxBNea4vD4TRNuKLgO+QAul9DHvKQzrlO6Rkb4rBRUMklCFLKIYEVfdqFIiZUDa3RgkqDBf0Twq+5LavVivfffx/PPPMMKioqAACPPPII9u3bxwM3cTjXGVxR8A3/a2wBOC0Pu2+Ez//IQYXBjPAAJSb1a43CCgOKq4wI81egf0L4Ne96+PPPP/HII49g//79LumdO3eGTqdz8bLI4XBaPlxR8AFENK+xZeC0TKr7TejZNhTeWp+qrKzEc889h7fffhuCIDjSO3bsiPfeew/Dhg3zUkscDqc5wRUFDqeRqe4boWdsiMdZgcIyvdf9JhARvvnmGzz22GPIz893pCsUCjz77LNYsGABVCrVNbXB4XCaL3x7JIfTiFT3jaA3WT36QrDnNVoEr/pNWLhwIcaOHeuiJAwbNgx//fUXnn/+ea4kcDjXOXxGgcNpRJx9I1ze9mjF53/kuCwxAMDhnFIYLWJESAAufhOuZVZhwoQJeOWVV2CxWBAZGYk33ngDEydOBGPczQeHw+GKAofTqFT3jQAAKrkY5Kk6lUYzglSXPR56y29C9+7dMW/ePJSWluI///kPQkNDr1yIw+FcN3BFgcNpRKr7RgAAg1lAoLKmC+RApRwGszijcDV+EzQaDRYsWIBu3brV8H/w8ssv8xkEDofjFq4ocDiNSHXfCAazgPhwP7chobtFB2H7iQsIVMngr5TV2W8CEeGzzz7DU089heLiYvj7+2PcuHGIi4tz5OFKAofD8QQ3ZuRwGhG7bwSlTFxuUMokbpUE57xqhRTFVUaoFVKPee1kZGRg6NChmDZtGoqLiwEAWq0W69at81mfOBxOy4LPKHA4jUx13wi1PfijQtR1MlzU6/V4+eWX8eqrr8Jsvmzv0LZtWyxbtgx33333tQnN4XCuG7iiwOG0MH788Uf885//xLlz5xxpUqkUTz75JJ5//nn4+/s3onQcDqe5wRWFRoAxFgggBoAfER1ubHk4vuVKDpW85USppKQEs2bNwpdffumSPnDgQKxYsQI9evS45r5wOJzrD26j0IAwxkYyxn4HUArgOID91a6HMMY22Q7+2tcCuJJDJW86UVKr1Th48KDjc2hoKFauXIlff/2VKwkcDueq4YpCA8EYWwhgM4D+zsnOeYioDIABwCgAqQ0nHcdXODtUOl5QgbOXtMjR6PD5Hzn4+cQFfP5HDnI0OqgVUjDGXJwo1Re1Wo3ly5cDAB544AFkZmbioYcegkTCv+YcDufqYUS1RzM+dOjQzQBuB9AHQCvw5Yp6U1lZGZibm9tJIpEIkZGR50NCQkpOnz7dw2q1yrt27XrIOW9FRUVQXl5ex4CAgNK2bdue81RnSyMrKytJr9cHBAYGIikpqbHF8RolWhNkEga92Qrnr5qVCIFKGSqNFkidtiYyBqjlUlgEqtU/gtVqRUlJCVq1alXjml6vh1p9bfEfOBxOw8IYg1KpRGBgIMLCwuqs4A8ZMgR79uwBgD1ENMQXsnl86B86dEgK4HHG2FSZTBYukUj8GWNyVHsL5lwZiUSiDAgIgFwut8hkskgAkUlJSTIigkql6uKc136jMMZCql9rySQkJKgEQYBUKm1sUbxKsFoOAhxul52R2GYQnGEAwFitXzKLxQKTyQQ/Pz9YLBbIZK5fY64kcDjNDyKCwWCAwWCAVqtFXFxck5kNrG12YLhEIpmuUCjiwsPDS4OCgvJUKpVJIpHUPgXBqcGRI0duACDt3r37calUKtjTLBaLvFu3bhnV8x8+fLgXETF311oqGRkZSVqtNsDPz69FzSiYLAIqDGaYLAIkDBAIUMgkCFLJoZBJHNeljEEiYRAEgpXIcd0Zo9GInJwcaLVaRxoRoVOnTi1OweJwrjcEQYBWq0VRURF0Oh1KSkoQERHR2GIBqF1RGCeTySJat25dHBERUdpgErVArFarTCKRWO1KAuf6wa4UlGiNsAoEqYS5KAH263qzBRarAJmUwV/uqiQIgoALFy6goKAAzkuFcrkcsbGxTeatg8PhXD0SiQSBgYEAgLy8PFRWVjYLRaGzRCLxDwkJKWgwaVooUqnUarFYZFarlUml0lpnZAwGg1wQBKlcLjc1lHwc36KQSeCnkLl8rn5dIXNvj1BZWYmcnBwYDAaX9MjISERHR9dYduBwOM0bu58To/Hagr15k9p+ZRQAJDKZjL8FXyNqtVpbWVkZXFZWFhweHl6rOXtRUVFrAPD3969qGOmub0wWwfY2T5BJGdRyWY0HuTfa0JksjhkFlVx6xTbMZjPy8vKg0Whc0v38/BAfH8+dJnE4LRR73JUrbTRoSPicZQMQHh5eDAAFBQUxJpPJo3JWUFAQWVxc3BoAIiIiihtKvusVu32AIAAyqQSCAIc9gdfbIDEstEB1ayM/P99FSZBIJIiLi0OXLl24ksDhtGCaYoA2Pm/ZAERERJSVlJSUVVRUhJw4caJrcHBwqSAIEgC4cOFCuMFgUJeXl4eYTCYlAISGhmpCQkIqG1fqlo/eLG5NFIhgMgsgAgQilGiNLksF14LOZIFA4g4HAJAwQMoY9GaLx+UGAIiJiUFpaSmsVitCQ0MRFxcHhaJu4aQ5HA7Hm3BFoYFITEw8l5WV1ba0tDRCo9FE2tNzc3PbOecLCwu71K5du9wGF/A6RFxukDiUBEB8kFsF70352Zcb7DAGSCQMFuvlGQWr1QoALjsX5HI52rZtC6lUipCQEK/Jw+FwOPWFKwoNhEQioQ4dOuRUVVVdLC4ujtBqtf4Wi0VOREwmk5n9/f2rIiIiNIGBgbrGlvV6QSYVtyM6LwXalwi8hX25QcJEJUFm2wIpk4ptlJeXIycnByEhIWjbtq1L2fDwcK/JweFwOFcLt1G4RmJiYnowxpKdD6VS2TsqKqrHXXfdlbBly5YA5/wBAQH6du3a5Xbr1i2zZ8+ef/Xq1etY9+7dM06dOlX64IMPRsbFxXVXq9U3+vv739i+fftukydPbrtv3746edApLS2VLFq0qPXAgQM7tWrV6ga5XN7b39//xsTExG733Xdf/KZNmwJ9MwrNE7VcBisRBCIA4lkhkyDMX4kgtdztEaiS4d23/ocBfXohMjQQwX4KtI1q5TF/mL8SCpkESpkEKrkUEsZgJYIMAs6ePYvTp0/DZDLh4sWLLv4R6sOQIUPAGMOqVau8Oj7egjHWIOuu06ZNa9BxmDVrFqRSKU6cONEg7V2v7Nu3D2lpaYiMjIRKpULHjh0xf/58lJeXX1O9P/30E9LS0hAVFQWFQoHWrVtj8ODBeO2112ott27dOgwfPhyRkZFQKpWIjo7GHXfcUeO+IyL06tULbdu2hU7XvN//uKLgJQYNGlQxZswYzZgxYzSDBg0qB4CtW7eGpqamJr3wwguRnsrp9Xo2ZsyYdikpKUlff/11uEKhoCFDhpTffPPNFRaLha1Zs6bVwIEDu86aNStGEDwbwK1fvz44ISGhxwsvvBB7+PDhgHbt2hlTUlJKBwwYUGGxWNhXX30VkZaW1mnEiBEJPuh+s8Tuw8C+3CBhcOvoyJnly5fjmWeeQX5+PkaOHImpU6di0qRJV25DAlisAhgDTJWlOJlxAqWll92TyGQymM1mr/avJbF7924wxjBkyJDGFgUAcPToUaxcuRJTpkxB165dG1ucFsvatWtx8803Y9OmTejUqRNGjx4Nk8mE//73v+jTpw8uXrxY7zoFQcDMmTORkpKCLVu2oFOnThg7diy6deuGkydPYsWKFW7LGQwGpKamYuLEifjtt9/Qs2dPjBkzBomJiThw4AC++OILl/yMMbz00kvIzc3Ff//736vqf1OBLz14iQULFhSlpqY6DBCNRiObMWNG3OrVq1u9/PLLsVOmTCnt0KGDy5NAEASkpqYm7Ny5MyQqKsq0cuXKrJEjR7psi1yzZk3wo48+2n7FihVt9Hq9ZNWqVTXsF1avXh38wAMPJBIRHn300aKXXnqpMDQ01EWrOHTokGrhwoXRWVlZKm/3vTlzJR8H1fnqq68c5zvuuKPObShkCmi1WuRk59R4uwgPD0dsbCzkcnk9pW8eZGQ0jIPRJUuW4Omnn0ZUVJTP25o3bx4EQcBzzz3n87auV/Ly8jBjxgwQETZt2oTRo0cDEF2YT5kyBevXr8fDDz+Mb775pl71Pvvss1i5ciX69euHL7/80mXJz2q14vDhw27LTZs2DVu2bEFqaio++eQTF2dIRqPR7cxSamoqevfujVdffRWzZs1CZKTHd8YmDZ9R8BFKpZJWrFiR6+/vL5jNZrZu3bp2WVlZcc7H4sWLO+3cuTMkICCAVq9eXdG1a9dQ5+sAMGnSpPLNmzefkslk9Omnn0Zu3LgxyLmdoqIi6SOPPNJeEAQ8//zzucuXL8+vriQAQHJysuGHH3449/rrr59vqDFoipgsAsr1JmiqjCjXm6AzWqAzWVBpMENnslxx22JurqindezYsc5tWiwW5OTkICMjw0VJUKlUSEpKQvv27VuskgAAnTt3RufOnX3eTlRUFDp37ozg4GCftnP8+HFs374dQ4YMQYcOHXza1vXMm2++Cb1ej6lTpzqUBECcffvggw8QFBSETZs21Wvp58SJE3jttdcQHh6O77//voZdkFQqRd++fWuU++mnn7B+/Xp06tQJGzZsqOExUalU4sYbb3Tb5vTp06HT6bBy5co6y9nU4IqCDwkICKB27doZAKCgoCBIo9FE2o/i4uLIlStXBgLA9OnTWWBgYITzdeedEbfccotu4sSJxQCwdOlSl9elV199tXVVVZU0KSlJv3DhwivOw40YMaJejpwqKiokzz33XOtevXp1DgwM7KVSqXrHxsb2GDFiRML69etdfpHtNhqe6rLbc5w8eVJRPb1r164BBQUF+PnnnzF06FCEhoaCMYZDhw4hJiYGjDEcPXrUo5zjxo0DY8wRZtkOEWHdunW48847ERERgUB/Nbp1SsSTcx7F2bNZuFRlhMVKV/RxYLcDyMrKAgC0b9/esfbuvDZJRPj8888xZMgQhIaGQqVSoX379njqqadQVFRkHyfExMSga9eu9gBgjjX8jz76CP369UNQUBAYYygrq3+46bpgNpvxzjvvONpSq9Xo0qULnn76aZSUlHgst3v3btx+++0ICgpCUFAQBg0ahG+//RbZ2dlgjKFdu3Y1yniyUSgoKMDs2bORmJgIlUoFPz8/tG3bFikpKfjggw8c+YYMGYKhQ4cCAPbs2eOor/pSxJVsFPbt24fJkycjPj4eSqUSERER6NOnD55//vkajq1q49133wUATJ061e31EydO4LnnnsPAgQMRHR0NhUKBVq1a4a677sKPP/7otsyqVavAGMO0adOg0Wjw2GOPoX379lAoFLjnnntc8ubm5uLxxx9HUlIS1Go1goKCcPPNN2PVqlVunfTk5ORgyZIlGDp0KOLi4qBUKhEWFoahQ4dizZo1de53Q7Np0yYAwOTJk2tcCwoKwqhRo1zy1YX33nsPVqsVM2bMqJd75HfeeQcAMHfuXCiVyjqXA4BJkyZBLpfj/fffR23Lx02ZJrf0kKPRyg9ml/oXa42yCH+lpU+7UG18uH+zXbytrKyUAkB0dHRFRESE3p7+119/yQoKCsIBYPLkyRqVSuVnMBjUUqnUEhoaWuNX66GHHir+/PPPWx06dCiguLhYGhERYQWArVu3hgDAhAkTNN72+X/q1CnF8OHDO2ZnZ6v8/PyE5OTkqqCgIGtBQYFiz549wRqNRj5+/PhrsyhyYvXq1fjyyy9x0003YcSIEcjNzYVcLsf999+PV155BatWrcIbb7xRo1xJSQm+++47KBQKTJw40ZFuNpsxYcIEbNy4EWq1Gr16J6NVq0hknDiOz1Z9jG83fYNPv/wWvXr3hgySWn0cpKSkoF27dtiwYQO0Wi3Gjh2LgADRTjUxMRGAqCRMmTIFa9asgVwux5AhQxAWFoa9e/diw4YN2L59Oz788EOMGjXK7Y/NnDlz8O677+Lmm29GamoqTp065RMjQIPBgBEjRmD37t3w8/PD0KFD4efnh19//RWvvPIK1q1bh507dyIhwdWc5YsvvsDUqVMhCAJ69+6NpKQkZGVl4Z577sFTTz1VLxkKCwuRnJyMoqIixMfHIyUlBUqlEvn5+UhPT0d2djZmzpwJQBx7lUqFn376Ca1bt0ZKSoqjnrrOVCxZsgT/+te/QETo1q0bBgwYgMrKSpw6dQqLFy/G0KFD62z/8O233wIAbr/9drfXX3/9dXz00Ufo0qULevbsiaCgIJw7dw5bt27F1q1b8b///Q9PPvmk27LFxcXo27cvysvLMXjwYPTp08dl98uuXbuQlpaG8vJyJCYmIiUlBVVVVUhPT8eDDz6InTt34rPPPnOp8/PPP8fChQvRoUMHdO7cGTfffDPy8vLw66+/Yvfu3di3bx/eeuutOvW9oaioqMDZs2cBwO0bvj199erV+PPPP+tc77Zt2wAAw4cPR1FREdauXYtTp05BrVajT58+GDNmDFQq19VZq9WKXbt2OcqdO3cO69evR3Z2tkNJGzVqlMfAbKGhoejduzf27duHw4cPo0+fPnWWt8lARG6PgwcPHvz777+1RHSwoY7s4qqj7+w8nf3ZH1lnNx7OPf3ZH1ln39l5Oju7uOpoQ8pRnyM6OtoIgL777ruT1a8dOHDgb6lUSjKZTDhz5oxLH954440sABQTE2O0p2k0mlOHDh2ynDx5sqR6XSaT6aDNnTZt3rz5pD3NFs2Ttm7dmunNflksloNdunTRAqBhw4aVXrx48U/n6yUlJYc3bdrk0mcAJN5StY9VZmbmMXfpUqmUVqxYQdXJyMggABQZGUlms7nG9eXLlxMAGjdunEv6ggULCADdcsstlJubS8WVBirTmaioXE8vv/o6AaD49gl0qrCUCsv0VFSupzKdiYorDTXasBMfH08AKCsry6McrVu3pr///tuRXlFRQRMnThTbi48ng8G1fvu4BQcH0759+zy27Ylbb72VANAnn3xSp/zz5s0jANS5c2fKy8tzpOt0OhozZgwBoP79+7uUycvLI39/fwJAH330kcu1jRs3ki2GCcXHx9doz+m+cPDCCy8QAHr44YdJEASXawaDgfbs2eOStmvXLgJAt956q8d+TZ061e04bNy4kQBQQEAAbd68uUa5/fv3U25ursd6nTl16hQBoLi4OI95du/e7fb+SE9Pp6CgIJLL5TXa++STTxzjdOedd1JFRUWN8gUFBRQaGkpSqZRWrVrlMm7nz5+nXr16ue3//v37Xe5H577ExcURAEpPT79Cz12xfw/qc7i7Nzxx9OhRAkAhISEe89j/r8nJyXWq02AwEGOMANA777xDgYGBNWSMi4ujgwcPupQ7efIkASCFQkHvvPMOKRSKGuW6d+9OZ8+e9dj23LlzCQC98sordZL1xIkTdOLEiTrltX//AewmD8/zaz2uaUbh2yP5Xl0MPJhT6m+yWCVq+eUoi3qzVfLhb1mKPvGhV7d3zAOje8V47U24OpcuXZLu2rXLf/78+W0FQcCSJUtyqxsyXrp0SQ4AERERjvSwsLAKi8Vy/vz58+0LCwtbRUVFXbJfk8vlCA4Otmo0GtmFCxdkAFBUVCSzT2VFRUV5ddZlzZo1IRkZGX7R0dGmTZs2nQsICHCZ0wwNDRVGjx7tVe+Ro0aNcvtW17lzZ/Tv3x/p6en44YcfcPfdd7tct081T5s2zZFWUlKCt99+GwEBAfjqq68QGRmJcr3J4Tdh+sxZ2LF9G3Zu/xG/7NiO1FGjavg4qC92y+bFixejW7dujvTAwEB8/PHH+O2335CTk4MNGza4nU6dP38+brrppqtqu67o9Xq89957AIC3334bMTExjmtqtRrvv/8+fvrpJ6Snp2Pv3r24+eabAYhLIlqtFsOGDcP06dNd6kxLS8PYsWPx5Zdf1lmOCxcuABBnC6rPmiiVStxyyy1X1T93vPDCCwDE/499utoZT2+s7jhy5AgAoEuXLh7z3HrrrW7T+/Xrh9mzZ+M///kPvv32W/zzn/+skcc+RW2PIujMm2++idLSUsyfP7/GskdcXBxWrlyJvn37YtmyZS7fBU/969ixIxYuXIiZM2diw4YN6Nevn8c+VWfcuHEoLq6fl/n6TPVXVYkrpLW5K7fP6FVW1u1nqLS01LE0M3fuXNx4441466230KNHD2RnZ+Nf//oXNm/ejLvuugvHjx93yGtfirNarZgzZw6GDx+OJUuWIDExERkZGZg7dy5+//13jBw5EkePHnXrQdW+M6Y+sx9NiSa19FBlMEuDVDKrc5pKJhEqDGb3czpNiFGjRnWqnqZQKOirr746PXbs2Iq61hMREVGam5vbTqPRuCgKQMMGCdm6dWsQAIwdO1ZTXUnwFfZ1aHdMmzYN6enpWLVqlYuikJGRgQMHDqBNmzYuU9K7du2CXq/HyJEjHZbGarnMFndB3Ao54OZB2Ln9R/z15wGMG3MPBIFgJYJ/PQ0LzWYzDhw4gOzsbEgkEgwfPrxGHpVKhcmTJ2Pp0qXYvXu3W0VhzJgx9Wr3ajh06BCqqqoce7+rExERgVGjRmHdunXYvXu3Q1HYs2cPAHjcCjpp0qR6KQo33XQT3n33XSxYsAAAcMcdd/gkhkVRURGOHj0KuVyOBx544Jrrs2/Hu5IzrMrKSmzZsgVHjhxBSUkJTCYxGOzp06cBAKdOnXJbrnfv3m7tPADghx9+AADce++9bq8nJycjICAAR44cgcFgcJlCNxgM+Omnn3DgwAFcunTJEZmwsLCwVnk8cSVfA9eK/bfOm0tvzvYBwcHB2L59u8PwtXv37vjmm29w44034tixY1i+fDmef/55l3JWqxVJSUn47rvvHFFb+/bti23btiExMRGZmZlYu3atW9uVsLAwAJcV5OZGk1IUAlRyq6HajILBIkgCVHJrbeWaAoMGDaqIjIw0ExEuXrwoP3jwYKDRaGQzZ85sn5SUlNm9e3eXmKEREREWACguLnZ5KkkkEmKMCfa4D3aMRiOrqKiQAkDr1q0tANCmTRuLRCKBIAgoLCyU9+zZ02txSfPy8pQA0LlzZ8OV8nqL2ra1TZgwAXPnzsWWLVug0WgcP9SffvopAGDKlCkua4Tnzp0DAGzZsuWKPzalGg0sVgEyKYO/vHY/Cs4QEYqLi5GXl+f4oY2IiEBZWRni4uJQ3WZzljBCAAAgAElEQVTEbiGfn5/vtr74+Pg6tXst2Ntu3769xzzu5LT/7UnG+sp+//33Y9u2bVizZg3S0tIglUrRvXt33HLLLZgwYQIGDhxYr/o8kZOTAwBo27Yt/Pz8rrk+u5OfoKAgj3m+/fZbTJ8+vVaj0IoK9+8OtY2j/Z6uywyIRqNxzBb98ccfuO+++5CXl1dveRoL+4yKfWbBHfZr7mZfaqsTAO67774au2MkEglmzpyJ2bNnY8eOHQ5Fwbnc9OnTa4R29/f3x5QpU/Daa69hx44dbhUF+/3iK+NkX3NNioK3p+97xYXovj9WGBykllkDlDKhymiRVOgt0tQbosqbukFjdT8KOTk58jvuuKPj6dOn1ZMmTWp/5MiRTOcHR79+/bQAkJ+fr8jPz5fFxMRYAECn0ykFQZBKpVIX5Sg9PV1tsViYRCJB//79dYA4TdmpUyd9ZmamOj093T8lJaVJh6a+ksVvbdbEwcHBSEtLw9q1a7FmzRrMmTMHgiA4nJw4T7UCl+MnJCUloX///i7XzE5xFuRSCQbdPADhAfWzZNbr9Th58qTjx8r+BiSRSNClS5caSoJzHk+o1XVywHlN1OVNrTY5PZWrryGtRCLB6tWr8cwzz+D777/H3r17sXfvXixbtgzLli3D9OnT8dFHH9WrzobAHnfD04M1Ly8PEydOhF6vx9NPP41JkyahXbt28Pf3h0QiwQcffICHH37Y4xjXdg/Y7+nx48fXMLirjv27pNPpkJaWhgsXLmDGjBmYNWsWEhMTERgYCIlEgm3btmH48OH1nq186qmnrmrpoa4zEfZZlbKyMlRUVLhVzOxblT3NwFQnMDAQ4eHh0Gg0HhVle7p9h1L1+utTzhn7/RIaGlonWZsaTWpGIT7c35x6Q1T5wexS/4uVRnmEv9JyS8dWTV5JcEd8fLx5/fr152666aauf/31l/+KFSvCHn30UccrRr9+/fTR0dGmgoICxYoVK8JffPHFCzqdTpmVlZUAAH5+fi4P/Q8//DACAHr37l3VqlUrhxKRkpJSlpmZqV63bl34okWLvDavFRsbawSAkydP1tlBk0wmI4vFwsrLyyXBwcEuWoHRaGR2u4yrZdq0aVi7di1WrVqFOXPmYPv27cjPz0efPn1cbAIAcc0WAHr06IEPPvwYerMFFqv4Y2gw24IwSZjDxXJ9OX36NKKjox2fY2NjAYhT055+dO1bK53tAhoau5x2WdzhTs7o6GicPHnS8YZenezs7KuSp3v37ujevTsAUZH84YcfMGnSJHz88ccYP3487rzzzquq1479DT03Nxd6vf6alTH7Mpan7ZTff/899Ho9xo4diyVLltS4fubMmatuOy4uDmfOnMHChQtr3O+e+OWXX3DhwgUkJyfjww8/9Jo8GzZs8HgveCI+Pr7OikJQUBA6dOiAs2fP4sCBAxg2bFiNPPv37wcAj/4L3NG7d29s377d4//PrvzY7R8AUcHo2LEjTp8+Xa9yztjLcYdLXiI+3N88Njm27OFbOhSPTY4ta45Kgp0bb7zRcP/9918CgCVLlsQfP3488eTJk4knT55MPHPmTOL06dOtAPDmm2/Gbt26tfuJEye66/V6P8YYRUVFFdjr+eWXX/zWrl0bAQALFiwodG5j3rx5FwMCAqwnT55UL168+Ip34Y8//uj+Tq5GSkpKBQBs2LAhXKfT1WmhMDIy0gwAR48eraFcfPPNN0FWq9WlHkEQmMlkkhGRI91qtUKn07k9Bg4ciJiYGBw+fBgHDhxwvHFOmjTJbV65XI7tP/+Mszm5MOgNMBkNqKjSQqvTw2IywmgwoLisEhWVWo9t2g+9Xo+LFy861prtMMYQFRWFYcOGISEhwWWWwxmz2ezYs96Ybojt69j5+fnYsWNHjesajQbfffcdAFc57caFa9eudVuvp/T6IJFIkJqa6nCu4+w3w24gZrFY6lVnmzZtcMMNN8BkMtXYNng19O7dGwA8OvmxLzfYFVVnjEYjvv7666tue8SIEQAuewetC7XJA+Cq/ShkZ2fX23K+vsqk/T5YvXp1jWsVFRWO+zQtLa3OddrtgNzd+wDw888/A0CNLYxXKmdP97T10X6/2O+f5kaTUxRaGi+99FKhv78/8vLyJOvXrw+urKx0HHfffbd60KBBqKqqwowZM5SHDx+GUqk0dOjQ4VRQUJAOEF04jx49uqPVamVTpky5NG7cOJc5z+joaMvy5cuzGWNYtGhR3OzZs2NKS0tr/F+PHTumHDVqVPsnnniibfVr7pg8eXJZ586d9QUFBYq0tLQEjUbjYlBaWloq+fbbb10WBwcNGlQBAIsWLYo2GAyOh//BgwdVTz31lKNdIoJOp1NqtVq1Tm9Q2+ezmVQGkkhhMFvdHiYr4b4Jop+Ed5a/i++//x4KhQJ3p42pkTc4LAIzHvoHysvKMG3KRJzIyECVwQyzRQCIYLRYUVpaitWfrcK58/ke23QcJgsqtXrn7iIgIABdu3ZFTEwMJBKJY2/8woULkZmZ6chntVoxf/585OTkID4+HuPGjavLv8AnqNVqPPLIIwCAxx9/3GHMBogGb7NmzUJVVRX69+/vMGQEgBkzZkCtVmP79u0OuxA7mzdvrtfDCwA+++wzt65yNRoN/vjjDwCu6/X22Y0zZ87UW1mwrzXPmzfPYRDozMGDB2tdv3cmISEBbdu2RV5entsydr8OX3/9tYvhmslkwpw5cxx2BlfDvHnzEBQUhP/85z9Yvny523FIT093+V/Y5dm5c6fLPSkIAhYvXoy9e/detTy+Zu7cuVCr1fj000+xefNmR7rFYsHDDz+MiooK3HPPPTVibezfv9+jN9CpU6ciNjYWBw8erDHjs2HDBqxevRpSqRSPPvqoy7XHH38cAQEB+Prrr2soLm+88QZ++eUX+Pv748EHH3TbF/s9XZvBdpPGk/bXGH4UmuNRmx8F+zFnzpwyABQbG2vOycnJKSgocBxnz57NscV3IADUoUMH/YgRI0qGDx9eEhsbawBAjDF66KGHiiwWi0c5vvjii9NBQUEWAKRUKoW+fftWpqamam6//fbShIQEvb3+1NRUTV37lpGRcaxt27YGAOTv72+95ZZbylJTUzW9e/euVKvV1r59+1Y65//777//CggIsACg6Oho4/Dhw0t69+5dKZfLhdGjR2vsY3XkyJGTpWXlJy+VlJ3TlFeejo6ONgOgzVu20omTZ6jKYPJ4HPn7uMv+5XvGjPGYt7RSS6PSxhJsPhq633AjjRh1D912Rwp17tbDsR96z77DtbZZZTBRld5E5VV6ioqOJgD0+x/ppNPpXPayC4Lg8JegUCho+PDhNGHCBEpISCAAFBoaSvv376fq2Ptytdj3USckJFC/fv08HocOHSIiIr1eT0OGDCHb/5VGjRpF9913H0VFRREAatu2rds94atWrXLsQ09OTqZJkybRgAEDCAA98cQTBIA6duxYp/6NHj2aIPoRoZEjR9LkyZMpJSXF4ath8ODBZDKZXMrceOONBIC6dOlCU6ZMoRkzZtCrr77quO7JjwLRZb8NAKhHjx40YcIEGjlyJCUmJhIA2rVrV53H+7HHHvPYjtlsdsgZGBhIo0aNonvvvZeio6PJ39/fUXbq1Kku5ex+FKqnV2fnzp0UFhZGACgqKopuv/12Gj9+PA0ePJiibffm+PHjXcqMGjWKbL8LlJKSQuPHj6eEhASSyWQ0f/78K/qnaEzWrFlDUqmUGGM0ePBgGj9+vMOHQ2JiIl24cKFGGbvPDU/fqT/++INsO7koKSmJxo4dS8nJyWT/rV22bJnbchs3biSZTEYAqFevXjRu3Djq2rWrY2w3btzotpxGoyG5XE5xcXFktVrr1O+m5keBKwoNoCiUlZUdDg8PNwOg119/Pdtdnh9++CEzLS2tOCYmxqhSqawqlcoaHx9vmDhx4qW9e/cer4ssxcXFfy5cuDC3X79+FWFhYWaZTCao1Wprhw4d9OPHj79Um4yejpKSksMLFizI69q1q9bPz8+qUqmsMTExxpEjR5Zs2LDhVPX8+/fv//u2224rDQwMtCiVSqFjx466F1988bzVanWM1cFDh89eKi0/W6k3ZVQZzCeio6NNDkXh1BnSGs21HgMGDnT8EHz9zbe15r1YoadP135Fd45IpcjWbUgul1NIaCh17NyFxk+aQp+sXk/5mooa5YrLKqm8Sl8jPa5tWwJAR479TeUVlaTT6cgZQRDos88+o8GDB1NwcDApFApq164dzZo1i86fP0/u8JaicKXD+WFoMpno7bffpr59+1JAQAAplUpKSkqi+fPnU3Fxsce2duzYQbfddhsFBARQQEAADRgwgDZs2EC//vorAaABAwbUqX+//PILPf7449S3b19q3bo1KRQKio6OpsGDB9NHH31UwykVEVFWVhbdd9991Lp1a4eDJ+cHXG2KAhHRr7/+Svfeey9FRUWRXC6niIgI6tu3Ly1atIg0Gk3tg+yE3QHYkCFD3F6vqKig+fPnU6dOnUipVFKbNm1owoQJlJmZ6VEhqKuiQERUWFhIzz77LPXs2dPxv4uPj6dbb72VlixZQmfOnHHJbzQaaenSpdStWzdSqVQUERFBo0aNovT09Do5smps0tPTafTo0RQREUEKhYI6dOhA8+bNo7KyMrf5r6QoEBGdO3eOpk+fTjExMY57YfTo0fTLL7/UKsuRI0cc96BcLqeoqCiaNGkSHTt2zGOZZcuWEQB66aWX6tZhanqKAiMPhleHDh06qFKpunTr1q1hQr+1YM6fPx8NAJGRkZdUKlWztbm4VgRBYFqtVm20CCqV2s9oezsFAORkZ8frDXo/Pz8/r24TtFgJWpMFZguBQDBZBMikEqjlUjAGCETwV8gcTpYEq4CLFy+itKwUgQGBiI2LdV8xAQaDHgqpBH5+6npb/bc0XnzxRTz33HOYPXs2li1b1tji+JyUlBRs27YNp0+f5oGhOLWSnJyMzMxMZGVl1dmY0R5xtTbHXnaGDBli93Oyh4iGXL2knrm+f90aiIsXL0ZdunSpjVKpvG6VBACwWCxSq0BSiVQqOCsJvkQmZfBXyGDf1aeQSaCWSwCITpeclYSKctG/fGlZKQCgsqoSVZUedpwyMdKcQFTvNfPmyvnz5906jPnhhx+wZMkSMMY8Bkpqabz66quQSCRYvHhxY4vCacJ8//33OHz4MObPn99sdzwATWx7ZEtFJpOZiUjiiwA/zQmr1SoViKRSmaxBHWjJpAwq+WVbTH+lq6NPk8mEosIiaHWuXsIDAgKgVHn2ryCVSmExWx3721s627Ztw8MPP4xevXohPj4eRISTJ0863n7+/e9/N8+AN1fBDTfcgJkzZ+L999/H/Pnz67xdkXP9QERYuHAhYmNjMW/evMYW55rgikIDoFartZWVlSFGo1F+Pc8qCILACIBEImkSsVaJCJpiDYo1xS6+D2QyGdq0aXNFj2+MSWC387keGDhwIKZOnYrffvsNO3bsgE6nQ1hYGO666y7MmjULqampjS1ig/Luu+86Qk5zONVhjDXb2A7V4YpCA9CmTZsLlZWVIbm5uTGJiYnZjS1Po0Lw2bSKxSpue7QKBKmEQSmTQiZlsFgJBrPVFuOBQSmTwGjQoaioqIZfhLCwMLRq1apONgeMMVwfKoJI165d8fHHHze2GBwOp4HhikIDEBwcXBUXF5edl5fXNjMzU9amTZsLgYGBWqlU2iTerBsSXz1Y7UaLEsYglTAIBGhNFqhkUhgsVhABEsZABJRV6XCp4LxLebVKjTZRba7oGpfD4XCuN7ii0AAcPnzY7mOUVVVVBZ85cyYYABhjQi12C9S7d+8jDSJgE0SwCtBqK2GweUXU63SwWkWjwcSkLlBUiwthtFhtigDBaCWAAIEAnckEuUQCiW2YrRYTzLoKSECwms2wqy4msuJiQT6CQkMREhIG3817cDgcTvOCKwoNgCAIbuexiUhyvaxv1xdtVSVyczzHI6iOfbnBriQAELc/CgRIBNg3+Bi1FagsuegoxxgDYwwWiwVVVZWoqqpEWYkGbdt1gETKNwVxOJyGpSk+E2pTFEwABIvFIpHJZNfdFLk36dKly/HGlqE5IpPJoFL7Qe3nB5lcjsK8XI957csNzmsbgiDAWFkKrUGLsCgx7LNCpUJwRGuo/fwREugPiS00tdViQVlpCS4WFUKn1aKoIA/RcXXyds3hcDhew64oNKVdcrUpCpmCICSUlZUFR0RElDaYRC0Qf39/Q2PL0NwIDApGoFO8eHM1o8PqKGVSaE0WCCTOJBh0WmjLih3LFfqKUoRFtII8IMjhZEkivfxFlMpkCG8VCcFqxaWLF1BeVoqo2Lgm9WXlcDgtH61W3KatVHremt3Q1Da3usFisRRfuHAhoqioKEKn0ykFQfDoyZFzmYyMjE6nTp1KaGw5mjX1fD7bHSuR1YKyS4Wo0BQ5lAQAsJhNsAg1nSxVR+XnB0DU6oXrxD8Ch8NpXIgIgiCgsrISRUVFAHDF7dkNSW0zCj8JgtDZaDROvXjxYnhxcXEsY0wObuZ1Reyhoo8fP35l/5vXEYIgSAhgjDFi1W4jk9mkgEAwVFUi68ypGmVJuKyg5uVkg0mq3YYEWKwW0Uuisy7LALlcDliMKMm/ss2D2Sy6uWCM4Xz2FSL9ESCQAAlj170LZw6H4z38/PwQFhbW2GI48KgoJCcnWwG8fujQoX1ms3kYgL4AWtVWhiOSkZGRLJVKLUlJSTxOhhN6vV5ltAhqhVJlqu7CubCwIN5oNPr5+/m7jfVgMplwOlMczsSkLlAqL29j1Ol1KCosgsHousITGhqKyMhIWAQGi+DZzEYQBJjNJpSXlaL44gUQESJaRaJ1VEyt/SEiGA16KOVS+NlmIjgcDudqYIxBqVQiMDAQYWFhTerl44oP/eTk5L0Amm7Q8iZInz59BADFRHR9+LOtI4yx7gG9RnQNuiktVx4a7WJ0ULTm3++TSZ/cOzkZH7z/QY2yF3Nz8MADDwAAvt21H+2j4wAAb739Fj7/7DOXvNHx7THpH3MxYdBNtcpzY7vwGu6XZTIZxk6ahvmLlogzEbVgsViQn/k3OkWF1Cl4C4fD4TRH+OwAp1kTG3s5uqNKpcKIe+/H0LvSIJVd+dYOb9UaVqsF2spKGAx6AMB990/H9EefuKKSwOFwONcLXFHgNGvS0tLw/XffITQsDPPnzUeOvu4P+J8PnAAgLiEUFeRh9cfvY83HK7Dlm6/wxgefo8+AQb4Sm8PhcJoNXFHwHcGMsWtxjE9ENONqCjLGJgGYBeAGAFIAmQA+AfAeEdXbJwZjLAXAkwD6AFABOAdgLYDXiMh4NTLWF5PJfTMSJsE7y5fD388fAJCTpal33YwxRMXE4amFLyE6Ng5Ln1uABbMfwne/HoKfrV4Oh8O5XuGKgu9QAZh6lWUZRNv9eisKjLHlAB4FYACwA4AZwDAA7wAYxhi7l4jqvO+PMTYfwCsArAB2AygFcCuAlwCkMsaGEZGuvnLWh/T0dLz4wvOOz9XtCvy9+DAfN2kaXn/5OVy6WIS9u3/GHXeN9lrdHA6H0xzhioLvMAP4oyEbZIyNhagkFAG4hYhO29JbA9gFIA3AbABv1bG+PgCWAtABuI2I9tnSAwBsAXALgJcBPOHdnogUFxfjf6//D9u3bQNZL0fn3vv7XiQm+cZ4UKFUIiQkDBcvFCI3J9snbXA4HE5zgisKvqOEiIY2cJvP2M4L7EoCABDRBcbYLIgzAk8zxpbVcQniaYizG6/YlQRbfVWMsQcBnAbwKGPsBSIq81YnrIIVGzduxDvL3oFWW1Xj+qBBvrMd0GmrUFJSDAB82YHD4XBQu2dGTjOCMRYLIBlijI6vql8noj0A8gG0AdC/DvUpAIywfVztpr5zEGdMFADuumrBq5FTVIIxE6bglaVLXZSEIbcOcfwtYe5v2+IqAzILK3E4pxSZhZUorqrpOdtisbgpeZkvPloBi83pUu+bBlxFDzgcDqdlwRWFloM9lPVxItJ7yHOgWt7aSALgB3Fm5KwX6qsTmvws5J8TPTOSYEV0VBu88spSPPTQZXONiopylJZoHIcgCCiuMmDfuRKYrQLeu78P/je+F15+cXENZSHttn5Y88n7yM3OconSlnX2NJY+vwDLX3sZADAsJRWdunTzVrc4HA6n2cKXHloO7W3nnFrynK+Wty71na8lT631McamAZjmlORvyP3bz3Qpx8AkEhfPjGQ23FCtNMyXspF9KRtPTHP193X/6DtcPncccCcsEhUEApw9Ox/ZuwOzHy1EkPrylsmcrLNY+twCLMUCMIkEUpkcgsUCQbhsIBncqg0qBDlmPjyzlq4DJAjQVZRCYSxDQEBArXk5HA7HFxw5csT+Z6Kv2uCKQsvB/qTS1pLHPpdfl2gj3qivHcQdEpexmkFmPeBh+eAydQ8+ln06E0xa039CafFFVBpclxpkIW0gmPQgkwEkWGAxGQEwQCqHRK6ERBUAvcQfR48eu2K7RAQyGWApK6izrBwOh+MjfBYcgisKLQf7u7S3wnt6o75sAHucPvtDKvdjcrW7GYVeIJKCMTOTq44BgDKu+1U1Wls5qUIN6VXV6gYiQLCqIe4KqU2hqi+9AAQDKAdw5Ap5Wzp8LC7Dx+IyfCwuMwCirZjPwt1yRcEHEFFj2H5U2s61zYHbr1XWksdr9RHRKgCr7J8ZY91Vcd09xHp4+n0y6ZOZXHWszaSlD9dBvkbHqq+Qlv+2JsFcnJNBREe9VS9jbDfEmZgjRDTEW/U2R/hYXIaPxWX4WFzGaSwO+aoNbszYcsi2nWuGXrxMXLW8damvrZfq43A4HE4zhCsKLYc/bedujDG1hzx9q+WtjUwAegBhjLEOHvLYwzPWpT4Oh8PhNEO4otBCIKJcAIchrlXdW/06Y+xWALEQvTZe0WMkEZkAbLV9nOymvgSIa2MmiF4a6yQmiAiCwK6ctRkgWJnLHksOh8NpgXBFoWWxxHZ+hTHm2CrDGIsE8K7t41Jnr4yMsdmMsUzG2Gdu6lsK0ZhxAWPsJqcyAQA+hnj/vFsPr4xWEixWwWxoEfedYNJLSLAKAGr34sThcDjNmBbxg80RIaINAN6D6H3xL8bYd4yxjRBdLXcFsAlicChnIiA6V6phi0BEByC6cfYD8DtjbBtj7EsAZyEaz+wD8K96iGgQDFUGS0lei/CNbC4+7ycYtQaIAbg4HA6nRcIVhRYGET0KcangMMSH+XAAZyAGgxpbn8iRtvpehejKeRdEG4dRAIoB/BvArfWMHFlq1uRVmEsK6uLHoclj1uQFWUoLKyBG1ORwOJwWCd8e2QIhojUA1tQx7yIAi66Q50cAP16zYECFpSSvyqotlejO7A/xS7zJa4GkGhp99pFAa2Wx0nTxXCXEvdwcDofTIuGKAqfBICJijBXoT6erwFg8mQ1SVbteZVJ1kM8chXgbwaiVGLKPBhty/wrXndl/HiQU1neWhsPhcJoTXFHgNChEdJExxvSn/oBgqIo0FpzsIA0I05PZKCerBbCY5FV/74xobDmrQ1azlIxauaWqxN+qLdUa8zJyrRUXs4mosLFl43A4HF/CFQVOg0NEFxhjJl3GLxpI5UGKiPgAq77iJ4COMav5gjbjl8YW0RUSAMFqFkw6nbk4t4DMhnIAGiLS+KjFVQB2gzuyAvhYOLMKfCzsrAIfCzur4OOxYHwbOKcxYYzJAAQBkAPeC8PgA6wQt0FWEJG5sYXhcDichoIrChwOh8PhcDzCt0dyOBwOh8PxCFcUOD6BMTaJMfYrY6ycMVbFGDvIGPsnY+yq7jnGWIrN4VMJY0zHGPubMfYvxpjS27J7G2+MBWNMwhgbyBh7yVZXHmPMxBi7wBj7gTF2jy/74C28fV9Uq3smY4xsR3XHYk0OH3xHpIyxhxljvzDGNIwxA2Ms1+Z4bZS35fcm3hwLxlgoY+w/jLG/GGNaxpiRMZbDGPucMdbLF/J7A8ZYEmPsccbYFzZvuYLtXh53jfVe+9gSET/44dUDwHKIrp/1AL4H8A2AClvaRgDSetY331bWAuBnAF8BuGhL+wOAX2P32ddjASDRVoYAaAD8BGAdgP1O6Z/AtpzYFA9v3xfV6o631SXY6nunsfvbkGMBIAyip1QCUAYx/so6AHsB6AB82Nh9boixgOhhNsdW9pKtvg0Qnc4RADNEx3ON3m83sr/p9F12PsY19tg2+uDwo2UdAMbabsJCAB2d0lsDOGG79ng96utj+/HXAujnlB4AYI+tvjcau9++HgsAHQDsAJBS/csN0QNnla2+Bxu73w1xX1Srm0FUIKsgWoA3aUXBB98RiU0hIAArAfhXux4AoHtj97uBxmKNrcwWOL1A2MZoke1aMQB5Y/fdjewPAXgVwH227/vua1EUvPr709iDw4+WdQA4aLsBH3Bz7VanG1dSx/o22Mo85+ZaAsTdCEYAIY3dd1+PxRXa+retvh2N3e+GHgsAs2zl5zg9DJqyouDt78jDtjK70YRnlBpoLAptZfq7uSaFOLtCALo2dt/r0JdrVRS8NrbcRoHjNRhjsQCSIYae/qr6dSLaAyAfYtCq/nWoTwExzgQArHZT3zmISw8KAHddteA+wNtjUQf+tJ1jvVCXV/HlWDDG2kN8C9uLmgHPmhw+GovZtvMrZHsKNAd8NBbGOuYrrmO+Zom3x5YrChxvcqPtfJyI9B7yHKiWtzaSIEauLCGis16oryHx9lhciY62c1P0FOmTsWCMMYjhzmUAZjSTh6RXx4Ix1gZAd4hr77sYYz0YY4sYY+/bDPruuHaRfYYv7gt7TJp/M8b87Im2e+U5AGoAm4noYn2FbWZ4dWy5Z0aON2lvO+fUkud8tbx1qe98LXnqU19D4u2x8IjtB/Ex28evr6UuH+GrsZgNYAiAp4no5FXI1Rh4eyxusJ2zATdUhM0AABEZSURBVCwE8AxEmw07zzDGfoFowNfU3qJ9cV/8G+KDbySAHMZYOsRZhp4QDV6/APBo/UVtdnh1bPmMAsebBNjO2lryVNnOdQk17e36GpKGlP1diF/2EwA+uMa6fIHXx4Ix1gHAEgCHALx29aI1ON4eizDbuT2AZwF8DqALRG+ntwHIAHALgC/rLanv8fp9YVOGbgPwKYAIAKkQjfoSAZwDsIeIKq9K2uaFV8eWKwocb2J/k/HWFLC362tIGkR2xthCAFMhhrq+j4jqukbbkHh1LJyWHBQAplPzit7p7fvC/hsug2jIOpWIMomokoh2AbgT4ta4oYyxW73Uprfw+neEMdYZor3OcAD3A4gCEAJgGMSH5krG2Mfeaq8J49Wx5YoCx5vYNfWAWvLYr9VFq/d2fQ2Jz2VnjD0JYDHEN4MRRHT8auppALw9Fo9BfEteQkTHrkWwRsBX3xHAzWwSEeVB3CoIiA/LpoRXx8IWN+ZriLMHY4joCyIqIqJyItoJ4A4AFwA8yBgbeg1yNwe8OrbcRoHjTbJt5/ha8sRVy1uX+tp6qb6GJNt29tZYuMAYmwPgfxDfFlOJ6I/61tGAZNvO3hqLNNv5Djdvye3seRhj3QFUEVFqHepsKLJtZ29/RwAgy0Mee3qbOtTXkGTbzt4ai34AugI45+77QEQljLGtAKYBuB3ArroK2gzJtp29MrZcUeB4E/sWvW6MMbUHa9u+1fLWRibEB2EYY6yDh50PN9WjvobE22PhgDH2TwBvAzAAuNu21akp46uxGFDLtWjbUV6P+hoCX3xHtAD8AYR7yBNhO1d5uN5YeHss7C8Utf3Py2znsFrytAS8OrZ86YHjNYgoF8BhiGvH91a/bnv7iwVQBNH/wZXqMwHYavs42U19CRAfFiZcnl5tEnh7LJzKPQLRX4ARwD1E9LNXBPYhPrgvhhARc3cAeMGWbbktLcR7Pbl2fDAWZoiueQE3SwuMMTnEZRpAdMDTZPDBd6TAdu7MGPP0f7f7DPA0+9Ii8PrYNrb3KX60rAPAOFz2+JXolB4J4DjcuA2FuM0tE8Bnburri8sunG9ySg/AZc9lTdWFs7fH4h+2sTAAuKux+9eYY1FLO4vQ9D0zevu+6AnRQ6kOwDCndCmA12315QFQN3bffTkWEB+K+bYyXwMIcromwWXvpWYAHRq773UYG/vvm0fPjBB3/mRCtNe55rH12E5jDwY/Wt4BcbseQVw2+A5i8JFyW9o3qBmrwP7jvttDfc5BobZB3Op1wZaWjqYdFMorYwGgFy4HPMqAGNPA3fFaY/e5oe4LD23YyzRZRcEXYwHRfbUAUWFIh+j6/CwuB4ka0Nh9boixgGiwaHfTXAxxRnIjxK2RZBuffzZ2nz2MQ2/b/85+2IM3nXJOr1ZmlS3PKm+MraeD2yhwvA4RPcoY+w3APyH6FJdC1Ho/BvAeEQn1rO9VxtgxAP8HcYZBBfGL/zbEB2NT3BIIwKtjEYLLW5462w535AB46uol9h3evi+aMz74jixjjP0F8X/fH+JDpxDiToglRJTtRfG9ijfHgoi2M8Z6AngSoj+FIbb6iiBG03yLiNK92wOvEQTRILM6Hd2k1QlvjS2zaR0cDofD4XA4NeDGjBwOh8PhcDzCFQUOh8PhcDge4YoCh8PhcDgcj3BFgcPhcDgcjke4osDhcDgcDscjXFHgcDgcDofjEa4ocDgcDofD8QhXFDicFgpjrDNjjBhjhsaWpSXAGCuyjWf/K+d2Wz7dVn6Ct2XjcHwJVxQ4nAaCMbbK9qC40jG3sWVtSBhjSz2MQyVj7ARj7B3GWGJjy+kJxlgiY2wRY2x2Y8viLRhjj3j4n+gYY1mMsXWMsdt81HaEbTz/7Yv6OfWHu3DmcBoeM4CSWq5rG0qQJoYVon9+O60AdLEd0xljE4hoc6NIJnIaYtyE6iF7EwE8D+AkxMiensiG6Iq7whfC+ZALTn+HAmhnO8YzxpYQ0bNebi8C4ngaAbzk5bo5VwFXFDichud3IhrS2EI0Qc4QkSOGBWNMCWAExMA2UQBWM8YSiOhSYwhHRIOvsXxzXHIwElEb+wfGmATADRDjrAwG8Axj7Gci2tlYAnJ8D1964HA4TRIiMhLRJgBTbUkBAKY0okjXPUQkENERAGkQZ1cA4IFGFInTAHBFgcNp4jDGYhljsxljWxljZ2zrxBWMsUOMsYWMsaCrrHcsY+xHxthFxpiZMaZhjGUyxlYzxsbVUi6NMfY9Y+wCY8xkO3/LGBt29b30DBFtB6CxfUx2I08IY+xFxthfjLEqm23DEcbYc4yxwFr6MYwx9g1jrMDWjzLG2GnG2EbG2ENu8tcwZmSMFUEMZQwASW7W9Cc45a1hzMgYe8mW9lttY8AYe9CWL8/2Vl/9+hDG2JeMsXxbXzSMsW21/R+vBSLSADhk+9jVg8z1vm8ZY+kQw6gDgNLNeD7tpkwCY+xd2/9Ob2vjAGPsKcaY2lt9vq5p7Bjc/ODH9XLgcuz43fUs972tnP0ohbieb/+cCaCNm3KdbdcNbq79r1qdFRDX3u2fs92UUQL4slq58mqfF1/FuCy196OWPEdteTa76WOuU/tVAHROn88CaO+mvjnV5NbayjrS3JQpsl3r75R2DKK9CQGw2PI4H6Od8qbb8k1wSutiSxMAxNfS/222fK9VS2cA3nTzPxGcPq+CLVJwPf4nj3i6d5zy7LDlOeSt+xbAFgCXnPJUH8/Z1fJPAGCo9n80OX0+DCC8Mb/3LeHgMwocTtMnE8AzEB8qaiL6//bOP9bLqo7jrw/ED7kgXhQQBkE6w1YwYW4tXCizNTYqRv4gsFbakgzLzKY2l5rNOXIrbFkyxUgMFdBkLZLcYvlj1GxiaUFmkZlokdEPo0HRuz8+59zvw3Of53u/33svl9v4vLazZ8/59ZzzPOd+z+ec8/l8bidwHHAusAOYAdzeamVmNgO4Mt1+AThJ0vHAKGAisATYWlF0FXAB8AJwITBa0ljgeOByfKL9vJktbreDLfDGdM3b3ZjZSOA7wBRgNzBf0migA1gA7AFOAR40s2GFcmOBL6XbO4ApkjpS2ZOAhbhA1COSZgHL0u0Lkk4uhc09lN+JC0GGT3rdMLMJQLYwWF9Kvhq4Alc4vBQ4IX2TjtSuP+FHN/1qSWNmJwJz0u1va7K1PW4lLcR1HyDpR5RCl7KomZ0F3AsMBW4GpkrqwMfxWbiQMBu4u6/9PeY52pJKhAjHSqCxo3CQ7iulHL7ZZp3jaazUJpfSKncU8DNlATvaeM7bUpm9+A9yVZ5c70/b7EPTHQXgfBorxBWF+I/l/gFvrig3G1/lC1hWiJ+X4v5CGyttKnYUUvyCZu0v5Ou2o5Dir07xP6spd3lV/bhQsz/1/4yasrmvfwSGttHXyh0F/Lj6DOCxwjdZ1Gq9fRm3FXU8lfJdVpN+Uuq3gJnttjFCI8SOQhAMPMPwlXtV6GynIrkFwE/wH/BWHQFl87zOtCpvhaxQuEnSSzV5NuA//HPMbFyL9daSzriXA3elqH3AtwtZ8vn7RknPl8tL2gFkc8oLC0m5/yOAPrezH7gPn8xmmVnVef/SdC3vJizBV+jb5AqG3ZD0GPAyMAG3VmiXEUk349Wkj/EvfDcgr/rvoPGOW6aX47YLM3srcCb+Le+qyiPpz8Cj6fZd7T4jaBDmkUEw8PxIbZpHmtlcfGt5LjAZ31ouM7nF6p7Ez7GnAU+a2TeARyW92KTM3HS9qIejBUthCs19RVQxw8xUk7YPeL+kvxbi8tb3tiZ1/hDX0J9TiPsl7tNgOrDdzL4GPFIlbAwEkl5KyozvxI8LuhwNmdk04B3p9r5S0fxN5qVJvI4sDE3FJ/l2mVgRdwi4RNI9zQr287gtkvs+CnjJzOryjU7Xqb14RpAIQSEIBjlmdh2HO575Dz5xHkz3Y4GRVP8Id0PSXjP7MHAPPoHemZ6zB1eaWyOprIU/KV3HpNATo1ppS4miwyXh2+q/xwWB1ZK6HP+Yzwx5Any5SZ1/SNfxOULSQTNbBjwEnAbclup8DVfQ+5akLb1of19YjwsKSykICunegKck/bpUJn+TUbT2vnvzTQ5IGglgZkNxXZFLgWuAr5jZM5J+XlWwv8dtidz3N1AtyJTpTd+DRBw9BMEgxszmAF9Mt1/GFcNGSBqnpOCFa5eDTygtIVeym46fRW/Cz98nAx8BHjezr5aK5N+K5ZKshfDjXnS3qAw4SdKpkuZLuqkoJFT0dUS7D5K0HTgV16u4F99hOBE/oviemT1cZYZ4BNmIe+w8xczeXoivO3aAxje5pcVvcn9fGijpkKTdkj6HCwDjgA1VJohHatwWyH3f3mLfP96LZwSJEBSCYHBzPv5DulnSVZJ2SfpvKU8rK6puSNonabWkCyRNws+w16bkT5b8IuSJutJmfqBJ7yAfbUxrknVKunbz5ihpv6R1kj4k6U24K+ZbU/Ii4OL+am9PyP0S/CDdLgVI+gqzcFPHByqKHc1vcjO+WzMD+FRF+hEbt4nc99MHWKA7JokXHASDmzzRVZ4tJ1O/M/vjQZKelXRx4VlnF5K3p+uitAU9GHg6Xec3yZPNCp9ukgcASb+RdDXwcIo6u1n+AnkC7M3KuEjeNViSJr+8m7BN0isV+fM3OTeNgwFD0gEaJqbXVDi26su4beV95r53Auc0bWzQZ0JQCILBzd/SdWZN+g245nvLmNnwHrLkf0td3NJfm67Tgat6qL8ty40+sCldF1VZC5jZbOB96XZDIb6n/ud/+tTqkUa2ojihxfx1bMYdBp2MCzjNjh0A7se/1WjglmYVH6Fvcje+q9MJrCil9WXc5vc5vM6zYrLyyJYetzaz3jGzjqIfjaB9QlAIgsFNNu86r+iS1swmmtkq3HHSa7Wlq/m0mW0xsw+YWdf2r5l1mtmNNLTsu5wupR/mr6fblWa2Kmnk57JjzGyBma0H1rXZnt6yDnfqMwT4rpmdk9piZvZu/Ax8KL6q3Vgot9jMnjCzS8ysSxs+TSifwJ1KQbXTqSp24YqYE8xsYW87I+mfNEwNV+I6FAeAB2vyvwpcn24vM3e9/ZacbmbHmdk8M1tNc8uQvrQ3j4nPmFlRYbDX4zb1K6c1O/5Zget1zAEeM3djPTQ9Z6iZzTSzG3DvnCe238Ogi4F02hAhwrEc6IULZ3z7tegK9xC+issuem/HV5YCri2VrXO4dG2hPgH/wD0eFuNuq2jLMGBNKd/fU9miy+Dvt/leenTh3KTs6fhZeX726/iqPN93c+GMe0As9mE/DTfMOTwEDCmVq3S4lNIeKJTdhytH/g54byFPpcOlUj3vKbejhXdwU+n9v576U3SXvLPN99qjC+eUbzwN199X9se4TeVXlsZnfp+XlfItSuk57wHccqboxlnAxKP99///HGJHIQgGMfJfw8W4ydyvcBMzgMdxb4PlLd9WWAssx1fZO/Ef8Q7gFfx8fqGkKyra8m9JH8XPhNfjpovDcRO3F/GV7wdxR0ADgqRd+Pb2zcAvaJxrPwvciHss3F0q9gjuQGod8BwuKIzBFR63AhcB56m78l0zLsEVIZ/Ht9SnpdCu6d9WDl9p1x07dCHpenxVvQYXjIak5+4BtuDfuk//IrvJs/fSOJb6rPm/Bu+PcXtdCs/hu0L5fR6miyG33jkNP3p5Bj+KGYsffTyRnj9D3a1mgjawJJUFQRAEQRB0I3YUgiAIgiCoJQSFIAiCIAhqCUEhCIIgCIJaQlAIgiAIgqCWEBSCIAiCIKglBIUgCIIgCGoJQSEIgiAIglpCUAiCIAiCoJYQFIIgCIIgqCUEhSAIgiAIaglBIQiCIAiCWv4HSDv6aUBbUxkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#####\n", "# 3. Make ROC curves to evaluate a model's overall useability.\n", "#####\n", "\n", "from sklearn.metrics import roc_curve, auc\n", "\n", "# a function to make 'pretty' ROC curves for this model\n", "def make_roc(name, clf, ytest, xtest, ax=None, labe=5, proba=True, skip=0):\n", " initial=False\n", " if not ax:\n", " ax=plt.gca()\n", " initial=True\n", " if proba:#for stuff like logistic regression\n", " fpr, tpr, thresholds=roc_curve(ytest, clf.predict_proba(xtest)[:,1])\n", " else:#for stuff like SVM\n", " fpr, tpr, thresholds=roc_curve(ytest, clf.decision_function(xtest))\n", " roc_auc = auc(fpr, tpr)\n", " if skip:\n", " l=fpr.shape[0]\n", " ax.plot(fpr[0:l:skip], tpr[0:l:skip], '.-', alpha=0.3, label='ROC curve for %s (area = %0.2f)' % (name, roc_auc))\n", " else:\n", " ax.plot(fpr, tpr, '.-', alpha=0.3, label='ROC curve for %s (area = %0.2f)' % (name, roc_auc))\n", " label_kwargs = {}\n", " label_kwargs['bbox'] = dict(\n", " boxstyle='round,pad=0.3', alpha=0.2,\n", " )\n", " if labe!=None:\n", " for k in range(0, fpr.shape[0],labe):\n", " #from https://gist.github.com/podshumok/c1d1c9394335d86255b8\n", " threshold = str(np.round(thresholds[k], 2))\n", " ax.annotate(threshold, (fpr[k], tpr[k]), **label_kwargs)\n", " if initial:\n", " ax.plot([0, 1], [0, 1], 'k--')\n", " ax.set_xlim([0.0, 1.0])\n", " ax.set_ylim([0.0, 1.05])\n", " ax.set_xlabel('False Positive Rate')\n", " ax.set_ylabel('True Positive Rate')\n", " ax.set_title('ROC')\n", " ax.legend(loc=\"lower right\")\n", " return ax\n", "\n", "\n", "sns.set_context(\"poster\")\n", "make_roc(\"Logistic\", logit, y_test, X_test, ax=None, labe=20, proba=True, skip=1);\n", " " ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Question**\n", "4. Use the ROC curves to select a threshold to balance the two types of errors.\n", "\n", "**Answer** It looks like based on the logistic model, a threshold of somwhere around 0.04 will give us a very high true positive rate (80\\% or so) before taking on \"too high\" of a false positive rate (just a tad over 50\\%). \n" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "## Part 4: Imputation\n", "\n", "In this part we ask that you explore the effects of imputation:\n", "1. Plot the histogram of `income`.\n", "2. Create a new variable `income_imp` that imputes the median or mean income for all the missing values and plot the histogram for this new variable.\n", "3. Compare the histograms above.\n", "\n", "\n", "4. Update your `poorhealth` prediction model(s) by incorporating `income_imp`. \n", "5. Compare the accuracy of this new model.\n", "\n", "\n", "And if there is time:\n", " \n", "6. Create a new variable `income_imp2` that imputes the value via a model.\n", "7. Update your `poorhealth` prediction model(s) by incorporating `income_imp2`. \n", "8. Compare the accuracy of this newest model." ] }, { "cell_type": "code", "execution_count": 96, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAELCAYAAADtIjDCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYNElEQVR4nO3df7BkZX3n8ffHweHHjuLq1sjGUX4IcUtiaggMgVjLjI6WP5aYsDgbxWyc/Niq5YdKReVHkfWPoHEgmmiiYO1myWRLqcoyKq5FaSnojBgJmUHHKEoJwgiyQrIiIwMoBr77xzmXadq+d7p77j3dd+77VXXr3D7nefo8/fBwP3O6n/N0qgpJkhba0ybdAEnS0mDgSJI6YeBIkjph4EiSOmHgSJI6cdCkGzBpSb4GHA3sAe6YcHMkabE4FlgB3FVVJwxTIUt9WnSSB4HDJ90OSVqkdlfVs4YpuOSvcGiubA4//PDDWb169aTbIkmLws6dO9m9ezc0f0OHYuA0b6M9b/Xq1WzdunXSbZGkRWHdunVs27YNRvgowkkDkqROGDiSpE4YOJKkThg4kqROGDiSpE4YOJKkThg4kqROeB/OYvbF907mvC+7eDLnlbSoeYUjSeqEgSNJ6oSBI0nqxMiBk+TQJBck2Z7kwSSPJLkryTVJXjqg/NOSnJtkR5I9SXYnuTHJG4c411lt2d1t3R3tcxmUkrTIjDRpIMnRwOdovgfhn4BtwE+Bo4DfAL4O/F1P+WXAJ4DXAT9u6x4MrAeuTnJqVb11lnN9GDgH+AlwA/Cztt6HgPVJNlTV46O0X5I0OUMHTpJ/BXweeCFwKXBpVf2s5/hzgOf0VTufJmy+Bby8qu5vyx4H3Ai8JckNVfWpvnOdSRM29wGnVdXt7f7nAl8EzgDOAz44/EuVJE3SKG9N/RFN2PyvqnpXb9gAVNUPq+o7M4/bq5sL2odnz4RNW/Z24ML24SUDzjUz7/bCmbBp690PnN0+vMi31iRp8RjqD3aS5cB/aR9uGvK5TwVWAt+vqi8NOH4Nzdtka5I8r+dcq4ATgcfaMk9RVduAe4EjgFOGbIskacKGfUvtRJq3y+6pqm8n+TXg9HbffcBnq+qmvjoz33G9fdATVtUjSW4FVrc/9/bVu7WqHp2lPduB57VlvzLka5AkTdCwgfOSdnt7ks3Am/uOvyvJx4H/3BMSR7fb783xvHfThM3RPfuGrddbVpI05YYNnGe329OAZcD7gI8AP2z3XQGcSTMT7ffasiva7cNzPO/Md2E/o2ffuPWelGQjsHGO+r1WD1lOkrQfhg2cmc96DgL+qqre2XPs/yT5v8A/AG9O8u6quhNIe7xGbNO49XodBazdj/qSpHk2bOA81PP7/+g/WFU7ktwCnASsA+7sqbOiv3yPmWO9zz9uvV67aO4RGsZq4PAhy0qSxjRs4Ozq+f2uWcrcRRM4R/TVOXKO533+gOcft96TqmozsHmO+k9KshWvhiRpwQ17H8tXe37vv7lzxr9ptzOfr8zUWTOocJLDgF9qH36t59DM78cnOXSWc63pKytJmnJDBU5V3Qvc3D5c3388yb8GfqV9uKPd3kSz/M2qJKcNeNoNwNOB7e3zz5zrHpqwWt6W6T/XWmAVzXTs/qnYkqQpNcqd+u9pt+9K8uTMriSHAFfSfA5yC20ItOuc/Wlb7MokK3vqHMfeG0hnnrfXzDeLXZbk2J56K2lmxAFsqqonRmi/JGmChl5Lrao+neR9wDuAm5PcTDMt+mTgF2hu3HxjVfXOLvtzmmnTv05zD88NNFc1rwAOAf6yfx219lxbklxJs4zNN5Jcz97FO58JXEuziKckaZEYaS2ydjr0f6RZEfolwGuBR4A/A07oXfesLf848JvAW4A7gFfRfEB/C/Cm2VaKbuueA7yJ5u21tW3dO2gW7TzTlaIlaXEZ6esJAKrqk8AnRyj/BM3VyMhXJFV1NXD1qPUkSdPH1ZYlSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdGDtwkvxJkmp/3jFHubOS3Jhkd5I9SXYkOTfJnOdO8uokn0vyQJJHknwzySVJDh63zZKkyRkrcJKsAS4Aah/lPgx8DDgJuBH4PPCLwIeALUmWzVLvAuAzwMuBrwLXASuBdwNbkxw2TrslSZMzcuC0VxibgfuBT81R7kzgHOA+4Jer6vSqOgM4Dvg2cAZw3oB6JwGbgEeAl1bVK6pqA3AM8CXgFOA9o7ZbkjRZ41zh/DHwYuC/ArvnKHdxu72wqm6f2VlV9wNntw8vGvDW2kVAgMuq6uaeenuA3wWeAM5J8qwx2i5JmpCRAifJrwJvB66uqk/PUW4VcCLwGHBN//Gq2gbcCxxBc8UyU2858Jr24ccG1LsTuAlYDrx2lLZLkiZr6MBJcgjwN8ADwNv2UfyEdntrVT06S5ntfWUBXgQcBjxQVd8doZ4kacodNELZ99AEwhuq6v/to+zR7fZ7c5S5u69s7+93M7tB9SRJU26owEnya8D5wLVV9bdDVFnRbh+eo8yedvuMeaj3FEk2AhvneI5eq4csJ0naD/sMnCSHAn8N/Jhm1tkw0m7nnDY9j/X6HQWs3c/nkCTNo2GucP6E5t6Z36uqHwz5vA+12xVzlJk59lDPvnHr9dsFbJvjeK/VwOFDlpUkjWmYwDmDZirym5O8ue/Yv2u3Zyc5Hbijqv6A5g8+wJFzPO/z2+2unn0zv79gxHpPUVWbae4V2qckW/FqSJIW3LCTBp7G3H+Uj2l/Zu6N+Vq7PT7JobPMVFvTVxbgNuBR4NlJXjjLTLWTB9STJE25fU6LrqqjqiqDfmimSQO8s923uq1zD82SNMuBDf3PmWQtsIpmFYKbes71GM2SNgBvGlDvGOBUmvt7rhvhdUqSJmwhV4t+b7u9LMmxMzuTrASuaB9uqqon+uptopk0cGGSk3vqrQCuomnzFVX14IK1XJI07xYscKpqC3AlzWoC30jy6SSfAG6nWRrnWppFPPvrbadZ3uYw4CvtitH/G/guzdt6NwOXLFS7JUkLY5QbP0dWVeck+TJwLk1YLKP5nOYq4MoBVzcz9S5P8o80y+isAQ4B7gT+AnhfVf10IdstSZp/+xU4VbWRfdxgWVVXA1eP8dyfBT47VsMkSVPHb/yUJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdcLAkSR1wsCRJHXCwJEkdeKgSTdgsTvqousmdu5dr5rYqSVpZAbOIvaBG74zkfOe/7KJnFbSIudbapKkThg4kqROGDiSpE4YOJKkThg4kqROGDiSpE4MFThJnp5kfZL3J/n7JD9I8liSe5NsSbJuH/XPSnJjkt1J9iTZkeTcJHOeP8mrk3wuyQNJHknyzSSXJDl4hNcoSZoCw17hrAWuB/4QOBK4Bfgk8ABwJvDFJH88qGKSDwMfA04CbgQ+D/wi8CFgS5Jls9S7APgM8HLgq8B1wErg3cDWJIcN2XZJ0hQYNnCeAD4OnFZV/7aqTq+q36qqlwBvAB4H/luSp9wSmORM4BzgPuCX23pnAMcB3wbOAM7rP1mSk4BNwCPAS6vqFVW1ATgG+BJwCvCe0V+uJGlShgqcqvpCVb2+qm4ccOxvgc3tw9/uO3xxu72wqm7vqXM/cHb78KIBb61dBAS4rKpu7qm3B/hdmgA8J8mzhmm/JGny5mvSwNfa7aqZHUlWAScCjwHX9Feoqm3AvcARNFcsM/WWA69pH35sQL07gZuA5cBr56f5kqSFNl+Bc1y7/UHPvhPa7a1V9egs9bb3lQV4EXAY8EBVfXeEepKkKbbfi3cmOQLY2D78eM+ho9vt9+aofndf2d7f72Z2g+r1tmljT5v2ZfWQ5SRJ+2G/AifJQcBHgcOBG6rq0z2HV7Tbh+d4ij3t9hnzUK/XUTQz6yRJU2J/r3A+AqwH7uHnJwyk3daIzzluvV67gG1Dll1NE5iSpAU0duAk+SDw+zRTntdX1X19RR5qtyuY3cyxh3r2jVvvSVW1mb0z5+aUZCteDUnSghtr0kCS9wNvBf6ZJmxuH1BsV7s9co6nen5f2d7fXzBiPUnSFBs5cJJcTrPiwA+BV1bVt2YpOjNV+vgkh85SZk1fWYDbgEeBZyd54Sz1Th5QT5I0xUYKnCSbgHcCP6IJm6/PVraq7qFZkmY5sGHAc62luW/nPpr7ambqPUazpA3AmwbUOwY4leb+nutGab8kaXKGDpwklwIXAg/ShM0wVxfvbbeXJTm257lWAle0DzdV1RN99TbRTBq4MMnJPfVWAFe17b6iqh4ctv2SpMkaatJAktcBf9Q+vAN4S5JBRW+rqk0zD6pqS5IraZax+UaS64Gf0cxseyZwLc0ink9RVduTXARcBnwlyRdogm4tzQKeNwOXDPUKF9j5B22ZdBMkaVEYdpbas3t+P6n9GWQbzdXJk6rqnCRfBs6lCYxlNJ/TXAVcOeDqZqbe5Un+EXg7zWc9hwB3An8BvK+qfjpk2yVJU2CowBllmvEs9a8Grh6j3meBz457XknS9PAbPyVJndjvtdS09Bx10WQmB+7a9B8mcl5J88MrHElSJwwcSVInDBxJUicMHElSJwwcSVInDBxJUicMHElSJwwcSVInDBxJUicMHElSJwwcSVInDBxJUicMHElSJwwcSVInDBxJUif8PhyN7PyDtkzozH4fjrSYeYUjSeqEgSNJ6oSBI0nqhIEjSeqEgSNJ6oSBI0nqhIEjSeqE9+FIQzjqousmct5dm7z3SAcOr3AkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnTBwJEmdMHAkSZ0wcCRJnfDGTy0ak7r5UtL88ApHktQJA0eS1AnfUpOGcP5BWyZ0ZtdS04HDKxxJUie8wtGiMbmrDEnzwSscSVInDBxJUicMHElSJwwcSVInDBxJUicMHElSJ5wWLU2xSa4ft2uTN512Zan8dzZwJA00qT+CBt2Ba+rfUktyVpIbk+xOsifJjiTnJpn6tkuS9prqK5wkHwbOAX4C3AD8DFgPfAhYn2RDVT0+wSZKC2qSqyt84F9eP7FzLzWTXUXDt9RIciZN2NwHnFZVt7f7nwt8ETgDOA/44MQaKR3AXLBU822a35a6uN1eOBM2AFV1P3B2+/Ai31qTpMVhKq9wkqwCTgQeA67pP15V25LcCzwPOAX4SrctlLRQlsqMraVoKgMHOKHd3lpVj85SZjtN4JyAgSMdMJbK5xlL0bQGztHt9ntzlLm7r+yTkmwENg55rlMBdu7cybp164asstf377pt5DqSptO1X143kfNO8u/IuK95586dM78eO2ydaQ2cFe324TnK7Gm3zxhw7Chg7Sgn3L17N9u2bRuliqQDzHfvvn/STejcPLzmFfsu0pjWwEm7rTHr7wKGTY8TgWXAA8AdI5xjNXA4sBvYuY+yathno7G/RmN/jW5/+uxYmrC5a9gK0xo4D7XbuZJz5thD/QeqajOweX6b9FRJttJcRe2sqnULea4DhX02GvtrNPbX6Lrus2mdUryr3R45R5nn95WVJE2xaQ2cr7Xb45McOkuZNX1lJUlTbCoDp6ruAb4KLAc29B9PshZYRbMKwU3dtk6SNI6pDJzWe9vtZUmenHaXZCVwRftwU1U90XnLJEkjm9ZJA1TVliRX0ixj840k17N38c5nAtfSLOIpSVoEpjZwAKrqnCRfBs6lmUmxDLgNuAq40qsbSVo8pjpwAKrqauDqSbdDkrR/pvkzHEnSAcTAkSR1YurfUptim4GteOPpKDZjn41iM/bXKDZjf41qMx32WarGXa5MkqTh+ZaaJKkTBo4kqRMGzhiSnJXkxiS7k+xJsiPJuUmWXH8m2Zyk5vgZ+M1SSZ7W9tmOtg93t336xq5fw3xL8qIkb0vy0SS3JXmi7YvXD1F3rLGV5NVJPpfkgSSPJPlmkkuSHDx/r2zhjNNn4469tu6iHX9Jnp5kfZL3J/n7JD9I8liSe5NsSbJuH/UnNsacNDCiJB8GzgF+AtzA3tUPPgSsT7Khqh6fYBMn5e8Y/H1CP+jfkWQZ8AngdcCPgc8BB9P049VJTq2qty5gWxfa2cDbRq007thKcgFwGfA4zQfAP6K5UfrdwOlJ1lfVI+O9lM6M1WetocceHBDjby3w+fb3+4BbaL6s8sXAmcCZSS6tqnf1V5z4GKsqf4b8af9jFs1APq5n/3OBb7XH3jbpdnbcJ5vb171xhDpvb+vcCjy3Z/9x7f9ABfzGpF/bfvTJHwCXA/8JeGH7P2gBr5/vsQWcBDzR/sH51Z79K2i+hLCAP590nyxQn4089g6E8Qe8HNgC/PsBx34L+Jf2Nbxs2sbYxDtvMf0AO9rO/Z0Bx9b2/Md82qTb2mGfjPQ/Pc3yRPe3dU4bcPzN7bF/mPRrm8c+GuaP51hjq/3DU8C7BtQ7huZfpD8FnjXpfliAPhs5cJbC+AP+qn0N/3PaxtiS+8xhXElW0Xwd9WPANf3Hq2obcC9wBHBKt61bVE4FVgLfr6ovDTh+Dc1l/pokz+u0ZRMy7thKshx4TfvwYwPq3Unz9R3LgdfOe8MXp6Uw/ma+I2zVzI5pGWMGzvBOaLe3VtWjs5TZ3ld2KXlZkj9L8t+TXJrkVbN8CDnTN9sHHKOa94FvbR+uXoiGTqFxx9aLgMOAB6rquyPUO9AMO/ZgaYy/49pt72dYUzHGnDQwvKPb7ffmKHN3X9ml5HcG7PtWkjdU1Td69g3bj6tZOv047tg6uu/YsPUONMOOPTjAx1+SI4CN7cOP9xyaijHmFc7wVrTbh+cos6fdPmOB2zJNdgJvBY6n6aNfAE4Hvk4za+b6vrcm7MefN26fLPW+HHXswQHcZ0kOAj4KHA7cUFWf7jk8FWPMK5zhpd26FlCPqvpA366HgeuSfJ5mBsspwMXAee1x+/HnjdsnS7ovxxh7cGD32UdopjjfA/x237GpGGNe4QzvoXa7Yo4yM8cemqPMklBVj7H3a8J7P0y0H3/euH1iXw4wx9iDA7TPknwQ+H2aad3rq+q+viJTMcYMnOHtardHzlHm+X1ll7qZO71739bY1W7tx712tdtR+2Tm9xeMWG8pGDT24AAcf0neT/PW4j/ThM3tA4rtarcTHWMGzvBmphoen+TQWcqs6Su71D2n3e7p2ffVdruGAZIcBvxS+3Cp9OO4Y+s24FHg2UleOEu9kwfUWwoGjT04wMZfksuBPwR+CLyyqr41S9GpGGMGzpCq6h6awboc2NB/PMlamnnv99HMS1dz1zg8dQrqTcA/AauSnDagzgbg6cD2qrp3gds3FcYdW+1bR59pH75pQL1jaO47eQy4bt4bPt0GjT04gMZfkk3AO2mWmXllVX19trJTM8YmfVfsYvoBXs/eu3GP7dm/kmbu/pJa2oZm6ujpwLK+/QfR/Kvr8bZPXtV3/B3sXVpkZc/+49q+neqlRcbop63s+675scYWzb9KZ5YdObln/4qe80790jaj9tm4Y+9AGX/ApW07fwScOGSdiY8xv4BtREmuoFlo8CfA9exd/O6ZwLU0/4MsicU7k/wm8EngAeA7wPdppka+hGaK6hPAxVV1eV+9ZW29X6dZPPEGmn9VvgI4BPjLmu7FE+eU5FeAK3p2vZimX26n6SsAquqUvnpjja2+hRW/ADxIs1TJSuBm4OU15Yt3jtpn4469tu6iHn9JXgd8qn24g703qva7rao29dWd7BibdFIvxh/gLJoVan9Mk/q3AOeyhNZQa/vhaOADwFdolsX4Cc37vbcDVzHHv7xo3s49j70r3f4Y+DJw1qRf1zz0yzqaf/XN+TOfYwt4Nc0Kwj9q/xvcClwCHDzp/liIPtufsbfYxx/NjZ377Ctg67SNMa9wJEmdcNKAJKkTBo4kqRMGjiSpEwaOJKkTBo4kqRMGjiSpEwaOJKkTBo4kqRMGjiSpEwaOJKkT/x89Dr8uVltHNQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#####\n", "# 1. Plot the histogram of `income`.\n", "# 2. Create a new variable `income_imp` that imputes the median or \n", "# mean income for all the missing values and plot the histogram for this new variable.\n", "#####\n", "\n", "\n", "\n", "# First create 'income_imp', and then add it into the predictor set\n", "income_imp = gssdata['income'].fillna(gssdata['income'].median())\n", "X['income_imp'] = income_imp\n", "\n", "# plot the original and the version with imputations\n", "plt.hist(gssdata['income'])\n", "plt.hist(X,alpha=0.5);\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Question:**\n", "3. Compare the histograms above.\n", "\n", "**Solution:** There is now a spike at the median compared to what was there before. They distributions are not all that similar in shape or spread (but center is very similar)." ] }, { "cell_type": "code", "execution_count": 173, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[316 130]\n", " [ 14 11]]\n", "[[316 130]\n", " [ 14 11]]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ageeducfemalemarital_never marriedmarital_no longer marriedrace_otherrace_whitesexornt_gaysexornt_heteropartyid_otherpartyid_repincome_imp
6942.014.010001011037.50
108232.015.011001000021.25
114322.014.001001000037.50
87124.012.000001010167.50
38040.012.01010000003.50
\n", "
" ], "text/plain": [ " age educ female marital_never married marital_no longer married \\\n", "69 42.0 14.0 1 0 0 \n", "1082 32.0 15.0 1 1 0 \n", "1143 22.0 14.0 0 1 0 \n", "871 24.0 12.0 0 0 0 \n", "380 40.0 12.0 1 0 1 \n", "\n", " race_other race_white sexornt_gay sexornt_hetero partyid_other \\\n", "69 0 1 0 1 1 \n", "1082 0 1 0 0 0 \n", "1143 0 1 0 0 0 \n", "871 0 1 0 1 0 \n", "380 0 0 0 0 0 \n", "\n", " partyid_rep income_imp \n", "69 0 37.50 \n", "1082 0 21.25 \n", "1143 0 37.50 \n", "871 1 67.50 \n", "380 0 3.50 " ] }, "execution_count": 173, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#####\n", "# 4. Update your `poorhealth` prediction model(s) by incorporating `income_imp`. \n", "# 5. Calculate and compare the accuracy of this new model.\n", "#####\n", "\n", "# re-use indices for splitting since now we added the imputed 'income' variable\n", "# Note: the response is unaffected so does not need to be redefined.\n", "X_train = X.loc[itrain]\n", "X_test = X.loc[itest]\n", "\n", "logit_imp1 = sk.linear_model.LogisticRegression(C=100000)\n", "logit_imp1.fit(X_train,y_train)\n", "\n", "knn15_imp1 = KNeighborsClassifier(15)\n", "knn15_imp1.fit(X_train,y_train)\n", "\n", "print(confusion_matrix(y_test,t_repredict(logit,0.07,X_test)))\n", "print(confusion_matrix(y_test,t_repredict(logit_imp1,0.07,X_test)))\n", "\n", "X_train.head()" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Question:**\n", "5. Compare the accuracies.\n", "\n", "**Answer:** Nothing has improved. The accuracies are identical for both the logistic and knn models (looking at varios different thresholds." ] }, { "cell_type": "code", "execution_count": 170, "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[330 116]\n", " [ 15 10]]\n", "[[322 124]\n", " [ 15 10]]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] } ], "source": [ "#####\n", "# And if there is time:\n", "# 6. Create a new variable `income_imp2` that imputes the value via a model.\n", "# 7. Update your `poorhealth` prediction model(s) by incorporating `income_imp2`. \n", "# 8. Calculate and compare the accuracy of this newest model.\n", "#####\n", "\n", "# the y_imp is observed incomes in train for building an imputation model to \n", "# impute values on training and test, and the X_imp is everything else in train \n", "# (we could use the true response health or poorhealth, but we decided not to here).\n", "\n", "income_train = gssdata['income'][itrain]\n", "\n", "# all the missing observations' predictors\n", "miss_index = gssdata['income'][gssdata['income'].isna()].index\n", "X_miss = X.loc[gssdata['income'].isna(),:]\n", "X_miss = X_miss.drop('income_imp',axis=1)\n", "\n", "# all the available observed incomes (and predictors) within train to be used to \n", "# build a model to predict income\n", "y_imp = income_train.dropna()\n", "X_imp = X_train.drop('income_imp',axis=1).loc[itrain]\n", "X_imp = X_imp.loc[income_train.isna()==False,:]\n", "\n", "# fit the model\n", "lm = sk.linear_model.LinearRegression()\n", "lm.fit(X_imp,y_imp)\n", "\n", "# do the predictions without noise, and turn it into a series for imputation\n", "y_miss = lm.predict(X_miss)\n", "y_miss_series = pd.Series(data = y_miss, index = miss_index)\n", "\n", "# create the imputed income variable for all observations\n", "income_imp2 = gssdata['income'].fillna(y_miss_series)\n", "\n", "# add income_imp into the train and test sets properly \n", "X_train['income_imp'] = income_imp2[itrain]\n", "X_test['income_imp'] = income_imp2[itest]\n", "\n", "# go back to the primary classifciatoin modeling...\n", "logit_imp2 = sk.linear_model.LogisticRegression(C=100000)\n", "logit_imp2.fit(X_train,y_train)\n", "\n", "knn15_imp2 = KNeighborsClassifier(15)\n", "knn15_imp2.fit(X_train,y_train)\n", "\n", "# quick peak at predictions to see if anything has changed \n", "print(confusion_matrix(y_test,t_repredict(logit,0.07,X_test)))\n", "print(confusion_matrix(y_test,t_repredict(logit_imp2,0.07,X_test)))\n", "\n" ] }, { "cell_type": "code", "execution_count": 155, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ageeducfemalemarital_never marriedmarital_no longer marriedrace_otherrace_whitesexornt_gaysexornt_heteropartyid_otherpartyid_rep
108232.015.0110010000
87124.012.0000010101
38040.012.0101000000
64843.016.0100010000
126461.010.0101010110
\n", "
" ], "text/plain": [ " age educ female marital_never married marital_no longer married \\\n", "1082 32.0 15.0 1 1 0 \n", "871 24.0 12.0 0 0 0 \n", "380 40.0 12.0 1 0 1 \n", "648 43.0 16.0 1 0 0 \n", "1264 61.0 10.0 1 0 1 \n", "\n", " race_other race_white sexornt_gay sexornt_hetero partyid_other \\\n", "1082 0 1 0 0 0 \n", "871 0 1 0 1 0 \n", "380 0 0 0 0 0 \n", "648 0 1 0 0 0 \n", "1264 0 1 0 1 1 \n", "\n", " partyid_rep \n", "1082 0 \n", "871 1 \n", "380 0 \n", "648 0 \n", "1264 0 " ] }, "execution_count": 155, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X_imp.head()" ] }, { "cell_type": "code", "execution_count": 171, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[319 127]\n", " [ 15 10]]\n", "[[324 122]\n", " [ 14 11]]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "//anaconda3/lib/python3.7/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] } ], "source": [ "######\n", "# Now to do imputation with uncertainty:\n", "# we can use the same imputation model from the previous part \n", "######\n", "\n", "# first add the standard y-hats just like before\n", "y_miss = lm.predict(X_miss)\n", "\n", "# we need to estimate the residual variance (MSE), sigma2_hat, from the observed incomes \n", "# that were used to train the model\n", "y_hat = lm.predict(X_imp)\n", "sigma2_hat = sk.metrics.mean_squared_error(y_imp,y_hat)\n", "\n", "# sample a residual from the assumed normal distribution\n", "e_miss = np.random.normal(loc=0,scale=np.sqrt(sigma2_hat),size=y_miss.shape[0])\n", "\n", "# create the income measurement with uncertainty to be imputed\n", "y_miss_series = pd.Series(data = y_miss+e_miss, index = miss_index)\n", "\n", "# imputed them properly into where they belong.\n", "income_imp3 = gssdata['income'].fillna(y_miss_series)\n", "\n", "# add income_imp into the train and test sets properly \n", "X_train['income_imp'] = income_imp3[itrain]\n", "X_test['income_imp'] = income_imp3[itest]\n", "\n", "# go back to the primary classifciatoin modeling...\n", "logit_imp3 = sk.linear_model.LogisticRegression(C=100000)\n", "logit_imp3.fit(X_train,y_train)\n", "\n", "knn15_imp3 = KNeighborsClassifier(15)\n", "knn15_imp3.fit(X_train,y_train)\n", "\n", "print(confusion_matrix(y_test,t_repredict(logit,0.07,X_test)))\n", "print(confusion_matrix(y_test,t_repredict(logit_imp3,0.07,X_test)))\n" ] }, { "cell_type": "code", "execution_count": 169, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AUC for logistic model when income was dropped: \n", " 0.6481614349775785\n", "AUC for logistic model when imputations were done with median imputation: \n", " 0.6481614349775785\n", "AUC for logistic model when imputation were done via linear regression: \n", " 0.6476233183856502\n", "AUC for logistic model when imputation were done via linear regression with uncertainty: \n", " 0.6590134529147982\n" ] } ], "source": [ "# Let's use AUC for evaluations \n", "\n", "fpr, tpr, thresholds=roc_curve(y_test, logit.predict_proba(X_test)[:,1])\n", "print(\"AUC for logistic model when income was dropped: \\n\",\n", " auc(fpr, tpr))\n", "\n", "fpr, tpr, thresholds=roc_curve(y_test, logit_imp1.predict_proba(X_test)[:,1])\n", "print(\"AUC for logistic model when imputations were done with median imputation: \\n\",\n", " auc(fpr, tpr))\n", "\n", "fpr, tpr, thresholds=roc_curve(y_test, logit_imp2.predict_proba(X_test)[:,1])\n", "print(\"AUC for logistic model when imputation were done via linear regression: \\n\",\n", " auc(fpr, tpr))\n", "fpr, tpr, thresholds=roc_curve(y_test, logit_imp3.predict_proba(X_test)[:,1])\n", "print(\"AUC for logistic model when imputation were done via linear regression with uncertainty: \\n\",\n", " auc(fpr, tpr))" ] }, { "cell_type": "markdown", "metadata": { "button": false, "new_sheet": false, "run_control": { "read_only": false } }, "source": [ "**Question:**\n", "8. Compare the accuracies.\n", "\n", "**Answer:** Things have improved! Using the uncertainty in the imputations has slightly improved the AUC for the classification model. But this is ONLY slightly!" ] } ], "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.9" } }, "nbformat": 4, "nbformat_minor": 1 }