001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.distribution;
018
019import org.apache.commons.statistics.distribution.ContinuousDistribution;
020import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
021import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
022import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
023import org.apache.commons.math4.legacy.exception.OutOfRangeException;
024import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
025import org.apache.commons.rng.UniformRandomProvider;
026import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
027import org.apache.commons.math4.core.jdkmath.JdkMath;
028
029/**
030 * Base class for probability distributions on the reals.
031 * Default implementations are provided for some of the methods
032 * that do not vary from distribution to distribution.
033 *
034 * <p>
035 * This base class provides a default factory method for creating
036 * a {@link org.apache.commons.statistics.distribution.ContinuousDistribution.Sampler
037 * sampler instance} that uses the
038 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
039 * inversion method</a> for generating random samples that follow the
040 * distribution.
041 * </p>
042 *
043 * @since 3.0
044 */
045public abstract class AbstractRealDistribution
046    implements ContinuousDistribution {
047    /** Default absolute accuracy for inverse cumulative computation. */
048    public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
049
050    /**
051     * For a random variable {@code X} whose values are distributed according
052     * to this distribution, this method returns {@code P(x0 < X <= x1)}.
053     *
054     * @param x0 Lower bound (excluded).
055     * @param x1 Upper bound (included).
056     * @return the probability that a random variable with this distribution
057     * takes a value between {@code x0} and {@code x1}, excluding the lower
058     * and including the upper endpoint.
059     * @throws NumberIsTooLargeException if {@code x0 > x1}.
060     *
061     * The default implementation uses the identity
062     * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
063     *
064     * @since 3.1
065     */
066    @Override
067    public double probability(double x0,
068                              double x1) {
069        if (x0 > x1) {
070            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
071                                                x0, x1, true);
072        }
073        return cumulativeProbability(x1) - cumulativeProbability(x0);
074    }
075
076    /**
077     * {@inheritDoc}
078     *
079     * The default implementation returns
080     * <ul>
081     * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
082     * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
083     * </ul>
084     */
085    @Override
086    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
087        /*
088         * IMPLEMENTATION NOTES
089         * --------------------
090         * Where applicable, use is made of the one-sided Chebyshev inequality
091         * to bracket the root. This inequality states that
092         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
093         * mu: mean, sig: standard deviation. Equivalently
094         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
095         * F(mu + k * sig) >= k^2 / (1 + k^2).
096         *
097         * For k = sqrt(p / (1 - p)), we find
098         * F(mu + k * sig) >= p,
099         * and (mu + k * sig) is an upper-bound for the root.
100         *
101         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
102         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
103         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
104         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
105         * F(mu - k * sig) <= 1 / (1 + k^2).
106         *
107         * For k = sqrt((1 - p) / p), we find
108         * F(mu - k * sig) <= p,
109         * and (mu - k * sig) is a lower-bound for the root.
110         *
111         * In cases where the Chebyshev inequality does not apply, geometric
112         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
113         * the root.
114         */
115        if (p < 0.0 || p > 1.0) {
116            throw new OutOfRangeException(p, 0, 1);
117        }
118
119        double lowerBound = getSupportLowerBound();
120        if (p == 0.0) {
121            return lowerBound;
122        }
123
124        double upperBound = getSupportUpperBound();
125        if (p == 1.0) {
126            return upperBound;
127        }
128
129        final double mu = getMean();
130        final double sig = JdkMath.sqrt(getVariance());
131        final boolean chebyshevApplies;
132        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
133                             Double.isInfinite(sig) || Double.isNaN(sig));
134
135        if (lowerBound == Double.NEGATIVE_INFINITY) {
136            if (chebyshevApplies) {
137                lowerBound = mu - sig * JdkMath.sqrt((1. - p) / p);
138            } else {
139                lowerBound = -1.0;
140                while (cumulativeProbability(lowerBound) >= p) {
141                    lowerBound *= 2.0;
142                }
143            }
144        }
145
146        if (upperBound == Double.POSITIVE_INFINITY) {
147            if (chebyshevApplies) {
148                upperBound = mu + sig * JdkMath.sqrt(p / (1. - p));
149            } else {
150                upperBound = 1.0;
151                while (cumulativeProbability(upperBound) < p) {
152                    upperBound *= 2.0;
153                }
154            }
155        }
156
157        final UnivariateFunction toSolve = new UnivariateFunction() {
158            /** {@inheritDoc} */
159            @Override
160            public double value(final double x) {
161                return cumulativeProbability(x) - p;
162            }
163        };
164
165        return UnivariateSolverUtils.solve(toSolve,
166                                           lowerBound,
167                                           upperBound,
168                                           getSolverAbsoluteAccuracy());
169    }
170
171    /**
172     * Returns the solver absolute accuracy for inverse cumulative computation.
173     * You can override this method in order to use a Brent solver with an
174     * absolute accuracy different from the default.
175     *
176     * @return the maximum absolute error in inverse cumulative probability estimates
177     */
178    protected double getSolverAbsoluteAccuracy() {
179        return SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
180    }
181
182    /**
183     * {@inheritDoc}
184     * <p>
185     * The default implementation simply computes the logarithm of {@code density(x)}.
186     */
187    @Override
188    public double logDensity(double x) {
189        return JdkMath.log(density(x));
190    }
191
192    /**
193     * Utility function for allocating an array and filling it with {@code n}
194     * samples generated by the given {@code sampler}.
195     *
196     * @param n Number of samples.
197     * @param sampler Sampler.
198     * @return an array of size {@code n}.
199     */
200    public static double[] sample(int n,
201                                  ContinuousDistribution.Sampler sampler) {
202        final double[] samples = new double[n];
203        for (int i = 0; i < n; i++) {
204            samples[i] = sampler.sample();
205        }
206        return samples;
207    }
208
209    /**{@inheritDoc} */
210    @Override
211    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
212        // Inversion method distribution sampler.
213        return InverseTransformContinuousSampler.of(rng, this::inverseCumulativeProbability)::sample;
214    }
215}