1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.distribution;
18
19 import java.io.Serializable;
20
21 import org.apache.commons.math.MathException;
22 import org.apache.commons.math.MathRuntimeException;
23 import org.apache.commons.math.special.Gamma;
24 import org.apache.commons.math.util.MathUtils;
25
26 /**
27 * Implementation for the {@link PoissonDistribution}.
28 *
29 * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
30 */
31 public class PoissonDistributionImpl extends AbstractIntegerDistribution
32 implements PoissonDistribution, Serializable {
33
34 /** Serializable version identifier */
35 private static final long serialVersionUID = -3349935121172596109L;
36
37 /** Distribution used to compute normal approximation. */
38 private NormalDistribution normal;
39
40 /**
41 * Holds the Poisson mean for the distribution.
42 */
43 private double mean;
44
45 /**
46 * Create a new Poisson distribution with the given the mean.
47 * The mean value must be positive; otherwise an
48 * <code>IllegalArgument</code> is thrown.
49 *
50 * @param p the Poisson mean
51 * @throws IllegalArgumentException if p ≤ 0
52 */
53 public PoissonDistributionImpl(double p) {
54 this(p, new NormalDistributionImpl());
55 }
56
57 /**
58 * Create a new Poisson distribution with the given the mean.
59 * The mean value must be positive; otherwise an
60 * <code>IllegalArgument</code> is thrown.
61 *
62 * @param p the Poisson mean
63 * @param z a normal distribution used to compute normal approximations.
64 * @throws IllegalArgumentException if p ≤ 0
65 * @since 1.2
66 */
67 public PoissonDistributionImpl(double p, NormalDistribution z) {
68 super();
69 setNormal(z);
70 setMean(p);
71 }
72
73 /**
74 * Get the Poisson mean for the distribution.
75 *
76 * @return the Poisson mean for the distribution.
77 */
78 public double getMean() {
79 return this.mean;
80 }
81
82 /**
83 * Set the Poisson mean for the distribution.
84 * The mean value must be positive; otherwise an
85 * <code>IllegalArgument</code> is thrown.
86 *
87 * @param p the Poisson mean value
88 * @throws IllegalArgumentException if p ≤ 0
89 */
90 public void setMean(double p) {
91 if (p <= 0) {
92 throw MathRuntimeException.createIllegalArgumentException(
93 "the Poisson mean must be positive ({0})",
94 p);
95 }
96 this.mean = p;
97 normal.setMean(p);
98 normal.setStandardDeviation(Math.sqrt(p));
99 }
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 }