[ https://issues.apache.org/jira/browse/RNG-91?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]
Alex D Herbert updated RNG-91:
------------------------------
Attachment: kemp.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
>
>
> 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)