Sequential Monte Carlo Methods (SMC) 序列蒙特卡洛/粒子滤波/Bootstrap Filtering

Problem Statement

我们考虑一个具有马尔可夫性质、非线性、非高斯的状态空间模型(State Space Model)：对于一个时间序列上的观测结果 $\{ y_t , t \in N \}$ ，我们认为每个观测结果 $y_t$ 的生成依赖于一个无法直接观察的隐变量 $x_t \in \{x_t , t \in N \}$ ，即： $p(y_t|x_{t}), t \in [1,N]$ ；我们假设隐变量具有一个先验的状态转移函数 $p(x_t | x_{t-1}), t \in [1,N]$ ，若给定一个初始分布 $p(x_0)$ ，那么这个模型就能由这三个参数描述： $p(x_0),p(x_t| x_{t-1}),p(y_t|x_t)$ 。这里，我们定义到时间 $t$ 的观察序列： $\mathbf{x_{0:t}}=\{ x_0, \dots, x_t\}$ 和对应的隐变量序列 $\mathbf{y_{1:t}}=\{ y_1, \dots, y_t\}$ 。

我们关心如何根据当前的观测序列来推断(infer)隐变量序列，即估计一个后验概率分布 $p(\mathbf{x_{0:t}| y_{1:t}})$ ，和它的边缘概率分布： $p(x_{t}| \mathbf{y_{1:t}})$ (这里通常被称为滤波, filtering) ，以及它对于某个函数 $f_t$ 的期望：
$I(f_t) = \mathbb{E}_{p(\mathbf{x_{0:t}| y_{1:t}})} f_t(\mathbf{x_{0;t}}) = \int p(\mathbf{x_{0:t}| y_{1:t}}) f_t(\mathbf{x_{0;t}}) d_{\mathbf{x_{1:t}}}$

在任何时间 $t$ ，我们可以推导 $p(\mathbf{x_{0:t+1}| y_{1:t+1}})$ 和 $p(\mathbf{x_{0:t}| y_{1:t}})$ 之间的关系：
$\begin{aligned} {p(\mathbf{x_{0:t+1}| y_{1:t+1}})} &= \frac {p(\mathbf{x_{0:t+1}, y_{1:t+1}})} {p(\mathbf{y_{1:t+1}})} \\ &= \frac{p (\mathbf{x_{0:t+1}},y_{t+1}|\mathbf{y}_t) p(\mathbf{y_{1:t}})}{p(\mathbf{y_{1:t+1}})} \\ &= \frac {p(\mathbf{x_{0:t}}|\mathbf{y_{1:t}})p(y_{t+1},x_{t+1}|\mathbf{x_{0:t}},\mathbf{y_{1:t}}) p(\mathbf{y_{1:t}}) }{p(\mathbf{y_{1:t}},y_{t+1})} \\ &= \frac {p(\mathbf{x_{0:t}}|\mathbf{y_{1:t}})p(x_{t+1}|\mathbf{y_{1:t}},\mathbf{x_{0:t}})p(y_{t+1}|\mathbf{x_{0:t}},\mathbf{y_{1:t}},x_{t+1}) p(\mathbf{y_{1:t}}) }{p(y_{t+1}| \mathbf{y_{1:t}}) p(\mathbf{y_{1:t}})} \\&= p(\mathbf{x_{0:t}}|\mathbf{y_{1:t}}) \frac {p(x_{t+1}|x_t) p(y_{t+1}|x_{t+1}) }{p(y_{t+1}|\mathbf{y_{1:t}})} \end{aligned}$

Monte Carlo Sampling

假设我们能够从 $p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})$ 生成N个独立同分布的随机样本（也被称为粒子particles)，那么 $f_t$ 的期望可以用以下表示：
$I_N(f_t) = \frac{1}{N} \sum_{i=1}^N f_t(\mathbf{x^{(i)}_{0:t}})$
当N足够大时，蒙特卡洛采样能够无偏估计 $f_t$ 的期望，然而我们通常是不知道先验分布 $p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})$ 的，因此采样的样本是随机采样的，这就会导致这种方法效率不会很高，尤其是当随机采样的分布与 $p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})$ 的分布相差较大时。

Importance Sampling (IS)

为了解决采样效率问题，我们引用一个重要性分布 $\pi(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})$ 来代替之前蒙特卡洛方法里的随机采样，那么采样权重(importance weight)可以表示为：
$w(\mathbf{x_{0:t}})= \frac { p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})} {\pi(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})}$

则 $f_t$ 的期望可以被写成：

$\begin{aligned} I(f_t) &= \int f_t(\mathbf{x_{0:t}}) p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}}) d_{\mathbf{x_{0:t}}} \\ &= \int \pi(\mathbf{x_{0:t}}|\mathbf{y_{0:t}}) \frac { p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})} {\pi(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})} f_t(\mathbf{x_{0:t}}) d_{\mathbf{x_{0:t}}} \\&= \mathbb{E}_{\pi(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})} [w(\mathbf{x_{0:t}})f_t(\mathbf{x_{0:t}})] \\ &= \frac{1}{N} \sum_{i=1}^N [w(\mathbf{x^{(i)}_{0:t}}) f_t(\mathbf{x^{(i)}_{0:t}})] \end{aligned}$

我们还可以推导另外一种写法：由于 $\int p(x)dx=1$ ，有 $\int \pi(x)w(x)dx=1$ ，因此 $\mathbb{E}_{\pi(x)}w(x)=1$ ，然后得到 $\frac{1}{N}w(x)=1$ ，即 $w(x)=N$ ，则上面的式子可以被这么写：

$\begin{aligned} I(f_t) &= \frac{1}{N} \sum_{i=1}^N [w_t^{(i)} f_t(\mathbf{x^{(i)}_{0:t}})] \\ & = \sum_{i=1}^N [\frac{w_t^{(i)}}{N} f_t(\mathbf{x^{(i)}_{0:t}})] \\ & = \sum_{i=1}^N [\frac{w_t^{(i)}}{\sum w_t^{(i)}} f_t(\mathbf{x^{(i)}_{0:t}})] \\ &= \sum_{i=1}^N \tilde{w_t}f_t(\mathbf{x^{(i)}_{0:t}}) \end{aligned}$
这里的 $\tilde{w_t}=\frac{w_t(x)}{\sum w_t(x)}$ 是一个归一化后的采样权重(normalized importance weight)，重要性采样IS是一个通用的蒙特卡洛方法，但是由于它的表现形式，它在时间序列上并不能进行递归估计，比如：为了计算 $p(\mathbf{x_{0:t}}|\mathbf{y_{0:t}})$ ，IS需要获得所有的 $\mathbf{y_{0:t}}$ 然后才能进行计算，然后当新的 $y_{t+1}$ 到来时，IS需要重新在整个观测序列 $\mathbf{y_{0:t+1}}$ 上来计算，时间复杂度高，因此我们接下来介绍在时间序列上的重要性采样。

Sequential Importance Sampling (SIS)

我们能够扩展IS方法，使它在 $t+1$ 时刻的计算能够使用 $t$ 时刻的计算结果。我们先将重要性函数 $\pi(\mathbf{x_{0:t}|y_{1:t}})$ 用 $t-1$ 时刻描述：

$\begin{aligned} \pi(\mathbf{x_{0:t}|y_{1:t}}) &= \pi(\mathbf{x_{0:t-1}|y_{1:t}})\pi(x_t| \mathbf{x_{0:t-1},y_{1:t}}) \\ &= \pi(\mathbf{x_{0:t-1}|y_{1:t-1}})\pi(x_t| \mathbf{x_{0:t-1},y_{1:t}}) \\ \end{aligned}$
于是我们有：
$\begin{aligned} \pi(\mathbf{x_{0:t}|y_{1:t}}) &= \pi(x_0) \prod_k^t \pi(x_k|\mathbf{x_{0:k-1},y_{1:k}}) \end{aligned}$
到这里，由于我们可以推导出 $p(\mathbf{x_{0:t}}|\mathbf{y_{1:t}})$ 和 $\pi(\mathbf{x_{0:t}}|\mathbf{y_{1:t}})$ 关于 $t-1$ 时刻的计算，则对于重要性权重 $w_t$ 我们同样有：

$\begin{aligned} w^{(i)}_t(\mathbf{x_{0:t}}) &= \frac { p(\mathbf{x_{0:t}}^{(i)}|\mathbf{y_{1:t}}^{(i)})} {\pi(\mathbf{x_{0:t}}^{(i)}|\mathbf{y_{1:t}}^{(i)})} \\ &= \frac { p(\mathbf{x_{0:t-1}}^{(i)}|\mathbf{y_{1:t-1}}^{(i)})} {\pi(\mathbf{x_{0:t-1}}^{(i)}|\mathbf{y_{1:t-1}}^{(i)})} \frac {p(x_{t}^{(i)}|x_{t-1}^{(i)}) p(y_{t}|x_{t}^{(i)}) }{p(y_{t}| \mathbf{y_{1:{t-1}}}) \pi(x_t|\mathbf{x_{0:t-1}}^{(i)}, \mathbf{y_{1:t}}) } \\ &= w^{(i)}_{t-1}(\mathbf{x_{0:t-1}}^{(i)}) \frac {p(x_{t}^{(i)}|x_{t-1}^{(i)}) p(y_{t}|x_{t}^{(i)}) }{p(y_{t}|\mathbf{y_{1:t-1}}) \pi(x_t^{(i)}|\mathbf{x_{0:t-1}}^{(i)},\mathbf{y_{1:t}}) } \end{aligned}$
对于当前 $t$ 时刻来说， $p(y_{t}|\mathbf{y_{1:t-1}})$ 是确定的，而且对于所有的采样样本 $\mathbf{x_{0:t}}^{(i)},i \in N$ 都一致，因此还我们可以这么表示 $\tilde w^{(i)}_t$ 的更新：
$\begin{aligned} \tilde w^{(i)}_t(\mathbf{x_{0:t}}) & \propto \tilde w^{(i)}_{t-1}(\mathbf{x_{0:t-1}}) \frac {p(x_{t}|x_{t-1}) p(y_{t}|x_{t}) }{\pi(x_t|\mathbf{x_{0:t-1},y_{1:t}}) } \end{aligned}$
当我们采用先验分布作为重要性分布时，我们的重要性函数变成了：
$\begin{aligned} \pi(\mathbf{x_{0:t}|y_{1:t}}) &= p(\mathbf{x_{0:t}}) =p(x_0) \prod_{k=1}^t \pi(x_k|x_{k-1}) \end{aligned}$
因此 $w$ 的更新被简化为： $\tilde w^{(i)}_t(\mathbf{x_{0:t}}) \propto \tilde w^{(i)}_{t-1}(\mathbf{x_{0:t-1}}) p (y_{t}|x_t^{(i)})$ 。SIS是一个不错的方法，但本质上就是一个带约束的重要性采样方法，但是它在高维度空间的计算并不是很有效，比如当 $t$ 非常大时。

Bootstrap Filtering = Resampling + SIS

SIS的问题是：当t增加时，重要性权重的分布的计算会越来越不准确，尤其是当时间t达到一定量时，通常只有很少的样本(粒子)是有非零权重的，而大部分样本(粒子)的权重都是0，因此会让这个采样分布无法正确描述出先验分布。为了处理这种退化(degeneracy)，我们需要加入一个额外的重采样(resampling)步骤。

重采样的目的是扩增(multiply)那些权值较高的粒子个数而减少(eliminate)权值低的粒子，具体来说，假设一个粒子具有 $\tilde w_t$ 权值，重采样会在此粒子附近再采样 $N \tilde w_t$ 个粒子，若 $\tilde w_t=0$ ，那么这个粒子会被淘汰(dead)，剩下的(serviving)粒子会进入到下一个时刻，注意重采样完后的每个粒子的权值都变成了： $\tilde w_t^{(i)}=\frac{1}{N}, i \in N$ 。

假设我们的重要性采样分布就是先验分布，那么Boostrap Filtering的过程可以被如下归纳：

初始化 $N$ 个粒子： $\{x^{(i)}_0=p(x_0),i\in N\}$ ，令 $t=1$ 。
采样 $N$ 个粒子： $\{\hat x_t^{(i)} \sim p(x_t|x_{t-1}^{(i)}),i \in N\}$ ，对于每个粒子，更新其采样轨迹： $\mathbf {\hat x_{0:t}}^{(i)} = (\mathbf {x_{0:t-1}}^{(i)},\hat x_t^{(i)})$ 。
估计每个粒子的权重： $w_t^{(i)}=p(y_t|\hat x_t^{(i)})$ 并归一化得到 $\tilde w_t^{(i)}$ （由于上一步的重采样使 $w_{t-1}^{(i)}$ 都相同，所以不需要参与权重计算）。
根据每个粒子的权重 $\tilde w_t^{(i)}$ ，从 $\{\hat x_{t}^{(i)},i \in N\}$ 中重新采样N个粒子： $\{ x_{t}^{(i)},i \in N\}$ ，更新 $\mathbf { x_{0:t}}^{(i)} = (\mathbf {x_{0:t-1}}^{(i)}, x_t^{(i)})$ ，此时 $\tilde w_t^{(i)}=\frac{1}{N}$ 。
重复2-4步骤，且 $t=t+1$

Bootstrap Filtering有几个优点：1. 易实现，计算简单；2. 扩展性好，当迁移至别的问题时，只需要更改重要性采样分布即可；3. 天然支持并行计算，4.重采样过程可以看成是一种黑匣子，输入当前的权值和粒子索引，输出每个粒子要重采样的次数，因此这种黑匣子有很多种实现方法，这样可以支持很多的复杂的模型在时间序列上进行推断(inference)

Conclusion

这种方法其实代表了一种时间序列上的推断框架，它有很多名字：粒子滤波 / Bootstrap Filtering / 序列化蒙特卡洛 / 序列化重要性采样 + 重采样，这种通用的思想可以被用到很多具体的推断问题上，特别是partially observed问题，如机器人定位，partially observed MDP (POMDPs)。借用一张图来表示Bootstrap Filtering：

在这里插入图片描述