root/mpdm/merg_pred.cpp @ 162

#include <estim/arx.h>
#include <stat/libEF.h>
#include <stat/loggers.h>
#include <stat/merger.h>
using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main() {
        // Setup model
        RV y ( "{y }" );
        RV u1 ( "{u1 }" );
        RV u2 ( "{u2 }" );

        // Full system
        vec thg ( "0.7 1 1" ); //Simulated system
        //y=a y_t-1 + u1 + u2
        double sqr=0.1;
        int ord = 1;

        // Estimated systems ARX(2)
        RV thri ( "{thr_i }",vec_1 ( 2+1 ) );
        RV thrg ( "{thr_g }",vec_1 ( 3+1 ) );
        // Setup values

        //ARX constructor
        mat V0 = 0.001*eye ( thri.count() ); V0 ( 0,0 ) *= 10; //
        mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; //
        double nu0 = ord+1;
        double frg = 0.99;

        ARX P1 ( thri, V0, nu0, frg );
        ARX P2 ( thri, V0, nu0, frg );
        ARX PG ( thrg, V0g, nu0, frg );
        //Test estimation
        int ndat = 10000;
        int t;

        // Logging
        dirfilelog L ( "exp/merg",ndat );

        int Eth1_log = L.add ( thri,"P1" );
        int Eth2_log = L.add ( thri,"P2" );
        int Ethg_log = L.add ( thrg,"PG" );
        int Data_log = L.add ( RV ( "{Y U1 U2 }" ), "" );
        int LL_log   = L.add ( RV ( "{1 2 G }" ), "LL" );

        L.init();

        vec Yt ( ndat );

        Yt.set_subvector ( 0,randn ( ord ) ); //initial values
        vec rgrg ( thg.length() );
        vec rgri ( thri.count() );

        for ( t=0; t<ndat; t++ ) {
                // True system
                if ( t>0 ) {
                        rgrg ( 0 ) =Yt ( t-1 );
                        rgrg ( 1 ) = pow(sin ( ( t/40.0 ) *pi ),3);
                        rgrg ( 2 ) = pow(cos ( ( t/40.0 ) *pi ),3);

                        Yt ( t ) = thg*rgrg + sqr * NorRNG();

                        // 1st
                        rgri ( 0 ) =Yt ( t );
                        rgri ( 1 ) =Yt ( t-1 );
                        rgri ( 2 ) =rgrg ( 1 );
                        P1.bayes ( rgri );
                        // 2nd
                        rgri ( 2 ) =rgrg ( 2 );
                        P2.bayes ( rgri );

                        //Global
                        PG.bayes ( concat ( Yt ( t ),rgrg ) );
                        
                        //Merger
                }
                L.logit ( Eth1_log, P1._epdf().mean() );
                L.logit ( Eth2_log, P2._epdf().mean() );
                L.logit ( Ethg_log, PG._epdf().mean() );
                L.logit ( Data_log, vec_3 ( Yt ( t ), rgrg ( 1 ), rgrg ( 2 ) ) );
                L.logit ( LL_log, vec_3 ( P1._ll(), P2._ll(), PG._ll() ) );
                L.step (  );
        }
        L.finalize( );
        L.itsave ( "merg.it" );
}