From 539c1b4331b27f27e1ad09d8d4053073cf9cc2ff Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Thu, 23 Apr 2015 14:11:47 +0300 Subject: [PATCH 1/4] playing around with moving stats --- src/Numerics/Numerics.csproj | 1 + src/Numerics/Statistics/MovingStatistics.cs | 251 ++++++++++++++++++ .../StatisticsTests/MovingStatisticsTests.cs | 66 +++++ src/UnitTests/UnitTests.csproj | 1 + 4 files changed, 319 insertions(+) create mode 100644 src/Numerics/Statistics/MovingStatistics.cs create mode 100644 src/UnitTests/StatisticsTests/MovingStatisticsTests.cs diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b783234e..f286723b 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -211,6 +211,7 @@ + diff --git a/src/Numerics/Statistics/MovingStatistics.cs b/src/Numerics/Statistics/MovingStatistics.cs new file mode 100644 index 00000000..6e97b720 --- /dev/null +++ b/src/Numerics/Statistics/MovingStatistics.cs @@ -0,0 +1,251 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Statistics +{ + /// + /// Running statistics over a window of data, allows updating by adding values. + /// + public class MovingStatistics + { + readonly double[] _oldValues; + readonly int _windowSize; + int _lastIndex; + double _m1; + double _m2; + double _m3; + double _m4; + double _max = double.NegativeInfinity; + double _min = double.PositiveInfinity; + + public MovingStatistics(int windowSize) + { + if (windowSize < 1) + { + throw new ArgumentException(string.Format(Resources.ArgumentMustBePositive), "windowSize"); + } + _windowSize = windowSize; + _oldValues = new double[_windowSize]; + } + + public MovingStatistics(int windowSize, IEnumerable values) : this(windowSize) + { + PushRange(values); + } + + public int WindowSize + { + get { return _windowSize; } + } + + /// + /// Gets the total number of samples. + /// + public long Count { get; private set; } + + /// + /// Returns the minimum value in the sample data. + /// Returns NaN if data is empty or if any entry is NaN. + /// + public double Minimum + { + get { return Count > 0 ? _min : double.NaN; } + } + + /// + /// Returns the maximum value in the sample data. + /// Returns NaN if data is empty or if any entry is NaN. + /// + public double Maximum + { + get { return Count > 0 ? _max : double.NaN; } + } + + /// + /// Evaluates the sample mean, an estimate of the population mean. + /// Returns NaN if data is empty or if any entry is NaN. + /// + public double Mean + { + get { return Count > 0 ? _m1 : double.NaN; } + } + + /// + /// Estimates the unbiased population variance from the provided samples. + /// On a dataset of size N will use an N-1 normalizer (Bessel's correction). + /// Returns NaN if data has less than two entries or if any entry is NaN. + /// + public double Variance + { + get { return Count < 2 ? double.NaN : _m2/(Count - 1); } + } + + /// + /// Evaluates the variance from the provided full population. + /// On a dataset of size N will use an N normalizer and would thus be biased if applied to a subset. + /// Returns NaN if data is empty or if any entry is NaN. + /// + public double PopulationVariance + { + get { return Count < 2 ? double.NaN : _m2/Count; } + } + + /// + /// Estimates the unbiased population standard deviation from the provided samples. + /// On a dataset of size N will use an N-1 normalizer (Bessel's correction). + /// Returns NaN if data has less than two entries or if any entry is NaN. + /// + public double StandardDeviation + { + get { return Count < 2 ? double.NaN : Math.Sqrt(_m2/(Count - 1)); } + } + + /// + /// Evaluates the standard deviation from the provided full population. + /// On a dataset of size N will use an N normalizer and would thus be biased if applied to a subset. + /// Returns NaN if data is empty or if any entry is NaN. + /// + public double PopulationStandardDeviation + { + get { return Count < 2 ? double.NaN : Math.Sqrt(_m2/Count); } + } + +/* /// + /// Estimates the unbiased population skewness from the provided samples. + /// Uses a normalizer (Bessel's correction; type 2). + /// Returns NaN if data has less than three entries or if any entry is NaN. + /// + public double Skewness + { + get { return Count < 3 ? double.NaN : (Count*_m3*Math.Sqrt(_m2/(Count - 1))/(_m2*_m2*(Count - 2)))*(Count - 1); } + } + + /// + /// Evaluates the population skewness from the full population. + /// Does not use a normalizer and would thus be biased if applied to a subset (type 1). + /// Returns NaN if data has less than two entries or if any entry is NaN. + /// + public double PopulationSkewness + { + get { return Count < 2 ? double.NaN : Math.Sqrt(Count)*_m3/Math.Pow(_m2, 1.5); } + } + + /// + /// Estimates the unbiased population kurtosis from the provided samples. + /// Uses a normalizer (Bessel's correction; type 2). + /// Returns NaN if data has less than four entries or if any entry is NaN. + /// + public double Kurtosis + { + get { return Count < 4 ? double.NaN : ((double) Count*Count - 1)/((Count - 2)*(Count - 3))*(Count*_m4/(_m2*_m2) - 3 + 6.0/(Count + 1)); } + } + + /// + /// Evaluates the population kurtosis from the full population. + /// Does not use a normalizer and would thus be biased if applied to a subset (type 1). + /// Returns NaN if data has less than three entries or if any entry is NaN. + /// + public double PopulationKurtosis + { + get { return Count < 3 ? double.NaN : Count*_m4/(_m2*_m2) - 3.0; } + }*/ + + /// + /// Update the running statistics by adding another observed sample (in-place). + /// + public void Push(double value) + { + if (Count < _windowSize) + { + _oldValues[Count] = value; + Count++; + var d = value - _m1; + var s = d/Count; +// var s2 = s * s; + var t = d*s*(Count - 1); + + _m1 += s; +// _m4 += t * s2 * (Count * Count - 3 * Count + 3) + 6 * s2 * _m2 - 4 * s * _m3; +// _m3 += t * s * (Count - 2) - 3 * s * _m2; + _m2 += t; + + if (value < _min || double.IsNaN(value)) + { + _min = value; + } + + if (value > _max || double.IsNaN(value)) + { + _max = value; + } + } + else + { + var oldValue = _oldValues[_lastIndex]; + var d = value - oldValue; + var s = d/Count; +// var s2 = s * s; + var oldM1 = _m1; + _m1 += s; + + var x = (value - _m1 + oldValue - oldM1); + var t = d*x; + _m2 += t; + +// _m4 += t * s2 * (Count * Count - 3 * Count + 3) + 6 * s2 * _m2 - 4 * s * _m3; +// _m3 += t * (x /(Count-1) - 3 * s * _m2; + + _oldValues[_lastIndex] = value; + _lastIndex++; + if (_lastIndex == WindowSize) + { + _lastIndex = 0; + } + _max = value > _max || double.IsNaN(value) ? value : _oldValues.Maximum(); + _min = value < _min || double.IsNaN(value)? value : _oldValues.Minimum(); + } + } + + /// + /// Update the running statistics by adding a sequence of observed sample (in-place). + /// + public void PushRange(IEnumerable values) + { + foreach (var value in values) + { + Push(value); + } + } + } +} diff --git a/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs new file mode 100644 index 00000000..60d77f07 --- /dev/null +++ b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs @@ -0,0 +1,66 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; +using System.Collections.Generic; +using System.Diagnostics; +using System.Runtime.InteropServices; +using MathNet.Numerics.Distributions; +using MathNet.Numerics.Random; +using MathNet.Numerics.Statistics; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.StatisticsTests +{ +#if !PORTABLE + [TestFixture, Category("Statistics")] + public class MovingStatisticsTests + { + [Test] + public void QuickTest() + { + var data = new double[1000000]; + (new Normal(50, 10, new SystemRandomSource(0))).Samples(data); + var ms = new MovingStatistics(5, data); + ms.PushRange(new[] { 11.11, 22.22, 33.33, 44.44, 55.55 }); + Assert.AreEqual(5, ms.Count); + Assert.AreEqual(11.11, ms.Minimum); + Assert.AreEqual(55.55, ms.Maximum); + + Assert.AreEqual(33.33, ms.Mean, 1e-11); + Assert.AreEqual(308.58025, ms.Variance, 1e-10); + + //AssertHelpers.AlmostEqualRelative(stats0.Mean, ms.Mean, 14); + //AssertHelpers.AlmostEqualRelative(stats0.Variance, ms.Variance, 14); + //AssertHelpers.AlmostEqualRelative(stats0.StandardDeviation, ms.StandardDeviation, 14); + } + } +#endif +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 3b9505a2..713d8740 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -377,6 +377,7 @@ + From ac29a58b3d3ab182da2ff1613c17d8d61c7da71a Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 25 Apr 2015 19:29:59 +0200 Subject: [PATCH 2/4] Statistics: MovingStatistics should handle NaN such that it only affects its window --- src/Numerics/Statistics/MovingStatistics.cs | 72 +++++++++++++------ .../StatisticsTests/MovingStatisticsTests.cs | 66 +++++++++++++---- 2 files changed, 104 insertions(+), 34 deletions(-) diff --git a/src/Numerics/Statistics/MovingStatistics.cs b/src/Numerics/Statistics/MovingStatistics.cs index 6e97b720..71c0ee44 100644 --- a/src/Numerics/Statistics/MovingStatistics.cs +++ b/src/Numerics/Statistics/MovingStatistics.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2015 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -41,7 +41,12 @@ namespace MathNet.Numerics.Statistics { readonly double[] _oldValues; readonly int _windowSize; + + long _count; + long _totalCountOffset; int _lastIndex; + int _lastNaNTimeToLive; + double _m1; double _m2; double _m3; @@ -72,7 +77,10 @@ namespace MathNet.Numerics.Statistics /// /// Gets the total number of samples. /// - public long Count { get; private set; } + public long Count + { + get { return _totalCountOffset + _count; } + } /// /// Returns the minimum value in the sample data. @@ -80,7 +88,7 @@ namespace MathNet.Numerics.Statistics /// public double Minimum { - get { return Count > 0 ? _min : double.NaN; } + get { return _count > 0 && _lastNaNTimeToLive == 0 ? _min : double.NaN; } } /// @@ -89,7 +97,7 @@ namespace MathNet.Numerics.Statistics /// public double Maximum { - get { return Count > 0 ? _max : double.NaN; } + get { return _count > 0 && _lastNaNTimeToLive == 0 ? _max : double.NaN; } } /// @@ -98,7 +106,7 @@ namespace MathNet.Numerics.Statistics /// public double Mean { - get { return Count > 0 ? _m1 : double.NaN; } + get { return _count > 0 && _lastNaNTimeToLive == 0 ? _m1 : double.NaN; } } /// @@ -108,7 +116,7 @@ namespace MathNet.Numerics.Statistics /// public double Variance { - get { return Count < 2 ? double.NaN : _m2/(Count - 1); } + get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : _m2/(_count - 1); } } /// @@ -118,7 +126,7 @@ namespace MathNet.Numerics.Statistics /// public double PopulationVariance { - get { return Count < 2 ? double.NaN : _m2/Count; } + get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : _m2/_count; } } /// @@ -128,7 +136,7 @@ namespace MathNet.Numerics.Statistics /// public double StandardDeviation { - get { return Count < 2 ? double.NaN : Math.Sqrt(_m2/(Count - 1)); } + get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : Math.Sqrt(_m2/(_count - 1)); } } /// @@ -138,7 +146,7 @@ namespace MathNet.Numerics.Statistics /// public double PopulationStandardDeviation { - get { return Count < 2 ? double.NaN : Math.Sqrt(_m2/Count); } + get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : Math.Sqrt(_m2/_count); } } /* /// @@ -186,26 +194,46 @@ namespace MathNet.Numerics.Statistics /// public void Push(double value) { - if (Count < _windowSize) + if (double.IsNaN(value)) + { + _totalCountOffset += _count + 1; + _count = 0; + _lastNaNTimeToLive = _windowSize; + + _m1 = 0.0; + _m2 = 0.0; + _m3 = 0.0; + _m4 = 0.0; + _max = double.NegativeInfinity; + _min = double.PositiveInfinity; + return; + } + + if (_lastNaNTimeToLive > 0) + { + _lastNaNTimeToLive--; + } + + if (_count < _windowSize) { - _oldValues[Count] = value; - Count++; + _oldValues[_count] = value; + _count++; var d = value - _m1; - var s = d/Count; + var s = d/_count; // var s2 = s * s; - var t = d*s*(Count - 1); + var t = d*s*(_count - 1); _m1 += s; // _m4 += t * s2 * (Count * Count - 3 * Count + 3) + 6 * s2 * _m2 - 4 * s * _m3; // _m3 += t * s * (Count - 2) - 3 * s * _m2; _m2 += t; - if (value < _min || double.IsNaN(value)) + if (value < _min) { _min = value; } - if (value > _max || double.IsNaN(value)) + if (value > _max) { _max = value; } @@ -214,7 +242,7 @@ namespace MathNet.Numerics.Statistics { var oldValue = _oldValues[_lastIndex]; var d = value - oldValue; - var s = d/Count; + var s = d/_count; // var s2 = s * s; var oldM1 = _m1; _m1 += s; @@ -232,8 +260,8 @@ namespace MathNet.Numerics.Statistics { _lastIndex = 0; } - _max = value > _max || double.IsNaN(value) ? value : _oldValues.Maximum(); - _min = value < _min || double.IsNaN(value)? value : _oldValues.Minimum(); + _max = value > _max ? value : _oldValues.Maximum(); + _min = value < _min ? value : _oldValues.Minimum(); } } diff --git a/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs index 60d77f07..3753d706 100644 --- a/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs @@ -3,9 +3,9 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// +// // Copyright (c) 2009-2015 Math.NET -// +// // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -14,10 +14,10 @@ // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: -// +// // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. -// +// // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND @@ -28,10 +28,6 @@ // OTHER DEALINGS IN THE SOFTWARE. // -using System; -using System.Collections.Generic; -using System.Diagnostics; -using System.Runtime.InteropServices; using MathNet.Numerics.Distributions; using MathNet.Numerics.Random; using MathNet.Numerics.Statistics; @@ -39,17 +35,18 @@ using NUnit.Framework; namespace MathNet.Numerics.UnitTests.StatisticsTests { -#if !PORTABLE [TestFixture, Category("Statistics")] public class MovingStatisticsTests { [Test] - public void QuickTest() + public void StabilityTest() { var data = new double[1000000]; - (new Normal(50, 10, new SystemRandomSource(0))).Samples(data); + Normal.Samples(new SystemRandomSource(0), data, 50, 10); + var ms = new MovingStatistics(5, data); ms.PushRange(new[] { 11.11, 22.22, 33.33, 44.44, 55.55 }); + Assert.AreEqual(5, ms.Count); Assert.AreEqual(11.11, ms.Minimum); Assert.AreEqual(55.55, ms.Maximum); @@ -61,6 +58,51 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests //AssertHelpers.AlmostEqualRelative(stats0.Variance, ms.Variance, 14); //AssertHelpers.AlmostEqualRelative(stats0.StandardDeviation, ms.StandardDeviation, 14); } + + [Test] + public void NaNTest() + { + var ms = new MovingStatistics(3); + Assert.That(ms.Minimum, Is.NaN); + Assert.That(ms.Maximum, Is.NaN); + Assert.That(ms.Mean, Is.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(1.0); + Assert.That(ms.Minimum, Is.Not.NaN); + Assert.That(ms.Maximum, Is.Not.NaN); + Assert.That(ms.Mean, Is.Not.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.NaN); + Assert.That(ms.Maximum, Is.Not.NaN); + Assert.That(ms.Mean, Is.Not.NaN); + Assert.That(ms.StandardDeviation, Is.Not.NaN); + + ms.Push(double.NaN); + Assert.That(ms.Minimum, Is.NaN); + Assert.That(ms.Maximum, Is.NaN); + Assert.That(ms.Mean, Is.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(1.0); + Assert.That(ms.Minimum, Is.NaN); + Assert.That(ms.Maximum, Is.NaN); + Assert.That(ms.Mean, Is.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(2.0); + Assert.That(ms.Minimum, Is.NaN); + Assert.That(ms.Maximum, Is.NaN); + Assert.That(ms.Mean, Is.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(3.0); + Assert.That(ms.Minimum, Is.Not.NaN); + Assert.That(ms.Maximum, Is.Not.NaN); + Assert.That(ms.Mean, Is.Not.NaN); + Assert.That(ms.StandardDeviation, Is.Not.NaN); + } } -#endif } From 1baaee0fd50a3e071fbb35f2417aec959dd3049f Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sat, 25 Apr 2015 22:19:44 +0200 Subject: [PATCH 3/4] Statistics: MovingAverage --- src/Numerics/Statistics/Statistics.cs | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/Numerics/Statistics/Statistics.cs b/src/Numerics/Statistics/Statistics.cs index 3fe30641..9a92fafe 100644 --- a/src/Numerics/Statistics/Statistics.cs +++ b/src/Numerics/Statistics/Statistics.cs @@ -931,5 +931,21 @@ namespace MathNet.Numerics.Statistics { return StreamingStatistics.Entropy(data.Where(d => d.HasValue).Select(d => d.Value)); } + + /// + /// Evaluates the sample mean over a moving window, for each samples. + /// Returns NaN if no data is empty or if any entry is NaN. + /// + /// The sample stream to calculate the mean of. + /// The number of last samples to consider. + public static IEnumerable MovingAverage(this IEnumerable samples, int windowSize) + { + var movingStatistics = new MovingStatistics(windowSize); + return samples.Select(sample => + { + movingStatistics.Push(sample); + return movingStatistics.Mean; + }); + } } } From 644ea702de0f5e76d86ee11aa8ee0606ef1309e6 Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Sun, 3 May 2015 16:33:38 +0300 Subject: [PATCH 4/4] added support for infinity --- src/Numerics/Statistics/MovingStatistics.cs | 218 ++++++++++++------ .../StatisticsTests/MovingStatisticsTests.cs | 155 ++++++++++++- 2 files changed, 300 insertions(+), 73 deletions(-) diff --git a/src/Numerics/Statistics/MovingStatistics.cs b/src/Numerics/Statistics/MovingStatistics.cs index 71c0ee44..44cf1662 100644 --- a/src/Numerics/Statistics/MovingStatistics.cs +++ b/src/Numerics/Statistics/MovingStatistics.cs @@ -46,11 +46,11 @@ namespace MathNet.Numerics.Statistics long _totalCountOffset; int _lastIndex; int _lastNaNTimeToLive; + int _lastPosInfTimeToLive; + int _lastNegInfTimeToLive; double _m1; double _m2; - double _m3; - double _m4; double _max = double.NegativeInfinity; double _min = double.PositiveInfinity; @@ -64,7 +64,8 @@ namespace MathNet.Numerics.Statistics _oldValues = new double[_windowSize]; } - public MovingStatistics(int windowSize, IEnumerable values) : this(windowSize) + public MovingStatistics(int windowSize, IEnumerable values) + : this(windowSize) { PushRange(values); } @@ -88,7 +89,20 @@ namespace MathNet.Numerics.Statistics /// public double Minimum { - get { return _count > 0 && _lastNaNTimeToLive == 0 ? _min : double.NaN; } + get + { + if (_lastNaNTimeToLive > 0) + { + return double.NaN; + } + + if (_lastNegInfTimeToLive > 0) + { + return double.NegativeInfinity; + } + + return (_count > 0 || _lastPosInfTimeToLive > 0) ? _min : double.NaN; + } } /// @@ -97,7 +111,20 @@ namespace MathNet.Numerics.Statistics /// public double Maximum { - get { return _count > 0 && _lastNaNTimeToLive == 0 ? _max : double.NaN; } + get + { + if (_lastNaNTimeToLive > 0) + { + return double.NaN; + } + + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } + + return (_count > 0 || _lastNegInfTimeToLive > 0) ? _max : double.NaN; + } } /// @@ -106,7 +133,25 @@ namespace MathNet.Numerics.Statistics /// public double Mean { - get { return _count > 0 && _lastNaNTimeToLive == 0 ? _m1 : double.NaN; } + get + { + if (_lastNaNTimeToLive > 0 || (_lastPosInfTimeToLive > 0 && _lastNegInfTimeToLive > 0)) + { + return double.NaN; + } + + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } + + if (_lastNegInfTimeToLive > 0) + { + return double.NegativeInfinity; + } + + return _count == 0 ? double.NaN : _m1; + } } /// @@ -116,7 +161,20 @@ namespace MathNet.Numerics.Statistics /// public double Variance { - get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : _m2/(_count - 1); } + get + { + if (_lastNaNTimeToLive > 0 || _lastNegInfTimeToLive > 0 || (_lastPosInfTimeToLive > 0 && _lastNegInfTimeToLive > 0)) + { + return double.NaN; + } + + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } + + return _count < 2 ? double.NaN : _m2 / (_count - 1); + } } /// @@ -126,7 +184,20 @@ namespace MathNet.Numerics.Statistics /// public double PopulationVariance { - get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : _m2/_count; } + get + { + if (_lastNaNTimeToLive > 0 || _lastNegInfTimeToLive > 0 || (_lastPosInfTimeToLive > 0 && _lastNegInfTimeToLive > 0)) + { + return double.NaN; + } + + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } + + return _count < 2 ? double.NaN : _m2 / _count; + } } /// @@ -136,7 +207,20 @@ namespace MathNet.Numerics.Statistics /// public double StandardDeviation { - get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : Math.Sqrt(_m2/(_count - 1)); } + get + { + if (_lastNaNTimeToLive > 0 || _lastNegInfTimeToLive > 0 || (_lastPosInfTimeToLive > 0 && _lastNegInfTimeToLive > 0)) + { + return double.NaN; + } + + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } + + return _count < 2 ? double.NaN : Math.Sqrt(_m2 / (_count - 1)); + } } /// @@ -146,72 +230,48 @@ namespace MathNet.Numerics.Statistics /// public double PopulationStandardDeviation { - get { return _count < 2 || _lastNaNTimeToLive > 0 ? double.NaN : Math.Sqrt(_m2/_count); } - } - -/* /// - /// Estimates the unbiased population skewness from the provided samples. - /// Uses a normalizer (Bessel's correction; type 2). - /// Returns NaN if data has less than three entries or if any entry is NaN. - /// - public double Skewness - { - get { return Count < 3 ? double.NaN : (Count*_m3*Math.Sqrt(_m2/(Count - 1))/(_m2*_m2*(Count - 2)))*(Count - 1); } - } + get + { + if (_lastNaNTimeToLive > 0 || _lastNegInfTimeToLive > 0 || (_lastPosInfTimeToLive > 0 && _lastNegInfTimeToLive > 0)) + { + return double.NaN; + } - /// - /// Evaluates the population skewness from the full population. - /// Does not use a normalizer and would thus be biased if applied to a subset (type 1). - /// Returns NaN if data has less than two entries or if any entry is NaN. - /// - public double PopulationSkewness - { - get { return Count < 2 ? double.NaN : Math.Sqrt(Count)*_m3/Math.Pow(_m2, 1.5); } - } + if (_lastPosInfTimeToLive > 0) + { + return double.PositiveInfinity; + } - /// - /// Estimates the unbiased population kurtosis from the provided samples. - /// Uses a normalizer (Bessel's correction; type 2). - /// Returns NaN if data has less than four entries or if any entry is NaN. - /// - public double Kurtosis - { - get { return Count < 4 ? double.NaN : ((double) Count*Count - 1)/((Count - 2)*(Count - 3))*(Count*_m4/(_m2*_m2) - 3 + 6.0/(Count + 1)); } + return _count < 2 ? double.NaN : Math.Sqrt(_m2 / _count); + } } - /// - /// Evaluates the population kurtosis from the full population. - /// Does not use a normalizer and would thus be biased if applied to a subset (type 1). - /// Returns NaN if data has less than three entries or if any entry is NaN. - /// - public double PopulationKurtosis - { - get { return Count < 3 ? double.NaN : Count*_m4/(_m2*_m2) - 3.0; } - }*/ - /// /// Update the running statistics by adding another observed sample (in-place). /// public void Push(double value) { + DecrementTimeToLive(); + if (double.IsNaN(value)) { - _totalCountOffset += _count + 1; - _count = 0; _lastNaNTimeToLive = _windowSize; + Reset(double.PositiveInfinity, double.NegativeInfinity); + return; + } - _m1 = 0.0; - _m2 = 0.0; - _m3 = 0.0; - _m4 = 0.0; - _max = double.NegativeInfinity; - _min = double.PositiveInfinity; + if (double.IsPositiveInfinity(value)) + { + _lastPosInfTimeToLive = _windowSize; + Reset(_min, double.NegativeInfinity); return; } - if (_lastNaNTimeToLive > 0) + if (double.IsNegativeInfinity(value)) { - _lastNaNTimeToLive--; + _lastNegInfTimeToLive = _windowSize; + Reset(double.PositiveInfinity, _max); + return; } if (_count < _windowSize) @@ -219,13 +279,10 @@ namespace MathNet.Numerics.Statistics _oldValues[_count] = value; _count++; var d = value - _m1; - var s = d/_count; -// var s2 = s * s; - var t = d*s*(_count - 1); + var s = d / _count; + var t = d * s * (_count - 1); _m1 += s; -// _m4 += t * s2 * (Count * Count - 3 * Count + 3) + 6 * s2 * _m2 - 4 * s * _m3; -// _m3 += t * s * (Count - 2) - 3 * s * _m2; _m2 += t; if (value < _min) @@ -242,18 +299,14 @@ namespace MathNet.Numerics.Statistics { var oldValue = _oldValues[_lastIndex]; var d = value - oldValue; - var s = d/_count; -// var s2 = s * s; + var s = d / _count; var oldM1 = _m1; _m1 += s; var x = (value - _m1 + oldValue - oldM1); - var t = d*x; + var t = d * x; _m2 += t; -// _m4 += t * s2 * (Count * Count - 3 * Count + 3) + 6 * s2 * _m2 - 4 * s * _m3; -// _m3 += t * (x /(Count-1) - 3 * s * _m2; - _oldValues[_lastIndex] = value; _lastIndex++; if (_lastIndex == WindowSize) @@ -275,5 +328,32 @@ namespace MathNet.Numerics.Statistics Push(value); } } + + private void DecrementTimeToLive() + { + if (_lastNaNTimeToLive > 0) + { + _lastNaNTimeToLive--; + } + + if (_lastPosInfTimeToLive > 0) + { + _lastPosInfTimeToLive--; + } + + if (_lastNegInfTimeToLive > 0) + { + _lastNegInfTimeToLive--; + } + } + + private void Reset(double min, double max) + { + _totalCountOffset += _count + 1; + _count = 0; + _m1 = 0; + _max = max; + _min = min; + } } } diff --git a/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs index 3753d706..867ab4a6 100644 --- a/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs +++ b/src/UnitTests/StatisticsTests/MovingStatisticsTests.cs @@ -28,6 +28,7 @@ // OTHER DEALINGS IN THE SOFTWARE. // +using System; using MathNet.Numerics.Distributions; using MathNet.Numerics.Random; using MathNet.Numerics.Statistics; @@ -38,6 +39,153 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests [TestFixture, Category("Statistics")] public class MovingStatisticsTests { + [Test] + public void PositiveInfinityTest() + { + var ms = new MovingStatistics(3); + + ms.Push(1.0); + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.EqualTo(double.PositiveInfinity)); + + ms.Push(double.PositiveInfinity); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(1.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(double.PositiveInfinity); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(3.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(4.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.EqualTo(double.PositiveInfinity)); + } + + [Test] + public void NegativeInfinityTest() + { + var ms = new MovingStatistics(3); + + ms.Push(1.0); + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.EqualTo(double.NegativeInfinity)); + + ms.Push(double.NegativeInfinity); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(1.0); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(double.NegativeInfinity); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(2.0); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(3.0); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(4.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.NaN); + } + + [Test] + public void MixedInfinityTest() + { + var ms = new MovingStatistics(3); + + ms.Push(1.0); + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.EqualTo(double.PositiveInfinity)); + + ms.Push(double.NegativeInfinity); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(1.0); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(double.PositiveInfinity); + Assert.That(ms.Minimum, Is.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.NaN); + Assert.That(ms.StandardDeviation, Is.NaN); + + ms.Push(2.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(3.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.EqualTo(double.PositiveInfinity)); + + ms.Push(4.0); + Assert.That(ms.Minimum, Is.Not.EqualTo(double.NegativeInfinity)); + Assert.That(ms.Maximum, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.Mean, Is.Not.EqualTo(double.PositiveInfinity)); + Assert.That(ms.StandardDeviation, Is.Not.EqualTo(double.PositiveInfinity)); + } + [Test] public void StabilityTest() { @@ -53,10 +201,9 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests Assert.AreEqual(33.33, ms.Mean, 1e-11); Assert.AreEqual(308.58025, ms.Variance, 1e-10); - - //AssertHelpers.AlmostEqualRelative(stats0.Mean, ms.Mean, 14); - //AssertHelpers.AlmostEqualRelative(stats0.Variance, ms.Variance, 14); - //AssertHelpers.AlmostEqualRelative(stats0.StandardDeviation, ms.StandardDeviation, 14); + Assert.AreEqual(17.5664524022354, ms.StandardDeviation, 1e-11); + Assert.AreEqual(246.8642, ms.PopulationVariance, 1e-10); + Assert.AreEqual(15.7119126779651, ms.PopulationStandardDeviation, 1e-10); } [Test]