【经典教程翻译】卡尔曼与贝叶斯滤波器:离散贝叶斯滤波器(上)
本期继续大神 Roger Labbe 的 Kalman and Bayesian Filters in Python
。本文为第二章离散贝叶斯过滤器,通过逐渐复杂的例子建立离散贝叶斯的直觉观念。
所有文章首发于 MyEncyclopedia公众号
【经典教程翻译】卡尔曼与贝叶斯滤波器:直觉理解滤波器背后的原理(上)
https://mp.weixin.qq.com/s/TTkmCgYimWhcJeMxBh05fg
【经典教程翻译】卡尔曼与贝叶斯滤波器:直觉理解滤波器背后的原理(下:滤波器的思考框架)
https://mp.weixin.qq.com/s/WObX_20GDr1Kc5vmWtae4Q
离散贝叶斯滤波器
卡尔曼滤波器属于称为贝叶斯滤波器
的滤波器家族。大多数卡尔曼滤波器的教科书介绍了贝叶斯公式,也许显示了它如何影响卡尔曼滤波器方程,但大多数情况下讨论的非常抽象。
这种方法需要对数学的几个领域有相当复杂的理解,尽管它仍然把理解和形成直觉的大部分工作留给了读者。
我将使用不同的方式来展开这个主题,我对 Dieter Fox 和 Sebastian Thrun 的工作感激不尽。他们通过跟踪穿过走廊的物体来建立贝叶斯统计如何工作的直觉——他们使用机器人,我使用狗。我喜欢狗,它们比机器人更难以预测,这给滤波带来了有趣的困难。我能找到的第一个发布的例子似乎是 Fox 1999 (Monte carlo localization: Efficient position estimation for mobile robots),Fox 2003 (Bayesian Filters for Location Estimation) 中有一个更完整的例子。Sebastian Thrun 在他出色的 Udacity 课程 Artificial Intelligence for Robotics
中也使用了这个公式。事实上,如果您喜欢看视频,我强烈建议您暂停阅读本书,转而阅读该课程的前几课,然后再返回本书以更深入地探讨该主题。
现在让我们使用一个简单的思想实验,就像我们对 g-h 滤波器所做的一样,看看我们如何推断使用概率进行滤波和跟踪。
追踪小狗
让我们从一个简单的问题开始。我们有一个欢迎狗的工作空间,所以人们带着他们的狗来工作。偶尔,这些狗会走出办公室,穿过大厅。我们希望能够追踪他们。因此,在黑客马拉松期间,有人发明了一种声纳传感器,可以挂在狗的项圈上。它发出信号,聆听回声,根据回声返回的速度,我们可以判断狗是否在敞开的门口前。它还会感知狗何时行走,并报告狗移动的方向。它通过 wifi 连接到网络并每秒发送一次更新。
我想追踪我的小狗 Simon ,所以我将设备系在它的项圈上,然后启动 Python,准备编写代码来追踪它穿过大楼。乍一看,这似乎是不可能的。如果我开始听 Simon 衣领上的传感器,我可能会读到door、hall、hall等等。我如何使用该信息确定 Simon 的位置?
为了使问题足够小以便于绘制,我们假设走廊中只有 10 个位置,我们将其编号为 0 到 9,其中 1 在 0 的右边。出于稍后将清楚的原因,我们还将假设走廊是圆形或矩形的。如果你从位置 9 向右移动,你将在位置 0。
当我开始听传感器时,我没有理由相信Simon在走廊的任何特定位置。从我的角度来看,他处于任何位置的可能性都是一样的。有 10 个位置,所以他在任何给定位置的概率是 1/10。
让我们在 NumPy 数组中表示我们对他的位置的信念。我可以使用 Python List,但 NumPy 数组提供了我们即将使用的功能。
In [3]:
在贝叶斯统计
中,这称为先验
。它是合并测量或其他信息之前的概率。更完整地说,这称为先验概率分布
。概率分布
是一个事件所有可能概率的集合。概率分布总和为 1,因为某些事情必须发生;该分布列出了所有可能的事件和每个事件的概率。
我相信您以前使用过概率——例如“今天下雨的概率是 30%”。最后一段听起来更多。但是贝叶斯统计是概率的一场革命,因为它将概率视为对单个事件的信念。让我们举个例子。我知道如果我无限次地掷一枚公平的硬币,我将得到 50% 的正面和 50% 的反面。这被称为频率派统计
,以区别于贝叶斯统计
,来计算基于事件发生的频率。
我又掷了一次硬币,让它落地。我认为它以哪种方式着陆?频率派概率论对此无话可说。它只会说 50% 的硬币以正面朝上。在某些方面,将概率分配给硬币的当前状态是没有意义的。它要么是正面,要么是反面,我们只是不知道是哪一个。贝叶斯将其视为对单个事件的信念——我相信或知道这个特定的抛硬币是正面的强度是 50%。有些人反对“信念”一词;信念可能意味着在没有证据的情况下认为某事是真实的。在本书中,它始终是我们知识强度的衡量标准。随着我们的进行,我们将了解更多相关信息。
贝叶斯统计考虑了过去的信息(先验信息)。我们观察到每 100 天下 4 次雨。由此我可以说明天下雨的可能性是 1/25。这不是天气预报的方式。如果我知道今天在下雨并且风暴前线停滞不前,那么明天很可能会下雨。天气预报是贝叶斯的。
在实践中,统计学家混合使用频率派和贝叶斯技术。有时很难或不可能找到先验,因而频率派占主导地位。在本书中我们可以找到先验。当我谈论某事的概率时,我指的是在给定过去事件的情况下某些特定事物为真的概率。当我这样做时,我采用了贝叶斯方法。
现在让我们创建走廊的地图。我们将把前两扇门靠在一起,然后将另一扇门放在更远的地方。我们将 1 用于门,0 用于墙:
In [4]:
我开始在网络上监听 Simon 的传输,我从传感器获得的第一个数据是door。目前假设传感器总是返回正确的答案。由此我断定他在一扇门前,但哪一扇呢?我没有理由相信他在第一、第二或第三扇门前。我能做的是为每扇门分配一个概率。所有门的可能性都一样,而且一共有三扇门,所以我为每扇门分配了 1/3 的概率。
In [5]:

这种分布称为分类分布
,它是描述观察概率的离散分布结果。这是一个
多峰分布
,因为我们对我们的狗的位置有多种看法。当然,我们并不是说我们认为他同时在三个不同的位置,只是我们将我们的知识范围缩小到这三个位置之一。我的(贝叶斯)信念是,有 33.3% 的机会在 0 号门,33.3% 的机会在 1 号门,33.3% 的机会在 8 号门。
这是两个方面的改进。我已经拒绝了一些不可能的走廊位置,我对剩余位置的信念强度从 10% 增加到 33%。随着我们知识的提高,概率将接近 100%。
关于分布的众数
的几句话。给定一个数字列表,例如 {1, 2, 2, 2, 3, 3, 4},众数
是最常出现的数字。对于这个集合,众数是 2。一个分布可以包含多个众数。列表 {1, 2, 2, 2, 3, 3, 4, 4, 4} 的众数是 2 和 4,因为它们都出现了三次。我们说前者是单峰
的,后者是多峰
的。
用于此分布的另一个术语是直方图
。直方图以图形方式描述了一组数字的分布。上面的条形图是一个直方图。
我在上面的代码中手工编码了belief
数组。我们将如何在代码中实现它?我们用 1 表示门,用 0 表示墙,所以我们将走廊变量乘以百分比,如下所示;
In [6]:
让我们把 Python 放在一边,稍微思考一下这个问题。假设我们要从 Simon 的传感器中读取以下内容:
门
向右移
门
我们能推断出Simon的位置吗?当然!考虑到走廊的布局,只有一个地方可以得到这个序列,那就是在左端。因此我们可以自信地说Simon在第二个门口。如果你不清楚,假设Simon是从第二或第三扇门开始的。向右移动后,他的传感器会返回“墙”。这与传感器读数不符,所以我们知道他不是从那里开始的。对于所有剩余的起始位置,我们可以继续使用该逻辑。唯一的可能,就是他现在就在第二道门前。我们的信念是:
In [7]:
我设计了走廊布局和传感器读数,以便快速为我们提供准确的答案。真正的问题并不是那么明确。但这应该触发你的直觉——第一个传感器读数只给了我们Simon位置的低概率 (0.333),但在位置更新和另一个传感器读数之后,我们更了解他在哪里。你可能会怀疑,如果你有一个很长的走廊,有很多门,那么在几次传感器读数和位置更新之后,我们要么能够知道Simon在哪里,要么将可能性缩小到少数可能性。当一组传感器读数仅匹配一个或几个起始位置时,这是可能的。
我们现在可以实施这个解决方案,但让我们考虑这个问题在现实世界中的复杂性。
有噪声的传感器
完美的传感器很少见。如果Simon坐在门前抓挠自己,传感器可能不会检测到门,或者如果他没有面朝走廊,传感器就会误读。因此,当我得到 door ,就不能使用 1/3 作为概率。我要给每扇门分配不到1/3,给每个空白墙位置分配一个小概率。就像是
乍一看,这似乎是无法克服的。如果传感器有噪音,就会对每条数据产生怀疑。如果我们总是不确定,我们怎么能得出结论呢?
对于上述问题,答案是有概率的。我们已经很乐意为狗的位置分配概率信念;现在我们必须考虑由传感器噪声引起的额外不确定性。
假设我们得到了door的读数,并假设测试表明传感器正确的可能性是错误的 3 倍。我们应该在有门的地方将概率分布缩放 3 倍。如果我们这样做,结果将不再是概率分布,但我们稍后会学习如何解决这个问题。
让我们在 Python 代码中看一下。这里我使用变量z
来表示测量。z
或者y
是测量文献中的习惯选择。作为一名程序员,我更喜欢有意义的变量名,但我希望您能够阅读文献或其他滤波代码,所以我现在将开始介绍这些缩写名称。
In [8]:

这不是概率分布,因为它的总和不等于 1.0。但代码基本上做对了——门被分配的数字 (0.3) 比墙 (0.1) 高 3 倍。我们需要做的就是对结果进行归一化,以便概率总和正确为 1.0。归一化是通过将每个元素除以列表中所有元素的总和来完成的。使用 NumPy 很容易:
Out[9]:
FilterPy 使用以下函数实现此normalize
功能:
说“正确的可能性是错误的 3 倍”有点奇怪。我们正在研究概率,所以让我们指定传感器正确的概率,并从中计算比例因子。方程式是
此外,for
循环很麻烦。作为一般规则,您将希望避免for
在 NumPy 代码中使用循环。NumPy 是在 C 和 Fortran 中实现的,因此如果您避免 for 循环,结果通常比等效循环运行速度快 100 倍。
我们如何摆脱这个for
循环?NumPy 允许您使用布尔数组索引数组。您使用逻辑运算符创建一个布尔数组。我们可以找到走廊中的所有门:
Out[10]:
当您使用布尔数组作为另一个数组的索引时,它只返回索引所在的元素True
。因此我们可以如此代替 for
教你 NumPy 超出了本书的范围。我将使用惯用的 NumPy 结构,并在我第一次展示它们时对其进行解释。如果您是 NumPy 的新手,有很多博客文章和视频介绍如何高效且地道地使用 NumPy。
这是我们的改进版本:
In [11]:

我们可以从输出中看到,总和现在是 1.0,门与墙的概率仍然是原来的三倍。结果也符合我们的直觉,门的概率必须小于 0.333,而墙的概率必须大于 0.0。最后,它应该符合我们的直觉,即我们还没有得到任何信息可以让我们区分任何给定的门或墙的位置,所以所有的门位置应该具有相同的值,墙的位置也应该如此。
这个结果称为后验
,是后验概率分布
的简称。它们都表示合并测量信息后的概率分布(在此上下文中后验表示“之后”)。回顾一下,先验
是包含测量信息之前的概率分布。
另一个术语是似然
。当我们计算时,belief[hall==z] *= scale
我们正在计算每个位置被测量的可能性。似然值不是概率分布,因为它总和不为一。
这些的组合给出了等式
当我们谈论滤波器的输出时,我们通常将执行预测后的状态称为先验
或预测
,并将更新后的状态称为后验
或状态估计
。
学习和内化这些术语非常重要,因为大多数文献都广泛使用它们。
scaled_update()
执行了此计算?确实如此。让我把它改写成这种形式:
In [12]:
这个功能并不完全通用。它包含有关走廊的知识,以及我们如何将测量值与其匹配。我们总是努力编写通用函数。这里我们将从函数中移除似然计算,并要求调用者自己计算似然。
这是该算法的完整实现:
似然的计算因问题而异。例如,传感器可能不会只返回 1 或 0,而是返回float
0 和 1 之间的值,表示在门前的概率。它可能会使用计算机视觉并报告一个斑点形状,然后您可以概率地将其与门匹配。它可能会使用声纳并返回距离读数。在每种情况下,似然的计算都会不同。我们将在本书中看到许多这样的例子,并学习如何进行这些计算。
FilterPy 实现update
. 这是完全通用形式的先前示例:
In [13]:
Out[13]:
结合运动测量值
回想一下,当我们结合一系列测量和运动更新时,我们能够以多快的速度找到一个精确的解决方案。然而,这发生在一个完美传感器的虚构世界中。我们能否找到带有噪声传感器的精确解决方案?
不幸的是,答案是否定的。即使传感器读数与极其复杂的走廊地图完美匹配,我们也不能 100% 确定狗在特定位置,毕竟,每个传感器读数都有错误的可能性很小!自然地,在更典型的情况下,大多数传感器读数都是正确的,我们可能接近 100% 确定我们的答案,但永远不会 100% 确定。这可能看起来很复杂,但让我们继续进行数学编程。
首先让我们处理简单的情况——假设运动传感器是完美的,它报告说狗向右移动了一个空间。我们将如何改变我们的belief
阵列?
我希望经过片刻的思考后,很明显我们应该将所有值向右移动一个空格。如果我们之前认为 Simon 有 50% 的机会在位置 3,那么在他向右移动一个位置后我们应该相信他有 50% 的机会在位置 4。走廊是圆形的,所以我们将使用模运算来执行移位。
In [14]:

我们可以看到我们正确地将所有值向右移动了一个位置,从数组的末尾回到开头。
下一个单元格对此进行动画处理,以便您可以看到它的运行情况。使用滑块及时向前和向后移动。这模拟了Simon在走廊里走来走去。它尚未包含新的测量值,因此概率分布不会改变形状,只会改变位置。
In [15]:
复习术语
让我们暂停一下来回顾一下术语。我在上一章介绍了这个术语,但让我们花点时间来帮助巩固你的知识。
系统
是我们试图建模或滤波的对象。这里的系统是我们的狗。状态
是其当前配置或值。在本章中,状态是我们狗的位置。我们很少知道实际状态,所以我们说我们的滤波器产生系统的状态估计
。在实践中,这通常被称为状态
,所以要小心理解上下文。
预测和更新测量的一个周期称为状态的或系统的迭代
,它是时间迭代
的简称。另一个术语是系统传播
。它指的是系统的状态如何随时间变化。对于滤波器,时间通常是离散的步长,例如 1 秒。对于我们的狗追踪器,系统状态是狗的位置,状态迭代
是经过离散时间后的位置。
我们使用过程模型
对系统行为进行建模。在这里,我们的过程模型是狗在每个时间步移动一个或多个位置。这不是一个关于狗行为的特别准确的模型。模型中的错误称为系统错误
或过程错误
。
预测是我们的新先验
。时间向前发展,我们在不知道测量结果的情况下做出了预测。
让我们举个例子。狗的当前位置是 17 m。我们的迭代时长为 2 秒,狗以 15 m/s 的速度行进。我们预测他将在两秒钟内到达哪里?
显然,
我在变量上使用条形图来表示它们是先验(预测)。我们可以这样写过程模型的等式:
是当前位置或状态。如果狗在 17 m 处,则
.
是 x 的状态传播函数。它描述了
在一个时间步内发生变化多少。对于我们的示例,它执行计算
,所以我们将其定义为
.
给预测增加不确定性
perfect_predict()
假设完美的测量,但所有传感器都有噪音。如果传感器报告我们的狗移动了一个空间,但他实际上移动了两个空间,或者零个空间怎么办?这听起来像是一个无法克服的问题,但让我们对其建模并看看会发生什么。
假设传感器的移动测量有 80% 的可能性是正确的,有 10% 的可能性向右超过一个位置,有 10% 的可能性低于向左的位置。也就是说,如果移动测量值为 4(即向右移动 4 个空格),则狗向右移动 4 个空格的可能性为 80%,向右移动 3 个空格的可能性为 10%,向右移动 5 个空格的可能性为 10%。
数组中的每个结果现在都需要包含 3 种不同情况的概率。例如,考虑测量值 2。如果我们 100% 确定狗从位置 3 开始,那么他有 80% 的机会在 5,10% 的机会在 4 或 6。让我们尝试编码:
In [16]:

它似乎工作正常。现在,当我们的信念不是 100% 确定时会发生什么?
In [17]:
Out[17]:

这里的结果比较复杂,但你应该仍然能够在脑海中计算出来。0.04 (新位置3)是由于 0.4 (原位置2)有10%过低测量导致的。0.38 (新位置4)是由于以下原因:80% 的机会我们从原位置2移动了 2 个位置(0.4$\times\times$0.1)。过高的测量值在这里不起作用,因为如果我们从 0.4 和 0.6 对应的位置得到过高测量值都会超过这个位置。我强烈建议做一些例子,直到所有这些都非常清楚,因为接下来的大部分内容都取决于对这一步的理解。
如果您在执行更新后查看概率,您可能会感到沮丧。在上面的示例中,我们从两个位置的概率 0.4 和 0.6 开始;执行更新后,概率不仅降低了,而且分散在地图上。
这不是巧合,也不是精心挑选的例子的结果——预测总是正确的。如果传感器有噪声,我们会丢失一些关于每个预测的信息。假设我们要进行无限次预测——结果会是什么?如果我们在每一步都丢失信息,我们最终肯定会完全没有信息,我们的概率将平均分布在数组中belief
。让我们尝试 100 次迭代。情节是动画的;使用滑块更改步数。
In [18]:
In [19]:
在 100 次迭代之后,我们几乎丢失了所有信息,即使我们 100% 确定我们从位置 0 开始。你可以随意修改数字以查看不同更新次数的效果。例如,在 100 次更新后留下少量信息,在 50 次更新后留下很多信息,但是经过 200 次迭代,基本上所有信息都丢失了。
而且,如果您在线查看此内容,这里是该输出的动画。

我不会在本书的其余部分生成这些独立的动画。请参阅前言以获取在 Web 上免费运行本书或在您的计算机上安装 IPython 。这将允许您运行所有单元格并查看动画。练习这段代码非常重要,而不仅仅是被动地阅读。
用卷积抽象
我们假设运动误差最多是一个位置。但错误可能是两个、三个或更多位置。作为程序员,我们总是希望抽象我们的代码,使其适用于所有情况。
这很容易用卷积
解决。卷积用另一个函数修改一个函数。在我们的例子中,我们正在使用传感器的误差函数修改概率分布。predict_move()
是一个卷积的实现,虽然我们没有这样称呼它。形式上,卷积定义为
是将 f 与 g 进行卷积的符号。这并不意味着倍增。
积分用于连续函数,但我们使用的是离散函数。我们用求和代替积分,用数组括号代替括号。
predict_move()
正在计算这个方程 - 它计算一系列乘法的总和。
可汗学院对卷积有很好的介绍,维基百科有一些关于卷积的优秀动画。但是大意已经很清楚了。您将一个称为内核
的数组滑过另一个数组,将当前单元格的邻居与第二个数组的值相乘。在上面的示例中,我们使用 0.8 表示移动到正确位置的概率,0.1 表示过低,0.1 表示过高。我们用 array 表示这个内核[0.1, 0.8, 0.1]
。我们需要做的就是编写一个循环遍历数组的每个元素,乘以内核,然后对结果求和。为了强调信念是一种概率分布,我将其命名为pdf
。
In [20]:
代码虽然实现了算法,但运行速度非常慢。SciPy在ndimage.filters
模块中提供了一个卷积方法 convolve()
。我们需要在卷积之前将 pdf 移动 offset
,np.roll()
就可以。这样,移动和预测可以用一行来实现:
FilterPy 使用 discrete_bayes.predict()
函数来实现。
In [21]:

除了中间的元素外,所有元素都没有变化。位置 4 和 6 的值应该是
位置5应该是
让我们确保对大于 1 的运动和不对称内核(下面例子的内核为 [0.05, 0.05, 0.6, 0.2, 0.1])能正确工作。
In [22]:

位置被正确地移动了 3 个位置,我们给过高测量给与更多权重,所以这看起来是正确的。
确保您了解我们在做什么。我们正在预测狗的移动位置,并对概率进行卷积以获得先验。
如果我们不使用概率,我们会使用我之前给出的这个等式:
先验概率是我们对狗在哪里的预测,是狗移动的量加上他当前的位置。例如,狗在 10 点,他移动了 5 米,所以他现在在 15 米处。再简单不过了。但是我们使用概率对此进行建模,所以我们的等式是:
我们正在将当前概率位置估计与我们认为狗移动了多少的概率估计进行卷积。这是相同的概念,但数学略有不同。 粗体表示它是一个数组。