doc/ipython-notebooks/regression/Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### By Saurabh Mahindre - <a href=\"https://github.com/Saurabh7\">github.com/Saurabh7</a> as a part of <a href=\"http://www.google-melange.com/gsoc/project/details/google/gsoc2014/saurabh7/5750085036015616\">Google Summer of Code 2014 project</a> mentored by - Heiko Strathmann - <a href=\"https://github.com/karlnapf\">github.com/karlnapf</a> - <a href=\"http://herrstrathmann.de/\">herrstrathmann.de</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates various regression methods provided in Shogun. Linear models like [Least Square regression](http://en.wikipedia.org/wiki/Ordinary_least_squares), [Ridge regression](http://en.wikipedia.org/wiki/Tikhonov_regularization), [Least Angle regression](http://en.wikipedia.org/wiki/Least-angle_regression), etc. and also [kernel based methods](http://en.wikipedia.org/wiki/Kernel_trick) like Kernel Ridge regression are discussed and applied to toy and real life data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Introduction](#Introduction)\n",
    "2. [Least Squares regression](#Least-Squares-regression)\n",
    "  1. [Prediction using Least Squares](#Prediction-using-Least Squares)\n",
    "  2. [Training and generating weights](#Training-and-generating-weights)\n",
    "3. [Ridge Regression](#Ridge-Regression)\n",
    "  1. [Weights and regularization](#Relationship-between-weights-and-regularization)\n",
    "4. [Least Angle Regression and LASSO](#Least-Angle-Regression-and-LASSO)\n",
    "5. [Kernel Ridge Regression](#Kernel-Ridge-Regression)\n",
    "6. [Support Vector Regression](#Support-Vector-Regression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression is a case of supervised learning where the goal is to learn a mapping from inputs $x\\in\\mathcal{X}$ to outputs $y\\in\\mathcal{Y}$, given a labeled set of input-output pairs $\\mathcal{D} = {(x_i,y_i)}^{\\text N}_{i=1} \\subseteq \\mathcal{X} \\times \\mathcal{Y}$. The response variable $y_i$ is continuous in regression analysis. Regression finds applications in many fields like for predicting stock prices or predicting consumption spending, etc. In linear regression, the mapping is a linear (straight-line) equation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least Squares regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Linear regression model can be defined as $\\text y =$  $\\bf {w} \\cdot \\bf{x} $ $+ b$. Here $\\text y$ is the predicted value, $\\text x$ the independent variable and $\\text w$ the so called weights.</br> We aim to find the linear function (line) that best explains the data, i.e. that minimises some measure of deviation to the training data $\\mathcal{D}$. One such measure is the sum of squared distances. The [Ordinary Least Sqaures method](http://en.wikipedia.org/wiki/Ordinary_least_squares) minimizes the sum of squared distances between the observed responses in the dataset and the responses predicted by the linear approximation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distances called residuals have to minimized. This can be represented as:$$E({\\bf{w}}) = \\sum_{i=1}^N(y_i-{\\bf w}\\cdot {\\bf x}_i)^2$$\n",
    "One can differentiate with respect to $\\bf w$ and equate to zero to determine the $\\bf w$ that minimises $E({\\bf w})$. This leads to solution of the form:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$${\\bf w} = \\left(\\sum_{i=1}^N{\\bf x}_i{\\bf x}_i^T\\right)^{-1}\\left(\\sum_{i=1}^N y_i{\\bf x}_i\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction using Least Squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression using Least squares is demonstrated below on toy data. Shogun provides the tool to do it using [CLeastSquaresRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLeastSquaresRegression.html) class. The data is a straight line with lot of noise and having slope 3. Comparing with the mathematical equation above we thus expect $\\text w$ to be around 3 for a good prediction. Once the data is converted to Shogun format, we are ready to train the machine. To label the training data [CRegressionLabels](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CRegressionLabels.html) are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "%matplotlib inline\n",
    "from cycler import cycler\n",
    "# import all shogun classes\n",
    "from modshogun import *\n",
    "slope = 3\n",
    "\n",
    "X_train = rand(30)*10\n",
    "y_train = slope*(X_train)+random.randn(30)*2+2\n",
    "y_true = slope*(X_train)+2\n",
    "X_test = concatenate((linspace(0,10, 50),X_train))\n",
    "\n",
    "#Convert data to shogun format features\n",
    "feats_train = RealFeatures(X_train.reshape(1,len(X_train)))\n",
    "feats_test = RealFeatures(X_test.reshape(1,len(X_test)))\n",
    "labels_train = RegressionLabels(y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training and generating weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[LeastSquaresRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLeastSquaresRegression.html) has to be initialised with the training features and training labels. Once that is done to learn from data we `train()` it. This also generates the $\\text w$ from the general equation described above. To access $\\text w$ use `get_w()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ls = LeastSquaresRegression(feats_train, labels_train)\n",
    "ls.train()\n",
    "w = ls.get_w()\n",
    "print 'Weights:'\n",
    "print w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This value of $\\text w$ is pretty close to 3, which certifies a pretty good fit for the training data. Now let's `apply` this trained machine to our test data to get the ouput values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "out = ls.apply(feats_test).get_labels()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an aid to visualisation, a plot of the output and also of the residuals is shown. The sum of the squares of these residuals is minimised. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "figure(figsize=(20,5))\n",
    "#Regression and true plot\n",
    "pl1 = subplot(131)\n",
    "title('Regression')\n",
    "_ = plot(X_train,labels_train, 'ro')\n",
    "_ = plot(X_test,out, color='blue')\n",
    "_ = plot(X_train, y_true, color='green')\n",
    "p1 = Rectangle((0, 0), 1, 1, fc=\"r\")\n",
    "p2 = Rectangle((0, 0), 1, 1, fc=\"b\")\n",
    "p3 = Rectangle((0, 0), 1, 1, fc=\"g\")\n",
    "pl1.legend((p1, p2, p3), [\"Training samples\", \"Predicted output\", \"True relationship\"], loc=2)\n",
    "xlabel('Samples (X)', fontsize=12)\n",
    "ylabel('Response variable (Y)', fontsize=12)\n",
    "\n",
    "#plot residues\n",
    "pl2 = subplot(132)\n",
    "title(\"Squared error and output\")\n",
    "_ = plot(X_test,out, linewidth=2)\n",
    "gray()\n",
    "_ = scatter(X_train,labels_train,c=ones(30) ,cmap=gray(), s=40)\n",
    "for i in range(50,80):\n",
    "    plot([X_test[i],X_test[i]],[out[i],y_train[i-50]] , linewidth=2, color='red')\n",
    "p1 = Rectangle((0, 0), 1, 1, fc=\"r\")\n",
    "p2 = Rectangle((0, 0), 1, 1, fc=\"b\")\n",
    "pl2.legend((p1, p2), [\"Error/residuals to be squared\", \"Predicted output\"], loc=2)\n",
    "xlabel('Samples (X)', fontsize=12)\n",
    "ylabel('Response variable (Y)', fontsize=12)\n",
    "jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function we choose should not only best fit the training data but also generalise well. If the coefficients/weights are unconstrained, they are susceptible to high variance and overfitting. To control variance, one has to regularize the coefficients i.e. control how large the coefficients grow. This is what is done in [Ridge regression](http://en.wikipedia.org/wiki/Tikhonov_regularization) which is L2 (sum of squared components of $\\bf w$) regularized form of least squares. A penalty is imposed on the size of coefficients. The error to be minimized is:\n",
    "\n",
    "\n",
    "$$E({\\bf{w}}) = \\sum_{i=1}^N(y_i-{\\bf w}\\cdot {\\bf x}_i)^2 + \\tau||{\\bf w}||^2$$\n",
    "\n",
    "Here $\\tau$ imposes a penalty on the weights.</br>\n",
    "By differentiating the regularised training error and equating to zero, we find the optimal $\\bf w$, given by:\n",
    "\n",
    "$${\\bf w} = \\left(\\tau {\\bf I}+ \\sum_{i=1}^N{\\bf x}_i{\\bf x}_i^T\\right)^{-1}\\left(\\sum_{i=1}^N y_i{\\bf x}_i\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ridge regression can be performed in Shogun using [CLinearRidgeRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLinearRidgeRegression.html) class. It takes the regularization constant $\\tau$ as an additional argument. Let us see the basic regression example solved using the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tau = 0.8\n",
    "rr = LinearRidgeRegression(tau, feats_train, labels_train)\n",
    "rr.train()\n",
    "w = rr.get_w()\n",
    "print w\n",
    "out = rr.apply(feats_test).get_labels()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(20,5))\n",
    "#Regression and true plot\n",
    "pl1 = subplot(131)\n",
    "title('Ridge Regression')\n",
    "_ = plot(X_train,labels_train, 'ro')\n",
    "_ = plot(X_test, out, color='blue')\n",
    "_ = plot(X_train, y_true, color='green')\n",
    "p1 = Rectangle((0, 0), 1, 1, fc=\"r\")\n",
    "p2 = Rectangle((0, 0), 1, 1, fc=\"b\")\n",
    "p3 = Rectangle((0, 0), 1, 1, fc=\"g\")\n",
    "pl1.legend((p1, p2, p3), [\"Training samples\", \"Predicted output\", \"True relationship\"], loc=2)\n",
    "xlabel('Samples (X)', fontsize=12)\n",
    "ylabel('Response variable (Y)', fontsize=12)\n",
    "jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relationship between weights and regularization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prediction in the basic regression example was simliar to that of least squares one. To actually see ridge regression's forte, we analyse how the weights change along with the regularization constant. Data with slightly higher dimensions is sampled to do this because overfitting  is more likely to occur in such data.  Here `set_tau()` method is used to set the necessary parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Generate Data\n",
    "def generate_data(N, D):\n",
    "    w = randn(D,1)\n",
    "    X = zeros((N,D))\n",
    "    y = zeros((N,1))\n",
    "    for i in range(N):\n",
    "        x = randn(1,D)\n",
    "        for j in range(D):\n",
    "            X[i][j] = x[0][j]\n",
    "    y = dot(X,w) + randn(N,1);\n",
    "    y.reshape(N,)\n",
    "    return X, y.T\n",
    "\n",
    "def generate_weights(taus, feats_train, labels_train):\n",
    "    preproc = PruneVarSubMean(True)\n",
    "    preproc.init(feats_train)\n",
    "    feats_train.add_preprocessor(preproc)\n",
    "    feats_train.apply_preprocessor()     \n",
    "    weights = []\n",
    "    rr = LinearRidgeRegression(tau, feats_train, labels_train)\n",
    "    \n",
    "    #vary regularization\n",
    "    for t in taus:\n",
    "        rr.set_tau(t)\n",
    "        rr.train()\n",
    "        weights.append(rr.get_w())\n",
    "    return weights, rr\n",
    "\n",
    "def plot_regularization(taus, weights):     \n",
    "    ax = gca()\n",
    "    ax.set_prop_cycle(cycler('color', ['b', 'r', 'g', 'c', 'k', 'y', 'm']))\n",
    "    ax.plot(taus, weights, linewidth=2)\n",
    "    xlabel('Tau', fontsize=12)\n",
    "    ylabel('Weights', fontsize=12)\n",
    "    ax.set_xscale('log')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean squared error (MSE) of an estimator measures the average of the squares of the errors. [CMeanSquaredError](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CMeanSquaredError.html) class is used to compute the MSE as :\n",
    "\n",
    "$$\\frac{1}{|L|} \\sum_{i=1}^{|L|} (L_i - R_i)^2$$\n",
    "\n",
    "Here $L$ is the vector of predicted labels and $R$ is the vector of real labels.\n",
    "\n",
    "We use 5-fold [cross-validation](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CCrossValidation.html) to compute MSE and have a look at how MSE varies with regularisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def xval_results(taus):\n",
    "    errors = []\n",
    "    for t in taus:\n",
    "        rr.set_tau(t)\n",
    "        splitting_strategy = CrossValidationSplitting(labels_train, 5)\n",
    "        # evaluation method\n",
    "        evaluation_criterium = MeanSquaredError()\n",
    "        # cross-validation instance\n",
    "        cross_validation = CrossValidation(rr, feats_train, labels_train, splitting_strategy, evaluation_criterium, False)\n",
    "        cross_validation.set_num_runs(100)\n",
    "        cross_validation.set_conf_int_alpha(0.05)\n",
    "        result = cross_validation.evaluate()\n",
    "        result = CrossValidationResult.obtain_from_generic(result)\n",
    "        errors.append(result.mean)\n",
    "    return errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data with dimension: 10 and number of samples: 200 is now sampled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n = 500\n",
    "taus = logspace(-6, 4, n)\n",
    "figure(figsize=(20,6))\n",
    "suptitle('Effect of Regularisation for 10-dimensional data with 200 samples', fontsize=12)\n",
    "\n",
    "matrix, y = generate_data(200,10)\n",
    "feats_train = RealFeatures(matrix.T)\n",
    "labels_train = RegressionLabels(y[0])\n",
    "weights, rr = generate_weights(taus, feats_train, labels_train)\n",
    "errors = xval_results(taus)\n",
    "\n",
    "p1=subplot(121)\n",
    "plot_regularization(taus, weights)\n",
    "\n",
    "p2 = subplot(122)\n",
    "plot(taus, errors)\n",
    "p2.set_xscale('log')\n",
    "xlabel('Tau', fontsize=12)\n",
    "ylabel('Error', fontsize=12)\n",
    "jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen from the plot of errors, regularisation doesn't seem to affect the errors significantly. One interpretation could be that this is beacuse there is less overfitting as we have large number of samples.  For a small sample size as compared to the dimensionality, the test set performance may be poor even. The reason for this is that the regression function will fit the noise too much, while the interesting part of the signal is too small. We now generate 10 samples of 10-dimensions to test this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(20,6))\n",
    "suptitle('Effect of Regularisation for 10-dimensional data with 10 samples', fontsize=12)\n",
    "matrix, y = generate_data(10,10)\n",
    "feats_train = RealFeatures(matrix.T)\n",
    "labels_train = RegressionLabels(y[0])\n",
    "weights, rr = generate_weights(taus, feats_train, labels_train)\n",
    "errors = xval_results(taus)\n",
    "\n",
    "p1 = subplot(121)\n",
    "plot_regularization(taus, weights)\n",
    "\n",
    "\n",
    "p2 = subplot(122)\n",
    "plot(taus, errors)\n",
    "p2.set_xscale('log')\n",
    "xlabel('Tau', fontsize=12)\n",
    "ylabel('Error', fontsize=12)\n",
    "jet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first plot is the famous ridge trace that is the signature of this technique. The plot is really very straight forward to read. It presents the standardized regression coefficients (weights) on the vertical axis and various values of tau (Regularisation constant) along the horizontal axis. Since the values of tau ($\\tau$) span several orders of magnitude, we adopt a logarithmic scale along this axis. As tau is increased, the values of the regression estimates change, often wildly at first. At some point, the coefficients seem to settle down and then gradually drift towards zero. Often the value of tau for which these coefficients are at their stable values is the best one. This should be supported by a low error value for that tau.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least Angle Regression and LASSO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LASSO (Least Absolute Shrinkage and Selection Operator) is another version of Least Squares regression, which uses a L1-norm of the parameter vector. This intuitively enforces sparse solutions, whereas L2-norm penalties usually result in smooth and dense solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\min \\|X^T\\beta - y\\|^2 + \\lambda\\|\\beta\\|_1$$\n",
    "\n",
    "In Shogun, following equivalent form is solved, where increasing $C$ selects more variables:\n",
    "\n",
    "$$\\min \\|X^T\\beta - y\\|^2 \\quad s.t. \\|\\beta\\|_1 \\leq C $$\n",
    "\n",
    "One way to solve this regularized form is by using [Least Angle Regression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLeastAngleRegression.html) (LARS)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " LARS is essentially [forward stagewise](http://en.wikipedia.org/wiki/Stepwise_regression) made fast. LARS can be briefly described as follows.\n",
    "\n",
    "1. Start with an empty set.\n",
    "2. Select $x_j$ that is most correlated with residuals.\n",
    "3. Proceed in the direction of $x_j$ until another variable $x_k$ is equally correlated with residuals.\n",
    "4. Choose equiangular direction between $x_j$ and $x_k$.\n",
    "5. Proceed until third variable enters the active set, etc.\n",
    "\n",
    "It should be noticed that instead of making tiny hops in the direction of one variable at a time, LARS makes optimally-sized leaps in optimal directions. These directions are chosen to make equal angles (equal correlations) with each of the variables currently in our set (equiangular). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shogun provides tools for Least angle regression (LARS) and lasso using [CLeastAngleRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLeastAngleRegression.html) class. As explained in the mathematical formaulation, LARS is just like [Stepwise Regression](http://en.wikipedia.org/wiki/Stepwise_regression) but increases the estimated variables in a direction equiangular to each one's correlations with the residual. The working of this is shown below by plotting the LASSO path. Data is generated in a similar way to the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#sample some data\n",
    "X=rand(10)*1.5\n",
    "for i in range(9):\n",
    "    x=random.standard_normal(10)*0.5\n",
    "    X=vstack((X, x))\n",
    "y=ones(10)\n",
    "\n",
    "feats_train=RealFeatures(X)\n",
    "labels_train=RegressionLabels(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[CLeastAngleRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLeastAngleRegression.html) requires the features to be normalized with a zero mean and unit norm. Hence we use two preprocessors: [PruneVarSubMean](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CPruneVarSubMean.html) and [NormOne](http://www.shogun-toolbox.org/doc/en/latest/classshogun_1_1CPruneVarSubMean.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Preprocess data\n",
    "preproc=PruneVarSubMean()\n",
    "preproc.init(feats_train)\n",
    "feats_train.add_preprocessor(preproc)\n",
    "feats_train.apply_preprocessor()\n",
    "\n",
    "preprocessor=NormOne()\n",
    "preprocessor.init(feats_train)\n",
    "feats_train.add_preprocessor(preprocessor)\n",
    "feats_train.apply_preprocessor()\n",
    "\n",
    "print \"(No. of attributes, No. of samples) of data:\"\n",
    "print feats_train.get_feature_matrix().shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we train on the data. Keeping in mind that we had 10 attributes/dimensions in our data, let us have a look at the size of LASSO path which is obtained readily using `get_path_size()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Train and generate weights\n",
    "la=LeastAngleRegression()\n",
    "la.set_labels(labels_train)\n",
    "la.train(feats_train)\n",
    "\n",
    "size=la.get_path_size()\n",
    "print (\"Size of path is %s\" %size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weights generated ($\\beta_i$) and their norm ($\\sum_i|\\beta_i|$) change with each step. This is when a new variable is added to path. To get the weights at each of these steps `get_w_for_var()` method is used. The argument is the index of the variable which should be in the range [0, path_size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#calculate weights\n",
    "weights=[]\n",
    "for i in range(size):\n",
    "    weights.append(la.get_w_for_var(i))  \n",
    "s = sum(abs(array(weights)), axis=1)\n",
    "print ('Max. norm is %s' %s[-1])\n",
    "\n",
    "figure(figsize(30,7))\n",
    "#plot 1\n",
    "ax=subplot(131)\n",
    "title('Lasso path')\n",
    "ax.plot(s, weights, linewidth=2)\n",
    "ymin, ymax = ylim()\n",
    "ax.vlines(s[1:-1], ymin, ymax, linestyle='dashed')\n",
    "xlabel(\"Norm\")\n",
    "ylabel(\"weights\")\n",
    "\n",
    "#Restrict norm to half for early termination\n",
    "la.set_max_l1_norm(s[-1]*0.5)\n",
    "la.train(feats_train)\n",
    "size=la.get_path_size()\n",
    "weights=[]\n",
    "for i in range(size):\n",
    "    weights.append(la.get_w_for_var(i))\n",
    "s = sum(abs(array(weights)), axis=1)\n",
    "\n",
    "#plot 2\n",
    "ax2=subplot(132)\n",
    "title('Lasso path with restricted norm')\n",
    "ax2.plot(s, weights, linewidth=2)\n",
    "ax2.vlines(s[1:-1], ymin, ymax, linestyle='dashed')\n",
    "xlabel(\"Norm\")\n",
    "ylabel(\"weights\")\n",
    "print ('Restricted norm is  %s' %(s[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each color in the plot represents a coefficient and the vertical lines denote steps. It is clear that the weights are piecewise linear function of the norm.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kernel Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kernel ridge regression (KRR) is a kernel-based regularized form of regression. The dual form of Ridge regression can be shown to be:\n",
    "$${\\bf \\alpha}=\\left({\\bf X}^T{\\bf X}+\\tau{\\bf I}\\right)^{-1}{\\bf y} \\quad \\quad(1)$$ \n",
    "It can be seen that the equation to compute $\\alpha$ only contains the vectors $\\bf X$ in inner products with each other. If a non-linear mapping\n",
    "$\\Phi : x \\rightarrow \\Phi(x) \\in \\mathcal F$  is used, the equation can be defined in terms of inner products  $\\Phi(x)^T \\Phi(x)$ instead. We can then use the [kernel trick](http://en.wikipedia.org/wiki/Kernel_method) where a kernel function, which can be evaluated efficiently, is choosen $K({\\bf x_i, x_j})=\\Phi({\\bf x_i})\\Phi({\\bf x_j})$. This is done because it is sufficient to know these inner products only, instead of the actual vectors $\\bf x_i$. Linear regression methods like above discussed Ridge Regression can then be carried out in the feature space by using a kernel function representing a non-linear map which amounts to nonlinear regression in original input space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "KRR can be performed in Shogun using [CKernelRidgeRegression](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CKernelRidgeRegression.html) class. Let us apply it on a non linear regression problem from the [Boston Housing Dataset](https://archive.ics.uci.edu/ml/datasets/Housing), where the task is to predict prices of houses by finding a relationship with the various attributes provided. The per capita crime rate attribute is used in this particular example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "feats = RealFeatures(CSVFile('../../../data/uci/housing/fm_housing.dat'))\n",
    "train_labels = RegressionLabels(CSVFile('../../../data/uci/housing/housing_label.dat'))\n",
    "mat = feats.get_feature_matrix()\n",
    "\n",
    "crime_rate = mat[0]\n",
    "feats_train = RealFeatures(crime_rate.reshape(1, len(mat[0])))\n",
    "\n",
    "preproc = RescaleFeatures()\n",
    "preproc.init(feats_train)\n",
    "feats_train.add_preprocessor(preproc)\n",
    "feats_train.apply_preprocessor(True)\n",
    "\n",
    "# Store preprocessed feature matrix.\n",
    "preproc_data = feats_train.get_feature_matrix()\n",
    "\n",
    "size=500\n",
    "x1=linspace(0, 1, size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "width=0.5\n",
    "tau=0.5\n",
    "kernel=GaussianKernel(feats_train, feats_train, width)\n",
    "krr=KernelRidgeRegression(tau, kernel, train_labels)\n",
    "krr.train(feats_train)\n",
    "\n",
    "feats_test=RealFeatures(x1.reshape(1,len(x1)))\n",
    "kernel.init(feats_train, feats_test)\n",
    "out = krr.apply().get_labels()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Visualization of regression\n",
    "fig=figure(figsize(6,6))\n",
    "#first plot with only one attribute\n",
    "title(\"Regression with 1st attribute\")\n",
    "_=scatter(preproc_data[0:], train_labels.get_labels(), c=ones(506), cmap=gray(), s=20)\n",
    "_=xlabel('Crime rate')\n",
    "_=ylabel('Median value of homes')\n",
    "\n",
    "_=plot(x1,out, linewidth=3)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen from the example KRR (using the kernel trick) can apply techniques for linear regression in the feature space to perform nonlinear regression in the input space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Support Vector Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Kernel Ridge Regression $(1)$  we have seen the result to be a dense solution. Thus all training examples are active which limits its usage to fewer number of training examples. Support Vector Regression (SVR) uses the  concept of support vectors as in [Support Vector Machines]() that leads to a sparse solution. In the SVM the penalty was paid for being on the wrong side of the discriminating plane. Here we do the same thing: we introduce a penalty for being far away from predicted line, but once you are close enough, i.e. in some “epsilon-tube” around this line, there is no penalty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are given a labeled set of input-output pairs $\\mathcal{D}=(x_i,y_i)^N_{i=1}\\subseteq \\mathcal{X} \\times \\mathcal{Y}$ where $x\\in\\mathcal{X}$ and $y\\in \\mathcal{Y}$ and the primary problem is as follows:\n",
    "$$\\arg\\min_{\\mathbf{w},\\mathbf{\\xi}, b } ({\\frac{1}{2} \\|\\mathbf{w}\\|^2 +C \\sum_{i=1}^n (\\xi_i+ {\\xi_i}^*)) }$$\n",
    "\n",
    "For the constraints:\n",
    "$$ {\\bf w}^T{\\bf x}_i+b-c_i-\\xi_i\\leq 0, \\, \\forall i=1\\dots N$$\n",
    "$$ -{\\bf w}^T{\\bf x}_i-b-c_i^*-\\xi_i^*\\leq 0, \\, \\forall i=1\\dots N $$\n",
    "with $c_i=y_i+ \\epsilon$ and $c_i^*=-y_i+ \\epsilon$\n",
    "\n",
    "The resulting dual optimaization problem is:\n",
    "$$ \\max_{{\\bf \\alpha},{\\bf \\alpha}^*} -\\frac{1}{2}\\sum_{i,j=1}^N(\\alpha_i-\\alpha_i^*)(\\alpha_j-\\alpha_j^*) {\\bf x}_i^T {\\bf x}_j-\\sum_{i=1}^N(\\alpha_i+\\alpha_i^*)\\epsilon - \\sum_{i=1}^N(\\alpha_i-\\alpha_i^*)y_i\\\\$$ $$ \\mbox{wrt}:$$\n",
    "$${\\bf \\alpha},{\\bf \\alpha}^*\\in{\\bf R}^N\\\\ \\mbox{s.t.}: 0\\leq \\alpha_i,\\alpha_i^*\\leq C,\\, \\forall i=1\\dots N\\\\ \\sum_{i=1}^N(\\alpha_i-\\alpha_i^*)y_i=0 $$\n",
    "\n",
    "This class also support the $\\nu$-SVR regression version of the problem, where $\\nu$ replaces the $\\epsilon$ parameter and represents an upper bound on the action of margin errors and a lower bound on the fraction of support vectors. The resulting problem generally takes a bit longer to solve. The details and comparison of these two versioins can be found in [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us try regression using Shogun's [LibSVR](http://shogun-toolbox.org/doc/en/latest/classshogun_1_1CLibSVR.html). The dataset from last section is used. The `svr_param` argument is  the $\\epsilon$-tube for the $\\epsilon$ version and is the $\\nu$ parameter in other case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Use different kernels\n",
    "gaussian_kernel=GaussianKernel(feats_train, feats_train, 0.1)\n",
    "#Polynomial kernel of degree 2\n",
    "poly_kernel=PolyKernel(feats_train, feats_train, 2, True)\n",
    "linear_kernel=LinearKernel(feats_train, feats_train)\n",
    "\n",
    "kernels=[linear_kernel, poly_kernel, gaussian_kernel]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "svr_param=1\n",
    "svr_C=10\n",
    "svr=LibSVR(svr_C, svr_param, gaussian_kernel, train_labels, LIBSVR_EPSILON_SVR)\n",
    "\n",
    "#Visualization of regression\n",
    "x1=linspace(0, 1, size)\n",
    "feats_test_=RealFeatures(x1.reshape(1,len(x1)))\n",
    "\n",
    "def svr_regress(kernels):\n",
    "    fig=figure(figsize(8,8))\n",
    "    for i, kernel in enumerate(kernels):\n",
    "        svr.set_kernel(kernel)\n",
    "        svr.train()\n",
    "        out=svr.apply(feats_test_).get_labels()\n",
    "        #subplot(1,len(kernels), i)\n",
    "        #first plot with only one attribute\n",
    "        title(\"Support Vector Regression\")\n",
    "        _=scatter(preproc_data[0:], train_labels.get_labels(), c=ones(506), cmap=gray(), s=20)\n",
    "        _=xlabel('Crime rate')\n",
    "        _=ylabel('Median value of homes')\n",
    "        _=plot(x1,out, linewidth=3)\n",
    "        ylim([0, 40])\n",
    "    p1 = Rectangle((0, 0), 1, 1, fc=\"r\")\n",
    "    p2 = Rectangle((0, 0), 1, 1, fc=\"b\")\n",
    "    p3 = Rectangle((0, 0), 1, 1, fc=\"g\")\n",
    "    _=legend((p1, p2, p3), [\"Gaussian Kernel\", \"Linear Kernel\", \"Polynomial Kernel\"], loc=1)\n",
    "\n",
    "    \n",
    "svr_regress(kernels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us do comparison of time taken for the 2 different models simliar to that done in section 6 of [1]. The [Boston Housing Dataset](https://archive.ics.uci.edu/ml/datasets/Housing) is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "gaussian_kernel=GaussianKernel(feats, feats, 13)\n",
    "nus=[0.2, 0.4, 0.6, 0.8]\n",
    "epsilons=[0.16, 0.09, 0.046, 0.0188]\n",
    "svr_C=10\n",
    "\n",
    "def compare_svr(nus, epsilons):\n",
    "    time_eps=[]\n",
    "    time_nus=[]\n",
    "    for i in range(len(epsilons)):\n",
    "        svr_param=1\n",
    "        svr=LibSVR(svr_C, epsilons[i], gaussian_kernel, train_labels, LIBSVR_EPSILON_SVR)\n",
    "        t_start=time.clock()\n",
    "        svr.train()\n",
    "        time_test=(time.clock() - t_start)\n",
    "        time_eps.append(time_test)\n",
    "        \n",
    "    for i in range(len(nus)):\n",
    "        svr_param=1\n",
    "        svr=LibSVR(svr_C, nus[i], gaussian_kernel, train_labels,  LIBSVR_NU_SVR)\n",
    "        t_start=time.clock()\n",
    "        svr.train()\n",
    "        time_test=(time.clock() - t_start)\n",
    "        time_nus.append(time_test)\n",
    "\n",
    "    print \"-\"*72      \n",
    "    print \"|\", \"%15s\" % 'Nu' ,\"|\", \"%15s\" % 'Epsilon',\"|\",\"%15s\" % 'Time (Nu)' ,\"|\", \"%15s\" % 'Time(Epsilon)' ,\"|\"\n",
    "    for i in range(len(nus)):\n",
    "        print \"-\"*72      \n",
    "        print \"|\", \"%15s\" % nus[i] ,\"|\", \"%15s\" %epsilons[i],\"|\",\"%15s\" %time_nus[i] ,\"|\", \"%15s\" %time_eps[i] ,\"|\"  \n",
    "    print \"-\"*72  \n",
    " \n",
    "title_='SVR Performance on Boston Housing dataset'\n",
    "print \"%50s\" %title_\n",
    "compare_svr(nus, epsilons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Training $\\nu$-Support Vector Regression: Theory and Algorithms by Chih-Chung Chang and Chih-Jen Lin"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}