Changeset 477 for library/bdm/stat/exp_family.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/exp_family.cpp (modified) (15 diffs)

library/bdm/stat/exp_family.cpp

@@ -4 +4 @@
 #include "exp_family.h"
 
-namespace bdm{
+namespace bdm {
 
 Uniform_RNG UniRNG;
@@ -12 +12 @@
 using std::cout;
 
-void BMEF::bayes ( const vec &dt ) {this->bayes ( dt,1.0 );};
+void BMEF::bayes ( const vec &dt ) {
+        this->bayes ( dt, 1.0 );
+};
 
 vec egiw::sample() const {
@@ -20 +22 @@
 
 double egiw::evallog_nn ( const vec &val ) const {
-        int vend = val.length()-1;
-
-        if ( dimx==1 ) { //same as the following, just quicker.
+        int vend = val.length() - 1;
+
+        if ( dimx == 1 ) { //same as the following, just quicker.
                 double r = val ( vend ); //last entry!
-                if (r<0) return -inf;
-                vec Psi ( nPsi+dimx );
+                if ( r < 0 ) return -inf;
+                vec Psi ( nPsi + dimx );
                 Psi ( 0 ) = -1.0;
-                Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest
-
-                double Vq=V.qform ( Psi );
-                return -0.5* ( nu*log ( r ) + Vq /r );
-        }
-        else {
-                mat Th= reshape ( val ( 0,nPsi*dimx-1 ),nPsi,dimx );
-                fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
-                double ldetR=R.logdet();
-                if (ldetR) return -inf;
-                mat Tmp=concat_vertical ( -eye ( dimx ),Th );
+                Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest
+
+                double Vq = V.qform ( Psi );
+                return -0.5* ( nu*log ( r ) + Vq / r );
+        } else {
+                mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
+                fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) );
+                double ldetR = R.logdet();
+                if ( ldetR ) return -inf;
+                mat Tmp = concat_vertical ( -eye ( dimx ), Th );
                 fsqmat iR ( dimx );
                 R.inv ( iR );
@@ -48 +49 @@
         const vec& D = V._D();
 
-        double m = nu - nPsi -dimx-1;
+        double m = nu - nPsi - dimx - 1;
 #define log2  0.693147180559945286226763983
 #define logpi 1.144729885849400163877476189
@@ -54 +55 @@
 #define Inf std::numeric_limits<double>::infinity()
 
-        double nkG = 0.5* dimx* ( -nPsi *log2pi + sum ( log ( D ( dimx,D.length()-1 ) ) ) );
+        double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1 ) ) ) );
         // temporary for lgamma in Wishart
-        double lg=0;
-        for ( int i =0; i<dimx;i++ ) {lg+=lgamma ( 0.5* ( m-i ) );}
-
-        double nkW = 0.5* ( m*sum ( log ( D ( 0,dimx-1 ) ) ) ) \
-                     - 0.5*dimx* ( m*log2 + 0.5* ( dimx-1 ) *log2pi )  - lg;
-
-        it_assert_debug ( ( ( -nkG-nkW ) >-Inf ) && ( ( -nkG-nkW ) <Inf ), "ARX improper" );
-        return -nkG-nkW;
+        double lg = 0;
+        for ( int i = 0; i < dimx; i++ ) {
+                lg += lgamma ( 0.5 * ( m - i ) );
+        }
+
+        double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \
+                     - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg;
+
+        it_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
+        return -nkG - nkW;
 }
 
 vec egiw::est_theta() const {
-        if ( dimx==1 ) {
+        if ( dimx == 1 ) {
                 const mat &L = V._L();
                 int end = L.rows() - 1;
 
-                mat iLsub = ltuinv ( L ( dimx,end,dimx,end ) );
-
-                vec L0 = L.get_col(0);
-
-                return iLsub * L0(1,end);
-        }
-        else {
-                it_error("ERROR: est_theta() not implemented for dimx>1");
+                mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
+
+                vec L0 = L.get_col ( 0 );
+
+                return iLsub * L0 ( 1, end );
+        } else {
+                it_error ( "ERROR: est_theta() not implemented for dimx>1" );
                 return 0;
         }
@@ -89 +91 @@
                 int end = D.length() - 1;
 
-                mat Lsub = L ( 1, end, 1, end);
-                mat Dsub = diag(D(1, end));
-
-                return inv( transpose(Lsub) * Dsub * Lsub );
-
-        }
-        else {
-                it_error("ERROR: est_theta_cov() not implemented for dimx>1");
+                mat Lsub = L ( 1, end, 1, end );
+                mat Dsub = diag ( D ( 1, end ) );
+
+                return inv ( transpose ( Lsub ) * Dsub * Lsub );
+
+        } else {
+                it_error ( "ERROR: est_theta_cov() not implemented for dimx>1" );
                 return 0;
         }
@@ -104 +105 @@
 vec egiw::mean() const {
 
-        if ( dimx==1 ) {
-                const vec &D= V._D();
-                int end = D.length()-1;
+        if ( dimx == 1 ) {
+                const vec &D = V._D();
+                int end = D.length() - 1;
 
                 vec m ( dim );
                 m.set_subvector ( 0, est_theta() );
-                m ( end ) = D ( 0 ) / ( nu -nPsi -2*dimx -2 );
+                m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
                 return m;
-        }
-        else {
+        } else {
                 mat M;
                 mat R;
-                mean_mat ( M,R );
-                return cvectorize ( concat_vertical ( M,R ) );
+                mean_mat ( M, R );
+                return cvectorize ( concat_vertical ( M, R ) );
         }
 
@@ -124 +124 @@
 vec egiw::variance() const {
 
-        if ( dimx==1 ) {
-                int l=V.rows();
-                const ldmat tmp(V,linspace(1,l-1));
-                ldmat itmp(l);
-                tmp.inv(itmp);
-                double cove = V._D() ( 0 ) / ( nu -nPsi -2*dimx -2 );
-                
-                vec var(l);
-                var.set_subvector(0,diag(itmp.to_mat())*cove);
-                var(l-1)=cove*cove/( nu -nPsi -2*dimx -2 );
+        if ( dimx == 1 ) {
+                int l = V.rows();
+                const ldmat tmp ( V, linspace ( 1, l - 1 ) );
+                ldmat itmp ( l );
+                tmp.inv ( itmp );
+                double cove = V._D() ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
+
+                vec var ( l );
+                var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
+                var ( l - 1 ) = cove * cove / ( nu - nPsi - 2 * dimx - 2 );
                 return var;
-        }
-        else {it_error("not implemented"); return vec(0);}
+        } else {
+                it_error ( "not implemented" );
+                return vec ( 0 );
+        }
 }
 
 void egiw::mean_mat ( mat &M, mat&R ) const {
-        const mat &L= V._L();
-        const vec &D= V._D();
-        int end = L.rows()-1;
-
-        ldmat ldR ( L ( 0,dimx-1,0,dimx-1 ), D ( 0,dimx-1 ) / ( nu -nPsi -2*dimx -2 ) ); //exp val of R
-        mat iLsub=ltuinv ( L ( dimx,end,dimx,end ) );
+        const mat &L = V._L();
+        const vec &D = V._D();
+        int end = L.rows() - 1;
+
+        ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R
+        mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
 
         // set mean value
-        mat Lpsi = L ( dimx,end,0,dimx-1 );
-        M= iLsub*Lpsi;
-        R= ldR.to_mat()  ;
+        mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+        M = iLsub * Lpsi;
+        R = ldR.to_mat()  ;
 }
 
 vec egamma::sample() const {
-        vec smp ( dim);
+        vec smp ( dim );
         int i;
 
-        for ( i=0; i<dim; i++ ) {
-                if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) {
-                        GamRNG.setup ( alpha ( i ),beta ( i ) );
-                }
-                else {
-                        GamRNG.setup ( alpha ( i ),std::numeric_limits<double>::epsilon() );
+        for ( i = 0; i < dim; i++ ) {
+                if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) {
+                        GamRNG.setup ( alpha ( i ), beta ( i ) );
+                } else {
+                        GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() );
                 }
 #pragma omp critical
@@ -190 +191 @@
         int i;
 
-        if (any(val<=0.)) return -inf;
-        if (any(beta<=0.)) return -inf;
-        for ( i=0; i<dim; i++ ) {
-                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
-        }
-        double tmp=res-lognc();;
-        it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+        if ( any ( val <= 0. ) ) return -inf;
+        if ( any ( beta <= 0. ) ) return -inf;
+        for ( i = 0; i < dim; i++ ) {
+                res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i );
+        }
+        double tmp = res - lognc();;
+        it_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
         return tmp;
 }
@@ -204 +205 @@
         int i;
 
-        for ( i=0; i<dim; i++ ) {
-                res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
+        for ( i = 0; i < dim; i++ ) {
+                res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ;
         }
 
@@ -211 +212 @@
 }
 
-void mgamma::set_parameters(double k0, const vec &beta0) {
+void mgamma::set_parameters ( double k0, const vec &beta0 ) {
         k = k0;
-        iepdf->set_parameters(k * ones(beta0.length()), beta0);
+        iepdf->set_parameters ( k * ones ( beta0.length() ), beta0 );
         dimc = e()->dimension();
 }
 
-ivec eEmp::resample (RESAMPLING_METHOD method) {
-        ivec ind=zeros_i ( n );
+ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+        ivec ind = zeros_i ( n );
         ivec N_babies = zeros_i ( n );
         vec cumDist = cumsum ( w );
         vec u ( n );
-        int i,j,parent;
+        int i, j, parent;
         double u0;
 
         switch ( method ) {
-                case MULTINOMIAL:
-                        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-                        for ( i = n - 2;i >= 0;i-- ) {
-                                u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-                        }
-
-                        break;
-
-                case STRATIFIED:
-
-                        for ( i = 0;i < n;i++ ) {
-                                u ( i ) = ( i + UniRNG.sample() ) / n;
-                        }
-
-                        break;
-
-                case SYSTEMATIC:
-                        u0 = UniRNG.sample();
-
-                        for ( i = 0;i < n;i++ ) {
-                                u ( i ) = ( i + u0 ) / n;
-                        }
-
-                        break;
-
-                default:
-                        it_error ( "PF::resample(): Unknown resampling method" );
+        case MULTINOMIAL:
+                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+
+                for ( i = n - 2; i >= 0; i-- ) {
+                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+                }
+
+                break;
+
+        case STRATIFIED:
+
+                for ( i = 0; i < n; i++ ) {
+                        u ( i ) = ( i + UniRNG.sample() ) / n;
+                }
+
+                break;
+
+        case SYSTEMATIC:
+                u0 = UniRNG.sample();
+
+                for ( i = 0; i < n; i++ ) {
+                        u ( i ) = ( i + u0 ) / n;
+                }
+
+                break;
+
+        default:
+                it_error ( "PF::resample(): Unknown resampling method" );
         }
 
@@ -259 +260 @@
         j = 0;
 
-        for ( i = 0;i < n;i++ ) {
+        for ( i = 0; i < n; i++ ) {
                 while ( u ( i ) > cumDist ( j ) ) j++;
 
@@ -270 +271 @@
 
         // find the first parent;
-        parent=0; while ( N_babies ( parent ) ==0 ) parent++;
+        parent = 0;
+        while ( N_babies ( parent ) == 0 ) parent++;
 
         // Build index
-        for ( i = 0;i < n;i++ ) {
+        for ( i = 0; i < n; i++ ) {
                 if ( N_babies ( i ) > 0 ) {
                         ind ( i ) = i;
                         N_babies ( i ) --; //this index was now replicated;
-                }
-                else {
+                } else {
                         // test if the parent has been fully replicated
                         // if yes, find the next one
-                        while ( ( N_babies ( parent ) ==0 ) || ( N_babies ( parent ) ==1 && parent>i ) ) parent++;
+                        while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
 
                         // Replicate parent
@@ -292 +293 @@
 
         // copy the internals according to ind
-        for ( i=0;i<n;i++ ) {
-                if ( ind ( i ) !=i ) {
-                        samples ( i ) =samples ( ind ( i ) );
-                }
-                w ( i ) = 1.0/n;
+        for ( i = 0; i < n; i++ ) {
+                if ( ind ( i ) != i ) {
+                        samples ( i ) = samples ( ind ( i ) );
+                }
+                w ( i ) = 1.0 / n;
         }
 
@@ -305 +306 @@
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         dim = epdf0->dimension();
-        w=w0;
-        w/=sum ( w0 );//renormalize
-        n=w.length();
+        w = w0;
+        w /= sum ( w0 );//renormalize
+        n = w.length();
         samples.set_size ( n );
         dim = epdf0->dimension();
 
-        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
+        for ( int i = 0; i < n; i++ ) {
+                samples ( i ) = epdf0->sample();
+        }
 }
 
 void eEmp::set_samples ( const epdf* epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
-        w=1;
-        w/=sum ( w );//renormalize
-
-        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
-}
-
-void migamma_ref::from_setting( const Setting &set ) 
-{                               
+        w = 1;
+        w /= sum ( w );//renormalize
+
+        for ( int i = 0; i < n; i++ ) {
+                samples ( i ) = epdf0->sample();
+        }
+}
+
+void migamma_ref::from_setting ( const Setting &set ) {
         vec ref;
-        UI::get( ref, set, "ref" , UI::compulsory);
-        set_parameters(set["k"],ref,set["l"]);
-}
-
-void mlognorm::from_setting( const Setting &set ) 
-{       
-        vec mu0;        
-        UI::get( mu0, set, "mu0", UI::compulsory);
-        set_parameters(mu0.length(),set["k"]);
-        condition(mu0);
+        UI::get ( ref, set, "ref" , UI::compulsory );
+        set_parameters ( set["k"], ref, set["l"] );
+}
+
+void mlognorm::from_setting ( const Setting &set ) {
+        vec mu0;
+        UI::get ( mu0, set, "mu0", UI::compulsory );
+        set_parameters ( mu0.length(), set["k"] );
+        condition ( mu0 );
 }