From 9af50fb47480dec7e625cfec9609148e3acbe352 Mon Sep 17 00:00:00 2001 From: ThearchyHelios Date: Wed, 21 Oct 2026 13:03:12 +0800 Subject: [PATCH] =?UTF-8?q?=E5=88=86=E6=9E=90=E7=AE=97=E5=AD=90=E7=9A=84?= =?UTF-8?q?=E7=BB=9F=E8=AE=A1=E5=87=BD=E6=95=B0=EF=BC=88app/analysis/stats?= =?UTF-8?q?.py=EF=BC=8C=E7=BA=AF=20Python=E3=80=81=E4=B8=8D=E5=8A=A0?= =?UTF-8?q?=E4=BE=9D=E8=B5=96=EF=BC=89=EF=BC=9A=E6=AD=A3=E6=80=81=E5=88=86?= =?UTF-8?q?=E5=B8=83=E5=87=BD=E6=95=B0=EF=BC=88math.erfc=EF=BC=89=E4=B8=8E?= =?UTF-8?q?=E5=88=86=E4=BD=8D=E6=95=B0=EF=BC=88Acklam=20=E9=80=BC=E8=BF=91?= =?UTF-8?q?=20+=20Halley=20=E4=BF=AE=E6=AD=A3=EF=BC=89=E3=80=81t=20?= =?UTF-8?q?=E5=88=86=E5=B8=83=EF=BC=88=E6=AD=A3=E5=88=99=E5=8C=96=E4=B8=8D?= =?UTF-8?q?=E5=AE=8C=E5=85=A8=20=CE=B2=20=E7=9A=84=E8=BF=9E=E5=88=86?= =?UTF-8?q?=E5=BC=8F=EF=BC=8C=E8=87=AA=E7=94=B1=E5=BA=A6=E5=BE=88=E5=A4=A7?= =?UTF-8?q?=E6=97=B6=E6=8C=89=E6=AD=A3=E6=80=81=EF=BC=89=E4=B8=8E=E5=88=86?= =?UTF-8?q?=E4=BD=8D=E6=95=B0=EF=BC=88=E4=BA=8C=E5=88=86=EF=BC=89=E3=80=81?= =?UTF-8?q?Wilson=EF=BC=88=E4=BB=8E=20eval/stats=20=E7=A7=BB=E8=BF=87?= =?UTF-8?q?=E6=9D=A5=EF=BC=8C=E8=AF=84=E6=B5=8B=E6=94=B9=E4=B8=BA=E5=BC=95?= =?UTF-8?q?=E7=94=A8=E3=80=81=E8=A1=8C=E4=B8=BA=E4=B8=8D=E5=8F=98=EF=BC=89?= =?UTF-8?q?=E3=80=81Newcombe=20=E5=8C=BA=E9=97=B4=E3=80=81=E5=90=88?= =?UTF-8?q?=E5=B9=B6=E6=AF=94=E4=BE=8B=20z=20=E6=A3=80=E9=AA=8C=E3=80=81Fi?= =?UTF-8?q?sher=20=E7=B2=BE=E7=A1=AE=E6=A3=80=E9=AA=8C=EF=BC=88lgamma?= =?UTF-8?q?=EF=BC=89=E3=80=81Welch=20t=20=E6=A3=80=E9=AA=8C=E4=B8=8E?= =?UTF-8?q?=E5=8C=BA=E9=97=B4=E3=80=81Holm=20=E6=A0=A1=E6=AD=A3=E3=80=81Ma?= =?UTF-8?q?nn=E2=80=93Kendall=EF=BC=88n=20=E2=89=A4=2010=20=E4=B8=94?= =?UTF-8?q?=E6=97=A0=E5=B9=B6=E5=88=97=E7=94=A8=E9=80=86=E5=BA=8F=E6=95=B0?= =?UTF-8?q?=E9=80=92=E6=8E=A8=E7=9A=84=E7=B2=BE=E7=A1=AE=E5=88=86=E5=B8=83?= =?UTF-8?q?=EF=BC=8C=E5=90=A6=E5=88=99=E6=AD=A3=E6=80=81=E8=BF=91=E4=BC=BC?= =?UTF-8?q?=E5=B8=A6=E5=B9=B6=E5=88=97=E4=BF=AE=E6=AD=A3=E5=92=8C=E8=BF=9E?= =?UTF-8?q?=E7=BB=AD=E6=80=A7=E6=A0=A1=E6=AD=A3=EF=BC=89=E3=80=81Theil?= =?UTF-8?q?=E2=80=93Sen=EF=BC=88=E7=B2=BE=E7=A1=AE=E7=9A=84=E6=95=B0?= =?UTF-8?q?=E5=9C=A8=E5=85=AC=E5=88=86=E6=AF=8D=E4=B8=8A=E6=8C=89=E6=95=B4?= =?UTF-8?q?=E6=95=B0=E6=AF=94=EF=BC=8C=E5=8F=AA=E5=AF=B9=E6=B5=AE=E7=82=B9?= =?UTF-8?q?=E7=9B=B8=E7=AD=89=E7=9A=84=E9=82=A3=E5=87=A0=E4=B8=AA=E7=B2=BE?= =?UTF-8?q?=E7=A1=AE=E6=AF=94=EF=BC=89=E3=80=81=E4=B8=AD=E4=BD=8D=E6=95=B0?= =?UTF-8?q?=20/=20MAD=20/=20=E5=90=8D=E6=AC=A1=EF=BC=9B=E6=B5=8B=E8=AF=95?= =?UTF-8?q?=E7=94=A8=E5=85=AC=E8=AE=A4=E8=A1=A8=E5=80=BC=E9=92=89=E4=BD=8F?= =?UTF-8?q?=EF=BC=88=E6=AD=A3=E6=80=81=200.975=20=E5=88=86=E4=BD=8D?= =?UTF-8?q?=E3=80=81t(10)=20=E4=B8=B4=E7=95=8C=E5=80=BC=E3=80=81=E5=93=81?= =?UTF-8?q?=E8=8C=B6=E8=A1=A8=E3=80=81Newcombe=20=E4=BE=8B=E5=AD=90?= =?UTF-8?q?=E3=80=81=E5=B0=8F=E6=A0=B7=E6=9C=AC=E7=B2=BE=E7=A1=AE=E5=88=86?= =?UTF-8?q?=E5=B8=83=E5=AF=B9=E7=85=A7=E6=9A=B4=E5=8A=9B=E6=9E=9A=E4=B8=BE?= =?UTF-8?q?=EF=BC=89?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- backend/app/analysis/stats.py | 467 +++++++++++++++++++++++++++ backend/app/eval/stats.py | 16 +- backend/tests/test_analysis_stats.py | 184 +++++++++++ 3 files changed, 654 insertions(+), 13 deletions(-) create mode 100644 backend/app/analysis/stats.py create mode 100644 backend/tests/test_analysis_stats.py diff --git a/backend/app/analysis/stats.py b/backend/app/analysis/stats.py new file mode 100644 index 00000000..79780b22 --- /dev/null +++ b/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 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", +] diff --git a/backend/app/eval/stats.py b/backend/app/eval/stats.py index 77f9dd22..0bfe7239 100644 --- a/backend/app/eval/stats.py +++ b/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 diff --git a/backend/tests/test_analysis_stats.py b/backend/tests/test_analysis_stats.py new file mode 100644 index 00000000..db125315 --- /dev/null +++ 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