root/library/tests/testSmp.cpp @ 440

Revision 394, 3.0 kB (checked in by mido, 15 years ago)

1) UI_File renamed to UIFile
2) part UI's documentation
3) stat/mixtures.h renamed to stat/emix.h and related changes applied

  • Property svn:eol-style set to native
Line 
1
2#include "stat/exp_family.h"
3#include "stat/emix.h"
4
5using namespace bdm;
6
7//These lines are needed for use of cout and endl
8using std::cout;
9using std::endl;
10
11void disp(const vec &tmu, const mat &tR,const mat &Smp){
12        int N = Smp.cols();
13        vec Emu = Smp*ones(N) /N ;
14        mat Er = (Smp*Smp.transpose())/N - outer_product(Emu,Emu);
15        cout << "True mu:" << tmu <<endl;
16        cout << "Emp  mu:" << Emu <<endl;
17       
18        cout << "True R:" << tR <<endl;
19        cout << "Emp  R:" << Er <<endl;
20}
21
22int main() {
23
24        RNG_randomize();
25
26        RV x("{x }","2");
27        RV y("{y }","2");
28        int N = 10000; //number of samples
29        vec mu0 = "1.5 1.7";
30        mat V0("1.2 0.3; 0.3 5");
31        ldmat R = ldmat(V0);
32       
33        cout << "====== ENorm ====== " <<endl;
34        enorm<ldmat> eN;
35        eN.set_parameters(mu0,R);
36        mat Smp = eN.sample(N);
37
38        disp(mu0,R.to_mat(),Smp);
39
40        cout << "====== MlNorm ====== " <<endl;
41        mat I = eye(2);
42        mlnorm<ldmat> ML;
43        ML.set_parameters(I,zeros(2),R);
44        Smp = ML.samplecond_m(mu0,N);
45       
46        disp(mu0,R.to_mat(),Smp);
47
48        cout << "====== EGamma ====== " <<endl; 
49        vec a = "100000,10000";
50        vec b = a/10.0;
51        egamma eG;
52        eG.set_parameters(a,b);
53       
54        cout << eG.evallog(a)<<endl;
55        Smp = eG.sample_m(N);
56
57        vec g_mu = elem_div(a,b);
58        vec g_var = elem_div(a,pow(b,2.0));
59        disp(g_mu,diag(g_var),Smp);
60
61        cout << "====== MGamma ====== " <<endl; 
62        mgamma mG;
63        double k = 10.0;
64        mG.set_parameters(k,mu0);
65
66        Smp=mG.samplecond_m(mu0,N);
67        disp(mu0,pow(mu0,2.0)/k,Smp);
68       
69        cout << "======= EMix ======== " << endl;
70        emix eMix;
71        Array<epdf*> Coms(2);
72        Coms(0) = &eG;
73        Coms(1) = &eN;
74       
75        eMix.set_parameters(vec_2(0.5,0.5), Coms);
76        vec smp = eMix.sample();
77        Smp = eMix.sample_m(N);
78        disp(eMix.mean(),zeros(2),Smp);
79
80        cout << "======= MEpdf ======== " << endl;
81        mepdf meMix(&eMix);
82       
83        Smp = meMix.samplecond_m(mu0,N);
84        disp(eMix.mean(),zeros(2),Smp);
85
86        cout << "======= MMix ======== " << endl;
87        mmix mMix;
88        Array<mpdf*> mComs(2);
89        mComs(0) = &mG;
90        eN.set_mu(vec_2(0.0,0.0));
91        mepdf mEnorm(&eN);
92        mComs(1) = &mEnorm;
93        mMix.set_parameters(vec_2(0.5,0.5),mComs);
94       
95        Smp = mMix.samplecond_m(mu0,N);
96        disp(mMix._epdf().mean(),zeros(2),Smp);
97
98        cout << "======= EProd ======== " << endl;
99        // we have to change eG.rv to y
100        eN.set_rv(x);
101        eG.set_rv(y);
102        //create array
103        Array<mpdf*> A(2);
104        mepdf meN(&eN);
105        mepdf meG(&eG);
106        A(0) = &meN;
107        A(1) = &meG;
108       
109        mprod eP(A);
110        mat epV=zeros(4,4);
111        epV.set_submatrix(0,0,V0);
112        epV.set_submatrix(2,2,diag(g_var));
113       
114        vec v0=vec(0);
115        Smp = eP.samplecond(v0,N);
116        disp(concat(eN.mean(),eG.mean()), epV,Smp);
117       
118        cout << "======= eWishart ======== " << endl;
119        mat wM="1.0 0.9; 0.9 1.0";
120        eWishartCh eW; eW.set_parameters(wM/100,100);
121        mat mea=zeros(2,2);
122        mat Ch(2,2);
123        for (int i=0;i<100;i++){Ch=eW.sample_mat(); mea+=Ch.T()*Ch;}
124        cout << mea /100 <<endl;
125       
126        cout << "======= rwiWishart ======== " << endl;
127        rwiWishartCh rwW; rwW.set_parameters(2,0.1,"1 1",0.9);
128        mea=zeros(2,2);
129        mat wMch=chol(wM);
130        for (int i=0;i<100;i++){ 
131                vec tmp=rwW.samplecond(vec(wMch._data(),4));
132                 copy_vector(4,tmp._data(), Ch._data()); 
133                 mea+=Ch.T()*Ch;
134        }
135        cout << mea /100 <<endl;
136        //Exit program:
137        return 0;
138
139}
Note: See TracBrowser for help on using the browser.