root/library/tests/testSmp.cpp @ 504

#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;
        shared_ptr<enorm<ldmat> > eN = new enorm<ldmat>();
        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;
        shared_ptr<egamma> eG = new egamma();
        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;
        shared_ptr<mgamma> mG = new mgamma();
        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;
        shared_ptr<emix> eMix = new emix();
        Array<shared_ptr<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 );

        // emix::set_parameters requires the first mpdf to be named
        mG->set_rv(x);
        mG->set_rvc(y);
        mComs ( 0 ) = mG;

        eN->set_mu ( vec_2 ( 0.0, 0.0 ) );
        shared_ptr<mepdf> mEnorm = new mepdf ( 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.get() );
        mepdf meG ( eG.get() );
        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;

}