Added Ensemble Kalman Filter chapter.
Had to rename and move directories around to allow the new chapter to fit in.
This commit is contained in:
1 parent
3c6abb8899
commit
a27ebac336
6 files changed
+955
-22
No files matched your search
@@ -0,0 +1,514 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:89143715003512c2e724530f782eafab4fe0ced8e56da4ce0c97b2064b321ad3"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "heading",
|
||||
"level": 1,
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Ensemble Kalman Filters"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"#format the book\n",
|
||||
"%matplotlib inline\n",
|
||||
"from __future__ import division, print_function\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0,'../code') # allow us to format the book\n",
|
||||
"import book_format\n",
|
||||
"book_format.load_style()"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"html": [
|
||||
"<style>\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
|
||||
"\n",
|
||||
" div.cell{\n",
|
||||
" width: 850px;\n",
|
||||
" margin-left: 0% !important;\n",
|
||||
" margin-right: auto;\n",
|
||||
" }\n",
|
||||
" div.text_cell code {\n",
|
||||
" background: transparent;\n",
|
||||
" color: #000000;\n",
|
||||
" font-weight: 600;\n",
|
||||
" font-size: 11pt;\n",
|
||||
" font-style: bold;\n",
|
||||
" font-family: 'Source Code Pro', Consolas, monocco, monospace;\n",
|
||||
" }\n",
|
||||
" h1 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
"\t}\n",
|
||||
"\t\n",
|
||||
" div.input_area {\n",
|
||||
" background: #F6F6F9;\n",
|
||||
" border: 1px solid #586e75;\n",
|
||||
" }\n",
|
||||
"\n",
|
||||
" .text_cell_render h1 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 30pt;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#c76c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 1em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: wrap;\n",
|
||||
" } \n",
|
||||
" h2 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h2 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 20pt;\n",
|
||||
" font-style: italic;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#c76c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" } \n",
|
||||
" h3 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h3 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-size: 18pt;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#d77c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 2em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" h4 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h4 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" color:#d77c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 0.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" h5 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h5 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-style: normal;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" margin-bottom: 0em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" div.text_cell_render{\n",
|
||||
" font-family: 'Arimo',verdana,arial,sans-serif;\n",
|
||||
" line-height: 135%;\n",
|
||||
" font-size: 125%;\n",
|
||||
" width:750px;\n",
|
||||
" margin-left:auto;\n",
|
||||
" margin-right:auto;\n",
|
||||
" text-align:justify;\n",
|
||||
" text-justify:inter-word;\n",
|
||||
" }\n",
|
||||
" div.output_subarea.output_text.output_pyout {\n",
|
||||
" overflow-x: auto;\n",
|
||||
" overflow-y: scroll;\n",
|
||||
" max-height: 50000px;\n",
|
||||
" }\n",
|
||||
" div.output_subarea.output_stream.output_stdout.output_text {\n",
|
||||
" overflow-x: auto;\n",
|
||||
" overflow-y: scroll;\n",
|
||||
" max-height: 50000px;\n",
|
||||
" }\n",
|
||||
" code{\n",
|
||||
" font-size: 70%;\n",
|
||||
" }\n",
|
||||
" .rendered_html code{\n",
|
||||
" background-color: transparent;\n",
|
||||
" }\n",
|
||||
" ul{\n",
|
||||
" margin: 2em;\n",
|
||||
" }\n",
|
||||
" ul li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.5em; \n",
|
||||
" }\n",
|
||||
" ul li li{\n",
|
||||
" padding-left: 0.2em; \n",
|
||||
" margin-bottom: 0.2em; \n",
|
||||
" margin-top: 0.2em; \n",
|
||||
" }\n",
|
||||
" ol{\n",
|
||||
" margin: 2em;\n",
|
||||
" }\n",
|
||||
" ol li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.5em; \n",
|
||||
" }\n",
|
||||
" ul li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.2em; \n",
|
||||
" }\n",
|
||||
" a:link{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" a:visited{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" }\n",
|
||||
" a:hover{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" }\n",
|
||||
" a:focus{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" a:active{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" .rendered_html :link {\n",
|
||||
" text-decoration: underline; \n",
|
||||
" }\n",
|
||||
" .rendered_html :hover {\n",
|
||||
" text-decoration: none; \n",
|
||||
" }\n",
|
||||
" .rendered_html :visited {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .rendered_html :focus {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .rendered_html :active {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .warning{\n",
|
||||
" color: rgb( 240, 20, 20 )\n",
|
||||
" } \n",
|
||||
" hr {\n",
|
||||
" color: #f3f3f3;\n",
|
||||
" background-color: #f3f3f3;\n",
|
||||
" height: 1px;\n",
|
||||
" }\n",
|
||||
" blockquote{\n",
|
||||
" display:block;\n",
|
||||
" background: #fcfcfc;\n",
|
||||
" border-left: 5px solid #c76c0c;\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" width:680px;\n",
|
||||
" padding: 10px 10px 10px 10px;\n",
|
||||
" text-align:justify;\n",
|
||||
" text-justify:inter-word;\n",
|
||||
" }\n",
|
||||
" blockquote p {\n",
|
||||
" margin-bottom: 0;\n",
|
||||
" line-height: 125%;\n",
|
||||
" font-size: 100%;\n",
|
||||
" }\n",
|
||||
"</style>\n",
|
||||
"<script>\n",
|
||||
" MathJax.Hub.Config({\n",
|
||||
" TeX: {\n",
|
||||
" extensions: [\"AMSmath.js\"]\n",
|
||||
" },\n",
|
||||
" tex2jax: {\n",
|
||||
" inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
|
||||
" displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
|
||||
" },\n",
|
||||
" displayAlign: 'center', // Change this to 'center' to center equations.\n",
|
||||
" \"HTML-CSS\": {\n",
|
||||
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
|
||||
" }\n",
|
||||
" });\n",
|
||||
"</script>\n"
|
||||
],
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 1,
|
||||
"text": [
|
||||
"<IPython.core.display.HTML at 0x7f5cc95696a0>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"> I am not well versed with Ensemble filters. I have implemented one for this book, and made it work, but I have not used one in real life. Different sources use slightly different forms of these equations. If I implement the equations given in the sources the filter does not work. It is possible that I am doing something wrong. However, in various places on the web I have seen comments by people stating that they do the kinds of things I have done in my filter to make it work. In short, I do not understand this topic well, but choose to present my lack of knowledge rather than to gloss over the subject. I hope to master this topic in the future and to author a more definitive chapter. At the end of the chapter I document my current confusion and questions. In any case if I got confused by the sources perhaps you will to, so documenting my confusion can help you avoid the same.\n",
|
||||
"\n",
|
||||
"The ensemble Kalman filter (EnKF) is very similar to the unscented Kalman filter (UKF) of the last chapter. If you recall, the UKF uses a set of deterministically chosen weighted sigma points passed through nonlinear state and measurement functions. After the sigma points are passed through the function, we find the mean and covariance of the points and use this as the filter's new mean and covariance. It is only an approximation of the true value, and thus suboptimal, but in practice the filter is highly accurate. It has the advantage of often producing more accurate estimates than the EKF does, and also does not require you to analytically derive the linearization of the state and measurement equations. \n",
|
||||
"\n",
|
||||
"The ensemble Kalman filter works in a similar way, except it uses a *Monte Carlo* method to choose the sigma points. The filter starts by randomly generating a large number of points distributed about the filter's initial state. This distribution is proportional to the filter's covariance $\\mathbf{P}$. In other words 68% of the points will be within one standard deviation of the mean, 95% percent within two standard deviations, and so on. Let's look at this in two dimensions. We will use `numpy.random.multivariate_normal()` function to randomly create points from a multivariate normal distribution drawn from the mean (5, 3) with the covariance\n",
|
||||
"\n",
|
||||
"$$\\begin{bmatrix}\n",
|
||||
"32 & 15 \\\\ 15 & 40\n",
|
||||
"\\end{bmatrix}$$\n",
|
||||
"\n",
|
||||
"I've drawn the covariance ellipse representing two standard deviations to illustrate how the points are distributed."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import numpy as np\n",
|
||||
"from numpy.random import multivariate_normal\n",
|
||||
"from stats import (covariance_ellipse, \n",
|
||||
" plot_covariance_ellipse)\n",
|
||||
"\n",
|
||||
"mean = (5, 3)\n",
|
||||
"P = np.array([[32, 15],\n",
|
||||
" [15., 40.]])\n",
|
||||
"\n",
|
||||
"x,y = multivariate_normal(mean=mean, cov=P, size=2500).T\n",
|
||||
"plt.scatter(x, y, alpha=0.3, marker='.')\n",
|
||||
"plt.axis('equal')\n",
|
||||
"\n",
|
||||
"plot_covariance_ellipse(mean=mean, cov=p,\n",
|
||||
" variance=2.**2,\n",
|
||||
" facecolor='none')"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAsQAAAF2CAYAAABpv45gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0XOd55/nvXaoKxMZVXACQAHcRBAESoCjBksWWbXmR\nKUu0FbXtVjyynbTT6nMkO5P2eGnPSVoz3Tk9nagjJ+5x2lQcJ3E6Q9labMWKbdkStZDiAoAkuG8A\nAXDfSSxVdZf54+WtKoAASRCQILF+n3N0BBRu3bp1i0f64eHzPq8VhmGIiIiIiEiessf6AkRERERE\nxpICsYiIiIjkNQViEREREclrCsQiIiIiktcUiEVEREQkrykQi4iIiEheUyAWERERkbw2aoH49ddf\nx7Zt1qxZA0A6nebLX/4ypaWlVFZWsnbt2tF6KRERERGRUeOOxkk8z+Mb3/gGixYtwrIsAJ566il2\n7NhBZ2cnzc3NrFq1isbGRioqKkbjJUVERERERsWoVIi/+93vsmrVKqZOnZp5bO3atTz++OOUlpay\ncuVKGhsbee6550bj5URERERERs2IA/GxY8f44Q9/yB/+4R/2e3zv3r0sXLiQRx55hH/6p3+iurqa\nPXv2jPTlRERERERG1YgD8R/90R/x7W9/m0Qi0e/x7u5uiouLaW1t5ciRI5SUlHDp0qWRvpyIiIiI\nyKgaUQ/xm2++yaFDh3j44YcBCMOQMAwBKCoqoru7m5aWFgCeeOIJSkpKrjjHK6+8MpJLEBERERG5\nLh/+8IcHfXxEgXjz5s2sX78e284WmtetW0draysLFixg165d1NfXA7Bz504eeOCBQc+zaNGikVzG\nqJg8eTI//elPWbly5VhfSl7S/R9buv9jR/d+bOn+jy3d/7GTj/d+165dQ/5sRC0TTzzxBEEQZP5Z\nuXIlP/jBD3jqqad4+OGHefrppzl//jyvvvoqGzZsYPXq1SN5ORERERGRUTcqY9cG87WvfY3du3cz\nc+ZMJk6cyDPPPEN5efk79XIiIiIiIjdkVAPxb3/72+yJXZc1a9ZkNup4P3gvtG7kM93/saX7P3Z0\n78eW7v/Y0v0fO7r3Wdq6OYf+YIwt3f+xpfs/dnTvx5bu/9jS/R87uvdZCsQiIiIiktcUiEVEREQk\nrykQi4iIiEheUyAWERERkbymQCwiIiIieU2BWERERETymgKxiIiIiOQ1BWIRERERyWsKxCIiIiKS\n1xSIRURERCSvKRCLiIiISF5TIBYRERGRvKZALCIiIiJ5TYFYRERERPKaArGIiIiI5DUFYhERERHJ\nawrEIiIiIpLXFIhFREREJK8pEIuIiIhIXlMgFhEREZG8pkAsIiIiInlNgVhERERE8poCsYiIiIjk\nNQViEREREclrCsQiIiIiktcUiEVEREQkrykQi4iIiEheUyAWERERkbymQCwiIiIieU2BWERERETy\nmgKxiIiIiOQ1BWIRERERyWsKxCIiIiKS1xSIRURERCSvKRCLiIiISF5TIBYRERGRvKZALCIiIiJ5\nbcSB+JFHHmHGjBmMHz+euro6XnzxRQDS6TRf/vKXKS0tpbKykrVr1474YkVERERERps70hN8/etf\nZ82aNSQSCX71q1+xatUqzpw5w1/91V+xY8cOOjs7aW5uZtWqVTQ2NlJRUTEa1y0iIiIiMipGXCGu\nra0lkUgQhiGpVIri4mIsy2Lt2rU8/vjjlJaWsnLlShobG3nuuedG45pFREREREbNiCvEAI899hjP\nPPMM48aN46WXXqKwsJC9e/eycOFCHnnkEe6//36qq6vZs2fPaLyciIiIiMioGZVA/L3vfY+nn36a\n73//+zzyyCPs3LmT7u5uiouLaW1tpaGhgZKSEjo6OgZ9/uTJk0fjMkYkFosB741ryUe6/2NL93/s\n6N6PLd3/saX7P3Z07/sblUAM4Lou//7f/3v+8i//kldeeYWioiK6u7tpaWkB4IknnqCkpGTQ5z75\n5JOZr++++25Wrlw5WpclIiIiInnotddeY926dZnv77nnniGPHbVAHAnDkDAMWbBgAbt27aK+vh6A\nnTt38sADDwz6nMcee6zf96dPnx7ty7qm6DeksXht0f0fa7r/Y0f3fmzp/o8t3f+xkw/3vqamhpqa\nmsz3u3btGvLYES2qO378OGvWrOHChQt4nsf3v/99Tpw4wQc+8AEefvhhnn76ac6fP8+rr77Khg0b\nWL169UheTkRERERk1I2oQuw4Dj/+8Y/5xje+QSqVorq6mhdffJFJkybxta99jd27dzNz5kwmTpzI\nM888Q3l5+Whdt4iIiIjIqBhRIJ4yZQqvvPLK4Cd2XdasWcOaNWtG8hIiIiIiIu8obd0sIiIiInlN\ngVhERERE8poCsYiIjCrfN/+IiLxfKBCLiMio8X3YsiXBli0JhWIRed9QIBYRERGRvDbqG3OIiEj+\nchxoaEhmvhYReT9QIBYRkVGlICwi7zdqmRARERGRvKZALCIiIiJ5TYFYRERERPKaArGIiIiI5DUF\nYhERERHJawrEIiIiN0A78oncPBSIRUREhkk78oncXBSIRURERCSvaWMOERGRYdKOfCI3FwViERGR\nG6AgLHLzUMuEiIiIiOQ1BWIRERERyWsKxCIiIiKS1xSIRURERCSvKRCLiIiISF5TIBYRERGRvKZA\nLCIiIiJ5TYFYRETGlO+j7Y9FZEwpEIuIyKgaTsD1fdiyJcGWLQmFYhEZMwrEIiIyahRwReT9SFs3\ni4jImHEcaGhIZr4WERkLCsQiIjJqHAeWLRtewFUQFpGxppYJEREZNb4Pzc0JmpvVMiEi7x8KxCIi\neUZTHURE+lPLhIhIHokWvYHp3R3tdgX1BIvI+5ECsYiIjCoFYRF5v1EgFhHJI+92BTdqzVBIFpH3\nMgViEZE8826F03e6PUNEZLRoUZ2IiLwnaLGfiIwVVYhFROQds2xZEtu+dlVa1WQRGUuqEIuIyKjI\nrfBGAbe5OTG2FyUich1GFIg9z+MLX/gCZWVlTJgwgQ996EPs3LkTgHQ6zZe//GVKS0uprKxk7dq1\no3LBIiLy3hMF4C1bbmxDjmixn6rDIjIWRhSIfd9n/vz5bN68mXPnzvGpT32KBx98EICnnnqKHTt2\n0NnZyY9+9CO+9KUv0dnZOSoXLSIi7203EnAd58YW/Kn3WERGakSBOJFI8J3vfIeysjIAHn30Ufbv\n38+pU6dYu3Ytjz/+OKWlpaxcuZLGxkaee+65UbloERF5bxksAN9owB2OkVamRURglBfVrV+/nvLy\nciZPnszevXtZuHAhjzzyCPfffz/V1dXs2bNnNF9ORETG0MAZwwPDr2YQi8j7xagF4vPnz/PVr36V\nP//zP8eyLLq7uykuLqa1tZWGhgZKSkro6OgY9LmTJ08ercu4YbFYDHhvXEs+0v0fW7r/Y8PzQizL\nwXXt992997yQ9evN142N4LrWsH5+rXPD9T/n3nuj44uv+zUi+rM/tnT/x47ufX+jEoiTySSrV6/m\ns5/9LL/zO78DQFFREd3d3bS0tADwxBNPUFJSMujzn3zyyczXd999NytXrhyNyxIRec+KAqPjwJ13\nBmN9Oe8Z/YN0eF2heDhhW0Tyx2uvvca6desy399zzz1DHjviQOz7Pp/73OdYsGABf/Inf5J5fMGC\nBezatYv6+noAdu7cyQMPPDDoOR577LF+358+fXqklzVs0W9IY/Haovs/1nT/332+D5cuJSguLiKd\nhvPnz4z1JQ3bwoXm3+fPZx/LbZMY7OfXEt0XgDNnrlyQN9ptGPqzP7Z0/8dOPtz7mpoaampqMt/v\n2rVryGNHHIi/8pWvYNs23/ve9/o9/vDDD/P000+zatUqmpub2bBhAz/84Q9H+nIiIjeFaBHapElF\n11XhfDf6cdNp8+/Lf5N6TYOF1RvdXCP3/TU0JEf9/CIiVzOiQNze3s4zzzxDYWEh48ePzzz+8ssv\n87WvfY3du3czc+ZMJk6cyDPPPEN5efmIL1hEZKyNVjh1nOv76/5rBcGB13Mj15dOw7PPFgHw0EPd\n1x2KR4OCroiMtREF4srKSoJg6N63NWvWsGbNmpG8hIjIe8p7LbwNvB54d68vN3xfrbo7Gt7p84tI\n/hrVsWsiIjL63o0gGIuZynD09fUY7JeDG7m+od7fYJVuBWEReScoEIuIDMNYVSmHeq3Brud6r29g\n4BwYhG+0NeRGnqd+YREZSwrEIiLD9F4LZwOv53qu73r6kjduND9fsWLwQDpYGB/tIGtZEATvzo53\n8N77bEXk3TGirZtFROTmFATQ1ubS1uZylaUi79j2zI4Dy5Yl8Txobn5nt2X2vFDbP4vkOVWIRUTy\n0LVaP2wb5s3zMl+P1nmHogqtiIwlBWIRkTx1tfDpOFBfP7xgO5Ke44FtFr5vKsOuC0uXjm4P8cDr\ndF1L0ytE8pwCsYhInrpWgL3eXuTIO7EILgyHV6G+lsHG1IGCsEi+UyAWEclD17P47WqB2ffB82Db\ntgRhaPp9b9RgbRaaOSwi7yYFYhGRm9RI+nKjwOz7sHx5st9INt+HpqYEe/e6OA7Mneth21cPsDdS\njR5Ohfp63+ONBG31N4vc/BSIRUTeJ4YTzK5VAb6eYOj7ZtKE62b7iXPZNlRWepke36EmNLxTM4Vv\n9LzDXeynecgiNz8FYhGR94HhBrMgMDN8w3DoY661qG758iSua84RBGaRW/T69fVJli5NYttkwvBo\nBUdVZEXk3aZALCLyHnYjc3GjCQ1Ru8ONBstYrH9lOPdahjN/eDhtCsMJ1u9Gn7F6mUXygwKxiMh7\nkO9fWZUdbjBznJFPaMjdgc69yv8xrhUc38nA+k5TEBa5+SkQi4iMMc8L8f0rtz+OQmjU9vLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f5ca4cf0198>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The Algorithm\n",
|
||||
"\n",
|
||||
"As I already stated, when the filter is initialized a large number of sigma points are drawn from the initial state ($\\mathbf{x}) and covariance ($\\mathbf{P}$. From there the algorithm proceeds very similarly to the UKF. During the prediction step the sigma points are passed through the state transition function, and then perturbed by adding a bit of noise to account for the process noise. During the update step the sigma points are translated into measurement space by passing them through the measurement function, they are perturbed by a small amount to account for the measurement noise. The Kalman gain is computed from the \n",
|
||||
"\n",
|
||||
"We already mentioned the main difference between the UKF and EnKF - the UKF choses the sigma points deterministically. There is another difference, implied by the algorithm above. With the UKF we generate new sigma points during each predict step, and after passing the points through the nonlinear function we reconstitute them into a mean and covariance by using the *unscented transform*. The EnKF just keeps propagating the originally created sigma points; we only need to compute a mean and covariance as outputs for the filter! \n",
|
||||
"\n",
|
||||
"Let's look at the equations for the filter. As usual, I will leave out the typical subscripts and superscripts; I am expressing an algorithm, not mathematical functions. Here $N$ is the number of sigma points, $\\chi$ is the set of sigma points.\n",
|
||||
"\n",
|
||||
"**Initialize Step**\n",
|
||||
"\n",
|
||||
"$$\\chi \\sim \\mathcal{N}(\\mathbf{x}_0, \\mathbf{P}_0)\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"This just says to select the sigma points from the filter's initial mean and covariance. In code this might look like\n",
|
||||
"\n",
|
||||
" N = 1000\n",
|
||||
" sigmas = multivariate_normal(mean=x, cov=P, size=N)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"** Predict Step**\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\chi &= f(\\chi, u) + v_Q \\\\\n",
|
||||
"x &= \\frac{1}{N} \\sum_1^N \\chi\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$\n",
|
||||
"That is short and sweet, but perhaps not entirely clear. The first line passes all of the sigma points through a use supplied state transition function and then adds some noise distributed according to the $\\mathbf{Q}$ matrix. In Python we might write"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
" for i, s in enumerate(sigmas):\n",
|
||||
" sigmas[i] = fx(x=s, dt=0.1, u=0.)\n",
|
||||
"\n",
|
||||
" sigmas += multivariate_normal(x, Q, N)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The second line computes the mean from the sigmas. In Python we will take advantage of `numpy.mean` to do this very concisely and quickly.\n",
|
||||
"\n",
|
||||
" x = np.mean(sigmas, axis=0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We can now optionally compute the covariance of the mean. The algorithm does not need to compute this value, but it is often useful for analysis. The equation is\n",
|
||||
"\n",
|
||||
"$$\\mathbf{P} = \\frac{1}{N-1}\\sum_1^N[\\chi-\\mathbf{x}^-][\\chi-\\mathbf{x}^-]^\\mathsf{T}$$\n",
|
||||
"\n",
|
||||
"$\\chi-\\mathbf{x}^-$ is a one dimensional vector, so we will use `numpy.outer` to compute the $[\\chi-\\mathbf{x}^-][\\chi-\\mathbf{x}^-]^\\mathsf{T}$ term. In Python we might write\n",
|
||||
"\n",
|
||||
" P = 0\n",
|
||||
" for s in sigmas:\n",
|
||||
" P += outer(s-x, s-x)\n",
|
||||
" P = P / (N-1)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**Update Step**\n",
|
||||
"\n",
|
||||
"In the update step we pass the sigma points through the measurement function, compute the mean and covariance of the sigma points, compute the Kalman gain from the covariance, and then update the Kalman state by scaling the residual by the Kalman gain. The equations are\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\chi_h &= h(\\chi, u)\\\\\n",
|
||||
"z_{mean} &= \\frac{1}{N}\\sum_1^N \\chi_h \\\\ \\\\\n",
|
||||
"P_{zz} &= \\frac{1}{N-1}\\sum_1^N [\\chi_h - z_{mean}][\\chi_h - z_{mean}]^\\mathsf{T} + R \\\\\n",
|
||||
"P_{xz} &= \\frac{1}{N-1}\\sum_1^N [\\chi - x^-][\\chi_h - z_{mean}]^\\mathsf{T} \\\\\n",
|
||||
"\\\\\n",
|
||||
"K &= P_{xz} P_{zz}^{-1}\\\\ \n",
|
||||
"x & = x + K[z - v_R + \\chi_h]\n",
|
||||
"\\end{aligned}\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This is very similar to the linear KF and the UKF. Let's just go line by line.\n",
|
||||
"\n",
|
||||
"The first line,\n",
|
||||
"\n",
|
||||
"$$\\chi_h = h(\\chi, u),$$\n",
|
||||
"\n",
|
||||
"just passes the sigma points through the measurement function $h$. We name the resulting points $\\chi_h$ to distinguish them from the sigma points. In Python we could write this as\n",
|
||||
"\n",
|
||||
" sigmas_h = h(sigmas, u)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The next line computes the mean of the measurement sigmas.\n",
|
||||
"$$z_{mean} = \\frac{1}{N}\\sum_1^N \\chi_h$$\n",
|
||||
"\n",
|
||||
"In Python we can write\n",
|
||||
"\n",
|
||||
" z_mean = np.mean(sigmas_h, axis=0)\n",
|
||||
" \n",
|
||||
"Now that we have the mean of the measurement sigmas we can compute the covariance for every measurement sigma point, and the *cross variance* for the measurement sigma points vs the sigma points. That is expressed by these two equations\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"P_{zz} &= \\frac{1}{N-1}\\sum_1^N [\\chi_h - z_{mean}][\\chi_h - z_{mean}]^\\mathsf{T} + R \\\\\n",
|
||||
"P_{xz} &= \\frac{1}{N-1}\\sum_1^N [\\chi - x^-][\\chi_h - z_{mean}]^\\mathsf{T}\n",
|
||||
"\\end{aligned}$$\n",
|
||||
"\n",
|
||||
"We can express this in Python with\n",
|
||||
"\n",
|
||||
" P_zz = 0\n",
|
||||
" for sigma in sigmas_h:\n",
|
||||
" s = sigma - z_mean\n",
|
||||
" P_zz += outer(s, s)\n",
|
||||
" P_zz = P_zz / (N-1) + R\n",
|
||||
"\n",
|
||||
" P_xz = 0\n",
|
||||
" for i in range(N):\n",
|
||||
" P_xz += outer(self.sigmas[i] - self.x, sigmas_h[i] - z_mean)\n",
|
||||
" P_xz /= N-1"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Computation of the Kalman gain is straightforward $K = P_{xz} P_{zz}^{-1}$.\n",
|
||||
"\n",
|
||||
"In Python this is the trivial\n",
|
||||
"\n",
|
||||
" K = np.dot(P_xz, inv(P_zz))\n",
|
||||
"\n",
|
||||
"Finally, we compute the new state \n",
|
||||
"x & = x + K[z - v_R + \\chi_h]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Outstanding Questions\n",
|
||||
"\n",
|
||||
"All of this should be considered as *my* questions, not lingering questions in the literature. However, I am copying equations directly from well known sources in the field, and they do not address the discrepencies.\n",
|
||||
"\n",
|
||||
"First, in Brown [1] we have all sums multipied by $\\frac{1}{N}$, as in \n",
|
||||
"\n",
|
||||
"$$ \\hat{x} = \\frac{1}{N}\\sum_{i=1}^N\\chi_k^{(i)}$$\n",
|
||||
"\n",
|
||||
"The same equation in Crassidis [2] reads (I'll use the same notation as in Brown, although Crassidis' is different)\n",
|
||||
"\n",
|
||||
"$$ \\hat{x} = \\frac{1}{N}\\sum_{i=1}^N\\chi_k^{(i)}$$\n",
|
||||
"\n",
|
||||
"The same is true in both sources for the sums in the computation for the covariances. Crassidis, in the context of talking about the filter's covariance, states that $N-1$ is used to ensure an unbiased estimate. Given the following standard equations for the mean and standard deviation (p.2 of Crassidis), this makes sense for the covariance.\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\begin{aligned}\n",
|
||||
"\\mu &= \\frac{1}{N}\\sum_{i=1}^N[\\tilde{z}(t_i) - \\hat{z}(t_i)] \\\\\n",
|
||||
" \\sigma^2 &= \\frac{1}{N-1}\\sum_{i=1}^N\\{[\\tilde{z}(t_i) - \\hat{z}(t_i)] - \\mu\\}^2\n",
|
||||
"\\end{aligned}\n",
|
||||
" $$\n",
|
||||
" \n",
|
||||
"However, I see no justification or reason to use $N-1$ to compute the mean. If I use $N-1$ in the filter for the mean the filter does not converge and the state essentially follows the measurements without any filtering. However, I do see a reason to use it for the covariance as in Crassidis, in contrast to Brown. Again, I support my decision emperically - $N-1$ works in the implementation of the filter, $N$ does not.\n",
|
||||
"\n",
|
||||
"My second question relates to the use of the $\\mathbf{R}$ matrix. In Brown $\\mathbf{R}$ is added to $\\mathbf{P}_{zz}$ whereas it isn't in Crassidis and other sources. I have read on the web notes by other implementers tthat adding R helps the filter, and it certainly seems reasonable and necessary to me, so this is what I do. \n",
|
||||
"\n",
|
||||
"My third question relates to the computation of the covariance $\\mathbf{P}$. Again, we have different equations in Crassidis and Brown. I have chosen the implementation given in Brown as it seems to give me the behavior that I expect (convergence of $\\mathbf{P}$ over time) and it closely compares to the form in the linear KF. In contrast I find the equations in Crassidis do not seem to converge much.\n",
|
||||
"\n",
|
||||
"I am not comfortable saying either book is wrong; it is quite likely that I missed some point that makes each set of equations work. I can say that when I implemented them as written I did not get a filter that worked. Between reading implementation notes on the web and reasoning about various issues I have chosen the implementation in this chapter, which does in fact seem to work correctly. I have yet to explore the significant amount of original literature that will likely definitively explain the discrepencies. I would like to leave this here in some form even if I do find an explanation that reconciles the various differences, as if I got confused by these books than probably others will as well.\n",
|
||||
"\n",
|
||||
"## References\n",
|
||||
"- [1] Crassidis, John L., and John L. Junkins. *Optimal estimation of dynamic systems*. CRC press, 2011.\n",
|
||||
"\n",
|
||||
"- [2] Brown, Robert Grover, and Patrick Y.C. Hwang. *Introduction to Random Signals and Applied Kalman Filtering, With MATLAB\u00ae excercises and solutions.* Wiley, 2012."
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
File renamed without changes.
File renamed without changes.
@@ -136,36 +136,41 @@ Kalman filter as covered only work for linear problems. Extended Kalman filters
|
||||
|
||||
* [**Chapter 10: Unscented Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/10_Unscented_Kalman_Filters/Unscented_Kalman_Filter.ipynb)
|
||||
|
||||
|
||||
Unscented Kalman filters (UKF) are a recent development in Kalman filter theory. They allow you to filter nonlinear problems without requiring a closed form solution like the Extended Kalman filter requires.
|
||||
|
||||
|
||||
* [**Chapter 11: Designing Nonlinear Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/11_Designing_Nonlinear_Kalman_Filters/Designing_Nonlinear_Kalman_Filters.ipynb)
|
||||
[**Chapter 11: Ensemble Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/11_Ensemble_Kalman_Filter/Ensemble_Kalman_Filter_Kalman_Filters.ipynb)
|
||||
|
||||
Discusses the ensemble Kalman Filter, which uses a Monte Carlo approach to deal with very large Kalman filter states in nonlinear systems.
|
||||
|
||||
|
||||
* [**Chapter 12: Designing Nonlinear Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/12_Designing_Nonlinear_Kalman_Filters/Designing_Nonlinear_Kalman_Filters.ipynb)
|
||||
|
||||
EKF and UKF are linear approximations of nonlinear problems. Unless programmed carefully, they are not numerically stable. We discuss some common approaches to this problem.
|
||||
|
||||
|
||||
* [**Chapter 12: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/12_HInfinity_Filters/HInfinity_Filters.ipynb)
|
||||
* [**Chapter 13: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/13_HInfinity_Filters/HInfinity_Filters.ipynb)
|
||||
|
||||
H-inifinity filters are a form of filter that is very robust in the presence of non-Gaussian noise. They do not perform as well as Kalman filters, but are less likely to diverge.
|
||||
|
||||
|
||||
* [**Chapter 13: Numerical Stability**](not implemented)
|
||||
|
||||
Not written yet.
|
||||
|
||||
|
||||
* [**Chapter 14: Smoothing**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/14_Smoothing/Smoothing.ipynb)
|
||||
|
||||
Kalman filters are recursive, and thus very suitable for real time filtering. However, they work well for post-processing data. We discuss some common approaches.
|
||||
|
||||
|
||||
|
||||
* [**Chapter XX: Numerical Stability**](not implemented)
|
||||
|
||||
Not written yet.
|
||||
|
||||
|
||||
* [**Chapter 15: Particle Filters**](not implemented)
|
||||
* [**Chapter XX: Particle Filters**](not implemented)
|
||||
|
||||
Not written yet
|
||||
|
||||
|
||||
* [**Chapter 16: Multihypothesis Tracking**](not implemented)
|
||||
* [**Chapter XX: Multihypothesis Tracking**](not implemented)
|
||||
|
||||
Not written yet.
|
||||
|
||||
|
||||
@@ -0,0 +1,415 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:11e1418ff66f54f922969f038a790dfe89ea35dd53a49764ef0cb06a9c32f67b"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"%matplotlib inline\n",
|
||||
"from __future__ import division, print_function\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import numpy as np\n",
|
||||
"from numpy.random import randn\n",
|
||||
"from numpy import array, asarray\n",
|
||||
"\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0, '../code') # allow us to import book_format\n",
|
||||
"import book_format\n",
|
||||
"book_format.load_style()\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"html": [
|
||||
"<style>\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
|
||||
"\n",
|
||||
" div.cell{\n",
|
||||
" width: 850px;\n",
|
||||
" margin-left: 0% !important;\n",
|
||||
" margin-right: auto;\n",
|
||||
" }\n",
|
||||
" div.text_cell code {\n",
|
||||
" background: transparent;\n",
|
||||
" color: #000000;\n",
|
||||
" font-weight: 600;\n",
|
||||
" font-size: 11pt;\n",
|
||||
" font-style: bold;\n",
|
||||
" font-family: 'Source Code Pro', Consolas, monocco, monospace;\n",
|
||||
" }\n",
|
||||
" h1 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
"\t}\n",
|
||||
"\t\n",
|
||||
" div.input_area {\n",
|
||||
" background: #F6F6F9;\n",
|
||||
" border: 1px solid #586e75;\n",
|
||||
" }\n",
|
||||
"\n",
|
||||
" .text_cell_render h1 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 30pt;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#c76c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 1em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: wrap;\n",
|
||||
" } \n",
|
||||
" h2 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h2 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 20pt;\n",
|
||||
" font-style: italic;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#c76c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" } \n",
|
||||
" h3 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h3 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-size: 18pt;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#d77c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 2em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" h4 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h4 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" color:#d77c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 0.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" h5 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h5 {\n",
|
||||
" font-weight: 300;\n",
|
||||
" font-style: normal;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" margin-bottom: 0em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" div.text_cell_render{\n",
|
||||
" font-family: 'Arimo',verdana,arial,sans-serif;\n",
|
||||
" line-height: 135%;\n",
|
||||
" font-size: 125%;\n",
|
||||
" width:750px;\n",
|
||||
" margin-left:auto;\n",
|
||||
" margin-right:auto;\n",
|
||||
" text-align:justify;\n",
|
||||
" text-justify:inter-word;\n",
|
||||
" }\n",
|
||||
" div.output_subarea.output_text.output_pyout {\n",
|
||||
" overflow-x: auto;\n",
|
||||
" overflow-y: scroll;\n",
|
||||
" max-height: 50000px;\n",
|
||||
" }\n",
|
||||
" div.output_subarea.output_stream.output_stdout.output_text {\n",
|
||||
" overflow-x: auto;\n",
|
||||
" overflow-y: scroll;\n",
|
||||
" max-height: 50000px;\n",
|
||||
" }\n",
|
||||
" code{\n",
|
||||
" font-size: 70%;\n",
|
||||
" }\n",
|
||||
" .rendered_html code{\n",
|
||||
" background-color: transparent;\n",
|
||||
" }\n",
|
||||
" ul{\n",
|
||||
" margin: 2em;\n",
|
||||
" }\n",
|
||||
" ul li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.5em; \n",
|
||||
" }\n",
|
||||
" ul li li{\n",
|
||||
" padding-left: 0.2em; \n",
|
||||
" margin-bottom: 0.2em; \n",
|
||||
" margin-top: 0.2em; \n",
|
||||
" }\n",
|
||||
" ol{\n",
|
||||
" margin: 2em;\n",
|
||||
" }\n",
|
||||
" ol li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.5em; \n",
|
||||
" }\n",
|
||||
" ul li{\n",
|
||||
" padding-left: 0.5em; \n",
|
||||
" margin-bottom: 0.5em; \n",
|
||||
" margin-top: 0.2em; \n",
|
||||
" }\n",
|
||||
" a:link{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" a:visited{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" }\n",
|
||||
" a:hover{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" }\n",
|
||||
" a:focus{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" a:active{\n",
|
||||
" font-weight: bold;\n",
|
||||
" color:#447adb;\n",
|
||||
" }\n",
|
||||
" .rendered_html :link {\n",
|
||||
" text-decoration: underline; \n",
|
||||
" }\n",
|
||||
" .rendered_html :hover {\n",
|
||||
" text-decoration: none; \n",
|
||||
" }\n",
|
||||
" .rendered_html :visited {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .rendered_html :focus {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .rendered_html :active {\n",
|
||||
" text-decoration: none;\n",
|
||||
" }\n",
|
||||
" .warning{\n",
|
||||
" color: rgb( 240, 20, 20 )\n",
|
||||
" } \n",
|
||||
" hr {\n",
|
||||
" color: #f3f3f3;\n",
|
||||
" background-color: #f3f3f3;\n",
|
||||
" height: 1px;\n",
|
||||
" }\n",
|
||||
" blockquote{\n",
|
||||
" display:block;\n",
|
||||
" background: #fcfcfc;\n",
|
||||
" border-left: 5px solid #c76c0c;\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" width:680px;\n",
|
||||
" padding: 10px 10px 10px 10px;\n",
|
||||
" text-align:justify;\n",
|
||||
" text-justify:inter-word;\n",
|
||||
" }\n",
|
||||
" blockquote p {\n",
|
||||
" margin-bottom: 0;\n",
|
||||
" line-height: 125%;\n",
|
||||
" font-size: 100%;\n",
|
||||
" }\n",
|
||||
"</style>\n",
|
||||
"<script>\n",
|
||||
" MathJax.Hub.Config({\n",
|
||||
" TeX: {\n",
|
||||
" extensions: [\"AMSmath.js\"]\n",
|
||||
" },\n",
|
||||
" tex2jax: {\n",
|
||||
" inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
|
||||
" displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
|
||||
" },\n",
|
||||
" displayAlign: 'center', // Change this to 'center' to center equations.\n",
|
||||
" \"HTML-CSS\": {\n",
|
||||
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
|
||||
" }\n",
|
||||
" });\n",
|
||||
"</script>\n"
|
||||
],
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 1,
|
||||
"text": [
|
||||
"<IPython.core.display.HTML at 0x7f7aee4a4978>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"pos = 1\n",
|
||||
"vel = 1\n",
|
||||
"p_pos = 1\n",
|
||||
"p_move = .5\n",
|
||||
"\n",
|
||||
"data = [(0,pos+randn()*p_pos) for i in range(1000)]\n",
|
||||
"data = np.asarray(data)\n",
|
||||
"\n",
|
||||
"pos += vel\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"plt.scatter(data[:,0], data[:,1], marker=',')\n",
|
||||
"plt.show()\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAF2CAYAAABDD3/PAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwlNX9x/HPsiy5E8gSIpsoSisWmqJVa6Vq4k96sSN2\nUodK7Yi3FpnBtF6mU8TeQEZHtFUrWAtIYiNlRp3SVJlqFYRAqtZ2tF5AinhBQiDgQoCQZDebPL8/\naNY8ogs+nN3nkLxfM8zwTbLPOc931vWbw/c5J+A4jiMAAAAAn8kgvycAAAAAHI8opAEAAAAPKKQB\nAAAADyikAQAAAA8opAEAAAAPKKQBAAAADyikAQAAAA+MF9IXXnihcnJyVFBQoIKCAl199dWmhwAA\nAAB8N9j0BQOBgB588EFdd911pi8NAAAAWCMtrR0clggAAID+Li2F9OzZs1VcXKxvfvOb2rRpUzqG\nAAAAAHxlvJD+zW9+o23btmnr1q0666yz9J3vfEeJRML0MAAAAICvAk4a+zAcx1FhYaFeeOEFlZeX\nS5JWr16druEAAAAAl0mTJqXt2sYfNvy4QCBwWM/0uHHj0j3sgBAOh7VixQpVVlb6PZV+gXyaQy7N\nIp9mkU+zyKc55NKscDisxsbGtI5htLVj3759evrppxWLxRSLxTR37lyVlJRo/PjxJocBAAAAfGd0\nRbqrq0s///nP9fbbbysUCumcc87RU089pWAwaHIYAAAAwHdGC+kRI0bolVdeMXlJHAFtMmaRT3PI\npVnk0yzyaRb5NIdcHl84Ivw4x39wZpFPc8ilWeTTLPJpFvk0h1weXyikAQAAAA8opAEAAAAPKKQB\nAAAADyikAQAAAA8opAEAAAAPKKQBAAAADyikAQAAAA8opAEAAAAPKKQBAAAADyikAQAAAA8opAEA\nAAAPKKQBAAAADyikAQAAAA8opAEAAAAPKKQBAAAADyikAQAAAA8opAEAAAAPBvs9AQDAp9u+PUuv\nvpqQJCUSuYpEulVaGvN5VgAAiUIaAKzW3BxUVVVeMq6vb1VpqY8TAgAk0doBAAAAeMCKNABYLBQK\naNasDlcMALADhTQAWKytLaD583OS8Zlndvk4GwBAX7R2AIDFsrJ6UsYAAP9QSAMAAAAe0NoBABaL\nRLq1cuVBSVIi0aVIpNvnGQEAelFIA4DFSktjmjAhX5IUje7zeTYAgL5o7QAAAAA8oJAGAAAAPKCQ\nBgAAADygkAYAAAA8oJAGAAAAPEhLIb1+/XoNGjRIS5cuTcflAQAAAN8ZL6QTiYRuvfVWjRs3ToFA\nwPTlAQAAACsYL6QXLFigyZMna+TIkaYvDQAAAFjDaCG9c+dOPfLII7rllltMXhYAAACwjtGTDX/6\n05/q5z//ubKyslL+XDgcNjnsgBUKhSSRT1PIpznk0izyaRb5NIt8mkMuzerNZzoZK6T/8Y9/6L33\n3tPll1+e/JrjOJ/4s/PmzUv+vaKiQpWVlaamAQAAgAGqoaFB69atkyQFg0FVVFSkdTxjhfS///1v\nvfjiixo06KNukYaGBm3YsEH33nuv62dnzpzpiqPRqKlpDCi9v7GSPzPIpznk0izyaRb5NIt8mkMu\nj115ebnKy8slHcpnY2NjWscz1iN94403qqenJ/mnsrJSDz/88GFFNAAAANAfcCALAAAA4IHRhw37\nWrNmTbouDQAAAPiOFWkAAADAg7StSAMAjt3mzblqauqSJDlOgU48sVtjx7b7PCsAgEQhDQBW++CD\noK6+Oj8Z//GPBzR2rI8TAgAkUUgDgMUKChzNmtXhigEAdqCQBgCLtbYGNH9+TjKuqUn4OBsAQF88\nbAgAFhsyJHUMAPAPhTQAWCwUclLGAAD/0NoBABY75ZQurVx5UJKUSHQpEun2eUYAgF6sSAMAAAAe\nsCINABZrbg6qqiovGdfXt6q01McJAQCSWJEGAAAAPKCQBgAAADygtQMALFZYKD36aJskyXEcFRb6\nPCEAQBKFNABYbN8+adq0j44Ir69v9XE2AIC+aO0AAAAAPKCQBgAAADygtQMALBaJdHMgCwBYikIa\nACxWWhrThAmHeqSj0X0+zwYA0BetHQAAAIAHFNIAAACAB7R2AIDFtm/P0quvJiRJiUSuIpFulZbG\nfJ4VAECikAYAqzU3B1VVlZeM6+tbVVrq44QAAEm0dgAAAAAesCINABYLhQKaNavDFQMA7EAhDQAW\na2sLaP78nGR85pldPs4GANAXrR0AAACABxTSAGCxWCyQMgYA+IfWDgCwWHFxjxYubHPFAAA7UEgD\ngMV27x6k6ur8ZFxXd8DH2QAA+qKQBgCLjRjhuFakR4xwfJwNAKAvCmkAsFhLS8C1Il1by4o0ANiC\nhw0BAAAAD4wX0ldeeaVGjRqlwsJCnX766XryySdNDwEAA0YwmDoGAPjHeGvHz372My1dulRZWVl6\n7rnnNHnyZO3du1e5ubmmhwKAfm/oUMd1suHQofRIA4AtjBfSEyZMkCQ5jqN4PK78/HwFAux7CgBe\nxGJynWz42GOcbAgAtkhLj/TMmTOVk5OjK6+8Uk899ZRycnKO/CIAwGGysnpSxgAA/wQcx0nLvxMm\nEgktWrRIv/3tb7Vx40ZlZ2dLklavXq3zzz8/HUMOOKFQSJLU1cUKlQnk0xxyac66dTE1NX205lFW\n1qOKiiwfZ3T84/1pFvk0h1yaFQqFtGbNGk2aNCltY6Rt+7vBgwfrhhtu0MKFC7V69Wpdcsklye/N\nmzcv+feKigpVVlamaxoAcFzbuTOg667jQBYAOBoNDQ1at26dJCkYDKqioiKt46V9H+lPWvCeOXOm\nK45Go+meRr8UDoclkT9TyKc55NKcRKLgYzF5PVa8P80in+aQy2NXXl6u8vJySYfy2djYmNbxjPZI\nt7S0aOnSpdq/f3+ytWPXrl2aOHGiyWEAAAAA3xldkQ4Gg1q+fLluvfVWxeNxjR8/Xk8++aSKiopM\nDgMAA0ZJifuI8JIStr8DAFsYLaRHjBih1atXm7wkAAxo0aj7iHB6pAHAHhwRDgAW6+5OHQMA/EMh\nDQAAAHiQ9l07AADelZQ4WrKEHmkAsBGFNABY7MMPA5o+nR5pALARrR0AYDF6pAHAXhTSAAAAgAe0\ndgCAxeiRBgB7UUgDgMVaWtw90rW19EgDgC1o7QAAiwWDqWMAgH8opAEAAAAPaO0AAIvl5zuaNavD\nFQMA7EAhDQAWa20NaP78nGRcU5PwcTYAgL5o7QAAiw0enDoGAPiHj2QAsFg47GjhwjZXDACwA4U0\nAFisq0vatu2jrTrKynp8nA0AoC8KaQCw2ODB0okndrtiAIAd+EgGAIvt2hVQdfVHB7LU1HAgCwDY\ngocNAQAAAA8opAHAYoFA6hgA4B9aOwDAYsXF7l07iovZtQMAbEEhDQAWi8fdu3ZEIuzaAQC2oJAG\nAIvt2+c+2bC2lpMNAcAW9EgDgMWCwdQxAMA/FNIAYLHu7tQxAMA/tHYAgMV42BAA7EUhDQAWGzRI\nyspyxwAAO1BIA4DFYjFpy5aPGqNHjGDXDgCwBWsbAGAxx0kdAwD8w4o0AFistdW9/V1NDdvfAYAt\nWJEGAItxRDgA2ItCGgAsNnhw6hgA4B8+kgHAYnl5jmbN6nDFAAA7GF2RTiQSuuqqqxSJRDRs2DBd\ndNFF2rhxo8khAGBAaWs71CPd+6etjd4OALCF0RXp7u5unXrqqbrrrrsUiUR0//33q6qqSps3bzY5\nDAAMGPn57hXp/HxWpAHAFkYL6aysLP3yl79Mxtdcc41uueUWRaNRhcNhk0MBwIDQuyLdq66OXTsA\nwBZpfdjwxRdfVGlpKUU0AHj08ZMMOdkQAOyRtocN9+3bp5tuukn33nvvYd+jsDYjFApJIp+mkE9z\nyKU5+fmxw1o7yOux4f1pFvk0h1ya1ZvPdEpLIR2LxfTd735X3//+9/W9733vsO/Pmzcv+feKigpV\nVlamYxoAcNzbv5/WDgA4Wg0NDVq3bp0kKRgMqqKiIq3jGS+ku7u7dcUVV2js2LGaO3fuJ/7MzJkz\nXXE0GjU9jQGh9zdW8mcG+TSHXJoTCBR8LCavx4r3p1nk0xxyeezKy8tVXl4u6VA+Gxsb0zqe8W67\nGTNmaNCgQfr9739v+tIAMOB0daWOAQD+MboivXXrVtXU1Cg3N1eFhYXJrz/zzDM677zzTA4FAAPC\nsGHu7e+GDWP7OwCwhdFCevTo0erp6TF5SQAY0Pbtc/dI19bSIw0AtuCIcACwWHGxo4UL21wxAMAO\nFNIAYLFgUMrKcscAADtQSAOAxXbsCGj69PxkXFNzQBMm+DghAEASZ2QBAAAAHlBIA4DFBg9OHQMA\n/MNHMgBYbOhQ9/Z3Q4fysCEA2IJCGgAstmePe/u7mhq2vwMAW1BIA4DFOJAFAOxFIQ0AFmtvd69I\nL1vGijQA2IKHDQHAYvF46hgA4B8KaQAAAMADWjsAwGIlJe4jwktK6JEGAFtQSAOAxaLRgKqrPzrZ\nsK7ugI+zAQD0RWsHAFhs0KDUMQDAP3wkA4DFeNgQAOxFawcAWIx9pAHAXhTSAGCx1lZONgQAW1FI\nA4DFTjjB0ZIlba4YAGAHCmkAsFhPjxSLuWMAgB0opAHAYrt2ube/q6lh+zsAsAWFNABYbNQod2vH\nqFG0dgCALSikAcBiO3YENH26e0V6wgQfJwQASGIfaQCwWCCQOgYA+IdCGgAAAPCA1g4AsFhxsaOF\nC9tcMQDADhTSAGCx3bvdu3bU1rJrBwDYgtYOALBYKJQ6BgD4hxVpALBYbq6jWbM6XDEAwA4U0gBg\nsdbWgObPz0nGNTUJH2cDAOiLQhoALDZypPthw5EjWZEGAFtQSAOAxXjYEADsxcOGAGAxx0Line truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f7aee499b70>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"from filterpy.kalman import KalmanFilter\n",
|
||||
"from numpy.random import multivariate_normal\n",
|
||||
"import matplotlib.cm as cm\n",
|
||||
"\n",
|
||||
"vel = 20.\n",
|
||||
"kf = KalmanFilter(dim_x=2, dim_z=1)\n",
|
||||
"kf.F = array([[1, 1], [0, 1]], dtype=float)\n",
|
||||
"kf.H = array([[1, 0]])\n",
|
||||
"kf.x = array([[0., vel]]).T\n",
|
||||
"kf.P *= 100\n",
|
||||
"kf.R *= 10\n",
|
||||
"kf.Q *= .0001\n",
|
||||
"\n",
|
||||
"steps = 6\n",
|
||||
"colors = cm.winter(np.linspace(0, .5, steps))\n",
|
||||
"pcolors = cm.Reds(np.linspace(.5, 1, steps))\n",
|
||||
"for i in range(steps):\n",
|
||||
" kf.predict()\n",
|
||||
" x,y = multivariate_normal(mean=kf.x.T[0], cov=kf.P, size=1000).T\n",
|
||||
" plt.scatter(x,y, marker=',', c=pcolors[i])\n",
|
||||
" \n",
|
||||
" z = (i+1)*vel\n",
|
||||
" kf.update(z)\n",
|
||||
" x,y = multivariate_normal(mean=kf.x.T[0], cov=kf.P, size=1000).T\n",
|
||||
" plt.scatter(x,y, marker='.', c=colors[i], alpha=0.7)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAF2CAYAAAB+nR6pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xt8VPWd//HXOWfOXDOTewgkIVwUFOll1V6ou6DYmy1I\ntYs/iwgFlVr8WUpru+62bO2y27rWrVu3+mirgFdqpRWv7W5bb6gtdq3YXxUQBbkkQCD3yWQu5/b7\n42QmM8kkDCRcEj7Px8PHI2cy58yZbwDf883n+/kqjuM4CCGEEEIIIQalnuwbEEIIIYQQYiSQ4CyE\nEEIIIUQBJDgLIYQQQghRAAnOQgghhBBCFECCsxBCCCGEEAWQ4CyEEEIIIUQBJDgLIYQQQghRgCEH\n53g8zpe+9CXKy8spLS3lhhtuAMAwDK655hoikQj19fVs2LBhyDcrhBBCCCHEyeIZ6gVWrlzJrl27\n2Lp1K1VVVWzbtg2AO+64g7feeouGhga2bNnCnDlzmDFjBrW1tUO+aSGEEEIIIU40ZSg7B8bjccrL\ny3n99dc566yzcr73oQ99iJUrV7JgwQIAPvnJTzJ37lxuvPHGod2xEEIIIYQQJ8GQSjV27NiBoihs\n3LiR6upqzjnnHB5//PHM96ZOncrChQv5xS9+wbRp03j77beH5aaFEEIIIYQ40YZUqtHZ2UkqleK9\n995j7969/OEPf2Du3Lm88847xGIxioqKePPNNznvvPMIh8Ps27dvuO5bCCGEEEKIE2pIwTkYDGJZ\nFl//+tfxer1ceOGFTJkyhc2bNxMKhYjFYrzxxhsArFixgnA43O8azz777FBuQQghhBBCiIJcfPHF\nQzp/SMF50qRJKIrS73HHcZgyZQrbtm3j3HPPBWDr1q3Mmzcv73XOPvvsodzGqFdeXs5jjz3GrFmz\nTvatnNJknAonY1UYGafCyVgVRsapMDJOhZOxKkx5eTkvv/zykK8zpBrn0tJSZs2axQ9/+ENM02TT\npk3s2LGDGTNmcMUVV3DnnXfS0dHBCy+8wObNm7nsssuGfMNCCCGEEEKcDENuR7d27Vq++MUvUlJS\nQm1tLQ899BDV1dWsXLmS7du3U1dXR2lpKWvXrqWmpmY47lkIIYQQQogTbsjBecKECbzwwgv9L+zx\nsGbNGtasWTPUlxBIOUuhZJwKJ2NVGBmnwslYFUbGqTAyToWTsTpxZMvtEUL+UhRGxqlwMlaFkXEq\nnIxVYWScCiPjVDgZqxNHgrMQQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUQAJzkIIIYQQQhRAgrMQ\nQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUQAJzkIIIYQQQhRAgrMQQgghhBAFkOAshBBCCCFEASQ4\nCyGEEEIIUQAJzkIIIYQQQhRAgrMQQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUQAJzkIIIYQQQhRA\ngrMQQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUQAJzkIIIYQQQhRAgrMQQgghhBAFkOAshBBCCCFE\nASQ4CyGEEEIIUQAJzkIIIYQQQhRAgrMQQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUQAJzkIIIYQQ\nQhRAgrMQQgghhBAFkOAshBBCCCFEASQ4CyGEEEIIUYAhB+cLL7yQQCBAOBwmHA6zePFiAAzD4Jpr\nriESiVBfX8+GDRuGfLNCCCGEEEKcLJ6hXkBRFO666y6WLl2a8/gdd9zBW2+9RUNDA1u2bGHOnDnM\nmDGD2traob6kEEIIIYQQJ9ywlGo4jtPvsQ0bNvCVr3yFSCTCrFmzmDFjBhs3bhyOlxNCCCGEEOKE\nG5bg/I//+I9UVlbyyU9+ku3btwOwY8cOpk6dysKFC/nFL37BtGnTePvtt4fj5YQQYlh4o20EGndm\n/vNG2072LQkhhDiFDblU4/bbb2f69OlYlsXq1au59NJL2bp1K7FYjKKiIt58803OO+88wuEw+/bt\ny3uN8vLyod7GqKbrOiDjdCQyToWTsXIZh/biX3dr5jixbBXhCWdkjmWcCidjVRgZp8LIOBVOxqow\n6XEaqiEH5/POOy/z9fe+9z3uuusutm3bRigUIhaL8cYbbwCwYsUKwuFw3musXr068/XMmTOZNWvW\nUG9LCCGEEEKcxl588UU2bdoEgKZpzJw5c8jXHHJw7ktRFBzHYcqUKWzbto1zzz0XgK1btzJv3ry8\n5yxfvjznuKWlZbhva0RLf4qUcRmcjFPhZKxcAcPAn3VsGAadWWMi41Q4GavCyDgVRsapcDJWA5s+\nfTrTp08H3HF6+eWXh3zNIQXnjo4O/vCHPzB79mwAvv/97zNmzBimTZvGFVdcwZ133smcOXPYsmUL\nmzdv5r777hvyDQshxHCxImVEl9yccyyEEEIMZEjB2TAMvvWtb/HOO++g6zof/vCHeeqpp/B4PKxc\nuZLt27dTV1dHaWkpa9eupaamZrjuWwghhiwVLoVw6cm+DSGEECPEkIJzRUUFr7/+ev4LezysWbOG\nNWvWDOUlhBBCCCGEOCXIlttCCCGEEEIUQIKzEEIIIYQQBZDgLIQQQgghRAEkOAshhBBCCFGAYe/j\nLIQQQpxq9NZm1NbDmWO7rBKjrOIk3pEQYiSS4CyEEGLUU1sPo6/+ZubYWHUbSHAWQhwlKdUQQggh\nhBCiABKchRBCCCGEKICUagghhBj17LJKtzwj61gIIY6WBGchhBCjnlFWITXNQoghk1INIYQQQggh\nCiAzzkKIUcUbbUPrbM0cW5EyUuHSk3hHQgghRgsJzkKIUUXrbCW87tbMcXTJzVBAcM4XuCkvPy73\nKIQQYmSS4CyEEAwQuEcgmXEXQojjR4KzEGJUUy2TQOPOzPGpECSPZ7g91hl3IYQQRybBWQgxqliR\nspzZYkXzEL73XzPHp0KQlHA7fGQrbSHEiSTBWQgxqqTCpTkhNHu2eTB9A7cVKRv2exPDT7bSFkKc\nSBKchRCC/oH7VDdQuYd8ABBCiONHgrMQYkQqtE74VAySw3FP2eUeyZlzMWsnZcbjVKjjFkKI0UiC\nsxBiRCq0TvhUnEke9nvy6ITW/yhzeDrVTMtW2kKIE0mCsxBCiBFLttIWQpxIEpyFEGIEyi73UJNx\nfCf5foQQ4nQgwVkIMSKdirXLaSdiE5Lscg9vtO2UHQshhBhNJDgLIUakU7F2Oe1E92k+lcdCCCFG\nE/Vk34AQQgghhBAjgcw4CyGOixNRriCGn/zchBBiYBKchRDHxem8rfSpXH99JKfzz00IIY5EgrMQ\nQgwzqTkeGr21GbX1cObYLqt0284JIcRJJsFZCHFSSWnAsRutY6e2HkZf/c3MsbHqNunVLIQ4JUhw\nFkIcE8NQANB1J+/3Cy1XkNKAY3c8xu5El5nkC/9CCHGqkuAshMg4UhjOft799xcBsHhxV97nD1au\nkB2W1GR8KLfc73owemZeT4YTXWaSL/zbJ+zVhRDi6EhwFkIAhYXh4ZIdlpKzLx/W64HMWo90dlml\nW56RdSyEEKeCYQvOL730ErNmzeKee+7hmmuuwTAMrr/+ejZs2EBpaSm333478+fPH66XE0KcRLru\nsHhxV+brITENYgtWYPsCwPH7Vf1onJUeyd07BmOUVUhNsxDilDQswdk0TW6++WbOPvtsFMX9Ve8d\nd9zBW2+9RUNDA1u2bGHOnDnMmDGD2tra4XhJIcQQJBL9fxl+tGF4uGakfZueIrXkZuI1k4/pfG+0\nraByj9E4Kz0auneM1vAvhBidhiU4/9d//Rdz5szht7/9beaxDRs2sHLlSiKRCLNmzWLGjBls3LiR\nG2+8cTheUghxjAxD4e673Q+4CxYoOQH4eJZnZBvOsKR1tuJp2JUp+TCr605a+DoVZ7VPlXsa6D5G\nQ/gXQpw+hhycDx48yH333cef/vSnnOC8Y8cOpk6dysKFC5k7dy7Tpk3j7bffHurLCSFGgeEOS75N\nT/Vee8nNxy0YHimEnoqz2qfKPaXvIzHjEmxUqB6Pz3dQejQLIUaUIQfnm266iW9961v4fL6cx2Ox\nGEVFRbz55pucd955hMNh9u3bl/ca5eXlQ72NUU3XdUDG6UhknAr3la+oAHg8I//X4sahvTnHuq4T\nzPNnwDC6SSxb1fu80oq8z+t7Lej9M2Uc2os/K4Qmlq0iPOEMjIP7oK0ZJ5V0Z75NA9+mpwa8lxOp\n0PE5Wqm9u7EPN+EkEzhGinikBBWHcHkV3vETBrwPGxUeWIMGaIDz3TuInDl1yPczUsi/U4WRcSqc\njFVh0uM0VEMKzq+88grvvfceV1xxBQCO4+A47q96Q6EQsViMN954A4AVK1YQDofzXmf16tWZr2fO\nnMmsWbOGcltCiCMIBDQADGMUNP4qrcgJxJTmn73Uq+ugum7QS6UDcJpVNgZ/bf2R76GtGf/Pev8d\nG45OIac6+3ATyndWoqSP5y+CDQ9gf/cOyBOchRDiRHvxxRfZtGkTAJqmMXPmzCFfc0jB+bXXXuOP\nf/wjqqpmHtu0aRNvvvkmU6ZMYdu2bZx77rkAbN26lXnz5uW9zvLly3OOW1pahnJbo076U6SMy+Bk\nnAo3qsZKD0LV+NzHjvF9BQ4dyClrSCxbhTFmXGacAoaBP+v5hmHQ2dLS73Gzuo7UkpuxAmE6j/Je\n0uUgqmWCaeAEwziWecy1yd5AGCO7nvwY7ikfn2mQb/7GNI28f64y99HegVbA80erUfV37ziScSqc\njNXApk+fzvTp0wF3nF5++eUhX3NIwXnFihWsWLEic3zRRRdx9dVXs3TpUn7wgx9w5513MmfOHLZs\n2cLmzZu57777hnq/QohhkkjYGIZywhYEjgaFLmq0fYFj7hLStyY5OftyfM89dsy1yalwKdLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f7ad12f6da0>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 99
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"def update(mu1, var1, mu2, var2):\n",
|
||||
" if var1 == 0.0:\n",
|
||||
" var1=1.e-80\n",
|
||||
" \n",
|
||||
" if var2 == 0:\n",
|
||||
" var2 = 1e-80\n",
|
||||
" \n",
|
||||
" mean = (var1*mu2 + var2*mu1) / (var1+var2)\n",
|
||||
" variance = 1. / (1/var1 + 1/var2)\n",
|
||||
" return (mean, variance)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"def predict(mu1, var1, mu2, var2):\n",
|
||||
" return (mu1+mu2, var1+var2)\n",
|
||||
"\n",
|
||||
"pos = 0\n",
|
||||
"pos_p = 100.\n",
|
||||
"move = 1\n",
|
||||
"move_p = .1\n",
|
||||
"sense_p = 8\n",
|
||||
"\n",
|
||||
"N = 10\n",
|
||||
"colors = cm.winter(np.linspace(0, .5, N))\n",
|
||||
"pcolors = cm.Reds(np.linspace(.5, 1, N))\n",
|
||||
"for i in range(N):\n",
|
||||
" pos, pos_p = predict(pos, pos_p, move, move_p)\n",
|
||||
" x,y = multivariate_normal(mean=(pos,i), cov=[[0.001, 0],[0, p_pos]], size=1000).T\n",
|
||||
" plt.scatter(x,y, marker=',', c=colors[i]) \n",
|
||||
" print(pos, pos_p)\n",
|
||||
" \n",
|
||||
" pos, pos_p = update(pos, pos_p, i, sense_p)\n",
|
||||
" x,y = multivariate_normal(mean=(i, pos), cov=asarray([[0.001, 0],[0, pos_p]]).T, size=1000).T\n",
|
||||
" plt.scatter(x,y, marker='.', c=pcolors[i], alpha=0.5) \n",
|
||||
" \n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"1 100.1\n",
|
||||
"1.0740055504162813 7.50795559666975\n",
|
||||
"2.0381768183201006 3.9730859395971154\n",
|
||||
"3.02550842348427 2.7546779733418045\n",
|
||||
"4.018974755764885 2.1491012228687625\n",
|
||||
"5.014956796940505"
|
||||
]
|
||||
},
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
" 1.7940228898505703\n",
|
||||
"6.012217081465885 1.5654022438192952\n",
|
||||
"7.010217725218009 1.409220211690132\n",
|
||||
"8.008687415099768 1.2981611057964606\n",
|
||||
"9.007474523188765 1.2169185743509554\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAF2CAYAAAB3dr7TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlwXed55/nv2e65Oy72HeAK7uImSiQlilqoxbJlS1ls\nJ3Ec2U56etJTdqZqYle5exJPp6q7UzUzSVzdVeOOl7TidDlxbMmRZFMSJVISJVIbKYo7uGLfcfft\nrPMHSFpyaFsUKQEHeD5VqMK9FID36L3A/Z33POd5Fd/3fYQQQgghhBDvoc72AIQQQgghhJiLJCgL\nIYQQQghxFRKUhRBCCCGEuAoJykIIIYQQQlyFBGUhhBBCCCGuQoKyEEIIIYQQVyFBWQghhBBCiKu4\n7qD8k5/8hG3bthEOh/nCF75w5flvfOMbGIZBIpEgkUiwZMmS6/1RQgghhBBCfGSuOyinUim++tWv\n8qUvfek9zyuKwu/8zu+Qz+fJ5/OcP3/+en+UEEIIIYQQHxn9er/Bzp07ATh06BClUunK877vI5v+\nCSGEEEKIoLphNcq/GIoVReHJJ5+koaGBjRs38tRTT92oHyWEEEIIIcSH7oYFZUVR3vP4M5/5DOfP\nn2dsbIw/+7M/47Of/Sy9vb036scJIYQQQgjxobru0ovLfnFFeeXKlVc+f+SRR7jzzjt55pln6Onp\nec9/9/zzz9+oIQghhBBCCPFL3XPPPdf039+woPyLK8rXYtWqVTdqGOIjUl9fz49//OMrNeoiWGT+\ngkvmLthk/oJN5i+46uvr2b9//zV/3XWXXnieR6VSwXEcXNelWq3iOA6PP/44mUwGz/N4+umnefHF\nF7n//vuv98cJIYQQQgjxkbjuFeXHHnuML37xi1cef//73+fP//zPOXHiBF/4whdwXZfly5fzj//4\nj/+q7EIIIYQQQoi56rqD8qOPPsqjjz56A4YigkZKZoJN5i+4ZO6CTeYv2GT+FhbZwlp8YPLHIthk\n/oJL5i7YZP6CTeZvYZGgLIQQQgghxFVIUBZCCCGEEOIqJCgLIYQQQghxFRKUhRBCCCGEuAoJykII\nIYQQQlyFBGUhhBBCCCGuQoKyEEIIIYQQVyFBWQghhBBCiKuQoCyEEEIIIcRVSFAWQgghhBDiKiQo\nCyGEEEIIcRUSlIUQQgghhLgKCcpCCCGEEEJchQRlIYQQQgghrkKCshBCCCGEEFchQVkIIYQQQoir\nkKAshBBCCCHEVUhQFkIIIYQQ4iokKAshhBBCCHEVEpSFEEIIIYS4CgnKQgghhBBCXIUEZSGEEEII\nIa5CgrIQQgghxALlOQ6e48z2MOYsCcpCCCGEEAuQ5ziceHI3J57cfSUsS3B+LwnKQgghhBDiqsF5\nodNnewBCCCGEEOKjp+o6qx964MrnEo7/NQnKQgghhBALlKrr7/n83cFZSFAWQgghhFjQhkoew8Wf\nrya3xXTaJSECEpSFEEIIIRa04aLDl5/tB6C/5PHEw4toj4ZmeVRzgwRlIYQQQoj3wbGs2R7Ch8Lw\nXW49f2DmwZJtszuYOUaCshBCCCHEr+E7Did/9BMAWu++A2We1fBubI7SENH52G0tpKs++0ctFid1\n2qMLu0Ha/JplIYQQQgjxvg2VPAqexqIHP0bF8Tgy7fA3bw0DSAkGEpSFEEIIIX4tRddZ9ZufAiCb\nz8/yaK6f5ziMlT1O5Xw+/9O+K89/7dZmttZrsziyuUWCshBCCCHE+3AubdE3nqfieFR8jfqIRmNU\noT2uzPbQrsnljUXGyi7963deeb4rqrIlpdB6/EkmSjbR2/5XQFaUhRBCCCHErzGWKfGjbz1ObU2E\ng42b2Lmilm2dJu3x4MapuOrTFVVpiyh8Pn2YyosKdbUJInqRmCEry8GdWSGEEEKIj5COy6L2FJGL\nJ7i5VCa55sHZHtIHcnljEWOyzPiTz/JnpkpdqoXcyT7yqJy5/RFenXL5WhW6HGdBbz6ysG9lFEII\nIYR4H3zHIf38HirpDPmWpeTt2R7R9blQ8Bh55VUSk0OYik/EruJaFm65QtJQWJVUmdzzLCee3L2g\nt7ZeuKcIQgghhBDvk++6JEKwqTNOeMddrPdVYuEQjdFg1SfDTKeLwYLDEyfGqUu0EevaxMa3nkbV\nNWoXdZF46zmSgyPU7ridQrkw28OdVRKUhRBCCCGuYqjgM5x3CeFS2vccdbUx7v+dj1GsVGZ7aNdl\nuOhwZNrmtUubi2woemyLRvB1g+rYCF61Skt3G7ll62iui0nphRBCCCGEeK/hvMvD/zDMv/3JOE+e\nKpLOltENY7aHdUO0xw36Sx4Ay2Iq0bpazGQcp1LFVjQmI3X89xNZCt7CvqFv4Z4iCCGEEEK8D/0F\nn9PLbqPxnnoc28Z3nMDvzBcPafz9A+0M795N+6hJJZcjfbGfRGsL4eU99L1+mPVZF25tZSG3iAv2\nLAshhBBCfIi64gptMYWtE4eY/OEU+8ol6pZ0s+QTDwQyLA+VPKouKEBE8VheH6EmHiU9PYZbreLa\nNomGBmobUtwS10mZOr05l57kwlxZDt4MCyGEEEJ8BDpiPv858TaFskv9iibUiSq6oYHvzfbQPrDh\nosNnnrxIV1Tlt9NH2GSNkM0VaFrVw5KuDnIj4ww8v4dYUxM1d93NUMVHrzoLNihLjbIQQgghxFU0\nmx4R1efoaJWzo3nGRyZwCgW8eJJ3JlzeGHF4Y8RhqODP9lA/kGQsjFUo4loW0bo6+l99nclTp4jU\n1mK2d/DtdyYxfIdEaGGGZJCgLIQQQogbzHcc/ID33vUdh4vP7MGIRFj90AN0pkxS0RB2vkhpZIxS\nxaXqKuy9UGE47872cK9Zf8lj2W23YoTDRBvqMcImruOgGiFali8jd+oUm175Mf7B/Rh+8I7vRrmu\noPyTn/yEbdu2EQ6H+cIXvnDledu2+dKXvkQymaS7u5sf/vCH1z1QIYQQQsx9vuMwuvtZRnc/G/yw\n7PmErBLL61UujubQelbjNLQyNTrF/71/mkPDFcJ6sPoot8V0nnh4EY892E0iYhBOJHCrVYaPHKN5\n9Qpa1q5i9Phxki2NxGsSOPjkLAnKH0gqleKrX/0qX/rSl97z/F/91V9x/PhxBgcHeeyxx/jiF7/I\n4ODgdQ1UCCGEEOKjFGmsxx4ZxDj8KkszZ1EPHyCUmSDZ3soXNyZY0Ri8bhDtUZUtjSFadJvivj2E\nEnG6btmMZobQzBAj7xwn0zeAGQ7T2tpAfNttVNAYKgW3Lvt6XFdQ3rlzJ4888gh1dXXvef6HP/wh\nX/7yl0kmk+zcuZNt27bx+OOPX9dAhRBCCDH3KbpOywP30fLAfYHsCvFuqq6jGSb5M2dRPBe7XCK1\nfBlaPMb0//wuSWxWNoaoCQerhtdzHKy338SdGKU4McHgobfJDY3gOR6xhjq0kEGmfwAlEqE3O7Oa\nPFwM9tWBD+qGvIJ9/71F7L29vaxYsYLPfe5zPPTQQ6xevZrTp0/fiB8lhBBCiDku6AH5suLwCKgK\n4dYWato6iMbDDL64H6uQZ+nqVQzmHP7kuUme+L02gtZIrFiq0rxuLeWxUUaPHEM1dPLjE4RTKWo6\nOyiMTVC79TZWmzGmyg4zDeUWnhsyq4ry3v95xWKReDzOsWPH2Lx5M4lEgoGBgV/69fX19TdiGOIj\nZFzamUjmLphk/oJL5i7Y5vv8nR4pMjZdxnMs8pZKLGbi+QqddQYrWmOzPbxr4lgWZiSCuWwJUxMZ\nMi/sIdbUSKy+jmhdLQXPoaNGpyuuEDWUwM1p1/Zt9H73O9R2thNraSZSU0NxYpzpiUlWPXAvoUiE\nc+M5apoM/uiZQZ76reWBO8Z3Mz7gjoofyopyLBajWCzy9ttvA/CVr3yFRCLxS7/+L/7iL658fscd\nd7Bz584bMSwhhBBCfERcy6I0nab3Bz8lkhnjbLSb+I67+fc/y/HUH7exonW2R3htzk8WmbBUouEI\nnpe5dGOigmbo2Pk83lg/rbvu5f9QD5F77h16Ew/R01Y728N+XyrFIkM/exqnWCA/PELnLTdjRiKc\nfzWPabuMnj5NU88KlB//Pe6KlTz3G5+gKWnO9rCv2YsvvshLL70EgKZp3HHHHdf8PT6UFeWenh5O\nnjzJpk2bADhx4gSf+tSnfunX//Ef//F7Hk9NTd2IYYkP0eWzSpmrYJL5Cy6Zu2Cbr/PnOw4jT+/G\nHh1npV2gbDiMJ0wu3/7l2HbgjnksazGac6gZPE3DimXUdt5P+vBhJo+fxEwlCbe08p1DOdyhCvEl\nEQany9Sbwbjh7eJ0mezIBK0b11McGaP/tTdIdbRT29VBJZNj4uw5SpNTxGqTlEdHKVSqrNZCgZvD\ntWvXsnbtWmDmd2///v3X/D2u62Y+z/OoVCo4joPrulSrVRzH4dOf/jTf/OY3yWaz7Nu3j4MHD/LI\nI49cz48Sc4TiuSjewm0TI4QQ4pfwPdxiAV1TibU00XXbLSxvC/OVu2pme2QfyEhZ4WTzGrxCnqnX\nX4NyET1sgu/iVW0MY2atcW/rdkbW3YVFcG7oO1dWyT/8+2RXbKJSKFCaTpMdGiY3OoYWMqhb0o1m\nhGi74w7Krd2U/YXb9eK6VpQvt3677Pvf/z7f+MY3+PrXv86pU6fo7OyktraW7373u7S3t1/3YMXs\nUjyX8Mm3APBqd6HqH6zeRwghxPwyXDLI3nwfodBBxo6coNw/zampl1n2iV1s7gzRVqMBwdq9Lqy4\ntI6dxS/kwNDwbAcjGqF22TIUw0AxTWpN6C/45Jzg7d82YqukfGhdsYypi/2Y8RhGKMTI0eOkOjuI\nNTfw0oGjPNOxlc/ZCsNFh/Zo8NrhXa/rCsqPPvoojz766FX/7Tvf+Q7f+c53rufbCyGEEGLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f7ac824d320>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 100
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
+11
-12
@@ -75,41 +75,40 @@
|
||||
"\n",
|
||||
"Unscented Kalman filters (UKF) are a recent development in Kalman filter theory. They allow you to filter nonlinear problems without requiring a closed form solution like the Extended Kalman filter requires.\n",
|
||||
"\n",
|
||||
"*Still early going.*\n",
|
||||
"[**Chapter 11: Ensemble Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/11_Ensemble_Kalman_Filter/Ensemble_Kalman_Filter_Kalman_Filters.ipynb)\n",
|
||||
"\n",
|
||||
"Discusses the ensemble Kalman Filter, which uses a Monte Carlo approach to deal with very large Kalman filter states in nonlinear systems.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"[**Chapter 11: Designing Nonlinear Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/11_Designing_Nonlinear_Kalman_Filters/Designing_Nonlinear_Kalman_Filters.ipynb)\n",
|
||||
"[**Chapter 12: Designing Nonlinear Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/12_Designing_Nonlinear_Kalman_Filters/Designing_Nonlinear_Kalman_Filters.ipynb)\n",
|
||||
"\n",
|
||||
"Works through some examples of the design of Kalman filters for nonlinear problems. *This is still very much a work in progress.*\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"[**Chapter 12: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/12_HInfinity_Filters/HInfinity_Filters.ipynb)\n",
|
||||
"[**Chapter 13: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/13_HInfinity_Filters/HInfinity_Filters.ipynb)\n",
|
||||
" \n",
|
||||
"Describes the $H_\\infty$ filter. \n",
|
||||
"\n",
|
||||
"*I have code that implements the filter, but no supporting text yet.*\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"[**Chapter 13: Numerical Stability**](not implemented)\n",
|
||||
"[**Chapter 14: Smoothing**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/14_Smoothing/Smoothing.ipynb)\n",
|
||||
" \n",
|
||||
"Kalman filters are recursive, and thus very suitable for real time filtering. However, they work well for post-processing data. We discuss some common approaches.\n",
|
||||
" \n",
|
||||
"[**Chapter XX: Numerical Stability**](not implemented)\n",
|
||||
"\n",
|
||||
"EKF and UKF are linear approximations of nonlinear problems. Unless programmed carefully, they are not numerically stable. We discuss some common approaches to this problem.\n",
|
||||
"\n",
|
||||
"*This chapter is not started. I'm likely to rearrange where this material goes - this is just a placeholder.*\n",
|
||||
"\n",
|
||||
"[**Chapter 14: Smoothing**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/14_Smoothing/Smoothing.ipynb)\n",
|
||||
" \n",
|
||||
"Kalman filters are recursive, and thus very suitable for real time filtering. However, they work well for post-processing data. We discuss some common approaches.\n",
|
||||
" \n",
|
||||
" \n",
|
||||
"[**Chapter 15: Particle Filters**](not implemented)\n",
|
||||
"[**Chapter XX: Particle Filters**](not implemented)\n",
|
||||
" \n",
|
||||
"Particle filters uses a Monte Carlo technique to filter. \n",
|
||||
"\n",
|
||||
"*This is not implemented, and I have not decided if I want to make it part of this book or not.*\n",
|
||||
" \n",
|
||||
"\n",
|
||||
"[**Chapter 16: Multihypothesis Tracking**](not implemented)\n",
|
||||
"[**Chapter XX: Multihypothesis Tracking**](not implemented)\n",
|
||||
" \n",
|
||||
"*Not implemented yet.*\n",
|
||||
"\n",
|
||||
|
||||
Reference in new issue
Block a user