root/tests/egiw_test.cpp @ 262

Revision 254, 2.2 kB (checked in by smidl, 16 years ago)

create namespace bdm

Line 
1#include <stat/libEF.h>
2using namespace bdm;
3
4//These lines are needed for use of cout and endl
5using std::cout;
6using std::endl;
7
8void Test ( const egiw &E ) {
9}
10
11int main() {
12        cout << "Testing eGiw(1,1)"<<endl;
13        {
14                // Setup model
15                double mu=1.1;
16                double s=0.1;
17
18                // TEST 1x1 EGIW
19                mat V ( 2,2 );
20                V ( 0,0 ) = pow ( mu,2 ) +s;
21                V ( 1,0 ) = mu;
22                V ( 0,1 ) = V ( 1,0 );
23                V ( 1,1 ) = 1.0;
24
25                double nu=10;
26
27                RV thr ( "{thr }", "2" );
28                egiw E ( thr,nu*V,nu );
29                cout << "egiw mean value:" << E.mean() <<endl;
30                cout << "egiw normalizing constant:" << E.lognc() <<endl;
31
32                int n=100;
33                vec t_val ( 2 );
34
35                mat pdf ( 2*n,n );
36                vec Mu ( 2*n );
37                vec Si ( n );
38
39                for ( int i=0;i<2*n;i++ ) {
40                        Mu ( i ) = -2+i* ( 1.0/ ( n ) ) *3.0;
41                        t_val ( 0 ) = Mu ( i );
42                        for ( int j=0;j<n;j++ ) {
43                                Si ( j ) = ( j+1 ) * ( 1.0/n ) *2;
44                                t_val ( 1 ) = Si ( j );
45
46                                pdf ( i,j ) =E.evallog ( t_val );
47                        }
48                }
49
50                mat Pdf=exp ( pdf );
51                vec fm=sum ( Pdf,2 );
52                vec fs=sum ( Pdf,1 );
53                cout << "Numerical mean: " << vec_2 ( Mu*fm/sum ( fm ), Si*fs/sum ( fs ) ) <<endl;
54                cout << "Numerical integral of pdf: "<<sumsum ( Pdf/n/n*3*2 ) <<endl;
55        }
56        cout << "Testing Egiw(1,2)"<<endl;
57        {
58                // Setup model
59                double mu=1.1; //unit step parametr
60                double b=3.0; // sequence of <1 -1 1 -1...>
61                double s=0.1;
62
63
64                // TEST 1x1 EGIW
65                mat V ( 3,3 );
66                V ( 0,0 ) = pow ( mu,2 ) +pow ( b ,2 )  +s;
67                V ( 1,0 ) = mu;
68                V ( 2,0 ) = b;
69
70                V ( 0,1 ) = V ( 1,0 );
71                V ( 1,1 ) = 1.0;
72                V ( 2,1 ) = 0.0;
73
74                V ( 0,2 ) = V ( 2,0 );
75                V ( 1,2 ) = V ( 2,1 );
76                V ( 2,2 ) = 1.0;
77
78
79                double nu=20;
80
81                RV thr ( "{thr }", "3" );
82                egiw E ( thr,nu*V,nu );
83                cout << "egiw mean value:" << E.mean() <<endl;
84                cout << "egiw normalizing constant:" << E.lognc() <<endl;
85
86                int n=100;
87                vec t_val ( 3 );
88
89                mat Tmp= zeros ( 2*n,n );
90
91                double summ=0.0;
92                for ( int k=0;k<n;k++ ) { // ALL b
93                        t_val ( 1 ) = 1 + k* ( 1.0/n ) * 4.0;
94                        for ( int i=0;i<2*n;i++ ) { //ALL mu
95                                t_val ( 0 ) = -2+i* ( 1.0/  n ) *3.0;
96                                for ( int j=0;j<n;j++ ) { // All sigma
97                                        t_val ( 2 ) = ( j+1 ) * ( 1.0/n ) *2.0;
98
99                                        Tmp ( i,j ) = E.evallog ( t_val );
100                                }
101                        }
102                        summ += sumsum ( exp ( Tmp ) ) /n/n/n*3.0*2.0*4.0;
103                }
104
105
106                cout << "Numerical integral of pdf: "<<summ <<endl;
107        }
108       
109       
110}
Note: See TracBrowser for help on using the browser.