Changeset 886 for library/bdm/stat/emix.cpp

Timestamp:

03/29/10 23:01:49 (14 years ago)

Author:

smidl

Message:

epdf and emix now have _base classes

Files:

• library/bdm/stat/emix.cpp (modified) (11 diffs)

library/bdm/stat/emix.cpp

@@ -3 +3 @@
 namespace bdm {
 
-void emix::validate (){
-        bdm_assert ( Coms.length() > 0, "There has to be at least one component." );
-
-        bdm_assert ( Coms.length() == w.length(), "It is obligatory to define weights of all the components." );
+void emix_base::validate (){
+        bdm_assert ( no_coms() > 0, "There has to be at least one component." );
+
+        bdm_assert ( no_coms() == w.length(), "It is obligatory to define weights of all the components." );
 
         double sum_w = sum ( w );
@@ -12 +12 @@
         w = w / sum_w;
 
-        dim = Coms ( 0 )->dimension();
-        RV rv_tmp = Coms ( 0 )->_rv() ;
-        bool isnamed = Coms ( 0 )->isnamed();
-        for ( int i = 1; i < Coms.length(); i++ ) {
-                bdm_assert ( dim == ( Coms ( i )->dimension() ), "Component sizes do not match!" );
-                isnamed &= Coms(i)->isnamed() & Coms(i)->_rv().equal(rv_tmp);
+        dim = component ( 0 )->dimension();
+        RV rv_tmp = component ( 0 )->_rv() ;
+        bool isnamed = component( 0 )->isnamed();
+        for ( int i = 1; i < no_coms(); i++ ) {
+                bdm_assert ( dim == ( component ( i )->dimension() ), "Component sizes do not match!" );
+                isnamed &= component(i)->isnamed() & component(i)->_rv().equal(rv_tmp);
         }
         if (isnamed)
@@ -23 +23 @@
 }
 
-void emix::from_setting ( const Setting &set ) {
-        UI::get ( Coms, set, "pdfs", UI::compulsory );
-
-        if ( !UI::get ( w, set, "weights", UI::optional ) ) {
-                int len = Coms.length();
-                w.set_length ( len );
-                w = 1.0 / len;
-        }
-}
-
-
-vec emix::sample() const {
+
+
+vec emix_base::sample() const {
         //Sample which component
         vec cumDist = cumsum ( w );
@@ -46 +37 @@
         }
 
-        return Coms ( i )->sample();
-}
-
-vec emix::mean() const {
+        return component ( i )->sample();
+}
+
+vec emix_base::mean() const {
         int i;
         vec mu = zeros ( dim );
         for ( i = 0; i < w.length(); i++ ) {
-                mu += w ( i ) * Coms ( i )->mean();
+                mu += w ( i ) * component ( i )->mean();
         }
         return mu;
 }
 
-vec emix::variance() const {
+vec emix_base::variance() const {
         //non-central moment
         vec mom2 = zeros ( dim );
         for ( int i = 0; i < w.length(); i++ ) {
-                mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+                mom2 += w ( i ) * ( component( i )->variance() + pow ( component ( i )->mean(), 2 ) );
         }
         //central moment
@@ -68 +59 @@
 }
 
-double emix::evallog ( const vec &val ) const {
+double emix_base::evallog ( const vec &val ) const {
         int i;
         double sum = 0.0;
         for ( i = 0; i < w.length(); i++ ) {
-                sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+                sum += w ( i ) * exp ( component ( i )->evallog ( val ) );
         }
         if ( sum == 0.0 ) {
@@ -82 +73 @@
 }
 
-vec emix::evallog_mat ( const mat &Val ) const {
+vec emix_base::evallog_mat ( const mat &Val ) const {
         vec x = zeros ( Val.cols() );
         for ( int i = 0; i < w.length(); i++ ) {
-                x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
+                x += w ( i ) * exp ( component( i )->evallog_mat ( Val ) );
         }
         return log ( x );
 };
 
-mat emix::evallog_coms ( const mat &Val ) const {
+mat emix_base::evallog_coms ( const mat &Val ) const {
         mat X ( w.length(), Val.cols() );
         for ( int i = 0; i < w.length(); i++ ) {
-                X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) );
+                X.set_row ( i, w ( i ) *exp ( component( i )->evallog_mat ( Val ) ) );
         }
         return X;
 }
 
-shared_ptr<epdf> emix::marginal ( const RV &rv ) const {
+shared_ptr<epdf> emix_base::marginal ( const RV &rv ) const {
         emix *tmp = new emix();
         shared_ptr<epdf> narrow ( tmp );
@@ -105 +96 @@
 }
 
-void emix::marginal ( const RV &rv, emix &target ) const {
+void emix_base::marginal ( const RV &rv, emix &target ) const {
         bdm_assert ( isnamed(), "rvs are not assigned" );
 
-        Array<shared_ptr<epdf> > Cn ( Coms.length() );
-        for ( int i = 0; i < Coms.length(); i++ ) {
-                Cn ( i ) = Coms ( i )->marginal ( rv );
+        Array<shared_ptr<epdf> > Cn ( no_coms() );
+        for ( int i = 0; i < no_coms(); i++ ) {
+                Cn ( i ) = component ( i )->marginal ( rv );
         }
 
@@ -118 +109 @@
 }
 
-shared_ptr<pdf> emix::condition ( const RV &rv ) const {
+shared_ptr<pdf> emix_base::condition ( const RV &rv ) const {
         bdm_assert ( isnamed(), "rvs are not assigned" );
         mratio *tmp = new mratio ( this, rv );
@@ -124 +115 @@
 }
 
-void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-        w = w0 / sum ( w0 );
-        int i;
-        for ( i = 0; i < w.length(); i++ ) {
-                bdm_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
-        }
-        if ( copy ) {
-                Coms.set_length ( Coms0.length() );
-                for ( i = 0; i < w.length(); i++ ) {
-                        bdm_error ( "Not implemented" );
-                        // *Coms ( i ) = *Coms0 ( i );
-                }
-                destroyComs = true;
-        } else {
-                Coms = Coms0;
-                destroyComs = false;
-        }
-}
-
-void    egiwmix::validate (){
-   dim = Coms ( 0 )->dimension();
-}
-
-vec egiwmix::sample() const {
-        //Sample which component
-        vec cumDist = cumsum ( w );
-        double u0;
-#pragma omp critical
-        u0 = UniRNG.sample();
-
-        int i = 0;
-        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
-                i++;
-        }
-
-        return Coms ( i )->sample();
-}
-
-vec egiwmix::mean() const {
-        int i;
-        vec mu = zeros ( dim );
-        for ( i = 0; i < w.length(); i++ ) {
-                mu += w ( i ) * Coms ( i )->mean();
-        }
-        return mu;
-}
-
-vec egiwmix::variance() const {
-        // non-central moment
-        vec mom2 = zeros ( dim );
-        for ( int i = 0; i < w.length(); i++ ) {
-                // pow is overloaded, we have to use another approach
-                mom2 += w ( i ) * ( Coms ( i )->variance() + elem_mult ( Coms ( i )->mean(), Coms ( i )->mean() ) );
-        }
-        // central moment
-        // pow is overloaded, we have to use another approach
-        return mom2 - elem_mult ( mean(), mean() );
-}
-
-shared_ptr<epdf> egiwmix::marginal ( const RV &rv ) const {
-        emix *tmp = new emix();
-        shared_ptr<epdf> narrow ( tmp );
-        marginal ( rv, *tmp );
-        return narrow;
-}
-
-void egiwmix::marginal ( const RV &rv, emix &target ) const {
-        bdm_assert_debug ( isnamed(), "rvs are not assigned" );
-
-        Array<shared_ptr<epdf> > Cn ( Coms.length() );
-        for ( int i = 0; i < Coms.length(); i++ ) {
-                Cn ( i ) = Coms ( i )->marginal ( rv );
-        }
-
-        target._w() = w;
-        target._Coms() = Cn;
-        target.validate();
-}
-
-egiw*   egiwmix::approx() {
-        // NB: dimx == 1 !!!
-        // The following code might look a bit spaghetti-like,
-        // consult Dedecius, K. et al.: Partial forgetting in AR models.
-
-        double sumVecCommon;                            // common part for many terms in eq.
-        int len = w.length();                           // no. of mix components
-        int dimLS = Coms ( 1 )->_V()._D().length() - 1;         // dim of LS
-        vec vecNu ( len );                                      // vector of dfms of components
-        vec vecD ( len );                                       // vector of LS reminders of comps.
-        vec vecCommon ( len );                          // vector of common parts
-        mat matVecsTheta;                               // matrix which rows are theta vects.
-
-        // fill in the vectors vecNu, vecD and matVecsTheta
-        for ( int i = 0; i < len; i++ ) {
-                vecNu.shift_left ( Coms ( i )->_nu() );
-                vecD.shift_left ( Coms ( i )->_V()._D() ( 0 ) );
-                matVecsTheta.append_row ( Coms ( i )->est_theta() );
-        }
-
-        // calculate the common parts and their sum
-        vecCommon = elem_mult ( w, elem_div ( vecNu, vecD ) );
-        sumVecCommon = sum ( vecCommon );
-
-        // LS estimator of theta
-        vec aprEstTheta ( dimLS );
-        aprEstTheta.zeros();
-        for ( int i = 0; i < len; i++ ) {
-                aprEstTheta +=  matVecsTheta.get_row ( i ) * vecCommon ( i );
-        }
-        aprEstTheta /= sumVecCommon;
-
-
-        // LS estimator of dfm
-        double aprNu;
-        double A = log ( sumVecCommon );                // Term 'A' in equation
-
-        for ( int i = 0; i < len; i++ ) {
-                A += w ( i ) * ( log ( vecD ( i ) ) - psi ( 0.5 * vecNu ( i ) ) );
-        }
-
-        aprNu = ( 1 + sqrt ( 1 + 2 * ( A - LOG2 ) / 3 ) ) / ( 2 * ( A - LOG2 ) );
-
-
-        // LS reminder (term D(0,0) in C-syntax)
-        double aprD = aprNu / sumVecCommon;
-
-        // Aproximation of cov
-        // the following code is very numerically sensitive, thus
-        // we have to eliminate decompositions etc. as much as possible
-        mat aprC = zeros ( dimLS, dimLS );
-        for ( int i = 0; i < len; i++ ) {
-                aprC += Coms ( i )->est_theta_cov().to_mat() * w ( i );
-                vec tmp = ( matVecsTheta.get_row ( i ) - aprEstTheta );
-                aprC += vecCommon ( i ) * outer_product ( tmp, tmp );
-        }
-
-        // Construct GiW pdf :: BEGIN
-        ldmat aprCinv ( inv ( aprC ) );
-        vec D = concat ( aprD, aprCinv._D() );
-        mat L = eye ( dimLS + 1 );
-        L.set_submatrix ( 1, 0, aprCinv._L() * aprEstTheta );
-        L.set_submatrix ( 1, 1, aprCinv._L() );
-        ldmat aprLD ( L, D );
-
-        egiw* aprgiw = new egiw ( 1, aprLD, aprNu );
-        return aprgiw;
-};
+void emix::from_setting ( const Setting &set ) {
+        UI::get ( Coms, set, "pdfs", UI::compulsory );
+        UI::get ( w, set, "weights", UI::compulsory );
+}
+
+void    emix::validate (){
+        emix_base::validate();
+        dim = Coms ( 0 )->dimension();
+}
+
 
 double mprod::evallogcond ( const vec &val, const vec &cond ) {
@@ -378 +232 @@
 }
 
-vec eprod::mean() const {
+vec eprod_base::mean() const {
         vec tmp ( dim );
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->mean();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor( i )->mean();
                 dls ( i )->pushup ( tmp, pom );
         }
@@ -387 +241 @@
 }
 
-vec eprod::variance() const {
+vec eprod_base::variance() const {
         vec tmp ( dim ); //second moment
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->mean();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor ( i )->mean();
                 dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
         }
         return tmp - pow ( mean(), 2 );
 }
-vec eprod::sample() const {
+vec eprod_base::sample() const {
         vec tmp ( dim );
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->sample();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor ( i )->sample();
                 dls ( i )->pushup ( tmp, pom );
         }
         return tmp;
 }
-double eprod::evallog ( const vec &val ) const {
+double eprod_base::evallog ( const vec &val ) const {
         double tmp = 0;
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-        }
-        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-        return tmp;
-}
-
-}
-// mprod::mprod ( Array<pdf*> mFacs, bool overlap) : pdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), pdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), irvcs_rvc ( n ) {
-//              int i;
-//              bool rvaddok;
-//              // Create rv
-//              for ( i = 0;i < n;i++ ) {
-//                      rvaddok=rv.add ( pdfs ( i )->_rv() ); //add rv to common rvs.
-//                      // If rvaddok==false, pdfs overlap => assert error.
-//                      epdfs ( i ) = & ( pdfs ( i )->posterior() ); // add pointer to epdf
-//              };
-//              // Create rvc
-//              for ( i = 0;i < n;i++ ) {
-//                      rvc.add ( pdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
-//              };
-//
-// //           independent = true;
-//              //test rvc of pdfs and fill rvinds
-//              for ( i = 0;i < n;i++ ) {
-//                      // find ith rv in common rv
-//                      rvsinrv ( i ) = pdfs ( i )->_rv().dataind ( rv );
-//                      // find ith rvc in common rv
-//                      rvcinrv ( i ) = pdfs ( i )->_rvc().dataind ( rv );
-//                      // find ith rvc in common rv
-//                      irvcs_rvc ( i ) = pdfs ( i )->_rvc().dataind ( rvc );
-//                      //
-// /*                   if ( rvcinrv ( i ).length() >0 ) {independent = false;}
-//                      if ( irvcs_rvc ( i ).length() >0 ) {independent = false;}*/
-//              }
-//      };
-
+        for ( int i = 0; i < no_factors(); i++ ) {
+                tmp += factor ( i )->evallog ( dls ( i )->pushdown ( val ) );
+        }
+        //bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+        return tmp;
+}
+
+}
+