Tally.java
package org.djutils.stats.summarizers;
import org.djutils.exceptions.Throw;
import org.djutils.stats.ConfidenceInterval;
import org.djutils.stats.DistNormalTable;
import org.djutils.stats.summarizers.quantileaccumulator.NoStorageAccumulator;
import org.djutils.stats.summarizers.quantileaccumulator.QuantileAccumulator;
/**
* The Tally class ingests a series of values and provides mean, standard deviation, etc. of the ingested values.
* <p>
* Copyright (c) 2002-2021 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
* for project information <a href="https://simulation.tudelft.nl/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
* project is distributed under a three-clause BSD-style license, which can be found at
* <a href="https://simulation.tudelft.nl/dsol/3.0/license.html" target="_blank">
* https://simulation.tudelft.nl/dsol/3.0/license.html</a>. <br>
* @author <a href="https://www.tudelft.nl/averbraeck" target="_blank"> Alexander Verbraeck</a>
* @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
* @author <a href="https://www.tudelft.nl/staff/p.knoppers/">Peter Knoppers</a>
*/
public class Tally implements TallyInterface
{
/** */
private static final long serialVersionUID = 20200228L;
/** The sum of this tally. */
private double sum = 0;
/** The mean of this tally. */
private double m1 = 0;
/** The summation for the second moment (variance). */
private double m2 = 0;
/** The summation for the third moment (skewness). */
private double m3 = 0;
/** The summation for the fourth moment (kurtosis). */
private double m4 = 0;
/** The minimum observed value of this tally. */
private double min = Double.NaN;
/** The maximum observed value of this tally. */
private double max = Double.NaN;
/** The number of measurements of this tally. */
private long n = 0;
/** The description of this tally. */
private final String description;
/** The quantile accumulator. */
private final QuantileAccumulator quantileAccumulator;
/** the synchronized lock. */
@SuppressWarnings("checkstyle:visibilitymodifier")
protected Object semaphore = new Object();
/**
* Constructs a new Tally.
* @param description String; the description of this tally
* @param quantileAccumulator QuantileAccumulator; the input series accumulator that can approximate or compute quantiles.
*/
public Tally(final String description, final QuantileAccumulator quantileAccumulator)
{
this.description = description;
this.quantileAccumulator = quantileAccumulator;
initialize();
}
/**
* Convenience constructor that uses a NoStorageAccumulator to estimate quantiles.
* @param description String; the description of this tally
*/
public Tally(final String description)
{
this(description, new NoStorageAccumulator());
}
/** {@inheritDoc} */
@Override
public final double getSampleMean()
{
if (this.n > 0)
{
return this.m1;
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getQuantile(final double probability)
{
return this.quantileAccumulator.getQuantile(this, probability);
}
/** {@inheritDoc} */
@Override
public final double[] getConfidenceInterval(final double alpha)
{
return this.getConfidenceInterval(alpha, ConfidenceInterval.BOTH_SIDE_CONFIDENCE);
}
/** {@inheritDoc} */
@Override
public final double[] getConfidenceInterval(final double alpha, final ConfidenceInterval side)
{
Throw.whenNull(side, "type of confidence level cannot be null");
Throw.when(alpha < 0 || alpha > 1, IllegalArgumentException.class,
"confidenceLevel should be between 0 and 1 (inclusive)");
synchronized (this.semaphore)
{
double sampleMean = getSampleMean();
if (Double.isNaN(sampleMean) || Double.valueOf(this.getSampleStDev()).isNaN())
{
return null; // TODO throw something
}
double level = 1 - alpha;
if (side.equals(ConfidenceInterval.BOTH_SIDE_CONFIDENCE))
{
level = 1 - alpha / 2.0;
}
double z = DistNormalTable.getInverseCumulativeProbability(0.0, 1.0, level);
double confidence = z * Math.sqrt(this.getSampleVariance() / this.n);
double[] result = {sampleMean - confidence, sampleMean + confidence};
if (side.equals(ConfidenceInterval.LEFT_SIDE_CONFIDENCE))
{
result[1] = sampleMean;
}
if (side.equals(ConfidenceInterval.RIGHT_SIDE_CONFIDENCE))
{
result[0] = sampleMean;
}
result[0] = Math.max(result[0], this.min);
result[1] = Math.min(result[1], this.max);
return result;
}
}
/** {@inheritDoc} */
@Override
public final String getDescription()
{
return this.description;
}
/** {@inheritDoc} */
@Override
public final double getMax()
{
return this.max;
}
/** {@inheritDoc} */
@Override
public final double getMin()
{
return this.min;
}
/** {@inheritDoc} */
@Override
public final long getN()
{
return this.n;
}
/** {@inheritDoc} */
@Override
public final double getSampleStDev()
{
synchronized (this.semaphore)
{
if (this.n > 1)
{
return Math.sqrt(getSampleVariance());
}
return Double.NaN;
}
}
/** {@inheritDoc} */
@Override
public final double getPopulationStDev()
{
synchronized (this.semaphore)
{
return Math.sqrt(getPopulationVariance());
}
}
/** {@inheritDoc} */
@Override
public final double getSum()
{
return this.sum;
}
/** {@inheritDoc} */
@Override
public final double getSampleVariance()
{
synchronized (this.semaphore)
{
if (this.n > 1)
{
return this.m2 / (this.n - 1);
}
return Double.NaN;
}
}
/** {@inheritDoc} */
@Override
public final double getPopulationVariance()
{
synchronized (this.semaphore)
{
if (this.n > 0)
{
return this.m2 / this.n;
}
return Double.NaN;
}
}
/** {@inheritDoc} */
@Override
public final double getSampleSkewness()
{
if (this.n > 2)
{
return getPopulationSkewness() * Math.sqrt(this.n * (this.n - 1)) / (this.n - 2);
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getPopulationSkewness()
{
if (this.n > 1)
{
return (this.m3 / this.n) / Math.pow(getPopulationVariance(), 1.5);
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getSampleKurtosis()
{
if (this.n > 3)
{
double sVar = getSampleVariance();
return this.m4 / (this.n - 1) / sVar / sVar;
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getPopulationKurtosis()
{
if (this.n > 2)
{
return (this.m4 / this.n) / (this.m2 / this.n) / (this.m2 / this.n);
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getSampleExcessKurtosis()
{
if (this.n > 3)
{
double g2 = getPopulationExcessKurtosis();
return (1.0 * (this.n - 1) / (this.n - 2) / (this.n - 3)) * ((this.n + 1) * g2 + 6.0);
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public final double getPopulationExcessKurtosis()
{
if (this.n > 2)
{
// convert kurtosis to excess kurtosis, shift by -3
return getPopulationKurtosis() - 3.0;
}
return Double.NaN;
}
/** {@inheritDoc} */
@Override
public void initialize()
{
synchronized (this.semaphore)
{
this.min = Double.NaN;
this.max = Double.NaN;
this.n = 0;
this.sum = 0.0;
this.m1 = 0.0;
this.m2 = 0.0;
this.m3 = 0;
this.m4 = 0;
this.quantileAccumulator.initialize();
}
}
/**
* Ingest an array of values.
* @param values double...; the values to ingest
*/
public void ingest(final double... values)
{
for (double value : values)
{
ingest(value);
}
}
/** {@inheritDoc} */
@Override
public double ingest(final double value)
{
Throw.when(Double.isNaN(value), IllegalArgumentException.class, "value may not be NaN");
synchronized (this.semaphore)
{
if (this.n == 0)
{
this.min = Double.MAX_VALUE;
this.max = -Double.MAX_VALUE;
}
this.n++;
double delta = value - this.m1;
double oldm2 = this.m2;
double oldm3 = this.m3;
// Eq 4 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
// Eq 1.1 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
this.m1 += delta / this.n;
// Eq 44 in https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
// Eq 1.2 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
this.m2 += delta * (value - this.m1);
// Eq 2.13 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
this.m3 += -3 * oldm2 * delta / this.n + (this.n - 1) * (this.n - 2) * delta * delta * delta / this.n / this.n;
// Eq 2.16 in https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
this.m4 += -4 * oldm3 * delta / this.n + 6 * oldm2 * delta * delta / this.n / this.n + (this.n - 1)
* (this.n * this.n - 3 * this.n + 3) * delta * delta * delta * delta / this.n / this.n / this.n;
this.sum += value;
if (value < this.min)
{
this.min = value;
}
if (value > this.max)
{
this.max = value;
}
this.quantileAccumulator.ingest(value);
}
return value;
}
/** {@inheritDoc} */
@Override
@SuppressWarnings("checkstyle:designforextension")
public String toString()
{
return "Tally [sum=" + this.sum + ", m1=" + this.m1 + ", m2=" + this.m2 + ", m3=" + this.m3 + ", m4=" + this.m4
+ ", min=" + this.min + ", max=" + this.max + ", n=" + this.n + ", description=" + this.description
+ ", quantileAccumulator=" + this.quantileAccumulator + "]";
}
}