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 * Create an instance. 58 * 59 * @param rng Generator of uniformly distributed random numbers. 60 * @param mean Mean. 61 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0} 62 */ 63 public SmallMeanPoissonSampler(UniformRandomProvider rng, 64 double mean) { 65 this(rng, mean, computeP0(mean)); 66 } 67 68 /** 69 * Instantiates a new small mean poisson sampler. 70 * 71 * @param rng Generator of uniformly distributed random numbers. 72 * @param mean Mean. 73 * @param p0 {@code Math.exp(-mean)}. 74 */ 75 private SmallMeanPoissonSampler(UniformRandomProvider rng, 76 double mean, 77 double p0) { 78 this.rng = rng; 79 this.p0 = p0; 80 // The returned sample is bounded by 1000 * mean 81 limit = (int) Math.ceil(1000 * mean); 82 } 83 84 /** 85 * @param rng Generator of uniformly distributed random numbers. 86 * @param source Source to copy. 87 */ 88 private SmallMeanPoissonSampler(UniformRandomProvider rng, 89 SmallMeanPoissonSampler source) { 90 this.rng = rng; 91 p0 = source.p0; 92 limit = source.limit; 93 } 94 95 /** {@inheritDoc} */ 96 @Override 97 public int sample() { 98 int n = 0; 99 double r = 1; 100 101 while (n < limit) { 102 r *= rng.nextDouble(); 103 if (r >= p0) { 104 n++; 105 } else { 106 break; 107 } 108 } 109 return n; 110 } 111 112 /** {@inheritDoc} */ 113 @Override 114 public String toString() { 115 return "Small Mean Poisson deviate [" + rng.toString() + "]"; 116 } 117 118 /** 119 * {@inheritDoc} 120 * 121 * @since 1.3 122 */ 123 @Override 124 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) { 125 return new SmallMeanPoissonSampler(rng, this); 126 } 127 128 /** 129 * Creates a new sampler for the Poisson distribution. 130 * 131 * @param rng Generator of uniformly distributed random numbers. 132 * @param mean Mean of the distribution. 133 * @return the sampler 134 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0}. 135 * @since 1.3 136 */ 137 public static SharedStateDiscreteSampler of(UniformRandomProvider rng, 138 double mean) { 139 return new SmallMeanPoissonSampler(rng, mean); 140 } 141 142 /** 143 * Compute {@code Math.exp(-mean)}. 144 * 145 * <p>This method exists to raise an exception before invocation of the 146 * private constructor; this mitigates Finalizer attacks 147 * (see SpotBugs CT_CONSTRUCTOR_THROW). 148 * 149 * @param mean Mean. 150 * @return the mean 151 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0} 152 */ 153 private static double computeP0(double mean) { 154 InternalUtils.requireStrictlyPositive(mean, "mean"); 155 final double p0 = Math.exp(-mean); 156 if (p0 > 0) { 157 return p0; 158 } 159 // This excludes NaN values for the mean 160 throw new IllegalArgumentException("No p(x=0) probability for mean: " + mean); 161 } 162 }