Removed 1D example.
We already do this in MKF chapter.
This commit is contained in:
Roger Labbe committed
1 parent
9ea7fe6ed0
commit
daed96898f
1 file changed
-156
@@ -938,162 +938,6 @@
|
||||
"The takeaway point is that *g* and *h* are specified fully by the variance and covariances of the measurement and predictions at time *n*. In other words, we are picking a point between the measurement and prediction by a scale factor determined by the quality of each of those two inputs. That is all the Kalman filter is. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Compare Performance to Univariate Kalman Filter"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's implement a simple 1-D kalman filter using the Kalman filter from this chapter, and compare its output to the kalman filter from the previous chapter by plotting it. We will use a simple model of tracking an object that starts at x=0 and moves by 1 at each step. We will assume the arbitrary value 5 for the measurement noise and .02 for the process noise.\n",
|
||||
"\n",
|
||||
"First, let's implement the filter from the *Univariate Filter* chapter."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from __future__ import division\n",
|
||||
"import numpy as np\n",
|
||||
"from numpy.random import randn\n",
|
||||
"from filterpy.kalman import KalmanFilter\n",
|
||||
"\n",
|
||||
"# 1-D Kalman filter equations\n",
|
||||
"def predict(pos, variance, movement, movement_variance):\n",
|
||||
" return (pos + movement, variance + movement_variance)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"def update (mu1, var1, mu2, var2):\n",
|
||||
" mean = (var1*mu2 + var2*mu1) / (var1+var2)\n",
|
||||
" variance = 1 / (1/var1 + 1/var2)\n",
|
||||
" return (mean, variance)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now, let's implement the same thing using the Kalman filter. I will implement it as a function that returns a KalmanFilter object so that you can run the analysis code several times with the KalmanFilter initialized to the same starting conditions each time."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from filterpy.kalman import KalmanFilter\n",
|
||||
"\n",
|
||||
"def mkf_filter(R, Q):\n",
|
||||
" f = KalmanFilter(dim_x=1, dim_z=1, dim_u=1)\n",
|
||||
" f.P = 500.\n",
|
||||
" f.H = np.array([[1.]])\n",
|
||||
" f.F = np.array([[1.]])\n",
|
||||
" f.B = np.array([[1.]])\n",
|
||||
" f.Q = Q\n",
|
||||
" f.R = R\n",
|
||||
" \n",
|
||||
" return f"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Finally, let's compare the two. I will plot the data from the 1-D Kalman filter as a blue line, and the output of the filter from this chapter as red dots. I wrote it as a function so you can easily modify the parameters and regenerate the plots."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuYAAAEjCAYAAAB+VMTJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNX9//HXzCSTScJkJQvZSICwr4LsBlBBBbFFiyxa\nRa3aFqtoKy5VWapYlLq1aJVfK3yrIlqXUrWKC6uArArITkISEhKWhCQkmWSSub8/AiNDEkjCJAR4\nPx8PHsK95557Jo9reOdw7ueYDMMwEBERERGR88p8vgcgIiIiIiIK5iIiIiIizYKCuYiIiIhIM6Bg\nLiIiIiLSDCiYi4iIiIg0AwrmIiIiIiLNgIK5iIiIiEgzoGAuItLMuFwu/vGPfzBs2DDCw8OxWq1E\nRkbSrVs3Jk2axKJFi9xt58+fj9lsZsaMGV659/79+zGbzQwbNswr/YmISN35nO8BiIjIT1wuFzfc\ncAOfffYZISEhXH/99cTFxeF0Ovnxxx/54IMPWLduHePGjfO4zmQyeeX+J/vxVn8iIlJ3CuYiIs3I\nu+++y2effUbPnj1Zvnw5drvd43x5eTmrVq2qdp23NnHWZtAiIuePlrKIiDQj3377LQCTJk2qFsoB\nrFYrV155pbvNnXfeCcCMGTMwm83uXytWrACgsLCQ559/niuvvJK4uDj8/PyIjIzkhhtuYM2aNR59\nz58/nzZt2gCwbNkyj/5OXyqzadMmxo8fT0xMDH5+fsTExHDbbbexb98+735BREQuIZoxFxFpRsLD\nwwHYtWvXWduOGTOGgoIC/vOf/zB06FCGDh3qPpeYmAjA9u3beeKJJxgyZAijR48mNDSU9PR0Fi9e\nzP/+9z8WL17MddddB0CvXr144IEHePnll0lMTGTSpEnu/k7t++2332bSpEnYbDZuuOEG4uPj2bNn\nDwsXLuS///0vy5Yto0ePHuf8tRARudSYDP27pYhIs/H999/Tr18/KioqmDhxIj//+c+57LLLSExM\nrHHd9/z587nzzjuZPn06Tz31VLXzhYWFVFRUEBYW5nE8KyuLvn37EhQUxI4dO9zH09PTSUpKYujQ\noXzzzTfV+tu7dy9du3YlPj6eFStW0KpVK/e55cuXc/XVV9OjRw82bNhwLl8GEZFLkpayiIg0Iz17\n9uStt94iKiqKt99+m7Fjx9K2bVtCQ0O5/vrrWbRoUb3WgQcFBVUL5QCxsbHcdNNN7Nq1iwMHDriP\nn63v1157jfLycl588UWPUA64Z+U3bdrkEfZFRKRutJRFRKSZGTt2LGPGjGHp0qV8++23bN68mVWr\nVvHZZ5/x2WefMX/+fBYvXoyvr2+d+vv22295+eWXWbNmDYcPH6a8vNzjfFZWFnFxcXXuC6pmx2ua\nFc/NzQVgx44ddOrUqU59iohIFQVzEZFmyMfHh+HDhzN8+HCgqoziBx98wJ133skXX3zBa6+9xv33\n33/Wfj766CN+8YtfEBAQwPDhw2nbti2BgYGYzWaWLl3K8uXLKSsrq/O4jh49CsBf/vKXWtuYTCaK\ni4vr3KeIiFRRMBcRuQCYzWbGjh3L1q1befrpp/nmm2/qFMyffPJJbDYbGzZsoEOHDh7nsrKyWL58\neb3GERwcjMlk4ujRo4SEhNTrWhEROTOtMRcRuYC0aNEC+GktuMViAaCysrLG9nv37qVz587VQrnL\n5aqxHvrZ+hs4cCCGYbjLMYqIiPcomIuINCMLFy7kq6++qvElzJycHObNmwdASkoK8FN5xfT09Br7\nS0pKYvfu3Rw8eNB9zDAMpk+fzo4dO6pVegkNDQUgMzOzxv7uu+8+rFYrv//972ss6VhRUcHSpUvP\n9jFFRKQGKpcoItKMPPjgg7z88stER0czePBgdz3ytLQ0Pv30UxwOBwMHDuTrr7/Gz8+PgoIC4uLi\ncDqd/PKXvyQhIQGTycRtt91GQkICb7zxBr/+9a+JjIzkxhtvxNfXl2+//ZYdO3Zw9dVXu+uOnwz6\nAIMHD2b16tVcf/319OrVC19fX4YMGcIVV1wBVO1Oescdd1BRUcG1115LcnIylZWVZGZm8u233+J0\nOsnLyzsfXz4RkQub4UXZ2dnGbbfdZkRERBg2m83o3LmzsXz5co8206ZNM2JiYgx/f39j6NChxo8/\n/ujNIYiIXNAOHDhgvPbaa8aNN95odOzY0QgODjZ8fX2N6OhoY/jw4cYbb7xhVFRUeFzz5ZdfGoMH\nDzbsdrthMpkMs9ns8b13/vz5Rs+ePY3AwEAjIiLCuPHGG41t27YZ06dPr9bWMAwjNTXVGDNmjNGy\nZUvDYrEYZrPZmDFjhkeb7du3G7/61a+MpKQkw8/PzwgNDTU6d+5s3Hnnncann37aeF8gEZGLmNdm\nzI8dO8Zll11GSkoK9913HxEREaSmptKqVSs6duwIwOzZs3nmmWdYsGAB7du3Z+bMmaxatYpdu3a5\n102KiIiIiFyKvBbMH3/8cVauXMnKlStrPG8YBjExMdx///089thjADgcDiIjI5kzZw733HOPN4Yh\nIiIiInJB8trLnx9//DF9+/Zl3LhxREVF0atXL+bOnes+n5aWRm5uLiNGjHAfs9lspKSksHr1am8N\nQ0RERETkguS1YJ6amsqrr75Ku3btWLJkCQ888ACPPvqoO5zn5OQAEBUV5XFdZGSk+5yIiIiIyKXK\naxsMuVwu+vbtyzPPPANAjx492LNnD3PnzmXy5MlnvPb0cl0FBQXeGpaIiIiISJMLDg6u9zVemzGP\niYmhc+fOHsc6duxIRkYGANHR0QDk5uZ6tMnNzXWfExERERG5VHktmA8aNIidO3d6HNu9e7e7Bm9S\nUhLR0dEsWbLEfd7hcLBq1SoGDhzorWGIiIiIiFyQvLaU5cEHH2TgwIHMmjWLm2++mc2bN/PXv/6V\nZ599FqharjJlyhRmzZpFx44dSU5O5umnn8ZutzNx4sRa+23IPwPIpWHDhg0A9OnT5zyPRJozPSdS\nF3pO5Gz0jEhdnOtybK8F8z59+vDxxx/z+OOP86c//YnWrVvz9NNP85vf/MbdZurUqZSWljJ58mTy\n8/Pp378/S5YsITAw0FvDEBERERG5IHktmAOMHDmSkSNHnrHNtGnTmDZtmjdvKyIiIiJywfPaGnMR\nEREREWk4BXMRERERkWZAwVxEREREpBnw6hpzERERkebOMAycTicul6vO17Ru3RqoKvUslyaz2Yyv\nr2+1jTG9ScFcRERELhmGYeBwOLBarfUKWTabrZFHJs2ZYRi4XC4cDgc2m63RwrmWsoiIiMglw+l0\nYrVasVgsjTrzKRcXk8mExWLBarXidDob7T4K5iIiInLJcLlcmM2KP9IwZrO5Xkug6t1/o/UsIiIi\n0gxpplwaqrGfHQVzEREREZFmQMFcRERERKQZUDAXERERkRpNmjSJpKSk83b/oUOHMmzYsPN2/6am\nYC4iIiJygZs+fTpms5lDhw7VeH7o0KF06tSp3v2aTKbzuib/XO7/6quvsmDBAi+PqHGpjrmIiIjI\nJaAhAXfevHkYhtEIo6mbL7/8ssHXvvrqq0RERHD77bd7cUSNSzPmIiIiIlIjHx8ffH19m/y+paWl\n7vv7+Fw688gK5iIiIiKXkGXLlmE2m3n33XeZNWsWcXFx+Pv7c/XVV7Nv3z6PtqevMe/atSspKSk1\n9pucnOyxHvwvf/kLgwcPJiIiAn9/f7p3784//vGPatclJiZy3XXX8fXXX9OvXz/8/f15/vnngZrX\nmNel38TERLZv387y5csxm82YzWaPz1FWVsaMGTNITk7GZrMRFxfHQw895P6B4Hy5dH4EERERERG3\n5557Dh8fH6ZOncqxY8d47rnnuOWWW1i7dq1Hu1OXwEyYMIGnnnqKrKwsYmNj3cc3btzIvn37ePjh\nh93HXnrpJUaPHs348eMxmUx8/PHH3H333VRUVHDvvfd69L93717Gjh3LPffcw913301CQkKN969r\nvy+//DK/+93vsNvt/PGPfwSgRYsWABiGwZgxY1ixYgX33HMPnTt3Zvv27bz66qv8+OOPfPHFF+f6\npW0wBXMRERGRS1BZWRnr1q1zLxUJDQ3lgQce4Mcff6RLly7udqeuMR8/fjxPPvkk77//PlOmTHEf\nX7RoEb6+vvziF79wH9uzZw82m83958mTJ3PNNdcwZ84cj2BuGAb79u1j8eLFXH/99Wcdd136/dnP\nfsYf//hHIiMjmThxosf1Cxcu5IsvvmDZsmVcccUV7uN9+vTh1ltv5csvv2T48OFnHUdjUDAXERER\nqcVnaxfy+XeLGq3/a/uNY2T/CY3W/5ncdtttHuu3Bw8eDEBaWppHMD9V27Zt6d27N4sWLfII5u+9\n9x5XXXUVYWFh7mMnw7PT6aSoqAiXy8XQoUP58ssvKSoqwm63u9vGx8fXKZTXt9+avPfee7Rv357O\nnTtz5MgR9/GUlBRMJhNLly5VMBcRERGRxnP6kpBTl4tA1Yw5QH5+/hn7GT9+PA8//DDp6em0bt2a\nNWvWkJGRwcyZMz3a/ec//+FPf/oTP/zwA5WVlR7jKCgo8AjQbdq0qfPnqE+/Ndm9eze7du0iIiKi\n2jmTycThw4frPBZvUzAXERERucCdnEWu7eXFkpISj+UfABaLpca2ZyuPOG7cOKZOncqiRYvc/7XZ\nbIwZM8bdZtWqVYwZM4aUlBRef/11YmJisFqtfPrpp7z44ou4XC6PPv39/c/6GRvSb01cLhddunTh\n5ZdfrvF8TExMncbSGBTMRURERGoxsv+E87bUpD5at24NwM6dO92/P6myspK9e/cydOhQr9wrLi6O\ngQMHsmjRIh5++GHef/99rr32Wo+Z6n//+98EBASwZMkSrFar+/jXX399TveuT7+11Line truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x8042e48>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtIAAAEjCAYAAAAfTOaLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW9x/HvmS0LywDZAwkBDFJ2NKBQC5S6VGkR2ypQ\nvSK3lWpFLdpXpdoWpFwp3lZrLVqlimhl0Vq1t1C3QoAUFFAICFSWAEFICAESkpBtZs79I2TIEEKS\nCZPJTD7v1ysvZ86c55zfnDzil8fnPMcwTdMUAAAAgGaxBLsAAAAAIBQRpAEAAAA/EKQBAAAAPxCk\nAQAAAD8QpAEAAAA/EKQBAAAAPxCkASCIDh48KIvFomnTpgW7FABAMxGkAbRbt99+uywWi5599tlG\n97311ltlsVj00ksvBaQWwzACclwAQOAYPJAFQHu1fv16jRkzRoMGDVJ2dnaD+xUUFKhHjx6Kjo7W\n0aNHFR0dfclqcLlcysnJkdPpVEJCwiU7LgAg8BiRBtBufe1rX1O/fv20Y8cObdq0qcH9lixZIpfL\npSlTplzSEC1JNptNffv2JUQDQAgiSANo16ZPny5JWrRoUYP7/PnPf5ZhGJo+fbrWrFmj6dOnq3//\n/nI6nYqOjtbAgQM1Z84cVVRU1Gs7Z84cWSwWLVmyRKtWrdLo0aPVuXNndevWTVLDc6T37NmjWbNm\nKSMjQ3FxcYqMjFRaWpruvvtuHT58uN55MjMzvcc5ePCgJk+erNjYWEVFRWn48OFauXJlg9/vjTfe\n0HXXXaeYmBjveW677Tb9+9//rrfvW2+9pWuvvVbdunVTZGSk+vbtq5///OcqKSlp8PgAEK4I0gDa\ntalTp8rhcGjFihUqLS2t9/natWu1d+9eXXHFFRo2bJgWLFigjz76SFdccYXuuece3X333YqIiNDc\nuXN1ww03yO12X/A8b775pm6++WbFxMTovvvu08SJE30+P3+O9N/+9je98MIL6tmzp26//XY98MAD\n6t+/v15++WUNHz5cR44cueB5Dh06pKuuukq5ubmaOnWqJk2apM8//1w333yzMjMzffY1TVN33XWX\nJk+erOzsbN1yyy16+OGHNWbMGG3evFlvvfWWz/4//vGPdeutt2rv3r367ne/qwceeEAJCQlasGCB\nrrnmmgtePwAIayYAtHN33HGHaRiG+eKLL9b77Pvf/77PZzk5ORc8xi9/+UvTMAxz+fLlPttnz55t\nGoZhWq1W8/3336/X7sCBA6ZhGOa0adN8th85csSsqqqqt/8HH3xgWq1W85577vHZvmbNGtMwDNMw\nDHPu3Lk+n73//vumYRjmTTfd5LP9hRdeMA3DMIcPH24WFRX5fObxeMwjR45437/22mumYRjmd7/7\nXbOiosJn33nz5pmGYZgPP/xwvXoBIJxxsyGAdi8rK0ujR49WRkaGz1zpkydPKjk5WRERETp69Kg6\ndOjQ4DFOnjyp2NhY/fd//7f+/Oc/e7fPmTNHc+fO1S233FJvhFeqmdrRu3dv3XXXXXr55ZebVO/g\nwYNVVlam/fv3e7dlZmZq3LhxSktL0/79++uNcPfs2VPl5eUqKCjwbhs0aJB27dqlTZs26corr7zo\nOa+88krt2LFDBQUF6tKli89nbrdbiYmJstlsysvLa9J3AIBwYAt2AQAQbNdcc42+8pWvaMuWLdq+\nfbsGDx4sSXrttddUVVWlqVOnekN0WVmZnnnmGb399tvas2ePSktLVXc8oqEpFyNGjGh2XX/5y1/0\nyiuvKDs7W0VFRT7TRiIiIi7YZujQoRdcSi8lJUWffPKJ931ZWZl27typ2NjYRkP0mTNntHXrVsXE\nxOj3v//9BfdxOBzKy8vTqVOn1LVr16Z8PQAIeQRpAFDNTYczZ87UokWLvOtKL1q0yHuToSRVV1dr\n3Lhx2rx5swYNGqQpU6YoLi5Odrtdpmnq8ccfV2Vl5QWPn5iY2Kx6Zs6cqWeeeUbJycm68cYb1b17\nd0VFRUmSFi9erNzc3Au2O3+0uJbNZpPH4/G+LyoqkiR179690VpOnTolSTpx4oTmzp3b4H6GYai0\ntJQgDaDdIEgDgKQ777xTs2bN0tKlS/W///u/2rp1q3bt2qVhw4Z5R2zfffddbd68WdOmTav3YJa8\nvDw9/vjjDR6/OQ9cKSgo0B/+8AcNGjRIGzZsqDel5PXXX2/GN7uw2sDd0Ah6XU6nU1LNlJJt27a1\n+NwAEC5YtQMAJHXt2lW33nqrTp06pTfffNO7HF7taLQk7du3T5L0ne98p177tWvXXrJacnJyZJqm\nrr/++noh+ssvv1ROTk6Lz9GhQwcNHDhQhYWF2rJly0X37dixowYOHKjdu3frxIkTLT43AIQLgjQA\nnFUbmp955hm98cYb6tixo26//Xbv57169ZIkrVmzxqddTk6OHnnkkUtWR+151q9f7zMdo7S0VHff\nfXeDS+w11wMPPCBJuvfee1VcXOzzmWmaPjcOPvzww6qurta0adO8Uz3qKikpuehDbQAgHDG1AwDO\nqr3p8LPPPpMk/eAHP1DHjh29n3/729/WZZddpqeeeko7duzQ0KFDlZubq5UrV+pb3/qWli9ffknq\nSEhI0OTJk7V8+XINHTpU1113nYqLi/Xhhx8qOjpaQ4cO9WuKxfmLNP3whz9UVlaWXn31VV122WWa\nOHGi4uLilJeXp7Vr1+rmm2/W008/Lalmve3PPvtMzz77rPr06aMbbrhBPXv2VFFRkQ4ePKj169fr\nhhtu0N/+9rdLcg0AIBQ0aUQ6Ly9PU6dOVXx8vKKiojRgwACtW7cu0LUBQKv70Y9+JEk+NxnWio6O\n1urVq/X9739fO3fu1LPPPqvPP/9cv/rVr/Taa69d8HiGYTRrfnStl156SY8++qjKy8v13HPP6cMP\nP9SECRP073//W06ns9nHbKiOV155Ra+//roGDhyov/71r3r66aeVmZmpESNGaNKkST77PvPMM/rn\nP/+p0aNHKzMzU08//bTeeecdFRYWasaMGZo9e3azvycAhLJG15EuKirSFVdcodGjR2vGjBmKi4tT\nTk6OkpKS1K9fv9aqEwAAAGhTGg3Sjz76qNavX6/169e3Vk0AAABAm9fo1I533nnH+7/4EhISNGzY\nMC1cuLA1agMAAADarEaDdE5Ojp577jlddtll+uCDD/Tggw9q1qxZhGkAAAC0a41O7XA4HBoxYoSy\nsrK82x577DG9/fbb2rVrl3fb+UsnAQAAAKGk9gFUTdXoiHRycrL69+/vs61fv34NPp4WAAAAaA8a\nDdJf/epX9Z///Mdn2549e5SWlhaomgAAAIA2r9EHssycOVOjRo3SE088odtuu01bt27Vs88+q/nz\n5zfYprnD4mgfah9DnJGREeRK0JbRT9AU9BM0Bf0ETdGS6cmNjkhnZGTonXfe0RtvvKFBgwbpl7/8\npebNm6d7773X75MCAAAAoa5Jjwi/6aabdNNNNwW6FgAAACBkNOkR4QAAAAB8EaQBAAAAPxCkAQAA\nAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8QpAEAAAA/EKQBAAAA\nPxCkAQAAAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8QpAEAAAA/\nEKQBAAAAPxCkAQAAAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8QpAEAAAA/EKQBAAAAPxCkAQAAAD8Q\npAEAAAA/EKQBAAAAPzQapOfMmSOLxeLzk5yc3Bq1AQAAAG2WrSk79evXT5mZmd73Vqs1UPUAAAAA\nIaFJQdpqtSo+Pj7QtQAAAAAho0lzpHNyctS9e3f17t1bU6ZM0YEDBy66v2mal6Q4AAAAoK1qNEhf\nffXVWrJkid5//30tWrRI+fn5GjVqlE6ePNlgm4oVK+QqKbmkhQIAAABtiWE2c/j4zJkz6tWrl2bN\nmqWZM2d6txcXF3tfd+7SRcWLFytn8GB5PJ5LVy0AAABwCaWnp3tfO53OZrVt9vJ30dHRGjBggPbt\n29fgPoakzrNnq1tVVXMPDwAAAISEJt1sWFdFRYV2796tcePGXXQ/S26uEisrlXb11X4Xh/CyZcsW\nSVJGRkaQK0FbRj9BU9BP0BT0EzRF3VkVzdXoiPRPf/pTrVu3TgcOHNAnn3yi733veyovL9fUqVP9\nPikAAAAQ6hodkT5y5IimTJmiwsJCxcXFaeTIkfr444+VkpJy0Xae1FQZvXpdskIBAACAtqTRIL1s\n2bJmH9SUVLlggSIbCdsAAABAqGr2zYZNUbFsmezjx8swjEAcHgAAAAi6Zt9s2BSRkyYRogEAABDW\nAjIibYonGwIAACC8BSRIF5eeCMRhAQAAgDYjIEH6eFFeIA4LAAAAtBkEaQAAAMAPBGkAAADADwEK\n0kcDcVgAAACgzQhIkC4szg/EYQEAAIA2IzBBuihfHtMTiEMDAAAAbUJAgnS1u0pFJSyBBwAAgPAV\nkCAtMU8aAAAA4S1gQZp50gAAAAhnjEgDAAAAfghYkC5gLWkAAACEscBN7SBIAwAAIIwFdI40S+AB\nAAAgXAUsSLvc1SoqKQzU4QEAAICgCliQlqTjTO8AAABAmCJIAwAAAH4IcJBmCTwAAACEp8AGaR7K\nAgAAgDDFiDQAAADgh4AG6cLifHk87kCeAgAAAAiKgATpTlFOSZLb7dKpUpbAAwAAQPgJSJCO65Ls\nfV1YxDxpAAAAhJ8ABekk7+sC5kkDAAAgDAUkSMfWCdKsJQ0AAIBwFPAR6UKCNAAAAMJQs4L0/Pnz\nZbFYdP/99190v7pzpBmRBgAAQDhqcpD++OOPtWjRIg0ePFiGYVx0X58R6dMsgQcAALine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x808fe10>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"\n",
|
||||
"def plot_kf_compare(x0, p0, R, Q, move):\n",
|
||||
" # storage for filter output\n",
|
||||
" x_uni, x_kf = [], []\n",
|
||||
" p1, p2 = [], []\n",
|
||||
"\n",
|
||||
" # initialize the filters\n",
|
||||
" f = mkf_filter(R, Q)\n",
|
||||
" f.x[0] = 0.\n",
|
||||
" f.P[0, 0] = p0\n",
|
||||
" pos = (x0, p0)\n",
|
||||
" for i in range(50):\n",
|
||||
" z = i*move + randn()\n",
|
||||
" pos = update(pos[0], pos[1], z, R)\n",
|
||||
" f.update(z)\n",
|
||||
"\n",
|
||||
" x_uni.append(pos[0])\n",
|
||||
" x_kf.append(f.x[0])\n",
|
||||
"\n",
|
||||
" p1.append(pos[1])\n",
|
||||
" p2.append(f.P[0, 0])\n",
|
||||
"\n",
|
||||
" u = move + randn()\n",
|
||||
" pos = predict(pos[0], pos[1], u, Q)\n",
|
||||
" f.predict(u=u)\n",
|
||||
"\n",
|
||||
" plt.scatter(range(len(x_kf)), x_kf, c='r',\n",
|
||||
" s=80, label='Kalman')\n",
|
||||
" plt.title('State')\n",
|
||||
" plt.plot(x_uni, label='Univariate')\n",
|
||||
" plt.xlim([0, 50])\n",
|
||||
" plt.legend()\n",
|
||||
"\n",
|
||||
" plt.figure()\n",
|
||||
" plt.plot(p1)\n",
|
||||
" plt.scatter(range(len(x_kf)), p2, c='r', s=80)\n",
|
||||
" plt.title('Variance')\n",
|
||||
" plt.xlim([0, 50])\n",
|
||||
" plt.show()\n",
|
||||
" \n",
|
||||
"plot_kf_compare(x0=0., p0=500., R=5., Q=.2, move=1.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"As you can see, both filters produce the same results. Feel free to vary the initial guess, the measurement noise, and the process noise; so long as you make the same changes to both filters the output should be the same. This is a solid demonstration, albeit not a rigorous proof, that both filters in fact implement the same math for the 1-D case. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
|
||||
Reference in new issue
Block a user