    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     */
    017    package org.apache.commons.math3.distribution;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.analysis.UnivariateFunction;
    022    import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
    023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    024    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    025    import org.apache.commons.math3.exception.OutOfRangeException;
    026    import org.apache.commons.math3.exception.util.LocalizedFormats;
    027    import org.apache.commons.math3.random.RandomGenerator;
    028    import org.apache.commons.math3.random.RandomDataImpl;
    029    import org.apache.commons.math3.util.FastMath;
    030
    031    /**
    032     * Base class for probability distributions on the reals.
    033     * Default implementations are provided for some of the methods
    034     * that do not vary from distribution to distribution.
    035     *
    036     * @version $Id: AbstractRealDistribution.java 1422195 2012-12-15 06:45:18Z psteitz $
    037     * @since 3.0
    038     */
    039    public abstract class AbstractRealDistribution
    040    implements RealDistribution, Serializable {
    041        /** Default accuracy. */
    042        public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    043        /** Serializable version identifier */
    044        private static final long serialVersionUID = -38038050983108802L;
    045         /**
    046          * RandomData instance used to generate samples from the distribution.
    047          * @deprecated As of 3.1, to be removed in 4.0. Please use the
    048          * {@link #random} instance variable instead.
    049          */
    050        @Deprecated
    051        protected RandomDataImpl randomData = new RandomDataImpl();
    052
    053        /**
    054         * RNG instance used to generate samples from the distribution.
    055         * @since 3.1
    056         */
    057        protected final RandomGenerator random;
    058
    059        /** Solver absolute accuracy for inverse cumulative computation */
    060        private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY;
    061
    062        /**
    063         * @deprecated As of 3.1, to be removed in 4.0. Please use
    064         * {@link #AbstractRealDistribution(RandomGenerator)} instead.
    065         */
    066        @Deprecated
    067        protected AbstractRealDistribution() {
    068            // Legacy users are only allowed to access the deprecated "randomData".
    069            // New users are forbidden to use this constructor.
    070            random = null;
    071        }
    072        /**
    073         * @param rng Random number generator.
    074         * @since 3.1
    075         */
    076        protected AbstractRealDistribution(RandomGenerator rng) {
    077            random = rng;
    078        }
    079
    080        /**
    081         * {@inheritDoc}
    082         *
    083         * The default implementation uses the identity
    084         * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
    085         *
    086         * @deprecated As of 3.1 (to be removed in 4.0). Please use
    087         * {@link #probability(double,double)} instead.
    088         */
    089        @Deprecated
    090        public double cumulativeProbability(double x0, double x1) throws NumberIsTooLargeException {
    091            return probability(x0, x1);
    092        }
    093
    094        /**
    095         * For a random variable {@code X} whose values are distributed according
    096         * to this distribution, this method returns {@code P(x0 < X <= x1)}.
    097         *
    098         * @param x0 Lower bound (excluded).
    099         * @param x1 Upper bound (included).
    100         * @return the probability that a random variable with this distribution
    101         * takes a value between {@code x0} and {@code x1}, excluding the lower
    102         * and including the upper endpoint.
    103         * @throws NumberIsTooLargeException if {@code x0 > x1}.
    104         *
    105         * The default implementation uses the identity
    106         * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
    107         *
    108         * @since 3.1
    109         */
    110        public double probability(double x0,
    111                                  double x1) {
    112            if (x0 > x1) {
    113                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
    114                                                    x0, x1, true);
    115            }
    116            return cumulativeProbability(x1) - cumulativeProbability(x0);
    117        }
    118
    119        /**
    120         * {@inheritDoc}
    121         *
    122         * The default implementation returns
    123         * <ul>
    124         * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
    125         * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
    126         * </ul>
    127         */
    128        public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
    129            /*
    130             * IMPLEMENTATION NOTES
    131             * --------------------
    132             * Where applicable, use is made of the one-sided Chebyshev inequality
    133             * to bracket the root. This inequality states that
    134             * P(X - mu >= k * sig) <= 1 / (1 + k^2),
    135             * mu: mean, sig: standard deviation. Equivalently
    136             * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
    137             * F(mu + k * sig) >= k^2 / (1 + k^2).
    138             *
    139             * For k = sqrt(p / (1 - p)), we find
    140             * F(mu + k * sig) >= p,
    141             * and (mu + k * sig) is an upper-bound for the root.
    142             *
    143             * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
    144             * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
    145             * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
    146             * P(X <= mu - k * sig) <= 1 / (1 + k^2),
    147             * F(mu - k * sig) <= 1 / (1 + k^2).
    148             *
    149             * For k = sqrt((1 - p) / p), we find
    150             * F(mu - k * sig) <= p,
    151             * and (mu - k * sig) is a lower-bound for the root.
    152             *
    153             * In cases where the Chebyshev inequality does not apply, geometric
    154             * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
    155             * the root.
    156             */
    157            if (p < 0.0 || p > 1.0) {
    158                throw new OutOfRangeException(p, 0, 1);
    159            }
    160
    161            double lowerBound = getSupportLowerBound();
    162            if (p == 0.0) {
    163                return lowerBound;
    164            }
    165
    166            double upperBound = getSupportUpperBound();
    167            if (p == 1.0) {
    168                return upperBound;
    169            }
    170
    171            final double mu = getNumericalMean();
    172            final double sig = FastMath.sqrt(getNumericalVariance());
    173            final boolean chebyshevApplies;
    174            chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
    175                                 Double.isInfinite(sig) || Double.isNaN(sig));
    176
    177            if (lowerBound == Double.NEGATIVE_INFINITY) {
    178                if (chebyshevApplies) {
    179                    lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
    180                } else {
    181                    lowerBound = -1.0;
    182                    while (cumulativeProbability(lowerBound) >= p) {
    183                        lowerBound *= 2.0;
    184                    }
    185                }
    186            }
    187
    188            if (upperBound == Double.POSITIVE_INFINITY) {
    189                if (chebyshevApplies) {
    190                    upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
    191                } else {
    192                    upperBound = 1.0;
    193                    while (cumulativeProbability(upperBound) < p) {
    194                        upperBound *= 2.0;
    195                    }
    196                }
    197            }
    198
    199            final UnivariateFunction toSolve = new UnivariateFunction() {
    200
    201                public double value(final double x) {
    202                    return cumulativeProbability(x) - p;
    203                }
    204            };
    205
    206            double x = UnivariateSolverUtils.solve(toSolve,
    207                                                       lowerBound,
    208                                                       upperBound,
    209                                                       getSolverAbsoluteAccuracy());
    210
    211            if (!isSupportConnected()) {
    212                /* Test for plateau. */
    213                final double dx = getSolverAbsoluteAccuracy();
    214                if (x - dx >= getSupportLowerBound()) {
    215                    double px = cumulativeProbability(x);
    216                    if (cumulativeProbability(x - dx) == px) {
    217                        upperBound = x;
    218                        while (upperBound - lowerBound > dx) {
    219                            final double midPoint = 0.5 * (lowerBound + upperBound);
    220                            if (cumulativeProbability(midPoint) < px) {
    221                                lowerBound = midPoint;
    222                            } else {
    223                                upperBound = midPoint;
    224                            }
    225                        }
    226                        return upperBound;
    227                    }
    228                }
    229            }
    230            return x;
    231        }
    232
    233        /**
    234         * Returns the solver absolute accuracy for inverse cumulative computation.
    235         * You can override this method in order to use a Brent solver with an
    236         * absolute accuracy different from the default.
    237         *
    238         * @return the maximum absolute error in inverse cumulative probability estimates
    239         */
    240        protected double getSolverAbsoluteAccuracy() {
    241            return solverAbsoluteAccuracy;
    242        }
    243
    244        /** {@inheritDoc} */
    245        public void reseedRandomGenerator(long seed) {
    246            random.setSeed(seed);
    247            randomData.reSeed(seed);
    248        }
    249
    250        /**
    251         * {@inheritDoc}
    252         *
    253         * The default implementation uses the
    254         * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
    255         * inversion method.
    256         * </a>
    257         */
    258        public double sample() {
    259            return inverseCumulativeProbability(random.nextDouble());
    260        }
    261
    262        /**
    263         * {@inheritDoc}
    264         *
    265         * The default implementation generates the sample by calling
    266         * {@link #sample()} in a loop.
    267         */
    268        public double[] sample(int sampleSize) {
    269            if (sampleSize <= 0) {
    270                throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES,
    271                        sampleSize);
    272            }
    273            double[] out = new double[sampleSize];
    274            for (int i = 0; i < sampleSize; i++) {
    275                out[i] = sample();
    276            }
    277            return out;
    278        }
    279
    280        /**
    281         * {@inheritDoc}
    282         *
    283         * @return zero.
    284         * @since 3.1
    285         */
    286        public double probability(double x) {
    287            return 0d;
    288        }
    289    }
    290