Chapter 16: The Fourth-Root Blend

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss

from kaggle_companion.activities._common import clean

SEED = 16
REPLICATES = 40                      # independent constructed datasets; every estimate is their mean
NOISE = np.array([1.0, 1.15, 1.3])   # how noisy each model's view of the hidden signal is: three similar-strength models
N_TRAIN, N_DEV, N_ASSESS = 600, 300, 20000
EPS = 1e-4


def generate(rng, n, shared):
    """A hidden signal drives the label. Each model sees a noisy copy; `shared` is the share of each copy's noise
    variance that is common to all three models, so it sets how alike their errors are."""
    signal = rng.normal(size=n)
    y = (rng.random(n) < 1 / (1 + np.exp(-1.6 * signal))).astype(int)
    common = rng.normal(size=n)
    views = [signal + sd * (np.sqrt(shared) * common + np.sqrt(1 - shared) * rng.normal(size=n)) for sd in NOISE]
    return np.column_stack(views), y


def fourth_root_weights(scores, baseline, maximize=True):
    """The chapter's heuristic: compress each positive margin over a measured baseline with a fourth root."""
    scores = np.asarray(scores, dtype=float)
    margin = scores - baseline if maximize else baseline - scores
    raw = np.maximum(margin, 0) ** 0.25
    return raw / raw.sum()


def greedy_weights(y, P, step=0.1, max_iterations=30):
    """Hill-climbing on development log loss: start from the best model, move weight to whichever model helps most."""
    loss = lambda q: log_loss(y, np.clip(q, EPS, 1 - EPS))
    weights = np.eye(P.shape[1])[int(np.argmin([loss(P[:, j]) for j in range(P.shape[1])]))]
    current = P @ weights
    for _ in range(max_iterations):
        trial = [loss((1 - step) * current + step * P[:, j]) for j in range(P.shape[1])]
        j = int(np.argmin(trial))
        if trial[j] >= loss(current) - 1e-12:
            break
        weights = weights * (1 - step)
        weights[j] += step
        current = P @ weights
    return weights


def run(shared):
    loss = lambda y, q: log_loss(y, np.clip(q, EPS, 1 - EPS))
    rows = []
    for rep in range(REPLICATES):
        rng = np.random.default_rng(SEED * 100 + rep)
        (X, y), (Xd, yd), (Xa, ya) = [generate(rng, n, shared) for n in (N_TRAIN, N_DEV, N_ASSESS)]
        models = [LogisticRegression().fit(X[:, [m]], y) for m in range(3)]      # one model per view
        Pd = np.column_stack([models[m].predict_proba(Xd[:, [m]])[:, 1] for m in range(3)])
        Pa = np.column_stack([models[m].predict_proba(Xa[:, [m]])[:, 1] for m in range(3)])
        dev_loss = np.array([loss(yd, Pd[:, m]) for m in range(3)])
        dummy = loss(yd, np.full(N_DEV, y.mean()))                              # measured baseline: a constant-prior predictor
        w_root = fourth_root_weights(dev_loss, dummy, maximize=False)
        w_greedy = greedy_weights(yd, Pd)
        logits = np.log(Pd / (1 - Pd))
        corr = np.corrcoef(logits.T)
        rows.append({
            "best": loss(ya, Pa[:, int(np.argmin(dev_loss))]),                  # strongest model, chosen on development rows
            "equal": loss(ya, Pa.mean(axis=1)), "root": loss(ya, Pa @ w_root), "greedy": loss(ya, Pa @ w_greedy),
            "equal_dev": loss(yd, Pd.mean(axis=1)), "greedy_dev": loss(yd, Pd @ w_greedy),
            "correlation": (corr[0, 1] + corr[0, 2] + corr[1, 2]) / 3, "root_spread": w_root.max() - w_root.min(),
            "w_root": w_root, "w_greedy": w_greedy,
        })
    col = lambda k: np.array([r[k] for r in rows])
    gain = {k: col("best") - col(k) for k in ("equal", "root", "greedy")}      # positive: the blend beats the best single model
    return clean({
        "shared_noise": shared, "correlation": float(col("correlation").mean()), "root_weight_spread": float(col("root_spread").mean()),
        "assessment_loss": {k: float(col(k).mean()) for k in ("best", "equal", "root", "greedy")},
        "gain_over_best": {k: float(v.mean()) for k, v in gain.items()},
        "gain_standard_error": {k: float(v.std() / np.sqrt(REPLICATES)) for k, v in gain.items()},
        "share_equal_beats_best": float((gain["equal"] > 0).mean()),
        "mean_weights": {"root": np.mean([r["w_root"] for r in rows], axis=0), "greedy": np.mean([r["w_greedy"] for r in rows], axis=0)},
        "greedy_vs_equal": {"development": float((col("equal_dev") - col("greedy_dev")).mean()),
                            "assessment": float((col("equal") - col("greedy")).mean())},
        "replicates": REPLICATES, "assessment_rows": N_ASSESS, "development_rows": N_DEV,
    })

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

VALUES = [0.0, 0.3, 0.7, 0.95]
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.ch16 import draw, NCOLS, HEIGHT

DEFAULT = 0.95
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()