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.special.Gamma;
024    import org.apache.commons.math.util.MathUtils;
025    
026    /**
027     * Implementation for the {@link PoissonDistribution}.
028     *
029     * @version $Revision: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $
030     */
031    public class PoissonDistributionImpl extends AbstractIntegerDistribution
032            implements PoissonDistribution, Serializable {
033    
034        /**
035         * Default maximum number of iterations for cumulative probability calculations.
036         * @since 2.1
037         */
038        public static final int DEFAULT_MAX_ITERATIONS = 10000000;
039    
040        /**
041         * Default convergence criterion.
042         * @since 2.1
043         */
044        public static final double DEFAULT_EPSILON = 1E-12;
045    
046        /** Serializable version identifier */
047        private static final long serialVersionUID = -3349935121172596109L;
048    
049        /** Distribution used to compute normal approximation. */
050        private NormalDistribution normal;
051    
052        /**
053         * Holds the Poisson mean for the distribution.
054         */
055        private double mean;
056    
057        /**
058         * Maximum number of iterations for cumulative probability.
059         *
060         * Cumulative probabilities are estimated using either Lanczos series approximation of
061         * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ.
062         */
063        private int maxIterations = DEFAULT_MAX_ITERATIONS;
064    
065        /**
066         * Convergence criterion for cumulative probability.
067         */
068        private double epsilon = DEFAULT_EPSILON;
069    
070        /**
071         * Create a new Poisson distribution with the given the mean. The mean value
072         * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
073         *
074         * @param p the Poisson mean
075         * @throws IllegalArgumentException if p &le; 0
076         */
077        public PoissonDistributionImpl(double p) {
078            this(p, new NormalDistributionImpl());
079        }
080    
081        /**
082         * Create a new Poisson distribution with the given mean, convergence criterion
083         * and maximum number of iterations.
084         *
085         * @param p the Poisson mean
086         * @param epsilon the convergence criteria for cumulative probabilites
087         * @param maxIterations the maximum number of iterations for cumulative probabilites
088         * @since 2.1
089         */
090        public PoissonDistributionImpl(double p, double epsilon, int maxIterations) {
091            setMean(p);
092            this.epsilon = epsilon;
093            this.maxIterations = maxIterations;
094        }
095    
096        /**
097         * Create a new Poisson distribution with the given mean and convergence criterion.
098         *
099         * @param p the Poisson mean
100         * @param epsilon the convergence criteria for cumulative probabilites
101         * @since 2.1
102         */
103        public PoissonDistributionImpl(double p, double epsilon) {
104            setMean(p);
105            this.epsilon = epsilon;
106        }
107    
108        /**
109         * Create a new Poisson distribution with the given mean and maximum number of iterations.
110         *
111         * @param p the Poisson mean
112         * @param maxIterations the maximum number of iterations for cumulative probabilites
113         * @since 2.1
114         */
115        public PoissonDistributionImpl(double p, int maxIterations) {
116            setMean(p);
117            this.maxIterations = maxIterations;
118        }
119    
120    
121        /**
122         * Create a new Poisson distribution with the given the mean. The mean value
123         * must be positive; otherwise an <code>IllegalArgument</code> is thrown.
124         *
125         * @param p the Poisson mean
126         * @param z a normal distribution used to compute normal approximations.
127         * @throws IllegalArgumentException if p &le; 0
128         * @since 1.2
129         * @deprecated as of 2.1 (to avoid possibly inconsistent state, the
130         * "NormalDistribution" will be instantiated internally)
131         */
132        @Deprecated
133        public PoissonDistributionImpl(double p, NormalDistribution z) {
134            super();
135            setNormalAndMeanInternal(z, p);
136        }
137    
138        /**
139         * Get the Poisson mean for the distribution.
140         *
141         * @return the Poisson mean for the distribution.
142         */
143        public double getMean() {
144            return mean;
145        }
146    
147        /**
148         * Set the Poisson mean for the distribution. The mean value must be
149         * positive; otherwise an <code>IllegalArgument</code> is thrown.
150         *
151         * @param p the Poisson mean value
152         * @throws IllegalArgumentException if p &le; 0
153         * @deprecated as of 2.1 (class will become immutable in 3.0)
154         */
155        @Deprecated
156        public void setMean(double p) {
157            setNormalAndMeanInternal(normal, p);
158        }
159        /**
160         * Set the Poisson mean for the distribution. The mean value must be
161         * positive; otherwise an <code>IllegalArgument</code> is thrown.
162         *
163         * @param z the new distribution
164         * @param p the Poisson mean value
165         * @throws IllegalArgumentException if p &le; 0
166         */
167        private void setNormalAndMeanInternal(NormalDistribution z,
168                                              double p) {
169            if (p <= 0) {
170                throw MathRuntimeException.createIllegalArgumentException(
171                        "the Poisson mean must be positive ({0})", p);
172            }
173            mean = p;
174            normal = z;
175            normal.setMean(p);
176            normal.setStandardDeviation(Math.sqrt(p));
177        }
178    
179        /**
180         * The probability mass function P(X = x) for a Poisson distribution.
181         *
182         * @param x the value at which the probability density function is
183         *            evaluated.
184         * @return the value of the probability mass function at x
185         */
186        public double probability(int x) {
187            double ret;
188            if (x < 0 || x == Integer.MAX_VALUE) {
189                ret = 0.0;
190            } else if (x == 0) {
191                ret = Math.exp(-mean);
192            } else {
193                ret = Math.exp(-SaddlePointExpansion.getStirlingError(x) -
194                      SaddlePointExpansion.getDeviancePart(x, mean)) /
195                      Math.sqrt(MathUtils.TWO_PI * x);
196            }
197            return ret;
198        }
199    
200        /**
201         * The probability distribution function P(X <= x) for a Poisson
202         * distribution.
203         *
204         * @param x the value at which the PDF is evaluated.
205         * @return Poisson distribution function evaluated at x
206         * @throws MathException if the cumulative probability can not be computed
207         *             due to convergence or other numerical errors.
208         */
209        @Override
210        public double cumulativeProbability(int x) throws MathException {
211            if (x < 0) {
212                return 0;
213            }
214            if (x == Integer.MAX_VALUE) {
215                return 1;
216            }
217            return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
218        }
219    
220        /**
221         * Calculates the Poisson distribution function using a normal
222         * approximation. The <code>N(mean, sqrt(mean))</code> distribution is used
223         * to approximate the Poisson distribution.
224         * <p>
225         * The computation uses "half-correction" -- evaluating the normal
226         * distribution function at <code>x + 0.5</code>
227         * </p>
228         *
229         * @param x the upper bound, inclusive
230         * @return the distribution function value calculated using a normal
231         *         approximation
232         * @throws MathException if an error occurs computing the normal
233         *             approximation
234         */
235        public double normalApproximateProbability(int x) throws MathException {
236            // calculate the probability using half-correction
237            return normal.cumulativeProbability(x + 0.5);
238        }
239    
240        /**
241         * Access the domain value lower bound, based on <code>p</code>, used to
242         * bracket a CDF root. This method is used by
243         * {@link #inverseCumulativeProbability(double)} to find critical values.
244         *
245         * @param p the desired probability for the critical value
246         * @return domain lower bound
247         */
248        @Override
249        protected int getDomainLowerBound(double p) {
250            return 0;
251        }
252    
253        /**
254         * Access the domain value upper bound, based on <code>p</code>, used to
255         * bracket a CDF root. This method is used by
256         * {@link #inverseCumulativeProbability(double)} to find critical values.
257         *
258         * @param p the desired probability for the critical value
259         * @return domain upper bound
260         */
261        @Override
262        protected int getDomainUpperBound(double p) {
263            return Integer.MAX_VALUE;
264        }
265    
266        /**
267         * Modify the normal distribution used to compute normal approximations. The
268         * caller is responsible for insuring the normal distribution has the proper
269         * parameter settings.
270         *
271         * @param value the new distribution
272         * @since 1.2
273         * @deprecated as of 2.1 (class will become immutable in 3.0)
274         */
275        @Deprecated
276        public void setNormal(NormalDistribution value) {
277            setNormalAndMeanInternal(value, mean);
278        }
279    }