View Javadoc
1   package org.djutils.stats.summarizers;
2   
3   import static org.junit.Assert.assertEquals;
4   import static org.junit.Assert.assertFalse;
5   import static org.junit.Assert.assertNull;
6   import static org.junit.Assert.assertTrue;
7   import static org.junit.Assert.fail;
8   
9   import java.util.ArrayList;
10  import java.util.List;
11  import java.util.Random;
12  
13  import org.djutils.stats.ConfidenceInterval;
14  import org.djutils.stats.DistNormalTable;
15  import org.djutils.stats.summarizers.quantileaccumulator.FullStorageAccumulator;
16  import org.djutils.stats.summarizers.quantileaccumulator.NoStorageAccumulator;
17  import org.djutils.stats.summarizers.quantileaccumulator.TDigestAccumulator;
18  import org.junit.Test;
19  
20  /**
21   * The TallyTest test the tally.
22   * <p>
23   * Copyright (c) 2002-2022 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
24   * for project information <a href="https://simulation.tudelft.nl/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
25   * project is distributed under a three-clause BSD-style license, which can be found at
26   * <a href="https://simulation.tudelft.nl/dsol/3.0/license.html" target="_blank">
27   * https://simulation.tudelft.nl/dsol/3.0/license.html</a>. <br>
28   * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
29   * @author <a href="https://www.tudelft.nl/staff/p.knoppers/">Peter Knoppers</a>
30   */
31  public class TallyTest
32  {
33      /** Test the tally. */
34      @SuppressWarnings("checkstyle:methodlength")
35      @Test
36      public void testTally()
37      {
38          String description = "THIS TALLY IS TESTED";
39          Tally tally = new Tally(description);
40  
41          // check the description
42          assertTrue(tally.toString().contains(description));
43          assertEquals(description, tally.getDescription());
44          assertTrue(tally.toString().startsWith("Tally"));
45  
46          // now we check the initial values
47          assertTrue(Double.valueOf(tally.getMin()).isNaN());
48          assertTrue(Double.valueOf(tally.getMax()).isNaN());
49          assertTrue(Double.valueOf(tally.getSampleMean()).isNaN());
50          assertTrue(Double.valueOf(tally.getPopulationMean()).isNaN());
51          assertTrue(Double.valueOf(tally.getSampleVariance()).isNaN());
52          assertTrue(Double.valueOf(tally.getPopulationVariance()).isNaN());
53          assertTrue(Double.valueOf(tally.getSampleStDev()).isNaN());
54          assertTrue(Double.valueOf(tally.getPopulationStDev()).isNaN());
55          assertTrue(Double.valueOf(tally.getSampleSkewness()).isNaN());
56          assertTrue(Double.valueOf(tally.getPopulationSkewness()).isNaN());
57          assertTrue(Double.valueOf(tally.getSampleKurtosis()).isNaN());
58          assertTrue(Double.valueOf(tally.getPopulationKurtosis()).isNaN());
59          assertTrue(Double.valueOf(tally.getSampleExcessKurtosis()).isNaN());
60          assertTrue(Double.valueOf(tally.getPopulationExcessKurtosis()).isNaN());
61          assertEquals(0, tally.getSum(), 0);
62          assertEquals(0L, tally.getN());
63          assertNull(tally.getConfidenceInterval(0.95));
64          assertNull(tally.getConfidenceInterval(0.95, ConfidenceInterval.LEFT_SIDE_CONFIDENCE));
65          assertNull(tally.getConfidenceInterval(0.95, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE));
66          assertNull(tally.getConfidenceInterval(0.95, ConfidenceInterval.BOTH_SIDE_CONFIDENCE));
67  
68          // Now we ingest some values
69          try
70          {
71              tally.register(1.1);
72              assertFalse("sample mean is now available", Double.isNaN(tally.getSampleMean()));
73              assertFalse("mean is now available", Double.isNaN(tally.getPopulationMean()));
74              assertEquals("smaple mean is 1.1", 1.1, tally.getSampleMean(), 0.0000001);
75              assertEquals("mean is 1.1", 1.1, tally.getPopulationMean(), 0.0000001);
76              assertTrue("sample variance is not available", Double.isNaN(tally.getSampleVariance()));
77              assertFalse("variance is not available", Double.isNaN(tally.getPopulationVariance()));
78              assertTrue("skewness is not available", Double.isNaN(tally.getPopulationSkewness()));
79              tally.register(1.2);
80              assertFalse("sample variance is now available", Double.isNaN(tally.getSampleVariance()));
81              assertTrue("sample skewness is not available", Double.isNaN(tally.getSampleSkewness()));
82              assertTrue("sample kurtosis is not available", Double.isNaN(tally.getSampleKurtosis()));
83              assertFalse("skewness is available", Double.isNaN(tally.getPopulationSkewness()));
84              assertTrue("kurtosis is not available", Double.isNaN(tally.getPopulationKurtosis()));
85              tally.register(1.3);
86              assertFalse("skewness is now available", Double.isNaN(tally.getSampleSkewness()));
87              assertFalse("kurtosis is now available", Double.isNaN(tally.getPopulationKurtosis()));
88              assertTrue("sample kurtosis is not available", Double.isNaN(tally.getSampleKurtosis()));
89              tally.register(1.4);
90              assertFalse("sample kurtosis is now available", Double.isNaN(tally.getSampleKurtosis()));
91              tally.register(1.5);
92              tally.register(1.6);
93              tally.register(1.7);
94              tally.register(1.8);
95              tally.register(1.9);
96              tally.register(2.0);
97              tally.register(1.0);
98          }
99          catch (Exception exception)
100         {
101             fail(exception.getMessage());
102         }
103 
104         // Now we check the tally
105         assertEquals(2.0, tally.getMax(), 1.0E-6);
106         assertEquals(1.0, tally.getMin(), 1.0E-6);
107         assertEquals(11, tally.getN());
108         assertEquals(16.5, tally.getSum(), 1.0E-6);
109         assertEquals(1.5, tally.getSampleMean(), 1.0E-6);
110         assertEquals(0.110000, tally.getSampleVariance(), 1.0E-6);
111         assertEquals(0.331662, tally.getSampleStDev(), 1.0E-6);
112         // From: https://atozmath.com/StatsUG.aspx
113         assertEquals(0.1, tally.getPopulationVariance(), 1.0E-6);
114         assertEquals(Math.sqrt(0.1), tally.getPopulationStDev(), 1.0E-6);
115         assertEquals(0.0, tally.getPopulationSkewness(), 1.0E-6);
116         assertEquals(1.78, tally.getPopulationKurtosis(), 1.0E-6);
117         assertEquals(-1.22, tally.getPopulationExcessKurtosis(), 1.0E-6);
118         assertEquals(0.0, tally.getSampleSkewness(), 1.0E-6);
119         assertEquals(-1.2, tally.getSampleExcessKurtosis(), 1.0E-6);
120         assertEquals(1.618182, tally.getSampleKurtosis(), 1.0E-6);
121 
122         assertEquals(1.304003602, tally.getConfidenceInterval(0.05)[0], 1E-05);
123         assertEquals(1.695996398, tally.getConfidenceInterval(0.05)[1], 1E-05);
124         assertEquals(1.335514637, tally.getConfidenceInterval(0.10)[0], 1E-05);
125         assertEquals(1.664485363, tally.getConfidenceInterval(0.10)[1], 1E-05);
126         assertEquals(1.356046853, tally.getConfidenceInterval(0.15)[0], 1E-05);
127         assertEquals(1.643953147, tally.getConfidenceInterval(0.15)[1], 1E-05);
128         assertEquals(1.432551025, tally.getConfidenceInterval(0.50)[0], 1E-05);
129         assertEquals(1.567448975, tally.getConfidenceInterval(0.50)[1], 1E-05);
130         assertEquals(1.474665290, tally.getConfidenceInterval(0.80)[0], 1E-05);
131         assertEquals(1.525334710, tally.getConfidenceInterval(0.80)[1], 1E-05);
132         assertEquals(1.493729322, tally.getConfidenceInterval(0.95)[0], 1E-05);
133         assertEquals(1.506270678, tally.getConfidenceInterval(0.95)[1], 1E-05);
134 
135         assertEquals(1.304003602, tally.getConfidenceInterval(0.05, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[0], 1E-05);
136         assertEquals(1.695996398, tally.getConfidenceInterval(0.05, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[1], 1E-05);
137         assertEquals(1.432551025, tally.getConfidenceInterval(0.50, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[0], 1E-05);
138         assertEquals(1.567448975, tally.getConfidenceInterval(0.50, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[1], 1E-05);
139         assertEquals(1.493729322, tally.getConfidenceInterval(0.95, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[0], 1E-05);
140         assertEquals(1.506270678, tally.getConfidenceInterval(0.95, ConfidenceInterval.BOTH_SIDE_CONFIDENCE)[1], 1E-05);
141 
142         assertEquals(1.304003602, tally.getConfidenceInterval(0.025, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[0], 1E-05);
143         assertEquals(1.500000000, tally.getConfidenceInterval(0.025, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[1], 1E-05);
144         assertEquals(1.432551025, tally.getConfidenceInterval(0.25, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[0], 1E-05);
145         assertEquals(1.500000000, tally.getConfidenceInterval(0.25, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[1], 1E-05);
146         assertEquals(1.474665290, tally.getConfidenceInterval(0.40, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[0], 1E-05);
147         assertEquals(1.500000000, tally.getConfidenceInterval(0.40, ConfidenceInterval.LEFT_SIDE_CONFIDENCE)[1], 1E-05);
148 
149         assertEquals(1.500000000, tally.getConfidenceInterval(0.025, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[0], 1E-05);
150         assertEquals(1.695996398, tally.getConfidenceInterval(0.025, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[1], 1E-05);
151         assertEquals(1.500000000, tally.getConfidenceInterval(0.25, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[0], 1E-05);
152         assertEquals(1.567448975, tally.getConfidenceInterval(0.25, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[1], 1E-05);
153         assertEquals(1.500000000, tally.getConfidenceInterval(0.40, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[0], 1E-05);
154         assertEquals(1.525334710, tally.getConfidenceInterval(0.40, ConfidenceInterval.RIGHT_SIDE_CONFIDENCE)[1], 1E-05);
155 
156         // we check the input of the confidence interval
157         try
158         {
159             tally.getConfidenceInterval(0.95, null);
160             fail("null is not defined as side of confidence level");
161         }
162         catch (Exception exception)
163         {
164             assertTrue(exception.getClass().equals(NullPointerException.class));
165         }
166         try
167         {
168             assertNull(tally.getConfidenceInterval(-0.95));
169             fail("should have reacted on wrong confidence level -0.95");
170         }
171         catch (Exception exception)
172         {
173             assertTrue(exception.getClass().equals(IllegalArgumentException.class));
174         }
175         try
176         {
177             assertNull(tally.getConfidenceInterval(1.14));
178             fail("should have reacted on wrong confidence level 1.14");
179         }
180         catch (Exception exception)
181         {
182             assertTrue(exception.getClass().equals(IllegalArgumentException.class));
183         }
184 
185         assertTrue(Math.abs(tally.getSampleMean() - 1.5) < 10E-6);
186 
187         // Let's compute the standard deviation
188         double varianceAccumulator = 0;
189         for (int i = 0; i < 11; i++)
190         {
191             varianceAccumulator = Math.pow(1.5 - (1.0 + i / 10.0), 2) + varianceAccumulator;
192         }
193 
194         assertEquals(varianceAccumulator / 10.0, tally.getSampleVariance(), 1.0E-6);
195         assertEquals(Math.sqrt(varianceAccumulator / 10.0), tally.getSampleStDev(), 1.0E-6);
196 
197         assertEquals(varianceAccumulator / 11.0, tally.getPopulationVariance(), 1.0E-6);
198         assertEquals(Math.sqrt(varianceAccumulator / 11.0), tally.getPopulationStDev(), 1.0E-6);
199 
200         // test confidence interval failure on uninitialized tally
201         Tally t = new Tally("unused");
202         assertNull(t.getConfidenceInterval(0.95));
203         t.register(1.0);
204         assertNull(t.getConfidenceInterval(0.95));
205     }
206 
207     /**
208      * Test the kurtosis example from https://en.wikipedia.org/wiki/Kurtosis where we assumed that the example for sample
209      * kurtosis actually calculates population kurtosis (!).
210      */
211     @Test
212     public void testKurtWikipedia()
213     {
214         Tally tally = new Tally("Wikipedia");
215         tally.register(0, 3, 4, 1, 2, 3, 0, 2, 1, 3, 2, 0, 2, 2, 3, 2, 5, 2, 3, 999);
216         assertEquals(18.05, tally.getPopulationKurtosis(), 0.01);
217         assertEquals(15.05, tally.getPopulationExcessKurtosis(), 0.01);
218     }
219 
220     /**
221      * Test skewness and kurtosis based on two Excel samples.
222      */
223     @Test
224     public void testSkewKurtExcel()
225     {
226         Tally tally1 = new Tally("Excel1");
227         tally1.register(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2);
228         assertEquals(-1.2, tally1.getSampleExcessKurtosis(), 0.01);
229         Tally tally2 = new Tally("Excel2");
230         tally2.register(2, 4, 6, 3, 2, 1, 2, 3, 4, 5, 9);
231         assertEquals(3.7272, tally2.getPopulationMean(), 0.01);
232         assertEquals(4.7438, tally2.getPopulationVariance(), 0.01);
233         assertEquals(5.2182, tally2.getSampleVariance(), 0.01);
234         assertEquals(2.1780, tally2.getPopulationStDev(), 0.01);
235         assertEquals(2.2843, tally2.getSampleStDev(), 0.01);
236         assertEquals(1.0903, tally2.getPopulationSkewness(), 0.01);
237         assertEquals(1.2706, tally2.getSampleSkewness(), 0.01);
238         assertEquals(1.7908, tally2.getSampleExcessKurtosis(), 0.01);
239     }
240 
241     /**
242      * Test Tally with the NoStorageAccumulator.
243      */
244     @Test
245     public void testNoStorageAccumulator()
246     {
247         Tally tally = new Tally("test with the NoStorageAccumulator", new NoStorageAccumulator());
248         assertTrue("mean of no data is NaN", Double.isNaN(tally.getSampleMean()));
249         try
250         {
251             tally.getQuantile(0.5);
252             fail("getQuantile of no data should have resulted in an IllegalArgumentException");
253         }
254         catch (IllegalArgumentException iae)
255         {
256             // Ignore expected exception
257         }
258 
259         tally.register(90.0);
260         assertEquals("mean of one value is that value", 90.0, tally.getSampleMean(), 0);
261         try
262         {
263             tally.getQuantile(0.5);
264             fail("getQuantile of one value should have resulted in an IllegalArgumentException");
265         }
266         catch (IllegalArgumentException iae)
267         {
268             // Ignore expected exception
269         }
270 
271         tally.register(110.0);
272         assertEquals("mean of two value", 100.0, tally.getSampleMean(), 0);
273         assertEquals("50% quantile", 100.0, tally.getQuantile(0.5), 0);
274         /*-
275         double sigma = tally.getSampleStDev();
276         double mu = tally.getSampleMean();
277         // For loop below makes painfully clear where the getQuantile method fails
278         // Values are from last table in https://en.wikipedia.org/wiki/Standard_normal_table
279         for (double probability : new double[] { 1 - 5.00000E-1, 1 - 1.58655E-1, 1 - 2.27501E-2, 1 - 1.34990E-3, 1 - 3.16712E-5,
280                 1 - 2.86652E-7, 1 - 9.86588E-10, 1 - 1.27981E-12, 1 - 6.22096E-16, 1 - 1.12859E-19, 1 - 7.61985E-24 })
281         {
282             double x = tally.getQuantile(probability);
283             System.out.println(String.format("probability=%19.16f 1-probability=%19.16f, x=%19.14f, sigmaCount=%19.16f",
284                     probability, 1 - probability, x, (x - mu) / sigma));
285         }
286         // Output shows that the inverse cumulative probability function works fine up to about 8 sigma
287          */
288 
289         assertEquals("84% is about one sigma", 1, DistNormalTable.getInverseCumulativeProbability(0, 1, 0.84), 0.01);
290         assertEquals("16% is about minus one sigma", -1, DistNormalTable.getInverseCumulativeProbability(0, 1, 0.16), 0.01);
291 
292         // Test for the problem that Peter Knoppers had in Tritapt where really small rounding errors caused sqrt(-1e-14).
293         double value = 166.0 / 25.0;
294         tally.initialize();
295         tally.register(value);
296         tally.register(value);
297         tally.register(value);
298         tally.register(value);
299         tally.initialize();
300         // Throw a lot of pseudo-randomly normally distributed values in and see if the expected mean and stddev come out
301         double mean = 123.456;
302         double stddev = 234.567;
303         Random random = new Random(123456);
304         for (int sample = 0; sample < 10000; sample++)
305         {
306             value = generateGaussianNoise(mean, stddev, random);
307             tally.register(value);
308         }
309         assertEquals("mean should approximately match", mean, tally.getSampleMean(), stddev / 10);
310         assertEquals("stddev should approximately match", stddev, tally.getSampleStDev(), stddev / 10);
311     }
312 
313     /**
314      * Test the FullStorageAccumulator.
315      */
316     @Test
317     public void testFullStorageAccumulator()
318     {
319         Tally tally = new Tally("Tally for FullStorageAccumulator test", new FullStorageAccumulator());
320         // Insert values from 0.0 .. 100.0 (step 1.0)
321         for (int step = 0; step <= 100; step++)
322         {
323             tally.register(1.0 * step);
324         }
325         for (double probability : new double[] {0.0, 0.01, 0.1, 0.49, 0.5, 0.51, 0.9, 0.99, 1.0})
326         {
327             double expected = 100 * probability;
328             double got = tally.getQuantile(probability);
329             assertEquals("quantile should match", expected, got, 0.00001);
330         }
331         try
332         {
333             tally.getQuantile(-0.01);
334             fail("negative probability should have thrown an exception");
335         }
336         catch (IllegalArgumentException iae)
337         {
338             // Ignore expected exception
339         }
340 
341         try
342         {
343             tally.getQuantile(1.01);
344             fail("Probability > 1 should have thrown an exception");
345         }
346         catch (IllegalArgumentException iae)
347         {
348             // Ignore expected exception
349         }
350 
351         assertTrue("toString returns something descriptive",
352                 new FullStorageAccumulator().toString().startsWith("FullStorageAccumulator"));
353     }
354 
355     /**
356      * Generate normally distributed values. Derived from https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
357      * @param mu double; mean
358      * @param sigma double; standard deviation
359      * @param random Random; entropy source
360      * @return double; one pseudo random value
361      */
362     double generateGaussianNoise(final double mu, final double sigma, final Random random)
363     {
364         final double epsilon = Double.MIN_VALUE;
365         final double twoPi = Math.PI * 2;
366 
367         double u1, u2;
368         do
369         {
370             u1 = random.nextDouble();
371             u2 = random.nextDouble();
372         }
373         while (u1 <= epsilon);
374 
375         return mu + sigma * Math.sqrt(-2.0 * Math.log(u1)) * Math.cos(twoPi * u2);
376     }
377 
378     /**
379      * Test the T-Digest accumulator.
380      */
381     @Test
382     public void testTDigestAccumulator()
383     {
384         Tally tally = new Tally("Tally for TDigestAccumulator test", new TDigestAccumulator());
385         // Insert values from 0.0 .. 100.0 (step 1.0)
386         for (int step = 0; step <= 100; step++)
387         {
388             tally.register(1.0 * step);
389         }
390         for (double probability : new double[] {0.0, 0.01, 0.1, 0.49, 0.5, 0.51, 0.9, 0.99, 1.0})
391         {
392             double expected = 100 * probability;
393             double got = tally.getQuantile(probability);
394             // System.out.println(String.format("probability %10.8f, expected %10.8f, got %10.8f", probability, expected, got));
395             assertEquals("quantile should match", expected, got, 1.0); // With 100 bins the error should be below 1%
396         }
397         // https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)
398         assertEquals("sample skewness should be 0", 0, tally.getSampleSkewness(), 0.0001);
399         assertEquals("skewness should be 0", 0, tally.getPopulationSkewness(), 0.0001);
400         assertEquals("sample excess kurtosis should be -1.2", -1.2, tally.getSampleExcessKurtosis(), 0.0001);
401         assertEquals("population excess kurtosis should be -1.2", -1.2, tally.getPopulationExcessKurtosis(), 0.01);
402         // System.out.println(String.format("%d uniformly distributed values: skewness %20.15f, kurtosis %20.15f", tally.getN(),
403         // tally.getSampleSkewness(), tally.getSampleKurtosis()));
404         tally.initialize();
405         // Insert values from 0.0 .. 100.0 (step 0.0001)
406         for (int step = 0; step <= 1000000; step++)
407         {
408             tally.register(0.0001 * step);
409         }
410         for (double probability : new double[] {0.0, 0.01, 0.1, 0.49, 0.5, 0.51, 0.9, 0.99, 1.0})
411         {
412             double expected = 100 * probability;
413             double got = tally.getQuantile(probability);
414             // System.out.println(String.format("probability %10.8f, expected %10.8f, got %10.8f", probability, expected, got));
415             assertEquals("quantile should match", expected, got, 0.01); // Uniformly distributed data yields very good estimates
416         }
417         // System.out.println(String.format("%d uniformly distributed values: skewness %20.15f, kurtosis %20.15f", tally.getN(),
418         // tally.getSampleSkewness(), tally.getSampleKurtosis()));
419         assertEquals("sample skewness should be 0", 0, tally.getSampleSkewness(), 0.0001);
420         assertEquals("skewness should be 0", 0, tally.getPopulationSkewness(), 0.0001);
421         assertEquals("sample kurtosis should be -1.2", -1.2, tally.getSampleExcessKurtosis(), 0.0001);
422         assertEquals("population excess kurtosis should be -1.2", -1.2, tally.getPopulationExcessKurtosis(), 0.0001);
423         tally = new Tally("Tally for TDigestAccumulator test", new TDigestAccumulator(4));
424         // Insert values from 0.0 .. 100.0 (step 0.0001)
425         for (int step = 0; step <= 1000000; step++)
426         {
427             tally.register(0.0001 * step);
428         }
429         for (double probability : new double[] {0.0, 0.01, 0.1, 0.49, 0.5, 0.51, 0.9, 0.99, 1.0})
430         {
431             double expected = 100 * probability;
432             double got = tally.getQuantile(probability);
433             // System.out.println(String.format("probability %10.8f, expected %10.8f, got %10.8f", probability, expected, got));
434             assertEquals("quantile should match", expected, got, 0.01); // Uniformly distributed data yields very good estimates
435         }
436         try
437         {
438             tally.getQuantile(-0.01);
439             fail("negative probability should have thrown an exception");
440         }
441         catch (IllegalArgumentException iae)
442         {
443             // Ignore expected exception
444         }
445 
446         try
447         {
448             tally.getQuantile(1.01);
449             fail("Probability > 1 should have thrown an exception");
450         }
451         catch (IllegalArgumentException iae)
452         {
453             // Ignore expected exception
454         }
455         assertTrue("toString returns something descriptive",
456                 new TDigestAccumulator().toString().startsWith("TDigestAccumulator"));
457         tally = new Tally("Tally for TDigestAccumulator test", new TDigestAccumulator());
458         // Throw a lot of pseudo-randomly normally distributed values in and see if the accumulator gets a good approximation
459         // of the distribution
460         double mean = 123.456;
461         double stddev = 234.567;
462         Random random = new Random(123456);
463         for (int sample = 0; sample < 10000; sample++)
464         {
465             double value = generateGaussianNoise(mean, stddev, random);
466             tally.register(value);
467         }
468         // Test that tally reports cumulative probabilities that roughly follow that of normal distribution
469         for (double probability : new double[] {0.01, 0.1, 0.25, 0.49, 0.5, 0.51, 0.75, 0.9, 0.99})
470         {
471             double expected = DistNormalTable.getInverseCumulativeProbability(mean, stddev, probability);
472             double got = tally.getQuantile(probability);
473             double margin = mean / 10 / Math.sqrt(Math.min(probability, 1 - probability));
474             // System.out.println(String.format("probability %12.7f, expected %12.7f, reasonable margin %12.7f, got %12.7f",
475             // probability, expected, margin, got));
476             assertEquals("quantile should match", expected, got, margin);
477         }
478     }
479 
480     /**
481      * Test skewness and kurtosis. Test data from http://web.ipac.caltech.edu/staff/fmasci/home/astro_refs/SkewStatSignif.pdf
482      */
483     @Test
484     public void testSkewnessAndKurtosis()
485     {
486         List<Double> testValues = new ArrayList<>();
487         for (int i = 5; --i >= 0;)
488         {
489             testValues.add(61.0);
490         }
491         for (int i = 18; --i >= 0;)
492         {
493             testValues.add(64.0);
494         }
495         for (int i = 42; --i >= 0;)
496         {
497             testValues.add(67.0);
498         }
499         for (int i = 27; --i >= 0;)
500         {
501             testValues.add(70.0);
502         }
503         for (int i = 8; --i >= 0;)
504         {
505             testValues.add(73.0);
506         }
507         int count = testValues.size();
508         Tally tally = new Tally("");
509         for (double value : testValues)
510         {
511             tally.register(value);
512         }
513         // System.out.println(tally);
514         // System.out.println(String.format("count %d mean %20.15f variance %20.15f skew %20.15f kurtosis %20.15f", count,
515         // tally.getSampleMean(), tally.getSampleVariance(), tally.getSampleSkewness(), tally.getSampleKurtosis()));
516         // Do the math the "classic" way (i.e. using two passes; the first pass gets the mean)
517         double mean = tally.getSampleMean();
518         double m2 = 0;
519         double m3 = 0;
520         double m4 = 0;
521         for (double value : testValues)
522         {
523             double delta = value - mean;
524             m2 += delta * delta;
525             m3 += delta * delta * delta;
526             m4 += delta * delta * delta * delta;
527         }
528         m2 /= count;
529         m3 /= count;
530         m4 /= count;
531         double g1 = m3 / Math.pow(m2, 1.5);
532         double sg1 = g1 * Math.sqrt(count * (count - 1)) / (count - 2);
533         double a4 = m4 / m2 / m2;
534         // System.out.println(String.format("m2 %20.15f, m3 %20.15f, m4 %20.15f", m2, m3, m4));
535         double g2 = a4 - 3;
536         double sg2 = 1.0 * (count - 1) / (count - 2) / (count - 3) * ((count + 1) * g2 + 6);
537         // System.out.println(String.format("g1 %20.15f sampleSkewness %20.15f, a4 %20.15f g2 %20.15f sampleKurtosis %20.15f",
538         // g1, sg1, a4, g2, sg2));
539         assertEquals("skew should match", sg1, tally.getSampleSkewness(), 0.0001);
540         assertEquals("kurtosis should match", sg2, tally.getSampleExcessKurtosis(), 0.0001);
541     }
542 
543 }