/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.special.Gamma;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.FastMath;
/**
* Implementation for the {@link PoissonDistribution}.
*
* @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
*/
public class PoissonDistributionImpl extends AbstractIntegerDistribution
implements PoissonDistribution, Serializable {
/**
* Default maximum number of iterations for cumulative probability calculations.
* @since 2.1
*/
public static final int DEFAULT_MAX_ITERATIONS = 10000000;
/**
* Default convergence criterion.
* @since 2.1
*/
public static final double DEFAULT_EPSILON = 1E-12;
/** Serializable version identifier */
private static final long serialVersionUID = -3349935121172596109L;
/** Distribution used to compute normal approximation. */
private NormalDistribution normal;
/**
* Holds the Poisson mean for the distribution.
*/
private double mean;
/**
* Maximum number of iterations for cumulative probability.
*
* Cumulative probabilities are estimated using either Lanczos series approximation of
* Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ.
*/
private int maxIterations = DEFAULT_MAX_ITERATIONS;
/**
* Convergence criterion for cumulative probability.
*/
private double epsilon = DEFAULT_EPSILON;
/**
* Create a new Poisson distribution with the given the mean. The mean value
* must be positive; otherwise an <code>IllegalArgument</code> is thrown.
*
* @param p the Poisson mean
* @throws IllegalArgumentException if p ≤ 0
*/
public PoissonDistributionImpl(double p) {
this(p, new NormalDistributionImpl());
}
/**
* Create a new Poisson distribution with the given mean, convergence criterion
* and maximum number of iterations.
*
* @param p the Poisson mean
* @param epsilon the convergence criteria for cumulative probabilites
* @param maxIterations the maximum number of iterations for cumulative probabilites
* @since 2.1
*/
public PoissonDistributionImpl(double p, double epsilon, int maxIterations) {
setMean(p);
this.epsilon = epsilon;
this.maxIterations = maxIterations;
}
/**
* Create a new Poisson distribution with the given mean and convergence criterion.
*
* @param p the Poisson mean
* @param epsilon the convergence criteria for cumulative probabilites
* @since 2.1
*/
public PoissonDistributionImpl(double p, double epsilon) {
setMean(p);
this.epsilon = epsilon;
}
/**
* Create a new Poisson distribution with the given mean and maximum number of iterations.
*
* @param p the Poisson mean
* @param maxIterations the maximum number of iterations for cumulative probabilites
* @since 2.1
*/
public PoissonDistributionImpl(double p, int maxIterations) {
setMean(p);
this.maxIterations = maxIterations;
}
/**
* Create a new Poisson distribution with the given the mean. The mean value
* must be positive; otherwise an <code>IllegalArgument</code> is thrown.
*
* @param p the Poisson mean
* @param z a normal distribution used to compute normal approximations.
* @throws IllegalArgumentException if p ≤ 0
* @since 1.2
* @deprecated as of 2.1 (to avoid possibly inconsistent state, the
* "NormalDistribution" will be instantiated internally)
*/
@Deprecated
public PoissonDistributionImpl(double p, NormalDistribution z) {
super();
setNormalAndMeanInternal(z, p);
}
/**
* Get the Poisson mean for the distribution.
*
* @return the Poisson mean for the distribution.
*/
public double getMean() {
return mean;
}
/**
* Set the Poisson mean for the distribution. The mean value must be
* positive; otherwise an <code>IllegalArgument</code> is thrown.
*
* @param p the Poisson mean value
* @throws IllegalArgumentException if p ≤ 0
* @deprecated as of 2.1 (class will become immutable in 3.0)
*/
@Deprecated
public void setMean(double p) {
setNormalAndMeanInternal(normal, p);
}
/**
* Set the Poisson mean for the distribution. The mean value must be
* positive; otherwise an <code>IllegalArgument</code> is thrown.
*
* @param z the new distribution
* @param p the Poisson mean value
* @throws IllegalArgumentException if p ≤ 0
*/
private void setNormalAndMeanInternal(NormalDistribution z,
double p) {
if (p <= 0) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.NOT_POSITIVE_POISSON_MEAN, p);
}
mean = p;
normal = z;
normal.setMean(p);
normal.setStandardDeviation(FastMath.sqrt(p));
}
/**
* The probability mass function P(X = x) for a Poisson distribution.
*
* @param x the value at which the probability density function is
* evaluated.
* @return the value of the probability mass function at x
*/
public double probability(int x) {
double ret;
if (x < 0 || x == Integer.MAX_VALUE) {
ret = 0.0;
} else if (x == 0) {
ret = FastMath.exp(-mean);
} else {
ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) -
SaddlePointExpansion.getDeviancePart(x, mean)) /
FastMath.sqrt(MathUtils.TWO_PI * x);
}
return ret;
}
/**
* The probability distribution function P(X <= x) for a Poisson
* distribution.
*
* @param x the value at which the PDF is evaluated.
* @return Poisson distribution function evaluated at x
* @throws MathException if the cumulative probability can not be computed
* due to convergence or other numerical errors.
*/
@Override
public double cumulativeProbability(int x) throws MathException {
if (x < 0) {
return 0;
}
if (x == Integer.MAX_VALUE) {
return 1;
}
return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
}
/**
* Calculates the Poisson distribution function using a normal
* approximation. The <code>N(mean, sqrt(mean))</code> distribution is used
* to approximate the Poisson distribution.
* <p>
* The computation uses "half-correction" -- evaluating the normal
* distribution function at <code>x + 0.5</code>
* </p>
*
* @param x the upper bound, inclusive
* @return the distribution function value calculated using a normal
* approximation
* @throws MathException if an error occurs computing the normal
* approximation
*/
public double normalApproximateProbability(int x) throws MathException {
// calculate the probability using half-correction
return normal.cumulativeProbability(x + 0.5);
}
/**
* Generates a random value sampled from this distribution.
*
* <p><strong>Algorithm Description</strong>:
* <ul><li> For small means, uses simulation of a Poisson process
* using Uniform deviates, as described
* <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
* The Poisson process (and hence value returned) is bounded by 1000 * mean.</li><
*
* <li> For large means, uses the rejection algorithm described in <br/>
* Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
* <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p>
*
* @return random value
* @since 2.2
* @throws MathException if an error occurs generating the random value
*/
@Override
public int sample() throws MathException {
return (int) FastMath.min(randomData.nextPoisson(mean), Integer.MAX_VALUE);
}
/**
* Access the domain value lower bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
* @return domain lower bound
*/
@Override
protected int getDomainLowerBound(double p) {
return 0;
}
/**
* Access the domain value upper bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* {@link #inverseCumulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
* @return domain upper bound
*/
@Override
protected int getDomainUpperBound(double p) {
return Integer.MAX_VALUE;
}
/**
* Modify the normal distribution used to compute normal approximations. The
* caller is responsible for insuring the normal distribution has the proper
* parameter settings.
*
* @param value the new distribution
* @since 1.2
* @deprecated as of 2.1 (class will become immutable in 3.0)
*/
@Deprecated
public void setNormal(NormalDistribution value) {
setNormalAndMeanInternal(value, mean);
}
/**
* Returns the lower bound of the support for the distribution.
*
* The lower bound of the support is always 0 no matter the mean parameter.
*
* @return lower bound of the support (always 0)
* @since 2.2
*/
public int getSupportLowerBound() {
return 0;
}
/**
* Returns the upper bound of the support for the distribution.
*
* The upper bound of the support is positive infinity,
* regardless of the parameter values. There is no integer infinity,
* so this method returns <code>Integer.MAX_VALUE</code> and
* {@link #isSupportUpperBoundInclusive()} returns <code>true</code>.
*
* @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity)
* @since 2.2
*/
public int getSupportUpperBound() {
return Integer.MAX_VALUE;
}
/**
* Returns the variance of the distribution.
*
* For mean parameter <code>p</code>, the variance is <code>p</code>
*
* @return the variance
* @since 2.2
*/
public double getNumericalVariance() {
return getMean();
}
}