diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java | 306 |
1 files changed, 306 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java b/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java new file mode 100644 index 0000000..b9e5bca --- /dev/null +++ b/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java @@ -0,0 +1,306 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.distribution; + +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils; +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.NumberIsTooLargeException; +import org.apache.commons.math3.exception.OutOfRangeException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.util.FastMath; + +import java.io.Serializable; + +/** + * Base class for probability distributions on the reals. Default implementations are provided for + * some of the methods that do not vary from distribution to distribution. + * + * @since 3.0 + */ +public abstract class AbstractRealDistribution implements RealDistribution, Serializable { + /** Default accuracy. */ + public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6; + + /** Serializable version identifier */ + private static final long serialVersionUID = -38038050983108802L; + + /** + * RandomData instance used to generate samples from the distribution. + * + * @deprecated As of 3.1, to be removed in 4.0. Please use the {@link #random} instance variable + * instead. + */ + @Deprecated + protected org.apache.commons.math3.random.RandomDataImpl randomData = + new org.apache.commons.math3.random.RandomDataImpl(); + + /** + * RNG instance used to generate samples from the distribution. + * + * @since 3.1 + */ + protected final RandomGenerator random; + + /** Solver absolute accuracy for inverse cumulative computation */ + private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY; + + /** + * @deprecated As of 3.1, to be removed in 4.0. Please use {@link + * #AbstractRealDistribution(RandomGenerator)} instead. + */ + @Deprecated + protected AbstractRealDistribution() { + // Legacy users are only allowed to access the deprecated "randomData". + // New users are forbidden to use this constructor. + random = null; + } + + /** + * @param rng Random number generator. + * @since 3.1 + */ + protected AbstractRealDistribution(RandomGenerator rng) { + random = rng; + } + + /** + * {@inheritDoc} + * + * <p>The default implementation uses the identity + * + * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)} + * + * @deprecated As of 3.1 (to be removed in 4.0). Please use {@link #probability(double,double)} + * instead. + */ + @Deprecated + public double cumulativeProbability(double x0, double x1) throws NumberIsTooLargeException { + return probability(x0, x1); + } + + /** + * For a random variable {@code X} whose values are distributed according to this distribution, + * this method returns {@code P(x0 < X <= x1)}. + * + * @param x0 Lower bound (excluded). + * @param x1 Upper bound (included). + * @return the probability that a random variable with this distribution takes a value between + * {@code x0} and {@code x1}, excluding the lower and including the upper endpoint. + * @throws NumberIsTooLargeException if {@code x0 > x1}. + * <p>The default implementation uses the identity {@code P(x0 < X <= x1) = P(X <= x1) - P(X + * <= x0)} + * @since 3.1 + */ + public double probability(double x0, double x1) { + if (x0 > x1) { + throw new NumberIsTooLargeException( + LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, x0, x1, true); + } + return cumulativeProbability(x1) - cumulativeProbability(x0); + } + + /** + * {@inheritDoc} + * + * <p>The default implementation returns + * + * <ul> + * <li>{@link #getSupportLowerBound()} for {@code p = 0}, + * <li>{@link #getSupportUpperBound()} for {@code p = 1}. + * </ul> + */ + public double inverseCumulativeProbability(final double p) throws OutOfRangeException { + /* + * IMPLEMENTATION NOTES + * -------------------- + * Where applicable, use is made of the one-sided Chebyshev inequality + * to bracket the root. This inequality states that + * P(X - mu >= k * sig) <= 1 / (1 + k^2), + * mu: mean, sig: standard deviation. Equivalently + * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2), + * F(mu + k * sig) >= k^2 / (1 + k^2). + * + * For k = sqrt(p / (1 - p)), we find + * F(mu + k * sig) >= p, + * and (mu + k * sig) is an upper-bound for the root. + * + * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and + * P(Y >= -mu + k * sig) <= 1 / (1 + k^2), + * P(-X >= -mu + k * sig) <= 1 / (1 + k^2), + * P(X <= mu - k * sig) <= 1 / (1 + k^2), + * F(mu - k * sig) <= 1 / (1 + k^2). + * + * For k = sqrt((1 - p) / p), we find + * F(mu - k * sig) <= p, + * and (mu - k * sig) is a lower-bound for the root. + * + * In cases where the Chebyshev inequality does not apply, geometric + * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket + * the root. + */ + if (p < 0.0 || p > 1.0) { + throw new OutOfRangeException(p, 0, 1); + } + + double lowerBound = getSupportLowerBound(); + if (p == 0.0) { + return lowerBound; + } + + double upperBound = getSupportUpperBound(); + if (p == 1.0) { + return upperBound; + } + + final double mu = getNumericalMean(); + final double sig = FastMath.sqrt(getNumericalVariance()); + final boolean chebyshevApplies; + chebyshevApplies = + !(Double.isInfinite(mu) + || Double.isNaN(mu) + || Double.isInfinite(sig) + || Double.isNaN(sig)); + + if (lowerBound == Double.NEGATIVE_INFINITY) { + if (chebyshevApplies) { + lowerBound = mu - sig * FastMath.sqrt((1. - p) / p); + } else { + lowerBound = -1.0; + while (cumulativeProbability(lowerBound) >= p) { + lowerBound *= 2.0; + } + } + } + + if (upperBound == Double.POSITIVE_INFINITY) { + if (chebyshevApplies) { + upperBound = mu + sig * FastMath.sqrt(p / (1. - p)); + } else { + upperBound = 1.0; + while (cumulativeProbability(upperBound) < p) { + upperBound *= 2.0; + } + } + } + + final UnivariateFunction toSolve = + new UnivariateFunction() { + /** {@inheritDoc} */ + public double value(final double x) { + return cumulativeProbability(x) - p; + } + }; + + double x = + UnivariateSolverUtils.solve( + toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy()); + + if (!isSupportConnected()) { + /* Test for plateau. */ + final double dx = getSolverAbsoluteAccuracy(); + if (x - dx >= getSupportLowerBound()) { + double px = cumulativeProbability(x); + if (cumulativeProbability(x - dx) == px) { + upperBound = x; + while (upperBound - lowerBound > dx) { + final double midPoint = 0.5 * (lowerBound + upperBound); + if (cumulativeProbability(midPoint) < px) { + lowerBound = midPoint; + } else { + upperBound = midPoint; + } + } + return upperBound; + } + } + } + return x; + } + + /** + * Returns the solver absolute accuracy for inverse cumulative computation. You can override + * this method in order to use a Brent solver with an absolute accuracy different from the + * default. + * + * @return the maximum absolute error in inverse cumulative probability estimates + */ + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** {@inheritDoc} */ + public void reseedRandomGenerator(long seed) { + random.setSeed(seed); + randomData.reSeed(seed); + } + + /** + * {@inheritDoc} + * + * <p>The default implementation uses the <a + * href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">inversion method. </a> + */ + public double sample() { + return inverseCumulativeProbability(random.nextDouble()); + } + + /** + * {@inheritDoc} + * + * <p>The default implementation generates the sample by calling {@link #sample()} in a loop. + */ + public double[] sample(int sampleSize) { + if (sampleSize <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, sampleSize); + } + double[] out = new double[sampleSize]; + for (int i = 0; i < sampleSize; i++) { + out[i] = sample(); + } + return out; + } + + /** + * {@inheritDoc} + * + * @return zero. + * @since 3.1 + */ + public double probability(double x) { + return 0d; + } + + /** + * Returns the natural logarithm of the probability density function (PDF) of this distribution + * evaluated at the specified point {@code x}. In general, the PDF is the derivative of the + * {@link #cumulativeProbability(double) CDF}. If the derivative does not exist at {@code x}, + * then an appropriate replacement should be returned, e.g. {@code Double.POSITIVE_INFINITY}, + * {@code Double.NaN}, or the limit inferior or limit superior of the difference quotient. Note + * that due to the floating point precision and under/overflow issues, this method will for some + * distributions be more precise and faster than computing the logarithm of {@link + * #density(double)}. The default implementation simply computes the logarithm of {@code + * density(x)}. + * + * @param x the point at which the PDF is evaluated + * @return the logarithm of the value of the probability density function at point {@code x} + */ + public double logDensity(double x) { + return FastMath.log(density(x)); + } +} |