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.exception.NotStrictlyPositiveException;
022 import org.apache.commons.math.exception.util.LocalizedFormats;
023 import org.apache.commons.math.special.Gamma;
024 import org.apache.commons.math.util.MathUtils;
025 import org.apache.commons.math.util.FastMath;
026
027 /**
028 * Implementation for the {@link PoissonDistribution}.
029 *
030 * @version $Id: PoissonDistributionImpl.java 1178295 2011-10-03 04:36:27Z psteitz $
031 */
032 public class PoissonDistributionImpl extends AbstractIntegerDistribution
033 implements PoissonDistribution, Serializable {
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 * Default convergence criterion.
041 * @since 2.1
042 */
043 public static final double DEFAULT_EPSILON = 1e-12;
044 /** Serializable version identifier. */
045 private static final long serialVersionUID = -3349935121172596109L;
046 /** Distribution used to compute normal approximation. */
047 private final NormalDistribution normal;
048 /** Mean of the distribution. */
049 private final double mean;
050 /**
051 * Maximum number of iterations for cumulative probability.
052 *
053 * Cumulative probabilities are estimated using either Lanczos series approximation of
054 * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ.
055 */
056 private final int maxIterations;
057 /**
058 * Convergence criterion for cumulative probability.
059 */
060 private final double epsilon;
061
062 /**
063 * Create a new Poisson distribution with the given the mean. The mean value
064 * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
065 *
066 * @param p the Poisson mean
067 * @throws NotStrictlyPositiveException if {@code p <= 0}.
068 */
069 public PoissonDistributionImpl(double p) {
070 this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
071 }
072
073 /**
074 * Create a new Poisson distribution with the given mean, convergence criterion
075 * and maximum number of iterations.
076 *
077 * @param p Poisson mean.
078 * @param epsilon Convergence criterion for cumulative probabilities.
079 * @param maxIterations the maximum number of iterations for cumulative
080 * probabilities.
081 * @since 2.1
082 */
083 public PoissonDistributionImpl(double p, double epsilon, int maxIterations) {
084 if (p <= 0) {
085 throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, p);
086 }
087 mean = p;
088 normal = new NormalDistributionImpl(p, FastMath.sqrt(p));
089 this.epsilon = epsilon;
090 this.maxIterations = maxIterations;
091 }
092
093 /**
094 * Create a new Poisson distribution with the given mean and convergence criterion.
095 *
096 * @param p Poisson mean.
097 * @param epsilon Convergence criterion for cumulative probabilities.
098 * @since 2.1
099 */
100 public PoissonDistributionImpl(double p, double epsilon) {
101 this(p, epsilon, DEFAULT_MAX_ITERATIONS);
102 }
103
104 /**
105 * Create a new Poisson distribution with the given mean and maximum number of iterations.
106 *
107 * @param p Poisson mean.
108 * @param maxIterations Maximum number of iterations for cumulative probabilities.
109 * @since 2.1
110 */
111 public PoissonDistributionImpl(double p, int maxIterations) {
112 this(p, DEFAULT_EPSILON, maxIterations);
113 }
114
115 /**
116 * {@inheritDoc}
117 */
118 public double getMean() {
119 return mean;
120 }
121
122 /**
123 * The probability mass function {@code P(X = x)} for a Poisson distribution.
124 *
125 * @param x Value at which the probability density function is evaluated.
126 * @return the value of the probability mass function at {@code x}.
127 */
128 public double probability(int x) {
129 double ret;
130 if (x < 0 || x == Integer.MAX_VALUE) {
131 ret = 0.0;
132 } else if (x == 0) {
133 ret = FastMath.exp(-mean);
134 } else {
135 ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) -
136 SaddlePointExpansion.getDeviancePart(x, mean)) /
137 FastMath.sqrt(MathUtils.TWO_PI * x);
138 }
139 return ret;
140 }
141
142 /**
143 * The probability distribution function {@code P(X <= x)} for a Poisson
144 * distribution.
145 *
146 * @param x Value at which the PDF is evaluated.
147 * @return the Poisson distribution function evaluated at {@code x}.
148 * due to convergence or other numerical errors.
149 */
150 @Override
151 public double cumulativeProbability(int x) {
152 if (x < 0) {
153 return 0;
154 }
155 if (x == Integer.MAX_VALUE) {
156 return 1;
157 }
158 return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
159 }
160
161 /**
162 * Calculates the Poisson distribution function using a normal
163 * approximation. The {@code N(mean, sqrt(mean))} distribution is used
164 * to approximate the Poisson distribution.
165 * The computation uses "half-correction" (evaluating the normal
166 * distribution function at {@code x + 0.5}).
167 *
168 * @param x Upper bound, inclusive.
169 * @return the distribution function value calculated using a normal
170 * approximation.
171 * approximation.
172 */
173 public double normalApproximateProbability(int x) {
174 // calculate the probability using half-correction
175 return normal.cumulativeProbability(x + 0.5);
176 }
177
178 /**
179 * Generates a random value sampled from this distribution.
180 * <br/>
181 * <strong>Algorithm Description</strong>:
182 * <ul>
183 * <li>For small means, uses simulation of a Poisson process
184 * using Uniform deviates, as described
185 * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here</a>.
186 * The Poisson process (and hence value returned) is bounded by 1000 * mean.
187 * </li>
188 * <li>For large means, uses the rejection algorithm described in
189 * <quote>
190 * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
191 * <strong>Computing</strong> vol. 26 pp. 197-207.
192 * </quote>
193 * </li>
194 * </ul>
195 *
196 * @return a random value.
197 * @since 2.2
198 */
199 @Override
200 public int sample() {
201 return (int) FastMath.min(randomData.nextPoisson(mean), Integer.MAX_VALUE);
202 }
203
204 /**
205 * Access the domain value lower bound, based on {@code p}, used to
206 * bracket a CDF root. This method is used by
207 * {@link #inverseCumulativeProbability(double)} to find critical values.
208 *
209 * @param p Desired probability for the critical value.
210 * @return the domain lower bound.
211 */
212 @Override
213 protected int getDomainLowerBound(double p) {
214 return 0;
215 }
216
217 /**
218 * Access the domain value upper bound, based on {@code p}, used to
219 * bracket a CDF root. This method is used by
220 * {@link #inverseCumulativeProbability(double)} to find critical values.
221 *
222 * @param p Desired probability for the critical value.
223 * @return the domain upper bound.
224 */
225 @Override
226 protected int getDomainUpperBound(double p) {
227 return Integer.MAX_VALUE;
228 }
229
230 /**
231 * {@inheritDoc}
232 *
233 * The lower bound of the support is always 0 no matter the mean parameter.
234 *
235 * @return lower bound of the support (always 0)
236 */
237 @Override
238 public int getSupportLowerBound() {
239 return 0;
240 }
241
242 /**
243 * {@inheritDoc}
244 *
245 * The upper bound of the support is positive infinity,
246 * regardless of the parameter values. There is no integer infinity,
247 * so this method returns <code>Integer.MAX_VALUE</code> and
248 * {@link #isSupportUpperBoundInclusive()} returns <code>true</code>.
249 *
250 * @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity)
251 */
252 @Override
253 public int getSupportUpperBound() {
254 return Integer.MAX_VALUE;
255 }
256
257 /**
258 * {@inheritDoc}
259 *
260 * For mean parameter <code>p</code>, the mean is <code>p</code>
261 *
262 * @return {@inheritDoc}
263 */
264 @Override
265 protected double calculateNumericalMean() {
266 return getMean();
267 }
268
269 /**
270 * {@inheritDoc}
271 *
272 * For mean parameter <code>p</code>, the variance is <code>p</code>
273 *
274 * @return {@inheritDoc}
275 */
276 @Override
277 protected double calculateNumericalVariance() {
278 return getMean();
279 }
280
281 /**
282 * {@inheritDoc}
283 */
284 @Override
285 public boolean isSupportUpperBoundInclusive() {
286 return true;
287 }
288 }