forked from tsai/mathnet-numerics
Christoph Ruegg
13 years ago
5 changed files with 451 additions and 1 deletions
@ -0,0 +1,240 @@ |
|||
// <copyright file="DescriptiveStatistics.cs" company="Math.NET">
|
|||
// 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-2014 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.
|
|||
// </copyright>
|
|||
//
|
|||
// Adapted from the old DescriptiveStatistics and inspired in design
|
|||
// among others by http://www.johndcook.com/skewness_kurtosis.html
|
|||
|
|||
using System; |
|||
using System.Collections.Generic; |
|||
|
|||
namespace MathNet.Numerics.Statistics |
|||
{ |
|||
/// <summary>
|
|||
/// Running statistics, allows updating by adding values,
|
|||
/// or combining
|
|||
/// </summary>
|
|||
public class RunningStatistics |
|||
{ |
|||
long _n; |
|||
double _m1; |
|||
double _m2; |
|||
double _m3; |
|||
double _m4; |
|||
double _min = Double.PositiveInfinity; |
|||
double _max = Double.NegativeInfinity; |
|||
|
|||
public RunningStatistics() |
|||
{ |
|||
} |
|||
|
|||
public RunningStatistics(IEnumerable<double> values) |
|||
{ |
|||
PushRange(values); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets the total number of samples.
|
|||
/// </summary>
|
|||
public long Count |
|||
{ |
|||
get { return _n; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the minimum value in the sample data.
|
|||
/// Returns NaN if data is empty or if any entry is NaN.
|
|||
/// </summary>
|
|||
public double Minimum |
|||
{ |
|||
get { return _n > 0 ? _min : double.NaN; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Returns the maximum value in the sample data.
|
|||
/// Returns NaN if data is empty or if any entry is NaN.
|
|||
/// </summary>
|
|||
public double Maximum |
|||
{ |
|||
get { return _n > 0 ? _max : double.NaN; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Evaluates the sample mean, an estimate of the population mean.
|
|||
/// Returns NaN if data is empty or if any entry is NaN.
|
|||
/// </summary>
|
|||
public double Mean |
|||
{ |
|||
get { return _n > 0 ? _m1 : double.NaN; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double Variance |
|||
{ |
|||
get { return _n < 2 ? double.NaN : _m2/(_n - 1); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double PopulationVariance |
|||
{ |
|||
get { return _n < 2 ? double.NaN : _m2/_n; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double StandardDeviation |
|||
{ |
|||
get { return _n < 2 ? double.NaN : Math.Sqrt(_m2/(_n - 1)); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double PopulationStandardDeviation |
|||
{ |
|||
get { return _n < 2 ? double.NaN : Math.Sqrt(_m2/_n); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double Skewness |
|||
{ |
|||
get { return _n < 3 ? double.NaN : (_n*_m3*Math.Sqrt(_m2/(_n - 1))/(_m2*_m2*(_n - 2)))*(_n - 1); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double PopulationSkewness |
|||
{ |
|||
get { return _n < 2 ? double.NaN : _m3*Math.Sqrt(_n*(_n - 1))*Math.Sqrt(_m2/(_n - 1))/(_m2*_m2); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double Kurtosis |
|||
{ |
|||
get { return _n < 4 ? double.NaN : ((double)_n*_n - 1)/((_n - 2)*(_n - 3))*(_n*_m4/(_m2*_m2) - 3 + 6.0/(_n + 1)); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// 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.
|
|||
/// </summary>
|
|||
public double PopulationKurtosis |
|||
{ |
|||
get { return _n < 3 ? double.NaN : (_m4*_n - 3*_m2*_m2)/(_m2*_m2); } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Update the running statistics by adding another observed sample (in-place).
|
|||
/// </summary>
|
|||
public void Push(double value) |
|||
{ |
|||
_n++; |
|||
double d = value - _m1; |
|||
double s = d/_n; |
|||
double s2 = s*s; |
|||
double t = d*s*(_n - 1); |
|||
|
|||
_m1 += s; |
|||
_m4 += t*s2*(_n*_n - 3*_n + 3) + 6*s2*_m2 - 4*s*_m3; |
|||
_m3 += t*s*(_n - 2) - 3*s*_m2; |
|||
_m2 += t; |
|||
|
|||
if (_min > value) |
|||
{ |
|||
_min = value; |
|||
} |
|||
if (_max < value) |
|||
{ |
|||
_max = value; |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Update the running statistics by adding a sequence of observed sample (in-place).
|
|||
/// </summary>
|
|||
public void PushRange(IEnumerable<double> values) |
|||
{ |
|||
foreach (double value in values) |
|||
{ |
|||
Push(value); |
|||
} |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Create a new running statistics over the combined samples of two existing running statistics.
|
|||
/// </summary>
|
|||
public static RunningStatistics Combine(RunningStatistics a, RunningStatistics b) |
|||
{ |
|||
long n = a._n + b._n; |
|||
double d = b._m1 - a._m1; |
|||
double d2 = d*d; |
|||
double d3 = d2*d; |
|||
double d4 = d2*d2; |
|||
|
|||
double m1 = (a._n*a._m1 + b._n*b._m1)/n; |
|||
double m2 = a._m2 + b._m2 + d2*a._n*b._n/n; |
|||
double m3 = a._m3 + b._m3 + d3*a._n*b._n*(a._n - b._n)/(n*n) |
|||
+ 3*d3*(a._n*b._m2 - b._n*a._m2)/n; |
|||
double m4 = a._m4 + b._m4 + d4*a._n*b._n*(a._n*a._n - a._n*b._n + b._n*b._n)/(n*n*n) |
|||
+ 6*d2*(a._n*a._n*b._m2 + b._n*b._n*a._m2)/(n*n) + 4*d*(a._n*b._m3 - b._n*a._m3)/n; |
|||
|
|||
return new RunningStatistics { _n = n, _m1 = m1, _m2 = m2, _m3 = m3, _m4 = m4 }; |
|||
} |
|||
|
|||
public static RunningStatistics operator +(RunningStatistics a, RunningStatistics b) |
|||
{ |
|||
return Combine(a, b); |
|||
} |
|||
} |
|||
} |
|||
@ -0,0 +1,170 @@ |
|||
// <copyright file="DescriptiveStatisticsTests.cs" company="Math.NET">
|
|||
// 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-2014 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.
|
|||
// </copyright>
|
|||
|
|||
namespace MathNet.Numerics.UnitTests.StatisticsTests |
|||
{ |
|||
#if !PORTABLE
|
|||
using System.Collections.Generic; |
|||
using NUnit.Framework; |
|||
using Statistics; |
|||
|
|||
/// <summary>
|
|||
/// Running statistics tests.
|
|||
/// </summary>
|
|||
/// <remarks>NOTE: this class is not included into Silverlight version, because it uses data from local files.
|
|||
/// In Silverlight access to local files is forbidden, except several cases.</remarks>
|
|||
[TestFixture, Category("Statistics")] |
|||
public class RunningStatisticsTests |
|||
{ |
|||
/// <summary>
|
|||
/// Statistics data.
|
|||
/// </summary>
|
|||
readonly IDictionary<string, StatTestData> _data = new Dictionary<string, StatTestData>(); |
|||
|
|||
/// <summary>
|
|||
/// Initializes a new instance of the DescriptiveStatisticsTests class.
|
|||
/// </summary>
|
|||
public RunningStatisticsTests() |
|||
{ |
|||
_data.Add("lottery", new StatTestData("./data/NIST/Lottery.dat")); |
|||
_data.Add("lew", new StatTestData("./data/NIST/Lew.dat")); |
|||
_data.Add("mavro", new StatTestData("./data/NIST/Mavro.dat")); |
|||
_data.Add("michelso", new StatTestData("./data/NIST/Michelso.dat")); |
|||
_data.Add("numacc1", new StatTestData("./data/NIST/NumAcc1.dat")); |
|||
_data.Add("numacc2", new StatTestData("./data/NIST/NumAcc2.dat")); |
|||
_data.Add("numacc3", new StatTestData("./data/NIST/NumAcc3.dat")); |
|||
_data.Add("numacc4", new StatTestData("./data/NIST/NumAcc4.dat")); |
|||
_data.Add("meixner", new StatTestData("./data/NIST/Meixner.dat")); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// <c>IEnumerable</c> Double.
|
|||
/// </summary>
|
|||
/// <param name="dataSet">Dataset name.</param>
|
|||
/// <param name="digits">Digits count.</param>
|
|||
/// <param name="skewness">Skewness value.</param>
|
|||
/// <param name="kurtosis">Kurtosis value.</param>
|
|||
/// <param name="median">Median value.</param>
|
|||
/// <param name="min">Min value.</param>
|
|||
/// <param name="max">Max value.</param>
|
|||
/// <param name="count">Count value.</param>
|
|||
[TestCase("lottery", 14, -0.09333165310779, -1.19256091074856, 522.5, 4, 999, 218)] |
|||
[TestCase("lew", 14, -0.050606638756334, -1.49604979214447, -162, -579, 300, 200)] |
|||
[TestCase("mavro", 11, 0.64492948110824, -0.82052379677456, 2.0018, 2.0013, 2.0027, 50)] |
|||
[TestCase("michelso", 11, -0.0185388637725746, 0.33968459842539, 299.85, 299.62, 300.07, 100)] |
|||
[TestCase("numacc1", 15, 0, double.NaN, 10000002, 10000001, 10000003, 3)] |
|||
[TestCase("numacc2", 13, 0, -2.003003003003, 1.2, 1.1, 1.3, 1001)] |
|||
[TestCase("numacc3", 9, 0, -2.003003003003, 1000000.2, 1000000.1, 1000000.3, 1001)] |
|||
[TestCase("numacc4", 7, 0, -2.00300300299913, 10000000.2, 10000000.1, 10000000.3, 1001)] |
|||
[TestCase("meixner", 8, -0.016649617280859657, 0.8171318629552635, -0.002042931016531602, -4.825626912281697, 5.3018298664184913, 10000)] |
|||
public void ConsistentWithNist(string dataSet, int digits, double skewness, double kurtosis, double median, double min, double max, int count) |
|||
{ |
|||
var data = _data[dataSet]; |
|||
var stats = new RunningStatistics(data.Data); |
|||
|
|||
AssertHelpers.AlmostEqualRelative(data.Mean, stats.Mean, 10); |
|||
AssertHelpers.AlmostEqualRelative(data.StandardDeviation, stats.StandardDeviation, digits); |
|||
AssertHelpers.AlmostEqualRelative(skewness, stats.Skewness, 8); |
|||
AssertHelpers.AlmostEqualRelative(kurtosis, stats.Kurtosis, 8); |
|||
Assert.AreEqual(stats.Minimum, min); |
|||
Assert.AreEqual(stats.Maximum, max); |
|||
Assert.AreEqual(stats.Count, count); |
|||
} |
|||
|
|||
[TestCase("lottery", 1e-8, -0.09268823, -0.09333165)] |
|||
[TestCase("lew", 1e-8, -0.0502263, -0.05060664)] |
|||
[TestCase("mavro", 1e-6, 0.6254181, 0.6449295)] |
|||
[TestCase("michelso", 1e-8, -0.01825961, -0.01853886)] |
|||
[TestCase("numacc1", 1e-8, 0, 0)] |
|||
//[TestCase("numacc2", 1e-20, 3.254232e-15, 3.259118e-15)] TODO: accuracy
|
|||
//[TestCase("numacc3", 1e-14, 1.747103e-09, 1.749726e-09)] TODO: accuracy
|
|||
//[TestCase("numacc4", 1e-13, 2.795364e-08, 2.799561e-08)] TODO: accuracy
|
|||
[TestCase("meixner", 1e-8, -0.01664712, -0.01664962)] |
|||
public void SkewnessConsistentWithR_e1071(string dataSet, double delta, double skewnessType1, double skewnessType2) |
|||
{ |
|||
var data = _data[dataSet]; |
|||
var stats = new RunningStatistics(data.Data); |
|||
|
|||
Assert.That(stats.Skewness, Is.EqualTo(skewnessType2).Within(delta), "Skewness"); |
|||
Assert.That(stats.PopulationSkewness, Is.EqualTo(skewnessType1).Within(delta), "PopulationSkewness"); |
|||
} |
|||
|
|||
[TestCase("lottery", -1.192781, -1.192561)] |
|||
[TestCase("lew", -1.48876, -1.49605)] |
|||
[TestCase("mavro", -0.858384, -0.8205238)] |
|||
[TestCase("michelso", 0.2635305, 0.3396846)] |
|||
[TestCase("numacc1", -1.5, double.NaN)] |
|||
[TestCase("numacc2", -1.999, -2.003003)] |
|||
[TestCase("numacc3", -1.999, -2.003003)] |
|||
[TestCase("numacc4", -1.999, -2.003003)] |
|||
[TestCase("meixner", 0.8161234, 0.8171319)] |
|||
public void KurtosisConsistentWithR_e1071(string dataSet, double kurtosisType1, double kurtosisType2) |
|||
{ |
|||
var data = _data[dataSet]; |
|||
var stats = new RunningStatistics(data.Data); |
|||
|
|||
Assert.That(stats.Kurtosis, Is.EqualTo(kurtosisType2).Within(1e-6), "Kurtosis"); |
|||
Assert.That(stats.PopulationKurtosis, Is.EqualTo(kurtosisType1).Within(1e-6), "PopulationKurtosis"); |
|||
} |
|||
|
|||
[Test] |
|||
public void ShortSequences() |
|||
{ |
|||
var stats0 = new RunningStatistics(new double[0]); |
|||
Assert.That(stats0.Skewness, Is.NaN); |
|||
Assert.That(stats0.Kurtosis, Is.NaN); |
|||
|
|||
var stats1 = new RunningStatistics(new[] { 1.0 }); |
|||
Assert.That(stats1.Skewness, Is.NaN); |
|||
Assert.That(stats1.Kurtosis, Is.NaN); |
|||
|
|||
var stats2 = new RunningStatistics(new[] { 1.0, 2.0 }); |
|||
Assert.That(stats2.Skewness, Is.NaN); |
|||
Assert.That(stats2.Kurtosis, Is.NaN); |
|||
|
|||
var stats3 = new RunningStatistics(new[] { 1.0, 2.0, -3.0 }); |
|||
Assert.That(stats3.Skewness, Is.Not.NaN); |
|||
Assert.That(stats3.Kurtosis, Is.NaN); |
|||
|
|||
var stats4 = new RunningStatistics(new[] { 1.0, 2.0, -3.0, -4.0 }); |
|||
Assert.That(stats4.Skewness, Is.Not.NaN); |
|||
Assert.That(stats4.Kurtosis, Is.Not.NaN); |
|||
} |
|||
|
|||
[Test] |
|||
public void ZeroVarianceSequence() |
|||
{ |
|||
var stats = new RunningStatistics(new[] { 2.0, 2.0, 2.0, 2.0 }); |
|||
Assert.That(stats.Skewness, Is.NaN); |
|||
Assert.That(stats.Kurtosis, Is.NaN); |
|||
} |
|||
} |
|||
#endif
|
|||
} |
|||
Loading…
Reference in new issue