AbstractRealDistribution.java
/*
* 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.math4.legacy.distribution;
import org.apache.commons.statistics.distribution.ContinuousDistribution;
import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
import org.apache.commons.math4.legacy.exception.OutOfRangeException;
import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
import org.apache.commons.math4.core.jdkmath.JdkMath;
/**
* 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.
*
* <p>
* This base class provides a default factory method for creating
* a {@link org.apache.commons.statistics.distribution.ContinuousDistribution.Sampler
* sampler instance} that uses the
* <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
* inversion method</a> for generating random samples that follow the
* distribution.
* </p>
*
* @since 3.0
*/
public abstract class AbstractRealDistribution
implements ContinuousDistribution {
/** Default absolute accuracy for inverse cumulative computation. */
public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* 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}.
*
* The default implementation uses the identity
* {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
*
* @since 3.1
*/
@Override
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}
*
* The default implementation returns
* <ul>
* <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
* <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
* </ul>
*/
@Override
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 = getMean();
final double sig = JdkMath.sqrt(getVariance());
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 * JdkMath.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 * JdkMath.sqrt(p / (1. - p));
} else {
upperBound = 1.0;
while (cumulativeProbability(upperBound) < p) {
upperBound *= 2.0;
}
}
}
final UnivariateFunction toSolve = new UnivariateFunction() {
/** {@inheritDoc} */
@Override
public double value(final double x) {
return cumulativeProbability(x) - p;
}
};
return UnivariateSolverUtils.solve(toSolve,
lowerBound,
upperBound,
getSolverAbsoluteAccuracy());
}
/**
* 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 SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
}
/**
* {@inheritDoc}
* <p>
* The default implementation simply computes the logarithm of {@code density(x)}.
*/
@Override
public double logDensity(double x) {
return JdkMath.log(density(x));
}
/**
* Utility function for allocating an array and filling it with {@code n}
* samples generated by the given {@code sampler}.
*
* @param n Number of samples.
* @param sampler Sampler.
* @return an array of size {@code n}.
*/
public static double[] sample(int n,
ContinuousDistribution.Sampler sampler) {
final double[] samples = new double[n];
for (int i = 0; i < n; i++) {
samples[i] = sampler.sample();
}
return samples;
}
/**{@inheritDoc} */
@Override
public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
// Inversion method distribution sampler.
return InverseTransformContinuousSampler.of(rng, this::inverseCumulativeProbability)::sample;
}
}