<jupyter_start><jupyter_text>--- Load diabetes data set.<jupyter_code>from sklearn import datasets
import pandas as pd

diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target
y -= y.mean()
features = "age sex bmi map tc ldl hdl tch ltg glu".split()
diabetesdf = pd.DataFrame(diabetes.data, columns=features)
diabetesdf.tail()<jupyter_output><empty_output><jupyter_text>--- OLS regression<jupyter_code>from sklearn.linear_model import LinearRegression

olsmodel = LinearRegression()
olsmodel.fit(X, y)
yfit = olsmodel.predict(X)

print("beta:", olsmodel.coef_)
print("MSE:", np.mean((y - yfit)**2))
print("R^2:", olsmodel.score(X, y))<jupyter_output>beta: [ -10.01219782 -239.81908937  519.83978679  324.39042769 -792.18416163
  476.74583782  101.04457032  177.06417623  751.27932109   67.62538639]
MSE: 2859.6903987680657
R^2: 0.5177494254132934<jupyter_text>--- Ridge regression with cross validationError to be minimized $= SSE + \lambda \sum_j \beta_j^2$$\lambda$ is refered to as a hyperparameter.We'll use cross validation to figure out what it should be.<jupyter_code>from sklearn.linear_model import RidgeCV

alphas = np.logspace(-5, -1, 50)  # 10^-5 ... 10^-1
ridgecvmodel = RidgeCV(alphas, store_cv_values=True)
ridgecvmodel.fit(X, y)
yfit = ridgecvmodel.predict(X)

print("alpha:", ridgecvmodel.alpha_)
print("beta:", ridgecvmodel.coef_)
print("MSE:", np.mean((y - yfit)**2))
print("R^2:", ridgecvmodel.score(X, y))

MSE_vs_alpha = ridgecvmodel.cv_values_.mean(axis=0)
idx = np.argmin(MSE_vs_alpha)
plt.plot(alphas, MSE_vs_alpha, 'o-')
plt.plot(alphas[idx], MSE_vs_alpha[idx], 'ro')
plt.xlabel('alpha')
plt.ylabel('RidgeCV MSE');
plt.semilogx();<jupyter_output>alpha: 0.004094915062380423
beta: [  -8.49508631 -237.25073982  521.05537865  322.45164462 -545.10254784
  280.74496815   -7.54438643  148.06050222  656.87377116   69.36047931]
MSE: 2862.0592638039234
R^2: 0.5173499463210277<jupyter_text>--- Let's do ridge regression with 10-fold cross validation by hand.<jupyter_code>from sklearn.linear_model import Ridge
from sklearn.model_selection import KFold

K = 10
kf = KFold(n_splits=K, shuffle=True, random_state=42)

alphas = np.logspace(-5, -1, 50)  # 10^-5 ... 10^-1
MSE_train_vs_alpha = []
MSE_val_vs_alpha = []
first = True
for alpha in alphas:
    # Perform 10-fold CV for each alpha and store the resulting average MSE for each of the 10 training and validation sets.
    MSE_train_cv = []
    MSE_val_cv = []
    for train_index, val_index in kf.split(X, y):
        # plot data split only on the first time
        if first:
            plt.figure()
            plt.bar(train_index, np.ones(train_index.shape), color='lightblue', label='training set')
            plt.bar(val_index, np.ones(val_index.shape), color='magenta', label='validation set')
            plt.xlabel('Patient #')
            plt.legend()
            first = False
        
        # split data into training and validation sets
        X_train = X[train_index, :]
        X_val = X[val_index, :]
        y_train = y[train_index]
        y_val = y[val_index]
        
        # new model for each split
        model = Ridge(alpha=alpha)
        
        # fit model to training set
        model.fit(X_train, y_train)
        
        # get error from training set
        y_train_fit = model.predict(X_train)
        MSE_train_cv.append(np.mean((y_train - y_train_fit)**2))
        
        # get error from test set
        y_val_fit = model.predict(X_val)
        MSE_val_cv.append(np.mean((y_val - y_val_fit)**2))
    
    MSE_train_vs_alpha.append(np.array(MSE_train_cv).mean())
    MSE_val_vs_alpha.append(np.array(MSE_val_cv).mean())

fig, ax = plt.subplots(2, 1, figsize=(8, 5))
idx = np.argmin(MSE_val_vs_alpha)
ax[0].plot(alphas, MSE_val_vs_alpha, 'o-')
ax[0].plot(alphas[idx], MSE_val_vs_alpha[idx], 'ro')
ax[0].set_xlabel('alpha')
ax[0].set_ylabel('Validation MSE');
ax[0].semilogx();

ax[1].plot(alphas, MSE_train_vs_alpha, 'o-')
ax[1].set_xlabel('alpha')
ax[1].set_ylabel('Training MSE');
ax[1].semilogx();<jupyter_output><empty_output><jupyter_text>OK, now we have selected a $\lambda$. What are the best weights for each of our features? For that, let's refit all of our data using our selected model (i.e. our selected $\lambda$).<jupyter_code>idx = np.argmin(MSE_val_vs_alpha)
best_alpha = alphas[idx]

ridgekfoldmodel = Ridge(alpha=best_alpha)
ridgekfoldmodel.fit(X, y)
yfit = ridgekfoldmodel.predict(X)

print("alpha:", best_alpha)
print("beta:", ridgekfoldmodel.coef_)
print("MSE:", np.mean((y - yfit)**2))
print("R^2:", ridgekfoldmodel.score(X, y))
plt.figure(figsize=(12, 5))
xticks = np.arange(len(features))
plt.bar(xticks-0.2, olsmodel.coef_, width=0.2, label='OLS')
plt.bar(xticks, ridgecvmodel.coef_, width=0.2, label='RidgeCV')
plt.bar(xticks+0.2, ridgekfoldmodel.coef_, width=0.2, label='Ridge K-fold CV')
plt.xticks(xticks, labels=features)
plt.ylabel('Weight')
plt.legend();<jupyter_output><empty_output><jupyter_text>--- K-fold cross validation wrapped up nice and neat for you.<jupyter_code>from sklearn.model_selection import cross_validate

model = Ridge(alpha=best_alpha)
results = cross_validate(model, X, y, cv=10, return_train_score=True)

results['test_score'], results['train_score']
model = Ridge(alpha=best_alpha)
results = cross_validate(model, X, y, cv=10, scoring='neg_mean_squared_error', return_train_score=True)

results['test_score'], results['train_score']
alphas = np.logspace(-5, -1, 50)  # 10^-5 ... 10^-1
MSE_train_vs_alpha = []
MSE_val_vs_alpha = []
for alpha in alphas:
    model = Ridge(alpha=alpha)
    results = cross_validate(model, X, y, cv=10, scoring='neg_mean_squared_error', return_train_score=True)
    MSE_train_vs_alpha.append(-results['train_score'].mean())
    MSE_val_vs_alpha.append(-results['test_score'].mean())
fig, ax = plt.subplots(2, 1, figsize=(8, 5))
idx = np.argmin(MSE_val_vs_alpha)
ax[0].plot(alphas, MSE_val_vs_alpha, 'o-')
ax[0].plot(alphas[idx], MSE_val_vs_alpha[idx], 'ro')
ax[0].set_ylabel('Validation MSE');
ax[0].semilogx();

ax[1].plot(alphas, MSE_train_vs_alpha, 'o-')
ax[1].set_xlabel('alpha')
ax[1].set_ylabel('Training MSE');
ax[1].semilogx();

fig.tight_layout();

alphas[idx]<jupyter_output><empty_output><jupyter_text>--- Compare with OLS.<jupyter_code>model = LinearRegression()
results = cross_validate(model, X, y, cv=10, scoring='neg_mean_squared_error', return_train_score=True)

print("train MSE:", -results['train_score'].mean())
print("test MSE:", -results['test_score'].mean())<jupyter_output>train MSE: 2852.5566887595296
test MSE: 3000.3810328390514<jupyter_text>--- More ways to wrap things up nicely for you.<jupyter_code>from sklearn.model_selection import GridSearchCV

params = {'alpha': np.logspace(-5, -1, 50)}
ridgemodel = Ridge()
cvmodel = GridSearchCV(ridgemodel, params, cv=10, scoring='neg_mean_squared_error')
cvmodel.fit(X, y)
cvmodel.best_params_['alpha']<jupyter_output>/anaconda3/lib/python3.7/site-packages/sklearn/model_selection/_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
  DeprecationWarning)
{'fraction_non_alphanumeric': 0.1444814198557959, 'fraction_numerical': 0.06239600665557404, 'mean_word_length': 4.89297385620915, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '234129', 'license': 'MIT', 'n_tokens_mistral': 3013, 'n_tokens_neox': 2713, 'n_words': 555, 'path': 'lecture_19_cross_validation/cross_validation.ipynb', 'repo_name': 'marcel-goldschen-ohm/CompNeuroSpring2019', 'size': '7212'}