粒子滤波算法

本文为《蒙特卡罗方法理论和应用》笔记

$\text{1.动态空间模型}$

动态空间模型包括过程模型和观测模型，分别由状态方程和观测方程描述，状态方程和观测方程为：
$$\begin{array}{rl}\mathbf{x}(t)& =a(\mathbf{x}(t-1),\mathbf{u}(t))\\ \mathbf{y}(t)& =b(\mathbf{x}(t),\mathbf{v}(t))\end{array}$$
当过程噪音和观测噪音可分离时，状态方程和观测方程为：
$$\begin{array}{rl}{\mathbf{x}}_k& =a\left({\mathbf{x}}_{k-1},{\mathbf{u}}_k\right)=a\left({\mathbf{x}}_{k-1}\right)+{\mathbf{u}}_k\\ {\mathbf{y}}_k& =b\left({\mathbf{x}}_k,{\mathbf{v}}_k\right)=b\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k\end{array}$$
式中，所有下标$k$表示时刻；$a({\mathbf{x}}_{k-1},{\mathbf{u}}_k)$和$a({\mathbf{x}}_{k-1})$为过程函数； $b({\mathbf{x}}_k,{\mathbf{v}}_k)$和$b({\mathbf{x}}_k)$为观测函数；${\mathbf{x}}_k$和${\mathbf{x}}_{k-1}$为状态，${\mathbf{y}}_k$为观测；${\mathbf{u}}_k$为过程噪声； ${\mathbf{v}}_k$为观测噪声。过程噪声的统计参数为${\mathbf{Q}}_k$ ，观测噪声的统计参数为 ${\mathbf{R}}_k$ ，对单变量模型${\mathbf{Q}}_k$ 和 ${\mathbf{R}}_k$表示方差，对多变量模型${\mathbf{Q}}_k$ 和 ${\mathbf{R}}_k$表示协方差。可用统计参数${\mathbf{Q}}_k$ 和 ${\mathbf{R}}_k$唯一地描述高斯噪声，但不能唯一地描述非高斯噪声。

$\text{2.隐马尔可夫随机过程}$

序贯随机过程的下一状态不仅与当前状态有关，而且还与以前所有状态都有关。序贯随机过程的概率分布为：
$$f(\mathbf{x})=f\left({\mathbf{x}}_1|{\mathbf{x}}_{1:0}\right)f\left({\mathbf{x}}_2|{\mathbf{x}}_{1:1}\right)f\left({\mathbf{x}}_3|{\mathbf{x}}_{1:2}\right)\cdots f\left({\mathbf{x}}_k|{\mathbf{x}}_{1:k-1}\right)\cdots f\left({\mathbf{x}}_m|{\mathbf{x}}_{1:m-1}\right)$$
隐马尔可夫随机过程的概率分布为：
$$f(\mathbf{x})=f\left({\mathbf{x}}_1|{\mathbf{x}}_0\right)f\left({\mathbf{x}}_2|{\mathbf{x}}_1\right)\cdots f\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1}\right)\cdots f\left({\mathbf{x}}_m|{\mathbf{x}}_{m-1}\right)$$
隐马尔可夫随机过程本身具有马尔可夫性，即无后效性，下一状态只与当前状态有关，与以前状态无关。${\mathbf{x}}_{k+1}$只与${\mathbf{x}}_k$ 有关，与${\mathbf{x}}_{0:k-1}$无关，${\mathbf{y}}_k$只与${\mathbf{x}}_k$有关，状态$\mathbf{x}$接受的过程噪声$\mathbf{u}$，观测$\mathbf{y}$观测噪声$\mathbf{v}$.

$\text{3.贝叶斯递推滤波}$

状态变量$\mathbf{X}$ 的先验概率分布为$f(\mathbf{x})$ ，观测变量$\mathbf{Y}$ 的概率分布为 $f(\mathbf{y}|\mathbf{x})$
故状态变量$\mathbf{X}$ 的后验概率分布为
$$f(\mathbf{x}|\mathbf{y})=\frac{f(\mathbf{y}|\mathbf{x})f(\mathbf{x})}{\int f(\mathbf{y}|\mathbf{x})f(\mathbf{x})d\mathbf{x}}$$
状态变量的估计值为
$$\hat{\mathbf{x}}=\int \mathbf{x}f(\mathbf{x}|\mathbf{y})d\mathbf{x}$$
如果后验概率分布$f(\mathbf{x}|\mathbf{y})$的分母积分和状态变量估计值的积分都能求解，则可得到状态变量的估计值，称为$\text{贝叶斯估计}$.

$\text{4.贝叶斯递推滤波描述}$

$0:k$时刻的状态值用${\mathbf{x}}_{0:k}$表示，
${\mathbf{x}}_{0:k}={\mathbf{x}}_0,{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_k$;
$1:k$时刻的观测值用${\mathbf{y}}_{1:k}$表示，
${\mathbf{y}}_{1:k}={\mathbf{y}}_1,{\mathbf{y}}_2,\cdots ,{\mathbf{y}}_k$.
贝叶斯递推滤波是给定所有观测值${\mathbf{y}}_{1:k}$构造状态${\mathbf{x}}_{0:k}$的后验概率分布，贝叶斯递推滤波分为$\underline{预测阶段}$和$\underline{更新阶段}$.

$\text{预测阶段}$，预测阶段主要基于${\mathbf{x}}_{k-1}$预测${\mathbf{x}}_k$ 预测方程为：$$f\left({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k-1}\right)=\int f\left({\mathbf{x}}_{0:k}|{\mathbf{x}}_{0:k-1}\right)f\left({\mathbf{x}}_{0:k-1}|{\mathbf{y}}_{1:k-1}\right)\mathrm{d}{\mathbf{x}}_{0:k-1}$$

$\text{更新阶段}$，更新阶段主要根据${\mathbf{y}}_k$更新${\mathbf{x}}_k$. 根据贝叶斯定理，后验概率分布为
$$f\left({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k}\right)=f\left({\mathbf{y}}_{1:k}|{\mathbf{x}}_{0:k}\right)f\left({\mathbf{x}}_{0:k}\right)/f\left({\mathbf{y}}_{1:k}\right)$$
式中， f ( y 1 : k ) = ∫ f ( y 1 : k ∣ x 0 : k ) f ( x 0 :