root/tests/testSmp.cpp @ 170

Revision 170, 2.5 kB (checked in by smidl, 16 years ago)

Mixtures of EF and related changes to libEF and BM

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
[162]26        RV x("{x }","2");
27        RV y("{y }","2");
[13]28        int N = 10000; //number of samples
29        vec mu0 = "1.5 1.7";
[22]30        mat V0("1.2 0.3; 0.3 5");
31        ldmat R = ldmat(V0);
[32]32       
33        cout << "====== ENorm ====== " <<endl;
[145]34        enorm<ldmat> eN(x);
[28]35        eN.set_parameters(mu0,R);
[13]36        mat Smp = eN.sample(N);
[32]37
38        disp(mu0,R.to_mat(),Smp);
39
40        cout << "====== MlNorm ====== " <<endl;
41        mat I = eye(2);
42        vec lik(N);
[145]43        mlnorm<ldmat> ML(x,x);
[32]44        ML.set_parameters(I,R);
45        Smp = ML.samplecond(mu0,lik,N);
[13]46       
[32]47        disp(mu0,R.to_mat(),Smp);
[13]48
[32]49        cout << "====== EGamma ====== " <<endl; 
50        vec a = "100000,10000";
51        vec b = a/10.0;
[145]52        egamma eG(x);
[32]53        eG.set_parameters(a,b);
54       
[104]55        cout << eG.eval(a)<<endl;
56        Smp = eG.sampleN(N);
[32]57
58        vec g_mu = elem_div(a,b);
59        vec g_var = elem_div(a,pow(b,2.0));
60        disp(g_mu,diag(g_var),Smp);
61
62        cout << "====== MGamma ====== " <<endl; 
[145]63        mgamma mG(x,x);
[32]64        double k = 10.0;
65        mG.set_parameters(k);
66
67        Smp=mG.samplecond(mu0,lik,N);
68        disp(mu0,pow(mu0,2.0)/k,Smp);
[104]69       
70        cout << "======= EMix ======== " << endl;
[145]71        emix eMix(x);
[104]72        Array<epdf*> Coms(2);
73        Coms(0) = &eG;
74        Coms(1) = &eN;
75       
76        eMix.set_parameters(vec_2(0.5,0.5), Coms);
77        vec smp = eMix.sample();
78        Smp = eMix.sampleN(N);
79        disp(eMix.mean(),zeros(2),Smp);
[32]80
[125]81        cout << "======= MEpdf ======== " << endl;
[145]82        mepdf meMix(eMix);
[104]83       
[125]84        Smp = meMix.samplecond(mu0,lik,N);
[104]85        disp(eMix.mean(),zeros(2),Smp);
86
[125]87        cout << "======= MMix ======== " << endl;
[145]88        mmix mMix(x,x);
[125]89        Array<mpdf*> mComs(2);
90        mComs(0) = &mG;
91        eN.set_mu(vec_2(0.0,0.0));
[145]92        mepdf mEnorm(eN);
[125]93        mComs(1) = &mEnorm;
94        mMix.set_parameters(vec_2(0.5,0.5),mComs);
95       
96        Smp = mMix.samplecond(mu0,lik,N);
97        disp(mMix._epdf().mean(),zeros(2),Smp);
98
[145]99        cout << "======= EProd ======== " << endl;
100        // we have to change eG.rv to y
[170]101        egamma eGy(x);
102        eGy.set_parameters(a,b);
[145]103        //create array
[162]104        Array<mpdf*> A(2);
105        mepdf meN(eN);
[170]106        mepdf meG(eGy);
[162]107        A(0) = &meN;
108        A(1) = &meG;
[145]109       
[162]110        mprod eP(A);
[145]111        mat epV=zeros(4,4);
112        epV.set_submatrix(0,0,V0);
113        epV.set_submatrix(2,2,diag(g_var));
114       
[162]115        vec v0=vec(0);
116        Smp = eP.samplecond(v0,lik,N);
[170]117        disp(concat(eN.mean(),eGy.mean()), epV,Smp);
[162]118       
[13]119        //Exit program:
120        return 0;
121
122}
Note: See TracBrowser for help on using the browser.