机器学习 [白板推导]（十二）[卡曼滤波、粒子滤波]

15. 线性动态系统（卡曼滤波，Kalman Filter）

15.1. 概述

15.1.1. 背景介绍

  变量随时间变化的系统叫做动态系统，其中隐变量取值离散的是隐马尔可夫模型（HMM），而隐变量取值连续的分为线性动态系统和粒子滤波（Particular Filter），前者各变量和变量各时刻满足线性关系和高斯分布，后者则不满足。

  线性关系：时刻之间：$z_t=A\cdot z_{t-1}+B+\epsilon$；变量之间：$x_t=C\cdot z_t+D+\delta$；其中 $\epsilon ,\delta$ 是噪声，在线性动态系统中服从高斯分布。

15.1.2. 模型定义

  设噪声 $\epsilon \sim \mathcal{N}(0,Q)$，$\delta \sim \mathcal{N}(0,R)$，同时设初始状态 $z_1\sim \mathcal{N}({\mu }_1,{\Sigma }_1)$，则可以通过高斯分布依次计算 $z_t|z_{t-1}\sim \mathcal{N}(A\cdot z_{t-1}+B,Q)$，$x_t|z_t\sim \mathcal{N}(C\cdot z_t+D,R)$ 得到所有变量的分布，因此模型参数为 $\theta =(A,B,C,D,Q,R,{\mu }_1,{\Sigma }_1)$.

  线性动态系统和HMM一样，也服从一阶齐次马尔可夫假设和观测独立假设（见上篇笔记机器学习 [白板推导]（十一））。

15.2. Filtering

  首先，filtering公式为 $P(z_t|x_{1:t})$，是一种Online算法，通过 t 时刻及之前的观测变量来推测当前状态，这既是模型求解（Learning）的必要信息，同时也是对prediction结果 $P(z_t|x_{1:t-1})$ 的修正（Correction）。

  filtering也可以递推求解，推导为：
$$\begin{array}{rl}P(z_t|x_{1:t})& =\frac{P(z_t,x_{1:t})}{P(x_{1:t})}=\frac{\stackrel{P(x_t|z_t)}{\overbrace{P(x_t|z_t,x_{1:t-1})}}\cdot P(z_t,x_{1:t-1})}{P(x_{1:t})}\\ & =\frac{P(x_t|z_t)\cdot P(z_t|x_{1:t-1})\cdot P(x_{1:t-1})}{P(x_{1:t})}\\ & =\underset{\text{Emission Probability}}{\underbrace{P(x_t|z_t)}}\cdot \underset{\text{Prediction of t-1}}{\underbrace{P(z_t|x_{1:t-1})}}\cdot \underset{\text{Constant}}{\underbrace{\frac{P(x_{1:t-1})}{P(x_{1:t})}}},\end{array}\text{\hspace{1em}(15.1)}$$
其中Prediction可以求解为：
$$\begin{array}{rl}P(z_t|x_{1:t-1})& ={\int }_{z_{t-1}}P(z_t,z_{t-1}|x_{1:t-1})d{z}_{t-1}\\ & ={\int }_{z_{t-1}}\underset{\underset{\text{State Transition Probability}}{\underbrace{P(z_t|z_{t-1})}}}{\underbrace{P(z_t|z_{t-1},x_{1:t-1})}}\cdot \underset{\text{Filtering of t-1}}{\underbrace{P(z_{t-1}|x_{1:t-1})}}d{z}_{t-1},\end{array}\text{\hspace{1em}(15.2)}$$
因此可以通过 $t-1$ 时刻的filtering计算 $t$ 时刻的prediction，进而计算t时刻的filtering，循环迭代进行求解。

  总结来说，Filtering的过程可以写作以下步骤：

• $t=1$，$P(z_1|x_1)=\frac{P(z_1,x_1)}{P(x_1)}=\frac{P(x_1|z_1)\cdot P(z_1)}{{\int }_{z_1}P(x_1|z_1)\cdot P(z_1)d{z}_1}$，并用上面的公式计算Prediction $P(z_2|x_1)$.
• $t=2$，根据上面的递推关系计算 $P(z_2|x_1,x_2)$（作为filtering的更新，update，以及prediction的修正，correction）和 $P(z_3|x_1,x_2)$（新的prediction）。
• ……
• $t=T$，计算 $P(z_T|x_{1:T-1})$.

  在整个计算过程中，Filtering和Prediction的计算结果均为高斯分布，其参数（均值和标准差）可以通过之前时刻的计算结果以及模型参数算出，推导过程参考机器学习 [白板推导]（一）中高斯分布的性质，这里不再详细介绍。

16. 粒子滤波（Particular Filtering）

16.1. 背景介绍

  在15. 线性动态系统（卡曼滤波）中，可以通过filtering和prediction交替迭代的方法求得隐状态序列，但应对非线性非高斯模型，无法求出解析解，只能通过采样模拟的方法求出近似解，如此便诞生了粒子滤波。

  粒子滤波应用范围非常之广，且已经产生了许多研究成果和变体，这里只介绍最基本的粒子滤波算法。

16.2. 算法介绍

16.2.1. 重要度采样及针对filtering的改进 —— 序列重要度采样（Sequential Importance Sampling，SIS）

  首先回顾重要度采样，给定一个复杂的分布 $p(x)$，需要计算期望 $E[f(x)]$ 时，直接从 $p(x)$ 中采样较为困难，此时可以设置一个简单的建议分布 $q(x)$，此时 $E[f(x)]={\int }_z\left[f(x)p(x)\right]dz={\int }_z\left[f(x)\cdot \frac{p(x)}{q(x)}\cdot q(x)\right]dz$，采样变得容易。

  将重要度采样应用于动态系统的filtering问题时，对每个时刻 $t$ 和每个样本 $x_t^{(i)}$，都需要计算重要度 $w_t^{(i)}=\frac{p(z_t^{(i)}|x_{1:t}^{(i)})}{q(z_t^{(i)}|x_{1:t}^{(i)})}$，而 $p(z_t^{(i)}|x_{1:t}^{(i)})$ 的计算是复杂的，还要计算N次（N个样本）。

  将filtering问题的计算目标进行调整，由 $p(z_t^{(i)}|x_{1:t}^{(i)})$ 改为 $p(z_{1:t}^{(i)}|x_{1:t}^{(i)})$（同样可以达到目标），此时重要度为 $w_t^{(i)}\propto \frac{p(z_{1:t}^{(i)}|x_{1:t}^{(i)})}{q(z_{1:t}^{(i)}|x_{1:t}^{(i)})}$，分别推导分子和分母的递推式，如果可以建立递推关系，则可以简化运算：

• 对于分子，可以推导如下
  $$\begin{array}{rl}p(z_{1:t}^{(i)}|x_{1:t}^{(i)})& =\frac{p(z_{1:t}^{(i)},x_{1:t}^{(i)})}{p(x_{1:t}^{(i)})}\\ & =\frac{\stackrel{\stackrel{\text{Emission Probability}}{\overbrace{p(x_t^{(i)}|z_t^{(i)})}}}{\overbrace{p(x_t^{(i)}|z_{1:t}^{(i)},x_{1:t-1}^{(i)})}}\cdot p(z_{1:t}^{(i)},x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ & =\frac{p(x_t^{(i)}|z_t^{(i)})\cdot \stackrel{\stackrel{\text{State Transition Probability}}{\overbrace{p(z_t^{(i)}|z_{t-1}^{(i)})}}}{\overbrace{p(z_t^{(i)}|z_{1:t-1}^{(i)},x_{1:t-1}^{(i)})}}\cdot p(z_{1:t-1}^{(i)},x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ & =\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})\cdot p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ & =\left[\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\cdot p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})\right]\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)}),\end{array}\text{\hspace{1em}(15.3)}$$
• 对于分母，可以设 $q(z_{1:t}^{(i)}|x_{1:t}^{(i)})$ 满足 $q(z_{1:t}^{(i)}|x_{1:t}^{(i)})=q(z_t^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})\cdot q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})$.
• 将分子分母带入得，
  $$\begin{array}{rl}{w}_t^{(i)}& \propto \frac{p(z_{1:t}^{(i)}|x_{1:t}^{(i)})}{q(z_{1:t}^{(i)}|x_{1:t}^{(i)})}\\ & =\frac{\left[\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\cdot p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})\right]\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}{q(z_t^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})\cdot q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}\\ & =\underset{\text{Constant}}{\underbrace{\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}}}\cdot \left[\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})}{\underset{q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})}{\underbrace{q(z_t^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})}}}\right]\cdot \underset{w_{t-1}^{(i)}}{\underbrace{\frac{p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}{q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}}}\\ & \propto \left[\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})}{q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})}\right]\cdot w_{t-1}^{(i)},\end{array}\text{\hspace{1em}(15.4)}$$
  如此便建立了递推计算关系，可以通过蒙特卡洛的重要度采样方法计算filtering问题。

  为了采样能收敛于真正的期望，需要对重要度进行归一化，即保证 ${\sum }_{i=1}^N{w}_t^{(i)}=1$，这也就是为什么递推关系中可以化简常数，只要得到正比的关系即可。

  上述算法即为基础的粒子滤波算法，其中对状态值的每一次抽样结果 $z_t^{(i)}\sim q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$ 被称为一个粒子。

16.2.2. 重采样

  在应用上述粒子滤波算法的过程中，很多时候存在粒子退化问题，即只有少数状态拥有真正有效的重要度 $w_t^{(i)}$，其他大多数状态的重要度接近0，但算法每次都需要对所有 $z_t^{(i)}$ 以 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$ 的概率抽样，而 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$ 与重要度 $w_t^{(i)}$ 之间没有关联，这意味着可能有很多情况 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$ 较大，$w_t^{(i)}$ 较小，抽样结果充斥着不重要的信息，或是反过来，有效信息在抽样结果中占比很低，这些情况都会使计算资源浪费，模型难以收敛，效果变差等等。

  重采样的思想用于解决上述问题，具体来说，当使用概率 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$ 完成抽样后，计算每个粒子的重要度，并对重要度计算CDF，再在0-1的均匀分布上抽取N个分裂因子 ${\mu }_t^{(i)}$，每个 ${\mu }_t^{(i)}$ 落在CDF的哪个粒子区间内，就把哪个粒子复制一次（有时也叫粒子分裂）。

  例如，抽样得到的结果为z1，z2，z3（从小到大排列），三者的重要度分别为0.1,0.1,0.8，计算CDF为0-0.1、0.1-0.2、0.8-1.0三段，则需要抽取三个分裂因子，假设为0.15,0.38,0.54，第一个落入z2区间（0.1-0.2），则z2复制一次，另外两个落入z3区间（0.2-0.8），则z3复制两次，通过这个方法尽可能让每次抽样的粒子拥有接近1/N的重要度（或者说重要度高的粒子分裂的次数也多，更加能被有效利用，而重要度低的粒子可以适当舍弃）。

  当样本的维度过高时粒子退化问题将更加可怕，此时重采样的必要性更高；

16.2.3. Sampling-Importance-Resampling（SIR）

  除了重采样外，还有一种方法可以解决上文中提到的粒子退化问题，就是设置一个更加合理的 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})$，有一个方法是令 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)})$，此时重要度的递推关系为 $w_t^{(i)}=p(x_t^{(i)}|z_t^{(i)})\cdot w_{t-1}^{(i)}$，这种改进版的粒子滤波再加上重采样，可以进一步提升算法性能，叫做SIR.

  为什么要令 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)})$，其实是逻辑上的考虑，filtering的目标为计算 $p(z_t^{(i)}|x_{1:t}^{(i)})$（或 $p(z_{1:t}^{(i)}|x_{1:t}^{(i)})$），而这个问题中粒子滤波的方案是要对 $z_t^{(i)}$ 进行采样来计算均值，一个很简单的思想是，在已知 $x_t^{(i)}$ 的前提下，$p(x_t^{(i)}|z_t^{(i)})$ 越大的 $z_t^{(i)}$ 很多时候越有可能成为真正的状态，因此应该赋予其越高的权重，这就是要设置一个特殊的 $q(z_t^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)})$ 的原因，就是为了凑出 $w_t^{(i)}=p(x_t^{(i)}|z_t^{(i)})\cdot w_{t-1}^{(i)}$，然而这是一个朴素的思想，后人对这个建议分布的取法进行了很多研究，这里只介绍最基本的粒子滤波，因此不做赘述。