root/library/tests/stresssuite/kalman_QR_stress.cpp @ 1137

Revision 1064, 2.2 kB (checked in by mido, 14 years ago)

astyle applied all over the library

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_QR_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    if ( !xxx ) {
25        bdm_error ( "kalman_stress.it not found" );
26    }
27    fin >> Dt;
28    fin.seek ( "A" );
29    fin >> A;
30    fin.seek ( "B" );
31    fin >> B;
32    fin.seek ( "C" );
33    fin >> C;
34    fin.seek ( "D" );
35    fin >> D;
36    fin.seek ( "R" );
37    fin >> R;
38    fin.seek ( "Q" );
39    fin >> Q;
40    fin.seek ( "P0" );
41    fin >> P0;
42    fin.seek ( "mu0" );
43    fin >> Mu0;
44    mu0 = Mu0.get_col ( 0 );
45
46    Ndat = Dt.cols();
47    XQRt = zeros ( 5, Ndat );
48    mat Xt = zeros ( 2, Ndat );
49
50//      cout << KF;
51    RV rx ( "{x }", "2" );
52    RV ru ( "{u }", "1" );
53    RV ry ( "{y }", "1" );
54    RV rQR ( "{Q,R }", "3" );
55    //
56    KFcondQR KF;
57    // KF with R unknown
58    KF.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) );
59    KF.set_est ( mu0, ldmat ( P0 ) );
60    //
61    Kalman<ldmat> KFtr;
62    // KF with R unknown
63    KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) );
64    KFtr.set_est ( mu0, ldmat ( P0 ) );
65
66    mgamma evolQR;
67    evolQR.set_parameters ( 10.0, "1 1 1" ); //sigma = 1/10 mu
68
69    MPF<KFcondQR> KF_QR;
70    KF_QR.set_parameters ( &evolQR, &evolQR, 100 );
71    evolQR.condition ( "1 1 1" );
72    egamma kfinit = egamma ( "10 10 10", "1 1 1" );
73    KF_QR.set_statistics ( egamma ( "10 10 10", "1 1 1" ), &KF );
74    const epdf& mpost = KF_QR.posterior();
75    const epdf& mposttr = KFtr.posterior();
76
77    XQRt.set_col ( 0, mpost.mean() );
78    Xt.set_col ( 0, mposttr.mean() );
79    for ( int t = 1; t < Ndat; t++ ) {
80        KF_QR.bayes ( Dt.get_col ( t ) );
81        KFtr.bayes ( Dt.get_col ( t ) );
82
83        XQRt.set_col ( t, mpost.mean() );
84        Xt.set_col ( t, mposttr.mean() );
85    }
86
87    it_file fou ( "kalman_stress_QR_res.it" );
88    fou << Name ( "xqrth" ) << XQRt;
89    fou << Name ( "xth" ) << Xt;
90}
Note: See TracBrowser for help on using the browser.