Uploaded image for project: 'Commons RNG'
  1. Commons RNG
  2. RNG-91

Kemp small mean poisson sampler



    • New Feature
    • Status: Closed
    • Minor
    • Resolution: Implemented
    • 1.3
    • 1.3
    • sampling
    • None


      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) algorithm requires 1 sample from the RNG and then accrues a cumulative probability using a recurrence relation to compute each successive Poisson probability:

      p(n+1) = p(n) * mean / (n+1)

      The full algorithm is here:

          mean = ...;
          final double p0 = Math.exp(-mean);
          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;
                  p = p * mean / x;
              return x;

      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.


        1. kemp.jpg
          18 kB
          Alex Herbert
        2. poisson-samplers.jpg
          19 kB
          Alex Herbert

        Issue Links



              aherbert Alex Herbert
              aherbert Alex Herbert
              0 Vote for this issue
              1 Start watching this issue



                Time Tracking

                  Original Estimate - Not Specified
                  Not Specified
                  Remaining Estimate - 0h
                  Time Spent - 1h 50m
                  1h 50m