mirror of
https://github.com/ThearchyHelios/AgentLab.git
synced 2026-11-15 06:02:01 +00:00
ThearchyHelios
185 lines
9.4 KiB
Python
185 lines
9.4 KiB
Python
"""分析用的统计函数（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
|