You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2019/10/02 23:52:40 UTC
[commons-rng] branch master updated: LargeMeanPoissonSampler: Added
comments of the algorithm steps.
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-rng.git
The following commit(s) were added to refs/heads/master by this push:
new eb79e2b LargeMeanPoissonSampler: Added comments of the algorithm steps.
eb79e2b is described below
commit eb79e2b0fe7fbb38ff3a5c8dd0b9397d9af1808b
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Oct 3 00:52:36 2019 +0100
LargeMeanPoissonSampler: Added comments of the algorithm steps.
The steps are taken from the reference paper by Devroye.
---
.../rng/sampling/distribution/LargeMeanPoissonSampler.java | 10 ++++++++++
1 file changed, 10 insertions(+)
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
index a8e4e54..ae16216 100644
--- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
@@ -240,8 +240,10 @@ public class LargeMeanPoissonSampler
double qr;
double qa;
while (true) {
+ // Step 1:
final double u = rng.nextDouble();
if (u <= p1) {
+ // Step 2:
final double n = gaussian.sample();
x = n * Math.sqrt(lambda + halfDelta) - 0.5d;
if (x > delta || x < -lambda) {
@@ -251,6 +253,7 @@ public class LargeMeanPoissonSampler
final double e = exponential.sample();
v = -e - 0.5 * n * n + c1;
} else {
+ // Step 3:
if (u > p1 + p2) {
y = lambda;
break;
@@ -259,21 +262,28 @@ public class LargeMeanPoissonSampler
y = Math.ceil(x);
v = -exponential.sample() - delta * (x + 1) / twolpd;
}
+ // The Squeeze Principle
+ // Step 4.1:
a = x < 0 ? 1 : 0;
t = y * (y + 1) / (2 * lambda);
+ // Step 4.2
if (v < -t && a == 0) {
y = lambda + y;
break;
}
+ // Step 4.3:
qr = t * ((2 * y + 1) / (6 * lambda) - 1);
qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
+ // Step 4.4:
if (v < qa) {
y = lambda + y;
break;
}
+ // Step 4.5:
if (v > qr) {
continue;
}
+ // Step 4.6:
if (v < y * logLambda - getFactorialLog((int) (y + lambda)) + logLambdaFactorial) {
y = lambda + y;
break;