[
https://issues.apache.org/jira/browse/RNG-91?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16819984#comment-16819984
]
Gilles commented on RNG-91:
---------------------------
{quote}
A recommendation would be to:
* \[...\] use the Kemp sampler \[...\]
* Maintain the switch point of 40 \[...\]
{quote}
+1
And deprecate {{SmallMeanPoissonSampler}}?
> 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: 40m
> 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)
