Changeset 204 for bdm

Timestamp:

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

Author:

smidl

Message:

merger is now in logarithms + new merge_test

Location:

bdm

Files:

• estim/arx.cpp (modified) (1 diff)
• estim/merger.cpp (modified) (7 diffs)
• stat/emix.h (modified) (9 diffs)
• stat/libEF.h (modified) (1 diff)

bdm/estim/arx.cpp

@@ -15 +15 @@
         nu+=w;
 
+        // log(sqrt(2*pi)) = 0.91893853320467
         if ( evalll ) {
                 lnc = est.lognc();
-                ll = lnc - last_lognc;
+                ll = lnc - last_lognc - 0.91893853320467;
                 last_lognc = lnc;
         }

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;
 
 }

bdm/stat/emix.h

@@ -47 +47 @@
         //!Default constructor. By default, the given epdf is not copied!
         //! It is assumed that this function will be used only temporarily.
-        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( rv,nom0->_rv().subt(rv) ), dl(rv,rvc,nom0->_rv()) {
+        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( rv,nom0->_rv().subt ( rv ) ), dl ( rv,rvc,nom0->_rv() ) {
                 if ( copy ) {
 //                      nom = nom0->_copy_();
-                        it_error("todo");
+                        it_error ( "todo" );
                         destroynom=true;
                 }
@@ -61 +61 @@
         };
         double evalcond ( const vec &val, const vec &cond ) {
-                vec nom_val(rv.count()+rvc.count());
-                dl.fill_val_cond(nom_val,val,cond);
+                vec nom_val ( rv.count() +rvc.count() );
+                dl.fill_val_cond ( nom_val,val,cond );
                 return exp ( nom->evalpdflog ( nom_val ) - den->evalpdflog ( cond ) );
         }
@@ -132 +132 @@
         //! Auxiliary function for taking ownership of the Coms()
         void ownComs() {destroyComs=true;}
-        
+
         //!access function
-        epdf* _Coms(int i){return Coms(i);}
+        epdf* _Coms ( int i ) {return Coms ( i );}
 };
 
@@ -175 +175 @@
                                                 // add logarithms
                                                 res += epdfs ( i )->evalpdflog ( dls ( i )->get_val ( val ) );*/
-                        res *= mpdfs ( i )->evalcond ( 
-                                                   dls ( i )->get_val ( val ),
-                        dls ( i )->get_cond ( val, cond )
-                                                                                 );
+                        res *= mpdfs ( i )->evalcond (
+                                   dls ( i )->get_val ( val ),
+                                   dls ( i )->get_cond ( val, cond )
+                               );
                 }
                 return res;
@@ -222 +222 @@
                         independent=rv.add ( epdfs ( i )->_rv() );
                         it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
-                        dls ( i ) = new datalink_e2e( epdfs ( i )->_rv() , rv );
+                }
+                for ( int i=0;i<epdfs.length();i++ ) {
+                        dls ( i ) = new datalink_e2e ( epdfs ( i )->_rv() , rv );
                 }
         }
@@ -230 +232 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls(i)->fill_val(tmp, pom );
+                        dls ( i )->fill_val ( tmp, pom );
                 }
                 return tmp;
@@ -238 +240 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->sample();
-                        dls(i)->fill_val(tmp, pom );
+                        dls ( i )->fill_val ( tmp, pom );
                 }
                 return tmp;
@@ -245 +247 @@
                 double tmp=0;
                 for ( int i=0;i<epdfs.length();i++ ) {
-                        tmp+=epdfs ( i )->evalpdflog ( dls(i)->get_val ( val ) );
+                        tmp+=epdfs ( i )->evalpdflog ( dls ( i )->get_val ( val ) );
                 }
                 return tmp;
@@ -251 +253 @@
         //!access function
         const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}
-        
+
         //!Destructor
-        ~eprod(){for(int i=0;i<epdfs.length();i++){delete dls(i);}}
+        ~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
 };
 

bdm/stat/libEF.h

@@ -48 +48 @@
         virtual double evalpdflog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
         //!Evaluate normalized log-probability
-        virtual double evalpdflog ( const vec &val ) const {return evalpdflog_nn ( val )-lognc();}
+        virtual double evalpdflog ( const vec &val ) const {double tmp;tmp= evalpdflog_nn ( val )-lognc();it_assert_debug(std::isfinite(tmp),"why?"); return tmp;}
         //!Evaluate normalized log-probability for many samples
         virtual vec evalpdflog ( const mat &Val ) const {