forked from Gitlink/gitlink-cli
309 lines
9.6 KiB
Python
309 lines
9.6 KiB
Python
"""statgate 统计核心：纯标准库实现的统计量计算。
|
||
|
||
只用 Python 标准库 math，不依赖 numpy/scipy。实现研究中最常用的统计检验与
|
||
效应量，并给出 p 值的解析近似或精确分布。覆盖：
|
||
|
||
- 描述统计：均值、标准差、中位数、四分位
|
||
- 正态性粗检：偏度/峰度（作为是否走参数检验的参考）
|
||
- 双样本：独立 t 检验（Welch）、配对 t 检验、Mann-Whitney U（秩和）
|
||
- 多组：单因素方差分析（One-way ANOVA）
|
||
- 分类：卡方独立性检验
|
||
- 效应量：Cohen's d、秩双列相关 r
|
||
|
||
p 值通过分布的累积函数近似计算（t 分布用数值积分，正态用 erf，卡方/F 用
|
||
不完全伽马/贝塔的级数展开）。精度足以支撑研究中的显著性判断。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from typing import Sequence
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 描述统计
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def mean(xs: Sequence[float]) -> float:
|
||
return sum(xs) / len(xs) if xs else 0.0
|
||
|
||
|
||
def variance(xs: Sequence[float], ddof: int = 1) -> float:
|
||
n = len(xs)
|
||
if n - ddof <= 0:
|
||
return 0.0
|
||
m = mean(xs)
|
||
return sum((x - m) ** 2 for x in xs) / (n - ddof)
|
||
|
||
|
||
def std(xs: Sequence[float], ddof: int = 1) -> float:
|
||
return math.sqrt(variance(xs, ddof))
|
||
|
||
|
||
def median(xs: Sequence[float]) -> float:
|
||
s = sorted(xs)
|
||
n = len(s)
|
||
if n == 0:
|
||
return 0.0
|
||
mid = n // 2
|
||
return s[mid] if n % 2 else (s[mid - 1] + s[mid]) / 2
|
||
|
||
|
||
def skewness(xs: Sequence[float]) -> float:
|
||
n = len(xs)
|
||
if n < 3:
|
||
return 0.0
|
||
m, sd = mean(xs), std(xs, ddof=1)
|
||
if sd == 0:
|
||
return 0.0
|
||
return (n / ((n - 1) * (n - 2))) * sum(((x - m) / sd) ** 3 for x in xs)
|
||
|
||
|
||
def kurtosis(xs: Sequence[float]) -> float:
|
||
"""超额峰度（正态为 0）。"""
|
||
n = len(xs)
|
||
if n < 4:
|
||
return 0.0
|
||
m, sd = mean(xs), std(xs, ddof=1)
|
||
if sd == 0:
|
||
return 0.0
|
||
g2 = sum(((x - m) / sd) ** 4 for x in xs)
|
||
return (n * (n + 1) / ((n - 1) * (n - 2) * (n - 3))) * g2 - 3 * (n - 1) ** 2 / ((n - 2) * (n - 3))
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 分布累积函数（用于 p 值）
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def _norm_cdf(z: float) -> float:
|
||
return 0.5 * (1 + math.erf(z / math.sqrt(2)))
|
||
|
||
|
||
def _betacf(a: float, b: float, x: float) -> float:
|
||
"""连分数展开（Numerical Recipes 思路，独立实现）。"""
|
||
MAXIT, EPS, FPMIN = 200, 3e-12, 1e-30
|
||
qab, qap, qam = a + b, a + 1, a - 1
|
||
c = 1.0
|
||
d = 1 - qab * x / qap
|
||
if abs(d) < FPMIN:
|
||
d = FPMIN
|
||
d = 1 / d
|
||
h = d
|
||
for m in range(1, MAXIT + 1):
|
||
m2 = 2 * m
|
||
aa = m * (b - m) * x / ((qam + m2) * (a + m2))
|
||
d = 1 + aa * d
|
||
if abs(d) < FPMIN:
|
||
d = FPMIN
|
||
c = 1 + aa / c
|
||
if abs(c) < FPMIN:
|
||
c = FPMIN
|
||
d = 1 / d
|
||
h *= d * c
|
||
aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
|
||
d = 1 + aa * d
|
||
if abs(d) < FPMIN:
|
||
d = FPMIN
|
||
c = 1 + aa / c
|
||
if abs(c) < FPMIN:
|
||
c = FPMIN
|
||
d = 1 / d
|
||
de = d * c
|
||
h *= de
|
||
if abs(de - 1) < EPS:
|
||
break
|
||
return h
|
||
|
||
|
||
def _betai(a: float, b: float, x: float) -> float:
|
||
"""正则化不完全贝塔函数 I_x(a,b)。"""
|
||
if x <= 0:
|
||
return 0.0
|
||
if x >= 1:
|
||
return 1.0
|
||
lbeta = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
|
||
bt = math.exp(lbeta + a * math.log(x) + b * math.log(1 - x))
|
||
if x < (a + 1) / (a + b + 2):
|
||
return bt * _betacf(a, b, x) / a
|
||
return 1 - bt * _betacf(b, a, 1 - x) / b
|
||
|
||
|
||
def _t_sf_two_sided(t: float, df: float) -> float:
|
||
"""t 分布双侧 p 值。"""
|
||
if df <= 0:
|
||
return float("nan")
|
||
x = df / (df + t * t)
|
||
return _betai(df / 2, 0.5, x)
|
||
|
||
|
||
def _f_sf(f: float, df1: float, df2: float) -> float:
|
||
"""F 分布上尾 p 值。"""
|
||
if f <= 0:
|
||
return 1.0
|
||
x = df2 / (df2 + df1 * f)
|
||
return _betai(df2 / 2, df1 / 2, x)
|
||
|
||
|
||
def _gammainc_upper_reg(s: float, x: float) -> float:
|
||
"""正则化上不完全伽马 Q(s,x)，用于卡方上尾 p。"""
|
||
if x <= 0:
|
||
return 1.0
|
||
if x < s + 1:
|
||
# 用下不完全的级数
|
||
term = 1.0 / s
|
||
summ = term
|
||
n = s
|
||
for _ in range(200):
|
||
n += 1
|
||
term *= x / n
|
||
summ += term
|
||
if abs(term) < abs(summ) * 1e-12:
|
||
break
|
||
p = summ * math.exp(-x + s * math.log(x) - math.lgamma(s))
|
||
return 1 - p
|
||
# 连分数
|
||
FPMIN = 1e-30
|
||
b = x + 1 - s
|
||
c = 1 / FPMIN
|
||
d = 1 / b
|
||
h = d
|
||
for i in range(1, 200):
|
||
an = -i * (i - s)
|
||
b += 2
|
||
d = an * d + b
|
||
if abs(d) < FPMIN:
|
||
d = FPMIN
|
||
c = b + an / c
|
||
if abs(c) < FPMIN:
|
||
c = FPMIN
|
||
d = 1 / d
|
||
de = d * c
|
||
h *= de
|
||
if abs(de - 1) < 1e-12:
|
||
break
|
||
return h * math.exp(-x + s * math.log(x) - math.lgamma(s))
|
||
|
||
|
||
def _chi2_sf(chi2: float, df: int) -> float:
|
||
return _gammainc_upper_reg(df / 2, chi2 / 2)
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 检验
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def welch_t_test(a: Sequence[float], b: Sequence[float]) -> dict[str, float]:
|
||
"""Welch 独立样本 t 检验（不假设方差齐）。"""
|
||
na, nb = len(a), len(b)
|
||
ma, mb = mean(a), mean(b)
|
||
va, vb = variance(a), variance(b)
|
||
se = math.sqrt(va / na + vb / nb)
|
||
if se == 0:
|
||
return {"t": 0.0, "df": na + nb - 2, "p_value": 1.0}
|
||
t = (ma - mb) / se
|
||
df = (va / na + vb / nb) ** 2 / ((va / na) ** 2 / (na - 1) + (vb / nb) ** 2 / (nb - 1))
|
||
return {"t": t, "df": df, "p_value": _t_sf_two_sided(t, df)}
|
||
|
||
|
||
def paired_t_test(a: Sequence[float], b: Sequence[float]) -> dict[str, float]:
|
||
"""配对 t 检验。"""
|
||
if len(a) != len(b):
|
||
raise ValueError("配对 t 检验要求两组样本长度相同")
|
||
diffs = [x - y for x, y in zip(a, b)]
|
||
n = len(diffs)
|
||
md, sd = mean(diffs), std(diffs)
|
||
se = sd / math.sqrt(n) if n else 0
|
||
if se == 0:
|
||
# 差值无变异：若均值也为 0，两组完全相同（p=1）；
|
||
# 若均值非 0（每对差值相同且非零），则为完全分离，视为极显著。
|
||
if md == 0:
|
||
return {"t": 0.0, "df": n - 1, "p_value": 1.0}
|
||
return {"t": float("inf") if md > 0 else float("-inf"), "df": n - 1, "p_value": 0.0}
|
||
t = md / se
|
||
return {"t": t, "df": n - 1, "p_value": _t_sf_two_sided(t, n - 1)}
|
||
|
||
|
||
def mann_whitney_u(a: Sequence[float], b: Sequence[float]) -> dict[str, float]:
|
||
"""Mann-Whitney U 检验（正态近似，含连续性校正）。"""
|
||
na, nb = len(a), len(b)
|
||
combined = [(v, 0) for v in a] + [(v, 1) for v in b]
|
||
combined.sort(key=lambda x: x[0])
|
||
# 秩（含并列均秩）
|
||
ranks = [0.0] * len(combined)
|
||
i = 0
|
||
while i < len(combined):
|
||
j = i
|
||
while j + 1 < len(combined) and combined[j + 1][0] == combined[i][0]:
|
||
j += 1
|
||
avg_rank = (i + j) / 2 + 1
|
||
for k in range(i, j + 1):
|
||
ranks[k] = avg_rank
|
||
i = j + 1
|
||
r1 = sum(ranks[k] for k in range(len(combined)) if combined[k][1] == 0)
|
||
u1 = r1 - na * (na + 1) / 2
|
||
u = min(u1, na * nb - u1)
|
||
mu = na * nb / 2
|
||
sigma = math.sqrt(na * nb * (na + nb + 1) / 12)
|
||
if sigma == 0:
|
||
return {"U": u, "z": 0.0, "p_value": 1.0}
|
||
z = (u - mu + 0.5) / sigma
|
||
p = 2 * _norm_cdf(z)
|
||
return {"U": u, "z": z, "p_value": min(1.0, p)}
|
||
|
||
|
||
def one_way_anova(*groups: Sequence[float]) -> dict[str, float]:
|
||
"""单因素方差分析。"""
|
||
k = len(groups)
|
||
all_vals = [v for g in groups for v in g]
|
||
n = len(all_vals)
|
||
grand = mean(all_vals)
|
||
ss_between = sum(len(g) * (mean(g) - grand) ** 2 for g in groups)
|
||
ss_within = sum((v - mean(g)) ** 2 for g in groups for v in g)
|
||
df_b, df_w = k - 1, n - k
|
||
if df_w <= 0 or ss_within == 0:
|
||
return {"F": 0.0, "df_between": df_b, "df_within": df_w, "p_value": 1.0}
|
||
ms_b, ms_w = ss_between / df_b, ss_within / df_w
|
||
f = ms_b / ms_w
|
||
return {"F": f, "df_between": df_b, "df_within": df_w, "p_value": _f_sf(f, df_b, df_w)}
|
||
|
||
|
||
def chi_square_test(table: Sequence[Sequence[float]]) -> dict[str, float]:
|
||
"""卡方独立性检验（列联表）。"""
|
||
rows = len(table)
|
||
cols = len(table[0])
|
||
total = sum(sum(r) for r in table)
|
||
row_sums = [sum(r) for r in table]
|
||
col_sums = [sum(table[i][j] for i in range(rows)) for j in range(cols)]
|
||
chi2 = 0.0
|
||
for i in range(rows):
|
||
for j in range(cols):
|
||
exp = row_sums[i] * col_sums[j] / total if total else 0
|
||
if exp > 0:
|
||
chi2 += (table[i][j] - exp) ** 2 / exp
|
||
df = (rows - 1) * (cols - 1)
|
||
return {"chi2": chi2, "df": df, "p_value": _chi2_sf(chi2, df)}
|
||
|
||
|
||
# ---------------------------------------------------------------------------
|
||
# 效应量
|
||
# ---------------------------------------------------------------------------
|
||
|
||
def cohens_d(a: Sequence[float], b: Sequence[float]) -> float:
|
||
"""Cohen's d（合并标准差）。"""
|
||
na, nb = len(a), len(b)
|
||
sp2 = ((na - 1) * variance(a) + (nb - 1) * variance(b)) / (na + nb - 2)
|
||
if sp2 <= 0:
|
||
return 0.0
|
||
return (mean(a) - mean(b)) / math.sqrt(sp2)
|
||
|
||
|
||
def interpret_d(d: float) -> str:
|
||
ad = abs(d)
|
||
if ad < 0.2:
|
||
return "可忽略"
|
||
if ad < 0.5:
|
||
return "小"
|
||
if ad < 0.8:
|
||
return "中等"
|
||
return "大"
|