root/applications/mpdm/SYSID09/merg_2a.cpp @ 1305

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

applications compile - but may be broken!

  • Property svn:eol-style set to native
Line 
1#include <estim/arx.h>
2#include <estim/merger.h>
3#include <stat/exp_family.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/*! \file
13Experiment for distributed identification using log-normal merging
14
15The 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.
16
17Lets assume that:
18\dot
19digraph compart{
20
21  edge [fontname="FreeSans",fontsize=10,labelfontname="FreeSans",labelfontsize=10];
22  node [fontname="FreeSans",fontsize=10,shape=record];
23
24rankdir=LR;
25
26 U [label="u",height=0.2,width=0.4,color="white", fillcolor="white", style="filled" fontcolor="black"];
27 AR1 [label="a,b,r",height=0.2,width=0.4,color="black", fillcolor="white", style="filled" fontcolor="black"]
28 U -> AR1 [color="midnightblue",style="solid"];
29 AR2 [label="a,c,r",height=0.2,width=0.4,color="black", fillcolor="white", style="filled" fontcolor="black"]
30 AR1 -> AR2 [color="midnightblue",style="solid",label="y"];
31 Z [label="z",height=0.2,width=0.4,color="white", fillcolor="white", style="filled" fontcolor="black"];
32 AR2 -> Z [color="midnightblue",style="solid"];
33
34 
35}
36\enddot
37*/
38
39int main() {
40        // Setup model
41        RV y ( "{y }" );
42        RV u ( "{u }" );
43        RV z ( "{z }" );
44        RV a ("{a }"); 
45        RV b ("{b }"); RV ab = a; ab.add(b);
46        RV c ("{c }"); RV ac = a; ac.add(c);
47        RV r ("{r }");
48       
49        double at = 0.9;
50        double bt = 5;
51        double ct = -0.5;
52        double rt = 0.50;
53        // Full system
54        vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term
55        vec thz =vec_2 ( at,ct ); //Simulated system - zero for constant term
56        vec Thy = concat ( thy, vec_1(rt) ); //Full parameter
57        vec Thz = concat ( thz, vec_1(rt) ); //Full parameter
58
59        //ARX constructor
60        mat V0 = 0.01*eye ( 3 ); V0 ( 0,0 ) = 1; //
61
62        ARX P1; P1.set_rv(concat(ab,r));
63        P1.set_statistics(1, V0, -1 );
64        P1.set_parameters(0.9);
65        ARX P2; P2.set_rv(concat(ac,r));
66        P2.set_statistics(1, V0, -1 );
67        P2.set_parameters(0.9);
68
69        //Test estimation
70        int ndat = 100;
71        int t;
72
73        // Logging
74        dirfilelog L ( "exp/merg_2a",100 );
75        int Li_Data = L.add ( RV ( "{U Y Z }" ), "" );
76//      int Li_LL   = L.add ( RV ( "{P1 P2 M1 M2 }" ), "LL" );
77        int Li_P1m   = L.add ( concat ( ab,r ), "P1m" );
78        int Li_P2m   = L.add ( concat ( ac,r ), "P2m" );
79        int Li_Mm   = L.add ( concat ( ab,concat(r,c) ), "Mm" );
80        int Li_Th   = L.add ( concat ( ab,concat(c,r) ), "T" );
81        L.init();
82
83        vec Ut ( ndat );
84        vec Yt ( ndat );
85        vec Zt ( ndat );
86        vec yt ( 1 );
87
88        //Proposal
89        enorm<ldmat> g0; g0.set_rv(ab); 
90        g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) );
91        egamma g1;g1.set_rv(r); 
92        g1.set_parameters ( "2  ", "2" );
93        enorm<ldmat> g2; g2.set_rv(c); 
94        g2.set_parameters ( "-1 ",mat ( "1" ) );
95
96        Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2; 
97        eprod G0; G0.set_parameters ( A );
98
99        epdf* proposal=&G0;
100       
101        vec rgry(2);
102        vec rgrz(2);
103        Yt(0) = 0.1;
104        Ut(0) = 0.0;
105        for ( t=1; t<ndat; t++ ) {
106                // True system
107                Ut ( t ) = 0.1*pow ( sin ( ( t/40.0 ) *pi ),3 );
108               
109                rgry(0) = Ut(t); rgry(1) = Ut(t-1);
110                Yt ( t ) = thy*rgry + rt * NorRNG();
111               
112                rgrz(0) = Yt(t); rgrz(1) = Yt(t-1);
113                Zt ( t ) = thz*rgrz + rt * NorRNG();
114
115                // Bayes for all
116                P1.bayes ( concat ( Yt ( t ),rgry ) );
117                P2.bayes ( concat ( Zt ( t ),rgrz ) );
118
119                if (t>50) {bt+=0.1;
120                        thy(1)=bt;}
121                if (t>50) {ct-=0.01;
122                        thz(1)=ct;}
123               
124                // Merge estimates
125                Array<pdf*> A ( 2 ); A ( 0 ) =&P1;A ( 1 ) =&P2;
126                merger_mix M ( A );
127                M.set_parameters ( 20 ,0.99); 
128                M.set_method(LOGNORMAL, 1.2);
129                M._Mix().set_method(QB);
130                //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB);
131/*              char fnm[100];
132                sprintf(fnm,"m2a_dbg%d.it",t);
133                M.debug_file(fnm);*/
134                M.set_support ( *proposal,100 );
135                M.merge();
136                //proposal = M.proposal();
137                //Likelihood
138                yt ( 0 ) = Yt ( t );
139
140                //Logger
141                L.logit(Li_Data, vec_3(Ut(t), Yt(t), Zt(t)));
142                L.logit(Li_P1m, P1._e()->mean());
143                L.logit(Li_P2m, P2._e()->mean());
144                L.logit(Li_Mm, M.mean());
145                L.logit(Li_Th, concat(thy,vec_2(ct,rt*rt)));
146                L.step ( );
147
148                cout << t << "," << endl;
149        }
150        L.finalize( );
151        L.itsave ( "merg_2a.it" );
152        cout << endl;
153}
Note: See TracBrowser for help on using the browser.