001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.rng.sampling.distribution; 018 019import org.apache.commons.rng.UniformRandomProvider; 020 021/** 022 * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution</a>. 023 * 024 * <ul> 025 * <li> 026 * For small means, a Poisson process is simulated using uniform deviates, as described in 027 * <blockquote> 028 * Knuth (1969). <i>Seminumerical Algorithms</i>. The Art of Computer Programming, 029 * Volume 2. Chapter 3.4.1.F.3 Important integer-valued distributions: The Poisson distribution. 030 * Addison Wesley. 031 * </blockquote> 032 * The Poisson process (and hence, the returned value) is bounded by {@code 1000 * mean}. 033 * </li> 034 * </ul> 035 * 036 * <p>This sampler is suitable for {@code mean < 40}. 037 * For large means, {@link LargeMeanPoissonSampler} should be used instead.</p> 038 * 039 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()} and requires on average 040 * {@code mean + 1} deviates per sample.</p> 041 * 042 * @since 1.1 043 */ 044public class SmallMeanPoissonSampler 045 implements SharedStateDiscreteSampler { 046 /** 047 * Pre-compute {@code Math.exp(-mean)}. 048 * Note: This is the probability of the Poisson sample {@code P(n=0)}. 049 */ 050 private final double p0; 051 /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */ 052 private final int limit; 053 /** Underlying source of randomness. */ 054 private final UniformRandomProvider rng; 055 056 /** 057 * Create an instance. 058 * 059 * @param rng Generator of uniformly distributed random numbers. 060 * @param mean Mean. 061 * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0} 062 */ 063 public SmallMeanPoissonSampler(UniformRandomProvider rng, 064 double mean) { 065 this(rng, mean, computeP0(mean)); 066 } 067 068 /** 069 * Instantiates a new small mean poisson sampler. 070 * 071 * @param rng Generator of uniformly distributed random numbers. 072 * @param mean Mean. 073 * @param p0 {@code Math.exp(-mean)}. 074 */ 075 private SmallMeanPoissonSampler(UniformRandomProvider rng, 076 double mean, 077 double p0) { 078 this.rng = rng; 079 this.p0 = p0; 080 // The returned sample is bounded by 1000 * mean 081 limit = (int) Math.ceil(1000 * mean); 082 } 083 084 /** 085 * @param rng Generator of uniformly distributed random numbers. 086 * @param source Source to copy. 087 */ 088 private SmallMeanPoissonSampler(UniformRandomProvider rng, 089 SmallMeanPoissonSampler source) { 090 this.rng = rng; 091 p0 = source.p0; 092 limit = source.limit; 093 } 094 095 /** {@inheritDoc} */ 096 @Override 097 public int sample() { 098 int n = 0; 099 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}