root/library/tests/stresssuite/kalman_QRexh_stress.cpp @ 1033

Revision 721, 1.6 kB (checked in by mido, 15 years ago)

stresssuite - halfway point

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