[ https://issues.apache.org/jira/browse/RNG91?page=com.atlassian.jira.plugin.system.issuetabpanels:commenttabpanel&focusedCommentId=16819878#comment16819878
]
Alex D Herbert commented on RNG91:

Original timings above were using the Kemp method that effectively computes the cumulative
probability of the Poisson distribution until the cumulative probability is above the value
of the random deviate.
The cumulative probability is computed *for each sample*. This is a lot of repeat computation.
However it should be noted that the value of the random deviate will be above 0.5 half of
the time. Therefore for repeat sampling it is worth computing and storing the cumulative probability
and the sample value that corresponds to approximately 0.5 in the cumulative probability function.
This will be known as the 50% hedge value.
This will *save half of the computation* for *half of the samples*. This is expected to produce
a speed increase of 50% and so *reduce runtime to 2/3* of that listed above.
I have tested an implementation that uses this 50% hedge.
Here is an output of the raw data for various means using the original Kemp algorithm and
the modified Kemp algorithm. The algorithm is compared to the current Small and Large mean
Poisson sampler.
Note that the XO_RO_SHI_RO_128_PLUS generator is the fastest currently in the library for
double floating point generation. The Small mean sampler is very dependent on the speed of
the generator so a slow algorithm is once again included to demonstrate this (WELL_44497B).
 Small  OriginalKemp  Kemp  Large  
meanXoSplitWellXoSplitWellXoSplitWellXoSplitWell
0.258.6310.0330.147.879.0722.168.139.0422.04   
0.512.3114.6037.7211.3413.0725.4211.0612.3925.09   
116.7119.9249.5315.2017.2529.5913.5815.4128.0491.78100.64203.57
221.3723.1467.8220.6421.0036.4016.7318.7031.3888.5197.72193.22
425.9428.17101.9326.9127.2243.4621.5223.5139.2481.6788.89172.74
835.4937.44167.9444.2644.5860.7233.0436.9149.2172.9679.27145.84
1652.6155.00300.2378.5579.2894.8153.5855.0368.8664.8969.99134.86
3287.1790.55562.73147.19147.15162.9290.8692.52107.1059.9462.67146.87
Sorry for the large table width but the raw data allows a lot of analysis to be done. Algorithms
are:
* Split = SPLIT_MIX_64
* Xo = XO_RO_SHI_RO_128_PLUS
* Well = WELL_44497_B
The LargeMean sampler is invalid at a mean below 1.
*Analysis*
The new hedge computation is very similar to the original Kemp sampler at very low mean. The
added complexity of the algorithm cancels the speed increase of the 50% hedge. It should also
be noted that at a mean of 0.25 the cumulative probability of a Poisson distribution is 0.779
for n=0. This high skew in the distribution effectively means the hedge is rarely used and
the complexity of the algorithm is a waste.
When the mean is large the new 50% hedge algorithm does show the approximate runtime at 2/3
of the original Kemp algorithm. So the hedge algorithm is working as expected,
When the mean is in the range 0.5 to 8 the new hedged Kemp sampler is the fastest.
Only at mean 16 and 32 is the current Small mean sampler faster and this is by less then 5%
and requires a very fast generator. When the generator is slow the new Kemp sampler is faster
by 5fold due to the small sampler requiring multiple double values from the RNG per poisson
sample.
The Large mean sampler is currently used in the library at mean 40. This table shows that
with a fast RNG the sampler starts outperforming the small samplers around a mean of 32. When
the sampler is slow it outperforms the small mean sampler at mean 8 but is still slower than
the new Kemp sampler when the mean is 32.
Note: The bound of 40 on the large mean poisson sampler may not be based on performance and
instead on the validity of the algorithm. It is based around approximating the poisson distribution
using a normal distribution with handling of the skew when the mean is low. This sampler may
not even be valid for very low means and the current units tests use a mean of 67.89 to
test the distribution.
This is perhaps best illustrated with a chart (this omits the original Kemp sampler and only
includes the new Kemp sampler with the 50% hedge):
!poissonsamplers.jpg!
The new code for the Kemp sampler and the performance test is submitted in a PR.
A recommendation would be to:
* Switch the code to use the Kemp sampler instead of the Small mean poisson sampler. Speed
is approximately equivalent with the fastest generators in the library and is much better
when the generator is slow as only 1 double from the RNG is used per sample.
* Maintain the switch point of 40 for the change to the large mean poisson sampler. This
is a valid point on the basis of performance and may be imposed anyway by the details of the
algorithm
> Kemp small mean poisson sampler
> 
>
> Key: RNG91
> URL: https://issues.apache.org/jira/browse/RNG91
> 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, poissonsamplers.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.85967654375977E305.
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.335439283623915E15. This is close to the expected
tolerance for floating point error (Note: 1  Math.nextDown(1) = 1.1102230246251565E16).
> Using a mean of 10 leaves u = 4.988586742717954E17. 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)
