Changeset 96

Timestamp:

05/09/08 17:25:10 (17 years ago)

Author:

smidl

Message:

changes in EF - introduction of the lognc() function

Location:

bdm/stat

Files:

• libBM.h (modified) (2 diffs)
• libEF.cpp (modified) (2 diffs)
• libEF.h (modified) (6 diffs)

bdm/stat/libBM.h

@@ -129 +129 @@
         //! Returns the required moment of the epdf
 //      virtual vec moment ( const int order = 1 );
-        //! Returns a sample from the density, \f$x \sim epdf(rv)\f$
+        //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
         virtual vec sample () const =0;
         //! Compute probability of argument \c val
@@ -224 +224 @@
         RV rv;
         //!Logarithm of marginalized data likelihood.
- double ll;
+        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;

bdm/stat/libEF.cpp

@@ -12 +12 @@
 using std::cout;
 
+vec egiw::sample() const {
+        it_warning("Function not implemented");
+        return vec_1(0.0);
+}
+
+double egiw::evalpdflog( const vec &val ) const {
+        int nPsi = rv.count()-1; // assuming 1dim y
+        double k = nu + nPsi + 2;
+        
+        double r = val(nPsi); //last entry!
+        vec Psi(nPsi+1);
+        Psi(0) = -1.0;
+        Psi.set_subvector(1,val); // fill the rest
+        
+        return -0.5*( k*log(r) + V.qform(Psi)) - lognc();
+}
+
+double egiw::lognc() const{
+        const vec& D = V._D();
+        int nPsi = D.length()-1; // assuming 1dim y
+
+// log(2) = 0.693147180559945286226763983
+// log(pi) = 1.144729885849400163877476189
+return lgamma(0.5*nu) + 0.5*((1.0-nu)*log(D(0)) - V.logdet() + (nu+nPsi)*0.693147180559945286226763983 + nPsi*1.144729885849400163877476189);
+}
+
+vec egiw::mean() const {
+        const mat &L= V._L();
+        const vec &D= V._D();
+        
+        int end = L.rows()-1;
+        vec L0 = L.get_col(0);
+        
+        vec m(D.length());
+        mat iLsub = ltuinv(L(1,end,1,end));
+        m.set_subvector(0,iLsub*L0(1,end));
+        m(end)= D(0)/(nu-2.0);
+        
+        return m;
+}
+
 vec egamma::sample() const {
         vec smp ( rv.count() );
@@ -40 +81 @@
 
 double egamma::evalpdflog ( const vec &val ) const {
+        double res = 0.0; //the rest will be added
+        int i;
+
+        for ( i=0; i<rv.count(); i++ ) {
+                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
+        }
+
+        return res-lognc();
+}
+
+double egamma::lognc() const {
         double res = 0.0; //will be added
         int i;
 
         for ( i=0; i<rv.count(); i++ ) {
-                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) + alpha ( i ) *std::log ( beta ( i ) ) - beta ( i ) *val ( i ) - lgamma ( alpha ( i ) );
+                res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
         }
 

bdm/stat/libEF.h

@@ -37 +37 @@
 
 class eEF : public epdf {
-
 public:
 //      eEF() :epdf() {};
         //! default constructor
         eEF ( const RV &rv ) :epdf ( rv ) {};
+        //! 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 ) {};
@@ -92 +93 @@
         double eval ( const vec &val ) const ;
         double evalpdflog ( const vec &val ) const;
+        double lognc () const;
         vec mean() const {return mu;}
 
@@ -106 +108 @@
         //! access method
         mat getR () {return R.to_mat();}
+};
+
+/*!
+* \brief Gauss-inverse-Wishart density stored in LD form
+
+* More?...
+*/
+class egiw : public eEF {
+protected:
+        //! Extended information matrix of sufficient statistics
+        ldmat V;
+        //! Number of data records (degrees of freedom) of sufficient statistics
+        double nu;
+public:
+        //!Default constructor
+        egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) {
+                it_assert_debug(rv.count()==V.rows(),"Incompatible V0.");
+        }
+
+        vec sample() const;
+        vec mean() const;
+        double evalpdflog ( const vec &val ) const;
+        double lognc () const;
+
+        //Access
+        //! returns a pointer to the internal statistics. Use with Care!
+        ldmat& _V() {return V;}
+        //! returns a pointer to the internal statistics. Use with Care!
+        double& _nu() {return nu;}
+
 };
 
@@ -132 +164 @@
         mat sample ( int N ) const;
         double evalpdflog ( const vec &val ) const;
+        double lognc () const;
         //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
         void _param ( vec* &a, vec* &b ) {a=αb=β};
@@ -255 +288 @@
 class mgamma_fix : public mgamma {
 protected:
+        //! parameter l
         double l;
+        //! reference vector
         vec refl;
 public:
@@ -367 +402 @@
 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
         // 1.83787706640935 = log(2pi)
-        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet() +R.invqform ( mu-val ) );
-};
-
+        return  -0.5* (  +R.invqform ( mu-val ) ) - lognc();
+};
+
+template<class sq_T>
+inline double enorm<sq_T>::lognc () const {
+        // 1.83787706640935 = log(2pi)
+        return -0.5* ( R.cols() * 1.83787706640935 +R.logdet());
+};
 
 template<class sq_T>