root/tests/egiw_test.cpp @ 305

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

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

void Test ( const egiw &E ) {
}

int main() {
        cout << "Testing eGiw(1,1)"<<endl;
        {
                // Setup model
                double mu=1.1;
                double s=0.1;

                // TEST 1x1 EGIW
                mat V ( 2,2 );
                V ( 0,0 ) = pow ( mu,2 ) +s;
                V ( 1,0 ) = mu;
                V ( 0,1 ) = V ( 1,0 );
                V ( 1,1 ) = 1.0;

                double nu=10;

                egiw E ( 1,nu*V,nu );
                cout << "egiw mean value:" << E.mean() <<endl;
                cout << "egiw normalizing constant:" << E.lognc() <<endl;

                int n=100;
                vec t_val ( 2 );

                mat pdf ( 2*n,n );
                vec Mu ( 2*n );
                vec Si ( n );

                for ( int i=0;i<2*n;i++ ) {
                        Mu ( i ) = -2+i* ( 1.0/ ( n ) ) *3.0;
                        t_val ( 0 ) = Mu ( i );
                        for ( int j=0;j<n;j++ ) {
                                Si ( j ) = ( j+1 ) * ( 1.0/n ) *2;
                                t_val ( 1 ) = Si ( j );

                                pdf ( i,j ) =E.evallog ( t_val );
                        }
                }

                mat Pdf=exp ( pdf );
                vec fm=sum ( Pdf,2 );
                vec fs=sum ( Pdf,1 );
                cout << "Numerical mean: " << vec_2 ( Mu*fm/sum ( fm ), Si*fs/sum ( fs ) ) <<endl;
                cout << "Numerical integral of pdf: "<<sumsum ( Pdf/n/n*3*2 ) <<endl;
        }
        cout << "Testing Egiw(1,2)"<<endl;
        {
                // Setup model
                double mu=1.1; //unit step parametr
                double b=3.0; // sequence of <1 -1 1 -1...>
                double s=0.1;


                // TEST 1x1 EGIW
                mat V ( 3,3 );
                V ( 0,0 ) = pow ( mu,2 ) +pow ( b ,2 )  +s;
                V ( 1,0 ) = mu;
                V ( 2,0 ) = b;

                V ( 0,1 ) = V ( 1,0 );
                V ( 1,1 ) = 1.0;
                V ( 2,1 ) = 0.0;

                V ( 0,2 ) = V ( 2,0 );
                V ( 1,2 ) = V ( 2,1 );
                V ( 2,2 ) = 1.0;


                double nu=20;

                egiw E ( 1,nu*V,nu );
                cout << "egiw mean value:" << E.mean() <<endl;
                cout << "egiw normalizing constant:" << E.lognc() <<endl;

                int n=100;
                vec t_val ( 3 );

                mat Tmp= zeros ( 2*n,n );

                double summ=0.0;
                for ( int k=0;k<n;k++ ) { // ALL b
                        t_val ( 1 ) = 1 + k* ( 1.0/n ) * 4.0;
                        for ( int i=0;i<2*n;i++ ) { //ALL mu
                                t_val ( 0 ) = -2+i* ( 1.0/  n ) *3.0;
                                for ( int j=0;j<n;j++ ) { // All sigma
                                        t_val ( 2 ) = ( j+1 ) * ( 1.0/n ) *2.0;

                                        Tmp ( i,j ) = E.evallog ( t_val );
                                }
                        }
                        summ += sumsum ( exp ( Tmp ) ) /n/n/n*3.0*2.0*4.0;
                }


                cout << "Numerical integral of pdf: "<<summ <<endl;
        }
        
        
}