Changeset 32 for tests

Timestamp:

03/03/08 13:00:32 (17 years ago)

Author:

smidl

Message:

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

Location:

tests

Files:

• CMakeLists.txt (modified) (3 diffs)
• test0.cpp (modified) (1 diff)
• testEpdf.cpp (added)
• testKF.cpp (modified) (3 diffs)
• testKF_QR.cpp (added)
• testPF.cpp (modified) (2 diffs)
• testResample.cpp (added)
• testSmp.cpp (modified) (2 diffs)

tests/CMakeLists.txt

@@ -6 +6 @@
 
 ## Save all needed libraries in variable BdmLibs
-SET(BdmLibs bdm itpp)
+SET(BdmLibs bdm itpp_debug)
+#SET(BdmLibs bdm itpp)
 
 IF(WIN32)
@@ -23 +24 @@
 add_executable (testPF testPF.cpp)
 add_executable (testSmp testSmp.cpp)
+add_executable (testEpdf testEpdf.cpp)
+add_executable (testResample testResample.cpp)
+add_executable (testKF_QR testKF_QR.cpp)
 
 # Link the executable to the Hello library.
@@ -29 +33 @@
 target_link_libraries (testPF ${BdmLibs})
 target_link_libraries (testSmp ${BdmLibs})
+target_link_libraries (testEpdf ${BdmLibs})
+target_link_libraries (testResample ${BdmLibs})
+target_link_libraries (testKF_QR ${BdmLibs})

tests/test0.cpp

@@ -12 +12 @@
 {
 
-        RV th = RV ( "1 2","{a b }","1 1","0 0","0 0" );
+        RV th = RV ( "1 2","{a b }","1 1","0 0");
         RV r = RV ( "3 4" );
 
         cout << th << r << endl;
 
-        ldmat ld = ldmat("1 0;-0.1 1","1.1 1.3");
+        ldmat ld = ldmat("1 0;0.5 1","1.1 1.3");
         vec v = "1 -0.1";
         

tests/testKF.cpp

@@ -18 +18 @@
         it_file fin( "testKF.it" );
 
-        mat Dt, Xt,Xt2,XtE;
+        mat Dt, Xt,Xt2,XtE,Xtf;
         int Ndat;
 
@@ -42 +42 @@
         Xt=zeros( 2,Ndat );
         Xt2=zeros( 2,Ndat );
+        Xtf=zeros( 2,Ndat );
         XtE=zeros( 2,Ndat );
 
 //      cout << KF;
-        RV rx("1","{x}","2","0","0");
-        RV ru("2","{u}","1","0","0");
-        RV ry("3","{y}","1","0","0");
+        RV rx("1","{x}","2","0");
+        RV ru("2","{u}","1","0");
+        RV ry("3","{y}","1","0");
         //
         Kalman<ldmat> KF(rx,ry,ru);
@@ -65 +66 @@
         KFE.set_est(mu0,P0);
 
-        Xt.set_col( 0,*((enorm<ldmat>*)(KF._epdf()))->_mu() );
+        epdf& KFep = KF._epdf();
+        epdf& KFfep = KFf._epdf();
+        epdf& KFEep = KFE._epdf();
+        
+        Xt.set_col( 0,KFep.mean() );
+        Xtf.set_col( 0,KFfep.mean() );
         Xt2.set_col( 0,KF2.mu );
-        XtE.set_col( 0,*((enorm<ldmat>*)(KFE._epdf()))->_mu() );
+        XtE.set_col( 0,KFEep.mean() );
         for ( int t=1;t<Ndat;t++ ) {
                 KFf.bayes( Dt.get_col( t ));
                 KF.bayes( Dt.get_col( t ));
                 KF2.bayes( Dt.get_col( t ));
-//              KFE.bayes( Dt.get_col( t ));
-                Xt.set_col(t,*((enorm<ldmat>*)(KF._epdf()))->_mu());
+                KFE.bayes( Dt.get_col( t ));
+                Xt.set_col( t,KFep.mean() );
+                Xtf.set_col( t,KFfep.mean() );
                 Xt2.set_col(t,KF2.mu);
-                XtE.set_col(t,*((enorm<ldmat>*)(KFE._epdf()))->_mu());
+                XtE.set_col( t,KFEep.mean() );
         }
 
         it_file fou( "testKF_res.it" );
         fou << Name("xth") << Xt;
+        fou << Name("xthf") << Xtf;
         fou << Name("xth2") << Xt2;
         fou << Name("xthE") << XtE;

tests/testPF.cpp

@@ -13 +13 @@
 
         RV x("1");
-        RV xm=x; xm.t(-1);
+        RV xm=x; xm.t(-1);const
         RV y("2");
         
@@ -19 +19 @@
         vec vR = "1";
         ldmat R(vR);
+                
+        eEmp emp(x,10);
+        euni eun(x);
+        eun.set_parameters("0","1");
+        emp.set_parameters(ones(10),&eun);
+        vec &v=emp._w();
+        Array<vec> &S=emp._samples();
         
+        for (int i=0;i<10;i++){ v(i) = exp(-0.5*sum(pow(S(i)-1,2.0))*10);}
+        v/=sum(v);
         
-        vec ptc = randn(10);
-        vec w = exp(-0.5*(pow(ptc,2))/0.2);
+        cout << "p:" << S << endl;
+        cout << "w:" << v << endl;
         
-        cout << "p:" << ptc << endl;
-        cout << "w:" << w << endl;
-        
-        PF pf(w);
-        ivec ind = pf.resample();
+        ivec ind = emp.resample();
         
         cout << ind << endl;
-        /*
-        mlnorm<ldmat> obs(x,xm,A,R);
-        mlnorm<ldmat> par(y,x,A,R);
         
-        TrivialPF TPF(obs,par,10);
-        */
         //Exit program:
         return 0;

tests/testSmp.cpp

@@ -8 +8 @@
 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();
+        RNG_randomize();
 
-        RV rv("1","{x }","2","0","0");
+        RV rv("1","{x }","2","0");
         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(rv);
         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(rv,rv);
+        ML.set_parameters(I,R);
+        Smp = ML.samplecond(mu0,lik,N);
         
-        mat Smp = eN.sample(N);
-        cout << "True:" <<endl;
-        cout << "mu:" << mu0 ;
-        cout << "R:" << R ;
-        cout << "R:" << R.to_mat() ;
+        disp(mu0,R.to_mat(),Smp);
+
+        cout << "====== EGamma ====== " <<endl; 
+        vec a = "100000,10000";
+        vec b = a/10.0;
+        egamma eG(rv);
+        eG.set_parameters(a,b);
         
-        vec Emu = Smp*ones(N) /N ;
-        cout << "Empirical:" <<endl;
-        cout << "mu:" << Emu;
-        cout << "R:" << (Smp*Smp.transpose())/N - Emu*Emu.transpose() <<endl;
+        Smp = eG.sample(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(rv,rv);
+        double k = 10.0;
+        mG.set_parameters(k);
+
+        Smp=mG.samplecond(mu0,lik,N);
+        disp(mu0,pow(mu0,2.0)/k,Smp);
 
         //Exit program:
@@ -35 +69 @@
 
 }
+