Changeset 283 for bdm/stat

Timestamp:

02/24/09 14:14:01 (15 years ago)

Author:

smidl

Message:

get rid of BMcond + adaptation in doprava/

Location:

bdm/stat

Files:

• libBM.h (modified) (7 diffs)
• libDS.h (modified) (5 diffs)
• libEF.cpp (modified) (2 diffs)
• libEF.h (modified) (7 diffs)

bdm/stat/libBM.h

@@ -126 +126 @@
         int length() const {return len;} ;
         int id ( int at ) const{return ids ( at );};
-        int size ( int at ) const {return RV_SIZES ( ids(at) );};
+        int size ( int at ) const {return RV_SIZES ( ids ( at ) );};
         int time ( int at ) const{return times ( at );};
-        std::string name ( int at ) const {return RV_NAMES ( ids(at) );};
+        std::string name ( int at ) const {return RV_NAMES ( ids ( at ) );};
         void set_time ( int at, int time0 ) {times ( at ) =time0;};
         //!@}
@@ -204 +204 @@
 
         //! access function
-        int _dimy() const{return dimy;}
+        int dimension() const{return dimy;}
 };
 
@@ -259 +259 @@
         //! return expected variance (not covariance!)
         virtual vec variance() const {it_error ( "not implemneted" );return vec ( 0 );};
+        //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
+        virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const {
+                vec mea=mean(); vec std=sqrt(variance()); 
+                lb = mea-2*std; ub=mea+2*std;
+        };
         //!@}
 
@@ -592 +597 @@
 /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
 
+This object represents exact or approximate evaluation of the Bayes rule:
+\f[
+f(\theta_t | d_1,\ldots,d_t) = \frac{f(y_t|\theta_t,\cdot) f(\theta_t|d_1,\ldots,d_{t-1})}{f(y_t|d_1,\ldots,d_{t-1})}
+\f]
+
+Access to the resulting posterior density is via function \c posterior().
+
+As a "side-effect" it also evaluates log-likelihood of the data, which can be accessed via function _ll().
+It can also evaluate predictors of future values of \f$y_t\f$, see functions epredictor() and predictor().
+
+Alternatively, it can evaluate posterior density conditioned by a known constant, \f$ c_t \f$:
+\f[
+f(\theta_t | c_t, d_1,\ldots,d_t) \propto  f(y_t,\theta_t|c_t,\cdot, d_1,\ldots,d_{t-1})
+\f]
+
+The value of \f$ c_t \f$ is set by function condition().
+
 */
 
@@ -606 +628 @@
         //! @{
 
-        BM () :ll ( 0 ),evalll ( true) {};
+        BM () :ll ( 0 ),evalll ( true ), LIDs ( 3 ), opt_L_bounds ( false ) {};
         BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
         //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
-        //! Prototype: \code BM* _copy_(){return new BM(*this);} \endcode
-        virtual BM* _copy_ () {return NULL;};
+        //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
+        virtual BM* _copy_ () const {return NULL;};
         //!@}
 
@@ -634 +656 @@
         //!@}
 
+        //! \name Extension to conditional BM
+        //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
+        //! Alternatively, it can be used for automated connection to DS when the condition is observed
+        //!@{
+
+        //! Name of extension variable
+        RV rvc;
+        //! access function
+        const RV& _rvc() const {return rvc;}
+
+        //! Substitute \c val for \c rvc.
+        virtual void condition ( const vec &val ) {it_error ( "Not implemented!" );};
+
+        //!@}
+
+
         //! \name Access to attributes
         //!@{
@@ -646 +684 @@
         //!@}
 
-};
-
-/*!
-\brief Conditional Bayesian Filter
-
-Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.
-
-This is an interface class used to assure that certain BM has operation \c condition .
-
-*/
-
-class BMcond :public bdmroot {
-protected:
-        //!dimension of the conditioning variable
-        int dimc;
-        //! Identificator of the conditioning variable
-        RV rvc;
-public:
-        //! Substitute \c val for \c rvc.
-        virtual void condition ( const vec &val ) =0;
-        //! Default constructor
-        BMcond ( ) :rvc ( ) {};
-        //! Destructor for future use
-        virtual ~BMcond() {};
-        //! access function
-        const RV& _rvc() const {return rvc;}
+        //! \name Logging of results
+        //!@{
+
+        //! Set boolean options from a string
+        void set_options ( const string &opt ) {
+                opt_L_bounds= ( opt.find ( "logbounds" ) !=string::npos );
+        }
+        //! IDs of storages in loggers
+        ivec LIDs;
+
+        //! Option for logging bounds
+        bool opt_L_bounds;
+        //! Add all logged variables to a logger
+        void log_add ( logger *L, const string &name="" ) {
+                // internal
+                RV r;
+                if ( posterior().isnamed() ) {r=posterior()._rv();}
+                else{r=RV ( "est", posterior().dimension() );};
+
+                // Add mean value
+                LIDs ( 0 ) =L->add ( r,name );
+                if ( opt_L_bounds ) {
+                        LIDs ( 1 ) =L->add ( r,name+"_lb" );
+                        LIDs ( 2 ) =L->add ( r,name+"_ub" );
+                }
+        }
+        void logit ( logger *L ) {
+                L->logit ( LIDs ( 0 ), posterior().mean() );
+                if ( opt_L_bounds ) {
+                        vec ub,lb;
+                        posterior().qbounds(lb,ub);
+                        L->logit ( LIDs ( 1 ), lb );
+                        L->logit ( LIDs ( 2 ), ub );
+                }
+        }
+        //!@}
 };
 

bdm/stat/libDS.h

@@ -28 +28 @@
 */
 class MemDS : public DS {
+        protected:
         //! internal matrix of data
         mat Data;
@@ -45 +46 @@
         void step();
         //!Default constructor
+        MemDS () {};
         MemDS ( mat &Dat, ivec &rowid, ivec &delays );
+};
+
+/*! Read Data Matrix from an IT file
+
+*/
+class FileDS: public MemDS {
+
+public:
+        FileDS ( const string &fname, const string &varname ) :MemDS() {
+                it_file it ( fname );
+                it << Name ( varname ); 
+                it >> Data;
+                time =0;
+                //rowid and delays are ignored
+        }
+        void getdata ( vec &dt ) {
+                it_assert_debug ( dt.length() ==Data.rows(),"" );
+                dt = Data.get_col(time);
+        };
+        void getdata ( vec &dt, const ivec &indeces ){
+                it_assert_debug ( dt.length() ==indeces.length(),"" );
+                vec tmp(indeces.length());
+                tmp = Data.get_col(time);
+                dt = tmp(indeces);
+        };
+        //! returns number of data in the file;
+        int ndat(){return Data.cols();}
 };
 
@@ -96 +125 @@
         { model.set_parameters ( Th0, mu0, sqR0 );};
         //! Set
-        void set_drv(RV &yrv, RV &urv, RV &rrv){
+        void set_drv ( RV &yrv, RV &urv, RV &rrv ) {
                 Rrv = rrv;
                 Urv = urv;
-                dt_size = yrv._dsize()+urv._dsize();
-                
-                RV drv = concat(yrv,urv);
+                dt_size = yrv._dsize() +urv._dsize();
+
+                RV drv = concat ( yrv,urv );
                 Drv = drv;
                 int td = rrv.mint();
-                H.set_size(drv._dsize()*(-td+1));
-                U.set_size(Urv._dsize());
-                for (int i=-1;i>=td;i--){
-                        drv.t(-1);
-                        Drv.add(drv); //shift u1
+                H.set_size ( drv._dsize() * ( -td+1 ) );
+                U.set_size ( Urv._dsize() );
+                for ( int i=-1;i>=td;i-- ) {
+                        drv.t ( -1 );
+                        Drv.add ( drv ); //shift u1
                 }
-                rgrlnk.set_connection(rrv,Drv);
-                
+                rgrlnk.set_connection ( rrv,Drv );
+
                 dtsize = Drv._dsize();
                 utsize = Urv._dsize();
@@ -121 +150 @@
         virtual void log_add ( logger &L ) {
                 //DS::log_add ( L ); too long!!
-                L_dt=L.add ( Drv(0,dt_size),"" );
+                L_dt=L.add ( Drv ( 0,dt_size ),"" );
                 L_ut=L.add ( Urv,"" );
 
@@ -131 +160 @@
         virtual void logit ( logger &L ) {
                 //DS::logit ( L );
-                L.logit( L_dt, H.left(dt_size));
-                L.logit(L_ut, U);
-                
+                L.logit ( L_dt, H.left ( dt_size ) );
+                L.logit ( L_ut, U );
+
                 mat &A =model._A();
                 mat R =model._R();

bdm/stat/libEF.cpp

@@ -183 +183 @@
 };
 
-ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+ivec eEmp::resample (RESAMPLING_METHOD method) {
         ivec ind=zeros_i ( n );
         ivec N_babies = zeros_i ( n );
@@ -268 +268 @@
 }
 
-void eEmp::set_parameters ( const vec &w0, const epdf* epdf0 ) {
+void eEmp::set_statistics ( const vec &w0, const epdf* epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         dim = epdf0->dimension();

bdm/stat/libEF.h

@@ -94 +94 @@
 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
 
-        BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
+        BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 
@@ -348 +348 @@
  \f]
 
+Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
+
  Inverse Gamma can be converted to Gamma using \f[
  x\sim iG(a,b) => 1/x\sim G(a,1/b)
@@ -354 +356 @@
  */
 
-class eigamma : public eEF {
-protected:
-        //!internal egamma
-        egamma eg;
-        //! Vector \f$\alpha\f$
-        vec α
-        //! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG)
-        vec β
+class eigamma : public egamma {
+protected:
 public :
         //! \name Constructors
-        //!@{
-        eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {};
-        eigamma ( const vec &a, const vec &b ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {eg.set_parameters ( a,b );};
-        void set_parameters ( const vec &a, const vec &b ) {eg.set_parameters ( a,b );};
-        //!@}
-
-        vec sample() const {return 1.0/eg.sample();};
-        //! TODO: is it used anywhere?
-//      mat sample ( int N ) const;
-        double evallog ( const vec &val ) const {return eg.evallog ( val );};
-        double lognc () const {return eg.lognc();};
+        //! All constructors are inherited
+        //!@{
+        //!@}
+
+        vec sample() const {return 1.0/egamma::sample();};
         //! Returns poiter to alpha and beta. Potentially dangerous: use with care!
         vec mean() const {return elem_div ( beta,alpha-1 );}
         vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
-        vec& _alpha() {return alpha;}
-        vec& _beta() {return beta;}
 };
 /*
@@ -490 +478 @@
         mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
         //!set mean function
-        void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->_dimy() ), R0 );}
+        void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
         void condition ( const vec &cond ) {mu=g->eval ( cond );};
 };
@@ -686 +674 @@
 
         //! Set samples and weights
-        void set_parameters ( const vec &w0, const epdf* pdf0 );
+        void set_statistics ( const vec &w0, const epdf* pdf0 );
+        //! Set samples and weights
+        void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
         //! Set sample
         void set_samples ( const epdf* pdf0 );
         //! Set sample
-        void set_n ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
+        void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
         //! Potentially dangerous, use with care.
         vec& _w()  {return w;};
@@ -700 +690 @@
         const Array<vec>& _samples() const {return samples;};
         //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
-        ivec resample ( RESAMPLING_METHOD method = SYSTEMATIC );
+        ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
         //! inherited operation : NOT implemneted
         vec sample() const {it_error ( "Not implemented" );return 0;}
@@ -714 +704 @@
                 for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
                 return pom-pow ( mean(),2 );
+        }
+        //! For this class, qbounds are minimum and maximum value of the population!
+        void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const {
+                // lb in inf so than it will be pushed below;
+                lb.set_size(dim);
+                ub.set_size(dim);
+                lb = std::numeric_limits<double>::infinity();
+                ub = -std::numeric_limits<double>::infinity();
+                int j;
+                for ( int i=0;i<n;i++ ) {
+                        for ( j=0;j<dim; j++ ) {
+                                if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
+                                if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
+                        }
+                }
         }
 };