Changeset 197

Timestamp:

11/04/08 14:54:33 (15 years ago)

Author:

smidl

Message:

opravy v bdm

Location:

bdm

Files:

• estim/merger.cpp (modified) (5 diffs)
• estim/mixef.cpp (modified) (2 diffs)
• estim/mixef.h (modified) (4 diffs)
• stat/libEF.cpp (modified) (1 diff)
• stat/libEF.h (modified) (3 diffs)

bdm/estim/merger.cpp

@@ -47 +47 @@
         mat V0=1e-8*eye ( rv.count() +1 );
         ARX A0 ( RV ( "{th_r  }", vec_1 ( rv.count() * ( rv.count() +1 ) ) ),\
-                 V0, rv.count() *rv.count() +3.0 ); //initial guess of Mix: zero mean, large variance
+                 V0, rv.count() *rv.count() +5.0 ); //initial guess of Mix: zero mean, large variance
+
+        Mix.init ( &A0, Smp_ex, Nc );
+        //Preserve initial mixture for repetitive estimation via flattening
+        MixEF Mix_init(Mix);
 
         // ============= MAIN LOOP ==================
@@ -59 +63 @@
                 //Re-estimate Mix
                 //Re-Initialize Mixture model
-                Mix.init ( &A0, Smp_ex, Nc );
+                Mix.flatten(&Mix_init);
                 Mix.bayesB ( Smp_ex, w*Ns );
                 Mpred = Mix.predictor ( rv ); // Allocation => must be deleted at the end!!
 
-                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,"Mpred_mean%d",niter );
+                dbg << Name ( str ) << Mpred->mean();
+
+                // Generate new samples
+                eSmp.set_samples ( Mpred );
+                for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
                 sprintf ( str,"Mpdf%d",niter );
@@ -109 +106 @@
                                 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 ) );
@@ -127 +123 @@
                                         lw_src ( k ) += mpdfs ( i )->_epdf().evalpdflog ( dls ( i )->get_val ( Smp ( k ) ) );
                                 }
-                                
+
                         }
                         lW.set_row ( i, lw_src ); // do not divide by mix
@@ -161 +157 @@
                 // ==== stopping rule ===
                 niter++;
-                converged = ( niter>4 );
+                converged = ( niter>9 );
         }
 

bdm/estim/mixef.cpp

@@ -31 +31 @@
                 int ind=floor(ndat*UniRNG.sample());
                 Coms ( i )->bayes ( Data.get_col ( ind ),1.0 );
+                //flatten back to oringinal
+                Coms(i)->flatten(Com0);
         }
 
@@ -139 +141 @@
         return tmp;
 }
+
+void MixEF::flatten(const BMEF* M2){
+        const MixEF* Mix2=dynamic_cast<const MixEF*>(M2);
+        it_assert_debug(Mix2->n==n,"Different no of coms");
+        //Flatten each component
+        for (int i=0; i<n;i++){
+                Coms(i)->flatten(Mix2->Coms(i));
+        }
+        //Faltten weights
+        weights.flatten(&(Mix2->weights));
+}

bdm/estim/mixef.h

@@ -39 +39 @@
 TODO: Extend BM to use rvc.
 */
-class MixEF: public BM {
+class MixEF: public BMEF {
 protected:
         //!Number of components
@@ -69 +69 @@
         //! Full constructor
         MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
-                        BM ( RV() ), n ( Coms0.length() ), Coms ( n ),
+                        BMEF ( RV() ), n ( Coms0.length() ), Coms ( n ),
                         weights ( RV ( "{w }", vec_1 ( n ) ),alpha0 ), method(QB) {
         //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
@@ -76 +76 @@
                 build_est();
         };
+        //! Constructor of empty mixture
         MixEF () :
-                        BM ( RV() ), n ( 0 ), Coms ( 0 ),
+                        BMEF ( RV() ), n ( 0 ), Coms ( 0 ),
                         weights ( RV ( "{w }", vec_1 ( 0 ) ),vec ( 0 ) ),method(QB) {build_est();}
+        //! Copy constructor
+        MixEF(const MixEF &M2): BMEF ( RV() ), n ( M2.n ), Coms ( n ),
+                   weights ( M2.weights ), method(M2.method) {
+        //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
+
+                           for ( int i=0;i<n;i++ ) {Coms ( i ) = M2.Coms ( i )->_copy_();}
+                           build_est();
+                   }
         //! Initializing the mixture by a random pick of centroids from data
         //! \param Com0 Initial component - necessary to determine its type.
@@ -98 +107 @@
         emix* predictor(const RV &rv);
         //! Flatten the density as if it was not estimated from the data
-        void flatten(double sumw=1.0);
+        void flatten(const BMEF* M2);
         //! Access function
         BMEF* _Coms(int i){return Coms(i);}

bdm/stat/libEF.cpp

@@ -98 +98 @@
         R= ldR.to_mat()  ;
 }
-
-void egiw::pow(double p){V*=p;nu*=p;}
 
 vec egamma::sample() const {

bdm/stat/libEF.h

@@ -93 +93 @@
         virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
         //!Flatten the posterior as if to keep nu0 data
-        virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+        
+        BMEF* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 
@@ -190 +192 @@
         //! returns a pointer to the internal statistics. Use with Care!
         double& _nu() {return nu;}
-        void pow ( double p );
+        void pow ( double p ){V*=p;nu*=p;};
 };
 
@@ -266 +268 @@
         }
         void flatten ( const BMEF* B ) {
-                const eDirich* E=dynamic_cast<const eDirich*> ( B );
+                const multiBM* E=dynamic_cast<const multiBM*> ( B );
                 // sum(beta) should be equal to sum(B.beta)
-                const vec &Eb=const_cast<eDirich*> ( E )->_beta();
-                est.pow ( sum ( beta ) /sum ( Eb ) );
+                const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
+                beta*= ( sum ( Eb) /sum ( beta ) );
                 if ( evalll ) {last_lognc=est.lognc();}
         }