Changeset 1009

Timestamp:

05/27/10 23:07:16 (14 years ago)

Author:

smidl

Message:

changes in bayes_batch

Location:

library

Files:

• bdm/base/bdmbase.cpp (modified) (1 diff)
• bdm/base/bdmbase.h (modified) (1 diff)
• bdm/estim/arx.cpp (modified) (3 diffs)
• bdm/estim/arx.h (modified) (4 diffs)
• bdm/estim/mixtures.cpp (modified) (7 diffs)
• bdm/estim/mixtures.h (modified) (3 diffs)
• bdm/stat/merger.cpp (modified) (1 diff)
• bdm/stat/merger.h (modified) (1 diff)
• tests/stresssuite/arx_elem_stress.cpp (modified) (1 diff)
• tests/stresssuite/mixtures_stress.cpp (modified) (1 diff)

library/bdm/base/bdmbase.cpp

@@ -575 +575 @@
 }
 
-void BM::bayes_batch ( const mat &Data, const vec &cond ) {
+double BM::bayes_batch ( const mat &Data, const vec &cond ) {
+        double levid=0.0;
         for ( int t = 0; t < Data.cols(); t++ ) {
                 bayes ( Data.get_col ( t ), cond );
-        }
-}
-
-void BM::bayes_batch ( const mat &Data, const mat &Cond ) {
+                levid+=ll;
+        }
+        return levid;
+}
+
+double BM::bayes_batch ( const mat &Data, const mat &Cond ) {
+        double levid=0.0;
         for ( int t = 0; t < Data.cols(); t++ ) {
                 bayes ( Data.get_col ( t ), Cond.get_col ( t ) );
         }
-}
-
-}
+        return levid;
+}
+
+}

library/bdm/base/bdmbase.h

@@ -1293 +1293 @@
         @param dt vector of input data
         */
-        virtual void bayes ( const vec &yt, const vec &cond = empty_vec ) = 0;
+        virtual void bayes ( const vec &yt, const vec &cond = empty_vec ) =0;
         //! Batch Bayes rule (columns of Dt are observations)
-        virtual void bayes_batch ( const mat &Dt, const vec &cond = empty_vec );
+        virtual double bayes_batch ( const mat &Dt, const vec &cond = empty_vec );
         //! Batch Bayes rule (columns of Dt are observations, columns of Cond are conditions)
-        virtual void bayes_batch ( const mat &Dt, const mat &Cond );
+        virtual double bayes_batch ( const mat &Dt, const mat &Cond );
         //! Evaluates predictive log-likelihood of the given data record
         //! I.e. marginal likelihood of the data with the posterior integrated out.
         //! This function evaluates only \f$ y_t \f$, condition is assumed to be the last used in bayes().
         //! See bdm::BM::predictor for conditional version.
-        virtual double logpred ( const vec &yt ) const NOT_IMPLEMENTED(0.0);
+        virtual double logpred ( const vec &yt, const vec &cond ) const NOT_IMPLEMENTED(0.0);
 
         //! Matrix version of logpred
-        vec logpred_mat ( const mat &Yt ) const {
+        vec logpred_mat ( const mat &Yt, const mat &Cond ) const {
                 vec tmp ( Yt.cols() );
                 for ( int i = 0; i < Yt.cols(); i++ ) {
-                        tmp ( i ) = logpred ( Yt.get_col ( i ) );
+                        tmp ( i ) = logpred ( Yt.get_col ( i ), Cond.get_col(i) );
                 }
                 return tmp;

library/bdm/estim/arx.cpp

@@ -3 +3 @@
 
 void ARX::bayes_weighted ( const vec &yt, const vec &cond, const double w ) {
-        bdm_assert_debug ( yt.length() == dimy, "ARX::bayes yt is of size "+num2str(yt.length())+" expected dimension is "+num2str(dimy) );
-        bdm_assert_debug ( cond.length() == rgrlen , "ARX::bayes cond is of size "+num2str(cond.length())+" expected dimension is "+num2str(rgrlen) );
+        bdm_assert_debug ( yt.length() == dimy, "BM::bayes yt is of size "+num2str(yt.length())+" expected dimension is "+num2str(dimy) );
+        bdm_assert_debug ( cond.length() == rgrlen , "BM::bayes cond is of size "+num2str(cond.length())+" expected dimension is "+num2str(rgrlen) );
         
         BMEF::bayes_weighted(yt,cond,w); //potential discount scheduling
@@ -46 +46 @@
 }
 
-double ARX::logpred ( const vec &yt ) const {
+double ARX::logpred ( const vec &yt, const vec &cond ) const {
         egiw pred ( est );
         ldmat &V = pred._V();
@@ -54 +54 @@
         vec dyad_p = dyad;
         dyad_p.set_subvector ( 0, yt );
-
+        dyad_p.set_subvector(dimy,cond);
+        
         if ( frg < 1.0 ) {
                 pred.pow ( frg );

library/bdm/estim/arx.h

@@ -88 +88 @@
                 bayes_weighted ( yt, cond, 1.0 );
         };
-        double logpred ( const vec &yt ) const;
+        double logpred ( const vec &yt, const vec &cond ) const;
         void flatten ( const BMEF* B );
         //! Conditioned version of the predictor
@@ -103 +103 @@
         //! Smarter structure estimation by Ludvik Tesar.\return indices of accepted regressors.
         ivec structure_est_LT ( const egiw &Eg0 );
-        //! reduce to
+        //! reduce structure to the given ivec of matrix V
         void reduce_structure(ivec &inds_in_V){
                 ldmat V = posterior()._V();
@@ -110 +110 @@
                 ldmat newV(V,inds_in_V);
                 est.set_parameters(dimy,newV, posterior()._nu());
+                
+                if (have_constant){
+                        ivec rgr_elem= find(inds_in_V<(V.rows()-1)); // < -- find non-constant 
+                        rgr = rgr.subselect(rgr_elem);
+                        rgrlen = rgr_elem.length();
+                } else{
+                        rgr = rgr.subselect(inds_in_V);
+                }
                 validate();
         }
@@ -184 +192 @@
                 
         } 
+        //! access function
+        RV & _rgr() {return rgr;}
+        bool _have_constant() {return have_constant;}
+        int _rgrlen() {return rgrlen;}
 };
 UIREGISTER ( ARX );

library/bdm/estim/mixtures.cpp

@@ -28 +28 @@
         for ( i = 0; i < Coms.length(); i++ ) {
                 //pick one datum
-                int ind = (int) floor ( ndat * UniRNG.sample() );
-                Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
+                if (ndat==Coms.length()){ //take the ith vector
+                        Coms ( i )->bayes ( Data.get_col ( i ), empty_vec );
+                } else { // pick at random
+                        int ind = (int) floor ( ndat * UniRNG.sample() );
+                        Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
+                }
                 //flatten back to oringinal
                 Coms ( i )->flatten ( Com0 );
@@ -36 +40 @@
 }
 
-void MixEF::bayes_batch ( const mat &data , const mat &cond, const vec &wData ) {
+double MixEF::bayes_batch ( const mat &data , const mat &cond, const vec &wData ) {
         int ndat = data.cols();
         int t, i, niter;
@@ -58 +62 @@
         int maxi;
         double maxll;
+        
+        double levid=0.0;
         //Estim
         while ( !converged ) {
+                levid=0.0;
                 // Copy components back to their initial values
                 // All necessary information is now in w and Coms0.
@@ -68 +75 @@
                         //#pragma omp parallel for
                         for ( i = 0; i < n; i++ ) {
-                                ll ( i ) = Coms ( i )->logpred ( data.get_col ( t ) );
+                                ll ( i ) = Coms ( i )->logpred ( data.get_col ( t ) , empty_vec);
                                 wtmp = 0.0;
                                 wtmp ( i ) = 1.0;
@@ -99 +106 @@
                 // !!!!    note  wData ==> this is extra weight of the data record
                 // !!!!    For typical cases wData=1.
+                vec logevid(n);
                 for ( t = 0; t < ndat; t++ ) {
                         //#pragma omp parallel for
                         for ( i = 0; i < n; i++ ) {
                                 Coms ( i )-> bayes_weighted ( data.get_col ( t ), empty_vec, W ( i, t ) * wData ( t ) );
+                                logevid(i) = Coms(i)->_ll();
                         }
                         weights.bayes ( W.get_col ( t ) * wData ( t ) );
                 }
-
+                levid += weights._ll()+log(weights.posterior().mean() * exp(logevid)); // inner product w*exp(evid)
+                
                 niter++;
                 //TODO better convergence rule.
@@ -116 +126 @@
                 delete Coms0 ( i );
         }
+        return levid;
 }
 
@@ -127 +138 @@
 
 
-double MixEF::logpred ( const vec &yt ) const {
+double MixEF::logpred ( const vec &yt, const vec &cond =empty_vec) const {
 
         vec w = weights.posterior().mean();
         double exLL = 0.0;
         for ( int i = 0; i < Coms.length(); i++ ) {
-                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt ) );
+                exLL += w ( i ) * exp ( Coms ( i )->logpred ( yt , cond ) );
         }
         return log ( exLL );

library/bdm/estim/mixtures.h

@@ -63 +63 @@
         MixEF_METHOD method;
 
+        //! maximum number of iterations
+        int max_niter;
 public:
         //! Full constructor
         MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
                         BMEF ( ), Coms ( Coms0.length() ),
-                        weights (), est(*this), method ( QB ) {
+                        weights (), est(*this), method ( QB ), max_niter(10) {
                 for ( int i = 0; i < Coms0.length(); i++ ) {
                         Coms ( i ) = ( BMEF* ) Coms0 ( i )->_copy();
@@ -78 +80 @@
         MixEF () :
                         BMEF ( ), Coms ( 0 ),
-                        weights (), est(*this), method ( QB ) {
+                        weights (), est(*this), method ( QB ), max_niter(10) {
         }
         //! Copy constructor
         MixEF ( const MixEF &M2 ) : BMEF ( ),  Coms ( M2.Coms.length() ),
-                        weights ( M2.weights ), est(*this), method ( M2.method ) {
+                        weights ( M2.weights ), est(*this), method ( M2.method ), max_niter(M2.max_niter) {
                 for ( int i = 0; i < M2.Coms.length(); i++ ) {
                         Coms ( i ) = (BMEF*) M2.Coms ( i )->_copy();
@@ -99 +101 @@
         void bayes ( const mat &yt, const vec &cond );
         //! batch weighted Bayes rule
-        void bayes_batch ( const mat &yt, const mat &cond, const vec &wData );
-        double logpred ( const vec &yt ) const;
+        double bayes_batch ( const mat &yt, const mat &cond, const vec &wData );
+        double logpred ( const vec &yt, const vec &cond ) const;
         //! return correctly typed posterior (covariant return)
         const eprod_mix& posterior() const {

library/bdm/stat/merger.cpp

@@ -258 +258 @@
                         sprintf ( dbg_str, "pdf%d", niter );
                         for ( int i = 0; i < Npoints; i++ ) {
-                                Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );
+                                Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ), empty_vec );
                         }
                         *dbg_file << Name ( dbg_str ) << Mix_pdf;

library/bdm/stat/merger.h

@@ -233 +233 @@
                 vec dtf = ones ( yt.length() + 1 );
                 dtf.set_subvector ( 0, yt );
-                return Mix.logpred ( dtf );
+                return Mix.logpred ( dtf, vec(0) );
         }
         //!@}

library/tests/stresssuite/arx_elem_stress.cpp

@@ -41 +41 @@
         mlstudent* Ap = Ar.predictor_student();
         vec Ap_x = Ap->evallogcond_mat ( X, empty_vec );
-        vec ll_x = Ar.logpred_mat ( X2 );
+        vec ll_x = Ar.logpred_mat ( X2 , empty_vec);
 
         cout << "normalize : " << xstep*sum ( exp ( Ap_x ) ) << endl;

library/tests/stresssuite/mixtures_stress.cpp

@@ -35 +35 @@
                 for ( double y = yb ( 0 ); y <= yb ( 1 ); y += ystep, j++ ) {
                         rgr ( 1 ) = y;
-                        PPdf ( i, j ) = exp ( M.logpred ( rgr ) );
+                        PPdf ( i, j ) = exp ( M.logpred ( rgr, empty_vec ) );
                 }
         }