[447] | 1 | #include "base/bdmbase.h" |
---|
| 2 | #include "base/user_info.h" |
---|
[516] | 3 | #include "stat/emix.h" |
---|
[447] | 4 | #include "itpp_ext.h" |
---|
| 5 | #include "mpdf_harness.h" |
---|
| 6 | #include "mat_checks.h" |
---|
| 7 | #include "UnitTest++.h" |
---|
| 8 | |
---|
| 9 | using namespace bdm; |
---|
| 10 | |
---|
[524] | 11 | static void check_mean(mpdf &distrib_obj, const vec &mu0, int nsamples, const vec &mean, double tolerance); |
---|
[516] | 12 | |
---|
[524] | 13 | static void check_covariance(mmix &distrib_obj, const vec &mu0, int nsamples, const mat &R, double tolerance); |
---|
[516] | 14 | |
---|
[477] | 15 | TEST ( test_mepdf ) { |
---|
| 16 | mpdf_harness::test_config ( "mepdf.cfg" ); |
---|
[462] | 17 | } |
---|
| 18 | |
---|
[477] | 19 | TEST ( test_mgamma ) { |
---|
| 20 | mpdf_harness::test_config ( "mgamma.cfg" ); |
---|
[447] | 21 | } |
---|
[459] | 22 | |
---|
[477] | 23 | TEST ( test_mlnorm ) { |
---|
| 24 | mpdf_harness::test_config ( "mlnorm.cfg" ); |
---|
[459] | 25 | } |
---|
[505] | 26 | |
---|
| 27 | TEST ( test_mprod ) { |
---|
| 28 | mpdf_harness::test_config ( "mprod.cfg" ); |
---|
| 29 | } |
---|
[516] | 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 | |
---|
[524] | 70 | // not using mpdf_harness because emix isn't configurable (yet?) |
---|
| 71 | TEST ( test_mepdf_emix ) { |
---|
| 72 | int N = 10000; //number of samples |
---|
| 73 | vec mu0 ( "1.5 1.7" ); |
---|
| 74 | mat V0 ( "1.2 0.3; 0.3 5" ); |
---|
| 75 | ldmat R = ldmat ( V0 ); |
---|
| 76 | |
---|
| 77 | shared_ptr<enorm<ldmat> > eN = new enorm<ldmat>(); |
---|
| 78 | eN->set_parameters ( mu0, R ); |
---|
| 79 | |
---|
| 80 | vec a = "100000,10000"; |
---|
| 81 | vec b = a / 10.0; |
---|
| 82 | shared_ptr<egamma> eG = new egamma(); |
---|
| 83 | eG->set_parameters ( a, b ); |
---|
| 84 | |
---|
| 85 | shared_ptr<emix> eMix = new emix(); |
---|
| 86 | Array<shared_ptr<epdf> > Coms ( 2 ); |
---|
| 87 | Coms ( 0 ) = eG; |
---|
| 88 | Coms ( 1 ) = eN; |
---|
| 89 | |
---|
| 90 | eMix->set_parameters ( vec_2 ( 0.5, 0.5 ), Coms ); |
---|
| 91 | mepdf meMix ( eMix ); |
---|
| 92 | check_mean ( meMix, mu0, N, eMix->mean(), 0.1 ); |
---|
| 93 | } |
---|
| 94 | |
---|
| 95 | static void check_mean(mpdf &distrib_obj, const vec &mu0, int nsamples, const vec &mean, double tolerance) { |
---|
[516] | 96 | int tc = 0; |
---|
| 97 | Array<vec> actual(CurrentContext::max_trial_count); |
---|
| 98 | do { |
---|
[524] | 99 | mat smp = distrib_obj.samplecond_m ( mu0, nsamples ); |
---|
[516] | 100 | vec emu = smp * ones ( nsamples ) / nsamples ; |
---|
| 101 | actual( tc ) = emu; |
---|
| 102 | ++tc; |
---|
| 103 | } while ( ( tc < CurrentContext::max_trial_count ) && |
---|
| 104 | !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) ); |
---|
| 105 | if ( ( tc == CurrentContext::max_trial_count ) && |
---|
| 106 | ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
---|
| 107 | UnitTest::MemoryOutStream stream; |
---|
| 108 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
---|
| 109 | stream << "Expected " << mean << " +/- " << tolerance << " but was " << actual; |
---|
| 110 | |
---|
| 111 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
---|
| 112 | } |
---|
| 113 | } |
---|
| 114 | |
---|
[524] | 115 | static void check_covariance(mmix &distrib_obj, const vec &mu0, int nsamples, const mat &R, double tolerance) { |
---|
[516] | 116 | int tc = 0; |
---|
| 117 | Array<mat> actual(CurrentContext::max_trial_count); |
---|
| 118 | do { |
---|
[524] | 119 | mat smp = distrib_obj.samplecond_m ( mu0, nsamples ); |
---|
[516] | 120 | vec emu = smp * ones ( nsamples ) / nsamples ; |
---|
| 121 | mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu ); |
---|
| 122 | actual( tc ) = er; |
---|
| 123 | ++tc; |
---|
| 124 | } while ( ( tc < CurrentContext::max_trial_count ) && |
---|
| 125 | !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) ); |
---|
| 126 | if ( ( tc == CurrentContext::max_trial_count ) && |
---|
| 127 | ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
---|
| 128 | UnitTest::MemoryOutStream stream; |
---|
| 129 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
---|
| 130 | stream << "Expected " << R << " +/- " << tolerance << " but was " << actual; |
---|
| 131 | |
---|
| 132 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
---|
| 133 | } |
---|
| 134 | } |
---|