Changeset 996 for library

Timestamp:

05/26/10 23:33:28 (14 years ago)

Author:

smidl

Message:

Scheduling of forgetting factor

Location:

library/bdm

Files:

• estim/arx.cpp (modified) (2 diffs)
• estim/arx.h (modified) (3 diffs)
• estim/arx_straux.h (modified) (1 diff)
• itpp_ext.cpp (modified) (1 diff)
• stat/exp_family.cpp (modified) (1 diff)
• stat/exp_family.h (modified) (3 diffs)

library/bdm/estim/arx.cpp

@@ -3 +3 @@
 
 void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
-
         bdm_assert_debug ( yt.length() == dimy, "ARX::bayes yt is of size "+num2str(yt.length())+" expected dimension is "+num2str(dimy) );
         bdm_assert_debug ( cond.length() == rgrlen , "ARX::bayes cond is of size "+num2str(cond.length())+" expected dimension is "+num2str(rgrlen) );
+        
+        BMEF::bayes_weighted(yt,cond,w); //potential discount scheduling
+        
         double lnc;
         //cache
@@ -216 +218 @@
 ivec ARX::structure_est_LT ( egiw est0 ) {
         //some stuff with beliefs etc.
-//      ivec ind = bdm::straux1(V,nu, est0._V(), est0._nu());
-        return ivec();//ind;
+        ivec belief = vec_1 ( 2 );        // default belief
+        int nbest = 1;           // nbest: how many regressors are returned
+        int nrep = 5;         // nrep: number of random repetions of structure estimation
+        double lambda   = 0.9;
+        int k = 2;
+        
+        Array<str_aux> o2;
+        
+        ivec ind = bdm::straux1(est._V(),est._nu(), est0._V(), est0._nu(), belief, nbest, nrep, lambda, k, o2);
+        
+        return ind;
 }
 

library/bdm/estim/arx.h

@@ -105 +105 @@
         //! Smarter structure estimation by Ludvik Tesar.\return indices of accepted regressors.
         ivec structure_est_LT ( egiw Eg0 );
+        //! reduce to
+        void reduce_structure(ivec &inds_in_V){
+                ldmat V = posterior()._V();
+                if (max(inds_in_V)>=V.rows()) {bdm_error("Incompatible structure");}
+                
+                ldmat newV(V,inds_in_V);
+                est.set_parameters(dimy,newV, posterior()._nu());
+                validate();
+        }
         //!@}
 
@@ -138 +147 @@
                 BMEF::validate();       
                 est.validate();
-                //if dimc not set set it from V
-                if(dimy>0) {//statistics is assigned
-                        if (posterior()._V().rows()>dimy) {//statistics is assigned
-                                rgrlen=posterior()._V().rows() - dimy - int ( have_constant == true );
-                        }
-                } else{
-                        bdm_error("No posterior or yrv assigned matrix assigned");
-                }
-                dimc =rgrlen;
-                
-                if(est._dimx()==0) { // no prior 
+                
+                // When statistics is defined, it has priority
+                if(posterior()._dimx()>0) {//statistics is assigned
+                        dimy = posterior()._dimx()>0;
+                        rgrlen=posterior()._V().rows() - dimy - int ( have_constant == true );
+                        dimc = rgrlen;
+                } else{  // statistics is not assigned - build it from dimy and rgrlen
+                        bdm_assert(dimy>0,"No way to validate egiw: empty statistics and empty dimy");
                         est.set_parameters(dimy, zeros(dimy+rgrlen+int(have_constant==true)));
                         set_prior_default(est);
@@ -213 +219 @@
 
 
+
 ////////////////////
 template<class sq_T>

library/bdm/estim/arx_straux.h

@@ -48 +48 @@
 
 
-//! Rplication of Ludvik Tesar original straux1 from mixtools straux1
+//! Replication of Ludvik Tesar original straux1 from mixtools straux1
 ivec straux1 ( ldmat Ld, double nu, ldmat Ld0, double nu0, ivec belief, int nbest, int max_nrep, double lambda, int order_k, Array<str_aux> &rgrsout );
 

library/bdm/itpp_ext.cpp

@@ -55 +55 @@
 }
 
+//! return vector of elements at positions given by indexlist
 vec get_vec ( const vec &v, const ivec &indexlist ) {
         int size = indexlist.size();

library/bdm/stat/exp_family.cpp

@@ -404 +404 @@
 }
 
-// mat egamma::sample ( int N ) const {
-//      mat Smp ( rv.count(),N );
-//      int i,j;
-//
-//      for ( i=0; i<rv.count(); i++ ) {
-//              GamRNG.setup ( alpha ( i ),beta ( i ) );
-//
-//              for ( j=0; j<N; j++ ) {
-//                      Smp ( i,j ) = GamRNG();
-//              }
-//      }
-//
-//      return Smp;
-// }
 
 double egamma::evallog ( const vec &val ) const {

library/bdm/stat/exp_family.h

@@ -81 +81 @@
         //! cached value of lognc() in the previous step (used in evaluation of \c ll )
         double last_lognc;
+        //! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
+        double frg_sched_factor;
 public:
         //! Default constructor (=empty constructor)
-        BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ) {}
+        BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
         //! Copy constructor
-        BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
+        BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
         //!get statistics from another model
         virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;
 
         //! Weighted update of sufficient statistics (Bayes rule)
-        virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {};
+        virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
+                if (frg_sched_factor>0) {frg = frg*(1-frg_sched_factor)+frg_sched_factor;}
+        };
         //original Bayes
         void bayes ( const vec &yt, const vec &cond = empty_vec );
@@ -102 +106 @@
                 BM::to_setting( set );
                 UI::save(frg, set, "frg");
+                UI::save( frg_sched_factor, set, "frg_sched_factor" );
         } 
 
@@ -108 +113 @@
                 if ( !UI::get ( frg, set, "frg" ) )
                         frg = 1.0;
+                UI::get ( frg_sched_factor, set, "frg_sched_factor",UI::optional );
         }