贝叶斯与卡尔曼滤波(3)--随机过程

前面的内容中，一般都是先验概率X和一个观测概率Y，求它的后验概率

而在实际中，我们更多的面对的是一个随机过程，它有不止一个随机变量\(X_0, X_1,X_2, ..., X_n\),以及不止一个观测值

有一个初值\(X_0\),同时又\(K\)个观测值\(y_1, y_2,...,y_k\),现在该怎么处理呢？

  1. 所有的\(X_0\backsim X_n\)的先验概率我们都靠猜来赋予初值。

    这个方法是可行的，但是它有缺点：过于依赖观测，等于是放弃了预测信息

  2. 至于\(X_0\)的概率是猜的，\(X_1\backsim X_n\)的先验概率是递推的。

    这个方法明显就比第一个方法要好。贝叶斯滤波实际上用的也是这个方法。

怎么做呢？一般使用状态方程。状态方程用来描述\(X_k\)与\(X_{k-1}\)之间的关系

比如：

\[X_k= \frac{1}{2}gt^2+Q
\]

\[\begin{equation}
\begin{aligned}
X_k&
= X_{k-1}+\dot{X}_{k-1}(t_k-t_{k-1})+Q\\&
=X_{k-1}+gt(t_k-t_{k-1})+Q

\end{aligned}
\end{equation}
\]

上式中，是用了泰勒展开，只展开一阶。

但是很遗憾，在实际中，大量的问题是无法写出上面这种状态方程的。一段随机波动的曲线，我们根本不知道它服从什么规律。这个时候我们可以随便写一个方程，比如\(X_k=X_{k-1}+Q\)，\(Q\backsim N(0, 1000)\)，状态方程可以建立的非常粗糙，这没关系。

  1. 状态方程：反映了\(X_k\)与\(X_{k-1}\)之间的关系：

    \[X_k=f(X_{k-1})+Q_k
    \]

    \(Q\)是预测噪声。

  2. 观测方程：反映了状态是如何引起传感器的读数

    比如测温度，状态：温度， 观测：温度，那么观测方程\(Y_K = X_k+R_k\)

    再比如，状态：位移， 观测：角度，那么观测方程可以是\(Y_k = arcsinX_k + R_k\)

    一般来说，观测方程为：

    \[Y_k = h(X_k)+R_k
    \]

    \(R_k\)是观测噪声

这样得到两个方程：

\[\begin{cases}
X_k=f(X_{k-1})+Q_k\\
Y_k = h(X_k)+R_k
\end{cases}
\]

那怎么递推呢？假设有一个随机过程\(X_k=2X_{k-1}\),假设无\(Q_k\),即预测完全准确。也没有观测。首先猜一个初值，\(X_0\backsim N(0, 1)\)，那么：

\[\begin{cases}
X_1\backsim N(0, 2^2)\\
X_2\backsim N(0, 2^4)\\
X_3\backsim N(0, 2^8)\\
...
\end{cases}
\]

可以看到，方差会越来越大，这很明显是不可行的。所以引入观测是十分有必要的。

下面我们来看一下真正的递推是什么样的。

\[X_0 \backsim N(0, 1) \xrightarrow[]{预测步} X_1^{-}\backsim N(0, 2^2) \xrightarrow[\uArr 观测y_1=0]{更新步/后验步} X_1^{+} \backsim N(0, 0.8) \xrightarrow[]{预测步} X_2^{-}\backsim N(0, 3.2) \xrightarrow[\uArr 观测y_2=0.1]{更新步} X_2^{+}\xrightarrow[]{预测步}

\]

以上就是完整的递推过程，可以总结为：

\[先验\xrightarrow[]{预测步}先验\xrightarrow[\uArr 观测]{更新步} 后验/先验\xrightarrow[]{预测步}先验\xrightarrow[\uArr 观测]{更新步}后验/先验 \xrightarrow[]{}
\]

所以如果没有观测，上一步的先验概率作为下一步的先验概率，方差会越来越大。有观测的话，上一步的后验概率作为下一步的先验概率，方差会一大一小，比较稳定。

那么具体该怎么做呢？

贝叶斯滤波算法的推导。

已知：

\[\begin{cases}
X_k=f(X_{k-1})+Q_k\\
Y_k = h(X_k)+R_k
\end{cases}
\]

其中：\(X_k, X_{k-1}, Y_k, Q_k, R_k\)都是随机变量

假设：\(X_0, Q_1,...Q_k,R_1,..., R_k相互独立\)。有观测值\(y_1, y_2,...,y_k\)。设初值\(X_0\)的概率密度为\(f_0(x)\), 它是一个定值。\(Q_k\)的概率密度为\(f_{Q_k}(x)\),\(R_k\)的概率密度为\(f_{R_k}(x)\)

此处引申出一个重要定理：条件概率里的条件可以做逻辑推导。

比如：

\[\begin{equation}
\begin{aligned}
P(X=1|Y=2,Z=3)&
= P(X+Y=3|Y=2,Z=3)=P(X+Y=3|Y=2,Z-Y=1)\\&
\not = P(X+Y=3|Z=3)

\end{aligned}
\end{equation}
\]

下面开始正式的推导。

首先预测步，利用化积分为求和

\[P(X_1<x) = \sum_{u=-\infty}^{x}P(X_1=u)
\]

根据全概率公式：

\[\begin{equation}
\begin{aligned}
P(X_1=u)&
=\sum_{v=-\infty}^{+\infty}P(X_1=u|X_0=v)P(X_0=v)\\&
=\sum_{v=-\infty}^{+\infty}P(X_1-f(X_0)=u-f(X_0)|X_0=v)P(X_0=v)\\&
=\sum_{v=-\infty}^{+\infty}P(Q_1=u-f(X_0)|X_0=v)P(X_0=v) \qquad\qquad 此处(X_1-f(X_0)=Q_1)

\end{aligned}
\end{equation}
\]

观察这个公式，\(Q_1\)和\(X_0\)是相互独立的，那么公式中的后面可以拿掉了。但是再往后递推的时候，\(Q_k\)和\(X_{k-1}\)是不是相互独立呢？答案是肯定的，他们肯定是相互独立的，后面我们再证明。继续化简：

\[\begin{equation}
\begin{aligned}
P(X_1=u)&
=\sum_{v=-\infty}^{+\infty}P(Q_1=u-f(X_0))P(X_0=v)\\&
=\lim_{\epsilon \to 0}\sum_{v=-\infty}^{+\infty}f_{Q_1}[u-f(v)]\cdot\epsilon\cdot f_0(v)
\cdot \epsilon \\&
=\lim_{\epsilon \to 0}\int_{-\infty}^{+\infty}f_{Q_1}[u-f(v)]\cdot f_0(v)dv
\cdot \epsilon
\end{aligned}
\end{equation}
\]

所以：

\[\begin{equation}
\begin{aligned}
P(X_1<x) &
=\sum_{u=-\infty}^{x}P(X_1=u)\\&
=\sum_{u=-\infty}^{x} \lim_{\epsilon \to 0}\int_{-\infty}^{+\infty}f_{Q_1}[u-f(v)]\cdot f_0(v)dv
\cdot \epsilon \\&
= \int_{-\infty}^{x}\int_{-\infty}^{+\infty}f_{Q_1}[u-f(v)]\cdot f_0(v)dvdu
\cdot \epsilon
\end{aligned}
\end{equation}
\]

所以先验概率密度：

\[\begin{equation}
\begin{aligned}
f_1^-(x)&
=\frac{d}{dx}(P(X_1<x))\\&
=\int_{-\infty}^{+\infty}f_{Q_1}[x-f(v)]f_0(v)dv

\end{aligned}
\end{equation}
\]

然后再看更新步，观测\(Y_1=y_1\)

上面说过，似然概率可以这么写：

\[\begin{equation}
\begin{aligned}
f(Y_1|X_1)(y_1|x)&
=\lim_{\epsilon \to 0}\frac{P(y_1<Y_1<y_1+\epsilon|X_1=x)}{\epsilon}\\&

=\lim_{\epsilon \to 0}\frac{P(y_1 - h(x)<Y_1-h(x)<y_1-h(x)+\epsilon|X_1=x)}{\epsilon}\\&
=\lim_{\epsilon\to 0}\frac{P(y_1-h(x)<R_1<y_1-h(x)+\epsilon|X_1=\epsilon)}{\epsilon}

\end{aligned}
\end{equation}
\]

这个时候，有个问题，\(R_1\)与\(X_1\)是相互独立的吗？他们是相互独立的，但是如何证明呢？后面再证明。继续化简：

\[\begin{equation}
\begin{aligned}
f(Y_1|X_1)(y_1|x)&

=\lim_{\epsilon\to 0}\frac{P(y_1-h(x)<R_1<y_1-h(x)+\epsilon)}{\epsilon}\\&
=f_{R_1}[y_1-h(x)]

\end{aligned}
\end{equation}
\]

根据贝叶斯公式：

\[f_1^+(x)=\eta\cdot f_{R_1}[y_1-h(x)]\cdot f_1^-(x)
\]

\[\eta = (\int_{-\infty}^{+\infty}f_{R_1}[y_1-h(x)]\cdot f_1^-(x)dx)^{-1}
\]

\[\begin{aligned} f_{0}(x) & \stackrel{\text { 预测 }}{\Longrightarrow} f_{1}^{-}(x)=\int_{-\infty}^{+\infty} f_{Q_{1}}[x-f(v)] f_{0}(v) d v \stackrel{\text { 观测更新 }}{\Longrightarrow} f_{1}^{+}(x)=\eta_{1} * f_{R_{1}}\left[y_{1}-h(x)\right] * f_{1}^{-}(x) \\ & \stackrel{\text { 预测 }}{\Longrightarrow} f_{2}^{-}(x)=\int_{-\infty}^{+\infty} f_{Q_{2}}[x-f(v)] f_{1}^{+}(x) d v \stackrel{\text { 观测更新 }}{\Longrightarrow} f_{2}^{+}(x)=\eta_{2} * f_{R_{2}}\left[y_{2}-h(x)\right] * f_{2}^{-}(x) \\ & \cdots \cdots \\ & \stackrel{\text { 预测 }}{\Longrightarrow} f_{k}^{-}(x)=\int_{-\infty}^{+\infty} f_{Q_{k}}[x-f(v)] f_{k-1}^{+}(x) d v \stackrel{\text { 观测更新 }}{\Longrightarrow} f_{k}^{+}(x)=\eta_{k} * f_{R_{k}}\left[y_{k}-h(x)\right] * f_{k}^{-}(x) \end{aligned}
\]

接下来证明\(Q_k与X_{K-1}\)独立，\(X_k与R_k\)独立。

因为：

\[X_{K-1}=f(X_{k-2})+Q_{k-1}=F_0(X_{k-2}, Q_{k-1})
\]

\[X_{K-2}=F_1(X_{k-3}, Q_{k-2})
\]

\[...
\]

\[X_1 =F_{k-1}(X_0, Q_1)
\]

所以\(X_k\)他就是一个关于\(X_0, Q_1, Q_2, ..., Q_k\)的函数：

\[X_k = F_K(X_0, Q_1, Q_2, ..., Q_k)
\]

\[X_{k-1} = F_{K-1}(X_0, Q_1, Q_2, ..., Q_{k-1})
\]

因为\(X_0, Q_1, Q_2, ..., Q_k， R_1, R_2, ...,R_k\)相互独立，所以

  • \(X_{k-1} = F_{K-1}(X_0, Q_1, Q_2, ..., Q_{k-1})\)与\(Q_k\)独立
  • \(X_k = F_K(X_0, Q_1, Q_2, ..., Q_k)\)与\(R_k\)独立

完整算法：

  1. 设置初值：\(X_0\)以及\(f_0(x)\)

  2. 预测步：

    \[f_k^{-}(x)=\int_{-\infty}^{+\infty}f_{Q_k}[x-f(v)]f_{k-1}^{+}(v)dv
    \]

  3. 更新步

    \[f_k^+(x)=\eta_k \cdot f_{R_k}[y_k-h(x)]\cdot f_k^-(x)
    \]

    \[\eta_k = (\int_{-\infty}^{+\infty}f_{R_k}[y_1-h(x)]\cdot f_k^-(x)dx)^{-1}
    \]

我们已经算出了后验概率密度函数，那么如何根据这个概率密度函数算出具体的数值呢？这其实就是一个期望的问题。根据贝叶斯公式，可以得到：

\[\hat x_k^{+}=\int_{-\infty}^{+\infty}xf_k^{+}(x)dx
\]

贝叶斯滤波的缺点：

\(f_{k-1}^{+}\to f_k^{-}\)、\(\eta\)、\(\hat x_k^{+}\)的计算都需要无穷积分，大多数情况下，是没有解析解的。而贝叶斯滤波的每一步几乎都需要\(f_{k-1}^{+}\to f_k^{-}\)、\(\eta\)的计算，只要有一个无法计算，那贝叶斯滤波基本就废了。

那这个问题如何解决呢？

  1. 做假设

    1. \(f_{k-1}\),\(h(x)\)是线性的，\(Q_k, R_k\)都是正态分布

      \(\implies\)Kalman Filter (\(KF\)) \(F\)和\(h\)都是常数

      \[\begin{cases}
      X_k=F\cdot X_{k-1}+Q_k\\
      Y_k = h\cdot X_k+R_k
      \end{cases}
      \]

    2. \(f_{k-1}\),\(h(x)\)是非线性的，\(Q_k, R_k\)都是正态分布

      \(\implies\)\(UKF\)、\(EKF\)

  2. 霸王硬上弓，直接对无穷积分做数值积分
    1. 高斯积分。不太常用
    2. 蒙特卡罗积分\(\implies\)粒子滤波 (\(PF\))
    3. 直方图滤波

贝叶斯滤波更多的是一种思想，实际使用中，更多的是使用卡尔曼滤波等衍生方法。