【超详细的贝叶斯滤波原理】（不看后悔）

贝叶斯公式

二维离散型随机变量的贝叶斯公式

对于二维离散型随机变量$(X,Y)$，由其条件概率质量函数与全概率公式，容易得到其贝叶斯公式:
f X ∣ Y ( x ∣ y ) = f X , Y ( x , y ) f Y ( y ) = f Y ∣ X ( y ∣ x ) f X ( x ) ∑ i = 1 ∞ f Y ∣ X ( y ∣ x i ) f X ( x i ) , ( x , y ) ∈ { x i , y j } , i , j = 1 , 2 , 3 , ⋯ f_{X \mid Y}(x \mid y)=\frac{f_{X, Y}(x, y)}{f_{Y}(y)}=\frac{f_{Y \mid X}(y \mid x) f_{X}(x)}{\sum_{i=1}^{\infty} f_{Y \mid X}\left(y \mid x_{i}\right) f_{X}\left(x_{i}\right)},(x, y) \in\left\{x_{i}, y_{j}\right\}, i, j=1,2,3, \cdots fX∣Y​(x∣y)=fY​(y)fX,Y​(x,y)​=∑i=1∞​fY∣X​(y∣xi​)fX​(xi​)fY∣X​(y∣x)fX​(x)​,(x,y)∈{xi​,yj​},i,j=1,2,3,⋯
二维离散型随机变量的贝叶斯公式可通过作图的方式轻松证得。

二维连续型随机变量的贝叶斯公式

结论

对于二维连续型随机变量 $(X,Y)$，由其条件概率密度函数与全概率公式，容易得到其贝叶斯公式：
$$f_{X\mid Y}(x\mid y)=\frac{f_{X,Y}(x,y)}{f_Y(y)}=\frac{f_{Y\mid X}(y\mid x)f_X(x)}{{\sum }_{i=1}^{\infty }{f}_{Y\mid X}\left(y\mid x_i\right)f_X\left(x_i\right)}$$

推导

二维连续型随机变量的贝叶斯公式无法通过作图的方式推得，下面进行公式推导，首先计算二维连续型随机变量的条件累积分布函数：
$$\begin{array}{rl}{F}_{X\mid Y}(x\mid y)& =P(X\le x\mid Y=y)\\ & =\sum _{u=-\infty }^{x}P(X=u\mid Y=y)\Rightarrow \text{ 化连续为离散无穷小的累加 }\\ & =\sum _{u=-\infty }^x\frac{P(Y=y\mid X=u)P(X=u)}{P(Y=y)}\Rightarrow \text{ 二维离散型随机变量的贝 }\\ & =\underset{ϵ\to 0}{\lim }\sum _{u=-\infty }^x\frac{P(y\le Y\le y+ϵ\mid X=u)P(u\le X\le u+ϵ)}{P(y\le Y\le y+ϵ)}\Rightarrow \text{ 化无穷小为极限形式 }\\ & =\underset{ϵ\to 0}{\lim }\sum _{u=-\infty }^x\frac{\left[F_{Y\mid X}(y+ϵ\mid u)-F_{Y\mid X}(y\mid u)\right]\left[F_X(u+ϵ)-F_X(u)\right]}{F_Y(y+ϵ)-F_Y(y)}\Rightarrow \text{分布函数性质}\\ & =\underset{ϵ\to 0}{\lim }\sum _{u=-\infty }^x\frac{\left[f_{Y\mid X}\left({\xi }_1\mid u\right)\cdot ϵ\right]\left[f_X\left({\xi }_2\right)\cdot ϵ\right]}{f_Y\left({\xi }_3\right)\cdot ϵ}\Rightarrow \text{ 拉格朗日中值定理, }{\xi }_1,{\xi }_3\in (y,y+ϵ),{\xi }_2\in (u,u+ϵ)\\ & =\underset{ϵ\to 0}{\lim }\sum _{u=-\infty }^x\frac{f_{Y\mid X}(y\mid u)f_X(u)}{f_Y(y)}\cdot ϵ\Rightarrow ϵ\to 0\text{ 时, }{\xi }_1\to y,{\xi }_2\to u,{\xi }_3\to y\\ & ={\int }_{-\infty }^x\frac{f_{Y\mid X}(y\mid u)f_X(u)}{f_Y(y)}\mathrm{d}u\Rightarrow \text{ 积分定义 }\\ & ={\int }_{-\infty }^x\frac{f_{Y\mid X}(y\mid x)f_X(x)}{f_Y(y)}\mathrm{d}x\Rightarrow \text{ 替换自变量符号 }u\text{ 为 }x\end{array}$$
故，二维连续型随机变量的条件概率密度函数为：
$$f_{X\mid Y}(x\mid y)=\frac{\mathrm{d}{F}_{X\mid Y}(x\mid y)}{\mathrm{d}x}=\frac{f_{Y\mid X}(y\mid x)f_X(x)}{f_Y(y)}$$
代入全概率公式：
$$f_{X\mid Y}(x\mid y)=\frac{f_{Y\mid X}(y\mid x)f_X(x)}{{\int }_{-\infty }^{+\infty }{f}_{Y\mid X}(y\mid x)f_X(x)\mathrm{d}x}$$
上式即为二维连续型随机变量的贝叶斯公式。

先验概率、似然概率与后验概率

在二维连续型随机变量的贝叶斯公式中，有如下定义：

• $f_X(x)$ 被称为先验概率密度（Prior Probability Density），表示根据以往的经验和分析，在本次试验或采样前便可获得的随机变量 $X$ 的概率密度；
• $f_{Y|X}(y|x)$ 被称为似然概率密度（Likelihood Probability Density），表示在状态随机变量 $X$ 取值为 $x$ 的条件下，观测随机变量 $Y$ 取值为 $y$ 的概率密度，状态为因，观测为果，即由因推果；
• $f_{X|Y}(x|y)$ 被成为后验概率密度（Posterior Probability Density），表示在观测随机变量 $Y$ 取值为 $y$ 的条件下，状态随机变量 $X$ 取值为 $x$ 的概率密度，状态为因，观测为果，即由果推因。
  此外，当 $y$ 为定值时，$\eta ={\left[{\int }_{-\infty }^{+\infty }{f}_{Y\mid X}(y\mid x)f_X(x)\mathrm{d}x\right]}^{-1}$ 为一常数，常被称为贝叶斯公式的归一化常数。
  因此，二维连续型随机变量的贝叶斯公式可表示为：
  $$后验概率密度=\eta \times 似然概率密度\times 先验概率密度$$

再谈似然概率

上文中提到，似然概率密度函数$f_{Y|X}(y|x)$表示在状态随机变量 $X$ 取值为 $x$ 的条件下，观测随机变量 $Y$ 取值为 $y$ 的概率密度。似然概率密度函数表征了传感器检测精度，对于给定的状态条件 $X=x$，观测结果 $Y=y$ 的概率分布通常有三种模型：

1. 等可能型
   观测值在状态量真值附近呈均匀分布，此时的似然概率密度函数为常数。
2. 阶梯型
   观测值在状态量真值附近呈阶梯分布，此时的似然概率密度函数为分段常数。
3. 正态分布型
   观测值在状态量真值附近呈高斯分布，此时的似然概率密度函数为高斯函数：
   $$f_{Y\mid X}(y\mid x)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(y-x{)}^2}{2{\sigma }^2}}$$
   若假定似然概率密度函数为高斯函数，此时，似然概率密度函数的均值 $x$ 代表状态量真值，$\sigma$ 代表传感器检测精度范围。若同时假定先验概率密度函数为高斯函数，即：
   $$f_X(x)\sim \mathcal{N}\left({\mu }_1,{\sigma }_1^2\right),f_{Y\mid X}(y\mid x)\sim \mathcal{N}\left({\mu }_2,{\sigma }_2^2\right)$$
   则
   $$f_{X\mid Y}(x\mid y)\sim \mathcal{N}\left(\frac{{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}{\mu }_1+\frac{{\sigma }_1^2}{{\sigma }_1^2+{\sigma }_2^2}{\mu }_2,\frac{{\sigma }_1^2{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}\right)$$
   由于
   $$\frac{{\sigma }_1^2{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}=\frac{{\sigma }_1^2}{1+\frac{{\sigma }_1^2}{{\sigma }_2^2}}<{\sigma }_1^2$$
   且
   $$\frac{{\sigma }_1^2{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}=\frac{{\sigma }_2^2}{1+\frac{{\sigma }_2^2}{{\sigma }_1^2}}<{\sigma }_2^2$$
   故，后验概率密度函数方差既小于先验概率密度函数方差，也小于似然概率密度函数方差，系统不确定度降低。
   若 ${\sigma }_1^2\gg {\sigma }_2^2$ ，则近似有：
   $$f_{X\mid Y}(x\mid y)\sim \mathcal{N}\left({\mu }_2,{\sigma }_2^2\right)$$
   此时，后验倾向于观测。
   若 ${\sigma }_1^2\ll {\sigma }_2^2$，则近似有：
   $$f_{X\mid Y}(x\mid y)\sim \mathcal{N}\left({\mu }_1,{\sigma }_1^2\right)$$
   此时，后验倾向于先验。

贝叶斯滤波推导

问题建模

1. 问题描述
   对于某状态量随机变量 $X$，从初始时刻 $0$ 开始，对其进行观测，得到 $0-k$ 时刻的观测值：
   $$y_0,y_1,y_2,\cdots ,y_k$$
   求解 $k$ 时刻状态量随机变量 $X_k$ 的最优估计 ${\hat{x}}_k$。
2. 求解思路
   以贝叶斯公式为求解方向，将问题转化为求解状态量随机变量 $X_k$ 后验概率密度函数的期望：
   $${\hat{x}}_k=E\left[f_{X_k}^{+}(x)\right]$$
   进而需要求解状态量随机变量 $X_k$ 的先验概率密度函数与似然概率密度函数。我们认为，$k$ 时刻的状态量随机变量 $X_k$ 与且仅与上一时刻的状态量随机变量 $X_{k-1}$ 有关，$k$ 时刻的观测量随机变量 $Y_k$ 与且仅与 $k$ 时刻的状态量随机变量 $X_k$ 有关，其中的数量关系我们分别称之为状态方程与观测方程：
   $$\{\begin{array}{l}{X}_k=f\left(X_{k-1}\right)+Q_k\Rightarrow \text{ 状态方程 }\\ Y_k=h\left(X_k\right)+R_k\Rightarrow \text{ 观测方程 }\end{array}$$
   $f(x)$ 被称为状态转移函数，$h(x)$ 被称为观测函数。
   对于 0 时刻的初始状态量随机变量 $X_0$，认为观测值 $y_0$ 即为其真值，其后验概率密度函数即为其先验概率密度函数。我们可以根据经验知识（建模精度和传感器精度）写出 0 时刻的初始状态量随机变量 $X_0$ 的后验概率密度函数 $f_{X_0}^{+}(x)$、$k$ 时刻过程噪声随机变量 $Q_k$ 的概率密度函数 $f_{Q_k}(x)$ 和 $k$ 时刻观测噪声随机变量 $R_k$ 的概率密度函数 $f_{R_k}(x)$。
3. 符号定义

• 各时刻的状态量随机变量
  $$X_0,X_1,X_2,\cdots ,X_k$$
• 各时刻的观测量随机变量
  $$Y_0,Y_1,Y_2,\cdots ,Y_k$$
• 各时刻的观测值
  $$y_0,y_1,y_2,\cdots ,y_k$$
• 各时刻的过程噪声随机变量
  $$Q_1,Q_2,\cdots ,Q_k$$
• 各时刻的观测噪声随机变量
  $$R_1,R_2,\cdots ,R_k$$
• 各时刻的过程噪声随机变量概率密度函数
  $$f_{Q_1}(x),f_{Q_2}(x),\cdots ,f_{Q_k}(x)$$
• 各时刻的观测噪声随机变量概率密度函数
  $$f_{R_1}(x),f_{R_2}(x),\cdots ,f_{R_k}(x)$$
• 各时刻的状态量随机变量先验概率密度函数
  $$f_{X_0}^{-}(x),f_{X_1}^{-}(x),f_{X_2}^{-}(x),\cdots ,f_{X_k}^{-}(x)$$
• 各时刻的状态量随机变量后验概率密度函数
  $$f_{X_0}^{-}(x),f_{X_1}^{+}(x),f_{X_2}^{+}(x),\cdots ,f_{X_k}^{+}(x)$$
• 各时刻状态量随机变量与观测量随机变量的似然概率密度函数
  $$f_{X_0}^{-}(x),f_{X_1}^{+}(x),f_{X_2}^{+}(x),\cdots ,f_{X_k}^{+}(x)$$

4. 重要假设

• $X_0$ 分别与 $Q_1,Q_2,\cdots ,Q_k$ 相互独立；
• $X_0$ 分别与 $R_1,R_2,\cdots ,R_k$ 相互独立；
• $X_{k-1}$ 与 $Q_k$ 相互独立；
• $X_k$ 与 $R_k$ 相互独立。

5. 重要定理
   条件概率里的条件可以作逻辑推导。例如：
   $$P(X=1\mid Y=2,Z=3)=P(X+Y=3\mid Y=2,Z=3)=P(X+Y=3\mid Y=2,Z-Y=1)$$

预测步推导

已知 0 时刻状态量随机变量 $X_0$ 的后验概率密度函数 $f_{X_0}^{+}(x)$，状态转移函数 $f(x)$，1 时刻过程噪声随机变量 $Q_1$ 的概率密度函数 $f_{Q_1}(x)$，求解 1 时刻状态量随机变量 X1 的先验概率密度函数 $f_{X_1}^{-}(x)$。
类似二维连续型随机变量贝叶斯公式的推导过程，我们从求解 X1 的先验累积分布函数 $F_{X_1}^{-}$ 入手。

$$=\sum _{u=-\infty }^x\sum _{v=-\infty }^{+\infty }P\left(X_1=u\mid X_0=v\right)P\left(X_0=v\right)$$

---

故，1 时刻状态量随机变量 $X_1$ 的先验概率密度函数为：
$$f_{X_1}^{-}(x)=\frac{\mathrm{d}{F}_{X_1}^{-}(x)}{\mathrm{d}x}={\int }_{-\infty }^{+\infty }{f}_{Q_1}[x-f(v)]f_{X_0}^{-}(v)\mathrm{d}v$$
推导完毕。可以发现，先验概率密度函数本质来源于状态方程。

更新步推导

已知 1 时刻观测量随机变量 $Y_1$ 的取值 $y_1$，求解 1 时刻状态量随机变量与观测量随机变量的似然概率密度函数 $f_{Y_1|X_1}(y_1|x)$，并联合预测步得到的 1 时刻状态量随机变量 $X_1$ 的先验概率密度函数 $f_{X_1}^{-}(x)$，求解 1 时刻状态量随机变量 $X_1$ 的后验概率密度函数 $f_{X_1}^{+}(x)$。
首先，求解似然概率密度函数 $f_{Y_1|X_1}(y_1|x)$：
KaTeX parse error: Expected 'EOF', got '&' at position 41: … \mid x\right) &̲=\lim _{\epsilo…

---

可以发现，似然概率密度函数本质来源于观测方程。
然后，联合预测步得到的 1 时刻状态量随机变量 $X_1$ 的先验概率密度函数 $f_{X_1}^{-}(x)$，求解 1 时刻状态量随机变量 $X_1$ 的后验概率密度函数$f_{X_1}^{+}(x)$：
$$f_{X_1}^{-}(x)={\eta }_1\cdot f_{Y_1\mid X_1}\left(y_1\mid x\right)\cdot f_{X_1}^{-}(x)={\eta }_1\cdot f_{R_1}\left[y_1-h(x)\right]\cdot f_{X_1}^{-}(x)$$
其中，归一化常数 ${\eta }_1$为：
$${\eta }_1={\left[{\int }_{-\infty }^{+\infty }{f}_{Y_1\mid X_1}\left(y_1\mid x\right)f_{X_1}^{-}(x)\mathrm{d}x\right]}^{-1}={\left\{{\int }_{-\infty }^{+\infty }{f}_{R_1}\left[y_1-h(x)\right]f_{X_1}^{-}(x)\mathrm{d}x\right\}}^{-1}$$

递推流程

由预测步和更新步的推导结果，可得到由 0 时刻状态量随机变量 $X_0$ 的后验概率密度函数 $f_{X_0}^{+}(x)$ 到 $k$ 时刻状态量随机变量 $X_k$ 的后验概率密度函数 $f_{X_k}^{+}(x)$ 的递推流程：
$$\begin{gathered} f_{X_0}^{-}(x)\stackrel{预测}{⟹}{f}_{X_1}^{-}(x)={\int }_{-\infty }^{-\infty }{f}_{Q_1}[x-f(v)]f_{X_0}^{+}(v)\mathrm{d}v\stackrel{观测更新}{⟹}{f}_{X_1}^{+}(x)={\eta }_1\cdot f_{R_1}\left[y_1-h(x)\right]\cdot f_{X_1}^{-}(x) \\ \stackrel{预测}{⟹}{f}_{X_2}^{-}(x)={\int }_{-\infty }^{-\infty }{f}_{Q_2}[x-f(v)]f_{X_1}^{+}(v)\mathrm{d}v\stackrel{观测更新}{⟹}{f}_{X_2}^{+}(x)={\eta }_2\cdot f_{R_2}\left[y_2-h(x)\right]\cdot f_{X_2}^{-}(x) \\ \cdots \\ \stackrel{预测}{⟹}{f}_{X_k}^{-}(x)={\int }_{-\infty }^{-\infty }{f}_{Q_k}[x-f(v)]f_{X_{k-1}}^{+}(v)\mathrm{d}v\stackrel{观测更新}{⟹}{f}_{X_k}^{+}(x)={\eta }_k\cdot f_{R_k}\left[y_k-h(x)\right]\cdot f_{X_k}^{-}(x) \end{gathered}$$

其中，归一化常数 ${\eta }_k$ 为：
$${\eta }_k={\left\{{\int }_{-\infty }^{+\infty }{f}_{R_k}\left[y_k-h(x)\right]f_{X_k}^{-}(x)\mathrm{d}x\right\}}^{-1}$$
最终，可得到 $k$ 时刻状态量随机变量 $X_k$ 的最优估计 ${\hat{x}}_k$：
$${\hat{x}}_k=E\left[f_{X_k}^{+}(x)\right]={\int }_{-\infty }^{-\infty }x{f}_{X_k}^{-}(x)\mathrm{d}x$$

完整算法框架

1. 设初值
   初始 0 时刻状态量随机变量 $X_0$ 的后验概率密度函数：
   $$f_{X_0}^{+}(x)$$
2. 预测步
   $k$ 时刻状态量随机变量 $X_k$ 的先验概率密度函数：
   $$f_{X_k}^{-}(x)={\int }_{-\infty }^{+\infty }{f}_{Q_{k}k}[x-f(v)]f_{X_{k-1}}^{+}(v)\mathrm{d}v$$
3. 更新步
   $k$ 时刻状态量随机变量 $X_k$ 的后验概率密度函数：
   $$f_{X_k}^{+}(x)={\eta }_k\cdot f_{R_k}\left[y_k-h(x)\right]\cdot f_{X_k}^{-}(x)$$
   归一化常数 ${\eta }_k$：
   $${\eta }_k={\left\{{\int }_{-\infty }^{+\infty }{f}_{R_k}\left[y_k-h(x)\right]f_{X_k}^{-}(x)\mathrm{d}x\right\}}^{-1}$$
4. 求解状态量后验估计
   $k$ 时刻状态量随机变量 $X_k$ 的后验估计：
   $${\hat{x}}_k^{-}=E\left[f_{X_k}^{+}(x)\right]={\int }_{-\infty }^{-\infty }x{f}_{X_k}^{-}(x)\mathrm{d}x$$

贝叶斯滤波的缺点及解决方法

缺点

从上文的推导及结论中可以发现，求解预测步中的先验概率密度函数 $f_{X_k}^{-}(x)$、更新步中的归一化常数 ${\eta }_k$、最终的最优估计 ${\hat{x}}^k$ 时均涉及到无穷积分，而大多数情况无法得到解析解，使得贝叶斯滤波算法的直接应用十分困难。

解决办法

为了解决贝叶斯滤波中的无穷积分问题，通常从两个角度出发:

1. 作理想假设

• 假设状态转移函数 $f(x)$ 和观测函数 $h(x)$ 均为线性函数，过程噪声随机变量 $Q_k$ 和 观测噪声随机变量 $R_k$ 均服从均值为 0 的正态分布——卡尔曼滤波（Kalman Filter）
• 假设状态转移函数 $f(x)$ 和（或）观测函数 $h(x)$ 为非线性函数，过程噪声随机变量 $Q_k$ 和 观测噪声随机变量 $R_k$ 均服从均值为 0 的正态分布——扩展卡尔曼滤波（Extended Kalman Filter）和无迹卡尔曼滤波（Unscented Kalman Filter）

2. 化连续为离散
   将无穷积分转化为数值积分，一般有以下方法：

• 高斯积分（不常用）
• 蒙特卡罗积分（粒子滤波，Particle Filter）
• 直方图滤波