[188] | 1 | #include <stat/libEF.h> |
---|
[254] | 2 | using namespace bdm; |
---|
[188] | 3 | |
---|
| 4 | //These lines are needed for use of cout and endl |
---|
| 5 | using std::cout; |
---|
| 6 | using std::endl; |
---|
| 7 | |
---|
| 8 | void Test ( const egiw &E ) { |
---|
| 9 | } |
---|
| 10 | |
---|
| 11 | int 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 | |
---|
[211] | 46 | pdf ( i,j ) =E.evallog ( t_val ); |
---|
[188] | 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 | |
---|
[211] | 99 | Tmp ( i,j ) = E.evallog ( t_val ); |
---|
[188] | 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 | } |
---|