From 15a56daaf2390362d58e4e949bbb95354f671828 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 21 Mar 2013 15:04:33 +0100 Subject: [PATCH] Statistics: Array and Streaming Variance --- src/Numerics/Statistics/ArrayStatistics.cs | 24 +++++- src/Numerics/Statistics/Statistics.cs | 75 ++----------------- .../Statistics/StreamingStatistics.cs | 35 ++++++++- .../StatisticsTests/StatisticsTests.cs | 12 +++ 4 files changed, 75 insertions(+), 71 deletions(-) diff --git a/src/Numerics/Statistics/ArrayStatistics.cs b/src/Numerics/Statistics/ArrayStatistics.cs index b304a5c2..b37818cc 100644 --- a/src/Numerics/Statistics/ArrayStatistics.cs +++ b/src/Numerics/Statistics/ArrayStatistics.cs @@ -80,7 +80,7 @@ namespace MathNet.Numerics.Statistics } /// - /// 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. /// /// Sample array, no sorting is assumed. @@ -97,5 +97,27 @@ namespace MathNet.Numerics.Statistics } return mean; } + + /// + /// 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. + /// + /// Sample array, no sorting is assumed. + 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); + } } } \ No newline at end of file diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index f003987d..e673b294 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -118,34 +118,10 @@ namespace MathNet.Numerics.Statistics /// The unbiased population variance of the sample. public static double Variance(this IEnumerable data) { - if (data == null) - { - throw new ArgumentNullException("data"); - } - - double variance = 0; - double t = 0; - ulong j = 0; - - using (IEnumerator 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); } /// @@ -155,47 +131,8 @@ namespace MathNet.Numerics.Statistics /// The population variance of the sample. public static double Variance(this IEnumerable data) { - if (data == null) - { - throw new ArgumentNullException("data"); - } - - double variance = 0; - double t = 0; - ulong j = 0; - - using (IEnumerator 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)); } /// diff --git a/src/Numerics/Statistics/StreamingStatistics.cs b/src/Numerics/Statistics/StreamingStatistics.cs index 52c676eb..5f717a2b 100644 --- a/src/Numerics/Statistics/StreamingStatistics.cs +++ b/src/Numerics/Statistics/StreamingStatistics.cs @@ -80,7 +80,7 @@ namespace MathNet.Numerics.Statistics } /// - /// 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. /// /// Sample stream, no sorting is assumed. @@ -98,5 +98,38 @@ namespace MathNet.Numerics.Statistics } return any ? mean : double.NaN; } + + /// + /// 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. + /// + /// Sample stream, no sorting is assumed. + public static double Variance(IEnumerable 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; + } } } diff --git a/src/UnitTests/StatisticsTests/StatisticsTests.cs b/src/UnitTests/StatisticsTests/StatisticsTests.cs index e462aa8d..f3a0c192 100644 --- a/src/UnitTests/StatisticsTests/StatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/StatisticsTests.cs @@ -83,10 +83,12 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.Throws(() => ArrayStatistics.Minimum(data)); Assert.Throws(() => ArrayStatistics.Maximum(data)); Assert.Throws(() => ArrayStatistics.Mean(data)); + Assert.Throws(() => ArrayStatistics.Variance(data)); Assert.Throws(() => StreamingStatistics.Minimum(data)); Assert.Throws(() => StreamingStatistics.Maximum(data)); Assert.Throws(() => StreamingStatistics.Mean(data)); + Assert.Throws(() => 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]