Spaces:
Paused
Paused
| """Anytime-valid sequential inference for paired mean differences. | |
| Why not repeated fixed-sample tests: a CLT/Wilson interval is valid at ONE | |
| pre-committed n. Recomputing it after every batch and stopping on the first | |
| significant look inflates type-I error without bound (the law of the iterated | |
| logarithm guarantees an eventual false crossing). Everything here is therefore | |
| ANYTIME-VALID: the guarantee holds simultaneously over all sample sizes and | |
| data-dependent stopping rules. | |
| Machinery (Waudby-Smith & Ramdas, "Estimating means of bounded random variables | |
| by betting", JRSS-B 2023): | |
| * `EmpiricalBernsteinCS` — predictable plug-in empirical-Bernstein confidence | |
| sequence for a mean of [0,1]-bounded observations. Variance-adaptive: for | |
| paired differences (variance << 1/4) it is far tighter than Hoeffding. | |
| * `OneSidedEProcess` — a betting supermartingale (e-process) against a | |
| one-sided null "mean >= m" or "mean <= m". Ville's inequality gives | |
| P(e-value ever >= 1/alpha) <= alpha under the null, at any stopping time. | |
| * `PairedSequentialTest` — combines both on paired differences d in [-1,1] | |
| (mapped to x=(d+1)/2), producing the mandatory THREE-OUTCOME decision | |
| against a practical-significance margin delta: | |
| REGRESSION e-process against H0: Delta >= -delta crossed 1/alpha_reg | |
| PASS e-process against H0: Delta <= -delta crossed 1/alpha_pass | |
| (an equivalence-style claim: the drop, if any, is | |
| confidently smaller than the margin) | |
| INSUFFICIENT_EVIDENCE neither — never interpreted as proven equivalence. | |
| * A Beta-Binomial posterior P(Delta < -delta | data) is computed alongside as | |
| an operational, human-readable complement (it is NOT the gate criterion). | |
| Observations fed in are CLUSTER means (see units.py), so temporal/spatial | |
| correlation is handled by construction rather than by post-hoc correction. | |
| """ | |
| from __future__ import annotations | |
| import math | |
| from typing import Dict, List, Optional | |
| from scipy import stats as _sstats | |
| DECISION_REGRESSION = "REGRESSION" | |
| DECISION_PASS = "PASS" | |
| DECISION_INSUFFICIENT = "INSUFFICIENT_EVIDENCE" | |
| class EmpiricalBernsteinCS: | |
| """Predictable plug-in empirical-Bernstein confidence sequence on [0,1]. | |
| Running-intersection bounds are kept (intersecting a CS over time preserves | |
| coverage), so the interval never widens as data accrue. | |
| """ | |
| def __init__(self, alpha: float = 0.05, truncation: float = 0.5): | |
| self.alpha = alpha | |
| self.truncation = truncation | |
| self.t = 0 | |
| self._sum_x = 0.0 | |
| self._mu_hat = 0.5 # mu_hat_0 = 1/2 (prior-regularized running mean) | |
| self._var_acc = 0.25 # 1/4 + sum_i (x_i - mu_hat_i)^2 | |
| self._sigma2_hat = 0.25 | |
| self._sum_lam = 0.0 | |
| self._sum_lam_x = 0.0 | |
| self._sum_v_psi = 0.0 | |
| self.lower = 0.0 | |
| self.upper = 1.0 | |
| def update(self, x: float) -> None: | |
| x = min(1.0, max(0.0, float(x))) | |
| t = self.t + 1 | |
| log2a = math.log(2.0 / self.alpha) | |
| lam = math.sqrt(2.0 * log2a / (max(self._sigma2_hat, 1e-12) * t * math.log1p(t))) | |
| lam = min(lam, self.truncation) | |
| v = 4.0 * (x - self._mu_hat) ** 2 | |
| psi = (-math.log1p(-lam) - lam) / 4.0 | |
| self._sum_lam += lam | |
| self._sum_lam_x += lam * x | |
| self._sum_v_psi += v * psi | |
| # advance predictable estimates AFTER betting with the old ones | |
| self.t = t | |
| self._sum_x += x | |
| self._mu_hat = (0.5 + self._sum_x) / (t + 1) | |
| self._var_acc += (x - self._mu_hat) ** 2 | |
| self._sigma2_hat = self._var_acc / (t + 1) | |
| if self._sum_lam > 0: | |
| center = self._sum_lam_x / self._sum_lam | |
| radius = (log2a + self._sum_v_psi) / self._sum_lam | |
| self.lower = max(self.lower, center - radius) | |
| self.upper = min(self.upper, center + radius) | |
| def interval(self) -> tuple: | |
| return (max(0.0, self.lower), min(1.0, self.upper)) | |
| class OneSidedEProcess: | |
| """Betting e-process for a one-sided mean null on [0,1] observations. | |
| side="below": H0: mu >= m, wealth grows when observations fall below m. | |
| side="above": H0: mu <= m, wealth grows when observations exceed m. | |
| Bet sizing: max of the empirical-Bernstein recipe (drives early growth | |
| while the sample variance estimate is still coarse) and an aGRAPA / | |
| approximate-Kelly bet lambda ~ drift / (sigma^2 + drift^2) (near-optimal | |
| asymptotic growth rate ~ KL divergence per observation once a real drift | |
| emerges). Both are predictable and clipped to the admissible range, so the | |
| wealth process remains a nonnegative supermartingale under the null — the | |
| sizing choice affects power only, never validity. | |
| """ | |
| def __init__(self, m: float, side: str, alpha_ref: float = 0.05, | |
| truncation: float = 0.75): | |
| assert side in ("below", "above") | |
| self.m = float(m) | |
| self.side = side | |
| self.alpha_ref = alpha_ref | |
| self.truncation = truncation | |
| self.t = 0 | |
| self.log_e = 0.0 | |
| self._sum_x = 0.0 | |
| self._mu_hat = 0.5 | |
| self._var_acc = 0.25 | |
| self._sigma2_hat = 0.25 | |
| def update(self, x: float) -> None: | |
| x = min(1.0, max(0.0, float(x))) | |
| t = self.t + 1 | |
| loga = math.log(1.0 / self.alpha_ref) | |
| sigma2 = max(self._sigma2_hat, 1e-12) | |
| lam_eb = math.sqrt(2.0 * loga / (sigma2 * t * math.log1p(t))) | |
| drift = (self.m - self._mu_hat) if self.side == "below" else (self._mu_hat - self.m) | |
| lam_kelly = max(drift, 0.0) / (sigma2 + drift * drift) | |
| lam = max(lam_eb, lam_kelly) | |
| if self.side == "below": | |
| lam = min(lam, self.truncation / max(1.0 - self.m, 1e-6)) | |
| factor = 1.0 - lam * (x - self.m) | |
| else: | |
| lam = min(lam, self.truncation / max(self.m, 1e-6)) | |
| factor = 1.0 + lam * (x - self.m) | |
| self.log_e += math.log(max(factor, 1e-12)) | |
| self.t = t | |
| self._sum_x += x | |
| self._mu_hat = (0.5 + self._sum_x) / (t + 1) | |
| self._var_acc += (x - self._mu_hat) ** 2 | |
| self._sigma2_hat = self._var_acc / (t + 1) | |
| def e_value(self) -> float: | |
| return math.exp(min(self.log_e, 700.0)) | |
| def p_analogue(self) -> float: | |
| """Anytime-valid p-value analogue: min(1, 1/e).""" | |
| return min(1.0, 1.0 / max(self.e_value, 1e-300)) | |
| class PairedSequentialTest: | |
| """Three-outcome anytime-valid test on paired differences d in [-1, 1]. | |
| Delta = E[d] (candidate minus baseline). Regression margin `delta` in | |
| absolute metric points: the H0 boundary sits at Delta = -delta. | |
| """ | |
| def __init__(self, delta: float = 0.005, alpha: float = 0.05, | |
| alpha_pass: Optional[float] = None): | |
| self.delta = float(delta) | |
| self.alpha = float(alpha) | |
| self.alpha_pass = float(alpha_pass if alpha_pass is not None else alpha) | |
| m = (1.0 - self.delta) / 2.0 # Delta = -delta <=> mu = (1-delta)/2 | |
| self.cs = EmpiricalBernsteinCS(alpha=self.alpha) | |
| self.e_reg = OneSidedEProcess(m, side="below", alpha_ref=self.alpha) | |
| self.e_pass = OneSidedEProcess(m, side="above", alpha_ref=self.alpha_pass) | |
| # object-level accounting (for reporting + Bayesian complement) | |
| self.n_objects = 0 | |
| self.n_clusters = 0 | |
| self.sum_baseline = 0.0 | |
| self.sum_candidate = 0.0 | |
| self.n01 = 0 # baseline correct, candidate wrong (regression pairs) | |
| self.n10 = 0 # candidate correct, baseline wrong (improvement pairs) | |
| self.decision = DECISION_INSUFFICIENT | |
| self.decided_at_n: Optional[int] = None | |
| self.trajectory: List[Dict] = [] | |
| # ---- updates ------------------------------------------------------- | |
| def update_clusters(self, cluster_means) -> None: | |
| """Feed per-cluster mean paired differences (each in [-1, 1]).""" | |
| for d in cluster_means: | |
| x = (float(d) + 1.0) / 2.0 | |
| self.cs.update(x) | |
| self.e_reg.update(x) | |
| self.e_pass.update(x) | |
| self.n_clusters += 1 | |
| def record_objects(self, b, c) -> None: | |
| """Track object-level paired outcomes for counts and the posterior.""" | |
| import numpy as np | |
| b = np.asarray(b, dtype=bool) | |
| c = np.asarray(c, dtype=bool) | |
| self.n_objects += int(b.size) | |
| self.sum_baseline += float(b.sum()) | |
| self.sum_candidate += float(c.sum()) | |
| self.n01 += int(np.sum(b & ~c)) | |
| self.n10 += int(np.sum(~b & c)) | |
| # ---- state --------------------------------------------------------- | |
| def delta_interval(self) -> tuple: | |
| lo, hi = self.cs.interval() | |
| return (2.0 * lo - 1.0, 2.0 * hi - 1.0) | |
| def delta_estimate(self) -> Optional[float]: | |
| if self.n_objects == 0: | |
| return None | |
| return (self.sum_candidate - self.sum_baseline) / self.n_objects | |
| def bayes_p_regression(self) -> Optional[float]: | |
| """Beta-Binomial posterior P(Delta < -delta | data). | |
| McNemar decomposition: among D = n01 + n10 discordant pairs let | |
| theta = P(pair is a 01 regression pair). With plug-in discordance rate | |
| nu = D/n, Delta = nu * (1 - 2*theta), so | |
| Delta < -delta <=> theta > (1 + delta*n/D) / 2. | |
| Prior theta ~ Beta(1,1). | |
| """ | |
| d_total = self.n01 + self.n10 | |
| if self.n_objects == 0 or d_total == 0: | |
| return None | |
| threshold = (1.0 + self.delta * self.n_objects / d_total) / 2.0 | |
| if threshold >= 1.0: | |
| return 0.0 | |
| return float(_sstats.beta.sf(threshold, 1 + self.n01, 1 + self.n10)) | |
| def evaluate(self, alpha_reg: Optional[float] = None, | |
| alpha_pass: Optional[float] = None) -> str: | |
| """Update the sticky three-outcome decision at the given levels.""" | |
| if self.decision != DECISION_INSUFFICIENT: | |
| return self.decision | |
| a_r = alpha_reg if alpha_reg is not None else self.alpha | |
| a_p = alpha_pass if alpha_pass is not None else self.alpha_pass | |
| if self.e_reg.e_value >= 1.0 / a_r: | |
| self.decision = DECISION_REGRESSION | |
| self.decided_at_n = self.n_objects | |
| elif self.e_pass.e_value >= 1.0 / a_p: | |
| self.decision = DECISION_PASS | |
| self.decided_at_n = self.n_objects | |
| return self.decision | |
| def snapshot(self) -> Dict: | |
| lo, hi = self.delta_interval() | |
| est = self.delta_estimate() | |
| return { | |
| "n": self.n_objects, | |
| "n_clusters": self.n_clusters, | |
| "delta_estimate": None if est is None else round(est, 5), | |
| "delta_ci": [round(lo, 5), round(hi, 5)], | |
| "e_regression": round(self.e_reg.e_value, 4), | |
| "e_pass": round(self.e_pass.e_value, 4), | |
| "p_analogue": round(self.e_reg.p_analogue(), 6), | |
| "bayes_p_regression": (None if self.bayes_p_regression() is None | |
| else round(self.bayes_p_regression(), 5)), | |
| "decision": self.decision, | |
| } | |
| def record_trajectory_point(self) -> None: | |
| """Append an evidence-vs-samples point (for the dashboard chart).""" | |
| self.trajectory.append({ | |
| "n": self.n_objects, | |
| "n_clusters": self.n_clusters, | |
| "delta_estimate": self.delta_estimate(), | |
| "delta_lower": round(self.delta_interval()[0], 5), | |
| "delta_upper": round(self.delta_interval()[1], 5), | |
| "log_e_regression": round(self.e_reg.log_e, 4), | |
| "log_e_pass": round(self.e_pass.log_e, 4), | |
| "decision": self.decision, | |
| }) | |
| def approx_mde(var_d: float, n_eff: float, delta: float, alpha: float, | |
| power: float = 0.9, anytime_inflation: float = 1.4) -> Optional[float]: | |
| """Approximate minimum detectable absolute effect at the current n_eff. | |
| Fixed-n normal approximation (z_{1-a} + z_{power}) * sqrt(var/n_eff), | |
| inflated for the anytime-valid overhead, plus the practical margin delta | |
| (effects must exceed the margin to be declarable).""" | |
| if n_eff <= 1 or var_d <= 0: | |
| return None | |
| z = _sstats.norm.ppf(1 - alpha) + _sstats.norm.ppf(power) | |
| return float(delta + anytime_inflation * z * math.sqrt(var_d / n_eff)) | |