Changeset 204 for bdm/estim/merger.cpp

Timestamp:

11/10/08 15:40:29 (16 years ago)

Author:

smidl

Message:

merger is now in logarithms + new merge_test

Files:

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

bdm/estim/merger.cpp

@@ -11 +11 @@
         switch ( nu ) {
                 case 2:
-                        coef= ( 1-0.5*sqrt ( ( 4*beta-3 ) /beta ) );
-                        return exp ( coef*sq2bl + mu );
+                        coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
+                        return  coef*sq2bl + mu ;
                         break;
                 case 3://Ratio of Bessel
                         coef = sqrt ( ( 3*beta-2 ) /3*beta );
-                        return elem_mult (
-                                   elem_div ( besselk ( 0,sq2bl*coef ), besselk ( 0,sq2bl ) ),
-                                   exp ( mu ) );
+                        return log( besselk ( 0,sq2bl*coef )) - log( besselk ( 0,sq2bl ) ) +  mu;
                         break;
                 case 4:
@@ -46 +44 @@
         vec lw_src ( Ns );
         vec lw_mix ( Ns );
+        vec lw ( Ns );
         mat lW=zeros ( n,Ns );
         vec vec0 ( 0 );
@@ -75 +74 @@
 //              dbg << Name ( str ) << Mpred->mean();
 
-                if ( niter<10 ) {
+                if ( 1 ) {
                         // Generate new samples
                         eSmp.set_samples ( Mpred );
-                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
-
-                }
+                        for ( int i=0;i<Ns;i++ ) {
+                                //////////// !!!!!!!!!!!!!
+                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = GamRNG(); }
+                                set_col_part ( Smp_ex,i,Smp ( i ) );
+                        }
+                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+                        cout << sum ( Smp_ex,2 ) /Ns <<endl;
+                }
+                else
+                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+
 //              sprintf ( str,"Mpdf%d",niter );
 //              for ( int i=0;i<Ns;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
@@ -127 +134 @@
                                                                         zdls ( i )->get_val ( Smp ( k ) ),
                                                                         zdls ( i )->get_cond ( Smp ( k ) ) ) );
+                                                if ( !std::isfinite ( lw_src ( k ) ) ) {
+                                                        cout << endl;
+                                                }
                                         }
                                         delete tmp_cond;
@@ -147 +157 @@
 //              dbg << Name ( str ) << lW;
 
-                w = lognorm_merge ( lW ); //merge
+                lw = lognorm_merge ( lW ); //merge
 
 //              sprintf ( str,"w%d",niter );
@@ -155 +165 @@
 
                 //Importance weighting
-                w /=exp ( lw_mix ); // hoping that it is not numerically sensitive...
+                lw -=  lw_mix; // hoping that it is not numerically sensitive...
+                w = exp(lw-max(lw));
                 //renormalize
-                w /=sum ( w );
-
+                double sumw =sum ( w );
+                if ( std::isfinite ( sumw ) ) {
+                        w = w/sumw;
+                }
+                else {
+                        
+                }
+                {
+                        double eff = 1./sum_sqr ( w );
+                        if ( eff<2 ) {
+                                int mi= max_index ( w );
+                                cout << w(mi) <<endl;
+                                cout << lW.get_col ( mi ) <<endl;
+                                mat mm(2,1);mm.set_col(0,lW.get_col(mi));
+                                cout << lognorm_merge(mm) <<endl;
+                                cout << lw_mix ( mi ) <<endl;
+                                cout << lw ( mi ) <<endl;
+                                cout << Smp ( mi ) <<endl;
+                                cout << Mix._Coms(0)->_e()->mean() <<endl;
+                        }
+                        cout << "Eff: " << eff <<endl;
+                }
 //              sprintf ( str,"w_is_%d",niter );
 //              dbg << Name ( str ) << w;
@@ -170 +201 @@
                 // ==== stopping rule ===
                 niter++;
-                converged = ( niter>5 );
+                converged = ( niter>20 );
         }
+        cout << endl;
 
 }