Comparing Machine Learning and Neural Network Performance on Heart Disease Data Set.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification using Heart Disease Data Set\n",
    "\n",
    "This notebook compares the different classification techniques of a logistic regression model, a support vector machine model, and a 3-layer neural network that uses a sigmoid activator in Python on a Heart Disease dataset that has the outcome being classified as having heart disease (indicated by a \"1\") or not (indicates by a \"0\"). \n",
    "\n",
    "Since many models are being compared, a train, validate, test split will be used with a respective split of 60%, 20%, and 20%. Performance will be assessed by comparing accuracy and F1 scores.\n",
    "\n",
    "More information about the dataset is under **The Data Source** section.Now to load the libraries and data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import researchpy as rp\n",
    "import numpy as np\n",
    "\n",
    "from sklearn import svm\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "from sklearn.model_selection import train_test_split # Fast way for splitting data into train, validate, and test sets\n",
    "from sklearn.metrics import accuracy_score, f1_score\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 303 entries, 0 to 302\n",
      "Data columns (total 14 columns):\n",
      "age         303 non-null int64\n",
      "sex         303 non-null int64\n",
      "cp          303 non-null int64\n",
      "trestbps    303 non-null int64\n",
      "chol        303 non-null int64\n",
      "fbs         303 non-null int64\n",
      "restecg     303 non-null int64\n",
      "thalach     303 non-null int64\n",
      "exang       303 non-null int64\n",
      "oldpeak     303 non-null float64\n",
      "slope       303 non-null int64\n",
      "ca          303 non-null int64\n",
      "thal        303 non-null int64\n",
      "target      303 non-null int64\n",
      "dtypes: float64(1), int64(13)\n",
      "memory usage: 33.2 KB\n"
     ]
    }
   ],
   "source": [
    "#df = pd.read_csv(\"C:\\\\Users\\\\CoreySSD\\\\Downloads\\\\heart.csv\")\n",
    "df = pd.read_csv(\"C:\\\\Users\\\\bryantcm\\\\Downloads\\\\heart.csv\")\n",
    "\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Data Source\n",
    "\n",
    "The data is from [Kaggle](https://www.kaggle.com/ronitf/heart-disease-uci) and was uploaded by the user <i>ronit</i>. The data contains attributes, that are a subset from a larger database, that are believed to predict the presence of heart disease. The variables in the data are:\n",
    "\n",
    "* *age* = patients age\n",
    "* *sex* = patients sex (1 = male, 0 = female)\n",
    "* *cp* = chest pain type\n",
    "    * 1 indicates typical angina\n",
    "    * 2 indicates atypical angina\n",
    "    * 3 indicates non-anginal pain\n",
    "    * 4 indicates asympotmatic\n",
    "* *trestbps* = resting blood pressure (measured in mm/Hg at time of hospital admission)\n",
    "* *chol* = serum cholestoral measured in mh/dl\n",
    "* *fbs* = fasting blood sugar (1 = fbs > 120 mg/dl, 0 = fbs &le; 120 mg/dl)\n",
    "* *restecg* = resting electrocardiographic results\n",
    "    * 0 indicates normal results\n",
    "    * 1 indicates having ST-T wave abnormality (T wave inversion and/or ST elevation or depression of > 0.05 mV)\n",
    "    * 2 indicates showing probable or definite left ventricular hypertrophy using the Estes' criteria\n",
    "* *thalach* = maximum heart rate achieved\n",
    "* *exang* = exercise induced angina (1 = yes, 0 = no)\n",
    "* *oldpeak* = ST depression induced by exercise relative to rest\n",
    "* *slope* = the slope of the peak exercise ST segment\n",
    "    * 1 indicates upsloping\n",
    "    * 2 indicates flat\n",
    "    * 3 indicate downsloping\n",
    "* *ca* = number of major vessels (0-3) colored by flourosopy\n",
    "* *thal* = (no description provided; 3 = normal, 6 = fixed defect, 7 = reversable defect)\n",
    "* *target* = presence of heart disease (1 = yes, 0 = no; 1 is the collapsed measure of any level of heart disease which typically ranges on a measure from 0-4 where 0 indicates no heart disease present)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n"
     ]
    },
    {
     "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>Variable</th>\n",
       "      <th>N</th>\n",
       "      <th>Mean</th>\n",
       "      <th>SD</th>\n",
       "      <th>SE</th>\n",
       "      <th>95% Conf.</th>\n",
       "      <th>Interval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>age</td>\n",
       "      <td>303.0</td>\n",
       "      <td>54.366337</td>\n",
       "      <td>9.082101</td>\n",
       "      <td>0.521753</td>\n",
       "      <td>53.339605</td>\n",
       "      <td>55.393069</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>trestbps</td>\n",
       "      <td>303.0</td>\n",
       "      <td>131.623762</td>\n",
       "      <td>17.538143</td>\n",
       "      <td>1.007540</td>\n",
       "      <td>129.641075</td>\n",
       "      <td>133.606450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chol</td>\n",
       "      <td>303.0</td>\n",
       "      <td>246.264026</td>\n",
       "      <td>51.830751</td>\n",
       "      <td>2.977599</td>\n",
       "      <td>240.404558</td>\n",
       "      <td>252.123495</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>thalach</td>\n",
       "      <td>303.0</td>\n",
       "      <td>149.646865</td>\n",
       "      <td>22.905161</td>\n",
       "      <td>1.315867</td>\n",
       "      <td>147.057435</td>\n",
       "      <td>152.236294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>oldpeak</td>\n",
       "      <td>303.0</td>\n",
       "      <td>1.039604</td>\n",
       "      <td>1.161075</td>\n",
       "      <td>0.066702</td>\n",
       "      <td>0.908344</td>\n",
       "      <td>1.170864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>ca</td>\n",
       "      <td>303.0</td>\n",
       "      <td>0.729373</td>\n",
       "      <td>1.022606</td>\n",
       "      <td>0.058747</td>\n",
       "      <td>0.613767</td>\n",
       "      <td>0.844979</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Variable      N        Mean         SD        SE   95% Conf.    Interval\n",
       "0       age  303.0   54.366337   9.082101  0.521753   53.339605   55.393069\n",
       "1  trestbps  303.0  131.623762  17.538143  1.007540  129.641075  133.606450\n",
       "2      chol  303.0  246.264026  51.830751  2.977599  240.404558  252.123495\n",
       "3   thalach  303.0  149.646865  22.905161  1.315867  147.057435  152.236294\n",
       "4   oldpeak  303.0    1.039604   1.161075  0.066702    0.908344    1.170864\n",
       "5        ca  303.0    0.729373   1.022606  0.058747    0.613767    0.844979"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Getting some information on the continuous variables\n",
    "\n",
    "rp.summary_cont(df[['age', 'trestbps', 'chol', 'thalach', 'oldpeak', 'ca']])"
   ]
  },
  {
   "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>Variable</th>\n",
       "      <th>Outcome</th>\n",
       "      <th>Count</th>\n",
       "      <th>Percent</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>sex</td>\n",
       "      <td>1</td>\n",
       "      <td>207</td>\n",
       "      <td>68.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td></td>\n",
       "      <td>0</td>\n",
       "      <td>96</td>\n",
       "      <td>31.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>cp</td>\n",
       "      <td>0</td>\n",
       "      <td>143</td>\n",
       "      <td>47.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td></td>\n",
       "      <td>2</td>\n",
       "      <td>87</td>\n",
       "      <td>28.71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>50</td>\n",
       "      <td>16.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td></td>\n",
       "      <td>3</td>\n",
       "      <td>23</td>\n",
       "      <td>7.59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>fbs</td>\n",
       "      <td>0</td>\n",
       "      <td>258</td>\n",
       "      <td>85.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>45</td>\n",
       "      <td>14.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>restecg</td>\n",
       "      <td>1</td>\n",
       "      <td>152</td>\n",
       "      <td>50.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td></td>\n",
       "      <td>0</td>\n",
       "      <td>147</td>\n",
       "      <td>48.51</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td></td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>1.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>exang</td>\n",
       "      <td>0</td>\n",
       "      <td>204</td>\n",
       "      <td>67.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>99</td>\n",
       "      <td>32.67</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>slope</td>\n",
       "      <td>2</td>\n",
       "      <td>142</td>\n",
       "      <td>46.86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>140</td>\n",
       "      <td>46.20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td></td>\n",
       "      <td>0</td>\n",
       "      <td>21</td>\n",
       "      <td>6.93</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>ca</td>\n",
       "      <td>0</td>\n",
       "      <td>175</td>\n",
       "      <td>57.76</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>65</td>\n",
       "      <td>21.45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td></td>\n",
       "      <td>2</td>\n",
       "      <td>38</td>\n",
       "      <td>12.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td></td>\n",
       "      <td>3</td>\n",
       "      <td>20</td>\n",
       "      <td>6.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td></td>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>1.65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>thal</td>\n",
       "      <td>2</td>\n",
       "      <td>166</td>\n",
       "      <td>54.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td></td>\n",
       "      <td>3</td>\n",
       "      <td>117</td>\n",
       "      <td>38.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td></td>\n",
       "      <td>1</td>\n",
       "      <td>18</td>\n",
       "      <td>5.94</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td></td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.66</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>target</td>\n",
       "      <td>1</td>\n",
       "      <td>165</td>\n",
       "      <td>54.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td></td>\n",
       "      <td>0</td>\n",
       "      <td>138</td>\n",
       "      <td>45.54</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Variable  Outcome  Count  Percent\n",
       "0       sex        1    207    68.32\n",
       "1                  0     96    31.68\n",
       "2        cp        0    143    47.19\n",
       "3                  2     87    28.71\n",
       "4                  1     50    16.50\n",
       "5                  3     23     7.59\n",
       "6       fbs        0    258    85.15\n",
       "7                  1     45    14.85\n",
       "8   restecg        1    152    50.17\n",
       "9                  0    147    48.51\n",
       "10                 2      4     1.32\n",
       "11    exang        0    204    67.33\n",
       "12                 1     99    32.67\n",
       "13    slope        2    142    46.86\n",
       "14                 1    140    46.20\n",
       "15                 0     21     6.93\n",
       "16       ca        0    175    57.76\n",
       "17                 1     65    21.45\n",
       "18                 2     38    12.54\n",
       "19                 3     20     6.60\n",
       "20                 4      5     1.65\n",
       "21     thal        2    166    54.79\n",
       "22                 3    117    38.61\n",
       "23                 1     18     5.94\n",
       "24                 0      2     0.66\n",
       "25   target        1    165    54.46\n",
       "26                 0    138    45.54"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Getting some information on the categorical variables\n",
    "rp.summary_cat(df[['sex', 'cp', 'fbs', 'restecg', 'exang', 'slope', \\\n",
    "                   'ca', 'thal', 'target']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a pretty even split of those with heart disease, 54.5%, and those without heart disease, 45.5%. Next, to split the data into our feature set and the outcome set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "features = df.iloc[:, :-1]\n",
    "target = df['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Normalization\n",
    "The features are measured on different scales and have a range of variation present. This is a good case to normalize the features to remove the scale units from the measurement. I will use feature standardization method, this can be expressed using the following formula:\n",
    "\n",
    "$$\n",
    "x'= \\frac{x - average(x)}{\\sigma}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feature_normalization(x):\n",
    "    x_prime = (x - x.mean()) / x.std()\n",
    "    \n",
    "    return x_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_norm = features.apply(feature_normalization)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Splitting the Data \n",
    "When splitting data it is best practice to split the data randomly. Scikit learn has a built-in method that allows the data to be split into train and test sets, but doesn't include the splitting into a validation set. It's an easy fix, just have to run the method twice.\n",
    "\n",
    "The data will be split into a 60% train, 20% validate, and 20% test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regular data set length= 303 \n",
      " Feature training set length= 181 which is 60% of the original \n",
      " Feature validate set length= 61 which is 20% of the original \n",
      " Feature test set length= 61 which is 20% of the original \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# The first spit - 60% train and 40% test\n",
    "features_train, features_test, target_train, target_test \\\n",
    "    = train_test_split(features_norm, target, test_size=0.4, random_state=1)\n",
    "\n",
    "# The second split - 20% validate and 20% test\n",
    "features_validate, features_test, target_validate, target_test \\\n",
    "    = train_test_split(features_test, target_test, \n",
    "                       test_size=0.5, random_state=1)\n",
    "\n",
    "print(\"Regular data set length=\", features.shape[0], \"\\n\", \n",
    "      \"Feature training set length=\", features_train.shape[0], f\"which is{features_train.shape[0]/features.shape[0]*100: .0f}% of the original\", \"\\n\",\n",
    "     \"Feature validate set length=\", features_validate.shape[0], f\"which is{features_validate.shape[0]/features.shape[0]*100: .0f}% of the original\", \"\\n\",\n",
    "     \"Feature test set length=\", features_test.shape[0], f\"which is{features_test.shape[0]/features.shape[0]*100: .0f}% of the original\", \"\\n\",)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression Model\n",
    "\n",
    "This logistic regression model will use the *newton-cg* solver with the default regularization penalty *l2*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using a Logistic Regression model, the accuracy achieved on the validation set was  80.3% and the f1 score was 0.8125.\n"
     ]
    }
   ],
   "source": [
    "classifier = LogisticRegression(random_state= 1, solver= 'newton-cg')\n",
    "classifier.fit(features_train, target_train)\n",
    "\n",
    "LogReg_y_predict = classifier.predict(features_test)\n",
    "LogReg_accuracy = accuracy_score(target_validate, classifier.predict(features_validate))\n",
    "LogReg_f1_score = f1_score(target_validate, classifier.predict(features_validate))\n",
    "\n",
    "print(\"Using a Logistic Regression model, the accuracy achieved on the validation set\",\n",
    "      f\"was {LogReg_accuracy*100: 0.1f}% and the f1 score was {LogReg_f1_score}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVM Model\n",
    "There are various kernals that can be used in the support vector machine algorithms. Scikit learn accepts the kernal options of: 'rbf' (the default), 'linear', 'poly', 'sigmoid', 'precomputed', or a callable function that the user created. With the outcome being binary, I will use the sigmoid function and test varying penalization parameters of C and gamma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the SVM Model\n",
    "There are a couple parameters that need to be decided on when using a support vector machine, the parameter *C and gamma*. Finding a good combination of these can be tricky, but they each effect the model's bias and variance. Easy way to test the combined effect of this is to use a for loop and compare accuracy metrics. In general: <br>\n",
    "\n",
    "**The C parameter** <br>\n",
    "This parameter is a penality parameter that is applied to the formula. A value of 1 means no penality\n",
    "> Large C value = lower bias and higher variance (more prone to overfitting) <br>\n",
    "> Small C value = higher bias and low variance (more prone to underfitting)\n",
    "\n",
    "**The gamma parameter** <br>\n",
    "By default, scikit learn svm uses $\\frac{1}{\\text{# of features * $\\sigma$}}$. This parameter effect the smoothness and overall shape of the peak of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using a SVM model, the accuracy achieved on the validation set was  80.3% and \n",
      "the f1 score was 0.8181818181818182.\n"
     ]
    }
   ],
   "source": [
    "clf = svm.SVC(gamma= 'scale', \n",
    "              C= 1, kernel= 'sigmoid', \n",
    "              max_iter= 3000)\n",
    "clf.fit(features_train, target_train)\n",
    "\n",
    "SVM_accuracy = accuracy_score(target_validate, clf.predict(features_validate))\n",
    "SVM_f1_score = f1_score(target_validate, clf.predict(features_validate))\n",
    "\n",
    "print(\"Using a SVM model, the accuracy achieved on the validation set\",\n",
    "      f\"was {SVM_accuracy*100: 0.1f}% and \\nthe f1 score was {SVM_f1_score}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial model achieved an accuracy of 80.3% with an F1 score of 0.82. Since the penalization parameter effects the cost of the function, I will test multiple C values and opt to use the best performing C value, in terms of accuracy and F1 score, in the final model selection. C values with an accuracy score $\\ge$ 80% will be displayed, and of those the best performing F1 score will be used to select the model for comparison against the logistic regression model and the 3-layer neural network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The C values that produced the largest accuracy score of  82.0% were: [24.84865656565657, 26.060737373737375, 27.272818181818185, 27.575838383838384, 27.878858585858588]. Of those C values, the best f1 score was 0.819672131147541 which was associated with the C value of 24.84865656565657\n"
     ]
    }
   ],
   "source": [
    "c = np.linspace(0.001, 30, num= 100)\n",
    "accuracy = {}\n",
    "f1 = {}\n",
    "\n",
    "for penality in c:\n",
    "    clf = svm.SVC(gamma= 'scale', \n",
    "                  C= penality, \n",
    "                  kernel= 'sigmoid', \n",
    "                  max_iter= 3000)\n",
    "    clf.fit(features_train, target_train)\n",
    "        \n",
    "    acc = accuracy_score(target_validate, clf.predict(features_validate))\n",
    "    f1_sc = f1_score(target_validate, clf.predict(features_validate))\n",
    "        \n",
    "    if acc >= 0.80000000:\n",
    "        accuracy.update({penality: acc})\n",
    "        f1.update({penality: f1_sc})\n",
    "            \n",
    "\n",
    "max_acc_value = max(accuracy.values())\n",
    "max_key = [key for key, value in accuracy.items() if value == max_acc_value]\n",
    "\n",
    "good_f1_values = [(value, key) for key, value in f1.items()]\n",
    "max_f1_value = max(good_f1_values)\n",
    "\n",
    "print(f\"The C values that produced the largest accuracy score of {max_acc_value*100: 0.1f}%\",\n",
    "      f\"were: {max_key}. Of those C values, the best f1 score was {max_f1_value[0]}\", \n",
    "      f\"which was associated with the C value of {max_f1_value[1]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several C values that produced accuracy scores that were greater than 80% - albeit not that much greater. More specifically, there were 5 C values that produced the maximum accuracy score of 82%. Of those C values, the best f1 score that was achieved corresponded to C = 24.85."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Network\n",
    "\n",
    "Here I will test a 3 layer neural network model on the data to see if it performs better than both the logistic regression and SVM models. In this neural network, I will use the sigmoid function as the activator again since the outcome is binary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = features_train.values # Converting Pandas Data Frame into Matrix\n",
    "x_train = np.append(np.ones([len(x_train), 1]), x_train, 1) # Adding bias paramter\n",
    "\n",
    "y_train = np.matrix(target_train) # Converting Pandas Data Frame into Matrix\n",
    "yVec_train = y_train.reshape((181,1)) # Reshape Marix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These functions are needed for the activation step (sigmoid) and\n",
    "#    for the gradient check in the neural network\n",
    "def sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z))\n",
    "\n",
    "def sigmoidGradient(z):\n",
    "    return np.multiply(sigmoid(z), (1 - sigmoid(z)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell is creating a function to randomly initialize $\\theta$'s to the correct input layer size with the corresponding output layer size. It can take a seed that will be used for the randomized number generation for reproducability purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Randomly Initializing Weights Function\n",
    "def randInitializeWeights(layer_size_in, layer_size_out, seed= 0):\n",
    "    \"\"\"\n",
    "    Randomly initializes the weights of a layer with a specified number \n",
    "        of incoming connections and outgoing connections. Where:\n",
    "            layer_size_in is the number of incoming connections, and\n",
    "            layer_size_out is the number of outgoing connections\n",
    "    \"\"\"\n",
    "    \n",
    "    np.random.seed(seed)\n",
    "    \n",
    "    epi = (6**1/2) / (layer_size_in + layer_size_out)**1/2\n",
    "    \n",
    "    W = np.random.rand(layer_size_out, layer_size_in + 1) * (2*epi) - epi\n",
    "    \n",
    "    return W\n",
    "\n",
    "initial_Theta1 = randInitializeWeights(features_norm.shape[1], 25)\n",
    "initial_Theta2 = randInitializeWeights(25, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next cell creates a function for a 3-layer neural network that can be regularized if desired. It calculates the forward propagation and backward error terms that will be used to train the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Neural Network Cost Function\n",
    "def nnCostFunction(Theta1, Theta2, num_labels, x, y, lmbda):\n",
    "    \n",
    "    Theta1_grad = 0\n",
    "    Theta2_grad = 0\n",
    "    \n",
    "    x = np.matrix(x)\n",
    "    yVec = np.matrix(y)\n",
    "    \n",
    "    m = x.shape[0]\n",
    "    J = 0\n",
    "    \n",
    "    \n",
    "    ## Forward Propogation\n",
    "    layer_z2 = x*Theta1.T\n",
    "    layer_a2 = sigmoid(layer_z2)\n",
    "\n",
    "    layer_a2 = np.append(np.ones([len(layer_a2), 1]), layer_a2, 1)\n",
    "\n",
    "    layer_z3 = layer_a2*Theta2.T\n",
    "    layer_a3 = sigmoid(layer_z3)\n",
    "\n",
    "    ## Cost\n",
    "    J = (1/m) * np.sum(np.sum(np.multiply((-1*yVec), np.log(layer_a3)) - \\\n",
    "                          np.multiply((1 - yVec), np.log(1 - layer_a3)))) \\\n",
    "    + lmbda/(2*m) * (np.sum(np.sum(np.power(Theta1[:, 1:], 2))) + \\\n",
    "                     np.sum(np.sum(np.power(Theta2[:, 1:], 2))))\n",
    "    \n",
    "    \n",
    "\n",
    "    # This computes the error between the last layer (contains prediction) and the known output \n",
    "    delta_3 = layer_a3 - yVec\n",
    "\n",
    "    # This computes the error between layer 2 and layer 3\n",
    "    delta_2 = np.multiply((delta_3 * Theta2), \\\n",
    "                          np.append(np.ones([len(sigmoidGradient(layer_z2)), 1]), \\\n",
    "                                    sigmoidGradient(layer_z2), 1))\n",
    "    delta_2 = delta_2[:, 1:]\n",
    "\n",
    "    Theta1_grad = Theta1_grad + (1/m) * (delta_2.T * x) + (lmbda/m) * np.append(np.zeros([len(Theta1), 1]), Theta1[:, 1:], 1)\n",
    "    Theta2_grad = Theta2_grad + (1/m) * (delta_3.T * layer_a2) + (lmbda/m) * np.append(np.zeros([len(Theta2), 1]), Theta2[:, 1:], 1)\n",
    "    \n",
    "    \n",
    "    return J, Theta1_grad, Theta2_grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient descent function in the next cell uses the previously created *nnCostFunction* to train the neural network using gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientDescentnn(x, y, Theta1, Theta2, alpha, num_iters,\n",
    "                     lmbda, num_labels):\n",
    "    \n",
    "    m = y.shape[0]\n",
    "    J_history = []\n",
    "    \n",
    "    for i in range(num_iters):\n",
    "        #nnCostFunction(theta1, theta2, num_labels, X, y, lmbda)\n",
    "        cost, Theta1_grad, Theta2_grad = nnCostFunction(Theta1, Theta2, \n",
    "                                                        num_labels, x, y, \n",
    "                                                        lmbda)\n",
    "        Theta1 = Theta1 - (alpha * Theta1_grad)\n",
    "        Theta2 = Theta2 - (alpha * Theta2_grad)\n",
    "        J_history.append(cost)\n",
    "    \n",
    "    return Theta1, Theta2, J_history"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to train the randomly initialized thetas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This works and runs\n",
    "Theta1_f, Theta2_f, J_history = \\\n",
    "    gradientDescentnn(x_train, yVec_train, \n",
    "                      initial_Theta1, initial_Theta2, \n",
    "                      0.8, 3000, 1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next step is to check that the algorithm learned through each iteration. If the algorithm learned properly there will be a clear decline in the cost function value (y-axis) as the number of iterations increased (x-axis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'Number of Iterations')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcXFWd9/HPt6rTaZYEggnIEghocEMFiYgjKm4YHRV93FAfhUcdBgUdt3HAQURQ3EZx9EERNeIysoxr8EERFRiUAdIsguAAIYJEtkDYIUkvv+ePcyp9u7q66nanq6u7832/XtV177nn3ntuVXX96pxz77mKCMzMzJqpdLoAZmY29TlYmJlZSw4WZmbWkoOFmZm15GBhZmYtOViYmVlLDhbjJOlCSe+epH29R9Jdkh6W9Lg27+tASddNdF4bG0k3SHp+G7a7o6TfS3pI0ucmevudJqlLUkha1OmybCpJP5B0fKfLUeNg0YSkWyQ9lr+k75L0HUlbj3Ebi/KHt2ucZZgFfAk4KCK2joh7C8uen8v2sKRH8n4eLjx2Hev+IuLCiHjaROcdD0mvkHRx/mK7Owfov9/Ebf5e0mFNlj+xwet4xabss0SZRnwpRMSTIuLiNuzuCOB2YG5E/MumbkzSu/Pr9cG69DslHbCp259o+f1/TNJOhbSlklaWXP9Tkk5vWwGnMAeL1l4dEVsDzwKeDRw7yfvfAegBRvyCj4iLcwDZGqh9aW9bS4uIvxbzS6pImhbvuaRDgLOAZcDOwI7AJ4HXTMb+C6/h1hGx72Tsc5LsBlwf47gat8kPnrXAMWP9ITUe4/3RVedRJv//eEwkVTtdhhEiwo9RHsAtwEsL818AfpGnLwTenacrpA/frcDdwPeAbfKyvwIBPJwfz22wn9nAl0m/+G7P07OBPYFHCuv/rklZF+V8XXXpvwdOBP4beCznezfwZ+Ah4ObaceT8LwVuKcyvBj4EXAs8AJwBzB5r3rz8GOBO4G/AP+TyLmpwLJWc54NNjrcCHFd4zU8n/VoG2BL4IXAvcD9wOTAf+BwwAKzLr+eXG2z3ienfouE+PwWcPlre/Fp/Ergkv7a/ArYrLH8BcGl+bW4D3g68F+gDNuQy/bTwWh6Yp3uArwB35NflS0B38T0APgqsyZ+fd4xS/u/X7evAktv+WH7fvtNgm+8m/S/8EvjXQvqdwAGF9+pjpM/aPcCZwLxGn6EGx/4p0o+GM/Jrehjw3Pw63p/L/RVgVs7fNdrnqvAeHZe3tXtOWwqsLOTZBfhpfj3/AhyZ01+VX7u+/PpdAbwMuKqw7oXAJYX5S4FX5emnARflcl8L/H0h3w+AU/Jn5pH83vwAOD4vnwv8F3AyoI58H3Zip9PlQSFYAAtJv+5PLHwoasHincBKYA9ga+AnwPfzskU0+BKv288J+UO1PbCA9GVzYtn1m+XL/xy3AE8BZuV/plfnsgp4MSmIPCPnbxQALgUeDzwOuLFw3GPJ+yrSF9lTgK1I//yjBYu98rKFTY738Lz93YE5wM/JX2bAkcDPgC2AKrAE2LrwehzWZLubGixuAhaTAtbFwKfyst1JX1Bvyu/BfGDvvGzjl0Lda3lgnj4pfyYW5M/IZcAnCu9BP/CJ/P6+hvRlM3eUYxi2r5LbPgnoBrZosL1asNiXVMPYNqcXg8VHgD+Qaog9wLcZ+v8oEyw2kD6zlfyePht4Tn4d98ifg6Ny/jLB4jBSgDk9p20MFvnzcjUpuHXn9/gW4CWjfAa2Iv34mJfz30kKYFsVlm2bl/2FFNRn5eN+GHhi4X25jxQIK6Qfiz8Ajs+fld76z8hkP6ZFk0SH/UzS/aQP2UWkf5x6bwO+FBGrIuJh0i/oQ8ZQZX4bcEJE3B0Ra0i/Tt8+AWWvWRYRf46Ivojoj4hzclkjIn4H/BZo1pn65Yi4M1J/yS+AvceR903At3M5HiEd42hqnfh3NMnzNuDfIuIvEfEQ6Z/7rbmZrY/0D/bEiBiIiN78vpQm6f7C4wNjWPXbEXFTRDwK/CdDx/+/gV9FxNn5PbgnIq4uuc23kb4o1kTE3aQfF8XPxzpSUOqLiOXAelKtdCK23Z+Xb4iIx0bbSERcQQoa/9xg8T8CH4uIv0XEOtIX4JvG0CT6+/yZHYyIxyJiRURcll/HVcBpwAtLbqvmJOB/SXpyXfr+pEB7Uj7mlaTgdkijjeTP8lWk/5/9gCtJtfjnAn9HavK7H3geKWB8Ib9PvyHVxorb/WlE/Hc+zvU5bWfS985/RMTxYzzGCTUR7X8z3WvzG9vMTqTmkJpbSa/tDiX30Wj9nUbJOx63FWckvQr4OOkXcIX0K3hFk/XvLEw/Cmw3jrw7kQJuwzLVqXXi79gkX6PXrJv0C/n0vPxsSXNJzS/HRkR/k30OExHbls1bp/74a+34C0nNMOOxIyOPdefC/D0RMTDKfjd123dFxIaS2/o4cImkL9el7wqcI2mwkBakmkwZ9Z/fJwNfJNVmtiT9r11Wcltp5xF3Svo66UfLdwqLdgN2zT8Qa6qkQDiai0jNRvfk6cdIwUt5HtLn8a+RqxFZ/Wvd6LP+GlKz5TebH1H7uWYxMW4nfchqdiX9IruL9E8xnvVvn7DSFcogaQvgR8BngB3yl+KvSR/sdrqD1BZcs7BJ3utJx//6JnkavWYbgDX5F+HxEfEU4ADgdaRf0FDu/RjNI6Qvp5rHj2Hd24AnjLKsVZnuYOSx/m0M+96UbZd+vSLiOuAcUs26aDXwsojYtvDoiYg7qXtNc228/vTw+jJ8A/gTqeY4l9QHMZ7P7+eAgxheU74NuKmurHMi4tWjlAWGgsUL8vRFpGDxQoaCxe3AQknFcpZ5rU8FLgD+n6QtGyyfNA4WE+MM4IOSds9nhJwEnJV/ya4BBkltq83WP1bSAknzSR/+H7SprLNJv8DXAAO5lvGSNu2r6GzgXZKelD/0Hx8tY0QMAh8Gjpd0qKS5+Uyu50s6NWc7A/hQPjV5DvBp4IyIGJT0Ykl75WaOB0nNUrVf3nfR/L1o5mrghZIWStoWOHoM6/4AWCrp9flagPmSnlmyTGcAx+V1FpBeu4n6fEz0to8nnbwwp5B2KnBS7VRuSdtLqp3V9j/AHEkvz6eJ1/pemplD+rX9iKSnkJq5xiwi1pJOJik2nf03sEHShyX1SKpKerqk2hlxdwGL6r70f0/qvN6H1Ol9DamvYwmp3wpSv1A/8GFJsyS9GHgl6f+iaTFJpzuvApZL6hnPsU4EB4uJsYzU1PFfpE6sdcD7AHLb9aeBP+T27/0brP8pUgfWNaSzJK7MaRMut59+kHS2x1rgDaS+hbaKiHOAr5Neo5tIHZ6Q2tcb5T8TeCvpi+d2UvPOCaSObEjV8rNI/4yrSJ3H/5SX7UQ6yeBB0kkJvyF9KUL6cnhLfi++NMbD+BXpdbuWdIbV8rIrRsRfSJ20/0J63a8Enp4Xfwt4pqT7JP2oweqfBP6Y93sNqcnlM2Ms+2gmdNu5jf8MhtfAvkR67X4r6SHSF+ezc/77SP8r3yX9yl7L8Ka8Rj4MHEp6z79B+hyM18kUftHnH3ivJPU/3EJqWvoG6Wwk8r66gbWSLs/rPER67a7J/RFB+nyszH135D6IVwMH521+BXhrRNzYqoB5e+8infX3U0mzN+F4x03Dm9DMJoekp5O+MGfnmoSZTWGuWdikkfQ6Sd1KQ5Z8Fvi5A4XZ9OBgYZPpSFIV/CZSU92RnS2OmZXlZigzM2vJNQszM2tpxlyUN3/+/Fi0aFGni2FmNq1cccUV90TEglb52hosJC0F/p10BeS3IuKzdctPBl6UZ7cEtq9dOSvpUIZGhvxURHy32b4WLVpEb2/vRBbfzGzGk3Rr61xtDBZ5iN1TSKMyrgZWSFoeEdfX8kTEBwv530e6qAVJ25EuzllCOgf6irzufe0qr5mZja6dfRb7kS5KWZXHljmTdEHKaN7C0IVTLwfOj4i1OUCcTxoZ0szMOqCdwWJnhg+MtZrhg2ZtJGk30hDOvxvLupIOl9QrqXfNmjUTUmgzMxupncGi0cBeo52newjwo8LImaXWjYjTImJJRCxZsKBl/4yZmY1TO4PFaoaPLLoLo4+keghDTVBjXdfMzNqsncFiBbA4j8TaTQoIIwZek/Qk0l2m/ruQfB5wkKR5kuaRhhE+r41lNTOzJtp2NlRE9Es6ivQlXyXdre06SScAvfmOXpA6ts8s3hQkItZKOpGhG/KckIcTNjOzDpgxw30sWbIkxnOdxQOP9nH6Jbfw4idvz9N32aYNJTMzm7okXRERS1rlmzFXcI+XKnDyb26kqyoHCzOzUWz2Y0PN7ZnFLvO24M93PNjpopiZTVmbfbAA2Gunbbhm9QOdLoaZ2ZTlYAE8a7dt+evaR7n34YZ3+DQz2+w5WACLd0j3ll91zyMdLomZ2dTUsoM7XyPxWmBRMX9EnNS+Yk2u3bZL95a/9d5Hefai7TpcGjOzqafM2VA/Jd0C8wpgoEXeaWmnbbcA4K4H13W4JGZmU1OZYLFbROzV9pJ0UM+sKlt2V7nvkQ2dLoqZ2ZRUps/iUklPbXtJOmzelt2sdbAwM2uoTM3iOcBVklYC60kjwkZEPKutJZtkj9u6m3sdLMzMGioTLF7b9lJMAXN6unhkfX+ni2FmNiW1bIaKiJuBLUi3R30Z0JPTZpQtu7t4ZMOM7L83M9tkLYNFHjn2bGDX/Dhb0nvbXbDJtlV31TULM7NRlGmGOhzYLyIeBpB0EnAJ8LV2FmyybTm7i0c3OFiYmTVS5mwoAX2F+T4a3/Z0Wks1CzdDmZk1UqZm8X3S6bM/zvOvA77bviJ1xpbdXTzWN8DAYFCtzLhYaGa2SVoGi4j4vKQLgOeTahRHRMSKFqtNO1t2VwF4rG+ArWdv9rf5MDMbZtRvRUlbRcQjkuYCN+RHbdnciJhRN4Do7kotchv6B2F2hwtjZjbFNPsJ/SPgFcB1QPHeq8rzu7axXJNudleqWWzoH+xwSczMpp5RO7gj4hX5eWFE7Fp4LIyIUoFC0lJJN0haKenoUfK8SdL1kq6T9MNC+oCkq/Nj+VgPbKyG1SzMzGyYMkOU/zoiDmqV1mC9KnAK6UK+1cAKScsj4vpCnsXAMcDzIuI+SdsXNvFYROw9hmPZJBuDxYDPiDIzq9esz6Ib6AF2kDSHodNl51KuCWo/YGVErMrbOxM4GLi+kOcfgFMi4j6AiLh7zEcwQbqrKVis63PNwsysXrPrLI4k9Vc8OT/XHucBp5bY9s7AbYX51TmtaE9gT0l/kHSppKWFZT2SenN6w/GpJB2e8/SuWbOmRJFGN3tWrWbhYGFmVm/UmkVEnAycLOkDEfHlcWy70cUKUTffBSwGDgR2AS6WtFdE3A/sGhG3S9oD+J2ka+vHpIqI04DTAJYsWVK/7TGZXXWfhZnZaMpcwb1O0ra1GUnzJB1eYr3VwMLC/C7A7Q3y/Dwi+iLiL6TTcxcDRMTt+XkVcCGwT4l9jps7uM3MRlcmWByRf+kDkPsX3lNivRXAYkm75/6PQ4D6s5p+BrwIQNJ8UrPUqhyQZhfSn8fwvo4J52BhZja6MpcqV4szkirArFYrRUR/HrH2vLyNZRFxnaQTgN6IWJ6XHSTpetL9vf85Iu6V9HfANyQNkgLaZ4tnUbVD7TqL9Q4WZmYjlAkW50s6g9SpHaRaxW/KbDwizgXOrUs7rjAdwIfyo5jnEuDpZfYxUWZVUxdLnzu4zcxGKBMs/hl4L/BBUqf1r4FvtLNQnTArd3D3D25SP7mZ2YxUZiDBAeCr+TFj1Uaa7XfNwsxshDJXcO8PfALYrZg/IvZsY7kmXVctWLhmYWY2QplmqO8AHwWuIHVCz0hdtWYo1yzMzEYoEywejIhz2l6SDqu6ZmFmNqoyweJ3kj4D/ARYX0uMiGvaVqoOqJ0NNeBgYWY2QplgcUDdM6RTaF8w8cXpHNcszMxGV+ZsqOdPRkE6bVal1mfhYGFmVq/M2VAfa5QeESdNfHE6p1IREvQPuoPbzKxemWao4hlQPcDfk4Yqn3FmVSpuhjIza6BMM9TnivOSPkcaAHDGqVbkU2fNzBooM+psvdnAEya6IFNBV1WuWZiZNVCmz+Iqhm5aVAV2BGZUf0VNV0Xu4DYza6BMn8UbCtP9wJ0RsX60zNNZV9V9FmZmjYzaDCXpYIB8K9M1EXFzRNw6UwMF1GoW7rMwM6vXrM/iE4XpC9tcjimhqypfwW1m1kCzYKFRpmesrkqFPgcLM7MRmvVZbCHp6aSA0pOnNwaNmTY2FKRmqAFflGdmNkKzYLEG+FqevqcwDTNwbChI11n0+WwoM7MRRg0WEzEmlKSlwL+TTrn9VkR8tkGeNwHHkwLQHyPirTn9UODYnO1TEfHdTS1PK7OqFfdZmJk1UObU2XGRVAVOAV4GrAZWSFoeEdcX8iwGjgGeFxH3Sdo+p29H6mBfQgoiV+R172tXeaFWs3AzlJlZvfFcwV3WfsDKiFgVERuAM4GD6/L8A3BKLQhExN05/eXA+RGxNi87H1jaxrICtT4L1yzMzOq1M1jsDNxWmF+d04r2BPaU9AdJl+Zmq7LrIulwSb2SetesWbPJBe6q+gpuM7NGWgYLSb8uk9Zo1QZp9d/EXcBi4EDgLcC3JG1bcl0i4rSIWBIRSxYsWFCiSM11VSoMhIOFmVm9UfssJHWThiTfQdIchr7A5wK7ltj2amBhYX4X4PYGeS6NiD7gL5JuIAWP1aQAUlz3whL73CSVigcSNDNrpFnN4kjSfSuenJ9rj/OAU0tsewWwWNLuOfAcAiyvy/Mz4EUAkuaTmqVW5X0cJGmepHnAQTmtrXydhZlZY81OnT0ZOFnSByLiy2PdcET0SzqK9CVfBZZFxHWSTgB6I2I5Q0HhetJNlv45Iu4FkHQiKeAAnBARa8dahrGqVoRPhjIzG6nMqbN/lTQnIh6SdDTwLOCkiLi61YoRcS5wbl3acYXpAD6UH/XrLgOWlSjfhHHNwsyssTJnQx2fA8XfAa8GzqJcM9S04z4LM7PGygSL2j24XwV8LSJ+TLpb3ozTVRGDDhZmZiOUaYa6Q9IpwCuAfXNndTuvz+iYqmsWZmYNlfnSfxNwEfDKfDX1fODotpaqQ3wFt5lZYy2DRUQ8DFwPHCjpCGBeRPyy7SXrANcszMwaK3MF91HA2aQL8XYFzpb03nYXrBOq7rMwM2uoTJ/F4cB+uYaBpJOASxh+f4sZoatScc3CzKyBMn0WAvoK833M0NusVuQ+CzOzRsrULL4PXCrpx3n+dUDbb0TUCV1V0e+L8szMRmgZLCLi85IuAJ5PqlEcERErWqw2LaU+i06Xwsxs6mk26uy+wPyIOC8HhxU5/VWS9omIqyarkJOlq+KahZlZI836LL4A3NQg/Ubg39pTnM6qSAwGPiPKzKxOs2CxICJW1SdGxI3Apt9paArqqqR+e98AycxsuGbBYosmy7ac6IJMBdVqDhauWZiZDdMsWPxO0ifrEyUdxyTcta4TNtYsHCzMzIZpdjbUh4Flkm4Eap3ZewPXAv+n3QXrhIpSsPCFeWZmwzW7U95DwBsl7Qk8LSd/PPdZzEi1moU7uM3MhitzncWNpDOgZrxqNbXKuWZhZjZcW+9LIWmppBskrcy3ZK1ffpikNZKuzo93F5YNFNKXt7OcNe6zMDNrrMxwH+MiqQqcArwMWA2skLQ8Iq6vy3pWRBzVYBOPRcTe7SpfI9WNfRa+MM/MrKjMEOWHNkj7VIlt7wesjIhVEbEBOBM4eOxFnDzVjX0WHS6ImdkUU6YZ6q2SDqnNSPoKsEuJ9XYGbivMr85p9V4v6RpJP5K0sJDeI6lX0qWSXltif5usq+qahZlZI2WaoV4HnCNpgHQf7kcj4rAS6zUaxry+M+Ac4IyIWJ/vwvdd4MV52a4RcbukPUjXfFwbETcP24F0OOl+G+y6664litRc1X0WZmYNjVqzkDRX0lxSQDkU+BjwKPCxnN7KaqBYU9gFuL2YISLujYj1efabwL6FZbfn51WkiwD3qd9BRJwWEUsiYsmCBZs+AknV11mYmTXUrBnqOuBP+XEJMJ/U51BLb2UFsFjS7pK6gUOAYWc1SdqxMPsa4M85fZ6k2Xl6PvA80n3A28o1CzOzxppdlLdQUoV0S9VLx7rhiOjP9+8+D6gCyyLiOkknAL0RsRx4v6TXAP3AWuCwvPpTgG9IGiQFtM82OItqwnV5bCgzs4aa9llExKCkLwP7j2fjEXEucG5d2nGF6WOAYxqsdwnw9PHsc1N4uA8zs8bKnA11vqQpfcrrROmqpJfDNQszs+HKnA11FLCNpPXAY6SznCIitmtryTrAfRZmZo2VCRbz216KKcJ9FmZmjZUZSHBA0iuBF+SkCyPiV+0tVmdUPNyHmVlDZYb7+DTwUWBVfny05HAf044HEjQza6xMM9SrgX0iYgBA0jLgSuDYdhasE9xnYWbWWNkhyotXbM9pR0GmAvdZmJk1VqZm8XngSkm/JZ0JdSBwXNM1pikP92Fm1liZDu4fSLoAeA4pWBwXEX9re8k6YOMQ5eFgYWZW1DJY5D6Ki4GLI2Jl+4vUObWL8voHHCzMzIrK9FmcCewOfDPfHvUsSUe2uVwdUXWfhZlZQ2WaoX4t6TfAs4CXAEeShhI/pc1lm3TuszAza6xMM9R5wDakIccvBvav3Wtiptl46qz7LMzMhinTDHUjaQjxxcCewBNr95qYaTZelDfgK7jNzIrKNEO9D0DSNsA7gO8D2wNbtLdok69ScTOUmVkjZZqhjgCeDzwbuAP4Hqk5asbxcB9mZo2VuShvHvA1YEVEbGhzeTrKfRZmZo2V6bN4UkT8oRgoJJ3eviJ1zlCfhYOFmVlRmWDxjOKMpCqpSWrGqbrPwsysoVGDhaR/kXQf8AxJayXdl+fXUHdf7SbbWCrphnwx39ENlh8maY2kq/Pj3YVlh0q6KT8OHcexjZkkKnKfhZlZvWZ9Fp8Hvgh8Btj4RV8bqryVXAM5BXgZsBpYIWl5RFxfl/WsiDiqbt3tgE8AS4AArsjr3ldm35uiq1Jxn4WZWZ1RaxaR9JNufPRG4F/yXfN2kbRviW3vB6yMiFW5v+NM4OCS5Xo5cH5ErM0B4nxgacl1N0m1ItcszMzqlOmz+CrwIuDtef5R4NQS6+0M3FaYX53T6r1e0jWSfiRp4VjWlXS4pF5JvWvWrClRpNaqFXkgQTOzOmWCxd9FxD8C6wAiYi3QXWI9NUir/xY+B1gUEc8AfgN8dwzrEhGnRcSSiFiyYMGCEkVqrVqRhyg3M6tTJlj0SaqQv6wlPQ4oMx7GamBhYX4XYNiYUhFxb0Ssz7PfJA1QWGrddumqiP5BD/dhZlZUJlicAvwYWCDpk8Dvgc+VWG8FsFjS7pK6gUOA5cUMknYszL4G+HOePg84SNI8SfOAg3Ja27nPwsxspDJjQ31P0hXAS0nNQ2+MiD+VWK9f0lGkL/kqsCwirpN0AtAbEcuB90t6DWmgwrXAYXndtZJOJAUcgBNy81fbuc/CzGykMsN9EBHXAdeNdeMRcS5112RExHGF6WOAY0ZZdxmwbKz73FTVinzqrJlZnVGDhaSHGOpUVmG6AsyKiBk7TLmboczMhhs1WETEnOK8pK2AI4D3Ar9oc7k6plKRh/swM6vTsoNb0lxJxwLXAguA50bEP7W9ZB3SVZEHEjQzq9OsGWo74IPA20j3sNh3Mobb6LSqh/swMxuhWQf3rcC9wLeBB4C3S0PXykXEV9pbtM5wn4WZ2UjNgsW/kzq1u0nNT5sF91mYmY3UrIP72MksyFSRaha+gtvMrKjMFdybFV/BbWY2koNFHfdZmJmN5GBRp+o+CzOzEZqdOvv+ZivO1LOhqhUx6GBhZjZMs7OhNpszoIq6XLMwMxuh2dlQH5/MgkwV7uA2Mxup5aizkmaThg5/GtBTS4+Iw9tXrM5xn4WZ2UhlOri/BywCXgVcBjyBfIvVmahaqbjPwsysTplgsWe+78TDEfFtYCmwV3uL1TnuszAzG6nUPbjz8/2SngLMAXZrX5E6y30WZmYjlblT3rfzfbA/QbpF6pbAcc1Xmb6qEv0e7sPMbJimwUJSFbgnD01+AbDrpJSqg6pVMeBYYWY2TNNmqIgYAD4w3o1LWirpBkkrJR3dJN8bJIWkJXl+kaTHJF2dH6eOtwxj5YEEzcxGKtMMdZ6kDwBnAY/UEiPiwWYr5VrJKcDLgNXACknLI+L6unxzgPeTzrQqujki9i5RvglVkTu4zczqlQkW/5ifP1xIC1o3Se0HrIyIVQCSzgQOBq6vy3ci8HngIyXK0nYeSNDMbKSWZ0NFxMIGjzJ9FzsDtxXmV+e0jSTtAyyMiF80WH93SVdJukjS8xvtQNLhknol9a5Zs6ZEkVpLfRYOFmZmRWWu4H5ro/SI+GGrVRutVthuBTiZdHV4vTuAXSPiXkn7Aj+T9LT6pq+IOA04DWDJkiUT8g3vmoWZ2UhlmqGKv+p7gBcDVwCtgsVqYGFhfhfg9sL8HNLFfRfme3s/Hlgu6TUR0QusB4iIKyTdDOwJ9JYo7yap5j6LiKB4z3Ezs81Zy2AREe8pzudrLk4vse0VwGJJuwN/Aw4BNtZSIuIBYH5huxcCH4mIXkkLgLURMSBpD2AxsKrEPjdZtZJa5gYDqo4VZmZAuZpFvYdIv/Kbioh+SUeRLuSrAssi4jpJJwC9EbG8yeovAE6Q1A8MAEdExNpxlHXMunKEGBgMqhVHCzMzKNdn8VOG+hoqpNFnf1Zm4xFxLnBuXVrDq78j4sDC9I+BH5fZx0SrBQj3W5iZDSlTs/i/hel+4NaIuKU9xem8au6nSEN+VDtbGDOzKaLMqbO/BW4ABiPiIuAOSVu1vWQdUqtZ+CJuM7MhLYOFpHcCy4Fv5aTdgJ+3s1CdVOuz8GCCZmZDygxR/n5gf+BBgIi4Edi+nYXqJPdZmJmNVCZYrIuIDbWZPObTjD1NaKjPwsHCzKymTLD4g6SPAj2SXkQaULDR8BwzgmsWZmYjlQkWHyVdW/E/wD8BvwX+tZ2F6qTidRZmZpaUuYJ7APh6fsx4FTdDmZmNMGqwkHQ+hYH/6kSl0ji0AAATpklEQVREvLw9Reqsrjzch2sWZmZDmtUsjm2QtoTULHVve4rTee6zMDMbadRgEREb71wn6XnAx4FtgKMi4pxJKFtHdDlYmJmN0LTPQtJLSEEigJMi4vxJKVUH1WoWvijPzGxIsz6LS0n3mPgCcHFOe0ZteURc0/bSdYCboczMRmpWs+gn3Rb1EODNDL8QL0jDiM84boYyMxupWZ/FAZNZkKnCNQszs5HKXJS3WaldlNfnYGFmtpGDRZ3ZXekeFuv7BjpcEjOzqcPBok7PrPSSrOv32VBmZjVlbqv6jAbJDwC3RcSM+0at1SzWuWZhZrZRmZrFt4ErgO8B3wd6gZ8CN+XrMEYlaamkGyStlHR0k3xvkBSSlhTSjsnr3SBp0oYWmZ1rFm6GMjMbUiZY3ATsGxF7R8QzgX2Bq4GXA18cbaV834tTgFcATwXeIumpDfLNId1gqXjF+FNJp+w+DVgKfC1vr+16ZuU+CzdDmZltVCZYPKV4AV5EXAs8KyJWtlhvP2BlRKzKN086Ezi4Qb4Tgc8D6wppBwNnRsT6iPgLsDJvr+163AxlZjZCmWBxs6SvSnpefnwFWClpNunCvdHsTLqor2Z1TttI0j7Awoiov5lSy3Xz+odL6pXUu2bNmhKH0tqsqqgI1vW5ZmFmVlMmWLyD9GV9NHAMcDtwKClQNOuzaHTr1Y0XL0iqACcDHx7ruhsTIk6LiCURsWTBggVNilKeJHpmVV2zMDMrKHPzo0eBz+VHvQearLoaWFiY34UUaGrmAHsBFyrdcOjxwHJJrymxblv1zKqyrt/Bwsyspsyps/sDnwB2K+aPiD1brLoCWCxpd+BvpA7rtxbWfwCYX9jPhcBHIqJX0mPADyV9CdgJWAxcXvKYNtnsroqboczMCloGC+A7pBseXQGU/rkdEf2SjgLOA6rAsoi4TtIJQG9ELG+y7nWSzgauJzV3HZlv7zop3AxlZjZcmWDx4HhvdhQR5wLn1qUdN0reA+vmPw18ejz73VSzuyo+ddbMrKBMsPidpM8APwHW1xJn6v0swDULM7N6ZYLFAXXPMIPvZwFpfKj17rMwM9uozNlQz5+MgkwlW3Z3cdeD61pnNDPbTDS7repbIuIMSe9vtDwivtK+YnXWnJ4uVt7d7HpDM7PNS7Oaxbz8PDFXu00jc3q6eGhdX6eLYWY2ZTS7rerX8vPHJ684U8Ocnlk8tK6fiCBfMGhmtlkrc1HefOCdwCKGX5R3ePuK1Vlze2bRPxis6xtki+5JGezWzGxKK3M21M+BS4HfM4aL8qazOT3pZXlwXZ+DhZkZ5YLFVhHRaLC/GasWLB5a18cOc3s6XBozs84rM+rsLyUd1PaSTCFzt5gFwIPrfEaUmRmUCxZHAL+S9LCktZLuk7S23QXrpLm5ZvHAYz4jyswMyjVDzW+dZWaZv/VsAO59eEOHS2JmNjU0uyhvcUTcRLoPdiMzdmyo7eekfoq7H/JV3GZm0LxmcTTwLuCUBstm9NhQW3RXmTO7i7sfXN86s5nZZqDZRXnvys+b3dhQAAvmznbNwswsK9NngaQnA08FNp5HGhE/bFehpoId5vRwl2sWZmZAuSu4jwUOAp5Muuvdy0kX6M3oYLFwuy244IY1nS6GmdmUUObU2TcDLwLuiIi3A8+kZI1kOnvCgq1Z89B6nz5rZka5YPFYvv91v6Q5wJ3AHu0tVuc9YcHWANy85uEOl8TMrPPKBIurJG0LLAN6gcuBK8tsXNJSSTdIWinp6AbLj5B0raSrJf1e0lNz+iJJj+X0qyWdOoZjmhBPevwcAP70twcme9dmZlNO0+YkpfG5j4+I+4FTJJ0HzI2IlsFCUpV02u3LgNXACknLI+L6QrYfRsSpOf9rgC8BS/OymyNi7zEf0QTZZd4W7LRND5etWss7nruoU8UwM5sSmtYsIiKAXxTmV5YJFNl+wMqIWBURG4AzgYPrtv9gYXYr0vUbU4Ik9t/jcVxy8z1s6Pf9uM1s81amGepySc8ax7Z3Bm4rzK/OacNIOlLSzcDngeItXHeXdJWkiyQ1vNZD0uGSeiX1rlkz8WcuvfqZO3Hfo32cf/1dE75tM7PpZNRgIanWRHUAKWDcIOnK/AVepnbR6BZzI2oOEXFKRDwB+Bfg2Jx8B7BrROwDfAj4oaS5DdY9LSKWRMSSBQsm/u6vL9hzAbvP34ov/voGHt3gEWjNbPPVrGZxeX5+LfAk4JXAG4E35OdWVgMLC/O7ALc3yX9m3hcRsT4i7s3TVwA3A3uW2OeEqlbEiQfvxS33PsJhy1Zw672PTHYRzMymhGYd3AKIiJvHue0VwGJJuwN/Aw4B3jpsB0ODFQL8PXBTTl8ArI2IAUl7AIuBVeMsxyY5YPF8Tn7z3hzzk2s58N8uZO+F27L3wm3ZZd6WPH5uD1v3dLFVd5Utu7vYsrtKd1eFroroqlboqopZlQrViphVle/nbWbTVrNgsUDSh0ZbGBFfarbhiOiXdBTpqu8qsCwirpN0AtAbEcuBoyS9FOgD7gMOzau/ADhBUj/pVq5HRETH7qFx8N4785zdH8cZl/+VP6y8hzMu/yvr+sbe6V2tKAWSHExmVZXTUkDpqqZl1Uot4NTmU55G87VtVivF7ZWbb7b9WdUKs6oVurvSerO7KnRXUzCspXV3VeiuVhwEzTYDSic8NVgg3QF8ncZ9D0TEJ9tYrjFbsmRJ9Pb2Tsq+IoL7H+3jzgfX8cj6fh7dMMCjG9Lzhv5B+geD/oH03DcQDAwO0jcQ9A8O0j8Qw5YPDJafH5pO2xo+n/ZTm0/LBhmchPPLuquF4FF75GAzu5A2q5rSi3mGPXcNrTN7VpWergo9s6rMzs/DpyvM7krPPbOqdFcrVCoOWmZjJemKiFjSKl+zmsUdEXHCBJZpxpDEvK26mbdVd6eL0tLgYBSCTgowfQOt5/sHUoDbMDDAhv5gw8AgG/prj4G8bJD1Oa2vuHyg7rl/kHV9gzy0rn/j/Pq8vLhe/yZGtu6uylCAmVWhp6txgJmdA0xP11C+2bNykOpKtafZXUOBbvaw5Xm+sKy7K9UMzWayln0WNr1VKqJ74xdZtaNlaWVwcCgAre8fYH3fIOv6Bljfn57X9aX0dfXp/UPLRq6Tph/d0M/aRwZZl/OsL6zTN7Dp1a9ZVdFdTTWi2YVa0lBzYt10tdYkWZhukL+a+726qqKi1EQokaaVpquVtKyi9H5vnFaDdSrDl1UqhenCfC2fJETxGUTaXqPpSp6mlp7XrRTyUttOXXperbCd4ftENEwfUZb68rqZdEI0CxYvmbRSmJG+6HoqqTYAsyZtv/0DqZazvm/oeX3/QA5aQ9O1GtH6voFCvlr6wMa8tZrUsGbHWtPiQPBof/+w5sRi02HfwODG5+J6k9GcOJNtDEYMDyQwFIjSdCE4jbLe8LzD822cL2y/GDiH8jZej8J6xQA4Yn/5T22tp+w4l//71vFcDldes5sfdaxD2WwypTPXKmw5hVsVBweDwQgGg/ycpgcGg6ibHqjlq19ncCjfYAQR5Lx5vcHh+QbztggIUv48u3Gf1NKpLUvptb7QKKxbS8+rDU8vbIecZ3AwCtsd2mat7FG3/RgtncblycUoTNeONc/HyHy59BvTN+ar287I7Rf2X5dvaH9D26e+3MO2X3gdc8Ku223Z8jO0qWb8UONmM0GlIipuGbYOKjPch5mZbeYcLMzMrCUHCzMza8nBwszMWnKwMDOzlhwszMysJQcLMzNrycHCzMxaGnXU2elG0hrg1k3YxHzgngkqTifNlOMAH8tUNVOOZaYcB2zasewWES1vNTpjgsWmktRbZpjeqW6mHAf4WKaqmXIsM+U4YHKOxc1QZmbWkoOFmZm15GAx5LROF2CCzJTjAB/LVDVTjmWmHAdMwrG4z8LMzFpyzcLMzFpysDAzs5Y2+2AhaamkGyStlHR0p8tThqRbJF0r6WpJvTltO0nnS7opP8/L6ZL0lXx810hq770XW5d9maS7Jf2pkDbmsks6NOe/SdKhU+Q4jpf0t/y+XC3plYVlx+TjuEHSywvpHf/8SVoo6QJJf5Z0naR/yunT8X0Z7Vim1XsjqUfS5ZL+mI/jkzl9d0mX5df3LEndOX12nl+Zly9qdXxjFvmWipvjA6gCNwN7AN3AH4GndrpcJcp9CzC/Lu3zwNF5+mjgc3n6lcAvSbfr3R+4rMNlfwHwLOBP4y07sB2wKj/Py9PzpsBxHA98pEHep+bP1mxg9/yZq06Vzx+wI/CsPD0HuDGXeTq+L6Mdy7R6b/Jru3WengVcll/rs4FDcvqpwHvy9HuBU/P0IcBZzY5vPGXa3GsW+wErI2JVRGwAzgQO7nCZxutg4Lt5+rvAawvp34vkUmBbSTt2ooAAEfFfQP393cda9pcD50fE2oi4DzgfWNr+0g8Z5ThGczBwZkSsj4i/ACtJn70p8fmLiDsi4so8/RDwZ2Bnpuf7MtqxjGZKvjf5tX04z87KjwBeDPwop9e/J7X36kfASySJ0Y9vzDb3YLEzcFthfjXNP1hTRQC/lnSFpMNz2g4RcQekfxhg+5w+HY5xrGWfysd0VG6aWVZrtmEaHUduvtiH9Et2Wr8vdccC0+y9kVSVdDVwNynw3gzcHxH9Dcq0sbx5+QPA45jA49jcg4UapE2Hc4mfFxHPAl4BHCnpBU3yTtdjhNHLPlWP6evAE4C9gTuAL+b0aXEckrYGfgx8ICIebJa1QdqUOp4GxzLt3puIGIiIvYFdSLWBpzQpU9uPY3MPFquBhYX5XYDbO1SW0iLi9vx8N/BT0gfprlrzUn6+O2efDsc41rJPyWOKiLvyP/gg8E2GqvtT/jgkzSJ9uf5HRPwkJ0/L96XRsUzn9yYi7gcuJPVZbCupq0GZNpY3L9+G1Ew6YcexuQeLFcDifIZBN6ljaHmHy9SUpK0kzalNAwcBfyKVu3b2yaHAz/P0cuAd+QyW/YEHak0LU8hYy34ecJCkebk54aCc1lF1fUGvI70vkI7jkHzGyu7AYuBypsjnL7dtfxv4c0R8qbBo2r0vox3LdHtvJC2QtG2e3gJ4Kan/5QLgDTlb/XtSe6/eAPwuUg/3aMc3dpPVuz9VH6QzO24ktQf+a6fLU6K8e5DObvgjcF2tzKT2yd8CN+Xn7WLorIpT8vFdCyzpcPnPIDUD9JF+9bxrPGUH3knqrFsJ/J8pchzfz+W8Jv+T7ljI/6/5OG4AXjGVPn/AAaSmiWuAq/PjldP0fRntWKbVewM8A7gql/dPwHE5fQ/Sl/1K4D+B2Tm9J8+vzMv3aHV8Y314uA8zM2tpc2+GMjOzEhwszMysJQcLMzNrycHCzMxacrAwM7OWHCxsypIUkr5YmP+IpOMnaNunS3pD65ybvJ835hFQL6hLX6Q8Yq2kvYujoE7APreV9N7C/E6SftRsHbNWHCxsKlsP/C9J8ztdkCJJ1TFkfxfw3oh4UZM8e5PO6R9LGbqaLN6WNAopkK74j4i2B0ab2RwsbCrrJ91b+IP1C+prBpIezs8HSrpI0tmSbpT0WUlvy/cGuFbSEwqbeamki3O+V+X1q5K+IGlFHnTuHwvbvUDSD0kXd9WX5y15+3+S9LmcdhzpIrFTJX2h0QHmq4NPAN6sdJ+FN+er9JflMlwl6eCc9zBJ/ynpHNJAkltL+q2kK/O+a6OifhZ4Qt7eF+pqMT2SvpPzXyXpRYVt/0TSr5TulfD5wutxej6uayWNeC9s89Ds14nZVHAKcE3ty6ukZ5IGXVtLuqfCtyJiP6Ub4bwP+EDOtwh4IWmAuQskPRF4B2n4imdLmg38QdKvc/79gL0iDfW8kaSdgM8B+wL3kb7IXxsRJ0h6Mek+Cr2NChoRG3JQWRIRR+XtnUQaruGdeciHyyX9Jq/yXOAZEbE21y5eFxEP5trXpZKWk+49sVekQehqo6/WHJn3+3RJT85l3TMv25s0Sut64AZJXyWNNLtzROyVt7Vt85feZirXLGxKizRi6PeA949htRWR7muwnjTMQe3L/lpSgKg5OyIGI+ImUlB5Mmk8o3coDQ19GWnIi8U5/+X1gSJ7NnBhRKyJNDz0f5BujjReBwFH5zJcSBrKYde87PyIqN1HQ8BJkq4BfkMaenqHFts+gDT0BRHxP8CtQC1Y/DYiHoiIdcD1wG6k12UPSV+VtBRoNhqtzWCuWdh08GXgSuA7hbR+8o+dPHhcd2HZ+sL0YGF+kOGf+fqxbmpDOr8vIoYNgCfpQOCRUcrXaBjoTSHg9RFxQ10ZnlNXhrcBC4B9I6JP0i2kwNJq26Mpvm4DQFdE3CfpmaQbGx0JvIk0/pNtZlyzsCkv/5I+m9RZXHMLqdkH0t3AZo1j02+UVMn9GHuQBlo7D3iP0jDXSNpTaXTfZi4DXihpfu78fgtw0RjK8RDpFqA15wHvy0EQSfuMst42wN05ULyIVBNotL2i/yIFGXLz066k424oN29VIuLHwMdJt5K1zZCDhU0XXwSKZ0V9k/QFfTlQ/4u7rBtIX+q/BI7IzS/fIjXBXJk7hb9Bixp4pOG5jyENH/1H4MqI+HmzdepcADy11sENnEgKftfkMpw4ynr/ASyR1EsKAP+Ty3Mvqa/lTw061r8GVCVdC5wFHJab60azM3BhbhI7PR+nbYY86qyZmbXkmoWZmbXkYGFmZi05WJiZWUsOFmZm1pKDhZmZteRgYWZmLTlYmJlZS/8fwOHBtEJfiQAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = pd.DataFrame(J_history).plot(title= \"Plot of Training Cost Function for Neural Network\",\n",
    "                             legend= False)\n",
    "\n",
    "ax.set_ylabel(\"Training Neural Network Cost Function\")\n",
    "ax.set_xlabel(\"Number of Iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Excellent, there is drastic decrease in the cost function and then a clear downward slope as each iteraction of the gradient descent neural network was completed; this indicates that the model was learning with each iteration.\n",
    "\n",
    "At the end of the 3,000 iterations, the cost of the neural network was $\\le$ 0.35. Next I test the trained theta values on the validation set for comparision to the other methods used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neural network's accuracy score=  80.3% \n",
      " Neural network's f1 score= 0.8181818181818182\n"
     ]
    }
   ],
   "source": [
    "x_validate = features_validate.values # Converting Pandas Data Frame into Matrix\n",
    "x_validate = np.append(np.ones([len(x_validate), 1]), x_validate, 1) # Adding bias paramter\n",
    "\n",
    "y_validate = np.matrix(target_validate) # Converting Pandas Data Frame into Matrix\n",
    "yVec_validate = y_validate.reshape((y_validate.shape[1],1)) \n",
    "\n",
    "\n",
    "def predict(Theta1, Theta2, x):\n",
    "    m = x.shape[0]\n",
    "    x = np.matrix(x)\n",
    "    \n",
    "    layer_z2 = x*Theta1.T # 5000x25, good\n",
    "    layer_a2 = sigmoid(layer_z2) # 5000x25, good\n",
    "\n",
    "    layer_a2 = np.append(np.ones([len(layer_a2), 1]), layer_a2, 1) # 5000x26, good\n",
    "\n",
    "    layer_z3 = layer_a2*Theta2.T # 5000x10\n",
    "    layer_a3 = sigmoid(layer_z3) # 5000x10\n",
    "    \n",
    "    return layer_a3\n",
    "\n",
    "\n",
    "predictions = predict(Theta1_f, Theta2_f, x_validate).round(0)\n",
    "\n",
    "print(f\"Neural network's accuracy score= {accuracy_score(yVec_validate, predictions) * 100: 0.1f}%\", \"\\n\",\n",
    "      f\"Neural network's f1 score= {f1_score(yVec_validate, predictions)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of Models on Image Classification Task\n",
    "\n",
    "Next is to create a nice table to easily compare the performance of the difference models on the cancer image classification task."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table>\n",
    "  <tr>\n",
    "    <th>Model Algorithm</th>\n",
    "    <th>Accuracy Score</th> \n",
    "    <th>F1 Score</th>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>Logistic Regression</td>\n",
    "    <td>80.3%</td>\n",
    "    <td>0.8125</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>SVM (3k iterations) with <br> C = 24.8, and <br> gamma= 'scale'</td>\n",
    "    <td>82.0%</td>\n",
    "    <td>0.8197</td>\n",
    "  </tr>\n",
    "  <tr>\n",
    "    <td>3 Layer Neural Network (3k iterations)<br> w/ sigmoid activation</td>\n",
    "    <td>80.3%</td>\n",
    "    <td>0.8182</td>\n",
    "  </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SVM model performed the best in terms of accuracy and F1 score. This model will be selected to be used on the final test set. The accuracy scores are not really desirable for this task given the severity of the problem, but this is a demonstration and is good enough for this purpose.\n",
    "\n",
    "Some steps that could be taken to improve performance are:\n",
    "* Increase the number of samples used\n",
    "    * The current data set had only 303 oservations which were then split into a training set (60% of observations), a validation set (20% of observations), and a test set (20% of observations). The overall sample size was small for this kind of application purpose.\n",
    "* Could increase the number of layers in the neural network to see if that increase classification performance\n",
    "* Could try different activation functions, lambda values, scaling parameters, and/or the method of normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The SVM Model on the Final Test Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using a SVM model, the accuracy achieved on the validation set was  67.2% and \n",
      "the f1 score was 0.7142857142857144.\n"
     ]
    }
   ],
   "source": [
    "clf = svm.SVC(gamma= 'scale', \n",
    "              C= 24.8, \n",
    "              kernel= 'sigmoid', \n",
    "              max_iter= 3000)\n",
    "clf.fit(features_test, target_test)\n",
    "\n",
    "SVM_accuracy = accuracy_score(target_test, clf.predict(features_test))\n",
    "SVM_f1_score = f1_score(target_test, clf.predict(features_test))\n",
    "\n",
    "print(\"Using a SVM model, the accuracy achieved on the validation set\",\n",
    "      f\"was {SVM_accuracy*100: 0.1f}% and \\nthe f1 score was {SVM_f1_score}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The selected model did not perform well on the final test set. Given the small number of observations, I would first seek to increase the sample size."
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}