02-Discrete Bayes Filter 离散贝叶斯滤波

02-Discrete Bayes Filter 离散贝叶斯滤波


本文翻译自 ISOEhy/Kalman-and-Bayesian-Filters-in-Python. (github.com),翻译采用的是google机翻,笔者翻译时对一些明显错误的翻译进行了更改,但是还出存在很多比较别的地方,好在保留了原文,两者结合起来可以帮助英文不太好的同学更好的理解文章的内容。部分翻译不算准确,可结合英文原文阅读。

这个系列的文章浅入深出,是学习kalman滤波非常好的参考资源,英文好的可以直接看英文,写得非常易于理解。笔者也会抽时间将这个系列更完。

%matplotlib inline
#format the book
import book_format
book_format.set_style()

The Kalman filter belongs to a family of filters called Bayesian filters. Most textbook treatments of the Kalman filter present the Bayesian formula, perhaps shows how it factors into the Kalman filter equations, but mostly keeps the discussion at a very abstract level.

That approach requires a fairly sophisticated understanding of several fields of mathematics, and it still leaves much of the work of understanding and forming an intuitive grasp of the situation in the hands of the reader.

I will use a different way to develop the topic, to which I owe the work of Dieter Fox and Sebastian Thrun a great debt. It depends on building an intuition on how Bayesian statistics work by tracking an object through a hallway - they use a robot, I use a dog. I like dogs, and they are less predictable than robots which imposes interesting difficulties for filtering. The first published example of this that I can find seems to be Fox 1999 [1], with a fuller example in Fox 2003 [2]. Sebastian Thrun also uses this formulation in his excellent Udacity course Artificial Intelligence for Robotics [3]. In fact, if you like watching videos, I highly recommend pausing reading this book in favor of first few lessons of that course, and then come back to this book for a deeper dive into the topic.

Let’s now use a simple thought experiment, much like we did with the g-h filter, to see how we might reason about the use of probabilities for filtering and tracking.

卡尔曼滤波器属于一个称为“贝叶斯滤波器”的滤波器家族。大多数关于卡尔曼滤波器的教科书处理都提出了贝叶斯公式,也许显示了它是如何
将因子纳入卡尔曼滤波器方程的,但大多将讨论保持在非常抽象的水平上。

这种方法需要对数学的几个领域有相当复杂的理解,而且它仍然把
理解和形成对情况的直观把握的大部分工作交给了读者。

我将用一种不同的方式来发展这个话题,我非常感谢迪特尔·福克斯和塞巴斯蒂安·特龙的工作。这取决于通过走廊跟踪物体来建立贝叶斯统计
如何工作的直觉——他们用的是机器人,我用的是狗。我喜欢狗,它们比机器人更不可预测,这给过滤带来了有趣的困难。我能找到的第一个出
版的例子似乎是Fox 1999[1],Fox 2003[2]中有一个更完整的例子。Sebastian Thrun在他出色的Udacity课程《机器人的人工智能》中
也使用了这个公式[3]。事实上,如果你喜欢看视频,我强烈建议你暂停阅读这本书,转而学习该课程的前几节课,然后再回到这本书中深入
探讨这个话题。

现在让我们使用一个简单的思维实验,就像我们对g-h滤波器所做的那样,看看我们如何推理使用概率进行滤波和跟踪。

Tracking a Dog

Let’s begin with a simple problem. We have a dog friendly workspace, and so people bring their dogs to work. Occasionally the dogs wander out of offices and down the halls. We want to be able to track them. So during a hackathon somebody invented a sonar sensor to attach to the dog’s collar. It emits a signal, listens for the echo, and based on how quickly an echo comes back we can tell whether the dog is in front of an open doorway or not. It also senses when the dog walks, and reports in which direction the dog has moved. It connects to the network via wifi and sends an update once a second.

I want to track my dog Simon, so I attach the device to his collar and then fire up Python, ready to write code to track him through the building. At first blush this may appear impossible. If I start listening to the sensor of Simon’s collar I might read door, hall, hall, and so on. How can I use that information to determine where Simon is?

To keep the problem small enough to plot easily we will assume that there are only 10 positions in the hallway, which we will number 0 to 9, where 1 is to the right of 0. For reasons that will be clear later, we will also assume that the hallway is circular or rectangular. If you move right from position 9, you will be at position 0.

When I begin listening to the sensor I have no reason to believe that Simon is at any particular position in the hallway. From my perspective he is equally likely to be in any position. There are 10 positions, so the probability that he is in any given position is 1/10.

Let’s represent our belief of his position in a NumPy array. I could use a Python list, but NumPy arrays offer functionality that we will be using soon.

让我们从一个简单的问题开始。我们有一个狗狗友好的工作空间,所以人们带着他们的狗来工作。这些狗偶尔会走出办公室,在大厅里游荡.
我们希望能够追踪他们。在黑客马拉松期间,有人发明了一个声纳传感器,安装在狗的项圈上。它发出一个信号,监听回声,根据回声返回的
速度,我们可以判断狗是否在敞开的门口。当狗走路时,它也能感知,并报告狗朝哪个方向移动。它通过wifi连接到网络,并每秒发送一次
更新。

我想跟踪我的狗西蒙,所以我把这个装置绑在他的项圈上,然后启动Python,准备写代码在大楼里跟踪他。乍一看,这似乎是不可能的。如
果我开始听西蒙项圈上的传感器,我可能会读到大厅大厅等等。我怎么用这些信息来确定西蒙在哪里?为了使问题足够
小,便于绘制,我们将假设走廊中只有10个位置,我们将其编号为0到9,其中1在0的右侧。出于稍后会清楚说明的原因,我们还将假设走廊
是圆形或矩形的。如果你从位置9向右移动,你将在位置0。当我开始听传感器时,我没有理由相信西蒙在走廊的任何特定位置。在我看来,他
同样有可能担任任何职位。有10个位置,所以他在任意位置的概率是1/10。

让我们在NumPy数组中表示我们对他位置的信念。我可以使用Python列表,但是NumPy数组提供了我们将很快使用的功能。

import numpy as np
belief = np.array([1/10]*10)
print(belief)
[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]

In Bayesian statistics this is called a prior. It is the probability prior to incorporating measurements or other information. More completely, this is called the prior probability distribution. A probability distribution is a collection of all possible probabilities for an event. Probability distributions always sum to 1 because something had to happen; the distribution lists all possible events and the probability of each.

I’m sure you’ve used probabilities before - as in “the probability of rain today is 30%”. The last paragraph sounds like more of that. But Bayesian statistics was a revolution in probability because it treats probability as a belief about a single event. Let’s take an example. I know that if I flip a fair coin infinitely many times I will get 50% heads and 50% tails. This is called frequentist statistics to distinguish it from Bayesian statistics. Computations are based on the frequency in which events occur.

贝叶斯统计中,这被称为先验。它是在纳入测量或其他信息之前的概率。更完整地说,这叫做先验概率分布。一个概率分布是一个事件的所有可能概率的集合。概率分布的和总是1因为总有事情发生;该分布列出了所有可能发生的事件及其发生的概率。

我相信你以前用过概率——比如“今天下雨的概率是30%”。最后一段听起来更像这样。但贝叶斯统计是概率论的一次革命,因为它将概率视为对单个事件的一种信念。让我们举个例子。我知道如果我抛一枚均匀的硬币无限次我将得到50%正面和50%反面。不同于贝叶斯统计,这被称为频率统计。计算是基于事件发生的频率。

I flip the coin one more time and let it land. Which way do I believe it landed? Frequentist probability has nothing to say about that; it will merely state that 50% of coin flips land as heads. In some ways it is meaningless to assign a probability to the current state of the coin. It is either heads or tails, we just don’t know which. Bayes treats this as a belief about a single event - the strength of my belief or knowledge that this specific coin flip is heads is 50%. Some object to the term “belief”; belief can imply holding something to be true without evidence. In this book it always is a measure of the strength of our knowledge. We’ll learn more about this as we go.

Bayesian statistics takes past information (the prior) into account. We observe that it rains 4 times every 100 days. From this I could state that the chance of rain tomorrow is 1/25. This is not how weather prediction is done. If I know it is raining today and the storm front is stalled, it is likely to rain tomorrow. Weather prediction is Bayesian.

In practice statisticians use a mix of frequentist and Bayesian techniques. Sometimes finding the prior is difficult or impossible, and frequentist techniques rule. In this book we can find the prior. When I talk about the probability of something I am referring to the probability that some specific thing is true given past events. When I do that I’m taking the Bayesian approach.

Now let’s create a map of the hallway. We’ll place the first two doors close together, and then another door further away. We will use 1 for doors, and 0 for walls:

我再抛一次硬币,让它落地。我相信它落在哪个方向?频率论的概率与此无关;它只会说明50%的硬币抛掷为正面。在某些方面,给硬币的当前状态分配一个概率是没有意义的。要么是正面,要么是反面,我们只是不知道是哪个。贝叶斯将其视为对单个事件的信念——我对这个特定的硬币投掷是正面的信念或知识的强度是50%。有些人反对“信仰”这个词;信仰可以意味着在没有证据的情况下坚持某事是正确的。在这本书中,它总是衡量我们知识的力量。我们以后会学到更多。

贝叶斯统计考虑了过去的信息(先验)。我们观察到每100天下雨4次。由此我可以说明天下雨的几率是1/25。天气预报不是这样做的。如果我知道今天下雨,风暴锋停止了,那么明天可能会下雨。天气预报是贝叶斯理论。

在实践中,统计学家混合使用频率论和贝叶斯技术。有时找到先验是困难的或不可能的,而频率论技术占主导地位。在这本书中我们可以找到先验。当我谈论某件事的概率时,我指的是在过去事件的前提下,某件事发生的概率。当我这样做的时候,我用的是贝叶斯方法。
现在我们来画一张走廊的地图。我们将把前两扇门放在一起,然后再把另一扇门放在远处。门用1,墙用0:

hallway = np.array([1, 1, 0, 0, 0, 0, 0, 0, 1, 0])

I start listening to Simon’s transmissions on the network, and the first data I get from the sensor is door. For the moment assume the sensor always returns the correct answer. From this I conclude that he is in front of a door, but which one? I have no reason to believe he is in front of the first, second, or third door. What I can do is assign a probability to each door. All doors are equally likely, and there are three of them, so I assign a probability of 1/3 to each door.

我开始听西蒙在网络上的传输,我从传感器得到的第一个数据是。目前假设传感器总是返回正确的答案。由此我得出结论,他在一扇门前,但哪一扇呢?我没有理由相信他在第一扇、第二扇或第三扇门前。我能做的就是给每扇门分配一个概率。所有门出现的概率都是等的,有三个门,所以我给每个门分配1/3的概率。

import kf_book.book_plots as book_plots
from kf_book.book_plots import figsize, set_figsize
import matplotlib.pyplot as plt

belief = np.array([1/3, 1/3, 0, 0, 0, 0, 0, 0, 1/3, 0])
book_plots.bar_plot(belief)


png

This distribution is called a categorical distribution, which is a discrete distribution describing the probability of observing n n n outcomes. It is a multimodal distribution because we have multiple beliefs about the position of our dog. Of course we are not saying that we think he is simultaneously in three different locations, merely that we have narrowed down our knowledge to one of these three locations. My (Bayesian) belief is that there is a 33.3% chance of being at door 0, 33.3% at door 1, and a 33.3% chance of being at door 8.

This is an improvement in two ways. I’ve rejected a number of hallway positions as impossible, and the strength of my belief in the remaining positions has increased from 10% to 33%. This will always happen. As our knowledge improves the probabilities will get closer to 100%.

这种分布称为分类分布,它是一种描述观察到 n n n个结果的概率的离散分布。这是一个多模态分布,因为我们对狗的位置有多个信念。当然,我们并不是说我们认为他同时在三个不同的地方,只是我们已经将我们的知识范围缩小到这三个地方中的一个。我(贝叶斯)认为有33.3%的概率出现在门0,33.3%的概率出现在门1,33.3%的概率出现在门8。

这在两个方面是一种改进。我已经拒绝了一些不可能的走廊位置,我对剩余位置的信心从10%增加到33%。这种情况总是会发生。随着我们知识的提高,这种可能性将接近100%。

A few words about the mode
of a distribution. Given a list of numbers, such as {1, 2, 2, 2, 3, 3, 4}, the mode is the number that occurs most often. For this set the mode is 2. A distribution can contain more than one mode. The list {1, 2, 2, 2, 3, 3, 4, 4, 4} contains the modes 2 and 4, because both occur three times. We say the former list is unimodal, and the latter is multimodal.

Another term used for this distribution is a histogram. Histograms graphically depict the distribution of a set of numbers. The bar chart above is a histogram.

I hand coded the belief array in the code above. How would we implement this in code? We represent doors with 1, and walls as 0, so we will multiply the hallway variable by the percentage, like so;

关于分布的模式的几句话。给定一个数字列表,例如{1,2,2,3,3,4},模式是出现频率最高的数字。对于这个设置,模式是2。一个分布可以包含多个模式。列表{1,2,2,2,2,3,3,4,4,4}包含模式2和4,因为它们都出现了三次。我们说前一个列表是单模态,后一个列表是多模态

这种分布的另一个术语是直方图。直方图以图形方式描述一组数字的分布。上面的条形图是直方图。我在上面的代码中手工编写了“信念”数组。我们如何在代码中实现它?我们用1表示门,用0表示墙,所以我们将走廊变量乘以百分比,像这样;

belief = hallway * (1/3)
print(belief)
[0.333 0.333 0.    0.    0.    0.    0.    0.    0.333 0.   ]

Extracting Information from Sensor Readings

Let’s put Python aside and think about the problem a bit. Suppose we were to read the following from Simon’s sensor:

  • door
  • move right
  • door

Can we deduce Simon’s location? Of course! Given the hallway’s layout there is only one place from which you can get this sequence, and that is at the left end. Therefore we can confidently state that Simon is in front of the second doorway. If this is not clear, suppose Simon had started at the second or third door. After moving to the right, his sensor would have returned ‘wall’. That doesn’t match the sensor readings, so we know he didn’t start there. We can continue with that logic for all the remaining starting positions. The only possibility is that he is now in front of the second door. Our belief is:

让我们把Python放在一边,思考一下这个问题。假设我们从西蒙的传感器读取以下信息:

  • 向右移动

我们能推断出西蒙的位置吗?当然!考虑到走廊的布局,你只能从一个地方获得这个序列,那就是左端。因此,我们可以自信地说,西蒙在第二个门口的前面。如果不清楚,假设西蒙从第二扇或第三扇门开始。向右移动后,他的传感器会返回“墙”。这和传感器的读数不符,所以我们知道他不是从那里开始的。我们可以对所有剩余的起始位置继续使用这个逻辑。唯一的可能是他现在在第二扇门前。我们的信念是:

belief = np.array([0., 1., 0., 0., 0., 0., 0., 0., 0., 0.])

I designed the hallway layout and sensor readings to give us an exact answer quickly. Real problems are not so clear cut. But this should trigger your intuition - the first sensor reading only gave us low probabilities (0.333) for Simon’s location, but after a position update and another sensor reading we know more about where he is. You might suspect, correctly, that if you had a very long hallway with a large number of doors that after several sensor readings and positions updates we would either be able to know where Simon was, or have the possibilities narrowed down to a small number of possibilities. This is possible when a set of sensor readings only matches one to a few starting locations.

We could implement this solution now, but instead let’s consider a real world complication to the problem.

我设计了走廊布局和传感器读数,以便快速给出准确的答案。真正的问题并不那么明确。但这应该会触发你的直觉——第一个传感器读数只给了我们Simon位置的低概率(0.333),但在位置更新和另一个传感器读数之后,我们更了解他在哪里。你可能会怀疑,如果你有一个很长的走廊,里面有很多扇门,在几次传感器读数和位置更新之后,我们要么能知道西蒙在哪里,要么把可能性缩小到少数几个。当一组传感器读数只匹配一个到几个起始位置时,这是可能的。

我们现在可以实现这个解决方案,但是让我们考虑一下这个问题在现实世界中的复杂情况。

Noisy Sensors

Perfect sensors are rare. Perhaps the sensor would not detect a door if Simon sat in front of it while scratching himself, or misread if he is not facing down the hallway. Thus when I get door I cannot use 1/3 as the probability. I have to assign less than 1/3 to each door, and assign a small probability to each blank wall position. Something like

完美的传感器是罕见的。也许,如果西蒙坐在门前挠痒痒,传感器就不会探测到门,或者如果他没有面向走廊,传感器就会误读。因此,当我得到时,我不能用1/3作为概率。我必须给每个门分配小于1/3的概率,给每个空墙的位置分配一个小概率。类似的

[.31, .31, .01, .01, .01, .01, .01, .01, .31, .01]

At first this may seem insurmountable. If the sensor is noisy it casts doubt on every piece of data. How can we conclude anything if we are always unsure?

The answer, as for the problem above, is with probabilities. We are already comfortable assigning a probabilistic belief to the location of the dog; now we have to incorporate the additional uncertainty caused by the sensor noise.

Say we get a reading of door, and suppose that testing shows that the sensor is 3 times more likely to be right than wrong. We should scale the probability distribution by 3 where there is a door. If we do that the result will no longer be a probability distribution, but we will learn how to fix that in a moment.

Let’s look at that in Python code. Here I use the variable z to denote the measurement. z or y are customary choices in the literature for the measurement. As a programmer I prefer meaningful variable names, but I want you to be able to read the literature and/or other filtering code, so I will start introducing these abbreviated names now.

乍一看,这似乎是无法克服的。如果传感器有噪声,就会对每一条数据产生怀疑。如果我们总是不确定,我们怎么能得出结论呢?对于上面的问题,答案是概率。我们已经习惯了给狗的位置分配一个概率信念;现在我们必须考虑由传感器噪声引起的额外不确定性。假设我们得到的读数,并假设测试表明传感器正确的可能性是错误的3倍。在有门的地方我们应该把概率分布乘以3。如果我们这样做,结果将不再是一个概率分布,但我们将学习如何解决这个问题。

让我们看一下Python代码。这里我用变量z来表示测量值。“z”或“y”是文献中测量的习惯选择。作为一名程序员,我更喜欢有意义的变量名,但我希望您能够阅读文献和/或其他过滤代码,所以我现在将开始介绍这些缩写名。

def update_belief(hall, belief, z, correct_scale):
    for i, val in enumerate(hall):
        if val == z:
            belief[i] *= correct_scale

belief = np.array([0.1] * 10)
reading = 1 # 1 is 'door'
update_belief(hallway, belief, z=reading, correct_scale=3.)
print('belief:', belief)
print('sum =', sum(belief))
plt.figure()
book_plots.bar_plot(belief)
belief: [0.3 0.3 0.1 0.1 0.1 0.1 0.1 0.1 0.3 0.1]
sum = 1.6000000000000003

png

This is not a probability distribution because it does not sum to 1.0. But the code is doing mostly the right thing - the doors are assigned a number (0.3) that is 3 times higher than the walls (0.1). All we need to do is normalize the result so that the probabilities correctly sum to 1.0. Normalization is done by dividing each element by the sum of all elements in the list. That is easy with NumPy:

这不是一个概率分布,因为它的和不等于1.0。但代码做的大多是正确的事情——门被分配了一个数字(0.3),比墙(0.1)高3倍。我们所需要做的就是将结果归一化,使概率正确地和为1.0。规范化是通过将每个元素除以列表中所有元素的和来完成的。这很容易用NumPy:

belief / sum(belief)
array([0.188, 0.188, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.188,
       0.062])

FilterPy implements this with the normalize function:

from filterpy.discrete_bayes import normalize
normalize(belief)

It is a bit odd to say “3 times as likely to be right as wrong”. We are working in probabilities, so let’s specify the probability of the sensor being correct, and compute the scale factor from that. The equation for that is

说“正确的可能性是错误的三倍”有点奇怪。我们用概率计算,所以让我们指定传感器正确的概率,并从中计算比例因子。它的方程是

s c a l e = p r o b c o r r e c t p r o b i n c o r r e c t = p r o b c o r r e c t 1 − p r o b c o r r e c t scale = \frac{prob_{correct}}{prob_{incorrect}} = \frac{prob_{correct}} {1-prob_{correct}} scale=probincorrectprobcorrect=1probcorrectprobcorrect

Also, the for loop is cumbersome. As a general rule you will want to avoid using for loops in NumPy code. NumPy is implemented in C and Fortran, so if you avoid for loops the result often runs 100x faster than the equivalent loop.

How do we get rid of this for loop? NumPy lets you index arrays with boolean arrays. You create a boolean array with logical operators. We can find all the doors in the hallway with:

此外,for循环也很麻烦。作为一般规则,您将希望避免在NumPy代码中使用’ for '循环。NumPy是用C和Fortran实现的,所以如果你避免for循环,结果通常比等效循环快100倍。

我们如何摆脱这个for循环?NumPy允许您使用布尔数组对数组进行索引。创建一个带有逻辑运算符的布尔数组。我们可以找到走廊里所有的门

hallway == 1
array([ True,  True, False, False, False, False, False, False,  True,
       False])

When you use the boolean array as an index to another array it returns only the elements where the index is True. Thus we can replace the for loop with

当你使用布尔数组作为另一个数组的索引时,它只返回索引为“True”的元素。因此,我们可以将for循环替换为

belief[hall==z] *= scale

and only the elements which equal z will be multiplied by scale.

Teaching you NumPy is beyond the scope of this book. I will use idiomatic NumPy constructs and explain them the first time I present them. If you are new to NumPy there are many blog posts and videos on how to use NumPy efficiently and idiomatically.

Here is our improved version:

只有等于’ z ‘的元素才会乘以’ scale '。
教你NumPy超出了本书的范围。我将使用惯用的NumPy结构,并在第一次介绍它们时进行解释。如果你是NumPy的新手,有很多博客文章和视频教你如何高效和习惯地使用NumPy。

这是我们改进后的版本:

from filterpy.discrete_bayes import normalize

def scaled_update(hall, belief, z, z_prob): 
    scale = z_prob / (1. - z_prob)
    belief[hall==z] *= scale
    normalize(belief)

belief = np.array([0.1] * 10)
scaled_update(hallway, belief, z=1, z_prob=.75)

print('sum =', sum(belief))
print('probability of door =', belief[0])
print('probability of wall =', belief[2])
book_plots.bar_plot(belief, ylim=(0, .3))
sum = 1.0
probability of door = 0.1875
probability of wall = 0.06249999999999999

png

We can see from the output that the sum is now 1.0, and that the probability of a door vs wall is still three times larger. The result also fits our intuition that the probability of a door must be less than 0.333, and that the probability of a wall must be greater than 0.0. Finally, it should fit our intuition that we have not yet been given any information that would allow us to distinguish between any given door or wall position, so all door positions should have the same value, and the same should be true for wall positions.

This result is called the posterior, which is short for posterior probability distribution. All this means is a probability distribution after incorporating the measurement information (posterior means ‘after’ in this context). To review, the prior is the probability distribution before including the measurement’s information.

Another term is the likelihood. When we computed belief[hall==z] *= scale we were computing how likely each position was given the measurement. The likelihood is not a probability distribution because it does not sum to one.

The combination of these gives the equation

我们可以从输出中看到,总和现在是1.0,门vs墙的概率仍然是3倍大。结果也符合我们的直觉,即门的概率必须小于0.333,而墙的概率必须大于0.0。最后,它应该符合我们的直觉,我们还没有得到任何信息,可以让我们区分任何给定的门或墙的位置,所以所有的门的位置应该有相同的值,墙的位置也应该是一样的。

这个结果被称为posterior,是后验概率分布 的缩写。所有这些都意味着包含测量信息的概率分布(后验在此上下文中表示“后”)。回顾一下,“先验”是包含测量信息之前的概率分布。

另一个术语是可能性 中文翻译叫似然。当我们计算belief[hall==z] *= scale时,我们是在计算每个位置被给定测量值的可能性。概率不是概率分布,因为它的和不等于1。

它们的组合给出了这个方程

p o s t e r i o r = l i k e l i h o o d × p r i o r n o r m a l i z a t i o n \mathtt{posterior} = \frac{\mathtt{likelihood} \times \mathtt{prior}}{\mathtt{normalization}} posterior=normalizationlikelihood×prior

When we talk about the filter’s output we typically call the state after performing the prediction the prior or prediction, and we call the state after the update either the posterior or the estimated state.

It is very important to learn and internalize these terms as most of the literature uses them extensively.

Does scaled_update() perform this computation? It does. Let me recast it into this form:

当我们讨论过滤器的输出时,我们通常将执行预测后的状态称为先验预测,而将更新后的状态称为后验估计状态

学习和内化这些术语是非常重要的,因为大多数文献都广泛使用它们。

scaled_update()执行此计算吗?它会。让我把它改写成这样:

def scaled_update(hall, belief, z, z_prob): 
    scale = z_prob / (1. - z_prob)
    likelihood = np.ones(len(hall))
    likelihood[hall==z] *= scale
    return normalize(likelihood * belief)

This function is not fully general. It contains knowledge about the hallway, and how we match measurements to it. We always strive to write general functions. Here we will remove the computation of the likelihood from the function, and require the caller to compute the likelihood themselves.

Here is a full implementation of the algorithm:

这个函数不是完全通用的。它包含了关于走廊的知识,以及我们如何将测量结果与之匹配。我们总是努力编写通用函数。这里,我们将从函数中删除可能性的计算,并要求调用者自己计算可能性。

以下是该算法的完整实现:

def update(likelihood, prior):
    return normalize(likelihood * prior)

Computation of the likelihood varies per problem. For example, the sensor might not return just 1 or 0, but a float between 0 and 1 indicating the probability of being in front of a door. It might use computer vision and report a blob shape that you then probabilistically match to a door. It might use sonar and return a distance reading. In each case the computation of the likelihood will be different. We will see many examples of this throughout the book, and learn how to perform these calculations.

FilterPy implements update. Here is the previous example in a fully general form:

每个问题的可能性计算方法不同。例如,传感器可能不只是返回1或0,而是返回一个介于0和1之间的“浮点数”,表示出现在门前的概率。它可能会使用计算机视觉并报告一个斑点形状,然后你可能会将其与门相匹配。它可能会使用声纳并返回距离读数。在每种情况下,可能性的计算将是不同的。我们将在本书中看到许多这样的例子,并学习如何执行这些计算。

FilterPy实现update。下面是前一个例子的完全一般形式:

from filterpy.discrete_bayes import update

def lh_hallway(hall, z, z_prob):
    """ compute likelihood that a measurement matches
    positions in the hallway."""
    
    try:
        scale = z_prob / (1. - z_prob)
    except ZeroDivisionError:
        scale = 1e8

    likelihood = np.ones(len(hall))
    likelihood[hall==z] *= scale
    return likelihood

belief = np.array([0.1] * 10)
likelihood = lh_hallway(hallway, z=1, z_prob=.75)
update(likelihood, belief)  
array([0.188, 0.188, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.188,
       0.062])

Incorporating Movement

Recall how quickly we were able to find an exact solution when we incorporated a series of measurements and movement updates. However, that occurred in a fictional world of perfect sensors. Might we be able to find an exact solution with noisy sensors?

Unfortunately, the answer is no. Even if the sensor readings perfectly match an extremely complicated hallway map, we cannot be 100% certain that the dog is in a specific position - there is, after all, a tiny possibility that every sensor reading was wrong! Naturally, in a more typical situation most sensor readings will be correct, and we might be close to 100% sure of our answer, but never 100% sure. This may seem complicated, but let’s go ahead and program the math.

First let’s deal with the simple case - assume the movement sensor is perfect, and it reports that the dog has moved one space to the right. How would we alter our belief array?

I hope that after a moment’s thought it is clear that we should shift all the values one space to the right. If we previously thought there was a 50% chance of Simon being at position 3, then after he moved one position to the right we should believe that there is a 50% chance he is at position 4. The hallway is circular, so we will use modulo arithmetic to perform the shift.

回想一下,当我们结合一系列测量和运动更新时,我们是如何快速找到精确的解决方案的。然而,这发生在一个拥有完美传感器的虚构世界。也许我们能找到一个带有噪声传感器的精确解决方案呢?

不幸的是,答案是否定的。即使传感器读数与极其复杂的走廊地图完全匹配,我们也不能100%确定狗在特定的位置——毕竟,每个传感器读数都有很小的可能性是错误的!当然,在更典型的情况下,大多数传感器读数都是正确的,我们可能接近100%确定我们的答案,但永远不会100%确定。这可能看起来很复杂,但让我们继续编写数学程序。

首先,让我们来处理一个简单的情况——假设运动传感器是完美的,它报告狗向右移动了一个空间。我们如何改变我们的“信念”数组?我希望经过片刻的思考之后,我们应该将所有的值向右移动一个空间。如果我们之前认为西蒙有50%的几率在位置3,那么在他向右移动一个位置之后我们应该相信他有50%的几率在位置4。走廊是圆形的,所以我们将使用模运算来执行移位。

def perfect_predict(belief, move):
    """ move the position by `move` spaces, where positive is 
    to the right, and negative is to the left
    """
    n = len(belief)
    result = np.zeros(n)
    for i in range(n):
        result[i] = belief[(i-move) % n]
    return result
        
belief = np.array([.35, .1, .2, .3, 0, 0, 0, 0, 0, .05])
plt.subplot(121)
book_plots.bar_plot(belief, title='Before prediction', ylim=(0, .4))

belief = perfect_predict(belief, 1)
plt.subplot(122)
book_plots.bar_plot(belief, title='After prediction', ylim=(0, .4))


png

We can see that we correctly shifted all values one position to the right, wrapping from the end of the array back to the beginning.

The next cell animates this so you can see it in action. Use the slider to move forwards and backwards in time. This simulates Simon walking around and around the hallway. It does not yet incorporate new measurements so the probability distribution does not change shape, only position.

我们可以看到,我们正确地将所有值向右移动了一个位置,从数组的末尾换行到开始。
下一个单元格将其动画化,这样你就可以看到它的动作。使用滑块按时间向前或向后移动。这模拟了西蒙在走廊里走来走去。它还没有纳入新的测量,所以概率分布不会改变形状,只会改变位置。

from ipywidgets import interact, IntSlider

belief = np.array([.35, .1, .2, .3, 0, 0, 0, 0, 0, .05])
perfect_beliefs = []

for _ in range(20):
    # Simon takes one step to the right
    belief = perfect_predict(belief, 1)
    perfect_beliefs.append(belief)

def simulate(time_step):
    book_plots.bar_plot(perfect_beliefs[time_step], ylim=(0, .4))
    plt.show()
    
interact(simulate, time_step=IntSlider(value=0, max=len(perfect_beliefs)-1));
interactive(children=(IntSlider(value=0, description='time_step', max=19), Output()), _dom_classes=('widget-in…

Terminology

Let’s pause a moment to review terminology. I introduced this terminology in the last chapter, but let’s take a second to help solidify your knowledge.

The system is what we are trying to model or filter. Here the system is our dog. The state is its current configuration or value. In this chapter the state is our dog’s position. We rarely know the actual state, so we say our filters produce the estimated state of the system. In practice this often gets called the state, so be careful to understand the context.

One cycle of prediction and updating with a measurement is called the state or system evolution, which is short for time evolution [7]. Another term is system propagation. It refers to how the state of the system changes over time. For filters, time is usually a discrete step, such as 1 second. For our dog tracker the system state is the position of the dog, and the state evolution is the position after a discrete amount of time has passed.

We model the system behavior with the process model. Here, our process model is that the dog moves one or more positions at each time step. This is not a particularly accurate model of how dogs behave. The error in the model is called the system error or process error.

The prediction is our new prior. Time has moved forward and we made a prediction without benefit of knowing the measurements.

Let’s work an example. The current position of the dog is 17 m. Our epoch is 2 seconds long, and the dog is traveling at 15 m/s. Where do we predict he will be in two seconds?

Clearly,

x ˉ = 17 + ( 15 ∗ 2 ) = 47 \begin{aligned} \bar x &= 17 + (15*2) \\ &= 47 \end{aligned} xˉ=17+(152)=47

I use bars over variables to indicate that they are priors (predictions). We can write the equation for the process model like this:

x ˉ k + 1 = f x ( ∙ ) + x k \bar x_{k+1} = f_x(\bullet) + x_k xˉk+1=fx()+xk

x k x_k xk is the current position or state. If the dog is at 17 m then x k = 17 x_k = 17 xk=17.

f x ( ∙ ) f_x(\bullet) fx() is the state propagation function for x. It describes how much the x k x_k xk changes over one time step. For our example it performs the computation 15 ⋅ 2 15 \cdot 2 152 so we would define it as

f x ( v x , t ) = v k t f_x(v_x, t) = v_kt fx(vx,t)=vkt

让我们暂停一下来回顾一下术语。我在上一章中介绍了这个术语,但让我们花点时间来巩固你的知识。“系统”是我们试图建模或过滤的东西。这个系统是我们的狗。state是它的当前配置或值。在本章中,状态是我们的狗的位置。我们很少知道实际状态,所以我们说我们的过滤器产生系统的“估计状态”。在实践中,这通常被称为状态,所以要小心理解上下文。

用一次测量进行预测和更新的一个周期称为状态或系统的“演化”,是“时间演化”的缩写[7]。另一个术语是“系统传播”。它指的是系统的状态如何随时间变化。对于滤波器,时间通常是一个离散的步骤,例如1秒。对于我们的狗跟踪器,系统状态是狗的位置,状态演化是一段离散时间后的位置。

我们用“流程模型”对系统行为进行建模。这里,我们的过程模型是狗在每个时间步移动一个或多个位置。这并不是一个特别准确的狗的行为模型。模型中的错误称为“系统错误”或“过程错误”。预测是我们的新“先验”。时间在流逝,我们在不知道测量结果的情况下做出了预测。我们来举个例子。狗的当前位置是17米。我们的时间是2秒,狗以15米/秒的速度运动。我们预测他两秒钟后会在哪里?

显然,

x ˉ = 17 + ( 15 ∗ 2 ) = 47 \begin{aligned} \bar x &= 17 + (15*2) \\ &= 47 \end{aligned} xˉ=17+(152)=47

我在变量上使用条形来表示它们是先验(预测)。

我们可以这样编写流程模型的方程:
x ˉ k + 1 = f x ( ∙ ) + x k \bar x_{k+1} = f_x(\bullet) + x_k xˉk+1=fx()+xk

x k x_k xk是当前位置或状态。

如果狗在17米,那么 x k = 17 x_k = 17 xk=17
f x ( ∙ ) f_x(\bullet) fx()是x的状态传播函数,它描述了 x k x_k xk在一个时间步长的变化。

对于我们的示例,它执行计算 15 ⋅ 2 15 \cdot 2 152,因此我们将其定义为 f x ( v x , t ) = v k t f_x(v_x, t) = v_k t fx(vx,t)=vkt

Adding Uncertainty to the Prediction

perfect_predict() assumes perfect measurements, but all sensors have noise. What if the sensor reported that our dog moved one space, but he actually moved two spaces, or zero? This may sound like an insurmountable problem, but let’s model it and see what happens.

Assume that the sensor’s movement measurement is 80% likely to be correct, 10% likely to overshoot one position to the right, and 10% likely to undershoot to the left. That is, if the movement measurement is 4 (meaning 4 spaces to the right), the dog is 80% likely to have moved 4 spaces to the right, 10% to have moved 3 spaces, and 10% to have moved 5 spaces.

Each result in the array now needs to incorporate probabilities for 3 different situations. For example, consider the reported movement of 2. If we are 100% certain the dog started from position 3, then there is an 80% chance he is at 5, and a 10% chance for either 4 or 6. Let’s try coding that:

perfect_predict()假设完美测量,但所有传感器都有噪声。如果传感器报告说我们的狗移动了一个空间,但它实际上移动了两个空间,或者零呢?这听起来像是一个无法克服的问题,但让我们建立一个模型,看看会发生什么。

假设传感器的运动测量值有80%可能是正确的,10%可能向右超调一个位置,10%可能向左过调一个位置。也就是说,如果移动测量值为4(即向右移动4个空间),狗有80%的可能向右移动了4个空间,10%的可能移动了3个空间,10%的可能移动了5个空间。

数组中的每个结果现在都需要包含3种不同情况的概率。例如,考虑报告的2的移动。如果我们100%确定狗从位置3开始,那么它有80%的概率在位置5,而4或6的概率是10%。让我们试着编码:

def predict_move(belief, move, p_under, p_correct, p_over):
    n = len(belief)
    prior = np.zeros(n)
    for i in range(n):
        prior[i] = (
            belief[(i-move) % n]   * p_correct +
            belief[(i-move-1) % n] * p_over +
            belief[(i-move+1) % n] * p_under)      
    return prior

belief = [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.]
prior = predict_move(belief, 2, .1, .8, .1)
book_plots.plot_belief_vs_prior(belief, prior)


png

It appears to work correctly. Now what happens when our belief is not 100% certain?

belief = [0, 0, .4, .6, 0, 0, 0, 0, 0, 0]
prior = predict_move(belief, 2, .1, .8, .1)
book_plots.plot_belief_vs_prior(belief, prior)
prior
array([0.  , 0.  , 0.  , 0.04, 0.38, 0.52, 0.06, 0.  , 0.  , 0.  ])


png

Here the results are more complicated, but you should still be able to work it out in your head. The 0.04 is due to the possibility that the 0.4 belief undershot by 1. The 0.38 is due to the following: the 80% chance that we moved 2 positions (0.4 × \times × 0.8) and the 10% chance that we undershot (0.6 × \times × 0.1). Overshooting plays no role here because if we overshot both 0.4 and 0.6 would be past this position. I strongly suggest working some examples until all of this is very clear, as so much of what follows depends on understanding this step.

If you look at the probabilities after performing the update you might be dismayed. In the example above we started with probabilities of 0.4 and 0.6 in two positions; after performing the update the probabilities are not only lowered, but they are strewn out across the map.

This is not a coincidence, or the result of a carefully chosen example - it is always true of the prediction. If the sensor is noisy we lose some information on every prediction. Suppose we were to perform the prediction an infinite number of times - what would the result be? If we lose information on every step, we must eventually end up with no information at all, and our probabilities will be equally distributed across the belief array. Let’s try this with 100 iterations. The plot is animated; use the slider to change the step number.

这里的结果更复杂,但你应该仍然能够在你的头脑中计算出来。0.04是由于0.4的信念比1低的可能性。0.38是由于以下原因:我们有80%的机会移动了2个位置(0.4美元乘以0.8美元)和10%的机会我们低于(0.6美元乘以0.1美元)。超调在这里不起作用,因为如果我们超调0.4和0.6就会超过这个位置。我强烈建议做一些例子,直到所有这些都非常清楚,因为下面的很多东西取决于理解这一步。

如果您在执行更新后查看概率,您可能会感到沮丧。在上面的例子中,我们从两个位置的0.4和0.6的概率开始;在执行更新后,概率不仅会降低,而且会分散在地图上。

这不是巧合,也不是一个精心挑选的例子的结果——预测总是正确的。如果传感器有噪声,我们会在每次预测中丢失一些信息。假设我们要进行无限次预测,结果会是什么?如果我们在每一步中都丢失信息,那么我们最终将完全没有信息,我们的概率将在“信念”数组中均匀分布。让我们尝试100次迭代。情节是动画的;使用滑块更改步长。

belief = np.array([1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
predict_beliefs = []
    
for i in range(100):
    belief = predict_move(belief, 1, .1, .8, .1)
    predict_beliefs.append(belief)

print('Final Belief:', belief)

# make interactive plot
def show_prior(step):
    book_plots.bar_plot(predict_beliefs[step-1])
    plt.title(f'Step {step}')
    plt.show()

interact(show_prior, step=IntSlider(value=1, max=len(predict_beliefs)));
Final Belief: [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]



interactive(children=(IntSlider(value=1, description='step'), Output()), _dom_classes=('widget-interact',))
print('Final Belief:', belief)
Final Belief: [0.104 0.103 0.101 0.099 0.097 0.096 0.097 0.099 0.101 0.103]

After 100 iterations we have lost almost all information, even though we were 100% sure that we started in position 0. Feel free to play with the numbers to see the effect of differing number of updates. For example, after 100 updates a small amount of information is left, after 50 a lot is left, but by 200 iterations essentially all information is lost.

在100次迭代之后,我们已经丢失了几乎所有的信息,即使我们100%确定我们从位置0开始。你可以随意摆弄这些数字,看看不同数量的更新所产生的效果。例如,在100次更新之后,只剩下少量的信息,在50次更新之后,剩下很多信息,但是在200次迭代之后,基本上所有的信息都丢失了。

And, if you are viewing this online here is an animation of that output.

I will not generate these standalone animations through the rest of the book. Please see the preface for instructions to run this book on the web, for free, or install IPython on your computer. This will allow you to run all of the cells and see the animations. It’s very important that you practice with this code, not just read passively.

我不会在本书的其余部分生成这些独立的动画。请参阅前言,了解在web上免费运行本书的说明,或者在您的计算机上安装IPython。这将允许您运行所有单元格并查看动画。练习这段代码是非常重要的,而不是被动地阅读。

Generalizing with Convolution

We made the assumption that the movement error is at most one position. But it is possible for the error to be two, three, or more positions. As programmers we always want to generalize our code so that it works for all cases.

This is easily solved with convolution. Convolution modifies one function with another function. In our case we are modifying a probability distribution with the error function of the sensor. The implementation of predict_move() is a convolution, though we did not call it that. Formally, convolution is defined as

我们假设运动误差最多为一个位置。但误差可能是两个、三个或更多个位置。作为程序员,我们总是希望泛化我们的代码,使其适用于所有情况。这很容易通过convolution解决。
卷积用一个函数修饰另一个函数。在我们的例子中,我们用传感器的误差函数修改一个概率分布。predict_move()的实现是一个卷积,尽管我们没有这样称呼它。形式上,卷积被定义为

( f ∗ g ) ( t ) = ∫ 0 t  ⁣ f ( τ )   g ( t − τ )   d τ (f \ast g) (t) = \int_0^t \!f(\tau) \, g(t-\tau) \, \mathrm{d}\tau (fg)(t)=0tf(τ)g(tτ)dτ

where f ∗ g f\ast g fg is the notation for convolving f by g. It does not mean multiply.

Integrals are for continuous functions, but we are using discrete functions. We replace the integral with a summation, and the parenthesis with array brackets.

其中 f ∗ g f\ast g fg是对f × g进行卷积的符号,它并不表示乘法。
积分适用于连续函数,但我们用的是离散函数。我们用求和代替积分,用数组括号代替圆括号。

( f ∗ g ) [ t ] = ∑ τ = 0 t  ⁣ f [ τ ]   g [ t − τ ] (f \ast g) [t] = \sum\limits_{\tau=0}^t \!f[\tau] \, g[t-\tau] (fg)[t]=τ=0tf[τ]g[tτ]

Comparison shows that predict_move() is computing this equation - it computes the sum of a series of multiplications.

Khan Academy [4] has a good introduction to convolution, and Wikipedia has some excellent animations of convolutions [5]. But the general idea is already clear. You slide an array called the kernel across another array, multiplying the neighbors of the current cell with the values of the second array. In our example above we used 0.8 for the probability of moving to the correct location, 0.1 for undershooting, and 0.1 for overshooting. We make a kernel of this with the array [0.1, 0.8, 0.1]. All we need to do is write a loop that goes over each element of our array, multiplying by the kernel, and summing the results. To emphasize that the belief is a probability distribution I have named it pdf.

比较表明predict_move()正在计算这个方程—它计算一系列乘法的和。
可汗学院[4]对卷积有很好的介绍,维基百科也有一些很棒的卷积动画[5]。但总体思路已经很清楚了。将一个称为“内核”的数组滑动到另一个数组,将当前单元格的邻居与第二个数组的值相乘。

在上面的例子中,我们使用0.8表示移动到正确位置的概率,0.1表示过冲,0.1表示过冲。我们用数组’[0.1,0.8,0.1]'创建一个内核。我们所需要做的就是编写一个循环,遍历数组的每个元素,乘以内核,然后对结果求和。为了强调信念是一个概率分布,我将其命名为pdf

def predict_move_convolution(pdf, offset, kernel):
    N = len(pdf)
    kN = len(kernel)
    width = int((kN - 1) / 2)

    prior = np.zeros(N)
    for i in range(N):
        for k in range (kN):
            index = (i + (width-k) - offset) % N
            prior[i] += pdf[index] * kernel[k]
    return prior

This illustrates the algorithm, but it runs very slow. SciPy provides a convolution routine convolve() in the ndimage.filters module. We need to shift the pdf by offset before convolution; np.roll() does that. The move and predict algorithm can be implemented with one line:

这说明了算法,但它运行得很慢。SciPy在nimage 中提供了一个卷积例程convolve() 。过滤器的模块。我们需要在卷积之前通过偏移量来移动pdf;np.roll() 就是这样做的。move和predict算法只需一行即可实现:

convolve(np.roll(pdf, offset), kernel, mode='wrap')

FilterPy implements this with discrete_bayespredict() function.

from filterpy.discrete_bayes import predict

belief = [.05, .05, .05, .05, .55, .05, .05, .05, .05, .05]
prior = predict(belief, offset=1, kernel=[.1, .8, .1])
book_plots.plot_belief_vs_prior(belief, prior, ylim=(0,0.6))


png

All of the elements are unchanged except the middle ones. The values in position 4 and 6 should be
( 0.1 × 0.05 ) + ( 0.8 × 0.05 ) + ( 0.1 × 0.55 ) = 0.1 (0.1 \times 0.05)+ (0.8 \times 0.05) + (0.1 \times 0.55) = 0.1 (0.1×0.05)+(0.8×0.05)+(0.1×0.55)=0.1

Position 5 should be ( 0.1 × 0.05 ) + ( 0.8 × 0.55 ) + ( 0.1 × 0.05 ) = 0.45 (0.1 \times 0.05) + (0.8 \times 0.55)+ (0.1 \times 0.05) = 0.45 (0.1×0.05)+(0.8×0.55)+(0.1×0.05)=0.45

Let’s ensure that it shifts the positions correctly for movements greater than one and for asymmetric kernels.

prior = predict(belief, offset=3, kernel=[.05, .05, .6, .2, .1])
book_plots.plot_belief_vs_prior(belief, prior, ylim=(0,0.6))


png

The position was correctly shifted by 3 positions and we give more weight to the likelihood of an overshoot vs an undershoot, so this looks correct.

Make sure you understand what we are doing. We are making a prediction of where the dog is moving, and convolving the probabilities to get the prior.

If we weren’t using probabilities we would use this equation that I gave earlier:

这个位置被正确地移动了3个位置,我们给了超调的可能性比过调的可能性更多的权重,所以这看起来是正确的。
确保你明白我们在做什么。我们正在预测狗的移动方向,并对概率进行卷积以得到先验。

如果我们不用概率我们就用我之前给出的这个方程

x ˉ k + 1 = x k + f x ( ∙ ) \bar x_{k+1} = x_k + f_{\mathbf x}(\bullet) xˉk+1=xk+fx()

The prior, our prediction of where the dog will be, is the amount the dog moved plus his current position. The dog was at 10, he moved 5 meters, so he is now at 15 m. It couldn’t be simpler. But we are using probabilities to model this, so our equation is:

先验,即我们对狗的位置的预测,是狗移动的量加上它当前的位置。狗在10米的地方,他移动了5米,所以他现在在15米的地方。再简单不过了。但是我们用概率来建模,所以我们的方程是:

x ˉ k + 1 = x k ∗ f x ( ∙ ) \bar{ \mathbf x}_{k+1} = \mathbf x_k \ast f_{\mathbf x}(\bullet) xˉk+1=xkfx()

We are convolving the current probabilistic position estimate with a probabilistic estimate of how much we think the dog moved. It’s the same concept, but the math is slightly different. x \mathbf x x is bold to denote that it is an array of numbers.

我们将当前的概率位置估计与我们认为狗移动了多少的概率估计进行卷积。这是相同的概念,但数学略有不同。 x \mathbf x x以粗体表示它是一个数字数组。

Integrating Measurements and Movement Updates

The problem of losing information during a prediction may make it seem as if our system would quickly devolve into having no knowledge. However, each prediction is followed by an update where we incorporate the measurement into the estimate. The update improves our knowledge. The output of the update step is fed into the next prediction. The prediction degrades our certainty. That is passed into another update, where certainty is again increased.

Let’s think about this intuitively. Consider a simple case - you are tracking a dog while he sits still. During each prediction you predict he doesn’t move. Your filter quickly converges on an accurate estimate of his position. Then the microwave in the kitchen turns on, and he goes streaking off. You don’t know this, so at the next prediction you predict he is in the same spot. But the measurements tell a different story. As you incorporate the measurements your belief will be smeared along the hallway, leading towards the kitchen. On every epoch (cycle) your belief that he is sitting still will get smaller, and your belief that he is inbound towards the kitchen at a startling rate of speed increases.

That is what intuition tells us. What does the math tell us?

We have already programmed the update and predict steps. All we need to do is feed the result of one into the other, and we will have implemented a dog tracker!!! Let’s see how it performs. We will input measurements as if the dog started at position 0 and moved right one position each epoch. As in a real world application, we will start with no knowledge of his position by assigning equal probability to all positions.

在预测过程中丢失信息的问题可能会让我们的系统看起来好像很快就会退化到毫无知识的地步。然而,每个预测之后都有一个更新,我们将测量合并到估计中。这个更新提高了我们的知识。更新步骤的输出被输入到下一个预测中。预测降低了我们的确定性。这被传递到另一个更新中,其中确定性再次增加。

让我们直观地思考一下。考虑一个简单的例子——你正在跟踪一只一动不动的狗。在每次预测中,你都预测他不会移动。你的过滤器迅速“收敛”到对他位置的准确估计。然后厨房里的微波炉打开了,他就开始裸奔了。你不知道这个,所以在下一次预测时,你预测他在同一个位置。但测量结果却告诉我们一个不同的故事。当你把这些尺寸结合在一起时,你的信念就会沿着走廊被涂抹,一直延伸到厨房。每过一个周期,你相信他坐着不动的信念就会减弱,而你相信他正以惊人的速度朝厨房冲去的信念就会增强。

这是直觉告诉我们的。数学告诉我们什么?

我们已经编写了更新程序并预测了步骤。我们所需要做的就是将其中一个的结果输入另一个,这样我们就实现了一个狗追踪器!!让我们看看它是如何执行的。我们将输入测量值,就好像狗从位置0开始,每个epoch向右移动一个位置。在现实世界的应用程序中,我们将从不知道他的位置开始,通过为所有位置分配相同的概率。

from filterpy.discrete_bayes import update

hallway = np.array([1, 1, 0, 0, 0, 0, 0, 0, 1, 0])
prior = np.array([.1] * 10)
likelihood = lh_hallway(hallway, z=1, z_prob=.75)
posterior = update(likelihood, prior)
book_plots.plot_prior_vs_posterior(prior, posterior, ylim=(0,.5))


png

After the first update we have assigned a high probability to each door position, and a low probability to each wall position.

在第一次更新之后,我们给每个门的位置分配了一个高概率,给每个墙的位置分配了一个低概率。

kernel = (.1, .8, .1)
prior = predict(posterior, 1, kernel)
book_plots.plot_prior_vs_posterior(prior, posterior, True, ylim=(0,.5))


png

The predict step shifted these probabilities to the right, smearing them about a bit. Now let’s look at what happens at the next sense.

预测步骤将这些概率向右移动,使它们模糊了一点。现在让我们看看下一种感觉会发生什么。

likelihood = lh_hallway(hallway, z=1, z_prob=.75)
posterior = update(likelihood, prior)
book_plots.plot_prior_vs_posterior(prior, posterior, ylim=(0,.5))


png

Notice the tall bar at position 1. This corresponds with the (correct) case of starting at position 0, sensing a door, shifting 1 to the right, and sensing another door. No other positions make this set of observations as likely. Now we will add an update and then sense the wall.

注意位置1上的高条。这对应于从位置0开始,感知一扇门,向右移动1,感知另一扇门的(正确)情况。没有其他立场能使这组观察结果如此可信。现在我们将添加一个更新,然后感知墙。

prior = predict(posterior, 1, kernel)
likelihood = lh_hallway(hallway, z=0, z_prob=.75)
posterior = update(likelihood, prior)
book_plots.plot_prior_vs_posterior(prior, posterior, ylim=(0,.5))


png

This is exciting! We have a very prominent bar at position 2 with a value of around 35%. It is over twice the value of any other bar in the plot, and is about 4% larger than our last plot, where the tallest bar was around 31%. Let’s see one more cycle.

这太令人兴奋了!我们在2号位置有一个非常突出的条,它的值大约是35%。它是该地块中任何其他柱形值的两倍多,并且比我们的上一个地块大约4%,其中最高的柱形约为31%。再看一个循环。

prior = predict(posterior, 1, kernel)
likelihood = lh_hallway(hallway, z=0, z_prob=.75)
posterior = update(likelihood, prior)
book_plots.plot_prior_vs_posterior(prior, posterior, ylim=(0,.5))


png

I ignored an important issue. Earlier I assumed that we had a motion sensor for the predict step; then, when talking about the dog and the microwave I assumed that you had no knowledge that he suddenly began running. I mentioned that your belief that the dog is running would increase over time, but I did not provide any code for this. In short, how do we detect and/or estimate changes in the process model if we aren’t directly measuring it?

For now I want to ignore this problem. In later chapters we will learn the mathematics behind this estimation; for now it is a large enough task just to learn this algorithm. It is profoundly important to solve this problem, but we haven’t yet built enough of the mathematical apparatus that is required, and so for the remainder of the chapter we will ignore the problem by assuming we have a sensor that senses movement.

我忽略了一个重要的问题。之前我假设我们有一个运动传感器用于预测步骤;然后,当谈到狗和微波炉的时候,我以为你不知道它突然跑了起来。我提到过,你认为狗在跑的信念会随着时间的推移而增强,但我没有提供任何代码。简而言之,如果我们不直接测量过程模型中的变化,我们如何检测和/或评估它呢?

现在我想忽略这个问题。在后面的章节中,我们将学习这种估计背后的数学原理;现在学习这个算法是一个很大的任务。解决这个问题是非常重要的,但我们还没有建立足够的数学设备,所以在本章的其余部分,我们将忽略这个问题,假设我们有一个传感器来感知运动。

The Discrete Bayes Algorithm

This chart illustrates the algorithm:

book_plots.predict_update_chart()


png

This filter is a form of the g-h filter. Here we are using the percentages for the errors to implicitly compute the g g g and h h h parameters. We could express the discrete Bayes algorithm as a g-h filter, but that would obscure the logic of this filter.

The filter equations are:

这个滤波器是g-h滤波器的一种形式。这里我们使用误差的百分比来隐式地计算参数g和h。我们可以将离散贝叶斯算法表示为g-h滤波器,但这会模糊该滤波器的逻辑。

滤波器方程为:

x ˉ = x ∗ f x ( ∙ )    Predict Step x = ∥ L ⋅ x ˉ ∥    Update Step \begin{aligned} \bar {\mathbf x} &= \mathbf x \ast f_{\mathbf x}(\bullet)\, \, &\text{Predict Step} \\ \mathbf x &= \|\mathcal L \cdot \bar{\mathbf x}\|\, \, &\text{Update Step}\end{aligned} xˉx=xfx()=LxˉPredict StepUpdate Step

L \mathcal L L is the usual way to write the likelihood function, so I use that. The ∥ ∥ \|\| ∥∥ notation denotes taking the norm. We need to normalize the product of the likelihood with the prior to ensure x x x is a probability distribution that sums to one.

We can express this in pseudocode.

L \mathcal L L是表示似然函数的常用方法,所以我用它。 ∥ ∥ \|\| ∥∥符号表示取范数。我们需要将可能性与先验的乘积归一化,以确保x是一个总和为1的概率分布。

我们可以用伪代码来表示。

Initialization

1. Initialize our belief in the state

Predict

1. Based on the system behavior, predict state for the next time step
2. Adjust belief to account for the uncertainty in prediction

Update

1. Get a measurement and associated belief about its accuracy
2. Compute how likely it is the measurement matches each state
3. Update state belief with this likelihood

初始化

1.初始化我们在状态中的信念

Predict

1.根据系统行为,预测下一个时间步骤2的状态。
2.调整信念以考虑预测中的不确定性

更新

1.获得一个测量值和与之相关的关于其准确性的信念
2.计算测量值与每个状态匹配的可能性。
3.用这个可能性更新状态信念

When we cover the Kalman filter we will use this exact same algorithm; only the details of the computation will differ.

Algorithms in this form are sometimes called predictor correctors. We make a prediction, then correct them.

Let’s animate this. First Let’s write functions to perform the filtering and to plot the results at any step. I’ve plotted the position of the doorways in black. Prior are drawn in orange, and the posterior in blue. I draw a thick vertical line to indicate where Simon really is. This is not an output of the filter - we know where Simon is only because we are simulating his movement.

当我们讲到卡尔曼滤波器时我们将使用完全相同的算法;只有计算的细节会有所不同。

这种形式的算法有时被称为“预测校正器”。我们做出预测,然后修正。

让我们把它做成动画。首先,让我们编写函数来执行过滤,并绘制出每一步的结果。我用黑色标出了门道的位置。先验用橙色表示,后验用蓝色表示。我画了一条粗竖线来表示西蒙的真实位置。这不是过滤器的输出——我们知道Simon在哪里只是因为我们模拟了他的运动。

def discrete_bayes_sim(prior, kernel, measurements, z_prob, hallway):
    posterior = np.array([.1]*10)
    priors, posteriors = [], []
    for i, z in enumerate(measurements):
        prior = predict(posterior, 1, kernel)
        priors.append(prior)

        likelihood = lh_hallway(hallway, z, z_prob)
        posterior = update(likelihood, prior)
        posteriors.append(posterior)
    return priors, posteriors


def plot_posterior(hallway, posteriors, i):
    plt.title('Posterior')
    book_plots.bar_plot(hallway, c='k')
    book_plots.bar_plot(posteriors[i], ylim=(0, 1.0))
    plt.axvline(i % len(hallway), lw=5)
    plt.show()
    
def plot_prior(hallway, priors, i):
    plt.title('Prior')
    book_plots.bar_plot(hallway, c='k')
    book_plots.bar_plot(priors[i], ylim=(0, 1.0), c='#ff8015')
    plt.axvline(i % len(hallway), lw=5)
    plt.show()

def animate_discrete_bayes(hallway, priors, posteriors):
    def animate(step):
        step -= 1
        i = step // 2    
        if step % 2 == 0:
            plot_prior(hallway, priors, i)
        else:
            plot_posterior(hallway, posteriors, i)
    
    return animate

Let’s run the filter and animate it.

# change these numbers to alter the simulation
kernel = (.1, .8, .1)
z_prob = 1.0
hallway = np.array([1, 1, 0, 0, 0, 0, 0, 0, 1, 0])

# measurements with no noise
zs = [hallway[i % len(hallway)] for i in range(50)]

priors, posteriors = discrete_bayes_sim(prior, kernel, zs, z_prob, hallway)
interact(animate_discrete_bayes(hallway, priors, posteriors), step=IntSlider(value=1, max=len(zs)*2));
interactive(children=(IntSlider(value=1, description='step'), Output()), _dom_classes=('widget-interact',))

Now we can see the results. You can see how the prior shifts the position and reduces certainty, and the posterior stays in the same position and increases certainty as it incorporates the information from the measurement. I’ve made the measurement perfect with the line z_prob = 1.0; we will explore the effect of imperfect measurements in the next section. Finally,

Another thing to note is how accurate our estimate becomes when we are in front of a door, and how it degrades when in the middle of the hallway. This should make intuitive sense. There are only a few doorways, so when the sensor tells us we are in front of a door this boosts our certainty in our position. A long stretch of no doors reduces our certainty.

现在我们可以看到结果了。您可以看到先验如何改变位置并降低确定性,而后验如何保持相同位置并增加确定性,因为它包含来自测量的信息。我用 z_prob = 1.0 行使测量变得完美;我们将在下一节探讨不完美测量的影响。

最后,另一件需要注意的事情是,当我们在门前时,我们的估计变得多么准确,而在走廊中间时,它又如何降低。这应该具有直觉意义。只有几个门口,所以当传感器告诉我们我们在门前时,这会增强我们对自己位置的确定性。长时间没有门降低了我们的确定性。

The Effect of Bad Sensor Data

You may be suspicious of the results above because I always passed correct sensor data into the functions. However, we are claiming that this code implements a filter - it should filter out bad sensor measurements. Does it do that?

To make this easy to program and visualize I will change the layout of the hallway to mostly alternating doors and hallways, and run the algorithm on 6 correct measurements:

错误传感器数据的影响

您可能会对上面的结果产生怀疑,因为我总是将正确的传感器数据传递到函数中。然而,我们声称此代码实现了一个过滤器 - 它应该过滤掉不良的传感器测量值。它这样做吗?

为了使其易于编程和可视化,我将走廊的布局更改为大部分交替的门和走廊,并在 6 个正确的测量值上运行算法:

hallway = np.array([1, 0, 1, 0, 0]*2)
kernel = (.1, .8, .1)
prior = np.array([.1] * 10)
zs = [1, 0, 1, 0, 0, 1]
z_prob = 0.75
priors, posteriors = discrete_bayes_sim(prior, kernel, zs, z_prob, hallway)
interact(animate_discrete_bayes(hallway, priors, posteriors), step=IntSlider(value=12, max=len(zs)*2));
interactive(children=(IntSlider(value=12, description='step', max=12), Output()), _dom_classes=('widget-intera…

We have identified the likely cases of having started at position 0 or 5, because we saw this sequence of doors and walls: 1,0,1,0,0. Now I inject a bad measurement. The next measurement should be 0, but instead we get a 1:

我们已经确定了从位置 0 或 5 开始的可能情况,因为我们看到了门和墙的顺序:1,0,1,0,0。现在我注入了一个错误的测量值。下一个测量值应该是 0,但我们得到的是 1:

measurements = [1, 0, 1, 0, 0, 1, 1]
priors, posteriors = discrete_bayes_sim(prior, kernel, measurements, z_prob, hallway);
plot_posterior(hallway, posteriors, 6)


png

That one bad measurement has significantly eroded our knowledge. Now let’s continue with a series of correct measurements.

那一次糟糕的测量极大地侵蚀了我们的知识。现在让我们继续进行一系列正确的测量。

with figsize(y=5.5):
    measurements = [1, 0, 1, 0, 0, 1, 1, 1, 0, 0]
    for i, m in enumerate(measurements):
        likelihood = lh_hallway(hallway, z=m, z_prob=.75)
        posterior = update(likelihood, prior)
        prior = predict(posterior, 1, kernel)
        plt.subplot(5, 2, i+1)
        book_plots.bar_plot(posterior, ylim=(0, .4), title=f'step {i+1}')
    plt.tight_layout()


png

We quickly filtered out the bad sensor reading and converged on the most likely positions for our dog.

Drawbacks and Limitations

Do not be mislead by the simplicity of the examples I chose. This is a robust and complete filter, and you may use the code in real world solutions. If you need a multimodal, discrete filter, this filter works.

With that said, this filter it is not used often because it has several limitations. Getting around those limitations is the motivation behind the chapters in the rest of this book.

The first problem is scaling. Our dog tracking problem used only one variable, p o s pos pos, to denote the dog’s position. Most interesting problems will want to track several things in a large space. Realistically, at a minimum we would want to track our dog’s ( x , y ) (x,y) (x,y) coordinate, and probably his velocity ( x ˙ , y ˙ ) (\dot{x},\dot{y}) (x˙,y˙) as well. We have not covered the multidimensional case, but instead of an array we use a multidimensional grid to store the probabilities at each discrete location. Each update() and predict() step requires updating all values in the grid, so a simple four variable problem would require O ( n 4 ) O(n^4) O(n4) running time per time step. Realistic filters can have 10 or more variables to track, leading to exorbitant computation requirements.

The second problem is that the filter is discrete, but we live in a continuous world. The histogram requires that you model the output of your filter as a set of discrete points. A 100 meter hallway requires 10,000 positions to model the hallway to 1cm accuracy. So each update and predict operation would entail performing calculations for 10,000 different probabilities. It gets exponentially worse as we add dimensions. A 100x100 m 2 ^2 2 courtyard requires 100,000,000 bins to get 1cm accuracy.

A third problem is that the filter is multimodal. In the last example we ended up with strong beliefs that the dog was in position 4 or 9. This is not always a problem. Particle filters, which we will study later, are multimodal and are often used because of this property. But imagine if the GPS in your car reported to you that it is 40% sure that you are on D street, and 30% sure you are on Willow Avenue.

A forth problem is that it requires a measurement of the change in state. We need a motion sensor to detect how much the dog moves. There are ways to work around this problem, but it would complicate the exposition of this chapter, so, given the aforementioned problems, I will not discuss it further.

With that said, if I had a small problem that this technique could handle I would choose to use it; it is trivial to implement, debug, and understand, all virtues.

缺点和限制

不要被我选择的示例的简单性误导。这是一个强大而完整的过滤器,您可以在实际解决方案中使用该代码。如果您需要一个多模式、离散的过滤器,这个过滤器就可以工作。
话虽如此,这个过滤器并不经常使用,因为它有几个限制。克服这些限制是本书其余章节的动机。

第一个问题是缩放。我们的狗追踪问题仅使用一个变量 p o s pos pos 来表示狗的位置。大多数有趣的问题都希望在一个大空间中跟踪多个事物。实际上,至少我们想要跟踪我们狗的 ( x , y ) (x,y) (x,y) 坐标,可能还有他的速度 ( x ˙ , y ˙ ) (\dot{x},\dot{y}) (x˙,y˙)。我们没有涉及多维情况,但是我们使用多维网格来存储每个离散位置的概率,而不是数组。每个 update()predict() 步骤都需要更新网格中的所有值,因此一个简单的四变量问题将需要 O ( n 4 ) O(n^4) O(n4) 运行时间每个时间步。现实的过滤器可能有 10 个或更多变量要跟踪,从而导致过高的计算要求。

第二个问题是滤波器是离散的,但我们生活在一个连续的世界中。直方图要求您将过滤器的输出建模为一组离散点。 100 米长的走廊需要 10,000 个位置才能将走廊建模到 1 厘米的精度。因此,每个更新和预测操作都需要对 10,000 个不同的概率进行计算。随着维度的增加,情况会呈指数级恶化。一个 100x100 m 2 ^2 2 的庭院需要 100,000,000 个垃圾桶才能达到 1 厘米的精度。

第三个问题是过滤器是多模式的。在最后一个例子中,我们最终坚信狗在位置 4 或 9。这并不总是一个问题。我们稍后将研究的粒子滤波器是多峰的,并且由于这个特性而经常被使用。但是想象一下,如果你车上的 GPS 向你报告说它有 40% 的把握你在 D 街,有 30% 的把握你在 Willow Avenue。第四个问题是它需要测量状态的变化。我们需要一个运动传感器来检测狗移动了多少。有一些方法可以解决这个问题,但这会使本章的阐述变得复杂,因此,鉴于上述问题,我将不再进一步讨论。

话虽如此,如果我遇到这种技术可以解决的小问题,我会选择使用它;实施、调试和理解所有优点都是微不足道的。

Tracking and Control

We have been passively tracking an autonomously moving object. But consider this very similar problem. I am automating a warehouse and want to use robots to collect all of the items for a customer’s order. Perhaps the easiest way to do this is to have the robots travel on a train track. I want to be able to send the robot a destination and have it go there. But train tracks and robot motors are imperfect. Wheel slippage and imperfect motors means that the robot is unlikely to travel to exactly the position you command. There is more than one robot, and we need to know where they all are so we do not cause them to crash.

So we add sensors. Perhaps we mount magnets on the track every few feet, and use a Hall sensor to count how many magnets are passed. If we count 10 magnets then the robot should be at the 10th magnet. Of course it is possible to either miss a magnet or to count it twice, so we have to accommodate some degree of error. We can use the code from the previous section to track our robot since magnet counting is very similar to doorway sensing.

But we are not done. We’ve learned to never throw information away. If you have information you should use it to improve your estimate. What information are we leaving out? We know what control inputs we are feeding to the wheels of the robot at each moment in time. For example, let’s say that once a second we send a movement command to the robot - move left 1 unit, move right 1 unit, or stand still. If I send the command ‘move left 1 unit’ I expect that in one second from now the robot will be 1 unit to the left of where it is now. This is a simplification because I am not taking acceleration into account, but I am not trying to teach control theory. Wheels and motors are imperfect. The robot might end up 0.9 units away, or maybe 1.2 units.

Now the entire solution is clear. We assumed that the dog kept moving in whatever direction he was previously moving. That is a dubious assumption for my dog! Robots are far more predictable. Instead of making a dubious prediction based on assumption of behavior we will feed in the command that we sent to the robot! In other words, when we call predict() we will pass in the commanded movement that we gave the robot along with a kernel that describes the likelihood of that movement.

跟踪和控制

我们一直在被动地跟踪一个自主移动的物体。但是请考虑这个非常相似的问题。我正在使仓库自动化,并希望使用机器人来收集客户订单的所有物品。也许最简单的方法是让机器人在火车轨道上行驶。我希望能够向机器人发送目的地并让它去那里。但是火车轨道和机器人电机并不完美。车轮打滑和不完善的电机意味着机器人不太可能准确地移动到您指定的位置。机器人不止一个,我们需要知道它们都在哪里,这样我们才不会导致它们崩溃。

所以我们添加传感器。也许我们每隔几英尺就在轨道上安装磁铁,并使用霍尔传感器来计算通过了多少磁铁。如果我们数 10 个磁铁,那么机器人应该位于第 10 个磁铁处。当然,有可能漏掉一块磁铁或者数了两次,所以我们必须容忍一定程度的误差。我们可以使用上一节中的代码来跟踪我们的机器人,因为磁铁计数与门口感应非常相似。

但我们还没有完成。我们学会了永远不要丢弃信息。如果您有信息,您应该使用它来改进您的估计。我们漏掉了什么信息?我们知道在每个时刻我们正在向机器人的轮子提供什么控制输入。例如,假设我们每秒向机器人发送一次移动命令 - 向左移动 1 个单位,向右移动 1 个单位,或者静止不动。如果我发送命令“向左移动 1 个单位”,我预计从现在开始的一秒钟内,机器人将比现在所在的位置向左移动 1 个单位。这是一种简化,因为我没有考虑加速度,但我也不想教授控制理论。轮子和马达不完美。机器人最终可能会离开 0.9 个单位,或者 1.2 个单位。

现在整个解决方案很清楚了。我们假设这只狗一直朝他之前移动的任何方向移动。这对我的狗来说是一个可疑的假设!机器人的可预测性要高得多。我们不会根据行为假设做出可疑的预测,而是输入我们发送给机器人的命令!换句话说,当我们调用 predict() 时,我们将传递给机器人的命令运动以及描述该运动可能性的内核。

Simulating the Train Behavior

We need to simulate an imperfect train. When we command it to move it will sometimes make a small mistake, and its sensor will sometimes return the incorrect value.

模拟火车行为

我们需要模拟不完美的火车。当我们命令它移动时,它有时会犯一个小错误,它的传感器有时会返回不正确的值。

class Train(object):

    def __init__(self, track_len, kernel=[1.], sensor_accuracy=.9):
        self.track_len = track_len
        self.pos = 0
        self.kernel = kernel
        self.sensor_accuracy = sensor_accuracy

    def move(self, distance=1):
        """ move in the specified direction
        with some small chance of error"""

        self.pos += distance
        # insert random movement error according to kernel
        r = random.random()
        s = 0
        offset = -(len(self.kernel) - 1) / 2
        for k in self.kernel:
            s += k
            if r <= s:
                break
            offset += 1
        self.pos = int((self.pos + offset) % self.track_len)
        return self.pos

    def sense(self):
        pos = self.pos
         # insert random sensor error
        if random.random() > self.sensor_accuracy:
            if random.random() > 0.5:
                pos += 1
            else:
                pos -= 1
        return pos

With that we are ready to write the filter. We will put it in a function so that we can run it with different assumptions. I will assume that the robot always starts at the beginning of the track. The track is implemented as being 10 units long, but think of it as a track of length, say 10,000, with the magnet pattern repeated every 10 units. A length of 10 makes it easier to plot and inspect.

这样我们就可以编写过滤器了。我们将把它放在一个函数中,这样我们就可以在不同的假设下运行它。我假设机器人总是从轨道的起点开始。轨道实施为 10 个单位长,但可以将其视为长度为 10,000 的轨道,并且每 10 个单位重复一次磁铁图案。长度为 10 使绘图和检查更容易。

def train_filter(iterations, kernel, sensor_accuracy, 
             move_distance, do_print=True):
    track = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    prior = np.array([.9] + [0.01]*9)
    posterior = prior[:]
    normalize(prior)
    
    robot = Train(len(track), kernel, sensor_accuracy)
    for i in range(iterations):
        # move the robot and
        robot.move(distance=move_distance)

        # peform prediction
        prior = predict(posterior, move_distance, kernel)       

        #  and update the filter
        m = robot.sense()
        likelihood = lh_hallway(track, m, sensor_accuracy)
        posterior = update(likelihood, prior)
        index = np.argmax(posterior)

        if do_print:
            print(f'time {i}: pos {robot.pos}, sensed {m}, at position {track[robot.pos]}')
            conf = posterior[index] * 100
            print(f'        estimated position is {index} with confidence {conf:.4f}%:')            

    book_plots.bar_plot(posterior)
    if do_print:
        print()
        print('final position is', robot.pos)
        index = np.argmax(posterior)
        conf = posterior[index]*100
        print(f'Estimated position is {index} with confidence {conf:.4f}')

Read the code and make sure you understand it. Now let’s do a run with no sensor or movement error. If the code is correct it should be able to locate the robot with no error. The output is a bit tedious to read, but if you are at all unsure of how the update/predict cycle works make sure you read through it carefully to solidify your understanding.

阅读代码并确保你理解它。现在让我们在没有传感器或运动错误的情况下运行。如果代码正确,它应该能够准确无误地定位机器人。输出读起来有点乏味,但如果您完全不确定更新/预测周期是如何工作的,请确保您仔细阅读它以巩固您的理解。

import random

random.seed(3)
np.set_printoptions(precision=2, suppress=True, linewidth=60)
train_filter(4, kernel=[1.], sensor_accuracy=.999,
             move_distance=4, do_print=True)
time 0: pos 4, sensed 4, at position 4
        estimated position is 4 with confidence 99.9900%:
time 1: pos 8, sensed 8, at position 8
        estimated position is 8 with confidence 100.0000%:
time 2: pos 2, sensed 2, at position 2
        estimated position is 2 with confidence 100.0000%:
time 3: pos 6, sensed 6, at position 6
        estimated position is 6 with confidence 100.0000%:

final position is 6
Estimated position is 6 with confidence 100.0000

png

We can see that the code was able to perfectly track the robot so we should feel reasonably confident that the code is working. Now let’s see how it fairs with some errors.

我们可以看到代码能够完美地跟踪机器人,因此我们应该有理由相信代码正在运行。现在让我们看看它如何处理一些错误。

random.seed(5)
train_filter(4, kernel=[.1, .8, .1], sensor_accuracy=.9,
         move_distance=4, do_print=True)
time 0: pos 4, sensed 4, at position 4
        estimated position is 4 with confidence 96.0390%:
time 1: pos 8, sensed 9, at position 8
        estimated position is 9 with confidence 52.1180%:
time 2: pos 3, sensed 3, at position 3
        estimated position is 3 with confidence 88.3993%:
time 3: pos 7, sensed 8, at position 7
        estimated position is 8 with confidence 49.3174%:

final position is 7
Estimated position is 8 with confidence 49.3174

png

There was a sensing error at time 1, but we are still quite confident in our position.

Now let’s run a very long simulation and see how the filter responds to errors.

时间 1 出现了感应错误,但我们对自己的位置仍然很有信心。

现在让我们运行一个很长的模拟,看看过滤器如何响应错误。

with figsize(y=5.5):
    for i in range (4):
        random.seed(3)
        plt.subplot(221+i)
        train_filter(148+i, kernel=[.1, .8, .1], 
                     sensor_accuracy=.8,
                     move_distance=4, do_print=False)
        plt.title (f'iteration {148 + i}')


png

We can see that there was a problem on iteration 149 as the confidence degrades. But within a few iterations the filter is able to correct itself and regain confidence in the estimated position.

我们可以看到,随着置信度下降,第 149 次迭代出现了问题。但在几次迭代中,过滤器能够自我纠正并重新获得对估计位置的信心。

Bayes Theorem and the Total Probability Theorem

We developed the math in this chapter merely by reasoning about the information we have at each moment. In the process we discovered Bayes’ Theorem and the Total Probability Theorem.

Bayes theorem tells us how to compute the probability of an event given previous information.

We implemented the update() function with this probability calculation:

p o s t e r i o r = l i k e l i h o o d × p r i o r n o r m a l i z a t i o n   f a c t o r \mathtt{posterior} = \frac{\mathtt{likelihood}\times \mathtt{prior}}{\mathtt{normalization\, factor}} posterior=normalizationfactorlikelihood×prior

We haven’t developed the mathematics to discuss Bayes yet, but this is Bayes’ theorem. Every filter in this book is an expression of Bayes’ theorem. In the next chapter we will develop the mathematics, but in many ways that obscures the simple idea expressed in this equation:

u p d a t e d   k n o w l e d g e = ∥ l i k e l i h o o d   o f   n e w   k n o w l e d g e × p r i o r   k n o w l e d g e ∥ updated\,knowledge = \big\|likelihood\,of\,new\,knowledge\times prior\, knowledge \big\| updatedknowledge= likelihoodofnewknowledge×priorknowledge

where ∥ ⋅ ∥ \| \cdot\| expresses normalizing the term.

We came to this with simple reasoning about a dog walking down a hallway. Yet, as we will see the same equation applies to a universe of filtering problems. We will use this equation in every subsequent chapter.

Likewise, the predict() step computes the total probability of multiple possible events. This is known as the Total Probability Theorem in statistics, and we will also cover this in the next chapter after developing some supporting math.

For now I need you to understand that Bayes’ theorem is a formula to incorporate new information into existing information.

我们仅通过推理我们在每个时刻所拥有的信息来发展本章中的数学。在这个过程中,我们发现了 贝叶斯定理 和 [总概率定理](https://en.wikipedia.org/wiki /Law_of_total_probability)。

贝叶斯定理告诉我们如何计算给定先前信息的事件的概率。

我们用这个概率计算实现了 update() 函数:

p o s t e r i o r = l i k e l i h o o d × p r i o r n o r m a l i z a t i o n   f a c t o r \mathtt{posterior} = \frac{\mathtt{likelihood}\times \mathtt{prior}}{\mathtt{normalization\, factor}} posterior=normalizationfactorlikelihood×prior

我们还没有发展出数学来讨论贝叶斯,但这是贝叶斯定理。本书中的每个过滤器都是贝叶斯定理的表达。在下一章中,我们将发展数学,但在许多方面模糊了这个等式中表达的简单思想:

更新   知识 = ∥ l i k e l i h o o d   o f   n e w   k n o w l e d g e × p r i o r   k n o w l e d g e ∥ 更新\,知识 = \big\|likelihood\,of\,new\,knowledge\times prior\, knowledge \big\| 更新知识= likelihoodofnewknowledge×priorknowledge

其中 ∥ ⋅ ∥ \| \cdot\| 表示规范化术语。

我们通过关于一只狗在走廊上行走的简单推理来得出这个结论。然而,正如我们将看到的,同样的等式适用于过滤问题的领域。我们将在后续的每一章中使用这个等式。

同样,“predict()”步骤计算多个可能事件的总概率。这在统计学中被称为总概率定理,我们还将在开发一些支持性数学之后在下一章中介绍它。

现在我需要你明白贝叶斯定理是一个将新信息合并到现有信息中的公式。

Summary

The code is very short, but the result is impressive! We have implemented a form of a Bayesian filter. We have learned how to start with no information and derive information from noisy sensors. Even though the sensors in this chapter are very noisy (most sensors are more than 80% accurate, for example) we quickly converge on the most likely position for our dog. We have learned how the predict step always degrades our knowledge, but the addition of another measurement, even when it might have noise in it, improves our knowledge, allowing us to converge on the most likely result.

This book is mostly about the Kalman filter. The math it uses is different, but the logic is exactly the same as used in this chapter. It uses Bayesian reasoning to form estimates from a combination of measurements and process models.

If you can understand this chapter you will be able to understand and implement Kalman filters. I cannot stress this enough. If anything is murky, go back and reread this chapter and play with the code. The rest of this book will build on the algorithms that we use here. If you don’t understand why this filter works you will have little success with the rest of the material. However, if you grasp the fundamental insight - multiplying probabilities when we measure, and shifting probabilities when we update leads to a converging solution - then after learning a bit of math you are ready to implement a Kalman filter.

代码很短,但结果令人印象深刻!我们已经实现了一种形式的贝叶斯过滤器。我们已经学会了如何从没有信息开始,并从有噪声的传感器中获取信息。尽管本章中的传感器非常嘈杂(例如,大多数传感器的准确率超过80%),但我们很快就会收敛到狗最可能的位置。我们已经了解到,预测步骤总是会降低我们的知识,但是添加另一个测量,即使它可能有噪声,也会提高我们的知识,使我们能够收敛到最可能的结果上。

这本书主要是关于卡尔曼滤波的。它使用的数学是不同的,但逻辑与本章使用的完全相同。它使用贝叶斯推理从测量和过程模型的组合中形成估计。

如果你能理解本章,你将能够理解和实现卡尔曼滤波器。我怎么强调都不为过。如果有任何不清楚的地方,请返回并重新阅读本章并使用代码。本书的其余部分将以我们在这里使用的算法为基础。如果你不明白为什么这个过滤器会起作用,你将很难成功地使用剩下的材料。然而,如果你掌握了基本的洞察力——当我们测量时乘以概率,当我们更新时移动概率会得到一个收敛的解决方案——那么在学习了一些数学知识之后,你就准备好实现卡尔曼滤波器了。

References

  • [1] D. Fox, W. Burgard, and S. Thrun. “Monte carlo localization: Efficient position estimation for mobile robots.” In Journal of Artifical Intelligence Research, 1999.

http://www.cs.cmu.edu/afs/cs/project/jair/pub/volume11/fox99a-html/jair-localize.html

  • [2] Dieter Fox, et. al. “Bayesian Filters for Location Estimation”. In IEEE Pervasive Computing, September 2003.

http://swarmlab.unimaas.nl/wp-content/uploads/2012/07/fox2003bayesian.pdf

  • [3] Sebastian Thrun. “Artificial Intelligence for Robotics”.

https://www.udacity.com/course/cs373

  • [4] Khan Acadamy. “Introduction to the Convolution”

https://www.khanacademy.org/math/differential-equations/laplace-transform/convolution-integral/v/introduction-to-the-convolution

  • [5] Wikipedia. “Convolution”

http://en.wikipedia.org/wiki/Convolution

  • [6] Wikipedia. “Law of total probability”

    http://en.wikipedia.org/wiki/Law_of_total_probability

  • [7] Wikipedia. “Time Evolution”

https://en.wikipedia.org/wiki/Time_evolution

  • [8] We need to rethink how we teach statistics from the ground up

http://www.statslife.org.uk/opinion/2405-we-need-to-rethink-how-we-teach-statistics-from-the-ground-up

离散卡尔曼滤波Discrete Kalman Filter)是一种常用的状态估计方法,用于从一系列观测数据中估计出系统的状态。在Simulink中,可以使用集成的KF模块来实现离散卡尔曼滤波。这个集成模块功能强大,但可能会比较复杂并且不易修改。 离散卡尔曼滤波的原理是基于系统的动力学模型和观测模型,通过对系统状态的预测和观测值的更新来进行状态估计。它通过估计协方差矩阵来考虑系统的不确定性,并通过最小化估计误差的平方和来优化估计结果。 使用Simulink集成的KF模块可以很方便地进行离散卡尔曼滤波,但如果需要对该模块进行定制化或修改,可能会比较困难。如果你希望实现一个简洁且易于修改的离散卡尔曼滤波模块,可以考虑自己编写代码来实现。在编写代码时,你可以参考已有的博客或文档中提供的初始化参数(如Para_cv1)以及离散卡尔曼滤波的原理。这样你就能够根据自己的需求进行修改和定制化了。 另外,除了离散卡尔曼滤波,还有其他一些滤波器可以用于系统状态估计,如连续时间卡尔曼滤波、混合时间卡尔曼滤波和无味卡尔曼滤波器(Unscented Kalman Filter,UKF)。每种滤波器都有其特点和适用的场景,你可以根据具体的需求选择合适的滤波器来进行状态估计。<span class="em">1</span><span class="em">2</span><span class="em">3</span> #### 引用[.reference_title] - *1* [离散时间卡尔曼滤波/SIMULINK](https://download.csdn.net/download/weixin_44044161/12609945)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_1"}}] [.reference_item style="max-width: 33.333333333333336%"] - *2* [连续/离散/混合卡尔曼滤波SIMULINK初始化](https://download.csdn.net/download/weixin_44044161/12609958)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_1"}}] [.reference_item style="max-width: 33.333333333333336%"] - *3* [卡尔曼滤波-Simulink仿真](https://blog.csdn.net/weixin_45696231/article/details/123972759)[target="_blank" data-report-click={"spm":"1018.2226.3001.9630","extra":{"utm_source":"vip_chatgpt_common_search_pc_result","utm_medium":"distribute.pc_search_result.none-task-cask-2~all~insert_cask~default-1-null.142^v93^chatsearchT3_1"}}] [.reference_item style="max-width: 33.333333333333336%"] [ .reference_list ]
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值