3. 随机数生成与随机抽样

蒙特卡洛仿真的全部随机性都来自一个源头：随机数。但计算机本身是确定性的，它无法凭空产生真正的随机性，只能通过确定性算法生成看起来像随机的数列，这就是伪随机数。理解伪随机数的生成机制、均匀分布抽样方法，以及如何把均匀随机数变换成任意目标分布的样本，是搭建一切蒙特卡洛模型的基础。本节从底层生成器讲到四种常用抽样方法，并给出可直接运行的实现思路。

一、伪随机数生成器：确定性的“随机”

伪随机数生成器（PRNG）本质上是一个状态机：给定一个初始种子（seed），它按固定递推公式不断产生新的状态，并把状态映射为 [0,1) 区间上的数值。只要种子相同，产生的序列就完全相同——这既是缺点（可预测），也是巨大优点（可复现）。仿真调试时，固定种子能让每次运行结果一致，便于定位问题。

最经典的生成器是线性同余发生器（LCG），递推式为：

X_{n+1} = (a * X_n + c) mod m
U_n = X_n / m        // 归一化到 [0,1)

其中 m 是模数，a 是乘子，c 是增量，X_0 是种子。例如取 m = 2^31、a = 1103515245、c = 12345，就是早期 C 标准库常用的参数。LCG 实现简单、速度快，但低位随机性差，且周期有限，高维上会出现规律性“晶格结构”，因此现代仿真更推荐使用梅森旋转（Mersenne Twister）或 PCG 等生成器。Python 的 random 模块底层就是梅森旋转，周期长达 2^19937−1，在绝大多数工程场景下足够可靠。

实践建议：仿真中不要自己手写 LCG 用于正式分析，优先调用成熟库（Python 的 random / numpy.random，C++ 的 <random>）。自写生成器容易在周期、均匀性和相关性上埋下隐患。

二、均匀分布抽样：一切抽样的起点

PRNG 直接输出的就是 [0,1) 上的均匀随机数 U。由 U 可以轻松得到任意区间 [a,b) 上的均匀样本：

X = a + (b - a) * U

例如要生成 [10, 20] 上的均匀随机数，取 U = 0.37，则 X = 10 + 10×0.37 = 13.7。均匀抽样是所有其他抽样方法的“原材料”，逆变换法和接受拒绝法都建立在它之上。

三、逆变换法：用累积分布函数的反函数抽样

逆变换法（Inverse Transform Method）是理论上最优雅的抽样方法。它的依据是：若 U ~ Uniform(0,1)，F 是某目标分布的累积分布函数（CDF），则 X = F⁻¹(U) 服从该分布。直观理解是：U 落在 [0,1) 的哪个位置，就对应 CDF 曲线上的哪个分位点，从而“搬运”到目标分布上。

以指数分布为例，其 CDF 为 F(x) = 1 − e^(−λx)，反解得到：

X = -ln(1 - U) / λ

由于 1−U 与 U 同为均匀分布，工程上常简写为 X = −ln(U)/λ。假设 λ = 0.5，取 U = 0.2，则 X = −ln(0.2)/0.5 ≈ 3.22。逆变换法的优点是每个样本只需一次计算，效率高；缺点是很多分布的 F⁻¹ 没有解析式（如正态分布），此时需要数值近似或改用其他方法。

四、接受拒绝法：用“试探—筛选”逼近任意分布

当目标分布 f(x) 难以直接抽样时，接受拒绝法（Acceptance-Rejection Method）提供了一条通用路径。它需要一个容易抽样的建议分布 g(x)，以及一个常数 M，使得对定义域内所有 x 都有 f(x) ≤ M·g(x)。算法流程如下：

  1. 从建议分布 g 中抽取一个候选样本 Y；
  2. 从 Uniform(0,1) 抽取 U；
  3. 若 U ≤ f(Y) / (M·g(Y))，则接受 Y 作为样本；否则拒绝，回到第 1 步。

接受率等于 1/M，因此 M 越接近 1，效率越高。若 M 过大，大量候选样本被丢弃，计算成本急剧上升。例如用 Uniform(0,1) 作为建议分布去抽样 Beta(2,3) 分布，其密度最大值为 1.5，取 M = 1.5，则平均接受率约 66.7%，效率尚可。若换成一个形状差异很大的建议分布，M 可能达到 10 以上，接受率不足 10%，此时应更换建议分布。

接受拒绝法的关键陷阱：M 必须严格满足 f(x) ≤ M·g(x) 对所有 x 成立。若 M 取值偏小，会错误地拒绝本应接受的样本，导致抽样分布偏离目标分布，且这种偏差往往不易察觉。

五、方法对比与实现要点

四种方法各有适用场景：PRNG 是底层引擎，均匀抽样是基础原料，逆变换法适合 CDF 可解析反解的分布（指数、威布尔、逻辑斯蒂等），接受拒绝法适合密度函数已知但难以反解的分布（Beta、截断正态等）。实际工程中，正态分布通常用 Box-Muller 变换或 Ziggurat 算法，它们比朴素接受拒绝法更快。

下面给出一个 Python 实现示例，同时演示逆变换法和接受拒绝法：

import numpy as np

rng = np.random.default_rng(seed=42)   # 固定种子，保证可复现

# 逆变换法：指数分布 lambda=0.5
def exp_inverse(n, lam):
    u = rng.random(n)
    return -np.log(u) / lam

# 接受拒绝法：Beta(2,3)，建议分布 Uniform(0,1)，M=1.5
def beta_rejection(n, M=1.5):
    samples = []
    while len(samples) < n:
        y = rng.random()
        u = rng.random()
        f = 12 * y * (1 - y)**2      # Beta(2,3) 密度
        if u <= f / M:
            samples.append(y)
    return np.array(samples)

exp_samples = exp_inverse(10000, 0.5)
beta_samples = beta_rejection(10000)
print("指数均值:", exp_samples.mean())   # 理论值 1/0.5 = 2
print("Beta均值:", beta_samples.mean())   # 理论值 2/(2+3) = 0.4

运行后，指数样本均值应接近 2，Beta 样本均值应接近 0.4。若偏差明显偏大，首先检查种子是否固定、样本量是否足够，再检查接受拒绝法的 M 是否满足上界条件。样本量较小时（如 n=100），均值波动属于正常统计涨落，不能据此判定算法错误。

随机数生成与抽样方法关系图 伪随机数生成器 种子 → 状态递推 均匀分布抽样 U ~ Uniform(0,1) 逆变换法 X = F⁻¹(U) 适用：CDF 可解析反解 接受拒绝法 U ≤ f(Y)/(M·g(Y)) 适用：密度已知难反解 目标分布样本（指数 / Beta / 正态 …） 关键参数 种子 seed / 接受率 1/M

总结来看，伪随机数生成器提供可复现的均匀随机源，均匀抽样是基础，逆变换法以解析反解换取高效率，接受拒绝法以通用性换取灵活性。选择哪种方法，取决于目标分布的 CDF 是否可反解、密度函数是否已知、以及对计算效率的要求。在后续的蒙特卡洛仿真中，这些抽样方法将直接决定输入随机变量的质量，进而影响整个稳健设计分析的可靠性。

← 上一章 📖 返回目录 下一章 →