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.rng.sampling.distribution; 18 19 import org.apache.commons.rng.UniformRandomProvider; 20 21 /** 22 * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson 23 * distribution</a>. 24 * 25 * <ul> 26 * <li> 27 * Kemp, A, W, (1981) Efficient Generation of Logarithmically Distributed 28 * Pseudo-Random Variables. Journal of the Royal Statistical Society. Vol. 30, No. 3, pp. 29 * 249-253. 30 * </li> 31 * </ul> 32 * 33 * <p>This sampler is suitable for {@code mean < 40}. For large means, 34 * {@link LargeMeanPoissonSampler} should be used instead.</p> 35 * 36 * <p>Note: The algorithm uses a recurrence relation to compute the Poisson probability 37 * and a rolling summation for the cumulative probability. When the mean is large the 38 * initial probability (Math.exp(-mean)) is zero and an exception is raised by the 39 * constructor.</p> 40 * 41 * <p>Sampling uses 1 call to {@link UniformRandomProvider#nextDouble()}. This method provides 42 * an alternative to the {@link SmallMeanPoissonSampler} for slow generators of {@code double}.</p> 43 * 44 * @see <a href="https://www.jstor.org/stable/2346348">Kemp, A.W. (1981) JRSS Vol. 30, pp. 45 * 249-253</a> 46 * @since 1.3 47 */ 48 public final class KempSmallMeanPoissonSampler 49 implements SharedStateDiscreteSampler { 50 /** Underlying source of randomness. */ 51 private final UniformRandomProvider rng; 52 /** 53 * Pre-compute {@code Math.exp(-mean)}. 54 * Note: This is the probability of the Poisson sample {@code p(x=0)}. 55 */ 56 private final double p0; 57 /** 58 * The mean of the Poisson sample. 59 */ 60 private final double mean; 61 62 /** 63 * @param rng Generator of uniformly distributed random numbers. 64 * @param p0 Probability of the Poisson sample {@code p(x=0)}. 65 * @param mean Mean. 66 */ 67 private KempSmallMeanPoissonSampler(UniformRandomProvider rng, 68 double p0, 69 double mean) { 70 this.rng = rng; 71 this.p0 = p0; 72 this.mean = mean; 73 } 74 75 /** {@inheritDoc} */ 76 @Override 77 public int sample() { 78 // Note on the algorithm: 79 // - X is the unknown sample deviate (the output of the algorithm) 80 // - x is the current value from the distribution 81 // - p is the probability of the current value x, p(X=x) 82 // - u is effectively the cumulative probability that the sample X 83 // is equal or above the current value x, p(X>=x) 84 // So if p(X>=x) > p(X=x) the sample must be above x, otherwise it is x 85 double u = rng.nextDouble(); 86 int x = 0; 87 double p = p0; 88 while (u > p) { 89 u -= p; 90 // Compute the next probability using a recurrence relation. 91 // p(x+1) = p(x) * mean / (x+1) 92 p *= mean / ++x; 93 // The algorithm listed in Kemp (1981) does not check that the rolling probability 94 // is positive. This check is added to ensure no errors when the limit of the summation 95 // 1 - sum(p(x)) is above 0 due to cumulative error in floating point arithmetic. 96 if (p == 0) { 97 return x; 98 } 99 } 100 return x; 101 } 102 103 /** {@inheritDoc} */ 104 @Override 105 public String toString() { 106 return "Kemp Small Mean Poisson deviate [" + rng.toString() + "]"; 107 } 108 109 /** {@inheritDoc} */ 110 @Override 111 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 112 return new KempSmallMeanPoissonSampler(rng, p0, mean); 113 } 114 115 /** 116 * Creates a new sampler for the Poisson distribution. 117 * 118 * @param rng Generator of uniformly distributed random numbers. 119 * @param mean Mean of the distribution. 120 * @return the sampler 121 * @throws IllegalArgumentException if {@code mean <= 0} or 122 * {@code Math.exp(-mean) == 0}. 123 */ 124 public static SharedStateDiscreteSampler of(UniformRandomProvider rng, 125 double mean) { 126 if (mean <= 0) { 127 throw new IllegalArgumentException("Mean is not strictly positive: " + mean); 128 } 129 130 final double p0 = Math.exp(-mean); 131 132 // Probability must be positive. As mean increases then p(0) decreases. 133 if (p0 > 0) { 134 return new KempSmallMeanPoissonSampler(rng, p0, mean); 135 } 136 137 // This catches the edge case of a NaN mean 138 throw new IllegalArgumentException("No probability for mean: " + mean); 139 } 140 }