Changeset 168

Timestamp:

09/18/08 19:54:09 (16 years ago)

Author:

smidl

Message:

Work on mixtures of EF, small changes

Files:

• CMakeLists.txt (modified) (1 diff)
• bdm/CMakeLists.txt (modified) (1 diff)
• bdm/estim/arx.h (modified) (1 diff)
• bdm/math/chmat.cpp (modified) (1 diff)
• bdm/math/chmat.h (modified) (1 diff)
• bdm/math/libDC.cpp (modified) (2 diffs)
• bdm/math/libDC.h (modified) (3 diffs)
• bdm/stat/emix.cpp (modified) (3 diffs)
• bdm/stat/emix.h (modified) (2 diffs)
• bdm/stat/libBM.cpp (modified) (1 diff)
• bdm/stat/libBM.h (modified) (2 diffs)
• bdm/stat/libEF.cpp (modified) (7 diffs)
• bdm/stat/libEF.h (modified) (2 diffs)
• bdm/stat/merger.h (modified) (3 diffs)
• tests/merger_test.cpp (modified) (1 diff)

CMakeLists.txt

@@ +1 @@
+# Compulsory line
+cmake_minimum_required(VERSION 2.6)
+
 # The name of our project is "BDM".  CMakeLists files in this project can
 # refer to the root source directory of the project as ${BDM_SOURCE_DIR} and

bdm/CMakeLists.txt

@@ -3 +3 @@
 SET(BdmMath math/libDC.cpp math/libDC.h math/chmat.cpp math/chmat.h)
 SET(BdmStat stat/libDS.cpp stat/libDS.h stat/libFN.cpp stat/libFN.h stat/libBM.cpp stat/libBM.h stat/libEF.cpp stat/libEF.h stat/loggers.cpp stat/loggers.h stat/emix.cpp stat/emix.h stat/merger.cpp)
-SET(BdmEstim estim/libKF.cpp estim/libKF.h estim/libPF.cpp estim/libPF.h estim/arx.cpp estim/arx.h)
+SET(BdmEstim estim/libKF.cpp estim/libKF.h estim/libPF.cpp estim/libPF.h estim/arx.cpp estim/arx.h estim/mixef.cpp estim/mixef.h)
 SET(BdmUI userinfo.cpp userinfo.h)
 

bdm/estim/arx.h

@@ -50 +50 @@
 public:
         //! Full constructor
-        ARX (RV &rv, mat &V0, double &nu0, 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) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();tll=0.0;};
         //! Set sufficient statistics
-        void set_parameters(mat &V0, double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;}
+        void set_parameters(const mat &V0, const double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;}
         //! Returns sufficient statistics
         void get_parameters(mat &V0, double &nu0){V0=est._V().to_mat(); nu0=est._nu();}

bdm/math/chmat.cpp

@@ -17 +17 @@
         Ch = R ( 0, Ch.rows()-1, 0, Ch.cols()-1 );
 };
-mat chmat::to_mat() {mat F=Ch.T() *Ch;return F;};
+mat chmat::to_mat() const {mat F=Ch.T() *Ch;return F;};
 void chmat::mult_sym ( const mat &C ) {
         it_error ( "not implemented" );

bdm/math/chmat.h

@@ -32 +32 @@
 
         void opupdt ( const vec &v, double w );
-        mat to_mat();
+        mat to_mat() const;
         void mult_sym ( const mat &C );
         void mult_sym ( const mat &C , chmat &U ) const;

bdm/math/libDC.cpp

@@ -7 +7 @@
 
 void fsqmat::opupdt ( const vec &v, double w ) {M+=outer_product ( v,v*w );};
-mat fsqmat::to_mat()  {return M;};
+mat fsqmat::to_mat() const {return M;};
 void fsqmat::mult_sym ( const mat &C) {M=C *M*C.T();};
 void fsqmat::mult_sym_t ( const mat &C) {M=C.T() *M*C;};
@@ -78 +78 @@
 }
 
-mat ldmat::to_mat() {
+mat ldmat::to_mat() const {
         int dim = D.length();
         mat V( dim, dim );

bdm/math/libDC.h

@@ -44 +44 @@
                 */
 
-                virtual mat to_mat() =0;
+                virtual mat to_mat() const =0;
 
                 /*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$
@@ -115 +115 @@
         public:
                 void opupdt ( const vec &v, double w );
-                mat to_mat() ;
+                mat to_mat() const;
                 void mult_sym ( const mat &C);
                 void mult_sym_t ( const mat &C);
@@ -198 +198 @@
 
                 void opupdt ( const vec &v, double w );
-                mat to_mat();
+                mat to_mat() const;
                 void mult_sym ( const mat &C);
                 void mult_sym_t ( const mat &C);

bdm/stat/emix.cpp

@@ -4 +4 @@
 
 void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0 ) {
-        w = w0;
+        w = w0/sum(w0);
         int i;
         for ( i=0;i<w.length();i++ ) {
@@ -38 +38 @@
                 };
                 
-                independent = true;
+//              independent = true;
                 //test rvc of mpdfs and fill rvinds
                 for ( i = 0;i < n;i++ ) {
@@ -48 +48 @@
                         rvcinds ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
                         //
-                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
-                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
+/*                      if ( rvcinrv ( i ).length() >0 ) {independent = false;}
+                        if ( rvcinds ( i ).length() >0 ) {independent = false;}*/
                 }
         };

bdm/stat/emix.h

@@ -82 +82 @@
         //! Indeces of rvc in common rvc
         Array<ivec> rvcinds;
-        //! Indicate independence of its factors
-        bool independent;
+//      //! Indicate independence of its factors
+//      bool independent;
 //      //! Indicate internal creation of mpdfs which must be destroyed
 //      bool intermpdfs;
 public:
-        /*!\brief Constructor from list of mFacs, 
+        /*!\brief Constructor from list of mFacs,
         Additional parameter overlap is left for future use. Do not set to true for mprod.
         */
-        mprod ( Array<mpdf*> mFacs, bool overlap=false ); 
+        mprod ( Array<mpdf*> mFacs, bool overlap=false );
 
         double evalpdflog ( const vec &val ) const {
@@ -122 +122 @@
                         set_subvector ( smp,rvinds ( i ), smpi );
                         // add ith likelihood contribution
-                        ll+=epdfs(i)->evalpdflog(smpi);
+                        ll+=epdfs ( i )->evalpdflog ( smpi );
                 }
                 return smp;
         }
         mat samplecond ( const vec &cond, vec &ll, int N ) {
-                mat Smp(rv.count(),N);
-                for(int i=0;i<N;i++){Smp.set_col(i,samplecond(cond,ll(i)));}
+                mat Smp ( rv.count(),N );
+                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond,ll ( i ) ) );}
                 return Smp;
         }
-//      vec mean() const {
-//              vec tmp ( rv.count() );
-//              if ( independent ) {
-//                      for ( int i=0;i<n;i++ ) {
-//                              vec pom = epdfs ( i )->mean();
-//                              set_subvector ( tmp,rvinds ( i ), pom );
-//                      }
-//              }
-//              else {
-//                      int N=50*rv.count();
-//                      it_warning ( "eprod.mean() computed by sampling" );
-//                      tmp = zeros ( rv.count() );
-//                      for ( int i=0;i<N;i++ ) { tmp += sample();}
-//                      tmp /=N;
-//              }
-//              return tmp;
-//      }
 
         ~mprod() {};
 };
+
+//! Product of independent epdfs. For dependent pdfs, use mprod.
+class eprod: public epdf {
+protected:
+        //! Components (epdfs)
+        Array<epdf*> epdfs;
+        //! Array of indeces
+        Array<ivec> rvinds;
+public:
+        eprod ( const Array<epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
+                bool independent=true;
+                for ( int i=0;i<epdfs.length();i++ ) {
+                        independent=rv.add ( epdfs ( i )->_rv() );
+                        it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
+                        rvinds ( i ) = ( epdfs ( i )->_rv() ).dataind ( rv );
+                }
+        }
+
+        vec mean() const {
+                vec tmp ( rv.count() );
+                for ( int i=0;i<epdfs.length();i++ ) {
+                        vec pom = epdfs ( i )->mean();
+                        set_subvector ( tmp,rvinds ( i ), pom );
+                }
+                return tmp;
+        }
+        vec sample() const {
+                vec tmp ( rv.count() );
+                for ( int i=0;i<epdfs.length();i++ ) {
+                        vec pom = epdfs ( i )->sample();
+                        set_subvector ( tmp,rvinds ( i ), pom );
+                }
+                return tmp;
+        }
+        double evalpdflog ( const vec &val ) const {
+                double tmp=0;
+                for ( int i=0;i<epdfs.length();i++ ) {
+                        tmp+=epdfs(i)->evalpdflog(val(rvinds(i)));
+                }
+                return tmp;
+        }
+};
+
 
 /*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type

bdm/stat/libBM.cpp

@@ -47 +47 @@
 }
 
-RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ), names ( 0 ), sizes ( 0 ), times ( 0 ) {};
+RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ),  sizes ( 0 ), times ( 0 ),names ( 0 ) {};
 
 bool RV::add ( const RV &rv2 ) {

bdm/stat/libBM.h

@@ -265 +265 @@
         //!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.
+        //!  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 computational time.
         bool evalll;
 public:
@@ -288 +288 @@
         //!access function
         double _ll() const {return ll;}
+        //!access function 
+        void set_evalll(bool evl0){evalll=evl0;}
 };
 

bdm/stat/libEF.cpp

@@ -13 +13 @@
 
 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{
+        it_warning ( "Function not implemented" );
+        return vec_1 ( 0.0 );
+}
+
+double egiw::evalpdflog ( const vec &val ) const {
+
+        if ( xdim==1 ) { //same as the following, just quicker.
+                double r = val ( val.length()-1 ); //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();
+        }
+        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 ) );
+                mat Tmp=concat_vertical ( -eye ( xdim ),Th );
+                fsqmat iR ( xdim );
+                R.inv ( iR );
+                
+                return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T()*V.to_mat() *Tmp ) ) + 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);
+
+        double m = nu - nPsi -xdim-1;
+#define log2  0.693147180559945286226763983
+#define logpi 1.144729885849400163877476189
+#define log2pi 1.83787706640935
+
+        double nkG = 0.5* xdim* ( -nPsi *log2pi + sum ( log ( D ( xdim,D.length()-1 ) ) ) );
+        // temporary for lgamma in Wishart
+        double lg=0;
+        for ( int i =0; i<xdim;i++ ) {lg+=lgamma ( 0.5* ( m-i ) );}
+
+        double nkW = 0.5* ( m*sum ( log ( D ( 0,xdim-1 ) ) ) ) \
+                     - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
+
+        return nkG+nkW;
 }
 
@@ -41 +61 @@
         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 ) );
         
-        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);
-        
+        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;
 }
@@ -59 +88 @@
         for ( i=0; i<rv.count(); i++ ) {
                 GamRNG.setup ( alpha ( i ),beta ( i ) );
-                #pragma omp critical
+#pragma omp critical
                 smp ( i ) = GamRNG();
         }
@@ -69 +98 @@
 //      mat Smp ( rv.count(),N );
 //      int i,j;
-// 
+//
 //      for ( i=0; i<rv.count(); i++ ) {
 //              GamRNG.setup ( alpha ( i ),beta ( i ) );
-// 
+//
 //              for ( j=0; j<N; j++ ) {
 //                      Smp ( i,j ) = GamRNG();
 //              }
 //      }
-// 
+//
 //      return Smp;
 // }
@@ -121 +150 @@
 
         switch ( method ) {
-        case MULTINOMIAL:
-                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-                for ( i = n - 2;i >= 0;i-- ) {
-                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-                }
-
-                break;
-
-        case STRATIFIED:
-
-                for ( i = 0;i < n;i++ ) {
-                        u ( i ) = ( i + UniRNG.sample() ) / n;
-                }
-
-                break;
-
-        case SYSTEMATIC:
-                u0 = UniRNG.sample();
-
-                for ( i = 0;i < n;i++ ) {
-                        u ( i ) = ( i + u0 ) / n;
-                }
-
-                break;
-
-        default:
-                it_error ( "PF::resample(): Unknown resampling method" );
+                case MULTINOMIAL:
+                        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+
+                        for ( i = n - 2;i >= 0;i-- ) {
+                                u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+                        }
+
+                        break;
+
+                case STRATIFIED:
+
+                        for ( i = 0;i < n;i++ ) {
+                                u ( i ) = ( i + UniRNG.sample() ) / n;
+                        }
+
+                        break;
+
+                case SYSTEMATIC:
+                        u0 = UniRNG.sample();
+
+                        for ( i = 0;i < n;i++ ) {
+                                u ( i ) = ( i + u0 ) / n;
+                        }
+
+                        break;
+
+                default:
+                        it_error ( "PF::resample(): Unknown resampling method" );
         }
 
@@ -173 +202 @@
                         ind ( i ) = i;
                         N_babies ( i ) --; //this index was now replicated;
-                } else {
+                }
+                else {
                         // test if the parent has been fully replicated
                         // if yes, find the next one
@@ -198 +228 @@
 
 void eEmp::set_parameters ( const vec &w0, epdf* epdf0 ) {
-//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
+        //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         w=w0;
         w/=sum ( w0 );//renormalize

bdm/stat/libEF.h

@@ -113 +113 @@
 * \brief Gauss-inverse-Wishart density stored in LD form
 
-* More?...
+* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
+*
 */
 class egiw : public eEF {
@@ -121 +122 @@
         //! Number of data records (degrees of freedom) of sufficient statistics
         double nu;
-public:
-        //!Default constructor
+        //! Dimension of the output
+        int xdim;
+        //! Dimension of the regressor
+        int nPsi;
+public:
+        //!Default constructor, assuming
         egiw(RV rv, mat V0, double nu0): eEF(rv), V(V0), nu(nu0) {
-                it_assert_debug(rv.count()==V.rows(),"Incompatible V0.");
+                xdim = rv.count()/V.rows();
+                it_assert_debug(rv.count()==xdim*V.rows(),"Incompatible V0.");
+                nPsi = V.rows()-xdim;
         }
 
         vec sample() const;
         vec mean() 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 lognc () const;

bdm/stat/merger.h

@@ -53 +53 @@
         //! sample from merged density
         //! weight w is a 
-        vec sample(double &w, vec smp0=zeros ( rv.count()+rvc.count())){
+/*      vec sample(double &w, const vec &smp0){
                 // result
                 vec smp=smp0;
-                vec cond=sample(condpdf); // Just like in samplecond, here it is not given!
+                vec cond=condpdf.sample(); // Just like in samplecond, here it is not given!
                 // temporary
                 mat smpi=zeros(rv.count()+rvc.count(), n);
@@ -63 +63 @@
                 
                 // Consider arithmetic mean as a proposal density
-                ll(i) = 0;
+                ll = 0;
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
                         if ( rvcinds ( i ).length() > 0 ) {
@@ -99 +99 @@
                 
                 // copied from mprod.sample
-        };
+        };*/
 //      project 
         //! for future use

tests/merger_test.cpp

@@ -26 +26 @@
 
         merger M(A);
-        eEmp res=M.merge(100);
+//      eEmp res=M.merge(100);
         
-        cout << res.mean() << end;
+//      cout << res.mean() << endl;
         //Exit program:
         return 0;