Various EM code improvements by Charles Vaske and Andy Nguyen
authorJoris Mooij <j.mooij@cs.ru.nl>
Tue, 29 Apr 2014 14:13:35 +0000 (16:13 +0200)
committerJoris Mooij <j.mooij@cs.ru.nl>
Tue, 29 Apr 2014 14:13:35 +0000 (16:13 +0200)
ChangeLog
include/dai/emalg.h
src/emalg.cpp

index ac57ef9..5257ffc 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,6 @@
 git HEAD
 --------
+* [Charles Vaske] Various improvements to EM code, including interface changes
 * Fixed regression in SWIG interface caused by new build system (thanks to Benjamin Mumm for reporting this)
 * Added 'inclusive' argument to DAG::ancestors() and DAG::descendants()
 * [Julien Rebetez] added FactorGraph::printDot() to SWIG interface
index 543cc67..9f99c63 100644 (file)
@@ -73,15 +73,22 @@ class ParameterEstimation {
             (*_registry)[method] = f;
         }
 
-        /// Estimate the factor using the accumulated sufficient statistics and reset.
-        virtual Prob estimate() = 0;
+        /// Estimate the factor using the provided expectations.
+        virtual Prob estimate(const Prob &p) { return parametersToFactor(parameters(p)); }
 
-        /// Accumulate the sufficient statistics for \a p.
-        virtual void addSufficientStatistics( const Prob &p ) = 0;
+        /// Convert a set of estimated parameters to a factor
+        virtual Prob parametersToFactor(const Prob &params) = 0;
 
-        /// Returns the size of the Prob that should be passed to addSufficientStatistics.
+        /// Return parameters for the estimated factor, in a format specific to the parameter estimation
+        virtual Prob parameters(const Prob &p) = 0;
+
+        /// Returns the size of the Prob that should be passed to estimate and parameters
         virtual size_t probSize() const = 0;
 
+        // Returns the name of the class of this parameter estimation
+        virtual const std::string& name() const = 0;
+
+        virtual const PropertySet& properties() const = 0;
     private:
         /// A static registry containing all methods registered so far.
         static std::map<std::string, ParamEstFactory> *_registry;
@@ -98,11 +105,14 @@ class CondProbEstimation : private ParameterEstimation {
     private:
         /// Number of states of the variable of interest
         size_t _target_dim;
-        /// Current pseudocounts
-        Prob _stats;
         /// Initial pseudocounts
         Prob _initial_stats;
 
+        static std::string _name; // = "CondProbEstimation";
+
+        /// PropertySet that allows reconstruction of this estimator
+        PropertySet _props;
+
     public:
         /// Constructor
         /** For a conditional probability \f$ P( X | Y ) \f$,
@@ -132,13 +142,18 @@ class CondProbEstimation : private ParameterEstimation {
         /** The format of the resulting Prob keeps all the values for
          *  \f$ P(X | Y=y) \f$ in sequential order in the array.
          */
-        virtual Prob estimate();
-
-        /// Accumulate sufficient statistics from the expectations in \a p
-        virtual void addSufficientStatistics( const Prob &p );
-
-        /// Returns the required size for arguments to addSufficientStatistics().
-        virtual size_t probSize() const { return _stats.size(); }
+        virtual Prob parameters(const Prob &p);
+
+        // For a discrete conditional probability distribution,
+        // the parameters are equivalent to the resulting factor
+        virtual Prob parametersToFactor(const Prob &p);
+
+        /// Returns the required size for arguments to estimate().
+        virtual size_t probSize() const { return _initial_stats.size(); }
+        
+        virtual const std::string& name() const { return _name; }
+        
+        virtual const PropertySet& properties() const { return _props; }
 };
 
 
@@ -171,6 +186,8 @@ class SharedParameters {
         ParameterEstimation *_estimation;
         /// Indicates whether \c *this gets ownership of _estimation
         bool _ownEstimation;
+        /// The accumulated expectations
+        Prob* _expectations;
 
         /// Calculates the permutation that permutes the canonical ordering into the desired ordering
         /** \param varOrder Desired ordering of variables
@@ -197,10 +214,11 @@ class SharedParameters {
         SharedParameters( std::istream &is, const FactorGraph &fg );
 
         /// Copy constructor
-        SharedParameters( const SharedParameters &sp ) : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _ownEstimation(sp._ownEstimation) {
+        SharedParameters( const SharedParameters &sp ) : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _ownEstimation(sp._ownEstimation), _expectations(NULL) {
             // If sp owns its _estimation object, we should clone it instead of copying the pointer
             if( _ownEstimation )
                 _estimation = _estimation->clone();
+            _expectations = new Prob(*sp._expectations);
         }
 
         /// Destructor
@@ -208,24 +226,39 @@ class SharedParameters {
             // If we own the _estimation object, we should delete it now
             if( _ownEstimation )
                 delete _estimation;
+            if( _expectations != NULL) 
+                delete _expectations;
         }
 
-        /// Collect the sufficient statistics from expected values (beliefs) according to \a alg
+        /// Collect the expected values (beliefs) according to \a alg
         /** For each of the relevant factors (that shares the parameters we are interested in),
          *  the corresponding belief according to \a alg is obtained and its entries are permuted
          *  such that their ordering corresponds with the shared parameters that we are estimating.
-         *  Then, the parameter estimation subclass method addSufficientStatistics() is called with
-         *  this vector of expected values of the parameters as input.
          */
-        void collectSufficientStatistics( InfAlg &alg );
+        void collectExpectations( InfAlg &alg );
+
+        /// Return the current accumulated expectations
+        const Prob& currentExpectations() const { return *_expectations; }
+
+        ParameterEstimation& getPEst() const { return *_estimation; }
 
         /// Estimate and set the shared parameters
-        /** Based on the sufficient statistics collected so far, the shared parameters are estimated
+        /** Based on the expectation statistics collected so far, the shared parameters are estimated
          *  using the parameter estimation subclass method estimate(). Then, each of the relevant
          *  factors in \a fg (that shares the parameters we are interested in) is set according 
          *  to those parameters (permuting the parameters accordingly).
          */
         void setParameters( FactorGraph &fg );
+
+        /// Return a reference to the vector of factor orientations
+        /** This is necessary for determing which variables were used
+         *  to estimate parameters, and analysis of expectations
+         *  after an Estimation step has been performed.
+         */
+        const FactorOrientations& getFactorOrientations() const { return _varorders; }
+
+        /// Reset the current expectations
+        void clear( ) { _expectations->fill(0); }
 };
 
 
@@ -254,6 +287,9 @@ class MaximizationStep {
 
         /// Using all of the currently added expectations, make new factors with maximized parameters and set them in the FactorGraph.
         void maximize( FactorGraph &fg );
+        
+        /// Clear the step, to be called at the begining of each step
+        void clear( );
 
     /// \name Iterator interface
     //@{
index e8958c2..68969f8 100644 (file)
@@ -34,6 +34,8 @@ ParameterEstimation* ParameterEstimation::construct( const std::string &method,
 }
 
 
+std::string CondProbEstimation::_name = "CondProbEstimation";
+
 ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
     size_t target_dimension = p.getStringAs<size_t>("target_dim");
     size_t total_dimension = p.getStringAs<size_t>("total_dim");
@@ -43,37 +45,39 @@ ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
     return new CondProbEstimation( target_dimension, Prob( total_dimension, pseudo_count ) );
 }
 
-
 CondProbEstimation::CondProbEstimation( size_t target_dimension, const Prob &pseudocounts )
-  : _target_dim(target_dimension), _stats(pseudocounts), _initial_stats(pseudocounts)
+  : _target_dim(target_dimension), _initial_stats(pseudocounts), _props()
 {
-    DAI_ASSERT( !(_stats.size() % _target_dim) );
+    DAI_ASSERT( !(_initial_stats.size() % _target_dim) );
+    _props.setAsString<size_t>("target_dim", _target_dim);
+    _props.setAsString<size_t>("total_dim", _initial_stats.size());
+    _props.setAsString<Real>("pseudo_count", _initial_stats.get(0));
 }
 
 
-void CondProbEstimation::addSufficientStatistics( const Prob &p ) {
-    _stats += p;
+Prob CondProbEstimation::parametersToFactor(const Prob& p) {
+    Prob result(p);
+    return result;
 }
 
-
-Prob CondProbEstimation::estimate() {
+// In the case of a conditional probability table, the 
+// parameters are identical to the estimated factor
+Prob CondProbEstimation::parameters(const Prob& p) {
+    Prob stats = p + _initial_stats;
     // normalize pseudocounts
-    for( size_t parent = 0; parent < _stats.size(); parent += _target_dim ) {
+    for( size_t parent = 0; parent < stats.size(); parent += _target_dim ) {
         // calculate norm
         size_t top = parent + _target_dim;
         Real norm = 0.0;
         for( size_t i = parent; i < top; ++i )
-            norm += _stats[i];
+            norm += stats[i];
         if( norm != 0.0 )
             norm = 1.0 / norm;
         // normalize
         for( size_t i = parent; i < top; ++i )
-            _stats.set( i, _stats[i] * norm );
+            stats.set( i, stats[i] * norm );
     }
-    // reset _stats to _initial_stats
-    Prob result = _stats;
-    _stats = _initial_stats;
-    return result;
+    return stats;
 }
 
 
@@ -87,6 +91,7 @@ void SharedParameters::setPermsAndVarSetsFromVarOrders() {
     if( _varorders.size() == 0 )
         return;
     DAI_ASSERT( _estimation != NULL );
+    _expectations = new Prob(_estimation->probSize(), 0);
 
     // Construct the permutation objects and the varsets
     for( FactorOrientations::const_iterator foi = _varorders.begin(); foi != _varorders.end(); ++foi ) {
@@ -99,7 +104,7 @@ void SharedParameters::setPermsAndVarSetsFromVarOrders() {
 
 
 SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool ownPE )
-  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE)
+  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE), _expectations(NULL)
 {
     // Calculate the necessary permutations and varsets
     setPermsAndVarSetsFromVarOrders();
@@ -107,7 +112,7 @@ SharedParameters::SharedParameters( const FactorOrientations &varorders, Paramet
 
 
 SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
-  : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true)
+  : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true), _expectations(NULL)
 {
     // Read the desired parameter estimation method from the stream
     std::string est_method;
@@ -163,7 +168,7 @@ SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
 }
 
 
-void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
+void SharedParameters::collectExpectations( InfAlg &alg ) {
     for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
         Permute &perm = i->second;
         VarSet &vs = _varsets[i->first];
@@ -172,13 +177,13 @@ void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
         Prob p( b.nrStates(), 0.0 );
         for( size_t entry = 0; entry < b.nrStates(); ++entry )
             p.set( entry, b[perm.convertLinearIndex(entry)] ); // apply inverse permutation
-        _estimation->addSufficientStatistics( p );
+        (*_expectations) += p;
     }
 }
 
 
 void SharedParameters::setParameters( FactorGraph &fg ) {
-    Prob p = _estimation->estimate();
+    Prob p = _estimation->estimate(this->currentExpectations());
     for( std::map<FactorIndex, Permute>::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
         Permute &perm = i->second;
         VarSet &vs = _varsets[i->first];
@@ -203,7 +208,7 @@ MaximizationStep::MaximizationStep( std::istream &is, const FactorGraph &fg_varl
 
 void MaximizationStep::addExpectations( InfAlg &alg ) {
     for( size_t i = 0; i < _params.size(); ++i )
-        _params[i].collectSufficientStatistics( alg );
+        _params[i].collectExpectations( alg );
 }
 
 
@@ -213,6 +218,11 @@ void MaximizationStep::maximize( FactorGraph &fg ) {
 }
 
 
+void MaximizationStep::clear( ) {
+    for( size_t i = 0; i < _params.size(); ++i )
+        _params[i].clear( );
+}
+
 const std::string EMAlg::MAX_ITERS_KEY("max_iters");
 const std::string EMAlg::LOG_Z_TOL_KEY("log_z_tol");
 const size_t EMAlg::MAX_ITERS_DEFAULT = 30;
@@ -266,6 +276,8 @@ Real EMAlg::iterate( MaximizationStep &mstep ) {
     Real logZ = 0;
     Real likelihood = 0;
 
+    mstep.clear();
+
     _estep.run();
     logZ = _estep.logZ();