| 29 | | //vec tmu = 0.5 * eN->mean() + 0.5 * mu0; |
| 30 | | |
| 31 | | /* |
| 32 | | RV x ( "{mmixx }", "2" ); |
| 33 | | RV y ( "{mmixy }", "2" ); |
| 34 | | int N = 10000; //number of samples |
| 35 | | vec mu0 ( "1.5 1.7" ); |
| 36 | | mat V0 ( "1.2 0.3; 0.3 5" ); |
| 37 | | ldmat R = ldmat ( V0 ); |
| 38 | | |
| 39 | | enorm_ldmat_ptr eN; |
| 40 | | eN->set_parameters ( mu0, R ); |
| 41 | | |
| 42 | | mgamma_ptr mG; |
| 43 | | double k = 10.0; |
| 44 | | mG->set_parameters ( k, mu0 ); |
| 45 | | |
| 46 | | mmix mMix; |
| 47 | | pdf_array mComs ( 2 ); |
| 48 | | |
| 49 | | // mmix::set_parameters requires the first pdf to be named |
| 50 | | mG->set_rv(x); |
| 51 | | mG->set_rvc(y); |
| 52 | | mComs ( 0 ) = mG; |
| 53 | | |
| 54 | | eN->set_mu ( vec_2 ( 0.0, 0.0 ) ); |
| 55 | | mComs ( 1 ) = eN; |
| 56 | | |
| 57 | | mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs ); |
| 58 | | |
| 59 | | vec tmu = 0.5 * eN->mean() + 0.5 * mu0; |
| 60 | | check_mean ( mMix, mu0, N, tmu, 0.1 ); |
| 61 | | |
| 62 | | mat observedR ( "1.27572 0.778247; 0.778247 3.33129" ); |
| 63 | | check_covariance( mMix, mu0, N, observedR, 0.2 ); |
| 64 | | */ |
| 70 | | /* |
| 71 | | int N = 10000; //number of samples |
| 72 | | vec mu0 ( "1.5 1.7" ); |
| 73 | | mat V0 ( "1.2 0.3; 0.3 5" ); |
| 74 | | ldmat R = ldmat ( V0 ); |
| 75 | | |
| 76 | | enorm_ldmat_ptr eN; |
| 77 | | eN->set_parameters ( mu0, R ); |
| 78 | | |
| 79 | | vec a = "100000,10000"; |
| 80 | | vec b = a / 10.0; |
| 81 | | egamma_ptr eG; |
| 82 | | eG->set_parameters ( a, b ); |
| 83 | | |
| 84 | | emix_ptr eMix; |
| 85 | | epdf_array Coms ( 2 ); |
| 86 | | Coms ( 0 ) = eG; |
| 87 | | Coms ( 1 ) = eN; |
| 88 | | |
| 89 | | eMix->set_parameters ( vec_2 ( 0.5, 0.5 ), Coms ); |
| 90 | | check_mean ( *eMix, mu0, N, eMix->mean(), 0.1 ); |
| 91 | | */ |
| 93 | | |
| 94 | | static void check_mean(pdf &distrib_obj, const vec &mu0, int nsamples, const vec &mean, double tolerance) { |
| 95 | | int tc = 0; |
| 96 | | Array<vec> actual(CurrentContext::max_trial_count); |
| 97 | | do { |
| 98 | | mat smp = distrib_obj.samplecond_mat ( mu0, nsamples ); |
| 99 | | vec emu = smp * ones ( nsamples ) / nsamples ; |
| 100 | | actual( tc ) = emu; |
| 101 | | ++tc; |
| 102 | | } while ( ( tc < CurrentContext::max_trial_count ) && |
| 103 | | !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) ); |
| 104 | | if ( ( tc == CurrentContext::max_trial_count ) && |
| 105 | | ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
| 106 | | UnitTest::MemoryOutStream stream; |
| 107 | | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
| 108 | | stream << "Expected " << mean << " +/- " << tolerance << " but was " << actual; |
| 109 | | |
| 110 | | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
| 111 | | } |
| 112 | | } |
| 113 | | |
| 114 | | static void check_covariance(mmix &distrib_obj, const vec &mu0, int nsamples, const mat &R, double tolerance) { |
| 115 | | int tc = 0; |
| 116 | | Array<mat> actual(CurrentContext::max_trial_count); |
| 117 | | do { |
| 118 | | mat smp = distrib_obj.samplecond_mat ( mu0, nsamples ); |
| 119 | | vec emu = smp * ones ( nsamples ) / nsamples ; |
| 120 | | mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu ); |
| 121 | | actual( tc ) = er; |
| 122 | | ++tc; |
| 123 | | } while ( ( tc < CurrentContext::max_trial_count ) && |
| 124 | | !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) ); |
| 125 | | if ( ( tc == CurrentContext::max_trial_count ) && |
| 126 | | ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
| 127 | | UnitTest::MemoryOutStream stream; |
| 128 | | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
| 129 | | stream << "Expected " << R << " +/- " << tolerance << " but was " << actual; |
| 130 | | |
| 131 | | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
| 132 | | } |
| 133 | | } |