动态线性模型
一元动态线性模型(DLM)的一般形式为
$\begin{array}{rl}
\pmb y_t=&\!\!\!\!\pmb F_t\pmb\theta_t+\pmb\nu_t,\quad\pmb\nu_t\sim N_m(\pmb 0,\pmb V_t),\\
\pmb\theta_t=&\!\!\!\!\pmb G_t\pmb\theta_{t-1}+\pmb\omega_t,\quad\pmb\omega_t\sim N_p(\pmb 0,\pmb W_t),
\end{array}$
其中 $\pmb y_t$ 称为测量向量,$\pmb\theta_t$ 称为状态向量,$\pmb x_t=\pmb F_t\pmb\theta_t$ 表示信号,$\pmb F_t$ 为回归矩阵,$\pmb G_t$ 为状态矩阵,而 $\pmb\nu_t$ 和 $\pmb\omega_t$ 为零均值的测量误差和状态信息(innovations)向量,且相互独立。
一元动态线性模型是线性状态空间(state space)模型,时间序列常见的趋势和季节性特征都可以以这种形式进行建模。
平稳的 AR(2) 模型的一般表达式为
$y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+\varepsilon_t,\quad\varepsilon_t\overset{\text{i.i.d.}}{\sim} N(0,\tau^2).$
定义状态向量 $\pmb\theta_t=(y_t,y_{t-1})’$,则状态转移方程为
$\pmb\theta_t=\left(\begin{array}{c}
y_t\\ y_{t-1}
\end{array}\right)=\left(\begin{array}{cc}
\phi_1 & \phi_2\\ 1 & 0\\
\end{array}\right)\left(\begin{array}{c}
y_{t-1}\\ y_{t-2}
\end{array}\right)+\left(\begin{array}{c}
\varepsilon_t\\ 0
\end{array}\right),$
且测量方程为
$y_t=(1,0)\pmb\theta_t+\nu_t.$
动态线性回归模型的一般表达式为
$\begin{array}{rl}
y_t=&\!\!\!\!\pmb F_t\pmb\theta_t+\nu_t,\quad\nu_t\sim N(0,V_t),\\
\pmb\theta_t=&\!\!\!\!\pmb\theta_{t-1}+\pmb\omega_t,\quad\pmb\omega_t\sim N_p(\pmb 0,\pmb W_t).
\end{array}$
记 $\mathcal D_{t-1}={y_1,y_2,\dots,y_{t-1}}$ 为在时间 $t-1$ 时刻时候的累积信息,并假设
$\pmb\theta_{t-1}\mid\mathcal D_{t-1}\sim N_p(\pmb m_{t-1},\pmb C_{t-1})$,且 $\pmb\theta_0\sim N_p(\pmb m_0,\pmb C_0).$
于是 $\pmb\theta_t$ 的先验分布为
$\pmb\theta_t\mid\mathcal D_{t-1}\sim N_p(\pmb m_{t-1},\pmb R_t),\quad\pmb R_t=\pmb C_{t-1}+\pmb W_t.$
进一步 $y_t$ 的一步预测分布为
$y_t\mid\mathcal D_{t-1}\sim N(f_t,Q_t),\quad f_t=\pmb F_t\pmb m_{t-1},~Q_t=\pmb F_t\pmb R_t\pmb F_t'+V_t.$
它们的联合分布为
$\left(\begin{array}{c}
\pmb\theta_t\\ y_t
\end{array}\right)\mid\mathcal D_{t-1}\sim N_{p+1}\left(\left(\begin{array}{c}
\pmb m_{t-1}\\ f_t
\end{array}\right),\left(\begin{array}{cc}
\pmb R_t & \pmb R_t\pmb F_t'\\ \pmb F_t\pmb R_t & Q_t
\end{array}\right)\right).$
所以 $\pmb\theta_t$ 在给定 $\mathcal D_t={\mathcal D_{t-1},y_t}$ 的后验分布为
$\pmb\theta_t\mid\mathcal D_t\sim N_p(\pmb m_t,\pmb C_t),\quad\pmb m_t=\pmb m_{t-1}+\pmb A_te_t,~\pmb C_t=\pmb R_t-\pmb A_tQ_t\pmb A_t',$
其中 $\pmb A_t=\pmb R_t\pmb F_t’Q_t^{-1}$,$e_t=y_t-f_t$。
下面定理给出了一般的结果。
定理(过滤递归). (1) 状态的一步向前预报密度可以基于过滤密度 $p(\pmb\theta_{t-1}\mid\mathcal D_{t-1})$ 由下式给出
$\displaystyle p(\pmb\theta_t\mid\mathcal D_{t-1})=\int p(\pmb\theta_t\mid\pmb\theta_{t-1})p(\pmb\theta_{t-1}\mid\mathcal D_{t-1})\text{d}\nu(\pmb\theta_{t-1}).$
(2) 观测的一步向前预报密度可以基于状态的预测密度给出
$\displaystyle f(\pmb y_t\mid\mathcal D_{t-1})=\int f(\pmb y_t\mid\pmb\theta_t)p(\pmb\theta_t\mid\mathcal D_{t-1})\text{d}\nu(\pmb\theta_t).$
(3) 过滤密度可以基于前述密度给出
$p(\pmb\theta_t\mid\mathcal D_t)=\dfrac{f(\pmb y_t\mid\pmb\theta_t)p(\pmb\theta_t\mid\mathcal D_{t-1})}{f(\pmb y_t\mid\mathcal D_{t-1})}.$
上述结果可以用来递归计算 $k$ 步向前预报密度($k\geq1$):
$\displaystyle p(\pmb\theta_{t+k}\mid\mathcal D_t)=\int p(\pmb\theta_{t+k}\mid\pmb\theta_{t+k-1})p(\pmb\theta_{t+k-1}\mid\mathcal D_t)\text{d}\nu(\pmb\theta_{t+k-1}),$
以及
$\displaystyle f(\pmb y_{t+k}\mid\mathcal D_t)=\int f(\pmb y_{t+k}\mid\pmb\theta_{t+k})p(\pmb\theta_{t+k}\mid\mathcal D_t)\text{d}\nu(\pmb\theta_{t+k}).$
卡尔曼滤波
卡尔曼滤波(Kalman filter)是一种均方误差意义下最优线性状态估计方法。应用范围非常广泛:
- 在空中交通管制系统中使用雷达角度和范围的测量值来确定飞机接近机场的三维位置和速度;
- 使用电池单元的外部电气和热条件的测量值确定混合动力电动汽车电池组的内部“充电状态”和“健康状态”;
- 使用对地球磁场,GPS 和加速度的测量来确定物体在三个维度上的位置、方向和速度;
- 使用金融市场数据预测交易性证券和商品的未来价格。
卡尔曼滤波的核心思想是利用预测+测量反馈来达到最优估计,它是一种递归过程,主要有两个更新过程:
- 时间更新:主要包括状态预测和协方差预测,主要是对系统的预测;
- 观测更新:主要包括计算卡尔曼增益、状态更新和协方差更新。
下面定理给出了卡尔曼滤波的结果,它的证明实际上通过多元正态分布的一些简单计算即可。
定理(卡尔曼滤波). 对一元动态线性模型
$\begin{array}{rl}
y_t=&\!\!\!\!\pmb F_t\pmb\theta_t+\nu_t,\quad\nu_t\sim N(0,V_t),\\
\pmb\theta_t=&\!\!\!\!\pmb G_t\pmb\theta_{t-1}+\pmb H_{t-1}\pmb u_{t-1}+\pmb\omega_t,\quad\pmb\omega_t\sim N_p(\pmb 0,\pmb W_t),
\end{array}$
其中 $\pmb u_t$ 为系统控制输入变量,$\pmb H_t$ 为对应的系数矩阵。记 $\mathcal D_{t-1}={y_1,\dots,y_{t-1},\pmb u_1,\dots,\pmb u_{t-1}}$ 为在时间 $t-1$ 时刻时候的累积信息,并假设
$\pmb\theta_{t-1}\mid\mathcal D_{t-1}\sim N_p(\pmb m_{t-1},\pmb C_{t-1})$,且 $\pmb\theta_0\sim N_p(\pmb m_0,\pmb C_0).$
于是 $\pmb\theta_t$ 的先验分布为
$\pmb\theta_t\mid\mathcal D_{t-1}\sim N_p(\pmb a_t,\pmb R_t),\quad\pmb a_t=\pmb G_t\pmb m_{t-1}+\pmb H_{t-1}\pmb u_{t-1},~\pmb R_t=\pmb G_t\pmb C_{t-1}\pmb G_t'+\pmb W_t.$
进一步 $y_t$ 的一步预测分布为
$y_t\mid\mathcal D_{t-1}\sim N(f_t,Q_t),\quad f_t=\pmb F_t\pmb a_t,~Q_t=\pmb F_t\pmb R_t\pmb F_t'+V_t.$
从而 $\pmb\theta_t$ 的后验分布为
$\pmb\theta_t\mid\mathcal D_t\sim N(\pmb m_t,\pmb C_t),\quad\pmb m_t=\pmb a_t+\pmb A_te_t,~\pmb C_t=\pmb R_t-\pmb A_tQ_t\pmb A_t',$
其中 $\pmb A_t=\pmb R_t\pmb F_t’Q_t^{-1}$,$e_t=y_t-f_t$。利用Lec10中的矩阵恒等式可知
$\pmb C_t=\pmb R_t-\pmb A_tQ_t\pmb A_t'=\pmb R_t-\pmb R_t\pmb F_t'(\pmb F_t\pmb R_t\pmb F_t'+V_t)^{-1}\pmb F_t\pmb R_t'=(\pmb R_t^{-1}+\pmb F_t'V_t^{-1}\pmb F_t)^{-1}.$
根据定理结果有
光滑
在时间序列中,经常对 $y_t$ 观察一定时期,$t=1,\dots,T$,记 $\mathcal D_T=(y_1,\dots,y_T)$,则可以从 $p(\pmb\theta_T\mid\mathcal D_T)$ 出发使用向后迭代算法计算条件密度 $\pmb\theta_t\mid\mathcal D_T$,$t<T$。
定理(光滑递归). (1) 给定 $\mathcal D_T$,状态序列 $(\pmb\theta_0,\dots,\pmb\theta_T)$ 具有向后转移概率
$p(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)=\dfrac{p(\pmb\theta_{t+1}\mid\pmb\theta_t)p(\pmb\theta_t\mid\mathcal D_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)},\quad t=T-1,T-2,\dots,0.$
(2) $\pmb\theta_t\mid\mathcal D_T$,$t=T-1,T-2,\dots,1$ 的光滑密度可以由下式计算
$\displaystyle p(\pmb\theta_t\mid\mathcal D_T)=p(\pmb\theta_t\mid\mathcal D_t)\int\dfrac{p(\pmb\theta_{t+1}\mid\pmb\theta_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)}p(\pmb\theta_{t+1}\mid\mathcal D_T)\text{d}\pmb\theta_{t+1}.$
证明: (1) 注意到 $\pmb\theta_{t+1}\bot(\pmb\theta_0,\dots,\pmb\theta_{t-1})\mid\pmb\theta_t,\mathcal D_T$,以及 $\pmb\theta_t\bot(y_{t+1},\dots,y_T)\mid\pmb\theta_{t+1}$,由贝叶斯公式,
$\begin{array}{rl}
p(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)=&\!\!\!\!p(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_t)\\
=&\!\!\!\!\dfrac{p(\pmb\theta_t\mid\mathcal D_t)p(\pmb\theta_{t+1}\mid\pmb\theta_t,\mathcal D_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)}\\
=&\!\!\!\!\dfrac{p(\pmb\theta_t\mid\mathcal D_t)p(\pmb\theta_{t+1}\mid\pmb\theta_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)},\quad t=T-1,T-2,\dots,0.
\end{array}$
(2) 对 $p(\pmb\theta_t,\pmb\theta_{t+1}\mid\mathcal D_T)$ 边际化有
$\begin{array}{rl}
p(\pmb\theta_t\mid\mathcal D_T)=&\!\!\!\!\displaystyle\int p(\pmb\theta_t,\pmb\theta_{t+1}\mid\mathcal D_T)\text{d}\pmb\theta_{t+1}\\
=&\!\!\!\!\displaystyle\int p(\pmb\theta_{t+1}\mid\mathcal D_T)p(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_t)\text{d}\pmb\theta_{t+1}\\
=&\!\!\!\!\displaystyle p(\pmb\theta_t\mid\mathcal D_t)\int\dfrac{p(\pmb\theta_{t+1}\mid\pmb\theta_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)}p(\pmb\theta_{t+1}\mid\mathcal D_T)\text{d}\pmb\theta_{t+1}.
\end{array}$
定理. 对一元动态线性模型
$\begin{array}{rl}
\pmb y_t=&\!\!\!\!\pmb F_t\pmb\theta_t+\pmb\nu_t,\quad\pmb\nu_t\sim N_m(\pmb 0,\pmb V_t),\\
\pmb\theta_t=&\!\!\!\!\pmb G_t\pmb\theta_{t-1}+\pmb\omega_t,\quad\pmb\omega_t\sim N_p(\pmb 0,\pmb W_t),
\end{array}$
如果 $\pmb\theta_{t+1}\mid\mathcal D_T\sim N_p(\pmb s_{t+1},\pmb S_{t+1})$,则 $\pmb\theta_t\mid\mathcal D_T\sim N_p(\pmb s_t,\pmb S_t)$,其中
$\begin{array}{rl}
\pmb s_t=&\!\!\!\!\pmb m_t+\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}(\pmb s_{t+1}-\pmb a_{t+1}),\\
\pmb S_t=&\!\!\!\!\pmb C_t+\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}(\pmb S_{t+1}-\pmb R_{t+1})\pmb R_{t+1}^{-1}\pmb G_{t+1}\pmb C_t.
\end{array}$
证明: 由多元正态分布的性质,易知 $\pmb\theta_t\mid\mathcal D_T$ 的条件密度为正态分布,因此只需要计算其期望与协方差矩阵即可。注意到
$\pmb\theta_{t+1}\mid\pmb\theta_t\sim N_p(\pmb G_{t+1}\pmb\theta_t,\pmb W_{t+1}),\quad\pmb\theta_t\mid\mathcal D_t\sim N_p(\pmb m_t,\pmb C_t).$
于是
$\begin{array}{rl}
p(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)=&\!\!\!\!\dfrac{p(\pmb\theta_{t+1}\mid\pmb\theta_t)p(\pmb\theta_t\mid\mathcal D_t)}{p(\pmb\theta_{t+1}\mid\mathcal D_t)}\\
\propto&\!\!\!\!\exp\left\{-\dfrac{1}{2}\left[(\pmb\theta_{t+1}-\pmb G_{t+1}\pmb\theta_t)'\pmb W_{t+1}^{-1}(\pmb\theta_{t+1}-\pmb G_{t+1}\pmb\theta_t)+(\pmb\theta_t-\pmb m_t)'\pmb C_t^{-1}(\pmb\theta_t-\pmb m_t)\right]\right\}\\
\propto&\!\!\!\!\exp\left\{-\dfrac{1}{2}\left[\pmb\theta_t'(\pmb G_{t+1}'\pmb W_{t+1}^{-1}\pmb G_{t+1}+\pmb C_t^{-1})\pmb\theta_t-2(\pmb\theta_{t+1}'\pmb W_{t+1}^{-1}\pmb G_{t+1}+\pmb m_t'\pmb C_t^{-1})\pmb\theta_t\right]\right\}
\end{array}$
所以
$\begin{array}{\rl}
\mathbb E(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)=&\!\!\!\!\pmb m_t+\pmb C_t\pmb G_{t+1}'(\pmb G_{t+1}\pmb C_t\pmb G_{t+1}'+\pmb W_{t+1})^{-1}(\pmb\theta_{t+1}-\pmb G_{t+1}\pmb m_t)\\
=&\!\!\!\!\pmb m_t+\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}(\pmb\theta_{t+1}-\pmb a_{t+1}),\\
\text{Var}(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)=&\!\!\!\!(\pmb G_{t+1}'\pmb W_{t+1}^{-1}\pmb G_{t+1}+\pmb C_t^{-1})^{-1}\\
=&\!\!\!\!\pmb C_t-\pmb C_t\pmb G_{t+1}'\left(\pmb W_{t+1}+\pmb G_{t+1}\pmb C_t\pmb G_{t+1}'\right)^{-1}\pmb G_{t+1}\pmb C_t\\
=&\!\!\!\!\pmb C_t-\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}\pmb G_{t+1}\pmb C_t.
\end{array}$
因此
$\begin{array}{rl}
\pmb s_t=&\!\!\!\!\mathbb E[\mathbb E(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)\mid\mathcal D_T]=\pmb m_t+\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}(\pmb s_{t+1}-\pmb a_{t+1}),\\
\pmb S_t=&\!\!\!\!\text{Var}(\pmb\theta_t\mid\mathcal D_T)=\text{Var}[\mathbb E(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)\mid\mathcal D_T]+\mathbb E[\text{Var}(\pmb\theta_t\mid\pmb\theta_{t+1},\mathcal D_T)\mid\mathcal D_T]\\
=&\!\!\!\!\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}\pmb S_{t+1}\pmb R_{t+1}^{-1}\pmb G_{t+1}\pmb C_t+\pmb C_t-\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}\pmb G_{t+1}\pmb C_t\\
=&\!\!\!\!\pmb C_t+\pmb C_t\pmb G_{t+1}'\pmb R_{t+1}^{-1}(\pmb S_{t+1}-\pmb R_{t+1})\pmb R_{t+1}^{-1}\pmb G_{t+1}\pmb C_t.
\end{array}$
预报
很多情形下,我们还关心 $k$ 步向前预报的结果,即根据 $\mathcal D_t$ 预报 $y_{t+k}$ 的结果。
对 $k\geq1$,定义
$\pmb a_t(k)=\mathbb E(\pmb\theta_{t+k}\mid\mathcal D_t),\quad\pmb R_t(k)=\text{Var}(\pmb\theta_{t+k}\mid\mathcal D_t),\\
f_t(k)=\mathbb E(y_{t+k}\mid\mathcal D_t),\quad Q_t(k)=\text{Var}(y_{t+k}\mid\mathcal D_t).$
下面定理给出了预报的结论,事实上,$k$ 步预报就是按一步预报进行不断地迭代即可。
定理(预报). 令 $\pmb a_t(0)=\pmb m_t$,$\pmb R_t(0)=\pmb C_t$,则对 $k\geq1$,下面结论成立:
(1) $\pmb\theta_{t+k}$ 给定 $\mathcal D_t$ 的分布是正态分布,均值和协方差矩阵分别为
$\pmb a_t(k)=\pmb G_{t+k}\pmb a_t(k-1),\quad\pmb R_t(k)=\pmb G_{t+k}\pmb R_t(k-1)\pmb G_{t+k}'+\pmb W_{t+k}.$
(2) $y_{t+k}$ 给定 $\mathcal D_t$ 的条件分布是正态分布,均值和方差分别为
$f_t(k)=\pmb F_{t+k}\pmb a_t(k),\quad Q_t(k)=\pmb F_{t+k}\pmb R_t(k)\pmb F_{t+k}'+V_{t+k}.$
证明: 同样地,我们只需要计算各自的均值和协方差矩阵(或方差)即可。
下面我们对 $k$ 进行归纳。当 $k=1$ 时,这就是我们前面证明的结论;当 $k>1$ 时,有
$\begin{array}{rl}
\pmb a_t(k)=&\!\!\!\!\mathbb E(\pmb\theta_{t+k}\mid\mathcal D_t)=\mathbb E[\mathbb E(\pmb\theta_{t+k}\mid\pmb\theta_{t+k-1},\mathcal D_t)\mid\mathcal D_t]\\
=&\!\!\!\!\mathbb E(\pmb G_{t+k}\pmb\theta_{t+k-1}\mid\mathcal D_t)=\pmb G_{t+k}\pmb a_t(k-1),\\
\pmb R_t(k)=&\!\!\!\!\text{Var}(\pmb\theta_{t+k}\mid\mathcal D_t)\\
=&\!\!\!\!\text{Var}[\mathbb E(\pmb\theta_{t+k}\mid\pmb\theta_{t+k-1},\mathcal D_t)\mid\mathcal D_t]+\mathbb E[\text{Var}(\pmb\theta_{t+k}\mid\pmb\theta_{t+k-1},\mathcal D_t)\mid\mathcal D_t]\\
=&\!\!\!\!\pmb G_{t+k}\pmb R_t(k-1)\pmb G_{t+k}'+\pmb W_{t+k};\\
f_t(k)=&\!\!\!\!\mathbb E(y_{t+k}\mid\mathcal D_t)=\mathbb E[\mathbb E(y_{t+k}\mid\pmb\theta_{t+k},\mathcal D_t)\mid\mathcal D_t]\\
=&\!\!\!\!\mathbb E(\pmb F_{t+k}\pmb\theta_{t+k}\mid\mathcal D_t)=\pmb F_{t+k}\pmb a_t(k),\\
Q_t(k)=&\!\!\!\!\text{Var}(y_{t+k}\mid\mathcal D_t)\\
=&\!\!\!\!\text{Var}[\mathbb E(y_{t+k}\mid\pmb\theta_{t+k},\mathcal D_t)\mid\mathcal D_t]+\mathbb E[\text{Var}(y_{t+k}\mid\pmb\theta_{t+k},\mathcal D_t)\mid\mathcal D_t]\\
=&\!\!\!\!\pmb F_{t+k}\pmb R_t(k)\pmb F_{t+k}'+V_{t+k};
\end{array}$
应用:卡车位置与速度
考虑一辆在无摩擦的直线导轨上的卡车。最初卡车固定在位置0,但其受到随机不受控制的外力抖动。我们每 $\Delta t$ 秒测量一次卡车的位置,但是这些测量结果并不精确,考虑构建如下刻画卡车位置和速度的模型。
记位置和速度通过线性状态空间 $\pmb x_k=(x,\dot{x})$ 表示,其中 $x$ 为位置,$\dot{x}$ 为速度。假设在 $k-1$ 和 $k$ 时刻之间,不可控制外力导致一个常数加速度 $a_k$,其服从均值为0、方差为 $\sigma_a^2$ 的正态分布,因此
$\pmb x_k=\pmb F\pmb x_{k-1}+\pmb Ga_k,$
其中
$\pmb F=\left(\begin{array}{cc}
1 & \Delta t\\ 0 & 1
\end{array}\right),\quad\pmb G=\left(\begin{array}{c}
\dfrac{1}{2}(\Delta t)^2\\ \Delta t
\end{array}\right).$
等价地,有
$\pmb x_k=\pmb F\pmb x_{k-1}+\pmb w_k,\quad\pmb w_k\sim N(\pmb 0,\pmb Q),$
其中
$\pmb Q=\pmb G\pmb G'\sigma_a^2=\left(\begin{array}{cc}
\dfrac{1}{4}(\Delta t)^4 & \dfrac{1}{2}(\Delta t)^3\\ \dfrac{1}{2}(\Delta t)^3 & (\Delta t)^2
\end{array}\right)\sigma_a^2.$
注意到矩阵 $\pmb Q$ 不是满秩的,因此 $N(\pmb 0,\pmb Q)$ 不存在概率密度,避免这种情况的表示为 $\pmb w_k\sim\pmb G\cdot N(0,\sigma_a^2)$。
下面构建测量方程,在每一个时间步,卡车真实位置的一个带噪观测为
$\pmb z_k=\pmb H\pmb x_k+\pmb v_k,$
其中 $\pmb H=(1,0)$,$\pmb v_k\sim N(0,\sigma_v^2)$ 为随机噪音。
现在初始位置 $\hat{\pmb x}{0\vert 0}=(0,0)$,如果我们知道初始位置真实的位置和速度,可以指定一个值为0的初始协方差矩阵 $\pmb P{0\vert 0}$,否则可以指定一个合适的协方差矩阵,于是卡尔曼滤波的步骤为
$\begin{array}{rl}
\hat{\pmb x}_{k\vert k}=&\!\!\!\!\pmb F_k\hat{\pmb x}_{k-1\vert k-1}+\pmb K_k\left[\pmb z_k-\pmb H_k\pmb F_k\hat{\pmb x}_{k-1\vert k-1}\right],\\
\pmb P_{k\vert k-1}=&\!\!\!\!\pmb F_k\pmb P_{k-1\vert k-1}\pmb F_k'+\pmb Q_k,\\
\pmb S_k=&\!\!\!\!\pmb R_k+\pmb H_k\pmb P_{k\vert k-1}\pmb H_k',\\
\pmb K_k=&\!\!\!\!\pmb P_{k\vert k-1}\pmb H_k'\pmb S_k^{-1},\\
\pmb P_{k\vert k}=&\!\!\!\!(\pmb I-\pmb K_k\pmb H_k)\pmb P_{k\vert k-1}.
\end{array}$
后记
本文内容参考自中国科学技术大学《贝叶斯分析》(2022Fall)课件。
如对本文内容有任何疑问,请联系 watthu@mail.ustc.edu.cn。