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.statistics.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
021
022/**
023 * Implementation of the exponential distribution.
024 *
025 * <p>The probability density function of \( X \) is:
026 *
027 * <p>\[ f(x; \mu) = \frac{1}{\mu} e^{-x / \mu} \]
028 *
029 * <p>for \( \mu &gt; 0 \) the mean and
030 * \( x \in [0, \infty) \).
031 *
032 * <p>This implementation uses the scale parameter \( \mu \) which is the mean of the distribution.
033 * A common alternative parameterization uses the rate parameter \( \lambda \) which is the reciprocal
034 * of the mean. The distribution can be be created using \( \mu  = \frac{1}{\lambda} \).
035 *
036 * @see <a href="https://en.wikipedia.org/wiki/Exponential_distribution">Exponential distribution (Wikipedia)</a>
037 * @see <a href="https://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a>
038 */
039public final class ExponentialDistribution extends AbstractContinuousDistribution {
040    /** Support lower bound. */
041    private static final double SUPPORT_LO = 0;
042    /** Support upper bound. */
043    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
044    /** The mean of this distribution. */
045    private final double mean;
046    /** The logarithm of the mean, stored to reduce computing time. */
047    private final double logMean;
048
049    /**
050     * @param mean Mean of this distribution.
051     */
052    private ExponentialDistribution(double mean) {
053        this.mean = mean;
054        logMean = Math.log(mean);
055    }
056
057    /**
058     * Creates an exponential distribution.
059     *
060     * @param mean Mean of this distribution. This is a scale parameter.
061     * @return the distribution
062     * @throws IllegalArgumentException if {@code mean <= 0}.
063     */
064    public static ExponentialDistribution of(double mean) {
065        if (mean <= 0) {
066            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, mean);
067        }
068        return new ExponentialDistribution(mean);
069    }
070
071    /** {@inheritDoc} */
072    @Override
073    public double density(double x) {
074        if (x < SUPPORT_LO) {
075            return 0;
076        }
077        return Math.exp(-x / mean) / mean;
078    }
079
080    /** {@inheritDoc} **/
081    @Override
082    public double logDensity(double x) {
083        if (x < SUPPORT_LO) {
084            return Double.NEGATIVE_INFINITY;
085        }
086        return -x / mean - logMean;
087    }
088
089    /** {@inheritDoc} */
090    @Override
091    public double cumulativeProbability(double x)  {
092        if (x <= SUPPORT_LO) {
093            return 0;
094        }
095        return -Math.expm1(-x / mean);
096    }
097
098    /** {@inheritDoc} */
099    @Override
100    public double survivalProbability(double x)  {
101        if (x <= SUPPORT_LO) {
102            return 1;
103        }
104        return Math.exp(-x / mean);
105    }
106
107    /**
108     * {@inheritDoc}
109     *
110     * <p>Returns {@code 0} when {@code p == 0} and
111     * {@link Double#POSITIVE_INFINITY} when {@code p == 1}.
112     */
113    @Override
114    public double inverseCumulativeProbability(double p) {
115        ArgumentUtils.checkProbability(p);
116        if (p == 1) {
117            return Double.POSITIVE_INFINITY;
118        }
119        // Subtract from zero to prevent returning -0.0 for p=-0.0
120        return 0 - mean * Math.log1p(-p);
121    }
122
123    /**
124     * {@inheritDoc}
125     *
126     * <p>Returns {@code 0} when {@code p == 1} and
127     * {@link Double#POSITIVE_INFINITY} when {@code p == 0}.
128     */
129    @Override
130    public double inverseSurvivalProbability(double p) {
131        ArgumentUtils.checkProbability(p);
132        if (p == 0) {
133            return Double.POSITIVE_INFINITY;
134        }
135        // Subtract from zero to prevent returning -0.0 for p=1
136        return 0 - mean * Math.log(p);
137    }
138
139    /** {@inheritDoc} */
140    @Override
141    public double getMean() {
142        return mean;
143    }
144
145    /**
146     * {@inheritDoc}
147     *
148     * <p>For mean \( \mu \), the variance is \( \mu^2 \).
149     */
150    @Override
151    public double getVariance() {
152        return mean * mean;
153    }
154
155    /**
156     * {@inheritDoc}
157     *
158     * <p>The lower bound of the support is always 0.
159     *
160     * @return 0.
161     */
162    @Override
163    public double getSupportLowerBound() {
164        return SUPPORT_LO;
165    }
166
167    /**
168     * {@inheritDoc}
169     *
170     * <p>The upper bound of the support is always positive infinity.
171     *
172     * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
173     */
174    @Override
175    public double getSupportUpperBound() {
176        return SUPPORT_HI;
177    }
178
179    /** {@inheritDoc} */
180    @Override
181    double getMedian() {
182        // Overridden for the probability(double, double) method.
183        // This is intentionally not a public method.
184        // ln(2) / rate = mean * ln(2)
185        return mean * Constants.LN_TWO;
186    }
187
188    /** {@inheritDoc} */
189    @Override
190    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
191        // Exponential distribution sampler.
192        return ZigguratSampler.Exponential.of(rng, getMean())::sample;
193    }
194}