Chapter 38: Calibration and Post-Processing

import numpy as np
from sklearn.isotonic import IsotonicRegression
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss, roc_auc_score
from sklearn.naive_bayes import GaussianNB

from kaggle_companion.activities._common import clean

SEED = 38
DRAWS = 30        # independent development sets per size; results are means over draws
EPS = 1e-3        # probabilities are kept inside [0.001, 0.999] so one confident miss cannot dominate log loss


def generate(n, rng):
    """Six noisy copies of one hidden signal plus two noise features. The true log odds bend upward for positive signal."""
    z = rng.normal(size=n)
    X = np.column_stack([z + rng.normal(0, 0.6, n) for _ in range(6)] + [rng.normal(size=n), rng.normal(size=n)])
    logit = 0.4 * z + 1.6 * np.maximum(z, 0) - 1.0
    y = (rng.random(n) < 1 / (1 + np.exp(-logit))).astype(int)
    return X, y


def to_logit(p):
    p = np.clip(p, EPS, 1 - EPS)
    return np.log(p / (1 - p))


def expected_calibration_error(y, p, bins=10):
    """The chapter's population-weighted ECE with equal-width bins; empty bins contribute nothing."""
    b = np.minimum((p * bins).astype(int), bins - 1)
    return float(sum((b == k).mean() * abs(y[b == k].mean() - p[b == k].mean()) for k in range(bins) if (b == k).any()))


def reliability(y, p, bins=10):
    b = np.minimum((p * bins).astype(int), bins - 1)
    return [[float(p[b == k].mean()), float(y[b == k].mean()), int((b == k).sum())] for k in range(bins) if (b == k).sum() >= 20]


def run(development_rows):
    rng = np.random.default_rng(SEED)
    X_base, y_base = generate(3000, rng)
    base = GaussianNB().fit(X_base, y_base)                    # frozen base model: never sees development or assessment labels
    X_a, y_a = generate(20000, rng)                            # assessment rows, used only for scoring
    p_a = np.clip(base.predict_proba(X_a)[:, 1], EPS, 1 - EPS)

    scores = {k: [] for k in ("none", "platt", "isotonic", "none_dev", "platt_dev", "isotonic_dev",
                              "auc_platt", "auc_isotonic")}
    for draw in range(DRAWS):
        X_d, y_d = generate(development_rows, rng)
        p_d = np.clip(base.predict_proba(X_d)[:, 1], EPS, 1 - EPS)
        platt = LogisticRegression(C=1e6).fit(to_logit(p_d)[:, None], y_d)
        iso = IsotonicRegression(out_of_bounds="clip", y_min=EPS, y_max=1 - EPS).fit(p_d, y_d)
        cal = {"platt": lambda p: platt.predict_proba(to_logit(p)[:, None])[:, 1], "isotonic": iso.predict}
        scores["none"].append(log_loss(y_a, p_a))
        scores["none_dev"].append(log_loss(y_d, p_d, labels=[0, 1]))
        for name, f in cal.items():
            scores[name].append(log_loss(y_a, f(p_a)))
            scores[name + "_dev"].append(log_loss(y_d, f(p_d), labels=[0, 1]))   # scored on the rows it was fitted to
            scores["auc_" + name].append(roc_auc_score(y_a, f(p_a)))
        if draw == 0:
            curves = {"none": reliability(y_a, p_a), "platt": reliability(y_a, cal["platt"](p_a)),
                      "isotonic": reliability(y_a, cal["isotonic"](p_a))}
            ece = {"none": expected_calibration_error(y_a, p_a),
                   "platt": expected_calibration_error(y_a, cal["platt"](p_a)),
                   "isotonic": expected_calibration_error(y_a, cal["isotonic"](p_a))}

    mean = {k: float(np.mean(v)) for k, v in scores.items()}
    return clean({
        "development_rows": development_rows,
        "assessment": {k: mean[k] for k in ("none", "platt", "isotonic")},
        "development": {k: mean[k + "_dev"] for k in ("none", "platt", "isotonic")},
        "isotonic_beats_platt": float(np.mean(np.array(scores["isotonic"]) < np.array(scores["platt"]))),
        "auc": {"none": roc_auc_score(y_a, p_a), "platt": mean["auc_platt"], "isotonic": mean["auc_isotonic"]},
        "ece_first_draw": ece, "reliability_first_draw": curves,
        "draws": DRAWS, "assessment_rows": len(y_a),
    })

import pandas as pd
from kaggle_companion.activities.ch38 import explain

VALUES = [50, 200, 1000, 5000]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch38 import draw, NCOLS, HEIGHT

DEFAULT = 200
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()