Changeset 205 for bdm/estim/merger.cpp

Timestamp:

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

Author:

smidl

Message:

merger posledni verze

Files:

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

bdm/estim/merger.cpp

@@ -16 +16 @@
                 case 3://Ratio of Bessel
                         coef = sqrt ( ( 3*beta-2 ) /3*beta );
-                        return log( besselk ( 0,sq2bl*coef )) - log( besselk ( 0,sq2bl ) ) +  mu;
+                        return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
                         break;
                 case 4:
@@ -31 +31 @@
         it_assert_debug ( rv.equal ( g0->_rv() ),"Incompatible g0" );
         //Empirical density - samples
-        eEmp eSmp ( rv,Ns );
         eSmp.set_parameters ( ones ( Ns ), g0 );
         Array<vec> &Smp = eSmp._samples(); //aux
@@ -62 +61 @@
         char str[100];
 
-        epdf* Mpred;
+        epdf* Mpred=Mix.predictor(rv);
         vec Mix_pdf ( Ns );
         while ( !converged ) {
@@ -69 +68 @@
                 Mix.flatten ( &Mix_init );
                 Mix.bayesB ( Smp_ex, w*Ns );
+                delete Mpred;
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
-
-//              sprintf ( str,"Mpred_mean%d",niter );
-//              dbg << Name ( str ) << Mpred->mean();
-
-                if ( 1 ) {
+                
+                if ( 1./sum_sqr ( w ) <0.5*Ns ) {
                         // Generate new samples
                         eSmp.set_samples ( Mpred );
                         for ( int i=0;i<Ns;i++ ) {
                                 //////////// !!!!!!!!!!!!!
-                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = GamRNG(); }
+                                if ( Smp ( i ) ( 1 ) <0 ) {Smp ( i ) ( 1 ) = 0.01; }
                                 set_col_part ( Smp_ex,i,Smp ( i ) );
                         }
-                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+                        {cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;}
                         cout << sum ( Smp_ex,2 ) /Ns <<endl;
+                        cout << Smp_ex*Smp_ex.T()/Ns << endl; 
                 }
-                else
-                        {cout<<"Eff. sample size=" << 1./sum_sqr ( w ) << endl;}
+//              sprintf ( str,"Mpred_mean%d",niter );
+//              dbg << Name ( str ) << Mpred->mean();
+
 
 //              sprintf ( str,"Mpdf%d",niter );
@@ -135 +134 @@
                                                                         zdls ( i )->get_cond ( Smp ( k ) ) ) );
                                                 if ( !std::isfinite ( lw_src ( k ) ) ) {
-                                                        cout << endl;
+                                                        lw_src ( k ) = -1e16; cout << "!";
                                                 }
                                         }
@@ -166 +165 @@
                 //Importance weighting
                 lw -=  lw_mix; // hoping that it is not numerically sensitive...
-                w = exp(lw-max(lw));
+                w = exp ( lw-max ( lw ) );
                 //renormalize
                 double sumw =sum ( w );
@@ -173 +172 @@
                 }
                 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;