root/mpdm/merg_pred.cpp @ 264

Revision 254, 4.3 kB (checked in by smidl, 16 years ago)

create namespace bdm

Line 
1#include <estim/arx.h>
2#include <estim/merger.h>
3#include <stat/libEF.h>
4#include <stat/loggers.h>
5//#include <stat/merger.h>
6using namespace bdm;
7
8//These lines are needed for use of cout and endl
9using std::cout;
10using std::endl;
11
12int main() {
13        // Setup model
14        RV y ( "{y }" );
15        RV u1 ( "{u1 }" );
16        RV u2 ( "{u2 }" );
17        RV ym=y; ym.t ( -1 );
18        RV yy = y; yy.add ( ym );
19        RV uu=u1; uu.add ( u2 );
20
21        // Full system
22        vec thg ( "0.7 1 1 0" ); //Simulated system - zero for constant term
23        //y=a y_t-1 + u1 + u2
24        double sqr=0.1;
25        int ord = 1;
26
27        // Estimated systems ARX(2)
28        RV thri ( "{thr_i }",vec_1 ( 3+1 ) );
29        RV thrg ( "{thr_g }",vec_1 ( 4+1 ) );
30        // Setup values
31
32        //ARX constructor
33        mat V0 = 0.001*eye ( thri.count() ); V0 ( 0,0 ) *= 10; //
34        mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; //
35        double nu0 = ord+6.0;
36        double frg = 0.95;
37
38        ARX P1 ( thri, V0, nu0, frg );
39        ARX P2 ( thri, V0, nu0, frg );
40        ARX PG ( thrg, V0g, nu0, frg );
41        //Test estimation
42        int ndat = 200;
43        int t;
44
45        // Logging
46        dirfilelog L ( "exp/merg",10 );
47
48        int Li_Eth1 = L.add ( thri,"P1" );
49        int Li_Eth2 = L.add ( thri,"P2" );
50        int Li_Ethg = L.add ( thrg,"PG" );
51        int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" );
52        int Li_LL   = L.add ( RV ( "{1 2 G }" ), "LL" );
53        int Li_Pred = L.add ( RV ( "{1 2 G ar ge }" ), "Pred" );
54
55        L.init();
56
57        vec Yt ( ndat );
58        vec yt ( 1 );
59
60        Yt.set_subvector ( 0,randn ( ord ) ); //initial values
61        vec rgrg ( thrg.count() -1 ); // constant terms are in!
62        vec rgr1 ( thri.count() -1 );
63        vec rgr2 ( thri.count() -1 );
64
65        vec PredLLs ( 5 );
66        vec PostLLs ( 3 );
67        vec PredLLs_m=zeros ( 5 );
68        ivec ind_r1 = "0 1 3";
69        ivec ind_r2 = "0 2 3";
70        for ( t=0; t<ndat; t++ ) {
71                // True system
72                if ( t>0 ) {
73                        rgrg ( 0 ) =Yt ( t-1 );
74                        rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3);
75                        rgrg ( 2 ) = pow(cos ( ( t/60.0 ) *pi ),3);
76                        rgrg ( 3 ) = 1.0; // constant term
77
78                        rgr1 ( 0 ) = rgrg ( 0 ); rgr1 ( 1 ) = rgrg ( 1 ); rgr1 ( 2 ) = rgrg ( 3 ); // no u2
79                        rgr2 ( 0 ) = rgrg ( 0 ); rgr2 ( 1 ) = rgrg ( 2 ); rgr2 ( 2 ) = rgrg ( 3 ); // no u1
80
81                        Yt ( t ) = thg*rgrg + sqr * NorRNG();
82
83                        // Test predictors
84                        if ( t>2 ) {
85                                        mlnorm<ldmat>* P1p = P1.predictor ( y,concat ( ym,u1 ) );
86                                        mlnorm<ldmat>* P2p = P2.predictor ( y,concat ( ym,u2 ) );
87                                        mlnorm<ldmat>* Pgp = PG.predictor ( y,concat ( ym,uu ) );
88
89                                        Array<mpdf*> A ( 2 ); A ( 0 ) =P1p;A ( 1 ) =P2p;
90                                        merger M ( A );
91                                        enorm<ldmat> g0 ( concat ( yy,uu ) ); g0.set_parameters ( "0 0 0 0 ",3*eye ( 4 ) );
92                                        M.set_parameters ( 1e8, 101,1 );
93                                        M.merge ( &g0 );
94                                yt ( 0 ) = Yt ( t );
95                                double P1pl = P1p->evallogcond ( yt,rgr1 );
96                                double P2pl = P2p->evallogcond ( yt,rgr2 );
97                                double PGpl = Pgp->evallogcond ( yt,rgrg );
98                                {
99                                        cout << "T: " << t << "yt: " << yt << endl;
100                                        cout << "yt_1: " << P1p->_epdf().mean() << endl;
101                                        cout << *P1p <<endl;
102                                        cout << "yt_2: " << P2p->_epdf().mean() << endl;
103                                        cout << "yt_G: " << P2p->_epdf().mean() << endl;
104                                }
105                                double cP1pl;
106                                double cP2pl;
107                                {
108                                        ARX* Apred = ( ARX* ) M._Mix()._Coms ( 0 );
109                                        enorm<ldmat>* MP= Apred->predictor ( concat ( yy,uu ) );
110                                        enorm<ldmat>* mP1p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u1 ) );
111                                        enorm<ldmat>* mP2p = ( enorm<ldmat>* ) MP->marginal ( concat ( yy,u2 ) );
112                                        mlnorm<ldmat>* cP1p = ( mlnorm<ldmat>* ) mP1p->condition ( y );
113                                        mlnorm<ldmat>* cP2p = ( mlnorm<ldmat>* ) mP2p->condition ( y );
114
115                                        cP1pl = cP1p->evallogcond ( yt,rgr1 );
116                                        cP2pl = cP2p->evallogcond ( yt,rgr2 );
117
118                                        cout << "ytm1: " << cP1p->_epdf().mean() << endl;
119                                        cout << "ytm2: " << cP2p->_epdf().mean() << endl;
120                                        cout << *cP2p << endl;
121                                }
122
123                                PredLLs *=frg;
124                                PredLLs += concat ( vec_3 ( P1pl, P2pl, PGpl ), vec_2 ( cP1pl,cP2pl ) );
125                                L.logit ( Li_Pred, PredLLs ); //log-normal
126
127                                delete P1p;
128                                delete P2p;
129                                delete Pgp;
130                        }
131
132                        if ( t<100 ) {
133                                // 1st
134                                P1.bayes ( concat ( Yt ( t ),rgr1 ) );
135                                // 2nd
136                                P2.bayes ( concat ( Yt ( t ),rgr2 ) );
137
138                                //Global
139                                PG.bayes ( concat ( Yt ( t ),rgrg ) );
140                        }
141                        //Merger
142                }
143                L.logit ( Li_Eth1, P1._epdf().mean() );
144                L.logit ( Li_Eth2, P2._epdf().mean() );
145                L.logit ( Li_Ethg, PG._epdf().mean() );
146                L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) );
147                PostLLs *= frg;
148                PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() );
149                L.logit ( Li_LL, PostLLs );
150                L.step ( );
151        }
152        L.finalize( );
153        L.itsave ( "merg.it" );
154}
Note: See TracBrowser for help on using the browser.