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 final double[] getConfidenceInterval(final double alpha)
102     {
103         return this.getConfidenceInterval(alpha, ConfidenceInterval.BOTH_SIDE_CONFIDENCE);
104     }
105 
106     /** {@inheritDoc} */
107     @Override
108     public final double[] getConfidenceInterval(final double alpha, final ConfidenceInterval side)
109     {
110         Throw.whenNull(side, "type of confidence level cannot be null");
111         Throw.when(alpha < 0 || alpha > 1, IllegalArgumentException.class,
112                 "confidenceLevel should be between 0 and 1 (inclusive)");
113         synchronized (this.semaphore)
114         {
115             double sampleMean = getSampleMean();
116             if (Double.isNaN(sampleMean) || Double.valueOf(this.getSampleStDev()).isNaN())
117             {
118                 return null; // TODO throw something
119             }
120             double level = 1 - alpha;
121             if (side.equals(ConfidenceInterval.BOTH_SIDE_CONFIDENCE))
122             {
123                 level = 1 - alpha / 2.0;
124             }
125             double z = DistNormalTable.getInverseCumulativeProbability(0.0, 1.0, level);
126             double confidence = z * Math.sqrt(this.getSampleVariance() / this.n);
127             double[] result = {sampleMean - confidence, sampleMean + confidence};
128             if (side.equals(ConfidenceInterval.LEFT_SIDE_CONFIDENCE))
129             {
130                 result[1] = sampleMean;
131             }
132             if (side.equals(ConfidenceInterval.RIGHT_SIDE_CONFIDENCE))
133             {
134                 result[0] = sampleMean;
135             }
136             result[0] = Math.max(result[0], this.min);
137             result[1] = Math.min(result[1], this.max);
138             return result;
139         }
140     }
141 
142     /** {@inheritDoc} */
143     @Override
144     public final String getDescription()
145     {
146         return this.description;
147     }
148 
149     /** {@inheritDoc} */
150     @Override
151     public final double getMax()
152     {
153         return this.max;
154     }
155 
156     /** {@inheritDoc} */
157     @Override
158     public final double getMin()
159     {
160         return this.min;
161     }
162 
163     /** {@inheritDoc} */
164     @Override
165     public final long getN()
166     {
167         return this.n;
168     }
169 
170     /** {@inheritDoc} */
171     @Override
172     public final double getSampleStDev()
173     {
174         synchronized (this.semaphore)
175         {
176             if (this.n > 1)
177             {
178                 return Math.sqrt(getSampleVariance());
179             }
180             return Double.NaN;
181         }
182     }
183 
184     /** {@inheritDoc} */
185     @Override
186     public final double getPopulationStDev()
187     {
188         synchronized (this.semaphore)
189         {
190             return Math.sqrt(getPopulationVariance());
191         }
192     }
193 
194     /** {@inheritDoc} */
195     @Override
196     public final double getSum()
197     {
198         return this.sum;
199     }
200 
201     /** {@inheritDoc} */
202     @Override
203     public final double getSampleVariance()
204     {
205         synchronized (this.semaphore)
206         {
207             if (this.n > 1)
208             {
209                 return this.m2 / (this.n - 1);
210             }
211             return Double.NaN;
212         }
213     }
214 
215     /** {@inheritDoc} */
216     @Override
217     public final double getPopulationVariance()
218     {
219         synchronized (this.semaphore)
220         {
221             if (this.n > 0)
222             {
223                 return this.m2 / this.n;
224             }
225             return Double.NaN;
226         }
227     }
228 
229     /** {@inheritDoc} */
230     @Override
231     public final double getSampleSkewness()
232     {
233         if (this.n > 2)
234         {
235             return getPopulationSkewness() * Math.sqrt(this.n * (this.n - 1)) / (this.n - 2);
236         }
237         return Double.NaN;
238     }
239 
240     /** {@inheritDoc} */
241     @Override
242     public final double getPopulationSkewness()
243     {
244         if (this.n > 1)
245         {
246             return (this.m3 / this.n) / Math.pow(getPopulationVariance(), 1.5);
247         }
248         return Double.NaN;
249     }
250 
251     /** {@inheritDoc} */
252     @Override
253     public final double getSampleKurtosis()
254     {
255         if (this.n > 3)
256         {
257             double sVar = getSampleVariance();
258             return this.m4 / (this.n - 1) / sVar / sVar;
259         }
260         return Double.NaN;
261     }
262 
263     /** {@inheritDoc} */
264     @Override
265     public final double getPopulationKurtosis()
266     {
267         if (this.n > 2)
268         {
269             return (this.m4 / this.n) / (this.m2 / this.n) / (this.m2 / this.n);
270         }
271         return Double.NaN;
272     }
273 
274     /** {@inheritDoc} */
275     @Override
276     public final double getSampleExcessKurtosis()
277     {
278         if (this.n > 3)
279         {
280             double g2 = getPopulationExcessKurtosis();
281             return (1.0 * (this.n - 1) / (this.n - 2) / (this.n - 3)) * ((this.n + 1) * g2 + 6.0);
282         }
283         return Double.NaN;
284     }
285 
286     /** {@inheritDoc} */
287     @Override
288     public final double getPopulationExcessKurtosis()
289     {
290         if (this.n > 2)
291         {
292             // convert kurtosis to excess kurtosis, shift by -3
293             return getPopulationKurtosis() - 3.0;
294         }
295         return Double.NaN;
296     }
297 
298     /** {@inheritDoc} */
299     @Override
300     public void initialize()
301     {
302         synchronized (this.semaphore)
303         {
304             this.min = Double.NaN;
305             this.max = Double.NaN;
306             this.n = 0;
307             this.sum = 0.0;
308             this.m1 = 0.0;
309             this.m2 = 0.0;
310             this.m3 = 0;
311             this.m4 = 0;
312             this.quantileAccumulator.initialize();
313         }
314     }
315 
316     /**
317      * Ingest an array of values.
318      * @param values double...; the values to ingest
319      */
320     public void ingest(final double... values)
321     {
322         for (double value : values)
323         {
324             ingest(value);
325         }
326     }
327 
328     /** {@inheritDoc} */
329     @Override
330     public double ingest(final double value)
331     {
332         Throw.when(Double.isNaN(value), IllegalArgumentException.class, "value may not be NaN");
333         synchronized (this.semaphore)
334         {
335             if (this.n == 0)
336             {
337                 this.min = Double.MAX_VALUE;
338                 this.max = -Double.MAX_VALUE;
339             }
340             this.n++;
341             double delta = value - this.m1;
342             double oldm2 = this.m2;
343             double oldm3 = this.m3;
344             // Eq 4 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
345             // Eq 1.1 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
346             this.m1 += delta / this.n;
347             // Eq 44 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
348             // Eq 1.2 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
349             this.m2 += delta * (value - this.m1);
350             // Eq 2.13 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
351             this.m3 += -3 * oldm2 * delta / this.n + (this.n - 1) * (this.n - 2) * delta * delta * delta / this.n / this.n;
352             // Eq 2.16 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
353             this.m4 += -4 * oldm3 * delta / this.n + 6 * oldm2 * delta * delta / this.n / this.n + (this.n - 1)
354                     * (this.n * this.n - 3 * this.n + 3) * delta * delta * delta * delta / this.n / this.n / this.n;
355             this.sum += value;
356             if (value < this.min)
357             {
358                 this.min = value;
359             }
360             if (value > this.max)
361             {
362                 this.max = value;
363             }
364             this.quantileAccumulator.ingest(value);
365         }
366         return value;
367     }
368 
369     /** {@inheritDoc} */
370     @Override
371     @SuppressWarnings("checkstyle:designforextension")
372     public String toString()
373     {
374         return "Tally [sum=" + this.sum + ", m1=" + this.m1 + ", m2=" + this.m2 + ", m3=" + this.m3 + ", m4=" + this.m4
375                 + ", min=" + this.min + ", max=" + this.max + ", n=" + this.n + ", description=" + this.description
376                 + ", quantileAccumulator=" + this.quantileAccumulator + "]";
377     }
378 
379 }