root/applications/mpdm/TR2244/merger_iter_cond.cpp @ 1439

Revision 702, 3.0 kB (checked in by smidl, 15 years ago)

applications compile - but may be broken!

  • Property svn:eol-style set to native
Line 
1
2#include <stat/exp_family.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        enorm<fsqmat>  f2; f2.set_rv( y );// f2.set_rvc( u2 );
25
26        //Differneces in constant term are essential
27        f1.set_parameters ( "1","1",mat ( "0.01" ) );
28        f2.set_parameters ( "2", mat ( "0.02" ) );
29
30        Array<pdf* > A ( 2 );
31        A ( 0 ) =&f1;
32        A ( 1 ) =&f2;
33
34        merger_mix M ( A );
35        M.set_debug_file("iter_cond_debug.it");
36        enorm<ldmat> g0; g0.set_rv ( all );
37        mat Cov=10*eye ( 2 );
38/*      Cov(1,2)=0.29;
39        Cov(2,1)=0.29;*/
40        g0.set_parameters ( vec ( "2 1" ),Cov);// +1*ones ( 3,3 ) );
41
42//      enorm<ldmat>* teste=g0.marginal(concat(u1,u2));
43//      mlnorm<ldmat>* testm=(mlnorm<ldmat>*)teste->condition(u2);
44       
45        M.set_parameters (1,0.5);
46        M.set_method(LOGNORMAL, 1.2);
47       
48        Array<vec> YUU(2); 
49        YUU(0)=linspace(1.5,2.5,20);
50        YUU(1)=linspace(0.5,1.5,20);
51       
52        YUU(0)=linspace(0,4,20);
53        YUU(1)=linspace(-1,3,20);
54//YUU(2)=linspace(-3,3,20);
55
56        //M.set_grid(YUU);
57       
58        int Ntrials=1;
59        vec A1s ( Ntrials );
60//      vec A2s ( Ntrials );
61        vec R1s ( Ntrials );
62        vec R2s ( Ntrials );
63        vec C1s ( Ntrials );
64        vec C2s ( Ntrials );
65
66        for ( int it=0;it<Ntrials;it++ ) {
67                M.set_support( g0,1000);
68                M.merge();
69               
70                MixEF &MM = M._Mix();
71                epdf* MP = MM._Coms ( 0 )->epredictor ( );
72                MP->set_rv(all);
73               
74                //Missing piece
75                enorm<ldmat>* miss=     (enorm<ldmat>*) MP->marginal(u1);
76                cout << "Missing: " << miss->_mu() << " +- " << miss->_R().to_mat() << endl;
77
78                RV yu1 = y; yu1.add ( u1 );
79//              RV yu2 = y; yu2.add ( u2 );
80                enorm<ldmat>* P1m= ( enorm<ldmat>* ) MP->marginal ( yu1 );
81                enorm<ldmat>* P2m= ( enorm<ldmat>* ) MP->marginal ( y );
82                mlnorm<ldmat>* P1c= ( mlnorm<ldmat>* ) ( P1m->condition ( y ) );
83//              mlnorm<ldmat>* P2c= ( mlnorm<ldmat>* ) ( P2m->condition ( y ) );
84
85                A1s(it) = P1c->_A()(0,0);
86//              A2s(it) = P2c->_A()(0,0);
87                R1s(it) = P1c->_R()(0,0);
88                R2s(it) = P2m->_R().to_mat()(0,0);
89                C1s(it) = P1c->_mu_const()(0);
90                C2s(it) = P2m->_mu()(0);
91
92                cout << "mean: " << MM._Coms(0)->_e()->mean() <<endl;
93        }
94       
95       
96       
97        cout << "A1s:" <<A1s<<endl;
98        cout << "C2s:" <<C2s<<endl;
99        double A1mean = sum(A1s)/Ntrials;
100//      double A2mean = sum(A2s)/Ntrials;
101        double C1mean = sum(C1s)/Ntrials;
102        double C2mean = sum(C2s)/Ntrials;
103        double R1mean = sum(R1s)/Ntrials;
104        double R2mean = sum(R2s)/Ntrials;
105        cout << "A1: " << A1mean << " +- " << 2*sqrt(sum_sqr(A1s)/Ntrials-A1mean*A1mean) <<endl;
106//      cout << "A2: " << A2mean << " +- " << 2*sqrt(sum_sqr(A2s)/Ntrials-A2mean*A2mean) <<endl;
107        cout << "C1: " << C1mean << " +- " << 2*sqrt(sum_sqr(C1s)/Ntrials-C1mean*C1mean) <<endl;
108        cout << "C2: " << C2mean << " +- " << 2*sqrt(sum_sqr(C2s)/Ntrials-C2mean*C2mean) <<endl;
109        cout << "R1: " << R1mean << " +- " << 2*sqrt(sum_sqr(R1s)/Ntrials-R1mean*R1mean) <<endl;
110        cout << "R2: " << R2mean << " +- " << 2*sqrt(sum_sqr(R2s)/Ntrials-R2mean*R2mean) <<endl;
111       
112       
113}
Note: See TracBrowser for help on using the browser.