Example of Regression Analysis Using the Boston Housing Data Set.¶
In [1]:
from sklearn.cross_validation import KFold
from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet, SGDRegressor
import numpy as np
import pylab as pl
In [2]:
from sklearn.datasets import load_boston
boston = load_boston()
In [3]:
print boston.DESCR
In [4]:
print boston.feature_names
In [5]:
print boston.data.shape
print boston.target.shape
In [6]:
np.set_printoptions(precision=2, linewidth=120, suppress=True, edgeitems=7)
In [7]:
print boston.data
In [8]:
# In order to do multiple regression we need to add a column of 1s for x0
x = np.array([np.concatenate((v,[1])) for v in boston.data])
y = boston.target
In [9]:
# First 10 elements of the data
print x[:10]
In [10]:
# First 10 elements of the response variable
print y[:10]
In [11]:
# Create linear regression object
linreg = LinearRegression()
# Train the model using the training sets
linreg.fit(x,y)
Out[11]:
In [12]:
# Let's see predictions for the first 10 instances
print linreg.predict(x[:10])
In [13]:
# Compute RMSE on training data
p = linreg.predict(x)
# Now we can constuct a vector of errors
err = abs(p-y)
# Let's see the error on the first 10 predictions
print err[:10]
In [14]:
# Dot product of error vector with itself gives us the sum of squared errors
total_error = np.dot(err,err)
# Compute RMSE
rmse_train = np.sqrt(total_error/len(p))
print rmse_train
In [15]:
# We can view the regression coefficients
print 'Regression Coefficients: \n', linreg.coef_
In [18]:
%matplotlib inline
def plot_coefficients(model, n_features, feature_names):
pl.barh(range(n_features), model.coef_[:-1], align='center')
pl.yticks(np.arange(n_features), feature_names)
pl.xlabel("Coefficient Value")
pl.ylabel("Feature")
pl.ylim(-1, n_features)
plot_coefficients(linreg, len(boston.feature_names), boston.feature_names)
In [19]:
print linreg.intercept_
In [20]:
# Plot outputs
%matplotlib inline
pl.plot(p, y,'ro')
pl.plot([0,50],[0,50], 'g-')
pl.xlabel('predicted')
pl.ylabel('real')
pl.show()
In [21]:
# Now let's compute RMSE using 10-fold x-validation
n = 10
kf = KFold(len(x), n_folds=n)
xval_err = 0
for train,test in kf:
linreg.fit(x[train],y[train])
p = linreg.predict(x[test])
e = p-y[test]
xval_err += np.sqrt(np.dot(e,e)/len(x[test]))
rmse_10cv = xval_err/n
In [22]:
method_name = 'Simple Linear Regression'
print('Method: %s' %method_name)
print('RMSE on training: %.4f' %rmse_train)
print('RMSE on 10-fold CV: %.4f' %rmse_10cv)
Let's try Ridge Regression:¶
In [23]:
# Create linear regression object with a ridge coefficient 0.5
ridge = Ridge(fit_intercept=True, alpha=0.5)
# Train the model using the training set
ridge.fit(x,y)
Out[23]:
In [24]:
# Compute RMSE on training data
p = ridge.predict(x)
err = p-y
total_error = np.dot(err,err)
rmse_train = np.sqrt(total_error/len(p))
# Compute RMSE using 10-fold x-validation
n = 10
kf = KFold(len(x), n_folds=n)
xval_err = 0
for train,test in kf:
ridge.fit(x[train],y[train])
p = linreg.predict(x[test])
e = p-y[test]
xval_err += np.sqrt(np.dot(e,e)/len(x[test]))
rmse_10cv = xval_err/n
method_name = 'Ridge Regression'
print('Method: %s' %method_name)
print('RMSE on training: %.4f' %rmse_train)
print('RMSE on 10-fold CV: %.4f' %rmse_10cv)
We can try different values of alpha and observe the impact on x-validation RMSE¶
In [25]:
print('Ridge Regression')
print('alpha\t RMSE_train\t RMSE_10cv\n')
alpha = np.linspace(.01,20,50)
t_rmse = np.array([])
cv_rmse = np.array([])
for a in alpha:
ridge = Ridge(alpha=a)
# computing the RMSE on training data
ridge.fit(x,y)
p = ridge.predict(x)
err = p-y
total_error = np.dot(err,err)
rmse_train = np.sqrt(total_error/len(p))
# computing RMSE using 10-fold cross validation
kf = KFold(len(x), n_folds=10)
xval_err = 0
for train, test in kf:
ridge.fit(x[train], y[train])
p = ridge.predict(x[test])
err = p - y[test]
xval_err += np.sqrt(np.dot(err,err)/len(x[test]))
rmse_10cv = xval_err/n
t_rmse = np.append(t_rmse, [rmse_train])
cv_rmse = np.append(cv_rmse, [rmse_10cv])
print('{:.3f}\t {:.4f}\t\t {:.4f}'.format(a,rmse_train,rmse_10cv))
In [26]:
pl.plot(alpha, t_rmse, label='RMSE-Train')
pl.plot(alpha, cv_rmse, label='RMSE_XVal')
pl.legend( ('RMSE-Train', 'RMSE_XVal') )
pl.ylabel('RMSE')
pl.xlabel('Alpha')
pl.show()
To make comparisons across methods easier, let's parametrize the regression methods:¶
In [27]:
a = 0.3
for name,met in [
('linear regression', LinearRegression()),
('lasso', Lasso(fit_intercept=True, alpha=a)),
('ridge', Ridge(fit_intercept=True, alpha=a)),
('elastic-net', ElasticNet(fit_intercept=True, alpha=a))
]:
met.fit(x,y)
p = met.predict(x)
e = p-y
total_error = np.dot(e,e)
rmse_train = np.sqrt(total_error/len(p))
kf = KFold(len(x), n_folds=10)
err = 0
for train,test in kf:
met.fit(x[train],y[train])
p = met.predict(x[test])
e = p-y[test]
err += np.dot(e,e)
rmse_10cv = np.sqrt(err/len(x))
print('Method: %s' %name)
print('RMSE on training: %.4f' %rmse_train)
print('RMSE on 10-fold CV: %.4f' %rmse_10cv)
print "\n"
Now let's try to do regression via Stochastic Gradient Descent.¶
In [28]:
# SGD is very senstitive to varying-sized feature values. So, first we need to do feature scaling.
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(x)
x_s = scaler.transform(x)
sgdreg = SGDRegressor(penalty='l2', alpha=0.1, n_iter=300)
# Compute RMSE on training data
sgdreg.fit(x_s,y)
p = sgdreg.predict(x_s)
err = p-y
total_error = np.dot(err,err)
rmse_train = np.sqrt(total_error/len(p))
# Compute RMSE using 10-fold x-validation
kf = KFold(len(x), n_folds=10)
xval_err = 0
for train,test in kf:
scaler = StandardScaler()
scaler.fit(x[train]) # Don't cheat - fit only on training data
xtrain_s = scaler.transform(x[train])
xtest_s = scaler.transform(x[test]) # apply same transformation to test data
sgdreg.fit(xtrain_s,y[train])
p = sgdreg.predict(xtest_s)
e = p-y[test]
xval_err += np.dot(e,e)
rmse_10cv = np.sqrt(xval_err/len(x))
method_name = 'Stochastic Gradient Descent Regression'
print('Method: %s' %method_name)
print('RMSE on training: %.4f' %rmse_train)
print('RMSE on 10-fold CV: %.4f' %rmse_10cv)
Using the regression implementation from Machine Learning in Action, Chapter 8:¶
In [29]:
def standRegres(xArr,yArr):
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T*xMat
if np.linalg.det(xTx) == 0.0:
print "This matrix is singular, cannot do inverse"
return
ws = xTx.I * (xMat.T*yMat)
return ws
In [30]:
w = standRegres(x,y)
In [31]:
print w
In [32]:
def ridgeRegres(xArr,yArr,lam=0.2):
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T*xMat
denom = xTx + np.eye(np.shape(xMat)[1])*lam
if np.linalg.det(denom) == 0.0:
print "This matrix is singular, cannot do inverse"
return
ws = denom.I * (xMat.T*yMat)
return ws
In [33]:
w_ridge = ridgeRegres(x,y,0.5)
print w_ridge
Now that we have the regression coefficients, we can compute the predictions:¶
In [34]:
xMat=np.mat(x)
yMat=np.mat(y)
yHat = xMat*w_ridge
In [35]:
yHat.shape
Out[35]:
In [36]:
print yHat[0:10]
In [37]:
print yMat.T[0:10]
Model evaluation and cross validation can be performed as before.¶
In [38]:
yHat.A.ravel()
Out[38]:
In [39]:
yHat.A.shape
Out[39]: