commons-issues mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Alex D Herbert (JIRA)" <j...@apache.org>
Subject [jira] [Created] (RNG-91) Kemp small mean poisson sampler
Date Fri, 12 Apr 2019 14:29:00 GMT
Alex D Herbert created RNG-91:
---------------------------------

             Summary: Kemp small mean poisson sampler
                 Key: RNG-91
                 URL: https://issues.apache.org/jira/browse/RNG-91
             Project: Commons RNG
          Issue Type: New Feature
          Components: sampling
    Affects Versions: 1.3
            Reporter: Alex D Herbert
            Assignee: Alex D Herbert
             Fix For: 1.3


The current algorithm for the {{SmallMeanPoissonSampler}} is used to generate Poisson samples
for any mean up to 40. The sampler requires approximately n samples from a RNG to generate
1 Poisson deviate for a mean of n.

The [Kemp (1981)|https://www.jstor.org/stable/2346348] algorithm requires 1 sample from the
RNG and then accrues a cumulative probability using a recurrence relation to compute each
successive Poisson probability:
{noformat}
p(n+1) = p(n) * mean / (n+1)
{noformat}
The full algorithm is here:
{code:java}
    mean = ...;
    final double p0 = Math.exp(-mean);

    @Override
    public int sample() {
        double u = rng.nextDouble();
        int x = 0;
        double p = p0;
        // The algorithm listed in Kemp (1981) does not check that the rolling probability
        // is positive. This check is added to ensure no errors when the limit of the summation
        // 1 - sum(p(x)) is above 0 due to cumulative error in floating point arithmetic.
        while (u > p && p != 0) {
            u -= p;
            x++;
            p = p * mean / x;
        }
        return x;
    }
{code}
The limit for the sampler is set by the ability to compute p0. This is approximately 744.440
when Math.exp(-mean) returns 0.

A conservative limit of 700 sets an initial probability p0 of 9.85967654375977E-305. When
run through the summation series for the limit (u initialised to 1) the result when the summation
ends (p is zero) leaves u = 3.335439283623915E-15. This is close to the expected tolerance
for floating point error (Note: 1 - Math.nextDown(1) = 1.1102230246251565E-16).

Using a mean of 10 leaves u = 4.988586742717954E-17. So smaller means have a similar error.
The error in the cumulative sum is expected to result in truncation of the long tail of the
Poisson distribution (which should be bounded at infinity).

This sampler should outperform the current {{SmallMeanPoissonSampler}} as it requires 1 uniform
deviate per sample.

Note that the \{[SmallMeanPoissonSampler}} uses a limit for the mean of Integer.MAX_VALUE
/ 2. This should be updated since it also relies on p0 being above zero.



--
This message was sent by Atlassian JIRA
(v7.6.3#76005)

Mime
View raw message