Changeset 286

Timestamp:

03/05/09 14:03:35 (15 years ago)

Author:

smidl

Message:

make mpdm work again

Files:

• CMakeLists.txt (modified) (1 diff)
• bdm/estim/arx.cpp (modified) (1 diff)
• bdm/estim/arx.h (modified) (3 diffs)
• bdm/estim/merger.cpp (modified) (3 diffs)
• bdm/estim/merger.h (modified) (2 diffs)
• bdm/estim/mixef.cpp (modified) (2 diffs)
• bdm/estim/mixef.h (modified) (5 diffs)
• bdm/stat/emix.h (modified) (12 diffs)
• bdm/stat/libBM.h (modified) (5 diffs)
• bdm/stat/libEF.h (modified) (2 diffs)

CMakeLists.txt

@@ -76 +76 @@
 add_subdirectory (tests)
 add_subdirectory (pmsm)
-#add_subdirectory (mpdm)
+add_subdirectory (mpdm)
 add_subdirectory (library)
 add_subdirectory (doprava)

bdm/estim/arx.cpp

@@ -57 +57 @@
 
 enorm<ldmat>* ARX::epredictor ( const vec &rgr ) const {
-        int dim=est.dimension();
+        int dim=dimx;//est.dimension();
         mat mu ( dim, V.rows()-dim );
         mat R ( dim,dim );

bdm/estim/arx.h

@@ -54 +54 @@
         //!@{
         ARX ( const double frg0=1.0 ) : BMEF ( frg0 ),est (), V ( est._V() ), nu ( est._nu() ) {};
-        ARX ( const ARX &A0 ) : BMEF (),est ( A0.est ), V ( est._V() ), nu ( est._nu() ) {};
+        ARX ( const ARX &A0 ) : BMEF (),est (), V ( est._V() ), nu ( est._nu() ) {
+                set_statistics ( A0.dimx,A0.V,A0.nu );
+                set_parameters(A0.frg);
+        };
         ARX* _copy_() const;
         void set_parameters ( double frg0 ) {frg=frg0;}
@@ -82 +85 @@
         enorm<ldmat>* epredictor ( const vec &rgr ) const;
         //! Predictor for empty regressor
-        enorm<ldmat>* epredictor() const {it_assert_debug ( est.dimension() ==V.rows()-1,"Regressor is not only 1" );return epredictor ( vec_1 ( 1.0 ) );}
+        enorm<ldmat>* epredictor() const {
+                it_assert_debug ( dimx==V.rows()-1,"Regressor is not only 1" );
+                return epredictor ( vec_1 ( 1.0 ) );
+        }
         //! conditional version of the predictor
         mlnorm<ldmat>* predictor() const;
@@ -103 +109 @@
                 if ( _yrv._dsize() !=dimx ) {
                         int i=0;
-                        while ( _yrv._dsize() <dimx ) {_yrv.add ( drv ( vec_1(i) ) );i++;}
+                        while ( _yrv._dsize() <dimx ) {_yrv.add ( drv ( vec_1 ( i ) ) );i++;}
                 }
                 //yrv should be ready by now

bdm/estim/merger.cpp

@@ -51 +51 @@
         // Initial component in the mixture model
         mat V0=1e-8*eye ( dim +1 );
-        ARX A0; 
-        A0.set_statistics(dim, V0, dim*dim +5.0 ); //initial guess of Mix: zero mean, large variance
+        ARX A0; A0.set_statistics(dim, V0); //initial guess of Mix: 
 
         Mix.init ( &A0, Smp_ex, Nc );
@@ -63 +62 @@
         char str[100];
 
-        epdf* Mpred=Mix.epredictor (  );
+        emix* Mpred=Mix.epredictor (  );
         vec Mix_pdf ( Ns );
         while ( !converged ) {
@@ -72 +71 @@
                 delete Mpred;
                 Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
+                Mpred->set_rv(rv); //the predictor predicts rv of this merger
 
                 // This will be active only later in iterations!!!

bdm/estim/merger.h

@@ -52 +52 @@
                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp() {
                 RV ztmp;
+                rv = getrv(false);
+                dim = rv._dsize();
                 // Extend rv by rvc!
                 RV rvc; setrvc ( rv,rvc );
@@ -57 +59 @@
                 for ( int i=0;i<n;i++ ) {
                         //Establich connection between mpdfs and merger
-                        dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
+                        dls ( i ) = new datalink_m2e;dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
                         // find out what is missing in each mpdf
                         ztmp= mpdfs ( i )->_rv();
                         ztmp.add ( mpdfs ( i )->_rvc() );
                         rvzs ( i ) =rv.subt ( ztmp );
-                        zdls ( i ) = new datalink_m2e ( rvzs ( i ), ztmp, rv ) ;
+                        zdls ( i ) = new datalink_m2e; zdls(i)->set_connection ( rvzs ( i ), ztmp, rv ) ;
                 };
                 //Set Default values of parameters

bdm/estim/mixef.cpp

@@ -15 +15 @@
         int ndat = Data.cols();
         //Estimate  Com0 from all data
-        Coms ( 0 ) = ( BMEF* ) Com0->_copy_();
+        Coms ( 0 ) = Com0->_copy_();
 //      Coms(0)->set_evalll(false);
         Coms ( 0 )->bayesB ( Data );
@@ -24 +24 @@
         for ( i=1;i<n;i++ ) {
                 //copy Com0 and create new rvs for them
-                Coms ( i ) = ( BMEF* ) Coms ( 0 )->_copy_ ( true );
+                Coms ( i ) =  Coms ( 0 )->_copy_ ( );
         }
         //Pick some data for each component and update it

bdm/estim/mixef.h

@@ -19 +19 @@
 #include "../stat/emix.h"
 
-namespace bdm{
+namespace bdm {
 
 enum MixEF_METHOD { EM = 0, QB = 1};
@@ -50 +50 @@
         eprod* est;
         ////!Indeces of component rvc in common rvc
-        
+
         //! Flag for a method that is used in the inference
         MixEF_METHOD method;
-        
+
         //! Auxiliary function for use in constructors
         void build_est() {
-                Array<const epdf*> epdfs ( n+1 );
-                for ( int i=0;i<Coms.length();i++ ) {
+                est = new eprod;
+                if ( n>0 ) {
+                        Array<const epdf*> epdfs ( n+1 );
+                        for ( int i=0;i<Coms.length();i++ ) {
 //                      it_assert_debug(!x,"MixEF::MixEF : Incompatible components");
-                        epdfs ( i ) =& ( Coms ( i )->posterior() );
+                                epdfs ( i ) =& ( Coms ( i )->posterior() );
+                        }
+                        // last in the product is the weight
+                        epdfs ( n ) =& ( weights.posterior() );
+                        est->set_parameters ( epdfs, false );
                 }
-                // last in the product is the weight
-                epdfs ( n ) =& ( weights.posterior() );
-                est = new eprod ( epdfs );
         }
 
@@ -69 +72 @@
         //! Full constructor
         MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
-                        BMEF (  ), n ( Coms0.length() ), Coms ( n ),
-                        weights (), method(QB) {
-        //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
+                        BMEF ( ), n ( Coms0.length() ), Coms ( n ),
+                        weights (), method ( QB ) {
+                //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
 
                 for ( int i=0;i<n;i++ ) {Coms ( i ) = ( BMEF* ) Coms0 ( i )->_copy_();}
@@ -79 +82 @@
         MixEF () :
                         BMEF ( ), n ( 0 ), Coms ( 0 ),
-                        weights (),method(QB) {build_est();}
+                        weights (),method ( QB ) {build_est();}
         //! Copy constructor
-        MixEF(const MixEF &M2): BMEF ( ), n ( M2.n ), Coms ( n ),
-                   weights ( M2.weights ), method(M2.method) {
-        //      it_assert_debug ( n>0,"MixEF::MixEF : Empty Component list" );
+        MixEF ( const MixEF &M2 ) : BMEF ( ), 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();
-                   }
+                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.
@@ -108 +111 @@
         emix* epredictor() const;
         //! Flatten the density as if it was not estimated from the data
-        void flatten(const BMEF* M2);
+        void flatten ( const BMEF* M2 );
         //! Access function
-        BMEF* _Coms(int i){return Coms(i);}
-        
+        BMEF* _Coms ( int i ) {return Coms ( i );}
+
         //!Set which method is to be used
-        void set_method(MixEF_METHOD M){method=M;}
+        void set_method ( MixEF_METHOD M ) {method=M;}
 };
 

bdm/stat/emix.h

@@ -18 +18 @@
 //#include <std>
 
-namespace bdm{
+namespace bdm {
 
 //this comes first because it is used inside emix!
@@ -47 +47 @@
         //!Default constructor. By default, the given epdf is not copied!
         //! It is assumed that this function will be used only temporarily.
-        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( rv,rvc,nom0->_rv() ) {
-                ep->set_rv(_rv());
-                set_rvc(nom0->_rv().subt ( rv ) );
-                if ( copy ) {
-//                      nom = nom0->_copy_();
-                        it_error ( "todo" );
-                        destroynom=true;
-                }
-                else {
-                        nom = nom0;
-                        destroynom = false;
-                }
-                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
+        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( ) {
+                // adjust rv and rvc
+                rvc = nom0->_rv().subt ( rv );
+                dimc = rvc._dsize();
+                ep = new epdf;
+                ep->set_parameters(rv._dsize());
+                ep->set_rv(rv);
+                
+                //prepare data structures
+                if ( copy ) {it_error ( "todo" ); destroynom=true; }
+                else { nom = nom0; destroynom = false; }
+                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );               
+                
+                // build denominator
                 den = nom->marginal ( rvc );
+                dl.set_connection(rv,rvc,nom0->_rv());
         };
         double evallogcond ( const vec &val, const vec &cond ) {
@@ -67 +69 @@
                 dl.pushup_cond ( nom_val,val,cond );
                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
-                it_assert_debug(std::isfinite(tmp),"Infinite value");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                 return tmp;
         }
@@ -96 +98 @@
 public:
         //!Default constructor
-        emix (  ) : epdf ( ) {};
+        emix ( ) : epdf ( ) {};
         //! Set weights \c w and components \c Coms
         //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
@@ -103 +105 @@
         vec sample() const;
         vec mean() const {
-                int i; vec mu = zeros ( dim);
+                int i; vec mu = zeros ( dim );
                 for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                 return mu;
@@ -109 +111 @@
         vec variance() const {
                 //non-central moment
-                vec mom2 = zeros(dim);
-                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow(Coms ( i )->mean(),2); }
+                vec mom2 = zeros ( dim );
+                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow ( Coms ( i )->mean(),2 ); }
                 //central moment
-                return mom2-pow(mean(),2);
+                return mom2-pow ( mean(),2 );
         }
         double evallog ( const vec &val ) const {
@@ -118 +120 @@
                 double sum = 0.0;
                 for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
-                if (sum==0.0){sum=std::numeric_limits<double>::epsilon();}
+                if ( sum==0.0 ) {sum=std::numeric_limits<double>::epsilon();}
                 double tmp=log ( sum );
-                it_assert_debug(std::isfinite(tmp),"Infinite");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                 return tmp;
         };
@@ -151 +153 @@
         //!access function
         epdf* _Coms ( int i ) {return Coms ( i );}
+        void set_rv(const RV &rv){
+                epdf::set_rv(rv);
+                for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+        }
 };
 
@@ -174 +180 @@
         mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf (), epdfs ( n ), dls ( n ) {
                 ep=&dummy;
-                RV rv=getrv(true);
-                set_rv(rv);dummy.set_parameters(rv._dsize());
+                RV rv=getrv ( true );
+                set_rv ( rv );dummy.set_parameters ( rv._dsize() );
                 setrvc ( ep->_rv(),rvc );
                 // rv and rvc established = > we can link them with mpdfs
                 for ( int i = 0;i < n;i++ ) {
-                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
+                        dls ( i ) = new datalink_m2m;
+                        dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
                 }
 
@@ -237 +244 @@
         Array<datalink*> dls;
 public:
-        eprod ( const Array<const epdf*> epdfs0 ) : epdf (  ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
+        eprod () : epdfs ( 0 ),dls ( 0 ) {};
+        void set_parameters ( const Array<const epdf*> &epdfs0, bool named=true ) {
+                epdfs=epdfs0;//.set_length ( epdfs0.length() );
+                dls.set_length ( epdfs.length() );
+
                 bool independent=true;
-                for ( int i=0;i<epdfs.length();i++ ) {
-                        independent=rv.add ( epdfs ( i )->_rv() );
-                        it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
-                }
-                for ( int i=0;i<epdfs.length();i++ ) {
-                        dls ( i ) = new datalink ( epdfs ( i )->_rv() , rv );
+                if ( named ) {
+                        for ( int i=0;i<epdfs.length();i++ ) {
+                                independent=rv.add ( epdfs ( i )->_rv() );
+                                it_assert_debug ( independent==true, "eprod:: given components are not independent." );
+                        }
+                        dim=rv._dsize();
+                }
+                else {
+                        dim =0; for ( int i=0;i<epdfs.length();i++ ) {
+                                dim+=epdfs ( i )->dimension();
+                        }
+                }
+                //
+                int cumdim=0;
+                int dimi=0;
+                int i;
+                for ( i=0;i<epdfs.length();i++ ) {
+                        dls ( i ) = new datalink;
+                        if ( named ) {dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );}
+                        else {
+                                dimi = epdfs ( i )->dimension();
+                                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim,cumdim+dimi-1 ) );
+                                cumdim+=dimi;
+                        }
                 }
         }
@@ -260 +289 @@
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls ( i )->pushup ( tmp, pow(pom,2) );
-                }
-                return tmp-pow(mean(),2);
+                        dls ( i )->pushup ( tmp, pow ( pom,2 ) );
+                }
+                return tmp-pow ( mean(),2 );
         }
         vec sample() const {
@@ -277 +306 @@
                         tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
                 }
-                it_assert_debug(std::isfinite(tmp),"Infinite");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                 return tmp;
         }

bdm/stat/libBM.h

@@ -261 +261 @@
         //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
         virtual void qbounds ( vec &lb, vec &ub, double percentage=0.95 ) const {
-                vec mea=mean(); vec std=sqrt(variance()); 
+                vec mea=mean(); vec std=sqrt ( variance() );
                 lb = mea-2*std; ub=mea+2*std;
         };
@@ -401 +401 @@
                 it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
         }
+        //! set connection using indeces
+        void set_connection ( int ds, int us, const ivec &upind ) {
+                downsize = ds;
+                upsize = us;
+                v2v_up= upind;
+
+                it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
+        }
         //! Get val for myself from val of "Up"
         vec pushdown ( const vec &val_up ) {
@@ -425 +433 @@
 
 public:
+        datalink_m2e() {};
         //! Constructor
-        datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
-                        datalink ( rv,rv_up ), condsize ( rvc._dsize() ) {
+        void set_connection ( const RV &rv,  const RV &rvc, const RV &rv_up ) {
+                datalink::set_connection ( rv,rv_up );
+                condsize=  rvc._dsize();
                 //establish v2c connection
                 rvc.dataind ( rv_up, v2c_lo, v2c_up );
@@ -454 +464 @@
 public:
         //! Constructor
-        datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
-                        datalink_m2e ( rv, rvc, rv_up ) {
+        datalink_m2m() {};
+        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
+                datalink_m2e::set_connection ( rv, rvc, rv_up );
                 //establish c2c connection
                 rvc.dataind ( rvc_up, c2c_lo, c2c_up );
@@ -714 +725 @@
                 if ( opt_L_bounds ) {
                         vec ub,lb;
-                        posterior().qbounds(lb,ub);
+                        posterior().qbounds ( lb,ub );
                         L->logit ( LIDs ( 1 ), lb );
                         L->logit ( LIDs ( 2 ), ub );

bdm/stat/libEF.h

@@ -94 +94 @@
 //      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
 
-        BMEF* _copy_ ( bool changerv=false ) const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
+        BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 
@@ -707 +707 @@
 class iW : public epdf{
         public:
-                vec sample(){}
+                vec sample(){return vec();}
 };