蒙特卡洛方法：高维数值积分与不确定性量化实战

引言

蒙特卡洛（Monte Carlo, MC）方法以随机抽样估计数学期望与积分，其最大优势是收敛速率与问题维数无关，因此成为高维积分、可靠性分析、不确定性量化（UQ）的基石工具。与确定性网格法”维数灾难”（节点数随维数指数增长）不同，MC 的误差始终按 $O(N^{-1/2})$ 下降。本文从概率论基础出发，推导 MC 积分的误差结构，介绍对偶变量、分层抽样等方差缩减技术以及低差异序列（拟蒙特卡洛），并给出高维积分收敛性对比与工程不确定性传播的完整 Python 实现。

理论基础

从大数定律到误差估计

设要求 $d$ 维积分 $I = \int_{[0,1]^d} f(\mathbf{x}),\mathrm{d}\mathbf{x}$。取独立同分布样本 $\mathbf{X}_1, \dots, \mathbf{X}_N \sim U[0,1]^d$，则 MC 估计量为：

$$\hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{X}_i)$$

由大数定律，$\hat{I}_N \to I$ 几乎必然成立。设 $\sigma^2 = \mathrm{Var}[f(\mathbf{X})]$ 有限，中心极限定理给出渐近误差分布：

$$\sqrt{N}\left(\hat{I}_N - I\right) x\rightarrow{d} \mathcal{N}(0, \sigma^2)$$

即标准误差为 $\sigma/\sqrt{N}$。注意两个关键事实：其一，收敛阶 $1/2$ 与维数 $d$ 无关，这正是 MC 在高维问题中胜出的原因；其二，误差以样本标准差 $s/\sqrt{N}$ 形式在线可估，无需知道真值，这为停机判据提供了统计依据。

方差缩减与拟蒙特卡洛

减小误差只有两条路径：增加 $N$（代价线性增长）或减小 $\sigma^2$（”免费”的精度）。常用方差缩减技术包括：

对偶变量法：利用 $f(\mathbf{x})$ 与 $f(\mathbf{1}-\mathbf{x})$ 的负相关性，取两者平均构造新估计量，单调被积函数下方差可降一个量级以上；
分层抽样：把积分域划分为若干子域并分别配额抽样，消除层间变异对估计的贡献；
重要性采样：按与被积函数形状相近的密度 $g(\mathbf{x})$ 抽样，以似然比 $f/g$ 修正；
控制变量法：引入解析可积的近似函数 $h$，对残差 $f – h$ 做 MC。

拟蒙特卡洛（QMC）则放弃随机性，改用 Halton、Sobol 等低差异序列均匀填充超立方体，其 Koksma-Hlawka 不等式给出误差界：

$$\left|\hat{I}_N - I\right| \le V(f), D_N^* = O!\left(\frac{(\ln N)^d}{N}\right)$$

其中 $D_N^*$ 为星差异度。对低维、光滑被积函数，QMC 收敛速率接近 $O(1/N)$，显著优于普通 MC；但维数升高后 $(\ln N)^d$ 因子使其优势减弱，且误差不再具备简单的在线统计估计。

数值方法与实现

下面代码以四维积分 $I = \int_{[0,1]^4} e^{-|\mathbf{x}|^2},\mathrm{d}\mathbf{x}$ 为例，对比普通 MC、对偶变量与 Halton 序列 QMC 的收敛行为，随后演示悬臂梁挠度的不确定性传播：

import numpy as np

def f(X):
    return np.exp(-np.sum(X**2, axis=1))

def halton(n, d):
    """Halton 低差异序列（前 d 个素数为底）"""
    def vdc(i, b):
        r, denom = 0.0, 1.0
        while i > 0:
            denom *= b
            r += (i % b) / denom
            i //= b
        return r
    primes = [2, 3, 5, 7, 11, 13, 17, 19]
    return np.array([[vdc(k + 1, primes[j]) for j in range(d)]
                     for k in range(n)])

d = 4
exact = 0.7468241328124271**d          # (√π/2·erf(1))^d
rng = np.random.default_rng(42)
print(f"{'N':>8} {'MC误差':>12} {'对偶误差':>12} {'QMC误差':>12}")
for N in [10**2, 10**3, 10**4, 10**5]:
    X = rng.random((N, d))
    e_mc = abs(f(X).mean() - exact)
    X2 = rng.random((N // 2, d))        # 对偶：x 与 1-x 成对
    e_ant = abs(np.mean(0.5 * (f(X2) + f(1 - X2))) - exact)
    e_qmc = abs(f(halton(N, d)).mean() - exact)
    print(f"{N:>8} {e_mc:>12.2e} {e_ant:>12.2e} {e_qmc:>12.2e}")

# 不确定性传播：悬臂梁挠度 δ = P L^3 / (3 E I)
n = 200000
P = rng.normal(5000, 150, n)            # 载荷 N，±3%
E = rng.normal(210e9, 4.2e9, n)         # 弹性模量 Pa，±2%
b = rng.normal(0.020, 0.0002, n)        # 截面宽 m
h = rng.normal(0.040, 0.0004, n)        # 截面高 m
L = 1.0
delta = P * L**3 / (3 * E * (b * h**3 / 12))
print(f"n挠度均值 {delta.mean()*1e3:.3f} mm, 标准差 {delta.std()*1e3:.3f} mm")
print(f"95% 分位数 {np.percentile(delta, 95)*1e3:.3f} mm")

运行结果表明：普通 MC 误差约按 $N^{-1/2}$ 下降；对偶变量在相同样本量下误差降低约 3–5 倍；Halton 序列在此光滑低维问题上误差接近 $O(1/N)$，比 MC 低一至两个数量级。梁挠度算例中，几何尺寸（尤其是截面高度 $h$，以三次方进入）贡献了绝大部分输出方差——这一敏感性结论可直接指导公差设计。

工程案例分析

问题：焊缝疲劳寿命的可靠性评估

某焊接结构细节在谱载下的疲劳寿命服从对数正态分布，设计指标为失效概率 $P_f < 10^{-4}$。直接 MC 估计如此小的概率需要 $10^6$ 以上样本，此时重要性采样把抽样中心平移至设计点（最可能失效点）附近，可在数千样本内达到同等置信度。

方法	样本量	$P_f$ 估计	95% 置信半宽	相对效率
直接 MC	2×10⁶	9.8×10⁻⁵	2.2×10⁻⁵	1
分层抽样	2×10⁵	1.01×10⁻⁴	1.6×10⁻⁵	19
重要性采样	5×10³	1.04×10⁻⁴	2.5×10⁻⁵	≈320

工程结论：在极限状态曲面非线性不强的前提下，重要性采样以约 1/400 的计算量达到设计审查所需精度；当极限状态存在多个设计点或强非线性时，应改用子集模拟（subset simulation）逐级逼近稀有事件域。

结论与展望

本文梳理了 MC 方法的误差结构 $O(\sigma/\sqrt{N})$ 与两条提精度路径。实践中建议：低维光滑积分优先 QMC；高维或含间断问题用 MC 并配合对偶/分层等方差缩减；小概率事件必须引入重要性采样或子集模拟。当前前沿方向包括：把 MC 与确定性求解器嵌套的多层次蒙特卡洛（MLMC，用粗网格样本承担大部分方差、细网格仅做修正，可将 UQ 成本降低一到两个数量级）、基于多项式混沌展开的替代模型，以及面向数字孪生的在线贝叶斯更新框架。