View Javadoc
1   package org.djutils.stats.summarizers;
2   
3   import org.djutils.exceptions.Throw;
4   import org.djutils.stats.ConfidenceInterval;
5   import org.djutils.stats.DistNormalTable;
6   import org.djutils.stats.summarizers.quantileaccumulator.NoStorageAccumulator;
7   import org.djutils.stats.summarizers.quantileaccumulator.QuantileAccumulator;
8   
9   /**
10   * The Tally class ingests a series of values and provides mean, standard deviation, etc. of the ingested values.
11   * <p>
12   * Copyright (c) 2002-2021 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
13   * for project information <a href="https://simulation.tudelft.nl/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
14   * project is distributed under a three-clause BSD-style license, which can be found at
15   * <a href="https://simulation.tudelft.nl/dsol/3.0/license.html" target="_blank">
16   * https://simulation.tudelft.nl/dsol/3.0/license.html</a>. <br>
17   * @author <a href="https://www.tudelft.nl/averbraeck" target="_blank"> Alexander Verbraeck</a>
18   * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
19   * @author <a href="https://www.tudelft.nl/staff/p.knoppers/">Peter Knoppers</a>
20   */
21  public class Tally implements TallyInterface
22  {
23      /** */
24      private static final long serialVersionUID = 20200228L;
25  
26      /** The sum of this tally. */
27      private double sum = 0;
28  
29      /** The mean of this tally. */
30      private double m1 = 0;
31  
32      /** The summation for the second moment (variance). */
33      private double m2 = 0;
34  
35      /** The summation for the third moment (skewness). */
36      private double m3 = 0;
37  
38      /** The summation for the fourth moment (kurtosis). */
39      private double m4 = 0;
40  
41      /** The minimum observed value of this tally. */
42      private double min = Double.NaN;
43  
44      /** The maximum observed value of this tally. */
45      private double max = Double.NaN;
46  
47      /** The number of measurements of this tally. */
48      private long n = 0;
49  
50      /** The description of this tally. */
51      private final String description;
52  
53      /** The quantile accumulator. */
54      private final QuantileAccumulator quantileAccumulator;
55  
56      /** the synchronized lock. */
57      @SuppressWarnings("checkstyle:visibilitymodifier")
58      protected Object semaphore = new Object();
59  
60      /**
61       * Constructs a new Tally.
62       * @param description String; the description of this tally
63       * @param quantileAccumulator QuantileAccumulator; the input series accumulator that can approximate or compute quantiles.
64       */
65      public Tally(final String description, final QuantileAccumulator quantileAccumulator)
66      {
67          this.description = description;
68          this.quantileAccumulator = quantileAccumulator;
69          initialize();
70      }
71  
72      /**
73       * Convenience constructor that uses a NoStorageAccumulator to estimate quantiles.
74       * @param description String; the description of this tally
75       */
76      public Tally(final String description)
77      {
78          this(description, new NoStorageAccumulator());
79      }
80  
81      /** {@inheritDoc} */
82      @Override
83      public final double getSampleMean()
84      {
85          if (this.n > 0)
86          {
87              return this.m1;
88          }
89          return Double.NaN;
90      }
91  
92      /** {@inheritDoc} */
93      @Override
94      public final double getQuantile(final double probability)
95      {
96          return this.quantileAccumulator.getQuantile(this, probability);
97      }
98  
99      /** {@inheritDoc} */
100     @Override
101     public double getCumulativeProbability(final double quantile)
102     {
103         return this.quantileAccumulator.getCumulativeProbability(this, quantile);
104     }
105 
106     /** {@inheritDoc} */
107     @Override
108     public final double[] getConfidenceInterval(final double alpha)
109     {
110         return this.getConfidenceInterval(alpha, ConfidenceInterval.BOTH_SIDE_CONFIDENCE);
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     public final double[] getConfidenceInterval(final double alpha, final ConfidenceInterval side)
116     {
117         Throw.whenNull(side, "type of confidence level cannot be null");
118         Throw.when(alpha < 0 || alpha > 1, IllegalArgumentException.class,
119                 "confidenceLevel should be between 0 and 1 (inclusive)");
120         synchronized (this.semaphore)
121         {
122             double sampleMean = getSampleMean();
123             if (Double.isNaN(sampleMean) || Double.valueOf(this.getSampleStDev()).isNaN())
124             {
125                 return null; // TODO throw something
126             }
127             double level = 1 - alpha;
128             if (side.equals(ConfidenceInterval.BOTH_SIDE_CONFIDENCE))
129             {
130                 level = 1 - alpha / 2.0;
131             }
132             double z = DistNormalTable.getInverseCumulativeProbability(0.0, 1.0, level);
133             double confidence = z * Math.sqrt(this.getSampleVariance() / this.n);
134             double[] result = {sampleMean - confidence, sampleMean + confidence};
135             if (side.equals(ConfidenceInterval.LEFT_SIDE_CONFIDENCE))
136             {
137                 result[1] = sampleMean;
138             }
139             if (side.equals(ConfidenceInterval.RIGHT_SIDE_CONFIDENCE))
140             {
141                 result[0] = sampleMean;
142             }
143             result[0] = Math.max(result[0], this.min);
144             result[1] = Math.min(result[1], this.max);
145             return result;
146         }
147     }
148 
149     /** {@inheritDoc} */
150     @Override
151     public final String getDescription()
152     {
153         return this.description;
154     }
155 
156     /** {@inheritDoc} */
157     @Override
158     public final double getMax()
159     {
160         return this.max;
161     }
162 
163     /** {@inheritDoc} */
164     @Override
165     public final double getMin()
166     {
167         return this.min;
168     }
169 
170     /** {@inheritDoc} */
171     @Override
172     public final long getN()
173     {
174         return this.n;
175     }
176 
177     /** {@inheritDoc} */
178     @Override
179     public final double getSampleStDev()
180     {
181         synchronized (this.semaphore)
182         {
183             if (this.n > 1)
184             {
185                 return Math.sqrt(getSampleVariance());
186             }
187             return Double.NaN;
188         }
189     }
190 
191     /** {@inheritDoc} */
192     @Override
193     public final double getPopulationStDev()
194     {
195         synchronized (this.semaphore)
196         {
197             return Math.sqrt(getPopulationVariance());
198         }
199     }
200 
201     /** {@inheritDoc} */
202     @Override
203     public final double getSum()
204     {
205         return this.sum;
206     }
207 
208     /** {@inheritDoc} */
209     @Override
210     public final double getSampleVariance()
211     {
212         synchronized (this.semaphore)
213         {
214             if (this.n > 1)
215             {
216                 return this.m2 / (this.n - 1);
217             }
218             return Double.NaN;
219         }
220     }
221 
222     /** {@inheritDoc} */
223     @Override
224     public final double getPopulationVariance()
225     {
226         synchronized (this.semaphore)
227         {
228             if (this.n > 0)
229             {
230                 return this.m2 / this.n;
231             }
232             return Double.NaN;
233         }
234     }
235 
236     /** {@inheritDoc} */
237     @Override
238     public final double getSampleSkewness()
239     {
240         if (this.n > 2)
241         {
242             return getPopulationSkewness() * Math.sqrt(this.n * (this.n - 1)) / (this.n - 2);
243         }
244         return Double.NaN;
245     }
246 
247     /** {@inheritDoc} */
248     @Override
249     public final double getPopulationSkewness()
250     {
251         if (this.n > 1)
252         {
253             return (this.m3 / this.n) / Math.pow(getPopulationVariance(), 1.5);
254         }
255         return Double.NaN;
256     }
257 
258     /** {@inheritDoc} */
259     @Override
260     public final double getSampleKurtosis()
261     {
262         if (this.n > 3)
263         {
264             double sVar = getSampleVariance();
265             return this.m4 / (this.n - 1) / sVar / sVar;
266         }
267         return Double.NaN;
268     }
269 
270     /** {@inheritDoc} */
271     @Override
272     public final double getPopulationKurtosis()
273     {
274         if (this.n > 2)
275         {
276             return (this.m4 / this.n) / (this.m2 / this.n) / (this.m2 / this.n);
277         }
278         return Double.NaN;
279     }
280 
281     /** {@inheritDoc} */
282     @Override
283     public final double getSampleExcessKurtosis()
284     {
285         if (this.n > 3)
286         {
287             double g2 = getPopulationExcessKurtosis();
288             return (1.0 * (this.n - 1) / (this.n - 2) / (this.n - 3)) * ((this.n + 1) * g2 + 6.0);
289         }
290         return Double.NaN;
291     }
292 
293     /** {@inheritDoc} */
294     @Override
295     public final double getPopulationExcessKurtosis()
296     {
297         if (this.n > 2)
298         {
299             // convert kurtosis to excess kurtosis, shift by -3
300             return getPopulationKurtosis() - 3.0;
301         }
302         return Double.NaN;
303     }
304 
305     /** {@inheritDoc} */
306     @Override
307     public void initialize()
308     {
309         synchronized (this.semaphore)
310         {
311             this.min = Double.NaN;
312             this.max = Double.NaN;
313             this.n = 0;
314             this.sum = 0.0;
315             this.m1 = 0.0;
316             this.m2 = 0.0;
317             this.m3 = 0;
318             this.m4 = 0;
319             this.quantileAccumulator.initialize();
320         }
321     }
322 
323     /**
324      * Ingest an array of values.
325      * @param values double...; the values to ingest
326      */
327     public void ingest(final double... values)
328     {
329         for (double value : values)
330         {
331             ingest(value);
332         }
333     }
334 
335     /** {@inheritDoc} */
336     @Override
337     public double ingest(final double value)
338     {
339         Throw.when(Double.isNaN(value), IllegalArgumentException.class, "value may not be NaN");
340         synchronized (this.semaphore)
341         {
342             if (this.n == 0)
343             {
344                 this.min = Double.MAX_VALUE;
345                 this.max = -Double.MAX_VALUE;
346             }
347             this.n++;
348             double delta = value - this.m1;
349             double oldm2 = this.m2;
350             double oldm3 = this.m3;
351             // Eq 4 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
352             // Eq 1.1 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
353             this.m1 += delta / this.n;
354             // Eq 44 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
355             // Eq 1.2 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
356             this.m2 += delta * (value - this.m1);
357             // Eq 2.13 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
358             this.m3 += -3 * oldm2 * delta / this.n + (this.n - 1) * (this.n - 2) * delta * delta * delta / this.n / this.n;
359             // Eq 2.16 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
360             this.m4 += -4 * oldm3 * delta / this.n + 6 * oldm2 * delta * delta / this.n / this.n + (this.n - 1)
361                     * (this.n * this.n - 3 * this.n + 3) * delta * delta * delta * delta / this.n / this.n / this.n;
362             this.sum += value;
363             if (value < this.min)
364             {
365                 this.min = value;
366             }
367             if (value > this.max)
368             {
369                 this.max = value;
370             }
371             this.quantileAccumulator.ingest(value);
372         }
373         return value;
374     }
375 
376     /** {@inheritDoc} */
377     @Override
378     @SuppressWarnings("checkstyle:designforextension")
379     public String toString()
380     {
381         return "Tally [sum=" + this.sum + ", m1=" + this.m1 + ", m2=" + this.m2 + ", m3=" + this.m3 + ", m4=" + this.m4
382                 + ", min=" + this.min + ", max=" + this.max + ", n=" + this.n + ", description=" + this.description
383                 + ", quantileAccumulator=" + this.quantileAccumulator + "]";
384     }
385 
386 }