Changeset 625

Timestamp:

09/18/09 00:17:05 (15 years ago)

Author:

smidl

Message:

ARX re-designed

Location:

library/bdm/estim

Files:

• arx.cpp (modified) (9 diffs)
• arx.h (modified) (5 diffs)
• kalman.h (modified) (1 diff)

library/bdm/estim/arx.cpp

@@ -5 +5 @@
         double lnc;
 
+                
         if ( frg < 1.0 ) {
                 est.pow ( frg );
@@ -11 +12 @@
                 }
         }
-        V.opupdt ( dt, w );
+        if (have_constant) {
+                _dt.set_subvector(0,dt);
+                V.opupdt ( _dt, w );
+        } else {
+                V.opupdt ( dt, w );
+        }
         nu += w;
 
@@ -78 +84 @@
 mlnorm<ldmat>* ARX::predictor ( ) const {
         int dim = est.dimension();
-        int dif = V.rows() - dim ;///<----------- TODO
-        bdm_assert_debug ( ( dif == 0 ) || ( dif == 1 ), "Give RVs do not match" );
-
+        
         mat mu ( dim, V.rows() - dim );
         mat R ( dim, dim );
@@ -90 +94 @@
         //correction for student-t  -- TODO check if correct!!
 
-        if ( dif == 0 ) { // no constant term
-                tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );
-        } else {
+        if ( have_constant) { // constant term
                 //Assume the constant term is the last one:
                 tmp->set_parameters ( mu.get_cols ( 0, mu.cols() - 2 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ) );
+        } else {
+                tmp->set_parameters ( mu, zeros ( dim ), ldmat ( R ) );
         }
         return tmp;
@@ -101 +105 @@
 mlstudent* ARX::predictor_student ( ) const {
         int dim = est.dimension();
-        int dif = V.rows() - est.dimension();//-------------TODO
-        bdm_assert_debug ( ( dif == 0 ) || ( dif == 1 ), "Give RVs do not match" );
 
         mat mu ( dim, V.rows() - dim );
@@ -117 +119 @@
 
 
-        if ( dif == 0 ) { // no constant term
-                tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam );
-        } else {
+        if ( have_constant) { // no constant term
                 //Assume the constant term is the last one:
                 if ( mu.cols() > 1 ) {
@@ -126 +126 @@
                         tmp->set_parameters ( mat ( dim, 0 ), mu.get_col ( mu.cols() - 1 ), ldmat ( R ), Lam );
                 }
+        } else {
+                // no constant term
+                tmp->set_parameters ( mu, zeros ( xdim ), ldmat ( R ), Lam );
         }
         return tmp;
@@ -203 +206 @@
 
 void ARX::from_setting ( const Setting &set ) {
-        shared_ptr<RV> yrv = UI::build<RV> ( set, "y", UI::compulsory );
+        shared_ptr<RV> yrv = UI::build<RV> ( set, "rv", UI::compulsory );
         shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );
         int ylen = yrv->_dsize();
+        // rgrlen - including constant!!!
         int rgrlen = rrv->_dsize();
-
+        
+        set_rv ( *yrv, *rrv );
+        
         string opt;
         if ( UI::get(opt, set,  "options", UI::optional) ) {
                 BM::set_options(opt);
         }
+        if (!UI::get(have_constant, set, "constant", UI::optional)){
+                have_constant=true;
+        }
+        if (have_constant) {rgrlen++;_dt=ones(rgrlen+ylen);}
 
         //init
         mat V0;
         vec dV0;
-        if ( !UI::get ( dV0, set, "dV0" ) )
-                dV0 = concat ( 1e-3 * ones ( ylen ), 1e-5 * ones ( rgrlen ) );
-        V0 = diag ( dV0 );
-
+        if (!UI::get(V0, set, "V0",UI::optional)){
+                if ( !UI::get ( dV0, set, "dV0" ) )
+                        dV0 = concat ( 1e-3 * ones ( ylen ), 1e-5 * ones ( rgrlen ) );
+                V0 = diag ( dV0 );
+        }
         double nu0;
         if ( !UI::get ( nu0, set, "nu0" ) )
@@ -230 +241 @@
         set_parameters ( frg );
         set_statistics ( ylen, V0, nu0 );
-        set_drv ( concat ( *yrv, *rrv ) );
-
+        
         //name results (for logging)
-        set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) );
-
-}
-
-}
+        shared_ptr<RV> rv_par=UI::build<RV>(set, "rv_param",UI::optional );
+        if (!rv_par){
+                est.set_rv ( RV ( "{theta r }", vec_2 ( ylen*rgrlen, ylen*ylen ) ) );
+        } else {
+                est.set_rv ( *rv_par );
+        }
+        validate();
+}
+
+}

library/bdm/estim/arx.h

@@ -48 +48 @@
         //! Do NOT access directly, only via \c get_yrv().
         RV _yrv;
+        //! rv of regressor
+        RV rgrrv;
         //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density
         egiw est;
@@ -54 +56 @@
         //! cached value of est.nu
         double ν
+        //! switch if constant is modelled or not
+        bool have_constant;
+        //! cached value of data vector for have_constant =true
+        vec _dt;
 public:
         //! \name Constructors
@@ -73 +79 @@
         //!@}
 
-//      //! Set parameters given by moments, \c mu (mean of theta), \c R (mean of R) and \c C (variance of theta)
-//      void set_parameters ( const vec &mu, const mat &R, const mat &C, double dfm){};
         //! Set sufficient statistics
         void set_statistics ( const BMEF* BM0 );
-//      //! Returns sufficient statistics
-//      void get_parameters ( mat &V0, double &nu0 ) {V0=est._V().to_mat(); nu0=est._nu();}
+
         //!\name Mathematical operations
         //!@{
@@ -121 +124 @@
         //!\name Connection
         //!@{
-        void set_drv ( const RV &drv0 ) {
-                drv = drv0;
+        void set_rv ( const RV &yrv0 , const RV &rgrrv0 ) {
+                _yrv = yrv0;
+                rgrrv=rgrrv0;
+                set_drv(concat(yrv0, rgrrv));
         }
 
@@ -140 +145 @@
         //!@}
 
-        // TODO dokumentace - aktualizovat
         /*! UI for ARX estimator
 
-        The ARX is constructed from a structure with fields:
         \code
-        estimator = {
-                class = "ARX";
-                y = {type="rv", ...}   // description of output variables
-                rgr = {type="rv", ...} // description of regressor variables
-                constant = true;       // boolean switch if the constant term is modelled or not
+        class = 'ARX';
+        rv    = RV({names_of_dt} )                 // description of output variables
+        rgr   = RV({names_of_regressors}, [-1,-2]} // description of regressor variables
+        constant = true;                           // boolean switch if the constant term is modelled or not
 
-                //optional fields
-                dV0 = [1e-3, 1e-5, 1e-5, 1e-5];
-                                                           // default: 1e-3 for y, 1e-5 for rgr
-                nu0 = 6;               // default: rgrlen + 2
-                frg = 1.0;             // forgetting, default frg=1.0
-        };
+        --- optional ---
+        V0  = [1 0;0 1];                           // Initial value of information matrix V
+          --- OR ---
+        dV0 = [1e-3, 1e-5, 1e-5, 1e-5];            // Initial value of diagonal of information matrix V
+                                                                                           // default: 1e-3 for rv, 1e-5 for rgr
+        nu0 = 6;                                                   // initial value of nu, default: rgrlen + 2
+        frg = 1.0;                                 // forgetting, default frg=1.0
+
+        rv_param   = RV({names_of_parameters}}     // description of parametetr names 
+                                                                                           // default: ["theta_i" and "r_i"]
         \endcode
-
-        The estimator will assign names of the posterior in the form ["theta_i" and "r_i"]
         */
         void from_setting ( const Setting &set );
 
+        void validate() {
+                bdm_assert(dimx == _yrv._dsize(), "RVs of parameters and regressor do not match");
+                
+        }
 };
 

library/bdm/estim/kalman.h

@@ -378 +378 @@
         public:
                 //! set up this object to match given mlnorm
-                void connect_mlnorm(const mlnorm<fsqmat > &ml){
-        //get ids of yrv                                
+                void connect_mlnorm(const mlnorm<fsqmat> &ml){
+                        //get ids of yrv                                
                         const RV &yrv = ml._rv();
                         //need to determine u_t - it is all in _rvc that is not in ml._rv()