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 described in
027 *   <blockquote>
028 *    Knuth (1969). <i>Seminumerical Algorithms</i>. The Art of Computer Programming,
029 *    Volume 2. Chapter 3.4.1.F.3 Important integer-valued distributions: The Poisson distribution.
030 *    Addison Wesley.
031 *   </blockquote>
032 *   The Poisson process (and hence, the returned value) is bounded by {@code 1000 * mean}.
033 *  </li>
034 * </ul>
035 *
036 * <p>This sampler is suitable for {@code mean < 40}.
037 * For large means, {@link LargeMeanPoissonSampler} should be used instead.</p>
038 *
039 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()} and requires on average
040 * {@code mean + 1} deviates per sample.</p>
041 *
042 * @since 1.1
043 */
044public class SmallMeanPoissonSampler
045    implements SharedStateDiscreteSampler {
046    /**
047     * Pre-compute {@code Math.exp(-mean)}.
048     * Note: This is the probability of the Poisson sample {@code P(n=0)}.
049     */
050    private final double p0;
051    /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */
052    private final int limit;
053    /** Underlying source of randomness. */
054    private final UniformRandomProvider rng;
055
056    /**
057     * Create an instance.
058     *
059     * @param rng  Generator of uniformly distributed random numbers.
060     * @param mean Mean.
061     * @throws IllegalArgumentException if {@code mean <= 0} or {@code Math.exp(-mean) == 0}
062     */
063    public SmallMeanPoissonSampler(UniformRandomProvider rng,
064                                   double mean) {
065        this(rng, mean, computeP0(mean));
066    }
067
068    /**
069     * Instantiates a new small mean poisson sampler.
070     *
071     * @param rng  Generator of uniformly distributed random numbers.
072     * @param mean Mean.
073     * @param p0 {@code Math.exp(-mean)}.
074     */
075    private SmallMeanPoissonSampler(UniformRandomProvider rng,
076                                    double mean,
077                                    double p0) {
078        this.rng = rng;
079        this.p0 = p0;
080        // The returned sample is bounded by 1000 * mean
081        limit = (int) Math.ceil(1000 * mean);
082    }
083
084    /**
085     * @param rng Generator of uniformly distributed random numbers.
086     * @param source Source to copy.
087     */
088    private SmallMeanPoissonSampler(UniformRandomProvider rng,
089                                    SmallMeanPoissonSampler source) {
090        this.rng = rng;
091        p0 = source.p0;
092        limit = source.limit;
093    }
094
095    /** {@inheritDoc} */
096    @Override
097    public int sample() {
098        int n = 0;
099        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}