| 30 | |
| 31 | // not using mpdf_harness because mmix isn't configurable (yet?) |
| 32 | TEST ( test_mmix ) { |
| 33 | RV x ( "{mmixx }", "2" ); |
| 34 | RV y ( "{mmixy }", "2" ); |
| 35 | int N = 10000; //number of samples |
| 36 | vec mu0 ( "1.5 1.7" ); |
| 37 | mat V0 ( "1.2 0.3; 0.3 5" ); |
| 38 | ldmat R = ldmat ( V0 ); |
| 39 | |
| 40 | shared_ptr<enorm<ldmat> > eN = new enorm<ldmat>(); |
| 41 | eN->set_parameters ( mu0, R ); |
| 42 | |
| 43 | shared_ptr<mgamma> mG = new mgamma(); |
| 44 | double k = 10.0; |
| 45 | mG->set_parameters ( k, mu0 ); |
| 46 | |
| 47 | mmix mMix; |
| 48 | Array<shared_ptr<mpdf> > mComs ( 2 ); |
| 49 | |
| 50 | // mmix::set_parameters requires the first mpdf to be named |
| 51 | mG->set_rv(x); |
| 52 | mG->set_rvc(y); |
| 53 | mComs ( 0 ) = mG; |
| 54 | |
| 55 | eN->set_mu ( vec_2 ( 0.0, 0.0 ) ); |
| 56 | shared_ptr<mepdf> mEnorm = new mepdf ( eN ); |
| 57 | mComs ( 1 ) = mEnorm; |
| 58 | |
| 59 | mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs ); |
| 60 | |
| 61 | double tolerance = 0.1; |
| 62 | |
| 63 | vec tmu = 0.5 * eN->mean() + 0.5 * mu0; |
| 64 | check_mean ( mMix, mu0, N, tmu, tolerance ); |
| 65 | |
| 66 | mat observedR ( "1.27572 0.778247; 0.778247 3.33129" ); |
| 67 | check_covariance( mMix, mu0, N, observedR, tolerance); |
| 68 | } |
| 69 | |
| 70 | static void check_mean(mmix &mMix, const vec &mu0, int nsamples, const vec &mean, double tolerance) { |
| 71 | int tc = 0; |
| 72 | Array<vec> actual(CurrentContext::max_trial_count); |
| 73 | do { |
| 74 | mat smp = mMix.samplecond_m ( mu0, nsamples ); |
| 75 | vec emu = smp * ones ( nsamples ) / nsamples ; |
| 76 | actual( tc ) = emu; |
| 77 | ++tc; |
| 78 | } while ( ( tc < CurrentContext::max_trial_count ) && |
| 79 | !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) ); |
| 80 | if ( ( tc == CurrentContext::max_trial_count ) && |
| 81 | ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
| 82 | UnitTest::MemoryOutStream stream; |
| 83 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
| 84 | stream << "Expected " << mean << " +/- " << tolerance << " but was " << actual; |
| 85 | |
| 86 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
| 87 | } |
| 88 | } |
| 89 | |
| 90 | static void check_covariance(mmix &mMix, const vec &mu0, int nsamples, const mat &R, double tolerance) { |
| 91 | int tc = 0; |
| 92 | Array<mat> actual(CurrentContext::max_trial_count); |
| 93 | do { |
| 94 | mat smp = mMix.samplecond_m ( mu0, nsamples ); |
| 95 | vec emu = smp * ones ( nsamples ) / nsamples ; |
| 96 | mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu ); |
| 97 | actual( tc ) = er; |
| 98 | ++tc; |
| 99 | } while ( ( tc < CurrentContext::max_trial_count ) && |
| 100 | !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) ); |
| 101 | if ( ( tc == CurrentContext::max_trial_count ) && |
| 102 | ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
| 103 | UnitTest::MemoryOutStream stream; |
| 104 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
| 105 | stream << "Expected " << R << " +/- " << tolerance << " but was " << actual; |
| 106 | |
| 107 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
| 108 | } |
| 109 | } |