【机器学习】粒子滤波原理推导（一）

前言

上一篇，我们已经对卡尔曼滤波进行数学推导和案例实现。可是，无论是隐马尔可夫模型还是卡尔曼滤波，它们或多或少的都存在着一些限制。对卡尔曼滤波而言，其要求噪声必须是高斯分布，并且必须是线性的。而对于粒子滤波而言，却没有这方面的限制。
数学基础：【概率论与数理统计知识复习-哔哩哔哩】

原理推导

对于粒子滤波，从贝叶斯的角度来看，它也是要求解Filter问题，即求解$P(z_t|x_1,x_2,\cdots ,x_t)$。z为隐变量，x为观测变量。我们对其进行一下定义
$$Z=\left(\begin{array}{c}{z}_1\\ ⋮\\ z_T\end{array}\right);X=\left(\begin{array}{c}{x}_1\\ ⋮\\ x_T\end{array}\right)$$
Z是隐序列，X为观测序列。

对于粒子滤波，依然满足两个假设，即

①齐次马尔可夫假设：当前隐状态仅和上一个时刻的隐状态有关。数学表达为
$$P(z_t|z_1,z_2,\cdots ,z_{t-1},x_1,\cdots .x_{t-1})=P(z_t|z_{t-1})$$
②观测独立性假设：当前观测状态仅和当前时刻的隐状态有关。数学表达为
$$P(x_t|x_1,x_2,\cdots ,x_{t-1},z_1,\cdots ,z_t)=P(x_t|z_t)$$
对于粒子滤波，其状态转移方程
$$x_t=f(x_{t-1},u)$$
$f$是一个函数的表达，其没有任何限制，可以是线性的，也可以不是线性的。$x_{t-1}$则是上一个时刻观测变量，而$u$是一个过程噪声，它也没有限制，可以是高斯的，也可以是非高斯的。

而观测方程
$$z_t=g(x_t,v)$$
同理对于$g$来说依然是一个没有限制的函数，$v$依然是一个噪声，可以服从任意概率分布。

流程图

问题

回到我们之前的问题，我们要求的是$P(z_t|x_1,x_2,\cdots ,x_t)$，这一次，不同于卡尔曼滤波可以对其进行求解，这一次因为它是非线性、非高斯的，我们压根就不能用上次的方法去求解。并且，如果你看过卡尔曼滤波中的对其的转化，你会发现要计算它是要进行积分运算的。积分一直都是一个难题，其复杂度太高。

那么，我们该如何去求解呢?实际上，对于Filter的应用问题，很多情况下我们并不是要求出的它的概率，而是更趋向于知道它的期望和方差是多少。比如上一个案例卡尔曼滤波的应用当中我们就仅仅是用到了它的期望卡尔曼滤波对运动轨迹优化。

所以，我们的问题就是要求解它的期望，得到期望之后自然能得方差。现在，我们定义它的期望为$\mu$

那么，如何去求这个期望呢？

蒙特卡罗

前面我们说过，积分问题是很难求解的，所以我们退而求其次，不去求它的概率分布，转而去求它的期望。可既然概率分布都求不出，又怎么才能去求出它的期望呢？

我们的方法就是采用蒙特卡罗的思想。即假设随机变量x的期望为
$$E(x)={\int }_{x}xP(x)dx=\frac{1}{n}\sum _{i=1}^n{x}^i$$
什么意思呢？我们知道，对于期望，我们是可以用均值去近似的，也就是对于随机变量$x$的期望，我们是可以从这个概率分布中抽样，不断地抽样出$x^1,\cdot ,x^n$，然后求出这些样本的均值。根据大数定律，只要抽样次数足够多，那么理论上我们的的均值就可以等于期望的值。

所以，要求出$(z_t|x_1,x_2,\cdots ,x_t)$这个随机变量的期望，我们就要从它的分布中去采样。得到样本。

但是很多情况下，由于这个分布的复杂性，我们往往无法从中去抽样。因此，我们会引入一个相对比较简单的分布，从这个相对简单的分布中去抽样，推导为
$$\begin{array}{c}\begin{array}{rl}E(x)=& {\int }_{x}xP(x)dx\\ =& {\int }_{x}x\frac{P(x)}{q(x)}q(x)dx\\ =& \frac{1}{n}\sum _{i=1}^n{x}^i\frac{P(x^i)}{q(x^i)}\end{array}\end{array}$$
其中，$\frac{P(x^i)}{q(x^i)}$又被称为权重，一般用$w^i$表示，即$w^i=\frac{P(x^i)}{q(x^i)}$

所以，对于随机变量$(z_t|x_1,x_2,\cdots ,x_t)$，我们也同理是要引入一个相对比较简单的概率分布，从中抽样，然后就可以求出它的期望。

SIS

前面，我们知道，我们要引入一个合适的分布，然后要从中采样，然后得到期望。但是，粒子滤波和卡尔曼滤波一样，在应用的时候，我们往往也都是需要求出$P(z_t|x_1,x_1,\cdots ,x_t)$的每一个时刻的值，也就是说，我们要求$P(z_1|x_0),P(z_2|x_1,x_2),\cdots ,P(z_t|x_1,x_1,\cdots ,x_t)$。我们在每一个时刻t，都必须要去采样，都要去计算每一个时刻的期望。

如果我们每一个时刻采样n次，那么我们就要采样$n\times T$次，如果我们每一次都要去计算这个权重值$w$，因为$P(z_t|x_1,x_2,\cdots ,x_t)$的复杂性，这无疑会大大增加计算量。即
$$\begin{array}{c}\begin{array}{rlr}& t=1:& \\ & & w_1^1,w_1^2,w_1^3,\cdots ,w_1^n\\ & t=2:& \\ & & w_2^1,w_2^2,w_2^3,\cdots ,w_2^n\\ & t=3:& \\ & & w_3^1,w_3^2,w_3^3,\cdots ,w_3^n\\ & ⋮& \\ & t=T:& \\ & & w_T^1,w_T^2,w_T^3,\cdots ,w_T^n\end{array}\end{array}$$
那我们就有一种想法，我们在$t$时刻计算出来的$w_t$权值，能不能在$t+1$时刻的时候复用，或者说，我们能不能找一出一个递推关系式，找出$w_t$和$w_{t+1}$的关系式，减少计算量？

事实上，它们还真的存在一定关系，推在递推关系

为了计算的简便，我们用不再去计算$P(z_t|x_1,x_2,\cdots ,x_t)$，转而去计算$P(z_1,z_2,\cdots ,z_t|x_1,x_2,\cdots ,x_t)$。

为了表达的简便，我们令$z_1,z_2,\cdots ,z_t=z_{1:t}$，同理得x的。所以我们就是要求出$P(z_{1:t}|x_{1:t})$。

可这样一来，我们求出的就是随机变量$z_{1:t}|x_{1:t}$所服从概率分布的期望，就不再就是原本的那个了。事实上，我们算出来的这个可以算是联合概率分布，而之前$P(z_t|x_1,x_2,\cdots ,x_t)$对应的概率分布则可以说是边缘分布。因此，我们只要在联合概率分布的基础上将其他时刻的$z$积分掉即可。

但是话说回来，积分会增加计算的复杂度。所以在实际的应用中，我们更趋向于计算出来的$P(z_{1:t}|x_{1:t})$的期望就近似等于我们Filter问题$P(z_t|x_1,x_2,\cdots ,x_t)$的期望。严格来说它们并不相等，只是为了计算的简便，我们一般认为它们近似相等

那么我们现在就去求出它们的权值的递推式吧

我们要求出$P(z_{1:t}|x_{1:t})$，要引进一个较为简单的概率分布$q(z_{1:t}|x_{1:t})$

我们先来看$P(z_{1:t}|x_{1:t})$
$$\begin{array}{c}\begin{array}{rl}P(z_{1:t}|x_{1:t})=& \frac{P(z_{1:t},x_{1:t})}{P(x_{1:t})}\\ =& \frac{P(x_t|z_{1:t},x_{1:t-1})P(z_{1:t},x_{1:t-1})}{P(x_{1:t})}\\ =& \frac{P(x_t|z_t)P(z_{1:t},x_{1:t-1})}{P(x_{1:t})}\\ =& \frac{P(x_t|z_t)P(z_t|z_{1:t-1},x_{1:t-1})P(z_{1:t-1},x_{1:t-1})}{P(x_{1:t})}\\ =& \frac{P(x_t|z_t)P(z_t|z_{t-1})P(z_{1:t-1},x_{1:t-1})}{P(x_{1:t})}\\ =& \frac{P(x_t|z_t)P(z_t|z_{t-1})P(z_{1:t-1}|x_{1:t-1})P(x_{1:t-1})}{P(x_{1:t})}\end{array}\end{array}$$

其中，里面用到了前面给到的两个假设。$P(x_{1:t})$和$P(x_{1:t-1})$是可以计算出来的，因为它们是给定的。因此
$$P(z_{1:t}|x_{1:t})\propto P(x_t|z_t)P(z_t|z_{t-1})P(z_{1:t-1}|x_{1:t-1})$$
再来看我们引入的一个比较简单的概率分布$q(z_{1:t}|x_{1:t})$，它是我们任意给定的，所以可以随意设定它。我们令
$$q(z_{1:t}|x_{1:t})=q(z_t|z_{1:t-1},x_{1:t})q(z_{1:t-1}|x_{1:t-1})\text{\hspace{1em}(1)}$$
请注意这个不是由贝叶斯定理得出的，是我们自己设定的。因为$q(z_{1:t}|x_{1:t})$本来就是引进来的，我们想令它等于何种形式就等于何种形式。

所以对于t时刻的权重，我们有
$$\begin{array}{c}\begin{array}{rl}{w}_t=& \frac{P(z_{1:t}|x_{1:t})}{q(z_{1:t}|x_{1:t})}\\ \propto & \frac{P(x_t|z_t)P(z_t|z_{t-1})P(z_{1:t-1}|x_{1:t-1})}{q(z_t|z_{1:t-1},x_{1:t})q(z_{1:t-1}|x_{1:t-1})}\\ =& \frac{P(x_t|z_t)P(z_t|z_{t-1})}{q(z_t|z_{1:t-1},x_{1:t})}{w}_{t-1}\end{array}\end{array}$$
好，现在，我们得到了权重$w$的递推式。但是现在还有一个问题，即便是每一次我们计算权值的时候不再需要计算复杂的$P(z_{1:t}|x_{1:t})$。但是我们抽样的时候也需要从里面去抽样，这是一个t维的分布，复杂度不用说，定然麻烦。所以我们能不能也像权重$w$一样，找一个递推式？

简单起见，我们就以$P(z_{1:2}|x_{1:2})$为例假设我们已经求出了$P(z_{1:t-1}|x_{1:t-1})$，我们用$P(z_{1:t-1}^{\ast }|x_{1:t-1})$作标记，要求出$P(z_{1:t}|x_{1:t})$

$$\begin{array}{c}\begin{array}{rl}P(z_{1:t}|x_{1:t})=& \frac{P(z_{1:t},x_{1:t})}{P(x_{1:t})}\\ \propto & P(z_{1:t},x_{1:t})\\ =& P(x_t|z_{1:t},x_{1:t-1})P(z_t|z_{1:t-1},x_{1:t-1})P(z_{1:t-1}|x_{1:t-1})P(x_{1:t-1})\\ =& P(x_t|x_{t-1})P(z_t|z_{t-1})P(z_{1:t-1}|x_{1:t-1})P(x_{1:t-1})\\ \propto & P(z_t|z_{t-1})P(z_{1:t-1}|x_{1:t-1})\end{array}\end{array}$$
我们知道，我们已经求出了前一个时刻的$P(z_{1:t-1}^{\ast }|x_{1:t-1})$，假设我们抽样是从$P(z_t|z_{t-1})$这个概率分布中去抽样的，得到的结果记作$P(z_t^{\ast }|z_{t-1}^{\ast })$，那么我们往回推，就能得到当前时刻的$P(z_{1:t}|x_{1:t})$。这很容易看出来是完全可行的。也就是说，我们每一次只需要从$P(z_t|z_{t-1})$抽样出来即可。

同理我们引入的概率分布中(式1)，那么同样我们也只需要从概率$q(z_t|z_{1:t-1},x_{1:t})$中去采样即可。因为按照我们上面的说法，t时刻的概率和t-1时刻的概率存在关系即可。

概率分布选定

所以，如何去选定$q(z_t|z_{1:t-1},x_{1:t})$的概率分布呢？一般情况下，为了简单起见，我们认为其同样符合齐次马尔可夫性质，所以$q(z_t|z_{1:t-1},x_{1:t})=q(z_t|z_{t-1},x_{1:t})$。在很多时候，选定$q(z_t|z_{1:t-1},x_{1:t})$这个概率分布的时候，我们认为最合适的概率分布是$P(z_t|z_{t-1})$（但请注意这个并不是最优的，只是为了简化运算）。

那么，我们最终的权重就可以简化为
$$\begin{array}{c}\begin{array}{rl}{w}_t\propto & \frac{P(x_t|z_t)P(z_t|z_{t-1})}{q(z_t|z_{1:t-1},x_{1:t})}{w}_{t-1}\\ =& \frac{P(x_t|z_t)P(z_t|z_{t-1})}{P(z_t|z_{t-1})}{w}_{t-1}\\ =& P(x_t|z_t)w_{t-1}\end{array}\end{array}$$
并且，一般情况下，我们都会对权重进行归一化。

另外，采样出来的样本，在粒子滤波中，实际上又被称为粒子。

存在问题

以上这种方法，看起来很好了，但并非是一点问题都没有，它会产生权值退化的问题。下一篇我们就来讲讲如何去解决。