Changeset 679

Show
Ignore:
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!!!

Files:
24 modified

Legend:

Unmodified
Added
Removed

  • applications/pmsm/pmsmDS.h

    r676 r679  
    285285        //connect 
    286286        shared_ptr<RV> drv = UI::build<RV>(root, "drv"); 
    287         set_drv(*drv); 
     287        set_yrv(*drv); 
    288288        shared_ptr<RV> rv = UI::build<RV>(root, "rv"); 
    289289        set_rv(*rv); 
    … …  
    363363        //connect 
    364364                shared_ptr<RV> drv = UI::build<RV>(root, "drv"); 
    365         set_drv(*drv); 
     365        set_yrv(*drv); 
    366366                shared_ptr<RV> rv = UI::build<RV>(root, "rv"); 
    367367        set_rv(*rv); 

  • library/bdm/base/bdmbase.h

    r678 r679  
    406406 
    407407        mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), dim( m.dim), rv( m.rv ) { } 
     408         
     409        //! copy of the current object - make sure to implement 
     410        virtual mpdf* _copy_() const {return new mpdf(*this);} 
    408411        //!@} 
    409412 
    … …  
    421424 
    422425        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently. 
    423         virtual double evallogcond ( const vec &dt, const vec &cond ) { 
     426        virtual double evallogcond ( const vec &yt, const vec &cond ) { 
    424427                bdm_error ( "Not implemented" ); 
    425428                return 0.0; 
    … …  
    427430 
    428431        //! Matrix version of evallogcond 
    429         virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) { 
    430                 vec v ( Dt.cols() ); 
    431                 for ( int i = 0; i < Dt.cols(); i++ ) { 
    432                         v ( i ) = evallogcond ( Dt.get_col ( i ), cond ); 
     432        virtual vec evallogcond_m ( const mat &Yt, const vec &cond ) { 
     433                vec v ( Yt.cols() ); 
     434                for ( int i = 0; i < Yt.cols(); i++ ) { 
     435                        v ( i ) = evallogcond ( Yt.get_col ( i ), cond ); 
    433436                } 
    434437                return v; 
    … …  
    436439 
    437440        //! Array<vec> version of evallogcond 
    438         virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) { 
     441        virtual vec evallogcond_m ( const Array<vec> &Yt, const vec &cond ) { 
    439442                bdm_error ( "Not implemented" ); 
    440443                return vec(); 
    … …  
    458461        } 
    459462 
     463        //! access function 
     464        void set_dim(int d) {dim=d;} 
     465        //! access function 
     466        void set_dimc(int d) {dimc=d;} 
    460467        //! Load from structure with elements: 
    461468        //!  \code 
    … …  
    506513                dim = dim0; 
    507514        } 
     515        epdf* _copy_() const {return new epdf(*this);} 
    508516        //!@} 
    509517 
    … …  
    10201028This object represents exact or approximate evaluation of the Bayes rule: 
    10211029\f[ 
    1022 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})} 
     1030f(\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})} 
    10231031\f] 
    1024  
     1032where: 
     1033 * \f$ y_t \f$ is the variable  
    10251034Access to the resulting posterior density is via function \c posterior(). 
    10261035 
    … …  
    10281037It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor(). 
    10291038 
    1030 Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$: 
     1039Alternatively, it can evaluate posterior density with rvc replaced by the given values, \f$ c_t \f$: 
    10311040\f[ 
    10321041f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1}) 
    10331042\f] 
    10341043 
    1035 The value of \f$ c_t \f$ is set by function condition(). 
    10361044 
    10371045*/ 
    … …  
    10401048protected: 
    10411049        //! Random variable of the data (optional) 
    1042         RV drv; 
     1050        RV yrv; 
     1051        //! size of the data record 
     1052        int dimy; 
     1053        //! Name of extension variable 
     1054        RV rvc; 
     1055        //! size of the conditioning vector 
     1056        int dimc; 
     1057         
    10431058        //!Logarithm of marginalized data likelihood. 
    10441059        double ll; 
    10451060        //!  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. 
    10461061        bool evalll; 
    1047          
     1062 
    10481063public: 
    10491064        //! \name Constructors 
    10501065        //! @{ 
    10511066 
    1052         BM() : ll ( 0 ), evalll ( true ) { }; 
    1053         BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {} 
     1067        BM() : yrv(),dimy(0),rvc(),dimc(0), ll ( 0 ), evalll ( true ) { }; 
     1068        BM ( const BM &B ) :  yrv ( B.yrv ), dimy(B.dimy), rvc ( B.rvc ), ll ( B.ll ), evalll ( B.evalll ) {} 
    10541069        //! \brief Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! 
    10551070        //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode 
    1056         virtual BM* _copy_() const { 
    1057                 return NULL; 
    1058         }; 
     1071        virtual BM* _copy_() const { return NULL; }; 
    10591072        //!@} 
    10601073 
    … …  
    10651078        @param dt vector of input data 
    10661079        */ 
    1067         virtual void bayes ( const vec &dt ) = 0; 
     1080        virtual void bayes ( const vec &yt, const vec &cond=empty_vec ) = 0; 
    10681081        //! Batch Bayes rule (columns of Dt are observations) 
    10691082        virtual void bayesB ( const mat &Dt ); 
    10701083        //! Evaluates predictive log-likelihood of the given data record 
    1071         //! I.e. marginal likelihood of the data with the posterior integrated out. 
    1072         virtual double logpred ( const vec &dt ) const { 
     1084        //! I.e. marginal likelihood of the data with the posterior integrated out.  
     1085        //! This function evaluates only \f$ y_t \f$, condition is assumed to be the last used in bayes().  
     1086        //! See bdm::BM::predictor for conditional version. 
     1087        virtual double logpred ( const vec &yt ) const { 
    10731088                bdm_error ( "Not implemented" ); 
    10741089                return 0.0; 
    … …  
    10761091 
    10771092        //! Matrix version of logpred 
    1078         vec logpred_m ( const mat &dt ) const { 
    1079                 vec tmp ( dt.cols() ); 
    1080                 for ( int i = 0; i < dt.cols(); i++ ) { 
    1081                         tmp ( i ) = logpred ( dt.get_col ( i ) ); 
     1093        vec logpred_m ( const mat &Yt ) const { 
     1094                vec tmp ( Yt.cols() ); 
     1095                for ( int i = 0; i < Yt.cols(); i++ ) { 
     1096                        tmp ( i ) = logpred ( Yt.get_col ( i ) ); 
    10821097                } 
    10831098                return tmp; 
    … …  
    10961111        //!@} 
    10971112 
    1098         //! \name Extension to conditional BM 
    1099         //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF). 
    1100         //! Alternatively, it can be used for automated connection to DS when the condition is observed 
    1101         //!@{ 
    1102  
    1103         //! Name of extension variable 
    1104         RV rvc; 
    1105         //! access function 
    1106         const RV& _rvc() const { 
    1107                 return rvc; 
    1108         } 
    1109  
    1110         //! Substitute \c val for \c rvc. 
    1111         virtual void condition ( const vec &val ) { 
    1112                 bdm_error ( "Not implemented!" ); 
    1113         } 
    1114  
    1115         //!@} 
    1116  
    11171113 
    11181114        //! \name Access to attributes 
    11191115        //!@{ 
     1116                //! access function 
     1117                const RV& _rvc() const { 
     1118                        return rvc; 
     1119                } 
     1120                //! access function 
     1121                int dimensionc() const { 
     1122                        return dimc; 
     1123                } 
     1124                //! access function 
     1125                int dimensiony() const { 
     1126                        return dimy; 
     1127                } 
     1128                //! access function 
     1129                int dimension() const { 
     1130                        return posterior().dimension(); 
     1131                } 
     1132                //! access function 
     1133        const RV& _yrv() const { 
     1134                return yrv; 
     1135        } 
    11201136        //! access function 
    1121         const RV& _drv() const { 
    1122                 return drv; 
    1123         } 
    1124         //! access function 
    1125         void set_drv ( const RV &rv ) { 
    1126                 drv = rv; 
     1137        void set_yrv ( const RV &rv ) { 
     1138                yrv = rv; 
    11271139        } 
    11281140        //! access to rv of the posterior 
    11291141        void set_rv ( const RV &rv ) { 
    1130                 const_cast<epdf&> ( posterior() ).set_rv ( rv ); 
     1142                const_cast<epdf&>(posterior()).set_rv ( rv ); 
     1143        } 
     1144        //! access function 
     1145        void set_dim ( int dim ) { 
     1146                const_cast<epdf&>(posterior()).set_dim ( dim ); 
    11311147        } 
    11321148        //! return internal log-likelihood of the last data vector 
    … …  
    11801196                shared_ptr<RV> r = UI::build<RV> ( set, "drv", UI::optional ); 
    11811197                if ( r ) { 
    1182                         set_drv ( *r ); 
     1198                        set_yrv ( *r ); 
    11831199                } 
    11841200                string opt; 
    … …  
    12151231 
    12161232template<class EPDF> 
    1217 double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) { 
     1233double mpdf_internal<EPDF>::evallogcond ( const vec &yt, const vec &cond ) { 
    12181234        double tmp; 
    12191235        condition ( cond ); 
    1220         tmp = iepdf.evallog ( dt ); 
     1236        tmp = iepdf.evallog ( yt ); 
    12211237        return tmp; 
    12221238} 
    12231239 
    12241240template<class EPDF> 
    1225 vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) { 
     1241vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Yt, const vec &cond ) { 
    12261242        condition ( cond ); 
    1227         return iepdf.evallog_m ( Dt ); 
     1243        return iepdf.evallog_m ( Yt ); 
    12281244} 
    12291245 
    12301246template<class EPDF> 
    1231 vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) { 
     1247vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Yt, const vec &cond ) { 
    12321248        condition ( cond ); 
    1233         return iepdf.evallog_m ( Dt ); 
     1249        return iepdf.evallog_m ( Yt ); 
    12341250} 
    12351251 

  • library/bdm/design/ctrlbase.cpp

    r620 r679  
    55void LQG::set_system(shared_ptr<StateSpace<fsqmat> > S0){ 
    66        S = S0; 
    7         dimx= S->_dimx(); 
    8         dimy= S->_dimy(); 
    9         dimu= S->_dimu(); 
     7        dimx= S->_A().rows(); 
     8        dimy= S->_C().rows(); 
     9        dimu= S->_B().cols(); 
    1010        pr=zeros(dimx+dimu+dimy,  dimu+dimx+dimu+dimy); 
    1111        pr.set_submatrix(dimx, dimu+dimx, eye(dimu+dimy)); 

  • library/bdm/estim/arx.cpp

    r665 r679  
    22namespace bdm { 
    33 
    4 void ARX::bayes ( const vec &dt, const double w ) { 
     4        void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) { 
    55        double lnc; 
    6                  
     6        //cache 
     7        ldmat &V=est._V();  
     8        double &nu=est._nu(); 
     9 
     10        dyad.set_subvector(0,yt); 
     11        dyad.set_subvector(dimy,cond); 
     12        // possible "1" is there from the beginning 
     13         
    714        if ( frg < 1.0 ) { 
    815                est.pow ( frg ); // multiply V and nu 
     16 
    917                 
    1018                //stabilize 
    1119                ldmat V0=alter_est._V(); //$ copy 
    1220                double &nu0=alter_est._nu(); 
     21                 
    1322                V0*=(1-frg); 
    1423                V += V0; //stabilization 
    … …  
    2029                } 
    2130        } 
    22         if (have_constant) { 
    23                 _dt.set_subvector(0,dt); 
    24                 V.opupdt ( _dt, w ); 
    25         } else { 
    26                 V.opupdt ( dt, w ); 
    27         } 
     31        V.opupdt ( dyad, w ); 
    2832        nu += w; 
    2933 
    … …  
    3640} 
    3741 
    38 double ARX::logpred ( const vec &dt ) const { 
     42double ARX::logpred ( const vec &yt ) const { 
    3943        egiw pred ( est ); 
    4044        ldmat &V = pred._V(); 
    … …  
    4246 
    4347        double lll; 
     48        vec dyad_p = dyad; 
     49        dyad_p.set_subvector(0,yt); 
    4450 
    4551        if ( frg < 1.0 ) { 
    … …  
    5359                } 
    5460 
    55         V.opupdt ( dt, 1.0 ); 
     61        V.opupdt ( dyad_p, 1.0 ); 
    5662        nu += 1.0; 
    5763        // log(sqrt(2*pi)) = 0.91893853320467 
    … …  
    6773        const ARX* A0 = dynamic_cast<const ARX*> ( B0 ); 
    6874 
    69         bdm_assert_debug ( V.rows() == A0->V.rows(), "ARX::set_statistics Statistics  differ" ); 
    70         set_statistics ( A0->dimx, A0->V, A0->nu ); 
     75        bdm_assert_debug ( dimension() == A0->dimension(), "Statistics of different dimensions" ); 
     76        set_statistics ( A0->dimensiony(), A0->posterior()._V(), A0->posterior()._nu() ); 
    7177} 
    7278 
    7379enorm<ldmat>* ARX::epredictor ( const vec &rgr ) const { 
    74         int dim = dimx;//est.dimension(); 
    75         mat mu ( dim, V.rows() - dim ); 
    76         mat R ( dim, dim ); 
     80        mat mu ( dimy, posterior()._V().rows() - dimy ); 
     81        mat R ( dimy, dimy ); 
    7782 
    7883        enorm<ldmat>* tmp; 
    7984        tmp = new enorm<ldmat> ( ); 
    8085        //TODO: too hackish 
    81         if ( drv._dsize() > 0 ) { 
     86        if ( yrv._dsize() > 0 ) { 
    8287        } 
    8388 
    … …  
    9196 
    9297mlnorm<ldmat>* ARX::predictor ( ) const { 
    93         int dim = est.dimension(); 
    94          
    95         mat mu ( dim, V.rows() - dim ); 
    96         mat R ( dim, dim ); 
     98        mat mu ( dimy, posterior()._V().rows() - dimy ); 
     99        mat R ( dimy, dimy ); 
    97100        mlnorm<ldmat>* tmp; 
    98101        tmp = new mlnorm<ldmat> ( ); 
    … …  
    106109                tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) ); 
    107110        } else { 
    108                 tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) ); 
     111                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) ); 
    109112        } 
    110113        return tmp; 
    … …  
    112115 
    113116mlstudent* ARX::predictor_student ( ) const { 
    114         int dim = est.dimension(); 
    115  
    116         mat mu ( dim, V.rows() - dim ); 
    117         mat R ( dim, dim ); 
     117        const ldmat &V = posterior()._V(); 
     118         
     119        mat mu ( dimy, V.rows() - dimy ); 
     120        mat R ( dimy, dimy ); 
    118121        mlstudent* tmp; 
    119122        tmp = new mlstudent ( ); 
    … …  
    122125        mu = mu.T(); 
    123126 
    124         int xdim = dimx; 
    125127        int end = V._L().rows() - 1; 
    126         ldmat Lam ( V._L() ( xdim, end, xdim, end ), V._D() ( xdim, end ) );  //exp val of R 
     128        ldmat Lam ( V._L() ( dimy, end, dimy, end ), V._D() ( dimy, end ) );  //exp val of R 
    127129 
    128130 
    … …  
    132134                        tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam ); 
    133135                } else { 
    134                         tmp->set_parameters ( mat ( dim, 0 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam ); 
     136                        tmp->set_parameters ( mat ( dimy, dimc ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam ); 
    135137                } 
    136138        } else { 
    137139                // no constant term 
    138                 tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam ); 
     140                tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ), Lam ); 
    139141        } 
    140142        return tmp; 
    … …  
    214216 
    215217void ARX::from_setting ( const Setting &set ) { 
    216         shared_ptr<RV> yrv = UI::build<RV> ( set, "rv", UI::compulsory ); 
     218        shared_ptr<RV> yrv_ = UI::build<RV> ( set, "rv", UI::compulsory ); 
    217219        shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory ); 
    218         int ylen = yrv->_dsize(); 
     220        dimy = yrv_->_dsize(); 
    219221        // rgrlen - including constant!!! 
    220         int rgrlen = rrv->_dsize(); 
    221          
    222         dimx = ylen; 
    223         set_rv ( *yrv, *rrv ); 
     222        dimc = rrv->_dsize(); 
     223         
     224        yrv = *yrv_; 
     225        rvc = *rrv; 
    224226         
    225227        string opt; 
    … …  
    233235                have_constant=constant>0; 
    234236        } 
    235         if (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);} 
     237        int rgrlen = dimc+int(have_constant==true); 
    236238 
    237239        //init 
    238240        shared_ptr<egiw> pri=UI::build<egiw>(set, "prior", UI::optional); 
    239241        if (pri) { 
    240                 bdm_assert(pri->_dimx()==ylen,"prior is not compatible"); 
    241                 bdm_assert(pri->_V().rows()==ylen+rgrlen,"prior is not compatible"); 
     242                bdm_assert(pri->_dimx()==dimy,"prior is not compatible"); 
     243                bdm_assert(pri->_V().rows()==dimy+rgrlen,"prior is not compatible"); 
    242244                est.set_parameters( pri->_dimx(),pri->_V(), pri->_nu()); 
    243245        }else{ 
    244                 est.set_parameters( ylen, zeros(ylen+rgrlen)); 
     246                est.set_parameters( dimy, zeros(dimy+rgrlen)); 
    245247                set_prior_default(est); 
    246248        } 
    … …  
    248250        shared_ptr<egiw> alt=UI::build<egiw>(set, "alternative", UI::optional); 
    249251        if (alt) { 
    250                 bdm_assert(alt->_dimx()==ylen,"alternative is not compatible"); 
    251                 bdm_assert(alt->_V().rows()==ylen+rgrlen,"alternative is not compatible"); 
     252                bdm_assert(alt->_dimx()==dimy,"alternative is not compatible"); 
     253                bdm_assert(alt->_V().rows()==dimy+rgrlen,"alternative is not compatible"); 
    252254                alter_est.set_parameters( alt->_dimx(),alt->_V(), alt->_nu()); 
    253255        } else { 
    … …  
    264266        shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional ); 
    265267        if (!rv_par){ 
    266                 est.set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) ); 
     268                est.set_rv ( RV ( "{theta r }", vec_2 ( dimy*rgrlen, dimy*dimy ) ) ); 
    267269        } else { 
    268270                est.set_rv ( *rv_par ); 
    … …  
    270272        validate(); 
    271273} 
    272  
    273 } 
     274} 
     275 

  • library/bdm/estim/arx.h

    r665 r679  
    4343class ARX: public BMEF { 
    4444protected: 
    45         //!size of output variable (needed in regressors) 
    46         int dimx; 
    47         //!description of modelled data \f$ y_t \f$ in the likelihood function 
    48         //! Do NOT access directly, only via \c get_yrv(). 
    49         RV _yrv; 
    50         //! rv of regressor 
    51         RV rgrrv; 
    52         //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density 
    53         egiw est; 
    54         //! cached value of est.V 
    55         ldmat &V; 
    56         //! cached value of est.nu 
    57         double &nu; 
    5845        //! switch if constant is modelled or not 
    5946        bool have_constant; 
    60         //! cached value of data vector for have_constant =true 
    61         vec _dt; 
     47        //! vector of dyadic update 
     48        vec dyad; 
     49        //! posterior density 
     50        egiw est; 
    6251        //! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy] 
    6352        egiw alter_est; 
    … …  
    6554        //! \name Constructors 
    6655        //!@{ 
    67         ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), est (), V ( est._V() ), nu ( est._nu() ), have_constant(true), _dt() {}; 
    68         ARX ( const ARX &A0 ) : BMEF (A0.frg), est (A0.est), V ( est._V() ), nu ( est._nu() ), have_constant(A0.have_constant), _dt(A0._dt) { 
    69                 dimx = A0.dimx; 
    70                 _yrv = A0._yrv; 
    71                 rgrrv = A0.rgrrv; 
    72                 set_drv(A0._drv()); 
    73         }; 
     56        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant(true), dyad(), est() {}; 
     57        ARX ( const ARX &A0 ) : BMEF (A0.frg),  have_constant(A0.have_constant), dyad(A0.dyad),est(est) { }; 
    7458        ARX* _copy_() const; 
    7559        void set_parameters ( double frg0 ) { 
    … …  
    7963                have_constant=const0; 
    8064        } 
    81         void set_statistics ( int dimx0, const ldmat V0, double nu0 = -1.0 ) { 
    82                 est.set_parameters ( dimx0, V0, nu0 ); 
     65        void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) { 
     66                est.set_parameters ( dimy0, V0, nu0 ); 
    8367                last_lognc = est.lognc(); 
    84                 dimx = dimx0; 
     68                dimy = dimy0; 
    8569        } 
    8670        //!@} 
    … …  
    9377 
    9478        //! Weighted Bayes \f$ dt = [y_t psi_t] \f$. 
    95         void bayes ( const vec &dt, const double w ); 
    96         void bayes ( const vec &dt ) { 
    97                 bayes ( dt, 1.0 ); 
     79        void bayes_weighted ( const vec &yt, const vec &cond=empty_vec, const double w=1.0 ); 
     80        void bayes( const vec &yt, const vec &cond=empty_vec ) { 
     81                bayes_weighted ( yt,cond, 1.0 ); 
    9882        }; 
    99         double logpred ( const vec &dt ) const; 
     83        double logpred ( const vec &yt ) const; 
    10084        void flatten ( const BMEF* B ) { 
    10185                const ARX* A = dynamic_cast<const ARX*> ( B ); 
    10286                // nu should be equal to B.nu 
    103                 est.pow ( A->nu / nu ); 
     87                est.pow ( A->posterior()._nu() / posterior()._nu() ); 
    10488                if ( evalll ) { 
    10589                        last_lognc = est.lognc(); 
    … …  
    11094        //! Predictor for empty regressor 
    11195        enorm<ldmat>* epredictor() const { 
    112                 bdm_assert_debug ( dimx == V.rows() - 1, "Regressor is not only 1" ); 
     96                bdm_assert_debug ( dimy == posterior()._V().rows() - 1, "Regressor is not only 1" ); 
    11397                return epredictor ( vec_1 ( 1.0 ) ); 
    11498        } 
    … …  
    127111                const egiw& posterior() const { 
    128112                return est; 
    129         } 
    130         //!@} 
    131  
    132         //!\name Connection 
    133         //!@{ 
    134         void set_rv ( const RV &yrv0 , const RV &rgrrv0 ) { 
    135                 _yrv = yrv0; 
    136                 rgrrv=rgrrv0; 
    137                 set_drv(concat(yrv0, rgrrv)); 
    138         } 
    139  
    140         RV& get_yrv() { 
    141                 //if yrv is not ready create it 
    142                 if ( _yrv._dsize() != dimx ) { 
    143                         int i = 0; 
    144                         while ( _yrv._dsize() < dimx ) { 
    145                                 _yrv.add ( drv ( vec_1 ( i ) ) ); 
    146                                 i++; 
    147                         } 
    148                 } 
    149                 //yrv should be ready by now 
    150                 bdm_assert_debug ( _yrv._dsize() == dimx, "incompatible drv" ); 
    151                 return _yrv; 
    152113        } 
    153114        //!@} 
    … …  
    174135 
    175136        void validate() { 
    176                 bdm_assert(dimx == _yrv._dsize(), "RVs of parameters and regressor do not match"); 
    177                  
     137                //if dimc not set set it from V 
     138                if (dimc==0){ 
     139                        dimc = posterior()._V().rows()-dimy-int(have_constant==true); 
     140                } 
     141                         
     142                if (have_constant) { 
     143                        dyad = ones(dimy+dimc+1); 
     144                } else {  
     145                        dyad = zeros(dimy+dimc); 
     146                } 
     147                         
    178148        } 
    179149        //! function sets prior and alternative density  

  • library/bdm/estim/kalman.cpp

    r653 r679  
    88 
    99 
    10 void KalmanFull::bayes ( const vec &dt ) { 
    11         bdm_assert_debug ( dt.length() == ( dimy + dimu ), "KalmanFull::bayes wrong size of dt" ); 
    12  
    13         vec u = dt.get ( dimy, dimy + dimu - 1 ); 
    14         vec y = dt.get ( 0, dimy - 1 ); 
    15         vec& mu = est->_mu(); 
    16         mat &P = est->_R(); 
     10void KalmanFull::bayes ( const vec &yt, const vec &cond ) { 
     11        bdm_assert_debug ( yt.length() == ( dimy ), "KalmanFull::bayes wrong size of dt" ); 
     12 
     13        const vec &u = cond; // in this case cond=ut 
     14        const vec &y = yt; 
     15         
     16        vec& mu = est._mu(); 
     17        mat &P = est._R(); 
    1718        mat& _Ry = fy._R(); 
    1819        vec& yp = fy._mu(); 
    … …  
    4243        phxu = phxu0; 
    4344 
    44         dimx = pfxu0->_dimx(); 
     45        set_dim( pfxu0->_dimx()); 
    4546        dimy = phxu0->dimension(); 
    46         dimu = pfxu0->_dimu(); 
    47         est->set_parameters( zeros(dimx), eye(dimx) ); 
    48  
    49         A.set_size ( dimx, dimx ); 
    50         C.set_size ( dimy, dimx ); 
     47        dimc = pfxu0->_dimu(); 
     48        est.set_parameters( zeros(dimension()), eye(dimension()) ); 
     49 
     50        A.set_size ( dimension(), dimension() ); 
     51        C.set_size ( dimy, dimension() ); 
    5152        //initialize matrices A C, later, these will be only updated! 
    52         pfxu->dfdx_cond ( est->_mu(), zeros ( dimu ), A, true ); 
     53        pfxu->dfdx_cond ( est._mu(), zeros ( dimc ), A, true ); 
    5354        B.clear(); 
    54         phxu->dfdx_cond ( est->_mu(), zeros ( dimu ), C, true ); 
     55        phxu->dfdx_cond ( est._mu(), zeros ( dimc ), C, true ); 
    5556        D.clear(); 
    5657 
    … …  
    5960} 
    6061 
    61 void EKFfull::bayes ( const vec &dt ) { 
    62         bdm_assert_debug ( dt.length() == ( dimy + dimu ), "EKFull::bayes wrong size of dt" ); 
    63  
    64         vec u = dt.get ( dimy, dimy + dimu - 1 ); 
    65         vec y = dt.get ( 0, dimy - 1 ); 
    66         vec &mu = est->_mu(); 
    67         mat &P = est->_R(); 
     62void EKFfull::bayes ( const vec &yt, const vec &cond) { 
     63        bdm_assert_debug ( yt.length() == ( dimy ), "EKFull::bayes wrong size of dt" ); 
     64        bdm_assert_debug ( cond.length() == ( dimc ), "EKFull::bayes wrong size of dt" ); 
     65         
     66        const vec &u = cond; 
     67        const vec &y = yt; //lazy to change it 
     68        vec &mu = est._mu(); 
     69        mat &P = est._R(); 
    6870        mat& _Ry = fy._R(); 
    6971        vec& yp = fy._mu(); 
    7072         
    71         pfxu->dfdx_cond ( mu, zeros ( dimu ), A, true ); 
    72         phxu->dfdx_cond ( mu, zeros ( dimu ), C, true ); 
     73        pfxu->dfdx_cond ( mu, zeros ( dimc ), A, true ); 
     74        phxu->dfdx_cond ( mu, zeros ( dimc ), C, true ); 
    7375 
    7476        //Time update 
    … …  
    9496        ( ( StateSpace<chmat>* ) this )->set_parameters ( A0, B0, C0, D0, Q0, R0 ); 
    9597         
    96         _K=zeros(dimx,dimy); 
    97         // Cholesky special! 
    98         initialize(); 
     98        _K=zeros(dimension(),dimy); 
    9999} 
    100100 
    101101void KalmanCh::initialize(){ 
    102         preA = zeros ( dimy + dimx + dimx, dimy + dimx ); 
     102        preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() ); 
    103103//      preA.clear(); 
    104104        preA.set_submatrix ( 0, 0, R._Ch() ); 
    105         preA.set_submatrix ( dimy + dimx, dimy, Q._Ch() ); 
    106 } 
    107  
    108 void KalmanCh::bayes ( const vec &dt ) { 
    109  
    110         vec u = dt.get ( dimy, dimy + dimu - 1 ); 
    111         vec y = dt.get ( 0, dimy - 1 ); 
     105        preA.set_submatrix ( dimy + dimension(), dimy, Q._Ch() ); 
     106} 
     107 
     108void KalmanCh::bayes ( const vec &yt, const vec &cond ) { 
     109        bdm_assert_debug(yt.length()==dimy, "yt mismatch"); 
     110        bdm_assert_debug(cond.length()==dimc, "yt mismatch"); 
     111         
     112        const vec &u = cond; 
     113        const vec &y = yt; 
    112114        vec pom ( dimy ); 
    113115 
    114         chmat &_P=est->_R(); 
    115         vec &_mu = est->_mu(); 
    116         mat _K(dimx,dimy); 
     116        chmat &_P=est._R(); 
     117        vec &_mu = est._mu(); 
     118        mat _K(dimension(),dimy); 
    117119        chmat &_Ry=fy._R(); 
    118120        vec &_yp = fy._mu(); 
    … …  
    130132 
    131133        ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 ); 
    132         _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimx - 1 ) ).T(); 
    133         ( _P._Ch() ) = postA ( dimy, dimy + dimx - 1, dimy, dimy + dimx - 1 ); 
     134        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimension() - 1 ) ).T(); 
     135        ( _P._Ch() ) = postA ( dimy, dimy + dimension() - 1, dimy, dimy + dimension() - 1 ); 
    134136 
    135137        _mu = A * ( _mu ) + B * u; 
    … …  
    154156        phxu = phxu0; 
    155157 
    156         dimx = pfxu0->dimension(); 
     158        set_dim( pfxu0->_dimx()); 
    157159        dimy = phxu0->dimension(); 
    158         dimu = pfxu0->_dimu(); 
    159          
    160         vec &_mu = est->_mu(); 
     160        dimc = pfxu0->_dimu(); 
     161         
     162        vec &_mu = est._mu(); 
    161163        // if mu is not set, set it to zeros, just for constant terms of A and C 
    162         if ( _mu.length() != dimx ) _mu = zeros ( dimx );  
    163         A = zeros ( dimx, dimx ); 
    164         C = zeros ( dimy, dimx ); 
    165         preA = zeros ( dimy + dimx + dimx, dimy + dimx ); 
     164        if ( _mu.length() != dimension() ) _mu = zeros ( dimension() );  
     165        A = zeros ( dimension(), dimension() ); 
     166        C = zeros ( dimy, dimension() ); 
     167        preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() ); 
    166168 
    167169        //initialize matrices A C, later, these will be only updated! 
    168         pfxu->dfdx_cond ( _mu, zeros ( dimu ), A, true ); 
     170        pfxu->dfdx_cond ( _mu, zeros ( dimc ), A, true ); 
    169171//      pfxu->dfdu_cond ( *_mu,zeros ( dimu ),B,true ); 
    170172        B.clear(); 
    171         phxu->dfdx_cond ( _mu, zeros ( dimu ), C, true ); 
     173        phxu->dfdx_cond ( _mu, zeros ( dimc ), C, true ); 
    172174//      phxu->dfdu_cond ( *_mu,zeros ( dimu ),D,true ); 
    173175        D.clear(); 
    … …  
    179181        preA.clear(); 
    180182        preA.set_submatrix ( 0, 0, R._Ch() ); 
    181         preA.set_submatrix ( dimy + dimx, dimy, Q._Ch() ); 
    182 } 
    183  
    184  
    185 void EKFCh::bayes ( const vec &dt ) { 
     183        preA.set_submatrix ( dimy + dimension(), dimy, Q._Ch() ); 
     184} 
     185 
     186 
     187void EKFCh::bayes ( const vec &yt, const vec &cond ) { 
    186188 
    187189        vec pom ( dimy ); 
    188         vec u = dt.get ( dimy, dimy + dimu - 1 ); 
    189         vec y = dt.get ( 0, dimy - 1 ); 
    190         vec &_mu = est->_mu(); 
    191         chmat &_P = est->_R(); 
     190        const vec &u = cond; 
     191        const vec &y = yt; 
     192        vec &_mu = est._mu(); 
     193        chmat &_P = est._R(); 
    192194        chmat &_Ry = fy._R(); 
    193195        vec &_yp = fy._mu(); 
    … …  
    211213 
    212214        ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 ); 
    213         _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimx - 1 ) ).T(); 
    214         ( _P._Ch() ) = postA ( dimy, dimy + dimx - 1, dimy, dimy + dimx - 1 ); 
     215        _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dimension() - 1 ) ).T(); 
     216        ( _P._Ch() ) = postA ( dimy, dimy + dimension() - 1, dimy, dimy + dimension() - 1 ); 
    215217 
    216218//      mat iRY = inv(_Ry->to_mat()); 
    … …  
    268270        //connect 
    269271        shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory ); 
    270         set_drv ( *drv ); 
     272        set_yrv ( *drv ); 
    271273        shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory ); 
    272274        set_rv ( *rv ); 
    … …  
    282284 
    283285        set_parameters ( A ); 
    284         set_drv ( A ( 0 )->_drv() ); 
     286        set_yrv ( A ( 0 )->_yrv() ); 
    285287        //set_rv(A(0)->_rv()); 
    286288 

  • library/bdm/estim/kalman.h

    r675 r679  
    3434{ 
    3535        protected: 
    36                 //! cache of rv.count() 
    37                 int dimx; 
    38                 //! cache of rvy.count() 
    39                 int dimy; 
    40                 //! cache of rvu.count() 
    41                 int dimu; 
    4236                //! Matrix A 
    4337                mat A; 
    … …  
    5347                sq_T R; 
    5448        public: 
    55                 StateSpace() : dimx (0), dimy (0), dimu (0), A(), B(), C(), D(), Q(), R() {} 
     49                StateSpace() : A(), B(), C(), D(), Q(), R() {} 
    5650                //!copy constructor 
    57                 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) {} 
     51                StateSpace(const StateSpace<sq_T> &S0) : A(S0.A), B(S0.B), C(S0.C), D(S0.D), Q(S0.Q), R(S0.R) {} 
    5852                //! set all matrix parameters 
    5953                void set_parameters (const mat &A0, const  mat &B0, const  mat &C0, const  mat &D0, const  sq_T &Q0, const sq_T &R0); 
    … …  
    8377                }                
    8478                //! access function 
    85                 int _dimx(){return dimx;} 
    86                 //! access function 
    87                 int _dimy(){return dimy;} 
    88                 //! access function 
    89                 int _dimu(){return dimu;} 
    90                 //! access function 
    9179                const mat& _A() const {return A;} 
    9280                //! access function 
    … …  
    10997                //! id of output 
    11098                RV yrv; 
    111                 //! id of input 
    112                 RV urv; 
    11399                //! Kalman gain 
    114100                mat  _K; 
    115101                //!posterior 
    116                 shared_ptr<enorm<sq_T> > est; 
     102                enorm<sq_T> est; 
    117103                //!marginal on data f(y|y) 
    118104                enorm<sq_T>  fy; 
    119105        public: 
    120                 Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(),urv(), _K(),  est(new enorm<sq_T>){} 
     106                Kalman<sq_T>() : BM(), StateSpace<sq_T>(), yrv(), _K(),  est(){} 
    121107                //! Copy constructor 
    122                 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){} 
     108                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){} 
    123109                //!set statistics of the posterior 
    124                 void set_statistics (const vec &mu0, const mat &P0) {est->set_parameters (mu0, P0); }; 
     110                void set_statistics (const vec &mu0, const mat &P0) {est.set_parameters (mu0, P0); }; 
    125111                //!set statistics of the posterior 
    126                 void set_statistics (const vec &mu0, const sq_T &P0) {est->set_parameters (mu0, P0); }; 
     112                void set_statistics (const vec &mu0, const sq_T &P0) {est.set_parameters (mu0, P0); }; 
    127113                //! return correctly typed posterior (covariant return) 
    128                 const enorm<sq_T>& posterior() const {return *est.get();} 
    129                 //! shared posterior 
    130                 shared_ptr<epdf> shared_posterior() {return est;} 
     114                const enorm<sq_T>& posterior() const {return est;} 
    131115                //! load basic elements of Kalman from structure 
    132116                void from_setting (const Setting &set) { 
    … …  
    139123                        // Initial values 
    140124                        UI::get (yrv, set, "yrv", UI::optional); 
    141                         UI::get (urv, set, "urv", UI::optional); 
    142                         set_drv(concat(yrv,urv)); 
     125                        UI::get (rvc, set, "urv", UI::optional); 
     126                        set_yrv(concat(yrv,rvc)); 
    143127                         
    144128                        validate(); 
    … …  
    147131                void validate() { 
    148132                        StateSpace<sq_T>::validate(); 
    149                         bdm_assert(est->dimension(), "Statistics and model parameters mismatch"); 
     133                        dimy = this->C.rows(); 
     134                        dimc = this->B.cols(); 
     135                        set_dim(this->A.rows()); 
     136                         
     137                        bdm_assert(est.dimension(), "Statistics and model parameters mismatch"); 
    150138                } 
    151139}; 
    … …  
    160148                KalmanFull() :Kalman<fsqmat>(){}; 
    161149                //! Here dt = [yt;ut] of appropriate dimensions 
    162                 void bayes (const vec &dt); 
     150                void bayes (const vec &yt, const vec &cond=empty_vec); 
    163151                BM* _copy_() const { 
    164152                        KalmanFull* K = new KalmanFull; 
    165153                        K->set_parameters (A, B, C, D, Q, R); 
    166                         K->set_statistics (est->_mu(), est->_R()); 
     154                        K->set_statistics (est._mu(), est._R()); 
    167155                        return K; 
    168156                } 
    … …  
    192180                        KalmanCh* K = new KalmanCh; 
    193181                        K->set_parameters (A, B, C, D, Q, R); 
    194                         K->set_statistics (est->_mu(), est->_R()); 
     182                        K->set_statistics (est._mu(), est._R()); 
    195183                        return K; 
    196184                } 
    … …  
    213201                Thus this object evaluates only predictors! Not filtering densities. 
    214202                */ 
    215                 void bayes (const vec &dt); 
     203                void bayes (const vec &yt, const vec &cond=empty_vec); 
    216204 
    217205                void from_setting(const Setting &set){ 
    218206                        Kalman<chmat>::from_setting(set); 
     207                        validate(); 
     208                } 
     209                void validate() { 
     210                        Kalman<chmat>::validate(); 
    219211                        initialize(); 
    220212                } 
    … …  
    244236 
    245237                //! Here dt = [yt;ut] of appropriate dimensions 
    246                 void bayes (const vec &dt); 
     238                void bayes (const vec &yt, const vec &cond=empty_vec); 
    247239                //! set estimates 
    248240                void set_statistics (const vec &mu0, const mat &P0) { 
    249                         est->set_parameters (mu0, P0); 
     241                        est.set_parameters (mu0, P0); 
    250242                }; 
    251243                //! access function 
    252244                const mat _R() { 
    253                         return est->_R().to_mat(); 
     245                        return est._R().to_mat(); 
    254246                } 
    255247                void from_setting (const Setting &set) { 
    … …  
    280272                        //connect 
    281273                        shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory ); 
    282                         set_drv ( *drv ); 
     274                        set_yrv ( *drv ); 
    283275                        shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory ); 
    284276                        set_rv ( *rv ); 
    … …  
    337329                        EKFCh* E = new EKFCh; 
    338330                        E->set_parameters (pfxu, phxu, Q, R); 
    339                         E->set_statistics (est->_mu(), est->_R()); 
     331                        E->set_statistics (est._mu(), est._R()); 
    340332                        return E; 
    341333                } 
    … …  
    344336 
    345337                //! Here dt = [yt;ut] of appropriate dimensions 
    346                 void bayes (const vec &dt); 
     338                void bayes (const vec &yt, const vec &cond=empty_vec); 
    347339 
    348340                void from_setting (const Setting &set); 
    … …  
    389381                        est.set_parameters (A (0)->posterior().mean(), A (0)->posterior()._R()); 
    390382                } 
    391                 void bayes (const vec &dt) { 
     383                void bayes (const vec &yt, const vec &cond=empty_vec) { 
    392384                        int n = Models.length(); 
    393385                        int i; 
    394386                        for (i = 0; i < n; i++) { 
    395                                 Models (i)->bayes (dt); 
     387                                Models (i)->bayes (yt); 
    396388                                _lls (i) = Models (i)->_ll(); 
    397389                        } 
    … …  
    475467                                Crv.add(urv.copy_t(t)); 
    476468                        } 
    477                          
    478                         this->dimx = xrv._dsize(); 
    479                         this->dimy = yrv._dsize(); 
    480                         this->dimu = urv._dsize(); 
    481                          
     469                                                 
    482470                        // get mapp 
    483471                        th2A.set_connection(xrv, ml._rvc()); 
    … …  
    486474 
    487475                        //set matrix sizes 
    488                         this->A=zeros(dimx,dimx); 
    489                         for (int j=1; j<dimx; j++){A(j,j-1)=1.0;} // off diagonal 
    490                                 this->B=zeros(dimx,1); 
     476                        this->A=zeros(xrv._dsize(),xrv._dsize()); 
     477                        for (int j=1; j<xrv._dsize(); j++){A(j,j-1)=1.0;} // off diagonal 
     478                                this->B=zeros(xrv._dsize(),1); 
    491479                                this->B(0) = 1.0; 
    492                                 this->C=zeros(1,dimx); 
     480                                this->C=zeros(1,xrv._dsize()); 
    493481                                this->D=zeros(1,urv._dsize()); 
    494                                 this->Q = zeros(dimx,dimx); 
     482                                this->Q = zeros(xrv._dsize(),xrv._dsize()); 
    495483                        // R is set by update 
    496484                         
    … …  
    538526template<class sq_T> 
    539527void StateSpace<sq_T>::validate(){ 
    540         dimx = A.rows(); 
    541         dimu = B.cols(); 
    542         dimy = C.rows(); 
    543         bdm_assert (A.cols() == dimx, "KalmanFull: A is not square"); 
    544         bdm_assert (B.rows() == dimx, "KalmanFull: B is not compatible"); 
    545         bdm_assert (C.cols() == dimx, "KalmanFull: C is not square"); 
    546         bdm_assert ( (D.rows() == dimy) && (D.cols() == dimu), "KalmanFull: D is not compatible"); 
    547         bdm_assert ( (Q.cols() == dimx) && (Q.rows() == dimx), "KalmanFull: Q is not compatible"); 
    548         bdm_assert ( (R.cols() == dimy) && (R.rows() == dimy), "KalmanFull: R is not compatible"); 
     528        bdm_assert (A.cols() == A.rows(), "KalmanFull: A is not square"); 
     529        bdm_assert (B.rows() == A.rows(), "KalmanFull: B is not compatible"); 
     530        bdm_assert (C.cols() == A.rows(), "KalmanFull: C is not compatible"); 
     531        bdm_assert ( (D.rows() == C.rows()) && (D.cols() == B.cols()), "KalmanFull: D is not compatible"); 
     532        bdm_assert ( (Q.cols() == A.rows()) && (Q.rows() == A.rows()), "KalmanFull: Q is not compatible"); 
     533        bdm_assert ( (R.cols() == C.rows()) && (R.rows() == C.rows()), "KalmanFull: R is not compatible"); 
    549534} 
    550535 

  • library/bdm/estim/mixtures.cpp

    r565 r679  
    3030                //pick one datum 
    3131                int ind = floor ( ndat * UniRNG.sample() ); 
    32                 Coms ( i )->bayes ( Data.get_col ( ind ), 1.0 ); 
     32                Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec ); 
    3333                //flatten back to oringinal 
    3434                Coms ( i )->flatten ( Com0 ); 
    … …  
    4141} 
    4242 
    43 void MixEF::bayesB ( const mat &data , const vec &wData ) { 
     43void MixEF::bayesB ( const mat &data , const mat &cond, const vec &wData ) { 
    4444        int ndat = data.cols(); 
    4545        int t, i, niter; 
    … …  
    106106                        //#pragma omp parallel for 
    107107                        for ( i = 0; i < n; i++ ) { 
    108                                 Coms ( i )-> bayes ( data.get_col ( t ), W ( i, t ) * wData ( t ) ); 
     108                                Coms ( i )-> bayes_weighted ( data.get_col ( t ), empty_vec, W ( i, t ) * wData ( t ) ); 
    109109                        } 
    110110                        weights.bayes ( W.get_col ( t ) * wData ( t ) ); 
    … …  
    122122} 
    123123 
    124 void MixEF::bayes ( const vec &data ) { 
     124void MixEF::bayes ( const vec &data, const vec &cond=empty_vec ) { 
    125125 
    126126}; 
    127127 
    128 void MixEF::bayes ( const mat &data ) { 
    129         this->bayesB ( data, ones ( data.cols() ) ); 
     128void MixEF::bayes ( const mat &data, const vec &cond=empty_vec ) { 
     129        this->bayesB ( data, cond, ones ( data.cols() ) ); 
    130130}; 
    131131 
    132132 
    133 double MixEF::logpred ( const vec &dt ) const { 
     133double MixEF::logpred ( const vec &yt ) const { 
    134134 
    135135        vec w = weights.posterior().mean(); 
    136136        double exLL = 0.0; 
    137137        for ( int i = 0; i < n; i++ ) { 
    138                 exLL += w ( i ) * exp ( Coms ( i )->logpred ( dt ) ); 
     138                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt ) ); 
    139139        } 
    140140        return log ( exLL ); 

  • library/bdm/estim/mixtures.h

    r660 r679  
    108108        } 
    109109        //! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006 
    110         void bayes ( const vec &dt ); 
     110        void bayes ( const vec &yt, const vec &cond ); 
    111111        //! EM algorithm 
    112         void bayes ( const mat &dt ); 
     112        void bayes ( const mat &yt, const vec &cond ); 
    113113        //! batch weighted Bayes rule  
    114         void bayesB ( const mat &dt, const vec &wData ); 
    115         double logpred ( const vec &dt ) const; 
     114        void bayesB ( const mat &yt, const mat &cond, const vec &wData ); 
     115        double logpred ( const vec &yt ) const; 
    116116        //! return correctly typed posterior (covariant return) 
    117117        const eprod& posterior() const { 

  • library/bdm/estim/particles.cpp

    r665 r679  
    4646} 
    4747 
    48 void PF::bayes ( const vec &dt ) { 
    49         vec yt=dt.left(obs->dimension()); 
    50         vec ut=dt.get(obs->dimension(),dt.length()-1); 
     48void PF::bayes( const vec &yt, const vec &cond ) { 
     49        const vec &ut=cond; //todo 
    5150         
    5251        int i; 
    … …  
    7473// } 
    7574 
    76 void MPF::bayes ( const vec &dt ) {      
     75void MPF::bayes ( const vec &yt, const vec &cond ) {     
    7776        // follows PF::bayes in most places!!!   
    7877        int i; 
    … …  
    8584        #pragma parallel for 
    8685        for ( i = 0; i < n; i++ ) { 
    87                 BMs(i) -> condition(pf->posterior()._sample(i)); 
    88                 BMs(i) -> bayes(dt); 
     86                BMs(i) -> bayes(this2bm.pushdown(yt), this2bm.get_cond(yt,cond)); 
    8987                lls ( i ) += BMs(i)->_ll(); 
    9088        } 

  • library/bdm/estim/particles.h

    r676 r679  
    9393                return eff < ( res_threshold*n ); 
    9494        } 
    95         void bayes ( const vec &dt ); 
     95        void bayes ( const vec &yt, const vec &cond ); 
    9696        //!access function 
    9797        vec& __w() { return _w; } 
    … …  
    131131                u.add(obs_u); // join both u, and check if they do not overlap 
    132132 
    133                 set_drv(concat(obs->_rv(),u) ); 
     133                set_yrv(concat(obs->_rv(),u) ); 
    134134        } 
    135135        //! auxiliary function reading parameter 'resmethod' from configuration file 
    … …  
    296296        mpfepdf jest; 
    297297 
    298         //! Log means of BMs 
    299         bool opt_L_mea; 
    300  
     298        //! datalink from global yt and cond (Up) to BMs yt and cond (Down) 
     299        datalink_m2m this2bm; 
     300        //! datalink from global yt and cond (Up) to PFs yt and cond (Down) 
     301        datalink_m2m this2pf; 
     302         
    301303public: 
    302304        //! Default constructor. 
    … …  
    320322        }; 
    321323 
    322         void bayes ( const vec &dt ); 
     324        void bayes ( const vec &yt, const vec &cond ); 
    323325        const epdf& posterior() const { 
    324326                return jest; 
    … …  
    328330        void set_options ( const string &opt ) { 
    329331                BM::set_options(opt); 
    330                 opt_L_mea = ( opt.find ( "logmeans" ) != string::npos ); 
    331332        } 
    332333 
    … …  
    352353 
    353354                pf = new PF; 
    354                 // rpior must be set before BM 
     355                // prior must be set before BM 
    355356                pf->prior_from_set(set); 
    356357                pf->resmethod_from_set(set); 
    357358                pf->set_model(par,obs); 
    358                  
     359                                 
    359360                shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory); 
    360361                set_BM(*BM0); 
    … …  
    367368                //find potential input - what remains in rvc when we subtract rv 
    368369                RV u = par->_rvc().remove_time().subt( par->_rv() );             
    369                 set_drv(concat(BM0->_drv(),u) ); 
     370                set_yrv(concat(BM0->_yrv(),u) ); 
    370371                validate(); 
    371372        } 
    … …  
    377378                } 
    378379                jest.read_parameters(); 
    379                 for ( int i = 0; i < pf->__w().length(); i++ ) { 
    380                         BMs ( i )->condition ( pf->posterior()._sample ( i ) ); 
    381                 } 
     380                this2bm.set_connection(BMs(0)->_yrv(), BMs(0)->_rvc(), yrv, rvc); 
     381                this2bm.set_connection(pf->_yrv(), pf->_rvc(), yrv, rvc); 
    382382        } 
    383383 

  • library/bdm/itpp_ext.cpp

    r662 r679  
    2626namespace itpp 
    2727{ 
     28        vec empty_vec = vec(0); 
     29         
    2830Array<int> to_Arr (const ivec &indices) 
    2931{ 

  • library/bdm/itpp_ext.h

    r662 r679  
    1919 
    2020namespace itpp { 
     21        extern vec empty_vec; 
     22         
    2123Array<int> to_Arr ( const ivec &indices ); 
    2224ivec linspace ( int from, int to ); 

  • library/bdm/mex/mex_BM.h

    r660 r679  
    7575                        est.from_setting(posterior); 
    7676                } 
    77                 void bayes(const vec &dt)  { 
     77                void bayes(const vec &yt, const vec &cond)  { 
    7878                //void bayes()  { 
    7979                        mxArray *tmp, *old; 

  • library/bdm/stat/exp_family.cpp

    r678 r679  
    1414        /////////// 
    1515 
    16 void BMEF::bayes ( const vec &dt ) { 
    17         this->bayes ( dt, 1.0 ); 
     16void BMEF::bayes( const vec &yt, const vec &cond ) { 
     17        this->bayes_weighted ( yt, cond, 1.0 ); 
    1818}; 
    1919 

  • library/bdm/stat/exp_family.h

    r678 r679  
    9696 
    9797                //! Weighted update of sufficient statistics (Bayes rule) 
    98                 virtual void bayes (const vec &data, const double w) {}; 
     98                virtual void bayes_weighted (const vec &data, const vec &cond=empty_vec, const double w=1.0) {}; 
    9999                //original Bayes 
    100                 void bayes (const vec &dt); 
     100                void bayes (const vec &yt, const vec &cond=empty_vec); 
    101101 
    102102                //!Flatten the posterior according to the given BMEF (of the same type!) 
    … …  
    436436                //! Sets sufficient statistics to match that of givefrom mB0 
    437437                void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;} 
    438                 void bayes (const vec &dt) { 
     438                void bayes (const vec &yt, const vec &cond=empty_vec) { 
    439439                        if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();} 
    440                         beta += dt; 
     440                        beta += yt; 
    441441                        if (evalll) {ll = est.lognc() - last_lognc;} 
    442442                } 
    443                 double logpred (const vec &dt) const { 
     443                double logpred (const vec &yt) const { 
    444444                        eDirich pred (est); 
    445445                        vec &beta = pred._beta(); 
    … …  
    452452                                else{lll = pred.lognc();} 
    453453 
    454                         beta += dt; 
     454                        beta += yt; 
    455455                        return pred.lognc() - lll; 
    456456                } 
    … …  
    832832                } 
    833833                void condition (const vec &cond) { 
    834                         iepdf._mu() = A * cond + mu_const; 
     834                        if (cond.length()>0) { 
     835                                iepdf._mu() = A * cond + mu_const; 
     836                        } else { 
     837                                iepdf._mu() =  mu_const; 
     838                        } 
    835839                        double zeta; 
    836840                        //ugly hack! 

  • library/bdm/stat/merger.cpp

    r565 r679  
    103103                //Re-Initialize Mixture model 
    104104                Mix.flatten ( &Mix_init ); 
    105                 Mix.bayesB ( Smp_ex, w*Npoints ); 
     105                Mix.bayesB ( Smp_ex,empty_vec, w*Npoints ); 
    106106                delete Mpred; 
    107107                Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!! 

  • library/bdm/stat/merger.h

    r569 r679  
    342342        } 
    343343        //! loglikelihood computed on mixture models 
    344         double evallog ( const vec &dt ) const { 
    345                 vec dtf = ones ( dt.length() + 1 ); 
    346                 dtf.set_subvector ( 0, dt ); 
     344        double evallog ( const vec &yt ) const { 
     345                vec dtf = ones ( yt.length() + 1 ); 
     346                dtf.set_subvector ( 0, yt ); 
    347347                return Mix.logpred ( dtf ); 
    348348        } 

  • library/tests/arx_elem_test.cpp

    r536 r679  
    99        ARX Ar; 
    1010        Ar.set_statistics ( 1, V0, -1.0 ); 
     11        Ar.set_constant(true); 
     12        Ar.validate(); 
    1113 
    1214        mat mu ( 1, 1 ); 
    … …  
    3638 
    3739        mlstudent* Ap = Ar.predictor_student(); 
    38         vec Ap_x = Ap->evallogcond_m ( X, vec_1 ( 1.0 ) ); 
     40        vec Ap_x = Ap->evallogcond_m ( X, empty_vec ); 
    3941        vec ll_x = Ar.logpred_m ( X2 ); 
    4042 

  • library/tests/arx_test.cpp

    r477 r679  
    2727        ARX Ar; 
    2828        Ar.set_statistics ( 1, V0, nu0 );               // Estimator 
     29        Ar.set_constant(false); 
     30        Ar.validate(); 
    2931        const epdf& f_thr = Ar.posterior();          // refrence to posterior of the estimator 
    3032 
    … …  
    3436        Yt.set_subvector ( 0, randn ( ord ) ); //initial values 
    3537        vec rgr ( ord );                // regressor 
    36         vec Psi ( ord + 1 );            // extended regressor 
    3738 
    3839        //print moments of the prior distribution 
    … …  
    4950                Yt ( t ) = th * rgr + sqr * NorRNG(); 
    5051 
    51                 Psi = concat ( Yt ( t ), rgr ); // Inefficient! Used for compatibility with Matlab! 
    52                 Ar.bayes ( Psi );             // Bayes rule 
     52                Ar.bayes ( vec_1(Yt(t)), rgr );             // Bayes rule 
    5353 
    5454                // Build predictor 

  • library/tests/dirfile-format.matrix

    r425 r679  
    1 r_0 RAW d 1 
    2 r_1 RAW d 1 
    3 th_alog RAW d 1 
    4 th_blog RAW d 1 
     1r.0 RAW d 1 
     2r.1 RAW d 1 
     3th.alog RAW d 1 
     4th.blog RAW d 1 

  • library/tests/mixtures_test.cpp

    r529 r679  
    9999        // Add ones for constant coefficients 
    100100        mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) ); 
    101         Post.bayes ( Data ); 
     101        Post.bayes ( Data , empty_vec); 
    102102 
    103103        cout << "Posterior mixture:" << endl; 
    … …  
    115115        mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );; 
    116116 
    117         RND.bayes ( Data ); 
     117        RND.bayes ( Data, empty_vec ); 
    118118        cout << endl << "== BAYES ==" << endl; 
    119119        disp_mix2D ( RND ); 

  • library/tests/tutorial/arx_simple.cpp

    r477 r679  
    99        ARX Ar; 
    1010        Ar.set_statistics ( 1, V0 ); //nu is default (set to have finite moments) 
     11        Ar.set_constant(true); 
     12        Ar.validate(); 
    1113        // forgetting is default: 1.0 
    1214        mat Data = concat_vertical ( randn ( 1, 100 ), ones ( 1, 100 ) ); 

  • library/tests/tutorial/kalman_simple.cpp

    r477 r679  
    2020        KF.set_parameters ( A, B, C, D,/*covariances*/ Q, R ); 
    2121        KF.set_statistics ( mu0, P0 ); 
     22        KF.validate(); 
    2223        // Estimation loop 
    2324        for ( int i = 0; i < 100; i++ ) { 
    24                 KF.bayes ( randn ( dx + du ) ); 
     25                KF.bayes ( randn ( dy), randn( du ) ); 
    2526        } 
    2627        //print results