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    }