lectures/Day3-ZeljkoIvezic/notebooks/classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "[The 1st IAA-CSIC Severo Ochoa School on Statistics, Data Mining and Machine Learning](https://www.granadacongresos.com/sostat)\n",
    "\n",
    "[Zeljko Ivezic, University of Washington](http://faculty.washington.edu/ivezic/) \n",
    "\n",
    "[This notebook](https://github.com/carmensg/IAA_School2019/tree/master/lectures/Day3-ZeljkoIvezic/notebooks/classification.ipynb)\n",
    "### November 6, 2019 \n",
    "\n",
    "# Lecture 3: Classification\n",
    "\n",
    "### o Introduction to classification\n",
    "### o Classification loss (risk)\n",
    "### o Receiver Operating Characteristic curves\n",
    "### o Generative Classification and examples\n",
    "### o Discriminative Classification and examples\n",
    "### o Decision Trees\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# first things first...\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats\n",
    "from scipy import integrate\n",
    "from scipy.stats import norm\n",
    "from scipy.stats import cauchy\n",
    "# astroML tools \n",
    "from astroML.plotting import hist\n",
    "from astroML.plotting.mcmc import convert_to_stdev\n",
    "from astroML.plotting import setup_text_plots\n",
    "from astroML.plotting.mcmc import convert_to_stdev\n",
    "setup_text_plots(fontsize=8, usetex=True)\n",
    "import warnings; warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Classification and density estimation\n",
    "\n",
    "In density estimation we estimate joint probability distributions from multivariate data sets to identify the inherent clustering. This is essentially <u> unsupervised classification </u>\n",
    "\n",
    "If we have labels for some of these data points (e.g., an object is tall, short, red, or blue) we can develop a relationship between the label and the properties of a source. This is <u> supervised classification </u>\n",
    "\n",
    "Here “supervised” means that there is prior information about the number and properties of clusters: for a training sample, we know the so-called “class labels” (for each data point in the training sample, we know to which cluster it belongs; characterizing these known clusters is typically easier than finding and characterizing unknown clusters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Classification, regression, and density estimation are all related. For example, the regression function  \n",
    "$ŷ =f(y|x⃗ )$ is the best estimated value of $y$ given a value of $x⃗$. \n",
    "\n",
    "In classification $y$ is categorical and $f(y|x⃗)$ is the called the discriminant function\n",
    "\n",
    "Using density estimation for classification is referred to as **generative classification** (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",
    "Classification that finds the decision boundary that separates classes is called **discriminative classification** \n",
    "\n",
    "Both have their place and role in astrophysical classification. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    " ## Classification loss: how well are we doing\n",
    "\n",
    "The first question we need to address is how we score (defined the success of our classification)\n",
    "\n",
    "We can define a _loss function_. A zero-one loss function assigns a value of one for a misclassification and zero for a correct classification (i.e. we will want to minimize the loss).\n",
    "\n",
    "If $\\hat{y}$ is the best guess value of $y$, the classification loss, $L(y,\\widehat{y})$, is\n",
    "\n",
    "> $L(y,\\widehat{y}) = \\delta(y \\neq \\widehat{y})$\n",
    "\n",
    "which means\n",
    "\n",
    ">$\\begin{eqnarray} L(y,\\hat{y}) & = & \\left\\{ \\begin{array}{cl} 1 & \\mbox{if $y\\neq\\hat{y}$}, \\\\  0 & \\mbox{otherwise.}       \t           \\end{array} \\right. \\end{eqnarray}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "The expectation (mean) value of the loss $\\mathbb{E} \\left[ L(y,\\hat{y}) \\right] = p(y\\neq \\hat{y})$  is called the <u>classification risk</u> \n",
    "\n",
    "This is related to regression loss functions: $L(y, \\hat{y}) = (y - \\hat{y})^2$ and <u>risk<u> $\\mathbb{E}[(y - \\hat{y})^2]$.\n",
    "\n",
    "We can then define:\n",
    "\n",
    "> $ {\\rm completeness} = \\frac{\\rm true\\ positives}\n",
    "  {\\rm true\\ positives + false\\ negatives}\n",
    "$\n",
    "\n",
    "> $  {\\rm contamination} = \\frac{\\rm false\\ positives}\n",
    "  {\\rm true\\ positives + false\\ positives}\n",
    "$\n",
    "\n",
    "or\n",
    "\n",
    "> $ {\\rm true\\ positive\\ rate} = \\frac{\\rm true\\ positives}\n",
    "  {\\rm true\\ positives + false\\ negatives}\n",
    "$\n",
    "\n",
    "> $  {\\rm false\\ positive\\ rate} = \\frac{\\rm false\\ positives}\n",
    "  {\\rm true\\ negatives + false\\ positives}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparing the performance of classifiers\n",
    "\n",
    "Best performance is a bit of a subjective topic (e.g. star-galaxy separation for correlation function studies or Galactic streams studies). We trade contamination as a function of completeness and this is science dependent.\n",
    "\n",
    "**ROC curves: Receiver Operating Characteristic curves**\n",
    "\n",
    "- Plot the true-positive vs the false-positive rate\n",
    "\n",
    "- Initially used to analyze radar results in WWII (a very productive era for statistics...).\n",
    "\n",
    "- One concern about ROC curves is that they are sensitive to the relative sample sizes (if there are many more background events than source events small false positive results can dominate a signal). For these cases we we can plot efficiency (1 - contamination) vs completeness\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"figures/ROC.png\" width=100%>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 9.17: \n",
    "### https://www.astroml.org/book_figures/chapter9/fig_ROC_curve.html \n",
    "from __future__ import print_function\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 sklearn.metrics import roc_curve\n",
    "from astroML.classification import GMMBayes\n",
    "from astroML.utils import split_samples, completeness_contamination\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "setup_text_plots(fontsize=16, usetex=True)\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",
    "# 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",
    "        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=(14, 7))\n",
    "fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15, top=0.9, wspace=0.25)\n",
    "\n",
    "\n",
    "# ax2 will show roc curves\n",
    "ax1 = plt.subplot(121)\n",
    "\n",
    "# ax1 will show completeness/efficiency\n",
    "ax2 = plt.subplot(122)\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",
    "\n",
    "    # add (0, 0) as first point\n",
    "    fpr = np.concatenate([[0], fpr])\n",
    "    tpr = np.concatenate([[0], tpr])\n",
    "\n",
    "    ax1.plot(fpr, tpr, label=labels[name])\n",
    "\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",
    "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",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Generative Classification\n",
    "\n",
    "In generative classifiers we model class-conditional densities $p_k(\\vec{x})$ given $p(\\vec{x}|y=y_k)$ \n",
    "\n",
    "$p(y=y_k)$, or $\\pi_k$ for short, is the probability of any point having class $k$ (equivalent to the prior probability of the class $k$).  \n",
    "\n",
    "Our goal is to estimate the $p_k$'s \n",
    "\n",
    "<u> The discriminative function </u>\n",
    "\n",
    "$\\hat{y} = f(y|\\vec{x})$ represents the best guess of $y$ given a value of $\\vec{x}$.\n",
    "\n",
    "$f(y|\\vec{x})$ is the _discriminant function_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<u>For a simple 2-class example</u>\n",
    "\n",
    "> $\\begin{eqnarray}\n",
    "g(\\vec{x}) & = &  \\int y \\, p(y|\\vec{x}) \\, dy \\\\\n",
    "%    & = & \\int y p(y|\\vec{x}) \\, dy \\\\\n",
    "       & = & 1 \\cdot p(y=1 | \\vec{x}) + 0 \\cdot p(y=0 | \\vec{x}) = p(y=1 | \\vec{x}).\n",
    "%     & = & p(y=1 | \\vec{x})\n",
    "\\end{eqnarray}\n",
    "$\n",
    "\n",
    "From Bayes rule\n",
    "\n",
    "> $\\begin{eqnarray} g(\\vec{x})     & = & \\frac{p(\\vec{x}|y=1) \\, p(y=1)}{p(\\vec{x}|y=1) \\, p(y=1)  + p(\\vec{x}|y=0) \\, p(y=0)} \\\\         & = & \\frac{\\pi_1 p_1(\\vec{x})}{\\pi_1 p_1(\\vec{x}) + \\pi_0 p_0(\\vec{x})} \\end{eqnarray}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<u> Bayes Classifier </u>\n",
    "\n",
    "The discriminant function gives a binary predictor called a Bayes classifier\n",
    "\n",
    ">$\\begin{eqnarray} \\widehat{y} & = & \\left\\{ \\begin{array}{cl}       \t           1 & \\mbox{if $g(\\vec{x}) > 1/2$}, \\\\       \t           0 & \\mbox{otherwise,}       \t           \\end{array}       \t   \\right. \\\\     & = & \\left\\{ \\begin{array}{cl}               1 & \\mbox{if $p(y=1|\\vec{x}) > p(y=0|\\vec{x})$}, \\\\               0 & \\mbox{otherwise,}               \\end{array}       \\right. \\\\     & = & \\left\\{ \\begin{array}{cl}               1 & \\mbox{if $\\pi_1 p_1(\\vec{x}) > \\pi_0 p_0(\\vec{x})$}, \\\\               0 & \\mbox{otherwise.}               \\end{array}       \\right.\\end{eqnarray}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<u> Decision Boundary </u>\n",
    "\n",
    "A set of $x$ values at which each class is equally likely;\n",
    "\n",
    ">$\n",
    "\\pi_1 p_1(\\vec{x}) = \\pi_2 p_2(\\vec{x});\n",
    "$\n",
    "\n",
    "> $g_1(\\vec{x}) = g_2(\\vec{x})$; $g_1(\\vec{x}) - g_2(\\vec{x}) = 0$;  $g(\\vec{x}) = 1/2$; in a two-class problem\n",
    " \n",
    "<img src=\"figures/boundary.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simplest classifier: Naive Bayes\n",
    "\n",
    "We want $p(x_1,x_2,x_3...x_n|y)$ but if we assume that all attributes are conditionally independent this simplifies 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^0,x^1,x^2,\\ldots,x^N}|y_k) = \\prod_i p(x^i|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_i p(x^i|y_k) p(y_k)}\n",
    "       {\\sum_j \\prod_i p(x^i|y_j) p(y_j)}.\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We calculate the most likely value of $y$ by maximizing over $y_k$,\n",
    "\n",
    "## $ \\widehat{y} = \\arg \\max_{y_k} \\frac{\\prod_i p(x^i|y_k) p(y_k)}\n",
    "        {\\sum_j \\prod_i p(x^i|y_j) p(y_j)},\n",
    "$\n",
    "\n",
    "or\n",
    "\n",
    "## $\\widehat{y} = \\arg \\max_{y_k} \\frac{\\prod_i p_k(x^i) \\pi_k}\n",
    "        {\\sum_j \\prod_i p_j(x^i) \\pi_j}.\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The rest of the class is now just estimating densities\n",
    "\n",
    "> <b> $p(\\vec{x}|y=y_k)$ and $\\pi_k$ are learned from a set of training data.\n",
    "- $\\pi_k$ is just the frequency of the class $k$ in the training set\n",
    "- $p(\\vec{x}|y=y_k)$ is just the density (probability) of a object with class $k$ having the attributes $x$\n",
    "</b>\n",
    "\n",
    "\n",
    "If the training set does not cover the full parameter space $p_k(x^i)$ can be $0$ for some value of $y_k$ and $x^i$.  The posterior probability is then $p(y_k|\\{x^i\\}) = 0/0$ which is unfortunate! The trick is _Laplace smoothing_ : an offset $\\alpha$ is added to the probability of each bin $p(\\vec{x}|y=y_k)$ for all $i, k$ (equivalent to the addition of a Bayesian prior to the naive Bayes classifier)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian Naive Bayes\n",
    "\n",
    "In Gaussian naive Bayes $p_k(x^i)$ are modeled as one-dimensional normal distributions, with means $\\mu^i_k$ and widths $\\sigma^i_k$. The naive Bayes estimator is then\n",
    "\n",
    "# $\\hat{y} = \\arg\\max_{y_k}\\left[\\ln \\pi_k - \\frac{1}{2}\\sum_{i=1}^N\\left(2\\pi(\\sigma^i_k)^2 + \\frac{(x^i - \\mu^i_k)^2}{(\\sigma^i_k)^2} \\right) \\right]$\n",
    "\n",
    "<b> Note: this is the log of the Bayes criterion with no normalization constant </b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import colors\n",
    "from astroML.datasets import fetch_imaging_sample\n",
    "from astroML.plotting.tools import draw_ellipse\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "def get_stars_and_galaxies(Nstars=10000, Ngals=10000):\n",
    "    \"\"\"Get the subset of star/galaxy data to plot\"\"\"\n",
    "    data = fetch_imaging_sample()\n",
    "\n",
    "    objtype = data['type']\n",
    "\n",
    "    stars = data[objtype == 6][:Nstars]\n",
    "    \n",
    "    galaxies = data[objtype == 3][:Ngals]\n",
    "\n",
    "    return np.concatenate([stars,galaxies]), np.concatenate([np.zeros(len(stars)),np.ones(len(galaxies))])\n",
    "\n",
    "data, y = get_stars_and_galaxies(Nstars=10000, Ngals=10000)\n",
    "\n",
    "# select r model mag and psf - model mag as columns\n",
    "X = np.column_stack((data['rRaw'], data['rRawPSF'] - data['rRaw']))\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 = (14, 23)\n",
    "ylim = (-1, 5)\n",
    "\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",
    "Z = Z[:, 1].reshape(xx.shape)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Oranges, zorder=2)\n",
    "\n",
    "ax.contour(xx, yy, Z, [0.5], linewidths=2., colors='blue')\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": [
    "## GMM and Bayes classification\n",
    "\n",
    "The natural extension to the Gaussian assumptions is to use GMM's to learn the density distribution. \n",
    "\n",
    "The number of Gaussian components $K$ must be chosen for each class independently\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 4.2: \n",
    "### https://www.astroml.org/book_figures/chapter4/fig_GMM_1D.html\n",
    "\n",
    "from sklearn.mixture import GaussianMixture\n",
    "setup_text_plots(fontsize=16, usetex=True)\n",
    "#------------------------------------------------------------\n",
    "# Set up the dataset.\n",
    "#  We'll create our dataset by drawing samples from Gaussians.\n",
    "random_state = np.random.RandomState(seed=1)\n",
    "X = np.concatenate([random_state.normal(-1, 1.5, 350),\n",
    "                    random_state.normal(0, 1, 500),\n",
    "                    random_state.normal(3, 0.5, 150)]).reshape(-1, 1)\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Learn the best-fit GaussianMixture models\n",
    "#  Here we'll use scikit-learn's GaussianMixture model. The fit() method\n",
    "#  uses an Expectation-Maximization approach to find the best\n",
    "#  mixture of Gaussians for the data\n",
    "\n",
    "# fit models with 1-10 components\n",
    "N = np.arange(1, 11)\n",
    "models = [None for i in range(len(N))]\n",
    "\n",
    "for i in range(len(N)):\n",
    "    models[i] = GaussianMixture(N[i]).fit(X)\n",
    "\n",
    "# compute the AIC and the BIC\n",
    "AIC = [m.aic(X) for m in models]\n",
    "BIC = [m.bic(X) for m in models]\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "#  We'll use three panels:\n",
    "#   1) data + best-fit mixture\n",
    "#   2) AIC and BIC vs number of components\n",
    "#   3) probability that a point came from each component\n",
    "\n",
    "fig = plt.figure(figsize=(15, 5.2))\n",
    "fig.subplots_adjust(left=0.12, right=0.97,\n",
    "                    bottom=0.21, top=0.9, wspace=0.5)\n",
    "\n",
    "\n",
    "# plot 1: data + best-fit mixture\n",
    "ax = fig.add_subplot(131)\n",
    "M_best = models[np.argmin(AIC)]\n",
    "\n",
    "x = np.linspace(-6, 6, 1000)\n",
    "logprob = M_best.score_samples(x.reshape(-1, 1))\n",
    "responsibilities = M_best.predict_proba(x.reshape(-1, 1))\n",
    "pdf = np.exp(logprob)\n",
    "pdf_individual = responsibilities * pdf[:, np.newaxis]\n",
    " \n",
    "ax.hist(X, 30, histtype='stepfilled', alpha=0.4, normed=True)\n",
    "ax.plot(x, pdf, '-k')\n",
    "ax.plot(x, pdf_individual, '--k')\n",
    "ax.text(0.04, 0.96, \"Best-fit Mixture\",\n",
    "        ha='left', va='top', transform=ax.transAxes)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$p(x)$')\n",
    "\n",
    "\n",
    "# plot 2: AIC and BIC\n",
    "ax = fig.add_subplot(132)\n",
    "ax.plot(N, AIC, '-k', label='AIC')\n",
    "ax.plot(N, BIC, '--k', label='BIC')\n",
    "ax.set_xlabel('n. components')\n",
    "ax.set_ylabel('information criterion')\n",
    "ax.legend(loc=2)\n",
    "\n",
    "\n",
    "# plot 3: posterior probabilities for each component\n",
    "ax = fig.add_subplot(133)\n",
    "\n",
    "p = responsibilities\n",
    "p = p[:, (1, 0, 2)]  # rearrange order so the plot looks better\n",
    "p = p.cumsum(1).T\n",
    "\n",
    "ax.fill_between(x, 0, p[0], color='gray', alpha=0.3)\n",
    "ax.fill_between(x, p[0], p[1], color='gray', alpha=0.5)\n",
    "ax.fill_between(x, p[1], 1, color='gray', alpha=0.7)\n",
    "ax.set_xlim(-6, 6)\n",
    "ax.set_ylim(0, 1)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel(r'$p({\\rm class}|x)$')\n",
    "\n",
    "ax.text(-5, 0.3, 'class 1', rotation='vertical')\n",
    "ax.text(0, 0.5, 'class 2', rotation='vertical')\n",
    "ax.text(3, 0.3, 'class 3', rotation='vertical')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 9.3: \n",
    "### https://www.astroml.org/book_figures/chapter9/fig_rrlyrae_naivebayes.html \n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\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",
    "setup_text_plots(fontsize=16, usetex=True)\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 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",
    "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",
    "    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=(15, 7.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.binary, 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()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## K-nearest neighbours\n",
    "\n",
    "As with density estimation (and kernel density estimation)  the intuitive justification is that $p(y|x) \\approx p(y|x')$ if $x'$ is very close to $x$.  \n",
    "\n",
    "The number of neighbors, $K$, regulates the complexity of the classification. In simplest form, a majority rule classification is adopted, where each of the $K$ points votes on the classification. Increasing $K$ decreases the variance in the classification but at the expense of an increase in the bias.  \n",
    "\n",
    "Weights can be assigned to  individual votes by weighting the vote by the distance to the nearest point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 9.7: \n",
    "### https://www.astroml.org/book_figures/chapter9/fig_rrlyrae_knn.html\n",
    "from __future__ import print_function\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",
    "setup_text_plots(fontsize=16, usetex=True)\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, 10]\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=(15, 7.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.binary, 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.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",
    "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": [
    "# Discriminative classification\n",
    "\n",
    "Rather than a probabilistc description we directly model the decision boundary between two or more classes\n",
    "\n",
    "For a two class example $\\{0,1\\}$, the discriminant function is given by $g(x) =  p(y=1 | x)$. Once we have it we can use the rule\n",
    "\n",
    "## $ \\widehat{y} = \\left\\{ \\begin{array}{cl}       \t           1 & \\mbox{if $g(x) > 1/2$}, \\\\       \t           0 & \\mbox{otherwise,}       \t           \\end{array}       \t   \\right.$\n",
    "\n",
    "to perform classification."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Logistic regression\n",
    "\n",
    "Assuming a binomial logistic regression and consider the linear model\n",
    "\n",
    "## $\n",
    "\\begin{eqnarray}\n",
    "  p(y=1|x) &=& \\frac{\\exp\\left[\\sum_j \\theta_j x^j\\right]}\n",
    "  {1 + \\exp\\left[\\sum_j \\theta_j x^j\\right]}\\nonumber\\\\\n",
    "  &=& p({\\theta}),\n",
    "\\end{eqnarray}\n",
    "$\n",
    "\n",
    "where we define\n",
    "\n",
    "## $\n",
    "\\mbox{logit}(p_i) = \\log\\left(\\frac{p_i}{1-p_i}\\right)\n",
    "= \\sum_j \\theta_j x_i^j.\n",
    "$\n",
    "\n",
    "The name logistic regression comes from the fact that the function $e^x/(1+e^x)$ is called the logistic function.  \n",
    "- useful for categorical regression as values cannot go above or below 1 and 0.\n",
    "\n",
    "Because $y$ is binary, it can be modeled as a Bernoulli distribution with (conditional) likelihood function\n",
    "\n",
    "## $\n",
    "L(\\beta) = \\prod_{i=1}^N p_i(\\beta)^{y_i} (1-p_i(\\beta))^{1-y_i}.\n",
    "$\n",
    "\n",
    "logistic regression is related to linear discriminant analysis (LDA). In LDA\n",
    "\n",
    "## \\begin{eqnarray}\n",
    "\\log\\left( \\frac{p(y=1|x)}{p(y=0|x)} \\right) & = &\n",
    "  - \\frac{1}{2}(\\mu_0+\\mu_1)^T \\Sigma^{-1} (\\mu_1-\\mu_0) \\\\ \\nonumber\n",
    "  & + & \\log\\left( \\frac{\\pi_0}{\\pi_1} \\right) + x^T \\Sigma^{-1} (\\mu_1-\\mu_0) \\\\ \\nonumber\n",
    "  & = & \\alpha_0 + \\alpha^T x.\n",
    "\\end{eqnarray}\n",
    "\n",
    "In logistic regression the model is by assumption\n",
    "\n",
    "## $\n",
    "\\log\\left( \\frac{p(y=1|x)}{p(y=0|x)} \\right) = \\beta_0 + \\beta^T x.\n",
    "$\n",
    "\n",
    "Logistic regression minimizes classification error rather than density estimation error.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Support Vector Machines\n",
    "<img style=\"float: left; width: 50%\" src=\"figures/svm_lines.png\"> \n",
    "<img style=\"float: left; width: 50%\" src=\"figures/svm.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Find the hyperplane that maximizes the distance of the closest point from either class. This distance is the margin (width of the line before it hits a point). We want the line that maximizes the margin (m).\n",
    "\n",
    "The points on the margin are called _support vectors_\n",
    "\n",
    "If we assume $y \\in \\{-1,1\\}$, (+1 is maximum margin, -1 is minimum, 0 is the decision boundary)\n",
    "\n",
    "The maximum is then just when $\\beta_0 + \\beta^T x_i = 1$ etc\n",
    "\n",
    "The hyperplane which maximizes the margin is given by finding\n",
    "\n",
    "> \\begin{equation}\n",
    "\\max_{\\beta_0,\\beta}(m) \\;\\;\\;\n",
    "  \\mbox{subject to} \\;\\;\\; \\frac{1}{||\\beta||} y_i ( \\beta_0 + \\beta^T x_i )\n",
    "  \\geq m \\,\\,\\, \\forall \\, i.\n",
    "\\end{equation}\n",
    "\n",
    "The constraints can be written as $y_i ( \\beta_0 + \\beta^T x_i ) \\geq m ||\\beta|| $. \n",
    "\n",
    "Thus the optimization problem is equivalent to minimizing\n",
    "> \\begin{equation}\n",
    "\\frac{1}{2} ||\\beta|| \\;\\;\\; \\mbox{subject to} \\;\\;\\; y_i\n",
    "  ( \\beta_0 + \\beta^T x_i ) \\geq 1 \\,\\,\\, \\forall \\, i.\n",
    "\\end{equation}\n",
    "\n",
    "This optimization  is a _quadratic programming_ problem (quadratic objective function with linear constraints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that because SVM uses a metric which maximizes the margin rather than a measure over all points in the data sets, it is similar in spirit to the rank-based estimators \n",
    "\n",
    "- The median of a distribution is unaffected by even large perturbations of outlying points, as long as those perturbations do not cross the boundary.\n",
    "- In the same way, once the support vectors are determined, changes to the positions or numbers of points beyond the margin will not change the decision boundary.  For this reason, SVM can be a very powerful tool for discriminative classification.\n",
    "\n",
    "- This is why there is a high completeness compared to the other methods: it does not matter that the background sources outnumber the RR Lyrae stars by a factor of $\\sim$200 to 1. It simply determines the best boundary between the small RR Lyrae clump and the large background clump.\n",
    "- This completeness, however, comes at the cost of a relatively large contamination level.\n",
    "\n",
    "- SVM is not scale invariant so it often worth rescaling the data to [0,1] or to whiten it to have a mean of 0 and variance 1 (remember to do this to the test data as well!)\n",
    "- The data dont need to be separable (we can put a constraint in minimizing the number of \"failures\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 9.10: \n",
    "### https://www.astroml.org/book_figures/chapter9/fig_rrlyrae_svm.html\n",
    "from __future__ import print_function\n",
    "from sklearn.svm import SVC\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",
    "setup_text_plots(fontsize=16, usetex=True)\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",
    "# SVM takes several minutes to run, and is order[N^2]\n",
    "#  truncating the dataset can be useful for experimentation.\n",
    "#X = X[::5]\n",
    "#y = y[::5]\n",
    "\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",
    "# Fit SVM\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "\n",
    "\n",
    "@pickle_results('SVM_rrlyrae.pkl')\n",
    "def compute_SVM(Ncolors):\n",
    "    classifiers = []\n",
    "    predictions = []\n",
    "\n",
    "    for nc in Ncolors:\n",
    "        # perform support vector classification\n",
    "        clf = SVC(kernel='linear', class_weight='balanced')\n",
    "        clf.fit(X_train[:, :nc], y_train)\n",
    "        y_pred = clf.predict(X_test[:, :nc])\n",
    "\n",
    "        classifiers.append(clf)\n",
    "        predictions.append(y_pred)\n",
    "\n",
    "    return classifiers, predictions\n",
    "\n",
    "classifiers, predictions = compute_SVM(Ncolors)\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",
    "w = clf.coef_[0]\n",
    "a = -w[0] / w[1]\n",
    "yy = np.linspace(-0.1, 0.4)\n",
    "xx = a * yy - clf.intercept_[0] / w[1]\n",
    "\n",
    "#----------------------------------------------------------------------\n",
    "# plot the results\n",
    "fig = plt.figure(figsize=(15, 7.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",
    "ax.plot(xx, yy, '-k')\n",
    "im = ax.scatter(X[-N_plot:, 1], X[-N_plot:, 0], c=y[-N_plot:],\n",
    "                s=4, lw=0, cmap=plt.cm.binary, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\n",
    "ax.set_xlim(0.7, 1.35)\n",
    "ax.set_ylim(-0.15, 0.4)\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, '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 = fig.add_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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 10.23: \n",
    "### https://www.astroml.org/book_figures/chapter10/fig_LINEAR_SVM.html\n",
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.model_selection import train_test_split\n",
    "from astroML.datasets import fetch_LINEAR_geneva\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=True)\n",
    "\n",
    "data = fetch_LINEAR_geneva()\n",
    "\n",
    "attributes = [('gi', 'logP'),\n",
    "              ('gi', 'logP', 'ug', 'iK', 'JK', 'amp', 'skew')]\n",
    "labels = ['$u-g$', '$g-i$', '$i-K$', '$J-K$',\n",
    "          r'$\\log(P)$', 'amplitude', 'skew']\n",
    "cls = 'LCtype'\n",
    "Ntrain = 3000\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Create attribute arrays\n",
    "X = []\n",
    "y = []\n",
    "\n",
    "for attr in attributes:\n",
    "    X.append(np.vstack([data[a] for a in attr]).T)\n",
    "    LCtype = data[cls].copy()\n",
    "\n",
    "    # there is no #3.  For a better color scheme in plots,\n",
    "    # we'll set 6->3\n",
    "    LCtype[LCtype == 6] = 3\n",
    "    y.append(LCtype)\n",
    "\n",
    "\n",
    "#@pickle_results('LINEAR_SVM.pkl')\n",
    "def compute_SVM_results(i_train, i_test):\n",
    "    classifiers = []\n",
    "    predictions = []\n",
    "    Xtests = []\n",
    "    ytests = []\n",
    "    Xtrains = []\n",
    "    ytrains = []\n",
    "\n",
    "    for i in range(len(attributes)):\n",
    "        Xtrain = X[i][i_train]\n",
    "        Xtest = X[i][i_test]\n",
    "        ytrain = y[i][i_train]\n",
    "        ytest = y[i][i_test]\n",
    "\n",
    "        clf = SVC(kernel='linear', class_weight=None)\n",
    "        clf.fit(Xtrain, ytrain)\n",
    "        y_pred = clf.predict(Xtest)\n",
    "\n",
    "        classifiers.append(clf)\n",
    "        predictions.append(y_pred)\n",
    "\n",
    "    return classifiers, predictions\n",
    "\n",
    "\n",
    "i = np.arange(len(data))\n",
    "i_train, i_test = train_test_split(i, random_state=0, train_size=2000)\n",
    "clfs, ypred = compute_SVM_results(i_train, i_test)\n",
    "\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "fig.subplots_adjust(hspace=0.1, wspace=0.1)\n",
    "\n",
    "class_labels = []\n",
    "\n",
    "for i in range(2):\n",
    "    Xtest = X[i][i_test]\n",
    "    ytest = y[i][i_test]\n",
    "    amp = data['amp'][i_test]\n",
    "\n",
    "    # Plot the resulting classifications\n",
    "    ax1 = fig.add_subplot(221 + 2 * i)\n",
    "    ax1.scatter(Xtest[:, 0], Xtest[:, 1],\n",
    "                c=ypred[i], edgecolors='none', s=4, linewidths=0)\n",
    "\n",
    "    ax1.set_ylabel(r'$\\log(P)$')\n",
    "\n",
    "    ax2 = plt.subplot(222 + 2 * i)\n",
    "    ax2.scatter(amp, Xtest[:, 1],\n",
    "                c=ypred[i], edgecolors='none', s=4, lw=0)\n",
    "\n",
    "    #------------------------------\n",
    "    # set axis limits\n",
    "    ax1.set_xlim(-0.6, 2.1)\n",
    "    ax2.set_xlim(0.1, 1.5)\n",
    "    ax1.set_ylim(-1.5, 0.5)\n",
    "    ax2.set_ylim(-1.5, 0.5)\n",
    "\n",
    "    ax2.yaxis.set_major_formatter(plt.NullFormatter())\n",
    "    if i == 0:\n",
    "        ax1.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "        ax2.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "    else:\n",
    "        ax1.set_xlabel(r'$g-i$')\n",
    "        ax2.set_xlabel(r'$A$')\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Second figure\n",
    "fig = plt.figure(figsize=(5, 5))\n",
    "fig.subplots_adjust(left=0.11, right=0.95, wspace=0.3)\n",
    "\n",
    "attrs = ['skew', 'ug', 'iK', 'JK']\n",
    "labels = ['skew', '$u-g$', '$i-K$', '$J-K$']\n",
    "ylims = [(-1.8, 2.2), (0.6, 2.9), (0.1, 2.6), (-0.2, 1.2)]\n",
    "\n",
    "for i in range(4):\n",
    "    ax = fig.add_subplot(221 + i)\n",
    "    ax.scatter(data['gi'][i_test], data[attrs[i]][i_test],\n",
    "               c=ypred[1], edgecolors='none', s=4, lw=0)\n",
    "    ax.set_xlabel('$g-i$')\n",
    "    ax.set_ylabel(labels[i])\n",
    "\n",
    "    ax.set_xlim(-0.6, 2.1)\n",
    "    ax.set_ylim(ylims[i])\n",
    "\n",
    "#------------------------------------------------------------\n",
    "# Save the results\n",
    "#\n",
    "# run the script as\n",
    "#\n",
    "#   >$ python fig_LINEAR_clustering.py --save\n",
    "#\n",
    "# to output the data file showing the cluster labels of each point\n",
    "import sys\n",
    "if len(sys.argv) > 1 and sys.argv[1] == '--save':\n",
    "    filename = 'cluster_labels_svm.dat'\n",
    "\n",
    "    print(\"Saving cluster labels to\", filename)\n",
    "\n",
    "    from astroML.datasets.LINEAR_sample import ARCHIVE_DTYPE\n",
    "    new_data = np.zeros(len(data),\n",
    "                        dtype=(ARCHIVE_DTYPE + [('2D_cluster_ID', 'i4'),\n",
    "                                                ('7D_cluster_ID', 'i4')]))\n",
    "\n",
    "    # switch the labels back 3->6\n",
    "    for i in range(2):\n",
    "        ypred[i][ypred[i] == 3] = 6\n",
    "\n",
    "    # need to put labels back in order\n",
    "    class_labels = [-999 * np.ones(len(data)) for i in range(2)]\n",
    "    for i in range(2):\n",
    "        class_labels[i][i_test] = ypred[i]\n",
    "\n",
    "    for name in data.dtype.names:\n",
    "        new_data[name] = data[name]\n",
    "    new_data['2D_cluster_ID'] = class_labels[0]\n",
    "    new_data['7D_cluster_ID'] = class_labels[1]\n",
    "\n",
    "    fmt = ('%.6f   %.6f   %.3f   %.3f   %.3f   %.3f   %.7f   %.3f   %.3f   '\n",
    "           '%.3f    %.2f     %i     %i      %s          %i              %i\\n')\n",
    "\n",
    "\n",
    "    F = open(filename, 'w')\n",
    "    F.write('#    ra           dec       ug      gi      iK      JK     '\n",
    "            'logP       Ampl    skew      kurt    magMed    nObs  LCtype  '\n",
    "            'LINEARobjectID  2D_cluster_ID   7D_cluster_ID\\n')\n",
    "    for line in new_data:\n",
    "        F.write(fmt % tuple(line[col] for col in line.dtype.names))\n",
    "    F.close()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Decision Trees\n",
    "\n",
    "The hierarchical application of decision boundaries lead to _decision trees_\n",
    "\n",
    "Tree structure:\n",
    "- top node contains the entire data set\n",
    "- at each branch the data are subdivided into two child nodes \n",
    "- split is based on a predefined decision boundary (usually axis aligned)\n",
    "- splitting repeats, recursively, until we reach a predefined stopping criteria "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"figures/decision.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The \"leaf nodes\" record the fraction of points that have one classification or the other\n",
    "\n",
    "Application of the tree to classification is simple (a series of binary decisions). The  fraction of points from the training set classified as one class or the other (in the leaf node) defines the class associated with that leaf node.\n",
    "\n",
    "Decision trees are simple to interpret (a set of questions)\n",
    "\n",
    "## Splitting Criteria\n",
    "\n",
    "In order to build a decision tree we must choose the feature and\n",
    "value on which we wish to split the data.\n",
    "\n",
    "\n",
    "## Information content or entropy of the data\n",
    "\n",
    "> $\n",
    "    E(x) = -\\sum_i p_i(x) \\ln (p_i(x)),\n",
    "$\n",
    "\n",
    "where $i$ is the class and $p_i(x)$ is the probability of that class\n",
    "given the training data. \n",
    "\n",
    "Information gain is the reduction in entropy due to the partitioning of the data \n",
    "\n",
    "For a binary split $IG(x)$ is\n",
    "\n",
    "> $ IG(x|x_i) = E(x) - \\sum_{i=0}^{1} \\frac{N_i}{N}  E(x_i), $\n",
    "\n",
    "where $N_i$ is the number of points, $x_i$, in the $i$th class,\n",
    "and $E(x_i)$ is the entropy associated with that class\n",
    "\n",
    "Search for the split is undertaken in a greedy fashion (one attribute at a time) where we sort the data on feature $i$ and maximize the information gain for a given split point, $s$,\n",
    "\n",
    ">$\n",
    "IG(x|s) = E(x) - \\arg\\max_s \\left( \\frac{N(x|x<s)}{N}  E(x|x<s) -\n",
    "  \\frac{N(x|x\\ge s)}{N}  E(x|x\\ge s)  \\right ).\n",
    "$\n",
    "\n",
    "To determine the depth of the tree we use cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Modeled after astroML book figure 9.13: \n",
    "### https://www.astroml.org/book_figures/chapter9/fig_rrlyrae_decisiontree.html\n",
    "from __future__ import print_function\n",
    "\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from astroML.datasets import fetch_rrlyrae_combined\n",
    "from astroML.utils import split_samples\n",
    "from astroML.utils import completeness_contamination\n",
    "setup_text_plots(fontsize=16, usetex=True)\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_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",
    "# Fit Decision tree\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "\n",
    "classifiers = []\n",
    "predictions = []\n",
    "Ncolors = np.arange(1, X.shape[1] + 1)\n",
    "depths = [7, 12]\n",
    "\n",
    "for depth in depths:\n",
    "    classifiers.append([])\n",
    "    predictions.append([])\n",
    "    for nc in Ncolors:\n",
    "        clf = DecisionTreeClassifier(random_state=0, max_depth=depth,\n",
    "                                     criterion='entropy')\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",
    "\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], 101),\n",
    "                     np.linspace(ylim[0], ylim[1], 101))\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=(15, 7.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.binary, zorder=2)\n",
    "im.set_clim(-0.5, 1)\n",
    "\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",
    "ax.text(0.02, 0.02, \"depth = %i\" % depths[1],\n",
    "        transform=ax.transAxes)\n",
    "\n",
    "# plot completeness vs Ncolors\n",
    "ax = fig.add_subplot(222)\n",
    "ax.plot(Ncolors, completeness[0], 'o-k', ms=6, label=\"depth=%i\" % depths[0])\n",
    "ax.plot(Ncolors, completeness[1], '^--k', ms=6, label=\"depth=%i\" % depths[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=\"depth=%i\" % depths[0])\n",
    "ax.plot(Ncolors, contamination[1], '^--k', ms=6, label=\"depth=%i\" % depths[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",
    "\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": [
    "## Alright, alright but what the @#%! should I use?\n",
    "\n",
    "A convenient cop-out: no single model can be known in advance to be the best classifier!\n",
    "\n",
    "In general the level of accuracy increases for parametric models as:\n",
    "- <b>naive Bayes</b>, \n",
    "- linear discriminant analysis (LDA),\n",
    "- logistic regression, \n",
    "- linear support vector machines, \n",
    "- quadratic discriminant analysis (QDA),\n",
    "- linear ensembles of linear models. \n",
    "\n",
    "For non-parametric models accuracy increases as:\n",
    "- decision trees\n",
    "- $K$-nearest-neighbor, \n",
    "- neural networks\n",
    "- kernel discriminant analysis,\n",
    "- <b> kernelized support vector machines</b>\n",
    "- <b> random forests</b>\n",
    "- boosting\n",
    "\n",
    "Naive Bayes and its variants are by far the easiest to compute. Linear support vector machines are more expensive, though several fast algorithms exist. Random forests can be easily parallelized. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected",
   "theme": "sky"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}