ml-basis-expansion-101/Basis Expansion.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b8a3206d7df375d88d2ccbc5123abca03ea42d8513444b68a23f253f854f9dfb"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt \n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Machine Learning, Basis Expansion, Kernels(1) and Kernels(2) \n",
      "\n",
      "2015-February-13\n",
      "Scott Hendrickson\n",
      "\n",
      "Attempt to draw a conceptual line connecting some machine learning techniques we have been using and talking about recently, a short discussion about Kernel choices for SVM and finding non-linear relationships in data..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Logistic Regression and SVM classifiers are linear\n",
      "\n",
      "These methods formulate an optimization problem that locates separating hyperplanes.\n",
      "\n",
      "###Logistic Regression\n",
      "\n",
      "[Following DrJ from last week] Cost function:\n",
      "\n",
      "$ C = - \\left( y \\: \\text{log}(h_{\\theta}(x)) + (1-y) \\: \\text{log}(1-h_{\\theta}) \\right) $\n",
      "\n",
      "The hypothesis function, $h$:\n",
      "\n",
      "$ h_{\\theta}(x) = \\frac{1}{1 + e^{-z}}$, where $ \\theta^T x = z $. \n",
      "\n",
      "This function takes an arbitrary input vector $x$ and creates a linear combination with the model vector $ \\theta $ to create an output on [0,1]."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sigmoid(z, z0=0, a=1):\n",
      "    return 1./(1. + np.exp(-a*(z-z0)))\n",
      "\n",
      "z = np.linspace(-8,8,50)\n",
      "plt.plot(z,sigmoid(z), color=\"magenta\")\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The common cost function formulation for Logistic Regression in terms of z:\n",
      "\n",
      "$ C = -y \\: \\text{log} \\frac{1}{1+exp(-z)} - (1-y) \\: \\text{log} \\left( 1-\\frac{1}{1+exp(-z)} \\right) $\n",
      "\n",
      "Five points to notice about the above forumlation:\n",
      "\n",
      "0. Optimizating $ C $ determines our model parameters $ \\theta $ from features $ x $ and targets $ y $.\n",
      "\n",
      "1. z is a scalar $ \\theta^T x = z $\n",
      "\n",
      "2. $ \\theta^T \\perp x $ implies $ z = 0 $ and this means $ h_{\\theta}(x) = \\frac{1}{1 + e^{0}} = \\frac{1}{2} $. This condition defines the direction of the separating hyperplane.\n",
      "\n",
      "3. The hyperplane passes through the origin.\n",
      "\n",
      "4. $ \\theta^T x = z $ is linear in features $x$\n",
      "\n",
      "Let's see it in action."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(n=100, sigma=1, mu=1):\n",
      "    # random normal vector\n",
      "    return sigma * np.random.randn(n) + mu\n",
      "    \n",
      "plt.scatter(f(),f())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Separate 2 noisy groups"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# make some 2d data\n",
      "x = f(200)\n",
      "y = np.append(f(100, mu=3), f(100, mu=0))\n",
      "t = np.append(np.ones(100), np.zeros(100))\n",
      "# plot it\n",
      "plt.scatter(x,y, s=30, c=t, cmap=plt.cm.Paired)\n",
      "X = np.array([x,y]).transpose()\n",
      "\n",
      "# meshgrid for shading and contours\n",
      "xx, yy = np.meshgrid(np.linspace(-4, 4, 500),\n",
      "                     np.linspace(-4, 6, 500))\n",
      "\n",
      "# make a model\n",
      "from sklearn.linear_model import LogisticRegression\n",
      "m = LogisticRegression()\n",
      "m.fit(X,t)\n",
      "# get the decision boundaries\n",
      "Z = m.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
      "Z = Z.reshape(xx.shape)\n",
      "plt.imshow(Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "\n",
      "# Prediction on the training set\n",
      "p = m.predict(X)\n",
      "#print(p)\n",
      "print(\"The hyperplane predicts {:.2%} of the training data.\".format(m.score(X,t)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###What about when a plane can't do the job?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# meshgrid for contours and shading\n",
      "xx, yy = np.meshgrid(np.linspace(-15, 15, 700),\n",
      "                     np.linspace(-15, 15, 700))\n",
      "\n",
      "# make some data\n",
      "x = np.append(f(100), f(100, sigma=5))\n",
      "y = np.append(f(100, mu=0), f(100, sigma=5))\n",
      "t = np.append(np.ones(100), np.zeros(100))\n",
      "plt.scatter(x,y, s=30, c=t, cmap=plt.cm.Paired)\n",
      "X = np.array([x,y]).transpose()\n",
      "# model\n",
      "m = LogisticRegression()\n",
      "m.fit(X,t)\n",
      "\n",
      "Z = m.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
      "Z = Z.reshape(xx.shape)\n",
      "plt.imshow(Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "# Prediction score\n",
      "p = m.predict(X)\n",
      "#print(p)\n",
      "print(\"The hyperplane predicts {:.2%} of the training data.\".format(m.score(X,t)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Maybe relax (3) and (4) a bit...\n",
      "\n",
      "We can make this learning method a little more flexible by relaxing (3) and (4).  As long as we don't change conditions (1-3), the Logistic Regression strategy still works (why?).\n",
      "\n",
      "Allow arbitrary translation from the origin by a constant and determined the constant as part of the model:\n",
      "\n",
      "$ z = \\theta^T x $\n",
      "\n",
      "$ \\theta = \\left( \\begin{array}{ccc}\n",
      "\\theta_1 \\\\\n",
      "\\theta_2 \\\\\n",
      "\\ldots \\\\\n",
      "\\theta_n \\\\\n",
      "\\theta_{\\text translate} \\end{array} \\right)$\n",
      "\n",
      "$ x = \\left( \\begin{array}{ccc}\n",
      "x_1 \\\\\n",
      "x_2 \\\\\n",
      "\\ldots \\\\\n",
      "x_n \\\\\n",
      "1 \\end{array} \\right)$\n",
      "\n",
      "$ z = \\theta^T x = \\theta_1 x_1 + \\theta_2 x_2 + \\ldots + \\theta_{\\text translate} $\n",
      "\n",
      "Or factoring out $\\theta_1$,\n",
      "\n",
      "$ z = \\theta^T x = \\theta_1 (x_1 + \\theta_2^\\prime x_2 + \\ldots + \\theta_{\\text translate}^\\prime) $\n",
      "\n",
      "where $ \\theta_2^\\prime = \\frac{\\theta_2}{\\theta_1} $ etc.\n",
      "\n",
      "Now, the scaling of the transition region, the hyperplane location and the hyperplane direction are determined by the model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# plot it in 1 dim\n",
      "plt.plot(z,sigmoid(z,2,3), color=\"cyan\", label=\"compress, shift\")\n",
      "plt.plot(z,sigmoid(z,-2), color=\"yellow\", label=\"shift\")\n",
      "plt.plot(z,sigmoid(z,0,0.5), color=\"green\", label=\"expand\")\n",
      "plt.legend(loc=2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###A second perspective...\n",
      "\n",
      "The strategy above for relaxing (4) can be seen as finding $ z(x) $ by embedding the original feature vector, $x$, in a space with 1 additional dimension. If the originial feature vectors has dimension $n$, then the new features space has dimension, $n+1$.\n",
      "\n",
      "We might have formulated the problem as a translation $x \\rightarrow x - x^0$, and expansion $ x - x^0 \\rightarrow \\alpha (x - x^0) $ giving the slightly less general formulation,\n",
      "\n",
      "$ z = \\theta^T \\alpha (x - x^0) = \\theta_1 \\alpha_1 (x_1 - x^0_1) + \\theta_2 \\alpha_2 (x_2 - x^0_2) + \\ldots + \\theta_n \\alpha_n (x_n - x^0_n) $\n",
      "\n",
      "Because we are going to let the optimizer find the model parameters $\\theta$, we usually absorb the $\\alpha$s into the $\\theta$s. In this limited formulation, $ z(x) $ is only a function of the distance to the hyperplane.\n",
      "\n",
      "This discussion suggests that one strategy for expanding the capabilties of our classifier is to think up favorable coordinate transformations. When we imagine a useful coordinate trasformation, we calculate $z$, \n",
      "\n",
      "$ z = \\theta^T \\mathbf{M} x = \\theta^T x^\\prime $.\n",
      "\n",
      "I propose to...Think up embedding/coordinate transformations that make the problem that is essentially a non-linear relationship in the original feature space into a linear relationship in some transformed feature space.\n",
      "\n",
      "Then our Logistic Regression will work to build models of non-linear relationships. Hold that thought."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Might inspire some creative geometric thinking, but...\n",
      "\n",
      "The question is \"In what coordinate system is this essentially non-linear relationship a linear relationship?\" Thinking up the right coordinate transformations from this geometric point of view might be practically impossible! "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### How about relaxing (5) directly?\n",
      "\n",
      "The \"absorb the $\\alpha$s\" comment above might push us to different interesting idea. What if we could find a useful representation of $z$ that:\n",
      "\n",
      "1. Is non-linear in input features $x$\n",
      "2. Lets the optimizer sort out the relative strenghts and types of the non-linearity at learning time\n",
      "\n",
      "From ESL,\n",
      "\n",
      "    The core idea of [basis expansion] is to augment/replace the vectors of inputs X with additional variables, which are transformations of X, and then use linear models in the new space of derived input features.  - ESL pg 139\n",
      "\n",
      "This functional approach changes our strategy toward finding $z$ as a functional expansion in $x$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Basis Expansion\n",
      "\n",
      "Maybe we can try a polynomial expansion of $z$? They are easy to write down and we are pretty confident you can represent any form of $z$ we need.\n",
      "\n",
      "$ z(x) = \\Sigma^\\infty_{m=1}\\beta_m x^n $\n",
      "\n",
      "Done. We can now represent a huge range of non-linear relationships between $z$ and $x$ in a way that allows the Logistic Regression method to work by finding a hyperplane. However, this is going to present some challenges for practical machine learning. We have to determine an inifite number of $\\beta$s in training!\n",
      "\n",
      "We need a finite expansion. \n",
      "\n",
      "One idea to assume the higher-order terms are less important and truncate the series after a couple of terms. Absorb the constants into the $\\theta$s.\n",
      "\n",
      "Try:\n",
      "\n",
      "$ z(x) = \\theta_1 x^2 + \\theta_2 x + \\theta_3 $"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# meshgrid for contours and shading\n",
      "xx, yy = np.meshgrid(np.linspace(-15, 15, 500),\n",
      "                     np.linspace(-15, 15, 500))\n",
      "\n",
      "plt.scatter(x,y, s=30, c=t, cmap=plt.cm.Paired)\n",
      "# quadratic terms, including cross terms between the dimensions x and y\n",
      "x1 = x*x\n",
      "y1 = y*y\n",
      "xy = x*y\n",
      "ext_X = np.array([x, y, x1, y1, xy]).transpose()\n",
      "ext_m = LogisticRegression()\n",
      "ext_m.fit(ext_X,t)\n",
      "\n",
      "gx, gy = xx.ravel(), yy.ravel()\n",
      "\n",
      "Z = ext_m.decision_function(np.array([gx, gy, gx*gx, gy*gy, gx*gy]).transpose())\n",
      "Z = Z.reshape(xx.shape)\n",
      "plt.imshow(Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "\n",
      "p = ext_m.predict(ext_X)\n",
      "#print(p)\n",
      "print(ext_m.score(ext_X,t))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Seems promising--we just found a boundary defined by a conical section. Maybe this is enough. If not, can we generalize the idea? \n",
      "\n",
      "Yes, bring on the basis functions.\n",
      "\n",
      "Rewrite $z$,\n",
      "\n",
      "$ z(x) = \\Sigma^M_{m=1}\\beta_m h_m(x) $\n",
      "\n",
      "Before we apply to machine learning, let's see how the function expansion might look.\n",
      "\n",
      "Let's choose a localized basis function with parameters that define the location and extent. One important example is the Raidal Basis Function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# raidal basis function\n",
      "x = np.linspace(-3,9,60)\n",
      "def b1(x, pars):\n",
      "    return np.exp(-(x-pars[0])*(x-pars[0])/(pars[1]*pars[1]))\n",
      "\n",
      "plt.plot(x, b1(x, (0,1)), color=\"green\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# target function\n",
      "def target(x):\n",
      "    return 0.1*np.power(x,3) - np.power(x,2) + 2*x + 5\n",
      "\n",
      "# point by point difference\n",
      "def c(x, pars):\n",
      "    return np.power(target(x) - pars[2]*b1(x,pars[:2]), 2)\n",
      "\n",
      "# cost function is some of diffs\n",
      "def cost(pars):\n",
      "    return np.sum(c(x,pars))\n",
      "\n",
      "plt.plot(x,target(x)        , color=\"red\")\n",
      "#plt.plot(x,  10*b1(x, (1,1)), color=\"green\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Best fit if we use only one basis function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.optimize import minimize\n",
      "a = minimize(cost,[1,1,3])\n",
      "print(a.x)\n",
      "plt.plot(x,target(x), color=\"red\")\n",
      "plt.plot(x, a.x[2]*b1(x, a.x), color=\"green\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Best fit if we use more than one basis function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# a fun python functional thing will be useful in a minute\n",
      "from functools import partial\n",
      "basetwo = partial(int, base=2)\n",
      "basetwo.__doc__ = 'Convert base 2 string to an int.'\n",
      "basetwo('10010')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To see how additional basis functions improve the fit, change the value of n in the line gx=parital(g,n) below. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def g(x, pars, n=1):\n",
      "    sum = 0.0\n",
      "    for i in range(n):\n",
      "        sum += pars[i*3 + 2]*b1(x, pars[i*3:i*3+2])\n",
      "    return sum\n",
      "\n",
      "from functools import partial\n",
      "plt.plot(x,  target(x), color=\"red\")\n",
      "\n",
      "gx = partial(g, n=4)\n",
      "\n",
      "def c(x, pars):\n",
      "    return np.power(target(x) - gx(x, pars), 2)\n",
      "\n",
      "\n",
      "a = minimize(cost,[ 1,1,3,\n",
      "                   -1,1,1,\n",
      "                   -1,2,1,\n",
      "                    1,2,1])\n",
      "print(a.x)\n",
      "\n",
      "plt.scatter(x, gx(x, a.x), color=\"green\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ok, n=4 gives a nice smooth fit.  Now we are ready to try a machine learning prolem in 2d. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def b2(xx, pars):\n",
      "    # xx is 2d\n",
      "    return np.array([np.exp(-(np.power(x[0]-pars[0],2) +\n",
      "                    np.power(x[1]-pars[1],2))\n",
      "                        /(pars[2]*pars[2])) for x in xx])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# meshgrid\n",
      "xx, yy = np.meshgrid(np.linspace(-6, 6, 500),\n",
      "                     np.linspace(-6, 6, 500))\n",
      "# some data\n",
      "np.random.seed(0)\n",
      "X = np.random.randn(300, 2)\n",
      "t = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)\n",
      "\n",
      "plt.scatter(X[:, 0], X[:, 1], s=30, c=t, cmap=plt.cm.Paired)\n",
      "\n",
      "m = LogisticRegression()\n",
      "m.fit(X,t)\n",
      "\n",
      "Z = m.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
      "Z = Z.reshape(xx.shape)\n",
      "plt.imshow(Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "\n",
      "plt.show()\n",
      "\n",
      "p = m.predict(X)\n",
      "#print(p)\n",
      "print(\"The hyperplane predicts {:.2%} of the training data.\".format(m.score(X,t)))\n",
      "#"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###XOR requires a non-linear model\n",
      "\n",
      "Lets try using the RBF.\n",
      "\n",
      "Two problems, one fundamental, one practical:\n",
      "\n",
      "1. How many basis functions do I need?\n",
      "2. We have non-muliplicative parameters\n",
      "\n",
      "The first is much like our problem with the infinite series.  Remember, we will always use regularization in practical applications, so have a few too many basis functions might still work fine, while having too few will underperform. A parameter search across the number of basis functions might prove useful.\n",
      "\n",
      "The second problem is a practical problem with common libraries. There is nothing in theory that won't let us minimize the cost function over the parameters of the functions, but additional care for stability, convergence etc will be required.\n",
      "\n",
      "Proposes a set of basis functions located in each of the quadrants."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Train on the XOR data using our basis funcitons"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make some basis functions\n",
      "v11 = partial(b2, pars=(1,1,1))\n",
      "v_11 = partial(b2, pars=(-1,1,1))\n",
      "v1_1 = partial(b2, pars=(1,-1,1))\n",
      "v_1_1 = partial(b2, pars=(-1,-1,1))\n",
      "v22 = partial(b2, pars=(2,2,1))\n",
      "v_22 = partial(b2, pars=(-2,2,1))\n",
      "v2_2 = partial(b2, pars=(2,-2,1))\n",
      "v_2_2 = partial(b2, pars=(-2,-2,1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xx, yy = np.meshgrid(np.linspace(-6, 6, 500),\n",
      "                     np.linspace(-6, 6, 500))\n",
      "#\n",
      "plt.scatter(X[:, 0], X[:, 1], s=30, c=t, cmap=plt.cm.Paired)\n",
      "\n",
      "ext_X = np.array([v11(X), v_11(X), v1_1(X), v_1_1(X),\n",
      "                  v22(X), v_22(X), v2_2(X), v_2_2(X)]).transpose()\n",
      "ext_m = LogisticRegression()\n",
      "ext_m.fit(ext_X,t)\n",
      "\n",
      "gr = np.c_[xx.ravel(), yy.ravel()]\n",
      "\n",
      "ext_Z = ext_m.decision_function(np.array([v11(gr)\n",
      "                                , v_11(gr)\n",
      "                                , v1_1(gr)\n",
      "                                , v_1_1(gr)\n",
      "                                , v22(gr)\n",
      "                                , v_22(gr)\n",
      "                                , v2_2(gr)\n",
      "                                , v_2_2(gr)]).transpose())\n",
      "ext_Z = ext_Z.reshape(xx.shape)\n",
      "plt.imshow(ext_Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, ext_Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "\n",
      "p = ext_m.predict(ext_X)\n",
      "#print(p)\n",
      "print(\"The hyperplane predicts {:.2%} of the training data.\".format(ext_m.score(ext_X,t)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Given this particular set of basis functions, we can also reach into the estimator above and find out what the relative weights were in fitting the model to this data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# nb: must be same order as the X array that went into the model fit (ext_X)\n",
      "feature_names = ['v11', 'v_11', 'v1_1', 'v_1_1', 'v22', 'v_22', 'v2_2', 'v_2_2']\n",
      "\n",
      "x = np.arange(len(feature_names))\n",
      "plt.bar(x, ext_m.coef_.ravel() )\n",
      "_ = plt.xticks(x + 0.5, feature_names, rotation=30)\n",
      "_ = plt.ylabel(\"model coefficients\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Kernel Smoothing\n",
      "\n",
      "Turn here for a language-driven diversion.\n",
      "\n",
      "Kernel smoothing is a class of regression techniques that achieve flexibility in estimating the regression function $f(x)$ by fitting a different but simple model separate at each query point $x_0$. This is done by using only the observations close to the target point $x_0$ to fit the simple model, and in such a way that the resulting estimated function $f(x)$ is smooth.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%run kernel.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##Kernel Methods\n",
      "\n",
      "Above, we looked a coordinate transformations as a way to find a space where the problem became linear. For the Logistic Regression, we noted a simple way to forumulate it in terms of only distances to the separating hyperplane. This strategy didn't look very promising because guessing the right coordinate transformation seemed impossible.\n",
      "\n",
      "If we had a machine learning technique that more natrually compatible with formuation of the cost function based only on inner products, we can use different trick to solve the problem.\n",
      "\n",
      "Remember last week when Josh talked about Support Vector Machines.  The cost function,\n",
      "\n",
      "$ \\text{min}_{\\theta} \n",
      "    C \\sum \\limits_{i=1}^{m} \\lbrack \n",
      "    y_i \\text{cost}_1(\\theta^T f) + \n",
      "    (1-y_i) \\text{cost}_0(\\theta^T f) \\rbrack + \n",
      "    \\frac{1}{2} \\sum \\limits_{j=1}^{n} \\theta_j^2 $\n",
      "\n",
      "with f defined by,\n",
      "\n",
      "$ f = \\text{exp} \\left( - \\frac{\\Vert x - l^i \\Vert ^2}{2 \\sigma^2} \\right) $\n",
      "\n",
      "This function depends on the square of the distance between the feature vectors.\n",
      "\n",
      "$ | x - x^\\prime |^2 $\n",
      "\n",
      "This can be written in terms of the inner product of the feature vectors. \n",
      "\n",
      "$ | x - x^\\prime |^2 = x \\cdot x^\\prime $"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###What's the trick?\n",
      "\n",
      "1. Data items are embeded into a vector space called the feature space.\n",
      "2. Linear relations are sought among the images of the data items in the feature space.\n",
      "3. The algorithms (e.g. SVM) are implemented in such a way that the coordinates of teh embedded points are not needed, only the pairwise inner products.\n",
      "4. The pariwise inner products can be computed efficiently direction from the origianl data items usein a kernel function.\n",
      "\n",
      "A kernel function is a function $\\kappa$ the for all $x$, $z$ $\\in$ $X$ satisfies,\n",
      "\n",
      "$ \\kappa(x,z) = <\\phi(x), \\phi(z)> $\n",
      "\n",
      "where $\\phi$ is a mapping from X to an inner product feature space F\n",
      "\n",
      "$\\phi:x \\rightarrow \\phi(x) \\in F$\n",
      "\n",
      "Since the inner product is expressed in the transformed space by the kernel function, we might make some headway by looking at the generic properties of the kernel functions and giving up on actually finding the transformation, $\\phi(x)$.\n",
      "\n",
      "Inner product spaces with the additional properties of being separable and complete are Hilbert spaces.  When $\\kappa$ satisfies  the finitely positive semi-definate property, we refer to it as a Reproducing kernel Hilbert Space (RKHS).\n",
      "\n",
      "It can be shown that once we have kernels $\\kappa_1$ and $\\kappa_2$ we can construct new kernels by:\n",
      "\n",
      "1. $\\kappa(x,z) = \\kappa_1(x,z) + \\kappa_2(x,z)$\n",
      "2. $\\kappa(x,z) = \\alpha \\kappa_1(x,z)$\n",
      "3. $\\kappa(x,z) = \\kappa_1(x,z) \\kappa_2(x,z)$\n",
      "4. $\\kappa(x,z) = f(x) f(z)$\n",
      "5. $\\kappa(x,z) = \\kappa_3(\\phi(x),\\phi(z))$\n",
      "6. $\\kappa(x,z) = x^\\prime B z$\n",
      "\n",
      "Now we just put our efforts into kernel construction.  \n",
      "\n",
      "## RBF and SVM\n",
      "It turns out that a useful kernel is the RBF."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn import svm\n",
      "\n",
      "xx, yy = np.meshgrid(np.linspace(-3, 3, 500),\n",
      "                     np.linspace(-3, 3, 500))\n",
      "np.random.seed(0)\n",
      "X = np.random.randn(300, 2)\n",
      "Y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0)\n",
      "\n",
      "# fit the model\n",
      "clf = svm.NuSVC()\n",
      "clf.fit(X, Y)\n",
      "\n",
      "# plot the decision function for each datapoint on the grid\n",
      "Z = clf.decision_function(np.c_[xx.ravel(), yy.ravel()])\n",
      "Z = Z.reshape(xx.shape)\n",
      "\n",
      "plt.imshow(Z, interpolation='nearest',\n",
      "           extent=(xx.min(), xx.max(), yy.min(), yy.max()), aspect='auto',\n",
      "           origin='lower', cmap=plt.cm.PuOr_r)\n",
      "contours = plt.contour(xx, yy, Z, levels=[0], linewidths=2,\n",
      "                       linetypes='--')\n",
      "plt.scatter(X[:, 0], X[:, 1], s=30, c=Y, cmap=plt.cm.Paired)\n",
      "plt.show()\n",
      "p = clf.predict(X)\n",
      "#print(p)\n",
      "print(\"The hyperplane predicts {:.2%} of the training data.\".format(clf.score(X,t)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Last Note: ESL points out that it is write the solution based no an expansion on RBF basis functions and based on the kernel smoothing propertes sketched above in a way that shows them to be isomorphic.  So the kernel(1) and kernel(2) concepts are not completely disconnected as they might at first appear. (See Chapters 5 & 6.)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}