Changeset 265
- Timestamp:
- 02/09/09 23:14:58 (16 years ago)
- Files:
-
- 4 added
- 6 modified
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
bdm/stat/libBM.cpp
r263 r265 52 52 RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ), sizes ( 0 ), times ( 0 ),names ( 0 ) {}; 53 53 54 RV::RV(string name, int id){ 55 it_assert(id>=RVcounter,"RV::This id is already taken"); 56 RVcounter=id; 54 RV::RV(string name, int id, int sz, int tm){ 55 if (id>RVcounter) {RVcounter=id;}; 57 56 Array<string> A(1); A(0)=name; 58 init(vec_1(id),A,vec_1( 1),vec_1(0));57 init(vec_1(id),A,vec_1(sz),vec_1(tm)); 59 58 } 60 59 -
bdm/stat/libBM.h
r263 r265 78 78 RV (); 79 79 //! Constructor of a single RV with given id 80 RV (string name, int id );80 RV (string name, int id, int sz=1, int tm=0); 81 81 82 82 //! Printing output e.g. for debugging. … … 490 490 L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) ); 491 491 } 492 }; 493 494 /*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities. 492 //!access function 493 virtual RV _drv() const {return concat(Drv,Urv);} 494 //!access function 495 const RV& _urv() const {return Urv;} 496 }; 497 498 /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities. 495 499 496 500 */ … … 500 504 //!Random variable of the posterior 501 505 RV rv; 506 //! Random variable of the data (optional) 507 RV drv; 502 508 //!Logarithm of marginalized data likelihood. 503 509 double ll; … … 537 543 //!access function 538 544 const RV& _rv() const {return rv;} 545 //!access function 546 const RV& _drv() const {return drv;} 547 //!set drv 548 void set_drv(const RV &rv){drv=rv;} 539 549 //!access function 540 550 double _ll() const {return ll;} -
bdm/stat/libDS.cpp
r263 r265 44 44 void ArxDS::step() { 45 45 double tmp; 46 //shift history 47 H.shift_right ( 0, Drv.count() +Urv.count() ); 48 49 H.set_subvector ( Drv.count(), U ); // write U after Drv 50 46 51 //get regressor 47 52 rgr = rgrlnk.get_val ( H ); 48 53 // Eval Y 49 Y = model.samplecond ( rgr,tmp);54 H.set_subvector ( 0, model.samplecond ( rgr,tmp ) ); 50 55 51 //shift history52 H.shift_right ( 0, Drv.count() +Urv.count() );53 //fill new data54 H.set_subvector ( Drv.count(),Y );55 56 //Leaving U.length() empty slots - these will be filled by write()57 56 } 58 57 … … 60 59 RV fullrgr ( const RV &drv0, const RV &urv0, const RV &rrv0 ) { 61 60 RV T ( urv0 ); 62 RV pom = concat (drv0, urv0);61 RV pom = concat ( drv0, urv0 ); 63 62 int mint = rrv0.mint(); 64 63 for ( int i=0; i>mint; i-- ) { … … 69 68 } 70 69 71 ArxDS::ArxDS ( RV &drv, RV &urv, RV &rrv ) : 72 DS (drv,urv), Rrv(rrv), Hrv( fullrgr ( drv,urv,rrv )), //RVs73 Y(drv.count()), H(Hrv.count()),rgr ( rrv.count() ), //tmp variables74 rgrlnk (Hrv ,rrv ) ,model ( drv,rrv ) {70 ArxDS::ArxDS ( RV &drv, RV &urv, RV &rrv ) : 71 DS ( drv,urv ), Rrv ( rrv ), Hrv ( concat(drv,fullrgr ( drv,urv,rrv )) ), //RVs 72 H ( Hrv.count() ) ,U ( urv.count() ),rgr ( rrv.count() ), //tmp variables 73 rgrlnk ( rrv, Hrv ) ,model ( drv,rrv ) { 75 74 } -
bdm/stat/libDS.h
r263 r265 53 53 */ 54 54 class ArxDS : public DS { 55 protected: 55 56 //! Rv of the regressor 56 57 RV Rrv; 57 58 //! Rv of the history (full regressor) 58 59 RV Hrv; 59 //! Internal storage of results 60 vec Y; 61 //! History, ordered as \f$[u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 60 //! History, ordered as \f$[y_t, u_t, y_{t-1 }, u_{t-1}, \ldots]\f$ 62 61 vec H; 62 //! (future) input 63 vec U; 63 64 //! temporary variable for regressor 64 65 vec rgr; 65 //! data link: val = y, cond = u; local = rgr;66 //! data link: H -> rgr 66 67 datalink_e2e rgrlnk; 67 //! model of Y - linear Gaussian 68 //! model of Y - linear Gaussian 68 69 mlnorm<chmat> model; 70 //! options 71 bool opt_L_theta; 72 //! loggers 73 int L_theta; 74 int L_R; 69 75 public: 70 void getdata ( vec &dt ){it_assert_debug(dt.length()==Y.length(),"ArxDS"); 71 dt=concat(Y,H.left(Urv.count()));}; 72 void getdata ( vec &dt, const ivec &indexes ){it_assert_debug(dt.length()==Y.length(),"ArxDS"); dt=Y(indexes);}; 73 void write ( vec &ut ){it_assert_debug(ut.length()==Urv.count(),"ArxDS"); H.set_subvector(0,ut);}; 74 void write ( vec &ut, const ivec &indexes ){it_assert_debug(ut.length()==Urv.count(),"ArxDS"); set_subvector(H,indexes,ut);}; 76 void getdata ( vec &dt ) { 77 it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" ); 78 dt=H.left ( Urv.count() +Drv.count() ); 79 }; 80 void getdata ( vec &dt, const ivec &indexes ) { 81 it_assert_debug ( dt.length() ==indeces.length(),"ArxDS" ); 82 dt=H ( indexes ); 83 }; 84 void write ( vec &ut ) { 85 it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" ); 86 U=ut; 87 }; 88 void write ( vec &ut, const ivec &indexes ) { 89 it_assert_debug ( ut.length() ==indeces.length(),"ArxDS" ); 90 set_subvector ( U, indexes,ut ); 91 }; 75 92 void step(); 76 93 //!Default constructor 77 ArxDS ( RV &drv, RV &urv, RV &rrv );94 ArxDS ( RV &drv, RV &urv, RV &rrv ); 78 95 //! Set parameters of the internal model 79 void set_parameters(const mat &Th0, const vec mu0, const chmat &sqR0) 80 {model.set_parameters(Th0, mu0, sqR0); }; 96 void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 ) 97 { model.set_parameters ( Th0, mu0, sqR0 ); }; 98 //! set options from a string 99 void set_options ( const string &s ) { 100 opt_L_theta= ( s.find ( "L_theta" ) !=string::npos ); 101 }; 102 virtual void log_add ( logger &L ) { 103 DS::log_add ( L ); 104 mat &A =model._A(); 105 mat R =model._R(); 106 if ( opt_L_theta ) {L_theta=L.add ( RV("{theta }", vec_1(A.rows() *A.cols())),"t" );} 107 if ( opt_L_theta ) {L_R=L.add ( RV("{R }", vec_1(R.rows() *R.cols())),"r" );} 108 } 109 virtual void logit ( logger &L ) { 110 DS::logit ( L ); 111 mat &A =model._A(); 112 mat R =model._R(); 113 if ( opt_L_theta ) {L.logit ( L_theta,vec ( A._data(), A.rows() *A.cols() ) );}; 114 if ( opt_L_theta ) {L.logit ( L_R, vec ( R._data(), R.rows() *R.rows() ) );}; 115 } 116 81 117 }; 82 118 119 class ARXDS : public ArxDS { 120 public: 121 ARXDS ( RV &drv, RV &urv, RV &rrv ) : ArxDS ( drv,urv,rrv ) {} 122 123 void getdata ( vec &dt ) {dt=H;} 124 void getdata ( vec &dt, const ivec &indeces ) {dt=H ( indeces );} 125 virtual RV _drv() const {return Hrv;} 126 127 }; 83 128 }; //namespace 84 129 -
matlab/mex/arx_example.cfg
r263 r265 18 18 r = [0.1, 0.0, 19 19 0.0, 1.0]; 20 opt="L_theta"; 20 21 }; 21 22 22 23 //store results 23 24 logger = { 24 type= " dirfilelog";25 dirname = "exp/arx_ui";26 maxlen = 1000; // 25 type= "mexlog"; 26 maxlen = 90; // 27 dirname = "exp/ax"; 27 28 }; 28 29 29 30 //estimation 30 31 estimator = { 31 type = "ARX"; 32 32 type = "ARXest"; 33 ychns = [1]; 34 rgrid = [ 1, 1, 1, 2]; 35 delays = [-1, -2, -3, -1]; 36 37 //optional fields 38 dV0 = [1.0, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr 39 nu0 = 8.; //default: rgrlen + 2 40 frg = 1.; // forgetting, default frg=1.0 33 41 }; 34 42 35 43 //experiment description 36 44 experiment:{ 37 ndat = 90 00;45 ndat = 90; 38 46 }; -
tests/UI/UIArxDS_test.cfg
r263 r265 18 18 r = [0.1, 0.0, 19 19 0.0, 1.0]; 20 opt="L_theta"; 20 21 }; 21 22 … … 29 30 //estimation 30 31 estimator = { 31 type = "ARX"; 32 32 type = "ARXest"; 33 ychns = [1]; 34 rgrid = [ 1, 1, 1, 2]; 35 delays = [-1, -2, -3, -1]; 36 37 //optional fields 38 dV0 = [1e-3, 1e-5, 1e-5, 1e-5, 1e-5]; //default: 1e-3 for y, 1e-5 for rgr 39 nu0 = 8.; //default: rgrlen + 2 40 frg = .990; // forgetting, default frg=1.0 33 41 }; 34 42 -
tests/UI/UIArxDS_test.cpp
r263 r265 1 1 /*! 2 \file 2 \file 3 3 \brief Test of UI builders for ARX 4 4 … … 9 9 10 10 #include <stat/libDS_ui.h> 11 #include <estim/arx_ui.h> 11 12 12 13 using namespace bdm; 13 14 int main() { 14 UIFile F ("UIArxDS_test.cfg");15 UIFile F ( "UIArxDS_test.cfg" ); 15 16 16 logger* L; 17 UIbuild (F.lookup("logger"),L);17 logger* L; 18 UIbuild ( F.lookup ( "logger" ),L ); 18 19 ArxDS * DS; 19 UIbuild (F.lookup("system"),DS);20 UIbuild ( F.lookup ( "system" ),DS ); 20 21 int Ndat; 21 F.lookupValue ("experiment.ndat",Ndat);22 F.lookupValue ( "experiment.ndat",Ndat ); 22 23 // SET SIMULATOR 24 BM* E; 25 UIbuild ( F.lookup ( "estimator" ),E ); 23 26 24 DS->log_add(*L); 27 DS->log_add ( *L ); 28 int L_est= L->add ( E->_epdf()._rv(), "est" ); // estimate 29 int L_lb = L->add ( E->_epdf()._rv(), "lb" ); // lower bound 30 int L_ub = L->add ( E->_epdf()._rv(), "ub" ); // upper bound 25 31 L->init(); 26 32 33 vec dt=zeros ( DS->_drv().count() ); //data variable 34 datalink_e2e dl ( E->_drv(),DS->_drv() ); //datalink between a datasource and estimator 35 27 36 for ( int tK=1;tK<Ndat;tK++ ) { 28 DS->step(); 29 DS->logit(*L); 37 DS->step(); // simulator step 38 DS->getdata ( dt ); // read data 39 E->bayes ( dl.get_val ( dt ) ); // update estimates 40 41 DS->logit ( *L ); 42 L->logit ( L_est, E->_epdf().mean() ); 43 L->logit ( L_lb, E->_epdf().mean()-2*sqrt ( E->_epdf().variance() ) ); 44 L->logit ( L_ub, E->_epdf().mean() +2*sqrt ( E->_epdf().variance() ) ); 45 30 46 L->step(); 31 47 } 32 33 L->finalize(); 34 48 49 L->finalize(); 50 35 51 delete L; 36 52 delete DS; 53 delete E; 37 54 return 0; 38 55 }