root/library/tests/emix_test.cpp @ 522

#include "shared_ptr.h"
#include "stat/exp_family.h"
#include "stat/emix.h"
#include "mat_checks.h"
#include "UnitTest++.h"
#include "test_util.h"

const double epsilon = 0.00001;

using namespace bdm;

static void check_mean ( emix &distrib_obj, int nsamples, const vec &mean, double tolerance );

static void check_covariance ( emix &distrib_obj, int nsamples, const mat &R, double tolerance);

TEST ( test_emix_1 ) {
        RV x ( "{emixx }" );
        RV y ( "{emixy }" );
        RV xy = concat ( x, y );
        vec mu0 ( "1.00054 1.0455" );

        shared_ptr<enorm<ldmat> > E1 = new enorm<ldmat>();
        E1->set_rv ( xy );
        E1->set_parameters ( mu0 , mat ( "0.740142 -0.259015; -0.259015 1.0302" ) );

        shared_ptr<enorm<ldmat> > E2 = new enorm<ldmat>();
        E2->set_rv ( xy );
        E2->set_parameters ( "-1.2 -0.1" , mat ( "1 0.4; 0.4 0.5" ) );

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

        emix M1;
        M1.set_rv ( xy );
        M1.set_parameters ( vec ( "1" ), A1 );

        // test if ARX and emix with one ARX are the same
        shared_ptr<epdf> Mm = M1.marginal ( y );
        shared_ptr<epdf> Am = E1->marginal ( y );
        shared_ptr<mpdf> Mc = M1.condition ( y );
        shared_ptr<mpdf> Ac = E1->condition ( y );

        mlnorm<ldmat> *wacnd = dynamic_cast<mlnorm<ldmat> *>( Ac.get() );
        CHECK(wacnd);
        if ( wacnd ) {
                CHECK_CLOSE ( mat ( "-0.349953" ), wacnd->_A(), epsilon );
                CHECK_CLOSE ( vec ( "1.39564" ), wacnd->_mu_const(), epsilon );
                CHECK_CLOSE ( mat ( "0.939557" ), wacnd->_R(), epsilon );
        }

        double same = -1.46433;
        CHECK_CLOSE ( same, Mm->evallog ( vec_1 ( 0.0 ) ), epsilon );
        CHECK_CLOSE ( same, Am->evallog ( vec_1 ( 0.0 ) ), epsilon );
        CHECK_CLOSE ( 0.145974, Mc->evallogcond ( vec_1 ( 0.0 ), vec_1 ( 0.0 ) ), epsilon );
        CHECK_CLOSE ( -1.92433, Ac->evallogcond ( vec_1 ( 0.0 ), vec_1 ( 0.0 ) ), epsilon );

        // mixture with two components
        Array<shared_ptr<epdf> > A2 ( 2 );
        A2 ( 0 ) = E1;
        A2 ( 1 ) = E2;

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


        // mixture normalization
        CHECK_CLOSE ( 1.0, normcoef ( &M2, vec ( "-3 3 " ), vec ( "-3 3 " ) ), 0.1 );

        int N = 3;
        mat Smp = M2.sample_m ( N );

        vec exp_ll ( "-5.0 -2.53563 -2.62171" );
        vec ll = M2.evallog_m ( Smp );
        CHECK_CLOSE ( exp_ll, ll, 5.0 );

        check_mean ( M2, N, mu0, 1.0 );

        mat observedR ( "0.740142 -0.259015; -0.259015 1.0302" );
        check_covariance ( M2, N, observedR, 2.0);

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

        // marginal mean
        CHECK_CLOSE ( vec ( "1.0" ), Mg->mean(), 0.1 );
}

TEST ( test_emix_2 ) {
        int N = 10000; //number of samples
        vec mu0 ( "1.5 1.7" );
        mat V0 ( "1.2 0.3; 0.3 5" );
        ldmat R = ldmat ( V0 );

        shared_ptr<enorm<ldmat> > eN = new enorm<ldmat>();
        eN->set_parameters ( mu0, R );

        vec a = "100000,10000";
        vec b = a / 10.0;
        shared_ptr<egamma> eG = new egamma();
        eG->set_parameters ( a, b );

        emix eMix;
        Array<shared_ptr<epdf> > Coms ( 2 );
        Coms ( 0 ) = eG;
        Coms ( 1 ) = eN;

        eMix.set_parameters ( vec_2 ( 0.5, 0.5 ), Coms );
        check_mean ( eMix, N, eMix.mean(), 0.1 );
}

static void check_mean ( emix &distrib_obj, int nsamples, const vec &mean, double tolerance ) {
        int tc = 0;
        Array<vec> actual(CurrentContext::max_trial_count);
        do {
                mat smp = distrib_obj.sample_m ( nsamples );
                vec emu = sum ( smp, 2 ) / nsamples;
                actual( tc ) = emu;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__);
                stream << "Expected " << mean << " +/- " << tolerance << " but was " << actual;

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
        }
}

static void check_covariance ( emix &distrib_obj, int nsamples, const mat &R, double tolerance) {
        int tc = 0;
        Array<mat> actual(CurrentContext::max_trial_count);
        do {
                mat smp = distrib_obj.sample_m ( nsamples );
                vec emu = sum ( smp, 2 ) / nsamples;
                mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
                actual( tc ) = er;
                ++tc;
        } while ( ( tc < CurrentContext::max_trial_count ) &&
                  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
        if ( ( tc == CurrentContext::max_trial_count ) &&
             ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
                UnitTest::MemoryOutStream stream;
                UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__);
                stream << "Expected " << R << " +/- " << tolerance << " but was " << actual;

                UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}