Re: Kalman filtering resources
- To: mathgroup at smc.vnet.net
- Subject: [mg15527] Re: Kalman filtering resources
- From: flip phillips <flip at mathematica-journal.com>
- Date: Mon, 18 Jan 1999 23:47:26 -0500
- Organization: the mathematica journal
- References: <774l8v$47j@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
oof! looks like igor beat me to the punch re: the time series package- also- i was recently able to implement a simple 1-d kalman filter for a little GPS project i was doing, I've attached the NB here for your use. It is really a trivial example, 1-D, no state change, no conrol input... still, it's great fun at parties. enjoy- fp CLERC Nicolas wrote: > Hi Mathematica Folks, > > I'm working on a financial problem that would require some Kalman > Filtering (a conditional CAPM with latent factors). > > I found nothing related to KF on the MathSource site. > > Does anyone know any relevant resources (URL or ftp sites) where I would > find some Mathematica code and/or advices for KF ? > > Thank you very much in advance for your support, > > Best, > > Nicolas -- flip phillips technical editor, the mathematica journal flip at mathematica-journal.com http://www.mathematica-journal.com filename="kalman.nb" (*********************************************************************** Mathematica-Compatible Notebook This notebook can be used on any computer system with Mathematica 3.0, MathReader 3.0, or any compatible application. The data for the notebook starts with the line of stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info at wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. ***********************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 20868, 705]*) (*NotebookOutlinePosition[ 21912, 741]*) (* CellTagsIndexPosition[ 21837, 735]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Simple Kalman Filter", "Title"], Cell[TextData[{ "Flip Phillips\n", ButtonBox["flip at skidmore.edu", ButtonData:>{ URL[ "mailto:flip at skidmore.edu"], None}, ButtonStyle->"Hyperlink"], "\n", ButtonBox["flip at mathematica-journal.com", ButtonData:>{ URL[ "mailto:flip at mathematica-journal.com"], None}, ButtonStyle->"Hyperlink"] }], "Subsubtitle"], Cell[TextData[{ "This is a pretty simple example of implementing a Kalman filter. There is \ no real tutorial as to the hows and whys of the Kalman filter, for a good \ tutorial see the ", ButtonBox["UNC Tech Report", ButtonData:>"ref:welch", ButtonStyle->"Hyperlink"], " cited at the end of the paper. In fact, this example is pretty much a \ straight forward duplication of one of their examples." }], "Text"], Cell[CellGroupData[{ Cell["Preliminaries", "Section"], Cell["\<\ So we can generate some nice noise\[Ellipsis]\ \>", "Text"], Cell[BoxData[ \(<< Statistics`NormalDistribution`\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Our friendly function", "Section"], Cell[TextData[{ "Let's say that we have a noisy system whose 'true' value, ", Cell[BoxData[ \(TraditionalForm\`z\)]], ", is some random constant:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(z = \(-0.6\)\)], "Input"], Cell[BoxData[ FormBox[ RowBox[{"-", StyleBox["0.6`", StyleBoxAutoDelete->True, PrintPrecision->1]}], TraditionalForm]], "Output"] }, Open ]], Cell["\<\ Now, we make a normally distributed noise function with a \[Mu] of \ 0 and \[Sigma] of some level:\ \>", "Text"], Cell[BoxData[ \(sd = 0.1; \nnDist = NormalDistribution[0, sd]; \n noise[] := Random[nDist]\)], "Input"], Cell["Finally, define the initial noisy function itself:", "Text"], Cell[BoxData[ \(zf[] := z + noise[]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["A sample of our system", "Section"], Cell[TextData[{ Cell[BoxData[ \(TraditionalForm\`d\)]], " contains ", StyleBox["iter", FontSlant->"Italic"], " samples of our function." }], "Text"], Cell[BoxData[ \(iter = 100; \nd = Table[zf[], {iter}]; \)], "Input"], Cell["\<\ And here we can take a look at it, along with the true value, \ plotted as a red line:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(g0 = ListPlot[d, Epilog \[Rule] {Hue[1], Line[{{0, z}, {iter, z}}]}]; \)\)], "Input"], Cell[GraphicsData["PostScript", "\<\ %! %%Creator: Mathematica %%AspectRatio: .61803 MathPictureStart /Mabs { Mgmatrix idtransform Mtmatrix dtransform } bind def /Mabsadd { Mabs 3 -1 roll add 3 1 roll add exch } bind def %% Graphics %%IncludeResource: font Courier %%IncludeFont: Courier /Courier findfont 10 scalefont setfont % Scaling calculations 0.0238095 0.00952381 1.19672 1.51654 [ [.21429 .12264 -6 -9 ] [.21429 .12264 6 0 ] [.40476 .12264 -6 -9 ] [.40476 .12264 6 0 ] [.59524 .12264 -6 -9 ] [.59524 .12264 6 0 ] [.78571 .12264 -6 -9 ] [.78571 .12264 6 0 ] [.97619 .12264 -9 -9 ] [.97619 .12264 9 0 ] [.01131 .28679 -24 -4.5 ] [.01131 .28679 0 4.5 ] [.01131 .43845 -24 -4.5 ] [.01131 .43845 0 4.5 ] [.01131 .5901 -24 -4.5 ] [.01131 .5901 0 4.5 ] [ 0 0 0 0 ] [ 1 .61803 0 0 ] ] MathScale % Start of Graphics 1 setlinecap 1 setlinejoin newpath 0 g ..25 Mabswid [ ] 0 setdash ..21429 .13514 m ..21429 .14139 L s [(20)] .21429 .12264 0 1 Mshowa ..40476 .13514 m ..40476 .14139 L s [(40)] .40476 .12264 0 1 Mshowa ..59524 .13514 m ..59524 .14139 L s [(60)] .59524 .12264 0 1 Mshowa ..78571 .13514 m ..78571 .14139 L s [(80)] .78571 .12264 0 1 Mshowa ..97619 .13514 m ..97619 .14139 L s [(100)] .97619 .12264 0 1 Mshowa ..125 Mabswid ..07143 .13514 m ..07143 .13889 L s ..11905 .13514 m ..11905 .13889 L s ..16667 .13514 m ..16667 .13889 L s ..2619 .13514 m ..2619 .13889 L s ..30952 .13514 m ..30952 .13889 L s ..35714 .13514 m ..35714 .13889 L s ..45238 .13514 m ..45238 .13889 L s ..5 .13514 m ..5 .13889 L s ..54762 .13514 m ..54762 .13889 L s ..64286 .13514 m ..64286 .13889 L s ..69048 .13514 m ..69048 .13889 L s ..7381 .13514 m ..7381 .13889 L s ..83333 .13514 m ..83333 .13889 L s ..88095 .13514 m ..88095 .13889 L s ..92857 .13514 m ..92857 .13889 L s ..25 Mabswid 0 .13514 m 1 .13514 L s ..02381 .28679 m ..03006 .28679 L s [(-0.6)] .01131 .28679 1 0 Mshowa ..02381 .43845 m ..03006 .43845 L s [(-0.5)] .01131 .43845 1 0 Mshowa ..02381 .5901 m ..03006 .5901 L s [(-0.4)] .01131 .5901 1 0 Mshowa ..125 Mabswid ..02381 .16547 m ..02756 .16547 L s ..02381 .1958 m ..02756 .1958 L s ..02381 .22613 m ..02756 .22613 L s ..02381 .25646 m ..02756 .25646 L s ..02381 .31712 m ..02756 .31712 L s ..02381 .34745 m ..02756 .34745 L s ..02381 .37778 m ..02756 .37778 L s ..02381 .40812 m ..02756 .40812 L s ..02381 .46878 m ..02756 .46878 L s ..02381 .49911 m ..02756 .49911 L s ..02381 .52944 m ..02756 .52944 L s ..02381 .55977 m ..02756 .55977 L s ..02381 .10481 m ..02756 .10481 L s ..02381 .07448 m ..02756 .07448 L s ..02381 .04415 m ..02756 .04415 L s ..02381 .01382 m ..02756 .01382 L s ..25 Mabswid ..02381 0 m ..02381 .61803 L s 0 0 m 1 0 L 1 .61803 L 0 .61803 L closepath clip newpath ..008 w ..03333 .03156 Mdot ..04286 .12178 Mdot ..05238 .0222 Mdot ..0619 .38764 Mdot ..07143 .29149 Mdot ..08095 .39042 Mdot ..09048 .17822 Mdot ..1 .4717 Mdot ..10952 .33344 Mdot ..11905 .3865 Mdot ..12857 .23656 Mdot ..1381 .36234 Mdot ..14762 .22694 Mdot ..15714 .08944 Mdot ..16667 .30609 Mdot ..17619 .10661 Mdot ..18571 .37695 Mdot ..19524 .23267 Mdot ..20476 .10122 Mdot ..21429 .16287 Mdot ..22381 .16472 Mdot ..23333 .31786 Mdot ..24286 .46501 Mdot ..25238 .24552 Mdot ..2619 .03828 Mdot ..27143 .42117 Mdot ..28095 .34721 Mdot ..29048 .01472 Mdot ..3 .24829 Mdot ..30952 .38052 Mdot ..31905 .12164 Mdot ..32857 .4457 Mdot ..3381 .24053 Mdot ..34762 .14739 Mdot ..35714 .44207 Mdot ..36667 .25755 Mdot ..37619 .29003 Mdot ..38571 .25914 Mdot ..39524 .19242 Mdot ..40476 .33753 Mdot ..41429 .25722 Mdot ..42381 .14494 Mdot ..43333 .11622 Mdot ..44286 .39171 Mdot ..45238 .16474 Mdot ..4619 .32785 Mdot ..47143 .32962 Mdot ..48095 .53362 Mdot ..49048 .25073 Mdot ..5 .06906 Mdot ..50952 .39791 Mdot ..51905 .30824 Mdot ..52857 .15923 Mdot ..5381 .47791 Mdot ..54762 .52189 Mdot ..55714 .12113 Mdot ..56667 .13565 Mdot ..57619 .46677 Mdot ..58571 .05859 Mdot ..59524 .37742 Mdot ..60476 .4017 Mdot ..61429 .22673 Mdot ..62381 .1401 Mdot ..63333 .41167 Mdot ..64286 .35077 Mdot ..65238 .23587 Mdot ..6619 .44248 Mdot ..67143 .02562 Mdot ..68095 .36632 Mdot ..69048 .26378 Mdot ..7 .24696 Mdot ..70952 .38528 Mdot ..71905 .07702 Mdot ..72857 .21008 Mdot ..7381 .40891 Mdot ..74762 .16632 Mdot ..75714 .1713 Mdot ..76667 .46003 Mdot ..77619 .03545 Mdot ..78571 .03593 Mdot ..79524 .21139 Mdot ..80476 .35172 Mdot ..81429 .1298 Mdot ..82381 .39469 Mdot ..83333 .60332 Mdot ..84286 .21631 Mdot ..85238 .31617 Mdot ..8619 .53767 Mdot ..87143 .08532 Mdot ..88095 .34952 Mdot ..89048 .57646 Mdot ..9 .17395 Mdot ..90952 .27613 Mdot ..91905 .05957 Mdot ..92857 .42349 Mdot ..9381 .23191 Mdot ..94762 .47184 Mdot ..95714 .23528 Mdot ..96667 .38204 Mdot ..97619 .19485 Mdot 1 0 0 r ..5 Mabswid ..02381 .28679 m ..97619 .28679 L s % End of Graphics MathPictureEnd \ \>"], "Graphics", ImageSize->{288, 177.938}, ImageMargins->{{43, 0}, {0, 0}}, ImageRegion->{{0, 1}, {0, 1}}, ImageCache->GraphicsData["Bitmap", "\<\ CF5dJ6E]HGAYHf4PAg9QL6QYHg<PAVmbKF5d0`40004P0000/A000`40O003h00Oog<`8G<`003oLc0Q Lc000?mc<25c<000og<`8G<`003oLc0QLc000?mc<25c<000og<`8G<`000KLc000`00Lc1c<03oLc03 Lc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<`Lc00og<`0g<`000KLc020013Lc02002nLc00 01]c<003001c<7<`049c<0800;ic<0006g<`00<007<`Lc0017<`0P00o7<`000KLc000`00Lc1c<004 Lc02002QLc02001ILc0001]c<004001c<000002VLc02001ILc0001]c<004001c<000003oLc02Lc00 01]c<003001c<7<`03]c<08008Ec<0800003Lc00000003]c<0006g<`00<007<`Lc00>g<`0P00QG<` 0P0000=c<0000000>g<`000KLc02003oLc04Lc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<` Lc00og<`0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`091c<08005Ic<08001Qc<000 6g<`00<007<`Lc00T7<`0P00EW<`0P0067<`000KLc000`00Lc1c<01jLc020026Lc0001]c<003001c <7<`07Yc<08008Ic<0006g<`0P00]7<`0P00CG<`000KLc000`00Lc1c<02cLc02001=Lc0001]c<003 001c<7<`0?mc<0=c<0006g<`00<007<`Lc00fg<`0P009G<`000KLc000`00Lc1c<00OLc02002jLc02 000ULc0001]c<003001c<7<`01mc<0800>5c<0006g<`00<007<`Lc00:W<`10000g<`0P00:W<`0`00 0g<`0P00:W<`0P0017<`0P00:W<`0P0017<`0P009g<`0`0017<`0P0017<`0P000g<`000KLc000`00 Lc1c<00[Lc030002Lc001000Lc1c<000:W<`00D007<`Lc1c<000009c<003001c<7<`02Ic<004001c <7<`0002Lc001000Lc1c<000:7<`00 at 007<`Lc00009c<004001c<7<`000WLc000`00Lc1c<002Lc00 1000Lc1c<0000W<`00 at 007<`Lc00009c<0006g<`0P009G<`0P001W<`0P000W<`00 at 007<`Lc0002Qc <0 at 0009c<004001c<7<`000XLc001000Lc1c<0000W<`00 at 007<`Lc0002Qc<004001c<7<`0002Lc00 1000Lc1c<0009g<`00<007<`Lc000W<`00 at 007<`Lc00009c<004001c<7<`0002Lc0001]c<003001c <7<`02Ac<08000Mc<004001c<7<`0002Lc000`00Lc1c<00VLc000`00Lc000003Lc001000Lc1c<000 :7<`0`000g<`00 at 007<`Lc0002Uc<08000=c<004001c<7<`000WLc000`00Lc1c<002Lc001000Lc1c <0000W<`00 at 007<`Lc00009c<0006g<`00<007<`Lc00:W<`00 at 007<`Lc00009c<004001c<7<`000X Lc000`00Lc000003Lc001000Lc1c<000:7<`00<007<`Lc000g<`00 at 007<`Lc0002Qc<004001c<7<` 0002Lc001000Lc1c<0009g<`00<007<`Lc000W<`00 at 007<`Lc00009c<004001c<7<`0002Lc0001]c <003001c<7<`02]c<08000Ac<08002Yc<08000Ac<08000=c<08002Ec<0<000=c<08002Yc<08000Ac <08002Mc<08000Ec<08000Ac<08000=c<0006g<`00 at 007<`Lc1c<08004Mc<08001ac<08001mc<080 07Mc<0006g<`00 at 007<`Lc1c<08004Mc<08003ec<08007Mc<0006g<`00<007<`Lc00c7<`0P00=7<` 000KLc000`00Lc1c<03<Lc02000dLc0001Ec<?l000T000=c<0006g<`00<007<`Lc002G<`00<007<` Lc002W<`00<007<`Lc002G<`00<007<`Lc002W<`00<007<`Lc002G<`00<007<`Lc002W<`00<007<` Lc002W<`00<007<`Lc002G<`00<007<`Lc002W<`00<007<`Lc002G<`00<007<`Lc002W<`00 at 007<` Lc1c<08000Ic<003001c<7<`00Ac<08000Ac<003001c<7<`00Uc<003001c<7<`00Yc<003001c<7<` 00Uc<003001c<7<`00Yc<003001c<7<`00Uc<003001c<7<`00Yc<003001c<7<`00Yc<003001c<7<` 00Mc<0006g<`00<007<`Lc00;g<`00<007<`Lc00;g<`00<007<`Lc000W<`0P00:g<`00<007<`Lc00 17<`0P00:G<`00<007<`Lc00<7<`00<007<`Lc001g<`000KLc000`00Lc1c<01BLc02000BLc02002J Lc0001]c<003001c<7<`059c<0800:ic<0006g<`00<007<`Lc00og<`0g<`000KLc000`00Lc1c<021 Lc02001oLc0001]c<003001c<7<`02ic<080055c<08007mc<0006g<`0P00;g<`0P0000=c<0000000 >W<`0P00C7<`0P00AG<`000KLc000`00Lc1c<00aLc02000jLc02001<Lc040013Lc0001]c<003001c <7<`0;ec<08002Ac<08001ec<0006g<`00<007<`Lc003W<`0P00dg<`0P007G<`000KLc000`00Lc1c <00>Lc02003bLc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<`Lc00og<`0g<`000KLc000`00 Lc1c<01NLc02002GLc020009Lc0001]c<08005mc<08009Mc<08000Uc<0006g<`00<007<`Lc00og<` 0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0;Ic<08004Yc<0006g<`00<007<`Lc00 ]W<`0P003g<`0P00>G<`000KLc000`00Lc1c<037Lc02000;Lc02000/Lc0001]c<003001c<7<`0=Ac <08002ac<0006g<`00<007<`Lc00og<`0g<`000KLc02000NLc02001iLc02001XLc0001]c<003001c <7<`01ec<08000Yc<08006ec<08005=c<08001=c<0006g<`00<007<`Lc0067<`0P003g<`0P00Mg<` 0P00BG<`0P000g<`0P003W<`000KLc000`00Lc1c<00HLc020028Lc02001>Lc02000>Lc0001]c<003 001c<7<`04mc<0800;5c<0006g<`00<007<`Lc00>7<`0P005G<`0P00GG<`0P00DW<`000KLc000`00 Lc1c<00hLc02000;Lc02000`Lc02000eLc02001BLc0001]c<003001c<7<`04Ec<080031c<08008Uc <0006g<`0P00F7<`0P000g<`0P001G<`0P00WG<`000KLc000`00Lc1c<01GLc020003Lc020005Lc02 002MLc0001]c<003001c<7<`0:ac<08005Ac<0006g<`00<007<`Lc00[7<`0P00E7<`000KLc000`00 Lc1c<03oLc03Lc0001]c<003001c<7<`0>Ec<08001]c<00027<`0P0017<`0P0017<`0P001G<`00<0 07<`Lc00iG<`0P006g<`0007Lc001000Lc1c<0000g<`0P000g<`00 at 007<`Lc0000Ac<003001c<7<` 0?mc<0=c<00000=c<00000000`0000Ec<000Lc1c<00000Qc<004001c<7<`0004Lc3lO009Lc0000Mc <004001c<7<`0008Lc030005Lc000`00Lc1c<009Lc02001>Lc02002WLc0000Mc<004001c<7<`0008 Lc000`00Lc1c<005Lc000`00Lc1c<009Lc02001>Lc02002WLc0000Qc<08000Yc<0<000Ac<003001c <7<`0?mc<0=c<0006g<`00<007<`Lc00og<`0g<`000KLc000`00Lc1c<00RLc02001KLc020021Lc00 01]c<003001c<7<`029c<08005]c<080085c<0006g<`00<007<`Lc00eW<`0P00:W<`000KLc02000d Lc02002QLc02000ZLc0001]c<003001c<7<`03=c<0800<ec<0006g<`00<007<`Lc00L7<`0P00T7<` 000KLc000`00Lc1c<01`Lc04002>Lc0001]c<003001c<7<`01=c<08005ec<08008ic<0006g<`00<0 07<`Lc004g<`0P00C7<`0P00Wg<`000KLc000`00Lc1c<01QLc02002OLc0001]c<003001c<7<`0?mc <0=c<0006g<`0P00 at G<`0P00W7<`0P008W<`000KLc000`00Lc1c<010Lc02001MLc02000YLc02000B Lc02000RLc0001]c<003001c<7<`09mc<08002Uc<08003Ic<0006g<`00<007<`Lc00og<`0g<`000K Lc000`00Lc1c<00JLc02003VLc0001]c<003001c<7<`01Yc<08008ec<08005Mc<0006g<`00<007<` Lc00ZG<`0P00Eg<`000KLc000`00Lc1c<00WLc02003ILc0001]c<08002Qc<08001mc<08004Uc<080 06ec<0006g<`00<007<`Lc00B7<`0P00BG<`0P00Gg<`0P0037<`000KLc000`00Lc1c<006Lc02000= Lc02002JLc020011Lc02000<Lc0001]c<003001c<7<`00Ic<08000=c<08000Qc<08005Ac<08004Ac <08004mc<0006g<`00<007<`Lc002g<`0P00GW<`0P00HW<`0P00<G<`000KLc000`00Lc1c<01lLc02 001ALc02000aLc0001]c<003001c<7<`07ac<08001Mc<08006]c<0006g<`00<007<`Lc00UG<`0P00 Jg<`000KLc02002iLc020018Lc0001]c<003001c<7<`09ec<08001Uc<08004Qc<0006g<`00<007<` Lc00WG<`0P00Hg<`000KLc000`00Lc1c<00mLc020033Lc0001]c<003001c<7<`03ec<0800:]c<080 01Ic<0006g<`00<007<`Lc00jW<`0P005W<`0008Lc020004Lc020003Lc030005Lc000`00Lc1c<03o Lc03Lc0000Mc<004001c<7<`0003Lc020006Lc000`00Lc1c<002Lc000`00Lc1c<03oLc03Lc000003 Lc00000000<00005Lc0007<`Lc00000;Lc000`00Lc1c<002Lc02003oLc04Lc0000Mc<004001c<7<` 0008Lc030005Lc000`00Lc1c<01DLc02001>Lc02001LLc0000Mc<004001c<7<`0008Lc000`00Lc1c <005Lc000`00Lc1c<01=Lc020005Lc02001>Lc02001LLc0000Qc<08000Uc<0 at 000Ac<003001c<7<` 04ec<0800;=c<0006g<`00<007<`Lc00og<`0g<`000KLc000`00Lc1c<030Lc020010Lc0001]c<003 001c<7<`0<1c<080041c<0006g<`00<007<`Lc00=W<`0P00EW<`0P00LW<`000KLc02000ALc02000T Lc02001FLc02001OLc02000ALc0001]c<003001c<7<`011c<0800=ec<080015c<0006g<`00<007<` Lc00Q7<`0P00O7<`000KLc000`00Lc1c<024Lc02001lLc0001]c<003001c<7<`0?mc<0=c<0006g<` 00<007<`Lc00og<`0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0?mc<0=c<0006g<` 0P00og<`17<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<` Lc00og<`0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<` Lc00QW<`0P00NW<`000KLc000`00Lc1c<026Lc02001jLc0001]c<0800?mc<0Ac<0006g<`00<007<` Lc00MG<`0P00Rg<`000KLc000`00Lc1c<01eLc02001RLc02000WLc0001]c<003001c<7<`0=Uc<080 02Mc<0006g<`00<007<`Lc00og<`0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0?mc <0=c<0006g<`00<007<`Lc00og<`0g<`000KLc02003oLc04Lc0001]c<003001c<7<`0?mc<0=c<000 6g<`00<007<`Lc00og<`0g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0>1c<080021c <0006g<`00<007<`Lc00h7<`0P0087<`0008Lc020004Lc020004Lc030004Lc000`00Lc1c<03oLc03 Lc0000Mc<004001c<7<`0003Lc020005Lc000`00Lc1c<003Lc000`00Lc1c<03oLc03Lc000003Lc00 000000<00005Lc0007<`Lc000008Lc040004Lc02003oLc04Lc0000Mc<004001c<7<`0008Lc000`00 Lc000005Lc000`00Lc1c<03oLc03Lc0000Mc<004001c<7<`0008Lc000`00Lc000005Lc000`00Lc1c <03oLc03Lc0000Qc<08000Yc<08000Ec<003001c<7<`0=5c<08002mc<0006g<`00<007<`Lc00dG<` 0P00;g<`000KLc000`00Lc1c<03oLc03Lc0001]c<003001c<7<`0?mc<0=c<0006g<`00<007<`Lc00 og<`0g<`003oLc0QLc000?mc<25c<000og<`8G<`003oLc0QLc000?mc<25c<000og<`8G<`003oLc0Q Lc000001\ \>"], ImageRangeCache->{{{0, 287}, {176.938, 0}} -> {-10.8859, -0.806954, 0.398872, 0.0025049}}] }, Open ]], Cell["Now, can we estimate the constant from our noisy data?", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["A filter", "Section"], Cell["\<\ In this simple Kalman filter the state does not change and there is \ no control input, thus we end up with a simple model, the time update \ equations are:\ \>", "Text"], Cell[TextData[{ Cell[BoxData[ \(TraditionalForm\`\(\(x\&^\)\_\(k + 1\)\%-\) = \(x\&^\)\)]], ",\n", Cell[BoxData[ \(TraditionalForm\`\(P\_\(k + 1\)\%-\) = P\_k + Q\)]], "," }], "Text"], Cell["and the measurement update equations are:", "Text"], Cell[TextData[Cell[BoxData[ \(TraditionalForm \`K\_k = \(P\_k\%-\)\/\(\(P\_k\%-\) + R\), \[IndentingNewLine]\(x\&^\)\_k = \(\(x\&^\)\_k\%-\) + K(z\_k - \(\(x\&^\)\_k\%-\)), \[IndentingNewLine]P\_k = \((1 - K\_k)\) \(P\_k\%-\)\)]]], "Text"], Cell[TextData[{ "Finally, here is some ", StyleBox["Mathematica", FontSlant->"Italic"], " code to filter the data, x0 is an seed value for the filter, Pk0 is the \ initial scaling paramter used for convergence, R the measurement error \ variance estimate and Q the process error estimates. It returns a list of \ pairs consisting of the current ", Cell[BoxData[ \(TraditionalForm\`P\_k\)]], " for the given time step and the respective filter estimate ", Cell[BoxData[ \(TraditionalForm\`\(x\&^\)\_k\)]], ". " }], "Text"], Cell[BoxData[ \(runfilter[data_, x0_, Pk0_, R_, Q_] := Module[{x = x0, Pk = Pk0, res}, \[IndentingNewLine]Kk[] := Pk\/\(Pk + R\); \[IndentingNewLine]Map[ \[IndentingNewLine]\((res = {Pk, x + Kk[] \((# - x)\)}; \[IndentingNewLine]x = res\[LeftDoubleBracket]2\[RightDoubleBracket]; \[IndentingNewLine]Pk = \((Pk + Q)\) \((1 - Kk[])\); \[IndentingNewLine]res)\)&, data]\[IndentingNewLine]]\)], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Demonstration", "Section"], Cell[TextData[{ "This strangely named function takes a packet of data, a measurement error \ variance and process error estimate, filters it, plots the resulting filter \ output, the set of ", Cell[BoxData[ \(TraditionalForm\`P\_k\)]], "s and finally a combination graph showing the initial data, true value and \ filter estimate. " }], "Text"], Cell[BoxData[ \(DoIt[data_, R_, Q_] := Module[{Pks, xs, g1, g2}, \[IndentingNewLine]\[IndentingNewLine]{Pks, xs} = Transpose[\[IndentingNewLine]runfilter[d, 0, 1.0, R, Q]]; \[IndentingNewLine]g1 = ListPlot[xs, PlotJoined \[Rule] True]; \[IndentingNewLine]ListPlot[Pks, PlotJoined \[Rule] True]; \[IndentingNewLine]Show[{g0, g1}]; \[IndentingNewLine]]\)], "Input"], Cell[CellGroupData[{ Cell["Some examples", "Subsection"], Cell[TextData[{ "These examples show the filter in operation with varying levels of ", Cell[BoxData[ \(TraditionalForm\`Q\)]], " and ", Cell[BoxData[ \(TraditionalForm\`R\)]], ". It is interesting to notice that the actual measurement error variance \ is ", Cell[BoxData[ \(TraditionalForm\`\[Sigma]\^2\)]], " or 0.01 in this example. We'd expect that to give us the best result:" }], "Text"], Cell[BoxData[ \(DoIt[d, 0.01, 10. \^\(-5\)]\)], "Input"], Cell[TextData[{ "and indeed it does a pretty good job. Luckily, of course, we ", StyleBox["knew", FontSlant->"Italic"], " the actual ", Cell[BoxData[ \(TraditionalForm\`\(\(\(\[Sigma]\^2\)\(.\)\)\(\ \)\)\)]], "What happens when we guess a little bit?" }], "Text"], Cell[BoxData[ \(DoIt[d, 0.0001, 10. \^\(-5\)]\)], "Input"], Cell["noisy fit...", "Text"], Cell[BoxData[ \(DoIt[d, 1, 10. \^\(-5\)]\)], "Input"], Cell[BoxData[ \(DoIt[d, .1, 10. \^\(-7\)]\)], "Input"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Conclusions", "Section"], Cell["\<\ Honestly, there are no useful conclusions to place here, this was \ just a self-imposed exercise to see if I could sucessfully implement this \ filter. \ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["References", "Section"], Cell["\<\ Welch,Gregory and Gary Bishop.1995. \"An Introduction to the Kalman \ Filter,\" University of North Carolina,Department of Computer Science,TR \ 95-041.\ \>", "Text", CellTags->"ref:welch"] }, Closed]] }, Open ]] }, FrontEndVersion->"Macintosh 3.0", ScreenRectangle->{{0, 1152}, {0, 850}}, WindowToolbars->"EditBar", WindowSize->{643, 705}, WindowMargins->{{1, Automatic}, {Automatic, 0}}, Magnification->1, StyleDefinitions -> "PastelColor.nb", MacintoshSystemPageSetup->"\<\ 00l0001804P000000_P2B?ofoo833P9F;085:0?l0 at 00005X0FP000003]P;J001 0 at 00I00100000@0200000BL?00400000000000000000030001000000000 at 0?o> okX?AP^^0201000000000000000001T1\>" ] (*********************************************************************** Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. ***********************************************************************) (*CellTagsOutline CellTagsIndex->{ "ref:welch"->{ Cell[20639, 696, 201, 5, 44, "Text", CellTags->"ref:welch"]} } *) (*CellTagsIndex CellTagsIndex->{ {"ref:welch", 21740, 728} } *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1731, 51, 37, 0, 68, "Title"], Cell[1771, 53, 342, 11, 58, "Subsubtitle"], Cell[2116, 66, 426, 9, 62, "Text"], Cell[CellGroupData[{ Cell[2567, 79, 32, 0, 50, "Section"], Cell[2602, 81, 69, 2, 26, "Text"], Cell[2674, 85, 66, 1, 35, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[2777, 91, 40, 0, 34, "Section"], Cell[2820, 93, 170, 5, 26, "Text"], Cell[CellGroupData[{ Cell[3015, 102, 45, 1, 35, "Input"], Cell[3063, 105, 168, 5, 34, "Output"] }, Open ]], Cell[3246, 113, 122, 3, 26, "Text"], Cell[3371, 118, 111, 2, 67, "Input"], Cell[3485, 122, 67, 0, 26, "Text"], Cell[3555, 124, 52, 1, 35, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[3644, 130, 41, 0, 34, "Section"], Cell[3688, 132, 166, 7, 28, "Text"], Cell[3857, 141, 72, 1, 51, "Input"], Cell[3932, 144, 110, 3, 26, "Text"], Cell[CellGroupData[{ Cell[4067, 151, 112, 2, 35, "Input"], Cell[4182, 155, 12292, 396, 194, 4589, 296, "GraphicsData", "PostScript", "Graphics"] }, Open ]], Cell[16489, 554, 70, 0, 26, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[16596, 559, 27, 0, 34, "Section"], Cell[16626, 561, 180, 4, 44, "Text"], Cell[16809, 567, 203, 7, 44, "Text"], Cell[17015, 576, 57, 0, 26, "Text"], Cell[17075, 578, 270, 5, 66, "Text"], Cell[17348, 585, 551, 14, 64, "Text"], Cell[17902, 601, 503, 10, 171, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[18442, 616, 32, 0, 34, "Section"], Cell[18477, 618, 356, 8, 62, "Text"], Cell[18836, 628, 421, 7, 147, "Input"], Cell[CellGroupData[{ Cell[19282, 639, 35, 0, 38, "Subsection"], Cell[19320, 641, 424, 12, 44, "Text"], Cell[19747, 655, 60, 1, 39, "Input"], Cell[19810, 658, 283, 8, 46, "Text"], Cell[20096, 668, 62, 1, 39, "Input"], Cell[20161, 671, 28, 0, 26, "Text"], Cell[20192, 673, 57, 1, 39, "Input"], Cell[20252, 676, 59, 1, 39, "Input"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[20360, 683, 30, 0, 50, "Section"], Cell[20393, 685, 177, 4, 44, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[20607, 694, 29, 0, 34, "Section"], Cell[20639, 696, 201, 5, 44, "Text", CellTags->"ref:welch"] }, Closed]] }, Open ]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)