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 */ 36 public class SmallMeanPoissonSampler 37 implements DiscreteSampler { 38 39 /** 40 * Pre-compute {@code Math.exp(-mean)}. 41 * Note: This is the probability of the Poisson sample {@code P(n=0)}. 42 */ 43 private final double p0; 44 /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */ 45 private final int limit; 46 /** Underlying source of randomness. */ 47 private final UniformRandomProvider rng; 48 49 /** 50 * @param rng Generator of uniformly distributed random numbers. 51 * @param mean Mean. 52 * @throws IllegalArgumentException if {@code mean <= 0}. 53 */ 54 public SmallMeanPoissonSampler(UniformRandomProvider rng, 55 double mean) { 56 this.rng = rng; 57 if (mean <= 0) { 58 throw new IllegalArgumentException(mean + " <= " + 0); 59 } 60 61 p0 = Math.exp(-mean); 62 // The returned sample is bounded by 1000 * mean or Integer.MAX_VALUE 63 limit = (int) Math.ceil(Math.min(1000 * mean, Integer.MAX_VALUE)); 64 } 65 66 /** {@inheritDoc} */ 67 @Override 68 public int sample() { 69 int n = 0; 70 double r = 1; 71 72 while (n < limit) { 73 r *= rng.nextDouble(); 74 if (r >= p0) { 75 n++; 76 } else { 77 break; 78 } 79 } 80 return n; 81 } 82 83 /** {@inheritDoc} */ 84 @Override 85 public String toString() { 86 return "Small Mean Poisson deviate [" + rng.toString() + "]"; 87 } 88 }