Changeset 162 for mpdm

Show
Ignore:
Timestamp:
09/04/08 20:27:01 (16 years ago)
Author:
smidl
Message:

opravy a dokumentace

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • mpdm/merg_pred.cpp

    r161 r162  
    22#include <stat/libEF.h> 
    33#include <stat/loggers.h> 
     4#include <stat/merger.h> 
    45using namespace itpp; 
    56 
    … …  
    1011int main() { 
    1112        // Setup model 
    12         RV y("1","{y }","1","0"); 
    13         RV u1("2","{u1 }","1","0"); 
    14         RV u2("3","{u2 }","1","0"); 
    15          
     13        RV y ( "{y }" ); 
     14        RV u1 ( "{u1 }" ); 
     15        RV u2 ( "{u2 }" ); 
     16 
    1617        // Full system 
    17         vec thg("0.7 1 1"); //Simulated system  
    18                                                 //y=a y_t-1 + u1 + u2 
     18        vec thg ( "0.7 1 1" ); //Simulated system 
     19        //y=a y_t-1 + u1 + u2 
    1920        double sqr=0.1; 
    2021        int ord = 1; 
    21          
     22 
    2223        // Estimated systems ARX(2) 
    23         RV thri("4","{thr }",vec_1(2+1),"0"); 
    24         RV thrg("5","{thr }",vec_1(3+1),"0"); 
     24        RV thri ( "{thr_i }",vec_1 ( 2+1 ) ); 
     25        RV thrg ( "{thr_g }",vec_1 ( 3+1 ) ); 
    2526        // Setup values 
    26          
     27 
    2728        //ARX constructor 
    28         mat V0 = 0.001*eye(thri.count()); V0(0,0)*= 10; // 
    29         mat V0g = 0.001*eye(thrg.count()); V0g(0,0)*= 10; // 
     29        mat V0 = 0.001*eye ( thri.count() ); V0 ( 0,0 ) *= 10; // 
     30        mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; // 
    3031        double nu0 = ord+1; 
    31         double frg = 0.94; 
    32          
    33         ARX P1(thri, V0, nu0, frg); 
    34         ARX P2(thri, V0, nu0, frg); 
    35         ARX PG(thrg, V0g, nu0, frg); 
     32        double frg = 0.99; 
     33 
     34        ARX P1 ( thri, V0, nu0, frg ); 
     35        ARX P2 ( thri, V0, nu0, frg ); 
     36        ARX PG ( thrg, V0g, nu0, frg ); 
    3637        //Test estimation 
    37         int ndat = 1000; 
     38        int ndat = 10000; 
    3839        int t; 
    39          
     40 
    4041        // Logging 
    41         dirfilelog L("exp/merg",1000); 
     42        dirfilelog L ( "exp/merg",ndat ); 
    4243 
    43         int Eth1_log = L.add(thri,"P1"); 
    44         int Eth2_log = L.add(thri,"P2"); 
    45         int Data_log = L.add(RV("10","{Y U1 U2 }","3","0"), ""); 
    46         int LL_log   = L.add(RV("11","{1 2 G }","3","0"), "LL"); 
    47          
    48         vec Yt(ndat); 
    49          
    50         Yt.set_subvector(0,randn(ord)); //initial values 
    51         vec rgrg(thg.length()); 
    52         vec rgri(thri.count()); 
    53          
    54         for (t=2; t<ndat; t++) { 
     44        int Eth1_log = L.add ( thri,"P1" ); 
     45        int Eth2_log = L.add ( thri,"P2" ); 
     46        int Ethg_log = L.add ( thrg,"PG" ); 
     47        int Data_log = L.add ( RV ( "{Y U1 U2 }" ), "" ); 
     48        int LL_log   = L.add ( RV ( "{1 2 G }" ), "LL" ); 
     49 
     50        L.init(); 
     51 
     52        vec Yt ( ndat ); 
     53 
     54        Yt.set_subvector ( 0,randn ( ord ) ); //initial values 
     55        vec rgrg ( thg.length() ); 
     56        vec rgri ( thri.count() ); 
     57 
     58        for ( t=0; t<ndat; t++ ) { 
    5559                // True system 
    56                 rgrg(0)=Yt(t-1);  
    57                 rgrg(1)= sin(t/(40)*pi); 
    58                 rgrg(2)= cos(t/(40)*pi); 
    59                  
    60                 Yt(t) = thg*rgrg + sqr * NorRNG(); 
    61                  
    62                 // 1st 
    63                 rgri(0)=Yt(t); 
    64                 rgri(1)=rgrg(1); 
    65                 P1.bayes(rgri); 
    66                 // 2nd 
    67                 rgri(1)=rgrg(2); 
    68                 P2.bayes(rgri); 
    69                  
    70                 //Global 
    71                 PG.bayes(concat(Yt(t),rgrg)); 
    72                 L.logit(Eth1_log, P1._epdf().mean()); 
    73                 L.logit(Eth2_log, P2._epdf().mean()); 
    74                 L.logit(Data_log, vec_3(Yt(t), rgrg(1), rgrg(2))); 
    75                 L.logit(LL_log, vec_3(P1._ll(), P2._ll(), PG._ll())); 
    76                 L.step(false); 
     60                if ( t>0 ) { 
     61                        rgrg ( 0 ) =Yt ( t-1 ); 
     62                        rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3); 
     63                        rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3); 
     64 
     65                        Yt ( t ) = thg*rgrg + sqr * NorRNG(); 
     66 
     67                        // 1st 
     68                        rgri ( 0 ) =Yt ( t ); 
     69                        rgri ( 1 ) =Yt ( t-1 ); 
     70                        rgri ( 2 ) =rgrg ( 1 ); 
     71                        P1.bayes ( rgri ); 
     72                        // 2nd 
     73                        rgri ( 2 ) =rgrg ( 2 ); 
     74                        P2.bayes ( rgri ); 
     75 
     76                        //Global 
     77                        PG.bayes ( concat ( Yt ( t ),rgrg ) ); 
     78                         
     79                        //Merger 
     80                } 
     81                L.logit ( Eth1_log, P1._epdf().mean() ); 
     82                L.logit ( Eth2_log, P2._epdf().mean() ); 
     83                L.logit ( Ethg_log, PG._epdf().mean() ); 
     84                L.logit ( Data_log, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) ); 
     85                L.logit ( LL_log, vec_3 ( P1._ll(), P2._ll(), PG._ll() ) ); 
     86                L.step (  ); 
    7787        } 
    78         L.step(true); 
    79          
     88        L.finalize( ); 
     89        L.itsave ( "merg.it" ); 
    8090}