PyMC 时间序列分布实战指南AR、GARCH11、随机游走与 Euler-Maruyama 的完整用法【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc导读时间序列建模是贝叶斯推断中最常见的应用场景之一。PyMC 在pymc.distributions.timeseries模块中内置了一组开箱即用的时间序列分布——AR自回归过程、GARCH11波动率模型、GaussianRandomWalk、MvGaussianRandomWalk、MvStudentTRandomWalk以及通用的RandomWalk基类还提供了基于 Euler-Maruyama 离散化的随机微分方程SDE分布EulerMaruyama。阅读本文后你将掌握这些分布的数学定义、完整参数语义、在pm.Model中的调用方式以及它们底层基于 PyTensorscan与符号随机变量的实现原理能够直接将其用于金融波动率、经济指标、生物过程等时间序列的贝叶斯建模。本文以 docs/source/api/distributions/timeseries.rst 中登记的全部 API 为主线结合 pymc/distributions/timeseries.py 的源码实现与 tests/distributions/test_timeseries.py 的测试用例展开讲解。一、API 总览时间序列分布家族timeseries.rst通过autosummary指令按模块pymc导出了 6 个核心类全部定义于 pymc/distributions/timeseries.py分布类数学模型典型场景AR带 p 阶滞后的自回归过程平稳经济序列、信号预测GARCH11条件异方差波动率模型金融收益率波动率建模GaussianRandomWalk正态增量随机游走状态空间模型、时变参数MvGaussianRandomWalk多元正态增量随机游走多维状态协同演化MvStudentTRandomWalk多元 Student-t 增量随机游走重尾多维随机游走EulerMaruyamaSDE 的欧拉-丸山离散化扩散过程、随机动力系统除此之外模块源码还额外导出了一个通用构建块RandomWalk见__all__列表pymc/distributions/timeseries.py。文档中列出的三个 RandomWalk 系列分布正是通过它派生的预定义实现理解了RandomWalk其余三个的用法便可一通百通。二、RandomWalk 通用基类从初始值 创新项构造任意随机游走2.1 核心思想RandomWalk的建模思路是把随机游走拆成两个部分pymc/distributions/timeseries.pyinit_dist初始值的分布必须是未命名分布即通过.dist()API 创建的分布innovation_dist每一步创新项increment的分布同样必须是.dist()API 创建的未命名分布steps随机游走的步数steps 0仅在未提供 shape 时必需。生成过程相当于x_t x_{t-1} ε_t其中ε_t ~ innovation_dist。init_dist与innovation_dist既可以是单变量分布也可以是多元分布例如Dirichlet、MvNormal两者支持维度必须一致。需要注意的两点约束由源码RandomWalk.dist强制校验pymc/distributions/timeseries.py传入的init_dist/innovation_dist必须是RandomVariable或SymbolicRandomVariable类型的分布变量否则抛出TypeError(init_dist must be a distribution variable)两个分布必须完全独立若其中一个出现在另一个的计算图中ancestors关系会抛出ValueError且它们会被克隆clone与传入的原始变量解耦。2.2 最小示例import pymc as pm with pm.Model(): init_dist pm.Normal.dist(0, 10) innovation_dist pm.Normal.dist(0, 1) rw pm.RandomWalk( rw, init_distinit_dist, innovation_distinnovation_dist, steps100 )2.3 steps 的推断规则steps并非必须显式给出。源码RandomWalk.get_stepspymc/distributions/timeseries.py会依次从steps、shape、dims、observed中推断支持维度长度dist()中若steps与shape都未给出则抛出ValueError(Must specify steps or shape parameter)。测试test_infer_steps验证了从 shape / dims / observed 三种来源推断 steps 的一致性test_inconsistent_steps_and_shape则验证了 steps 与 shape 冲突时会抛出断言错误tests/distributions/test_timeseries.py。2.4 底层实现RandomWalkRVRandomWalk对应一个符号随机变量RandomWalkRVpymc/distributions/timeseries.py。其rv_op的实现值得注意先将init_dist缩放到形状(1, B, S)、将innovation_dist缩放到(steps, B, S)即把步数维放到最前通过moveaxis将步数维换到支持维位置拼接后沿时间轴做cumsum累计求和得到(B, T, S)形状的随机游走样本由于 logp 只能对直接叠加在 RandomVariable 之上的 dimshuffle 正确求导实现特意先换轴、再拼接并在注册的_logprob中沿-1轴对 logp 求和pymc/distributions/timeseries.py从而把时间维度折叠掉。此外RandomWalkRV还注册了支持点support point计算random_walk_support_point用 init_dist 与 innovation_dist 支持点的累计和和尺寸变换change_random_walk_size这些保证了它能在pm.Model中与 NUTS 等采样器、change_dist_sizeAPI 正常工作。三、三种预定义随机游走Gaussian / MvGaussian / MvStudentT三个预定义分布通过抽象元类PredefinedRandomWalk复用RandomWalk的构造逻辑pymc/distributions/timeseries.py每个类只需实现get_dists()返回(init_dist, innovation_dist, kwargs)即可把预置的创新项分布注入随机游走。3.1 GaussianRandomWalk数学形式x_t x_{t-1} ε_tε_t ~ Normal(mu, sigma)。参数pymc/distributions/timeseries.pymu创新项漂移tensor_like of float默认0sigma创新项标准差sigma 0默认1init_dist初始值的未命名单变量分布可选steps步数可选可用 shape 替代。重要默认行为若未显式传入init_dist源码会发出UserWarning并默认使用Normal.dist(0, 100)宽先验提示如下Initial distribution not specified, defaulting to Normal.dist(0, 100). You can specify an init_dist manually to suppress this warning.示例import pymc as pm with pm.Model(): mu pm.Normal(mu, 0, 1) sigma pm.HalfNormal(sigma, 1) y pm.GaussianRandomWalk( y, mumu, sigmasigma, steps100, init_distpm.Normal.dist(0, 10), # 显式指定初始值分布避免默认宽先验告警 )测试test_gaussian_inference给出了完整的端到端验证用已知mu2, sigma1生成 1000 步观测在Uniform先验下用GaussianRandomWalk建模并sample(chains1)后验均值能以 0.2 的容差恢复真实参数tests/distributions/test_timeseries.py。3.2 MvGaussianRandomWalk数学形式x_t x_{t-1} ε_tε_t ~ MvNormal(mu, Σ)。参数pymc/distributions/timeseries.pymu多元创新漂移tensor_like of floatcov协方差矩阵正定tau精度矩阵协方差逆矩阵正定chol协方差矩阵的 Cholesky 分解lowerchol是否为下三角矩阵默认Trueinit_dist初始值的未命名多元分布可选默认MvNormal.dist(0, I*100)并告警steps步数可选。cov、tau、chol三种参数化只需提供其一源码注释明确说明 Only one of cov, tau or chol is required内部统一交给MvNormal.dist处理。测试test_mvgaussian对三种参数化分别验证了其底层创新分布确为MvNormaltests/distributions/test_timeseries.py。值得注意的进阶用法chol可以是模型内的随机变量。测试test_mvgaussian_with_chol_cov_rv展示了与LKJCholeskyCov的组合——先构造相关矩阵的 Cholesky 因子再让MvGaussianRandomWalk的协方差来自该随机变量从而对多维随机游走的相关结构做贝叶斯推断import pymc as pm with pm.Model() as model: mu pm.Normal(mu, 0.0, 1.0, shape3) sd_dist pm.Exponential.dist(1.0, shape3) chol, corr, stds pm.LKJCholeskyCov( chol_cov, n3, eta2, sd_distsd_dist, compute_corrTrue ) mv pm.MvGaussianRandomWalk( mv, mu, cholchol, shape(10, 7, 3) # (steps, batch, dim) )该测试还断言draw(mv, draws5)的形状为(5, 10, 7, 3)直观说明了抽样数、步数、批次、维度的排布tests/distributions/test_timeseries.py。3.3 MvStudentTRandomWalk数学形式x_t x_{t-1} ε_tε_t ~ MvStudentT(nu, mu, Σ)适合重尾outlier 较多的多元增量。参数pymc/distributions/timeseries.pynu自由度intmu多元创新漂移scale尺度矩阵正定tau精度矩阵cholCholesky 分解lower是否下三角默认Trueinit_dist初始值分布可选steps步数可选。同样地scale、tau、chol三者只需提供一个cov亦可传入源码会从kwargs中弹出并转交MvStudentT.dist。测试test_mvstudentt验证了nu4时底层创新分布为MvStudentTtests/distributions/test_timeseries.py。四、AR带 p 阶滞后的自回归过程4.1 数学定义与参数AR的模型pymc/distributions/timeseries.py$$ x_t \rho_0 \rho_1 x_{t-1} \ldots \rho_p x_{t-p} \epsilon_t, \quad \epsilon_t \sim N(0,\sigma^2) $$创新项既可用标准差也可用精度参数化二者关系为 $\tau 1/\sigma^2$。参数说明rho自回归系数张量最后一个维度的第 n 个元素是第 n 阶滞后的系数sigma创新标准差sigma 0默认1与tau二选一tau创新精度tau 0可选constant布尔默认False。为True时rho的第一个元素被用作常数项截距此时有效自回归系数个数为rho.shape[-1] - 1init_dist初始值分布标量或向量形状应为(*shape[:-1], ar_order)若不匹配会被自动缩放默认Normal.dist(0, 100, shape...)并发出告警ar_orderAR 阶数可选。未指定时从rho最后一维长度推断ar_order rho.shape[-1] if constant else rho.shape[-1] - 1stepsAR 过程的步数steps 0可选可用 shape 替代。4.2 ar_order 的推断与约束源码_get_ar_orderpymc/distributions/timeseries.py通过constant_fold对rho.shape[-1]做常量折叠来推断阶数。若rho来自无静态形状的随机变量如Normal(size(5, 3))也可被识别为 3但完全动态的形状无法推断会抛出ValueError提示显式传入ar_order推断出的阶数小于 1 时同样报错。测试test_batched_sigma中就用pytensor.shared动态参数配合显式ar_orderar_order来构造模型。4.3 官方示例三阶 AR 带常数项文档内附的示例直接继承自 pymc/distributions/timeseries.pyimport pymc as pm # Create an AR of order 3, with a constant term with pm.Model() as AR3: # The first coefficient will be the constant term coefs pm.Normal(coefs, 0, size4) # We need one init variable for each lag, hence size3 init pm.Normal.dist(5, size3) ar3 pm.AR(ar3, coefs, sigma1.0, init_distinit, constantTrue, steps500)注意两点细节constantTrue时coefs长度 4 常数项 1 三阶滞后系数 3init_dist的 size 必须等于ar_order3 个滞后初始值init变量会在内部被克隆。4.4 logp 验证与批处理能力测试test_order1_logp与test_order2_logp用平凡的Normal模型对照验证了AR的 logp 数值正确性对一阶、二阶 AR其似然应等价于把每个时刻的期望写成muphi[0]*x[t-1]...的独立Normal的乘积tests/distributions/test_timeseries.py。test_batched_size、test_batched_rhos、test_batched_sigma进一步验证了rho、sigma、init_dist分别作为批维张量时与逐条建模的 logp 完全一致np.testing.assert_allclose这为在分层模型中复用同一组 AR 系数提供了保证。AR底层由AutoRegressiveRVpymc/distributions/timeseries.py实现通过pytensor.scan以tapsrange(-ar_order, 0)的延迟窗口迭代生成constant_term作为编译期属性保存在 Op 上其 logp 通过将系数与观测值卷积得到期望项再与Normal.dist(0, sigma)对比求和同时叠加init_dist对前ar_order个初始值的对数似然pymc/distributions/timeseries.py。五、GARCH11条件异方差波动率模型5.1 数学定义GARCH11用于刻画波动率随时间变化的时间序列如金融收益率模型为pymc/distributions/timeseries.py$$ y_t \sim N(0, \sigma_t^2) $$$$ \sigma_t^2 \omega \alpha_1 y_{t-1}^2 \beta_1 \sigma_{t-1}^2 $$其中误差方差 $\sigma_t^2$ 服从一个 ARMA(1,1) 结构。参数omegaomega 0基础方差mean variancealpha_1alpha_1 0自回归项系数ARCH 项beta_1beta_1 0且alpha_1 beta_1 1移动平均项系数GARCH 项initial_volinitial_vol 0初始波动率 $\sigma_0$。约束条件alpha_1 beta_1 1保证了条件方差过程的平稳性。5.2 使用示例import pymc as pm with pm.Model() as model: omega pm.Uniform(omega, 0, 1) alpha_1 pm.Uniform(alpha_1, 0, 1) beta_1 pm.Uniform(beta_1, 0, 1) # 需满足 alpha_1 beta_1 1 的平稳性约束可通过变换实现 y pm.GARCH11( y, omegaomega, alpha_1alpha_1, beta_1beta_1, initial_vol0.1, observedreturns, # 收益率观测序列 )steps与shape二选一提供支持长度源码在GARCH11.dist中要求二者必居其一否则抛出ValueError(Must specify steps or shape parameter)pymc/distributions/timeseries.py。5.3 实现与验证要点GARCH11的随机采样由GARCH11RV通过pytensor.scan完成pymc/distributions/timeseries.py每一步根据prev_y、prev_sigma递推new_sigma sqrt(omega alpha_1 * prev_y^2 beta_1 * prev_sigma^2)再从Normal(0, new_sigma)抽取new_y。其 logp 则用scan反向递推波动率序列后把每一步视为Normal(0, sigma_t)的独立观测并求和pymc/distributions/timeseries.py。测试test_logp手工用循环递推波动率验证GARCH11的 logp 与逐点Normal(0, vol)完全一致float64 精度 7 位小数tests/distributions/test_timeseries.pytest_batched_size验证了omega、alpha_1、beta_1、initial_vol四个参数各自支持批维扩展。六、EulerMaruyama随机微分方程的离散化分布6.1 原理与参数EulerMaruyama将随机微分方程SDE$$ dX_t f(X_t, \theta) , dt g(X_t, \theta) , dW_t $$用 Euler-Maruyama 格式离散化为pymc/distributions/timeseries.py$$ x_{t1} x_t \Delta t \cdot f(x_t, \theta) \sqrt{\Delta t} \cdot g(x_t, \theta) \cdot \epsilon_t, \quad \epsilon_t \sim N(0,1) $$参数dt离散化时间步长floatsde_fn可调用对象返回 SDE 的漂移项与扩散项系数(f, g)签名形如sde_fn(x, *sde_pars) - (f, g)sde_parsSDE 参数元组会以*args方式传入sde_fninit_dist初始值的标量分布形状应为(*shape[:-1],)不匹配会被自动缩放默认Normal.dist(0, 100, shape...)并告警steps离散化步数可选可用 shape 替代。6.2 官方测试中的经典用法以下取自测试test_linear_model的完整建模流程演示了用EulerMaruyama做状态空间反演tests/distributions/test_timeseries.py。它模拟 Ornstein-Uhlenbeck 型线性 SDEdx lam*x*dt sig2*dW观测带噪声然后用EulerMaruyama作为隐状态先验、Normal作为观测模型进行推断import numpy as np import pymc as pm lam, sig2, N, dt -0.78, 5e-3, 300, 1e-1 rng np.random.default_rng(42) # 生成 SDE 路径 def _gen_sde_path(sde, pars, dt, n, x0): xs [x0] wt rng.normal(sizen) for i in range(n): f, g sde(xs[-1], *pars) xs.append(xs[-1] f * dt np.sqrt(dt) * g * wt[i]) return np.array(xs) sde lambda x, lam: (lam * x, sig2) x _gen_sde_path(sde, (lam,), dt, N, 5.0) z x rng.standard_normal(sizex.size) * sig2 # 含噪观测 with pm.Model() as model: lamh pm.Flat(lamh) # 漂移系数的无信息先验 xh pm.EulerMaruyama( xh, dt, sde, (lamh,), stepsN, initvalx, init_distpm.Normal.dist(0, 10), ) pm.Normal(zh, muxh, sigmasig2, observedz) # 观测似然 with model: trace pm.sample(chains1, random_seedrng) ppc pm.sample_posterior_predictive(trace, modelmodel, random_seedrng)测试断言漂移系数lamh的 95% 后验区间包含真实值-0.78且后验预测区间覆盖观测值的比例大于 95%。测试中还出现了另一类 SDE 示例test_change_dist_size2其漂移项来自群体遗传学模型def sde2(p, s): N 500.0 return s * p * (1 - p) / (1 s * p), pm.math.sqrt(p * (1 - p) / N)可见sde_fn的返回值中扩散项g不必是常量可以依赖当前状态x与参数且支持用pm.math构造符号表达式。6.3 实现要点EulerMaruyamaRV同样基于pytensor.scanpymc/distributions/timeseries.py每步计算mu prev_y dt * f、sigma sqrt(dt) * g后从Normal抽取。其 logp 实现pymc/distributions/timeseries.py值得注意由于sde_fn是用户函数、不一定支持时间维广播源码将各 SDE 参数先扩展出[..., None]维度再对相邻两个时间点(x[..., :-1], x[..., 1:])计算条件正态对数似然并沿时间轴求和叠加init_dist对初始值的对数似然。七、贯穿所有分布的三条共性实现机制从源码可以看出时间序列家族共享同一套 PyMC 5 的现代分布注册机制理解它们有助于排查使用问题SymbolicRandomVariable 子图封装RandomWalkRV、AutoRegressiveRV、GARCH11RV、EulerMaruyamaRV都继承自SymbolicRandomVariable将一段由pytensor.scan构建的计算子图封装为一个带extended_signature的虚拟分布节点。它们统一使用pt.random.shared_rng(seedNone)管理随机数状态并通过update(node)方法返回 RNG 更新映射保证采样可复现、可沿时间维正确推进随机流。_logprob/_support_point/_change_dist_size三重注册每个分布都注册了自定义的 logp用于 MCMC 似然计算、support point用于初始点与后验矩估计以及尺寸变换函数配合change_dist_size与draw使用。测试中大量使用assert_support_point_is_expected与change_dist_size来验证这些注册行为。init_dist 的克隆与未命名约束所有分布都要求init_dist必须是通过.dist()API 创建的未命名分布未注册进模型且传入后会被克隆与外部变量解耦。若把已注册的分布如pm.Normal(init, ...)传入会触发check_dist_not_registered检查并抛错测试test_dists_not_registered_check验证了这一点tests/distributions/test_timeseries.py。八、常见问题与使用建议忘记传init_dist除GARCH11其初始值分布固定为Normal(0, initial_vol)外其余分布都会回退到Normal.dist(0, 100)之类的宽先验并发出UserWarning。生产建模中建议总是显式传入避免告警并让初始值先验更贴合领域知识。steps与shape冲突二者同时给出但长度不一致时会在图求值阶段触发 support_shape does not match respective shape dimension 断言建议只显式提供其一另一个交给框架推断。AR 的动态rho当rho的形状无法静态推断例如来自pytensor.shared或某些随机变量时ar_order推断会失败此时必须显式传入ar_order。GARCH 平稳性务必保证alpha_1 beta_1 1若用Uniform(0,1)等先验直接采样建议通过变换如 Dirichlet 构造把约束编码进模型否则生成的波动率过程可能发散。EulerMaruyama 的步长dt越小离散化误差越小但所需steps越大、计算成本越高sde_fn必须返回(漂移, 扩散)二元组扩散项建议保持非负。结语pymc.distributions.timeseries覆盖了贝叶斯时间序列建模中最高频的三类需求可定制创新项分布的随机游走RandomWalk及其三个预定义变体、经典自回归过程AR、条件异方差波动率GARCH11以及通过 Euler-Maruyama 接入任意 SDE 的通用方案EulerMaruyama。它们共享的SymbolicRandomVariable子图封装与 logp/support-point/size 注册机制让这些分布可以无缝参与pm.Model的构建、NUTS 采样与后验预测。结合本文给出的参数语义、官方示例与测试对照你可以直接把这些分布落地到自己的时序建模项目中。【免费下载链接】pymcBayesian Modeling and Probabilistic Programming in Python项目地址: https://gitcode.com/GitHub_Trending/py/pymc创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考