diff --git a/doc/samples/vector_binary_redux.cpp b/doc/samples/vector_binary_redux.cpp
index abc36412..eae4488c 100644
--- a/doc/samples/vector_binary_redux.cpp
+++ b/doc/samples/vector_binary_redux.cpp
@@ -19,5 +19,6 @@ int main () {
         v1 (i) = v2 (i) = i;
 
     std::cout << inner_prod (v1, v2) << std::endl;
+    std::cout << covariance (v1, v2) << std::endl;
 }
 
diff --git a/doc/samples/vector_unary_redux.cpp b/doc/samples/vector_unary_redux.cpp
index 7de1918b..ec2684fc 100644
--- a/doc/samples/vector_unary_redux.cpp
+++ b/doc/samples/vector_unary_redux.cpp
@@ -19,6 +19,10 @@ int main () {
         v (i) = i;
 
     std::cout << sum (v) << std::endl;
+    std::cout << mean (v) << std::endl;
+    std::cout << mean_iterative (v) << std::endl;
+    std::cout << variance (v) << std::endl;
+    std::cout << variance_iterative (v) << std::endl;
     std::cout << norm_1 (v) << std::endl;
     std::cout << norm_2 (v) << std::endl;
     std::cout << norm_inf (v) << std::endl;
diff --git a/include/boost/numeric/ublas/detail/config.hpp b/include/boost/numeric/ublas/detail/config.hpp
index 9e674105..5ec91e10 100644
--- a/include/boost/numeric/ublas/detail/config.hpp
+++ b/include/boost/numeric/ublas/detail/config.hpp
@@ -58,6 +58,12 @@
 
 #endif
 
+// Oracle Developer Studio compiler
+#if defined (__SUNPRO_CC) && ! defined (BOOST_STRICT_CONFIG)
+#if __SUNPRO_CC >= 0x5130
+#define BOOST_UBLAS_USEFUL_ARRAY_PLACEMENT_NEW
+#endif
+#endif
 
 // GNU Compiler Collection
 #if defined (__GNUC__) && ! defined (BOOST_STRICT_CONFIG)
diff --git a/include/boost/numeric/ublas/functional.hpp b/include/boost/numeric/ublas/functional.hpp
index f54b441f..a03ce5d4 100644
--- a/include/boost/numeric/ublas/functional.hpp
+++ b/include/boost/numeric/ublas/functional.hpp
@@ -371,6 +371,202 @@ namespace boost { namespace numeric { namespace ublas {
         }
     };
 
+template<class V>
+    struct vector_mean: 
+        public vector_scalar_unary_functor<V> {
+        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_unary_functor<V>::result_type result_type;
+
+        template<class E>
+        static BOOST_UBLAS_INLINE
+        result_type apply (const vector_expression<E> &e) { 
+            result_type t = result_type (0);
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) 
+                t += e () (i);
+            return t / size;
+        }
+        // Dense case
+        template<class D, class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (D size, I it) { 
+            result_type t = result_type (0);
+            D i (0);
+            while (++ i <= size) {
+                t += *it;
+		++ it;
+            }
+            return t / size; 
+        }
+        // Sparse case
+        template<class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I it, const I &it_end) {
+            typedef typename I::difference_type vector_difference_type;
+            result_type t = result_type (0);
+            vector_difference_type size (0);
+            while (it != it_end) {
+                t += *it;
+                ++ it;
+                ++ size;
+            }
+            return t / size; 
+        }
+    };
+
+template<class V>
+    struct vector_mean_iterative: 
+        public vector_scalar_unary_functor<V> {
+        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_unary_functor<V>::result_type result_type;
+
+        template<class E>
+        static BOOST_UBLAS_INLINE
+        result_type apply (const vector_expression<E> &e) { 
+            result_type t = result_type (0);
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) 
+                t += (e () (i) - t) / (i + 1);
+            return t;
+        }
+        // Dense case
+        template<class D, class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (D size, I it) { 
+            result_type t = result_type (0);
+            D i (0);
+            while (++ i <= size) {
+                t += (*it - t) / i;
+		++ it;
+	    }
+            return t; 
+        }
+        // Sparse case
+        template<class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I it, const I &it_end) {
+            typedef typename I::difference_type vector_difference_type;
+            result_type t = result_type (0);
+            vector_difference_type size (1);
+            while (it != it_end) {
+                t += (*it - t) / size;
+                ++ it;
+                ++ size;
+            }
+            return t; 
+        }
+    };
+
+template<class V>
+    struct vector_variance: 
+        public vector_scalar_unary_functor<V> {
+        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_unary_functor<V>::result_type result_type;
+
+        template<class E>
+        static BOOST_UBLAS_INLINE 
+        result_type apply (const vector_expression<E> &e) { 
+            result_type sum = result_type (0);
+	    result_type sumsq = result_type (0);
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
+                sum += e () (i);
+                sumsq += e () (i) * e () (i);
+	    }
+            return (sumsq - (sum * sum) / size) / size ;
+        }
+        // Dense case
+        template<class D, class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (D size, I it) {
+            result_type sum = result_type (0);
+	    result_type sumsq = result_type (0);
+	    result_type n = result_type (0);
+            while (++ n <= size) {
+                sum += *it;
+                sumsq += *it * *it;
+                ++ it;
+            }
+            return (sumsq - (sum * sum) / size) / size ; 
+        }
+        // Sparse case
+        template<class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I it, const I &it_end) {
+            typedef typename I::difference_type vector_difference_type;
+            result_type sum = result_type (0);
+	    result_type sumsq = result_type (0);
+	    vector_difference_type n (0);
+            while (it != it_end) { 
+                ++ n;
+	        sum += *it;
+                sumsq += *it * *it;
+                ++ it;
+	    }  
+            return (sumsq - (sum * sum) / n) / n ;
+        }
+    };
+
+template<class V>
+    struct vector_variance_iterative: 
+        public vector_scalar_unary_functor<V> {
+        typedef typename vector_scalar_unary_functor<V>::value_type value_type;
+        typedef typename vector_scalar_unary_functor<V>::result_type result_type;
+
+        template<class E>
+        static BOOST_UBLAS_INLINE 
+        result_type apply (const vector_expression<E> &e) { 
+            result_type mean = result_type (0);
+	    result_type var = result_type (0);
+	    result_type del;
+            typedef typename E::size_type vector_size_type;
+            vector_size_type size (e ().size ());
+            for (vector_size_type i = 0; i < size; ++ i) {
+                del = e () (i) - mean;
+                mean += del / (i + 1);
+                var += del * (e () (i) - mean);
+	    }
+            return var / size;
+        }
+        // Dense case
+        template<class D, class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (D size, I it) { 
+            result_type mean = result_type (0);
+	    result_type n = result_type (0);
+	    result_type var = result_type (0);
+	    result_type del;
+            while (++ n <= size){
+                del = *it - mean;
+                mean += del / n;
+                var += del * (*it - mean);
+                ++ it;
+	    }
+            return var / size; 
+        }
+        // Sparse case
+        template<class I>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I it, const I &it_end) {
+            typedef typename I::difference_type vector_difference_type;
+            result_type mean = result_type (0);
+	    vector_difference_type n (0);
+	    result_type var = result_type (0);
+	    result_type del;
+            while (it != it_end)  {
+                ++ n;
+                del = *it - mean;
+                mean += del / n;
+                var += del * (*it - mean);
+                ++ it;
+	    }
+            return var / n; 
+        }
+    };
+
     // Unary returning real scalar 
     template<class V>
     struct vector_scalar_real_unary_functor {
@@ -814,6 +1010,202 @@ namespace boost { namespace numeric { namespace ublas {
         }
     };
 
+    template<class V1, class V2, class TV>
+    struct vector_covariance:
+        public vector_scalar_binary_functor<V1, V2, TV> {
+        typedef typename vector_scalar_binary_functor<V1, V2, TV>::value_type value_type;
+        typedef typename vector_scalar_binary_functor<V1, V2, TV>::result_type result_type;
+
+        template<class C1, class C2>
+        static BOOST_UBLAS_INLINE
+        result_type apply (const vector_container<C1> &c1,
+                           const vector_container<C2> &c2) {
+#ifdef BOOST_UBLAS_USE_SIMD
+            using namespace raw;
+            typedef typename C1::size_type vector_size_type;
+            vector_size_type size (BOOST_UBLAS_SAME (c1 ().size (), c2 ().size ()));
+            const typename V1::value_type *data1 = data_const (c1 ());
+            const typename V1::value_type *data2 = data_const (c2 ());
+            vector_size_type s1 = stride (c1 ());
+            vector_size_type s2 = stride (c2 ());
+            result_type t = result_type (0);
+	    result_type mean1 = result_type (0);
+            result_type mean2 = result_type (0);
+	    result_type del1, del2;
+            if (s1 == 1 && s2 == 1) {
+                for (vector_size_type i = 0; i < size; ++ i) {
+		    del1 = (data1[i] - mean1) / (i + 1);
+		    mean1 += del1;
+                    del2 = (data2[i] - mean2) / (i + 1);
+                    mean2 += del2;
+                    t += i * del1 * del2 - t / (i + 1);
+		}
+            } else if (s2 == 1) {
+                for (vector_size_type i = 0, i1 = 0; i < size; ++ i, i1 += s1) {
+		    del1 = (data1[i1] - mean1) / (i + 1);
+		    mean1 += del1;
+                    del2 = (data2[i] - mean2) / (i + 1);
+                    mean2 += del2;
+                    t += i * del1 * del2 - t / (i + 1);
+		}
+            } else if (s1 == 1) {
+                for (vector_size_type i = 0, i2 = 0; i < size; ++ i, i2 += s2) {
+		    del1 = (data1[i] - mean1) / (i + 1);
+		    mean1 += del1;
+                    del2 = (data2[i2] - mean2) / (i + 1);
+                    mean2 += del2;
+                    t += i * del1 * del2 - t / (i + 1);
+		}
+            } else {
+                for (vector_size_type i = 0, i1 = 0, i2 = 0; i < size; ++ i, i1 += s1, i2 += s2) {
+		    del1 = (data1[i1] - mean1) / (i + 1);
+		    mean1 += del1;
+                    del2 = (data2[i2] - mean2) / (i + 1);
+                    mean2 += del2;
+                    t += i * del1 * del2 - t / (i + 1);
+		}
+            }
+            return size / (size - 1) * t;
+#else
+            return apply (static_cast<const vector_expression<C1> > (c1), static_cast<const vector_expression<C2> > (c2));
+#endif
+        }
+        template<class E1, class E2>
+        static BOOST_UBLAS_INLINE
+        result_type apply (const vector_expression<E1> &e1,
+                           const vector_expression<E2> &e2) {
+            typedef typename E1::size_type vector_size_type;
+            vector_size_type size (BOOST_UBLAS_SAME (e1 ().size (), e2 ().size ()));
+            result_type t = result_type (0);
+	    result_type mean1 = result_type (0);
+            result_type mean2 = result_type (0);
+	    result_type del1, del2;
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+            for (vector_size_type i = 0; i < size; ++ i) {
+                del1 = (e1 () (i) - mean1) / (i + 1);
+                mean1 += del1;
+                del2 = (e2 () (i) - mean2) / (i + 1);
+                mean2 += del2;
+                t += i * del1 * del2 - t / (i + 1);
+	    }
+#else
+            vector_size_type i (0);
+            DD (size, 4, r, del1 = (e1 () (i) - mean1) / (i + 1),
+                mean1 += del1,
+                del2 = (e2 () (i) - mean2) / (i + 1),
+                mean2 += del2,
+                t += i * del1 * del2 - t / (i + 1), ++ i);
+#endif
+            return size / (size - 1) * t;
+        }
+        // Dense case
+        template<class D, class I1, class I2>
+        static BOOST_UBLAS_INLINE
+        result_type apply (D size, I1 it1, I2 it2) {
+            result_type t = result_type (0);
+            result_type i = result_type (0);
+	    result_type mean1 = result_type (0);
+            result_type mean2 = result_type (0);
+	    result_type del1, del2;
+#ifndef BOOST_UBLAS_USE_DUFF_DEVICE
+            while (++ i <= size) {
+                del1 = (*it1 - mean1) / i;
+                mean1 += del1;
+                del2 = (*it2 - mean2) / i;
+                mean2 += del2;
+                t += (i - 1) * del1 * del2 - t / i;
+                ++ it1;
+		++ it2;
+            }
+#else
+            DD (size, 4, r, del1 = (*it1 - mean1) / (i + 1),
+                mean1 += del1,
+                del2 = (*it2 - mean2) / (i + 1),
+                mean2 += del2,
+                t += i * del1 * del2 - t / (i + 1),
+                ++ it1, ++ it2, ++ i);
+#endif
+            return size / (size - 1) * t;
+        }
+        // Packed case
+        template<class I1, class I2>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end) {
+            result_type t = result_type (0);
+            result_type i = result_type (0);
+	    result_type mean1 = result_type (0);
+            result_type mean2 = result_type (0);
+	    result_type del1, del2;
+            typedef typename I1::difference_type vector_difference_type;
+            vector_difference_type it1_size (it1_end - it1);
+            vector_difference_type it2_size (it2_end - it2);
+            vector_difference_type diff (0);
+            if (it1_size > 0 && it2_size > 0)
+                diff = it2.index () - it1.index ();
+            if (diff != 0) {
+                vector_difference_type size = (std::min) (diff, it1_size);
+                if (size > 0) {
+                    it1 += size;
+                    it1_size -= size;
+                    diff -= size;
+                }
+                size = (std::min) (- diff, it2_size);
+                if (size > 0) {
+                    it2 += size;
+                    it2_size -= size;
+                    diff += size;
+                }
+            }
+            vector_difference_type size ((std::min) (it1_size, it2_size));
+            while (++ i <= size) {
+                del1 = (*it1 - mean1) / i;
+                mean1 += del1;
+                del2 = (*it2 - mean2) / i;
+                mean2 += del2;
+                t += (i - 1) * del1 * del2 - t / i;
+                ++ it1;
+		++ it2;
+            }
+            return size / (size - 1) * t;
+        }
+        // Sparse case
+        template<class I1, class I2>
+        static BOOST_UBLAS_INLINE
+        result_type apply (I1 it1, const I1 &it1_end, I2 it2, const I2 &it2_end, sparse_bidirectional_iterator_tag) {
+            typedef typename I1::difference_type vector_difference_type;
+            result_type t = result_type (0);
+            vector_difference_type i (1);
+	    result_type mean1 = result_type (0);
+            result_type mean2 = result_type (0);
+	    result_type del1, del2;
+            if (it1 != it1_end && it2 != it2_end) {
+                for (;;) {
+                    if (it1.index () == it2.index ()) {
+                        del1 = (*it1 - mean1) / i;
+                        mean1 += del1;
+                        del2 = (*it2 - mean2) / i;
+                        mean2 += del2;
+                        t += (i - 1) * del1 * del2 - t / i;
+                        ++ i;
+                        ++ it1;
+                        ++ it2;
+                        if (it1 == it1_end || it2 == it2_end)
+                            break;
+                    } else if (it1.index () < it2.index ()) {
+                        increment (it1, it1_end, it2.index () - it1.index ());
+                        if (it1 == it1_end)
+                            break;
+                    } else if (it1.index () > it2.index ()) {
+                        increment (it2, it2_end, it1.index () - it2.index ());
+                        if (it2 == it2_end)
+                            break;
+                    }
+                }
+            }
+            return i / (i - 1) * t;
+        }
+    };
+
     // Matrix functors
 
     // Binary returning vector
diff --git a/include/boost/numeric/ublas/vector_expression.hpp b/include/boost/numeric/ublas/vector_expression.hpp
index 6f952a46..afbba101 100644
--- a/include/boost/numeric/ublas/vector_expression.hpp
+++ b/include/boost/numeric/ublas/vector_expression.hpp
@@ -1589,6 +1589,42 @@ namespace boost { namespace numeric { namespace ublas {
         return expression_type (e ());
     }
 
+    // mean v - lazy mean
+    template<class E>
+    BOOST_UBLAS_INLINE
+    typename vector_scalar_unary_traits<E, vector_mean<E> >::result_type
+    mean (const vector_expression<E> &e) {
+        typedef typename vector_scalar_unary_traits<E, vector_mean<E> >::expression_type expression_type;
+        return expression_type (e ());
+    }
+
+    // mean v - iterative mean
+    template<class E>
+    BOOST_UBLAS_INLINE
+    typename vector_scalar_unary_traits<E, vector_mean_iterative<E> >::result_type
+    mean_iterative (const vector_expression<E> &e) {
+        typedef typename vector_scalar_unary_traits<E, vector_mean_iterative<E> >::expression_type expression_type;
+        return expression_type (e ());
+    }
+
+    // variance v = variance (v [i])
+    template<class E>
+    BOOST_UBLAS_INLINE
+    typename vector_scalar_unary_traits<E, vector_variance<E> >::result_type
+    variance (const vector_expression<E> &e) {
+        typedef typename vector_scalar_unary_traits<E, vector_variance<E> >::expression_type expression_type;
+        return expression_type (e ());
+    }
+
+    // variance v = variance (v [i])
+    template<class E>
+    BOOST_UBLAS_INLINE
+    typename vector_scalar_unary_traits<E, vector_variance_iterative<E> >::result_type
+    variance_iterative (const vector_expression<E> &e) {
+        typedef typename vector_scalar_unary_traits<E, vector_variance_iterative<E> >::expression_type expression_type;
+        return expression_type (e ());
+    }
+
     // real: norm_1 v = sum (abs (v [i]))
     // complex: norm_1 v = sum (abs (real (v [i])) + abs (imag (v [i])))
     template<class E>
@@ -1757,6 +1793,21 @@ namespace boost { namespace numeric { namespace ublas {
         return expression_type (e1 (), e2 ());
     }
 
+    // covariance (v1, v2)
+    template<class E1, class E2>
+    BOOST_UBLAS_INLINE
+    typename vector_scalar_binary_traits<E1, E2, vector_covariance<E1, E2,
+                                                                   typename promote_traits<typename E1::value_type,
+                                                                                           typename E2::value_type>::promote_type> >::result_type
+    covariance (const vector_expression<E1> &e1,
+                const vector_expression<E2> &e2) {
+        typedef typename vector_scalar_binary_traits<E1, E2, vector_covariance<E1, E2,
+                                                                   typename promote_traits<typename E1::value_type,
+                                                                                           typename E2::value_type>::promote_type> >::expression_type expression_type;
+        return expression_type (e1 (), e2 ());
+    }
+
+
 }}}
 
 #endif
diff --git a/test/test11.cpp b/test/test11.cpp
index 2d3d2752..0c82b896 100644
--- a/test/test11.cpp
+++ b/test/test11.cpp
@@ -119,6 +119,14 @@ struct test_my_vector {
         initialize_vector (v1);
         t = ublas::sum (v1);
         std::cout << "sum (v1) = " << t << std::endl;
+        t = ublas::mean (v1);
+        std::cout << "mean (v1) = " << t << std::endl;
+        t = ublas::mean_iterative (v1);
+        std::cout << "mean_iterative (v1) = " << t << std::endl;
+        t = ublas::variance (v1);
+        std::cout << "variance (v1) = " << t << std::endl;
+        t = ublas::variance_iterative (v1);
+        std::cout << "variance_iterative (v1) = " << t << std::endl;
         n = ublas::norm_1 (v1);
         std::cout << "norm_1 (v1) = " << n << std::endl;
         n = ublas::norm_2 (v1);
@@ -134,6 +142,8 @@ struct test_my_vector {
         initialize_vector (v2);
         t = ublas::inner_prod (v1, v2);
         std::cout << "inner_prod (v1, v2) = " << t << std::endl;
+        t = ublas::covariance (v1, v2);
+        std::cout << "covariance (v1, v2) = " << t << std::endl;
 
         // Scalar and Binary vector expression resulting in a vector
         initialize_vector (v1);
diff --git a/test/test31.cpp b/test/test31.cpp
index 77622872..5e5cdc59 100644
--- a/test/test31.cpp
+++ b/test/test31.cpp
@@ -105,6 +105,14 @@ struct test_my_vector {
             initialize_vector (v1);
             t = ublas::sum (v1);
             std::cout << "sum (v1) = " << t << std::endl;
+            t = ublas::mean (v1);
+            std::cout << "mean (v1) = " << t << std::endl;
+            t = ublas::mean_iterative (v1);
+            std::cout << "mean_iterative (v1) = " << t << std::endl;
+            t = ublas::variance (v1);
+            std::cout << "variance (v1) = " << t << std::endl;
+            t = ublas::variance_iterative (v1);
+            std::cout << "variance_iterative (v1) = " << t << std::endl;
             n = ublas::norm_1 (v1);
             std::cout << "norm_1 (v1) = " << n << std::endl;
             n = ublas::norm_2 (v1);
@@ -120,6 +128,8 @@ struct test_my_vector {
             initialize_vector (v2);
             t = ublas::inner_prod (v1, v2);
             std::cout << "inner_prod (v1, v2) = " << t << std::endl;
+            t = ublas::covariance (v1, v2);
+            std::cout << "covariance (v1, v2) = " << t << std::endl;
         }
     }
     void operator () () const {
diff --git a/test/test71.cpp b/test/test71.cpp
index 89e9ccf7..8f2ab346 100644
--- a/test/test71.cpp
+++ b/test/test71.cpp
@@ -89,6 +89,14 @@ struct test_my_vector {
             initialize_vector (v1);
             t = ublas::sum (v1);
             std::cout << "sum (v1) = " << t << std::endl;
+            t = ublas::mean (v1);
+            std::cout << "mean (v1) = " << t << std::endl;
+            t = ublas::mean_iterative (v1);
+            std::cout << "mean_iterative (v1) = " << t << std::endl;
+            t = ublas::variance (v1);
+            std::cout << "variance (v1) = " << t << std::endl; 
+            t = ublas::variance_iterative (v1);
+            std::cout << "variance_iterative (v1) = " << t << std::endl;
             n = ublas::norm_1 (v1);
             std::cout << "norm_1 (v1) = " << n << std::endl;
             n = ublas::norm_2 (v1);
@@ -104,6 +112,8 @@ struct test_my_vector {
             initialize_vector (v2);
             t = ublas::inner_prod (v1, v2);
             std::cout << "inner_prod (v1, v2) = " << t << std::endl;
+            t = ublas::covariance (v1, v2);
+            std::cout << "covariance (v1, v2) = " << t << std::endl;
         }
     }
     void operator () () const {