粒子滤波学习笔记（二）

接上一篇，现在我们介绍蒙特卡罗积分和几种采样方法。目前我只是对粒子滤波有一个比较基础的认识，基本原理上只涉及到了蒙特卡罗采样和Important Sampling,Sequential Important Sampling,Sampling-Importance Resampling三种，所以就只介绍这几种。

1 蒙特卡罗积分

蒙特卡罗方法又被称为随机采样法或统计试验法，其基本思路就是将实际求解的问题描述成某种随机变量，然后从已知的概率分布采样，建立各种估计量，得到所求的解。所谓的蒙特卡罗积分也是这个思想。

假设$D$是$n$维空间$\mathcal R^n$的一个区域，考虑多重积分

$$I=\int_D G(\mathbf x)d\mathbf x$$ 式中， $D$为积分区域。

蒙特卡罗积分就是将积分值$I$看作是某种随机变量的数学期望，于是有下式

$$I=\int_D f(\mathbf x)p(\mathbf x)d\mathbf x$$ 假设 $G(\mathbf x)=f(\mathbf x)p(\mathbf x)$, $p(\mathbf x)$可以看作是 $f(\mathbf x)$的数学期望，于是自然就有要求 $p(\mathbf x)\ge0,\int_{R_n}p(\mathbf x)d\mathbf x=1$。

假设从概率密度$p(\mathbf x)$中采样得到$N_p$个独立同分布的样本值$[\mathbf x^{1},...,\mathbf x^{N_p}]$,那么$G(\mathbf x)$的蒙特卡罗估计可以为

$$\hat{f_{N_p}}=\frac{1}{N_p}\sum_{i=0}^{N_p}f(\mathbf x^{i})$$

由大数定理可以知道，当$[\mathbf x^{1},...,\mathbf x^{N_p}]$相互独立而且$N_p\to\infty$的时候，我们上面得到的估计$\hat{f_{N_p}}$就会以概率１收敛与积分的真值，也就是

$$p(\lim_{N_p\to \infty}\hat{f_{N_p}}=I)=1$$
由中心极限定理可以知道，收敛的速率为 $$\sqrt{N_p}(\hat{f_{N_p}}-\mathbb E(f))\sim \mathscr N(0, \sigma^2)$$ $\sigma^2$是 $f(\mathbf x)$的方差。
从以上推理我们基本可以知道，估计的误差率是 $\mathscr O({N_p}^{-1/2})$，而且估计的方差和采样的个数成反比。

2 几种采样方法

2.1 重要采样

重要采样是采样分布的“重要”区域，这样的话可以提高计算效率。一般来说先验分布$p(\mathbf x)$难以采样，我们间接的选取一个重要分布$q(\mathbf x)$来取代原来真实的分布,于是有下面的公式：

$$\int f(\mathbf x)p(\mathbf x)d\mathbf x=\int f(\mathbf x)\frac{p(\mathbf x)}{q(\mathbf x)}q(\mathbf x)d\mathbf x$$
蒙特卡罗重要性采样就是从重要分布 $q(\mathbf x)$采大量的(比如说 $\mathbf N_p$)独立样本以此来获得加权的求和项来估计 $f(\mathbf x)$。
$$\hat f=\frac 1{N_p}\sum_{i=1}^{N_p}W(\mathbf x^{i})f(\mathbf x^{i})$$
其中 $W(\mathbf x^{i})=\frac{p(\mathbf x^{i})}{q(\mathbf x^{i})}$叫做重要性权值，但是如果归一化因子 $p(\mathbf x)$不知道，我们只能说权值正比于一个常数，也就是 $W(\mathbf x^{i})\propto\frac{p(\mathbf x^{i})}{q(\mathbf x^{i})}$，为了保证 $\sum_{i=1}^{N_p}W(\mathbf x^{i})=1$,我们通常对权值做归一化，也就是 $$\hat f=\frac{\frac 1{N_p}\sum_{i=1}^{N_p}W(\mathbf x^{i})f(\mathbf x^{i})}{\frac 1{N_p}\sum_{i=1}^{N_p}W(\mathbf x^{i})}$$ 其中 $\hat W(\mathbf x^{(i)})=\frac {W(\mathbf x^{i})}{\sum_{j=1}^{N_p}W(\mathbf x^{j})}$叫做重要权重。

2.2 序贯重要采样

一个好的重要分布对于重要采样的效果是很重要的，因此怎么选择重要分布是成功应用重要采样的关键，但是一般来说寻找一个重要分布尤其是在高维空间中寻找是很难的。一个很自然的想法就是我们可以按照顺序来建立重要分布，这就是所谓的序贯重要采样（SIS）的基本思想。

如果我们的重要分布式按照一个因式分解的形式选择的，如下所示，那么重要采样就可以采用迭代的方式来进行。

$$q(\mathbf x_{0:n}|\mathbf y_{0:n})=q(\mathbf x_0)\prod_{t=1}^n{q(\mathbf x_t|\mathbf x_{0:t-1},\mathbf y_{0:t})}$$

为了便于理解，我们考虑一个简化的例子，也就是不考虑条件，按照条件概率我们可以有

$$p(\mathbf x_{0:n})=p(\mathbf x_0)p(\mathbf x_1|\mathbf x_0)...p(\mathbf x_n|\mathbf x_0,...,\mathbf x_{n-1})$$
$$q(\mathbf x_{0:n})=q_0(\mathbf x_0)q_1(\mathbf x_1|\mathbf x_0)...q_n(\mathbf x_n|\mathbf x_0,...,\mathbf x_{n-1})$$
因此重要权重就可以看作 $$W(\mathbf x_{0:n})=\frac{p(\mathbf x_0)p(\mathbf x_1|\mathbf x_0)...p(\mathbf x_n|\mathbf x_0,...,\mathbf x_{n-1})}{q_0(\mathbf x_0)q_1(\mathbf x_1|\mathbf x_0)...q_n(\mathbf x_n|\mathbf x_0,...,\mathbf x_{n-1})}$$
于是也就可以用迭代的形式来计算 $$W(\mathbf x_{0:n})=W_{n-1}(\mathbf x_{0:n-1}) \frac{p(\mathbf x_n|\mathbf x_{0:n-1})}{q_n(\mathbf x_n|\mathbf x_{0:n-1})}$$

2.3 重采样

重采样的目的可以基于以下几点：
(1)如果重要权重不是均匀分布的，那么计算一些“琐碎”的权重很没有必要，是对计算能力的浪费；
(2)当重要权值不准确的时候，重采样可以提供选择重要样本的机会。
注意：
(1)一般来说，重采样是在两次重要采样之间执行的。在重采样的时候，粒子和相关重要权重$\{\mathbf x^{(i)},\hat W^{(i)}\}$被新的等权重的粒子替换。重采样可以每一步都执行也可以在需要的时候执行。
(2)重采样不会提高当前状态估计的精度，因为重采样也会带来额外的蒙特卡罗方差。

在目前我只用到了一种重要采样算法，我在这里给出：
(1)产生一个均匀分布$u\sim \mathcal U(0,1)$,建立一个权重的累积概率，也就是$s_i=\sum_{j=1}^i\hat W^{(j)}$;
(2)寻找一个$s_i$使得$s_{i-1}\le u \le s_i$,索引为$i$的粒子被选择；
(3)给定$\{\mathbf x^{(j)},\hat W^{(j)}\}\ \ j=1,...,N_p$,按照相关的权重$\hat W^{(j)}$复制粒子$\mathbf x^{(j)}$,进而得到新的粒子$\mathbf x^{(i)}$;
(4)$W^{(i)}=1/N_p$.