root/tests/enorm_test.cpp @ 179

Revision 179, 1.9 kB (checked in by smidl, 16 years ago)

Tests of the new methods for emorm and mlnorm

RevLine 
[179]1#include <stat/libEF.h>
2#include <stat/emix.h>
3using namespace itpp;
4
5//These lines are needed for use of cout and endl
6using std::cout;
7using std::endl;
8
9double normcoef ( const epdf* ep,const vec &xb, const vec &yb, int Ngr=100 ) {
10        mat PPdf ( Ngr+1,Ngr+1 );
11        vec rgr ( 2 );
12
13        int i=0,j=0;
14        double xstep= ( xb ( 1 )-xb ( 0 ) ) /Ngr;
15        double ystep= ( yb ( 1 )-yb ( 0 ) ) /Ngr;
16
17        for ( double x=xb ( 0 );x<=xb ( 1 );x+= xstep,i++ ) {
18                rgr ( 0 ) =x;j=0;
19                for ( double y=yb ( 0 );y<=yb ( 1 );y+=ystep,j++ ) {
20                        rgr ( 1 ) =y;
21                        PPdf ( i,j ) =exp ( ep->evalpdflog ( rgr ) );
22                }
23        }
24        return sumsum ( PPdf ) *xstep*ystep;
25}
26
27int main() {
28        RNG_randomize();
29
30        cout << "Testing enorm(2,2)"<<endl;
31        // Setup model
32        vec mu ( "1.1 -1" );
33        ldmat R ( mat ( "1 -0.5; -0.5 2" ) );
34
35        RV x ( "{x }" );
36        RV y ( "{y }" );
37
38        enorm<ldmat> E ( concat ( x,y ) );
39        E.set_parameters ( mu,R );
40        cout << "enorm mean value:" << E.mean() <<endl;
41        cout << "enorm mean variance:" << R.to_mat() <<endl;
42        cout << "enorm normalizing constant:" << E.lognc() <<endl;
43        cout << " Numerically integrates to [1.0]: " << normcoef ( &E, vec ( "-5 5" ), vec ( "-5 5" ) ) << endl;
44
45        int N=1000;
46        vec ll ( N );
47        mat Smp=E.sample ( 1000 );
48        vec Emu = sum ( Smp,2 ) /N;
49        mat Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );
50
51        cout << "original empirical mean: " <<Emu <<endl;
52        cout << "original empirical variance: " <<Er <<endl;
53
54        enorm<ldmat>* Mg = E.marginal ( y );
55        mlnorm<ldmat>* Cn = E.condition ( x );
56
57
58        cout<< "marginal mean: " << Mg->mean() <<endl;
59
60        cout << "========== putting them back together ======= "<<endl;
61        mepdf mMg ( *Mg );
62        Array<mpdf*> A ( 2 );
63        A ( 0 ) = Cn;
64        A ( 1 ) = &mMg;
65        mprod mEp ( A );
66
67        Smp=mEp.samplecond ( vec ( 0 ),ll,1000 );
68        Emu = sum ( Smp,2 ) /N;
69        Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );
70
71        cout << "composite mean: " <<Emu <<endl;
72        cout << "composite variance: " <<Er <<endl;
73}
Note: See TracBrowser for help on using the browser.