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

这种递归是通过以下方式实现的：<EstimateLikelion>=PredictiveLikelion<-maxWeight+log（sumWeights）-log（noParticles）LogLikelion<-LogLikelion+predictiveLikelihood后验分布的估计通过将对数似然的估计插入（12）中得出<estimateState>时间t处的潜在状态xt可通过（15）进行估计，给定观测值y1:（t+1），其对应于估计器bxnt=NNXi=1x（i）t，（24）通过以下方式实现：<estimateState>=xHatFiltered[t]<-均值（粒子[，t]）。注意，这对应于权重w（i）t=1/N，其对应于根据（23）中的表达式，所有粒子都是相同的。然而，这并不是PMH在非线性动力学模型粒子数（N）10 20 50 100 200 500 1000对数偏差推断中的情况-3.70-3.96-4.57-4.85-5.19-5.67-6.08log-MSE-6.94-7.49-8.72-9.29-9.91-10.87-11.67表1：不同N的过滤状态的对数偏差和对数MSE。一般来说，估算器（24）是在粒子过滤器的传播和称重步骤中使用最佳选择的结果。这种选择对应于一种特殊类型的算法，称为完全自适应粒子过滤器（faPF；Pitt和Shephard，1999），更多细节请参见第3.3节。在faPF中，我们在重采样步骤和构造πt（xt）的经验近似时使用单独的权重。在第5节中，我们介绍了另一种类型的粒子过滤器，以及BXNT的另一种估计器，该估计器可用于大多数TSSM<returnEstimates>particleFilter的输出由以下内容返回：<returnEstimates>=list（xHatFiltered=xHatFiltered，loglikelion=loglikelion）3.2。数值说明我们利用粒子过滤器的实现来估计过滤状态，并研究有限N的这种估计的性质。

完整的实现可以在脚本/函数示例1中找到。我们使用图2所示的数据，该数据是在本节开头生成的。将粒子滤波器的估计值与卡尔曼滤波器的相应估计值进行比较。后者源自于利用高斯分布的特性精确求解（5），这仅适用于LGSS模型。在图3的中间，我们展示了卡尔曼滤波器的最佳状态估计与使用N=20个粒子的粒子滤波器的估计之间的差异。两种可供选择的精度度量是状态估计的偏差（绝对误差）和均方误差（MSE）。这些是根据toBias计算的bxNt=TTXt=1bxNt-通道管理系统, MSEbxNt=TTXt=1bxNt-通道管理系统,其中bxt表示卡尔曼滤波器获得的最佳状态估计。在表1和图3的底部，我们给出了N不同值的偏差和MSE的对数。我们注意到，当N增加时，偏差和MSE迅速减少。因此，我们得出结论，对于这个模型，N不需要很大，BxNt才能准确。3.3. 展望：设计选择和概括在本教程继续实施PMHalgorithm之前，我们做一个简短的回顾。在本节中，我们将为颗粒过滤器的重要改进和进一步研究提供一些参考。请注意，有关粒子过滤器的科学文献非常丰富，这只是对主题和参考文献的一小部分自然偏颇的选择。Johan Dahlin，Thomas B.Sch"on 150 50 100 150 200-6-4-2 0 2 6时间观测0 50 100 150 200-0.10-0.05 0.00 0.10状态估计中的时间误差0 200 400 600 800 1000-7-6-5-4-3无。颗粒（N）对数bias0 200 400 600 800 1000-12-11-10-9-8-7-6个。

粒子（N）对数均方误差图3：顶部和中部：LGSS模型的一组模拟观测值（顶部），以及使用N=20的粒子过滤器进行的潜在状态估计误差（中部）。下图：当改变粒子数N时，粒子过滤器的估计对数偏差（左）和对数均方误差（右）。16开始在非线性动力学模型中推断PMH。Doucet和Johansen（2011）、Cappé等人（2007）给出了更全面的粒子过滤器调查和教程，Arulampalam等人（2002年）和Ala Luhtala等人（2016年）也提供了伪代码。这些文件中讨论的更改可以很容易地合并到本节中开发的实现中。在第3.1节中，我们使用faPF（Pitt and Shephard 1999）来估计过滤状态和可能性。该算法利用了粒子建议和加权函数的最佳选择，对应于能够从p（xt+1 | xt，yt+1）中采样并计算p（yt+1 | xt）。这对于许多SSM来说是不可能的。一种可能性是创建所需量的高斯近似，见Doucet等人（2001）和Pitt等人（2012）。然而，这些方法依赖于拟牛顿优化，当N较大时，该优化具有计算上的可行性。有关讨论和伪代码，请参见Arulampalam等人（2002）。另一种可能性是在粒子过滤器中混合使用提议和加权函数，如Kronander和Sch"on（2014）所述。这类过滤器基于多重重要性抽样，通常用于计算机图形学等领域。一种不同的方法是使用一个额外的粒子过滤器来近似完全适应的提案，从而形成一个嵌套结构，其中一个粒子过滤器与另一个粒子过滤器一起使用，以构建该粒子过滤器的提案分布。

结果结构称为嵌套顺序蒙特卡罗（SMC），详情见Naesseth等人（2015）。嵌套SMC结构使得考虑Nx较大的状态空间成为可能，这是其他类型的粒子过滤器难以解决的问题。引导式粒子过滤器（bPF）是默认选择，因为它可以用于大多数TSSMS。bPF的一个缺点是，当状态空间nx的维数超过约10时，其统计精度较差。在此设置中，例如，状态推断需要faPF或nestedSMC。然而，在某些情况下，有一些方法可能会提高bPF的性能，应考虑有效实施。这些包括：（i）并行实现，（ii）更好的重采样方案，（iii）利用SSM的线性子结构，以及（iv）使用准蒙特卡罗。粒子过滤器（i）的并行实现是一个正在进行的研究主题，但Lee等人（2010年）和Lee and Whiteley（2016年）报告了一些令人鼓舞的结果。在本教程中，我们将在粒子过滤器中使用多项式重采样。交替采样方案（ii）可用于降低估计值的方差，有关一些比较，请参见Holet al.（2006）和Douc and Cappé（2005）。一般来说，建议对颗粒过滤进行系统重采样。另一个可能的改进是卡尔曼滤波和粒子滤波（iii）的结合，如果模型在某些状态下是条件线性和高斯的，这是可能的。然后利用卡尔曼滤波来估计这些状态，并对剩余状态进行粒子滤波，同时保持线性状态与卡尔曼估计值一致。这些类型的模型在工程中很常见，像这样的Rao Blackwellization方案可以大幅降低方差。例如，见Doucet等人。

（2000），Chen and Liu（2000）和Sch"on等人（2005），以获取更多信息和比较。准蒙特卡罗（iv）基于所谓的准随机数，由确定性递归生成，与标准随机数相比，它能更好地填充空间。这些在标准蒙特卡罗中很有用，可以减少估计中的方差。Gerber和Chopin（2015）讨论了准蒙特卡罗在粒子滤波中的应用。最后，粒子滤波是SMC方法的一个实例（Del Moral等人，2006年），Johan Dahlin，Thomas B.Sch"on 17代表了一类基于重要性抽样的通用算法。SMC可以用于许多统计模型的推理，是MCMC的补充/替代。一个特别有趣的成员是SMC算法（肖邦等人，2013；Fulop和Li2013），它基本上是一个两级粒子过滤器，类似于嵌套SMC。外部粒子过滤器维持以πθ（θ）为目标的粒子系统。内部粒子过滤器以πt（xt）为目标，并针对外部过滤器中的每个粒子运行，因为这包含θ值的假设。内部过滤器的似然估计用于计算外部过滤器的权重。因此，SMC是PMH参数推断的替代方案，请参见Svensson和Sch"on（2016）进行比较。4.估计线性高斯SSM中的参数本教程现在继续进行主事件，PMH算法根据第2节提供的大纲实施。第3.1节中的粒子过滤器实现用于近似PMH算法中的接受概率（13）。我们将LGSS模型作为一个玩具示例来说明实现，并在本节末尾提供PMH文献的概述。先验知识要求能够进行贝叶斯参数推断。

本节的目的是估计θ=φ的后验参数，同时保持{σv，σe}固定在其真值上。因此，我们选择先验的p（φ）=tn(-1,1）（φ；0,0.5），以确保系统稳定（即，XT的值对每t有界）。区间z内平均μ、标准差σ的截断高斯分布∈ [a，b]由N[a，b]（z；u，σ）=I（a<z<b）N（z；u，σ）定义，其中I（s）表示指示函数，如果s为真，则值为1，否则为零。4.1. 实现粒子Metropolis Hastings函数particleMetropolis中提供了PMH算法实现的完整源代码。它的代码框架由以下公式给出：particleMetropolisHastings<-函数（y，initialPhi，sigmav，sigmae，noParticles，initialState，noIterations，stepSize）{<initialization>for（k in 2:noIterations）{<proposeParameters><computeAcceptProbability><acceptRejectStep>}phi}这个函数有输入：y（带有T+1个观测值的向量），initialPhi（φ（φ（0）φ的初始值），sigmav，sigmae（参数{σv，σe}）、noParticles、initialState、noIterations（编号PMH迭代K）和stepSize（建议书中的步长）。函数particleMetropolishings的输出为φ（1:K），即根据参数后验πθ近似分布的相关样本。18从PMH开始在非线性动力学模型中进行推理<初始化>我们分配一些变量来存储马尔可夫链Phi的当前状态，建议的状态Phi，提出了当前对数似然对数似然和建议的对数似然对数似然。此外，如果建议的参数被接受，我们将二进制变量ProposedPhiAccept分配为1，否则分配为0。最后，我们使用参数{φ（0）、σv、σe}运行粒子滤波器来估计初始可能性。

初始化通过以下方式实现：<初始化>=phi<-矩阵（0，nrow=noIterations，ncol=1）建议的phi<-矩阵（0，nrow=noIterations，ncol=1）对数似然<-矩阵（0，nrow=noIterations，ncol=1）建议的logLikelihood<-矩阵（0，nrow=noIterations，ncol=1）phi[1]<-初始phitheta<-c（phi[1]，sigmav，sigmae）outputPF<-particleFilter（y，θ，noParticles，initialState）对数似然[1]<-outputPF$Log似然<proposeParameters>如第2节所述，提案分布有很多选择。在本教程中，我们使用byq给出的高斯随机游走方案θθ（k）-1)= Nθ; θ（k）-1), , （25）在哪里 > 0表示随机游动的步长，即增量的标准偏差。建议步骤通过以下方式实现：<proposeParameters>=phiProposed[k]<-phi[k-1]+步长*rForm（1）<computeAcceptProbability>由于（25）在θ（即q）上是对称的，因此可以简化接受概率（13）θθ（k）-1)= Qθ（k）-1)θ,结果（13）可以表示为α（θ，θ（k-1） ）=min（1，explog“p（θ）p（θ（k-1） ）#+log“bpNθ（y1:T）bpNθ（k-1） （y1:T）#！）。对数似然和对数先验用于避免数值精度的损失。其前提是α（θ，θ（k-1） ）=当|φ|>1时为0。因此，在这种情况下，粒子过滤器不会运行以节省计算，而接受概率设置为零以确保θ被弹出。接受概率的计算通过以下方式实现：<computeAcceptProbability>=if（abs（phiaproved[k]）<1.0）{theta<-c（phiaproved[k]，sigmav，sigmae）outputPF<-particleFilter（y，theta，noParticles，initialState）Johan Dahlin，Thomas B。

Sch"on 19logLikelihoodProposed[k]<-outputPF$Loglikelion}priorPart<-dnorm（phiProposed[k]，log=TRUE）priorPart<-priorPart-dnorm（phi[k-1]，log=TRUE）likelihoodDifference<-LoglikelihoodProposition[k]Loglikelikelization[k-1]acceptProbability<-exp（priorPart+likelihoodDifference）acceptProbability<-acceptProbability*（abs）（abs（phiProposed[k]）<1.0）<Acce，我们需要决定是否接受提议的参数。这是通过内置的R命令runif模拟[0,1]上的均匀随机变量ω来实现的。如果ω<α，我们接受θθ、 θ（k-1)将其及其对应的对数似然值存储为马尔可夫链的当前状态。否则，我们将保留状态的当前值和上一次迭代的对数可能性。接受/拒绝步骤由以下步骤实现：<acceptRejectStep>=uniformRandomVariable<-runif（1）if（uniformRandomVariable<acceptProbability）{phi[k]<-phi建议的[k]对数似然[k]<-logLikelihood建议的[k]建议的[k]phi接受的[k]<-1]对数似然[k]<-loglikelike[k]建议的phi接受的[k]<-0}4.2。数值说明我们利用图2中的数据，利用粒子群估计θ=φ。完整的实现和代码可在函数/脚本示例2_lgss中找到。我们将马尔可夫链初始化为θ=0.5，并利用K=5000次迭代（无迭代），将第一次Kb=1000次迭代（无迭代）作为磨合丢弃。也就是说，我们只使用最后4000个样本来构造参数后验的经验近似值，以确保马尔可夫链实际上已经达到其平稳状态，有关更多信息，请参见第6节。在图4中，使用不同的步长显示了三次PMH运行 = 0.01（左）、0.10（中）和0.50（右）。

该病例的后验平均值估计值为isbφ=0.66 = 0.10. 这是通过自相关函数（ACF）计算得出的R>mean（φ[Noburniterations:noIterations]），我们看到 影响马尔可夫链中的相关性。一个不错的选择 因此，获得高效的算法非常重要。我们将在第6.3.20节“PMH入门”中继续讨论这个问题，以便在非线性动力学模型φ后验估计0中进行推理。50 0.60 0.70 0.800 2 4 6 8 10 12φ后验估计0。50 0.60 0.70 0.800 2 4 6 8 10 12φ后验估计0。50 0.60 0.70 0.800 2 4 6 8 10 121000 1200 1400 1800 20000.4 0.5 0.6 0.7 0.8迭代φ1000 1200 1200 20000.4 0.5 0.6 0.7 0.8迭代φ1000 1200 1600 1800 20000.4 0.5 0.6 0.7 0.8迭代φ0 10 20 30 40 50-0.2 0.0 0.2 0.4 0.6 0.8迭代φ10 20 20 20 20 40 50-0.0.0.2 0.0 0 0.0 0 0 0.2 0 0 0.2 0.0 0 0.2 0.0 0 0.4 0.0 0 0 0 0 0.4 0.6 0.0-0.0.0.0迭代图1使用PMH算法估计LGSS模型中的πθ[φ]，使用三种不同的步长： = 0.01（左）、0.10（中）和0.50（右）。上图：πθ的估计以直方图和核密度估计（实线）表示。中间：老化后1000次迭代的马尔可夫链状态。底部：马尔科夫链的估计ACF。上图和中图中的虚线表示对后验概率的估计。底图中的虚线表示ACFcoe效率的95%置信区间。约翰·达林，托马斯·B。

Sch"on 21观察次数（T）1020501000200 500估计后验均值0.596 0.794 0.765 0.727 0.696 0.719估计后验方差0.040 0 0.013 0.009 0.006 0.003 0.001表2：变化T时的估计后验均值和方差。我们注意到，参数估计值与真实值0.75略有不同，并且参数后验估计值的不确定性相当大。这是因为样本量T（和a fite K）相对较小。根据Bayesianestimator的渐近理论，我们知道后部质量倾向于集中在真实参数T（和K）周围。我们在表2中举例说明了这一点，当T增加时，我们使用相同的设置来估计后验均值和方差。这项小型研究表明，真实参数确实是通过极限内的后验平均估计恢复的。然而，这种收敛的速度是由模型决定的，因此不可能给出任何关于T需要多大才能达到一定精度的通用指南。4.3. 展望：概括我们用这一节为PMH算法的重要改进和进一步研究提供一些参考。关于粒子过滤器，与PMI有关的科学文献数量巨大且增长迅速。因此，这只是对近期发展的一小部分有偏见的选择。有关SSM中参数推断的更广泛观点，请参见Kantas等人（2015）和Sch"on等人（2015）的调查。如第2节所述，PMH算法属于精确近似或伪边缘算法家族（Andrieu和Roberts，2009）。在这里，粒子过滤器用于估计目标，但也可以为此目的使用重要性抽样等方法。PMH也是所谓的粒子MCMC（PMCMC；Andrieuet al。