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

Timestamp:

08/05/09 14:40:03 (15 years ago)

Author:

mido

Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Files:

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

library/bdm/stat/emix.cpp

@@ -1 +1 @@
 #include "emix.h"
 
-namespace bdm{
+namespace bdm {
 
 void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
-        w = w0/sum ( w0 );
-        dim = Coms0(0)->dimension();
-        bool isnamed = Coms0(0)->isnamed();
+        w = w0 / sum ( w0 );
+        dim = Coms0 ( 0 )->dimension();
+        bool isnamed = Coms0 ( 0 )->isnamed();
         int i;
         RV tmp_rv;
-        if (isnamed) tmp_rv=Coms0(0)->_rv();
-        
-        for ( i=0;i<w.length();i++ ) {
-                it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" ); 
-                it_assert_debug ( !isnamed || tmp_rv.equal( Coms0 ( i )->_rv() ),"Component RVs do not match!" );
+        if ( isnamed ) tmp_rv = Coms0 ( 0 )->_rv();
+
+        for ( i = 0; i < w.length(); i++ ) {
+                it_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
+                it_assert_debug ( !isnamed || tmp_rv.equal ( Coms0 ( i )->_rv() ), "Component RVs do not match!" );
         }
         if ( copy ) {
-                Coms.set_length(Coms0.length());
-                for ( i=0;i<w.length();i++ ) {it_error("Not imp...");
-                        *Coms ( i ) =*Coms0 ( i );}
-                destroyComs=true;
-        }
-        else {
+                Coms.set_length ( Coms0.length() );
+                for ( i = 0; i < w.length(); i++ ) {
+                        it_error ( "Not imp..." );
+                        *Coms ( i ) = *Coms0 ( i );
+                }
+                destroyComs = true;
+        } else {
                 Coms = Coms0;
-                destroyComs=false;
-        }
-        if (isnamed) epdf::set_rv(tmp_rv); //coms aer already OK, no need for set_rv
+                destroyComs = false;
+        }
+        if ( isnamed ) epdf::set_rv ( tmp_rv ); //coms aer already OK, no need for set_rv
 }
 
@@ -32 +33 @@
         vec cumDist = cumsum ( w );
         double u0;
-        #pragma omp critical
+#pragma omp critical
         u0 = UniRNG.sample();
 
-        int i=0;
-        while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+        int i = 0;
+        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+                i++;
+        }
 
         return Coms ( i )->sample();
 }
 
-emix* emix::marginal(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-                        
-        Array<epdf*> Cn(Coms.length());
-        for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+emix* emix::marginal ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+        Array<epdf*> Cn ( Coms.length() );
+        for ( int i = 0; i < Coms.length(); i++ ) {
+                Cn ( i ) = Coms ( i )->marginal ( rv );
+        }
         emix* tmp = new emix();
-        tmp->set_parameters(w,Cn,false);
+        tmp->set_parameters ( w, Cn, false );
         tmp->ownComs();
         return tmp;
 }
 
-mratio* emix::condition(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-        return new mratio(this,rv);
+mratio* emix::condition ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+        return new mratio ( this, rv );
 };
 
 void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-        w = w0/sum ( w0 );
-        dim = Coms0(0)->dimension();
+        w = w0 / sum ( w0 );
+        dim = Coms0 ( 0 )->dimension();
         int i;
-        for ( i=0;i<w.length();i++ ) {
-                it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" );
+        for ( i = 0; i < w.length(); i++ ) {
+                it_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++ ) {it_error("Not imp...");
-                        *Coms ( i ) =*Coms0 ( i );}
-                destroyComs=true;
-        }
-        else {
+                Coms.set_length ( Coms0.length() );
+                for ( i = 0; i < w.length(); i++ ) {
+                        it_error ( "Not imp..." );
+                        *Coms ( i ) = *Coms0 ( i );
+                }
+                destroyComs = true;
+        } else {
                 Coms = Coms0;
-                destroyComs=false;
+                destroyComs = false;
         }
 }
@@ -80 +86 @@
         vec cumDist = cumsum ( w );
         double u0;
-        #pragma omp critical
+#pragma omp critical
         u0 = UniRNG.sample();
 
-        int i=0;
-        while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+        int i = 0;
+        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+                i++;
+        }
 
         return Coms ( i )->sample();
@@ -90 +98 @@
 
 vec egiwmix::mean() const {
-        int i; vec mu = zeros ( dim );
-        for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
+        int i;
+        vec mu = zeros ( dim );
+        for ( i = 0; i < w.length(); i++ ) {
+                mu += w ( i ) * Coms ( i )->mean();
+        }
         return mu;
 }
@@ -98 +109 @@
         // non-central moment
         vec mom2 = zeros ( dim );
-        for ( int i = 0;i < w.length();i++ ) {
+        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() )); 
+                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() );
-}
-
-emix* egiwmix::marginal(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-                        
-        Array<epdf*> Cn(Coms.length());
-        for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+        return mom2 - elem_mult ( mean(), mean() );
+}
+
+emix* egiwmix::marginal ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+        Array<epdf*> Cn ( Coms.length() );
+        for ( int i = 0; i < Coms.length(); i++ ) {
+                Cn ( i ) = Coms ( i )->marginal ( rv );
+        }
         emix* tmp = new emix();
-        tmp->set_parameters(w,Cn,false);
+        tmp->set_parameters ( w, Cn, false );
         tmp->ownComs();
         return tmp;
@@ -124 +137 @@
 
         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
+        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()  );
+        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);
+        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 );
+        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) );
+        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 ) );
 
 
@@ -167 +181 @@
         // 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);
+        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);
+        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;
 };