DistNormalTrunc.java
package org.opentrafficsim.core.dsol;
import static nl.tudelft.simulation.jstats.distributions.DistNormal.CUMULATIVE_NORMAL_PROBABILITIES;
import java.util.Locale;
import nl.tudelft.simulation.jstats.distributions.DistContinuous;
import nl.tudelft.simulation.jstats.streams.MersenneTwister;
import nl.tudelft.simulation.jstats.streams.StreamInterface;
/**
* The Normal Truncated distribution. For more information on the normal distribution see
* <a href="http://mathworld.wolfram.com/NormalDistribution.html"> http://mathworld.wolfram.com/NormalDistribution.html </a>
* <p>
* This version of the normal distribution uses the numerically approached inverse cumulative distribution.
* <p>
* (c) copyright 2002-2018 <a href="http://www.simulation.tudelft.nl">Delft University of Technology </a>, the Netherlands. <br>
* See for project information <a href="http://www.simulation.tudelft.nl"> www.simulation.tudelft.nl </a> <br>
* License of use: <a href="http://www.gnu.org/copyleft/lesser.html">Lesser General Public License (LGPL) </a>, no warranty.
* @author <a href="mailto:a.verbraeck@tudelft.nl"> Alexander Verbraeck </a> <br>
* <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
*/
public class DistNormalTrunc extends DistContinuous
{
/** */
private static final long serialVersionUID = 1L;
/** mu refers to the mean of the normal distribution. */
private final double mu;
/** mu refers to the mean of the normal distribution. */
private final double sigma;
/** minimum x-value of the distribution. */
private final double min;
/** maximum x-value of the distribution. */
private final double max;
/** cumulative distribution value of the minimum. */
private final double cumulMin;
/** cumulative distribution value difference between max and min. */
private final double cumulDiff;
/** factor on probability density to normalize to 1. */
private final double probDensFactor;
/**
* constructs a normal distribution with mu=0 and sigma=1. Errors of various types, e.g., in the impact point of a bomb;
* quantities that are the sum of a large number of other quantities by the virtue of the central limit theorem.
* @param stream StreamInterface; the numberstream
* @param min double; minimum x-value of the distribution
* @param max double; maximum x-value of the distribution
*/
public DistNormalTrunc(final StreamInterface stream, final double min, final double max)
{
this(stream, 0.0, 1.0, min, max);
}
/**
* constructs a normal distribution with mu and sigma.
* @param stream StreamInterface; the numberstream
* @param mu double; the medium
* @param sigma double; the standard deviation
* @param min double; minimum x-value of the distribution
* @param max double; maximum x-value of the distribution
*/
public DistNormalTrunc(final StreamInterface stream, final double mu, final double sigma, final double min,
final double max)
{
super(stream);
if (max < min)
{
throw new IllegalArgumentException("Error Normal Truncated - max < min");
}
this.mu = mu;
this.sigma = sigma;
this.min = min;
this.max = max;
this.cumulMin = getCumulativeProbabilityNotTruncated(min);
this.cumulDiff = getCumulativeProbabilityNotTruncated(max) - this.cumulMin;
this.probDensFactor = 1.0 / this.cumulDiff;
}
/** {@inheritDoc} */
@Override
public double draw()
{
return getInverseCumulativeProbabilityNotTruncated(this.cumulMin + this.cumulDiff * this.stream.nextDouble());
}
/**
* returns the cumulative probability of the x-value.
* @param x double; the observation x
* @return double the cumulative probability
*/
public double getCumulativeProbability(final double x)
{
if (x < this.min)
{
return 0.0;
}
if (x > this.max)
{
return 1.0;
}
return (getCumulativeProbabilityNotTruncated(x) - this.cumulMin) * this.probDensFactor;
}
/**
* returns the cumulative probability of the x-value.
* @param x double; the observation x
* @return double the cumulative probability
*/
private double getCumulativeProbabilityNotTruncated(final double x)
{
double z = (x - this.mu) / this.sigma * 100;
double absZ = Math.abs(z);
int intZ = (int) absZ;
double f = 0.0;
if (intZ >= 1000)
{
intZ = 999;
f = 1.0;
}
else
{
f = absZ - intZ;
}
if (z >= 0)
{
return (1 - f) * CUMULATIVE_NORMAL_PROBABILITIES[intZ] + f * CUMULATIVE_NORMAL_PROBABILITIES[intZ + 1];
}
return 1 - ((1 - f) * CUMULATIVE_NORMAL_PROBABILITIES[intZ] + f * CUMULATIVE_NORMAL_PROBABILITIES[intZ + 1]);
}
/**
* returns the x-value of the given cumulativePropability.
* @param cumulativeProbability double; reflects cum prob
* @return double the inverse cumulative probability
*/
public double getInverseCumulativeProbability(final double cumulativeProbability)
{
if (cumulativeProbability < 0 || cumulativeProbability > 1)
{
throw new IllegalArgumentException("1<cumulativeProbability<0 ?");
}
/*
* For extreme cases we return the min and max directly. The method getInverseCumulativeProbabilityNotTruncated() can
* only return values from "mu - 10*sigma" to "mu + 10*sigma". If min or max is beyond these values, those values would
* result erroneously. For any cumulative probability that is slightly above 0.0 or slightly below 1.0, values in the
* range from "mu - 10*sigma" to "mu + 10*sigma" will always result.
*/
if (cumulativeProbability == 0.0)
{
return this.min;
}
if (cumulativeProbability == 1.0)
{
return this.max;
}
return getInverseCumulativeProbabilityNotTruncated(this.cumulMin + cumulativeProbability * this.cumulDiff);
}
/** {@inheritDoc} */
@Override
public double probDensity(final double x)
{
if (x < this.min || x > this.max)
{
return 0.0;
}
return this.probDensFactor / (Math.sqrt(2 * Math.PI * Math.pow(this.sigma, 2)))
* Math.exp(-1 * Math.pow(x - this.mu, 2) / (2 * Math.pow(this.sigma, 2)));
}
/** {@inheritDoc} */
@Override
public String toString()
{
return "NormalTrunc(" + this.mu + "," + this.sigma + "," + this.min + "," + this.max + ")";
}
/**
* Test.
* @param args String[]; args
*/
public static void main(final String[] args)
{
StreamInterface stream = new MersenneTwister();
double mu = 2.0;
double sigma = 3.0;
double min = -5.0;
double max = 4.0;
DistNormalTrunc dist = new DistNormalTrunc(stream, mu, sigma, min, max);
System.out.println("<< probability density >>");
double sum = 0.0;
double step = (max - min) / 96;
for (double x = min - 2 * step; x <= max + 2 * step; x += step)
{
double p = dist.probDensity(x);
System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", x, p));
sum += p * step;
}
System.out.println(String.format(Locale.GERMAN, "Approx. sum = %.8f", sum));
System.out.println("");
System.out.println("<< cumulative density >>");
for (double x = min - 2 * step; x <= max + 2 * step; x += step)
{
double c = dist.getCumulativeProbability(x);
System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", x, c));
}
System.out.println("");
System.out.println("<< inverse cumulative density >>");
for (double c = 0.0; c < 1.005; c += 0.01)
{
double x = dist.getInverseCumulativeProbability(Math.min(c, 1.0)); // want to include 1.0, also if 1.0000000000001
System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", c, x));
}
System.out.println("");
System.out.println("<< 10000 random numbers. >>");
for (int i = 1; i < 10000; i++)
{
System.out.println(String.format(Locale.GERMAN, "%,8f", dist.draw()));
}
}
/**
* returns the x-value of the given cumulativePropability.
* @param cumulativeProbability double; reflects cum prob
* @return double the inverse cumulative probability
*/
private double getInverseCumulativeProbabilityNotTruncated(final double cumulativeProbability)
{
if (cumulativeProbability < 0 || cumulativeProbability > 1)
{
throw new IllegalArgumentException("1<cumulativeProbability<0 ?");
}
boolean located = false;
double prob = cumulativeProbability;
if (cumulativeProbability < 0.5)
{
prob = 1 - cumulativeProbability;
}
int i = 0;
double f = 0.0;
while (!located)
{
if (CUMULATIVE_NORMAL_PROBABILITIES[i] < prob && CUMULATIVE_NORMAL_PROBABILITIES[i + 1] >= prob)
{
located = true;
if (CUMULATIVE_NORMAL_PROBABILITIES[i] < CUMULATIVE_NORMAL_PROBABILITIES[i + 1])
{
f = (prob - CUMULATIVE_NORMAL_PROBABILITIES[i])
/ (CUMULATIVE_NORMAL_PROBABILITIES[i + 1] - CUMULATIVE_NORMAL_PROBABILITIES[i]);
}
}
else
{
i++;
}
}
if (cumulativeProbability < 0.5)
{
return this.mu - ((f + i) / 100.0) * this.sigma;
}
return ((f + i) / 100.0) * this.sigma + this.mu;
}
}