root/library/tests/pdf_test.cpp @ 694

Revision 693, 3.8 kB (checked in by mido, 15 years ago)

mpdf renamed to pdf in the whole library

Line 
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
9using namespace bdm;
10
11static void check_mean(pdf &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 ( mepdf_test ) {
16        pdf_harness::test_config ( "mepdf.cfg" );
17}
18
19TEST ( mgamma_test ) {
20        pdf_harness::test_config ( "mgamma.cfg" );
21}
22
23TEST ( mlnorm_test ) {
24        pdf_harness::test_config ( "mlnorm.cfg" );
25}
26
27TEST ( mprod_test ) {
28        pdf_harness::test_config ( "mprod.cfg" );
29}
30
31// not using pdf_harness because mmix isn't configurable (yet?)
32TEST ( 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?)
68TEST ( 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
91static 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
111static 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}
Note: See TracBrowser for help on using the browser.