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.special.Gamma; 025 import org.apache.commons.math.util.MathUtils; 026 import org.apache.commons.math.util.FastMath; 027 028 /** 029 * Implementation for the {@link PoissonDistribution}. 030 * 031 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 032 */ 033 public class PoissonDistributionImpl extends AbstractIntegerDistribution 034 implements PoissonDistribution, Serializable { 035 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 /** 043 * Default convergence criterion. 044 * @since 2.1 045 */ 046 public static final double DEFAULT_EPSILON = 1E-12; 047 048 /** Serializable version identifier */ 049 private static final long serialVersionUID = -3349935121172596109L; 050 051 /** Distribution used to compute normal approximation. */ 052 private NormalDistribution normal; 053 054 /** 055 * Holds the Poisson mean for the distribution. 056 */ 057 private double mean; 058 059 /** 060 * Maximum number of iterations for cumulative probability. 061 * 062 * Cumulative probabilities are estimated using either Lanczos series approximation of 063 * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ. 064 */ 065 private int maxIterations = DEFAULT_MAX_ITERATIONS; 066 067 /** 068 * Convergence criterion for cumulative probability. 069 */ 070 private double epsilon = DEFAULT_EPSILON; 071 072 /** 073 * Create a new Poisson distribution with the given the mean. The mean value 074 * must be positive; otherwise an <code>IllegalArgument</code> is thrown. 075 * 076 * @param p the Poisson mean 077 * @throws IllegalArgumentException if p ≤ 0 078 */ 079 public PoissonDistributionImpl(double p) { 080 this(p, new NormalDistributionImpl()); 081 } 082 083 /** 084 * Create a new Poisson distribution with the given mean, convergence criterion 085 * and maximum number of iterations. 086 * 087 * @param p the Poisson mean 088 * @param epsilon the convergence criteria for cumulative probabilites 089 * @param maxIterations the maximum number of iterations for cumulative probabilites 090 * @since 2.1 091 */ 092 public PoissonDistributionImpl(double p, double epsilon, int maxIterations) { 093 setMean(p); 094 this.epsilon = epsilon; 095 this.maxIterations = maxIterations; 096 } 097 098 /** 099 * 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 }