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