随机过程模拟

模拟随机过程的样本路径——包括离散马尔可夫链、连续时间过程、随机微分方程和 MCMC 采样器——并进行收敛诊断、方差缩减技术和轨迹可视化。

适用场景

• 需要从随机过程生成用于估计、预测或可视化的样本路径
• 解析解不可得，模拟是唯一可行的方法
• 正在运行蒙特卡洛估计，需要收敛保证和不确定性量化
• 想要对比解析结果（平稳分布、命中时间）与经验模拟进行验证
• 需要使用 MCMC 从复杂后验分布中采样
• 在承诺进行完整解析处理之前对随机模型进行原型验证

输入

必需

process_type

"dtmc"

"ctmc"

"random_walk"

"brownian_motion"

"sde"

"mcmc"

parameters

n_paths

n_steps

可选

initial_state

dt

seed

burn_in

n_steps / 10

thinning

variance_reduction

"none"

"none"

"antithetic"

"stratified"

"control_variate"

target_function

步骤

第 1 步：定义过程模型和参数

1.1. 识别过程类型并收集所有必需参数：

• DTMC：转移矩阵

  P

  和状态空间。验证

  P

  是行随机的。
• CTMC：速率矩阵

  Q

  。验证行和为 0 且非对角线元素非负。
• 随机游走：步长分布（例如

  {-1, +1}

  等概率），边界（如有）。
• 布朗运动：漂移

  mu

  ，波动率

  sigma

  ，维度

  d

  。
• SDE（伊藤）：漂移函数

  a(x,t)

  ，扩散函数

  b(x,t)

  。
• MCMC：目标对数密度，提议机制（随机游走 Metropolis、Hamilton、Gibbs 分量）。

1.2. 验证参数一致性：

• 矩阵维度与状态空间大小匹配。
• SDE 系数满足增长和 Lipschitz 条件（至少非正式地）以适应所选求解器。
• MCMC 提议在目标分布的支撑上有良好定义。

1.3. 设置随机种子以确保可重复性。

预期结果： 一个完全指定的随机模型，参数已验证，随机状态可重复。

失败处理： 如果参数不一致（例如非随机矩阵），在继续之前修正。如果 SDE 系数具有病理性，考虑使用不同的离散化方案。

第 2 步：选择模拟方法

2.1. 根据过程类型选择适当的算法：

dt

2.2. 对于 SDE 方法，选择足够小的

dt

dt = 0.01

2.3. 对于 MCMC，调整提议尺度以达到大约以下接受率：

• 高维随机游走 Metropolis：23.4%
• 一维目标：57.4%
• HMC：65-90%（取决于轨迹长度）

2.4. 如果请求方差缩减，进行配置：

• 对偶变量：对于每条具有随机增量

  Z

  的路径，同时模拟

  -Z

  。
• 分层采样：将概率空间分区并在每层中采样。
• 控制变量：识别一个具有已知期望的相关量以减少方差。

预期结果： 选定的模拟算法匹配过程类型，具有适当的调参。

失败处理： 如果所选方法不稳定（例如 Euler-Maruyama 发散），切换到隐式方法或减小

dt

第 3 步：实现并运行模拟

3.1. 为

n_paths

n_steps

3.2. 对每条路径

i = 1, ..., n_paths

DTMC / 随机游走：

• 设置

  x[0] = initial_state

• 对

  t = 1, ..., n_steps

  ：根据

  x[t-1]

  从转移分布中采样

  x[t]

CTMC (Gillespie)：

• 设置

  x[0] = initial_state

  ，

  time = 0

• 当

  time < T_max

  时：
  • 计算总速率

    lambda = -Q[x, x]

  • 采样停留时间

    tau ~ Exp(lambda)

  • 从转移概率

    Q[x, j] / lambda

    （

    j != x

    ）采样下一状态
  • 更新

    time += tau

    ，记录转移

SDE (Euler-Maruyama)：

• 设置

  x[0] = initial_state

• 对

  t = 1, ..., n_steps

  ：

  • dW = sqrt(dt) * N(0, I)

    （Wiener 增量）

  • x[t] = x[t-1] + a(x[t-1], t*dt) * dt + b(x[t-1], t*dt) * dW

MCMC (Metropolis-Hastings)：

• 设置

  x[0] = initial_state

• 对

  t = 1, ..., n_steps

  ：
  • 提议

    x' ~ q(x' | x[t-1])

  • 计算接受率

    alpha = min(1, p(x') * q(x[t-1]|x') / (p(x[t-1]) * q(x'|x[t-1])))

  • 以概率

    alpha

    接受：如果接受则

    x[t] = x'

    ，否则

    x[t] = x[t-1]

  • 记录接受决策

3.3. 如果提供了

target_function

3.4. 应用稀疏化：保留每第

thinning

3.5. 丢弃每条路径开头的

burn_in

预期结果：

n_paths

失败处理： 如果模拟产生 NaN 或 Inf 值，对 SDE 方法减小

dt

第 4 步：应用收敛诊断

4.1. 轨迹图：为路径子集绘制每个分量随时间的值。对平稳性进行目视检查（无趋势、方差稳定）。

4.2. Gelman-Rubin 诊断 (R-hat)：对于多链 MCMC：

• 计算链内方差

  W

  和链间方差

  B

  。

• R_hat = sqrt((n-1)/n + B/(n*W))

• R_hat < 1.01

  （严格）或

  R_hat < 1.1

  （宽松）表示收敛。

4.3. 有效样本量 (ESS)：

• 估计递增滞后处的自相关。

• ESS = n_samples / (1 + 2 * sum(autocorrelations))

• 经验法则：

  ESS > 400

  用于可靠的后验总结。

4.4. Geweke 诊断：比较每条链前 10% 和后 50% 的均值。z 分数应在 [-2, 2] 范围内表示收敛。

4.5. 非 MCMC 过程：验证时间平均统计量（均值、方差）随路径长度增加而稳定。绘制运行平均值。

4.6. 报告汇总表：

预期结果： 所有收敛诊断通过其阈值。轨迹图显示稳定、良好混合的链。

失败处理： 如果 R-hat > 1.1，运行更长的链或改进提议。如果 ESS 非常低，增加稀疏化或切换到更好的采样器（例如 HMC）。如果 Geweke 失败，延长预热期。

第 5 步：计算带置信区间的汇总统计量

5.1. 对每个感兴趣的量（状态占用、函数期望、命中时间）：

• 计算点估计为路径间的样本均值（预热和稀疏化之后）。
• 使用有效样本量计算标准误：

  SE = SD / sqrt(ESS)

  。

5.2. 构造置信区间：

• 正态近似：

  estimate +/- z_{alpha/2} * SE

• 对偏斜分布，使用百分位 bootstrap 或批次均值。

5.3. 如果应用了方差缩减，计算方差缩减因子：

• VRF = Var(朴素估计量) / Var(缩减估计量)

• 报告有效加速比。

5.4. 对蒙特卡洛积分估计：

• 报告估计值、标准误、95% CI、ESS 和函数评估次数。

5.5. 对分布估计：

• 计算经验分位数（中位数、2.5 百分位、97.5 百分位）。
• 连续量的核密度估计。

5.6. 将所有汇总统计量及其不确定性制表。

预期结果： 点估计及其标准误和置信区间。方差缩减（如已应用）产生 VRF > 1。

失败处理： 如果置信区间太宽，增加

n_paths

n_steps

第 6 步：可视化轨迹和分布

6.1. 轨迹图：绘制样本路径的代表性子集（5-20 条路径）随时间的变化。使用透明度处理重叠路径。

6.2. 集合统计：叠加所有路径的平均轨迹和逐点 95% 置信带。

6.3. 边际分布：在选定的时间点，绘制跨路径的状态分布直方图或密度估计。

6.4. 平稳分布比较：如果有解析平稳分布可用，将其叠加在最终时间切片的经验直方图上。

6.5. 自相关图：对 MCMC，绘制每个分量到合理滞后的自相关函数 (ACF)。

6.6. 诊断仪表板：将轨迹图、ACF 图、运行均值图和边际密度组合成单一多面板图形，进行综合评估。

6.7. 将所有图形保存为矢量格式（PDF/SVG）和光栅格式（PNG），用于文档记录。

预期结果： 出版质量的图形，展示轨迹行为、分布收敛和诊断摘要。解析解（如可用）与经验结果匹配。

失败处理： 如果可视化揭示了模型未预期的非平稳性或多模态，回到第 1-2 步检查参数或方法错误。如果图形杂乱，减少显示的路径数量或增加图形尺寸。

验证清单

• 所有模拟轨迹保持在有效状态空间内（无越界值、无 NaN/Inf）
• 对 DTMC/CTMC：经验平稳分布收敛到解析分布（在预期蒙特卡洛误差范围内）
• 对 SDE：将

  dt

  减半不会定性地改变结果（收敛阶检查）
• 对 MCMC：R-hat < 1.01，ESS > 400，Geweke z 分数在 [-2, 2] 范围内
• 置信区间宽度按

  1/sqrt(n_paths)

  比例减小（中心极限定理）
• 方差缩减技术产生 VRF > 1（估计改善而非恶化）
• 可重复性：使用相同种子重新运行产生相同结果

常见问题

• MCMC 预热不足：从不良初始状态开始需要长时间预热，样本才能代表目标分布。始终检查轨迹图并使用收敛诊断，而不是猜测预热长度
• 刚性 SDE 的 Euler-Maruyama 不稳定性：如果漂移项有大梯度，显式 Euler-Maruyama 可能发散。切换到隐式方法或使用自适应步长
• 混淆 SDE 的强收敛和弱收敛：强收敛衡量逐路径误差（对单个轨迹重要）；弱收敛衡量分布误差（对期望足够）。Euler-Maruyama 弱阶 1.0 但强阶 0.5
• 伪随机数生成器质量：对于非常长的模拟，低质量 RNG 可能产生相关样本。使用经过良好测试的生成器（Mersenne Twister、PCG 或 Xoshiro）并验证独立性
• 忽略 MCMC 中的自相关：将自相关的 MCMC 样本视为独立会低估不确定性。始终使用有效样本量而非原始样本数来计算标准误
• 非单调函数的对偶变量：对偶采样仅在估计量是底层均匀分布的单调函数时减少方差。对非单调函数，它可能增加方差
• 大规模模拟的内存：存储许多长路径的所有时间步可能耗尽内存。当不需要完整轨迹用于可视化时，使用在线统计（运行均值、方差）

相关技能

• 马尔可夫链建模 — 提供模拟验证的转移矩阵和解析解
• 拟合隐马尔可夫模型 — 从拟合的 HMM 进行模拟实现后验预测检验和合成数据生成