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 */ 17package org.apache.commons.math.distribution; 18 19import java.io.Serializable; 20 21import org.apache.commons.math.MathException; 22import org.apache.commons.math.MathRuntimeException; 23import org.apache.commons.math.exception.util.LocalizedFormats; 24import org.apache.commons.math.special.Gamma; 25import org.apache.commons.math.util.MathUtils; 26import org.apache.commons.math.util.FastMath; 27 28/** 29 * Implementation for the {@link PoissonDistribution}. 30 * 31 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 32 */ 33public class PoissonDistributionImpl extends AbstractIntegerDistribution 34 implements PoissonDistribution, Serializable { 35 36 /** 37 * Default maximum number of iterations for cumulative probability calculations. 38 * @since 2.1 39 */ 40 public static final int DEFAULT_MAX_ITERATIONS = 10000000; 41 42 /** 43 * Default convergence criterion. 44 * @since 2.1 45 */ 46 public static final double DEFAULT_EPSILON = 1E-12; 47 48 /** Serializable version identifier */ 49 private static final long serialVersionUID = -3349935121172596109L; 50 51 /** Distribution used to compute normal approximation. */ 52 private NormalDistribution normal; 53 54 /** 55 * Holds the Poisson mean for the distribution. 56 */ 57 private double mean; 58 59 /** 60 * Maximum number of iterations for cumulative probability. 61 * 62 * Cumulative probabilities are estimated using either Lanczos series approximation of 63 * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ. 64 */ 65 private int maxIterations = DEFAULT_MAX_ITERATIONS; 66 67 /** 68 * Convergence criterion for cumulative probability. 69 */ 70 private double epsilon = DEFAULT_EPSILON; 71 72 /** 73 * Create a new Poisson distribution with the given the mean. The mean value 74 * must be positive; otherwise an <code>IllegalArgument</code> is thrown. 75 * 76 * @param p the Poisson mean 77 * @throws IllegalArgumentException if p ≤ 0 78 */ 79 public PoissonDistributionImpl(double p) { 80 this(p, new NormalDistributionImpl()); 81 } 82 83 /** 84 * Create a new Poisson distribution with the given mean, convergence criterion 85 * and maximum number of iterations. 86 * 87 * @param p the Poisson mean 88 * @param epsilon the convergence criteria for cumulative probabilites 89 * @param maxIterations the maximum number of iterations for cumulative probabilites 90 * @since 2.1 91 */ 92 public PoissonDistributionImpl(double p, double epsilon, int maxIterations) { 93 setMean(p); 94 this.epsilon = epsilon; 95 this.maxIterations = maxIterations; 96 } 97 98 /** 99 * Create a new Poisson distribution with the given mean and convergence criterion. 100 * 101 * @param p the Poisson mean 102 * @param epsilon the convergence criteria for cumulative probabilites 103 * @since 2.1 104 */ 105 public PoissonDistributionImpl(double p, double epsilon) { 106 setMean(p); 107 this.epsilon = epsilon; 108 } 109 110 /** 111 * Create a new Poisson distribution with the given mean and maximum number of iterations. 112 * 113 * @param p the Poisson mean 114 * @param maxIterations the maximum number of iterations for cumulative probabilites 115 * @since 2.1 116 */ 117 public PoissonDistributionImpl(double p, int maxIterations) { 118 setMean(p); 119 this.maxIterations = maxIterations; 120 } 121 122 123 /** 124 * Create a new Poisson distribution with the given the mean. The mean value 125 * must be positive; otherwise an <code>IllegalArgument</code> is thrown. 126 * 127 * @param p the Poisson mean 128 * @param z a normal distribution used to compute normal approximations. 129 * @throws IllegalArgumentException if p ≤ 0 130 * @since 1.2 131 * @deprecated as of 2.1 (to avoid possibly inconsistent state, the 132 * "NormalDistribution" will be instantiated internally) 133 */ 134 @Deprecated 135 public PoissonDistributionImpl(double p, NormalDistribution z) { 136 super(); 137 setNormalAndMeanInternal(z, p); 138 } 139 140 /** 141 * Get the Poisson mean for the distribution. 142 * 143 * @return the Poisson mean for the distribution. 144 */ 145 public double getMean() { 146 return mean; 147 } 148 149 /** 150 * Set the Poisson mean for the distribution. The mean value must be 151 * positive; otherwise an <code>IllegalArgument</code> is thrown. 152 * 153 * @param p the Poisson mean value 154 * @throws IllegalArgumentException if p ≤ 0 155 * @deprecated as of 2.1 (class will become immutable in 3.0) 156 */ 157 @Deprecated 158 public void setMean(double p) { 159 setNormalAndMeanInternal(normal, p); 160 } 161 /** 162 * Set the Poisson mean for the distribution. The mean value must be 163 * positive; otherwise an <code>IllegalArgument</code> is thrown. 164 * 165 * @param z the new distribution 166 * @param p the Poisson mean value 167 * @throws IllegalArgumentException if p ≤ 0 168 */ 169 private void setNormalAndMeanInternal(NormalDistribution z, 170 double p) { 171 if (p <= 0) { 172 throw MathRuntimeException.createIllegalArgumentException( 173 LocalizedFormats.NOT_POSITIVE_POISSON_MEAN, p); 174 } 175 mean = p; 176 normal = z; 177 normal.setMean(p); 178 normal.setStandardDeviation(FastMath.sqrt(p)); 179 } 180 181 /** 182 * The probability mass function P(X = x) for a Poisson distribution. 183 * 184 * @param x the value at which the probability density function is 185 * evaluated. 186 * @return the value of the probability mass function at x 187 */ 188 public double probability(int x) { 189 double ret; 190 if (x < 0 || x == Integer.MAX_VALUE) { 191 ret = 0.0; 192 } else if (x == 0) { 193 ret = FastMath.exp(-mean); 194 } else { 195 ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) - 196 SaddlePointExpansion.getDeviancePart(x, mean)) / 197 FastMath.sqrt(MathUtils.TWO_PI * x); 198 } 199 return ret; 200 } 201 202 /** 203 * The probability distribution function P(X <= x) for a Poisson 204 * distribution. 205 * 206 * @param x the value at which the PDF is evaluated. 207 * @return Poisson distribution function evaluated at x 208 * @throws MathException if the cumulative probability can not be computed 209 * due to convergence or other numerical errors. 210 */ 211 @Override 212 public double cumulativeProbability(int x) throws MathException { 213 if (x < 0) { 214 return 0; 215 } 216 if (x == Integer.MAX_VALUE) { 217 return 1; 218 } 219 return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations); 220 } 221 222 /** 223 * Calculates the Poisson distribution function using a normal 224 * approximation. The <code>N(mean, sqrt(mean))</code> distribution is used 225 * to approximate the Poisson distribution. 226 * <p> 227 * The computation uses "half-correction" -- evaluating the normal 228 * distribution function at <code>x + 0.5</code> 229 * </p> 230 * 231 * @param x the upper bound, inclusive 232 * @return the distribution function value calculated using a normal 233 * approximation 234 * @throws MathException if an error occurs computing the normal 235 * approximation 236 */ 237 public double normalApproximateProbability(int x) throws MathException { 238 // calculate the probability using half-correction 239 return normal.cumulativeProbability(x + 0.5); 240 } 241 242 /** 243 * Generates a random value sampled from this distribution. 244 * 245 * <p><strong>Algorithm Description</strong>: 246 * <ul><li> For small means, uses simulation of a Poisson process 247 * using Uniform deviates, as described 248 * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a> 249 * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li>< 250 * 251 * <li> For large means, uses the rejection algorithm described in <br/> 252 * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i> 253 * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p> 254 * 255 * @return random value 256 * @since 2.2 257 * @throws MathException if an error occurs generating the random value 258 */ 259 @Override 260 public int sample() throws MathException { 261 return (int) FastMath.min(randomData.nextPoisson(mean), Integer.MAX_VALUE); 262 } 263 264 /** 265 * Access the domain value lower bound, based on <code>p</code>, used to 266 * bracket a CDF root. This method is used by 267 * {@link #inverseCumulativeProbability(double)} to find critical values. 268 * 269 * @param p the desired probability for the critical value 270 * @return domain lower bound 271 */ 272 @Override 273 protected int getDomainLowerBound(double p) { 274 return 0; 275 } 276 277 /** 278 * Access the domain value upper bound, based on <code>p</code>, used to 279 * bracket a CDF root. This method is used by 280 * {@link #inverseCumulativeProbability(double)} to find critical values. 281 * 282 * @param p the desired probability for the critical value 283 * @return domain upper bound 284 */ 285 @Override 286 protected int getDomainUpperBound(double p) { 287 return Integer.MAX_VALUE; 288 } 289 290 /** 291 * Modify the normal distribution used to compute normal approximations. The 292 * caller is responsible for insuring the normal distribution has the proper 293 * parameter settings. 294 * 295 * @param value the new distribution 296 * @since 1.2 297 * @deprecated as of 2.1 (class will become immutable in 3.0) 298 */ 299 @Deprecated 300 public void setNormal(NormalDistribution value) { 301 setNormalAndMeanInternal(value, mean); 302 } 303 304 /** 305 * Returns the lower bound of the support for the distribution. 306 * 307 * The lower bound of the support is always 0 no matter the mean parameter. 308 * 309 * @return lower bound of the support (always 0) 310 * @since 2.2 311 */ 312 public int getSupportLowerBound() { 313 return 0; 314 } 315 316 /** 317 * Returns the upper bound of the support for the distribution. 318 * 319 * The upper bound of the support is positive infinity, 320 * regardless of the parameter values. There is no integer infinity, 321 * so this method returns <code>Integer.MAX_VALUE</code> and 322 * {@link #isSupportUpperBoundInclusive()} returns <code>true</code>. 323 * 324 * @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity) 325 * @since 2.2 326 */ 327 public int getSupportUpperBound() { 328 return Integer.MAX_VALUE; 329 } 330 331 /** 332 * Returns the variance of the distribution. 333 * 334 * For mean parameter <code>p</code>, the variance is <code>p</code> 335 * 336 * @return the variance 337 * @since 2.2 338 */ 339 public double getNumericalVariance() { 340 return getMean(); 341 } 342 343} 344