Changeset 170 for bdm/stat/libEF.h

Timestamp:

09/24/08 13:07:50 (16 years ago)

Author:

smidl

Message:

Mixtures of EF and related changes to libEF and BM

Files:

• bdm/stat/libEF.h (modified) (10 diffs)

bdm/stat/libEF.h

@@ -41 +41 @@
         //! default constructor
         eEF ( const RV &rv ) :epdf ( rv ) {};
-        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$ 
-        virtual double lognc()const =0;
+        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
+        virtual double lognc() const =0;
         //!TODO decide if it is really needed
-        virtual void tupdate ( double phi, mat &vbar, double nubar ) {};
-        //!TODO decide if it is really needed
-        virtual void dupdate ( mat &v,double nu=1.0 ) {};
+        virtual void dupdate ( mat &v ) {it_error ( "Not implemneted" );};
+        //!Evaluate normalized log-probability
+        virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemneted" );return 0.0;};
+        //!Evaluate normalized log-probability
+        virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
+        //!Evaluate normalized log-probability for many samples
+        virtual vec evalpdflog ( const mat &Val ) const {
+                vec x ( Val.cols() );
+                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog_nn ( Val.get_col ( i ) ) ;}
+                return x-lognc();
+        }
+        //!Power of the density, used e.g. to flatten the density
+        virtual void pow ( double p ) {it_error ( "Not implemented" );};
 };
 
@@ -60 +70 @@
         //! Default constructor
         mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
+};
+
+//! Estimator for Exponential family
+class BMEF : public BM {
+protected:
+        //! forgetting factor
+        double frg;
+        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
+        double last_lognc;
+public:
+        //! Default constructor
+        BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
+        //! Copy constructor
+        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
+        //!get statistics from another model
+        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
+        //! Weighted update of sufficient statistics (Bayes rule)
+        virtual void bayes ( const vec &data, const double w ) {};
+        //original Bayes
+        void bayes ( const vec &dt );
+        //!Flatten the posterior
+        virtual void flatten ( BMEF * B) {it_error ( "Not implemented" );}
 };
 
@@ -99 +131 @@
         //! returns a pointer to the internal mean value. Use with Care!
         vec& _mu() {return mu;}
-        
+
         //! access function
-        void set_mu(const vec mu0) { mu=mu0;}
+        void set_mu ( const vec mu0 ) { mu=mu0;}
 
         //! returns pointers to the internal variance and its inverse. Use with Care!
@@ -128 +160 @@
 public:
         //!Default constructor, assuming
-        egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) {
-                xdim = rv.count()/V.rows();
-                it_assert_debug(rv.count()==xdim*V.rows(),"Incompatible V0.");
+        egiw ( RV rv, mat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
+                xdim = rv.count() /V.rows();
+                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
+                nPsi = V.rows()-xdim;
+        }
+        //!Full constructor for V in ldmat form
+        egiw ( RV rv, ldmat V0, double nu0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
+                xdim = rv.count() /V.rows();
+                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
                 nPsi = V.rows()-xdim;
         }
@@ -136 +174 @@
         vec sample() const;
         vec mean() const;
+        void mean_mat ( mat &M, mat&R ) const;
         //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
-        double evalpdflog ( const vec &val ) const;
+        double evalpdflog_nn ( const vec &val ) const;
         double lognc () const;
 
@@ -145 +184 @@
         //! returns a pointer to the internal statistics. Use with Care!
         double& _nu() {return nu;}
-
+        void pow ( double p );
+};
+
+/*! \brief Dirichlet posterior density
+Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
+\f[
+f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n]x_i^(\beta_i-1)
+\f]
+where \f$\gamma=\sum_i beta_i\f$.
+*/
+class eDirich: public eEF {
+protected:
+        //!sufficient statistics
+        vec beta;
+public:
+        //!Default constructor
+        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
+        //! Copy constructor
+        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
+        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
+        vec mean() const {return beta/sum ( beta );};
+        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
+        double evalpdflog_nn ( const vec &val ) const {return ( beta-1 ) *log ( val );};
+        double lognc () const {
+                double gam=sum ( beta );
+                double lgb=0.0;
+                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
+                return lgb-lgamma ( gam );
+        };
+        //!access function
+        vec& _beta() {return beta;}
+};
+
+//! Estimator for Multinomial density
+class multiBM : public BMEF {
+protected:
+        //! Conjugate prior and posterior
+        eDirich est;
+        vec &beta;
+public:
+        //!Default constructor
+        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {last_lognc=est.lognc();}
+        //!Copy constructor
+        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
+
+        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
+        void bayes ( const vec &dt ) {
+                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
+                beta+=dt;
+                if ( evalll ) {ll=est.lognc()-last_lognc;}
+        }
+        double logpred ( const vec &dt ) const {
+                eDirich pred ( est );
+                vec &beta = pred._beta();
+                
+                double lll;
+                if ( frg<1.0 )
+                        {beta*=frg;lll=pred.lognc();}
+                else
+                        if ( evalll ) {lll=last_lognc;}
+                        else{lll=pred.lognc();}
+
+                beta+=dt;
+                return pred.lognc()-lll;
+        }
+        void flatten (BMEF* B ) {
+                eDirich* E=dynamic_cast<eDirich*>(B);
+                // sum(beta) should be equal to sum(B.beta)
+                const vec &Eb=E->_beta();
+                est.pow ( sum(beta)/sum(Eb) );
+                if(evalll){last_lognc=est.lognc();}
+        }
+        const epdf& _epdf() const {return est;};
+        //!access funct
 };
 
@@ -151 +263 @@
  \brief Gamma posterior density
 
- Multvariate Gamma density as product of independent univariate densities.
+ Multivariate Gamma density as product of independent univariate densities.
  \f[
  f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
@@ -213 +325 @@
         double evalpdflog ( const vec &val ) const  {return lnk;}
         vec sample() const {
-                vec smp ( rv.count() ); 
-                #pragma omp critical
+                vec smp ( rv.count() );
+#pragma omp critical
                 UniRNG.sample_vector ( rv.count(),smp );
-                return low+elem_mult(distance,smp);
+                return low+elem_mult ( distance,smp );
         }
         //! set values of \c low and \c high
@@ -408 +520 @@
 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
         // 1.83787706640935 = log(2pi)
-        return  -0.5* (  +R.invqform ( mu-val ) ) - lognc();
+        return  -0.5* ( +R.invqform ( mu-val ) ) - lognc();
 };
 
@@ -414 +526 @@
 inline double enorm<sq_T>::lognc () const {
         // 1.83787706640935 = log(2pi)
-        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet());
-};
-
-template<class sq_T>
-mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu(epdf._mu()) { ep =&epdf;
+        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
+};
+
+template<class sq_T>
+mlnorm<sq_T>::mlnorm ( RV &rv0,RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
+        ep =&epdf;
 }