week-5-lasso-assignment-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression Week 5: Feature Selection and LASSO (Interpretation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, you will use LASSO to select features, building on a pre-implemented solver for LASSO (using GraphLab Create, though you can use other solvers). You will:\n",
    "* Run LASSO with different L1 penalties.\n",
    "* Choose best L1 penalty using a validation set.\n",
    "* Choose best L1 penalty using a validation set, with additional constraint on the size of subset.\n",
    "\n",
    "In the second notebook, you will implement your own LASSO solver, using coordinate descent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fire up graphlab create"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import graphlab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load in house sales data\n",
    "\n",
    "Dataset is from house sales in King County, the region where the city of Seattle, WA is located."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO] This non-commercial license of GraphLab Create is assigned to kcbighuge@gmail.com and will expire on October 25, 2016. For commercial licensing options, visit https://dato.com/buy/.\n",
      "\n",
      "[INFO] Start server at: ipc:///tmp/graphlab_server-13031 - Server binary: /Users/jjl/.graphlab/anaconda/lib/python2.7/site-packages/graphlab/unity_server - Server log: /tmp/graphlab_server_1454301823.log\n",
      "[INFO] GraphLab Server Version: 1.8.1\n"
     ]
    }
   ],
   "source": [
    "sales = graphlab.SFrame('kc_house_data.gl/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create new features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in Week 2, we consider features that are some transformations of inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import log, sqrt\n",
    "sales['sqft_living_sqrt'] = sales['sqft_living'].apply(sqrt)\n",
    "sales['sqft_lot_sqrt'] = sales['sqft_lot'].apply(sqrt)\n",
    "sales['bedrooms_square'] = sales['bedrooms']*sales['bedrooms']\n",
    "\n",
    "# In the dataset, 'floors' was defined with type string, \n",
    "# so we'll convert them to float, before creating a new feature.\n",
    "sales['floors'] = sales['floors'].astype(float) \n",
    "sales['floors_square'] = sales['floors']*sales['floors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Squaring bedrooms will increase the separation between not many bedrooms (e.g. 1) and lots of bedrooms (e.g. 4) since 1^2 = 1 but 4^2 = 16. Consequently this variable will mostly affect houses with many bedrooms.\n",
    "* On the other hand, taking square root of sqft_living will decrease the separation between big house and small house. The owner may not be exactly twice as happy for getting a house that is twice as big."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learn regression weights with L1 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us fit a model with all the features available, plus the features we just created above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17\n"
     ]
    }
   ],
   "source": [
    "all_features = ['bedrooms', 'bedrooms_square',\n",
    "            'bathrooms',\n",
    "            'sqft_living', 'sqft_living_sqrt',\n",
    "            'sqft_lot', 'sqft_lot_sqrt',\n",
    "            'floors', 'floors_square',\n",
    "            'waterfront', 'view', 'condition', 'grade',\n",
    "            'sqft_above',\n",
    "            'sqft_basement',\n",
    "            'yr_built', 'yr_renovated']\n",
    "print len(all_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying L1 penalty requires adding an extra parameter (`l1_penalty`) to the linear regression call in GraphLab Create. (Other tools may have separate implementations of LASSO.)  Note that it's important to set `l2_penalty=0` to ensure we don't introduce an additional L2 penalty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PROGRESS: Linear regression:\n",
      "PROGRESS: --------------------------------------------------------\n",
      "PROGRESS: Number of examples          : 21613\n",
      "PROGRESS: Number of features          : 17\n",
      "PROGRESS: Number of unpacked features : 17\n",
      "PROGRESS: Number of coefficients    : 18\n",
      "PROGRESS: Starting Accelerated Gradient (FISTA)\n",
      "PROGRESS: --------------------------------------------------------\n",
      "PROGRESS: +-----------+----------+-----------+--------------+--------------------+---------------+\n",
      "PROGRESS: | Iteration | Passes   | Step size | Elapsed Time | Training-max_error | Training-rmse |\n",
      "PROGRESS: +-----------+----------+-----------+--------------+--------------------+---------------+\n",
      "PROGRESS: Tuning step size. First iteration could take longer than subsequent iterations.\n",
      "PROGRESS: | 1         | 2        | 0.000002  | 1.523467     | 6962915.603493     | 426631.749026 |\n",
      "PROGRESS: | 2         | 3        | 0.000002  | 1.566143     | 6843144.200219     | 392488.929838 |\n",
      "PROGRESS: | 3         | 4        | 0.000002  | 1.611693     | 6831900.032123     | 385340.166783 |\n",
      "PROGRESS: | 4         | 5        | 0.000002  | 1.669655     | 6847166.848958     | 384842.383767 |\n",
      "PROGRESS: | 5         | 6        | 0.000002  | 1.713882     | 6869667.895833     | 385998.458623 |\n",
      "PROGRESS: | 6         | 7        | 0.000002  | 1.755608     | 6847177.773672     | 380824.455891 |\n",
      "PROGRESS: +-----------+----------+-----------+--------------+--------------------+---------------+\n",
      "PROGRESS: TERMINATED: Iteration limit reached.\n",
      "PROGRESS: This model may not be optimal. To improve it, consider increasing `max_iterations`.\n"
     ]
    }
   ],
   "source": [
    "model_all = graphlab.linear_regression.create(sales, target='price', features=all_features,\n",
    "                                              validation_set=None, \n",
    "                                              l2_penalty=0., l1_penalty=1e10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find what features had non-zero weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------------+-------+-------+--------+\n",
      "|       name      | index | value | stderr |\n",
      "+-----------------+-------+-------+--------+\n",
      "|     bedrooms    |  None |  0.0  |  None  |\n",
      "| bedrooms_square |  None |  0.0  |  None  |\n",
      "|     sqft_lot    |  None |  0.0  |  None  |\n",
      "|  sqft_lot_sqrt  |  None |  0.0  |  None  |\n",
      "|      floors     |  None |  0.0  |  None  |\n",
      "|  floors_square  |  None |  0.0  |  None  |\n",
      "|    waterfront   |  None |  0.0  |  None  |\n",
      "|       view      |  None |  0.0  |  None  |\n",
      "|    condition    |  None |  0.0  |  None  |\n",
      "|  sqft_basement  |  None |  0.0  |  None  |\n",
      "|     yr_built    |  None |  0.0  |  None  |\n",
      "|   yr_renovated  |  None |  0.0  |  None  |\n",
      "+-----------------+-------+-------+--------+\n",
      "[12 rows x 4 columns]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model_all.coefficients[model_all.coefficients['value']==0].print_rows(num_rows=len(model_all.coefficients))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\"><table frame=\"box\" rules=\"cols\">\n",
       "    <tr>\n",
       "        <th style=\"padding-left: 1em; padding-right: 1em; text-align: center\">name</th>\n",
       "        <th style=\"padding-left: 1em; padding-right: 1em; text-align: center\">index</th>\n",
       "        <th style=\"padding-left: 1em; padding-right: 1em; text-align: center\">value</th>\n",
       "        <th style=\"padding-left: 1em; padding-right: 1em; text-align: center\">stderr</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">(intercept)</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">274873.05595</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">bathrooms</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">8468.53108691</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">sqft_living</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">24.4207209824</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">sqft_living_sqrt</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">350.060553386</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">grade</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">842.068034898</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">sqft_above</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">20.0247224171</td>\n",
       "        <td style=\"padding-left: 1em; padding-right: 1em; text-align: center; vertical-align: top\">None</td>\n",
       "    </tr>\n",
       "</table>\n",
       "[6 rows x 4 columns]<br/>\n",
       "</div>"
      ],
      "text/plain": [
       "Columns:\n",
       "\tname\tstr\n",
       "\tindex\tstr\n",
       "\tvalue\tfloat\n",
       "\tstderr\tfloat\n",
       "\n",
       "Rows: 6\n",
       "\n",
       "Data:\n",
       "+------------------+-------+---------------+--------+\n",
       "|       name       | index |     value     | stderr |\n",
       "+------------------+-------+---------------+--------+\n",
       "|   (intercept)    |  None |  274873.05595 |  None  |\n",
       "|    bathrooms     |  None | 8468.53108691 |  None  |\n",
       "|   sqft_living    |  None | 24.4207209824 |  None  |\n",
       "| sqft_living_sqrt |  None | 350.060553386 |  None  |\n",
       "|      grade       |  None | 842.068034898 |  None  |\n",
       "|    sqft_above    |  None | 20.0247224171 |  None  |\n",
       "+------------------+-------+---------------+--------+\n",
       "[6 rows x 4 columns]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_all.coefficients[model_all.coefficients['value']!=0][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that a majority of the weights have been set to zero. So by setting an L1 penalty that's large enough, we are performing a subset selection. \n",
    "\n",
    "***QUIZ QUESTION***:\n",
    "According to this list of weights, which of the features have been chosen? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Selecting an L1 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find a good L1 penalty, we will explore multiple values using a validation set. Let us do three way split into train, validation, and test sets:\n",
    "* Split our sales data into 2 sets: training and test\n",
    "* Further split our training data into two sets: train, validation\n",
    "\n",
    "Be *very* careful that you use seed = 1 to ensure you get the same answer!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(training_and_validation, testing) = sales.random_split(.9,seed=1) # initial train/test split\n",
    "(training, validation) = training_and_validation.random_split(0.5, seed=1) # split training into train and validate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we write a loop that does the following:\n",
    "* For `l1_penalty` in [10^1, 10^1.5, 10^2, 10^2.5, ..., 10^7] (to get this in Python, type `np.logspace(1, 7, num=13)`.)\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list.\n",
    "    * Compute the RSS on VALIDATION data (here you will want to use `.predict()`) for that `l1_penalty`\n",
    "* Report which `l1_penalty` produced the lowest RSS on validation data.\n",
    "\n",
    "When you call `linear_regression.create()` make sure you set `validation_set = None`.\n",
    "\n",
    "Note: you can turn off the print out of `linear_regression.create()` with `verbose = False`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "l1_search = np.logspace(1, 7, num=13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def model_l1(train_data=training, valid_data=validation, l1_penalty=10., output_name='price', features_list=all_features):\n",
    "    model = graphlab.linear_regression.create(train_data, target = output_name, features = features_list, \\\n",
    "                                              validation_set = None, verbose = False, l1_penalty=l1_penalty, l2_penalty=0.)\n",
    "    e = compute_rss(model.predict(valid_data[features_list]), valid_data[output_name])\n",
    "    print 'L1 ', l1_penalty, ':', e\n",
    "    return e\n",
    "\n",
    "def compute_rss(pred, actual):\n",
    "    return sum((pred-actual)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L1  10.0 : 6.25766285142e+14\n",
      "L1  31.6227766017 : 6.25766285362e+14\n",
      "L1  100.0 : 6.25766286058e+14\n",
      "L1  316.227766017 : 6.25766288257e+14\n",
      "L1  1000.0 : 6.25766295212e+14\n",
      "L1  3162.27766017 : 6.25766317206e+14\n",
      "L1  10000.0 : 6.25766386761e+14\n",
      "L1  31622.7766017 : 6.25766606749e+14\n",
      "L1  100000.0 : 6.25767302792e+14\n",
      "L1  316227.766017 : 6.25769507644e+14\n",
      "L1  1000000.0 : 6.25776517727e+14\n",
      "L1  3162277.66017 : 6.25799062845e+14\n",
      "L1  10000000.0 : 6.25883719085e+14\n"
     ]
    }
   ],
   "source": [
    "y = []\n",
    "for i in l1_search:\n",
    "    y.append(model_l1(l1_penalty=i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGa1JREFUeJzt3X20XHV97/H3J4mBJARiQEBBEgMEhBgQELRinSJdAcoV\nL4IgDxUqYqvYrrbegiwvOVdZpdd121UBRQIhlWsR5akEKA+Wy0ADRM4hxMNDggkQAklXTCSwIBDN\nw/f+sScwnMyZxz2zZ/b5vNaaxey9f/Pb3x8TPtn89sMoIjAzs/walXUBZmbWXg56M7Occ9CbmeWc\ng97MLOcc9GZmOeegNzPLuY4HvaS5ktZIGqyj7ackPS5pk6STh2y7W9J6SfObqKEg6QlJT0l6YJg2\nP5G0VNKgpGsljS6t/7SkVyUtKr2+XVo/vdTnotI/X5P0l82Or7R9oqSXJF3e6BjNzLbJ4oh+HjCr\nzrYvAl8C/rXCtu8BZ1X7sKQXKqzbBfgBcGJEzABOHebjP4mIAyNiJjAeOK9s20MRcVjpdSlARPw6\nIj4aEYcBhwMbgNuqD6/q+AC+CzxYow8zs6o6HvQRsQBYX75O0rTSEXq/pAclTS+1XRkRTwHb3dUV\nEQ8Ab9TaXYV1ZwC3RMSqUj/rhqnznrLFx4C9y0uusd9jgeci4iVobnySDgd2B+6rsS8zs6q6ZY5+\nDnBBRHwM+B/AVSn1WymQpwOTJT1QCt6zq3YgjQHOBsqD/xOSFku6S9JBFT52GvDTsuWGxidJwP8B\nvjnMGMzM6jYm6wIkTQD+ALipFHAA72mhvyuBT5YW3y9pUen9TRFxGcmYDwOOASYAj0p6NCKWD9Pl\nD4EHI+Lh0vLjwD4R8aak44F/I/nLY9v+3wN8FriohfF9DbgrIlaXPuKwN7OmZR70JP9Xsb40t92y\niLhg23tJz1fo92VgXURsBDZKegg4BNgu6CVdAuwWEeeX9f9G2fu7Jf1Q0uSIeKW0+njg8YhY28L4\nPgEcLelrwETgPZJej4iLG+jDzAyoY+pG0t6S/p+kpyU9OdyVJJIul7SsNKVxaK1uSy8i4nXgBUmn\nlPU1c5jPDNtPlf0MdTtJiI6WNB44Cliy3Qel80hOGn9xyPo9yt4fCags5Cm1f3vappnxRcRZETE1\nIqaRTN9c75A3s6ZFRNUXsCdwaOn9TsCzwIFD2hxPMtUASXAurNLfDcBq4HfASuBcYApwN7AYeAr4\ndqntEcBLwOvAWuDJsn4eAtaQXN2yEvjjCvt6fpgavgk8DQwC3yhbfxewZ+n9JmAZ8ASwqKymr5dq\nfAJ4BDiq7PPjS3VOHLK/hsdX9tkvAZfX+p788ssvv4Z7KaKxxxRL+jfgioi4v2zdj4AHIuJnpeUl\nQCEi1jTUuZmZpa6hq24kTQUOBX45ZNNeJEem26wqrTMzs4zVfTJW0k7AzcBfRdkJyUZI8q+cmJk1\nISKavvquriP60rXkNwP/NyJur9BkFfDBsuW9S+u20+pc0+zZs1tuV2lbrXVDt29brtamm8dXbbme\n9+0cX6Nj65bxteu7S2N8vfRnc6SNr55saVW9UzfXAc9ExPeH2T4f+FMASR8HXo02zc8XCoWW21Xa\nVmvd0O3bluutp16dGl+15WrjblU9/TU6tkrrsxhfu767SuvzNL5G/7zmbXwdyZZaf3uR3Hy0heSK\nkW1XoBwHfBU4v6zdlSTXov8KOGyYviLPZs+enXUJbeXx9a48jy0i/+MrZWfT/7dSc44+kjtCR9fR\n7oJabfIu7aPfbuPx9a48jw3yP75WNXx5ZUs7k6KT+zMzywNJRLtPxpqZWe9y0JuZ5ZyD3sws5xz0\nZmY556A3M8s5B72ZWc456M3Mcs5Bb2aWcw56M7Occ9CbmeWcg97MLOcc9GZmOeegNzPLOQe9mVnO\nOejNzHLOQW9mlnMOejOznHPQm5nlnIPezCznHPRmZl1s+fLW+3DQm5l1sSuuaL0PB72ZWRfr72+9\nj5pBL2mupDWSBofZvrOk+ZIWS3pS0jmtl2VmZps2wWDF5G1MPUf084BZVbZ/HXg6Ig4F/gj4R0lj\nWi/NzGxke/pp2Gef1vupGfQRsQBYX60JMLH0fiLw24jY3HppZmYjW38/fOxjrfeTxhz9lcBBklYD\nvwL+KoU+zcxGvLSCPo0pllnAExFxjKR9gV9ImhkRb1Rq3NfX9/b7QqFAoVBIoQQzs/woFosUi0Xm\nz4etW1vvTxFRu5E0BbgjImZW2HYncFlEPFxavh+4MCIGKrSNevZnZjbSvfUW7LorvPIKjBsnIkLN\n9lXv1I1Kr0peBI4FkLQHMB14vtmCzMwMFi+GD38Ydtyx9b5qTt1IugEoALtKWgnMBsYCERFzgEuB\nfym7/PLvIuKV1kszMxu50pqfhzqCPiLOqLH9v6h++aWZmTWovx/SOoXpO2PNzLpQmkf0dZ2MTYtP\nxpqZ1fbaa7DXXvDqqzBmDEidORlrZmYdsmgRHHJIEvJpcNCbmXWZNKdtwEFvZtZ1HPRmZjnnoDcz\ny7G1a5OTsPvtl16fDnozsy7S3w9HHAGjUkxnB72ZWRdJe9oGHPRmZl3FQW9mlmMRDnozs1x76SWQ\nYO+90+3XQW9m1iW2Hc2r6YcdVOagNzPrEu2YtgEHvZlZ19h2aWXa/PRKM7MusHUrTJ4My5bB+973\n7m1+eqWZWQ4sXw6TJm0f8mlw0JuZdYF2zc+Dg97MrCs46M3Mcq6dQe+TsWZmGdu8OZmfX70adt55\n++0+GWtm1uOefho++MHKIZ8GB72ZWcbaOW0DdQS9pLmS1kgarNKmIOkJSU9JeiDdEs3M8i3zoAfm\nAbOG2yhpF+AHwIkRMQM4NaXazMxGhMyDPiIWAOurNDkDuCUiVpXar0upNjOz3Nu4EZYuhUMOad8+\n0pijnw5MlvSApH5JZ6fQp5nZiLB4MRxwAIwb1759jEmpj8OAY4AJwKOSHo2I5ZUa9/X1vf2+UChQ\nKBRSKMHMrDdVmrYpFosUi8XU9lHXdfSSpgB3RMTMCtsuBHaMiP9VWr4WuDsibqnQ1tfRm5mV+dKX\n4Oij4StfGb5Np66jV+lVye3A0ZJGSxoPHAUsabYgM7ORpN0nYqGOqRtJNwAFYFdJK4HZwFggImJO\nRCyVdC8wCGwB5kTEM22s2cwsF15/HVauhIMPbu9+/AgEM7OMFItw8cXwyCPV2/kRCGZmPaoT0zbg\noDczy4yD3sws5xz0ZmY5tnYtrF8P++/f/n056M3MMjAwAIcfDqM6kMIOejOzDHRq2gYc9GZmmXDQ\nm5nlWEQS9Ecc0Zn9OejNzDrs5Zdh61bYZ5/O7M9Bb2bWYdumbdT0va6NcdCbmXXYwEDn5ufBQW9m\n1nGdPBELfqiZmVlHRcDkyfDss7D77vV9xg81MzPrIcuXw8471x/yaXDQm5l1UKenbcBBb2bWUQ56\nM7OcyyLofTLWzKxDNm+GSZNg1SrYZZf6P+eTsWZmPeKZZ2DvvRsL+TQ46M3MOqSTz7cp56A3M+uQ\nLObnwUFvZtYxWQW9T8aamXXAxo3JHbHr1sH48Y19tu0nYyXNlbRG0mCNdh+TtEnSyc0WY2aWV4OD\nMH164yGfhnqmbuYBs6o1kDQK+Afg3jSKMjPLm6ymbaCOoI+IBcD6Gs2+AdwM/CaNoszM8qarg74W\nSR8APhcRVwEdeoy+mVlvyTLox6TQxz8DF5YtVw37vr6+t98XCgUKhUIKJZiZda/XX4cVK2DGjPra\nF4tFisViavuv66obSVOAOyJiZoVtz297C+wGbADOj4j5Fdr6qhszG3EefBAuuggefbS5z7d61U29\nR/RimCP1iJhWVsw8kr8Qtgt5M7ORKstpG6gj6CXdABSAXSWtBGYDY4GIiDlDmvtw3cxsiP5+OPHE\n7PbvG6bMzNps2jS46y748Ieb+7yfXmlm1sXWrUteBxyQXQ0OejOzNhoYgMMPh1EZpq2D3sysjbI+\nEQsOejOztnLQm5nl3MCAg97MLLdWrYJNm2DKlGzrcNCbmbXJtmkbZfwUMAe9mVmbdMP8PDjozcza\npluC3nfGmpm1QQTsuissWQJ77NFaX74z1sysCz33HOy0U+shnwYHvZlZG3TLtA046M3M2qK/H444\nIusqEg56M7M26KYjep+MNTNL2ebNMGkSvPxy8s9W+WSsmVmXWbIEPvCBdEI+DQ56M7OUddO0DTjo\nzcxS1w0PMivnoDczS1m3HdH7ZKyZWYp+9zuYPBnWroXx49Pp0ydjzcy6yOAg7LdfeiGfBge9mVmK\num3aBhz0Zmap6smglzRX0hpJg8NsP0PSr0qvBZI+kn6ZZma9oSeDHpgHzKqy/XngDyPiEOBS4Jo0\nCjMz6zVvvAEvvAAzZmRdybuNqdUgIhZIGvYXDyNiYdniQmCvNAozM+s1ixYlIT92bNaVvFvac/Tn\nAXen3KeZWU/oxmkbqOOIvl6S/gg4Fzi6Wru+vr633xcKBQqFQlolmJllqr8fjj++9X6KxSLFYrH1\njkrqumGqNHVzR0TMHGb7TOAW4LiIeK5KP75hysxya9994Y474KCD0u23UzdMqfSqVMA+JCF/drWQ\nNzPLs9/+Nrkb9oADsq5kezWnbiTdABSAXSWtBGYDY4GIiDnA/wQmAz+UJGBTRBzZvpLNzLrPwAAc\ndhiMHp11Jdur56qbM2ps/wrwldQqMjPrQd32xMpyvjPWzCwF3XrFDTjozcxS4aA3M8ux1auTxxNP\nnZp1JZU56M3MWrTtaF5NXwDZXg56M7MWdfO0DTjozcxa1t8PRxyRdRXD808Jmpm1IAJ22w2eegre\n//727MM/JWhmlqHnn4dx49oX8mlw0JuZtaDb5+fBQW9m1hIHvZlZzvVC0PtkrJlZk7ZsgUmTYOVK\neO9727cfn4w1M8vI0qWw557tDfk0OOjNzJrUC9M24KA3M2uag97MLOd6Jeh9MtbMrAm//31yInbt\nWpgwob378slYM7MMDA4mPwbe7pBPg4PezKwJvTJtAw56M7OmOOjNzHKul4LeJ2PNzBq0YQO8733w\n6qswdmz79+eTsWZmHbZoEcyY0ZmQT0PNoJc0V9IaSYNV2lwuaZmkxZIOTbdEM7Pu0kvTNlDfEf08\nYNZwGyUdD+wbEfsDXwV+lFJtZmZdKXdBHxELgPVVmpwEXF9q+0tgF0l7pFOemVn3GRjIWdDXYS/g\npbLlVaV1Zma5s349rFkDBx6YdSX1G9PpHfb19b39vlAoUCgUOl2CmVnTBgbgox+F0aPbt49isUix\nWEytv7our5Q0BbgjImZW2PYj4IGI+FlpeSnw6YhYU6GtL680s561aRN85jNw0knwt3/buf126vJK\nlV6VzAf+tFTMx4FXK4W8mVmv+9a3kmfb/PVfZ11JY2pO3Ui6ASgAu0paCcwGxgIREXMi4t8lnSBp\nObABOLedBZuZZeHmm+GWW5Kpm1E9dgeS74w1M6th6VL41Kfgnnvg8MM7v3/fGWtm1kZvvAGf/zz8\n/d9nE/Jp8BG9mdkwIuCMM2DcOJg7F9T0MXVrWj2i7/jllWZmveKKK5Jpm0ceyS7k0+AjejOzCh55\nBD73OVi4EKZNy7YWz9GbmaVszRo47TS47rrsQz4NDnozszKbN8Ppp8M558CJJ2ZdTTo8dWNmVubC\nC+GJJ+Duu9v7mING+GSsmVlKbrsNbrwRHn+8e0I+DT6iNzMDli2DT34S7rwTjjwy62rezSdjzcxa\ntGEDnHwyfOc73RfyafARvZmNaBFw9tnJ82t+/OPuvF7ec/RmZi246ip48kl49NHuDPk0+IjezEas\nhQvhs59Nbo7ab7+sqxme5+jNzJqwdi184QtwzTXdHfJp8BG9mY04W7bArFnJD3xfdlnW1dTmI3oz\nswZdcklyEva73826ks7wyVgzG1Hmz4frr09uihozQhJwhAzTzAyWL4fzzoPbb4fdd8+6ms7x1I2Z\njQhvvpn8UtQll8AnPpF1NZ3lk7FmlnsRcO65sGkT/OQnvXe9vG+YMjOrYc4cGBiAX/6y90I+DT6i\nN7Nc6++HE06Ahx+G6dOzrqY5vrzSzGwY69bBKafA1Vf3bsinoa6gl3ScpKWSfi3pwgrbd5Y0X9Ji\nSU9KOif1Ss3MGrBlC5x5ZvKTgCefnHU12ao5dSNpFPBr4DPAaqAfOD0ilpa1+Rawc0R8S9JuwLPA\nHhGxeUhfnroxs4645BJ46CH4j//o/evlO3Ey9khgWUS8WNrhjcBJwNKyNgFMLL2fCPx2aMibmXXK\nXXclP+w9MND7IZ+Gev4V7AW8VLb8Mkn4l7sSmC9pNbATcFo65ZmZNeaFF+DP/gxuvRX23DPrarpD\nWn/XzQKeiIhjJO0L/ELSzIh4Y2jDvr6+t98XCgUKhUJKJZjZSPfWW8lNURdfnPwsYK8qFosUi8XU\n+qtnjv7jQF9EHFdavgiIiPjfZW3uBC6LiIdLy/cDF0bEwJC+PEdvZm3z5S8nPwv405/m63r5TszR\n9wP7SZoC/BdwOvDFIW1eBI4FHpa0BzAdeL7ZoszMGjV3bvIrUY89lq+QT0PNoI+ILZIuAO4juRxz\nbkQskfTVZHPMAS4F/kXSYOljfxcRr7StajOzMo8/DhddBP/5n7DTTllX0318Z6yZ9bTnnoNjj4Xv\nfQ9OPTXratrDd8aa2YizaVNyVc2sWXDUUfD1r+c35NPgK0zNrGesWJH8xut118H++8P55yfPlt9x\nx6wr624OejPraps2wZ13Js+rGRiAs86C+++Hgw7KurLe4aA3s6607eh93jzYd9/k6P2222DcuKwr\n6z0OejPrGtuO3ufMSR4vfNZZ8ItfwMEHZ11Zb3PQm1nmVqyAa69N5t6nTUuO3m+91UfvaXHQm1km\nNm9+Z+79sceSo/f77oMZM7KuLH8c9GbWUS+++M7R+4c+5KP3TnDQm1nbbd6cPDr46quT320980y4\n914fvXeKg97M2ubFF5Nn0MydC1OnJkfvN98M48dnXdnI4qA3s9REwG9+AwsXJlfOLFyYHL3fcw98\n5CNZVzdyOejNrG4RyQ9ur1iR/MDHihXvvF54ITmCnzAhuRzynHPgppt89N4N/FAzM3tbBLzyyvBB\nvmIF7LBDchJ16tR3XtuWp0yBiROH7d6a1OpDzRz0ZiNIBLz6avUgHz26epDvsktm5Y9YDnqznNu6\nFTZuhDffTF4bNrzzvtq68vWvvQYrVyZhvnVrEtzDhfmkSdmO17bXiV+YStUb2/2KbH508u+wevdV\nT7tG2wxtP9y2Vt93+rVly7tfmzc3ttxIm82bk983rRXWGzYkIb/DDsnc9/jx27+GWz958jvvJ05M\njsanToX3vte/wDTSdPyIfsKEfB/Rd/I/oHr3VU+7RtsMbT/ctlbfd/I1evQ7rzFjqi+n0aZaSJev\nHzcORvmXI0Y0T92YmeWcf2HKzMyqctCbmeWcg97MLOcc9GZmOVdX0Es6TtJSSb+WdOEwbQqSnpD0\nlKQH0i2zNxSLxaxLaCuPr3fleWyQ//G1qmbQSxoFXAnMAg4GvijpwCFtdgF+AJwYETOAU9tQa9fL\n+x82j6935XlskP/xtaqeI/ojgWUR8WJEbAJuBE4a0uYM4JaIWAUQEevSLfMd9X6h1dpV2lZr3dDt\n25bT/gPWqfFVW6427lbV01+jY6u0Povxteu7q7Q+T+Nr9M9r3sbXiWypJ+j3Al4qW365tK7cdGCy\npAck9Us6O60Ch8rzl9FIfw766uvzFBSV1udpfA76YsVtaY6t5g1Tkj4PzIqI80vLZwFHRsRflrW5\nAjgcOAaYADwKnBARy4f05bulzMya0O5n3awC9ilb3ru0rtzLwLqI2AhslPQQcAjwrqBvpVAzM2tO\nPVM3/cB+kqZIGgucDswf0uZ24GhJoyWNB44ClqRbqpmZNaPmEX1EbJF0AXAfyV8McyNiiaSvJptj\nTkQslXQvMAhsAeZExDNtrdzMzOrS0YeamZlZ5/nOWDOznHPQm5nlXOZBL+lDkq6V9POsa2kHSSdJ\nmiPpp5L+OOt60iTpQElXSfq5pD/Pup52kDS+dG/ICVnXkjZJn5b0UOk7/MOs60mbEpdKuryd9/Zk\nQdLRpe/tGkkLarXPPOgj4oWIOC/rOtolIm4v3YPwF8AXsq4nTRGxNCL+AjgN+IOs62mTC4GfZV1E\nmwTwOrADySXSeXMSyeXgvydn44uIBaX/9u4EflyrfepBL2mupDWSBoesr/lgtF7Qwvi+TfI8oK7V\nzNgk/TeSP2z/3slam9Ho+CQdCzwDrAW6/h6QRscXEQ9FxJ8AFwHf6XS9jWriz+cBwMMR8U3gax0t\ntkEt5MoZwA01dxARqb6Ao4FDgcGydaNIbp6aArwHWAwcOORzN6VdSztezYwP+AfgmKxrb9d3V2p3\nZ9b1pz0+4FLgn4B7gduyrr9d3x8wFvh51vW34fs7Ezil9P7GrOtP+7sDPghcXU//qR/RR8QCYP2Q\n1cM+GE3SZElXAYf2wpF+E+P7BvAZ4BRJ53e02AY1MbZPS/q+pB8Bd3W22sY1Or6I+HZE/A3wr8A1\nHS22CU18f/+99N39mOQJtV2t0fEBtwLHSfo+8GDnKm1cE2MD+DIwr57+63kEQhoqPRjtSICIeIVk\n/rqXVRvfFcAVWRSVkmpje5Au/w+oDsOOb5uIuL6jFaWr2vd3G3BbFkWlqNr43gJ6+fxf1T+bEdFX\nb0eZn4w1M7P26lTQ1/NgtF6W5/HleWzg8fW6PI8vtbG1K+jFu69SqOfBaL0kz+PL89jA4/P4ulf7\nxtaGs8c3AKuB3wErgXNL648HngWWARdlfZbb4xtZY/P4PL5ufrV7bH6omZlZzvlkrJlZzjnozcxy\nzkFvZpZzDnozs5xz0JuZ5ZyD3sws5xz0ZmY556A3M8u5/w9FIsbmK0x+vQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111c2bd90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.plot(l1_search, y)\n",
    "plt.xscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*** QUIZ QUESTIONS ***\n",
    "1. What was the best value for the `l1_penalty`?\n",
    "2. What is the RSS on TEST data of the model with the best `l1_penalty`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.0\n"
     ]
    }
   ],
   "source": [
    "best = np.argmin(y)\n",
    "print l1_search[best]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.56983602382e+14\n"
     ]
    }
   ],
   "source": [
    "test_model = graphlab.linear_regression.create(training, target = 'price', features = all_features, \\\n",
    "                                               validation_set = None, verbose = False, l1_penalty=10.0, l2_penalty=0.)\n",
    "print sum((test_model.predict(testing[all_features]) - testing['price'])**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTION***\n",
    "Also, using this value of L1 penalty, how many nonzero weights do you have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------------------+-------+------------------+--------+\n",
      "|       name       | index |      value       | stderr |\n",
      "+------------------+-------+------------------+--------+\n",
      "|   (intercept)    |  None |  18993.4272128   |  None  |\n",
      "|     bedrooms     |  None |  7936.96767903   |  None  |\n",
      "| bedrooms_square  |  None |  936.993368193   |  None  |\n",
      "|    bathrooms     |  None |  25409.5889341   |  None  |\n",
      "|   sqft_living    |  None |  39.1151363797   |  None  |\n",
      "| sqft_living_sqrt |  None |  1124.65021281   |  None  |\n",
      "|     sqft_lot     |  None | 0.00348361822299 |  None  |\n",
      "|  sqft_lot_sqrt   |  None |  148.258391011   |  None  |\n",
      "|      floors      |  None |   21204.335467   |  None  |\n",
      "|  floors_square   |  None |  12915.5243361   |  None  |\n",
      "|    waterfront    |  None |  601905.594545   |  None  |\n",
      "|       view       |  None |  93312.8573119   |  None  |\n",
      "|    condition     |  None |  6609.03571245   |  None  |\n",
      "|      grade       |  None |  6206.93999188   |  None  |\n",
      "|    sqft_above    |  None |  43.2870534193   |  None  |\n",
      "|  sqft_basement   |  None |  122.367827534   |  None  |\n",
      "|     yr_built     |  None |  9.43363539372   |  None  |\n",
      "|   yr_renovated   |  None |  56.0720034488   |  None  |\n",
      "+------------------+-------+------------------+--------+\n",
      "[18 rows x 4 columns]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "test_model.get(\"coefficients\").print_rows(num_rows=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Limit the number of nonzero weights\n",
    "\n",
    "What if we absolutely wanted to limit ourselves to, say, 7 features? This may be important if we want to derive \"a rule of thumb\" --- an interpretable model that has only a few features in them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, you are going to implement a simple, two phase procedure to achive this goal:\n",
    "1. Explore a large range of `l1_penalty` values to find a narrow region of `l1_penalty` values where models are likely to have the desired number of non-zero weights.\n",
    "2. Further explore the narrow region you found to find a good value for `l1_penalty` that achieves the desired sparsity.  Here, we will again use a validation set to choose the best value for `l1_penalty`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "max_nonzeros = 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the larger range of values to find a narrow range with the desired sparsity\n",
    "\n",
    "Let's define a wide range of possible `l1_penalty_values`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "l1_penalty_values = np.logspace(8, 10, num=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, implement a loop that search through this space of possible `l1_penalty` values:\n",
    "\n",
    "* For `l1_penalty` in `np.logspace(8, 10, num=20)`:\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list. When you call `linear_regression.create()` make sure you set `validation_set = None`\n",
    "    * Extract the weights of the model and count the number of nonzeros. Save the number of nonzeros to a list.\n",
    "        * *Hint: `model['coefficients']['value']` gives you an SArray with the parameters you learned.  If you call the method `.nnz()` on it, you will find the number of non-zero parameters!* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def model_coeffs(train_data=training, valid_data=validation, l1_penalty=10., output_name='price', features_list=all_features):\n",
    "    model = graphlab.linear_regression.create(train_data, target = output_name, features = features_list, \\\n",
    "                                              validation_set = None, verbose = False, l1_penalty=l1_penalty, l2_penalty=0.)\n",
    "    e = compute_rss(model.predict(valid_data[features_list]), valid_data[output_name])\n",
    "    coeffs = model['coefficients']['value']\n",
    "    print 'L1 ', l1_penalty, ':', coeffs.nnz() , e\n",
    "    return e, coeffs, coeffs.nnz()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L1  100000000.0 : 18 6.27492659875e+14\n",
      "L1  127427498.57 : 18 6.28210516771e+14\n",
      "L1  162377673.919 : 18 6.29176689541e+14\n",
      "L1  206913808.111 : 18 6.30650082719e+14\n",
      "L1  263665089.873 : 17 6.32940229287e+14\n",
      "L1  335981828.628 : 17 6.3626814023e+14\n",
      "L1  428133239.872 : 17 6.41261198311e+14\n",
      "L1  545559478.117 : 17 6.48983455376e+14\n",
      "L1  695192796.178 : 17 6.60962217696e+14\n",
      "L1  885866790.41 : 16 6.77261520728e+14\n",
      "L1  1128837891.68 : 15 7.01046815867e+14\n",
      "L1  1438449888.29 : 15 7.37850622829e+14\n",
      "L1  1832980710.83 : 13 7.9616310964e+14\n",
      "L1  2335721469.09 : 12 8.69018172894e+14\n",
      "L1  2976351441.63 : 10 9.66925692362e+14\n",
      "L1  3792690190.73 : 6 1.08186759232e+15\n",
      "L1  4832930238.57 : 5 1.24492736032e+15\n",
      "L1  6158482110.66 : 3 1.38416149024e+15\n",
      "L1  7847599703.51 : 1 1.23079472046e+15\n",
      "L1  10000000000.0 : 1 1.22915716064e+15\n"
     ]
    }
   ],
   "source": [
    "y = []\n",
    "nonzeros = []\n",
    "coeffs = []\n",
    "for i in l1_penalty_values:\n",
    "    rss_error, weights, nz_params = model_coeffs(l1_penalty=i)\n",
    "    y.append(rss_error)\n",
    "    nonzeros.append(nz_params)\n",
    "    coeffs.append(weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out of this large range, we want to find the two ends of our desired narrow range of `l1_penalty`.  At one end, we will have `l1_penalty` values that have too few non-zeros, and at the other end, we will have an `l1_penalty` that has too many non-zeros.  \n",
    "\n",
    "More formally, find:\n",
    "* The largest `l1_penalty` that has more non-zeros than `max_nonzero` (if we pick a penalty smaller than this value, we will definitely have too many non-zero weights)\n",
    "    * Store this value in the variable `l1_penalty_min` (we will use it later)\n",
    "* The smallest `l1_penalty` that has fewer non-zeros than `max_nonzero` (if we pick a penalty larger than this value, we will definitely have too few non-zero weights)\n",
    "    * Store this value in the variable `l1_penalty_max` (we will use it later)\n",
    "\n",
    "\n",
    "*Hint: there are many ways to do this, e.g.:*\n",
    "* Programmatically within the loop above\n",
    "* Creating a list with the number of non-zeros for each value of `l1_penalty` and inspecting it to find the appropriate boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "l1_penalty_min = np.max(l1_penalty_values[[i for i in range(len(nonzeros)) if nonzeros[i] >max_nonzeros]])\n",
    "l1_penalty_max = np.min(l1_penalty_values[[i for i in range(len(nonzeros)) if nonzeros[i] <max_nonzeros]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2976351441.63\n",
      "3792690190.73\n"
     ]
    }
   ],
   "source": [
    "print l1_penalty_min\n",
    "print l1_penalty_max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTIONS***\n",
    "\n",
    "What values did you find for `l1_penalty_min` and`l1_penalty_max`? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the narrow range of values to find the solution with the right number of non-zeros that has lowest RSS on the validation set \n",
    "\n",
    "We will now explore the narrow region of `l1_penalty` values we found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l1_penalty_values = np.linspace(l1_penalty_min,l1_penalty_max,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* For `l1_penalty` in `np.linspace(l1_penalty_min,l1_penalty_max,20)`:\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list. When you call `linear_regression.create()` make sure you set `validation_set = None`\n",
    "    * Measure the RSS of the learned model on the VALIDATION set\n",
    "\n",
    "Find the model that the lowest RSS on the VALIDATION set and has sparsity *equal* to `max_nonzero`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L1  2976351441.63 : 10 9.66925692362e+14\n",
      "L1  3019316638.95 : 10 9.74019450085e+14\n",
      "L1  3062281836.27 : 10 9.81188367942e+14\n",
      "L1  3105247033.59 : 10 9.89328342459e+14\n",
      "L1  3148212230.92 : 10 9.98783211266e+14\n",
      "L1  3191177428.24 : 10 1.00847716702e+15\n",
      "L1  3234142625.56 : 10 1.01829878055e+15\n",
      "L1  3277107822.88 : 10 1.02824799221e+15\n",
      "L1  3320073020.2 : 8 1.03461690923e+15\n",
      "L1  3363038217.52 : 8 1.03855473594e+15\n",
      "L1  3406003414.84 : 8 1.04323723787e+15\n",
      "L1  3448968612.16 : 7 1.04693748875e+15\n",
      "L1  3491933809.48 : 7 1.05114762561e+15\n",
      "L1  3534899006.81 : 7 1.05599273534e+15\n",
      "L1  3577864204.13 : 7 1.06079953176e+15\n",
      "L1  3620829401.45 : 6 1.0657076895e+15\n",
      "L1  3663794598.77 : 6 1.06946433543e+15\n",
      "L1  3706759796.09 : 6 1.07350454959e+15\n",
      "L1  3749724993.41 : 6 1.07763277558e+15\n",
      "L1  3792690190.73 : 6 1.08186759232e+15\n"
     ]
    }
   ],
   "source": [
    "y = []\n",
    "nonzeros = []\n",
    "coeffs = []\n",
    "for i in l1_penalty_values:\n",
    "    rss_error, weights, nz_params = model_coeffs(l1_penalty=i)\n",
    "    y.append(rss_error)\n",
    "    nonzeros.append(nz_params)\n",
    "    coeffs.append(weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTIONS***\n",
    "1. What value of `l1_penalty` in our narrow range has the lowest RSS on the VALIDATION set and has sparsity *equal* to `max_nonzeros`?\n",
    "2. What features in this model have non-zero coefficients?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2976351441.63\n"
     ]
    }
   ],
   "source": [
    "best = np.argmin(y)\n",
    "print l1_penalty_values[best]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_max7 = min(np.array(y)[[i for i in range(len(nonzeros)) if nonzeros[i] == max_nonzeros]])\n",
    "y.index(best_max7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3448968612.16\n"
     ]
    }
   ],
   "source": [
    "print l1_penalty_values[y.index(best_max7)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('(intercept)', 222253.19254432785),\n",
       " ('bedrooms', 661.7227177822587),\n",
       " ('bedrooms_square', 0.0),\n",
       " ('bathrooms', 15873.957259267981),\n",
       " ('sqft_living', 32.41022145125964),\n",
       " ('sqft_living_sqrt', 690.1147733133256),\n",
       " ('sqft_lot', 0.0),\n",
       " ('sqft_lot_sqrt', 0.0),\n",
       " ('floors', 0.0),\n",
       " ('floors_square', 0.0),\n",
       " ('waterfront', 0.0),\n",
       " ('view', 0.0),\n",
       " ('condition', 0.0),\n",
       " ('grade', 2899.4202697498786),\n",
       " ('sqft_above', 30.011575302201045),\n",
       " ('sqft_basement', 0.0),\n",
       " ('yr_built', 0.0),\n",
       " ('yr_renovated', 0.0)]"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zip(test_model.get(\"coefficients\")['name'], coeffs[y.index(best_max7)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}