root/mpdm/merger_iter_cond.cpp @ 289

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

props

  • Property svn:eol-style set to native
Line 
1
2#include <stat/libEF.h>
3#include <estim/merger.h>
4
5using namespace bdm;
6
7//These lines are needed for use of cout and endl
8using std::cout;
9using std::endl;
10
11int main() {
12
13        RNG_randomize();
14
15        RV y ( "{y }","1" );
16        RV u1 ( "{u1 }","1" );
17        RV u2 ( "{u2 }","1" );
18
19        RV all = y;
20        all.add ( u1 );
21        all.add ( u2 );
22
23        mlnorm<fsqmat> f1; f1.set_rv( y ); f1.set_rvc( u1 );
24        mlnorm<fsqmat> f2; f2.set_rv( y ); f2.set_rvc( u2 );
25
26        //Differneces in constant term are essential
27        f1.set_parameters ( "0.4","1",mat ( "0.04" ) );
28        f2.set_parameters ( "0.2","-1",mat ( "0.08" ) );
29
30        Array<mpdf* > A ( 2 );
31        A ( 0 ) =&f1;
32        A ( 1 ) =&f2;
33
34        merger M ( A );
35        enorm<fsqmat> g0; g0.set_rv ( all );
36        g0.set_parameters ( vec ( "1 -1 1" ),3*eye ( 3 ));// +1*ones ( 3,3 ) );
37
38        M.set_parameters ( 1.2,1000,1 );
39
40        int Ntrials=100;
41        vec A1s ( Ntrials );
42        vec A2s ( Ntrials );
43        vec R1s ( Ntrials );
44        vec R2s ( Ntrials );
45        vec C1s ( Ntrials );
46        vec C2s ( Ntrials );
47
48        for ( int it=0;it<Ntrials;it++ ) {
49                M.merge ( &g0 );
50
51                MixEF &MM = M._Mix();
52                epdf* MP = MM._Coms ( 0 )->epredictor ( );
53
54                RV yu1 = y; yu1.add ( u1 );
55                RV yu2 = y; yu2.add ( u2 );
56                enorm<ldmat>* P1m= ( enorm<ldmat>* ) MP->marginal ( yu1 );
57                enorm<ldmat>* P2m= ( enorm<ldmat>* ) MP->marginal ( yu2 );
58                mlnorm<ldmat>* P1c= ( mlnorm<ldmat>* ) ( P1m->condition ( y ) );
59                mlnorm<ldmat>* P2c= ( mlnorm<ldmat>* ) ( P2m->condition ( y ) );
60
61                A1s(it) = P1c->_A()(0,0);
62                A2s(it) = P2c->_A()(0,0);
63                R1s(it) = P1c->_R()(0,0);
64                R2s(it) = P2c->_R()(0,0);
65                C1s(it) = P1c->_mu_const()(0);
66                C2s(it) = P2c->_mu_const()(0);
67
68        }
69        double A1mean = sum(A1s)/Ntrials;
70        double A2mean = sum(A2s)/Ntrials;
71        double C1mean = sum(C1s)/Ntrials;
72        double C2mean = sum(C2s)/Ntrials;
73        double R1mean = sum(R1s)/Ntrials;
74        double R2mean = sum(R2s)/Ntrials;
75        cout << "A1: " << A1mean << " +- " << 2*(sum_sqr(A1s)/Ntrials-A1mean*A1mean) <<endl;
76        cout << "A2: " << A2mean << " +- " << 2*(sum_sqr(A2s)/Ntrials-A2mean*A2mean) <<endl;
77        cout << "C1: " << C1mean << " +- " << 2*(sum_sqr(C1s)/Ntrials-C1mean*C1mean) <<endl;
78        cout << "C2: " << C2mean << " +- " << 2*(sum_sqr(C2s)/Ntrials-C2mean*C2mean) <<endl;
79        cout << "R1: " << R1mean << " +- " << 2*(sum_sqr(R1s)/Ntrials-R1mean*R1mean) <<endl;
80        cout << "R2: " << R2mean << " +- " << 2*(sum_sqr(R2s)/Ntrials-R2mean*R2mean) <<endl;
81}
Note: See TracBrowser for help on using the browser.