Particle Filtering for Partially Observed Gaussian State Space Models

Andrieu, Christophe; Doucet, Arnaud

doi:10.1111/1467-9868.00363

总结

当后验分布通过数据积累随时间演化时，解决贝叶斯估计问题在动态模型中有许多应用。为了解决这些问题，最近提出了大量基于粒子滤波方法的算法，也称为序列蒙特卡罗算法。我们提出了一种特殊的粒子滤波方法，该方法使用正态分布的随机混合来表示部分观测到的高斯状态空间模型的后验分布。该算法基于边缘化思想来提高效率，与标准算法相比，可以带来显著的收益。它不同于以往仅适用于条件线性高斯状态空间模型的算法。对动态tobit和probit模型进行了计算机仿真，以评估所提算法的性能。

贝叶斯估计,过滤,广义线性时间序列,重要性抽样,顺序蒙特卡罗采样,状态空间模型

1.简介

1.1. 背景

许多数据分析任务涉及在只有部分或不准确观测数据可用时估计动态模型的状态（West和Harrison，1997). 除了在一些特殊情况下，包括线性高斯状态空间模型，在线状态估计是一个不允许闭式解的问题。由于大多数现实世界模型都是非线性和非高斯的，因此开发有效的计算方法来数值求解这一所谓的贝叶斯滤波问题是非常有意义的。

为了克服这个问题，人们提出了许多近似方案，如扩展卡尔曼滤波器。然而，在许多实际场景中，这些近似方法是不可靠的，并且很难在线诊断故障。最近，对于非线性或非高斯时间序列分析的序贯蒙特卡罗（SMC）方法（Doucet等。,2001). 这些方法始于戈登等。(1993)，使用基于随机样本（或粒子）的后验概率分布表示。

1.2. 一般问题

对于任何序列我_t吨，我们定义

我_{我 : j个} ≜ (我_{我}, 我_{我 + 1}, \dots, 我_{j个}) .

在本文中，我们将集中讨论以下一类状态空间模型。让t吨=1,2，…表示离散时间：则

{x个}_{t吨} = {A类}_{t吨} {x个}_{t吨 - 1} + B_{t吨} {v（v）}_{t吨} + {F类}_{t吨} {u个}_{t吨}, {x个}_{0} \sim N个 ({\hat{x个}}_{0}, {P（P）}_{0}),

(1)

年_{t吨} = {C类}_{t吨} {x个}_{t吨} + 天_{t吨} ε_{t吨} + {G公司}_{t吨} {u个}_{t吨},

(2)

{z（z）}_{t吨} \sim 对 ({z（z）}_{t吨} ∣ 年_{t吨}),

（3）

哪里u个_t吨∈ℝ^n个u个是一个外生过程x个_t吨∈ℝ^n个x个和年_t吨∈ℝ^n个年是未观察到的过程。序列

{v（v）}_{t吨} \overset{IID公司}{\sim} N个 (0, 我_{{n个}_{v（v）}}) \in ℝ^{{n个}_{v（v）}}

和

ε_{t吨} \overset{IID公司}{\sim} N个 (0, 我_{{n个}_{ε}}) \in ℝ^{{n个}_{ε}}

是独立同分布（IID）高斯。我们假设P（P）₀>0;x个₀,v（v）_t吨和w个_t吨都是相互独立的t吨、和模型参数

λ ≜ ({\hat{x个}}_{0}, {P（P）}_{0}, {A类}_{t吨}, B_{t吨}, {C类}_{t吨}, 天_{t吨}, {F类}_{t吨}, {G公司}_{t吨}; t吨 = 1, 2, \dots)

已知。过程(x个_t吨)和(年_t吨)定义一个标准的线性高斯状态空间模型。我们没有观察到(年_t吨)就我们而言，但是(z（z）_t吨). 这个观察(z（z）_t吨)在给定过程的条件下是独立的(x个_t吨)和(年_t吨)并根据对(z（z）_t吨|年_t吨); 假设对(z（z）_t吨|年_t吨)可以按点计算，直到得到一个规范化常数。通常情况下对(z（z）_t吨|年_t吨)属于指数族。或者z（z）_t吨可以是的审查或量化版本年_t吨。这类部分观测高斯状态空间模型有许多应用；许多例子都被讨论过，例如在德容(1997)、曼里克和谢泼德(1998)West和Harrison(1997).

我们想估算一下按时间顺序后验分布的一些特征对(x个_0:t吨|z（z）_1:t吨). 通常，我们对计算感兴趣E类(x个_t吨|z（z）_1:t吨)（过滤），E类(x个_{t吨+L（左）}|z（z）_1:t吨)（预测）和E类(x个_{t吨−L（左）}|z（z）_1:t吨)（固定滞后平滑），其中L（左）是一个正整数。这些估计一般不允许使用解析表达式，我们必须采用数值方法。

1.3. 分辨率

SMC方法松散地说，是重要采样和重采样方法的组合，它使我们能够随着时间的推移有效地传播一组大致按照对(x个_0:t吨|z（z）_1:t吨). 我们可以应用标准的SMC方法，如引导过滤器（Gordon等。,1993)估计对(x个_0:t吨,年_1:t吨|z（z）_1:t吨)因此对(x个_0:t吨|z（z）_1:t吨). 然而，在其标准形式中，该算法并没有使用模型的所有显著结构。我们的算法基于边缘化技术，通常称为Rao–Blackwellization方法（Gelfand and Smith，1990)，这提高了程序的效率。它侧重于对对(年_1:t吨|z（z）_1:t吨)而不是接缝密度对(x个_0:t吨,年_1:t吨|z（z）_1:t吨). 这个过程(x个_t吨)通过分析进行整合。一次对(年_1:t吨|z（z）_1:t吨)我们可以得到E类(x个_t吨|z（z）_1:t吨),E类(x个_{t吨+L（左）}|z（z）_1:t吨)和E类(x个_{t吨−L（左）}|z（z）_1:t吨)通过进一步讨论的卡尔曼滤波器。在马尔可夫链蒙特卡罗框架中，德容(1997)提出了所谓的扫描取样器对(年_1:t吨|z（z）_1:t吨)在类似类别的状态空间模型中；参见Manrique和Shephard(1998)对于某些应用程序。

1.4. 计划

论文的其余部分组织如下。章节2展示了如何限制我们自己进行估计对(年_1:t吨|z（z）_1:t吨)而不是对(x个_0:t吨,年_1:t吨|z（z）_1:t吨)从而提高了蒙特卡罗效率。然后详细描述了粒子滤波算法。章节三通过对动态tobit和probit模型的计算机仿真，验证了该算法的性能。

2.Rao–Blackwellized粒子过滤

2.1. 边缘化

考虑由方程定义的状态空间模型(1)–(三). 我们有

对 ({x个}_{0 : t吨} ∣ {z（z）}_{1 : t吨}) = {¦Β}^{​} 对 ({x个}_{0 : t吨} ∣ 年_{1 : t吨}) 对 (年_{1 : t吨} ∣ {z（z）}_{1 : t吨}) d日 年_{1 : t吨} .

因此，如果我们（通过进一步描述的SMC方法）获得与密度相关的概率分布的近似值对(年_1:t吨|z（z）_1:t吨)表单的

{\hat{对}}_{N个} (d日 年_{1 : t吨} ∣ {z（z）}_{1 : t吨}) = {\sum^{​}}_{我 = 1}^{N个 ​} {w个}_{t吨}^{(我)} δ_{{\tilde{年}}_{1 : t吨}^{(我)}} (d日 年_{1 : t吨}), {w个}_{t吨}^{(我)} ⩾ 0, {\sum^{​}}_{我 = 1}^{N个} {w个}_{t吨}^{(我)} = 1,

然后对(x个_0:t吨|z（z）_1:t吨)可以通过使用

{\hat{对}}_{N个} ({x个}_{0 : t吨} ∣ {z（z）}_{1 : t吨}) = {\sum^{​}}_{我 = 1}^{N个} {w个}_{t吨}^{(我)} 对 ({x个}_{0 : t吨} ∣ {\tilde{年}}_{1 : t吨}^{(我)}),

即高斯密度的混合。根据这样的近似值，我们可以估计E类(x个_t吨|z（z）_1:t吨)和E类(x个_{t吨−L（左）}|z（z）_1:t吨). 例如，估计E类(x个_t吨|z（z）_1:t吨)由提供

{\hat{E类}}_{N个} ({x个}_{t吨} ∣ {z（z）}_{1 : t吨}) = {¦Β}^{​} {x个}_{t吨} {\hat{对}}_{N个} ({x个}_{0 : t吨} ∣ {z（z）}_{1 : t吨}) d日 {x个}_{0 : t吨} = {\sum^{​}}_{我 = 1}^{N个} {w个}_{t吨}^{(我)} E类 ({x个}_{t吨} ∣ {\tilde{年}}_{1 : t吨}^{(我)}),

哪里 $E类 ({x个}_{t吨} ∣ {\tilde{年}}_{1 : t吨}^{(我)})$ 通过与方程定义的线性高斯状态空间模型相关联的卡尔曼滤波器进行计算(1)和(2). 使用方差分解公式，很明显，对于任何函数小时(·)

无功功率，无功功率 {小时 ({x个}_{t吨}) ∣ {z（z）}_{1 : t吨}} ⩾ 无功功率，无功功率 [E类 {小时 ({x个}_{t吨}) ∣ 年_{1 : t吨}, {z（z）}_{1 : t吨}} ∣ {z（z）}_{1 : t吨}],

这表明估计对(年_1:t吨|z（z）_1:t吨)只有这样才能更有效率。

获得边缘后部密度的SMC近似值对(年_1:t吨|z（z）_1:t吨)，我们需要能够逐点估计这个“目标”密度，直到一个规范化常数。我们有

对 (年_{1 : t吨} ∣ {z（z）}_{1 : t吨}) ⑪ {\prod^{​}}_{k个 = 1}^{我} 对 ({z（z）}_{k个} ∣ 年_{k个}) 对 (年_{k个} ∣ 年_{1 : k个 - 1}),

(4)

哪里

对 (年_{1} ∣ 年_{1 : 0}) ≜ 对 (年_{1}) .

作为对(z（z）_k个|年_k个)假设已知到标准化常数，则只需估计对(年_k个|年_1:k个−1)达到规范化常数。可以使用卡尔曼滤波器计算预测密度。

卡尔曼滤波方程如下。设置 ${x个}_{0 ∣ 0} = {\hat{x个}}_{0}$ 和P（P）_0|0=P（P）₀; 然后针对t吨=1,2，…计算

\begin{matrix} {x个}_{t吨 ∣ t吨 - 1} = {A类}_{t吨} {x个}_{t吨 - 1 ∣ t吨 - 1} + {F类}_{t吨} {u个}_{t吨}, \\ {P（P）}_{t吨 ∣ t吨 - 1} = {A类}_{t吨} {P（P）}_{t吨 - 1 ∣ t吨 - 1} {A类}_{t吨}^{T型} + B_{t吨} B_{t吨}^{T型}, \\ 年_{t吨 ∣ t吨 - 1} = {C类}_{t吨} {x个}_{t吨 ∣ t吨 - 1} + {G公司}_{t吨} {u个}_{t吨}, \\ {S公司}_{t吨} = {C类}_{t吨} {P（P）}_{t吨 ∣ t吨 - 1} {C类}_{t吨}^{T型} + 天_{t吨} 天_{t吨}^{T型}, \\ {x个}_{t吨 ∣ t吨} = {x个}_{t吨 ∣ t吨 - 1} + {P（P）}_{t吨 ∣ t吨 - 1} {C类}_{t吨}^{T型} {S公司}_{t吨}^{- 1} (年_{t吨} - 年_{t吨 ∣ t吨 - 1}), \\ {P（P）}_{t吨 ∣ t吨} = {P（P）}_{t吨 ∣ t吨 - 1} - {P（P）}_{t吨 ∣ t吨 - 1} {C类}_{t吨}^{T型} {S公司}_{t吨}^{- 1} {C类}_{t吨} {P（P）}_{t吨 ∣ t吨 - 1}, \end{matrix}}

（5）

哪里

\begin{array}{c} {x个}_{t吨 ∣ t吨 - 1} ≜ E类 ({x个}_{t吨} ∣ 年_{1 : t吨 - 1}) \\ {x个}_{t吨 ∣ t吨} ≜ E类 ({x个}_{t吨} ∣ 年_{1 : t吨}) \\ 年_{t吨 ∣ t吨 - 1} ≜ E类 (年_{t吨} ∣ 年_{1 : t吨 - 1}) \\ {P（P）}_{t吨 ∣ t吨 - 1} ≜ 覆盖（cov） ({x个}_{t吨} ∣ 年_{1 : t吨 - 1}) \\ {P（P）}_{t吨 ∣ t吨} ≜ 冠状病毒 ({x个}_{t吨} ∣ 年_{1 : t吨}) \end{array}

和

{S公司}_{t吨} ≜ 覆盖（cov） (年_{t吨} ∣ 年_{1 : t吨 - 1}) .

我们获得对(年_k个|年_1:k个−1)=𝒩(年_k个;年_{k个|k个−1},S公司_k个)，其中𝒩(年_k个;年_{k个|k个−1},S公司_k个)是参数的高斯分布年_k个，平均值年_{k个|k个−1}和协方差S公司_k个.

2.2. 粒子过滤

2.2.1. 顺序重要性抽样和重抽样

我们在这里简要描述了如何应用顺序重要性抽样-重抽样（SISR）方法从对(年_1:t吨|z（z）_1:t吨); 参见Doucet等。(2001)了解更多详细信息。

时间t吨−1，假设我们已经说过N个粒子 ${{\tilde{年}}_{1 : t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ 大致按照对(年_1:t吨−1|z（z）_1:t吨−1)我们想获得N个粒子 ${{\tilde{年}}_{1 : t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ 分配依据对(年_1:t吨|z（z）_1:t吨). 时间t吨，我们“延伸”每个粒子 ${\tilde{年}}_{1 : t吨 - 1}^{(我)}$ 通过采样 ${\tilde{年}}_{t吨}^{(我)}$ 根据条件密度 ${q个}_{t吨} (年_{t吨} ∣ {\tilde{年}}_{1 : t吨 - 1}^{(我)}, {z（z）}_{1 : t吨})$ ⁠因此，每个粒子 ${\tilde{年}}_{1 : t吨}^{(我)}$ 根据对(年_1:t吨−1|z（z）_1:t吨−1)q个_t吨(年_t吨|年_1:t吨−1,z（z）_1:t吨). 更正以下内容之间的差异对(年_1:t吨−1|z（z）_1:t吨−1)×q个_t吨(年_t吨|年_1:t吨−1,z（z）_1:t吨)和对(年_1:t吨|z（z）_1:t吨)，我们使用重要性抽样，以便与密度相关的分布对(年_1:t吨|z（z）_1:t吨)近似值为

{\hat{对}}_{N个} (d日 年_{1 : t吨} ∣ {z（z）}_{1 : t吨}) = \frac{\sum_{我 = 1}^{N个} w个 ({\tilde{年}}_{1 : t吨}^{(我)}) δ_{{\tilde{年}}_{1 : t吨}^{(我)}} (d日 年_{1 : t吨})}{\sum_{我 = 1}^{N个} w个 ({\tilde{年}}_{1 : t吨}^{(我)})} = \sum_{我 = 1}^{N个} {w个}_{t吨}^{(我)} δ_{{\tilde{年}}_{1 : t吨}^{(我)}} (d日 年_{1 : t吨}),

（6）

其中，使用公式(4)，我们有重要的权重

w个 (年_{1 : t吨}) ⑪ \frac{对 (年_{1 : t吨} ∣ {z（z）}_{1 : t吨})}{对 (年_{1 : t吨 - 1} ∣ {z（z）}_{1 : t吨 - 1}) {q个}_{t吨} (年_{t吨} ∣ 年_{1 : t吨 - 1}, {z（z）}_{1 : t吨})} ⑪ \frac{对 ({z（z）}_{t吨} ∣ 年_{t吨}) 对 (年_{t吨} ∣ 年_{1 : t吨 - 1})}{{q个}_{t吨} (年_{t吨} ∣ 年_{1 : t吨 - 1}, {z（z）}_{1 : t吨})} .

算法的性能取决于重要性密度q个_t吨(年_t吨|年_1:t吨−1,z（z）_1:t吨). 我们可以选择对(年_t吨|年_1:t吨−1)因为它是高斯密度。在这种情况下，相关的重要性权重等于w个(年_1:t吨)∝对(z（z）_t吨|年_t吨). 请注意，“最佳”重要性密度，即最小化权重条件方差的密度年_1:t吨−1（Doucet等。,2000)，是

对 (年_{t吨} ∣ 年_{1 : t吨 - 1}, {z（z）}_{1 : t吨}) ⑪ 对 ({z（z）}_{t吨} ∣ 年_{t吨}) 对 (年_{t吨} ∣ 年_{1 : t吨 - 1}),

关联的重要性权重为

w个 (年_{1 : t吨}) ⑪ 对 ({z（z）}_{t吨} ∣ 年_{1 : t吨 - 1}) = {¦Β}^{​} 对 ({z（z）}_{t吨} ∣ 年_{t吨}) 对 (年_{t吨} ∣ 年_{1 : t吨 - 1}) d日 年_{t吨} .

(7)

根据以下条件，可能计算或不计算该重量对(z（z）_t吨|年_t吨).

最后，我们得到N个粒子 ${年_{1 : t吨}^{(我)}}_{我 = 1}^{N个}$ 近似分布依据对(年_1:t吨|z（z）_1:t吨)通过从等式中给出的加权经验分布重采样(6). 文献中提供了几种重采样程序。我们在这里采用了Kitagawa中描述的分层抽样方案(1996).

替代SMC方法可用于估算对(年_1:t吨|z（z）_1:t吨). 特别是Pitt和Shephard的辅助粒子滤波（APF）技术(1999)可以使用。APF背后的理念是扩展现有粒子 ${\tilde{年}}_{1 : t吨 - 1}^{(我)}$ 从预测可能性的意义上来说，这是最有希望的 $对 ({z（z）}_{t吨} ∣ {\tilde{年}}_{1 : t吨 - 1}^{(我)})$ 很大。何时 $对 ({z（z）}_{t吨} ∣ {\tilde{年}}_{1 : t吨 - 1}^{(我)})$ APF提出了一个分析近似值。在这种情况下，APF和SISR有显著差异。然而，当 $对 ({z（z）}_{t吨} ∣ {\tilde{年}}_{1 : t吨 - 1}^{(我)})$ 可以解析计算，然后APF使用最佳重要性密度。这被称为“完美适应”（皮特和谢泼德，1999). 在这种特殊情况下，APF和SISR基本上相似，不同之处在于APF颠倒了采样和重采样步骤的顺序；这是可能的，因为重要性权重与年_t吨这显然效率更高。

2.2.2. 算法

我们的演示仅限于重要密度的标准选择，其中q个_t吨(年_t吨|年_1:t吨−1,z（z）_1:t吨)取决于(年_1:t吨−1,z（z）_1:t吨)仅通过z（z）_t吨和低维充分统计集x个_{t吨|t吨−1}和P（P）_{t吨|t吨−1}.我们将写信

{q个}_{t吨} (年_{t吨} ∣ {x个}_{t吨 ∣ t吨 - 1}, {P（P）}_{t吨 ∣ t吨 - 1}, {z（z）}_{t吨}) ≜ {q个}_{t吨} (年_{t吨} ∣ 年_{1 : t吨 - 1}, {z（z）}_{1 : t吨})

这类密度包括

对 (年_{t吨} ∣ {x个}_{t吨 ∣ t吨 - 1}, {P（P）}_{t吨 ∣ t吨 - 1}) ≜ 对 (年_{t吨} ∣ 年_{1 : t吨 - 1}) = N个 (年_{t吨}; 年_{t吨 ∣ t吨 - 1}, {S公司}_{t吨}),

哪里年_{t吨|t吨−1}和S公司_t吨是的确定函数x个_{t吨|t吨−1}和P（P）_{t吨|t吨−1}作为一个典型的关注边缘特征的人对(年_t吨|z（z）_1:t吨)，只需要存储在内存中 ${年_{t吨}^{(我)}, {x个}_{t吨 ∣ t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ 和P（P）_{t吨|t吨−1}而不是 ${年_{1 : t吨}^{(我)}}_{我 = 1}^{N个}$ ⁠与Chen和Liu提出的算法相反(2000)和Doucet等。(2000)，我们指出我们不必计算N个“完整”卡尔曼滤波递归，因为大多数计算只需进行一次。更准确地说，我们注意到 ${P（P）}_{t吨 ∣ t吨 - 1}^{(我)} = {P（P）}_{t吨 ∣ t吨 - 1}$ 和 ${S公司}_{t吨 ∣ t吨}^{(我)} = {S公司}_{t吨] t吨}$ 对于任何我∈ {1, … ,N个}.

鉴于N个粒子 ${年_{t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ 时间t吨−1大约根据对(年_t吨−1|z（z）_1:t吨−1)以及相关的足够统计数据 ${{x个}_{t吨 ∣ t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ 和P（P）_{t吨|t吨−1}，粒子过滤器在以下时间进行t吨。在顺序重要性抽样步骤:

（a）
对于我=1, … ,N个，套
${\tilde{x个}}_{t吨 ∣ t吨 - 1}^{(我)} ≜ {x个}_{t吨 ∣ t吨 - 1}^{(我)}$
和样品
${\tilde{年}}_{t吨}^{(我)} ~ {q个}_{t吨} (年_{t吨} ∣ {\tilde{x个}}_{t吨 ∣ t吨 - 1}^{(我)}, {P（P）}_{t吨 ∣ t吨 - 1}, {z（z）}_{t吨});$
（b）
对于我=1, … ,N个、评估并规范重要性权重
${w个}_{t吨}^{(我)} ⑪ \frac{对 ({z（z）}_{t吨} ∣ {\tilde{年}}_{t吨}^{(我)}) 对 ({\tilde{年}}_{t吨}^{(我)} ∣ {\tilde{x个}}_{t吨 ∣ t吨 - 1}^{(我)}, {P（P）}_{t吨 ∣ t吨 - 1})}{{q个}_{t吨} ({\tilde{年}}_{t吨}^{(我)} ∣ {\tilde{x个}}_{t吨 ∣ t吨 - 1}^{(我)}, {P（P）}_{t吨 ∣ t吨 - 1}, {z（z）}_{t吨})}, \sum_{我 = 1}^{N个} {w个}_{t吨}^{(我)} = 1$
(8)

在重采样步骤：繁殖或丢弃粒子 ${{\tilde{年}}_{t吨, 我}^{(我)}, {\tilde{x个}}_{t吨 ∣ {t吨}_{- 1}}^{(我)}}_{我 = 1}^{N个}$ 关于高或低重要性权重 ${w个}_{t吨}^{(我)}$ 以获得N个粒子 ${年_{t吨}^{(我)}, {x个}_{t吨 ∣ t吨 - 1}^{(我)}}_{我 = 1}^{N个}$ ⁠.

在更新步骤:

（a）
计算P（P）_t吨+1|t吨鉴于P（P）_{t吨|t吨−1}使用卡尔曼递归的一步（5）；
（b）
对于我=1, … ,N个，使用卡尔曼递归（5）的一个步骤计算 ${x个}_{t吨 + 1 ∣ t吨}^{(我)}$ 鉴于 $年_{t吨}^{(我)}, {x个}_{t吨 ∣ t吨 - 1}^{(我)}$ 和P（P）_{t吨|t吨−1}.

该算法在每个时间步长的计算复杂度为𝒪(N个). 如果方程给出的（未规范化的）重要性权重(8)是上界的，则渐近收敛(N个→ ∞ ) 可以确保蒙特卡罗估计接近其真实值（Crisan，2001).

2.3. 扩展

模型和算法都有许多潜在的扩展。

2.3.1. 模型

对于线性高斯模型（1）–（2），我们可以很容易地考虑以下情况v（v）_t吨和w个_t吨相关和/或添加非线性项ϕ(年_1:t吨−1)到等式的右边(2). 也可以应用上述边缘化方法进行分析整合(x个_t吨)当模型（1）–（2）不是线性高斯，而是Shephard中描述的条件线性高斯时(1994); 这种扩展允许我们考虑高斯分布的有限或连续混合。

另一个有趣的扩展包括部分观察到的隐马尔可夫模型：(x个_t吨)被建模为有限状态空间马尔可夫链

对 (年_{1 : t吨}, {z（z）}_{1 : t吨} ∣ {x个}_{1 : t吨}) = {\prod^{​}}_{k个 = 1}^{t吨} 对 (年_{k个} ∣ {x个}_{k个}, 年_{1 : k个 - 1}) 对 ({z（z）}_{k个} ∣ 年_{k个}) .

我们可以整合(x个_t吨)并计算对(年_k个|年_1:k个−1)使用隐马尔可夫模型滤波器代替卡尔曼滤波器，作为上一节中开发的方法的一部分。

2.3.2. 算法

当重要性权重的分布 ${w个}_{t吨}^{(我)}$ 倾斜，粒子 ${{\tilde{年}}_{1 : t吨}^{(我)}}_{我 = 1}^{N个}$ 多次选择重要性较高的权重；这将导致样品“耗尽”为大量颗粒 $年_{1 : t吨}^{(我)} 和年_{1 : t吨}^{(j个)}$ 实际上等于我≠j个为了执行样本“再生”，可以使用最近提出的基于马尔可夫链蒙特卡罗步骤的方法（Gilks和Berzuini，2001). 它包括应用于每个粒子 $年_{1 : t吨}^{(我)}$ 一个（可能是非遍历的）转换核 ${K（K）}_{t吨} (年_{1 : t吨} ∣ 年_{1 : t吨}^{(我)})$ 密度不变对(年_1:t吨|z（z）_1:t吨). 这个内核有无数种可能的选择。一种可能是及时更新t吨价值观年_{t吨−M（M）+1:t吨}(M（M）>0）使用高效扫描取样器（de Jong，1997)，其计算复杂度为𝒪(M（M）). 虽然这一步对于确保算法的理论收敛性来说不是必需的，但它可以改善结果。

3.模拟

3.1. 动态tobit模型

让我们考虑以下tobit模型（Manrique和Shephard，1998):

\begin{matrix} {x个}_{t吨 + 1} = ϕ {x个}_{t吨} + σ_{v（v）} {v（v）}_{t吨 + 1}, {x个}_{0} \sim N个 {0, σ_{v（v）}^{2} / (1 - ϕ^{2})}, {v（v）}_{t吨} \overset{IID公司}{\sim} N个 (0, 1), \\ 年_{t吨} = {x个}_{t吨} + σ_{ε} ε_{t吨}, ε_{t吨} \overset{IID公司}{\sim} N个 (0, 1), \\ {z（z）}_{t吨} = 最大值 (年_{t吨}, 0) . \end{matrix}

很明显，该模型的形式为（1）-（3）。我们选择“最佳”密度作为重要性密度对(年_t吨|年_1:t吨−1,z（z）_t吨). （我们无法使用对(年_t吨|年_1:t吨−1)对于重要性密度，当z（z）_t吨=0.在这种情况下，重要性权重不存在。）如果z（z）_t吨>0，那么年_t吨=z（z）_t吨如果z（z）_t吨=0，则

对于重要性权重，我们使用以下公式获得(7)

w个 (年_{1 : t吨}) ⑪ {\begin{array}{l} Φ (- 年_{t吨 ∣ t吨 - 1} / \sqrt {S公司}_{t吨}) \\ N个 ({z（z）}_{t吨}; 年_{t吨 ∣ t吨 - 1}, {S公司}_{t吨}) \end{array} \begin{matrix} 如果 {z（z）}_{t吨} = 0, \\ 如果 {z（z）}_{t吨} > 0, \end{matrix}

其中Φ（·）是标准正态分布的累积函数。

我们模拟了T型=200个具有已知超参数的观测值φ=0.99, $σ_{v（v）}^{2} = 0.05$ 和 $σ_{w个}^{2} = 0.30$ ⁠.对于不同数量的颗粒N个，我们生成了K（K）=100种不同的Rao–Blackwellized（RB）滤波器实现和标准算法估计对(x个_0:t吨,年_1:t吨|z（z）_1:t吨)使用重要性密度对(x个_t吨,年_t吨|x个_t吨−1,年_t吨−1,z（z）_t吨)作为对(x个_t吨,年_t吨|x个_t吨−1,年_t吨−1)无法使用。我们的比较是根据平方误差SE的平均值和方差进行的，计算如下：

\begin{matrix} 东南方 (我) = \sum_{t吨 = 1}^{T型} {{x个}_{t吨} - {\hat{E类}}_{N个, 我} ({x个}_{t吨} ∣ {z（z）}_{1 : t吨})}^{2}, \\ 米 (东南方) = \frac{1}{K（K）} \sum_{我 = 1}^{K（K）} 东南方 (我), \\ σ^{2} (东南方) = \frac{1}{K（K）} \sum_{我 = 1}^{K（K）} {东南方 (我) - 米 (东南方)}^{2} \end{matrix}

哪里 ${\hat{E类}}_{N个, 我} ({x个}_{t吨} ∣ {z（z）}_{1 : t吨})$ 使用我粒子滤波的实现。我们在表格中列出1和2标准和RB滤波器的性能。对于固定数量的粒子，RB过滤器的计算强度不比标准过滤器高，并且性能显著更好。

表1

米（SE）用于标准过滤器和RB过滤器

算法	以下N值的结果：
	100	250	500	1000	2500	5000	10000	25000
标准过滤器	33.70	33.64	33.90	33.41	33.45	33.55	33.54	33.52
RB过滤器	33.52	33.49	33.51	33.50	33.50	33.49	33.51	33.50

新标签中打开

表1

米（SE）用于标准过滤器和RB过滤器

算法	以下N值的结果：
	100	250	500	1000	2500	5000	10000	25000
标准过滤器	33.70	33.64	33.90	33.41	33.45	33.55	33.54	33.52
RB过滤器	33.52	33.49	33.51	33.50	33.50	33.49	33.51	33.50

新标签中打开

表2

10σ（SE）用于标准过滤器和RB过滤器

算法	以下N值的结果：
	100	250	500	1000	2500	5000	10000	25000
标准过滤器	3.76	2.19	1.61	1.20	0.83	0.62	0.37	0.26
RB滤波器	2.50	1.30	1.21	0.99	0.51	0.38	0.29	0.14

新标签中打开

表2

10σ（SE）用于标准过滤器和RB过滤器

算法	以下N值的结果：
	100	250	500	1000	2500	5000	10000	25000
标准过滤器	3.76	2.19	1.61	1.20	0.83	0.62	0.37	0.26
RB过滤器	2.50	1.30	1.21	0.99	0.51	0.38	0.29	0.14

新标签中打开

3.2. 动态probit模型

我们在这里分析了一个非平稳的二进制时间序列，更具体地说，分析了（聚合的）东京降雨量数据集（Knorr-Held，1999). 它包括T型=366个观察值z（z）_t吨=1表示下雨t吨年的第天和z（z）_t吨否则=0。我们建模z（z）_t吨通过使用动态probit模型，即。

公共关系 ({z（z）}_{t吨} = 1 ∣ α_{t吨}) = Φ (α_{t吨}),

其中(α_t吨)通过使用二阶随机行走进行建模

α_{t吨} = 2 α_{t吨 - 1} - α_{t吨 - 2} + σ_{v（v）} {v（v）}_{t吨}, {v（v）}_{t吨} \overset{IID公司}{\sim} N个 (0, 1) .

(9)

让我们介绍一个人工的潜在过程(年_t吨)这样的话

年_{t吨} = α_{t吨} + ε_{t吨}, ε_{t吨} \overset{IID公司}{\sim} N个 (0, 1),

(10)

并定义

(11)

很容易检查我们是否有

公共关系 ({z（z）}_{t吨} = 1 ∣ α_{t吨}) = 公共关系 (年_{t吨} > 0 ∣ α_{t吨}) = 公共关系 (ε_{t吨} > - α_{t吨}) = Φ (α_{t吨}) .

我们可以很容易地改写方程式(9)–(11)在形式为（1）-（3）的状态空间模型中，通过定义

{x个}_{t吨} ≜ (α_{t吨}, α_{t吨 - 1}) .

在这种情况下，作为(年_t吨)是人为的，有必要将我们的程序与应用于估计对(x个_0:t吨|z（z）_1:t吨)而不是与对(x个_0:t吨,年_1:t吨|z（z）_1:t吨). 引入的动机(年_t吨)这是因为可以使用最佳密度作为重要性密度

这是一个截断高斯分布，并且，使用公式(7)，我们获得

w个 (年_{1 : t吨}) ⑪ 对 ({z（z）}_{t吨} ∣ 年_{1 : t吨 - 1}) = {1 - Φ (- \frac{年_{t吨 ∣ t吨 - 1}}{\sqrt {S公司}_{t吨}})}^{{z（z）}_{t吨}} Φ {(- \frac{年_{t吨 ∣ t吨 - 1}}{\sqrt {S公司}_{t吨}})}^{1 - {z（z）}_{t吨}} .

如果我们考虑对(x个_0:t吨|z（z）_1:t吨)，由于相关的重要性权重，因此无法使用最佳密度对(z（z）_t吨|x个_t吨−1)不接受分析表达式。请注意(年_t吨)已经用于开发高效的马尔可夫链蒙特卡罗采样器（Albert和Chib，1993).

超参数 $σ_{v（v）}^{2}$ 设置为0.01。在图中。1，我们显示观察结果z（z）_t吨,E类{Φ(α_t吨)|年_1:t吨}和E类｛Φ(α_t吨)|年_1:t吨}±√变量{Φ(α_t吨)|年_1:t吨}；估计值是通过使用N个=1000个粒子。为了获得类似的结果，引导过滤器（Gordon等。,1993)需要多达N个=5000个粒子。

图1

（a）二进制观测值（zt）和（b）E{Φ（αt）|y1:t}（--）和E{φ（αt

新标签中打开下载幻灯片

（a）二进制观测(z（z）_t吨)和（b）E类{Φ(α_t吨)|年_1:t吨} (——)和E类{Φ(α_t吨)|年_1:t吨}±√变量{Φ(α_t吨)|年_1:t吨} (——)

4.结论

本文提出了一种部分观测高斯模型的递推状态估计方法。我们的算法是一种基于边缘化的SMC方法。这种边缘化是通过卡尔曼滤波方法实现的。我们的仿真表明，我们的方法可以显著优于标准SMC方法。

在本文中，模型参数λ假设已知。然而，结合此处开发的粒子滤波方法和Liu和West最近提出的方法，可以对这些参数进行批量和递归估计(2001).

致谢

我们感谢Leonhard Knorr-Held提供东京降雨量数据集，感谢Sanjeev Arulampalam、Nando de Freitas、Arthur Gretton和Elena Punskaya提出宝贵意见。

参考文献

1

艾伯特

,

J。

和

芯片

,

美国。

(

1993

)

二进制和多光子响应数据的贝叶斯分析

.

《美国统计杂志》。助理。

,

88

,

669

——

679

.

谷歌学者

交叉参考

书目数据库

2

陈

,

R。

和

线路接口单元

,

J.S.公司。

(

2000

)

混合卡尔曼滤波器

.

J.R.统计。Soc公司

.B中，

62

,

493

——

508

.

谷歌学者

交叉参考

书目数据库

三

克里桑

,

D。

(

2001

)颗粒过滤器——理论观点。在

序贯蒙特卡罗方法在实践中的应用

（编辑

答：。

水龙头

,

J·F·G。

德弗里塔斯

和

新泽西州。

戈登

)，第页。

17

——

41

纽约：

施普林格

.

4

水龙头

,

答：。

,

De Freitas公司

,

J·F·G。

和

戈登

,

新泽西州。

（编辑）(

2001

)

序列蒙特卡罗方法在实践中的应用

纽约：

施普林格

.

5

水龙头

,

答：。

,

戈德斯基

,

S.J.公司。

和

安德烈厄

,

C、。

(

2000

)

贝叶斯滤波的序贯蒙特卡罗抽样方法

.

统计师。计算。

,

10

,

197

——

208

.

谷歌学者

交叉参考

书目数据库

6

盖尔芬德

,

答：。

和

史密斯

,

A.F.M.公司。

(

1990

)

基于抽样的边际密度计算方法

.

《美国统计杂志》。助理。

,

85

,

398

——

409

.

谷歌学者

交叉参考

书目数据库

7

吉尔克斯

,

W.R.公司。

和

贝尔祖伊尼

,

C、。

(

2001

)

动态贝叶斯模型的跟踪运动目标蒙特卡罗推断

.

J.R.统计。Soc公司

.B、，

63

,

127

——

146

.

谷歌学者

交叉参考

书目数据库

8

戈登

,

新泽西州。

,

鲑鱼

,

D.J.博士。

和

史密斯

,

A.F.M.公司。

(

1993

)

非线性/非高斯贝叶斯状态估计的新方法

.

IEE程序

.F、，

140

,

107

——

113

.

谷歌学者

OpenURL占位符文本

书目数据库

9

德容

,

第页。

(

1997

)

扫描取样器

.

生物特征

,

84

,

929

——

937

.

谷歌学者

交叉参考

书目数据库

10

北川

,

G.公司。

(

1996

)

非高斯非线性状态空间模型的蒙特卡罗滤波和平滑器

.

J.计算。图表。统计师。

,

5

,

1

——

25

.

谷歌学者

OpenURL占位符文本

书目数据库

11

Knorr-Held公司

,

L。

(

1999

)

动态模型中的条件优先建议

.

扫描。J.统计学家。

,

26

,

129

——

144

.

谷歌学者

交叉参考

书目数据库

12

线路接口单元

,

J。

和

西部

,

M。

(

2001

)基于仿真的滤波中的参数和状态联合估计

序贯蒙特卡罗方法在实践中的应用

（编辑

答：。

水龙头

,

J·F·G。

德弗里塔斯

和

新泽西州。

戈登

)，第页。

197

——

223

纽约：

施普林格

.

13

曼里克

,

答：。

和

谢泼德

,

N。

(

1998

)

有限相依过程的模拟似然推理

.

计量经济学。J。

,

1

,

174

——

202

.

谷歌学者

交叉参考

书目数据库

14

皮特

,

M.K.博士。

和

谢泼德

,

N。

(

1999

)

通过模拟过滤：辅助粒子过滤器

.

《美国统计杂志》。助理。

,

94

,

590

——

599

.

谷歌学者

交叉参考

书目数据库

15

谢泼德

,

N。

(

1994

)

部分非高斯时间序列模型

.

生物特征

,

81

,

115

——

131

.

谷歌学者

交叉参考

书目数据库

16

西部

,

M。

和

哈里森

,

P.J.公司。

(

1997

)

贝叶斯预测和动态模型

，第2版。纽约：

施普林格

.

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

下载所有幻灯片

月份：	总浏览次数：
2023年4月	11
2023年5月	6
2023年6月	7
2023年7月	10
2023年8月	22
2023年9月	10
2023年10月	9
2023年11月	7
2023年12月	32
2024年1月	18
2024年2月	10
2024年3月	13
2024年4月	47
2024年5月	11

文章内容

部分观测高斯状态空间模型的粒子滤波

总结

1.简介

1.1. 背景

1.2. 一般问题

1.3. 分辨率

1.4. 计划

2.Rao–Blackwellized粒子过滤

2.1. 边缘化

2.2. 粒子过滤

2.2.1. 顺序重要性抽样和重抽样

2.2.2. 算法

2.3. 扩展

2.3.1. 模型

2.3.2. 算法

3.模拟

3.1. 动态tobit模型

3.2. 动态probit模型

4.结论

致谢

参考文献

引文

视图

海拔高度

电子邮件警报

通过引用文章

最新的

阅读次数最多

引用次数最多

文章内容

部分观测高斯状态空间模型的粒子滤波

总结

1.简介

1.1. 背景

1.2. 一般问题

1.3. 分辨率

1.4. 计划

2.Rao–Blackwellized粒子过滤

2.1. 边缘化

2.2. 粒子过滤

2.2.1. 顺序重要性抽样和重抽样

2.2.2. 算法

2.3. 扩展

2.3.1. 模型

2.3.2. 算法

3.模拟

3.1. 动态tobit模型

3.2. 动态probit模型

4.结论

致谢

参考文献

引文

视图

海拔高度

电子邮件警报

通过引用文章

最新的

阅读次数最多

引用次数最多

此功能仅对订阅服务器可用