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 < 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 < 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 < 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 < <i>lower bound</i>) < <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 < <i>upper bound</i>) > <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 < μ) > .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 < μ) > .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 }