root/library/tests/emix_test.cpp @ 477

#include "shared_ptr.h"
#include "stat/exp_family.h"
#include "stat/emix.h"

using namespace bdm;

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

double normcoef ( const epdf* ep, const vec &xb, const vec &yb, int Ngr = 100 ) {
        mat PPdf ( Ngr + 1, Ngr + 1 );
        vec rgr ( 2 );

        int i = 0, j = 0;
        double xstep = ( xb ( 1 ) - xb ( 0 ) ) / Ngr;
        double ystep = ( yb ( 1 ) - yb ( 0 ) ) / Ngr;

        for ( double x = xb ( 0 ); x <= xb ( 1 ); x += xstep, i++ ) {
                rgr ( 0 ) = x;
                j = 0;
                for ( double y = yb ( 0 ); y <= yb ( 1 ); y += ystep, j++ ) {
                        rgr ( 1 ) = y;
                        PPdf ( i, j ) = exp ( ep->evallog ( rgr ) );
                }
        }
        return sumsum ( PPdf ) *xstep*ystep;
}

int main() {
        RNG_randomize();

        RV x ( "{x }" );
        RV y ( "{y }" );
        RV xy = concat ( x, y );

        enorm<ldmat> E1;
        E1.set_rv ( xy );

        E1.set_parameters ( "1.00054 1.0455" , mat ( "0.740142 -0.259015; -0.259015 1.0302" ) );
        enorm<ldmat> E2;
        E2.set_rv ( xy );
        E2.set_parameters ( "-1.2 -0.1" , mat ( "1 0.4; 0.4 0.5" ) );

        Array<epdf*> A1 ( 1 );
        A1 ( 0 ) = &E1;

        emix M1;
        M1.set_rv ( xy );
        M1.set_parameters ( vec ( "1" ), A1, false );
        cout << "======== test if ARX and emix with one ARX are the same ===" << endl;

        epdf* Mm = M1.marginal ( y );
        epdf* Am = E1.marginal ( y );
        mpdf* Mc = M1.condition ( y );
        mpdf* Ac = E1.condition ( y );

        cout << * ( ( mlnorm<ldmat>* ) Ac ) << endl;

        cout << "Mix marg at 0: " << Mm->evallog ( vec_1 ( 0.0 ) ) << endl;
        cout << "ARX marg at 0: " << Am->evallog ( vec_1 ( 0.0 ) ) << endl;
        cout << "Mix cond at 0,0: " << Mc->evallogcond ( vec_1 ( 0.0 ), vec_1 ( 0.0 ) ) << endl;
        cout << "ARX cond at 0,0: " << Ac->evallogcond ( vec_1 ( 0.0 ), vec_1 ( 0.0 ) ) << endl;

        cout << "======== Mixture with two components ===" << endl;
        Array<epdf*> A2 ( 2 );
        A2 ( 0 ) = &E1;
        A2 ( 1 ) = &E2;

        emix M2;
        M2.set_rv ( xy );
        M2.set_parameters ( vec ( "1" ), A2, false );


        cout << "Mixture normalization: " << normcoef ( &M2, vec ( "-3 3 " ), vec ( "-3 3 " ) ) << endl;

        int N = 3;
        vec ll ( N );
        vec ll2 ( N );
        mat Smp = M2.sample_m ( N );
        ll = M2.evallog_m ( Smp );

        vec Emu = sum ( Smp, 2 ) / N;
        mat Er = ( Smp * Smp.transpose() ) / N - outer_product ( Emu, Emu );

        cout << "original empirical mean: " << Emu << endl;
        cout << "original empirical variance: " << Er << endl;

        shared_ptr<epdf> Mg ( dynamic_cast<epdf *> ( M2.marginal ( y ) ) );
        mpdf* Cn = ( mpdf* ) M2.condition ( x );

        it_file it ( "emix_test.it" );
        it << Name ( "Smp" ) << Smp;

        cout << "marginal mean: " << Mg->mean() << endl;

        cout << "========== putting them back together ======= " << endl;
        mepdf mMg ( Mg );
        Array<mpdf*> AA ( 2 );
        AA ( 0 ) = Cn;
        AA ( 1 ) = &mMg;
        mprod mEp ( AA );

        for ( int j = 0; j < N; j++ ) {
                ll2 ( j ) = mEp.evallogcond ( Smp.get_col ( j ), vec ( 0 ) );
        }
        it << Name ( "ll" ) << ll;
        it << Name ( "ll2" ) << ll2;


}