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