最近看到一本深入浅出的机器学习开源书，是大神 Roger Labbe 的 Kalman and Bayesian Filters in Python
这本书的特点除了写得循序渐进之外，所有的动手环节，都用 python Jupyter Notebook 来实现。另外，大家也可通过 nbviewer 来边阅读边实践体验书籍。
g-h 滤波器
在我们开始之前，请确保你了解如何使用 Jupyter Notebooks，并熟悉 SciPy、NumPy和 Matplotlib 包，因为它们在整本书中都被使用。前言包含对这些包的介绍。
通过思想实验建立直觉
想象一下，我们生活在一个没有秤的世界——你可以站在秤上称重自己的设备。有一天在工作中，一位同事跑到你面前，向你宣布她发明了一个“体重秤”。在她解释之后，你迫不及待地站在上面宣布结果：“172 磅”。你欣喜若狂——这是你一生中第一次知道自己的体重。更重要的是，当你想象将此设备卖给世界各地的减肥诊所时，美元在你的眼中跳动！这是太棒了！
另一位同事听到了骚动，过来想知道是什么让你如此兴奋。你解释了这项发明，然后再次站在秤上，自豪地宣布结果：“161 磅”。然后你犹豫，困惑。
“几秒钟前读数为 172 磅”，你向你的同事抱怨。
“我从未说过它是准确的，”她回答道。
传感器不准确。这就是大量滤波器发明背后的动机，而解决这个问题就是本书的主题。我只能提供过去半个世纪以来发展的解决方案，但这些解决方案是通过对我们所知道的以及我们如何知道的本质提出非常基本的问题而获得。在我们尝试数学之前，让我们跟随那个发现之旅，看看它是否能告诉我们关于滤波器的直觉。
尝试另一种体重秤
有什么方法可以改进这个结果吗？显而易见，首先要尝试的是获得更好的传感器。不幸的是，你的同事告诉你她已经制作了 10 个秤，而且它们的准确度都差不多。你让她拿出另一台秤，你用一台秤称体重，然后用另一台称重。第一个秤 (A) 读数为“160 磅”，第二个秤 (B) 读数为“170 磅”。关于你的体重，我们可以得出什么结论？
那么，我们的选择是什么？
- 我们可以选择只相信 A，并将 160 磅分配给我们的体重估计。
- 我们可以选择只相信 B，并将 170 磅分配给我们的体重。
- 我们可以选择一个小于 A 和 B 的数字。
- 我们可以选择一个大于 A 和 B 的数字。
- 我们可以在 A 和 B 之间选择一个数字。
前两个选择是合理的，但我们没有理由偏向于其中一种。为什么我们会选择相信 A 而不是 B？我们没有理由这样相信。第三个和第四个选择是不合理的。诚然，这些比例不是很准确，但根本没有理由选择一个超出他们测量范围的数字。最后的选择是唯一合理的。如果两种秤都不准确，并且可能给出高于或低于我的实际体重的结果，则答案通常介于 A 和 B 之间。
在数学中，这个概念被形式化为期望值，稍后我们将深入介绍它。现在问问自己，如果我们读取一百万次读数，“通常”会发生什么。有时两个秤的读数都太低，有时都太高，其余时间它们会覆盖实际重量。如果它们覆盖实际重量，那么我们当然应该在 A 和 B 之间选择一个数字。如果它们没有覆盖，那么我们不知道它们是太高还是太低，但是通过选择 A 和 B 之间的一个数字，我们可以至少减轻最差测量的影响。例如，假设我们的实际体重是 180 磅。160 磅是一个很大的错误。但如果我们选择 160 磅到 170 磅之间的重量，我们的估计值将优于 160 磅。如果两个秤返回的值都大于实际重量，则同样的论点成立。
我们稍后会更正式地处理这个问题，但现在我希望很清楚我们的最佳估计是 A 和 B 的平均值。
我们可以图形化地看一下。我绘制了 A 和 B 的测量值，假设误差为±8 磅。测量值介于 160 到 170 之间，因此唯一有意义的重量必须在 160 到 170 磅之间。
In [3]:
import kf_book.book_plots as book_plots
from kf_book.book_plots import plot_errorbars
plot_errorbars([(160, 8, 'A'), (170, 8, 'B')], xlims=(150, 180))
所以 165 磅看起来是一个合理的估计，但这里有更多信息我们或许可以利用。唯一可能的重量位于 A 和 B 的误差线之间的交叉点。例如，161 磅的重量是不可能的，因为秤 B 无法给出 170 磅的读数，最大误差为 8 磅。同样，169 磅的重量也是不可能的，因为秤 A 不能给出 160 磅的读数，最大误差为 8 磅。在这个例子中，唯一可能的重量在 162 到 168 磅的范围内。
这还不能让我们找到更好的重量估计，但让我们再开启“假设”。如果我们现在被告知 A 比 B 准确三倍怎么办？考虑我们上面列出的 5 个选项。选择 A 和 B 范围之外的数字仍然没有意义，因此我们不会考虑这些。选择 A 作为我们的估计似乎更有说服力——毕竟，我们知道它更准确，为什么不用它代替 B？B 是否可以单独提高我们对 A 的了解？
答案可能与直觉相反：是的，它可以。首先，让我们看一下 A=160 和 B=170 的相同测量值，但误差为 A ±3 磅 B 的误差是它的 3 倍，±9 磅。
In [4]:
plot_errorbars([(160, 3, 'A'), (170, 9, 'B')], xlims=(150, 180))
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。我们应该估计自己的体重是多少？让我们以图形方式看一下。
In [5]:
plot_errorbars([(160, 1, 'A'), (170, 9, 'B')], xlims=(150, 180))
在这里我们可以看到唯一可能的重量是 161 磅。这是一个重要的结果。使用两个相对不准确的传感器，我们能够推断出极其准确的结果。
所以两个传感器，即使一个不如另一个准确，也比一个好。我将在本书的其余部分反复强调这一点。我们从不丢弃信息，无论信息多么糟糕。我们将开发数学和算法，使我们能够包含所有可能的信息来源，以形成可能的最佳估计。
然而，我们偏离了我们的问题。没有客户会想要购买多个秤，此外，我们最初假设所有秤都同样（不）准确。这种不考虑准确性而使用所有测量值的见解将在以后发挥重要作用，所以不要忘记它。
如果我有一个体重秤，但我称了很多次体重怎么办？我们得出结论，如果我们有两个精度相同的秤，我们应该对它们的测量结果进行平均。如果我用一台秤称自己 10,000 次会怎样？我们已经说过，秤返回一个太大的数字和返回一个太小的数字的可能性是一样的。证明大量权重的平均值会非常接近实际权重这并不难，让我们写一个模拟来确认。我将使用 NumPy，它是数值计算 SciPy 生态系统的一部分。
In [6]:
import numpy as np
measurements = np.random.uniform(160, 170, size=10000)
mean = measurements.mean()
print(f'Average of measurements is {mean:.4f}')
Average of measurements is 164.9353
打印出确切数字取决于你的随机数生成器，但它应该非常接近 165。
这段代码做出了一个可能不正确的假设——对于 165 磅的真实重量，秤的读数可能为 160 和 165。这几乎从来都不是真的。真正的传感器更有可能获得更接近真实值的读数，并且越来越不可能获得距离真实值越远的读数。我们将在高斯章节中详细介绍这一点。现在，我将不做进一步解释地使用该numpy.random.normal()函数，它将在接近 165 磅的地方产生更多的值，而在更远的地方产生更少的值。现在请相信，这将产生类似于真实秤工作方式的有噪声测量结果。
In [7]:
mean = np.random.normal(165, 5, size=10000).mean()
print(f'Average of measurements is {mean:.4f}')
Average of measurements is 164.9762
同样，答案非常接近 165。
太好了，我们的传感器问题有了答案！但这不是一个非常实际的答案。没有人有耐心给自己称一万次，或者十几次。
那么，让我们开启“假设”。如果你每天测量一次体重，读数分别为 170、161 和 169，结果会怎样？你的体重是增加了还是减轻了，或者这只是有噪声的测量结果？
我们真的不能说。第一次测量是 170，最后一次是 169，意味着减重 1 磅。但如果秤只精确到 10 磅，那可以用噪声来解释。我本来可以增重的；也许我第一天的体重是 165 磅，第三天是 172 磅。可以通过体重增加获得这些体重读数。我的体重秤告诉我我正在减肥，实际上我正在增加体重！让我们在图表中看一下。我已经绘制了测量值和误差条，然后绘制了一些可能的体重增加/减少，这些增加/减少可以用绿色虚线中的这些测量值来解释。
In [8]:
mean = np.random.normal(165, 5, size=10000).mean()
print(f'Average of measurements is {mean:.4f}')
正如我们所见，这三个测量值可以解释体重变化的极端范围。事实上，有无数种选择。我们要放弃吗？至少我不会！回想一下，我们正在谈论测量一个人的体重。没有合理的方法可以让一个人在第 1 天体重达到 180 磅，在第 3 天体重达到 160 磅，或者在一天内减掉 30 磅却在第二天又恢复回来（我们假设此人没有截肢或其他创伤发生）。
我们正在测量的物理系统的行为应该影响我们如何解释测量结果。如果我们每天称重一块石头，我们会将所有差异归因于噪声。如果我们正在称量一个用于做家务的雨水蓄水池，我们可能会相信这种重量变化是真实的。
假设我使用不同的秤，得到以下测量值：169、170、169、171、170、171、169、170、169、170。你的直觉告诉你什么？例如，你可能每天增加 1 磅，而有噪声的测量结果恰好看起来你保持相同的体重。同样，你可以每天减掉 1 磅并获得相同的读数。但这可能吗？抛硬币连续出现 10 次正面朝上的可能性有多大？不太可能。我们不能仅根据这些读数来证明这一点，但我的体重似乎很可能保持稳定。在下面的图表中，我用误差条绘制了测量值，并用绿色虚线标出了可能的真实重量。这条虚线并不意味着是这个问题的“正确”答案，它只是一个合理的并且可以通过测量来解释的答案。
In [9]:
gh.plot_hypothesis2()
另一个假设：如果读数是 158.0、164.2、160.3、159.9、162.1、164.6、169.6、167.4、166.4、171.0 怎么办？让我们看一下图表，再回答问题。
In [10]:
gh.plot_hypothesis3()
我是否“看起来”减肥了，而这些只是噪声非常大的数据？并不太可能。或者看起来我有可能保持相同的体重吗？也并不太可能。该数据随时间呈上升趋势；不明显，但绝对是向上的趋势。我们不能确定，但这看起来像是体重增加，而且是显着的体重增加。让我们用更多的图来测试这个假设。与表格相比，图表中的数据通常更容易一眼看出来。
那么让我们来看看两个假设。首先，假设我们的体重没有改变。为了得到这个数字，我们认同应该平均测量值。。
In [11]:
gh.plot_hypothesis4()
这看起来不太有说服力。事实上，我们可以看到并不能在所有误差线内绘制一条水平线。
现在，假设我们体重增加了。多少？我不知道，但 NumPy 知道！我们想通过看起来“大约”正确的测量值画一条线。NumPy 具有称为最小二乘拟合的函数。我们不用担心该绘制函数的细节（如果你有兴趣，我会使用 polyfit()），只需绘制结果即可。
In [12]:
gh.plot_hypothesis5()
这看起来好多了，至少在我看来是这样的。现在请注意，假设的线离每个测量值都非常近，而在之前的图中，假设的线离很多测量值都很远。我体重增加的可能性似乎比我没有体重增加的可能性大得多。我真的增加了 13 磅吗？谁能说？这似乎无法回答。
“但这不可能吗？” 大家也许会问。
让我们尝试一些疯狂的事情。让我们假设我知道我每天增加大约一磅。现在我是怎么知道的并不重要，假设我知道它是大致正确的。也许我每天进食 6000 卡路里，这会导致体重增加。或者也许还有另一种方法来估计体重增加。这是一个思想实验，细节并不重要。让我们看看我们是否可以利用此类信息（如果可用的话）。
第一次测量是 158。我们无法知道有什么不同，所以让我们接受它作为我们的估计。如果我们今天的体重是 158，明天会是多少？好吧，我们认为我们的体重增加了 1 磅/天，所以我们的预测是 159，如下所示：
In [13]:
gh.plot_estimate_chart_1()
好吧，但这有什么好处呢？当然，我们可以假设 1 磅/天是准确的，并预测接下来 10 天的体重，但如果我们不合并其读数，为什么还要使用体重秤呢？那么让我们看看下一个测量。我们再次踏上体重秤，它显示 164.2 磅。
In [14]:
gh.plot_estimate_chart_2()
现在出现问题了。我们的预测与测量不符。但是，这就是我们所期望的，对吧？如果预测总是与测量完全相同，那么它就无法向滤波器添加任何信息。当然，也就没有必要去测量，因为我们的预测是完美的。
整本书的关键见解在下一段。仔细阅读！
那么我们该怎么办？如果我们仅从测量中形成估计，那么预测将不会影响结果。如果我们只根据预测形成估计，那么测量将被忽略。如果这是可行的，我们需要将预测和测量进行某种混合（我已将关键点加粗）。
混合两个值——这听起来很像之前的两个秤的问题。使用与之前相同的推理，我们可以看到唯一有意义的是在预测和测量之间选择一个数字。例如，165 的估计值没有意义，157 也没有意义。我们的估计值应该介于 159（预测值）和 164.2（测量值）之间。
再说一遍，这里太重要了。我们一致认为，当出现两个有误差的值时，我们应该在这两个值之间形成一个估计值。这些值是如何生成的并不重要。在本章开头我们有两个测量值，但现在我们有一个测量值和一个预测值。在这两种情况下，推理和数学都是相同的。我们从不丢弃信息。我是认真的。我看到很多商业软件会丢弃噪声数据。不要这样做！我们对体重增加的预测可能不是很准确，但只要有一些信息，我们就应该使用它。
我必须坚持让你停下来认真考虑一下。我所做的只是用基于人体生理学的不准确体重预测代替不准确的体重秤。它仍然是数据。数学不知道数据是来自测量还是预测。我们有两份带有一定噪声的数据，我们想把它们结合起来。在本书的其余部分，我们将开发一些相当复杂的数学来执行此计算，但数学从不关心数据来自何处，它只根据这些值和其准确性进行计算。
估计值应该介于测量值和预测值之间吗？也许吧，但总的来说，我们似乎知道我们的预测与测量结果相比或多或少是准确的。我们预测的准确性可能与测量的准确性不同。回想一下我们在体重秤 A 比体重秤 B 准确得多时所做的事情——我们将答案调整为更接近 A 而不是 B。让我们在图表中看一下。
In [15]:
gh.plot_estimate_chart_3()
现在让我们尝试随机选择一个数字来缩放我们的估计：4/10. 我们的估计将是测量值的十分之四，其余部分将来自预测。换句话说，我们在这里表达了一种信念，一种相信预测比测量更可能正确的信念。我们计算为
测量值和预测值之间的差异称为残差，在上图中用黑色垂直线表示。这将成为稍后使用的重要值，因为它是测量值与滤波器输出之间差异的精确计算。较小的残差意味着更好的性能。
让我们编写代码，并在我们根据上面的一系列权重对其进行测试时查看结果。我们必须考虑另一个因素。体重增加的单位是磅/时间，所以一般来说我们需要添加一个时间步长 ，我们将其设置为 1（天）。
我手工生成的体重数据对应于 160 磅的真实起始体重，以及每天增加 1 磅的体重。换句话说，第一天（第零天）真实体重是 160 磅，第二天（第一天，称重的第一天）真实体重是 161 磅，依此类推。
我们需要对初始权重进行猜测。现在谈论初始化策略还为时过早，所以现在我假设 160 磅。
In [16]:
from kf_book.book_plots import figsize
import matplotlib.pyplot as plt
weights = [158.0, 164.2, 160.3, 159.9, 162.1, 164.6,
169.6, 167.4, 166.4, 171.0, 171.2, 172.6]
time_step = 1.0 # day
scale_factor = 4.0/10
def predict_using_gain_guess(estimated_weight, gain_rate, do_print=False):
# storage for the filtered results
estimates, predictions = [estimated_weight], []
# most filter literature uses 'z' for measurements
for z in weights:
# predict new position
predicted_weight = estimated_weight + gain_rate * time_step
# update filter
estimated_weight = predicted_weight + scale_factor * (z - predicted_weight)
# save and log
estimates.append(estimated_weight)
predictions.append(predicted_weight)
if do_print:
gh.print_results(estimates, predicted_weight, estimated_weight)
return estimates, predictions
initial_estimate = 160.
estimates, predictions = predict_using_gain_guess(
estimated_weight=initial_estimate, gain_rate=1, do_print=True)
previous estimate: 160.00, prediction: 161.00, estimate 159.80
previous estimate: 159.80, prediction: 160.80, estimate 162.16
previous estimate: 162.16, prediction: 163.16, estimate 162.02
previous estimate: 162.02, prediction: 163.02, estimate 161.77
previous estimate: 161.77, prediction: 162.77, estimate 162.50
previous estimate: 162.50, prediction: 163.50, estimate 163.94
previous estimate: 163.94, prediction: 164.94, estimate 166.80
previous estimate: 166.80, prediction: 167.80, estimate 167.64
previous estimate: 167.64, prediction: 168.64, estimate 167.75
previous estimate: 167.75, prediction: 168.75, estimate 169.65
previous estimate: 169.65, prediction: 170.65, estimate 170.87
previous estimate: 170.87, prediction: 171.87, estimate 172.16
In [17]:
# plot results
book_plots.set_figsize(10)
gh.plot_gh_results(weights, estimates, predictions, [160, 172])
weights
Out[17]:
[158.0,
164.2,
160.3,
159.9,
162.1,
164.6,
169.6,
167.4,
166.4,
171.0,
171.2,
172.6]
结果不错！图中数据很多，下面我们就来说说如何解读。蓝色粗线显示滤波器的估计值。它从第 0 天开始，初始猜测为 160 磅。红线显示根据前一天的体重做出的预测。所以，在第一天，之前的体重是 160 磅，体重增加了 1 磅，所以第一个预测是 161 磅。第一天的估计值介于预测值和测量值之间，为 159.8 磅。图表下方是打印出的前一天体重、预测体重和每天的新估计体重。最后，细黑线表示被称重者的实际体重增加。
观察一遍每天的点，确保你了解每一步的预测和估计是如何形成的。请注意估计值如何始终介于测量值和预测值之间。
估计值不是一条直线，但它们比测量值更直，并且有点接近我们创建的趋势线。此外，它似乎随着时间的推移变得更好。
滤波器的结果可能会让你觉得很愚蠢；当然，如果我们假设结论是我们的体重增加大约 1 磅/天，那么数据看起来会很好！如果我们的初始猜测是错误的，让我们看看滤波器会做什么。让我们预测每天体重减轻 1 磅：
In [18]:
e, p = predict_using_gain_guess(initial_estimate, -1.)
gh.plot_gh_results(weights, e, p, [160, 172])
那不是那么令人印象深刻。估计很快就偏离了测量结果。显然，要求我们正确猜测变化率的滤波器不是很有用。即使我们最初的猜测是正确的，只要变化率发生变化，滤波器就会失效。如果我停止暴饮暴食，滤波器将很难适应这种变化。注意是调整！尽管我们告诉它我们每天减掉 1 磅，但估计值仍在攀升。它只是不能足够快地调整。
但是，“如果”呢？如果我们不是将体重增加保持在 1 磅（或其他）的初始猜测值，而是根据现有的测量值和估计值进行计算，会怎样呢？在第一天，我们对体重的估计是：
第二天，我们测出 164.2，这意味着体重增加了 4.4 磅（因为 164.2 - 159.8 = 4.4），而不是 1。我们能否以某种方式使用此信息？这似乎是合理的。毕竟，体重测量本身是基于我们体重的真实世界测量，因此有有用的信息。我们对体重增加的估计可能并不完美，但肯定比仅仅猜测我们的体重增加了 1 磅要好。数据比猜测更好，即使它有噪声。
人们在这一点上真的很犹豫，所以要确保你同意。两次有噪声体重测量结果为我们提供了隐含的体重增加/减少。如果测量不准确，该估计将非常不准确，但此计算中仍有信息。想象一下，用精确到 1 磅的秤称一头牛，它显示这头牛增加了 10 磅。根据误差，这头奶牛可能增加了 8 磅到 12 磅，但我们知道它增加了重量，并且大致增加了多少。这是信息。永远不要扔掉它！
回到我的体重。我们应该将新的增重/天设置为 4.4 磅吗？昨天我们认为体重增加了 1 磅，今天我们认为是 4.4 磅。我们有两个数字，想以某种方式将它们组合起来。嗯，听起来又是我们同样的问题。让我们使用相同的工具，也是迄今为止唯一的工具——在两者之间选择一个值。这次我将使用另一个任意选择的数字，。该等式与体重估计相同，只是我们必须考虑时间，因为这是一个比率（增重/天）：
In [19]:
weight = 160. # initial guess
gain_rate = -1.0 # initial guess
time_step = 1.
weight_scale = 4./10
gain_scale = 1./3
estimates = [weight]
predictions = []
for z in weights:
# prediction step
weight = weight + gain_rate*time_step
gain_rate = gain_rate
predictions.append(weight)
# update step
residual = z - weight
gain_rate = gain_rate + gain_scale * (residual/time_step)
weight = weight + weight_scale * residual
estimates.append(weight)
gh.plot_gh_results(weights, estimates, predictions, [160, 172])
我认为这开始看起来非常好。由于最初对体重增加的猜测很差 -1，滤波器需要几天时间才能准确预测体重，但一旦做到这一点，它就会开始准确跟踪体重。我们没有使用任何方法来选择我们的比例因子 4/10 和 1/3（实际上，对于这个问题，它们是糟糕的选择），但除此之外，所有数学都是从非常合理的假设得出的。回想一下，你可以将参数的值更改为time_step更大的值，然后重新运行单元格，如果你希望看到逐步绘制的绘图。
在我们继续之前的最后一点。在预测步骤中，我写了这行
gain_rate = gain_rate这显然没有效果，可以去掉。我写这篇文章是为了强调在预测步骤中，你需要预测所有变量的下一个值，包括weight和gain_rate。这将很快变得相关。在这种情况下，我们假设增量不变，但当我们推广此算法时，我们将去除该假设。
如果大家觉得不错并且想看到下一期连载，请一键三连哦。更多AI干货，欢迎关注 MyEncyclopedia