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 ≤ 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 ≤ 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 ≤ 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 ≤ 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 }