Changeset 170 for bdm

Timestamp:

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

Author:

smidl

Message:

Mixtures of EF and related changes to libEF and BM

Location:

bdm

Files:

• estim/arx.cpp (modified) (4 diffs)
• estim/arx.h (modified) (2 diffs)
• estim/libKF.h (modified) (2 diffs)
• estim/libPF.h (modified) (5 diffs)
• math/libDC.cpp (modified) (2 diffs)
• stat/emix.h (modified) (2 diffs)
• stat/libBM.cpp (modified) (1 diff)
• stat/libBM.h (modified) (7 diffs)
• stat/libEF.cpp (modified) (4 diffs)
• stat/libEF.h (modified) (10 diffs)

bdm/estim/arx.cpp

@@ -3 +3 @@
 using namespace itpp;
 
-void ARX::bayes ( const vec &dt ) {
+void ARX::bayes ( const vec &dt, const double w ) {
         double lnc;
 
-        if (frg<1.0) {
-                V*=frg;
-                nu*=frg;
-                if (evalll) {
+        if ( frg<1.0 ) {
+                est.pow ( frg );
+                if ( evalll ) {
                         last_lognc = est.lognc();
                 }
         }
-        V.opupdt ( dt,1.0);
-        nu+=1.0;
+        V.opupdt ( dt,w );
+        nu+=w;
 
         if ( evalll ) {
@@ -20 +19 @@
                 ll = lnc - last_lognc;
                 last_lognc = lnc;
-                tll +=ll;
         }
 }
 
+double ARX::logpred ( const vec &dt ) const {
+        egiw pred ( est );
+        ldmat &V=pred._V();
+        double &nu=pred._nu();
+
+        double lll;
+
+        if ( frg<1.0 ) {
+                pred.pow ( frg );
+                lll = pred.lognc();
+        }
+        else 
+                if(evalll){lll=last_lognc;}else{lll=pred.lognc();}
+
+        V.opupdt ( dt,1.0 );
+        nu+=1.0;
+
+        return pred.lognc()-lll;
+}
+
+BM* ARX::_copy_ ( bool changerv) {
+        ARX* Tmp=new ARX ( *this );
+        if ( changerv ) {Tmp->rv.newids(); Tmp->est._renewrv(Tmp->rv);}
+        return Tmp;
+}
+
+void ARX::set_statistics ( const BMEF* B0 ) {
+        const ARX* A0=dynamic_cast<const ARX*> ( B0 );
+
+        it_assert_debug ( V.rows() ==A0->V.rows(),"ARX::set_statistics Statistics  differ" );
+        set_parameters ( A0->V,A0->nu );
+}
 /*! \brief Return the best structure
 @param Eg a copy of GiW density that is being examined
@@ -48 +78 @@
 
 
-cout << "bb:(" << indeces <<") ll=" << Egll <<endl;
+        cout << "bb:(" << indeces <<") ll=" << Egll <<endl;
 
         //try to remove only one rv
@@ -62 +92 @@
                 tmpll = Eg.lognc()-Eg0.lognc(); // likelihood is difference of norm. coefs.
 
-cout << "i=(" << i <<") ll=" << tmpll <<endl;
-                
+                cout << "i=(" << i <<") ll=" << tmpll <<endl;
+
                 //
                 if ( tmpll > Egll ) { //increase of the likelihood

bdm/estim/arx.h

@@ -34 +34 @@
 Extension for time-variant parameters \f$\theta_t,r_t\f$ may be achived using exponential forgetting (Kulhavy and Zarrop, 1993). In such a case, the forgetting factor \c frg \f$\in <0,1>\f$ should be given in the constructor. Time-invariant parameters are estimated for \c frg = 1.
 */
-class ARX: public BM {
+class ARX: public BMEF {
 protected:
         //! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density
@@ -42 +42 @@
         //! cached value of est.nu
         double &nu;
-        //! forgetting factor 
-        double frg;
-        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
-        double last_lognc;
-        //! total likelihood
-        double tll;
 public:
         //! Full constructor
-        ARX (const RV &rv, const mat &V0, const double &nu0, const double frg0=1.0) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();tll=0.0;};
+        ARX ( const RV &rv, const mat &V0, const double &nu0, const double frg0=1.0 ) : BMEF ( rv,frg0 ),est ( rv,V0,nu0 ), V ( est._V() ), nu ( est._nu() )
+        {last_lognc=est.lognc();};
+
+        //!Copy constructor
+        ARX ( const ARX &A0 ) : BMEF ( A0),est ( rv,A0.V,A0.nu ), V ( est._V() ), nu ( est._nu() ) {};
+
+        //!Auxiliary function
+        BM* _copy_(bool changerv=false);
+        
+//      //! Set parameters given by moments, \c mu (mean of theta), \c R (mean of R) and \c C (variance of theta)
+//      void set_parameters ( const vec &mu, const mat &R, const mat &C, double dfm){};
         //! Set sufficient statistics
-        void set_parameters(const mat &V0, const double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;}
+        void set_parameters ( const ldmat &V0, const double &nu0 )
+        {est._V() =V0;est._nu() =nu0;last_lognc=est.lognc();}
+        void set_statistics ( const BMEF* BM0 );
         //! Returns sufficient statistics
-        void get_parameters(mat &V0, double &nu0){V0=est._V().to_mat(); nu0=est._nu();}
+        void get_parameters ( mat &V0, double &nu0 ) {V0=est._V().to_mat(); nu0=est._nu();}
         //! Here \f$dt = [y_t psi_t] \f$.
-        void bayes ( const vec &dt );
-        epdf& _epdf() {return est;}
+        void bayes ( const vec &dt, const double w );
+        void bayes ( const vec &dt ) {bayes ( dt,1.0 );};
+        const epdf& _epdf() const {return est;}
+        double logpred ( const vec &dt ) const;
+        void flatten (BMEF* B ) {
+                ARX* A=dynamic_cast<ARX*>(B);
+                // nu should be equal to B.nu
+                est.pow ( A->nu/nu);
+                if(evalll){last_lognc=est.lognc();}
+        }
+
         //! Brute force structure estimation.\return indeces of accepted regressors.
-        ivec structure_est(egiw Eg0);
-        //!access function
-        double _tll(){return tll;}
+        ivec structure_est ( egiw Eg0 );
 };
 

bdm/estim/libKF.h

@@ -121 +121 @@
         void bayes ( const vec &dt );
         //!access function
-        epdf& _epdf() {return est;}
+        const epdf& _epdf() const {return est;}
         //!access function
         mat& __K() {return _K;}
@@ -186 +186 @@
         void set_est (vec mu0, mat P0){mu=mu0;P=P0;};
         //!dummy!
-        epdf& _epdf(){return E;};
+        const epdf& _epdf()const{return E;};
 };
 

bdm/estim/libPF.h

@@ -71 +71 @@
                 eEmp &E;
                 vec &_w;
-                Array<epdf*> Coms;
+                Array<const epdf*> Coms;
         public:
                 mpfepdf ( eEmp &E0, const RV &rvc ) :
@@ -80 +80 @@
                 };
 
-                void set_elements ( int &i, double wi, epdf* ep )
+                void set_elements ( int &i, double wi, const epdf* ep )
                 {_w ( i ) =wi; Coms ( i ) =ep;};
 
@@ -112 +112 @@
                 for ( int i=0;i<n;i++ ) {
                         Bms[i] = new BM_T ( BMcond0 ); //copy constructor
-                        epdf& pom=Bms[i]->_epdf();
+                        const epdf& pom=Bms[i]->_epdf();
                         jest.set_elements ( i,1.0/n,&pom );
                 }
@@ -121 +121 @@
 
         void bayes ( const vec &dt );
-        epdf& _epdf() {return jest;}
+        const epdf& _epdf() const {return jest;}
         //! Set postrior of \c rvc to samples from epdf0. Statistics of Bms are not re-computed! Use only for initialization!
         void set_est ( const epdf& epdf0 ) {
@@ -197 +197 @@
                                 delete Bms[i];
                                 Bms[i] = new BM_T ( *Bms[ind ( i ) ] ); //copy constructor
-                                epdf& pom=Bms[i]->_epdf();
+                                const epdf& pom=Bms[i]->_epdf();
                                 jest.set_elements ( i,1.0/n,&pom );
                         }

bdm/math/libDC.cpp

@@ -157 +157 @@
         double x = 0.0, sum;
         int i,j;
-        vec s(v.length());
-        vec S=L*v;
 
         for ( i=0; i<D.length(); i++ ) { //rows of L
@@ -164 +162 @@
                 for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );}
                 x +=D( i )*sum*sum;
-                s(i)=sum;
         };
         return x;

bdm/stat/emix.h

@@ -139 +139 @@
 protected:
         //! Components (epdfs)
-        Array<epdf*> epdfs;
+        Array<const epdf*> epdfs;
         //! Array of indeces
         Array<ivec> rvinds;
 public:
-        eprod ( const Array<epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
+        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
                 bool independent=true;
                 for ( int i=0;i<epdfs.length();i++ ) {
@@ -175 +175 @@
                 return tmp;
         }
+        //!access function
+        const epdf* operator () (int i) const {it_assert_debug(i<epdfs.length(),"wrong index");return epdfs(i);}
 };
 

bdm/stat/libBM.cpp

@@ -180 +180 @@
         return pom;
 }
+
+void RV::newids(){ids=linspace ( RVcounter+1, RVcounter+len ),RVcounter+=len;}
+
+void BM::bayesB(const mat &Data){
+        for(int t=0;t<Data.cols();t++){bayes(Data.get_col(t));}
+}

bdm/stat/libBM.h

@@ -103 +103 @@
         //!access function
         std::string name ( int at ) {return names ( at );};
+        //!Assign unused ids to this rv 
+        void newids();
 };
 
@@ -144 +146 @@
         epdf ( const RV &rv0 ) :rv ( rv0 ) {};
 
-        //! Returns the required moment of the epdf
+//      //! Returns the required moment of the epdf
 //      virtual vec moment ( const int order = 1 );
+        
         //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
         virtual vec sample () const =0;
@@ -156 +159 @@
         virtual double evalpdflog ( const vec &val ) const =0;
 
+        //! Compute log-probability of multiple values argument \c val
+        virtual vec evalpdflog ( const mat &Val ) const {
+                vec x ( Val.cols() );
+                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;}
+                return x;
+        }
+
         //! return expected value
         virtual vec mean() const =0;
@@ -162 +172 @@
         virtual ~epdf() {};
         //! access function, possibly dangerous!
-        RV& _rv() {return rv;}
+        const RV& _rv() const {return rv;}
+        //! modifier function - useful when copying epdfs
+        void _renewrv(const RV &in_rv){rv=in_rv;}
 };
 
@@ -270 +282 @@
 
         //!Default constructor
-        BM ( const RV &rv0 ) :rv ( rv0 ), ll ( 0 ),evalll ( true ) {//Fixme: test rv
+        BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0) {//Fixme: test rv
         };
+        //!Copy constructor
+        BM (const BM &B) : rv(B.rv), ll(B.ll), evalll(B.evalll) {}
 
         /*! \brief Incremental Bayes rule
@@ -278 +292 @@
         virtual void bayes ( const vec &dt ) = 0;
         //! Batch Bayes rule (columns of Dt are observations)
-        void bayes ( mat Dt );
+        virtual void bayesB (const mat &Dt );
         //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
-        virtual epdf& _epdf() =0;
-
+        virtual const epdf& _epdf() const =0;
+
+        //! Evaluates predictive log-likelihood of the given data record
+        //! I.e. marginal likelihood of the data with the posterior integrated out.
+        virtual double logpred(const vec &dt)const{it_error("Not implemented");return 0.0;}
+        
         //! Destructor for future use;
         virtual ~BM() {};
@@ -290 +308 @@
         //!access function 
         void set_evalll(bool evl0){evalll=evl0;}
+        
+        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
+        //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; }
+        virtual BM* _copy_(bool changerv=false){it_error("function _copy_ not implemented for this BM"); return NULL;};
 };
 

bdm/stat/libEF.cpp

@@ -12 +12 @@
 using std::cout;
 
+void BMEF::bayes ( const vec &dt ) {this->bayes ( dt,1.0 );};
+
 vec egiw::sample() const {
         it_warning ( "Function not implemented" );
@@ -17 +19 @@
 }
 
-double egiw::evalpdflog ( const vec &val ) const {
+double egiw::evalpdflog_nn ( const vec &val ) const {
+        int vend = val.length()-1;
 
         if ( xdim==1 ) { //same as the following, just quicker.
-                double r = val ( val.length()-1 ); //last entry!
+                double r = val ( vend ); //last entry!
                 vec Psi ( nPsi+xdim );
                 Psi ( 0 ) = -1.0;
-                Psi.set_subvector ( 1,val ); // fill the rest
-
-                return -0.5* ( nu*log ( r ) + V.qform ( Psi ) /r ) + lognc();
+                Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest
+
+                double Vq=V.qform ( Psi );
+                return -0.5* ( nu*log ( r ) + Vq /r );
         }
         else {
-                int vend = val.length()-1;
                 mat Th= reshape ( val ( 0,nPsi*xdim-1 ),nPsi,xdim );
                 fsqmat R ( reshape ( val ( nPsi*xdim,vend ),xdim,xdim ) );
@@ -34 +37 @@
                 fsqmat iR ( xdim );
                 R.inv ( iR );
-                
-                return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T()*V.to_mat() *Tmp ) ) + lognc();
+
+                return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) );
         }
 }
@@ -55 +58 @@
                      - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
 
-        return nkG+nkW;
+        return -nkG-nkW;
 }
 
 vec egiw::mean() const {
+
+        if ( xdim==1 ) {
+                const mat &L= V._L();
+                const vec &D= V._D();
+                int end = L.rows()-1;
+
+                vec m ( rv.count() );
+                mat iLsub = ltuinv ( L ( xdim,end,xdim,end ) );
+
+                vec L0 = L.get_col ( 0 );
+
+                m.set_subvector ( 0,iLsub*L0 ( 1,end ) );
+                m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2 );
+                return m;
+        } else {
+                mat M;
+                mat R;
+                mean_mat(M,R);
+                return cvectorize (concat_vertical(M,R));
+        }
+
+}
+void egiw::mean_mat(mat &M, mat&R) const {
         const mat &L= V._L();
         const vec &D= V._D();
         int end = L.rows()-1;
-
-        vec m(rv.count());
-        mat iLsub = ltuinv ( L ( xdim,end,xdim,end ) );
+                
+        ldmat ldR ( L ( 0,xdim-1,0,xdim-1 ), D ( 0,xdim-1 ) / ( nu -nPsi -2*xdim -2 ) ); //exp val of R
+        mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) );
         
-        if ( xdim==1 ) {
-                vec L0 = L.get_col ( 0 );
-
-                m.set_subvector ( 0,iLsub*L0 ( 1,end ) );
-                m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2);
-        }
-        else {
-                ldmat ldR(L(0,xdim,0,xdim), D(0,xdim)/( nu -nPsi -2*xdim -2)); //exp val of R
-                mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) );
-                mat M = iLsub*L(xdim,end,0,xdim);
-                
-                m = cvectorize(concat_vertical(M,ldR.to_mat()));
-        }
-        return m;
-}
+        // set mean value
+        M= iLsub*L ( xdim,end,0,xdim-1 );
+        R= ldR.to_mat()  ;
+}
+
+void egiw::pow(double p){V*=p;nu*=p;}
 
 vec egamma::sample() const {

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