分析算子的统计函数（app/analysis/stats.py，纯 Python、不加依赖）：正态分布函数（math.erfc）与分位数（Acklam 逼近 + Halley 修正）、t 分布（正则化不完全 β 的连分式，自由度很大时按正态）与分位数（二分）、Wilson（从 eval/stats 移过来，评测改为引用、行为不变）、Newcombe 区间、合并比例 z 检验、Fisher 精确检验（lgamma）、Welch t 检验与区间、Holm 校正、Mann–Kendall（n ≤ 10 且无并列用逆序数递推的精确分布，否则正态近似带并列修正和连续性校正）、Theil–Sen（精确的数在公分母上按整数比，只对浮点相等的那几个精确比）、中位数 / MAD / 名次；测试用公认表值钉住（正态 0.975 分位、t(10) 临界值、品茶表、Newcombe 例子、小样本精确分布对照暴力枚举）

backend/app/analysis/stats.py

@@ -0,0 +1,467 @@
+"""分析用的统计函数（docs/design/analysis-round/01-算子与判定.md §3–§5，总方案 D13）：纯 Python，不加依赖、不用 numpy。
+
+为什么不装 scipy：几十 MB 的依赖只为几个函数，本机装包慢；评测本来就有一套纯 Python 的统计（Wilson 区间、精确 McNemar），
+两边要是同一套。所以 Wilson 从 eval/stats 搬到这里，评测改为引用（行为不变）。每个函数都用公认的表值钉住
+（tests/test_analysis_stats.py：正态 0.975 分位 1.959964、t(10) 双侧 0.05 临界值 2.228139、Fisher 品茶表双侧 0.4857……）。
+
+- 正态：分布函数用 math.erfc；分位数用 Acklam 的有理逼近再按分布函数做两步 Halley 修正（到浮点精度）。
+- t 分布：分布函数用正则化不完全 β 函数（连分式，Lentz 法）；自由度很大（> 1e7）时和正态差在 1e-7 以内，直接按正态算；
+  分位数按分布函数二分求逆。
+- 比率：Wilson 区间（缺省按评测的写法四舍五入到 4 位，places=None 不取整）、两比率之差的 Newcombe 区间（混合 score
+  法，即 Newcombe 1998 的方法 10）、合并比例的 z 检验、Fisher 精确检验（超几何概率用 math.lgamma 累加，双侧按「不比
+  观测到的那张表更可能」的表求和，比较时留 1e-7 的相对余量，同 R 和 scipy）。
+- 均值：Welch t 检验，自由度按 Welch–Satterthwaite。
+- 多重比较：Holm 逐步校正。
+- 趋势：Mann–Kendall（n ≤ 10 且没有并列时用精确分布：S 的分布由逆序数的分布递推算出，不查表；否则正态近似，方差按
+  并列修正，带连续性校正）、Theil–Sen 斜率（两两斜率的中位数）。
+- 中位数、中位绝对偏差、名次。
+
+能精确算的照精确算：中位数、MAD、Theil–Sen 对 Fraction 输入给 Fraction，Mann–Kendall 的精确 p 是 Fraction。浮点只出现在
+要开方、取对数、查分布的地方，结论里的数再按 analysis/numbers 统一成规范的十进制字符串。
+"""
+from __future__ import annotations
+
+import math
+from collections import Counter
+from collections.abc import Sequence
+from fractions import Fraction
+from typing import Any, NamedTuple
+
+#: 95% 区间的 z（正态 0.975 分位）
+Z95 = 1.959963984540054
+_SQRT2 = math.sqrt(2.0)
+#: t 分布自由度大到这个数以上时按正态算（相对误差约 1/自由度）
+_T_AS_NORMAL_DF = 1e7
+
+# --------------------------------------------------------------------------
+# 正态分布
+# --------------------------------------------------------------------------
+
+
+def normal_cdf(x: float) -> float:
+    """标准正态分布函数 Φ(x)。"""
+    return 0.5 * math.erfc(-x / _SQRT2)
+
+
+def normal_two_sided(z: float) -> float:
+    """标准正态的双侧尾概率 P(|Z| ≥ |z|)。"""
+    return math.erfc(abs(z) / _SQRT2)
+
+
+_A = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02, 1.383577518672690e+02,
+      -3.066479806614716e+01, 2.506628277459239e+00)
+_B = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02, 6.680131188771972e+01,
+      -1.328068155288572e+01)
+_C = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00, -2.549732539343734e+00,
+      4.374664141464968e+00, 2.938163982698783e+00)
+_D = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00)
+
+
+def normal_ppf(p: float) -> float:
+    """标准正态的 p 分位数 Φ⁻¹(p)。p ≤ 0 为负无穷，p ≥ 1 为正无穷。"""
+    if p <= 0.0:
+        return -math.inf
+    if p >= 1.0:
+        return math.inf
+    low = 0.02425
+    if p < low:
+        q = math.sqrt(-2 * math.log(p))
+        x = (((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5]) / \
+            ((((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1)
+    elif p <= 1 - low:
+        q = p - 0.5
+        r = q * q
+        x = (((((_A[0] * r + _A[1]) * r + _A[2]) * r + _A[3]) * r + _A[4]) * r + _A[5]) * q / \
+            (((((_B[0] * r + _B[1]) * r + _B[2]) * r + _B[3]) * r + _B[4]) * r + 1)
+    else:
+        q = math.sqrt(-2 * math.log1p(-p))
+        x = -(((((_C[0] * q + _C[1]) * q + _C[2]) * q + _C[3]) * q + _C[4]) * q + _C[5]) / \
+            ((((_D[0] * q + _D[1]) * q + _D[2]) * q + _D[3]) * q + 1)
+    for _ in range(2):
+        # Halley 修正：误差按分布函数现算，两步就到浮点精度
+        e = (normal_cdf(x) - p) if p < 0.5 else -(0.5 * math.erfc(x / _SQRT2) - (1 - p))
+        u = e * math.sqrt(2 * math.pi) * math.exp(x * x / 2)
+        x = x - u / (1 + x * u / 2)
+    return x
+
+
+# --------------------------------------------------------------------------
+# 正则化不完全 β 函数与 t 分布
+# --------------------------------------------------------------------------
+
+_CF_MAX_ITER = 20_000
+_CF_EPS = 1e-15
+_CF_TINY = 1e-300
+
+
+def _beta_cf(a: float, b: float, x: float) -> float:
+    """不完全 β 函数的连分式（Lentz 法，Numerical Recipes 的 betacf）。"""
+    qab, qap, qam = a + b, a + 1.0, a - 1.0
+    c = 1.0
+    d = 1.0 - qab * x / qap
+    d = 1.0 / (d if abs(d) >= _CF_TINY else _CF_TINY)
+    h = d
+    for m in range(1, _CF_MAX_ITER + 1):
+        m2 = 2 * m
+        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
+        d = 1.0 + aa * d
+        d = 1.0 / (d if abs(d) >= _CF_TINY else _CF_TINY)
+        c = 1.0 + aa / c
+        c = c if abs(c) >= _CF_TINY else _CF_TINY
+        h *= d * c
+        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
+        d = 1.0 + aa * d
+        d = 1.0 / (d if abs(d) >= _CF_TINY else _CF_TINY)
+        c = 1.0 + aa / c
+        c = c if abs(c) >= _CF_TINY else _CF_TINY
+        step = d * c
+        h *= step
+        if abs(step - 1.0) < _CF_EPS:
+            break
+    return h
+
+
+def betainc(a: float, b: float, x: float) -> float:
+    """正则化不完全 β 函数 I_x(a, b)，a、b > 0，0 ≤ x ≤ 1。"""
+    if x <= 0.0:
+        return 0.0
+    if x >= 1.0:
+        return 1.0
+    log_front = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b) + a * math.log(x) + b * math.log1p(-x)
+    if x < (a + 1.0) / (a + b + 2.0):
+        return math.exp(log_front) * _beta_cf(a, b, x) / a
+    return 1.0 - math.exp(log_front) * _beta_cf(b, a, 1.0 - x) / b
+
+
+def t_two_sided(t: float, df: float) -> float:
+    """自由度 df 的 t 分布的双侧尾概率 P(|T| ≥ |t|)。"""
+    if math.isinf(t):
+        return 0.0
+    if df > _T_AS_NORMAL_DF or math.isinf(df):
+        return normal_two_sided(t)
+    return min(1.0, betainc(df / 2.0, 0.5, df / (df + t * t)))
+
+
+def t_cdf(t: float, df: float) -> float:
+    """自由度 df 的 t 分布的分布函数。"""
+    tail = t_two_sided(t, df) / 2.0
+    return 1.0 - tail if t >= 0 else tail
+
+
+def t_ppf(p: float, df: float) -> float:
+    """自由度 df 的 t 分布的 p 分位数：按分布函数二分求逆。"""
+    if p <= 0.0:
+        return -math.inf
+    if p >= 1.0:
+        return math.inf
+    if p == 0.5:
+        return 0.0
+    if df > _T_AS_NORMAL_DF or math.isinf(df):
+        return normal_ppf(p)
+    target = min(p, 1.0 - p)
+    lo, hi = -1.0, 0.0
+    while t_cdf(lo, df) > target:
+        lo, hi = lo * 2.0, lo
+        if lo < -1e12:
+            break
+    for _ in range(200):
+        mid = (lo + hi) / 2.0
+        if t_cdf(mid, df) > target:
+            hi = mid
+        else:
+            lo = mid
+        if hi - lo <= 1e-15 * max(1.0, abs(mid)):
+            break
+    x = (lo + hi) / 2.0
+    return x if p < 0.5 else -x
+
+
+# --------------------------------------------------------------------------
+# 比率
+# --------------------------------------------------------------------------
+
+
+def wilson(passed: int, n: int, z: float = Z95, *, places: int | None = 4) -> tuple[float, float] | None:
+    """Wilson 区间（缺省 95%）。n 为 0 时没有区间。places 是四舍五入的位数（缺省 4 位，评测一直这样显示）；None 不取整。"""
+    if n <= 0:
+        return None
+    p = passed / n
+    denom = 1 + z * z / n
+    centre = (p + z * z / (2 * n)) / denom
+    half = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / denom
+    low, high = max(0.0, centre - half), min(1.0, centre + half)
+    if places is None:
+        return low, high
+    return round(low, places), round(high, places)
+
+
+def newcombe(x1: int, n1: int, x2: int, n2: int, z: float = Z95) -> tuple[float, float] | None:
+    """两比率之差 p1 − p2 的 Newcombe 区间（混合 score 法：两边各自的 Wilson 区间合成）。有一边样本为 0 时没有区间。"""
+    w1, w2 = wilson(x1, n1, z, places=None), wilson(x2, n2, z, places=None)
+    if w1 is None or w2 is None:
+        return None
+    p1, p2 = x1 / n1, x2 / n2
+    (l1, u1), (l2, u2) = w1, w2
+    diff = p1 - p2
+    delta = math.sqrt((p1 - l1) ** 2 + (u2 - p2) ** 2)
+    epsilon = math.sqrt((u1 - p1) ** 2 + (p2 - l2) ** 2)
+    return diff - delta, diff + epsilon
+
+
+def pooled_z(x1: int, n1: int, x2: int, n2: int) -> tuple[float, float] | None:
+    """两比率相等的合并比例 z 检验：(z, 双侧 p)。有一边样本为 0 时为 None；合并比例是 0 或 1（两边全是或全不是）时 z 为 0、
+    p 为 1。"""
+    if n1 <= 0 or n2 <= 0:
+        return None
+    pooled = (x1 + x2) / (n1 + n2)
+    se = math.sqrt(pooled * (1 - pooled) * (1 / n1 + 1 / n2))
+    if se == 0:
+        return 0.0, 1.0
+    z = (x1 / n1 - x2 / n2) / se
+    return z, normal_two_sided(z)
+
+
+def min_expected(x1: int, n1: int, x2: int, n2: int) -> float:
+    """2×2 表（成功、失败 × 两组）各格的期望频数里最小的那个。"""
+    total = n1 + n2
+    if total <= 0:
+        return 0.0
+    successes, failures = x1 + x2, total - x1 - x2
+    return min(r * c / total for r in (n1, n2) for c in (successes, failures))
+
+
+def _log_comb(n: int, k: int) -> float:
+    return math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)
+
+
+def fisher_exact(a: int, b: int, c: int, d: int) -> float:
+    """2×2 表 [[a, b], [c, d]] 的 Fisher 精确检验，双侧 p：边际固定时，不比观测到的表更可能出现的那些表的概率之和。"""
+    r1, r2, c1 = a + b, c + d, a + c
+    n = r1 + r2
+    if n == 0:
+        return 1.0
+    lo, hi = max(0, c1 - r2), min(r1, c1)
+    base = _log_comb(n, c1)
+
+    def log_p(x: int) -> float:
+        return _log_comb(r1, x) + _log_comb(r2, c1 - x) - base
+
+    observed = log_p(a) + 1e-7          # 相对余量 1 + 1e-7（对数里约等于加 1e-7），同 R 和 scipy
+    total = sum(math.exp(lp) for lp in (log_p(x) for x in range(lo, hi + 1)) if lp <= observed)
+    return min(1.0, total)
+
+
+# --------------------------------------------------------------------------
+# 均值
+# --------------------------------------------------------------------------
+
+
+class Welch(NamedTuple):
+    t: float
+    df: float
+    p: float
+
+
+def welch(m1: float, v1: float, n1: int, m2: float, v2: float, n2: int) -> Welch | None:
+    """两组均值相等的 Welch t 检验。v 是样本方差（n − 1 作分母）。有一组样本不到 2 个时为 None（方差算不了）。两组方差都是 0：
+    均值相同 p 为 1，不同 p 为 0。"""
+    if n1 < 2 or n2 < 2:
+        return None
+    s1, s2 = v1 / n1, v2 / n2
+    se2 = s1 + s2
+    if se2 <= 0:
+        return Welch(0.0, math.inf, 1.0) if m1 == m2 else Welch(math.copysign(math.inf, m1 - m2), math.inf, 0.0)
+    t = (m1 - m2) / math.sqrt(se2)
+    denom = (s1 * s1 / (n1 - 1) if s1 else 0.0) + (s2 * s2 / (n2 - 1) if s2 else 0.0)
+    df = se2 * se2 / denom if denom > 0 else math.inf
+    return Welch(t, df, t_two_sided(t, df))
+
+
+def welch_interval(m1: float, v1: float, n1: int, m2: float, v2: float, n2: int,
+                   level: float = 0.95) -> tuple[float, float] | None:
+    """m1 − m2 的 Welch 区间（缺省 95%）。算不了检验时也没有区间。"""
+    test = welch(m1, v1, n1, m2, v2, n2)
+    if test is None:
+        return None
+    se = math.sqrt(v1 / n1 + v2 / n2)
+    if se == 0:
+        return m1 - m2, m1 - m2
+    q = t_ppf(0.5 + level / 2.0, test.df)
+    return m1 - m2 - q * se, m1 - m2 + q * se
+
+
+# --------------------------------------------------------------------------
+# 多重比较
+# --------------------------------------------------------------------------
+
+
+def holm(pvalues: Sequence[float | None]) -> list[float | None]:
+    """Holm 逐步校正：按原始 p 从小到大，第 i 个（从 1 数）乘以 m − i + 1，再取到它为止的最大值，不超过 1。None（这一对算不了）
+    不计入次数、原样留着。"""
+    present = sorted((p, i) for i, p in enumerate(pvalues) if p is not None)
+    m = len(present)
+    out: list[float | None] = [None] * len(pvalues)
+    running = 0.0
+    for rank, (p, i) in enumerate(present):
+        running = max(running, min(1.0, (m - rank) * p))
+        out[i] = running
+    return out
+
+
+# --------------------------------------------------------------------------
+# 趋势
+# --------------------------------------------------------------------------
+
+
+def _sign(x: Any) -> int:
+    return (x > 0) - (x < 0)
+
+
+class MannKendall(NamedTuple):
+    s: int
+    var: float               # S 的方差（并列修正后）
+    z: float | None          # 正态近似的 z（精确分布时为 None）
+    p: float | Fraction      # 双侧 p；精确分布时是 Fraction
+    exact: bool
+
+
+#: 精确分布最多到几期：再多 n! 太大，正态近似已经足够
+EXACT_MAX_N = 10
+
+
+def _inversion_counts(n: int) -> list[int]:
+    """1..n 的排列里，逆序数为 k 的排列有几个（k = 0 … n(n−1)/2）。递推：第 m 个数插进 m 个位置，逆序数各加 0 … m−1。"""
+    counts = [1]
+    for m in range(2, n + 1):
+        nxt = [0] * (len(counts) + m - 1)
+        for k, c in enumerate(counts):
+            for j in range(m):
+                nxt[k + j] += c
+        counts = nxt
+    return counts
+
+
+def _scaled(values: Sequence[Any]) -> tuple[list[int], int] | None:
+    """一组精确的数（整数、Fraction）乘上公分母换成整数：(整数, 公分母)。有浮点的为 None。两两相减、比大小在整数上做，
+    比 Fraction 快得多，结果一样。"""
+    if not all(isinstance(v, (int, Fraction)) and not isinstance(v, bool) for v in values):
+        return None
+    scale = math.lcm(*(Fraction(v).denominator for v in values)) if values else 1
+    return [int(Fraction(v) * scale) for v in values], scale
+
+
+def mann_kendall(values: Sequence[Any]) -> MannKendall:
+    """Mann–Kendall 趋势检验：S = Σ_{i<j} sgn(y_j − y_i)，双侧 p。n ≤ EXACT_MAX_N 且没有并列时用精确分布（没有并列时
+    S = N − 2·逆序数，N = n(n−1)/2，逆序数的分布递推算出）；否则正态近似：Var(S) = [n(n−1)(2n+5) − Σt(t−1)(2t+5)] / 18
+    （t 是每组并列的个数），z = (S − sgn S) / √Var（连续性校正）。"""
+    n = len(values)
+    scaled = _scaled(values)
+    ys = scaled[0] if scaled is not None else list(values)
+    s = sum(_sign(ys[j] - ys[i]) for i in range(n) for j in range(i + 1, n))
+    ties = [t for t in Counter(ys).values() if t > 1]
+    var = (n * (n - 1) * (2 * n + 5) - sum(t * (t - 1) * (2 * t + 5) for t in ties)) / 18.0
+    if n < 2:
+        return MannKendall(0, 0.0, None, Fraction(1), True)
+    if not ties and n <= EXACT_MAX_N:
+        counts = _inversion_counts(n)
+        pairs = n * (n - 1) // 2
+        hit = sum(c for k, c in enumerate(counts) if abs(pairs - 2 * k) >= abs(s))
+        return MannKendall(s, var, None, Fraction(hit, math.factorial(n)), True)
+    if var <= 0:
+        return MannKendall(s, var, 0.0, 1.0, False)
+    z = (s - _sign(s)) / math.sqrt(var)
+    return MannKendall(s, var, z, normal_two_sided(z), False)
+
+
+def theil_sen(xs: Sequence[Any], ys: Sequence[Any]) -> Any:
+    """Theil–Sen 斜率：所有两点（x 不同）连线斜率的中位数。不到两个点为 None。xs、ys 都是精确的数时结果是 Fraction：斜率
+    在公分母上写成整数之比，按 int / int（正确舍入、单调）的浮点排，只把浮点相等的那几个按精确值比，n = 261 的三万多个
+    斜率也不逐个造 Fraction。"""
+    sx, sy = _scaled(xs), _scaled(ys)
+    if sx is None or sy is None:
+        slopes = [(ys[j] - ys[i]) / (xs[j] - xs[i]) for i in range(len(xs)) for j in range(i + 1, len(xs))
+                  if xs[j] != xs[i]]
+        return median(slopes) if slopes else None
+    (X, dx), (Y, dy) = sx, sy
+    pairs = []
+    for i in range(len(X)):
+        for j in range(i + 1, len(X)):
+            den = (X[j] - X[i]) * dy
+            if den:
+                num = (Y[j] - Y[i]) * dx
+                pairs.append((num, den) if den > 0 else (-num, -den))
+    if not pairs:
+        return None
+    keys = [num / den for num, den in pairs]
+    order = sorted(range(len(pairs)), key=keys.__getitem__)
+
+    def kth(k: int) -> Fraction:
+        pivot = keys[order[k]]
+        lo, hi = k, k + 1
+        while lo > 0 and keys[order[lo - 1]] == pivot:
+            lo -= 1
+        while hi < len(order) and keys[order[hi]] == pivot:
+            hi += 1
+        block = sorted(Fraction(*pairs[t]) for t in order[lo:hi])
+        return block[k - lo]
+
+    m = len(pairs)
+    return kth(m // 2) if m % 2 else (kth(m // 2 - 1) + kth(m // 2)) / 2
+
+
+# --------------------------------------------------------------------------
+# 位置与离散
+# --------------------------------------------------------------------------
+
+
+#: 超过这么多个值时中位数先按浮点排（Theil–Sen 的两两斜率有几万个，逐个按 Fraction 比太慢）
+_FAST_MEDIAN = 512
+
+
+def _kth(values: Sequence[Any], k: int) -> Any:
+    """从小到大第 k 个（从 0 数），精确：Fraction 换成浮点是正确舍入的、单调的——浮点小的精确值一定小，排错的只可能是浮点
+    相等的那几个，所以先按浮点排，再只把和第 k 个浮点相等的那一段按精确值排。"""
+    keyed = sorted((float(v), i) for i, v in enumerate(values))
+    pivot = keyed[k][0]
+    lo, hi = k, k + 1
+    while lo > 0 and keyed[lo - 1][0] == pivot:
+        lo -= 1
+    while hi < len(keyed) and keyed[hi][0] == pivot:
+        hi += 1
+    block = sorted(values[i] for _, i in keyed[lo:hi])
+    return block[k - lo]
+
+
+def median(values: Sequence[Any]) -> Any:
+    """中位数（偶数个时取中间两个的平均）。Fraction 进 Fraction 出。没有值为 None。"""
+    n = len(values)
+    if n == 0:
+        return None
+    mid = n // 2
+    if n > _FAST_MEDIAN:
+        upper = _kth(values, mid)
+        return upper if n % 2 else (_kth(values, mid - 1) + upper) / 2
+    data = sorted(values)
+    if n % 2:
+        return data[mid]
+    return (data[mid - 1] + data[mid]) / 2
+
+
+def mad(values: Sequence[Any], center: Any | None = None) -> Any:
+    """中位绝对偏差：各值离中位数（或给定的中心）的绝对距离的中位数。"""
+    if not values:
+        return None
+    c = median(values) if center is None else center
+    return median([abs(v - c) for v in values])
+
+
+def rank_desc(x: Any, others: Sequence[Any]) -> int:
+    """x 在「others + x」里从大到小的名次（并列取最好的名次：比它大的有几个再加 1）。"""
+    return 1 + sum(1 for o in others if o > x)
+
+
+__all__ = [
+    "EXACT_MAX_N", "MannKendall", "Welch", "Z95", "betainc", "fisher_exact", "holm", "mad", "mann_kendall", "median",
+    "min_expected", "newcombe", "normal_cdf", "normal_ppf", "normal_two_sided", "pooled_z", "rank_desc", "t_cdf",
+    "t_ppf", "t_two_sided", "theil_sen", "welch", "welch_interval", "wilson",
+]

backend/app/eval/stats.py

@@ -22,25 +22,15 @@ from collections import Counter
 from collections.abc import Iterable, Mapping, Sequence
 from typing import Any
 
-#: 95% 区间的 z
-Z95 = 1.959963984540054
+#: Wilson 区间和 95% 区间的 z 在分析的统计里（总方案 D13：分析和评测共用一套纯 Python 的统计），这里引用、行为不变
+from app.analysis.stats import Z95, wilson
+
 #: 问数据页等复核最多这么久（毫秒，同 answer.REVIEW_PAGE_WAIT_S）
 REVIEW_PAGE_WAIT_MS = 25_000
 TIERS = ("formal", "degraded", "withheld", "explore")
 TIMING_KEYS = ("build_ms", "run_wall_ms", "run_active_ms", "review_ms", "total_ms")
 
 
-def wilson(passed: int, n: int, z: float = Z95) -> tuple[float, float] | None:
-    """Wilson 区间（缺省 95%）。n 为 0 时没有区间。"""
-    if n <= 0:
-        return None
-    p = passed / n
-    denom = 1 + z * z / n
-    centre = (p + z * z / (2 * n)) / denom
-    half = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / denom
-    return round(max(0.0, centre - half), 4), round(min(1.0, centre + half), 4)
-
-
 def mcnemar_exact(only_a: int, only_b: int) -> float:
     """精确 McNemar 双侧 p 值：不一致的那几题里，「只 A 对」服从 Binomial(b + c, 0.5)。没有不一致的题时为 1。"""
     n = only_a + only_b

backend/tests/test_analysis_stats.py

@@ -0,0 +1,184 @@
+"""分析用的统计函数（app/analysis/stats.py）：用公认的表值钉住，外加边界——并列、n 很小、p 接近 0 / 1、t 自由度很大趋近正态。
+Wilson 从 eval/stats 搬过来，评测引用的是同一个函数、结果和以前一样。"""
+from __future__ import annotations
+
+import itertools
+import math
+import random
+from fractions import Fraction
+
+import pytest
+
+from app.analysis import stats
+from app.eval import stats as eval_stats
+
+
+def test_normal_quantiles_and_distribution_match_the_tables():
+    assert stats.normal_ppf(0.975) == pytest.approx(1.959964, abs=1e-6)
+    assert stats.normal_ppf(0.975) == pytest.approx(stats.Z95, rel=1e-15)
+    assert stats.normal_ppf(0.5) == 0.0
+    assert stats.normal_ppf(0.95) == pytest.approx(1.644854, abs=1e-6)
+    assert stats.normal_ppf(0.995) == pytest.approx(2.575829, abs=1e-6)
+    assert stats.normal_ppf(1e-10) == pytest.approx(-6.361341, abs=1e-6)       # 很靠近 0 的尾巴
+    assert stats.normal_ppf(0) == -math.inf and stats.normal_ppf(1) == math.inf
+    for p in (1e-12, 0.001, 0.02425, 0.3, 0.7, 0.97575, 0.999, 1 - 1e-9):
+        assert stats.normal_cdf(stats.normal_ppf(p)) == pytest.approx(p, rel=1e-9)
+    assert stats.normal_cdf(1.959963984540054) == pytest.approx(0.975, abs=1e-15)
+    assert stats.normal_two_sided(1.959963984540054) == pytest.approx(0.05, abs=1e-15)
+
+
+def test_t_quantiles_match_the_tables():
+    # 双侧 0.05 的临界值（t 分布表）
+    table = {1: 12.706205, 2: 4.302653, 5: 2.570582, 10: 2.228139, 30: 2.042272, 100: 1.983972}
+    for df, value in table.items():
+        assert stats.t_ppf(0.975, df) == pytest.approx(value, abs=1e-6), df
+    assert stats.t_two_sided(2.228139, 10) == pytest.approx(0.05, abs=1e-6)
+    # 自由度 2 有闭式解：F(t) = 1/2 + t / (2√(2 + t²))
+    for t in (0.3, 1.5, 4.0, 25.0):
+        assert stats.t_cdf(t, 2) == pytest.approx(0.5 + t / (2 * math.sqrt(2 + t * t)), rel=1e-12)
+    assert stats.t_cdf(0.0, 7) == 0.5 and stats.t_two_sided(0.0, 7) == 1.0
+    assert stats.t_cdf(-2.0, 10) == pytest.approx(1 - stats.t_cdf(2.0, 10), rel=1e-12)
+
+
+def test_t_with_a_huge_df_approaches_the_normal():
+    assert stats.t_ppf(0.975, 1e6) == pytest.approx(stats.Z95, abs=1e-5)
+    assert stats.t_ppf(0.975, 1e9) == pytest.approx(stats.Z95, abs=1e-9)
+    assert stats.t_two_sided(1.96, 5e6) == pytest.approx(stats.normal_two_sided(1.96), abs=1e-6)
+    # 非整数自由度（Welch）也行
+    assert stats.t_ppf(0.975, 41.784) == pytest.approx(2.018, abs=1e-3)
+
+
+def test_regularized_incomplete_beta():
+    assert stats.betainc(2, 3, 0.4) == pytest.approx(0.5248, abs=1e-12)            # 多项式，手算
+    assert stats.betainc(0.5, 0.5, 0.5) == pytest.approx(0.5, abs=1e-12)
+    assert stats.betainc(1, 1, 0.3) == pytest.approx(0.3, abs=1e-12)
+    assert stats.betainc(3, 2, 0.0) == 0.0 and stats.betainc(3, 2, 1.0) == 1.0
+
+
+def test_wilson_moved_here_is_the_same_function_the_eval_uses():
+    assert eval_stats.wilson is stats.wilson
+    assert eval_stats.Z95 == stats.Z95
+    assert stats.wilson(8, 10) == (0.4902, 0.9433)
+    assert stats.wilson(0, 0) is None
+    assert stats.wilson(0, 5)[0] == 0.0 and stats.wilson(5, 5)[1] == 1.0
+    low, high = stats.wilson(8, 10, places=None)
+    assert round(low, 4) == 0.4902 and round(high, 4) == 0.9433 and low != 0.4902
+
+
+def test_newcombe_matches_the_published_example():
+    # Newcombe（1998）的例子：56/70 对 48/80，混合 score 法（方法 10）0.0524 到 0.3339
+    low, high = stats.newcombe(56, 70, 48, 80)
+    assert low == pytest.approx(0.0524, abs=1e-4) and high == pytest.approx(0.3339, abs=1e-4)
+    assert stats.newcombe(1, 0, 2, 5) is None
+
+
+def test_pooled_z_and_min_expected():
+    z, p = stats.pooled_z(56, 70, 48, 80)
+    assert z == pytest.approx(2.6502, abs=1e-4) and p == pytest.approx(0.00805, abs=1e-5)
+    assert stats.pooled_z(0, 10, 0, 20) == (0.0, 1.0)               # 两边全不是：没有差异可言
+    assert stats.pooled_z(3, 0, 1, 5) is None
+    assert stats.min_expected(2, 40, 5, 230) == pytest.approx(7 * 40 / 270)
+
+
+def test_fisher_exact_two_sided():
+    # Fisher 的品茶实验
+    assert stats.fisher_exact(3, 1, 1, 3) == pytest.approx(0.4857142857, abs=1e-9)
+    assert stats.fisher_exact(1, 9, 11, 3) == pytest.approx(0.002759456, abs=1e-8)
+    assert stats.fisher_exact(0, 10, 10, 0) == pytest.approx(2 / math.comb(20, 10), rel=1e-9)   # p 接近 0
+    assert stats.fisher_exact(5, 5, 5, 5) == pytest.approx(1.0)                                # p 接近 1
+    assert stats.fisher_exact(0, 0, 0, 0) == 1.0
+
+
+def test_welch_t_test():
+    test = stats.welch(10, 4, 20, 8, 9, 25)
+    # t = 2 / √(4/20 + 9/25) = 2.6726；自由度按 Welch–Satterthwaite
+    assert test.t == pytest.approx(2.672612, abs=1e-6)
+    assert test.df == pytest.approx(0.56 ** 2 / (0.2 ** 2 / 19 + 0.36 ** 2 / 24), rel=1e-12)
+    assert test.p == pytest.approx(stats.t_two_sided(test.t, test.df))
+    assert stats.welch(1, 1, 1, 2, 1, 5) is None                    # 一组只有一个样本：方差算不了
+    assert stats.welch(5, 0, 3, 5, 0, 4).p == 1.0                    # 没有方差、均值一样
+    assert stats.welch(5, 0, 3, 6, 0, 4).p == 0.0
+    low, high = stats.welch_interval(10, 4, 20, 8, 9, 25)
+    assert low < 2 < high and (low + high) / 2 == pytest.approx(2)
+
+
+def test_holm_step_down():
+    assert stats.holm([0.01, 0.04, 0.03, 0.005]) == pytest.approx([0.03, 0.06, 0.06, 0.02])
+    assert stats.holm([0.01, None, 0.5]) == [0.02, None, 0.5]
+    assert stats.holm([0.6, 0.9]) == [1.0, 1.0]                     # 不超过 1
+    assert stats.holm([]) == []
+
+
+def _exact_two_sided(values: list[int]) -> Fraction:
+    """暴力枚举全部排列，算 Mann–Kendall 的精确双侧 p（对照用）。"""
+    def s_of(seq: tuple[int, ...]) -> int:
+        return sum((b > a) - (b < a) for i, a in enumerate(seq) for b in seq[i + 1:])
+
+    observed = abs(s_of(tuple(values)))
+    perms = list(itertools.permutations(values))
+    return Fraction(sum(1 for p in perms if abs(s_of(p)) >= observed), len(perms))
+
+
+def test_mann_kendall_small_samples_use_the_exact_distribution():
+    mk = stats.mann_kendall([1, 2, 3, 4])
+    assert mk.exact and mk.s == 6 and mk.p == Fraction(1, 12)              # 4! 里只有两头两个排列
+    assert stats.mann_kendall([1, 2, 3, 4, 5]).p == Fraction(1, 60)
+    for values in ([3, 1, 4, 2, 6, 5], [1, 3, 2, 5, 4, 7, 6], [2, 1, 4, 3, 6, 5, 8]):
+        assert stats.mann_kendall(values).p == _exact_two_sided(values)
+    # n = 10 的单调序列：2 / 10!
+    assert stats.mann_kendall(list(range(10))).p == Fraction(2, math.factorial(10))
+
+
+def test_mann_kendall_with_ties_or_more_points_uses_the_corrected_normal_approximation():
+    values = [1, 2, 2, 3, 4, 5, 5, 6, 7, 8, 9, 10]
+    mk = stats.mann_kendall(values)
+    n = len(values)
+    var = (n * (n - 1) * (2 * n + 5) - 2 * (2 * 1 * 9)) / 18          # 两组各两个并列
+    assert not mk.exact and mk.var == pytest.approx(var)
+    assert mk.z == pytest.approx((mk.s - 1) / math.sqrt(var))         # 连续性校正
+    assert mk.p == pytest.approx(stats.normal_two_sided(mk.z))
+    tied = stats.mann_kendall([4, 4, 4, 4, 4, 4])                      # 全是并列：方差为 0
+    assert tied.s == 0 and tied.p == 1.0
+    small_tied = stats.mann_kendall([1, 1, 2, 3])                      # n ≤ 10 但有并列：也走正态近似
+    assert not small_tied.exact
+
+
+def test_theil_sen_median_mad_and_rank():
+    xs = [Fraction(i) for i in range(4)]
+    assert stats.theil_sen(xs, [Fraction(1), Fraction(3), Fraction(2), Fraction(6)]) == Fraction(19, 12)
+    assert stats.theil_sen([Fraction(0)], [Fraction(1)]) is None
+    assert stats.median([Fraction(3), Fraction(1), Fraction(2)]) == 2
+    assert stats.median([Fraction(1), Fraction(2), Fraction(3), Fraction(10)]) == Fraction(5, 2)
+    assert stats.median([]) is None
+    assert stats.mad([Fraction(v) for v in (1, 1, 2, 2, 4, 6, 9)]) == 1
+    assert stats.rank_desc(5, [1, 7, 5, 9]) == 3 and stats.rank_desc(10, [1, 2]) == 1
+
+
+def test_theil_sen_on_many_points_matches_the_plain_definition():
+    rng = random.Random(9)
+    for n in (2, 3, 11, 60):
+        xs = [Fraction(i) for i in range(n)]
+        ys = [Fraction(rng.randint(-500, 500), rng.choice([1, 2, 3, 7, 100])) for _ in range(n)]
+        slopes = sorted((ys[j] - ys[i]) / (xs[j] - xs[i]) for i in range(n) for j in range(i + 1, n))
+        m = len(slopes)
+        expected = slopes[m // 2] if m % 2 else (slopes[m // 2 - 1] + slopes[m // 2]) / 2
+        assert stats.theil_sen(xs, ys) == expected
+    assert stats.theil_sen([0.0, 1.0, 2.0], [1.0, 3.0, 2.0]) == 0.5           # 浮点照样算
+    # 浮点分不开的斜率：精确地比
+    tiny = Fraction(1, 10 ** 25)
+    ys = [Fraction(0), Fraction(1) + tiny, Fraction(2) - tiny, Fraction(3)]
+    xs = [Fraction(i) for i in range(4)]
+    slopes = sorted((ys[j] - ys[i]) / (xs[j] - xs[i]) for i in range(4) for j in range(i + 1, 4))
+    assert stats.theil_sen(xs, ys) == (slopes[2] + slopes[3]) / 2
+
+
+def test_the_fast_median_of_many_values_is_still_exact():
+    rng = random.Random(5)
+    for size in (513, 514, 2001):
+        # 浮点分不开的值：1/3 加上远小于浮点精度的扰动
+        values = [Fraction(1, 3) + Fraction(rng.randint(-3, 3), 10 ** 30) for _ in range(size)]
+        values += [Fraction(rng.randint(-1000, 1000), rng.randint(1, 260)) for _ in range(size)]
+        ordered = sorted(values)
+        n = len(ordered)
+        expected = ordered[n // 2] if n % 2 else (ordered[n // 2 - 1] + ordered[n // 2]) / 2
+        assert stats.median(values) == expected