Skip to content

Commit 38bbbe8

Browse files
cdrnetcdrnet
authored
Merge pull request #865 from rca22/rca22_weightedStats
Add classes generalizing running and descriptive statistics to (reliability) weighted samples
2 parents 97a905c + 3f74901 commit 38bbbe8

6 files changed

+1562
-34
lines changed

‎src/Numerics.Tests/StatisticsTests/RunningWeightedStatisticsTests.cs

+324
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,324 @@
1+
// <copyright file="DescriptiveStatisticsTests.cs" company="Math.NET">
2+
// Math.NET Numerics, part of the Math.NET Project
3+
// http://numerics.mathdotnet.com
4+
// http://github.com/mathnet/mathnet-numerics
5+
//
6+
// Copyright (c) 2009-2016 Math.NET
7+
//
8+
// Permission is hereby granted, free of charge, to any person
9+
// obtaining a copy of this software and associated documentation
10+
// files (the "Software"), to deal in the Software without
11+
// restriction, including without limitation the rights to use,
12+
// copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
// copies of the Software, and to permit persons to whom the
14+
// Software is furnished to do so, subject to the following
15+
// conditions:
16+
//
17+
// The above copyright notice and this permission notice shall be
18+
// included in all copies or substantial portions of the Software.
19+
//
20+
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21+
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
22+
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
23+
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
24+
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
25+
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26+
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
27+
// OTHER DEALINGS IN THE SOFTWARE.
28+
// </copyright>
29+
30+
using System.Linq;
31+
using System.Collections.Generic;
32+
using MathNet.Numerics.Distributions;
33+
using MathNet.Numerics.Random;
34+
using MathNet.Numerics.Statistics;
35+
using NUnit.Framework;
36+
37+
namespace MathNet.Numerics.UnitTests.StatisticsTests
38+
{
39+
/// <summary>
40+
/// Running statistics tests.
41+
/// </summary>
42+
/// <remarks>NOTE: this class is not included into Silverlight version, because it uses data from local files.
43+
/// In Silverlight access to local files is forbidden, except several cases.</remarks>
44+
[TestFixture, Category("Statistics")]
45+
public class RunningWeightedStatisticsTests
46+
{
47+
/// <summary>
48+
/// Statistics data.
49+
/// </summary>
50+
readonly IDictionary<string, StatTestData> _data = new Dictionary<string, StatTestData>();
51+
52+
/// <summary>
53+
/// Initializes a new instance of the DescriptiveStatisticsTests class.
54+
/// </summary>
55+
public RunningWeightedStatisticsTests()
56+
{
57+
_data.Add("lottery", new StatTestData("NIST.Lottery.dat"));
58+
_data.Add("lew", new StatTestData("NIST.Lew.dat"));
59+
_data.Add("mavro", new StatTestData("NIST.Mavro.dat"));
60+
_data.Add("michelso", new StatTestData("NIST.Michelso.dat"));
61+
_data.Add("numacc1", new StatTestData("NIST.NumAcc1.dat"));
62+
_data.Add("numacc2", new StatTestData("NIST.NumAcc2.dat"));
63+
_data.Add("numacc3", new StatTestData("NIST.NumAcc3.dat"));
64+
_data.Add("numacc4", new StatTestData("NIST.NumAcc4.dat"));
65+
_data.Add("meixner", new StatTestData("NIST.Meixner.dat"));
66+
}
67+
68+
/// <summary>
69+
/// <c>IEnumerable</c> Double.
70+
/// </summary>
71+
/// <param name="dataSet">Dataset name.</param>
72+
/// <param name="digits">Digits count.</param>
73+
/// <param name="skewness">Skewness value.</param>
74+
/// <param name="kurtosis">Kurtosis value.</param>
75+
/// <param name="median">Median value.</param>
76+
/// <param name="min">Min value.</param>
77+
/// <param name="max">Max value.</param>
78+
/// <param name="count">Count value.</param>
79+
[TestCase("lottery", 14, -0.09333165310779, -1.19256091074856, 522.5, 4, 999, 218)]
80+
[TestCase("lew", 14, -0.050606638756334, -1.49604979214447, -162, -579, 300, 200)]
81+
[TestCase("mavro", 11, 0.64492948110824, -0.82052379677456, 2.0018, 2.0013, 2.0027, 50)]
82+
[TestCase("michelso", 11, -0.0185388637725746, 0.33968459842539, 299.85, 299.62, 300.07, 100)]
83+
[TestCase("numacc1", 15, 0, double.NaN, 10000002, 10000001, 10000003, 3)]
84+
[TestCase("numacc2", 13, 0, -2.003003003003, 1.2, 1.1, 1.3, 1001)]
85+
[TestCase("numacc3", 9, 0, -2.003003003003, 1000000.2, 1000000.1, 1000000.3, 1001)]
86+
[TestCase("numacc4", 7, 0, -2.00300300299913, 10000000.2, 10000000.1, 10000000.3, 1001)]
87+
[TestCase("meixner", 8, -0.016649617280859657, 0.8171318629552635, -0.002042931016531602, -4.825626912281697, 5.3018298664184913, 10000)]
88+
public void ConsistentWithNist(string dataSet, int digits, double skewness, double kurtosis, double median, double min, double max, int count)
89+
{
90+
var data = _data[dataSet];
91+
var stats = new RunningWeightedStatistics(data.Data.Select(x => System.Tuple.Create(1.0, x)));
92+
93+
AssertHelpers.AlmostEqualRelative(data.Mean, stats.Mean, 10);
94+
AssertHelpers.AlmostEqualRelative(data.StandardDeviation, stats.StandardDeviation, digits);
95+
AssertHelpers.AlmostEqualRelative(skewness, stats.Skewness, 8);
96+
AssertHelpers.AlmostEqualRelative(kurtosis, stats.Kurtosis, 8);
97+
Assert.AreEqual(stats.Minimum, min);
98+
Assert.AreEqual(stats.Maximum, max);
99+
Assert.AreEqual(stats.Count, count);
100+
}
101+
102+
[TestCase("lottery", 1e-8, -0.09268823, -0.09333165)]
103+
[TestCase("lew", 1e-8, -0.0502263, -0.05060664)]
104+
[TestCase("mavro", 1e-6, 0.6254181, 0.6449295)]
105+
[TestCase("michelso", 1e-8, -0.01825961, -0.01853886)]
106+
[TestCase("numacc1", 1e-8, 0, 0)]
107+
//[TestCase("numacc2", 1e-20, 3.254232e-15, 3.259118e-15)] TODO: accuracy
108+
//[TestCase("numacc3", 1e-14, 1.747103e-09, 1.749726e-09)] TODO: accuracy
109+
//[TestCase("numacc4", 1e-13, 2.795364e-08, 2.799561e-08)] TODO: accuracy
110+
[TestCase("meixner", 1e-8, -0.01664712, -0.01664962)]
111+
public void SkewnessConsistentWithR_e1071(string dataSet, double delta, double skewnessType1, double skewnessType2)
112+
{
113+
var data = _data[dataSet];
114+
var stats = new RunningWeightedStatistics(data.Data.Select(x => System.Tuple.Create(1.0, x)));
115+
116+
Assert.That(stats.Skewness, Is.EqualTo(skewnessType2).Within(delta), "Skewness");
117+
Assert.That(stats.PopulationSkewness, Is.EqualTo(skewnessType1).Within(delta), "PopulationSkewness");
118+
}
119+
120+
[TestCase("lottery", -1.192781, -1.192561)]
121+
[TestCase("lew", -1.48876, -1.49605)]
122+
[TestCase("mavro", -0.858384, -0.8205238)]
123+
[TestCase("michelso", 0.2635305, 0.3396846)]
124+
[TestCase("numacc1", -1.5, double.NaN)]
125+
[TestCase("numacc2", -1.999, -2.003003)]
126+
[TestCase("numacc3", -1.999, -2.003003)]
127+
[TestCase("numacc4", -1.999, -2.003003)]
128+
[TestCase("meixner", 0.8161234, 0.8171319)]
129+
public void KurtosisConsistentWithR_e1071(string dataSet, double kurtosisType1, double kurtosisType2)
130+
{
131+
var data = _data[dataSet];
132+
var stats = new RunningWeightedStatistics(data.Data.Select(x => System.Tuple.Create(1.0, x)));
133+
134+
Assert.That(stats.Kurtosis, Is.EqualTo(kurtosisType2).Within(1e-6), "Kurtosis");
135+
Assert.That(stats.PopulationKurtosis, Is.EqualTo(kurtosisType1).Within(1e-6), "PopulationKurtosis");
136+
}
137+
138+
[Test]
139+
public void NegativeWeightsThrow()
140+
{
141+
Assert.That(() => new RunningWeightedStatistics(new[] { System.Tuple.Create(-1.0, 1.0) }), Throws.TypeOf<System.ArgumentOutOfRangeException>());
142+
var stats0 = new RunningWeightedStatistics(new System.Tuple<double, double>[0]);
143+
Assert.That(() => stats0.Push(-1.0, 1.0), Throws.TypeOf<System.ArgumentOutOfRangeException>());
144+
Assert.That(() => stats0.PushRange(new[] { System.Tuple.Create(-1.0, 1.0) }), Throws.TypeOf<System.ArgumentOutOfRangeException>());
145+
}
146+
147+
[Test]
148+
public void ShortSequences()
149+
{
150+
var stats0 = new RunningWeightedStatistics(new System.Tuple<double, double>[0]);
151+
Assert.That(stats0.Skewness, Is.NaN);
152+
Assert.That(stats0.Kurtosis, Is.NaN);
153+
154+
var stats1 = new RunningWeightedStatistics(new[] { System.Tuple.Create(1.0, 1.0) });
155+
Assert.That(stats1.Skewness, Is.NaN);
156+
Assert.That(stats1.Kurtosis, Is.NaN);
157+
158+
var stats2 = new RunningWeightedStatistics(new[] { System.Tuple.Create(1.0, 1.0), System.Tuple.Create(1.0, 2.0) });
159+
Assert.That(stats2.Skewness, Is.NaN);
160+
Assert.That(stats2.Kurtosis, Is.NaN);
161+
162+
var stats3 = new RunningWeightedStatistics(new[] { System.Tuple.Create(1.0, 1.0), System.Tuple.Create(1.0, 2.0), System.Tuple.Create(1.0, -3.0) });
163+
Assert.That(stats3.Skewness, Is.Not.NaN);
164+
Assert.That(stats3.Kurtosis, Is.NaN);
165+
166+
var stats4 = new RunningWeightedStatistics(new[] { System.Tuple.Create(1.0, 1.0), System.Tuple.Create(1.0, 2.0), System.Tuple.Create(1.0, -3.0), System.Tuple.Create(1.0, -4.0) });
167+
Assert.That(stats4.Skewness, Is.Not.NaN);
168+
Assert.That(stats4.Kurtosis, Is.Not.NaN);
169+
}
170+
171+
[Test]
172+
public void ZeroVarianceSequence()
173+
{
174+
var stats = new RunningWeightedStatistics(new[] { System.Tuple.Create(1.0, 2.0), System.Tuple.Create(1.0, 2.0), System.Tuple.Create(1.0, 2.0), System.Tuple.Create(1.0, 2.0) });
175+
Assert.That(stats.Skewness, Is.NaN);
176+
Assert.That(stats.Kurtosis, Is.NaN);
177+
}
178+
179+
[Test]
180+
public void CombineUnweighted()
181+
{
182+
var rnd = new SystemRandomSource(10);
183+
var a = Generate.Random(200, new Erlang(2, 0.2, rnd)).Select(datum => System.Tuple.Create(1.0, datum)).ToArray();
184+
var b = Generate.Random(100, new Beta(1.2, 1.4, rnd)).Select(datum => System.Tuple.Create(1.0, datum)).ToArray();
185+
var c = Generate.Random(150, new Rayleigh(0.8, rnd)).Select(datum => System.Tuple.Create(1.0, datum)).ToArray();
186+
187+
var d = a.Concat(b).Concat(c);
188+
var direct = d.Select(datum => datum.Item2).ToArray();
189+
190+
var x = new RunningWeightedStatistics(d);
191+
192+
var y = new RunningWeightedStatistics(a);
193+
y.PushRange(b);
194+
y.PushRange(c);
195+
196+
var za = new RunningWeightedStatistics(a);
197+
var zb = new RunningWeightedStatistics(b);
198+
var zc = new RunningWeightedStatistics(c);
199+
var z = za + zb + zc;
200+
201+
Assert.That(x.Mean, Is.EqualTo(direct.Mean()).Within(1e-12), "Mean Reference");
202+
Assert.That(y.Mean, Is.EqualTo(x.Mean).Within(1e-12), "Mean y");
203+
Assert.That(z.Mean, Is.EqualTo(x.Mean).Within(1e-12), "Mean z");
204+
205+
Assert.That(x.Variance, Is.EqualTo(direct.Variance()).Within(1e-12), "Variance Reference");
206+
Assert.That(y.Variance, Is.EqualTo(x.Variance).Within(1e-12), "Variance y");
207+
Assert.That(z.Variance, Is.EqualTo(x.Variance).Within(1e-12), "Variance z");
208+
209+
Assert.That(x.PopulationVariance, Is.EqualTo(direct.PopulationVariance()).Within(1e-12), "PopulationVariance Reference");
210+
Assert.That(y.PopulationVariance, Is.EqualTo(x.PopulationVariance).Within(1e-12), "PopulationVariance y");
211+
Assert.That(z.PopulationVariance, Is.EqualTo(x.PopulationVariance).Within(1e-12), "PopulationVariance z");
212+
213+
Assert.That(x.StandardDeviation, Is.EqualTo(direct.StandardDeviation()).Within(1e-12), "StandardDeviation Reference");
214+
Assert.That(y.StandardDeviation, Is.EqualTo(x.StandardDeviation).Within(1e-12), "StandardDeviation y");
215+
Assert.That(z.StandardDeviation, Is.EqualTo(x.StandardDeviation).Within(1e-12), "StandardDeviation z");
216+
217+
Assert.That(x.PopulationStandardDeviation, Is.EqualTo(direct.PopulationStandardDeviation()).Within(1e-12), "PopulationStandardDeviation Reference");
218+
Assert.That(y.PopulationStandardDeviation, Is.EqualTo(x.PopulationStandardDeviation).Within(1e-12), "PopulationStandardDeviation y");
219+
Assert.That(z.PopulationStandardDeviation, Is.EqualTo(x.PopulationStandardDeviation).Within(1e-12), "PopulationStandardDeviation z");
220+
221+
Assert.That(x.Skewness, Is.EqualTo(direct.Skewness()).Within(1e-12), "Skewness Reference (not independent!)");
222+
Assert.That(y.Skewness, Is.EqualTo(x.Skewness).Within(1e-12), "Skewness y");
223+
Assert.That(z.Skewness, Is.EqualTo(x.Skewness).Within(1e-12), "Skewness z");
224+
225+
Assert.That(x.PopulationSkewness, Is.EqualTo(direct.PopulationSkewness()).Within(1e-12), "PopulationSkewness Reference (not independent!)");
226+
Assert.That(y.PopulationSkewness, Is.EqualTo(x.PopulationSkewness).Within(1e-12), "PopulationSkewness y");
227+
Assert.That(z.PopulationSkewness, Is.EqualTo(x.PopulationSkewness).Within(1e-12), "PopulationSkewness z");
228+
229+
Assert.That(x.Kurtosis, Is.EqualTo(direct.Kurtosis()).Within(1e-12), "Kurtosis Reference (not independent!)");
230+
Assert.That(y.Kurtosis, Is.EqualTo(x.Kurtosis).Within(1e-12), "Kurtosis y");
231+
Assert.That(z.Kurtosis, Is.EqualTo(x.Kurtosis).Within(1e-12), "Kurtosis z");
232+
233+
Assert.That(x.PopulationKurtosis, Is.EqualTo(direct.PopulationKurtosis()).Within(1e-12), "PopulationKurtosis Reference (not independent!)");
234+
Assert.That(y.PopulationKurtosis, Is.EqualTo(x.PopulationKurtosis).Within(1e-12), "PopulationKurtosis y");
235+
Assert.That(z.PopulationKurtosis, Is.EqualTo(x.PopulationKurtosis).Within(1e-12), "PopulationKurtosis z");
236+
}
237+
238+
[Test]
239+
/// Tests that combination of data via + / Combine is consistent with the incremental approach.
240+
public void CombineWeighted()
241+
{
242+
var rnd = new SystemRandomSource(10);
243+
var wa = Generate.Random(200, new ContinuousUniform(1.0, 10.0));
244+
var a = Generate.Random(200, new Erlang(2, 0.2, rnd)).Select((datum, i) => System.Tuple.Create(wa[i], datum)).ToArray();
245+
var wb = Generate.Random(100, new ContinuousUniform(1.0, 10.0));
246+
var b = Generate.Random(100, new Beta(1.2, 1.4, rnd)).Select((datum, i) => System.Tuple.Create(wb[i], datum)).ToArray();
247+
var wc = Generate.Random(150, new ContinuousUniform(1.0, 10.0));
248+
var c = Generate.Random(150, new Rayleigh(0.8, rnd)).Select((datum, i) => System.Tuple.Create(wc[i], datum)).ToArray();
249+
250+
var d = a.Concat(b).Concat(c);
251+
252+
var x = new RunningWeightedStatistics(d);
253+
254+
var y = new RunningWeightedStatistics(a);
255+
y.PushRange(b);
256+
y.PushRange(c);
257+
258+
var za = new RunningWeightedStatistics(a);
259+
var zb = new RunningWeightedStatistics(b);
260+
var zc = new RunningWeightedStatistics(c);
261+
var z = za + zb + zc;
262+
263+
Assert.That(y.Mean, Is.EqualTo(x.Mean).Within(1e-12), "Mean y");
264+
Assert.That(z.Mean, Is.EqualTo(x.Mean).Within(1e-12), "Mean z");
265+
266+
Assert.That(y.Variance, Is.EqualTo(x.Variance).Within(1e-12), "Variance y");
267+
Assert.That(z.Variance, Is.EqualTo(x.Variance).Within(1e-12), "Variance z");
268+
269+
Assert.That(y.PopulationVariance, Is.EqualTo(x.PopulationVariance).Within(1e-12), "PopulationVariance y");
270+
Assert.That(z.PopulationVariance, Is.EqualTo(x.PopulationVariance).Within(1e-12), "PopulationVariance z");
271+
272+
Assert.That(y.StandardDeviation, Is.EqualTo(x.StandardDeviation).Within(1e-12), "StandardDeviation y");
273+
Assert.That(z.StandardDeviation, Is.EqualTo(x.StandardDeviation).Within(1e-12), "StandardDeviation z");
274+
275+
Assert.That(y.PopulationStandardDeviation, Is.EqualTo(x.PopulationStandardDeviation).Within(1e-12), "PopulationStandardDeviation y");
276+
Assert.That(z.PopulationStandardDeviation, Is.EqualTo(x.PopulationStandardDeviation).Within(1e-12), "PopulationStandardDeviation z");
277+
278+
Assert.That(y.Skewness, Is.EqualTo(x.Skewness).Within(1e-12), "Skewness y");
279+
Assert.That(z.Skewness, Is.EqualTo(x.Skewness).Within(1e-12), "Skewness z");
280+
281+
Assert.That(y.PopulationSkewness, Is.EqualTo(x.PopulationSkewness).Within(1e-12), "PopulationSkewness y");
282+
Assert.That(z.PopulationSkewness, Is.EqualTo(x.PopulationSkewness).Within(1e-12), "PopulationSkewness z");
283+
284+
Assert.That(y.Kurtosis, Is.EqualTo(x.Kurtosis).Within(1e-12), "Kurtosis y");
285+
Assert.That(z.Kurtosis, Is.EqualTo(x.Kurtosis).Within(1e-12), "Kurtosis z");
286+
287+
Assert.That(y.PopulationKurtosis, Is.EqualTo(x.PopulationKurtosis).Within(1e-12), "PopulationKurtosis y");
288+
Assert.That(z.PopulationKurtosis, Is.EqualTo(x.PopulationKurtosis).Within(1e-12), "PopulationKurtosis z");
289+
}
290+
291+
[TestCase("lottery")]
292+
[TestCase("lew")]
293+
[TestCase("mavro")]
294+
[TestCase("michelso")]
295+
[TestCase("numacc1")]
296+
[TestCase("meixner")]
297+
/// Generates samples with weightings that are integral and compares that to the unweighted statistics result. Doesn't correspond with the
298+
/// higher order sample statistics because our weightings represent reliability weights, *not* frequency weights, and the Bessel correction is
299+
/// calculated appropriately - so don't let the construction of the test mislead you.
300+
public void ConsistentWithUnweighted (string dataSet)
301+
{
302+
var data = _data[dataSet].Data.ToArray();
303+
var gen = new DiscreteUniform(1, 5);
304+
var weights = new int[data.Length];
305+
gen.Samples(weights);
306+
307+
var stats = new RunningWeightedStatistics( data.Select((x, i) => System.Tuple.Create((double)weights[i], x)) );
308+
var stats2 = new RunningStatistics();
309+
for (int i = 0; i < data.Length; ++i)
310+
for(int j = 0; j < weights[i] ; ++j)
311+
stats2.Push(data[i]);
312+
var sumWeights = weights.Sum();
313+
Assert.That(stats.TotalWeight, Is.EqualTo(sumWeights), "TotalWeight");
314+
Assert.That(stats.Count, Is.EqualTo(weights.Length), "Count");
315+
Assert.That(stats2.Minimum, Is.EqualTo(stats.Minimum), "Minimum");
316+
Assert.That(stats2.Maximum, Is.EqualTo(stats.Maximum), "Maximum");
317+
Assert.That(stats2.Mean, Is.EqualTo(stats.Mean).Within(1e-8), "Mean");
318+
Assert.That(stats2.PopulationVariance, Is.EqualTo(stats.PopulationVariance).Within(1e-9), "PopulationVariance");
319+
Assert.That(stats2.PopulationStandardDeviation, Is.EqualTo(stats.PopulationStandardDeviation).Within(1e-9), "PopulationStandardDeviation");
320+
Assert.That(stats2.PopulationSkewness, Is.EqualTo(stats.PopulationSkewness).Within(1e-8), "PopulationSkewness");
321+
Assert.That(stats2.PopulationKurtosis, Is.EqualTo(stats.PopulationKurtosis).Within(1e-8), "PopulationKurtosis");
322+
}
323+
}
324+
}

0 commit comments

Comments
 (0)