root/tests/testSmp.cpp @ 125

Revision 125, 2.1 kB (checked in by smidl, 16 years ago)

drobnosti...

RevLine 
[13]1#include <itpp/itbase.h>
[22]2#include <stat/libEF.h>
[104]3#include <stat/emix.h>
[13]4
5using namespace itpp;
6
7//These lines are needed for use of cout and endl
8using std::cout;
9using std::endl;
10
[32]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
[13]22int main() {
23
[32]24        RNG_randomize();
[13]25
[32]26        RV rv("1","{x }","2","0");
[13]27        int N = 10000; //number of samples
28        vec mu0 = "1.5 1.7";
[22]29        mat V0("1.2 0.3; 0.3 5");
30        ldmat R = ldmat(V0);
[32]31       
32        cout << "====== ENorm ====== " <<endl;
[28]33        enorm<ldmat> eN(rv);
34        eN.set_parameters(mu0,R);
[13]35        mat Smp = eN.sample(N);
[32]36
37        disp(mu0,R.to_mat(),Smp);
38
39        cout << "====== MlNorm ====== " <<endl;
40        mat I = eye(2);
41        vec lik(N);
42        mlnorm<ldmat> ML(rv,rv);
43        ML.set_parameters(I,R);
44        Smp = ML.samplecond(mu0,lik,N);
[13]45       
[32]46        disp(mu0,R.to_mat(),Smp);
[13]47
[32]48        cout << "====== EGamma ====== " <<endl; 
49        vec a = "100000,10000";
50        vec b = a/10.0;
51        egamma eG(rv);
52        eG.set_parameters(a,b);
53       
[104]54        cout << eG.eval(a)<<endl;
55        Smp = eG.sampleN(N);
[32]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(rv,rv);
63        double k = 10.0;
64        mG.set_parameters(k);
65
66        Smp=mG.samplecond(mu0,lik,N);
67        disp(mu0,pow(mu0,2.0)/k,Smp);
[104]68       
69        cout << "======= EMix ======== " << endl;
70        emix eMix(rv);
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.sampleN(N);
78        disp(eMix.mean(),zeros(2),Smp);
[32]79
[125]80        cout << "======= MEpdf ======== " << endl;
81        mepdf meMix(rv,rv,&eMix);
[104]82       
[125]83        Smp = meMix.samplecond(mu0,lik,N);
[104]84        disp(eMix.mean(),zeros(2),Smp);
85
[125]86        cout << "======= MMix ======== " << endl;
87        mmix mMix(rv,rv);
88        Array<mpdf*> mComs(2);
89        mComs(0) = &mG;
90        eN.set_mu(vec_2(0.0,0.0));
91        mepdf mEnorm(rv,rv,&eN);
92        mComs(1) = &mEnorm;
93        mMix.set_parameters(vec_2(0.5,0.5),mComs);
94       
95        Smp = mMix.samplecond(mu0,lik,N);
96        disp(mMix._epdf().mean(),zeros(2),Smp);
97
[13]98        //Exit program:
99        return 0;
100
101}
Note: See TracBrowser for help on using the browser.