Added animation to PF chapter.
This commit is contained in:
Roger Labbe committed
1 parent
41d8d246e0
commit
5794485341
3 files changed
+538
-445
No files matched your search
+91
-103
@@ -16,7 +16,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -262,7 +262,7 @@
|
||||
"<IPython.core.display.HTML object>"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -281,29 +281,29 @@
|
||||
"source": [
|
||||
"## Motivation\n",
|
||||
"\n",
|
||||
"Here is our problem. We have objects moving in a space, and we want to track them. Maybe the objects are fighter jets and missiles in the sky, or maybe we are tracking people playing cricket in a field. It doesn't really matter. Which of the filters that we have learned can handle this problem? Well, none of them are ideal. Let's think about the characteristics of this problem. \n",
|
||||
"Here is our problem. We have moving objects that we want to track. Maybe the objects are fighter jets and missiles, or maybe we are tracking people playing cricket in a field. It doesn't really matter. Which of the filters that we have learned can handle this problem? Unfortunately, none of them are ideal. Let's think about the characteristics of this problem. \n",
|
||||
"\n",
|
||||
"1. **multi-modal**: We want to track zero, one, or more than one object simultaneously.\n",
|
||||
"* **multimodal**: We want to track zero, one, or more than one object simultaneously.\n",
|
||||
"\n",
|
||||
"2. **occlusions**: One object can hide another, causing there to be only one measurement for multiple objects.\n",
|
||||
"* **occlusions**: One object can hide another, resulting in one measurement for multiple objects.\n",
|
||||
"\n",
|
||||
"3. **nonlinear behavior**: Aircraft are buffeted by winds, balls move in parabolas, and people collide into each other.\n",
|
||||
"* **nonlinear behavior**: Aircraft are buffeted by winds, balls move in parabolas, and people collide into each other.\n",
|
||||
"\n",
|
||||
"4. **nonlinear measurements**: Radar gives us the distance to an object. Converting that to an (x,y,z) coordinate requires a square root, which is nonlinear.\n",
|
||||
"* **nonlinear measurements**: Radar gives us the distance to an object. Converting that to an (x,y,z) coordinate requires a square root, which is nonlinear.\n",
|
||||
"\n",
|
||||
"5. **non-Gaussian noise:** as objects move across a background the computer vision can mistake part of the background for the object. \n",
|
||||
"* **non-Gaussian noise:** as objects move across a background the computer vision can mistake part of the background for the object. \n",
|
||||
"\n",
|
||||
"6. **continuous:** the object's position and velocity (i.e. the state space) can smoothly vary over time.\n",
|
||||
"* **continuous:** the object's position and velocity (i.e. the state space) can smoothly vary over time.\n",
|
||||
"\n",
|
||||
"7. **multivariate**: we want to track several attributes, such as position, velocity, turn rates, etc.\n",
|
||||
"* **multivariate**: we want to track several attributes, such as position, velocity, turn rates, etc.\n",
|
||||
"\n",
|
||||
"None of the filters we have learned work well with all of these constraints. \n",
|
||||
"\n",
|
||||
"* **Discrete Bayes filter**: This has most of the attributes. It is multimodal, can handle nonlinear measurements, and can be made to work with nonlinear behavior (in this book we only handled the linear case). However, it is discrete, not continuous, and it is *univariate*, not multivariate.\n",
|
||||
"* **Discrete Bayes filter**: This has most of the attributes. It is multimodal, can handle nonlinear measurements, and can be extended to work with nonlinear behavior. However, it is discrete and univariate.\n",
|
||||
"\n",
|
||||
"* **Kalman filter**: The Kalman filter produces optimal estimates for unimodal linear systems with Gaussian noise. None of these are true for our problem.\n",
|
||||
"\n",
|
||||
"* **Unscented Kalman filter**: The UKF handles nonlinear, continuous, multivariate problems. However, it is neither multimodal nor does it handle occlusions. It can handle noise that is modestly non-Gaussian, but does not do well with distributions that are very non-Gaussian or problems that are very nonlinear.\n",
|
||||
"* **Unscented Kalman filter**: The UKF handles nonlinear, continuous, multivariate problems. However, it is not multimodal nor does it handle occlusions. It can handle noise that is modestly non-Gaussian, but does not do well with distributions that are very non-Gaussian or problems that are very nonlinear.\n",
|
||||
"\n",
|
||||
"* **Extended Kalman filter**: The EKF has the same strengths and limitations as the UKF, except that is it even more sensitive to strong nonlinearities and non-Gaussian noise."
|
||||
]
|
||||
@@ -319,16 +319,16 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuoAAADaCAYAAADwgmciAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXts29eV5z+09eBDpGRJpCiJskkpsmVbkePYVizX4ya2\nY092imSyGWy27aDtIkD+ycw2EywCBFhMgkE3g2KLIChm2gVm0WkmmSKd3c40md2dRJ06beLYqfwc\nSlYkKxIpmZQoUU9SfEmyuH8I9+ZHitTLelr3AxiQ+fvxx/v7UTr33HPP+R5dMplMolAoFAqFQqFQ\nKDYVOzZ6AAqFQqFQKBQKhWI+ylFXKBQKhUKhUCg2IcpRVygUCoVCoVAoNiHKUVcoFAqFQqFQKDYh\nylFXKBQKhUKhUCg2IcpRVygUCoVCoVAoNiHKUVcoFAqFQqFQKDYhylFXKBSKbcxf/uVfcuzYMQoL\nC7HZbDz55JPcunVr3nmvvfYalZWVGI1GHnvsMdrb21OOJxIJ/vRP/xSr1UpBQQFPPfUUfr9/vW5D\noVAo7kuUo65QKBTbmN/+9rf8yZ/8CZcvX+bChQvk5ORw9uxZxsbG5Dnf//73eeONN/irv/orrly5\ngs1m4/HHH2dyclKe8+KLL/KP//iPvPvuu3zyySeEQiG+9rWvMTs7uxG3pVAoFPcFOtWZVKFQKBSC\nSCRCYWEh7733Hn/wB39AMpmkoqKC//yf/zOvvPIKAPF4HJvNxg9+8AOef/55JiYmsNls/PSnP+Xr\nX/86AD6fjz179vAv//IvnDt3biNvSaFQKLYsKqKuUCgUCkkoFGJ2dpZdu3YB4PF4GBwcTHG29Xo9\np06d4tKlSwBcu3aN6enplHMcDgf79++X5ygUCoVi+ShHXaFQKBSS7373uxw+fJimpiYAAoEAAGVl\nZSnn2Ww2eSwQCLBz505KSkpSzikrK2NwcHAdRq1QKBT3JzkbPYBMTExMbPQQFAqFYl0oLCzc6CFI\nXnrpJS5dusTFixfR6XSLnr+Uc9JR9l2hUGwXVsO+q4i6QqFQKPizP/szfv7zn3PhwgWcTqd83W63\nA8yLjA8ODspjdrudu3fvMjIyknJOIBCQ5ygUCoVi+ShHXaFQKLY53/3ud6WTvnfv3pRjLpcLu91O\nc3OzfC0ej3Px4kVOnDgBwJEjR8jNzU05x+fz0dHRIc9RKBQKxfLZlKkvWjbTtrBCoVCsBpsp/eOF\nF17gnXfe4Ze//CWFhYUy79xsNmMymdDpdLz44ou8/vrr1NXVUVtby/e+9z3MZjPf+MY3gDk7/dxz\nz/Hyyy9js9koLi7mpZde4tChQ5w9ezbrZyv7fv8SCoW4ePEiwWAQo9GIyWQCwGQycfjwYSwWy7zz\ngXmvrzebZRyKrctq2/dN76grFAqFYu348Y9/jE6n48yZMymvv/baa/z5n/85AC+//DKxWIwXXniB\nsbExjh8/TnNzs3S+AN58801ycnJ49tlnicVinD17lnfeeWdFeeyK+wOTyYTJZKKmpgaA7u7urOdu\ntGO8EQ76Qp+pFgwKwabUUdeuRlTERaFQ3G9sZxu3ne99u5HubG5W5zMUCnHjxg2AjNH+9f7MjRiP\nYvVYbRunctQVCoVCoVCsOhaLJcXJTHfYBaFQaN5r24FIJEIkEtnoYSg2OSr1RaFQKBQKxZqTHikW\nr4mUmI2KHlssFjmezRC93mzjUWwsylFXKBQKhUKx5oRCISKRCCaTSTroIqKsrXfYCDbCIV7onpWD\nrhAoR12hUCgUCsWa8vnnn9PW1iaLS4Ujqv3/VnROV5p3r6LmiqWiHHWFQqFQKBSrQibH1efz8f77\n7+Pz+Th16pR0ytMd1c1abJqNbEWf2vsQ+ffahchWu0/FxqIcdYVCoVAoFPfMQmoler0eh8NBfX29\nPHcrKp0s5mRr76Ompga3201XVxe1tbWcPHkSYEvcp2LzoBx1hUKhUCgUKyaTYovWoXU4HDzzzDPy\n/1vVUc20mBAa8VvpPhRbi1WTZ/zLv/xLjh07RmFhITabjSeffJJbt27NO++1116jsrISo9HIY489\nRnt7+2oNQaHY9GxXGTKFQnF/IpxXrQNbU1NDd3c3N27ckPbO4XDgcDjk+4aGhvD7/SnXqqmp2VLO\nuyiIdbvd+Hw+4Mvcc5F3f/LkSZ555hkaGhpSjm+l+1RsLKvmqP/2t7/lT/7kT7h8+TIXLlwgJyeH\ns2fPMjY2Js/5/ve/zxtvvMFf/dVfceXKFWw2G48//jiTk5OrNQyFYtOindCUs65QKO5HFisKtVgs\n2Gw2uru7aW5uxufzSdu4UOfSzUAmJzsSidDV1YXb7ZaBGOHAi8WLxWKhu7ubixcv4vP55j0j7fvU\n3KBIZ9VSXz744IOU/7/99tsUFhZy6dIl/uAP/oBkMsmbb77JK6+8wtNPPw3AW2+9hc1m42c/+xnP\nP//8ag1FoVAoFArFOpCpKFS8ls3pNJvNGAwGYrEY4XB43SLLq1HEmd7ASUTKF5OcFA59ps+/ceOG\nbH5kMpk4efKkirYrJGuWox4KhZidnWXXrl0AeDweBgcHOXfunDxHr9dz6tQpLl26pBx1xX2PkuNS\nKBT3I9nsmbaRkfZch8PBuXPnaGtrY2hoiMrKyjW3jWtVrOpwOOZdK5PkZCaHHpA57pFIhJ6eHvR6\nPQ0NDWqOUEjWzFH/7ne/y+HDh2lqagIgEAgAUFZWlnKezWajv79/rYahUGwqlPFVKBTbjUzdRysr\nKxkaGtrgka0OWruebcGRyaEX54ndBxGBV/OEQsuaOOovvfQSly5d4uLFi+h0ukXPX8o5CoVi9VA6\nvgqFYi3R7iAu5fhaK8Gs145munZ6+jFBpnQh7c8KhWDVHfU/+7M/4x/+4R/46KOPcDqd8nW73Q7A\n4OBgSuX34OCgPKZQ3E9sVmd4rfWKN+t9KxSKtSPT3/1ikeb0BkALXWs1WK3rLTS+dPua7bMzOfTK\nZioysWqqLzCX7vLzn/+cCxcusHfv3pRjLpcLu91Oc3OzfC0ej3Px4kVOnDixmsNQKDac7arwsl3v\nW6HYzmj/7oWKSzrZnHJxTKipAIvakI1UR1mOjVvouShbqVgqqxZRf+GFF3jnnXf45S9/SWFhocxJ\nN5vNmEwmdDodL774Iq+//jp1dXXU1tbyve99D7PZzDe+8Y3VGoZCoVgEVdSqUCjWgkgkgtvtxmQy\nzdutW2wnTxthFqopmdjsHUwzpfws9FwUisVYNUf9xz/+MTqdjjNnzqS8/tprr/Hnf/7nALz88svE\nYjFeeOEFxsbGOH78OM3NzSkSRortyf2WLrHZneH0CTT9tXu57ma+b4VCsfpoCyJXqoUuHHCtTOFm\nZCEbl8mWZnsuylYqlooumUwmN3oQ6UxMTMifCwsLN3AkivVgs0dI7mdCoRAXL14EUNq968h2tnHb\n+d63A4vlby90TDjqQEr0Of196xXYWc7n+Hy+BaPm91swSpGd1bZxaybPqFAoVpe1MPShUEg24VDa\nvQqF4l5ZrCvpQsfSU0aEk54eyFkPO7WcAFIoFMLtdtPV1UVFRcWiii8KxXJQjrpiw0k30JmM3GZD\nW62/2SaN5WCxWKitrZU/L2c8y32PQqHYHqzUPmxle2IymaioqCAWi+F2u1PmhkzPQ9lQxVJRjrpi\nU5AtcrIZEekiXV1d1NbWbumUEYvFwsmTJ+XPS2GrfE8KhWL9Wcg+rMQ53ahc7vTPXWjs4ly/38/7\n77+Px+PB5XJlndeUDVUsB+WoKxRbgLWcrO71eioypFAoMpEuR+h2u4Hl18NslG3RRsQXc6yFAx6N\nRhkYGKCtrY3Kysp1Ha/i/kQ56opNw1apghdRaJHTvV5j3SzPZL07CioUiq1DugKMKA6NRCL09PSg\n1+uz1sNs9UW/xWKhvr4eo9EoVWuyzWs1NTXrOn8oti7KUVdsKraK0cpkYLf6JLMcFmpeolAotjeZ\nbKDJZKK6uhqTyZRVyWWzLvqzOduZbH5TUxNNTU0pc8RCevIKxWIoR11x37MeDvRmnmTWkq2yC6JQ\nKNaXTCouWoS9FJHljWaxeSJTYCbd4d6Oc4Bi7VGOuuK+Zrs60OuJeqYKhSIT2WyDcIrTO3Yup5HQ\napJtnljO54qOqos1alosuLGddmYVS0M56grFKqAiy/NZaeMThUJxf6L9u8/UsXOhRkiwvsGWxT43\nfcdA3MdSdggWWsCowJIiHeWoK+5r7tWBXo5DKfIul6oDvxmc1bUaQzZJMoGajBSK7UU2JzS9qHIj\n7OJK54n0MYv8e4ViNVGOuuK+Z6WGc7nRjeVsn65G5GS527JLHcNadUDV5qQqFIrtTaYcb/FaJBKh\noaEBh8OxbruVC0XMlxIhX40xqp1ZRSaUo67Y9GyGyPNKWa2tzHtRSFjKuUuJdi8nlSVbs5BMxxQK\nxfbAZrNhNpvn2QQtkUiErq4u4Et1rY3WUc+G1u6t1hiVTVSkoxx1xaZmI3P2lutQrlUEZjWfQfrE\nkp4nmi3avVi3wUzHtOco51yh2L6Ibs6tra3Y7XbOnDmDw+GQ9kbruDc0NAAsWpS53mwlOUnF/YVy\n1BUKDZkiw9nyzjNFmJezfbrcnPZIJCLPXe6ioKamhnA4LB1yMbHcy3bvcllMAUJNdArF/YU2ai4a\nHnV3d2M0GmlqapKdSoEU23Ty5Elg89iEhZxyrV1WKNaCHRs9AMX9j3BGV4JwHNcjYiGM8Y0bN+R4\nM7220OvZ7mEhNYPFriEcbZibzJa73Soi5h6PR3YJTL++9hkv5MAvpH6wku9oOc9RsXZ8/PHHPPnk\nkzgcDnbs2MFbb72Vcvw73/kOO3bsSPl34sSJlHMSiQR/+qd/itVqpaCggKeeegq/3Line truncated
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuoAAADaCAYAAADwgmciAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWtsXNd1739DSRzOk6RIDodP8RFKJDWmrNCSRZtRbMmW\n6iR1kLqAb/KhzkXu9Yfr9sY1CgMGitoocl0E1zCMoEk+tEjj2nCdpnETF2lSNlL8YEnrLVMUH6LJ\nISkOOS++5s2HOPcD794+M5zhQ6Qokto/wLB45sw5+5yZWXudtdf6L108Ho+jUCgUCoVCoVAothQZ\nd3sACoVCoVAoFAqFYinKUVcoFAqFQqFQKLYgylFXKBQKhUKhUCi2IMpRVygUCoVCoVAotiDKUVco\nFAqFQqFQKLYgylFXKBQKhUKhUCi2IMpRVygUCoVCoVAotiDKUVcoFIp7mL/5m7/hyJEjZGdnY7PZ\nePLJJ7l+/fqS/V555RVKSkowGo08+uijdHV1Jbw+MzPDn/3Zn1FQUIDZbObrX/86Lpdrsy5DoVAo\ndiTKUVcoFIp7mA8//JA//dM/pb29nbNnz7J7924ee+wxJicn5T7f//73ef311/nbv/1bLly4gM1m\n4/HHHycUCsl9nn/+ed577z3effddPv74YwKBAF/72tdYWFi4G5elUCgUOwKd6kyqUCgUCkE4HCY7\nO5tf/epXfPWrXyUej1NcXMz//t//m5deegmAWCyGzWbjtdde49lnn2V6ehqbzcZPf/pTvvnNbwIw\nMjLCvn37+M1vfsOpU6fu5iUpFArFtkVF1BUKhUIhCQQCLCwskJubC4DT6cTj8SQ421lZWRw/fpy2\ntjYALl26xNzcXMI+paWl1NXVyX0UCoVCsXaUo65QKBQKyXe/+10OHz5MU1MTAG63G4DCwsKE/Ww2\nm3zN7Xaza9cu8vLyEvYpLCzE4/FswqgVCoViZ7L7bg8gFdPT03d7CAqFQrEpZGdn3+0hSF544QXa\n2tpobW1Fp9OtuP9q9klG2XeFQnGvsBH2XUXUFQqFQsGf//mf87Of/YyzZ89SUVEht9vtdoAlkXGP\nxyNfs9vt3Lp1i/Hx8YR93G633EehUCgUa0c56gqFQnGP893vflc66fv37094rbKyErvdTktLi9wW\ni8VobW3loYceAqCxsZE9e/Yk7DMyMkJPT4/cR6FQKBRrZ0umvmjZSsvCCoVCsRFspfSP5557jrff\nfptf/vKXZGdny7xzi8WCyWRCp9Px/PPP8+qrr1JbW0tNTQ3f+973sFgsfOtb3wIW7fR3vvMdXnzx\nRWw2G3v37uWFF17g0KFDPPbYY2nPrez71iYQCHDlyhUADh8+jNVqXfP7A4EAVqtVvld7zOrqavr7\n+5c9fiAQAFjzuTcD7dhu516t9/4qtiYbbd+3vKOuUCgUijvHj3/8Y3Q6HSdPnkzY/sorr/BXf/VX\nALz44otEo1Gee+45JicnOXbsGC0tLZhMJrn/G2+8we7du3n66aeJRqM89thjvP3227eVx67YGlit\nVg4fPiz/fTus5IivZgx3m3QPC+sZmzjmeu+vYuezJXXUtU8jKuKiUCh2GveyjbuXr/1eQ0SMw+Ew\nDQ0NlJaWyu3weSRa/Hsrspao92qvRUXSdzYbbeNURF2hUCgUCsWGY7Vaqa6upqOjg/7+fpkCo3VM\ntf8Wjq74v3DstwvK4VbcCZSjrlAoFAqFIiXrjXhbrdaEFKnlznPlyhW8Xi8DAwNkZWXx1FNP3XVn\nfSPSf1JRXV295KFFoUiFctQVCoVCoVAsYSNSNO6Uo7uZaAthtX/fDsn3VKFYCeWoKxQ7kK2e96lQ\nKO4dVmOHtA79Vkx9WetDy2pssLLTitWgHHWFYoehCpUUCsVGsNnRcHGO7W6zlrPB4p4GAoF1K+Io\n7g2Uo65QKBQKhSIlyoFcJPmh5Xaj4SqKrlgrSp5RodiBqMlga3Mv27h7+doVS9kutmqtzY1SXVe6\n/PStfu2KtbHRNi5j3Uf4//zN3/wNR44cITs7G5vNxpNPPsn169eX7PfKK69QUlKC0Wjk0Ucfpaur\na6OGoFAo/j9KTUChUGwVRIfSVNuvXLnClStXUr6+VbidcSbb4FT3QNlpxWrYMEf9ww8/5E//9E9p\nb2/n7Nmz7N69m8cee4zJyUm5z/e//31ef/11/vZv/5YLFy5gs9l4/PHHCYVCGzUMheKeJd1kqFAo\nFHeLVE6u1laFw2HC4fDdHOKaEWkwq80tF/egv7+f6upqlZOuWBMblqP+29/+NuHvt956i+zsbNra\n2vjqV79KPB7njTfe4KWXXuIb3/gGAG+++SY2m4133nmHZ599dqOGolDcc6gCUoVCsR3Q2qrq6uo1\nvxc2P1UkVVHtenTllX1WrIU7VkwaCARYWFggNzcXAKfTicfj4dSpU3KfrKwsjh8/Tltbm3LUFQqF\nQqHYYaQrwhQkN0NK54zf7WDEes63E7TkFXePO+aof/e73+Xw4cM0NTUB4Ha7ASgsLEzYz2azMTo6\neqeGoVDcE6iJQKFQ3E2Wi3ankicU/9ZGqEdGRujo6MBkMu24lcGddC2KzeWOOOovvPACbW1ttLa2\notPpVtx/NfsoFIrl2Q4TwXZReFAoFKtnrdHuZMddHKOjo4O+vj5qampSvkcFIxT3IhtWTCr48z//\nc372s59x9uxZKioq5Ha73Q6Ax+NJ2N/j8cjXFArFzmW7KDwoFIqNYy1F7iaTiZqaGhoaGtJG5pWT\nrrjX2NCI+ne/+11+/vOf8/vf/579+/cnvFZZWYndbqelpYXGxkYAYrEYra2tvPbaaxs5DIViy6Ei\nyQqFYqeSLtqdqnA0nbO9XSPmyrYr7jQb5qg/99xzvP322/zyl78kOztb5qRbLBZMJhM6nY7nn3+e\nV199ldraWmpqavje976HxWLhW9/61kYNQ7HD2Y5G8W4XQW0VtutErFAoVkabwqL9GxYlGNvb2xkd\nHaW4uJimpiZKS0vTHmO7oGy7YjPYMEf9xz/+MTqdjpMnTyZsf+WVV/irv/orAF588UWi0SjPPfcc\nk5OTHDt2jJaWliVV3wpFKpRR3P6oz0yh2LmkstGHDx+W+eexWIz+/n5MJtM9n8ayHYNOirvDhjnq\nCwsLq9rv5Zdf5uWXX96o0yoUW56tEElWk4JCobgbCIfcarVSWVmJ0+ncMcG5VLZ9NbZWBZ0Ua+GO\nyTMqFBvNVnB4b5e7Od7bmRSUY69QKNbKSjbaYrHIQtF0tmW72Z50OfmrsbXhcJhAILBtrlVxd9hw\n1ReF4k5yry+Xbga3q86yFnUHhUKxM0llowOBAK2trfziF7+go6Mjra0YGRmhtbV1yytDrdfWWa1W\nWVzb39+/pa9VcfdREXWFYoezXJRro6JXailXoVAsRzgcJhaLEQ6HUzY1WklHfauwnK1bSyqM1Wrd\nMSlAijuLctQVinuAVI5zuglnO6QYbbflcYXiXiP5N2oymaiqqsLhcOD1elO+ZyUd9e2A1kEPBAL0\n9/cDq3PqFYpUKEddoWBrO353Y2xrPddmTjoqeq9QbG2Sf6Ow6ISbTCZKSkooKSkBlnYo3Q6Oa/I4\nU9lncf3hcBggbeR8K1+nYuugHHXFPc9Wdvzu5Ng2emLcSvdNoVDcPQKBAF6vVzqoq7U1W9WGJDvj\n2qj5cvbZZDJRXV2taqsU60I56grFFmIrR/a3Ctsl8qZQ3IuIXPOBgQGqqqrk9u36W72dYImyUYqN\nRDnqinuejTSq63G0U00I6xnbSmPZyisJK7GdxqpQKBZJlyaSvG2l92wVlrPPW3G8iu2JctQVW47t\nkJOdijvl+Io8yLXo7W5nJ1yhUGxfrFYrDQ0NhMPhhNzsVDZpJTu1FezYSsGSVNtGRkYAKC0tvbOD\nU9wTKEddsaXYCob5TrHSA0i6CeFOPgCo5VmFQrGRCDtns9nu8kg2jrXYx5GREX7xi18A8NRTTyln\nXbFulKOu2HZs1aXQlfTK0znb2uvZbCd8Pefbqp+Dlu0wRoVip6C1c8lFlKls0mqi1SqYoLjXUY66\nYkuxkmHe6OjyRjtyt5NHvtL1rOaepHptubGs97q3w8rHdhijQrFTSRV4uJ2Ga3frd7taG6ndT6Qn\nPvXUU4BKfVFsDMpRV2w5Nsswb5YjJwz5eiJDG1kQqhxYhUJxJxBBBWHzlnN2t7IdWu3YklcQ0jU3\nUijWg3LUFduKtS6FblTEPN1x1quscjtLuxtxTaIRx+1itVqprq5OOY6tkm6ils0VirtDf39/QrOf\n5RzXcDi8pkL5zSR5bKuxbeFwGJfLteJ+CsVqybjbA1AoVotW+WQtkeMrV65IA6tFOHKpcsa1+6c7\njnb7yMhIynOshrXkpiePJd01pHrf7Y4v3fH6+/vp7+9f1b26W6hGI6vjo48+4sknn6S0tJSMjAze\nfPPNhNe//e1vk5GRkfDfQw89lLDPzMwMf/Znf0ZBQQFms5mvf/3r0mFRKJLRPuwn25G7TaqxpbJt\nWvtbWlpKdXU14XCYlpYWWltbt9Q1KbYvKqKu2BasZinydiK5qRz0tS5fhsNhOjo6UkaONmMFYDU5\nlKnaea/3vJvFVh7bTiEcDtPQ0MAzzzzDn/zJn6DT6RJe1+l0PP7447z11ltyW2ZmZsI+zz//PO+/\n/z7vvvsue/fu5YUXXuBrX/saly5dIiNDxYTuFZLTX5Z7WLZarUts0VZhtWNLtvfbybYqtgfKUVfs\nCIQzKhyO0tLSNTnJ2vdDoiOb7jjaCUk496nYSO3z9aRziGVcEfkRx1hOqUGMK/l8K92T2xlfKrZy\nHutO4oknnuCJJ54AFqPnycTjcTIzM9NK7k1PT/OTn/yEn/70p5w8eRKAt956i3379Line truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x4707b38>"
|
||||
"<matplotlib.figure.Figure at 0x13c7710>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
@@ -344,29 +344,25 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"On the left I created 3000 points normally distributed based on the Gaussian\n",
|
||||
"The left plot shows 3,000 points normally distributed based on the Gaussian\n",
|
||||
"\n",
|
||||
"$$\\mu = \\begin{bmatrix}0\\\\0\\end{bmatrix},\\, \\, \\, \\Sigma = \\begin{bmatrix}32&15\\\\15&40\\end{bmatrix}$$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"and passed them through the system\n",
|
||||
"The right plots shows these points passed through this set of equations:\n",
|
||||
"\n",
|
||||
"$$\\begin{aligned}x&=x+y\\\\\n",
|
||||
"y &= 0.1x^2 + y^2\\end{aligned}$$ \n",
|
||||
"\n",
|
||||
"I then showed how poorly the EKF did at estimating the new mean and covariance compared to the UKF.\n",
|
||||
"Using a finite number of randomly sampled points to compute a result is called a *Monte Carlo* (MC) method. The idea is simple. Generate enough points to get a representative sample of the problem, run the points through the system you are modeling, and then compute the results on the transformed points. \n",
|
||||
"\n",
|
||||
"This technique of using a finite number of randomly sampled to compute a result is called a *Monte Carlo* (MC) method. The idea is simple. Generate enough points to get a representative sample of the problem, run the points through the system you are modeling, and then compute the results on the transformed points. \n",
|
||||
"\n",
|
||||
"In a nutshell this is what particle filtering is. A bit later I'll demonstrate how MC can integrate over a probability distribution, but you don't need that formality to understand how they work. It's very simple. It is the Bayesian filter algorithm we have been using throughout the book applied to thousands of particles, where each particle represents a *possible* state for the system. We extract the estimated state from the thousands of particles using weighted statistics of the particles.\n",
|
||||
"In a nutshell this is what particle filtering does. The Bayesian filter algorithm we have been using throughout the book is applied to thousands of particles, where each particle represents a *possible* state for the system. We extract the estimated state from the thousands of particles using weighted statistics of the particles.\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"**Generic Particle Filter Algorithm**\n",
|
||||
"## Generic Particle Filter Algorithm\n",
|
||||
"\n",
|
||||
"1. **Randomly generate a bunch of particles**\n",
|
||||
" \n",
|
||||
" Particles can have position, heading, and/or whatever other state variable you need to estimate.\n",
|
||||
" Each has a weight indicating how likely it matches the actual state of the system.\n",
|
||||
" Particles can have position, heading, and/or whatever other state variable you need to estimate. Each has a weight (probability) indicating how likely it matches the actual state of the system.\n",
|
||||
" \n",
|
||||
"2. **Predict next state of the particles**\n",
|
||||
"\n",
|
||||
@@ -374,18 +370,17 @@
|
||||
" \n",
|
||||
"3. **Update**\n",
|
||||
"\n",
|
||||
" Update the weighting of the particles based on a measurement.\n",
|
||||
" Update the weighting of the particles based on measurements.\n",
|
||||
" \n",
|
||||
"4. **Resample**\n",
|
||||
"\n",
|
||||
" Discard highly improbable particle and replace them with copies of more probable particles\n",
|
||||
" Discard highly improbable particle and replace them with copies of more probable particles.\n",
|
||||
" \n",
|
||||
" Optionally, compute mean and covariance of the set of particles to get the most likely current state.\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" \n",
|
||||
" \n",
|
||||
"This naive algorithm runs into some practical difficulties which we will need to overcome, but this is the general idea. To demonstrate this I wrote a particle filter that performs robot localization. This is the same robot localization problem used in the UKF and EKF chapter. In this problem we have a robot that we can control. It has sensors that measure its distance to various landmarks, and we can steer the robot. Both the sensors and control mechanism have noise in them, and we need to estimate the robot's position.\n",
|
||||
"5. **Compute Estimate**\n",
|
||||
"\n",
|
||||
" Optionally, compute weighted mean and covariance of the set of particles to get a state estimate.\n",
|
||||
" \n",
|
||||
"This naive algorithm has practical difficulties which we will need to overcome, but this is the general idea. Let's see an example. I wrote a particle filter for the robot localization problem from the UKF and EKF chapters. The robot has steering and velocity control inputs. It has sensors that measures distance to visible landmarks. Both the sensors and control mechanism have noise in them, and we need to track the robot's position.\n",
|
||||
"\n",
|
||||
"Here I run a particle filter and plotted the positions of the particles. The plot on the left is after one iteration, and on the right is after 10. The red 'X' shows the actual position of the robot, and the large circle is the computed weighted mean position."
|
||||
]
|
||||
@@ -416,7 +411,21 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"After the first iteration the particles are still largely randomly scattered around the map, but you can see that some have already collected near the robot's position. The computed mean is quite close to the robot's position. This is because each particle is weighted based on how closely it matches the measurement. The measurement (with noise) says the robot is near (1,1), so particles that are near (1, 1) will have a high weight, and particles at, say, (15, 19.5) will have a very low weight. The estimated is computed as the weighted mean of all of the particle's positions, so the estimate is quite accurate. This is partially the result of the random number generation - other runs could produce a poor estimate if no particles were particularly near the robot.\n",
|
||||
"If you are viewing this in a browser, this animation shows the entire sequence:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<img src='animations/particle_filter_anim.gif'>"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"After the first iteration the particles are still largely randomly scattered around the map, but you can see that some have already collected near the robot's position. The computed mean is quite close to the robot's position. This is because each particle is weighted based on how closely it matches the measurement. The robot is near (1,1), so particles that are near (1, 1) will have a high weight because they closely match the measurements. Particles that are far from the robot will not match the measurements, and thus have a very low weight. The estimated position is computed as the weighted mean of positions of the particles. Particles near the robot contribute more to the computation so the estimate is quite accurate.\n",
|
||||
"\n",
|
||||
"Several iterations later you can see that all the particles have clustered around the robot. This is due to the *resampling* step. Resampling discards particles that are very improbable (very low weight) and replaces them with particles with higher probability. There are multiple algorithms for this, and we will discuss them later. \n",
|
||||
"\n",
|
||||
@@ -428,15 +437,15 @@
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Probability distributions via MC\n",
|
||||
"Suppose we want to know the area under the curve $y= \\mathrm{e}^{\\sin(x)}$ in the interval [0, $\\pi$]. We can compute the area with the definite integral $\\int_0^\\pi \\mathrm{e}^{\\sin(x)}\\, \\mathrm{d}x$. As an exercise, go ahead and compute the answer; I'll wait. If you are wise you did not take that challenge; $y= \\mathrm{e}^{\\sin(x)}$ cannot be integrated analytically.\n",
|
||||
"Suppose we want to know the area under the curve $y= \\mathrm{e}^{\\sin(x)}$ in the interval [0, $\\pi$]. Area is computed with the definite integral $\\int_0^\\pi \\mathrm{e}^{\\sin(x)}\\, \\mathrm{d}x$. As an exercise, go ahead and find the answer; I'll wait. \n",
|
||||
"\n",
|
||||
"The world is filled with equations which we cannot integrate. We cannot integrate something as simple as $\\int_a^b x^x\\, \\mathrm{d}x$. For a realistic example from physics, think of how we calculate luminosity of an object. An object reflects light. Some of that light bounces off of other objects and restrikes the original object, increasing the luminosity. This creates a *recursive integral*. Good luck with that one.\n",
|
||||
"If you are wise you did not take that challenge; $\\mathrm{e}^{\\sin(x)}$ cannot be integrated analytically. The world is filled with equations which we cannot integrate. For a realistic example from physics, think of how we calculate the luminosity of an object. An object reflects some of the light that strike it. Some of the reflected light bounces off of other objects and restrikes the original object, increasing the luminosity. This creates a *recursive integral*. Good luck with that one.\n",
|
||||
"\n",
|
||||
"However, these are trivial to compute using a Monte Carlo technique. Create a bounding box that contains the curve in the desired interval, generate random pairs $(x,y)$, and count how many fall under the curve. Multiply the area of the bounding box by the ratio of points that were under the curve vs the total number of points and you will have computed the are under the curve. As you tend towards infinite points you can achieve any arbitrary precision. In practice, a few thousand points will give you a fairly accurate result.\n",
|
||||
"However, integrals are trivial to compute using a Monte Carlo technique. To find the area under a curve create a bounding box that contains the curve in the desired interval, generate random pairs $(x,y)$, and count how many fall under the curve. Multiply the area of the bounding box by the ratio of points that were under the curve vs the total number of points and you will have computed the area under the curve. As you tend towards infinite points you can achieve any arbitrary precision. In practice, a few thousand points will give you a fairly accurate result.\n",
|
||||
"\n",
|
||||
"Think of how powerful this technique is. You can use it to numerically integrate any function. The function can be of any arbitrary difficulty, including being non integrable and noncontinuous, and yet we can trivially find the answer. \n",
|
||||
"You can use this technique to numerically integrate a function of any arbitrary difficulty. this includes non-integrable and noncontinuous functions. This technique was invented by Stanley Ulam at Los Alamos National Laboratory to allow him to perform computations for nuclear reactions which were unsolvable on paper.\n",
|
||||
"\n",
|
||||
"Let's use it to compute $\\pi$ by computing the area of a circle. This is very easy. We will create a circle with radius 1, and bound it in a square. The means that the length of a side of the box will be 2, and hence the area of the box is 4. We will generate a set of uniformly distributed random points within the box, and count how many fall inside the circle. The area of the circle is then the area of the box times the ratio of points inside the circle over the total number number of points. Finally, we know that $A = \\pi r^2$, so $\\pi = A / r^2$.\n",
|
||||
"Let's compute $\\pi$ by finding the area of a circle. We will define a circle with a radius of 1, and bound it in a square. The side of the square has length 2, so the area is 4. We generate a set of uniformly distributed random points within the box, and count how many fall inside the circle. The area of the circle is computed as the area of the box times the ratio of points inside the circle vs. the total number of points. Finally, we know that $A = \\pi r^2$, so we compute $\\pi = A / r^2$.\n",
|
||||
"\n",
|
||||
"We start by creating the points.\n",
|
||||
"\n",
|
||||
@@ -445,19 +454,19 @@
|
||||
"ps = uniform(-1, 1, (N, 2))\n",
|
||||
"```\n",
|
||||
"\n",
|
||||
"A point is inside the circle if it's distance from the center point (0, 0) is less than or equal to one. We can compute the distance by using `numpy.linalg.norm`. If you aren't familiar with this, the norm is the magnitude of a vector. Since vectors start at (0, 0) calling norm will compute the point's distance from the origin.\n",
|
||||
"A point is inside a circle if its distance from the center of the circle is less than or equal to the radius. We compute the distance with `numpy.linalg.norm`. If you aren't familiar with this, the norm is the magnitude of a vector. Since vectors start at (0, 0) calling norm will compute the point's distance from the origin.\n",
|
||||
"\n",
|
||||
"```python\n",
|
||||
"dist = np.linalg.norm(ps, axis=1)\n",
|
||||
"```\n",
|
||||
"\n",
|
||||
"Next we compute which of this distances fit the criteria with this code, which returns a bool array that contains `True` if it meets the condition\n",
|
||||
"Next we compute which of this distances fit the criteria. This code returns a bool array that contains `True` if it meets the condition `dist <= 1`:\n",
|
||||
"\n",
|
||||
"```python\n",
|
||||
"in_circle = dist <= 1\n",
|
||||
"```\n",
|
||||
"\n",
|
||||
"All that is left is to count the points inside the circle, compute pi, and plot the results. I've put it all in one cell so you can experiment with alternative values for `N`."
|
||||
"All that is left is to count the points inside the circle, compute pi, and plot the results. I've put it all in one cell so you can experiment with alternative values for `N`, the number of points."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -492,13 +501,15 @@
|
||||
"from numpy.random import uniform \n",
|
||||
"\n",
|
||||
"N = 20000\n",
|
||||
"radius = 1\n",
|
||||
"area = (2*radius)**2\n",
|
||||
"\n",
|
||||
"ps = uniform(-1, 1, (N, 2))\n",
|
||||
"dist = np.linalg.norm(ps, axis=1)\n",
|
||||
"in_circle = dist <= 1\n",
|
||||
"\n",
|
||||
"in_circle_count = np.count_nonzero(in_circle)\n",
|
||||
"area = 2*2\n",
|
||||
"pi = area*(in_circle_count / N)\n",
|
||||
"pts_in_circle = np.count_nonzero(in_circle)\n",
|
||||
"pi = area * (pts_in_circle / N)\n",
|
||||
"\n",
|
||||
"plt.scatter(ps[in_circle,0], ps[in_circle,1], \n",
|
||||
" marker=',', edgecolor='k', s=1)\n",
|
||||
@@ -514,15 +525,38 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I show this to you to help illustrate that you really are computing a numerical integration. An alternative way to find the area of a circle is by integrating the infinitely many circles that lie inside the circle:\n",
|
||||
"I show this to you to help illustrate that you are performing a numerical integration. An alternative way to find the area of a circle is by integrating the infinitely many circles that lie inside the circle:\n",
|
||||
"\n",
|
||||
"$$A = \\int_0^{2\\pi r} 2\\pi t\\, dt = \\pi r^2 $$\n",
|
||||
"\n",
|
||||
"The Monte Carlo technique does something very similar by summing the infinite number of points inside the circle. Computers can't sum infinities, so we settle for a representative sampling of points and less than infinite precision in the results. Using Monte Carlo for this problem is remarkably unnecessary, but imagine instead trying to compute the recursive integral for luminosity that I mentioned. As far as I know there is no analytic solution to that problem, yet we could use the code above to find the answer. MC is a powerful tool to reach for when analysis fails. \n",
|
||||
"The Monte Carlo technique does something very similar by summing the infinite number of points inside the circle. Computers can't sum infinities, so we settle for a representative sampling of points and less than infinite precision in the results. Using Monte Carlo for this trivial problem is overkill, but imagine trying to compute the recursive integral for luminosity. As far as I know there is no analytic solution to that problem, yet we could use the code above to find the answer. \n",
|
||||
"\n",
|
||||
"This insight leads us to the realization that we can use Monte Carlo to compute the probability density of any probability distribution. For example, suppose we have a Gaussian. The Gaussian has a probability distribution as shown below."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from filterpy.stats import plot_gaussian_pdf\n",
|
||||
"plot_gaussian_pdf(mean=2, variance=3);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The *probability density function* (PDF) is the probability that the random value falls between 2 values. For example, we may want to know the probability of x being between 0 and 2 in the graph above. This is a continuous function, so we need to take the integral to find the area under the curve, as the area is equal to the probability for that range of values to occur. \n",
|
||||
" \n",
|
||||
"$$P[a \\le X \\le b] = \\int_a^b f_X(x) \\, dx$$\n",
|
||||
"\n",
|
||||
"It is easy to compute this integral for a Gaussian. But real life is not so easy. For example, the plot below shows a possible and realistic probability distribution. There is no way to analytically describe an arbitrary curve, let alone integrate it analytically."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
@@ -532,53 +566,9 @@
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAADBCAYAAADGrth2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczXX///EnxjLGICXZkkJUspdCSfuOIoQWRBSVtqtN\nqUjyvdJySYtUSqWir0pX+4VScSF79mRJyjYmY7bz++P5m+8pITEzn3PO53G/3c5tzkyWFw2e8zqv\n9+tdJBKJRAQAAACEQNGgCwAAAAAKC+EXAAAAoUH4BQAAQGgQfgEAABAahF8AAACEBuEXAAAAoUH4\nBQAAQGgQfgEAABAahF8AAACEBuEXAAAAoUH4BQAAQGgQfgEAABAahF8AAACEBuEXAAAAoUH4BQAA\nQGgQfgEAABAahF8AAACEBuEXAAAAoUH4BQAAQGgQfgEAABAahF8AAACERlLQBQAAAOzNjz9KX3zh\nx/r1UmqqVKaMVK6cVLy4VKqUVLasdPzx0jnnSMWKBV0xYh3hFwAAxIRIRFqyRJo+3Y9p06RVq/b/\n+1erJl17rdS1q1S7dsHVifhWJBKJRIIuAgAAhNvcudLAgdJnn/3x42XLSq1aSa1bO9BmZEjbt/uR\nnS3t3Clt3ixNmSItXx79fvXqSZdcInXsKDVqJBUpUqi/HMQwwi8AAAjMunXSvfdKY8e681uunHTu\nuVKLFlLLllKDBvs3yhCJeDTixRelyZOlrVuj/61BA+npp/1jAoRfAABQ6LZskR55RHriCXdzk5Kk\nG25wEK5Q4eB+7Kwsj028/bb05pvSpk3u/N56qzR4sOeEEV6EXwAAUGgyM6VHH5VGjIh2Zy+/XBoy\npGDmdDMypAcfdNDOzZXq1nWX+eST8//nQnwg/AIAgEKxcaOD7vTpfr9NG4fSZs0K/uf+5hvpmmuk\nxYulokXdZR4yREpJKfifG7GF8AsAAArc1187+K5bJ1WpIr38snTmmYVbQ0aGNGiQu845OVKdOtJr\nr0lNmhRuHQgW4RcAABSYSEQaNUq66SbP4rZoIb31lnTEEcHVNGeO1K2btHChdwWPHu2uMMKBG94A\nAECB2LFDuvJKqV8/B98bb/QqsyCDr+TVZzNnevQhK8u7gR97LNiaUHjo/AIAgHy3eLF02WV+m5Ii\nPfec1Llz0FX92VNPOZRL3gRx773B1oOCR/gFAAD56osvpIsuktLTfdnE22/7bax6+WWPPeTmSvfd\nJz3wQNAVoSAx9gAAAPLNxx9LF1zg4Nuxo/Ttt7EdfCWpe3dp3DhvgRg82A8kLsIvAADIF5MmueO7\nc6fUo4c3KZQpE3RV+6dz52gAHjRIGjYs6IpQUAi/AADgoL35pleZZWZK/ftLzz67f9cSx5LOnaWX\nXvJtcHfeKY0cGXRFKAiEXwAAcFBee83BMSdHuusu6fHH3UGNR127evWZ5PVsY8YEWw/yHwfeAADA\nAXvxRalnz+hhsfvvd+c03uVtgShRQvryS6lp06ArQn6J06/LAABA0IYP947c3FzpwQe9JSERgq/k\nHcDXX+8xjg4dpM2bg64I+YXOLwAA+NuGDJHuvtvPR470nG+i2bVLatlSmjVLuuQSH+hLlHAfZoRf\nAADwt4weLfXp4yD40ku+KjhRrV4tNWwobdsmjRgh3XJL0BXhYBF+AQDAfhs/3lcWRyIOwdddF3RF\nBW/iRKl9e2+vmDZNOuWUoCvCwWDmFwAA7JdJk9zljUSkhx4KR/CVpHbt3PHNyZG6dHEXGPGLzi8A\nAPhLU6ZIbdv6ANg99/iAW5hkZrrjO3u2b4R76aWgK8KBIvwCAIB9mj1batVK+u03acAA6Z//DOfB\nryVLpEaNpIyM6G5jxB/GHgAAwF4tXCidf76Db7du4Q2+klS3rn/9ktS7tw/DIf7Q+QUAAHu0cKF0\nxhnSpk3S2WdLkydLJUsGXVWwIhHpsst8CO7ccz0OEtYvBuIV4RcAAPzJhg2+1Wz9eoe8iROl5OSg\nq4oNGzdK9epJW7ZIb7whdewYdEX4Oxh7AAAAf5CZKV1xhYNvq1YE391VqiQNHernt9zikRDED8Iv\nAAD4P5GIr/WdNk2qUkWaMIHguye9evnw27p10iOPBF0N/g7GHgAAwP/Ju7Y4OVn64gvppJOCrih2\nTZ/uzniJEtKCBVLt2kFXhP1B5xcAAEiS3nvPO3yLFPFNbgTffWvZUrrmGo+JcO1x/KDzCwAAtHy5\nD7ht2yY9/LB0111BVxQffvrJHd8dO6QPP/ThQMQ2Or8AAIRcWppvb9u2zVf5/uMfQVcUP444Qrrv\nPj/v14/Db/GA8AsAQIhFItK113qnb7160tix7K39u266SapfX1qxQho2LOhq8FcYewAAIMRGjnR4\nS02VZs2S6tQJuqL4NG2adNppUunS0uLF0pFHBl0R9obOLwAAIfXtt9Jtt/n5mDEE34PRqpXUoYPH\nHgYODLoa7AudXwAAQmjrVu+pXb1a6t/fHWAcnLVr/QXEzp1eg9aiRdAVYU/o/AIAEDKRiNSnj4Nv\nkybS8OFBV5QYqlWLdn1vuknKzQ22HuwZnV8AAEJm7Fjvp01JkebM4XKG/JSe7u7v+vXeldypU9AV\nYXd0fgEACJFly6QbbvDzp58m+Oa3lBTpgQf8/J57fAEGYgvhFwCAkMjKkrp2dXfyiiuk7t2Drigx\nXX21dOyxXn323HNBV4PdMfYAAEBI3Huv9NBDXsP13XdS+fJBV5S4Jk6U2reXKlaUli7l9zqW0PkF\nACAEvvpKGjLEF1i88gphrKC1bSu1bClt2uTfd8QOOr8AACS49HSpQQO/DH/HHdIjjwRdUTjMmiU1\nayaVLOkb9I45JuiKINH5BQAg4d15p4PviSdGD2Oh4DVtKnXrJu3aJQ0Y4BVzCB6dXwAAEtgXX0hn\nnCElJUkzZ0oNGwZdUbj89JMPv23fLr3zjtSuXdAVgc4vAAAJKj1d6tHDz++5h+AbhCOOkB5+2M/7\n9/f1xwgW4RcAgAQ1aJC0cqXnfe+6K+hqwuv6632V9Nq13KYXCxh7AAAgAc2cKTVv7ufffOP5UwRn\n6lTp9NOl5GRp8WKpRo2gKwovOr8AACSY7GypZ08pN1e6+WaCbyw47TRfdbxzpzRwYNDVhBudXwAA\nEszw4dLtt0tHHeUVW6VLB10RJOnHH6W6dT33++mnUps2QVcUTnR+AQBIICtXetZXkkaNIvjGkurV\no7PXN97o66ZR+Ai/AAAkkP79/dJ6ly7SeecFXQ12N3CgVKuWtGiR9PjjQVcTTow9AACQIN57T7r4\nYik1VVq61Gu2EHs+/FA6/3x35Rct4vBbYaPzCwBAAvjtN6lfPz9/4AGCbyw77zypQwf/P+vbl5vf\nChvhFwCABDBsmLRmjS+y6N8/6GrwV0aOlMqVkz74QJowIehqwoWxBwAA4tycOdLJJ/sA1dSpUqtW\nQVeE/TF6tNSnj1SpkscfKlQIuqJwoPMLAEAcy8mRrr3WwbdfP4JvPOnVy/t/N270PmYUDjq/AADE\nsSeekAYM8BqtJUtYbRZvli2TTjxRysiQpkxhQ0dhoPMLAECcWrpUuvNOP3/ySYJvPKpdWxo82M/7\n9JG2bw+2njAg/AIAEIdyc6XrrvNO365dpUsvDboiHKibb5YaN5Z++EG66aagq0l8jD0AABCHnnlG\nuv56qWJFafFi6dBDg64IB2PxYqlRI2nXLsYfChqdXwAA4szGjdIdd/j5008TfBNBvXrR8YdevaTN\nm4OtJ5ERfgEAiCORiLc6bN8uXXCBdPnlQVeE/HLLLdIpp0hr17L9oSAx9gAAIRKJeEY0O1sqXlxK\nSvKjSJGgK8P+mjRJatdOKlNGmj9fOuqooCtCflq6VGrQwNsfPvjA1yAjfxF+ASAORCLSr79Kq1ZJ\nGzZI69b57a+/ugOYk+NAu3OnlJbml0zT0x1qs7J8jWpmpt9mZ//xxy5aVEpOlooV87cvXdpzpOXK\nSeXL+23ZslJqqlS5sk+n16rl0FWiRCC/HaG1YYNvcPv5Z684u/HGoCtCQXjsMem227y+bv58/xlE\n/iH8AkAMycz0rtb58/2YM8d7QNev90GY/FCqlLu9WVkOwjk5B/bjFC0qVa3qEFytmnTccT6x3qiR\nQzLyV2am1Lq1NGOG1KaN9PHH/n+AxJOdLZ16qjRzptSpk/Taa7w6k58IvwAQoC1bpE8+kaZPl775\nxmE3M3PP37ZsWQfN6tUdLqtU8UGn8uUdZvM6uKmpviY1JcXfLynJ3dySJf12925tXsc4LwRv3+6O\n8rZt0tatfpuW5o+vWSMtXy6tWOHnubl7rrV6dalZM3cpmzXzrWN59eDA3HGH9Oij/kLjv/+VDj88\n6IpQkL7/XmraVNqxQ3r9demKK4KuKHEQfgGgEEUi0rx50vvve55vxow/B8hataT69f1o0EA64QR3\nWGMtPGZm+mDO6tUOwt995/A+Z86fF/UnJ3t108UXe4bxiCMCKTlu/e//eo9vsWLS1KnuCiLxPfus\n1Lu3VKmSXxEqXz7oihID4RcAClhOjru7Eyc69K5dG/1vSUlSy5bSmWdKJ5/sLmm8/wOXm+t/qGfN\nciCePl369ts/fptTT5XatpUuu0w6+uhg6owXS5e6A5iWJj3ySHTFGRJfbq50+un+M9Srl8MwDh7h\nFwAKyI8/Si+8II0Z4+d5Kld29/PCC6WzzvI4Q6Jbu1Z69113uz/7zCfZ85xyirtbHTu6Q4yo7dv9\n+7NokdShg/TGG8x+hs2iRR4fysqSPvxQOvfcoCuKf4RfAMhHubk+iPSvf0nvvRcdaTjmGF9Be/HF\nPhQW5gCzY4dvsJo4UZo82e9LnlXu3t07bOvVC7bGWJCb6+745Mn+/fj663B8oYQ/GLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.axes._subplots.AxesSubplot at 0x82837f0>"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAADaCAYAAADnoPfrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVPX+P/DXzLANCMi+b7IpoqgsKuKuKJqZpqiZpall\nmbm0eO1n9+Y3r7eya6aldr/1LVMzSk3LcMFd1FREEAUBRTZlRjbZmYFhfn9go3NBQAUOy+v5ePAA\n3udzZt6Hcxjec/ic9xGp1Wo1iIiIiIioVYiFToCIiIiIqDNhAU5ERERE1IpYgBMRERERtSIW4ERE\nRERErYgFOBERERFRK2IBTkRERETUiliAExERERG1oiYV4Js2bYKbmxukUikCAgIQHR39yLEnTpzA\nxIkTYW9vDyMjI/j5+eG7776rM0YsFtf5SElJebqtISIiIiJq4xotwCMiIrBkyRKsXLkScXFxCA4O\nRlhYGLKysuodf+7cOfj5+WH37t24du0aXn/9dbz66qvYuXNnnbGJiYmQyWSaDw8Pj6ffIiIiIiKi\nNkzU2J0w+/fvjz59+uDrr7/WxLy8vDBlyhSsWbOmSU8ybdo0qFQq7Nq1C0DtGfARI0YgNzcXFhYW\nT5E+EREREVH70uAZcKVSidjYWISGhmrFQ0NDcfbs2SY/SVFREczNzevEAwICYG9vj1GjRuHEiRNN\nfjwiIiIiovZKp6GFeXl5UKlUsLGx0YpbW1tDJpM16Qn279+PY8eOaRXs9vb22LJlCwIDA6FQKLBt\n2zaMHDkSJ0+eREhIyBNsBhERERFR+9BgAf60zpw5g5kzZ2Ljxo0ICAjQxL28vODl5aX5fsCAAUhP\nT8fatWu1CvCioqKWTI+IiIiIqEWZmprWiTU4BcXS0hISiQRyuVwrLpfLYWdn1+CTRUdHY9y4cfjo\no4/w2muvNZpcUFAQUlNTGx1HRERERNSeNViA6+npwd/fH4cPH9aKR0VFITg4+JHrnTp1CuPGjcOq\nVavw1ltvNSmRuLg42NvbN2ksEREREVF71egUlGXLlmHWrFkICgpCcHAwtmzZAplMhgULFgAAVqxY\ngYsXL+LIkSMAajucjB8/Hm+++SZmzJihmSsukUhgZWUFAFi/fj3c3Nzg4+MDpVKJ7du3Y9++fdiz\nZ88j86jv9D09nZiYGADQmh5EnQf3P/EYIB4DnRv3f8tpbBp1owV4eHg48vPzsXr1auTk5KBXr16I\njIyEk5MTAEAmkyEtLU0zfuvWraisrMTatWuxdu1aTdzV1VUzrqqqCu+++y6ys7MhlUrh6+uLyMhI\njB079ok2koiIiIiovWi0D7iQHn73wDPgzY/vfDs37n/iMUA8Bjo37v+W01gN26Rb0RMRERERUfNg\nAU5ERERE1IpYgBMRERERtSIW4ERERERErYgFOBERERFRK2IBTkRERETUiliAExERERG1okZvxENE\nRO2HqkaFguK7yL13B/dK81FcVoiS8iIUlxeiorIUSpUSVVUKVFUrUV5RDpFIjMNJRhCLJdCR6MJQ\n3whGBsYwkhrD0MAYZl0sYWFqA3MTa5h1sYREwj8bRERPi6+kRETtVFFpATLv3kDW3ZvIzr2FuwXZ\nyCuWo6ZG9ViPU1JZ0KRxIpEYVl3tYG/pAgdLN9hbusDZxgOmRuZPkj4RUafFApyIqB1Qq9XIvZeD\n1OwEpGZfxc3b11BU1rTCuflyqMHdwtu4W3gbcalnNXELUxu42/ugm30PeDj4wqqrHUQiUavmRkTU\nnrAAJyJqo5TVCqRkXsHVWxdwLT0WRaX5TVrP1MgcVmb2sDC2hrGRGYwNTWFiaAYjA2Po6epDV0cP\nejr6SLyWCDXU6OHTA6qaalRVV6G8sgTlilKUVZSgtKIIBcW5yC+WI79YjuKywnqfL79IjvwiOS4k\nHQcAWJrawse1H3xc/eHh6As9Hf1m+5kQEXUELMCJiNoQhbICV9LOIy71LK5nxqGqWvnIsXo6+nC0\n6gYnG3c4WbvD3tIFVl3toa9r0KTnypLKAAB2Fs5NGq+sUiAnPxN38tJxOy8dt3NvIVN+A1Uq7Rzz\nimQ4FR+JU/GR0NPRR0+3APT1HAQfN38W40REYAFORCQ4laoa1zPjEHP9JBLSLkBZrah3nIGeIdwd\nfODp2Auejr5wsHSFWCxptTz1dPXhYusJF1tPTaxaVYWsuzeRdicJN25fQ2r2VSirKjXLldUKXE49\ng8upZ6Cna4BeboEI8hkBb2c/iEVsxEVEnRMLcCIigRQU5+LctcM4d/UIisvrn95hY+6IXm5B8O0W\nBBdbT0haseBuCh2JLtzsusPNrjtG+k9CVXUV0u4kIjH9Eq6lX8LdwtuascqqSlxKOY1LKadhbmyF\n/j4j0d9nJMxNrATcAiKi1scCnIioFdWoa5CUHoszCYdwLf0S1OqaOmPsLJzh7z0EfTyCYW1mL0CW\nT05XRxfezn7wdvbDpCGv4E5eRu0Z8JRo3L13RzOuoCQXB87/hIPnI9CzWyCG950ADwdfXrxJRJ0C\nC3AiolZQVa3ExesncOzSXq1C9C8mRmYI7D4MAd5D4WDl2voJthB7SxfYW7pg3IAZuJ13C+cTj+Hi\n9ZMorywBAKihxtW0C7iadgGOVt0wrO8E9PMKgY5EV+DMiYhaDgtwIqIWVK4oxZkrh3Aybn+900y8\nnf0Q0mssfN0CO/RNbkQiERytusFxaDc8O+hlJKSdx7mrUUjOiteMyc5Nw/bDX+C36B8wwn8iBvUa\n2+QLSomI2pOO+2pPRCSgckUpjsfuw4m4/VAoK7SWSfUMMdB3NIJ9x7S7KSbNQVdHF/28QtDPKwSy\ngiycvLwfF64f13R8KS4vxN7T3+NIzK8Y0W8iBvcOg76eVOCsiYiaDwtwIqJmpFBW4ETcfhyL3YsK\nRZnWMtMuFhjedwIG9gyFVN9QoAzbFltzJ0wb+TqeCZ6JM1cP43R8pOYGQ6UVRfjtzA84eulXjOj3\nHIb2eQZ6umxjSETtHwtwIqJmUFWtxOkrBxAVsxtlFcVay2zMHDHSfxICug/h3OZHMJKaIDRwCob3\nnYjziUcRFbMbhSW5AICyyhL8fnYbTsX/gbABM9DfZ0Sb6wZDRPQ4WIATET0FtVqNuBtn8Vv0D8gv\nlmstszK1Q9iA6ejnFdKq/brbM10dXYT0HosBPUfiQtJxHL64CwXFdwEARWUF+OnoVzhx+TdMGDQL\nvm6B7JpCRO0SC3AioieUKb+BPae+RdqdJK24mbEVxvafhqAew3mm9gnpSHQR7BuK/j1G4M/Eozhw\n/icUl9VexCoryML//r4GHg498fzQ+R2qawwRdQ4swImIHlNx2T38dmYrLiQd14obGhgjrP80BPuO\nga4Op5o0B4lEB4N6jUFA96E4cfk3HLn0q+ai1hu3r2HtzmUY7DcOYQOmw1C/i8DZEhE1DQtwIqIm\nqlHX4GzCYfx+5gdUKMs1cbFYgiF+4zE2KByGBiwCW4K+rgHGBIUj2DcUhy78jOiEQ6ipUaFGXYOT\ncftxKfk0Jgyahf4+I3iLeyJq81iAExE1QdbdNPx8bDMy5Kla8V7dgjAxZHanbCcoBGPDrpgy7FWE\n9A7DrhP/i5SsKwBqO6bsPPIlzl2NwoxRC2Fn4SxwpkREj9ak0wSbNm2Cm5sbpFIpAgICEB0d/cix\nJ06cwMSJE2Fvbw8jIyP4+fnhu+++qzPu5MmT8Pf3h1Qqhbu7O77++usn3woiohZSqazAnpPf4rOf\n3tEqvi1NbfH6c//A/Anvs/gWgK25ExZOWoVXxr0Hsy6Wmni6LBmf/rgMkX/uRFV1lYAZEhE9WqMF\neEREBJYsWYKVK1ciLi4OwcHBCAsLQ1ZWVr3jz507Bz8/P+zevRvXrl3D66+/jldffRU7d+7UjLl1\n6xbGjRuHkJAQxMXFYcWKFVi0aBH27NnTfFtGRPSUkjPj8a/tb+FE3O9Qq2sA1M5JHtt/Gla8uAE9\nXPoKnGHnJhKJ0MczGO+/9CVCA6dCIq79p66qphoHz0fg051LkXbnusBZEhHVJVKr1eqGBvTv3x99\n+vTROkPt5eWFKVOmYM2aNU16kmnTpkGlUmHXrl0AgOXLl2Pv3r1ITk7WjJk/fz6uXbuGs2fPamJF\nRUWar01NTZu2RdRkMTExAICAgACBMyEhcP8/WoWiHPuiv8fZq4e14l5OvRE+/DVYmzkIlFnz6mjH\nQE5+FnYe/RLpOQ/+toggQkjvMDw7aBbvplmPjnYM0OPh/m85jdWwDZ4BVyqViI2NRWhoqFY8NDRU\nq1BuShLm5uaa78+dO1fvY8bExEClUjX5cYmImltyZjw+3rFYq/g2MjDGrDFLsXDSqg5TfHdEdhZO\nWDJlDaYMmw99XQMAgBpqnL4SiY9/XFKnXSQRkVAavAgzLy8PKpUKNjY2WnFra2vIZLImPcH+/ftx\n7NgxrYJdLpfXeUwbGxtUV1cjLy+vzjLgwbs0an782XZu3P+1qlRKxNyKQqr8slbc2dwb/d3DICo1\nwqVLlwTKrmV1tGPAEDYY7zcf528ewO3CGwCA/CI51v+yAj0dgtHHeYhmugrV6mjHAD0e7v/m5+np\n2eDyFn0FOnPmDGbOnImNGzfy3xtE1GblldzB6ZRfUVJZqInp60gR1G0sXC19eLfFdqiLvilG9JiG\nW7lXcT7tIKpUCgDAtdtncafwBkK8JsLMqO7JHiKi1tBgAW5paQmJRAK5XPv2ynK5HHZ2dg0+cHR0\nNMaPH4+PPvoIr732mtYyW1vbOmfQ5XI5dHR0YGlpifqwgG9+nPvVuXH/AzU1Khy59CsOXt2JmpoH\n0996uw9A+PAFMDHqKmB2La8zHAOBCMTokmewI2qjpmVhYfldRF75DuMHvoAR/s916Line truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x7eb4780>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from filterpy.stats import plot_gaussian_pdf\n",
|
||||
"plot_gaussian_pdf(mean=2, variance=3)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The *probability density function* (PDF) is the probability that the random value falls between 2 values. For example, we may want to know the probability of x being between 0 and 2 in the graph above. This is a continuous function, so we need to take the integral to find the area under the curve, as the area is equal to the probability for that range of values to occur. \n",
|
||||
" \n",
|
||||
"$$P[a \\le X \\le b] = \\int_a^b f_X(x) \\, dx$$\n",
|
||||
"\n",
|
||||
"It is easy to compute this integral for a Gaussian. But real life is not so easy. For example, the plot below shows a possible and realistic probability distribution. There is no way to analytically describe the curve, let alone integrate it symbolically."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAADBCAYAAADGrth2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8DnT6//H3cRxkpxRlKTtFIZElMhOmBY3SNy3apk0y\nTZs0JdVMU1NKy5QWk2kZU31VUlFCypqiw0k4SISsx3Y46/374/3zPZkknOVzL6/n43E/3Ec4l5Lz\nvq/7+lyfpEgkEhEAAACQAEqFLgAAAAAoKYRfAAAAJAzCLwAAABIG4RcAAAAJg/ALAACAhEH4BQAA\nQMIg/AIAgKC2bpUuv1xKSip4VKvm73vrLWnLlp//nNxcacUK6eOPpaefltq3L/i59epJ//qXlJ9f\n8r8XRL8k9vwCAIBQ5syRLrxQWr1aKldOuv56qU8fqUMHKSXl4H+dSET68ENpyBBp4UJ/X6tW0sMP\nS7/9rUMxIBF+AQBAAJGI9Nxz0uDBUk6O1K6d9MorUqNGhft18/Kk116Thg6VfvjB39etmzRmjFS7\nduHrRuwj/AIAgBK1bZs0cKBDqiQNGiQ9+qhUpkzRfY7MTOnJJ935zciQ6tSRJk+WGjcuus+B2ET4\nBQAAJWbyZOmKK9yVLV9eeuEFqX//4vt8mzdLvXpJM2dKNWtK06ZJTZoU3+dD9CP8AgCAEjFunHTR\nRT6sdtppHkVo2rT4P++uXdJ550lTp0pHH+1Dci1bFv/nRXRi2wMAACh2H35YEHz/9Cd3Yksi+EpS\nhQrSe+9JZ50lbdggde0qzZ1bMp8b0YfOLwAAKFazZvnQ2Z49Dr6PPhpm+8KePVK/fg7ClSo5kHfs\nWPJ1ICzCLwAAKDbLljlgbtwoXX21Z3xDrh3LyfH+4LFj3RGeNIkAnGgIvwAAoFh8953UubO0Zo3U\ns6c7rqVLh67K69CuuEJ69VWpcmXpiy/YApFICL8AAKDIbdwodekiLV7szurEiVLFiqGrKpCb6xGI\nt9/24bc5c3zJBuIfB94AAECRysjw4bLFi6UTT5QmTIiu4Cu5Az1mjNSwoZSaKt19d+iKUFLo/AIA\ngCITiUh9+7qj2rix9Omn3q8brebO9VXKeXnSRx85tCO+0fkFAABF5qmnHHwrV5Y++CC6g6/kfcPD\nhvn5pZdKa9eGrQfFj84vAAAoErNnS2ec4Y0Kb73lDnAsyMuTuneXpkzxSrbJk8NupEDxovMLAAAK\nbeNG6cILHXwHD46d4CtJycnS669LRx3lADxyZOiKUJzo/AIAgELJz5fOPdeXRnTo4GuEy5QJXdWh\ne+cd6fzzpZQUacYMqW3b0BWhOND5BQAAhfLEEw6+1av78ohYDL6S1KePdNNN7l737eurkBF/6PwC\nAIDDNmuW53xzc9057d07dEWFk5Ulde3q+eXzz5fGjQtdEYoa4RcAAByW9eul1q2ldeukW26RRowI\nXVHRWL1aat5c2rnTl3P06BG6IhQlxh4AAMAhy8uTLrrIwbdzZ+nhh0NXVHTq1JHuucfPb7nFXW3E\nD8IvAAA4ZI89Jk2fLtWqJb3xhg+JxZPBg6UGDXxL3bPPhq4GRYmxBwAAcEgWLpROPVXKzvZFFr/7\nXeiKise77/oQXNWq0rffSsccE7oiFAU6vwAA4KBlZUmXXebge+218Rt8JalXL6lnTykjQ7r99tDV\noKjQ+QUAAAftzjulRx6R6teXvv5aqlgxdEXFa/lyH37LzvaYR+fOoStCYdH5BQAAB2XqVOnvf5dK\nlZJeeSX+g6/kud8hQ/z8mmukPXvC1oPCI/wCAIBftWmTdOmlUiQi/fnPvsktUQwdKjVrJi1dKj34\nYOhqUFiMPQAAgAOKRHzhw7vvSh07StOmSaVLh66qZM2cKXXqJCUnS/PmSSefHLoiHC46vwAA4IBG\njXLwrVJFevXVxAu+kjvdN97onb9XXsnu31hG5xcAAPyitDSvNduzRxo71hdbJKqdO6WTTpJWrfLu\n3+uvD10RDgfhFwAA7Nfu3Q6+33wjDRggvfxy6IrCe/NNqV8/6aijpGXLvAMYsYWxBwAAsF9Dhjj4\nNm0qPfNM6GqiwwUXSGec4QOAw4eHrgaHg84vAAD4malTpW7dPN87e7bUpk3oiqLHggX+95GU5F3H\nJ54YuiIcCjq/AABgHzt3Sldf7ef33EPw/W+nnOLb7fLypEGDvA0DsYPOLwAA2MfAgdI//uGQN3eu\nlJISuqLos3mz1KSJv331VemSS0JXhINF+AUAAP/no4+kHj087sA+2wMbPdod8ho1pCVLpGrVQleE\ng8HYAwAAkCRt3ChdfrmfDx9O8P01V17pw28bN0p33RW6GhwsOr8AAECRiNS7t/Tee1KXLtInn/g2\nMxxYWprUqpWUk+Ob77p0CV0Rfg2dXwAAoH/8w8G3alVpzBiC78E68URp6FA/v/pqKTMzbD34dYRf\nAAAS3MKF0m23+fkLL0j16oWtJ9YMHSq1aCEtXy795S+hq8GvYewBAIAEtmOHb3FbulS66irppZdC\nVxSbZs2SOnZ0x3z+fF+DjOhE5xcAgAQ2eLCDb4sW0lNPha4mdp1+unTddVJurvSHP0j5+aErwi+h\n8wsAQIIaP96H3MqWlb78kpvKCmvbNqlZM2ndOnfQr7oqdEXYHzq/AAAkoIwM6frr/fyhhwi+RaFK\nFenvf/fz22+XNmwIWw/2j/ALAECCiUSkm292h/L00/0cRaN/f6l7d2nLFumOO0JXg/1h7AEAgATz\n4oueSy1f3re4NWsWuqL4kp7uTnp2tjRjhtShQ+iK8FN0fgEASCDffisNGuTnzz1H8C0ODRt67EGS\nBg70IThED8IvAAAJIjdXGjBA2rNHuuwyP1A87rrL+5IXLJD+9rfQ1eCnGHsAACBBDBsm3X+/VLu2\ntGiRD2ih+Hz8sed/y5SRvvlGatAgdEWQ6PwCAJAQZs6UHnhASkry9cUE3+J31lnS5Zd79veWW0JX\ng73o/AIAEOd275ZOOcWXWQwZ4tVmKBnr1klNmvgmvfffl84+O3RFoPMLAECcGz7cwbd5c+m++0JX\nk1hq1fK4ieTDb7t3h60HdH4BAIhrX30lnXaar9udOVNq3z50RYknN1dq3VpauFB68EHp7rtDV5TY\n6PwCABCncnKkq6+W8vJ8kQXBN4zSpaUnnvDzv/5VWrMmbD2JjvALAECcGjHCq7aOP176y19CV5PY\nunWT+vaVMjM5/BYaYw8AAMShFSukk07yjOnEiVKPHqErwvffe+561y5p3Djp/PNDV5SY6PwCABBn\nIhHphhscfPv3J/hGi7p1PfYgSTfd5A0QKHmEXwAA4sxrr0kffSRVqyY9/njoavBTN93kA4hr1/rC\nEZQ8xh4AAIgjW7ZITZtKGzdKo0dLV14ZuiL8ty+/lNq2lZKTvQGiadPQFSUWOr8AAMSRIUMcfLt0\nka64InQ12J82baRrrvEKtD/+0WMqKDl0fgEAiBMzZ0odO0opKdLXX0vNmoWuCL9kwwapcWNp2zbp\nzTelCy4IXVHioPMLAEAcyM31ITdJuuMOgm+0O/rogmumb7vNK9BQMgi/AADEgaeeklJTpRNOkIYO\nDV0NDsa110otW0qrVkl/+1voahIHYw8AAMS4deukJk28OmvCBOmcc0JXhIP1+edS585SmTLS4sVS\n/fqhK4p/dH4BAIhxt9zi4NurF8E31nTqJF12mZSdLd11V+hqEgOdXwAAYtgHHzjwli8vpaX5KmPE\nltWrffhtzx5p9mypXbvQFcU3Or8AAMSozExp4EA/f+ABgm+sqlPH3XtJuv12Vp8VNzq/AADEqLvv\n9nW5LVv64oTSpUNXhMO1bZvUoIG0ebP0/vvS2WeHrih+0fkFACAGpadLjz7q588+S/CNdVWqFMz8\n3n23lJ8ftp54RvgFACAG3XabD0kNGCB16BC6GhSFG2+UjjtOWrBAeuut0NXEL8YeAACIMZMmST17\nShUqSMuWSbVqha4IRWXUKOn666VGjaRvvqGjXxzo/AIAEENyc6Vbb/Xze+8l+Mabq66SGjb0i5pX\nXgldTXwi/AIAEENGjfJKsxNOkAYPDl0NilpKijRsmJ/ff79HW1C0CL8AAMSIrVvd7ZWkxx6TypYN\nWw+Kx8UXS82aSd99J738cuhq4g/hFwCAGPHww9KWLVLXrlKfPqGrQXFJTi54kfPQQ1JOTth64g0H\n3gAgzmRlSevWOSRlZUk7d0q7dknbt0ubNkkbNvgq3Lw8f1HNzvaPy8ryDVPJyVLFir4xrFw5H7hJ\nSZHKlJGqVfNFCscfLx17rOdNU1JC/44Tw+LFUuvW/m80d67Utm3oilCc8vKk5s2lpUulf/3LVyCj\naBB+ASCGZGZKy5dLK1ZI338vrVolrV0rrV9f8Ni6teTqSU6Watd2GD7pJF/L2rkzN40Vtdxc6fTT\npXnzpCuukP75z9AVoST8858+ANesmbRokVSK9+uLBOEXAKLMnj3SkiUOOgsWSCtXSj/+KK1Z43D7\na0qXlmrWlI46yjOhFSr4UamSVKOGdPTRUuXKDq57O7plyxY88vLcKd61y93g3Fw/srLcOV650qF7\n3To/9vdVpG5dqUsXB7YzznAHKymp6P9dJYqnnpJuvtnX4C5a5P9+iH85Od788P330ptvShdcELqi\n+ED4BYDA0tOl6dOlmTOl+fOlhQt/ecavdGmpfn1fg1q3rjusxx0nHXOMA2+tWtKRR5Zchygry1+Y\nV66UvvpKmjVL+vRTX9X6U3XqSP37u4vVuHHJ1BYv1q5152/7dumdd6TevUNXhJL0j39IAwd65GXe\nPF5EFgXCLwCUsJwc6fPPHWQ++MDh96eSkhwQW7aUTj1VatLE4fa44xxuo33pfV6eVLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x879c860>"
|
||||
"<matplotlib.figure.Figure at 0x9016550>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
@@ -593,7 +583,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"But we have already stated that we can use Monte Carlo methods to compute any integral. The PDF is computed with an integral, hence we can compute the PDF of this curve using Monte Carlo. "
|
||||
"We can use Monte Carlo methods to compute any integral. The PDF is computed with an integral, hence we can compute the PDF of this curve using Monte Carlo. "
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -602,23 +592,21 @@
|
||||
"source": [
|
||||
"## The Particle Filter\n",
|
||||
"\n",
|
||||
"All of this brings us to the particle filter. Consider tracking a robot or a car in an urban environment. For consistency I will use the robot localization problem from the ends of the EKF and UKF chapters. In this problem we tracked a robot that had a sensor that could detect the range and bearing to landmarks. \n",
|
||||
"All of this brings us to the particle filter. Consider tracking a robot or a car in an urban environment. For consistency I will use the robot localization problem from the EKF and UKF chapters. In this problem we tracked a robot that had a sensor measures the range and bearing to landmarks. \n",
|
||||
"\n",
|
||||
"Before I continue let me point out that the filter I am going to develop next is a specific form of particle filter that is intuitive to grasp and relates to the problems we have studied in this book. In other words, particle filters are a family of algorithms. A more academic approach would give you a better grasp of this family and the interrelationships between the algorithms, but as in previous chapters I choose to develop concepts based on examples and basic reasoning, and only once that is done go back and develop the theory. This will leave a few of the steps seeming a bit 'magical' since I haven't offered full explanation for their use. Full explanations will follow later in the chapter.\n",
|
||||
"Before I continue let me point out that particle filters are a family of algorithms. I'm presenting a specific form of a particle filter that is intuitive to grasp and relates to the problems we have studied in this book. A more academic approach would start by defining the characteristics and interrelationships between the algorithms. I choose to develop concepts based on examples and basic reasoning, and only once that is done go back and develop the theory. This will leave a few of the steps seeming a bit 'magical' since I haven't offered full explanation for their use. Full explanations will follow later in the chapter.\n",
|
||||
"\n",
|
||||
"Taking insight from the discussion in the previous section we start by creating several thousand *particles*. Each particle has a position that represents a possible belief of where the robot is in the scene, and perhaps a heading and velocity. Suppose that when we initialize the filter we have no knowledge of the location of the robot. We would want to scatter the particles uniformly over the entire scene. If there was a large clump of particles near a specific location that would imply that we were more certain that the robot is there. If you think of all of the particles representing a probability distribution, locations where there are more particles represent a higher belief, and locations with fewer particles represents a lower belief.\n",
|
||||
"\n",
|
||||
"Think back to the *Discrete Bayes* chapter. In that chapter we modeled positions in a hallway as discrete and uniformly spaced. This is very similar except the particles are randomly distributed.\n",
|
||||
"Taking insight from the discussion in the previous section we start by creating several thousand *particles*. Each particle has a position that represents a possible belief of where the robot is in the scene, and perhaps a heading and velocity. Suppose that we have no knowledge of the location of the robot. We would want to scatter the particles uniformly over the entire scene. If you think of all of the particles representing a probability distribution, locations where there are more particles represent a higher belief, and locations with fewer particles represents a lower belief. If there was a large clump of particles near a specific location that would imply that we were more certain that the robot is there.\n",
|
||||
"\n",
|
||||
"Each particle needs a weight - ideally the probability that it represents the true position of the robot. This probability is rarely computable, so we only require it be *proportional* to that probability, which is computable. At initialization we have no reason to favor one particle over another, so we would typically assign a weight of $1/n$, where $n$ is the number of particles. We did the same thing in the discrete Bayes chapter. When we initialized the filter we assigned a probability of 1/N to each hallway position. The reason for $1/n$ is simple - the sum of all probabilities must equal one, and of course $\\sum\\limits_{i=1}^n \\frac{1}{n} = 1$.\n",
|
||||
"\n",
|
||||
"We will be operating on many particles so it is imperative that we are as efficient in our programming as possible. That means storing the particles in an NumPy `array` and making full use of vectorization and ufuncs.\n",
|
||||
"The combination of particles and weights forms the *probability distribution* for our problem. Think back to the *Discrete Bayes* chapter. In that chapter we modeled positions in a hallway as discrete and uniformly spaced. This is very similar except the particles are randomly distributed in a continuous space rather than constrained to discrete locations.\n",
|
||||
"\n",
|
||||
"To track our robot we need to maintain states for x, y, and heading. We will store `N` particles in a `(N, 3)` shaped array. The three columns contain x, y, and heading, in that order.\n",
|
||||
"\n",
|
||||
"In this problem the robot can move on a plane of some arbitrary dimension, with the lower right corner at (0,0).\n",
|
||||
"\n",
|
||||
"Here is a partial implementation of this problem. The particles are initially randomly distributed over the space, and are given a random heading. "
|
||||
"Here is a partial implementation of this problem. The particles are initially randomly distributed over the space with a random heading."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -658,11 +646,11 @@
|
||||
"source": [
|
||||
"### Predict Step\n",
|
||||
"\n",
|
||||
"The predict step in the Bayes algorithm predicts the movement of the robot and to update our belief in the robot's position by that movement. How would we do that with particles? Each particle represents a possible position for the robot. Suppose we send a command to the robot move 0.1 meters while turning 0.0073 radians. We would want to update each particle by moving it using the same command. If we did that exactly we would soon run into a problem. The robot's controls are not perfect so it will not move exactly as commanded. Therefore we need to add noise to the particle's movements to have a reasonable chance of capturing the actual movement of the robot. If you do not model the uncertainty in the system the particle filter will not correctly model the probability distribution of our belief in the robot's position. \n",
|
||||
"The predict step in the Bayes algorithm uses the process model to update it's belief in the system state. How would we do that with particles? Each particle represents a possible position for the robot. Suppose we send a command to the robot to move 0.1 meters while turning by 0.007 radians. We could move each particle by this amount. If we did that we would soon run into a problem. The robot's controls are not perfect so it will not move exactly as commanded. Therefore we need to add noise to the particle's movements to have a reasonable chance of capturing the actual movement of the robot. If you do not model the uncertainty in the system the particle filter will not correctly model the probability distribution of our belief in the robot's position. \n",
|
||||
"\n",
|
||||
"If this wasn't a robot, but something you are tracking passively, then you would want to include things like velocity in the particle state and use that estimate to make the prediction. When using particle filters I try to minimize the dimensionality of the problem; if I add 2 random variables I will probably have to increase the number of particles to ensure I have enough samples to correctly sample from the probability distribution of each. This can quickly become intractable. Thus I've excluded velocity from this problem.\n",
|
||||
"\n",
|
||||
"I cannot easily modify a class definition in Jupyter Notebook, so I will show you the method in text. The entire class is in pf_internal, and we will be importing when we run the filter. \n",
|
||||
"The predict step for the class ParticleFilter is listed below: \n",
|
||||
"\n",
|
||||
"```python\n",
|
||||
"def predict(self, u, std, dt=1.):\n",
|
||||
|
||||
+447
-342
@@ -1,342 +1,447 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/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: 16pt;\n",
|
||||
" font-style: italic;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#c76c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 1.5em;\n",
|
||||
" display: inline;\n",
|
||||
" white-space: wrap;\n",
|
||||
" } \n",
|
||||
" h3 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h3 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 14pt;\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: 100;\n",
|
||||
" font-size: 14pt;\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: 200;\n",
|
||||
" font-style: normal;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" margin-bottom: 0em;\n",
|
||||
" margin-top: 0.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: 125%;\n",
|
||||
" font-size: 120%;\n",
|
||||
" width:740px;\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",
|
||||
" div.output_wrapper{\n",
|
||||
" margin-top:0.2em;\n",
|
||||
" margin-bottom:0.2em;\n",
|
||||
"}\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",
|
||||
" scale:85,\n",
|
||||
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
|
||||
" }\n",
|
||||
" });\n",
|
||||
"</script>\n"
|
||||
],
|
||||
"text/plain": [
|
||||
"<IPython.core.display.HTML object>"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from __future__ import division, print_function\n",
|
||||
"%matplotlib inline\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0,'..') # allow us to format the book\n",
|
||||
"sys.path.insert(0,'../code') \n",
|
||||
"# use same formatting as rest of book so that the plots are\n",
|
||||
"# consistant with that look and feel.\n",
|
||||
"import book_format\n",
|
||||
"book_format.load_style(directory='..')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD9CAYAAABA3unhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztfd+LXUdyf2ssWTJfs/4yEBCbyNpEEC7Lvo/zkIcQo2Wd\nPyCIsA9aMiGwYHFhH2YHQvIUTDBIoAhBBmUC3qwiyYmDSEb2TBQNyUOM3gaRYELQ3CtlQeTBsTX3\nx8w9Pre/D/fbM337VnVXdVefc8a+Bctad87p011dXV31qerqE1prreY0pznNCaCFujswpznNqbk0\nVxBzmtOcUJoriDnNaU4ozRXEnOY0J5TmCmJOc5oTSiexP3z55ZdV9mNOc5pTA+iNN96Y+vfcgpjT\nnOaE0lxBzGlOc8JJI/TFF18c/u/s2VJfu9bXrVahHzwY6KIodKezrzudfV2WJdaE1lrrTmdfnz1b\n6uXlyf8rpbVSWp89W+rd3aF+8GCgz54t9dmzpX7wYAC2V5blzHO7u8OZ9trtyW8bG/1gv2KpLMuZ\nsWN9sf/98GFvqo/23zqdffRbPv4Y3kLfMXO2uFiSvnP//tF7778/0K1WgT6fSlC/DU+pvMHGYdrh\nyGgVBMkN5W8p36OsLa2n17pLJAWxuzvU3e4wagBm0hcXJ4LndpjKHPc5lwH37vX06upAt9tDvbVF\nVxCcycGY7v6+sdHXW1v9w3/fv9/TrVahldJ6cXGycCkTF+oftKBcxWQUUug73e5Qt9tDvby8rxcX\ny+DiTBFqKh9Dfaa02XSi9pvLb46yTVYQsZ0sy4mV8PDhZIG0WoV++LCnd3eHIpNn+rOz0z9cgEpp\n3WoVemenT1I6HKGyld3y8r5utyeK0+6LrcDMv7vdaauB2r8QFUUxtfPbisgIBfU7HF5ILEZMlmIV\nT6r1URdR+k3ht8u3yhUEVyjc5z/6qKefPOnpoihwbkWSzQzODs0Vqk5nX7daxZQldP9+77B9n9Cn\n7Iy+NlutQrfbEyU8Go2SFi51cTZxMWJ9ilE4OUx+br85z9jyZTbhp08HU1ZsVheDOpDQoNrtYRbT\nz2YQx8fnjqksj3x8SBB9izNWULE261wQTVQQEK+Kopj6bWurr3d3/a6y2879+3k2NV+/Ka6kzW/M\njd/Y6JOggWwKwjCb6h8boDJFmEKm6c5On7xwYnb2p08HWQA2iHxtSn0vxjK03cYm+fs+M5tqWWLA\nb87xhRQ0FayGAgEUmcjiYrhAHLRj2s+///6ABH75iLJLd7vDKb8c2kns97gg5dZWf0rQTMQEEqxU\nnMGnBIpiYkq228PD6FLMd7B+UwT1/v2e7nZl8KQcZI+Nungwy5crs9Juiq89My+Q9by7Owy2nQWk\ndIE3DFyR3G18CwbywwwYKrXbhkBKN6qSOmZsTK7CSzGDId6srfXASFATXQsfUV1PW65d4Pfq1b5e\nXR2gShP77sbG7CaSk6DN8dq1PimiJ6YgbOIIi5RfHGtyU/vqi0SY/261Cr28vK+Xl/dncgWwqAp3\nIbkC6048FK2QsspMvsu1a/1D5WdIEveoikzfut3h1KL1hVZHo5F++LCnV1YG+ubN3sw7IYLyYuyd\nPCe/uOFqrTMpCIlQF0ZY2zGgHbWvFBdqNBrpe/eOBObePXjnTtlpob5SErFScZ3Hj/d0uz2cSqza\n2ekH++Zz3yRJYlFxNp2yLA/xLMhi9BGEgxle5lw32HhCrm4WBaF1vsy10GLHQMqUKALme04zGo5g\nQHyJFQJsgt3f7IxJCSEL7Xr22CAgkJIbwiFXtnIkE2kN89sg/zs7fb20NELD2hh1u8MpnMq2xnK7\naTGubjYFgXXq/v1JzkNsQlQsE4tikoC0s9Nn++QhBdFqFXpr62VWtwrrR7c7m5IunUoc6zcbt8td\nRCmWhYu9bG+/nLFuILwr5nvue1tb/Sk+fPgh350zYHa7PZzJ7K0Cx+G6utkVBLa4KIIGLSQzaXYS\nUGjBp5puPhfDCOnKykBfvz6dQk2JfHAWMrZQq/DzY7EiKDcEsnrMzhxq37ZKbMXji4SFFl4oCmD+\nJuXOxVq6kkRVRrUoiOXlfdIuizHLRZOpadApmhkCKV2A8OrVvr5yZUjyR31Yyu7uUO/s9GfChL7d\nh9rvmLGmELSoXAXRatHn07QHhSaxZLtUDMrXjrQ7VxWoSx135S4GdooQ819TIg+GcplumHVEMcEx\nd8G2ENxQFHccMTuS9C5GAS+xDFSoLZNrAoUmqTkaVACbOpY6FXAKUfqRXUGYjuzuDvWjR3tToSHb\nRKag87EKAvIlQ2m1FIL6geUJUN6FTG87EccH+lG/EVKMKeY4Rpir6DPdoX7a479yZajX1+m7N2QB\ncnzxlPG771flRkhQJQrC0AQo7OnHj/d0p+PX4hD45uIQ3J2x08Fj3jHk9oOTEUpViraCsHdQV8lC\nJK0gcgk3tV23b7EnX2PQfEqbHAwlJxCZ2kebKlMQMWYeFTziTKj0BJXlZOc3SPri4mRRh7LrIKzB\nBSGhbDdst/UBulIuRirvUudTSkFh1lsVVkFdCiKWd5UpiDp2Jm4/Ysn032QZhsbhG68PpMT6v7s7\nRMFLSZAyZQeXmmMJ/11aBurAMXL20aZGKAitm4fexrSLnRZ1KUVAof4/ezYdYr1+fTYVWoLsb7da\nxVTmaEjQ6zStfeOQkIEYTMzwyQVsqcewc/fRULKCoC7qOsAZTOlggFUO0DJGQYT64/4dyuLc2ekF\n+5tiXUDKcHv7JerLN0lB2OOQsKxSZBviy8rKgJzjQ6XaXIwYoDC3lWC+RU3BdTM9Y48puxgCBiD6\n8jhiJjJGQaQqbEiw2+0hCtQeN/TepVD/JTGxO3emD99J8akWkFIS7JNUHim7eWx1KxNl8CUy2XgF\ntEvEALa+/P5U/vjGikVwsOI/OTaIqjadXBaQy8f1dbnTuBLkUxCV3YsxHo/V5uaBeuutU+qtt06p\nzc0DNR6Pwee63QPV7cJ/j/nul1+W6tKlkVpcPGqv1zuhLl8+rZ4/L1jtPX9eqB/+8Iy6enXyvx/+\n8MxMG8+fF+ry5dPqs89OqqtXz6jf+73X1C9+UXrb1Vp7+fPtb59U3/nOWF26NFKXLo3Ud74zVt/+\nNnoxmggtLCyoixdPq08+2VeXLo3Un/7pafX55xORef11rdbXD9S5c6dm3jl//rQ6f/60WlhIFy+q\n3NRNPrk1fPz000L927+N1K/+qlZffHFMrqShaBUJkzEU4eh09g+Reur3KCahm+FpDhbFVrei7DIU\n/MHtd6gAT6fjr0UBkVTVKUkXjUtV4hqxLhL3Pe4xgtxUGUjpI59ZbZvjd+/u6QsXiqlkKg6YF/qm\nD2SjEEUY3GegrE6336GFQFkoPuQ8tfhqldiSTSnRg9QQKzU0GaPE6uInRJVmUmKELSyzM9qn9m7d\n6unFxZJ0wIerIHyFdjljCbVhnqFmdXKtodDf4ROWPTT3QpIkhZ+zO0uCpJy2qrZypBVLIxSE1vDg\nOp198FCOAfe4prpUaEqKOMLDDX36voPx1LhanNvHOH3NwXPqopBcqFBbWMJYVXKW6zuNURAQlSW2\n2/WDB3y4ZncdphwFj8jhwrnWl3vC1j3/wemDT1DrzIXIrSBWVgZRlqwU5eJtoxWE1jhogwmimQxq\nVmOd5FtMuU1i40NDpxmNgojpg09Q61QQZSlXTdpt69q1vr5+vZ98dUMKfWMVhNa8jEhfKnDdgA9E\nGPDls5BidiQfD7EDYjFCF4pI1eXWlSW/2I6PoArREpc/xdLXysXIZXJBprTERbhVEAVIjN3VKd+G\nDojFIvA5Mg5TSXqHhdrLdX0klb4WIGUuTVeW08eum+pWYBTCCezITlVmOuTeUcJ7dWM7EEnzzZXj\npt8kBhFlnkTyIDhMySHg7mSZpKemuhUQQXyB8jyqUhCGp9AN4ZyCwU2hHBuTWym9iYoRIyo/RA5r\npZ5eSxVwLOwU6lOTJpQ6YVX58RhPoTL2dfOOSjlzMNyS+E3fnKjrUOSwFmeR5xDwHL5yboKElQPG\n+v4tQZiC8N1jmUpNUtghouSWNNm9bayC0FpeEKTDcbmJ018K6JcDtDQl/V0Xg1KBOjbKUnfiGoeO\nu4Ko3cWoejfgfu+4JO2Ens0Jvi0tjfSHH06DlKPRyJviHrvQzTi4d11SxyQti8fdxdB6OuUfq/Iu\nDlLau0+TmVUUBemyXYhSBS5VQdjgJXTSE7o3M6Zv0OU0ZszY+GMVVqcDX9OXGinJaZlU4erlphB/\nRMOc5mPHwdwyAsk5Iq21jMCluBjuTrWx0debm9Np0ylJQBQFQX2fM/dlied+hHjSlMNSpm/d7vDw\nOEBViiJWOYX4I6ogzMdiBKtqStnpJMbGmVD7Wfgqu95UVl9Kgpi9+NxsVEp6cooC5V6gI+GiSVJZ\n8u4tkfxuLM9rURDuxaopxUhymWypvnJdyg+LLpjfFhdLUul9H9m+6eZmj11OP3beqHPCmYNcLgY0\nxk4HPn2cWz6o/MAiZ1ldDNv/2t2doN2tVqFbrUnFolgzqwpU+zii7dD37QIwkgKKCV6Vfj32DOf7\nRjal6l1g32+ygsD6HHKJRKIYrl8skXZa907tIw5AlsMCwnaCTkf2FCs2B02YG66LJpEWHeJxk10M\naM6ePh0EgwkieRA5tGYd4JJP4LiLvQ5LQzqC5NspfXPTNHQf6i/3YJXNC1P7AYocNRWkhHiwvf0y\nuMYaqyCqXGChb8X0JaeCC/mSkHtn3tndHbJubqL4rfaOHIq61BHyhubCPZrNWWA//elA3727d1jw\n175mIGYjqUJ5QnO2ujqr6MQVRE4hqIp5ocUcs9hzKYiYXd19J/Ywmz0fRVFMYU6mLLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x5a2cda0>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from random import uniform\n",
|
||||
"import matplotlib.pylab as pylab\n",
|
||||
"from matplotlib import animation\n",
|
||||
"from gif_animate import animate\n",
|
||||
"import matplotlib as mpl\n",
|
||||
"from math import pi\n",
|
||||
"\n",
|
||||
"def make_particle():\n",
|
||||
" return [uniform(0,100), uniform(0,100), uniform(0, 2), uniform(0, 2*pi)]\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"def move(particles, t=1.):\n",
|
||||
" vx = particles[:, 2]*np.cos(particles[:, 3])\n",
|
||||
" vy = particles[:, 2]*np.sin(particles[:, 3])\n",
|
||||
" particles[:, 0] = (particles[:, 0] + vx*t) % 100\n",
|
||||
" particles[:, 1] = (particles[:, 1] + vy*t) % 100\n",
|
||||
" \n",
|
||||
"particles = np.asarray([make_particle() for i in range(1000)])\n",
|
||||
"\n",
|
||||
"def panimate(frame):\n",
|
||||
" move(particles, t=1.)\n",
|
||||
"\n",
|
||||
" plt.cla()\n",
|
||||
" plt.scatter(particles[:, 0], particles[:, 1])\n",
|
||||
" plt.xlim(0, 100)\n",
|
||||
" plt.ylim(0, 100)\n",
|
||||
" plt.gca().get_xaxis().set_ticks([])\n",
|
||||
" plt.gca().get_yaxis().set_ticks([])\n",
|
||||
"\n",
|
||||
" \n",
|
||||
"animate('13_particle_move.gif', panimate, \n",
|
||||
" frames=500, interval=50, figsize=(4, 4))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"source": [
|
||||
"<img src='13_particle_move.gif'>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.4.3"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/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",
|
||||
"@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
|
||||
"\n",
|
||||
".CodeMirror pre {\n",
|
||||
" font-family: 'Source Code Pro', Consolas, monocco, monospace;\n",
|
||||
"}\n",
|
||||
" div.cell{\n",
|
||||
" width: 900px;\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",
|
||||
" text-align: left;\n",
|
||||
" } \n",
|
||||
" h2 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" text-align: left;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h2 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 16pt;\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: wrap;\n",
|
||||
" text-align: left;\n",
|
||||
" } \n",
|
||||
" h3 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h3 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-size: 14pt;\n",
|
||||
" line-height: 100%;\n",
|
||||
" color:#d77c0c;\n",
|
||||
" margin-bottom: 0.5em;\n",
|
||||
" margin-top: 2em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: wrap;\n",
|
||||
" text-align: left;\n",
|
||||
" }\n",
|
||||
" h4 {\n",
|
||||
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
|
||||
" }\n",
|
||||
" .text_cell_render h4 {\n",
|
||||
" font-weight: 100;\n",
|
||||
" font-size: 14pt;\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",
|
||||
"\n",
|
||||
" .text_cell_render h5 {\n",
|
||||
" font-weight: 200;\n",
|
||||
" font-style: normal;\n",
|
||||
" color: #1d3b84;\n",
|
||||
" font-size: 16pt;\n",
|
||||
" margin-bottom: 0em;\n",
|
||||
" margin-top: 0.5em;\n",
|
||||
" display: block;\n",
|
||||
" white-space: nowrap;\n",
|
||||
" }\n",
|
||||
" div.text_cell_render{\n",
|
||||
" font-family: 'Fira sans', verdana,arial,sans-serif;\n",
|
||||
" line-height: 150%;\n",
|
||||
" font-size: 110%;\n",
|
||||
" font-weight: 400;\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",
|
||||
" div.output_wrapper{\n",
|
||||
" margin-top:0.2em;\n",
|
||||
" margin-bottom:0.2em;\n",
|
||||
"}\n",
|
||||
"\n",
|
||||
" code{\n",
|
||||
" font-size: 6pt;\n",
|
||||
"\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",
|
||||
" equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
|
||||
" },\n",
|
||||
" tex2jax: {\n",
|
||||
" inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
|
||||
" displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
|
||||
" },\n",
|
||||
" displayAlign: 'center', // Change this to 'center' to center equations.\n",
|
||||
" \"HTML-CSS\": {\n",
|
||||
" scale:95,\n",
|
||||
" availableFonts: [],\n",
|
||||
" preferredFont:null,\n",
|
||||
" webFont: \"TeX\",\n",
|
||||
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
|
||||
" }\n",
|
||||
" });\n",
|
||||
"</script>\n"
|
||||
],
|
||||
"text/plain": [
|
||||
"<IPython.core.display.HTML object>"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from __future__ import division, print_function\n",
|
||||
"%matplotlib inline\n",
|
||||
"import sys\n",
|
||||
"sys.path.insert(0,'..') # allow us to format the book\n",
|
||||
"sys.path.insert(0,'../code') \n",
|
||||
"# use same formattibng as rest of book so that the plots are\n",
|
||||
"# consistant with that look and feel.\n",
|
||||
"import book_format\n",
|
||||
"book_format.load_style(directory='..')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import numpy as np\n",
|
||||
"from numpy.random import randn, random, uniform, seed\n",
|
||||
"import scipy.stats\n",
|
||||
"\n",
|
||||
"class ParticleFilter(object):\n",
|
||||
"\n",
|
||||
" def __init__(self, N, x_dim, y_dim):\n",
|
||||
" self.particles = np.empty((N, 3)) # x, y, heading\n",
|
||||
" self.N = N\n",
|
||||
" self.x_dim = x_dim\n",
|
||||
" self.y_dim = y_dim\n",
|
||||
"\n",
|
||||
" # distribute particles randomly with uniform weight\n",
|
||||
" self.weights = np.empty(N)\n",
|
||||
" self.weights.fill(1./N)\n",
|
||||
" self.particles[:, 0] = uniform(0, x_dim, size=N)\n",
|
||||
" self.particles[:, 1] = uniform(0, y_dim, size=N)\n",
|
||||
" self.particles[:, 2] = uniform(0, 2*np.pi, size=N)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" def predict(self, u, std):\n",
|
||||
" \"\"\" move according to control input u with noise std\"\"\"\n",
|
||||
"\n",
|
||||
" self.particles[:, 2] += u[0] + randn(self.N) * std[0]\n",
|
||||
" self.particles[:, 2] %= 2 * np.pi\n",
|
||||
"\n",
|
||||
" d = u[1] + randn(self.N)\n",
|
||||
" self.particles[:, 0] += np.cos(self.particles[:, 2]) * d\n",
|
||||
" self.particles[:, 1] += np.sin(self.particles[:, 2]) * d\n",
|
||||
"\n",
|
||||
" self.particles[:, 0:2] += u + randn(self.N, 2) * std\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" def weight(self, z, var):\n",
|
||||
" dist = np.sqrt((self.particles[:, 0] - z[0])**2 +\n",
|
||||
" (self.particles[:, 1] - z[1])**2)\n",
|
||||
"\n",
|
||||
" # simplification assumes variance is invariant to world projection\n",
|
||||
" n = scipy.stats.norm(0, np.sqrt(var))\n",
|
||||
" prob = n.pdf(dist)\n",
|
||||
"\n",
|
||||
" # particles far from a measurement will give us 0.0 for a probability\n",
|
||||
" # due to floating point limits. Once we hit zero we can never recover,\n",
|
||||
" # so add some small nonzero value to all points.\n",
|
||||
" prob += 1.e-12\n",
|
||||
" self.weights += prob\n",
|
||||
" self.weights /= sum(self.weights) # normalize\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" def neff(self):\n",
|
||||
" return 1. / np.sum(np.square(self.weights))\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" def resample(self):\n",
|
||||
" p = np.zeros((self.N, 3))\n",
|
||||
" w = np.zeros(self.N)\n",
|
||||
"\n",
|
||||
" cumsum = np.cumsum(self.weights)\n",
|
||||
" for i in range(self.N):\n",
|
||||
" index = np.searchsorted(cumsum, random())\n",
|
||||
" p[i] = self.particles[index]\n",
|
||||
" w[i] = self.weights[index]\n",
|
||||
"\n",
|
||||
" self.particles = p\n",
|
||||
" self.weights = w / np.sum(w)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" def estimate(self):\n",
|
||||
" \"\"\" returns mean and variance \"\"\"\n",
|
||||
" pos = self.particles[:, 0:2]\n",
|
||||
" mu = np.average(pos, weights=self.weights, axis=0)\n",
|
||||
" var = np.average((pos - mu)**2, weights=self.weights, axis=0)\n",
|
||||
"\n",
|
||||
" return mu, var"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 30,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAETCAYAAAD6a4mGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXtcVNX6/z8z3GFsQBFBxyugHLAGDD2KIaOmHe8GkXY5\nhfrLc1JE0JNlmXdTzvmaCYiWmVanTA0sj5euXBwMTY/DlHq4itqUiKiDjgJyWb8/aO/mxnUPl8Hn\n/XrxYth77bWfPcz6zFrPetazRIwxBoIgCAGIO9oAgiCsHxISgiAEQ0JCEIRgSEgIghAMCQlBEIIh\nISEIQjAkJITVo1AoIBbTR7kj6dLvfl1dHXbt2oWxY8eiR48esLe3h4eHBx5++GFERUVh3759BuX3\n7NkDsViMNWvWdJDFjaPRaLBhwwZERkbCx8cHYrEYYrEY+fn5DV7zwQcfYObMmfDx8YFUKoVEIoG/\nvz/mz5/f6HUthWvMjf209n3l6r5y5YrZ8yKRCCKRSIj5FicjIwNisRhz5szpaFPaBduONqCtqKur\nw/Tp03H06FG4urpi6tSpkMlkqK6uxvnz55GSkoIff/wRs2bNMrm2s30oOc6cOYM333wTYrEYAwcO\nhKurK8rLyxu95pNPPkFJSQlGjRoFT09PiMVinDt3Drt378ZHH32EgwcPYtKkSRazMSoqCgMGDDB7\nTqFQtLrexv4nH330ESoqKlpdd1vSWT9LlqbLCslnn32Go0ePIjAwEJmZmejWrZvB+fv37yMrK8vs\ntZ012Hf48OFQKpWQy+WQSCRQKBQ4fvx4o9ccPXoUDg4OJse/++47TJw4EbGxsRYXkjFjxlisPg7G\nWIP/l759+1r8fpais36WLA7roixYsICJRCK2devWZpV/8cUXmUgkMvuTmZlpUPbzzz9n48ePZ25u\nbszBwYH5+vqy1157jd2+fduk3v79+zORSMSqqqrY8uXLWf/+/ZmDgwPz8fFh69atY/fv32/1M4aF\nhTGRSMTy8vJadb2rqyuztbVt9f3N2WL8XjXGoUOH2Pjx45mXlxdzcHBgXl5ebPTo0eytt97iyzT0\nPxkwYIDJvfVJT09nIpGIRUVFscLCQhYREcG6d+/OunXrxiZOnMjOnTvHGGOstLSUvfTSS8zLy4s5\nOjqy4OBglp6ebmLrr7/+ytasWcNCQkJYr169mL29Pevduzd75pln2Pnz5w3Krlq1qkG79+zZY1A2\nLS2NTZs2jbm7uzN7e3vWv39/tmDBAlZSUtLs97Ez0GV7JD169AAA5OXlNav8k08+ifLycnz55ZdQ\nKBQG3fD+/fvzrxcsWIAdO3agb9++iIiIgJubG7KzsxEfH4+jR4/ixIkTkEgkJvU/9dRTOHv2LJ56\n6inY2dnh4MGDWLlyJf773//i4MGDwh62FWRlZaG8vByBgYEm56KiovDRRx9h1apVWLVqVYvqZc38\nBn7vvffw97//HZ6enpg6dSo8PDxQVlaG8+fP491338Xy5csBAKtWrcKePXtw+fJlxMbGwtXVFQD4\n3xwNDSEuXbqEkSNHQi6XY968eTh//jyOHTuGsWPH4sSJE/jLX/4Cd3d3PPPMM/j1119x4MABTJo0\nCfn5+QY9nePHjyM+Ph7jxo3DsGHDIJFIkJ+fj5SUFBw6dAhZWVn8ezl27FhcvnwZH374IQIDAzFz\n5ky+nqCgIP51fHw8li9fjh49emDKlCnw9PSEWq3G9u3bcejQIZw8eRJ9+vRp1vvZ4XS0krUVKpWK\n2dvbM7FYzJ5//nn2+eefs4sXL7K6uroGr9m9ezcTiURszZo1Zs9//PHHTCQSsYiICFZZWWlwbv36\n9UwkErGlS5caHOd6JEOGDGFarZY/XlFRwUaMGMFEIhHbu3dvq56xJT2SAwcOsFWrVrFly5axmTNn\nMnt7e9azZ0928uRJk7JRUVFMLBY3+D40ZktUVBRbtWqV2Z9Lly7x5YcNG8YcHR1ZaWmpSV03btww\nW/fly5cbvLdYLDY4xvVIRCIRe/vttw3OzZ8/n4lEIta9e3cWGxtrcO6tt95iIpGIxcXFGRwvLS1l\nOp3O5N5qtZpJJBL2l7/8xeB4RkYGE4lEbM6cOWZtzszMZCKRiIWEhLDy8nKDc/qfM2uhywoJY4zt\n37+feXl5GXQtpVIpmzJlCvvss89MRKUpIRk2bBizs7Njt27dMjlXU1PD3N3dmaenp8FxTkj+/e9/\nm1zzzTffMJFIxCZMmNCq52uJkMyePdvgffDz82NnzpwxW/bq1assLy+PlZWVtdiWhn7EYrHBsGfY\nsGHMxcWF3bx5s9l1t0ZIvL29TcorlUomEolYt27d2L179wzOXblyhYlEIjZu3LjmPDZjjLFp06Yx\nR0dHVlNTY3L/hoQkPDyciUQi9vPPP5s9HxQUxGxtbc2KV2ekyw5tACAyMhJPPvkk0tPTceLECahU\nKmRlZeHo0aM4evQo9uzZg0OHDsHOzq7Juu7duweVSoUePXrgnXfeMVvG3t4eV69exc2bN9G9e3f+\nuEgkQlhYmEn50NBQAEBOTk4rn7D57N27F3v37oVOp8O5c+ewbt06jB49GsnJyZg7d65BWU9PT3h6\nerbqPhkZGc1ytj7//PNYunQp/P39MWvWLISGhiIkJAReXl6tum9DyOVyk2Pcs/n6+sLJycngHHd/\njUZjct2RI0ewY8cOnDlzBjdu3EBNTQ1/TiQSoaysDL169WqWXSdOnICtrS1SUlLw+eefm5yvqqpC\nbW0t8vLyMGzYsGbV2ZF0aSEBAFtbW0yYMAETJkwAUD8tnJKSgrlz5+Lrr7/G9u3bERMT02Q9t27d\nAgDcuHEDa9eubbCcSCTC3bt3DYQEgNkPmKOjIx566KEmp3AtiUQiwciRI/Hll18iODgYCxYswOOP\nP45+/fq1mw0AEBcXBw8PD2zfvh3btm1DQkICAGDkyJHYuHGjWeFtDVKp1OSYra1tk+eqq6sNjm/d\nuhVxcXHo3r07JkyYgH79+sHZ2RkikQgHDx6EWq1GVVVVs+26ceMGamtrG42t4T5L1kCXFxJjxGIx\nIiMj8fPPP2P9+vVIS0trlpBwH7pHHnmkVT2IkpISk2nKyspK3L59G+7u7i2uTyi2trYYO3Ysfvrp\nJ5w8ebLdhQQAnnvuOTz33HPQ6XTIzs7GoUOHsHPnTkyaNAlqtRq+vr7tbpM5ampqsHr1anh5eeHs\n2bMmXwonTpxocZ1SqRQ1NTXQarWWMrND6dKRrY3BzawwvVkGGxsbAEBtba3Z8kOHDsX//vc/3Lhx\no0X3YowhMzPT5DgXA6LvyW9Pfv31VwB/fAt3FBKJBBMmTEBiYiKWLl2KyspKHDt2jD/f2P+lPSgr\nK0N5eTlCQkJMRESn0+Hs2bMms0ZN2RwSEoLbt2/jp59+ahuj25kuKyR79+7Fd999Z3Y6sqSkBDt3\n7gQAg/E81zO4fPmy2TqXLl2K6upqzJkzhx/q6HPnzh38+OOPZq9dt26dwRCmoqICK1asAADBYdTm\nnhEAbt68iYsXL5o9d/jwYRw8eBAuLi4mEaclJSXIzc1tsWC2hLS0NLPHr169CgBwcXHhj3FT+Q39\nX9oaDw8PODs748yZMwZDjerqaixevNjs+9SUzUuWLAEAzJ8/nxd0fSorKxsMmOyMdNmhzY8//oit\nW7fC09MTjz32GB+2XVxcjCNHjqCyshIhISGIjo7mrwkJCYGLiws+++wz2NnZoV+/fhCJRHjhhRfQ\nr18/vPjiizh79iwSExPh7e2NJ554Av3794dWq8WlS5egVCrxxBNPIDU11cSeP/3pTwgICMBTTz0F\nGxsbfPHFFyguLsbMmTMxe/bsZj9XVFQU/zo3NxcAsHz5cjz00EMA6uNhZsyYAQC4cuUKHn30UQwf\nPhyDBw9Gnz59oNVqkZOTg1OnTsHOzg47d+408ee89tprrY4j2b17d4Mi8ac//YlfkhAeHs77a/r3\n7w+RSIQff/wRWVlZ8PHxwdNPP81fN3HiRHz++ed46aWXEB4ejm7dusHNzQ0LFy7kyzQkppZALBYj\nJiYGmzZtwsMPP4zp06fj/v37SE9Ph1arxdixY5Genm5wjZ+fH/r27QulUonnn38evr6+sLGxwYwZ\nM/Dwww9DoVDgX//6F1599VX4+vpi8uTJGDhwICoqKnDlyhUcP34cgwYNwtmzZ9vsuSxKY1M6mZmZ\nbNq0aaxPnz5mo/IYq4/i6927N3NycmIKhcIkyq+j0Gg0bPv27Sw8PJz5+fkxqVTK7OzsmKenJ5sw\nYQJ77733DKbrOL799lv22GOPsW7dupmdtmSMsa+++orNmDGDeXp6Mnt7e9arVy/26KOPsmXLlrGc\nnByDsv3792disZiPbB0wYABzcHBg3t7ebN26day6urpFz8XZxP3mfri/9aeub926xVasWMFCQ0OZ\nl5cXs7e3ZxKJhPn7+7OXX36ZXbhwwew9WhNHolAoeDsa+nnyySf58jt27GDh4eHM29ububi4MFdX\nVyaXy9maNWtMpoTr6urYypUrmY+PD7O3t2cikYgNHDjQ5N76NDb9WlxczEQiERs7dqzZZzGun7H6\n6f23336b+fv7MycnJ+bl5cVeeOEFduXKFf79Mp6ePnv2LJswYQJzdXXl/08ffvihQZmTJ0+yZ555\nhslkMmZvb8/c3d3ZI488wqKjo5lSqWzkHe9ciBhrWMqPHTuGEydOICgoCC+88AK2b9+OF154gT8f\nHx+PDRs24MMPP8TgwYOxdu1aZGVlIS8vz2x054PIgAED8Msvv3TY+J4g2oXmKo5EIjFQ07q6Oubp\n6WmwLqKiooJ169aNvfvuu5YUO6uG65EQRFem1c7W4uJiXLt2DRMnTuSPOTo6YsyYMLine truncated
|
||||
"text/plain": [
|
||||
"<matplotlib.figure.Figure at 0x8a1a908>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from pf_internal import plot_pf\n",
|
||||
"\n",
|
||||
"seed(1234)\n",
|
||||
"N = 3000\n",
|
||||
"pf = ParticleFilter(N, 20, 20)\n",
|
||||
"xs = np.linspace (1, 10, 20)\n",
|
||||
"ys = np.linspace (1, 10, 20)\n",
|
||||
"zxs = xs + randn(20)\n",
|
||||
"zys = xs + randn(20)\n",
|
||||
"\n",
|
||||
"def animatepf(i):\n",
|
||||
" if i == 0:\n",
|
||||
" plot_pf(pf, 10, 10, weights=False)\n",
|
||||
" \n",
|
||||
" idx = int((i-1) / 3)\n",
|
||||
" x, y = xs[idx], ys[idx]\n",
|
||||
" z = [x + randn()*0.2, y + randn()*0.2]\n",
|
||||
"\n",
|
||||
" step = (i % 3) + 1\n",
|
||||
" if step == 2:\n",
|
||||
" pf.predict((0.5, 0.5), (0.2, 0.2))\n",
|
||||
" pf.weight(z=z, var=.6)\n",
|
||||
" plot_pf(pf, 10, 10, weights=False)\n",
|
||||
" plt.title('Step {}: Predict'.format(idx+1))\n",
|
||||
" elif step == 3:\n",
|
||||
" pf.resample()\n",
|
||||
" plot_pf(pf, 10, 10, weights=False)\n",
|
||||
" plt.title('Step {}: Resample'.format(idx+1))\n",
|
||||
"\n",
|
||||
" else:\n",
|
||||
" mu, var = pf.estimate()\n",
|
||||
" plot_pf(pf, 10, 10, weights=False)\n",
|
||||
" plt.scatter(mu[0], mu[1], color='g', s=100, label='PF')\n",
|
||||
" plt.scatter(x, y, marker='x', color='r', s=180, lw=3, label='Robot')\n",
|
||||
" plt.title('Step {}: Estimate'.format(idx+1))\n",
|
||||
" #plt.scatter(mu[0], mu[1], color='g', s=100, label=\"PF\")\n",
|
||||
" #plt.scatter([x+1], [x+1], marker='x', color='r', s=180, label=\"True\", lw=3)\n",
|
||||
" plt.legend(scatterpoints=1, loc=2)\n",
|
||||
" plt.tight_layout()\n",
|
||||
"\n",
|
||||
"from gif_animate import animate\n",
|
||||
"animate('particle_filter_anim.gif', animatepf, \n",
|
||||
" frames=40, interval=800, figsize=(4, 4))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<img src='particle_filter_anim.gif'>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.5.0"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 316 KiB |
Reference in new issue
Block a user