13.1 把用户稳稳地分成了两组。现在问题变成：两组差多少才算真的差？ 这是整章最容易出错、也最容易被运营当成玄学的一节。很多团队的做法是「实验跑一周，看哪个数高就上哪个」——这几乎必然会把噪声当成信号。本节要把「怎么判」变成能算、能复现的流程。
做法上有一条硬约束：不引入 scipy 或 statsmodels。正态 CDF、两比例 z 检验、Welch t 检验、样本量公式全部手写，随机数用与 13.1 同一套的确定性 PRNG。原因很实际——这些公式是实验流程的地基，把它们当黑箱调库，就永远说不清「p 值到底是怎么来的」，也就看不出下面那些陷阱。
假设检验的四个数
一次两样本检验，绕不开四个数：
- α（第一类错误率）：原假设为真时，却判「有差异」的概率。本套取 50‰（0.05）。它控制「假阳性」。
- power（功效）：真有差异时，能检出来的概率。本套取 800‰（0.8）。它控制「漏报」。
- p 值：在「两组真的一样」的前提下，观测到当前或更极端差异的概率。它不是「有差异的概率」——这是最常被误读的一点。
- 样本量 n：由 α、power 和「想检出的差异」共同决定。n 不够，再好的设计也检不出真实效应。
手写实现：不依赖统计库
下面这段是 S 层统计函数库，把上面四个数需要的公式全部手写：正态 CDF 用 Abramowitz & Stegun 7.1.26 的 erf 近似，t 分布 CDF 用正则化不完全 beta 的连分式，样本量用标准公式，随机数用向量化的 mulberry32（与 sim/n13_split.py 同算法）：
# file: sim/n13_stats.py
# tier: S
# verify: python3 sim/n13_stats.py
"""手写统计函数库（13.2）：两比例 z 检验、Welch t 检验、样本量公式、向量化 mulberry32。
不引入 scipy/statsmodels：正态 CDF 用 A&S 7.1.26 的 erf 近似，t 分布 CDF 用正则化不完全 beta 连分式。
被 sim/n13_pitfalls.py（陷阱仿真）与 sim/n13_stats_bench.py（B 基准）复用。
"""
from __future__ import annotations
import math
import numpy as np
def norm_cdf(x):
"""标准正态 CDF：erf 用 A&S 7.1.26 近似（绝对误差 < 1.5e-7）。"""
x = np.asarray(x, dtype=float)
ax = np.abs(x) / math.sqrt(2.0)
t = 1.0 / (1.0 + 0.3275911 * ax)
poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741
+ t * (-1.453152027 + t * 1.061405429))))
return 0.5 * (1.0 + np.sign(x) * (1.0 - poly * np.exp(-ax * ax)))
def norm_ppf(p: float) -> float:
"""标准正态分位数：在 CDF 上二分（手写，不查表）。"""
lo, hi = -10.0, 10.0
for _ in range(200):
mid = 0.5 * (lo + hi)
lo, hi = (mid, hi) if float(norm_cdf(mid)) < p else (lo, mid)
return 0.5 * (lo + hi)
def two_prop_z(c1, n1, c2, n2):
"""两比例 z 检验（合并方差），返回 z 与双侧 p 值。"""
c1, c2 = np.asarray(c1, float), np.asarray(c2, float)
p = (c1 + c2) / (n1 + n2)
se = np.sqrt(p * (1.0 - p) * (1.0 / n1 + 1.0 / n2))
z = np.divide(c1 / n1 - c2 / n2, se, out=np.zeros_like(se), where=se > 0)
return z, 2.0 * (1.0 - norm_cdf(np.abs(z)))
def _betacf(a: float, b: float, x: float) -> float:
"""正则化不完全 beta 的连分式（Numerical Recipes 风格）。"""
MAXIT, EPS, FPMIN = 300, 3e-14, 1e-300
qab, qap, qam = a + b, a + 1.0, a - 1.0
c, d = 1.0, 1.0 - qab * x / qap
d = 1.0 / (d if abs(d) > FPMIN else FPMIN)
h = d
for m in range(1, MAXIT + 1):
m2 = 2 * m
for aa in (m * (b - m) * x / ((qam + m2) * (a + m2)),
-(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))):
d = 1.0 + aa * d
c = 1.0 + aa / c
d = 1.0 / (d if abs(d) > FPMIN else FPMIN)
c = c if abs(c) > FPMIN else FPMIN
h *= d * c
if abs(d * c - 1.0) < EPS:
break
return h
def _betai(a: float, b: float, x: float) -> float:
if x <= 0.0:
return 0.0
if x >= 1.0:
return 1.0
bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
+ a * math.log(x) + b * math.log(1.0 - x))
if x < (a + 1.0) / (a + b + 2.0):
return bt * _betacf(a, b, x) / a
return 1.0 - bt * _betacf(b, a, 1.0 - x) / b
def student_t_cdf(t: float, df: float) -> float:
"""学生 t 分布 CDF（df 自由度）。"""
ib = _betai(df / 2.0, 0.5, df / (df + t * t))
return 1.0 - 0.5 * ib if t > 0 else 0.5 * ib
def welch_t(m1, s1, n1, m2, s2, n2):
"""Welch 两样本 t 检验，返回 t、df、双侧 p。"""
v1, v2 = s1 * s1 / n1, s2 * s2 / n2
t = (m1 - m2) / math.sqrt(v1 + v2)
df = (v1 + v2) ** 2 / (v1 * v1 / (n1 - 1) + v2 * v2 / (n2 - 1))
return t, df, 2.0 * (1.0 - student_t_cdf(abs(t), df))
def sample_size_prop(p1: float, p2: float, alpha: float, power: float) -> int:
"""两比例每组样本量：n = (z_{1-α/2}+z_β)^2 (p1(1-p1)+p2(1-p2)) / (p1-p2)^2。"""
return math.ceil((norm_ppf(1 - alpha / 2) + norm_ppf(power)) ** 2
* (p1 * (1 - p1) + p2 * (1 - p2)) / (p1 - p2) ** 2)
def sample_size_mean(sigma: float, delta: float, alpha: float, power: float) -> int:
"""两均值每组样本量：n = 2σ²(z_{1-α/2}+z_β)^2 / δ²。"""
return math.ceil(2.0 * sigma * sigma * (norm_ppf(1 - alpha / 2) + norm_ppf(power)) ** 2 / delta ** 2)
def mulberry32_stream(seed: int, n: int) -> np.ndarray:
"""生成 n 个 [0,1) 均匀数：把 mulberry32 状态序列向量化，与标量版逐位一致。"""
state = (np.uint64(seed & 0xFFFFFFFF)
+ np.uint64(0x6D2B79F5) * np.arange(1, n + 1, dtype=np.uint64)) & np.uint64(0xFFFFFFFF)
t = ((state ^ (state >> np.uint64(15))) * (state | np.uint64(1))) & np.uint64(0xFFFFFFFF)
inner = ((t ^ (t >> np.uint64(7))) * (t | np.uint64(61))) & np.uint64(0xFFFFFFFF)
t = (t ^ (t + inner)) & np.uint64(0xFFFFFFFF)
return ((t ^ (t >> np.uint64(14))) & np.uint64(0xFFFFFFFF)).astype(np.float64) / 4294967296.0
if __name__ == "__main__":
print("norm_cdf(1.96) =", round(float(norm_cdf(1.96)), 6))
print("norm_ppf(0.975) =", round(norm_ppf(0.975), 4))
print("t_cdf(1.96, 1e6) =", round(student_t_cdf(1.96, 1e6), 6))
print("sample_size_prop(0.08,0.096,0.05,0.8) =", sample_size_prop(0.08, 0.096, 0.05, 0.8))
print("first 3 uniforms(seed=20261012) =", [round(float(v), 9) for v in mulberry32_stream(20261012, 3)])
本机自测输出：norm_cdf(1.96) = 0.975002（解析值是 0.975），norm_ppf(0.975) = 1.96，t_cdf(1.96, 1e6) = 0.975002（自由度极大时 t 分布收敛到正态，对得上），sample_size_prop(0.08, 0.096, 0.05, 0.8) = 4918。手写实现的第一步，是拿已知值校准——正态 CDF 在 1.96 处必须是 0.975，这是它没写错的证据。
样本量：想检出 80‰→96‰ 要多少人
把主实验的假设代进公式：对照组 80‰、实验组 96‰、α=0.05、power=0.8。库里的 sample_size_prop 算出每组 4918 人。这意味着：
- 每组少于 4918 人，实验大概率检不出这个效应（功效不足，容易漏报）。
- 这个数随效应变小而平方级增长。差异从 16‰ 减半到 8‰，每组样本量要涨到约 4 倍（
sim/n13_pitfalls.py的输出里，均值型指标sampleSizeMeanPerArm = 14128，同样是被方差和效应一起推高的）。
样本量必须在开实验之前算好，而不是「跑一周看看」。跑完再看数据决定「还要不要继续收」，就掉进了下面的 peeking 陷阱。
陷阱一：peeking（边跑边看）
真实运营最常见的动作是：实验上线后每天盯报表，一旦某个指标显著了，就宣布成功、停止实验。这个动作有个学名——peeking（偷看），它会把假阳性率推高。sim/n13_pitfalls.py 模拟了「在 10 个中间点各查一次，任一次显著就算成功」：
| 指标 | 实测值 |
|---|---|
| 单次检验的第一类错误率 | 0.046 |
| 偷看 10 次后的「假阳性率」 | 0.1885 |
| 膨胀倍数 | 4.098 |
同一批数据、同一个 α=0.05，只因为多看了 9 次，假阳性率从 4.6% 涨到 18.9%——膨胀约 4 倍。 如果 10 次偷看是独立的，理论上界是 1 − 0.95^10 = 0.401；实际 0.1885 更低，因为中间点的数据是嵌套的（后一次偷看包含了前一次的样本），高度相关，所以没到理论上界。但即便如此，4 倍膨胀已经足以让「每天看一眼」的流程变成假阳性制造机。
正确做法是预先固定样本量或固定实验时长，中途不看，到点再看。要做序贯检验（sequential test）也行，但必须用专门的 α 消耗函数，不能拿普通 z 检验反复查。
陷阱二：多重比较
peeking 是「同一实验看多次」，多重比较是「同时跑多个实验，挑一个显著的报」。sim/n13_pitfalls.py 模拟 20 个独立的 H0 检验，看「至少一个显著」的概率：
| 量 | 值 |
|---|---|
| 蒙特卡洛实测 FWER | 0.6358 |
解析值 1 − 0.95^20 | 0.6415 |
只要跑 20 个「其实都没用」的改动，至少有 64% 的概率能凑出一个「显著」的。 实测 0.6358 与解析 0.6415 吻合，说明独立假设基本成立。对策是多重比较校正：Bonferroni 把显著性阈值从 α 改成 α/m（20 个检验就用 0.0025），或者用 FDR 控制。「我跑了 20 个实验，终于找到一个有效的」不是惊喜，是数学必然。
陷阱三：SRM（样本比例失衡）
SRM（Sample Ratio Mismatch）是分流的体检指标：分的时候说好 50/50，跑完发现是 51/49。这通常意味着分流本身有 bug（比如某个端的分桶逻辑不一致、某个变体加载失败导致用户没进组），此时任何指标比较都不可信。sim/n13_pitfalls.py 对观测到的 [51000, 49000] 做卡方检验：
期望 [50000, 50000]，观测 [51000, 49000]
chi2 = (1000²/50000) × 2 = 40.0
df = 1，p ≈ 5×10⁻¹⁰ ≪ 0.05 → SRM 成立
卡方 40、p 值远小于 0.05，判定 SRM 成立。SRM 检验应该是每次实验的第一个动作——先确认两组人数对不对，再看指标。人数都不对，后面的 p 值没有意义。
陷阱四：辛普森悖论
最反直觉的一个：分层看每组都是 A 好，合起来看却是 B 好。sim/n13_pitfalls.py 构造了两层用户：
| 分层 | 对照组（‰） | 实验组（‰） | 谁好 |
|---|---|---|---|
| 新玩家 | 120 | 100 | 对照 |
| 老玩家 | 60 | 50 | 对照 |
| 聚合 | 72 | 90 | 实验 |
分层里对照组两层都赢，一聚合实验组反而赢了。原因是分层权重变了：新玩家整体转化率高（120‰），而实验组里新玩家占了绝大多数（4000/5000）。聚合时高转化的新玩家把实验组的平均值拉了上去——聚合掩盖了分层内的真实方向。这就是为什么实验必须事先定义分层（按新老、按渠道、按付费档），并在分层内看结论；「一锅烩」的聚合数字经常是错的。
连续指标与 t 检验
转化率是二值的，用 z 检验；ARPU 这类连续指标用 Welch t 检验。sim/n13_pitfalls.py 对一个「ARPU 从 300 分涨到 330 分」的结果做检验（标准差 900/910，每组 5000 人）：
t = −1.6574，df ≈ 9996.78，p = 0.097463
p = 0.097 > 0.05，不显著。 均值看着涨了 10%，但 ARPU 的标准差（900）远大于均值差（30），在 5000 人的样本量下这点差异完全可能是噪声。这条与转化率实验形成对照：指标方差越大，同样的人数越检不出差异。连续指标的样本量公式里那个 σ² 就是这么进来的。
B 基准：把模拟值和理论值对上
上面所有结论都建立在「手写公式是对的」这个前提上。下面这段是 B 层基准脚本（软告警不硬失败），用蒙特卡洛把实测值和理论值逐项对照：
# file: sim/n13_stats_bench.py
# tier: B
# verify: python3 sim/n13_stats_bench.py --config config/ --out out/
"""统计公式基准（B 层，软告警不硬失败）：把蒙特卡洛实测值与理论值对照。
1. 第一类错误率的实测覆盖率：H0 为真时拒绝率应 ≈ α。
2. p 值均匀性：H0 下 p 值应均匀，用 KS 风格的最大 CDF 偏差检验。
3. 功效：样本量公式算出的 n 下，实测功效应 ≈ 目标功效。
产出 out/n13_stats_bench.json。
"""
from __future__ import annotations
import argparse
import json
import os
import numpy as np
from n13_stats import mulberry32_stream, sample_size_prop, two_prop_z
def main() -> int:
here = os.path.dirname(os.path.abspath(__file__))
ap = argparse.ArgumentParser()
ap.add_argument("--config", default=os.path.join(here, "..", "config"))
ap.add_argument("--out", default=os.path.join(here, "..", "out"))
args = ap.parse_args()
os.makedirs(args.out, exist_ok=True)
with open(os.path.join(args.config, "n13_experiment.json"), encoding="utf-8") as fh:
cfg = json.load(fh)
seed, ex, mc = cfg["seed"], cfg["experiment"], cfg["montecarlo"]
alpha, target_power = ex["alphaPermille"] / 1000.0, ex["powerPermille"] / 1000.0
p_c, p_t = ex["controlRatePermille"] / 1000.0, ex["treatRatePermille"] / 1000.0
trials, n_arm, h0 = mc["trials"], mc["nPerArm"], mc["h0RatePermille"] / 1000.0
# 1) H0 下的 p 值与第一类错误率
u = mulberry32_stream(seed, trials * 2 * n_arm).reshape(trials, 2, n_arm)
_, pv = two_prop_z((u[:, 0] < h0).sum(axis=1), n_arm, (u[:, 1] < h0).sum(axis=1), n_arm)
emp_alpha = float((pv < alpha).mean())
# 2) p 值均匀性：经验 CDF 与均匀 CDF 的最大偏差（KS 风格；离散 p 值下比卡方更稳）
grid = np.linspace(0.01, 0.99, 99)
ks = float(np.max(np.abs(np.searchsorted(np.sort(pv), grid, side="right") / trials - grid)))
crit = 1.358 / np.sqrt(trials)
# 3) 功效：公式 n 下的实测功效
n_prop = sample_size_prop(p_c, p_t, alpha, target_power)
pt = min(trials, 2000)
ub = mulberry32_stream((seed ^ 0x2222) & 0xFFFFFFFF, pt * 2 * n_prop).reshape(pt, 2, n_prop)
_, pvb = two_prop_z((ub[:, 0] < p_c).sum(axis=1), n_prop, (ub[:, 1] < p_t).sum(axis=1), n_prop)
emp_power = float((pvb < alpha).mean())
checks = {
"type1": {"empirical": round(emp_alpha, 5), "nominal": alpha, "coverage": round(emp_alpha / alpha, 4),
"tol": 0.01, "pass": abs(emp_alpha - alpha) <= 0.01},
"pValueUniformity": {"ks": round(ks, 5), "critical": round(float(crit), 5), "pass": bool(ks < crit)},
"power": {"empirical": round(emp_power, 4), "target": target_power, "nPerArm": n_prop,
"pass": abs(emp_power - target_power) <= 0.05},
}
report = {"version": 1, "kind": "baseline", "trials": trials, "nPerArm": n_arm, "checks": checks,
"pValueBins": np.histogram(pv, bins=10, range=(0.0, 1.0))[0].astype(int).tolist(),
"pass": all(c["pass"] for c in checks.values())}
with open(os.path.join(args.out, "n13_stats_bench.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
print(json.dumps(report, ensure_ascii=False))
return 0
if __name__ == "__main__":
raise SystemExit(main())
本机实跑 out/n13_stats_bench.json：
| 对照项 | 实测值 | 理论值 | 判定 |
|---|---|---|---|
| 第一类错误率（覆盖率） | 0.046 | 0.05（覆盖率 0.92） | 通过 |
| p 值均匀性（KS 统计量） | 0.01375 | 临界值 0.02147 | 通过 |
| 功效（n=4918/组） | 0.8 | 目标 0.8 | 通过 |
三项全部通过，pass = true。三条对照分别验证了三件事：α 控制住了（实测 0.046 与名义 0.05 相差在容差 0.01 内）、p 值分布是均匀的（KS 统计量 0.01375 小于 5% 临界值 0.02147，说明「p 值在 H0 下均匀分布」这一前提成立）、样本量公式给出的 n 真能达到目标功效（实测 0.8 精确命中目标 0.8）。B 基准的意义就在这里：它不是再算一遍，而是验证前面所有结论所依赖的公式本身没写错。
小结
本节把「怎么判」变成了可跑流程：手写 z 检验 / t 检验 / 样本量公式，用 B 基准验证公式无误，再量化四个陷阱——peeking 把假阳性率抬 4 倍、20 次多重比较有 64% 概率凑出显著、SRM 靠卡方一票否决、辛普森悖论让聚合方向反转。有了这些，「实验跑出来了」才不等于「结论可信」。
但即便统计上显著，还有一个问题没解决：显著就等于该上线吗？ 下一节把 p 值、置信区间和业务门槛合起来，给出一个能拍板的决策规则。
阅读导航：上一节：13.1 实验设计与分流 · 下一节：13.3 从实验到决策 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。