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.AhrensDieterExponentialSampler;
021
022/**
023 * Implementation of the <a href="http://en.wikipedia.org/wiki/Exponential_distribution">exponential distribution</a>.
024 */
025public class ExponentialDistribution extends AbstractContinuousDistribution {
026    /** Support lower bound. */
027    private static final double SUPPORT_LO = 0;
028    /** Support upper bound. */
029    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
030    /** The mean of this distribution. */
031    private final double mean;
032    /** The logarithm of the mean, stored to reduce computing time. */
033    private final double logMean;
034
035    /**
036     * Creates a distribution.
037     *
038     * @param mean Mean of this distribution.
039     * @throws IllegalArgumentException if {@code mean <= 0}.
040     */
041    public ExponentialDistribution(double mean) {
042        if (mean <= 0) {
043            throw new DistributionException(DistributionException.NEGATIVE, mean);
044        }
045        this.mean = mean;
046        logMean = Math.log(mean);
047    }
048
049    /** {@inheritDoc} */
050    @Override
051    public double density(double x) {
052        final double logDensity = logDensity(x);
053        return logDensity == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logDensity);
054    }
055
056    /** {@inheritDoc} **/
057    @Override
058    public double logDensity(double x) {
059        if (x < SUPPORT_LO ||
060            x >= SUPPORT_HI) {
061            return Double.NEGATIVE_INFINITY;
062        }
063
064        return -x / mean - logMean;
065    }
066
067    /**
068     * {@inheritDoc}
069     *
070     * The implementation of this method is based on:
071     * <ul>
072     * <li>
073     * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
074     * Exponential Distribution</a>, equation (1).</li>
075     * </ul>
076     */
077    @Override
078    public double cumulativeProbability(double x)  {
079        if (x <= SUPPORT_LO) {
080            return 0;
081        }
082
083        return 1 - Math.exp(-x / mean);
084    }
085
086    /**
087     * {@inheritDoc}
088     *
089     * Returns {@code 0} when {@code p= = 0} and
090     * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
091     */
092    @Override
093    public double inverseCumulativeProbability(double p) {
094        double ret;
095
096        if (p < 0 ||
097            p > 1) {
098            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
099        } else if (p == 1) {
100            ret = Double.POSITIVE_INFINITY;
101        } else {
102            ret = -mean * Math.log(1 - p);
103        }
104
105        return ret;
106    }
107
108    /** {@inheritDoc} */
109    @Override
110    public double getMean() {
111        return mean;
112    }
113
114    /**
115     * {@inheritDoc}
116     *
117     * For mean parameter {@code k}, the variance is {@code k^2}.
118     */
119    @Override
120    public double getVariance() {
121        return mean * mean;
122    }
123
124    /**
125     * {@inheritDoc}
126     *
127     * The lower bound of the support is always 0 no matter the mean parameter.
128     *
129     * @return lower bound of the support (always 0)
130     */
131    @Override
132    public double getSupportLowerBound() {
133        return SUPPORT_LO;
134    }
135
136    /**
137     * {@inheritDoc}
138     *
139     * The upper bound of the support is always positive infinity
140     * no matter the mean parameter.
141     *
142     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
143     */
144    @Override
145    public double getSupportUpperBound() {
146        return SUPPORT_HI;
147    }
148
149    /**
150     * {@inheritDoc}
151     *
152     * The support of this distribution is connected.
153     *
154     * @return {@code true}
155     */
156    @Override
157    public boolean isSupportConnected() {
158        return true;
159    }
160
161    /**
162     * {@inheritDoc}
163     *
164     * <p>Sampling algorithm uses the
165     *  <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html">
166     *   inversion method</a> to generate exponentially distributed
167     *  random values from uniform deviates.
168     * </p>
169     */
170    @Override
171    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
172        // Exponential distribution sampler.
173        return new AhrensDieterExponentialSampler(rng, mean)::sample;
174    }
175}