Changeset 189 for bdm/estim/mixef.cpp

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

Files:

• bdm/estim/mixef.cpp (modified) (7 diffs)

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;