root/tests/emix_test.cpp @ 309

Revision 278, 2.7 kB (checked in by smidl, 16 years ago)

props

  • Property svn:eol-style set to native
Line 
1#include <stat/libEF.h>
2#include <stat/emix.h>
3using namespace bdm;
4
5//These lines are needed for use of cout and endl
6using std::cout;
7using std::endl;
8
9double normcoef ( const epdf* ep,const vec &xb, const vec &yb, int Ngr=100 ) {
10        mat PPdf ( Ngr+1,Ngr+1 );
11        vec rgr ( 2 );
12
13        int i=0,j=0;
14        double xstep= ( xb ( 1 )-xb ( 0 ) ) /Ngr;
15        double ystep= ( yb ( 1 )-yb ( 0 ) ) /Ngr;
16
17        for ( double x=xb ( 0 );x<=xb ( 1 );x+= xstep,i++ ) {
18                rgr ( 0 ) =x;j=0;
19                for ( double y=yb ( 0 );y<=yb ( 1 );y+=ystep,j++ ) {
20                        rgr ( 1 ) =y;
21                        PPdf ( i,j ) =exp ( ep->evallog ( rgr ) );
22                }
23        }
24        return sumsum ( PPdf ) *xstep*ystep;
25}
26
27int main() {
28        RNG_randomize();
29
30        RV x ( "{x }" );
31        RV y ( "{y }" );
32        RV xy=concat(x,y);
33       
34        enorm<ldmat> E1; E1.set_rv ( xy );
35               
36        it_file iti("xxx");
37        mat R; vec mu;
38        iti >> Name("R") >>R;
39        iti >> Name("mu") >>mu;
40
41        cout << mu << R <<endl;
42        E1.set_parameters( mu, R);//"1.00054 1.0455" , mat ( "0.740142 -0.259015; -0.259015 1.0302" ));
43        enorm<ldmat> E2;E2.set_rv ( xy );
44        E2.set_parameters( "-1.2 -0.1" , mat ( "1 0.4; 0.4 0.5" ));
45
46        Array<epdf*> A1(1);
47        A1(0) = &E1;
48       
49        emix M1; M1.set_rv(xy);
50        M1.set_parameters(vec("1"), A1, false);
51        cout << "======== test if ARX and emix with one ARX are the same ==="<<endl;
52       
53        epdf* Mm = M1.marginal(y);
54        epdf* Am = E1.marginal(y);
55        mpdf* Mc = M1.condition(y);
56        mpdf* Ac = E1.condition(y);
57       
58        cout << *((mlnorm<ldmat>*)Ac) <<endl;
59       
60        cout << "Mix marg at 0: " << Mm->evallog(vec_1(0.0)) <<endl;
61        cout << "ARX marg at 0: " << Am->evallog(vec_1(0.0)) <<endl;
62        cout << "Mix cond at 0,0: " << Mc->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
63        cout << "ARX cond at 0,0: " << Ac->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
64       
65        cout << "======== Mixture with two components ==="<<endl;
66        Array<epdf*> A2(2);
67        A2(0) = &E1;
68        A2(1) = &E2;
69       
70        emix M2; M2.set_rv(xy);
71        M2.set_parameters(vec("1"), A2, false);
72       
73       
74        cout << "Mixture normalization: " << normcoef(&M2, vec("-3 3 "), vec("-3 3 ")) <<endl;
75       
76        int N=3;
77        vec ll ( N );
78        vec ll2 ( N );
79        mat Smp=M2.sample_m ( N );
80        ll = M2.evallog_m(Smp);
81       
82        vec Emu = sum ( Smp,2 ) /N;
83        mat Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );
84
85        cout << "original empirical mean: " <<Emu <<endl;
86        cout << "original empirical variance: " <<Er <<endl;
87
88        epdf* Mg = (epdf*)M2.marginal ( y );
89        mpdf* Cn = (mpdf*)M2.condition ( x );
90
91        it_file it("emix_test.it");
92        it << Name("Smp") << Smp;
93
94        cout<< "marginal mean: " << Mg->mean() <<endl;
95
96        cout << "========== putting them back together ======= "<<endl;
97        mepdf mMg ( Mg );
98        Array<mpdf*> AA ( 2 );
99        AA ( 0 ) = Cn;
100        AA ( 1 ) = &mMg;
101        mprod mEp ( AA );
102
103        for (int j=0; j<N; j++){
104                ll2(j)=mEp.evallogcond (Smp.get_col(j), vec ( 0 ));
105        }
106        it << Name("ll") << ll;
107        it << Name("ll2") << ll2;
108       
109       
110}
Note: See TracBrowser for help on using the browser.