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 * </ul>
031 *
032 * @since 1.1
033 *
034 * This sampler is suitable for {@code mean < 40}.
035 * For large means, {@link LargeMeanPoissonSampler} should be used instead.
036 */
037public class SmallMeanPoissonSampler
038    implements DiscreteSampler {
039    /** Upper bound to avoid truncation. */
040    private static final double MAX_MEAN = 0.5 * Integer.MAX_VALUE;
041    /**
042     * Pre-compute {@code Math.exp(-mean)}.
043     * Note: This is the probability of the Poisson sample {@code P(n=0)}.
044     */
045    private final double p0;
046    /** Pre-compute {@code 1000 * mean} as the upper limit of the sample. */
047    private final int limit;
048    /** Underlying source of randomness. */
049    private final UniformRandomProvider rng;
050
051    /**
052     * @param rng  Generator of uniformly distributed random numbers.
053     * @param mean Mean.
054     * @throws IllegalArgumentException if {@code mean <= 0}.
055     */
056    public SmallMeanPoissonSampler(UniformRandomProvider rng,
057                                   double mean) {
058        this.rng = rng;
059        if (mean <= 0) {
060            throw new IllegalArgumentException(mean + " <= " + 0);
061        }
062        if (mean > MAX_MEAN) {
063            throw new IllegalArgumentException(mean + " > " + MAX_MEAN);
064        }
065
066        p0 = Math.exp(-mean);
067        // The returned sample is bounded by 1000 * mean or Integer.MAX_VALUE
068        limit = (int) Math.ceil(Math.min(1000 * mean, Integer.MAX_VALUE));
069    }
070
071    /** {@inheritDoc} */
072    @Override
073    public int sample() {
074        int n = 0;
075        double r = 1;
076
077        while (n < limit) {
078            r *= rng.nextDouble();
079            if (r >= p0) {
080                n++;
081            } else {
082                break;
083            }
084        }
085        return n;
086    }
087
088    /** {@inheritDoc} */
089    @Override
090    public String toString() {
091        return "Small Mean Poisson deviate [" + rng.toString() + "]";
092    }
093}