spb/anomaly-atlas Public License

1# =============================================================================2#  Project   : anomaly-atlas3#  File      : src/anomaly_atlas/stats/multiple_testing.py4#  Purpose   : FDR/Bonferroni corrections (White RC / SPA / DSR arrive in expF)5#  Author    : Simon-Pierre Boucher6#  Contact   : contact@spboucher.ai7#  Data src  : hfmarketdata.io (sole data source)8#  Created   : 2026-08-129#  Modified  : 2026-08-1210#  Platform  : macOS / Apple Silicon (arm64)11#  License   : All rights reserved (research code)12# =============================================================================13"""Multiple-testing corrections for the scan stage (Level-0 triage).1415Scans use Benjamini–Hochberg FDR (the right error rate when a controlled16fraction of false leads into the next stage is acceptable). Level-117promotion uses SPA/StepM against artifact nulls — implemented with expF.18Validated on synthetic ground truth (§8.1 gate).19"""2021from __future__ import annotations2223import numpy as np242526def benjamini_hochberg(pvals: np.ndarray, alpha: float = 0.05) -> np.ndarray:27    """BH step-up FDR procedure. Returns a boolean rejection mask.2829    NaN p-values are never rejected but still COUNT toward m (conservative:30    an unevaluable test is a spent test, not a free one).31    """32    p = np.asarray(pvals, dtype=float)33    m = len(p)34    if m == 0:35        return np.zeros(0, dtype=bool)36    finite = np.where(np.isfinite(p))[0]37    reject = np.zeros(m, dtype=bool)38    if len(finite) == 0:39        return reject40    order = finite[np.argsort(p[finite])]41    thresholds = alpha * (np.arange(1, len(order) + 1) / m)42    passed = np.where(p[order] <= thresholds)[0]43    if len(passed):44        reject[order[: passed.max() + 1]] = True45    return reject464748def bonferroni(pvals: np.ndarray, alpha: float = 0.05) -> np.ndarray:49    """Bonferroni FWE mask (reported alongside FDR for reference)."""50    p = np.asarray(pvals, dtype=float)51    return np.isfinite(p) & (p <= alpha / max(len(p), 1))525354def bootstrap_pvalue(samples: np.ndarray, null_value: float = 0.0) -> float:55    """Two-sided percentile-bootstrap p-value of a statistic vs a null value.5657    p = 2 * min(P(boot <= null), P(boot >= null)), with the +1/(B+1)58    correction so p is never exactly 0. Triage-grade inference for scans —59    Level-1 promotion re-tests with SPA machinery.60    """61    s = np.asarray(samples, dtype=float)62    s = s[np.isfinite(s)]63    b = len(s)64    if b < 50:65        return float("nan")66    lo = (np.sum(s <= null_value) + 1) / (b + 1)67    hi = (np.sum(s >= null_value) + 1) / (b + 1)68    return float(min(1.0, 2.0 * min(lo, hi)))69