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