HW05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f33c522c",
   "metadata": {},
   "source": [
    "# HW05: Practice with algorithm selection, assessment, hyperparameter tuning, multiclass and one-class classification, and imbalanced data.\n",
    "\n",
    "[Kara Conrad, 9082326472]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98d20274",
   "metadata": {},
   "source": [
    "##### Hello Students:\n",
    " Start by downloading HW05.ipynb from this folder. Then develop it into your solution.\n",
    "- Write code where you see \"... your code here ...\" below.\n",
    "  (You are welcome to use more than one cell.)\n",
    "- If you have questions, please ask them in class, office hours, or piazza. Our TA\n",
    "  and I are very happy to help with the programming (provided you start early\n",
    "  enough, and provided we are not helping so much that we undermine your learning).\n",
    "- When you are done, run these Notebook commands:\n",
    "  - Shift-L (once, so that line numbers are visible)\n",
    "  - Kernel > Restart and Run All (run all cells from scratch)\n",
    "  - Esc S (save)\n",
    "  - File > Download as > HTML\n",
    "- Turn in:\n",
    "  - HW03.ipynb to Canvas's HW03.ipynb assignment\n",
    "  - HW03.html to Canvas's HW03.html assignment\n",
    "  - As a check, download your files from Canvas to a new 'junk' folder. Try 'Kernel > Restart\n",
    "  and Run All' on the '.ipynb' file to make sure it works. Glance through the '.html' file.\n",
    "- Turn in partial solutions to Canvas before the deadline. e.g. Turn in part 1,\n",
    "  then parts 1 and 2, then your whole solution. That way we can award partial credit\n",
    "  even if you miss the deadline. We will grade your last submission before the deadline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fe65b353",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: imbalanced-learn in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (0.10.1)\n",
      "Requirement already satisfied: numpy>=1.17.3 in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (from imbalanced-learn) (1.21.5)\n",
      "Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (from imbalanced-learn) (2.2.0)\n",
      "Requirement already satisfied: joblib>=1.1.1 in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (from imbalanced-learn) (1.2.0)\n",
      "Requirement already satisfied: scipy>=1.3.2 in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (from imbalanced-learn) (1.7.3)\n",
      "Requirement already satisfied: scikit-learn>=1.0.2 in /Users/karaconrad/opt/anaconda3/lib/python3.9/site-packages (from imbalanced-learn) (1.0.2)\n"
     ]
    }
   ],
   "source": [
    " !pip install -U imbalanced-learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c8fd2f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn import mixture\n",
    "\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn import svm, linear_model, datasets\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "\n",
    "from sklearn.metrics import (confusion_matrix, precision_score, recall_score,\n",
    "                             accuracy_score, roc_auc_score, RocCurveDisplay)\n",
    "\n",
    "from sklearn.datasets import make_classification\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from imblearn.over_sampling import RandomOverSampler"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b5a813d",
   "metadata": {},
   "source": [
    "## 1. Algorithm selection for multiclass classification by optical recognition of handwritten digits\n",
    "\n",
    "The [digits](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_digits.html) dataset\n",
    "  has 1797 labeled images of hand-written digits.\n",
    "  - $X$ = `digits.data` has shape (1797, 64).\n",
    "    - Each image $\\mathbf{x}_i$ is represented as the $i$th row of 64 pixel values in the 2D\n",
    "      `digits.data` array that corresponds to an 8x8 photo of a handwritten digit.\n",
    "  - $y$ = `digits.target` has shape (1797,). Each $y_i$ is a number from 0 to 9 indicating\n",
    "    the handwritten digit that was photographed and stored in $\\mathbf{x}_i$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e018ed28",
   "metadata": {},
   "source": [
    "### 1(a) Load the digits dataset and split it into training, validation, and test sets as I did in the lecture example code [07ensemble.html](https://pages.stat.wisc.edu/~jgillett/451/burkov/07/07ensemble.html).\n",
    "This step does not need to display any output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e13337b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = datasets.load_digits()\n",
    "X = digits.data\n",
    "y = digits.target\n",
    "\n",
    "X_train, X_tmp, y_train, y_tmp = train_test_split(X, y, test_size=.2, random_state=0, stratify=y)\n",
    "\n",
    "X_valid, X_test, y_valid, y_test = train_test_split(X_tmp, y_tmp, test_size=.5, random_state=0, stratify=y_tmp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1610740",
   "metadata": {},
   "source": [
    "## 1(b) Use algorithm selection on training and validation data to choose a best classifier.\n",
    "Loop through these four classifiers and corresponding parameters, doing a grid search\n",
    "to find the best hyperparameter setting. Use only the training data for the grid search.\n",
    "- SVM:\n",
    "  - Try all values of `kernel` in 'linear', 'rbf'.\n",
    "  - Try all values of `C` in 0.01, 1, 100.\n",
    "- logistic regression:\n",
    "  - Use `max_iter=5000` to avoid a nonconvergence warning.\n",
    "  - Try all values of `C` in 0.01, 1, 100.\n",
    "- ID3 decision tree:\n",
    "  - Use `criterion='entropy` to get our ID3 tree.\n",
    "  - Try all values of `max_depth` in 1, 3, 5, 7.\n",
    "- kNN:\n",
    "  - (Use the default Euclidean distance).\n",
    "  - Try all values of `n_neighbors` in 1, 2, 3, 4.\n",
    "\n",
    "Hint:\n",
    "- Make a list of the four classifiers without setting any hyperparameters.\n",
    "- Make a list of four corresponding parameter dictionaries.\n",
    "- Loop through 0, 1, 2, 3:\n",
    "  - Run grid search on the $i$th classifier with the $i$th parameter dictionary on the\n",
    "    training data. (The grid search does its own cross-validation using the training data.)\n",
    "  - Use the $i$th classifier with its best hyperparameter settings (just `clf` from\n",
    "    `clf = GridSearchCV(...)`) to find the accuracy of the model on the validation data, i.e.\n",
    "    find `clf.score(X_valid, y_valid)`.\n",
    "- Keep track, as your loop progresses, of:\n",
    "  - the index $i$ of the best classifier (initialize it to `-1` or some other value)\n",
    "  - the best accuracy score on validation data (initialize it to `-np.Inf`)\n",
    "  - the best classifier with its hyperparameter settings, that is the best `clf` from\n",
    "    `clf = GridSearchCV(...)` (initialize it to `None` or some other value)\n",
    "\n",
    "I needed about 30 lines of code to do this. It took a minute to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9d0208f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best classifier is classifier 0\n",
      "Best accuracy is 0.9888888888888889\n",
      "Best classifier hyperparameters are {'C': 1, 'kernel': 'rbf'}\n"
     ]
    }
   ],
   "source": [
    "classifiers = [svm.SVC(),  linear_model.LogisticRegression(max_iter=5000), \n",
    "               DecisionTreeClassifier(criterion='entropy'), KNeighborsClassifier(metric='euclidean')]\n",
    "\n",
    "params = [{'kernel': ['linear', 'rbf'], 'C': [0.01, 1, 100]},\n",
    "          {'C': [0.01, 1, 100]},\n",
    "          {'max_depth': [1, 3, 5, 7]},\n",
    "          {'n_neighbors': [1, 2, 3, 4]}]\n",
    "\n",
    "best_classifier_idx = -1\n",
    "best_accuracy = -np.Inf\n",
    "best_clf = None\n",
    "\n",
    "# Loop through the classifiers and corresponding parameter grids\n",
    "for i in range(len(classifiers)):\n",
    "    clf = GridSearchCV(classifiers[i], params[i], cv=5)\n",
    "    clf.fit(X_train, y_train)\n",
    "    accuracy = clf.score(X_valid, y_valid)\n",
    "    \n",
    "    if accuracy > best_accuracy:\n",
    "        best_accuracy = accuracy\n",
    "        best_i = i\n",
    "        best_clf = clf\n",
    "        \n",
    "print(\"Best classifier is classifier\", best_i)\n",
    "print(\"Best accuracy is\", best_accuracy)\n",
    "print(\"Best classifier hyperparameters are\", best_clf.best_params_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bf500da",
   "metadata": {},
   "source": [
    "### 1(c) Use the test data to evaluate your already-fit best classifier and its hyperparameter settings from 1(b).\n",
    "- Well, there are two tied for 'best'. Please use the first of these two.\n",
    "- Report the result of calling `.score(X_test, y_test)` on your best classifier/hyperparameters.\n",
    "- Show a confusion matrix from the true `y_test` values and the corresponding $\\hat{y}$ values\n",
    "  predicted by your best classifier/hyperparameters on `X_test`.\n",
    "- For each of the wrong predictions (where `y_test` and your $\\hat{y}$ values disagree), show:\n",
    "  - The index $i$ in the test data of that example $\\mathbf{x}$\n",
    "  - The correct label $y_i$\n",
    "  - Your incorrect prediction $\\hat{y}_i$\n",
    "  - A plot of that image (to see whether the confusion was reasonable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c015c000",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 0.9777777777777777\n",
      "Confusion matrix:\n",
      " [[18  0  0  0  0  0  0  0  0  0]\n",
      " [ 0 18  0  0  0  0  0  0  0  0]\n",
      " [ 0  0 18  0  0  0  0  0  0  0]\n",
      " [ 0  0  0 18  0  0  0  1  0  0]\n",
      " [ 0  0  0  0 17  0  0  0  1  0]\n",
      " [ 0  0  0  0  0 18  0  0  0  0]\n",
      " [ 0  0  0  0  0  0 18  0  0  0]\n",
      " [ 0  0  0  0  0  0  0 18  0  0]\n",
      " [ 0  1  0  0  0  0  0  0 16  0]\n",
      " [ 0  0  0  0  0  1  0  0  0 17]]\n",
      "Index: 61\n",
      "Correct label: 9\n",
      "Incorrect prediction: 5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAKx0lEQVR4nO3d24td9RnG8efpqLQeB1pbJBMyChIIhcxICEhATGxLrOLMRS8SUEgo5EpxaEG0d/0HZHpRhCFqBFOljQoiVitosEJrzWGnNZlY0mDJNNooZTwFGhLfXsxOiXbsrL33Os2b7weCc9jk927061qzZ6/1c0QIQB5fa3oAAOUiaiAZogaSIWogGaIGkrmkir/UdsqX1NesWVPreqdPn65trRMnTtS21rlz52pbK7OI8GJfdxW/0soadafTSbve1NRUbWvNz8/XtlZmXxU1p99AMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDKFora92fY7to/ZfrDqoQD0b8mobQ9J+qWk2yWtkbTVdr1vggZQWJEj9XpJxyLieESckfS0pIlqxwLQryJRr5B04SU8c92vfYHtHbb32d5X1nAAelfk0svFrgT5n6uwImJG0oyU9yotYDkocqSek7Tygs9HJJ2sZhwAgyoS9VuSbrR9ve3LJG2R9Hy1YwHo15Kn3xFx1va9kl6WNCTpsYg4XPlkAPpS6HZGEfGipBcrngVACXhHGZAMUQPJEDWQDFEDyRA1kAxRA8kQNZBMJdvu1Glior4LxtauXVvbWpI0PDxc21rbtm2rba3p6ena1roYcaQGkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiCZIjt0PGb7lO236xgIwGCKHKl3Sdpc8RwASrJk1BHxuqR/1TALgBKUdpWW7R2SdpT19wHoT2lRs+0O0A68+g0kQ9RAMkV+pfWUpD9IWm17zvaPqx8LQL+K7KW1tY5BAJSD028gGaIGkiFqIBmiBpIhaiAZogaSIWogmWW/7c5HH33U9AiV2bVrV21rTU1N1bYW2+5UiyM1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJFLlH2Urbr9metX3Y9v11DAagP0Xe+31W0k8j4oDtqyTtt/1KRBypeDYAfSiy7c57EXGg+/EnkmYlrah6MAD96ekqLdujksYlvbnI99h2B2iBwlHbvlLSM5KmIuLjL3+fbXeAdij06rftS7UQ9O6IeLbakQAMosir35b0qKTZiHi4+pEADKLIkXqDpHskbbLd6f75YcVzAehTkW133pDkGmYBUALeUQYkQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMo4o/9qLOi/oGB0drWsp7d27t7a1JGnVqlW1rleX7du317ZWnfuR1S0iFn1TGEdqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiCZIjce/LrtP9k+1N125+d1DAagP0Xu+/1vSZsi4tPurYLfsP3biPhjxbMB6EORGw+GpE+7n17a/cPN+oGWKnoz/yHbHUmnJL0SEYtuu2N7n+19Jc8IoAeFoo6IcxExJmlE0nrb313kMTMRsS4i1pU8I4Ae9PTqd0TMS9oraXMVwwAYXJFXv6+1Pdz9+BuSvifpaMVzAehTkVe/r5P0hO0hLfxP4NcR8UK1YwHoV5FXv/+shT2pASwDvKMMSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSW/bY7dRoeHq51vbGxsdrWmpycrG2tW2+9NeVakjQ/P1/bWmy7A1wkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSKZw1N0b+h+0zU0HgRbr5Uh9v6TZqgYBUI6i2+6MSLpD0s5qxwEwqKJH6mlJD0j6/KsewF5aQDsU2aHjTkmnImL//3sce2kB7VDkSL1B0l2235X0tKRNtp+sdCoAfVsy6oh4KCJGImJU0hZJr0bE3ZVPBqAv/J4aSKbIBnn/FRF7tbCVLYCW4kgNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJMO2O6hdp9Opba3p6ena1pKkXbt21bYW2+4AFwmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSKXQ7o+6dRD+RdE7SWW4DDLRXL/co2xgRH1Y2CYBScPoNJFM06pD0O9v7be9Y7AFsuwO0Q9HT7w0RcdL2tyW9YvtoRLx+4QMiYkbSjMSll0CTCh2pI+Jk95+nJD0naX2VQwHoX5EN8q6wfdX5jyX9QNLbVQ8GoD9FTr+/I+k52+cf/6uIeKnSqQD0bcmoI+K4pLU1zAKgBPxKC0iGqIFkiBpIhqiBZIgaSIaogWSIGkiGbXd6MDExUet64+Pjta118ODB2tbauHFjbWtNTk7WtpYkjY6O1rYW2+4AFwmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSKRS17WHbe2wftT1r++aqBwPQn6L3/f6FpJci4ke2L5N0eYUzARjAklHbvlrSLZK2SVJEnJF0ptqxAPSryOn3DZI+kPS47YO2d3bv//0FbLsDtEORqC+RdJOkRyJiXNJnkh788oMiYiYi1rHNLdCsIlHPSZqLiDe7n+/RQuQAWmjJqCPifUknbK/ufuk2SUcqnQpA34q++n2fpN3dV76PS9pe3UgABlEo6ojoSOJnZWAZ4B1lQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSTDXlo96HQ6TY9QmTr3gLrmmmtqW+vQoUO1rSVJY2Njta3FXlrARYKogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkhmyahtr7bdueDPx7anapgNQB+WvEdZRLwjaUySbA9J+oek56odC0C/ej39vk3S3yLi71UMA2BwRW8RfN4WSU8t9g3bOyTtGHgiAAMpfKTu3vP7Lkm/Wez7bLsDtEMvp9+3SzoQEf+sahgAg+sl6q36ilNvAO1RKGrbl0v6vqRnqx0HwKCKbrtzWtI3K54FQAl4RxmQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyVS17c4Hknq9PPNbkj4sfZh2yPrceF7NWRUR1y72jUqi7oftfVmv8Mr63Hhe7cTpN5AMUQPJtCnqmaYHqFDW58bzaqHW/EwNoBxtOlIDKAFRA8m0Imrbm22/Y/uY7QebnqcMtlfafs32rO3Dtu9veqYy2R6yfdD2C03PUibbw7b32D7a/Xd3c9Mz9arxn6m7GwT8VQu3S5qT9JakrRFxpNHBBmT7OknXRcQB21dJ2i9pcrk/r/Ns/0TSOklXR8SdTc9TFttPSPp9ROzs3kH38oiYb3isnrThSL1e0rGIOB4RZyQ9LWmi4ZkGFhHvRcSB7sefSJqVtKLZqcphe0TSHZJ2Nj1LmWxfLekWSY9KUkScWW5BS+2IeoWkExd8Pqck//GfZ3tU0rikNxsepSzTkh6Q9HnDc5TtBkkfSHq8+6PFTttXND1Ur9oQtRf5Wprfs9m+UtIzkqYi4uOm5xmU7TslnYqI/U3PUoFLJN0k6ZGIGJf0maRl9xpPG6Kek7Tygs9HJJ1saJZS2b5UC0Hvjogst1feIOku2+9q4UelTbafbHak0sxJmouI82dUe7QQ+bLShqjfknSj7eu7L0xskfR8wzMNzLa18LPZbEQ83PQ8ZYmIhyJiJCJGtfDv6tWIuLvhsUoREe9LOmF7dfdLt0ladi9s9rpBXuki4qzteyW9LGlI0mMRcbjhscqwQdI9kv5iu9P92s8i4sXmRkIB90na3T3AHJe0veF5etb4r7QAlKsNp98ASkTUQDJEDSRD1EAyRA0kQ9RAMkQNJPMfrdijXuhcWvEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index: 65\n",
      "Correct label: 4\n",
      "Incorrect prediction: 8\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAKpElEQVR4nO3d72ud9RnH8c9nUdn8GdjckKasKqMggzVSClIQU7dRp9g82IMWFCqDPlKUDUT3bP+AZA+GEKpOsFO2+gMRpxM0OGFztjXbrKmjK45m1UWZ0epgpXrtQU5HdXH5nnPuX7n6fkExOTnkex3r2/vOyTn31xEhAHl8oe0BAFSLqIFkiBpIhqiBZIgaSOasOr6p7ZRPqV9++eWNrnf8+PHG1lpYWGhsLVQjIrzc7a7jV1pZo37iiScaXW9mZqaxtaamphpbC9X4vKg5/QaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkimK2vZW22/YPmz7rrqHAjC4FaO2PSLpZ5Kuk3SFpB22r6h7MACDKTlSb5J0OCKORMQJSY9I2lbvWAAGVRL1GklHT/t8vnfbp9jeZXuf7X1VDQegfyVvvVzunSD/8y6siJiWNC3lfZcWsBqUHKnnJa097fMxScfqGQfAsEqifkXSN2xfavscSdslPVnvWAAGteLpd0SctH2rpGcljUi6PyIO1j4ZgIEUXc4oIp6W9HTNswCoAK8oA5IhaiAZogaSIWogGaIGkiFqIBmiBpJZ9Tt0bNvW3BvGmt6hY3x8vLG1FhcXG1uryZ1HJicnG1tLkmZnZxtbix06gDMEUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyZTs0HG/7QXbrzUxEIDhlBypfy5pa81zAKjIilFHxIuS/tnALAAqUHQ10RK2d0naVdX3AzCYyqJm2x2gG3j2G0iGqIFkSn6l9bCk30lab3ve9g/qHwvAoEr20trRxCAAqsHpN5AMUQPJEDWQDFEDyRA1kAxRA8kQNZBMZa/9bsvExERja73//vuNrSU1u4VLk1vhjI6ONrZWk/8Ou4IjNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyZRco2yt7Rdsz9k+aPv2JgYDMJiS136flPSjiDhg+wJJ+20/FxGv1zwbgAGUbLvzVkQc6H18XNKcpDV1DwZgMH29S8v2Oknjkl5e5mtsuwN0QHHUts+X9KikOyLig89+nW13gG4oevbb9tlaCnpPRDxW70gAhlHy7Lcl3SdpLiLuqX8kAMMoOVJvlnSzpC22Z3t/vlfzXAAGVLLtzkuS3MAsACrAK8qAZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSGbV76XV5F5JF110UWNrSdLi4mJjazX52CYnJxtb60zEkRpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSKbkwoNftP0H23/sbbvzkyYGAzCYkpeJ/lvSloj4sHep4Jds/zoifl/zbAAGUHLhwZD0Ye/Ts3t/uFg/0FGlF/MfsT0raUHScxGx7LY7tvfZ3lfxjAD6UBR1RHwcERskjUnaZPuby9xnOiI2RsTGimcE0Ie+nv2OiEVJM5K21jEMgOGVPPt9se3R3sdfkvRtSYdqngvAgEqe/b5E0oO2R7T0P4FfRsRT9Y4FYFAlz37/SUt7UgNYBXhFGZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJeOmdlRV/UzvlWzN37tzZ6HpTU1ONrTUzM9PYWmy7U42I8HK3c6QGkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiCZ4qh7F/R/1TYXHQQ6rJ8j9e2S5uoaBEA1SrfdGZN0vaTd9Y4DYFilR+opSXdK+uTz7sBeWkA3lOzQcYOkhYjY///ux15aQDeUHKk3S7rR9puSHpG0xfZDtU4FYGArRh0Rd0fEWESsk7Rd0vMRcVPtkwEYCL+nBpIp2SDvvyJiRktb2QLoKI7UQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJsu9OH0dHRRtd77733GltrYmKisbWa3OInM7bdAc4QRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJFN0OaPelUSPS/pY0kkuAwx0Vz/XKJuIiHdrmwRAJTj9BpIpjTok/cb2ftu7lrsD2+4A3VB6+r05Io7Z/qqk52wfiogXT79DRExLmpbyvvUSWA2KjtQRcaz3zwVJj0vaVOdQAAZXskHeebYvOPWxpO9Keq3uwQAMpuT0+2uSHrd96v6/iIhnap0KwMBWjDoijkj6VgOzAKgAv9ICkiFqIBmiBpIhaiAZogaSIWogGaIGkmHbnT5MTU01ut6GDRsaW+uaa65pbC1Ug213gDMEUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRRFbXvU9l7bh2zP2b6q7sEADKb0ut8/lfRMRHzf9jmSzq1xJgBDWDFq2xdKulrSTkmKiBOSTtQ7FoBBlZx+XybpHUkP2H7V9u7e9b8/hW13gG4oifosSVdKujcixiV9JOmuz94pIqYjYiPb3ALtKol6XtJ8RLzc+3yvliIH0EErRh0Rb0s6ant976ZrJb1e61QABlb67Pdtkvb0nvk+IumW+kYCMIyiqCNiVhI/KwOrAK8oA5IhaiAZogaSIWogGaIGkiFqIBmiBpIhaiCZ0leUQc3vNzU7O9voesiBIzWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kMyKUdteb3v2tD8f2L6jgdkADGDFl4lGxBuSNkiS7RFJf5f0eL1jARhUv6ff10r6a0T8rY5hAAyv3zd0bJf08HJfsL1L0q6hJwIwlOIjde+a3zdK+tVyX2fbHaAb+jn9vk7SgYj4R13DABheP1Hv0OecegPojqKobZ8r6TuSHqt3HADDKt1251+SvlzzLAAqwCvKgGSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkjGEVH9N7XfkdTv2zO/IundyofphqyPjcfVnq9HxMXLfaGWqAdhe1/Wd3hlfWw8rm7i9BtIhqiBZLoU9XTbA9Qo62PjcXVQZ36mBlCNLh2pAVSAqIFkOhG17a2237B92PZdbc9TBdtrbb9ge872Qdu3tz1TlWyP2H7V9lNtz1Il26O299o+1Pu7u6rtmfrV+s/UvQ0C/qKlyyXNS3pF0o6IeL3VwYZk+xJJl0TEAdsXSNovaXK1P65TbP9Q0kZJF0bEDW3PUxXbD0r6bUTs7l1B99yIWGx5rL504Ui9SdLhiDgSESckPSJpW8szDS0i3oqIA72Pj0uak7Sm3amqYXtM0vWSdrc9S5VsXyjpakn3SVJEnFhtQUvdiHqNpKOnfT6vJP/xn2J7naRxSS+3PEpVpiTdKemTlueo2mWS3pH0QO9Hi922z2t7qH51IWovc1ua37PZPl/So5LuiIgP2p5nWLZvkLQQEfvbnqUGZ0m6UtK9ETEu6SNJq+45ni5EPS9p7Wmfj0k61tIslbJ9tpaC3hMRWS6vvFnSjbbf1NKPSltsP9TuSJWZlzQfEafOqPZqKfJVpQtRvyLpG7Yv7T0xsV3Sky3PNDTb1tLPZnMRcU/b81QlIu6OiLGIWKelv6vnI+KmlseqRES8Lemo7fW9m66VtOqe2Ox3g7zKRcRJ27dKelbSiKT7I+Jgy2NVYbOkmyX92fZs77YfR8TT7Y2EArdJ2tM7wByRdEvL8/St9V9pAahWF06/AVSIqIFkiBpIhqiBZIgaSIaogWSIGkjmP9cam9d+sg6SAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index: 152\n",
      "Correct label: 8\n",
      "Incorrect prediction: 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAKpElEQVR4nO3d34tc9RnH8c+nq6H1N7S2SDZ0FSQQCt1ICEhAbGxLrKK56EUCiiuFXClKC6K96z8g6UURlmgimCptVBSxWkFXK7TWJK6tycaSxg3ZRhulrL8KDdGnFzuBaNfud2bOr318vyC4Ozvsecb49pw5M3O+jggByOMrbQ8AoFpEDSRD1EAyRA0kQ9RAMmfV8Uttc0q9AitWrGhsW2vWrGlsW2+99VZj23r//fcb21bTIsKL3e46XtIi6mqMjY01tq3p6enGtnXLLbc0tq0nnniisW017Yui5vAbSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkimKGrbm2y/afuw7bvrHgrA4JaM2vaIpF9JulbSGklbbTf3RmEAfSnZU6+XdDgijkTESUmPSLqx3rEADKok6pWSjp3x/Vzvts+wvc32Xtt7qxoOQP9KPnq52CdB/udTWBExKWlS4lNaQJtK9tRzklad8f2opOP1jANgWCVRvyrpctuX2l4haYukJ+sdC8Cgljz8johTtm+T9KykEUkPRMSB2icDMJCiyxlFxNOSnq55FgAV4B1lQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDKs0NGHJlfMaHp7ExMTjW2rycd19dVXN7atprFCB/AlQdRAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDIlK3Q8YPuE7TeaGAjAcEr21Lskbap5DgAVWTLqiHhJ0r8amAVABYquJlrC9jZJ26r6fQAGU1nULLsDdANnv4FkiBpIpuQlrYcl/VHSattztn9S/1gABlWyltbWJgYBUA0Ov4FkiBpIhqiBZIgaSIaogWSIGkiGqIFkKnvvd1vGx8cb29bU1FRj25Kk+fn5xrbV5GNr+t/jlw17aiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkim5Rtkq2y/YnrF9wPYdTQwGYDAl7/0+JelnEbHf9vmS9tl+LiIO1jwbgAGULLvzdkTs7339oaQZSSvrHgzAYPr6lJbtMUlrJb2yyM9YdgfogOKobZ8n6VFJd0bEB5//OcvuAN1QdPbb9tlaCHp3RDxW70gAhlFy9tuS7pc0ExH31j8SgGGU7Kk3SLpZ0kbb070/P6p5LgADKll252VJbmAWABXgHWVAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJLPs19LavHlz2yPUpsl1wnbt2tXYtlhLq17sqYFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZEouPPhV23+2/Xpv2Z1fNDEYgMGUvE30P5I2RsRHvUsFv2z7dxHxp5pnAzCAkgsPhqSPet+e3fvDxfqBjiq9mP+I7WlJJyQ9FxGLLrtje6/tvRXPCKAPRVFHxCcRMS5pVNJ6299Z5D6TEbEuItZVPCOAPvR19jsi5iVNSdpUxzAAhldy9vti2xf1vv6apO9LOlTzXAAGVHL2+xJJD9oe0cL/BH4TEU/VOxaAQZWc/f6LFtakBrAM8I4yIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpJZ9svubN++ve0RajM7O9vYti688MLGtpV5qaQuYE8NJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyxVH3Luj/mm0uOgh0WD976jskzdQ1CIBqlC67MyrpOkk76h0HwLBK99TbJd0l6dMvugNraQHdULJCx/WSTkTEvv93P9bSArqhZE+9QdINtmclPSJpo+2Hap0KwMCWjDoi7omI0YgYk7RF0vMRcVPtkwEYCK9TA8n0dTmjiJjSwlK2ADqKPTWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQjCOi+l9qV/9Lv4QmJiYa29bOnTsb29aLL77Y2LaaXuJnfn6+sW1FhBe7nT01kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJFF3OqHcl0Q8lfSLpFJcBBrqrn2uUfS8i3qttEgCV4PAbSKY06pD0e9v7bG9b7A4suwN0Q+nh94aIOG77m5Kes30oIl468w4RMSlpUuKjl0CbivbUEXG8988Tkh6XtL7OoQAMrmSBvHNtn3/6a0k/lPRG3YMBGEzJ4fe3JD1u+/T9fx0Rz9Q6FYCBLRl1RByR9N0GZgFQAV7SApIhaiAZogaSIWogGaIGkiFqIBmiBpLp56OXaNjs7Gxj2zp69Ghj2xofH0+5LUmamppqdHuLYU8NJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRVHbvsj2HtuHbM/YvrLuwQAMpvS937+U9ExE/Nj2Cknn1DgTgCEsGbXtCyRdJWlCkiLipKST9Y4FYFAlh9+XSXpX0k7br9ne0bv+92ew7A7QDSVRnyXpCkn3RcRaSR9Luvvzd4qIyYhYxzK3QLtKop6TNBcRr/S+36OFyAF00JJRR8Q7ko7ZXt276RpJB2udCsDASs9+3y5pd+/M9xFJt9Y3EoBhFEUdEdOSeK4MLAO8owxIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZFhLq8OaXEurSU2uN9WFta2axp4aSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkhmyahtr7Y9fcafD2zf2cBsAAaw5NtEI+JNSeOSZHtE0j8kPV7vWAAG1e/h9zWS/h4RR+sYBsDw+v1AxxZJDy/2A9vbJG0beiIAQyneU/eu+X2DpN8u9nOW3QG6oZ/D72sl7Y+If9Y1DIDh9RP1Vn3BoTeA7iiK2vY5kn4g6bF6xwEwrNJld/4t6es1zwKgAryjDEiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkHBHV/1L7XUn9fjzzG5Leq3yYbsj62Hhc7fl2RFy82A9qiXoQtvdm/YRX1sfG4+omDr+BZIgaSKZLUU+2PUCNsj42HlcHdeY5NYBqdGlPDaACRA0k04mobW+y/abtw7bvbnueKtheZfsF2zO2D9i+o+2ZqmR7xPZrtp9qe5Yq2b7I9h7bh3p/d1e2PVO/Wn9O3Vsg4G9auFzSnKRXJW2NiIOtDjYk25dIuiQi9ts+X9I+SZuX++M6zfZPJa2TdEFEXN/2PFWx/aCkP0TEjt4VdM+JiPmWx+pLF/bU6yUdjogjEXFS0iOSbmx5pqFFxNsRsb/39YeSZiStbHeqatgelXSdpB1tz1Il2xdIukrS/ZIUESeXW9BSN6JeKenYGd/PKcl//KfZHpO0VtIrLY9Sle2S7pL0actzVO0ySe9K2tl7arHD9rltD9WvLkTtRW5L8zqb7fMkPSrpzoj4oO15hmX7ekknImJf27PU4CxJV0i6LyLWSvpY0rI7x9OFqOckrTrj+1FJx1uapVK2z9ZC0LsjIsvllTdIusH2rBaeKm20/VC7I1VmTtJcRJw+otqjhciXlS5E/aqky21f2jsxsUXSky3PNDTb1sJzs5mIuLfteaoSEfdExGhEjGnh7+r5iLip5bEqERHvSDpme3XvpmskLbsTm/0ukFe5iDhl+zZJz0oakfRARBxoeawqbJB0s6S/2p7u3fbziHi6vZFQ4HZJu3s7mCOSbm15nr61/pIWgGp14fAbQIWIGkiGqIFkiBpIhqiBZIgaSIaogWT+C697jzONoTA9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index: 177\n",
      "Correct label: 3\n",
      "Incorrect prediction: 7\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAK1ElEQVR4nO3d24td9RnG8efpRGnigUBri2ZCJ4IMSKETCQEJSBrbEqtoLnqRgEKkkCslQwuivUr/AbUXRRiiRjBV2nhAxGoFDVZorUmctsaJJQ27ZBptlDJ4qCRE317MCkQ7dtbee53mzfcDwTkxv3eTfF1r9qy9fo4IAcjjK20PAKBaRA0kQ9RAMkQNJEPUQDLL6vimtlM+pb58+fJG17viiisaW2vFihWNrXXq1KnG1jp+/Hhja0nSJ5980thaEeGFPl5L1FmNj483ut6uXbsaW2tiYqKxtXq9XmNrTU5ONraWJE1PTze63kI4/QaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkikVte3Ntt+2fdT23XUPBWBwi0Zte0TSLyXdIOlqSdtsX133YAAGU+ZIvV7S0Yg4FhGnJT0u6ZZ6xwIwqDJRr5J07ktdZouPfY7tHbYP2D5Q1XAA+lfmVVoLvbzrf15aGRFTkqakvC+9BJaCMkfqWUmrz3l/VNKJesYBMKwyUb8u6Srba2xfKGmrpGfqHQvAoBY9/Y6IM7bvkPSCpBFJD0XE4donAzCQUnc+iYjnJD1X8ywAKsAVZUAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyrmPT+Sav/R4bG2tqKT399NONrSU1u9vDnj17Gltr+/btja3V5L8PSdq4cWNja33ZtjscqYFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSKbMDh0P2T5p+80mBgIwnDJH6j2SNtc8B4CKLBp1RLwi6d8NzAKgAqXuJlqG7R2SdlT1/QAMprKo2XYH6Aae/QaSIWogmTK/0npM0h8kjduetf3j+scCMKgye2lta2IQANXg9BtIhqiBZIgaSIaogWSIGkiGqIFkiBpIprJrv88HvV6v0fV27drV2FpNPrYtW7Y0tlaTWxd1BUdqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSKXOPstW2X7Y9Y/uw7Z1NDAZgMGWu/T4j6acRccj2JZIO2n4xIt6qeTYAAyiz7c47EXGoePtDSTOSVtU9GIDB9PUqLdtjktZKem2Bz7HtDtABpaO2fbGkJyRNRsQHX/w82+4A3VDq2W/bF2g+6L0R8WS9IwEYRplnvy3pQUkzEXFv/SMBGEaZI/UGSbdJ2mR7uvjzw5rnAjCgMtvuvCrJDcwCoAJcUQYkQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMo6o/rUXvKCjGhMTE42tNTk52dhaTWr6cc3NzTW2VkQseFEYR2ogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIJkyNx78qu0/2f5zse3Oz5sYDMBgytz3+5SkTRHxUXGr4Fdt/zYi/ljzbAAGUObGgyHpo+LdC4o/XNsNdFTZm/mP2J6WdFLSixGx4LY7tg/YPlDxjAD6UCrqiPg0IiYkjUpab/vbC3zNVESsi4h1Fc8IoA99PfsdEXOS9kvaXMcwAIZX5tnvy2yvLN5eLul7ko7UPBeAAZV59vtySY/YHtH8/wR+HRHP1jsWgEGVefb7L5rfkxrAEsAVZUAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kU+aKMhSa3sLlvvvua3S9pqxZs6axtZrcBqcrOFIDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZBM6aiLG/q/YZubDgId1s+ReqekmboGAVCNstvujEq6UdLuescBMKyyR+r7Jd0l6bMv+wL20gK6ocwOHTdJOhkRB//f17GXFtANZY7UGyTdbLsn6XFJm2w/WutUAAa2aNQRcU9EjEbEmKStkl6KiFtrnwzAQPg9NZBMX7czioj9mt/KFkBHcaQGkiFqIBmiBpIhaiAZogaSIWogGaIGkmHbnT7s37+/0fVsN7ZWr9drbK2VK1c2ttb5iCM1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJlLpMtLiT6IeSPpV0htsAA93Vz7Xf342I92ubBEAlOP0GkikbdUj6ne2Dtncs9AVsuwN0Q9nT7w0RccL2NyS9aPtIRLxy7hdExJSkKUmyHRXPCaCkUkfqiDhR/PekpKckra9zKACDK7NB3kW2Lzn7tqQfSHqz7sEADKbM6fc3JT1V3IVjmaRfRcTztU4FYGCLRh0RxyR9p4FZAFSAX2kByRA1kAxRA8kQNZAMUQPJEDWQDFEDybDtTh+mp6cbXa/J7Wnm5uZSrnU+4kgNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAypaK2vdL2PttHbM/YvrbuwQAMpuy137+Q9HxE/Mj2hZJW1DgTgCEsGrXtSyVdJ2m7JEXEaUmn6x0LwKDKnH5fKek9SQ/bfsP27uL+35/DtjtAN5SJepmkayQ9EBFrJX0s6e4vflFETEXEOra5BdpVJupZSbMR8Vrx/j7NRw6ggxaNOiLelXTc9njxoeslvVXrVAAGVvbZ7zsl7S2e+T4m6fb6RgIwjFJRR8S0JH5WBpYArigDkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIJklv5fWxMREY2tt3LixsbUkaXJysrG1du7c2dhavV6vsbXORxypgWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkFo3a9rjt6XP+fGB7soHZAAxg0ctEI+JtSROSZHtE0j8lPVXvWAAG1e/p9/WS/h4R/6hjGADD6/cFHVslPbbQJ2zvkLRj6IkADKX0kbq45/fNkn6z0OfZdgfohn5Ov2+QdCgi/lXXMACG10/U2/Qlp94AuqNU1LZXSPq+pCfrHQfAsMpuu/MfSV+reRYAFeCKMiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaScURU/03t9yT1+/LMr0t6v/JhuiHrY+NxtedbEXHZQp+oJepB2D6Q9RVeWR8bj6ubOP0GkiFqIJkuRT3V9gA1yvrYeFwd1JmfqQFUo0tHagAVIGogmU5EbXuz7bdtH7V9d9vzVMH2atsv256xfdj2zrZnqpLtEdtv2H627VmqZHul7X22jxR/d9e2PVO/Wv+Zutgg4G+av13SrKTXJW2LiLdaHWxIti+XdHlEHLJ9iaSDkrYs9cd1lu2fSFon6dKIuKnteapi+xFJv4+I3cUddFdExFzLY/WlC0fq9ZKORsSxiDgt6XFJt7Q809Ai4p2IOFS8/aGkGUmr2p2qGrZHJd0oaXfbs1TJ9qWSrpP0oCRFxOmlFrTUjahXSTp+zvuzSvKP/yzbY5LWSnqt5VGqcr+kuyR91vIcVbtS0nuSHi5+tNht+6K2h+pXF6L2Ah9L83s22xdLekLSZER80PY8w7J9k6STEXGw7VlqsEzSNZIeiIi1kj6WtOSe4+lC1LOSVp/z/qikEy3NUinbF2g+6L0RkeX2yhsk3Wy7p/kflTbZfrTdkSozK2k2Is6eUe3TfORLSheifl3SVbbXFE9MbJX0TMszDc22Nf+z2UxE3Nv2PFWJiHsiYjQixjT/d/VSRNza8liViIh3JR23PV586HpJS+6JzX43yKtcRJyxfYekFySNSHooIg63PFYVNki6TdJfbU8XH/tZRDzX3kgo4U5Je4sDzDFJt7c8T99a/5UWgGp14fQbQIWIGkiGqIFkiBpIhqiBZIgaSIaogWT+C9wikF6nIA4GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "best_clf = classifiers[best_i]\n",
    "best_params = params[best_i]\n",
    "\n",
    "clf = GridSearchCV(best_clf, best_params, cv=5)\n",
    "clf.fit(X_train, y_train)\n",
    "test_acc = clf.score(X_test, y_test)\n",
    "print(\"Test accuracy:\", test_acc)\n",
    "\n",
    "\n",
    "y_pred = clf.predict(X_test)\n",
    "\n",
    "\n",
    "conf_mat = confusion_matrix(y_test, y_pred)\n",
    "print(\"Confusion matrix:\\n\", conf_mat)\n",
    "\n",
    "\n",
    "for i in range(len(y_test)):\n",
    "    if y_test[i] != y_pred[i]:\n",
    "        print(\"Index:\", i)\n",
    "        print(\"Correct label:\", y_test[i])\n",
    "        print(\"Incorrect prediction:\", y_pred[i])\n",
    "        plt.imshow(X_test[i].reshape(8, 8), cmap='gray')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "442c16b1",
   "metadata": {},
   "source": [
    "## 2. One-class classification (outlier detection)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54a8dbf8",
   "metadata": {},
   "source": [
    "### 2(a) There is an old gradebook at [http://pages.stat.wisc.edu/~jgillett/451/data/midtermGrades.txt](http://pages.stat.wisc.edu/~jgillett/451/data/midtermGrades.txt).\n",
    "Use `pd.read_table()` to read it into a DataFrame.\n",
    "\n",
    "Hint: `pd.read_table()` has many parameters. Check its documentation to find three parameters to:\n",
    "- Read from the given URL\n",
    "- Use the separator '\\s+', which means 'one or more whitespace characters'\n",
    "- Skip the first 12 rows, as they are a note to students and not part of the gradebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb5daf79",
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"http://pages.stat.wisc.edu/~jgillett/451/data/midtermGrades.txt\"\n",
    "df = pd.read_table(url, sep='\\s+', skiprows=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54e6f536",
   "metadata": {},
   "source": [
    "### 2(b) Use `clf = mixture.GaussianMixture(n_components=1)` to make a one-class Gaussian model to decide which $\\mathbf{x}=(\\text{Exam1}, \\text{Exam2})$ are outliers:\n",
    "\n",
    "- Set a matrix X to the first two columns, Exam1 and Exam.\n",
    "- These exams were worth 125 points each. Transform scores to percentages in $[0, 100]$.\n",
    "\n",
    "  Hint: I tried the MinMaxScaler() first, but it does the wrong thing if there aren't scores\n",
    "  of 0 and 125 in each column. So, instead, I just multiplied the whole matrix by 100 / 125.\n",
    "- Fit your classifier to X.\n",
    "  \n",
    "  Hint:\n",
    "  - The reference page for `mixture.GaussianMixture` includes a `fit(X, y=None)` method\n",
    "    with the comment that y is ignored (as this is an unsupervised learning algorithm--there\n",
    "    is no $y$) but present for API consistency. So we can fit with just X.\n",
    "  - I got a warning about \"KMeans ... memory leak\". You may ignore this\n",
    "    warning if you see it. I still got satisfactory results.\n",
    "- Print the center $\\mathbf{\\mu}$ and covariance matrix $\\mathbf{\\Sigma}$ from the two-variable\n",
    "  $N_2(\\mathbf{\\mu}, \\mathbf{\\Sigma})$ distribution you estimated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5072d5b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df.iloc[:,:2] * 100 / 125\n",
    "\n",
    "clf = mixture.GaussianMixture(n_components=1)\n",
    "clf.fit(X)\n",
    "\n",
    "print(\"Mean:\\n\", clf.means_)\n",
    "print(\"\\nCovariance Matrix:\\n\", clf.covariances_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78955408",
   "metadata": {},
   "source": [
    "### 2(c) Here I have given you code to make a contour plot of the negative log likelihood $-\\ln f_{\\mathbf{\\mu}, \\mathbf{\\Sigma}}(\\mathbf{x})$ for $\\mathbf{X} \\sim N_2(\\mathbf{\\mu}, \\mathbf{\\Sigma})$, provided you have set `clf`.\n",
    "\n",
    "```\n",
    "# make contour plot of log-likelihood of samples from clf.score_samples()\n",
    "margin = 10\n",
    "x = np.linspace(0 - margin, 100 + margin)\n",
    "y = np.linspace(0 - margin, 100 + margin)\n",
    "grid_x, grid_y = np.meshgrid(x, y)\n",
    "two_column_grid_x_grid_y = np.array([grid_x.ravel(), grid_y.ravel()]).T\n",
    "negative_log_pdf_values = -clf.score_samples(two_column_grid_x_grid_y)\n",
    "grid_z = negative_log_pdf_values\n",
    "grid_z = grid_z.reshape(grid_x.shape)\n",
    "plt.contour(grid_x, grid_y, grid_z, levels=10) # X, Y, Z\n",
    "plt.title('(Exam1, Exam2) pairs')\n",
    "```\n",
    "\n",
    "Paste my code into your code cell below and add more code:\n",
    "- Add black $x$- and $y$- axes. Label them Exam1 and Exam2.\n",
    "- Plot the data points in blue.\n",
    "- Plot $\\mathbf{\\mu}=$ `clf.means_` as a big lime dot.\n",
    "- Overplot (i.e. plot again) in red the 8 outliers determined by a threshold consisting\n",
    "  of the 0.02 quantile of the pdf values $f_{\\mathbf{\\mu}, \\mathbf{\\Sigma}}(\\mathbf{x})$\n",
    "  for each $\\mathbf{x}$ in X.\n",
    "  \n",
    "  Hint: `clf.score_samples(X)` gives log likelihood, so `np.exp(clf.score_samples(X))`\n",
    "  gives the required $f_{\\mathbf{\\mu}, \\mathbf{\\Sigma}}(\\mathbf{x})$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32ad915c",
   "metadata": {},
   "outputs": [],
   "source": [
    "margin = 10\n",
    "x = np.linspace(0 - margin, 100 + margin)\n",
    "y = np.linspace(0 - margin, 100 + margin)\n",
    "grid_x, grid_y = np.meshgrid(x, y)\n",
    "two_column_grid_x_grid_y = np.array([grid_x.ravel(), grid_y.ravel()]).T\n",
    "negative_log_pdf_values = -clf.score_samples(two_column_grid_x_grid_y)\n",
    "grid_z = negative_log_pdf_values\n",
    "grid_z = grid_z.reshape(grid_x.shape)\n",
    "plt.contour(grid_x, grid_y, grid_z, levels=10)\n",
    "plt.xlabel('Exam1')\n",
    "plt.ylabel('Exam2')\n",
    "plt.title('(Exam1, Exam2) pairs')\n",
    "\n",
    "plt.scatter(X['Exam1'], X['Exam2'], color='blue')\n",
    "plt.scatter(clf.means_[:, 0], clf.means_[:, 1], color='lime', s=150)\n",
    "\n",
    "threshold = np.quantile(np.exp(clf.score_samples(X)), 0.02)\n",
    "pdf_values = np.exp(clf.score_samples(X))\n",
    "outliers = X[pdf_values < threshold]\n",
    "plt.scatter(outliers['Exam1'], outliers['Exam2'], color='red', s=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31574f59",
   "metadata": {},
   "source": [
    "### What characterizes 7 of these 8 outliers? Write your answer in a markdown cell."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09581fc2",
   "metadata": {},
   "source": [
    "If you print the outliers df, you can see that for 7 of the outliers, at lease one of the exam scores for Exam 1 and 2 is a 0.0. Only Student 357 does not have a zero at least one of the exam grades."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56532eab",
   "metadata": {},
   "source": [
    "### 2(d) Write a little code to report whether, by the 0.02 quantile criterion, $\\mathbf{x}=$ (Exam1=50, Exam2=100) is an outlier.\n",
    "\n",
    "Hint: Compare $f_{\\mathbf{\\mu}, \\mathbf{\\Sigma}}(\\mathbf{x})$ to your threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aaf8764",
   "metadata": {},
   "outputs": [],
   "source": [
    "two_d = pd.DataFrame([[50,100]])\n",
    "threshold = np.quantile(np.exp(clf.score_samples(X)), 0.02)\n",
    "np.exp(clf.score_samples(two_d)) < threshold\n",
    "\n",
    "# Yes these exam scores together would be an outlier"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c849af",
   "metadata": {},
   "source": [
    "## 3. Explore the fact that accuracy can be misleading for imbalanced data.\n",
    "Here I make a fake imbalanced data set by randomly sampling $y$ from a distribution with $P(y = 0) = 0.980$ and $P(y = 1) = 0.020$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1153fda6",
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = make_classification(n_samples=1000, n_features=4, n_classes=2, weights=[0.980, 0.020],\n",
    "                           n_clusters_per_class=1, flip_y=0.01, random_state=0)\n",
    "print(f'np.bincount(y)={np.bincount(y)}; we expect about 980 zeros and 20 ones.')\n",
    "print(f'np.mean(y)={np.mean(y)}; we expect the proportion of ones to be about 0.020.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66eca00d",
   "metadata": {},
   "source": [
    "Here I split the data into 50% training and 50% testing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f04c1e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.5,\n",
    "                                                    random_state=0, stratify=y)\n",
    "print(f'np.bincount(y_train)={np.bincount(y_train)}')\n",
    "print(f'np.mean(y_train)={np.mean(y_train)}.')\n",
    "print(f'np.bincount(y_test)={np.bincount(y_test)}.')\n",
    "print(f'np.mean(y_test)={np.mean(y_test)}.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d57e1786",
   "metadata": {},
   "source": [
    "### 3a. Train and assess a gradient boosting model.\n",
    "- Train on the training data.\n",
    "- Use 100 trees of maximum depth 1 and learning rate $\\alpha = 0.25$.\n",
    "- Use `random_state=0` (to give us all a chance of getting the same results).\n",
    "- Display the accuracy, precision, recall, and AUC on the test data. Use 3 decimal places.\n",
    "  Use a labeled print statement with 3 decimal places so the reader can easily find each metric.\n",
    "- Make an ROC curve from your classifier and the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e12227b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = GradientBoostingClassifier(n_estimators=100, max_depth=1, learning_rate=0.25, random_state=0)\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "y_pred = clf.predict(X_test)\n",
    "acc = accuracy_score(y_test, y_pred)\n",
    "prec = precision_score(y_test, y_pred)\n",
    "rec = recall_score(y_test, y_pred)\n",
    "auc = roc_auc_score(y_test, clf.predict_proba(X_test)[:,1])\n",
    "print(f\"Accuracy: {acc:.3f}, Precision: {prec:.3f}, Recall: {rec:.3f}, AUC: {auc:.3f}\")\n",
    "\n",
    "RocCurveDisplay.from_estimator(clf, X_test, y_test)\n",
    "plt.plot([0, 1], [0, 1], 'k--', label='Random Guess')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('ROC Curve')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e67008",
   "metadata": {},
   "source": [
    "Note the high accuracy but lousy precision, recall, and AUC.\n",
    "\n",
    "Note that since the data have about 98% $y = 0$, we could get about 98% accuracy\n",
    "by just always predicting $\\hat{y} = 0$. High accuracy alone is not necessarily helpful."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16399469",
   "metadata": {},
   "source": [
    "### 3b. Now oversample the data to get a balanced data set.\n",
    "- Use the `RandomOverSampler(random_state=0)` to oversample and get a balanced data set.\n",
    "- Repeat my `train_test_split()` block from above.\n",
    "- Repeat your train/assess block from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e14447b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = RandomOverSampler(random_state=0).fit_resample(X, y)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.5,\n",
    "                                                    random_state=0, stratify=y)\n",
    "\n",
    "#from last cell\n",
    "clf = GradientBoostingClassifier(n_estimators=100, max_depth=1, learning_rate=0.25, random_state=0)\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "y_pred = clf.predict(X_test)\n",
    "acc = accuracy_score(y_test, y_pred)\n",
    "prec = precision_score(y_test, y_pred)\n",
    "rec = recall_score(y_test, y_pred)\n",
    "auc = roc_auc_score(y_test, clf.predict_proba(X_test)[:,1])\n",
    "print(f\"Accuracy: {acc:.3f}, Precision: {prec:.3f}, Recall: {rec:.3f}, AUC: {auc:.3f}\")\n",
    "\n",
    "RocCurveDisplay.from_estimator(clf, X_test, y_test)\n",
    "plt.plot([0, 1], [0, 1], 'k--', label='Random Guess')\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('ROC Curve')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d218737",
   "metadata": {},
   "source": [
    "Note that we traded a little accuracy for much improved precision, recall, and AUC.\n",
    "\n",
    "If you do classification in your project and report accuracy, please\n",
    "also report the proportions of $y = 0$ and $y = 1$ in your test data so that\n",
    "we get insight into whether your model improves upon always guessing $\\hat{y} = 0$\n",
    "or always guessing $\\hat{y} = 1$."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}