Changeset 1189 for library

Timestamp:

09/16/10 13:22:11 (14 years ago)

Author:

smidl

Message:

new tests of egiw

Location:

library

Files:

• bdm/estim/arx.h (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (2 diffs)
• bdm/stat/exp_family.h (modified) (9 diffs)
• tests/epdf_harness.cpp (modified) (2 diffs)
• tests/epdf_harness.h (modified) (1 diff)
• tests/testsuite/egiw.cfg (modified) (4 diffs)
• tests/testsuite/egiw_support_gen.m (added)
• tests/testsuite/epdf_test.cpp (modified) (1 diff)
• tests/testsuite/rectangular_support_test.cpp (modified) (1 diff)

library/bdm/estim/arx.h

@@ -221 +221 @@
         
         void bayes(const vec &dt, const vec &psi){
-          ldmat &Pbeta = est.P;
-          ldmat &Palpha = alternative.P;
+//        ldmat &Pbeta = est.P;
+//        ldmat &Palpha = alternative.P;
           
           

library/bdm/stat/exp_family.cpp

@@ -49 +49 @@
     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
+    eWishartCh Omega; //iWishart, same degrees of freedom as in Wishard in terms of delta
+        // delta = nu - npsi- p -1
+    Omega.set_parameters ( iChLam, nu - nPsi -dimx -1);
     Omega.validate();
 
-    mat OmChi;
+    chmat OmChi;
     mat Z ( M.rows(), M.cols() );
 
     mat Mi;
-    mat RChiT;
+    chmat Ri;
     mat tmp ( dimension(), n );
     M=M.T();// ugly hack == decide what to do with M.
     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() ) ) );
+        OmChi.inv(Ri );
+                if (M._datasize()>0){
+                        Z = randn ( M.rows(), M.cols() );
+                        Mi = M + Ri._Ch().T() * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
+                }
+        tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( Ri.to_mat() ) ) );
     }
     return tmp;
@@ -79 +81 @@
     factorize ( M, Vz, Lam );
 
-    chmat Ch;
-    Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) );
+    chmat Ch(Lam.to_mat());
     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
+    Omega.set_parameters ( iCh, nu - nPsi - dimx -1 ); // 2*nPsi is there to match numercial simulations - check if analytically correct
     Omega.validate();
 
     chmat Omi;
-    Omi.setCh ( Omega.sample_mat() );
-
+    Omi= Omega.sample_mat();
+        Omi.inv ( Ri );
+        
     if (M._datasize()>0) {
+                M = M.T();
         mat Z = randn ( M.rows(), M.cols() );
-        Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
-    }
-    Omi.inv ( Ri );
+        Mi = M + Ri._Ch().T() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+    }
 }
 

library/bdm/stat/exp_family.h

@@ -388 +388 @@
 * \brief Gauss-inverse-Wishart density stored in LD form
 
-* For \f$p\f$-variate densities, given rv.count() should be \f$p    imes\f$ V.rows().
+* For \f$p\f$-variate densities, 
+* 
+* \f$ M,R \sim GiW( V_t, \nu_t) \propto |R|^{0.5\nu}\exp(-1/2 tr(R^{-1}[I,M] V_t [I;M'])) \f$
 *
+* Factorizes as:
+* \f$ M|R \sim N( \hat{M}, R \otimes Vz^{-1}) \propto |R|^{-0.5dim(psi)}\exp(-1/2 tr((M-\hat{M})R^{-1}(M-\hat{M})Vz) \f$
+* \f$ R \sim iW( \Lambda,\delta) \propto |R|^{-0.5(\nu - dim(psi))} |\Lambda|^{}\exp(-1/2 tr(R^{-1}\Lambda) \f$
+* where in standard notation \f$ |R|^{-0.5(\delta + p +1)}\f$, i.e. 
+* \f$ \delta = \nu-dim(psi) -p-1 \f$
 */
 class egiw : public eEF {
@@ -426 +433 @@
 
     void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
-    //! LS estimate of \f$    heta\f$
+    //! LS estimate of \f$\theta\f$
     vec est_theta() const;
 
@@ -1612 +1619 @@
 
 /*! \brief Inverse Wishart density defined on Choleski decomposition
+ * 
+ * here \f$\Omega \sim |\Sigma|^{-0.5(\delta}|\Omega|^{0.5(\delta -p -1)} \exp(-1/2 tr(\Omega\Sigma^{-1}))\$f
 */
 class eWishartCh : public epdf {
 protected:
     //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
-    chmat Y;
+    chmat Sigma;
     //! dimension of matrix  \f$ \Psi \f$
     int p;
@@ -1624 +1633 @@
     //! Set internal structures
     void set_parameters ( const mat &Y0, const double delta0 ) {
-        Y = chmat ( Y0 );
+        Sigma = chmat ( Y0 );
         delta = delta0;
-        p = Y.rows();
+        p = Sigma.rows();
     }
     //! Set internal structures
     void set_parameters ( const chmat &Y0, const double delta0 ) {
-        Y = Y0;
+        Sigma = Y0;
         delta = delta0;
-        p = Y.rows();
+        p = Sigma.rows();
     }
 
@@ -1640 +1649 @@
     }
 
-    //! Sample matrix argument
-    mat sample_mat() const {
-        mat X = zeros ( p, p );
+    //! Sample matrix argument - lower 
+    //! Using Bartlet decomposition: W=L A A^T L^T where A is lower triag and L is choleski factor of Sigma.
+    chmat sample_mat() const {
+        mat A_T = zeros ( p, p ); // A transpose
 
         //sample diagonal
@@ -1648 +1658 @@
             GamRNG.setup ( 0.5* ( delta - i ) , 0.5 );   // no +1 !! index if from 0
 #pragma omp critical
-            X ( i, i ) = sqrt ( GamRNG() );
+            A_T ( i, i ) = sqrt ( GamRNG() );
         }
         //do the rest
@@ -1654 +1664 @@
             for ( int j = i + 1; j < p; j++ ) {
 #pragma omp critical
-                X ( i, j ) = NorRNG.sample();
+                A_T ( i, j ) = NorRNG.sample();
             }
         }
-        return X*Y._Ch();// return upper triangular part of the decomposition
+        chmat tmp;
+                tmp._Ch()=A_T*Sigma._Ch();
+        return tmp;// return upper triangular part of the decomposition
     }
 
     vec sample () const {
-        return vec ( sample_mat()._data(), p*p );
+        return vec ( sample_mat().to_mat()._data(), p*p );
     }
 
@@ -1673 +1685 @@
     //! fast access function y0 will be copied into Y.Ch.
     void setY ( const mat &Ch0 ) {
-        copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
+        copy_vector ( dim, Ch0._data(), Sigma._Ch()._data() );
     }
 
     //! fast access function y0 will be copied into Y.Ch.
     void _setY ( const vec &ch0 ) {
-        copy_vector ( dim, ch0._data(), Y._Ch()._data() );
+        copy_vector ( dim, ch0._data(), Sigma._Ch()._data() );
     }
 
     //! access function
     const chmat& getY() const {
-        return Y;
+        return Sigma;
     }
 };
@@ -1715 +1727 @@
     vec sample() const {
         mat iCh;
-        iCh = inv ( W.sample_mat() );
+        iCh = inv ( W.sample_mat()._Ch() );
         return vec ( iCh._data(), dim );
     }

library/tests/epdf_harness.cpp

@@ -28 +28 @@
     UI::get ( variance, set, "variance", UI::compulsory );
 
-    support = UI::build<rectangular_support> ( set, "support", UI::optional );
+    support = UI::build<support_base> ( set, "support", UI::optional );
     UI::get ( nsamples, set, "nsamples", UI::optional );
     UI::get ( R, set, "R", UI::optional );
@@ -45 +45 @@
     if ( support ) { // support is given
         grid_fnc ep_disc;
-        ep_disc.set_support ( *support );
+        ep_disc.set_support ( support );
         ep_disc.set_values ( *hepdf );
         // ep_disc is discretized at support points
 
-        double point_volume = prod ( support->_steps() );
-        CHECK_CLOSE ( 1.0, sum ( ep_disc._values() ) *point_volume, 0.01 );
+        CHECK_CLOSE ( 1.0,  ep_disc.integrate(), 0.01 );
 
         vec pdf = ep_disc._values();

library/tests/epdf_harness.h

@@ -29 +29 @@
     vec mean;
     vec variance;
-    shared_ptr<rectangular_support> support;
+    shared_ptr<support_base> support;
     int nsamples;
     mat R;

library/tests/testsuite/egiw.cfg

@@ -5 +5 @@
     class = "egiw";
     fV = ( "matrix", 2, 2, [ 13.1, 11.0, 11.0, 10.0 ] );
+// m=1.1; var_e = 0.1;
     nu = 10.0;
     dimx = 1;
@@ -28 +29 @@
     class = "egiw";
     fV = ( "matrix", 2, 2, [ 200.0, 80.0, 80.0, 40.0 ] );
+        // m=2; var_e = 1;
     nu = 40.0;
     dimx = 1;
@@ -45 +47 @@
         };
   tolerance = 0.01;
-});
-
-
-broken = 
+},
 {
   class = "epdf_harness";
   epdf = {
     class = "egiw";
-    V = ( "matrix", 2, 2, [ 10.0, 5.0, 5.0, 15.0 ] );
-    nu = 10.0;
+    fV = ( "matrix", 2, 2, [ 10.0, 5.0, 5.0, 15.0 ] );
+    nu = 16.0;
     dimx = 2;
     rv : 
@@ -62 +61 @@
     };
   };
-  mean = [1.0, 0.5, 0.5, 1.5];
-  variance = [1.0, 1.0];
-  R = ( "matrix", 2, 2, [ 1.0, 0.0, 0.0, 1.0 ] );
+  mean = [1.0, 0.5, 0.5, 1.5];// V/(delta-2-1)=V/(nu-2-1-2-1)
+  variance = [0.2500, 0.2045, 0.2045, 0.5625]; // m=13; var(i,j)=(12Vij^2+10Vii*Vjj)/(11*100*8)
+  R = ("matrix",4,4,[0.2500, 0.0, 0.0, 0.0, 0.0, 0.2045, 0.0, 0.0, 0.0, 0.0, 0.2045, 0.0, 0.0, 0.0, 0.0, 0.5625]); 
   support = {
-        class = "rectangular_support";
-        ranges = ( [ 0.0, 4.0 ], [-2.0, 2.0], [-2.0, 2.0], [0.0, 4.0 ] );
-        gridsizes = [ 10, 10, 10, 10];
+        class = "discrete_support";
+        points = "s@egiw_support.cfg";
+        volumes = "v@egiw_support.cfg";
         };
-  tolerance = 0.01;
-};
+  tolerance = 0.03;
+}
+);

library/tests/testsuite/epdf_test.cpp

@@ -53 +53 @@
     mat wM = "1.1 0.9; 0.9 1.0";
     eWishartCh eW;
-    eW.set_parameters ( wM / 100, 100 );
+    eW.set_parameters ( wM , 10 );
     eW.validate();
     mat mea = zeros ( 2, 2 );
-    mat Ch;
+    chmat Ch;
     for ( int i = 0; i < 100; i++ ) {
         Ch = eW.sample_mat();
-        mea += Ch.T() * Ch;
+        mea += Ch.to_mat();
     }
 
     mat actual = mea / 100;
-    CHECK_CLOSE ( wM, actual, 0.1 );
+    CHECK_CLOSE ( 10*wM, actual, 0.3 );
 }
 

library/tests/testsuite/rectangular_support_test.cpp

@@ -9 +9 @@
 TEST ( rectangular_support_test ) {
     rectangular_support rs;
-    CHECK_EQUAL ( rs.first_vec(), vec ( 0 ) );
+    //CHECK_EQUAL ( rs.first_vec(), vec ( 0 ) );
 
     Array<vec> range ( 2 );