[ https://issues.apache.org/jira/browse/RNG-91?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]
Alex D Herbert updated RNG-91: ------------------------------ Attachment: poisson-samplers.jpg > 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 > Priority: Minor > Fix For: 1.3 > > Attachments: kemp.jpg, poisson-samplers.jpg > > Time Spent: 20m > Remaining Estimate: 0h > > 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)