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的过程可以被如下归纳：

1. 初始化$N$个粒子：$\{x^{(i)}_0=p(x_0),i\in N\}$，令$t=1$。
2. 采样$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)})$。
3. 估计每个粒子的权重：$w_t^{(i)}=p(y_t|\hat x_t^{(i)})$ 并归一化得到$\tilde w_t^{(i)}$（由于上一步的重采样使$w_{t-1}^{(i)}$都相同，所以不需要参与权重计算）。
4. 根据每个粒子的权重$\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}$。
5. 重复2-4步骤，且$t=t+1$

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

Conclusion

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

Reference

https://www.stats.ox.ac.uk/~doucet/doucet_defreitas_gordon_smcbookintro.pdf