Changeset 205 for mpdm/merg_giw.cpp

Show
Ignore:
Timestamp:
11/10/08 15:40:30 (16 years ago)
Author:
smidl
Message:

merger posledni verze

Files:
1 modified

Legend:

Unmodified
Added
Removed
  • mpdm/merg_giw.cpp

    r204 r205  
    1717        RV uu=u1; uu.add ( u2 ); 
    1818 
     19        double a1t = 1.5; 
     20        double a2t = 0.8; 
     21        double sqr=0.10; 
    1922        // 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 
    2224        vec Th = concat ( thg, sqr ); //Full parameter 
    2325 
     
    2830        RV all =a1; all.add ( a2 ); all.add ( r ); 
    2931        RV allj =a1; allj.add ( r ); allj.add ( a2 ); 
    30         vec Thj("1 1 1"); 
     32        vec Thj=vec_3 ( a1t,sqr,a2t ); 
    3133        // Setup values 
    3234 
     
    3739        ARX P1 ( concat ( a1,r ), V0, -1 ); 
    3840        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 ); 
    4043 
    4144        //Test estimation 
     
    4447 
    4548        // Logging 
    46         dirfilelog L ( "exp/merg_giw",1 ); 
     49        dirfilelog L ( "exp/merg_giw",ndat ); 
    4750        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" ); 
    4954        int Li_Gm   = L.add ( RV ( "{a1 a2 r }" ), "G" ); 
    5055        int Li_Mm   = L.add ( RV ( "{a1 r a2 }" ), "M" ); 
     
    5459        vec yt ( 1 ); 
    5560 
    56         vec LLs ( 4 ); 
     61        vec LLs ( 2 ); 
    5762        vec rgrg ( 2 ); 
    58          
     63 
    5964        //Proposal 
    60         enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1  ",mat("1") ); 
     65        enorm<ldmat> g0 ( a1 ); g0.set_parameters ( "1  ",mat ( "1" ) ); 
    6166        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" ) ); 
    6368 
    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 
    6772        for ( t=0; t<ndat; t++ ) { 
    6873                // True system 
    6974                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 ); 
    7176 
    7277                Yt ( t ) = thg*rgrg + sqr * NorRNG(); 
    7378 
    7479                // 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 ) ) ) ); 
    7782                PG.bayes ( concat ( Yt ( t ),rgrg ) ); 
    7883 
     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 
    7993                // Merge estimates 
    80                 mepdf eG1(P1._e()); 
    81                 mepdf eG2(P2._e()); 
     94                mepdf eG1 ( P1._e() ); 
     95                mepdf eG2 ( P2._e() ); 
    8296                Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2; 
    8397                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); 
    85100                M.merge ( &G0 ); 
     101                //M2.merge ( &G0 ); 
    86102 
    87103                //Likelihood 
    88104                yt ( 0 ) = Yt ( t ); 
    89105 
    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)) ); 
    94111                L.logit ( Li_LL, LLs ); //log-normal 
    95112 
    96113                //Logger 
    97114                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() ); 
    98117                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 
    102120                L.step ( ); 
     121 
     122                cout << "Vg: " << PG._e()->_V().to_mat() <<endl; 
     123                vec mea = M.mean(); 
     124                cout << "Ve: " << M.variance() <<endl; 
    103125        } 
    104126        L.finalize( ); 
    105127        L.itsave ( "merg_egiw.it" ); 
     128        cout << endl; 
    106129}