Changeset 964

Timestamp:

05/21/10 00:44:04 (14 years ago)

Author:

smidl

Message:

Corrections in ARX and PF

Files:

• applications/bdmtoolbox/tutorial/userguide/arx_basic_example.m (modified) (1 diff)
• applications/bdmtoolbox/tutorial/userguide/frg_estim.m (modified) (3 diffs)
• library/bdm/estim/arx.cpp (modified) (4 diffs)
• library/bdm/estim/arx.h (modified) (6 diffs)
• library/bdm/estim/kalman.h (modified) (1 diff)
• library/bdm/estim/particles.h (modified) (3 diffs)
• library/bdm/stat/emix.cpp (modified) (1 diff)
• library/bdm/stat/exp_family.h (modified) (1 diff)

applications/bdmtoolbox/tutorial/userguide/arx_basic_example.m

@@ -10 +10 @@
 
 A1.class = 'ARX';
-A1.rv = y;
+A1.yrv = y;
 A1.rgr = RVtimes([y,u],[-3,-1]) ; % correct structure is {y,y}
 A1.log_level = 'logbounds,logevidence';

applications/bdmtoolbox/tutorial/userguide/frg_estim.m

@@ -10 +10 @@
 
 A1.class = 'ARXfrg';
-A1.rv = y;
+A1.yrv = y;
+A1.rv = RV({'theta','r'},[3,1]);
 A1.rgr = RVtimes([y,u],[-3,-1]) ; 
 A1.log_level = 'logbounds';
@@ -20 +21 @@
 phi_pdf.class = 'mDirich';         % random walk on coefficient phi
 phi_pdf.rv    = RV({'phi','1_phi'});       % 2D random walk - frg is the first element
-phi_pdf.k     = 0.01;              % width of the random walk
+phi_pdf.k     = 0.001;              % width of the random walk
 phi_pdf.betac = [0.1 0.1];         % stabilizing elememnt of random walk
 
@@ -32 +33 @@
 E.class = 'PF';
 E.particle = p;                    % ARX is the analytical part
-E.res_threshold = 0.90;             % resampling parameter
-E.n = 100;                          % number of particles
+E.res_threshold = 0.0;             % resampling parameter
+E.n = 10;                          % number of particles
 E.prior.class = 'eDirich';         % prior on non-linear part
 E.prior.beta  = [2 1]; % 
-E.log_level = 'logbounds,logweights,logmeans';
+E.log_level = 'logbounds,logweights,logmeans,logvars';
 E.name = 'MPF';
 
-[M,Str]=estimator(DS,{E});
+A2=A1;
+A2.class='ARX';
+A2.frg=1.0;
+A2.name = 'MPFx';
+
+[M,Str]=estimator(DS,{E,A2});
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

library/bdm/estim/arx.cpp

@@ -4 +4 @@
 void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
 
-        bdm_assert_debug ( yt.length() >= dimy, "ARX::bayes yt is smaller then dimc" );
-        bdm_assert_debug ( cond.length() >= dimc, "ARX::bayes cond is smaller then dimc" );
+        bdm_assert_debug ( yt.length() == dimy, "ARX::bayes yt is smaller then dimensiony()" );
+        bdm_assert_debug ( (cond.length() == rgrlen - int(have_constant==true)) , "ARX::bayes dimension of cond is not rgrlen" );
         double lnc;
         //cache
@@ -217 +217 @@
         BMEF::from_setting(set);
         
-        shared_ptr<RV> yrv_ = UI::build<RV> ( set, "rv", UI::compulsory );
-        shared_ptr<RV> rrv = UI::build<RV> ( set, "rgr", UI::compulsory );
-        dimy = yrv_->_dsize();
+        RV rrv; 
+        UI::get (rrv, set, "rgr", UI::compulsory );
+        
+        dimy = yrv._dsize();
+        bdm_assert(dimy>0,"ARX::yrv should not be empty");
         // rgrlen - including constant!!!
-        dimc = rrv->_dsize();
-
-        yrv = *yrv_;
-        rvc = *rrv;
+        rgrlen = rrv._dsize();
 
         int constant;
@@ -232 +231 @@
                 have_constant = constant > 0;
         }
-        int rgrlen = dimc + int ( have_constant == true );
+        dimc = rgrlen;
+        rvc = rrv;
+        rgrlen += int ( have_constant == true );
 
         //init
@@ -259 +260 @@
         // frg handled by BMEF
 
-        //name results (for logging)
-        shared_ptr<RV> rv_par = UI::build<RV> ( set, "rv_param", UI::optional );
-        if ( !rv_par ) {
-                est.set_rv ( RV ( "{theta r }", vec_2 ( dimy*rgrlen, dimy*dimy ) ) );
-        } else {
-                est.set_rv ( *rv_par );
-        }       
-}
-
-}
-
+}
+
+}
+

library/bdm/estim/arx.h

@@ -47 +47 @@
         //! vector of dyadic update
         vec dyad;
+        //! length of the regressor
+        int rgrlen;
         //! posterior density
         egiw est;
@@ -54 +56 @@
         //! \name Constructors
         //!@{
-        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant ( true ), dyad(), est(), alter_est() {};
-        ARX ( const ARX &A0 ) : BMEF ( A0 ),  have_constant ( A0.have_constant ), dyad ( A0.dyad ), est ( A0.est ), alter_est ( A0.alter_est ) { };
+        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ),  have_constant ( true ), dyad(), rgrlen(),est(), alter_est() {};
+        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 ) { };
 
         ARX* _copy() const;
@@ -115 +117 @@
         \code
         class = 'ARX';
-        rv    = RV({names_of_dt} )                 // description of output variables
+        yrv   = RV({names_of_dt} )                 // description of output variables
         rgr   = RV({names_of_regressors}, [-1,-2]} // description of regressor variables
         constant = 1;                              // 0/1 switch if the constant term is modelled or not
@@ -125 +127 @@
         frg = 1.0;                                 // forgetting, default frg=1.0
 
-        rv_param   = RV({names_of_parameters}}     // description of parametetr names
-                                                                                           // default: ["theta_i" and "r_i"]
+        rv  = RV({names_of_parameters}}            // description of parametetr names
+                                                                                           // default: [""]
         \endcode
         */
@@ -133 +135 @@
         void validate() {
                 BMEF::validate();       
-
+                est.validate();
                 //if dimc not set set it from V
                 if ( dimc == 0 ) {
@@ -192 +194 @@
 
         void bayes ( const vec &val, const vec &cond ) {
-                frg = cond ( dimc - 1 ); //  last in cond is phi
-                ARX::bayes ( val, cond );
+                int arx_cond_size=rgrlen -int(have_constant==true);
+                bdm_assert_debug(cond.size()>arx_cond_size, "ARXfrg: Insufficient conditioning, frg not given.");
+                frg = cond ( arx_cond_size ); // the first part after rgrlen
+                ARX::bayes ( val, cond.left(arx_cond_size) );
         }
         void validate() {

library/bdm/estim/kalman.h

@@ -147 +147 @@
                 R     = [];                   // Matrix R
                 prior = struct('class','epdf_offspring');    // Prior density - will be converted to gaussian
-                rvy   = RV('some_names');     // Description of required observations
+                yrv   = RV('some_names');     // Description of required observations
                 rvc   = RV('some_names');     // Description of required inputs
                 \endcode

library/bdm/estim/particles.h

@@ -323 +323 @@
                 //set drv
 
-                set_yrv ( bm0->_rv() );
                 rvc = bm0->_rvc();
+                dimc = bm0->dimensionc();
                 BM::set_rv(bm0->_rv());
                 yrv=bm0->_yrv();
+                dimy = bm0->dimensiony();
         }
         
@@ -339 +340 @@
                         }
                 }
+                if (log_level[logvars]){
+                        for (int i=0; i<particles.length(); i++){
+                                L.add_vector( log_level, logvars, RV ( particles(i)->dimension() ), prefix , i);
+                        }
+                }
         };
         void log_write ( ) const {
@@ -347 +353 @@
                 if (log_level[logmeans]){
                         for (int i=0; i<particles.length(); i++){
-                                 log_level.store( logmeans, particles(i)->posterior().mean(), i);
+                                log_level.store( logmeans, particles(i)->posterior().mean(), i);
+                        }
+                }
+                if (log_level[logvars]){
+                        for (int i=0; i<particles.length(); i++){
+                                log_level.store( logvars, particles(i)->posterior().variance(), i);
                         }
                 }

library/bdm/stat/emix.cpp

@@ -53 +53 @@
         //non-central moment
         vec mom2 = zeros ( dim );
+        vec mom1 = zeros ( dim );
+        
         for ( int i = 0; i < w.length(); i++ ) {
-                mom2 += w ( i ) * ( component( i )->variance() + pow ( component ( i )->mean(), 2 ) );
+                vec vi=component( i )->variance();
+                vec mi=component ( i )->mean();
+                mom2 += w ( i ) * ( vi + pow ( mi, 2 ) );
+                mom1 += w(i) * mi;
         }
         //central moment
-        return mom2 - pow ( mean(), 2 );
+        return mom2 - pow ( mom1, 2 );
 }
 

library/bdm/stat/exp_family.h

@@ -449 +449 @@
         */
         void from_setting ( const Setting &set );
+        //! see egiw::from_setting
         void to_setting ( Setting& set ) const;
         void validate();