Exercises: Regression Algorithms¶

These exercises build sequentially — each one adds a layer of complexity on top of the previous. Work through them in order. The warm-up exercises establish your baseline. The main exercises build toward a complete modelling workflow. The stretch exercises require genuine debugging and thought.

Exercise 1 — Warm-Up: Fit a Linear Regression Baseline¶

Dataset: fetch_california_housing() from scikit-learn (or any housing CSV you have)

Tasks:

Load the dataset and inspect it: how many features, how many rows, what is the target's range?
Split into 80% train / 20% test using random_state=42
Fit a LinearRegression wrapped in a Pipeline with StandardScaler
Report: MAE, RMSE, and R² on the test set
Print the coefficients as a sorted pd.Series and identify the top 2 most influential features (by absolute coefficient value)

Expected output (approximate):

MAE:  0.533
RMSE: 0.745
R²:   0.576
Top positive coefficient: MedInc (0.855)
Top negative coefficient: Latitude (-0.899)

Show answer

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

print(f"Dataset: {X.shape[0]} rows, {X.shape[1]} features")
print(f"Target range: {y.min():.2f} to {y.max():.2f} ($100k units)")

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("model", LinearRegression())
])
pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)

mae  = mean_absolute_error(y_test, y_pred)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2   = r2_score(y_test, y_pred)

print(f"\nMAE:  {mae:.3f}")
print(f"RMSE: {rmse:.3f}")
print(f"R²:   {r2:.3f}")

coefs = pd.Series(
    pipe.named_steps["model"].coef_,
    index=X.columns
).sort_values()

print("\nCoefficients (sorted):")
print(coefs.round(3))

Exercise 2 — Warm-Up: Build and Read a Residual Plot¶

Using the model from Exercise 1:

Compute the residuals (y_test - y_pred)
Create a 2-panel plot:
Left panel: residuals vs. predicted values, with a red horizontal line at 0
Right panel: histogram of residuals
Answer these questions in a comment in your code:
Is there a visible pattern in the residuals vs. fitted plot?
Are the residuals approximately symmetric?
Does the spread of residuals seem consistent across the range of predictions, or does it fan out?

Show answer

import matplotlib.pyplot as plt

residuals = y_test.values - y_pred

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].scatter(y_pred, residuals, alpha=0.3, s=10, color="steelblue")
axes[0].axhline(0, color="red", linewidth=1.5)
axes[0].set_xlabel("Predicted Value (units of $100k)")
axes[0].set_ylabel("Residual")
axes[0].set_title("Residuals vs. Fitted — Linear Regression")

axes[1].hist(residuals, bins=50, edgecolor="black", color="steelblue")
axes[1].axvline(0, color="red", linewidth=1.5)
axes[1].set_xlabel("Residual")
axes[1].set_title("Residual Distribution")

plt.tight_layout()
plt.show()

# Observations for California Housing + Linear Regression:
# - The residuals vs. fitted plot shows a fan shape: spread increases
#   for higher predicted values — heteroscedasticity is present.
# - The distribution is approximately bell-shaped but has a right tail
#   (the model underestimates some high-value houses).
# - There is a clear "ceiling" effect around predicted values of 5.0
#   ($500k) — the dataset caps house values at $500k, creating
#   a cluster of residuals that linear regression cannot handle.

Exercise 3 — Main: Compare Five Models¶

Dataset: fetch_california_housing()

Tasks:

Split data (80/20, random_state=42)
Fit and evaluate all five models below on the test set:
LinearRegression (with StandardScaler)
RidgeCV (with StandardScaler, search alphas with np.logspace(-3, 4, 50))
LassoCV (with StandardScaler)
RandomForestRegressor (n_estimators=200, random_state=42)
GradientBoostingRegressor (n_estimators=200, learning_rate=0.1, max_depth=4, random_state=42)
Print a comparison table: model name, MAE, RMSE, R²
Identify the best model by RMSE
Print the feature importances for the best tree-based model (sorted descending)

Expected output format:

Model                    MAE      RMSE     R²
LinearRegression         0.533    0.745    0.576
RidgeCV                  0.533    0.745    0.576
LassoCV                  0.533    0.745    0.576
RandomForest             0.328    0.502    0.805
GradientBoosting         0.299    0.462    0.833

Show answer

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

models = {
    "LinearRegression": Pipeline([
        ("scaler", StandardScaler()),
        ("model", LinearRegression())
    ]),
    "RidgeCV": Pipeline([
        ("scaler", StandardScaler()),
        ("model", RidgeCV(alphas=np.logspace(-3, 4, 50)))
    ]),
    "LassoCV": Pipeline([
        ("scaler", StandardScaler()),
        ("model", LassoCV(cv=5, max_iter=10_000, random_state=42))
    ]),
    "RandomForest": RandomForestRegressor(
        n_estimators=200, n_jobs=-1, random_state=42
    ),
    "GradientBoosting": GradientBoostingRegressor(
        n_estimators=200, learning_rate=0.1, max_depth=4, random_state=42
    ),
}

results = []
fitted_models = {}
for name, model in models.items():
    model.fit(X_train, y_train)
    fitted_models[name] = model
    preds = model.predict(X_test)
    results.append({
        "Model": name,
        "MAE":  round(mean_absolute_error(y_test, preds), 3),
        "RMSE": round(np.sqrt(mean_squared_error(y_test, preds)), 3),
        "R²":   round(r2_score(y_test, preds), 3),
    })

comparison = pd.DataFrame(results).set_index("Model")
print(comparison.to_string())

best_model_name = comparison["RMSE"].idxmin()
print(f"\nBest model by RMSE: {best_model_name}")

# Feature importance for Gradient Boosting (best tree-based model)
gb_model = fitted_models["GradientBoosting"]
feat_importance = pd.Series(
    gb_model.feature_importances_, index=X.columns
).sort_values(ascending=False)
print("\nGradient Boosting Feature Importances:")
print(feat_importance.round(3))

Exercise 4 — Main: Regularisation Strength Analysis¶

Goal: Understand how alpha controls model complexity in Ridge regression.

Tasks:

Using fetch_california_housing(), create a Ridge pipeline (with StandardScaler)
Train it for these alpha values: [0.001, 0.01, 0.1, 1, 10, 100, 1000]
For each alpha, record: train R², test R², and the number of non-zero coefficients (Ridge never zeros out, so count coefficients with |coef| > 0.001)
Print a comparison table
Plot train R² and test R² vs. log(alpha) on the same chart
Answer: at which alpha does overfitting start? How can you tell?

Show answer

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

alphas = [0.001, 0.01, 0.1, 1, 10, 100, 1000]
rows = []

for alpha in alphas:
    pipe = Pipeline([
        ("scaler", StandardScaler()),
        ("model", Ridge(alpha=alpha))
    ])
    pipe.fit(X_train, y_train)

    train_r2 = r2_score(y_train, pipe.predict(X_train))
    test_r2  = r2_score(y_test, pipe.predict(X_test))
    coef_values = pipe.named_steps["model"].coef_
    large_coefs = np.sum(np.abs(coef_values) > 0.001)

    rows.append({
        "Alpha": alpha,
        "Train R²": round(train_r2, 4),
        "Test R²":  round(test_r2, 4),
        "Large Coefficients": large_coefs
    })

df_alpha = pd.DataFrame(rows).set_index("Alpha")
print(df_alpha.to_string())

# Plot
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(alphas, df_alpha["Train R²"], marker="o", label="Train R²", color="steelblue")
ax.semilogx(alphas, df_alpha["Test R²"],  marker="s", label="Test R²", color="orange")
ax.set_xlabel("Alpha (log scale)")
ax.set_ylabel("R²")
ax.set_title("Ridge: Train vs. Test R² across Alpha Values")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Answer: With Ridge on California Housing (well-conditioned), the test R² is
# quite stable across a wide range of alpha. The gap between train and test R²
# is small throughout — this dataset does not overfit heavily with linear models.
# On high-dimensional datasets, you would see train R² staying high while
# test R² drops at low alpha (overfitting) and drops again at very high alpha
# (underfitting — model is too constrained).

Exercise 5 — Main: Diagnose a Failing Model¶

Goal: Use residual plots to identify specific weaknesses in a model.

Setup: Fit a LinearRegression on fetch_california_housing().

Tasks:

Fit the model and compute residuals
Create the following four diagnostic plots in a 2x2 grid:
Residuals vs. Fitted values
Actual vs. Predicted (with the perfect-prediction diagonal)
Residual histogram
Absolute residuals vs. Fitted (scale-location plot)
Investigate: do high-predicted-value homes have larger or smaller residuals on average?
Check: is the mean residual close to zero?
Write three bullet points describing specific weaknesses the plots reveal

Show answer

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

pipe = Pipeline([("scaler", StandardScaler()), ("model", LinearRegression())])
pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)
residuals = y_test.values - y_pred

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Residuals vs. Fitted
axes[0, 0].scatter(y_pred, residuals, alpha=0.2, s=8, color="steelblue")
axes[0, 0].axhline(0, color="red", linewidth=1.5)
axes[0, 0].set_xlabel("Predicted")
axes[0, 0].set_ylabel("Residual")
axes[0, 0].set_title("Residuals vs. Fitted")

# 2. Actual vs. Predicted
min_val, max_val = y_test.min(), y_test.max()
axes[0, 1].scatter(y_test, y_pred, alpha=0.2, s=8, color="steelblue")
axes[0, 1].plot([min_val, max_val], [min_val, max_val],
                color="red", linestyle="--", linewidth=1.5)
axes[0, 1].set_xlabel("Actual")
axes[0, 1].set_ylabel("Predicted")
axes[0, 1].set_title("Actual vs. Predicted")

# 3. Residual histogram
axes[1, 0].hist(residuals, bins=60, edgecolor="black", color="steelblue")
axes[1, 0].axvline(0, color="red", linewidth=1.5)
axes[1, 0].set_xlabel("Residual")
axes[1, 0].set_title("Residual Distribution")

# 4. Scale-location
axes[1, 1].scatter(y_pred, np.abs(residuals), alpha=0.2, s=8, color="steelblue")
axes[1, 1].set_xlabel("Predicted")
axes[1, 1].set_ylabel("|Residual|")
axes[1, 1].set_title("Scale-Location (Heteroscedasticity Check)")

plt.suptitle("Linear Regression Diagnostics — California Housing", y=1.02)
plt.tight_layout()
plt.show()

# Analysis
high_mask = y_pred > 3.0
print(f"Mean residual (all):        {np.mean(residuals):.4f}")
print(f"Mean residual (high-pred):  {np.mean(residuals[high_mask]):.4f}")
print(f"Residual std (all):         {np.std(residuals):.4f}")
print(f"Residual std (high-pred):   {np.std(residuals[high_mask]):.4f}")

# Expected output:
# Mean residual (all):        0.0000
# Mean residual (high-pred): -0.32   (model underestimates expensive houses)
# Residual std (all):         0.745
# Residual std (high-pred):   0.94   (larger errors for high-value homes)

# Weaknesses revealed by the plots:
# 1. Heteroscedasticity: the scale-location plot fans outward at higher
#    predicted values — the model is less accurate for expensive houses.
# 2. Systematic underprediction at the top: actual vs. predicted shows
#    many points BELOW the diagonal for high actual values — the model
#    consistently undershoots expensive homes.
# 3. Non-normality: the residual histogram has a right tail, indicating
#    some large positive residuals (true price much higher than predicted).

Exercise 6 — Stretch: Tune a Random Forest with Cross-Validation¶

Goal: Find optimal hyperparameters for a Random Forest regressor using grid search.

Dataset: fetch_california_housing()

Tasks:

Split data (80/20, random_state=42)
Define a hyperparameter grid:
n_estimators: [100, 200]
max_depth: [None, 10, 20]
min_samples_leaf: [1, 4, 8]
max_features: ["sqrt", 0.5]
Use RandomizedSearchCV with n_iter=15, cv=5, scoring="neg_root_mean_squared_error", random_state=42
Print the best parameters and the corresponding cross-validated RMSE
Evaluate the best estimator on the test set
Compare to the default Random Forest (no tuning) — how much did tuning help?

Show answer

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, RandomizedSearchCV
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Default model baseline
default_rf = RandomForestRegressor(n_estimators=100, n_jobs=-1, random_state=42)
default_rf.fit(X_train, y_train)
default_preds = default_rf.predict(X_test)
default_rmse = np.sqrt(mean_squared_error(y_test, default_preds))
print(f"Default RF RMSE: {default_rmse:.4f}")

# Hyperparameter search
param_grid = {
    "n_estimators":    [100, 200],
    "max_depth":       [None, 10, 20],
    "min_samples_leaf": [1, 4, 8],
    "max_features":    ["sqrt", 0.5],
}

rf_search = RandomizedSearchCV(
    estimator=RandomForestRegressor(n_jobs=-1, random_state=42),
    param_distributions=param_grid,
    n_iter=15,
    cv=5,
    scoring="neg_root_mean_squared_error",
    random_state=42,
    n_jobs=-1,
    verbose=1
)
rf_search.fit(X_train, y_train)

print(f"\nBest parameters: {rf_search.best_params_}")
best_cv_rmse = -rf_search.best_score_
print(f"Best CV RMSE:    {best_cv_rmse:.4f}")

# Evaluate best estimator on test set
best_preds = rf_search.best_estimator_.predict(X_test)
best_rmse = np.sqrt(mean_squared_error(y_test, best_preds))
best_r2   = r2_score(y_test, best_preds)
best_mae  = mean_absolute_error(y_test, best_preds)

print(f"\nBest Model Test Performance:")
print(f"  MAE:  {best_mae:.4f}")
print(f"  RMSE: {best_rmse:.4f}")
print(f"  R²:   {best_r2:.4f}")
print(f"\nImprovement from tuning: {(default_rmse - best_rmse) / default_rmse * 100:.1f}% RMSE reduction")

Exercise 7 — Stretch: Apply a Log Transform to Handle a Skewed Target¶

Context: Many real-world regression targets (prices, incomes, durations) are right-skewed. Predicting in log-space and back-transforming can substantially improve model accuracy and fix heteroscedasticity.

Dataset: Use fetch_california_housing(), but treat the target as if prices are skewed (apply np.log1p as a demonstration of the technique).

Tasks:

Inspect the distribution of the original target (y.hist())
Apply y_log = np.log1p(y) and plot the transformed distribution — is it more symmetric?
Fit a RandomForestRegressor on both:
Original target
Log-transformed target (back-transform predictions with np.expm1)
Compare MAE and RMSE for both approaches
Plot residuals for both models side by side
Explain in a comment: when would the log-transform approach help most?

Show answer

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

housing = fetch_california_housing(as_frame=True)
X, y = housing.data, housing.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Inspect original target distribution
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].hist(y, bins=50, edgecolor="black")
axes[0].set_title("Original Target Distribution")
axes[0].set_xlabel("Median House Value ($100k)")

y_log = np.log1p(y)
axes[1].hist(y_log, bins=50, edgecolor="black")
axes[1].set_title("Log-Transformed Target Distribution")
axes[1].set_xlabel("log(1 + Median House Value)")

plt.tight_layout()
plt.show()

# Model 1: Original target
rf_original = RandomForestRegressor(n_estimators=200, n_jobs=-1, random_state=42)
rf_original.fit(X_train, y_train)
preds_original = rf_original.predict(X_test)

# Model 2: Log-transformed target
y_train_log = np.log1p(y_train)
rf_log = RandomForestRegressor(n_estimators=200, n_jobs=-1, random_state=42)
rf_log.fit(X_train, y_train_log)
preds_log = np.expm1(rf_log.predict(X_test))  # back-transform

# Compare
mae_orig  = mean_absolute_error(y_test, preds_original)
rmse_orig = np.sqrt(mean_squared_error(y_test, preds_original))
r2_orig   = r2_score(y_test, preds_original)

mae_log   = mean_absolute_error(y_test, preds_log)
rmse_log  = np.sqrt(mean_squared_error(y_test, preds_log))
r2_log    = r2_score(y_test, preds_log)

print(f"{'':30} {'MAE':>8} {'RMSE':>8} {'R²':>8}")
print(f"{'Random Forest (original)':30} {mae_orig:>8.3f} {rmse_orig:>8.3f} {r2_orig:>8.3f}")
print(f"{'Random Forest (log target)':30} {mae_log:>8.3f}  {rmse_log:>8.3f} {r2_log:>8.3f}")

# Residual comparison
res_orig = y_test.values - preds_original
res_log  = y_test.values - preds_log

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].scatter(preds_original, res_orig, alpha=0.2, s=8)
axes[0].axhline(0, color="red", linewidth=1)
axes[0].set_title("Residuals — Original Target")

axes[1].scatter(preds_log, res_log, alpha=0.2, s=8)
axes[1].axhline(0, color="red", linewidth=1)
axes[1].set_title("Residuals — Log-Transformed Target")

plt.tight_layout()
plt.show()

# When to use log-transform:
# - Target is right-skewed (long right tail of high values)
# - Residuals fan out for high predicted values (heteroscedasticity)
# - Errors as a percentage of the true value are more meaningful
#   than absolute errors (e.g. being off by $50k on a $100k house
#   is worse than being off by $50k on a $2M house)
# - Log-transform compresses the scale and makes the model treat
#   proportional errors equally, which usually reduces RMSE on
#   right-skewed targets.

05-regression-metrics | ../Day-02-Part-1-Classification-Algorithms/00-agenda