Changeset 193

Timestamp:

10/22/08 10:46:42 (16 years ago)

Author:

smidl

Message:

oprava merger_iter a jeho casti

Files:

• bdm/estim/merger.cpp (modified) (5 diffs)
• bdm/estim/mixef.h (modified) (1 diff)
• bdm/stat/emix.h (modified) (12 diffs)
• bdm/stat/libBM.h (modified) (1 diff)
• bdm/stat/libEF.cpp (modified) (1 diff)
• bdm/stat/libEF.h (modified) (1 diff)
• matlab/merger_iter_debug.m (modified) (3 diffs)
• tests/CMakeLists.txt (modified) (1 diff)
• tests/emix_test.cpp (added)
• tests/merger_iter_test.cpp (modified) (3 diffs)

bdm/estim/merger.cpp

@@ -60 +60 @@
                 //Re-Initialize Mixture model
                 Mix.init ( &A0, Smp_ex, Nc );
-                Mix.bayesB ( Smp_ex, ones ( Ns ) );//w*Ns );
+                Mix.bayesB ( Smp_ex, w*Ns );
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
 
-                // Generate new samples
-                eSmp.set_samples ( Mpred );
-                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+                if ( 1 ) {
+                        // Generate new samples
+                        eSmp.set_samples ( Mpred );
+                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+                }
+                else {
+                        for ( int ii=0;ii<10;ii++ ) {
+                                for ( int jj=0; jj<10; jj++ ) {
+                                        Smp ( ii+jj*10 ) =vec_2 ( -1.0+6*ii/10.0, -1.0+6*jj/10.0 );
+                                }
+                        }
+                        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+                }
 
                 sprintf ( str,"Mpdf%d",niter );
@@ -98 +108 @@
                                 // Compute likelihood of the missing variable
                                 if ( rv.count() > ( mpdfs ( i )->_rv().count() + mpdfs ( i )->_rvc().count() ) ) {
+                                        ///////////////
+                                        cout << Mpred->mean() <<endl;
                                         // There are variales unknown to mpdfs(i) : rvzs
                                         mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
                                         // Compute likelihood
+                                        vec lw_dbg=lw_src;
                                         for ( int k= 0; k<Ns; k++ ) {
                                                 lw_src ( k ) += log (
@@ -111 +124 @@
                                 // Compute likelihood of the partial source
                                 for ( int k= 0; k<Ns; k++ ) {
-                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp(k) ) );
-                                        lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( dls ( i )->get_val ( Smp(k) ) );
+                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
+                                        lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( dls ( i )->get_val ( Smp ( k ) ) );
                                 }
+                                
                         }
                         lW.set_row ( i, lw_src ); // do not divide by mix
@@ -139 +153 @@
                 dbg << Name ( str ) << w;
 
-                eSmp.resample(); // So that it can be used in bayes
-                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+//              eSmp.resample(); // So that it can be used in bayes
+//              for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
                 sprintf ( str,"Smp_res%d",niter );
@@ -147 +161 @@
                 // ==== stopping rule ===
                 niter++;
-                converged = ( niter>20 );
+                converged = ( niter>4 );
         }
 

bdm/estim/mixef.h

@@ -99 +99 @@
         //! Flatten the density as if it was not estimated from the data
         void flatten(double sumw=1.0);
+        //! Access function
+        BMEF* _Coms(int i){return Coms(i);}
         
         //!Set which method is to be used

bdm/stat/emix.h

@@ -42 +42 @@
         //!flag for destructor
         bool destroynom;
+        //!datalink between conditional and nom
+        datalink_m2e dl;
 public:
         //!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,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");
                         destroynom=true;
                 }
@@ -54 +57 @@
                         destroynom = false;
                 }
-                rvc = nom->_rv().subt ( rv );
                 it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
                 den = nom->marginal ( rvc );
         };
         double evalcond ( const vec &val, const vec &cond ) {
-                return exp ( nom->evalpdflog ( concat ( val,cond ) ) - den->evalpdflog ( 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 ) );
         }
         //! Object takes ownership of nom and will destroy it
@@ -105 +109 @@
         };
         vec evalpdflog_m ( const mat &Val ) const {
-                vec x=ones ( Val.cols() );
-                vec y ( Val.cols() );
+                vec x=zeros ( Val.cols() );
                 for ( int i = 0; i < w.length(); i++ ) {
-                        y = w ( i ) *exp ( Coms ( i )->evalpdflog_m ( Val ) );
-                        elem_mult_inplace ( y,x ); //result is in x
+                        x+= w ( i ) *exp ( Coms ( i )->evalpdflog_m ( Val ) );
                 }
                 return log ( x );
@@ -130 +132 @@
         //! Auxiliary function for taking ownership of the Coms()
         void ownComs() {destroyComs=true;}
+        
+        //!access function
+        epdf* _Coms(int i){return Coms(i);}
 };
 
@@ -153 +158 @@
                 // rv and rvc established = > we can link them with mpdfs
                 for ( int i = 0;i < n;i++ ) {
-                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs(i)->_rvc(), rv, rvc );
+                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv, rvc );
                 }
 
@@ -163 +168 @@
         double evalcond ( const vec &val, const vec &cond ) {
                 int i;
-                double res = 0.0;
+                double res = 1.0;
                 for ( i = n - 1;i >= 0;i-- ) {
-                        if ( mpdfs(i)->_rvc().count() >0) {
-                                mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
-                        }
-                        // add logarithms
-                        res += epdfs ( i )->evalpdflog ( dls ( i )->get_val ( val ) );
-                }
-                return exp ( res );
+                        /*                      if ( mpdfs(i)->_rvc().count() >0) {
+                                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
+                                                }
+                                                // 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 )
+                                                                                 );
+                }
+                return res;
         }
         vec samplecond ( const vec &cond, double &ll ) {
@@ -180 +189 @@
                 // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
-                        if ( mpdfs(i)->_rvc().count() ) {
+                        if ( mpdfs ( i )->_rvc().count() ) {
                                 mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
                         }
                         smpi = epdfs ( i )->sample();
                         // copy contribution of this pdf into smp
-                        dls(i)->fill_val(smp, smpi);
+                        dls ( i )->fill_val ( smp, smpi );
                         // add ith likelihood contribution
                         ll+=epdfs ( i )->evalpdflog ( smpi );
@@ -206 +215 @@
         Array<const epdf*> epdfs;
         //! Array of indeces
-        Array<ivec> rvinds;
-public:
-        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
+        Array<datalink_e2e*> dls;
+public:
+        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
                 bool independent=true;
                 for ( int i=0;i<epdfs.length();i++ ) {
                         independent=rv.add ( epdfs ( i )->_rv() );
                         it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
-                        rvinds ( i ) = ( epdfs ( i )->_rv() ).dataind ( rv );
+                        dls ( i ) = new datalink_e2e( epdfs ( i )->_rv() , rv );
                 }
         }
@@ -221 +230 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        set_subvector ( tmp,rvinds ( i ), pom );
+                        dls(i)->fill_val(tmp, pom );
                 }
                 return tmp;
@@ -229 +238 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->sample();
-                        set_subvector ( tmp,rvinds ( i ), pom );
+                        dls(i)->fill_val(tmp, pom );
                 }
                 return tmp;
@@ -236 +245 @@
                 double tmp=0;
                 for ( int i=0;i<epdfs.length();i++ ) {
-                        tmp+=epdfs ( i )->evalpdflog ( val ( rvinds ( i ) ) );
+                        tmp+=epdfs ( i )->evalpdflog ( dls(i)->get_val ( val ) );
                 }
                 return tmp;
@@ -242 +251 @@
         //!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);}}
 };
 

bdm/stat/libBM.h

@@ -296 +296 @@
                 return tmp;
         }
+        void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
+                it_assert_debug ( valsize==val.length(),"Wrong val" );
+                it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
+                set_subvector ( val_up, v2v_up, val );
+                set_subvector ( val_up, v2c_up, cond );
+        }
 };
 //!DataLink is a connection between mpdf and its superordinate (Up)

bdm/stat/libEF.cpp

@@ -82 +82 @@
                 mat R;
                 mean_mat(M,R);
-                return cvectorize (concat_vertical(M,R));
+                return cvectorize (concat_vertical(M.T(),R));
         }
 

bdm/stat/libEF.h

@@ -620 +620 @@
 
         //fixme - could this be done in general for all sq_T?
-        mat S=R.to_mat();
+        mat S=Rn.to_mat();
         //fixme
         int n=rvn.count()-1;

matlab/merger_iter_debug.m

@@ -7 +7 @@
 
 figure(1);
-for it=0:20
+for it=0:4
     figure(it+1)
     subplot(2,2,1);
@@ -24 +24 @@
     subplot(2,2,3);
     eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(2,:)))']);
+    %eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lw_cond' si '))']);
     title('Second source');
     set(gca,'XLim',XL);
@@ -30 +31 @@
     
     subplot(2,2,4);
-    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ')']);
+    hold off
+    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ','':'')']);
     title('Merged density'); 
     set(gca,'XLim',XL);
     set(gca,'YLim',YL);
+    hold on
+    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w_is_' si ')']);
     
 end

tests/CMakeLists.txt

@@ -34 +34 @@
 TEST(enorm_test)
 TEST(egiw_test)
+TEST(emix_test)
 TEST(test0)
 TEST(testResample)

tests/merger_iter_test.cpp

@@ -11 +11 @@
 int main() {
 
-        RNG_randomize();
+        //RNG_randomize();
         
         RV x ( "{x }","1" );
@@ -21 +21 @@
         enorm<fsqmat> f2 ( x );
                 
-        f1.set_parameters ( "0 0",mat ( "0.5 0; 0 0.5" ) );
+        f1.set_parameters ( "3 2",mat ( "0.5 0; 0 0.5" ) );
         f2.set_parameters ( "3",mat ( "0.5" ) );
         
@@ -44 +44 @@
         g0.set_parameters(vec("1 1"),mat("1 0; 0 1"));
         
-        M.set_parameters(1.2,200,2);
+        M.set_parameters(1.2,100,1);
         M.merge(&g0);