Changeset 1009 for library/bdm/estim/mixtures.cpp

Timestamp:

05/27/10 23:07:16 (14 years ago)

Author:

smidl

Message:

changes in bayes_batch

Files:

• library/bdm/estim/mixtures.cpp (modified) (7 diffs)

library/bdm/estim/mixtures.cpp

@@ -28 +28 @@
         for ( i = 0; i < Coms.length(); i++ ) {
                 //pick one datum
-                int ind = (int) floor ( ndat * UniRNG.sample() );
-                Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
+                if (ndat==Coms.length()){ //take the ith vector
+                        Coms ( i )->bayes ( Data.get_col ( i ), empty_vec );
+                } else { // pick at random
+                        int ind = (int) floor ( ndat * UniRNG.sample() );
+                        Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
+                }
                 //flatten back to oringinal
                 Coms ( i )->flatten ( Com0 );
@@ -36 +40 @@
 }
 
-void MixEF::bayes_batch ( const mat &data , const mat &cond, const vec &wData ) {
+double MixEF::bayes_batch ( const mat &data , const mat &cond, const vec &wData ) {
         int ndat = data.cols();
         int t, i, niter;
@@ -58 +62 @@
         int maxi;
         double maxll;
+        
+        double levid=0.0;
         //Estim
         while ( !converged ) {
+                levid=0.0;
                 // Copy components back to their initial values
                 // All necessary information is now in w and Coms0.
@@ -68 +75 @@
                         //#pragma omp parallel for
                         for ( i = 0; i < n; i++ ) {
-                                ll ( i ) = Coms ( i )->logpred ( data.get_col ( t ) );
+                                ll ( i ) = Coms ( i )->logpred ( data.get_col ( t ) , empty_vec);
                                 wtmp = 0.0;
                                 wtmp ( i ) = 1.0;
@@ -99 +106 @@
                 // !!!!    note  wData ==> this is extra weight of the data record
                 // !!!!    For typical cases wData=1.
+                vec logevid(n);
                 for ( t = 0; t < ndat; t++ ) {
                         //#pragma omp parallel for
                         for ( i = 0; i < n; i++ ) {
                                 Coms ( i )-> bayes_weighted ( data.get_col ( t ), empty_vec, W ( i, t ) * wData ( t ) );
+                                logevid(i) = Coms(i)->_ll();
                         }
                         weights.bayes ( W.get_col ( t ) * wData ( t ) );
                 }
-
+                levid += weights._ll()+log(weights.posterior().mean() * exp(logevid)); // inner product w*exp(evid)
+                
                 niter++;
                 //TODO better convergence rule.
@@ -116 +126 @@
                 delete Coms0 ( i );
         }
+        return levid;
 }
 
@@ -127 +138 @@
 
 
-double MixEF::logpred ( const vec &yt ) const {
+double MixEF::logpred ( const vec &yt, const vec &cond =empty_vec) const {
 
         vec w = weights.posterior().mean();
         double exLL = 0.0;
         for ( int i = 0; i < Coms.length(); i++ ) {
-                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt ) );
+                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt , cond ) );
         }
         return log ( exLL );