Changeset 1013

Timestamp:

05/27/10 23:07:49 (15 years ago)

Author:

smidl

Message:

Flatten has an extra argument

Files:

• applications/bdmtoolbox/mex/mixef_init.cpp (modified) (1 diff)
• applications/bdmtoolbox/tutorial/userguide/mixef_basic.m (modified) (4 diffs)
• library/bdm/estim/arx.cpp (modified) (1 diff)
• library/bdm/estim/arx.h (modified) (1 diff)
• library/bdm/estim/mixtures.cpp (modified) (5 diffs)
• library/bdm/estim/mixtures.h (modified) (2 diffs)
• library/bdm/stat/exp_family.cpp (modified) (2 diffs)
• library/bdm/stat/exp_family.h (modified) (2 diffs)

applications/bdmtoolbox/mex/mixef_init.cpp

@@ -54 +54 @@
         MixEF mix;
         mix.init(&(*com1),Data,no_com);
-        mix.bayes_batch(Data,zeros(0,Data.cols()),ones(Data.cols()));
+//      mix.bayes_batch(Data,zeros(0,Data.cols()));
         shared_ptr<epdf> p=mix.epredictor(vec(0));
         UIFile prf;

applications/bdmtoolbox/tutorial/userguide/mixef_basic.m

@@ -5 +5 @@
 com.rgr = RV({},[],[]);
 
-Data = randn(2,20);
+Data = [randn(2,50) randn(2,50)+10*ones(2,50)];
 
 [Mix0,P0]=mixef_init(Data,com);
-
 % show predictor
 
 Pred = bm_epredictor(Mix0);
+figure(1);
+hold off
+plot(Data(1,:), Data(2,:),'.');
+hold on
 epdf_2dplot(Pred);
 
@@ -29 +32 @@
     % normalize weights
     w = exp(log_w_nn-max(log_w_nn));
-    w = w/sum(w);
+    w = w/sum(w)
     
     for c=1:n
         yt = Data(:,t);
-        [MixQB.Coms{c}]=bm_bayes(MixQB.Coms{c}, yt);
+        [MixQB.Coms{c}]=bm_bayesweighted(MixQB.Coms{c}, yt, [], w(c));
     end
+    MixQB.weights = bm_bayes(MixQB.weights, w);
 end
 toc
@@ -40 +44 @@
 tic
 % should be the same as:
-% MixQBc = bm_bayes_batch(MixQB,Data);
+MixQBc = bm_bayes_batch(MixQB,Data);
 toc
 
 %% display results
 
+PredQB = bm_epredictor(MixQB);
 figure(2);
-PredQB = bm_epredictor(MixQB);
-epdf_2dplot(Pred);
+hold off
+plot(Data(1,:), Data(2,:),'.');
+hold on
+epdf_2dplot(PredQB);
 
 %% relations between components
@@ -56 +63 @@
     end
 end
+figure(3)
 imagesc(Batta_dist);

library/bdm/estim/arx.cpp

@@ -72 +72 @@
 }
 
-void ARX::flatten ( const BMEF* B ) {
+void ARX::flatten ( const BMEF* B , double weight =1.0) {
         const ARX* A = dynamic_cast<const ARX*> ( B );
         // nu should be equal to B.nu
-        est.pow ( A->posterior()._nu() / posterior()._nu() );
+        est.pow ( A->posterior()._nu() / posterior()._nu() *weight);
         if ( evalll ) {
                 last_lognc = est.lognc();

library/bdm/estim/arx.h

@@ -89 +89 @@
         };
         double logpred ( const vec &yt, const vec &cond ) const;
-        void flatten ( const BMEF* B );
+        void flatten ( const BMEF* B , double weight );
         //! Conditioned version of the predictor
         enorm<ldmat>* epredictor ( const vec &rgr ) const;

library/bdm/estim/mixtures.cpp

@@ -9 +9 @@
         Coms.set_size ( c );
         weights.set_parameters ( ones ( c ) ); //assume at least one observation in each comp.
+        multiBM weights0(weights);
         //est will be done at the end
         //
@@ -18 +19 @@
         Coms ( 0 )->bayes_batch ( Data );
         // Flatten it to its original shape
-        Coms ( 0 )->flatten ( Com0 );
+        shared_ptr<BMEF> SharpCom((BMEF*)Coms(0)->_copy());
+        Coms ( 0 )->flatten ( Com0 ); 
 
         //Copy it to the rest
@@ -32 +34 @@
                 } else { // pick at random
                         int ind = (int) floor ( ndat * UniRNG.sample() );
-                        Coms ( i )->bayes ( Data.get_col ( ind ), empty_vec );
+                        Coms ( i )->bayes_weighted ( Data.get_col ( ind ), empty_vec, ndat/Coms.length() );
                 }
-                //flatten back to oringinal
-                Coms ( i )->flatten ( Com0 );
+                //sharpen to the sharp component
+                //Coms ( i )->flatten ( SharpCom.get(), 1.0/Coms.length() );
         }
-
 }
 
-double MixEF::bayes_batch ( const mat &data , const mat &cond, const vec &wData ) {
+double MixEF::bayes_batch_weighted ( const mat &data , const mat &cond, const vec &wData ) {
         int ndat = data.cols();
         int t, i, niter;
@@ -133 +134 @@
 };
 
-void MixEF::bayes ( const mat &data, const vec &cond = empty_vec ) {
-        this->bayes_batch ( data, cond, ones ( data.cols() ) );
-};
-
-
 double MixEF::logpred ( const vec &yt, const vec &cond =empty_vec) const {
 
@@ -162 +158 @@
 }
 
-void MixEF::flatten ( const BMEF* M2 ) {
+void MixEF::flatten ( const BMEF* M2, double weight=1.0 ) {
         const MixEF* Mix2 = dynamic_cast<const MixEF*> ( M2 );
         bdm_assert_debug ( Mix2->Coms.length() == Coms.length(), "Different no of coms" );
         //Flatten each component
         for ( int i = 0; i < Coms.length(); i++ ) {
-                Coms ( i )->flatten ( Mix2->Coms ( i ) );
+                Coms ( i )->flatten ( Mix2->Coms ( i ) , weight);
         }
         //Flatten weights = make them equal!!

library/bdm/estim/mixtures.h

@@ -98 +98 @@
         //! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006
         void bayes ( const vec &yt, const vec &cond );
-        //! EM algorithm
-        void bayes ( const mat &yt, const vec &cond );
         //! batch weighted Bayes rule
-        double bayes_batch ( const mat &yt, const mat &cond, const vec &wData );
+        double bayes_batch_weighted ( const mat &yt, const mat &cond, const vec &wData );
+        double bayes_batch ( const mat &yt, const mat &cond){
+                return bayes_batch_weighted(yt,cond,ones(yt.cols()));
+        };
         double logpred ( const vec &yt, const vec &cond ) const;
         //! return correctly typed posterior (covariant return)
@@ -110 +111 @@
         emix* epredictor(const vec &cond=vec()) const;
         //! Flatten the density as if it was not estimated from the data
-        void flatten ( const BMEF* M2 );
+        void flatten ( const BMEF* M2, double weight );
         //! Access function
         BMEF* _Coms ( int i ) {

library/bdm/stat/exp_family.cpp

@@ -377 +377 @@
         return pred.lognc() - lll;
 }
-void multiBM::flatten ( const BMEF* B ) {
+void multiBM::flatten ( const BMEF* B, double weight=1.0 ) {
         const multiBM* E = dynamic_cast<const multiBM*> ( B );
         // sum(beta) should be equal to sum(B.beta)
         const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
-        beta *= ( sum ( Eb ) / sum ( beta ) );
+        beta *= ( sum ( Eb ) / sum ( beta ) * weight);
         if ( evalll ) {
                 last_lognc = est.lognc();
@@ -768 +768 @@
                 dimc = _beta.length();
         }
-
-};
+}

library/bdm/stat/exp_family.h

@@ -99 +99 @@
 
         //!Flatten the posterior according to the given BMEF (of the same type!)
-        virtual void flatten ( const BMEF * B ) NOT_IMPLEMENTED_VOID;
+        virtual void flatten ( const BMEF * B, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;
 
 
@@ -620 +620 @@
         double logpred ( const vec &yt ) const;
 
-        void flatten ( const BMEF* B );
+        void flatten ( const BMEF* B , double weight);
 
         //! return correctly typed posterior (covariant return)