classification/1-Classification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from astroML.datasets import fetch_sdss_corrected_spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification: part 1\n",
    "\n",
    "The goal of a classification problem is to train a model that will successfully predict the classes of new objects that we feed in to the model. We do this a lot in astronomy! \n",
    "\n",
    "* Type 1 / Type 2 anything\n",
    "* RR lyrae: RRab / RRc / RRd\n",
    "* Seyfert 1 / 2 galaxies\n",
    "* Mini-neptunes / Super Earths / etc.\n",
    "* Blue cloud / red sequence / green valley galaxies\n",
    "* ...the list goes on\n",
    "\n",
    "Sometimes our classes are observationally determined, in which case we then try to make causal relationships to theory. Other times theory predicts clear separation of objects. In either case, we often want to be able to classify new objects, or select objects of a particular class. Sometimes it makes sense to draw a hard boundary and do hard classification, other times we prefer to have probabilities of belonging to classes, i.e. soft classification. There are reasons to do both in science, and there are many well-studied and used algorithms for classification problems.\n",
    "\n",
    "Deciding on a classifier or classification method is generally very context dependent, and depends on whether you want soft/hard classifications, have many or few classes, have many or few data points, have a training set or not. Even within a given model, we often then have to choose a way to \"score\" the performance of the model. As with other machine learning methods, you also have to decide on what features to compute or use from the data or things you want to classify.\n",
    "\n",
    "In choosing a classifier, it's important to weight several things:\n",
    "* Do I care about population statistics? Does this classifier create a really complex selection function?\n",
    "* How fast do I need the prediction to be?\n",
    "* How much data do I have?\n",
    "* How many features do I have? (i.e. dimensionality)\n",
    "\n",
    "There are\n",
    "\n",
    "To name a few classifiers:\n",
    "* Nearest neighbor (discriminative):\n",
    "    * +Very simple\n",
    "    * +Can compare well to very complicated methods if you have a lot of data\n",
    "    * +Non-parametric representation of complex decision boundaries \"for free\"\n",
    "    * -Doesn't scale well to high dimensions\n",
    "    * -Have to weight features or normalize in some way to compute distance\n",
    "    * -Slow prediction (no real training)\n",
    "* (Naive) Bayes (generative):\n",
    "    * +Approximation to Bayesian hierarchical inference\n",
    "    * +Simple, and therefore fast\n",
    "    * +Interpretable, especially if your density distributions are probability distributions\n",
    "    * +Can handle missing data\n",
    "    * -Naive Bayes can't handle correlated features (generic Bayes is fine though)\n",
    "* Decision trees / random forest (discriminative):\n",
    "    * +Non-parametric\n",
    "    * +Handles high-dimensional feature spaces\n",
    "    * +Will tell you what features are most discriminating\n",
    "    * -Have to rebuild full tree any time you get new data\n",
    "    * -Susceptible to over-fitting\n",
    "* Neural networks (discriminative or generative):\n",
    "    * +Non-parametric\n",
    "    * +Can extend to handle convolutional relationships between features\n",
    "    * +All the rage right now, so good tools exist\n",
    "    * -Designing good network architectures is an art, often done by intuition, not theory\n",
    "    * -Generally require a lot of data and time to train\n",
    "    * -Not interpretable\n",
    "    \n",
    "This week, we'll discuss nearest neighbors and simple decision trees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data munging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = fetch_sdss_corrected_spectra()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need do some pre-processing to condense our labels a bit. In the original data, the value of `lineindex_cln` is:\n",
    "* 4 if HII region\n",
    "* 5 if AGN\n",
    "* 0, 1, 2, 3, 6, 9 else (other things)\n",
    "\n",
    "We'll make a new label array, `y`, that is:\n",
    "* 0 if else\n",
    "* 1 if HII\n",
    "* 2 if AGN\n",
    "\n",
    "These classes were determined using the line ratios, but we're going to use the spectral fluxes themselves as features to try to predict these classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nspectra = len(data['lineindex_cln'])\n",
    "y = np.zeros(nspectra)\n",
    "y[data['lineindex_cln'] == 4] = 1\n",
    "y[data['lineindex_cln'] == 5] = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the data look like in the line ratio space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6,5))\n",
    "cs = ax.scatter(data['log_NII_Ha'], data['log_OIII_Hb'], \n",
    "                c=data['z'], vmin=0, vmax=0.2,\n",
    "                marker='o', s=20, linewidth=0, alpha=0.4)\n",
    "ax.set_xlim(-1.5, 0.5)\n",
    "ax.set_ylim(-1.5, 1.5)\n",
    "cb = fig.colorbar(cs)\n",
    "cb.set_label('redshift, $z$')\n",
    "ax.set_xlabel(r'$\\log \\left( [{\\rm NII}] / {\\rm H}\\alpha \\right)$')\n",
    "ax.set_ylabel(r'$\\log \\left( [{\\rm OIII}] / {\\rm H}\\beta \\right)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to prepare two feature sets that we'll use in various applications: \n",
    "1. A small section of the spectrum (~40 pixels or features)\n",
    "2. The full spectra themselves (they have been put onto the same wavelength grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at a few of the spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rnd = np.random.RandomState(seed=42)\n",
    "\n",
    "fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True, sharey=True)\n",
    "\n",
    "for ax, i in zip(axes, rnd.choice(nspectra, size=len(axes), replace=False)):\n",
    "    ax.plot(data['spectra'][i], marker='', drawstyle='steps-mid')\n",
    "    ax.set_ylim(-5, 25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's pack the spectra into feature matrices - in fact, Jake VdP has done that for us already! The full spectra have 1000 pixels (wavelength values), and our mini spectra will be a 40 pixel subset of the full spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_spec_mini = data['spectra'][:, 200:240]\n",
    "X_spec = data['spectra']\n",
    "X_spec_mini.shape, X_spec.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectra have already been cleaned of any NaN or Inf values (set to 0).\n",
    "\n",
    "Let's now take a look at a few applications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# K Nearest Neighbors\n",
    "\n",
    "We'll see how well we can do at predicting whether an object is an AGN vs. an HII region with KNN, after optimizing over the number of neighbors, `K`. For this example, we'll just use the line ratio data.\n",
    "\n",
    "Convention is to define the classifier object as \"`clf`\". For the simplest version of this algorithm, it has 1 tunable parameter: the number of neighbors K to consider. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = KNeighborsClassifier(n_neighbors=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by fitting the classifier to all but the last 128 objects, then predict the classes of the last data points. In the case of KNN, what exactly is KNN doing when we \"fit\" the model? There are no parameters to fit for! Here it's just constructing a tree based on the features to help make the prediction step more efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X_spec_mini[:-128], y[:-128])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = clf.predict(X_spec_mini[-128:])\n",
    "y_true = y[-128:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Predicted:\", y_pred[:16])\n",
    "print(\"True:     \", y_true[:16])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice, our predictions are pretty good! But how do we evaluate how well our classifier is doing? We need to compute some kind of \"score\" for our classifier. In machine learning, this is often the most subtle challenge in the methodology. Be aware that scikit-learn often lets you specify or define your own score functions, but for now we're going to table this idea and stick to some simple defaults.\n",
    "\n",
    "For discrete classification, a common score function is the \"misclassification error.\" This is just the fractional number of data points that we incorrectly classified with whatever classification method we used. We can compute this fairly simply, and then take the compliment of the error to compute the \"accuracy\" of the classifier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "misclassification_error = np.sum(y_true != y_pred) / len(y_true)\n",
    "accuracy = 1 - misclassification_error\n",
    "accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a pretty good accuracy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn also provides a function to compute this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy_score(y_true, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: \n",
    "\n",
    "How does the accuracy change as you increase the number of neighbors (\"K\") from 1-10?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns = np.arange(1, 15+1)\n",
    "accs = []\n",
    "for N in Ns:\n",
    "    clf = KNeighborsClassifier(n_neighbors=N)\n",
    "    clf.fit(X_spec_mini[:-1000], y[:-1000])\n",
    "    \n",
    "    y_pred = clf.predict(X_spec_mini[-1000:])\n",
    "    y_true = y[-1000:]\n",
    "    acc = accuracy_score(y_true, y_pred)\n",
    "    accs.append(acc)\n",
    "    \n",
    "    print(\"n={0}, accuracy={1:.2f}\".format(N, acc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(Ns, accs, drawstyle='steps-mid', marker='')\n",
    "plt.xlabel('$K$ neighbors')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.xlim(0, 15.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Cross-validation\n",
    "\n",
    "Cross-validation provides a way to tune or set hyperparameters of machine learning algorithms to try to avoid overfitting to training data. The way these methods generally work for supervised learning problems (like classification) is by training the model on some subset of the full training set, predicting the labels of the held-back data, computing an accuracy metric by comparing to the true labels, then optimizing over the accuracy metric. With jargon, the labeled dataset is typically split into (at least two) subsets called the _training data_ and _test data_. The models are trained on the training data, and then evaluated on the test data, which the model has never seen before. Provided you separate the train/test data in a sensible way, this gives you an unbiased way to validate the models. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now use Sciki-learn's cross-validation score function to automatically do the train-test splitting (multiple times). How does the cross-validation score depend on `k`, the number of neighbors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vals = np.arange(2, 100, 5)\n",
    "mean_scores = []\n",
    "for n_neighbors in vals:\n",
    "    clf = KNeighborsClassifier(n_neighbors=n_neighbors)\n",
    "    scores = cross_val_score(clf, X_spec_mini, y)\n",
    "    mean_scores.append(np.mean(scores))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(vals, mean_scores)\n",
    "plt.xlabel('$K$ neighbors')\n",
    "plt.ylabel('Cross-validation accuracy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question**: Why does the accuracy decrease when considering more neighbors? What is your expectation for the asymptotic behavior of this accuracy?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn also has a utility that will automatically do the optimization to find the model parameter that produces the highest-accuracy predictions (see [the documentation for some tips](http://scikit-learn.org/stable/modules/grid_search.html#grid-search-tips)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "knn_clf = KNeighborsClassifier()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercise__: Split the above datasets into train and test sets. We'll need to split both the features (`X`) and the labels (`y`) in the same way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = np.random.choice(len(X_spec_mini), size=len(X_spec_mini), replace=False)\n",
    "\n",
    "# 85% training data, 15% test data\n",
    "split = int(0.85 * len(X_spec_mini))\n",
    "\n",
    "train_X = X_spec_mini[idx[:split]]\n",
    "train_y = y[idx[:split]]\n",
    "\n",
    "test_X = X_spec_mini[idx[split:]]\n",
    "test_y = y[idx[split:]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use GridSearchCV to find the optimal number of neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'n_neighbors': np.arange(2, 30, 1)}\n",
    "cv_clf = GridSearchCV(knn_clf, param_grid=params)\n",
    "cv_clf.fit(X_spec_mini, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_clf.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_clf.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Naive Bayes\n",
    "\n",
    "With Gaussians!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.naive_bayes import GaussianNB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = GaussianNB()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X_spec_mini, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = clf.predict(X_spec_mini)\n",
    "accuracy_score(y, y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question**: Why does Gaussian Naive Bayes perform so much worse than nearest neighbors?\n",
    "\n",
    "A demonstration with a pathological data set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = np.random.multivariate_normal([1, 1], np.array([[0.5, 0.499], [0.499, 0.5]])**2, size=1000)\n",
    "x2 = np.random.multivariate_normal([1, 1], np.array([[0.5**2, -0.499**2], [-0.499**2, 0.5**2]]), size=1000)\n",
    "X_patho = np.vstack((x1, x2))\n",
    "y_patho = np.concatenate((np.zeros(len(x1)), np.ones(len(x2))))\n",
    "\n",
    "# shuffle\n",
    "idx = np.random.choice(X_patho.shape[0], size=X_patho.shape[0], replace=False)\n",
    "X_patho = X_patho[idx]\n",
    "y_patho = y_patho[idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
    "ax.scatter(x1[:, 0], x1[:, 1])\n",
    "ax.scatter(x2[:, 0], x2[:, 1])\n",
    "ax.set_xlim(-1, 3)\n",
    "ax.set_ylim(-1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gnb = GaussianNB()\n",
    "gnb.fit(X_patho[:-128], y_patho[:-128])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_patho_pred = gnb.predict(X_patho[-128:])\n",
    "accuracy_score(y_patho[-128:], y_patho_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Decision trees\n",
    "\n",
    "Make binary decisions for each feature with decision boundaries perpendicular to feature axes. Need to specify the depth/height of the tree, $D$ - this partitions the feature space into $2^D$ regions. Can set $D$ with cross-validation!\n",
    "\n",
    "Can improve bias and variance using boosting, bagging.\n",
    "\n",
    "For now, we'll just use a vanilla decision tree, but using the full spectrum as our feature set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = DecisionTreeClassifier(max_depth=8, \n",
    "                             min_samples_split=2, \n",
    "                             min_samples_leaf=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X_spec[:-1024], y[:-1024])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = clf.predict(X_spec[-1024:])\n",
    "accuracy_score(y[-1024:], y_pred[-1024:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Use cross-validation to find the optimal `max_depth`. How does this compare to nearest neighbors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}