ExponentialDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  19. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  20. import org.apache.commons.math3.exception.OutOfRangeException;
  21. import org.apache.commons.math3.exception.util.LocalizedFormats;
  22. import org.apache.commons.math3.random.RandomGenerator;
  23. import org.apache.commons.math3.random.Well19937c;
  24. import org.apache.commons.math3.util.CombinatoricsUtils;
  25. import org.apache.commons.math3.util.FastMath;
  26. import org.apache.commons.math3.util.ResizableDoubleArray;

  28. /**
  29.  * Implementation of the exponential distribution.
  30.  *
  31.  * @see <a href="http://en.wikipedia.org/wiki/Exponential_distribution">Exponential distribution (Wikipedia)</a>
  32.  * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a>
  33.  */
  34. public class ExponentialDistribution extends AbstractRealDistribution {
  35.     /**
  36.      * Default inverse cumulative probability accuracy.
  37.      * @since 2.1
  38.      */
  39.     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
  40.     /** Serializable version identifier */
  41.     private static final long serialVersionUID = 2401296428283614780L;
  42.     /**
  43.      * Used when generating Exponential samples.
  44.      * Table containing the constants
  45.      * q_i = sum_{j=1}^i (ln 2)^j/j! = ln 2 + (ln 2)^2/2 + ... + (ln 2)^i/i!
  46.      * until the largest representable fraction below 1 is exceeded.
  47.      *
  48.      * Note that
  49.      * 1 = 2 - 1 = exp(ln 2) - 1 = sum_{n=1}^infty (ln 2)^n / n!
  50.      * thus q_i -> 1 as i -> +inf,
  51.      * so the higher i, the closer to one we get (the series is not alternating).
  52.      *
  53.      * By trying, n = 16 in Java is enough to reach 1.0.
  54.      */
  55.     private static final double[] EXPONENTIAL_SA_QI;
  56.     /** The mean of this distribution. */
  57.     private final double mean;
  58.     /** The logarithm of the mean, stored to reduce computing time. **/
  59.     private final double logMean;
  60.     /** Inverse cumulative probability accuracy. */
  61.     private final double solverAbsoluteAccuracy;

  63.     /**
  64.      * Initialize tables.
  65.      */
  66.     static {
  67.         /**
  68.          * Filling EXPONENTIAL_SA_QI table.
  69.          * Note that we don't want qi = 0 in the table.
  70.          */
  71.         final double LN2 = FastMath.log(2);
  72.         double qi = 0;
  73.         int i = 1;

  75.         /**
  76.          * ArithmeticUtils provides factorials up to 20, so let's use that
  77.          * limit together with Precision.EPSILON to generate the following
  78.          * code (a priori, we know that there will be 16 elements, but it is
  79.          * better to not hardcode it).
  80.          */
  81.         final ResizableDoubleArray ra = new ResizableDoubleArray(20);

  83.         while (qi < 1) {
  84.             qi += FastMath.pow(LN2, i) / CombinatoricsUtils.factorial(i);
  85.             ra.addElement(qi);
  86.             ++i;
  87.         }

  89.         EXPONENTIAL_SA_QI = ra.getElements();
  90.     }

  92.     /**
  93.      * Create an exponential distribution with the given mean.
  94.      * <p>
  95.      * <b>Note:</b> this constructor will implicitly create an instance of
  96.      * {@link Well19937c} as random generator to be used for sampling only (see
  97.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  98.      * needed for the created distribution, it is advised to pass {@code null}
  99.      * as random generator via the appropriate constructors to avoid the
  100.      * additional initialisation overhead.
  101.      *
  102.      * @param mean mean of this distribution.
  103.      */
  104.     public ExponentialDistribution(double mean) {
  105.         this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  106.     }

  108.     /**
  109.      * Create an exponential distribution with the given mean.
  110.      * <p>
  111.      * <b>Note:</b> this constructor will implicitly create an instance of
  112.      * {@link Well19937c} as random generator to be used for sampling only (see
  113.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  114.      * needed for the created distribution, it is advised to pass {@code null}
  115.      * as random generator via the appropriate constructors to avoid the
  116.      * additional initialisation overhead.
  117.      *
  118.      * @param mean Mean of this distribution.
  119.      * @param inverseCumAccuracy Maximum absolute error in inverse
  120.      * cumulative probability estimates (defaults to
  121.      * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  122.      * @throws NotStrictlyPositiveException if {@code mean <= 0}.
  123.      * @since 2.1
  124.      */
  125.     public ExponentialDistribution(double mean, double inverseCumAccuracy) {
  126.         this(new Well19937c(), mean, inverseCumAccuracy);
  127.     }

  129.     /**
  130.      * Creates an exponential distribution.
  131.      *
  132.      * @param rng Random number generator.
  133.      * @param mean Mean of this distribution.
  134.      * @throws NotStrictlyPositiveException if {@code mean <= 0}.
  135.      * @since 3.3
  136.      */
  137.     public ExponentialDistribution(RandomGenerator rng, double mean)
  138.         throws NotStrictlyPositiveException {
  139.         this(rng, mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  140.     }

  142.     /**
  143.      * Creates an exponential distribution.
  144.      *
  145.      * @param rng Random number generator.
  146.      * @param mean Mean of this distribution.
  147.      * @param inverseCumAccuracy Maximum absolute error in inverse
  148.      * cumulative probability estimates (defaults to
  149.      * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  150.      * @throws NotStrictlyPositiveException if {@code mean <= 0}.
  151.      * @since 3.1
  152.      */
  153.     public ExponentialDistribution(RandomGenerator rng,
  154.                                    double mean,
  155.                                    double inverseCumAccuracy)
  156.         throws NotStrictlyPositiveException {
  157.         super(rng);

  159.         if (mean <= 0) {
  160.             throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
  161.         }
  162.         this.mean = mean;
  163.         logMean = FastMath.log(mean);
  164.         solverAbsoluteAccuracy = inverseCumAccuracy;
  165.     }

  167.     /**
  168.      * Access the mean.
  169.      *
  170.      * @return the mean.
  171.      */
  172.     public double getMean() {
  173.         return mean;
  174.     }

  176.     /** {@inheritDoc} */
  177.     public double density(double x) {
  178.         final double logDensity = logDensity(x);
  179.         return logDensity == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logDensity);
  180.     }

  182.     /** {@inheritDoc} **/
  183.     @Override
  184.     public double logDensity(double x) {
  185.         if (x < 0) {
  186.             return Double.NEGATIVE_INFINITY;
  187.         }
  188.         return -x / mean - logMean;
  189.     }

  191.     /**
  192.      * {@inheritDoc}
  193.      *
  194.      * The implementation of this method is based on:
  195.      * <ul>
  196.      * <li>
  197.      * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
  198.      * Exponential Distribution</a>, equation (1).</li>
  199.      * </ul>
  200.      */
  201.     public double cumulativeProbability(double x)  {
  202.         double ret;
  203.         if (x <= 0.0) {
  204.             ret = 0.0;
  205.         } else {
  206.             ret = 1.0 - FastMath.exp(-x / mean);
  207.         }
  208.         return ret;
  209.     }

  211.     /**
  212.      * {@inheritDoc}
  213.      *
  214.      * Returns {@code 0} when {@code p= = 0} and
  215.      * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
  216.      */
  217.     @Override
  218.     public double inverseCumulativeProbability(double p) throws OutOfRangeException {
  219.         double ret;

  221.         if (p < 0.0 || p > 1.0) {
  222.             throw new OutOfRangeException(p, 0.0, 1.0);
  223.         } else if (p == 1.0) {
  224.             ret = Double.POSITIVE_INFINITY;
  225.         } else {
  226.             ret = -mean * FastMath.log(1.0 - p);
  227.         }

  229.         return ret;
  230.     }

  232.     /**
  233.      * {@inheritDoc}
  234.      *
  235.      * <p><strong>Algorithm Description</strong>: this implementation uses the
  236.      * <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html">
  237.      * Inversion Method</a> to generate exponentially distributed random values
  238.      * from uniform deviates.</p>
  239.      *
  240.      * @return a random value.
  241.      * @since 2.2
  242.      */
  243.     @Override
  244.     public double sample() {
  245.         // Step 1:
  246.         double a = 0;
  247.         double u = random.nextDouble();

  249.         // Step 2 and 3:
  250.         while (u < 0.5) {
  251.             a += EXPONENTIAL_SA_QI[0];
  252.             u *= 2;
  253.         }

  255.         // Step 4 (now u >= 0.5):
  256.         u += u - 1;

  258.         // Step 5:
  259.         if (u <= EXPONENTIAL_SA_QI[0]) {
  260.             return mean * (a + u);
  261.         }

  263.         // Step 6:
  264.         int i = 0; // Should be 1, be we iterate before it in while using 0
  265.         double u2 = random.nextDouble();
  266.         double umin = u2;

  268.         // Step 7 and 8:
  269.         do {
  270.             ++i;
  271.             u2 = random.nextDouble();

  273.             if (u2 < umin) {
  274.                 umin = u2;
  275.             }

  277.             // Step 8:
  278.         } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1

  280.         return mean * (a + umin * EXPONENTIAL_SA_QI[0]);
  281.     }

  283.     /** {@inheritDoc} */
  284.     @Override
  285.     protected double getSolverAbsoluteAccuracy() {
  286.         return solverAbsoluteAccuracy;
  287.     }

  289.     /**
  290.      * {@inheritDoc}
  291.      *
  292.      * For mean parameter {@code k}, the mean is {@code k}.
  293.      */
  294.     public double getNumericalMean() {
  295.         return getMean();
  296.     }

  298.     /**
  299.      * {@inheritDoc}
  300.      *
  301.      * For mean parameter {@code k}, the variance is {@code k^2}.
  302.      */
  303.     public double getNumericalVariance() {
  304.         final double m = getMean();
  305.         return m * m;
  306.     }

  308.     /**
  309.      * {@inheritDoc}
  310.      *
  311.      * The lower bound of the support is always 0 no matter the mean parameter.
  312.      *
  313.      * @return lower bound of the support (always 0)
  314.      */
  315.     public double getSupportLowerBound() {
  316.         return 0;
  317.     }

  319.     /**
  320.      * {@inheritDoc}
  321.      *
  322.      * The upper bound of the support is always positive infinity
  323.      * no matter the mean parameter.
  324.      *
  325.      * @return upper bound of the support (always Double.POSITIVE_INFINITY)
  326.      */
  327.     public double getSupportUpperBound() {
  328.         return Double.POSITIVE_INFINITY;
  329.     }

  331.     /** {@inheritDoc} */
  332.     public boolean isSupportLowerBoundInclusive() {
  333.         return true;
  334.     }

  336.     /** {@inheritDoc} */
  337.     public boolean isSupportUpperBoundInclusive() {
  338.         return false;
  339.     }

  341.     /**
  342.      * {@inheritDoc}
  343.      *
  344.      * The support of this distribution is connected.
  345.      *
  346.      * @return {@code true}
  347.      */
  348.     public boolean isSupportConnected() {
  349.         return true;
  350.     }
  351. }
Created with JaCoCo 0.7.2.201409121644