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 27 * described <a href="http://mathaa.epfl.ch/cours/PMMI2001/interactive/rng7.htm">here</a>. 28 * The Poisson process (and hence, the returned value) is bounded by 1000 * mean. 29 * </li> 30 * </ul> 31 * 32 * @since 1.1 33 * 34 * This sampler is suitable for {@code mean < 40}. 35 * For large means, {@link LargeMeanPoissonSampler} should be used instead. 36 */ 37 public class SmallMeanPoissonSampler 38 implements DiscreteSampler { 39 /** Upper bound to avoid truncation. */ 40 private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE; 41 /** 42 * Pre-compute {@code Math.exp(-mean)}. 43 * Note: This is the probability of the Poisson sample {@code P(n=0)}. 44 */ 45 private final double p0; 46 /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */ 47 private final int limit; 48 /** Underlying source of randomness. */ 49 private final UniformRandomProvider rng; 50 51 /** 52 * @param rng Generator of uniformly distributed random numbers. 53 * @param mean Mean. 54 * @throws IllegalArgumentException if {@code mean <= 0}. 55 */ 56 public SmallMeanPoissonSampler(UniformRandomProvider rng, 57 double mean) { 58 this.rng = rng; 59 if (mean <= 0) { 60 throw new IllegalArgumentException(mean + " <= " + 0); 61 } 62 if (mean > MAX_MEAN) { 63 throw new IllegalArgumentException(mean + " > " + MAX_MEAN); 64 } 65 66 p0 = Math.exp(-mean); 67 // The returned sample is bounded by 1000 * mean or Integer.MAX_VALUE 68 limit = (int) Math.ceil(Math.min(1000 * mean, Integer.MAX_VALUE)); 69 } 70 71 /** {@inheritDoc} */ 72 @Override 73 public int sample() { 74 int n = 0; 75 double r = 1; 76 77 while (n < limit) { 78 r *= rng.nextDouble(); 79 if (r >= p0) { 80 n++; 81 } else { 82 break; 83 } 84 } 85 return n; 86 } 87 88 /** {@inheritDoc} */ 89 @Override 90 public String toString() { 91 return "Small Mean Poisson deviate [" + rng.toString() + "]"; 92 } 93 }