| 1 | #include "base/bdmbase.h" |
|---|
| 2 | #include "base/user_info.h" |
|---|
| 3 | #include "stat/emix.h" |
|---|
| 4 | #include "itpp_ext.h" |
|---|
| 5 | #include "pdf_harness.h" |
|---|
| 6 | #include "mat_checks.h" |
|---|
| 7 | #include "UnitTest++.h" |
|---|
| 8 | |
|---|
| 9 | using namespace bdm; |
|---|
| 10 | |
|---|
| 11 | static void check_mean(pdf &distrib_obj, const vec &mu0, int nsamples, const vec &mean, double tolerance); |
|---|
| 12 | |
|---|
| 13 | static void check_covariance(mmix &distrib_obj, const vec &mu0, int nsamples, const mat &R, double tolerance); |
|---|
| 14 | |
|---|
| 15 | TEST ( mepdf_test ) { |
|---|
| 16 | pdf_harness::test_config ( "mepdf.cfg" ); |
|---|
| 17 | } |
|---|
| 18 | |
|---|
| 19 | TEST ( mgamma_test ) { |
|---|
| 20 | pdf_harness::test_config ( "mgamma.cfg" ); |
|---|
| 21 | } |
|---|
| 22 | |
|---|
| 23 | TEST ( mlnorm_test ) { |
|---|
| 24 | pdf_harness::test_config ( "mlnorm.cfg" ); |
|---|
| 25 | } |
|---|
| 26 | |
|---|
| 27 | TEST ( mprod_test ) { |
|---|
| 28 | pdf_harness::test_config ( "mprod.cfg" ); |
|---|
| 29 | } |
|---|
| 30 | |
|---|
| 31 | // not using pdf_harness because mmix isn't configurable (yet?) |
|---|
| 32 | TEST ( mmix_test ) { |
|---|
| 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 | enorm_ldmat_ptr eN; |
|---|
| 41 | eN->set_parameters ( mu0, R ); |
|---|
| 42 | |
|---|
| 43 | mgamma_ptr mG; |
|---|
| 44 | double k = 10.0; |
|---|
| 45 | mG->set_parameters ( k, mu0 ); |
|---|
| 46 | |
|---|
| 47 | mmix mMix; |
|---|
| 48 | pdf_array mComs ( 2 ); |
|---|
| 49 | |
|---|
| 50 | // mmix::set_parameters requires the first pdf 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 | mComs ( 1 ) = eN; |
|---|
| 57 | |
|---|
| 58 | mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs ); |
|---|
| 59 | |
|---|
| 60 | vec tmu = 0.5 * eN->mean() + 0.5 * mu0; |
|---|
| 61 | check_mean ( mMix, mu0, N, tmu, 0.1 ); |
|---|
| 62 | |
|---|
| 63 | mat observedR ( "1.27572 0.778247; 0.778247 3.33129" ); |
|---|
| 64 | check_covariance( mMix, mu0, N, observedR, 0.2 ); |
|---|
| 65 | } |
|---|
| 66 | |
|---|
| 67 | // not using pdf_harness because emix isn't configurable (yet?) |
|---|
| 68 | TEST ( mepdf_emix_test ) { |
|---|
| 69 | int N = 10000; //number of samples |
|---|
| 70 | vec mu0 ( "1.5 1.7" ); |
|---|
| 71 | mat V0 ( "1.2 0.3; 0.3 5" ); |
|---|
| 72 | ldmat R = ldmat ( V0 ); |
|---|
| 73 | |
|---|
| 74 | enorm_ldmat_ptr eN; |
|---|
| 75 | eN->set_parameters ( mu0, R ); |
|---|
| 76 | |
|---|
| 77 | vec a = "100000,10000"; |
|---|
| 78 | vec b = a / 10.0; |
|---|
| 79 | egamma_ptr eG; |
|---|
| 80 | eG->set_parameters ( a, b ); |
|---|
| 81 | |
|---|
| 82 | emix_ptr eMix; |
|---|
| 83 | epdf_array Coms ( 2 ); |
|---|
| 84 | Coms ( 0 ) = eG; |
|---|
| 85 | Coms ( 1 ) = eN; |
|---|
| 86 | |
|---|
| 87 | eMix->set_parameters ( vec_2 ( 0.5, 0.5 ), Coms ); |
|---|
| 88 | check_mean ( *eMix, mu0, N, eMix->mean(), 0.1 ); |
|---|
| 89 | } |
|---|
| 90 | |
|---|
| 91 | static void check_mean(pdf &distrib_obj, const vec &mu0, int nsamples, const vec &mean, double tolerance) { |
|---|
| 92 | int tc = 0; |
|---|
| 93 | Array<vec> actual(CurrentContext::max_trial_count); |
|---|
| 94 | do { |
|---|
| 95 | mat smp = distrib_obj.samplecond_m ( mu0, nsamples ); |
|---|
| 96 | vec emu = smp * ones ( nsamples ) / nsamples ; |
|---|
| 97 | actual( tc ) = emu; |
|---|
| 98 | ++tc; |
|---|
| 99 | } while ( ( tc < CurrentContext::max_trial_count ) && |
|---|
| 100 | !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) ); |
|---|
| 101 | if ( ( tc == CurrentContext::max_trial_count ) && |
|---|
| 102 | ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
|---|
| 103 | UnitTest::MemoryOutStream stream; |
|---|
| 104 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
|---|
| 105 | stream << "Expected " << mean << " +/- " << tolerance << " but was " << actual; |
|---|
| 106 | |
|---|
| 107 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
|---|
| 108 | } |
|---|
| 109 | } |
|---|
| 110 | |
|---|
| 111 | static void check_covariance(mmix &distrib_obj, const vec &mu0, int nsamples, const mat &R, double tolerance) { |
|---|
| 112 | int tc = 0; |
|---|
| 113 | Array<mat> actual(CurrentContext::max_trial_count); |
|---|
| 114 | do { |
|---|
| 115 | mat smp = distrib_obj.samplecond_m ( mu0, nsamples ); |
|---|
| 116 | vec emu = smp * ones ( nsamples ) / nsamples ; |
|---|
| 117 | mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu ); |
|---|
| 118 | actual( tc ) = er; |
|---|
| 119 | ++tc; |
|---|
| 120 | } while ( ( tc < CurrentContext::max_trial_count ) && |
|---|
| 121 | !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) ); |
|---|
| 122 | if ( ( tc == CurrentContext::max_trial_count ) && |
|---|
| 123 | ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) { |
|---|
| 124 | UnitTest::MemoryOutStream stream; |
|---|
| 125 | UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), __LINE__); |
|---|
| 126 | stream << "Expected " << R << " +/- " << tolerance << " but was " << actual; |
|---|
| 127 | |
|---|
| 128 | UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() ); |
|---|
| 129 | } |
|---|
| 130 | } |
|---|
