root/tests/mixef_test.cpp @ 265

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

create namespace bdm

Line 
1#include <estim/mixef.h>
2#include <estim/arx.h>
3#include <stat/libEF.h>
4using namespace bdm;
5
6//These lines are needed for use of cout and endl
7using std::cout;
8using std::endl;
9
10void disp_mix2D ( const MixEF &Mi ) {
11        const eprod &ep=((const eprod&)(Mi._epdf()));
12        int i;
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;
17        }
18        cout <<"Weights: "<< ep ( i )->mean() <<endl;
19}
20
21mat grid2D ( const MixEF &M,const vec &xb, const vec &yb, int Ngr=100 ) {
22        mat PPdf ( Ngr+1,Ngr+1 );
23        vec rgr ( 3 );
24
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;j=0;
32                for ( double y=yb(0);y<=yb(1);y+=ystep,j++ ) {
33                        rgr ( 1 ) =y;
34                        PPdf ( i,j ) =exp ( M.logpred ( rgr ) );
35                }
36        }
37        return PPdf;
38}
39
40int main() {
41        RV x ( "{x }","2" );
42
43       
44        cout << "Test for estimation of a 2D Gaussian mixture" << endl;
45        ///////////////////////////////////
46        // Setup model
47
48        vec m1 ( "-2 -2" );
49        vec m2 ( "2 2" );
50
51        fsqmat V1 ( mat ( "1 0.5; 0.5 1" ) );
52        fsqmat V2 ( mat ( "2 -0.1; -0.1 2" ) );
53
54        enorm<fsqmat> C1 ( x );
55        C1.set_parameters ( m1,V1 );
56        enorm<fsqmat> C2 ( x );
57        C2.set_parameters ( m2,V2 );
58
59        Array<epdf*> Sim ( 2 ); Sim ( 0 ) =&C1;Sim ( 1 ) =&C2;
60        emix Simul ( x );
61        Simul.set_parameters ( "0.5 0.6", Sim, false );
62
63        // Sample parameters
64        int ndat = 100;
65        mat Smp =Simul.sample_m ( ndat );
66
67        cout << "Simulated means: " << m1 <<" and " << m2 <<endl;
68        cout << "Simulated covariances: " << endl<<V1 <<" and " <<endl<< V2 <<endl;
69
70        //////////////////////////////
71        // Prepare estimator manually
72
73        // Initial values of components
74        mat V0=1e-4*eye ( 3 ); V0 ( 0,0 ) *=100;V0 ( 1,1 ) *=100;
75        mat Vp="0 0 1; 0 0 1; 1 1 0";
76        Vp*=5*1e-5;
77
78        ARX M1 ( RV ( "{Th1 }","6" ),V0 - Vp,8 );
79        ARX M2 ( RV ( "{Th2 }","6" ),V0 + Vp,8 );
80       
81        // Build mixture model
82        Array<BMEF*> A ( 2 ); A ( 0 ) =&M1; A ( 1 ) =&M2;
83        MixEF Post ( A,"0.5 0.5" );
84        cout << "Initial mixture:"<<endl;
85        disp_mix2D(Post);
86       
87        // Add ones for constant coefficients
88        mat Data = concat_vertical ( Smp, ones ( 1,Smp.cols() ) );
89        Post.bayes ( Data );
90
91        cout << "Posterior mixture:" <<endl;
92        disp_mix2D(Post);
93
94        //TEST random initialization
95        RNG_randomize();
96        ARX Aflat=ARX(RV("{Th3 }",vec_1(6)),V0,7);
97        MixEF RND;
98        RND.init(&Aflat,Data,10); // already initialized!
99        cout << endl<< "== Randomly initialized mixture ==" <<endl;
100        cout << endl<< "== INIT ==" <<endl;
101        disp_mix2D(RND);
102        mat PPdf_I=grid2D(RND, vec_2(-5.0,5.0), vec_2(-5.0,5.0));;
103       
104        RND.bayes(Data);
105        cout << endl<< "== BAYES ==" <<endl;
106        disp_mix2D(RND);
107
108       
109        it_file of ( "mixef_test.it",true );
110        of<< Name ( "Smp" ) <<Smp;
111        of<<Name ( "PPdf" ) <<grid2D(Post, vec_2(-5.0,5.0), vec_2(-5.0,5.0));
112        of<<Name ( "PPdf_I" ) <<PPdf_I;
113        of<<Name ( "PPdf_RND" ) <<grid2D(RND, vec_2(-5.0,5.0), vec_2(-5.0,5.0));
114
115        cout << endl <<"Variables written to file mixef_test.it:"<<endl;
116        cout << "Smp  : data sampled from the simulator " <<endl;
117        cout << "PPdf : posterior density on a grid <-5,5><-5,5>"<<endl;
118        cout << "PPdf_RND : posterior density of randomly initilaized mixture on a grid <-5,5><-5,5>"<<endl;
119}
Note: See TracBrowser for help on using the browser.