Chapter 20: Horse Health: Small Data, Big Discipline

import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold

from kaggle_companion.activities._common import clean

SEED = 20
REPLICATES = 14      # independent constructed datasets per development size
N_NEW = 4000         # fresh rows from the same generator, never used for fitting or tuning
GRID = np.exp(np.linspace(np.log(0.5), np.log(2.0), 13))
FACTORS = np.array([(1.0, a, b) for a in GRID for b in GRID])   # divide class 0, 1, 2 probabilities by these, then take argmax


def generate(n, rng, weights, bias):
    """Three-class task with moderate imbalance (about 45% / 33% / 22%) and a few nonlinear terms."""
    X = rng.normal(size=(n, 8))
    z = X @ weights + bias + 0.5 * np.column_stack([np.sin(2 * X[:, 0]), X[:, 1] * X[:, 2], np.zeros(n)])
    p = np.exp(z - z.max(1, keepdims=True))
    p /= p.sum(1, keepdims=True)
    return X, np.minimum((rng.random(n)[:, None] > p.cumsum(1)).sum(1), 2)


def model():
    return HistGradientBoostingClassifier(learning_rate=0.1, max_leaf_nodes=6, max_iter=40, early_stopping=False, random_state=0)


def tune(proba, y):
    """Grid-search the factors that maximize accuracy on these rows; return them and that accuracy."""
    accuracy = (np.argmax(proba[:, None, :] / FACTORS[None], 2) == y[:, None]).mean(0)
    best = int(np.argmax(accuracy))
    return FACTORS[best], accuracy[best]


def gain(proba, y, factors):
    """Accuracy of factor-adjusted argmax minus accuracy of plain argmax, on these rows."""
    return float((np.argmax(proba / factors, 1) == y).mean() - (proba.argmax(1) == y).mean())


def one_dataset(n, seed):
    rng = np.random.default_rng(seed)
    weights, bias = rng.normal(0, 0.55, (8, 3)), np.array([0.5, 0.0, -0.5])
    X, y = generate(n, rng, weights, bias)
    X_new, y_new = generate(N_NEW, rng, weights, bias)
    oof = np.zeros((n, 3))
    for fit, out in StratifiedKFold(5, shuffle=True, random_state=1).split(X, y):
        oof[out] = model().fit(X[fit], y[fit]).predict_proba(X[out])
    new = model().fit(X, y).predict_proba(X_new)         # the deployed model, refitted on every development row
    prior = np.bincount(y, minlength=3) / n
    half = rng.permutation(n)
    halves = (half[: n // 2], half[n // 2:])
    out = {}
    # "plain": the model's own probabilities. "rebalanced": divided by training class frequency and renormalized,
    # which distorts them toward the rare classes in the way class-weighted training does.
    for name, (p_oof, p_new) in {"plain": (oof, new), "rebalanced": (oof / prior, new / prior)}.items():
        p_oof, p_new = p_oof / p_oof.sum(1, keepdims=True), p_new / p_new.sum(1, keepdims=True)
        factors, tuned_accuracy = tune(p_oof, y)
        # Same rows: tuned on the OOF rows and scored on those same rows. Cross-half: tuned on one half, scored on the other.
        same_rows = tuned_accuracy - (p_oof.argmax(1) == y).mean()
        cross = np.mean([gain(p_oof[b], y[b], tune(p_oof[a], y[a])[0]) for a, b in (halves, halves[::-1])])
        out[name] = (same_rows, cross, gain(p_new, y_new, factors), (p_oof.argmax(1) == y).mean())
    return out


def run(development_rows):
    sets = [one_dataset(development_rows, SEED * 1000 + i) for i in range(REPLICATES)]
    models = {}
    for name in ("plain", "rebalanced"):
        v = np.array([s[name] for s in sets])
        models[name] = {"same_rows": float(v[:, 0].mean()), "cross_half": float(v[:, 1].mean()), "new_rows": float(v[:, 2].mean()),
                        "new_rows_se": float(v[:, 2].std() / np.sqrt(REPLICATES)), "argmax_accuracy": float(v[:, 3].mean())}
    return clean({"development_rows": development_rows, "new_rows": N_NEW, "replicates": REPLICATES, "models": models})

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

VALUES = [250, 500, 988, 3000]
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.ch20 import draw, NCOLS, HEIGHT

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