<jupyter_start><jupyter_text>Description and Details: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:* Polynomial regression with a small penalty:    * Use the pre-built implementation of regression  to run polynomial regression    * Use matplotlib to visualize polynomial regressions on four datasets    * Use a pre-built implementation of regression to run polynomial regression, this time with L2_small_penalty    * Observe overfitting: Use matplotlib to visualize polynomial regressions under L2_small_penalty* Polynomial regression with an appropriate penalty:    * Use matplotlib to visualize polynomial regressions under L2_large_penalty    * Choose best L2 penalty by building hands-on cross-validation function    * Another option to choose best L2 penalty: using sklearn.ridgeCV    * Assess the final fit using test data.<jupyter_code>import pandas as pd
import numpy as np
from sklearn import linear_model
import matplotlib.pyplot as plt
%matplotlib inline
import glob<jupyter_output><empty_output><jupyter_text>Load Data<jupyter_code>dtype_dict = {'bathrooms':float, 'waterfront':int, 'sqft_above':int, 'sqft_living15':float, 'grade':int, 'yr_renovated':int, 'price':float, 'bedrooms':float, 'zipcode':str, 'long':float, 'sqft_lot15':float, 'sqft_living':float, 'floors':float, 'condition':int, 'lat':float, 'date':str, 'sqft_basement':int, 'yr_built':int, 'id':str, 'sqft_lot':int, 'view':int}

sales = pd.read_csv('kc_house_data.csv', dtype=dtype_dict)
sales = sales.sort_values(['sqft_living','price'])<jupyter_output><empty_output><jupyter_text>Polynomial Regression with L2 Regularization<jupyter_code>def polynomial(feature, degree):
    poly_dataframe = pd.DataFrame()
    poly_dataframe['power_1'] = feature

    if degree > 1:
        for power in range(2, degree+1): 
            name = 'power_' + str(power)
            tmp = feature.apply(lambda x: x**power)
            poly_dataframe[name] = tmp

    return poly_dataframe
poly15_data = polynomial(sales['sqft_living'], 15) 
poly15_data
l2_small_penalty = 1.5e-5
model = linear_model.Ridge(alpha=l2_small_penalty, normalize=True)
model.fit(poly15_data, sales['price'])
print(model.coef_)<jupyter_output>[ 1.24873306e+02 -4.77376011e-02  3.01446238e-05 -2.44419942e-09
 -1.94153675e-13  8.54085686e-18  1.51142121e-21  8.27979094e-26
  6.52603100e-31 -3.27895017e-34 -3.87962315e-38 -2.72437650e-42
 -1.07790800e-46  3.78242694e-51  1.39790296e-54]<jupyter_text>Observe Overfitting<jupyter_code>set_1 = pd.read_csv('wk3_kc_house_set_1_data.csv')#, dtype=dtype_dict)
set_2 = pd.read_csv('wk3_kc_house_set_2_data.csv', dtype=dtype_dict)
set_3 = pd.read_csv('wk3_kc_house_set_3_data.csv', dtype=dtype_dict)
set_4 = pd.read_csv('wk3_kc_house_set_4_data.csv', dtype=dtype_dict)
print(set_1.shape)
print(set_2.shape)
print(set_3.shape)
print(set_4.shape)
l2_small_penalty = 1.5e-5
model = linear_model.Ridge(alpha=l2_small_penalty, normalize=True)<jupyter_output><empty_output><jupyter_text>for set_1:<jupyter_code>poly1_data = polynomial(set_1['sqft_living'], 15)
features = poly1_data.columns
poly1_data['price'] = set_1['price']

print(poly1_data.min())
print(poly1_data.shape)

model.fit(poly1_data[features], poly1_data['price'])
print(model.intercept_)

#plt.plot(poly1_data['power_1'], poly1_data['price'],'.')
#plt.plot(poly1_data['power_1'], model.predict(poly1_data[features]), ',', color = 'red', alpha = 1)<jupyter_output>power_1     4.300000e+02
power_2     1.849000e+05
power_3     7.950700e+07
power_4     3.418801e+10
power_5     1.470084e+13
power_6     6.321363e+15
power_7     2.718186e+18
power_8     1.168820e+21
power_9     5.025926e+23
power_10    2.161148e+26
power_11    9.292937e+28
power_12    3.995963e+31
power_13    1.718264e+34
power_14    7.388536e+36
power_15    3.177070e+39
price       8.000000e+04
dtype: float64
(5404, 16)
288998.10746196454<jupyter_text>Loop over set_1-4 :<jupyter_code>print('all csv files in data directory:', glob.glob('wk3_kc_house_set_*_data.csv'))
for filename in glob.glob('wk3_kc_house_set_*_data.csv'):      # loop over datasets
    raw_data = pd.read_csv(filename)
    poly_data = polynomial(raw_data['sqft_living'], 15)
    features = poly_data.columns
    poly_data['price'] = raw_data['price']
    print(poly_data.shape)
    print(poly_data.iloc[:2,:2])
    model = linear_model.Ridge(alpha=l2_small_penalty, normalize=True)
    model.fit(poly_data[features], poly_data['price'])
    print('coefficient of power_1 is:', model.intercept_)
    plt.plot(poly_data['power_1'],poly_data['price'],'.',
        poly_data['power_1'], model.predict(poly_data[features]),'-')
    print(filename)
    plt.show()<jupyter_output>(5404, 16)
   power_1  power_2
0      430   184900
1      460   211600
coefficient of power_1 is: 288998.10746196454
wk3_kc_house_set_1_data.csv<jupyter_text>Ridge Regresion 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.<jupyter_code>l2_large_penalty = 1.23e2
for filename in glob.glob('wk3_kc_house_set_*_data.csv'):      # loop over datasets
    raw_data = pd.read_csv(filename)
    poly_data = polynomial(raw_data['sqft_living'], 15)
    features = poly_data.columns
    poly_data['price'] = raw_data['price']
    print(poly_data.shape)
    print(poly_data.iloc[:2,:2])
    model = linear_model.Ridge(alpha=l2_large_penalty, normalize=True)
    model.fit(poly_data[features], poly_data['price'])
    print('coefficient of power_1 is:', model.intercept_)
    plt.plot(poly_data['power_1'],poly_data['price'],'.',
        poly_data['power_1'], model.predict(poly_data[features]),'-')
    print(filename)
    plt.show()<jupyter_output>(5404, 16)
   power_1  power_2
0      430   184900
1      460   211600
coefficient of power_1 is: 539103.0310577955
wk3_kc_house_set_1_data.csv<jupyter_text>Selecting an L2 penalty via hands-on cross-validation<jupyter_code>train_valid_shuffled = pd.read_csv('wk3_kc_house_train_valid_shuffled.csv', dtype=dtype_dict)
test = pd.read_csv('wk3_kc_house_test_data.csv', dtype=dtype_dict)
len(train_valid_shuffled)
n = len(train_valid_shuffled)
k = 10  # number of folds

for i in range(10):
    start = round((n*i)/k)
    end = round((n*(i+1))/k)
    print(i, (start, end), len(range(start, end)))
poly_data = polynomial(train_valid_shuffled['sqft_living'], 15)
features = poly_data.columns                       # get the name of the features
poly_data['price'] = train_valid_shuffled['price'] # add price to the data since 
poly_data
def k_fold_cross_validation(k, l2_penalty, data, output):
    poly_data = polynomial(data['sqft_living'], 15)
    features = poly_data.columns                       
    poly_data['price'] = data['price']     
    RSS_final = np.zeros((k, len(l2_penalty)))
    n = len(data)
    for i in range(k):
        start = round((n*i)/k)
        end = round((n*(i+1))/k)
        validation_set = poly_data[start:end]
        train_set = poly_data[0:start].append(poly_data[end:n+1])
        for j in range(len(l2_penalty)):
            model = linear_model.Ridge(alpha = l2_penalty[j], normalize=True)
            model.fit(train_set[features], train_set[output])
            predictions = model.predict(validation_set[features])
            RSS = np.sum((validation_set[output] - predictions)**2)
            RSS_final[i][j] = RSS
    return RSS_final.mean(axis = 0)         
a = np.logspace(-5, 1, num = 20)
cv_error = k_fold_cross_validation(k = 10, l2_penalty = a, data = train_valid_shuffled, output = 'price')
result = pd.DataFrame({'a': a, 'error': cv_error})
result.set_index('a', inplace = True)   # 把a设为index

bestlam = result.index[result['error'] == result['error'].min()].tolist()
bestlam
plt.plot(a,cv_error)
plt.xlabel('$\ell_2$ penalty')
plt.ylabel('cross validation error')
plt.xscale('log')
plt.yscale('log')
print ('The penalties are:','\n', a, '\n')
print ('The average validation error are:', '\n', cv_error)<jupyter_output>The penalties are: 
 [1.00000000e-05 2.06913808e-05 4.28133240e-05 8.85866790e-05
 1.83298071e-04 3.79269019e-04 7.84759970e-04 1.62377674e-03
 3.35981829e-03 6.95192796e-03 1.43844989e-02 2.97635144e-02
 6.15848211e-02 1.27427499e-01 2.63665090e-01 5.45559478e-01
 1.12883789e+00 2.33572147e+00 4.83293024e+00 1.00000000e+01] 

The average validation error are: 
 [2.12721097e+14 4.68515068e+14 6.94153920e+14 7.63463191e+14
 6.91513195e+14 5.57978406e+14 4.20863037e+14 2.99383684e+14
 2.04776499e+14 1.48592223e+14 1.26044160e+14 1.20453079e+14
 1.19850480e+14 1.20133800e+14 1.21125799e+14 1.24331935e+14
 1.34245911e+14 1.58319145e+14 1.93852057e+14 2.23690752e+14]<jupyter_text>Another option to choose L2 penalty : using sklearn.RidgeCV<jupyter_code>from sklearn.linear_model import RidgeCV
poly_data = polynomial(train_valid_shuffled['sqft_living'], 15)  # 产生1-15次项
features = poly_data.columns                       
poly_data['price'] = train_valid_shuffled['price']
poly_data
a = np.logspace(-5,1, num = 20)
clf = RidgeCV(alphas = a, cv = 10, normalize = True).fit(poly_data[features], 
                                                         poly_data['price'])
print('The coefficients are:', '\n', clf.coef_, '\n')
print('The best lambda is:', '\n', clf.alpha_, '\n')
bestlam_auto = clf.alpha_
if bestlam_auto == bestlam:
    print('Pass!')<jupyter_output>Pass!<jupyter_text>Retrain the model on full training set using selected lambda<jupyter_code>final_data = polynomial(train_valid_shuffled['sqft_living'], 15)
features = final_data.columns                       
final_data['price'] = train_valid_shuffled['price']  
final_data
final_model = linear_model.Ridge(alpha = bestlam, normalize=True)
final_model.fit(final_data[features], final_data['price'])
test_data = pd.read_csv('wk3_kc_house_test_data.csv')
test_data.head()
final_test_data = polynomial(test_data['sqft_living'], 15)
features = final_test_data.columns                       
final_test_data['price'] = test_data['price']  
final_test_data
final_prediction = final_model.predict(final_test_data[features])
RSS_final = np.sum((final_test_data['price'] - final_prediction)**2)
RSS_final<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12337782721542455, 'fraction_numerical': 0.0967741935483871, 'mean_word_length': 4.8699673558215455, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '609971', 'license': 'MIT', 'n_tokens_mistral': 4564, 'n_tokens_neox': 3929, 'n_words': 797, 'path': 'week4/HW4_1.ipynb', 'repo_name': 'jessie-sun927/UW-Machine-Learning-Regression', 'size': '10788'}