root/tests/enorm_test.cpp @ 179

#include <stat/libEF.h>
#include <stat/emix.h>
using namespace itpp;

//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->evalpdflog ( rgr ) );
                }
        }
        return sumsum ( PPdf ) *xstep*ystep;
}

int main() {
        RNG_randomize();

        cout << "Testing enorm(2,2)"<<endl;
        // Setup model
        vec mu ( "1.1 -1" );
        ldmat R ( mat ( "1 -0.5; -0.5 2" ) );

        RV x ( "{x }" );
        RV y ( "{y }" );

        enorm<ldmat> E ( concat ( x,y ) );
        E.set_parameters ( mu,R );
        cout << "enorm mean value:" << E.mean() <<endl;
        cout << "enorm mean variance:" << R.to_mat() <<endl;
        cout << "enorm normalizing constant:" << E.lognc() <<endl;
        cout << " Numerically integrates to [1.0]: " << normcoef ( &E, vec ( "-5 5" ), vec ( "-5 5" ) ) << endl;

        int N=1000;
        vec ll ( N );
        mat Smp=E.sample ( 1000 );
        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;

        enorm<ldmat>* Mg = E.marginal ( y );
        mlnorm<ldmat>* Cn = E.condition ( x );


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

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

        Smp=mEp.samplecond ( vec ( 0 ),ll,1000 );
        Emu = sum ( Smp,2 ) /N;
        Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );

        cout << "composite mean: " <<Emu <<endl;
        cout << "composite variance: " <<Er <<endl;
}