mixpp: arx.h Source File

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  doxygen 1.4.7