root/tests/testSmp.cpp @ 196

#include <itpp/itbase.h>
#include <stat/libEF.h>
#include <stat/emix.h>

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

void disp(const vec &tmu, const mat &tR,const mat &Smp){
        int N = Smp.cols();
        vec Emu = Smp*ones(N) /N ;
        mat Er = (Smp*Smp.transpose())/N - outer_product(Emu,Emu);
        cout << "True mu:" << tmu <<endl;
        cout << "Emp  mu:" << Emu <<endl;
        
        cout << "True R:" << tR <<endl;
        cout << "Emp  R:" << Er <<endl;
}

int main() {

        RNG_randomize();

        RV x("{x }","2");
        RV y("{y }","2");
        int N = 10000; //number of samples
        vec mu0 = "1.5 1.7";
        mat V0("1.2 0.3; 0.3 5");
        ldmat R = ldmat(V0);
        
        cout << "====== ENorm ====== " <<endl;
        enorm<ldmat> eN(x);
        eN.set_parameters(mu0,R);
        mat Smp = eN.sample(N);

        disp(mu0,R.to_mat(),Smp);

        cout << "====== MlNorm ====== " <<endl;
        mat I = eye(2);
        vec lik(N);
        mlnorm<ldmat> ML(x,x);
        ML.set_parameters(I,zeros(2),R);
        Smp = ML.samplecond(mu0,lik,N);
        
        disp(mu0,R.to_mat(),Smp);

        cout << "====== EGamma ====== " <<endl; 
        vec a = "100000,10000";
        vec b = a/10.0;
        egamma eG(x);
        eG.set_parameters(a,b);
        
        cout << eG.eval(a)<<endl;
        Smp = eG.sampleN(N);

        vec g_mu = elem_div(a,b);
        vec g_var = elem_div(a,pow(b,2.0));
        disp(g_mu,diag(g_var),Smp);

        cout << "====== MGamma ====== " <<endl; 
        mgamma mG(x,x);
        double k = 10.0;
        mG.set_parameters(k);

        Smp=mG.samplecond(mu0,lik,N);
        disp(mu0,pow(mu0,2.0)/k,Smp);
        
        cout << "======= EMix ======== " << endl;
        emix eMix(x);
        Array<epdf*> Coms(2);
        Coms(0) = &eG;
        Coms(1) = &eN;
        
        eMix.set_parameters(vec_2(0.5,0.5), Coms);
        vec smp = eMix.sample();
        Smp = eMix.sampleN(N);
        disp(eMix.mean(),zeros(2),Smp);

        cout << "======= MEpdf ======== " << endl;
        mepdf meMix(eMix);
        
        Smp = meMix.samplecond(mu0,lik,N);
        disp(eMix.mean(),zeros(2),Smp);

        cout << "======= MMix ======== " << endl;
        mmix mMix(x,x);
        Array<mpdf*> mComs(2);
        mComs(0) = &mG;
        eN.set_mu(vec_2(0.0,0.0));
        mepdf mEnorm(eN);
        mComs(1) = &mEnorm;
        mMix.set_parameters(vec_2(0.5,0.5),mComs);
        
        Smp = mMix.samplecond(mu0,lik,N);
        disp(mMix._epdf().mean(),zeros(2),Smp);

        cout << "======= EProd ======== " << endl;
        // we have to change eG.rv to y
        egamma eGy(x);
        eGy.set_parameters(a,b);
        //create array 
        Array<mpdf*> A(2);
        mepdf meN(eN);
        mepdf meG(eGy);
        A(0) = &meN;
        A(1) = &meG;
        
        mprod eP(A);
        mat epV=zeros(4,4);
        epV.set_submatrix(0,0,V0);
        epV.set_submatrix(2,2,diag(g_var));
        
        vec v0=vec(0);
        Smp = eP.samplecond(v0,lik,N);
        disp(concat(eN.mean(),eGy.mean()), epV,Smp);
        
        //Exit program:
        return 0;

}