root/mpdm/merg_2a.cpp @ 287

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

props

  • Property svn:eol-style set to native
Line 
1#include <estim/arx.h>
2#include <estim/merger.h>
3#include <stat/libEF.h>
4#include <stat/loggers.h>
5//#include <stat/merger.h>
6using namespace bdm;
7
8//These lines are needed for use of cout and endl
9using std::cout;
10using std::endl;
11
12/*! The purpose of this experiment is to test merging fragmental pdfs between two participants sharing the same parameter, a. However, this parameter has a different role in each model.
13
14Lets assume that:
15 u -> |a,r| -> y -> |a,r| -> z
16*/
17
18int main() {
19        // Setup model
20        RV y ( "{y }" );
21        RV u ( "{u }" );
22        RV z ( "{z }" );
23        RV a ("{a }"); 
24        RV b ("{b }"); RV ab = a; ab.add(b);
25        RV c ("{c }"); RV ac = a; ac.add(c);
26        RV r ("{r }");
27       
28        double at = 1.5;
29        double bt = 0.5;
30        double ct = -0.5;
31        double rt = 0.30;
32        // Full system
33        vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term
34        vec thz =vec_2 ( at,ct ); //Simulated system - zero for constant term
35        vec Thy = concat ( thy, vec_1(rt) ); //Full parameter
36        vec Thz = concat ( thz, vec_1(rt) ); //Full parameter
37
38        //ARX constructor
39        mat V0 = 0.001*eye ( 3 ); V0 ( 0,0 ) = 1; //
40
41        ARX P1;
42        P1.set_statistics(2, V0, -1 );
43        ARX P2;
44         P2.set_statistics(2, V0, -1 );
45
46        //Test estimation
47        int ndat = 100;
48        int t;
49
50        // Logging
51        dirfilelog L ( "exp/merg_2a",3 );
52        int Li_Data = L.add ( RV ( "{U Y Z }" ), "" );
53//      int Li_LL   = L.add ( RV ( "{P1 P2 M1 M2 }" ), "LL" );
54        int Li_P1m   = L.add ( concat ( ab,r ), "P1m" );
55        int Li_P2m   = L.add ( concat ( ac,r ), "P2m" );
56        int Li_Mm   = L.add ( concat ( ab,concat(r,c) ), "Mm" );
57        int Li_Th   = L.add ( concat ( ab,concat(c,r) ), "T" );
58        L.init();
59
60        vec Ut ( ndat );
61        vec Yt ( ndat );
62        vec Zt ( ndat );
63        vec yt ( 1 );
64
65        //Proposal
66        enorm<ldmat> g0; 
67        g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) );
68        egamma g1; 
69        g1.set_parameters ( "2  ", "2" );
70        enorm<ldmat> g2; 
71        g2.set_parameters ( "1 ",mat ( "1" ) );
72
73        Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2; 
74        eprod G0 ( A );
75
76        vec rgru(2);
77        vec rgry(2);
78        Yt(0) = 0.1;
79        Ut(0) = 0.0;
80        for ( t=1; t<ndat; t++ ) {
81                // True system
82                Ut ( t ) = pow ( sin ( ( t/40.0 ) *pi ),3 );
83                rgru(0) = Ut(t); rgru(1) = Ut(t-1);
84                Yt ( t ) = thy*rgru + rt * NorRNG();
85                rgry(0) = Yt(t); rgry(1) = Yt(t-1);
86                Zt ( t ) = thz*rgry + rt * NorRNG();
87
88                // Bayes for all
89                P1.bayes ( concat ( Yt ( t ),rgru ) );
90                P2.bayes ( concat ( Zt ( t ),rgry ) );
91
92                // Merge estimates
93                mepdf eG1 ( P1._e() );
94                mepdf eG2 ( P2._e() );
95                Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2;
96                merger M ( A );
97                M.set_parameters ( 10, 100,3 ); //M._Mix().set_method(QB);
98                //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
99                M.merge ( &G0 );
100
101                //Likelihood
102                yt ( 0 ) = Yt ( t );
103
104                //Logger
105                L.logit(Li_Data, vec_3(Ut(t), Yt(t), Zt(t)));
106                L.logit(Li_P1m, P1._e()->mean());
107                L.logit(Li_P2m, P2._e()->mean());
108                L.logit(Li_Mm, M.mean());
109                L.logit(Li_Th, concat(thy,vec_2(ct,rt)));
110                L.step ( );
111
112        }
113        L.finalize( );
114        L.itsave ( "merg_2a.it" );
115        cout << endl;
116}
Note: See TracBrowser for help on using the browser.