概述
文章目录
- 贝叶斯滤波(*Bayesian filtering equations*)
- 贝叶斯滤波方程
- Kalman滤波
- 贝叶斯平滑(*Bayesian smoothing*)
- 贝叶斯最优平滑方程
- Rauch-Tung-Striebel smoother
贝叶斯滤波(Bayesian filtering equations)
贝叶斯滤波的目的是在观测到信号
y
1
:
k
boldsymbol y_{1:k}
y1:k后,得到当前第
k
k
k个时刻的状态
x
k
boldsymbol{ x}_k
xk的后验边际分布:
p
(
x
k
∣
y
1
:
k
)
p(boldsymbol x_k| boldsymbol y_{1:k})
p(xk∣y1:k)
贝叶斯滤波方程
定理(Bayesian filtering equations)
The recursive equations (the Bayesian filter) for computing the predicted distribution
p
(
x
k
∣
y
1
:
k
−
1
)
p(boldsymbol x_k| boldsymbol y_{1:k-1})
p(xk∣y1:k−1) and the filtering distribution
p
(
x
k
∣
y
1
:
k
)
p(boldsymbol x_k| boldsymbol y_{1:k})
p(xk∣y1:k) at the time
k
k
k are given by the following Bayesian filtering eqautions.
-
Initialization: 迭代从先验分布 p ( x 0 ) p(boldsymbol{x}_0) p(x0)开始
-
Prediction: 给定动态模型(dynamic model) p ( x k ∣ x k − 1 ) p(boldsymbol x_k| boldsymbol x_{k-1}) p(xk∣xk−1),状态 x k boldsymbol{x}_k xk的预测分布(predictive distribution)为:
p ( x k ∣ y 1 : k − 1 ) = ∫ p ( x k ∣ x k − 1 ) p ( x k − 1 ∣ y 1 : k − 1 ) d x k − 1 (1) p(boldsymbol x_{k} | boldsymbol y_{1:k-1} ) = int p(boldsymbol x_{k}|boldsymbol x_{k-1}) p(boldsymbol x_{k-1}| boldsymbol y_{1:k-1}) text{d} boldsymbol x_{k-1} tag{1} p(xk∣y1:k−1)=∫p(xk∣xk−1)p(xk−1∣y1:k−1)dxk−1(1) -
Update: 给定前 k k k次观测向量 y 1 : k boldsymbol{ y}_{1:k} y1:k,状态 x k boldsymbol{ x}_k xk的后验分布为
p ( x k ∣ y 1 : k ) = 1 Z k p ( y k ∣ x k ) p ( x k ∣ y 1 : k − 1 ) (2) p(boldsymbol x_{k} | boldsymbol y_{1:k} ) =frac{1}{Z_k} p(boldsymbol y_k| boldsymbol x_k) p(boldsymbol x_{k} | boldsymbol y_{1:k-1} )tag{2} p(xk∣y1:k)=Zk1p(yk∣xk)p(xk∣y1:k−1)(2)其中 Z k Z_k Zk是归一化系数:
Z k = ∫ 1 Z k p ( y k ∣ x k ) p ( x k ∣ y 1 : k − 1 ) d x k (3) Z_k = int frac{1}{Z_k} p(boldsymbol y_k| boldsymbol x_k) p(boldsymbol x_{k} | boldsymbol y_{1:k-1} ) text{d} boldsymbol x_ktag{3} Zk=∫Zk1p(yk∣xk)p(xk∣y1:k−1)dxk(3)
根据之前的博客概率状态空间模型的马尔可夫性质和观测向量的条件独立性,不难证明上述等式,这里省略。
Kalman滤波
在经典Kalman滤波中,动态模型(dynamic model)和测量模型(measurement model)都是线性高斯的:
dynamic model:
x
k
=
A
k
−
1
x
k
−
1
+
q
k
−
1
measurement model:
y
k
=
H
k
x
k
+
r
k
(4)
begin{aligned} text{dynamic model: }& boldsymbol x_k = boldsymbol A_{k-1} boldsymbol x_{k-1} + boldsymbol q_{k-1} \ text{measurement model: }& boldsymbol y_k = boldsymbol H_k boldsymbol x_k + boldsymbol r_k tag{4} end{aligned}
dynamic model: measurement model: xk=Ak−1xk−1+qk−1yk=Hkxk+rk(4)
其中
- x k ∈ R n boldsymbol x_k in mathbb{R}^{n} xk∈Rn 是 k k k时刻的状态
- y k ∈ R m boldsymbol y_k in mathbb{R}^{m} yk∈Rm 是 k k k时刻的测量向量
- q k − 1 ∼ N ( 0 , Q k − 1 ) boldsymbol q_{k-1} sim mathcal{N}(boldsymbol{0}, boldsymbol{Q}_{k-1}) qk−1∼N(0,Qk−1)是过程噪声(process noise)
- r k ∼ N ( 0 , R k ) boldsymbol r_k sim mathcal{N}(boldsymbol{0}, boldsymbol{R}_{k}) rk∼N(0,Rk)是测量噪声
- x 0 ∼ N ( m 0 , P 0 ) boldsymbol x_0 sim mathcal{N}(boldsymbol m_0, boldsymbol P_0) x0∼N(m0,P0)是状态的先验分布
- A k − 1 boldsymbol A_{k-1} Ak−1是动态模型中的变换矩阵(transition matrix)
- H k boldsymbol H_{k} Hk是测量矩阵
用概率的形式表示模型(4)为:
p
(
x
k
∣
x
k
−
1
)
=
N
(
x
k
∣
A
k
−
1
x
k
−
1
,
Q
k
−
1
)
p
(
y
k
∣
x
k
)
=
N
(
x
k
∣
H
k
x
k
,
R
k
)
(5)
begin{aligned} p(boldsymbol x_{k} | boldsymbol x_{k-1} ) & = mathcal{N} (boldsymbol x_k| boldsymbol A_{k-1} boldsymbol x_{k-1}, boldsymbol{Q}_{k-1}) \ p(boldsymbol y_{k} | boldsymbol x_{k} ) & = mathcal{N} (boldsymbol x_k| boldsymbol H_{k} boldsymbol x_{k}, boldsymbol{R}_{k}) end{aligned} tag{5}
p(xk∣xk−1)p(yk∣xk)=N(xk∣Ak−1xk−1,Qk−1)=N(xk∣Hkxk,Rk)(5)
定理(Kalman Filter)
The Bayesian filtering equations for the linear filtering model (4) can be evaluated in closed form and the resulting distributions are Gaussian:
p
(
x
k
∣
y
1
:
k
−
1
)
=
N
(
x
k
∣
m
k
−
,
P
k
−
)
p
(
x
k
∣
y
1
:
k
)
=
N
(
x
k
∣
m
k
,
P
k
)
p
(
y
k
∣
y
1
:
k
−
1
)
=
N
(
y
k
∣
H
k
m
k
−
,
S
k
)
(6)
begin{aligned} p(boldsymbol x_k| boldsymbol y_{1:k-1}) &= mathcal{N}(boldsymbol x_k| boldsymbol m^{-}_k, boldsymbol P^{-}_k) \ p(boldsymbol x_k| boldsymbol y_{1:k}) &= mathcal{N}(boldsymbol x_k| boldsymbol m^{}_k, boldsymbol P_k) \ p(boldsymbol y_k | boldsymbol y_{1:k-1})&= mathcal{N}(boldsymbol y_k| boldsymbol H_k boldsymbol m^{-}_k, boldsymbol S_k) end{aligned} tag{6}
p(xk∣y1:k−1)p(xk∣y1:k)p(yk∣y1:k−1)=N(xk∣mk−,Pk−)=N(xk∣mk,Pk)=N(yk∣Hkmk−,Sk)(6)
上述分布的参数的迭代可以被描述为
-
prediction step:
m k − = A k − 1 m k − 1 P k − = A k − 1 P k − 1 A k − 1 T + Q k − 1 (7) begin{aligned} boldsymbol m^{-}_k & = boldsymbol A_{k-1} boldsymbol m_{k-1} \ boldsymbol P^{-}_k &=boldsymbol A_{k-1} boldsymbol P_{k-1} boldsymbol A^T_{k-1} + boldsymbol Q_{k-1} end{aligned} tag{7} mk−Pk−=Ak−1mk−1=Ak−1Pk−1Ak−1T+Qk−1(7) -
update step
v k = y k − H k m k − S k = H k P k − 1 H k T + R k K k = P k − 1 H k T + R k m k = m k − + K k v k P k = P k − 1 − K k S k K k T (8) begin{aligned} boldsymbol v_k &= boldsymbol y_k - boldsymbol H_k boldsymbol m^{-}_k \ boldsymbol S_k &= boldsymbol H_k boldsymbol P^{-1}_k boldsymbol H^T_k + boldsymbol R_k \ boldsymbol{ K}_k &= boldsymbol P^{-1}_k boldsymbol H^T_k + boldsymbol R_k \ boldsymbol m_{k} &= boldsymbol m^{-}_k + boldsymbol{ K}_k boldsymbol v_k \ boldsymbol P_k &= boldsymbol P^{-1}_k - boldsymbol{ K}_k boldsymbol S_k boldsymbol{ K}^T_k tag{8} end{aligned} vkSkKkmkPk=yk−Hkmk−=HkPk−1HkT+Rk=Pk−1HkT+Rk=mk−+Kkvk=Pk−1−KkSkKkT(8)
上述迭代过程从先验分布的均值
m
0
boldsymbol{ m}_0
m0和协方差
P
0
boldsymbol{P}_0
P0开始。
证明:
(1)prediction相关
联合高斯分布与条件高斯分布的相关性质(之前写的博客)总结了联合高斯分布的表达式,据此,我们可以得到
p
(
x
k
−
1
,
x
k
∣
y
1
:
k
−
1
)
=
p
(
x
k
∣
x
k
−
1
)
p
(
x
k
−
1
∣
y
1
:
k
−
1
)
=
N
(
x
k
∣
A
k
−
1
x
k
−
1
,
Q
k
−
1
)
⋅
N
(
x
k
−
1
∣
m
k
−
1
,
P
k
−
1
)
=
N
(
[
x
k
−
1
x
k
]
∣
m
′
,
P
′
)
begin{aligned} & p(boldsymbol x_{k-1} ,boldsymbol x_k | boldsymbol y_{1:k-1}) \ & = p(boldsymbol x_{k} | boldsymbol x_{k-1} ) p(boldsymbol x_{k-1} | boldsymbol y_{1:k-1}) \ & = mathcal{N}(boldsymbol x_k| boldsymbol A_{k-1} boldsymbol x_{k-1}, boldsymbol Q_{k-1}) cdot mathcal{N}(boldsymbol x_{k-1}| boldsymbol m_{k-1}, boldsymbol P_{k-1}) \ & = mathcal{N} left( left[ begin{array}{c} boldsymbol{x}_{k-1}\ boldsymbol{x}_k\ end{array} right] |boldsymbol{m}^{prime},boldsymbol{P}^{prime} right) end{aligned}
p(xk−1,xk∣y1:k−1)=p(xk∣xk−1)p(xk−1∣y1:k−1)=N(xk∣Ak−1xk−1,Qk−1)⋅N(xk−1∣mk−1,Pk−1)=N([xk−1xk]∣m′,P′)
其中
m
′
=
[
m
k
−
1
A
k
−
1
m
k
−
1
]
P
′
=
[
P
k
−
1
P
k
−
1
A
k
−
1
T
A
k
−
1
P
k
−
1
A
k
−
1
P
k
−
1
A
k
−
1
T
+
Q
k
−
1
]
begin{aligned} boldsymbol{m}^{prime} &= left[ begin{array}{c} boldsymbol{m}_{k-1}\ boldsymbol{A}_{k-1}boldsymbol{m}_{k-1}\ end{array} right] \ boldsymbol P^{prime} & = left[ begin{matrix} boldsymbol{P}_{k-1}& boldsymbol{P}_{k-1}boldsymbol{A}_{k-1}^{T}\ boldsymbol{A}_{k-1}boldsymbol{P}_{k-1}& boldsymbol{A}_{k-1}boldsymbol{P}_{k-1}boldsymbol{A}_{k-1}^{T}+boldsymbol{Q}_{k-1}\ end{matrix} right] end{aligned}
m′P′=[mk−1Ak−1mk−1]=[Pk−1Ak−1Pk−1Pk−1Ak−1TAk−1Pk−1Ak−1T+Qk−1]
根据上述博客总结中的边际概率表达式,我们可以写出
p
(
x
k
∣
y
1
:
k
−
1
)
=
N
(
x
k
∣
m
k
−
,
P
k
−
)
p(boldsymbol x_k| boldsymbol y_{1:k-1}) = mathcal{N}(boldsymbol x_k| boldsymbol m^{-}_k, boldsymbol P^{-}_k)
p(xk∣y1:k−1)=N(xk∣mk−,Pk−)
其中
m
k
−
=
A
k
−
1
m
k
−
1
,
P
k
−
=
A
k
−
1
P
k
−
1
A
k
−
1
T
+
Q
k
−
1
boldsymbol m^{-}_k = boldsymbol{A}_{k-1}boldsymbol{m}_{k-1}, boldsymbol{P}^{-}_{k} = boldsymbol{A}_{k-1}boldsymbol{P}_{k-1}boldsymbol{A}_{k-1}^{T}+boldsymbol{Q}_{k-1}
mk−=Ak−1mk−1, Pk−=Ak−1Pk−1Ak−1T+Qk−1
(2)Update相关
类似地,我们可以得到
p
(
x
k
,
y
k
∣
y
1
:
k
−
1
)
=
p
(
y
k
∣
x
k
)
p
(
x
k
∣
y
1
:
k
−
1
)
=
N
(
y
k
∣
H
k
x
k
,
R
k
)
⋅
N
(
x
k
∣
m
k
−
,
P
k
−
)
=
N
(
[
x
k
y
k
]
∣
m
′
′
,
P
′
′
)
begin{aligned} & p(boldsymbol x_{k} ,boldsymbol y_k | boldsymbol y_{1:k-1}) \ & = p(boldsymbol y_{k} | boldsymbol x_{k} ) p(boldsymbol x_{k} | boldsymbol y_{1:k-1}) \ & = mathcal{N}(boldsymbol y_k| boldsymbol H_{k} boldsymbol x_{k}, boldsymbol R_{k}) cdot mathcal{N}(boldsymbol x_{k}| boldsymbol m^{-}_{k}, boldsymbol P^{-}_{k}) \ & = mathcal{N} left( left[ begin{array}{c} boldsymbol{x}_{k}\ boldsymbol{y}_k\ end{array} right] |boldsymbol{m}^{ prime prime},boldsymbol{P}^{prime prime} right) end{aligned}
p(xk,yk∣y1:k−1)=p(yk∣xk)p(xk∣y1:k−1)=N(yk∣Hkxk,Rk)⋅N(xk∣mk−,Pk−)=N([xkyk]∣m′′,P′′)
其中
m
′
′
=
[
m
k
−
H
k
m
k
−
]
P
′
′
=
[
P
k
−
P
k
−
H
k
T
H
k
P
k
−
H
k
P
k
−
H
k
T
+
R
k
]
begin{aligned} boldsymbol{m}^{prime prime} &= left[ begin{array}{c} boldsymbol{m}^{-}_{k}\ boldsymbol{H}_{k}boldsymbol{m}^{-}_{k}\ end{array} right] \ boldsymbol P^{prime prime} & = left[ begin{matrix} boldsymbol{P}^{-}_{k}& boldsymbol{P}^{-}_{k}boldsymbol{H}_{k}^{T}\ boldsymbol{H}_{k}boldsymbol{P}^{-}_{k}& boldsymbol{H}_{k}boldsymbol{P}^{-}_{k} boldsymbol{H}_{k}^{T} +boldsymbol{R}_{k}\ end{matrix} right] end{aligned}
m′′P′′=[mk−Hkmk−]=[Pk−HkPk−Pk−HkTHkPk−HkT+Rk]
进一步,关于
p
(
x
k
,
y
k
∣
y
1
:
k
−
1
)
p(boldsymbol x_{k} ,boldsymbol y_k | boldsymbol y_{1:k-1})
p(xk,yk∣y1:k−1)的边际概率分布为
p
(
x
k
∣
y
k
,
y
1
:
k
−
1
)
=
p
(
x
k
∣
y
1
:
k
)
=
N
(
x
k
∣
m
k
,
P
k
)
begin{aligned} p(boldsymbol x_{k} |boldsymbol y_k , boldsymbol y_{1:k-1}) & = p(boldsymbol x_{k} |boldsymbol y_{1:k}) \ & = mathcal{N} (boldsymbol x_k| boldsymbol m_k ,boldsymbol P_k) end{aligned}
p(xk∣yk,y1:k−1)=p(xk∣y1:k)=N(xk∣mk,Pk)
其中
m
k
=
m
k
−
+
P
k
−
H
k
T
(
H
k
P
k
−
H
k
T
+
R
k
)
−
1
(
y
k
−
H
k
m
k
−
1
)
P
k
=
P
k
−
−
P
k
−
H
k
T
(
H
k
P
k
−
H
k
T
+
R
k
)
−
1
H
k
P
k
−
begin{aligned} boldsymbol m_k &= boldsymbol m^{-}_k + boldsymbol P^{-}_k boldsymbol H^T_k left ( boldsymbol H_k boldsymbol P^{-}_k boldsymbol H^T_k + boldsymbol R_k right)^{-1} left ( boldsymbol y_k - boldsymbol H_k boldsymbol m^{-1}_k right) \ boldsymbol P_k & = boldsymbol P^{-}_k - boldsymbol P^{-}_k boldsymbol H^T_k left ( boldsymbol H_k boldsymbol P^{-}_k boldsymbol H^T_k + boldsymbol R_k right)^{-1} boldsymbol H_k boldsymbol P^{-}_k end{aligned}
mkPk=mk−+Pk−HkT(HkPk−HkT+Rk)−1(yk−Hkmk−1)=Pk−−Pk−HkT(HkPk−HkT+Rk)−1HkPk−
上述迭代过程,也可以被写为式(7,8)的形式。完成证明。
贝叶斯平滑(Bayesian smoothing)
贝叶斯平滑的目的是在观测到信号
y
1
:
T
boldsymbol y_{1:T}
y1:T后,计算之前第
k
k
k个时刻的状态
x
k
boldsymbol{ x}_k
xk的后验边际概率:
p
(
x
k
∣
y
1
:
T
)
,
T
>
k
p(boldsymbol x_k| boldsymbol y_{1:T}), T >k
p(xk∣y1:T), T>k
贝叶斯最优平滑方程
定理(Bayesian optimal smoothing equations)
The backward recursive equations (the Bayesian smoother) for computing the smoothed distributions
p
(
x
k
∣
y
1
:
T
)
p(boldsymbol x_k| boldsymbol y_{1:T})
p(xk∣y1:T) for any
k
<
T
k < T
k<T are given by the following Bayesian (fixed-interval) smooth equations:
p
(
x
k
+
1
∣
y
1
:
k
)
=
∫
p
(
x
k
+
1
∣
x
k
)
p
(
x
k
∣
y
1
:
k
)
d
x
k
p
(
x
k
∣
y
1
:
T
)
=
p
(
x
k
∣
y
1
:
k
)
∫
[
p
(
x
k
+
1
∣
x
K
)
p
(
x
k
+
1
∣
y
1
:
T
)
p
(
x
k
+
1
∣
y
1
:
k
)
]
d
x
k
+
1
(9)
begin{aligned} p(boldsymbol x_{k+1} | boldsymbol y_{1:k} ) &= int p(boldsymbol x_{k+1}|boldsymbol x_k) p(boldsymbol x_k| boldsymbol y_{1:k}) text{d} boldsymbol x_{k} \ p(boldsymbol x_k| boldsymbol y_{1:T}) &= p(boldsymbol x_k | boldsymbol y_{1:k}) int left [ frac{p(boldsymbol x_{k+1}| boldsymbol x_K) p(boldsymbol x_{k+1}| boldsymbol y_{1:T})}{p(boldsymbol x_{k+1}| boldsymbol y_{1:k})} right] text{d} boldsymbol x_{k+1} tag{9} end{aligned}
p(xk+1∣y1:k)p(xk∣y1:T)=∫p(xk+1∣xk)p(xk∣y1:k)dxk=p(xk∣y1:k)∫[p(xk+1∣y1:k)p(xk+1∣xK)p(xk+1∣y1:T)]dxk+1(9)
证明:注意式(1)的第一个等式,积分项里边本应为: p ( x k + 1 ∣ x k , y 1 : k ) p ( x k ∣ y 1 : k ) p(boldsymbol x_{k+1}|boldsymbol x_k,boldsymbol y_{1:k}) p(boldsymbol x_k| boldsymbol y_{1:k}) p(xk+1∣xk,y1:k)p(xk∣y1:k),但是因为马尔可夫性,所以 p ( x k + 1 ∣ x k , y 1 : k ) = p ( x k + 1 ∣ x k ) p(boldsymbol x_{k+1}|boldsymbol x_k,boldsymbol y_{1:k}) = p(boldsymbol x_{k+1}|boldsymbol x_k) p(xk+1∣xk,y1:k)=p(xk+1∣xk)。
回顾上一篇博客概率状态空间模型的马尔可夫性质,我们知道
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
k
)
p(boldsymbol x_k|boldsymbol x_{k+1},boldsymbol y_{1:T}) = p(boldsymbol x_k| boldsymbol x_{k+1},boldsymbol y_{1:k})
p(xk∣xk+1,y1:T)=p(xk∣xk+1,y1:k),借助Bayes公式,我们有
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
k
)
=
p
(
x
k
,
x
k
+
1
∣
y
1
:
k
)
p
(
x
k
+
1
∣
y
1
:
k
)
=
p
(
x
k
+
1
∣
x
k
,
y
1
:
k
)
p
(
x
k
∣
y
1
:
k
)
p
(
x
k
+
1
∣
y
1
:
k
)
=
p
(
x
k
+
1
∣
x
k
)
p
(
x
k
∣
y
1
:
k
)
p
(
x
k
+
1
∣
y
1
:
k
)
begin{aligned} p(boldsymbol x_k|boldsymbol x_{k+1},boldsymbol y_{1:T}) &= p(boldsymbol x_k| boldsymbol x_{k+1},boldsymbol y_{1:k}) \ & = frac{ p(boldsymbol x_k ,boldsymbol x_{k+1}|boldsymbol y_{1:k}) } {p(boldsymbol x_{k+1}|boldsymbol y_{1:k})} \ & = frac{p(boldsymbol x_{k+1} |boldsymbol x_{k},boldsymbol y_{1:k}) p(boldsymbol x_{k}|boldsymbol y_{1:k})} {p(boldsymbol x_{k+1}|boldsymbol y_{1:k})} \ & = frac{p(boldsymbol x_{k+1} |boldsymbol x_{k}) p(boldsymbol x_{k}|boldsymbol y_{1:k})} {p(boldsymbol x_{k+1}|boldsymbol y_{1:k})} end{aligned}
p(xk∣xk+1,y1:T)=p(xk∣xk+1,y1:k)=p(xk+1∣y1:k)p(xk,xk+1∣y1:k)=p(xk+1∣y1:k)p(xk+1∣xk,y1:k)p(xk∣y1:k)=p(xk+1∣y1:k)p(xk+1∣xk)p(xk∣y1:k)
给定
y
1
:
T
boldsymbol y_{1:T}
y1:T,
x
k
boldsymbol x_k
xk和
x
k
+
1
boldsymbol x_{k+1}
xk+1的联合分布为:
p
(
x
k
,
x
k
+
1
∣
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
p
(
x
k
+
1
∣
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
k
)
p
(
x
k
+
1
∣
y
1
:
T
)
=
p
(
x
k
+
1
∣
x
k
)
p
(
x
k
∣
y
1
:
k
)
p
(
x
k
+
1
∣
y
1
:
T
)
p
(
x
k
+
1
∣
y
1
:
k
)
begin{aligned} p(boldsymbol x_k,boldsymbol x_{k+1}|boldsymbol y_{1:T}) &= p(boldsymbol x_k|boldsymbol x_{k+1},boldsymbol y_{1:T}) p(boldsymbol x_{k+1}|boldsymbol y_{1:T}) \ &= p(boldsymbol x_k|boldsymbol x_{k+1},boldsymbol y_{1:k}) p(boldsymbol x_{k+1}|boldsymbol y_{1:T}) \ & = frac{p(boldsymbol x_{k+1} |boldsymbol x_{k}) p(boldsymbol x_{k}|boldsymbol y_{1:k}) p(boldsymbol x_{k+1}|boldsymbol y_{1:T})} {p(boldsymbol x_{k+1}|boldsymbol y_{1:k})} end{aligned}
p(xk,xk+1∣y1:T)=p(xk∣xk+1,y1:T)p(xk+1∣y1:T)=p(xk∣xk+1,y1:k)p(xk+1∣y1:T)=p(xk+1∣y1:k)p(xk+1∣xk)p(xk∣y1:k)p(xk+1∣y1:T)
其中 p ( x k + 1 ∣ y 1 : T ) p(boldsymbol x_{k+1}|boldsymbol y_{1:T}) p(xk+1∣y1:T)是时刻 k + 1 k+1 k+1的平滑分布(smoothed distribution),关于 x k boldsymbol{x}_k xk的边际分布就是对上式关于 x k + 1 boldsymbol{x}_{k+1} xk+1积分,得到式(1)。
Rauch-Tung-Striebel smoother
Rauch-Tung-Striebel smoother(RTSS)的平滑解:
p
(
x
k
∣
y
1
:
T
)
=
N
(
x
k
∣
m
k
s
,
P
k
s
)
(10)
p(boldsymbol x_k | boldsymbol y_{1:T}) = mathcal{N}(boldsymbol x_k | boldsymbol m^s_k, boldsymbol P^s_k) tag{10}
p(xk∣y1:T)=N(xk∣mks,Pks)(10)
定理(RTS smoother)
The backward recursion equations for the (fixed interval) Rauch-Tung-Striebel smoother are given as
m
k
+
1
−
=
A
k
m
k
P
k
+
1
−
=
A
k
P
k
A
k
T
+
Q
k
G
k
=
P
k
A
k
T
[
P
k
+
1
−
]
−
1
m
k
s
=
m
k
+
G
k
[
m
k
+
1
s
−
m
k
+
1
−
]
P
k
s
=
P
k
+
G
k
[
P
k
+
1
s
−
P
k
+
1
−
]
G
k
T
(11)
begin{aligned} boldsymbol m^{-}_{k+1} & = boldsymbol A_{k} boldsymbol m_{k} \ boldsymbol P^{-}_{k+1} &=boldsymbol A_{k} boldsymbol P_{k} boldsymbol A^T_{k} + boldsymbol Q_{k} \ boldsymbol G_k &= boldsymbol P_{k} boldsymbol A^T_{k}[boldsymbol P^{-}_{k+1}]^{-1} \ boldsymbol m^s_{k} &= boldsymbol m^{}_k + boldsymbol{ G}_k [boldsymbol m^s_{k+1} - boldsymbol m^{-}_{k+1} ]\ boldsymbol P^{s}_{k} &= boldsymbol P^{}_k + boldsymbol{ G}_k [boldsymbol P^{s}_{k+1} - boldsymbol P^{-}_{k+1}]boldsymbol G^T_k tag{11} end{aligned}
mk+1−Pk+1−GkmksPks=Akmk=AkPkAkT+Qk=PkAkT[Pk+1−]−1=mk+Gk[mk+1s−mk+1−]=Pk+Gk[Pk+1s−Pk+1−]GkT(11)
其中 m k boldsymbol{m}_k mk和 P k boldsymbol P_k Pk由Kalman滤波得到。迭代过程从时刻 T T T开始,且 m T s = m T boldsymbol{m}^s_T = boldsymbol{m}_T mTs=mT, P T s = P T boldsymbol P^s_T = boldsymbol P_T PTs=PT
联合高斯分布与条件高斯分布的相关性质(之前写的博客)总结了联合高斯分布的表达式,据此,给定
y
1
:
k
boldsymbol{y}_{1:k}
y1:k,
x
k
boldsymbol{x}_k
xk和
x
k
+
1
boldsymbol{x}_{k+1}
xk+1的联合概率密度函数为
p
(
x
k
,
x
k
+
1
∣
y
1
:
k
)
=
p
(
x
k
+
1
∣
x
k
)
p
(
x
k
∣
y
1
:
k
)
=
N
(
x
k
+
1
∣
A
k
x
k
,
Q
k
)
⋅
N
(
x
k
∣
m
k
,
P
k
)
=
N
(
[
x
k
x
k
+
1
]
∣
m
~
1
,
P
~
1
)
begin{aligned} & p(boldsymbol x_{k} ,boldsymbol x_{k+1} | boldsymbol y_{1:k}) \ & = p(boldsymbol x_{k+1} | boldsymbol x_{k} ) p(boldsymbol x_{k} | boldsymbol y_{1:k}) \ & = mathcal{N}(boldsymbol x_{k+1}| boldsymbol A_{k} boldsymbol x_{k}, boldsymbol Q_{k}) cdot mathcal{N}(boldsymbol x_{k}| boldsymbol m_{k}, boldsymbol P_{k}) \ & = mathcal{N} left( left[ begin{array}{c} boldsymbol{x}_{k}\ boldsymbol{x}_{k+1}\ end{array} right] | tilde{ boldsymbol{m}}_{1},tilde{ boldsymbol{P}}_{1}right) end{aligned}
p(xk,xk+1∣y1:k)=p(xk+1∣xk)p(xk∣y1:k)=N(xk+1∣Akxk,Qk)⋅N(xk∣mk,Pk)=N([xkxk+1]∣m~1,P~1)
其中
m
~
1
=
[
m
k
A
k
m
k
]
P
~
1
=
[
P
k
P
k
A
k
T
A
k
P
k
A
k
P
k
A
k
T
+
Q
k
]
begin{aligned} tilde{ boldsymbol{m}}_{1} &= left[ begin{array}{c} boldsymbol{m}_{k}\ boldsymbol{A}_{k}boldsymbol{m}_{k}\ end{array} right] \ tilde{ boldsymbol{P}}_{1} & = left[ begin{matrix} boldsymbol{P}_{k}& boldsymbol{P}_{k}boldsymbol{A}_{k}^{T}\ boldsymbol{A}_{k}boldsymbol{P}_{k}& boldsymbol{A}_{k}boldsymbol{P}_{k}boldsymbol{A}_{k}^{T}+boldsymbol{Q}_{k}\ end{matrix} right] end{aligned}
m~1P~1=[mkAkmk]=[PkAkPkPkAkTAkPkAkT+Qk]
根据状态的马尔可夫性,我们有
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
k
)
p(boldsymbol x_{k} | boldsymbol x_{k+1} , boldsymbol y_{1:T}) = p(boldsymbol x_{k} | boldsymbol x_{k+1} , boldsymbol y_{1:k})
p(xk∣xk+1,y1:T)=p(xk∣xk+1,y1:k)
根据上述博客总结中的边际概率表达式,我们可以写出
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
k
)
=
N
(
x
k
∣
m
~
2
,
P
~
2
)
begin{aligned} p(boldsymbol x_{k} | boldsymbol x_{k+1} , boldsymbol y_{1:T}) &= p(boldsymbol x_{k} | boldsymbol x_{k+1} , boldsymbol y_{1:k}) \ & = mathcal{N} (boldsymbol x_k| tilde {boldsymbol m}_2, tilde {boldsymbol P}_2) end{aligned}
p(xk∣xk+1,y1:T)=p(xk∣xk+1,y1:k)=N(xk∣m~2,P~2)
其中
G
k
=
P
k
A
k
T
(
A
k
P
k
A
k
T
+
Q
k
)
−
1
m
~
2
=
m
k
+
G
k
(
x
k
+
1
−
A
k
m
k
)
P
~
2
=
P
k
−
G
k
(
A
k
P
k
A
k
T
+
Q
k
)
G
k
T
begin{aligned} boldsymbol G_k & = boldsymbol P_k boldsymbol A^T_k left ( boldsymbol A_k boldsymbol P_k boldsymbol A^T_k + boldsymbol Q_k right)^{-1} \ tilde { boldsymbol m}_2 & = boldsymbol m_k + boldsymbol G_k left ( boldsymbol x_{k+1} - boldsymbol A_k boldsymbol m_k right) \ tilde { boldsymbol P }_2 &= boldsymbol P_k - boldsymbol G_k left ( boldsymbol A_k boldsymbol P_k boldsymbol A^T_k + boldsymbol Q_k right) boldsymbol G^T_k end{aligned}
Gkm~2P~2=PkAkT(AkPkAkT+Qk)−1=mk+Gk(xk+1−Akmk)=Pk−Gk(AkPkAkT+Qk)GkT
给定所有数据
y
1
:
T
boldsymbol{y}_{1:T}
y1:T,
x
k
boldsymbol{x}_k
xk和
x
k
+
1
boldsymbol{x}_{k+1}
xk+1的联合分布为
p
(
x
k
+
1
,
x
k
∣
y
1
:
T
)
=
p
(
x
k
∣
x
k
+
1
,
y
1
:
T
)
p
(
x
k
+
1
∣
y
1
:
T
)
=
N
(
x
k
∣
m
~
2
,
P
~
2
)
N
(
x
k
+
1
∣
m
k
+
1
s
,
P
k
+
1
s
)
=
N
(
[
x
k
+
1
x
k
]
∣
m
~
3
,
P
~
3
)
begin{aligned} p(boldsymbol{x}_{k+1}, boldsymbol{x_k}|boldsymbol{y}_{1:T}) &= p( boldsymbol{x}_k | boldsymbol{x}_{k+1},boldsymbol{y}_{1:T}) p (boldsymbol{x}_{k+1} | boldsymbol{y}_{1:T}) \ &= mathcal{N}(boldsymbol{x}_k| tilde{boldsymbol{m}}_2, tilde{ boldsymbol{P}}_2) mathcal{N}(boldsymbol{x}_{k+1}| boldsymbol{m}^s_{k+1},boldsymbol{P}^s_{k+1}) \ & = mathcal{N} left( left[ begin{array}{c} boldsymbol{x}_{k+1}\ boldsymbol{x}_{k}\ end{array} right] | tilde{ boldsymbol{m}}_{3},tilde{ boldsymbol{P}}_{3}right) end{aligned}
p(xk+1,xk∣y1:T)=p(xk∣xk+1,y1:T)p(xk+1∣y1:T)=N(xk∣m~2,P~2)N(xk+1∣mk+1s,Pk+1s)=N([xk+1xk]∣m~3,P~3)
其中
m
~
3
=
[
m
k
+
1
s
m
k
+
G
k
(
m
k
+
1
s
−
A
k
m
k
)
]
P
~
3
=
[
P
k
+
1
s
P
k
+
1
s
G
k
T
G
k
P
k
+
1
s
G
k
P
k
+
1
s
G
k
T
+
P
~
2
]
tilde{ boldsymbol{m}}_{3} = left[ begin{array}{c} boldsymbol{m}^s_{k+1}\ boldsymbol{m}_{k} + boldsymbol G_k (boldsymbol m^s_{k+1} - boldsymbol A_k boldsymbol m_k) end{array} right] \ tilde{ boldsymbol P}_{3} = left[ begin{matrix} boldsymbol{P}^{s}_{k+1}& boldsymbol{P}^{s}_{k+1}boldsymbol{G}_{k}^{T}\ boldsymbol{G}_{k}boldsymbol{P}^{s}_{k+1}& boldsymbol{G}_{k}boldsymbol{P}^{s}_{k+1} boldsymbol{G}_{k}^{T} +tilde{boldsymbol{P}}_{2}\ end{matrix} right]
m~3=[mk+1smk+Gk(mk+1s−Akmk)]P~3=[Pk+1sGkPk+1sPk+1sGkTGkPk+1sGkT+P~2]
进一步,关于
x
k
boldsymbol{x}_k
xk的边际概率分布为
p
(
x
k
∣
y
1
:
T
)
=
N
(
x
k
∣
m
k
s
,
P
k
s
)
p(boldsymbol{x}_k | boldsymbol{y}_{1:T}) = mathcal{N}(boldsymbol{x}_k | boldsymbol{m}^s_k, boldsymbol{P}^s_k)
p(xk∣y1:T)=N(xk∣mks,Pks)
其中
m
k
s
=
m
k
+
G
k
(
m
k
+
1
s
−
A
k
m
k
)
P
k
s
=
P
k
+
G
k
(
P
k
s
−
A
k
P
k
A
k
T
−
Q
k
)
G
k
T
begin{aligned} boldsymbol{m}^s_k & = boldsymbol{m}_k + boldsymbol{G}_k (boldsymbol{m}^s_{k+1} - boldsymbol{A}_k boldsymbol{m}_k) \ boldsymbol{P}^s_k &= boldsymbol{P}_k + boldsymbol{G}_k ( boldsymbol{P}^s_k - boldsymbol{A}_k boldsymbol{P}_k boldsymbol{A}^T_k - boldsymbol{Q}_k ) boldsymbol{G}^T_k end{aligned}
mksPks=mk+Gk(mk+1s−Akmk)=Pk+Gk(Pks−AkPkAkT−Qk)GkT
这样也就完成了对式(11),关于RTS平滑的证明。
最后
以上就是乐观裙子为你收集整理的贝叶斯滤波和贝叶斯平滑(Kalman滤波,RTS平滑)贝叶斯滤波(Bayesian filtering equations)贝叶斯平滑(Bayesian smoothing)的全部内容,希望文章能够帮你解决贝叶斯滤波和贝叶斯平滑(Kalman滤波,RTS平滑)贝叶斯滤波(Bayesian filtering equations)贝叶斯平滑(Bayesian smoothing)所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复