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

Here are some JMH timing results:
meanrandomSourceNamesamplerTypeMethodScoreErrorMedianRuntime
   baseline2.0450.0042.045 
1SPLIT_MIX_64KempSmallMeanPoissonSamplersample17.1550.24417.14315.098
1SPLIT_MIX_64SmallMeanPoissonSamplersample19.5860.16919.58017.535
2SPLIT_MIX_64KempSmallMeanPoissonSamplersample22.1670.14822.15820.113
2SPLIT_MIX_64SmallMeanPoissonSamplersample24.2560.13824.26422.219
4SPLIT_MIX_64KempSmallMeanPoissonSamplersample32.1970.24232.17830.132
4SPLIT_MIX_64SmallMeanPoissonSamplersample28.1650.22928.17126.126
8SPLIT_MIX_64KempSmallMeanPoissonSamplersample56.9590.10956.96254.917
8SPLIT_MIX_64SmallMeanPoissonSamplersample39.53317.95937.42935.384
16SPLIT_MIX_64KempSmallMeanPoissonSamplersample107.8360.157107.832105.787
16SPLIT_MIX_64SmallMeanPoissonSamplersample58.02326.62254.88252.837
32SPLIT_MIX_64KempSmallMeanPoissonSamplersample211.4230.365211.453209.408
32SPLIT_MIX_64SmallMeanPoissonSamplersample89.9350.37989.90087.855
1WELL_44497_BKempSmallMeanPoissonSamplersample31.9370.22231.91629.871
1WELL_44497_BSmallMeanPoissonSamplersample49.3060.61349.25647.211
2WELL_44497_BKempSmallMeanPoissonSamplersample36.7350.21236.73534.690
2WELL_44497_BSmallMeanPoissonSamplersample67.6690.55367.61565.570
4WELL_44497_BKempSmallMeanPoissonSamplersample47.5961.53247.41445.369
4WELL_44497_BSmallMeanPoissonSamplersample101.4870.586101.51599.470
8WELL_44497_BKempSmallMeanPoissonSamplersample72.0940.32172.12470.079
8WELL_44497_BSmallMeanPoissonSamplersample167.5690.761167.561165.515
16WELL_44497_BKempSmallMeanPoissonSamplersample123.3200.615123.283121.238
16WELL_44497_BSmallMeanPoissonSamplersample298.6981.452298.780296.734
32WELL_44497_BKempSmallMeanPoissonSamplersample226.7990.678226.761224.716
32WELL_44497_BSmallMeanPoissonSamplersample562.2601.813562.136560.091
!kemp.jpg!
The plot shows the Runtime (which is the median minus the median of the baseline).
For reference here are the two algorithms:
{code:java}
public int sampleKemp() {
double u = rng.nextDouble();
int x = 0;
double p = p0;
while (u > p) {
u = p;
p = p * mean / ++x;
if (p == 0) {
break;
}
}
return x;
}
public int sampleSmallMean() {
int n = 0;
double r = 1;
while (n < limit) {
r *= rng.nextDouble();
if (r >= p0) {
n++;
} else {
break;
}
}
return n;
}
{code}
Conclusions:
* The speed of the RNG is very important for the the {{SmallMean}}. This is obvious as it
uses a uniform deviate within the sample loop. What is not obvious is that with a fast RNG
it can outperform the Kemp method.
* The Kemp method runtime is largely independent of the RNG. This is obvious as it only uses
a single uniform deviate.
* The Kemp method always outperforms the SmallMean method when the RNG is slow.
* The Kemp method always outperforms the SmallMean method when the RNG is fast and the mean
is <=2.
* The SmallMean method outperforms the Kemp method when the RNG is fast and the mean is above
2.
This creates some dilemmas about what to do with the Kemp method. It is a viable generator
of Poisson samples. Speed will scale the same irrespective of the generator. But is does not
provide maximum performance for the expected use case range of 0 to 40 (Note: above a mean
of 40 the LargeMeanPoissonSampler is recommended).
It may be worth using in the library for means below 1. For example a Poisson sample with
a fractional mean is used in the LargeMeanPoissonSampler.
> 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
>
>
> 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)
