5.1 积分与估计
前四节的输出都是一个确定的像素值：求交给出 $t$，着色给出 $(r,g,b)$，深度限制把递归截断。从本节起这件事变了——要算的量变成一个积分，而积分的输出是一个带方差的估计值，判据也从"相等"换成"收敛"。实时 GI 那篇把这件事压在屏幕空间里做近似（见 实时全局光照 ）；本卷要把估计量本身拆到能逐位复现的粒度：估计量怎么写、无偏性靠什么成立、随机数的哪一位决定了它还能不能信。
从积分到估计量
着色方程里到处都是这种量，$f$ 是被积函数（沿某方向的入射辐亮度乘 BRDF 与余弦项），$\Omega$ 是定义域：
$$ I = \int_{\Omega} f(x), dx $$
除了极少数构型，它没有闭式解。蒙特卡洛的做法是把它换成一个估计量，$x_i$ 按概率密度 $p$ 采样：
$$ \langle I \rangle = \frac{1}{N}\sum_{i=1}^{N} \frac{f(x_i)}{p(x_i)} $$
这个式子的全部内容就是：在 $p$ 大的地方少放样本、在 $p$ 小的地方多放大权重。当 $p$ 恰好正比于 $f$ 时每一项都相等，方差为零；当 $p$ 是均匀分布时它就是朴素的平均值。估计量要无偏，只需一个条件：
$$ p(x) > 0 \quad \text{只要} \quad f(x) \neq 0 $$
一旦漏掉一块区域，那个区域对积分的贡献就永久丢失——不是"少算一点"，是"再也算不回来"。这就是偏差：$N$ 加到无穷，均值收敛到的是 $\int f,dx$ 的一个固定错误值。方差是随机项，多采样就能压小；偏差是系统项，多采样只会让它更精确地停在错的地方。本节的实跑就围绕这条分界展开：估 $\int_0^1 x^2,dx = 1/3$，正确写法用 $p(x)=1$，把 pdf 故意写错成常数 $1.1$ 就得到一个有偏版本。
SplitMix64：掩码必须显式写出来
采样之前先要有随机数。本卷全程用 SplitMix64——一个纯整数、全 64 位掩码的确定性函数，两侧（Python 的 int 与 TypeScript 的 BigInt）逐位一致。它的核心是四步：state = (state + GOLDEN) & MASK；z = state；两次 z = ((z ^ (z >> k)) * MIX) & MASK；最后 z = z ^ (z >> 31)。
这里有一个语言层面的陷阱：Python 的 int 是任意精度，BigInt 也是，它们都不会自动截断到 64 位。& MASK 必须手写，漏掉一个，乘法就把位宽撑到 191 位。脚本里专门留了一个漏掉掩码的对照版本 splitmix_unmasked，并用它和正确版本跑同一批状态：
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_05_1_estimator.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_05_1_estimator.py
"""5.1 积分与估计：无偏估计量、SplitMix64 的掩码、拒绝采样与取模的偏置。
I = ∫ f(x) dx 没有闭式解时，用 N 个样本的均值去估它：<I> = (1/N) Σ f(x_i)/p(x_i)。
三件事做成可断言的实跑：(1) 掩码——每步混合都要显式 `& RNG_MASK`，Python 的 int 与 JS
的 BigInt 都不自动截断，漏掉掩码后输出 z 涨到 191 位，`z >> 48` 返回 143 位整数，Q16.16
采样当场失效；(2) 无偏性——估 ∫_0^1 x² dx = 1/3，pdf 写错成常数 1.1 后均值收敛到
(1/3)/1.1 并停在那里；(3) 取模的偏置——`% n` 在 n 不整除值域宽度 R 时相对超配
(n − R mod n)/R，R = 2^64 时不可测，把随机数截窄到 8 位、4 位则升到 0.78%、12.5%。
产物：config/cg_trace_05_1_rng.json（采样器定义 + 逐位样本头 + 全量计数，P 层读它复现）
out/cg_trace_05_1_estimator.json（估计量表、掩码位宽表、直方图与卡方）
"""
import json
import math
import os
import sys
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 # 子流起点间隔，远大于任何一条子流的样本数
P_BIAS = 72090 # 错的 pdf：1.1 * 65536 = 72090.6 -> 72090
N_LIST = (1, 16, 64, 256, 1024, 4096)
TRIALS = 1024
M_HIST = 2_000_000
N_BINS = 3
RANGES = (4, 8, 12, 16, 64) # 64 表示「不截窄，直接用 64 位拒绝采样」
SIGMA_X2 = math.sqrt(1 / 5 - 1 / 9) # x 均匀时 x² 的标准差 = 0.298142…
def stream(index):
"""第 index 条独立子流的起点（index 从 1 开始）。"""
return (SEED0 + index * STREAM_STRIDE) & cgtrace.RNG_MASK
def splitmix_unmasked(state):
"""**故意漏掉掩码**的 SplitMix64：只保留 `state += GOLDEN`，三次乘法不截断。"""
state = state + cgtrace.RNG_GOLDEN
z = state
z = (z ^ (z >> 30)) * cgtrace.RNG_MIX1
z = (z ^ (z >> 27)) * cgtrace.RNG_MIX2
return state, z ^ (z >> 31)
def mask_table(seed, steps=6):
"""逐步对比「有掩码 / 无掩码」：位宽，以及 `z >> 48` 是否还在 Q16.16 的值域内。"""
rows, sm, su = [], seed, seed
for k in range(1, steps + 1):
sm, zm = cgtrace.rng_next(sm)
su, zu = splitmix_unmasked(su)
rows.append({"step": k, "maskedBits": zm.bit_length(),
"unmaskedBits": zu.bit_length(), "maskedQ16": zm >> 48,
"unmaskedQ16Bits": (zu >> 48).bit_length(),
"unmaskedInRange": (zu >> 48) < ONE})
return rows
def estimator_table():
"""无偏 vs 有偏：**同一批样本**喂给两个估计量，故比值不受统计噪声影响。"""
exact, rows = ONE / 3.0, []
for N in N_LIST:
su = sb = 0
for t in range(TRIALS):
state = stream(t + 1)
acc_u = acc_b = 0
for _ in range(N):
state, q = cgtrace.rng_q16(state)
f = cglib.fmul(q, q) # f(x) = x²，Q16.16
acc_u += f
acc_b += cglib.fdiv(f, P_BIAS) # 用错的 pdf
su += acc_u // N
sb += acc_b // N
mu_u, mu_b = su / TRIALS, sb / TRIALS
# 无偏均值的标准误（相对值）：σ(x²)/√(N·T) / 真值
rel_se = SIGMA_X2 / math.sqrt(N * TRIALS) / (1 / 3)
rows.append({"N": N, "unbiasedMeanQ16": int(mu_u + 0.5),
"unbiasedRelBias": round((mu_u - exact) / exact, 6),
"unbiasedStdErr": round(rel_se, 6),
"unbiasedSigma": round(((mu_u - exact) / exact) / rel_se, 3),
"biasedMeanQ16": int(mu_b + 0.5), "ratio": round(mu_b / mu_u, 6),
"expectedRatio": round(ONE / P_BIAS, 6)})
return rows
def histogram(state, n, m, bits=64):
"""bits = 64 时走 `rng_below` 拒绝采样；否则取低 bits 位再取模。"""
mask, counts, head = (1 << bits) - 1, [0] * n, []
for i in range(m):
if bits == 64:
state, v = cgtrace.rng_below(state, n)
else:
state, z = cgtrace.rng_next(state)
v = (z & mask) % n
counts[v] += 1
if i < 256:
head.append(v)
return counts, head, state
def chi2(counts, m):
"""卡方统计量：Σ (obs - exp)² / exp，exp = M/n。"""
exp = m / len(counts)
return sum((c - exp) ** 2 / exp for c in counts)
def chi2_theory(n, r, m):
"""理论卡方（非中心度）：值域 R = n*q + r 时 r 个大箱各占 (q+1)/R。"""
q, rem, u = *divmod(r, n), 1.0 / n
return m * (rem * ((q + 1) / r - u) ** 2 / u + (n - rem) * (q / r - u) ** 2 / u)
def chi2_crit(df, z=3.0902):
"""Wilson–Hilferty 近似：卡方分布的 99.9% 分位点。"""
return df * (1 - 2 / (9 * df) + z * math.sqrt(2 / (9 * df))) ** 3
def main():
os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)
mrows, est_rows = mask_table(SEED0), estimator_table()
head64, s = [], SEED0 # 前 8 个 64 位原始输出，P 层靠它比对
for _ in range(8):
s, z = cgtrace.rng_next(s)
head64.append(z)
# 五个值域宽度、同一分箱数、各一条独立子流
hist = {}
for i, bits in enumerate(RANGES):
counts, head, _ = histogram(stream(100 + i), N_BINS, M_HIST, bits)
hist[bits] = {"counts": counts, "head": head, "chi2": chi2(counts, M_HIST),
"theory": chi2_theory(N_BINS, 1 << bits, M_HIST),
"over": (N_BINS - ((1 << bits) % N_BINS)) / (1 << bits)}
crit = chi2_crit(N_BINS - 1)
# 可检出样本数：M* = 临界值 / (理论卡方每样本)。偏置小于噪声时要多少样本才够。
for b in RANGES:
h = hist[b]
h["detectN"] = None if h["theory"] <= 0 else int(crit * M_HIST / h["theory"])
# 一份 report 同时落盘到 config/ 与 out/：P 层读前者复现，读者看后者核对。
# 64 位整数超过 JS 的 2^53，一律以字符串落盘，TS 侧用 BigInt 还原。
report = {"seed": str(SEED0), "streamStride": STREAM_STRIDE, "M": M_HIST,
"n": N_BINS, "ranges": list(RANGES), "exactQ16": ONE / 3.0,
"maskTable": mrows, "head64": [str(v) for v in head64],
"counts": {str(b): hist[b]["counts"] for b in RANGES},
"heads3": {str(b): hist[b]["head"] for b in RANGES},
"estimator": {"nList": list(N_LIST), "trials": TRIALS, "pBias": P_BIAS},
"estimatorRows": est_rows,
"estimatorMeans": {"unbiased": [r["unbiasedMeanQ16"] for r in est_rows],
"biased": [r["biasedMeanQ16"] for r in est_rows]},
"hist": {"M": M_HIST, "n": N_BINS, "crit": round(crit, 4),
"rows": [{"bits": b, "counts": hist[b]["counts"],
"chi2": hist[b]["chi2"], "chi2Theory": hist[b]["theory"],
"over": hist[b]["over"], "detectN": hist[b]["detectN"]}
for b in RANGES]}}
for name in ("out/cg_trace_05_1_estimator.json", "config/cg_trace_05_1_rng.json"):
with open(os.path.join(ROOT, name), "w", encoding="utf-8") as fh:
json.dump(report, fh, ensure_ascii=False, indent=2)
print("== 1. 掩码：SplitMix64 的输出位宽与 `z >> 48`（同一条 64 位状态流）")
for r in mrows:
print(f" 步 {r['step']}：掩码 {r['maskedBits']} 位 / 无掩码 {r['unmaskedBits']} 位"
f"｜掩码 q16 = {r['maskedQ16']}｜无掩码 >>48 = {r['unmaskedQ16Bits']} 位"
f"｜在 [0,65536) 内 = {r['unmaskedInRange']}")
print(f" 掩码后前 3 个 64 位输出：{[str(v) for v in head64[:3]]}")
print("== 2. 估计量：∫_0^1 x² dx = 1/3 = 0.333333（每条试验 1024 次）")
for r in est_rows:
print(f" N = {r['N']:>4}：无偏 {r['unbiasedMeanQ16'] / ONE:.6f}"
f"（相对偏差 {r['unbiasedRelBias']:+.6f}，标准误 {r['unbiasedStdErr']:.6f}，"
f"{r['unbiasedSigma']:+.3f}σ）｜有偏 {r['biasedMeanQ16'] / ONE:.6f}"
f"｜有偏/无偏 = {r['ratio']:.6f}")
print(f" 期望比值 1/1.1 = {est_rows[0]['expectedRatio']:.6f}")
print("== 3. 直方图：n = 3，M = {:,}，期望每箱 {:.1f}".format(M_HIST, M_HIST / N_BINS))
print(f" χ² 的 99.9% 临界值（df = 2）= {crit:.3f}")
for b in RANGES:
h = hist[b]
tag = "拒绝采样（R = 2^64）" if b == 64 else f"取低 {b:>2} 位"
dn = "不可检出" if h["detectN"] is None else f"{h['detectN']:,}"
print(f" {tag}: 计数 {h['counts']}｜实测 χ² = {h['chi2']:9.3f}"
f"｜理论 χ² = {h['theory']:11.4g}｜相对超配 = {h['over'] * 100:.6f}%"
f"｜M* = {dn}")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一遍（9.98 秒，单线程），输出节选：
== 1. 掩码：SplitMix64 的输出位宽与 `z >> 48`（同一条 64 位状态流）
步 1：掩码 64 位 / 无掩码 191 位｜掩码 q16 = 54087｜无掩码 >>48 = 143 位｜在 [0,65536) 内 = False
步 2：掩码 64 位 / 无掩码 192 位｜掩码 q16 = 58932｜无掩码 >>48 = 144 位｜在 [0,65536) 内 = False
步 3：掩码 64 位 / 无掩码 192 位｜掩码 q16 = 59937｜无掩码 >>48 = 144 位｜在 [0,65536) 内 = False
步 4：掩码 64 位 / 无掩码 193 位｜掩码 q16 = 45597｜无掩码 >>48 = 145 位｜在 [0,65536) 内 = False
步 5：掩码 63 位 / 无掩码 193 位｜掩码 q16 = 26635｜无掩码 >>48 = 145 位｜在 [0,65536) 内 = False
步 6：掩码 64 位 / 无掩码 193 位｜掩码 q16 = 57670｜无掩码 >>48 = 145 位｜在 [0,65536) 内 = False
掩码后前 3 个 64 位输出：['15224344366058555842', '16587940323163529587', '16870881601195667877']
== 3. 直方图：n = 3，M = 2,000,000，期望每箱 666666.7
χ² 的 99.9% 临界值（df = 2）= 14.133
取低 4 位: 计数 [750139, 624774, 625087]｜实测 χ² = 15677.242｜理论 χ² = 1.562e+04｜相对超配 = 12.500000%｜M* = 1,808
取低 8 位: 计数 [672568, 662987, 664445]｜实测 χ² = 79.952｜理论 χ² = 61.04｜相对超配 = 0.781250%｜M* = 463,099
取低 12 位: 计数 [667342, 667220, 665438]｜实测 χ² = 3.408｜理论 χ² = 0.2384｜相对超配 = 0.048828%｜M* = 118,553,572
取低 16 位: 计数 [666572, 665535, 667893]｜实测 χ² = 4.190｜理论 χ² = 0.0009313｜相对超配 = 0.003052%｜M* = 30,349,714,502
拒绝采样（R = 2^64）: 计数 [666653, 667332, 666015]｜实测 χ² = 1.301｜理论 χ² = 0｜相对超配 = 0.000000%｜M* = 不可检出
读这张表要注意两列。掩码位数恒在 63–64：状态加一个常数，位数本来就可能掉到 63（第 5 步就是），这不是 bug。无掩码位数从 191 涨到 193：三次 64 位乘法每次贡献约 64 位，越乘越长。真正致命的是最后一列——无掩码时 z >> 48 返回的是 143–145 位的整数，而 rng_q16 约定它落在 $[0, 65536)$。于是 fmul(q, q) 的中间量直接炸到 $2^{290}$ 量级，Q16.16 的定点约定当场作废。掩码不是优化，是语义的一部分。
估计量：无偏与有偏的实跑
每条试验（共 1024 条）取一条独立子流——子流起点间隔 1 << 40，远大于任何一条子流的样本数；否则相邻试验的样本会重叠，“1024 条试验的均值"就退化成一次估计。Q16.16 定点下真值是 $65536/3 = 21845.33$：
| $N$ | 无偏均值 | 相对偏差 | 标准误 | 偏离/标准误 | 有偏均值 | 有偏/无偏 |
|---|---|---|---|---|---|---|
| 1 | 0.325241 | −0.024259 | 0.027951 | −0.868σ | 0.295670 | 0.909062 |
| 16 | 0.332352 | −0.002946 | 0.006988 | −0.422σ | 0.302124 | 0.909061 |
| 64 | 0.331161 | −0.006536 | 0.003494 | −1.871σ | 0.301041 | 0.909060 |
| 256 | 0.332504 | −0.002482 | 0.001747 | −1.421σ | 0.302261 | 0.909060 |
| 1024 | 0.333054 | −0.000826 | 0.000873 | −0.946σ | 0.302765 | 0.909061 |
| 4096 | 0.333328 | −0.000006 | 0.000437 | −0.013σ | 0.303024 | 0.909060 |
两列行为完全不同。无偏那列：相对偏差从 −0.024259 掉到 −0.000006，而偏离/标准误始终在 ±2σ 内游走——它没有系统性错误，剩下的全是噪声；标准误严格按 $1/\sqrt{N}$ 缩（0.027951 → 0.000437，$N$ 涨 4096 倍，标准误降 64 倍）。有偏那列：均值从 0.295670 爬到 0.303024，停在那里。加 4096 倍样本只让它更靠近 $0.303030 = (1/3)/1.1$，永远到不了 0.333333。同一批样本给出的比值是 0.909060–0.909062，与理论值 $65536/72090 = 0.909086$ 只差 2.6e-5——差值来自 Q16.16 的截断与逐样本取整，不是统计噪声。这就是"多采样救不了偏差"的量化版本。
拒绝采样为什么不能换成取模
rng_below(state, n) 要返回 $[0,n)$ 的整数。最省事的写法是 z % n，但这是错的。设随机数的值域宽度为 $R$（取低 $k$ 位就是 $R = 2^k$），把 $R$ 个等概率残数分给 $n$ 个箱：前 $R \bmod n$ 个箱各多分到一个残数，于是
$$ \text{相对超配} = \frac{n - (R \bmod n)}{R} $$
$R = 2^{64}$、$n = 3$ 时这个量是 $2^{-63}$ 量级——理论上存在，实际上测不出来。问题在于人们常常先把随机数截窄再用：z & 0xFF 之后取模，$R$ 就从 $2^{64}$ 掉到 $2^8$，相对超配立刻变成 0.78%。脚本用 3 分箱直方图 + 卡方把这条曲线量了出来，$M = 2{,}000{,}000$ 个样本、每条曲线一条独立子流：
| 值域 | 相对超配（理论） | 实测 $\chi^2$ | 理论 $\chi^2$ | 判定（临界 14.133） |
|---|---|---|---|---|
| $2^4$ | 12.500% | 15677.242 | 1.562e+04 | 拒绝 |
| $2^8$ | 0.781% | 79.952 | 61.04 | 拒绝 |
| $2^{12}$ | 0.0488% | 3.408 | 0.2384 | 接受 |
| $2^{16}$ | 0.00305% | 4.190 | 0.0009313 | 接受 |
| $2^{64}$（拒绝采样） | 0 | 1.301 | 0 | 接受 |
理论 $\chi^2$ 是"非中心度”：值域 $R = n q + r$ 时，$r$ 个大箱各占 $(q+1)/R$、其余占 $q/R$，把偏离代入 $\sum (p_k - 1/n)^2/(1/n)$ 再乘 $M$。实测与理论严格对应：$2^4$ 那一行 15677 对 15620，$2^8$ 是 79.952 对 61.04。$2^{64}$ 那一行理论值恰好为 0——$2^{64} \bmod 3 = 1$，但那"多分一个残数"摊到 $2^{64}$ 个残数上是 $2^{-63}$，落在浮点精度之外。最后一列 M* 是可检出样本数：$M^* = \text{临界值} \cdot M / \chi^2_{\text{理论}}$，$2^4$ 只要 1808 个样本就能测出偏置，$2^{16}$ 要 3.03e10 个，$2^{64}$ 则永远测不出来。结论不是"取模有害"，而是：取模的偏置与值域宽度成反比，所以危险不在取模，在你把值域截窄了。 拒绝采样（lim = MASK - (MASK % n)，超出 lim 就重抽）在 64 位值域上把偏置压到不可测，代价是期望重抽次数小于 2 次。第 7 章讲重要性采样时会再遇到这一取舍。
P 层：两侧逐位一致
TypeScript 侧读同一份 config（seed 与 64 位输出都以字符串落盘——它们超过 Number 的 2^53 精度），用同一套整数运算复现全部五组直方图与估计量均值：
// snippet: 取自 lantern/core/src/lantern/cg_trace_05_1_estimator.ts 的掩码对照与取模直方图（逐字摘录，完整文件另含 chi2 与 estimateMeanSquare）
/** **故意漏掉掩码**的 SplitMix64：只保留 `state += GOLDEN`，三次乘法不截断。 */
export function rngNextUnmasked(state: bigint): [bigint, bigint] {
const s = state + RNG_GOLDEN;
let z = s;
z = (z ^ (z >> 30n)) * RNG_MIX1;
z = (z ^ (z >> 27n)) * RNG_MIX2;
z = z ^ (z >> 31n);
return [s, z];
}
/** 第 index 条独立子流的起点（index 从 1 开始）。 */
export function streamStart(seed0: bigint, index: number, stride: bigint): bigint {
return (seed0 + BigInt(index) * stride) & RNG_MASK;
}
/** 取低 `bits` 位再取模：值域 2^bits。返回计数、前 256 个样本、末状态。 */
export function histTruncated(
state: bigint, n: number, bits: number, m: number,
): { counts: number[]; head: number[]; state: bigint } {
const mask = (1n << BigInt(bits)) - 1n;
const nb = BigInt(n);
const counts = new Array<number>(n).fill(0);
const head: number[] = [];
let s = state;
for (let i = 0; i < m; i += 1) {
const r = rngNext(s);
s = r[0];
const v = Number((r[1] & mask) % nb);
counts[v] += 1;
if (i < 256) head.push(v);
}
return { counts, head, state: s };
}
八条断言全绿，其中三条是硬口径：前 8 个 64 位输出逐位相同（'15224344366058555842' 起）、五组直方图计数数组逐项相同（[750139, 624774, 625087] 等）、六个估计量均值逐位相同（[21315, 21781, 21703, 21791, 21827, 21845]）。一条负例把掩码摘掉，断言 z >> 48 > 65535 立刻成立，证明这个判据确实有区分力。整套 node --test 跑完 2.8 秒。到这里，“算一个积分"已被拆成三件独立可验证的事：估计量怎么写、随机数怎么保证语义、偏差从哪来——三件事都没涉及"估得准不准”，那是下一节的问题。
阅读导航：上一节：4.3 递归的光线树 · 下一节：5.2 方差与收敛 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。