Statistics: Array and Streaming Variance

src/Numerics/Statistics/ArrayStatistics.cs

@@ -80,7 +80,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// Returns the arithmetic sample mean from the unsorted data array.
+        /// Estimates the arithmetic sample mean from the unsorted data array.
         /// Returns NaN if data is empty or any entry is NaN.
         /// </summary>
         /// <param name="data">Sample array, no sorting is assumed.</param>
@@ -97,5 +97,27 @@ namespace MathNet.Numerics.Statistics
             }
             return mean;
         }
+
+        /// <summary>
+        /// Estimates the unbiased population or sample variance from the unsorted data array.
+        /// On a dataset of size N will use an N-1 normalizer
+        /// Returns NaN if data is empty or any entry is NaN.
+        /// </summary>
+        /// <param name="data">Sample array, no sorting is assumed.</param>
+        public static double Variance(double[] data)
+        {
+            if (data == null) throw new ArgumentNullException("data");
+            if (data.Length <= 1) return double.NaN;
+
+            double variance = 0;
+            double t = data[0];
+            for (int i = 1; i < data.Length; i++)
+            {
+                t += data[i];
+                double diff = ((i + 1) * data[i]) - t;
+                variance += (diff * diff) / ((i + 1) * i);
+            }
+            return variance/(data.Length - 1);
+        }
     }
 }

src/Numerics/Statistics/Statistics.cs

@@ -118,34 +118,10 @@ namespace MathNet.Numerics.Statistics
         /// <returns>The unbiased population variance of the sample.</returns>
         public static double Variance(this IEnumerable<double> data)
         {
-            if (data == null)
-            {
-                throw new ArgumentNullException("data");
-            }
-
-            double variance = 0;
-            double t = 0;
-            ulong j = 0;
-
-            using (IEnumerator<double> iterator = data.GetEnumerator())
-            {
-                if (iterator.MoveNext())
-                {
-                    j++;
-                    t = iterator.Current;
-                }
-
-                while (iterator.MoveNext())
-                {
-                    j++;
-                    double xi = iterator.Current;
-                    t += xi;
-                    double diff = (j * xi) - t;
-                    variance += (diff * diff) / (j * (j - 1));
-                }
-            }
-
-            return j > 1 ? variance/(j - 1) : double.NaN;
+            var array = data as double[];
+            return array != null
+                ? ArrayStatistics.Variance(array)
+                : StreamingStatistics.Variance(data);
         }
 
         /// <summary>
@@ -155,47 +131,8 @@ namespace MathNet.Numerics.Statistics
         /// <returns>The population variance of the sample.</returns>
         public static double Variance(this IEnumerable<double?> data)
         {
-            if (data == null)
-            {
-                throw new ArgumentNullException("data");
-            }
-
-            double variance = 0;
-            double t = 0;
-            ulong j = 0;
-
-            using (IEnumerator<double?> iterator = data.GetEnumerator())
-            {
-                while (true)
-                {
-                    bool hasNext = iterator.MoveNext();
-                    if (!hasNext)
-                    {
-                        break;
-                    }
-
-                    if (iterator.Current.HasValue)
-                    {
-                        j++;
-                        t = iterator.Current.Value;
-                        break;
-                    }
-                }
-
-                while (iterator.MoveNext())
-                {
-                    if (iterator.Current.HasValue)
-                    {
-                        j++;
-                        double xi = iterator.Current.Value;
-                        t += xi;
-                        double diff = (j * xi) - t;
-                        variance += (diff * diff) / (j * (j - 1));
-                    }
-                }
-            }
-
-            return j > 1 ? variance/(j - 1) : double.NaN;
+            if (data == null) throw new ArgumentNullException("data");
+            return StreamingStatistics.Variance(data.Where(d => d.HasValue).Select(d => d.Value));
         }
 
         /// <summary>

src/Numerics/Statistics/StreamingStatistics.cs

@@ -80,7 +80,7 @@ namespace MathNet.Numerics.Statistics
         }
 
         /// <summary>
-        /// Returns the arithmetic sample mean from the enumerable, in a single pass without memoization.
+        /// Estimates the arithmetic sample mean from the enumerable, in a single pass without memoization.
         /// Returns NaN if data is empty or any entry is NaN.
         /// </summary>
         /// <param name="stream">Sample stream, no sorting is assumed.</param>
@@ -98,5 +98,38 @@ namespace MathNet.Numerics.Statistics
             }
             return any ? mean : double.NaN;
         }
+
+        /// <summary>
+        /// Estimates the unbiased population or sample variance from the enumerable, in a single pass without memoization.
+        /// On a dataset of size N will use an N-1 normalizer
+        /// Returns NaN if data is empty or any entry is NaN.
+        /// </summary>
+        /// <param name="stream">Sample stream, no sorting is assumed.</param>
+        public static double Variance(IEnumerable<double> stream)
+        {
+            if (stream == null) throw new ArgumentNullException("stream");
+
+            double variance = 0;
+            double t = 0;
+            ulong j = 0;
+            using (var iterator = stream.GetEnumerator())
+            {
+                if (iterator.MoveNext())
+                {
+                    j++;
+                    t = iterator.Current;
+                }
+
+                while (iterator.MoveNext())
+                {
+                    j++;
+                    double xi = iterator.Current;
+                    t += xi;
+                    double diff = (j * xi) - t;
+                    variance += (diff * diff) / (j * (j - 1));
+                }
+            }
+            return j > 1 ? variance / (j - 1) : double.NaN;
+        }
     }
 }

src/UnitTests/StatisticsTests/StatisticsTests.cs

@@ -83,10 +83,12 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Minimum(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Maximum(data));
             Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Mean(data));
+            Assert.Throws<ArgumentNullException>(() => ArrayStatistics.Variance(data));
 
             Assert.Throws<ArgumentNullException>(() => StreamingStatistics.Minimum(data));
             Assert.Throws<ArgumentNullException>(() => StreamingStatistics.Maximum(data));
             Assert.Throws<ArgumentNullException>(() => StreamingStatistics.Mean(data));
+            Assert.Throws<ArgumentNullException>(() => StreamingStatistics.Variance(data));
         }
 
         [Test]
@@ -117,10 +119,12 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.DoesNotThrow(() => ArrayStatistics.Minimum(data));
             Assert.DoesNotThrow(() => ArrayStatistics.Maximum(data));
             Assert.DoesNotThrow(() => ArrayStatistics.Mean(data));
+            Assert.DoesNotThrow(() => ArrayStatistics.Variance(data));
 
             Assert.DoesNotThrow(() => StreamingStatistics.Minimum(data));
             Assert.DoesNotThrow(() => StreamingStatistics.Maximum(data));
             Assert.DoesNotThrow(() => StreamingStatistics.Mean(data));
+            Assert.DoesNotThrow(() => StreamingStatistics.Variance(data));
         }
 
         [TestCase("lottery")]
@@ -231,8 +235,10 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             AssertHelpers.AlmostEqual(2d, Statistics.StandardDeviation(gaussian.Samples().Take(10000)), 2);
 
             AssertHelpers.AlmostEqual(1e+9, ArrayStatistics.Mean(gaussian.Samples().Take(10000).ToArray()), 11);
+            AssertHelpers.AlmostEqual(4d, ArrayStatistics.Variance(gaussian.Samples().Take(10000).ToArray()), 1);
 
             AssertHelpers.AlmostEqual(1e+9, StreamingStatistics.Mean(gaussian.Samples().Take(10000)), 11);
+            AssertHelpers.AlmostEqual(4d, StreamingStatistics.Variance(gaussian.Samples().Take(10000)), 1);
         }
 
         [Test]
@@ -278,6 +284,12 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.That(Statistics.Variance(new double[0]), Is.NaN);
             Assert.That(Statistics.Variance(new[] { 2d }), Is.NaN);
             Assert.That(Statistics.Variance(new[] { 2d, 3d }), Is.Not.NaN);
+            Assert.That(ArrayStatistics.Variance(new double[0]), Is.NaN);
+            Assert.That(ArrayStatistics.Variance(new[] { 2d }), Is.NaN);
+            Assert.That(ArrayStatistics.Variance(new[] { 2d, 3d }), Is.Not.NaN);
+            Assert.That(StreamingStatistics.Variance(new double[0]), Is.NaN);
+            Assert.That(StreamingStatistics.Variance(new[] { 2d }), Is.NaN);
+            Assert.That(StreamingStatistics.Variance(new[] { 2d, 3d }), Is.Not.NaN);
         }
 
         [Test]