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 | | } |