diff --git a/examples/Makefile.am b/examples/Makefile.am
index 123bafb..8a51549 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -20,7 +20,8 @@ noinst_PROGRAMS =           \
     exPValue                \
     exRand                  \
     exRootFinder            \
-    exThreadPool
+    exThreadPool            \
+    exVarBenchmark
 
 exCompiledDoubleFunction_SOURCES  = exCompiledDoubleFunction.cpp
 exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
@@ -78,4 +79,8 @@ exThreadPool_SOURCES              = exThreadPool.cpp
 exThreadPool_CXXFLAGS             = $(COM_CXXFLAGS)
 exThreadPool_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
+exVarBenchmark_SOURCES            = exVarBenchmark.cpp
+exVarBenchmark_CXXFLAGS           = $(COM_CXXFLAGS)
+exVarBenchmark_LDFLAGS            = -L../lib/.libs -lLatAnalyze
+
 ACLOCAL_AMFLAGS = -I .buildutils/m4
diff --git a/examples/exVarBenchmark.cpp b/examples/exVarBenchmark.cpp
new file mode 100644
index 0000000..4c0ea9c
--- /dev/null
+++ b/examples/exVarBenchmark.cpp
@@ -0,0 +1,119 @@
+#include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Statistics/Random.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+
+using namespace std;
+using namespace Latan;
+
+constexpr Index  size          = 1000;
+constexpr Index  nSample       = 2000;
+
+int main(void)
+{
+    random_device rd;
+    DMat          var(size, size);
+    DVec          mean(size);
+    DMatSample    sample(nSample, size, 1), sample2(nSample, size, 1);
+
+    cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
+    var   = DMat::Random(size, size);
+    var  *= var.adjoint();
+    mean  = DVec::Random(size); 
+    RandomNormal mgauss(mean, var, rd());
+    sample[central] = mgauss();
+    FOR_STAT_ARRAY(sample, s)
+    {
+        sample[s] = mgauss();
+    }
+    sample2[central] = mgauss();
+    FOR_STAT_ARRAY(sample, s)
+    {
+        sample2[s] = mgauss();
+    }
+
+    cout << "-- check new routines" << endl;
+
+    DMat v, vo;
+
+    cout << "var" << endl;
+    auto start = chrono::high_resolution_clock::now();
+    vo = sample.varianceOld();
+    auto end = chrono::high_resolution_clock::now();
+    chrono::duration<double> diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.variance();
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    cout << "cov" << endl;
+    start = chrono::high_resolution_clock::now();
+    vo = sample.covarianceOld(sample2);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.covariance(sample2);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    cout << "mean" << endl;
+    start = chrono::high_resolution_clock::now();
+    vo = sample.meanOld(3, 5);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.mean(3, 5);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    cout << "varmat" << endl;
+    start = chrono::high_resolution_clock::now();
+    vo = sample.varianceMatrixOld();
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.varianceMatrix();
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    cout << "covarmat" << endl;
+    start = chrono::high_resolution_clock::now();
+    vo = sample.covarianceMatrixOld(sample2);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.covarianceMatrix(sample2);
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    cout << "corrmat" << endl;
+    start = chrono::high_resolution_clock::now();
+    vo = sample.correlationMatrixOld();
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    start = chrono::high_resolution_clock::now();
+    v = sample.correlationMatrix();
+    end = chrono::high_resolution_clock::now();
+    diff = end - start;
+    cout << "time " << diff.count() << endl;
+    cout << "diff " << (v - vo).norm() << endl;
+
+    return EXIT_SUCCESS;
+}
diff --git a/lib/Statistics/MatSample.hpp b/lib/Statistics/MatSample.hpp
index a48db9d..a103215 100644
--- a/lib/Statistics/MatSample.hpp
+++ b/lib/Statistics/MatSample.hpp
@@ -103,6 +103,10 @@ public:
                      const Index nCol);
     // resize all matrices
     void resizeMat(const Index nRow, const Index nCol);
+    // covariance matrix
+    Mat<T> covarianceMatrix(const MatSample<T> &sample) const;
+    Mat<T> varianceMatrix(void) const;
+    Mat<T> correlationMatrix(void) const;
 };
 
 // non-member operators
@@ -379,6 +383,78 @@ void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
     }
 }
 
+// covariance matrix ///////////////////////////////////////////////////////////
+template <typename T>
+Mat<T> MatSample<T>::covarianceMatrix(const MatSample<T> &sample) const
+{
+    if (((*this)[central].cols() != 1) or (sample[central].cols() != 1))
+    {
+        LATAN_ERROR(Size, "samples have more than one column");
+    }
+
+    Index  n1 = (*this)[central].rows(), n2 = sample[central].rows();
+    Index  nSample = this->size();
+    Mat<T> tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2);
+    Mat<T> s1(n1, 1), s2(n2, 1), one(nSample, 1);
+
+    one.fill(1.);
+    s1.fill(0.);
+    s2.fill(0.);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s1          += (*this)[s];
+        tmp1.col(s)  = (*this)[s];
+    }
+    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s2          += sample[s];
+        tmp2.col(s)  = sample[s];
+    }
+    tmp2 -= s2*one.transpose()/static_cast<double>(nSample);
+    res   = tmp1*tmp2.transpose()/static_cast<double>(nSample - 1);
+
+    return res;
+}
+
+template <typename T>
+Mat<T> MatSample<T>::varianceMatrix(void) const
+{
+    if ((*this)[central].cols() != 1)
+    {
+        LATAN_ERROR(Size, "samples have more than one column");
+    }
+
+    Index  n1 = (*this)[central].rows();
+    Index  nSample = this->size();
+    Mat<T> tmp1(n1, nSample), res(n1, n1);
+    Mat<T> s1(n1, 1), one(nSample, 1);
+
+    one.fill(1.);
+    s1.fill(0.);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s1          += (*this)[s];
+        tmp1.col(s)  = (*this)[s];
+    }
+    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
+    res   = tmp1*tmp1.transpose()/static_cast<double>(nSample - 1);
+
+    return res;
+}
+
+template <typename T>
+Mat<T> MatSample<T>::correlationMatrix(void) const
+{
+    Mat<T> res = varianceMatrix();
+    Mat<T> invDiag(res.rows(), 1);
+
+    invDiag = res.diagonal();
+    invDiag = invDiag.cwiseInverse().cwiseSqrt();
+    res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
+
+    return res;
+}
 
 END_LATAN_NAMESPACE
 
diff --git a/lib/Statistics/StatArray.hpp b/lib/Statistics/StatArray.hpp
index 9637089..15179e1 100644
--- a/lib/Statistics/StatArray.hpp
+++ b/lib/Statistics/StatArray.hpp
@@ -51,14 +51,17 @@ public:
     const T & operator[](const Index s) const;
     // statistics
     void bin(Index binSize);
+    T    sum(const Index pos = 0, const Index n = -1) const;
+    T    meanOld(const Index pos = 0, const Index n = -1) const;
     T    mean(const Index pos = 0, const Index n = -1) const;
-    T    covariance(const StatArray<T, os> &array, const Index pos = 0,
-                    const Index n = -1) const;
-    T    covarianceMatrix(const StatArray<T, os> &array, const Index pos = 0,
+    T    covarianceOld(const StatArray<T, os> &array) const;
+    T    covariance(const StatArray<T, os> &array) const;
+    T    covarianceMatrixOld(const StatArray<T, os> &array, const Index pos = 0,
                           const Index n = -1) const;
-    T    variance(const Index pos = 0, const Index n = -1) const;
-    T    varianceMatrix(const Index pos = 0, const Index n = -1) const;
-    T    correlationMatrix(const Index pos = 0, const Index n = -1) const;
+    T    varianceOld(void) const;
+    T    variance(void) const;
+    T    varianceMatrixOld(const Index pos = 0, const Index n = -1) const;
+    T    correlationMatrixOld(const Index pos = 0, const Index n = -1) const;
     // IO type
     virtual IoType getType(void) const;
 public:
@@ -66,7 +69,7 @@ public:
 };
 
 // reduction operations
-namespace ReducOp
+namespace StatOp
 {
     // general templates
     template <typename T>
@@ -149,42 +152,82 @@ void StatArray<T, os>::bin(Index binSize)
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::mean(const Index pos, const Index n) const
+T StatArray<T, os>::meanOld(const Index pos, const Index n) const
 {
     T           result = T();
     const Index m = (n >= 0) ? n : size();
 
     if (m)
     {
-        result = this->segment(pos+os, m).redux(&ReducOp::sum<T>);
+        result = this->segment(pos+os, m).redux(&StatOp::sum<T>);
     }
     return result/static_cast<double>(m);
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::covariance(const StatArray<T, os> &array, const Index pos,
-                               const Index n) const
+T StatArray<T, os>::sum(const Index pos, const Index n) const
+{
+    T           result;
+    const Index m = (n >= 0) ? n : size();
+
+    result = (*this)[pos];
+    for (Index i = pos + 1; i < pos + m; ++i)
+    {
+        result += (*this)[i];
+    }
+
+    return result;
+}
+
+template <typename T, Index os>
+T StatArray<T, os>::mean(const Index pos, const Index n) const
+{
+    const Index m = (n >= 0) ? n : size();
+
+    return sum(pos, n)/static_cast<double>(m);
+}
+
+template <typename T, Index os>
+T StatArray<T, os>::covarianceOld(const StatArray<T, os> &array) const
 {
     T           s1, s2, prs, res = T();
-    const Index m = (n >= 0) ? n : size();
+    const Index m = size();
     
     if (m)
     {
-        auto arraySeg = array.segment(pos+os, m);
-        auto thisSeg  = this->segment(pos+os, m);
+        auto arraySeg = array.segment(os, m);
+        auto thisSeg  = this->segment(os, m);
         
-        s1  = thisSeg.redux(&ReducOp::sum<T>);
-        s2  = arraySeg.redux(&ReducOp::sum<T>);
-        prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod<T>)
-                     .redux(&ReducOp::sum<T>);
-        res = prs - ReducOp::prod(s1, s2)/static_cast<double>(m);
+        s1  = thisSeg.redux(&StatOp::sum<T>);
+        s2  = arraySeg.redux(&StatOp::sum<T>);
+        prs = thisSeg.binaryExpr(arraySeg, &StatOp::prod<T>)
+                     .redux(&StatOp::sum<T>);
+        res = prs - StatOp::prod(s1, s2)/static_cast<double>(m);
     }
     
     return res/static_cast<double>(m - 1);
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
+T StatArray<T, os>::covariance(const StatArray<T, os> &array) const
+{
+    T s1, s2, res;
+
+    s1  = array.sum();
+    s2  = this->sum();
+    res = StatOp::prod<T>(array[0], (*this)[0]);
+    for (Index i = 1; i < size(); ++i)
+    {
+        res += StatOp::prod<T>(array[i], (*this)[i]);
+    }
+    res -= StatOp::prod<T>(s1, s2)/static_cast<double>(size());
+    res /= static_cast<double>(size() - 1);
+
+    return res;
+}
+
+template <typename T, Index os>
+T StatArray<T, os>::covarianceMatrixOld(const StatArray<T, os> &array,
                                      const Index pos, const Index n) const
 {
     T           s1, s2, prs, res = T();
@@ -195,32 +238,38 @@ T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
         auto arraySeg = array.segment(pos+os, m);
         auto thisSeg  = this->segment(pos+os, m);
 
-        s1  = thisSeg.redux(&ReducOp::sum<T>);
-        s2  = arraySeg.redux(&ReducOp::sum<T>);
-        prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd<T>)
-                     .redux(&ReducOp::sum<T>);
-        res = prs - ReducOp::tensProd(s1, s2)/static_cast<double>(m);
+        s1  = thisSeg.redux(&StatOp::sum<T>);
+        s2  = arraySeg.redux(&StatOp::sum<T>);
+        prs = thisSeg.binaryExpr(arraySeg, &StatOp::tensProd<T>)
+                     .redux(&StatOp::sum<T>);
+        res = prs - StatOp::tensProd(s1, s2)/static_cast<double>(m);
     }
     
     return res/static_cast<double>(m - 1);
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::variance(const Index pos, const Index n) const
+T StatArray<T, os>::variance(void) const
 {
-    return covariance(*this, pos, n);
+    return covariance(*this);
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) const
+T StatArray<T, os>::varianceOld(void) const
 {
-    return covarianceMatrix(*this, pos, n);
+    return covarianceOld(*this);
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
+T StatArray<T, os>::varianceMatrixOld(const Index pos, const Index n) const
 {
-    T res = varianceMatrix(pos, n);
+    return covarianceMatrixOld(*this, pos, n);
+}
+
+template <typename T, Index os>
+T StatArray<T, os>::correlationMatrixOld(const Index pos, const Index n) const
+{
+    T res = varianceMatrixOld(pos, n);
     T invDiag(res.rows(), 1);
 
     invDiag = res.diagonal();
@@ -231,8 +280,14 @@ T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
 }
 
 // reduction operations ////////////////////////////////////////////////////////
-namespace ReducOp
+namespace StatOp
 {
+    template <typename T>
+    inline void zero(T &a)
+    {
+        a = 0.;
+    }
+
     template <typename T>
     inline T sum(const T &a, const T &b)
     {