I think there is a mistake in The formula of Total Probability Theorem in 02-Discrete-Bayes.ipynb.

P(X^{t−1}_{i}) should be P(X^{t−1}_{j}).

When executing

#format the book
%matplotlib inline
%load_ext autoreload
%autoreload 2  
from __future__ import division, print_function
import sys
sys.path.insert(0,'./code')
from book_format import load_style
load_style()

from any page, I get the following warning:

D:\Chad\Documents\WinPython-64bit-3.4.3.6\python-3.4.3.amd64\lib\site-packages\matplotlib\__init__.py:876: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

This Python distro is using Matplotlib 1.5.0rc3. I have confirmed that changing the axes.color_cycle line in 538.json to

  "axes.prop_cycle": "cycler('color', ['#6d904f','#013afe', '#202020','#fc4f30','#e5ae38','#A60628','#30a2da','#008080','#7A68A6','#CF4457','#188487','#E24A33'])",

fixes the warning. But I am unsure how to make the json file conditional based on Matplotlib version, so as to be backwards compatible with older distributions.

Recently I'm bump into your comment in the filterpy implementation of sigma point calculation:

If your method returns a triangular matrix it must be upper
triangular. Do not use numpy.linalg.cholesky - for historical
reasons it returns a lower triangular matrix.

In my opinion it will be great add note about requirement upper triangular matrix instead of lower after words:

SciPy provides cholesky() method in scipy.linalg.
If your language of choice is Fortran, C, or C++,
libraries such as LAPACK provide this routine.
Matlab provides chol().

P.S.

Can you explain where did this requirements from?

I've tried use both of matrix view. When I use lower triangular matrix the covariance matrix P is growing up (after each step my confidence decreases).

Is it expected behaviour?

I try to found out how to find upper triangular matrix using Cholesky decomposition. The most of sources say that Cholesky decomposition is:

A = L * L'

The few sources say about another form:

A = U' * U

So it's looks like U = L', am I right? Does the upper triangular root matrix just transposed lower triangular root matrix?

Here's a couple typos I noticed while reading chapter 1:

• "Now, let's assume we we gained weight."

  "we" is duplicated.

• "For example, for the 100 kg weight weight our estimate might be 99.327 kg due to sensor errors."

  "weight" is duplicated.

Hi,

I am trying to work with this ipython files, but they dont work.
It appears this error, after calling any ipython file, any advice?
thanks in advance, Jaime

cd Kalman-and-Bayesian-Filters-in-Python/
jaimebayes@jaimebayes-OptiPlex-755:~/Kalman-and-Bayesian-Filters-in-Python$ ipython notebook
2015-06-02 13:34:11.222 [NotebookApp] Using existing profile dir: u'/home/jaimebayes/.ipython/profile_default'
2015-06-02 13:34:11.320 [NotebookApp] Using MathJax from CDN: https://cdn.mathjax.org/mathjax/latest/MathJax.js
2015-06-02 13:34:11.519 [NotebookApp] Serving notebooks from local directory: /home/jaimebayes/Kalman-and-Bayesian-Filters-in-Python
2015-06-02 13:34:11.520 [NotebookApp] 0 active kernels
2015-06-02 13:34:11.520 [NotebookApp] The IPython Notebook is running at: http://localhost:8888/
2015-06-02 13:34:11.520 [NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation).

(process:4174): GLib-CRITICAL **: g_slice_set_config: assertion 'sys_page_size == 0' failed
2015-06-02 13:34:27.679 [NotebookApp] WARNING | Unreadable Notebook: /home/jaimebayes/Kalman-and-Bayesian-Filters-in-Python/00_Preface.ipynb global name 'NBFormatError' is not defined
WARNING:tornado.access:400 GET /api/notebooks/00_Preface.ipynb?_=1433270067150 (127.0.0.1) 26.97ms referer=http://localhost:8888/notebooks/00_Preface.ipynb

A reader stated that they had difficulty using Python 2, and had to switch to 3. It is time to test 2.7, 3.3, 3.4, 3.5.

In http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/05_Kalman_Filters/Kalman_Filters.ipynb

Hi,
I am rather enjoying going though your book. (Currently working though and translating bits of it into julia for my understanding of both Julia and Kalman filters.)
I appreciate it is a work in progress.

at In[89]:

    def update(self, z):
        self.x = (self.P * z + self.x * self.R) / (self.P + self.R)
        self.P = 1. / (1./self.P + 1./self.R)

My instincts are telling me that this should be:

    def update(self, z):
        self.x = (self.R * z + self.x * self.P) / (self.P + self.R)
        self.P = 1. / (1./self.P + 1./self.R)

Because P is the variance (uncertainty) of the state,
and R is the variance of the measurement.

I'm not to certain though.
If I am wrong, perhaps the section could be written to make it more clear why this should be the case.

At the very end of the Tracking and Control section of chapter 2, the last sentence appears to be truncated:

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

(Maybe the likelihood of distance/destination?)

The license on the readme (updated Mar 9, 2015) is CC BY, but the license called out in license.html (May 16, 2014) is CC BY-NC-SA. Legally, the BY takes precedence over the BY-NC-SA and is irrevocable, but they should be the same to avoid confusion.

Hi Roger.
After a long pause on the topic I came back to download the latest version and noticed the book PDF is no longer present in the zip (and in the repository).

edit: I just read in one of the recent commits' comment that you actually removed the PDF.
I'm often offline and was relying on the PDF to continue reading...

I get a good number of errors when trying to build the PDF on my laptop... Would you please consider adding back the PDF to the repo?

Kind regards,
Mario

It might be an operating system/browser problem, but I cannot see the equations properly for chapter 10 (I haven't tried other chapters). Attached is a screenshot of what I see using Firefox Mozilla v. 42.0 on Ubuntu 12.04 LTS.

Hi Roger,
first, it's really great that didactically skilled people like you share their knowledge with dumb people like me in such a nice way. I really learned a lot from your book!

I found some issues I want to address. I don't know if this is the right place to do, but I interpreted from your preface that it is.
Anyway: I implemented my own Kalman filter and it looked quite ugly, so I took the input parameters that you used in chaper 6 when implementing the first full KF. It still didn't look the same and I found that you used for the state transition matrix
F = 1 1
0 1
The top right index is supposed to be the step, right? So, here you use dt=1 and further down for the calculation of the Q matrix in this example you use dt=0.1.

Is that a bug or am I missing something? I haven't read past chapter 6 yet, but to my feeling this is wrong, or? After I changed in my version the state transition to your version, the KF performed like yours.

Cheers
Marko

In Chapter 7, iterative least squares for sensor fusion section there is inconsistent wording and use of the - and + superscripts. This occurs to a lesser extent in the rest of the book.

PDF is floating around the web with no obvious link back to the github account, and the creative common license is buried in the preface. I'd like a page directly after the title page giving the license and the links to github.

Hi there,

In binder if I try to compile this listing:

import book_plots
from book_plots import interactive_plot
import gh_internal as gh
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 # day
scale_factor = 4/10

def predict_using_gain_guess(weight, gain_rate, do_print=True, sim_rate=0): 
    # store the filtered results
    estimates, predictions = [weight], []

    # most filter literature uses 'z' for measurements
    for z in weights: 
        # predict new position
        prediction = weight + gain_rate * time_step

        # update filter 
        weight = prediction + scale_factor * (z - prediction)

        # save
        estimates.append(weight)
        predictions.append(prediction)
        if do_print:
            gh.print_results(estimates, prediction, weight)

    # plot results
    gh.plot_gh_results(weights, estimates, predictions, sim_rate)

initial_guess = 160.
with interactive_plot():
    predict_using_gain_guess(weight=initial_guess, gain_rate=1)   

The following exception is thrown:

ImportError                               Traceback (most recent call last)
<ipython-input-1-1dc46de6cfa4> in <module>()
----> 1 import book_plots
      2 from book_plots import interactive_plot
      3 import gh_internal as gh
      4 import matplotlib.pyplot as plt
      5 

ImportError: No module named 'book_plots'

If I change the import from import book_plots to import code.book_plots it seems to behave better but then the exception is thrown about filterpy:

/home/main/notebooks/book_format.py in test_filterpy_version()
     45 def test_filterpy_version():
     46 
---> 47     import filterpy
     48     from distutils.version import LooseVersion
     49 

ImportError: No module named 'filterpy'

I am noticing that in Firefox Ubuntu the div.inner_cell is a few pixels too short, cutting off the end of the row of words and is visually distracting. I notice that other notebooks (from other books) don't have this issue, even though they seem to share the same .css files! But, your font is different than the linked example, so something must be different when you compile it. I can't tell from a first glance of the HTML what is different though.

Looks like a great book otherwise. Gotta tear into it :)

In 06-Multivariate-Kalman-Filters.ipynb you write

Your job as a designer will be to design the state (x,P), the process (x,P), the measurement (z,R), and the unexplained measurement function H. If the system has control inputs, such as a robot, you will also design B and u.

I think you mean .... the process(F,Q) .... ?

hi,

Is it suitable to use Kalman filter for predicting spikes ?
F.e. let say I monitor CPU usage where the values go from 0 to 100 ... will the Kalman filter be "agile" enough to detect and predict when the usage go above say 80% ... with say better tha 60-70% hit rate.

Or you would use some other technique.

Early in the book I use 'state transtition function' everywhere. Which is a valid term, but it is not so clear maybe that this is just the process model. I do use 'process noise', so why not make the connection more obvious, especially since i mostly use process model later in the book.

In /06-Multivariate-Kalman-Filters.ipynb you are explaining ndarray and state

However, NumPy recognizes 1D arrays as vectors, so I can simplify this line to use a 1D array.
x = np.array([10.0, 4.5])
x

Since before you have explained how X should be a nx1 matrix having said that

x = np.array([[10.0, 4.5]]).T

I assume you want to say that you could write x = np.array([10.0, 4.5]).T (i.e. with transpose) to mean the same thing as double brackets, transpose.

I wrote to you on hackernews (but know I'll forget to check answer, so I copy/past here)

Sorry about the notebook format change. Did you loose work? Conversion should be handled for you, if not it is a bug, we can fix it on 3.1.
We would be happy to get more feedback on your writing process and your need, feel free to directly contact the team (IPython-dev at scipy.org, or issue on main IPython repo is fine).
As for concept of "book" or collection of notebook, we are working on that (integration with sphinx)

If you loose things in migration it's not normal, we can fix that, please open an issues. We would also love to get more details feedback of your pain points in writing. We are having our bi-annual IPython dev meeting this week so it is the time where we will sketch the future for the next 6 to 12 month on IPython now.

I see you are in SF Bay Area, we are meeting in UC Berkeley, I guess you could even come and say Hi, if you prefer in person face-to-face feedback/complaint :-)

Chapter 2

• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/02-Discrete-Bayes.ipynb#Terminology
  In the 4th paragraph of the section, "how dog's behave" should be "how dogs behave".
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/02-Discrete-Bayes.ipynb#The-Discrete-Bayes-Algorithm
  In the update step 3, "Determine whether whether the measurement" should be "Determine whether the measurement"

Chapter 3

• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/03-Gaussians.ipynb#Nomenclature
  In 2nd to last paragraph of the section, "forcing us to use computationally expensive alternative" should be "forcing us to use a computationally expensive alternative" or "forcing us to use computationally expensive alternatives" or something.

Chapter 4

• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/04-One-Dimensional-Kalman-Filters.ipynb#Understanding-Gaussian-Multiplication
  After two code samples, "because I know have two readings" should be "because I now have two readings".
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/04-One-Dimensional-Kalman-Filters.ipynb#Exercise:-Modify-Variance-Values
  The second sentence in the section, "Which has a larger effect on the variance convergence?." has a period after a question mark when it shouldn't.
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/04-One-Dimensional-Kalman-Filters.ipynb#KF-Animation
  In the last paragraph of the section, "how much it changed for the firs noisy measurement" should be "how much it changed for the first noisy measurement".
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/04-One-Dimensional-Kalman-Filters.ipynb#Introduction-to-Designing-a-Filter
  In the first paragraph, "lets get a feel for it" should be "let's get a feel for it".
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/04-One-Dimensional-Kalman-Filters.ipynb#Example:-Extreme-Amounts-of-Noise
  In the first paragraph, "standard deviation of 300 meters?." has a period after a question mark when it shouldn't.

Chapter 5

• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/05-Multivariate-Gaussians.ipynb#Multivariate-Normal-Distribution-Equation
  After the first code sample, "probability density of for any value" has "of" and "for".
• http://nbviewer.jupyter.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/05-Multivariate-Gaussians.ipynb#Multivariate-Normal-Distribution-Equation
  In the next paragraph, "If we think of this as a the Gaussian" has "a" and "the".

The chapter based on example when dog Simon tracking in hallway. The hallway is described as one-dimensional array with length 10.

To keep the problem small, we will assume that there are only 10 positions in a single hallway to consider, which we will number 0 to 9, where 1 is to the right of 0, 2 is to the right of 1, and so on.

And all sections use this data domain. But two sections use different data: Incorporating Movement Data and Adding Noise to the Prediction.

In the Incorporating Movement Data section is used array with length 4 (but it is still talk about the Simon in the hallway).

At the beginning of the Adding Noise to the Prediction section is used array with length 8, but later used common data domain - array with length 10.

This is not error! All charts matched with corresponding code in Python and description, but, in my opinion this is some kind of drawback.

I have a few suggestion to improve problem domain description of para "Tracking an Airplane" of "Unscented Kalman Filter" chapter. First of all, this improvements not related with the main topic of your book, so them is not important. Also improving of problem domain description can pollute the main idea. So you should decide accept some of them or reject, but anyway I'd like notice about them.

You wrote:

We will track one dimension on the ground and the altitude of the aircraft.

But later you wrote:

By timing how long it takes for the reflected signal to return it can compute the slant distance and bearing to the target.

We can’t compute bearing to the target by timing.

The angular determination of the target is determined by the directivity of the antenna.
link

If we measure altitude then more appropriate use the "elevation angle" instead of "bearing".

The elevation angle is the angle between the horizontal plane and the line of sight, measured in the vertical plane.
link

Also such radars called height finders.

A height finder is a ground based aircraft altitude measuring device.
link

You wrote:

Radars also provide velocity measurements via the Doppler effect, but we will be ignoring this complication for the moment.

Velocity of flying object can be calculated using radar measurement, but strictly speaking radar can't measure velocity, instead radar can measure radial velocity.

Only the radial component of the speed is relevant. When the reflector is moving at right angle to the radar beam, it has no relative velocity. Vehicles and weather moving parallel to the radar beam produce the maximum Doppler frequency shift.
link

You wrote:

We compute the (x,y) position of the aircraft from the slant distance and bearing as illustrated by this diagram:

This is very subjective judgement but in my opinion using (x,z) instead of (x,y) will be better (by the way, in code snippet 15 you use dz as delta of altitude).

In your chart ekf_internal.show_radar_chart:
You used theta as elevation angle usually theta used for bearing but epsilon used for elevation angle.

You wrote:

A typical radar might update only once every 12 seconds so we will use that for our epoch period.

As far as I'm concerned this is true if we talk about surveillance radars, but for height finder this time usually smaller (about 3 seconds).

I hope that this will be helpful for You!

I have a question about Chapter 10 "Unscented Kalman Filter".

In paragraph "Implementation of the UKF: Sigma Points" in the sigma point equation under square root is placed lambda. Later (in the text) instead of lambda used kappa. Later at the source code snippet:

U = scipy.linalg.cholesky((n+kappa)*P)

used kappa too.

In paragraph "Choosing the Sigma Parameters" used lambda again.

But in the description of the filter by link used kappa.

So my question is: what is the correct equation for the sigma points? Is this error or I have missed some thing?

We now only need 1d arrays. code is corrected, but (maybe) not all the text.

I have a question about Chapter 10 "Unscented Kalman Filter".

In paragraph "Implementation of the UKF: Choosing Sigma Points" specified that sum of weights of means and covariances must be equal one.

Later In paragraph "Implementation of the UKF: Van der Merwe's Scaled Sigma Point Algorithm: Weight Computation" specified the equation for weights computation.

When I try to compute weights using the algorithm with suggested parameters (for two dimensions):

• alpha = 0 ... 1;
• beta = 2;
• kappa = 1.

I will get residual with constraint for weights of covariances. For different values of alpha the sum of weights lies between 3 and 4.

Changing number of dimensions or kappa has no effect. I have found out that to reduce the residual beta must be equal -2.

So my questions:

• do you have the similar results of computations (sum of weights of covariances)?
• if yes, should I always respect the constraints?
• If yes, should I normalize weights myself?
• If no, maybe instead adding beta should be used subtraction?

I'm getting questions that show me that I have not talked about tracking situations and how the order of the filter affects the performance. Yes, this belongs in a IMM chapter, but the regular chapters can use this info as well in shorter form.

Hi,

Just tried to use binder to read your book. (Nice work, well done).

But when trying to run some code, such as the first cell in 01-g-h-filter.ipynb I get the following error which implies that filerpy has not been installed, though I note that filerpy is listed in requirements.txt.

Here is the error message I get.

ImportError Traceback (most recent call last)
in ()
2 get_ipython().magic('matplotlib inline')
3 from future import division, print_function
----> 4 from book_format import load_style
5 load_style()

/home/main/notebooks/book_format.py in ()
59 # called when this module is imported at the top of each book
60 # chapter so the reader can see that they need to update FilterPy.
---> 61 test_filterpy_version()
62
63

/home/main/notebooks/book_format.py in test_filterpy_version()
39
40 def test_filterpy_version():
---> 41 import filterpy
42 min_version = [0,0,28]
43 v = filterpy.version

ImportError: No module named 'filterpy'

Hi,

Nice work on the book and on providing a python implementation.

I recently came across a paper that discussed estimating the measurement noise covariance matrix in an online fashion as follows:

                 \\(\hat{R(k)} = C_v + H \bar{P}(k) H^T\\) 

where

            \\(\bar{P}(k) \\) is the posteriori error covariance matrix and \\(C_v\\) , at time \\(k\\), is computed from an average moving window of size \\(m\\), as

                      \\(C_v (k) =1/ m  ∑ v(k − 1)v(k − 1)^T\\)

The thing is when I implement this, I find my estimates are really smooth but it takes a long time to reach convergence.

Wondering if you have any idea why this is so. The paper in question is:

J. Wang, “Stochastic modeling for real-time kinematic gps/glonass position,” Navigation, vol. 46, pp. 297–305, 2000.

Looking forward to your reply.

Thank you!

I just cloned the repo and ran jypyter notebook from the root directory. Now I see /tree in my browser and double click 01-g-h-filter.ipynb. However, when I try to run

import gh_internal as gh
gh.plot_errorbar1()

I get the error below

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-7-dd29f703c030> in <module>()
----> 1 import gh_internal as gh
      2 gh.plot_errorbar1()

ImportError: No module named 'gh_internal'

I'm sure this isn't an issue with the code but my lack of understanding for how the gh_internal.py is being imported from the code directory to the notebook.

Also, I notice the same error on mybinder as well: http://app.mybinder.org/1459191046/notebooks/01-g-h-filter.ipynb

How should this be done?

I sort of have this with the notation, but it needs to be expanded. There are so many different terms in the literature for the same thing - dynamic matrix, state transition matrix, it goes on and on. It is just confusing for everyone.

I'm thinking a table with all the different terms for a concept grouped together, along with the typical math symbols for the corresponding matrix/vector, if applicable.

This is not an issue but a question.

In chapter 6 "Multivariate Kalman Filters" you mention the measurement space. Mostly the reason of using the measurement space is clear, but I have a several questions.

1. First of all I'd like to know how called the opposite of the measurement space? I guess that it might be called the state space, am I right?

2. Some members (Residual, Measurement and Covariance) of equations are explicit belong to the measurement space, but other not. So my question is: which members are belong to the measurement space and which members are not:

   • State (obviously no).
   • Covariance (yes);
   • State transition function (looks like no);
   • Process noise (looks like yes);
   • Control input (looks like no);
   • Control function (looks like no);
   • Residual (yes);
   • Measurement (yes);
   • Measurement function (?);
   • Measurement noise (looks like yes);
   • System uncertainty (looks like yes);
   • Kalman gain (?!).

   Or other words: to respect your example about measurement of temperature, what units (Volt, Celsius, Volt/Celsius, etc) must have each member?

3. This question is consequence of previous question. The most obscure for me how we can update State by adding previous state estimation (which is not in the measurement space) and Residual (which is in the measurement space)? Does Kalman gain performs transformation of Residual?

Thanks!

P.S. Sorry if answers on this questions already are in your book but I have not found their.

A lot of people are interested in this book due to the robotics/UAV use of the filters. It is hard to do so when there are no examples of how to account for the control inputs.

Chapter five, near the bottom: "the aircraft's state includes things such as as heading, engine RPM, weight, color, the first name of the pilot, and so on. We cannot sense these directly using the position sensor so they are not observed."

It seems to me like the aircraft's heading would be an observable, hidden variable like velocity; you could estimate it, and use it to improve your results, the same way you did with velocity.

According to description of the example weight have to gain of 1 lb per day.

I hand generated the weight data to correspond to a true starting weight of 160 lbs, and a weight gain of 1 lb per day. In other words one day one the true weight is 160lbs, on the second day the true weight is 161 lbs, and so on.

It means that when weight equal 160 lb at 1-st day it will equal 171 lb at 12-th day:
1 - 160
2 - 161
3 - 162
4 - 163
5 - 164
6 - 165
7 - 166
8 - 167
9 - 168
10 - 169
11 - 170
12 - 171

But on the chart weight at 12-th day equals 172 lb.

It looks like off-by-one error.

In my opinion correct values of weight should be:
0 - 160 (today)
1 - 161 (tomorrow)
2 - 162
3 - 163
4 - 164
5 - 165
6 - 166
7 - 167
8 - 168
9 - 169
10 - 170
11 - 171
12 - 172

In chapter 3, under the heading "Product of Gaussians", there's the formula for the product's mean. I think the second term in the numerator uses the wrong variance. The formula is also in chapter 4, and it looks correct there.

I put a try/except block around code in chapter 5 because it depends on code in FilterPy that is not released yet. Remove block once FilterPy is updated.

Hi (again),
In http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/09_Extended_Kalman_Filters/Extended_Kalman_Filters.ipynb#Designing-Q
We get some details on how to find Q,
It would be great if this concluded by substiting in values for $w_alt$ and $w_vel$, to produce the final $Q$.

To may knowledge it goes:
$w_{vel}^2=5$, $v_{vel}^2=10$, and $w_{alt} w_{vel}=0$ (as they are independent and thus uncorrelated)

This gives a final matrix of:

Q=
[ 0.000208333  0.00625  0.0
 0.00625      0.25     0.0
 0.0          0.0      0.5]

However that produces awful results, so i think I must have got the math wrong.
It would be great to have this section more complete so I can better understand my mistake.

PS: Is it ok if I raise issues like this? If you would prefer not to receive such feedback at this stage, I can be quiet (and come back when you are, if you would like. :-) )

I get questions about papers - how would I implement this paper in your code. I think there would be a lot of pedagogical value in taking a few open source papers, working through them, and implementing them in FilterPy.

With my recent edits the g-h chapter is pretty close to the KF chapter. The KF chapter intro should remind us of what was done in the g-h chapter to lay the land for the rest of the chapter.

I hate how it shows 3 scatter points in the legend, clamp it to one.

The formula for total probability theorem is wrong here. Probability of being at position i at time t,
P(x_i^t) is given by the sum of probabilities of being at position j at time t-1 multiplied by probability of moving from position j to i at time t. Therefore, it will be,
P(x_i^t) = for all j, sum of ( P(x_j^t-1) * P(x_i|x_j) )

While reading I stumbled upon the following minor issues. I didn't want to open a new issue for each of them, so here they are:

• 06-Multivariate-Kalman-Filters
  • [..]So, if K is large then (1−cK) is small, and P will be made smaller than it was. If K is small, then (1−cK) is large, and P will be made larger than it was. [..]
    • Doesn't P get smaller in both cases, just smaller in the first one?
• 07-Kalman-Filter-Math
  • [..] division is scalar's analogous operation for matrix inversion [..]
    • Reciprocal is the analogous operation, it's only somewhat division here because of the subsequent multiplication
  • integration routines such as scipy.integrate.ode. These routines are robust, but [end of paragraph]
• 08-Designing-Kalman-Filters
  • [..] a Kalman filter with R=0.5 as a thin blue line, and a Kalman filter with R=10 as a thick red line [..]
    • On binder at least R=0.5 is a thin green line, and R=10 is a thick blue line
  • Tracking Noisy Data
    author's note: not sure if I want to cover this. If I do, not sure I want to cover this here.
    • Ok, not really a bug. I just hope you do because I would really love that chapter :)

From emailed report:

In EKF:
def residual(a,b):
""" compute residual between two measurement containing [range, bearing].
Bearing is normalized to [0, 360)"""

Testing looks like [-180, 180]
print('a:', np.rad2deg(residual(np.array([[0], [np.deg2rad(195.0)]]),
np.array([[1], [np.deg2rad(355.0)]]))[1]))
print('b:', np.rad2deg(residual(np.array([[0], [np.deg2rad(355.0)]]),
np.array([[1], [np.deg2rad(5.0)]]))[1]))

Most of the simulation have no process noise, so their behavior is a bit too ideal. Add a bit of noise to the vel/acc.

The Non-linear filter notebook opens with the example,

There can be nonlinearity in the process model. Suppose we wanted to track the motion of a weight on a spring, such as an automobile's suspension. The equation for the motion with $m$ being the mass, $k$ the spring constant, and $c$ the damping force is

$m\frac{d^2x}{dt^2} + c\frac{dx}{dt} +kx = 0$

There is no linear solution for $x(t)$ for this second order differential equation, and therefore we cannot design a Kalman filter using the theory that we have learned.

which is a linear system described by the state variables x, $\frac{dx}{dt}$.

Hello,

As I am still a beginner when it comes to kalman filters, I was wondering if someone could provide me with help on how to solve my problem.
At the moment I am trying to estimate the position, velocity and acceleration of an e-puck robot. The robot comes equipped with encoders for the motors' position and an accelerometer giving a 3D acceleration vector in m/s^2.
Here is how I initialized my Unscented Kalman Filter:

    points = sigPts(9, alpha=0.1, beta=2.0, kappa=0.0)
    self._ukf = ukf(dim_x=9, dim_z=5, dt=TIME_STEP, hx=self.measurements, fx=EpuckNode.transition, points=points)
    self._ukf.x = np.array([float(init_pos_x), 0.0, 0.0, float(init_pos_y), 0.0, 0.0, float(init_pos_z), 0.0, 0.0])
    self._ukf.P = np.identity(9, dtype=float)
    self._ukf.Q[0:3, 0:3] = Q_discrete_white_noise(3, dt=TIME_STEP, var=0.02)
    self._ukf.Q[3:6, 3:6] = Q_discrete_white_noise(3, dt=TIME_STEP, var=0.02)
    self._ukf.Q[6:9, 6:9] = Q_discrete_white_noise(3, dt=TIME_STEP, var=0.02)
    self._ukf.R = np.diag([0.2 ** 2, 0.2 ** 2, 0.2 ** 2, 0.03 ** 2, 0.03 ** 2])  # Sensors' sensitivity

This is the transition function:

def transition(state, dt):
    # Define the transition matrix
    f = np.array([
        [1, dt, (dt ** 2)/2, 0, 0, 0, 0, 0, 0],
        [0, 1, dt, 0, 0, 0, 0, 0, 0],
        [0, 0, 1, 0, 0, 0, 0, 0, 0],
        [0, 0, 0, 1, dt, (dt ** 2)/2, 0, 0, 0],
        [0, 0, 0, 0, 1, dt, 0, 0, 0],
        [0, 0, 0, 0, 0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, dt, (dt ** 2)/2],
        [0, 0, 0, 0, 0, 0, 0, 1, dt],
        [0, 0, 0, 0, 0, 0, 0, 0, 1]
    ], dtype=float)
   # Return the array containing all the transforms
    return np.dot(f, state)

And finally the measurement function:

def measurements(self, measures):
        """
        Define the computation required for the measurements
        :param measures: An array containing: [encoder_L, encoder_R, accel_X, accel_Y]
        :return: 
        """

        # Compute the position and heading
        left_encoder = measures[0]
        right_encoder = measures[3]

        # Compute the number of steps each wheel executed
        left_steps_diff = left_encoder * MOT_STEP_DIST - self._leftStepsPrev    # Expressed in meters.
        right_steps_diff = right_encoder * MOT_STEP_DIST - self._rightStepsPrev  # Expressed in meters.

        # Compute the rotation and step differences
        delta_theta = (right_steps_diff - left_steps_diff) / WHEEL_DISTANCE # Expressed in radiant.
        delta_steps = (right_steps_diff + left_steps_diff) / 2  # Expressed in meters.

        # Extract the robot's position and orientation
        pos_x = self._pos_x + delta_steps*math.cos(self._pos_z + delta_theta/2)  # Expressed in meters.
        pos_y = self._pos_y + delta_steps*math.sin(self._pos_z + delta_theta/2)  # Expressed in meters.
        pos_z = self._pos_z + delta_theta   # Expressed in radiant.

        # Filter the measurements obtained from the accelerometer and the encoders

        # Keep track of the number of steps for both wheels
        self._leftStepsPrev = left_encoder * MOT_STEP_DIST  # Expressed in meters.
        self._rightStepsPrev = right_encoder * MOT_STEP_DIST    # Expressed in meters.

        # Return the array containing all the measurements
        return np.array([pos_x, pos_y, pos_z, measures[2], measures[5]])

However each time I try to execute the filter through the following test code:

# Read the sensors' value
left_encoder = self.getLeftEncoder()
right_encoder = self.getRightEncoder()
accels = self._accelerometer.getValues()

# Predict the next state
self._ukf.predict()
self._ukf.update([left_encoder, right_encoder, accels[0], accels[1]])

# Update the robot's state
self._pos_x = self._ukf.x[0]
self._speed_x = self._ukf.x[1]
self._accel_x = self._ukf.x[2]
self._pos_y = self._ukf.x[3]
self._speed_y = self._ukf.x[4]
self._accel_y = self._ukf.x[5]
self._pos_z = self._ukf.x[6]
self._speed_z = self._ukf.x[7]
self._accel_z = self._ukf.x[8]

print(self._ukf.x)

I get an error:

Traceback (most recent call last):
self._ukf.update([left_encoder, right_encoder, accels[0], accels[1]])
File "/usr/local/lib/python2.7/dist-packages/filterpy/kalman/UKF.py", line 341, in update
y = self.residual_z(z, zp) #residual
ValueError: operands could not be broadcast together with shapes (4,) (5,) 

Could you please explain or at least point me in the direction of my mistake?
Thank you very much for your help.

Sincerely.

I'm using the latest version of Anaconda (still python 2.7) along with Chrome under Ubuntu but I am not able to load any of the pages when I click on the notebooks. I get an error message "Bad request". Any ideas on how to get around this?

It takes them many minutes to slow down or speed up significantly. So, if I know that the train is at kilometer marker 23km at time t and moving at 60 kph, I can be extremely confident in predicting its position at time t + 1 second. And why is that important? Suppose we can only measure its position with an accuracy of 500 meters. So at t+1 sec the measurement could be anywhere from 22.5 km to 23 km. But the train is moving at 60 kph, which is 16.6 meters/second. So if the next measurement says the position is at 23.4 we know that must be wrong. Even if at time t the engineer slammed on the brakes the train will still be very close to 23.0166 km because a train cannot slow down very much in 1 second.

I think it should be
"So at t+1 sec the measurement could be anywhere from 22.516 km to 23.484 km. But the train is moving at 60 kph, which is 16.6 meters/second. "
since the statement 'with an accuracy of 500 meters' means both + and - ways.
so the next 'So if the next measurement says the position is at 23.4 we know that must be wrong. ' also need a fix!