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.math.distribution;
018    
019    import java.io.Serializable;
020    
021    import org.apache.commons.math.MathException;
022    import org.apache.commons.math.MathRuntimeException;
023    import org.apache.commons.math.exception.util.LocalizedFormats;
024    import org.apache.commons.math.util.FastMath;
025    
026    /**
027     * The default implementation of {@link ExponentialDistribution}.
028     *
029     * @version $Revision: 1055914 $ $Date: 2011-01-06 16:34:34 +0100 (jeu. 06 janv. 2011) $
030     */
031    public class ExponentialDistributionImpl extends AbstractContinuousDistribution
032        implements ExponentialDistribution, Serializable {
033    
034        /**
035         * Default inverse cumulative probability accuracy
036         * @since 2.1
037         */
038        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
039    
040        /** Serializable version identifier */
041        private static final long serialVersionUID = 2401296428283614780L;
042    
043        /** The mean of this distribution. */
044        private double mean;
045    
046        /** Inverse cumulative probability accuracy */
047        private final double solverAbsoluteAccuracy;
048    
049        /**
050         * Create a exponential distribution with the given mean.
051         * @param mean mean of this distribution.
052         */
053        public ExponentialDistributionImpl(double mean) {
054            this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
055        }
056    
057        /**
058         * Create a exponential distribution with the given mean.
059         * @param mean mean of this distribution.
060         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
061         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
062         * @since 2.1
063         */
064        public ExponentialDistributionImpl(double mean, double inverseCumAccuracy) {
065            super();
066            setMeanInternal(mean);
067            solverAbsoluteAccuracy = inverseCumAccuracy;
068        }
069    
070        /**
071         * Modify the mean.
072         * @param mean the new mean.
073         * @throws IllegalArgumentException if <code>mean</code> is not positive.
074         * @deprecated as of 2.1 (class will become immutable in 3.0)
075         */
076        @Deprecated
077        public void setMean(double mean) {
078            setMeanInternal(mean);
079        }
080        /**
081         * Modify the mean.
082         * @param newMean the new mean.
083         * @throws IllegalArgumentException if <code>newMean</code> is not positive.
084         */
085        private void setMeanInternal(double newMean) {
086            if (newMean <= 0.0) {
087                throw MathRuntimeException.createIllegalArgumentException(
088                      LocalizedFormats.NOT_POSITIVE_MEAN, newMean);
089            }
090            this.mean = newMean;
091        }
092    
093        /**
094         * Access the mean.
095         * @return the mean.
096         */
097        public double getMean() {
098            return mean;
099        }
100    
101        /**
102         * Return the probability density for a particular point.
103         *
104         * @param x The point at which the density should be computed.
105         * @return The pdf at point x.
106         * @deprecated - use density(double)
107         */
108        @Deprecated
109        public double density(Double x) {
110            return density(x.doubleValue());
111        }
112    
113        /**
114         * Return the probability density for a particular point.
115         *
116         * @param x The point at which the density should be computed.
117         * @return The pdf at point x.
118         * @since 2.1
119         */
120        @Override
121        public double density(double x) {
122            if (x < 0) {
123                return 0;
124            }
125            return FastMath.exp(-x / mean) / mean;
126        }
127    
128        /**
129         * For this distribution, X, this method returns P(X &lt; x).
130         *
131         * The implementation of this method is based on:
132         * <ul>
133         * <li>
134         * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
135         * Exponential Distribution</a>, equation (1).</li>
136         * </ul>
137         *
138         * @param x the value at which the CDF is evaluated.
139         * @return CDF for this distribution.
140         * @throws MathException if the cumulative probability can not be
141         *            computed due to convergence or other numerical errors.
142         */
143        public double cumulativeProbability(double x) throws MathException{
144            double ret;
145            if (x <= 0.0) {
146                ret = 0.0;
147            } else {
148                ret = 1.0 - FastMath.exp(-x / mean);
149            }
150            return ret;
151        }
152    
153        /**
154         * For this distribution, X, this method returns the critical point x, such
155         * that P(X &lt; x) = <code>p</code>.
156         * <p>
157         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
158         *
159         * @param p the desired probability
160         * @return x, such that P(X &lt; x) = <code>p</code>
161         * @throws MathException if the inverse cumulative probability can not be
162         *            computed due to convergence or other numerical errors.
163         * @throws IllegalArgumentException if p < 0 or p > 1.
164         */
165        @Override
166        public double inverseCumulativeProbability(double p) throws MathException {
167            double ret;
168    
169            if (p < 0.0 || p > 1.0) {
170                throw MathRuntimeException.createIllegalArgumentException(
171                      LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
172            } else if (p == 1.0) {
173                ret = Double.POSITIVE_INFINITY;
174            } else {
175                ret = -mean * FastMath.log(1.0 - p);
176            }
177    
178            return ret;
179        }
180    
181        /**
182         * Generates a random value sampled from this distribution.
183         *
184         * <p><strong>Algorithm Description</strong>: Uses the <a
185         * href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html"> Inversion
186         * Method</a> to generate exponentially distributed random values from
187         * uniform deviates. </p>
188         *
189         * @return random value
190         * @since 2.2
191         * @throws MathException if an error occurs generating the random value
192         */
193        @Override
194        public double sample() throws MathException {
195            return randomData.nextExponential(mean);
196        }
197    
198        /**
199         * Access the domain value lower bound, based on <code>p</code>, used to
200         * bracket a CDF root.
201         *
202         * @param p the desired probability for the critical value
203         * @return domain value lower bound, i.e.
204         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
205         */
206        @Override
207        protected double getDomainLowerBound(double p) {
208            return 0;
209        }
210    
211        /**
212         * Access the domain value upper bound, based on <code>p</code>, used to
213         * bracket a CDF root.
214         *
215         * @param p the desired probability for the critical value
216         * @return domain value upper bound, i.e.
217         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
218         */
219        @Override
220        protected double getDomainUpperBound(double p) {
221            // NOTE: exponential is skewed to the left
222            // NOTE: therefore, P(X < &mu;) > .5
223    
224            if (p < .5) {
225                // use mean
226                return mean;
227            } else {
228                // use max
229                return Double.MAX_VALUE;
230            }
231        }
232    
233        /**
234         * Access the initial domain value, based on <code>p</code>, used to
235         * bracket a CDF root.
236         *
237         * @param p the desired probability for the critical value
238         * @return initial domain value
239         */
240        @Override
241        protected double getInitialDomain(double p) {
242            // TODO: try to improve on this estimate
243            // TODO: what should really happen here is not derive from AbstractContinuousDistribution
244            // TODO: because the inverse cumulative distribution is simple.
245            // Exponential is skewed to the left, therefore, P(X < &mu;) > .5
246            if (p < .5) {
247                // use 1/2 mean
248                return mean * .5;
249            } else {
250                // use mean
251                return mean;
252            }
253        }
254    
255        /**
256         * Return the absolute accuracy setting of the solver used to estimate
257         * inverse cumulative probabilities.
258         *
259         * @return the solver absolute accuracy
260         * @since 2.1
261         */
262        @Override
263        protected double getSolverAbsoluteAccuracy() {
264            return solverAbsoluteAccuracy;
265        }
266    
267        /**
268         * Returns the lower bound of the support for the distribution.
269         *
270         * The lower bound of the support is always 0, regardless of the mean.
271         *
272         * @return lower bound of the support (always 0)
273         * @since 2.2
274         */
275        public double getSupportLowerBound() {
276            return 0;
277        }
278    
279        /**
280         * Returns the upper bound of the support for the distribution.
281         *
282         * The upper bound of the support is always positive infinity,
283         * regardless of the mean.
284         *
285         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
286         * @since 2.2
287         */
288        public double getSupportUpperBound() {
289            return Double.POSITIVE_INFINITY;
290        }
291    
292        /**
293         * Returns the mean of the distribution.
294         *
295         * For mean parameter <code>k</code>, the mean is
296         * <code>k</code>
297         *
298         * @return the mean
299         * @since 2.2
300         */
301        public double getNumericalMean() {
302            return getMean();
303        }
304    
305        /**
306         * Returns the variance of the distribution.
307         *
308         * For mean parameter <code>k</code>, the variance is
309         * <code>k^2</code>
310         *
311         * @return the variance
312         * @since 2.2
313         */
314        public double getNumericalVariance() {
315            final double m = getMean();
316            return m * m;
317        }
318    
319    }