In this notebook you will compare different regression models in order to assess which model fits best. We will be using polynomial regression as a means to examine this topic. In particular you will:
We will continue to use the House data from previous notebooks.
In [1]:
import graphlab
Next we're going to write a polynomial function that takes an SArray and a maximal degree and returns an SFrame with columns containing the SArray to all the powers up to the maximal degree.
The easiest way to apply a power to an SArray is to use the .apply() and lambda x: functions. For example to take the example array and compute the third power we can do as follows: (note running this cell the first time may take longer than expected since it loads graphlab)
In [2]:
tmp = graphlab.SArray([1., 2., 3.])
tmp_cubed = tmp.apply(lambda x: x**3)
print tmp
print tmp_cubed
We can create an empty SFrame using graphlab.SFrame() and then add any columns to it with ex_sframe['column_name'] = value. For example we create an empty SFrame and make the column 'power_1' to be the first power of tmp (i.e. tmp itself).
In [3]:
ex_sframe = graphlab.SFrame()
ex_sframe['power_1'] = tmp
print ex_sframe
Using the hints above complete the following function to create an SFrame consisting of the powers of an SArray up to a specific degree:
In [16]:
def polynomial_sframe(feature, degree):
# assume that degree >= 1
# initialize the SFrame:
poly_sframe = graphlab.SFrame()
# then loop over the remaining degrees:
for power in range(1, degree+1):
# then assign poly_sframe[name] to the appropriate power of feature
poly_sframe['power_' + str(power)] = feature.apply(lambda x: x**power)
return poly_sframe
To test your function consider the smaller tmp variable and what you would expect the outcome of the following call:
In [17]:
print polynomial_sframe(tmp, 3)
Let's use matplotlib to visualize what a polynomial regression looks like on some real data.
In [20]:
sales = graphlab.SFrame('../Data/kc_house_data.gl/')
sales.head()
Out[20]:
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 [21]:
sales = sales.sort(['sqft_living', 'price'])
Let's start with a degree 1 polynomial using 'sqft_living' (i.e. a line) to predict 'price' and plot what it looks like.
In [25]:
poly1_data = polynomial_sframe(sales['sqft_living'], 1)
print poly1_data
poly1_data['price'] = sales['price'] # add price to the data since it's the target
print poly1_data
NOTE: for all the models in this notebook use validation_set = None to ensure that all results are consistent across users.
In [26]:
model1 = graphlab.linear_regression.create(poly1_data, target = 'price', features = ['power_1'], validation_set = None)
In [27]:
#let's take a look at the weights before we plot
model1.get("coefficients")
Out[27]:
In [28]:
import matplotlib.pyplot as plt
%matplotlib inline
In [29]:
plt.plot(poly1_data['power_1'],poly1_data['price'],'.',
poly1_data['power_1'], model1.predict(poly1_data),'-')
Out[29]:
Let's unpack that plt.plot() command. The first pair of SArrays we passed are the 1st power of sqft and the actual price we then ask it to print these as dots '.'. The next pair we pass is the 1st power of sqft and the predicted values from the linear model. We ask these to be plotted as a line '-'.
We can see, not surprisingly, that the predicted values all fall on a line, specifically the one with slope 280 and intercept -43579. What if we wanted to plot a second degree polynomial?
In [31]:
poly2_data = polynomial_sframe(sales['sqft_living'], 2)
print poly2_data
my_features = poly2_data.column_names() # get the name of the features
print my_features
poly2_data['price'] = sales['price'] # add price to the data since it's the target
print poly2_data
model2 = graphlab.linear_regression.create(poly2_data, target = 'price', features = my_features, validation_set = None)
In [32]:
model2.get("coefficients")
Out[32]:
In [33]:
plt.plot(poly2_data['power_1'],poly2_data['price'],'.',
poly2_data['power_1'], model2.predict(poly2_data),'-')
Out[33]:
The resulting model looks like half a parabola. Try on your own to see what the cubic looks like:
In [34]:
poly3_data = polynomial_sframe(sales['sqft_living'], 3)
print poly3_data
my_features3 = poly3_data.column_names() # get the name of the features
print my_features3
poly3_data['price'] = sales['price'] # add price to the data since it's the target
print poly3_data
model3 = graphlab.linear_regression.create(poly3_data, target = 'price', features = my_features3, validation_set = None)
In [35]:
model3.get("coefficients")
Out[35]:
In [36]:
plt.plot(poly3_data['power_1'],poly3_data['price'],'.',
poly3_data['power_1'], model3.predict(poly3_data),'-')
Out[36]:
Now try a 15th degree polynomial:
In [37]:
poly15_data = polynomial_sframe(sales['sqft_living'], 15)
print poly15_data
my_features15 = poly15_data.column_names() # get the name of the features
print my_features15
poly15_data['price'] = sales['price'] # add price to the data since it's the target
print poly15_data
model15 = graphlab.linear_regression.create(poly15_data, target = 'price', features = my_features15, validation_set = None)
In [38]:
model15.get("coefficients")
Out[38]:
In [39]:
plt.plot(poly15_data['power_1'],poly15_data['price'],'.',
poly15_data['power_1'], model15.predict(poly15_data),'-')
Out[39]:
What do you think of the 15th degree polynomial? Do you think this is appropriate? If we were to change the data do you think you'd get pretty much the same curve? Let's take a look.
We're going to split the sales data into four subsets of roughly equal size. Then you will estimate a 15th degree polynomial model on all four subsets of the data. Print the coefficients (you should use .print_rows(num_rows = 16) to view all of them) and plot the resulting fit (as we did above). The quiz will ask you some questions about these results.
To split the sales data into four subsets, we perform the following steps:
.random_split(0.5, seed=0). .random_split(0.5, seed=0).We set seed=0 in these steps so that different users get consistent results.
You should end up with 4 subsets (set_1, set_2, set_3, set_4) of approximately equal size.
In [41]:
set_a, set_b = sales.random_split(0.5,seed=0)
set_1, set_2 = set_a.random_split(0.5,seed=0)
set_3, set_4 = set_b.random_split(0.5,seed=0)
In [42]:
print len(set_1), len(set_2), len(set_3), len(set_4)
Fit a 15th degree polynomial on set_1, set_2, set_3, and set_4 using sqft_living to predict prices. Print the coefficients and make a plot of the resulting model.
In [44]:
set_1_15_data = polynomial_sframe(set_1['sqft_living'], 15)
set_2_15_data = polynomial_sframe(set_2['sqft_living'], 15)
set_3_15_data = polynomial_sframe(set_3['sqft_living'], 15)
set_4_15_data = polynomial_sframe(set_4['sqft_living'], 15)
#
my_features_x_15 = set_1_15_data.column_names() # get the name of the features
#
set_1_15_data['price'] = set_1['price'] # add price to the data since it's the target
set_2_15_data['price'] = set_2['price'] # add price to the data since it's the target
set_3_15_data['price'] = set_3['price'] # add price to the data since it's the target
set_4_15_data['price'] = set_4['price'] # add price to the data since it's the target
#
model_1_15 = graphlab.linear_regression.create(set_1_15_data, target='price', features=my_features_x_15, validation_set=None)
model_2_15 = graphlab.linear_regression.create(set_2_15_data, target='price', features=my_features_x_15, validation_set=None)
model_3_15 = graphlab.linear_regression.create(set_3_15_data, target='price', features=my_features_x_15, validation_set=None)
model_4_15 = graphlab.linear_regression.create(set_4_15_data, target='price', features=my_features_x_15, validation_set=None)
In [45]:
model_1_15.get("coefficients").print_rows(num_rows = 16)
model_2_15.get("coefficients").print_rows(num_rows = 16)
model_3_15.get("coefficients").print_rows(num_rows = 16)
model_4_15.get("coefficients").print_rows(num_rows = 16)
In [47]:
plt.plot(set_1_15_data['power_1'],set_1_15_data['price'],'.',set_1_15_data['power_1'],model_1_15.predict(set_1_15_data),'-')
Out[47]:
In [48]:
plt.plot(set_2_15_data['power_1'],set_2_15_data['price'],'.',set_2_15_data['power_1'],model_2_15.predict(set_2_15_data),'-')
Out[48]:
In [49]:
plt.plot(set_3_15_data['power_1'],set_3_15_data['price'],'.',set_3_15_data['power_1'],model_3_15.predict(set_3_15_data),'-')
Out[49]:
In [50]:
plt.plot(set_4_15_data['power_1'],set_4_15_data['price'],'.',set_4_15_data['power_1'],model_4_15.predict(set_4_15_data),'-')
Out[50]:
Whenever we have a "magic" parameter like the degree of the polynomial there is one well-known way to select these parameters: validation set. (We will explore another approach in week 4).
We split the sales dataset 3-way into training set, test set, and validation set as follows:
training_and_validation and testing. Use random_split(0.9, seed=1).training and validation. Use random_split(0.5, seed=1).Again, we set seed=1 to obtain consistent results for different users.
In [51]:
training_and_validation, testing = sales.random_split(0.9,seed=1)
training, validation = training_and_validation.random_split(0.5,seed=1)
Next you should write a loop that does the following:
(Note you can turn off the print out of linear_regression.create() with verbose = False)
In [53]:
def get_residual_sum_of_squares(model, data, outcome):
# First get the predictions
predictions = model.predict(data)
# Then compute the residuals/errors
residuals = predictions - outcome
# Then square and add them up
RSS = sum(pow(residuals,2))
return(RSS)
In [78]:
def minimize_rss( training_data, validation_data, degrees ):
degree_rss = {}
for degree in range(1,degrees+1):
poly_degree = polynomial_sframe(training_data['sqft_living'], degree)
poly_features = poly_degree.column_names()
poly_degree['price'] = training_data['price']
poly_model = graphlab.linear_regression.create(
poly_degree, target='price', features=poly_features, validation_set=None, verbose=False)
poly_validation_data = polynomial_sframe(validation_data['sqft_living'], degree)
degree_rss[degree] = get_residual_sum_of_squares(poly_model, poly_validation_data, validation_data['price'])
print (degree,degree_rss[degree])
min_value = min(degree_rss.values())
min_key = [key for key, value in degree_rss.iteritems() if value == min_value]
return( min_key[0], min_value )
In [79]:
print minimize_rss( training, validation, 15 )
Now that you have chosen the degree of your polynomial using validation data, compute the RSS of this model on TEST data. Report the RSS on your quiz.
In [80]:
print minimize_rss( training, testing, 15 )
In [ ]: