@@ -18,7 +18,7 @@ kernelspec:
18
18
</div>
19
19
```
20
20
21
- # 卡尔曼滤波初探
21
+ # 卡尔曼滤波器的初步介绍
22
22
23
23
``` {index} single: 卡尔曼滤波
24
24
```
@@ -38,17 +38,17 @@ tags: [hide-output]
38
38
39
39
## 概述
40
40
41
- 本讲座为卡尔曼滤波提供了一个简单直观的介绍 ,适合以下读者:
41
+ 本讲座为卡尔曼滤波器提供了一个简单直观的介绍 ,适合以下读者:
42
42
43
- * 听说过卡尔曼滤波但不知道它如何工作的人 ,或者
44
- * 知道卡尔曼滤波方程但不知道这些方程从何而来的人
43
+ * 听说过卡尔曼滤波器但不知道它如何工作的人 ,或者
44
+ * 知道卡尔曼滤波的方程但不知道这些方程从何而来的人
45
45
46
46
关于卡尔曼滤波的更多(进阶)阅读材料,请参见:
47
47
48
48
* {cite}` Ljungqvist2012 ` ,第2.7节
49
49
* {cite}` AndersonMoore2005 `
50
50
51
- 第二个参考文献对卡尔曼滤波进行了全面的阐述 。
51
+ 第二个参考文献对卡尔曼滤波器进行了全面的阐述 。
52
52
53
53
所需知识:熟悉矩阵运算、多元正态分布、协方差矩阵等。
54
54
@@ -73,7 +73,7 @@ from scipy.linalg import eigvals
73
73
74
74
## 基本概念
75
75
76
- 卡尔曼滤波在经济学中有许多应用 ,但现在让我们假装我们是火箭科学家。
76
+ 卡尔曼滤波器在经济学中有许多应用 ,但现在让我们假装我们是火箭科学家。
77
77
78
78
一枚导弹从Y国发射,我们的任务是追踪它。
79
79
@@ -89,7 +89,7 @@ from scipy.linalg import eigvals
89
89
90
90
密度 $p$ 被称为随机变量 $x$ 的* 先验* 。
91
91
92
- 为了使我们的例子便于处理,我们假设我们的先验是高斯分布 。
92
+ 为了使我们的例子便于处理,我们假设我们的先验分布是高斯分布 。
93
93
94
94
特别地,我们采用
95
95
@@ -127,7 +127,7 @@ p = N(\hat x, \Sigma)
127
127
---
128
128
tags: [output_scroll]
129
129
---
130
- # 设置高斯先验密度 p
130
+ # 设定高斯先验分布 p
131
131
Σ = [[0.4, 0.3], [0.3, 0.45]]
132
132
Σ = np.matrix(Σ)
133
133
x_hat = np.matrix([0.2, -0.2]).T
@@ -142,7 +142,7 @@ Q = 0.3 * Σ
142
142
# y 的观测值
143
143
y = np.matrix([2.3, -1.9]).T
144
144
145
- # 设置绘图网格
145
+ # 设定绘图网格
146
146
x_grid = np.linspace(-1.5, 2.9, 100)
147
147
y_grid = np.linspace(-3.1, 1.7, 100)
148
148
X, Y = np.meshgrid(x_grid, y_grid)
@@ -227,8 +227,7 @@ plt.show()
227
227
228
228
坏消息是我们的传感器并不精确。
229
229
230
- 具体来说,我们应该将传感器的输出理解为不是
231
- $y=x$,而是
230
+ 具体来说,我们不应该将传感器的输出理解为$y=x$,而是
232
231
233
232
``` {math}
234
233
:label: kl_measurement_model
@@ -238,9 +237,9 @@ y = G x + v, \quad \text{where} \quad v \sim N(0, R)
238
237
239
238
这里 $G$ 和 $R$ 是 $2 \times 2$ 矩阵,其中 $R$ 是正定矩阵。两者都被假定为已知,且噪声项 $v$ 被假定与 $x$ 独立。
240
239
241
- 那么,我们应该如何将我们的先验 $p(x) = N(\hat x, \Sigma)$ 和这个新信息 $y$ 结合起来,以改进我们对导弹位置的理解呢 ?
240
+ 那么,我们应该如何将我们的先验分布 $p(x) = N(\hat x, \Sigma)$ 和这个新信息 $y$ 结合起来,以提高我们对导弹位置的了解呢 ?
242
241
243
- 正如你可能已经猜到的 ,答案是使用贝叶斯定理,它告诉我们通过以下方式将先验 $p(x)$ 更新为 $p(x \, |\, y)$:
242
+ 你可能已经猜到了 ,答案是使用贝叶斯定理,它告诉我们通过以下方式将先验分布 $p(x)$ 更新为 $p(x \, |\, y)$:
244
243
245
244
$$
246
245
p(x \,|\, y) = \frac{p(y \,|\, x) \, p(x)} {p(y)}
257
256
258
257
由于我们处在线性和高斯框架中,可以通过计算总体线性回归来得到更新后的密度。
259
258
260
- 具体来说,已知解 [ ^ f1 ] 为
259
+ 具体来说,我们可以得出解 [ ^ f1 ] 为
261
260
262
261
$$
263
262
p(x \,|\, y) = N(\hat x^F, \Sigma^F)
@@ -307,13 +306,13 @@ plt.show()
307
306
308
307
到目前为止我们取得了什么成果?
309
308
310
- 我们已经获得了基于先验和当前信息的状态 (导弹)当前位置的概率。
309
+ 我们在给定先验分布和当前信息的情况下,已经获得了状态 (导弹)当前位置的概率。
311
310
312
311
这被称为"滤波"而不是预测,因为我们是在过滤噪声而不是展望未来。
313
312
314
313
* $p(x \, |\, y) = N(\hat x^F, \Sigma^F)$ 被称为* 滤波分布*
315
314
316
- 但现在假设我们有另一个任务:预测导弹在一个时间单位 (无论是什么单位)后的位置 。
315
+ 但现在假设我们有另一个任务:预测导弹在一个时间单位后 (无论是什么单位)的位置 。
317
316
318
317
为此我们需要一个状态演化的模型。
319
318
@@ -422,9 +421,9 @@ plt.show()
422
421
423
422
然后我们使用当前测量值$y$更新为$p(x \, |\, y)$。
424
423
425
- 最后,我们使用$\{ x_t\} $的运动方程{eq}` kl_xdynam ` 更新为 $p_ {new}(x)$。
424
+ 最后,我们使用$\{ x_t\} $的运动方程{eq}` kl_xdynam ` 将其更新为 $p_ {new}(x)$。
426
425
427
- 如果我们现在进入下一个周期,我们就可以再次循环,将$p_ {new}(x)$作为当前先验 。
426
+ 如果我们现在进入下一个周期,我们就可以再次循环,将$p_ {new}(x)$作为当前的先验分布 。
428
427
429
428
将符号$p_t(x)$替换为$p(x)$,将$p_ {t+1}(x)$替换为$p_ {new}(x)$,完整的递归程序为:
430
429
@@ -473,15 +472,15 @@ plt.show()
473
472
474
473
这是一个关于 $\Sigma_t$ 的非线性差分方程。
475
474
476
- {eq}` kalman_sdy ` 的不动点是满足以下条件的常数矩阵 $\Sigma$:
475
+ {eq}` kalman_sdy ` 的固定点是满足以下条件的常数矩阵 $\Sigma$:
477
476
478
477
``` {math}
479
478
:label: kalman_dare
480
479
481
480
\Sigma = A \Sigma A' - A \Sigma G' (G \Sigma G' + R)^{-1} G \Sigma A' + Q
482
481
```
483
482
484
- 方程 {eq}` kalman_sdy ` 被称为离散时间里卡提差分方程 。
483
+ 方程 {eq}` kalman_sdy ` 被称为离散时间黎卡提差分方程 。
485
484
486
485
方程 {eq}` kalman_dare ` 被称为[ 离散时间代数黎卡提方程] ( https://en.wikipedia.org/wiki/Algebraic_Riccati_equation ) 。
487
486
@@ -498,10 +497,10 @@ plt.show()
498
497
``` {index} single: Kalman Filter; Programming Implementation
499
498
```
500
499
501
- 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类实现了卡尔曼滤波器
500
+ [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类实现了卡尔曼滤波器
502
501
503
502
* 实例数据包括:
504
- * 当前先验的矩 $(\hat x_t, \Sigma_t)$
503
+ * 当前先验分布的矩 $(\hat x_t, \Sigma_t)$
505
504
* 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 的 [ LinearStateSpace] ( https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py ) 类的一个实例
506
505
507
506
后者表示形式如下的线性状态空间模型
516
515
517
516
其中冲击项 $w_t$ 和 $v_t$ 是独立同分布的标准正态分布。
518
517
519
- 为了与本讲座的符号保持一致 ,我们设定
518
+ 为了与本章节的符号保持一致 ,我们设定
520
519
521
520
$$
522
521
Q := CC' \quad \text{和} \quad R := HH'
523
522
$$
524
523
525
- * 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类有许多方法,其中一些我们会等到后续讲座中学习更高级的应用时再使用 。
524
+ * [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类有许多方法,其中一些我们会等到后续章节中学习更高级的应用时再使用 。
526
525
* 与本讲座相关的方法有:
527
526
528
- * ` prior_to_filtered ` ,将 $(\hat x_t, \Sigma_t)$ 更新为 $(\hat x_t^F, \Sigma_t^F)$
527
+ * ` prior_to_filtered ` ,将 $(\hat x_t, \Sigma_t)$ 更新为 $(\hat x_t^F, \Sigma_t^F)$
529
528
* ` filtered_to_forecast ` ,将滤波分布更新为预测分布 -- 成为新的先验分布 $(\hat x_ {t+1}, \Sigma_ {t+1})$
530
529
* ` update ` ,结合上述两种方法
531
530
* ` stationary_values ` ,计算{eq}` kalman_dare ` 的解和相应的(稳态)卡尔曼增益
574
573
A, C, G, H = 1, 0, 1, 1
575
574
ss = LinearStateSpace(A, C, G, H, mu_0=θ)
576
575
577
- # 设置先验 ,初始化卡尔曼滤波器
576
+ # 设定先验分布 ,初始化卡尔曼滤波器
578
577
x_hat_0, Σ_0 = 8, 1
579
578
kalman = Kalman(ss, x_hat_0, Σ_0)
580
579
@@ -583,7 +582,7 @@ N = 5
583
582
x, y = ss.simulate(N)
584
583
y = y.flatten()
585
584
586
- # 设置图形
585
+ # 设定图形
587
586
fig, ax = plt.subplots(figsize=(10,8))
588
587
xgrid = np.linspace(θ - 5, θ + 2, 200)
589
588
614
613
z_t := 1 - \int_{\theta - \epsilon}^{\theta + \epsilon} p_t(x) dx
615
614
$$
616
615
617
- 对于 $t = 0, 1, 2, \ldots, T$。
616
+ 其中 $t = 0, 1, 2, \ldots, T$。
618
617
619
- 绘制 $z_t$ 与 $T$ 的关系图,设置 $\epsilon = 0.1$ 和 $T = 600$。
618
+ 绘制 $z_t$ 与 $T$ 的关系图,设定 $\epsilon = 0.1$ 和 $T = 600$。
620
619
621
620
你的图应该显示误差不规则地下降,类似这样
622
621
@@ -670,21 +669,21 @@ plt.show()
670
669
:label: kalman_ex3
671
670
```
672
671
673
- 如{ref}` 上文所述 <kalman_convergence> ` ,如果冲击序列 $\{ w_t\} $ 不是退化的,那么在 $t-1$ 时刻通常无法无误地预测 $x_t$(即使我们能观察到 $x_ {t-1}$ 也是如此 )。
672
+ 如{ref}` 上文所述 <kalman_convergence> ` ,如果冲击序列 $\{ w_t\} $ 不是退化的,那么在 $t-1$ 时刻通常无法无误地预测 $x_t$(即使我们能观察到 $x_ {t-1}$ ,情况也是如此 )。
674
673
675
- 让我们现在将卡尔曼滤波得到的预测值 $\hat x_t$ 与一个** 被允许** 观察 $x_ {t-1}$ 的竞争者进行比较。
674
+ 让我们现在将在卡尔曼滤波器得到的预测值 $\hat x_t$ 与一个** 被允许** 观察 $x_ {t-1}$ 的竞争者进行比较。
676
675
677
676
这个竞争者将使用条件期望 $\mathbb E[ x_t \, |\, x_ {t-1}] $,在这种情况下等于 $A x_ {t-1}$。
678
677
679
678
条件期望被认为是在最小化均方误差方面的最优预测方法。
680
679
681
- (更准确地说,关于 $g$ 的 $ \mathbb E \, \| x_t - g(x_ {t-1}) \| ^2$ 的最小化器是 $g^* (x_ {t-1}) := \mathbb E[ x_t \, |\, x_ {t-1}] $)
680
+ (更准确地说, $ \mathbb E \, \| x_t - g(x_ {t-1}) \| ^2$ 关于 $g$ 的最小值是 $g^* (x_ {t-1}) := \mathbb E[ x_t \, |\, x_ {t-1}] $)
682
681
683
- 因此,我们是在将卡尔曼滤波与一个拥有更多信息(在能够观察潜在状态的意义上 )的竞争者进行比较,并且
682
+ 因此,我们是在将卡尔曼滤波器与一个拥有更多信息(能够观察潜在状态 )的竞争者进行比较,并且
684
683
685
684
在最小化平方误差方面表现最优。
686
685
687
- 我们的对比竞赛将以平方误差来评估 。
686
+ 我们的赛马式竞争将以平方误差来评估 。
688
687
689
688
具体来说,你的任务是生成一个图表,绘制 $\| x_t - A x_ {t-1} \| ^2$ 和 $\| x_t - \hat x_t \| ^2$ 对 $t$ 的观测值,其中 $t = 1, \ldots, 50$。
690
689
702
701
\right)
703
702
$$
704
703
705
- 要初始化先验密度 ,设定
704
+ 要初始化先验分布 ,设定
706
705
707
706
$$
708
707
\Sigma_0
@@ -741,10 +740,10 @@ A = [[0.5, 0.4],
741
740
[0.6, 0.3]]
742
741
C = np.sqrt(0.3) * np.identity(2)
743
742
744
- # 设置状态空间模型 ,初始值 x_0 设为零
743
+ # 设定状态空间模型 ,初始值 x_0 设为零
745
744
ss = LinearStateSpace(A, C, G, H, mu_0 = np.zeros(2))
746
745
747
- # 定义先验密度
746
+ # 定义先验分布
748
747
Σ = [[0.9, 0.3],
749
748
[0.3, 0.9]]
750
749
Σ = np.array(Σ)
@@ -759,7 +758,7 @@ print(eigvals(A))
759
758
760
759
# 打印平稳 Σ
761
760
S, K = kn.stationary_values()
762
- print("平稳预测误差方差 :")
761
+ print("平稳的预测误差方差 :")
763
762
print(S)
764
763
765
764
# 生成图表
@@ -789,7 +788,7 @@ plt.show()
789
788
``` {exercise}
790
789
:label: kalman_ex4
791
790
792
- 尝试上下调整系数 $0.3$ (在 $ Q = 0.3 I$ 中) 。
791
+ 尝试上下调整$ Q = 0.3 I$ 中的系数 $0.3$ 。
793
792
794
793
观察平稳解 $\Sigma$ (参见 {eq}`kalman_dare`) 中的对角线值如何随这个系数增减而变化。
795
794
0 commit comments