# Regression Week 4: Ridge Regression (interpretation)

In this notebook, we will run ridge regression multiple times with different L2 penalties to see which one produces the best fit. We will revisit the example of polynomial regression as a means to see the effect of L2 regularization. In particular, we will:

• Use a pre-built implementation of regression (GraphLab Create) to run polynomial regression
• Use matplotlib to visualize polynomial regressions
• Use a pre-built implementation of regression (GraphLab Create) to run polynomial regression, this time with L2 penalty
• Use matplotlib to visualize polynomial regressions under L2 regularization
• Choose best L2 penalty using cross-validation.
• Assess the final fit using test data.

We will continue to use the House data from previous notebooks. (In the next programming assignment for this module, you will implement your own ridge regression learning algorithm using gradient descent.)

# Fire up graphlab create

``````

In [1]:

import graphlab

``````

# Polynomial regression, revisited

We build on the material from Week 3, where we wrote the function to produce an SFrame with columns containing the powers of a given input. Copy and paste the function `polynomial_sframe` from Week 3:

``````

In [2]:

def polynomial_sframe(feature, degree):
# assume that degree >= 1
# initialize the SFrame:
poly_sframe = graphlab.SFrame()
# and set poly_sframe['power_1'] equal to the passed feature
poly_sframe['power_1'] = feature
# first check if degree > 1
if degree > 1:
# then loop over the remaining degrees:
# range usually starts at 0 and stops at the endpoint-1. We want it to start at 2 and stop at degree
for power in range(2, degree+1):
# first we'll give the column a name:
name = 'power_' + str(power)
name_left = 'power_' + str(power-1)
# then assign poly_sframe[name] to the appropriate power of feature
poly_sframe[name] = feature * poly_sframe[name_left]

return poly_sframe

``````

Let's use matplotlib to visualize what a polynomial regression looks like on the house data.

``````

In [3]:

import matplotlib.pyplot as plt
%matplotlib inline

``````
``````

This non-commercial license of GraphLab Create for academic use is assigned to gonadarush@gmail.com and will expire on July 07, 2017.

[INFO] graphlab.cython.cy_server: GraphLab Create v2.1 started. Logging: /tmp/graphlab_server_1471725437.log

``````
``````

In [4]:

sales = graphlab.SFrame('kc_house_data.gl/')

``````

As in Week 3, we will use the sqft_living variable. For plotting purposes (connecting the dots), you'll need to sort by the values of sqft_living. For houses with identical square footage, we break the tie by their prices.

``````

In [5]:

sales = sales.sort(['sqft_living','price'])

``````

Let us revisit the 15th-order polynomial model using the 'sqft_living' input. Generate polynomial features up to degree 15 using `polynomial_sframe()` and fit a model with these features. When fitting the model, use an L2 penalty of `1e-5`:

``````

In [6]:

l2_small_penalty = 1e-5

``````

Note: When we have so many features and so few data points, the solution can become highly numerically unstable, which can sometimes lead to strange unpredictable results. Thus, rather than using no regularization, we will introduce a tiny amount of regularization (`l2_penalty=1e-5`) to make the solution numerically stable. (In lecture, we discussed the fact that regularization can also help with numerical stability, and here we are seeing a practical example.)

With the L2 penalty specified above, fit the model and print out the learned weights.

Hint: make sure to add 'price' column to the new SFrame before calling `graphlab.linear_regression.create()`. Also, make sure GraphLab Create doesn't create its own validation set by using the option `validation_set=None` in this call.

``````

In [8]:

poly_sframe = polynomial_sframe(sales['sqft_living'], 15)
my_features = poly_sframe.column_names()
poly_sframe['price'] = sales['price']
model = graphlab.linear_regression.create(poly_sframe, 'price', features=my_features,
validation_set=None, l2_penalty=1e-5)

``````
``````

Linear regression:

--------------------------------------------------------

Number of examples          : 21613

Number of features          : 15

Number of unpacked features : 15

Number of coefficients    : 16

Starting Newton Method

--------------------------------------------------------

+-----------+----------+--------------+--------------------+---------------+

| Iteration | Passes   | Elapsed Time | Training-max_error | Training-rmse |

+-----------+----------+--------------+--------------------+---------------+

| 1         | 2        | 1.027601     | 2662555.737333     | 245656.462164 |

+-----------+----------+--------------+--------------------+---------------+

SUCCESS: Optimal solution found.

``````
``````

In [12]:

model.get('coefficients')

``````
``````

Out[12]:

name
index
value
stderr

(intercept)
None
167924.863999
nan

power_1
None
103.090931763
nan

power_2
None
0.134604574374
nan

power_3
None
-0.000129071378069
nan

power_4
None
5.18929005625e-08
nan

power_5
None
-7.77169405643e-12
nan

power_6
None
1.71144993357e-16
nan

power_7
None
4.51177796874e-20
nan

power_8
None
-4.78838623813e-25
nan

power_9
None
-2.33343500089e-28
nan

[16 rows x 4 columns]Note: Only the head of the SFrame is printed.You can use print_rows(num_rows=m, num_columns=n) to print more rows and columns.

``````

QUIZ QUESTION: What's the learned value for the coefficient of feature `power_1`?

# Observe overfitting

Recall from Week 3 that the polynomial fit of degree 15 changed wildly whenever the data changed. In particular, when we split the sales data into four subsets and fit the model of degree 15, the result came out to be very different for each subset. The model had a high variance. We will see in a moment that ridge regression reduces such variance. But first, we must reproduce the experiment we did in Week 3.

First, split the data into split the sales data into four subsets of roughly equal size and call them `set_1`, `set_2`, `set_3`, and `set_4`. Use `.random_split` function and make sure you set `seed=0`.

``````

In [13]:

(semi_split1, semi_split2) = sales.random_split(.5,seed=0)
(set_1, set_2) = semi_split1.random_split(0.5, seed=0)
(set_3, set_4) = semi_split2.random_split(0.5, seed=0)

``````

Next, fit a 15th degree polynomial on `set_1`, `set_2`, `set_3`, and `set_4`, using 'sqft_living' to predict prices. Print the weights and make a plot of the resulting model.

Hint: When calling `graphlab.linear_regression.create()`, use the same L2 penalty as before (i.e. `l2_small_penalty`). Also, make sure GraphLab Create doesn't create its own validation set by using the option `validation_set = None` in this call.

``````

In [17]:

def fitAndPlot(data, degree, l2):
poly_data = polynomial_sframe(data['sqft_living'], degree)
my_features = poly_data.column_names() # get the name of the features
poly_data['price'] = data['price'] # add price to the data since it's the target
model = graphlab.linear_regression.create(poly_data, target = 'price', l2_penalty=l2, verbose=False,
features = my_features, validation_set = None)
plt.plot(poly_data['power_1'],poly_data['price'],'.',
poly_data['power_1'], model.predict(poly_data),'-')
model.get("coefficients").print_rows(num_rows = 16)

``````
``````

In [18]:

for data in [set_1, set_2, set_3, set_4]:
fitAndPlot(data, 15, l2_small_penalty)

``````
``````

+-------------+-------+--------------------+-------------------+
|     name    | index |       value        |       stderr      |
+-------------+-------+--------------------+-------------------+
| (intercept) |  None |   9306.46480693    |   695453.031576   |
|   power_1   |  None |   585.865810528    |   2868.03758336   |
|   power_2   |  None |  -0.397305881003   |   4.90807132554   |
|   power_3   |  None | 0.000141470892661  |  0.00463227884253 |
|   power_4   |  None | -1.52945968443e-08 | 2.71604166164e-06 |
|   power_5   |  None | -3.79756581598e-13 | 1.06336760596e-09 |
|   power_6   |  None | 5.97481767382e-17  | 2.93741971421e-13 |
|   power_7   |  None | 1.06888530756e-20  | 5.94742259737e-17 |
|   power_8   |  None | 1.59343815508e-25  | 8.91670367395e-21 |
|   power_9   |  None | -6.92834892798e-29 | 1.00025123147e-24 |
|   power_10  |  None | -6.83813298819e-33 | 8.77679004471e-29 |
|   power_11  |  None | -1.6268619866e-37  | 4.22510722907e-33 |
|   power_12  |  None |  2.8511863672e-41  |        nan        |
|   power_13  |  None | 3.79998248991e-45  |        nan        |
|   power_14  |  None | 1.52652624269e-49  | 7.46943531098e-46 |
|   power_15  |  None | -2.33807325958e-53 | 2.64974348654e-50 |
+-------------+-------+--------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+--------------------+-------------------+
|     name    | index |       value        |       stderr      |
+-------------+-------+--------------------+-------------------+
| (intercept) |  None |   -25115.8941917   |        nan        |
|   power_1   |  None |   783.493762459    |        nan        |
|   power_2   |  None |  -0.767759249298   |        nan        |
|   power_3   |  None | 0.000438766331686  |        nan        |
|   power_4   |  None | -1.15169153433e-07 |        nan        |
|   power_5   |  None | 6.84281171027e-12  |        nan        |
|   power_6   |  None | 2.51195165245e-15  |        nan        |
|   power_7   |  None | -2.06440494281e-19 |        nan        |
|   power_8   |  None | -4.59673153569e-23 |        nan        |
|   power_9   |  None | -2.71278478359e-29 |        nan        |
|   power_10  |  None | 6.21818485319e-31  | 1.44613409189e-27 |
|   power_11  |  None | 6.51741444991e-35  |        nan        |
|   power_12  |  None | -9.41317420769e-40 |        nan        |
|   power_13  |  None | -1.02421358996e-42 |        nan        |
|   power_14  |  None | -1.00391102292e-46 |        nan        |
|   power_15  |  None | 1.30113368153e-50  |  4.1518530229e-48 |
+-------------+-------+--------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+--------------------+-------------------+
|     name    | index |       value        |       stderr      |
+-------------+-------+--------------------+-------------------+
| (intercept) |  None |   462426.577972    |   1492853.45808   |
|   power_1   |  None |   -759.251889293   |   7591.23648924   |
|   power_2   |  None |   1.02867011342    |   16.3156355215   |
|   power_3   |  None | -0.000528264572927 |  0.0196235956464  |
|   power_4   |  None | 1.15422924267e-07  | 1.47351150492e-05 |
|   power_5   |  None | -2.26096191111e-12 | 7.28022751943e-09 |
|   power_6   |  None | -2.08214287633e-15 |  2.4207180876e-12 |
|   power_7   |  None | 4.08770804373e-20  | 5.41506001666e-16 |
|   power_8   |  None | 2.57079137663e-23  | 7.95708136404e-20 |
|   power_9   |  None | 1.24311232005e-27  | 7.34819960116e-24 |
|   power_10  |  None | -1.72025897349e-31 | 4.04571966188e-28 |
|   power_11  |  None | -2.96760975077e-35 |        nan        |
|   power_12  |  None | -1.06574910679e-39 |        nan        |
|   power_13  |  None | 2.42635690906e-43  | 1.97115191803e-40 |
|   power_14  |  None | 3.55598631717e-47  | 2.01733173599e-44 |
|   power_15  |  None | -2.8577741082e-51  | 4.89169891569e-49 |
+-------------+-------+--------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+--------------------+-------------------+
|     name    | index |       value        |       stderr      |
+-------------+-------+--------------------+-------------------+
| (intercept) |  None |   -170240.043013   |   1305887.37827   |
|   power_1   |  None |   1247.59037346    |   7944.94142547   |
|   power_2   |  None |   -1.22460914458   |   19.8015393014   |
|   power_3   |  None | 0.000555254626863  |  0.0262699173506  |
|   power_4   |  None | -6.38262295229e-08 | 1.96119990139e-05 |
|   power_5   |  None | -2.20216028435e-11 | 5.98716339152e-09 |
|   power_6   |  None | 4.81834737761e-15  |        nan        |
|   power_7   |  None | 4.21461720665e-19  |        nan        |
|   power_8   |  None | -7.99881032162e-23 |        nan        |
|   power_9   |  None | -1.32365892877e-26 |        nan        |
|   power_10  |  None | 1.60197996533e-31  |        nan        |
|   power_11  |  None | 2.39904325812e-34  | 6.51519869368e-31 |
|   power_12  |  None | 2.33354506252e-38  | 8.53881404629e-35 |
|   power_13  |  None | -1.79874100684e-42 | 1.05086749241e-38 |
|   power_14  |  None | -6.02862621607e-46 | 7.21432423682e-43 |
|   power_15  |  None | 4.39472650429e-50  | 1.92575917591e-47 |
+-------------+-------+--------------------+-------------------+
[16 rows x 4 columns]

``````
``````

In [ ]:

power_1   |  None |   1247.59037346    |   7944.94142547
power_1   |  None |   -759.251889293   |   7591.2364892
power_1   |  None |   783.493762459    |        nan
power_1   |  None |   585.865810528    |   2868.03758336

``````

The four curves should differ from one another a lot, as should the coefficients you learned.

QUIZ QUESTION: For the models learned in each of these training sets, what are the smallest and largest values you learned for the coefficient of feature `power_1`? (For the purpose of answering this question, negative numbers are considered "smaller" than positive numbers. So -5 is smaller than -3, and -3 is smaller than 5 and so forth.)

# Ridge regression comes to rescue

Generally, whenever we see weights change so much in response to change in data, we believe the variance of our estimate to be large. Ridge regression aims to address this issue by penalizing "large" weights. (Weights of `model15` looked quite small, but they are not that small because 'sqft_living' input is in the order of thousands.)

With the argument `l2_penalty=1e5`, fit a 15th-order polynomial model on `set_1`, `set_2`, `set_3`, and `set_4`. Other than the change in the `l2_penalty` parameter, the code should be the same as the experiment above. Also, make sure GraphLab Create doesn't create its own validation set by using the option `validation_set = None` in this call.

``````

In [19]:

for data in [set_1, set_2, set_3, set_4]:
fitAndPlot(data, 15, 1e5)

``````
``````

+-------------+-------+-------------------+-------------------+
|     name    | index |       value       |       stderr      |
+-------------+-------+-------------------+-------------------+
| (intercept) |  None |   530317.024516   |   1046571.66157   |
|   power_1   |  None |   2.58738875673   |   4316.04540175   |
|   power_2   |  None |  0.00127414400592 |   7.38604640294   |
|   power_3   |  None | 1.74934226932e-07 |  0.00697101248391 |
|   power_4   |  None | 1.06022119097e-11 | 4.08731014987e-06 |
|   power_5   |  None | 5.42247604482e-16 | 1.60023804873e-09 |
|   power_6   |  None | 2.89563828343e-20 | 4.42045701357e-13 |
|   power_7   |  None | 1.65000666351e-24 | 8.95014280934e-17 |
|   power_8   |  None | 9.86081528409e-29 | 1.34185472722e-20 |
|   power_9   |  None | 6.06589348254e-33 | 1.50525563305e-24 |
|   power_10  |  None |  3.7891786887e-37 | 1.32079943911e-28 |
|   power_11  |  None | 2.38223121312e-41 | 6.35826905955e-33 |
|   power_12  |  None | 1.49847969215e-45 |        nan        |
|   power_13  |  None | 9.39161190285e-50 |        nan        |
|   power_14  |  None | 5.84523161981e-54 | 1.12405855888e-45 |
|   power_15  |  None | 3.60120207203e-58 | 3.98753951386e-50 |
+-------------+-------+-------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+-------------------+-------------------+
|     name    | index |       value       |       stderr      |
+-------------+-------+-------------------+-------------------+
| (intercept) |  None |   519216.897383   |        nan        |
|   power_1   |  None |   2.04470474182   |        nan        |
|   power_2   |  None |  0.0011314362684  |        nan        |
|   power_3   |  None | 2.93074277549e-07 |        nan        |
|   power_4   |  None | 4.43540598453e-11 |        nan        |
|   power_5   |  None | 4.80849112204e-15 |        nan        |
|   power_6   |  None | 4.53091707826e-19 |        nan        |
|   power_7   |  None | 4.16042910575e-23 |        nan        |
|   power_8   |  None | 3.90094635128e-27 |        nan        |
|   power_9   |  None |  3.7773187602e-31 |        nan        |
|   power_10  |  None | 3.76650326842e-35 | 1.99309029163e-27 |
|   power_11  |  None | 3.84228094754e-39 |        nan        |
|   power_12  |  None | 3.98520828414e-43 |        nan        |
|   power_13  |  None | 4.18272762394e-47 |        nan        |
|   power_14  |  None | 4.42738332878e-51 |        nan        |
|   power_15  |  None | 4.71518245412e-55 | 5.72216504581e-48 |
+-------------+-------+-------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+-------------------+-------------------+
|     name    | index |       value       |       stderr      |
+-------------+-------+-------------------+-------------------+
| (intercept) |  None |   522911.518048   |   2081057.27089   |
|   power_1   |  None |   2.26890421877   |   10582.2830804   |
|   power_2   |  None |  0.00125905041842 |   22.7442095856   |
|   power_3   |  None | 2.77552918155e-07 |  0.0273555493206  |
|   power_4   |  None |  3.2093309779e-11 | 2.05409433488e-05 |
|   power_5   |  None | 2.87573572364e-15 | 1.01487325035e-08 |
|   power_6   |  None | 2.50076112671e-19 |  3.3745127157e-12 |
|   power_7   |  None | 2.24685265906e-23 | 7.54866457855e-16 |
|   power_8   |  None | 2.09349983135e-27 |  1.1092275627e-19 |
|   power_9   |  None | 2.00435383296e-31 | 1.02434864756e-23 |
|   power_10  |  None | 1.95410800249e-35 | 5.63978619118e-28 |
|   power_11  |  None | 1.92734119456e-39 |        nan        |
|   power_12  |  None | 1.91483699013e-43 |        nan        |
|   power_13  |  None | 1.91102277046e-47 | 2.74781158783e-40 |
|   power_14  |  None | 1.91246242302e-51 | 2.81218685883e-44 |
|   power_15  |  None | 1.91699558035e-55 | 6.81909234989e-49 |
+-------------+-------+-------------------+-------------------+
[16 rows x 4 columns]

+-------------+-------+-------------------+-------------------+
|     name    | index |       value       |       stderr      |
+-------------+-------+-------------------+-------------------+
| (intercept) |  None |   513667.087087   |   1726876.91195   |
|   power_1   |  None |   1.91040938244   |   10506.2168015   |
|   power_2   |  None |  0.00110058029175 |   26.1851225533   |
|   power_3   |  None | 3.12753987879e-07 |  0.0347387642353  |
|   power_4   |  None | 5.50067886825e-11 | 2.59344786219e-05 |
|   power_5   |  None | 7.20467557825e-15 | 7.91729394199e-09 |
|   power_6   |  None | 8.24977249384e-19 |        nan        |
|   power_7   |  None | 9.06503223498e-23 |        nan        |
|   power_8   |  None | 9.95683160453e-27 |        nan        |
|   power_9   |  None | 1.10838127982e-30 |        nan        |
|   power_10  |  None | 1.25315224143e-34 |        nan        |
|   power_11  |  None | 1.43600781402e-38 | 8.61555627852e-31 |
|   power_12  |  None |  1.662699678e-42  | 1.12915409685e-34 |
|   power_13  |  None |  1.9398172453e-46 | 1.38964419166e-38 |
|   power_14  |  None |  2.2754148577e-50 | 9.54006460834e-43 |
|   power_15  |  None | 2.67948784897e-54 | 2.54658182183e-47 |
+-------------+-------+-------------------+-------------------+
[16 rows x 4 columns]

``````

These curves should vary a lot less, now that you applied a high degree of regularization.

QUIZ QUESTION: For the models learned with the high level of regularization in each of these training sets, what are the smallest and largest values you learned for the coefficient of feature `power_1`? (For the purpose of answering this question, negative numbers are considered "smaller" than positive numbers. So -5 is smaller than -3, and -3 is smaller than 5 and so forth.)

# Selecting an L2 penalty via cross-validation

Just like the polynomial degree, the L2 penalty is a "magic" parameter we need to select. We could use the validation set approach as we did in the last module, but that approach has a major disadvantage: it leaves fewer observations available for training. Cross-validation seeks to overcome this issue by using all of the training set in a smart way.

We will implement a kind of cross-validation called k-fold cross-validation. The method gets its name because it involves dividing the training set into k segments of roughtly equal size. Similar to the validation set method, we measure the validation error with one of the segments designated as the validation set. The major difference is that we repeat the process k times as follows:

Set aside segment 0 as the validation set, and fit a model on rest of data, and evalutate it on this validation set
Set aside segment 1 as the validation set, and fit a model on rest of data, and evalutate it on this validation set
...
Set aside segment k-1 as the validation set, and fit a model on rest of data, and evalutate it on this validation set

After this process, we compute the average of the k validation errors, and use it as an estimate of the generalization error. Notice that all observations are used for both training and validation, as we iterate over segments of data.

To estimate the generalization error well, it is crucial to shuffle the training data before dividing them into segments. GraphLab Create has a utility function for shuffling a given SFrame. We reserve 10% of the data as the test set and shuffle the remainder. (Make sure to use `seed=1` to get consistent answer.)

``````

In [20]:

(train_valid, test) = sales.random_split(.9, seed=1)
train_valid_shuffled = graphlab.toolkits.cross_validation.shuffle(train_valid, random_seed=1)

``````

Once the data is shuffled, we divide it into equal segments. Each segment should receive `n/k` elements, where `n` is the number of observations in the training set and `k` is the number of segments. Since the segment 0 starts at index 0 and contains `n/k` elements, it ends at index `(n/k)-1`. The segment 1 starts where the segment 0 left off, at index `(n/k)`. With `n/k` elements, the segment 1 ends at index `(n*2/k)-1`. Continuing in this fashion, we deduce that the segment `i` starts at index `(n*i/k)` and ends at `(n*(i+1)/k)-1`.

With this pattern in mind, we write a short loop that prints the starting and ending indices of each segment, just to make sure you are getting the splits right.

``````

In [21]:

n = len(train_valid_shuffled)
k = 10 # 10-fold cross-validation

for i in xrange(k):
start = (n*i)/k
end = (n*(i+1))/k-1
print i, (start, end)

``````
``````

0 (0, 1938)
1 (1939, 3878)
2 (3879, 5817)
3 (5818, 7757)
4 (7758, 9697)
5 (9698, 11636)
6 (11637, 13576)
7 (13577, 15515)
8 (15516, 17455)
9 (17456, 19395)

``````

Let us familiarize ourselves with array slicing with SFrame. To extract a continuous slice from an SFrame, use colon in square brackets. For instance, the following cell extracts rows 0 to 9 of `train_valid_shuffled`. Notice that the first index (0) is included in the slice but the last index (10) is omitted.

``````

In [22]:

train_valid_shuffled[0:10] # rows 0 to 9

``````
``````

Out[22]:

id
date
price
bedrooms
bathrooms
sqft_living
sqft_lot
floors
waterfront

2780400035
2014-05-05 00:00:00+00:00
665000.0
4.0
2.5
2800.0
5900
1
0

1703050500
2015-03-21 00:00:00+00:00
645000.0
3.0
2.5
2490.0
5978
2
0

5700002325
2014-06-05 00:00:00+00:00
640000.0
3.0
1.75
2340.0
4206
1
0

0475000510
2014-11-18 00:00:00+00:00
594000.0
3.0
1.0
1320.0
5000
1
0

0844001052
2015-01-28 00:00:00+00:00
365000.0
4.0
2.5
1904.0
8200
2
0

2658000373
2015-01-22 00:00:00+00:00
305000.0
4.0
2.0
1610.0
6250
1
0

3750603471
2015-03-27 00:00:00+00:00
239950.0
3.0
2.5
1560.0
4800
2
0

2114700540
2014-10-21 00:00:00+00:00
366000.0
3.0
2.5
1320.0
4320
1
0

2596400050
2014-07-30 00:00:00+00:00
375000.0
3.0
1.0
1960.0
7955
1
0

4140900050
2015-01-26 00:00:00+00:00
440000.0
4.0
1.75
2180.0
10200
1
0

view
condition
grade
sqft_above
sqft_basement
yr_built
yr_renovated
zipcode
lat

0
3
8
1660
1140
1963
0
98115
47.68093246

0
3
9
2490
0
2003
0
98074
47.62984888

0
5
7
1170
1170
1917
0
98144
47.57587004

0
4
7
1090
230
1920
0
98107
47.66737217

0
5
7
1904
0
1999
0
98010
47.31068733

0
4
7
1610
0
1952
0
98118
47.52930128

0
4
7
1560
0
1974
0
98001
47.26533057

0
3
6
660
660
1918
0
98106
47.53271982

0
4
7
1260
700
1963
0
98177
47.76407345

2
3
8
2000
180
1966
0
98028
47.76382378

long
sqft_living15
sqft_lot15

-122.28583258
2580.0
5900.0

-122.02177564
2710.0
6629.0

-122.28796
1360.0
4725.0

-122.36472902
1700.0
5000.0

-122.0012452
1560.0
12426.0

-122.27097145
1310.0
6000.0

-122.28506088
1510.0
12240.0

-122.34716948
1190.0
4200.0

-122.36361517
1850.0
8219.0

-122.27022456
2590.0
10445.0

[10 rows x 21 columns]

``````

Now let us extract individual segments with array slicing. Consider the scenario where we group the houses in the `train_valid_shuffled` dataframe into k=10 segments of roughly equal size, with starting and ending indices computed as above. Extract the fourth segment (segment 3) and assign it to a variable called `validation4`.

``````

In [25]:

validation4 = train_valid_shuffled[5818:7758]

``````

To verify that we have the right elements extracted, run the following cell, which computes the average price of the fourth segment. When rounded to nearest whole number, the average should be \$536,234.

``````

In [26]:

print int(round(validation4['price'].mean(), 0))

``````
``````

536234

``````

After designating one of the k segments as the validation set, we train a model using the rest of the data. To choose the remainder, we slice (0:start) and (end+1:n) of the data and paste them together. SFrame has `append()` method that pastes together two disjoint sets of rows originating from a common dataset. For instance, the following cell pastes together the first and last two rows of the `train_valid_shuffled` dataframe.

``````

In [27]:

n = len(train_valid_shuffled)
first_two = train_valid_shuffled[0:2]
last_two = train_valid_shuffled[n-2:n]
print first_two.append(last_two)

``````
``````

+------------+---------------------------+-----------+----------+-----------+
|     id     |            date           |   price   | bedrooms | bathrooms |
+------------+---------------------------+-----------+----------+-----------+
| 2780400035 | 2014-05-05 00:00:00+00:00 |  665000.0 |   4.0    |    2.5    |
| 1703050500 | 2015-03-21 00:00:00+00:00 |  645000.0 |   3.0    |    2.5    |
| 4139480190 | 2014-09-16 00:00:00+00:00 | 1153000.0 |   3.0    |    3.25   |
| 7237300290 | 2015-03-26 00:00:00+00:00 |  338000.0 |   5.0    |    2.5    |
+------------+---------------------------+-----------+----------+-----------+
+-------------+----------+--------+------------+------+-----------+-------+------------+
| sqft_living | sqft_lot | floors | waterfront | view | condition | grade | sqft_above |
+-------------+----------+--------+------------+------+-----------+-------+------------+
|    2800.0   |   5900   |   1    |     0      |  0   |     3     |   8   |    1660    |
|    2490.0   |   5978   |   2    |     0      |  0   |     3     |   9   |    2490    |
|    3780.0   |  10623   |   1    |     0      |  1   |     3     |   11  |    2650    |
|    2400.0   |   4496   |   2    |     0      |  0   |     3     |   7   |    2400    |
+-------------+----------+--------+------------+------+-----------+-------+------------+
+---------------+----------+--------------+---------+-------------+
| sqft_basement | yr_built | yr_renovated | zipcode |     lat     |
+---------------+----------+--------------+---------+-------------+
|      1140     |   1963   |      0       |  98115  | 47.68093246 |
|       0       |   2003   |      0       |  98074  | 47.62984888 |
|      1130     |   1999   |      0       |  98006  | 47.55061236 |
|       0       |   2004   |      0       |  98042  | 47.36923712 |
+---------------+----------+--------------+---------+-------------+
+---------------+---------------+-----+
|      long     | sqft_living15 | ... |
+---------------+---------------+-----+
| -122.28583258 |     2580.0    | ... |
| -122.02177564 |     2710.0    | ... |
| -122.10144844 |     3850.0    | ... |
| -122.12606473 |     1880.0    | ... |
+---------------+---------------+-----+
[4 rows x 21 columns]

``````

Extract the remainder of the data after excluding fourth segment (segment 3) and assign the subset to `train4`.

``````

In [29]:

train4 = train_valid_shuffled[0:5818].append(train_valid_shuffled[7758:])

``````

To verify that we have the right elements extracted, run the following cell, which computes the average price of the data with fourth segment excluded. When rounded to nearest whole number, the average should be \$539,450.

``````

In [30]:

print int(round(train4['price'].mean(), 0))

``````
``````

539450

``````

Now we are ready to implement k-fold cross-validation. Write a function that computes k validation errors by designating each of the k segments as the validation set. It accepts as parameters (i) `k`, (ii) `l2_penalty`, (iii) dataframe, (iv) name of output column (e.g. `price`) and (v) list of feature names. The function returns the average validation error using k segments as validation sets.

• For each i in [0, 1, ..., k-1]:
• Compute starting and ending indices of segment i and call 'start' and 'end'
• Form validation set by taking a slice (start:end+1) from the data.
• Form training set by appending slice (end+1:n) to the end of slice (0:start).
• Train a linear model using training set just formed, with a given l2_penalty
• Compute validation error using validation set just formed
``````

In [42]:

def k_fold_cross_validation(k, l2_penalty, data, output_name, features_list):
n = len(data)
RSS_total = 0.
for i in xrange(k):
start = (n*i)/k
end = (n*(i+1))/k
validation = data[start:end]
train = data[:start].append(data[end:])
model = graphlab.linear_regression.create(train, target=output_name, l2_penalty=l2_penalty,
verbose=False, features=features_list, validation_set=None)
predictions = model.predict(validation)
errors = predictions - validation[output_name]
RSS = (errors * errors).sum()
RSS_total = RSS_total + RSS
return RSS_total / k

``````

Once we have a function to compute the average validation error for a model, we can write a loop to find the model that minimizes the average validation error. Write a loop that does the following:

• We will again be aiming to fit a 15th-order polynomial model using the `sqft_living` input
• For `l2_penalty` in [10^1, 10^1.5, 10^2, 10^2.5, ..., 10^7] (to get this in Python, you can use this Numpy function: `np.logspace(1, 7, num=13)`.)
• Run 10-fold cross-validation with `l2_penalty`
• Report which L2 penalty produced the lowest average validation error.

Note: since the degree of the polynomial is now fixed to 15, to make things faster, you should generate polynomial features in advance and re-use them throughout the loop. Make sure to use `train_valid_shuffled` when generating polynomial features!

``````

In [40]:

train_data = polynomial_sframe(train_valid_shuffled['sqft_living'], 15)
feature_list = train_data.column_names()
train_data['price'] = train_valid_shuffled['price']

``````
``````

In [37]:

import numpy as np

``````
``````

In [46]:

validation_error = float('inf')
l2_min = 0
errors = []
l2_list = np.logspace(1, 7, num=13)
for l2_penalty in l2_list:
new_error = k_fold_cross_validation(10, l2_penalty, train_data, 'price', feature_list)
errors.append(new_error)
if new_error < validation_error:
l2_min = l2_penalty
validation_error = new_error
print l2_min

``````
``````

1000.0

``````

QUIZ QUESTIONS: What is the best value for the L2 penalty according to 10-fold validation?

You may find it useful to plot the k-fold cross-validation errors you have obtained to better understand the behavior of the method.

``````

In [49]:

# Plot the l2_penalty values in the x axis and the cross-validation error in the y axis.
# Using plt.xscale('log') will make your plot more intuitive.
plt.xscale('log')
plt.plot(np.logspace(1, 7, num=13), errors)

``````
``````

Out[49]:

[<matplotlib.lines.Line2D at 0x12ef69390>]

``````

Once you found the best value for the L2 penalty using cross-validation, it is important to retrain a final model on all of the training data using this value of `l2_penalty`. This way, your final model will be trained on the entire dataset.

``````

In [51]:

final_model = graphlab.linear_regression.create(train_data, target='price', l2_penalty=1000,
verbose=False, features=feature_list, validation_set=None)

``````

QUIZ QUESTION: Using the best L2 penalty found above, train a model using all training data. What is the RSS on the TEST data of the model you learn with this L2 penalty?

``````

In [52]:

test_data = polynomial_sframe(test['sqft_living'], 15)
test_data['price'] = test['price']
predictions = final_model.predict(test_data)

``````
``````

In [53]:

test_error = predictions - test['price']

``````
``````

In [54]:

RSS = (test_error * test_error).sum()

``````
``````

In [55]:

RSS

``````
``````

Out[55]:

128780855058449.28

``````
``````

In [ ]:

``````