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

当N时，可以证明经验近似收敛于真实分布→ ∞ 在一些规律性假设下。然后，通过将经验滤波分布（10）插入（4），通过分解（3）估计可能性。这导致估计量bpNθ（y1:T）=bpNθ（y）TYt=2bpNθ（yt | y1:T）-1） =TYt=1“NNXi=1v（i）t#，（11）其中，非归一化参数后验值的逐点估计由bγNθ（θ）=p（θ）bpNθ（y1:t）给出。（12）支持MH的理论细节被故意忽略。为了使遍历定理成立，马尔可夫链需要是遍历的，这意味着马尔可夫链应该能够探索整个空间，而不是陷入其中的某些部分。感兴趣的读者可以在中找到更多关于理论细节的信息，例如MH的Tierney（1994）和Meyn and Tweedie（2009）以及PMH的Andrieu等人（2010）。我们再次省略了许多关于粒子过滤器的重要理论细节。有关这些以及与该算法相关的许多其他重要主题的更多信息，请参阅Crisan and Doucet（2002）和Doucet and Johansen（2011）。Johan Dahlin，Thomas B.Sch"on 7在这里，假设先验可以通过闭合形式表达式逐点计算。验收概率（7）可以通过插入（12）来近似计算，从而得出αθ、 θ（k-1)= min（1，bγNθ（θ）bγNθ（θ（k）-1） ）qθ（k）-1)θQθθ（k）-1))= min（1，p（θ）p（θ）k-1） bpNθ（y1:T）bpNθ（k）-1） （y1:T）qθ（k）-1)θQθθ（k）-1)). （13） 只有当bpnθ（y1:T）是似然的无偏估计量时，这种近似才有效。幸运的是，当使用粒子滤波器作为其似然估计量时，这种情况是无偏的（Del Moral 2013；Pitt et al。

2012）对于任何N≥ 1.与参数推理问题一样，人们的兴趣通常在于计算一个行为良好的测试函数的期望值→ 关于πt（xt）的Rn~n，由πt[~n]给出，Eh~n（xt）y1:ti=ZXа（xt）πt（xt）dxt，其中再次选择测试函数а（xt）=xt对应于计算边缘过滤分布的平均值。由于πt（xt）未知，这一预期难以实现。同样，我们用经验近似（10）提供的估计值代替它，结果是：inbπNt[~n]，ZX~n（xt）bπNt（dxt）=NXi=1w（i）tZX~n（xt）δx（i）t（dxt）=NXi=1w（i）t~nx（i）t, （14） 对于大约0≤ T≤ T由狄拉克δ函数的性质决定。在接下来的章节中，我们特别感兴趣的是边际滤波分布平均值的估计值bxNt，bπNt[x]=NXi=1w（i）tx（i）t，（15），这被称为滤波状态估计。在某些假设下，bπNt[~n]的性质与PMH算法中的估值器类似，更多信息请参见Crisan和Doucet（2002）以及Doucet和Johansen（2011）。因此，我们证明了当N→ ∞) 误差为高斯分布，方差减少1/N。8开始使用PMH进行非线性动力学模型2的推理。5.展望：伪边缘算法在前面关于PMH的讨论中采用的观点是，它是一种MH算法，使用粒子过滤器来近似其他难以解决的接受概率。另一个更普遍的观点是将PMH视为伪边缘算法（Andrieuand Roberts，2009）。

在这种类型的算法中，会引入一些辅助变量来简化计算，但在算法运行期间，这些辅助变量会被边缘化。这导致PMH算法可以被视为在非标准扩展空间Θ×U上运行的标准MH算法，其中Θ和U分别表示参数空间和辅助变量空间。由此产生的扩展目标由πθ，u（θ，u）=bπNθ（θ| u）m（u）给出。（16） 这里，参数“后验”增加了u∈ U，表示密度为m（U）的多变量变量。从上面的讨论中，我们知道u可以通过（12）的粒子滤波来构造bπNθ（θ| u）的逐点估计。基于粒子滤波器的似然估计的无偏性给出了hbγNθ（θ| u）i=ZUbγNθ（θ| u）m（u）du=γθ（θ）。（17） 这意味着，可以通过边缘化所有辅助变量u来恢复非标准化参数后验γθ（θ）。在实现中，这导致可以忽略u的样本值。因此，我们不会将其存储在后续实施中，而只保留θ（和xt）的样本。我们通过推导（16）中的伪边缘算法得出结论。θ和u的建议是用formq选择的θ、 uθ（k）-1） ，u（k）-1)= qθθθ（k）-1） ，u（k）-1)曲Uu（k）-1), （18） 这是选择asqθ的两个方案的乘积θθ（k）-1） ，u（k）-1)= qθθθ（k）-1), 曲Uu（k）-1)= m（u）。（19） 这对应于u的独立提案和不包括deu的θ提案。其他选项是当前研究的主题，请参见第4.3节了解一些参考资料。

选择目标和方案的结果接受概率由α（θ，u，θ（k）给出-1） ，u（k）-1） ）=min（1，πθ，u（θ，u）πθ，u（θ（k）-1） ，u（k）-1） q（θ（k）-1） ，u（k）-1） |θ，u）q（θ，u |θ（k）-1） ，u（k）-1） ）=min（1，bγNθ（θ| u）bγNθ（θ（k）-1） |u（k）-1） qθ（θ（k）-1） |θ）qθ（θ|θ（k）-1） ）），注意，得到的接受概率与（13）中的相同。这些论据提供了PMH从targetdistribution生成样本的证据，并由Flury和Shephard（2011）首次提出。有关PMH有效性的更正式的处理和证明，请参见Andrieu和Roberts（2009）和Andrieuet al.（2010）。Johan Dahlin，Thomas B.Sch"on 90 50 100 150 200 250-4-2 0 2 4 6时间潜伏期xt0 50 100 150 200 250-4-2 0 2 4时间观测yt0 5 10 15 20 25-0.2 0.0 0.2 0.4 0.6 0.8 1.0 YTGF图2：LGSS模型的模拟数据，具有潜伏期（橙色）、观测值（绿色）和观测值的自相关函数（浅绿色）。在线性高斯SSMWe中估计状态现在准备好开始基于第2节中的材料实施PMH算法。为了简化演示，本节讨论了如何估计线性高斯状态空间（LGSS）模型的滤波分布{πt（xt）}Tt=0。这些分布可用于计算xn0:TandbpNθ（y1:T），即过滤状态的估计值和可能性的估计值。参数推断问题将在下一节中讨论。byx给出了所考虑的特定LGSS模型~ δx（x），xt | xt-1.~ Nxt；φxt-1，σv, yt | xt~ Nyt；xt，σe, （20） 其中参数用θ={φ，σv，σe}表示。φ ∈ (-1，1）决定状态的持续性，而σv，σe∈ R+分别表示状态转移噪声和观测噪声的标准差。高斯密度用N（x；u，σ）表示，平均u和标准偏差σ>0。

图2显示了使用θ={0.75,1.00,0.10}和x=0的T=250个观测值的模型模拟数据记录。数据生成的完整代码可在generateData中找到。3.1. 实现粒子过滤器功能particleFilter中提供了实现粒子过滤器的完整源代码。其代码框架由以下公式给出：particleFilter<-函数（y，θ，noParticles，initialState）{<initialization><initializeStates>for（t in 2:t）{<resamplepoticles><propagateParticles><weightsparticles>10开始用PMH在非线性动力学模型中进行推理<normalizeWeights><estimatelikelibility><estate state><returnEstimates>}函数particles Filter有输入：y（带有T+1观测值的向量）、θ（参数{φ，σv，σe}）、无粒子和初始状态。输出将是过滤状态的估计值Bxn0:T和可能性bpn（θ| y1:T）的估计值。注意particleFilter在t=2，T，对应于时间指数T=1，T-1in（20），因为向量元素的编号从R中的1开始。迭代在时间索引T处终止-1因为每次迭代都需要未来的观测yt+1。也就是说，我们假设新的观测值在传播和加权步骤之间到达算法。因此，如果可能的话，在加权步骤中使用信息是有意义的<初始化>通过为变量分配内存来初始化粒子过滤器：粒子、索引、权重和归一化权重。它们分别对应于粒子、它们的祖先、非标准化权重和标准化权重。此外，两个变量xHatFiltered和loglikelike分别分配给storebxN0:TandlogbpN（θ| y1:T）。

这一切都是通过以下方式实现的：<初始化>=T<-长度（y）-1粒子<-矩阵（0，nrow=noParticles，ncol=T+1）取消存储索引<-矩阵（0，nrow=noParticles，ncol=T+1）权重<-矩阵（1，nrow=noParticles，ncol=T+1）归一化权重<-矩阵（0，nrow=noParticles，ncol=T+1）xHatFiltered<-矩阵（0，nrow=T，logncol=1）似然<-0（20） ，我们有μθ（x）=δ（x），所以所有粒子最初都设置为x（1:N）=x=0，所有权重都设置为w（1:N）=1/N（因为所有粒子都是相同的）。该操作通过以下方式实现：<initializeStates>=AncestorIndexs[，1]<-1:NoParticlesPaticles[，1]<-InitializeStateXhatFiltered[，1]<-InitializeStateNormalizedWeights[，1]=1/noParticles<resampleParticles>过滤器中的粒子对应于潜伏期数值的许多不同假设。权重表示给定粒子在模型下产生观测的概率。在重采样步骤中，该信息用于将粒子过滤器的计算效果集中到状态空间的相关部分，即重量较大的粒子。这是通过随机复制重量较大的粒子并丢弃重量较小的粒子来实现的，这样粒子总数始终保持不变。请注意，Johan Dahlin，Thomas B.Sch"on 11重采样步骤是无偏的，因为重采样粒子的预期比例由粒子权重给出。这一步很重要，因为粒子系统在经过几次迭代后，将只由一个唯一的粒子组成。这将导致（14）中出现较大差异。在我们的实现中，我们使用多项式重采样，这也被称为带替换的加权自举。

这个过程的输出是每个粒子i的所谓的指数a（i）t，可以解释为时间t时粒子i的父指数。对于每个i=1，N、 祖先指数从多项式（分类）分布中取样，PHA（i）t=ji=w（j）t-1，j=1，N.重采样步骤是通过调用函数sample来实现的：<重采样粒子>=新祖先<-sample（noParticles，replace=TRUE，prob=normalizedWeights[，t-1]）ancestorindex[，1:（t-1）]<-ancestorindex[new祖先，1:（t-1）]ancestorindex[，t]<-newAncestors其中，生成的祖先索引a（1:N）作为新祖先返回，并存储在存储条件中记账。请注意，粒子谱系也会与粒子同时重新采样。所有这些都是为了随着时间的推移追踪粒子系统的谱系。也就是说，粒子系统的整个历史<传播粒子>通过使用模型向前传播粒子系统，更新关于状态的假设。这相当于从粒子建议分布中取样，以生成新粒子x（i）tbyx（i）t | xa（i）tt-1.~ pθx（i）txa（i）tt-1，yt, （21）其中来自先前粒子xa（i）tt的信息-1和电流测量值可包括在内。粒子方案有几种不同的选择，最常见的选择是利用模型本身，即pθ（xt | xt）-1，yt）=fθ（xt | xt-1). 然而，对于LGSS模型，可以利用高斯分布的性质得出最优方案，从而得到inoptθx（i）txa（i）tt-1，yt∝ gθ（yt）xt）fθxtxa（i）tt-1.= Nx（i）t；σhσ-2eyt+σ-2vφxa（i）tt-1i，σ, （22）带σ-2= σ-2v+σ-2e。这个粒子方案是最优的，因为它使当前时间步的增量粒子权重的方差最小化。

对于大多数其他人来说，这种无偏性对于PMH算法的有效性至关重要。因此，在PMH中不能使用仅在某些迭代中对粒子系统重新采样的自适应方法（Doucet and Johansen 2011，第3.5节）。然而，目前尚不清楚这究竟如何影响整个粒子系统，也就是说，这是否是全球最佳选择。12从PMH开始，对于非线性动力学模型的推理，最优方案难以解决，而使用状态转移模型。传播步骤由以下步骤实现：<传播粒子>=part1<-（sigmav^（-2）+sigmae^（-2））^（-1）part2<-sigmae^（-2）*y[t]part2<-part2+sigmav^（-2）*phi*粒子[new祖先，t-1]粒子[，t]<-part1*part2+rnorm（noParticles，0，sqrt（part1）），从（22）可以看出噪声方差之间的比率决定了提案的形状。本质上，有两种不同的极端情况（i）σe/σv 1和（ii）σe/σv 1.在第一个极端（i）中，提案的位置基本上由yt决定，标度主要由σe决定。这相当于基本上模拟gθ（yt | xt）中的粒子。当σeis较小时，通常允许仅使用少量颗粒运行颗粒过滤器。在第二个极端（ii）中，该方案基本上模拟了fθ（xt | xt）-1） 并且没有考虑到观察结果。总之，最优方案的性能和特点因此与模型的噪声方差及其相关大小有关<权重粒子>和<normalizeWeights>在这些步骤中，计算采样步骤所需的权重。可以使用不同的所谓权重函数计算权重，标准选择是观测模型，即vt（xt）=gθ（yt | xt）。

然而，对于LGSS模型，我们可以通过再次应用高斯分布的特性来获得最优加权函数，从而得到v（i）t，p（yt+1 | x（i）t）=Zgθyt+1xt+1fθxt+1x（i）tdxt+1=Nyt+1；φx（i）t，σv+σe.请记住，在这个实现中，新的观测yt+1首次在传播和加权步骤之间的粒子滤波器中引入。这就是为什么时间指数的选择有点复杂的原因。然而，在这一步中，yt+1中包含的信息是有益的，因为它提高了过滤器的准确性。也就是说，我们保留了传播后可能导致观测yt+1的粒子。在许多情况下，产生的权重很小，因此使用移位的原木权重可以避免数值精度问题。这是通过应用转换ev（i）t=log v（i）t来实现的- vmax，对于i=1，N、 其中，vmax表示logv（1:N）t中的最大元素。然后通过w（i）t=exp对权重进行归一化（确保它们的总和为1）ev（i）tPNj=1expev（j）t=经验(-vmax）exp对数v（i）t经验(-vmax）PNj=1exp对数v（j）t=v（i）tPNj=1v（j）t，（23）其中-VMAX取消，不影响砝码的相对大小。因此，（23）中的第一个和第三个表达式是等效的，但第一个表达式具有更好的数值精度。权重的计算由：Johan Dahlin，Thomas B。

Sch"on 13<weightParticles>=yhatMean<-phi*粒子[t]yhatVariance<-sqrt（sigmae^2+sigmav^2）权重[，t]<-dnorm（y[t+1]，yhatMean，yhatVariance，log=TRUE）<NormalizeWeghts>=maxWeight<-max（weights[，t]）weights[，t]-exp（weights[，MaxWeights）sumWeights<-sumWeights<-sum（weights[，sumWeights[，t]）和（weights[，t]）weights[，t]，-weights[，t]-weights[，t]，/Sumwei粒子[，t]由于R中的索引数组的约定，对应于y和x（1:N）t-1分别在模型中。这也是评论particleFilter中for循环的原因。现在我们看到权重取决于下一次观察，这就是为什么循环需要运行toT- 1而不是T<EstimateLibability>可通过将未规范化权重插入（11）中来估计似然p（θ| y1:t）。然而，由于可能性通常很小，因此最好是估算对数可能性，以获得更好的数值精度。这导致重写（11）以获得递归logbpNθ（y1:t）=logbpNθ（y1:t）-1） +（vmax+log“NXi=1expv（i）t#- 对数N），其中使用对数偏移的权重来提高数值精度。请注意，在估计对数可能性时，需要补偿Vmax的移位。