arx.h
Go to the documentation of this file.00001 00013 #ifndef AR_H 00014 #define AR_H 00015 00016 #include "../math/functions.h" 00017 #include "../stat/exp_family.h" 00018 #include "../base/user_info.h" 00019 //#include "../estim/kalman.h" 00020 #include "arx_straux.h" 00021 00022 namespace bdm { 00023 00043 class ARX: public BMEF { 00044 protected: 00046 bool have_constant; 00048 vec dyad; 00050 RV rgr; 00052 int rgrlen; 00054 egiw est; 00056 egiw alter_est; 00057 public: 00060 ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), have_constant ( true ), dyad(), rgrlen(),est(), alter_est() {}; 00061 ARX ( const ARX &A0 ) : BMEF ( A0 ), have_constant ( A0.have_constant ), dyad ( A0.dyad ),rgrlen(A0.rgrlen), est ( A0.est ), alter_est ( A0.alter_est ) { }; 00062 00063 ARX* _copy() const; 00064 00065 void set_frg ( double frg0 ) { 00066 frg = frg0; 00067 } 00068 void set_constant ( bool const0 ) { 00069 have_constant = const0; 00070 } 00071 void set_statistics ( int dimy0, const ldmat V0, double nu0 = -1.0 ) { 00072 est.set_parameters ( dimy0, V0, nu0 ); 00073 est.validate(); 00074 last_lognc = est.lognc(); 00075 dimy = dimy0; 00076 } 00078 00080 void set_statistics ( const BMEF* BM0 ); 00081 00084 00086 void bayes_weighted ( const vec &yt, const vec &cond = empty_vec, const double w = 1.0 ); 00087 void bayes ( const vec &yt, const vec &cond = empty_vec ) { 00088 bayes_weighted ( yt, cond, 1.0 ); 00089 }; 00090 double logpred ( const vec &yt, const vec &cond ) const; 00091 vec samplepred ( const vec &cond ) ; 00092 void flatten ( const BMEF* B , double weight ); 00094 enorm<ldmat>* epredictor ( const vec &rgr ) const; 00095 estudent<ldmat>* epredictor_student ( const vec &rgr ) const; 00097 template<class sq_T> 00098 shared_ptr<mlnorm<sq_T> > ml_predictor() const; 00100 template<class sq_T> 00101 void ml_predictor_update ( mlnorm<sq_T> &pred ) const; 00102 mlstudent* predictor() const; 00104 ivec structure_est ( const egiw &Eg0 ); 00106 ivec structure_est_LT ( const egiw &Eg0 ); 00108 void reduce_structure(ivec &inds_in_V) { 00109 ldmat V = posterior()._V(); 00110 if (max(inds_in_V)>=V.rows()) { 00111 bdm_error("Incompatible structure"); 00112 } 00113 00114 ldmat newV(V,inds_in_V); 00115 est.set_parameters(dimy,newV, posterior()._nu()); 00116 00117 if (have_constant) { 00118 ivec rgr_elem= find(inds_in_V<(V.rows()-1)); // < -- find non-constant 00119 rgr = rgr.subselect(rgr_elem); 00120 rgrlen = rgr_elem.length(); 00121 } else { 00122 rgr = rgr.subselect(inds_in_V); 00123 } 00124 validate(); 00125 } 00127 00131 const egiw& posterior() const { 00132 return est; 00133 } 00135 00155 void from_setting ( const Setting &set ); 00156 00157 void validate() { 00158 BMEF::validate(); 00159 est.validate(); 00160 00161 // When statistics is defined, it has priority 00162 if(posterior()._dimx()>0) {//statistics is assigned 00163 dimy = posterior()._dimx(); 00164 rgrlen=posterior()._V().rows() - dimy - int ( have_constant == true ); 00165 dimc = rgrlen; 00166 } else { // statistics is not assigned - build it from dimy and rgrlen 00167 bdm_assert(dimy>0,"No way to validate egiw: empty statistics and empty dimy"); 00168 est.set_parameters(dimy, zeros(dimy+rgrlen+int(have_constant==true))); 00169 set_prior_default(est); 00170 } 00171 if (alter_est.dimension()==0) alter_est=est; 00172 00173 dyad = ones ( est._V().rows() ); 00174 } 00176 void set_prior ( const epdf *prior ); 00178 void set_prior_default ( egiw &prior ) { 00179 //assume 00180 vec dV0 ( prior._V().rows() ); 00181 dV0.set_subvector ( 0, prior._dimx() - 1, 1.0 ); 00182 if (dV0.length()>prior._dimx()) 00183 dV0.set_subvector ( prior._dimx(), dV0.length() - 1, 1e-5 ); 00184 00185 prior.set_parameters ( prior._dimx(), ldmat ( dV0 ) ); 00186 prior.validate(); 00187 } 00188 00189 void to_setting ( Setting &set ) const 00190 { 00191 BMEF::to_setting( set ); // takes care of rv, yrv, rvc 00192 UI::save(rgr, set, "rgr"); 00193 int constant = have_constant ? 1 : 0; 00194 UI::save(constant, set, "constant"); 00195 UI::save(&alter_est, set, "alternative"); 00196 UI::save(&posterior(), set, "posterior"); 00197 00198 } 00200 RV & _rgr() { 00201 return rgr; 00202 } 00203 bool _have_constant() { 00204 return have_constant; 00205 } 00206 int _rgrlen() { 00207 return rgrlen; 00208 } 00209 }; 00210 UIREGISTER ( ARX ); 00211 SHAREDPTR ( ARX ); 00212 00214 class ARXls : public BMEF{ 00215 public: 00216 egw_ls<ldmat> est; 00217 00218 egw_ls<ldmat> alternative; 00219 00220 const egw_ls<ldmat>& posterior() {return est;}; 00221 00222 void bayes(const vec &dt, const vec &psi){ 00223 // ldmat &Pbeta = est.P; 00224 // ldmat &Palpha = alternative.P; 00225 00226 00227 } 00228 }; 00229 00235 class ARXfrg : public ARX { 00236 public: 00237 ARXfrg() : ARX() {}; 00239 ARXfrg ( const ARXfrg &A0 ) : ARX ( A0 ) {}; 00240 virtual ARXfrg* _copy() const { 00241 ARXfrg *A = new ARXfrg ( *this ); 00242 return A; 00243 } 00244 00245 void bayes ( const vec &val, const vec &cond ) { 00246 bdm_assert_debug(cond.size()>rgrlen, "ARXfrg: Insufficient conditioning, frg not given."); 00247 frg = cond ( rgrlen); // the first part after rgrlen 00248 ARX::bayes ( val, cond.left(rgrlen) ); 00249 } 00250 void validate() { 00251 ARX::validate(); 00252 rvc.add ( RV ( "{phi }", vec_1 ( 1 ) ) ); 00253 dimc += 1; 00254 } 00255 }; 00256 UIREGISTER ( ARXfrg ); 00257 00263 class ARXmaxent : public ARX { 00264 protected: 00265 double maxent_frg; 00266 public: 00267 ARXmaxent() : ARX() {}; 00269 ARXmaxent ( const ARXmaxent &A0 ) : ARX ( A0 ),maxent_frg(A0.maxent_frg) {}; 00270 virtual ARXmaxent* _copy() const { 00271 ARXmaxent *A = new ARXmaxent ( *this ); 00272 return A; 00273 } 00274 00275 void bayes ( const vec &val, const vec &cond ) { 00276 ARX::bayes_weighted ( val, cond, maxent_frg ); 00277 } 00278 void from_setting(const Setting &set){ 00279 ARX::from_setting(set); 00280 maxent_frg=frg; 00281 frg = 1.0; 00282 } 00283 }; 00284 00285 UIREGISTER ( ARXmaxent ); 00286 00291 class ARXwin : public ARX { 00292 protected: 00293 mat Y; 00294 mat Cond; 00295 00296 int win_length; 00297 ldmat V0; 00298 double nu0; 00299 public: 00300 ARXwin() : ARX() {}; 00302 void set_parameters(const int win_length0){ 00303 win_length=win_length0; 00304 } 00305 ARXwin ( const ARXwin &A0 ) : ARX(A0), Y( A0.Y), Cond(A0.Cond) {}; 00306 00307 virtual ARXwin* _copy() const { 00308 ARXwin *A = new ARXwin ( *this ); 00309 return A; 00310 } 00311 00312 void bayes ( const vec &val, const vec &cond ) { 00313 00314 if(cond.size()>0) 00315 { 00316 // fill window 00317 Y.append_col(val); 00318 Cond.append_col(cond); 00319 if (Y.cols()>win_length){ 00320 // shift the buffer 00321 Y=Y.get_cols(1,Y.cols()-1); 00322 Cond=Cond.get_cols(1,Cond.cols()-1); 00323 } 00324 00325 est._V()=V0; 00326 est._nu()=nu0; 00327 for ( int t = 0; t < Y.cols(); t++ ) { 00328 ARX::bayes ( Y.get_col ( t ), Cond.get_col ( t ) ); 00329 } 00330 } 00331 else 00332 { 00333 Y.append_col(val); 00334 00335 if (Y.cols()>win_length){ 00336 // shift the buffer 00337 Y=Y.get_cols(1,Y.cols()-1); 00338 } 00339 00340 est._V()=V0; 00341 est._nu()=nu0; 00342 00343 for ( int t = 0; t < Y.cols(); t++ ) { 00344 ARX::bayes (Y.get_col ( t )); 00345 } 00346 } 00347 00348 } 00349 00350 void from_setting(const Setting &set){ 00351 ARX::from_setting(set); 00352 UI::get(win_length,set,"win_length",UI::compulsory); 00353 } 00354 void validate() { 00355 ARX::validate(); 00356 V0=est._V(); 00357 nu0=est._nu(); 00358 } 00359 }; 00360 00361 UIREGISTER ( ARXwin ); 00362 00363 00364 00385 class ARXpartialforg : public ARX { 00386 public: 00387 ARXpartialforg() : ARX(1.0) {}; 00389 ARXpartialforg ( const ARXpartialforg &A0 ) : ARX ( A0 ) {}; 00390 virtual ARXpartialforg* _copy() const { 00391 ARXpartialforg *A = new ARXpartialforg ( *this ); 00392 return A; 00393 } 00394 00395 void bayes ( const vec &val, const vec &cond ); 00396 00397 void validate() { 00398 ARX::validate(); 00399 int philen = 1 << (est._V().cols() - 1); 00400 rvc.add ( RV ( "{phi }", vec_1(philen) ) ); // pocet 2^parametru 00401 dimc += philen; 00402 } 00403 }; 00404 UIREGISTER ( ARXpartialforg ); 00405 00406 00408 template<class sq_T> 00409 shared_ptr< mlnorm<sq_T> > ARX::ml_predictor ( ) const { 00410 shared_ptr< mlnorm<sq_T> > tmp = new mlnorm<sq_T> ( ); 00411 tmp->set_rv ( yrv ); 00412 tmp->set_rvc ( _rvc() ); 00413 00414 ml_predictor_update ( *tmp ); 00415 tmp->validate(); 00416 return tmp; 00417 } 00418 00419 template<class sq_T> 00420 void ARX::ml_predictor_update ( mlnorm<sq_T> &pred ) const { 00421 mat mu ( dimy, posterior()._V().rows() - dimy ); 00422 mat R ( dimy, dimy ); 00423 00424 est.mean_mat ( mu, R ); //mu = 00425 mu = mu.T(); 00426 //correction for student-t -- TODO check if correct!! 00427 00428 if ( have_constant ) { // constant term 00429 //Assume the constant term is the last one: 00430 pred.set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), sq_T ( R ) ); 00431 } else { 00432 pred.set_parameters ( mu, zeros ( dimy ), sq_T ( R ) ); 00433 } 00434 } 00435 00436 }; 00437 #endif // AR_H 00438 00439
Generated on 2 Dec 2013 for mixpp by
1.4.7
