Exercises: Model Evaluation

These exercises require you to apply the full evaluation workflow — not just run a function and read a number. The stretch exercises require you to interpret results and make decisions, which is the part that matters in practice.

Setup

Run this block first. All exercises use these datasets and imports.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.datasets import (
    fetch_california_housing,
    load_breast_cancer,
    make_classification
)
from sklearn.model_selection import (
    train_test_split,
    cross_val_score,
    cross_validate,
    StratifiedKFold,
    GridSearchCV,
    RandomizedSearchCV,
    learning_curve,
    TimeSeriesSplit
)
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor, GradientBoostingClassifier
from sklearn.dummy import DummyClassifier, DummyRegressor
from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    r2_score,
    confusion_matrix,
    classification_report,
    roc_auc_score,
    average_precision_score,
    roc_curve,
    precision_recall_curve,
    f1_score,
    precision_score,
    recall_score,
    ConfusionMatrixDisplay
)
from scipy.stats import randint, uniform
import warnings
warnings.filterwarnings("ignore")

np.random.seed(42)

# --- Regression dataset ---
housing = fetch_california_housing()
X_reg, y_reg = housing.data, housing.target
X_reg_train, X_reg_test, y_reg_train, y_reg_test = train_test_split(
    X_reg, y_reg, test_size=0.2, random_state=42
)

# --- Balanced classification dataset ---
cancer = load_breast_cancer()
X_clf, y_clf = cancer.data, cancer.target
X_clf_train, X_clf_test, y_clf_train, y_clf_test = train_test_split(
    X_clf, y_clf, test_size=0.2, random_state=42, stratify=y_clf
)

# --- Imbalanced classification dataset: 5% positives ---
X_imb, y_imb = make_classification(
    n_samples=3000,
    n_features=20,
    weights=[0.95, 0.05],
    flip_y=0.01,
    random_state=42
)
X_imb_train, X_imb_test, y_imb_train, y_imb_test = train_test_split(
    X_imb, y_imb, test_size=0.2, random_state=42, stratify=y_imb
)

print("Setup complete.")
print(f"  Regression:             {X_reg_train.shape[0]} train / {X_reg_test.shape[0]} test")
print(f"  Classification:         {X_clf_train.shape[0]} train / {X_clf_test.shape[0]} test")
print(f"  Imbalanced (5%+):       {X_imb_train.shape[0]} train / {X_imb_test.shape[0]} test")
print(f"  Imbalanced positive %:  {y_imb_train.mean():.1%}")

---

Warm-Up Exercises

These reinforce the core mechanics. If you find any of these difficult, re-read the relevant note before moving to the main exercises.

WU-1: Baseline Comparison

Train a DummyRegressor and a RandomForestRegressor on the California housing data. Print a comparison table with MAE, RMSE, and R² for both models.

Expected output format:

Model                   MAE       RMSE         R²
--------------------------------------------------
Baseline (mean)       0.534      0.736      0.000
Random Forest         0.329      0.505      0.805

Show answer

# Baseline
baseline_reg = DummyRegressor(strategy="mean")
baseline_reg.fit(X_reg_train, y_reg_train)
y_base_pred = baseline_reg.predict(X_reg_test)

# Random Forest
rf_reg = Pipeline([
    ("scaler", StandardScaler()),
    ("rf", RandomForestRegressor(n_estimators=100, random_state=42))
])
rf_reg.fit(X_reg_train, y_reg_train)
y_rf_pred = rf_reg.predict(X_reg_test)

print(f"{'Model':<25} {'MAE':>8} {'RMSE':>8} {'R²':>8}")
print("-" * 55)
for name, preds in [("Baseline (mean)", y_base_pred), ("Random Forest", y_rf_pred)]:
    mae  = mean_absolute_error(y_reg_test, preds)
    rmse = np.sqrt(mean_squared_error(y_reg_test, preds))
    r2   = r2_score(y_reg_test, preds)
    print(f"{name:<25} {mae:>8.3f} {rmse:>8.3f} {r2:>8.3f}")

WU-2: Reading Cross-Validation Output

Run 5-fold stratified cross-validation on the breast cancer dataset using a LogisticRegression pipeline with StandardScaler. Use cross_validate with return_train_score=True. Print the mean and std for both train and validation F1 scores. Identify whether the model is overfitting, underfitting, or well-fitted.

Show answer

lr_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("lr", LogisticRegression(max_iter=1000, random_state=42))
])

cv_results = cross_validate(
    lr_pipeline,
    X_clf_train,
    y_clf_train,
    cv=StratifiedKFold(n_splits=5, shuffle=True, random_state=42),
    scoring="f1",
    return_train_score=True
)

train_f1 = cv_results["train_score"]
val_f1   = cv_results["test_score"]

print(f"Train F1: {train_f1.mean():.3f} ± {train_f1.std():.3f}")
print(f"Val F1:   {val_f1.mean():.3f} ± {val_f1.std():.3f}")
print(f"Gap:      {train_f1.mean() - val_f1.mean():.3f}")

gap = train_f1.mean() - val_f1.mean()
if val_f1.mean() < 0.7:
    verdict = "Underfitting — both scores are low."
elif gap > 0.05:
    verdict = "Overfitting — large gap between train and validation."
else:
    verdict = "Well-fitted — small gap, good validation score."

print(f"\nVerdict: {verdict}")

WU-3: Confusion Matrix and Classification Report

Train a RandomForestClassifier on the imbalanced dataset. Print the confusion matrix and full classification report. Identify the number of false negatives (missed positives).

Show answer

rf_imb = Pipeline([
    ("scaler", StandardScaler()),
    ("rf", RandomForestClassifier(n_estimators=100, random_state=42))
])
rf_imb.fit(X_imb_train, y_imb_train)
y_imb_pred = rf_imb.predict(X_imb_test)

cm = confusion_matrix(y_imb_test, y_imb_pred)
print("Confusion Matrix:")
print(cm)

tn, fp, fn, tp = cm.ravel()
print(f"\nTrue Negatives  (correctly cleared): {tn}")
print(f"False Positives (false alarms):       {fp}")
print(f"False Negatives (missed positives):   {fn}")
print(f"True Positives  (caught positives):   {tp}")
print(f"\nMissed {fn} out of {fn + tp} actual positives ({fn / (fn + tp):.1%} miss rate)")

print("\nClassification Report:")
print(classification_report(y_imb_test, y_imb_pred, target_names=["Negative", "Positive"]))

---

Main Exercises

These require combining multiple concepts. Expect to write 20–50 lines of code each.

M-1: Plot Learning Curves and Diagnose

Plot learning curves for two models on the breast cancer dataset: - LogisticRegression (may underfit) - RandomForestClassifier with max_depth=None (may overfit)

For each model, plot training and validation F1 score against training set size. Summarise in one sentence what each learning curve tells you about the model's bias-variance behaviour.

Show answer

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
train_sizes = np.linspace(0.1, 1.0, 8)

models_lc = [
    (Pipeline([("scaler", StandardScaler()), ("lr", LogisticRegression(max_iter=1000))]),
     "Logistic Regression", axes[0]),
    (Pipeline([("scaler", StandardScaler()), ("rf", RandomForestClassifier(max_depth=None, random_state=42))]),
     "Random Forest (unbounded depth)", axes[1]),
]

for model, title, ax in models_lc:
    sizes, train_scores, val_scores = learning_curve(
        model,
        X_clf,
        y_clf,
        train_sizes=train_sizes,
        cv=StratifiedKFold(n_splits=5, shuffle=True, random_state=42),
        scoring="f1",
        n_jobs=-1
    )
    t_mean = train_scores.mean(axis=1)
    v_mean = val_scores.mean(axis=1)
    t_std  = train_scores.std(axis=1)
    v_std  = val_scores.std(axis=1)

    ax.plot(sizes, t_mean, label="Train F1",      color="steelblue",  lw=2)
    ax.plot(sizes, v_mean, label="Validation F1", color="darkorange", lw=2)
    ax.fill_between(sizes, t_mean - t_std, t_mean + t_std, alpha=0.15, color="steelblue")
    ax.fill_between(sizes, v_mean - v_std, v_mean + v_std, alpha=0.15, color="darkorange")
    ax.set_title(title)
    ax.set_xlabel("Training set size")
    ax.set_ylabel("F1 Score")
    ax.set_ylim([0.7, 1.05])
    ax.legend()

plt.suptitle("Learning Curves — Breast Cancer Dataset", y=1.02)
plt.tight_layout()
plt.savefig("learning_curves_exercise.png", dpi=150)
plt.show()

# Diagnosis
print("Logistic Regression:")
print("  Both curves converge to ~0.95 with a small gap — well-fitted with low variance.")
print("  More data does not help much; the model has reached its capacity.")
print()
print("Random Forest (unbounded):")
print("  Training F1 ≈ 1.00, validation F1 ≈ 0.95. Gap closes slowly with more data.")
print("  Signs of overfitting — the model memorises training data.")

M-2: Tune a Classifier with RandomizedSearchCV

Use RandomizedSearchCV on the imbalanced dataset to tune a GradientBoostingClassifier. Search over at least 4 hyperparameters with distributions (not lists). Use n_iter=30, cv=5, and score on F1.

Report: 1. The best hyperparameters 2. The best CV F1 score 3. The test set F1, precision, and recall at the default 0.5 threshold

Show answer

gb_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("gb", GradientBoostingClassifier(random_state=42))
])

param_distributions = {
    "gb__n_estimators":      randint(50, 300),
    "gb__max_depth":         randint(2, 7),
    "gb__learning_rate":     uniform(0.01, 0.29),
    "gb__subsample":         uniform(0.6, 0.4),
    "gb__min_samples_leaf":  randint(1, 15)
}

random_search = RandomizedSearchCV(
    gb_pipeline,
    param_distributions,
    n_iter=30,
    cv=StratifiedKFold(n_splits=5, shuffle=True, random_state=42),
    scoring="f1",
    refit=True,
    n_jobs=-1,
    random_state=42,
    verbose=0
)

random_search.fit(X_imb_train, y_imb_train)

print(f"Best CV F1:      {random_search.best_score_:.3f}")
print(f"Best parameters: {random_search.best_params_}")

best_gb = random_search.best_estimator_
y_pred_gb = best_gb.predict(X_imb_test)

print(f"\nTest set results (threshold = 0.5):")
print(f"  F1:        {f1_score(y_imb_test, y_pred_gb):.3f}")
print(f"  Precision: {precision_score(y_imb_test, y_pred_gb):.3f}")
print(f"  Recall:    {recall_score(y_imb_test, y_pred_gb):.3f}")

M-3: ROC and PR Curves Side by Side

For the imbalanced dataset, compare a LogisticRegression and the best GradientBoostingClassifier from M-2. Plot both their ROC curves and PR curves on the same figures. Print AUC and PR-AUC for each.

Discuss in a comment: which model is better, and does ROC-AUC or PR-AUC give a more informative comparison for this dataset?

Show answer

# Train logistic regression for comparison
lr_imb = Pipeline([
    ("scaler", StandardScaler()),
    ("lr", LogisticRegression(max_iter=1000, random_state=42))
])
lr_imb.fit(X_imb_train, y_imb_train)

models_compare = {
    "Logistic Regression": lr_imb,
    "Gradient Boosting":   best_gb
}
colors = {"Logistic Regression": "steelblue", "Gradient Boosting": "darkorange"}

fig, (ax_roc, ax_pr) = plt.subplots(1, 2, figsize=(14, 5))
baseline_prev = y_imb_test.mean()

for model_name, model in models_compare.items():
    y_proba = model.predict_proba(X_imb_test)[:, 1]
    color   = colors[model_name]

    # ROC
    fpr, tpr, _ = roc_curve(y_imb_test, y_proba)
    auc = roc_auc_score(y_imb_test, y_proba)
    ax_roc.plot(fpr, tpr, color=color, lw=2, label=f"{model_name} (AUC={auc:.3f})")

    # PR
    prec, rec, _ = precision_recall_curve(y_imb_test, y_proba)
    pr_auc = average_precision_score(y_imb_test, y_proba)
    ax_pr.plot(rec, prec, color=color, lw=2, label=f"{model_name} (PR-AUC={pr_auc:.3f})")

    print(f"{model_name}:  ROC-AUC = {auc:.3f},  PR-AUC = {pr_auc:.3f}")

ax_roc.plot([0, 1], [0, 1], "k--", lw=1, label="Random")
ax_roc.set_xlabel("False Positive Rate")
ax_roc.set_ylabel("True Positive Rate")
ax_roc.set_title("ROC Curves")
ax_roc.legend()

ax_pr.axhline(y=baseline_prev, color="gray", linestyle="--", lw=1,
              label=f"Baseline (prevalence={baseline_prev:.2f})")
ax_pr.set_xlabel("Recall")
ax_pr.set_ylabel("Precision")
ax_pr.set_title("Precision-Recall Curves")
ax_pr.legend()

plt.tight_layout()
plt.savefig("roc_pr_comparison.png", dpi=150)
plt.show()

# Discussion:
# On this imbalanced dataset (5% positive), ROC-AUC values for both models
# are high and look similar. PR-AUC shows a larger gap between models because
# it focuses entirely on the positive class — where the models actually differ.
# PR-AUC is the more informative metric here.

M-4: Threshold Optimisation

Using the best model from M-2, scan thresholds from 0.1 to 0.9 in steps of 0.05 on the test set predictions. Create a table of precision, recall, and F1 at each threshold.

Then answer: if the business says "we cannot afford to miss more than 20% of positives" (i.e., recall >= 0.80), what is the highest precision you can achieve while meeting that constraint?

Show answer

y_proba_gb = best_gb.predict_proba(X_imb_test)[:, 1]

threshold_results = []
for thresh in np.arange(0.05, 0.95, 0.05):
    y_thresh = (y_proba_gb >= thresh).astype(int)
    prec = precision_score(y_imb_test, y_thresh, zero_division=0)
    rec  = recall_score(y_imb_test, y_thresh)
    f1   = f1_score(y_imb_test, y_thresh, zero_division=0)
    threshold_results.append({
        "threshold": round(thresh, 2),
        "precision": round(prec, 3),
        "recall":    round(rec, 3),
        "f1":        round(f1, 3)
    })

df_thresh = pd.DataFrame(threshold_results)
print(df_thresh.to_string(index=False))

# Business constraint: recall >= 0.80
valid = df_thresh[df_thresh["recall"] >= 0.80]
if not valid.empty:
    best_constrained = valid.loc[valid["precision"].idxmax()]
    print(f"\nWith recall >= 0.80 constraint:")
    print(f"  Best threshold: {best_constrained['threshold']}")
    print(f"  Precision:      {best_constrained['precision']}")
    print(f"  Recall:         {best_constrained['recall']}")
    print(f"  F1:             {best_constrained['f1']}")
else:
    print("\nNo threshold achieves recall >= 0.80 on this dataset.")

# Note: the threshold decision must be documented.
# In a production setting, you select the threshold on validation data,
# not test data. Here we do it on test data for illustration only.

---

Stretch Exercises

These require you to go beyond the notes and think critically. There is no single correct answer.

S-1: Nested Cross-Validation vs Simple GridSearch

On the breast cancer dataset, compare two estimates of model performance:

1. Simple GridSearchCV: GridSearchCV with cv=5 on X_clf_train, then report best_score_.
2. Nested CV: outer KFold(n_splits=5) wrapping an inner GridSearchCV(cv=3).

Run both, print both scores, and explain in a comment why they differ and which one is the honest estimate.

Show answer

from sklearn.model_selection import KFold

rf_clf = Pipeline([
    ("scaler", StandardScaler()),
    ("rf", RandomForestClassifier(random_state=42))
])

param_grid = {
    "rf__n_estimators": [50, 100],
    "rf__max_depth": [None, 5, 10]
}

# --- Simple GridSearchCV ---
simple_grid = GridSearchCV(rf_clf, param_grid, cv=5, scoring="f1", refit=True)
simple_grid.fit(X_clf_train, y_clf_train)
simple_score = simple_grid.best_score_

# --- Nested CV ---
inner_cv  = KFold(n_splits=3, shuffle=True, random_state=42)
outer_cv  = KFold(n_splits=5, shuffle=True, random_state=42)
inner_search = GridSearchCV(rf_clf, param_grid, cv=inner_cv, scoring="f1")
nested_scores = cross_val_score(inner_search, X_clf_train, y_clf_train,
                                cv=outer_cv, scoring="f1")

print(f"Simple GridSearchCV best CV F1: {simple_score:.3f}")
print(f"Nested CV F1:                   {nested_scores.mean():.3f} ± {nested_scores.std():.3f}")
print(f"Optimism (simple - nested):     {simple_score - nested_scores.mean():.3f}")

# Explanation:
# The simple GridSearchCV score is optimistic because the same folds were used
# both to select the best hyperparameters and to estimate performance.
# The model that "wins" is the one that happened to do best on those folds —
# which includes some selection noise.
# Nested CV removes this selection bias by using separate folds for tuning (inner)
# and evaluation (outer). The nested CV score is the honest estimate of
# how the best-tuned model will perform on truly unseen data.

S-2: Write a Model Evaluation Report

Train the best pipeline you can on the California housing dataset (use RandomizedSearchCV on a RandomForestRegressor). Then write a structured model evaluation report as a printed summary covering:

• Problem statement
• Baseline performance (MAE and R²)
• Model and tuning approach
• Best hyperparameters
• Cross-validated performance
• Test set performance
• Residual analysis: is the error distribution symmetric? Are errors larger for high-value houses?
• Recommendation: is the model good enough to deploy?

There is no fixed answer — the quality lies in how thoroughly you investigate and explain.

Show answer (structure)

from scipy.stats import randint, uniform

# Baseline
baseline_r = DummyRegressor(strategy="mean")
baseline_r.fit(X_reg_train, y_reg_train)
base_mae = mean_absolute_error(y_reg_test, baseline_r.predict(X_reg_test))
base_r2  = r2_score(y_reg_test, baseline_r.predict(X_reg_test))

# Tuned Random Forest
rf_reg_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("rf", RandomForestRegressor(random_state=42))
])
param_dist = {
    "rf__n_estimators":     randint(100, 400),
    "rf__max_depth":        randint(5, 20),
    "rf__min_samples_leaf": randint(1, 10),
    "rf__max_features":     uniform(0.3, 0.7)
}
rs = RandomizedSearchCV(rf_reg_pipeline, param_dist, n_iter=30,
                        cv=5, scoring="neg_mean_absolute_error",
                        refit=True, n_jobs=-1, random_state=42)
rs.fit(X_reg_train, y_reg_train)

y_reg_pred = rs.predict(X_reg_test)
test_mae  = mean_absolute_error(y_reg_test, y_reg_pred)
test_rmse = np.sqrt(mean_squared_error(y_reg_test, y_reg_pred))
test_r2   = r2_score(y_reg_test, y_reg_pred)
residuals = y_reg_test - y_reg_pred

# Residual analysis
high_value_mask = y_reg_test > np.percentile(y_reg_test, 75)
mae_high = mean_absolute_error(y_reg_test[high_value_mask], y_reg_pred[high_value_mask])
mae_low  = mean_absolute_error(y_reg_test[~high_value_mask], y_reg_pred[~high_value_mask])

print("=" * 60)
print("MODEL EVALUATION REPORT")
print("=" * 60)
print(f"\nProblem:         Predict median house value (in $100k)")
print(f"Dataset:         California Housing ({X_reg_train.shape[0]} train, {X_reg_test.shape[0]} test rows)")
print(f"\nBaseline (predict mean):")
print(f"  MAE:  {base_mae:.3f}  (${base_mae * 100_000:,.0f})")
print(f"  R²:   {base_r2:.3f}")
print(f"\nModel: RandomForestRegressor via RandomizedSearchCV (30 iterations, 5-fold CV)")
print(f"  Best hyperparameters: {rs.best_params_}")
print(f"  Best CV MAE: {-rs.best_score_:.3f}  (${-rs.best_score_ * 100_000:,.0f})")
print(f"\nTest Set Performance:")
print(f"  MAE:  {test_mae:.3f}  (${test_mae * 100_000:,.0f})")
print(f"  RMSE: {test_rmse:.3f} (${test_rmse * 100_000:,.0f})")
print(f"  R²:   {test_r2:.3f}")
print(f"\nResidual Analysis:")
print(f"  Mean residual:  {residuals.mean():.4f}  (near zero = unbiased)")
print(f"  Skewness:       {pd.Series(residuals).skew():.3f}")
print(f"  MAE (top 25%):  {mae_high:.3f}  (${mae_high * 100_000:,.0f})")
print(f"  MAE (bot 75%):  {mae_low:.3f}   (${mae_low * 100_000:,.0f})")
print(f"\nRecommendation:")
if test_r2 > 0.80 and test_mae < 0.40:
    print("  Model is strong (R² > 0.80, MAE < $40k). Suitable for deployment")
    print("  with monitoring. Note higher errors on high-value properties.")
else:
    print("  Model needs improvement before deployment.")
print("=" * 60)

S-3: Time Series Cross-Validation

Generate 3 years of monthly sales data with a trend and seasonality component. Fit a Ridge regression predicting next-month sales from lagged features (previous 3 months). Evaluate using TimeSeriesSplit(n_splits=6) and compare to a standard KFold split. Show why the K-Fold MAE is optimistic.

Show answer

from sklearn.linear_model import Ridge

# Generate monthly sales data: trend + seasonality + noise
np.random.seed(42)
n_months = 48  # 4 years
months = np.arange(n_months)
sales = (
    1000                                         # base
    + months * 5                                 # trend
    + 150 * np.sin(2 * np.pi * months / 12)     # seasonality
    + np.random.normal(0, 30, n_months)          # noise
)

# Build lag features: predict month t using months t-1, t-2, t-3
lags = 3
X_ts = np.array([sales[i:n_months - lags + i] for i in range(lags)]).T
y_ts = sales[lags:]

ridge = Ridge(alpha=1.0)

# TimeSeriesSplit — respects chronological order
tscv = TimeSeriesSplit(n_splits=6)
ts_scores = cross_val_score(ridge, X_ts, y_ts,
                            cv=tscv, scoring="neg_mean_absolute_error")
ts_mae = -ts_scores.mean()

# Standard KFold — shuffles the data (wrong for time series)
kfold = KFold(n_splits=6, shuffle=True, random_state=42)    # type: ignore
kf_scores = cross_val_score(ridge, X_ts, y_ts,
                            cv=kfold, scoring="neg_mean_absolute_error")
kf_mae = -kf_scores.mean()

print(f"TimeSeriesSplit MAE: {ts_mae:.1f}")
print(f"KFold (shuffled) MAE: {kf_mae:.1f}")
print(f"Optimism (KFold - TS): {ts_mae - kf_mae:.1f}")
print()
print("KFold MAE is lower because it uses future data to predict the past.")
print("The model 'knows' about the trend direction in its training fold.")
print("TimeSeriesSplit gives the honest estimate — training only uses past data.")

---

Self-Check Questions

Answer these without code. If you cannot answer them confidently, go back to the relevant note.

1. A model achieves ROC-AUC = 0.97 on a dataset where 99% of cases are class 0. What additional information do you need before concluding the model is good?

2. Your cross-validation result is F1 = 0.84 ± 0.09. Your colleague's model achieves F1 = 0.87 ± 0.02. Which model would you prefer to deploy and why?

3. You run GridSearchCV with cv=5 on 20 hyperparameter combinations. The best combination achieves validation F1 = 0.93. Is 0.93 an unbiased estimate of how the model will perform on truly unseen data? Why or why not?

4. Explain in one sentence why scaling inside a Pipeline is different from scaling before splitting.

5. You are building a model to detect a rare disease (0.1% prevalence). Your manager asks for "accuracy". What do you say?

---

Previous: Model Selection and Tuning | Next: Intro to Deep Learning