Changeset 213 for bdm/estim/merger.cpp

Timestamp:

11/25/08 12:24:42 (15 years ago)

Author:

smidl

Message:

Merging - new experiment

Files:

• bdm/estim/merger.cpp (modified) (6 diffs)

bdm/estim/merger.cpp

@@ -6 +6 @@
         int nu=lW.rows();
         vec mu = sum ( lW ) /nu; //mean of logs
-        vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 );
+        // vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
+        vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
         double coef=0.0;
         vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
@@ -61 +62 @@
         char str[100];
 
-        epdf* Mpred=Mix.predictor(rv);
+        epdf* Mpred=Mix.predictor ( rv );
         vec Mix_pdf ( Ns );
         while ( !converged ) {
@@ -70 +71 @@
                 delete Mpred;
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
-                
+
+                // This will be active only later in iterations!!!
                 if ( 1./sum_sqr ( w ) <0.5*Ns ) {
                         // Generate new samples
@@ -76 +78 @@
                         for ( int i=0;i<Ns;i++ ) {
                                 //////////// !!!!!!!!!!!!!
-                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = 0.01; }
+                                if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
                                 set_col_part ( Smp_ex,i,Smp ( i ) );
                         }
-                        {cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;}
+                        if(0){cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
                         cout << sum ( Smp_ex,2 ) /Ns <<endl;
-                        cout << Smp_ex*Smp_ex.T()/Ns << endl; 
+                        cout << Smp_ex*Smp_ex.T() /Ns << endl;}
                 }
 //              sprintf ( str,"Mpred_mean%d",niter );
@@ -172 +174 @@
                 }
                 else {
-
+                        it_file itf ( "merg_err.it" );
+                        itf << Name ( "w" ) << w;
                 }
 
@@ -188 +191 @@
                 converged = ( niter>20 );
         }
+        delete Mpred;
         cout << endl;