Java tutorial
/* * 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(); } }