root/library/tests/enorm_test.cpp @ 394

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

        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;
        E.set_rv(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;

        epdf* Mg = E.marginal ( y );
        mpdf* 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 ),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;
        
        cout << endl << " test of pdflog at zero "<<endl;
        cout << "original: " << E.evallog(vec("0 0")) <<endl;
        cout << "composite: " << mEp.evallogcond(vec("0 0"),vec(0)) << endl;
}