🔑 Solutions: Machine Learning Basics

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

# Generate synthetic data
n_samples = 100
X = np.linspace(-2, 2, n_samples).reshape(-1, 1)
y = X**2 + np.random.normal(scale=0.5, size=X.shape)

# TODO: Statsmodels
poly = PolynomialFeatures(degree=2)
X_poly = poly.fit_transform(X)
model = sm.OLS(y, X_poly) 
fit = model.fit()
print(fit.summary())

# TODO: Numpy
coefs = np.polyfit(X.ravel(), y.ravel(), deg=2) 
y_np = np.polyval(coefs, X.ravel()) # highest power first: [b2, b1, b0]
print("\nNumpy coefficients:", coefs)

# TODO: Scikit-learn
X_poly = poly.fit_transform(X)
model = LinearRegression()
model.fit(X_poly, y)

# TODO: Plot the predictions
fig, ax = plt.subplots()
ax.scatter(X, y, label='Data')
ax.plot(X, model.predict(X_poly), label='Model', color='red')
ax.set(xlabel='X', ylabel='y', title='Polynomial Regression')
ax.legend();

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

# Generate synthetic data
np.random.seed(55)
n_samples = 100
X = np.linspace(0, 2*np.pi, n_samples).reshape(-1, 1)  # Reshape for sklearn
y = np.sin(X) + np.random.normal(scale=0.5, size=X.shape)

# 1. Visualize the data
plt.scatter(X, y, label='Data', color='blue')

# 2. Split the data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 3. Fit polynomial regression models for degrees 1 to 12
degrees = range(1, 13)
train_errors = []
test_errors = []

for degree in degrees:
    # Transform the features into polynomial features
    poly = PolynomialFeatures(degree=degree)
    X_train_poly = poly.fit_transform(X_train)
    X_test_poly = poly.transform(X_test)

    # Fit a linear regression model on the polynomial features
    model = LinearRegression()
    model.fit(X_train_poly, y_train)

    # Predict on both training and test sets
    y_train_pred = model.predict(X_train_poly)
    y_test_pred = model.predict(X_test_poly)

    # Calculate mean squared errors
    train_errors.append(mean_squared_error(y_train, y_train_pred))
    test_errors.append(mean_squared_error(y_test, y_test_pred))

# 4. Plot the training and test MSE against polynomial degree
fig, ax = plt.subplots()
ax.plot(degrees, train_errors, marker='o', label='Training MSE')
ax.plot(degrees, test_errors, marker='s', label='Test MSE')
ax.set(xlabel="Polynomial Degree", ylabel="Mean Squared Error", title="Training vs. Test Error")
ax.legend();

# 5. Repeat the split with different values for random_state
fig, ax = plt.subplots()

for random_state in [0, 1, 2, 3, 4]:
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=random_state)
    test_errors = []

    for degree in degrees:
        poly = PolynomialFeatures(degree=degree)
        model = LinearRegression().fit(poly.fit_transform(X_train), y_train)
        test_errors.append(mean_squared_error(y_test, model.predict(poly.transform(X_test))))

    best = degrees[np.argmin(test_errors)]
    print(f"random_state={random_state}: best degree = {best}, test MSE = {min(test_errors):.3f}")
    ax.plot(degrees, test_errors, marker='o', label=f"random_state={random_state}")

ax.set(xlabel="Polynomial Degree", ylabel="Test MSE", title="Test error for different splits")
ax.legend();

# Bonus: Estimate bias² and variance by simulating many datasets
rng = np.random.default_rng(0)
f_true = np.sin(X).ravel()  # the true function at each x
sigma = 0.5                 # the noise level used to generate the data
n_sims = 500

bias2, variance = [], []

for degree in degrees:
    poly = PolynomialFeatures(degree=degree)
    X_poly = poly.fit_transform(X)

    # Fit one model per simulated dataset (same x, new noise) and store its predictions
    preds = np.empty((n_sims, len(X)))
    for s in range(n_sims):
        y_sim = f_true + rng.normal(0, sigma, size=len(X))
        preds[s] = LinearRegression().fit(X_poly, y_sim).predict(X_poly)

    bias2.append(np.mean((preds.mean(axis=0) - f_true) ** 2))
    variance.append(np.mean(preds.var(axis=0)))

total = np.array(bias2) + np.array(variance) + sigma**2

fig, ax = plt.subplots()
ax.plot(degrees, bias2, marker='o', label='Bias²')
ax.plot(degrees, variance, marker='o', label='Variance')
ax.axhline(sigma**2, color='grey', ls=':', label='Irreducible error $\\sigma^2$')
ax.plot(degrees, total, marker='o', color='k', label='Expected test MSE')
ax.set(xlabel="Polynomial Degree", ylabel="Error", yscale="log", title="Bias-variance decomposition")
ax.legend();

from sklearn.datasets import fetch_california_housing

data = fetch_california_housing(as_frame=True)
df = data.frame

df.head()

df.hist(figsize=(12, 10), bins=30, edgecolor="black");

# TODO: Prepare data
X, y = fetch_california_housing(return_X_y=True)

from sklearn.model_selection import KFold
from sklearn.linear_model import LinearRegression 
from sklearn.model_selection import cross_val_score 

# TODO: Implement CV
k_fold = KFold(n_splits=5, shuffle=True, random_state=42)  # shuffle, since the rows are ordered
model = LinearRegression()

neg_mse = cross_val_score(model, X, y, cv=k_fold, scoring='neg_mean_squared_error')
mse = -neg_mse

print(f"Average MSE:       {mse.mean()}")
print(f"MSE for each fold: {mse}")

# 3. Visualize the folds with and without shuffling and compare the fold MSEs
import numpy as np
import matplotlib.pyplot as plt

kf = KFold(n_splits=5)
kf_shuffled = KFold(n_splits=5, shuffle=True, random_state=42)

# Mark the validation rows of each fold with 1 (training rows stay 0)
folds = np.zeros((5, len(y)))
for fold, (train_idx, val_idx) in enumerate(kf.split(X)):
    folds[fold, val_idx] = 1

folds_shuffled = np.zeros((5, len(y)))
for fold, (train_idx, val_idx) in enumerate(kf_shuffled.split(X)):
    folds_shuffled[fold, val_idx] = 1

# Plot the folds
fig, ax = plt.subplots(2, 1, figsize=(10, 4), sharex=True)

ax[0].imshow(folds, aspect="auto", cmap="Blues", interpolation="nearest")
ax[0].set(title="Without shuffling", ylabel="Fold")

ax[1].imshow(folds_shuffled, aspect="auto", cmap="Blues", interpolation="nearest")
ax[1].set(title="With shuffling", xlabel="Row in the dataset", ylabel="Fold")

plt.tight_layout()

# Compare the MSE of the individual folds
mse = -cross_val_score(LinearRegression(), X, y, cv=kf, scoring='neg_mean_squared_error')
mse_shuffled = -cross_val_score(LinearRegression(), X, y, cv=kf_shuffled, scoring='neg_mean_squared_error')

print(f"Without shuffling: {np.round(mse, 3)}   mean: {mse.mean():.3f}")
print(f"With shuffling:    {np.round(mse_shuffled, 3)}   mean: {mse_shuffled.mean():.3f}")

# TODO: Implement LOOCV
import numpy as np
from matplotlib import pyplot as plt
from sklearn.model_selection import LeaveOneOut

# 1. Random subset of 5000 observations
rng = np.random.default_rng(42)
idx = rng.choice(len(y), size=5000, replace=False)
X_sub, y_sub = X[idx], y[idx]

model = LinearRegression()

# LOOCV
mse = -cross_val_score(model, X_sub, y_sub, cv=LeaveOneOut(), scoring='neg_mean_squared_error')

# 5-fold CV on the same subset for comparison
k_fold = KFold(n_splits=5, shuffle=True, random_state=42)
mse_kfold = -cross_val_score(model, X_sub, y_sub, cv=k_fold, scoring='neg_mean_squared_error')

# 2. and 4. Get MSE metrics
print(f"Average MSE (LOOCV):  {mse.mean():.3f}")
print(f"Average MSE (5-fold): {mse_kfold.mean():.3f}")
print(f"Minimum MSE:          {mse.min():.2e}")
print(f"Maximum MSE:          {mse.max():.1f}")
print(f"Median MSE:           {np.median(mse):.3f}")

# 3. Plot MSE
plt.hist(mse, bins=50, range=(0, 6));

# Which observation causes the largest error?
feature_names = fetch_california_housing().feature_names
worst = np.argmax(mse)
print(dict(zip(feature_names, X_sub[worst].round(1))))