    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.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.special.Gamma;
    025    import org.apache.commons.math.util.MathUtils;
    026    import org.apache.commons.math.util.FastMath;
    027
    028    /**
    029     * Implementation for the {@link PoissonDistribution}.
    030     *
    031     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
    032     */
    033    public class PoissonDistributionImpl extends AbstractIntegerDistribution
    034            implements PoissonDistribution, Serializable {
    035
    036        /**
    037         * Default maximum number of iterations for cumulative probability calculations.
    038         * @since 2.1
    039         */
    040        public static final int DEFAULT_MAX_ITERATIONS = 10000000;
    041
    042        /**
    043         * Default convergence criterion.
    044         * @since 2.1
    045         */
    046        public static final double DEFAULT_EPSILON = 1E-12;
    047
    048        /** Serializable version identifier */
    049        private static final long serialVersionUID = -3349935121172596109L;
    050
    051        /** Distribution used to compute normal approximation. */
    052        private NormalDistribution normal;
    053
    054        /**
    055         * Holds the Poisson mean for the distribution.
    056         */
    057        private double mean;
    058
    059        /**
    060         * Maximum number of iterations for cumulative probability.
    061         *
    062         * Cumulative probabilities are estimated using either Lanczos series approximation of
    063         * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ.
    064         */
    065        private int maxIterations = DEFAULT_MAX_ITERATIONS;
    066
    067        /**
    068         * Convergence criterion for cumulative probability.
    069         */
    070        private double epsilon = DEFAULT_EPSILON;
    071
    072        /**
    073         * Create a new Poisson distribution with the given the mean. The mean value
    074         * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
    075         *
    076         * @param p the Poisson mean
    077         * @throws IllegalArgumentException if p &le; 0
    078         */
    079        public PoissonDistributionImpl(double p) {
    080            this(p, new NormalDistributionImpl());
    081        }
    082
    083        /**
    084         * Create a new Poisson distribution with the given mean, convergence criterion
    085         * and maximum number of iterations.
    086         *
    087         * @param p the Poisson mean
    088         * @param epsilon the convergence criteria for cumulative probabilites
    089         * @param maxIterations the maximum number of iterations for cumulative probabilites
    090         * @since 2.1
    091         */
    092        public PoissonDistributionImpl(double p, double epsilon, int maxIterations) {
    093            setMean(p);
    094            this.epsilon = epsilon;
    095            this.maxIterations = maxIterations;
    096        }
    097
    098        /**
    099         * Create a new Poisson distribution with the given mean and convergence criterion.
    100         *
    101         * @param p the Poisson mean
    102         * @param epsilon the convergence criteria for cumulative probabilites
    103         * @since 2.1
    104         */
    105        public PoissonDistributionImpl(double p, double epsilon) {
    106            setMean(p);
    107            this.epsilon = epsilon;
    108        }
    109
    110        /**
    111         * Create a new Poisson distribution with the given mean and maximum number of iterations.
    112         *
    113         * @param p the Poisson mean
    114         * @param maxIterations the maximum number of iterations for cumulative probabilites
    115         * @since 2.1
    116         */
    117        public PoissonDistributionImpl(double p, int maxIterations) {
    118            setMean(p);
    119            this.maxIterations = maxIterations;
    120        }
    121
    122
    123        /**
    124         * Create a new Poisson distribution with the given the mean. The mean value
    125         * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
    126         *
    127         * @param p the Poisson mean
    128         * @param z a normal distribution used to compute normal approximations.
    129         * @throws IllegalArgumentException if p &le; 0
    130         * @since 1.2
    131         * @deprecated as of 2.1 (to avoid possibly inconsistent state, the
    132         * "NormalDistribution" will be instantiated internally)
    133         */
    134        @Deprecated
    135        public PoissonDistributionImpl(double p, NormalDistribution z) {
    136            super();
    137            setNormalAndMeanInternal(z, p);
    138        }
    139
    140        /**
    141         * Get the Poisson mean for the distribution.
    142         *
    143         * @return the Poisson mean for the distribution.
    144         */
    145        public double getMean() {
    146            return mean;
    147        }
    148
    149        /**
    150         * Set the Poisson mean for the distribution. The mean value must be
    151         * positive; otherwise an <code>IllegalArgument</code> is thrown.
    152         *
    153         * @param p the Poisson mean value
    154         * @throws IllegalArgumentException if p &le; 0
    155         * @deprecated as of 2.1 (class will become immutable in 3.0)
    156         */
    157        @Deprecated
    158        public void setMean(double p) {
    159            setNormalAndMeanInternal(normal, p);
    160        }
    161        /**
    162         * Set the Poisson mean for the distribution. The mean value must be
    163         * positive; otherwise an <code>IllegalArgument</code> is thrown.
    164         *
    165         * @param z the new distribution
    166         * @param p the Poisson mean value
    167         * @throws IllegalArgumentException if p &le; 0
    168         */
    169        private void setNormalAndMeanInternal(NormalDistribution z,
    170                                              double p) {
    171            if (p <= 0) {
    172                throw MathRuntimeException.createIllegalArgumentException(
    173                        LocalizedFormats.NOT_POSITIVE_POISSON_MEAN, p);
    174            }
    175            mean = p;
    176            normal = z;
    177            normal.setMean(p);
    178            normal.setStandardDeviation(FastMath.sqrt(p));
    179        }
    180
    181        /**
    182         * The probability mass function P(X = x) for a Poisson distribution.
    183         *
    184         * @param x the value at which the probability density function is
    185         *            evaluated.
    186         * @return the value of the probability mass function at x
    187         */
    188        public double probability(int x) {
    189            double ret;
    190            if (x < 0 || x == Integer.MAX_VALUE) {
    191                ret = 0.0;
    192            } else if (x == 0) {
    193                ret = FastMath.exp(-mean);
    194            } else {
    195                ret = FastMath.exp(-SaddlePointExpansion.getStirlingError(x) -
    196                      SaddlePointExpansion.getDeviancePart(x, mean)) /
    197                      FastMath.sqrt(MathUtils.TWO_PI * x);
    198            }
    199            return ret;
    200        }
    201
    202        /**
    203         * The probability distribution function P(X <= x) for a Poisson
    204         * distribution.
    205         *
    206         * @param x the value at which the PDF is evaluated.
    207         * @return Poisson distribution function evaluated at x
    208         * @throws MathException if the cumulative probability can not be computed
    209         *             due to convergence or other numerical errors.
    210         */
    211        @Override
    212        public double cumulativeProbability(int x) throws MathException {
    213            if (x < 0) {
    214                return 0;
    215            }
    216            if (x == Integer.MAX_VALUE) {
    217                return 1;
    218            }
    219            return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
    220        }
    221
    222        /**
    223         * Calculates the Poisson distribution function using a normal
    224         * approximation. The <code>N(mean, sqrt(mean))</code> distribution is used
    225         * to approximate the Poisson distribution.
    226         * <p>
    227         * The computation uses "half-correction" -- evaluating the normal
    228         * distribution function at <code>x + 0.5</code>
    229         * </p>
    230         *
    231         * @param x the upper bound, inclusive
    232         * @return the distribution function value calculated using a normal
    233         *         approximation
    234         * @throws MathException if an error occurs computing the normal
    235         *             approximation
    236         */
    237        public double normalApproximateProbability(int x) throws MathException {
    238            // calculate the probability using half-correction
    239            return normal.cumulativeProbability(x + 0.5);
    240        }
    241
    242        /**
    243         * Generates a random value sampled from this distribution.
    244         *
    245         * <p><strong>Algorithm Description</strong>:
    246         * <ul><li> For small means, uses simulation of a Poisson process
    247         * using Uniform deviates, as described
    248         * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
    249         * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li><
    250         *
    251         * <li> For large means, uses the rejection algorithm described in <br/>
    252         * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
    253         * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p>
    254         *
    255         * @return random value
    256         * @since 2.2
    257         * @throws MathException if an error occurs generating the random value
    258         */
    259        @Override
    260        public int sample() throws MathException {
    261            return (int) FastMath.min(randomData.nextPoisson(mean), Integer.MAX_VALUE);
    262        }
    263
    264        /**
    265         * Access the domain value lower bound, based on <code>p</code>, used to
    266         * bracket a CDF root. This method is used by
    267         * {@link #inverseCumulativeProbability(double)} to find critical values.
    268         *
    269         * @param p the desired probability for the critical value
    270         * @return domain lower bound
    271         */
    272        @Override
    273        protected int getDomainLowerBound(double p) {
    274            return 0;
    275        }
    276
    277        /**
    278         * Access the domain value upper bound, based on <code>p</code>, used to
    279         * bracket a CDF root. This method is used by
    280         * {@link #inverseCumulativeProbability(double)} to find critical values.
    281         *
    282         * @param p the desired probability for the critical value
    283         * @return domain upper bound
    284         */
    285        @Override
    286        protected int getDomainUpperBound(double p) {
    287            return Integer.MAX_VALUE;
    288        }
    289
    290        /**
    291         * Modify the normal distribution used to compute normal approximations. The
    292         * caller is responsible for insuring the normal distribution has the proper
    293         * parameter settings.
    294         *
    295         * @param value the new distribution
    296         * @since 1.2
    297         * @deprecated as of 2.1 (class will become immutable in 3.0)
    298         */
    299        @Deprecated
    300        public void setNormal(NormalDistribution value) {
    301            setNormalAndMeanInternal(value, mean);
    302        }
    303
    304        /**
    305         * Returns the lower bound of the support for the distribution.
    306         *
    307         * The lower bound of the support is always 0 no matter the mean parameter.
    308         *
    309         * @return lower bound of the support (always 0)
    310         * @since 2.2
    311         */
    312        public int getSupportLowerBound() {
    313            return 0;
    314        }
    315
    316        /**
    317         * Returns the upper bound of the support for the distribution.
    318         *
    319         * The upper bound of the support is positive infinity,
    320         * regardless of the parameter values. There is no integer infinity,
    321         * so this method returns <code>Integer.MAX_VALUE</code> and
    322         * {@link #isSupportUpperBoundInclusive()} returns <code>true</code>.
    323         *
    324         * @return upper bound of the support (always <code>Integer.MAX_VALUE</code> for positive infinity)
    325         * @since 2.2
    326         */
    327        public int getSupportUpperBound() {
    328            return Integer.MAX_VALUE;
    329        }
    330
    331        /**
    332         * Returns the variance of the distribution.
    333         *
    334         * For mean parameter <code>p</code>, the variance is <code>p</code>
    335         *
    336         * @return the variance
    337         * @since 2.2
    338         */
    339        public double getNumericalVariance() {
    340            return getMean();
    341        }
    342
    343    }