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.special.Gamma; 024 import org.apache.commons.math.util.MathUtils; 025 026 /** 027 * Implementation for the {@link PoissonDistribution}. 028 * 029 * @version $Revision: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $ 030 */ 031 public class PoissonDistributionImpl extends AbstractIntegerDistribution 032 implements PoissonDistribution, Serializable { 033 034 /** 035 * Default maximum number of iterations for cumulative probability calculations. 036 * @since 2.1 037 */ 038 public static final int DEFAULT_MAX_ITERATIONS = 10000000; 039 040 /** 041 * Default convergence criterion. 042 * @since 2.1 043 */ 044 public static final double DEFAULT_EPSILON = 1E-12; 045 046 /** Serializable version identifier */ 047 private static final long serialVersionUID = -3349935121172596109L; 048 049 /** Distribution used to compute normal approximation. */ 050 private NormalDistribution normal; 051 052 /** 053 * Holds the Poisson mean for the distribution. 054 */ 055 private double mean; 056 057 /** 058 * Maximum number of iterations for cumulative probability. 059 * 060 * Cumulative probabilities are estimated using either Lanczos series approximation of 061 * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ. 062 */ 063 private int maxIterations = DEFAULT_MAX_ITERATIONS; 064 065 /** 066 * Convergence criterion for cumulative probability. 067 */ 068 private double epsilon = DEFAULT_EPSILON; 069 070 /** 071 * Create a new Poisson distribution with the given the mean. The mean value 072 * must be positive; otherwise an <code>IllegalArgument</code> is thrown. 073 * 074 * @param p the Poisson mean 075 * @throws IllegalArgumentException if p ≤ 0 076 */ 077 public PoissonDistributionImpl(double p) { 078 this(p, new NormalDistributionImpl()); 079 } 080 081 /** 082 * Create a new Poisson distribution with the given mean, convergence criterion 083 * and maximum number of iterations. 084 * 085 * @param p the Poisson mean 086 * @param epsilon the convergence criteria for cumulative probabilites 087 * @param maxIterations the maximum number of iterations for cumulative probabilites 088 * @since 2.1 089 */ 090 public PoissonDistributionImpl(double p, double epsilon, int maxIterations) { 091 setMean(p); 092 this.epsilon = epsilon; 093 this.maxIterations = maxIterations; 094 } 095 096 /** 097 * Create a new Poisson distribution with the given mean and convergence criterion. 098 * 099 * @param p the Poisson mean 100 * @param epsilon the convergence criteria for cumulative probabilites 101 * @since 2.1 102 */ 103 public PoissonDistributionImpl(double p, double epsilon) { 104 setMean(p); 105 this.epsilon = epsilon; 106 } 107 108 /** 109 * Create a new Poisson distribution with the given mean and maximum number of iterations. 110 * 111 * @param p the Poisson mean 112 * @param maxIterations the maximum number of iterations for cumulative probabilites 113 * @since 2.1 114 */ 115 public PoissonDistributionImpl(double p, int maxIterations) { 116 setMean(p); 117 this.maxIterations = maxIterations; 118 } 119 120 121 /** 122 * Create a new Poisson distribution with the given the mean. The mean value 123 * must be positive; otherwise an <code>IllegalArgument</code> is thrown. 124 * 125 * @param p the Poisson mean 126 * @param z a normal distribution used to compute normal approximations. 127 * @throws IllegalArgumentException if p ≤ 0 128 * @since 1.2 129 * @deprecated as of 2.1 (to avoid possibly inconsistent state, the 130 * "NormalDistribution" will be instantiated internally) 131 */ 132 @Deprecated 133 public PoissonDistributionImpl(double p, NormalDistribution z) { 134 super(); 135 setNormalAndMeanInternal(z, p); 136 } 137 138 /** 139 * Get the Poisson mean for the distribution. 140 * 141 * @return the Poisson mean for the distribution. 142 */ 143 public double getMean() { 144 return mean; 145 } 146 147 /** 148 * Set the Poisson mean for the distribution. The mean value must be 149 * positive; otherwise an <code>IllegalArgument</code> is thrown. 150 * 151 * @param p the Poisson mean value 152 * @throws IllegalArgumentException if p ≤ 0 153 * @deprecated as of 2.1 (class will become immutable in 3.0) 154 */ 155 @Deprecated 156 public void setMean(double p) { 157 setNormalAndMeanInternal(normal, p); 158 } 159 /** 160 * Set the Poisson mean for the distribution. The mean value must be 161 * positive; otherwise an <code>IllegalArgument</code> is thrown. 162 * 163 * @param z the new distribution 164 * @param p the Poisson mean value 165 * @throws IllegalArgumentException if p ≤ 0 166 */ 167 private void setNormalAndMeanInternal(NormalDistribution z, 168 double p) { 169 if (p <= 0) { 170 throw MathRuntimeException.createIllegalArgumentException( 171 "the Poisson mean must be positive ({0})", p); 172 } 173 mean = p; 174 normal = z; 175 normal.setMean(p); 176 normal.setStandardDeviation(Math.sqrt(p)); 177 } 178 179 /** 180 * The probability mass function P(X = x) for a Poisson distribution. 181 * 182 * @param x the value at which the probability density function is 183 * evaluated. 184 * @return the value of the probability mass function at x 185 */ 186 public double probability(int x) { 187 double ret; 188 if (x < 0 || x == Integer.MAX_VALUE) { 189 ret = 0.0; 190 } else if (x == 0) { 191 ret = Math.exp(-mean); 192 } else { 193 ret = Math.exp(-SaddlePointExpansion.getStirlingError(x) - 194 SaddlePointExpansion.getDeviancePart(x, mean)) / 195 Math.sqrt(MathUtils.TWO_PI * x); 196 } 197 return ret; 198 } 199 200 /** 201 * The probability distribution function P(X <= x) for a Poisson 202 * distribution. 203 * 204 * @param x the value at which the PDF is evaluated. 205 * @return Poisson distribution function evaluated at x 206 * @throws MathException if the cumulative probability can not be computed 207 * due to convergence or other numerical errors. 208 */ 209 @Override 210 public double cumulativeProbability(int x) throws MathException { 211 if (x < 0) { 212 return 0; 213 } 214 if (x == Integer.MAX_VALUE) { 215 return 1; 216 } 217 return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations); 218 } 219 220 /** 221 * Calculates the Poisson distribution function using a normal 222 * approximation. The <code>N(mean, sqrt(mean))</code> distribution is used 223 * to approximate the Poisson distribution. 224 * <p> 225 * The computation uses "half-correction" -- evaluating the normal 226 * distribution function at <code>x + 0.5</code> 227 * </p> 228 * 229 * @param x the upper bound, inclusive 230 * @return the distribution function value calculated using a normal 231 * approximation 232 * @throws MathException if an error occurs computing the normal 233 * approximation 234 */ 235 public double normalApproximateProbability(int x) throws MathException { 236 // calculate the probability using half-correction 237 return normal.cumulativeProbability(x + 0.5); 238 } 239 240 /** 241 * Access the domain value lower bound, based on <code>p</code>, used to 242 * bracket a CDF root. This method is used by 243 * {@link #inverseCumulativeProbability(double)} to find critical values. 244 * 245 * @param p the desired probability for the critical value 246 * @return domain lower bound 247 */ 248 @Override 249 protected int getDomainLowerBound(double p) { 250 return 0; 251 } 252 253 /** 254 * Access the domain value upper bound, based on <code>p</code>, used to 255 * bracket a CDF root. This method is used by 256 * {@link #inverseCumulativeProbability(double)} to find critical values. 257 * 258 * @param p the desired probability for the critical value 259 * @return domain upper bound 260 */ 261 @Override 262 protected int getDomainUpperBound(double p) { 263 return Integer.MAX_VALUE; 264 } 265 266 /** 267 * Modify the normal distribution used to compute normal approximations. The 268 * caller is responsible for insuring the normal distribution has the proper 269 * parameter settings. 270 * 271 * @param value the new distribution 272 * @since 1.2 273 * @deprecated as of 2.1 (class will become immutable in 3.0) 274 */ 275 @Deprecated 276 public void setNormal(NormalDistribution value) { 277 setNormalAndMeanInternal(value, mean); 278 } 279 }