root/applications/mpdm/TR2244/merger_iter_cond.cpp @ 379

#include <stat/libEF.h>
#include <estim/merger.h>

using namespace bdm;

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

int main() {

        RNG_randomize();

        RV y ( "{y }","1" );
        RV u1 ( "{u1 }","1" );
//      RV u2 ( "{u2 }","1" );

        RV all = y;
        all.add ( u1 );
//      all.add ( u2 );

        mlnorm<fsqmat> f1; f1.set_rv( y ); f1.set_rvc( u1 );
        enorm<fsqmat>  f2; f2.set_rv( y );// f2.set_rvc( u2 );

        //Differneces in constant term are essential
        f1.set_parameters ( "1","1",mat ( "0.01" ) );
        f2.set_parameters ( "2", mat ( "0.02" ) );

        mepdf mf2(&f2);
        
        Array<mpdf* > A ( 2 );
        A ( 0 ) =&f1;
        A ( 1 ) =&mf2;

        merger_mix M ( A );
        M.set_debug_file("iter_cond_debug.it");
        enorm<ldmat> g0; g0.set_rv ( all );
        mat Cov=10*eye ( 2 );
/*      Cov(1,2)=0.29;
        Cov(2,1)=0.29;*/
        g0.set_parameters ( vec ( "2 1" ),Cov);// +1*ones ( 3,3 ) );

//      enorm<ldmat>* teste=g0.marginal(concat(u1,u2));
//      mlnorm<ldmat>* testm=(mlnorm<ldmat>*)teste->condition(u2);
        
        M.set_parameters (1,0.5);
        M.set_method(LOGNORMAL, 1.2);
        
        Array<vec> YUU(2); 
        YUU(0)=linspace(1.5,2.5,20);
        YUU(1)=linspace(0.5,1.5,20);
        
        YUU(0)=linspace(0,4,20);
        YUU(1)=linspace(-1,3,20);
//YUU(2)=linspace(-3,3,20);

        //M.set_grid(YUU);
        
        int Ntrials=1;
        vec A1s ( Ntrials );
//      vec A2s ( Ntrials );
        vec R1s ( Ntrials );
        vec R2s ( Ntrials );
        vec C1s ( Ntrials );
        vec C2s ( Ntrials );

        for ( int it=0;it<Ntrials;it++ ) {
                M.set_support( g0,1000);
                M.merge();
                
                MixEF &MM = M._Mix();
                epdf* MP = MM._Coms ( 0 )->epredictor ( );
                MP->set_rv(all);
                
                //Missing piece
                enorm<ldmat>* miss=     (enorm<ldmat>*) MP->marginal(u1);
                cout << "Missing: " << miss->_mu() << " +- " << miss->_R().to_mat() << endl;

                RV yu1 = y; yu1.add ( u1 );
//              RV yu2 = y; yu2.add ( u2 );
                enorm<ldmat>* P1m= ( enorm<ldmat>* ) MP->marginal ( yu1 );
                enorm<ldmat>* P2m= ( enorm<ldmat>* ) MP->marginal ( y );
                mlnorm<ldmat>* P1c= ( mlnorm<ldmat>* ) ( P1m->condition ( y ) );
//              mlnorm<ldmat>* P2c= ( mlnorm<ldmat>* ) ( P2m->condition ( y ) );

                A1s(it) = P1c->_A()(0,0);
//              A2s(it) = P2c->_A()(0,0);
                R1s(it) = P1c->_R()(0,0);
                R2s(it) = P2m->_R().to_mat()(0,0);
                C1s(it) = P1c->_mu_const()(0);
                C2s(it) = P2m->_mu()(0);

                cout << "mean: " << MM._Coms(0)->_e()->mean() <<endl;
        }
        
        
        
        cout << "A1s:" <<A1s<<endl;
        cout << "C2s:" <<C2s<<endl;
        double A1mean = sum(A1s)/Ntrials;
//      double A2mean = sum(A2s)/Ntrials;
        double C1mean = sum(C1s)/Ntrials;
        double C2mean = sum(C2s)/Ntrials;
        double R1mean = sum(R1s)/Ntrials;
        double R2mean = sum(R2s)/Ntrials;
        cout << "A1: " << A1mean << " +- " << 2*sqrt(sum_sqr(A1s)/Ntrials-A1mean*A1mean) <<endl;
//      cout << "A2: " << A2mean << " +- " << 2*sqrt(sum_sqr(A2s)/Ntrials-A2mean*A2mean) <<endl;
        cout << "C1: " << C1mean << " +- " << 2*sqrt(sum_sqr(C1s)/Ntrials-C1mean*C1mean) <<endl;
        cout << "C2: " << C2mean << " +- " << 2*sqrt(sum_sqr(C2s)/Ntrials-C2mean*C2mean) <<endl;
        cout << "R1: " << R1mean << " +- " << 2*sqrt(sum_sqr(R1s)/Ntrials-R1mean*R1mean) <<endl;
        cout << "R2: " << R2mean << " +- " << 2*sqrt(sum_sqr(R2s)/Ntrials-R2mean*R2mean) <<endl;
        
        
}