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 distribution</a>. 23 * 24 * <ul> 25 * <li> 26 * For small means, a Poisson process is simulated using uniform deviates, as described in 27 * <blockquote> 28 * Knuth (1969). <i>Seminumerical Algorithms</i>. The Art of Computer Programming, 29 * Volume 2. Chapter 3.4.1.F.3 Important integer-valued distributions: The Poisson distribution. 30 * Addison Wesley. 31 * </blockquote> 32 * The Poisson process (and hence, the returned value) is bounded by {@code 1000 * mean}. 33 * </li> 34 * </ul> 35 * 36 * <p>This sampler is suitable for {@code mean < 40}. 37 * For large means, {@link LargeMeanPoissonSampler} should be used instead.</p> 38 * 39 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()} and requires on average 40 * {@code mean + 1} deviates per sample.</p> 41 * 42 * @since 1.1 43 */ 44 public class SmallMeanPoissonSampler 45 implements SharedStateDiscreteSampler { 46 /** 47 * Pre-compute {@code Math.exp(-mean)}. 48 * Note: This is the probability of the Poisson sample {@code P(n=0)}. 49 */ 50 private final double p0; 51 /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */ 52 private final int limit; 53 /** Underlying source of randomness. */ 54 private final UniformRandomProvider rng; 55 56 /** 57 * @param rng Generator of uniformly distributed random numbers. 58 * @param mean Mean. 59 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0} 60 */ 61 public SmallMeanPoissonSampler(UniformRandomProvider rng, 62 double mean) { 63 this.rng = rng; 64 if (mean <= 0) { 65 throw new IllegalArgumentException("mean is not strictly positive: " + mean); 66 } 67 p0 = Math.exp(-mean); 68 if (p0 > 0) { 69 // The returned sample is bounded by 1000 * mean 70 limit = (int) Math.ceil(1000 * mean); 71 } else { 72 // This excludes NaN values for the mean 73 throw new IllegalArgumentException("No p(x=0) probability for mean: " + mean); 74 } 75 } 76 77 /** 78 * @param rng Generator of uniformly distributed random numbers. 79 * @param source Source to copy. 80 */ 81 private SmallMeanPoissonSampler(UniformRandomProvider rng, 82 SmallMeanPoissonSampler source) { 83 this.rng = rng; 84 p0 = source.p0; 85 limit = source.limit; 86 } 87 88 /** {@inheritDoc} */ 89 @Override 90 public int sample() { 91 int n = 0; 92 double r = 1; 93 94 while (n < limit) { 95 r *= rng.nextDouble(); 96 if (r >= p0) { 97 n++; 98 } else { 99 break; 100 } 101 } 102 return n; 103 } 104 105 /** {@inheritDoc} */ 106 @Override 107 public String toString() { 108 return "Small Mean Poisson deviate [" + rng.toString() + "]"; 109 } 110 111 /** 112 * {@inheritDoc} 113 * 114 * @since 1.3 115 */ 116 @Override 117 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 118 return new SmallMeanPoissonSampler(rng, this); 119 } 120 121 /** 122 * Creates a new sampler for the Poisson distribution. 123 * 124 * @param rng Generator of uniformly distributed random numbers. 125 * @param mean Mean of the distribution. 126 * @return the sampler 127 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0}. 128 * @since 1.3 129 */ 130 public static SharedStateDiscreteSampler of(UniformRandomProvider rng, 131 double mean) { 132 return new SmallMeanPoissonSampler(rng, mean); 133 } 134 }