Particle Filter Tutorial 粒子滤波：从推导到应用（二）

二、蒙特卡洛采样

假设我们能从一个目标概率分布p(x)中采样到一系列的样本（粒子） $x_{1},...,x_{N}$ ，（至于怎么生成服从p(x)分布的样本，这个问题先放一放），那么就能利用这些样本去估计这个分布的某些函数的期望值。譬如：

$E\left ( f\left ( x \right ) \right )=\int_{a}^{b}f\left ( x \right )p\left ( x \right )dx$

$Var\left ( f\left ( x \right ) \right )=E\left ( f\left ( x \right ) - E\left (f \left ( x \right ) \right )\right )^{2}=\int_{a}^{b}\left ( f\left ( x \right ) - E\left (f \left ( x \right ) \right )\right )^{2}p\left ( x \right )dx$

上面的式子其实都是计算期望的问题，只是被积分的函数不同。

蒙特卡洛采样的思想就是用平均值来代替积分，求期望：

$E\left ( f\left ( x \right ) \right )\approx \frac{f\left ( x_{1} \right )+...+f\left ( x_{N}\right )}{N}$

这可以从大数定理的角度去理解它。我们用这种思想去指定不同的f(x)以便达到估计不同东西的目的。比如：要估计一批同龄人的体重，不分男女，在大样本中男的有100个，女的有20个，为了少做事，我们按比例抽取10个男的，2个女的，测算这12个人的体重求平均就完事了。注意这里的按比例抽取，就可以看成从概率分布p(x)中进行抽样。

下面再看一个稍微学术一点的例子：

假设有一粒质地均匀的骰子。规定在一次游戏中，连续四次抛掷骰子，至少出现一次6个点朝上就算赢。现在来估计赢的概率。我们用 $x_{k}^{\left ( n \right )}$ 来表示在第n次游戏中，第k次投掷的结果，k=1...4。对于分布均匀的骰子，每次投掷服从均匀分布,即：

$x_{k}\sim u\left ( 1,6 \right )$

这里的区间是取整数，1,2,3,4,5,6，代表6个面。由于每次投掷都是独立同分布的，所以这里的目标分布p(x)也是一个均匀分布 $\chi =\left \{ 1,...6\right \}^{4}$ 。一次游戏就是 $\chi$ 空间中的一个随机点。

为了估计取胜的概率，在第n次游戏中定义一个指示函数:

$f\left ( x^{\left ( n \right )} \right ) = \mathbb{I}\left \{ 0<\sum_{k=1}^{4}\mathbb{I}\left \{ x_{k}^{(n)}=6 \right \} \right \}$

其中，指示函数I{x }是指，若x的条件满足，则结果为1，不满足结果为0。回到这个问题，这里函数 f()的意义就是单次游戏中，若四次投掷中只要有一个6朝上，f()的结果就会是1。由此，就可以估计在这样的游戏中取胜的期望，也就是取胜的概率：

$\theta =E\left ( f(x) \right )\approx \frac{1}{N}\sum_{n=1}^{N}f(x^{(n)})$

当抽样次数N足够大的时候，上式就逼近真实取胜概率了，看上面这种估计概率的方法，是通过蒙特卡洛方法的角度去求期望达到估计概率的目的。是不是就跟我们抛硬币的例子一样，抛的次数足够多就可以用来估计正面朝上或反面朝上的概率了。

当然可能有人会问，这样估计的误差有多大，对于这个问题，有兴趣的请去查看我最下面列出的参考文献2。（啰嗦一句：管的太多太宽，很容易让我们忽略主要问题。博主就是在看文献过程中，这个是啥那个是啥，都去查资料，到头来粒子滤波是干嘛完全不知道了，又重新看资料。个人感觉有问题还是先放一放，主要思路理顺了再关注细节。）

接下来，回到我们的主线上，在滤波中蒙特卡洛又是怎么用的呢？

由上面我们知道，它可以用来估计概率，而在上一节中，贝叶斯后验概率的计算里要用到积分，为了解决这个积分难的问题，可以用蒙特卡洛采样来代替计算后验概率。

假设可以从后验概率中采样到N个样本，那么后验概率的计算可表示为：

$\hat{p}(x_{n}|y_{1:k})=\frac{1}{N}\sum_{i=1}^{N}\delta (x_{n}-x_{n}^{(i)})\approx p(x_{n}|y_{1:k})$

其中，在这个蒙特卡洛方法中，我们定义 $f(x)=\delta (x_{n}-x_{n}^{(i)})$ ,是狄拉克函数(dirac delta function)，跟上面的指示函数意思差不多。

看到这里，既然用蒙特卡洛方法能够用来直接估计后验概率，现在估计出了后验概率，那到底怎么用来做图像跟踪或者滤波呢？要做图像跟踪或者滤波，其实就是想知道当前状态的期望值：

$E[f(x_{n})]\approx \int f(x_{n})\hat{p}(x_{n}|y_{1:k})dx_{n}$

$=\frac{1}{N}\sum_{i=1}^{N}\int f(x_{n})\delta (x_{n}-x_{n}^{(i)})dx_{n}$

$=\frac{1}{N}\sum_{i=1}^{N}f(x_{n}^{(i)})$ (1)

也就是用这些采样的粒子的状态值直接平均就得到了期望值，也就是滤波后的值，这里的 f(x) 就是每个粒子的状态函数。这就是粒子滤波了，只要从后验概率中采样很多粒子，用它们的状态求平均就得到了滤波结果。

思路看似简单，但是要命的是，后验概率不知道啊，怎么从后验概率分布中采样！所以这样直接去应用是行不通的，这时候得引入重要性采样这个方法来解决这个问题。

三、重要性采样

无法从目标分布中采样，就从一个已知的可以采样的分布里去采样如 q(x|y)，这样上面的求期望问题就变成了：

$E\left [ f\left ( x_{k} \right ) \right ] = \int f(x_{k})\frac{p(x_{k}|y_{1:k})}{q(x_{k}|y_{1:k})}q(x_{k}|y_{1:k})dx_{k}$

$= \int f(x_{k})\frac{p(y_{1:k}|x_{k})p(x_{k})}{p(y_{1:k})q(x_{k}|y_{1:k})}q(x_{k}|y_{1:k})dx_{k}$

$= \int f(x_{k})\frac{W_k(x_{k})}{p(y_{1:k})}q(x_{k}|y_{1:k})dx_{k}$ (2)式

其中

$W_k(x_k) = \frac{p(y_{1:k}|x_k)p(x_k)}{q(x_k|y_{1:k})}$ $\propto \frac{p(x_k|y_{1:k})}{q(x_k|y_{1:k})}$

由于：

$p(y_{1:k})=\int p(y_{1:k}|x_{k})p(x_{k})dx_k$

所以(2)式可以进一步写成：

$E[f(x_k)]=\frac{1}{p(y_{1:k})}\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k$

$=\frac{\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k}{\int p(y_{1:k}|x_{k})p(x_k)dx_k}$

$=\frac{\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k}{\int W_k(x_k)q(x_k|y_{1:k})dx_k}$

$=\frac{E_{q(x_k|y_{1:k})} [W_k(x_k)f(x_k)]}{E_{q(x_k|y_{1:k})} [W_k(x_k)]}$ (3)式

上面的期望计算都可以通过蒙特卡洛方法来解决它，也就是说，通过采样N个样本 $\left \{ x_k^{(i)} \right \}\sim q(x_k|y_{1:k})$ ,用样本的平均来求它们的期望，所以上面的（3）式可以近似为：

$E\left [ f(x_k) \right ]\approx \frac{\frac{1}{N}\sum_{i=1}^{N}W_k(x_k^{(i)})f(x_k^{(i)})}{\frac{1}{N}\sum_{i=1}^{N}W_k(x_k^{(i)})}$

(4)式

其中：

$\tilde{W}_k(x_k^{(i)})=\frac{W_k(x_k^{(i)})}{\sum_{i=1}^{N}W_k(x_k^{(i)})}$

这就是归一化以后的权重，而之前在(2)式中的那个权重是没有归一化的。

注意上面的(4)式，它不再是（1）式中所有的粒子状态直接相加求平均了，而是一种加权和的形式。不同的粒子都有它们相应的权重，如果粒子权重大，说明信任该粒子比较多。

到这里已经解决了不能从后验概率直接采样的问题，但是上面这种每个粒子的权重都直接计算的方法，效率低，因为每增加一个采样，p( x(k) |y(1:k))都得重新计算，并且还不好计算这个式子。所以求权重时能否避开计算p(x(k)|y(1:k))？而最佳的形式是能够以递推的方式去计算权重，这就是所谓的序贯重要性采样（SIS），粒子滤波的原型。

下面开始权重w递推形式的推导：

假设重要性概率密度函数，这里x的下标是0:k，也就是说粒子滤波是估计过去所有时刻的状态的后验。假设它可以分解为：

后验概率密度函数的递归形式可以表示为：

$\propto p(y_k|x_k)p(x_k|x_{k-1})p(x_{0:k-1}|Y_{k-1})$

其中，为了表示方便，将 y(1:k) 用 Y(k) 来表示，注意 Y 与 y 的区别。同时，上面这个式子和上一节贝叶斯滤波中后验概率的推导是一样的，只是之前的x(k)变成了这里的x(0:k)，就是这个不同，导致贝叶斯估计里需要积分，而这里后验概率的分解形式却不用积分。

粒子权值的递归形式可以表示为:

( 5)式

注意，这种权重递推形式的推导是在前面（2）式的形式下进行推导的，也就是没有归一化。而在进行状态估计的公式为这个公式中的的权重是归一化以后的，所以在实际应用中，递推计算出w(k)后,要进行归一化，才能够代入(4)式中去计算期望。同时，上面(5)式中的分子是不是很熟悉，在上一节贝叶斯滤波中我们都已经做了介绍，p( y|x ),p( x(k)|x(k-1) )的形状实际上和状态方程中噪声的概率分布形状是一样的，只是均值不同了。因此这个递推的(5)式和前面的非递推形式相比，公式里的概率都是已知的，权重的计算可以说没有编程方面的难度了。权重也有了以后，只要进行稍微的总结就可以得到SIS Filter。

四、Sequential Importance Sampling (SIS) Filter

在实际应用中我们可以假设重要性分布q()满足：

$q(x_{k}|x_{0:k-1},y_{1:k})=q(x_{k}|x_{k-1},y_{k})$

这个假设说明重要性分布只和前一时刻的状态x(k-1)以及测量y(k)有关了，那么(5)式就可以转化为：

$w_k^{(i)}\propto w_{k-1}^{(i)}\frac{p(y_k|x_k^{(i)})p(x_k^{(i)}|x_{k-1}^{(i)})}{q(x_{k}^{(i)}|x_{k-1}^{(i)},y_{k})}$

在做了这么多假设和为了解决一个个问题以后，终于有了一个像样的粒子滤波算法了，他就是序贯重要性采样滤波。

下面用伪代码的形式给出这个算法：

----------------------pseudo code-----------------------------------

For i=1:N

(1)采样： $x_k^{(i)}\sim q(x_{k}^{(i)}|x_{k-1}^{(i)},y_{k})$ ；

(2)根据 $w_k^{(i)}\propto w_{k-1}^{(i)}\frac{p(y_k|x_k^{(i)})p(x_k^{(i)}|x_{k-1}^{(i)})}{q(x_{k}^{(i)}|x_{k-1}^{(i)},y_{k})}$ 递推计算各个粒子的权重；

End For

粒子权值归一化。粒子有了，粒子的权重有了，就可以由(4)式,对每个粒子的状态进行加权去估计目标的状态了。

-----------------------end -----------------------------------------------

这个算法就是粒子滤波的前身了。只是在实际应用中，又发现了很多问题，如粒子权重退化的问题，因此就有了重采样( resample )，就有了基本的粒子滤波算法。还有就是重要性概率密度q()的选择问题，等等。都留到下一章去解决。

在这一章中，我们是用的重要性采样这种方法去解决的后验概率无法采样的问题。实际上，关于如何从后验概率采样，也就是如何生成特定概率密度的样本，有很多经典的方法（如拒绝采样，Markov Chain Monte Carlo，Metropolis-Hastings 算法，Gibbs采样），这里面可以单独作为一个课题去学习了，有兴趣的可以去看看《统计之都的一篇博文》，强烈推荐，参考文献里的前几个也都不错。

（转载请注明作者和出处：http://blog.csdn.net/heyijia0327 未经允许请勿用于商业用途）

reference：

1. Gabriel A. Terejanu 《Tutorial on Monte Carlo Techniques》

2. Taylan Cemgil 《A Tutorial Introduction to Monte Carlo methods, Markov Chain Monte Carlo and Particle Filtering》

3. M. Sanjeev Arulampalam 《A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking》

4. ZHE CHEN 《Bayesian Filtering: From Kalman Filters to Particle Filters, and Beyond》

5.百度文库《粒子滤波理论》

6. Haykin 《Neural Networks and learning Machines 》Chapter 14

7. 统计之都 <LDA-math-MCMC 和 Gibbs Sampling>

Particle Filter Tutorial 粒子滤波：从推导到应用（二）相关推荐

Particle Filter Tutorial 粒子滤波：从推导到应用（四）
六.Sampling Importance Resampling Filter (SIR) SIR滤波器很容易由前面的基本粒子滤波推导出来,只要对粒子的重要性概率密度函数做出特定的选择即可.在SIR中 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（一）
前言: 早几年前写博客的时候,编辑公式都是在线网站编辑的,有时候这个网站不稳定,贴的公式看不到.如果有需要,pdf版下载地址点击打开链接,给您阅读带来的不便表示抱歉,祝大家好运. 博主在自主学习粒子滤 ...
Particle Filter Tutorial 粒子滤波：从推导到应用（三）
五.重采样在应用SIS 滤波的过程中,存在一个退化的问题.就是经过几次迭代以后,很多粒子的权重都变得很小,可以忽略了,只有少数粒子的权重比较大.并且粒子权值的方差随着时间增大,状态空间中的有效粒子数 ...
Particle Filter Tutorial 粒子滤波：从推导到应用
转载自作者白巧克力亦唯心的文章:粒子滤波文章目录 1.贝叶斯滤波 2.蒙特卡洛采样 3.重要性采样 4.SIS Filter(Sequential Importance Sampling) 5.重采 ...
粒子滤波的推导到应用
本文转自 http://blog.csdn.net/heyijia0327, Particle Filter Tutorial 粒子滤波:从推导到应用系列一.贝叶斯滤波假设有一个系统,我们知道它的 ...
学习笔记：粒子滤波的推导过程
前面的文章我们了解了到对于一个线性动态系统(比如说卡尔曼滤波)中的filtering问题可以用预测和跟新两个步骤求解,因为它们的后验概率和条件概率都服从高斯分布,所以比较容易求出解析解,但是对于非线性 ...
视频跟踪算法之粒子滤波
1. 写在前面最近在看视频跟踪方面的一些硕博士毕业论文,几乎看到的每一篇都会涉及到粒子滤波算法,所以这段时间花了很多时间在看相关的内容. 浏览了大量的博客和文章,跟着不停推导公式,感觉还是无法完全掌 ...
双通道粒子滤波(PF)原理详解
1.粒子滤波是动态系统非线性非高斯的情况 2.粒子滤波无法像卡尔曼滤波求得解析解,需要用蒙特卡洛采样 3.采取重要性采样方法近似计算也无法直接求解粒子滤波问题 4.序列重要性采样得到t时刻与t-1时刻 ...
[转]粒子滤波（particle filtering）的思路发展过程及应用（详细深度好文）
粒子滤波作为视觉SLAM中后端进行状态估计的主要算法之一,很好的完成了扩展卡尔曼滤波无法有效处理的复杂状态方程下的状态估计任务.这篇文章详细地描述了粒子滤波的思想历程,即如何一步步从简单的状态估计.采 ...

Particle Filter Tutorial 粒子滤波：从推导到应用（二）

（转载请注明作者和出处：http://blog.csdn.net/heyijia0327 未经允许请勿用于商业用途）

Particle Filter Tutorial 粒子滤波：从推导到应用（二）相关推荐

最新文章

热门文章