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  import org.apache.commons.rng.sampling.distribution.InternalUtils.FactorialLog;
21  
22  /**
23   * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution</a>.
24   *
25   * <ul>
26   *  <li>
27   *   For large means, we use the rejection algorithm described in
28   *   <blockquote>
29   *    Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i><br>
30   *    <strong>Computing</strong> vol. 26 pp. 197-207.
31   *   </blockquote>
32   *  </li>
33   * </ul>
34   *
35   * @since 1.1
36   *
37   * This sampler is suitable for {@code mean >= 40}.
38   */
39  public class LargeMeanPoissonSampler
40      implements DiscreteSampler {
41  
42      /** Class to compute {@code log(n!)}. This has no cached values. */
43      private static final InternalUtils.FactorialLog NO_CACHE_FACTORIAL_LOG;
44  
45      static {
46          // Create without a cache.
47          NO_CACHE_FACTORIAL_LOG = FactorialLog.create();
48      }
49  
50      /** Underlying source of randomness. */
51      private final UniformRandomProvider rng;
52      /** Exponential. */
53      private final ContinuousSampler exponential;
54      /** Gaussian. */
55      private final ContinuousSampler gaussian;
56      /** Local class to compute {@code log(n!)}. This may have cached values. */
57      private final InternalUtils.FactorialLog factorialLog;
58  
59      // Working values
60  
61      /** Algorithm constant: {@code Math.floor(mean)}. */
62      private final double lambda;
63      /** Algorithm constant: {@code mean - lambda}. */
64      private final double lambdaFractional;
65      /** Algorithm constant: {@code Math.log(lambda)}. */
66      private final double logLambda;
67      /** Algorithm constant: {@code factorialLog((int) lambda)}. */
68      private final double logLambdaFactorial;
69      /** Algorithm constant: {@code Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1))}. */
70      private final double delta;
71      /** Algorithm constant: {@code delta / 2}. */
72      private final double halfDelta;
73      /** Algorithm constant: {@code 2 * lambda + delta}. */
74      private final double twolpd;
75      /**
76       * Algorithm constant: {@code a1 / aSum} with
77       * <ul>
78       *  <li>{@code a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(c1)}</li>
79       *  <li>{@code aSum = a1 + a2 + 1}</li>
80       * </ul>
81       */
82      private final double p1;
83      /**
84       * Algorithm constant: {@code a1 / aSum} with
85       * <ul>
86       *  <li>{@code a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd)}</li>
87       *  <li>{@code aSum = a1 + a2 + 1}</li>
88       * </ul>
89       */
90      private final double p2;
91      /** Algorithm constant: {@code 1 / (8 * lambda)}. */
92      private final double c1;
93  
94      /** The internal Poisson sampler for the lambda fraction. */
95      private final DiscreteSampler smallMeanPoissonSampler;
96  
97      /**
98       * @param rng  Generator of uniformly distributed random numbers.
99       * @param mean Mean.
100      * @throws IllegalArgumentException if {@code mean <= 0}.
101      */
102     public LargeMeanPoissonSampler(UniformRandomProvider rng,
103                                    double mean) {
104         this.rng = rng;
105         if (mean <= 0) {
106             throw new IllegalArgumentException(mean + " <= " + 0);
107         }
108 
109         gaussian = new ZigguratNormalizedGaussianSampler(rng);
110         exponential = new AhrensDieterExponentialSampler(rng, 1);
111         // Plain constructor uses the uncached function.
112         factorialLog = NO_CACHE_FACTORIAL_LOG;
113 
114         // Cache values used in the algorithm
115         lambda = Math.floor(mean);
116         lambdaFractional = mean - lambda;
117         logLambda = Math.log(lambda);
118         logLambdaFactorial = factorialLog((int) lambda);
119         delta = Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1));
120         halfDelta = delta / 2;
121         twolpd = 2 * lambda + delta;
122         c1 = 1 / (8 * lambda);
123         final double a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(c1);
124         final double a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd);
125         final double aSum = a1 + a2 + 1;
126         p1 = a1 / aSum;
127         p2 = a2 / aSum;
128 
129         // The algorithm requires a Poisson sample from the remaining lambda fraction.
130         smallMeanPoissonSampler = (lambdaFractional < Double.MIN_VALUE) ?
131             null : // Not used.
132             new SmallMeanPoissonSampler(rng, lambdaFractional);
133     }
134 
135     /** {@inheritDoc} */
136     @Override
137     public int sample() {
138 
139         final int y2 = (smallMeanPoissonSampler == null) ?
140             0 : // No lambda fraction
141             smallMeanPoissonSampler.sample();
142 
143         double x = 0;
144         double y = 0;
145         double v = 0;
146         int a = 0;
147         double t = 0;
148         double qr = 0;
149         double qa = 0;
150         while (true) {
151             final double u = rng.nextDouble();
152             if (u <= p1) {
153                 final double n = gaussian.sample();
154                 x = n * Math.sqrt(lambda + halfDelta) - 0.5d;
155                 if (x > delta || x < -lambda) {
156                     continue;
157                 }
158                 y = x < 0 ? Math.floor(x) : Math.ceil(x);
159                 final double e = exponential.sample();
160                 v = -e - 0.5 * n * n + c1;
161             } else {
162                 if (u > p1 + p2) {
163                     y = lambda;
164                     break;
165                 }
166                 x = delta + (twolpd / delta) * exponential.sample();
167                 y = Math.ceil(x);
168                 v = -exponential.sample() - delta * (x + 1) / twolpd;
169             }
170             a = x < 0 ? 1 : 0;
171             t = y * (y + 1) / (2 * lambda);
172             if (v < -t && a == 0) {
173                 y = lambda + y;
174                 break;
175             }
176             qr = t * ((2 * y + 1) / (6 * lambda) - 1);
177             qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
178             if (v < qa) {
179                 y = lambda + y;
180                 break;
181             }
182             if (v > qr) {
183                 continue;
184             }
185             if (v < y * logLambda - factorialLog((int) (y + lambda)) + logLambdaFactorial) {
186                 y = lambda + y;
187                 break;
188             }
189         }
190 
191         return (int) Math.min(y2 + (long) y, Integer.MAX_VALUE);
192     }
193 
194     /**
195      * Compute the natural logarithm of the factorial of {@code n}.
196      *
197      * @param n Argument.
198      * @return {@code log(n!)}
199      * @throws IllegalArgumentException if {@code n < 0}.
200      */
201     private double factorialLog(int n) {
202         return factorialLog.value(n);
203     }
204 
205     /** {@inheritDoc} */
206     @Override
207     public String toString() {
208         return "Large Mean Poisson deviate [" + rng.toString() + "]";
209     }
210 }