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
023 * distribution</a>.
024 *
025 * <ul>
026 *   <li>
027 *     Kemp, A, W, (1981) Efficient Generation of Logarithmically Distributed
028 *     Pseudo-Random Variables. Journal of the Royal Statistical Society. Vol. 30, No. 3, pp.
029 *     249-253.
030 *   </li>
031 * </ul>
032 *
033 * <p>This sampler is suitable for {@code mean < 40}. For large means,
034 * {@link LargeMeanPoissonSampler} should be used instead.</p>
035 *
036 * <p>Note: The algorithm uses a recurrence relation to compute the Poisson probability
037 * and a rolling summation for the cumulative probability. When the mean is large the
038 * initial probability (Math.exp(-mean)) is zero and an exception is raised by the
039 * constructor.</p>
040 *
041 * <p>Sampling uses 1 call to {@link UniformRandomProvider#nextDouble()}. This method provides
042 * an alternative to the {@link SmallMeanPoissonSampler} for slow generators of {@code double}.</p>
043 *
044 * @see <a href="https://www.jstor.org/stable/2346348">Kemp, A.W. (1981) JRSS Vol. 30, pp.
045 * 249-253</a>
046 * @since 1.3
047 */
048public final class KempSmallMeanPoissonSampler
049    implements SharedStateDiscreteSampler {
050    /** Underlying source of randomness. */
051    private final UniformRandomProvider rng;
052    /**
053     * Pre-compute {@code Math.exp(-mean)}.
054     * Note: This is the probability of the Poisson sample {@code p(x=0)}.
055     */
056    private final double p0;
057    /**
058     * The mean of the Poisson sample.
059     */
060    private final double mean;
061
062    /**
063     * @param rng Generator of uniformly distributed random numbers.
064     * @param p0 Probability of the Poisson sample {@code p(x=0)}.
065     * @param mean Mean.
066     */
067    private KempSmallMeanPoissonSampler(UniformRandomProvider rng,
068                                        double p0,
069                                        double mean) {
070        this.rng = rng;
071        this.p0 = p0;
072        this.mean = mean;
073    }
074
075    /** {@inheritDoc} */
076    @Override
077    public int sample() {
078        // Note on the algorithm:
079        // - X is the unknown sample deviate (the output of the algorithm)
080        // - x is the current value from the distribution
081        // - p is the probability of the current value x, p(X=x)
082        // - u is effectively the cumulative probability that the sample X
083        //   is equal or above the current value x, p(X>=x)
084        // So if p(X>=x) > p(X=x) the sample must be above x, otherwise it is x
085        double u = rng.nextDouble();
086        int x = 0;
087        double p = p0;
088        while (u > p) {
089            u -= p;
090            // Compute the next probability using a recurrence relation.
091            // p(x+1) = p(x) * mean / (x+1)
092            p *= mean / ++x;
093            // The algorithm listed in Kemp (1981) does not check that the rolling probability
094            // is positive. This check is added to ensure no errors when the limit of the summation
095            // 1 - sum(p(x)) is above 0 due to cumulative error in floating point arithmetic.
096            if (p == 0) {
097                return x;
098            }
099        }
100        return x;
101    }
102
103    /** {@inheritDoc} */
104    @Override
105    public String toString() {
106        return "Kemp Small Mean Poisson deviate [" + rng.toString() + "]";
107    }
108
109    /** {@inheritDoc} */
110    @Override
111    public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
112        return new KempSmallMeanPoissonSampler(rng, p0, mean);
113    }
114
115    /**
116     * Creates a new sampler for the Poisson distribution.
117     *
118     * @param rng Generator of uniformly distributed random numbers.
119     * @param mean Mean of the distribution.
120     * @return the sampler
121     * @throws IllegalArgumentException if {@code mean <= 0} or
122     * {@code Math.exp(-mean) == 0}.
123     */
124    public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
125                                                double mean) {
126        if (mean <= 0) {
127            throw new IllegalArgumentException("Mean is not strictly positive: " + mean);
128        }
129
130        final double p0 = Math.exp(-mean);
131
132        // Probability must be positive. As mean increases then p(0) decreases.
133        if (p0 > 0) {
134            return new KempSmallMeanPoissonSampler(rng, p0, mean);
135        }
136
137        // This catches the edge case of a NaN mean
138        throw new IllegalArgumentException("No probability for mean: " + mean);
139    }
140}