Session 4 - Ensembles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](http://datascience.uci.edu/wp-content/uploads/sites/2/2014/09/data_science_logo_with_image1.png 'UCI_data_science')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Predictive Modeling with Python \n",
    "## Session #2: Ensembles of Classifiers\n",
    "Author: [Eric Nalisnick](http://www.ics.uci.edu/~enalisni/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schedule for Today\n",
    "\n",
    "|Start Time | Session |\n",
    "|-----------|---------|\n",
    "|8:30am     | Check In|\n",
    "|9:00am     | Feature Engineering I |\n",
    "|10:30am    | Break|\n",
    "|10:45am    | Feature Engineering II |\n",
    "|12:30pm    | Lunch |\n",
    "|1:00pm    | Feature Engineering III |\n",
    "|2:30pm    | Break |\n",
    "|2:45pm    | **Ensembling** |\n",
    "|5:00pm    | End |\n",
    "\n",
    "### Goals of this Lesson\n",
    "- Understand the Bias-Variance Tradeoff\n",
    "   \n",
    "- Ensembling Methods\n",
    "    - Bagging\n",
    "    - Voting \n",
    "    - Stacking\n",
    "\n",
    "### References \n",
    "- Chapter 8 of [*Elements of Statistical Learning* by Hastie, Tibshirani, Friedman](http://statweb.stanford.edu/~tibs/ElemStatLearn/printings/ESLII_print10.pdf)\n",
    "- [A Few Useful Things to Know about Machine Learning](https://homes.cs.washington.edu/~pedrod/papers/cacm12.pdf)\n",
    "- [SciKit-Learn's documentation on ensemble methods](http://scikit-learn.org/stable/modules/ensemble.html)\n",
    "\n",
    "## 0.  Preliminaries\n",
    "First we need to import Numpy, Pandas, MatPlotLib..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import time\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we need functions for shuffling the data and calculating classification errrors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### function for shuffling the data and labels\n",
    "def shuffle_in_unison(features, labels):\n",
    "    rng_state = np.random.get_state()\n",
    "    np.random.shuffle(features)\n",
    "    np.random.set_state(rng_state)\n",
    "    np.random.shuffle(labels)\n",
    "    \n",
    "### calculate classification errors\n",
    "# return a percentage: (number misclassified)/(total number of datapoints)\n",
    "def calc_classification_error(predictions, class_labels):\n",
    "    n = predictions.size\n",
    "    num_of_errors = 0.\n",
    "    for idx in xrange(n):\n",
    "        if (predictions[idx] >= 0.5 and class_labels[idx]==0) or (predictions[idx] < 0.5 and class_labels[idx]==1):\n",
    "            num_of_errors += 1\n",
    "    return num_of_errors/n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0.1  Load the dataset of paintings\n",
    "We are going to use the Bob Ross paintings dataset throughout this session.  Let's again load the data and run PCA..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "# load the 403 x 360,000 matrix\n",
    "br_paintings = np.load(open('./data/bob_ross/bob_ross_paintings.npy','rb'))\n",
    "\n",
    "# perform PCA again\n",
    "pca = PCA(n_components=400)\n",
    "start_time = time.time()\n",
    "pca_paintings = pca.fit_transform(br_paintings)\n",
    "end_time = time.time()\n",
    "\n",
    "# remove the br_paintings from memory\n",
    "br_paintings = None\n",
    "\n",
    "print \"Training took a total of %.2f seconds.\" %(end_time-start_time)\n",
    "print \"Preserved percentage of original variance: %.2f%%\" %(pca.explained_variance_ratio_.sum() * 100) \n",
    "print \"Dataset is now of size: %d x %d\"%(pca_paintings.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to formulate a binary classification problem.  In the data folder there's a file that has labels denoting what is in each painting (tree, mountain, etc.).  Let's load it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "br_labels_data = pd.read_csv('./data/bob_ross/elements-by-episode.csv')\n",
    "br_labels_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Let's make two classes: 1 = 'painting contains hill or mountain', 0 = 'doesn't contain hill/mountain': "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "labels = (br_labels_data['HILLS'] + br_labels_data['MOUNTAIN'] + br_labels_data['MOUNTAINS'] + br_labels_data['SNOWY_MOUNTAIN'] > 0).astype('int8').as_matrix()\n",
    "print \"Contains mountain?: \"+str(bool(labels[5]))\n",
    "recon_img = pca.inverse_transform(pca_paintings[5,:])\n",
    "plt.imshow(np.reshape(recon_img, (300, 400, 3)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Make training and test split..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# set the random number generator for reproducability\n",
    "np.random.seed(182)\n",
    "\n",
    "# shuffle data\n",
    "N = pca_paintings.shape[0]\n",
    "shuffle_in_unison(pca_paintings, labels)\n",
    "\n",
    "# split into train and test sets\n",
    "train_features = pca_paintings[:int(.8*N), :]\n",
    "test_features = pca_paintings[int(.8*N):, :]\n",
    "train_labels = labels[:int(.8*N)]\n",
    "test_labels = labels[int(.8*N):]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 0.2  Run a baseline classifier\n",
    "In order to see the improvements that ensembling provides, let's train a baseline logistic regression classifier for later comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(182)\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# initialize and train a logistic regression model\n",
    "lr_model = LogisticRegression()\n",
    "lr_model.fit(train_features, train_labels)\n",
    "\n",
    "# compute error on test data\n",
    "lr_predictions = lr_model.predict(test_features)\n",
    "one_model_test_error_rate = calc_classification_error(lr_predictions, test_labels)\n",
    "\n",
    "print \"Classification error on test set: %.2f%%\" %(one_model_test_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %((sum(test_labels)*100.)/len(test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.  The Bias-Variance Tradeoff\n",
    "\n",
    "When faced with an important decision, its common to ask multiple people for their advice.  Why should a classification decision be any different?  If computer power is not a limiting factor--which is usually the case--why not train multiple classifiers and combine their predictions?  This is exactly what *ensembling* classifiers does.  In this section we'll cover three methods for combining classifiers: bagging, averaging, and stacking.  But first, let's examine why one classifier is usually not enough.  It can be formalized as a tradeoff between *bias* and *variance*.  \n",
    "\n",
    "Recall the squared loss function: $$\\mathcal{L} = \\sum_{i}^{N} (y_{i} - f(\\mathbf{x}_{i}))^{2}. $$  This loss is over a particular training set {$\\mathbf{X}, \\mathbf{y}$} but we are really interested in the loss over all possible datasets we could have observed, $\\{\\mathbf{X}, \\mathbf{y}\\} \\sim p(\\mathcal{D})$: $$\\mathbb{E}_{p(\\mathcal{D})}[\\mathcal{L}] = \\mathbb{E}_{p(\\mathcal{D})}[(y_{i} - f(\\mathbf{x}_{i}))^{2}]. $$  After some [algebraic manipulations](https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff#Derivation), we can re-write the expected loss as $$\\mathbb{E}_{p(\\mathcal{D})}[\\mathcal{L}] = (f^{*}(\\mathbf{x}_{i}) - \\mathbb{E}[f(\\mathbf{x}_{i})])^{2} + \\text{Var}[f(\\mathbf{x}_{i})] + \\text{constant (error)}. $$  The first term, $(f^{*}(\\mathbf{x}_{i}) - \\mathbb{E}[f(\\mathbf{x}_{i})])^{2}$, is the squared difference between the expected value of the classifier $f$ and the **perfect, true** classifier $f^{*}$.  This difference is known as the *bias* of a classifier.  For instance, a linear model has a strong bias since its functional form is rather simple (unless the optimial classifier is also a linear function).  The second term, $\\text{Var}[f(\\mathbf{x}_{i})]$, is the variance of our classifier.  Basically, this term captures the variability in outputs.  The main point is that if a classifier has *low* bias, meaning it is a very powerful function, then it will usually have high *variance* since this power allows it to generate a wide range of outputs.  And vice versa.  What I just said can be represented graphically as \n",
    "![bias_variance_pic](./graphics/bias-variance.png)\n",
    "Ensembling classifiers all but always produces better performance because it **reduces variance without incurring bias**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.  Bootstrap Aggregating ('Bagging')\n",
    "In section 1, when I mentioned using multiple classifiers, you probably thought I was talking about training and combining several different kinds of classifiers.  We will do that.  But first we'll discuss something simpler: training the same classifier on multiple datasets.  \n",
    "\n",
    "### 2.1 Bootstrap Resampling\n",
    "We (always) want more training data, but unfortunately, it is not available.  We can use the training data we do have and resample it *with replacement* to generate additional 'fake' datasets.  Formally, our original dataset is $$\\{\\mathbf{y}, \\mathbf{X}\\} \\sim p(\\mathcal{D}),$$ where $p(\\mathcal{D})$ is the unknown population distribution.  We then treat the original data as a substitute for the population, writing $$\\{\\mathbf{\\tilde y}, \\mathbf{ \\tilde X}\\}_{1} , \\{\\mathbf{\\tilde y}, \\mathbf{ \\tilde X}\\}_{2}, ... \\sim \\{\\mathbf{y}, \\mathbf{X}\\}.$$  $\\{\\mathbf{\\tilde y}, \\mathbf{ \\tilde X}\\}$ are called bootstrap (re)samples.  Usually, they contain the same number of instances as the original training set.  A diagram showing sampling with replacement is below    \n",
    "![bootstrap_diagram](./graphics/bootstrap_graphic.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's write some Python code to generate Bootstrap samples from a given dataset..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### function for bootstrap resampling\n",
    "def bootstrap_resample(features, labels, n_resamples):\n",
    "    N = features.shape[0]\n",
    "    idxs = np.arange(N)\n",
    "    # numpy's choice() handles the sampling w/ replacement\n",
    "    resampled_idxs = np.random.choice(idxs, size=(N,n_resamples), replace=True)\n",
    "    boot_samps_x = []\n",
    "    boot_samps_y = []\n",
    "    for i in xrange(n_resamples):\n",
    "        boot_samps_x.append(features[resampled_idxs[:,i],:])\n",
    "        boot_samps_y.append(labels[resampled_idxs[:,i]])\n",
    "    return boot_samps_x, boot_samps_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Training on Bootstrap Samples\n",
    "Given the Bootstrap datasets, we next train a classifier on each dataset and then combine their predictions.  These classifiers can all be instances of the model, as is the case in the code below, or different ones.  Once the models are trained, we can combine them in two way: by averaging the probabilities or the predictions.  Both methods are shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fit_and_predict_on_bootstrap_samples(model, bootstrapped_features, bootstrapped_labels, test_features, n_bootstrap_samples):\n",
    "    n_test = test_features.shape[0]\n",
    "    ensemb_probs = np.zeros((n_test,))\n",
    "    ensemb_preds = np.zeros((n_test,))\n",
    "    for idx in xrange(n_bootstrap_samples):\n",
    "        print \"training model #%d\" %(idx+1)\n",
    "        model.fit(bootstrapped_features[idx], bootstrapped_labels[idx])\n",
    "        ensemb_probs += model.predict_proba(test_features)[:,1]\n",
    "        ensemb_preds += model.predict(test_features)\n",
    "    ensemb_probs /= n_bootstrap_samples\n",
    "    ensemb_preds /= n_bootstrap_samples\n",
    "    ensemb_probs = np.around(ensemb_probs)\n",
    "    ensemb_preds = np.around(ensemb_preds)\n",
    "    return ensemb_probs, ensemb_preds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now run the function..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(182)\n",
    "\n",
    "n_bootstrap_samples = 7\n",
    "bootstrapped_features, bootstrapped_labels = bootstrap_resample(train_features, train_labels, n_bootstrap_samples)\n",
    "\n",
    "ensembled_probs, ensembled_preds = fit_and_predict_on_bootstrap_samples(lr_model, bootstrapped_features, bootstrapped_labels, test_features, n_bootstrap_samples)\n",
    "\n",
    "print\n",
    "print \"Averaging probabilities: classification error on test set is %.2f%%\" %(calc_classification_error(ensembled_probs, test_labels)*100)\n",
    "print \"Averaging predictions: classification error on test set is %.2f%%\" %(calc_classification_error(ensembled_preds, test_labels)*100)\n",
    "print\n",
    "print \"One logistic regression model error: %.2f%%\"%(one_model_test_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline error: %.2f%%\" %((sum(test_labels)*100.)/len(test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.  Voting\n",
    "Now we'll consider combining a diverse set of classifers, each trained on an identitical copy of the data.  But first, we need to quickly introduce three new types of classifiers.  Unfortunately, we don't have enough time to cover each in detail and will have to use them as somewhat of a black-box.\n",
    "\n",
    "### 2.1 Overview of Three Classifiers: Decision Tree, k-Nearest Neighbors, and Naive Bayes\n",
    "\n",
    "![three_classifiers](./graphics/classifiers_diagram.png)\n",
    "\n",
    "SciKit-Learn Documentation:\n",
    "- [Decision Tree](http://scikit-learn.org/stable/modules/tree.html#classification)\n",
    "- [k-Nearest Neighbors](http://scikit-learn.org/stable/modules/neighbors.html#nearest-neighbors-classification)\n",
    "- [Naive Bayes](http://scikit-learn.org/stable/modules/naive_bayes.html#gaussian-naive-bayes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Voting via Averaging Predictions\n",
    "\n",
    "$$ \\hat f(\\mathbf{x}_{i}) = \\frac{1}{4} \\hat y_{\\text{DT}} + \\frac{1}{4} \\hat y_{\\text{kNN}} + \\frac{1}{4} \\hat y_{\\text{NB}} + \\frac{1}{4} \\hat y_{\\text{LogReg}}.$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(182)\n",
    "\n",
    "# import the three new classifiers\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "# initialize models\n",
    "d_tree_model = DecisionTreeClassifier()\n",
    "knn_model = KNeighborsClassifier(n_neighbors=5)\n",
    "nb_model = GaussianNB()\n",
    "\n",
    "# fit models\n",
    "d_tree_model.fit(train_features, train_labels)\n",
    "knn_model.fit(train_features, train_labels)\n",
    "nb_model.fit(train_features, train_labels)\n",
    "\n",
    "# predict on test data\n",
    "tree_predictions = d_tree_model.predict(test_features)\n",
    "knn_predictions = knn_model.predict(test_features)\n",
    "nb_predictions = nb_model.predict(test_features)\n",
    "\n",
    "# average predictions\n",
    "# add in the logistic regression predictions calcuated previously\n",
    "avg_predictions = np.around((tree_predictions + knn_predictions + nb_predictions + lr_predictions)/4.)\n",
    "\n",
    "print \"Averaging predictions: classification error on test set is %.2f%%\" %(calc_classification_error(avg_predictions, test_labels)*100)\n",
    "print\n",
    "print \"One logistic regression model error: %.2f%%\"%(one_model_test_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline error: %.2f%%\" %((sum(test_labels)*100.)/len(test_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.  Stacking Models\n",
    "When we performed the averaging above, we used this formula: $$ \\hat f(\\mathbf{x}_{i}) = \\frac{1}{4} \\hat y_{\\text{DT}} + \\frac{1}{4} \\hat y_{\\text{kNN}} + \\frac{1}{4} \\hat y_{\\text{NB}} + \\frac{1}{4} \\hat y_{\\text{LogReg}}.$$  That is, we gave each classifier equal weighting.  While this approach is reasonable, probably it would be better to give an unequal weighting to each classifer, allowing the 'smartest' model to contribute most to the decision.  We can accomplish this by training a second-level logistic regression classifier on the predicted probabilities: $$ \\hat f(\\mathbf{x}_{i}) = \\sigma ( \\alpha_{1} f_{\\text{DT}}(\\mathbf{x}_{i}) + \\alpha_{2}f_{\\text{kNN}}(\\mathbf{x}_{i}) + \\alpha_{3}f_{\\text{NB}}(\\mathbf{x}_{i}) + \\alpha_{4}f_{\\text{LogReg}}(\\mathbf{x}_{i})) \\text{ where } \\sigma (\\cdot) \\text{ is the logistic function}.$$  A depiction of the pipeline is below:\n",
    "![stacking_diagram](./graphics/stacking_diagram.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(182)\n",
    "\n",
    "### TRAINING\n",
    "# calculate probabilities on the training data\n",
    "tree_probs = d_tree_model.predict_proba(train_features)[:,1][np.newaxis].T\n",
    "knn_probs = knn_model.predict_proba(train_features)[:,1][np.newaxis].T\n",
    "nb_probs = nb_model.predict_proba(train_features)[:,1][np.newaxis].T\n",
    "logReg_probs = lr_model.predict_proba(train_features)[:,1][np.newaxis].T\n",
    "\n",
    "# combine into a new 'feature' matrix\n",
    "train_probs_matrix = np.hstack([tree_probs, knn_probs, nb_probs, logReg_probs])\n",
    "\n",
    "# train logistic regression\n",
    "meta_classifier = LogisticRegression()\n",
    "meta_classifier.fit(train_probs_matrix, train_labels)\n",
    "\n",
    "# plot the weights learned for each classifier\n",
    "f,ax = plt.subplots()\n",
    "ticks = np.arange(4)\n",
    "ax.bar(ticks, meta_classifier.coef_[0])\n",
    "ax.set_xticks(ticks+.4)\n",
    "ax.set_xticklabels(['Decision Tree', 'kNN', 'Naive Bayes', 'Log. Regression'])\n",
    "ax.set_title('Weights Learned for Each Classifier')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### TESTING\n",
    "# calculate probabilities on the test data\n",
    "tree_probs = d_tree_model.predict_proba(test_features)[:,1][np.newaxis].T\n",
    "knn_probs = knn_model.predict_proba(test_features)[:,1][np.newaxis].T\n",
    "nb_probs = nb_model.predict_proba(test_features)[:,1][np.newaxis].T\n",
    "logReg_probs = lr_model.predict_proba(test_features)[:,1][np.newaxis].T\n",
    "\n",
    "# combine into a new 'feature' matrix\n",
    "test_probs_matrix = np.hstack([tree_probs, knn_probs, nb_probs, logReg_probs])\n",
    "\n",
    "stacked_predictions = meta_classifier.predict(test_probs_matrix)\n",
    "\n",
    "print \"Averaging predictions: classification error on test set is %.2f%%\" %(calc_classification_error(stacked_predictions, test_labels)*100)\n",
    "print\n",
    "print \"One logistic regression model error: %.2f%%\"%(one_model_test_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline error: %.2f%%\" %((sum(test_labels)*100.)/len(test_labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "train_stacked_predictions = meta_classifier.predict(train_probs_matrix)\n",
    "\n",
    "print \"Stacking train error: %.2f%%\" %(calc_classification_error(train_stacked_predictions, train_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah!  The training error is zero percent--a dead giveaway that the added power of stacking caused us to overfit.  One way of prevent this overfitting is to split the training data and fit the base learners and meta-classifier on different subsets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## <span style=\"color:red\">STUDENT ACTIVITY (until end of session)</span> \n",
    "\n",
    "\n",
    "## 4.  Mini Competition\n",
    "\n",
    "For the remainder of the session, we'll have a mini predictive modeling competition on the [Labeled Faces in the Wild](http://vis-www.cs.umass.edu/lfw/) dataset.  Your task to get as low of an error on the test set as possible by empolying all the techniques we covered today: dimensionality reduction (PCA), Bagging, stacking, etc.  Be mindful to save some data as a validation set.\n",
    "\n",
    "The code below will load the data, display an image, convert it to feature vectors, and train a logistic regression classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_lfw_pairs\n",
    "lfw_train = fetch_lfw_pairs(subset='train')\n",
    "lfw_test = fetch_lfw_pairs(subset='test')\n",
    "lfw_train_pairs = lfw_train['pairs']\n",
    "lfw_train_targets = lfw_train['target']\n",
    "lfw_test_pairs = lfw_test['pairs']\n",
    "lfw_test_targets = lfw_test['target']\n",
    "\n",
    "print \"The training data is of size: %d instances x %d faces x %d pixels x %d pixels\" %(lfw_train_pairs.shape)\n",
    "print \"The test data is of size: %d instances x %d faces x %d pixels x %d pixels\" %(lfw_test_pairs.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's visualize the images..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "face_idx=0\n",
    "\n",
    "# subplot containing first image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(lfw_train_pairs[face_idx,0,:,:],cmap='Greys_r')\n",
    "\n",
    "# subplot containing second image\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(lfw_train_pairs[face_idx,1,:,:],cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a last step of pre-processing, let's flatten the data tensor..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "train_x = np.reshape(lfw_train_pairs, (2200, 5828))\n",
    "train_y = lfw_train_targets\n",
    "test_x = np.reshape(lfw_test_pairs, (1000, 5828))\n",
    "test_y = lfw_test_targets\n",
    "\n",
    "# print the shapes just to check its what we expect\n",
    "print train_x.shape\n",
    "print train_y.shape\n",
    "print test_x.shape\n",
    "print test_y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And just to get you started, here's code to train a logistic regression classifier..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classification error on test set: 47.70%\n",
      "Baseline Error: 50.00%\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(182)\n",
    "\n",
    "# initialize and train a logistic regression model\n",
    "lr_model = LogisticRegression()\n",
    "lr_model.fit(train_x, train_y)\n",
    "\n",
    "# compute error on test data\n",
    "lr_predictions = lr_model.predict(test_x)\n",
    "one_model_test_error_rate = calc_classification_error(lr_predictions, test_y)\n",
    "\n",
    "print \"Classification error on test set: %.2f%%\" %(one_model_test_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %((sum(test_y)*100.)/len(test_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some perspective, state of the art performance on this dataset is [around 5%](http://udrc.eng.ed.ac.uk/sites/udrc.eng.ed.ac.uk/files/publications/class_specific2014.pdf).  [This paper from 1991](https://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf) used PCA and a distance metric to get around 40% error."
   ]
  }
 ],
 "metadata": {
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}