""" ============================================================================= CMU 17-630 · 前置讲 P3 共享模块：stats_kit.py —— 提示实验需要的那点统计学 ============================================================================= 【这个模块是做什么的】 "提示 B 的准确率是 66%，A 是 62%，B 更好"——这句话可能是真的，也可能只是运气。 这个模块提供判断它需要的工具，全部只用 Python 标准库，每个函数都写清楚了公式和适用范围： 区间 wald_ci, wilson_ci 一个准确率的 95% 置信区间 bootstrap_ci 任意统计量的自助法（bootstrap）区间 比较 two_proportion_test 两组**不同**的题上比较两个准确率（不配对） mcnemar_exact **同一组**题上比较两个提示（配对） paired_diff_ci 配对差值的置信区间（正态近似） paired_bootstrap_ci 配对差值的 bootstrap 区间 多重 bonferroni, holm 试了很多个提示时，校正 p 值 规划 n_for_ci_halfwidth 区间想窄到 ±h，要多少道题 n_paired_mcnemar, n_unpaired 想分辨出 δ 的差别（80% 功效），要多少道题 聚类 cluster_se 每道题采样了好几次时，正确的标准误 以及一个"模拟的提示实验世界"（ItemBank）：题目有难度，提示有"能力"，同一道题上不同提示的 结果是**相关的**（都受题目难度影响）——这正是真实评测的样子，也是配对检验更灵敏的原因。 【运行】python3 stats_kit.py # 自检：几个公式和模拟结果的对照 【依赖】只用 Python 标准库（math、random、statistics）。 """ import math import random from dataclasses import dataclass from statistics import NormalDist Z = NormalDist() # 标准正态分布：Z.cdf(x) 是 P(X ≤ x)，Z.inv_cdf(q) 是它的反函数 # --------------------------------------------------------------------------- # 1. 一个准确率的置信区间 # --------------------------------------------------------------------------- def z_for(conf=0.95): """双侧置信水平 conf 对应的正态分位数，例如 0.95 → 1.96。""" return Z.inv_cdf(0.5 + conf / 2) def wald_ci(k, n, conf=0.95): """ 最朴素的区间：p̂ ± z · sqrt(p̂(1−p̂)/n)。 【边界】p̂ = 0 或 1 时宽度为 0（"绝对不可能出错"）；n 小、p 靠近 0 或 1 时覆盖率明显不足 95%（01 号脚本演示 2）。 结果截断到 [0, 1]。n = 0 时返回 (0, 1)。 """ if n == 0: return 0.0, 1.0 p = k / n h = z_for(conf) * math.sqrt(p * (1 - p) / n) return max(0.0, p - h), min(1.0, p + h) def wilson_ci(k, n, conf=0.95): """ Wilson 区间：把"p 落在 p̂ ± z·sqrt(p(1−p)/n) 里"这个不等式对 p 解出来，而不是先把 p 换成 p̂。 【公式】center = (p̂ + z²/2n)/(1 + z²/n)，half = z·sqrt(p̂(1−p̂)/n + z²/4n²)/(1 + z²/n) 【为什么更好】它的覆盖率在小样本、极端 p 时也接近名义水平，而且永远在 [0, 1] 里。推荐默认用它。 """ if n == 0: return 0.0, 1.0 z = z_for(conf) p = k / n d = 1 + z * z / n c = (p + z * z / (2 * n)) / d h = z * math.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / d return max(0.0, c - h), min(1.0, c + h) def bootstrap_ci(values, stat=None, B=2000, conf=0.95, seed=0): """ 自助法区间：从数据里**有放回地**重抽同样多个，算一次统计量；重复 B 次，取中间 conf 那一段（百分位法）。 【参数】values：数据列表（例如每道题 0/1）；stat：统计量函数（默认均值）；B：重抽次数 【思路】我们不知道"再出一份题目"会得到什么结果，就用手上这份题目代替总体，模拟"再出一份"。 它不依赖正态近似，可以用于任何统计量（中位数、F1、两个指标的差……）。 【边界】数据很少（< 20 个）时，bootstrap 区间往往偏窄；B 太小时区间端点不稳定。 """ rng = random.Random(seed) stat = stat or (lambda xs: math.fsum(xs) / len(xs)) n = len(values) reps = sorted(stat([values[rng.randrange(n)] for _ in range(n)]) for _ in range(B)) lo = reps[int((1 - conf) / 2 * B)] hi = reps[min(B - 1, int((1 + conf) / 2 * B))] return lo, hi # --------------------------------------------------------------------------- # 2. 比较两个提示 # --------------------------------------------------------------------------- def two_proportion_test(k1, n1, k2, n2): """ 不配对的两比例 z 检验：两组**不同**的题（或者当作不同的题）上，准确率 k1/n1 和 k2/n2 有没有差别。 【返回】(差值 p2 − p1, z 统计量, 双侧 p 值) 【公式】合并比例 p̄ = (k1+k2)/(n1+n2)，SE = sqrt(p̄(1−p̄)(1/n1 + 1/n2))，z = (p2 − p1)/SE 【边界】它假设两组结果相互独立。两个提示在**同一组**题上测出来的结果并不独立（都受题目难度影响）， 这时用它会**浪费**统计功效（02 号脚本）——应该用 mcnemar_exact 或 paired_diff_ci。 """ p1, p2 = k1 / n1, k2 / n2 pbar = (k1 + k2) / (n1 + n2) se = math.sqrt(pbar * (1 - pbar) * (1 / n1 + 1 / n2)) if se == 0: return p2 - p1, 0.0, 1.0 z = (p2 - p1) / se return p2 - p1, z, 2 * (1 - Z.cdf(abs(z))) def mcnemar_exact(a, b): """ McNemar 精确检验：同一组题上，提示 A 和 B 的 0/1 结果有没有系统性的差别。 【参数】a, b：等长的 0/1 列表（a[i] 是 A 在第 i 道题上对不对） 【返回】(n01：A 错 B 对的题数, n10：A 对 B 错的题数, 双侧 p 值) 【思路】两个提示都对、都错的题，对"谁更好"没有任何信息；只看**不一致**的题。 如果 A、B 一样好，每道不一致的题是"A 赢"还是"B 赢"就像抛硬币：n10 ~ Binomial(n01 + n10, 1/2)。 p 值 = 两侧的二项尾概率。 【边界】不一致的题为 0 时 p = 1。用的是精确的二项分布，不需要大样本近似。 """ n01 = sum(1 for x, y in zip(a, b) if x == 0 and y == 1) n10 = sum(1 for x, y in zip(a, b) if x == 1 and y == 0) n = n01 + n10 if n == 0: return n01, n10, 1.0 m = min(n01, n10) tail = math.fsum(math.comb(n, i) for i in range(m + 1)) / 2 ** n return n01, n10, min(1.0, 2 * tail) def paired_diff_ci(a, b, conf=0.95): """ 配对差值的正态近似区间：d_i = b_i − a_i，区间 = mean(d) ± z · sd(d)/sqrt(n)。 【返回】(差值的估计, 下界, 上界) 【为什么比不配对的区间窄】Var(b − a) = Var(a) + Var(b) − 2·Cov(a, b)。同一道题上两个提示的结果正相关 （难题大家都错、简单题大家都对），协方差项把方差减掉了一大块。 """ n = len(a) d = [y - x for x, y in zip(a, b)] m = math.fsum(d) / n var = math.fsum((x - m) ** 2 for x in d) / (n - 1) if n > 1 else 0.0 h = z_for(conf) * math.sqrt(var / n) return m, m - h, m + h def paired_bootstrap_ci(a, b, B=2000, conf=0.95, seed=0): """配对 bootstrap：每次重抽**题目**（A、B 的结果跟着题目一起走），算差值的均值。返回 (下界, 上界)。""" pairs = list(zip(a, b)) return bootstrap_ci(pairs, stat=lambda ps: math.fsum(y - x for x, y in ps) / len(ps), B=B, conf=conf, seed=seed) # --------------------------------------------------------------------------- # 3. 多重比较 # --------------------------------------------------------------------------- def bonferroni(pvals, alpha=0.05): """Bonferroni：每个检验用 α/m 的门槛（m 是检验的个数）。返回每个检验是否显著的列表。""" m = len(pvals) return [p <= alpha / m for p in pvals] def holm(pvals, alpha=0.05): """ Holm 逐步校正：把 p 值从小到大排序，第 i 小的（i 从 0 开始）和 α/(m − i) 比较， 一旦有一个不显著，后面的全部不显著。 【为什么】它和 Bonferroni 一样保证"至少有一个假阳性"的概率 ≤ α（族错误率，FWER）， 但比 Bonferroni 更容易发现真正的差别。 【返回】和输入顺序一致的布尔列表 """ m = len(pvals) order = sorted(range(m), key=lambda i: pvals[i]) out = [False] * m for rank, i in enumerate(order): if pvals[i] <= alpha / (m - rank): out[i] = True else: break return out # --------------------------------------------------------------------------- # 4. 样本量 # --------------------------------------------------------------------------- def n_for_ci_halfwidth(p, h, conf=0.95): """想让准确率的置信区间半宽 ≤ h（例如 0.02 = ±2 个百分点），大约要多少道题：n ≈ z²·p(1−p)/h²。""" z = z_for(conf) return math.ceil(z * z * p * (1 - p) / (h * h)) def n_unpaired(p1, p2, alpha=0.05, power=0.8): """ 不配对的比较：每组要多少道题，才能以 power 的概率在 α 水平上发现 p1 和 p2 的差别。 【公式】n = (z_{1−α/2}·sqrt(2·p̄q̄) + z_{power}·sqrt(p1q1 + p2q2))² / (p1 − p2)²，p̄ = (p1+p2)/2 """ za, zb = z_for(1 - alpha), Z.inv_cdf(power) pbar = (p1 + p2) / 2 num = za * math.sqrt(2 * pbar * (1 - pbar)) + zb * math.sqrt(p1 * (1 - p1) + p2 * (1 - p2)) return math.ceil(num * num / (p1 - p2) ** 2) def n_paired_mcnemar(delta, psi, alpha=0.05, power=0.8): """ 配对的比较（McNemar）：要多少道题，才能以 power 的概率发现差值 delta。 【参数】delta：两个提示准确率的差（= p10 − p01，B 的准确率减 A 的）； psi：不一致的题占的比例（= p01 + p10），通常从一个小规模的预实验里估计 【公式】（Connor 1987 的近似）n = (z_{1−α/2}·sqrt(ψ) + z_{power}·sqrt(ψ − δ²))² / δ² 【直觉】需要的题数大致和 ψ/δ² 成正比：两个提示在越多的题上"意见不同"，同样的平均差就越难分辨。 """ za, zb = z_for(1 - alpha), Z.inv_cdf(power) num = za * math.sqrt(psi) + zb * math.sqrt(max(0.0, psi - delta * delta)) return math.ceil(num * num / (delta * delta)) # --------------------------------------------------------------------------- # 5. 每道题采样多次：聚类标准误 # --------------------------------------------------------------------------- def naive_se(outcomes): """ 把所有样本当成独立的：outcomes 是 [[每道题的 k 次结果], ...]，SE = sqrt(p̂(1−p̂)/(n·k))。 【边界】同一道题的 k 次结果并不独立（难题 k 次都容易错），所以这个 SE **偏小**。只用来做对照。 """ flat = [x for item in outcomes for x in item] p = math.fsum(flat) / len(flat) return math.sqrt(p * (1 - p) / len(flat)) def cluster_se(outcomes): """ 聚类（按题目）标准误：先算每道题的平均正确率 m_i，再算 m_i 的均值的标准误 sd(m)/sqrt(n)。 【为什么】我们关心的是"换一批题目"时结果会怎样变，题目才是独立的单位。 Miller（2024，"Adding Error Bars to Evals"）建议评测 LLM 时都这样报告。 """ means = [math.fsum(item) / len(item) for item in outcomes] n = len(means) m = math.fsum(means) / n var = math.fsum((x - m) ** 2 for x in means) / (n - 1) return math.sqrt(var / n) # --------------------------------------------------------------------------- # 6. 模拟的提示实验世界 # --------------------------------------------------------------------------- @dataclass class ItemBank: """ 一个"题目总体"：每道题有一个难度 d ~ N(0, sd_difficulty²)；每个提示在每道题上还有一点自己的"手感" ε ~ N(0, sd_interaction²)（同一个提示在不同题上的表现不完全由难度决定）。 【贪心解码】提示（能力 θ）答对第 i 题 ⟺ θ + ε_i − d_i > 0。于是它在总体上的真实准确率 = P(θ + ε − d > 0) = Φ(θ / sqrt(sd_d² + sd_ε²))——可以精确算出来（true_accuracy）。 【采样解码】温度 > 0 时，每次采样以概率 σ(θ + ε_i − d_i) 答对（σ 是 sigmoid），同一道题的多次采样共享 ε_i、d_i。 【为什么这样设计】 - 同一道题上，两个提示的结果都受 d_i 影响 → 正相关 → 配对检验更灵敏（02 号脚本）； - ε 让两个提示在一部分题上"意见不同" → McNemar 检验有不一致的题可用； - 难度差别大 → 每道题的正确率很不一样 → 多次采样不能代替更多的题目（04 号脚本）。 """ sd_difficulty: float = 1.5 sd_interaction: float = 0.7 @property def scale(self): return math.sqrt(self.sd_difficulty ** 2 + self.sd_interaction ** 2) def theta_for(self, acc): """真实准确率为 acc 的提示的能力 θ（true_accuracy 的反函数）。""" return self.scale * Z.inv_cdf(acc) def true_accuracy(self, theta): """能力 θ 的提示在题目总体上的真实（贪心）准确率。""" return Z.cdf(theta / self.scale) def draw_items(self, n, rng): """抽 n 道题（只需要它们的难度）。""" return [rng.gauss(0.0, self.sd_difficulty) for _ in range(n)] def greedy(self, theta, items, rng): """ 一个提示在这些题上的贪心结果（0/1 列表）。每次调用都为这个提示重新抽一份 ε—— 所以对两个不同的提示分别调用，就得到了两个"手感"不同、但面对同样难度的提示。 """ return [1 if theta + rng.gauss(0.0, self.sd_interaction) - d > 0 else 0 for d in items] def sampled(self, theta, items, k, rng): """温度 > 0：每道题采样 k 次，返回 [[k 个 0/1], ...]。""" out = [] for d in items: p = 1 / (1 + math.exp(-(theta + rng.gauss(0.0, self.sd_interaction) - d))) out.append([1 if rng.random() < p else 0 for _ in range(k)]) return out if __name__ == "__main__": print(f" Wilson(7/10) = {tuple(round(x, 3) for x in wilson_ci(7, 10))}，Wald(7/10) = {tuple(round(x, 3) for x in wald_ci(7, 10))}") print(f" McNemar：A 对 B 错 2 题、A 错 B 对 10 题 → p = {mcnemar_exact([1]*2 + [0]*10, [0]*2 + [1]*10)[2]:.4f}") print(f" ±2 个百分点（p = 0.8）需要 {n_for_ci_halfwidth(0.8, 0.02)} 道题；" f"发现 0.70 → 0.73 的差别：不配对每组 {n_unpaired(0.70, 0.73)} 道，配对（ψ = 0.2）{n_paired_mcnemar(0.03, 0.2)} 道") bank = ItemBank() rng = random.Random(0) th = bank.theta_for(0.7) items = bank.draw_items(20000, rng) print(f" ItemBank：θ = {th:.3f} 的真实准确率 {bank.true_accuracy(th):.3f}，2 万道题上的观测 {math.fsum(bank.greedy(th, items, rng)) / 20000:.3f}")