5.2 方差与收敛
5.1 只问了估计量「对不对」——无偏性是均值层面的性质，说的是 $N \to \infty$ 时它落在哪儿。本节问的是「多准」：给定 $N$ 个样本，误差有多大，多喂样本能把它压多快。这是全卷判据从「相等」换成「收敛」之后的第一件实事。实时 GI 那篇靠时间上的重投影把噪声糊掉（见 实时全局光照 ），离线光追没有帧与帧之间的免费信息，只能把样本堆到误差可接受——堆多少，取决于误差怎么随 $N$ 缩。
从「相等」到「收敛」
固定 $N$ 时，$\langle I \rangle$ 是一个随机变量，不是一个数。它的期望是 $I$，方差则完全由被积函数与 pdf 的比值决定：
$$ \mathrm{Var},\langle I \rangle = \frac{1}{N}\left(\int_\Omega \frac{f(x)^2}{p(x)},dx - I^2\right) \equiv \frac{\sigma^2}{N} $$
括号里那个量与 $N$ 无关，只与 $p$ 有关；于是整条收敛律就是一句话：误差的均方根按 $N^{-1/2}$ 缩。想要精度翻倍，样本数得乘 4；想要精度提高 10 倍，样本数得乘 100。这是离线渲染里最贵的那个常数，也是本章所有优化的靶子——5.3 讲的分层与低差异序列，改的不是这条指数，而是括号里那个 $\sigma^2$。
有一点必须说清：$N^{-1/2}$ 描述的是方差的期望行为，不是某一次运行的表现。单跑一次，误差是随机的；只有把很多条独立子流的估计值聚在一起，方差本身才成为一个可测量的量。下面的脚本就是这么做的——每档 $N$ 跑 2048 条独立子流（起点间隔 1 << 40），由这 2048 个估计值反推方差。
逐档实测
真值仍然是 $\int_0^1 x^2,dx = 1/3$（Q16.16 下 $65536/3 = 21845.33$）。同一次内层循环里同时累两个估计量：无偏版用 $p(x)=1$，有偏版把 pdf 写成常数 $1.1$，因此两个序列共享同一批样本，比值不受统计噪声影响。
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_05_2_variance.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_05_2_variance.py
"""5.2 方差与收敛：相对误差按 N^(-1/2) 下降，斜率用 log-log 最小二乘量出来。
<I>_N = (1/N) Σ f(x_i)/p(x_i) 的方差是 σ²/N，故误差以 1/√N 收敛。三件事做成可断言
的实跑：(1) N = 1…4096 每档 2048 条独立子流，实测方差、相对标准差、相对 RMS 误差；
(2) 对 log N ~ log(相对 RMS 误差) 做最小二乘，无偏版斜率必须落在 -0.5 ± 0.03 内；
(3) 对照——pdf 写错成 1.1 后误差里多出一个不随 N 缩的常数，log-log 曲线在 N 大时被
压平，斜率从 -0.500 掉到 -0.27；而「相对标准差」那条仍是 -0.5：偏差不进方差，只进
总误差。SVG 用纯 xml.etree 手写，写完读回解析、核对圆点数。
产物：config/cg_trace_05_2_variance.json（逐档整数和 + 拟合结果，P 层读它复现）
out/cg_trace_05_2_convergence.svg（log-log 双曲线 + 两条拟合直线）
"""
import json
import math
import os
import sys
import xml.etree.ElementTree as ET
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
ONE = cglib.ONE
SEED0 = 0x00C0FFEE_12345678
STREAM_STRIDE = 1 << 40
N_LIST = (1, 4, 16, 64, 256, 1024, 4096)
TRIALS = 2048
P_BIAS = 72090 # 错的 pdf：1.1 * 65536 = 72090.6 -> 72090
TRUTH = ONE / 3.0 # ∫_0^1 x² dx = 1/3，Q16.16
W, H = 760, 440
COLORS = {"unbiased": "#1f77b4", "biased": "#d62728"}
Y_TICKS = ((0.0, "1"), (-0.5, "0.32"), (-1.0, "0.1"), (-1.5, "0.032"), (-2.0, "0.01"))
LEGEND = (("unbiased", "无偏 pdf = 1"), ("biased", "有偏 pdf = 1.1"))
def stream(index):
"""第 index 条独立子流的起点（index 从 1 开始）。"""
return (SEED0 + index * STREAM_STRIDE) & cgtrace.RNG_MASK
def sweep():
"""每档 N 跑 TRIALS 条独立子流；同一轮循环里同时累无偏与有偏两个估计量。"""
out = {}
for n in N_LIST:
su = suq = sb = sbq = 0
for t in range(TRIALS):
state = stream(t + 1)
au = ab = 0
for _ in range(n):
state, q = cgtrace.rng_q16(state)
f = cglib.fmul(q, q)
au += f
ab += cglib.fdiv(f, P_BIAS)
u, b = au // n, ab // n
su += u
suq += u * u
sb += b
sbq += b * b
out[n] = {"sum": su, "sumSq": suq, "biasSum": sb, "biasSumSq": sbq}
return out
def moments(sumv, sumsq, trials):
"""由整数和还原均值与样本方差：Var = (T·Σx² − (Σx)²) / (T(T−1))，分子全程整数。"""
return sumv / trials, (sumsq * trials - sumv * sumv) / (trials * (trials - 1))
def rows(data, ksum, ksq):
"""逐档 N 的相对标准差、相对 RMS 误差、相对偏差与「偏离/标准误」。"""
out = []
for n in N_LIST:
mu, var = moments(data[n][ksum], data[n][ksq], TRIALS)
sd = math.sqrt(var)
out.append({"n": n, "meanQ16": int(mu + 0.5), "stdQ16": round(sd, 4),
"relStd": round(sd / TRUTH, 6),
"relRms": round(math.sqrt(var + (mu - TRUTH) ** 2) / TRUTH, 6),
"relBias": round((mu - TRUTH) / TRUTH, 6),
"sigma": round((mu - TRUTH) / math.sqrt(var / TRIALS), 3)})
return out
def fit(rs, key):
"""最小二乘直线 log10(y) = a·log10(N) + b；返回斜率、截距与最大残差。"""
xs = [math.log10(r["n"]) for r in rs]
ys = [math.log10(r[key]) for r in rs]
mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
sxx = sum((x - mx) ** 2 for x in xs)
a = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sxx
b = my - a * mx
return a, b, max(abs(y - (a * x + b)) for x, y in zip(xs, ys))
def svg(series, fits, path):
"""纯 xml.etree 手写 log-log 折线图；返回落盘后解析出的圆点数。"""
l, r, t, b = 82, 24, 38, 56
x0, x1, y0, y1 = -0.25, 3.9, -2.05, 0.15
def px(x):
return round(l + (x - x0) / (x1 - x0) * (W - l - r), 2)
def py(y):
return round(H - b - (y - y0) / (y1 - y0) * (H - t - b), 2)
def el(tag, **kw):
# 关键字参数里的下划线还原成 SVG 属性名的连字符：stroke_width -> stroke-width
return ET.SubElement(root, tag, {k.replace("_", "-"): str(v) for k, v in kw.items()})
root = ET.Element("svg", {"xmlns": "http://www.w3.org/2000/svg", "width": str(W),
"height": str(H), "viewBox": f"0 0 {W} {H}"})
el("rect", x=0, y=0, width=W, height=H, fill="#ffffff")
for k, lab in Y_TICKS:
el("line", x1=l, y1=py(k), x2=px(x1), y2=py(k), stroke="#e6e6e6", stroke_width=1)
el("text", x=l - 8, y=py(k) + 4, fill="#555", font_size=12, text_anchor="end").text = lab
for k in (0, 1, 2, 3):
el("line", x1=px(k), y1=py(y0), x2=px(k), y2=py(y1), stroke="#e6e6e6", stroke_width=1)
el("text", x=px(k), y=py(y0) + 18, fill="#555", font_size=12,
text_anchor="middle").text = f"1e{k}"
el("line", x1=l, y1=py(y0), x2=px(x1), y2=py(y0), stroke="#333", stroke_width=1.5)
el("line", x1=l, y1=py(y0), x2=l, y2=py(y1), stroke="#333", stroke_width=1.5)
el("text", x=(l + W - r) / 2, y=20, fill="#111", font_size=14,
text_anchor="middle").text = "相对 RMS 误差 vs 样本数 N（双对数）"
for name in ("unbiased", "biased"):
pts = " ".join(f"{px(math.log10(pt['n']))},{py(math.log10(pt['relRms']))}"
for pt in series[name])
el("polyline", points=pts, fill="none", stroke=COLORS[name], stroke_width=2)
a, bb, _ = fits[name]["relRms"]
el("line", x1=px(x0), y1=py(a * x0 + bb), x2=px(x1), y2=py(a * x1 + bb),
stroke=COLORS[name], stroke_width=1, stroke_dasharray="6 4")
for pt in series[name]:
el("circle", cx=px(math.log10(pt["n"])), cy=py(math.log10(pt["relRms"])), r=3.5,
fill=COLORS[name])
for i, (name, lab) in enumerate(LEGEND):
y = t + 12 + i * 18
el("line", x1=W - 214, y1=y, x2=W - 190, y2=y, stroke=COLORS[name], stroke_width=3)
el("text", x=W - 184, y=y + 4, fill="#333", font_size=12).text = lab
ET.ElementTree(root).write(path, encoding="utf-8", xml_declaration=True)
return sum(1 for _ in ET.parse(path).getroot().iter("{http://www.w3.org/2000/svg}circle"))
def main():
os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)
data = sweep()
series = {"unbiased": rows(data, "sum", "sumSq"), "biased": rows(data, "biasSum", "biasSumSq")}
fits = {name: {k: fit(rs, k) for k in ("relRms", "relStd")} for name, rs in series.items()}
svg_path = os.path.join(ROOT, "out/cg_trace_05_2_convergence.svg")
circles = svg(series, fits, svg_path)
report = {"seed": str(SEED0), "streamStride": STREAM_STRIDE, "trials": TRIALS,
"nList": list(N_LIST), "pBias": P_BIAS, "exactQ16": TRUTH, "circles": circles,
"sums": {str(n): {k: str(v) for k, v in data[n].items()} for n in N_LIST},
"rows": series,
"fit": {name: {k: [round(v, 8) for v in fits[name][k]] for k in ("relRms", "relStd")}
for name in fits}}
with open(os.path.join(ROOT, "config/cg_trace_05_2_variance.json"), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
print("== 1. 逐档实测：每档 2048 条独立子流，真值 1/3 = 0.333333")
for r in series["unbiased"]:
print(f" N = {r['n']:>4}：σ/真值 = {r['relStd']:.6f}｜RMS/真值 = {r['relRms']:.6f}"
f"｜相对偏差 = {r['relBias']:+.6f}｜{r['sigma']:+.3f}σ")
print("== 2. log-log 最小二乘斜率（-0.5 即 1/√N 定律，容差 ±0.03）")
for name in ("unbiased", "biased"):
for k in ("relRms", "relStd"):
a, b, resid = fits[name][k]
print(f" {name:>8} · {k:<6}：斜率 = {a:+.4f}｜截距 = {b:+.4f}｜最大残差 = {resid:.4f}")
print("== 3. 有偏对照：常数偏差把总误差曲线压平")
for r in series["biased"]:
print(f" N = {r['n']:>4}：相对偏差 = {r['relBias']:+.6f}｜RMS/真值 = {r['relRms']:.6f}"
f"｜σ/真值 = {r['relStd']:.6f}")
print(f"== 4. SVG 落盘并解析回读：{circles} 个数据点圆")
return 0
if __name__ == "__main__":
sys.exit(main())
每档 $N$ 的相对标准差、相对 RMS 误差与相对偏差（真值 0.333333）：
| $N$ | $\sigma$/真值 | RMS/真值 | 相对偏差 | 偏离/$\sigma$ |
|---|---|---|---|---|
| 1 | 0.884147 | 0.884147 | −0.000124 | −0.006σ |
| 4 | 0.443527 | 0.443574 | −0.006420 | −0.655σ |
| 16 | 0.222989 | 0.223024 | −0.003947 | −0.801σ |
| 64 | 0.109941 | 0.109974 | −0.002692 | −1.108σ |
| 256 | 0.057667 | 0.057669 | −0.000367 | −0.288σ |
| 1024 | 0.028301 | 0.028301 | −0.000020 | −0.032σ |
| 4096 | 0.013867 | 0.013867 | +0.000063 | +0.204σ |
$N$ 从 1 涨到 4096（4096 倍），$\sigma$ 从 0.884147 掉到 0.013867（63.8 倍）——$\sqrt{4096} = 64$。这就是 $N^{-1/2}$ 的手感：样本数乘 4096，误差只除 64。同时注意最后一列：无偏版的偏离始终在 ±1.2σ 内游走，没有系统性方向，这正是无偏的统计表现。
脚本后半段把这件事做成了可断言的数，并给出有偏对照：
== 2. log-log 最小二乘斜率（-0.5 即 1/√N 定律，容差 ±0.03）
unbiased · relRms：斜率 = -0.4978｜截距 = -0.0527｜最大残差 = 0.0124
unbiased · relStd：斜率 = -0.4978｜截距 = -0.0527｜最大残差 = 0.0124
biased · relRms：斜率 = -0.2638｜截距 = -0.2418｜最大残差 = 0.1575
biased · relStd：斜率 = -0.4978｜截距 = -0.0941｜最大残差 = 0.0124
== 3. 有偏对照：常数偏差把总误差曲线压平
N = 1：相对偏差 = -0.091050｜RMS/真值 = 0.808907｜σ/真值 = 0.803766
N = 4：相对偏差 = -0.096775｜RMS/真值 = 0.414656｜σ/真值 = 0.403205
N = 16：相对偏差 = -0.094527｜RMS/真值 = 0.223672｜σ/真值 = 0.202716
N = 64：相对偏差 = -0.093386｜RMS/真值 = 0.136785｜σ/真值 = 0.099946
N = 256：相对偏差 = -0.091273｜RMS/真值 = 0.105257｜σ/真值 = 0.052424
N = 1024：相对偏差 = -0.090958｜RMS/真值 = 0.094526｜σ/真值 = 0.025727
N = 4096：相对偏差 = -0.090883｜RMS/真值 = 0.091753｜σ/真值 = 0.012606
log-log 斜率：把 1/√N 量成一个小数
「误差按 $N^{-1/2}$ 缩」等价于说：把 $\log N$ 当横轴、$\log(\text{误差})$ 当纵轴，数据点应该落在一条斜率 −0.5 的直线上。所以判据不需要眼睛看曲线，只需要最小二乘：
$$ \log_{10} E(N) = a \log_{10} N + b, \qquad a \overset{?}{=} -0.5 $$
实测 $a = -0.4978$，最大残差 0.0124。与 $-0.5$ 差 0.0022，落在 $\pm 0.03$ 的容差里。这个容差不是随手定的：$N = 1$ 那一档只有 2048 个估计值，$\sigma$ 的估计本身带约 $1/\sqrt{2 \times 2048} = 1.6%$ 的相对误差，折算到 $\log_{10}$ 轴上约 0.007，经过 7 点最小二乘后对斜率的贡献约 0.002。容差取 0.03 已经留了十几倍余量，跑偏到这个数以外就不可能是噪声。
脚本用纯 xml.etree 手写了这张双对数图（out/cg_trace_05_2_convergence.svg）：两条实测折线、两条虚线拟合直线、7 档刻度、14 个数据点圆。写完之后立刻用 ET.parse() 读回来数圆点，确认落盘的是一份能被解析的 SVG，而不是一段拼错的文本——这是「产物可验证」的最低门槛。
斜率跑偏说明什么
有偏版的 relRms 斜率是 −0.2638，明显偏离 $-0.5$。原因在表里：相对偏差从 $N=1$ 的 −0.091050 到 $N=4096$ 的 −0.090883，几乎不动。总误差是偏差与噪声的合成 $E = \sqrt{\text{bias}^2 + \sigma^2/N}$：小 $N$ 时噪声主导，曲线还贴着 $-0.5$；大 $N$ 时噪声被压到 0.012606，而偏差始终是 0.090909，于是曲线被这条水平的地板压平。
最能说明问题的是同一张表里的第三行：有偏版的 relStd 斜率仍是 −0.4978，与无偏版一模一样。偏差不进方差。 一个估计量的方差可以完美服从 $N^{-1/2}$，同时它的期望永远错着——5.1 的「多采样救不了偏差」在这里有了第二个侧面：方差曲线正常，不代表估计量正常。
于是 log-log 斜率成了一个诊断工具：
| 实测斜率 | 含义 |
|---|---|
| $\approx -0.5$ | 无偏、独立采样，一切正常 |
| $\approx 0$ | 误差被一条不随 $N$ 缩的地板钉住，几乎一定是偏差 |
| 明显大于 $-0.5$（如 $-0.26$） | 偏差与噪声在同一个量级，正在从噪声主导过渡到偏差主导 |
| 明显小于 $-0.5$（如 $-1$） | 样本不再独立，或分层数与样本数不匹配 |
最后一行是 5.3 的伏笔。分层采样（stratified sampling）把定义域切成 $S$ 块、每块各抽 $N/S$ 个样本，方差按 $S^{-1}$ 量级的因子下降；但如果 $N$ 不是 $S$ 的整数倍，或者写代码时把「层数」当成了「每层样本数」，某些层会被抽空、另一些层被抽重，此时各层的权重不再等于 $1/S$，误差曲线要么提前触底、要么比 $-0.5$ 更陡。斜率是这类实现错误的探针，而不是装饰。
P 层：整数和落盘，方差两侧各自还原
跨语言断言这里有个选择：落盘方差本身，还是落盘它的分子？方差是浮点量，两侧最后一位可能不同；$\sum x$ 与 $\sum x^2$ 是整数，只要 RNG 与定点乘除一致就必然逐位相同。TypeScript 侧因此重算全部 7 档 × 2048 条子流，累出 BigInt 整数和，再各自还原方差与斜率：
// snippet: 取自 lantern/core/src/lantern/cg_trace_05_2_variance.ts 的方差还原与斜率拟合（逐字摘录，完整文件另含 sweep 与 streamStart）
/**
* 由整数和还原均值与样本方差：Var = (T·Σx² − (Σx)²) / (T(T−1))。
* 分子分母都远在 2^53 之下（N = 1 时最大，约 1.6e15），故 Number() 转换不丢位。
*/
export function moments(
sumv: bigint, sumsq: bigint, trials: number,
): { mu: number; varr: number } {
const tb = BigInt(trials);
return {
mu: Number(sumv) / trials,
varr: Number(sumsq * tb - sumv * sumv) / Number(tb * (tb - 1n)),
};
}
/** 最小二乘直线 y = a·x + b；返回斜率、截距与最大残差。 */
export function fitSlope(xs: number[], ys: number[]): [number, number, number] {
const n = xs.length;
const mx = xs.reduce((s, x) => s + x, 0) / n;
const my = ys.reduce((s, y) => s + y, 0) / n;
let sxx = 0;
let sxy = 0;
for (let i = 0; i < n; i += 1) {
sxx += (xs[i] - mx) ** 2;
sxy += (xs[i] - mx) * (ys[i] - my);
}
const a = sxy / sxx;
const b = my - a * mx;
let resid = 0;
for (let i = 0; i < n; i += 1) resid = Math.max(resid, Math.abs(ys[i] - (a * xs[i] + b)));
return [a, b, resid];
}
八条断言全绿，其中三条是硬口径：7 档的 $\sum x$ 与 $\sum x^2$ 逐位相同（BigInt 精确整数比对，连有偏序列的四个和也一样）、斜率与最大残差复现到 $10^{-6}$ 内、SVG 解析出 14 个数据点圆（$2 \times 7$）。断言里还钉了两条反向判据：有偏版的总误差斜率必须与 $-0.5$ 差出 0.2 以上，而它的方差斜率必须仍在 $\pm 0.03$ 内——这条配对断言才是本节的核心，它把「偏差不进方差」从一句结论变成了一个会失败的检查。到这里，收敛速度已经是一个能被拟合出来的数了；下一节去改分母上那个 $\sigma^2$。
阅读导航：上一节：5.1 积分与估计 · 下一节：5.3 分层与低差异序列 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。