Changeset 180 for bdm/estim

Timestamp:

10/15/08 19:08:06 (16 years ago)

Author:

smidl

Message:

Modifications of BDM to reflect changes in basics

Location:

bdm/estim

Files:

• arx.cpp (modified) (1 diff)
• arx.h (modified) (2 diffs)
• merger.cpp (modified) (2 diffs)
• merger.h (modified) (2 diffs)
• mixef.cpp (modified) (3 diffs)
• mixef.h (modified) (3 diffs)

bdm/estim/arx.cpp

@@ -54 +54 @@
         set_parameters ( A0->V,A0->nu );
 }
+
+enorm<ldmat>* ARX::predictor(const RV &rv){
+        mat mu(rv.count(), 1);
+        mat R(rv.count(),rv.count());
+        enorm<ldmat>* tmp; 
+        tmp=new enorm<ldmat>(rv);
+        
+        est.mean_mat(mu,R);
+        tmp->set_parameters(mu.get_col(0),ldmat(R));
+        return tmp;
+}
+
 /*! \brief Return the best structure
 @param Eg a copy of GiW density that is being examined

bdm/estim/arx.h

@@ -66 +66 @@
         const epdf& _epdf() const {return est;}
         double logpred ( const vec &dt ) const;
-        void flatten (BMEF* B ) {
-                ARX* A=dynamic_cast<ARX*>(B);
+        void flatten (const BMEF* B ) {
+                const ARX* A=dynamic_cast<const ARX*>(B);
                 // nu should be equal to B.nu
                 est.pow ( A->nu/nu);
@@ -73 +73 @@
         }
 
+        enorm<ldmat>* predictor(const RV &rv); 
         //! Brute force structure estimation.\return indeces of accepted regressors.
         ivec structure_est ( egiw Eg0 );

bdm/estim/merger.cpp

@@ -18 +18 @@
                         break;
         }
-        return vec(0);
+        return vec ( 0 );
 }
 
 void merger::merge ( const epdf* g0 ) {
-        mat Smp ( rv.count(),Ns );
-        mat lW=zeros ( n,Ns );
+        it_file dbg ( "merger_debug.it" );
 
+        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
+        vec &w = eSmp._w(); //aux
+
+        mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples The last row is ones for the ARX model
+        for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+
+        dbg << Name ( "Smp0" ) << Smp_ex;
+
+        // Stuff for merging
         vec lw_src ( Ns );
         vec lw_mix ( Ns );
+        mat lW=zeros ( n,Ns );
         vec vec0 ( 0 );
 
-        //Make sure g0 is compatible
-        it_assert_debug ( rv.equal(g0->_rv()),"Incompatible g0" );
+        // Initial component in the mixture model
+        mat V0=1e-8*eye ( rv.count() +1 );
+        ARX A0 ( RV ( "{th_r  }", vec_1 ( rv.count() * ( rv.count() +1 ) ) ),\
+                 V0, 4.0 ); //initial guess of Mix: zero mean, large variance
 
-        // Initial component in the mixture model from g0
-        mat V0=1e-8*eye ( rv.count() +1 );
-        ARX A0 ( rv, V0, 3.0 ); //initial guess of Mix: zero mean, large variance
 
-        //generate samples from the proposal g0
-        for ( int i =0; i<Ns; i++ ) {Smp.set_col ( i,g0->sample() );}
 
-        //Initialize Mixture model
-        Mix.init ( &A0, Smp, Nc );
+        // ============= MAIN LOOP ==================
+        bool converged=false;
+        int niter = 0;
+        char str[100];
+        
+        vec Mix_pdf(Ns);
+        while ( !converged ) {
+                //Re-estimate Mix
+                //Re-Initialize Mixture model
+                Mix.init ( &A0, Smp_ex, Nc );
+                Mix.bayesB ( Smp_ex );
 
-        bool converged=false;
-        while ( !converged ) {
-                Mix.bayesB ( Smp );
-                //Generate weight for each particle
+                // Generate new samples
+                eSmp.set_samples ( Mix.predictor ( rv ) );
+                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
+
+                sprintf ( str,"Mpdf%d",niter );
+                for (int i=0;i<Ns;i++){Mix_pdf(i) = Mix.logpred(Smp_ex.get_col(i));}
+                dbg << Name ( str ) << Mix_pdf;
+
+                sprintf ( str,"Smp%d",niter );
+                dbg << Name ( str ) << Smp_ex;
+
+                //Importace weighting
                 for ( int i=0;i<n;i++ ) {
                         //Split according to dependency in rvs
@@ -51 +78 @@
                         if ( rvsinrv ( i ).length() ==rv.count() ) {
                                 // no need for conditioning or marginalization
-                                lw_src=mpdfs ( i )->_epdf().evalpdflog_m ( Smp );
-                                lw_mix=Mix.logpred_m ( Smp );
+                                for ( int j=0;j<Ns;j++ ) {
+                                        lw_src ( j ) =mpdfs ( i )->_epdf().evalpdflog ( Smp ( j ) );
+                                }
                         }
-                        lW.set_row ( i, lw_src-lw_mix );
+                        lW.set_row ( i, lw_src ); // do not divide by mix
                 }
-                converged = true;
+                //Importance of the mixture
+                for ( int j=0;j<Ns;j++ ) {
+                        lw_mix ( j ) =Mix.logpred ( Smp_ex.get_col ( j ) );
+                }
+                sprintf ( str,"lW%d",niter );
+                dbg << Name ( str ) << lW;
+
+                w = lognorm_merge ( lW ); //merge
+
+                sprintf ( str,"w%d",niter );
+                dbg << Name ( str ) << w;
+
+                //Importance weighting
+                w /=exp ( lw_mix ); // hoping that it is not numerically sensitive...
+                //renormalize
+                w /=sum ( 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 ) );}
+
+                sprintf ( str,"Smp_res%d",niter );
+                dbg << Name ( str ) << Smp;
+
+                // ==== stopping rule ===
+                niter++;
+                converged = ( niter>3 );
         }
 

bdm/estim/merger.h

@@ -41 +41 @@
                         compositepdf ( S ), epdf ( getrv ( false ) ),
                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) )
-        { beta=2.0; Ns=100; Nc=10;}
+        {setindices(rv); beta=2.0; Ns=100; Nc=10;}
         //! Set internal parameters used in approximation
         void set_parameters ( double beta0, int Ns0, int Nc0 ) {        beta=beta0;Ns=Ns0;Nc=Nc0;}
@@ -60 +60 @@
         //! weight w is a
         vec sample ( )const { return Mix._epdf().sample();}
-        double evalpdflog ( const vec &dt ) const{ return Mix._epdf().evalpdflog ( dt );}
+        double evalpdflog ( const vec &dt ) const{ vec dtf=zeros(dt.length()+1); dtf.set_subvector(0,dt); return Mix.logpred ( dtf );}
         vec mean()const {return Mix._epdf().mean();}
         //! for future use

bdm/estim/mixef.cpp

@@ -9 +9 @@
         Coms.set_size ( c );
         n=c;
-        weights.set_parameters ( 1/ ( double ) c*ones ( c ) ); //assume at least one observation in each comp.
+        weights.set_parameters ( ones ( c ) ); //assume at least one observation in each comp.
         //est will be done at the end
         //
@@ -30 +30 @@
                 //pick one datum
                 int ind=ndat*UniRNG.sample();
-                Coms ( i )->bayes ( Data.get_col ( ind ),ndat/n );
+                Coms ( i )->bayes ( Data.get_col ( ind ),1.0 );
         }
 
@@ -110 +110 @@
         return log ( exLL );
 }
+
+emix* MixEF::predictor(const RV &rv){
+        Array<epdf*> pC(n);
+        for(int i=0;i<n;i++){pC(i)=Coms(i)->predictor(rv);}
+        emix* tmp;
+        tmp = new emix(rv);
+        tmp->set_parameters(weights._epdf().mean(), pC, false);
+        tmp->ownComs();
+        return tmp;
+}

bdm/estim/mixef.h

@@ -33 +33 @@
 The characteristic feature of this model is that if the exact values of the latent variable were known, estimation of the parameters can be handled by a single model. For example, for the case of mixture models, posterior density for each component parameters would be a BayesianModel from Exponential Family.
 
-This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence approximate estimation project this class itself belongs to the exponential family.
+This class uses EM-style type algorithms for estimation of its parameters. Under this simplification, the posterior density is a product of exponential family members, hence under EM-style approximate estimation this class itself belongs to the exponential family.
 
 TODO: Extend BM to use rvc.
@@ -48 +48 @@
         eprod* est;
         ////!Indeces of component rvc in common rvc
-
+        
         //! Auxiliary function for use in constructors
         void build_est() {
@@ -91 +91 @@
         double logpred ( const vec &dt ) const;
         const epdf& _epdf() const {return *est;}
+        emix* predictor(const RV &rv);
+        //! Flatten the density as if it was not estimated from the data
+        void flatten(double sumw=1.0);
 };