4、随机变量抽样：逆变换法、接受-拒绝法、Box-Muller变换、常用分布抽样实现
蒙特卡洛模拟的核心在于用随机数驱动模型演化，而随机数的质量与抽样效率直接决定模拟结果的精度与计算成本。在工程实践中，我们通常先获得均匀分布 U(0,1) 的随机数（由线性同余发生器或梅森旋转算法产生），再通过特定变换将其映射为目标分布。本章围绕四种经典抽样方法展开：逆变换法（Inverse Transform）、接受-拒绝法（Acceptance-Rejection）、Box-Muller 变换以及常用分布的直接抽样实现。理解这些方法的适用边界与数值稳定性，是构建高效蒙特卡洛引擎的基础能力。
4.1 逆变换法：从 CDF 反函数出发
核心原理
逆变换法基于概率积分变换定理：若随机变量 X 的累积分布函数为 F(x)，则 Y = F(X) 服从 U(0,1) 均匀分布。反之，若 U ~ U(0,1)，则 X = F⁻¹(U) 服从分布 F。该方法要求 F(x) 严格单调且可解析求逆，或至少能通过数值方法高效求解。
操作流程
- 推导 CDF：明确目标分布的累积分布函数 F(x) 及其定义域。
- 求逆函数：解析求解 x = F⁻¹(u)，若无法解析则准备数值求根方案（如二分法、牛顿法）。
- 生成均匀随机数：调用随机数发生器产生 u ~ U(0,1)。
- 映射：计算 x = F⁻¹(u)，即为目标分布的一个样本。
案例拆解：指数分布抽样
指数分布的概率密度函数为 f(x) = λe^(-λx)，x ≥ 0，其 CDF 为 F(x) = 1 - e^(-λx)。令 u = 1 - e^(-λx)，解得 x = -ln(1-u)/λ。由于 1-u 与 u 同分布，可简化为 x = -ln(u)/λ。这一简洁形式使指数分布成为逆变换法最经典的教科书案例。
常见误区
- 忽略定义域：逆变换法要求 F 在定义域上严格单调，对于分段函数需逐段处理。
- 数值求逆的精度：当 F 接近 0 或 1 时，数值求逆的相对误差可能放大，导致尾部样本失真。
- 误用离散分布：离散分布的 CDF 是阶梯函数，逆变换需通过查找表实现，而非简单解析求逆。
4.2 接受-拒绝法：用包络分布逼近目标
核心原理
当目标分布 f(x) 的 CDF 难以求逆时，接受-拒绝法提供了一种通用替代方案。其思想是：选择一个容易抽样的建议分布 g(x)，并找到常数 c ≥ 1 使得 c·g(x) ≥ f(x) 对所有 x 成立。然后从 g(x) 中抽取候选样本 Y，再从 U(0,1) 中抽取 u，若 u ≤ f(Y)/(c·g(Y)) 则接受 Y，否则拒绝并重试。
判断方法与效率分析
| 关键指标 | 公式/含义 | 工程建议 |
|---|---|---|
| 接受概率 | P(accept) = 1/c | c 越接近 1，效率越高；c 过大则大量样本被浪费 |
| 包络条件 | c·g(x) ≥ f(x), ∀x | 需在目标分布支撑集上逐点验证 |
| 建议分布选择 | g(x) 应形状接近 f(x) | 常用均匀分布、正态分布或 t 分布作为包络 |
| 尾部行为 | c·g(x) 与 f(x) 的尾部比值 | 尾部若包络不足，会产生系统性偏差 |
案例拆解：从半正态分布抽样
假设目标分布为截断正态分布 f(x) ∝ exp(-x²/2)，x ∈ [0, 3]。选择建议分布 g(x) = U(0,3) 的均匀分布，则需找到 c 使得 c/3 ≥ f(x) 在 [0,3] 上成立。f(x) 在 x=0 处取最大值 1/√(2π) ≈ 0.399，因此 c = 3 × 0.399 ≈ 1.197。接受概率约为 1/c ≈ 0.835，效率尚可。若将区间扩展到 [0, 10]，则 c 需增大至 3.99，接受概率降至 0.25，此时应改用指数分布作为建议分布以提升效率。
常见误区
- c 值选取过小：若 c·g(x) 在某些区域小于 f(x)，则会产生偏差样本，且该偏差难以通过增大样本量消除。
- 忽略建议分布的支撑集：g(x) 的支撑集必须完全覆盖 f(x) 的支撑集，否则目标分布尾部将无法被采样。
- 效率评估缺失：未计算接受概率就投入大规模模拟，可能导致计算资源浪费。
4.3 Box-Muller 变换：高效生成正态分布
核心原理
Box-Muller 变换利用二维均匀分布到二维正态分布的映射关系。设 U₁, U₂ 为两个独立的 U(0,1) 随机数，则：
Z₁ = √(-2·ln U₁) · cos(2π·U₂)
Z₂ = √(-2·ln U₁) · sin(2π·U₂)得到的 Z₁, Z₂ 为两个独立的标准正态分布 N(0,1) 样本。该方法的数学基础是极坐标变换：将均匀分布的点映射到高斯分布的径向与角度分量。
操作流程与数值稳定性
- 生成两个独立均匀随机数 u₁, u₂。
- 计算 r = √(-2·ln u₁)，θ = 2π·u₂。
- 输出 z₁ = r·cos θ，z₂ = r·sin θ。
实际工程中，更推荐使用极坐标形式的改进版（Marsaglia polar method），它避免了三角函数调用，计算速度更快且数值更稳定：
do {
x = 2·u₁ - 1; y = 2·u₂ - 1;
s = x² + y²;
} while (s ≥ 1 || s == 0);
z₁ = x · √(-2·ln s / s);
z₂ = y · √(-2·ln s / s);案例拆解：金融资产价格路径模拟
在期权定价中，几何布朗运动模型 dS = μS dt + σS dW 的离散化需要生成标准正态随机数。假设模拟 10,000 条路径、每条路径 252 个时间步，共需 2,520,000 个正态样本。使用 Box-Muller 变换，每次调用生成 2 个样本，仅需 1,260,000 次 RNG 调用。若使用逆变换法（正态 CDF 无解析逆），需调用数值求逆函数，单次耗时约为 Box-Muller 的 3-5 倍。实测中，Box-Muller 配合 Marsaglia 极坐标法可将路径模拟总耗时降低约 40%。
常见误区
- u₁ = 0 导致无穷大：部分 RNG 可能产生 0 值，需做保护处理（如取 max(u₁, 1e-12)）。
- 忽略相关性检验：若 RNG 的周期不够长，Box-Muller 生成的两个样本可能呈现微弱负相关，影响模拟精度。
- 误用单精度：在计算 √(-2·ln u₁) 时，单精度浮点数在 u₁ 接近 1 时精度不足，建议使用 double。
4.4 常用分布的直接抽样实现
核心原理与实现策略
除上述通用方法外，许多常用分布拥有专用抽样算法，这些算法在效率、数值稳定性和实现复杂度之间取得了平衡。下表总结了典型分布的推荐抽样方法：
| 分布 | 推荐方法 | 关键公式/步骤 | 注意事项 |
|---|---|---|---|
| 指数分布 Exp(λ) | 逆变换法 | x = -ln(u)/λ | λ 较大时注意浮点精度 |
| 正态分布 N(μ, σ²) | Box-Muller / Marsaglia | z = √(-2·ln u₁)·cos(2πu₂) | 避免 u₁=0；建议使用极坐标法 |
| 泊松分布 Pois(λ) | Knuth 算法（λ 较小） | 累乘均匀数直到乘积 < e^(-λ) | λ 较大时改用正态近似或 PTRD 算法 |
| 二项分布 Bin(n, p) | BTPE 算法 | 基于接受-拒绝的专用算法 | n 很大时可用正态近似 |
| Gamma 分布 | Marsaglia-Tsang 方法 | 基于接受-拒绝的改进 | 形状参数 k 较小时需特殊处理 |
| 贝塔分布 Beta(α, β) | Johnk 方法 / 比值法 | 利用 Gamma 样本比值 | α, β 均小于 1 时需谨慎 |
案例拆解：泊松分布抽样的工程实现
在排队论模拟中，顾客到达过程通常建模为泊松过程。假设平均到达率 λ = 5 人/分钟，需要生成 10,000 分钟的到达人数。Knuth 算法实现如下：
function poissonSample(lambda) {
double L = Math.exp(-lambda);
double p = 1.0;
int k = 0;
do {
k++;
p *= Math.random();
} while (p > L);
return k - 1;
}当 λ = 5 时，平均每次调用需生成约 6 个均匀随机数，效率尚可。但当 λ = 1000 时，Knuth 算法平均需生成 1001 个随机数，效率急剧下降。此时应改用 PTRD（Poisson-Triangle-Review-Draw）算法或正态近似（λ > 1000 时近似精度已足够）。
常见误区与总结复盘
- 一刀切使用逆变换法：正态分布、Gamma 分布等 CDF 无解析逆，强行数值求逆会损失效率。
- 忽略参数边界：如泊松分布 λ 极大时，Knuth 算法退化；Beta 分布参数小于 1 时，Johnk 方法可能失效。
- 未验证样本质量：任何抽样算法实现后，都应进行分布拟合检验（如 KS 检验、卡方检验）和相关性检验。
- 随机数质量被忽视：底层 RNG 的周期、均匀性和独立性直接影响所有抽样方法的效果。建议使用梅森旋转（MT19937）或 PCG 等高质量 RNG。