hack-house/tests/test_sor_stats.py

"""Pre-registered §6 statistics — calibrated on synthetic ground truth only.

These tests validate the inference toolkit (BCa bootstrap CIs, Miller-Madow
entropy, Spearman, Holm) against distributions whose answer is known by
construction — never against confirmatory-cell data (there is none). They also
pin the §6 reproducibility contract: a fixed resample seed yields an identical
CI, and the three-seed spot-check agrees to Monte-Carlo error.
"""

import math

from cmd_chat.sor.analysis.detectors import auc, shannon_entropy_bits
from cmd_chat.sor.analysis.stats import (
    bootstrap_ci,
    holm_bonferroni,
    miller_madow_entropy_bits,
    spearman,
    two_sample_diff_ci,
)


# --------------------------------------------------------------------------- #
# Bootstrap CI — one sample.
# --------------------------------------------------------------------------- #
def _labelled_scores(sep: float, n: int = 60):
    """n linked ('pos') and n unlinked ('neg') scores separated by `sep`, from a
    fixed LCG so the fixture is deterministic. Returns units = (score, is_pos)."""
    state = 0x1234_5678
    def nxt():
        nonlocal state
        state = (1103515245 * state + 12345) & 0x7FFFFFFF
        return state / 0x7FFFFFFF
    units = []
    for _ in range(n):
        units.append((sep + nxt(), True))
        units.append((nxt(), False))
    return units


def _auc_stat(units):
    pos = [s for s, is_pos in units if is_pos]
    neg = [s for s, is_pos in units if not is_pos]
    return auc(pos, neg)


def test_bootstrap_ci_excludes_half_for_separated_scores():
    units = _labelled_scores(sep=0.9)
    ci = bootstrap_ci(units, _auc_stat, n_resamples=2000, seed=1)
    assert ci.point > 0.6
    assert ci.excludes(0.5) and ci.strictly_greater(0.5)


def test_bootstrap_ci_includes_half_for_overlapping_scores():
    units = _labelled_scores(sep=0.0)  # pos/neg drawn from the same distribution
    ci = bootstrap_ci(units, _auc_stat, n_resamples=2000, seed=1)
    assert not ci.excludes(0.5)


def test_bootstrap_ci_is_reproducible_from_seed():
    units = _labelled_scores(sep=0.7)
    a = bootstrap_ci(units, _auc_stat, n_resamples=1500, seed=42)
    b = bootstrap_ci(units, _auc_stat, n_resamples=1500, seed=42)
    assert (a.lo, a.hi, a.point, a.method) == (b.lo, b.hi, b.point, b.method)


def test_bootstrap_three_seed_spot_check_agrees_to_mc_error():
    units = _labelled_scores(sep=0.8)
    cis = [bootstrap_ci(units, _auc_stat, n_resamples=2000, seed=s) for s in (1, 2, 3)]
    los = [c.lo for c in cis]
    his = [c.hi for c in cis]
    assert max(los) - min(los) < 0.05
    assert max(his) - min(his) < 0.05


# --------------------------------------------------------------------------- #
# Two-sample difference CI (ΔH shape).
# --------------------------------------------------------------------------- #
def test_two_sample_diff_ci_detects_positive_shift():
    a = [5.0 + (i % 3) * 0.1 for i in range(40)]
    b = [3.0 + (i % 3) * 0.1 for i in range(40)]
    ci = two_sample_diff_ci(a, b, lambda xs: sum(xs) / len(xs),
                            n_resamples=2000, seed=7)
    assert ci.point > 1.8
    assert ci.strictly_greater(0.0)


def test_two_sample_diff_ci_spans_zero_for_equal_arms():
    a = [1.0 + (i % 5) * 0.2 for i in range(40)]
    b = [1.0 + (i % 5) * 0.2 for i in range(40)]
    ci = two_sample_diff_ci(a, b, lambda xs: sum(xs) / len(xs),
                            n_resamples=2000, seed=7)
    assert not ci.excludes(0.0)


# --------------------------------------------------------------------------- #
# Miller-Madow entropy.
# --------------------------------------------------------------------------- #
def test_miller_madow_equals_log2_n_at_the_limit_and_corrects_upward():
    counts = [100] * 8  # 8 equiprobable senders
    plug = shannon_entropy_bits(counts)
    mm = miller_madow_entropy_bits(counts)
    assert math.isclose(plug, 3.0, abs_tol=1e-9)  # log2(8)
    assert mm > plug  # bias correction adds (K-1)/(2 N ln2)
    assert mm - plug < 0.02  # small at N=800


def test_miller_madow_zero_for_empty():
    assert miller_madow_entropy_bits([]) == 0.0
    assert miller_madow_entropy_bits([0, 0]) == 0.0


# --------------------------------------------------------------------------- #
# Spearman.
# --------------------------------------------------------------------------- #
def test_spearman_monotone_and_antitone():
    xs = [1, 2, 3, 4, 5, 6]
    assert math.isclose(spearman(xs, [2, 4, 6, 8, 10, 12]), 1.0, abs_tol=1e-9)
    assert math.isclose(spearman(xs, [12, 10, 8, 6, 4, 2]), -1.0, abs_tol=1e-9)


def test_spearman_degenerate_inputs():
    assert spearman([], []) == 0.0
    assert spearman([1, 2, 3], [5, 5, 5]) == 0.0  # zero variance


# --------------------------------------------------------------------------- #
# Holm-Bonferroni with a frozen family larger than the reported subset.
# --------------------------------------------------------------------------- #
def test_holm_uses_full_family_size_not_reported_subset():
    # 4 reported tests embedded in a frozen family of 7: the k-th smallest uses
    # multiplier (7 - k + 1) = 7,6,5,4 rather than 4,3,2,1.
    ps = {"RQ1-P1": 0.001, "RQ1-P2": 0.004, "RQ2-P1": 0.02, "RQ2-P3": 0.30}
    res = holm_bonferroni(ps, alpha=0.05, family_size=7)
    by_name = {r.name: r for r in res}
    assert by_name["RQ1-P1"].multiplier == 7
    assert by_name["RQ1-P2"].multiplier == 6
    assert by_name["RQ2-P1"].multiplier == 5
    assert by_name["RQ2-P3"].multiplier == 4
    # 0.001*7 = 0.007 rejected; 0.30*4 = 1.0 not.
    assert by_name["RQ1-P1"].reject
    assert not by_name["RQ2-P3"].reject


def test_holm_is_more_conservative_than_reported_only():
    ps = {"a": 0.01, "b": 0.02}
    full = {r.name: r.p_adjusted for r in holm_bonferroni(ps, family_size=7)}
    subset = {r.name: r.p_adjusted for r in holm_bonferroni(ps, family_size=2)}
    assert full["a"] > subset["a"]  # larger denominator -> larger adjusted p


def test_holm_adjusted_p_is_monotone_and_capped():
    ps = {"a": 0.01, "b": 0.2, "c": 0.9}
    res = holm_bonferroni(ps, family_size=3)
    adj = [r.p_adjusted for r in res]
    assert adj == sorted(adj)  # non-decreasing in rank
    assert all(p <= 1.0 for p in adj)


def test_holm_rejects_family_size_smaller_than_reported():
    try:
        holm_bonferroni({"a": 0.1, "b": 0.2}, family_size=1)
    except ValueError:
        return
    raise AssertionError("expected ValueError for family_size < reported count")