diff --git a/include/boost/math/tools/univariate_statistics.hpp b/include/boost/math/tools/univariate_statistics.hpp
index f88c8a724..32e3a9bb3 100644
--- a/include/boost/math/tools/univariate_statistics.hpp
+++ b/include/boost/math/tools/univariate_statistics.hpp
@@ -148,9 +148,9 @@ auto skewness(ForwardIterator first, ForwardIterator last)
         {
             Real delta21 = *it - M1;
             Real tmp = delta21/n;
-            M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
-            M2 = M2 + tmp*(n-1)*delta21;
-            M1 = M1 + tmp;
+            M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
+            M2 += tmp*(n-1)*delta21;
+            M1 += tmp;
             n += 1;
         }
 
diff --git a/test/univariate_statistics_test.cpp b/test/univariate_statistics_test.cpp
index 6bcbc5fe0..4896f52cc 100644
--- a/test/univariate_statistics_test.cpp
+++ b/test/univariate_statistics_test.cpp
@@ -228,6 +228,23 @@ void test_variance()
     }
     Real m2 = boost::math::tools::variance(v);
     BOOST_TEST(abs(m1 - m2) < tol*abs(m1));
+
+    // Wikipedia example for a variance of N sided die:
+    // https://en.wikipedia.org/wiki/Variance
+    for (size_t j = 16; j < 2048; j *= 2)
+    {
+        v.resize(1024);
+        Real n = v.size();
+        for (size_t i = 0; i < v.size(); ++i)
+        {
+            v[i] = i + 1;
+        }
+
+        sigma_sq = boost::math::tools::variance(v);
+
+        BOOST_TEST(abs(sigma_sq - (n*n-1)/Real(12)) <= tol*sigma_sq);
+    }
+
 }
 
 template<class Z>