teaching/fall-18/static/evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Applied Machine Learning (HOUSECS 59-03)\n",
    "# Fall 2018, Duke University\n",
    "# Evaluation of Machine Learning Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import datasets, linear_model, svm\n",
    "from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix, roc_curve, auc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin by revisiting the dataset that looks at housing prices in Boston: load the data into a Pandas dataframe, and plot it with respect to a single feature, LSTAT (% lower status of the population)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Boston House Prices dataset\n",
      "===========================\n",
      "\n",
      "Notes\n",
      "------\n",
      "Data Set Characteristics:  \n",
      "\n",
      "    :Number of Instances: 506 \n",
      "\n",
      "    :Number of Attributes: 13 numeric/categorical predictive\n",
      "    \n",
      "    :Median Value (attribute 14) is usually the target\n",
      "\n",
      "    :Attribute Information (in order):\n",
      "        - CRIM     per capita crime rate by town\n",
      "        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.\n",
      "        - INDUS    proportion of non-retail business acres per town\n",
      "        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)\n",
      "        - NOX      nitric oxides concentration (parts per 10 million)\n",
      "        - RM       average number of rooms per dwelling\n",
      "        - AGE      proportion of owner-occupied units built prior to 1940\n",
      "        - DIS      weighted distances to five Boston employment centres\n",
      "        - RAD      index of accessibility to radial highways\n",
      "        - TAX      full-value property-tax rate per $10,000\n",
      "        - PTRATIO  pupil-teacher ratio by town\n",
      "        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town\n",
      "        - LSTAT    % lower status of the population\n",
      "        - MEDV     Median value of owner-occupied homes in $1000's\n",
      "\n",
      "    :Missing Attribute Values: None\n",
      "\n",
      "    :Creator: Harrison, D. and Rubinfeld, D.L.\n",
      "\n",
      "This is a copy of UCI ML housing dataset.\n",
      "http://archive.ics.uci.edu/ml/datasets/Housing\n",
      "\n",
      "\n",
      "This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.\n",
      "\n",
      "The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic\n",
      "prices and the demand for clean air', J. Environ. Economics & Management,\n",
      "vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics\n",
      "...', Wiley, 1980.   N.B. Various transformations are used in the table on\n",
      "pages 244-261 of the latter.\n",
      "\n",
      "The Boston house-price data has been used in many machine learning papers that address regression\n",
      "problems.   \n",
      "     \n",
      "**References**\n",
      "\n",
      "   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.\n",
      "   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.\n",
      "   - many more! (see http://archive.ics.uci.edu/ml/datasets/Housing)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# load dataset\n",
    "boston = datasets.load_boston()\n",
    "print(boston[\"DESCR\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CRIM</th>\n",
       "      <th>ZN</th>\n",
       "      <th>INDUS</th>\n",
       "      <th>CHAS</th>\n",
       "      <th>NOX</th>\n",
       "      <th>RM</th>\n",
       "      <th>AGE</th>\n",
       "      <th>DIS</th>\n",
       "      <th>RAD</th>\n",
       "      <th>TAX</th>\n",
       "      <th>PTRATIO</th>\n",
       "      <th>B</th>\n",
       "      <th>LSTAT</th>\n",
       "      <th>MEDV</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00632</td>\n",
       "      <td>18.0</td>\n",
       "      <td>2.31</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.538</td>\n",
       "      <td>6.575</td>\n",
       "      <td>65.2</td>\n",
       "      <td>4.0900</td>\n",
       "      <td>1.0</td>\n",
       "      <td>296.0</td>\n",
       "      <td>15.3</td>\n",
       "      <td>396.90</td>\n",
       "      <td>4.98</td>\n",
       "      <td>24.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.02731</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.07</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.469</td>\n",
       "      <td>6.421</td>\n",
       "      <td>78.9</td>\n",
       "      <td>4.9671</td>\n",
       "      <td>2.0</td>\n",
       "      <td>242.0</td>\n",
       "      <td>17.8</td>\n",
       "      <td>396.90</td>\n",
       "      <td>9.14</td>\n",
       "      <td>21.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.02729</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.07</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.469</td>\n",
       "      <td>7.185</td>\n",
       "      <td>61.1</td>\n",
       "      <td>4.9671</td>\n",
       "      <td>2.0</td>\n",
       "      <td>242.0</td>\n",
       "      <td>17.8</td>\n",
       "      <td>392.83</td>\n",
       "      <td>4.03</td>\n",
       "      <td>34.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.03237</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.18</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.458</td>\n",
       "      <td>6.998</td>\n",
       "      <td>45.8</td>\n",
       "      <td>6.0622</td>\n",
       "      <td>3.0</td>\n",
       "      <td>222.0</td>\n",
       "      <td>18.7</td>\n",
       "      <td>394.63</td>\n",
       "      <td>2.94</td>\n",
       "      <td>33.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.06905</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.18</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.458</td>\n",
       "      <td>7.147</td>\n",
       "      <td>54.2</td>\n",
       "      <td>6.0622</td>\n",
       "      <td>3.0</td>\n",
       "      <td>222.0</td>\n",
       "      <td>18.7</td>\n",
       "      <td>396.90</td>\n",
       "      <td>5.33</td>\n",
       "      <td>36.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      CRIM    ZN  INDUS  CHAS    NOX     RM   AGE     DIS  RAD    TAX  \\\n",
       "0  0.00632  18.0   2.31   0.0  0.538  6.575  65.2  4.0900  1.0  296.0   \n",
       "1  0.02731   0.0   7.07   0.0  0.469  6.421  78.9  4.9671  2.0  242.0   \n",
       "2  0.02729   0.0   7.07   0.0  0.469  7.185  61.1  4.9671  2.0  242.0   \n",
       "3  0.03237   0.0   2.18   0.0  0.458  6.998  45.8  6.0622  3.0  222.0   \n",
       "4  0.06905   0.0   2.18   0.0  0.458  7.147  54.2  6.0622  3.0  222.0   \n",
       "\n",
       "   PTRATIO       B  LSTAT  MEDV  \n",
       "0     15.3  396.90   4.98  24.0  \n",
       "1     17.8  396.90   9.14  21.6  \n",
       "2     17.8  392.83   4.03  34.7  \n",
       "3     18.7  394.63   2.94  33.4  \n",
       "4     18.7  396.90   5.33  36.2  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.DataFrame(data=boston.data, columns=boston.feature_names)\n",
    "df[\"MEDV\"] = boston.target\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1a1e61bdd8>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnX+MJMd1379vZmdszq0E8novASFqe6XAiHMKDIV3cRzIEQSdE8iMETuBLMhYEQfIwOqGDnCBkcTnLBDbAfYPUkgs/UMxG1nyhTOxrDhx5D8EOARNIwYCyFnKlESBocnEuxtZhKi7I6GcVsidbyt/TPeod7arurq7+lfN9wM0dq+3u/p138y3Xr969UqUUiCEENJ9ek0bQAghxA0UdEII8QQKOiGEeAIFnRBCPIGCTgghnkBBJ4QQT6CgE0KIJ1DQCSHEEyjohBDiCSt1XmxtbU1tbGzUeUlCCOk8zz///A2l1Lms42oV9I2NDezt7dV5SUII6TwicmBzHEMuhBDiCRR0QgjxBAo6IYR4AgWdEEI8gYJOCCGeYCXoIrIvIl8TkRdEZC/ad1ZEnhGRV6KfD1Rh4HQ6xcbGBkQEKysrEBFsbGxgOp2eOO6xxx5Dv9+HiJzY+v0+3vWud83PXVlZwWOPPTZve21tbX7s2toaptPp/Jq9Xi/1WkXst21Ld7xLmwghnqKUytwA7ANYW9j3BIBr0e/XADye1c6FCxdUHiaTiRqNRgrAqW00GqnJZKKUUmo8HqceY9ouXbqkBoPBqf39fl8Nh0Pttcrab2pLd/x4PM7VDiHELwDsKQutFmWxBJ2I7AO4qJS6kdj3MoD3KaVeE5EHAfyhUuqvmtq5ePGiypOHvrGxgYMDffplGIbY39/HysoK7t27Z91uEeJr5UFnv64t3fH9fj/1/orYRAjpHiLyvFLqYuZxloL+ZwDewMw7/LdKqV0ReVMpdX/imDeUUqfCLiKyBWALANbX1y+YBHqRXq8Hk30iguPjY4iIdZtFia+VB539uray7teFTYSQ7mEr6LaDou9RSj0M4CcB/IKIvNfWEKXUrlLqolLq4rlzmTNXT7C+vm71936/n6vdImTZkuecvPt191fEJkKIv1gJulLqm9HP1wH8LoAfBfCtKNSC6Ofrro3b2dnBaDRK/dtoNMLOzg4AYGtrK3fbly5dwmAwOLW/3+9jOBxqr5WHNPtNbemO39raytUOIWRJyQqyAzgD4C2J3/87gA8A+DhODoo+kdVW3kFRpWYDhWEYzgcsAagwDE8NCI7HY9Xr9U4NcvZ6PXX+/Pn5uf1+X43H43nbQRDMjw2CQE0mk/k1ReTEtXT7bey3PcfltQkhfgDLQVEbQX8ngK9E29cBbEf7AwDPAngl+nk2q60igu6KsoKYN2OFEEJcYSvoVoOirsib5eKK6XSKra0tHB0dzfeJCK5cuYInn3zSqo28GSuEEOIK14OinWZ7e/uEmAOzN5OnnnrKeoLO4eFhrv2EEFI3SyHoOtFVSmF7e9uqjbyZKYQQUjdLIegm0bX1sPNmrBBCSN0shaDv7OxoJx/Zetibm5vY3d1FGIYQEYRhiN3dXWxubro0lRBCCrMUg6LArHjXU089dWIm5mg0oigTQloPB0UXePLJJ/H000/TwyaEeMvSCDohhPjOStMG1MViLvrBwcG8ZAC9dEKIDyyNh56Wi350dGSdtkgIIW3HK0E3rerDiUGEEN/xRtDjkMrBwQGUUvOQSizqnBhECPEdb2LoWSGV27dvnzqHE4MIIT7hTR66abWf0Wh0SuyDIMAnP/lJDogSQlrP0uWhm1b7WRRzAFhdXaWYE0K8whtB19Va0S0efXBwcGrglBBCuow3gq6rtRKGofacxYFTQgjpMt7E0HWkLW6xCBepIIS0maWLoetIeu46mItOCPEB7wQ9bXLR5uYm9vf3taLOXHRCiA94Jehpk4s++tGPYm1tDb1eD7dv38ZwODxxDnPRCSG+4JWgp00uunPnDm7evAml1PxnEAQsoUsI8Q5vZooCs6yVLO7evYvV1VXcuHGjBosIIaQ+vPLQ+/2+1XF1DYKaioURQohrvPLQdZOIFqljEJT11wkhdeOVh25KTYwZDAa1DIKy/johpG68EvSdnZ1TWSyLvPWtb63FQ2b9dUJI3Xgl6AC0FRdjbt26VYsdrL9OCKkbrwR9e3sbd+/eNR5Tl6DqioUx550QUhWdFfS0DJKscEadgqorFsYBUUJIZSilatsuXLigXDCZTNRoNFIA5ttoNFJBEJzYl9zCMFSTySS1rTAMlYhojyGEkCYBsKcsNLaTaYu6DJL77rvv1OpEo9FI6xkztZAQ4hOdDLnoQiu3bt3KFeZgaiEhxCc6WQ99Y2MjdZp/3rrmunVIRQTHx8dlTCSEEGd4XQ/dVQYJUwsJIT7RSUF3lUHC1EJCiE9Yh1xEpA9gD8CfK6V+SkTeAeBzAM4C+DKAR5VSd0xtNLEEXRbT6RTb29s4PDzE+vo6dnZ2OCBKCGkVVYRcrgJ4KfHvxwH8ulLqhwC8AeDn85nYDuLVjI6Pj7G/v1+JmLPqIiGkDqwEXUQeAvD3AXw6+rcAeD+A34kOuQ7gZ6owsOukraK0tbVFUSeEOMfWQ/8EgH8OIE79CAC8qZT6i+jf3wDwtrQTRWRLRPZEZO/b3/52KWO7CFMjCSF1kSnoIvJTAF5XSj2f3J1yaGowXim1q5S6qJS6eO7cuYJmdhdWXSSE1IWNh/4eAP9ARPYxGwR9P2Ye+/0iEs80fQjANyuxsOMwNZIQUheZgq6U+mWl1ENKqQ0AHwbwB0qpTQDPAfhgdNhlAF+ozErH1DlIydRIQkhdlMlD/yUAvygir2IWU/8NNyZVS92DlKy6SAipi05O/S+Dq7IBhBBSF15P/S+DbjDy4OCAueKEkE6zdIKuG4wUEeaKE0I6jXeCnjXgmTZIKSKnqi4yV5wQ0jW8EvS0Ac9HH30UIjIX97RBSt04wsHBAb10Qkhn8GpQVDfgGaNbvch03mAwwGc/+1lmpRBCGmMpB0WzZl/qwihpYZiYu3fv4urVq07sI4SQKvFK0G1mX6aJfhyG0XHz5s1SdhFCSB14I+jT6RS3b9/OPE4n+nWHVFhSlxDiGi8EPR4MzfKks6bcB0GQur/X61kLro1Q552tSvEnhFihlKptu3DhgqqCMAwVZtUeT2xBEKgwDJWIqDAM1WQyMbYzmUzUcDhMbWs0GlmdPxqNtOdNJhOtrQBUGIa52ySE+A+APWWhsV5kufR6vdTUQxHB8fFxyhl6ptMpLl++jHv37p36W1Z5AFNZgZ2dHWxtbZ2qjZ5lL0sVEEKWKsulTInaOJwhIlhZWcFHPvKRVDEHsrNoTLXP0xa6sLGX9dQJIbZ4IehFS9QmY9kAtEIeo+sg4k5B97azvr6eKcA6e1lPnRBijU1cxtVWVQxdqe/Hp23j5UrpY+9pmy5unRbjTjsvK3aus5cxdEIILGPonRH0IoKd1YatmCMaYE27pq1QlxFmF/dOCOkuXgm6Cy81y5Mu6qWLSOqxIpJqA4WZEJIXW0HvRJaLi0yPrDovtixec21tLTX/PQgC3Lhxo/T1CCHEqywXF5kerrJCmF1CCGkrnRB0F5keebNCRMSqnVu3bqUep9tfFM4WJYRk0QlBL5qWmNXGYDDAcDg8dWwQBLhy5YrVNetIK6x7YWtCSEexCbS72qrMcrEZcEw7xnSeTZvj8fjUwKjrtEJdJk1aqQBCiH/ApyyXLPJkwUwmExUEwal0RNsMlGROeVqGi4io8Xjs9P7yZNIQQvxjqQTd1oOdTCZqMBgUTlG0TX3M8pzzpi/SQydkuVkqQbf1YPNOJorFNhbffr9vdZ7Jcy6SU8/ZooQsN0sl6LYerE74szz1vOeYPOei3jYnJRGyvNgKeieyXLLQZcE88sgjJ1L9zp49m6vdfr+fWSFxEVP2zXQ61U5uyspv39zcxP7+Po6Pj7G/v89Fqwkhp1hp2gAXxOJ29erV+azNXq+HT3/607h79y4A4ODgAMPhEL1ez6pG+mg0yi3mvV7vxELUSdGNUw91sHoiIaQsXnjoMd/73vfmv9++fXsu5jF37tzBAw88oF1qLp5MFIYhdnd3EYah8Xr9fh/ALG99OBzOO4q0PHFTPfS8OfWEEJKGN4Jus4AEANy8eRO3bt1CGIYYj8cIwxAigjAM8fTTT0MpNa/VYlp0OgxDXL9+HWEY4ubNm7hz586Jvyc9dcAcUtnd3WUIhRBSmk4U57JBtwydidFohN3dXQCzDuHw8BDr6+t45JFHcP36daNHffnyZeMxwMkl5biUHCGkKF4V57KhSAz66OgIH/vYx/Doo4+emFb/qU99SivUcTjmi1/8Yq4l5VyULyCEEBPeCHqaYA6HQ22RrZjvfve71p69iMwzTPIuKbe5uTmPy8chnjpCLSzqRcgSYZPb6Gqrcgk6pdJztYvknuu2ZK540SXl6qTJCUnMmyfEHVimiUUm8s4O1W2LQtiF2ZtNlQzowrMhpEs4E3QAPwjgjwF8BcDXAfxatP8dAL4E4BUAvw1gmNVWE4LuYum5fr9vXb2xjJ2uPdqminqx9gwhbnEp6AJgNfp9EIn4jwH4PIAPR/ufAjDOaqsJQVeq3OLQRdYuzVOiNz6nCo+2KWFldUhC3FJJyAXACMCXAfwtADcArET7/zaA3886vylBX8RW1PN6ypPJRA2HwxNtrKysnKrwuCjWVQlvWkcRi22VcW166IS4xamgA+gDeAHAbQCPA1gD8Gri728H8GJWO3ULus4zthksLSI+yTrrNp1FTJUeral+u81bQJFQEGPohLilKg/9fgDPAfg7KYL+Nc05WwD2AOytr6/XdPvZomIS23iRirxClieUkxTrOjzaItcoI8zMciHEHZUI+qxd/AqAf9b2kEuWgGV500WELI+g9/v9udiNx+PKPdoibwEMnRDSDmwFPXNikYicE5H7o9/vA/ATAF6KPPUPRoddBvCFrLbqRDfx5/DwENPpFN/5zneM5y/OAl2szZKGruhXGvfu3YNSs5mp169fx+XLlyuddFRkMWvTMySEtJAsxQfwIwD+BMBXAbwI4F9G+9+JWTrjqwD+I4AfyGqrLR56mdx0U+aKzRJ3ulWPbLzeMmGMIuETeuiEtAMs+8Si8XisHQQsO3t0MBicymaJ2zSFckTEKvSh6yzKhmXydggc3CSkHSy1oOvS9cbjsVLK3exRk3AX9dBNqYZNeMsc3CSkeZZa0LNCBS5mj5YR9cUt6fXm7WzKpDZSrAnpBkst6HnCGlWKuqn9ZJZLUkjzhoOKeugMpxDSHWwF3ZvyuUlsMjriRZfDjGXm0hgMBhgOh8Zj+v2+MRvk+Pg4dcHnvHXdb9++XagkbtoKTzaZPL7B8sLEK2xU39XWZAxd530W8YiTWS7IGV5BhmddJBxUxLNmvRW+pZDugGUOuShlHx/OE3axqbqoG/TMIxpFwkF5Qy/LlpKY9nlYtmdAusvSC7otaV7acDjU5pNniXGWt57MtrHBVtDzetbL5J3q7tXVsySkaijoOZhMJifyx4MgUOPx2Oht6yYY2RTo0nn6adh66kW8ymXJctE9wzKTvAipEwp6DvJ6cCZPPm3SkY2nb6oMmWVLU551VzoE01vTsrylkG5DQc9BXg8uaztz5oyVZ23Ki08Ky2Qy0dqSx9t3SZdCNqZYeVc6JbLcUNBzkMeDs9100/WTWxyrtRmca1tWSpcGFLvU+RCSBgU9BzYeXF5BT3rfWbFaG7Fum4C2rYPJgp446TIU9BzYeHC6bJgs79umfRuxbpuX2bYOhhCfoaDnxGYh58VMmMV9JmEztW8r1m3yMtvWwRDiMxT0kiTFMwiCzDrnppBNWpt5/rZYCriOhZ5taFMHQ4jPUNAjXC1yXHQbDodzL77IIs3j8djYflvEnRBSHRR0VSwsYBrErGLLijkXLcNbFfTKCakfCrrKP3BXR530NA9bZ0vZ7JpF4S0rxl2M9RPiAxR0lT+1zmbafh0eepmOJa4Vs3jvaTNY83r0XczGIcQHKOgqn4c+mUxqF3MROTEbNPZqy4R8giDIFaYJgsDam+5ivjwhPkBBV/m8xSLhDRebzs4i22g0Kv2WYfKmuzijldQDw2zVQkGPsPmgNeGdxyKXVdXRtMU1Y5L3VmSxDZNALz6nrA6SHvrywTBb9VDQLWliIHRR1MueHwTBPNTiIkNn0ZvWrc4UT65aPGZxGw6HTr/ctt4gvcZ6YCdePRR0S5oKtTS1JRen1oVnkhUcTR1e7IVldYqDwcCZmObJtKHXWA8Ms1UPBd2SpgW2qS1exCNLrLM6vDAMrTpFV96arTdIr7E++Kyrh4JugauYcx3b6uqq1cIZebbBYFAqhg/MvDCbZ5iVb28bFjFdy2YsgV6je/g2VD0UdAu6EG6JB06zSgAU3WIxLXN+UQ+9iBDYXMuU7UOvsRo4XlEtFHQLuuKdVznhKf4CFjk3GUM3vT3kTRU1ia7tIHYQBPQaiTdQ0C3ogodeV4dhW00yDs+EYajG47FVjF0nokXDIklvUHfdeNIWvUbiAxR0C5pOWWzTNhwOrd5YYmG0eXamfHZTR5B1XlKgmwqtsLMgdUJBt0SXY+1i60pIJ95sywYMBgPrMNCi4GV1BLqwSNp5g8FArayspNpXpcByEJDUDQW9AEUrHOq28XjcSMGvMtuZM2cqadcmDdLk6eb5fwmCoNLPCdP0SN1Q0EtgEg9brzspKq47CtttMplUJtBFNlPcOytunudtp+rURKZEkrqxFfQeyCl2dnYwGo1O7Q+CAFeuXEEYhhARhGGIS5cuQUROHCci+NCHPjT/9+bmJvb39xGGYeW2J9ne3sbR0VGt1zRxeHiI9fX11L/p9tv+veixRSh6D4RUjo3qu9q64qErdXpR6HhLi5Wm1R9fPM40gFfl1sSbQa/XS90f15yxeaZp/x82MXTXdWNsbWEMnVQJXIVcALwdwHMAXgLwdQBXo/1nATwD4JXo5wNZbXVJ0JVyN828qWyauCZLE4Ozi/c7HA5TUyNj8bfJFFnMLBmPx6farHpAVGcLxZxUiUtBfxDAw9HvbwHwpwDOA3gCwLVo/zUAj2e11TVBt42VZh3XZL67UqqR6yafyZkzZ6zeTvIueN21wUl2AqQozgT91AnAFwD8XQAvA3hQfV/0X846t2uC7spDbyp9MauuSZs3mxBG3YOTZQSZYRpShkoEHcAGgEMAbwXw5sLf3sg6v2uC7qpUaxMeum21xDZvyTK+aZg6UtfecFlB7trbBGkXzgUdwCqA5wH8o+jfVoIOYAvAHoC99fX1eu7eIS4WU9DF0HWDhy6FsE3eeZFxBJNo6kQ2rSxwWW+4rCDbvE0wJEN0OBV0AAMAvw/gFxP7vA+5uCQtE8Z2un3ebTwez6/bJg/90qVLhe7X5HEn30LiOjO6csBlvGGT3S7i/QzJEBPOBB2AAPj3AD6xsP/jODko+kRWW8ss6LovdFpVwLIivzipqS31asrc1+K5yQFU00Idi+e4/v+zFd6iYTmGZIhSbgX9x6MP11cBvBBtjwAIADyLWdriswDOZrW1zIJueuVOS8crK+o6L7aO7fz5863sKMqIY1bHaNO2KaTC2afEhDNBd7kts6Dr0vZ0dUdcLGix6DnWIeo2JXWb2FyELyaTibFTKQM9dGKCgt4ydILe6/W0QuNiZmkQBHOvsI66Lm0J7yyKbXJcoQxVCS9j6MQEBb1lmMICui9umzJU0jbd4GOZNUrLiLbp76483SqFl1ku1dH1Z0tBbxlZYYg0wakzdBHH7m3FM2vCUp2dUbzYtemNxmUsuuviYIsv9+nD2w8FvWVkDaqlCU7aOcnsDpdiHl8vy7tOhi9cdzhFPfvhcJiZ6VJ1jXTf8EEEY3wYn6CgtxCTYOo+XCYvKUtQTQs3L35JbdIbF2PRLlMik3YU8e6zOgNXgt4mrzXLljK2+iCCMT5kEFHQW4pLzydLUNNy3OPl4xa/5LYDsItfaBcpkYti46KDMH15iwpdm7zWLFvK2uqDCMb40DlR0FuMSy8vK5XO5lqmNmy/0EVEVicwRToIWw+9jNC1SRiybClra5vutSxt6oiLQkHvGE2+HucR0LQ2i4RJTIW3TO2lvXXo6rcsbqYceZtn1SavNcuWsraWEcE2haXabFMeKOgdoqwHUdXrdZUetW4QOP7Sra6uaq9vU9dFd82szsr0hW+T11q1h65UMRH0wRtuIxT0DlHHly/t73ni36YvdNGCW4v2pa1yFMf7gyA4MTEqCIJC9uhCM1lLCOps9DWGXpQ2dXo+QUHvEFW/yqd9uQeDgVUWjM2Sbnk99Di2b9NG3JGkLV+nWz80yx7bYmi68FJbXt2rzHIpSpvCUj5BQe8QVXs1RUIieeyYTCZWnUO8pU3DNx1vysDJG9OPhS0pdKaOp8uYBL0qsaeHXg0U9A5R9etx2VmbNsJmm/aY5vEXzT2PtzTSKlbqaroUEaE2eeppmD5TVZcvaEtYyico6B2jSoGo2kNXKl+nsdheWfuCIEjNrb906dKpY13ExrsgWqZOqmovuu2dXRehoHtOni9NmgDZFtbKyiyJySPKix6/67ovcRpjng4qeX+6DiLrXtsUVjDFshnn7h4UdI8p4iEuCrKpPntaNkzW9dKOsR1sdFEmeHEzrddqEi6be+2CIDbpoRP3UNA9xsUXMs+ixbae7mKnYbtYcxWCbtqSdtt2dMlzuiCIdcbQGWKpHgq6x7jwELNEKe1LX+R6pvz3rCyTqjZTrrbNvXYhhq5UPVkuRcYfmhD/rnc6FHSPceEhZn0RbYS2iEeaR0Tj8I9LMV9dXZ3bkufNIOttxJVAdE148nwWm+oIu9IBm6Cge4yrD6hJPLIGKvOWJoivk6fmeZxm6KpEb3Ii0mRiX5DM9ZffVLqgSeEp0pnkeVtsKlSVFTas8vm66qAp6J5TtSfn6kuQxyPXiWnWakS2HnbSZtP96QaGqw5RNBmbL9qZ5LG5qcFkl85JHlx20BR0UgpXH0YXIZN+v59azyWPmC9i+pIvpiyaOiVTTZk8zyPuLKoUPFOnVLQzyfM5aaOHXqUNLu+Xgk5K48Irtckxr3JR6aL1XpLilPV2YFPvJsbUjqlkcFmyhNdkk03bNp+TNsXQq+o0k7jsoCnopBXohDPpdacV3nK1LYptMrSRVhqg6HVsva4ibbsQdJO3OJmYa9+4pOksl7L/f3mgh068w0U2TdktmS2jE/GydmTVd49nn7pqOy8mb1F372lVMbtOnW8JjKETLymTTWMreGXOjz0mnbClrZKU5XWVHQw2tV2EIrF7AKWv20bqfEtglgtZKsp6xnEWTNl2lMqeXanzsNNi6K7ePJKLfPiYXUPsoKCTTlDEk02GSbLqydhs/X5/3sZ4PJ4P0vb7/VPldheFXZflUvStYTAYnFilaXF8QXfvts86zVtsOv+dZENBbyldmwnogqx7nkwmxkyXxfizKVVw8Vq2cev4XFfCZuOhJzuO5H3a2u5SdIt8Lpfxs9wUFPQWsoyekO09Tyb6ZeZsi3zZXl8XWjDF0IuIXZ63BdslAXW2180yfpabhILeQpYxVpnnnnXhjLLPzZSqmBQh2zBJns5kMcslb/mDrK2pkr3L+FluEgp6C+lCHW3XuLhnl8+tyGxJ18LlckGPpgR0GT/LTWIr6D2Q2lhfX8+1v2tMp1NsbGyg1+thY2MD0+nUyT27fG6bm5vY39/H8fEx9vf3sbm5Of/bzs4ORqORVTsHBwdYW1s7ca+25LE7CAKEYQgAEJETfxuNRtjZ2bFuyyW+f5Y7i43qu9qW3UP3Oe6ou7cy8e+stquaDFJkUDVv9Umb+HraqlBtGYT0+bPcRsCQSztp05fSJaaYqot7zkonrIo8g5t5Bk+zMnv6/X6jGSw25/n6WW4jzgQdwGcAvA7gxcS+swCeAfBK9PMBm4tR0P0lT0w1rxA07Q1OJva10/PYaGrXZYXFIs+uimfODqA4LgX9vQAeXhD0JwBci36/BuBxm4tR0P3FNuuhiFC0IaMiz4BpHhtt1jCtwu6s9l0/86Y75a7jNOQCYGNB0F8G8GD0+4MAXrZph4LuL7Zf2CJC0YaMiqKzULNsdCF0RWrlZNnl+pm3oVPuMlUL+psLf3/Dph0Kut/YvFIXEYosMajrVX4xr9yFh57Wbp7aLUWrWdbtobehU+4yrRF0AFsA9gDsra+v13LzpL0UEYqsollNvcpnhWHqyOax6ezaEEOnh16OqgWdIRdSiDICk+aFNykUafdSpnhWVeGoqrJc8sAYejmqFvSP4+Sg6BM27VDQiVJuhaLKV3kbO5u+ly55vk2lnvqAM0EH8FsAXgNwF8A3APw8gADAs5ilLT4L4KzNxSjoxDUuBC1NlJvwKF2Ho3T3ZrrvqqCHXg6nHrqrjYJOXFNWKMbjcWrBLlfphHlEU5dJYyoXbLpGm8YemniT8CnvnYJOloYyMeK8hbLyhHKKiOZkkr4yUhGxNYloFQJbRfpkGVt8eiOgoBOSQVaWSlnBazpl0CSirgW2qvTJonRpbMEGCjohGZi887SFofN6eE1P6qnTQ68qfbIovuW9U9AJyUAnQiIyjzOXicHW6aHnHdh1LbCu0yebevZthYJOSAa6QcgzZ840mnud97ws4Xad5ZJ2nq2A2qaCdqnkch1Q0AmxYHHa/eK6pnVWGCxaAqBOb1QnlDZ176us96Ozteosl7oyaSjohOSkydf0Mh5lnfFi0zPKEjfb59uV+HedbwG2gi6zY+vh4sWLam9vr7brEZKHXq+HtO+DiOD4+LjSa29sbODg4ODU/jAMsb+/X9m5eSnzjGzPrfN+ylCnnSLyvFLqYtZxXFOUkIgm18k8PDzMtT9J2lqoVa03avOM0taWtT0XqPd+ylDm/6wybNx4VxtDLqTNtLFyo224p65Yrk2pAReZNV2Y5VlniA6MoROSn6aEpEtZGaZnZJOP3nahtqWNMXQKOiEtoUqxq0tIuzKg6Yq2ZblwUJQQz5lOp9ja2sLR0dF832g0wu7uLjY3N51eqysDml2Dg6KEEACYNPKVAAAFJklEQVTA9vb2CTEHgKOjI2xvbzu/VlcGNH2Fgk5IDegyP+qgzmyMzc1N7O7uIgxDiAjCMMTu7i4ANHb/S4VNXMbVxhg6WUaaHvBsuq5J0/fvA7CModNDJ6Ri6gx5pNF0GKTp+18mKOiEVEzTE1B0YRDXA6I6mr7/ZWKlaQMI8Z319fXUzI86ZqDGbG5u1ibgi7Th/pcFeuiEVEzTIY+mWfb7rxMKOiEV03TIo2mW/f7rhBOLCCGk5XBiESGELBkUdEII8QQKOiGEeAIFnRBCPIGCTgghnlBrlouIfBvAdwHcqO2ixVkD7XRJF+zsgo0A7XRNF+wMlVLnsg6qVdABQET2bNJvmoZ2uqULdnbBRoB2uqYrdtrAkAshhHgCBZ0QQjyhCUHfbeCaRaCdbumCnV2wEaCdrumKnZnUHkMnhBBSDQy5EEKIJ9Qm6CLyARF5WUReFZFrdV03LyKyLyJfE5EXRKQ1lcRE5DMi8rqIvJjYd1ZEnhGRV6KfDzRpY2RTmp2/KiJ/Hj3TF0TkkSZtjGx6u4g8JyIvicjXReRqtL9Vz9RgZ6ueqYj8oIj8sYh8JbLz16L97xCRL0XP87dFZNhSO39TRP4s8Tzf3aSdhbFZp67sBqAP4H8BeCeAIYCvADhfx7UL2LoPYK1pO1Lsei+AhwG8mNj3BIBr0e/XADzeUjt/FcA/bdq2BTsfBPBw9PtbAPwpgPNte6YGO1v1TAEIgNXo9wGALwH4MQCfB/DhaP9TAMYttfM3AXyw6edYdqvLQ/9RAK8qpf63UuoOgM8B+Omaru0FSqn/BuDWwu6fBnA9+v06gJ+p1agUNHa2DqXUa0qpL0e//18ALwF4G1r2TA12tgo143b0z0G0KQDvB/A70f42PE+dnV5Ql6C/DcD/Sfz7G2jhhzJCAfivIvK8iGw1bUwGf1kp9Row++ID+EsN22PiH4vIV6OQTOOhoSQisgHgb2DmrbX2mS7YCbTsmYpIX0ReAPA6gGcweyt/Uyn1F9EhrfjeL9qplIqf5070PH9dRH6gQRMLU5egS8q+tvaK71FKPQzgJwH8goi8t2mDPOBTAP4KgHcDeA3Av27WnO8jIqsA/hOAf6KU+k7T9uhIsbN1z1QpdU8p9W4AD2H2Vv7X0g6r16oUAxbsFJG/DuCXAfwwgL8J4CyAX2rQxMLUJejfAPD2xL8fAvDNmq6dC6XUN6OfrwP4Xcw+mG3lWyLyIABEP19v2J5UlFLfir5ExwD+HVryTEVkgJlITpVS/zna3bpnmmZnW58pACil3gTwh5jFpu8XkXgx+lZ97xN2fiAKbSml1P8D8Fm06HnmoS5B/x8Afiga8R4C+DCA36vp2taIyBkReUv8O4C/B+BF81mN8nsALke/XwbwhQZt0RILZMQ/RAueqYgIgN8A8JJS6t8k/tSqZ6qzs23PVETOicj90e/3AfgJzOL9zwH4YHRYG55nmp3/M9GJC2Zx/sY/o0WobWJRlFb1CcwyXj6jlGrdkt8i8k7MvHIAWAHwH9pip4j8FoD3YVYZ7lsAfgXAf8Esi2AdwCGAn1VKNTogqbHzfZiFBhRmWUQfi+PUTSEiPw7gjwB8DcBxtPtfYBafbs0zNdj5c2jRMxWRH8Fs0LOPmaP4eaXUv4q+U5/DLIzxJwA+EnnBbbPzDwCcwyw8/AKAK4nB087AmaKEEOIJnClKCCGeQEEnhBBPoKATQognUNAJIcQTKOiEEOIJFHRCCPEECjohhHgCBZ0QQjzh/wN34doI2VevNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot LSTAT vs. MEDV\n",
    "plt.scatter(df[\"LSTAT\"], df[\"MEDV\"], color=\"black\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "We can create a second-degree regularized (via ridge regression) multivariate model, taking three of the features available (RM, DIS, and LSTAT) as input to predict the value of a home (MEDV). We can use the first 300 points as training, and test on the remaining 206, making sure to normalize first. Take a look at the mean absolute error, mean squared error (MSE), and $R^2$ value of the model, and a plot of how we did with respect to LSTAT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = sklearn.preprocessing.PolynomialFeatures(degree=2, include_bias=True)\n",
    "poly_feat = poly.fit_transform(df[[\"RM\", \"DIS\", \"LSTAT\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ridge(alpha=0, copy_X=True, fit_intercept=True, max_iter=None, normalize=True,\n",
       "   random_state=None, solver='auto', tol=0.001)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create and fit regularized model\n",
    "regularized_regr = sklearn.linear_model.Ridge(normalize=True, alpha=0)\n",
    "regularized_regr.fit(poly_feat[:300], df.loc[:299, \"MEDV\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean absolute error 4.471556587178042\n",
      "mean squared error 41.708742334633925\n",
      "r^2 score 0.322279932371771\n"
     ]
    }
   ],
   "source": [
    "# predict model\n",
    "regularized_pred = regularized_regr.predict(poly_feat[300:])\n",
    "print(\"mean absolute error\", sklearn.metrics.mean_absolute_error(regularized_pred, df.loc[300:, \"MEDV\"]))\n",
    "print(\"mean squared error\", sklearn.metrics.mean_squared_error(regularized_pred, df.loc[300:, \"MEDV\"]))\n",
    "print(\"r^2 score\", sklearn.metrics.r2_score(df.loc[300:, \"MEDV\"], regularized_pred)) # y_true, y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1a1e67e6d8>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnX+MJGeZ379P987A9RhktnaJLJbphoASrNPhwxsCASHCkhO3hy6XCCSWIfIdiEnah+MTScCwUpKLNCfZKAFEgslwGDh1n4FccuF0IrkggxVfooBmDx8/5CO+3EyvHSPWu9g5746x1zNP/qiq7erqt6re+l1d/f1Ir3q6un689U73933qeZ/3eUVVQQghpD106q4AIYSQYqGwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIyzhS5cWOHTumg8GgyksSQsjCc+7cuYuqetx2/0qFfTAYYGdnp8pLEkLIwiMikzT70xVDCCEtg8JOCCEtg8JOCCEtg8JOCCEtg8JOCCEtg8JOCCEtg8JOCCEtg8JOCCEtY6mFfTwGBgOg03Ffx+O6a0QIIfmpdOZpkxiPgc1NYH/ffT+ZuO8BYGOjvnoRQkheltZiP3t2Kuo++/vudkIIWWSWVtjPn0+3nRBCFoWlFfb19XTbCSFkUVhaYd/aAnq92W29nrudEEIWmaUV9o0NYHsb6PcBEfd1e5sDp4SQxWdpo2IAV8Qp5ISQtrG0FjshhLQVCjshhLQMCjshhLQMCjshhLQMK2EXkT0R+Z6IPCgiO962oyLydRF52Ht9UblVBcbjMQaDATqdDo4dO4Zjx46h0+lgMBjg1ltvvfbZYDDAuMGJX4L3Eaxr1HZCCEmFqiYWAHsAjoW23QXgDu/vOwDcmXSem2++WbMyGo201+spAKvS6/V0NBplvl5ZmO6j1+vpcDg0bm/iPRBCqgXAjlpotV/EPSYeEdkDcFJVLwa2/RDAm1X1RyJyA4D7VfWvxZ3n5MmTurOzY9PfzDEYDDCZTFId0+/3sbe3l+l6ZRF1H91uFwcHB3Pbm3gPhJBqEZFzqnrSen9LYd8F8ARcS/Lfq+q2iDypqtcH9nlCVefcMSKyCWATANbX129OK84+nU4HNnUNXRuHh4eZrlcWae+jifdACKmWtMJuO3j6BlV9DYBfBPDrIvIm2wuo6raqnlTVk8ePH7c9bI71DElcshxTNlF16na7qfYnhJAorIRdVR/zXi8A+H0ArwXwY88FA+/1QlmVBICtrS30wsldYuj1ethqYOIX0330ej1sbm4atzfxHgghzSZR2EVkTURe4P8N4BcAfB/AHwC4xdvtFgBfLauSALCxsYHt7W30+32ICBzHgeM4EBH0+30Mh8Nrn/X7fWxvb2OjgfkCwvfh1/XTn/60cXsT74EQ0mwSfewi8nK4Vjrg5pb5XVXdEhEHwFcArAM4D+CdqvqTuHPlGTwlhJBlJa2PPTEJmKr+BYBXG7ZfAnAqXfUIIYSUDWeeEkJIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy7AWdhHpish3ROQPvfcvE5FvicjDIvJlEVktr5qEEEJsSWOx3w7gocD7OwF8XFVfCeAJAO8rsmKEEEKyYSXsInICwC8B+G3vvQB4C4Df83b5IoBfKaOChBBC0mFrsX8CwIcAHHrvHQBPqupz3vtHAbyk4LoRQgjJQKKwi8jbAVxQ1XPBzYZdNeL4TRHZEZGdxx9/PGM1CSGE2GJjsb8BwC+LyB6AL8F1wXwCwPUicsTb5wSAx0wHq+q2qp5U1ZPHjx8voMqEEELiSBR2Vf2Iqp5Q1QGAdwH4hqpuAPgmgHd4u90C4Kul1ZIQQog1eeLYPwzggyLy53B97p8rpkqEEELykErYVfV+VX279/dfqOprVfUVqvpOVX2mnCpOGY+BwQDodNzX8bjsKxJCyOJxJHmXZjAeA5ubwP6++34ycd8DwMZGffUihJCmsTApBc6enYq6z/6+u50QQsiUhRH28+fTbSeEkGVlYYR9fT3ddkIIWVYWRti3toBeb3Zbr+duJ4QQMmVhhH1jA9jeBvp9QMR93d7mwCkhhIRZmKgYwBVxCjkhhMSzMBZ7nTB+nhCySFDYE/Dj5ycTQNV9fe97n8OxY/8YnU4Hg8EAYyo9IaRBUNgTMMXPP/vsEVy69EGoKiaTCTY3NynuhJDGQGFPIDpOfhpnub+/j7OcKUUIaQgU9gSi4+RnFf88Z0oRQhoChT0BU/w8cAXAR2e2rHOmFCGkIVDYEwjHzzvOZaysfADAvdf26fV62OJMKUJIQ6CwW7CxAeztAYeHwMWL1+Hzn38r+v0+RAT9fh/b29vYYIA9IaQhiKpxqdJSOHnypO7s7FR2PUIIaQMick5VT9ruT4udEEJaRmuEPe/sUM4uJYS0hYXKFRNF3tWVbr0V+Mxn3JmlWY4nhJAm0QqLPc/qSuPxrKibjh+PxxgMBtdSCNx66x/TuieENJZWCLvt6kphgR6Pxzh7dl7Ug8ePx2Nsbm5iMpl4KQT+Fu6+++dncsdsblLcCSHNoRVRMYOBK7Bh+n03TBGYCvR+wLTv9XrY378MQIzn7fcBYIDJzMl3AQxir0UIIUWylFExNqsrnT17dkbUATfHS6fzaMRZFVtbplQB5hmmzChACGkKrRB2m9WVXIE+A9fiPvBez+Dw8MNwUwQEOcTa2u9gY8OUKsCs4MwoQAhpCq0QdmB2duje3nw0y9GjHwDwWbhulI73+lnv0/cD2ANw6L2+B/v7vwYA2NraQm/mceCjCHcEIsDp08XdCyGE5KE1wp7MbwFYC21bQ6dzJ9y8Ly8D0PVe771mqW9sbGB7ezuQQuB/4tSpRyABt7wq8MUvcgCVENIMWjF4akOnExX9ouj1rpsbVI3L/2IzWEsIIUWxlIOnSYzH48hB0n5fQhZ5clIv2/BKQgipg9YLux/meHDwIYR9437kzMbGBvb29nB4eIi9vb05UQ+nGzh61HwtDqASQppAK1IKxHH77d/C/v4egGPelgMAgm73MWxvn0hMGWBKV7C6CqysAFevTvcLh1cSQkhdtNpiH4+BS5fuAnAc7iQkgTtAehUHBx+2ygNjXswaeOEL48MrCSGkLlptsbu5Xp5v+OR56HbvtDpHlN/8Jz8BLl7MWjNCCCmPRItdRJ4vIt8WkT8VkR+IyG96218mIt8SkYdF5Msislp+ddMRN5h5cPASq3NE+c3pTyeENBUbV8wzAN6iqq8GcBOAt4nI6wDcCeDjqvpKAE8AeF951czG0aOXIz/r9835YcIDpadPJ6crMB1XZ0y7KdkZIWSJUFXrAqAH4E8A/E0AFwEc8ba/HsAfJR1/8803a5U4zm0KPK1uBHuwPKOj0fz+o5Fqrze7b6+nOhyq9vuqIu5r+Nio40zXKJvRaKS9Xk8BXCu9Xk9HdVSGEFIIAHY0hVZbTVASkS6AcwBeAeDfAfgYgP+lqq/wPn8pgP+iqj9rOHYTwCYArK+v3zwxzewpiU6nA9V3AfgkplExFwH8BlTnrdioiUfdrjuzNGpwtEkTlgaDcDZKvy597HH2FCELSSkTlFT1QFVvAnACwGsBvMq0W8Sx26p6UlVPHj9+3LZeheCmBbgXwIvh3moHwIvR7/8P4/5RPvmDg/ic602asDSfjTJ+OyGkfaQKd1TVJwHcD+B1AK4XET+q5gSAx4qtWn7mE3i56QK2IgLO4wZE41Zkijqu06ne1z6fjTJ+OyGkfdhExRwXkeu9v38GwFsBPATgmwDe4e12C4CvllXJrMwn8IpPF2DK6x4kyuiNOi7J0i+DtJ0ZIaSFJDnhAfwcgO8A+C6A7wP45972lwP4NoA/B/AfADwv6VxVD55mYTRS7XbDg61u6feLO240Gmm/31cR0X6/X+jgZpnnJoRUD1IOnqaKislbFkHYVbNFuYxGIwUOjMIuMr8vI1cIIbakFfZWpxQAssWX26zINHsNN9GY7epKUcv0nY1y4ueEce2ELBlpeoG8pWqLvez48tHIdbO4lvquAp9S4HLi9URkxlr3i4RN+0LqyKcDQhYd0BUzxRXd+RLnK7dlOHxARa6Ezn3ZE/fda2Jv0s9+v28U9n4RFavxWk2G4w5kkaGwBxAxC3tew3g0GqnIxHhuV9TjxbNKKzrt04H/FBI1y3YR4VMLWXQo7AGiLHbHSU4REPe5awWbB0rd7cnCUW5UzLT+3e4jCpyxstiblBqhSPjUQhYdCnsAk1CtrqqurESLl8nFEhY31wrejbTYbYS6LMvYdM+ui2gq7lGdTpmuqzqpckyDkDKgsIcIC6jjRItXnIvFcabndC3AMxoeKBW5YpEwbOQlJzMPsuYV/Chx7nYfSXw6KMt1VTe02MmiQ2FPIEq8pm6Uw8jPfT2c+mzPqD9QKjLRU6ceinVlTI/bjew88rpC8ohzWy12+tjJokNhTyBKvGxKUOCmlrcr7I7zVOzTgHvtvicsUf75/MKaR5zb6mNXZVQMWWwo7AmYfdD2xXe1+FawzTG+tTz19e6mumYaV0hecW5jVAwhi05aYW/9zNMw4Vmlabn77mnudVW7Y/yZp9MMix8FcMX6mmkSM6adNWs6fm8PODx0X7lANyGLx9IJOzArXv1+udcKLqN3+vQIwOMAxnAXozoAzGnsZ44/fTpdWgSKMyHLzVIKe5CkVL1pWVv7KTqdn8AfpxP5KQBXjD/3uTfCXclJvNKNPZcI8PrXu6s3TSbuE8JkUn0qYELIYrGUwh5MinX27ACvf/2fZXLLzKO4cuWPcXj4PPjifeXK8/He9z6H228Hnn3WdIxEXlsVuP9+d5GPIHGLfhBCyNIJu5+JcTKZQFUxmUzwjW/0rP3l8wQPFABvAbA2s8ezzx7BpUsxZ4i59sGBeXvZK90xIyQh2ciSUbZojiTv0i5MKXNVT1gerXDFO+o9kKWv7PeBRx+NFnETZa5053d+fjtNJhMvLTEiV58ihLgivrk5fcr2XadAtWNdS2exmxd1tjF/DzEv4vb+G8cBVlfnt6+suH7+OFEPjwGIuAOqZVF1vnhC2sLZs81wnS6dsJsXdTaFHz4DN4LlEMBeyqsczrzrdg/wyU8C99zjCrzP2hrwwhcC73mPAnjOeKZ+H7jlltnQTFV3QLWsRzxz5xe9nRDiEvUTqfqns3TCblrsGbgXwPvhCrgv5L8G4MVwI1deBpFHI854GHr/U7hhjFN8Ud7YAC5edIV5NHJfXd+7wPWKzTrb/VDJr31t3g9vawVk8feZO7/o7Xlogj+SkKKI+omU6To1kmY2U97ShJmnqrPTy7vd7kwOEVNxZ4yesVpYA7hgNaU/OrXBVQUOtNt95Nqsz7gZrvHrsKabheq3y/Se53OrFDk136Z+ttdjygDSBMpKywGmFEiHm9HRnNZ1vkyTfrmv83nObRa0Ho2ihdrP5x5MKes4T8WkG7iiw+EDxntLkzfGlCjLbxdfKItOppVUP9vrMckXaRJlpOWgsIewseSGw2EKcU8S/qsJYpWUq2b3mpj6mNL8zpbHr91X8H7jOplwuziOY7ynYD2KTn+blInS9npMy0vaDoU9QBpLLuiGyC7qZvFdXb16rdeOzy7pLogRtpR9V1B0SuFDdZzbDPe7a9zfcZ6aa5c4N5RP0QtWJFnsttfjQhqk7aQV9lYPnlYbtvdbCE9MAoBO5xD33ONOFxgMpgnE5lEAn4fIl9weF9P48aNHj8Id4I06WHDp0gcN9zsf7eOOG390rl1czgDYhTv4uwvgzMyAadGDqqZ0DsHcOrbXq3Kwl5CFIE0vkLdUv5h1vCUXN1iYvkTnWHccd0m+aEvdL1fU5Ld3HCewsEeU1X4QsOx3dToO4A7uBv195ns1PXFc1uHwAatB1azE+SPpYyfEBXTFTInzvZrEwLb44jYbUWN2e6QvuyGxdUXacZ7StbX3a1TUjbvQh8kXf1lFNmbGGMztEn3epEHVMmFUDCEU9hniLLl8/nRT5xDtY09XDmLOd6DAf53b7vvwo6Nnpp1Fr9fT4XA413HEPQnMW/buMVyIg5BqoLCHiLLk8rpewu4cEfEsZnNUjH25qq77JOo8BxqMne90zlvEux9q0MXT7/d1dfWzardE326oA5g9puil87iCEyHzUNgtKcJiN2G7nunqqurKSpwQJ4nttC5+Z2UTceMe826La/j1+JTaPIk4jltsBdkk4G1ec5W0hzqMDwq7JXETchzHUcdxPCvc0dXV1TlhdxzH6Mftdh+JtLQ7nUtzQtbtJncCZqt9vpOxjZHvdM5bXue5FPvaC3KUgCctBk5I3dRlfFDYUxD0tfsDoaaBt9FoZJzAE468GI1G3gCnyTf+qWudR3Amp62FP1uuanD2q+8WGg4f0E7nksb5y93OzMZa9+ttu6+9IMc/WcwXhqOTppBmNneRFC7sAF4K4JsAHgLwAwC3e9uPAvg6gIe91xclnatpwq7qivHKyq9qMERwZeVX58Q9aXbj/EDq9HymEMaVlRXvSWA3pWCGhfayOs5tOhw+YOhQZoubf2ak9uMAV3Vt7enMwg6Y2jv9ObpdumNIM0iaLV0WZQj7DQBe4/39AgD/G8CNAO4CcIe3/Q4Adyadq4nCHhUi6Di3zeyXFBOf3WefJprGbD07zlMxLiC3BB8X0z0lHKbcf7bMxqXHu4ocJ/pz+tpJE2iNxT53APBVAH8HwA8B3KBT8f9h0rFNFPZoi3l3Zr8ki30q/DaJwkzi7h+TNaomSnwPtdOZFc9sVvhhZMcSVxzHbb/RaBTb+fghm0njDlUNVjEuPj9tjHBqpY8dwADuckMvBPBk6LMnko5vprBHCeLBzH5Jsxtd4TfP3nT97rYWfHrruNuNG7SdF2NXOLP4zv9fpvpN2y6681lb+4LXzsnnK/uHxJms+WlzhFOromIAXAfgHIC/7723EnYAmwB2AOysr6+X3wIpiZrU4zhPze0bZ8W5ybomEWK0Z4ysMbt2os4RX8w+9jjxziLs2QZSp087cU8jB3rqVFJUT/Cc5X0nmC0yP3W5LIK06YmhFGEHsALgjwB8MLCtVldMUf+00ch1AwS/fMFsjGmIi0ZxHGemUyjG5z77Yzl16iHN4w8vQ9ivu041eXZrtrp0u494uWyK/QHHTV6j1W5HXYOMPm17Yihc2OGu2/Y7AD4R2v6x0ODpXUnnKkrYi/6nFSUM0WF8u3MpZJPF3RdDP7TxOaMwrqxM65s2jLBsUT9yJG4SVlHlWe12nyv0Bxz3v6FLxo66Lfa6r180ZQj7G70v9XcBPOiV0wAcAPfBDXe8D8DRpHMVJexN/ae5/mHT8nln5h7j7ZOQmSx4dyBT5OLM6klxS+jlE/S0TwFu3aImHFVR8nwXkv43dMkkU7fFXPcTQ9GUOniatxQl7E39p7mRH/9ATVExw+EwYv+kNVd3I8RrV+cHcG1EOuqzpxX46dwP0fXdp7HW3TQErjsjiygX564JMhxOo226Xfd90v8yehxkQdWhYur0cTfV+MvKUgh7U/9ps4/ws2GPpklPqjbJyKKjdsIWZHJKgajBy6vXJjmZfojpxdY9X7Z0CcUVv/7DofnzJHEvahC1TYN4i0KRY2dNYCmEve7HvChmY9mTJz2p2kxs2o0Qrl01WZCuO+iCQYwvRwp0kgHq+vdNdYgW/PCPqo4SjN83lW43/r6LCHts6ne17djOKF8UlkLYVZtpBU1FOlqMw5jEY3V1NZCE7DaDSAYzNUYt7jw/USoq1j3JAAX+rUHEDxXYTxTOusU9qSSRd6JSU58uTTTxN5WVtoWsLo2wN5HkiTgHkcfFiYf/g3PFdFfDs1nD/vsoS3M4fMDgqrl8bSHsKNwfiZ8j/lCnOeOTQjMPLePSywzRjC42Fnve2adVjgflqW/bnizatsA5hb1mRqNRZKpb06SnNKSxQqJ+5G4GyPMaHtyNczHER4mc0SjfvZt0LM5y9zsqu5zvRZcbb4weUC1q9mlVFnve+i7Sk4UNtNgp7IVT1sBNXiskTxifKcVxfEjmZQXefa09zPnXbwudo8gJTHHlUG+80fyZL+5R7qwsA6dVWMJ5hayMJ4s6XTvz3/Uz3qzuw4V0M1HYG4LpS53fX5vvxxsXtRPM6x51L258+kTNic3mz+f4GcBi2mP2x1eVsEcP7Popgs0D0E9f66zSUIXA5e30i7bYm+DamRoj71aRK7XWJS8U9oZSTIRF8jniOo8sUTvmEMrZwduoEhR2072EnwDcRUKyirVtp5C8X/xEr7+0/n9VSd5Ov2ghTttRlNn5tcHNRGFvKMXFRMcnIgsL/8rKyrUIm6kLZdf4RTeNAcSlSUgS9ugngOjB3dXVLILtr81qsrJ/qm5WyqKeBg4T49/Lxu7pJ1toZlHimsa1U7Z139QJjWmgsDeUKkbpk2Pi/RKdPne+3lECd5B4nahOK66TGw4f8Nw9vkvnU9d8o9FiezVwnqlLyHGeCsQyFyHq0+I49TzKx4lg1Tnk4zqCNFZy2RY1LXYKe2lUMUqfPIvVL2ah63YfMdQ7Stx2Y68RZy3GdXLxom/ujDqdz0ReO/u6ssnF7/SCojYazS7KnbUDiBLNpohUkpWdxgov26Jugr8/LxT2hlLF4g32Fnt8FMtsved/FO5A1LvVcZzARKrp30nWYpx4x3VOIqKrq5/VYDz9qVMPaXixccdxZq6fnEMnf3Hz6kRnsxwOZ8XacdxisnbjhKgpbgWbDsbWtVNFZ7Xok68o7A0my6NymmNsM0a6vnZzKJ/pekX/KOI6uaBA2xTXNz+M7TTNA8DB8rQeORKVNsG+5Jlp63cMcZ1Qv1+soOahyA6mDRZ12VDYW0QWKz8ozI7jzK3cFCeESQJZ1D2Fo2GmHchIV1ZWUgk7AO10OsbtQTdXUOzW1p72InAOZhbrqDPNsC+KSZ8X6QLJQxnhkYtsUZcNhb1FFOGXj7L4TdvLHgewWzc2najHlSwD0+aIEzvRLzs3jv9vKGrQMg+2HQgFuxgo7C2i6nwXZV8vqeOwH/y1K0UOTKsGJ2rNW9dJPva8xcbqdidWmYtv7RfrUos/H10sxUFhbxFVRNIEcaf472o4j0ye6wWfDOJEWLV4i9124la2+zKL2mikurZWrKjbiHDSOILjVC+yTYngaQMU9hZRRSTN9FqmHOqXU+WwDounyWdvKl0vzaL9coHJZTadQXXtOL1mMaJuW8W4QVc3L0/1IhtVnzomBlUZ41/GtSjsLaOqL2SUMNhmpDSJZxrXimlQNXnZQHMJi3bVTz4+7rKCz2QW9ZiMDHPEDbzWESYZd83qY+6rNJDKuRaFnWQi7w8/jxvFcRzjYiNra2vW54jr/LKOHRSTtC0YVnpBw+vKRpW0bpKkMMmqLfao+vi+/iqpsmMv61oUdpKJvP7QrAOfvV4vdey6qYQJinKU5R/3YyvC8jLXNSj0B+pOtDr0/nZz2mQZ2Ixa19UvKys6l4enTB973BNE1VQZhFDWtSjsJBN5IxiiLJXwFz2YlMy3gouIhpm9l2RffZJI57G84kJHbZ48bAgP3tqEYzpOdaGHTRo4pcVOYV9q8oTDRWdsHM65M8IzTKMmGNmW8EBplIXe6XSs3Sp53Dd5BoBtLLvkmbTmUuWgZZNCHeljp7CTHNj4pLPOLvUFPDyTdnV1dWbSVZKo2i8Vl83yyhuy2e12EzufOH96XCnLWo4L/WzK5CRGxVDYSYlkFT7f6on70dicOyjMaXPb21heRU6y8s8VrltS6gFTKctabpJl3mYo7A2g6rzYi0SS8Jks7nC2xqznDro68q5GFUXWjivJHTWbesEs3kEfelz2yCJpki+9zVDYa6aOyTBNwUYI44Sv701qCsaxD1MsV5TGYo/a13GcXJ1ykZOsoureJCu5KWmE2w6FvWaqHIFvErYdWpSPfXV1NXd2ySRRDZ7L1mWSpVO2CbXMUoIDq03xX9NirwYKe81UnbirKaTp0KIWxiiiUwynLY5a/CONyyRPp1ykz70u4yCuE2nS00ObobDXzLJY7GG3i42VmUSVnWJal0lW90yaDsR3A/n3nPfJoQhshLspTw9Nq0uRUNhrZhl87GnywqTp0KruFMOdk+0M2CLdQ34Jhm2a6lbX92eRXC1tfnqgsDeApvwoyyJKgPNamaPRKDZOvWzSzIJNM8gaN2kqeL4i7yPL98903CINji5SJ5SWwoUdwD0ALgD4fmDbUQBfB/Cw9/oim4u1Sdjb+shnQ5z4hYUhjciYBlZXVlYq7Rht3SZprfgky70od1PWJ8ao4xznqcxiWbWBs0idUFrKEPY3AXhNSNjvAnCH9/cdAO60uVhbhL3Nj3w22LpM0opME8Yn0vjE09YzznIv6h6ztmF0+Odtmb7rdbgkabGndMUAGISE/YcAbvD+vgHAD23O0xZhb/MXyAbbH21akbEZPC3bCswTh26X5yW/4MW1QdYB6Ljjsjyd1tFJt9ngqkrYnwx9/oTNedoi7G1+5LPFRmDTikySGFRlBYbvbTgcFmKxh88PTBcTse2kktqgaIs9qxDXFfbbVhdp44QdwCaAHQA76+vrFTRB+Sy7xW5LWrEoS7TKvJc8VneWTqqszq/oTrMJbrU2UZWwL7Urps2PfEWSRSzKcDPY1jXuCSQuxDOLSyir8JXprirSzVV3hFPbqErYP4bZwdO7bM7TFmFXbe8jX9EUKRZlWYFp0iEUdS9ZO6m4NmhSmG0TIpzaROHCDuBeAD8CcBXAowDeB8ABcB/ccMf7ABy1uVibhJ1UTx53QZxfu8gOw1Zc41w7ccdFtYFNnp0qhb8uV0yTOrciKcViL6pQ2Elesvxw4yJdiowtT9PxpElYZtMGdQw81+U2i6tPW2d9U9hPy2e0AAAG5UlEQVQJCZE08FlUbHmWweIky92WJCEt2oJu4kB3mwdsKeyEhLBJE1CEpVdGDLktSaJWtAVdxxNCEm3OrJpW2DsgpOWsr6/Hft7v97G9vY1+vw8RufZ+Y2OjkOskXT/LcePxGIPBAJ1OB4PBAKdPn0av15vZp9frYWtrK1fdojh//nzs9o2NjVRtGr6f8Xicuk5F3+NCk6YXyFtosZM6iPNn5wlZtLlOGTHkcQOoUT7vIgae0/j0444tqt3KOk8TAV0xhMxjioopI+e57eCuaXar7aBwVl9yUQPPaaJwykg9UfQ9NvEaYSjshFhQ50BbXsuySl9yXDslCZxtGy+Sb7yup4K0wi7uMdVw8uRJ3dnZqex6hETR6XRg+u6LCA4PD0u99mAwwGQymdve7/ext7dX+vFpyNNOtsdWeT95qauuInJOVU/a7s/BU7KU1DnQljTwmMTW1lbsQGmR2LRT1MCnbRtXeT95yfu/q4w05n3eQlcMaQp1DrQV4Qaqys+b1E5xn6edsFW13zoLdbnwQB87IXbUJSZVdCpFJ/SKOpdNPPsiCLYti+Jjp7ATUgNlCl6V4rNIA59FsQhRMRw8JaRlVDnAt0gDn4sMB08JaQBFzKTMSpUDfHEDn3W2wbJDYSekYMbjMTY3NzGZTKCqmEwm2NzcrEzYqoz4iUodAKDWNlh26IohpGDqdk/4Hcv+/v61bb1eL1P+m6zU3QZtg64YQmqm7ljntAm4yqDuNlh2jtRdAULaxvr6utFarTLL4MbGRqVCHqYJbbDM0GInpGAWaSZlWbAN6oXCTkjBNMEVUjdsg3rh4CkhhDQcDp4SQsiSQ2EnhJCWQWEnhJCWQWEnhJCWQWEnhJCWUWlUjIg8DmB+1kI9HANwse5KWMK6Fs+i1BNgXctgUeoJuHVdU9XjtgdUKuxNQkR20oQP1QnrWjyLUk+AdS2DRaknkK2udMUQQkjLoLATQkjLWGZh3667AilgXYtnUeoJsK5lsCj1BDLUdWl97IQQ0laW2WInhJBWspTCLiJ7IvI9EXlQRBqVlUxE7hGRCyLy/cC2oyLydRF52Ht9UZ119Opkque/FJH/67XrgyJyus46+ojIS0XkmyLykIj8QERu97Y3ql1j6tm4dhWR54vIt0XkT726/qa3/WUi8i2vTb8sIqsNrusXRGQ30K431V1XABCRroh8R0T+0Hufuk2XUtg9/raq3tTAkKcvAHhbaNsdAO5T1VcCuM97XzdfwHw9AeDjXrvepKpfq7hOUTwH4J+o6qsAvA7Ar4vIjWheu0bVE2heuz4D4C2q+moANwF4m4i8DsCdcOv6SgBPAHhfjXX0iaorAPyzQLs+WF8VZ7gdwEOB96nbdJmFvZGo6n8H8JPQ5r8L4Ive318E8CuVVspARD0biar+SFX/xPv7Kbg/mpegYe0aU8/GoS6XvbcrXlEAbwHwe9722tsUiK1r4xCREwB+CcBve+8FGdp0WYVdAfw3ETknIpt1V8aCv6KqPwLcHz+AF9dcnzg+ICLf9Vw1tbuMwojIAMDPA/gWGtyuoXoCDWxXz2XwIIALAL4O4P8AeFJVn/N2eRQN6ZjCdVVVv123vHb9uIg8r8Yq+nwCwIcAHHrvHWRo02UV9jeo6msA/CLcx9031V2hlnA3gL8K93H3RwD+db3VmUVErgPwHwH8hqr+Zd31icJQz0a2q6oeqOpNAE4AeC2AV5l2q7ZWZsJ1FZGfBfARAH8dwN8AcBTAh2usIkTk7QAuqOq54GbDroltupTCrqqPea8XAPw+3C9lk/mxiNwAAN7rhZrrY0RVf+z9gA4BfBYNalcRWYErlmNV/U/e5sa1q6meTW5XAFDVJwHcD3dc4HoROeJ9dALAY3XVy0Sgrm/zXF+qqs8A+Dzqb9c3APhlEdkD8CW4LphPIEObLp2wi8iaiLzA/xvALwD4fvxRtfMHAG7x/r4FwFdrrEskvkh6/D00pF09P+XnADykqv8m8FGj2jWqnk1sVxE5LiLXe3//DIC3wh0T+CaAd3i71d6mQGRd/yzQqQtcv3Wt7aqqH1HVE6o6APAuAN9Q1Q1kaNOlm6AkIi+Ha6UDwBEAv6uqjVk6XUTuBfBmuBndfgzgXwD4zwC+AmAdwHkA71TVWgcuI+r5ZrjuAgWwB+Af+j7sOhGRNwJ4AMD3MPVdfhSu/7ox7RpTzzNoWLuKyM/BHcjrwjUQv6Kq/8r7fX0JrmvjOwDe41nEtRFT128AOA7X3fEggH8UGGStFRF5M4B/qqpvz9KmSyfshBDSdpbOFUMIIW2Hwk4IIS2Dwk4IIS2Dwk4IIS2Dwk4IIS2Dwk4IIS2Dwk4IIS2Dwk4IIS3j/wMlBdtWFkS60gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot predictions w.r.t. LSTAT\n",
    "plt.scatter(df.loc[300:, [\"LSTAT\"]], df.loc[300:, \"MEDV\"], color=\"black\") # true labels\n",
    "plt.scatter(df.loc[300:, [\"LSTAT\"]], regularized_pred, color=\"blue\") # predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a neighbor, discuss the interpretation of the mean absolute error, mean squared error, and $R^2$ value of the model. How well did our model perform, and does the plot corroborate our results?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Perhaps we now want to determine the optimal coefficient $\\alpha$ in a fourth-degree polynomial model. Let's use 10-fold cross-validation to determine such a value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = sklearn.preprocessing.PolynomialFeatures(degree=4, include_bias=False)\n",
    "poly_feat = poly.fit_transform(df[[\"RM\", \"DIS\", \"LSTAT\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best alpha \t average validation MSE \t test MSE\n",
      "0.2 \t\t 24.9 \t\t\t\t 11.7\n",
      "0.2 \t\t 25.7 \t\t\t\t 7.0\n",
      "0.2 \t\t 25.1 \t\t\t\t 9.6\n",
      "0.2 \t\t 24.6 \t\t\t\t 27.1\n",
      "0.2 \t\t 24.4 \t\t\t\t 17.6\n",
      "0.2 \t\t 24.6 \t\t\t\t 21.9\n",
      "0.2 \t\t 24.9 \t\t\t\t 11.8\n",
      "0.2 \t\t 16.4 \t\t\t\t 98.0\n",
      "0.2 \t\t 26.1 \t\t\t\t 21.3\n",
      "0.2 \t\t 24.4 \t\t\t\t 14.5\n"
     ]
    }
   ],
   "source": [
    "alpha_grid = np.linspace(0, 10, 51) # hyperparameter grid\n",
    "kf_test = sklearn.model_selection.KFold(n_splits=10)\n",
    "kf_val = sklearn.model_selection.KFold(n_splits=9)\n",
    "\n",
    "print(\"best alpha\", \"\\t\", \"average validation MSE\", \"\\t\", \"test MSE\")\n",
    "for train_val_index, test_index in kf_test.split(poly_feat): # outer (testing) loop\n",
    "    X_train_val, y_train_val = poly_feat[train_val_index], np.array(df[\"MEDV\"][train_val_index])\n",
    "    X_test, y_test = poly_feat[test_index], np.array(df[\"MEDV\"][test_index])\n",
    "    \n",
    "    alpha_mse_grid = np.array([0]*len(alpha_grid))\n",
    "    for train_index, val_index in kf_val.split(X_train_val): # inner (validation) loop\n",
    "        X_train, y_train = X_train_val[train_index], y_train_val[train_index]\n",
    "        X_val, y_val = X_train_val[val_index], y_train_val[val_index]\n",
    "        \n",
    "        for n, alpha in enumerate(alpha_grid): # hyperparameter grid search (loop)\n",
    "            val_regr = sklearn.linear_model.Ridge(normalize=True, alpha=alpha)\n",
    "            val_regr.fit(X_train, y_train)\n",
    "            val_mse = sklearn.metrics.mean_squared_error(val_regr.predict(X_val), y_val)\n",
    "            alpha_mse_grid[n] += val_mse\n",
    "    alpha_mse_grid = alpha_mse_grid/kf_val.n_splits # compute average\n",
    "    best_alpha, best_avg_val_mse = alpha_grid[np.argmin(alpha_mse_grid)], np.min(alpha_mse_grid)\n",
    "\n",
    "    test_regr = sklearn.linear_model.Ridge(normalize=True, alpha=best_alpha)\n",
    "    test_regr.fit(X_train_val, y_train_val)\n",
    "    test_mse = sklearn.metrics.mean_squared_error(test_regr.predict(X_test), y_test)\n",
    "    print(best_alpha, \"\\t\\t\", round(best_avg_val_mse, 1), \"\\t\\t\\t\\t\", round(test_mse, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There seems to be consensus that the best value of $\\alpha$ is 0.2. Let's see how our model performs, with $\\alpha = 0.2$ and degree = $4$. For comparison, change the value of $\\alpha$ a bit, and see how the MSE adjusts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean absolute error 4.919735838436536\n",
      "mean squared error 46.464565629175716\n",
      "r^2 score 0.24500316245755627\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1a1e99d8d0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXuMJNd13r/TvdOUelYCtbWbgBE93VIgJBICixE3igAKgsJ1DFlWHDmQDS9HAQULnqRlEwycxKK8QBIHmACkEYmOBEYeW6QYzFgPy3FkGAEShhITxn9ImA1piQLD0PbMrBgKWj5DLoekltMnf1TVdnX1vVW3Xl3V1d8PKHR3TT1u3en+7rnnnnuuqCoIIYS0h07dBSCEEFIuFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZx+Z5s5MnT+pwOJznLQkhZOE5f/7806p6yvX4uQr7cDjE7u7uPG9JCCELj4gcZDmerhhCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZFHZCCGkZSy3sOzvAcAh0Ov7rzk7dJSKEkOLMdeZpk9jZATY2gMND//PBgf8ZANbX6ysXIYQUZWkt9nPnJqIecnjo7yeEkEVmaYX9woVs+wkhZFFYWmFfW8u2nxBCFoWlFfbNTaDfn97X7/v7CSFkkVlaYV9fB7a2gMEAEPFft7Y4cEoIWXychF1E9kXkuyLysIjsBvtOiMh9IvJ48PqmaosK7OzsYDgcotPp4OTJkzh58iQ6nQ6GwyE+8YlPXPnbcDjEjkPs4vo6sL8PjMf+67xEPfoc0bLa9hNCSCZUNXUDsA/gZGzfHQBuC97fBuD2tOtcf/31mpft7W3t9/sKwGnr9/u6vb2d+35VYXqOfr+vo9HIuL+Jz0AImS8AdtVBq8NN/HOSEZF9AKdV9enIvscAvF9VfyAi1wB4QFX/WtJ1Tp8+rXlXUBoOhzg4yLSICAaDAfb393Pdrypsz9HtdnF0dDSzv4nPQAiZLyJyXlVPOx/vKOx7AJ6Db0n+tqpuicjzqnp15JjnVHXGHSMiGwA2AGBtbe36rOIc0ul04FLW2L0xHo9z3a8qsj5HE5+BEDJfsgq76+DpDar6LgA/BeCXReR9rjdQ1S1VPa2qp0+dcl6LdYa1HHGIec6pGluZut1upuMJIcSGk7Cr6pPB60UAfwjg3QB+GLhgELxerKqQALC5uYl+PD4xgX6/j80Gxi6anqPf72NjY8O4v4nPQAhpNqnCLiKrIvKG8D2AnwTwCIA/AnBzcNjNAL5eVSEBYH19HVtbWxgMBhAReJ4Hz/MgIhgMBhiNRlf+NhgMsLW1hfUGxi7GnyMs61133WXc38RnIIQ0m1Qfu4i8Fb6VDvhJw35PVTdFxAPwVQBrAC4A+DlVfTbpWkUGTwkhZFnJ6mNPze6oqn8B4J2G/c8AOJOteIQQQqpmaWeeEkJIW6GwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy6CwE0JIy3AWdhHpishDIvLHwee3iMi3RORxEfmKiPSqKyYhhBBXsljstwJ4NPL5dgCfUdW3AXgOwMfLLBghhJB8OAm7iFwL4KcB/G7wWQDcCOBrwSH3AvhwFQUkhBCSDVeL/U4AvwZgHHz2ADyvqq8Fn58A8OaSy0YIISQHqcIuIh8CcFFVz0d3Gw5Vy/kbIrIrIrtPPfVUzmISQghxxcVivwHAz4jIPoAvw3fB3AngahE5FhxzLYAnTSer6paqnlbV06dOnSqhyIQQQpJIFXZV/ZSqXquqQwC/AOAbqroO4JsAPhIcdjOAr1dWSkIIIc4UiWP/JIBfFZE/g+9z/0I5RSKEEFKEY+mHTFDVBwA8ELz/CwDvLr9IhBBCisCZp4QQ0jIo7IQQ0jIo7IQQ0jIo7IQQ0jIo7IQQ0jIo7IQQ0jIo7IQQ0jIo7IQQ0jIo7DnY2dnBcDhEp9PBcDjEzs5O3UUihJArZJp5SnxR39jYwOHhIQDg4OAAGxsbAID19fU6i0YIIQBosWfm3LlzV0Q95PDwEOfOnaupRIQQMg2F3YGdHWA4BDod4ODgAQBnZ465cOHCvItFCCFG6IpJYWcH2NgAJkb6EMDvBO+/dOW4tbW1+RaMEEIsLJTFHrWch0P/c9WcOxcV9ZBVAP/myqd+v4/Nzc3qC0MIIQ4sjLCHlvPBAaDqv25sVC/udg/LGkQEg8EAW1tbHDglhDQGUTUuVVoJp0+f1t3d3VznDoe+mMcZDID9/ULFauR9CSEkRETOq+pp1+MXxmK3Wc5Vj1lubgL9/vS+ft/fTwghTWRhhN02Nln1mOX6OrC15VvoIv7r1pa/nxBCmsjCCHudlvP6uu92GY/9V4o6IaTJLIywV2k5p0XbMIUAIWSRWJjB06qYjVP3ewJhoxFPIQAAKysfwxvf+Fk8++xxrK35vQZa8YSQqmjt4GkZmCxvU5z64aEfvw6YUgicxeXLn8Mzzxyfa9glIYS4sjQWu8ny7vf7ODy8BEBmjhfxfeqdTgfTdbQHf/bpNAx/JIRUBS12C7bkXZ3OE8bjT5y4BMCUKsAchsNUMYSQprA0wu4n6ToL3+I+Cl7PYjz+JICXYkeP8corfwAA2NzcRH8qHMes4EwVQwhpCksj7CdO/Ar85F1D+I89xCSZ1z0AxpGjO3jppY9gZ8fPsb61tYXBYAARged9Gr3ea1PX5oQlQkiTWBph95N2rcb2raLTuR3AhzBbFatXBlDX19exv7+P8XiMp5/+d7j77mPwvMmRr399ZYUmhJDMLE3a3mefPW7cPx5fC8A8gJzkN3/55cn7Z57xI2MAhj0SQuqnNRZ70iSjnR1/v4nBQOB5M3l5Adj95mkhkoQQUietsNjjk4zC2PKQX/zF13B0NPuoE9/4ceMkJZvfvK6EZIQQ4kIrhD3Jgr506RJ+9KNZN0ynM8bWVmfKdXLunC/OabNJ19bMqXwZGUMIaQKtmKDU6fiLb8QRAVTHMHucxlDN54lKS0NACCFlspQTlJJT+tr8I/n9JkzlSwhpMqnCLiKvE5Fvi8ifisj3ROQ3gv1vEZFvicjjIvIVEelVX1wzSSl9e73/hukYdQB4CZ73aefrmwZmmcqXENJUXCz2VwHcqKrvBHAdgA+IyHsA3A7gM6r6NgDPAfh4dcVMxmZB/8mf/E/86EdnMf2YYwD34rd+6287XbuutVYJISQvqcKuPpeCjyvBpgBuBPC1YP+9AD5cSQkdMVnQW1tDzE5K6qDT+XvOi0+7hjam5XSfJ8wfT8hy4+RjF5GuiDwM4CKA+wD8OYDnVTWcW/8EgDdXU8T8HB39FeP+8dhe1FCgRYBjx8zRL8B0aGOTrPowi+XBwQFUFQcHB9jY2KC4E7JEOAm7qh6p6nUArgXwbgBvNx1mOldENkRkV0R2n3rqqfwlzUG3+2Sm/VGBBoCjI/u1owO2TZqwZMtieY6zpwhZGjJFxajq8wAeAPAeAFeLSBgHfy0Ao1qq6paqnlbV06dOnSpS1sxsbOxjNnPjS8H+WUwCbSI+ealJE5YuWG5q2x+nSS4lQkg+XKJiTonI1cH71wP4CQCPAvgmgI8Eh90M4OtVFTIvd931XoxGD6HbfQLAGN3uExiNHsJdd73XeLyL9plCG23hlp3O/IVxNn988v4oTXIpEUIKoKqJG4AfB/AQgO8AeATAvwj2vxXAtwH8GYDfB3BV2rWuv/56bTKDgaovaeZtMDCft72t2u+bz+n3/b/PnrOtg8FARUQHg4Fumw7Kwfb2tvb7fYXvGlMA2u/3na5ve37bcxNC5gOAXU3R1+jmfGAZW9OFPY9Ah4xGDypw2UkYi4iv23PkazREzM8uUkqxCCE5ySrsrZh5WhbReHgA6Hb9V9vM0kkEjeLzn1+DzbMVd/FUPcAZzR+/HyzE6hL+mDyDlxCyMGRpBYpuTbfYs2C27o+cLHYRmbLWw01KMo23t/17iqh63ou6svIxp96B6ZnSeiqEkOoBLfZ0yoj8uPXWS4YImg7i6QtM6X+LDHCmER8AfeaZ47h8+XPw13v1sfUO2pwDh5O2yFKRpRUoujXBYnexSqMW72Awa7Fub29brXN//54CR9rtft86cFqVj90+ALxXSe9gEah6TIOQqgEHT5NJi/wwCb+I6mgUvcYgEEq7gKYJR1VRMbYBUL/BmQjbYIlCXfz/16zra5nqgCw2WYV96VwxaZOJTJOUVIHPf37isvEn+/w6TJOfgHMYDAbY2tpKyUezDmAfvutmP/hcHLs3Z/Lg/X4fm7bloebMPCZEFZ20RcjCkaUVKLotgsUOjC0W7+SYiQV49orbBdhTz7vFqQxJ7qA0N1Cea/d6l9Xzbim9d1CUeQ3W0mIniw7oikkmWVS3VeTAKuwi/nGe92Ig5heD7UhFDnQ0etDJP9/tfj/xHkWFrmjjMC/mNSGKPnay6FDYY5hEziZ8vmV31jow6nn2CUyA6sqKaq9nF+aJwNgGXucjdE1hnhOiqhrTIGQeZBX2Vqx5aiPL2qQ7O8BHP7oPYA3AZQA9ADJ1jL+GavZyDAZ+jvjhcIiDgwMAewCGzueL+Hnm28ZwaE6LHNYXIcRnKdc8tZFlkYyNDcAX2w6AqxAXdSCfqAOTgdmDgxvgp7QfAOYsx0baOvMzaUlDQkh+Wi3srul0b73VLV1vXtbWwmiPewCcgt9ohA1HssD3esAHP9jOVLptnhBFSJ0cSz9kcVlbM3f1oxbwzg7wzDPVlUHEF2a/l2Ba71vgi/tsDwEAVlaAe++dNDxhKl2gHQK4vt6O5yCkSbTaYnfp6hfLu2Wztif7VX1hti2x52MWdQB46aXmrM5ECFkMWi3stq4+MMkbcnBQZFTSJsjT+w8PJ5kiTQwGk4ySrnBuDSHERquFHfDFfX/fjyrxIy2mF3uOzsi0M0aWwU4TR0e+vzzOyorfg9jc9N/H/+Z55uudOFGoOKkwaRYhi0vrhT3ObC50U2qAVwE8hcl0/48CsPlS4ha/uQcwGAB33z0t1J4H3HPPxMcssQ6ACPDzP29uEF54obpB1J2d6cbv4OAAGxsbFHdCFoUsQe9FtybMPDXnQp9ODeB/Nh1zKTaZ5pICn42cux98fjl23Cupsz+TZmF6Xv6JS3lmodY5BX9RZs0SMk/AmafJ2ETLZTtz5guBAI8DEZ9uAERERdYVeGVKgI8de82Svlcj1zPnqBFJytiYLHx5c7FUvRBI2eUlpO1Q2FMw5Q1x2Y4fPz51ndXVVYtV77ruaXJ6guh5SYtsJwlf3lwsdVnsXEybEDMUdgeieUO63W6qqIcW7GAw0NFoZBE+k6tm2vKOkiTWcdHe3vYzNNqO8zzzc2bNxRLWS/SZwy1MmlVlzpWkXPJp92MuGNJmKOwZ8TM6ml0P2ba9RJH2vGnfsf3YsYa+/mjyMH/dUltK4bHRas9iAU/usaeTsYabrjRooahXmSXRZfUn0/2YvZG0HQp7DBdLbjQalSDuSRkbx9rtJon5rIhF3R6THoK98fC8F2ee1/NumbH0ba4bz7tFTYPD0RzzZbpobFk3Z91TlzQ+lhG/H/Otk7ZDYY+QxZKLuiGqsNjdtrECn51y/Uz3KM6q3Wo/Mj7vysrH1PNenIkyiTd4SUv9hfjlmI0gyjqo6rLQiL//spoileL3q2uwl5B5QWGPkMWSKy7syT52921v6pr+wh9RcbtoPW+ST35aeOPPax5Ati/OHZJk1WcJU8yz7mzUcqfFTpYNCnuENEsuabAwv7jbBzndtnFwnaS4ebO4+j5x0zk3TdXLRAijjYC5J+B5L+po9GCw6pPdx59l5ae0Qd0kXzt97GQZobBHSLLk8oY9pjcC2VZHMm+vKPCaVdziVvnKysdSltzbnxpjmLhU0noYYz1z5lGH48ybzWBOs9jt0TFj9bwXLXMCksdSOPGJLDIU9ghJllwxt4vfOJivsWcVunAGadKEo/TtaOp+nuddEbEkizp0Y/T7ffU8T+0unenzOp0LuctqcnH7a8bOHhu18NNCQeO9gfR1Zv1lC6PXWFmhuJPFgcIew2bJFXW9iIix4eh2/6HOWrhHeubMo5EyucWxm7eLM8IeknzNiY96dfWXEhqB6BYOXuYrq+ukLM+LDuraUyiYru0yW9V2PVv8PyFNg8LuSBkWu+psw+Fbw7MDmJ5FRbIL/CsaD/8LG6vRKK03sKdpvYp4g5Q32sfkY883aGre0vzx0UYl6TqELAIUdkdM1nZoxbtY81EXSJSkc81hlu5iNivQ0TEDl+uEbhwXa92/j2/dZ/exdzp+QzNdN/lE2rYlzQ2IuoHyCjv98qQpUNgzMPG1n1U/qdeRdrvfD5J9Tazw0WgUWOLTQm2KvEjqCXS73SmX0OT+N0WiTmyDptFtMjs1jCN3E8W9oDFzidwZK3BRR6MH9cyZLzieM7tFxT3/oGn2LWqx53HFMCEZaRIU9oyMRg+qKURwNHpw6jiXWOnt7W1jA2DaVlZWtNfrzTQU7tb0pKx+qGPaeWG8+bZm8Zv3epeD6+eL0+92o3VtP87z3HzreRoUP9dO/LmKxdoTMk9KF3YAPwbgmwAeBfA9ALcG+08AuA/A48Hrm9KuVaawl9VNtoUIdrvfnzrOJSY+b/jktFVvC1lMEsUXE84bz9RP9ntcDoT9s5kahXALSetVdLuz0St5t3jUS9bvS9YEannvQ4gLVQj7NQDeFbx/A4D/A+AdAO4AcFuw/zYAt6ddqyxhL7Ob7DLjUjXdYi86GDvZbsrhc1e1T0763IzQ+L2U7D0DtxDJWbH2/2duPYXV1fJcMp3O5H008sYFum9Ik6jcFQPg6wD+LoDHAFyjE/F/LO3csoS9zG6yq8WeNrsxadDUJTVwuHmeFwjvRWfx7XbDhuVs5LyxAv9PgVenjhXx3RTHj8dXeXLZsjYGqmfO+D0KX9SLzsotvsUHdE1sb6uKmBuhJGGn+4ZURaXCDmAIf/XnNwJ4Pva35yznbADYBbC7trZWykPm7SabcPWxqybPbrRZ7N1uV0ejkbObZnV1NXJsNGwyWRQnaXfT/eAiaplRmibc2YT9zBlTpE72xiFpe93rsh0vMh0zb8owmZT7Puk7Vub3kpAolQk7gOMAzgP4B8FnJ2GPbk202FU1kgvFj4oxiXoaST72fr9/ZYGOSTbFLO6ZqCVuf27fMnYTOP953RsOf3vNWZgHg+zhiysr2QdQRY702LFs54QibhpQTbt/mB7ZhO15u1363EkxKhF2ACsA/guAX43sq80V01Rfpp+HZV3jk5OivvgQd3FPjkaJPnc23/RR7D5p/u+4oNvXaQ17PNnK89oV/7/IS5mEenX15SlBTruvSN4InEtX8vKY///pcwmYyoDkoXRhByAA/gOAO2P7fxPTg6d3pF2riVExZbK9rQYR9qfyx3OD2yz8WV/9nlUk4s+dxULudJ7R6QbINjDqGlsfbn4EjT8L13bMbCPxjnc8EfHFhz0DV7fN0cz/IW6Nx+stu6iHkUHJ6YCj30tbA1NGKoMmfv9JdVQh7O8NROY7AB4Otg8C8ADcDz/c8X4AJ9Ku1cQ49jKxu0L2jGLgFvduj9qJ+/ldZ7H2eqrdblysX1Y/XcFkX7+fJ4LG7wn4KzhlEU7bc7rce29G6PwB4tljwx5Ocv3E3VLTqzi5LuCRdI8iNLXHSqqj0sHTolvbhT1JhEeWcIx0l8yeVcxC/31c3K+66gWDIPo51z3Pbkmurr5stAKzC/vZYDwhi7AnbaHQmsox1l7vv6euyBR/pqRwRj8RWdh72IuI+mRMwsVKrkrYGX2zfFDYayRJhG0LP6TnpbEtuDGxIM3+e1MislsSLXqbIZqcDth0nYMggqgsYfddLbYQzU7H3KB2u3bxdUnlO3GXmQev06xkW+MRxteHuW6yulIYfeOGy3rHiwKFvUbMS8eNA1GYXaJO1W0QNVy31LT+p8k1YG8sbA3PRGBMJPvfzftDa9ltslXywG2ncyF4pmwNTJr4uvip0wZzk6zkNF9/XJRdYuxVabG70LZVtSjsNTKJJTeFJs4uUReeE/8C9no99TxvxtJwyVeTdFySgCYJoN9gxa3llxME3/dTh+Juz8I4Dhob03J/0WM+G5Q/3wQnF8GziXxaryNLioGkbJThtVx0hz72dFx/K4sChb1m/Knz+8Yfbnw2a/Qcly7jaDQyflnj/nubtWIf3L2cGLs/abD2dNJruEnTQjHDmG+7rz0acpm0Xmy44EfeyU1HiXVrE0p/4Di5N5FFJ1zGHFyvZ55c1R7XQ1HScjstGhT2BpDkMihCFivE9CO3zbIFzqZ2U+PXmxZk2/MeBeW2Cdle7FmKrhdrK8fkPqFbKyqI9vI9lVgmFys5i8UO5PeTt831UBRa7BT20qnKB1rUCrENqmb50pvXi92zCqp/jtki9l080etkTzI2vV00+MOjA82zPYzkMQB7D8ElqVieRVTyxrhXKWSLGDPftoaOwt4Akrr2RbrKRX+8aRE4aQ2EPW2COXJndfWXpupkOs78QUMMv6uwmyYvhQJ+VkUOdOK/jw40mxug5DEA899cxC4pxUDZwl6V62GR/fltck1R2BuCSciKWhAuVoh7orJZyz2tgUiO4HFf59XeQLi5YnxhHOvE9x4XcL+xmy2vfRJUXLx8y9/c0Hiem9glhSWWHbJo640VtdgZgdMMKOwNpayucpJwmwRzZWXlSoSN53m6srKiNgs7LfmZy1qwLtairS7yLDJiWtw7vPdsfexZxPQg6E1NN8SmbJn+ilJuYmcTRc9LttrzuDuyZCnNAmPmmwGFvaHMY5TeJSa+1+sFrgqT4JgzF04GGM3WsW2zNVr2BuImw3T+l9UfxBxPLZwxvV203jvaEPo+fbMP3jYA7Z+zp8CRet6LQa72dLHzZ6/OHrOy4hbbntXdkRSWWcQvvmgW+yKOB7hAYW8o8xild7eo3VaNUrUNAE7PfDVtSW6mpLqICmk8O6ZdCMep9540TmYXTpYGNk3sRiOz+Gdd17XMUMq8fvFF8rHPq6x1NB4U9oYyj1F6F4s9ySURRrFMX9MmFnvqed6Um8c0qSprXSQ3TvYw0uhArOd5DsnRzGkZXH60SQKSZNGH13UV9iydOZf0DXltiEWxgufRu3BpPKqoLwp7g8kzSp/lHNcFtf2c8bP+WM+7xRCvbssHU01dJGe7NA9mXnXVC4mNZlocfXisLTzRFNqYd6aqSyx7HkFyCa2syi/eFOGfx3hAWuNRVa+Bwt4i8lj50z5lT3u93sz5o9FoZibpysrHjEv42fzx+a2/5MFff3DXJuxn1Zza4Cbj8aEVbreSx1MrZiWJsuuPM8kiT/pbfPwgS+6YkNEofVC2bJrkqpmHxZ7WeFRVBgp7iyjDL28TUtN+8/3Ozkz6ye+vTW6o3FxJ9glWsz0TCa5rF7vo86S5SVyq3cUlYtpMy/tlqec0i92fR1G+Zd2kwdV5NDJpz1tVr4HC3iLmne8iKVqlDEFIa6iyhlOmbVG/eZqbYjAonvDLdq+ieeld6jwtKmY0qkb0koSsDhdN9J7hYHX8fZGypDUetNhJKvOIpKn6ftGeQZIIJ90/7xZ180wWLLePGaQ1AK7VYBK0vJa8qwgnNR5J9y/6VUqK1Tc1cFndS3lx6cEUEXdbg0UfO0ll3vkuit4v7t4x+exNW7fbtd4/7xad9Tp93b1EgbPFnxdxP/kN1k2Jed1dtiQRThsfsP2tSOcvqa5sYZ221MRlW/dVRgmlMTEijqbGcIpAYW8Z8853kfd+JlHO4lqJ+vi73e7Ua9Yt3hhN9wTMicBco16y1uV0nUzy2MQt7F5vdkUnmwibypZmodoGVYuEQCZFECX1IOL3rMLKdXN/Fcu2aqIqY4zCTmqhiBvF8zzjYiOrq6vO10hqjGYbGLe1S4s2qkmuLd/anayK5Xm3XElrkCSIaTH0SUJWpnimuXaSJmLFewlVuIlcLHbb+ghFqMp9SmEntZB34NNfACQpdt1tSyLPj60MyyutzPE6S4qnD0U4TQST/l6muyMt+iNJ2OPVXkUkSfqAuXlFs6JUFfBAYSe1YBPP+Bc9mpQstILLiIaJE4/nj8fHp4l0EcvLHjrq1vPwr2EW4TQRLOLWyCL8ecP+gNnrVjWwO2kIj9Sf3HZRy8x8aYIW+4LTlNl2TcFm4Y5Goxl3hp9ca2KldzqdQqIeHyg19QBs68jayGt5FR0ATru+iwjm+W5mbRDyhv2ZMjlXHX8+zyAE+tgXmCbNtmsSLj7p9NmlyQIen0nb6/WmJl0liaotP7yJvJZXkbEGwB8wzrM+a1VRJEmPm9SAjEbm69nCHas2lOYZhFDFvSjsc6BJs+3KZB69kLzCN/E/uy4kYt5cFyXJa3mVOcnKdr8q/k/lL/zRzt9IXVDY50C6n3PxluSaVy8kTfhMFnc8W2Pea0ct7qKrUdnI2nCFZU4K7ZzHd6hsIZ5HQq5lgsI+B5IjDxZzEd0yftguQpgkfINgUlM0jn2UYaqii6hO8seYjy06yGV2B8Xz29w0JeYu5a76O1R2w06LvVwo7HMgORytGsGomqIWlmuDZvOx93o940zVrDNf0wYuw/9DknVftLcVbeBsKZLjGSmz9DaqokwXD8ehyoXCPifs4WjzTdxVFkUtrCwNWjxyJXS1lNEoRq9jixNPKm/SOXlIWtQk7d6L9h2Kw8ix8qCw18yiWuzZw93iC3IUF6OyG8Wsg6O2+6+urua24pOWIbR9T+qy2ElzobDXzKL62FXdLawsophFjObdKLo2TvEty//TTxvgZrFnGdglywWFvQEsYlRMFqpyY2xvbyfGqVdNFnH3PM/pf7y9rdrrXY6J+uxi4PHnLPIdKpJdlhOlAAAIoUlEQVTIrc3f20WmdGEHcDeAiwAeiew7AeA+AI8Hr29yudmyCHvbyTLwmEUsTAOrKysrcxOY7e1tZ2HPYsVvb+uVNK62VZ+yTJ5Ke4Y81n4VvQQ2FOVRhbC/D8C7YsJ+B4Dbgve3Abjd5WYU9nbg6jLJKhZNGJ8okpAsqZxpETtlDYzmrcOy657upHKpxBUDYBgT9scAXBO8vwbAYy7XobC3A9cfbVaxaEJEUZFcLy55ZGwTkfJE/pgs4bx1WHbdN6GRbhPzEvbnY39/zuU6FPb24NLNzioWLmIwj+59/B6j0aiwxR69dpHZtWmNalMs9iY00m2iccIOYAPALoDdtbW1OVQBaQpZxSJNtOrs3qcNrGadSGVy+bhcI61Om+Jjp8VeLvMSdrpiSCp5xCJPkq95iEVSiGeenkPeZ3GxhJsQFUMfe7nMS9h/E9ODp3e4XIfCvnyUKRZVdu9dUw7X/SxJDULTolCK5P0h05Qu7AC+BOAHAC4DeALAxwF4AO6HH+54P4ATLjejsJMilGGxm8SvDusy77PYyuqSZ2eewk+LvVwqsdjL2ijspAhFxCLJr20Lcczj4nEVzyKuHdM9qvK956Uut1nTei1lQWEnrSbPDzdvCGNWF09W8XRNWOZCmmunCqGtIuyyCG3uJVDYCYmRFtFi27KKXp2hhmnXqCLBWhVhl0VocyQOhZ2QGGnpcD3PK8XSm+fkIFOs/TyFtmmuH9V2x85T2AmJkWSxu6ylWvQ+ZVvsSQOotmcoW2jLDruss/4XAQo7ITFsPvbV1VX1PK+0gbZ5TQ4qElGTRzzzDNZmuWdZjQ597BR2smRkdV2UdZ+0gdNoeVxFd54uh6LhlVXkFEorb9VRMXVE3lDYCXGgzm57UctynmVPuleawLmWc5F843X1CijshDhQp5gUFeZ5iksRX7prHS+Sb7yuslLYCXGgTjEpo1GZlzugSPRLFj/8ovjG6zIIKOyEONDETJFNtFCLxKtnqeNFmTFKi53CThpOXWIyj0al7EyNtmulWbCLItiu0MdOYSfESpWCN0/xWaTeR1ksQlSM+OfMh9OnT+vu7u7c7kfIMjIcDnFwcDCzfzAYYH9/v9R77ezsYGNjA4eHh1f29ft9bG1tYX19vdR7LTMicl5VT7se36myMIQsKzs7OxgOh+h0OhgOh9jZ2ZnbvS9cuJBpfxHW19extbWFwWAAEcFgMLgi6nXWwdKTxbwvutEVQ5aBuqM8muAeqbsO2gYyumJosRNSMufOnZtyTQDA4eEhzp07N5f7b25uot/vT+3r9/vY3Nycy/2B+utg2aGwE1Iy83SFmEhyj8yLuutg2TlWdwEIaRtra2vGwcu1tbW5lWF9fb3Wwcsm1MEyQ4udkJJpgiukblgH9UJhJ6RkmuAKqRvWQb0wjp0QQhoO49gJIWTJobATQkjLoLATQkjLoLATQkjLoLATQkjLmGtUjIg8BWB21kI9nATwdN2FcIRlLZ9FKSfAslbBopQT8Mu6qqqnXE+Yq7A3CRHZzRI+VCcsa/ksSjkBlrUKFqWcQL6y0hVDCCEtg8JOCCEtY5mFfavuAmSAZS2fRSknwLJWwaKUE8hR1qX1sRNCSFtZZoudEEJayVIKu4jsi8h3ReRhEWlUVjIRuVtELorII5F9J0TkPhF5PHh9U51lDMpkKue/EpH/G9TrwyLywTrLGCIiPyYi3xSRR0XkeyJya7C/UfWaUM7G1auIvE5Evi0ifxqU9TeC/W8RkW8FdfoVEek1uKxfFJG9SL1eV3dZAUBEuiLykIj8cfA5c50upbAH/B1Vva6BIU9fBPCB2L7bANyvqm8DcH/wuW6+iNlyAsBngnq9TlX/85zLZOM1AP9UVd8O4D0AfllE3oHm1autnEDz6vVVADeq6jsBXAfgAyLyHgC3wy/r2wA8B+DjNZYxxFZWAPjnkXp9uL4iTnErgEcjnzPX6TILeyNR1f8B4NnY7r8P4N7g/b0APjzXQhmwlLORqOoPVPV/Be9fhP+jeTMaVq8J5WwcwRrLl4KPK8GmAG4E8LVgf+11CiSWtXGIyLUAfhrA7wafBTnqdFmFXQH8VxE5LyIbdRfGgb+sqj8A/B8/gL9Uc3mS+BUR+U7gqqndZRRHRIYA/iaAb6HB9RorJ9DAeg1cBg8DuAjgPgB/DuB5VX0tOOQJNKRhipdVVcN63Qzq9TMiclWNRQy5E8CvARgHnz3kqNNlFfYbVPVdAH4Kfnf3fXUXqCX8ewB/FX539wcA/m29xZlGRI4D+AMA/0RVX6i7PDYM5WxkvarqkapeB+BaAO8G8HbTYfMtlZl4WUXkbwD4FIC/DuBvATgB4JM1FhEi8iEAF1X1fHS34dDUOl1KYVfVJ4PXiwD+EP6Xssn8UESuAYDg9WLN5TGiqj8MfkBjAL+DBtWriKzAF8sdVf2Pwe7G1aupnE2uVwBQ1ecBPAB/XOBqETkW/OlaAE/WVS4TkbJ+IHB9qaq+CuAe1F+vNwD4GRHZB/Bl+C6YO5GjTpdO2EVkVUTeEL4H8JMAHkk+q3b+CMDNwfubAXy9xrJYCUUy4GfRkHoN/JRfAPCoqn468qdG1autnE2sVxE5JSJXB+9fD+An4I8JfBPAR4LDaq9TwFrW/x1p1AW+37rWelXVT6nqtao6BPALAL6hquvIUadLN0FJRN4K30oHgGMAfk9VG7N0uoh8CcD74Wd0+yGAfwngPwH4KoA1ABcA/Jyq1jpwaSnn++G7CxTAPoB/FPqw60RE3gvgQQDfxcR3+evw/deNqdeEcp5Fw+pVRH4c/kBeF76B+FVV/dfB7+vL8F0bDwH4aGAR10ZCWb8B4BR8d8fDAP5xZJC1VkTk/QD+map+KE+dLp2wE0JI21k6VwwhhLQdCjshhLQMCjshhLQMCjshhLQMCjshhLQMCjshhLQMCjshhLQMCjshhLSM/w8CXT+qTRzkjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "poly = sklearn.preprocessing.PolynomialFeatures(degree=4, include_bias=False)\n",
    "poly_feat = poly.fit_transform(df[[\"RM\", \"DIS\", \"LSTAT\"]])\n",
    "\n",
    "# create and fit regularized model\n",
    "regularized_regr = sklearn.linear_model.Ridge(normalize=True, alpha=0.2)\n",
    "regularized_regr.fit(poly_feat[:300], df.loc[:299, \"MEDV\"])\n",
    "\n",
    "# predict model\n",
    "regularized_pred = regularized_regr.predict(poly_feat[300:])\n",
    "print(\"mean absolute error\", sklearn.metrics.mean_absolute_error(regularized_pred, df.loc[300:, \"MEDV\"]))\n",
    "print(\"mean squared error\", sklearn.metrics.mean_squared_error(regularized_pred, df.loc[300:, \"MEDV\"]))\n",
    "print(\"r^2 score\", sklearn.metrics.r2_score(df.loc[300:, \"MEDV\"], regularized_pred)) # y_true, y_pred\n",
    "\n",
    "# plot predictions w.r.t. LSTAT\n",
    "plt.scatter(df.loc[300:, [\"LSTAT\"]], df.loc[300:, \"MEDV\"], color=\"black\") # true labels\n",
    "plt.scatter(df.loc[300:, [\"LSTAT\"]], regularized_pred, color=\"blue\") # predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there are other hyperparameters that can be tuned, including the degree of the polynomial. A remaining exercise is to tune additional hyperparameters, including degree."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Let's explore a new dataset that looks at handwritten digits, and load the data into matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optical Recognition of Handwritten Digits Data Set\n",
      "===================================================\n",
      "\n",
      "Notes\n",
      "-----\n",
      "Data Set Characteristics:\n",
      "    :Number of Instances: 5620\n",
      "    :Number of Attributes: 64\n",
      "    :Attribute Information: 8x8 image of integer pixels in the range 0..16.\n",
      "    :Missing Attribute Values: None\n",
      "    :Creator: E. Alpaydin (alpaydin '@' boun.edu.tr)\n",
      "    :Date: July; 1998\n",
      "\n",
      "This is a copy of the test set of the UCI ML hand-written digits datasets\n",
      "http://archive.ics.uci.edu/ml/datasets/Optical+Recognition+of+Handwritten+Digits\n",
      "\n",
      "The data set contains images of hand-written digits: 10 classes where\n",
      "each class refers to a digit.\n",
      "\n",
      "Preprocessing programs made available by NIST were used to extract\n",
      "normalized bitmaps of handwritten digits from a preprinted form. From a\n",
      "total of 43 people, 30 contributed to the training set and different 13\n",
      "to the test set. 32x32 bitmaps are divided into nonoverlapping blocks of\n",
      "4x4 and the number of on pixels are counted in each block. This generates\n",
      "an input matrix of 8x8 where each element is an integer in the range\n",
      "0..16. This reduces dimensionality and gives invariance to small\n",
      "distortions.\n",
      "\n",
      "For info on NIST preprocessing routines, see M. D. Garris, J. L. Blue, G.\n",
      "T. Candela, D. L. Dimmick, J. Geist, P. J. Grother, S. A. Janet, and C.\n",
      "L. Wilson, NIST Form-Based Handprint Recognition System, NISTIR 5469,\n",
      "1994.\n",
      "\n",
      "References\n",
      "----------\n",
      "  - C. Kaynak (1995) Methods of Combining Multiple Classifiers and Their\n",
      "    Applications to Handwritten Digit Recognition, MSc Thesis, Institute of\n",
      "    Graduate Studies in Science and Engineering, Bogazici University.\n",
      "  - E. Alpaydin, C. Kaynak (1998) Cascading Classifiers, Kybernetika.\n",
      "  - Ken Tang and Ponnuthurai N. Suganthan and Xi Yao and A. Kai Qin.\n",
      "    Linear dimensionalityreduction using relevance weighted LDA. School of\n",
      "    Electrical and Electronic Engineering Nanyang Technological University.\n",
      "    2005.\n",
      "  - Claudio Gentile. A New Approximate Maximal Margin Classification\n",
      "    Algorithm. NIPS. 2000.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "digits = datasets.load_digits()\n",
    "print(digits[\"DESCR\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = digits.data, digits.target%2 == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though we haven't formally introduced a classification model in the course yet, we'll use a model called support vector machines. The input will be an 8x8 image of integer pixels (between 0 and 16), and the output will be whether the digit is even (binary classification). We'll compute various classification metrics over our model. We'll use the first 1200 instances for training, and the last 597 for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = X[:1200], y[:1200]\n",
    "X_test, y_test = X[1200:], y[1200:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,\n",
       "  decision_function_shape='ovr', degree=3, gamma='auto', kernel='rbf',\n",
       "  max_iter=-1, probability=True, random_state=None, shrinking=True,\n",
       "  tol=0.001, verbose=False)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "svm_model = svm.SVC(probability=True)\n",
    "svm_model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy 0.7202680067001676\n",
      "misclassification error 0.27973199329983245\n"
     ]
    }
   ],
   "source": [
    "y_pred = svm_model.predict(X_test)\n",
    "print(\"accuracy\", accuracy_score(y_test, y_pred))\n",
    "print(\"misclassification error\", 1-accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the testing set, we can calculate the accuracy and misclassification error. We can also compute other important metrics, including precision, recall, and the F1-score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "precision score 1.0\n",
      "recall score 0.4358108108108108\n",
      "F1-score 0.6070588235294118\n"
     ]
    }
   ],
   "source": [
    "print(\"precision score\", precision_score(y_test, y_pred))\n",
    "print(\"recall score\", recall_score(y_test, y_pred))\n",
    "print(\"F1-score\", f1_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also of importance are confusion matrices, which provide a concise presentation of the predictive power of a model, and ROC curves, both of which are easy to produce using scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "confusion matrix\n",
      "TP 301 \t FP 0\n",
      "FN 167 \t TN 129\n"
     ]
    }
   ],
   "source": [
    "conf_matrix = confusion_matrix(y_test, y_pred)\n",
    "print(\"confusion matrix\")\n",
    "print(\"TP\", conf_matrix[0][0], \"\\t\", \"FP\", conf_matrix[0][1])\n",
    "print(\"FN\", conf_matrix[1][0], \"\\t\", \"TN\", conf_matrix[1][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_roc(fpr, tpr, roc_auc):\n",
    "    plt.figure()\n",
    "    lw = 2\n",
    "    plt.plot(fpr, tpr, color='darkorange', lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "    plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\n",
    "    plt.xlim([0.0, 1.0])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.xlabel('False Positive Rate')\n",
    "    plt.ylabel('True Positive Rate')\n",
    "    plt.title('ROC curve')\n",
    "    plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcjXX7wPHPNYtZ7GOQLI0tS9ZIpMeaJYSoB0mbEqLQT3hoQz2SiCxDJWl5VEpEFFJKFLIkhKSM7MvYZsYs1++Pc4wxZs6cYc6cOTPX+/Wal3Pf516uuc2ca77f731fX1FVjDHGmPT4eTsAY4wxOZslCmOMMS5ZojDGGOOSJQpjjDEuWaIwxhjjkiUKY4wxLlmiMMYY45IlCpPriMg+EYkRkbMickhE5ohIgVTb3CYi34jIGRGJFpEvRKR6qm0KicjrIvK381h7nMvh2fsdGeNdlihMbnWXqhYA6gB1gREX3xCRRsDXwELgeqA8sAVYIyIVnNvkA1YCNwFtgULAbcBxoIGnghaRAE8d25irZYnC5Gqqegj4CkfCuGg8MFdVJ6vqGVU9oaqjgHXAC85tHgDKAXer6nZVTVLVI6o6RlW/TOtcInKTiCwXkRMiclhE/uNcP0dExqbYrpmIRKVY3iciw0RkK3BOREaJyPxUx54sIlOcrwuLyNsiclBEDojIWBHxv8ZLZUy6LFGYXE1EygB3Anucy6E4WgafpLH5x0Ar5+s7gGWqetbN8xQEVgDLcLRSKuFokbirB9AeKAK8B7QTkULOY/sD/wY+dG77LpDgPEddoDXwaCbOZUymWKIwudXnInIG2A8cAZ53rg/D8XN/MI19DgIXxx+KpbNNejoAh1T1NVWNdbZUfsrE/lNUdb+qxqjqX8AvQGfney2A86q6TkRK4kh8g1T1nKoeASYB3TNxLmMyxRKFya06q2pBoBlQlUsJ4CSQBJRKY59SwDHn6+PpbJOessAfVxWpw/5Uyx/iaGUA3Mel1sQNQCBwUEROicgpYCZQ4hrObYxLlihMrqaq3wFzgAnO5XPAWuDeNDb/N5e6i1YAbUQkv5un2g9UTOe9c0BoiuXr0go11fInQDNn19ndXEoU+4E4IFxVizi/CqnqTW7GaUymWaIwecHrQCsRuTigPRx4UESeFJGCIlLUOdjcCHjRuc17OD6UPxWRqiLiJyLFROQ/ItIujXMsBq4TkUEiEuQ87q3O9zbjGHMIE5HrgEEZBayqR4FvgXeAP1V1h3P9QRx3bL3mvH3XT0QqikjTq7guxrjFEoXJ9ZwfunOBZ53LPwBtgC44xiH+wjEofLuq7nZuE4djQHsnsBw4DfyMowvrirEHVT2DYyD8LuAQsBto7nz7PRy33+7D8SH/kZuhf+iM4cNU6x8A8gHbcXSlzSdz3WTGZIrYxEXGGGNcsRaFMcYYlyxRGGOMcckShTHGGJcsURhjjHHJ5wqQhYeHa0REhLfDMMYYn7Jx48Zjqlr8avb1uUQRERHBhg0bvB2GMcb4FBH562r3ta4nY4wxLlmiMMYY45IlCmOMMS5ZojDGGOOSJQpjjDEuWaIwxhjjkscShYjMFpEjIrItnfdFRKaIyB4R2SoiN3sqFmOMMVfPk89RzAGm4ijvnJY7gcrOr1uBGc5/TW6VlAAHfoCEGG9HYkyecuFC0jXt77FEoaqrRSTCxSadgLnqqHO+TkSKiEgp58Qs5lrEnYb14yH2hLcjudyWGd6OwJg8Z+gXrdj0z7VNV+LNJ7NLc/k8wVHOdVckChHpA/QBKFeuXLYEd4WkRDj+G6gbmXn3Z7B1Jvjl83xcaTkb5Z3zZkZEW29HYEyeUKNmOFPWRFzTMbyZKCSNdWnOoqSqs4BZAPXr1/fMTEvnj8Fvc9LvFvnxOY+c1qPCa0Ltvt6O4nLBYVC5C/h7KYkak8tt336UX345yP331wLggS5K06ejKV9+9FUf05uJIgoom2K5DPCPx88aHwNf9rzyr+5D693bX/whvEbG2/kFQqtICCmR+RizQkAQhHrp3MaYbHf+fDxjx67m1Vd/xN9faNiwDJUqhSEiREQUuaZjezNRLAIGiMg8HIPY0Vk2PnFyD0R9d+X6c4dgzSjX+xavDRXvSvu9QuWhxsMgaTWGjDHGO5Yu3c0TT3zJn3+eAqB373oUKxaSZcf3WKIQkf8BzYBwEYkCngcCAVQ1EvgSaAfsAc4DD2fqBMe2wbqxkBB75Xt/LMx4/4i2cNuLl68LKgJhN2YqDGOM8ZYDB04zaNBXzJ+/HYBatUoSGdmeRo3KZrBn5njyrqceGbyvwBNXfYKts+D3j1xvU6U7BIZevk78oWZvKGV34hpjfNsTT3zJwoW/ExoayOjRzXjqqYYEBGT943E+Nx9FsqR4x7+1+kDEnVe+H34TFK2cvTEZY4yHJSQkJSeDV165g8BAf157rTXlyhX22Dl9M1Hs+AC2RDpeF68DlTt7Nx5jjPGw6OhYRo36hl27TrBsWU9EhCpVwvnkk3s9fm7fTBQ75116HV7Te3EYY4yHqSqffLKdQYOWcfDgWfz9hc2bD1G37rU9RJcZvpkoLrpzLpS53dtRGGOMR/zxxwkGDFjKsmV7AGjUqAyRkR2oVatktsbh24ki6NruDTbGmJxqwoQfefbZVcTGJlCkSDCvvHIHjz56M35+2X97vm8nCmOMyaXOn48nNjaBXr1qMWFCa0qUyO+1WCxRGGNMDnD06Dl+//04t9/uqGc3bFhjmjWLoEmTG7wcmU1cZIwxXpWUpLz11i9UqTKVLl0+4sQJR725oKCAHJEkwFoUxhjjNdu2HaFv38WsWeMopN2qVQXOn48nLCzrym9kBUsUxhiTzc6du8Do0d8xceI6EhKSKFkyP6+/3pZu3W5CcmAtOUsUxhiTze655xOWLduDCPTvX5+XXmpJkSLB3g4rXb6XKBIvwN7F3o7CGGOu2rBhjTl8+CwzZrTn1lvLeDucDPleoji979Jr/5ybgY0xBhy1md544yf27TvF5MmOunTNmkWwYUMfrzwTcTV8L1EkJTj+rdgRyjb1bizGGOPCzz8f4PHHF7N58yEA+vSpx003OSYU85UkAb58e+xtL9p0msaYHOnUqVj6919Cw4ZvsXnzIW64oTBffNEjOUn4Gt9rURhjTA42b942Bg1axuHD5wgI8OPppxvx7LNNyJ/fd/+wtURhjDFZ6Ouv/+Dw4XM0blyWGTPaU7Nm9hbw8wRLFMYYcw3i4hI4cOAMFSoUBWD8+Fb861/lePDBOj41DuGK745RGGOMl33zzZ/UqhVJ+/YfcuFCIgDh4aE8/HDdXJMkwBKFMcZk2uHDZ+nVawEtW85l167jAERFnfZyVJ5jXU/GGOOmpCTlzTc3Mnz4Sk6diiU4OIBRo/7F0KGNyZfP39vheYwlCmOMcdPdd3/EokW/A9CmTUWmTWtHxYphXo7K86zryRhj3NSlS1Wuu64AH310D0uX9swTSQKsRWGMMelatOh3oqJO07//LQA88EBtunSpRsGCQV6OLHtZojDGmFT+/juaJ59cysKFvxMU5E/btpWoUKEoIpLnkgRYojDGmGTx8YlMmfITzz//LefOxVOwYD7Gjm3BDTcU9nZoXmWJwhhjgHXronj88cVs3XoYgHvvrc6kSW0oXbqQlyPzPksUxhgDPPvsKrZuPUz58kWYOrUd7dpV9nZIOYYlCmNMnqSqnDlzgUKFHGMOU6feydy5Wxg5sgmhoYFeji5nsdtjjTF5zu+/H+OOO96jS5ePUFUAqlQJ56WXWlqSSIO1KIwxeUZsbAL//e/3jBu3hgsXEilWLIR9+05RvnxRb4eWo1miMMbkCcuX/0H//l+yZ88JAB55pA7jx7eiWLFQL0eW83m060lE2orI7yKyR0SGp/F+ORFZJSKbRGSriLTzZDzGmLxHVXnkkYW0bv0+e/acoHr14qxe/RBvv93JkoSbPNaiEBF/YBrQCogC1ovIIlXdnmKzUcDHqjpDRKoDXwIRnorJGJP3iAgREUUICQngueeaMmRIo1xdwM8TPNn11ADYo6p7AURkHtAJSJkoFLh4k3Jh4B8PxmOMySM2bz7EwYNnuPNOxy2uw4Y1plevWjYWcZU82fVUGtifYjnKuS6lF4D7RSQKR2tiYFoHEpE+IrJBRDYkJCR4IlZjTC5w5kwcQ4Z8Rb16s3jwwc85cSIGgKCgAEsS18CTiSKt6Z001XIPYI6qlgHaAe+JyBUxqeosVa2vqvUDAmz83RhzOVVlwYIdVK8+nUmT1gFw3301CQy0JwCygic/daOAsimWy3Bl11JvoC2Aqq4VkWAgHDjiwbiMMbnIX3+dYsCApSxevAuA+vWvZ+bMDtx8cykvR5Z7eDLdrgcqi0h5EckHdAcWpdrmb6AlgIhUA4KBox6MyRiTi6gqXbt+zOLFuyhUKIipU+9k3breliSymMdaFKqaICIDgK8Af2C2qv4mIqOBDaq6CHgaeFNEBuPolnpILz4maYwx6UhKUvz8BBFhwoTWREZuYNKkNpQqVdDboeVK4mufy/XLh+qGATHQaxOUqOPtcIwx2ej48fMMH74CgDff7OjlaHyLiGxU1fpXs6+N9BhjcjxV5d13N1O16jTeemsTc+duJSrqtLfDyjPsFiJjTI62Y8dR+vVbwnff/QVAs2YRzJjRnjJlbJ6I7GKJwhiTI6kqzz23ildeWUN8fBLh4aG89lprevWqhUhad98bT7FEYYzJkUSEAwfOEB+fxGOP3cy4cXcQFhbi7bDyJEsUxpgc459/znDs2Hlq1SoJwPjxrejduy6NG5fzcmR5mw1mG2O8LjExialTf6ZatWl07z6fCxcSAQgPD7UkkQNYi8IY41W//HKQxx9fzIYNjsINTZrcwOnTcYSHWwnwnMKtROF8srqcqu7xcDzGmDzi9Ok4nn32G6ZOXU9SklKmTCGmTGlL585VbbA6h8kwUYhIe2AikA8oLyJ1gOdV9W5PB2eMyZ1UlSZN3mHLlsP4+wtDhjTkhReaUbBgkLdDM2lwZ4xiNHArcApAVTcDlTwZlDEmdxMRBg9uSIMGpdmwoQ+vvdbGkkQO5k7XU7yqnkrVFPRe3Y/EOOcLa5oa4ysuXEhk4sS1+PsLQ4c2BuCBB2pz//218Pe3e2pyOncSxQ4R+TfgJyLlgaeAdZ4NywVNgrLNIKyK10Iwxrjv++//om/fJWzffpSgIH8eeKA2JUsWQETw97c/+HyBO6l8AFAPSAI+A2JxJAvv6fAJBAR7NQRjjGvHjp3nkUcW0qTJHLZvP0rlymEsXnwfJUsW8HZoJpPcaVG0UdVhwLCLK0SkC46k4R12R4QxOZaqMmfOZoYOXc7x4zHky+fPiBG3M3z47QQH2x35vsidFsWoNNaNzOpAjDG5x/vv/8rx4zG0aFGerVv78sILzSxJ+LB0/+dEpA2OaUpLi8jEFG8VwtENZYwxAJw/H090dCylShVERJg+vR3r1/9Dz5417ZmIXMBVij8CbMMxJvFbivVngOGeDMoY4zuWLt3NE098SYUKRVm+vBciQpUq4VSpEu7t0EwWSTdRqOomYJOIfKCqsdkYkzHGBxw4cJpBg75i/vztABQsGMTx4zFWeiMXcqfTsLSIvARUB5JvNVLVGz0WlTEmx0pMTGLatPWMGvUNZ85cIH/+QEaPbs6TT95KQIA9E5EbuZMo5gBjgQnAncDD2BiFMXlSUpLStOkc1qzZD0DnzlWZPLkt5coV9nJkxpPcSf+hqvoVgKr+oaqjgOaeDcsYkxP5+QmtW1ekbNlCLFzYnQULulmSyAPcaVHEieO2hT9EpC9wACjh2bCMMTmBqvLxx78REOBH167VARg2rDFDhjSiQIF8Xo7OZBd3EsVgoADwJPASUBh4xJNBGWO8748/TtC//5d8/fUfFC8eSosW5SlaNISgoACCrH5fnpJholDVn5wvzwC9AESkjCeDMsZ4T1xcAq+++iMvvfQ9sbEJFC0azEsvtaBwYSubk1e5TBQicgtQGvhBVY+JyE04Snm0ACxZGJPLfPvtPvr1W8LOnccA6NWrFhMmtKZEifxejsx4U7qD2SLyX+ADoCewTERGAquALYDdGmtMLpOYmET//o4kUaVKMb755gHmzr3bkoRx2aLoBNRW1RgRCQP+cS7/nj2hGWM8LSlJiY1NIDQ0EH9/P2bMaM/q1X/xzDONCQqy2kzGwdVPQqyqxgCo6gkR2WlJwpjc49dfD9O37xKqVi3G2293AqBp0wiaNo3wbmAmx3GVKCqIyMVS4gJEpFhGVbt4NDJjjEecO3eB0aO/Y+LEdSQkJPHnnyc5eTKGokVDvB2ayaFcJYquqZanejIQY4znffHF7wwYsJS//45GBPr3r89LL7WkSBG7o8mkz1VRwJXZGYgxxnMSEpLo1m0+n322A4A6da5j5swONGhQ2suRGV9go1XG5AEBAX4ULhxEgQL5GDOmOQMGNLACfsZtoqqeO7hIW2Ay4A+8parj0tjm38ALgAJbVPU+V8esX1Z0w65jEFLMAxEbk3v89FMUALfe6njk6fjx88TEJFCmTCFvhmW8REQ2qmr9q9nX7RaFiASpalwmtvcHpgGtgChgvYgsUtXtKbapDIwAGqvqSRGxGlLGXKNTp2IZMWIFM2dupGrVcDZv7ku+fP4UK2bzRJirk2HbU0QaiMivwG7ncm0RecONYzcA9qjqXlW9AMzD8WxGSo8B01T1JICqHslU9MaYZKrKhx/+StWqU4mM3Ii/vx8dO1YhMdFmBTDXxp0WxRSgA/A5gKpuERF3yoyXBvanWI4Cbk21zY0AIrIGR/fUC6q6zI1jG2NS2L37OP37f8mKFXsBaNy4LJGRHahRwxrp5tq5kyj8VPWvVBOkJ7qxX1ozqqceEAkAKgPNcNSO+l5EaqjqqcsOJNIH6ANQzypMGXOZ+PhEWrSYS1TUacLCQhg//g4efrgufn5p/Qoak3nuJIr9ItIAUOe4w0Bglxv7RQFlUyyXwVEGJPU261Q1HvhTRH7HkTjWp9xIVWcBs8AxmO3GuY3J9VQVESEw0J+XXmrBqlX7GD/+DooXt9pMJmu5c39cP2AIUA44DDR0rsvIeqCyiJQXkXxAd2BRqm0+xzlbnoiE4+iK2ute6MbkTYcPn6VXrwWMHbs6ed0DD9TmnXc6WZIwHuFOiyJBVbtn9sCqmiAiA4CvcIw/zFbV30RkNLBBVRc532stIttxdGcNVdXjmT2XMXlBUpLy5psbGT58JadOxVKkSDCDBjWkYEGbRch4VobPUYjIH8DvwEfAZ6p6JjsCS489R2Hyoi1bDtG37xLWrXM8G9G2bSWmTWtHhQpFvRyZ8RUefY5CVSuKyG04uo5eFJHNwDxVnXc1JzTGuC8+PpERI1by+uvrSExUSpUqwOTJbbnnnuqkusHEGI9x6xl+Vf1RVZ8EbgZO45jQyBjjYQEBfmzadIikJGXgwAbs2PEE9957kyUJk60ybFGISAEcD8p1B6oBC4HbPByXMXnW339Hk5iYRPnyRRERIiPbEx0dR/3613s7NJNHuTOYvQ34Ahivqt97OB5j8qz4+EQmT/6J55//lkaNyrB8eS9EhMqVbTzOeJc7iaKCqloNAGM8aO3a/fTtu4StWw8DEBYWwvnz8eTPn8/LkRnjIlGIyGuq+jTwqciVD7nZDHfGXLuTJ2MYPnwFs2b9AkD58kWYNq0dd95Z2cuRGXOJqxbFR85/bWY7YzwgLi6BOnVm8vff0QQG+jF06G2MHNmE0NBAb4dmzGVczXD3s/NlNVW9LFk4H6SzGfCMuQZBQQH07l2XlSv/ZMaM9lSvXtzbIRmTJnceuPtFVW9OtW6Tqtb1aGTpsAfujK+KjU3gv//9nipVwrnvvpqAY4pSf3+x212Nx3nkgTsR6YbjltjyIvJZircKAqfS3ssYk5bly/+gf/8v2bPnBCVK5Ofuu6sSEhJo05Ean+BqjOJn4DiOqq/TUqw/A2zyZFDG5BaHDp1lyJCv+N//tgFw003FiYzsQEiIjUMY3+FqjOJP4E9gRfaFY0zukJiYxMyZG/nPf1YSHR1HSEgAzz/flMGDG5Evn7+3wzMmU1x1PX2nqk1F5CSXTzgkgKpqmMejM8ZHJSYqb7zxM9HRcbRrV5mpU++kfHkr4Gd8k6uup4vTnYZnRyDG+LozZ+JITFSKFAkmXz5/3nzzLg4fPkuXLtVssNr4tHRH0lI8jV0W8FfVRKAR8Dhgs6MY46SqfPbZDqpVm8bTT3+VvP7228vRtatVeTW+z51bLj7HMQ1qRWAujsKAH3o0KmN8xL59p+jYcR5du37MgQNn2LbtKLGxCd4Oy5gs5U6iSHLOad0FeF1VBwKlPRuWMTlbfHwir7zyA9WrT2Px4l0UKhTE1Kl38uOPjxAc7E4JNWN8h1tToYrIvUAvoLNznd3bZ/Ks8+fjadjwLX799QgA3bvXYOLE1pQqVdDLkRnjGe4kikeA/jjKjO8VkfLA/zwbljE5V2hoIPXrX8/58/FMn96e1q0rejskYzwqwxIeACISAFRyLu5RVa91wloJD5PdVJW5c7dQsWIYt99eDoDo6Fjy5fO3B+eMz/DonNki8i/gPeAAjmcorhORXqq65mpOaIwv2bHjKP36LeG77/6iWrVwNm/uS758/hQuHOzt0IzJNu50PU0C2qnqdgARqYYjcVxVZjLGF8TExPPSS98zfvwa4uOTKF48lBEjbicw0GozmbzHnUSR72KSAFDVHSJi026ZXGvZsj088cSX7N17EoDHHruZcePuICwsxMuRGeMd7iSKX0RkJo5WBEBPrCigyaXOnr1Ar14LOHbsPDVqlCAysj2NG5fzdljGeJU7iaIv8CTwDI4xitXAG54MypjslJiYRFKSEhjoT4EC+Zg8uS1RUacZPLghgYFWwM8Yl3c9iUhNoCLwm6ruzraoXLC7nkxW2rjxHx5/fDGdOlXh2WebejscYzzmWu56SndkTkT+g6N8R09guYg8cpXxGZPjnD4dx1NPLaVBg7fYuPEg7723lfj4RG+HZUyO5KrrqSdQS1XPiUhx4EtgdvaEZYxnqCrz52/nqaeWcfDgWfz9hSFDGvLii82tm8mYdLhKFHGqeg5AVY+KiN0XaHzamTNxdOs2n6VL9wBw662liYzsQJ0613k5MmNyNleJokKKubIFqJhy7mxV7eLRyIzJYgUK5CMuLpHChYMYN+4O+vSph5+flQA3JiOuEkXXVMtTPRmIMZ6wevVflCpVgMqViyEizJ7dkeDgAEqWLODt0IzxGa7mzF6ZnYEYk5WOHTvPM88s5513NtOyZXmWL++FiHDDDUW8HZoxPscK55tcJSlJmTNnM0OHLufEiRjy5fPnX/8qR2KiEhBg3UzGXA2PDlCLSFsR+V1E9ojIcBfb3SMiKiJWP8pctd9+O0KzZnPo3XsRJ07E0LJleX79tR/PP9+MgAC7F8OYq+V2i0JEglQ1LhPb+wPTgFZAFLBeRBalrBvl3K4gjie/f3L32MakFh0dS8OGb3P27AVKlMjPxImtue++mjZftTFZIMM/s0SkgYj8Cux2LtcWEXdKeDTAMXfFXlW9AMwDOqWx3RhgPBDrftjGOFysLFC4cDDDhjWmb9967Nz5BD171rIkYUwWcac9PgXoABwHUNUtQHM39isN7E+xHEWqubZFpC5QVlUXuzqQiPQRkQ0issGN85o84MCB09xzz8e8//7W5HUjR/6LGTM6ULSoVXk1Jiu5kyj8VPWvVOvcqXWQ1p9zyYWlnA/wTQKezuhAqjpLVetfbZ0Sk3skJCQxefI6qladxqef7uD5578lMTEJwFoQxniIO2MU+0WkAaDOcYeBwC439osCyqZYLgP8k2K5IFAD+Nb5C34dsEhEOqqqtRzMFdavP0Dfvkv45ZeDAHTuXJUpU9ri728D1cZ4kjuJoh+O7qdywGFghXNdRtYDlUWkPI5pVLsD9118U1WjgfCLyyLyLfB/liRMaufOXWDYsBVMn74eVShXrjBvvHEnHTtW8XZoxuQJGSYKVT2C40M+U1Q1QUQGAF8B/sBsVf1NREYDG1R1UaajNXlSQIAfK1bsxc9PGDKkEc8/35T8+W2SRWOyi8v5KABE5E1SjC1cpKp9PBWUKzYfRd7wxx8nKFIkmGLFQgFHt1NwcAA1a5b0cmTG+CaPzEeRwgpgpfNrDVACcPt5CmMyIy4ugbFjV1OjxgyGDVuRvP6WW0pbkjDGS9zpevoo5bKIvAcs91hEJs/69tt99Ou3hJ07jwGOO5wSE5NssNoYL7uaWk/lgRuyOhCTdx05co6hQ5czd+4WAKpUKcaMGe1p3ry8lyMzxoAbiUJETnJpjMIPOAGkW7fJmMw4duw81apN48SJGIKC/Bk58l8880xjgoKsXqUxOYXL30ZxPOBQG8ftrQBJmtHotzGZEB4eSqdOVYiKOs306e2pVCnM2yEZY1JxmShUVUVkgarWy66ATO527twFRo/+jvbtb6RJE0cP5vTp7QkK8rcnq43JodwZJfxZRG72eCQm1/vii9+pXn0648f/SP/+S0hKcjROg4MDLEkYk4Ol26IQkQBVTQBuBx4TkT+AczhqOKmqWvIwbtm/P5qnnlrGggU7Aahb9zpmzuxg81Ub4yNcdT39DNwMdM6mWEwuk5CQxJQpP/Hcc6s4dy6eAgXyMXZsc554ooFNJGSMD3GVKARAVf/IplhMLnP6dBz//e8PnDsXT9eu1Xj99baUKVPI22EZYzLJVaIoLiJD0ntTVSd6IB7j406diiUkJICgoADCwkKYObMDQUH+tG9/o7dDM8ZcJVftf3+gAI5y4Gl9GZNMVfnww1+pUmUq48evSV7fpUs1SxLG+DhXLYqDqjo62yIxPmvXruP077+ElSv/BGD16r9RVbuTyZhcIsMxCmPSExubwCuv/MDLL//AhQuJhIWF8OqrrXjooTqWJIzJRVwlipbZFoXxOYcOnaVJk3fYvfsEAA89VIdXX21FeHiolyMzxmS1dBOFqp7IzkCMbylZMj9lyxYmIMCPGTPa07RphLdDMsZ4iFVeM25JSlLefHMjzZuX58YbiyEifPhhF4oWDSHz4buVAAAYNklEQVRfPn9vh2eM8SB76slkaMuWQzRuPJu+fZfQv/8SLtaFLFmygCUJY/IAa1GYdJ09e4EXXviW119fR2Kicv31Benb96pmUjTG+DBLFCZNn3++k4EDlxIVdRo/P2HgwAaMHduCQoWCvB2aMSabWaIwVzhw4DTdu88nLi6RevVKERnZgfr1r/d2WMYYL7FEYQCIj08kIMAPEaF06UK89FIL8uXzp3//W2zOamPyOPsEMPz4437q1ZvF++9vTV739NO3MXDgrZYkjDGWKPKyEydiePzxL2jceDa//nqE6dM3YDPdGmNSs66nPEhVef/9rTz99NccPXqewEA/nnmmMSNH/stKbxhjrmCJIo85fPgsPXp8yqpV+wBo2vQGZsxoT7Vqxb0bmDEmx7JEkccUKRLMwYNnCQ8PZcKEVjzwQG1rRRhjXLJEkQcsX/4HN99cimLFQgkKCuCTT+6lVKkCFCtmBfyMMRmzwexc7ODBM/To8SmtW7/PsGErktfXqFHCkoQxxm3WosiFEhOTmDlzIyNGrOT06ThCQgKoUqWYTSZkjLkqlihymV9+OUjfvotZv/4fANq3r8zUqe2IiCji5ciMMb7KEkUusm/fKRo0eJPERKV06YJMmXInd99d1VoRxphr4tFEISJtgcmAP/CWqo5L9f4Q4FEgATgKPKKqf3kyptwsIqIIDz9ch4IFg3jxxWYULGgF/Iwx185jg9ki4g9MA+4EqgM9RKR6qs02AfVVtRYwHxjvqXhyo337TnHXXf/ju+/2Ja+bNesuJk5sY0nCGJNlPNmiaADsUdW9ACIyD+gEbL+4gaquSrH9OuB+D8aTa8THJzJx4lpefPE7YmISOHbsPGvX9gawbiZjTJbz5O2xpYH9KZajnOvS0xtYmtYbItJHRDaIyIYsjM8n/fDD39StO5Phw1cSE5NA9+41+Oyzf3s7LGNMLubJFkVaf9qmWXFORO4H6gNN03pfVWcBswDql5U8WbXu5MkYhg5dzttvbwKgYsWiTJ/entatK3o5MmNMbufJRBEFlE2xXAb4J/VGInIHMBJoqqpxHozHpyUlKQsX/k5goB/Dh9/OiBG3ExIS6O2wjDF5gCcTxXqgsoiUBw4A3YH7Um4gInWBmUBbVT3iwVh80s6dxyhfvghBQQEUKxbKBx90oVy5wlStGu7t0IwxeYjHxihUNQEYAHwF7AA+VtXfRGS0iHR0bvYqUAD4REQ2i8giT8XjS86fj2fkyJXUqjWD8ePXJK9v3bqiJQljTLbz6HMUqvol8GWqdc+leH2HJ8/vi5Yt20P//kv4889TABw7dt7LERlj8jp7MjuH+OefMwwatIxPPnHcPVyzZgkiIztw221lM9jTGGM8yxJFDrBr13Hq15/FmTMXCA0N5IUXmjJoUEMCA/29HZoxxvhgosh/HQSEeDuKLFW5chi33FKa/PkDeeONO7nhBivgZ4zJOXwvURQoDYG+PZfC6dNxPPfcKvr3v4UbbyyGiLBoUXfy58/n7dCMMeYKvpcofJiqMn/+dp56ahkHD55l585jLFvmqFpiScIYk1NZosgme/eeZMCAL1m6dA8ADRuW4ZVX7KYvY0zOZ4nCwy5cSGTChB8ZM2Y1sbEJFCkSzLhxLXnssXr4+VkBP2NMzmeJwsP2749m9OjviItLpGfPmrz2WmtKlizg7bCMMcZtlig84OTJGIoUCUZEqFgxjMmT21KpUhgtW1bwdmjGGJNpniwznuckJSmzZ2+iUqU3eP/9rcnrH3+8viUJY4zPskSRRX777QjNms2hd+9FnDgRkzxobYwxvs66nq7R+fPxjBnzHRMmrCUhIYkSJfIzaVIbevSo4e3QjDEmS1iiuAa7dh2nTZv32bfvFCLQt289Xn65JUWL5q4nx40xeZslimtwww2FCQ4OoHbtkkRGdqBhwzLeDsnkIPHx8URFRREbG+vtUEweEhwcTJkyZQgMzLqJzSxRZEJCQhKRkRvo0aMGxYqFEhQUwLJlPSlduhABATbcYy4XFRVFwYIFiYiIQMSemTGep6ocP36cqKgoypcvn2XHtU83N/388wEaNHiTgQOXMmzYiuT1N9xQxJKESVNsbCzFihWzJGGyjYhQrFixLG/FWosiA9HRsYwc+Q3Tp69HFcqVK0ynTlW8HZbxEZYkTHbzxM+cJYp0qCofffQbgwd/xaFDZwkI8GPIkIY891xTK+BnjMlTrM8kHVu2HKZHj085dOgst91Wll9+6cMrr7SyJGF8ir+/P3Xq1KFGjRrcddddnDp1Kvm93377jRYtWnDjjTdSuXJlxowZg6omv7906VLq169PtWrVqFq1Kv/3f//njW/BpU2bNvHoo496O4x0HT9+nObNm1OgQAEGDBiQ7nYnTpygVatWVK5cmVatWnHy5EnA8Qfrk08+SaVKlahVqxa//PILAEePHqVt27bZ8j2AJYrLJCYmJb+uU+c6Bg9uyJtv3sX33z9MzZolvRiZMVcnJCSEzZs3s23bNsLCwpg2bRoAMTExdOzYkeHDh7Nr1y62bNnCjz/+yPTp0wHYtm0bAwYM4P3332fHjh1s27aNChWytrpAQkLCNR/j5ZdfZuDAgdl6zswIDg5mzJgxTJgwweV248aNo2XLluzevZuWLVsybtw4wJGsd+/eze7du5k1axb9+vUDoHjx4pQqVYo1a9Z4/HsA63pKtmrVn/Tv/yUzZ3agSZMbAJg4sY2XozK5xmseGqt4WjPexqlRo0Zs3eooLfPhhx/SuHFjWrduDUBoaChTp06lWbNmPPHEE4wfP56RI0dStWpVAAICAujfv/8Vxzx79iwDBw5kw4YNiAjPP/88Xbt2pUCBApw9exaA+fPns3jxYubMmcNDDz1EWFgYmzZtok6dOixYsIDNmzdTpIhjVsdKlSqxZs0a/Pz86Nu3L3///TcAr7/+Oo0bN77s3GfOnGHr1q3Url0bgJ9//plBgwYRExNDSEgI77zzDlWqVGHOnDksWbKE2NhYzp07xzfffMOrr77Kxx9/TFxcHHfffTcvvvgiAJ07d2b//v3Exsby1FNP0adPH7evb1ry58/P7bffzp49ris1LFy4kG+//RaABx98kGbNmvHKK6+wcOFCHnjgAUSEhg0bcurUKQ4ePEipUqXo3LkzH3zwwRXXxRPyfKI4cuQcQ4cuZ+7cLQBMnLg2OVEYk1skJiaycuVKevfuDTi6nerVq3fZNhUrVuTs2bOcPn2abdu28fTTT2d43DFjxlC4cGF+/fVXgOQuE1d27drFihUr8Pf3JykpiQULFvDwww/z008/ERERQcmSJbnvvvsYPHgwt99+O3///Tdt2rRhx44dlx1nw4YN1KhxqQJC1apVWb16NQEBAaxYsYL//Oc/fPrppwCsXbuWrVu3EhYWxtdff83u3bv5+eefUVU6duzI6tWradKkCbNnzyYsLIyYmBhuueUWunbtSrFixS477+DBg1m1atUV31f37t0ZPnx4ht9/Wg4fPkypUqUAKFWqFEeOHAHgwIEDlC1bNnm7MmXKcODAAUqVKkX9+vUZNWrUVZ0vs/JsokhKUt5++xeGDVvByZOxBAX5M2pUE4YOvc3boZncKBN/+WelmJgY6tSpw759+6hXrx6tWrUCHH3f6d0dk5m7ZlasWMG8efOSl4sWLZrhPvfeey/+/v4AdOvWjdGjR/Pwww8zb948unXrlnzc7du3J+9z+vRpzpw5Q8GCBZPXHTx4kOLFiycvR0dH8+CDD7J7925EhPj4+OT3WrVqRVhYGABff/01X3/9NXXr1gUcraLdu3fTpEkTpkyZwoIFCwDYv38/u3fvviJRTJo0yb2LkwVSjhlddPH/p0SJEvzzzz/ZEkeeTBR//nmS++9fwI8/7gegdeuKTJvWjkqVwrwcmTFZ6+IYRXR0NB06dGDatGk8+eST3HTTTaxevfqybffu3UuBAgUoWLAgN910Exs3bkzu1klPegkn5brU9/Tnz58/+XWjRo3Ys2cPR48e5fPPP0/+CzkpKYm1a9cSEpJ+OZyQkJDLjv3ss8/SvHlzFixYwL59+2jWrFma51RVRowYweOPP37Z8b799ltWrFjB2rVrCQ0NpVmzZmk+j+CJFkXJkiWTu5QOHjxIiRIlAEcLYv/+/cnbRUVFcf311wOO6+rq+mSlPDmYXahQELt2Hee66wowb15Xli3raUnC5GqFCxdmypQpTJgwgfj4eHr27MkPP/zAihWOh0djYmJ48skneeaZZwAYOnQoL7/8Mrt27QIcH9wTJ0684ritW7dm6tSpycsXu55KlizJjh07kruW0iMi3H333QwZMoRq1aol//We+ribN2++Yt9q1apd1vcfHR1N6dKlAZgzZ06652zTpg2zZ89OHkM5cOAAR44cITo6mqJFixIaGsrOnTtZt25dmvtPmjSJzZs3X/F1tUkCoGPHjrz77rsAvPvuu3Tq1Cl5/dy5c1FV1q1bR+HChZO7qHbt2nVZ15sn5ZlE8dVXe4iLc9zxUKxYKIsWdWfnzifo1q2GPRRl8oS6detSu3Zt5s2bR0hICAsXLmTs2LFUqVKFmjVrcssttyTfwlmrVi1ef/11evToQbVq1ahRowYHDx684pijRo3i5MmT1KhRg9q1ayf/pT1u3Dg6dOhAixYtkj/Y0tOtWzfef//95G4ngClTprBhwwZq1apF9erViYyMvGK/qlWrEh0dzZkzZwB45plnGDFiBI0bNyYxMTHd87Vu3Zr77ruPRo0aUbNmTe655x7OnDlD27ZtSUhIoFatWjz77LM0bNgw44vqhoiICIYMGcKcOXMoU6ZMcpfao48+yoYNGwAYPnw4y5cvp3Llyixfvjw56bRr144KFSpQqVIlHnvsseS70gBWrVpF+/btsyTGjEhafWA5Wf369fXixXXH/v3RPPnkMj7/fCdjxjRn1KgmHozOmEt27NhBtWrVvB1GrjZp0iQKFiyYo5+l8JQmTZqwcOHCNMeF0vrZE5GNqlr/as6Va1sUCQlJTJy4lmrVpvH55zspUCAfYWFW/tuY3KRfv34EBQV5O4xsd/ToUYYMGeLWzQNZIVcOZq9bF0XfvovZsuUwAF27VmPy5LaULl3Iy5EZY7JScHAwvXr18nYY2a548eJ07tw5286X6xLFTz9Fcdttb6MKERFFmDr1Ttq3v9HbYZk8ytVtqMZ4gieGE3JdomjQoDRt2lSibt3rGDWqCaGhWTd5hzGZERwczPHjx63UuMk2F+ejCA4OztLj+nyi2L37OIMHf8XEiW248UbHL+SSJffh52e/mMa7ypQpQ1RUFEePHvV2KCYPuTjDXVby2UQRF5fAuHE/8N///kBcXCLBwQHMn/9vAEsSJkcIDAzM0lnGjPEWj971JCJtReR3EdkjIlc8jSIiQSLykfP9n0Qkwp3jrly5l1q1Innhhe+Ii0vk4YfrEBnZIavDN8YYgwdbFCLiD0wDWgFRwHoRWaSq21Ns1hs4qaqVRKQ78ArQ7cqjXfLnn6e44473AKhWLZzIyA5WxM8YYzzIky2KBsAeVd2rqheAeUCnVNt0At51vp4PtJQMRv1OnowhODiAl19uwebNfS1JGGOMh3nsyWwRuQdoq6qPOpd7Abeq6oAU22xzbhPlXP7Duc2xVMfqA1wsDF8D2OaRoH1POHAsw63yBrsWl9i1uMSuxSVVVLVgxptdyZOD2Wm1DFJnJXe2QVVnAbMARGTD1T6GntvYtbjErsUldi0usWtxiYi4X/soFU92PUUBZVMslwFSF09P3kZEAoDCwAkPxmSMMSaTPJko1gOVRaS8iOQDugOLUm2zCHjQ+foe4Bv1tSqFxhiTy3ms60lVE0RkAPAV4A/MVtXfRGQ0sEFVFwFvA++JyB4cLYnubhx6lqdi9kF2LS6xa3GJXYtL7FpcctXXwufKjBtjjMleubbMuDHGmKxhicIYY4xLOTZReKr8hy9y41oMEZHtIrJVRFaKSK59CjGja5Fiu3tEREUk194a6c61EJF/O382fhORD7M7xuzixu9IORFZJSKbnL8n7bwRp6eJyGwROeJ8Ri2t90VEpjiv01YRudmtA6tqjvvCMfj9B1AByAdsAaqn2qY/EOl83R34yNtxe/FaNAdCna/75eVr4dyuILAaWAfU93bcXvy5qAxsAoo6l0t4O24vXotZQD/n6+rAPm/H7aFr0QS4GdiWzvvtgKU4nmFrCPzkznFzaovCI+U/fFSG10JVV6nqeefiOhzPrORG7vxcAIwBxgOx2RlcNnPnWjwGTFPVkwCqeiSbY8wu7lwLBS5OcVmYK5/pyhVUdTWun0XrBMxVh3VAEREpldFxc2qiKA3sT7Ec5VyX5jaqmgBEA8WyJbrs5c61SKk3jr8YcqMMr4WI1AXKquri7AzMC9z5ubgRuFFE1ojIOhFpm23RZS93rsULwP0iEgV8CQzMntBynMx+ngA5dz6KLCv/kQu4/X2KyP1AfaCpRyPyHpfXQkT8gEnAQ9kVkBe583MRgKP7qRmOVub3IlJDVU95OLbs5s616AHMUdXXRKQRjue3aqhqkufDy1Gu6nMzp7YorPzHJe5cC0TkDmAk0FFV47IptuyW0bUoiKNo5Lcisg9HH+yiXDqg7e7vyEJVjVfVP4HfcSSO3Mada9Eb+BhAVdcCwTgKBuY1bn2epJZTE4WV/7gkw2vh7G6ZiSNJ5NZ+aMjgWqhqtKqGq2qEqkbgGK/pqKpXXQwtB3Pnd+RzHDc6ICLhOLqi9mZrlNnDnWvxN9ASQESq4UgUeXGO2kXAA867nxoC0ap6MKOdcmTXk3qu/IfPcfNavAoUAD5xjuf/raodvRa0h7h5LfIEN6/FV0BrEdkOJAJDVfW496L2DDevxdPAmyIyGEdXy0O58Q9LEfkfjq7GcOd4zPNAIICqRuIYn2kH7AHOAw+7ddxceK2MMcZkoZza9WSMMSaHsERhjDHGJUsUxhhjXLJEYYwxxiVLFMYYY1yyRGFyHBFJFJHNKb4iXGwbkV6lzEye81tn9dEtzpIXVa7iGH1F5AHn64dE5PoU770lItWzOM71IlLHjX0GiUjotZ7b5F2WKExOFKOqdVJ87cum8/ZU1do4ik2+mtmdVTVSVec6Fx8Crk/x3qOquj1LorwU53Tci3MQYInCXDVLFMYnOFsO34vIL86v29LY5iYR+dnZCtkqIpWd6+9PsX6miPhncLrVQCXnvi2dcxj86qz1H+RcP04uzQEywbnuBRH5PxG5B0fNrQ+c5wxxtgTqi0g/ERmfIuaHROSNq4xzLSkKuonIDBHZII65J150rnsSR8JaJSKrnOtai8ha53X8REQKZHAek8dZojA5UUiKbqcFznVHgFaqejPQDZiSxn59gcmqWgfHB3WUs1xDN6Cxc30i0DOD898F/CoiwcAcoJuq1sRRyaCfiIQBdwM3qWotYGzKnVV1PrABx1/+dVQ1JsXb84EuKZa7AR9dZZxtcZTpuGikqtYHagFNRaSWqk7BUcunuao2d5byGAXc4byWG4AhGZzH5HE5soSHyfNinB+WKQUCU5198ok46halthYYKSJlgM9UdbeItATqAeud5U1CcCSdtHwgIjHAPhxlqKsAf6rqLuf77wJPAFNxzHXxlogsAdwuaa6qR0Vkr7POzm7nOdY4j5uZOPPjKFeRcoayf4tIHxy/16VwTNCzNdW+DZ3r1zjPkw/HdTMmXZYojK8YDBwGauNoCV8xKZGqfigiPwHtga9E5FEcZZXfVdURbpyjZ8oCgiKS5vwmztpCDXAUmesODABaZOJ7+Qj4N7ATWKCqKo5PbbfjxDGL2zhgGtBFRMoD/wfcoqonRWQOjsJ3qQmwXFV7ZCJek8dZ15PxFYWBg875A3rh+Gv6MiJSAdjr7G5ZhKMLZiVwj4iUcG4TJu7PKb4TiBCRSs7lXsB3zj79wqr6JY6B4rTuPDqDo+x5Wj4DOuOYI+Ej57pMxamq8Ti6kBo6u60KAeeAaBEpCdyZTizrgMYXvycRCRWRtFpnxiSzRGF8xXTgQRFZh6Pb6Vwa23QDtonIZqAqjikft+P4QP1aRLYCy3F0y2RIVWNxVNf8RER+BZKASBwfuoudx/sOR2sntTlA5MXB7FTHPQlsB25Q1Z+d6zIdp3Ps4zXg/1R1C475sX8DZuPozrpoFrBURFap6lEcd2T9z3medTiulTHpsuqxxhhjXLIWhTHGGJcsURhjjHHJEoUxxhiXLFEYY4xxyRKFMcYYlyxRGGOMcckShTHGGJf+H+ThgwIOMLDrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fpr, tpr, thresholds = roc_curve(y_test, svm_model.predict_proba(X_test)[:, 1])\n",
    "roc_auc = auc(fpr, tpr)\n",
    "plot_roc(fpr, tpr, roc_auc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's all! A remaining exercise is to perform k-fold cross-validation to tune a hyperparameter for support vector machines, if time remains. (One possible such hyperparameter is the kernel used.)"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}