root/library/tests/stresssuite/mixtures_stress.cpp @ 722

Revision 722, 3.7 kB (checked in by mido, 15 years ago)

astyler run over all test sources
general_suite added
cleanup of \test directory finished

Line 
1#include "estim/mixtures.h"
2#include "estim/arx.h"
3#include "stat/exp_family.h"
4#include "../mat_checks.h"
5
6using namespace bdm;
7
8//These lines are needed for use of cout and endl
9using std::cout;
10using std::endl;
11
12void 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
17                ( ( egiw* ) ep ( i ) )->mean_mat ( M, Var );
18                cout << "Mean: " << M << endl << "Variance: " << endl << Var << endl;
19        }
20        cout << "Weights: " << ep ( i )->mean() << endl;
21}
22
23mat 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
43TEST ( mixtures_stress ) {
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 );
69        Simul.set_parameters ( "0.5 0.6", Sim );
70
71        // Sample parameters
72        int ndat = 100;
73        mat Smp = Simul.sample_mat ( ndat );
74
75        cout << "Simulated means: " << m1 << " and " << m2 << endl;
76        cout << "Simulated covariances: " << endl << V1 << " and " << endl << V2 << endl;
77
78        //////////////////////////////
79        // Prepare estimator manually
80
81        // Initial values of components
82        mat V0 = 1e-4 * eye ( 3 );
83        V0 ( 0, 0 ) *= 100;
84        V0 ( 1, 1 ) *= 100;
85        mat Vp = "0 0 1; 0 0 1; 1 1 0";
86        Vp *= 5 * 1e-5;
87
88        ARX M1;
89        M1.set_statistics ( 2, V0 - Vp, 8 );
90        ARX M2;
91        M2.set_statistics ( 2, V0 + Vp, 8 );
92
93        // Build mixture model
94        Array<BMEF*> A ( 2 );
95        A ( 0 ) = &M1;
96        A ( 1 ) = &M2;
97        MixEF Post ( A, "0.5 0.5" );
98        cout << "Initial mixture:" << endl;
99        disp_mix2D ( Post );
100
101        // Add ones for constant coefficients
102        mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) );
103        Post.bayes ( Data , empty_vec );
104
105        cout << "Posterior mixture:" << endl;
106        disp_mix2D ( Post );
107
108        //TEST random initialization
109        RNG_randomize();
110        ARX Aflat;
111        Aflat.set_statistics ( 2, V0, 7 );
112        MixEF RND;
113        RND.init ( &Aflat, Data, 10 ); // already initialized!
114        cout << endl << "== Randomly initialized mixture ==" << endl;
115        cout << endl << "== INIT ==" << endl;
116        disp_mix2D ( RND );
117        mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );;
118
119        RND.bayes ( Data, empty_vec );
120        cout << endl << "== BAYES ==" << endl;
121        disp_mix2D ( RND );
122
123
124        it_file of ( "mixef_test.it", true );
125        of << Name ( "Smp" ) << Smp;
126        of << Name ( "PPdf" ) << grid2D ( Post, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );
127        of << Name ( "PPdf_I" ) << PPdf_I;
128        of << Name ( "PPdf_RND" ) << grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );
129
130        cout << endl << "Variables written to file mixef_test.it:" << endl;
131        cout << "Smp  : data sampled from the simulator " << endl;
132        cout << "PPdf : posterior density on a grid <-5,5><-5,5>" << endl;
133        cout << "PPdf_RND : posterior density of randomly initilaized mixture on a grid <-5,5><-5,5>" << endl;
134}
Note: See TracBrowser for help on using the browser.