All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math.distribution.PoissonDistributionImpl Maven / Gradle / Ivy

There is a newer version: 2024.11.18751.20241128T090041Z-241100
Show newest version
/*
 * 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 IllegalArgument 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 IllegalArgument 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 IllegalArgument 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 IllegalArgument 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 N(mean, sqrt(mean)) distribution is used
     * to approximate the Poisson distribution.
     * 

* The computation uses "half-correction" -- evaluating the normal * distribution function at x + 0.5 *

* * @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. * *

Algorithm Description: *

  • For small means, uses simulation of a Poisson process * using Uniform deviates, as described * here. * The Poisson process (and hence value returned) is bounded by 1000 * mean.
  • < * *
  • For large means, uses the rejection algorithm described in
    * Devroye, Luc. (1981).The Computer Generation of Poisson Random Variables * Computing vol. 26 pp. 197-207.

* * @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 p, 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 p, 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 Integer.MAX_VALUE and * {@link #isSupportUpperBoundInclusive()} returns true. * * @return upper bound of the support (always Integer.MAX_VALUE for positive infinity) * @since 2.2 */ public int getSupportUpperBound() { return Integer.MAX_VALUE; } /** * Returns the variance of the distribution. * * For mean parameter p, the variance is p * * @return the variance * @since 2.2 */ public double getNumericalVariance() { return getMean(); } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy