    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 org.apache.commons.math3.exception.NotStrictlyPositiveException;
    020    import org.apache.commons.math3.exception.OutOfRangeException;
    021    import org.apache.commons.math3.exception.util.LocalizedFormats;
    022    import org.apache.commons.math3.util.FastMath;
    023    import org.apache.commons.math3.util.ArithmeticUtils;
    024    import org.apache.commons.math3.util.ResizableDoubleArray;
    025    import org.apache.commons.math3.random.RandomGenerator;
    026    import org.apache.commons.math3.random.Well19937c;
    027
    028    /**
    029     * Implementation of the exponential distribution.
    030     *
    031     * @see <a href="http://en.wikipedia.org/wiki/Exponential_distribution">Exponential distribution (Wikipedia)</a>
    032     * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a>
    033     * @version $Id: ExponentialDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    034     */
    035    public class ExponentialDistribution extends AbstractRealDistribution {
    036        /**
    037         * Default inverse cumulative probability accuracy.
    038         * @since 2.1
    039         */
    040        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    041        /** Serializable version identifier */
    042        private static final long serialVersionUID = 2401296428283614780L;
    043        /**
    044         * Used when generating Exponential samples.
    045         * Table containing the constants
    046         * q_i = sum_{j=1}^i (ln 2)^j/j! = ln 2 + (ln 2)^2/2 + ... + (ln 2)^i/i!
    047         * until the largest representable fraction below 1 is exceeded.
    048         *
    049         * Note that
    050         * 1 = 2 - 1 = exp(ln 2) - 1 = sum_{n=1}^infty (ln 2)^n / n!
    051         * thus q_i -> 1 as i -> +inf,
    052         * so the higher i, the closer to one we get (the series is not alternating).
    053         *
    054         * By trying, n = 16 in Java is enough to reach 1.0.
    055         */
    056        private static final double[] EXPONENTIAL_SA_QI;
    057        /** The mean of this distribution. */
    058        private final double mean;
    059        /** Inverse cumulative probability accuracy. */
    060        private final double solverAbsoluteAccuracy;
    061
    062        /**
    063         * Initialize tables.
    064         */
    065        static {
    066            /**
    067             * Filling EXPONENTIAL_SA_QI table.
    068             * Note that we don't want qi = 0 in the table.
    069             */
    070            final double LN2 = FastMath.log(2);
    071            double qi = 0;
    072            int i = 1;
    073
    074            /**
    075             * ArithmeticUtils provides factorials up to 20, so let's use that
    076             * limit together with Precision.EPSILON to generate the following
    077             * code (a priori, we know that there will be 16 elements, but it is
    078             * better to not hardcode it).
    079             */
    080            final ResizableDoubleArray ra = new ResizableDoubleArray(20);
    081
    082            while (qi < 1) {
    083                qi += FastMath.pow(LN2, i) / ArithmeticUtils.factorial(i);
    084                ra.addElement(qi);
    085                ++i;
    086            }
    087
    088            EXPONENTIAL_SA_QI = ra.getElements();
    089        }
    090
    091        /**
    092         * Create an exponential distribution with the given mean.
    093         * @param mean mean of this distribution.
    094         */
    095        public ExponentialDistribution(double mean) {
    096            this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    097        }
    098
    099        /**
    100         * Create an exponential distribution with the given mean.
    101         *
    102         * @param mean Mean of this distribution.
    103         * @param inverseCumAccuracy Maximum absolute error in inverse
    104         * cumulative probability estimates (defaults to
    105         * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    106         * @throws NotStrictlyPositiveException if {@code mean <= 0}.
    107         * @since 2.1
    108         */
    109        public ExponentialDistribution(double mean, double inverseCumAccuracy) {
    110            this(new Well19937c(), mean, inverseCumAccuracy);
    111        }
    112
    113        /**
    114         * Creates an exponential distribution.
    115         *
    116         * @param rng Random number generator.
    117         * @param mean Mean of this distribution.
    118         * @param inverseCumAccuracy Maximum absolute error in inverse
    119         * cumulative probability estimates (defaults to
    120         * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    121         * @throws NotStrictlyPositiveException if {@code mean <= 0}.
    122         * @since 3.1
    123         */
    124        public ExponentialDistribution(RandomGenerator rng,
    125                                       double mean,
    126                                       double inverseCumAccuracy)
    127            throws NotStrictlyPositiveException {
    128            super(rng);
    129
    130            if (mean <= 0) {
    131                throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
    132            }
    133            this.mean = mean;
    134            solverAbsoluteAccuracy = inverseCumAccuracy;
    135        }
    136
    137        /**
    138         * Access the mean.
    139         *
    140         * @return the mean.
    141         */
    142        public double getMean() {
    143            return mean;
    144        }
    145
    146        /** {@inheritDoc} */
    147        public double density(double x) {
    148            if (x < 0) {
    149                return 0;
    150            }
    151            return FastMath.exp(-x / mean) / mean;
    152        }
    153
    154        /**
    155         * {@inheritDoc}
    156         *
    157         * The implementation of this method is based on:
    158         * <ul>
    159         * <li>
    160         * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
    161         * Exponential Distribution</a>, equation (1).</li>
    162         * </ul>
    163         */
    164        public double cumulativeProbability(double x)  {
    165            double ret;
    166            if (x <= 0.0) {
    167                ret = 0.0;
    168            } else {
    169                ret = 1.0 - FastMath.exp(-x / mean);
    170            }
    171            return ret;
    172        }
    173
    174        /**
    175         * {@inheritDoc}
    176         *
    177         * Returns {@code 0} when {@code p= = 0} and
    178         * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
    179         */
    180        @Override
    181        public double inverseCumulativeProbability(double p) throws OutOfRangeException {
    182            double ret;
    183
    184            if (p < 0.0 || p > 1.0) {
    185                throw new OutOfRangeException(p, 0.0, 1.0);
    186            } else if (p == 1.0) {
    187                ret = Double.POSITIVE_INFINITY;
    188            } else {
    189                ret = -mean * FastMath.log(1.0 - p);
    190            }
    191
    192            return ret;
    193        }
    194
    195        /**
    196         * {@inheritDoc}
    197         *
    198         * <p><strong>Algorithm Description</strong>: this implementation uses the
    199         * <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html">
    200         * Inversion Method</a> to generate exponentially distributed random values
    201         * from uniform deviates.</p>
    202         *
    203         * @return a random value.
    204         * @since 2.2
    205         */
    206        @Override
    207        public double sample() {
    208            // Step 1:
    209            double a = 0;
    210            double u = random.nextDouble();
    211
    212            // Step 2 and 3:
    213            while (u < 0.5) {
    214                a += EXPONENTIAL_SA_QI[0];
    215                u *= 2;
    216            }
    217
    218            // Step 4 (now u >= 0.5):
    219            u += u - 1;
    220
    221            // Step 5:
    222            if (u <= EXPONENTIAL_SA_QI[0]) {
    223                return mean * (a + u);
    224            }
    225
    226            // Step 6:
    227            int i = 0; // Should be 1, be we iterate before it in while using 0
    228            double u2 = random.nextDouble();
    229            double umin = u2;
    230
    231            // Step 7 and 8:
    232            do {
    233                ++i;
    234                u2 = random.nextDouble();
    235
    236                if (u2 < umin) {
    237                    umin = u2;
    238                }
    239
    240                // Step 8:
    241            } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1
    242
    243            return mean * (a + umin * EXPONENTIAL_SA_QI[0]);
    244        }
    245
    246        /** {@inheritDoc} */
    247        @Override
    248        protected double getSolverAbsoluteAccuracy() {
    249            return solverAbsoluteAccuracy;
    250        }
    251
    252        /**
    253         * {@inheritDoc}
    254         *
    255         * For mean parameter {@code k}, the mean is {@code k}.
    256         */
    257        public double getNumericalMean() {
    258            return getMean();
    259        }
    260
    261        /**
    262         * {@inheritDoc}
    263         *
    264         * For mean parameter {@code k}, the variance is {@code k^2}.
    265         */
    266        public double getNumericalVariance() {
    267            final double m = getMean();
    268            return m * m;
    269        }
    270
    271        /**
    272         * {@inheritDoc}
    273         *
    274         * The lower bound of the support is always 0 no matter the mean parameter.
    275         *
    276         * @return lower bound of the support (always 0)
    277         */
    278        public double getSupportLowerBound() {
    279            return 0;
    280        }
    281
    282        /**
    283         * {@inheritDoc}
    284         *
    285         * The upper bound of the support is always positive infinity
    286         * no matter the mean parameter.
    287         *
    288         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
    289         */
    290        public double getSupportUpperBound() {
    291            return Double.POSITIVE_INFINITY;
    292        }
    293
    294        /** {@inheritDoc} */
    295        public boolean isSupportLowerBoundInclusive() {
    296            return true;
    297        }
    298
    299        /** {@inheritDoc} */
    300        public boolean isSupportUpperBoundInclusive() {
    301            return false;
    302        }
    303
    304        /**
    305         * {@inheritDoc}
    306         *
    307         * The support of this distribution is connected.
    308         *
    309         * @return {@code true}
    310         */
    311        public boolean isSupportConnected() {
    312            return true;
    313        }
    314    }