Particle filter 推导理解——第二部分

前情回顾

上一部分主要讲述了隐马尔可夫模型、滤波与边沿似然。
状态估计模型可以抽象为隐马尔可夫模型，因为系统真实值未知，观测值可知，并且能够假设系统当前状态仅与前一时刻的状态有关。以HMM（隐马尔可夫）模型为基础，将系统估计抽象成为一个贝叶斯模型。并且得到了系统的先验概率和似然函数（若有疑问可以查看Particle filter 推导与理解——第一部分）。后验估计根据贝叶斯模型进行展开。
滤波与边沿似然基于马尔科夫模型确定的系统状态和观测的公式，和贝叶斯公式结合，得到了系统第$n$时刻的后验分布和第$n-1$时刻的后验分布之间的关系，滤波的目标为“当前的状态的估计”，所以将第1到$n$时刻的状态值的联合分布$p(x_{1:n}|y_{1:n})$转化成了第$n$时刻的边沿分布$p(x_n|y_{1:n})$。但是教程中说明了一下，并不是所有的Particle filter都是这样表示的。

本部分内容简介

介绍蒙特卡罗采样思想是如何作用于粒子滤波原理的，以及粒子滤波算法的构成。主要包括蒙特卡罗样本估计简介、重要性采样、顺序重要采样、重采样。

蒙特卡罗采样估计

什么是蒙特卡罗采样估计？
学过概率论的童鞋对“样本估计整体”这句话一定不陌生，蒙特卡罗模型在这里依旧是发挥了一个用样本估计整体的角色。
假设有一关于序列$x_{1:n}$的未知分布$${\pi }_n(x_{1:n})=\frac{{\gamma }_n(x_{1:n})}{Z_n}\text{\hspace{1em}(1)}$$
其中${\gamma }_n$是唯一我们能够得到的值（我认为可以理解为一个符合${\pi }_n$分布的实例），$Z_n$是归一化常数：
$$Z_n={\int }_{}{\gamma }_n(x_{1:n})d{x}_{1:n}\text{\hspace{1em}(2)}$$
蒙特卡罗估计就是从分布${\pi }_n(x_{1:n})$中抽取$N$个独立同分布的样本${X_n}^i,i=1,...,N$，抽取的目的是用样本估计整体的情况。则通过蒙特卡罗估计得：
$$\hat{\pi }(x_{1:n})=\frac{1}{N}\sum _{i=1}^N{\delta }_{{x_{1:n}}^i}(x_{1:n})\text{\hspace{1em}(3)}$$
以式对于任意一个边沿分布的估计有：
$$\hat{\pi }(x_k)=\frac{1}{N}\sum _{i=1}^N{\delta }_{{x_k}^i}(x_k)\text{\hspace{1em}(4)}$$
其中$\delta$表示delta函数。
上面几个公式表明了蒙特卡罗采样估计的一般公式表达。
需要说明一点：蒙特卡罗估计是对参数的无偏估计——也就是估计值的方差随着采样点数的增大最终收敛到0（也就是估计值最终收敛到真实值）。

为什么要提到蒙特卡罗采样估计？
我们在上一部分已经得到了关于此状态估计模型的后验分布表达式：
$$p(x_n|y_{1:n})=\frac{g(y_n|x_n)p(x_n|y_{1:n-1})}{p(y_n|y_{1:n-1})}\text{\hspace{1em}(5)}$$
我们希望对该后验分布进行估计，从而达到滤波的目的。那么怎么对它进行估计呢?蒙特卡罗告诉我们只要在该后验分布上采样，用样本估计整体就可以了。但是，在该状态估计模型中，后验分布是未知的，无法直接进行采样。需要说明的是，在单纯的蒙特卡罗估计中，存在高维度样本采样困难的问题，同时也存在计算复杂度大的问题。蒙特卡罗模型通常使用重要性采样解决采样困难的问题。同理，对于我们目前遇到的后验概率分布未知的难题，同样可以使用重要性采样的方式进行解决。

什么是重要性采样？
我们状态估计遇到的问题是：需要估计的值的概率分布未知。重要性采样的思想是：针对概率未分布知的情况，我们建立一个容易采样的概率分布模型（例如，高斯分布，均匀分布；选取的分布类型会对估计产生一定影响，但是如何影响的我还没有仿真验证），并且对该概率分布模型进行采样，通过对样本设置不同的权重，来近似未知分布的采样样本，从而达到对未知分布进行估计的目的。
即，对于未知的${\pi }_n(x_{1:n})$，有重要性分布$q_n(x_{1:n})$满足${\pi }_n(x_{1:n})>0=>q_n(x_{1:n})>0$
所以，
$${\pi }_n(x_{1:n})=\frac{w_n(x_{1:n})q_n(x_{1:n})}{Z_n}\text{\hspace{1em}(6)}$$
$$Z_n={\int }_{}{w}_n(x_{1:n})q_n(x_{1:n})d{x}_{1:n}\text{\hspace{1em}(7)}$$
其中$w_n(x_{1:n})$是未归一化的权重函数$w_n(x_{1:n})=\frac{{\gamma }_n(x_{1:n})}{q_n(x_{1:n})}$

为什么重要性采样？
这个问题显而易见。我们的状态估计模型后验分布未知，自然是希望通过重要性采样的方式估计后验概率。在particle filter 中重要性采样的意义在于，将无法采样的后验分布问题转化成了一个对"容易采样的分布的样本的权值分布"问题。权值的计算依赖于对应的观测值——每个样本对应的观测值与系统的观测值比较，距离小的分配较大的权值。
但是，上述重要性采样面临一个问题：得到一个测试函数的最优估计比较容易，但是我们对于一个特定的测试函数不感兴趣，而对一系列测试函数的均值感兴趣。即便是对一个特定的函数感兴趣，第$n$时刻的最优估计对于我们查看整个过程的状态是有害的（意思是在$n-1$时刻得到的最优估计${\hat{\varphi }}_{n-1}(x_{1:n-1})$与在$n$时刻的估计值${\hat{\varphi }}_n(x_{1:n})$中包含的$x_{1:n-1}$的边沿估计不相似——前一时刻的最优估计和后一时刻的最优估计没什么联系）。

顺序重要性采样（sequential importance sampling）

SIS是IS的一种，重要性分布满足：
$$q_n(x_{1:n})=q_{n-1}(x_{1:n-1})q_n(x_n|x_{1:n-1})=q_1(x_1)\prod _{k=2}^n{q}_k(x_k|x_{k-1})\text{\hspace{1em}(8)}$$
$$w_n(x_{1:n})=\frac{{\gamma }_n(x_{1:n})}{q_n(x_{1:n})}=\frac{{\gamma }_{n-1}(x_{1:n-1})}{q_{n-1}(x_{1:n-1})}\frac{{\gamma }_n(x_{1:n})}{{\gamma }_{n-1}(x_{1:n-1})q_n(x_n|x_{1:n-1})}\text{\hspace{1em}(9)}$$
将式（9）简写为：
$$w_n(x_{1:n})=w_{n-1}(x_{1:n-1}){\alpha }_n(x_{1:n})\text{\hspace{1em}(10)}$$
这样就提供了一个迭代方式，解决了计算量的问题。但是，利用SIS对系统状态进行估计仍然面临估计方差随着时间的推移产生指数性增长的事实。不加证明地说明，增加采样个数减小估计方差的方式产生的代价很大。
针对上述的问题，提出了重采样

什么是重采样？

例如，IS方法在$n$时刻对未知分布${\pi }_n(x_{1:n})$，通过对重要性分布$q_n(x_{1:n})$加权估计的结果为${\hat{\pi }}_n(x_{1:n})$，下一时刻我们在分布${\hat{\pi }}_n(x_{1:n})$上依据$q_n(x_{1:n})$的样本的权值$W_n^i$采样。重采样之后各个样本的权值一致。在particle filter 迭代中：计算重要性权值–>重采样–>计算重要性权值……。

上述知识如何应用于Particle filter？

Particle filter 的目的是对$p(x_n|y_{1:n})$进行估计——通过对系统的观测，估计系统的状态。而蒙特卡罗通过采样的方式进行估计整体，所以希望通过对系统后验$p(x_n|y_{1:n})$采样，估计系统的真实状态。我们遇到的问题是：后验分布$p(x_n|y_{1:n})$未知。所以要通过重要性采样的方式解决后验分布未知的难题。

假设${\gamma }_n(x_{1:n})=p(x_{1:n},y_{1:n})$，${\pi }_n(x_{1:n})=p(x_{1:n}|y_{1:n})$，$Z_n=p(y_{1:n})$。利用蒙特卡罗估计唯一需要做的就是选择重要性分布$q_n(x_n|x_{1:n-1})$且$q_n^{opt}(x_n|x_{1:n-1})={\pi }_n(x_n|x_{1:n-1})$（注意这是理想的情况，是不能达到的，我们能做的就是尽量逼近上述的结果）且有
$${\pi }_n(x_n|x_{1:n-1})=\frac{g(y_n|x_n)f(x_n|x_{n-1})}{p(y_n|x_{n-1})}\text{\hspace{1em}(11)}$$
式（11）的推导暂时有点疑问。
在上述的理想情况下结合${\alpha }_n(x_{1:n})=\frac{{\gamma }_n(x_{1:n})}{{\gamma }_{n-1}(x_{1:n-1})q_n(x_n|x_{1:n-1})}$ ,${\gamma }_n(x_{1:n})=p(x_{1:n},y_{1:n})$得到：
$${\alpha }_n(x_{1:n})=p(y_n|x_{n-1})\text{\hspace{1em}(12)}$$
根据已有经验，通常选择$q_n(x_n|x_{1:n-1})=q(x_n|y_n,x_{n-1})$也就是说这种重要性采样方式致使我们得不到任何关于$y_{1:n-1}$和$x_{1:n-2}$的任何信息。所以:
$${\alpha }_n(x_{1:n})={\alpha }_n(x_{n-1:n})=\frac{g(y_n|x_n)f(x_n|x_{n-1})}{q(x_n|y_n,x_{n-1})}\text{\hspace{1em}(13)}$$
其中，$q(x_n|y_n,x_{n-1})$是对$p(x_n|y_n,x_{n-1})$的近似。
最终我们在时间$n$得到系统状态的估计值：
$$\hat{p}(x_{1:n}|y_{1:n})=\sum _{i=1}^n{W}_n^i{\delta }_{x_{1:n}^i}(x_{1:n})\text{\hspace{1em}(14)}$$

算法截图如下：

总结起来particle filter 通过蒙特卡罗样本估计整体的思想为指导，通过重要性采样和重采样迭代逼近后验分布，最终达到滤波目的。