4. 蒙特卡洛方法基础
蒙特卡洛方法是不确定度量化中最通用、最容易落地的一类数值手段。它的核心思想非常朴素：把待求量写成某个随机变量的期望，然后通过大量独立抽样，用样本均值去逼近这个期望。设待求量为
θ = E[g(X)] = ∫ g(x) p(x) dx
其中 X 是输入随机向量，p(x) 是其联合概率密度。若能从 p(x) 中抽取独立同分布样本 X₁, X₂, …, XN，则估计量为
θ̂N = (1/N) Σi=1N g(Xi)
这一形式对 g 几乎没有光滑性、线性或低维要求，这正是蒙特卡洛方法在复杂仿真模型中不可替代的原因。
随机抽样：从分布到样本
抽样的第一步是把均匀随机数转换为目标分布的样本。最基本的方法是逆变换法：若 U ~ Uniform(0,1)，F 为目标分布函数且存在反函数，则 X = F⁻¹(U) 服从目标分布。例如指数分布 F(x)=1−e−λx，可得 X = −ln(1−U)/λ。对正态分布等无解析反函数的分布，常用 Box–Muller 变换：由两个独立均匀数 U₁、U₂ 生成 Z = √(−2lnU₁)·cos(2πU₂)，得到标准正态样本。
对多维相关输入，通常先用 Cholesky 分解相关矩阵 Σ = LLᵀ，再令 X = μ + LZ，其中 Z 为标准正态独立样本，从而得到均值 μ、协方差 Σ 的相关正态样本。若输入分布为任意边缘分布，可先按目标相关矩阵生成高斯 copula 样本，再逐维做逆变换。抽样质量的关键在于：随机数发生器周期足够长、样本之间无可见相关性、低差异序列（如 Sobol、Halton）在需要加速收敛时使用。
大数定律与收敛性
大数定律保证 θ̂N 依概率收敛到 θ：当 N → ∞ 时，对任意 ε > 0，有 P(|θ̂N − θ| > ε) → 0。若 g(X) 的方差 σ² = Var[g(X)] 有限，中心极限定理进一步给出
√N (θ̂N − θ) / σ → N(0,1)
这意味着蒙特卡洛误差以 O(N−1/2) 的速度下降，与输入维数无关——这是它相对确定性网格法在高维问题中的根本优势。但收敛速度也意味着：要把误差缩小 10 倍，样本量需增加约 100 倍。因此，样本量确定和方差缩减技术是实际应用中的关键环节。
误差估计
蒙特卡洛估计的误差由样本标准差决定。记样本方差 s² = (1/(N−1)) Σ (g(Xi) − θ̂N)²，则标准误为 SE = s/√N。在 95% 置信水平下，误差半宽约为 1.96·SE。需要注意，这里的误差是统计误差，不包含模型结构误差、输入分布设定误差和数值求解误差。实际报告中应把统计误差与其它不确定度来源分开列出，避免把“样本量不够”误判为“模型不准”。
当 g(X) 的方差很大或存在重尾时，样本均值的收敛可能很慢，甚至出现“看起来收敛但实际未收敛”的假象。此时应检查样本直方图、运行均值轨迹和自相关函数，必要时采用分层抽样、重要性抽样、对偶变量或控制变量等方差缩减手段。
样本量确定
样本量确定通常分两步：先做小规模预采样（如 N₀ = 200～500）估计 s，再按目标误差半宽 δ 计算所需样本量 N ≥ (z1−α/2 · s / δ)²。例如预采样得 s = 0.5，要求 95% 置信下误差不超过 0.02，则 N ≥ (1.96×0.5/0.02)² ≈ 2401，取 2500。若模型单次评估耗时较长，可先设定可接受的计算预算，再反推能达到的误差水平，并在报告中如实说明。
基础实现
下面给出一个可直接运行的最小实现，用于估计标准正态输入下 g(X)=X² 的期望（理论值为 1），并输出估计值、标准误和 95% 置信区间。
import numpy as np
rng = np.random.default_rng(20240501)
N = 10000
x = rng.standard_normal(N) # 标准正态抽样
g = x**2 # 模型评估
theta_hat = g.mean() # 样本均值
s = g.std(ddof=1) # 样本标准差
se = s / np.sqrt(N) # 标准误
ci = (theta_hat - 1.96*se, theta_hat + 1.96*se)
print(f"theta_hat = {theta_hat:.4f}")
print(f"SE = {se:.4f}")
print(f"95% CI = [{ci[0]:.4f}, {ci[1]:.4f}]")
一次典型运行可能输出 θ̂ ≈ 0.995、SE ≈ 0.014、95% CI ≈ [0.968, 1.022]，区间覆盖理论值 1，说明估计合理。若把 N 提高到 100000，SE 应约缩小到原来的 1/√10 ≈ 0.316 倍，即约 0.0045。实际使用时建议固定随机种子以便复现，同时用不同种子重复若干次，观察估计值的波动是否与 SE 一致。