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 }