diff --git a/src/doCMASA.h b/src/doCMASA.h
index cdf446c8f..aa8f56927 100644
--- a/src/doCMASA.h
+++ b/src/doCMASA.h
@@ -33,6 +33,7 @@
 #include "doEstimator.h"
 #include "doModifierMass.h"
 #include "doSampler.h"
+#include "doStats.h"
 
 using namespace boost::numeric::ublas;
 
@@ -106,6 +107,7 @@ public:
 
 	eoPop< EOT > selected_pop;
 
+
 	//-------------------------------------------------------------
 	// Temporary solution used by plot to enumerate iterations in
 	// several files.
@@ -158,19 +160,16 @@ public:
 	    D distrib = _estimator(pop);
 
 	    double size = distrib.size();
-
 	    assert(size > 0);
 
-	    double vacc = 1;
+	    doHyperVolume hv;
 
 	    for (int i = 0; i < size; ++i)
 		{
-		    vacc *= sqrt(distrib.variance()[i]);
+		    hv.update( distrib.variance()[i] );
 		}
 
-	    assert(vacc <= std::numeric_limits<double>::max());
-
-	    ofs_params_var << vacc << std::endl;
+	    ofs_params_var << hv << std::endl;
 	}
 
 	do
@@ -297,30 +296,27 @@ public:
 		// dicreasement
 		//-------------------------------------------------------------
 
-		{
-		    // double size = distrib.size();
+		// {
+		//     double size = distrib.size();
+		//     assert(size > 0);
 
-		    // assert(size > 0);
+		//     vector< double > vmin(size);
+		//     vector< double > vmax(size);
 
-		    // vector< double > vmin(size);
-		    // vector< double > vmax(size);
+		//     std::copy(distrib.param(0).begin(), distrib.param(0).end(), vmin.begin());
+		//     std::copy(distrib.param(1).begin(), distrib.param(1).end(), vmax.begin());
 
-		    // std::copy(distrib.param(0).begin(), distrib.param(0).end(), vmin.begin());
-		    // std::copy(distrib.param(1).begin(), distrib.param(1).end(), vmax.begin());
+		//     vector< double > vrange = vmax - vmin;
 
-		    // vector< double > vrange = vmax - vmin;
+		//     doHyperVolume hv;
 
-		    // double vacc = 1;
+		//     for (int i = 0, size = vrange.size(); i < size; ++i)
+		//     	{
+		// 	    hv.update( vrange(i) );
+		//     	}
 
-		    // for (int i = 0, size = vrange.size(); i < size; ++i)
-		    // 	{
-		    // 	    vacc *= vrange(i);
-		    // 	}
-
-		    // assert(vacc <= std::numeric_limits<double>::max());
-
-		    // ofs_params << vacc << std::endl;
-		}
+		//     ofs_params << hv << std::endl;
+		// }
 
 		//-------------------------------------------------------------
 
@@ -332,19 +328,16 @@ public:
 
 		{
 		    double size = distrib.size();
-
 		    assert(size > 0);
 
-		    double vacc = 1;
+		    doHyperVolume hv;
 
 		    for (int i = 0; i < size; ++i)
 			{
-			    vacc *= sqrt(distrib.variance()[i]);
+			    hv.update( distrib.variance()[i] );
 			}
 
-		    assert(vacc <= std::numeric_limits<double>::max());
-
-		    ofs_params_var << vacc << std::endl;
+		    ofs_params_var << hv << std::endl;
 		}
 
 		//-------------------------------------------------------------
diff --git a/src/doEstimatorNormal.h b/src/doEstimatorNormal.h
index 133b98032..84430f350 100644
--- a/src/doEstimatorNormal.h
+++ b/src/doEstimatorNormal.h
@@ -19,14 +19,16 @@ public:
 	unsigned int dimsize = pop[0].size();
 	assert(dimsize > 0);
 
-	std::vector< Var > var(dimsize);
+	//std::vector< doVar > var(dimsize);
+	doCovMatrix cov(dimsize);
 
 	for (unsigned int i = 0; i < popsize; ++i)
 	    {
-		for (unsigned int d = 0; d < dimsize; ++d)
-		    {
-			var[d].update(pop[i][d]);
-		    }
+		cov.update(pop[i]);
+		// for (unsigned int d = 0; d < dimsize; ++d)
+		//     {
+		// 	var[d].update(pop[i][d]);
+		//     }
 	    }
 
 	EOT mean(dimsize);
@@ -34,9 +36,12 @@ public:
 
 	for (unsigned int d = 0; d < dimsize; ++d)
 	    {
-		mean[d] = var[d].get_mean();
-		variance[d] = var[d].get_var();
+		// mean[d] = var[d].get_mean();
+		// variance[d] = var[d].get_var();
 		//variance[d] = var[d].get_std(); // perhaps I should use this !?!
+
+		mean[d] = cov.get_mean(d);
+		variance[d] = cov.get_var(d);
 	    }
 
 	return doNormal< EOT >(mean, variance);
diff --git a/src/doStats.h b/src/doStats.h
index c46a72932..27a4b9c48 100644
--- a/src/doStats.h
+++ b/src/doStats.h
@@ -18,109 +18,96 @@
 #ifndef _doStats_h
 #define _doStats_h
 
-#include <cassert>
-#include <vector>
+#include <eoPrintable.h>
 
-class Mean {
+class doStats : public eoPrintable
+{
+public:
+    doStats();
 
-    double n;
-    double mean;
+    virtual void printOn(std::ostream&) const;
 
-    public:
-    Mean() : n(0), mean(0) {}
-
-    void update(double v) {
-	n++;
-	double d = v - mean;
-	mean += 1/n * d;
-    }
-
-    double get_mean() const { return mean; }
+protected:
+    double _n;
 };
 
-class Var {
-    double n;
-    double mean;
-    double sumvar;
+class doMean : public doStats
+{
+public:
+    doMean();
 
-    public:
-    Var() : n(0), mean(0), sumvar(0) {}
+    virtual void update(double);
+    virtual void printOn(std::ostream&) const;
 
-    void update(double v) {
-	n++;
-	double d = v - mean;
-	mean += 1/n * d;
-	sumvar += (n-1)/n * d * d;
-    }
+    double get_mean() const;
 
-    double get_mean() const { return mean;           }
-    double get_var()  const { return sumvar / (n-1); }
-    double get_std()  const { return sqrt(get_var());  }
+protected:
+    double _mean;
+};
+
+class doVar : public doMean
+{
+public:
+    doVar();
+
+    virtual void update(double);
+    virtual void printOn(std::ostream&) const;
+
+    double get_var() const;
+    double get_std() const;
+
+protected:
+    double _sumvar;
 };
 
 /** Single covariance between two variates */
-class Cov {
-    double n;
-    double meana;
-    double meanb;
-    double sumcov;
+class doCov : public doStats
+{
+public:
+    doCov();
 
-    public:
-    Cov() : n(0), meana(0), meanb(0), sumcov(0) {}
+    virtual void update(double, double);
+    virtual void printOn(std::ostream&) const;
 
-    void update(double a, double b) {
-	++n;
-	double da = a - meana;
-	double db = b - meanb;
+    double get_meana() const;
+    double get_meanb() const;
+    double get_cov()   const;
 
-	meana += 1/n * da;
-	meanb += 1/n * db;
-
-	sumcov += (n-1)/n * da * db;
-    }
-
-    double get_meana() const { return meana; }
-    double get_meanb() const { return meanb; }
-    double get_cov()   const { return sumcov / (n-1); }
+protected:
+    double _meana;
+    double _meanb;
+    double _sumcov;
 };
 
-class CovMatrix {
-    double n;
-    std::vector<double> mean;
-    std::vector< std::vector<double> > sumcov;
+class doCovMatrix : public doStats
+{
+public:
+    doCovMatrix(unsigned dim);
 
-    public:
-    CovMatrix(unsigned dim) : n(0), mean(dim), sumcov(dim , std::vector<double>(dim)) {}
+    virtual void update(const std::vector<double>&);
 
-    void update(const std::vector<double>& v) {
-	assert(v.size() == mean.size());
+    double get_mean(int) const;
+    double get_var(int) const;
+    double get_std(int) const;
+    double get_cov(int, int) const;
 
-	n++;
+protected:
+    std::vector< double > _mean;
+    std::vector< std::vector< double > > _sumcov;
+};
 
-	for (unsigned i = 0; i < v.size(); ++i) {
-	    double d = v[i] - mean[i];
-	    mean[i] += 1/n * d;
+class doHyperVolume : public doStats
+{
+public:
+    doHyperVolume();
 
-	    sumcov[i][i] += (n-1)/n * d * d;
+    virtual void update(double);
+    virtual void printOn(std::ostream&) const;
 
-	    for (unsigned j = i; j < v.size(); ++j) {
-		double e = v[j] - mean[j]; // mean[j] is not updated yet
-
-		double upd = (n-1)/n * d * e;
-
-		sumcov[i][j] += upd;
-		sumcov[j][i] += upd;
-
-	    }
-	}
-
-    }
-
-    double get_mean(int i) const { return mean[i]; }
-    double get_var(int i ) const { return sumcov[i][i] / (n-1); }
-    double get_std(int i) const  { return sqrt(get_var(i)); }
-    double get_cov(int i, int j) const { return sumcov[i][j] / (n-1); }
+    double get_hypervolume() const;
 
+protected:
+    double _hv;
 };
 
 #endif // !_doStats_h
diff --git a/src/test.cpp b/src/test.cpp
deleted file mode 100644
index 97e5ab2cc..000000000
--- a/src/test.cpp
+++ /dev/null
@@ -1,4 +0,0 @@
-namespace DO
-{
-    void	test(){}
-}