Changeset 737 for library/bdm/stat/exp_family.cpp

Timestamp:

11/25/09 12:14:38 (15 years ago)

Author:

mido

Message:

ASTYLER RUN OVER THE WHOLE LIBRARY, JUPEE

Files:

• library/bdm/stat/exp_family.cpp (modified) (8 diffs)

library/bdm/stat/exp_family.cpp

@@ -12 +12 @@
 using std::cout;
 
-        ///////////
-
-void BMEF::bayes( const vec &yt, const vec &cond ) {
+///////////
+
+void BMEF::bayes ( const vec &yt, const vec &cond ) {
         this->bayes_weighted ( yt, cond, 1.0 );
 };
 
-void egiw::set_parameters (int dimx0, ldmat V0, double nu0) {
+void egiw::set_parameters ( int dimx0, ldmat V0, double nu0 ) {
         dimx = dimx0;
         nPsi = V0.rows() - dimx;
-        dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
-        
+        dim = dimx * ( dimx + nPsi ); // size(R) + size(Theta)
+
         V = V0;
-        if (nu0 < 0) {
+        if ( nu0 < 0 ) {
                 nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
                 // terms before that are sufficient for finite normalization
@@ -35 +35 @@
         mat M;
         chmat R;
-        sample_mat(M,R);
-        
-        return concat (cvectorize(M),cvectorize(R.to_mat()));
-}
-
-mat egiw::sample_mat(int n) const {
+        sample_mat ( M, R );
+
+        return concat ( cvectorize ( M ), cvectorize ( R.to_mat() ) );
+}
+
+mat egiw::sample_mat ( int n ) const {
         // TODO - correct approach - convert to product of norm * Wishart
         mat M;
         ldmat Vz;
         ldmat Lam;
-        factorize(M,Vz,Lam);
-        
-        chmat ChLam(Lam.to_mat());
+        factorize ( M, Vz, Lam );
+
+        chmat ChLam ( Lam.to_mat() );
         chmat iChLam;
-        ChLam.inv(iChLam);
-        
-        eWishartCh Omega; //inverse Wishart, result is R, 
-        Omega.set_parameters(iChLam,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
-        
+        ChLam.inv ( iChLam );
+
+        eWishartCh Omega; //inverse Wishart, result is R,
+        Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+
         mat OmChi;
-        mat Z(M.rows(),M.cols());
+        mat Z ( M.rows(), M.cols() );
 
         mat Mi;
         mat RChiT;
-        mat tmp(dimension(), n);
-        for (int i=0; i<n;i++){
-                OmChi=Omega.sample_mat();
-                RChiT=inv(OmChi);
-                Z=randn(M.rows(), M.cols());
-                Mi = M + RChiT * Z * inv(Vz._L().T() *diag(sqrt(Vz._D())));
-        
-                tmp.set_col(i,concat (cvectorize(Mi),cvectorize(RChiT*RChiT.T())));
+        mat tmp ( dimension(), n );
+        for ( int i = 0; i < n; i++ ) {
+                OmChi = Omega.sample_mat();
+                RChiT = inv ( OmChi );
+                Z = randn ( M.rows(), M.cols() );
+                Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
+
+                tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) );
         }
         return tmp;
 }
 
-void egiw::sample_mat(mat &Mi, chmat &Ri)const{
-        
+void egiw::sample_mat ( mat &Mi, chmat &Ri ) const {
+
         // TODO - correct approach - convert to product of norm * Wishart
         mat M;
         ldmat Vz;
         ldmat Lam;
-        factorize(M,Vz,Lam);
-        
+        factorize ( M, Vz, Lam );
+
         chmat Ch;
-        Ch.setCh(Lam._L()*diag(sqrt(Lam._D())));
+        Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) );
         chmat iCh;
-        Ch.inv(iCh);
-        
-        eWishartCh Omega; //inverse Wishart, result is R, 
-        Omega.set_parameters(iCh,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
-        
+        Ch.inv ( iCh );
+
+        eWishartCh Omega; //inverse Wishart, result is R,
+        Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+
         chmat Omi;
-        Omi.setCh(Omega.sample_mat());
-        
-        mat Z=randn(M.rows(), M.cols());
-        Mi = M+ Omi._Ch() * Z * inv(Vz._L()*diag(sqrt(Vz._D())));
-        Omi.inv(Ri);
+        Omi.setCh ( Omega.sample_mat() );
+
+        mat Z = randn ( M.rows(), M.cols() );
+        Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+        Omi.inv ( Ri );
 }
 
@@ -140 +140 @@
 
 //      bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
-        if (-nkG - nkW==Inf){
-                cout << "??" <<endl;
+        if ( -nkG - nkW == Inf ) {
+                cout << "??" << endl;
         }
         return -nkG - nkW;
@@ -162 +162 @@
 }
 
-void egiw::factorize(mat &M, ldmat &Vz, ldmat &Lam) const{      
+void egiw::factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const {
         const mat &L = V._L();
         const vec &D = V._D();
         int end = L.rows() - 1;
-        
-        Vz=ldmat(L ( dimx, end, dimx, end ), D(dimx,end));
-        mat iLsub = ltuinv ( Vz._L());
+
+        Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) );
+        mat iLsub = ltuinv ( Vz._L() );
         // set mean value
         mat Lpsi = L ( dimx, end, 0, dimx - 1 );
         M = iLsub * Lpsi;
-        
-        Lam =ldmat ( L (0, dimx-1, 0, dimx-1 ), D (0, dimx-1 ) );  //exp val of R
-         if (1){ // test with Peterka
-                 mat VF=V.to_mat();
-                 mat Vf=VF(0,dimx-1,0, dimx-1);
-                 mat Vzf = VF(dimx,end,0,dimx-1);
-                 mat VZ = VF(dimx,end,dimx,end);
-                 
-                 mat Lam2 = Vf-Vzf.T()*inv(VZ)*Vzf;
-         }
+
+        Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) );  //exp val of R
+        if ( 1 ) { // test with Peterka
+                mat VF = V.to_mat();
+                mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 );
+                mat Vzf = VF ( dimx, end, 0, dimx - 1 );
+                mat VZ = VF ( dimx, end, dimx, end );
+
+                mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf;
+        }
 }
 
@@ -193 +193 @@
 //              mat Dsub = diag ( D ( 1, end ) );
 
-                ldmat LD(inv(Lsub).T(), 1.0/D(1,end));
+                ldmat LD ( inv ( Lsub ).T(), 1.0 / D ( 1, end ) );
                 return LD;
 
@@ -217 +217 @@
                 mat R;
                 mean_mat ( M, R );
-                return concat( cvectorize (  M),cvectorize( R ) );
+                return concat ( cvectorize ( M ), cvectorize ( R ) );
         }
 
@@ -229 +229 @@
         ldmat itmp ( l );
         tmp.inv ( itmp );
-        
+
         // following Wikipedia notation
-        // m=nu-nPsi-dimx-1, p=dimx 
-        double mp1p=nu-nPsi-2*dimx; // m-p+1
-        double mp1m=mp1p-2;         // m-p-1
-        
+        // m=nu-nPsi-dimx-1, p=dimx
+        double mp1p = nu - nPsi - 2 * dimx; // m-p+1
+        double mp1m = mp1p - 2;     // m-p-1
+
         if ( dimx == 1 ) {
                 double cove = V._D() ( 0 ) / mp1m ;
@@ -240 +240 @@
                 vec var ( l );
                 var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
-                var ( l - 1 ) = cove * cove / (mp1m-2);
+                var ( l - 1 ) = cove * cove / ( mp1m - 2 );
                 return var;
         } else {
-                ldmat Vll( V, linspace(0,dimx-1)); // top-left part of V
-                mat Y=Vll.to_mat();
-                mat varY(Y.rows(), Y.cols()); 
-                
-                double denom = (mp1p-1)*mp1m*mp1m*(mp1m-2);         // (m-p)(m-p-1)^2(m-p-3)
-                
-                int i,j;
-                for ( i=0; i<Y.rows(); i++){
-                        for ( j=0; j<Y.cols(); j++){
-                                varY(i,j) = (mp1p*Y(i,j)*Y(i,j) + mp1m * Y(i,i)* Y(j,j)) /denom;
+                ldmat Vll ( V, linspace ( 0, dimx - 1 ) ); // top-left part of V
+                mat Y = Vll.to_mat();
+                mat varY ( Y.rows(), Y.cols() );
+
+                double denom = ( mp1p - 1 ) * mp1m * mp1m * ( mp1m - 2 );         // (m-p)(m-p-1)^2(m-p-3)
+
+                int i, j;
+                for ( i = 0; i < Y.rows(); i++ ) {
+                        for ( j = 0; j < Y.cols(); j++ ) {
+                                varY ( i, j ) = ( mp1p * Y ( i, j ) * Y ( i, j ) + mp1m * Y ( i, i ) * Y ( j, j ) ) / denom;
                         }
                 }
-                vec mean_dR = diag(Y)/mp1m; // corresponds to cove
-                vec var_th=diag ( itmp.to_mat() );
-                vec var_Th ( mean_dR.length()*var_th.length() );
+                vec mean_dR = diag ( Y ) / mp1m; // corresponds to cove
+                vec var_th = diag ( itmp.to_mat() );
+                vec var_Th ( mean_dR.length() *var_th.length() );
                 // diagonal of diag(mean_dR) \kron diag(var_th)
-                for (int i=0; i<mean_dR.length(); i++){
-                        var_Th.set_subvector(i*var_th.length(), var_th*mean_dR(i));  
-                }
-                
-                return concat(var_Th, cvectorize(varY));
+                for ( int i = 0; i < mean_dR.length(); i++ ) {
+                        var_Th.set_subvector ( i*var_th.length(), var_th*mean_dR ( i ) );
+                }
+
+                return concat ( var_Th, cvectorize ( varY ) );
         }
 }