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.math3.distribution;
018
019 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020 import org.apache.commons.math3.exception.OutOfRangeException;
021 import org.apache.commons.math3.exception.util.LocalizedFormats;
022 import org.apache.commons.math3.util.FastMath;
023 import org.apache.commons.math3.util.ArithmeticUtils;
024 import org.apache.commons.math3.util.ResizableDoubleArray;
025 import org.apache.commons.math3.random.RandomGenerator;
026 import org.apache.commons.math3.random.Well19937c;
027
028 /**
029 * Implementation of the exponential distribution.
030 *
031 * @see <a href="http://en.wikipedia.org/wiki/Exponential_distribution">Exponential distribution (Wikipedia)</a>
032 * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential distribution (MathWorld)</a>
033 * @version $Id: ExponentialDistribution.java 1416643 2012-12-03 19:37:14Z tn $
034 */
035 public class ExponentialDistribution extends AbstractRealDistribution {
036 /**
037 * Default inverse cumulative probability accuracy.
038 * @since 2.1
039 */
040 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
041 /** Serializable version identifier */
042 private static final long serialVersionUID = 2401296428283614780L;
043 /**
044 * Used when generating Exponential samples.
045 * Table containing the constants
046 * q_i = sum_{j=1}^i (ln 2)^j/j! = ln 2 + (ln 2)^2/2 + ... + (ln 2)^i/i!
047 * until the largest representable fraction below 1 is exceeded.
048 *
049 * Note that
050 * 1 = 2 - 1 = exp(ln 2) - 1 = sum_{n=1}^infty (ln 2)^n / n!
051 * thus q_i -> 1 as i -> +inf,
052 * so the higher i, the closer to one we get (the series is not alternating).
053 *
054 * By trying, n = 16 in Java is enough to reach 1.0.
055 */
056 private static final double[] EXPONENTIAL_SA_QI;
057 /** The mean of this distribution. */
058 private final double mean;
059 /** Inverse cumulative probability accuracy. */
060 private final double solverAbsoluteAccuracy;
061
062 /**
063 * Initialize tables.
064 */
065 static {
066 /**
067 * Filling EXPONENTIAL_SA_QI table.
068 * Note that we don't want qi = 0 in the table.
069 */
070 final double LN2 = FastMath.log(2);
071 double qi = 0;
072 int i = 1;
073
074 /**
075 * ArithmeticUtils provides factorials up to 20, so let's use that
076 * limit together with Precision.EPSILON to generate the following
077 * code (a priori, we know that there will be 16 elements, but it is
078 * better to not hardcode it).
079 */
080 final ResizableDoubleArray ra = new ResizableDoubleArray(20);
081
082 while (qi < 1) {
083 qi += FastMath.pow(LN2, i) / ArithmeticUtils.factorial(i);
084 ra.addElement(qi);
085 ++i;
086 }
087
088 EXPONENTIAL_SA_QI = ra.getElements();
089 }
090
091 /**
092 * Create an exponential distribution with the given mean.
093 * @param mean mean of this distribution.
094 */
095 public ExponentialDistribution(double mean) {
096 this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
097 }
098
099 /**
100 * Create an exponential distribution with the given mean.
101 *
102 * @param mean Mean of this distribution.
103 * @param inverseCumAccuracy Maximum absolute error in inverse
104 * cumulative probability estimates (defaults to
105 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
106 * @throws NotStrictlyPositiveException if {@code mean <= 0}.
107 * @since 2.1
108 */
109 public ExponentialDistribution(double mean, double inverseCumAccuracy) {
110 this(new Well19937c(), mean, inverseCumAccuracy);
111 }
112
113 /**
114 * Creates an exponential distribution.
115 *
116 * @param rng Random number generator.
117 * @param mean Mean of this distribution.
118 * @param inverseCumAccuracy Maximum absolute error in inverse
119 * cumulative probability estimates (defaults to
120 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
121 * @throws NotStrictlyPositiveException if {@code mean <= 0}.
122 * @since 3.1
123 */
124 public ExponentialDistribution(RandomGenerator rng,
125 double mean,
126 double inverseCumAccuracy)
127 throws NotStrictlyPositiveException {
128 super(rng);
129
130 if (mean <= 0) {
131 throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
132 }
133 this.mean = mean;
134 solverAbsoluteAccuracy = inverseCumAccuracy;
135 }
136
137 /**
138 * Access the mean.
139 *
140 * @return the mean.
141 */
142 public double getMean() {
143 return mean;
144 }
145
146 /** {@inheritDoc} */
147 public double density(double x) {
148 if (x < 0) {
149 return 0;
150 }
151 return FastMath.exp(-x / mean) / mean;
152 }
153
154 /**
155 * {@inheritDoc}
156 *
157 * The implementation of this method is based on:
158 * <ul>
159 * <li>
160 * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">
161 * Exponential Distribution</a>, equation (1).</li>
162 * </ul>
163 */
164 public double cumulativeProbability(double x) {
165 double ret;
166 if (x <= 0.0) {
167 ret = 0.0;
168 } else {
169 ret = 1.0 - FastMath.exp(-x / mean);
170 }
171 return ret;
172 }
173
174 /**
175 * {@inheritDoc}
176 *
177 * Returns {@code 0} when {@code p= = 0} and
178 * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
179 */
180 @Override
181 public double inverseCumulativeProbability(double p) throws OutOfRangeException {
182 double ret;
183
184 if (p < 0.0 || p > 1.0) {
185 throw new OutOfRangeException(p, 0.0, 1.0);
186 } else if (p == 1.0) {
187 ret = Double.POSITIVE_INFINITY;
188 } else {
189 ret = -mean * FastMath.log(1.0 - p);
190 }
191
192 return ret;
193 }
194
195 /**
196 * {@inheritDoc}
197 *
198 * <p><strong>Algorithm Description</strong>: this implementation uses the
199 * <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html">
200 * Inversion Method</a> to generate exponentially distributed random values
201 * from uniform deviates.</p>
202 *
203 * @return a random value.
204 * @since 2.2
205 */
206 @Override
207 public double sample() {
208 // Step 1:
209 double a = 0;
210 double u = random.nextDouble();
211
212 // Step 2 and 3:
213 while (u < 0.5) {
214 a += EXPONENTIAL_SA_QI[0];
215 u *= 2;
216 }
217
218 // Step 4 (now u >= 0.5):
219 u += u - 1;
220
221 // Step 5:
222 if (u <= EXPONENTIAL_SA_QI[0]) {
223 return mean * (a + u);
224 }
225
226 // Step 6:
227 int i = 0; // Should be 1, be we iterate before it in while using 0
228 double u2 = random.nextDouble();
229 double umin = u2;
230
231 // Step 7 and 8:
232 do {
233 ++i;
234 u2 = random.nextDouble();
235
236 if (u2 < umin) {
237 umin = u2;
238 }
239
240 // Step 8:
241 } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1
242
243 return mean * (a + umin * EXPONENTIAL_SA_QI[0]);
244 }
245
246 /** {@inheritDoc} */
247 @Override
248 protected double getSolverAbsoluteAccuracy() {
249 return solverAbsoluteAccuracy;
250 }
251
252 /**
253 * {@inheritDoc}
254 *
255 * For mean parameter {@code k}, the mean is {@code k}.
256 */
257 public double getNumericalMean() {
258 return getMean();
259 }
260
261 /**
262 * {@inheritDoc}
263 *
264 * For mean parameter {@code k}, the variance is {@code k^2}.
265 */
266 public double getNumericalVariance() {
267 final double m = getMean();
268 return m * m;
269 }
270
271 /**
272 * {@inheritDoc}
273 *
274 * The lower bound of the support is always 0 no matter the mean parameter.
275 *
276 * @return lower bound of the support (always 0)
277 */
278 public double getSupportLowerBound() {
279 return 0;
280 }
281
282 /**
283 * {@inheritDoc}
284 *
285 * The upper bound of the support is always positive infinity
286 * no matter the mean parameter.
287 *
288 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
289 */
290 public double getSupportUpperBound() {
291 return Double.POSITIVE_INFINITY;
292 }
293
294 /** {@inheritDoc} */
295 public boolean isSupportLowerBoundInclusive() {
296 return true;
297 }
298
299 /** {@inheritDoc} */
300 public boolean isSupportUpperBoundInclusive() {
301 return false;
302 }
303
304 /**
305 * {@inheritDoc}
306 *
307 * The support of this distribution is connected.
308 *
309 * @return {@code true}
310 */
311 public boolean isSupportConnected() {
312 return true;
313 }
314 }