Index: library/bdm/stat/discrete.cpp
===================================================================
--- library/bdm/stat/discrete.cpp (revision 723)
+++ library/bdm/stat/discrete.cpp (revision 730)
@@ -122,3 +122,12 @@
 }
 
+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() ));
+		}
+	}
 }
+
+}
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 728)
+++ library/bdm/stat/exp_family.h (revision 730)
@@ -221,4 +221,5 @@
 
 		vec sample() const;
+		mat sample_mat(int n) const;
 		vec mean() const;
 		vec variance() const;
Index: library/bdm/stat/discrete.h
===================================================================
--- library/bdm/stat/discrete.h (revision 706)
+++ library/bdm/stat/discrete.h (revision 730)
@@ -62,5 +62,8 @@
 	//! 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,5 +143,8 @@
 	//! 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,4 +155,5 @@
 		return vec_1 ( nearest_val ( val ) );
 	}
+	const vec & _values() const { return values; }
 };
 UIREGISTER ( grid_fnc );
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 725)
+++ library/bdm/stat/exp_family.cpp (revision 730)
@@ -38,4 +38,35 @@
 	
 	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;
 }
 
Index: library/tests/epdf_harness.cpp
===================================================================
--- library/tests/epdf_harness.cpp (revision 725)
+++ library/tests/epdf_harness.cpp (revision 730)
@@ -28,6 +28,5 @@
 	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,17 +43,26 @@
 	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,21 +120,4 @@
 }
 
-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() {
Index: library/tests/epdf_harness.h
===================================================================
--- library/tests/epdf_harness.h (revision 722)
+++ library/tests/epdf_harness.h (revision 730)
@@ -20,4 +20,5 @@
 #include "base/user_info.h"
 #include "stat/emix.h"
+#include <stat/discrete.h>
 
 namespace bdm {
@@ -28,6 +29,5 @@
 	vec mean;
 	vec variance;
-	mat support;
-	ivec nbins;
+	shared_ptr<rectangular_support> support;
 	int nsamples;
 	mat R;
Index: library/tests/testsuite/egiw.cfg
===================================================================
--- library/tests/testsuite/egiw.cfg (revision 725)
+++ library/tests/testsuite/egiw.cfg (revision 730)
@@ -14,5 +14,10 @@
   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,6 +38,10 @@
   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;
 } );
Index: library/tests/testsuite/user_info_matrix.cfg
===================================================================
--- library/tests/testsuite/user_info_matrix.cfg (revision 721)
+++ library/tests/testsuite/user_info_matrix.cfg (revision 730)
@@ -21,5 +21,5 @@
   manufacturer = "noname";
   electricLights = false;
-  matr = ( 0, 0, [ ] );
+  matr = ( "matrix", 0, 0, [ ] );
 };
 pepikovo : 
Index: library/tests/testsuite/enorm.cfg
===================================================================
--- library/tests/testsuite/enorm.cfg (revision 717)
+++ library/tests/testsuite/enorm.cfg (revision 730)
@@ -14,5 +14,9 @@
   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,5 +40,9 @@
   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 ];
