Changeset 679 for library

Timestamp:

10/23/09 00:05:25 (15 years ago)

Author:

smidl

Message:

Major changes in BM -- OK is only test suite and tests/tutorial -- the rest is broken!!!

Location:

library

Files:

• bdm/base/bdmbase.h (modified) (14 diffs)
• bdm/design/ctrlbase.cpp (modified) (1 diff)
• bdm/estim/arx.cpp (modified) (16 diffs)
• bdm/estim/arx.h (modified) (7 diffs)
• bdm/estim/kalman.cpp (modified) (10 diffs)
• bdm/estim/kalman.h (modified) (17 diffs)
• bdm/estim/mixtures.cpp (modified) (4 diffs)
• bdm/estim/mixtures.h (modified) (1 diff)
• bdm/estim/particles.cpp (modified) (3 diffs)
• bdm/estim/particles.h (modified) (8 diffs)
• bdm/itpp_ext.cpp (modified) (1 diff)
• bdm/itpp_ext.h (modified) (1 diff)
• bdm/mex/mex_BM.h (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (1 diff)
• bdm/stat/exp_family.h (modified) (4 diffs)
• bdm/stat/merger.cpp (modified) (1 diff)
• bdm/stat/merger.h (modified) (1 diff)
• tests/arx_elem_test.cpp (modified) (2 diffs)
• tests/arx_test.cpp (modified) (3 diffs)
• tests/dirfile-format.matrix (modified) (1 diff)
• tests/mixtures_test.cpp (modified) (2 diffs)
• tests/tutorial/arx_simple.cpp (modified) (1 diff)
• tests/tutorial/kalman_simple.cpp (modified) (1 diff)

library/bdm/base/bdmbase.h

@@ -406 +406 @@
 
         mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), dim( m.dim), rv( m.rv ) { }
+        
+        //! copy of the current object - make sure to implement
+        virtual mpdf* _copy_() const {return new mpdf(*this);}
         //!@}
 
@@ -421 +424 @@
 
         //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
-        virtual double evallogcond ( const vec &dt, const vec &cond ) {
+        virtual double evallogcond ( const vec &yt, const vec &cond ) {
                 bdm_error ( "Not implemented" );
                 return 0.0;
@@ -427 +430 @@
 
         //! Matrix version of evallogcond
-        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {
-                vec v ( Dt.cols() );
-                for ( int i = 0; i < Dt.cols(); i++ ) {
-                        v ( i ) = evallogcond ( Dt.get_col ( i ), cond );
+        virtual vec evallogcond_m ( const mat &Yt, const vec &cond ) {
+                vec v ( Yt.cols() );
+                for ( int i = 0; i < Yt.cols(); i++ ) {
+                        v ( i ) = evallogcond ( Yt.get_col ( i ), cond );
                 }
                 return v;
@@ -436 +439 @@
 
         //! Array<vec> version of evallogcond
-        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
+        virtual vec evallogcond_m ( const Array<vec> &Yt, const vec &cond ) {
                 bdm_error ( "Not implemented" );
                 return vec();
@@ -458 +461 @@
         }
 
+        //! access function
+        void set_dim(int d) {dim=d;}
+        //! access function
+        void set_dimc(int d) {dimc=d;}
         //! Load from structure with elements:
         //!  \code
@@ -506 +513 @@
                 dim = dim0;
         }
+        epdf* _copy_() const {return new epdf(*this);}
         //!@}
 
@@ -1020 +1028 @@
 This object represents exact or approximate evaluation of the Bayes rule:
 \f[
-f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
+f(\theta_t | y_1,\ldots,y_t, u_1,\ldots,u_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
 \f]
-
+where:
+ * \f$ y_t \f$ is the variable 
 Access to the resulting posterior density is via function \c posterior().
 
@@ -1028 +1037 @@
 It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor().
 
-Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$:
+Alternatively, it can evaluate posterior density with rvc replaced by the given values, \f$ c_t \f$:
 \f[
 f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1})
 \f]
 
-The value of \f$ c_t \f$ is set by function condition().
 
 */
@@ -1040 +1048 @@
 protected:
         //! Random variable of the data (optional)
-        RV drv;
+        RV yrv;
+        //! size of the data record
+        int dimy;
+        //! Name of extension variable
+        RV rvc;
+        //! size of the conditioning vector
+        int dimc;
+        
         //!Logarithm of marginalized data likelihood.
         double ll;
         //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
         bool evalll;
-        
+
 public:
         //! \name Constructors
         //! @{
 
-        BM() : ll ( 0 ), evalll ( true ) { };
-        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
+        BM() : yrv(),dimy(0),rvc(),dimc(0), ll ( 0 ), evalll ( true ) { };
+        BM ( const BM &B ) :  yrv ( B.yrv ), dimy(B.dimy), rvc ( B.rvc ), ll ( B.ll ), evalll ( B.evalll ) {}
         //! \brief Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
         //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
-        virtual BM* _copy_() const {
-                return NULL;
-        };
+        virtual BM* _copy_() const { return NULL; };
         //!@}
 
@@ -1065 +1078 @@
         @param dt vector of input data
         */
-        virtual void bayes ( const vec &dt ) = 0;
+        virtual void bayes ( const vec &yt, const vec &cond=empty_vec ) = 0;
         //! Batch Bayes rule (columns of Dt are observations)
         virtual void bayesB ( const mat &Dt );
         //! Evaluates predictive log-likelihood of the given data record
-        //! I.e. marginal likelihood of the data with the posterior integrated out.
-        virtual double logpred ( const vec &dt ) const {
+        //! I.e. marginal likelihood of the data with the posterior integrated out. 
+        //! This function evaluates only \f$ y_t \f$, condition is assumed to be the last used in bayes(). 
+        //! See bdm::BM::predictor for conditional version.
+        virtual double logpred ( const vec &yt ) const {
                 bdm_error ( "Not implemented" );
                 return 0.0;
@@ -1076 +1091 @@
 
         //! Matrix version of logpred
-        vec logpred_m ( const mat &dt ) const {
-                vec tmp ( dt.cols() );
-                for ( int i = 0; i < dt.cols(); i++ ) {
-                        tmp ( i ) = logpred ( dt.get_col ( i ) );
+        vec logpred_m ( const mat &Yt ) const {
+                vec tmp ( Yt.cols() );
+                for ( int i = 0; i < Yt.cols(); i++ ) {
+                        tmp ( i ) = logpred ( Yt.get_col ( i ) );
                 }
                 return tmp;
@@ -1096 +1111 @@
         //!@}
 
-        //! \name Extension to conditional BM
-        //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
-        //! Alternatively, it can be used for automated connection to DS when the condition is observed
-        //!@{
-
-        //! Name of extension variable
-        RV rvc;
-        //! access function
-        const RV& _rvc() const {
-                return rvc;
-        }
-
-        //! Substitute \c val for \c rvc.
-        virtual void condition ( const vec &val ) {
-                bdm_error ( "Not implemented!" );
-        }
-
-        //!@}
-
 
         //! \name Access to attributes
         //!@{
+                //! access function
+                const RV& _rvc() const {
+                        return rvc;
+                }
+                //! access function
+                int dimensionc() const {
+                        return dimc;
+                }
+                //! access function
+                int dimensiony() const {
+                        return dimy;
+                }
+                //! access function
+                int dimension() const {
+                        return posterior().dimension();
+                }
+                //! access function
+        const RV& _yrv() const {
+                return yrv;
+        }
         //! access function
-        const RV& _drv() const {
-                return drv;
-        }
-        //! access function
-        void set_drv ( const RV &rv ) {
-                drv = rv;
+        void set_yrv ( const RV &rv ) {
+                yrv = rv;
         }
         //! access to rv of the posterior
         void set_rv ( const RV &rv ) {
-                const_cast<epdf&> ( posterior() ).set_rv ( rv );
+                const_cast<epdf&>(posterior()).set_rv ( rv );
+        }
+        //! access function
+        void set_dim ( int dim ) {
+                const_cast<epdf&>(posterior()).set_dim ( dim );
         }
         //! return internal log-likelihood of the last data vector
@@ -1180 +1196 @@
                 shared_ptr<RV> r = UI::build<RV> ( set, "drv", UI::optional );
                 if ( r ) {
-                        set_drv ( *r );
+                        set_yrv ( *r );
                 }
                 string opt;
@@ -1215 +1231 @@
 
 template<class EPDF>
-double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
+double mpdf_internal<EPDF>::evallogcond ( const vec &yt, const vec &cond ) {
         double tmp;
         condition ( cond );
-        tmp = iepdf.evallog ( dt );
+        tmp = iepdf.evallog ( yt );
         return tmp;
 }
 
 template<class EPDF>
-vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
+vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Yt, const vec &cond ) {
         condition ( cond );
-        return iepdf.evallog_m ( Dt );
+        return iepdf.evallog_m ( Yt );
 }
 
 template<class EPDF>
-vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
+vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Yt, const vec &cond ) {
         condition ( cond );
-        return iepdf.evallog_m ( Dt );
+        return iepdf.evallog_m ( Yt );
 }
 

library/bdm/design/ctrlbase.cpp

@@ -5 +5 @@
 void LQG::set_system(shared_ptr<StateSpace<fsqmat> > S0){
         S = S0;
-        dimx= S->_dimx();
-        dimy= S->_dimy();
-        dimu= S->_dimu();
+        dimx= S->_A().rows();
+        dimy= S->_C().rows();
+        dimu= S->_B().cols();
         pr=zeros(dimx+dimu+dimy,  dimu+dimx+dimu+dimy);
         pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy));

library/bdm/estim/arx.cpp

@@ -2 +2 @@
 namespace bdm {
 
-void ARX::bayes ( const vec &dt, const double w ) {
+        void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
         double lnc;
-                
+        //cache
+        ldmat &V=est._V(); 
+        double &nu=est._nu();
+
+        dyad.set_subvector(0,yt);
+        dyad.set_subvector(dimy,cond);
+        // possible "1" is there from the beginning
+        
         if ( frg < 1.0 ) {
                 est.pow ( frg ); // multiply V and nu
+
                 
                 //stabilize
                 ldmat V0=alter_est._V(); //$ copy
                 double &nu0=alter_est._nu();
+                
                 V0*=(1-frg);
                 V += V0; //stabilization
@@ -20 +29 @@
                 }
         }
-        if (have_constant) {
-                _dt.set_subvector(0,dt);
-                V.opupdt ( _dt, w );
-        } else {
-                V.opupdt ( dt, w );
-        }
+        V.opupdt ( dyad, w );
         nu += w;
 
@@ -36 +40 @@
 }
 
-double ARX::logpred ( const vec &dt ) const {
+double ARX::logpred ( const vec &yt ) const {
         egiw pred ( est );
         ldmat &V = pred._V();
@@ -42 +46 @@
 
         double lll;
+        vec dyad_p = dyad;
+        dyad_p.set_subvector(0,yt);
 
         if ( frg < 1.0 ) {
@@ -53 +59 @@
                 }
 
-        V.opupdt ( dt, 1.0 );
+        V.opupdt ( dyad_p, 1.0 );
         nu += 1.0;
         // log(sqrt(2*pi)) = 0.91893853320467
@@ -67 +73 @@
         const ARX* A0 = dynamic_cast<const ARX*> ( B0 );
 
-        bdm_assert_debug ( V.rows() == A0->V.rows(), "ARX::set_statistics Statistics  differ" );
-        set_statistics ( A0->dimx, A0->V, A0->nu );
+        bdm_assert_debug ( dimension() == A0->dimension(), "Statistics of different dimensions" );
+        set_statistics ( A0->dimensiony(), A0->posterior()._V(), A0->posterior()._nu() );
 }
 
 enorm<ldmat>* ARX::epredictor ( const vec &rgr ) const {
-        int dim = dimx;//est.dimension();
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        mat mu ( dimy, posterior()._V().rows() - dimy );
+        mat R ( dimy, dimy );
 
         enorm<ldmat>* tmp;
         tmp = new enorm<ldmat> ( );
         //TODO: too hackish
-        if ( drv._dsize() > 0 ) {
+        if ( yrv._dsize() > 0 ) {
         }
 
@@ -91 +96 @@
 
 mlnorm<ldmat>* ARX::predictor ( ) const {
-        int dim = est.dimension();
-        
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        mat mu ( dimy, posterior()._V().rows() - dimy );
+        mat R ( dimy, dimy );
         mlnorm<ldmat>* tmp;
         tmp = new mlnorm<ldmat> ( );
@@ -106 +109 @@
                 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );
         } else {
-                tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );
+                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) );
         }
         return tmp;
@@ -112 +115 @@
 
 mlstudent* ARX::predictor_student ( ) const {
-        int dim = est.dimension();
-
-        mat mu ( dim, V.rows() - dim );
-        mat R ( dim, dim );
+        const ldmat &V = posterior()._V();
+        
+        mat mu ( dimy, V.rows() - dimy );
+        mat R ( dimy, dimy );
         mlstudent* tmp;
         tmp = new mlstudent ( );
@@ -122 +125 @@
         mu = mu.T();
 
-        int xdim = dimx;
         int end = V._L().rows() - 1;
-        ldmat Lam ( V._L() ( xdim, end, xdim, end ), V._D() ( xdim, end ) );  //exp val of R
+        ldmat Lam ( V._L() ( dimy, end, dimy, end ), V._D() ( dimy, end ) );  //exp val of R
 
 
@@ -132 +134 @@
                         tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
                 } else {
-                        tmp->set_parameters ( mat ( dim, 0 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
+                        tmp->set_parameters ( mat ( dimy, dimc ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
                 }
         } else {
                 // no constant term
-                tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam );
+                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ), Lam );
         }
         return tmp;
@@ -214 +216 @@
 
 void ARX::from_setting ( const Setting &set ) {
-        shared_ptr<RV> yrv = UI::build<RV> ( set, "rv", UI::compulsory );
+        shared_ptr<RV> yrv_ = UI::build<RV> ( set, "rv", UI::compulsory );
         shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );
-        int ylen = yrv->_dsize();
+        dimy = yrv_->_dsize();
         // rgrlen - including constant!!!
-        int rgrlen = rrv->_dsize();
-        
-        dimx = ylen;
-        set_rv ( *yrv, *rrv );
+        dimc = rrv->_dsize();
+        
+        yrv = *yrv_;
+        rvc = *rrv;
         
         string opt;
@@ -233 +235 @@
                 have_constant=constant>0;
         }
-        if (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);}
+        int rgrlen = dimc+int(have_constant==true);
 
         //init
         shared_ptr<egiw> pri=UI::build<egiw>(set, "prior", UI::optional);
         if (pri) {
-                bdm_assert(pri->_dimx()==ylen,"prior is not compatible");
-                bdm_assert(pri->_V().rows()==ylen+rgrlen,"prior is not compatible");
+                bdm_assert(pri->_dimx()==dimy,"prior is not compatible");
+                bdm_assert(pri->_V().rows()==dimy+rgrlen,"prior is not compatible");
                 est.set_parameters( pri->_dimx(),pri->_V(), pri->_nu());
         }else{
-                est.set_parameters( ylen, zeros(ylen+rgrlen));
+                est.set_parameters( dimy, zeros(dimy+rgrlen));
                 set_prior_default(est);
         }
@@ -248 +250 @@
         shared_ptr<egiw> alt=UI::build<egiw>(set, "alternative", UI::optional);
         if (alt) {
-                bdm_assert(alt->_dimx()==ylen,"alternative is not compatible");
-                bdm_assert(alt->_V().rows()==ylen+rgrlen,"alternative is not compatible");
+                bdm_assert(alt->_dimx()==dimy,"alternative is not compatible");
+                bdm_assert(alt->_V().rows()==dimy+rgrlen,"alternative is not compatible");
                 alter_est.set_parameters( alt->_dimx(),alt->_V(), alt->_nu());
         } else {
@@ -264 +266 @@
         shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional );
         if (!rv_par){
-                est.set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) );
+                est.set_rv ( RV ( "{theta r }", vec_2 ( dimy*rgrlen, dimy*dimy ) ) );
         } else {
                 est.set_rv ( *rv_par );
@@ -270 +272 @@
         validate();
 }
-
-}
+}
+

library/bdm/estim/arx.h

@@ -43 +43 @@
 class ARX: public BMEF {
 protected:
-        //!size of output variable (needed in regressors)
-        int dimx;
-        //!description of modelled data \f$ y_t \f$ in the likelihood function
-        //! Do NOT access directly, only via \c get_yrv().
-        RV _yrv;
-        //! rv of regressor
-        RV rgrrv;
-        //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density
-        egiw est;
-        //! cached value of est.V
-        ldmat &V;
-        //! cached value of est.nu
-        double ν
         //! switch if constant is modelled or not
         bool have_constant;
-        //! cached value of data vector for have_constant =true
-        vec _dt;
+        //! vector of dyadic update
+        vec dyad;
+        //! posterior density
+        egiw est;
         //! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy]
         egiw alter_est;
@@ -65 +54 @@
         //! \name Constructors
         //!@{
-        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), est (), V ( est._V() ), nu ( est._nu() ), have_constant(true), _dt() {};
-        ARX ( const ARX &A0 ) : BMEF (A0.frg), est (A0.est), V ( est._V() ), nu ( est._nu() ), have_constant(A0.have_constant), _dt(A0._dt) {
-                dimx = A0.dimx;
-                _yrv = A0._yrv;
-                rgrrv = A0.rgrrv;
-                set_drv(A0._drv());
-        };
+        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant(true), dyad(), est() {};
+        ARX ( const ARX &A0 ) : BMEF (A0.frg),  have_constant(A0.have_constant), dyad(A0.dyad),est(est) { };
         ARX* _copy_() const;
         void set_parameters ( double frg0 ) {
@@ -79 +63 @@
                 have_constant=const0;
         }
-        void set_statistics ( int dimx0, const ldmat V0, double nu0 = -1.0 ) {
-                est.set_parameters ( dimx0, V0, nu0 );
+        void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) {
+                est.set_parameters ( dimy0, V0, nu0 );
                 last_lognc = est.lognc();
-                dimx = dimx0;
+                dimy = dimy0;
         }
         //!@}
@@ -93 +77 @@
 
         //! Weighted Bayes \f$ dt = [y_t psi_t] \f$.
-        void bayes ( const vec &dt, const double w );
-        void bayes ( const vec &dt ) {
-                bayes ( dt, 1.0 );
+        void bayes_weighted ( const vec &yt, const vec &cond=empty_vec, const double w=1.0 );
+        void bayes( const vec &yt, const vec &cond=empty_vec ) {
+                bayes_weighted ( yt,cond, 1.0 );
         };
-        double logpred ( const vec &dt ) const;
+        double logpred ( const vec &yt ) const;
         void flatten ( const BMEF* B ) {
                 const ARX* A = dynamic_cast<const ARX*> ( B );
                 // nu should be equal to B.nu
-                est.pow ( A->nu / nu );
+                est.pow ( A->posterior()._nu() / posterior()._nu() );
                 if ( evalll ) {
                         last_lognc = est.lognc();
@@ -110 +94 @@
         //! Predictor for empty regressor
         enorm<ldmat>* epredictor() const {
-                bdm_assert_debug ( dimx == V.rows() - 1, "Regressor is not only 1" );
+                bdm_assert_debug ( dimy == posterior()._V().rows() - 1, "Regressor is not only 1" );
                 return epredictor ( vec_1 ( 1.0 ) );
         }
@@ -127 +111 @@
                 const egiw& posterior() const {
                 return est;
-        }
-        //!@}
-
-        //!\name Connection
-        //!@{
-        void set_rv ( const RV &yrv0 , const RV &rgrrv0 ) {
-                _yrv = yrv0;
-                rgrrv=rgrrv0;
-                set_drv(concat(yrv0, rgrrv));
-        }
-
-        RV& get_yrv() {
-                //if yrv is not ready create it
-                if ( _yrv._dsize() != dimx ) {
-                        int i = 0;
-                        while ( _yrv._dsize() < dimx ) {
-                                _yrv.add ( drv ( vec_1 ( i ) ) );
-                                i++;
-                        }
-                }
-                //yrv should be ready by now
-                bdm_assert_debug ( _yrv._dsize() == dimx, "incompatible drv" );
-                return _yrv;
         }
         //!@}
@@ -174 +135 @@
 
         void validate() {
-                bdm_assert(dimx == _yrv._dsize(), "RVs of parameters and regressor do not match");
-                
+                //if dimc not set set it from V
+                if (dimc==0){
+                        dimc = posterior()._V().rows()-dimy-int(have_constant==true);
+                }
+                        
+                if (have_constant) {
+                        dyad = ones(dimy+dimc+1);
+                } else { 
+                        dyad = zeros(dimy+dimc);
+                }
+                        
         }
         //! function sets prior and alternative density 

library/bdm/estim/kalman.cpp

@@ -8 +8 @@
 
 
-void KalmanFull::bayes ( const vec &dt ) {
-        bdm_assert_debug ( dt.length() == ( dimy + dimu ), "KalmanFull::bayes wrong size of dt" );
-
-        vec u = dt.get ( dimy, dimy + dimu - 1 );
-        vec y = dt.get ( 0, dimy - 1 );
-        vec& mu = est->_mu();
-        mat &P = est->_R();
+void KalmanFull::bayes ( const vec &yt, const vec &cond ) {
+        bdm_assert_debug ( yt.length() == ( dimy ), "KalmanFull::bayes wrong size of dt" );
+
+        const vec &u = cond; // in this case cond=ut
+        const vec &y = yt;
+        
+        vec& mu = est._mu();
+        mat &P = est._R();
         mat& _Ry = fy._R();
         vec& yp = fy._mu();
@@ -42 +43 @@
         phxu = phxu0;
 
-        dimx = pfxu0->_dimx();
+        set_dim( pfxu0->_dimx());
         dimy = phxu0->dimension();
-        dimu = pfxu0->_dimu();
-        est->set_parameters( zeros(dimx), eye(dimx) );
-
-        A.set_size ( dimx, dimx );
-        C.set_size ( dimy, dimx );
+        dimc = pfxu0->_dimu();
+        est.set_parameters( zeros(dimension()), eye(dimension()) );
+
+        A.set_size ( dimension(), dimension() );
+        C.set_size ( dimy, dimension() );
         //initialize matrices A C, later, these will be only updated!
-        pfxu->dfdx_cond ( est->_mu(), zeros ( dimu ), A, true );
+        pfxu->dfdx_cond ( est._mu(), zeros ( dimc ), A, true );
         B.clear();
-        phxu->dfdx_cond ( est->_mu(), zeros ( dimu ), C, true );
+        phxu->dfdx_cond ( est._mu(), zeros ( dimc ), C, true );
         D.clear();
 
@@ -59 +60 @@
 }
 
-void EKFfull::bayes ( const vec &dt ) {
-        bdm_assert_debug ( dt.length() == ( dimy + dimu ), "EKFull::bayes wrong size of dt" );
-
-        vec u = dt.get ( dimy, dimy + dimu - 1 );
-        vec y = dt.get ( 0, dimy - 1 );
-        vec &mu = est->_mu();
-        mat &P = est->_R();
+void EKFfull::bayes ( const vec &yt, const vec &cond) {
+        bdm_assert_debug ( yt.length() == ( dimy ), "EKFull::bayes wrong size of dt" );
+        bdm_assert_debug ( cond.length() == ( dimc ), "EKFull::bayes wrong size of dt" );
+        
+        const vec &u = cond;
+        const vec &y = yt; //lazy to change it
+        vec &mu = est._mu();
+        mat &P = est._R();
         mat& _Ry = fy._R();
         vec& yp = fy._mu();
         
-        pfxu->dfdx_cond ( mu, zeros ( dimu ), A, true );
-        phxu->dfdx_cond ( mu, zeros ( dimu ), C, true );
+        pfxu->dfdx_cond ( mu, zeros ( dimc ), A, true );
+        phxu->dfdx_cond ( mu, zeros ( dimc ), C, true );
 
         //Time update
@@ -94 +96 @@
         ( ( StateSpace<chmat>* ) this )->set_parameters ( A0, B0, C0, D0, Q0, R0 );
         
-        _K=zeros(dimx,dimy);
-        // Cholesky special!
-        initialize();
+        _K=zeros(dimension(),dimy);
 }
 
 void KalmanCh::initialize(){
-        preA = zeros ( dimy + dimx + dimx, dimy + dimx );
+        preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() );
 //      preA.clear();
         preA.set_submatrix ( 0, 0, R._Ch() );
-        preA.set_submatrix ( dimy + dimx, dimy, Q._Ch() );
-}
-
-void KalmanCh::bayes ( const vec &dt ) {
-
-        vec u = dt.get ( dimy, dimy + dimu - 1 );
-        vec y = dt.get ( 0, dimy - 1 );
+        preA.set_submatrix ( dimy + dimension(), dimy, Q._Ch() );
+}
+
+void KalmanCh::bayes ( const vec &yt, const vec &cond ) {
+        bdm_assert_debug(yt.length()==dimy, "yt mismatch");
+        bdm_assert_debug(cond.length()==dimc, "yt mismatch");
+        
+        const vec &u = cond;
+        const vec &y = yt;
         vec pom ( dimy );
 
-        chmat &_P=est->_R();
-        vec &_mu = est->_mu();
-        mat _K(dimx,dimy);
+        chmat &_P=est._R();
+        vec &_mu = est._mu();
+        mat _K(dimension(),dimy);
         chmat &_Ry=fy._R();
         vec &_yp = fy._mu();
@@ -130 +132 @@
 
         ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 );
-        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimx - 1 ) ).T();
-        ( _P._Ch() ) = postA ( dimy, dimy + dimx - 1, dimy, dimy + dimx - 1 );
+        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimension() - 1 ) ).T();
+        ( _P._Ch() ) = postA ( dimy, dimy + dimension() - 1, dimy, dimy + dimension() - 1 );
 
         _mu = A * ( _mu ) + B * u;
@@ -154 +156 @@
         phxu = phxu0;
 
-        dimx = pfxu0->dimension();
+        set_dim( pfxu0->_dimx());
         dimy = phxu0->dimension();
-        dimu = pfxu0->_dimu();
-        
-        vec &_mu = est->_mu();
+        dimc = pfxu0->_dimu();
+        
+        vec &_mu = est._mu();
         // if mu is not set, set it to zeros, just for constant terms of A and C
-        if ( _mu.length() != dimx ) _mu = zeros ( dimx ); 
-        A = zeros ( dimx, dimx );
-        C = zeros ( dimy, dimx );
-        preA = zeros ( dimy + dimx + dimx, dimy + dimx );
+        if ( _mu.length() != dimension() ) _mu = zeros ( dimension() ); 
+        A = zeros ( dimension(), dimension() );
+        C = zeros ( dimy, dimension() );
+        preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() );
 
         //initialize matrices A C, later, these will be only updated!
-        pfxu->dfdx_cond ( _mu, zeros ( dimu ), A, true );
+        pfxu->dfdx_cond ( _mu, zeros ( dimc ), A, true );
 //      pfxu->dfdu_cond ( *_mu,zeros ( dimu ),B,true );
         B.clear();
-        phxu->dfdx_cond ( _mu, zeros ( dimu ), C, true );
+        phxu->dfdx_cond ( _mu, zeros ( dimc ), C, true );
 //      phxu->dfdu_cond ( *_mu,zeros ( dimu ),D,true );
         D.clear();
@@ -179 +181 @@
         preA.clear();
         preA.set_submatrix ( 0, 0, R._Ch() );
-        preA.set_submatrix ( dimy + dimx, dimy, Q._Ch() );
-}
-
-
-void EKFCh::bayes ( const vec &dt ) {
+        preA.set_submatrix ( dimy + dimension(), dimy, Q._Ch() );
+}
+
+
+void EKFCh::bayes ( const vec &yt, const vec &cond ) {
 
         vec pom ( dimy );
-        vec u = dt.get ( dimy, dimy + dimu - 1 );
-        vec y = dt.get ( 0, dimy - 1 );
-        vec &_mu = est->_mu();
-        chmat &_P = est->_R();
+        const vec &u = cond;
+        const vec &y = yt;
+        vec &_mu = est._mu();
+        chmat &_P = est._R();
         chmat &_Ry = fy._R();
         vec &_yp = fy._mu();
@@ -211 +213 @@
 
         ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 );
-        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimx - 1 ) ).T();
-        ( _P._Ch() ) = postA ( dimy, dimy + dimx - 1, dimy, dimy + dimx - 1 );
+        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimension() - 1 ) ).T();
+        ( _P._Ch() ) = postA ( dimy, dimy + dimension() - 1, dimy, dimy + dimension() - 1 );
 
 //      mat iRY = inv(_Ry->to_mat());
@@ -268 +270 @@
         //connect
         shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory );
-        set_drv ( *drv );
+        set_yrv ( *drv );
         shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory );
         set_rv ( *rv );
@@ -282 +284 @@
 
         set_parameters ( A );
-        set_drv ( A ( 0 )->_drv() );
+        set_yrv ( A ( 0 )->_yrv() );
         //set_rv(A(0)->_rv());
 

library/bdm/estim/kalman.h

@@ -34 +34 @@
 {
         protected:
-                //! cache of rv.count()
-                int dimx;
-                //! cache of rvy.count()
-                int dimy;
-                //! cache of rvu.count()
-                int dimu;
                 //! Matrix A
                 mat A;
@@ -53 +47 @@
                 sq_T R;
         public:
-                StateSpace() : dimx (0), dimy (0), dimu (0), A(), B(), C(), D(), Q(), R() {}
+                StateSpace() :  A(), B(), C(), D(), Q(), R() {}
                 //!copy constructor
-                StateSpace(const StateSpace<sq_T> &S0) : dimx (S0.dimx), dimy (S0.dimy), dimu (S0.dimu), A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {}
+                StateSpace(const StateSpace<sq_T> &S0) :  A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {}
                 //! set all matrix parameters
                 void set_parameters (const mat &A0, const  mat &B0, const  mat &C0, const  mat &D0, const  sq_T &Q0, const sq_T &R0);
@@ -83 +77 @@
                 }               
                 //! access function
-                int _dimx(){return dimx;}
-                //! access function
-                int _dimy(){return dimy;}
-                //! access function
-                int _dimu(){return dimu;}
-                //! access function
                 const mat& _A() const {return A;}
                 //! access function
@@ -109 +97 @@
                 //! id of output
                 RV yrv;
-                //! id of input
-                RV urv;
                 //! Kalman gain
                 mat  _K;
                 //!posterior
-                shared_ptr<enorm<sq_T> > est;
+                enorm<sq_T> est;
                 //!marginal on data f(y|y)
                 enorm<sq_T>  fy;
         public:
-                Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(),urv(), _K(),  est(new enorm<sq_T>){}
+                Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(), _K(),  est(){}
                 //! Copy constructor
-                Kalman<sq_T>(const Kalman<sq_T> &K0) : BM(K0), StateSpace<sq_T>(K0), yrv(K0.yrv),urv(K0.urv), _K(K0._K),  est(new enorm<sq_T>(*K0.est)), fy(K0.fy){}
+                Kalman<sq_T>(const Kalman<sq_T> &K0) : BM(K0), StateSpace<sq_T>(K0), yrv(K0.yrv), _K(K0._K),  est(K0.est), fy(K0.fy){}
                 //!set statistics of the posterior
-                void set_statistics (const vec &mu0, const mat &P0) {est->set_parameters (mu0, P0); };
+                void set_statistics (const vec &mu0, const mat &P0) {est.set_parameters (mu0, P0); };
                 //!set statistics of the posterior
-                void set_statistics (const vec &mu0, const sq_T &P0) {est->set_parameters (mu0, P0); };
+                void set_statistics (const vec &mu0, const sq_T &P0) {est.set_parameters (mu0, P0); };
                 //! return correctly typed posterior (covariant return)
-                const enorm<sq_T>& posterior() const {return *est.get();}
-                //! shared posterior
-                shared_ptr<epdf> shared_posterior() {return est;}
+                const enorm<sq_T>& posterior() const {return est;}
                 //! load basic elements of Kalman from structure
                 void from_setting (const Setting &set) {
@@ -139 +123 @@
                         // Initial values
                         UI::get (yrv, set, "yrv", UI::optional);
-                        UI::get (urv, set, "urv", UI::optional);
-                        set_drv(concat(yrv,urv));
+                        UI::get (rvc, set, "urv", UI::optional);
+                        set_yrv(concat(yrv,rvc));
                         
                         validate();
@@ -147 +131 @@
                 void validate() {
                         StateSpace<sq_T>::validate();
-                        bdm_assert(est->dimension(), "Statistics and model parameters mismatch");
+                        dimy = this->C.rows();
+                        dimc = this->B.cols();
+                        set_dim(this->A.rows());
+                        
+                        bdm_assert(est.dimension(), "Statistics and model parameters mismatch");
                 }
 };
@@ -160 +148 @@
                 KalmanFull() :Kalman<fsqmat>(){};
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
                 BM* _copy_() const {
                         KalmanFull* K = new KalmanFull;
                         K->set_parameters (A, B, C, D, Q, R);
-                        K->set_statistics (est->_mu(), est->_R());
+                        K->set_statistics (est._mu(), est._R());
                         return K;
                 }
@@ -192 +180 @@
                         KalmanCh* K = new KalmanCh;
                         K->set_parameters (A, B, C, D, Q, R);
-                        K->set_statistics (est->_mu(), est->_R());
+                        K->set_statistics (est._mu(), est._R());
                         return K;
                 }
@@ -213 +201 @@
                 Thus this object evaluates only predictors! Not filtering densities.
                 */
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
 
                 void from_setting(const Setting &set){
                         Kalman<chmat>::from_setting(set);
+                        validate();
+                }
+                void validate() {
+                        Kalman<chmat>::validate();
                         initialize();
                 }
@@ -244 +236 @@
 
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
                 //! set estimates
                 void set_statistics (const vec &mu0, const mat &P0) {
-                        est->set_parameters (mu0, P0);
+                        est.set_parameters (mu0, P0);
                 };
                 //! access function
                 const mat _R() {
-                        return est->_R().to_mat();
+                        return est._R().to_mat();
                 }
                 void from_setting (const Setting &set) {
@@ -280 +272 @@
                         //connect
                         shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory );
-                        set_drv ( *drv );
+                        set_yrv ( *drv );
                         shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory );
                         set_rv ( *rv );
@@ -337 +329 @@
                         EKFCh* E = new EKFCh;
                         E->set_parameters (pfxu, phxu, Q, R);
-                        E->set_statistics (est->_mu(), est->_R());
+                        E->set_statistics (est._mu(), est._R());
                         return E;
                 }
@@ -344 +336 @@
 
                 //! Here dt = [yt;ut] of appropriate dimensions
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
 
                 void from_setting (const Setting &set);
@@ -389 +381 @@
                         est.set_parameters (A (0)->posterior().mean(), A (0)->posterior()._R());
                 }
-                void bayes (const vec &dt) {
+                void bayes (const vec &yt, const vec &cond=empty_vec) {
                         int n = Models.length();
                         int i;
                         for (i = 0; i < n; i++) {
-                                Models (i)->bayes (dt);
+                                Models (i)->bayes (yt);
                                 _lls (i) = Models (i)->_ll();
                         }
@@ -475 +467 @@
                                 Crv.add(urv.copy_t(t));
                         }
-                        
-                        this->dimx = xrv._dsize();
-                        this->dimy = yrv._dsize();
-                        this->dimu = urv._dsize();
-                        
+                                                
                         // get mapp
                         th2A.set_connection(xrv, ml._rvc());
@@ -486 +474 @@
 
                         //set matrix sizes
-                        this->A=zeros(dimx,dimx);
-                        for (int j=1; j<dimx; j++){A(j,j-1)=1.0;} // off diagonal
-                                this->B=zeros(dimx,1);
+                        this->A=zeros(xrv._dsize(),xrv._dsize());
+                        for (int j=1; j<xrv._dsize(); j++){A(j,j-1)=1.0;} // off diagonal
+                                this->B=zeros(xrv._dsize(),1);
                                 this->B(0) = 1.0;
-                                this->C=zeros(1,dimx);
+                                this->C=zeros(1,xrv._dsize());
                                 this->D=zeros(1,urv._dsize());
-                                this->Q = zeros(dimx,dimx);
+                                this->Q = zeros(xrv._dsize(),xrv._dsize());
                         // R is set by update
                         
@@ -538 +526 @@
 template<class sq_T>
 void StateSpace<sq_T>::validate(){
-        dimx = A.rows();
-        dimu = B.cols();
-        dimy = C.rows();
-        bdm_assert (A.cols() == dimx, "KalmanFull: A is not square");
-        bdm_assert (B.rows() == dimx, "KalmanFull: B is not compatible");
-        bdm_assert (C.cols() == dimx, "KalmanFull: C is not square");
-        bdm_assert ( (D.rows() == dimy) && (D.cols() == dimu), "KalmanFull: D is not compatible");
-        bdm_assert ( (Q.cols() == dimx) && (Q.rows() == dimx), "KalmanFull: Q is not compatible");
-        bdm_assert ( (R.cols() == dimy) && (R.rows() == dimy), "KalmanFull: R is not compatible");
+        bdm_assert (A.cols() == A.rows(), "KalmanFull: A is not square");
+        bdm_assert (B.rows() == A.rows(), "KalmanFull: B is not compatible");
+        bdm_assert (C.cols() == A.rows(), "KalmanFull: C is not compatible");
+        bdm_assert ( (D.rows() == C.rows()) && (D.cols() == B.cols()), "KalmanFull: D is not compatible");
+        bdm_assert ( (Q.cols() == A.rows()) && (Q.rows() == A.rows()), "KalmanFull: Q is not compatible");
+        bdm_assert ( (R.cols() == C.rows()) && (R.rows() == C.rows()), "KalmanFull: R is not compatible");
 }
 

library/bdm/estim/mixtures.cpp

@@ -30 +30 @@
                 //pick one datum
                 int ind = floor ( ndat * UniRNG.sample() );
-                Coms ( i )->bayes ( Data.get_col ( ind ), 1.0 );
+                Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
                 //flatten back to oringinal
                 Coms ( i )->flatten ( Com0 );
@@ -41 +41 @@
 }
 
-void MixEF::bayesB ( const mat &data , const vec &wData ) {
+void MixEF::bayesB ( const mat &data , const mat &cond, const vec &wData ) {
         int ndat = data.cols();
         int t, i, niter;
@@ -106 +106 @@
                         //#pragma omp parallel for
                         for ( i = 0; i < n; i++ ) {
-                                Coms ( i )-> bayes ( data.get_col ( t ), W ( i, t ) * wData ( t ) );
+                                Coms ( i )-> bayes_weighted ( data.get_col ( t ), empty_vec, W ( i, t ) * wData ( t ) );
                         }
                         weights.bayes ( W.get_col ( t ) * wData ( t ) );
@@ -122 +122 @@
 }
 
-void MixEF::bayes ( const vec &data ) {
+void MixEF::bayes ( const vec &data, const vec &cond=empty_vec ) {
 
 };
 
-void MixEF::bayes ( const mat &data ) {
-        this->bayesB ( data, ones ( data.cols() ) );
+void MixEF::bayes ( const mat &data, const vec &cond=empty_vec ) {
+        this->bayesB ( data, cond, ones ( data.cols() ) );
 };
 
 
-double MixEF::logpred ( const vec &dt ) const {
+double MixEF::logpred ( const vec &yt ) const {
 
         vec w = weights.posterior().mean();
         double exLL = 0.0;
         for ( int i = 0; i < n; i++ ) {
-                exLL += w ( i ) * exp ( Coms ( i )->logpred ( dt ) );
+                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt ) );
         }
         return log ( exLL );

library/bdm/estim/mixtures.h

@@ -108 +108 @@
         }
         //! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006
-        void bayes ( const vec &dt );
+        void bayes ( const vec &yt, const vec &cond );
         //! EM algorithm
-        void bayes ( const mat &dt );
+        void bayes ( const mat &yt, const vec &cond );
         //! batch weighted Bayes rule 
-        void bayesB ( const mat &dt, const vec &wData );
-        double logpred ( const vec &dt ) const;
+        void bayesB ( const mat &yt, const mat &cond, const vec &wData );
+        double logpred ( const vec &yt ) const;
         //! return correctly typed posterior (covariant return)
         const eprod& posterior() const {

library/bdm/estim/particles.cpp

@@ -46 +46 @@
 }
 
-void PF::bayes ( const vec &dt ) {
-        vec yt=dt.left(obs->dimension());
-        vec ut=dt.get(obs->dimension(),dt.length()-1);
+void PF::bayes( const vec &yt, const vec &cond ) {
+        const vec &ut=cond; //todo
         
         int i;
@@ -74 +73 @@
 // }
 
-void MPF::bayes ( const vec &dt ) {     
+void MPF::bayes ( const vec &yt, const vec &cond ) {    
         // follows PF::bayes in most places!!!  
         int i;
@@ -85 +84 @@
         #pragma parallel for
         for ( i = 0; i < n; i++ ) {
-                BMs(i) -> condition(pf->posterior()._sample(i));
-                BMs(i) -> bayes(dt);
+                BMs(i) -> bayes(this2bm.pushdown(yt), this2bm.get_cond(yt,cond));
                 lls ( i ) += BMs(i)->_ll();
         }

library/bdm/estim/particles.h

@@ -93 +93 @@
                 return eff < ( res_threshold*n );
         }
-        void bayes ( const vec &dt );
+        void bayes ( const vec &yt, const vec &cond );
         //!access function
         vec& __w() { return _w; }
@@ -131 +131 @@
                 u.add(obs_u); // join both u, and check if they do not overlap
 
-                set_drv(concat(obs->_rv(),u) );
+                set_yrv(concat(obs->_rv(),u) );
         }
         //! auxiliary function reading parameter 'resmethod' from configuration file
@@ -296 +296 @@
         mpfepdf jest;
 
-        //! Log means of BMs
-        bool opt_L_mea;
-
+        //! datalink from global yt and cond (Up) to BMs yt and cond (Down)
+        datalink_m2m this2bm;
+        //! datalink from global yt and cond (Up) to PFs yt and cond (Down)
+        datalink_m2m this2pf;
+        
 public:
         //! Default constructor.
@@ -320 +322 @@
         };
 
-        void bayes ( const vec &dt );
+        void bayes ( const vec &yt, const vec &cond );
         const epdf& posterior() const {
                 return jest;
@@ -328 +330 @@
         void set_options ( const string &opt ) {
                 BM::set_options(opt);
-                opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
         }
 
@@ -352 +353 @@
 
                 pf = new PF;
-                // rpior must be set before BM
+                // prior must be set before BM
                 pf->prior_from_set(set);
                 pf->resmethod_from_set(set);
                 pf->set_model(par,obs);
-                
+                                
                 shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory);
                 set_BM(*BM0);
@@ -367 +368 @@
                 //find potential input - what remains in rvc when we subtract rv
                 RV u = par->_rvc().remove_time().subt( par->_rv() );            
-                set_drv(concat(BM0->_drv(),u) );
+                set_yrv(concat(BM0->_yrv(),u) );
                 validate();
         }
@@ -377 +378 @@
                 }
                 jest.read_parameters();
-                for ( int i = 0; i < pf->__w().length(); i++ ) {
-                        BMs ( i )->condition ( pf->posterior()._sample ( i ) );
-                }
+                this2bm.set_connection(BMs(0)->_yrv(), BMs(0)->_rvc(), yrv, rvc);
+                this2bm.set_connection(pf->_yrv(), pf->_rvc(), yrv, rvc);
         }
 

library/bdm/itpp_ext.cpp

@@ -26 +26 @@
 namespace itpp
 {
+        vec empty_vec = vec(0);
+        
 Array<int> to_Arr (const ivec &indices)
 {

library/bdm/itpp_ext.h

@@ -19 +19 @@
 
 namespace itpp {
+        extern vec empty_vec;
+        
 Array<int> to_Arr ( const ivec &indices );
 ivec linspace ( int from, int to );

library/bdm/mex/mex_BM.h

@@ -75 +75 @@
                         est.from_setting(posterior);
                 }
-                void bayes(const vec &dt)  {
+                void bayes(const vec &yt, const vec &cond)  {
                 //void bayes()  {
                         mxArray *tmp, *old;

library/bdm/stat/exp_family.cpp

@@ -14 +14 @@
         ///////////
 
-void BMEF::bayes ( const vec &dt ) {
-        this->bayes ( dt, 1.0 );
+void BMEF::bayes( const vec &yt, const vec &cond ) {
+        this->bayes_weighted ( yt, cond, 1.0 );
 };
 

library/bdm/stat/exp_family.h

@@ -96 +96 @@
 
                 //! Weighted update of sufficient statistics (Bayes rule)
-                virtual void bayes (const vec &data, const double w) {};
+                virtual void bayes_weighted (const vec &data, const vec &cond=empty_vec, const double w=1.0) {};
                 //original Bayes
-                void bayes (const vec &dt);
+                void bayes (const vec &yt, const vec &cond=empty_vec);
 
                 //!Flatten the posterior according to the given BMEF (of the same type!)
@@ -436 +436 @@
                 //! Sets sufficient statistics to match that of givefrom mB0
                 void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
-                void bayes (const vec &dt) {
+                void bayes (const vec &yt, const vec &cond=empty_vec) {
                         if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
-                        beta += dt;
+                        beta += yt;
                         if (evalll) {ll = est.lognc() - last_lognc;}
                 }
-                double logpred (const vec &dt) const {
+                double logpred (const vec &yt) const {
                         eDirich pred (est);
                         vec &beta = pred._beta();
@@ -452 +452 @@
                                 else{lll = pred.lognc();}
 
-                        beta += dt;
+                        beta += yt;
                         return pred.lognc() - lll;
                 }
@@ -832 +832 @@
                 }
                 void condition (const vec &cond) {
-                        iepdf._mu() = A * cond + mu_const;
+                        if (cond.length()>0) {
+                                iepdf._mu() = A * cond + mu_const;
+                        } else {
+                                iepdf._mu() =  mu_const;
+                        }
                         double zeta;
                         //ugly hack!

library/bdm/stat/merger.cpp

@@ -103 +103 @@
                 //Re-Initialize Mixture model
                 Mix.flatten ( &Mix_init );
-                Mix.bayesB ( Smp_ex, w*Npoints );
+                Mix.bayesB ( Smp_ex,empty_vec, w*Npoints );
                 delete Mpred;
                 Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!

library/bdm/stat/merger.h

@@ -342 +342 @@
         }
         //! loglikelihood computed on mixture models
-        double evallog ( const vec &dt ) const {
-                vec dtf = ones ( dt.length() + 1 );
-                dtf.set_subvector ( 0, dt );
+        double evallog ( const vec &yt ) const {
+                vec dtf = ones ( yt.length() + 1 );
+                dtf.set_subvector ( 0, yt );
                 return Mix.logpred ( dtf );
         }

library/tests/arx_elem_test.cpp

@@ -9 +9 @@
         ARX Ar;
         Ar.set_statistics ( 1, V0, -1.0 );
+        Ar.set_constant(true);
+        Ar.validate();
 
         mat mu ( 1, 1 );
@@ -36 +38 @@
 
         mlstudent* Ap = Ar.predictor_student();
-        vec Ap_x = Ap->evallogcond_m ( X, vec_1 ( 1.0 ) );
+        vec Ap_x = Ap->evallogcond_m ( X, empty_vec );
         vec ll_x = Ar.logpred_m ( X2 );
 

library/tests/arx_test.cpp

@@ -27 +27 @@
         ARX Ar;
         Ar.set_statistics ( 1, V0, nu0 );               // Estimator
+        Ar.set_constant(false);
+        Ar.validate();
         const epdf& f_thr = Ar.posterior();          // refrence to posterior of the estimator
 
@@ -34 +36 @@
         Yt.set_subvector ( 0, randn ( ord ) ); //initial values
         vec rgr ( ord );                // regressor
-        vec Psi ( ord + 1 );            // extended regressor
 
         //print moments of the prior distribution
@@ -49 +50 @@
                 Yt ( t ) = th * rgr + sqr * NorRNG();
 
-                Psi = concat ( Yt ( t ), rgr ); // Inefficient! Used for compatibility with Matlab!
-                Ar.bayes ( Psi );             // Bayes rule
+                Ar.bayes ( vec_1(Yt(t)), rgr );             // Bayes rule
 
                 // Build predictor

library/tests/dirfile-format.matrix

@@ -1 @@
-r_0 RAW d 1
-r_1 RAW d 1
-th_alog RAW d 1
-th_blog RAW d 1
+r.0 RAW d 1
+r.1 RAW d 1
+th.alog RAW d 1
+th.blog RAW d 1

library/tests/mixtures_test.cpp

@@ -99 +99 @@
         // Add ones for constant coefficients
         mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) );
-        Post.bayes ( Data );
+        Post.bayes ( Data , empty_vec);
 
         cout << "Posterior mixture:" << endl;
@@ -115 +115 @@
         mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );;
 
-        RND.bayes ( Data );
+        RND.bayes ( Data, empty_vec );
         cout << endl << "== BAYES ==" << endl;
         disp_mix2D ( RND );

library/tests/tutorial/arx_simple.cpp

@@ -9 +9 @@
         ARX Ar;
         Ar.set_statistics ( 1, V0 ); //nu is default (set to have finite moments)
+        Ar.set_constant(true);
+        Ar.validate();
         // forgetting is default: 1.0
         mat Data = concat_vertical ( randn ( 1, 100 ), ones ( 1, 100 ) );

library/tests/tutorial/kalman_simple.cpp

@@ -20 +20 @@
         KF.set_parameters ( A, B, C, D,/*covariances*/ Q, R );
         KF.set_statistics ( mu0, P0 );
+        KF.validate();
         // Estimation loop
         for ( int i = 0; i < 100; i++ ) {
-                KF.bayes ( randn ( dx + du ) );
+                KF.bayes ( randn ( dy), randn( du ) );
         }
         //print results