- Timestamp:
- 11/10/08 15:40:30 (16 years ago)
- Location:
- mpdm
- Files:
-
- 1 added
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
mpdm/merg_giw.cpp
r204 r205 17 17 RV uu=u1; uu.add ( u2 ); 18 18 19 double a1t = 1.5; 20 double a2t = 0.8; 21 double sqr=0.10; 19 22 // Full system 20 vec thg ( "1 1" ); //Simulated system - zero for constant term 21 double sqr=1.0; 23 vec thg =vec_2 ( a1t,a2t ); //Simulated system - zero for constant term 22 24 vec Th = concat ( thg, sqr ); //Full parameter 23 25 … … 28 30 RV all =a1; all.add ( a2 ); all.add ( r ); 29 31 RV allj =a1; allj.add ( r ); allj.add ( a2 ); 30 vec Thj ("1 1 1");32 vec Thj=vec_3 ( a1t,sqr,a2t ); 31 33 // Setup values 32 34 … … 37 39 ARX P1 ( concat ( a1,r ), V0, -1 ); 38 40 ARX P2 ( concat ( a2,r ), V0, -1 ); 39 ARX PG ( all, V0g, -1 ); 41 ARX PG ( all, V0g, -1 ); //or -1? 42 // ARX PGk ( all, V0g, -1 ); 40 43 41 44 //Test estimation … … 44 47 45 48 // Logging 46 dirfilelog L ( "exp/merg_giw", 1);49 dirfilelog L ( "exp/merg_giw",ndat ); 47 50 int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" ); 48 int Li_LL = L.add ( RV ( "{1 2 G M }" ), "LL" ); 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" ); 49 54 int Li_Gm = L.add ( RV ( "{a1 a2 r }" ), "G" ); 50 55 int Li_Mm = L.add ( RV ( "{a1 r a2 }" ), "M" ); … … 54 59 vec yt ( 1 ); 55 60 56 vec LLs ( 4);61 vec LLs ( 2 ); 57 62 vec rgrg ( 2 ); 58 63 59 64 //Proposal 60 enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1 ",mat ("1") );65 enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1 ",mat ( "1" ) ); 61 66 egamma g1 ( r ); g1.set_parameters ( "2 ", "2" ); 62 enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1 ",mat ("1") );67 enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1 ",mat ( "1" ) ); 63 68 64 Array<const epdf*> A (3); A(0) = &g0; A(1)=&g1; A(2) = &g2;65 eprod G0 (A);66 69 Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ( 2 ) = &g2; 70 eprod G0 ( A ); 71 67 72 for ( t=0; t<ndat; t++ ) { 68 73 // True system 69 74 rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 ); 70 rgrg ( 1 ) = pow ( cos ( ( t/40.0 ) *pi ),3 );75 rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 ); 71 76 72 77 Yt ( t ) = thg*rgrg + sqr * NorRNG(); 73 78 74 79 // Bayes for all 75 P1.bayes ( concat ( Yt ( t ),vec_1 (rgrg ( 0 ) )) );76 P2.bayes ( concat ( Yt ( t ),vec_1 (rgrg ( 1 ) )) );80 P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) ); 81 P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) ); 77 82 PG.bayes ( concat ( Yt ( t ),rgrg ) ); 78 83 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 79 93 // Merge estimates 80 mepdf eG1 (P1._e());81 mepdf eG2 (P2._e());94 mepdf eG1 ( P1._e() ); 95 mepdf eG2 ( P2._e() ); 82 96 Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2; 83 97 merger M ( A ); 84 M.set_parameters ( 10.0, 100,2 ); 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); 85 100 M.merge ( &G0 ); 101 //M2.merge ( &G0 ); 86 102 87 103 //Likelihood 88 104 yt ( 0 ) = Yt ( t ); 89 105 90 LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") ); 91 LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2") ); 92 LLs ( 2 ) = PG._e()->evalpdflog (Th ); 93 LLs ( 3 ) = M._Mix().logpred ( concat(Thj, vec_1(1.0)) ); 106 // LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") ); 107 // LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2") ); 108 LLs ( 0 ) = PG._e()->evalpdflog ( Th ); 109 LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) ); 110 // LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) ); 94 111 L.logit ( Li_LL, LLs ); //log-normal 95 112 96 113 //Logger 97 114 L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) ); 115 L.logit ( Li_P1m, P1._e()->mean() ); 116 L.logit ( Li_P2m, P2._e()->mean() ); 98 117 L.logit ( Li_Gm, PG._e()->mean() ); 99 emix *tm =M._Mix().predictor(allj); 100 L.logit ( Li_Mm, tm->mean() ); 101 delete tm; 118 L.logit ( Li_Mm, M.mean() ); 119 102 120 L.step ( ); 121 122 cout << "Vg: " << PG._e()->_V().to_mat() <<endl; 123 vec mea = M.mean(); 124 cout << "Ve: " << M.variance() <<endl; 103 125 } 104 126 L.finalize( ); 105 127 L.itsave ( "merg_egiw.it" ); 128 cout << endl; 106 129 }