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 }