- Timestamp:
- 11/25/08 12:24:42 (16 years ago)
- Location:
- mpdm
- Files:
-
- 3 modified
-
CMakeLists.txt (modified) (1 diff)
-
merg_2a.cpp (modified) (3 diffs)
-
merg_pred.cpp (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
mpdm/CMakeLists.txt
r204 r213 7 7 ## Save all needed libraries in variable BdmLibs 8 8 #SET(BdmLibs bdm itpp_debug) 9 SET(BdmLibs bdm itpp)10 SET(ITppLibs itpp)11 9 12 10 IF(WIN32) 13 SET(BdmLibs ${BdmLibs} libacml_dll) 14 SET(ITppLibs ${ITppLibs} libacml_dll) 11 SET(AddLib libacml_dll) 15 12 ENDIF(WIN32) 13 14 # Define macro for testing a file 15 MACRO(TEST FN) 16 add_executable (${FN} ${FN}.cpp) 17 target_link_libraries (${FN} debug itpp_debug) 18 target_link_libraries (${FN} optimized itpp) 19 target_link_libraries (${FN} bdm ${AddLib}) 20 ENDMACRO(TEST) 16 21 17 22 # Add executable called "helloDemo" that is built from the source files 18 23 # "demo.cxx" and "demo_b.cxx". The extensions are automatically found. 19 add_executable (merg_pred merg_pred.cpp) 20 target_link_libraries (merg_pred ${BdmLibs}) 21 22 add_executable (merg_giw merg_giw.cpp) 23 target_link_libraries (merg_giw ${BdmLibs}) 24 25 add_executable (merger_iter_cond merger_iter_cond.cpp) 26 target_link_libraries (merger_iter_cond ${BdmLibs}) 24 TEST(merg_pred) 25 TEST(merg_giw) 26 TEST(merger_iter_cond) 27 TEST(merg_2a) -
mpdm/merg_2a.cpp
r211 r213 10 10 using std::endl; 11 11 12 /*! The purpose of this experiment is to test merging fragmental pdfs between two participants sharing the same parameter, a. However, this parameter has a different role in each model. 13 14 Lets assume that: 15 u -> |a,r| -> y -> |a,r| -> z 16 */ 17 12 18 int main() { 13 19 // Setup model 14 20 RV y ( "{y }" ); 15 21 RV u ( "{u }" ); 16 RV um = u; um.t(-1);17 22 RV z ( "{z }" ); 18 RV a ("{a }"); 19 RV b ("{b }"); 20 RV c ("{c }"); 23 RV a ("{a }"); 24 RV b ("{b }"); RV ab = a; ab.add(b); 25 RV c ("{c }"); RV ac = a; ac.add(c); 21 26 RV r ("{r }"); 22 27 23 28 double at = 1.5; 24 double bt = 0. 8;25 double ct = 0.50;26 double sig = 0.10;29 double bt = 0.5; 30 double ct = -0.5; 31 double rt = 0.30; 27 32 // Full system 28 33 vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term 29 vec th u=vec_2 ( at,ct ); //Simulated system - zero for constant term30 vec Thy = concat ( thy, vec_1( sig) ); //Full parameter31 vec Th u = concat ( thu, vec_1(sig) ); //Full parameter34 vec thz =vec_2 ( at,ct ); //Simulated system - zero for constant term 35 vec Thy = concat ( thy, vec_1(rt) ); //Full parameter 36 vec Thz = concat ( thz, vec_1(rt) ); //Full parameter 32 37 33 38 //ARX constructor 34 39 mat V0 = 0.001*eye ( 3 ); V0 ( 0,0 ) = 1; // 35 40 36 ARX P1 ( concat ( a , concat(b,r)), V0, -1 );37 ARX P2 ( concat ( a , concat(c,r)), V0, -1 );41 ARX P1 ( concat ( ab,r ), V0, -1 ); 42 ARX P2 ( concat ( ac,r ), V0, -1 ); 38 43 39 44 //Test estimation … … 42 47 43 48 // Logging 44 dirfilelog L ( "exp/merg_2a",ndat ); 45 int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" ); 46 int Li_LL = L.add ( RV ( "{G M }" ), "LL" ); 47 int Li_P1m = L.add ( RV ( "{a b r }" ), "P1" ); 48 int Li_P2m = L.add ( RV ( "{a c r }" ), "P2" ); 49 dirfilelog L ( "exp/merg_2a",3 ); 50 int Li_Data = L.add ( RV ( "{U Y Z }" ), "" ); 51 // int Li_LL = L.add ( RV ( "{P1 P2 M1 M2 }" ), "LL" ); 52 int Li_P1m = L.add ( concat ( ab,r ), "P1m" ); 53 int Li_P2m = L.add ( concat ( ac,r ), "P2m" ); 54 int Li_Mm = L.add ( concat ( ab,concat(r,c) ), "Mm" ); 49 55 L.init(); 50 56 57 vec Ut ( ndat ); 51 58 vec Yt ( ndat ); 59 vec Zt ( ndat ); 52 60 vec yt ( 1 ); 53 61 54 62 //Proposal 55 enorm<ldmat> g0 ( a 1 ); g0.set_parameters ( "1 ",mat ( "1" ) );63 enorm<ldmat> g0 ( ab ); g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) ); 56 64 egamma g1 ( r ); g1.set_parameters ( "2 ", "2" ); 57 enorm<ldmat> g2 ( a2 ); g2.set_parameters ( "1",mat ( "1" ) );65 enorm<ldmat> g2 ( c ); g2.set_parameters ( "1 ",mat ( "1" ) ); 58 66 59 Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A ( 2 ) = &g2;67 Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2; 60 68 eprod G0 ( A ); 61 69 62 for ( t=0; t<ndat; t++ ) { 70 vec rgru(2); 71 vec rgry(2); 72 Yt(0) = 0.1; 73 Ut(0) = 0.0; 74 for ( t=1; t<ndat; t++ ) { 63 75 // True system 64 rgrg ( 0 ) = pow ( sin ( ( t/40.0 ) *pi ),3 ); 65 rgrg ( 1 ) = pow ( cos ( ( t/40.0 +0.1 ) *pi ),3 ); 66 67 Yt ( t ) = thg*rgrg + sqr * NorRNG(); 76 Ut ( t ) = pow ( sin ( ( t/40.0 ) *pi ),3 ); 77 rgru(0) = Ut(t); rgru(1) = Ut(t-1); 78 Yt ( t ) = thy*rgru + rt * NorRNG(); 79 rgry(0) = Yt(t); rgry(1) = Yt(t-1); 80 Zt ( t ) = thz*rgry + rt * NorRNG(); 68 81 69 82 // Bayes for all 70 P1.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 0 ) ) ) ); 71 P2.bayes ( concat ( Yt ( t ),vec_1 ( rgrg ( 1 ) ) ) ); 72 PG.bayes ( concat ( Yt ( t ),rgrg ) ); 73 74 75 // crippling PGk by substituting zeros. 76 /* ldmat &Vk=const_cast<egiw*>(PGk._e())->_V(); //PG ldmat does not like 0! 77 mat fVk=PG._e()->_V().to_mat(); 78 fVk(1,2) = 0.0; 79 fVk(2,1) = 0.0; 80 Vk = ldmat(fVk); 81 */ //PGk is now krippled 83 P1.bayes ( concat ( Yt ( t ),rgru ) ); 84 P2.bayes ( concat ( Zt ( t ),rgry ) ); 82 85 83 86 // Merge estimates … … 89 92 //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB); 90 93 M.merge ( &G0 ); 91 //M2.merge ( &G0 );92 94 93 95 //Likelihood 94 96 yt ( 0 ) = Yt ( t ); 95 97 96 // LLs ( 0 ) = P1._e()->evalpdflog ( get_vec(Th, "1 2") );97 // LLs ( 1 ) = P2._e()->evalpdflog ( get_vec(Th, "3 2") );98 LLs ( 0 ) = PG._e()->evalpdflog ( Th );99 LLs ( 1 ) = M._Mix().logpred ( concat ( Thj, vec_1 ( 1.0 ) ) );100 // LLs ( 3 ) = M2._Mix().logpred ( concat(Thj, vec_1(1.0)) );101 L.logit ( Li_LL, LLs ); //log-normal102 103 98 //Logger 104 L.logit ( Li_Data, vec_3 ( Yt ( t ), rgrg ( 0 ), rgrg ( 1 ) ) ); 105 L.logit ( Li_P1m, P1._e()->mean() ); 106 L.logit ( Li_P2m, P2._e()->mean() ); 107 L.logit ( Li_Gm, PG._e()->mean() ); 108 L.logit ( Li_Mm, M.mean() ); 109 99 L.logit(Li_Data, vec_3(Ut(t), Yt(t), Zt(t))); 100 L.logit(Li_P1m, P1._e()->mean()); 101 L.logit(Li_P2m, P2._e()->mean()); 102 L.logit(Li_Mm, M.mean()); 110 103 L.step ( ); 111 104 112 cout << "Vg: " << PG._e()->_V().to_mat() <<endl;113 vec mea = M.mean();114 cout << "Ve: " << M.variance() <<endl;115 105 } 116 106 L.finalize( ); 117 L.itsave ( "merg_ egiw.it" );107 L.itsave ( "merg_2a.it" ); 118 108 cout << endl; 119 109 } -
mpdm/merg_pred.cpp
r211 r213 15 15 RV u1 ( "{u1 }" ); 16 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);17 RV ym=y; ym.t ( -1 ); 18 RV yy = y; yy.add ( ym ); 19 RV uu=u1; uu.add ( u2 ); 20 20 21 21 // Full system … … 44 44 45 45 // Logging 46 dirfilelog L ( "exp/merg", ndat);46 dirfilelog L ( "exp/merg",10 ); 47 47 48 48 int Li_Eth1 = L.add ( thri,"P1" ); … … 56 56 57 57 vec Yt ( ndat ); 58 vec yt (1);58 vec yt ( 1 ); 59 59 60 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 );61 vec rgrg ( thrg.count() -1 ); // constant terms are in! 62 vec rgr1 ( thri.count() -1 ); 63 vec rgr2 ( thri.count() -1 ); 64 64 65 vec PredLLs (5);66 vec PostLLs (3);67 vec PredLLs_m=zeros (5);65 vec PredLLs ( 5 ); 66 vec PostLLs ( 3 ); 67 vec PredLLs_m=zeros ( 5 ); 68 68 ivec ind_r1 = "0 1 3"; 69 69 ivec ind_r2 = "0 2 3"; … … 73 73 rgrg ( 0 ) =Yt ( t-1 ); 74 74 rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3); 75 rgrg ( 2 ) = pow(cos ( ( t/ 40.0 ) *pi ),3);76 rgrg ( 3) = 1.0; // constant term77 78 rgr1 (0) = rgrg(0); rgr1(1) = rgrg(1); rgr1(2) = rgrg(3); // no u279 rgr2 (0) = rgrg(0); rgr2(1) = rgrg(2); rgr2(2) = rgrg(3); // no u175 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 80 81 81 Yt ( t ) = thg*rgrg + sqr * NorRNG(); 82 82 83 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(10000000.0, 100,1); 93 M.merge(&g0); 94 95 yt(0) = Yt(t); 96 double P1pl = P1p->evallogcond(yt,rgr1); 97 double P2pl = P2p->evallogcond(yt,rgr2); 98 double PGpl = Pgp->evallogcond(yt,rgrg); 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 ); 99 98 { 100 cout << " yt: " << yt << endl;99 cout << "T: " << t << "yt: " << yt << endl; 101 100 cout << "yt_1: " << P1p->_epdf().mean() << endl; 101 cout << *P1p <<endl; 102 102 cout << "yt_2: " << P2p->_epdf().mean() << endl; 103 103 cout << "yt_G: " << P2p->_epdf().mean() << endl; … … 106 106 double cP2pl; 107 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);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 114 115 cP1pl = cP1p->evallogcond (yt,rgr1);116 cP2pl = cP2p->evallogcond (yt,rgr2);117 115 cP1pl = cP1p->evallogcond ( yt,rgr1 ); 116 cP2pl = cP2p->evallogcond ( yt,rgr2 ); 117 118 118 cout << "ytm1: " << cP1p->_epdf().mean() << endl; 119 119 cout << "ytm2: " << cP2p->_epdf().mean() << endl; 120 cout << *cP2p << endl; 120 121 } 121 122 122 123 PredLLs *=frg; 123 PredLLs += log(concat(vec_3(P1pl, P2pl, PGpl), vec_2(cP1pl,cP2pl)));124 L.logit (Li_Pred, PredLLs); //log-normal125 124 PredLLs += concat ( vec_3 ( P1pl, P2pl, PGpl ), vec_2 ( cP1pl,cP2pl ) ); 125 L.logit ( Li_Pred, PredLLs ); //log-normal 126 126 127 delete P1p; 127 128 delete P2p; 128 129 delete Pgp; 129 130 } 130 131 // 1st132 P1.bayes ( concat(Yt(t),rgr1) );133 // 2nd134 P2.bayes ( concat(Yt(t),rgr2) );135 131 136 //Global 137 PG.bayes ( concat ( Yt ( t ),rgrg ) ); 138 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 } 139 141 //Merger 140 142 } … … 146 148 PostLLs += vec_3 ( P1._ll(), P2._ll(), PG._ll() ); 147 149 L.logit ( Li_LL, PostLLs ); 148 L.step ( 150 L.step ( ); 149 151 } 150 152 L.finalize( );
