001 package com.softnetConsult.utils.random;
002
003 import java.util.ArrayList;
004 import java.util.List;
005 import java.util.Random;
006
007 import com.softnetConsult.utils.math.MathTools;
008
009 /**
010 * Creates Poisson-distributed pseudorandom numbers.
011 *
012 * <p style="font-size:smaller;">This product includes software developed by the
013 * <strong>SoftNet-Consult Java Utility Library</strong> project and its contributors.<br />
014 * (<a href="http://java-tools.sourceforge.net" target="_blank">http://java-tools.sourceforge.net</a>)<br />
015 * Copyright (c) 2007-2008 SoftNet-Consult.<br />
016 * Copyright (c) 2007-2008 G. Paperin.<br />
017 * All rights reserved.
018 * </p>
019 * <p style="font-size:smaller;">File: Poisson.java<br />
020 * Library API version: {@value com.softnetConsult.utils.APIProperties#apiVersion}<br />
021 * Java compliance version: {@value com.softnetConsult.utils.APIProperties#javaComplianceVersion}
022 * </p>
023 * <p style="font-size:smaller;">Redistribution and use in source and binary forms, with or
024 * without modification, are permitted provided that the following terms and conditions are met:
025 * </p>
026 * <p style="font-size:smaller;">1. Redistributions of source code must retain the above
027 * acknowledgement of the SoftNet-Consult Java Utility Library project, the above copyright
028 * notice, this list of conditions and the following disclaimer.<br />
029 * 2. Redistributions in binary form must reproduce the above acknowledgement of the
030 * SoftNet-Consult Java Utility Library project, the above copyright notice, this list of
031 * conditions and the following disclaimer in the documentation and/or other materials
032 * provided with the distribution.<br />
033 * 3. All advertising materials mentioning features or use of this software or any derived
034 * software must display the following acknowledgement:<br />
035 * <em>This product includes software developed by the SoftNet-Consult Java Utility Library
036 * project and its contributors.</em>
037 * </p>
038 * <p style="font-size:smaller;">THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY
039 * OF ANY KIND, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
040 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
041 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
042 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
043 * IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
044 * </p>
045 * @author Greg Paperin (<a href="http://www.paperin.org" target="_blank">http://www.paperin.org</a>)
046 * @version {@value com.softnetConsult.utils.APIProperties#apiVersion}
047 */
048 public class Poisson {
049
050
051 /**
052 * Default value for epsilon when no other value is specified.
053 */
054 public static final double DEFAULT_EPSILON = Math.pow(10, -10);
055
056
057 /**
058 * Specifies the completeness of the CDF table - the table is calculated until the
059 * approximated CDF value reaches {@code 1.0 - epsilon}.
060 */
061 private double epsilon = DEFAULT_EPSILON;
062
063
064 /**
065 * The lambda parameter of this Poisson distribution.
066 * Both, the mean and the variance of a Poisson-distributed random variable are equal to lambda.
067 */
068 private double lambda = 0;
069
070
071 /**
072 * The CDF table.
073 */
074 private double [] table = null;
075
076
077 /**
078 * The source of uniformly independently distributed random numbers.
079 */
080 private Random rndSource = null;
081
082
083 /**
084 * Creates a generator for Poisson distributed random numbers with the specified lambda-parameter.
085 * Calling this constructor is equivalent to {@code Poisson(lambda, new Random(), Poisson.DEFAULT_EPSILON)}.
086 *
087 * @param lambda The LAMBDA-parameter to the Poisson distribution.
088 * @see #Poisson(double, Random, double)
089 */
090 public Poisson(double lambda) {
091 this(lambda, new Random(), DEFAULT_EPSILON);
092 }
093
094 /**
095 * Creates a generator for Poisson distributed random numbers with the specified lambda-parameter
096 * and the specified source of random numbers.
097 *
098 * Calling this constructor is equivalent to {@code Poisson(lambda, rndSource, Poisson.DEFAULT_EPSILON)}.
099 *
100 * @param lambda The LAMBDA-parameter to the Poisson distribution.
101 * @param rndSource The source of uniformly independently distributed (pseudo-) random numbers to use for
102 * this distribution.
103 * @see #Poisson(double, Random, double)
104 */
105 public Poisson(double lambda, Random rndSource) {
106 this(lambda, rndSource, DEFAULT_EPSILON);
107 }
108
109
110 /**
111 * Creates a generator for Poisson distributed random numbers.
112 *
113 * @param lambda The LAMBDA parameter for this Poisson distribution.
114 * Both, the mean and the variance of a Poisson-distributed random variable are equal to {@code lambda}.
115 * @param rndSource The source of uniformly independently distributed (pseudo-) random numbers to use for
116 * this distribution.
117 * @param epsilon Completeness threshold used during the computation of the cumulative distribution function
118 * table. The CDF table is calculated until the approximated CDF value reaches {@code 1.0 - epsilon}.
119 */
120 public Poisson(double lambda, Random rndSource, double epsilon) {
121
122 if (1 > lambda)
123 throw new IllegalArgumentException("Poisson parameter lambda must be >= 1.0");
124
125 if (null == rndSource)
126 throw new NullPointerException("Cannot use a null random source");
127
128 if (epsilon < 0.0)
129 throw new IllegalArgumentException("Epsilon must be a very small positive number");
130
131 if (epsilon > 0.05)
132 throw new IllegalArgumentException("The specified value of epsilon is too large");
133
134 this.lambda = lambda;
135 this.rndSource = rndSource;
136 this.epsilon = epsilon;
137
138 List<Double> tab = new ArrayList<Double>();
139 double cdf = 0, pmf = 0;
140 int k = 0;
141 do {
142 pmf = MathTools.poisson(lambda, k++);
143 if (Double.isInfinite(pmf) || Double.isNaN(pmf)) {
144 throw new IllegalArgumentException("Poisson parameter lambda (=" + lambda + ") is too large,"
145 + " or epsilon (=" + epsilon + ") is too small:"
146 + " the cumulative distribution table cannot be computed with"
147 + " sufficient accuracy");
148 }
149 cdf += pmf;
150 tab.add(cdf);
151 } while (cdf < 1. - epsilon);
152
153 this.table = new double[tab.size()];
154 for (k = 0; k < this.table.length; k++)
155 this.table[k] = tab.get(k);
156
157 //printTable();
158 }
159
160 /*
161 private void printTable() {
162 // For debug only:
163 for (int x = 0; x < table.length; x++) {
164 System.out.println("PoissonCDF(lambda="+lambda+", x="+x+") = " + table[x]);
165 }
166 }
167 */
168
169 /**
170 * Returns expected mean and variance of this Poisson distribution.
171 *
172 * @return The LAMBDA parameter of this Poisson distribution.
173 */
174 public double getLambda() {
175 return lambda;
176 }
177
178 /**
179 * The epsilon value used during the calculation of the CDF table.
180 *
181 * @return Completeness threshold used during the computation of the cumulative distribution function
182 * table. The CDF table is calculated until the approximated CDF value reaches {@code 1.0 - epsilon}.
183 */
184 public double getEpsilon() {
185 return epsilon;
186 }
187
188 /**
189 * Generates a new (pseudo-) random integer number from this Poisson-distribution.
190 *
191 * @return A new Poisson-distributed integer.
192 */
193 public int nextRandom() {
194 double dice = rndSource.nextDouble();
195 return tableSearch(dice);
196 }
197
198
199 /**
200 * Used internally to look up the value of the inverse of the cumulative distribution function
201 * of the Poisson distribution.
202 *
203 * @param dice A number in interval [0, 1]
204 * @return The smallest integer <em>k</em> such that <em>{@code dice} <= CDF(k)</em>,
205 * where <em>CDF</em> is the cumulative distribution function of this Poisson distribution.
206 */
207 private int tableSearch(double dice) {
208 int left = 0;
209 int right = table.length - 1;
210
211 while (left <= right) {
212 int pivot = (left + right) >>> 1;
213 double pivotVal = table[pivot];
214
215 if (pivotVal < dice)
216 left = pivot + 1;
217 else if (pivotVal > dice)
218 right = pivot - 1;
219 else
220 return pivot;
221 }
222 return left;
223 }
224
225 } // public class Poisson