Changeset 665

Timestamp:

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

Author:

smidl

Message:

Compilation and minor extensions

Location:

library/bdm

Files:

• base/bdmbase.h (modified) (5 diffs)
• base/datasources.cpp (modified) (1 diff)
• base/user_info.h (modified) (1 diff)
• bdmroot.h (modified) (1 diff)
• estim/arx.cpp (modified) (4 diffs)
• estim/arx.h (modified) (2 diffs)
• estim/particles.cpp (modified) (5 diffs)
• estim/particles.h (modified) (5 diffs)
• math/functions.h (modified) (1 diff)
• stat/exp_family.cpp (modified) (1 diff)
• stat/exp_family.h (modified) (5 diffs)

library/bdm/base/bdmbase.h

@@ -149 +149 @@
         friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
 
-        string to_string() {ostringstream o; o << this; return o.str();}
+        string to_string() const {ostringstream o; o << *this; return o.str();}
         
         //! total size of a random variable
@@ -303 +303 @@
         //! Length of the output vector
         int dimy;
+        //! Length of the input vector
+        int dimc;
 public:
         //!default constructor
@@ -317 +319 @@
         int dimension() const {
                 return dimy;
+        }
+        //! access function
+        int dimensionc() const {
+                return dimc;
         }
 };
@@ -671 +677 @@
                 return downsize;
         }
+        //! for future use
+        virtual ~datalink(){}
 };
 
@@ -1219 +1227 @@
                 }
                 string opt;
-                if ( !UI::get ( opt, set, "options", UI::optional ) ) {
+                if ( UI::get ( opt, set, "options", UI::optional ) ) {
                         set_options ( opt );
                 }

library/bdm/base/datasources.cpp

@@ -118 +118 @@
         time = 0;
         //rowid and delays are ignored
-
+        rowid = linspace(0,Data.rows()-1);
         set_drv ( *rvtmp, RV() );
 }

library/bdm/base/user_info.h

@@ -575 +575 @@
         //!@}
 
+        //! for future use
+        virtual ~UI(){}
 };
 

library/bdm/bdmroot.h

@@ -33 +33 @@
 
         //! This method returns a basic info about the current instance
-        virtual string to_string() {
+        virtual string to_string() const {
                 return "";
         }

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();
         }

library/bdm/math/functions.h

@@ -84 +84 @@
         vec eval ( const vec &cond ) {
                 bdm_assert_debug ( cond.length() == ( dimx + dimu ), "linfn::eval Wrong cond." );
-                return eval ( cond ( 0, dimx - 1 ), cond ( dimx, dimx + dimu - 1 ) );//-1 = end (in matlab)
+                if (dimu>0){
+                        return eval ( cond ( 0, dimx - 1 ), cond ( dimx, dimx + dimu - 1 ) );//-1 = end (in matlab)
+                } else {
+                        return eval ( cond ( 0, dimx - 1 ), vec(0) );//-1 = end (in matlab)
+                }
+                                
         }
 

library/bdm/stat/exp_family.cpp

@@ -82 +82 @@
 
 //      bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
-if (-nkG - nkW==Inf){
-        cout << "??" <<endl;
-}
+        if (-nkG - nkW==Inf){
+                cout << "??" <<endl;
+        }
         return -nkG - nkW;
 }

library/bdm/stat/exp_family.h

@@ -244 +244 @@
                 double& _nu()  {return nu;}
                 const double& _nu() const {return nu;}
+                const int & _dimx() const {return dimx;}
                 /*! Create Gauss-inverse-Wishart density 
                 \f[ f(rv) = GiW(V,\nu) \f]
@@ -250 +251 @@
                 class = 'egiw';
                 V     = [];               // square matrix
+                dV    = [];               // vector of diagonal of V (when V not given)
                 nu    = [];               // scalar \nu ((almost) degrees of freedom)
                                                                   // when missing, it will be computed to obtain proper pdf
@@ -260 +262 @@
                 void from_setting (const Setting &set) {
                         epdf::from_setting(set);
-                        if (!UI::get (nu, set, "nu", UI::compulsory)) {nu=-1;}
                         UI::get (dimx, set, "dimx", UI::compulsory);
+                        if (!UI::get (nu, set, "nu", UI::optional)) {nu=-1;}
                         mat V;
-                        UI::get (V, set, "V", UI::compulsory);
-                        set_parameters (dimx, V, nu);
+                        if (!UI::get (V, set, "V", UI::optional)){
+                                vec dV;
+                                UI::get (dV, set, "dV", UI::compulsory);
+                                set_parameters (dimx, ldmat(dV), nu);
+                                
+                        } else {
+                                set_parameters (dimx, V, nu);
+                        }
+                }
+                void validate(){
+                        // check sizes, rvs etc.
                 }
                 //!@}
@@ -594 +605 @@
 
                 double evallog (const vec &val) const  {
-                        if (any (val < low) && any (val > high)) {return inf;}
+                        if (any (val < low) && any (val > high)) {return -inf;}
                         else return lnk;
                 }
@@ -632 +643 @@
 UIREGISTER(euni);
 
+//! Uniform density with conditional mean value
+class mguni : public mpdf_internal<euni>{
+        //! function of the mean value
+        shared_ptr<fnc> mean;
+        //! distance from mean to both sides
+        vec delta;
+        public:
+                void condition(const vec &cond){
+                        vec mea=mean->eval(cond);
+                        iepdf.set_parameters(mea-delta,mea+delta);
+                }
+                //! load from
+                void from_setting(const Setting &set){
+                        mpdf::from_setting(set); //reads rv and rvc
+                        UI::get(delta,set,"delta",UI::compulsory);
+                        mean = UI::build<fnc>(set,"mean",UI::compulsory);
+                        
+                        iepdf.set_parameters(-delta,delta);
+                        dimc = mean->dimensionc();
+                        validate();
+                }
+};
+UIREGISTER(mguni);
 /*!
  \brief Normal distributed linear function with linear function of mean value;