[204] | 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> |
---|
[254] | 6 | using namespace bdm; |
---|
[204] | 7 | |
---|
| 8 | //These lines are needed for use of cout and endl |
---|
| 9 | using std::cout; |
---|
| 10 | using std::endl; |
---|
| 11 | |
---|
[288] | 12 | /*! \file |
---|
| 13 | Experiment for distributed identification using log-normal merging |
---|
[213] | 14 | |
---|
[288] | 15 | 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. |
---|
| 16 | |
---|
[213] | 17 | Lets assume that: |
---|
[288] | 18 | \dot |
---|
| 19 | digraph compart{ |
---|
| 20 | |
---|
| 21 | edge [fontname="FreeSans",fontsize=10,labelfontname="FreeSans",labelfontsize=10]; |
---|
| 22 | node [fontname="FreeSans",fontsize=10,shape=record]; |
---|
| 23 | |
---|
| 24 | rankdir=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 |
---|
[213] | 37 | */ |
---|
| 38 | |
---|
[204] | 39 | int main() { |
---|
| 40 | // Setup model |
---|
| 41 | RV y ( "{y }" ); |
---|
[211] | 42 | RV u ( "{u }" ); |
---|
| 43 | RV z ( "{z }" ); |
---|
[213] | 44 | RV a ("{a }"); |
---|
| 45 | RV b ("{b }"); RV ab = a; ab.add(b); |
---|
| 46 | RV c ("{c }"); RV ac = a; ac.add(c); |
---|
[211] | 47 | RV r ("{r }"); |
---|
| 48 | |
---|
[321] | 49 | double at = 0.9; |
---|
| 50 | double bt = 5; |
---|
[213] | 51 | double ct = -0.5; |
---|
[321] | 52 | double rt = 0.50; |
---|
[204] | 53 | // Full system |
---|
[211] | 54 | vec thy =vec_2 ( at,bt ); //Simulated system - zero for constant term |
---|
[213] | 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 |
---|
[204] | 58 | |
---|
| 59 | //ARX constructor |
---|
[321] | 60 | mat V0 = 0.01*eye ( 3 ); V0 ( 0,0 ) = 1; // |
---|
[204] | 61 | |
---|
[288] | 62 | ARX P1; P1.set_rv(concat(ab,r)); |
---|
| 63 | P1.set_statistics(1, V0, -1 ); |
---|
[321] | 64 | P1.set_parameters(0.9); |
---|
[288] | 65 | ARX P2; P2.set_rv(concat(ac,r)); |
---|
[321] | 66 | P2.set_statistics(1, V0, -1 ); |
---|
| 67 | P2.set_parameters(0.9); |
---|
[204] | 68 | |
---|
| 69 | //Test estimation |
---|
| 70 | int ndat = 100; |
---|
| 71 | int t; |
---|
| 72 | |
---|
| 73 | // Logging |
---|
[321] | 74 | dirfilelog L ( "exp/merg_2a",100 ); |
---|
[213] | 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" ); |
---|
[217] | 80 | int Li_Th = L.add ( concat ( ab,concat(c,r) ), "T" ); |
---|
[204] | 81 | L.init(); |
---|
| 82 | |
---|
[213] | 83 | vec Ut ( ndat ); |
---|
[204] | 84 | vec Yt ( ndat ); |
---|
[213] | 85 | vec Zt ( ndat ); |
---|
[204] | 86 | vec yt ( 1 ); |
---|
| 87 | |
---|
| 88 | //Proposal |
---|
[288] | 89 | enorm<ldmat> g0; g0.set_rv(ab); |
---|
[270] | 90 | g0.set_parameters ( "1 1 ",mat ( "1 0; 0 1" ) ); |
---|
[288] | 91 | egamma g1;g1.set_rv(r); |
---|
[270] | 92 | g1.set_parameters ( "2 ", "2" ); |
---|
[288] | 93 | enorm<ldmat> g2; g2.set_rv(c); |
---|
[321] | 94 | g2.set_parameters ( "-1 ",mat ( "1" ) ); |
---|
[204] | 95 | |
---|
[213] | 96 | Array<const epdf*> A ( 3 ); A ( 0 ) = &g0; A ( 1 ) =&g1; A(2) = &g2; |
---|
[288] | 97 | eprod G0; G0.set_parameters ( A ); |
---|
[205] | 98 | |
---|
[321] | 99 | epdf* proposal=&G0; |
---|
| 100 | |
---|
[213] | 101 | vec rgry(2); |
---|
[321] | 102 | vec rgrz(2); |
---|
[213] | 103 | Yt(0) = 0.1; |
---|
| 104 | Ut(0) = 0.0; |
---|
| 105 | for ( t=1; t<ndat; t++ ) { |
---|
[204] | 106 | // True system |
---|
[321] | 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(); |
---|
[204] | 114 | |
---|
| 115 | // Bayes for all |
---|
[321] | 116 | P1.bayes ( concat ( Yt ( t ),rgry ) ); |
---|
| 117 | P2.bayes ( concat ( Zt ( t ),rgrz ) ); |
---|
[204] | 118 | |
---|
[321] | 119 | if (t>50) {bt+=0.1; |
---|
| 120 | thy(1)=bt;} |
---|
| 121 | if (t>50) {ct-=0.01; |
---|
| 122 | thz(1)=ct;} |
---|
| 123 | |
---|
[204] | 124 | // Merge estimates |
---|
[205] | 125 | mepdf eG1 ( P1._e() ); |
---|
| 126 | mepdf eG2 ( P2._e() ); |
---|
[204] | 127 | Array<mpdf*> A ( 2 ); A ( 0 ) =&eG1;A ( 1 ) =&eG2; |
---|
| 128 | merger M ( A ); |
---|
[321] | 129 | M.set_parameters ( 2, 100,20 ,0.99); M._Mix().set_method(QB); |
---|
[205] | 130 | //M2.set_parameters ( 100.0, 1000,3 ); //M2._Mix().set_method(QB); |
---|
[321] | 131 | /* char fnm[100]; |
---|
| 132 | sprintf(fnm,"m2a_dbg%d.it",t); |
---|
| 133 | M.debug_file(fnm);*/ |
---|
| 134 | M.merge ( proposal ); |
---|
[204] | 135 | |
---|
[321] | 136 | //proposal = M.proposal(); |
---|
[204] | 137 | //Likelihood |
---|
| 138 | yt ( 0 ) = Yt ( t ); |
---|
| 139 | |
---|
| 140 | //Logger |
---|
[213] | 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()); |
---|
[288] | 145 | L.logit(Li_Th, concat(thy,vec_2(ct,rt*rt))); |
---|
[204] | 146 | L.step ( ); |
---|
[205] | 147 | |
---|
[321] | 148 | cout << t << "," << endl; |
---|
[204] | 149 | } |
---|
| 150 | L.finalize( ); |
---|
[213] | 151 | L.itsave ( "merg_2a.it" ); |
---|
[205] | 152 | cout << endl; |
---|
[204] | 153 | } |
---|