3. 蒙特卡洛积分原理与误差分析

高维积分之所以困难，根源在于确定性网格方法所需的采样点数随维数指数增长。蒙特卡洛方法换了一条思路：不再试图在规则网格上“铺满”积分区域，而是用随机采样点的函数值平均来逼近积分值。这一节我们从无偏性、方差与收敛速率三个层面推导它的数学基础，再讨论随机数质量对精度的影响，最后用高维球体积计算验证误差阶。

估计量的构造与无偏性

设积分目标为 I = ∫_Ω f(x) dx，其中 Ω 是 d 维区域。把积分改写为期望形式：若在 Ω 上按均匀分布抽取随机点 X，则

I = V · E[f(X)]，  V = ∫_Ω 1 dx 为区域体积

于是自然构造估计量：独立抽取 N 个均匀样本 X₁,…,X_N，令

Î_N = (V / N) · Σ_{i=1}^{N} f(X_i)

无偏性来自期望的线性性：E[Î_N] = (V/N)·Σ E[f(X_i)] = V·E[f(X)] = I。注意这里没有对 f 做任何光滑性假设，只要求 f 可积且方差有限。这正是蒙特卡洛相对确定性求积的核心优势——收敛性不依赖被积函数的正则性，只依赖其统计特性。

方差与收敛速率

由于样本独立，估计量的方差可直接相加：

Var(Î_N) = (V²/N²) · Σ Var(f(X_i)) = V²·σ²/N
σ² = Var(f(X)) = ∫_Ω f²(x)dx / V − (I/V)²

标准差为 σ_Î = V·σ/√N，即误差量级为 O(N^{-1/2})。这个速率与维数 d 无关，是高维问题的关键优势；代价是它比一维梯形公式的 O(N^{-2}) 慢得多。因此蒙特卡洛并非在所有场合都优于确定性方法，而是在维数足够高、被积函数不够光滑时胜出。粗略估算：若一维网格法误差为 h² ~ N^{-2/d}，当 2/d < 1/2 即 d > 4 时，蒙特卡洛的收敛阶开始占优。

实际中 σ 未知，用样本标准差估计：

s² = (1/(N−1)) · Σ (f(X_i) − f̄)²，  f̄ = (1/N)Σ f(X_i)
误差标准误 ≈ V·s/√N

这给出一个可操作的置信区间：Î_N ± 1.96·V·s/√N 约对应 95% 置信水平。工程上常据此判断“再增加多少样本才能把误差压到目标值”——样本量需按 N ∝ 1/ε² 增长，误差减半意味着样本量翻四倍。

随机数质量的影响

上述推导全部建立在“样本独立同分布且服从均匀分布”的前提上。伪随机数发生器（PRNG）若周期过短、低位相关性明显或种子选择不当，会破坏独立性，使实际误差偏离 N^{-1/2} 的理论阶。常见问题包括：

• 周期耗尽：样本数接近发生器周期时出现重复序列，方差估计失真；
• 格点结构：线性同余发生器在高维下样本点落在少数超平面上，导致某些方向的积分被系统性低估或高估；
• 低差异序列误用：拟蒙特卡洛（如 Sobol 序列）可把误差阶提升到接近 N^{-1}，但其误差界含与维数相关的常数，且需要正确的方向数配置，否则效果可能反而不如纯随机。

实践中应优先使用经过统计检验的现代发生器（如 PCG、xoshiro 系列），并在高维场景下检查样本的投影分布是否均匀。一个简单有效的自检手段是：用同一被积函数、不同种子重复计算 M 次，观察 M 个估计值的样本标准差是否与 V·s/√N 吻合；若明显偏大，说明随机源或实现存在问题。

案例：高维球体积计算与误差阶验证

考虑 d 维单位球 B_d = {x : ‖x‖₂ ≤ 1}，其体积有解析式 V_d = π^{d/2} / Γ(d/2 + 1)。在包围盒 [−1,1]^d 内均匀采样，落入球内的点比例乘以盒体积 2^d 即为估计值。取 d = 10，解析体积约为 2.5502。下面是一段可直接运行的 Python 实现：

import numpy as np

def mc_ball_volume(d, N, rng):
    pts = rng.uniform(-1.0, 1.0, size=(N, d))
    inside = np.sum(pts**2, axis=1) <= 1.0
    return (2.0**d) * inside.mean(), inside.std() / np.sqrt(N) * (2.0**d)

rng = np.random.default_rng(20240501)
for N in [10**3, 10**4, 10**5, 10**6]:
    est, se = mc_ball_volume(10, N, rng)
    print(f"N={N:>7d}  est={est:.4f}  se={se:.4f}")

示例运行结果（属于教学演示数据，非严格实测）大致呈现如下规律：N 从 10³ 增至 10⁶，标准误按 √N 比例缩小约 31.6 倍，估计值围绕 2.5502 波动的幅度同步收窄。若把每个 N 下的重复实验做 50 次，统计估计值的样本标准差，会发现它与公式给出的 V·s/√N 在同一量级，从而验证了 O(N^{-1/2}) 的误差阶。

从图中可以直观看到：理论误差曲线与示例估计值的波动点基本贴合，说明实现正确时误差阶成立；而右下角的警示框提醒我们，一旦随机源质量不达标，实际误差会系统性偏离这条曲线。因此在高维积分工程中，验证随机数质量与验证积分结果本身同等重要。