Changeset 189

Timestamp:

10/22/08 10:46:36 (16 years ago)

Author:

smidl

Message:

extend MixEF to allow for EM algorithm and alow estimation of weighted empirical density

Location:

bdm

Files:

• estim/mixef.cpp (modified) (7 diffs)
• estim/mixef.h (modified) (6 diffs)
• stat/emix.h (modified) (1 diff)

bdm/estim/mixef.cpp

@@ -21 +21 @@
         Coms ( 0 )->flatten ( Com0 );
 
-        //Copy it to the rest 
+        //Copy it to the rest
         for ( i=1;i<n;i++ ) {
                 //copy Com0 and create new rvs for them
@@ -38 +38 @@
 
 }
-void MixEF::bayesB ( const mat &Data ) {
-        this->bayes ( Data );
-}
 
-void MixEF::bayes ( const vec &data ) {
-
-};
-
-void MixEF::bayes ( const mat &data ) {
+void MixEF::bayesB ( const mat &data , const vec &wData ) {
         int ndat=data.cols();
         int t,i,niter;
@@ -63 +56 @@
         // tmp for weights
         vec wtmp = zeros ( n );
+        int maxi;
+        double maxll;
         //Estim
         while ( !converged ) {
@@ -75 +70 @@
                                 ll ( i ) += weights.logpred ( wtmp );
                         }
-                        w = exp ( ll-max ( ll ) );
-                        W.set_col ( t, w/sum ( w ) );
+                        
+                        maxll = max(ll,maxi);
+                        switch (method) {
+                                case QB:
+                                        w = exp ( ll-maxll );
+                                        w/=sum(w);
+                                        break;
+                                case EM:
+                                        w = 0.0;
+                                        w(maxi) = 1.0;
+                                        break;
+                        }
+                        
+                        W.set_col ( t, w );
                 }
 
+                // copy initial statistics
                 for ( i=0;i<n;i++ ) {
                         Coms ( i )-> set_statistics ( Coms0 ( i ) );
@@ -84 +92 @@
                 weights.set_statistics ( &weights0 );
 
+                // Update statistics
+                // !!!!    note  wData ==> this is extra weight of the data record
+                // !!!!    For typical cases wData=1.
                 for ( t=0;t<ndat;t++ ) {
                         for ( i=0;i<n;i++ ) {
-                                Coms ( i )-> bayes ( data.get_col ( t ),W ( i,t ) );
+                                Coms ( i )-> bayes ( data.get_col ( t ),W ( i,t ) * wData ( t ) );
                         }
-                        weights.bayes ( W.get_col ( t ) );
+                        weights.bayes ( W.get_col ( t ) * wData ( t ) );
                 }
 
@@ -98 +109 @@
         //Clean Coms0
         for ( i=0;i<n;i++ ) {delete Coms0 ( i );}
+}
+
+void MixEF::bayes ( const vec &data ) {
+
+};
+
+void MixEF::bayes ( const mat &data ) {
+        this->bayesB ( data, ones ( data.cols() ) );
 };
 
@@ -111 +130 @@
 }
 
-emix* MixEF::predictor(const RV &rv){
-        Array<epdf*> pC(n);
-        for(int i=0;i<n;i++){pC(i)=Coms(i)->predictor(rv);}
+emix* MixEF::predictor ( const RV &rv ) {
+        Array<epdf*> pC ( n );
+        for ( int i=0;i<n;i++ ) {pC ( i ) =Coms ( i )->predictor ( rv );}
         emix* tmp;
-        tmp = new emix(rv);
-        tmp->set_parameters(weights._epdf().mean(), pC, false);
+        tmp = new emix ( rv );
+        tmp->set_parameters ( weights._epdf().mean(), pC, false );
         tmp->ownComs();
         return tmp;

bdm/estim/mixef.h

@@ -20 +20 @@
 
 using namespace itpp;
+
+enum MixEF_METHOD { EM = 0, QB = 1};
 
 /*!
@@ -49 +51 @@
         ////!Indeces of component rvc in common rvc
         
+        //! Flag for a method that is used in the inference
+        MixEF_METHOD method;
+        
         //! Auxiliary function for use in constructors
         void build_est() {
@@ -65 +70 @@
         MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
                         BM ( RV() ), n ( Coms0.length() ), Coms ( n ),
-                        weights ( RV ( "{w }", vec_1 ( n ) ),alpha0 ) {
+                        weights ( RV ( "{w }", vec_1 ( n ) ),alpha0 ), method(QB) {
         //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
 
@@ -73 +78 @@
         MixEF () :
                         BM ( RV() ), n ( 0 ), Coms ( 0 ),
-                        weights ( RV ( "{w }", vec_1 ( 0 ) ),vec ( 0 ) ) {build_est();}
+                        weights ( RV ( "{w }", vec_1 ( 0 ) ),vec ( 0 ) ),method(QB) {build_est();}
         //! Initializing the mixture by a random pick of centroids from data
         //! \param Com0 Initial component - necessary to determine its type.
@@ -88 +93 @@
         //! EM algorithm
         void bayes ( const mat &dt );
-        void bayesB ( const mat &dt );
+        void bayesB ( const mat &dt, const vec &wData );
         double logpred ( const vec &dt ) const;
         const epdf& _epdf() const {return *est;}
@@ -94 +99 @@
         //! Flatten the density as if it was not estimated from the data
         void flatten(double sumw=1.0);
+        
+        //!Set which method is to be used
+        void set_method(MixEF_METHOD M){method=M;}
 };
 

bdm/stat/emix.h

@@ -112 +112 @@
                 return log(x);
         };
+        mat evalpdflog_M ( const mat &Val ) const {
+                mat X(w.length(), Val.cols());
+                for (int i = 0; i < w.length(); i++ ) {
+                        X.set_row(i, w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) ));
+                }
+                return X;
+        };
 
         emix* marginal ( const RV &rv ) const;