@@ -18,7 +18,7 @@ kernelspec:
1818</div>
1919```
2020
21- # 卡尔曼滤波初探
21+ # 初见卡尔曼滤波器
2222
2323``` {index} single: 卡尔曼滤波
2424```
@@ -38,17 +38,17 @@ tags: [hide-output]
3838
3939## 概述
4040
41- 本讲座为卡尔曼滤波提供了一个简单直观的介绍 ,适合以下读者:
41+ 本讲座为卡尔曼滤波器提供了一个简单直观的介绍 ,适合以下读者:
4242
43- * 听说过卡尔曼滤波但不知道它如何工作的人 ,或者
44- * 知道卡尔曼滤波方程但不知道这些方程从何而来的人
43+ * 听说过卡尔曼滤波器但不知道它如何运作的人 ,或者
44+ * 知道卡尔曼滤波的方程但不知道这些方程从何而来的人
4545
4646关于卡尔曼滤波的更多(进阶)阅读材料,请参见:
4747
4848* {cite}` Ljungqvist2012 ` ,第2.7节
4949* {cite}` AndersonMoore2005 `
5050
51- 第二个参考文献对卡尔曼滤波进行了全面的阐述 。
51+ 第二个参考文献对卡尔曼滤波器进行了全面的阐述 。
5252
5353所需知识:熟悉矩阵运算、多元正态分布、协方差矩阵等。
5454
@@ -73,7 +73,7 @@ from scipy.linalg import eigvals
7373
7474## 基本概念
7575
76- 卡尔曼滤波在经济学中有许多应用 ,但现在让我们假装我们是火箭科学家。
76+ 卡尔曼滤波器在经济学中有许多应用 ,但现在让我们假装我们是火箭科学家。
7777
7878一枚导弹从Y国发射,我们的任务是追踪它。
7979
@@ -83,13 +83,15 @@ from scipy.linalg import eigvals
8383
8484总结我们知识的一种方式是点预测 $\hat x$
8585
86- * 但如果总统想知道导弹目前在日本海上空的概率呢?
87- * 那么用二元概率密度 $p$ 来总结我们的初始认知会更好
88- * $\int_E p(x)dx$ 表示我们认为导弹在区域 E 内的概率。
86+ 然而,点预测可能不够用。例如,我们可能需要回答"导弹目前在日本海上空的概率是多少"这样的问题。
87+
88+ 为了回答这类问题,我们需要用二元概率密度函数 $p$ 来描述我们对导弹位置的认知。
89+
90+ 对于任意区域 $E$,积分 $\int_E p(x)dx$ 给出了我们认为导弹在该区域内的概率。
8991
9092密度 $p$ 被称为随机变量 $x$ 的* 先验* 。
9193
92- 为了使我们的例子便于处理,我们假设我们的先验是高斯分布 。
94+ 为了使我们的例子便于处理,我们假设我们的先验分布是高斯分布 。
9395
9496特别地,我们采用
9597
@@ -127,7 +129,7 @@ p = N(\hat x, \Sigma)
127129---
128130tags: [output_scroll]
129131---
130- # 设置高斯先验密度 p
132+ # 设定高斯先验分布 p
131133Σ = [[0.4, 0.3], [0.3, 0.45]]
132134Σ = np.matrix(Σ)
133135x_hat = np.matrix([0.2, -0.2]).T
@@ -142,7 +144,7 @@ Q = 0.3 * Σ
142144# y 的观测值
143145y = np.matrix([2.3, -1.9]).T
144146
145- # 设置绘图网格
147+ # 设定绘图网格
146148x_grid = np.linspace(-1.5, 2.9, 100)
147149y_grid = np.linspace(-3.1, 1.7, 100)
148150X, Y = np.meshgrid(x_grid, y_grid)
@@ -186,7 +188,7 @@ def bivariate_normal(x, y, σ_x=1.0, σ_y=1.0, μ_x=0.0, μ_y=0.0, σ_xy=0.0):
186188
187189def gen_gaussian_plot_vals(μ, C):
188190 "用于绘制二元高斯 N(μ, C) 的 Z 值"
189- m_x, m_y = float(μ[0]) , float(μ[1])
191+ m_x, m_y = float(μ[0].item()) , float(μ[1].item() )
190192 s_x, s_y = np.sqrt(C[0, 0]), np.sqrt(C[1, 1])
191193 s_xy = C[0, 1]
192194 return bivariate_normal(X, Y, s_x, s_y, m_x, m_y, s_xy)
@@ -227,20 +229,19 @@ plt.show()
227229
228230坏消息是我们的传感器并不精确。
229231
230- 具体来说,我们应该将传感器的输出理解为不是
231- $y=x$,而是
232+ 具体来说,我们不应该将传感器的输出理解为$y=x$,而是
232233
233234``` {math}
234235:label: kl_measurement_model
235236
236- y = G x + v, \quad \text{where } \quad v \sim N(0, R)
237+ y = G x + v, \quad \text{且 } \quad v \sim N(0, R)
237238```
238239
239240这里 $G$ 和 $R$ 是 $2 \times 2$ 矩阵,其中 $R$ 是正定矩阵。两者都被假定为已知,且噪声项 $v$ 被假定与 $x$ 独立。
240241
241- 那么,我们应该如何将我们的先验 $p(x) = N(\hat x, \Sigma)$ 和这个新信息 $y$ 结合起来,以改进我们对导弹位置的理解呢 ?
242+ 那么,我们应该如何将我们的先验分布 $p(x) = N(\hat x, \Sigma)$ 和这个新信息 $y$ 结合起来,以提高我们对导弹位置的了解呢 ?
242243
243- 正如你可能已经猜到的 ,答案是使用贝叶斯定理,它告诉我们通过以下方式将先验 $p(x)$ 更新为 $p(x \, |\, y)$:
244+ 你可能已经猜到了 ,答案是使用贝叶斯定理,它告诉我们通过以下方式将先验分布 $p(x)$ 更新为 $p(x \, |\, y)$:
244245
245246$$
246247p(x \,|\, y) = \frac{p(y \,|\, x) \, p(x)} {p(y)}
257258
258259由于我们处在线性和高斯框架中,可以通过计算总体线性回归来得到更新后的密度。
259260
260- 具体来说,已知解 [ ^ f1 ] 为
261+ 具体来说,我们可以得出解 [ ^ f1 ] 为
261262
262263$$
263264p(x \,|\, y) = N(\hat x^F, \Sigma^F)
@@ -293,7 +294,7 @@ new_Z = gen_gaussian_plot_vals(x_hat_F, Σ_F)
293294cs2 = ax.contour(X, Y, new_Z, 6, colors="black")
294295ax.clabel(cs2, inline=1, fontsize=10)
295296ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
296- ax.text(float(y[0]) , float(y[1]), "$y$", fontsize=20, color="black")
297+ ax.text(float(y[0].item()) , float(y[1].item() ), "$y$", fontsize=20, color="black")
297298
298299plt.show()
299300```
@@ -307,13 +308,13 @@ plt.show()
307308
308309到目前为止我们取得了什么成果?
309310
310- 我们已经获得了基于先验和当前信息的状态 (导弹)当前位置的概率。
311+ 我们在给定先验分布和当前信息的情况下,已经获得了状态 (导弹)当前位置的概率。
311312
312313这被称为"滤波"而不是预测,因为我们是在过滤噪声而不是展望未来。
313314
314315* $p(x \, |\, y) = N(\hat x^F, \Sigma^F)$ 被称为* 滤波分布*
315316
316- 但现在假设我们有另一个任务:预测导弹在一个时间单位 (无论是什么单位)后的位置 。
317+ 但现在假设我们有另一个任务:预测导弹在一个时间单位后 (无论是什么单位)的位置 。
317318
318319为此我们需要一个状态演化的模型。
319320
@@ -322,7 +323,7 @@ plt.show()
322323``` {math}
323324:label: kl_xdynam
324325
325- x_{t+1} = A x_t + w_{t+1}, \quad \text{where } \quad w_t \sim N(0, Q)
326+ x_{t+1} = A x_t + w_{t+1}, \quad \text{且 } \quad w_t \sim N(0, Q)
326327```
327328
328329我们的目标是将这个运动定律和我们当前的分布 $p(x \, |\, y) = N(\hat x^F, \Sigma^F)$ 结合起来,得出一个新的一个时间单位后位置的* 预测* 分布。
@@ -406,7 +407,7 @@ new_Z = gen_gaussian_plot_vals(new_x_hat, new_Σ)
406407cs3 = ax.contour(X, Y, new_Z, 6, colors="black")
407408ax.clabel(cs3, inline=1, fontsize=10)
408409ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
409- ax.text(float(y[0]) , float(y[1]), "$y$", fontsize=20, color="black")
410+ ax.text(float(y[0].item()) , float(y[1].item() ), "$y$", fontsize=20, color="black")
410411
411412plt.show()
412413```
@@ -422,9 +423,9 @@ plt.show()
422423
423424然后我们使用当前测量值$y$更新为$p(x \, |\, y)$。
424425
425- 最后,我们使用$\{ x_t\} $的运动方程{eq}` kl_xdynam ` 更新为 $p_ {new}(x)$。
426+ 最后,我们使用$\{ x_t\} $的运动方程{eq}` kl_xdynam ` 将其更新为 $p_ {new}(x)$。
426427
427- 如果我们现在进入下一个周期,我们就可以再次循环,将$p_ {new}(x)$作为当前先验 。
428+ 如果我们现在进入下一个周期,我们就可以再次循环,将$p_ {new}(x)$作为当前的先验分布 。
428429
429430将符号$p_t(x)$替换为$p(x)$,将$p_ {t+1}(x)$替换为$p_ {new}(x)$,完整的递归程序为:
430431
@@ -473,17 +474,17 @@ plt.show()
473474
474475这是一个关于 $\Sigma_t$ 的非线性差分方程。
475476
476- {eq}` kalman_sdy ` 的不动点是满足以下条件的常数矩阵 $\Sigma$:
477+ {eq}` kalman_sdy ` 的固定点是满足以下条件的常数矩阵 $\Sigma$:
477478
478479``` {math}
479480:label: kalman_dare
480481
481482\Sigma = A \Sigma A' - A \Sigma G' (G \Sigma G' + R)^{-1} G \Sigma A' + Q
482483```
483484
484- 方程 {eq}` kalman_sdy ` 被称为离散时间里卡提差分方程 。
485+ 方程 {eq}` kalman_sdy ` 被称为离散时间黎卡提差分方程 。
485486
486- 方程 {eq}` kalman_dare ` 被称为[ 离散时间代数黎卡提方程] ( https://en.wikipedia.org/wiki/Algebraic_Riccati_equation ) 。
487+ 方程 {eq}` kalman_dare ` 被称为[ 离散时间代数黎卡提方程] ( https://zhuanlan.zhihu.com/p/692283143 ) 。
487488
488489关于固定点存在的条件以及序列 $\{ \Sigma_t\} $ 收敛到该固定点的条件在 {cite}` AHMS1996 ` 和 {cite}` AndersonMoore2005 ` 第4章中有详细讨论。
489490
@@ -498,10 +499,10 @@ plt.show()
498499``` {index} single: Kalman Filter; Programming Implementation
499500```
500501
501- 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类实现了卡尔曼滤波器
502+ [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类实现了卡尔曼滤波器
502503
503504* 实例数据包括:
504- * 当前先验的矩 $(\hat x_t, \Sigma_t)$
505+ * 当前先验分布的矩 $(\hat x_t, \Sigma_t)$
505506 * 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 的 [ LinearStateSpace] ( https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py ) 类的一个实例
506507
507508后者表示形式如下的线性状态空间模型
516517
517518其中冲击项 $w_t$ 和 $v_t$ 是独立同分布的标准正态分布。
518519
519- 为了与本讲座的符号保持一致 ,我们设定
520+ 为了与本章节的符号保持一致 ,我们设定
520521
521522$$
522523Q := CC' \quad \text{和} \quad R := HH'
523524$$
524525
525- * 来自 [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类有许多方法,其中一些我们会等到后续讲座中学习更高级的应用时再使用 。
526+ * [ QuantEcon.py] ( http://quantecon.org/quantecon-py ) 包的 ` Kalman ` 类有许多方法,其中一些我们会等到后续章节中学习更高级的应用时再使用 。
526527* 与本讲座相关的方法有:
527528
528- * ` prior_to_filtered ` ,将 $(\hat x_t, \Sigma_t)$ 更新为 $(\hat x_t^F, \Sigma_t^F)$
529+ * ` prior_to_filtered ` ,将 $(\hat x_t, \Sigma_t)$ 更新为 $(\hat x_t^F, \Sigma_t^F)$
529530 * ` filtered_to_forecast ` ,将滤波分布更新为预测分布 -- 成为新的先验分布 $(\hat x_ {t+1}, \Sigma_ {t+1})$
530531 * ` update ` ,结合上述两种方法
531532 * ` stationary_values ` ,计算{eq}` kalman_dare ` 的解和相应的(稳态)卡尔曼增益
574575A, C, G, H = 1, 0, 1, 1
575576ss = LinearStateSpace(A, C, G, H, mu_0=θ)
576577
577- # 设置先验 ,初始化卡尔曼滤波器
578+ # 设定先验分布 ,初始化卡尔曼滤波器
578579x_hat_0, Σ_0 = 8, 1
579580kalman = Kalman(ss, x_hat_0, Σ_0)
580581
@@ -583,7 +584,7 @@ N = 5
583584x, y = ss.simulate(N)
584585y = y.flatten()
585586
586- # 设置图形
587+ # 设定图形
587588fig, ax = plt.subplots(figsize=(10,8))
588589xgrid = np.linspace(θ - 5, θ + 2, 200)
589590
614615z_t := 1 - \int_{\theta - \epsilon}^{\theta + \epsilon} p_t(x) dx
615616$$
616617
617- 对于 $t = 0, 1, 2, \ldots, T$。
618+ 其中 $t = 0, 1, 2, \ldots, T$。
618619
619- 绘制 $z_t$ 与 $T$ 的关系图,设置 $\epsilon = 0.1$ 和 $T = 600$。
620+ 绘制 $z_t$ 与 $T$ 的关系图,设定 $\epsilon = 0.1$ 和 $T = 600$。
620621
621622你的图应该显示误差不规则地下降,类似这样
622623
@@ -647,7 +648,7 @@ y = y.flatten()
647648
648649for t in range(T):
649650 # 记录当前预测的均值和方差并绘制其密度
650- m, v = [float(temp) for temp in (kalman.x_hat, kalman.Sigma)]
651+ m, v = [float(temp.item() ) for temp in (kalman.x_hat, kalman.Sigma)]
651652
652653 f = lambda x: norm.pdf(x, loc=m, scale=np.sqrt(v))
653654 integral, error = quad(f, θ - ϵ, θ + ϵ)
@@ -670,21 +671,21 @@ plt.show()
670671:label: kalman_ex3
671672```
672673
673- 如{ref}` 上文所述 <kalman_convergence> ` ,如果冲击序列 $\{ w_t\} $ 不是退化的,那么在 $t-1$ 时刻通常无法无误地预测 $x_t$(即使我们能观察到 $x_ {t-1}$ 也是如此 )。
674+ 如{ref}` 上文所述 <kalman_convergence> ` ,如果冲击序列 $\{ w_t\} $ 不是退化的,那么在 $t-1$ 时刻通常无法无误地预测 $x_t$(即使我们能观察到 $x_ {t-1}$ ,情况也是如此 )。
674675
675- 让我们现在将卡尔曼滤波得到的预测值 $\hat x_t$ 与一个** 被允许** 观察 $x_ {t-1}$ 的竞争者进行比较。
676+ 让我们现在将在卡尔曼滤波器得到的预测值 $\hat x_t$ 与一个** 被允许** 观察 $x_ {t-1}$ 的竞争者进行比较。
676677
677678这个竞争者将使用条件期望 $\mathbb E[ x_t \, |\, x_ {t-1}] $,在这种情况下等于 $A x_ {t-1}$。
678679
679680条件期望被认为是在最小化均方误差方面的最优预测方法。
680681
681- (更准确地说,关于 $g$ 的 $ \mathbb E \, \| x_t - g(x_ {t-1}) \| ^2$ 的最小化器是 $g^* (x_ {t-1}) := \mathbb E[ x_t \, |\, x_ {t-1}] $)
682+ (更准确地说, $ \mathbb E \, \| x_t - g(x_ {t-1}) \| ^2$ 关于 $g$ 的最小值是 $g^* (x_ {t-1}) := \mathbb E[ x_t \, |\, x_ {t-1}] $)
682683
683- 因此,我们是在将卡尔曼滤波与一个拥有更多信息(在能够观察潜在状态的意义上 )的竞争者进行比较,并且
684+ 因此,我们是在将卡尔曼滤波器与一个拥有更多信息(能够观察潜在状态 )的竞争者进行比较,并且
684685
685686在最小化平方误差方面表现最优。
686687
687- 我们的对比竞赛将以平方误差来评估 。
688+ 我们的赛马式竞争将以平方误差来评估 。
688689
689690具体来说,你的任务是生成一个图表,绘制 $\| x_t - A x_ {t-1} \| ^2$ 和 $\| x_t - \hat x_t \| ^2$ 对 $t$ 的观测值,其中 $t = 1, \ldots, 50$。
690691
702703 \right)
703704$$
704705
705- 要初始化先验密度 ,设定
706+ 要初始化先验分布 ,设定
706707
707708$$
708709\Sigma_0
@@ -741,10 +742,10 @@ A = [[0.5, 0.4],
741742 [0.6, 0.3]]
742743C = np.sqrt(0.3) * np.identity(2)
743744
744- # 设置状态空间模型 ,初始值 x_0 设为零
745+ # 设定状态空间模型 ,初始值 x_0 设为零
745746ss = LinearStateSpace(A, C, G, H, mu_0 = np.zeros(2))
746747
747- # 定义先验密度
748+ # 定义先验分布
748749Σ = [[0.9, 0.3],
749750 [0.3, 0.9]]
750751Σ = np.array(Σ)
@@ -759,7 +760,7 @@ print(eigvals(A))
759760
760761# 打印平稳 Σ
761762S, K = kn.stationary_values()
762- print("平稳预测误差方差 :")
763+ print("平稳的预测误差方差 :")
763764print(S)
764765
765766# 生成图表
@@ -789,12 +790,12 @@ plt.show()
789790``` {exercise}
790791:label: kalman_ex4
791792
792- 尝试上下调整系数 $0.3$ (在 $ Q = 0.3 I$ 中) 。
793+ 尝试上下调整$ Q = 0.3 I$ 中的系数 $0.3$ 。
793794
794795观察平稳解 $\Sigma$ (参见 {eq}`kalman_dare`) 中的对角线值如何随这个系数增减而变化。
795796
796797这说明 $x_t$ 运动规律中的随机性越大,会导致预测中的(永久性)不确定性越大。
797798```
798799
799- [ ^ f1 ] : 例如,参见 {cite}` Bishop2006 ` 第93页。要从他的表达式得到上面使用的表达式,你还需要应用 [ Woodbury矩阵恒等式] ( https://en.wikipedia.org/wiki/Woodbury_matrix_identity ) 。
800+ [ ^ f1 ] : 例如,参见 {cite}` Bishop2006 ` 第93页。要从他的表达式得到上面使用的表达式,你还需要应用 [ Woodbury矩阵恒等式] ( https://zhuanlan.zhihu.com/p/388027547 ) 。
800801
0 commit comments