root/library/tests/test_kalman_QRexh.cpp @ 550

Revision 477, 1.5 kB (checked in by mido, 15 years ago)

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

  • Property svn:eol-style set to native
Line 
1
2#include <estim/kalman.h>
3#include <estim/particles.h>
4
5using namespace bdm;
6
7//These lines are needed for use of cout and endl
8using std::cout;
9using std::endl;
10
11int main() {
12        // Klaman filter
13        mat A, B, C, D, R, Q, P0;
14        vec mu0;
15        mat Mu0;// read from matlab
16        // input from Matlab
17        it_file fin ( "testKF.it" );
18
19        mat Dt, XQRt, eR, eQ;
20        int Ndat;
21
22        bool xxx = fin.seek ( "d" );
23
24        if ( !xxx ) {
25                it_error ( "testKF.it not found" );
26        }
27
28        fin >> Dt;
29
30        fin.seek ( "A" );
31        fin >> A;
32        fin.seek ( "B" );
33        fin >> B;
34        fin.seek ( "C" );
35        fin >> C;
36        fin.seek ( "D" );
37        fin >> D;
38        fin.seek ( "R" );
39        fin >> R;
40        fin.seek ( "Q" );
41        fin >> Q;
42        fin.seek ( "P0" );
43        fin >> P0;
44        fin.seek ( "mu0" );
45        fin >> Mu0;
46        mu0 = Mu0.get_col ( 0 );
47
48        vec vQ1 = "0.01:0.04:1";
49        vec vQ2 = "0.01:0.04:1";
50
51        mat LL = zeros ( vQ1.length(),  vQ2.length() );
52
53        Ndat = Dt.cols();
54
55        RV rx ( "{x}", "2" );
56        RV ru ( "{u}", "1" );
57        RV ry ( "{y}", "2" );
58        //
59        //
60        Kalman<ldmat> KFtr;
61
62        for ( int i = 0; i < vQ1.length(); i++ ) {
63
64                for ( int j = 0; j < vQ2.length(); j++ ) {
65                        // KF with R unknown
66                        mat Qj = Q;
67                        Qj ( 0, 0 ) = vQ1 ( i );
68                        Qj ( 1, 1 ) = vQ2 ( j );
69                        KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Qj ) );
70                        KFtr.set_est ( mu0, ldmat ( P0 ) );
71
72                        for ( int t = 1; t < Ndat; t++ ) {
73                                KFtr.bayes ( Dt.get_col ( t ) );
74                                LL ( i, j ) += KFtr._ll();
75                        }
76                }
77        }
78
79        it_file fou ( "testKF_QR_exh.it" );
80
81        fou << Name ( "LL" ) << LL;
82        fou << Name ( "Q1" ) << vQ1;
83        fou << Name ( "Q2" ) << vQ2;
84        //Exit program:
85        return 0;
86
87}
Note: See TracBrowser for help on using the browser.