Index: /library/tests/mat_checks.h
===================================================================
--- /library/tests/mat_checks.h (revision 477)
+++ /library/tests/mat_checks.h (revision 481)
@@ -16,4 +16,5 @@
 #include "../bdm/itpp_ext.h"
 #include "UnitTest++.h"
+#include <string>
 
 namespace UnitTest {
@@ -52,8 +53,12 @@
 }
 
-/*! CHECK_EQUAL_EX and CHECK_CLOSE_EX macros should be used only in
-  blocks having an instance of this class (which sets the globals for
-  error reporting). */
+/*! CHECK_*_EX macros should be used only in blocks having an instance
+  of this class (which sets the globals for error reporting). */
 class CurrentContext {
+public:
+	// how many times to repeat a failing test before reporting
+	// failure
+	static const int max_trial_count = 3;
+
 private:
 	static const char *config_name;
@@ -65,4 +70,10 @@
 	~CurrentContext();
 
+	/* Should be called only in blocks having an instance of
+	   CurrentContext. The argument, when not default, should be
+	   __LINE__ (and it is included in the returned string).
+	*/
+	static std::string format_context( int ln = -1 );
+
 	template<typename Expected, typename Actual>
 	static void CheckEqualEx ( UnitTest::TestResults& results,
@@ -72,5 +83,5 @@
 		if ( ! ( expected == actual ) ) {
 			UnitTest::MemoryOutStream stream;
-			stream << "error at " << config_name << '[' << index << "]: expected " << expected << " but was " << actual;
+			stream << format_context() << "expected " << expected << " but was " << actual;
 
 			results.OnTestFailure ( details, stream.GetText() );
@@ -86,5 +97,5 @@
 		if ( !UnitTest::AreClose ( expected, actual, tolerance ) ) {
 			UnitTest::MemoryOutStream stream;
-			stream << "error at " << config_name << '[' << index << "]: expected " << expected << " +/- " << tolerance << " but was " << actual;
+			stream << format_context() << "expected " << expected << " +/- " << tolerance << " but was " << actual;
 
 			results.OnTestFailure ( details, stream.GetText() );
Index: /library/tests/mat_checks.cpp
===================================================================
--- /library/tests/mat_checks.cpp (revision 477)
+++ /library/tests/mat_checks.cpp (revision 481)
@@ -1,3 +1,4 @@
 #include "mat_checks.h"
+#include <sstream>
 
 namespace UnitTest {
@@ -66,2 +67,14 @@
 	index = -1;
 }
+
+std::string CurrentContext::format_context(int ln) {
+	std::stringstream ss;
+	ss << "error at " << config_name << '[' << index << ']';
+
+	if ( ln >= 0 ) {
+    		ss << ", harness line " << ln;
+	}
+
+	ss << ": ";
+	return ss.str();
+}
Index: /library/tests/epdf_harness.h
===================================================================
--- /library/tests/epdf_harness.h (revision 477)
+++ /library/tests/epdf_harness.h (revision 481)
@@ -21,4 +21,5 @@
 #include "base/bdmbase.h"
 #include "base/user_info.h"
+#include "stat/emix.h"
 
 namespace bdm {
@@ -53,4 +54,25 @@
 		return tolerance;
 	}
+
+private:
+	// assumes support is valid
+	void check_support_mean();
+
+	// assumes support is valid
+	void check_support_integral();
+
+	// assumes R is valid - it probably shouldn't, but calling it
+	// unconditionally leads to a lot of not implemented
+	// functions...
+	void check_sample_mean();
+
+	// assumes R is valid
+	void check_covariance();
+
+	// assumes marginal_rv is valid
+	void check_cond_mean( mprod &mep );
+
+	// assumes marginal_rv and R are valid
+	void check_cond_covariance( mprod &mep );
 };
 
Index: /library/tests/epdf_harness.cpp
===================================================================
--- /library/tests/epdf_harness.cpp (revision 480)
+++ /library/tests/epdf_harness.cpp (revision 481)
@@ -3,5 +3,4 @@
 #include "base/user_info.h"
 #include "stat/exp_family.h"
-#include "stat/emix.h"
 #include "mat_checks.h"
 #include "test_util.h"
@@ -50,7 +49,4 @@
 
 	if ( support.rows() == 2 ) {
-		vec xb = support.get_row ( 0 );
-		vec yb = support.get_row ( 1 );
-
 		int old_size = nbins.size();
 		if ( old_size < 2 ) {
@@ -63,18 +59,11 @@
 		}
 
-		CHECK_CLOSE_EX ( mean, num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) ), tolerance );
-		CHECK_CLOSE_EX ( 1.0, normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) ), tolerance );
+		check_support_mean();
+		check_support_integral();
 	}
 
 	if ( R.rows() > 0 ) {
-		mat smp = hepdf->sample_m ( nsamples );
-		vec emu = smp * ones ( nsamples ) / nsamples;
-		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
-
-		// simplify overloading for Visual Studio
-		vec delta = sqrt ( variance ) / sqrt ( static_cast<double> ( nsamples ) );
-		CHECK_CLOSE_EX ( mean, emu, delta );
-
-		CHECK_CLOSE_EX ( R, er, tolerance );
+		check_sample_mean();
+		check_covariance();
 	}
 
@@ -90,14 +79,7 @@
 		mprod mEp ( aa );
 
-		mat smp = mEp.samplecond ( vec ( 0 ), nsamples );
-		vec emu = sum ( smp, 2 ) / nsamples;
-
-		// simplify overloading for Visual Studio
-		vec delta = sqrt ( variance ) / sqrt ( static_cast<double> ( nsamples ) );
-		CHECK_CLOSE_EX ( mean, emu, delta );
-
+		check_cond_mean(mEp);
 		if ( R.rows() > 0 ) {
-			mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
-			CHECK_CLOSE_EX ( R, er, tolerance );
+			check_cond_covariance( mEp );
 		}
 
@@ -113,3 +95,142 @@
 }
 
-}
+void epdf_harness::check_support_mean() {
+	vec xb = support.get_row ( 0 );
+	vec yb = support.get_row ( 1 );
+
+	int tc = 0;
+	Array<vec> actual(CurrentContext::max_trial_count);
+	do {
+		vec emu = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
+		actual( tc ) = emu;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+void epdf_harness::check_support_integral() {
+	vec xb = support.get_row ( 0 );
+	vec yb = support.get_row ( 1 );
+
+	int tc = 0;
+	Array<double> actual(CurrentContext::max_trial_count);
+	do {
+		double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
+		actual( tc ) = nc;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( 1.0, actual( tc - 1 ), tolerance ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( 1.0, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+void epdf_harness::check_sample_mean() {
+	// simplify overloading for Visual Studio
+	vec delta = 2 * ( sqrt ( variance ) / sqrt ( static_cast<double> ( nsamples ) ) );
+
+	int tc = 0;
+	Array<vec> actual(CurrentContext::max_trial_count);
+	do {
+		mat smp = hepdf->sample_m ( nsamples );
+		vec emu = smp * ones ( nsamples ) / nsamples;
+		actual( tc ) = emu;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+
+void epdf_harness::check_covariance() {
+	int tc = 0;
+	Array<mat> actual(CurrentContext::max_trial_count);
+	do {
+		mat smp = hepdf->sample_m ( nsamples );
+		vec emu = smp * ones ( nsamples ) / nsamples;
+		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
+		actual( tc ) = er;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << R << " +/- " << tolerance << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+void epdf_harness::check_cond_mean( mprod &mep ) {
+	// simplify overloading for Visual Studio
+	vec delta = 2 * ( sqrt ( variance ) / sqrt ( static_cast<double> ( nsamples ) ) );
+
+	int tc = 0;
+	Array<vec> actual(CurrentContext::max_trial_count);
+	do {
+		mat smp = mep.samplecond ( vec ( 0 ), nsamples );
+		vec emu = sum ( smp, 2 ) / nsamples;
+		actual( tc ) = emu;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+void epdf_harness::check_cond_covariance( mprod &mep ) {
+	int tc = 0;
+	Array<mat> actual(CurrentContext::max_trial_count);
+	do {
+		mat smp = mep.samplecond ( vec ( 0 ), nsamples );
+		vec emu = sum ( smp, 2 ) / nsamples;
+		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
+		actual( tc ) = er;
+		++tc;
+	} while ( ( tc < CurrentContext::max_trial_count ) &&
+		  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
+	if ( ( tc == CurrentContext::max_trial_count ) &&
+	     ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
+		UnitTest::MemoryOutStream stream;
+		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;
+
+		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);
+
+		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
+       }
+}
+
+}
