[204] | 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> |
---|
[254] | 6 | using namespace bdm; |
---|
[204] | 7 | |
---|
| 8 | //These lines are needed for use of cout and endl |
---|
| 9 | using std::cout; |
---|
| 10 | using std::endl; |
---|
| 11 | |
---|
| 12 | int main() { |
---|
| 13 | // Setup model |
---|
| 14 | RV y ( "{y }" ); |
---|
| 15 | RV u1 ( "{u1 }" ); |
---|
| 16 | RV u2 ( "{u2 }" ); |
---|
| 17 | RV uu=u1; uu.add ( u2 ); |
---|
| 18 | |
---|
[205] | 19 | double a1t = 1.5; |
---|
| 20 | double a2t = 0.8; |
---|
| 21 | double sqr=0.10; |
---|
[204] | 22 | // Full system |
---|
[205] | 23 | vec thg =vec_2 ( a1t,a2t ); //Simulated system - zero for constant term |
---|
[204] | 24 | vec Th = concat ( thg, sqr ); //Full parameter |
---|
| 25 | |
---|
| 26 | // Estimated systems ARX(2) |
---|
| 27 | RV a1 ( "{a1 }" ); |
---|
| 28 | RV a2 ( "{a2 }" ); |
---|
| 29 | RV r ( "{r }" ); |
---|
| 30 | RV all =a1; all.add ( a2 ); all.add ( r ); |
---|
| 31 | RV allj =a1; allj.add ( r ); allj.add ( a2 ); |
---|
[205] | 32 | vec Thj=vec_3 ( a1t,sqr,a2t ); |
---|
[204] | 33 | // Setup values |
---|
| 34 | |
---|
| 35 | //ARX constructor |
---|
| 36 | mat V0 = 0.001*eye ( 2 ); V0 ( 0,0 ) = 1; // |
---|
| 37 | mat V0g = 0.001*eye ( 3 ); V0g ( 0,0 ) = 1; // |
---|
| 38 | |
---|
[270] | 39 | ARX P1; P1.set_statistics(2, V0, -1 ); |
---|
| 40 | ARX P2; P2.set_statistics(2, V0, -1 ); |
---|
| 41 | ARX PG; PG.set_statistics(3, V0g, -1 ); //or -1? |
---|
[205] | 42 | // ARX PGk ( all, V0g, -1 ); |
---|
[204] | 43 | |
---|
| 44 | //Test estimation |
---|
| 45 | int ndat = 100; |
---|
| 46 | int t; |
---|
| 47 | |
---|
| 48 | // Logging |
---|
[205] | 49 | dirfilelog L ( "exp/merg_giw",ndat ); |
---|
[204] | 50 | int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" ); |
---|
[205] | 51 | int Li_LL = L.add ( RV ( "{G M }" ), "LL" ); |
---|
| 52 | int Li_P1m = L.add ( RV ( "{a1 r }" ), "P1" ); |
---|
| 53 | int Li_P2m = L.add ( RV ( "{a2 r }" ), "P2" ); |
---|
[204] | 54 | int Li_Gm = L.add ( RV ( "{a1 a2 r }" ), "G" ); |
---|
| 55 | int Li_Mm = L.add ( RV ( "{a1 r a2 }" ), "M" ); |
---|
| 56 | L.init(); |
---|
| 57 | |
---|
| 58 | vec Yt ( ndat ); |
---|
| 59 | vec yt ( 1 ); |
---|
| 60 | |
---|
[205] | 61 | vec LLs ( 2 ); |
---|
[204] | 62 | vec rgrg ( 2 ); |
---|
[205] | 63 | |
---|
[204] | 64 | //Proposal |
---|
[270] | 65 | enorm<ldmat> g0; g0.set_rv( a1 ); g0.set_parameters ( "1 ",mat ( "1" ) ); |
---|
| 66 | egamma g1; g1.set_rv ( r ); g1.set_parameters ( "2 ", "2" ); |
---|
| 67 | enorm<ldmat> g2; g2.set_rv ( a2 ); g2.set_parameters ( "1 ",mat ( "1" ) ); |
---|
[204] | 68 | |
---|
[205] | 69 | Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ( 2 ) = &g2; |
---|
| 70 | eprod G0 ( A ); |
---|
| 71 | |
---|
[204] | 72 | for ( t=0; t<ndat; t++ ) { |
---|
| 73 | // True system |
---|
| 74 | rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 ); |
---|
[205] | 75 | rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 ); |
---|
[204] | 76 | |
---|
| 77 | Yt ( t ) = thg*rgrg + sqr * NorRNG(); |
---|
| 78 | |
---|
| 79 | // Bayes for all |
---|
[205] | 80 | P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) ); |
---|
| 81 | P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) ); |
---|
[204] | 82 | PG.bayes ( concat ( Yt ( t ),rgrg ) ); |
---|
| 83 | |
---|
[205] | 84 | |
---|
| 85 | // crippling PGk by substituting zeros. |
---|
| 86 | /* ldmat &Vk=const_cast<egiw*>(PGk._e())->_V(); //PG ldmat does not like 0! |
---|
| 87 | mat fVk=PG._e()->_V().to_mat(); |
---|
| 88 | fVk(1,2) = 0.0; |
---|
| 89 | fVk(2,1) = 0.0; |
---|
| 90 | Vk = ldmat(fVk); |
---|
| 91 | */ //PGk is now krippled |
---|
| 92 | |
---|
[204] | 93 | // Merge estimates |
---|
[205] | 94 | mepdf eG1 ( P1._e() ); |
---|
| 95 | mepdf eG2 ( P2._e() ); |
---|
[204] | 96 | Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2; |
---|
| 97 | merger M ( A ); |
---|
[205] | 98 | M.set_parameters ( 1.5, 100,3 ); //M._Mix().set_method(QB); |
---|
| 99 | //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB); |
---|
[204] | 100 | M.merge ( &G0 ); |
---|
[205] | 101 | //M2.merge ( &G0 ); |
---|
[204] | 102 | |
---|
| 103 | //Likelihood |
---|
| 104 | yt ( 0 ) = Yt ( t ); |
---|
| 105 | |
---|
[211] | 106 | // LLs ( 0 ) = P1._e()->evallog ( get_vec(Th, "1 2") ); |
---|
| 107 | // LLs ( 1 ) = P2._e()->evallog ( get_vec(Th, "3 2") ); |
---|
| 108 | LLs ( 0 ) = PG._e()->evallog ( Th ); |
---|
[205] | 109 | LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) ); |
---|
| 110 | // LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) ); |
---|
[204] | 111 | L.logit ( Li_LL, LLs ); //log-normal |
---|
| 112 | |
---|
| 113 | //Logger |
---|
| 114 | L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) ); |
---|
[205] | 115 | L.logit ( Li_P1m, P1._e()->mean() ); |
---|
| 116 | L.logit ( Li_P2m, P2._e()->mean() ); |
---|
[204] | 117 | L.logit ( Li_Gm, PG._e()->mean() ); |
---|
[205] | 118 | L.logit ( Li_Mm, M.mean() ); |
---|
| 119 | |
---|
[204] | 120 | L.step ( ); |
---|
[205] | 121 | |
---|
| 122 | cout << "Vg: " << PG._e()->_V().to_mat() <<endl; |
---|
| 123 | vec mea = M.mean(); |
---|
| 124 | cout << "Ve: " << M.variance() <<endl; |
---|
[204] | 125 | } |
---|
| 126 | L.finalize( ); |
---|
| 127 | L.itsave ( "merg_egiw.it" ); |
---|
[205] | 128 | cout << endl; |
---|
[204] | 129 | } |
---|