随机过程模拟

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

适用场景

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

输入

必需

输入	类型	描述
process_type	string	过程类型："dtmc"、"ctmc"、"random_walk"、"brownian_motion"、"sde"、"mcmc"
parameters	dict	过程特定参数（转移矩阵、漂移/扩散系数、目标密度等）
n_paths	integer	独立样本路径数量
n_steps	integer	每条路径的时间步数（或 MCMC 总迭代数）

可选

输入	类型	默认值	描述
initial_state	scalar/vector	过程特定	每条路径的起始状态或分布
dt	float	0.01	连续时间离散化的时间步长
seed	integer	random	用于可重复性的随机种子
burn_in	integer	n_steps / 10	要丢弃的初始步数（MCMC）
thinning	integer	1	保留每第 k 个样本以减少自相关
variance_reduction	string	"none"	方法："none"、"antithetic"、"stratified"、"control_variate"
target_function	callable	none	沿路径评估的用于蒙特卡洛估计的函数

步骤

第 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. 根据过程类型选择适当的算法：

过程	方法	关键性质
DTMC	从转移行直接采样	精确
CTMC	Gillespie 算法 (SSA)	精确，事件驱动
CTMC（近似）	Tau-leaping	近似，高速率时更快
随机游走	直接采样增量	精确
布朗运动	高斯增量的累积和	固定 dt 时精确
SDE（一般）	Euler-Maruyama	强收敛 0.5 阶，弱收敛 1.0 阶
SDE（高阶）	Milstein	强收敛 1.0 阶（标量噪声）
SDE（刚性）	隐式 Euler-Maruyama	刚性漂移稳定
MCMC（一般）	Metropolis-Hastings	渐近精确
MCMC（梯度）	Hamilton 蒙特卡洛 (HMC)	高维混合更好
MCMC（条件）	Gibbs 采样器	条件分布可得时精确

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（或对事件驱动方法如 Gillespie 动态分配）。

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 个样本（主要用于 MCMC）。

预期结果： n_paths 条完整轨迹存储在内存中，带有可选的函数评估。MCMC 接受率在目标范围内。

失败处理： 如果模拟产生 NaN 或 Inf 值，对 SDE 方法减小 dt 或检查参数有效性。如果 MCMC 接受率接近 0% 或 100%，调整提议尺度。

第 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.01	...
ESS（最小）	...	> 400	...
Geweke z（最大绝对值）	...	< 2.0	...
接受率	...	0.15-0.50	...

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

失败处理： 如果 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。如果方差缩减恶化了估计（VRF < 1），禁用它——控制变量或对偶方案可能不适合该问题。

第 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 进行模拟实现后验预测检验和合成数据生成