[757] | 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> |
---|
| 6 | using namespace bdm; |
---|
| 7 | |
---|
| 8 | //These lines are needed for use of cout and endl |
---|
| 9 | using std::cout; |
---|
| 10 | using std::endl; |
---|
| 11 | |
---|
| 12 | /*! \file |
---|
| 13 | Experiment for distributed identification using log-normal merging |
---|
| 14 | |
---|
| 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 | |
---|
| 17 | Lets assume that: |
---|
| 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 |
---|
| 37 | */ |
---|
| 38 | |
---|
| 39 | int 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 | } |
---|