Changeset 28 for bdm/stat

Timestamp:

02/22/08 16:40:12 (16 years ago)

Author:

smidl

Message:

prelozitelna verze

Location:

bdm/stat

Files:

• libBM.h (modified) (5 diffs)
• libEF.h (modified) (7 diffs)

bdm/stat/libBM.h

@@ -84 +84 @@
         //! access function
         int _dimy()const{return dimy;}
+
+                //! Destructor for future use;
+                virtual ~fnc(){};
 };
 
-//! Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
-class BM {
-public:
-        //!Logarithm of marginalized data likelihood.
-        double ll;
-
-        //!Default constructor
-        BM(){ll=0;};
-        
-        /*! \brief Incremental Bayes rule
-        @param dt vector of input data
-        @param evall If true, the filter will compute likelihood of the data record and store it in \c ll
-        */
-        virtual void bayes ( const vec &dt, bool evall=true ) = 0;
-        //! Batch Bayes rule (columns of Dt are observations)
-        void bayes ( mat Dt );
-};
 
 //! Probability density function with numerical statistics, e.g. posterior density.
@@ -108 +94 @@
         RV rv;
 public:
+        //!default constructor
+        epdf():rv(ivec(0)){};
         //! Returns the required moment of the epdf
 //      virtual vec moment ( const int order = 1 );
@@ -113 +101 @@
         virtual vec sample ()=0;
         //! Compute probability of argument \c val
-        virtual double eval(const vec &val)
-        {
-                return 0.0;
-        };
+        virtual double eval(const vec &val){return 0.0;};
+        //! Compute log-probability of argument \c val
+        virtual double evalpdflog(const vec &val){return 0.0;};
+
+                //! Destructor for future use;
+                virtual ~epdf(){};
 };
 
@@ -130 +120 @@
 //      virtual fnc moment ( const int order = 1 );
         //! Returns a sample from the density conditioned on \c cond, $x \sim epdf(rv|cond)$
-        virtual vec samplecond (vec &cond, double lik){};
+        virtual vec samplecond (vec &cond, double lik){return vec(0);};
         virtual void condition (vec &cond){};
+
+                //! Destructor for future use;
+                virtual ~mpdf(){};
 };
 
@@ -164 +157 @@
         //! Moves from $t$ to $t+1$, i.e. perfroms the actions and reads response of the system. 
         void step();
+                
 };
 
+/*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.
+        
+*/
+class BM {
+public:
+        //!Logarithm of marginalized data likelihood.
+        double ll;
+        //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save time.
+        bool evalll;
+
+        //!Default constructor
+        BM():ll(0),evalll(true){};
+        
+        /*! \brief Incremental Bayes rule
+        @param dt vector of input data
+        */
+        virtual void bayes ( const vec &dt) = 0;
+        //! Batch Bayes rule (columns of Dt are observations)
+        void bayes ( mat Dt );
+        //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
+        virtual epdf* _epdf()=0;
+
+                //! Destructor for future use;
+                virtual ~BM(){};
+};
 
 #endif // BM_H

bdm/stat/libEF.h

@@ -26 +26 @@
 * More?...
 */
+
 class eEF : public epdf {
 
 public:
-        virtual void tupdate( double phi, mat &vbar, double nubar ) {};
-        virtual void dupdate( mat &v,double nu=1.0 ) {};
+        //! default constructor
+        eEF() :epdf() {};
+
+        virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
+
+        virtual void dupdate ( mat &v,double nu=1.0 ) {};
 };
 
@@ -45 +50 @@
 */
 template<class sq_T>
+
 class enorm : public eEF {
+protected:
+        //! mean value
+        vec mu;
+        //! Covariance matrix in decomposed form
+        sq_T R;
+        //! Cache:  _iR = inv(R);
+        sq_T _iR;
+        //! indicator if \c _iR is chached
+        bool cached;
+        //! Random number generator for sampling; Fixme: what about *static* correct ?
+        Normal_RNG RNG;
+        //! dimension (redundant from rv.count() for easier coding )
         int dim;
-        vec mu;
-        sq_T R;
-public:
-        enorm( RV &rv, vec &mu, sq_T &R );
-        enorm();
+public:
+        enorm() :eEF() {};
+
+        enorm ( RV &rv );
+        void set_parameters ( const vec &mu,const sq_T &R );
         //! tupdate in exponential form (not really handy)
-        void tupdate( double phi, mat &vbar, double nubar );
-        void dupdate( mat &v,double nu=1.0 );
-        //!tupdate used in KF
-        void tupdate();
-        //!dupdate used in KF
-        double dupdate();
-         
+        void tupdate ( double phi, mat &vbar, double nubar );
+        void dupdate ( mat &v,double nu=1.0 );
+
         vec sample();
-        mat sample(int N);
-        double eval( const vec &val );
-        Normal_RNG RNG;
+        mat sample ( int N );
+        double eval ( const vec &val );
+        double evalpdflog ( const vec &val );
+
+//Access methods
+        //! returns a pointer to the internal mean value. Use with Care!
+        vec* _mu() {return &mu;}
+
+        //! returns pointers to the internal variance and its inverse. Use with Care!
+        void _R ( sq_T* &pR, sq_T* &piR ) {
+                pR=&R;
+                piR=&_iR;
+        }
+
+        //! set cache as inconsistent
+        void _cached ( bool what ) {cached=what;}
 };
 
@@ -70 +97 @@
 */
 template<class sq_T>
+
 class mlnorm : public mEF {
         enorm<sq_T> epdf;
@@ -75 +103 @@
 public:
         //! Constructor
-        mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R );
+        mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R );
         //!Generate one sample of the posterior
-        vec samplecond( vec &cond, double &lik );
-        mat samplecond( vec &cond, vec &lik, int n );
-        void condition( vec &cond );
+        vec samplecond ( vec &cond, double &lik );
+        mat samplecond ( vec &cond, vec &lik, int n );
+        void condition ( vec &cond );
 };
 
@@ -85 +113 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
-        dim = rv.count();
+enorm<sq_T>::enorm ( RV &rv ) :eEF(), mu ( rv.count() ),R(rv.count()),_iR(rv.count()),cached ( false ),dim ( rv.count() ) {};
+
+template<class sq_T>
+void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
+//Fixme test dimensions of mu0 and R0;
         mu = mu0;
         R = R0;
-};
-
-template<class sq_T>
-void enorm<sq_T>::dupdate( mat &v, double nu ) {
+        if(_iR.rows()!=R.rows()) _iR=R; // memory allocation!
+        R.inv ( _iR ); //update cache
+        cached=true;
+};
+
+template<class sq_T>
+void enorm<sq_T>::dupdate ( mat &v, double nu ) {
         //
 };
 
 template<class sq_T>
-void enorm<sq_T>::tupdate( double phi, mat &vbar, double nubar ) {
+void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
         //
 };
@@ -103 +137 @@
 template<class sq_T>
 vec enorm<sq_T>::sample() {
-        vec x( dim );
-        RNG.sample_vector( dim,x );
-        vec smp = R.sqrt_mult( x );
+        vec x ( dim );
+        RNG.sample_vector ( dim,x );
+        vec smp = R.sqrt_mult ( x );
 
         smp += mu;
@@ -112 +146 @@
 
 template<class sq_T>
-mat enorm<sq_T>::sample( int N ) {
-        mat X( dim,N );
-        vec x( dim );
+mat enorm<sq_T>::sample ( int N ) {
+        mat X ( dim,N );
+        vec x ( dim );
         vec pom;
         int i;
+
         for ( i=0;i<N;i++ ) {
-                RNG.sample_vector( dim,x );
-                pom = R.sqrt_mult( x );
+                RNG.sample_vector ( dim,x );
+                pom = R.sqrt_mult ( x );
                 pom +=mu;
-                X.set_col( i, pom);
+                X.set_col ( i, pom );
         }
+
         return X;
 };
 
 template<class sq_T>
-double enorm<sq_T>::eval( const vec &val ) {
-        return 0.0;     
-};
-
-
-template<class sq_T>
-enorm<sq_T>::enorm() {};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R ) {
+double enorm<sq_T>::eval ( const vec &val ) {
+        return exp ( evalpdflog ( val ) );
+};
+
+template<class sq_T>
+double enorm<sq_T>::evalpdflog ( const vec &val ) {
+        if ( !cached ) {R.inv ( _iR );cached=true;}
+
+        return -0.5* ( R.cols() *0.79817986835811504957 +R.logdet() +_iR.qform ( mu-val ) );
+};
+
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( RV &rv,RV &rvc, mat &A, sq_T &R ) {
         int dim = rv.count();
-        vec mu( dim );
-
-        epdf = enorm<sq_T>( rv,mu,R );
-}
-
-template<class sq_T>
-vec mlnorm<sq_T>::samplecond( vec &cond, double &lik ) {
-        this->condition( cond );
+        vec mu ( dim );
+
+        epdf = enorm<sq_T> ( rv,mu,R );
+}
+
+template<class sq_T>
+vec mlnorm<sq_T>::samplecond ( vec &cond, double &lik ) {
+        this->condition ( cond );
         vec smp = epdf.sample();
-        lik = epdf.eval( smp );
+        lik = epdf.eval ( smp );
         return smp;
 }
 
 template<class sq_T>
-mat mlnorm<sq_T>::samplecond( vec &cond, vec &lik, int n ) {
+mat mlnorm<sq_T>::samplecond ( vec &cond, vec &lik, int n ) {
         int i;
         int dim = rv.count();
-        mat Smp( dim,n );
-        vec smp( dim );
-        this->condition( cond );
+        mat Smp ( dim,n );
+        vec smp ( dim );
+        this->condition ( cond );
+
         for ( i=0; i<dim; i++ ) {
                 smp = epdf.sample();
-                lik( i ) = epdf.eval( smp );
-                Smp.set_col( i ,smp );
+                lik ( i ) = epdf.eval ( smp );
+                Smp.set_col ( i ,smp );
         }
+
         return Smp;
 }
 
 template<class sq_T>
-void mlnorm<sq_T>::condition( vec &cond ) {
+void mlnorm<sq_T>::condition ( vec &cond ) {
         epdf.mu = A*cond;
 //R is already assigned;
 }
 
+///////////
+
+//#ifdef HAVE_EXTERN_TEMPLATE
+
+//  extern template class enorm<fsqmat>;
+//  extern template class enorm<ldmat>;
+//#endif
+
 #endif //EF_H