机器学习笔记之卡尔曼滤波(二)滤波思想的推导过程

引言

上一节介绍了动态模型，本节将介绍卡尔曼滤波模型中的滤波问题。

回顾：隐马尔可夫模型 VS 卡尔曼滤波

动态模型(Dynamic Model)的局部概率图模型表示如下：

这里并非单独比较隐马尔可夫模型和卡尔曼滤波两种模型，而是对描述隐马尔可夫模型和卡尔曼滤波的相关性质 进行比较：
(这里从模型参数的角度进行比较)
'非线性、非高斯动态模型'的代表(Non-Linear,Non-Gaussian Dynamic Model)——粒子滤波(Particle Filter)在后续介绍时再进行归纳,这里仅归纳2种模型。

离散状态动态模型(Discrete State Dynamic Model)

具有代表性的模型——隐马尔可夫模型。

• 状态转移概率$\mathcal{P}(i_t\mid i_{t-1})$
  由于离散状态动态模型中的隐变量是离散型随机变量，因此 $\mathcal{P}(i_t\mid i_{t-1})$通过查找状态转移矩阵$\mathcal{A}$得到对应结果：
  $$\begin{array}{rl}\mathcal{A}& =[a_{ij}{]}_{\mathcal{K}\times \mathcal{K}}={\left[\begin{array}{c}{a}_{11},a_{12},\cdots ,a_{1\mathcal{K}}\\ a_{21},a_{22},\cdots ,a_{2\mathcal{K}}\\ ⋮\\ a_{\mathcal{K}1},a_{\mathcal{K}2},\cdots ,a_{\mathcal{K}\mathcal{K}}\end{array}\right]}_{\mathcal{K}\times \mathcal{K}}\\ a_{ij}& =\mathcal{P}(i_t=q_j\mid i_{t-1}=q_i)\end{array}$$
  其中$q_i,q_j$均是 隐变量取值的离散集合$\mathcal{Q}$中的元素：
  q i , q j ∈ Q = { q 1 , q 2 , ⋯   , q K } qi,qj∈Q={q1,q2,⋯,qK} qi​,qj​∈Q={q1​,q2​,⋯,qK​}​
• 发射概率$\mathcal{P}(o_t\mid i_t)$
  离散状态动态模型中对观测变量 O = { o 1 , o 2 , ⋯   , o T } \mathcal O =\{o_1,o_2,\cdots,o_T\} O={o1​,o2​,⋯,oT​}没有具体要求，它可以是离散型随机变量，也可以是连续型随机变量。这里为容易表达起见，设定$\mathcal{O}$是离散型随机变量。因此 $\mathcal{P}(o_t\mid i_t)$通过查找发射矩阵$\mathcal{B}$得到相应结果：
  $$\begin{array}{rl}\mathcal{B}& =[b_j(k){]}_{\mathcal{K}\times \mathcal{M}}={\left[\begin{array}{c}{b}_1(1),b_1(2),\cdots ,b_1(\mathcal{M})\\ b_2(1),b_2(2),\cdots ,b_2(\mathcal{M})\\ ⋮\\ b_{\mathcal{K}}(1),b_{\mathcal{K}}(2),\cdots ,b_{\mathcal{K}}(\mathcal{M})\end{array}\right]}_{\mathcal{K}\times \mathcal{M}}\\ b_j(k)& =\mathcal{P}(o_t=v_k\mid i_t=q_j)\end{array}$$
  而$v_k$表示 观测变量取值的离散集合$\mathcal{V}$中的元素：
  v k ∈ V = { v 1 , v 2 , ⋯   , v M } v_k \in \mathcal V = \{v_1,v_2,\cdots,v_{\mathcal M}\} vk​∈V={v1​,v2​,⋯,vM​}
• 初始概率$\mathcal{P}(i_1)$
  在隐马尔可夫模型中介绍过，初始概率分布使用$\pi$进行表示：
  $$\mathcal{P}(i_1)=\pi$$
  综上，离散状态动态模型需要求解的模型参数具体表示如下：
  $$\lambda =(\pi ,\mathcal{A},\mathcal{B})$$

线性高斯动态模型(Linear Gaussian Dynamic Model)

具有代表性的模型——卡尔曼滤波。
相比于离散状态动态模型，该模型更突出的是线性：隐变量与观测变量均是连续型随机变量。

• 状态转移概率$\mathcal{P}(i_t\mid i_{t-1})$
  线性高斯动态模型中隐变量之间服从线性关系，且对应噪声服从高斯分布：
  $$\begin{array}{r}{i}_t=\mathcal{A}\cdot i_{t-1}+\mathcal{B}+ϵϵ\sim \mathcal{N}(0,\mathcal{Q})\\ \mathcal{P}(i_t\mid i_{t-1})\sim \mathcal{N}(\mathcal{A}\cdot i_{t-1}+\mathcal{B},\mathcal{Q})\end{array}$$
  其中$\mathcal{A},\mathcal{B}$表示线性关系的模型参数；$\mathcal{Q}$表示转移过程高斯分布噪声的协方差信息。
• 发射概率$\mathcal{P}(o_t\mid i_t)$
  同理，隐变量与观测变量之间同样服从线性关系，对应噪声服从高斯分布：
  $$\begin{array}{r}{o}_t=\mathcal{C}\cdot i_t+\mathcal{D}+\delta \delta \sim \mathcal{N}(0,\mathcal{R})\\ \mathcal{P}(o_t\mid i_t)\sim \mathcal{N}(\mathcal{C}\cdot i_t+\mathcal{D},\mathcal{R})\end{array}$$
  这里的$\mathcal{C},\mathcal{D}$表示线性关系的模型参数；$\mathcal{R}$表示发射过程噪声高斯分布的协方差信息。
• 初始概率$\mathcal{P}(i_1)$
  不同于$\pi$这种具体的概率值结果，线性高斯动态模型的初始概率同样是高斯分布：
  $$\mathcal{P}(i_1)\sim \mathcal{N}({\mu }_1,{\Sigma }_1)$$
  综上，线性高斯动态模型需要求解的模型参数表示如下：
  $$\lambda =(\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{D},\mathcal{Q},\mathcal{R},{\mu }_1,{\Sigma }_1)$$

滤波问题思想推导

公式推导过程

我们需要解决的滤波问题具体表示如下：
$$\mathcal{P}(i_t\mid o_t,o_{t-1},\cdots ,o_1)$$
类似于求值问题$\mathcal{P}(\mathcal{O}\mid \lambda )$，我们希望通过迭代方式表示 $t$时刻滤波结果与其他时刻滤波结果之间的关联关系：

• 首先，滤波问题本身是一个条件概率。根据条件概率的定义，改写为如下形式：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)=\frac{\mathcal{P}(i_t,o_1,\cdots ,o_t)}{\mathcal{P}(o_1,\cdots ,o_t)}$$
  由于$\mathcal{P}(o_1,\cdots ,o_t)$是初始时刻到$t$时刻观测变量的联合概率分布，而观测变量是给定的数据集合，因此$\mathcal{P}(o_1,\cdots ,o_t)$是可求的。令${\mathcal{C}}_1=\mathcal{P}(o_1,\cdots ,o_t)$，则有：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)=\frac{1}{{\mathcal{C}}_1}\mathcal{P}(i_t,o_1,\cdots ,o_t)$$
• 将联合概率分布使用条件概率公式展开，展开为$o_t$作为后验的条件概率的乘积形式：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)=\frac{1}{{\mathcal{C}}_1}\left[\mathcal{P}(o_t\mid o_1,\cdots ,o_{t-1},i_t)\cdot \mathcal{P}(o_1,\cdots ,o_{t-1},i_t)\right]$$
  观察中括号中的第一项，可以使用观测独立性假设改写成如下形式：
  观测独立性假设是‘隐马尔可夫模型’中介绍的，需要的去复习一下~隐马尔可夫模型介绍-传送门
  $$\mathcal{P}(o_t\mid o_1,\cdots ,o_{t-1},i_t)=\mathcal{P}(o_t\mid i_t)$$
  从而最终改写成如下形式：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)=\frac{1}{{\mathcal{C}}_1}\left[\mathcal{P}(o_t\mid i_t)\cdot \mathcal{P}(o_1,\cdots ,o_{t-1},i_t)\right]$$
• 将括号中的$\mathcal{P}(o_1,\cdots ,o_{t-1},i_t)$通过条件概率公式，展开成以$i_t$为后验的条件概率乘积形式：
  $$\mathcal{P}(o_1,\cdots ,o_{t-1},i_t)=\mathcal{P}(i_t\mid o_1,\cdots ,o_{t-1})\cdot \mathcal{P}(o_1,\cdots ,o_{t-1})$$
  其中$\mathcal{P}(o_1,\cdots ,o_{t-1})$同样也是观测变量的联合概率分布，是可求的。因此定义${\mathcal{C}}_2=\mathcal{P}(o_1,\cdots ,o_{t-1})$，从而有：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)=\frac{{\mathcal{C}}_2}{{\mathcal{C}}_1}\left[\mathcal{P}(o_t\mid i_t)\cdot \mathcal{P}(i_t\mid o_1,\cdots ,o_{t-1})\right]$$
• 继续观察中括号中的第二项：这明显是一个预测问题，通常采用方法是 通过积分，引入隐变量$i_{t-1}$：
  其中红色框表示‘条件项’与‘后验项’;需要引入中间变量(蓝色框)将它们关联起来。
不要忘记，隐变量$i_t$是连续型随机变量，其对应的积分是${\int }_{i_t}$。
  
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_{t-1})={\int }_{i_{t-1}}\mathcal{P}(i_t,i_{t-1}\mid o_1,\cdots ,o_{t-1})d{i}_{t-1}$$
  再根据条件概率的推导式，转化为如下格式：
  $${\int }_{i_{t-1}}\mathcal{P}(i_t\mid i_{t-1},o_1,\cdots ,o_{t-1})\cdot \mathcal{P}(i_{t-1}\mid o_1,\cdots ,o_{t-1})d{i}_{t-1}$$
  其中第一项使用齐次马尔可夫假设将其简化为$\mathcal{P}(i_t\mid i_{t-1})$；第二项正是$t-1$时刻的滤波问题。

至此，$t$时刻与$t-1$时刻滤波结果的关联关系。最终结果整理如下：
$$\begin{array}{rl}\mathcal{P}(i_t\mid o_1,\cdots ,o_t)& =\frac{{\mathcal{C}}_2}{{\mathcal{C}}_1}\left[\mathcal{P}(o_t\mid i_t){\int }_{i_{t-1}}\mathcal{P}(i_t\mid i_{t-1})\cdot \mathcal{P}(i_{t-1}\mid o_1,\cdots ,o_{t-1})d{i}_{t_1}\right]\\ & =\frac{\mathcal{P}(o_1,\cdots ,o_{t-1})}{\mathcal{P}(o_1,\cdots ,o_t)}\cdot \left[\mathcal{P}(o_t\mid i_t){\int }_{i_{t-1}}\mathcal{P}(i_t\mid i_{t-1})\cdot \mathcal{P}(i_{t-1}\mid o_1,\cdots ,o_{t-1})d{i}_{t-1}\right]\end{array}$$

滤波问题求解步骤

卡尔曼滤波在处理滤波问题时，是使用在线算法(On-line Algorithm)。即 执行到某时刻时，才能够计算出该时刻关于隐变量的后验信息。

它的求解步骤是一个迭代过程。每一次迭代均包含2个步骤：

• 更新步骤(Update)：根据给定的观测变量结果(从初始时刻 $to$ 当前时刻)，求出当前时刻隐变量的后验概率分布：
  $$\mathcal{P}(i_t\mid o_1,\cdots ,o_t)$$
• 预测步骤(Prediction)：根据给定的观测变量结果(从初始时刻 $to$ 当前时刻)，求出下一时刻隐变量的后验概率分布：
  $$\mathcal{P}(i_{t+1}\mid o_1,\cdots ,o_t)$$

具体过程表示如下：
已知条件：

• 隐变量初始时刻的概率分布$\mathcal{P}(i_1)$：
  $$\mathcal{P}(i_1)\sim \mathcal{N}({\mu }_1,{\Sigma }_1)$$
• 基于观测独立性假设，观测变量$o_t$在给定对应时刻隐变量$i_t$的条件概率 $\mathcal{P}(o_t\mid i_t)$：
  $$\mathcal{P}(o_t\mid i_t)\sim \mathcal{N}(\mathcal{C}\cdot i_t+\mathcal{D},\mathcal{R})$$
• 基于齐次马尔可夫假设，隐变量$i_t$在给定上一时刻隐变量$i_{t-1}$的条件概率 $\mathcal{P}(i_t\mid i_{t-1})$：
  $$\mathcal{P}(i_t\mid i_{t-1})\sim \mathcal{N}(\mathcal{A}\cdot i_{t-1}+\mathcal{B},\mathcal{Q})$$

相关公式介绍：
本质上，卡尔曼滤波待求解的模型参数有很多：
$$\lambda =(\mathcal{A},\mathcal{B},\mathcal{C},\mathcal{D},\mathcal{Q},\mathcal{R},{\mu }_1,{\Sigma }_1)$$
但实际上，这些参数都是用来描述正态分布 的参数。因此这里给出出现条件概率、积分情况下 概率分布的变化：
这属于高斯分布的常用计算公式范畴，给大家推荐一篇相关推导文章。PRML笔记-高斯分布-传送门

给定 变量$\mathcal{X}$的边缘概率分布$\mathcal{P}(\mathcal{X})$ 与 给定$\mathcal{X}$条件下，变量$\mathcal{Y}$的条件概率分布$\mathcal{P}(\mathcal{Y}\mid \mathcal{X})$ 如下：
这里假设'协方差矩阵'$\Lambda ,\mathcal{L}$是‘正定矩阵’,它们均可以求逆。
$$\{\begin{array}{l}\mathcal{P}(\mathcal{X}):x\in \mathcal{P}(\mathcal{X}),x\sim \mathcal{N}(\mu ,{\Lambda }^{-1})\\ \mathcal{P}(\mathcal{Y}\mid \mathcal{X}):y\in \mathcal{P}(\mathcal{Y}\mid \mathcal{X}),y\sim \mathcal{N}(\mathcal{A}\cdot x+\mathcal{B},{\mathcal{L}}^{-1})\end{array}$$
则变量$\mathcal{Y}$的边缘概率分布$\mathcal{P}(\mathcal{Y})$可表示为：
$$\begin{array}{rl}\mathcal{P}(\mathcal{Y})& ={\int }_{\mathcal{X}}\mathcal{P}(\mathcal{X},\mathcal{Y})d\mathcal{X}={\int }_{\mathcal{X}}\mathcal{P}(\mathcal{X})\cdot \mathcal{P}(\mathcal{Y}\mid \mathcal{X})d\mathcal{X}\\ & \to \mathcal{P}(\mathcal{Y}):y\in \mathcal{P}(\mathcal{Y}),y\sim \mathcal{N}(\mathcal{A}\cdot \mu +\mathcal{B},{\mathcal{L}}^{-1}+\mathcal{A}{\Lambda }^{-1}{\mathcal{A}}^T)\end{array}$$
给定$\mathcal{Y}$条件下，变量$\mathcal{X}$的条件概率分布$\mathcal{P}(\mathcal{X}\mid \mathcal{Y})$可表示为：
P ( X ∣ Y ) = P ( Y ∣ X ) ⋅ P ( X ) P ( Y ) → P ( X ∣ Y ) : x ∈ P ( X ∣ Y ) , x ∼ N ( Σ { A T L ( y − B ) + A μ } , Σ ) Σ = Λ + A T L A − 1 P(X∣Y)=P(Y∣X)⋅P(X)P(Y)→P(X∣Y):x∈P(X∣Y),x∼N(Σ{ATL(y−B)+Aμ},Σ)Σ=Λ+ATLA−1 P(X∣Y)​=P(Y)P(Y∣X)⋅P(X)​→P(X∣Y):x∈P(X∣Y),x∼N(Σ{ATL(y−B)+Aμ},Σ)Σ=Λ+ATLA−1​

具体过程：

• 初始步骤$(t=1)$：
  • 隐变量$i_1$的初始化作为$i_1$的 更新步骤(Update)：
    $$\mathcal{P}(i_1\mid o_1)=\mathcal{P}(i_1)\sim \mathcal{N}({\mu }_1,{\Sigma }_1)$$
  • 预测步骤：基于$\mathcal{P}(i_1\mid o_1)$，求解下一时刻隐变量$i_2$基于$o_1$的条件概率$\mathcal{P}(i_2\mid o_1)$：
    $$\begin{array}{r}\mathcal{P}(i_2\mid o_1)={\int }_{i_1}\mathcal{P}(i_2\mid i_1)\cdot \mathcal{P}(i_1\mid o_1)\end{array}$$
    根据已知条件，$\mathcal{P}(i_2\mid i_1)\sim \mathcal{N}(\mathcal{A}\cdot i_1+\mathcal{B},\mathcal{Q})$，结合初始概率分布，$\mathcal{P}(i_2\mid o_1)$的概率分布表示如下：
    $$\mathcal{P}(i_2\mid o_1)\sim \mathcal{N}(\mathcal{A}\cdot {\mu }_1+\mathcal{B},\mathcal{Q}+\mathcal{A}{\Sigma }_1{\mathcal{A}}^T)$$
    至此，我们通过预测步骤求解出$\mathcal{P}(i_2\mid o_1)$的概率分布：
    注意，这仅是一个‘预测结果’。
    $$\{\begin{array}{l}{\mu }_2=\mathcal{A}\cdot {\mu }_1+\mathcal{B}\\ {\Sigma }_2=\mathcal{Q}+\mathcal{A}{\Sigma }_1{\mathcal{A}}^T\end{array}$$
• $t=2$时刻：
  通过 预测步骤 得到了$\mathcal{P}(i_2\mid o_1)$的概率分布结果，结合观测独立性概率$\mathcal{P}(o_2\mid i_2)$，求解$\mathcal{P}(i_2\mid o_1,o_2)$。
  • 更新步骤：
    已知：
    $$\begin{array}{rl}\mathcal{P}(o_2\mid i_2)& \sim \mathcal{N}(\mathcal{C}\cdot i_2+\mathcal{D},\mathcal{R})\\ \mathcal{P}(i_2\mid o_1)& \sim \mathcal{N}({\mu }_2,{\Sigma }_2)\end{array}$$
    则有：
    $$\begin{array}{rl}\mathcal{P}(i_2\mid o_1,o_2)& =\frac{\mathcal{P}(o_1)}{\mathcal{P}(o_1,o_2)}\left[\mathcal{P}(o_2\mid i_2)\cdot \mathcal{P}(i_2\mid o_1)\right]\\ & \propto \mathcal{P}(o_2\mid i_2)\cdot \mathcal{P}(i_2\mid o_1)\end{array}$$
    将结果代入上述公式中：
    P ( i 2 ∣ o 1 , o 2 ) ∼ N ( Σ { C T R − 1 ( o 2 − D ) + C ⋅ μ 2 } , Σ ) Σ = Σ 2 − 1 + C R − 1 C − 1 { μ 2 ∗ = ( Σ 2 − 1 + C R − 1 C − 1 ) { C T R − 1 ( o 2 − D ) + C ⋅ μ 2 } Σ 2 ∗ = Σ 2 − 1 + C R − 1 C − 1 P(i2∣o1,o2)∼N(Σ{CTR−1(o2−D)+C⋅μ2},Σ)Σ=Σ−12+CR−1C−1{μ∗2=(Σ−12+CR−1C−1){CTR−1(o2−D)+C⋅μ2}Σ∗2=Σ−12+CR−1C−1 P(i2​∣o1​,o2​)​∼N(Σ{CTR−1(o2​−D)+C⋅μ2​},Σ)Σ=Σ2−1​+CR−1C−1{μ2∗​=(Σ2−1​+CR−1C−1){CTR−1(o2​−D)+C⋅μ2​}Σ2∗​=Σ2−1​+CR−1C−1​​
• 后续时刻以此类推。

迭代过程总结

观察上面的求解步骤，它明显包含两个步骤：

• 预测步骤(Prediction)：对下一时刻的隐变量进行一个预测：
  $$\begin{gathered} \mathcal{P}(i_2\mid o_1)\sim \mathcal{N}({\mu }_2,{\Sigma }_2) \\ \{\begin{array}{l}{\mu }_2=\mathcal{A}\cdot {\mu }_1+\mathcal{B}\\ {\Sigma }_2=\mathcal{Q}+\mathcal{A}{\Sigma }_1{\mathcal{A}}^T\end{array} \end{gathered}$$
• 更新步骤(Update)：在上一时刻预测的基础上，对当前时刻隐变量进行更新。同时对下一时刻隐变量进行预测。
  P ( i 2 ∣ o 1 , o 2 ) ∼ N ( Σ { C T R − 1 ( o 2 − D ) + C ⋅ μ 2 } , Σ ) Σ = Σ 2 − 1 + C R − 1 C − 1 { μ 2 ∗ = ( Σ 2 − 1 + C R − 1 C − 1 ) { C T R − 1 ( o 2 − D ) + C ⋅ μ 2 } Σ 2 ∗ = Σ 2 − 1 + C R − 1 C − 1 P(i2∣o1,o2)∼N(Σ{CTR−1(o2−D)+C⋅μ2},Σ)Σ=Σ−12+CR−1C−1{μ∗2=(Σ−12+CR−1C−1){CTR−1(o2−D)+C⋅μ2}Σ∗2=Σ−12+CR−1C−1 P(i2​∣o1​,o2​)​∼N(Σ{CTR−1(o2​−D)+C⋅μ2​},Σ)Σ=Σ2−1​+CR−1C−1{μ2∗​=(Σ2−1​+CR−1C−1){CTR−1(o2​−D)+C⋅μ2​}Σ2∗​=Σ2−1​+CR−1C−1​​
• 重复执行上述两个步骤。
  $({\mu }_2,{\Sigma }_2)\to ({\mu }_2^{\ast },{\Sigma }_2^{\ast })$，是一个明显的“先预测，再对预测修正”的过程。

至此，卡尔曼滤波部分介绍结束，下一节将介绍粒子滤波(Particle Filter)。

相关参考：
【PRML】高斯分布
机器学习-白板推导系列(十五)-线性动态系统-卡曼滤波（Kalman Filter）笔记
机器学习-线性动态系统2-Filtering问题
机器学习-线性动态系统3-Filtering问题求解