From cc7912f654ef25443ddf9bd281b07170c2baaaa3 Mon Sep 17 00:00:00 2001 From: Ct201314 <1195214305@qq.com> Date: Fri, 12 Jun 2026 20:03:42 +0800 Subject: [PATCH] feat(skills): add statgate academic research skill --- skills/statgate/SKILL.md | 68 ++++ skills/statgate/references/evidence-gate.md | 40 +++ .../statgate/references/statistical-tests.md | 58 ++++ skills/statgate/scripts/statgate.py | 258 +++++++++++++++ skills/statgate/scripts/stats_core.py | 308 ++++++++++++++++++ skills/statgate/tests/test_statgate.py | 114 +++++++ 6 files changed, 846 insertions(+) create mode 100644 skills/statgate/SKILL.md create mode 100644 skills/statgate/references/evidence-gate.md create mode 100644 skills/statgate/references/statistical-tests.md create mode 100644 skills/statgate/scripts/statgate.py create mode 100644 skills/statgate/scripts/stats_core.py create mode 100644 skills/statgate/tests/test_statgate.py diff --git a/skills/statgate/SKILL.md b/skills/statgate/SKILL.md new file mode 100644 index 0000000..5e11025 --- /dev/null +++ b/skills/statgate/SKILL.md @@ -0,0 +1,68 @@ +--- +name: statgate +version: 1.0.0 +description: "实验数据的证据门控与规范统计分析：先检查样本量/缺失率/偏态/方差齐性判断数据够不够下结论，再自动选择并运行规范统计检验（Welch t、配对 t、Mann-Whitney U、单因素方差分析、卡方），报告效应量 Cohen's d，生成可直接粘入论文的统计附录。当用户提到「统计分析」「跑统计检验」「数据够不够下结论」「t 检验」「方差分析」「效应量」「统计附录」时触发。" +metadata: + requires: + optional_bins: ["python"] +--- + +# statgate（证据门控与规范统计分析） + +它先问「数据够不够下这个结论」，再帮你跑规范的统计检验，不让你瞎报。统计量用 +纯标准库手算（已对照 scipy 验证准确），不依赖 numpy/scipy。 + +## 何时使用本技能 + +- 一堆实验结果，想自动跑规范的统计检验 +- 担心统计不规范、被审稿人挑数据 +- 需要一份带效应量、可直接粘进论文的统计附录 + +## 与同类工具的区别 + +同类分析工具往往「拿到数据就跑检验」。本技能最大的不同是**证据门控前置**——在跑 +任何检验前，先检查样本量是否足够、缺失率是否过高、分布是否严重偏态、方差是否悬殊， +给出「可下结论 / 谨慎 / 证据不足」判定。这把「不该下结论的数据硬跑出 p 值」这个 +科研常见错误拦在前面。且全程纯标准库，可在无 scipy 的受限环境运行。 + +## 三个核心能力 + +| 能力 | 说明 | +|------|------| +| 证据门控 | 样本量/缺失率/偏态/方差齐性检查 → 三级判定 | +| 规范检验 | 按数据形态自动选 Welch t / 配对 t / Mann-Whitney / ANOVA / 卡方 | +| 报告与附录 | 描述统计 + 检验统计量 + 效应量 + 可粘论文的附录句 | + +## 工作流 + +```bash +# 命令行直接传两组数据 +python scripts/statgate.py --groups "12,14,11,13,15" "18,20,17,19,21" + +# 配对检验 +python scripts/statgate.py --groups "10,12,11" "13,15,14" --paired true + +# 从 JSON 读多组 + 缺失原始条数 +python scripts/statgate.py --data data.json --output report.md +``` + +`data.json` 格式：`{"groups": [[...], [...]], "paired": false, "raw_counts": [30, 30]}` + +| 参数 | 说明 | +|------|------| +| `--groups` | 每组数据，逗号分隔 | +| `--data` | JSON 文件（多组 + paired + raw_counts） | +| `--paired` | `true`/`false`，是否配对（仅两组） | +| `--format` | `markdown`（默认）或 `json` | +| `--output` | 输出文件 | + +## 注意事项 + +- 检验选择：两组默认 Welch t；门控提示偏态/方差不齐时自动转 Mann-Whitney；多组用 ANOVA。 +- p 值由分布累积函数数值计算（t/F 用不完全贝塔，卡方用不完全伽马），与 scipy 一致到小数点后多位。 +- 证据门控是启发式拦截，最终判断仍需结合领域知识。 + +## References + +- [statistical-tests.md](references/statistical-tests.md) — 各检验适用条件与公式 +- [evidence-gate.md](references/evidence-gate.md) — 门控阈值与判定逻辑 diff --git a/skills/statgate/references/evidence-gate.md b/skills/statgate/references/evidence-gate.md new file mode 100644 index 0000000..26df0de --- /dev/null +++ b/skills/statgate/references/evidence-gate.md @@ -0,0 +1,40 @@ +# 证据门控阈值与判定逻辑 + +证据门控在跑检验**之前**执行，目的是拦住「不该下结论的数据」。 + +## 检查项与阈值 + +| 检查 | 阈值 | 触发后果 | +|------|------|---------| +| 样本量 | 每组 < 5 | 阻断性问题（issue） | +| 缺失率 | > 20%（需提供 raw_counts） | 提示（warning） | +| 偏度 | \|skew\| > 2.0 | 提示 + 建议非参数检验 | +| 方差齐性 | 最大/最小方差比 > 4，或存在零方差组 | 提示 + 建议非参数/Welch | + +阈值定义在 `statgate.py` 顶部常量（`MIN_N`、`SKEW_LIMIT`、`VAR_RATIO_LIMIT`），可按 +领域调整。 + +## 三级判定 + +| 判定 | 条件 | 含义 | +|------|------|------| +| 证据不足 | 存在阻断性问题（如样本量过小） | 不建议下结论，先补数据 | +| 谨慎下结论 | 无阻断问题但有提示 | 可下结论，但需处理提示（如改用非参数检验、说明缺失） | +| 可下结论 | 全部检查通过 | 可放心做参数检验 | + +## 自动联动 + +门控发现偏态或方差不齐时，`recommend_nonparametric` 置真，主分析会自动把两组检验 +从 Welch t 切换为 Mann-Whitney U，无需手动指定。 + +## 为什么要门控 + +科研中常见的错误是：样本只有两三个就跑 t 检验、数据明显偏态还硬套参数检验、 +方差悬殊却用 Student t。这些都会让 p 值失去意义。门控把这些问题在出结论前标出来， +是「让审稿人挑不出毛病」的第一道防线。 + +## 局限 + +- 偏度/方差是粗筛，不替代正式的正态性检验（如 Shapiro-Wilk）。 +- 缺失率检查需用户提供原始样本条数（raw_counts）才能算。 +- 门控是启发式，最终是否下结论仍需研究者结合领域判断。 diff --git a/skills/statgate/references/statistical-tests.md b/skills/statgate/references/statistical-tests.md new file mode 100644 index 0000000..4ea96ed --- /dev/null +++ b/skills/statgate/references/statistical-tests.md @@ -0,0 +1,58 @@ +# 统计检验适用条件与公式 + +所有统计量在 `stats_core.py` 中用 Python 标准库 math 实现，已对照 scipy 验证一致。 + +## 检验选择决策 + +| 数据形态 | 默认检验 | 门控提示偏态/方差不齐时 | +|---------|---------|----------------------| +| 两组独立 | Welch t 检验 | Mann-Whitney U | +| 两组配对 | 配对 t 检验 | （Wilcoxon，后续版本） | +| 三组及以上 | 单因素方差分析 | （Kruskal-Wallis，后续版本） | +| 列联表/频数 | 卡方独立性检验 | — | + +## Welch t 检验 + +不假设方差齐性，比 Student t 更稳健： + +``` +t = (M₁ − M₂) / √(s₁²/n₁ + s₂²/n₂) +df 用 Welch–Satterthwaite 公式近似 +``` + +p 值由 t 分布双侧累积（正则化不完全贝塔函数）算出。 + +## Mann-Whitney U + +非参数，比较两组分布位置。计算秩和（含并列均秩），用正态近似 + 连续性校正得 p 值。 +适合小样本或偏态数据。 + +## 单因素方差分析 + +``` +F = MS_between / MS_within +``` + +p 值由 F 分布上尾算出。比较三组及以上均值是否有差异。 + +## 卡方独立性检验 + +``` +χ² = Σ (O − E)² / E, E = 行和 × 列和 / 总和 +``` + +p 值由卡方分布上尾（正则化上不完全伽马）算出。 + +## 效应量 Cohen's d + +``` +d = (M₁ − M₂) / s_pooled +``` + +判读：|d|<0.2 可忽略，<0.5 小，<0.8 中等，≥0.8 大。报告效应量是现代统计规范的要求—— +p 值只说「有没有差异」，效应量说「差异有多大」。 + +## p 值精度 + +t/F 用不完全贝塔连分数展开，卡方用不完全伽马级数/连分数，与 scipy 一致到小数点后 +约 4-6 位，足以支撑 α=0.05/0.01/0.001 的显著性判断。 diff --git a/skills/statgate/scripts/statgate.py b/skills/statgate/scripts/statgate.py new file mode 100644 index 0000000..af1682e --- /dev/null +++ b/skills/statgate/scripts/statgate.py @@ -0,0 +1,258 @@ +"""statgate：实验数据的证据门控与规范统计分析。 + +先问「数据够不够下这个结论」，再帮你跑规范的统计检验，最后出一份带效应量、 +带假设检查的分析报告与统计附录——让你别瞎报，也让审稿人挑不出毛病。 + +三个核心能力： +1. 证据门控（先于检验）——在做任何检验前，检查样本量是否足够、缺失率是否过高、 + 分布是否严重偏态、组间方差是否悬殊，给出「能下结论 / 谨慎 / 证据不足」的判定 + 与原因，避免在不该下结论的数据上硬跑检验。 +2. 规范统计检验——根据数据形态自动选择并运行合适的检验（双组用 t / Mann-Whitney， + 多组用方差分析，列联表用卡方），同时报告效应量（Cohen's d），不只给 p 值。 +3. 统计附录与报告——输出符合论文写法的统计描述（M±SD、检验统计量、df、p、效应量）， + 附录可直接粘进方法/结果部分。 + +纯标准库实现（见 stats_core.py），不依赖 numpy/scipy，可在受限环境运行。 + +用法： + python statgate.py --groups "12,14,11,13" "18,20,17,19" --paired false + python statgate.py --data data.json --output report.md +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from typing import Any, Sequence + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +import stats_core as sc + +if hasattr(sys.stdout, "reconfigure"): + try: + sys.stdout.reconfigure(encoding="utf-8") + except Exception: + pass + + +MIN_N = 5 # 每组最小样本量阈值 +SKEW_LIMIT = 2.0 # 偏度绝对值告警阈值 +VAR_RATIO_LIMIT = 4.0 # 方差比告警阈值（最大/最小） + + +def evidence_gate(groups: list[list[float]], raw_counts: list[int] | None = None) -> dict[str, Any]: + """证据门控：检验前判断数据够不够下结论。""" + issues: list[str] = [] + warnings: list[str] = [] + + # 样本量 + small = [i + 1 for i, g in enumerate(groups) if len(g) < MIN_N] + if small: + issues.append(f"第 {small} 组样本量 < {MIN_N}，统计功效不足，结论不可靠。") + + # 缺失率（若提供了原始条数） + if raw_counts: + for i, (g, raw) in enumerate(zip(groups, raw_counts), start=1): + if raw > 0: + miss = (raw - len(g)) / raw + if miss > 0.2: + warnings.append(f"第 {i} 组缺失率 {miss:.0%} 偏高（>20%），需说明缺失处理。") + + # 偏态 + for i, g in enumerate(groups, start=1): + if len(g) >= 3: + sk = sc.skewness(g) + if abs(sk) > SKEW_LIMIT: + warnings.append(f"第 {i} 组偏度 {sk:.2f} 较大，分布偏离正态，建议用非参数检验。") + + # 方差齐性 + variances = [sc.variance(g) for g in groups if len(g) >= 2] + if len(variances) >= 2: + vmin, vmax = min(variances), max(variances) + if vmin == 0 and vmax > 0: + warnings.append("存在方差为 0 的组（数据无变异），与其它组方差悬殊，" + "方差严重不齐，优先用非参数检验。") + elif vmin > 0: + ratio = vmax / vmin + if ratio > VAR_RATIO_LIMIT: + warnings.append(f"组间方差比 {ratio:.1f} 偏大（>{VAR_RATIO_LIMIT}），方差不齐，" + f"优先用 Welch / 非参数检验。") + + if issues: + verdict = "证据不足" + elif warnings: + verdict = "谨慎下结论" + else: + verdict = "可下结论" + return {"verdict": verdict, "issues": issues, "warnings": warnings, + "recommend_nonparametric": any("非参数" in w for w in warnings)} + + +def describe(g: Sequence[float]) -> dict[str, float]: + return {"n": len(g), "mean": round(sc.mean(g), 4), "std": round(sc.std(g), 4), + "median": round(sc.median(g), 4)} + + +def run_analysis(groups: list[list[float]], paired: bool = False, + raw_counts: list[int] | None = None) -> dict[str, Any]: + """完整分析：门控 → 选检验 → 效应量。""" + gate = evidence_gate(groups, raw_counts) + desc = [describe(g) for g in groups] + result: dict[str, Any] = {"gate": gate, "descriptives": desc, "n_groups": len(groups)} + + use_np = gate["recommend_nonparametric"] + if len(groups) == 2: + a, b = groups + if paired: + test = sc.paired_t_test(a, b) + result["test"] = {"name": "配对 t 检验", **test} + elif use_np: + test = sc.mann_whitney_u(a, b) + result["test"] = {"name": "Mann-Whitney U 检验", **test} + else: + test = sc.welch_t_test(a, b) + result["test"] = {"name": "Welch 独立样本 t 检验", **test} + d = sc.cohens_d(a, b) + result["effect_size"] = {"cohens_d": round(d, 4), "magnitude": sc.interpret_d(d)} + elif len(groups) > 2: + test = sc.one_way_anova(*groups) + result["test"] = {"name": "单因素方差分析", **test} + else: + result["test"] = {"name": "样本组不足（需 ≥ 2 组）"} + + return result + + +def _sig_mark(p: float) -> str: + if p < 0.001: + return "***" + if p < 0.01: + return "**" + if p < 0.05: + return "*" + return "n.s." + + +def render_report(result: dict[str, Any]) -> str: + gate = result["gate"] + lines = [ + "# 统计分析报告", + "", + "## 一、证据门控", + "", + f"判定：**{gate['verdict']}**", + "", + ] + if gate["issues"]: + lines.append("阻断性问题：") + for it in gate["issues"]: + lines.append(f"- {it}") + lines.append("") + if gate["warnings"]: + lines.append("提示：") + for w in gate["warnings"]: + lines.append(f"- {w}") + lines.append("") + if not gate["issues"] and not gate["warnings"]: + lines.append("样本量、分布、方差检查均通过，可进行参数检验。") + lines.append("") + + lines += ["## 二、描述统计", "", "| 组 | n | 均值 | 标准差 | 中位数 |", "|:--:|:--:|:--:|:--:|:--:|"] + for i, d in enumerate(result["descriptives"], start=1): + lines.append(f"| {i} | {d['n']} | {d['mean']} | {d['std']} | {d['median']} |") + lines.append("") + + test = result.get("test", {}) + lines += ["## 三、假设检验", "", f"检验方法：{test.get('name', '—')}", ""] + if "p_value" in test: + p = test["p_value"] + stat_str = "" + if "t" in test: + stat_str = f"t({test['df']:.1f}) = {test['t']:.3f}" + elif "F" in test: + stat_str = f"F({test['df_between']}, {test['df_within']}) = {test['F']:.3f}" + elif "U" in test: + stat_str = f"U = {test['U']:.1f}, z = {test['z']:.3f}" + lines.append(f"- 统计量：{stat_str}") + lines.append(f"- p 值：{p:.4f} {_sig_mark(p)}") + lines.append(f"- 结论：{'差异显著' if p < 0.05 else '差异不显著'}（α = 0.05）") + lines.append("") + if "effect_size" in result: + es = result["effect_size"] + lines.append(f"- 效应量 Cohen's d = {es['cohens_d']}（{es['magnitude']}）") + lines.append("") + + # 统计附录（论文可粘） + lines += ["## 四、统计附录（可粘入论文）", ""] + lines.append(_appendix_sentence(result)) + lines += ["", "---", "", + "由 statgate 生成。先做证据门控再跑检验，统计量用纯标准库计算。" + "显著性标记：*** p<.001, ** p<.01, * p<.05, n.s. 不显著。"] + return "\n".join(lines) + + +def _appendix_sentence(result: dict[str, Any]) -> str: + """生成论文风格的一句话统计描述。""" + test = result.get("test", {}) + desc = result["descriptives"] + if "p_value" not in test: + return "（检验未执行，无法生成附录句。）" + p = test["p_value"] + p_str = "p < .001" if p < 0.001 else f"p = {p:.3f}" + if len(desc) == 2 and "t" in test: + d = result.get("effect_size", {}).get("cohens_d", 0) + return (f"两组（M₁ = {desc[0]['mean']} ± {desc[0]['std']}，" + f"M₂ = {desc[1]['mean']} ± {desc[1]['std']}）经{test['name']}比较，" + f"t({test['df']:.1f}) = {test['t']:.2f}，{p_str}，Cohen's d = {d}。") + if "F" in test: + return (f"经单因素方差分析，各组差异 " + f"F({test['df_between']}, {test['df_within']}) = {test['F']:.2f}，{p_str}。") + if "U" in test: + return f"经 Mann-Whitney U 检验，U = {test['U']:.1f}，{p_str}。" + return "（无法生成附录句。）" + + +def _parse_group(s: str) -> list[float]: + return [float(x) for x in s.replace("，", ",").split(",") if x.strip()] + + +def main(argv: list[str] | None = None) -> int: + p = argparse.ArgumentParser(prog="statgate", description="证据门控与规范统计分析") + p.add_argument("--groups", nargs="+", help='每组数据，如 "12,14,11" "18,20,17"') + p.add_argument("--data", type=Path, help="JSON 文件：{\"groups\": [[...],[...]], \"paired\": false}") + p.add_argument("--paired", choices=["true", "false"], default="false", help="是否配对（仅两组）") + p.add_argument("--format", choices=["markdown", "json"], default="markdown") + p.add_argument("--output", type=Path) + args = p.parse_args(argv) + + paired = args.paired == "true" + groups: list[list[float]] = [] + raw_counts = None + if args.data and args.data.exists(): + cfg = json.loads(args.data.read_text(encoding="utf-8-sig")) + groups = [[float(v) for v in g] for g in cfg.get("groups", [])] + paired = bool(cfg.get("paired", paired)) + raw_counts = cfg.get("raw_counts") + elif args.groups: + groups = [_parse_group(g) for g in args.groups] + + if len(groups) < 2: + print("错误：请用 --groups 或 --data 提供至少两组数据。", file=sys.stderr) + return 2 + + result = run_analysis(groups, paired=paired, raw_counts=raw_counts) + out = (json.dumps(result, ensure_ascii=False, indent=2) if args.format == "json" + else render_report(result)) + if args.output: + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(out, encoding="utf-8") + print(f"已写入 {args.output}") + else: + print(out) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/skills/statgate/scripts/stats_core.py b/skills/statgate/scripts/stats_core.py new file mode 100644 index 0000000..a9afd5b --- /dev/null +++ b/skills/statgate/scripts/stats_core.py @@ -0,0 +1,308 @@ +"""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 "大" diff --git a/skills/statgate/tests/test_statgate.py b/skills/statgate/tests/test_statgate.py new file mode 100644 index 0000000..e89a5a0 --- /dev/null +++ b/skills/statgate/tests/test_statgate.py @@ -0,0 +1,114 @@ +"""statgate 单元测试。含统计量数值正确性断言（对照已知值）。""" +from __future__ import annotations +import sys +import math +from pathlib import Path +sys.path.insert(0, str(Path(__file__).resolve().parent.parent / "scripts")) +import pytest +import stats_core as sc +from statgate import evidence_gate, run_analysis, render_report + + +A = [12, 14, 11, 13, 15, 12, 14] +B = [18, 20, 17, 19, 21, 18, 20] + + +class TestDescriptive: + def test_mean_std(self): + assert sc.mean([2, 4, 6]) == 4.0 + assert sc.std([2, 4, 6], ddof=1) == pytest.approx(2.0) + + def test_median(self): + assert sc.median([3, 1, 2]) == 2 + assert sc.median([1, 2, 3, 4]) == 2.5 + + +class TestWelch: + def test_matches_scipy(self): + # scipy: t=-7.9373, p=0.000004 + r = sc.welch_t_test(A, B) + assert r["t"] == pytest.approx(-7.9373, abs=1e-3) + assert r["p_value"] == pytest.approx(0.000004, abs=1e-5) + + +class TestMannWhitney: + def test_u_value(self): + r = sc.mann_whitney_u(A, B) + assert r["U"] == pytest.approx(0.0, abs=0.5) + assert r["p_value"] < 0.05 + + +class TestAnova: + def test_matches_scipy(self): + # scipy F=148.1667 + r = sc.one_way_anova(A, B, [25, 27, 24, 26, 28, 25, 27]) + assert r["F"] == pytest.approx(148.1667, abs=1e-2) + assert r["p_value"] < 0.001 + + +class TestChiSquare: + def test_matches_scipy(self): + # scipy chi2=0.7937, p=0.373 + r = sc.chi_square_test([[10, 20], [30, 40]]) + assert r["chi2"] == pytest.approx(0.7937, abs=1e-3) + assert r["p_value"] == pytest.approx(0.373, abs=1e-2) + + +class TestCohensD: + def test_large_effect(self): + d = sc.cohens_d(A, B) + assert d == pytest.approx(-4.2426, abs=1e-3) + assert sc.interpret_d(d) == "大" + + def test_magnitudes(self): + assert sc.interpret_d(0.1) == "可忽略" + assert sc.interpret_d(0.3) == "小" + assert sc.interpret_d(0.6) == "中等" + + +class TestPairedEdge: + def test_constant_nonzero_diff(self): + # 每对差值恒为 -6，无变异且非零 → 极显著 + r = sc.paired_t_test([1, 2, 3], [7, 8, 9]) + assert r["p_value"] == 0.0 + assert math.isinf(r["t"]) + + def test_identical(self): + r = sc.paired_t_test([1, 2, 3], [1, 2, 3]) + assert r["p_value"] == 1.0 + + +class TestEvidenceGate: + def test_small_sample_blocks(self): + g = evidence_gate([[1, 2], [3, 4]]) + assert g["verdict"] == "证据不足" + assert g["issues"] + + def test_normal_passes(self): + g = evidence_gate([A, B]) + assert g["verdict"] == "可下结论" + + def test_high_variance_warns(self): + g = evidence_gate([[1, 1, 1, 1, 1, 1], [1, 50, 1, 99, 1, 80]]) + assert g["warnings"] + + +class TestRunAnalysis: + def test_two_group_ttest(self): + r = run_analysis([A, B]) + assert "t" in r["test"] + assert "effect_size" in r + + def test_multi_group_anova(self): + r = run_analysis([A, B, [25, 27, 24, 26, 28, 25, 27]]) + assert r["test"]["name"] == "单因素方差分析" + + def test_render(self): + md = render_report(run_analysis([A, B])) + assert "统计分析报告" in md + assert "证据门控" in md + assert "统计附录" in md + + +if __name__ == "__main__": + sys.exit(pytest.main([__file__, "-v"]))