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: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
030     */
031    public class PoissonDistributionImpl extends AbstractIntegerDistribution
032            implements PoissonDistribution, Serializable {
033    
034        /** Serializable version identifier */
035        private static final long serialVersionUID = -3349935121172596109L;
036    
037        /** Distribution used to compute normal approximation. */
038        private NormalDistribution normal;
039        
040        /**
041         * Holds the Poisson mean for the distribution.
042         */
043        private double mean;
044    
045        /**
046         * Create a new Poisson distribution with the given the mean.
047         * The mean value must be positive; otherwise an 
048         * <code>IllegalArgument</code> is thrown.
049         * 
050         * @param p the Poisson mean
051         * @throws IllegalArgumentException if p &le; 0
052         */
053        public PoissonDistributionImpl(double p) {
054            this(p, new NormalDistributionImpl());
055        }
056    
057        /**
058         * Create a new Poisson distribution with the given the mean.
059         * The mean value must be positive; otherwise an 
060         * <code>IllegalArgument</code> is thrown.
061         * 
062         * @param p the Poisson mean
063         * @param z a normal distribution used to compute normal approximations.
064         * @throws IllegalArgumentException if p &le; 0
065         * @since 1.2
066         */
067        public PoissonDistributionImpl(double p, NormalDistribution z) {
068            super();
069            setNormal(z);
070            setMean(p);
071        }
072    
073        /**
074         * Get the Poisson mean for the distribution.
075         * 
076         * @return the Poisson mean for the distribution.
077         */
078        public double getMean() {
079            return this.mean;
080        }
081    
082        /**
083         * Set the Poisson mean for the distribution.
084         * The mean value must be positive; otherwise an 
085         * <code>IllegalArgument</code> is thrown.
086         * 
087         * @param p the Poisson mean value
088         * @throws IllegalArgumentException if p &le; 0
089         */
090        public void setMean(double p) {
091            if (p <= 0) {
092                throw MathRuntimeException.createIllegalArgumentException(
093                      "the Poisson mean must be positive ({0})",
094                      p);
095            }
096            this.mean = p;
097            normal.setMean(p);
098            normal.setStandardDeviation(Math.sqrt(p));
099        }
100    
101        /**
102         * The probability mass function P(X = x) for a Poisson distribution.
103         * 
104         * @param x the value at which the probability density function is evaluated.
105         * @return the value of the probability mass function at x
106         */
107        public double probability(int x) {
108            if (x < 0 || x == Integer.MAX_VALUE) {
109                return 0;
110            }
111            return Math.pow(getMean(), x) / 
112                MathUtils.factorialDouble(x) * Math.exp(-mean);
113        }
114        
115        /**
116         * The probability distribution function P(X <= x) for a Poisson distribution.
117         * 
118         * @param x the value at which the PDF is evaluated.
119         * @return Poisson distribution function evaluated at x
120         * @throws MathException if the cumulative probability can not be
121         *            computed due to convergence or other numerical errors.
122         */
123        @Override
124        public double cumulativeProbability(int x) throws MathException {
125            if (x < 0) {
126                return 0;
127            }
128            if (x == Integer.MAX_VALUE) {
129                return 1;
130            }
131            return Gamma.regularizedGammaQ((double)x + 1, mean, 
132                    1E-12, Integer.MAX_VALUE);
133        }
134    
135        /**
136         * Calculates the Poisson distribution function using a normal
137         * approximation.  The <code>N(mean, sqrt(mean))</code>
138         * distribution is used to approximate the Poisson distribution.
139         * <p>
140         * The computation uses "half-correction" -- evaluating the normal
141         * distribution function at <code>x + 0.5</code></p>
142         * 
143         * @param x the upper bound, inclusive
144         * @return the distribution function value calculated using a normal approximation
145         * @throws MathException if an error occurs computing the normal approximation
146         */
147        public double normalApproximateProbability(int x) throws MathException {
148            // calculate the probability using half-correction
149            return normal.cumulativeProbability(x + 0.5);
150        }
151    
152        /**
153         * Access the domain value lower bound, based on <code>p</code>, used to
154         * bracket a CDF root.  This method is used by
155         * {@link #inverseCumulativeProbability(double)} to find critical values.
156         * 
157         * @param p the desired probability for the critical value
158         * @return domain lower bound
159         */
160        @Override
161        protected int getDomainLowerBound(double p) {
162            return 0;
163        }
164    
165        /**
166         * Access the domain value upper bound, based on <code>p</code>, used to
167         * bracket a CDF root.  This method is used by
168         * {@link #inverseCumulativeProbability(double)} to find critical values.
169         * 
170         * @param p the desired probability for the critical value
171         * @return domain upper bound
172         */
173        @Override
174        protected int getDomainUpperBound(double p) {
175            return Integer.MAX_VALUE;
176        }
177        
178        /**
179         * Modify the normal distribution used to compute normal approximations.
180         * The caller is responsible for insuring the normal distribution has the
181         * proper parameter settings.
182         * @param value the new distribution
183         * @since 1.2
184         */
185        public void setNormal(NormalDistribution value) {
186            normal = value;
187        }
188        
189    }