总结
针对一类部分可观测的多元扩散,我们引入了一种新的粒子滤波方案。我们考虑了各种观测方案,包括误差观测的扩散、多元扩散分量子集的观测以及强度为扩散(Cox过程)已知函数的泊松过程的到达时间。与当前可用的方法不同,我们的粒子过滤器不需要使用时间离散来近似过渡和/或观测密度。相反,他们建立在最近的方法论基础上,用于精确模拟扩散过程和无偏估计过渡密度。我们引入了广义泊松估计,它推广了Beskos及其同事的泊松估计。给出了粒子滤波方案的中心极限定理。
1.简介
在许多不同的科学学科中,人们对使用扩散过程来模拟连续时间现象非常感兴趣。这些过程可用于直接建模观测数据和/或描述分层模型中未观测到的过程。本文着重于估计给定部分信息的扩散路径。我们开发了一种新的粒子滤波器,用于分析在一组离散时间点部分观测到的一类多元扩散。
粒子滤波方法是分析部分观测到的离散时间动态模型的标准蒙特卡罗方法(Doucet等。,2001). 它们包括通过加权粒子群估计感兴趣的过滤密度。近似误差随着粒子数的增加而减小,N个,增加。然而,扩散过程的滤波比离散时间马尔可夫模型的滤波困难得多,因为除了少数特殊情况外,扩散的转移密度在所有情况下都不可用。在许多情况下,即使是观测密度也是难以处理的。因此,无法常规应用粒子滤波算法中的标准传播-加权-重采样步骤。
为了避免这些复杂情况,建议采用基于扩散时间离散化的进一步近似方法(参见Crisan示例等。(1999)和德尔·道德等。(2001)). 每个粒子从一个观测时间到下一个观测时间的传播是通过将时间增量拆分为M(M)比如说,曲目和表演M(M)根据适当的高斯分布进行中间模拟。作为M(M)这种高斯近似收敛到真正的扩散动力学。在这个框架中,算法的计算成本是有序的M(M)×N个,并获得真实的滤波分布M(M)和N个增加。
我们的方法不依赖于时间离散化,而是建立在扩散模拟的精确算法(EA)的最新研究基础上(Beskos和Roberts,2005; 贝斯科斯等。,2006a、,2005a) 扩散跃迁密度的无偏估计(Beskos等。,2006b、,2005b) 。该算法可以用多种方式来避免滤波问题中的时间离散化。Beskos的讨论中提出了EA在过滤问题中的潜力等。(2006b) ;请参阅肖邦、昆施,尤其是罗塞特和多塞特的贡献,他们也建议在这种情况下使用随机加权粒子滤波器。
一种可能性是在Gordon的引导粒子过滤器的实现中使用EA来传播粒子等。(1993)从而完全避免了M(M)每对观测时间之间的中间近似模拟。我们称之为精确传播粒子滤波器(EPPF)。在可能的情况下,更好的方法是通过使用拒绝采样使EA直接从(粒子近似)过滤密度进行模拟;我们称之为精确模拟粒子滤波器(ESPF)。
然而,我们偏爱的方法走向了不同的方向。我们在皮特和谢泼德的辅助粒子过滤器框架内工作(1999),其中粒子根据用户特定的密度从每个观测时间传播到下一个观测时间,然后适当加权以提供新滤波分布的一致估计值。由于过渡密度不可用,因此与每个粒子关联的权重很难处理。然而,我们的方法是给每个粒子分配一个随机的正权重,它是真实权重的无偏估计量。我们称之为随机加权粒子滤波器(RWPF)。我们的算法得到了滤波分布的一致估计。在比本文所考虑的更一般的情况下,用正无偏估计量替换粒子滤波器中的权重是一种有趣的可能性。确实,在第节中3.2我们表明,这种方法相当于用辅助变量方便地增加状态。
权的无偏估计的构造是本文的主要贡献之一,它具有独立的意义。这是基于泊松估计量Beskos的(PE)等。(2006b) ,我们称之为广义泊松估计(GPE)。该估计器保证返回正估计值(与PE不同),其效率(在方差和计算成本方面)可以比PE高出几个数量级。PE和GPE的最佳实现从理论上和通过仿真进行了深入研究。
我们引入的所有三个无时间离散化粒子过滤器都易于实现,RWPF是最容易和最灵活的,可以适应比这里考虑的更通用的上下文。仿真研究表明,RWPF比ESPF效率高得多,ESPF比EPPF效率更高。我们还提供了一个理论结果,该结果表明,与时间离散化方法相比,我们的滤波器可以具有显著的计算优势。我们建立了一个中心极限定理(CLT),用于使用EPPF、ESPF或RWPF估计滤波分布的期望值。这是肖邦成果的延伸(2004). CLT表明,对于固定的计算成本K(K),滤波分布的粒子近似误差随着K(K)−1/2在我们的方法中,虽然已知速率为K(K)−1/3或时间离散化方法较慢。
这里提出的方法的主要局限性是要求指定潜在扩散过程的随机微分方程(SDE)可以转换为具有正交扩散矩阵和梯度漂移的随机微分方程式。尽管该框架排除了一些重要类型的模型(如随机波动率模型),但它包含了一系列可以成功模拟许多物理过程的过程。然而,我们的方法可以处理各种离散时间观测方案。本文考虑了三种方案:扩散过程的噪声观测、多元扩散分量子集的观测以及强度随机且由已知扩散函数给出的泊松过程的到达时间。
论文组织如下。章节2介绍了潜在扩散模型和必要的符号,我们考虑的观测方案,以及我们测试所提方法的模拟数据集。章节三介绍了RWPF并说明了CLT。章节4介绍了构建RWPF所需的主要工具:GPE。建立了GPE的几个理论结果,并进行了仿真研究以评估其性能。章节5致力于对我们介绍的各种粒子滤波器的性能进行实证研究。还讨论了几个实施问题。章节6最后讨论了该方法的扩展,附录中包含了技术结果和证明。
2.信号、数据和假设
让信号由d日-维扩散过程
1
我们在整篇论文中假设漂移是已知的。我们的方法需要一些假设,我们在本段中总结了这些假设:
在这三个条件中(a)和(c)是弱的,而最严格的是(b),在遍历情况下对应于X(X)是时间可逆扩散。
过程(1)的跃迁密度通常很难处理,但有一个有用的表达式(参见Beskos示例等。(2006b) 以及达库纳-卡斯特尔和佛罗伦萨-兹米罗(1986))以下为:
2
在这个表达式中𝒩t吨(u个)表示d日-带平均值的维正态分布0和方差t吨我d日评估时间:u个∈R(右)d日,并对布朗桥进行期望,周秒,秒∈ [0,t吨],使用周0=x个0和周t吨=x个t吨。通常无法评估此公式中的期望值。
数据由部分观测值组成年1,年2,…,年n个,在离散时间点0t吨1<t吨2<…<t吨n个。我们考虑了三种可能的观察制度。
- (i)
观察到有误差的扩散:观察结果年我与当时的信号有关t吨我通过已知的密度函数(f)(年我|x个t吨我). 该模型扩展了一般状态空间模型,允许信号在时间上不断演化。有一个广泛的应用程序适合于这个框架;参见Doucet等。(2001)供参考。
- (ii)
部分信息:时间t吨我我们观察到年我=ζ(X(X)t吨我)对于某些不可逆的已知函数ζ(·). 例如,我们可以观察到d日-维度扩散。在此型号中(f)(年我|x个t吨我)=1代表所有x个t吨我对于其中ζ(x个t吨我)=年我.
- (iii)
斯过程:在这种情况下,数据由观测时间组成t吨我它们是随机的,被假定为泊松速率过程的到达ν(X(X)秒),对于某些已知函数ν这样的模型在保险中很流行(Dassios和Jang,2005)和金融(恩格尔,2000; 达菲和辛格尔顿,1999)它们最近被用来分析单分子实验的数据(寇等。,2005). 这一观察制度与前两个制度之间存在着重大差异。为了使符号与(i)和(ii)一致,我们让年我=t吨我表示我观察并定义可能性(f)(年我|x个t吨我−1,x个t吨我)为下一次观测后的概率密度t吨我−1时间为t吨我此密度可通过积分获得三
关于(X(X)秒,秒∈ (t吨我−1,t吨我))有条件地X(X)t吨我−1=x个t吨我−1和X(X)t吨我=x个t吨我这个条件过程的分布具有已知的关于布朗桥测度的密度,它在Beskos引理1中给出等。(2006b) 。因此,我们可以证明感兴趣的密度是4
其中期望是关于布朗桥定律的x个t吨我−1到x个t吨我.
我们采用贝叶斯方法并假设X(X)0我们的兴趣在于在线计算滤波密度:信号在时间上的后验密度t吨我根据最新观察结果t吨我,每个1我n个虽然这些密度很难处理,但我们提出了一种粒子滤波方案来递归估计每个观测时间点的这些密度。正如我们在第节中指出的那样6,我们的方法允许估计连续时间路径的滤波分布(X(X)秒,t吨我−1<秒<t吨我).
信号的一个更灵活的模型是扩散过程Z轴它解出了一个比我们在过程(1)中假设的方程更一般的SDE:
5
与过程(1)相比,过程(5)允许扩散系数依赖于状态。我们的方法直接适用于所有此类过程,前提是存在显式转换Z轴秒↦η(Z轴秒)=:X(X)秒,其中X(X)求解类型(1)的SDE;隐含漂移α可以很容易地用b条和Σ它必须满足我们已经规定的条件。在模型(i)中,可能性变为(f){年我|η−1(X(X)t吨我)},在模型(ii)中,数据为年我=ζ{η−1(X(X)t吨我)}在模型(iii)中,泊松强度为ν{η−1(X(X)秒)},其中η−1表示逆变换。因此,当d日=1; 在温和条件下,过程(5)可以通过以下方式转化为过程(1),因为一些武断u个*在扩散的状态空间中。此外,转换过程的漂移通常会满足我们指定的三个条件。然而,在更高的维度中,扩展更加困难。当d日>这可能是难以解决的,甚至是不可能的(阿伊特萨赫勒,2004). 即使这种转换是明确的,也可能意味着X(X)违反条件(b)。然而,许多物理系统可以用扩散成功地建模,扩散可以转化为过程(1)。
我们的粒子滤波方法将在两组模拟数据上进行说明,如下所示。
2.1. 示例1:观测到的正弦扩散有误差
信号满足
6
并且数据由噪声观测组成,年我∼N个(X(X)t吨我,σ2). 图。1(a) 在第节中5.1显示了此模型的模拟σ=0.2. 在这种情况下
7
这个过程与圆上的布朗运动密切相关。作为一个示例,它很方便,因为通过使用最基本的EA形式(Beskos中的EA1),可以很容易地从该过程中模拟离散时间骨架等。(2006a) ;R代码可根据作者的要求提供)。
2.2. 示例2:Ornstein–Uhlenbeck驱动的Cox过程
第二个数据集由泊松过程的到达时间组成,年我=t吨我,其强度由ν(X(X)秒),秒0,其中
和X(X)是Ornstein–Uhlenbeck(OU)流程,
OU过程是平稳的,具有高斯边缘分布,N个(0,1/2ρ). 因此,该模型的解释是X(X)增加泊松强度,而一对应于强度,当X(X)处于平均水平。示例数据集如图所示。三在第节中5.2,我们已经采取了虽然OU过程的转变密度是众所周知的,
观测密度(f)(年我+1|x个t吨我,x个t吨我+1)很难对付。
实施例1和实施例2分别是观察制度(i)和(iii)的实施例。我们将证明观测制度(ii)可以以与制度(i)类似的方式处理,因此我们没有包括附带的示例。
3.随机加权粒子滤波器
如第节所示2我们将表示观测时间t吨我通过年我、和第页t吨(·|·)表示系统随时间的跃迁密度t吨(见方程式(2)). 我们将写Δ我=t吨我+1−t吨我和过滤密度第页(x个t吨我|年1:我)将用表示π我(x个t吨我),其中根据标准惯例年1:我=(年1,…,年我). 为了简化符号,当我们在下面引入加权粒子时,我们将粒子和权重都下标为我而不是t吨我.
我们的目的是计算过滤密度π我(x个t吨我)递归地。基本概率计算得出这些密度的以下标准过滤递归:
8
粒子过滤器近似值π我(x个t吨我)由离散分布表示,其支持是一组N个颗粒,,具有关联的概率权重.替换对于π我(x个t吨我)式(8)中给出了(连续密度)近似值π我+1(x个t吨我+1),
9
粒子滤波算法的一次迭代的目的是构建进一步的粒子(离散分布)近似.
我们可以通过重要性抽样获得这样的粒子近似,而皮特和谢泼德的辅助粒子滤波器给出了实现这一点的一般框架(1999). 我们选择表格的提案密度
10
选择合适的提案,即选择s和q个,将在第节中对我们的具体应用进行分析5.
在某个时间模拟新粒子t吨我+1我们
- (a)
模拟粒子时间我,其中k个是一个离散随机变量的实现,其值为j个∈ {1,2,…,N个}有可能、和
- (b)
在time中模拟新粒子t吨我+1从.
指定给这对粒子的权重与…成比例
11
这是重复的N个时间产生加权粒子集的时间t吨我+1,,它提供了一个重要的采样近似值π我+1(x个t吨我+1). 可以重新规范权重,但不会对方法或其准确性产生重大影响。可以改进步骤(a)中的独立抽样:参见Carpenter的分层抽样思想等。(1999). 由此产生的粒子过滤器具有良好的理论特性,包括一致性(Crisan,2001)和CLT用于估计后验矩(Del Moral和Miclo,2000; 肖邦,2004; 昆施,2005),作为N个→∞. 在与初始条件指数遗忘相关的条件下,粒子滤波误差稳定为n个→∞ (德尔·道德和吉奥内特,2001; 昆施,2005).
当信号X(X)扩散过程是指过渡密度由于方程中的期望项,对于大多数感兴趣的扩散来说,表达式(11)中出现的是难以处理的(2). 此外,对于观测模型(iii)(也适用于更一般的模型),似然项表达式(4)中给出的不能进行解析计算。
我们通过给每个新粒子分配一个随机权重来规避这些问题,该权重是平均值为表达式(11)的随机变量的实现。该随机变量的构造和模拟见第节4,它基于方程中过渡密度的特定表达式(2). 在比本文考虑的更一般的情况下,用正无偏估计量替换权重是一种有趣的可能性。事实上,在第节3.2我们表明,这种方法相当于用辅助变量方便地增加状态。
3.1. 重量模拟
在所有模型中,与该对相关联的重量等于
12
哪里小时我+1是已知函数,对于0<u个<t吨,
期望是针对d日-维布朗桥周,开始于u个从周u个=x个并按时完成t吨在周t吨=z(z).
3.1.1. 模型(i)和(ii)
对于型号(i)和(ii)
和.在模型类型(ii)中,提案分发应选择仅建议以下值x个t吨我+1这样的话ζ(x个t吨我+1)=年我+1; 然后(f)(年我+1|x个t吨我+1)=1.
3.1.2. 模型(iii)
表达式(2)、(4)和(11)的合成,给出
章节4显示如何为每对(x、 z(z))和时间(u个,t吨),使用u个<t吨,其他辅助的变量V(V)、和函数第页(v(v),x、 z(z),u个,t吨)0,具有辅助变量根据适当的条件分布进行模拟(·|x、 z(z),u个,t吨)、和第页易于评估。我们的方法在权重中取代了难处理的术语及其无偏估计第页.
3.1.3. 随机加权粒子滤波器
对于我=0,…,n个−1和j个=1,…,N个执行以下步骤。
第2步:计算,; 如果是ESS<C类,对于某些固定常数C类,模拟k个我,j个从,k个=1,…,N个、和设置; 否则设置k个我,j个=j个和.
步骤3:模拟从.
第4步:模拟.
请注意,此算法包含关于是否在步骤2中传播之前对粒子重新采样的决定,该决定基于.常数C类可以解释为可接受的最小有效样本量。(见刘和陈(1998)以了解基于这种条件的重采样的基本原理。)是否进行重采样将影响赋予新粒子集的权重,这由以下不同的值来解释在步骤2中。最理想的情况是,步骤2的重新采样将在N个样本,例如Carpenter的分层抽样方案等。(1999)或刘和陈的剩余采样(1998).
3.2. 通过增强状态的等效公式
在上一节中,我们描述了一个通用的序贯蒙特卡罗方案,其中滤波分布的重要性抽样近似中的精确权重被正无偏估计量取代。我们现在证明,该方案相当于将普通辅助粒子滤波器应用于具有更丰富潜在结构的模型。我们证明了模型类型(i)和(ii)的这种等价表示,因为对参数的明显修改建立了模型类型的等价性(iii)。
根据我们的建设,有条件地X(X)t吨我,X(X)t吨我+1,t吨我和t吨我+1,V(V)我+1独立于V(V)j个和X(X)t吨j个对于任何j个这与我或我+1.此外,它很容易遵循第页有条件地X(X)t吨我=x个,第页(v(v)我+1,x个,x个t吨我+1,t吨我,t吨我+1)是的概率密度函数(X(X)t吨我+1,V(V)我+1)关于产品测量Leb(dz(z))×(d)v(v)|x、 z(z),t吨我,t吨我+1),其中Leb表示Lebesgue度量。
现在考虑一个具有未观测状态的替代离散时间模型(Z轴我,V(V)我),我=1,…,n个,Z轴我∈R(右)d日,具有非齐次马尔可夫转移密度
(该密度相对于Leb×问克)和观测数据年我具有观测密度(f)(年我+1|z(z)我,z(z)我+1). 通过构造Z轴我在这个模型中π我(x个t吨我)即表达式(8)中的过滤密度。考虑一个应用于此模型的辅助粒子滤波器,其中我们以概率进行选择每个现有粒子,并根据提案生成新粒子
哪里q个与表达式(10)中的建议密度相同。在这个离散时间模型中,与每个粒子相关联的权重是可处理的,并由方程式给出(13). 因此,加权样本是精确的粒子近似π我+1(x个t吨我+1),RWPF等价于该离散时间模型上的辅助粒子滤波器,其潜在结构已被辅助变量增强V(V)我.
这种等效表示揭示了我们方法的许多方面。首先,很明显,将正无偏估计量的多个实现平均化是低效的每粒子。相反,只需为每对粒子模拟一个估计器实现,就可以更有效地生成更多粒子。
其次,说明了RWPF结合了自举和辅助粒子滤波的优点。虽然很容易从概率分布进行模拟(如第节所述4),很难推导出其密度(相对于适当的参考测量值)。自从V(V)我根据该测度传播s,避免了其计算。这是bootstrap过滤器的一个吸引人的特性,它可以传播粒子,而不需要分析系统跃迁密度。然而Z轴我s通过用户指定的密度完成,该密度将信息纳入数据中。
第三,它表明RWPF将与应用于离散时间模型的辅助粒子滤波器具有相似的理论特性。这在第节中进行了探讨3.3.
3.3. 理论性质
考虑某些函数的后验均值的估计ϕ当时的状态t吨我,[ϕ(x个t吨我)|年1:我]. 研究粒子滤波器有效性的一种自然方法是将算法的极限行为视为N个→∞。对于标准辅助颗粒过滤器,肖邦(2004)引入了一种CLT来估计这类期望。本CLT直接适用于EPPF和ESPF。
在附录F我们扩展了肖邦的结果(2004)并给出了一个关于RWPF中随机权重的进一步必要条件,在该条件下,CLT仍然成立。这个额外条件是条件2。估计量方差的表达式[ϕ(x个t吨我)|年1:我]用RWPF得到的结果与标准情况下(即当权重已知时)的表达式不同,因为权重的随机性导致了一个额外的项(见方程式(27)–(29)以及对中定理3的评述附录F更多详细信息)。观察到RWPF可以重新表示为增强状态的标准粒子滤波器,这有助于肖邦方法的快速适应(参见第节3.2).
这种CLT的一个重要结果是,估计误差[ϕ(x个t吨我)|年1:我]符合规定N个−1/2。我们考虑的扩散问题的先前过滤方法基于
参见Crisan的示例等。(1999). Del Moral给出了给出此类方案中错误顺序的结果等。(2001). 对于(i)和(ii)等模型,误差是有序的N个−1/2前提是中间时间步数M(M)每次观察之间都会以一定的速度增加N个1/2因此,对于固定计算成本K(K)∝明尼苏达州错误以一定速度减少K(K)−1/3对于像(iii)这样的模型,如果可能性取决于两次连续观测之间的状态路径,则误差减少的速度将较慢,例如。K(K)−1/4(德尔·莫拉尔等。,2001)或K(K)−1/6(Crisan的定理1.1等。(1999))。
4.广义泊松估计
我们已经激发了对正无偏估计量的模拟需求
14
期望是针对d日-维布朗桥周在本节中,我们介绍了一种推导此类估计量的方法,并提供了有关建议估计量方差的理论和模拟结果。这些结果在粒子滤波之外具有独立的意义,因此我们以通用的方式介绍我们的方法,其中克是一个任意函数,假设它只在上是连续的R(右)d日。我们假设周0=x个和周t吨=z(z),用于任意x、 z(z)∈R(右)d日和t吨>通过布朗桥的时间同质性,我们的方法扩展到积分极限变为u个和u个+t吨,对于任何u个>0.
贝斯科斯等。(2006b) 提出了期望(14)的无偏估计量,PE:
15
κ是具有平均值的泊松随机变量λt吨,的ψj个s均匀分布于[0,t吨]、和c(c)∈R(右)和λ>0是任意常数。(这里和下面我们假设空产品,即κ=0,取值1。)PE的两个主要缺点是它可能返回负估计值,并且它的方差不一定是有限的。当克有界。然而,在我们的背景下,这是一个非常有限制性的假设。因此,在这里,我们引入了一组无偏和正估计[E类](14)推广了PE。我们考虑的方法允许c(c)和λ依赖周和许可证κ具有一般离散分布。首先,我们需要能够模拟随机变量L(左)周和U型周具有
16
并且能够模拟周秒在任何时候秒给定不等式(16)所隐含的条件。对于无界克这是非平凡的。然而,由于在Beskos中引入了一种有效的算法,这两种模拟都变得可行等。(2005a) ●●●●。结构概述见附录A.
让U型周和L(左)周满足条件(16)和ψj个,j个1是[0,t吨]. 然后,估计器(14)可以重新表示为
17
哪里κ是具有条件概率的离散随机变量P(P)(κ=k个|U型周,L(左)周)=第页(k个|U型周,L(左)周). 上述论点中的第二个等式是通过使用支配收敛和Fubini定理(由和的正性所支持)获得的。
我们可以推导出[E类](14)通过指定第页(·|U型周,L(左)周). 所有此类估算器系列将被称为GPE:
18
以下定理(在附录B)为以下项目提供了最佳选择第页(·|U型周,L(左)周).
定理1。给定GPE的条件二阶矩U型周和L(左)周是
19
如果
20
然后通过选择将第二个力矩最小化
21
最小二阶矩由
22
虽然表达式(21)的右侧无法进行分析评估,但它可以指导选择合适的第页(·|U型周,L(左)周). 如果周已知,最优方案为平均泊松
23
我们将讨论表达式(23)用于选择好提案的两种可能方式。
保守的方法需要第页(·|U型周,L(左)周)为平均泊松(U型周−L(左)周)t吨(的上界λ周). 我们称此估计器为GPE-1。GPE-1的一个优点是它的二阶矩有界于[出口(-2L(左)周t吨)]. 因此,在温和明确的条件下克,包含在以下定理中(在附录C),保证估计量的方差是有限的。
定理2。GPE-1具有有限方差的一个充分条件是
自λ周是随机的,另一种方法是引入(外生)随机平均值,并假设第页(·|U型周,L(左)周)是具有此随机平均值的泊松。对于可控制性,我们选择随机平均值以获得伽马分布,当第页(·······································································U型周,L(左)周)变为负二项分布(估计量GPE-2):
24
哪里γ周和β分别表示负二项分布的平均值和离散参数。由于负二项分布的尾部比PE重,因此只要PE有有限方差,GPE-2就会有有限方差。然而,如果γ周选择为近似值[λ周|U型周,L(左)周]. 有各种各样的特别的可以提供此期望的粗略估计的方法。应用Jensen不等式交换方程中的平方幂积分(23),然后近似通过,建议服用
25
模拟研究(部分内容见第节4.1)表明这种选择在实践中非常有效,GPE-2的方差比PE或GPE-1小几个数量级。积分通常很容易计算;否则,可以使用粗略的近似值。
我们的介绍仅限于表达式(14)中的期望是关于布朗桥测度的情况。然而,正如Beskos指出的那样等。(2006b) 当对任意扩散桥测度取期望值时,只要能从该测度模拟出精确的骨架,就可以用完全相同的方法构造PE。GPE也可以在这个更广泛的框架中实施,前提是周可以构造为满足条件(16)。
4.1. 模拟研究
我们考虑一个光滑有界测试函数={sin(u个)2+科斯(u个)+1}/2. 这是根据示例1选择的。功能克是周期性的,周期为2π.In[0,2π]它在0和2处有局部极小值π,全局最小值为π最大值为π/3和5π/3.自克以9/8为界,我们可以通过设置c(c)9/8. 在这种约束下,Beskos等。(2006b) 认为一个好的选择是c(c)=λ=9/8. 仿真实验表明,GPE-2的性能对色散参数的选择具有很强的鲁棒性β。我们在示例中将其修正为β=10.表格1总结了基于10的估计量方差估计4模拟值。我们看到GPE-2可以比PE更有效,特别是考虑到[κ]. 一般来说,PE的性能对c(c)和λGPE-1的效率通常低于GPE-2。表1也给出了var的值(E类)它的值比PE、GPE-1或GPE-2中的任何一个都要小得多(几个数量级),说明了这些辅助变量构造的绝对效率成本。
. | 估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
方差 | 体育课 | 0.202 | 0.200 | 0.027 |
| GPE-1级 | 4.21×10−3 | 0.208 | 0.034 |
| GPE-2级 | 2.08×10−3 | 0.220 | 0.033 |
| 无功功率,无功功率(E类) | 3.74×10−5 | 3.27×10−3 | 4.72×10−3 |
[κ] | 体育课 | 1.118 | 1.126 | 1.121 |
| GPE-1级 | 0.130 | 1.091 | 0.744 |
| GPE-2级 | 0.119 | 0.329 | 0.735 |
. | 估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
方差 | 体育课 | 0.202 | 0.200 | 0.027 |
| GPE-1级 | 4.21×10−3 | 0.208 | 0.034 |
| GPE-2级 | 2.08×10−3 | 0.220 | 0.033 |
| 无功功率,无功功率(E类) | 3.74×10−5 | 3.27×10−3 | 4.72×10−3 |
[κ] | 体育课 | 1.118 | 1.126 | 1.121 |
| GPE-1级 | 0.130 | 1.091 | 0.744 |
| GPE-2级 | 0.119 | 0.329 | 0.735 |
. | 估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
方差 | 体育课 | 0.202 | 0.200 | 0.027 |
| GPE-1级 | 4.21×10−3 | 0.208 | 0.034 |
| GPE-2级 | 2.08×10−3 | 0.220 | 0.033 |
| 无功功率,无功功率(E类) | 3.74×10−5 | 3.27×10−3 | 4.72×10−3 |
[κ] | 体育课 | 1.118 | 1.126 | 1.121 |
| GPE-1级 | 0.130 | 1.091 | 0.744 |
| GPE-2级 | 0.119 | 0.329 | 0.735 |
. | 估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
方差 | 体育课 | 0.202 | 0.200 | 0.027 |
| GPE-1级 | 4.21×10−3 | 0.208 | 0.034 |
| GPE-2级 | 2.08×10−3 | 0.220 | 0.033 |
| 无功功率,无功功率(E类) | 3.74×10−5 | 3.27×10−3 | 4.72×10−3 |
[κ] | 体育课 | 1.118 | 1.126 | 1.121 |
| GPE-1型 | 0.130 | 1.091 | 0.744 |
| GPE-2级 | 0.119 | 0.329 | 0.735 |
我们还研究了PE和GPE-2的效率如何随时间增量而变化t吨尤其是对于小型t吨(未显示结果)。这些经验结果表明,PE和GPE-2的误差变异系数为O(运行)(t吨δ)对一些人来说δ>0,但其值为δ这两个估计值不同。在我们调查的案例中,GPE-2似乎比PE具有更快的收敛速度。
此模拟研究的结果已被其他功能验证克(结果未显示)。我们已经对可微性进行了实验(例如。)和不可微(例如。=|u个|)无限函数。在这些情况下,不可能设计一个概率为1且返回正估计值的PE。我们再次发现,GPE-2的表现明显优于PE。
值得一提的是,存在其他蒙特卡罗方法,这些方法产生一致但有偏差的表达式估计值(14)。通过将表达式(14)中的时间积分替换为基于数字的黎曼近似,可以得到一个这样的估计量,M(M)比如说,中间点。该技术用于构建Nicolau的转移密度估计器(2002)并有效地支持Durham和Gallant的转移密度估计(2002)(当扩散过程具有恒定的扩散系数时)。达勒姆和格兰特的做法(2002)已用于马尔可夫链蒙特卡罗和滤波应用(Golightly和Wilkinson,2006; Chib公司等。,2006; 离子,2003). 在过滤上下文中,它提供了一种替代RWPF的方法,其中权重是近似值。本文的目的不是仔细比较RWPF与此类变体。然而,作为一个例子,我们在估算跃迁密度的背景下进行了一个非常小的比较,第页t吨(z(z)|x个),表达式(6)的t吨=1和x个和z(z)如表所示1。我们比较了四种方法。其中两个基于方程式(2)并使用PE和GPE-2生成期望的估计值。另外两个DG-1和DG-5是Durham和Gallant的两个实现(2002)估计器,分别有一个和五个中间点。我们根据它们的均方根误差除以真实值(即变异系数)来比较这些方法。作为真实值,我们使用了GPE-2的估计值。比较结果见表2注意,DG-1和DG-5模拟的变量比GPE-2多得多,以构建它们的估计。
估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
体育课 | 1.25 | 0.93 | 0.17 |
GPE-2级 | 0.13 | 0.78 | 0.2 |
DG-1公司 | 0.5 | 0.45 | 0.3 |
DG-5型 | 0.28 | 0.19 | 0.22 |
估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
体育课 | 1.25 | 0.93 | 0.17 |
GPE-2级 | 0.13 | 0.78 | 0.2 |
DG-1型 | 0.5 | 0.45 | 0.3 |
DG-5型 | 0.28 | 0.19 | 0.22 |
估算员. | 以下(x,z)对的结果:. |
---|
x个=0时,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
体育课 | 1.25 | 0.93 | 0.17 |
GPE-2型 | 0.13 | 0.78 | 0.2 |
DG-1公司 | 0.5 | 0.45 | 0.3 |
DG-5型 | 0.28 | 0.19 | 0.22 |
估算员. | 以下(x,z)对的结果:. |
---|
x个=0,z(z)=0. | x个=0,z(z)=π. | x个=π,z(z)=π. |
---|
体育课 | 1.25 | 0.93 | 0.17 |
GPE-2级 | 0.13 | 0.78 | 0.2 |
DG-1公司 | 0.5 | 0.45 | 0.3 |
DG-5型 | 0.28 | 0.19 | 0.22 |
5.粒子滤波器对模拟数据的比较
现在,我们在第节中介绍的两个示例中演示了各种粒子过滤器的性能2.
5.1. 正弦扩散分析
我们首先考虑分析示例1的正弦扩散。模拟数据如图所示。1(a) ●●●●。我们比较了粒子滤波器的四种实现,每种实现都使用基于EA的方法来模拟扩散,从而避免了时间离散化:
- (a)
EPPF,使用EA实现引导过滤器;
- (b)
ESPF,它使EA通过过滤密度的拒绝采样进行模拟;
- (c)
RWPF1,使用PE实现RWPF(见表1)模拟重量;
- (d)
RWPF2是使用GPE-2模拟权重的RWPF实现。
EPPF和ESPF的实施细节见附录D.
在这个简单的例子中,ESPF比EPPF更有效,因为它具有相同的计算成本,但它是从最优方案分布中提出的。然而,我们已经有效地实现了ESPF,利用了这个简单模型的几个细节,特别是高斯似然和漂移有界的事实。在更一般的模型中,ESPF的实施可能要困难得多,并且由于接受概率较小,与EPPF的比较也不太有利。
在这种情况下,其中φ是有界的,我们可以通过替换来加快GPE-2的实现,而实际上不会损失效率U型周在表达式(24)和(25)中乘以9/8,这是φ在这种情况下,无需模拟U型周和L(左)周我们在RWPF2中实现了这种简化。
算法EPPF、RWPF1和RWPF2使用Carpenter的分层重采样算法等。(1999),每次迭代都会重新采样。对于RWPF1和RWPF2,我们根据最佳建议分布为新粒子选择建议分布,如果正弦扩散由Ozaki离散方案近似,则获得最佳建议分布(详细信息见附录E). 对于EPPF,我们选择了s为从该近似值中获得的值。
设置每个算法中使用的粒子数,以便每个过滤器具有可比较的中央处理器单元(CPU)成本,从而为每个算法使用500、500、910和1000个粒子。对于这些数量的粒子,EPPF和ESPF平均需要1360个粒子的提案,并且需要在算法每次迭代的接受-拒绝步骤(c)中进行675个布朗桥模拟。通过比较,RWPF1和RWPF2分别模拟了910和1000个粒子,平均需要1025和850个布朗桥模拟来生成每次迭代的随机权重。
请注意,四种算法的比较CPU成本,尤其是EPPF和ESPF与RWPF1和RWPF2的比较,取决于基础扩散路径。EPPF和ESPF内的接受概率取决于和x个t吨我+1,当这两个值都接近0(mod(2π)). (从长远来看,扩散将不经常访问这些区域,并将在那里停留很短的时间。)因此,在状态空间的这个区域花费更多或更少时间的模拟路径将导致EPPF和ESPF分别具有更大或更小的CPU成本。
我们比较了四个滤波器,基于每个滤波器500次独立运行期间状态滤波分布平均值估计值的可变性。结果如图所示。2,一次运行RWPF2的输出如图所示。1图中的比较结果。2用于估计每次迭代时滤波分布的平均值(对滤波分布的各个分位数获得了类似的结果)。他们显示RWPF2表现最好,平均效率比RWPF1提高15%,比ESPF提高50%,比EPPF提高200%。对这些结果的解释表明(例如)ESPF需要与N个=750(取该数据集CPU成本的1.5倍),以获得与RWPF2相当的精度。
改变模型的参数和算法的实现将影响算法的相对性能。尤其是增加或减少σ2测量误差的方差将分别增加或降低EPPF相对于其他滤波器的相对效率。类似的结果出现在Δ我分别减少或增加。其他三种算法的相对性能似乎对此类变化更为稳健。我们考虑实施EPPF,并使用正弦扩散的欧拉近似而非Ozaki近似来构建RWPF1和RWPF2的建议分布,但这些变化对方法的性能没有任何显著影响。我们还考虑减少重采样频率,设置C类=N个/RWPF算法第2步中的4(因此,当s小于N个/4) 这大大降低了算法的性能(RWPF1和RWPF2的系数为2)。
我们还研究了增加观测之间的时间量Δ的影响。为此,我们使用了上述数据
时间点。
为了测量这些不同场景下过滤器的性能,我们使用了Carpenter的有效样本大小ESS等。(1999). ESS是根据滤波器独立运行期间后验均值估计值的方差来计算的,但将此方差与后验方差进行比较,以确定从后验值中独立提取多少个参数可以产生相同精度水平的估计值。我们着重于估计观测时间状态的后验平均值;如果秒2是粒子滤波器估计值的样本方差[X(X)t吨我|年1:我]100次独立运行,以及是var的估计值(X(X)t吨我|年1:我),则ESS为注意,通过ESS比较滤波器相当于根据估计量的方差比较滤波器。
表三给出了各种Δ值的ESS值。我们看到,随着Δ的增加,ESS值急剧下降,并且当Δ=20时,滤波器效率低下。这种性能下降是由于在这种情况下随机权重的巨大可变性。这些权重的可变性是由于
过滤器. | 以下值的结果Δ:. |
---|
Δ=10. | Δ=20. |
---|
RWPF2型 | 73 | 5 |
离散化 | 80 | 12 |
伪RWPF2 | 923 | 933 |
过滤器. | 以下值的结果Δ:. |
---|
Δ=10. | Δ=20. |
---|
RWPF2型 | 73 | 5 |
离散化 | 80 | 12 |
伪RWPF2 | 923 | 933 |
过滤器. | 以下值的结果Δ:. |
---|
Δ=10. | Δ=20. |
---|
RWPF2型 | 73 | 5 |
离散化 | 80 | 12 |
伪RWPF2 | 923 | 933 |
过滤器. | 以下值的结果Δ:. |
---|
Δ=10. | Δ=20. |
---|
RWPF2型 | 73 | 5 |
离散化 | 80 | 12 |
伪RWPF2 | 923 | 933 |
为了评估由(a)引起的数量,我们尝试了一种粒子过滤器,该过滤器通过在一组离散时间点上模拟布朗桥来数值估计表达式(26)(在本例中,我们每半个时间单位采样一次值),然后使用这些值来数值计算积分。这种方法与达勒姆和加兰特的重要抽样方法密切相关(2002)和尼古拉(2002); 参见第节4.1表中也给出了该过滤器的结果三(请注意,ESS值忽略了通过此数值近似引入的任何偏差),我们再次看到较小的ESS值,尤其是对于Δ=20。该滤波器的性能与RWPF非常相似,这表明(b)中的蒙特卡罗可变性是导致这种情况下RWPF性能不佳的一个小原因。
最后,我们尝试在当前没有观测到的所有整数时间间隔引入伪观测。然后按上述方式运行RWPF,但在存在这些非信息性观察的时间点,对权重没有可能的贡献。想法是现在Δ=1,这样随机权重的方差表现良好,但我们仍然可以在观察之前调整单位时间间隔内的扩散路径,以考虑该观察中的信息。结果再次显示在表中三ESS值非常高(接近最佳值,即粒子数1000)。由于布朗桥模拟的总计算成本不变,因此,通过添加这些额外的伪观测,计算成本仅大致增加了一倍。这些结果对于选择引入这些无信息观察的频率是相当可靠的(结果未显示)。
5.2. 考克斯过程分析
我们现在考虑将RWPF应用于第节中的示例22,OU驱动的Cox过程。我们分析的数据如图所示。三(a) ●●●●。要使另外两个基于EA-的粒子过滤器(EPPF和ESPF)适应此问题,要么是不可能的,要么是困难的。例如,我们不能实现EPPF,因为似然函数是不可处理的。所以我们只关注RWPF在估计滤波分布时的效率|X(X)t吨|.
我们对RWPF的实施是基于从先前分布中提出粒子,所以和就是OU跃迁密度.我们用GPE-2模拟了随机权重。我们计算了每个观测时间以及所选的56个伪观测时间的滤波密度,因此计算滤波密度的两个连续时间之间的最大时间差为0.1。这是必要的,以避免需要布朗桥模拟来模拟权重对于长时间的互观测而言太大的次数,并控制随机权重的方差(见上文)。这些非观测时间的似然函数通过删除ν(x个t吨我)来自表达式(4)。
我们将粒子数设置为1000,并在s小于100(C类=N个/第节算法步骤2中的10三). 虽然正弦扩散的结果表明,这将导致算法重采样次数太少,但我们选择了一个较低的阈值,以便通过粒子过滤器权重的ESS如何随时间衰减来监测粒子过滤器的性能。该滤波器一次运行的结果如图所示。三(a) ●●●●。该方法的计算效率可以通过图来衡量。三(b) 其中的ESSs随时间绘制。
6.讨论
我们描述了如何在粒子滤波器中使用最新的方法来精确模拟扩散和无偏估计扩散指数泛函,从而避免了时间离散化的需要。在我们介绍的方法中,特别关注了RWPF,它实现了一个辅助粒子过滤器,但模拟了分配给每个粒子的权重。我们表明,这种方法相当于应用于适当扩展模型的辅助粒子滤波器。我们预计,这种方法将在不同于本文所考虑的模型的模型中有有趣的应用,然而,这些模型涉及棘手的动力学或可能性。
我们专注于滤波问题,根据迄今为止的观察结果估计当前状态。然而,对预测的扩展是微不足道的,只需要从状态方程进行模拟的能力,这可以通过EA算法实现。使用北川的想法也很简单(1996),其中每个粒子存储其轨迹的历史,以获得平滑密度的近似值(给定到目前为止的观测值,过去某个时间的状态密度)。
尽管粒子只存储每个观测时间的状态值,但填充这些时间之间的扩散路径以在任何时候产生关于状态的推断都是简单的。分布的粒子近似(X(X)秒,t吨我−1<秒<t吨我)以数据为条件年1:我可以通过使用当前的一组加权粒子来构建带砝码,如下所示。首先我们需要介绍一些符号;我们用表示粒子在时间上的值t吨我−1从中j个第th个粒子t吨我下降。粒子近似由一组加权路径给出带砝码。每条路径都是一个扩散桥,从并于完成可以使用EA进行模拟,如Beskos所述等。(2006年a,2005a) ●●●●。在观测区域(i)和(ii)中,应用EA模拟具有密度的扩散桥,其布朗桥测度由下式给出,而在区域(iii)中,相应的密度为。可以直接利用此表示来推断X(X)在观察时间之间。
致谢
作者感谢联合编辑和裁判的宝贵意见。
参考文献
附录A:分层布朗运动
Beskos中提出的算法等。(2005a) 首先创建示例空间的分区周对于给定的周0=x个和周t吨=年.书写x个=(x个1,…,x个d日),对于用户特定的常量一>√(t吨/3) ,的子集序列R(右)d日形成为一个j个={u个=(u个1,…,u个d日):分钟(x个t吨我,年我)−青年成就组织<u个我 最大(x个t吨我,年我)+青年成就组织},j个0,其中j个 一个j个=R(右)d日。此序列定义窗体示例空间的分区,其中路径属于D类j个当且仅当路径超过由一个j个−1但不是由一个j个.在Beskos等。(2005a) 演示了如何模拟随机变量,该变量决定了D类j个秒周属于,以及如何模拟周在这个随机变量的任何时间集合上:分层布朗桥结构。自克假设是连续的,知道周∈D类j个足以确定U型周和L(左)周满足条件(16)。事实上,在简化的设置中克是有界的,正如示例1中的正弦扩散一样,可以避免分层布朗桥结构,因为它易于选择U型周和L(左)周独立于周.
附录B:定理1的证明
然后,表达式(19)建立如下:
上面使用了Fubini定理和支配收敛性(这是有效的,因为被积函数几乎肯定是正的)。通过使用以下结果(可以通过使用Jensen不等式容易地证明)来获得表达式(22)。让(f)我>0用于我=1,2,…. 然后是第页我s最小化在约束∑下第页我=1由下式给出第页我=√(f)我/Σ √(f)我.
附录C:定理证明2
GPE-1估计值小于或等于exp(−L(左)周t吨)因此,如果[经验(−L(左)周t吨)]<∞,其中期望是关于d日-维布朗桥x个时间为0到年时间t吨然而,
哪里M(M)我=支持0秒t吨|周我|利用定理2中的增长界。此外,
因此,仍需约束d日这个表达式右边的积分。然而,根据Bachelier-Levy公式计算布朗运动和桥梁的撞击时间,
等等
它后退得像w个−k个日志(w个)作为w个→∞, 从而得出结论。
附录D:精确传播粒子滤波器和精确模拟粒子滤波器,例如1
EPPF根据以下程序生成新粒子:
- (a)
选择当前粒子之一,其中粒子j个被概率选择;
- (b)
建议x个t吨我+1从平均值的正态分布和方差Δ我;
- (c)
以概率exp{−cos接受此建议(x个t吨我+1)−1}; 如果提案被拒绝,则返回(a);
- (d)
其中期望是关于布朗桥定律的和周Δ我=x个t吨我+1、和φ在方程式中给出(7). 如果提案被拒绝,则返回(a);否则x个t吨我+1是当时的新粒子t吨我+1带重量w个我+1=(f)(年我+1|x个我+1).
步骤(d)通过使用Beskos中所述的回顾性抽样进行等。(2006a) ●●●●。ESPF如上所述进行,但步骤(a)和(b)替换为步骤
重复该算法,直到N个的值x个t吨我+1已验收,每个重量为1/N个.
附录E:示例1的提案分发
考虑满足SDE(1)的扩散,其中d日=1表示简单。该SDE的Ozaki近似是基于某个值附近漂移的一阶泰勒级数展开x个对于例1的正弦扩散,我们得到了近似的SDE
所以是示例2中定义的OU进程ρ=cos(x个)和σ=1.计算我们计算了由Ozaki近似给出的关于和似然函数(f)(年我+1|x个t吨我+1). 定义
和
我们得到了平均值正常(ησ2+年我+1τ2)/τ2σ2和方差η2τ2/(η2+τ2). 此外,我们计算.
附录F:中心极限定理
为了简化符号,我们考虑了粒子滤波器的一个特例,该特例与肖邦的粒子滤波器相似(2004). 我们为时间选择提案密度t吨我+1拥有; 我们假设在步骤2中。肖邦的粒子滤波器(2004)每次拆分模拟粒子t吨我+1进入之内
- (a)
粒子的时间重采样t吨我和
- (b)
这些粒子随时间的传播t吨我+1.
我们的独立同分布抽样假设等价于肖邦的多项式重抽样情况(2004). (如果刘和陈的剩余采样方法相同,则CLT的条件相同(1998),但差异不同。)为了简单起见,我们考虑观察模型(i)或(ii),尽管结果很容易扩展到观察模型(iii)。
让,其中k个我,j个是在模拟j个-时间粒子t吨我和是j个第th个粒子t吨我和粒子一起t吨我−1它是从那里传下来的。还让θ我表示给定的条件期望θ我类似地,让μ我(θ我)=(x个我−1,x个我,t吨我−1,t吨我),并表示为R(右)我无偏估计量μ我(θ我),即。[R(右)我]=μ我(θ我). 一个重要的数量是.
我们定义我[ϕ]和var我(ϕ)为任意函数的后验均值和方差ϕ(θ)时间我,并考虑粒子滤波估计我[ϕ]. 让是密度第页(x个t吨我−1|年1:我−1)q个(x个t吨我|x个t吨我−1). 最后定义q个我[ϕ]和varq个我(ϕ)是的条件期望和方差的简写ϕ(θ我)关于q个(x个t吨我|x个t吨我−1)(这些是的功能x个t吨我−1). 我们将‖·‖表示为欧几里德范数,并递归地定义Φ我是一组可测量的函数ϕ这样,对一些人来说δ>0,,并且该函数x个t吨我−1↦q个我[小时我μ我ϕ]为Φ我−1.
定理3。考虑一个函数ϕ; 定义通过归纳,
27
28
29
那么如果,为了所有人我,
对于任何ϕ∈ Φ我,我[ϕ]和V(V)我(ϕ)是有限的,并且我们在分布上收敛为粒子数,N个,趋于∞:
备注1。方程式(27)–(29)指迭代时由于传播、加权和重采样阶段而导致加权粒子方差的变化我.仅方程式(28)与肖邦的各自成绩不同(2004)这是由于右侧的第二项,它表示由于权重的随机性而导致的方差增加。条件(a)取自肖邦(2004)并适用于标准颗粒过滤器;条件(b)和(c)是新的,是限制随机权重方差的条件,确保V(V)我(ϕ)是有限的。
证明。我们改编了肖邦的感应证明(2004),依次考虑传播、加权和重采样步骤。我们的滤波器仅在加权步长方面与标准粒子滤波器不同,因此我们只需要调整肖邦引理A2的结果(2004). 事实上,方程式(27)和(29)与肖邦的对应量相同(2004); 因此,仍需显示等式(28). 我们定义常数和ϕ*=R(右)我小时我(ϕ−E类我[ϕ])/K(K)在扩大的信号空间框架内,我们可以应用方程式(4)肖邦的(2004),给予
现在我们可以计算q个我[ϕ*]首先对辅助变量取期望值(条件是θ我). 这给了q个我[ϕ*]=q个我[μ国际卫生组织我(ϕ−E类我[ϕ])/K(K)]. 同样,我们获得
30
31
(这里是方程中的期望和方差(30)与辅助变量有关。)将这些结果结合起来得出方程式(28). 正则性条件(a)-(c)也直接转换。
©2008皇家统计学会