Changeset 730

Timestamp:

11/17/09 20:55:03 (14 years ago)

Author:

smidl

Message:

Reworked epdf_harness - testing is now more general.

Corrections of existing tests.

Files:

• applications/bdmtoolbox/mex/epdf_sample_mat.cpp (modified) (1 diff)
• library/bdm/stat/discrete.cpp (modified) (1 diff)
• library/bdm/stat/discrete.h (modified) (3 diffs)
• library/bdm/stat/exp_family.cpp (modified) (1 diff)
• library/bdm/stat/exp_family.h (modified) (1 diff)
• library/tests/epdf_harness.cpp (modified) (3 diffs)
• library/tests/epdf_harness.h (modified) (2 diffs)
• library/tests/testsuite/egiw.cfg (modified) (2 diffs)
• library/tests/testsuite/enorm.cfg (modified) (2 diffs)
• library/tests/testsuite/user_info_matrix.cfg (modified) (1 diff)

applications/bdmtoolbox/mex/epdf_sample_mat.cpp

@@ -27 +27 @@
 
         int nos=10;
-        if (n_input>1) {nos = (int)mexGetPr(input[1]);}
+        if (n_input>1) {nos = (int)(*mxGetPr(input[1]));}
         
         UImxArray Cfg(input[0]);

library/bdm/stat/discrete.cpp

@@ -122 +122 @@
 }
 
+void grid_fnc::set_values (const epdf &ep ) {
+        if ( sup.points() > 0 ) {
+                values ( 0 ) =  exp(ep.evallog( sup.first_vec() ));
+                for ( int j = 1; j < sup.points(); j++ ) {
+                        values ( j ) = exp(ep.evallog( sup.next_vec() ));
+                }
+        }
 }
+
+}

library/bdm/stat/discrete.h

@@ -62 +62 @@
         //! Get next active vector, call ONLY after first_vector()!
         const vec& next_vec();
-
+        
+        //! return active vector, call ONLY after first_vector()!
+        const vec& act_vec(){return actvec;};
+        
         //! \todo to je asi navic .. v predkovi!
         ivec nearest_point ( const vec &val );
@@ -140 +143 @@
         //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
         void set_values ( double ( *evalptr ) ( const vec& ) );
-
+        
+        //! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
+        void set_values ( const epdf &ep);
+        
         //! get value nearest to the given point
         double nearest_val ( const vec &val ) {
@@ -149 +155 @@
                 return vec_1 ( nearest_val ( val ) );
         }
+        const vec & _values() const { return values; }
 };
 UIREGISTER ( grid_fnc );

library/bdm/stat/exp_family.cpp

@@ -38 +38 @@
         
         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());
+        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
+        
+        mat OmChi;
+        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())));
+        }
+        return tmp;
 }
 

library/bdm/stat/exp_family.h

@@ -221 +221 @@
 
                 vec sample() const;
+                mat sample_mat(int n) const;
                 vec mean() const;
                 vec variance() const;

library/tests/epdf_harness.cpp

@@ -28 +28 @@
         UI::get ( variance, set, "variance", UI::compulsory );
 
-        UI::get ( support, set, "support", UI::optional );
-        UI::get ( nbins, set, "nbins", UI::optional );
+        support = UI::build<rectangular_support> ( set, "support", UI::optional );
         UI::get ( nsamples, set, "nsamples", UI::optional );
         UI::get ( R, set, "R", UI::optional );
@@ -44 +43 @@
         CHECK_CLOSE_EX ( variance, hepdf->variance(), tolerance );
 
-        if ( support.rows() == 2 ) {
-                int old_size = nbins.size();
-                if ( old_size < 2 ) {
-                        ivec new_nbins ( "100 100" );
-                        for ( int i = 0; i < old_size; ++i ) {
-                                new_nbins ( i ) = nbins ( i );
-                        }
-
-                        nbins = new_nbins;
-                }
-
-                check_support_mean();
-                check_support_integral();
+        if ( support ) { // support is given
+                grid_fnc ep_disc;
+                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);
+
+                vec pdf=ep_disc._values();
+                pdf /=sum(pdf); // normalize
+                
+                vec mea=pdf(0) * support->first_vec();
+                mat Remp=pdf(0) * outer_product(support->act_vec(), support->act_vec());
+
+                // run through all points
+                for (int i=1; i<support->points(); i++){
+                        mea += pdf(i)*support->next_vec();
+                        Remp += pdf(i) * outer_product(support->act_vec(), support->act_vec());
+                }
+                CHECK_CLOSE(mean, mea, tolerance);
+                CHECK_CLOSE(R, Remp-outer_product(mea,mea), tolerance);
         }
 
@@ -112 +120 @@
 }
 
-void epdf_harness::check_support_mean() {
-        vec xb = support.get_row ( 0 );
-        vec yb = support.get_row ( 1 );
-
-        vec actual;
-        actual = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
-
-        CHECK_CLOSE(mean, actual, tolerance);
-}
-
-void epdf_harness::check_support_integral() {
-        vec xb = support.get_row ( 0 );
-        vec yb = support.get_row ( 1 );
-
-        double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
-        CHECK_CLOSE(1.0,nc,0.01);
-}
 
 void epdf_harness::check_sample_mean() {

library/tests/epdf_harness.h

@@ -20 +20 @@
 #include "base/user_info.h"
 #include "stat/emix.h"
+#include <stat/discrete.h>
 
 namespace bdm {
@@ -28 +29 @@
         vec mean;
         vec variance;
-        mat support;
-        ivec nbins;
+        shared_ptr<rectangular_support> support;
         int nsamples;
         mat R;

library/tests/testsuite/egiw.cfg

@@ -14 +14 @@
   mean = [ 1.1, 0.2 ];
   variance = [ 0.02, 0.01333 ];
-  support = ( "matrix", 2, 2, [ -2.0, 4.0, 0.01, 2.0 ] );
+  R = ( "matrix", 2, 2, [ 0.02, 0.0, 0.0, 0.022 ] );
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ -2.0, 4.0 ], [0.01, 2.0 ] );
+        gridsizes = [ 100, 100 ];
+        };
   nbins = [ 100, 200 ];
   tolerance = 0.01;
@@ -33 +38 @@
   mean = [2.0, 1.14286];
   variance = [0.0285714, 0.0395795];
-  support = ( "matrix", 2, 2, [ 0.0, 4.0, 0.01, 3.0 ] );
-  nbins = [ 100, 200 ];
+  R = ( "matrix", 2, 2, [ 0.028, 0.0, 0.0, 0.078 ] );
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ 0.0, 4.0 ], [0.01, 3.0 ] );
+        gridsizes = [ 100, 100 ];
+        };
   tolerance = 0.01;
 } );

library/tests/testsuite/enorm.cfg

@@ -14 +14 @@
   mean = [ 1.1, -1.0 ];
   variance = [ 1, 2 ];
-  support = ( "matrix", 2, 2, [ -5.0, 5.0, -5.0, 5.0 ] );
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ -5.0, 5.0 ], [-5.0, 5.0 ] );
+        gridsizes = [ 100, 100 ];
+        };
   R = ( "matrix", 2, 2, [ 1.0, -0.5, -0.5, 2.0 ] );
   marginal_rv =     {
@@ -36 +40 @@
   mean = [ 0.0, 0.0 ];
   variance = [ 5.0, 5.2 ];
-  support = ( "matrix", 2, 2, [ -10.0, 10.0, -10.0, 10.0 ] );
+  support = {
+        class = "rectangular_support";
+        ranges = ( [ -10.0, 10.0 ], [-10.0, 10.0 ] );
+        gridsizes = [ 100, 100 ];
+        };
   R = ( "matrix", 2, 2, [ 5.0, -0.05, -0.05, 5.2 ] );
   nbins = [ 200, 200 ];

library/tests/testsuite/user_info_matrix.cfg

@@ -21 +21 @@
   manufacturer = "noname";
   electricLights = false;
-  matr = ( 0, 0, [ ] );
+  matr = ( "matrix", 0, 0, [ ] );
 };
 pepikovo :