Refinement trials and the paired test

import lensml
import numpy as np
import pandas as pd
from lensml import (
    FEATURES,
    MAX_TESTED_TRIALS,
    TRIALS,
    Register,
    at_base_rate,
    recall_vs_calls,
    setup,
    show,
    show_p,
    verdict,
)
from scipy.stats import t as t_dist

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_colwidth", 70)

S = setup()  # the cache, the test members, the transforms (chapters 1-2)
train, test = S.train, S.test
X_tr, y_tr = train[FEATURES].values, train["wet"].values
g_tr = train["member_id"].values
print(
    f"train: {S.n_train_members} members / {len(train)} winters   "
    f"test: {test.member_id.nunique()} members / {len(test)} winters   "
    f"N_FOLDS = {S.n_folds}"
)

# Adapted from scikit-learn's example "Statistical comparison of models
# using grid search" (examples/model_selection/plot_grid_search_stats.py).
# Copyright (c) 2007-2026 The scikit-learn developers. BSD-3-Clause; full
# text in LICENSES/scikit-learn-BSD-3-Clause.txt.


def corrected_std(differences, n_train, n_test):
    """Corrects standard deviation using Nadeau and Bengio's approach."""
    kr = len(differences)
    corrected_var = np.var(differences, ddof=1) * (1 / kr + n_test / n_train)
    return np.sqrt(corrected_var)


def compute_corrected_ttest(differences, df, n_train, n_test):
    """Computes right-tailed paired t-test with corrected variance."""
    mean = np.mean(differences)
    std = corrected_std(differences, n_train, n_test)
    if std == 0:  # added: identical scores in every fold (guards 0/0)
        return 0.0, 1.0
    t_stat = mean / std
    # changed from the example: t_stat, not np.abs(t_stat), so the test is
    # right-tailed as the docstring says and a LOSS prints p near 1
    p_val = t_dist.sf(t_stat, df)
    return t_stat, p_val

show(verdict, Register.record)

reg = Register(S)
for name, make_model, what_changed in TRIALS[:3]:  # rungs 1-3
    reg.record(name, make_model, what_changed)
print(
    f"incumbent after the ladder ({len(reg.rows)} trials): {reg.incumbent_name}"
)
reg.bonferroni_line()
reg.table()

show(at_base_rate)
reg.record(*TRIALS[3])  # logistic@base-rate
print(f"incumbent after {len(reg.rows)} trials: {reg.incumbent_name}")
reg.bonferroni_line()
reg.table().tail(1)

reg.record(*TRIALS[4])  # boosting@base-rate
print(f"incumbent after {len(reg.rows)} trials: {reg.incumbent_name}")
reg.bonferroni_line()
recall_vs_calls(reg, "boosting@base-rate")
reg.table().tail(1)

YOUR_TRIALS = [
    # ('tree@base-rate', at_base_rate('tree'), 'tree, threshold lowered to the base rate'),
]
for t in YOUR_TRIALS:
    reg.record(*t)  # a name already recorded is skipped
left = MAX_TESTED_TRIALS - len(reg.rows)
print(
    f"incumbent after {len(reg.rows)} trials: {reg.incumbent_name}   "
    f"(tested rows left: {left})"
)
reg.bonferroni_line()
reg.table()

promoted = [r for r in reg.rows if r["verdict"].startswith("promoted")]
last = promoted[-1] if promoted else reg.rows[-1]
PAIR_A, PAIR_B = (
    last["incumbent"],
    last["candidate"],
)  # any two names in the register
promoted_on = (
    last["verdict"].split(":")[1].split()[0]
    if last["verdict"].startswith("promoted")
    else None
)
HAND_METRIC = promoted_on or "precision"
print(
    f"default pair: register row {last['trial']} ({last['candidate']} vs "
    f"{last['incumbent']}, verdict '{last['verdict']}') → metric {HAND_METRIC}"
    + ("" if promoted_on else "  (nothing promoted: precision by default)")
)

n_test = len(train) // S.n_folds
n_train = len(train) - n_test  # exactly as compare() sets them
a, b = reg.scores[PAIR_A][HAND_METRIC], reg.scores[PAIR_B][HAND_METRIC]
diff = b - a  # one paired difference per fold: candidate minus incumbent
k = len(diff)
df_hand = k - 1
sd_hand = diff.std(ddof=1)  # the plain per-fold SD
se_hand = np.sqrt(sd_hand**2 * (1 / k + n_test / n_train))  # corrected SE
if se_hand == 0:  # identical scores in every fold
    t_hand, p_hand = 0.0, 1.0
else:
    t_hand = diff.mean() / se_hand
    p_hand = t_dist.sf(t_hand, df_hand)  # one-sided: the RIGHT tail
p_loss = t_dist.cdf(t_hand, df_hand)  # the LEFT tail: small when PAIR_B lost

# the hand arithmetic == the functions above == lensml's, on this pair
assert np.isclose(se_hand, corrected_std(diff, n_train, n_test))
ref = lensml.compare(
    reg.scores[PAIR_A], reg.scores[PAIR_B], len(train), S.n_folds, HAND_METRIC
)
assert np.allclose(
    [t_hand, p_hand],
    compute_corrected_ttest(diff, df_hand, n_train, n_test),
)
assert np.allclose(
    [t_hand, p_hand], [ref[f"{HAND_METRIC}_t"], ref[f"{HAND_METRIC}_p"]]
)
for r in reg.rows:  # ... and on every register row, so the two cannot drift
    inc, cand = reg.scores[r["incumbent"]], reg.scores[r["candidate"]]
    for m in ("recall", "precision"):
        d = cand[m] - inc[m]
        ref = lensml.compare(inc, cand, len(train), S.n_folds, m)
        assert np.allclose(
            compute_corrected_ttest(d, len(d) - 1, n_train, n_test),
            [ref[f"{m}_t"], ref[f"{m}_p"]],
        )

print(
    f"{PAIR_B} − {PAIR_A} on {HAND_METRIC}: {len(diff)} paired fold scores   "
    f"n_train = {n_train}   n_test = {n_test}"
)
print(
    f"means {b.mean():.3f} vs {a.mean():.3f}   gap {diff.mean():+.4f}   per-fold SD of "
    f"the gap {sd_hand:.4f}   corrected standard error "
    f"{se_hand:.4f}  (the spread the test uses)"
)
if p_hand < 0.05:
    reading = (
        f"p < 0.05: a measured GAIN for {PAIR_B} under the uncorrected rule; "
        "read it against the Bonferroni line below"
    )
elif p_loss < 0.05:
    reading = (
        f"a measured LOSS for {PAIR_B}: {abs(t_hand):.1f} corrected standard errors "
        f"(left-tail p = {show_p(p_loss)})"
    )
else:
    reading = "NOT a measured difference (neither tail under 0.05)"
print(f"t = {t_hand:.2f}   p = {show_p(p_hand)}   →  {reading}")
reg.bonferroni_line()