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.util.LocalizedFormats; 021 import org.apache.commons.math3.special.Gamma; 022 import org.apache.commons.math3.util.MathUtils; 023 import org.apache.commons.math3.util.ArithmeticUtils; 024 import org.apache.commons.math3.util.FastMath; 025 import org.apache.commons.math3.random.RandomGenerator; 026 import org.apache.commons.math3.random.Well19937c; 027 028 /** 029 * Implementation of the Poisson distribution. 030 * 031 * @see <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution (Wikipedia)</a> 032 * @see <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution (MathWorld)</a> 033 * @version $Id: PoissonDistribution.java 1416643 2012-12-03 19:37:14Z tn $ 034 */ 035 public class PoissonDistribution extends AbstractIntegerDistribution { 036 /** 037 * Default maximum number of iterations for cumulative probability calculations. 038 * @since 2.1 039 */ 040 public static final int DEFAULT_MAX_ITERATIONS = 10000000; 041 /** 042 * Default convergence criterion. 043 * @since 2.1 044 */ 045 public static final double DEFAULT_EPSILON = 1e-12; 046 /** Serializable version identifier. */ 047 private static final long serialVersionUID = -3349935121172596109L; 048 /** Distribution used to compute normal approximation. */ 049 private final NormalDistribution normal; 050 /** Distribution needed for the {@link #sample()} method. */ 051 private final ExponentialDistribution exponential; 052 /** Mean of the distribution. */ 053 private final double mean; 054 055 /** 056 * Maximum number of iterations for cumulative probability. Cumulative 057 * probabilities are estimated using either Lanczos series approximation 058 * of {@link Gamma#regularizedGammaP(double, double, double, int)} 059 * or continued fraction approximation of 060 * {@link Gamma#regularizedGammaQ(double, double, double, int)}. 061 */ 062 private final int maxIterations; 063 064 /** Convergence criterion for cumulative probability. */ 065 private final double epsilon; 066 067 /** 068 * Creates a new Poisson distribution with specified mean. 069 * 070 * @param p the Poisson mean 071 * @throws NotStrictlyPositiveException if {@code p <= 0}. 072 */ 073 public PoissonDistribution(double p) throws NotStrictlyPositiveException { 074 this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS); 075 } 076 077 /** 078 * Creates a new Poisson distribution with specified mean, convergence 079 * criterion and maximum number of iterations. 080 * 081 * @param p Poisson mean. 082 * @param epsilon Convergence criterion for cumulative probabilities. 083 * @param maxIterations the maximum number of iterations for cumulative 084 * probabilities. 085 * @throws NotStrictlyPositiveException if {@code p <= 0}. 086 * @since 2.1 087 */ 088 public PoissonDistribution(double p, double epsilon, int maxIterations) 089 throws NotStrictlyPositiveException { 090 this(new Well19937c(), p, epsilon, maxIterations); 091 } 092 093 /** 094 * Creates a new Poisson distribution with specified mean, convergence 095 * criterion and maximum number of iterations. 096 * 097 * @param rng Random number generator. 098 * @param p Poisson mean. 099 * @param epsilon Convergence criterion for cumulative probabilities. 100 * @param maxIterations the maximum number of iterations for cumulative 101 * probabilities. 102 * @throws NotStrictlyPositiveException if {@code p <= 0}. 103 * @since 3.1 104 */ 105 public PoissonDistribution(RandomGenerator rng, 106 double p, 107 double epsilon, 108 int maxIterations) 109 throws NotStrictlyPositiveException { 110 super(rng); 111 112 if (p <= 0) { 113 throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, p); 114 } 115 mean = p; 116 this.epsilon = epsilon; 117 this.maxIterations = maxIterations; 118 119 // Use the same RNG instance as the parent class. 120 normal = new NormalDistribution(rng, p, FastMath.sqrt(p), 121 NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 122 exponential = new ExponentialDistribution(rng, 1, 123 ExponentialDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 124 } 125 126 /** 127 * Creates a new Poisson distribution with the specified mean and 128 * convergence criterion. 129 * 130 * @param p Poisson mean. 131 * @param epsilon Convergence criterion for cumulative probabilities. 132 * @throws NotStrictlyPositiveException if {@code p <= 0}. 133 * @since 2.1 134 */ 135 public PoissonDistribution(double p, double epsilon) 136 throws NotStrictlyPositiveException { 137 this(p, epsilon, DEFAULT_MAX_ITERATIONS); 138 } 139 140 /** 141 * Creates a new Poisson distribution with the specified mean and maximum 142 * number of iterations. 143 * 144 * @param p Poisson mean. 145 * @param maxIterations Maximum number of iterations for cumulative 146 * probabilities. 147 * @since 2.1 148 */ 149 public PoissonDistribution(double p, int maxIterations) { 150 this(p, DEFAULT_EPSILON, maxIterations); 151 } 152 153 /** 154 * Get the mean for the distribution. 155 * 156 * @return the mean for the distribution. 157 */ 158 public double getMean() { 159 return mean; 160 } 161 162 /** {@inheritDoc} */ 163 public double probability(int x) { 164 double ret; 165 if (x < 0 || x == Integer.MAX_VALUE) { 166 ret = 0.0; 167 } else if (x == 0) { 168 ret = FastMath.exp(-mean); 169 } else { 170 ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) - 171 SaddlePointExpansion.getDeviancePart(x, mean)) / 172 FastMath.sqrt(MathUtils.TWO_PI * x); 173 } 174 return ret; 175 } 176 177 /** {@inheritDoc} */ 178 public double cumulativeProbability(int x) { 179 if (x < 0) { 180 return 0; 181 } 182 if (x == Integer.MAX_VALUE) { 183 return 1; 184 } 185 return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, 186 maxIterations); 187 } 188 189 /** 190 * Calculates the Poisson distribution function using a normal 191 * approximation. The {@code N(mean, sqrt(mean))} distribution is used 192 * to approximate the Poisson distribution. The computation uses 193 * "half-correction" (evaluating the normal distribution function at 194 * {@code x + 0.5}). 195 * 196 * @param x Upper bound, inclusive. 197 * @return the distribution function value calculated using a normal 198 * approximation. 199 */ 200 public double normalApproximateProbability(int x) { 201 // calculate the probability using half-correction 202 return normal.cumulativeProbability(x + 0.5); 203 } 204 205 /** 206 * {@inheritDoc} 207 * 208 * For mean parameter {@code p}, the mean is {@code p}. 209 */ 210 public double getNumericalMean() { 211 return getMean(); 212 } 213 214 /** 215 * {@inheritDoc} 216 * 217 * For mean parameter {@code p}, the variance is {@code p}. 218 */ 219 public double getNumericalVariance() { 220 return getMean(); 221 } 222 223 /** 224 * {@inheritDoc} 225 * 226 * The lower bound of the support is always 0 no matter the mean parameter. 227 * 228 * @return lower bound of the support (always 0) 229 */ 230 public int getSupportLowerBound() { 231 return 0; 232 } 233 234 /** 235 * {@inheritDoc} 236 * 237 * The upper bound of the support is positive infinity, 238 * regardless of the parameter values. There is no integer infinity, 239 * so this method returns {@code Integer.MAX_VALUE}. 240 * 241 * @return upper bound of the support (always {@code Integer.MAX_VALUE} for 242 * positive infinity) 243 */ 244 public int getSupportUpperBound() { 245 return Integer.MAX_VALUE; 246 } 247 248 /** 249 * {@inheritDoc} 250 * 251 * The support of this distribution is connected. 252 * 253 * @return {@code true} 254 */ 255 public boolean isSupportConnected() { 256 return true; 257 } 258 259 /** 260 * {@inheritDoc} 261 * <p> 262 * <strong>Algorithm Description</strong>: 263 * <ul> 264 * <li>For small means, uses simulation of a Poisson process 265 * using Uniform deviates, as described 266 * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here</a>. 267 * The Poisson process (and hence value returned) is bounded by 1000 * mean. 268 * </li> 269 * <li>For large means, uses the rejection algorithm described in 270 * <quote> 271 * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i> 272 * <strong>Computing</strong> vol. 26 pp. 197-207. 273 * </quote> 274 * </li> 275 * </ul> 276 * </p> 277 * 278 * @return a random value. 279 * @since 2.2 280 */ 281 @Override 282 public int sample() { 283 return (int) FastMath.min(nextPoisson(mean), Integer.MAX_VALUE); 284 } 285 286 /** 287 * @param meanPoisson Mean of the Poisson distribution. 288 * @return the next sample. 289 */ 290 private long nextPoisson(double meanPoisson) { 291 final double pivot = 40.0d; 292 if (meanPoisson < pivot) { 293 double p = FastMath.exp(-meanPoisson); 294 long n = 0; 295 double r = 1.0d; 296 double rnd = 1.0d; 297 298 while (n < 1000 * meanPoisson) { 299 rnd = random.nextDouble(); 300 r = r * rnd; 301 if (r >= p) { 302 n++; 303 } else { 304 return n; 305 } 306 } 307 return n; 308 } else { 309 final double lambda = FastMath.floor(meanPoisson); 310 final double lambdaFractional = meanPoisson - lambda; 311 final double logLambda = FastMath.log(lambda); 312 final double logLambdaFactorial = ArithmeticUtils.factorialLog((int) lambda); 313 final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional); 314 final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1)); 315 final double halfDelta = delta / 2; 316 final double twolpd = 2 * lambda + delta; 317 final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / 8 * lambda); 318 final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd); 319 final double aSum = a1 + a2 + 1; 320 final double p1 = a1 / aSum; 321 final double p2 = a2 / aSum; 322 final double c1 = 1 / (8 * lambda); 323 324 double x = 0; 325 double y = 0; 326 double v = 0; 327 int a = 0; 328 double t = 0; 329 double qr = 0; 330 double qa = 0; 331 for (;;) { 332 final double u = random.nextDouble(); 333 if (u <= p1) { 334 final double n = random.nextGaussian(); 335 x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d; 336 if (x > delta || x < -lambda) { 337 continue; 338 } 339 y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x); 340 final double e = exponential.sample(); 341 v = -e - (n * n / 2) + c1; 342 } else { 343 if (u > p1 + p2) { 344 y = lambda; 345 break; 346 } else { 347 x = delta + (twolpd / delta) * exponential.sample(); 348 y = FastMath.ceil(x); 349 v = -exponential.sample() - delta * (x + 1) / twolpd; 350 } 351 } 352 a = x < 0 ? 1 : 0; 353 t = y * (y + 1) / (2 * lambda); 354 if (v < -t && a == 0) { 355 y = lambda + y; 356 break; 357 } 358 qr = t * ((2 * y + 1) / (6 * lambda) - 1); 359 qa = qr - (t * t) / (3 * (lambda + a * (y + 1))); 360 if (v < qa) { 361 y = lambda + y; 362 break; 363 } 364 if (v > qr) { 365 continue; 366 } 367 if (v < y * logLambda - ArithmeticUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) { 368 y = lambda + y; 369 break; 370 } 371 } 372 return y2 + (long) y; 373 } 374 } 375 }