🎯 What You'll Learn
- Use probability scores to evaluate a model's ranking ability, independent of any threshold
- Build and interpret ROC curves: True Positive Rate vs False Positive Rate at all thresholds
- Understand AUC as the probability that the model ranks a random positive higher than a random negative
- Perform cost-benefit analysis on the ROC curve to find the optimal business threshold
- Implement Stratified K-Fold cross-validation to prevent data leakage and unreliable estimates
- Test whether one model's cross-validation score is really better than another's with a paired t-test and confidence interval, not just eyeballing the averages
- Check whether a model's predicted probabilities are calibrated — and fix them with
CalibratedClassifierCVwhen they aren't
1 Beyond Hard Labels: Probability Scores
In Lesson 19, we evaluated models using hard 0/1 predictions. But most classifiers — logistic regression, SVM with probability=True, random forests — output a probability score before thresholding. Evaluating these scores directly, rather than their thresholded labels, reveals more about the model's true capability.
The probability score tells us how confident the model is. A sample with P(y=1|x) = 0.95 and a sample with P(y=1|x) = 0.55 both become the same "1" label when we threshold at 0.5 — but they represent very different levels of certainty. Metrics that use raw probabilities can capture this distinction.
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
cancer = load_breast_cancer()
X, y = cancer.data, cancer.target # 0=malignant, 1=benign
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.25, random_state=42, stratify=y)
pipe = Pipeline([
('scaler', StandardScaler()),
('lr', LogisticRegression(max_iter=1000, random_state=42))
])
pipe.fit(X_tr, y_tr)
# Raw probability scores — shape (n_samples, n_classes)
y_proba = pipe.predict_proba(X_te)
y_scores = y_proba[:, 1] # P(y=1|x) for each test sample
# Inspect the distribution of scores
print("Distribution of probability scores:")
for threshold in [0.1, 0.3, 0.5, 0.7, 0.9]:
n_positive = (y_scores >= threshold).sum()
print(f" Scores >= {threshold}: {n_positive}/{len(y_scores)} samples predicted positive")
# Show score vs actual label for first 10 test samples
print("\nFirst 10 test samples: score vs true label")
print(f"{'Score':>8} {'Predicted':>10} {'Actual':>8}")
for score, actual in zip(y_scores[:10], y_te[:10]):
pred = "positive" if score >= 0.5 else "negative"
actual_str = "positive" if actual == 1 else "negative"
print(f"{score:>8.4f} {pred:>10} {actual_str:>8}")2 The Precision-Recall Tradeoff
Changing the classification threshold directly changes the tradeoff between precision and recall. Lower threshold → predict positive more often → higher recall, lower precision. Higher threshold → predict positive only when very confident → lower recall, higher precision. sklearn provides precision_recall_curve() to sweep all thresholds at once.
from sklearn.metrics import precision_recall_curve, average_precision_score
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
import matplotlib.pyplot as plt
cancer = load_breast_cancer()
X, y = cancer.data, 1 - cancer.target # flip: malignant=1 (the "positive" case we care about)
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.25, random_state=42, stratify=y)
pipe = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
pipe.fit(X_tr, y_tr)
y_scores = pipe.predict_proba(X_te)[:, 1]
# Precision-Recall curve: arrays of (precision, recall, threshold) at every threshold
precisions, recalls, thresholds = precision_recall_curve(y_te, y_scores)
# Average Precision (area under PR curve) — useful summary for imbalanced data
ap = average_precision_score(y_te, y_scores)
print(f"Average Precision (AP): {ap:.4f}")
# Plot the tradeoff
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# Left: PR curve
axes[0].plot(recalls, precisions, 'b-', linewidth=2)
axes[0].set_xlabel('Recall')
axes[0].set_ylabel('Precision')
axes[0].set_title(f'Precision-Recall Curve (AP={ap:.3f})')
axes[0].grid(True, alpha=0.3)
# Right: Precision and Recall vs threshold
axes[1].plot(thresholds, precisions[:-1], 'b-', label='Precision')
axes[1].plot(thresholds, recalls[:-1], 'r-', label='Recall')
axes[1].set_xlabel('Classification Threshold')
axes[1].set_ylabel('Score')
axes[1].set_title('Precision and Recall vs Threshold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Find threshold that maximises F1
f1_scores = 2 * precisions[:-1] * recalls[:-1] / (precisions[:-1] + recalls[:-1] + 1e-10)
best_idx = np.argmax(f1_scores)
print(f"Best threshold for F1: {thresholds[best_idx]:.3f}")
print(f" Precision: {precisions[best_idx]:.4f}")
print(f" Recall: {recalls[best_idx]:.4f}")
print(f" F1: {f1_scores[best_idx]:.4f}")3 The ROC Curve
The ROC (Receiver Operating Characteristic) curve plots the True Positive Rate (TPR) against the False Positive Rate (FPR) at every possible classification threshold. By sweeping from threshold=1.0 (predict nothing positive) to threshold=0.0 (predict everything positive), we trace out the curve.
- TPR = Recall = TP/(TP+FN): Fraction of actual positives that are correctly detected.
- FPR = FP/(FP+TN): Fraction of actual negatives that are incorrectly flagged as positive (the "false alarm rate").
Key points on the ROC curve:
- (0, 0): Threshold=1.0 — predict nothing positive. TPR=0, FPR=0.
- (1, 1): Threshold=0.0 — predict everything positive. TPR=1, FPR=1.
- (0, 1): Perfect classifier — all positives detected, no false alarms.
- Diagonal (0,0)→(1,1): Random classifier — no better than chance.
The ROC curve is built by sweeping the classification threshold from 1.0 down to 0.0. Every threshold yields one (FPR, TPR) point — and one confusion matrix. Drag the slider below to pick a threshold on a synthetic 50-sample dataset (25 positives, 25 negatives) with overlapping score distributions: watch the marker slide along the curve, and the confusion matrix update, in lockstep.
At threshold = 0.50, the operating point sits here on the ROC curve.
Confusion matrix at the current threshold.
from sklearn.metrics import roc_curve, roc_auc_score, RocCurveDisplay
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt
import numpy as np
cancer = load_breast_cancer()
X, y = cancer.data, cancer.target
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.25, random_state=42, stratify=y)
# Train two models — RandomForest is a Lesson 23 model, used here only
# as "a second classifier to compare against" (black box is fine for now)
lr_pipe = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
rf_pipe = Pipeline([('s', StandardScaler()),
('rf', RandomForestClassifier(n_estimators=100, random_state=42))])
lr_pipe.fit(X_tr, y_tr)
rf_pipe.fit(X_tr, y_tr)
lr_scores = lr_pipe.predict_proba(X_te)[:, 1]
rf_scores = rf_pipe.predict_proba(X_te)[:, 1]
# ROC curves
lr_fpr, lr_tpr, lr_thresholds = roc_curve(y_te, lr_scores)
rf_fpr, rf_tpr, rf_thresholds = roc_curve(y_te, rf_scores)
lr_auc = roc_auc_score(y_te, lr_scores)
rf_auc = roc_auc_score(y_te, rf_scores)
# Plot
fig, ax = plt.subplots(figsize=(7, 5))
ax.plot(lr_fpr, lr_tpr, 'b-', linewidth=2, label=f'Logistic Regression (AUC={lr_auc:.3f})')
ax.plot(rf_fpr, rf_tpr, 'g-', linewidth=2, label=f'Random Forest (AUC={rf_auc:.3f})')
ax.plot([0, 1], [0, 1], 'k--', label='Random classifier (AUC=0.5)')
ax.scatter([0], [1], s=100, color='gold', zorder=5, label='Perfect classifier')
ax.set_xlabel('False Positive Rate (FPR)')
ax.set_ylabel('True Positive Rate (TPR = Recall)')
ax.set_title('ROC Curves: Logistic Regression vs Random Forest')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Logistic Regression AUC: {lr_auc:.4f}")
print(f"Random Forest AUC: {rf_auc:.4f}")The ROC curve is the standard choice when classes are approximately balanced. For severely imbalanced data (e.g., 1% positive), the ROC curve can look optimistic because FPR = FP/(FP+TN) — with many negatives, even a high FP count gives a low FPR. In this case, use the Precision-Recall curve instead, which better reflects performance on the minority class.
4 AUC: Area Under the Curve
The AUC (Area Under the ROC Curve) summarises the ROC curve as a single number between 0 and 1. It has a powerful probabilistic interpretation:
"AUC is the probability that the model assigns a higher score to a randomly chosen positive sample than to a randomly chosen negative sample."
An AUC of 0.85 means: pick any positive sample and any negative sample at random. The model correctly ranks the positive higher 85% of the time. AUC = 0.5 means the model is no better than random ranking. AUC = 1.0 means the model perfectly separates all positives from negatives.
AUC is threshold-independent — it evaluates the model's ranking ability regardless of where you set the decision threshold. This makes it useful for comparing models before you've decided on an operational threshold.
from sklearn.metrics import roc_auc_score
import numpy as np
# Verify the probabilistic interpretation manually
def auc_by_simulation(y_true, y_scores, n_simulations=100_000):
"""
Manually estimate AUC by randomly sampling pairs.
AUC = P(score_positive > score_negative)
"""
positives = y_scores[y_true == 1]
negatives = y_scores[y_true == 0]
rng = np.random.default_rng(42)
pos_samples = rng.choice(positives, size=n_simulations, replace=True)
neg_samples = rng.choice(negatives, size=n_simulations, replace=True)
# Handle ties: count as 0.5
correct = (pos_samples > neg_samples).sum()
tied = (pos_samples == neg_samples).sum()
return (correct + 0.5 * tied) / n_simulations
# Generate scores for demonstration
np.random.seed(42)
y_true = np.array([1] * 50 + [0] * 50)
y_scores = np.concatenate([
np.random.beta(8, 3, 50), # positives: scores concentrated near 1
np.random.beta(3, 8, 50) # negatives: scores concentrated near 0
])
sklearn_auc = roc_auc_score(y_true, y_scores)
manual_auc = auc_by_simulation(y_true, y_scores)
print(f"sklearn roc_auc_score: {sklearn_auc:.4f}")
print(f"Simulation (100K pairs): {manual_auc:.4f}")
print(f"Interpretation: The model ranks a random positive above a random negative {sklearn_auc:.1%} of the time")
# What different AUC ranges mean
print("\nAUC interpretation guide:")
for auc, label in [(0.5, "Random / useless"), (0.6, "Poor"),
(0.7, "Fair"), (0.8, "Good"), (0.9, "Excellent"), (0.99, "Near-perfect")]:
bar = '█' * int((auc - 0.5) * 40)
print(f" AUC={auc:.2f} {label:20} {bar}")5 Cost-Benefit Analysis Using the ROC Curve
The ROC curve becomes directly actionable when you assign business costs to false positives and false negatives. For each threshold on the ROC curve, you can compute the expected value of using that threshold, then choose the one that maximises it (or minimises expected cost).
Consider a fraud detection system:
- True Positive: Fraud caught → save $500 in chargeback losses
- False Positive: Legitimate transaction blocked → cost $10 in customer service + friction
- False Negative: Fraud missed → lose $500
- True Negative: No action needed → $0 cost
from sklearn.metrics import roc_curve
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
import matplotlib.pyplot as plt
# Simulate a fraud detection dataset (5% fraud rate)
X, y = make_classification(n_samples=10000, weights=[0.95, 0.05],
n_features=20, n_informative=10,
random_state=42)
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.3, stratify=y, random_state=42)
pipe = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
pipe.fit(X_tr, y_tr)
y_scores = pipe.predict_proba(X_te)[:, 1]
fpr_arr, tpr_arr, thresholds = roc_curve(y_te, y_scores)
# Business parameters
n_test = len(y_te)
n_positives = y_te.sum()
n_negatives = n_test - n_positives
cost_fp = 10 # cost per false positive (blocked legit transaction)
cost_fn = 500 # cost per false negative (fraud loss)
benefit_tp = 0 # we're minimizing cost, so TP avoidance is implicit in FN cost
# Compute expected total cost at each threshold
def expected_cost(tpr, fpr, n_pos, n_neg, cost_fn, cost_fp):
fn = (1 - tpr) * n_pos # expected false negatives
fp = fpr * n_neg # expected false positives
return fn * cost_fn + fp * cost_fp
costs = [expected_cost(tpr, fpr, n_positives, n_negatives, cost_fn, cost_fp)
for tpr, fpr in zip(tpr_arr, fpr_arr)]
# Find optimal threshold
best_idx = np.argmin(costs)
best_threshold = thresholds[best_idx]
best_cost = costs[best_idx]
best_tpr = tpr_arr[best_idx]
best_fpr = fpr_arr[best_idx]
print(f"=== Cost-Benefit Analysis ===")
print(f"FP cost: ${cost_fp} per false positive")
print(f"FN cost: ${cost_fn} per false negative")
print(f"\nOptimal threshold: {best_threshold:.4f}")
print(f"At this threshold:")
print(f" TPR (recall): {best_tpr:.4f} — catches {best_tpr:.0%} of frauds")
print(f" FPR: {best_fpr:.4f} — {best_fpr:.0%} of legit transactions flagged")
print(f" Expected total cost on test set: ${best_cost:,.0f}")
# Compare with default threshold 0.5
default_idx = np.argmin(np.abs(thresholds - 0.5))
default_cost = costs[default_idx]
print(f"\nCost at default threshold (0.5): ${default_cost:,.0f}")
print(f"Savings from optimal threshold: ${default_cost - best_cost:,.0f}")
# Plot cost vs threshold
plt.figure(figsize=(8, 4))
plt.plot(thresholds, costs, 'b-', linewidth=2)
plt.axvline(x=best_threshold, color='red', linestyle='--',
label=f'Optimal threshold={best_threshold:.3f}')
plt.axvline(x=0.5, color='gray', linestyle='--', label='Default threshold=0.5')
plt.xlabel('Classification Threshold')
plt.ylabel('Expected Total Cost ($)')
plt.title('Expected Cost vs Threshold — Fraud Detection')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()6 K-Fold Cross-Validation
A single train/test split is noisy. The test set might be unusually easy or unusually hard by chance. K-fold cross-validation solves this by evaluating the model K times on different subsets of the data and averaging the results.
The procedure: divide the dataset into K equal-sized folds. For each fold i: train on all K-1 folds, evaluate on fold i. Average the K evaluation scores. This gives a more stable and unbiased estimate of generalization performance.
from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
cancer = load_breast_cancer()
X, y = cancer.data, cancer.target
# Single train/test split — noisy estimate
from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)
pipe_single = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
pipe_single.fit(X_tr, y_tr)
single_acc = pipe_single.score(X_te, y_te)
print(f"Single split accuracy (test_size=0.2, seed=42): {single_acc:.4f}")
# Try different random seeds to show variability
print("\nAccuracy with 5 different random seeds:")
for seed in [0, 1, 2, 3, 4]:
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=seed)
pipe = Pipeline([('s', StandardScaler()), ('lr', LogisticRegression(max_iter=1000))])
pipe.fit(X_tr, y_tr)
print(f" seed={seed}: {pipe.score(X_te, y_te):.4f}")
# K-fold cross-validation — stable estimate
pipe_cv = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
for k in [3, 5, 10]:
kf = KFold(n_splits=k, shuffle=True, random_state=42)
scores = cross_val_score(pipe_cv, X, y, cv=kf, scoring='accuracy')
print(f"\n{k}-Fold CV accuracy: {scores.mean():.4f} ± {scores.std():.4f}")
print(f" Per-fold scores: {scores.round(4)}")7 Stratified K-Fold
Standard K-fold splits data randomly. For imbalanced datasets, this can produce folds with wildly different class proportions — or even folds with zero positive samples, causing metrics like precision and AUC to be undefined or misleading.
Stratified K-fold ensures that each fold has approximately the same proportion of each class as the overall dataset. This is the correct choice for classification problems and should be used by default over standard K-fold.
sklearn's cross_val_score with a classifier automatically uses stratified splitting when you pass an integer to cv=. But if you create a KFold object explicitly, it does NOT stratify. To be explicit and safe: always use StratifiedKFold for classification tasks, especially with imbalanced data.
The diagram below shows 20 samples (16 majority-class circles, 4 minority-class triangles — an 80/20 split) divided into 4 folds two ways. With standard KFold, samples are assigned to folds in whatever order a random shuffle happens to produce, so the minority class can land unevenly — even entirely missing a fold. With StratifiedKFold, every fold gets the same 80/20 ratio:
Same 20 samples, two splitting strategies. Standard KFold's random assignment can starve a fold of minority-class examples (here, Fold 1 gets none). StratifiedKFold guarantees each fold reflects the dataset's overall class balance.
from sklearn.model_selection import KFold, StratifiedKFold, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import numpy as np
# Severely imbalanced dataset: 95% class 0, 5% class 1
X, y = make_classification(n_samples=1000, weights=[0.95, 0.05],
n_features=20, random_state=42)
print(f"Dataset: {(y==0).sum()} class 0, {(y==1).sum()} class 1")
print(f"Class 1 rate: {y.mean():.2%}")
# Show class proportions per fold
print("\n=== Standard KFold (no stratification) ===")
kf = KFold(n_splits=5, shuffle=True, random_state=42)
for fold, (tr_idx, te_idx) in enumerate(kf.split(X, y)):
fold_rate = y[te_idx].mean()
print(f" Fold {fold+1}: {len(te_idx)} samples, class 1 rate = {fold_rate:.2%}")
print("\n=== StratifiedKFold ===")
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
for fold, (tr_idx, te_idx) in enumerate(skf.split(X, y)):
fold_rate = y[te_idx].mean()
print(f" Fold {fold+1}: {len(te_idx)} samples, class 1 rate = {fold_rate:.2%}")
# Compare AUC estimates: standard vs stratified
pipe = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))])
scores_kf = cross_val_score(pipe, X, y, cv=kf, scoring='roc_auc')
scores_skf = cross_val_score(pipe, X, y, cv=skf, scoring='roc_auc')
print(f"\nKFold AUC: {scores_kf.mean():.4f} ± {scores_kf.std():.4f}")
print(f"StratKFold AUC: {scores_skf.mean():.4f} ± {scores_skf.std():.4f}")
print(f"\nStratifiedKFold has lower std (more stable estimates) on imbalanced data")8 Cross-Validation Pitfalls and Pipelines
The most dangerous mistake in cross-validation is data leakage through preprocessing. If you fit a StandardScaler (or imputer, or feature selector) on the entire dataset before cross-validation, the test fold's statistics (mean, std) contaminate the scaler's fit. The model has effectively "seen" the test data through the scaler. This inflates performance estimates.
The solution is always to wrap preprocessing and model inside a Pipeline. sklearn's cross-validation runs the Pipeline's fit on the training fold and transform on the test fold in each iteration — exactly what you want.
from sklearn.model_selection import cross_val_score, StratifiedKFold, cross_validate
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.datasets import load_breast_cancer
import numpy as np
cancer = load_breast_cancer()
X, y = cancer.data, cancer.target
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# ---- WRONG: leakage ----
# Fitting scaler on ALL data before cross-validation
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X) # test fold data leaks into scaler here!
lr = LogisticRegression(max_iter=1000)
wrong_scores = cross_val_score(lr, X_scaled, y, cv=skf, scoring='accuracy')
# ---- CORRECT: no leakage ----
# Scaler is inside the pipeline; it's fit only on training folds
pipe = Pipeline([
('scaler', StandardScaler()),
('lr', LogisticRegression(max_iter=1000))
])
correct_scores = cross_val_score(pipe, X, y, cv=skf, scoring='accuracy')
print("=== Data Leakage Demonstration ===")
print(f"WRONG (leaky): {wrong_scores.mean():.5f} ± {wrong_scores.std():.5f}")
print(f"CORRECT (pipeline): {correct_scores.mean():.5f} ± {correct_scores.std():.5f}")
print(f"Leakage inflates accuracy by: {(wrong_scores.mean() - correct_scores.mean())*100:.3f} pp")
# cross_validate: multiple metrics at once
from sklearn.metrics import make_scorer, roc_auc_score
metrics = {
'accuracy': 'accuracy',
'roc_auc': 'roc_auc',
'f1': 'f1',
'precision': 'precision',
'recall': 'recall'
}
cv_results = cross_validate(pipe, X, y, cv=skf, scoring=metrics, return_train_score=True)
print("\n=== Multiple Metrics via cross_validate ===")
for metric in metrics:
test_mean = cv_results[f'test_{metric}'].mean()
train_mean = cv_results[f'train_{metric}'].mean()
print(f" {metric:12} Train={train_mean:.4f} Test={test_mean:.4f}")Every preprocessing step that looks at data to compute statistics — StandardScaler, MinMaxScaler, SimpleImputer, SelectKBest, PCA — must be inside a Pipeline when used with cross-validation. Steps that apply a fixed transformation (like adding a polynomial feature) don't need to be in a Pipeline, but it's good practice to always use one anyway.
9 Is Model B Actually Better? Statistical Significance
Section 6 gave you K-Fold cross-validation scores — say, Model A averages 0.842 accuracy and Model B averages 0.851. Model B looks better. But is an 0.9-point gap a real, reproducible improvement, or could it just as easily be noise from which rows happened to land in which fold? This is exactly the question the Central Limit Theorem (Lesson 7, Section 9) was building toward: comparing two averages requires knowing how much those averages would wobble on a different random split.
The standard tool is a paired t-test on the per-fold scores — "paired" because both models were evaluated on the same folds, so you're testing whether the fold-by-fold differences are consistently positive, not just comparing two unrelated averages:
import numpy as np
from scipy import stats
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
X, y = make_classification(n_samples=500, n_features=20, n_informative=8, random_state=42)
skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
# Evaluate BOTH models on the exact same 10 folds -- this is what makes
# the comparison "paired" rather than two independent experiments.
scores_a = cross_val_score(LogisticRegression(max_iter=1000), X, y, cv=skf, scoring='accuracy')
scores_b = cross_val_score(RandomForestClassifier(n_estimators=200, random_state=42), X, y, cv=skf, scoring='accuracy')
print(f"Model A (Logistic Regression): {scores_a.mean():.4f} ± {scores_a.std():.4f}")
print(f"Model B (Random Forest): {scores_b.mean():.4f} ± {scores_b.std():.4f}")
print(f"Per-fold differences (B - A): {(scores_b - scores_a).round(4)}")
# Paired t-test: is the mean of (scores_b - scores_a) significantly != 0?
t_stat, p_value = stats.ttest_rel(scores_b, scores_a)
print(f"\nPaired t-test: t={t_stat:.3f}, p={p_value:.4f}")
alpha = 0.05
if p_value < alpha:
print(f"p < {alpha} -- the difference is statistically significant. Model B is reliably better.")
else:
print(f"p >= {alpha} -- cannot reject the null hypothesis that A and B perform equally.")
print("The observed gap may just be fold-to-fold noise, not a real improvement.")The p-value answers a precise question: "if the two models truly performed identically, how likely would we be to see a difference this large (or larger) by chance alone?" A small p-value (conventionally below 0.05) means the observed gap would be unusual under pure chance — evidence the difference is real. It is not the probability that Model B is better, a common and important misreading.
A confidence interval complements the p-value by quantifying the gap itself, not just whether it's non-zero:
import numpy as np
from scipy import stats
diffs = scores_b - scores_a
n = len(diffs)
mean_diff = diffs.mean()
se_diff = diffs.std(ddof=1) / np.sqrt(n) # standard error, per the CLT
# 95% confidence interval for the true mean difference
ci_low, ci_high = stats.t.interval(0.95, df=n-1, loc=mean_diff, scale=se_diff)
print(f"Mean improvement: {mean_diff:.4f}")
print(f"95% CI: [{ci_low:.4f}, {ci_high:.4f}]")
# If the interval excludes 0, that agrees with a significant p-value --
# and the WIDTH of the interval tells you how precisely you've pinned down
# the improvement, which a p-value alone never communicates.With enough folds or enough data, even a genuinely tiny, operationally meaningless improvement (0.851 vs 0.850) can become "statistically significant." Always report the confidence interval alongside the p-value, and ask whether the size of the improvement — not just its existence — justifies the cost of shipping a more complex or expensive model. This same logic scales up directly into A/B testing a model change in production (Lesson 69), where "statistically significant" and "worth shipping" are explicitly treated as two separate questions.
10 Probability Calibration: Can You Trust the Score Itself?
Section 1 introduced probability scores and Section 3–4 used them to sweep a threshold and trace an ROC curve — but nothing so far checked whether a predicted probability of 0.8 actually means "80% of the time, this really is the positive class." A model can rank examples perfectly (great AUC) while its raw probability outputs are systematically over- or under-confident. That property — do the numbers mean what they claim to mean — is called calibration, and it's a separate axis of quality from everything measured so far in this lesson.
import numpy as np
from sklearn.calibration import calibration_curve
from sklearn.svm import SVC
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
X, y = make_classification(n_samples=3000, n_features=20, n_informative=10,
weights=[0.5, 0.5], random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
models = {
'Logistic Regression': LogisticRegression(max_iter=1000),
'SVM (decision_function-based)': SVC(probability=True, random_state=42),
}
for name, model in models.items():
model.fit(X_train, y_train)
probs = model.predict_proba(X_test)[:, 1]
# Bin predictions into 10 buckets by predicted probability, compare each
# bucket's AVERAGE predicted probability to its ACTUAL fraction positive
frac_positive, mean_predicted = calibration_curve(y_test, probs, n_bins=10)
print(f"\n{name}:")
for mp, fp in zip(mean_predicted, frac_positive):
gap = fp - mp
flag = " <- overconfident" if gap < -0.05 else (" <- underconfident" if gap > 0.05 else "")
print(f" predicted~{mp:.2f} actual={fp:.2f}{flag}")A perfectly calibrated model's points would sit exactly on the diagonal: among all the times it said "70% confident," the positive class should show up roughly 70% of the time. Some model families are calibrated well out of the box (logistic regression, by construction, since it directly optimizes log loss); others systematically distort their scores — SVMs push probabilities away from the boundary (overconfident), and boosted trees are frequently underconfident near the extremes.
import numpy as np
from sklearn.calibration import CalibratedClassifierCV
from sklearn.svm import SVC
from sklearn.metrics import brier_score_loss
# Wrap an uncalibrated model with CalibratedClassifierCV -- it internally
# cross-validates, fitting a correction (Platt scaling / sigmoid, or
# isotonic regression) mapping raw scores to calibrated probabilities.
raw_svm = SVC(probability=True, random_state=42)
calibrated_svm = CalibratedClassifierCV(raw_svm, method='sigmoid', cv=5)
raw_svm.fit(X_train, y_train)
calibrated_svm.fit(X_train, y_train)
raw_probs = raw_svm.predict_proba(X_test)[:, 1]
calibrated_probs = calibrated_svm.predict_proba(X_test)[:, 1]
# Brier score: mean squared error between predicted probability and the
# actual 0/1 outcome -- lower is better, and unlike AUC it directly
# penalizes miscalibration, not just ranking mistakes.
print(f"Raw SVM Brier score: {brier_score_loss(y_test, raw_probs):.4f}")
print(f"Calibrated SVM Brier score: {brier_score_loss(y_test, calibrated_probs):.4f}")
# Calibration typically leaves AUC nearly unchanged (ranking is preserved)
# while meaningfully improving the Brier score (the numbers become trustworthy).If you only ever threshold a model's output to get a hard label, miscalibration is often harmless — Section 5's cost-benefit threshold search still finds the right cutoff even if the raw numbers are distorted. But the moment you use the probability value itself — displaying "73% risk of default" to a loan officer, feeding a score into a downstream expected-value calculation, or combining probabilities from multiple models — calibration becomes essential. Two of the clearest use cases: risk scoring shown to a human decision-maker, and Lesson 30's GMM responsibilities, which are only meaningful as genuine probabilities in the first place.
Real-World Spotlight: Credit Scoring and Optimal Fraud Thresholds
Banks do not use 0.5 as their fraud detection threshold. They compute the expected cost at every point on the ROC curve and find the threshold that minimises total expected losses. The cost calculation considers: cost of fraud investigation (FP), value of fraud prevented (TP), and cost of chargebacks (FN). Different banks make different tradeoffs based on their customer service budget and fraud loss tolerance.
The entire pipeline is then validated with Stratified K-Fold cross-validation using roc_auc as the evaluation metric (because AUC is threshold-independent and gives a complete view of the model's discrimination ability before the threshold decision). The threshold is chosen separately from model selection.
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import roc_curve, roc_auc_score
from sklearn.datasets import make_classification
import numpy as np
import matplotlib.pyplot as plt
# Simulate credit card transactions (3% fraud rate)
np.random.seed(42)
X, y = make_classification(
n_samples=20000, weights=[0.97, 0.03],
n_features=30, n_informative=15, n_redundant=5,
random_state=42)
# Step 1: Model selection via StratifiedKFold + AUC
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
lr_pipe = Pipeline([('s', StandardScaler()), ('lr', LogisticRegression(max_iter=2000))])
gbc_pipe = Pipeline([('s', StandardScaler()),
('gbc', GradientBoostingClassifier(n_estimators=100, random_state=42))])
lr_auc = cross_val_score(lr_pipe, X, y, cv=skf, scoring='roc_auc')
gbc_auc = cross_val_score(gbc_pipe, X, y, cv=skf, scoring='roc_auc')
print("=== Model Selection (5-Fold Stratified CV) ===")
print(f"Logistic Regression AUC: {lr_auc.mean():.4f} ± {lr_auc.std():.4f}")
print(f"Gradient Boosting AUC: {gbc_auc.mean():.4f} ± {gbc_auc.std():.4f}")
print(f"Winner: {'GradientBoosting' if gbc_auc.mean() > lr_auc.mean() else 'LogisticRegression'}")
# Step 2: Train best model on 80% holdout, evaluate on 20%
from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, stratify=y, random_state=42)
best_pipe = gbc_pipe if gbc_auc.mean() > lr_auc.mean() else lr_pipe
best_pipe.fit(X_tr, y_tr)
y_scores = best_pipe.predict_proba(X_te)[:, 1]
fpr_arr, tpr_arr, thresholds = roc_curve(y_te, y_scores)
final_auc = roc_auc_score(y_te, y_scores)
print(f"\nHoldout test AUC: {final_auc:.4f}")
# Step 3: Cost-benefit threshold selection
n_pos = y_te.sum()
n_neg = len(y_te) - n_pos
cost_fp = 50 # manual review cost
cost_fn = 500 # fraud loss
costs = [(1-tpr)*n_pos*cost_fn + fpr*n_neg*cost_fp
for tpr, fpr in zip(tpr_arr, fpr_arr)]
best_i = np.argmin(costs)
print(f"\n=== Optimal Threshold (cost-benefit) ===")
print(f"Threshold: {thresholds[best_i]:.4f}")
print(f"TPR (recall): {tpr_arr[best_i]:.4f}")
print(f"FPR: {fpr_arr[best_i]:.4f}")
print(f"Expected cost: ${costs[best_i]:,.0f}")
print(f"vs default 0.5: ${costs[np.argmin(np.abs(thresholds - 0.5))]:,.0f}")
print(f"Annual savings (if 20K transactions/day): ${(costs[np.argmin(np.abs(thresholds-0.5))] - costs[best_i])*365:,.0f}")Quick Check
✍️ Practice Exercises
- Train a Logistic Regression and a KNN classifier on the breast cancer dataset. Plot both ROC curves on the same axes. Report the AUC for each and identify which model has better discriminative ability.
- Using a severely imbalanced dataset (
make_classification(weights=[0.98, 0.02])), demonstrate that the ROC curve can look optimistic while the Precision-Recall curve reveals the model's weakness. Plot both side-by-side. - Simulate the cost-benefit analysis: on any binary classification dataset, set FP_cost=100 and FN_cost=200. Plot expected cost vs threshold. Find the optimal threshold. Then repeat with FN_cost=1000. How does the optimal threshold shift?
- Demonstrate data leakage: on the breast cancer dataset, scale before cross-validation and then inside a Pipeline. Report the accuracy for both approaches and the magnitude of the leakage inflation.
- Use
cross_validatewithStratifiedKFold(n_splits=10)to evaluate a LogisticRegression on load_digits across 5 metrics simultaneously (accuracy, roc_auc with ovr, f1_macro, precision_macro, recall_macro). Report mean ± std for each.
▶ cross_validate multi-metric hint
from sklearn.model_selection import cross_validate, StratifiedKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_digits
X, y = load_digits(return_X_y=True)
pipe = Pipeline([('s', StandardScaler()),
('lr', LogisticRegression(max_iter=5000, multi_class='ovr'))])
skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
scoring = {
'accuracy': 'accuracy',
'roc_auc_ovr': 'roc_auc_ovr',
'f1_macro': 'f1_macro',
'precision_macro': 'precision_macro',
'recall_macro': 'recall_macro',
}
results = cross_validate(pipe, X, y, cv=skf, scoring=scoring)
for metric in scoring:
m = results[f'test_{metric}']
print(f"{metric:20} {m.mean():.4f} ± {m.std():.4f}")📚 Primary Sources
sklearn: ROC and AUC — official documentation with worked examples and multi-class ROC.
sklearn: Cross-Validation — comprehensive guide to all CV strategies and pitfalls.