Changeset 979

Timestamp:

05/25/10 20:55:06 (14 years ago)

Author:

smidl

Message:

estimator returns the array of posterior estimators as second argument

Files:

• applications/bdmtoolbox/mex/estimator.cpp (modified) (4 diffs)
• applications/bdmtoolbox/tutorial/userguide/arx_basic_example.m (modified) (1 diff)
• library/bdm/base/bdmbase.cpp (modified) (1 diff)
• library/bdm/estim/arx.cpp (modified) (2 diffs)
• library/bdm/estim/arx.h (modified) (2 diffs)
• library/bdm/estim/particles.h (modified) (6 diffs)
• library/bdm/stat/emix.h (modified) (2 diffs)

applications/bdmtoolbox/mex/estimator.cpp

@@ -38 +38 @@
 Execute command:
 \code
->> estimator('config_file.cfg');
+>> estimator
 \endcode
-when using loggers storing results on hard drives, and
-\code
->> Res=estimator('config_file.cfg');
-\endcode
-when using logger of the type \c "mex_logger". The results will be stored in structure \c M.
+and follow the help there.
 
  */
@@ -69 +65 @@
         // Check the number of inputs and output arguments
         if ( n_input<2 ) mexErrMsgTxt ( "Usage:\n"
-                                                "result=estimator(system, estimators, experiment, logger)\n"
+                                                "[Res,estimators,Res2]=estimator(system, estimators, experiment, logger)\n"
                                                 "  system     = struct('class','datasource',...);  % Estimated system\n"
                                                 "  estimators = {struct('class','estimator',...),  % Estimators\n"
@@ -76 +72 @@
                                                 "  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
                                                 "  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
+                                                "Output:\n"
+                                                "  Res          Matlab structure with logged results, \n"  
+                                                "  estimators   Array of estimators updated with data \n"
+                                                "  Res2         When logfull log_level is on, this structure is filled with structures of posterior densities\n\n"
                                                 "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM." );
 
@@ -206 +206 @@
                 if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
                 output[0] = mL->toCell();
-                if (n_output>1) {
+                if (n_output>1){ // write estimators
+                        UImxArray Ests;
+                        UI::save(Es, Ests,"estimators");
+                        output[1]=UImxArray::create_mxArray(Ests);
+                }
+                if (n_output>2) {
                         mL->_setting_conf().setAutoConvert(true);
-                        output[1]= UImxArray::create_mxArray(mL->_setting_conf().getRoot());
+                        output[2]= UImxArray::create_mxArray(mL->_setting_conf().getRoot());
                 }
         }
 #endif
         for (int i=0;i<Dls.length(); i++){delete Dls(i); delete Dlsc(i);}
+        UIFile F;
+        UI::save(Es, F.getRoot(),"estimators");
+        F.writeFile("estim_out.cfg");
 }

applications/bdmtoolbox/tutorial/userguide/arx_basic_example.m

@@ -15 +15 @@
 A1.log_level = 'logbounds,logevidence';
 
-M=estimator(DS,{A1});
+[M,Apost]=estimator(DS,{A1});
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

library/bdm/base/bdmbase.cpp

@@ -421 +421 @@
         root::log_register ( L, prefix );
 
+        if (dimension()==0) return;
+        
         RV r;
         if ( isnamed() ) {

library/bdm/estim/arx.cpp

@@ -223 +223 @@
         BMEF::from_setting(set);
         
-        RV rrv; 
-        UI::get (rrv, set, "rgr", UI::compulsory );
+        UI::get (rgr, set, "rgr", UI::compulsory );
         
         dimy = yrv._dsize();
         bdm_assert(dimy>0,"ARX::yrv should not be empty");
         // rgrlen - including constant!!!
-        rgrlen = rrv._dsize();
+        rgrlen = rgr._dsize();
 
         int constant;
@@ -238 +237 @@
         }
         dimc = rgrlen;
-        rvc = rrv;
+        rvc = rgr;
         rgrlen += int ( have_constant == true );
 

library/bdm/estim/arx.h

@@ -47 +47 @@
         //! vector of dyadic update
         vec dyad;
+        //! RV of regressor
+        RV rgr;
         //! length of the regressor
         int rgrlen;
@@ -165 +167 @@
         {                       
                 BMEF::to_setting( set ); // takes care of rv, yrv, rvc
+                UI::save(rgr, set, "rgr");
                 int constant = have_constant ? 1 : 0;
                 UI::save(constant, set, "constant");
-                UI::save(&est, set, "prior");
                 UI::save(&alter_est, set, "alternative");
-                
+                UI::save(&posterior(), set, "posterior");
                 
         } 

library/bdm/estim/particles.h

@@ -465 +465 @@
 UIREGISTER ( PF );
 
+/*! Marginalized particle for state-space models with unknown parameters of distribuution of residues on \f$v_t\f$. 
+
+\f{eqnarray*}{
+        x_t &=& g(x_{t-1}) + v_t,\\
+        y_t &\sim &fy(x_t),
+        \f}
+        
+        This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance.
+        
+        */
+class NoiseParticleX : public MarginalizedParticleBase{
+        protected:
+                //! function transforming xt, ut -> x_t+1
+                shared_ptr<fnc> g; // pdf for non-linear part
+                //! function transforming xt,ut -> yt
+                shared_ptr<pdf> fy; // pdf for non-linear part
+                
+                RV rvx;
+                RV rvxc;
+                RV rvyc;
+                
+                //!link from condition to f
+                datalink_part cond2g;
+                //!link from condition to h
+                datalink_part cond2fy;
+                //!link from xt to f
+                datalink_part x2g;
+                //!link from xt to h
+                datalink_part x2fy;
+        
+        public:
+                BM* _copy() const{return new NoiseParticleX(*this);};
+                void bayes(const vec &dt, const vec &cond){
+                        shared_ptr<epdf> pred_v=bm->epredictor();
+                        
+                        vec vt=pred_v->sample();
+                        
+                        //new sample
+                        vec &xtm=est_emp.point;
+                        vec g_args(g->dimensionc());
+                        x2g.filldown(xtm,g_args);
+                        cond2g.filldown(cond,g_args);
+                        vec xt = g->eval(g_args) + vt;
+                        est_emp.point=xt;
+
+                        // the vector [v_t] updates bm, 
+                        bm->bayes(vt);
+                        
+                        // residue of observation
+                        vec fy_args(fy->dimensionc());
+                        x2fy.filldown(xt,fy_args);
+                        cond2fy.filldown(cond,fy_args);
+                        
+                        ll=bm->_ll() + fy->evallogcond(dt,fy_args);
+                }
+                void from_setting(const Setting &set){
+                        MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc...
+                                
+                        g=UI::build<fnc>(set,"g",UI::compulsory);
+                        fy=UI::build<pdf>(set,"fy",UI::compulsory);
+                        UI::get(rvx,set,"rvx",UI::compulsory);
+                        est_emp.set_rv(rvx);
+                        
+                        UI::get(rvxc,set,"rvxc",UI::compulsory);
+                        UI::get(rvyc,set,"rvyc",UI::compulsory);
+                        
+                }
+                void validate(){
+                        MarginalizedParticleBase::validate();
+                        
+                        dimy = fy->dimension();
+                        bm->set_yrv(rvx);
+                        
+                        est_emp.set_rv(rvx);
+                        est_emp.set_dim(rvx._dsize());
+                        est.validate();
+                        // 
+                        //check dimensions
+                        rvc = rvxc.subt(rvx.copy_t(-1));
+                        rvc.add( rvyc);
+                        rvc=rvc.subt(rvx);
+                        
+                        bdm_assert(g->dimension()==rvx._dsize(),"rvx is not described");
+                        bdm_assert(g->dimensionc()==rvxc._dsize(),"rvxc is not described");
+                        bdm_assert(fy->dimensionc()==rvyc._dsize(),"rvyc is not described");
+                        
+                        bdm_assert(bm->dimensiony()==g->dimension(), 
+                                           "Incompatible noise estimator of dimension " +
+                                                num2str(bm->dimensiony()) + " does not match dimension of g , " + 
+                                                num2str(g->dimension()));
+                                                
+                        dimc = rvc._dsize();
+                                        
+                        //establish datalinks
+                        x2g.set_connection(rvxc, rvx.copy_t(-1));
+                        cond2g.set_connection(rvxc, rvc);
+                        
+                        x2fy.set_connection(rvyc, rvx);
+                        cond2fy.set_connection(rvyc, rvc);
+                }
+};
+UIREGISTER(NoiseParticleX);
+
 /*! Marginalized particle for state-space models with unknown parameters of residues distribution
 
-\f[
-\begin{eqnarray}
-x_t = g(x_{t-1}) + v_t,\\
-z_t = h(x_{t-1}) + w_t,\\
-\end{eqnarray}
-\f]
-
-This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance.
-
-*/
+\f{eqnarray*}{
+        x_t &=& g(x_{t-1}) + v_t,\\
+        z_t &= &h(x_{t-1}) + w_t,
+        \f}
+        
+        This particle is a only a shell creating the residues calling internal estimator of their parameters. The internal estimator can be of any compatible type, e.g. ARX for Gaussian residues with unknown mean and variance.
+        
+        */
 class NoiseParticle : public MarginalizedParticleBase{
         protected:
@@ -496 +597 @@
                 //!link from xt to h
                 datalink_part x2h;
-                
+        
         public:
                 BM* _copy() const{return new NoiseParticle(*this);};
@@ -524 +625 @@
                 void from_setting(const Setting &set){
                         MarginalizedParticleBase::from_setting(set); //reads bm, yrv,rvc, bm_rv, etc...
-                                
+                        
                         UI::get(g,set,"g",UI::compulsory);
                         UI::get(h,set,"h",UI::compulsory);
@@ -532 +633 @@
                         UI::get(rvxc,set,"rvxc",UI::compulsory);
                         UI::get(rvyc,set,"rvyc",UI::compulsory);
+                        
                 }
                 void validate(){
@@ -554 +656 @@
                         bdm_assert(bm->dimensiony()==g->dimension()+h->dimension(), 
                                            "Incompatible noise estimator of dimension " +
-                                                num2str(bm->dimensiony()) + " does not match dimension of g and h, " + 
-                                                num2str(g->dimension())+" and "+ num2str(h->dimension()) );
-                                                
-                        dimc = rvc._dsize();
-                                        
-                        //establish datalinks
-                        x2g.set_connection(rvxc, rvx.copy_t(-1));
-                        cond2g.set_connection(rvxc, rvc);
-                        
-                        x2h.set_connection(rvyc, rvx);
-                        cond2h.set_connection(rvyc, rvc);
+                                           num2str(bm->dimensiony()) + " does not match dimension of g and h, " + 
+                                           num2str(g->dimension())+" and "+ num2str(h->dimension()) );
+                                           
+                                           dimc = rvc._dsize();
+                                           
+                                           //establish datalinks
+                                           x2g.set_connection(rvxc, rvx.copy_t(-1));
+                                           cond2g.set_connection(rvxc, rvc);
+                                           
+                                           x2h.set_connection(rvyc, rvx);
+                                           cond2h.set_connection(rvyc, rvc);
                 }
 };
@@ -570 +672 @@
 
 
-
 }
 #endif // KF_H

library/bdm/stat/emix.h

@@ -278 +278 @@
         Array<datalink*> dls;
         //! interface for a factor
+public:
         virtual const epdf* factor(int i) const NOT_IMPLEMENTED(NULL); 
         //!number of factors
-        virtual const int no_factors() const =0;
-public:
+        virtual const int no_factors() const NOT_IMPLEMENTED(0);
         //! Default constructor
         eprod_base () :  dls (0) {};
@@ -394 +394 @@
 UIREGISTER ( mmix );
 
+
+//! Base class for all BM running as parallel update of internal BMs
+
+class ProdBMBase : public BM {
+        protected :
+                Array<vec_from_vec> bm_yt;
+                Array<vec_from_2vec> bm_cond;
+                class eprod_bm : public eprod_base {
+                        ProdBMBase & pb;
+                        public :
+                        eprod_bm(ProdBMBase &pb0): pb(pb0){}
+                        const epdf* factor(int i ) const {return &(pb.bm(i)->posterior());}
+                        const int no_factors() const {return pb.no_bms();}
+                } est;
+        public:
+                ProdBMBase():est(*this){}
+                virtual BM* bm(int i) NOT_IMPLEMENTED(NULL);
+                virtual int no_bms() const {return 0;}
+                const epdf& posterior() const {return est;}
+                void set_prior(const epdf *pri){
+                        const eprod_base* ep=dynamic_cast<const eprod_base*>(pri);
+                        if (ep){
+                                bdm_assert(ep->no_factors()!=no_bms() , "Given prior has "+ num2str(ep->no_factors()) + " while this ProdBM has "+
+                                        num2str(no_bms()) + "BMs");
+                                for (int i=0; i<no_bms(); i++){
+                                        bm(i)->set_prior(ep->factor(i));
+                                }
+                        }
+                }
+                
+                void validate() {
+                        est.validate();
+                        BM::validate();
+                        // set links
+                        bm_yt.set_length(no_bms());
+                        bm_cond.set_length(no_bms());
+                        
+                        // 
+                        
+                        for (int i=0; i<no_bms(); i++){
+                                yrv.add(bm(i)->_yrv());
+                                rvc.add(bm(i)->_rvc());
+                        }
+                        rvc=rvc.subt(yrv);
+                        
+                        dimy = yrv._dsize();
+                        dimc = rvc._dsize();
+                        
+                        for (int i=0; i<no_bms(); i++){
+                                bm_yt(i).connect(bm(i)->_yrv(), yrv);
+                                bm_cond(i).connect(bm(i)->_rvc(), yrv, rvc);
+                        }
+                }
+                void bayes(const vec &dt, const vec &cond){
+                        ll=0;
+                        for(int i=0;i<no_bms(); i++){
+                                bm_yt(i).update(dt);
+                                bm_cond(i).update(dt,cond);
+                                bm(i)->bayes(bm_yt(i), bm_cond(i));
+                        }
+                }
+                
+};
+
+class ProdBM: public ProdBMBase{
+        protected:
+                Array<shared_ptr<BM> > BMs;
+        public:
+                virtual BM* bm(int i) {return BMs(i).get();}
+                virtual int no_bms() const {return BMs.length();}
+                void from_setting(const Setting &set){
+                        BM::from_setting(set);
+                        UI::get(BMs,set,"BMs");
+                }
+                void to_setting(Setting &set) const{
+                        BM::to_setting(set);
+                        UI::save(BMs,set,"BMs");
+                }
+                
+};
+UIREGISTER(ProdBM);
+
 }
 #endif //MX_H