[386] | 1 | #include "estim/mixtures.h" |
---|
| 2 | #include "estim/arx.h" |
---|
| 3 | #include "stat/exp_family.h" |
---|
[254] | 4 | using namespace bdm; |
---|
[176] | 5 | |
---|
| 6 | //These lines are needed for use of cout and endl |
---|
| 7 | using std::cout; |
---|
| 8 | using std::endl; |
---|
| 9 | |
---|
| 10 | void disp_mix2D ( const MixEF &Mi ) { |
---|
[477] | 11 | const eprod &ep = ( ( const eprod& ) ( Mi.posterior() ) ); |
---|
[176] | 12 | int i; |
---|
[477] | 13 | mat Var ( 2, 2 ), M ( 1, 2 ); |
---|
| 14 | for ( i = 0; i < ep._rv().length() - 1; i++ ) { // -1 => last one is the weight |
---|
| 15 | ( ( egiw* ) ep ( i ) )->mean_mat ( M, Var ); |
---|
| 16 | cout << "Mean: " << M << endl << "Variance: " << endl << Var << endl; |
---|
[176] | 17 | } |
---|
[477] | 18 | cout << "Weights: " << ep ( i )->mean() << endl; |
---|
[176] | 19 | } |
---|
| 20 | |
---|
[477] | 21 | mat grid2D ( const MixEF &M, const vec &xb, const vec &yb, int Ngr = 100 ) { |
---|
| 22 | mat PPdf ( Ngr + 1, Ngr + 1 ); |
---|
[176] | 23 | vec rgr ( 3 ); |
---|
| 24 | |
---|
[477] | 25 | rgr ( 2 ) = 1; |
---|
| 26 | int i = 0, j = 0; |
---|
| 27 | double xstep = ( xb ( 1 ) - xb ( 0 ) ) / Ngr; |
---|
| 28 | double ystep = ( yb ( 1 ) - yb ( 0 ) ) / Ngr; |
---|
| 29 | |
---|
| 30 | for ( double x = xb ( 0 ); x <= xb ( 1 ); x += xstep, i++ ) { |
---|
| 31 | rgr ( 0 ) = x; |
---|
| 32 | j = 0; |
---|
| 33 | for ( double y = yb ( 0 ); y <= yb ( 1 ); y += ystep, j++ ) { |
---|
| 34 | rgr ( 1 ) = y; |
---|
| 35 | PPdf ( i, j ) = exp ( M.logpred ( rgr ) ); |
---|
[176] | 36 | } |
---|
| 37 | } |
---|
| 38 | return PPdf; |
---|
| 39 | } |
---|
| 40 | |
---|
| 41 | int main() { |
---|
[477] | 42 | RV x ( "{x }", "2" ); |
---|
[176] | 43 | |
---|
[477] | 44 | |
---|
[176] | 45 | cout << "Test for estimation of a 2D Gaussian mixture" << endl; |
---|
| 46 | /////////////////////////////////// |
---|
| 47 | // Setup model |
---|
| 48 | |
---|
| 49 | vec m1 ( "-2 -2" ); |
---|
| 50 | vec m2 ( "2 2" ); |
---|
| 51 | |
---|
| 52 | fsqmat V1 ( mat ( "1 0.5; 0.5 1" ) ); |
---|
| 53 | fsqmat V2 ( mat ( "2 -0.1; -0.1 2" ) ); |
---|
| 54 | |
---|
[529] | 55 | enorm_fsqmat_ptr C1; |
---|
[504] | 56 | C1->set_rv ( x ); |
---|
| 57 | C1->set_parameters ( m1, V1 ); |
---|
[529] | 58 | enorm_fsqmat_ptr C2; |
---|
[504] | 59 | C2->set_rv ( x ); |
---|
| 60 | C2->set_parameters ( m2, V2 ); |
---|
[176] | 61 | |
---|
[529] | 62 | epdf_array Sim ( 2 ); |
---|
[504] | 63 | Sim ( 0 ) = C1; |
---|
| 64 | Sim ( 1 ) = C2; |
---|
[477] | 65 | emix Simul; |
---|
| 66 | Simul.set_rv ( x ); |
---|
[504] | 67 | Simul.set_parameters ( "0.5 0.6", Sim ); |
---|
[176] | 68 | |
---|
| 69 | // Sample parameters |
---|
| 70 | int ndat = 100; |
---|
[477] | 71 | mat Smp = Simul.sample_m ( ndat ); |
---|
[176] | 72 | |
---|
[477] | 73 | cout << "Simulated means: " << m1 << " and " << m2 << endl; |
---|
| 74 | cout << "Simulated covariances: " << endl << V1 << " and " << endl << V2 << endl; |
---|
[176] | 75 | |
---|
| 76 | ////////////////////////////// |
---|
| 77 | // Prepare estimator manually |
---|
| 78 | |
---|
| 79 | // Initial values of components |
---|
[477] | 80 | mat V0 = 1e-4 * eye ( 3 ); |
---|
| 81 | V0 ( 0, 0 ) *= 100; |
---|
| 82 | V0 ( 1, 1 ) *= 100; |
---|
| 83 | mat Vp = "0 0 1; 0 0 1; 1 1 0"; |
---|
| 84 | Vp *= 5 * 1e-5; |
---|
[176] | 85 | |
---|
[477] | 86 | ARX M1; |
---|
| 87 | M1.set_statistics ( 2, V0 - Vp, 8 ); |
---|
| 88 | ARX M2; |
---|
| 89 | M2.set_statistics ( 2, V0 + Vp, 8 ); |
---|
| 90 | |
---|
[176] | 91 | // Build mixture model |
---|
[477] | 92 | Array<BMEF*> A ( 2 ); |
---|
| 93 | A ( 0 ) = &M1; |
---|
| 94 | A ( 1 ) = &M2; |
---|
| 95 | MixEF Post ( A, "0.5 0.5" ); |
---|
| 96 | cout << "Initial mixture:" << endl; |
---|
| 97 | disp_mix2D ( Post ); |
---|
| 98 | |
---|
[176] | 99 | // Add ones for constant coefficients |
---|
[477] | 100 | mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) ); |
---|
[679] | 101 | Post.bayes ( Data , empty_vec); |
---|
[176] | 102 | |
---|
[477] | 103 | cout << "Posterior mixture:" << endl; |
---|
| 104 | disp_mix2D ( Post ); |
---|
[176] | 105 | |
---|
| 106 | //TEST random initialization |
---|
| 107 | RNG_randomize(); |
---|
[270] | 108 | ARX Aflat; |
---|
[477] | 109 | Aflat.set_statistics ( 2, V0, 7 ); |
---|
[176] | 110 | MixEF RND; |
---|
[477] | 111 | RND.init ( &Aflat, Data, 10 ); // already initialized! |
---|
| 112 | cout << endl << "== Randomly initialized mixture ==" << endl; |
---|
| 113 | cout << endl << "== INIT ==" << endl; |
---|
| 114 | disp_mix2D ( RND ); |
---|
| 115 | mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );; |
---|
[176] | 116 | |
---|
[679] | 117 | RND.bayes ( Data, empty_vec ); |
---|
[477] | 118 | cout << endl << "== BAYES ==" << endl; |
---|
| 119 | disp_mix2D ( RND ); |
---|
[176] | 120 | |
---|
[477] | 121 | |
---|
| 122 | it_file of ( "mixef_test.it", true ); |
---|
| 123 | of << Name ( "Smp" ) << Smp; |
---|
| 124 | of << Name ( "PPdf" ) << grid2D ( Post, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); |
---|
| 125 | of << Name ( "PPdf_I" ) << PPdf_I; |
---|
| 126 | of << Name ( "PPdf_RND" ) << grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); |
---|
| 127 | |
---|
| 128 | cout << endl << "Variables written to file mixef_test.it:" << endl; |
---|
| 129 | cout << "Smp : data sampled from the simulator " << endl; |
---|
| 130 | cout << "PPdf : posterior density on a grid <-5,5><-5,5>" << endl; |
---|
| 131 | cout << "PPdf_RND : posterior density of randomly initilaized mixture on a grid <-5,5><-5,5>" << endl; |
---|
[176] | 132 | } |
---|