Changeset 1196 for library

Timestamp:

09/24/10 21:36:43 (14 years ago)

Author:

smidl

Message:

monitor Neff + test enorm

Location:

library

Files:

• bdm/estim/particles.cpp (modified) (1 diff)
• bdm/estim/particles.h (modified) (6 diffs)
• bdm/math/functions.h (modified) (2 diffs)
• bdm/stat/exp_family.h (modified) (2 diffs)
• tests/testsuite/enorm.cfg (modified) (1 diff)

library/bdm/estim/particles.cpp

@@ -46 +46 @@
                 log_level.store( logweights, w);
         }
+        if (log_level[logNeff]) {
+                log_level.store( logNeff, 1.0 / ( w * w ));
+        }
         
     if ( do_resampling() ) {

library/bdm/estim/particles.h

@@ -258 +258 @@
     //! \var log_level_enums logmeans
     //! means of particles will be logged
-    LOG_LEVEL(PF,logweights,logmeans,logvars);
+    LOG_LEVEL(PF,logweights,logmeans,logvars,logNeff);
 
     class pf_mix: public emix_base {
@@ -392 +392 @@
     void log_register ( bdm::logger& L, const string& prefix ) {
         BM::log_register(L,prefix);
-        if (log_level[logweights]) {
-            L.add_vector( log_level, logweights, RV ( particles.length()), prefix);
-        }
-        if (log_level[logmeans]) {
+                if (log_level[logweights]) {
+                        L.add_vector( log_level, logweights, RV ( particles.length()), prefix);
+                }
+                if (log_level[logNeff]) {
+                        L.add_vector( log_level, logNeff, RV ( 1), prefix);
+                }
+                if (log_level[logmeans]) {
             for (int i=0; i<particles.length(); i++) {
                 L.add_vector( log_level, logmeans, RV ( particles(i)->dimension() ), prefix , i);
@@ -548 +551 @@
         //shared_ptr<epdf> pred_v=bm->epredictor();
 
-        //vec vt=pred_v->sample();
-                vec vt = bm->samplepred();
+vec xt;
+vec xthat;
+//vec vt=pred_v->sample();
+//              vec vt = bm->samplepred();
+if (1){
+                //vec vt=pred_v->sample();
+                int dim_x = bm->dimensiony();
+                int dim_y = dimensiony();
+                
+                vec post_mean =bm->posterior().mean();
+                mat SigV;
+                if (post_mean.length()==dim_x){
+                        SigV=diag(vec(post_mean._data(), dim_x));
+                } else {
+                        SigV=mat(post_mean._data(), dim_x, dim_x);
+                }
+
+                vec &xtm=est_emp.point;
+                vec g_args(g->dimensionc());
+                x2g.filldown(xtm,g_args);
+                cond2g.filldown(cond,g_args);
+                xthat = g->eval(g_args);
+
+                mat IC = eye(dim_x+dim_y);
+                IC.set_submatrix(0,0,eye(dim_x));
+                mat SigJo=eye(dim_x + dim_y);
+                SigJo.set_submatrix(0,0,SigV);
+                mlnorm<ldmat>* mlfy=dynamic_cast<mlnorm<ldmat>* >(fy.get());
+                
+                SigJo.set_submatrix(dim_x,dim_x, mlfy->_R());
+                IC.set_submatrix(dim_x,0,-mlfy->_A());
+
+                mat iIC=inv(IC);
+                enorm<chmat> Jo;
+                Jo._mu()=iIC*concat(xthat, zeros(dim_y));
+                Jo._R()=iIC*SigJo*iIC.T();
+                Jo.set_rv(concat(bm->_yrv(),yrv));
+                Jo.validate();
+                shared_ptr<pdf> Cond=Jo.condition(bm->_yrv());
+                //new sample
+
+                xt = Cond->samplecond(dt);//  + vt;
+} else{
+
 
         //new sample
@@ -556 +601 @@
         x2g.filldown(xtm,g_args);
         cond2g.filldown(cond,g_args);
-        vec xt = g->eval(g_args) + vt;
+                xthat = g->eval(g_args);
+        
+                xt = xthat + bm->samplepred();
+                
+}
         est_emp.point=xt;
 
         // the vector [v_t] updates bm,
-        bm->bayes(vt);
+        bm->bayes(xt-xthat);
 
         // residue of observation
@@ -635 +684 @@
         */
 class NoiseParticleXY : public BM {
+protected:
         //! discrte particle
         dirac est_emp;
@@ -804 +854 @@
 UIREGISTER(NoiseParticleXY);
 
+class NoiseParticleXYprop: public NoiseParticleXY{
+public:
+        BM* _copy() const {
+                return new NoiseParticleXYprop(*this);
+        };
+        void bayes(const vec &dt, const vec &cond) {
+                int dim_x = bmx->dimensiony();
+                int dim_y = bmy->dimensiony();
+                
+                //vec vt=pred_v->sample();
+                mat SigV=mat(bmx->posterior().mean()._data(), dim_x, dim_x);
+                mat SigW=mat(bmy->posterior().mean()._data(), dim_y, dim_y);
+                
+                vec &xtm=est_emp.point;
+                vec g_args(g->dimensionc());
+                x2g.filldown(xtm,g_args);
+                cond2g.filldown(cond,g_args);
+                vec xthat = g->eval(g_args);
+                
+                mat IC = eye(dim_x+dim_y);
+                IC.set_submatrix(0,0,eye(dim_x));
+                mat SigJo=eye(dim_x + dim_y);
+                SigJo.set_submatrix(0,0,SigV);
+                SigJo.set_submatrix(dim_x,dim_x, 10*SigW);
+                
+                diffbifn* hb=dynamic_cast<bilinfn*>(h.get());
+                mat C=zeros(dim_y,dim_x);
+                if (hb){
+                        hb->dfdx_cond(xtm, cond, C,true);
+                        IC.set_submatrix(dim_x,0,-C);
+                }
+                
+                mat iIC=inv(IC);
+                enorm<chmat> Jo;
+                Jo._mu()=iIC*concat(xthat, zeros(dim_y));
+                Jo._R()=iIC*SigJo*iIC.T();
+                Jo.set_rv(concat(bmx->_yrv(), bmy->_yrv()));
+                Jo.validate();
+                shared_ptr<pdf> Cond=Jo.condition(bmx->_yrv());
+                //new sample
+                
+                vec xt = Cond->samplecond(dt);//  + vt;
+                
+                est_emp.point=xt;
+                
+                // the vector [v_t] updates bm,
+                bmx->bayes(xt-xthat);
+                
+                // residue of observation
+                vec h_args(h->dimensionc());
+                x2h.filldown(xt,h_args);
+                cond2h.filldown(cond,h_args);
+                
+                vec z_y =h->eval(h_args)-dt;
+                //                      ARX *abm = dynamic_cast<ARX*>(bmy.get());
+                /*                      double ll2;
+                 i f (abm){ //ARX*
+                 shared_ptr<epdf> pr_y(abm->epredictor_student(empty_vec));
+                 ll2=pr_y->evallog(z_y);
+        } else{
+                shared_ptr<epdf> pr_y(bmy->epredictor(empty_vec));
+                ll2=pr_y->evallog(z_y);
+        }*/
+                
+                bmy->bayes(z_y);
+                // test _lls
+                ll= bmy->_ll();
+        }
+};
+UIREGISTER(NoiseParticleXYprop);
+
 /*! Marginalized particle for state-space models with unknown parameters of residues distribution
 

library/bdm/math/functions.h

@@ -184 +184 @@
         bdm_assert_debug ( x0.length() == dimx, "bilinfn::eval Wrong xcond." );
         bdm_assert_debug ( u0.length() == dimu, "bilinfn::eval Wrong ucond." );
-        return A*x0 + B*u0;
+        return dimu>0 ?  A*x0 + B*u0 : A*x0;
     }
 
@@ -196 +196 @@
         if ( full ) F = B;      //else : nothing has changed no need to regenerate
     }
-    //!@}
-};
+    void from_setting(const Setting &set) {
+                diffbifn::from_setting(set);
+                UI::get(A,set,"A",UI::compulsory);
+                UI::get(B,set,"B",UI::compulsory);
+        }
+        void validate() {
+                dimy = A.rows();
+                dimx = A.cols();
+                dimu = B.cols();
+                dimc =dimx+dimu;
+        }
+        //!@}
+};
+UIREGISTER(bilinfn);
 
 } //namespace

library/bdm/stat/exp_family.h

@@ -384 +384 @@
 UIREGISTER2(estudent,ldmat);
 UIREGISTER2(estudent,chmat);
+
+// template < class sq_T, template <typename> class TEpdf = estudent >
+// class mstudent : public pdf_internal< TEpdf<sq_T> > {
+// protected:
+//      //! Internal epdf that arise by conditioning on \c rvc
+//      mat A;
+//      //! Constant additive term
+//      vec mu_const;
+//      
+// public:
+//      void condition( const vec &cond ) {
+//              iepdf._mu()=A*val*mu_const;
+//      }
+//              
+// };
 
 /*!
@@ -1136 +1151 @@
 class mlnorm : public pdf_internal< TEpdf<sq_T> > {
 protected:
-    //! Internal epdf that arise by conditioning on \c rvc
-    mat A;
-    //! Constant additive term
-    vec mu_const;
-//            vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
+        //! Internal epdf that arise by conditioning on \c rvc
+        mat A;
+        //! Constant additive term
+        vec mu_const;
+        //            vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
 public:
     //! \name Constructors

library/tests/testsuite/enorm.cfg

@@ -25 +25 @@
     };
   tolerance = 0.3;
+},
+{
+  class = "epdf_harness";
+  epdf = {
+    class = "enorm<ldmat>";
+    mu = [ 1.1, -2.0, 3.0];
+    R = ( "matrix", 3, 3, [ 1.0, -0.5, 0.8, -0.5, 2.0, 0.3, 0.8, 0.3, 3.0  ] );
+    rv : 
+    {
+      class = "RV";
+      names = ( "x", "y", "z" );
+    };
+  };
+  mean = [ 1.1, -2.0, 3.0 ];
+  variance = [ 1, 2, 3 ];
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ -5.0, 5.0 ], [-5.0, 5.0 ], [-1.0, 9.0] );
+        gridsizes = [ 50, 50, 50 ];
+        };
+ R = ( "matrix", 3, 3, [ 1.0, -0.5, 0.8, -0.5, 2.0, 0.3, 0.8, 0.3, 3.0  ] );
+  marginal_rv =     {
+      class = "RV";
+      names = ( "z" );
+    };
+  tolerance = 0.1;
 },
 {