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 }