lectures/Lecture_19.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Classification: I\n",
    "\n",
    "*S. R. Taylor (2021)*\n",
    "\n",
    "This lecture and notebook are based on the \"Classification\" and \"Classification2\" lectures of of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), which in turn are based on materials from Andy Connolly, and Ivezic et al. Chapter 9, Andy Connolly's [blog](http://connolly.github.io/introAstroML/blog/regression.html), and Aurelien Geron's [book](https://github.com/ageron/handson-ml2). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 9.\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [Generative vs Discriminative Classification](#one)\n",
    "* [Comparing the performance of classifiers](#two)\n",
    "* [Generative classifiers](#three)\n",
    "* [Simplest classifier: naive Bayes](#four)\n",
    "* [Gaussian naive Bayes](#five)\n",
    "* [Linear & Quadratic discriminant analysis](#six)\n",
    "* [GMM & Bayes classification](#seven)\n",
    "* [$K$-nearest neighbor classifier](#eight)\n",
    "    \n",
    "---\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>. Points will be given this week mostly for producing and understanding the notebook plots. But there are also discussions to be had amongst your breakout rooms.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Several weeks ago we looked at density estimation and clustering, which are **unsupervised** forms of classification. Let's now look at **supervised** classification. This is where we actually know the \"truth\" for some of our objects and can use that information to help guide the classification of unknown objects.\n",
    "\n",
    "\n",
    "## Generative vs. Discriminative Classification <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "We will talk about two different types of classification where each has a slightly different approach.  As an example, if you are trying to determine whether your neighbor is speaking Spanish or Portuguese, you could 1) learn both Spanish and Portuguese so that you'd know exactly what they are saying or 2) learn the keys rules about the differences between the languages.\n",
    "\n",
    "- If we find ourselves asking which category is most likely to generate the observed result, then we are using using **density estimation** for classification and this is referred to as **generative classification**. Here we have a full model of the density for each class or we have a model which describes how data could be generated from each class. \n",
    "\n",
    "\n",
    "- But if we don't care about the full distribution, then we are doing something more like clustering, where we don't need to map the distribution, we just need to define *boundaries*.  Classification that finds the **decision boundary** that separates classes is called **discriminative classification**.  For high-dimensional data, this may be a better choice.\n",
    "\n",
    "For example, in the figure below, to classify a new object with $x=1$, it would suffice to know that either \n",
    "1. model 1 is a better fit than model 2 (***generative classification***), or \n",
    "2. that the decision boundary is at $x=1.4$ (***discriminative classification***).\n",
    "\n",
    "![Ivezic, Figure 9.1](http://www.astroml.org/_images/fig_bayes_DB_1.png)\n",
    "\n",
    "\n",
    "### Assessing Your Results\n",
    "\n",
    "Before even introducing different schemes, let's talk about how we define the success of our classification. We'll bring back some concepts from the start of the course. \n",
    "\n",
    "Let's first consider **binary classification** where each observation is assigned to either **class=+1** or **class=-1** (or 0 depending on your preference). Let's ignore for now the fact that each observation can be assigned a probability of belonging to class 1 or -1 and we'll only allow those two possibilities (e.g., in the figure above $x=2$ is class +1 and $x=-0.5$ is class -1). \n",
    "\n",
    "In that case, there are 2 types of errors (where our goal is to identify class +1 objects):\n",
    "* a **[False Positive](https://en.wikipedia.org/wiki/False_positives_and_false_negatives#False_positive_error)**, where we have assigned a *true* class label when it is really false. ---> ***Type-1 error***.\n",
    "* a **[False Negative](https://en.wikipedia.org/wiki/False_positives_and_false_negatives#False_positive_error)**, where we have assigned a *false* class label when it is really true. ---> ***Type-II error***.\n",
    "\n",
    "All 4 [possibilities](https://en.wikipedia.org/wiki/Sensitivity_and_specificity) are (if you want apples [$g_2$], but not oranges [$g_1$]):\n",
    "- **True Positive** = ***correctly identified***  (apple identified as apple)\n",
    "- **True Negative** = ***correctly rejected***  (orange rejected as orange)\n",
    "- **False Positive** = ***incorrectly identified***  (orange identified as apple)\n",
    "- **False Negative** = ***incorrectly rejected***  (apple rejected as orange)\n",
    "\n",
    "Based on these, we usually define either of the following pairs of terms.  Which is used is largely a matter of preference in different fields, but we'll see that there are some key differences.\n",
    "\n",
    "$$ \n",
    "{\\rm completeness} = \\frac{\\rm true\\ positives}{\\rm true\\ positives + false\\ negatives} = {\\rm true\\ positive\\ rate\\ /\\ sensitivity\\ /\\ recall}\n",
    "$$\n",
    "\n",
    "$$  \n",
    "{\\rm contamination} = \\frac{\\rm false\\ positives}{\\rm true\\ positives + false\\ positives} = {\\rm false\\ discovery\\ rate}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$ \n",
    "{\\rm true\\ positive\\ rate} = \\frac{\\rm true\\ positives} {\\rm true\\ positives + false\\ negatives}\n",
    "$$\n",
    "\n",
    "$$  \n",
    "{\\rm false\\ positive\\ rate} = \\frac{\\rm false\\ positives} {\\rm true\\ negatives + false\\ positives} = {\\rm Type1\\ error}\n",
    "$$\n",
    "\n",
    "Similarly \n",
    " \n",
    "$$\n",
    "{\\rm efficiency} = 1 - {\\rm contamination} = {\\rm precision}. \n",
    "$$\n",
    "\n",
    "Depending on your goals, you may want to maximize the completeness or the efficiency, or a combination of both.\n",
    "\n",
    "To better understand the differences between these measures, let's take the needle in a haystack problem.  You have 100,000 stars and 1000 quasars.  If you correctly identify 900 quasars as such and mistake 1000 stars for quasars, then we have:\n",
    "- TP = 900\n",
    "- FN  = 100\n",
    "- TN = 99,000\n",
    "- FP = 1000\n",
    "\n",
    "Which gives\n",
    "\n",
    "$$\n",
    "{\\rm true\\ positive\\ rate} = \\frac{900}{900 + 100} = 0.9 = {\\rm completeness}\n",
    "$$\n",
    "\n",
    "$$\n",
    "{\\rm false\\ positive\\ rate} = \\frac{1000}{99000 + 1000} = 0.01\n",
    "$$\n",
    "\n",
    "Not bad right?  Well, sort of.  The FPR isn't bad, but there are *lots* of stars, so the contamination rate isn't so great. <font color='red'>Compute and write this below.</font>\n",
    "\n",
    "<!--- \n",
    "$$\n",
    "{\\rm contamination} = \\frac{1000}{900 + 1000} = 0.53\n",
    "$$\n",
    "-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparing the performance of classifiers <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "\"Best\" performance is subjective; we trade contamination versus completeness depending on the science. Before we even talk about different classification schemes, let's first talk about how we can quantify which of the methods is \"best\".  (N.B.  We have skipped ahead to Ivezic $\\S$ 9.9).  \n",
    "\n",
    "We do this with a [**Receiver Operating Characteristic (ROC)**](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) curve. **A ROC curve simply plots the true-positive vs. the false-positive rate.** \n",
    "\n",
    "One concern about ROC curves is that they are sensitive to the relative sample sizes. As we already demonstrated above, **if there are many more background events than source events, small false positive results can dominate a signal.** For these cases ***we can plot completeness versus efficiency (1 - contamination).***\n",
    "\n",
    "Here is a comparison of the two types of plots:\n",
    "\n",
    "![Ivezic, Figure 9.17](http://www.astroml.org/_images/fig_ROC_curve_1.png)\n",
    "\n",
    "Here we see that to get higher completeness, you could actually suffer significantly in terms of efficiency, but your FPR might not go up that much if there are lots of true negatives.\n",
    "\n",
    "Note that you ***choose*** the completeness and efficiency that you want by choosing a **threshold (decision boundary)**. The curves show you what your possible choices are (depending on where you set the threshold). Generally, you want to chose a decision boundary that maximizes the area under the ROC curve.\n",
    "\n",
    "Below is the code that makes these plots.  We'll talk about the data that goes into it in a bit.  For now, we'll concentrate on how to generate the ROC and completeness-contamination plots. We'll be comparing 7 different classifiers (with a generic `clf` object), making training and test sets with `split_samples`, then using these tools to generate our plots:\n",
    "\n",
    "- [sklearn.metrics.roc_curve(y_test, y_prob)](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html)\n",
    "- [sklearn.metrics.precision_recall_curve(y_test, y_prob)](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.precision_recall_curve.html)\n",
    "- astroML.utils.completeness_contamination(y_pred, y_test)\n",
    "\n",
    "Note that the [`sklearn.metrics` algorithms](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics) take `y_test`, which are classes, and `y_prob`, which are not class predictions, but rather probabilities, whereas the AstroML algorithm wants `y_pred` (which we get by converting `y_prob` into discrete predictions as a function of the probability).\n",
    "\n",
    "**NOTE!!** `sklearn-0.23.2` works with the following code. Newer versions may break. If so, revert to 0.23.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 9.17, edits by GTR\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.discriminant_analysis import (LinearDiscriminantAnalysis,\n",
    "                                           QuadraticDiscriminantAnalysis)\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from astroML.classification import GMMBayes\n",
    "\n",
    "from sklearn.metrics import roc_curve\n",
    "from sklearn.metrics import precision_recall_curve\n",
    "\n",
    "from astroML.utils import split_samples, completeness_contamination\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined()\n",
    "y = y.astype(int)\n",
    "(X_train, X_test), (y_train, y_test) = split_samples(X, y, [0.75, 0.25],\n",
    "                                                     random_state=0)\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Fit all the models to the training data\n",
    "def compute_models(*args):\n",
    "    names = []\n",
    "    probs = []\n",
    "    for classifier, kwargs in args:\n",
    "        print(classifier.__name__)\n",
    "        clf = classifier(**kwargs)\n",
    "        clf.fit(X_train, y_train)\n",
    "        \n",
    "        #Note that we are outputing the probabilities [of class 1], not the classes\n",
    "        y_probs = clf.predict_proba(X_test)[:, 1]\n",
    "\n",
    "        names.append(classifier.__name__)\n",
    "        probs.append(y_probs)\n",
    "\n",
    "    return names, probs\n",
    "\n",
    "\n",
    "names, probs = compute_models((GaussianNB, {}),\n",
    "                              (LinearDiscriminantAnalysis, {}),\n",
    "                              (QuadraticDiscriminantAnalysis, {}),\n",
    "                              (LogisticRegression,\n",
    "                               dict(class_weight='balanced')),\n",
    "                              (KNeighborsClassifier,\n",
    "                               dict(n_neighbors=10)),\n",
    "                              (DecisionTreeClassifier,\n",
    "                               dict(random_state=0, max_depth=12,\n",
    "                                    criterion='entropy')),\n",
    "                              (GMMBayes, dict(n_components=3, tol=1E-5,\n",
    "                                              covariance_type='full')))\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot ROC curves and completeness/efficiency\n",
    "fig = plt.figure(figsize=(15, 5))\n",
    "fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15, top=0.9, wspace=0.25)\n",
    "\n",
    "# ax1 will show roc curves\n",
    "ax1 = plt.subplot(131)\n",
    "\n",
    "# ax2 will show completeness/efficiency\n",
    "ax2 = plt.subplot(132)\n",
    "\n",
    "# ax3 will show precision/recall\n",
    "ax3 = plt.subplot(133)\n",
    "\n",
    "labels = dict(GaussianNB='GNB',\n",
    "              LinearDiscriminantAnalysis='LDA',\n",
    "              QuadraticDiscriminantAnalysis='QDA',\n",
    "              KNeighborsClassifier='KNN',\n",
    "              DecisionTreeClassifier='DT',\n",
    "              GMMBayes='GMMB',\n",
    "              LogisticRegression='LR')\n",
    "\n",
    "thresholds = np.linspace(0, 1, 1001)[:-1]\n",
    "\n",
    "# iterate through and show results\n",
    "for name, y_prob in zip(names, probs):\n",
    "    fpr, tpr, thresh = roc_curve(y_test, y_prob)\n",
    "    precision, recall, thresh2 = precision_recall_curve(y_test, y_prob)\n",
    "\n",
    "    # add (0, 0) as first point\n",
    "    fpr = np.concatenate([[0], fpr])\n",
    "    tpr = np.concatenate([[0], tpr])\n",
    "    # Here we add (1,0) \n",
    "    precision = np.concatenate([[0], precision])\n",
    "    recall = np.concatenate([[1], recall])\n",
    "    thresh2 = np.concatenate([[0], thresh2])\n",
    "\n",
    "    ax1.plot(fpr, tpr, label=labels[name])\n",
    "\n",
    "    #See note above about astroML vs. sklearn\n",
    "    #Note that the range of threshhold values here is 0% to 100% (0.0 to 1.0)\n",
    "    comp = np.zeros_like(thresholds)\n",
    "    cont = np.zeros_like(thresholds)\n",
    "    for i, t in enumerate(thresholds):\n",
    "        y_pred = (y_prob >= t)\n",
    "        comp[i], cont[i] = completeness_contamination(y_pred, y_test)\n",
    "    ax2.plot(1 - cont, comp, label=labels[name])\n",
    "    \n",
    "    ax3.plot(precision, recall, label=labels[name])\n",
    "\n",
    "ax1.set_xlim(0, 0.04)\n",
    "ax1.set_ylim(0, 1.02)\n",
    "ax1.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "ax1.set_xlabel('false positive rate')\n",
    "ax1.set_ylabel('true positive rate')\n",
    "ax1.legend(loc=4)\n",
    "\n",
    "ax2.set_xlabel('efficiency')\n",
    "ax2.set_ylabel('completeness')\n",
    "ax2.set_xlim(0, 1.0)\n",
    "ax2.set_ylim(0.2, 1.02)\n",
    "\n",
    "ax3.set_xlabel('precision')\n",
    "ax3.set_ylabel('recall')\n",
    "ax3.set_xlim(0, 1.0)\n",
    "ax3.set_ylim(0.2, 1.02)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that I've plotted both recall-precision and completeness-efficiency just to demonstrate that they are the same thing.\n",
    "\n",
    "The plot below shows the values of precision and recall as a function of the threshold value, highlighting the value where precision is 90%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#From Geron\n",
    "#Uses the values of precisions, recalls, and thresholds from above\n",
    "def plot_precision_recall_vs_threshold(precisions, recalls, thresholds):\n",
    "    plt.plot(thresholds, precisions[:-1], \"b--\", label=\"Precision\", linewidth=2)\n",
    "    plt.plot(thresholds, recalls[:-1], \"g-\", label=\"Recall\", linewidth=2)\n",
    "    plt.legend(loc=\"best\", fontsize=16)\n",
    "    plt.xlabel(\"Threshold\", fontsize=16)\n",
    "    plt.grid(True)                \n",
    "    plt.axis([0, 1, 0, 1])             \n",
    "\n",
    "recall_90_precision = recall[np.argmax(precision >= 0.90)]\n",
    "threshold_90_precision = thresh2[np.argmax(precision >= 0.90)]\n",
    "\n",
    "plt.figure(figsize=(8, 4))                                                                  \n",
    "plot_precision_recall_vs_threshold(precision, recall, thresh2)\n",
    "plt.plot([threshold_90_precision, threshold_90_precision], [0., 0.9], \"r:\")                 \n",
    "plt.plot([0, threshold_90_precision], [0.9, 0.9], \"r:\")                                \n",
    "plt.plot([0, threshold_90_precision], [recall_90_precision, recall_90_precision], \"r:\")\n",
    "plt.plot([threshold_90_precision], [0.9], \"ro\")                                             \n",
    "plt.plot([threshold_90_precision], [recall_90_precision], \"ro\")                             \n",
    "#save_fig(\"precision_recall_vs_threshold_plot\")                                              \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we'll look at the same in terms of completeness and contamination.  While these are the same as recall and precision, their calculation is slightly different (using either sklearn or astroML)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eff = 1-cont\n",
    "def plot_comp_eff_vs_threshold(eff, comp, thresholds):\n",
    "    plt.plot(thresholds, eff, \"b--\", label=\"Efficiency\", linewidth=2)\n",
    "    plt.plot(thresholds, comp, \"g-\", label=\"Completeness\", linewidth=2)\n",
    "    plt.legend(loc=\"best\", fontsize=16)\n",
    "    plt.xlabel(\"Threshold\", fontsize=16)\n",
    "    plt.grid(True)                \n",
    "    plt.axis([0, 1, 0, 1])             \n",
    "\n",
    "comp_90_eff = comp[np.argmax(eff >= 0.90)]\n",
    "threshold_90_eff = thresholds[np.argmax(eff >= 0.90)]\n",
    "\n",
    "plt.figure(figsize=(8, 4))                                                                  \n",
    "plot_comp_eff_vs_threshold(eff, comp, thresholds)\n",
    "plt.plot([threshold_90_eff, threshold_90_eff], [0., 0.9], \"r:\")                 \n",
    "plt.plot([0, threshold_90_eff], [0.9, 0.9], \"r:\")                                \n",
    "plt.plot([0, threshold_90_eff], [comp_90_eff, comp_90_eff], \"r:\")\n",
    "plt.plot([threshold_90_eff], [0.9], \"ro\")                                             \n",
    "plt.plot([threshold_90_eff], [comp_90_eff], \"ro\")                                                                          \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, here's how to visualize a constraint on precision (in terms of what it means for recall) in our precision-recall plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def plot_precision_vs_recall(precisions, recalls):\n",
    "    plt.plot(recalls, precisions, \"b-\", linewidth=2)\n",
    "    plt.xlabel(\"Recall\", fontsize=16)\n",
    "    plt.ylabel(\"Precision\", fontsize=16)\n",
    "    plt.axis([0, 1, 0, 1])\n",
    "    plt.grid(True)\n",
    "\n",
    "    \n",
    "recall_90_precision = recall[np.argmax(precision >= 0.90)]\n",
    "threshold_90_precision = thresh2[np.argmax(precision >= 0.90)]\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plot_precision_vs_recall(precision, recall)\n",
    "#plt.plot([0.4368, 0.4368], [0., 0.9], \"r:\")\n",
    "plt.plot([recall_90_precision,recall_90_precision], [0., 0.9], \"r:\")\n",
    "plt.plot([0.0, recall_90_precision], [0.9, 0.9], \"r:\")\n",
    "plt.plot(recall_90_precision, [0.9], \"ro\")\n",
    "#save_fig(\"precision_vs_recall_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So far we have only looked at binary classification.  However, you may instead want:\n",
    "\n",
    "* [Multiclass](https://scikit-learn.org/stable/modules/multiclass.html): to distinguish between more than 2 classes (e.g., MNIST digits)\n",
    "\n",
    "* [Multilabel](https://scikit-learn.org/stable/modules/multiclass.html): if you want to allow multiple class labels for each object\n",
    "\n",
    "But those are beyond the scope of what we'll see this week."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Generative Classification <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "With some assessment criteria defined, we can talk about classification itself.  **We can use Bayes' theorem to relate the labels to the features in an $N\\times D$ data set $X$.**  The $j$th feature of the $i$th point is $x^j_i$ and there are $k$ classes giving discrete labels $y_k$.  We have\n",
    "\n",
    "$$p(y_k|x_i) = \\frac{p(x_i|y_k)p(y_k)}{\\sum_i p(x_i|y_k)p(y_k)},$$\n",
    "\n",
    "where $x_i$ is assumed to be a vector with $j$ components.\n",
    "\n",
    "$p(y=y_k)$ is the probability of any point having class $k$ (equivalent to the prior probability of the class $k$). \n",
    "\n",
    "In generative classifiers we model class-conditional densities $p(x|y=y_k)$. \n",
    "\n",
    "Before we get into the generative classification algortithms, we'll first discuss 3 general concepts:\n",
    "- (a) Discriminant Functions\n",
    "- (b) Bayes Classifiers\n",
    "- (c) Decision Boundaries\n",
    "\n",
    "### (a) The Discriminant Function\n",
    "\n",
    "**We can relate classification to density estimation and regression**. $\\hat{y} = f(y|x)$ represents the best guess of $y$ given $x$. **So classification is just regression with discrete $y$ values, e.g., $y=\\{0,1\\}$.**\n",
    "\n",
    "In classification we refer to $f(y|x)$ as the [**discriminant function**](https://en.wikipedia.org/wiki/Discriminant_function_analysis).\n",
    "\n",
    "For a simple 2-class example, where $y=\\{0,1\\}$:\n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "g(x) = f(y|x) & = &  \\int y \\, p(y|x) \\, dy \\\\\n",
    "%    & = & \\int y p(y|x) \\, dy \\\\\n",
    "       & = & 1 \\cdot p(y=1 | x) + 0 \\cdot p(y=0 | x) = p(y=1 | x).\n",
    "%     & = & p(y=1 | x)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "and then using Bayes' rule:\n",
    "\n",
    "$$g(x) = \\frac{p(x|y=1) \\, p(y=1)}{p(x|y=1) \\, p(y=1)  + p(x|y=0) \\, p(y=0)}$$\n",
    "\n",
    "The first equation is just the expectation value of y.\n",
    "\n",
    "### (b) Bayes Classifier\n",
    "\n",
    "If the discriminant function gives a binary prediction, we call it a **Bayes classifier**, formulated as\n",
    "\n",
    "$$\\begin{eqnarray} \\widehat{y} & = & \\left\\{ \\begin{array}{cl}       \t           1 & \\mbox{if $g(x) > 1/2$}, \\\\       \t           0 & \\mbox{otherwise,}       \t           \\end{array}     \t   \\right. \\\\     & = & \\left\\{\n",
    "\\begin{array}{cl}               1 & \\mbox{if $p(y=1|x) > p(y=0|x)$}, \\\\               0 & \\mbox{otherwise.}               \\end{array}       \\right.\\end{eqnarray}$$\n",
    "\n",
    "Where this can be generalized to any number of classes, $k$, and not just two.\n",
    "\n",
    "### (c) Decision Boundary\n",
    "\n",
    "A **decision boundary** is just set of $x$ values at which each class is equally likely:\n",
    "\n",
    "$$\n",
    "p(x|y=1)p(y=1)  =  p(x|y=0)p(y=0);\n",
    "$$\n",
    "\n",
    "$$g_1(x) = g_2(x) \\; {\\rm or}\\; g(x) = 1/2$$\n",
    "\n",
    "Below is an example of a decision boundary in 1-D, where each class is equally likely so we can just look at $p(x)$.  In short, we assign classifications according to which pdf is higher at every given $x$.\n",
    "\n",
    "![Ivezic, Figure 9.1](http://www.astroml.org/_images/fig_bayes_DB_1.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simplest Classifier: Naive Bayes <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "In practice, classification is complicated as the data are generally multi-dimensional (that is we don't just have $x$, we have $x^{j=0},x^1,x^2,x^3...x^N$, so we want $p(x^0,x^1,x^2,x^3...x^N|y)$.\n",
    "\n",
    "Naive Bayes classification comes to the rescue here: note that it's not very clearly in the textbook section that this is for discrete features $\\{x\\}$. If we **assume** that ***all attributes are conditionally independent*** (which is not always true, but is often close enough), then multi-dimensional pdfs simplify to\n",
    "\n",
    "$$ p(x^i,x^j|y_k) = p(x^i|y)p(x^j|y_k)$$\n",
    "  \n",
    "which can be written as\n",
    "\n",
    "$$ p({x^{j=0},x^1,x^2,\\ldots,x^N}|y_k) = \\prod_j p(x^j|y_k).$$\n",
    "\n",
    "From Bayes' rule and conditional independence we get\n",
    "\n",
    "$$\n",
    "  p(y_k | {x^0,x^1,\\ldots,x^N}) =\n",
    "  \\frac{\\prod_j p(x^j|y_k) p(y_k)}\n",
    "       {\\sum_l \\prod_j p(x^j|y_l) p(y_l)}.\n",
    "$$\n",
    "\n",
    "We calculate the most likely value of $y$ by maximizing this over choices of $y_k$:\n",
    "\n",
    "$$\n",
    "\\hat{y} = \\arg \\max_{y_k} \\frac{\\prod_j p(x^j|y_k) p(y_k)}\n",
    "        {\\sum_l \\prod_j p(x^j|y_l) p(y_l)},\n",
    "$$\n",
    "\n",
    "\n",
    "From there the process is just estimating densities: $p(x|y=y_k)$ and $p(y=y_k)$ are learned from a set of training data, where\n",
    "- **$p(y=y_k)$ is just the frequency of the class $k$ in the training set**\n",
    "- **$p(x|y=y_k)$ is just the density (probability) of an object with class $k$ having the attributes $x$**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian Naive Bayes <a class=\"anchor\" id=\"five\"></a>\n",
    "\n",
    "One way to handle continuous values for $X$ is to model $p(x^j|y=y_k)$ as one-dimensional normal distributions, with means $\\mu^j_k$ and widths $\\sigma^j_k$. The naive Bayes estimator is then\n",
    "\n",
    "$$\\hat{y} = \\arg\\max_{y_k}\\left[\\ln p(y=y_k) - \\frac{1}{2}\\sum_{j=1}^N\\left(2\\pi(\\sigma^j_k)^2 + \\frac{(x^j - \\mu^j_k)^2}{(\\sigma^j_k)^2} \\right) \\right]$$\n",
    "\n",
    "In Scikit-Learn [`Gaussian Naive Bayes`](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html) classification is implemented as follows, with a simple example given in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "X = np.random.random((100,2))\n",
    "y = (X[:,0] + X[:,1] > 1).astype(int) #Class 1 if sum of both features is >1\n",
    "\n",
    "gnb = GaussianNB()\n",
    "gnb.fit(X,y)\n",
    "\n",
    "y_pred = gnb.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 9.2, edits by GTR and SRT\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "from matplotlib import colors\n",
    "\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Simulate some data\n",
    "np.random.seed(0)\n",
    "mu1 = [1, 1]\n",
    "cov1 = 0.3 * np.eye(2)\n",
    "\n",
    "mu2 = [5, 3]\n",
    "cov2 = np.eye(2) * np.array([0.4, 0.1])\n",
    "\n",
    "X = np.concatenate([np.random.multivariate_normal(mu1, cov1, 100),\n",
    "                    np.random.multivariate_normal(mu2, cov2, 100)])\n",
    "y = np.zeros(200)\n",
    "y[100:] = 1\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Fit the Naive Bayes classifier\n",
    "clf = GaussianNB()\n",
    "clf.fit(X, y)\n",
    "\n",
    "# predict the classification probabilities on a grid\n",
    "xlim = (-1, 8)\n",
    "ylim = (-1, 5)\n",
    "xx, yy = np.meshgrid(np.linspace(xlim[0], xlim[1], 71),\n",
    "                     np.linspace(ylim[0], ylim[1], 81))\n",
    "Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])\n",
    "# Gives probability for both class 1 and class 2 for each grid point\n",
    "# As these are degenerate, take just one and then\n",
    "# re-shape it to the grid pattern needed for contour plotting\n",
    "Z = Z[:, 1].reshape(xx.shape)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.binary, zorder=2)\n",
    "\n",
    "\n",
    "# Add the decision boundary, which is just the contour where\n",
    "# the probability exceeds some threshold, here 0.5.\n",
    "ax.contour(xx, yy, Z, [0.5], colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**And now an example using real data-- this uses Gaussian Naive Bayes classification to separate RR Lyrae stars from non-variable main squence stars.** Here we have a 4-D $X$ and we are going to take them 1-D at a time to see how much improvement comes from adding each new dimension of $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 9.3, edits by GTR and SRT\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import split_samples\n",
    "from astroML.utils import completeness_contamination\n",
    "# Added by GTR\n",
    "from sklearn.metrics import precision_recall_curve, roc_curve\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined() # X is a 4-D color-color-color-color space\n",
    "X = X[:, [1, 0, 2, 3]]  # rearrange columns for better 1-color results\n",
    "(X_train, X_test), (y_train, y_test) = split_samples(X, y, [0.75, 0.25],\n",
    "                                                     random_state=0)\n",
    "\n",
    "N_tot = len(y)\n",
    "N_st = np.sum(y == 0)\n",
    "N_rr = N_tot - N_st\n",
    "N_train = len(y_train)\n",
    "N_test = len(y_test)\n",
    "N_plot = 5000 + N_rr\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# perform Naive Bayes\n",
    "classifiers = []\n",
    "predictions = []\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "\n",
    "order = np.array([1, 0, 2, 3])\n",
    "\n",
    "y_prob = np.array([])\n",
    "\n",
    "for nc in Ncolors:\n",
    "    clf = GaussianNB()\n",
    "    clf.fit(X_train[:, :nc], y_train)\n",
    "    y_pred = clf.predict(X_test[:, :nc])\n",
    "    \n",
    "    # Added by GTR to be able to compute precision, recall, fpr, and tpr\n",
    "    # Gives the probability for both classes, take just one\n",
    "    y_prob = np.append(y_prob,clf.predict_proba(X_test[:, :nc])[:,1])\n",
    "\n",
    "    classifiers.append(clf)\n",
    "    predictions.append(y_pred)\n",
    "\n",
    "completeness, contamination = completeness_contamination(predictions, y_test)\n",
    "\n",
    "print(\"completeness\", completeness)\n",
    "print(\"contamination\", contamination)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the decision boundary\n",
    "clf = classifiers[1]\n",
    "xlim = (0.7, 1.35)\n",
    "ylim = (-0.15, 0.4)\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xlim[0], xlim[1], 81),\n",
    "                     np.linspace(ylim[0], ylim[1], 71))\n",
    "\n",
    "Z = clf.predict_proba(np.c_[yy.ravel(), xx.ravel()])\n",
    "Z = Z[:, 1].reshape(xx.shape)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(10, 5))\n",
    "fig.subplots_adjust(bottom=0.15, top=0.95, hspace=0.0,\n",
    "                    left=0.1, right=0.95, wspace=0.2)\n",
    "\n",
    "# left plot: data and decision boundary\n",
    "ax = fig.add_subplot(121)\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.Oranges, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "im = ax.imshow(Z, origin='lower', aspect='auto',\n",
    "               cmap=plt.cm.binary, zorder=1,\n",
    "               extent=xlim + ylim)\n",
    "im.set_clim(0, 1.5)\n",
    "ax.contour(xx, yy, Z, [0.5], colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$u-g$')\n",
    "ax.set_ylabel('$g-r$')\n",
    "\n",
    "# Plot completeness vs Ncolors\n",
    "ax = plt.subplot(222)\n",
    "ax.plot(Ncolors, completeness, 'o-k', ms=6)\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "\n",
    "ax.set_ylabel('completeness')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "# Plot contamination vs Ncolors\n",
    "ax = plt.subplot(224)\n",
    "ax.plot(Ncolors, contamination, 'o-k', ms=6)\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%i'))\n",
    "\n",
    "ax.set_xlabel('N colors')\n",
    "ax.set_ylabel('contamination')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- In the **left panel**, peach points represent non-variable sources, while dark grey points are variable sources. The decision boundary is a curved black line, with the classification porbability shown as a shaded background.\n",
    "- In the **right panel**, more and more stellar color features are used to train the classification scheme. \n",
    "\n",
    "<font color='red'>If you shifted the decision boundary \"up\" by hand, discuss with your breakout room what would happen to the completeness, contamination, and false positive rate? What happens if you change the fraction of objects in the training set?</font>\n",
    "\n",
    "The \"*naive*\" in Naive Bayes refers to the fact that we are assuming that all of the variables are conditionally independent. If we relax that assumption and allow for covariances, then we have a **Gaussian Bayes classifier**.  But note that this comes with a large jump in computational cost! Recall the form of a multivariate Gaussian distribution, where $\\vec{x}$ and $\\vec{\\mu}_k$ are $D$-dimensional vectors over data features, and $\\Sigma_k$ is a $D\\times D$ symmetric covariance matrix:\n",
    "\n",
    "$$ p_k(\\vec{x}) = \\frac{1}{\\sqrt{\\mathrm{det}(2\\pi\\Sigma_k)}}\\exp\\left[ -\\frac{1}{2}(\\vec{x}-\\vec{\\mu}_k)^T \\Sigma_k^{-1}(\\vec{x}-\\vec{\\mu}_k) \\right].$$\n",
    "\n",
    "**The Gaussian Bayes classifier version of $\\hat{y}$ is then**\n",
    "\n",
    "$$\\hat{y} = \\arg\\max_{y_k}\\left[\\ln p(y=y_k) - \\frac{1}{2}\\ln\\mathrm{det}(\\Sigma_k) - \\frac{1}{2}(\\vec{x}-\\vec{\\mu}_k)^T \\Sigma_k^{-1}(\\vec{x}-\\vec{\\mu}_k)\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear & Quadratic Discriminant Analysis <a class=\"anchor\" id=\"six\"></a>\n",
    "\n",
    "In **[Linear Discriminant Analysis (LDA)](https://en.wikipedia.org/wiki/Linear_discriminant_analysis)** we simplify the Gaussian Bayes classifier by assuming that the class distributions have ***identical***\n",
    "covariances for all $k$ classes (all classes are a set of shifted Gaussians). \n",
    "\n",
    "**If we ignore terms involving only the data, $X$, and the class-dependent covariances, then this cancels out terms that have quadratic dependence on $X$.** The log of the class posteriors\n",
    "\n",
    "$$ g_k(\\vec{x}) = \\vec{x}^T \\Sigma^{-1} \\vec{\\mu_k} - \\frac{1}{2}\\vec{\\mu_k}^T \\Sigma^{-1} \\vec{\\mu_k} + \\ln p(y=y_k),$$\n",
    "\n",
    "with $\\vec{\\mu}_k$ the mean of class $k$ and $\\Sigma_k$ the covariance of the Gaussians.\n",
    "\n",
    "<!--- ** note different from book --->\n",
    "\n",
    "***The Bayes classifier is, therefore, linear with respect to $X$***, and  discriminant boundary between classes is the line that minimizes the overlap between Gaussians.\n",
    "\n",
    "<!--- > $  g_k(\\vec{x}) - g_\\ell(\\vec{x}) = \\vec{x}^T \\Sigma^{-1} (\\mu_k-\\mu_\\ell)  - \\frac{1}{2}(\\mu_k - \\mu_\\ell)^T \\Sigma^{-1}(\\mu_k -\\mu_\\ell)  + \\log (\\frac{\\pi_k}{\\pi_\\ell}) = 0. $ --->\n",
    "\n",
    "**Relaxing the requirement that the covariances of the Gaussians are equal, the discriminant function becomes quadratic in $X$.**\n",
    "\n",
    "$$ g(\\vec{x}) = -\\frac{1}{2} \\ln\\mathrm{det}(\\Sigma_k) - \\frac{1}{2}(\\vec{x}-\\mu_k)^T \\Sigma_k^{-1}(\\vec{x}-\\mu_k) + \\ln p(y=y_k). $$\n",
    "\n",
    "This is sometimes known as **[Quadratic Discriminant Analysis (QDA)](https://en.wikipedia.org/wiki/Quadratic_classifier#Quadratic_discriminant_analysis)**.\n",
    "\n",
    "[`LDA`](http://scikit-learn.org/0.16/modules/generated/sklearn.lda.LDA.html) and [`QDA`](http://scikit-learn.org/0.16/modules/generated/sklearn.qda.QDA.html#sklearn.qda.QDA) are implemented in `Scikit-Learn` as follows and an example using the same data as above is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\n",
    "from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis as QDA\n",
    "\n",
    "X = np.random.random((100,2))\n",
    "y = (X[:,0] + X[:,1] > 1).astype(int)\n",
    "\n",
    "lda = LDA()\n",
    "lda.fit(X,y)\n",
    "y_pred = lda.predict(X)\n",
    "\n",
    "qda = QDA()\n",
    "qda.fit(X,y)\n",
    "y_pred = qda.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic, Figures 9.4 and 9.5, spliced together by GTR\n",
    "\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import split_samples\n",
    "from astroML.utils import completeness_contamination\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined()\n",
    "X = X[:, [1, 0, 2, 3]]  # rearrange columns for better 1-color results\n",
    "(X_train, X_test), (y_train, y_test) = split_samples(X, y, [0.75, 0.25], \n",
    "                                                     random_state=0)\n",
    "\n",
    "N_tot = len(y)\n",
    "N_stars = np.sum(y == 0)\n",
    "N_rrlyrae = N_tot - N_stars\n",
    "N_train = len(y_train)\n",
    "N_test = len(y_test)\n",
    "N_plot = 5000 + N_rrlyrae\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# perform LDA\n",
    "lda = LDA()\n",
    "lda.fit(X_train[:, :2], y_train)\n",
    "y_predLDA = lda.predict(X_test[:, :2])\n",
    "\n",
    "# perform QDA\n",
    "qda = QDA()\n",
    "qda.fit(X_train[:, :2], y_train)\n",
    "y_predQDA = qda.predict(X_test[:, :2])\n",
    "    \n",
    "completenessLDA, contaminationLDA = completeness_contamination(y_predLDA, y_test)\n",
    "completenessQDA, contaminationQDA = completeness_contamination(y_predQDA, y_test)\n",
    "\n",
    "print(\"completeness\", completenessLDA, completenessQDA)\n",
    "print(\"contamination\", contaminationLDA, contaminationQDA)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the decision boundary\n",
    "xlim = (0.7, 1.35)\n",
    "ylim = (-0.15, 0.4)\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xlim[0], xlim[1], 71),\n",
    "                     np.linspace(ylim[0], ylim[1], 81))\n",
    "\n",
    "Z_LDA = lda.predict_proba(np.c_[yy.ravel(), xx.ravel()])\n",
    "Z_LDA = Z_LDA[:, 1].reshape(xx.shape)\n",
    "Z_QDA = qda.predict_proba(np.c_[yy.ravel(), xx.ravel()])\n",
    "Z_QDA = Z_QDA[:, 1].reshape(xx.shape)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(10, 5))\n",
    "fig.subplots_adjust(bottom=0.15, top=0.95, hspace=0.0,\n",
    "                    left=0.1, right=0.95, wspace=0.2)\n",
    "\n",
    "# left plot: data and decision boundary\n",
    "ax = fig.add_subplot(121)\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.Oranges, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "im = ax.imshow(Z_LDA, origin='lower', aspect='auto',\n",
    "               cmap=plt.cm.binary, zorder=1,\n",
    "               extent=xlim + ylim)\n",
    "im.set_clim(0, 1.5)\n",
    "\n",
    "ax.contour(xx, yy, Z_LDA, [0.5], linewidths=2., colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$u-g$')\n",
    "ax.set_ylabel('$g-r$')\n",
    "\n",
    "# right plot: qda\n",
    "ax = fig.add_subplot(122)\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.Oranges, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "im = ax.imshow(Z_QDA, origin='lower', aspect='auto',\n",
    "               cmap=plt.cm.binary, zorder=1,\n",
    "               extent=xlim + ylim)\n",
    "im.set_clim(0, 1.5)\n",
    "\n",
    "ax.contour(xx, yy, Z_QDA, [0.5], linewidths=2., colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$u-g$')\n",
    "ax.set_ylabel('$g-r$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "If it is obvious from looking at your data that a linear or quadratic boundary will work well, then great.  But what if that is not the case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## GMM and Bayes Classification <a class=\"anchor\" id=\"seven\"></a>\n",
    "\n",
    "Our classifications so far have made some restrictive assumptions. Either\n",
    "- conditional independence, or \n",
    "- Gaussianity of the distributions.  \n",
    "\n",
    "However, a more flexible model might improve the completeness and efficiency of the classification. For that we can look to the **density estimation techniques** from Chapter 6.\n",
    "\n",
    "The natural extension of the Gaussian assumptions is to use GMM's (Gaussian Mixture Models) to determine the density distribution, i.e., a **GMM Bayes Classifier**. Note that the number of Gaussian components, $K$, must be chosen for each class, $k$, independently.\n",
    "\n",
    "`astroML` implements GMM Bayes classification as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from astroML.classification import GMMBayes\n",
    "#from astroML.classification import GaussianMixture as GMMBayes\n",
    "\n",
    "X = np.random.random((100,2))\n",
    "y = (X[:,0] + X[:,1] > 1).astype(int)\n",
    "\n",
    "gmmb = GMMBayes(3) # 3 clusters per class\n",
    "gmmb.fit(X,y)\n",
    "\n",
    "y_pred = gmmb.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We now apply the GMM Bayes classifier to the real stellar data from above. ***With just one component, we get results that are similar to those from Naive Bayes. But with 5 components (and all 4 attributes), we do pretty well.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 9.6, edits by GTR and SRT\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "from astroML.classification import GMMBayes\n",
    "from astroML.utils.decorators import pickle_results\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import split_samples\n",
    "from astroML.utils import completeness_contamination\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=12, usetex=False)\n",
    "\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined()\n",
    "X = X[:, [1, 0, 2, 3]]  # rearrange columns for better 1-color results\n",
    "\n",
    "# GMM-bayes takes several minutes to run, and is order[N^2]\n",
    "#  truncating the dataset can be useful for experimentation.\n",
    "#X = X[::10]\n",
    "#y = y[::10]\n",
    "\n",
    "(X_train, X_test), (y_train, y_test) = split_samples(X, y, [0.75, 0.25],\n",
    "                                                     random_state=0)\n",
    "N_tot = len(y)\n",
    "N_st = np.sum(y == 0)\n",
    "N_rr = N_tot - N_st\n",
    "N_train = len(y_train)\n",
    "N_test = len(y_test)\n",
    "N_plot = 5000 + N_rr\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# perform GMM Bayes\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "Ncomp = [1, 5]\n",
    "\n",
    "\n",
    "@pickle_results('GMMbayes_rrlyrae.pkl')\n",
    "def compute_GMMbayes(Ncolors, Ncomp):\n",
    "    classifiers = []\n",
    "    predictions = []\n",
    "\n",
    "    for ncm in Ncomp:\n",
    "        classifiers.append([])\n",
    "        predictions.append([])\n",
    "        for nc in Ncolors:\n",
    "            clf = GMMBayes(ncm, tol=1E-5, covariance_type='full')\n",
    "            clf.fit(X_train[:, :nc], y_train)\n",
    "            y_pred = clf.predict(X_test[:, :nc])\n",
    "\n",
    "            classifiers[-1].append(clf)\n",
    "            predictions[-1].append(y_pred)\n",
    "\n",
    "    return classifiers, predictions\n",
    "\n",
    "classifiers, predictions = compute_GMMbayes(Ncolors, Ncomp)\n",
    "\n",
    "completeness, contamination = completeness_contamination(predictions, y_test)\n",
    "\n",
    "print(\"completeness\", completeness)\n",
    "print(\"contamination\", contamination)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the decision boundary\n",
    "clf = classifiers[1][1]\n",
    "xlim = (0.7, 1.35)\n",
    "ylim = (-0.15, 0.4)\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xlim[0], xlim[1], 71),\n",
    "                     np.linspace(ylim[0], ylim[1], 81))\n",
    "\n",
    "Z = clf.predict_proba(np.c_[yy.ravel(), xx.ravel()])\n",
    "Z = Z[:, 1].reshape(xx.shape)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(10, 5))\n",
    "fig.subplots_adjust(bottom=0.15, top=0.95, hspace=0.0,\n",
    "                    left=0.1, right=0.95, wspace=0.2)\n",
    "\n",
    "# left plot: data and decision boundary\n",
    "ax = fig.add_subplot(121)\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.Oranges, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "im = ax.imshow(Z, origin='lower', aspect='auto',\n",
    "               cmap=plt.cm.binary, zorder=1,\n",
    "               extent=xlim + ylim)\n",
    "im.set_clim(0, 1.5)\n",
    "\n",
    "ax.contour(xx, yy, Z, [0.1], colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$u-g$')\n",
    "ax.set_ylabel('$g-r$')\n",
    "\n",
    "# plot completeness vs Ncolors\n",
    "ax = fig.add_subplot(222)\n",
    "ax.plot(Ncolors, completeness[0], '^--k', ms=6, label='N=%i' % Ncomp[0])\n",
    "ax.plot(Ncolors, completeness[1], 'o-k', ms=6, label='N=%i' % Ncomp[1])\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "\n",
    "ax.set_ylabel('completeness')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "# plot contamination vs Ncolors\n",
    "ax = fig.add_subplot(224)\n",
    "ax.plot(Ncolors, contamination[0], '^--k', ms=6, label='N=%i' % Ncomp[0])\n",
    "ax.plot(Ncolors, contamination[1], 'o-k', ms=6, label='N=%i' % Ncomp[1])\n",
    "ax.legend(loc='lower right',\n",
    "          bbox_to_anchor=(1.0, 0.78))\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%i'))\n",
    "\n",
    "ax.set_xlabel('N colors')\n",
    "ax.set_ylabel('contamination')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**NOTE:** We can take this to the extreme by having one mixture component at each training point. We also don't have to restrict ourselves to a Gaussian kernel, we can use any kernel that we like. The resulting ***non-parametric*** Bayes classifier is referred to as **Kernel Discriminant Analysis (KDA)**. \n",
    "\n",
    "It seems like this would be a *lot* more computationally intensive, but now we don't have to optimize the locations of the components...we just need to determine the bandwidth of the kernel. In the end, it can result in better classification.\n",
    "\n",
    "One of the tricks to speed things up is that we don't need to know the actually class probability, we just need to know which is larger.  This is explained in more detail in [Riegel, Gray, & Richards 2008](http://epubs.siam.org/doi/abs/10.1137/1.9781611972788.19), and it is implemented in a series of papers starting with [Richards et al. 2004](http://adsabs.harvard.edu/abs/2004ApJS..155..257R).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## K-Nearest Neighbor Classifier <a class=\"anchor\" id=\"eight\"></a>\n",
    "\n",
    "If we did KDA with a variable bandwidth that depended only on the distance of the nearest neighbor, then we'd have what we call a **Nearest-Neighbor** classifier. Here if $x$ is close to $x'$, then $p(y|x) \\approx p(y|x')$, i.e. **we use the class label of the nearest point**. Note that we have not assumed anything about the conditional density distribution, so it is completely non-parametric.\n",
    "\n",
    "The number of neighbors, $K$, regulates the complexity of the classification, where a larger $K$ decreases the variance in the classification but leads to an increase in the bias. In the simplest form, a majority rule classification is adopted, where each of the $K$ neighbors votes on the classification (N.B., the 3rd different use of $K$ or $k$ in this notebook!)\n",
    "\n",
    "The distance measure is usually N-dimensional Euclidean. However, if the attributes have very different properties, then normalization, weighting, etc. may be needed.\n",
    "\n",
    "Scikit-learn implements **[`K-Nearest Neighbors`](http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsClassifier.html) classification** as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "\n",
    "X = np.random.random((100,2))\n",
    "y = (X[:,0] + X[:,1] > 1).astype(int)\n",
    "\n",
    "knc = KNeighborsClassifier(5) # use 5 nearest neighbors\n",
    "knc.fit(X,y)\n",
    "\n",
    "y_pred = knc.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Implementing it for the same example as above shows that it isn't all that great for this particular case.  See below.  We probably need more training data to reduce the variance for it to work better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Ivezic v2, Figure 9.7, edits by GTR\n",
    "# Author: Jake VanderPlas\n",
    "# License: BSD\n",
    "#   The figure produced by this code is published in the textbook\n",
    "#   \"Statistics, Data Mining, and Machine Learning in Astronomy\" (2013)\n",
    "#   For more information, see http://astroML.github.com\n",
    "#   To report a bug or issue, use the following forum:\n",
    "#    https://groups.google.com/forum/#!forum/astroml-general\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import split_samples\n",
    "from astroML.utils import completeness_contamination\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# This function adjusts matplotlib settings for a uniform feel in the textbook.\n",
    "# Note that with usetex=True, fonts are rendered with LaTeX.  This may\n",
    "# result in an error if LaTeX is not installed on your system.  In that case,\n",
    "# you can set usetex to False.\n",
    "if \"setup_text_plots\" not in globals():\n",
    "    from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined()\n",
    "X = X[:, [1, 0, 2, 3]]  # rearrange columns for better 1-color results\n",
    "(X_train, X_test), (y_train, y_test) = split_samples(X, y, [0.75, 0.25],\n",
    "                                                     random_state=0)\n",
    "\n",
    "N_tot = len(y)\n",
    "N_st = np.sum(y == 0)\n",
    "N_rr = N_tot - N_st\n",
    "N_train = len(y_train)\n",
    "N_test = len(y_test)\n",
    "N_plot = 5000 + N_rr\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# perform Classification\n",
    "\n",
    "classifiers = []\n",
    "predictions = []\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "kvals = [1, 8]\n",
    "\n",
    "for k in kvals:\n",
    "    classifiers.append([])\n",
    "    predictions.append([])\n",
    "    for nc in Ncolors:\n",
    "        clf = KNeighborsClassifier(n_neighbors=k)\n",
    "        clf.fit(X_train[:, :nc], y_train)\n",
    "        y_pred = clf.predict(X_test[:, :nc])\n",
    "\n",
    "        classifiers[-1].append(clf)\n",
    "        predictions[-1].append(y_pred)\n",
    "\n",
    "completeness, contamination = completeness_contamination(predictions, y_test)\n",
    "\n",
    "print(\"completeness\", completeness)\n",
    "print(\"contamination\", contamination)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Compute the decision boundary\n",
    "clf = classifiers[1][1]\n",
    "xlim = (0.7, 1.35)\n",
    "ylim = (-0.15, 0.4)\n",
    "\n",
    "xx, yy = np.meshgrid(np.linspace(xlim[0], xlim[1], 71),\n",
    "                     np.linspace(ylim[0], ylim[1], 81))\n",
    "\n",
    "Z = clf.predict(np.c_[yy.ravel(), xx.ravel()])\n",
    "Z = Z.reshape(xx.shape)\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(10, 5))\n",
    "fig.subplots_adjust(bottom=0.15, top=0.95, hspace=0.0,\n",
    "                    left=0.1, right=0.95, wspace=0.2)\n",
    "\n",
    "# left plot: data and decision boundary\n",
    "ax = fig.add_subplot(121)\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.Oranges, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "im = ax.imshow(Z, origin='lower', aspect='auto',\n",
    "               cmap=plt.cm.binary, zorder=1,\n",
    "               extent=xlim + ylim)\n",
    "im.set_clim(0, 2)\n",
    "\n",
    "ax.contour(xx, yy, Z, [0.1], colors='k')\n",
    "\n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "\n",
    "ax.set_xlabel('$u-g$')\n",
    "ax.set_ylabel('$g-r$')\n",
    "\n",
    "ax.text(0.02, 0.02, \"k = %i\" % kvals[1],\n",
    "        transform=ax.transAxes)\n",
    "\n",
    "# plot completeness vs Ncolors\n",
    "ax = fig.add_subplot(222)\n",
    "\n",
    "ax.plot(Ncolors, completeness[0], 'o-k', ms=6, label='k=%i' % kvals[0])\n",
    "ax.plot(Ncolors, completeness[1], '^--k', ms=6, label='k=%i' % kvals[1])\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "\n",
    "ax.set_ylabel('completeness')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "# plot contamination vs Ncolors\n",
    "ax = fig.add_subplot(224)\n",
    "ax.plot(Ncolors, contamination[0], 'o-k', ms=6, label='k=%i' % kvals[0])\n",
    "ax.plot(Ncolors, contamination[1], '^--k', ms=6, label='k=%i' % kvals[1])\n",
    "ax.legend(loc='lower right',\n",
    "          bbox_to_anchor=(1.0, 0.79))\n",
    "\n",
    "ax.xaxis.set_major_locator(plt.MultipleLocator(1))\n",
    "ax.yaxis.set_major_locator(plt.MultipleLocator(0.2))\n",
    "ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%i'))\n",
    "ax.set_xlabel('N colors')\n",
    "ax.set_ylabel('contamination')\n",
    "ax.set_xlim(0.5, 4.5)\n",
    "ax.set_ylim(-0.1, 1.1)\n",
    "ax.grid(True)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$K=1$ is clearly terrible, but even $K=8$ is not that great. <font color='red'>Where do you think the decision to try $K=8$ came from? Before moving on, discuss with your breakout room the strategy you would use.</font>\n",
    "\n",
    "Regardless of whether this is the best algorithm or not, we can choose $K$ to minimize the classification error rate by using cross-validation. See below for how this is computed (may take a few minutes to run)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import completeness_contamination\n",
    "\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# get data and split into training & testing sets\n",
    "X, y = fetch_rrlyrae_combined()\n",
    "X = X[:, [1, 0, 2, 3]]  # rearrange columns for better 1-color results\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# perform Classification\n",
    "scores = []\n",
    "kvals = np.arange(1,20)\n",
    "for k in kvals:\n",
    "    clf = KNeighborsClassifier(n_neighbors=k)\n",
    "    CVpredk = cross_val_predict(clf, X, y)\n",
    "    scores.append(accuracy_score(y, CVpredk)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(\"max score is for k={:d}\".format(kvals[np.argmax(scores)]))\n",
    "\n",
    "# Plot number of neighbors vs score\n",
    "fig = plt.figure(figsize=(7, 3))\n",
    "u = np.arange(len(scores))+1\n",
    "plt.plot(u,scores)\n",
    "plt.xlabel('$K$ nearest neighbors')\n",
    "plt.ylabel('Classification accuracy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Ah, that's where $K=8$ came from. \n",
    "\n",
    "We can also use the [`metrics` module](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics) in sklearn to compute some statistics for us. Remember the definitions of **precision** and **recall** from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import metrics \n",
    "\n",
    "#print(k, nc)\n",
    "print(\"accuracy:\", metrics.accuracy_score(y_test, y_pred))\n",
    "print(\"precision:\", metrics.precision_score(y_test, y_pred))\n",
    "print(\"recall:\", metrics.recall_score(y_test, y_pred))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:astr8070] *",
   "language": "python",
   "name": "conda-env-astr8070-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}