Changeset 665 for library/bdm/estim

Timestamp:

10/19/09 22:24:45 (15 years ago)

Author:

smidl

Message:

Compilation and minor extensions

Location:

library/bdm/estim

Files:

• arx.cpp (modified) (4 diffs)
• arx.h (modified) (2 diffs)
• particles.cpp (modified) (5 diffs)
• particles.h (modified) (5 diffs)

library/bdm/estim/arx.cpp

@@ -6 +6 @@
                 
         if ( frg < 1.0 ) {
-                est.pow ( frg );
+                est.pow ( frg ); // multiply V and nu
                 
                 //stabilize
-                ldmat V0(eye(V.rows()));
-                V0*=(1-frg)*1e-3;
+                ldmat V0=alter_est._V(); //$ copy
+                double &nu0=alter_est._nu();
+                V0*=(1-frg);
                 V += V0; //stabilization
-                nu +=(1-frg)*(0.1 + V.rows() + 1* dimx + 2);
-                
-                // recompute loglikelihood of "prior"
+                nu +=(1-frg)*nu0;
+                
+                // recompute loglikelihood of new "prior"
                 if ( evalll ) {
                         last_lognc = est.lognc();
@@ -219 +220 @@
         int rgrlen = rrv->_dsize();
         
+        dimx = ylen;
         set_rv ( *yrv, *rrv );
         
@@ -234 +236 @@
 
         //init
-        mat V0;
-        vec 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" ) )
-                nu0 = rgrlen + ylen + 2;
+        shared_ptr<egiw> pri=UI::build<egiw>(set, "prior", UI::optional);
+        if (pri) {
+                bdm_assert(pri->_dimx()==ylen,"prior is not compatible");
+                bdm_assert(pri->_V().rows()==ylen+rgrlen,"prior is not compatible");
+                est.set_parameters( pri->_dimx(),pri->_V(), pri->_nu());
+        }else{
+                est.set_parameters( ylen, zeros(ylen+rgrlen));
+                set_prior_default(est);
+        }
+                
+        shared_ptr<egiw> alt=UI::build<egiw>(set, "alternative", UI::optional);
+        if (alt) {
+                bdm_assert(alt->_dimx()==ylen,"alternative is not compatible");
+                bdm_assert(alt->_V().rows()==ylen+rgrlen,"alternative is not compatible");
+                alter_est.set_parameters( alt->_dimx(),alt->_V(), alt->_nu());
+        } else {
+                alter_est = est;
+        }
 
         double frg;
@@ -250 +260 @@
 
         set_parameters ( frg );
-        set_statistics ( ylen, V0, nu0 );
         
         //name results (for logging)

library/bdm/estim/arx.h

@@ -162 +162 @@
 
         --- 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
+        prior = {class='egiw',...};                // Prior density, when given default is used instead
+        alternative = {class='egiw',...};          // Alternative density in stabilized estimation, when not given prior is used
+        
         frg = 1.0;                                 // forgetting, default frg=1.0
 
@@ -178 +176 @@
                 bdm_assert(dimx == _yrv._dsize(), "RVs of parameters and regressor do not match");
                 
+        }
+        //! function sets prior and alternative density 
+        void set_prior(const RV &drv, egiw &prior){
+                //TODO check ranges in RV and build prior
+        };
+        //! build default prior and alternative when all values are set
+        void set_prior_default(egiw &prior){
+                //assume 
+                vec dV0(prior._V().rows());
+                dV0.set_subvector(0,prior._dimx()-1, 1.0);
+                dV0.set_subvector(prior._dimx(),dV0.length()-1, 1e-5);
+                
+                prior.set_parameters(prior._dimx(),ldmat(dV0));
         }
 };

library/bdm/estim/particles.cpp

@@ -5 +5 @@
 using std::endl;
 
-void PF::bayes_gensmp(){
-        #pragma parallel for
-        for (int i = 0; i < n; i++ ) {
-                _samples ( i ) = par->samplecond ( _samples ( i ) );
-                lls(i) = 0;
+void PF::bayes_gensmp(const vec &ut){
+        if (ut.length()>0){
+                vec cond(par->dimensionc());
+                cond.set_subvector(par->dimension(),ut);
+                #pragma parallel for
+                for (int i = 0; i < n; i++ ) {
+                        cond.set_subvector(0,_samples ( i ));
+                        _samples ( i ) = par->samplecond ( cond );
+                        lls(i) = 0;
+                }
+        } else {                
+                #pragma parallel for
+                for (int i = 0; i < n; i++ ) {
+                        _samples ( i ) = par->samplecond ( _samples ( i ) );
+                        lls(i) = 0;
+                }
         }
 }
@@ -28 +39 @@
                 }
                 sw = sum(_w);
-                if (!std::isfinite(sw)) {
+                if (!std::isfinite(sw) || sw==0.0) {
                         bdm_error("Particle filter is lost; no particle is good enough.");
                 }
@@ -36 +47 @@
 
 void PF::bayes ( const vec &dt ) {
+        vec yt=dt.left(obs->dimension());
+        vec ut=dt.get(obs->dimension(),dt.length()-1);
+        
         int i;
         // generate samples - time step
-        bayes_gensmp();
+        bayes_gensmp(ut);
         // weight them - data step
         for ( i = 0; i < n; i++ ) {
-                lls ( i ) += obs->evallogcond ( dt, _samples ( i ) ); //+= because lls may have something from gensmp!
+                lls ( i ) += obs->evallogcond ( yt, _samples ( i ) ); //+= because lls may have something from gensmp!
         }
 
@@ -61 +75 @@
 
 void MPF::bayes ( const vec &dt ) {     
-        // follows PF::bayes in most places!!!
-        
+        // follows PF::bayes in most places!!!  
         int i;
         int n=pf->__w().length();
@@ -68 +81 @@
         
         // generate samples - time step
-        pf->bayes_gensmp();
+        pf->bayes_gensmp(vec(0));
         // weight them - data step
         #pragma parallel for

library/bdm/estim/particles.h

@@ -93 +93 @@
         }
         //! bayes I - generate samples and add their weights to lls
-        virtual void bayes_gensmp();
+        virtual void bayes_gensmp(const vec &ut);
         //! bayes II - compute weights of the 
         virtual void bayes_weights();
@@ -124 +124 @@
         */
         void from_setting(const Setting &set){
+                BM::from_setting(set);
                 par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
                 obs = UI::build<mpdf>(set,"observation_pdf",UI::compulsory);
@@ -137 +138 @@
                 
                 u.add(obs_u); // join both u, and check if they do not overlap
-                
+
                 set_drv(concat(obs->_rv(),u) );
         }
@@ -163 +164 @@
                         res_threshold=0.5;
                 }
+                validate();
         }
         //! load prior information from set and set internal structures accordingly
@@ -172 +174 @@
                         est=*test_emp;
                 } else {
-                        int n;
+                        //int n;
                         if (!UI::get(n,set,"n",UI::optional)){n=10;}
                         // sample from prior
                         set_statistics(ones(n)/n, *pri);
                 }
+                n = est._w().length();
                 //validate();
         }