pylab_2d.ipynb

{
  "metadata": {
    "name": "Simple Plotting of Classifier Behavior"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Simple Plotting of Classifier Behavior"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This example runs a number of classifiers on a simple 2D dataset and plots the\ndecision surface of each classifier.\n\n", 
            "First compose some sample data -- no PyMVPA involved."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "import numpy as np\n\n# set up the labeled data\n# two skewed 2-D distributions\nnum_dat = 200\ndist = 4\n# Absolute max value allowed. Just to assure proper plots\nxyamax = 10\nfeat_pos=np.random.randn(2, num_dat)\nfeat_pos[0, :] *= 2.\nfeat_pos[1, :] *= .5\nfeat_pos[0, :] += dist\nfeat_pos = feat_pos.clip(-xyamax, xyamax)\nfeat_neg=np.random.randn(2, num_dat)\nfeat_neg[0, :] *= .5\nfeat_neg[1, :] *= 2.\nfeat_neg[0, :] -= dist\nfeat_neg = feat_neg.clip(-xyamax, xyamax)\n\n# set up the testing features\nnpoints = 101\nx1 = np.linspace(-xyamax, xyamax, npoints)\nx2 = np.linspace(-xyamax, xyamax, npoints)\nx,y = np.meshgrid(x1, x2);\nfeat_test = np.array((np.ravel(x), np.ravel(y)))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now load PyMVPA and convert the data into a proper\n", 
            "[Dataset](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.html#mvpa2-datasets-base-dataset)."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "from mvpa2.suite import *\n\n# create the pymvpa dataset from the labeled features\npatternsPos = dataset_wizard(samples=feat_pos.T, targets=1)\npatternsNeg = dataset_wizard(samples=feat_neg.T, targets=0)\nds_lin = vstack((patternsPos, patternsNeg))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Let's add another dataset: XOR. This problem is not linear separable\nand therefore need a non-linear classifier to be solved. The dataset is\nprovided by the PyMVPA dataset warehouse."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# 30 samples per condition, SNR 2\nds_nl = pure_multivariate_signal(30, 2)\nl1 = ds_nl.sa['targets'].unique[1]\n\ndatasets = {'linear': ds_lin, 'non-linear': ds_nl}"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This demo utilizes a number of classifiers. The instantiation of a\nclassifier involves almost no runtime costs, so it is easily possible\ncompile a long list, if necessary."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# set up classifiers to try out\nclfs = {\n        'Ridge Regression': RidgeReg(),\n        'Linear SVM': LinearNuSVMC(probability=1,\n                      enable_ca=['probabilities']),\n        'RBF SVM': RbfNuSVMC(probability=1,\n                      enable_ca=['probabilities']),\n        'SMLR': SMLR(lm=0.01),\n        'Logistic Regression': PLR(criterion=0.00001),\n        '3-Nearest-Neighbour': kNN(k=3),\n        '10-Nearest-Neighbour': kNN(k=10),\n        'GNB': GNB(common_variance=True),\n        'GNB(common_variance=False)': GNB(common_variance=False),\n        'LDA': LDA(),\n        'QDA': QDA(),\n        }\n\n# How many rows/columns we need\nnx = int(ceil(np.sqrt(len(clfs))))\nny = int(ceil(len(clfs)/float(nx)))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we are ready to run the classifiers. The following loop trains\nand queries each classifier to finally generate a nice plot showing\nthe decision surface of each individual classifier, both for the linear and\nthe non-linear dataset."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "for id, ds in datasets.iteritems():\n    # loop over classifiers and show how they do\n    fig = 0\n\n    # make a new figure\n    pl.figure(figsize=(nx*4, ny*4))\n\n    print \"Processing %s problem...\" % id\n\n    for c in sorted(clfs):\n        # tell which one we are doing\n        print \"Running %s classifier...\" % (c)\n\n        # make a new subplot for each classifier\n        fig += 1\n        pl.subplot(ny, nx, fig)\n\n        # select the clasifier\n        clf = clfs[c]\n\n        # enable saving of the estimates used for the prediction\n        clf.ca.enable('estimates')\n\n        # train with the known points\n        clf.train(ds)\n\n        # run the predictions on the test values\n        pre = clf.predict(feat_test.T)\n\n        # if ridge, use the prediction, otherwise use the values\n        if c == 'Ridge Regression':\n            # use the prediction\n            res = np.asarray(pre)\n        elif 'Nearest-Ne' in c:\n            # Use the dictionaries with votes\n            res = np.array([e[l1] for e in clf.ca.estimates]) \\\n                  / np.sum([e.values() for e in clf.ca.estimates], axis=1)\n        elif c == 'Logistic Regression':\n            # get out the values used for the prediction\n            res = np.asarray(clf.ca.estimates)\n        elif c in ['SMLR']:\n            res = np.asarray(clf.ca.estimates[:, 1])\n        elif c in ['LDA', 'QDA'] or c.startswith('GNB'):\n            # Since probabilities are logprobs -- just for\n            # visualization of trade-off just plot relative\n            # \"trade-off\" which determines decision boundaries if an\n            # alternative log-odd value was chosen for a cutoff\n            res = np.asarray(clf.ca.estimates[:, 1]\n                             - clf.ca.estimates[:, 0])\n            # Scale and position around 0.5\n            res = 0.5 + res/max(np.abs(res))\n        else:\n            # get the probabilities from the svm\n            res = np.asarray([(q[1][1] - q[1][0] + 1) / 2\n                    for q in clf.ca.probabilities])\n\n        # reshape the results\n        z = np.asarray(res).reshape((npoints, npoints))\n\n        # plot the predictions\n        pl.pcolor(x, y, z, shading='interp')\n        pl.clim(0, 1)\n        pl.colorbar()\n        # plot decision surfaces at few levels to emphasize the\n        # topology\n        pl.contour(x, y, z, [0.1, 0.4, 0.5, 0.6, 0.9],\n                   linestyles=['dotted', 'dashed', 'solid', 'dashed', 'dotted'],\n                   linewidths=1, colors='black', hold=True)\n\n        # plot the training points\n        pl.plot(ds.samples[ds.targets == 1, 0],\n               ds.samples[ds.targets == 1, 1],\n               \"r.\")\n        pl.plot(ds.samples[ds.targets == 0, 0],\n               ds.samples[ds.targets == 0, 1],\n               \"b.\")\n\n        pl.axis('tight')\n        # add the title\n        pl.title(c)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }
      ], 
      "metadata": {}
    }
  ]
}