001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.sampling.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution</a>.
023 *
024 * <ul>
025 *  <li>
026 *   For small means, a Poisson process is simulated using uniform deviates, as
027 *   described <a href="http://mathaa.epfl.ch/cours/PMMI2001/interactive/rng7.htm">here</a>.
028 *   The Poisson process (and hence, the returned value) is bounded by 1000 * mean.
029 *  </li>
030 *  <li>
031 *   For large means, we use the rejection algorithm described in
032 *   <blockquote>
033 *    Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i><br>
034 *    <strong>Computing</strong> vol. 26 pp. 197-207.
035 *   </blockquote>
036 *  </li>
037 * </ul>
038 */
039public class PoissonSampler
040    extends SamplerBase
041    implements DiscreteSampler {
042    /** Value for switching sampling algorithm. */
043    private static final double PIVOT = 40;
044    /** Mean of the distribution. */
045    private final double mean;
046    /** Exponential. */
047    private final ContinuousSampler exponential;
048    /** Gaussian. */
049    private final ContinuousSampler gaussian;
050    /** {@code log(n!)}. */
051    private final InternalUtils.FactorialLog factorialLog;
052
053    /**
054     * @param rng Generator of uniformly distributed random numbers.
055     * @param mean Mean.
056     * @throws IllegalArgumentException if {@code mean <= 0}.
057     */
058    public PoissonSampler(UniformRandomProvider rng,
059                          double mean) {
060        super(rng);
061        if (mean <= 0) {
062            throw new IllegalArgumentException(mean + " <= " + 0);
063        }
064
065        this.mean = mean;
066
067        gaussian = new BoxMullerGaussianSampler(rng, 0, 1);
068        exponential = new AhrensDieterExponentialSampler(rng, 1);
069        factorialLog = mean < PIVOT ?
070            null : // Not used.
071            InternalUtils.FactorialLog.create().withCache((int) Math.min(mean, 2 * PIVOT));
072    }
073
074    /** {@inheritDoc} */
075    @Override
076    public int sample() {
077        return (int) Math.min(nextPoisson(mean), Integer.MAX_VALUE);
078    }
079
080    /** {@inheritDoc} */
081    @Override
082    public String toString() {
083        return "Poisson deviate [" + super.toString() + "]";
084    }
085
086    /**
087     * @param meanPoisson Mean.
088     * @return the next sample.
089     */
090    private long nextPoisson(double meanPoisson) {
091        if (meanPoisson < PIVOT) {
092            double p = Math.exp(-meanPoisson);
093            long n = 0;
094            double r = 1;
095
096            while (n < 1000 * meanPoisson) {
097                r *= nextDouble();
098                if (r >= p) {
099                    n++;
100                } else {
101                    break;
102                }
103            }
104            return n;
105        } else {
106            final double lambda = Math.floor(meanPoisson);
107            final double lambdaFractional = meanPoisson - lambda;
108            final double logLambda = Math.log(lambda);
109            final double logLambdaFactorial = factorialLog((int) lambda);
110            final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
111            final double delta = Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1));
112            final double halfDelta = delta / 2;
113            final double twolpd = 2 * lambda + delta;
114            final double a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(1 / (8 * lambda));
115            final double a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd);
116            final double aSum = a1 + a2 + 1;
117            final double p1 = a1 / aSum;
118            final double p2 = a2 / aSum;
119            final double c1 = 1 / (8 * lambda);
120
121            double x = 0;
122            double y = 0;
123            double v = 0;
124            int a = 0;
125            double t = 0;
126            double qr = 0;
127            double qa = 0;
128            while (true) {
129                final double u = nextDouble();
130                if (u <= p1) {
131                    final double n = gaussian.sample();
132                    x = n * Math.sqrt(lambda + halfDelta) - 0.5d;
133                    if (x > delta || x < -lambda) {
134                        continue;
135                    }
136                    y = x < 0 ? Math.floor(x) : Math.ceil(x);
137                    final double e = exponential.sample();
138                    v = -e - 0.5 * n * n + c1;
139                } else {
140                    if (u > p1 + p2) {
141                        y = lambda;
142                        break;
143                    } else {
144                        x = delta + (twolpd / delta) * exponential.sample();
145                        y = Math.ceil(x);
146                        v = -exponential.sample() - delta * (x + 1) / twolpd;
147                    }
148                }
149                a = x < 0 ? 1 : 0;
150                t = y * (y + 1) / (2 * lambda);
151                if (v < -t && a == 0) {
152                    y = lambda + y;
153                    break;
154                }
155                qr = t * ((2 * y + 1) / (6 * lambda) - 1);
156                qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
157                if (v < qa) {
158                    y = lambda + y;
159                    break;
160                }
161                if (v > qr) {
162                    continue;
163                }
164                if (v < y * logLambda - factorialLog((int) (y + lambda)) + logLambdaFactorial) {
165                    y = lambda + y;
166                    break;
167                }
168            }
169            return y2 + (long) y;
170        }
171    }
172
173    /**
174     * Compute the natural logarithm of the factorial of {@code n}.
175     *
176     * @param n Argument.
177     * @return {@code log(n!)}
178     * @throws IllegalArgumentException if {@code n < 0}.
179     */
180    private double factorialLog(int n) {
181        return factorialLog.value(n);
182    }
183}