1 | /*
|
---|
2 | * Licensed to the Apache Software Foundation (ASF) under one or more
|
---|
3 | * contributor license agreements. See the NOTICE file distributed with
|
---|
4 | * this work for additional information regarding copyright ownership.
|
---|
5 | * The ASF licenses this file to You under the Apache License, Version 2.0
|
---|
6 | * (the "License"); you may not use this file except in compliance with
|
---|
7 | * the License. You may obtain a copy of the License at
|
---|
8 | *
|
---|
9 | * http://www.apache.org/licenses/LICENSE-2.0
|
---|
10 | *
|
---|
11 | * Unless required by applicable law or agreed to in writing, software
|
---|
12 | * distributed under the License is distributed on an "AS IS" BASIS,
|
---|
13 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
---|
14 | * See the License for the specific language governing permissions and
|
---|
15 | * limitations under the License.
|
---|
16 | */
|
---|
17 | package agents.org.apache.commons.math.distribution;
|
---|
18 |
|
---|
19 | import java.io.Serializable;
|
---|
20 |
|
---|
21 | import agents.org.apache.commons.math.MathException;
|
---|
22 | import agents.org.apache.commons.math.MathRuntimeException;
|
---|
23 | import agents.org.apache.commons.math.exception.util.LocalizedFormats;
|
---|
24 | import agents.org.apache.commons.math.util.FastMath;
|
---|
25 |
|
---|
26 | /**
|
---|
27 | * The default implementation of {@link ExponentialDistribution}.
|
---|
28 | *
|
---|
29 | * @version $Revision: 1055914 $ $Date: 2011-01-06 16:34:34 +0100 (jeu. 06 janv. 2011) $
|
---|
30 | */
|
---|
31 | public class ExponentialDistributionImpl extends AbstractContinuousDistribution
|
---|
32 | implements ExponentialDistribution, Serializable {
|
---|
33 |
|
---|
34 | /**
|
---|
35 | * Default inverse cumulative probability accuracy
|
---|
36 | * @since 2.1
|
---|
37 | */
|
---|
38 | public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
|
---|
39 |
|
---|
40 | /** Serializable version identifier */
|
---|
41 | private static final long serialVersionUID = 2401296428283614780L;
|
---|
42 |
|
---|
43 | /** The mean of this distribution. */
|
---|
44 | private double mean;
|
---|
45 |
|
---|
46 | /** Inverse cumulative probability accuracy */
|
---|
47 | private final double solverAbsoluteAccuracy;
|
---|
48 |
|
---|
49 | /**
|
---|
50 | * Create a exponential distribution with the given mean.
|
---|
51 | * @param mean mean of this distribution.
|
---|
52 | */
|
---|
53 | public ExponentialDistributionImpl(double mean) {
|
---|
54 | this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
|
---|
55 | }
|
---|
56 |
|
---|
57 | /**
|
---|
58 | * Create a exponential distribution with the given mean.
|
---|
59 | * @param mean mean of this distribution.
|
---|
60 | * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
|
---|
61 | * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
|
---|
62 | * @since 2.1
|
---|
63 | */
|
---|
64 | public ExponentialDistributionImpl(double mean, double inverseCumAccuracy) {
|
---|
65 | super();
|
---|
66 | setMeanInternal(mean);
|
---|
67 | solverAbsoluteAccuracy = inverseCumAccuracy;
|
---|
68 | }
|
---|
69 |
|
---|
70 | /**
|
---|
71 | * Modify the mean.
|
---|
72 | * @param mean the new mean.
|
---|
73 | * @throws IllegalArgumentException if <code>mean</code> is not positive.
|
---|
74 | * @deprecated as of 2.1 (class will become immutable in 3.0)
|
---|
75 | */
|
---|
76 | @Deprecated
|
---|
77 | public void setMean(double mean) {
|
---|
78 | setMeanInternal(mean);
|
---|
79 | }
|
---|
80 | /**
|
---|
81 | * Modify the mean.
|
---|
82 | * @param newMean the new mean.
|
---|
83 | * @throws IllegalArgumentException if <code>newMean</code> is not positive.
|
---|
84 | */
|
---|
85 | private void setMeanInternal(double newMean) {
|
---|
86 | if (newMean <= 0.0) {
|
---|
87 | throw MathRuntimeException.createIllegalArgumentException(
|
---|
88 | LocalizedFormats.NOT_POSITIVE_MEAN, newMean);
|
---|
89 | }
|
---|
90 | this.mean = newMean;
|
---|
91 | }
|
---|
92 |
|
---|
93 | /**
|
---|
94 | * Access the mean.
|
---|
95 | * @return the mean.
|
---|
96 | */
|
---|
97 | public double getMean() {
|
---|
98 | return mean;
|
---|
99 | }
|
---|
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 | }
|
---|