总结

针对一类部分可观测的多元扩散,我们引入了一种新的粒子滤波方案。我们考虑了各种观测方案,包括误差观测的扩散、多元扩散分量子集的观测以及强度为扩散(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日-维扩散过程

d日X(X)=α(X(X))d日+d日B类,[0,t吨].
1

我们在整篇论文中假设漂移是已知的。我们的方法需要一些假设,我们在本段中总结了这些假设:

  • (a)

    α在其所有自变量中都是连续可微的,

  • (b)

    有一个函数一个:R(右)d日R(右)这样的话α(u个)=∇一个(u个)和

  • (c)

    有一个>−∞这样φ(u个):={‖α(u个)‖2+∇2一个(u个)}/2−0

在这三个条件中(a)和(c)是弱的,而最严格的是(b),在遍历情况下对应于X(X)是时间可逆扩散。

过程(1)的跃迁密度通常很难处理,但有一个有用的表达式(参见Beskos示例等。(2006b) 以及达库纳-卡斯特尔和佛罗伦萨-兹米罗(1986))以下为:

第页t吨(x个t吨x个0)=N个t吨(x个t吨x个0)经验{一个(x个t吨)一个(x个0)t吨}E类[经验{0t吨ϕ()d日}].
2

在这个表达式中𝒩t吨(u个)表示d日-带平均值的维正态分布0和方差t吨d日评估时间:u个R(右)d日,并对布朗桥进行期望,,∈ [0,t吨],使用0=x个0t吨=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吨1t吨ν(X(X))d日},
    关于(X(X),∈ (t吨−1,t吨))有条件地X(X)t吨−1=x个t吨−1X(X)t吨=x个t吨这个条件过程的分布具有已知的关于布朗桥测度的密度,它在Beskos引理1中给出等。(2006b) 。因此,我们可以证明感兴趣的密度是
    ν(x个t吨)N个t吨t吨1(x个t吨x个t吨1)第页t吨t吨1(x个t吨x个t吨1)经验{一个(x个t吨)一个(x个t吨1)}E类[经验[t吨1t吨{ϕ()+ν()}d日]],
    4
    其中期望是关于布朗桥定律的x个t吨−1x个t吨.

我们采用贝叶斯方法并假设X(X)0我们的兴趣在于在线计算滤波密度:信号在时间上的后验密度t吨根据最新观察结果t吨,每个1n个虽然这些密度很难处理,但我们提出了一种粒子滤波方案来递归估计每个观测时间点的这些密度。正如我们在第节中指出的那样6,我们的方法允许估计连续时间路径的滤波分布(X(X),t吨−1<<t吨).

信号的一个更灵活的模型是扩散过程Z轴它解出了一个比我们在过程(1)中假设的方程更一般的SDE:

d日Z轴=b条(Z轴)d日+Σ(Z轴)d日B类,[0,t吨].
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)η(Z轴)=u个*Z轴(z(z))1d日z(z),因为一些武断u个*在扩散的状态空间中。此外,转换过程的漂移通常会满足我们指定的三个条件。然而,在更高的维度中,扩展更加困难。d日>这可能是难以解决的,甚至是不可能的(阿伊特萨赫勒,2004). 即使这种转换是明确的,也可能意味着X(X)违反条件(b)。然而,许多物理系统可以用扩散成功地建模,扩散可以转化为过程(1)。

我们的粒子滤波方法将在两组模拟数据上进行说明,如下所示。

2.1. 示例1:观测到的正弦扩散有误差

信号满足

d日X(X)=(X(X))d日+d日B类,
6

并且数据由噪声观测组成,N个(X(X)t吨,σ2). 图。1(a) 在第节中5.1显示了此模型的模拟σ=0.2. 在这种情况下

ϕ(u个)={(u个)2+科斯(u个)+1}/2
7
图1

(a) 正弦扩散的实现(图解的)于[0100](图解的,单位时间间隔100个观测值;图解的,RWPF2获得的观测时间扩散的滤波分布平均值N个=1000个颗粒)和(b)观测数据与过滤平均值之间的差异(图解的)和90%可信区间(图解的)来自RWPF2(尽管为了清晰起见,它们始终显示,但可信区间仅在观测时间计算)

这个过程与圆上的布朗运动密切相关。作为一个示例,它很方便,因为通过使用最基本的EA形式(Beskos中的EA1),可以很容易地从该过程中模拟离散时间骨架等。(2006a) ;R代码可根据作者的要求提供)。

2.2. 示例2:Ornstein–Uhlenbeck驱动的Cox过程

第二个数据集由泊松过程的到达时间组成,=t吨,其强度由ν(X(X)),0,其中

ν(x个)=+β|x个|,

X(X)是Ornstein–Uhlenbeck(OU)流程,

d日X(X)=ρX(X) d日+d日B类.

OU过程是平稳的,具有高斯边缘分布,N个(0,1/2ρ). 因此,该模型的解释是X(X)增加泊松强度,而对应于强度,当X(X)处于平均水平。示例数据集如图所示。在第节中5.2,我们已经采取了=0,β=20ρ=12虽然OU过程的转变密度是众所周知的,

X(X)t吨X(X)0=x个0~N个[经验(ρt吨)x个0,12ρ{1经验(2ρt吨)}],

观测密度(f)(+1|x个t吨,x个t吨+1)很难对付。

实施例1和实施例2分别是观察制度(i)和(iii)的实施例。我们将证明观测制度(ii)可以以与制度(i)类似的方式处理,因此我们没有包括附带的示例。

3.随机加权粒子滤波器

如第节所示2我们将表示观测时间t吨通过、和第页t吨(·|·)表示系统随时间的跃迁密度t吨(见方程式(2)). 我们将写Δ=t吨+1t吨和过滤密度第页(x个t吨|1:)将用表示π(x个t吨),其中根据标准惯例1:=(1,…,). 为了简化符号,当我们在下面引入加权粒子时,我们将粒子和权重都下标为而不是t吨.

我们的目的是计算过滤密度π(x个t吨)递归地。基本概率计算得出这些密度的以下标准过滤递归:

π+1(x个t吨+1)(f)(+1x个t吨,x个t吨+1)第页Δ(x个t吨+1x个t吨)π(x个t吨)d日x个t吨.
8

粒子过滤器近似值π(x个t吨)由离散分布表示π^(X(X)t吨),其支持是一组N个颗粒,{x个(j个)}j个=1N个,具有关联的概率权重{w个(j个)}j个=1N个.替换π^(X(X)t吨)对于π(x个t吨)式(8)中给出了(连续密度)近似值π+1(x个t吨+1),

π˜+1(x个t吨+1)j个=1N个w个(j个)(f)(+1x个(j个),x个t吨+1)第页Δ(x个t吨+1x个(j个)).
9

粒子滤波算法的一次迭代的目的是构建进一步的粒子(离散分布)近似π˜+1(x个t吨+1).

我们可以通过重要性抽样获得这样的粒子近似,而皮特和谢泼德的辅助粒子滤波器给出了实现这一点的一般框架(1999). 我们选择表格的提案密度

j个=1N个β(j个)q个(x个t吨+1x个(j个),t吨+1).
10

选择合适的提案,即选择β(j个)s和q个,将在第节中对我们的具体应用进行分析5.

在某个时间模拟新粒子t吨+1我们

  • (a)

    模拟粒子x个(k个)时间,其中k个是一个离散随机变量的实现,其值为j个∈ {1,2,…,N个}有可能β(j个)、和

  • (b)

    在time中模拟新粒子t吨+1q个(x个t吨+1x个(k个),+1).

指定给这对粒子的权重(x个(k个),x个t吨+1)与…成比例

w个(k个)(f)(+1x个(k个),x个t吨+1)第页Δ(x个t吨+1x个(k个))β(k个)q个(x个t吨+1x个(k个),+1).
11

这是重复的N个时间产生加权粒子集的时间t吨+1,t吨+1,{(x个+1(j个),w个+1(j个))}j个=1N个,它提供了一个重要的采样近似值π+1(x个t吨+1). 可以重新规范权重,但不会对方法或其准确性产生重大影响。可以改进步骤(a)中的独立抽样:参见Carpenter的分层抽样思想等。(1999). 由此产生的粒子过滤器具有良好的理论特性,包括一致性(Crisan,2001)和CLT用于估计后验矩(Del Moral和Miclo,2000; 肖邦,2004; 昆施,2005),作为N个→∞. 在与初始条件指数遗忘相关的条件下,粒子滤波误差稳定为n个→∞ (德尔·道德和吉奥内特,2001; 昆施,2005).

当信号X(X)扩散过程是指过渡密度第页Δ(x个t吨+1x个(k个))由于方程中的期望项,对于大多数感兴趣的扩散来说,表达式(11)中出现的是难以处理的(2). 此外,对于观测模型(iii)(也适用于更一般的模型),似然项(f)(+1x个(k个),x个t吨+1)表达式(4)中给出的不能进行解析计算。

我们通过给每个新粒子分配一个随机权重来规避这些问题,该权重是平均值为表达式(11)的随机变量的实现。该随机变量的构造和模拟见第节4,它基于方程中过渡密度的特定表达式(2). 在比本文考虑的更一般的情况下,用正无偏估计量替换权重是一种有趣的可能性。事实上,在第节3.2我们表明,这种方法相当于用辅助变量方便地增加状态。

3.1. 重量模拟

在所有模型中,与该对相关联的重量(X(X)(k个),X(X)t吨+1)等于

小时+1(x个(k个),x个t吨+1,+1)μ(x个(k个),x个t吨+1,t吨,t吨+1)
12

哪里小时+1是已知函数,对于0<u个<t吨,

μ(x个,z(z),u个,t吨):=E类[经验{u个t吨()d日}].

期望是针对d日-维布朗桥,开始于u个u个=x个并按时完成t吨t吨=z(z).

3.1.1. 模型(i)和(ii)

对于型号(i)和(ii)

小时+1(x个(k个),x个t吨+1,+1)=w个(k个)(f)(+1x个t吨+1)N个Δ(x个t吨+1x个(k个))经验{一个(x个t吨+1)一个(x个(k个))}β(k个)q个(x个t吨+1x个(k个),+1),

=ϕ.在模型类型(ii)中,提案分发q个(x个t吨+1x个(k个),+1)应选择仅建议以下值x个t吨+1这样的话ζ(x个t吨+1)=+1; 然后(f)(+1|x个t吨+1)=1.

3.1.2. 模型(iii)

表达式(2)、(4)和(11)的合成=ϕ+ν,给出

小时+1(x个(k个),x个t吨+1,+1)=w个(k个)ν(x个t吨+1)N个Δ(x个t吨+1x个(k个))经验{一个(x个t吨+1)一个(x个(k个))}β(k个)q个(x个t吨+1x个(k个),+1).

章节4显示如何为每对(x、 z(z))和时间(u个,t吨),使用u个<t吨,其他辅助的变量V(V)、和函数第页(v(v),x、 z(z),u个,t吨)0,具有E类[第页(V(V),x个,z(z),u个,t吨)x个,z(z)]=μ(x个,z(z),u个,t吨)辅助变量根据适当的条件分布进行模拟(·|x、 z(z),u个,t吨)、和第页易于评估。我们的方法在权重中取代了难处理的术语μ及其无偏估计第页.

3.1.3. 随机加权粒子滤波器

  • 第1步:模拟样本x个0(1),,x个0(N个)第页(x个0),并设置w个0(j个)=1/N个.

对于=0,…,n个−1和j个=1,…,N个执行以下步骤。

  • 第2步:计算{β(k个)},嵌入式安全子系统={Σk个=1N个(β(k个))2}1; 如果是ESS<C类,对于某些固定常数C类,模拟k个,j个第页(k个)=β(k个),k个=1,…,N个、和设置δ+1(j个)=1; 否则设置k个,j个=j个δ+1(j个)=β(j个).

  • 步骤3:模拟x个+1(j个)q个(x个t吨+1x个(k个,j个),+1).

  • 第4步:模拟v(v)+1~q个(x个(k个,j个),x个t吨+1,t吨,t吨+1).

  • 第5步:指定粒子x个+1(j个)砝码
    w个+1(j个)=δ+1(j个)小时+1(x个(k个,j个),x个+1(j个),+1)第页(v(v)+1,x个(k个,j个),x个t吨+1,t吨,t吨+1).
    13

请注意,此算法包含关于是否在步骤2中传播之前对粒子重新采样的决定,该决定基于β(j个).常数C类可以解释为可接受的最小有效样本量。(见刘和陈(1998)以了解基于这种条件的重采样的基本原理。)是否进行重采样将影响赋予新粒子集的权重,这由以下不同的值来解释δ+1(j个)在步骤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日,具有非齐次马尔可夫转移密度

第页+1(z(z)+1,v(v)+1z(z),v(v))=第页(v(v)+1,z(z),z(z)+1,t吨,t吨+1)

(该密度相对于Leb×)和观测数据具有观测密度(f)(+1|z(z),z(z)+1). 通过构造Z轴在这个模型中π(x个t吨)即表达式(8)中的过滤密度。考虑一个应用于此模型的辅助粒子滤波器,其中我们以概率进行选择β(j个)每个现有粒子(z(z)(j个),v(v)(j个)),并根据提案生成新粒子

(z(z)+1,v(v)+1)~q个(z(z)+1z(z)(k个),+1)( d日v(v)+1z(z)(k个),z(z)+1,t吨,t吨+1)勒布(d日z(z)+1),

哪里q个与表达式(10)中的建议密度相同。在这个离散时间模型中,与每个粒子相关联的权重是可处理的,并由方程式给出(13). 因此,加权样本{(z(z)+1(j个),w个+1(j个))}j个=1N个是精确的粒子近似π+1(x个t吨+1),RWPF等价于该离散时间模型上的辅助粒子滤波器,其潜在结构已被辅助变量增强V(V).

这种等效表示揭示了我们方法的许多方面。首先,很明显,将正无偏估计量的多个实现平均化是低效的μ每粒子。相反,只需为每对粒子模拟一个估计器实现,就可以更有效地生成更多粒子。

其次,说明了RWPF结合了自举和辅助粒子滤波的优点。虽然很容易从概率分布进行模拟(如第节所述4),很难推导出其密度(相对于适当的参考测量值)。自从V(V)根据该测度传播s,避免了其计算。这是bootstrap过滤器的一个吸引人的特性,它可以传播粒子,而不需要分析系统跃迁密度。然而Z轴s通过用户指定的密度完成,该密度将信息纳入数据中。

第三,它表明RWPF将与应用于离散时间模型的辅助粒子滤波器具有相似的理论特性。这在第节中进行了探讨3.3.

3.3. 理论性质

考虑某些函数的后验均值的估计ϕ当时的状态t吨,E类[ϕ(x个t吨)|1:]. 研究粒子滤波器有效性的一种自然方法是将算法的极限行为视为N个→∞。对于标准辅助颗粒过滤器,肖邦(2004)引入了一种CLT来估计这类期望。本CLT直接适用于EPPF和ESPF。

附录F我们扩展了肖邦的结果(2004)并给出了一个关于RWPF中随机权重的进一步必要条件,在该条件下,CLT仍然成立。这个额外条件是条件2。估计量方差的表达式E类[ϕ(x个t吨)|1:]用RWPF得到的结果与标准情况下(即当权重已知时)的表达式不同,因为权重的随机性导致了一个额外的项(见方程式(27)–(29)以及对中定理3的评述附录F更多详细信息)。观察到RWPF可以重新表示为增强状态的标准粒子滤波器,这有助于肖邦方法的快速适应(参见第节3.2).

这种CLT的一个重要结果是,估计误差E类[ϕ(x个t吨)|1:]符合规定N个−1/2。我们考虑的扩散问题的先前过滤方法基于

  • (a)

    离散时间及其引入M(M)每个观测时间之间的中间时间点,

  • (b)

    使用欧拉或更高阶的扩散近似(1),以及

  • (c)

    对这个近似的离散时间模型应用粒子或其他过滤器。

参见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.广义泊松估计

我们已经激发了对正无偏估计量的模拟需求

E类[E类]  E类=:经验{0t吨()d日},
14

期望是针对d日-维布朗桥在本节中,我们介绍了一种推导此类估计量的方法,并提供了有关建议估计量方差的理论和模拟结果。这些结果在粒子滤波之外具有独立的意义,因此我们以通用的方式介绍我们的方法,其中是一个任意函数,假设它只在上是连续的R(右)d日。我们假设0=x个t吨=z(z),用于任意x、 z(z)R(右)d日t吨>通过布朗桥的时间同质性,我们的方法扩展到积分极限变为u个u个+t吨,对于任何u个>0.

贝斯科斯等。(2006b) 提出了期望(14)的无偏估计量,PE:

经验{(λc(c))t吨}λκj个=1κ{c(c)(ψj个)};
15

κ是具有平均值的泊松随机变量λt吨,的ψj个s均匀分布于[0,t吨]、和c(c)R(右)λ>0是任意常数。(这里和下面我们假设空产品,即κ=0,取值1。)PE的两个主要缺点是它可能返回负估计值,并且它的方差不一定是有限的。有界。然而,在我们的背景下,这是一个非常有限制性的假设。因此,在这里,我们引入了一组无偏和正估计E类[E类](14)推广了PE。我们考虑的方法允许c(c)λ依赖和许可证κ具有一般离散分布。首先,我们需要能够模拟随机变量L(左)U型具有

L(左)()U型,为所有人[0,t吨]
16

并且能够模拟在任何时候给定不等式(16)所隐含的条件。对于无界这是非平凡的。然而,由于在Beskos中引入了一种有效的算法,这两种模拟都变得可行等。(2005a) ●●●●。结构概述见附录A.

U型L(左)满足条件(16)和ψj个,j个1是[0,t吨]. 然后,估计器(14)可以重新表示为

E类[经验(U型t吨)经验[0t吨{U型()}d日]]                                          =E类[经验(U型t吨)k个=01k个![0t吨{U型()}d日]k个]=E类[经验(U型t吨)E类[k个=0t吨k个k个!j个=1k个{U型(ψj个)}U型,L(左)]]=E类[经验(U型t吨)t吨κκ!第页(κU型,L(左))j个=1κ{U型(ψj个)}],
17

哪里κ是具有条件概率的离散随机变量P(P)(κ=k个|U型,L(左))=第页(k个|U型,L(左)). 上述论点中的第二个等式是通过使用支配收敛和Fubini定理(由和的正性所支持)获得的。

我们可以推导出E类[E类](14)通过指定第页(·|U型,L(左)). 所有此类估算器系列将被称为GPE:

经验(U型t吨)t吨κκ!第页(κU型,L(左))j个=1κ{U型(ψj个)}.
18

以下定理(在附录B)为以下项目提供了最佳选择第页(·|U型,L(左)).

定理1。给定GPE的条件二阶矩U型L(左)

经验(2U型t吨)k个=0t吨k个第页(k个U型,L(左))k个!2E类[[0t吨{U型()}2 d日]k个U型,L(左)].
19

如果

k个=0t吨k个/2k个!E类[[0t吨{U型()}2 d日]k个U型,L(左)]1/2<,
20

然后通过选择将第二个力矩最小化

第页(k个U型,L(左))t吨k个/2k个!E类[[0t吨{U型()}2 d日)k个U型,L(左)]1/2,
21

最小二阶矩由

(经验(U型t吨)k个=0t吨k个/2k个![[0t吨{U型()}2 d日]k个U型,L(左)]1/2)2<,                                                                                                       对于几乎全部的U型L(左).
22

虽然表达式(21)的右侧无法进行分析评估,但它可以指导选择合适的第页(·|U型,L(左)). 如果已知,最优方案为平均泊松

λ:=[t吨0t吨{U型()}2 d日]1/2.
23

我们将讨论表达式(23)用于选择好提案的两种可能方式。

保守的方法需要第页(·|U型,L(左))为平均泊松(U型L(左))t吨(的上界λ). 我们称此估计器为GPE-1。GPE-1的一个优点是它的二阶矩有界于E类[出口(-2L(左)t吨)]. 因此,在温和明确的条件下,包含在以下定理中(在附录C),保证估计量的方差是有限的。

定理2。GPE-1具有有限方差的一个充分条件是

(u个1,,u个d日)δ=1d日(1+|u个|),  为所有人u个R(右),1d日,δ0

λ是随机的,另一种方法是引入(外生)随机平均值,并假设第页(·|U型,L(左))是具有此随机平均值的泊松。对于可控制性,我们选择随机平均值以获得伽马分布,当第页(·······································································U型,L(左))变为负二项分布(估计量GPE-2):

经验(U型t吨)t吨κΓ(β)(β+γ)β+κΓ(β+κ)ββγκj个=1κ{U型(ψj个)},
24

哪里γβ分别表示负二项分布的平均值和离散参数。由于负二项分布的尾部比PE重,因此只要PE有有限方差,GPE-2就会有有限方差。然而,如果γ选择为近似值E类[λ|U型,L(左)]. 有各种各样的特别的可以提供此期望的粗略估计的方法。应用Jensen不等式交换方程中的平方幂积分(23),然后近似E类[()|U型w个,L(左)w个]通过(E类[]),建议服用

γ=t吨U型0t吨(x个t吨t吨+t吨)d日>0
25

模拟研究(部分内容见第节4.1)表明这种选择在实践中非常有效,GPE-2的方差比PE或GPE-1小几个数量级。积分通常很容易计算;否则,可以使用粗略的近似值。

我们的介绍仅限于表达式(14)中的期望是关于布朗桥测度的情况。然而,正如Beskos指出的那样等。(2006b) 当对任意扩散桥测度取期望值时,只要能从该测度模拟出精确的骨架,就可以用完全相同的方法构造PE。GPE也可以在这个更广泛的框架中实施,前提是可以构造为满足条件(16)。

4.1. 模拟研究

我们考虑一个光滑有界测试函数(u个)={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更有效,特别是考虑到E类[κ]. 一般来说,PE的性能对c(c)λGPE-1的效率通常低于GPE-2。1也给出了var的值(E类)它的值比PE、GPE-1或GPE-2中的任何一个都要小得多(几个数量级),说明了这些辅助变量构造的绝对效率成本。

表1

表达式(14)四个估计量方差的蒙特卡罗估计,其中(u个)={罪恶(u个)2+科斯(u个)+1}/2

估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
方差体育课0.2020.2000.027
GPE-1级4.21×10−30.2080.034
GPE-2级2.08×10−30.2200.033
无功功率,无功功率(E类)3.74×10−53.27×10−34.72×10−3
E类[κ]体育课1.1181.1261.121
GPE-1级0.1301.0910.744
GPE-2级0.1190.3290.735
估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
方差体育课0.2020.2000.027
GPE-1级4.21×10−30.2080.034
GPE-2级2.08×10−30.2200.033
无功功率,无功功率(E类)3.74×10−53.27×10−34.72×10−3
E类[κ]体育课1.1181.1261.121
GPE-1级0.1301.0910.744
GPE-2级0.1190.3290.735

为了进行比较,我们还提供了var(E类). 我们还报告了E类[κ]. 我们考虑三对不同的起点和终点(x个,z(z))和时间增量t吨=1.估计值来自10个样本4实现。

表1

表达式(14)四个估计量方差的蒙特卡罗估计,其中(u个)={罪恶(u个)2+科斯(u个)+1}/2

估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
方差体育课0.2020.2000.027
GPE-1级4.21×10−30.2080.034
GPE-2级2.08×10−30.2200.033
无功功率,无功功率(E类)3.74×10−53.27×10−34.72×10−3
E类[κ]体育课1.1181.1261.121
GPE-1级0.1301.0910.744
GPE-2级0.1190.3290.735
估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
方差体育课0.2020.2000.027
GPE-1级4.21×10−30.2080.034
GPE-2级2.08×10−30.2200.033
无功功率,无功功率(E类)3.74×10−53.27×10−34.72×10−3
E类[κ]体育课1.1181.1261.121
GPE-1型0.1301.0910.744
GPE-2级0.1190.3290.735

为了进行比较,我们还提供了var(E类). 我们还报告了E类[κ]. 我们考虑三对不同的起点和终点(x个,z(z))和时间增量t吨=1.估计值来自10个样本4实现。

我们还研究了PE和GPE-2的效率如何随时间增量而变化t吨尤其是对于小型t吨(未显示结果)。这些经验结果表明,PE和GPE-2的误差变异系数为O(运行)(t吨δ)对一些人来说δ>0,但其值为δ这两个估计值不同。在我们调查的案例中,GPE-2似乎比PE具有更快的收敛速度。

此模拟研究的结果已被其他功能验证(结果未显示)。我们已经对可微性进行了实验(例如。(u个)=u个)和不可微(例如。(u个)=|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多得多,以构建它们的估计。

表2

蒙特卡洛估计基于104平方根误差除以四个估计量真值的实现第页t吨(z(z)|x个),表达式(6)的t吨=1和各种x个z(z)

估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
体育课1.250.930.17
GPE-2级0.130.780.2
DG-1公司0.50.450.3
DG-5型0.280.190.22
估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
体育课1.250.930.17
GPE-2级0.130.780.2
DG-1型0.50.450.3
DG-5型0.280.190.22

作为真值,我们取GPE-2给出的估计值的平均值。对于DG-1和DG-5,每个估计器使用的中间点数量分别为1和5;表中给出了PE和GPE-2的布朗桥模拟次数1(E类[κ]).

表2

蒙特卡洛估计基于104平方根误差除以四个估计量真值的实现第页t吨(z(z)|x个),表达式(6)的t吨=1和各种x个z(z)

估算员以下(x,z)对的结果:
x个=0时,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
体育课1.250.930.17
GPE-2型0.130.780.2
DG-1公司0.50.450.3
DG-5型0.280.190.22
估算员以下(x,z)对的结果:
x个=0,z(z)=0x个=0,z(z)=πx个=π,z(z)=π
体育课1.250.930.17
GPE-2级0.130.780.2
DG-1公司0.50.450.3
DG-5型0.280.190.22

作为真值,我们取GPE-2给出的估计值的平均值。对于DG-1和DG-5,每个估计器使用的中间点数量分别为1和5;表中给出了PE和GPE-2的布朗桥模拟次数1(E类[κ]).

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,我们选择了β(k个)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个k个,j个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

四种粒子滤波算法估计滤波平均值的相对效率E类[X(X)t吨|1:](与RWPF2相比,每行给出了一种算法的相对效率(图解的); 详见正文):图解的,RWPF1;图解的ESPF;图解的、EPPF

改变模型的参数和算法的实现将影响算法的相对性能。尤其是增加或减少σ2测量误差的方差将分别增加或降低EPPF相对于其他滤波器的相对效率。类似的结果出现在Δ分别减少或增加。其他三种算法的相对性能似乎对此类变化更为稳健。我们考虑实施EPPFβ(k个)=w个(k个),并使用正弦扩散的欧拉近似而非Ozaki近似来构建RWPF1和RWPF2的建议分布,但这些变化对方法的性能没有任何显著影响。我们还考虑减少重采样频率,设置C类=N个/RWPF算法第2步中的4(因此,当β(j个)s小于N个/4) 这大大降低了算法的性能(RWPF1和RWPF2的系数为2)。

我们还研究了增加观测之间的时间量Δ的影响。为此,我们使用了上述数据

  • (a)

    每10或

  • (b)

    每20天

时间点。

为了测量这些不同场景下过滤器的性能,我们使用了Carpenter的有效样本大小ESS等。(1999). ESS是根据滤波器独立运行期间后验均值估计值的方差来计算的,但将此方差与后验方差进行比较,以确定从后验值中独立提取多少个参数可以产生相同精度水平的估计值。我们着重于估计观测时间状态的后验平均值;如果2是粒子滤波器估计值的样本方差E类[X(X)t吨|1:]100次独立运行,以及σ^2是var的估计值(X(X)t吨|1:),则ESS为σ^2/2注意,通过ESS比较滤波器相当于根据估计量的方差比较滤波器。

给出了各种Δ值的ESS值。我们看到,随着Δ的增加,ESS值急剧下降,并且当Δ=20时,滤波器效率低下。这种性能下降是由于在这种情况下随机权重的巨大可变性。这些权重的可变性是由于

表3

观测值之间不同时间间隔的滤波器平均ESS值比较(Δ)

过滤器以下值的结果Δ:
Δ=10Δ=20
RWPF2型735
离散化8012
伪RWPF2923933
过滤器以下值的结果Δ:
Δ=10Δ=20
RWPF2型735
离散化8012
伪RWPF2923933

结果是使用GPE-2(RWPF2)的RWPF,这是一种通过离散扩散过程(离散化)数值近似权重的滤波器,以及在单位时间间隔引入非信息性观测值后的RWPF2。

表3

观测值之间不同时间间隔的滤波器平均ESS值比较(Δ)

过滤器以下值的结果Δ:
Δ=10Δ=20
RWPF2型735
离散化8012
伪RWPF2923933
过滤器以下值的结果Δ:
Δ=10Δ=20
RWPF2型735
离散化8012
伪RWPF2923933

结果是使用GPE-2(RWPF2)的RWPF,这是一种通过离散扩散过程(离散化)数值近似权重的滤波器,以及在单位时间间隔引入非信息性观测值后的RWPF2。

  • (a)
    的可变性
    经验{t吨t吨+1()d日},
    26
    跨越不同的扩散路径,以及
  • (b)

    估计给定路径的蒙特卡罗可变性。

为了评估由(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吨|.

图3

(a) 示例2中Cox过程的模拟和RWPF分析的结果(图解的,潜在扩散的绝对路径;图解的,观察到的到达时间;图解的,过滤RWPF中的估计值;图解的,扩散绝对值的90%可信区间(尽管为清晰起见,它们始终显示,可信区间仅在显示滤波器估计值的时间计算并适用),以及(b)RWPF权重的ESS,定义为(Σj个=1N个w个(j个))2/Σj个=1N个w个(J型)2,随着时间的推移(ESS的急剧增加对应于重采样时间)

我们对RWPF的实施是基于从先前分布中提出粒子,所以β(k个)=w个(k个)q个(x个t吨+1x个j个,+1)就是OU跃迁密度第页(x个t吨+1x个j个).我们用GPE-2模拟了随机权重。我们计算了每个观测时间以及所选的56个伪观测时间的滤波密度,因此计算滤波密度的两个连续时间之间的最大时间差为0.1。这是必要的,以避免需要布朗桥模拟来模拟权重对于长时间的互观测而言太大的次数,并控制随机权重的方差(见上文)。这些非观测时间的似然函数通过删除ν(x个t吨)来自表达式(4)。

我们将粒子数设置为1000,并在β(j个)s小于100(C类=N个/第节算法步骤2中的10). 虽然正弦扩散的结果表明,这将导致算法重采样次数太少,但我们选择了一个较低的阈值,以便通过粒子过滤器权重的ESS如何随时间衰减来监测粒子过滤器的性能。该滤波器一次运行的结果如图所示。(a) ●●●●。该方法的计算效率可以通过图来衡量。(b) 其中的ESSw个(j个)s随时间绘制。

6.讨论

我们描述了如何在粒子滤波器中使用最新的方法来精确模拟扩散和无偏估计扩散指数泛函,从而避免了时间离散化的需要。在我们介绍的方法中,特别关注了RWPF,它实现了一个辅助粒子过滤器,但模拟了分配给每个粒子的权重。我们表明,这种方法相当于应用于适当扩展模型的辅助粒子滤波器。我们预计,这种方法将在不同于本文所考虑的模型的模型中有有趣的应用,然而,这些模型涉及棘手的动力学或可能性。

我们专注于滤波问题,根据迄今为止的观察结果估计当前状态。然而,对预测的扩展是微不足道的,只需要从状态方程进行模拟的能力,这可以通过EA算法实现。使用北川的想法也很简单(1996),其中每个粒子存储其轨迹的历史,以获得平滑密度的近似值(给定到目前为止的观测值,过去某个时间的状态密度)。

尽管粒子只存储每个观测时间的状态值,但填充这些时间之间的扩散路径以在任何时候产生关于状态的推断都是简单的。分布的粒子近似(X(X),t吨−1<<t吨)以数据为条件1:可以通过使用当前的一组加权粒子来构建{(x个1(j个),x个(j个))}j个=1N个带砝码{w个(j个)},如下所示。首先我们需要介绍一些符号;我们用表示x个1|(j个)粒子在时间上的值t吨−1从中j个第th个粒子t吨下降。粒子近似由一组加权路径给出{(x个S公司,t吨1<<t吨)(j个)}j个=1N个带砝码{w个(j个)}。每条路径都是一个扩散桥,从x个1|(j个)并于完成x个(j个)可以使用EA进行模拟,如Beskos所述等。(2006年a,2005a) ●●●●。在观测区域(i)和(ii)中,应用EA模拟具有密度的扩散桥,其布朗桥测度由下式给出经验{t吨1t吨ϕ(X(X))d日},而在区域(iii)中,相应的密度为经验[t吨1t吨{ϕ(X(X))+ν(X(X))}d日]。可以直接利用此表示来推断X(X)在观察时间之间。

致谢

作者感谢联合编辑和裁判的宝贵意见。

参考文献

1

阿伊特萨赫勒
,
年。
(
2004
)
多元扩散的闭式似然展开
.工作文件.
普林斯顿大学本德海姆金融中心
,普林斯顿大学。(可从http://www.princeton.edu/~yacine/research.htm.)

2

贝斯科斯
,
答:。
,
帕帕斯皮利奥普洛斯
,
O。
罗伯茨
,
总干事。
(
2005年a
)
扩散测度的一种新分解和有限样本路径构造
。待发布。

贝斯科斯
,
答:。
,
帕帕斯皮利奥普洛斯
,
O。
罗伯茨
,
总干事。
(
2005年b
)
离散观测扩散过程的Monte Carlo极大似然估计
。待发布。

4

贝斯科斯
,
答:。
,
帕帕斯皮利奥普洛斯
,
O。
罗伯茨
,
总干事。
(
2006年a
)
扩散样品路径的回顾性精确模拟及其应用
.
伯努利
,
12
,
1077
1098
.

5

贝斯科斯
,
答:。
,
帕帕斯皮利奥普洛斯
,
O。
,
罗伯茨
,
总干事。
羽毛头
,
第页。
(
2006年b
)
离散观测扩散过程的精确且计算效率高的似然估计(附讨论)
.
J.R.统计。Soc.B公司
,
68
,
333
382
.

6

贝斯科斯
,
答:。
罗伯茨
,
总干事。
(
2005
)
扩散的精确模拟
.
附录申请。普罗巴伯。
,
15
,
2422
2444
.

7

木匠
,
J。
,
克利福德
,
第页。
羽毛头
,
第页。
(
1999
)
非线性问题的改进粒子滤波
.
IEE程序。雷达、声纳导航
,
146
,
2
7
.

8

芯片
,
美国。
,
皮特
,
M.K.博士。
谢泼德
,
N。
(
2006
)
扩散驱动状态空间模型的似然推理
.
工作文件
.
纳菲尔德学院
牛津大学。(可从http://www.nuff.ox.ac.uk/users/shephard网站/.)

9

肖邦
,
N。
(
2004
)
序列蒙特卡罗方法的中心极限定理及其在贝叶斯推理中的应用
.
Ann.Statist公司。
,
32
,
2385
2411
.

10

克里桑
,
D。
(
2001
)颗粒过滤器——理论观点。
序贯蒙特卡罗方法在实践中的应用
(编辑
答:。
水龙头
,
N。
德弗里塔斯
N。
戈登
),第页。
17
41
纽约:
施普林格
.

11

克里桑
,
D。
,
德尔莫拉尔
,
第页。
里昂
,
T·J。
(
1999
)
Kushner-Stratonovich方程的相互作用粒子系统近似
.
高级申请。普罗巴伯。
,
31
,
819
838
.

12

达昆哈城堡
,
D。
佛罗伦萨-兹米罗
,
D。
(
1986
)
从离散观测值估计扩散系数
.
随机性
,
19
,
263
284
.

13

达西奥斯
,
答:。
,
高-宽。
(
2005
)
具有散粒噪声强度的Cox过程驱动线性系统的卡尔曼-布基滤波及其在再保险合同定价中的应用
.
J.应用。普罗巴伯。
,
42
,
93
107
.

14

德尔莫拉尔
,
第页。
吉奥内
,
答:。
(
2001
)
交互过程的稳定性及其在滤波和遗传算法中的应用
.
Ann.Inst.H.Poincare概率。统计师。
,
37
,
155
194
.

15

德尔莫拉尔
,
第页。
,
加科德
,
J。
普罗特
,
第页。
(
2001
)
离散时间观测滤波的Monte-Carlo方法
.
普罗巴伯。理论相关Flds
,
120
,
346
368
.

16

德尔莫拉尔
,
第页。
米克罗
,
L。
(
2000
)
分支和相互作用粒子系统:Feymann-Kac公式的近似及其在非线性滤波中的应用
.
半概率。斯特拉斯。
,
34
,
1
145
.

17

水龙头
,
答:。
,
德弗里塔斯
,
N。
戈登
,
N。
(
2001
)序贯蒙特卡罗方法简介。
序贯蒙特卡罗方法在实践中的应用
,第页。
14
纽约:
施普林格
.

18

达菲
,
D。
辛格尔顿
,
K·J。
(
1999
)
可违约债券的期限结构建模
.
财务版次。螺柱。
,
12
,
687
720
.

19

达勒姆
,
G.B.公司。
盖朗
,
A.R.公司。
(
2002
)
连续时间扩散过程最大似然估计的数值技术(附注释)
.
J.公交车。经济学。统计师。
,
20
,
297
338
.

20

恩格尔
,
无线电频率。
(
2000
)
超高频数据的计量经济学
.
计量经济学
,
68
,
1
22
.

21

Golightly公司
,
答:。
威尔金森
,
D.J.博士。
(
2006
)
非线性多元扩散的贝叶斯序列推断
.
统计人员。计算。
,
16
,
323
338
.

22

戈登
,
N。
,
鲑鱼
,
D。
史密斯
,
A.F.M.公司。
(
1993
)
非线性/非高斯贝叶斯状态估计的新方法
.
IEE程序。F类
,
140
,
107
113
.

23

离子
,
E.公司。
(
2003
)
部分观测扩散过程的序贯蒙特卡罗推断与滤波
.工作文件.
密歇根大学
安娜堡。(可从http://www.stat.lsa.umich.edu/~ionides/pubs/工作纸过滤器.pdf.)

24

北川
,
G.公司。
(
1996
)
非高斯非线性状态空间模型的蒙特卡罗滤波和平滑器
.
J.计算图表。统计师。
,
5
,
1
25
.

25

,
南卡罗来纳州。
,
,
X·S。
线路接口单元
,
J.S.公司。
(
2005
)
单分子实验数据的贝叶斯分析(附讨论)
.
申请。统计师。
,
54
,
469
506
.

26

昆施
,
H.R.公司。
(
2005
)
蒙特卡罗滤波器:算法和理论分析
.
Ann.Statist公司。
,
33
,
1983
2021
.

27

线路接口单元
,
J.S.公司。
,
R。
(
1998
)
动态系统的序贯蒙特卡罗方法
.
《美国统计杂志》。助理。
,
93
,
1032
1044
.

28

尼古拉
,
J。
(
2002
)
一种模拟随机微分方程可能性的新方法
.
计量经济学。J。
,
5
,
91
103
.

29

皮特
,
M.K.博士。
谢泼德
,
N。
(
1999
)
通过模拟进行过滤:辅助粒子过滤器
.
《美国统计杂志》。助理。
,
94
,
590
599
.

附录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日。此序列定义窗体示例空间的分区j个=1D类j个,其中路径属于D类j个当且仅当路径超过由一个j个−1但不是由一个j个.在Beskos等。(2005a) 演示了如何模拟随机变量,该变量决定了D类j个属于,以及如何模拟在这个随机变量的任何时间集合上:分层布朗桥结构。假设是连续的,知道D类j个足以确定U型L(左)满足条件(16)。事实上,在简化的设置中是有界的,正如示例1中的正弦扩散一样,可以避免分层布朗桥结构,因为它易于选择U型L(左)独立于.

附录B:定理1的证明

:=t吨κκ!第页(κU型,L(左))j个=1κ{U型(ψj个)}.

然后,表达式(19)建立如下:

E类[2U型,L(左)]=E类[E类[2κ,]]=E类[t吨2κ{κ!第页(κU型,L(左))}2[0t吨{U型()}2t吨 d日]κ]=E类[t吨κ{κ!第页(k个U型,L(左))}2E类[[0t吨{U型()}2 d日]κU型,L(左),κ]]=k个=0t吨k个第页(k个U型,L(左))k个!2E类[[0t吨{U型()}2 d日]k个U型,L(左)].

上面使用了Fubini定理和支配收敛性(这是有效的,因为被积函数几乎肯定是正的)。通过使用以下结果(可以通过使用Jensen不等式容易地证明)来获得表达式(22)。(f)>0用于=1,2,…. 然后是第页s最小化=0(f)/第页在约束∑下第页=1由下式给出第页=√(f)/Σ √(f).

附录C:定理证明2

GPE-1估计值小于或等于exp(−L(左)t吨)因此,如果E类[经验(−L(左)t吨)]<∞,其中期望是关于d日-维布朗桥x个时间为0到时间t吨然而,

E类[经验(L(左)t吨)]=0P(P){经验(L(左))>w个}d日w个                                  =0P(P){L(左)<日志(w个)}d日w个0P(P){δ=1d日(1+M(M))>日志(w个)}d日w个

哪里M(M)=支持0t吨||利用定理2中的增长界。此外,

0P(P){δ=1d日(1+M(M))>日志(w个)}d日w个0=1d日P(P){δ(1+M(M))>d日1日志(w个)}d日w个                                                                    =0=1d日P(P){M(M)>(d日δ)1日志(w个)1}d日w个.

因此,仍需约束d日这个表达式右边的积分。然而,根据Bachelier-Levy公式计算布朗运动和桥梁的撞击时间,

P(P)(M(M)>v(v))经验[2(v(v)最大值{x个,})2/t吨]+经验[2(最小值{x个,+v(v)})2/t吨]

等等

P(P){M(M)>(d日δ)1日志(w个)1}经验([2{(d日δ)1日志(w个)1}最大值{x个,}]2/t吨)+经验(2[n个{x个,+(d日δ)1日志(w个)1}]2/t吨)

它后退得像w个k个日志(w个)作为w个→∞, 从而得出结论。

附录D:精确传播粒子滤波器和精确模拟粒子滤波器,例如1

EPPF根据以下程序生成新粒子:

  • (a)

    选择当前粒子之一x个(k个,j个),其中粒子j个被概率选择β(j个);

  • (b)

    建议x个t吨+1从平均值的正态分布x个(k个,j个)和方差Δ;

  • (c)

    以概率exp{−cos接受此建议(x个t吨+1)−1}; 如果提案被拒绝,则返回(a);

  • (d)
    有可能接受这个建议
    E类[经验{0Δϕ()d日}],

    其中期望是关于布朗桥定律的0=x个(k个,j个)Δ=x个t吨+1、和φ在方程式中给出(7). 如果提案被拒绝,则返回(a);否则x个t吨+1是当时的新粒子t吨+1带重量w个+1=(f)(+1|x个+1).

步骤(d)通过使用Beskos中所述的回顾性抽样进行等。(2006a) ●●●●。ESPF如上所述进行,但步骤(a)和(b)替换为步骤

  • (a)提议(x个(k个,j个),x个t吨+1)根据密度成比例
    经验{科斯(x个(k个,j个))(+1x个(k个,j个))22(σ2+Δ)}经验{(x个t吨+1η)22τ2}

    哪里η=(σ2+Δ)1(x个(k个,j个)σ2+Δ+1)、和τ=σ2Δ/(σ2).

重复该算法,直到N个的值x个t吨+1已验收,每个重量为1/N个.

附录E:示例1的提案分发

考虑满足SDE(1)的扩散,其中d日=1表示简单。该SDE的Ozaki近似是基于某个值附近漂移的一阶泰勒级数展开x个对于例1的正弦扩散,我们得到了近似的SDE

d日X(X)˜=科斯(x个){x个棕褐色的(x个)X(X)˜}d日+d日B类.

所以X(X)˜{x个棕褐色的(x个)}是示例2中定义的OU进程ρ=cos(x个)和σ=1.计算q个(x个t吨+1x个(j个),+1)我们计算了由Ozaki近似给出的关于x个=x个(j个)和似然函数(f)(+1|x个t吨+1). 定义

τ2=[1经验{2科斯(x个(j个))Δ}]/2科斯(x个(j个))

η=x个(j个)棕褐色的(x个(j个))[1经验{科斯(x个(j个))Δ}]

我们得到了q个(x个t吨+1x个(j个),+1)平均值正常(ησ2++1τ2)/τ2σ2和方差η2τ2/(η2+τ2). 此外,我们计算β(j个)w个(j个)N个τ2+σ2(+1η).

附录F:中心极限定理

为了简化符号,我们考虑了粒子滤波器的一个特例,该特例与肖邦的粒子滤波器相似(2004). 我们为时间选择提案密度t吨+1拥有βj个=w个(j个); 我们假设X(X)t吨(j个)在步骤2中。肖邦的粒子滤波器(2004)每次拆分模拟粒子t吨+1进入之内

  • (a)

    粒子的时间重采样t吨

  • (b)

    这些粒子随时间的传播t吨+1.

我们的独立同分布抽样假设等价于肖邦的多项式重抽样情况(2004). (如果刘和陈的剩余采样方法相同,则CLT的条件相同(1998),但差异不同。)为了简单起见,我们考虑观察模型(i)或(ii),尽管结果很容易扩展到观察模型(iii)。

θ(j个)=(x个t吨(j个),x个t吨1(k个,j个)),其中k个,j个是在模拟j个-时间粒子t吨θ(j个)j个第th个粒子t吨和粒子一起t吨−1它是从那里传下来的。还让E类θ表示给定的条件期望θ类似地,让μ(θ)=μ(x个−1,x个,t吨−1,t吨),并表示为R(右)无偏估计量μ(θ),即。E类[R(右)]=μ(θ). 一个重要的数量是σ2(θ)=无功功率,无功功率(R(右)).

我们定义E类[ϕ]和var(ϕ)为任意函数的后验均值和方差ϕ(θ)时间,并考虑粒子滤波估计E类[ϕ]. π˜(θ)是密度第页(x个t吨−1|1:−1)q个(x个t吨|x个t吨−1). 最后定义E类q个[ϕ]和varq个(ϕ)是的条件期望和方差的简写ϕ(θ)关于q个(x个t吨|x个t吨−1)(这些是的功能x个t吨−1). 我们将‖·‖表示为欧几里德范数,并递归地定义Φ是一组可测量的函数ϕ这样,对一些人来说δ>0,E类,并且该函数x个t吨−1E类q个[小时μϕ]为Φ−1.

定理3。考虑一个函数ϕ; 定义V(V)˜0=无功功率,无功功率π˜(x个0)(φ)通过归纳,

V(V)˜(φ)=V(V)^1(E类q个[φ])+E类1[无功功率,无功功率q个(φ)],>0,
27
V(V)(φ)=V(V)˜{μ小时(φE类[φ])}+E类π˜[(φE类[φ])2σ2小时2]E类π˜[μ小时]2,0,
28
V(V)^(φ)=V(V)(ϕ)+无功功率,无功功率(ϕ),0
29

那么如果,为了所有人,

  • (a)

    x个t吨1属于Φ,

  • (b)

    E类π˜[小时2σ2]<

  • (c)

    E类π˜[σφ小时]2+δ<对一些人来说δ>0

对于任何ϕ∈ Φ,E类[ϕ]和V(V)(ϕ)是有限的,并且我们在分布上收敛为粒子数,N个,趋于∞:

N个1/2{j个=1N个w个(j个)φ(x个t吨(j个))j个=1N个w个(j个)E类[φ]}N个{0,V(V)(φ)}.

备注1。方程式(27)–(29)指迭代时由于传播、加权和重采样阶段而导致加权粒子方差的变化.仅方程式(28)与肖邦的各自成绩不同(2004)这是由于右侧的第二项,它表示由于权重的随机性而导致的方差增加。条件(a)取自肖邦(2004)并适用于标准颗粒过滤器;条件(b)和(c)是新的,是限制随机权重方差的条件,确保V(V)(ϕ)是有限的。

证明。我们改编了肖邦的感应证明(2004),依次考虑传播、加权和重采样步骤。我们的滤波器仅在加权步长方面与标准粒子滤波器不同,因此我们只需要调整肖邦引理A2的结果(2004). 事实上,方程式(27)和(29)与肖邦的对应量相同(2004); 因此,仍需显示等式(28). 我们定义常数K(K)=E类π˜[R(右)小时]ϕ*=R(右)小时(ϕE类[ϕ])/K(K)在扩大的信号空间框架内,我们可以应用方程式(4)肖邦的(2004),给予

V(V)(φ)=V(V)˜(φ*)=V(V)^1{E类q个[φ*]}+E类1[无功功率,无功功率q个(φ*)].

现在我们可以计算E类q个[ϕ*]首先对辅助变量取期望值(条件是θ). 这给了E类q个[ϕ*]=E类q个[μ国际卫生组织(ϕE类[ϕ])/K(K)]. 同样,我们获得

无功功率,无功功率q个(φ*)=无功功率,无功功率q个{E类[R(右)小时(φE类[φ])/K(K)]}+E类q个[无功功率,无功功率{R(右)小时(φE类[φ])/K(K)}]
30
=无功功率,无功功率q个{μ小时(φE类[φ])/K(K)}+E类q个[σ2小时2{φE类[φ]}2/K(K)2].
31

(这里是方程中的期望和方差(30)与辅助变量有关。)将这些结果结合起来得出方程式(28). 正则性条件(a)-(c)也直接转换。

本文根据牛津大学出版社标准期刊出版模式的条款出版和发行(https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)