View Javadoc
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       * Create an instance.
58       *
59       * @param rng  Generator of uniformly distributed random numbers.
60       * @param mean Mean.
61       * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0}
62       */
63      public SmallMeanPoissonSampler(UniformRandomProvider rng,
64                                     double mean) {
65          this(rng, mean, computeP0(mean));
66      }
67  
68      /**
69       * Instantiates a new small mean poisson sampler.
70       *
71       * @param rng  Generator of uniformly distributed random numbers.
72       * @param mean Mean.
73       * @param p0 {@code Math.exp(-mean)}.
74       */
75      private SmallMeanPoissonSampler(UniformRandomProvider rng,
76                                      double mean,
77                                      double p0) {
78          this.rng = rng;
79          this.p0 = p0;
80          // The returned sample is bounded by 1000 * mean
81          limit = (int) Math.ceil(1000 * mean);
82      }
83  
84      /**
85       * @param rng Generator of uniformly distributed random numbers.
86       * @param source Source to copy.
87       */
88      private SmallMeanPoissonSampler(UniformRandomProvider rng,
89                                      SmallMeanPoissonSampler source) {
90          this.rng = rng;
91          p0 = source.p0;
92          limit = source.limit;
93      }
94  
95      /** {@inheritDoc} */
96      @Override
97      public int sample() {
98          int n = 0;
99          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 }