Lasso and Least angle regression

LassoRegression.ipynb

{
 "metadata": {
  "name": "LassoRegression"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Biblio\n",
      "\n",
      " - The [lasso/lars page](http://www-stat.stanford.edu/~tibs/lasso.html) of Tibshirani and the [original Lasso paper](http://www-stat.stanford.edu/~tibs/lasso/lasso.pdf);\n",
      " - The [Annals of statistics 2004](http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.aos/1083178935) paper (with blurred figures!)\n",
      " - An [earlier pre-print](http://www.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdf) with nice vectorial figures\n",
      " - [Slides](http://www.stanford.edu/~hastie/TALKS/bradfest.pdf) from hastie"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Example on real data sets\n",
      "\n",
      "Classic [Boston housing data set](http://archive.ics.uci.edu/ml/datasets/Housing)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!wget http://lib.stat.cmu.edu/datasets/boston_corrected.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "--2013-05-20 23:05:32--  http://lib.stat.cmu.edu/datasets/boston_corrected.txt\r\n",
        "Resolving lib.stat.cmu.edu (lib.stat.cmu.edu)... "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "128.2.234.14\r\n",
        "Connecting to lib.stat.cmu.edu (lib.stat.cmu.edu)|128.2.234.14|:80... "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "connected.\r\n",
        "HTTP request sent, awaiting response... "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "200 OK\r\n",
        "Length: 62136 (61K) [text/plain]\r\n",
        "Saving to: `boston_corrected.txt'\r\n",
        "\r\n",
        "\r",
        " 0% [                                       ] 0           --.-K/s              "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "13% [====>                                  ] 8,398       35.7K/s              "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "54% [====================>                  ] 33,900      76.9K/s              "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        "97% [====================================>  ] 60,380      92.2K/s              \r",
        "100%[======================================>] 62,136      93.4K/s   in 0.6s    \r\n",
        "\r\n",
        "2013-05-20 23:05:35 (93.4 KB/s) - `boston_corrected.txt' saved [62136/62136]\r\n",
        "\r\n"
       ]
      }
     ],
     "prompt_number": 100
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "! head -n 15 boston_corrected.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "This file contains the Harrison and Rubinfeld (1978) data corrected for a few minor errors and augmented with the latitude and longitude of the observations. This file appears under boston in the statlib index. One can obtain matlab and spreadsheet versions of the information below from www.finance.lsu/re under spatial statistics links.\r\n",
        "\r\n",
        "Harrison, David, and Daniel L. Rubinfeld, \ufffdHedonic Housing Prices and the Demand for Clean Air,\ufffd Journal of Environmental Economics and Management, Volume 5, (1978), 81-102. Original data.\r\n",
        "\r\n",
        "Gilley, O.W., and R. Kelley Pace, \ufffdOn the Harrison and Rubinfeld Data,\ufffd Journal of Environmental Economics and Management, 31 (1996), 403-405. Provided corrections and examined censoring.\r\n",
        "\r\n",
        "Pace, R. Kelley, and O.W. Gilley, \ufffdUsing the Spatial Configuration of the Data to Improve Estimation,\ufffd  Journal of the Real Estate Finance and Economics 14 (1997), 333-340. Added georeferencing and spatial estimation.\r\n",
        "\r\n",
        "\r\n",
        "OBS.\tTOWN\tTOWN#\tTRACT\tLON\tLAT\tMEDV\tCMEDV\tCRIM\tZN\tINDUS\tCHAS\tNOX\tRM\tAGE\tDIS\tRAD\tTAX\tPTRATIO\tB\tLSTAT\r\n",
        "1\tNahant\t0\t2011\t-70.955000\t42.255000\t24.0\t24.0\t0.00632\t18.0\t2.31\t0\t0.538\t6.575\t65.2\t4.0900\t1\t296\t15.3\t396.90\t4.98\r\n",
        "2\tSwampscott\t1\t2021\t-70.950000\t42.287500\t21.6\t21.6\t0.02731\t0.0\t7.07\t0\t0.469\t6.421\t78.9\t4.9671\t2\t242\t17.8\t396.90\t9.14\r\n",
        "3\tSwampscott\t1\t2022\t-70.936000\t42.283000\t34.7\t34.7\t0.02729\t0.0\t7.07\t0\t0.469\t7.185\t61.1\t4.9671\t2\t242\t17.8\t392.83\t4.03\r\n",
        "4\tMarblehead\t2\t2031\t-70.928000\t42.293000\t33.4\t33.4\t0.03237\t0.0\t2.18\t0\t0.458\t6.998\t45.8\t6.0622\t3\t222\t18.7\t394.63\t2.94\r\n",
        "5\tMarblehead\t2\t2032\t-70.922000\t42.298000\t36.2\t36.2\t0.06905\t0.0\t2.18\t0\t0.458\t7.147\t54.2\t6.0622\t3\t222\t18.7\t396.90\t5.33\r\n"
       ]
      }
     ],
     "prompt_number": 101
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import pandas as pd\n",
      "boston = pd.read_csv('boston_corrected.txt', sep='\\t',\n",
      "    skiprows=9, index_col='OBS.')\n",
      "boston.head()[['TOWN', 'LAT', 'LON', 'CRIM', 'AGE', 'B', 'TAX', 'INDUS', 'MEDV']]\n",
      "print boston.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(507, 20)\n"
       ]
      }
     ],
     "prompt_number": 123
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.subplot?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 187
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "f, (ax1, ax2) = plt.subplots(1,2)\n",
      "ax1.hist(boston['MEDV'], bins=20)\n",
      "\n",
      "ax2.hist(boston['MEDV'], bins=np.logspace(0, 3, 20))\n",
      "ax2.set_xscale('log')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEECAYAAAA8tB+vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVNX5B/BvIOm0NESSQjaYRZdC3gnZLRjsjJSNYWPf\niEEwEhR3SFItjh2pjAidURL6E5aqUwPW0So6a5wSM1pD6kAaKGykdjRYEmsnpmHaTfPC7lbYoAnB\nhiTn9wdlZZN9f8ndvfl+ZnZm976cPHf37pOz5557TowQQoCIiGRjhtQBEBFRaDGxExHJDBM7EZHM\nMLETEckMEzsRkcwwsRMRyYzXxF5TU4Pc3FwsWbIENTU1AAC73Q6dTof09HQUFRXh4sWLYQ+UKBC9\nvb0oKChATk4OlixZgv379wMAqqqqoFQqodFooNFocPToUcc+e/fuRVpaGjIzM9Hc3CxV6EQBi/HU\nj/3vf/87ysrKcPr0acTFxeH73/8+XnzxRbz00kuYO3cutm/fjn379mFgYAAGg2Eq4ybyidVqhdVq\nhVqtxtDQEJYtW4aGhgbU19dj9uzZePTRR5227+jowMaNG3H69Gn09/dj9erV6OrqwowZ/HFL0cPj\n2drZ2YkVK1bg61//OmbOnIlVq1bh7bffRmNjI/R6PQBAr9ejoaFhSoIl8ldKSgrUajUAID4+HllZ\nWejv7wcAuKrTHD58GGVlZYiLi4NKpcLixYvR2to6pTETBctjYl+yZAlOnToFu92O4eFhHDlyBH19\nfbDZbFAoFAAAhUIBm802JcESBaO7uxttbW249dZbAQAHDhxAXl4eKioqHM2J586dg1KpdOyjVCod\n/wiIokWsp5WZmZl4/PHHUVRUhG9+85tQq9WYOXOm0zYxMTGIiYlxub+75USh4uuIGENDQ1i/fj1q\namoQHx+PLVu24MknnwQAPPHEE9i2bRsOHjzoct+J5zHPa5oKwYz24rXhsLy8HB999BFaWlqQmJiI\n9PR0KBQKWK1WAIDFYkFycrLH4IJ57Nq1S/IyIiEGljH54asrV65g3bp1uO+++1BSUgIASE5OdlRK\nKisrHc0tqamp6O3tdezb19eH1NTUkJ7Xvhy7u21cLZ+47PrX3p4H8zkEcxxyOhZ/j8OX+IPlNbH/\n5z//AQD09PTg97//PTZu3Iji4mIYjUYAgNFodHxZiCKNEAIVFRXIzs7G1q1bHcstFovj+TvvvIPc\n3FwAQHFxMerq6jAyMgKz2YyzZ88iPz8/pDFptdqAt3G1fOKy61/78jxQwRyHu3XReCz+HsfE19ee\nh+I4HIQXK1euFNnZ2SIvL0+cOHFCCCHEhQsXRGFhoUhLSxM6nU4MDAy43NeH4r3atWuX5GVEQgws\nYzJfzq9Tp06JmJgYkZeXJ9RqtVCr1eLIkSNi06ZNIjc3VyxdulTceeedwmq1OvZ56qmnxKJFi0RG\nRoZoamoK6O9Gi1B8DpFCTscS7DkW1jM0FF+AkydPSl5GJMTAMiaTKsHKKbGH4nOIFHI6lmDPMY/9\n2IMVExMTkvYiIlekOr94XlO4BXuO8a4LIiKZYWInIpIZJnYiIplhYicikhkmdiIimWFiJyKSGSZ2\nIiKZYWInIpIZJnYiIplhYicikhkmdiIimWFinyAhIckxTrerR0JCktQhEkmK35HIx0HAJrg6O46n\nmKPvmOSKg4BJg9+R8OMgYERE5ISJnYhIZrwm9r179yInJwe5ubnYuHEj/vvf/8Jut0On0yE9PR1F\nRUWOGd6JiEh6HhN7d3c3Xn75ZZw5cwaffPIJxsbGUFdXB4PBAJ1Oh66uLhQWFsJgMExVvERE5IXH\nxJ6QkIC4uDgMDw9jdHQUw8PDuPHGG9HY2Ai9Xg8A0Ov1aGhomJJgiYjIO4+JPSkpCdu2bcNNN92E\nG2+8EXPmzIFOp4PNZoNCoQAAKBQK2Gy2KQmWiIi8i/W08p///Ceee+45dHd344YbbsDdd9+NN954\nw2mba31X3amqqnI812q10Gq1QQVM05fJZILJZJI6DKKI57Ef+5tvvoljx47hlVdeAQDU1tbigw8+\nwIkTJ3Dy5EmkpKTAYrGgoKAAnZ2dkwuPwv6+7KMbPdiPXRr8joRfWPuxZ2Zm4oMPPsDly5chhMDx\n48eRnZ2NNWvWwGg0AgCMRiNKSkoCDoCIiELL652nv/rVr2A0GjFjxgx85zvfwSuvvILBwUGUlpai\np6cHKpUK9fX1mDNnzuTCo7Bmw9pI9GCNXRr8joRfsOcYhxSYgCdt9GBilwa/I+HHIQWIiMgJEzsR\nkcwwsRMRyQwTOxGRzDCxExHJDBM7EZHMMLETEckMEzsRkcwwsRMRyQwTOxGRzDCxExHJDBM7yVpv\nby8KCgqQk5ODJUuWYP/+/QDgcd7evXv3Ii0tDZmZmWhubpYqdKKAcRCwCTjAUfTw5fyyWq2wWq1Q\nq9UYGhrCsmXL0NDQgNdeew1z587F9u3bsW/fPgwMDMBgMKCjowMbN27E6dOn0d/fj9WrV6Orqwsz\nZnxVB4rG8zqU+B0JPw4CRuRBSkoK1Go1ACA+Ph5ZWVno7+93O2/v4cOHUVZWhri4OKhUKixevBit\nra2SxU8UCI9T45H/EhKSMDg44HLd7NmJ+OIL+xRHRNd0d3ejra0NK1ascDtv77lz53Drrbc69lEq\nlejv759UFqd8pFAK9bSPTOwhdjWpu/4JNTjofm5YCq+hoSGsW7cONTU1mD17ttM6b/P2ulp3fWIn\nCtbEykF1dXVQ5XltivnHP/4BjUbjeNxwww3Yv3+/x4tPRJHkypUrWLduHTZt2uSYxlGhUMBqtQIA\nLBYLkpOTAQCpqano7e117NvX14fU1NSpD5ooCF4Te0ZGBtra2tDW1oa//vWvmDVrFtauXQuDwQCd\nToeuri4UFhbCYDBMRbwRINZRw3P1oMgihEBFRQWys7OxdetWx/Li4mKX8/YWFxejrq4OIyMjMJvN\nOHv2LPLz8yWJnShgwg9//OMfxW233SaEECIjI0NYrVYhhBAWi0VkZGRM2t7P4iMCAAEID49g1kff\n+xHJfHk/T506JWJiYkReXp5Qq9VCrVaLo0ePigsXLojCwkKRlpYmdDqdGBgYcOzz1FNPiUWLFomM\njAzR1NQU0N+VM1++AxScYN9Dv7o7lpeXY/ny5XjooYeQmJiIgYGBa/8ckJSU5Hh9TUxMDHbt2uV4\nHQ0XmXzpyhX4enYDC8bEC0zV1dWc81QC7O4YflM2mfXIyAhSU1PR0dGBefPmOSV2AEhKSoLd7tzj\nIxq/AEzs0YOTWUuDiT38pqwf+9GjR7Fs2TLMmzcPgPuLT0REJC2fE/uhQ4dQVlbmeO3u4hMREUnL\np6aYS5cu4eabb4bZbHb0Abbb7SgtLUVPTw9UKhXq6+sxZ84c58Kj8Ccrm2KiB5tipMGmmPCbsjb2\ngAqPwi8AE3v0YGKXBhN7+HGsGCIicsLETkQkM0zsREQyw8RORCQzTOxERDLDxE5EJDNM7EREMsPE\nTkQkM0zsREQyw8RORCQzTOxERDLDxE5EJDNM7EREMsPETkQkM0zsREQyw8RORCQzPiX2ixcvYv36\n9cjKykJ2djY+/PBD2O126HQ6pKeno6ioCBcvXgx3rERE5AOfEvsjjzyCH/7wh/j000/xt7/9DZmZ\nmTAYDNDpdOjq6kJhYSEMBkO4YyUiIh94nRrv888/h0ajwb/+9S+n5ZmZmWhpaYFCoYDVaoVWq0Vn\nZ6dz4VE4hRinxosenBpPGpwaL/yCPcdivW1gNpsxb948bN68GR9//DGWLVuG5557DjabDQqFAgCg\nUChgs9lc7l9VVeV4rtVqodVqAw6WpjeTyQSTySR1GEQRz2uN/aOPPsJ3v/td/OUvf8Ett9yCrVu3\nYvbs2Xj++ecxMDDg2C4pKQl2u9258Cis2bDGHj1YY5cGa+zhF/bJrJVKJZRKJW655RYAwPr163Hm\nzBmkpKTAarUCACwWC5KTkwMOgoiIQsdrYk9JScGCBQvQ1dUFADh+/DhycnKwZs0aGI1GAIDRaERJ\nSUl4IyUiIp94bYoBgI8//hiVlZUYGRnBokWL8Nprr2FsbAylpaXo6emBSqVCfX095syZ41x4FP5k\nZVNM9GBTjDTYFBN+wZ5jPiX2gAuPwi8AE3v0YGKXBhN7+IW9jZ2IiKILEzsRkcwwsRMRyQwTOxGR\nzDCxExHJDBM7yVp5eTkUCgVyc3Mdy6qqqqBUKqHRaKDRaHD06FHHur179yItLQ2ZmZlobm6WImSi\noLG74wTs7hg9fDm/Tp06hfj4eNx///345JNPAADV1dWYPXs2Hn30UadtOzo6sHHjRpw+fRr9/f1Y\nvXo1urq6MGOGc/0nGs/rUGJ3x/Bjd0ciD1auXInExMRJy119aQ4fPoyysjLExcVBpVJh8eLFaG1t\nnYowiUKKiZ2mpQMHDiAvLw8VFRWOSWLOnTsHpVLp2EapVKK/v1+qEIkC5nXYXiK52bJlC5588kkA\nwBNPPIFt27bh4MGDLre92uwwGYejplAK9ZDUTOw07Vw/EmllZSXWrFkDAEhNTUVvb69jXV9fH1JT\nU12WcX1iJwrWxMpBdXV1UOWxKYamHYvF4nj+zjvvOHrMFBcXo66uDiMjIzCbzTh79izy8/OlCpMo\nYKyxk6yVlZWhpaUF58+fx4IFC1BdXQ2TyYT29nbExMRg4cKFeOmllwAA2dnZKC0tRXZ2NmJjY/HC\nCy+4bYohimTs7jgBuztGD47uKA12dww/dnckIiInPjXFqFQqJCQkYObMmYiLi0Nrayvsdjvuuece\n/Pvf/3Y70QYREU09n2rsMTExMJlMaGtrc9ywYTAYoNPp0NXVhcLCQhgMhrAGSkREvvG5KWZie09j\nYyP0ej0AQK/Xo6GhIbSRERFRQHyusa9evRrLly/Hyy+/DACw2WxQKBQAAIVCAZvNFr4oiYjIZz61\nsb///vuYP38+PvvsM+h0OmRmZjqtj4mJ4R16FHahvjuPSK787u5YXV2N+Ph4vPzyyzCZTEhJSYHF\nYkFBQQE6OzudC4/CbmHs7hg92N1RGuzuGH5h7+44PDyMwcFBAMClS5fQ3NyM3NxcFBcXw2g0AgCM\nRiNKSkoCDoKIiELHa43dbDZj7dq1AIDR0VHce++92LlzJ+x2O0pLS9HT0+O2u2M01mxYY48erLFL\ngzX28Av2HOOdpxMwsUcPJnZpMLGHH+88JSIiJ0zsREQyw8RORCQzTOxERDLDxE5EJDNM7EREMsPE\nTkQkM0zsREQyw8RORCQzTOxERDLDxE5EJDNM7EREMsPETkQkM0zsREQyw8RORCQz0zKxJyQkOeZp\nnfggIop2PiX2sbExaDQarFmzBgBgt9uh0+mQnp6OoqIiXLx4MaxBhtrg4ACuThTg6kFEFN18Suw1\nNTXIzs521GgNBgN0Oh26urpQWFgIg8EQ1iCJiMh3XhN7X18fjhw5gsrKSsdUTY2NjdDr9QAAvV6P\nhoaG8EZJREQ+85rYf/7zn+Ppp5/GjBlfbWqz2aBQKAAACoUCNpstfBESEZFfYj2tfPfdd5GcnAyN\nRgOTyeRyG28XHauqqhzPtVottFptIHH6JSEh6X/t6CQnJpPJ7XlIRF+JER6mwv7FL36B2tpaxMbG\n4ssvv8QXX3yBu+66C6dPn4bJZEJKSgosFgsKCgrQ2dk5uXAJZ5H3Nou6+/XB7Ou9bM7eHjpSnl/T\n+XP05fs1nd+fUAj2HPPYFLNnzx709vbCbDajrq4Ot99+O2pra1FcXAyj0QgAMBqNKCkpCTgAIiIK\nLb/6sV9rctmxYweOHTuG9PR0nDhxAjt27AhLcERE5D+PTTFBF86mmEnr+BM1dNgUIw02xYRfWJti\niIgo+jCxk6yVl5dDoVAgNzfXsczTndN79+5FWloaMjMz0dzcLEXIkvM05AaH3YgOUZvYOd4L+WLz\n5s1oampyWubuzumOjg68+eab6OjoQFNTEx566CGMj49LEbakPA+5wSaWaBC1iZ3jvZAvVq5cicTE\nRKdl7u6cPnz4MMrKyhAXFweVSoXFixejtbV1ymMmClbUJnaiQLm7c/rcuXNQKpWO7ZRKJfr7+yWJ\nkSgYHu88JZI7b8137tZJcUc1yVeo76pmYqdpR6FQwGq1Ou6cTk5OBgCkpqait7fXsV1fXx9SU1Nd\nlnF9YicK1sTKQXV1dVDlsSmGph13d04XFxejrq4OIyMjMJvNOHv2LPLz86UMlSggrLGTrJWVlaGl\npQXnz5/HggULsHv3buzYsQOlpaU4ePAgVCoV6uvrAQDZ2dkoLS1FdnY2YmNj8cILL7CXFUWlqL3z\n1PPdb7zzdDrgnafhEdyd21fXy/n9mQq885SIiJwwsRMRyQwTOxGRzDCxExHJDBM7EZHMMLETEcmM\nx8T+5ZdfYsWKFVCr1cjOzsbOnTsBeB72lIimu1iPw/4mJCRJHaDsee3HPjw8jFmzZmF0dBS33XYb\nnnnmGTQ2NmLu3LnYvn079u3bh4GBAcfQp06Fsx/7pHXs3xs67MceHqHox85+7sEJez/2WbNmAQBG\nRkYwNjaGxMREt8OeEhGR9Lwm9vHxcajVaigUChQUFCAnJ8ftsKdERCQ9r2PFzJgxA+3t7fj8889x\nxx134OTJk07rvQ17yuFNKVRCPbQpkVz5NVbML3/5S3zjG9/AK6+8ApPJ5Bj2tKCgAJ2dnZMLZxv7\npHVsWwwdtrGHB9vYpRfWNvbz5887erxcvnwZx44dg0ajcTvsKRERSc9jU4zFYoFer8f4+DjGx8ex\nadMmFBYWQqPRuBz2lIiIpMdhe0O6r/ey+RM0dNgUEx5sipEeh+0lIiInTOxERDLDxE5EJDNM7ERE\nMsPETkQkM0zsREQyw8RORCQzTOxERDLDxE5EJDNM7EREMsPETkQkM0zsREQyw8QeQRISkjgBMBEF\njaM7hnRf72V7ej+8HRNHxHPG0R3Dg6M7So+jOxIRkRMmdiIimfGa2Ht7e1FQUICcnBwsWbIE+/fv\nBwDY7XbodDqkp6ejqKjIMYUeERFJy2sbu9VqhdVqhVqtxtDQEJYtW4aGhga89tprmDt3LrZv3459\n+/ZhYGAABoPBuXC2sU9axzb20GEbe3iwjV16YW9jT0lJgVqtBgDEx8cjKysL/f39aGxshF6vBwDo\n9Xo0NDQEHAQREYWOx8msJ+ru7kZbWxtWrFgBm80GhUIBAFAoFLDZbC73qaqqcjzXarXQarUBB0vT\nm8lkgslkkjoMoojnc3fHoaEhrFq1Ck888QRKSkqQmJiIgYEBx/qkpCTY7XbnwtkUM2kdm2JCh00x\n4cGmGOlNSXfHK1euYN26ddi0aRNKSkoAXK2lW61WAIDFYkFycnLAQRBJQaVSYenSpdBoNMjPzwfA\nTgEkD14TuxACFRUVyM7OxtatWx3Li4uLYTQaAQBGo9GR8IlCxdOduFdrlcGJiYmByWRCW1sbWltb\nAQAGgwE6nQ5dXV0oLCyc1CGAKBp4bYr585//jO9973tYunSp48u0d+9e5Ofno7S0FD09PVCpVKiv\nr8ecOXOcC2dTzKR1bIrxnS9NAsG8JwsXLsRHH32Eb33rW45lmZmZaGlpcfwi1Wq16OzsnBSXnD8L\nNsVIL9hzLGKHFEhISMLg4ICXraZTYo8DMOp239mzE/HFF3a366NRuBP7t7/9bdxwww2YOXMmHnzw\nQfzkJz9xunYkhEBSUpLTtaRrce3atcvxWm6dApjYp97EjgHV1dXyTOzBnVxyTOzT78sS7sRusVgw\nf/58fPbZZ9DpdDhw4ACKi4sl7RQQCZjYpcexYogCNH/+fADAvHnzsHbtWrS2trJTAMkCEztNS8PD\nwxgcHAQAXLp0Cc3NzcjNzWWnAJIFNsWEdF/vZbMpxnfhbIoxm81Yu3YtAGB0dBT33nsvdu7cCbvd\nLmmngEjAphjpyfbiKRO7v39XfhdXw93GHigmdib2cAv2HPNrSAGKZKPw9GUaHAy+3zcRRQe2sRMR\nyQwTOxGRzDCxExHJDBP7tBHrdsyVhIQkqYMjohDixdNpw/3FVV5YJZIX1tiJiGSGiZ2ISGaY2ImI\nZEayxB7uSRQik/sLmPI9ZqKJPH8PeDE/eJINKRDe25Yjd0iBSI0rEm/x5pAC0piKIQUi8XONJGEf\ntre8vBwKhQK5ubmOZZwXkogocnlN7Js3b0ZTU5PTMs4LSRS9pmcz6PTiU1NMd3c31qxZg08++QSA\nb/NCAmyKiZx9vZcdiT992RQTHpHQ1BKJn2skkWR0R5vNBoVCAQBQKBSw2Wxut62qqnI8l9vckDTV\nTP97EJEnQd956u3n2/WJnSg42v89rqmWJgyiCBdQYr/WBJOSkuJ1Xsh58xYGHBwREfkvoMR+bV7I\nxx9/3Ou8kOfPn3CxtBfAqkD+NBEReeH14mlZWRlaWlpw/vx5KBQK7N69G3feeafXeSEBTxdpzAC+\n7WadY+8g1kfuRcpIjSsSL1bx4ml48OJp5Iv4OU+Z2CNhX+9lR+IXiYk9PJjYI1/Yb1AiIqLowsRO\nRCQzTOwEDspEJC9M7ISvZldy/RgcHJAwNvIXhwwgTo1HJDNX/xF7u3hJcsYaOxGRzDCxExHJDBM7\nBc1zm+7XeGGWaIqxjZ2C5rlN1/PNKIODbO8lCjXW2ImIZIaJnYhIZpjYiYhkhomdiEhmmNiJiGSG\nvWJIYrG8zZ0oxIKqsTc1NSEzMxNpaWnYt29fqGKawBQBZURCDJFThskUfBlfxeFpnBppTM15HSlM\nUgcQMqE5L+Uh4MQ+NjaGhx9+GE1NTejo6MChQ4fw6aefhjK2/zFFQBmREIOUZTiP/lhQUBCCQaUC\niSP8pu68DlxoB/kyhSNESTCxfyXgxN7a2orFixdDpVIhLi4OGzZswOHDh0MZG0WMibXqXYiEmnU4\nTMV57UsCcreNyWS67oawa4+TcP95mHx4HihfyvC0jat1zsuufx98eR6oYD8Tb8s8vb72PJT/mAJO\n7P39/ViwYIHjtVKpRH9//6Tt4uK2TXrExj4V6J8lCitfz+tt27a5fezevdvj3wh1EvE9ebp7Hihf\nyvC0jat1zsuY2AMkAvTWW2+JyspKx+va2lrx8MMPO20DT4N888FHCB6hxvOaj0h5BCPgXjGpqano\n7e11vO7t7YVSqXTaRkzzCWkp+vC8JjkIuClm+fLlOHv2LLq7uzEyMoI333wTxcXFoYyNaMrxvCY5\nCLjGHhsbi+effx533HEHxsbGUFFRgaysrFDGRjTleF6TLATVkOPGzTffLHJzc4VarRa33HKLT/ts\n3rxZJCcniyVLljiWXbhwQaxevVqkpaUJnU4nBgYG/C5j165dIjU1VajVaqFWq8XRo0c9ltHT0yO0\nWq3Izs4WOTk5oqamxu9Y3JXhTyyXL18W+fn5Ii8vT2RlZYkdO3b4HYe7Mvx9T0ZHR4VarRY//vGP\n/Y7BXRn+xuDqnAokDqLpICyJXaVSiQsXLvi1z3vvvSfOnDnjlJQfe+wxsW/fPiGEEAaDQTz++ON+\nl1FVVSWeffZZn+OwWCyira1NCCHE4OCgSE9PFx0dHX7F4q4Mf2O5dOmSEEKIK1euiBUrVohTp075\n/Z64KsPfOJ599lmxceNGsWbNGiGE/5+LqzL8jcHVORVIHETTQdjGihF+XmBauXIlEhMTnZY1NjZC\nr9cDAPR6PRoaGvwuw99YUlJSoFarAQDx8fHIyspCf3+/X7G4K8PfWGbNmgUAGBkZwdjYGBITE/1+\nT1yV4U8cfX19OHLkCCorKx37+BuDqzLE1UqFTzFcM3F7f+MIl0uXLkGv1+OBBx7A7373O0liCBWz\n2YzKykrcfffdUocStMOHD+OBBx7Ahg0bcOzYManDCVhnZye2bNmC0tJSHDx40LedwvHfYuHChUKt\nVotly5aJ3/72tz7vZzabnWrbc+bMcTwfHx93eu1rGVVVVeLmm28WS5cuFeXl5X79XDebzeKmm24S\nX3zxRUCxXF/G4OCg37GMjY2JvLw8ER8fLx577DEhhP/viasy/Ilj/fr14syZM8JkMjmaUfyNwVUZ\n/r4Xrs6pQD+TUHv99dfFu+++K4QQ4p577pEkhlBbv3691CGEzMDAgKioqJA6jKCNjY2Ju+++26dt\nw1Jjf//999HW1oajR4/iN7/5DU6dOhV0mYHevr5lyxaYzWa0t7dj/vz52LZtm0/7DQ0NYd26daip\nqcHs2bMDimVoaAjr169HTU0N4uPj/Y5lxowZaG9vR19fH9577z2cPHnS7zgmlmEymXyO491330Vy\ncjI0Go3b2rW3GNyV4e974e2cCm54g8nKy8uhUCiQm5vrtNzVODLX39Q0c+bMkMUQKv4cS6QL5Fj+\n7//+Dw8//PBUhumVv8fxhz/8AT/60Y+wYcMG3/5AWP/FiKs1s2eeecanbSfWtjMyMoTFYhFCCHHu\n3DmRkZHhdxm+rrveyMiIKCoqEr/+9a8DjsVVGYHEcs3u3bvF008/HdB7MrEMX+PYuXOnUCqVQqVS\niZSUFDFr1ixx3333+RWDqzI2bdrkcwyuXDungnkvvHF1vWZ0dFQsWrRImM1mMTIyIvLy8kRHR4eo\nra111Ng3bNgQshhCxZ9juSZSa+z+HMv4+LjYvn27OH78uIQRuxbIZyKEEMXFxT6VH/Ia+/DwMAYH\nBwFcbXtsbm6e9F/JV8XFxTAajQAAo9GIkpISv8uwWCyO5++8847XWIQQqKioQHZ2NrZu3RpQLO7K\n8CeW8+fP4+LFiwCAy5cv49ixY9BoNH7F4a4Mq9XqUxx79uxBb28vzGYz6urqcPvtt6O2ttavGFyV\n8frrr/v1Xrg7p0Jxfrjj6nqNu3Fk7rrrLrz99tt46KGHIrLPuz/HYrfb8dOf/hTt7e0RWYv351ie\nf/55/OlPf8Jbb72Fl156SaKIXfPnOFpaWvDII4/gwQcfREFBgU/lh3w8dpvNhrVr1wIARkdHce+9\n96KoqMjrfmVlZWhpacH58+exYMEC7N69Gzt27HBcMFCpVKivr/erjOrqaphMJrS3tyMmJgYLFy70\n+gG///5Mnp6XAAABoElEQVT7eOONN7B06VJoNBoAwN69e/2KxVUZe/bswaFDh3yOxWKxQK/XY3x8\nHOPj49i0aRMKCwuh0Wh8jsNdGffff79f78k115o6/P1crhFCOMrYvn07Pv74Y59icHdOLV++PKA4\nAuVqHJkPP/wQs2bNwquvvhrWvx1q7o4lKSkJL774ooSR+c/dsRw4cAA/+9nPJIzMP+6OY9WqVVi1\napVfZYU8sS9cuBDt7e1+73fo0CGXy48fPx5UGeXl5X7Fcdttt2F8fDyoWNyV8YMf/MDnOHJzc3Hm\nzJlJy5OSknyOw10Zr7/+us9xXHP9yeVPDNfTarXQarUAgNraWp/3c3dOBRpHoOQ0IQiPJfKE8jg4\nNR6Rj3wZRyZa8FgiTyiPg4mdyEdyGkeGxxJ5QnocobvOSyQfGzZsEPPnzxdf+9rXhFKpFK+++qoQ\nQogjR46I9PR0sWjRIrFnzx6Jo/QNjyXyhPs4YoTgGKRERHLCphgiIplhYicikhkmdiIimWFiJyKS\nGSZ2IiKZYWInIpIZJnYiIplhYicikhkmdiIimfl/JOXHN3wlHlwAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 192
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[`sklearn.linear_model.Lasso`](http://scikit-learn.org/dev/modules/generated/sklearn.linear_model.Lasso.html) and  [`sklearn.grid_search.GridSearchCV`](http://scikit-learn.org/dev/modules/generated/sklearn.grid_search.GridSearchCV.html)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "boston = boston.dropna(axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 129
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn import cross_validation\n",
      "\n",
      "explanatory_cols = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT']\n",
      "X_train, X_test, y_train, y_test = cross_validation.train_test_split(\n",
      "    boston[explanatory_cols], boston.ix[:,'MEDV'], test_size=0.4, random_state=0)\n",
      "from sklearn.preprocessing import StandardScaler\n",
      "scaler = StandardScaler()\n",
      "X_train = scaler.fit_transform(X_train)\n",
      "X_test = scaler.transform(X_test)\n",
      "\n",
      "print boston.shape, X_train.shape, X_test.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(506, 20) (303, 13) (203, 13)\n"
       ]
      }
     ],
     "prompt_number": 253
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import itertools\n",
      "second_order_combinations = list(itertools.combinations_with_replacement(range(X_train.shape[1]), 2))\n",
      "\n",
      "X_train = np.transpose([X_train[:,i] * X_train[:,j]\n",
      "    for i,j in second_order_combinations])\n",
      "X_test = np.transpose([X_test[:,i] * X_test[:,j]\n",
      "    for i,j in second_order_combinations])\n",
      "print X_train.shape, y_train.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(303, 91) (303,)\n"
       ]
      }
     ],
     "prompt_number": 251
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.linear_model import LinearRegression \n",
      "ols = LinearRegression()\n",
      "ols_fit = clf.fit(X=X_train, y=y_train)\n",
      "\n",
      "print '{} out of {} non zero coefficients with ||w||_2 = {}'.format(\n",
      "    np.sum(np.abs(ols_fit.coef_)>=.01),\n",
      "    len(ols_fit.coef_),\n",
      "    np.linalg.norm(ols_fit.coef_, ord=2))\n",
      "\n",
      "print ols_fit.score(X=X_train, y=y_train), ols_fit.score(X=X_test, y=y_test) \n",
      "\n",
      "def get_significant_coef(coefs, explanatory_cols, second_order_combinations, threshold=.5):\n",
      "    return ['{} * {} = {}'.format(explanatory_cols[i], explanatory_cols[j], c) for (i,j), c in zip(second_order_combinations, coefs) if abs(c) > threshold]\n",
      "print get_significant_coef(ols_fit.coef_, explanatory_cols, second_order_combinations)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3 out of 13 non zero coefficients with ||w||_2 = 4.83142803599\n",
        "0.689113551553 0.625487738575\n",
        "['CRIM * RM = 2.45410544226', 'CRIM * PTRATIO = -1.8409232126', 'CRIM * LSTAT = -3.73243420177']\n"
       ]
      }
     ],
     "prompt_number": 254
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for p in np.linspace(.1, .9, 10):\n",
      "    x_t, a, y_t, b = cross_validation.train_test_split(\n",
      "        boston[explanatory_cols], boston.ix[:,'MEDV'], test_size=p, random_state=0) \n",
      "    ols_fit = LinearRegression().fit(X=x_t, y=y_t)\n",
      "    print 'x_t.shape = {},  {} out of {} non zero coefficients with ||w||_2 = {}'.format(\n",
      "        x_t.shape,\n",
      "        np.sum(np.abs(ols_fit.coef_)>=.01),\n",
      "        len(ols_fit.coef_),\n",
      "        np.linalg.norm(ols_fit.coef_, ord=2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x_t.shape = (455, 13),  11 out of 13 non zero coefficients with ||w||_2 = 16.7136620181\n",
        "x_t.shape = (410, 13),  10 out of 13 non zero coefficients with ||w||_2 = 17.3710283206\n",
        "x_t.shape = (365, 13),  10 out of 13 non zero coefficients with ||w||_2 = 16.4377194102\n",
        "x_t.shape = (320, 13),  11 out of 13 non zero coefficients with ||w||_2 = 18.0964699833\n",
        "x_t.shape = (275, 13),  11 out of 13 non zero coefficients with ||w||_2 = 20.3958815354\n",
        "x_t.shape = (230, 13),  11 out of 13 non zero coefficients with ||w||_2 = 17.7868462582\n",
        "x_t.shape = (185, 13),  12 out of 13 non zero coefficients with ||w||_2 = 15.3978473822\n",
        "x_t.shape = (140, 13),  10 out of 13 non zero coefficients with ||w||_2 = 16.6725552627\n",
        "x_t.shape = (95, 13),  11 out of 13 non zero coefficients with ||w||_2 = 12.4778921797\n",
        "x_t.shape = (50, 13),  11 out of 13 non zero coefficients with ||w||_2 = 9.8257665043\n"
       ]
      }
     ],
     "prompt_number": 255
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.linear_model import Lasso, Ridge\n",
      "ridge_regressor = Ridge(alpha=1, max_iter=10000)\n",
      "ridge_fit = ridge_regressor.fit(X=X_train, y=y_train)\n",
      "\n",
      "from sklearn.grid_search import GridSearchCV\n",
      "ridge_cv = GridSearchCV(Ridge(), param_grid={'alpha':np.logspace(0,3,100)}, cv=3)\n",
      "ridge_cv.fit(X=X_train, y=y_train)\n",
      "\n",
      "ridge_cv.best_params_\n",
      "\n",
      "plt.plot(ridge_cv.param_grid['alpha'], [score for params, score, fold_scores in ridge_cv.grid_scores_]);\n",
      "\n",
      "print '{} out of {} non zero coefficients with ||w||_2 = {}'.format(\n",
      "    np.sum(ridge_cv.best_estimator_.coef_>=.01),\n",
      "    len(ridge_cv.best_estimator_.coef_),\n",
      "    np.linalg.norm(ridge_cv.best_estimator_.coef_, ord=2))\n",
      "\n",
      "plt.xscale('log')\n",
      "print ridge_cv.best_estimator_.score(X=X_train, y=y_train), ridge_cv.best_estimator_.score(X=X_test, y=y_test)\n",
      "\n",
      "print get_significant_coef(ridge_cv.best_estimator_.coef_, explanatory_cols, second_order_combinations)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "6 out of 13 non zero coefficients with ||w||_2 = 7.00985085438\n",
        "0.766331309749 0.688536093212\n",
        "['CRIM * CRIM = -0.856731829886', 'CRIM * ZN = 1.26193356082', 'CRIM * CHAS = 0.657234629535', 'CRIM * NOX = -2.02862669449', 'CRIM * RM = 2.38964596169', 'CRIM * DIS = -3.09571278554', 'CRIM * RAD = 1.86044309443', 'CRIM * TAX = -1.59886442926', 'CRIM * PTRATIO = -2.38502341411', 'CRIM * LSTAT = -3.86919632853']\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEECAYAAADAoTRlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X90VOWdx/H3wEQhiLCRqGEmGiQxGSShwQSlgAQLBKsN\n8qs7SCtigBxOUTlr1x+traR1tXHxx0pWzLGltSghFdtGFAKGOrrgmlFQqwU1CrHTbEUQwUjUkGH2\njyv5AclkJpnJncl8Xuc8Z37c57nzjVe+97nPvfe5Fp/P50NERGJKP7MDEBGR3qfkLyISg5T8RURi\nkJK/iEgMUvIXEYlBSv4iIjGoy+RfVVVFRkYGaWlplJSUnLZ81apVZGdnk52dTWZmJlarlSNHjgCQ\nkpJCVlYW2dnZjBs3LvTRi4hIt1j8Xefv9XpJT0+nuroam81Gbm4u5eXlOByODus/99xzPPzww1RX\nVwMwYsQIdu3aRUJCQniiFxGRbvHb83e73aSmppKSkkJcXBxOp5PKyspO669fv5758+e3+073kImI\nRB6/yb++vp7k5OSWz3a7nfr6+g7rNjY2snXrVubMmdPyncViYerUqeTk5PD444+HKGQREekpq7+F\nFosl4BVt2rSJiRMnMnTo0Jbvdu7cSVJSEgcPHmTatGlkZGQwadKkbv+GiIi06snIit+ev81mw+Px\ntHz2eDzY7fYO627YsOG0IZ+kpCQAEhMTmTVrFm63u8O2Pp/P9HL33Xebvq5g2gVSt6s6nS0P5vtQ\n/neLhG0XLdsv2GWRuu2icftFyr+9nvKb/HNycqitraWuro6mpiYqKiooKCg4rd7Ro0d5+eWXmTlz\nZst3jY2NNDQ0AHDs2DG2bdtGZmZmjwMOl7y8PNPXFUy7QOp2Vaez5cF+b7ZQxxUN2y/YZZG67SD6\ntl+f+bfn68LmzZt9F198sW/kyJG+e++91+fz+XyPPfaY77HHHmup87vf/c43f/78du327dvnGzNm\njG/MmDG+Sy65pKXtqQIIQSLY3XffbXYI0k3adtGtp7nT76WevcFisYTkEEbM4XK5IrpXKZ3Ttotu\nPc2dSv4iIlGop7lT0zuIiMQgJX8RkRik5C8iEoOU/EVEYpCSv4hIDFLyFxGJQUr+IiIxSMlfRCQG\nKfmLiMQgJX8RkRik5C8iEoOU/EVEYpCSv4hIDFLyFxGJQUr+IiIxSMlfRCQGKfmLiMQgq9kBAJSX\nw4kTRvH5On49dXmgdTr7Ltg2ncXl77WjAh1/FwiLpfW1benXr/W1o9K/f+ur1dr6emqJi2stZ5zR\nvpx5ZmsZMKD1dcAAGDjQKPHxxmv//qH/f0REQisikv+zz7ZPYG2T2snS9vOp7y2W1gTXUb22y/r3\n77hNV+voKK7O4j01OZ9aoOPv/Olsx3Hqzsbrbb/T8npbX0+W5majtH1//LhRGhuhqcl4f/L1669P\nL1991Vq+/NIojY3G6xlnGDuC+Hg46ywYNMh4HTy49fXss1vLkCFGGTq0tfzLvxjfaUciEh56hq+E\nlM9n7ByOHWtfvvjCKA0NreXzz+HoUeP1yBHj/WeftZaGBmMHcM45Rhk2rLUkJsJ558G55xqv559v\nvI+LM/u/gEjv0APcpc/yeo2dwOHD8OmncOiQUQ4eNMonn8CBA63l4EHjqCEpCWy21pKc3FouuMA4\n+hCJdkr+It/weo2dwz//CfX1rcXjaS1//7sxHHXhhZCSAiNHwkUXGSUtzdg5aKhJokHYk39VVRUr\nVqzA6/WyePFibr/99nbLV61axVNPPQVAc3Mze/fu5dChQwwdOrTLtqH4A0SC4fMZRwh1dbB/v1E+\n/NAotbXGshEjID0dMjKM4nDAqFHGuQqRSBHW5O/1eklPT6e6uhqbzUZubi7l5eU4HI4O6z/33HM8\n/PDDVFdXB9xWyV8iSWOjsSN4/33YuxfefRf27DFezz0XRo+GrCwYM8YoaWk6UhBz9DR3+r3ax+12\nk5qaSkpKCgBOp5PKyspOk//69euZP39+t9qKRIL4eMjMNEpbXq9xlPD22/DXv0JFBfzkJ8a5hm99\nCy69FHJyYNw4Y4fQT3fQSITzm/zr6+tJTk5u+Wy326mpqemwbmNjI1u3buXRRx8Nuq1IpOvfH1JT\njTJrVuv3R47AG2/Arl3w/PPw858bVy3l5sKECUa57DKdZJbI4zf5WwK5AP0bmzZtYuLEiQwdOjTo\ntitXrmx5n5eXR15eXsBtRcw0dChMmWKUkw4cgFdfhVdeMXYGb74Jl1zSWm/CBJ0/kOC5XC5cLlfI\n1uc3+dtsNjweT8tnj8eD3W7vsO6GDRtahnyCbds2+YtEu/POg5kzjQLGjW9uN7z4Itx3n3GUMHYs\nTJ8O+fnGe503kK6c2jEuLi7u0fr8nvBtbm4mPT2d7du3M3z4cMaNG9fhSdujR49y0UUX8Y9//IOB\nAwcG1VYnfCXWNDbCyy/Dtm2wdatxhdF3vwvXXGPsEM4+2+wIJRqE9YSv1WqltLSU/Px8vF4vhYWF\nOBwOysrKACgqKgLgz3/+M/n5+S2J319bkVgXHw8zZhgFjMtOn38e1q6FwkKYOBFmz4aCAuNOZpFw\n0E1eIhGkoQE2b4ZnnjGOCsaNg/nzjZ3BN6fTRADd4SvSZzU2GkcE5eWwfTtMnQo33GAcMWgOI1Hy\nF4kBR47Axo3w298aN6H94AewdClcfLHZkYlZepo7dSuKSBQYOhQWL4adO42TxVYrTJoE3/kOPP20\nMfW2SDDU8xeJUl9/DX/6E6xZYxwNLF9uHA0kJJgdmfQG9fxFYtSZZ4LTCS+9BJs2GXMRjRwJP/qR\nMRWFiD9K/iJ9QHY2PPGEsQM4+2xjnqHrrzcmpRPpiJK/SB9y/vnGXcT79hlTUU+ZAv/6r8ZOQaQt\nJX+RPmjIELjzTuNcQHY2XHGFcSSwb5/ZkUmkUPIX6cPOOgvuuAM++MB4WlluLtx6q/FoTIltSv4i\nMWDIEFi50jgH0NhoPKHsoYd0iWgs06WeIjFozx74t38znmm8erVxv4BEF93hKyLd4vNBZSWsWGHM\nIfTQQ2CzmR2VBErX+YtIt1gscO21xlFAerrxOMqyMjhxwuzIpDeo5y8iALzzjjGFxBlnwOOPGzsE\niVzq+YtISIwebcwdNHeu8ajJRx7RUUBfpp6/iJymtta4L2DQIGMm0eRksyOSU6nnLyIhl5YG//M/\ncOWVcOmlxnTS0reo5y8ifr3+ujFFxPTp8OCD0OZprWIi9fxFJKxycmD3bvjsM7j8cnjvPbMjklBQ\n8heRLg0ZYjxO8kc/Mh4i8+c/mx2R9JSGfUQkKG63cUXQ9ddDcTH07292RLFJd/iKSK/75BPjPMDA\ngcYRwZAhZkcUezTmLyK97txz4YUXjJlCv/1tPTksGnWZ/KuqqsjIyCAtLY2SkpIO67hcLrKzsxk9\nejR5eXkt36ekpJCVlUV2djbjxo0LWdAiYj6rFUpLYdkyYwewc6fZEUkw/A77eL1e0tPTqa6uxmaz\nkZubS3l5OQ6Ho6XOkSNHmDBhAlu3bsVut3Po0CGGDRsGwIgRI9i1axcJfp4orWEfkehXVWWcA1i9\n2hgOkvAL67CP2+0mNTWVlJQU4uLicDqdVFZWtquzfv165syZg91uB2hJ/CcpsYv0fTNmQHW18aCY\nRx4xOxoJhN/kX19fT3Kb+7rtdjv19fXt6tTW1nL48GGmTJlCTk4O69ata1lmsViYOnUqOTk5PP74\n4yEOXUQiSVaWMfTz6KPG08PU74tsVn8LLRZLlys4fvw4u3fvZvv27TQ2NjJ+/Hguv/xy0tLS2LFj\nB8OHD+fgwYNMmzaNjIwMJk2adNo6Vq5c2fI+Ly+v3XkDEYkeF14IO3bANddAYaExO6guBQ0Nl8uF\ny+UK2fr8Jn+bzYbH42n57PF4WoZ3TkpOTmbYsGEMHDiQgQMHcsUVV/DWW2+RlpbG8OHDAUhMTGTW\nrFm43e4uk7+IRLdhw2D7dpg5ExYsgHXrIC7O7Kii36kd4+Li4h6tz++wT05ODrW1tdTV1dHU1ERF\nRQUFBQXt6sycOZMdO3bg9XppbGykpqaGUaNG0djYSENDAwDHjh1j27ZtZGZm9ihYEYkOgwbBc8/B\nF1/AvHnw9ddmRySn8tvzt1qtlJaWkp+fj9frpbCwEIfDQVlZGQBFRUVkZGQwY8YMsrKy6NevH0uW\nLGHUqFHs27eP2bNnA9Dc3MyCBQuYPn16+P8iEYkIAwbAH/8IP/gBFBQYU0JoUrjIoTt8RSSsmpuN\ny0APHzZ2AAMGmB1R36A7fEUkolmt8PvfG1NAzJmjIaBIoeQvImFntcKTTxrDPvPmQVOT2RGJkr+I\n9Iq4OGMSuH79jPMAXq/ZEcU2JX8R6TVxcbBhA3z6qTEnkE73mUfJX0R61YABxonfN9+EO+80O5rY\npeQvIr1u8GDYsgU2bYL77zc7mtjk9zp/EZFwOecc2LbNmA7abofrrjM7otii5C8iprHZ4Pnn4cor\nISkJpkwxO6LYoWEfETHV6NFQUWE8B+Cdd8yOJnYo+YuI6aZMgYcegquvhn/+0+xoYoOGfUQkIixY\nAPv2GbOBvvSS5gEKN83tIyIRw+drPfG7fj0E8EiRmKW5fUSkz7BYYO1a4wjgnnvMjqZv07CPiESU\ngQONm8AuuwwcDpg71+yI+iYN+4hIRNq9G/Lz4cUXjSuCpD0N+4hInzR2LDzwAMyaBUeOmB1N36Oe\nv4hEtJtvhv37obLSmBFUDOr5i0if9sADcPQo/OIXZkfSt6jnLyIR7+OPIScHfvMb4zyAqOcvIjHg\n/PPhqafghhvgH/8wO5q+QclfRKLC5MnG+L/TCcePmx1N9NOwj4hEjRMn4JprjEs/Y/05AD3NnUr+\nIhJVDh0yLgN99FFjRxCrwj7mX1VVRUZGBmlpaZSUlHRYx+VykZ2dzejRo8nLywuqrYhIMIYNM8b/\nFy/WDKA94bfn7/V6SU9Pp7q6GpvNRm5uLuXl5TgcjpY6R44cYcKECWzduhW73c6hQ4cYNmxYQG1B\nPX8R6Z6774ZXXoGtW2Pz+v+w9vzdbjepqamkpKQQFxeH0+mksrKyXZ3169czZ84c7HY7AMOGDQu4\nrYhId/3sZ/Dll7BqldmRRCe/E7vV19eTnJzc8tlut1NTU9OuTm1tLcePH2fKlCk0NDRwyy238MMf\n/jCgtietXLmy5X1eXl67oSMRkY5YrcbwT26u8TCY3FyzIwovl8uFy+UK2fr8Jn9LAJNpHz9+nN27\nd7N9+3YaGxsZP348l19+eUBtT2qb/EVEAnXhhVBaajwI5o03YNAgsyMKn1M7xsXFxT1an99hH5vN\nhsfjafns8XhahndOSk5OZvr06QwcOJBzzjmHK664grfeeiugtiIiPfX97xvTP992m9mRRBe/yT8n\nJ4fa2lrq6upoamqioqKCgoKCdnVmzpzJjh078Hq9NDY2UlNTw6hRowJqKyISCqtXw6ZNxslfCYzf\nYR+r1UppaSn5+fl4vV4KCwtxOByUlZUBUFRUREZGBjNmzCArK4t+/fqxZMkSRo0aBdBhWxGRUBs6\nFH77W1i4EN56C845x+yIIp9u8hKRPmPFCuPa/4oKsyMJP03sJiLyjfvug7/+FZ5+2uxIIp96/iLS\np7z6Klx7Lbz9NiQmmh1N+GhuHxGRU9x2G3z0Ud8e/tGwj4jIKYqLjRO/GzeaHUnkUs9fRPqk//1f\nmD3bOAfQF4d/NOwjItKJW2+FAwfgySfNjiT0lPxFRDpx7BhkZsJ//zdcdZXZ0YSWxvxFRDoxaBCU\nlcGyZfDFF2ZHE1nU8xeRPm/hQkhIgIceMjuS0NGwj4hIFz79FC65BCorjUng+gIN+4iIdOGcc+DB\nB2HJEjh+3OxoIoOSv4jEhPnz4fzz4ZFHzI4kMmjYR0RiRm0tjB8Pu3fDBReYHU3PaNhHRCRAaWlw\n881wyy1mR2I+JX8RiSm33w579sCzz5odibk07CMiMWf7digshL/9LXqf+6thHxGRIH3nO/Dtb8O9\n95odiXnU8xeRmPR//wdZWcYEcGlpZkcTPPX8RUS6YfhwuOMO4+RvLPY/lfxFJGbdfDPs3w+bNpkd\nSe/TsI+IxLTqali61Dj5O3Cg2dEETsM+IiI9MHUqXHop3H+/2ZH0LvX8RSTmffQRjB0Lb74Jyclm\nRxOYsPf8q6qqyMjIIC0tjZKSktOWu1wuhgwZQnZ2NtnZ2fzyl79sWZaSkkJWVhbZ2dmMGzeu20GK\niITThRfC8uXGDWCxwupvodfrZfny5VRXV2Oz2cjNzaWgoACHw9Gu3uTJk3m2g9vlLBYLLpeLhISE\n0EYtIhJit90GGRmwYwdMnGh2NOHnt+fvdrtJTU0lJSWFuLg4nE4nlZWVp9Xzd+ihIR0RiQaDBkFJ\nCaxYASdOmB1N+Pnt+dfX15PcZgDMbrdTU1PTro7FYuGVV15hzJgx2Gw2Vq1axahRo1qWTZ06lf79\n+1NUVMSSJUs6/J2VK1e2vM/LyyMvL6+bf46ISPfNn2887/eJJ2DRIrOjac/lcuFyuUK2Pr/J32Kx\ndLmCsWPH4vF4iI+PZ8uWLVx77bW8//77AOzcuZOkpCQOHjzItGnTyMjIYNKkSaeto23yFxExi8UC\nDz8MM2fC3LkweLDZEbU6tWNcXFzco/X5Hfax2Wx4PJ6Wzx6PB7vd3q7O4MGDiY+PB+Cqq67i+PHj\nHD58GICkpCQAEhMTmTVrFm63u0fBioiEW26ucflnX7/002/yz8nJoba2lrq6OpqamqioqKCgoKBd\nnQMHDrSM67vdbnw+HwkJCTQ2NtLQ0ADAsWPH2LZtG5mZmWH6M0REQuc//gMefRTa9H37HL/DPlar\nldLSUvLz8/F6vRQWFuJwOCgrKwOgqKiIjRs3smbNGqxWK/Hx8WzYsAGAjz/+mNmzZwPQ3NzMggUL\nmD59epj/HBGRnktOhmXL4Kc/hd//3uxowkM3eYmIdKChAS6+2Jj3JyfH7GhOp+kdRETCYPBg+MUv\n4NZb++asn0r+IiKduPFGOHwYOri9Kepp2EdExI/Nm43e/9tvg9XvWdLepWEfEZEwuuoqSEqCtWvN\njiS01PMXEenC668bN369/37kPPBdPX8RkTDLyYFJk+Chh8yOJHTU8xcRCcCHH8Jll8HevZCYaHY0\nPc+dSv4iIgG6+Wbj9ZFHzI0DlPxFRHrNwYPGnP+vvw4jRpgbi8b8RUR6SWIi3HQT9IWJiNXzFxEJ\nwuefQ1oabN8Oo0ebF4d6/iIivejss+GOO4xJ36KZev4iIkH66itj0reKChg/3pwY1PMXEellAwYY\n4/533BG9k74p+YuIdMP118OBA/DCC2ZH0j1K/iIi3WC1QnEx3HVXdPb+lfxFRLpp3jz4+mvjgS/R\nRslfRKSb+vWDX/4SfvYzOHHC7GiCo+QvItID3/senHkmbNxodiTB0aWeIiI9tG2bMe/PO+/03gNf\ndKmniIjJpk2Dc8+Fp54yO5LAqecvIhICL70EhYXw7ru90/tXz19EJAJMngwXXADr1pkdSWC6TP5V\nVVVkZGSQlpZGSUnJactdLhdDhgwhOzub7Oxs7rnnnoDbioj0JcXFxtU/x4+bHUnX/A77eL1e0tPT\nqa6uxmazkZubS3l5OQ6Ho6WOy+XiwQcf5Nlnnw26LWjYR0T6lqlTYf58YwgonMI67ON2u0lNTSUl\nJYW4uDicTieVlZWn1esogEDbioj0JcXFcM890NRkdiT++T0tUV9fT3Jycstnu91OTU1NuzoWi4VX\nXnmFMWPGYLPZWLVqFaNGjQqo7Ukr2zwZIS8vj7y8vG78KSIi5pswwZjv/3e/g6VLQ7del8uFy+UK\n2fr8Jn+LxdLlCsaOHYvH4yE+Pp4tW7Zw7bXX8v777wcVxMq+8FgcEZFvFBeD0wk33ABnnBGadZ7a\nMS4uLu7R+vwO+9hsNjweT8tnj8eD3W5vV2fw4MHEx8cDcNVVV3H8+HEOHz6M3W7vsq2ISF80frzx\nrN8nnjA7ks75Tf45OTnU1tZSV1dHU1MTFRUVFBQUtKtz4MCBljF/t9uNz+cjISEhoLYiIn3Vz38O\n994buVf++B32sVqtlJaWkp+fj9frpbCwEIfDQVlZGQBFRUVs3LiRNWvWYLVaiY+PZ8OGDX7biojE\nggkTYORI47r/G280O5rT6Q5fEZEwefllWLQI3nsv9Hf96g5fEZEIdcUVxl2/kTjnj3r+IiJh9OKL\nxiWfe/eGtvevnr+ISATLy4OkJPjmdGjEUM9fRCTMtm2DFSuM+f77hajLrZ6/iEiEmzYNzjoL/vQn\nsyNppeQvIhJmFgvcdZcx50+kDHQo+YuI9IJrrjEe8r55s9mRGJT8RUR6Qb9+8NOfRk7vX8lfRKSX\nzJkDn31mXP5pNiV/EZFe0r8/3Hmn0fs3m5K/iEgvuu46+PBDePVVc+NQ8hcR6UVxcXDbbXDffebG\noZu8RER62ZdfwogRUF0No0d3bx26yUtEJMoMHGjc8furX5kXg3r+IiImOHrUmO/f7YaLLgq+vXr+\nIiJRaMgQKCqC++835/fV8xcRMcnBg5CeDn/7mzHzZzDU8xcRiVKJibBgAfzXf/X+b6vnLyJioro6\nuPRS2LfPGAoKlHr+IiJRLCUFvvtdWLOmd39XPX8REZO9/TZMnw7798OAAYG1Uc9fRCTKZWZCTg48\n8UTv/aZ6/iIiEWDHDli4EN57L7AHvYe9519VVUVGRgZpaWmUlJR0Wu+1117DarXyzDPPtHyXkpJC\nVlYW2dnZjBs3rttBioj0dRMnGpd7tkmhYeW35+/1eklPT6e6uhqbzUZubi7l5eU4HI7T6k2bNo34\n+HgWLVrEnDlzABgxYgS7du0iISGh8wDU8xcRAWDTJiguhtdeMx796E9Ye/5ut5vU1FRSUlKIi4vD\n6XRSWVl5Wr3Vq1czd+5cEhMTT1umxC4iEpirr4bGxt552IvfkaX6+nqSk5NbPtvtdmpqak6rU1lZ\nyV/+8hdee+01LG12VxaLhalTp9K/f3+KiopYsmRJh7+zcuXKlvd5eXnk5eV1408REYlu/frBv/+7\nMeXDlVe2X+ZyuXC5XCH7Lb/J39LVcQewYsUKfvWrX7UcgrTt6e/cuZOkpCQOHjzItGnTyMjIYNKk\nSaeto23yFxGJZdddB3fdBW+9BWPGtH5/ase4uLi4R7/jN/nbbDY8Hk/LZ4/Hg91ub1dn165dOJ1O\nAA4dOsSWLVuIi4ujoKCApG8mq0hMTGTWrFm43e4Ok7+IiBjOPBNuuQX+8z/hySfD9zt+x/xzcnKo\nra2lrq6OpqYmKioqKCgoaFdn37597N+/n/379zN37lzWrFlDQUEBjY2NNDQ0AHDs2DG2bdtGZmZm\n+P4SEZE+oqgItmyBjz4K32/47flbrVZKS0vJz8/H6/VSWFiIw+GgrKzsmwCLOm378ccfM3v2bACa\nm5tZsGAB06dPD2HoIiJ905AhUFgIDz4YvknfdJOXiEgEqq837vz94APo6Gp5Te8gItIH2WxQUACP\nPRae9avnLyISoU5O+FZXZ5wIbks9fxGRPiozE771LXjqqdCvWz1/EZEItn073HQTvPOOcRPYSer5\ni4j0YVdeaczxv2VLaNer5C8iEsEsFvjxj2HVqtCuV8lfRCTCzZtnPOP39ddDt04lfxGRCBcXZ0z5\n8MADoVunTviKiESBzz+HESPgjTfgggt0wldEJCacfTYsWhS66R7U8xcRiRJ//ztkZxvj/0OHqucv\nIhITLrgA8vPh17/u+brU8xcRiSKvvw6zZ4PHo56/iEjMyMkxTvz2lJK/iEiUue++nq9Dwz4iIlFI\nl3qKiEjQlPxFRGKQkr+ISAxS8hcRiUFK/iIiMUjJX0QkBnWZ/KuqqsjIyCAtLY2SkpJO67322mtY\nrVaeeeaZoNtK9HK5XGaHIN2kbRfb/CZ/r9fL8uXLqaqqYs+ePZSXl7N3794O691+++3MmDEj6LYS\n3ZRAope2XWzzm/zdbjepqamkpKQQFxeH0+mksrLytHqrV69m7ty5JCYmBt02UoTyH0J31xVMu0Dq\ndlWns+XBfm+2UMcVDdsv2GWRuu0g+rZfX/m35zf519fXk5yc3PLZbrdTX19/Wp3KykqWLVsGGHed\nBdo2kij5d/97s0Vb8gi0rpJ/764v1pI/Pj82btzoW7x4ccvndevW+ZYvX96uzty5c32vvvqqz+fz\n+RYuXOjbuHFjwG2/mVpCRUVFRaUbpSes+GGz2fB4PC2fPR4Pdru9XZ1du3bhdDoBOHToEFu2bCEu\nLi6gthjR+wtBRETCwG/yz8nJoba2lrq6OoYPH05FRQXl5eXt6uzbt6/l/aJFi/je975HQUEBzc3N\nXbYVERFz+E3+VquV0tJS8vPz8Xq9FBYW4nA4KCsrA6CoqCjotiIiYj7Tp3QWEZHepzt8RURiUMQl\n/2PHjrFw4UKWLl3K+vXrzQ5HgrR//34WL17MvHnzzA5FglRZWcnSpUtxOp288MILZocjQXr33XdZ\ntmwZ3//+9/nNb37TZf2IG/ZZt24dCQkJXH311TidTjZs2GB2SNIN8+bN4+mnnzY7DOmGI0eO8OMf\n/5hf//rXZoci3XDixAmcTid/+MMf/NbrlZ7/jTfeyHnnnUdmZma77zua+6ftzWH9+/fvjfCkC8Fs\nP4ks3dl299xzD8uXL+/NMKUTwW6/TZs2tXScu9SjuwQC9PLLL/t2797tGz16dMt3zc3NvpEjR/r2\n79/va2pq8o0ZM8a3Z88e37p163zPPfecz+fz+ZxOZ2+EJ10IZvudNHfuXDNClVMEs+1OnDjhu+22\n23zV1dUmRixtdeffns/n8xUUFHS5br+XeobKpEmTqKura/dd27l/gJa5f26++WaWL1/O888/T0FB\nQW+EJ10IZvudd955/OQnP+HNN9+kpKSE22+/vfcDlhbBbLvq6mq2b9/O559/zgcffOD3Um7pHcFs\nv08++YQ//vGPfPXVV0yZMqXLdfdK8u9IR3P/1NTUEB8fz9q1a80KSwLU2fZLSEjgscceMzEy6Upn\n22716tWnvD7WAAAA1ElEQVTcdNNNJkYmgehs+02ePJnJkycHvB7TrvY5OQGcRCdtv+ilbRfdQrX9\nTEv+gc79I5FJ2y96adtFt1BtP9OSf9t5g5qamqioqNAYfxTR9ote2nbRLWTbL+SnpzvgdDp9SUlJ\nvjPOOMNnt9t9a9eu9fl8Pt/mzZt9F198sW/kyJG+e++9tzdCkW7Q9ote2nbRLZzbL+Ju8hIRkfCL\nuOkdREQk/JT8RURikJK/iEgMUvIXEYlBSv4iIjFIyV9EJAYp+YuIxCAlfxGRGKTkLyISg/4fl/Jh\nyoxrDRsAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 256
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.linear_model import LassoCV\n",
      "lasso_cv = LassoCV(alphas=np.logspace(-1, 0, 100), cv=10)\n",
      "lasso_regressor = lasso_cv.fit(X_train, y_train)\n",
      "\n",
      "print '{} out of {} non zero coefficients with ||w||_2 = {}'.format(\n",
      "    np.sum(lasso_regressor.coef_>=.01),\n",
      "    len(lasso_regressor.coef_),\n",
      "    np.linalg.norm(lasso_regressor.coef_, ord=2))\n",
      "\n",
      "print lasso_regressor.score(X=X_train, y=y_train), lasso_regressor.score(X=X_test, y=y_test)\n",
      "\n",
      "print get_significant_coef(lasso_regressor.coef_, explanatory_cols, second_order_combinations)\n",
      "\n",
      "plt.plot(lasso_regressor.alphas_, np.mean(lasso_regressor.mse_path_, axis=1))\n",
      "plt.xscale('log')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "5 out of 13 non zero coefficients with ||w||_2 = 6.32070740393\n",
        "0.761425081673 0.681285149439\n",
        "['CRIM * CRIM = -0.610457952752', 'CRIM * ZN = 0.996977254132', 'CRIM * CHAS = 0.640027183162', 'CRIM * NOX = -1.71312404649', 'CRIM * RM = 2.4446399333', 'CRIM * DIS = -2.64752441419', 'CRIM * RAD = 0.77895668647', 'CRIM * TAX = -0.656158990126', 'CRIM * PTRATIO = -2.29294815659', 'CRIM * LSTAT = -3.98054863351']\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEECAYAAADXg6SsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHCZJREFUeJzt3Xl01OW9x/H3oPEiBougECVUUaROEkJCWBXCRDaxFKOA\niBc3wtXidoNLFYUrWA2KYiC44oJeLAhUCXgtaRWYQFAhagjgwiKkjkBahaKMRYjkd/94MIrAkGVm\nfjO/+bzO+Z0TM8kzXw7Hz3n4zrO4LMuyEBERR2hkdwEiIhI8CnUREQdRqIuIOIhCXUTEQRTqIiIO\nolAXEXGQgKHu8/nIysoiOTmZlJQUCgoKACgvL6dHjx6kpqYyePBg9u7dG5ZiRUQkMFegdeqVlZVU\nVlaSlpaG3+8nIyODwsJCrr32Wp544gl69erFrFmz2LZtGw8++GA46xYRkaMIOFNPSEggLS0NgPj4\neNxuN9u3b2fz5s306tULgL59+/L666+HvlIRETmuWvfUKyoqKCsro1u3biQnJ7No0SIAFixYgM/n\nC1mBIiJSB1Yt7N2718rIyLAWLlxoWZZlffbZZ1b//v2tjIwMa9KkSVaLFi2O+nuAHj169Oip49MQ\nx/3tAwcOWP3797fy8/OP+vrGjRutrl27HjPUJXI88MADdpcQNtHyZ42EOsNVQyjfJ5hjB2OshozR\n0NwM2H6xLIucnBySkpLIzc2t+f5XX30FQHV1NQ899BBjxowJNIxECI/HY3cJYRMtf9ZIqDNcNYTy\nfYI5djDGsvPvNeDql5KSEjIzM0lNTcXlcgGQl5fH5s2beeqppwAYMmQIeXl5Rx/c5SLA8CIi8gsN\nzc2Aod5QCnURkbppaG5qR6mIiIMo1EVEHEShLiLiIAp1EREHUaiLiDiIQl1ExEEU6iIiDqJQFxFx\nEIW6iIiDKNRFRBxEoS4i4iAKdRERB1Goi4g4iEJdRMRBFOoiIg6iUBcRcRCFuoiIgyjURUQcRKEu\nIuIgCnUREQdRqIuIOIhCXUQkAuzYASNHNnwchbqIiI0OHIApUyA1FX7964aPp1AXEbGBZcH8+ZCU\nBCtWwPvvQ15ew8c9seFDiIhIXRQXw913w8GD8Nxz0KdP8MZWqIuIhMnq1TBhAmzZAg8/DMOHQ6Mg\n90vUfhERCbHychg8GIYONc/GjTBiRPADHRTqIiIhs2kTXHUVDBhgWiybN8ONN0JcXOjeM2Co+3w+\nsrKySE5OJiUlhYKCAgDWrFlD165dSU9Pp0uXLpSWloauQhGRKLNlC4waBRdeaFa1bNkC//3f0Lhx\n6N/bZVmWdawXKysrqaysJC0tDb/fT0ZGBoWFhYwZM4Zx48YxYMAAlixZwpQpU1i+fPmRg7tcBBhe\nRMRRPvvM9MqXLIFbboHcXDjttLqN0dDcDDhTT0hIIC0tDYD4+Hjcbjfbt2/nzDPP5JtvvgFgz549\ntG7dut4FiIhEu02bzMahzEy44AL4/HOYNKnugR4MAWfqP1dRUUHv3r35+OOP2bVrFz179sTlclFd\nXc17771HmzZtjhxcM3URcbDPP4c//hHeesvMym+/HZo2bdiYDc3NWi1p9Pv9DB06lOnTpxMfH092\ndjYFBQVcfvnlLFiwgFGjRvH2228f9XcnTpxY87XH48Hj8dS7WBGRSPDFF/DQQ/DGG3DrraZn/qtf\n1W8sr9eL1+sNWm3HnalXVVUxaNAgBg4cSG5uLgCnnnoq3377LQCWZdGsWbOadsxhg2umLiIO4vOZ\nLf1z5phVLHfdBS1aBPc9QtpTtyyLnJwckpKSagIdoF27dhQXFwOwbNky2rdvX+8CREQi3ZYtMHo0\ndOxoVrB88glMnhz8QA+GgDP1kpISMjMzSU1NxeVyAZCXl8cZZ5zBLbfcwv79+zn55JN5+umnSU9P\nP3JwzdRFJIqVlsJjj8Hy5XDzzaZnHuogb2hu1vqD0noNrlAXkShz8CD85S8wdSps2wZjx0JOTsM/\nAK2tsHxQKiLidF9/DS++CM88A61amc1Cw4aFdvdnKOiYABGJaRs3wk03wfnnw6efwoIF5uCtq6+O\nvkAHzdRFJAZZFni9kJ9vzjEfM8aEe8uWdlfWcAp1EYkZ330Hf/oTzJhheue33w6vvQZNmthdWfAo\n1EXE8TZsMJdRzJkDPXvCtGlw8cVwaFGfoyjURcSR/H5zXdwLL5gdoDk5UFYWnHtAI5mWNIqIY1gW\nlJTAyy+bLfy9e5sjcC+9FE6MkimsljSKSMzbvBlefRVmz4aTT4brrjMrWRIS7K4s/BTqIhKV/v53\n016ZN8+cyTJihFmO2KmTM3vltaX2i4hEBcsys+/CQvNs3QqXX26ui+vdO3raK8ejYwJExLH27YPi\nYrNtf8kS2L8fsrPN06tXdG4OOh6Fuog4xrffwpo1sGIFrFwJH3wA6ekwcKB5OnZ0fmtFoS4iUenr\nr6G83Dwffmgen8+EeGameS68EE491e5Kw0uhLiIRbf9+c4fn+vWwbp15ysvN7s7UVDP77tQJOncG\nt9s5vfH6UqiLSET55htYvNg869aZjT9t20Jysgnw1FTznH2281sp9aFQFxHbVVfDm2/C88+bfnhW\nllmZ0rUrtGsHJ51kd4XRQ6EuIrbZtw9eeQWeeAKaNYPbboPLLou9PngwaUepiISdzwdPP20uleje\n3Zyv0quX2imRQJdkiEitbdkCV15peuP79sGqVaZ3npmpQI8UCnUROa7vv4cHHzSz8k6dzBb9adPM\nbUESWdR+EZGA3nvPHJCVkgIffeT8o2ujnUJdRI5p7lxzAfPMmWZrvkQ+hbqIHMGy4KGHzAegS5dC\nhw52VyS1pVAXkcNUVcF//Rd8/LG5lPnMM+2uSOpCoS4iNfx+GDrUnH7o9cIpp9hdkdSVVr+ICAD/\n+Ad4PNCmDSxcqECPVgp1EaGiAi66CAYNMh+KxvqhWtFMf3UiMW7bNnNWy113wa232l2NNJRCXSSG\nbd0KF18M99wDY8bYXY0EQ8D2i8/nIysri+TkZFJSUigoKABg+PDhpKenk56eTtu2bUlPTw9LsSIS\nPJ9/bmbo48Yp0J0k4Ew9Li6O/Px80tLS8Pv9ZGRk0K9fP+bNm1fzM3fddRfNmjULeaEiEjyVldCv\nH9x3H9x0k93VSDAFDPWEhAQSEhIAiI+Px+12s2PHDtxuNwCWZTF//nyWL18e+kpFJCj8fvOB6PXX\nK9CdqNY99YqKCsrKyujWrVvN91auXEmrVq0477zzjvl7EydOrPna4/Hg8XjqVaiINFxVFQwbZu4B\nnTDB7moEwOv14vV6gzZerS7J8Pv9eDwexo8fT/bPDoAYM2YM7du3Z+zYsUcfXJdkiEQMy4LRo03r\nZdEiLVuMVCG/JKOqqoohQ4YwcuTIwwL9hx9+YOHChXz00Uf1fnMRCQ/Lgrvvhg0bzFkuCnTnCvhX\na1kWOTk5JCUlkZube9hr77zzDm63m7POOiukBYpIwz30EPztb2brf3y83dVIKAVc0rhq1SpeffVV\nli9fXrOEsaioCIB58+YxYsSIsBQpIvU3fTrMnm1CvXlzu6uRUNPF0yIONnMm5OXBihW63CJa6OJp\nETmqp56Cxx4zPXQFeuxQqIs40LRpUFBgeujnnGN3NRJOCnURh5kyxbRdvF7N0GORQl3EISwL7r/f\nnIXu9UJiot0ViR0U6iIOcPAg3HILfPghrFwJp59ud0ViF4W6SJSrqoJrrjE3Fy1dCqeeandFYieF\nukgUq642B3N9+y0sWQKNG9tdkdhNoS4SpSzL3FT05ZdQVKRAF0OhLhKl7r8fSktNy+Xkk+2uRiKF\nQl0kCj3yiDlpsbhYPXQ5nEJdJMo8/ji89JJZtqhVLvJLCnWRKDJtGjz7rAl0HZAqR6NQF4kSTz5p\nTlwsLtbGIjk2hbpIFJgxA6ZO1dZ/OT6FukiEmzbNzNB1OJfUhkJdJII98YRpu3i9cPbZdlcj0UCh\nLhKhHn/8pw9F1XKR2lKoi0SgKVPg+ed12qLUnUJdJMI8+ii88IIJ9Nat7a5Goo1CXSSCPPLITxuL\nFOhSHwp1kQjx4IMwZ442FknDKNRFbGZZ8D//Y24sKi6GVq3srkiimUJdxEaWBffea47OXb4czjjD\n7ook2inURWxSXW3OQ1+zBpYtgxYt7K5InEChLmKDqiq44Qbw+Uyg6/hcCRaFukiY7dsHw4ebmXpR\nkS64kOBqZHcBIrFk927o3x+aNjUfjCrQJdgU6iJh8ve/Q8+e0L07zJ4NcXF2VyROFDDUfT4fWVlZ\nJCcnk5KSQkFBQc1rM2bMwO12k5KSwj333BPyQkWiWXk5XHQR3HQTPPYYNNJ0SkIkYE89Li6O/Px8\n0tLS8Pv9ZGRk0K9fPyorK1m8eDHr1q0jLi6Or776Klz1ikSdt96C66+Hp5+GYcPsrkacLmCoJyQk\nkJCQAEB8fDxut5vt27fz/PPPM27cOOIO/fvxDC2uFTmCZZlz0KdMgcWLoUcPuyuSWFDr1S8VFRWU\nlZXRrVs37r77blasWMF9991H48aNefzxx+ncufNRf2/ixIk1X3s8HjweT0NrFol4VVVw222wahW8\n957OQpdj83q9eL3eoI3nsizLOt4P+f1+PB4P48ePJzs7mw4dOnDxxRczffp0SktLGT58OFu3bj1y\ncJeLWgwv4ij/+AdceaVZ4TJnjtagS900NDeP+3FNVVUVQ4YMYeTIkWRnZwOQmJjIFVdcAUCXLl1o\n1KgRu3btqncRIk5RWgpdukDv3qblokCXcAsY6pZlkZOTQ1JSErm5uTXfz87OZtmyZQBs2rSJAwcO\n0EJ7nCXGvfwy/Pa3UFBgTlzUChexQ8D2S0lJCZmZmaSmpuJyuQCYPHkyffr0YdSoUaxdu5aTTjqJ\nqVOnHrVXrvaLxIKqKrjzTvjrX6GwENxuuyuSaNbQ3KxVT73egyvUxeG++sr0z5s0gT/9CZo1s7si\niXYh76mLyNH92D/v0cP0zxXoEgl0oJdIHVkWPPecudji2Wfh0JoBkYigUBepg+++gzFjoKwMSkqg\nfXu7KxI5nNovIrW0fr1pt7hc8P77CnSJTAp1keOwLHjmGbj4YnP13CuvwCmn2F2VyNGp/SISwDff\nwOjRsGWL2fKv2blEOs3URY5h7Vro3BlatjTntyjQJRoo1EWO4oUXoF8/szP0qaegcWO7KxKpHbVf\nRH5m/364/XZYuRJWrNDuUIk+CnWRQ3buhKFDTbtl9WpzyqJItFH7RQRYswa6djWXQr/+ugJdopdm\n6hLzZs+GO+6A55+HQ6dLi0QthbrErB9+MOvOCwvB64XkZLsrEmk4hbrEpF274OqrobratF6aN7e7\nIpHgUE9dYk5pKWRkQGoqLFmiQBdn0UxdYoZlmb75+PE6XVGcS6EuMWHvXnO64tq1Zg36b35jd0Ui\noaH2izheWZlpt5x8sumfK9DFyRTq4liWBTNmmLXnEyea1kuTJnZXJRJaar+II1VWwqhR8PXX8O67\ncP75dlckEh6aqYvjLF4M6emm5bJqlQJdYotm6uIYe/bA2LFQXAwLFkDPnnZXJBJ+mqmLIxQVQYcO\n5sPQdesU6BK7NFOXqLZ7N9x1FyxdCrNmQd++dlckYi/N1CUqWRbMn2/OaznlFNiwQYEuApqpSxSq\nqIDbboOtW+GNN6BHD7srEokcmqlL1DhwACZPNveGdu8OH32kQBf5Jc3UJSq88465Zq5tW7Mr9Nxz\n7a5IJDIFnKn7fD6ysrJITk4mJSWFgoICACZOnEhiYiLp6emkp6dTVFQUlmIl9mzdCpdfDjfeCA8/\nDP/3fwp0kUBclmVZx3qxsrKSyspK0tLS8Pv9ZGRkUFhYyPz582natCl33HFH4MFdLgIML3JM334L\njzwCzz1nbiW6805o3NjuqkRCr6G5GbD9kpCQQEJCAgDx8fG43W62b98OoLCWkPjhB3NGy6RJcMkl\nUF4OiYl2VyUSPWr9QWlFRQVlZWV0794dgBkzZtCxY0dycnLYs2dPyAqU2GBZZiVLhw7w5z+byyte\nflmBLlJXAdsvP/L7/Xg8HsaPH092djb//Oc/OeOMMwCYMGECO3fu5MUXXzxycJeLBx54oOa/PR4P\nHo8neNWLIyxdCvfdB/v3Q14eDBwILpfdVYmEh9frxev11vz3pEmTGtQJOW6oV1VVMWjQIAYOHEhu\nbu4Rr1dUVPC73/2O9evXHzm4euoSwMqV8MAD8MUX8Mc/wvDh0EiLbCXGNTQ3A/4vZFkWOTk5JCUl\nHRboO3furPl64cKFdOjQod4FSOxZtQr69YNrr4WRI+HTT2HECAW6SDAEnKmXlJSQmZlJamoqrkP/\nHs7Ly2Pu3LmsXbsWl8tF27Ztee6552jVqtWRg2umLodYlmmzPPyw2RE6bhxcfz2cdJLdlYlElobm\nZq166vUeXKEe8w4eNOebP/oofPONCfMRIyAuzu7KRCJTSJc0itTX99/D//4vTJ0Kv/oV/OEPZhPR\nCSfYXZmIsynUJah27oRnnjGbhrp0gZkzITNTq1lEwkUfTUmDWRa8/z5cc405CnfXLlixwmzp791b\ngS4STuqpS739+9/w2mvw1FPmKrkxYyAnB047ze7KRKKXPiiVsCsvN22VuXPhoovg5pthwAAtSRQJ\nBn1QKmHxr3+ZEJ81CyorYfRoE+5t2thdmYj8nGbqckxVVfC3v8Hs2eYslksugRtuMBuHtIpFJDTU\nfpGgsixYvRrmzIF58+C888yuz6uugubN7a5OxPnUfpEGsyxYv960V157zezyvPpqePddE+oiEj0U\n6jHKsmDDBliwAObPN5uFhg83x9+mpWkZoki0UvslhlRXQ2mpCe6FC81Rt8OGwZVXmo1CCnIR+6n9\nIgHt2wfLlpnzV95802zZv+IK02rp1ElBLuI0CnUH+uIL+Mtf4K23oLgY0tNh8GDweqF9e7urE5FQ\nUvslylkWbNtmzigvLjbPnj1m+eFvfwv9+2vVikg00ZLGGLJvH3zyiVmpsm4dlJXB2rUQHw/du5tz\nVnr3NuevaHenSHRSqDvQv/4FmzbBxo3mVqBPPzVh7vOZ9kmHDuZJTzfPoetiRcQBFOpRyrJgxw4z\n2/7kExPgGzeaMN+3z4T3b34DbjckJZmnXTtdLiHidAr1KLFzp1lOWFoKH3wAH31klhimp0NKignw\nH4P8zDO1KkUkVinUI9DevWYGvnq1edasge++M2vBf3w6dYLWrRXeInI4hbrNqqtNz3vlSrMC5YMP\nzJLClBTo2tV8gNmtm9lurwAXkeNRqNvgyy+hqMicXOj1mkshevY0Z4t37Wr63+p9i0h9KNTDwLJM\nO6Ww0Dw7dpj135dcAn37wlln2V2hiDiFQj1EDhwwG3kWLzZPXBxcfjlkZ5uWis4TF5FQ0NkvQXTg\ngLkUYs4cs83e7Tbb6996y2zoUU9cRCJdzM/UDx40H3LOmwd//rNZUvif/2kOvWrVyu7qRCTWaKZe\nD9XV8N57PwV5y5bmLPHSUjjnHLurExGpv5gJdcsy68XnzTMXQ5x6qgny5cvN7FxExAkcHeqWZXZu\nzptnbvf5j/8wQV5UZHrkIiJOE/AsP5/PR1ZWFsnJyaSkpFBQUHDY61OnTqVRo0bs3r07pEXWRXW1\nuVvzzjvh3HNNiJ94IixaBJ99Bg8+qEAXEecKOFOPi4sjPz+ftLQ0/H4/GRkZ9OvXD7fbjc/n4+23\n3+bss88OV63H9P33sHSpCe4334QWLWDIELOmPDVVq1ZEJHYEDPWEhAQSEhIAiI+Px+12s2PHDtxu\nN3fccQdTpkzhsssuC0uhv/Tllz/d7uP1QseOcNll8Ic/mNMMRURiUa176hUVFZSVldGtWzcWLVpE\nYmIiqampoaztMPv3Q0mJ6Yf/9a+wfbvZ0Tl8OLz0kpmdi4jEulqFut/vZ+jQoUyfPp1GjRqRl5fH\n22+/XfN6oDWVEydOrPna4/Hg8XhqVdjBg1BeDu+8Y1or775rDsm65BKYOdOcdKhdnSIS7bxeL16v\nN2jjHXfzUVVVFYMGDWLgwIHk5uayfv16+vbtS5MmTQD48ssvad26NWvWrKFly5aHD16HRfQHD5pr\n2oqLzTLDFSvM5p++faFPH/B4oFmz+v0hRUSiRUjPfrEsi+uuu44WLVqQn59/1J9p27YtH374Ic2P\ncrtxoOL+/W9zTO2qVWZH57vvmsshMjMhK8uE+KF2vohIzAhpqJeUlJCZmUlqaiquQ0tI8vLyGDhw\nYM3PnHvuuXzwwQcBQ92yoKIC3n//p2fDBtNO6dEDevUyzy8m+iIiMSfiT2m89FKL0lKzVrxHD3Nh\nRPfu0LkzHOrgiIjIIREf6m+8YdGli65uExGpjYgP9Ug/pVFEJJI0NDcDHhMgIiLRRaEuIuIgCnUR\nEQdRqIuIOIhCXUTEQRTqIiIOolAXEXEQhbqIiIMo1EVEHEShLiLiIAp1EREHUaiLiDiIQl1ExEEU\n6iIiDqJQFxFxEIW6iIiDKNRFRBxEoS4i4iAKdRERB1Goi4g4iEJdRMRBFOoiIg6iUBcRcRCFuoiI\ngyjURUQcRKEuIuIgAUPd5/ORlZVFcnIyKSkpFBQUADBhwgQ6duxIWloaffr0wefzhaVYEREJzGVZ\nlnWsFysrK6msrCQtLQ2/309GRgaFhYUkJibStGlTAGbMmEF5eTkvvPDCkYO7XAQYXkREfqGhuRlw\npp6QkEBaWhoA8fHxuN1uduzYURPoAH6/n9NPP73eBUj4eL1eu0sIm2j5s0ZCneGqIZTvE8yxgzGW\nnX+vte6pV1RUUFZWRrdu3QC4//77+fWvf80rr7zCvffeG7ICJXgiIUDCJVr+rJFQp0I9+GPZ+vdq\n1cLevXutjIwMa+HChUe8NnnyZOv6668/6u8BevTo0aOnjk9DBOypA1RVVTFo0CAGDhxIbm7uEa9/\n8cUXXHrppWzYsCHQMCIiEgYB2y+WZZGTk0NSUtJhgb558+aarxctWkR6enroKhQRkVoLOFMvKSkh\nMzOT1NRUXC4XAHl5ebz44ots3LiRE044gfPOO49nnnmGli1bhq1oERE5uuO2X0REJHpoR6mIiIOE\nPdS3bdvG6NGjGTZsWLjfWkQk6nz33Xdcd9113HjjjcyZM+e4Px/2UG/btu1Rd5+KiMiR3njjDa68\n8kpmzpzJ4sWLj/vz9Q71UaNG0apVKzp06HDY94uKirjgggs4//zzefTRR+s7vIiIY9UlP7dv306b\nNm0AOOGEE447dr1D/YYbbqCoqOiw7x08eJBbb72VoqIiPvnkE+bOncunn37K7NmzGTt2LDt27Kjv\n24mIOEZd8jMxMbHm0MTq6urjjl3vUO/VqxennXbaYd9bs2YN7dq145xzziEuLo6rrrqKRYsWcc01\n15Cfn89ZZ53F7t27+f3vf8/atWs1kxeRmFSX/Lziiit4/fXXufnmmxk8ePBxxz4xmIX+/J8JAImJ\niaxevfqwn2nevDnPPvtsMN9WRCTqHSs/mzRpwksvvVTrcYL6QemPG5RERKRugpWfQQ311q1bH3Zh\nhs/nIzExMZhvISLiSMHKz6CGeufOndm8eTMVFRUcOHCAefPm1aoHJCIS64KVn/UO9REjRnDhhRey\nadMm2rRpw6xZszjxxBN58sknGTBgAElJSQwfPhy3213ftxARcaRQ5qfOfhERcRCd/SIi4iAKdRER\nB1Goi4g4iEJdRMRBFOoiIg6iUBcRcRCFuoiIgyjURUQcRKEuIuIg/w8gHzOnvzDkLgAAAABJRU5E\nrkJggg==\n"
      }
     ],
     "prompt_number": 258
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}