root/tests/testSmp.cpp @ 89

Revision 32, 1.4 kB (checked in by smidl, 16 years ago)

test KF : estimation of R in KF is not possible! Likelihood of y_t is growing when R -> 0

Line 
1#include <itpp/itbase.h>
2#include <stat/libEF.h>
3
4using namespace itpp;
5
6//These lines are needed for use of cout and endl
7using std::cout;
8using std::endl;
9
10void disp(const vec &tmu, const mat &tR,const mat &Smp){
11        int N = Smp.cols();
12        vec Emu = Smp*ones(N) /N ;
13        mat Er = (Smp*Smp.transpose())/N - outer_product(Emu,Emu);
14        cout << "True mu:" << tmu <<endl;
15        cout << "Emp  mu:" << Emu <<endl;
16       
17        cout << "True R:" << tR <<endl;
18        cout << "Emp  R:" << Er <<endl;
19}
20
21int main() {
22
23        RNG_randomize();
24
25        RV rv("1","{x }","2","0");
26        int N = 10000; //number of samples
27        vec mu0 = "1.5 1.7";
28        mat V0("1.2 0.3; 0.3 5");
29        ldmat R = ldmat(V0);
30       
31        cout << "====== ENorm ====== " <<endl;
32        enorm<ldmat> eN(rv);
33        eN.set_parameters(mu0,R);
34        mat Smp = eN.sample(N);
35
36        disp(mu0,R.to_mat(),Smp);
37
38        cout << "====== MlNorm ====== " <<endl;
39        mat I = eye(2);
40        vec lik(N);
41        mlnorm<ldmat> ML(rv,rv);
42        ML.set_parameters(I,R);
43        Smp = ML.samplecond(mu0,lik,N);
44       
45        disp(mu0,R.to_mat(),Smp);
46
47        cout << "====== EGamma ====== " <<endl; 
48        vec a = "100000,10000";
49        vec b = a/10.0;
50        egamma eG(rv);
51        eG.set_parameters(a,b);
52       
53        Smp = eG.sample(N);
54
55        vec g_mu = elem_div(a,b);
56        vec g_var = elem_div(a,pow(b,2.0));
57        disp(g_mu,diag(g_var),Smp);
58
59        cout << "====== MGamma ====== " <<endl; 
60        mgamma mG(rv,rv);
61        double k = 10.0;
62        mG.set_parameters(k);
63
64        Smp=mG.samplecond(mu0,lik,N);
65        disp(mu0,pow(mu0,2.0)/k,Smp);
66
67        //Exit program:
68        return 0;
69
70}
Note: See TracBrowser for help on using the browser.