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

学习particle filter 有一段日子了，今日进行总结讲解，宏观理解与仿真均没有太大问题，但是公式的推导和物理意义的理解还存在疏漏。借此机会总结梳理一番，本文分为上下两个部分，第一部分是关于particle filter推导的前期准备，第二部分涉及Monte Carlo 思想如何作用于particle filter。

Particle filter 宏观理解

总结起来，Particle filter 通过对系统真实值进行采样、估计。先在一个已知分布上进行采样，通过比较样本与系统观测值的“距离”确定各个样本的权值。然后通过重采样“复制”权重较大的样本，“丢弃”权值小的样本，并对样本的权值进行归一化。通过迭代，在真实值附近的样本个数数量越来越多，也就达到了particle filter 追踪的目的。

隐马尔可夫模型

什么是隐马尔可夫模型？
最主要的特点：“状态值未知，观测值可知”，且当前状态仅依赖于前一时刻的状态。推荐阅读http://www.cnblogs.com/skyme/p/4651331.html。
为什么要提到隐马尔可夫模型？
Particle filter 假设需要估计的系统是一个“状态值未知，观测值可知”的模型，且假设系统的当前状态仅依赖于上一时刻的状态。以此为基础进行公式的推导。
Particle filter 系统的状态转移模型$$f(x_n|x_{n-1})\text{\hspace{1em}(1)}$$其中$x_n$表示系统当前时刻的状态，$x_{n-1}$表示上一时刻的状态，$f$代表条件概率。Particle filter 系统的观测模型$$g(y_n|x_n)\text{\hspace{1em}(2)}$$其中$y_n$表示系统当前的观测值，是可知的。
马尔可夫模型有什么特点？
假设第一时刻系统状态$x_1$符合分布$u(x_1)$，$f(x_k|x_{k-1})$表示马尔可夫链状态转移概率，那么根据上述“状态模型和观测模型”我们能够得到一个贝叶斯模型： { X n } n > = 1 \{X_n\}_{n>=1} {Xn​}n>=1​的先验分布：
$$p(x_{1:n})=u(x_1)\prod _{k=2}^{n}f(x_k|x_{k-1})\text{\hspace{1em}(3)}$$
其中$f(x_k|x_{k-1})$表示马尔可夫链中的状态转移模型。
和似然函数：
$$p(y_{1:n}|x_{1:n})=\prod _{k=1}^{n}g(y_k|x_k)\text{\hspace{1em}(4)}$$
其中观测值 { Y n } n > = 1 \{Y_n\}_n>=1 {Yn​}n​>=1之间统计独立。（这是公式（4）可以连乘的原因）
而
$$p(x_{1:n}|y_{1:n})=\frac{p(x_{1:n})p(y_{1:n}|x_{1:n})}{p(y_{1:n})}\text{\hspace{1em}(5)}$$
其中
$$p(y_{1:n})={\int }_{}^{}p(y_{1:n},x_{1:n})d{x}_{1:n}\text{\hspace{1em}(6)}$$
到这里，我们回顾一下前面都做了些什么。我们将动态系统抽象成了一个隐马尔可夫模型，这样抽象的目的是这个系统符合HMM模型的特点。通过对系统的抽象，我们得到了系统的状态值与系统第一时刻的状态的链式关系式（3），观测值与系统状态的关系如式（4）。
接下来对需要做的就是对$p(x_{1:n}|y_{1:n})$的估计。（思考$p(x_{1:n}|y_{1:n})$和$p(x_n|y_{1:n})$的区别与联系）。在有限状态空间的尔可夫模型中，式（6）中的积分可以转化为求和的形式。在线性高斯环境下的估计也可以通过著名的Kalman解决。对于非线性非高斯模型，Particle filter 提供了一种灵活有效的估计方式——它提供后验分布$p(x_{1:n}|y_{1:n})$的近似分布并且促进对$p(y_{1:n})$的近似计算 （思考促进对$p(y_{1:n})$的近似计算在particle filter中的作用）。

$p(x_{1:n}|y_{1:n})$和$p(x_n|y_{1:n})$的区别与联系
$p(x_{1:n}|y_{1:n})$表示序列近似分布的联合分布，而$p(x_n|y_{1:n})$表示边沿分布。Particle filter 中主要对边沿分布$p(x_n|y_{1:n})$进行估计。这就意味着在时刻1，我们希望得到（估计）$p(x_1|y_1)$和$p(y_1)$，在时刻2我们希望得到（估计）$p(x_2|y_{1:2})$和$p(y_{1:2})$。

滤波和边沿似然

后验分布$p(x_{1:n},y_{1:n})$满足
$$p(x_{1:n},y_{1:n})=p(x_{1:n})p(y_{1:n}|x_{1:n})$$
结合式（3）和（4）
得到
$$p(x_{1:n},y_{1:n})=p(x_{1:n-1},y_{1:n-1})f(x_n|x_{n-1})g(y_n|x_n)\text{\hspace{1em}(7)}$$
从而式（5）后验分布满足：
$$p(x_{1:n}|y_{1:n})=p(x_{1:n-1}|y_{1:n-1})\frac{f(x_n|x_{n-1})g(y_n|x_n)}{p(y_n|y_{1:n-1})}\text{\hspace{1em}(8)}$$
其中
$$p(y_n|y_{1:n-1})={\int }_{}p(x_{n-1}|y_{1:n-1})f(x_n|x_{n-1})g(y_n|x_n)d{x}_{n-1:n}\text{\hspace{1em}(9)}$$
式（9）的推导参考式（7）
对式（8）中关于$x_{1:n-1}$做积分，得到
$$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}(10)}$$
其中
$$p(x_n|y_{1:n-1})={\int }_{}f(x_n|x_{n-1})p(x_{n-1}|y_{1:n-1})d{x}_{n-1}\text{\hspace{1em}(11)}$$
式（10）就是常提到的系统的更新步骤，而式（11）是预测步骤。
需要提到的是，边沿分布$p(y_{1:n})$可以通过迭代得到：
$$p(y_{1:n})=p(y_1)\prod _{k=2}^{n}p(y_k|y_{1:k-1})\text{\hspace{1em}(12)}$$

$p(y_k|y_{1:k-1})$的形式如式（9）

现在总结一下我们做了些什么：通过马尔科夫过程的特点得到了后验联合分布关于状态方程和观测方程的表达形式式（7），然后通过式（7）得到了后验分布的条件概率形式。通过对前$n-1$个时刻的状态进行积分，得到了当前时刻的状态值在给定前$n$个观测值的条件下的概率分布$p(x_n|y_{1:n})$。在particle filter的当前教程中，particle filter的估计目标就是$p(x_n|y_{1:n})$——通过对系统的观测估计系统的状态真实值。

至此，第一部分总结完毕。本部分是particle filter 利用蒙特卡罗模型近似估计系统状态的基础，下一部分介绍蒙特卡罗如何作用于particle filter，遗留问题也会在第二部分更新。