Changeset 204

Timestamp:

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

Author:

smidl

Message:

merger is now in logarithms + new merge_test

Files:

• bdm/estim/arx.cpp (modified) (1 diff)
• bdm/estim/merger.cpp (modified) (7 diffs)
• bdm/stat/emix.h (modified) (9 diffs)
• bdm/stat/libEF.h (modified) (1 diff)
• mpdm/CMakeLists.txt (modified) (1 diff)
• mpdm/merg_giw.cpp (added)
• mpdm/merg_pred.cpp (modified) (7 diffs)

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 {

mpdm/CMakeLists.txt

@@ -20 +20 @@
 target_link_libraries (merg_pred ${BdmLibs})
 
+add_executable (merg_giw merg_giw.cpp)
+target_link_libraries (merg_giw ${BdmLibs})
+
 add_executable (merger_iter_cond merger_iter_cond.cpp)
 target_link_libraries (merger_iter_cond ${BdmLibs})

mpdm/merg_pred.cpp

@@ -34 +34 @@
         mat V0g = 0.001*eye ( thrg.count() ); V0g ( 0,0 ) *= 10; //
         double nu0 = ord+6.0;
-        double frg = 1.0;//0.99;
+        double frg = 0.95;
 
         ARX P1 ( thri, V0, nu0, frg );
@@ -40 +40 @@
         ARX PG ( thrg, V0g, nu0, frg );
         //Test estimation
-        int ndat = 500;
+        int ndat = 200;
         int t;
 
@@ -51 +51 @@
         int Li_Data = L.add ( RV ( "{Y U1 U2 }" ), "" );
         int Li_LL   = L.add ( RV ( "{1 2 G }" ), "LL" );
-        int Li_Pred = L.add ( RV ( "{1 2 G ar ge ln}" ), "Pred" );
+        int Li_Pred = L.add ( RV ( "{1 2 G ar ge }" ), "Pred" );
 
         L.init();
@@ -66 +66 @@
         vec PostLLs(3);
         vec PredLLs_m=zeros(5);
-        vec PostLLs_m=zeros(3);
         ivec ind_r1 = "0 1 3";
         ivec ind_r2 = "0 2 3";
@@ -84 +83 @@
                         // Test predictors
                         if (t>2){
-                                mlnorm<ldmat>* P1p = P1.predictor_student(y,concat(ym,u1));
-                                mlnorm<ldmat>* P2p = P2.predictor_student(y,concat(ym,u2));
-                                mlnorm<ldmat>* Pgp = PG.predictor_student(y,concat(ym,uu));
+                                mlnorm<ldmat>* P1p = P1.predictor(y,concat(ym,u1));
+                                mlnorm<ldmat>* P2p = P2.predictor(y,concat(ym,u2));
+                                mlnorm<ldmat>* Pgp = PG.predictor(y,concat(ym,uu));
                                 
                                 Array<mpdf*> A(2); A(0)=P1p;A(1)=P2p;
                                 merger M(A);
                                 enorm<ldmat> g0(concat(yy,uu)); g0.set_parameters("0 0 0 0 ",3*eye(4));
-                                M.set_parameters(1.2, 100,1);
+                                M.set_parameters(10000000.0, 100,1);
                                 M.merge(&g0);
                                 
@@ -98 +97 @@
                                 double P2pl = P2p->evalcond(yt,rgr2);   
                                 double PGpl = Pgp->evalcond(yt,rgrg);   
-                                PredLLs(0) = log(P1pl); PredLLs(1) = log(P2pl); PredLLs(2) = log(PGpl);
+                                {
+                                        cout << "yt: " << yt << endl;
+                                        cout << "yt_1: " << P1p->_epdf().mean() << endl;
+                                        cout << "yt_2: " << P2p->_epdf().mean() << endl;
+                                        cout << "yt_G: " << P2p->_epdf().mean() << endl;
+                                }
+                                double cP1pl;
+                                double cP2pl;
                                 {
                                         ARX* Apred = (ARX*)M._Mix()._Coms(0);
@@ -107 +113 @@
                                         mlnorm<ldmat>* cP2p = (mlnorm<ldmat>*)mP2p->condition(y);
 
-                                        double cP1pl = cP1p->evalcond(yt,rgr1); 
-                                        double cP2pl = cP2p->evalcond(yt,rgr2); 
-                                        PredLLs(3) = log(cP1pl);
-                                        PredLLs(4) = log(cP2pl);
+                                        cP1pl = cP1p->evalcond(yt,rgr1);        
+                                        cP2pl = cP2p->evalcond(yt,rgr2);        
+                                
+                                        cout << "ytm1: " << cP1p->_epdf().mean() << endl;
+                                        cout << "ytm2: " << cP2p->_epdf().mean() << endl;
                                 }
 
-                                L.logit(Li_Pred, PredLLs_m*frg+ PredLLs); //log-normal
-                                PredLLs_m *=frg;
-                                PredLLs_m += PredLLs;
+                                PredLLs *=frg;
+                                PredLLs += log(concat(vec_3(P1pl, P2pl, PGpl), vec_2(cP1pl,cP2pl)));
+                                L.logit(Li_Pred, PredLLs); //log-normal
                                 
                                 delete P1p;