Changeset 679 for library/bdm/estim
- Timestamp:
- 10/23/09 00:05:25 (15 years ago)
- Location:
- library/bdm/estim
- Files:
-
- 8 modified
Legend:
- Unmodified
- Added
- Removed
-
library/bdm/estim/arx.cpp
r665 r679 2 2 namespace bdm { 3 3 4 void ARX::bayes ( const vec &dt, const double w ) {4 void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) { 5 5 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 7 14 if ( frg < 1.0 ) { 8 15 est.pow ( frg ); // multiply V and nu 16 9 17 10 18 //stabilize 11 19 ldmat V0=alter_est._V(); //$ copy 12 20 double &nu0=alter_est._nu(); 21 13 22 V0*=(1-frg); 14 23 V += V0; //stabilization … … 20 29 } 21 30 } 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 ); 28 32 nu += w; 29 33 … … 36 40 } 37 41 38 double ARX::logpred ( const vec & dt ) const {42 double ARX::logpred ( const vec &yt ) const { 39 43 egiw pred ( est ); 40 44 ldmat &V = pred._V(); … … 42 46 43 47 double lll; 48 vec dyad_p = dyad; 49 dyad_p.set_subvector(0,yt); 44 50 45 51 if ( frg < 1.0 ) { … … 53 59 } 54 60 55 V.opupdt ( d t, 1.0 );61 V.opupdt ( dyad_p, 1.0 ); 56 62 nu += 1.0; 57 63 // log(sqrt(2*pi)) = 0.91893853320467 … … 67 73 const ARX* A0 = dynamic_cast<const ARX*> ( B0 ); 68 74 69 bdm_assert_debug ( V.rows() == A0->V.rows(), "ARX::set_statistics Statistics differ" );70 set_statistics ( A0->dim x, 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() ); 71 77 } 72 78 73 79 enorm<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 ); 77 82 78 83 enorm<ldmat>* tmp; 79 84 tmp = new enorm<ldmat> ( ); 80 85 //TODO: too hackish 81 if ( drv._dsize() > 0 ) {86 if ( yrv._dsize() > 0 ) { 82 87 } 83 88 … … 91 96 92 97 mlnorm<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 ); 97 100 mlnorm<ldmat>* tmp; 98 101 tmp = new mlnorm<ldmat> ( ); … … 106 109 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) ); 107 110 } else { 108 tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );111 tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ) ); 109 112 } 110 113 return tmp; … … 112 115 113 116 mlstudent* 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 ); 118 121 mlstudent* tmp; 119 122 tmp = new mlstudent ( ); … … 122 125 mu = mu.T(); 123 126 124 int xdim = dimx;125 127 int end = V._L().rows() - 1; 126 ldmat Lam ( V._L() ( xdim, end, xdim, end ), V._D() ( xdim, end ) ); //exp val of R128 ldmat Lam ( V._L() ( dimy, end, dimy, end ), V._D() ( dimy, end ) ); //exp val of R 127 129 128 130 … … 132 134 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam ); 133 135 } 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 ); 135 137 } 136 138 } else { 137 139 // no constant term 138 tmp->set_parameters ( mu, zeros ( xdim), ldmat ( R ), Lam );140 tmp->set_parameters ( mu, zeros ( dimy ), ldmat ( R ), Lam ); 139 141 } 140 142 return tmp; … … 214 216 215 217 void 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 ); 217 219 shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory ); 218 int ylen = yrv->_dsize();220 dimy = yrv_->_dsize(); 219 221 // 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; 224 226 225 227 string opt; … … 233 235 have_constant=constant>0; 234 236 } 235 i f (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);}237 int rgrlen = dimc+int(have_constant==true); 236 238 237 239 //init 238 240 shared_ptr<egiw> pri=UI::build<egiw>(set, "prior", UI::optional); 239 241 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"); 242 244 est.set_parameters( pri->_dimx(),pri->_V(), pri->_nu()); 243 245 }else{ 244 est.set_parameters( ylen, zeros(ylen+rgrlen));246 est.set_parameters( dimy, zeros(dimy+rgrlen)); 245 247 set_prior_default(est); 246 248 } … … 248 250 shared_ptr<egiw> alt=UI::build<egiw>(set, "alternative", UI::optional); 249 251 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"); 252 254 alter_est.set_parameters( alt->_dimx(),alt->_V(), alt->_nu()); 253 255 } else { … … 264 266 shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional ); 265 267 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 ) ) ); 267 269 } else { 268 270 est.set_rv ( *rv_par ); … … 270 272 validate(); 271 273 } 272 273 } 274 } 275 -
library/bdm/estim/arx.h
r665 r679 43 43 class ARX: public BMEF { 44 44 protected: 45 //!size of output variable (needed in regressors)46 int dimx;47 //!description of modelled data \f$ y_t \f$ in the likelihood function48 //! Do NOT access directly, only via \c get_yrv().49 RV _yrv;50 //! rv of regressor51 RV rgrrv;52 //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density53 egiw est;54 //! cached value of est.V55 ldmat &V;56 //! cached value of est.nu57 double ν58 45 //! switch if constant is modelled or not 59 46 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; 62 51 //! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy] 63 52 egiw alter_est; … … 65 54 //! \name Constructors 66 55 //!@{ 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) { }; 74 58 ARX* _copy_() const; 75 59 void set_parameters ( double frg0 ) { … … 79 63 have_constant=const0; 80 64 } 81 void set_statistics ( int dim x0, const ldmat V0, double nu0 = -1.0 ) {82 est.set_parameters ( dim x0, V0, nu0 );65 void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) { 66 est.set_parameters ( dimy0, V0, nu0 ); 83 67 last_lognc = est.lognc(); 84 dim x = dimx0;68 dimy = dimy0; 85 69 } 86 70 //!@} … … 93 77 94 78 //! 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 ); 98 82 }; 99 double logpred ( const vec & dt ) const;83 double logpred ( const vec &yt ) const; 100 84 void flatten ( const BMEF* B ) { 101 85 const ARX* A = dynamic_cast<const ARX*> ( B ); 102 86 // nu should be equal to B.nu 103 est.pow ( A-> nu / nu);87 est.pow ( A->posterior()._nu() / posterior()._nu() ); 104 88 if ( evalll ) { 105 89 last_lognc = est.lognc(); … … 110 94 //! Predictor for empty regressor 111 95 enorm<ldmat>* epredictor() const { 112 bdm_assert_debug ( dim x == V.rows() - 1, "Regressor is not only 1" );96 bdm_assert_debug ( dimy == posterior()._V().rows() - 1, "Regressor is not only 1" ); 113 97 return epredictor ( vec_1 ( 1.0 ) ); 114 98 } … … 127 111 const egiw& posterior() const { 128 112 return est; 129 }130 //!@}131 132 //!\name Connection133 //!@{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 it142 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 now150 bdm_assert_debug ( _yrv._dsize() == dimx, "incompatible drv" );151 return _yrv;152 113 } 153 114 //!@} … … 174 135 175 136 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 178 148 } 179 149 //! function sets prior and alternative density -
library/bdm/estim/kalman.cpp
r653 r679 8 8 9 9 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(); 10 void 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(); 17 18 mat& _Ry = fy._R(); 18 19 vec& yp = fy._mu(); … … 42 43 phxu = phxu0; 43 44 44 dimx = pfxu0->_dimx();45 set_dim( pfxu0->_dimx()); 45 46 dimy = phxu0->dimension(); 46 dim u= pfxu0->_dimu();47 est ->set_parameters( zeros(dimx), eye(dimx) );48 49 A.set_size ( dim x, dimx);50 C.set_size ( dimy, dim x);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() ); 51 52 //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 ); 53 54 B.clear(); 54 phxu->dfdx_cond ( est ->_mu(), zeros ( dimu), C, true );55 phxu->dfdx_cond ( est._mu(), zeros ( dimc ), C, true ); 55 56 D.clear(); 56 57 … … 59 60 } 60 61 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(); 62 void 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(); 68 70 mat& _Ry = fy._R(); 69 71 vec& yp = fy._mu(); 70 72 71 pfxu->dfdx_cond ( mu, zeros ( dim u), A, true );72 phxu->dfdx_cond ( mu, zeros ( dim u), C, true );73 pfxu->dfdx_cond ( mu, zeros ( dimc ), A, true ); 74 phxu->dfdx_cond ( mu, zeros ( dimc ), C, true ); 73 75 74 76 //Time update … … 94 96 ( ( StateSpace<chmat>* ) this )->set_parameters ( A0, B0, C0, D0, Q0, R0 ); 95 97 96 _K=zeros(dimx,dimy); 97 // Cholesky special! 98 initialize(); 98 _K=zeros(dimension(),dimy); 99 99 } 100 100 101 101 void KalmanCh::initialize(){ 102 preA = zeros ( dimy + dim x + dimx, dimy + dimx);102 preA = zeros ( dimy + dimension() + dimension(), dimy + dimension() ); 103 103 // preA.clear(); 104 104 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 108 void 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; 112 114 vec pom ( dimy ); 113 115 114 chmat &_P=est ->_R();115 vec &_mu = est ->_mu();116 mat _K(dim x,dimy);116 chmat &_P=est._R(); 117 vec &_mu = est._mu(); 118 mat _K(dimension(),dimy); 117 119 chmat &_Ry=fy._R(); 118 120 vec &_yp = fy._mu(); … … 130 132 131 133 ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 ); 132 _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dim x- 1 ) ).T();133 ( _P._Ch() ) = postA ( dimy, dimy + dim x - 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 ); 134 136 135 137 _mu = A * ( _mu ) + B * u; … … 154 156 phxu = phxu0; 155 157 156 dimx = pfxu0->dimension();158 set_dim( pfxu0->_dimx()); 157 159 dimy = phxu0->dimension(); 158 dim u= pfxu0->_dimu();159 160 vec &_mu = est ->_mu();160 dimc = pfxu0->_dimu(); 161 162 vec &_mu = est._mu(); 161 163 // if mu is not set, set it to zeros, just for constant terms of A and C 162 if ( _mu.length() != dim x ) _mu = zeros ( dimx);163 A = zeros ( dim x, dimx);164 C = zeros ( dimy, dim x);165 preA = zeros ( dimy + dim x + 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() ); 166 168 167 169 //initialize matrices A C, later, these will be only updated! 168 pfxu->dfdx_cond ( _mu, zeros ( dim u), A, true );170 pfxu->dfdx_cond ( _mu, zeros ( dimc ), A, true ); 169 171 // pfxu->dfdu_cond ( *_mu,zeros ( dimu ),B,true ); 170 172 B.clear(); 171 phxu->dfdx_cond ( _mu, zeros ( dim u), C, true );173 phxu->dfdx_cond ( _mu, zeros ( dimc ), C, true ); 172 174 // phxu->dfdu_cond ( *_mu,zeros ( dimu ),D,true ); 173 175 D.clear(); … … 179 181 preA.clear(); 180 182 preA.set_submatrix ( 0, 0, R._Ch() ); 181 preA.set_submatrix ( dimy + dim x, 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 187 void EKFCh::bayes ( const vec &yt, const vec &cond ) { 186 188 187 189 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(); 192 194 chmat &_Ry = fy._R(); 193 195 vec &_yp = fy._mu(); … … 211 213 212 214 ( _Ry._Ch() ) = postA ( 0, dimy - 1, 0, dimy - 1 ); 213 _K = inv ( A ) * ( postA ( 0, dimy - 1 , dimy, dimy + dim x- 1 ) ).T();214 ( _P._Ch() ) = postA ( dimy, dimy + dim x - 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 ); 215 217 216 218 // mat iRY = inv(_Ry->to_mat()); … … 268 270 //connect 269 271 shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory ); 270 set_ drv ( *drv );272 set_yrv ( *drv ); 271 273 shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory ); 272 274 set_rv ( *rv ); … … 282 284 283 285 set_parameters ( A ); 284 set_ drv ( A ( 0 )->_drv() );286 set_yrv ( A ( 0 )->_yrv() ); 285 287 //set_rv(A(0)->_rv()); 286 288 -
library/bdm/estim/kalman.h
r675 r679 34 34 { 35 35 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;42 36 //! Matrix A 43 37 mat A; … … 53 47 sq_T R; 54 48 public: 55 StateSpace() : dimx (0), dimy (0), dimu (0),A(), B(), C(), D(), Q(), R() {}49 StateSpace() : A(), B(), C(), D(), Q(), R() {} 56 50 //!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) {} 58 52 //! set all matrix parameters 59 53 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 } 84 78 //! access function 85 int _dimx(){return dimx;}86 //! access function87 int _dimy(){return dimy;}88 //! access function89 int _dimu(){return dimu;}90 //! access function91 79 const mat& _A() const {return A;} 92 80 //! access function … … 109 97 //! id of output 110 98 RV yrv; 111 //! id of input112 RV urv;113 99 //! Kalman gain 114 100 mat _K; 115 101 //!posterior 116 shared_ptr<enorm<sq_T>> est;102 enorm<sq_T> est; 117 103 //!marginal on data f(y|y) 118 104 enorm<sq_T> fy; 119 105 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(){} 121 107 //! 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){} 123 109 //!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); }; 125 111 //!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); }; 127 113 //! 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;} 131 115 //! load basic elements of Kalman from structure 132 116 void from_setting (const Setting &set) { … … 139 123 // Initial values 140 124 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)); 143 127 144 128 validate(); … … 147 131 void validate() { 148 132 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"); 150 138 } 151 139 }; … … 160 148 KalmanFull() :Kalman<fsqmat>(){}; 161 149 //! Here dt = [yt;ut] of appropriate dimensions 162 void bayes (const vec & dt);150 void bayes (const vec &yt, const vec &cond=empty_vec); 163 151 BM* _copy_() const { 164 152 KalmanFull* K = new KalmanFull; 165 153 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()); 167 155 return K; 168 156 } … … 192 180 KalmanCh* K = new KalmanCh; 193 181 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()); 195 183 return K; 196 184 } … … 213 201 Thus this object evaluates only predictors! Not filtering densities. 214 202 */ 215 void bayes (const vec & dt);203 void bayes (const vec &yt, const vec &cond=empty_vec); 216 204 217 205 void from_setting(const Setting &set){ 218 206 Kalman<chmat>::from_setting(set); 207 validate(); 208 } 209 void validate() { 210 Kalman<chmat>::validate(); 219 211 initialize(); 220 212 } … … 244 236 245 237 //! Here dt = [yt;ut] of appropriate dimensions 246 void bayes (const vec & dt);238 void bayes (const vec &yt, const vec &cond=empty_vec); 247 239 //! set estimates 248 240 void set_statistics (const vec &mu0, const mat &P0) { 249 est ->set_parameters (mu0, P0);241 est.set_parameters (mu0, P0); 250 242 }; 251 243 //! access function 252 244 const mat _R() { 253 return est ->_R().to_mat();245 return est._R().to_mat(); 254 246 } 255 247 void from_setting (const Setting &set) { … … 280 272 //connect 281 273 shared_ptr<RV> drv = UI::build<RV> ( set, "drv", UI::compulsory ); 282 set_ drv ( *drv );274 set_yrv ( *drv ); 283 275 shared_ptr<RV> rv = UI::build<RV> ( set, "rv", UI::compulsory ); 284 276 set_rv ( *rv ); … … 337 329 EKFCh* E = new EKFCh; 338 330 E->set_parameters (pfxu, phxu, Q, R); 339 E->set_statistics (est ->_mu(), est->_R());331 E->set_statistics (est._mu(), est._R()); 340 332 return E; 341 333 } … … 344 336 345 337 //! Here dt = [yt;ut] of appropriate dimensions 346 void bayes (const vec & dt);338 void bayes (const vec &yt, const vec &cond=empty_vec); 347 339 348 340 void from_setting (const Setting &set); … … 389 381 est.set_parameters (A (0)->posterior().mean(), A (0)->posterior()._R()); 390 382 } 391 void bayes (const vec & dt) {383 void bayes (const vec &yt, const vec &cond=empty_vec) { 392 384 int n = Models.length(); 393 385 int i; 394 386 for (i = 0; i < n; i++) { 395 Models (i)->bayes ( dt);387 Models (i)->bayes (yt); 396 388 _lls (i) = Models (i)->_ll(); 397 389 } … … 475 467 Crv.add(urv.copy_t(t)); 476 468 } 477 478 this->dimx = xrv._dsize(); 479 this->dimy = yrv._dsize(); 480 this->dimu = urv._dsize(); 481 469 482 470 // get mapp 483 471 th2A.set_connection(xrv, ml._rvc()); … … 486 474 487 475 //set matrix sizes 488 this->A=zeros( dimx,dimx);489 for (int j=1; j< dimx; j++){A(j,j-1)=1.0;} // off diagonal490 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); 491 479 this->B(0) = 1.0; 492 this->C=zeros(1, dimx);480 this->C=zeros(1,xrv._dsize()); 493 481 this->D=zeros(1,urv._dsize()); 494 this->Q = zeros( dimx,dimx);482 this->Q = zeros(xrv._dsize(),xrv._dsize()); 495 483 // R is set by update 496 484 … … 538 526 template<class sq_T> 539 527 void 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"); 549 534 } 550 535 -
library/bdm/estim/mixtures.cpp
r565 r679 30 30 //pick one datum 31 31 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 ); 33 33 //flatten back to oringinal 34 34 Coms ( i )->flatten ( Com0 ); … … 41 41 } 42 42 43 void MixEF::bayesB ( const mat &data , const vec &wData ) {43 void MixEF::bayesB ( const mat &data , const mat &cond, const vec &wData ) { 44 44 int ndat = data.cols(); 45 45 int t, i, niter; … … 106 106 //#pragma omp parallel for 107 107 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 ) ); 109 109 } 110 110 weights.bayes ( W.get_col ( t ) * wData ( t ) ); … … 122 122 } 123 123 124 void MixEF::bayes ( const vec &data ) {124 void MixEF::bayes ( const vec &data, const vec &cond=empty_vec ) { 125 125 126 126 }; 127 127 128 void MixEF::bayes ( const mat &data ) {129 this->bayesB ( data, ones ( data.cols() ) );128 void MixEF::bayes ( const mat &data, const vec &cond=empty_vec ) { 129 this->bayesB ( data, cond, ones ( data.cols() ) ); 130 130 }; 131 131 132 132 133 double MixEF::logpred ( const vec & dt ) const {133 double MixEF::logpred ( const vec &yt ) const { 134 134 135 135 vec w = weights.posterior().mean(); 136 136 double exLL = 0.0; 137 137 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 ) ); 139 139 } 140 140 return log ( exLL ); -
library/bdm/estim/mixtures.h
r660 r679 108 108 } 109 109 //! 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 ); 111 111 //! EM algorithm 112 void bayes ( const mat & dt);112 void bayes ( const mat &yt, const vec &cond ); 113 113 //! 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; 116 116 //! return correctly typed posterior (covariant return) 117 117 const eprod& posterior() const { -
library/bdm/estim/particles.cpp
r665 r679 46 46 } 47 47 48 void PF::bayes ( const vec &dt ) { 49 vec yt=dt.left(obs->dimension()); 50 vec ut=dt.get(obs->dimension(),dt.length()-1); 48 void PF::bayes( const vec &yt, const vec &cond ) { 49 const vec &ut=cond; //todo 51 50 52 51 int i; … … 74 73 // } 75 74 76 void MPF::bayes ( const vec & dt) {75 void MPF::bayes ( const vec &yt, const vec &cond ) { 77 76 // follows PF::bayes in most places!!! 78 77 int i; … … 85 84 #pragma parallel for 86 85 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)); 89 87 lls ( i ) += BMs(i)->_ll(); 90 88 } -
library/bdm/estim/particles.h
r676 r679 93 93 return eff < ( res_threshold*n ); 94 94 } 95 void bayes ( const vec & dt);95 void bayes ( const vec &yt, const vec &cond ); 96 96 //!access function 97 97 vec& __w() { return _w; } … … 131 131 u.add(obs_u); // join both u, and check if they do not overlap 132 132 133 set_ drv(concat(obs->_rv(),u) );133 set_yrv(concat(obs->_rv(),u) ); 134 134 } 135 135 //! auxiliary function reading parameter 'resmethod' from configuration file … … 296 296 mpfepdf jest; 297 297 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 301 303 public: 302 304 //! Default constructor. … … 320 322 }; 321 323 322 void bayes ( const vec & dt);324 void bayes ( const vec &yt, const vec &cond ); 323 325 const epdf& posterior() const { 324 326 return jest; … … 328 330 void set_options ( const string &opt ) { 329 331 BM::set_options(opt); 330 opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );331 332 } 332 333 … … 352 353 353 354 pf = new PF; 354 // rpior must be set before BM355 // prior must be set before BM 355 356 pf->prior_from_set(set); 356 357 pf->resmethod_from_set(set); 357 358 pf->set_model(par,obs); 358 359 359 360 shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory); 360 361 set_BM(*BM0); … … 367 368 //find potential input - what remains in rvc when we subtract rv 368 369 RV u = par->_rvc().remove_time().subt( par->_rv() ); 369 set_ drv(concat(BM0->_drv(),u) );370 set_yrv(concat(BM0->_yrv(),u) ); 370 371 validate(); 371 372 } … … 377 378 } 378 379 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); 382 382 } 383 383