MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

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.
***********************************************************************)



  • Prev by Date: Re: ContourPlot problem
  • Next by Date: Re: Using Mathematica over Telnet
  • Previous by thread: Re: Kalman filtering resources
  • Next by thread: a tricky limit