g-h 过滤器(α-β过滤器) 上
在开始之前你要明确如是使用 jupyter Notebooks, 还要熟悉 Scipy、NumPy、Matplotilb 软件包,在这本书中你会经常使用到它们。前言部分介绍了这些软件包。(本书前言部分目前我没有进行翻译,并且这些软件包目前先知道就行,后面跟着用就慢慢会了,再不行就看看这些软件包的介绍就可以了)。
通过思维实验建立直觉
想象我们生活在一个没有秤的世界中–秤就是你站在上面称体重的装置。一天工作时,一位同事跑到你面前,她告诉你她发明了“秤”。在她告诉你怎么用之后,你站到秤上测量了一下,秤显示你目前重 172 磅。你特别的兴奋,因为你第一次知道了你自己的体重。更重要的是,当你想象着把这个设备卖给世界各地的减肥诊所时,你的眼睛里就会出现金钱符号!这太棒了!
(注:一磅大概是 0.45 公斤)
另一个同事听到了骚动,过来看看是什么让你这么兴奋。你解释了自己的发明,然后再次站到秤上,自豪地宣布结果:“161磅。”然后你犹豫并困惑了。
“前几秒读数还是 172 磅呢”,你向你的同事抱怨。
她解释道:“我从来没有说过这玩意是精准的呀!”。
传感器都是不精准的,这就是大量对滤波算法研究的意义,并且本书的意义就是解决这些问题的。我可以提供过去半个世纪发展起来的解决方案,但这些解决方案是通过问一些非常基本的问题来发展的,这些问题是关于我们所知的以及我们是如何知道的。在我们尝试讨论数学相关之前,让我们跟随这个探索之旅,看看它是否做到了符合我们直觉的过滤。
尝试使用另外的秤
我们有什么办法可以改进这个结果吗?显而易见,首先要做的就是找一个更好的传感器。不幸的是,你的同事告诉你,她已经做了10个秤,而且它们的操作精度都差不多。你让她拿出另一个秤,你用一个秤称体重,然后用另一个秤称。第一个秤(A)显示“160磅”,第二个秤(B)显示“170磅”。关于你的体重我们能得出什么结论?
好的,以下就是我们的选择:
- 我们选择只相信秤(A),然后给出我们的评估值 160 磅。
- 我们选择只相信秤(B),然后给出我们的评估值 170 磅。
- 我们只相信显示数值小的秤(哈哈哈,很符合大多数人的心理)。
- 我们只相信显示数值小的秤。
- 我们在 A 和 B 之间选择一个值。
前两种选择似乎是合理的,但我们没有理由偏爱其中一种。为什么我们会选择相信A而不是B呢?我们没有理由相信这一点。第三和第四个选择是非理性的。诚然,这两种测量并不十分准确,但完全没有理由选择一个超出它们所测量范围的数字。最后的选择是唯一合理的选择。如果两种体重秤都不准确,而且给出的结果可能高于实际体重,也可能低于实际体重,那么答案往往介于a和B之间。
在数学中,这个概念被形式化为期望值,我们将在后面深入讨论它。现在问问你自己,如果我们读了一百万次,“通常”会发生什么。有时两个秤的读数都过低,有时两个秤的读数都过高,而其他时候它们会跨越实际重量。如果它们超过了实际的重量,那么我们当然应该在a和B之间选择一个值。如果它们不超越,那么我们不知道它们是太高还是太低,但是通过在a和B之间选择一个数字,我们至少可以减轻最坏均值的影响。
后面我们会讨论出更正式的的结论,但是现在下面的就是清晰的给出最好的评估值,那就是 A 和 B 取平均值。
160
+
170
2
=
165
\frac{160+170}{2}=165
2160+170=165
我们可以从图形上看。我绘制了A和B的测量值,假设误差为±8磅。测量值在160到170磅之间,所以唯一有意义的体重必须在160到170磅之间。
说说我是如何生成这张图的。我从kf_book子目录中的模块book_plots导入代码。生成这个图需要大量的Python样板文件,读起来并不有趣。我在书中经常用这个图钉。当单元格运行时,调用plot_errorbars()并将情节插入到书中。
如果这是你第一次使用Jupyter Notebook,那么上面的代码位于单元格中。文本“In[2]:”将此单元格标记为可以输入输入的单元格,括号中的数字表示该单元格是第二个运行的单元格。要运行单元格,请用鼠标单击它,使其具有焦点,然后按键盘上的CTRL+ENTER。随着我们的继续,您将能够更改单元格内的代码并重新运行它们。尝试将值“160”,“170”和“8”更改为其他值并运行单元格。打印的输出应该根据您输入的内容而变化。
如果您想查看plot_error的代码可以在编辑器中打开它,或者创建一个新单元格,并键入函数名后跟两个问号。按Ctrl+Enter,浏览器将打开一个显示源代码的窗口。这是Jupyter Notebooks的一个特性。如果您只想查看该函数的文档,请执行相同的操作,但要加上一个问号。
plot_errorbars??
或者
plot_errorbars?
所以165磅看起来是一个合理的估计,但是这里有更多的信息我们可以利用。唯一可能的重量是在A和B误差条的交叉点上。例如,161磅的重量是不可能的,因为秤B不能在最大误差为8磅的情况下给出170磅的读数。同样,169磅的重量是不可能的,因为秤a不能在最大误差8磅的情况下给出160磅的读数。在这个例子中,唯一可能的重量在162到168磅之间。
这还不能让我们找到一个更好的重量的评估值,但让我们再玩一下“如果”。如果我们现在被告知 A 的准确率是 B 的三倍呢? 考虑我们上面列出的5个选项。选择a和B范围之外的数字仍然没有意义,所以我们不考虑这些。选择A作为我们获取到的评估值似乎更有说服力——毕竟,我们知道它更准确,为什么不用它来代替 B 呢? 再使用 B 的值能比 A 单独提高准确性吗?
答案可能与直觉相反,是可以的。首先,让我们看看 A=160 和 B=170 的相同测量值,但误差为 A±3 磅,B的误差是 3 倍,±9磅。
A 和 B 误差条的重叠是唯一可能的真正的重量。这个重叠小于单独 A 的误差。更重要的是,在这种情况下,我们可以看到重叠不包括 160 磅或 165 磅。如果我们只使用 A 的测量值,因为它比B 更 准确,我们会给出160磅的估计值。如果取 A 和 B 的平均值,得到165磅。考虑到我们对天平精度的了解,这两种重量都不可能。通过包括B的测量,我们将给出一个介于 161 磅和 163 磅之间的估计值,这是两个误差条相交的极限。
让我们来看看极限。假设我们知道秤 A 精确到1磅。换句话说,如果我们真的重 170 磅,它可能会报告 169、170 或 171 磅。我们还知道秤 B 精确到 9 磅。我们在每个磅秤上称重,得到a =160, B=170。我们应该估计自己的体重是多少?让我们来图解一下。
这里我们可以看到唯一可能的重量是161磅。这是一个重要的结果。用两个相对不准确的传感器,我们可以推断出极其准确的结果。
所以两个传感器,即使一个不如另一个精确,也比只用一个传感器好。 在这本书的其余部分,我将反复强调这一点。我们从不丢弃信息,不管它多么看起来没用。我们将利用数学的相关算法,使我们能够包括所有可能的信息来源,以形成尽可能最好的估计值。
然而,我们偏离了我们的问题。没有客户会想要购买好几个秤,此外,我们最初假设所有的秤都是同样准确的。这种使用所有测量而不考虑准确性的见解将在后面发挥重要作用,所以不要忘记它。
如果我只有一个体重秤,但我要称很多次呢? 我们得出结论,如果我们有两个精度相等的刻度,我们就应该取它们测量结果的平均值。如果我用一个秤称一万次体重呢?我们已经说过,刻度返回一个太大的数字和返回一个太小的数字的可能性是一样的。证明大量重量的平均值将非常接近实际重量并不难,所以现在让我们编写一个模拟。我将使用NumPy,这是用于数值计算的SciPy生态系统的一部分。
打印的确切数字取决于你的随机数生成器的散列情况,但它应该非常接近165。
这段代码做了一个可能不正确的假设——对于165磅的真实体重,体重秤上的读数可能是 160 和 165 但这几乎都不是真实的值。真正的传感器更有可能获得接近真实值的读数,并且越来越不可能获得远离真实值的读数。我们将在高斯(Gaussian )的章节中详细讨论这一点。现在,我将不作进一步解释地使用 numpy.random.normal()
函数,它将在 165 磅附近产生更多的值,在更远的地方产生更少的值。现在请相信,这将产生类似于实际规模的噪声测量结果。
再一次的结果在 165 附近。
太好了,我们有传感器问题的答案了!但这不是一个非常实际的答案。没有人有耐心称自己一万次,甚至十几次。
所以,让我们在换一种情况。如果你每天量一次体重,得到的读数分别是170、161和169。你是增重了,还是减重了,还是这些都只是嘈杂的测量?
我们真的不能说。第一次测量是170,最后一次是169,这意味着减掉了1磅。但如果体重秤只能精确到10磅,那就可以用噪声来解释了。我本来可以增重的;也许我第一天的体重是165磅,第三天是172磅。
有可能在体重增加的情况下得到这些体重读数。我的体重秤告诉我我在减肥,而实际上我在增重!让我们用图表来看一下。我将测量结果与误差条一起绘制出来,然后用绿线虚线表示一些可能的体重增加/减少。
正如我们所看到的,体重变化的极端范围可以用这三种测量来解释。事实上,有无数的选择。我们应该放弃吗? 不! 回想一下,我们正在谈论测量一个人的体重。一个人在第一天体重180磅,在第三天体重160磅是不可能的。或者一天减掉30磅,第二天又重回来(我们假设这个人没有截肢或其他创伤)。
我们正在测量的物理系统的行为直接影响我们如何解释测量结果。如果我们每天给一块石头称重,我们会把所有的变化都归因于噪声。如果我们给一个收集雨水用来做家务的水箱称重,我们可能会相信这种重量变化是真实的。
假设我取一个不同的刻度,我得到以下的测量值:169 170 169 171 171 171 169 170 170 169 169 170。你的直觉告诉你什么?例如,有可能你每天增重1磅,而含有噪声的测量结果恰好看起来你的体重保持不变。同样地,你可以一天减掉1磅,得到同样的读数。但这可能吗? 抛硬币连续得到10次正面的可能性有多大? 这实际是不太可能的。我们不能仅仅根据这些读数来证明这一点,但我的体重似乎很可能保持稳定。在下面的图表中,我用误差条绘制了测量结果,绿色虚线表示可能的真实重量。这条虚线并不意味着是这个问题的“正确”答案,只是一个合理的,可以通过测量来解释的答案。
另一个假设是:如果读数是158.0 164.2 160.3 159.9 162.1 164.6 169.6 167.4 166.4 171.0?让我们看一个图表,然后解释一些问题。
我是不是“看起来”很可能瘦了,而这只是一个充满噪声污染的数据? 不是真的。我的体重看起来和你一样重吗? 再一次,没有。随着时间的推移,这些数据呈上升趋势;不是均匀的,但肯定是上升的。我们不能确定,但这看起来像是体重增加了,而且是明显的体重增加。让我们用更多的图来验证这个假设。在图表中“观察”数据通常比在表格中更容易。
我们来看两个假设。首先,让我们假设体重没有变化。为了得到这个数字,我们同意对测量结果取平均值。我们来看看。
这看起来不太令人信服。事实上,我们可以看到,在所有的误差条内,我们无法画出一条水平线。
现在,假设我们的体重增加了。增加了多少?我不知道,但是NumPy知道! 我们想通过测量画一条看起来“大致”正确的线。NumPy有一些函数会根据“最小二乘拟合”的规则来做这件事。让我们不要担心计算的细节(如果您感兴趣,我使用polyfit())
,只需绘制结果。
这个看起来好多了,至少在我看来是这样。注意现在假设离每个测量值都很近,而在之前的图中假设离测量值很远。我体重增加的可能性比我没有体重增加的可能性要大得多。我真的增重了13磅吗?谁知道呢?这个问题似乎无法回答。
“但这是不可能的吗?”一位同事问道。
让我们来尝试一些疯狂的东西。让我们假设我知道我每天增加一磅。我现在怎么知道并不重要,假设我知道它是近似正确的。也许我每天吃6000卡路里的食物,这会导致体重增加。或者也许有另一种方法来估计体重增加。这是一个思想实验,细节并不重要。让我们看看是否可以利用这些信息,如果有的话。
第一次测量是158。我们没有办法知道有什么不同,所以我们就接受这个估计吧。如果我们今天的体重是158,明天会是多少?好吧,我们认为我们每天体重增加1磅,所以我们的预测是159磅,像这样:
好吧,但这有什么用呢? 当然,我们可以假设每天1磅是准确的,并预测我们未来10天的体重,但如果我们不考虑它的读数,为什么要使用秤呢?我们来看下一个测量。我们再次站到秤上,它显示164.2磅。
现在我们有麻烦了。我们的预测与我们的测量结果不符。但是,这是我们所期望的,对吧?如果预测总是与测量完全相同,它将无法向过滤器添加任何信息。当然,没有理由去测量,因为我们的预测是完美的。
整本书的关键洞见在下一段。仔细阅读!
所以我们要怎么做?如果我们只是从测量来估计的话那么这样预估出来的值是和结果没有关系的。如果我们仅从预判去评估的话,那么又缺少测量的值。只有把预估和测量值两者进行混合才是有效的(我会加深关键字)。
混合两个值,这样的作法看起来特别像前面提到的两个称的问题。使用与之前相同的推理,我们可以看到唯一有意义的是在预测和测量之间选择一个数字。例如,估计值 165 没有意义,157 也没有意义。我们的估计值应介于 159(预测)和 164.2(测量)之间。
再次强调,这是非常重要的。我们认同当出现两个含有误差的值的时候,我们应该形成一个评估值,这个值应该在两个误差值之间。暂时先不关心这些值是如何通过相关的算法产生的。
在这张开始的时候我们有两次测量值,但是现在我们有一个测量值和一个预估值。两种情况下的推理和数学原理都是一样的。我们从来没有丢弃数据信息。我其实是知道的,我看到很多的商业软件直接将噪声数据丢弃。请不要这样做。我们对重量增长的预估可能没有非常的精准,但是我们应该利用好所有的数据信息。
我坚持让你停下来并把这件事请想清楚。我所做的就是用基于人体生理学的不准确体重预测来取代不准确的体重秤。仍是这些数据,数学上的算法不知道这些数据是来自于体重秤还是来自预估。我们有两条带有一定噪声的数据,让后我们想将它们两个进行结合。在这本书的剩下的部分,我们将设计出一些相当复杂的数学算法去实现这些计算,但是这些算法完全不会关心这些数据的来源。这些算法只会基于一些给定的值计算出他们的精确值。
评估值应该介于测量值和预测值之间吗?也许是的,因为在生成这些值的时候我们也许知道我们的预测值是更多还是更少的接近准确值。我们预测的准确性可能与秤测量出来的准确性是不同的。回想一下当 秤A 比 秤B 准确得多时我们所做的事情 - 我们将答案调整为更接近 A 而不是 B。让我们在图表中看一下。
现在让我们尝试随机选择的数字来衡量我们的估计: 4/10 。 我们预估的值的十分之四来自于测量值剩下的十分之六来自于与预测值。通过这些我们表达一些想法,就是相信预测比测量更可能正确。我们计算:
在实际测量值和预测值之间的差被称作残差(residual),由上图中的黑色垂直线表示。后面我们会利用这个非常重要的值,因为它是测量值与滤波器输出之间差异的精确计算。 残差越小意味着性能越好。
让我们对其进行编码,并在针对上面的一系列重量进行测试同时查看一下输出的结果。另外我们还要引入另外一个因子。体重增加的单位为 磅/时间,因此一般来说我们需要添加一个时间增量 t,我们将其设置为 1(磅每天)。
我手工生成了体重数据,对应于 160 磅的真实起始体重和每天 1 磅的体重增加。 换句话说,第一天(第 0 天)的真实重量是 160 磅,第二天(第一天,称重的第一天)真实重量是 161 磅,依此类推。
我们需要猜测初始重量。 现在谈论初始化阶段的策略还为时过早,所以现在我直接假设以 160 磅作为初始值。
太好了,我们现在已经有了很多的数据,现在我们来讨论如何来解释这些数据。蓝色粗线条体现了过滤器的测量值,它开始从第0天我们假设的 160 磅。红色的线条代表预测值,它们是一句前一天的重量得来的,第一天的前一天是 160 磅,重量增长 1 磅,所以第一个预测值是 161 磅。第一天的预估值在测量值于预测值之间,为 159.8 磅。图表下方打印了之前的体重、预测体重和每天的新估计体重。 最后,细黑线显示被称重者的实际体重增加。
每天完成此操作,确保你了解每个步骤的预测和估计是如何形成的。 请注意估计值总是落在测量值和预测值之间。
预估值不是一条直线,但是预估值仍比测量值更直一些,有点接近我们常见的趋势线了。
过滤器的结果可能会让你觉得很愚蠢; 当然,如果我们假设我们的体重增加约为 1 磅/天,那么数据看起来会很好! 让我们看看如果我们最初的猜测是错误的,过滤器会做什么。 我们假设每天减重 1 磅:
这看起来就不是很好了,估计值很快的就和测量值偏离了。显然,要求我们正确猜测变化率的过滤器并不是很有用。即使我们在初始阶段猜测是正确的,过滤器仍会在速率变化后很快发生错误。如果我停止“暴饮暴食”那么过滤器会很难对变化做出调整。请注意,这里是调整,尽管我们告诉它我们每天减重 1 磅,但估计值仍在攀升。 它只是调整得不够快。
但是,‘如果’呢? 如果我们不将体重增加保留在 1 磅(或其他)的初始猜测值,而是根据现有的测量和估计来计算,结果会怎样呢? 第一天我们对重量的估计是:
在第二天我们测量为 164.2,这暗示总量增长了 4.4磅(164.2 - 159.8 = 4.4),而不是 1 磅。我们可以使用某种办法使用这些信息吗?这看起来似乎是很有道理的。毕竟,对重量本身的测量是依据我们现实中测量我们的体重,所以这是很重要的信息。我们对体重增长的估计可能不是完美的,但是肯定比只是猜测我们的体重总是增长1磅要好。数据肯定是好于猜测,哪怕是个噪声。
人们大多在这一点上非常的由于,所以确保你是认同的。两次有噪声的体重测量结果为我们提供了隐含的体重增加/减少。如果测量值是不精确的那么评估值也是非常不精确的,但是这里任然有一些信息在这个计算之中。想象一下用精确到 1 磅的秤给一头牛称重,然后秤显示牛增长了 10 磅。由于误差的存在这个牛可能增长了 8 磅或者 12 磅。 但我们知道它体重增加了,以及大概增加了多少。 这是信息。 我们用信息做什么? 永远不要扔掉它!
再回到我们“暴饮暴食”的例子中。我们应该设置新的每天增加为 4.4 磅?之前我们认为重量应该增加 1 磅,而现在我们认为是 4.4 磅。现在我们有两个数值,而且想要将它们合成。额~看起来相似的问题又来了。让我们使用相同的工具,也是我们迄今为止唯一的工具 - 在两者之间选择一个值。 这次我将使用另一个任意选择的数字,13。 该方程与重量估计相同,只是我们必须考虑时间,因为这是一个比率(增益/天):
我认为这看起来不错,由于最初对体重增加的猜测较差,为 -1,过滤器需要几天时间才能准确预测体重,但一旦做到这一点,它就开始准确跟踪体重。我们没有使用方法去选取系数,4/10 和 1/3(事实上,在这个问题上它们很汗选择)。但除此之外,所有的数学计算都遵循非常合理的假设。回忆一下,你可以改变参数 time_step 到一个很大的值然后重新运行这个代码然后你就看到新的图表了。
在我们继续之前还有最后一点。 在预测步骤中我写了这一行
gain_rate = gain_rate
这个看起来一点用都没有,甚至可以被移除。我写这篇文章是为了强调,在预测步骤中,您需要预测所有变量(权重和增益率)的下一个值。这很快就会变得有意义。 在这种情况下,我们假设增益不变,但当我们推广该算法时,我们将删除该假设。
g-h Filter 上半部分的翻译结束,如果有错误欢迎大家指正!
这本书的名字叫做 《Kalman Filters and Random Signals in Python》 原著的作者 Roger Labbe 是一名在航天航空领域工作近 20 年的资深软件工程师。本书尝试不使用复杂的数学理论去教会让大家了解滤波器的细节。本书中的代码都是使用 Jupyter Notebook 编写的。 使得作者可以更好的将文本、数学、Python 和 Python 输出组合到一处。并且本书有配套的习题和答案,更好的帮你理解滤波器逻辑。
本系列博客可能更新较慢,大家如果感兴趣非常欢迎直接去读原版。
Kalman Filters and Random Signals in Python 原文 Link
因为原文采用了 CC BY-NC-SA 4.0 DEED(署名-非商业性使用-相同方式共享 4.0 国际) 的 license。所以我发翻译也是合乎规则的。
我本身对于上大学学的数学知识也是忘得 7788 了。并且我会以我的角度在未来翻译的时候增加一些名词,数学基础的讲解,这样做的目的也是能将自己学习过程进行记录,如果能有幸帮助被大家读到觉得有帮助的话,还请点赞支持!
CC BY-NC-SA 4.0 DEED License 中文
CC BY-NC-SA 4.0 DEED License 英文