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 * <li>
35 * For large means, we use the rejection algorithm described in
36 * <blockquote>
37 * Devroye, Luc. (1981). <i>The Computer Generation of Poisson Random Variables</i><br>
38 * <strong>Computing</strong> vol. 26 pp. 197-207.
39 * </blockquote>
40 * </li>
41 * </ul>
42 *
43 * <p>Sampling uses:</p>
44 *
45 * <ul>
46 * <li>{@link UniformRandomProvider#nextDouble()}
47 * <li>{@link UniformRandomProvider#nextLong()} (large means only)
48 * </ul>
49 *
50 * @since 1.0
51 */
52 public class PoissonSampler
53 extends SamplerBase
54 implements SharedStateDiscreteSampler {
55
56 /**
57 * Value for switching sampling algorithm.
58 *
59 * <p>Package scope for the {@link PoissonSamplerCache}.
60 */
61 static final double PIVOT = 40;
62 /** The internal Poisson sampler. */
63 private final SharedStateDiscreteSampler poissonSamplerDelegate;
64
65 /**
66 * This instance delegates sampling. Use the factory method
67 * {@link #of(UniformRandomProvider, double)} to create an optimal sampler.
68 *
69 * @param rng Generator of uniformly distributed random numbers.
70 * @param mean Mean.
71 * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean > 0.5 *}
72 * {@link Integer#MAX_VALUE}.
73 */
74 public PoissonSampler(UniformRandomProvider rng,
75 double mean) {
76 super(null);
77
78 // Delegate all work to specialised samplers.
79 poissonSamplerDelegate = of(rng, mean);
80 }
81
82 /** {@inheritDoc} */
83 @Override
84 public int sample() {
85 return poissonSamplerDelegate.sample();
86 }
87
88 /** {@inheritDoc} */
89 @Override
90 public String toString() {
91 return poissonSamplerDelegate.toString();
92 }
93
94 /**
95 * {@inheritDoc}
96 *
97 * @since 1.3
98 */
99 @Override
100 public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
101 // Direct return of the optimised sampler
102 return poissonSamplerDelegate.withUniformRandomProvider(rng);
103 }
104
105 /**
106 * Creates a new Poisson distribution sampler.
107 *
108 * @param rng Generator of uniformly distributed random numbers.
109 * @param mean Mean.
110 * @return the sampler
111 * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean > 0.5 *}
112 * {@link Integer#MAX_VALUE}.
113 * @since 1.3
114 */
115 public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
116 double mean) {
117 // Each sampler should check the input arguments.
118 return mean < PIVOT ?
119 SmallMeanPoissonSampler.of(rng, mean) :
120 LargeMeanPoissonSampler.of(rng, mean);
121 }
122 }