[262] | 1 | |
---|
[386] | 2 | #include "stat/exp_family.h" |
---|
[394] | 3 | #include "stat/emix.h" |
---|
[13] | 4 | |
---|
[254] | 5 | using namespace bdm; |
---|
[13] | 6 | |
---|
| 7 | //These lines are needed for use of cout and endl |
---|
| 8 | using std::cout; |
---|
| 9 | using std::endl; |
---|
| 10 | |
---|
[477] | 11 | void disp ( const vec &tmu, const mat &tR, const mat &Smp ) { |
---|
[32] | 12 | int N = Smp.cols(); |
---|
[477] | 13 | vec Emu = Smp * ones ( N ) / N ; |
---|
| 14 | mat Er = ( Smp * Smp.transpose() ) / N - outer_product ( Emu, Emu ); |
---|
| 15 | cout << "True mu:" << tmu << endl; |
---|
| 16 | cout << "Emp mu:" << Emu << endl; |
---|
| 17 | |
---|
| 18 | cout << "True R:" << tR << endl; |
---|
| 19 | cout << "Emp R:" << Er << endl; |
---|
[32] | 20 | } |
---|
| 21 | |
---|
[13] | 22 | int main() { |
---|
| 23 | |
---|
[32] | 24 | RNG_randomize(); |
---|
[13] | 25 | |
---|
[477] | 26 | RV x ( "{x }", "2" ); |
---|
| 27 | RV y ( "{y }", "2" ); |
---|
[13] | 28 | int N = 10000; //number of samples |
---|
| 29 | vec mu0 = "1.5 1.7"; |
---|
[477] | 30 | mat V0 ( "1.2 0.3; 0.3 5" ); |
---|
| 31 | ldmat R = ldmat ( V0 ); |
---|
| 32 | |
---|
| 33 | cout << "====== ENorm ====== " << endl; |
---|
[270] | 34 | enorm<ldmat> eN; |
---|
[477] | 35 | eN.set_parameters ( mu0, R ); |
---|
| 36 | mat Smp = eN.sample_m ( N ); |
---|
[32] | 37 | |
---|
[477] | 38 | disp ( mu0, R.to_mat(), Smp ); |
---|
[32] | 39 | |
---|
[477] | 40 | cout << "====== MlNorm ====== " << endl; |
---|
| 41 | mat I = eye ( 2 ); |
---|
[270] | 42 | mlnorm<ldmat> ML; |
---|
[477] | 43 | ML.set_parameters ( I, zeros ( 2 ), R ); |
---|
| 44 | Smp = ML.samplecond_m ( mu0, N ); |
---|
[13] | 45 | |
---|
[477] | 46 | disp ( mu0, R.to_mat(), Smp ); |
---|
| 47 | |
---|
| 48 | cout << "====== EGamma ====== " << endl; |
---|
[32] | 49 | vec a = "100000,10000"; |
---|
[477] | 50 | vec b = a / 10.0; |
---|
[270] | 51 | egamma eG; |
---|
[477] | 52 | eG.set_parameters ( a, b ); |
---|
[32] | 53 | |
---|
[477] | 54 | cout << eG.evallog ( a ) << endl; |
---|
| 55 | Smp = eG.sample_m ( N ); |
---|
[32] | 56 | |
---|
[477] | 57 | vec g_mu = elem_div ( a, b ); |
---|
| 58 | vec g_var = elem_div ( a, pow ( b, 2.0 ) ); |
---|
| 59 | disp ( g_mu, diag ( g_var ), Smp ); |
---|
| 60 | |
---|
| 61 | cout << "====== MGamma ====== " << endl; |
---|
[270] | 62 | mgamma mG; |
---|
[32] | 63 | double k = 10.0; |
---|
[477] | 64 | mG.set_parameters ( k, mu0 ); |
---|
[32] | 65 | |
---|
[477] | 66 | Smp = mG.samplecond_m ( mu0, N ); |
---|
| 67 | disp ( mu0, pow ( mu0, 2.0 ) / k, Smp ); |
---|
| 68 | |
---|
[104] | 69 | cout << "======= EMix ======== " << endl; |
---|
[270] | 70 | emix eMix; |
---|
[477] | 71 | Array<epdf*> Coms ( 2 ); |
---|
| 72 | Coms ( 0 ) = &eG; |
---|
| 73 | Coms ( 1 ) = &eN; |
---|
| 74 | |
---|
| 75 | eMix.set_parameters ( vec_2 ( 0.5, 0.5 ), Coms ); |
---|
[104] | 76 | vec smp = eMix.sample(); |
---|
[477] | 77 | Smp = eMix.sample_m ( N ); |
---|
| 78 | disp ( eMix.mean(), zeros ( 2 ), Smp ); |
---|
[32] | 79 | |
---|
[125] | 80 | cout << "======= MEpdf ======== " << endl; |
---|
[477] | 81 | mepdf meMix ( &eMix ); |
---|
[104] | 82 | |
---|
[477] | 83 | Smp = meMix.samplecond_m ( mu0, N ); |
---|
| 84 | disp ( eMix.mean(), zeros ( 2 ), Smp ); |
---|
| 85 | |
---|
[125] | 86 | cout << "======= MMix ======== " << endl; |
---|
[270] | 87 | mmix mMix; |
---|
[477] | 88 | Array<mpdf*> mComs ( 2 ); |
---|
| 89 | mComs ( 0 ) = &mG; |
---|
| 90 | eN.set_mu ( vec_2 ( 0.0, 0.0 ) ); |
---|
| 91 | mepdf mEnorm ( &eN ); |
---|
| 92 | mComs ( 1 ) = &mEnorm; |
---|
| 93 | mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs ); |
---|
[125] | 94 | |
---|
[477] | 95 | Smp = mMix.samplecond_m ( mu0, N ); |
---|
| 96 | disp ( mMix.e()->mean(), zeros ( 2 ), Smp ); |
---|
| 97 | |
---|
[145] | 98 | cout << "======= EProd ======== " << endl; |
---|
| 99 | // we have to change eG.rv to y |
---|
[477] | 100 | eN.set_rv ( x ); |
---|
| 101 | eG.set_rv ( y ); |
---|
| 102 | //create array |
---|
| 103 | Array<mpdf*> A ( 2 ); |
---|
| 104 | mepdf meN ( &eN ); |
---|
| 105 | mepdf meG ( &eG ); |
---|
| 106 | A ( 0 ) = &meN; |
---|
| 107 | A ( 1 ) = &meG; |
---|
| 108 | |
---|
| 109 | mprod eP ( A ); |
---|
| 110 | mat epV = zeros ( 4, 4 ); |
---|
| 111 | epV.set_submatrix ( 0, 0, V0 ); |
---|
| 112 | epV.set_submatrix ( 2, 2, diag ( g_var ) ); |
---|
| 113 | |
---|
| 114 | vec v0 = vec ( 0 ); |
---|
| 115 | Smp = eP.samplecond ( v0, N ); |
---|
| 116 | disp ( concat ( eN.mean(), eG.mean() ), epV, Smp ); |
---|
| 117 | |
---|
[294] | 118 | cout << "======= eWishart ======== " << endl; |
---|
[477] | 119 | mat wM = "1.0 0.9; 0.9 1.0"; |
---|
| 120 | eWishartCh eW; |
---|
| 121 | eW.set_parameters ( wM / 100, 100 ); |
---|
| 122 | mat mea = zeros ( 2, 2 ); |
---|
| 123 | mat Ch ( 2, 2 ); |
---|
| 124 | for ( int i = 0; i < 100; i++ ) { |
---|
| 125 | Ch = eW.sample_mat(); |
---|
| 126 | mea += Ch.T() * Ch; |
---|
| 127 | } |
---|
| 128 | cout << mea / 100 << endl; |
---|
| 129 | |
---|
[294] | 130 | cout << "======= rwiWishart ======== " << endl; |
---|
[477] | 131 | rwiWishartCh rwW; |
---|
| 132 | rwW.set_parameters ( 2, 0.1, "1 1", 0.9 ); |
---|
| 133 | mea = zeros ( 2, 2 ); |
---|
| 134 | mat wMch = chol ( wM ); |
---|
| 135 | for ( int i = 0; i < 100; i++ ) { |
---|
| 136 | vec tmp = rwW.samplecond ( vec ( wMch._data(), 4 ) ); |
---|
| 137 | copy_vector ( 4, tmp._data(), Ch._data() ); |
---|
| 138 | mea += Ch.T() * Ch; |
---|
[294] | 139 | } |
---|
[477] | 140 | cout << mea / 100 << endl; |
---|
[13] | 141 | //Exit program: |
---|
| 142 | return 0; |
---|
| 143 | |
---|
| 144 | } |
---|