<jupyter_start><jupyter_text>Chapter 6 Question 9Regularisation on the `College` data set<jupyter_code>import statsmodels.api as sm
import sklearn.model_selection
import sklearn.linear_model
import sklearn.metrics
import numpy as np
import matplotlib.pyplot as plt
import sklearn.decomposition
import sklearn.pipeline
import sklearn.cross_decomposition
import pandas as pd
college = sm.datasets.get_rdataset("College", "ISLR").data
college["Private"] = college["Private"] == "Yes"
college.head()<jupyter_output><empty_output><jupyter_text>(a) Split the data set into a training set and a test set<jupyter_code>X = college.drop(columns=["Apps"])
y = college["Apps"]
X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(X, y, test_size=0.33)<jupyter_output><empty_output><jupyter_text>(b) Fit a linear model using least squares, and report the test error<jupyter_code>least_squares = sklearn.linear_model.LinearRegression()
least_squares.fit(X_train,y_train)
y_pred = least_squares.predict(X_test)
least_squares_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(least_squares_mse)
# The above is vast!<jupyter_output><empty_output><jupyter_text>(c) Fit a ridge regression model, choosing $\lambda$ by CV. Report the test error<jupyter_code>ridge = sklearn.linear_model.RidgeCV(alphas= np.linspace(0.001,10,num=1000))
ridge.fit(X_train,y_train)
y_pred = ridge.predict(X_test)
ridge_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(ridge_mse)<jupyter_output>1540869.3506382525<jupyter_text>(d) Fit a lasso model, choosing $\lambda$ by CV. Report the test error, along with the number of non-zero coefficients.<jupyter_code># Use the LassoCV
lasso_model = sklearn.linear_model.LassoCV(cv=5, max_iter=1e6)
lasso_model.fit(X_train,y_train)

mses = list(map(np.mean,lasso_model.mse_path_))
alphas = lasso_model.alphas_
plt.plot(alphas,np.log(mses))
plt.ylabel("log(mse)")
plt.xlabel("alpha")
plt.show()
print(lasso_model.coef_)
print(lasso_model.intercept_)
y_pred = lasso_model.predict(X_test)
lasso_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
lasso_mse
plt.scatter(college["Accept"], college["Apps"])  # Lasso suggests a roughly one-to-one mapping between these<jupyter_output><empty_output><jupyter_text>(e) Fit a PCR model, with M chosen by cross-validation. Report the test error, along with the value of M selected.<jupyter_code># Standardise each predictor:
scaler = sklearn.preprocessing.StandardScaler()
regressor = sklearn.linear_model.LinearRegression()
pca = sklearn.decomposition.PCA()

pipe = sklearn.pipeline.Pipeline(steps=[("scaling", scaler), ("pca", pca), ("linear regression", regressor)])

p = len(X_train.columns)
params = {"pca__n_components": list(range(1,p+1))}

search = sklearn.model_selection.GridSearchCV(pipe, params, cv=5, return_train_score=True)
search.fit(X_train, y_train)
pca.fit(X_train_scaled)

fig, (ax0, ax1) = plt.subplots(nrows=2, sharex=True, figsize=(6, 6))
ax0.plot(pca.explained_variance_ratio_, linewidth=2)
ax0.set_ylabel('PCA explained variance')

ax0.axvline(search.best_estimator_.named_steps['pca'].n_components,
            linestyle=':', label='n_components chosen')
ax0.legend(prop=dict(size=12))

# For each number of components, find the best classifier results
results = pd.DataFrame(search.cv_results_)
components_col = 'param_pca__n_components'
best_clfs = results.groupby(components_col).apply(
    lambda g: g.nlargest(1, 'mean_test_score'))

best_clfs.plot(x=components_col, y='mean_test_score', yerr='std_test_score',
               legend=False, ax=ax1)
ax1.set_ylabel('Classification accuracy (val)')
ax1.set_xlabel('n_components')

plt.tight_layout()
plt.show()
# The above graph suggests PCR with 5 components.

pipe.set_params(pca__n_components=5).fit(X_train, y_train)
y_pred = pipe.predict(X_test)
pcr_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
pcr_mse<jupyter_output>/home/will/.local/anaconda3/lib/python3.6/site-packages/sklearn/preprocessing/data.py:617: DataConversionWarning: Data with input dtype bool, int64, float64 were all converted to float64 by StandardScaler.
  return self.partial_fit(X, y)
/home/will/.local/anaconda3/lib/python3.6/site-packages/sklearn/base.py:465: DataConversionWarning: Data with input dtype bool, int64, float64 were all converted to float64 by StandardScaler.
  return self.fit(X, y, **fit_params).transform(X)
/home/will/.local/anaconda3/lib/python3.6/site-packages/sklearn/pipeline.py:331: DataConversionWarning: Data with input dtype bool, int64, float64 were all converted to float64 by StandardScaler.
  Xt = transform.transform(Xt)<jupyter_text>(f) Fit a PLS model, with M chosen by cross-validation. Report the test error, along with the value of M selected.<jupyter_code># Standardise each predictor:
# scaler = sklearn.preprocessing.StandardScaler()
pls = sklearn.cross_decomposition.PLSRegression()

# pipe = sklearn.pipeline.Pipeline(steps=[("scaling", scaler), ("pca", pca), ("linear regression", regressor)])

p = len(X_train.columns)
params = {"n_components": list(range(1,p+1))}

search = sklearn.model_selection.GridSearchCV(pls, params, cv=5, return_train_score=True)
search.fit(X_train, y_train)
best_pls = search.best_estimator_
print(f"Number of components: {best_pls.n_components}")
y_pred = best_pls.predict(X_test)
pls_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
pls_mse
fig, ax = plt.subplots() 
results = pd.DataFrame(search.cv_results_)
results.head()

plt.errorbar(results.param_n_components, results.mean_test_score, yerr=results.std_test_score)
# Actually, M=6 looks best here.
pls = sklearn.cross_decomposition.PLSRegression(n_components=6)
pls.fit(X_train, y_train)
y_pred = pls.predict(X_test)
pls_mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
pls_mse<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12020636285468615, 'fraction_numerical': 0.016680997420464316, 'mean_word_length': 7.332378223495702, 'pattern_counts': {'":': 2, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '502077', 'license': 'MIT', 'n_tokens_mistral': 2196, 'n_tokens_neox': 2099, 'n_words': 416, 'path': 'Chapter6/question9.ipynb', 'repo_name': 'wllgrnt/islr-python-examples', 'size': '5815'}