# commons-issues mailing list archives

##### Site index · List index
Message view
Top
From "Alex D Herbert (JIRA)" <j...@apache.org>
Subject [jira] [Updated] (RNG-91) Kemp small mean poisson sampler
Date Fri, 12 Apr 2019 15:46:00 GMT
```
[ 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)

```
Mime
View raw message