Session 1 - Features_I.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 #1: Feature Engineering I\n",
    "Author: [Eric Nalisnick](http://www.ics.uci.edu/~enalisni/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminaries\n",
    "\n",
    "1.  If you haven't done so already, download and **install the [Anaconda Scientific Python Distribution version 2.7](https://store.continuum.io/cshop/anaconda/)**.  If it offers to make itself your default Python distribution, allow it.\n",
    "1. Whether you've just installed Anaconda, or you have done so previously, you should now **update Anaconda** to the latest version of the distribution.  It changes a lot so do this today even if you did recently.\n",
    " 1. Open a terminal or command prompt.\n",
    " 1. Type ```conda update conda``` and press enter or return.  Confirm that you'd like it to make any changes that it offers.\n",
    "1. **Download the code repository**.  \n",
    " 1. Go to [bit.ly/uci_predictive2](http://bit.ly/uci_predictive2) and click the \"download zip\" button on the right to download a zip file containing this entire repository.\n",
    " 1. Unzip that file into a directory you know how to find; you'll need it several times throughout the day.  \n",
    "1. **Start an Jupyter notebook server**.\n",
    " 1. Open a terminal and type ```jupyter notebook```.  Navigate to the directory where you unzipped this repository.\n",
    " 1. Open \"Test Notebook.ipynb\".\n",
    " 1. Click \"Cell\" at the top of the opened notebook, followed by \"Run All\" and ensure that 1) there are no errors and that 2) the output from the first cell is the same as that in the second.  If it doesn't match, raise your hand.\n",
    " 1. If everything looks good, close the browser tab containing the test notebook but keep open the tab listing all the other notebooks.\n",
    " \n",
    "### 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",
    "- Feature Transformations\n",
    "    - Standard Normal Transform\n",
    "    - Domain-Specific Transform\n",
    "    - Mystery Transform\n",
    "    \n",
    "- Dimensionality Reduction\n",
    "    - PCA: Model and Learning\n",
    "    - PCA for Images\n",
    "    - PCA for Visualization\n",
    "\n",
    "### References \n",
    "- Chapter 14 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 data preprocessing](http://scikit-learn.org/stable/modules/preprocessing.html#preprocessing)\n",
    "- [SciKit-Learn's documentation on dimensionality reduction](http://scikit-learn.org/stable/modules/decomposition.html#decompositions)\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": [
    "I've created a function that we'll use later to create visualizations.  It is a bit messy and not essential to the material so don't worry about understanding it.  I'll be happy to explain it to anyone interested during a break or after the session.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# Another messy looking function to make pretty plots of basketball courts\n",
    "def visualize_court(log_reg_model, coord_type='cart', court_image = './data/nba/nba_court.jpg'):\n",
    "    two_class_cmap = ListedColormap(['#FFAAAA', '#AAFFAA']) # light red for miss, light green for make\n",
    "    x_min, x_max = 0, 50 #width (feet) of NBA court\n",
    "    y_min, y_max = 0, 47 #length (feet) of NBA half-court\n",
    "    grid_step_size = 0.2\n",
    "    grid_x, grid_y = np.meshgrid(np.arange(x_min, x_max, grid_step_size), np.arange(y_min, y_max, grid_step_size))\n",
    "    features = np.c_[grid_x.ravel(), grid_y.ravel()]\n",
    "    # change coordinate system\n",
    "    if coord_type == 'polar':\n",
    "        features = np.c_[grid_x.ravel(), grid_y.ravel()]\n",
    "        hoop_location = np.array([25., 0.])\n",
    "        features -= hoop_location\n",
    "        dists = np.sqrt(np.sum(features**2, axis=1))\n",
    "        angles = np.arctan2(features[:,1], features[:,0])\n",
    "        features = np.hstack([dists[np.newaxis].T, angles[np.newaxis].T])\n",
    "        \n",
    "    grid_predictions = log_reg_model.predict(features)\n",
    "    grid_predictions = grid_predictions.reshape(grid_x.shape)\n",
    "    fig, ax = plt.subplots()\n",
    "    court_image = plt.imread(court_image)\n",
    "    ax.imshow(court_image, interpolation='bilinear', origin='lower',extent=[x_min,x_max,y_min,y_max])\n",
    "    ax.imshow(grid_predictions, cmap=two_class_cmap, interpolation = 'nearest',\n",
    "              alpha = 0.60, origin='lower',extent=[x_min,x_max,y_min,y_max])\n",
    "    plt.xlim(x_min, x_max)\n",
    "    plt.ylim(y_min, y_max)\n",
    "    plt.title( \"Make / Miss Prediction Boundaries\" )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also 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": [
    "## 1.  Warm-up\n",
    "\n",
    "Let's start with a warm-up exercise.  In the data directory you'll find a dataset of recent movies and their ratings according to several popular websites.  Let's load it with Pandas...  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load a dataset of recent movies and their ratings across several websites\n",
    "movie_data = pd.read_csv('./data/movie_ratings.csv')\n",
    "# reduce it to just the ratings categories\n",
    "movie_data = movie_data[['FILM','RottenTomatoes','RottenTomatoes_User','Metacritic','Metacritic_User','Fandango_Ratingvalue', 'IMDB']]\n",
    "movie_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "movie_data.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression Review\n",
    "\n",
    "_**Data**_\n",
    "\n",
    "We observe pairs $(\\mathbf{x}_{i},y_{i})$ where\n",
    "\\begin{eqnarray*}\n",
    "y_{i} \\in \\{ 0, 1\\} &:& \\mbox{class label} \\\\\n",
    "\\mathbf{x}_{i} = (x_{i,1}, \\dots, x_{i,D}) &:& \\mbox{set of $D$ explanatory variables (aka features)} \n",
    "\\end{eqnarray*}\n",
    "\n",
    "_** Parameters**_\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\mathbf{\\beta}^{T} = (\\beta_{0}, \\dots, \\beta_{D}) : \\mbox{values encoding the relationship between the features and label}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "_** Transformation Function **_\n",
    "\n",
    "\\begin{equation*}\n",
    "f(z_{i}=\\mathbf{x}_{i} \\mathbf{\\beta} ) = (1+e^{-\\mathbf{x}_{i} \\mathbf{\\beta} })^{-1}\n",
    "\\end{equation*}\n",
    "\n",
    "_**Error Function**_\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\mathcal{L} = \\sum_{i=1}^{N} -y_{i} \\log f(\\mathbf{x}_{i} \\mathbf{\\beta} ) - (1-y_{i}) \\log (1-f(\\mathbf{x}_{i} \\mathbf{\\beta} ))\n",
    "\\end{eqnarray*}\n",
    "\n",
    "_** Learning $\\beta$ **_\n",
    "- Randomly initialize $\\beta$\n",
    "- Until $\\alpha || \\nabla \\mathcal{L} || < tol $:\n",
    "    - $\\mathbf{\\beta}_{t+1} = \\mathbf{\\beta}_{t} - \\alpha \\nabla_{\\mathbf{\\beta}} \\mathcal{L}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\">STUDENT ACTIVITY (10 MINS)</span> \n",
    "Let's run a logistic regression classifier on this data via SciKit-Learn.  If you need a refresher, [check out the notebook from the first course](https://github.com/UCIDataScienceInitiative/PredictiveModeling_withPython/blob/master/Solutions/Session%203%20-%20Solutions.ipynb) and the [SciKit-Learn documentation on logistic regression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html).  The goal is to predict if the *IMDB* rating is under or over 7/10, using the other ratings as features.  I've started the code.  You just need to fill-in the lines for training and computing the error.  Note there is no test set yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# set the random number generator for reproducability\n",
    "np.random.seed(123)\n",
    "\n",
    "# let's try to predict the IMDB rating from the others\n",
    "features = movie_data[['RottenTomatoes','RottenTomatoes_User','Metacritic','Metacritic_User','Fandango_Ratingvalue']].as_matrix()\n",
    "# create classes: more or less that 7/10 rating\n",
    "labels = (movie_data['IMDB'] >= 7.).astype('int').tolist()\n",
    "shuffle_in_unison(features, labels)\n",
    "\n",
    "### Your Code Goes Here ###\n",
    "\n",
    "# initialize and train a logistic regression model\n",
    "\n",
    "# compute error on training data\n",
    "\n",
    "###########################\n",
    "\n",
    "print \"Classification error on training set: %.2f%%\" %(train_error_rate*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %((sum(labels)*100.)/len(labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.  Feature Transformations\n",
    "Good features are crucial for training well-performing classifiers: 'garbage in, garbage out.'  In this section we introduce several transformations that are commonly applied to data as a preprocessing step before training a classifier.\n",
    "\n",
    "### 2.1 Normal Standardization\n",
    "Recall the formula of the standard linear model: $$\\hat Y = f(\\beta^{T} \\mathbf{X}) $$ where $\\hat Y$ is the predictions, $f(\\cdot)$ is the transformation function, $\\beta$ is the weights (parameters), and $X$ is the $N \\times D$ matrix of features.  For simplicity, assume there are just two features: $$ \\beta^{T} \\mathbf{x}_{i} = \\beta_{1}x_{i,1} + \\beta_{2}x_{i,2}.$$  Usually $x_{i,1}$ and $x_{i,2}$ will be measured in different units.  For instance, in the movie ratings data, the Rotten Tomatoes dimension is on a $0-100$ scale and the Fandango ratings are on $0-5$.  The difference in scale causes one dimension to dominate the inner product.  Linear models can learn to cope with this imbalance by changing the scales of the weights accordingly, but this makes optimization harder because gradient steps are unequal across dimensions.\n",
    "\n",
    "One way to get rid of hetergeneous scales is to standardize the data so that the values in each dimension are distributed according to the standard Normal distribution.  In math, this means we'll transform the data like so: $$\\mathbf{X}_{std} = \\frac{\\mathbf{X} - \\boldsymbol{\\mu}_{X}}{\\boldsymbol{\\sigma}_{X}}. $$  This is also called 'z-score scaling.'  Let's examine the affect of this transformation on training error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# perform z-score scaling\n",
    "features_mu = np.mean(features, axis=0)\n",
    "features_sigma = np.std(features, axis=0)\n",
    "std_features = (features - features_mu)/features_sigma\n",
    "\n",
    "# re-train model\n",
    "lm.fit(std_features, labels)\n",
    "\n",
    "### compute error on training data\n",
    "predictions = lm.predict(std_features)\n",
    "print \"Classification error on training set: %.3f%%\" %(calc_classification_error(predictions, labels)*100)\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %((sum(labels)*100.)/len(labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standard Normal scaling is a common and usually default first step, especially when you know the data in measured in different units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Domain-Specific Transformations\n",
    "Sometimes the data calls for a specific transformation.  We'll demonstrate this on the NBA dataset used in our first workshop.  Let's load it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nba_shot_data = pd.read_csv('./data/nba/NBA_xy_features.csv')\n",
    "nba_shot_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's run logistic regression on the data just as we did before..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# split data into train and test\n",
    "train_set_size = int(.80*len(nba_shot_data))\n",
    "train_features = nba_shot_data.ix[:train_set_size,['x_Coordinate','y_Coordinate']].as_matrix()\n",
    "test_features = nba_shot_data.ix[train_set_size:,['x_Coordinate','y_Coordinate']].as_matrix()\n",
    "train_class_labels = nba_shot_data.ix[:train_set_size,['shot_outcome']].as_matrix()\n",
    "test_class_labels = nba_shot_data.ix[train_set_size:,['shot_outcome']].as_matrix()\n",
    "\n",
    "#Train logistic regression model\n",
    "start_time = time.time()\n",
    "lm.fit(train_features, np.ravel(train_class_labels))\n",
    "end_time = time.time()\n",
    "print \"Training ended after %.2f seconds.\" %(end_time-start_time)\n",
    "\n",
    "# compute the classification error on training data\n",
    "predictions = lm.predict(test_features)\n",
    "print \"Classification Error on the Test Set: %.2f%%\" %(calc_classification_error(predictions, np.array(test_class_labels)) * 100)\n",
    "\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %(np.sum(test_class_labels)/len(test_class_labels)*100)\n",
    "\n",
    "# visualize the boundary on the basketball court\n",
    "visualize_court(lm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's transform the Cartesian coordinates into *polar* coordinates: (x,y) --> (radius, angle)...  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Transform coordinate system\n",
    "\n",
    "# radius coordinate: calculate distance from point to hoop\n",
    "hoop_location = np.array([25.5, 0.])\n",
    "train_features -= hoop_location\n",
    "test_features -= hoop_location\n",
    "train_dists = np.sqrt(np.sum(train_features**2, axis=1))\n",
    "test_dists = np.sqrt(np.sum(test_features**2, axis=1))\n",
    "\n",
    "# angle coordinate: use arctan2 function\n",
    "train_angles = np.arctan2(train_features[:,1], train_features[:,0])\n",
    "test_angles = np.arctan2(test_features[:,1], test_features[:,0])\n",
    "\n",
    "# combine vectors into polar coordinates\n",
    "polar_train_features = np.hstack([train_dists[np.newaxis].T, train_angles[np.newaxis].T])\n",
    "polar_test_features = np.hstack([test_dists[np.newaxis].T, test_angles[np.newaxis].T])\n",
    "\n",
    "pd.DataFrame(polar_train_features, columns=[\"Radius\",\"Angle\"]).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Train model\n",
    "start_time = time.time()\n",
    "lm.fit(polar_train_features, np.ravel(train_class_labels))\n",
    "end_time = time.time()\n",
    "print \"Training ended after %.2f seconds.\" %(end_time-start_time)\n",
    "\n",
    "# compute the classification error on test data\n",
    "predictions = lm.predict(polar_test_features)\n",
    "print \"Classification Error on the Test Set: %.2f%%\" %(calc_classification_error(predictions, np.array(test_class_labels)) * 100)\n",
    "\n",
    "# compute the baseline error since the classes are imbalanced\n",
    "print \"Baseline Error: %.2f%%\" %(np.sum(test_class_labels)/len(test_class_labels)*100)\n",
    "\n",
    "# visualize the boundary on the basketball court\n",
    "visualize_court(lm, coord_type='polar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\">STUDENT ACTIVITY (10 mins)</span> \n",
    "\n",
    "## 2.3  Mystery Data\n",
    "The data folder contains some mysterious data that can't be modeled well with a linear function.  Running the code below, we see the squared error is over 70.  However, the error can be driven to zero using one of two transformations.  See if you can find one or both.  The transformations are common ones you surely know.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "# load (x,y) where y is the mystery data\n",
    "x = np.arange(0, 30, .2)[np.newaxis].T\n",
    "y = np.load(open('./data/mystery_data.npy','rb'))\n",
    "\n",
    "### transformation goes here ###\n",
    "\n",
    "#x = ?\n",
    "\n",
    "################################\n",
    "\n",
    "# initialize regression model\n",
    "lm = LinearRegression()\n",
    "lm.fit(x,y)\n",
    "y_hat = lm.predict(x)\n",
    "squared_error = np.sum((y - y_hat)**2)\n",
    "\n",
    "if not np.isclose(squared_error,0):\n",
    "    print \"The squared error should be zero!  Yours is %.8f.\" %(squared_error)\n",
    "else:\n",
    "    print \"You found the secret transformation!  Your squared error is %.8f.\" %(squared_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.  Dimensionality Reduction\n",
    "Sometimes the data calls for more aggressive transformations.  High-dimensional data is usually hard to model because classifiers are likely to [overfit](https://en.wikipedia.org/wiki/Overfitting).  Regularization is one way to combat high dimensionality, but often it can not be enough.  This section will cover *dimensionality reduction*--a technique for reducing the number of features while still preserving curcial information.  This is a form of unsupervised learning since we use no class information.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Image Dataset: [Bob Ross](https://en.wikipedia.org/wiki/Bob_Ross) Paintings\n",
    "![alt text](http://i2.wp.com/espnfivethirtyeight.files.wordpress.com/2014/04/bob-ross1.jpg?quality=70&strip=all&w=600&ssl=1 'Bob_Ross_Banner')\n",
    "In this section and throughout the next session, we'll use a dataset of [Bob Ross'](https://en.wikipedia.org/wiki/Bob_Ross) paintings.  Images are a type of data that notoriously have redundant features and whose dimensionality can be reduced significantly, without much loss of information.  We'll explore this phenomenom via 403 $400 \\times 300$ full-color images of natural landscape paintings.\n",
    "\n",
    "Before we load the data, let's take a minute to review how image data is stored on a computer.  Of course, all the computer sees are numbers ranging from 0 to 255.  Each pixel takes on one of these values.  Furthermore, there are three layers to color images--one for red, blue, and green values.  Therefore, the painting we'll examine are represented as $300 \\times 400 \\times 3$-dimensional tensors (multi-dimensional arrays).  This layering is depicted below.\n",
    "\n",
    "![image_data_format](./graphics/image_data_format.png)\n",
    "\n",
    "While images need to be represented with three dimensions to be visualized, the learning algorithms we'll consider don't need any notion of color values so I've already flattened the images into vector form, i.e. to create a matrix of size $403 \\times 360000$.  Let's load the dataset..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# un-zip the paintings file\n",
    "import zipfile\n",
    "zipper = zipfile.ZipFile('./data/bob_ross/bob_ross_paintings.npy.zip')\n",
    "zipper.extractall('./data/bob_ross/')\n",
    "\n",
    "# load the 403 x 360,000 matrix\n",
    "br_paintings = np.load(open('./data/bob_ross/bob_ross_paintings.npy','rb'))\n",
    "print \"Dataset size: %d x %d\"%(br_paintings.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then visualize two of the images..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# subplot containing first image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "br_painting = br_paintings[70,:]\n",
    "ax1.imshow(np.reshape(br_painting, (300, 400, 3)))\n",
    "\n",
    "# subplot containing second image\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "br_painting = br_paintings[33,:]\n",
    "ax2.imshow(np.reshape(br_painting, (300, 400, 3)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 3.2 Principal Component Analysis\n",
    "\n",
    "As we've seen, the dataset has many, many, many more features (columns) than examples (rows).  Simple Lasso or Ridge regularization probably won't be enough to prevent overfitting so we have to do something more drastic.  In this section, we'll cover [Principal Component Analysis](https://en.wikipedia.org/wiki/Principal_component_analysis), a popular technique for reducing the dimensionality of data.\n",
    "\n",
    "#### Unsupervised Learning\n",
    "PCA does not take into consideration labels, only the input features.  We can think of PCA as performing unsupervised 'inverse' prediction.  Our goal is: for a datapoint $\\mathbf{x}_{i}$, find a lower-dimensional representation $\\mathbf{h}_{i}$ such that $\\mathbf{x}_{i}$ can be 'predicted' from $\\mathbf{h}_{i}$ using a linear transformation.  In math, this statement can be written as $$\\mathbf{\\tilde x}_{i} = \\mathbf{h}_{i} \\mathbf{W}^{T} \\text{ where } \\mathbf{h}_{i} = \\mathbf{x}_{i} \\mathbf{W}. $$  $\\mathbf{W}$ is a $D \\times K$ matrix of parameters that need to be learned--much like the $\\beta$ vector in regression models.  $D$ is the dimensionality of the original data, and $K$ is the dimensionality of the compressed representation $\\mathbf{h}_{i}$.  The graphic below reiterates the above described PCA pipline: \n",
    "\n",
    "![pca_pipeline](./graphics/pca_pipeline.png)\n",
    "\n",
    "#### Optimization\n",
    "Having defined the PCA model, we look to write learning as an optimization process.  Recall that we wish to make a reconstruction of the data, denoted $\\mathbf{\\tilde x}_{i}$, as close as possible to the original input: $$\\mathcal{L} = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{\\tilde x}_{i})^{2}.$$  We can make a substitution for $\\mathbf{\\tilde x}_{i}$ from the equation above: $$ = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{h}_{i}\\mathbf{W}^{T})^{2}.$$  And we can make another substitution for $\\mathbf{h}_{i}$, bringing us to the final form of the loss function: $$ = \\sum_{i=1}^{N} (\\mathbf{x}_{i} - \\mathbf{x}_{i}\\mathbf{W}\\mathbf{W}^{T})^{2}.$$   \n",
    "\n",
    "We could perform gradient descent on $\\mathcal{L}$, just like we do for logistic regression models, but there exists a deterministic solution.  We won't show the derivation here, but you can find it [here](http://stats.stackexchange.com/questions/32174/pca-objective-function-what-is-the-connection-between-maximizing-variance-and-m).  $\\mathbf{W}$ is optimal when it contains the eigenvectors of the data's covariance matrix, and thus we can use a standard eigen decomposition to learn the transform: $$ \\boldsymbol{\\Sigma}_{\\mathbf{X}} = \\mathbf{W} \\boldsymbol{\\Lambda} \\boldsymbol{W}^{T} $$ where $\\boldsymbol{\\Sigma}_{\\mathbf{X}}$ is the data's empirical covariance matrix and $\\boldsymbol{\\Lambda}$ is a diagonal matrix of eigenvalues.  Eigen decompositions can be performed effeciently by any scientific computing library, including [numpy](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig).   \n",
    "\n",
    "#### Intuition\n",
    "The connection to the data's (co-)variance becomes a little more clear when the intuitions behind PCA are examined.  The PCA transformation projects the data onto linear subspaces oriented in the directions of highest variance.  To elaborate, assume the data resides in two dimensions according to the following scatter plot.  The columns of $\\mathbf{W}$--the $K=2$ principal components--would be the green lines below:\n",
    "![alt text](http://weigend.com/files/teaching/stanford/2008/stanford2008.wikispaces.com/file/view/pca_example.gif)\n",
    "'PCA 1st Dimension' is the direction of greatest variance, and if the data is projected down to one dimension, the new representations would be produced by collapsing the data onto that line.\n",
    "\n",
    "#### Principal Component Analysis (PCA) Overview\n",
    "\n",
    "_*Data*_\n",
    "\n",
    "We observe $\\mathbf{x}_{i}$ where\n",
    "\\begin{eqnarray*}\n",
    "\\mathbf{x}_{i} = (x_{i,1}, \\dots, x_{i,D}) &:& \\mbox{set of $D$ explanatory variables (aka features).  No labels.} \n",
    "\\end{eqnarray*}\n",
    "\n",
    "_* Parameters*_\n",
    "\n",
    "$\\mathbf{W}$: Matrix with dimensionality $D \\times K$, where $D$ is the dimensionality of the original data and $K$ the dimensionality of the new features. The matrix encodes the transformation between the original and new feature spaces.\n",
    "\n",
    "_*Error Function*_\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\mathcal{L} = \\sum_{i=1}^{N} ( \\mathbf{x}_{i} - \\mathbf{x}_{i} \\mathbf{W} \\mathbf{W}^{T})^{2}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "#### PCA on Bob Ross dataset\n",
    "\n",
    "Now let's run PCA on the Bob Ross paintings dataset...\n",
    "#### <span style=\"color:red\">Caution: Running PCA on this dataset can take from 30 seconds to several minutes, depending on your computer's processing power.</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA(n_components=400)\n",
    "start_time = time.time()\n",
    "reduced_paintings = pca.fit_transform(br_paintings)\n",
    "end_time = time.time()\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\"%(reduced_paintings.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Let's visualize two of the paintings..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_idx = 70\n",
    "reconstructed_img = pca.inverse_transform(reduced_paintings[img_idx,:])\n",
    "original_img = br_paintings[70,:]\n",
    "\n",
    "# subplot for original image\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.imshow(np.reshape(original_img, (300, 400, 3)))\n",
    "ax1.set_title(\"Original Painting\")\n",
    "\n",
    "# subplot for reconstruction\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.imshow(np.reshape(reconstructed_img, (300, 400, 3)))\n",
    "ax2.set_title(\"Reconstruction\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "We can also visualize the transformation matrix $\\mathbf{W}^{T}$.  It's rows act as 'filters' or 'feature detectors'..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get the transformation matrix\n",
    "transformation_mat = pca.components_ # This is the W^T matrix\n",
    "# two components to show\n",
    "comp1 = 13\n",
    "comp2 = 350\n",
    "\n",
    "# subplot \n",
    "ax1 = plt.subplot(1,2,1)\n",
    "filter1 = transformation_mat[comp1-1,:]\n",
    "ax1.imshow(np.reshape(filter1, (300, 400, 3)))\n",
    "ax1.set_title(\"%dth Principal Component\"%(comp1))\n",
    "\n",
    "# subplot \n",
    "ax2 = plt.subplot(1,2,2)\n",
    "filter2 = transformation_mat[comp2-1,:]\n",
    "ax2.imshow(np.reshape(filter2, (300, 400, 3)))\n",
    "ax2.set_title(\"%dth Principal Component\"%(comp2))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 PCA for Visualization\n",
    "PCA can also be done for visualization purposes.  Let's perform PCA on the movie ratings dataset and see if any semblence of the class structure can be seen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get the movie features\n",
    "movie_features = movie_data[['RottenTomatoes','RottenTomatoes_User','Metacritic','Metacritic_User','Fandango_Ratingvalue']].as_matrix()\n",
    "\n",
    "# perform standard scaling again but via SciKit-Learn\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "z_scaler = StandardScaler()\n",
    "movie_features = z_scaler.fit_transform(movie_features)\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "start_time = time.time()\n",
    "movie_2d_proj = pca.fit_transform(movie_features)\n",
    "end_time = time.time()\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\"%(movie_2d_proj.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "labels = movie_data['FILM'].tolist()\n",
    "classes = movie_data['IMDB'].tolist()\n",
    "\n",
    "# color the points by IMDB ranking\n",
    "labels_to_show = []\n",
    "colors = []\n",
    "for idx, c in enumerate(classes):\n",
    "    if c > 7.25:\n",
    "        colors.append('g')\n",
    "        if c > 8.:\n",
    "            labels_to_show.append(labels[idx])\n",
    "    else:\n",
    "        colors.append('r')\n",
    "        if c < 4.75:\n",
    "            labels_to_show.append(labels[idx])\n",
    "\n",
    "# plot data\n",
    "plt.scatter(movie_2d_proj[:, 0], movie_2d_proj[:, 1], marker = 'o', c = colors, s = 150, alpha = .6)\n",
    "\n",
    "# add movie title annotations\n",
    "for label, x, y in zip(labels, movie_2d_proj[:, 0].tolist(), movie_2d_proj[:, 1].tolist()):\n",
    "    if label not in labels_to_show:\n",
    "        continue\n",
    "    if x < 0:\n",
    "        text_x = -20\n",
    "    else:\n",
    "        text_x = 150\n",
    "    plt.annotate(label.decode('utf-8'),xy = (x, y), xytext = (text_x, 40),\n",
    "        textcoords = 'offset points', ha = 'right', va = 'bottom',\n",
    "        arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3, rad=0'),\n",
    "        bbox = dict(boxstyle = 'round,pad=0.5', fc = 'b', alpha = 0.2))\n",
    "    \n",
    "plt.title('PCA Projection of Movies')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\">STUDENT ACTIVITY (until end of session)</span> \n",
    "Your task is to reproduce the above PCA examples on a new dataset of images.  Let's load it..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_olivetti_faces\n",
    "\n",
    "faces_dataset = fetch_olivetti_faces(shuffle=True)\n",
    "faces = faces_dataset.data # 400 flattened 64x64 images\n",
    "person_ids = faces_dataset.target # denotes the identity of person (40 total)\n",
    "\n",
    "print \"Dataset size: %d x %d\" %(faces.shape)\n",
    "print \"And the images look like this...\"\n",
    "plt.imshow(np.reshape(faces[200,:], (64, 64)), cmap='Greys_r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset contains 400 64x64 pixel images of 40 people each exhibiting 10 facial expressions.  The images are in gray-scale, not color, and therefore flattened vectors contain 4096 dimensions.\n",
    "\n",
    "### <span style=\"color:red\">Subtask 1: Run PCA</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# train PCA model on 'faces'\n",
    "\n",
    "###########################\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\"%(faces_reduced.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 2: Reconstruct an image</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Use learned transformation matrix to project back to the original 4096-dimensional space\n",
    "# Remember you need to use np.reshape() \n",
    "\n",
    "###########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your output should look something like what's below (although could be a different face):\n",
    "![reconstruction_solution](./graphics/face_reconstruction_soln.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 3: Visualize one or more components of the transformation matrix (W)</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Now visualize one of the principal components\n",
    "# Again, remember you need to use np.reshape() \n",
    "\n",
    "###########################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your output should look something like what's below (although could have differently ranked components):\n",
    "![reconstruction_solution](./graphics/face_components_soln.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:red\">Subtask 4: Generate a 2D scatter plot</span> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "### Your code goes here ###\n",
    "\n",
    "# Train another PCA model to project the data into two dimensions\n",
    "# Bonus: color the scatter plot according to the person_ids to see if any structure can be seen\n",
    "\n",
    "# Run PCA for 2 components\n",
    "\n",
    "# Generate plot\n",
    "\n",
    "###########################"
   ]
  }
 ],
 "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
}