Labs/Lab6_Classification_PCA/classification-and-pca-lab.ipynb

{
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3", 
      "name": "python3", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "3.6.1", 
      "pygments_lexer": "ipython3", 
      "codemirror_mode": {
        "version": 3, 
        "name": "ipython"
      }
    }
  }, 
  "nbformat": 4, 
  "nbformat_minor": 1, 
  "cells": [
    {
      "source": [
        "# Classification and PCA Lab"
      ], 
      "cell_type": "markdown", 
      "metadata": {
        "hide": true
      }
    }, 
    {
      "source": [
        "$$\n", 
        "\\renewcommand{\\like}{{\\cal L}}\n", 
        "\\renewcommand{\\loglike}{{\\ell}}\n", 
        "\\renewcommand{\\err}{{\\cal E}}\n", 
        "\\renewcommand{\\dat}{{\\cal D}}\n", 
        "\\renewcommand{\\hyp}{{\\cal H}}\n", 
        "\\renewcommand{\\Ex}[2]{E_{#1}[#2]}\n", 
        "\\renewcommand{\\x}{{\\mathbf x}}\n", 
        "\\renewcommand{\\v}[1]{{\\mathbf #1}}\n", 
        "$$"
      ], 
      "cell_type": "markdown", 
      "metadata": {
        "variables": {
          "\\cal D": {}, 
          "\\cal E": {}, 
          "\\mathbf #1": {}, 
          "\\mathbf x": {}, 
          "\\cal L": {}, 
          "\\cal H": {}, 
          "\\ell": {}
        }, 
        "hide": true
      }
    }, 
    {
      "execution_count": 1, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline\n", 
        "import numpy as np\n", 
        "import scipy as sp\n", 
        "import matplotlib as mpl\n", 
        "import matplotlib.cm as cm\n", 
        "import matplotlib.pyplot as plt\n", 
        "import pandas as pd\n", 
        "pd.set_option('display.width', 500)\n", 
        "pd.set_option('display.max_columns', 100)\n", 
        "pd.set_option('display.notebook_repr_html', True)\n", 
        "import seaborn.apionly as sns\n", 
        "sns.set_style(\"whitegrid\")\n", 
        "#from PIL import Image"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true, 
        "hide": true
      }
    }, 
    {
      "source": [
        "### Setting up some code\n", 
        "\n", 
        "In doing homework so far you have probably seen strange behaviours when you run and rerun code. This happens in the jupyter notebook because one is generally using global variables, and you might change a value 10 cells down and then rerun a cell 10 cells before. \n", 
        "\n", 
        "To work around such behavior. encapsulate code withon functions and minimize your use of global variables!"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 2, 
      "cell_type": "code", 
      "source": [
        "c0=sns.color_palette()[0]\n", 
        "c1=sns.color_palette()[1]\n", 
        "c2=sns.color_palette()[2]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "source": [
        "A function to plot the points on the training and test set, and the prediction regions associated with a classifier that has 2 features. Adapted from scikit-learn examples."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 3, 
      "cell_type": "code", 
      "source": [
        "from matplotlib.colors import ListedColormap\n", 
        "cmap_light = ListedColormap(['#FFAAAA', '#AAFFAA', '#AAAAFF'])\n", 
        "cmap_bold = ListedColormap(['#FF0000', '#00FF00', '#0000FF'])\n", 
        "cm = plt.cm.RdBu\n", 
        "cm_bright = ListedColormap(['#FF0000', '#0000FF'])\n", 
        "\n", 
        "def points_plot(ax, Xtr, Xte, ytr, yte, clf, mesh=True, colorscale=cmap_light, cdiscrete=cmap_bold, alpha=0.1, psize=10, zfunc=False, predicted=False):\n", 
        "    h = .02\n", 
        "    X=np.concatenate((Xtr, Xte))\n", 
        "    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5\n", 
        "    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5\n", 
        "    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 100),\n", 
        "                         np.linspace(y_min, y_max, 100))\n", 
        "\n", 
        "    #plt.figure(figsize=(10,6))\n", 
        "    if zfunc:\n", 
        "        p0 = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 0]\n", 
        "        p1 = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]\n", 
        "        Z=zfunc(p0, p1)\n", 
        "    else:\n", 
        "        Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])\n", 
        "    ZZ = Z.reshape(xx.shape)\n", 
        "    if mesh:\n", 
        "        plt.pcolormesh(xx, yy, ZZ, cmap=cmap_light, alpha=alpha, axes=ax)\n", 
        "    if predicted:\n", 
        "        showtr = clf.predict(Xtr)\n", 
        "        showte = clf.predict(Xte)\n", 
        "    else:\n", 
        "        showtr = ytr\n", 
        "        showte = yte\n", 
        "    ax.scatter(Xtr[:, 0], Xtr[:, 1], c=showtr-1, cmap=cmap_bold, s=psize, alpha=alpha,edgecolor=\"k\")\n", 
        "    # and testing points\n", 
        "    ax.scatter(Xte[:, 0], Xte[:, 1], c=showte-1, cmap=cmap_bold, alpha=alpha, marker=\"s\", s=psize+10)\n", 
        "    ax.set_xlim(xx.min(), xx.max())\n", 
        "    ax.set_ylim(yy.min(), yy.max())\n", 
        "    return ax,xx,yy"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true, 
        "hide": true
      }
    }, 
    {
      "source": [
        "A function to add contours to such a plot. I use it while showing predictions as opposed to the default \"actual test values\" in `points_plot`."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 4, 
      "cell_type": "code", 
      "source": [
        "def points_plot_prob(ax, Xtr, Xte, ytr, yte, clf, colorscale=cmap_light, cdiscrete=cmap_bold, ccolor=cm, psize=10, alpha=0.1):\n", 
        "    ax,xx,yy = points_plot(ax, Xtr, Xte, ytr, yte, clf, mesh=False, colorscale=colorscale, cdiscrete=cdiscrete, psize=psize, alpha=alpha, predicted=True) \n", 
        "    Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]\n", 
        "    Z = Z.reshape(xx.shape)\n", 
        "    plt.contourf(xx, yy, Z, cmap=ccolor, alpha=.2, axes=ax)\n", 
        "    cs2 = plt.contour(xx, yy, Z, cmap=ccolor, alpha=.6, axes=ax)\n", 
        "    plt.clabel(cs2, fmt = '%5.4f', colors = 'k', fontsize=14, axes=ax)\n", 
        "    return ax "
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true, 
        "hide": true
      }
    }, 
    {
      "source": [
        "### Digits dataset: constructing a classification dataset\n", 
        "\n", 
        "This problem is adapted from Jake Vanderpls's tutorial at pydata:  (http://nbviewer.jupyter.org/github/jakevdp/sklearn_pydata2015/blob/master/notebooks/02.2-Basic-Principles.ipynb)\n", 
        "\n", 
        "The classification problem there is a multiway digit classification problem"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 5, 
      "cell_type": "code", 
      "source": [
        "from sklearn import datasets\n", 
        "digits = datasets.load_digits()\n", 
        "digits.images.shape"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 6, 
      "cell_type": "code", 
      "source": [
        "print(digits.DESCR)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "This code, taken from Jake's notebook above, plots the images against the targets so that we can see what we are dealing with."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 7, 
      "cell_type": "code", 
      "source": [
        "fig, axes = plt.subplots(10, 10, figsize=(8, 8))\n", 
        "fig.subplots_adjust(hspace=0.1, wspace=0.1)\n", 
        "\n", 
        "for i, ax in enumerate(axes.flat):\n", 
        "    ax.imshow(digits.images[i], cmap='binary', interpolation='nearest')\n", 
        "    ax.text(0.05, 0.05, str(digits.target[i]),\n", 
        "            transform=ax.transAxes, color='green')\n", 
        "    ax.set_xticks([])\n", 
        "    ax.set_yticks([])"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 8, 
      "cell_type": "code", 
      "source": [
        "d2d = digits.images.reshape(1797,64,)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 9, 
      "cell_type": "code", 
      "source": [
        "d2d[0].shape, d2d.shape"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 10, 
      "cell_type": "code", 
      "source": [
        "df = pd.DataFrame(d2d)\n", 
        "df['target'] = digits.target\n", 
        "df.head()"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 11, 
      "cell_type": "code", 
      "source": [
        "df.groupby('target').count()"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "> YOUR CODE HERE: To create a stripped down problem for this lab, let us take 2 numbers and try and distinguish them between images. Lets take 8 and 9. make a dtaframe called `dftwo` for this"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 12, 
      "cell_type": "code", 
      "source": [
        "# your code here\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 13, 
      "cell_type": "code", 
      "source": [
        "dftwo.shape"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "### Logistic Regression"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 42, 
      "cell_type": "code", 
      "source": [
        "from sklearn.linear_model import LogisticRegression\n", 
        "from sklearn.model_selection import train_test_split"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 43, 
      "cell_type": "code", 
      "source": [
        "itrain, itest = train_test_split(range(dftwo.shape[0]), train_size=0.6)\n", 
        "set1={}\n", 
        "set1['Xtrain'] = dftwo[list(range(64))].iloc[itrain, :]\n", 
        "set1['Xtest'] = dftwo[list(range(64))].iloc[itest, :]\n", 
        "set1['ytrain'] = dftwo.target.iloc[itrain]==8\n", 
        "set1['ytest'] = dftwo.target.iloc[itest]==8"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "source": [
        "> YOUR TURN HERE: Carry out an unregularized logistic regression and calculate the score on the `set1` test set."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 44, 
      "cell_type": "code", 
      "source": [
        "# your code here\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "### Logistic Regression using Cross Validation and Regularization\n", 
        "\n", 
        "A function to grid search on parameters while doing cross-validation. Note we return the grid-search meta estimator. Be default `GridSearchCV` will refit on the entire training set. Note the use of `scoring`, which will allow for a use of a different scoring function on the cross-validation set than the loss used to train the model on the training set. (Kevin talked about this in class...and the default in `sklearn` is to use the 1-0 loss for scoring on the validation sets, and the log-loss for example in `LogisticRegression`, for training and parameter estimation.\n", 
        "\n", 
        "I keel these separate in my head as **estimation** and **decision** losses. After all, classification requires you to make a decision as to what threshold you will choose."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 45, 
      "cell_type": "code", 
      "source": [
        "from sklearn.model_selection import GridSearchCV\n", 
        "def cv_optimize(clf, parameters, Xtrain, ytrain, n_folds=5, scoring=None):\n", 
        "    if not scoring:\n", 
        "        gs = GridSearchCV(clf, param_grid=parameters, cv=n_folds)\n", 
        "    else:\n", 
        "        gs = GridSearchCV(clf, param_grid=parameters, cv=n_folds, scoring=scoring)\n", 
        "    gs.fit(Xtrain, ytrain)\n", 
        "    print(\"BEST PARAMS\", gs.best_params_)\n", 
        "    return gs"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "source": [
        "`do_classify` is an omnibus function which will take a dataframe, a set of column names to use as features, a name for the target, and do the entire machine learning process for you. For the reason of comparing classifiers, it can take an existing testing and training set as well. If you ask it to, it will standardize as well."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "This was what I had earlier and refactored. What more could you do?"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "def do_classify(clf, parameters, indf, featurenames, targetname, target1val, n_folds=5, standardize=False, train_size=0.6, sets=None, scoring=None):\n", 
        "    if sets:\n", 
        "        Xtrain, Xtest, ytrain, ytest = sets['Xtrain'], sets['Xtest'], sets['ytrain'], sets['ytest']\n", 
        "    else:       \n", 
        "        subdf=indf[featurenames]\n", 
        "        y=(indf[targetname].values==target1val)*1\n", 
        "        # do it stratified? TODO\n", 
        "        itrain, itest = train_test_split(range(subdf.shape[0]), train_size=train_size)\n", 
        "        if standardize:\n", 
        "            dftrain=(subdf.iloc[itrain] - subdf.iloc[itrain].mean())/subdf.iloc[itrain].std()\n", 
        "            dftest=(subdf.iloc[itest] - subdf.iloc[itest].mean())/subdf.iloc[itest].std()\n", 
        "        else:\n", 
        "            dftrain=subdf.iloc[itrain]\n", 
        "            dftest=subdf.iloc[itest]\n", 
        "        Xtrain, Xtest, ytrain, ytest = dftrain.values, dftest.values, y[itrain], y[itest]\n", 
        "    clf = cv_optimize(clf, parameters, Xtrain, ytrain, n_folds=n_folds, scoring=scoring)\n", 
        "    training_score = clf.score(Xtrain, ytrain)\n", 
        "    test_score = clf.score(Xtest, ytest)\n", 
        "    print(\"Score on training data: %0.2f\" % (training_score))\n", 
        "    print(\"Score on test data:     %0.2f\" % (test_score))\n", 
        "    return clf, Xtrain, ytrain, Xtest, ytest"
      ], 
      "cell_type": "raw", 
      "metadata": {
        "collapsed": true, 
        "hide": true
      }
    }, 
    {
      "source": [
        "So we refactor:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 46, 
      "cell_type": "code", 
      "source": [
        "def classify_with_sets(clf, parameters, sets, n_folds = 5, scoring=None):\n", 
        "    Xtrain, Xtest, ytrain, ytest = sets['Xtrain'], sets['Xtest'], sets['ytrain'], sets['ytest']\n", 
        "    gs = cv_optimize(clf, parameters, Xtrain, ytrain, n_folds=n_folds, scoring=scoring)\n", 
        "    training_score = gs.score(Xtrain, ytrain)\n", 
        "    test_score = gs.score(Xtest, ytest)\n", 
        "    print(\"Score on training data: %0.2f\" % (training_score))\n", 
        "    print(\"Score on test data:     %0.2f\" % (test_score))\n", 
        "    return gs"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 47, 
      "cell_type": "code", 
      "source": [
        "def classify_from_dataframe(clf, parameters, indf, featurenames, targetname, target1val, n_folds=5, standardize=False, train_size=0.6, scoring=None):\n", 
        "    subdf=indf[featurenames]\n", 
        "    y=(indf[targetname].values==target1val)*1\n", 
        "    itrain, itest = train_test_split(range(subdf.shape[0]), train_size=train_size)\n", 
        "    inset = {}\n", 
        "    if standardize:\n", 
        "        Xtr = (subdf.iloc[itrain] - subdf.iloc[itrain].mean())/subdf.iloc[itrain].std()\n", 
        "        inset['Xtrain'] = Xtr.values\n", 
        "        Xte = (subdf.iloc[itest] - subdf.iloc[itest].mean())/subdf.iloc[itest].std()\n", 
        "        inset['Xtest'] = Xte.values\n", 
        "    else:\n", 
        "        inset['Xtrain'] = subdf.iloc[itrain].values\n", 
        "        inset['Xtest'] = subdf.iloc[itest].values\n", 
        "    inset['ytrain'] = y[itrain]\n", 
        "    inset['ytest'] = y[itest]\n", 
        "    clf = classify_with_sets(clf, parameters, inset, n_folds=n_folds, scoring=scoring)\n", 
        "    return clf, inset['Xtrain'], inset['ytrain'], inset['Xtest'], inset['ytest']\n", 
        "    \n"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 48, 
      "cell_type": "code", 
      "source": [
        "cvals = [1e-20, 1e-15, 1e-10, 1e-5, 1e-3, 1e-1, 1, 10, 100, 10000, 100000]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 49, 
      "cell_type": "code", 
      "source": [
        "digitstwo_log_set1 = classify_with_sets(\n", 
        "    LogisticRegression(), \n", 
        "    {\"C\": cvals},  \n", 
        "    set1,\n", 
        "n_folds=5) "
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "### The confusion matrix: comparing classifiers"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "We have written two classifiers. A classifier will get some samples right, and some wrong. Generally we see which ones it gets right and which ones it gets wrong on the test set. There,\n", 
        "\n", 
        "- the samples that are +ive and the classifier predicts as +ive are called True Positives (TP)\n", 
        "- the samples that are -ive and the classifier predicts (wrongly) as +ive are called False Positives (FP)\n", 
        "- the samples that are -ive and the classifier predicts as -ive are called True Negatives (TN)\n", 
        "- the samples that are +ive and the classifier predicts as -ive are called False Negatives (FN)\n", 
        "\n", 
        "A classifier produces a confusion matrix from these which lookslike this:\n", 
        "\n", 
        "![hwimages](./images/confusionmatrix.png)\n", 
        "\n", 
        "\n", 
        "IMPORTANT NOTE: In sklearn, to obtain the confusion matrix in the form above, always have the observed `y` first, i.e.: use as `confusion_matrix(y_true, y_pred)`"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 50, 
      "cell_type": "code", 
      "source": [
        "from sklearn.metrics import confusion_matrix\n", 
        "confusion_matrix(set1['ytest'], clf.predict(set1['Xtest']))"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "> YOUR TURN NOW: Calculate the confusion matrix for the regularized logistic regression"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 51, 
      "cell_type": "code", 
      "source": [
        "confusion_matrix(set1['ytest'], digitstwo_log_set1.predict(set1['Xtest']))"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "> YOUR TURN NOW: As an exercise to do this completely with a new train/test split done directly on a dataframe. Call your classifier/estimator `digitstwo_log` and your training/test sets dictionary as `set2`. Compute the confusion matrix for this `set2`"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 52, 
      "cell_type": "code", 
      "source": [
        "# your code here\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 53, 
      "cell_type": "code", 
      "source": [
        "confusion_matrix(set2['ytest'], digitstwo_log_set2.predict(set2['Xtest']))"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "From the department of not-kosher things to do, (why?) we calculate the performance of this classifier on `set1`."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 54, 
      "cell_type": "code", 
      "source": [
        "confusion_matrix(set1['ytest'], digitstwo_log_set2.predict(set1['Xtest']))"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "### Plotting scores against hyper-parameters\n", 
        "\n", 
        "Finally `plot_scores` takes the output of a grid search on one parameter, and plots for you a graph showing the test performance against the parameter. You could augment this with a training set diagram if you like."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 55, 
      "cell_type": "code", 
      "source": [
        "def plot_scores(fitmodel, pname):\n", 
        "    params = [d[pname] for d in fitmodel.cv_results_['params']]\n", 
        "    scores = fitmodel.cv_results_['mean_test_score']\n", 
        "    stds = fitmodel.cv_results_['std_test_score']\n", 
        "    plt.plot(params, scores,'.-');\n", 
        "    plt.fill_between(params, scores - stds, scores + stds, alpha=0.3);"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 58, 
      "cell_type": "code", 
      "source": [
        "plot_scores(digitstwo_log_set2, 'C')\n", 
        "plt.xscale('log')\n", 
        "plt.ylim(0.6,1)"
      ], 
      "outputs": [], 
      "metadata": {
        "scrolled": true
      }
    }, 
    {
      "execution_count": 59, 
      "cell_type": "code", 
      "source": [
        "plot_scores(digitstwo_log_set1, 'C')\n", 
        "plt.xscale('log')\n", 
        "plt.ylim(0.6,1)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "## Feature engineering"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "Our images here are relatively small, but in general you will have as many features as pizels multiplied by the color channels. This is a lot of features! Having too many features can lead to overfitting.\n", 
        "\n", 
        "Indeed, it is possible to have more features than data points. Thus there is a high chance that a few attributes will correlate with $y$ purely coincidentally!\n", 
        "[^Having lots of images, or \"big-data\" helps in combatting overfitting!]\n", 
        "\n", 
        "We need to do something similar to what happened in the regularized regression or classification here! We will engage in some *a-priori* feature selection that will reduce the dimensionality of the problem. The idea we'll use here is something called **Principal Components Analysis**, or PCA.\n", 
        "\n", 
        "PCA is an unsupervized learning technique. The basic idea behind PCA is to rotate the co-ordinate axes of the feature space. We first find the direction in which the data varies the most. We set up one co-ordinate axes along this direction, which is called the first principal component. We then look for a perpendicular direction in which the data varies the second most. This is the second principal component. The diagram illustrates this process. There are as many principal components as the feature dimension: all we have done is a rotation.\n", 
        "\n", 
        "![pcanim](https://i.stack.imgur.com/Q7HIP.gif)\n", 
        "\n", 
        "(diagram taken from http://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues which also has nice discussions)\n", 
        "\n", 
        "How does this then achieve feature selection? We decide on a threshold of variation; once the variation in a particular direction falls below a certain number, we get rid of all the co-ordinate axes after that principal component. For example, if the variation falls below 10% after the third axes, and we decide that 10% is an acceptable cutoff, we remove all dimensions from the fourth dimension onwards. In other words, we took our higher dimensional problem and projected it onto a 3 dimensional **subspace**.\n", 
        "\n", 
        "These two ideas illustrate one of the most important reasons that learning is even feasible: we believe that **most datasets, in either their unsupervized form $\\{\\v{x\\}}$, or their supervized form $\\{y, \\v{x}\\}$, live on a lower dimensional subspace**. If we can find this subspace, we can then hope to find a method which respectively separates or fits the data."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 61, 
      "cell_type": "code", 
      "source": [
        "from sklearn.decomposition import PCA"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "source": [
        "The explained variance ratio `pca.explained_variance_ratio_` tells us how much of the variation in the features is explained by these 60 features. When we sum it up over the features, we see that 94% is explained: good enough to go down to a 60 dimensional space from a 136452 dimensional one!\n", 
        "\n", 
        "We can see the individual variances as we increase the dimensionality:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "The first dimension accounts for 35% of the variation, the second 6%, and it goes steadily down from there.\n", 
        "\n", 
        "Let us create a dataframe with these 16 features labelled pc1,pc2...,pc60 and the labels of the sample:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "Lets see what these principal components look like:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 62, 
      "cell_type": "code", 
      "source": [
        "def normit(a):\n", 
        "    a=(a - a.min())/(a.max() -a.min())\n", 
        "    a=a*256\n", 
        "    return np.round(a)\n", 
        "def getNC(pc, j):\n", 
        "    size=8*8\n", 
        "    g=pc.components_[j][0:size]\n", 
        "    g=normit(g)\n", 
        "    return g\n", 
        "def display_component(pc, j):\n", 
        "    g = getNC(pc,j)\n", 
        "    print(g.shape)\n", 
        "    plt.imshow(g.reshape(8,8))\n", 
        "    plt.xticks([])\n", 
        "    plt.yticks([])"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "source": [
        "You might be a bit confused: we needed to use 16 components to explain 90% of the variation in the features, but only 1 or 2 components to separate checks from dollars? This is because PCA is unsupervised: the only variation we are explaining is the variation in the 64 dimensional feature space. We are not explaining the variation in the $y$ or the label, and it might turn out, as it does in this case, that with the additional information in $y$, the dimensionality needed for classification is much lower.\n", 
        "\n", 
        "We could thus choose just the first few principal components to make our classifier. For the purposes of this lab, since two components can be easily visualized (even though adding some more features may leads to better separability), we'll go with learning a 2-dimensional classifier in the `pc1` and `pc2` dimensions! "
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 63, 
      "cell_type": "code", 
      "source": [
        "pca_digits = PCA(n_components=16)\n", 
        "X2 = pca_digits.fit_transform(dftwo[list(range(64))].values)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 64, 
      "cell_type": "code", 
      "source": [
        "X2"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 65, 
      "cell_type": "code", 
      "source": [
        "print(pca_digits.explained_variance_ratio_.sum())"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 66, 
      "cell_type": "code", 
      "source": [
        "100*pca_digits.explained_variance_ratio_"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 67, 
      "cell_type": "code", 
      "source": [
        "display_component(pca_digits, 0)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 68, 
      "cell_type": "code", 
      "source": [
        "display_component(pca_digits, 1)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 69, 
      "cell_type": "code", 
      "source": [
        "display_component(pca_digits, 2)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 70, 
      "cell_type": "code", 
      "source": [
        "dfpca = pd.DataFrame({\"target\":dftwo.target})\n", 
        "for i in range(pca_digits.explained_variance_ratio_.shape[0]):\n", 
        "    dfpca[\"pc%i\" % (i+1)] = X2[:,i]\n", 
        "dfpca.head()"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 71, 
      "cell_type": "code", 
      "source": [
        "colors = [c0, c2]\n", 
        "for label, color in zip(dfpca['target'].unique(), colors):\n", 
        "    mask = dfpca['target']==label\n", 
        "    plt.scatter(dfpca[mask]['pc1'], dfpca[mask]['pc2'], c=color, label=label, alpha=0.5)\n", 
        "plt.legend()"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "> YOUR CODE NOW: Do a regularized Logistic Regression using the first two principal components. Store the classifier in `digitspca_log2` and the sets in `setf`"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 72, 
      "cell_type": "code", 
      "source": [
        "# your code here\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "We use points plot to see misclassification and the decision boundary:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 73, 
      "cell_type": "code", 
      "source": [
        "plt.figure()\n", 
        "ax=plt.gca()\n", 
        "points_plot(ax, setf['Xtrain'], setf['Xtest'], setf['ytrain'], setf['ytest'], digitspca_log2, alpha=0.5, psize=20);"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "And a probability contour plot to see probabilities"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 75, 
      "cell_type": "code", 
      "source": [
        "plt.figure()\n", 
        "ax=plt.gca()\n", 
        "points_plot_prob(ax,  setf['Xtrain'], setf['Xtest'], setf['ytrain'], setf['ytest'], digitspca_log2, alpha=0.5, psize=20);"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "And a scores plot to see the hyper-parameter landscape."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 76, 
      "cell_type": "code", 
      "source": [
        "plot_scores(digitspca_log2, 'C')\n", 
        "plt.xscale('log')"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "And a confusion matrix..."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 77, 
      "cell_type": "code", 
      "source": [
        "from sklearn.metrics import confusion_matrix\n", 
        "confusion_matrix(setf['ytest'], digitspca_log2.predict(setf['Xtest']), )"
      ], 
      "outputs": [], 
      "metadata": {}
    }
  ]
}