开始使用粒子Metropolis Hastings进行推理

在本教程中，我们利用L=100滞后来通过[IACT（θ（1:K））=1+2Xτ=1bρKτ来估计IACTθ（1:K）,其中，bρKτ=COR（φ（θ（K）），φ（θ（K+τ）））表示对φ的滞后-τ自相关的估计。对于图5所示的结果，这对应于三个参数中每一个的IACT{135，86，66}。调整提案可以通过调整提案分布来改善混合，以更好地适应后验分布。这需要对后验协方差进行估计，其作用类似于预处理矩阵P。通过试验运行byR>var（静态）来估计协方差非常简单，其中静态表示老化被丢弃后PMH生成的马尔可夫链的轨迹。对于第5节中考虑的问题，这将导致以下预处理矩阵xBP=10-4.1371-16 15-16 5 -1015-10 31,可用于形成byq给出的新的改进方案θθ（k）-1)= Nθ；θ（k）-1） ，2.562bP！。（28）Sherlock et al.（2015）提出了协方差矩阵的这种缩放，以在从高斯目标分布采样时最小化其影响。在图6中，我们展示了（26）中每对参数的边缘后验概率（颜色），以及（28）中的原始（顶部）和优化建议（底部）。TunedPropos更好地适应了后验结果，预计这将改善MarkovJohan Dahlin、Thomas B.Sch"on 29μφ-1.0-0.50.0.5 1.00.90 0.95 1.00 1.05μσv-1.0-0.5 0.0 1.00.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35φσv0的混合效果。90 0.95 1.00 1.050.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35μφ-1.0-0.50.0 0.5 1.00.90 0.95 1.00 1.05μσv-1.0-0.50.0.5 1.00.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35φσv0。90 0.95 1.00 1.050.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35图6：从PMH获得的u和φ（紫色）、u和σv（品红）以及φ和σv（绿色）的估计边缘后验值。

深色轮廓线（顶部/底部）分别表示第5节中的原始提案和从试点运行中估算的新（改进）提案。这两个建议都以后验平均值估计为中心。链条在图7中，我们展示了在使用调整后的方案时，通过第5节的实现获得的结果。此代码在函数示例4_sv中提供。由此得出的IACT估计值为{32,32,28}，与第5.3节中获得的结果相比，这在混合u和σV方面是一个明显的改进。其结果是，K（和计算成本）可以减少135/32=4.2，同时在估计中保持相同的方差。请注意，最大IACT是比较的，因为这些是限制因素。我们还可以比较图5和图7中后验估计的支持度。根据（27），改进的混合应该减少估计的渐近方差。然而，在这起案件中没有看到这样的改善。这可能是因为后验估计的差异很大程度上源于数据中的有限信息量。重新参数化模型另一种改进混合的方法是重新参数化模型，以获得无约束的参数，这些参数可以假定实线上的任何值。在第5节中，我们将φ和σvt约束到区域|φ|<1和σv>0，以获得稳定有效的SSM。

这导致PMH在非线性动力学模型100 200 300 400 500-4-2 0 2 4timelog-returns0 100 200 300 400-2-1 0 1 2timelog波动率估计μ后验估计-1.0-0.50.0 0 0.0 1.0 1.0 52500 3000 3500 4000-1.0-0.0 0 0 0 0 0.5 1.0迭代μ0 20 40 60 80-0.2 0.2 0 0 1.6 1.0。88 0.92 0.96 1.000 10 20 30 40 502500 3000 3500 40000.88 0.92 0.96 1.00迭代φ0 20 40 60 80 100-0.2 0.2 0.6 1.0迭代φσV后验估计0。0.1 0.2 0.3 0.40 5 10 15 202500 3000 3500 40000.0 0.1 0.2 0.3 0.4迭代σv0 20 40 60 80 100-0.2 0.2 0.6 1.0迭代σvFigure 7：上图：2012年1月2日至2014年1月2日期间，纳斯达克OMXS30指数95%置信区间内的每日对数收益率（深绿色）和估计对数波动率（橙色）。底部：从PMH获得的u（紫色）、φ（品红）和σv（绿色）的后验估计值（左）、马尔科夫链轨迹（中）和相应的ACF（右）。左图和中图中的虚线和实线分别表示参数后验平均值和参数前验值。Johan Dahlin，Thomas B.Sch"on 31如果提出了许多违反这些约束的候选参数，这会增加自相关性。通过对φ=tanh给出的模型进行重新参数化，可以缓解这个问题ψ, σv=exp, （29）使得ψ，∈ R是无约束参数。因此，PMH的目标改变为θ={u，ψ，}的后部。这种参数转换可以补偿，直到从θ的后验点获得样本。这是通过考虑（29）的雅可比矩阵来实现的，雅可比矩阵由ψtanh-1(ψ) =1 - ψ,log（）=-1.

（30）根据α（θ，θk）计算得到的接受概率-1） =min（1，p（θ）p（θ）k-1） ）bpNθ（y1:T|u）bpNθk-1（y1:T|uk-1)1.- (φ)1 - φk-1.σvσv，k-1.), （31）原始参数用于计算先验值和估计可能性。因此，该提议对θ进行操作，以提出新的候选参数θ，但应用（30）来获得θ，该θ用于计算接受概率（31）。在运行该实现之后，可以通过将（30）应用于θ后部的样本来获得θ后部的样本。为了实现这一点，我们需要改变对粒子过滤器的调用和接受概率的计算。完整的实现和代码在functionexample5_sv中提供。由此得出的参数后验平均估计为{-0.16,0.96,0.17}，标准偏差{0.24,0.02,0.06}和IACT{21,29,17}。请注意，maximumIACT甚至比调整提案时更小，这导致与第5节中的实施相比，速度提高了4.7。选择颗粒数量颗粒过滤器在PMH算法中起着重要作用，对混合有很大影响。如果对数似然估计值有噪声（粒子滤波器中的N太小），链往往会在多次迭代中卡住，这会导致性能不佳。这是一个事实的结果，有时pθ（y1:T）bpNθ（y1:T），这是由于估计量的随机性。因此，在合理的计算成本下平衡N和K以获得良好的性能是一个重要的问题。一个小N可能会导致我们需要把K取大，反之亦然。皮特等人（2012年）和杜塞特等人（2015年）对这个问题进行了调查。一个简单的经验法则是选择N，使得对数似然估计的标准偏差在1.0到1.7之间。

这些结果是在某些严格的假设下得出的，这些假设在实践中可能不成立。在这里，我们想在实际环境中研究这一经验法则。我们在调整方案时使用了该实现，并针对不同的N选择重新运行它。我们还使用后验平均值的估计作为参数，对每种N选择的粒子滤波器的对数似然估计值的标准偏差进行估计，并进行1000次独立的蒙特卡罗运行。在PMH开始之前，建议的分布使用一个试验运行32进行调整，以推断非线性动力学模型中的粒子数（N）50 100 200 300 400 500 bpnθ（y1:T）的标准偏差2.6 1.8 1.2 1.0 0 0.9 0.7接受概率（%）6 14 22 28 29 33最大IACT 165 139 128 92 114 97每次PMH迭代时间0.03 0.05 0.09 0.12 0.15 0.1917 18表3：对数似然估计的标准偏差、可接受概率、θ的最大IACT、PMH每次迭代的计算时间和变化N的基准量（见正文）。马尔可夫链在后验平均值的估计中初始化。IACT的计算采用L=100滞后。表3给出了与颗粒过滤器和PMH相关的一些数量，因为N是变化的（保持所有其他参数不变）。根据Doucet等人（2015）提出的经验法则，N的最佳选择在100到300之间，取决于提案的选择。用户的目标通常是最小化PMH的计算时间，以便从后验数据中获得一定数量的不相关样本。因此，合适的基准量是最大IACT乘以一次迭代的计算时间。记住，IACT是两个不相关样本之间的估计迭代次数。

因此，基准数量（表3的最后一行）是对每个不相关样本的计算时间的估计。结果有点不确定，这是因为IACT的估计非常困难，因此有噪声。然而，经验法则似乎对这种型号有效，将N设置为100似乎是一个不错的选择。6.4. 在实际应用中，包括几何信息在内，当参数p的数量超过5时，或者如前面所讨论的，当链初始化远离高后验概率区域时，我们通常会遇到性能问题。即使使用Sherlock等人（2015）的经验法则来调整提案，也会出现这个问题。原因是，问题在于使用随机游走，这是已知的对高维空间探索较差的方法。为了缓解这个问题，可以考虑关于后验概率的几何信息。Girolami和Calderhead（2011）展示了如何利用对数后验概率的梯度和hessian将马尔可夫链引导到后验概率较高的区域。在本文中，这是由黎曼流形上的扩散过程驱动的。然而，在优化中可以找到一个更简单的类比，我们可以在方案中利用噪声梯度或牛顿更新。梯度信息有助于引导马尔科夫链在磨合期间朝向感兴趣的区域，也有助于在磨合阶段后将其保持在该区域。Hessian信息可用于重新缩放参数，使其更接近各向同性，这大大简化了采样。在Dahlin等人（2015b）中，作者展示了如何在PMH算法中使用这种类型的建议。我们指的是仅使用梯度信息asPMH1（第一阶）的提案。

利用梯度和海森信息的方案称为PMH2（二阶）。这里的挑战是获得梯度和Hessian的良好估计，对于非线性SSM，这两种方法在分析上都是难以解决的。FurJohan Dahlin，Thomas B.Sch"on 33thermore，这些量的计算通常需要使用计算成本高的粒子平滑器，这在PMH算法中是禁止的。为了缓解这一问题，他们建议利用更快但更不准确的固定滞后粒子平滑，并在非正定义时将Hessian估计正规化。在Dahlin等人（2015c）中，作者描述了一种基于噪声拟牛顿更新的拟牛顿PMH2方案（qPMH2），该方案不需要任何Hessian信息，但基于梯度信息构造了Hessian的局部近似。Nemeth等人（2016年）从理论上研究了PMH1和类似算法，该算法根据后验协方差的估计，提供了一个调整步长的Tumb规则。6.5. 控制变量控制变量是标准MonteCarlo中一种常见且有用的方差缩减技术，也可应用于PMH等MCMC算法。Mira等人（2013年）、Papamarkou等人（2014年）、Mijatovi'c和Vogrinc（2018年）以及Dellaportas和Kontoyiannis（2012年）在这些方面为MH做了一些有趣的工作，这些工作对于PMH的实施应该是相当前瞻性的。然而，据本文作者所知，尚未对PMH的控制变量进行适当调查，但Dahlin（2016）第4.3章根据Papamarkou等人（2014）给出了一些令人鼓舞的初步结果。相关软件与当前教程相关的软件包有很多，它们实现（i）PMH和/或（ii）粒子过滤和SMC。

这些可用于在新SSM中快速执行Bayesian参数推断，或作为希望创建自己的PMH实现的用户的构建块。软件LibBi（Murray 2015）提供了一个在SSM中使用串行和并行硬件进行贝叶斯推理的平台，特别关注高性能计算。它是用C++编写的，允许用户使用类似于JAGS（Plummer 2003）或BUGS（Lunn等人，2009）的建模语言定义新模型。这使得用户能够在许多不同类型的SSM中使用PMH和SMC快速解决参数推断问题。在第5节中介绍的SV模型中使用LibBi进行推理的示例可通过主页获得，这有助于学习如何使用该软件。此外，R和Octave的接口（Eaton等人，2017年）可分别通过RBi通过CRAN和OctBi通过Libbi主页获得。该软件适用于希望利用多核计算能力在大型数据集上快速进行推理的用户。另一种选择是Biips（Todeschini et al.2014），它是在C++中实现的，通过包RBIPS和工具箱Matbiips分别与R和MATLAB接口。Biips的功能与LibBi相似，关键区别在于后者的并行计算。因此，Biips允许用户使用BUGSlanguage定义模型，并使用PMH进行推理。连接到Biips的主页包含一些实现第5节模型的示例代码。http://libbi.org/https://biips.github.io/34开始使用PMH进行非线性动力学模型的推理。R package pomp（King等人，2016）包括一系列令人印象深刻的SSM不同推理方法。其中包括粒子滤波、PMH和近似贝叶斯计算（ABC；Marin等人。

2012). 对于经常使用R的读者来说，这个软件是一个很好的选择，他们想探索一些不稳定计算的最新发展。它还包含许多预先定义的流行病学模型。pomp中的模型规格与LibBi和Biips相比要复杂一些，这对一些用户来说可能是一个缺点。概率编程语言（PPL）是统计建模软件的最新贡献。他们用随机分支、循环和递归扩展图形模型。允许在SSM中使用PMH进行推理的三种流行PPL是英国国教（Tolpin等人，2016年）、Venture（Mansinghka等人，2014年）和概率C（Paige and Wood，2014年）。圣公会是这三种语言中最发达的一种，它作为一个独立的应用程序运行。在英国国教中实现PMH需要一些额外的编码，例如LibBi。然而，它是一个更通用的框架，用于对更大类别的模型进行推理。SMCTC（Johansen 2009）提供了用于粒子过滤和更通用SMC算法的C++模板。该模板通过软件包CPPSMC（Eddelbuettel and Johansen 2018）与R的接口进行补充，并通过软件vSMC（Zhou 2015）扩展到C++中的并行计算。这些模板需要用户进行额外的实现才能进行推断。它们不允许简单地指定新模型，或直接使用PMH进行推断。然而，它们为用户提供了对粒子过滤器的完全控制，从而允许添加扩展并根据特定问题定制算法。Lee等人（2010）提供了一个使用CUDA的SMC并行实现，源代码可在线获取。Casarin等人（2015）开发了一个用于并行运行SMC算法的MATLAB工具箱，名为DeCo。

DeCo利用SMC组合与经济学预测相关的密度，但可能会推广到其他应用中。最后，我们想简要介绍一些附加软件。Python Library PyParticleTest（Nordh 2017）提供了使用不同类型的粒子过滤器进行状态估计的功能。它还包括通过期望最大化（EM；Dempster等人，1977；McLachlan和Krishnan，2008）算法使用最大似然法进行参数估计。Bos（2011）在Ox（Doornik 2007）中实施MH，以估计第5节中SV模型的状态和参数。在此设置中，MH用于估计状态和参数，因此不需要粒子过滤器。Wilkinson（2011）随附的R软件包smfsb（Wilkinson 2018）包含在应用中使用PMH的演示代码。本教程第一作者的GitHub Repository包含PMH1、PMH2、qPMH2（Dahlin等人2015b，c）和伪边际MH的Python代码，以及相关随机变量（Dahlin等人2015a）。这段代码与本教程提供的Python代码非常相似，对于想要尝试更高级的PMH实现的感兴趣的读者来说，这是一个很好的起点。http://www.oxford-man.ox.ac.uk/gpuss/cuda_mc.htmlhttps://www.github.com/compops.Johan达林，托马斯·B·舍恩358。结论我们描述了非线性SSMs中贝叶斯参数推断的PMH算法。这包括粒子过滤器，因为它在PMH中起着重要作用，并提供了非负无偏的可能性估计。此外，我们还利用合成数据和真实数据，在LGSS模型和SV模型中应用PMH进行推理。我们还确定并讨论了与PMH有关的各种实际问题。