R7/analysis: pre-registered §6 inference toolkit (BCa bootstrap, Miller-Madow, Spearman, Holm)

Turnkey implementation of the frozen prereg §6 analysis plan, written before any
confirmatory data exists (analysis-precedes-data, rigor-standards §Statistics).
Pure stdlib, no I/O, no engine, no traffic — calibrated on synthetic ground truth only:

- bootstrap_ci / two_sample_diff_ci: BCa 95% CIs (10k resamples default) with a
  percentile fallback when bias/acceleration terms are degenerate; seeded and
  reproducible (CIResult carries method+seed for the §6 three-seed spot-check).
- miller_madow_entropy_bits: plug-in Shannon entropy + Miller-Madow bias
  correction (§3 estimator).
- spearman: rank correlation for RQ2-P3.
- holm_bonferroni: step-down multiplicity correction with explicit family_size so
  a lead paper reporting a subset of the frozen 7-test family still corrects
  against the full family (never re-optimised to the reported subset).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

cmd_chat/sor/analysis/stats.py

@@ -0,0 +1,326 @@
+"""Pre-registered inferential statistics for the RQ1/RQ2 confirmatory analysis.
+
+This module is the *turnkey* implementation of the frozen prereg §6 analysis plan
+(`sor-consent-prereg.md`), written **before** any confirmatory data exists so the
+analysis precedes the data (`rigor-standards §Statistics`). It computes nothing
+study-specific by itself — it is a small, stdlib-only inference toolkit:
+
+  * :func:`bootstrap_ci` / :func:`two_sample_diff_ci` — BCa (bias-corrected and
+    accelerated) bootstrap 95% CIs, 10,000 resamples by default, with a percentile
+    fallback when the acceleration/bias terms are degenerate (§6: "Bootstrap:
+    10,000 resamples, BCa intervals; seed spot-check");
+  * :func:`miller_madow_entropy_bits` — plug-in Shannon entropy with the
+    Miller–Madow finite-sample bias correction (§3 estimator [APPROVAL]);
+  * :func:`spearman` — Spearman rank correlation (RQ2-P3 mechanism test);
+  * :func:`holm_bonferroni` — Holm step-down multiplicity correction over the
+    frozen confirmatory family, with an explicit ``family_size`` so a lead paper
+    that reports a subset of the 7 pre-registered tests still corrects against the
+    full family (§6; never re-optimised to the reported subset).
+
+It performs no I/O, moves no traffic, and spawns no engine. Resampling is seeded
+(``random.Random(seed)``) and the seed is returned in every :class:`CIResult`, so
+a CI is reproducible and the §6 three-seed spot-check is mechanical.
+"""
+
+from __future__ import annotations
+
+import math
+import random
+from dataclasses import dataclass
+from statistics import NormalDist
+from typing import Callable, List, Sequence, TypeVar
+
+T = TypeVar("T")
+
+_NORM = NormalDist(0.0, 1.0)
+DEFAULT_RESAMPLES = 10_000
+DEFAULT_ALPHA = 0.05
+
+
+@dataclass(frozen=True)
+class CIResult:
+    """A point estimate with a bootstrap confidence interval and full provenance
+    (method actually used, resample count, seed, alpha) so it is reproducible and
+    auditable. ``excludes(v)`` is the pre-registered gate primitive: True iff the
+    whole interval lies on one side of ``v`` (e.g. RQ1-P1's "CI excludes 0.5")."""
+
+    point: float
+    lo: float
+    hi: float
+    alpha: float
+    n_resamples: int
+    method: str  # "bca" | "percentile"
+    seed: int
+
+    def excludes(self, v: float) -> bool:
+        return (self.lo > v) or (self.hi < v)
+
+    def strictly_greater(self, v: float) -> bool:
+        return self.lo > v
+
+    def strictly_less(self, v: float) -> bool:
+        return self.hi < v
+
+    def as_dict(self) -> dict:
+        return {
+            "point": self.point,
+            "ci_lo": self.lo,
+            "ci_hi": self.hi,
+            "alpha": self.alpha,
+            "n_resamples": self.n_resamples,
+            "method": self.method,
+            "seed": self.seed,
+        }
+
+
+def _percentile(sorted_vals: Sequence[float], q: float) -> float:
+    """Linear-interpolation percentile of an already-sorted sequence, ``q`` in
+    [0, 1]. Empty -> NaN; clamps out-of-range q to the endpoints."""
+    n = len(sorted_vals)
+    if n == 0:
+        return float("nan")
+    if q <= 0:
+        return float(sorted_vals[0])
+    if q >= 1:
+        return float(sorted_vals[-1])
+    pos = q * (n - 1)
+    lo = int(math.floor(pos))
+    hi = int(math.ceil(pos))
+    if lo == hi:
+        return float(sorted_vals[lo])
+    frac = pos - lo
+    return float(sorted_vals[lo]) * (1 - frac) + float(sorted_vals[hi]) * frac
+
+
+def _bca_endpoints(
+    thetas: Sequence[float], theta_hat: float, jack: Sequence[float], alpha: float
+) -> tuple[float, float, str]:
+    """Return (q_lo, q_hi, method) adjusted percentiles for a BCa interval, or the
+    plain (alpha/2, 1-alpha/2, "percentile") pair when the bias/acceleration terms
+    are degenerate (all resamples equal, or zero jackknife spread)."""
+    b = len(thetas)
+    n_less = sum(1 for t in thetas if t < theta_hat)
+    # Bias correction z0. If every resample is on one side, z0 is undefined ->
+    # fall back to the percentile interval rather than emit a garbage bound.
+    if n_less == 0 or n_less == b:
+        return alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
+    z0 = _NORM.inv_cdf(n_less / b)
+
+    # Acceleration from the jackknife leave-one-out estimates.
+    jbar = sum(jack) / len(jack)
+    diffs = [jbar - j for j in jack]
+    denom = 6.0 * (sum(d * d for d in diffs) ** 1.5)
+    if denom == 0.0:
+        return alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
+    a = sum(d ** 3 for d in diffs) / denom
+
+    z_lo = _NORM.inv_cdf(alpha / 2.0)
+    z_hi = _NORM.inv_cdf(1.0 - alpha / 2.0)
+
+    def adjust(z: float) -> float:
+        num = z0 + z
+        return _NORM.cdf(z0 + num / (1.0 - a * num))
+
+    q_lo = adjust(z_lo)
+    q_hi = adjust(z_hi)
+    if not (0.0 < q_lo < q_hi < 1.0):
+        return alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
+    return q_lo, q_hi, "bca"
+
+
+def bootstrap_ci(
+    units: Sequence[T],
+    statistic: Callable[[Sequence[T]], float],
+    *,
+    n_resamples: int = DEFAULT_RESAMPLES,
+    alpha: float = DEFAULT_ALPHA,
+    seed: int = 0,
+    method: str = "bca",
+) -> CIResult:
+    """One-sample bootstrap CI of ``statistic`` over ``units`` (the unit of
+    analysis — a circuit-pair for RQ1, a circuit for RQ2). Resamples ``units`` with
+    replacement ``n_resamples`` times. ``method="bca"`` applies bias-correction +
+    acceleration (falling back to percentile if degenerate); ``"percentile"``
+    forces the plain interval."""
+    n = len(units)
+    if n == 0:
+        return CIResult(float("nan"), float("nan"), float("nan"),
+                        alpha, n_resamples, "empty", seed)
+    theta_hat = float(statistic(units))
+    rng = random.Random(seed)
+    thetas: List[float] = []
+    for _ in range(n_resamples):
+        sample = [units[rng.randrange(n)] for _ in range(n)]
+        thetas.append(float(statistic(sample)))
+    thetas.sort()
+
+    if method == "bca" and n > 1:
+        jack = [float(statistic([units[j] for j in range(n) if j != i])) for i in range(n)]
+        q_lo, q_hi, used = _bca_endpoints(thetas, theta_hat, jack, alpha)
+    else:
+        q_lo, q_hi, used = alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
+
+    return CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
+                    alpha, n_resamples, used, seed)
+
+
+def two_sample_diff_ci(
+    units_a: Sequence[T],
+    units_b: Sequence[T],
+    statistic: Callable[[Sequence[T]], float],
+    *,
+    n_resamples: int = DEFAULT_RESAMPLES,
+    alpha: float = DEFAULT_ALPHA,
+    seed: int = 0,
+    method: str = "bca",
+) -> CIResult:
+    """Bootstrap CI for the difference ``statistic(A) - statistic(B)`` of two
+    independent arms (RQ2-P1: ΔH = H(federated) − H(single-house, matched N)).
+    Each arm is resampled independently. BCa uses a combined leave-one-out
+    jackknife across both arms (each point dropped from its own arm)."""
+    na, nb = len(units_a), len(units_b)
+    if na == 0 or nb == 0:
+        return CIResult(float("nan"), float("nan"), float("nan"),
+                        alpha, n_resamples, "empty", seed)
+    theta_hat = float(statistic(units_a)) - float(statistic(units_b))
+    rng = random.Random(seed)
+    thetas: List[float] = []
+    for _ in range(n_resamples):
+        sa = [units_a[rng.randrange(na)] for _ in range(na)]
+        sb = [units_b[rng.randrange(nb)] for _ in range(nb)]
+        thetas.append(float(statistic(sa)) - float(statistic(sb)))
+    thetas.sort()
+
+    if method == "bca" and na > 1 and nb > 1:
+        stat_b_full = float(statistic(units_b))
+        stat_a_full = float(statistic(units_a))
+        jack: List[float] = []
+        for i in range(na):
+            jack.append(float(statistic([units_a[j] for j in range(na) if j != i])) - stat_b_full)
+        for i in range(nb):
+            jack.append(stat_a_full - float(statistic([units_b[j] for j in range(nb) if j != i])))
+        q_lo, q_hi, used = _bca_endpoints(thetas, theta_hat, jack, alpha)
+    else:
+        q_lo, q_hi, used = alpha / 2.0, 1.0 - alpha / 2.0, "percentile"
+
+    return CIResult(theta_hat, _percentile(thetas, q_lo), _percentile(thetas, q_hi),
+                    alpha, n_resamples, used, seed)
+
+
+def mean(xs: Sequence[float]) -> float:
+    """Arithmetic mean; 0.0 for an empty sequence (a bootstrap resample is never
+    empty, but degenerate jackknife folds can be)."""
+    xs = list(xs)
+    return sum(xs) / len(xs) if xs else 0.0
+
+
+# --------------------------------------------------------------------------- #
+# RQ2 — Miller–Madow bias-corrected entropy (§3 estimator [APPROVAL]).
+# --------------------------------------------------------------------------- #
+def miller_madow_entropy_bits(counts: Sequence[float]) -> float:
+    """Shannon entropy in bits with the Miller–Madow bias correction.
+
+    H_MM = H_plugin + (K − 1) / (2 N ln 2), where K is the number of observed
+    (non-zero) categories and N the total count. This corrects the systematic
+    downward bias of the plug-in (MLE) estimator at finite N. Empty/all-zero -> 0.
+    """
+    vals = [c for c in counts if c > 0]
+    total = float(sum(vals))
+    if total <= 0:
+        return 0.0
+    h = 0.0
+    for c in vals:
+        p = c / total
+        h -= p * math.log2(p)
+    k = len(vals)
+    return h + (k - 1) / (2.0 * total * math.log(2.0))
+
+
+# --------------------------------------------------------------------------- #
+# RQ2-P3 — Spearman rank correlation.
+# --------------------------------------------------------------------------- #
+def _rankdata(values: Sequence[float]) -> List[float]:
+    """Fractional ranks (ties get the average of the ranks they span)."""
+    order = sorted(range(len(values)), key=lambda i: values[i])
+    ranks = [0.0] * len(values)
+    i = 0
+    n = len(values)
+    while i < n:
+        j = i
+        while j + 1 < n and values[order[j + 1]] == values[order[i]]:
+            j += 1
+        avg = (i + j) / 2.0 + 1.0  # 1-based average rank over the tie block
+        for k in range(i, j + 1):
+            ranks[order[k]] = avg
+        i = j + 1
+    return ranks
+
+
+def _pearson(xs: Sequence[float], ys: Sequence[float]) -> float:
+    n = len(xs)
+    if n == 0 or n != len(ys):
+        return 0.0
+    mx = sum(xs) / n
+    my = sum(ys) / n
+    num = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
+    dx = math.sqrt(sum((x - mx) ** 2 for x in xs))
+    dy = math.sqrt(sum((y - my) ** 2 for y in ys))
+    if dx == 0.0 or dy == 0.0:
+        return 0.0
+    return num / (dx * dy)
+
+
+def spearman(xs: Sequence[float], ys: Sequence[float]) -> float:
+    """Spearman ρ = Pearson correlation of the fractional ranks. Returns 0.0 for
+    empty, length-mismatched, or zero-variance inputs."""
+    if len(xs) != len(ys) or not xs:
+        return 0.0
+    return _pearson(_rankdata(xs), _rankdata(ys))
+
+
+# --------------------------------------------------------------------------- #
+# Multiplicity — Holm–Bonferroni over the frozen confirmatory family.
+# --------------------------------------------------------------------------- #
+@dataclass(frozen=True)
+class HolmResult:
+    name: str
+    p: float
+    p_adjusted: float
+    reject: bool
+    rank: int  # 1-based ascending rank among the reported tests
+    multiplier: int  # Holm denominator actually used (family_size − rank + 1)
+
+
+def holm_bonferroni(
+    pvalues: dict, alpha: float = DEFAULT_ALPHA, family_size: int | None = None
+) -> List[HolmResult]:
+    """Holm step-down correction.
+
+    ``pvalues`` maps test-name -> raw p. ``family_size`` is the size of the frozen
+    confirmatory family (default = number of tests supplied). When a lead paper
+    reports a **subset** of the pre-registered family (e.g. the 4 RQ1/RQ2 tests of
+    a 7-test frozen family), pass ``family_size=7``: the k-th smallest reported p
+    is then tested against ``alpha / (family_size − k + 1)`` — i.e. the reported
+    tests are treated as occupying the *smallest* slots of the full family, giving
+    the largest (most conservative) Holm multipliers. This is strictly no less
+    stringent than the true embedded correction and can never re-optimise the
+    family down to the reported subset (which would inflate the false-rejection
+    rate and constitute p-hacking).
+
+    Adjusted p-values are made monotone non-decreasing in rank (standard Holm) and
+    capped at 1.0. ``reject`` is ``p_adjusted <= alpha``.
+    """
+    items = sorted(pvalues.items(), key=lambda kv: kv[1])
+    m = family_size if family_size is not None else len(items)
+    if m < len(items):
+        raise ValueError(f"family_size {m} < number of reported tests {len(items)}")
+
+    out: List[HolmResult] = []
+    running = 0.0
+    for k, (name, p) in enumerate(items, start=1):
+        mult = m - k + 1
+        adj = min(1.0, mult * p)
+        running = max(running, adj)  # enforce step-down monotonicity
+        out.append(HolmResult(name=name, p=p, p_adjusted=running,
+                              reject=running <= alpha, rank=k, multiplier=mult))
+    return out

tests/test_sor_stats.py

@@ -0,0 +1,164 @@
+"""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")