root/library/tests/mpdf_test.cpp @ 625

Revision 529, 3.8 kB (checked in by vbarta, 15 years ago)

defined *_ptr wrappers of shared pointers

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