root/library/tests/testSmp.cpp @ 488

#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;

void disp ( const vec &tmu, const mat &tR, const mat &Smp ) {
        int N = Smp.cols();
        vec Emu = Smp * ones ( N ) / N ;
        mat Er = ( Smp * Smp.transpose() ) / N - outer_product ( Emu, Emu );
        cout << "True mu:" << tmu << endl;
        cout << "Emp  mu:" << Emu << endl;

        cout << "True R:" << tR << endl;
        cout << "Emp  R:" << Er << endl;
}

int main() {

        RNG_randomize();

        RV x ( "{x }", "2" );
        RV y ( "{y }", "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 );

        cout << "====== ENorm ====== " << endl;
        enorm<ldmat> eN;
        eN.set_parameters ( mu0, R );
        mat Smp = eN.sample_m ( N );

        disp ( mu0, R.to_mat(), Smp );

        cout << "====== MlNorm ====== " << endl;
        mat I = eye ( 2 );
        mlnorm<ldmat> ML;
        ML.set_parameters ( I, zeros ( 2 ), R );
        Smp = ML.samplecond_m ( mu0, N );

        disp ( mu0, R.to_mat(), Smp );

        cout << "====== EGamma ====== " << endl;
        vec a = "100000,10000";
        vec b = a / 10.0;
        egamma eG;
        eG.set_parameters ( a, b );

        cout << eG.evallog ( a ) << endl;
        Smp = eG.sample_m ( N );

        vec g_mu = elem_div ( a, b );
        vec g_var = elem_div ( a, pow ( b, 2.0 ) );
        disp ( g_mu, diag ( g_var ), Smp );

        cout << "====== MGamma ====== " << endl;
        mgamma mG;
        double k = 10.0;
        mG.set_parameters ( k, mu0 );

        Smp = mG.samplecond_m ( mu0, N );
        disp ( mu0, pow ( mu0, 2.0 ) / k, Smp );

        cout << "======= EMix ======== " << endl;
        emix eMix;
        Array<epdf*> Coms ( 2 );
        Coms ( 0 ) = &eG;
        Coms ( 1 ) = &eN;

        eMix.set_parameters ( vec_2 ( 0.5, 0.5 ), Coms );
        vec smp = eMix.sample();
        Smp = eMix.sample_m ( N );
        disp ( eMix.mean(), zeros ( 2 ), Smp );

        cout << "======= MEpdf ======== " << endl;
        mepdf meMix ( &eMix );

        Smp = meMix.samplecond_m ( mu0, N );
        disp ( eMix.mean(), zeros ( 2 ), Smp );

        cout << "======= MMix ======== " << endl;
        mmix mMix;
        Array<shared_ptr<mpdf> > mComs ( 2 );
        mComs ( 0 ) = &mG;
        eN.set_mu ( vec_2 ( 0.0, 0.0 ) );
        mepdf mEnorm ( &eN );
        mComs ( 1 ) = &mEnorm;
        mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs );

        Smp = mMix.samplecond_m ( mu0, N );
        disp ( 0.5*eN.mean()+0.4*eG.mean(), zeros ( 2 ), Smp );

        cout << "======= EProd ======== " << endl;
        // we have to change eG.rv to y
        eN.set_rv ( x );
        eG.set_rv ( y );
        //create array
        Array<mpdf*> A ( 2 );
        mepdf meN ( &eN );
        mepdf meG ( &eG );
        A ( 0 ) = &meN;
        A ( 1 ) = &meG;

        mprod eP ( A );
        mat epV = zeros ( 4, 4 );
        epV.set_submatrix ( 0, 0, V0 );
        epV.set_submatrix ( 2, 2, diag ( g_var ) );

        vec v0 = vec ( 0 );
        Smp = eP.samplecond ( v0, N );
        disp ( concat ( eN.mean(), eG.mean() ), epV, Smp );

        cout << "======= eWishart ======== " << endl;
        mat wM = "1.0 0.9; 0.9 1.0";
        eWishartCh eW;
        eW.set_parameters ( wM / 100, 100 );
        mat mea = zeros ( 2, 2 );
        mat Ch ( 2, 2 );
        for ( int i = 0; i < 100; i++ ) {
                Ch = eW.sample_mat();
                mea += Ch.T() * Ch;
        }
        cout << mea / 100 << endl;

        cout << "======= rwiWishart ======== " << endl;
        rwiWishartCh rwW;
        rwW.set_parameters ( 2, 0.1, "1 1", 0.9 );
        mea = zeros ( 2, 2 );
        mat wMch = chol ( wM );
        for ( int i = 0; i < 100; i++ ) {
                vec tmp = rwW.samplecond ( vec ( wMch._data(), 4 ) );
                copy_vector ( 4, tmp._data(), Ch._data() );
                mea += Ch.T() * Ch;
        }
        cout << mea / 100 << endl;
        //Exit program:
        return 0;

}