Changeset 165

Timestamp:

09/16/08 10:46:21 (16 years ago)

Author:

smidl

Message:

Switch from eprod to mprod, merger devel

Files:

• bdm/stat/emix.cpp (modified) (1 diff)
• bdm/stat/emix.h (modified) (3 diffs)
• bdm/stat/libBM.cpp (modified) (2 diffs)
• bdm/stat/libBM.h (modified) (1 diff)
• bdm/stat/merger.h (modified) (1 diff)
• tests/merger_test.cpp (modified) (1 diff)

bdm/stat/emix.cpp

@@ -22 +22 @@
         return Coms ( i )->sample();
 }
+
+mprod::mprod ( Array<mpdf*> mFacs, bool overlap) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), rvcinds ( n ) {
+                int i;
+                bool rvaddok;
+                // Create rv
+                for ( i = 0;i < n;i++ ) {
+                        rvaddok=rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
+                        // If rvaddok==false, mpdfs overlap => assert error.
+                        it_assert_debug(rvaddok||overlap,"mprod::mprod() input mpdfs overlap in rv!");
+                        epdfs ( i ) = & ( mpdfs ( i )->_epdf() ); // add pointer to epdf
+                };
+                // Create rvc
+                for ( i = 0;i < n;i++ ) {
+                        rvc.add ( mpdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
+                };
+                
+                independent = true;
+                //test rvc of mpdfs and fill rvinds
+                for ( i = 0;i < n;i++ ) {
+                        // find ith rv in common rv
+                        rvinds ( i ) = mpdfs ( i )->_rv().dataind ( rv );
+                        // find ith rvc in common rv
+                        rvcinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
+                        // find ith rvc in common rv
+                        rvcinds ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
+                        //
+                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
+                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
+                }
+        };

bdm/stat/emix.h

@@ -87 +87 @@
 //      bool intermpdfs;
 public:
-        //!Constructor from list of eFacs or list of mFacs
-        mprod ( Array<mpdf*> mFacs ) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), rvcinds ( n ) {
-                int i;
-                // Create rv
-                for ( i = 0;i < n;i++ ) {
-                        rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
-                        epdfs ( i ) = & ( mpdfs ( i )->_epdf() ); // add pointer to epdf
-                };
-                // Create rvc
-                for ( i = 0;i < n;i++ ) {
-                        rvc.add ( mpdfs ( i )->_rv().subt ( rv ) ); //add rv to common rvs.
-                };
-
-                independent = true;
-                //test rvc of mpdfs and fill rvinds
-                for ( i = 0;i < n;i++ ) {
-                        // find ith rv in common rv
-                        rvinds ( i ) = mpdfs ( i )->_rv().dataind ( rv );
-                        // find ith rvc in common rv
-                        rvcinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
-                        // find ith rvc in common rv
-                        rvcinds ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
-                        //
-                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
-                        if ( rvcinds ( i ).length() >0 ) {independent = false;}
-                }
-        };
+        /*!\brief Constructor from list of mFacs, 
+        Additional parameter overlap is left for future use. Do not set to true for mprod.
+        */
+        mprod ( Array<mpdf*> mFacs, bool overlap=false ); 
 
         double evalpdflog ( const vec &val ) const {
@@ -129 +106 @@
                 vec smp=zeros ( rv.count() );
                 vec condi;
+                vec smpi;
+                ll = 0;
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
                         if ( rvcinds ( i ).length() > 0 ) {
@@ -139 +118 @@
                                 mpdfs ( i )->condition ( condi );
                         }
+                        smpi = epdfs ( i )->sample();
                         // copy contribution of this pdf into smp
-                        set_subvector ( smp,rvinds ( i ), epdfs ( i )->sample() );
+                        set_subvector ( smp,rvinds ( i ), smpi );
+                        // add ith likelihood contribution
+                        ll+=epdfs(i)->evalpdflog(smpi);
                 }
                 return smp;

bdm/stat/libBM.cpp

@@ -70 +70 @@
                 tsize = sum ( sizes );
                 len = ids.length();
-                return ( index.length() <rv2.len );
+                return ( index.length() ==rv2.len );
         }
         else { //rv2 is empty
@@ -141 +141 @@
 
 ivec RV::dataind ( RV rv2 ) const {
-        str str2 = rv2.tostr();
         ivec res ( 0 );
-        ivec part;
-        int i;
-        for ( i = 0;i < len;i++ ) {
-                part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
-                res = concat ( res, part );
+        if ( rv2.count()>0 ) {
+                str str2 = rv2.tostr();
+                ivec part;
+                int i;
+                for ( i = 0;i < len;i++ ) {
+                        part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
+                        res = concat ( res, part );
+                }
         }
         return res;
+
 }
 
 RV RV::subt ( const RV rv2 ) const {
-        cout << *this <<endl;
-        cout << rv2 <<endl;
-
         ivec res = this->findself ( rv2 ); // nonzeros
         ivec valid = itpp::find ( res == -1 ); //-1 => value not found => it remains

bdm/stat/libBM.h

@@ -77 +77 @@
         //! Compare if \c rv2 is identical to this \c RV
         bool equal (const RV &rv2 ) const;
-        //! Add (concat) another variable to the current one, \return 0 if all rv2 were added, 1 if rv2 is in conflict
+        //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
         bool add ( const RV &rv2 );
         //! Subtract  another variable from the current one

bdm/stat/merger.h

@@ -28 +28 @@
 */
 
-class merger {
+class merger : public mprod{
 protected:      
-        mat Smp;
-        mprod joint;
+        //! Additional pdf on the part in condition (if undefined);
+        enorm<fsqmat> condpdf;
+        //! Find potential overlaps in rv
+        Array<ivec> overlaps;
 public:
 //!Default constructor
-        merger (const Array<mpdf*> &in_sources) : joint(in_sources) {
-                RV rvc=joint._rvc();
+        merger (const Array<mpdf*> &in_sources) : mprod(in_sources,true), condpdf(rvc), overlaps(rv.count()) {
+                if (rvc.count()>0) {
+                        vec mu = zeros(rvc.count());
+                        mat R = 100*eye(rvc.count());
+                        condpdf.set_parameters(mu,R);
+                }
+                for (int i=0;i<rv.count();i++){ // cycle over rv
+                        ivec tmp(0);
+                        for (int j=0;j<n;j++){//cycle over rvinds
+                                ivec id=itpp::find(rvinds(j)==i);
+                                tmp=concat(tmp,id); // add i to tmp if found
+                        }
+                        overlaps(i)=tmp;
+                }
         }
+        //! sample from merged density
+        //! weight w is a 
+        vec sample(double &w, vec smp0=zeros ( rv.count()+rvc.count())){
+                // result
+                vec smp=smp0;
+                vec cond=sample(condpdf); // Just like in samplecond, here it is not given!
+                // temporary
+                mat smpi=zeros(rv.count()+rvc.count(), n);
+                vec condi;
+                vec ll(n);
+                
+                // Consider arithmetic mean as a proposal density
+                ll(i) = 0;
+                for ( int i = ( n - 1 );i >= 0;i-- ) {
+                        if ( rvcinds ( i ).length() > 0 ) {
+                                condi = zeros ( rvcinrv.length() + rvcinds.length() );
+                                // copy data in condition
+                                set_subvector ( condi,rvcinds ( i ), cond );
+                                // copy data in already generated sample
+                                set_subvector ( condi,rvcinrv ( i ), smp );
 
+                                mpdfs ( i )->condition ( condi );
+                        }
+                        smpi.set_col(i, epdfs ( i )->sample());
+                        // add ith likelihood contribution
+                        ll+=epdfs(i)->evalpdflog(smpi);
+                }
+                // Now lets analyze samples smpi
+                for (int i=0;i<rv.count();i++){
+                        // number of components
+                        int noc = overlaps(i).length();
+                        if (noc<2) { // only one sample in this dimension, take it
+                                smp(i) = smpi(i,overlaps(i)(0));
+                        }
+                        else {
+                                //pick dimension
+                                double unis = UniRNG.sample();
+                                dim = 0;
+                                while(unis<(dim+1)/noc) {dim++;} // 
+                                //in dim we now have randomly picked dimension
+                                smp(i) = smpi(i,overlaps(i)(dim));
+                        }
+                }
+                
+                // now we have sample how to evaluate weight?
+                return smp;
+                
+                // copied from mprod.sample
+        };
 //      project 
         //! for future use

tests/merger_test.cpp

@@ -26 +26 @@
 
         merger M(A);
-
+        eEmp res=M.merge(100);
+        
+        cout << res.mean() << end;
         //Exit program:
         return 0;