Changeset 477 for library/tests/testSmp.cpp
- Timestamp:
- 08/05/09 14:40:03 (15 years ago)
- Files:
-
- 1 modified
-
library/tests/testSmp.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
library/tests/testSmp.cpp
r461 r477 9 9 using std::endl; 10 10 11 void disp (const vec &tmu, const mat &tR,const mat &Smp){11 void disp ( const vec &tmu, const mat &tR, const mat &Smp ) { 12 12 int N = Smp.cols(); 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;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; 20 20 } 21 21 … … 24 24 RNG_randomize(); 25 25 26 RV x ("{x }","2");27 RV y ("{y }","2");26 RV x ( "{x }", "2" ); 27 RV y ( "{y }", "2" ); 28 28 int N = 10000; //number of samples 29 29 vec mu0 = "1.5 1.7"; 30 mat V0 ("1.2 0.3; 0.3 5");31 ldmat R = ldmat (V0);32 33 cout << "====== ENorm ====== " << endl;30 mat V0 ( "1.2 0.3; 0.3 5" ); 31 ldmat R = ldmat ( V0 ); 32 33 cout << "====== ENorm ====== " << endl; 34 34 enorm<ldmat> eN; 35 eN.set_parameters (mu0,R);36 mat Smp = eN.sample_m (N);35 eN.set_parameters ( mu0, R ); 36 mat Smp = eN.sample_m ( N ); 37 37 38 disp (mu0,R.to_mat(),Smp);38 disp ( mu0, R.to_mat(), Smp ); 39 39 40 cout << "====== MlNorm ====== " << endl;41 mat I = eye (2);40 cout << "====== MlNorm ====== " << endl; 41 mat I = eye ( 2 ); 42 42 mlnorm<ldmat> ML; 43 ML.set_parameters(I,zeros(2),R); 44 Smp = ML.samplecond_m(mu0,N); 45 46 disp(mu0,R.to_mat(),Smp); 43 ML.set_parameters ( I, zeros ( 2 ), R ); 44 Smp = ML.samplecond_m ( mu0, N ); 47 45 48 cout << "====== EGamma ====== " <<endl; 46 disp ( mu0, R.to_mat(), Smp ); 47 48 cout << "====== EGamma ====== " << endl; 49 49 vec a = "100000,10000"; 50 vec b = a /10.0;50 vec b = a / 10.0; 51 51 egamma eG; 52 eG.set_parameters(a,b); 53 54 cout << eG.evallog(a)<<endl; 55 Smp = eG.sample_m(N); 52 eG.set_parameters ( a, b ); 56 53 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); 54 cout << eG.evallog ( a ) << endl; 55 Smp = eG.sample_m ( N ); 60 56 61 cout << "====== MGamma ====== " <<endl; 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; 62 62 mgamma mG; 63 63 double k = 10.0; 64 mG.set_parameters (k,mu0);64 mG.set_parameters ( k, mu0 ); 65 65 66 Smp =mG.samplecond_m(mu0,N);67 disp (mu0,pow(mu0,2.0)/k,Smp);68 66 Smp = mG.samplecond_m ( mu0, N ); 67 disp ( mu0, pow ( mu0, 2.0 ) / k, Smp ); 68 69 69 cout << "======= EMix ======== " << endl; 70 70 emix eMix; 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);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 ); 76 76 vec smp = eMix.sample(); 77 Smp = eMix.sample_m (N);78 disp (eMix.mean(),zeros(2),Smp);77 Smp = eMix.sample_m ( N ); 78 disp ( eMix.mean(), zeros ( 2 ), Smp ); 79 79 80 80 cout << "======= MEpdf ======== " << endl; 81 mepdf meMix (&eMix);82 83 Smp = meMix.samplecond_m (mu0,N);84 disp (eMix.mean(),zeros(2),Smp);81 mepdf meMix ( &eMix ); 82 83 Smp = meMix.samplecond_m ( mu0, N ); 84 disp ( eMix.mean(), zeros ( 2 ), Smp ); 85 85 86 86 cout << "======= MMix ======== " << endl; 87 87 mmix mMix; 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);94 95 Smp = mMix.samplecond_m (mu0,N);96 disp (mMix.e()->mean(), zeros(2), Smp);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 ); 94 95 Smp = mMix.samplecond_m ( mu0, N ); 96 disp ( mMix.e()->mean(), zeros ( 2 ), Smp ); 97 97 98 98 cout << "======= EProd ======== " << endl; 99 99 // we have to change eG.rv to y 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 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 118 118 cout << "======= eWishart ======== " << endl; 119 mat wM="1.0 0.9; 0.9 1.0"; 120 eWishartCh eW; eW.set_parameters(wM/100,100); 121 mat mea=zeros(2,2); 122 mat Ch(2,2); 123 for (int i=0;i<100;i++){Ch=eW.sample_mat(); mea+=Ch.T()*Ch;} 124 cout << mea /100 <<endl; 125 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 126 130 cout << "======= rwiWishart ======== " << endl; 127 rwiWishartCh rwW; rwW.set_parameters(2,0.1,"1 1",0.9); 128 mea=zeros(2,2); 129 mat wMch=chol(wM); 130 for (int i=0;i<100;i++){ 131 vec tmp=rwW.samplecond(vec(wMch._data(),4)); 132 copy_vector(4,tmp._data(), Ch._data()); 133 mea+=Ch.T()*Ch; 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; 134 139 } 135 cout << mea / 100 <<endl;140 cout << mea / 100 << endl; 136 141 //Exit program: 137 142 return 0;
