From 2c8a114f766d05929e908fd79c5e4baf5a3841ae Mon Sep 17 00:00:00 2001 From: Phil Steitz Date: Mon, 12 Oct 2009 02:04:06 +0000 Subject: [PATCH] Implemented alternative algorithm for generating poisson deviates when the mean is large. JIRA: MATH-294. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@824214 13f79535-47bb-0310-9956-ffa450edef68 --- .../commons/math/random/RandomDataImpl.java | 149 ++++++++---------- src/site/xdoc/changes.xml | 4 + .../commons/math/random/RandomDataTest.java | 22 ++- 3 files changed, 83 insertions(+), 92 deletions(-) diff --git a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java index 18097bcca..d5d2474f7 100644 --- a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java +++ b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java @@ -322,30 +322,17 @@ public class RandomDataImpl implements RandomData, Serializable { /** * {@inheritDoc} *

- * Algorithm Description: For small means, uses simulation - * of a Poisson process using Uniform deviates, as described here. - *

- *

- * The Poisson process (and hence value returned) is bounded by 1000 * mean. - *

+ * Algorithm Description: + *

* - *

- * References: - *

- *

- * - * @param mean - * mean of the Poisson distribution. + * @param mean mean of the Poisson distribution. * @return the random Poisson value. */ public long nextPoisson(double mean) { @@ -356,7 +343,7 @@ public class RandomDataImpl implements RandomData, Serializable { final RandomGenerator generator = getRan(); - double pivot = 6.0; + final double pivot = 40.0d; if (mean < pivot) { double p = Math.exp(-mean); long n = 0; @@ -374,68 +361,70 @@ public class RandomDataImpl implements RandomData, Serializable { } return n; } else { - double mu = Math.floor(mean); - double delta = Math.floor(pivot + (mu - pivot) / 2.0); // integer - // between 6 - // and mean - double mu2delta = 2.0 * mu + delta; - double muDeltaHalf = mu + delta / 2.0; - double logMeanMu = Math.log(mean / mu); + final double lambda = Math.floor(mean); + final double lambdaFractional = mean - lambda; + final double logLambda = Math.log(lambda); + final double logLambdaFactorial = MathUtils.factorialLog((int) lambda); + final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional); + final double delta = Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1)); + final double halfDelta = delta / 2; + final double twolpd = 2 * lambda + delta; + final double a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(1 / 8 * lambda); + final double a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd); + final double aSum = a1 + a2 + 1; + final double p1 = a1 / aSum; + final double p2 = a2 / aSum; + final double c1 = 1 / (8 * lambda); - double muFactorialLog = MathUtils.factorialLog((int) mu); - - double c1 = Math.sqrt(Math.PI * mu / 2.0); - double c2 = c1 + - Math.sqrt(Math.PI * muDeltaHalf / - (2.0 * Math.exp(1.0 / mu2delta))); - double c3 = c2 + 2.0; - double c4 = c3 + Math.exp(1.0 / 78.0); - double c = c4 + 2.0 / delta * mu2delta * - Math.exp(-delta / mu2delta * (1.0 + delta / 2.0)); - - double y = 0.0; - double x = 0.0; - double w = Double.POSITIVE_INFINITY; - - boolean accept = false; - while (!accept) { - double u = nextUniform(0.0, c); - double e = nextExponential(mean); - - if (u <= c1) { - double z = nextGaussian(0.0, 1.0); - y = -Math.abs(z) * Math.sqrt(mu) - 1.0; - x = Math.floor(y); - w = -z * z / 2.0 - e - x * logMeanMu; - if (x < -mu) { - w = Double.POSITIVE_INFINITY; + double x = 0; + double y = 0; + double v = 0; + int a = 0; + double t = 0; + double qr = 0; + double qa = 0; + for (;;) { + final double u = nextUniform(0.0, 1); + if (u <= p1) { + final double n = nextGaussian(0d, 1d); + x = n * Math.sqrt(lambda + halfDelta) - 0.5d; + if (x > delta || x < -lambda) { + continue; } - } else if (c1 < u && u <= c2) { - double z = nextGaussian(0.0, 1.0); - y = 1.0 + Math.abs(z) * Math.sqrt(muDeltaHalf); - x = Math.ceil(y); - w = (-y * y + 2.0 * y) / mu2delta - e - x * logMeanMu; - if (x > delta) { - w = Double.POSITIVE_INFINITY; + y = x < 0 ? Math.floor(x) : Math.ceil(x); + final double e = nextExponential(1d); + v = -e - (n * n / 2) + c1; + } else { + if (u > p1 + p2) { + y = lambda; + break; + } else { + x = delta + (twolpd / delta) * nextExponential(1d); + y = Math.ceil(x); + v = -nextExponential(1d) - delta * (x + 1) / twolpd; } - } else if (c2 < u && u <= c3) { - x = 0.0; - w = -e; - } else if (c3 < u && u <= c4) { - x = 1.0; - w = -e - logMeanMu; - } else if (c4 < u) { - double v = nextExponential(mean); - y = delta + v * 2.0 / delta * mu2delta; - x = Math.ceil(y); - w = -delta / mu2delta * (1.0 + y / 2.0) - e - x * logMeanMu; } - accept = w <= x * Math.log(mu) - - MathUtils.factorialLog((int) (mu + x)) / muFactorialLog; + a = x < 0 ? 1 : 0; + t = y * (y + 1) / (2 * lambda); + if (v < -t && a == 0) { + y = lambda + y; + break; + } + qr = t * ((2 * y + 1) / (6 * lambda) - 1); + qa = qr - (t * t) / (3 * (lambda + a * (y + 1))); + if (v < qa) { + y = lambda + y; + break; + } + if (v > qr) { + continue; + } + if (v < y * logLambda - MathUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) { + y = lambda + y; + break; + } } - // cast to long is acceptable because both x and mu are whole - // numbers. - return (long) (x + mu); + return y2 + (long) y; } } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 56e03fe2e..a7c896045 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -39,6 +39,10 @@ The type attribute can be add,update,fix,remove. + + Fixed implementation of RandomDataImpl#nextPoisson by implementing an alternative + algorithm for large means. + Added support for multidimensional interpolation using the robust microsphere algorithm. diff --git a/src/test/java/org/apache/commons/math/random/RandomDataTest.java b/src/test/java/org/apache/commons/math/random/RandomDataTest.java index 678ed0605..d75636476 100644 --- a/src/test/java/org/apache/commons/math/random/RandomDataTest.java +++ b/src/test/java/org/apache/commons/math/random/RandomDataTest.java @@ -228,14 +228,20 @@ public class RandomDataTest extends RetryTestCase { } public void testNextPoissionConistency() throws Exception { - // TODO: once MATH-294 is resolved, increase upper bounds on test means - for (int i = 1; i < 6; i++) { + // Small integral means + for (int i = 1; i < 100; i++) { checkNextPoissonConsistency(i); } + // non-integer means RandomData randomData = new RandomDataImpl(); for (int i = 1; i < 10; i++) { - checkNextPoissonConsistency(randomData.nextUniform(1, 6)); + checkNextPoissonConsistency(randomData.nextUniform(1, 1000)); } + // large means + // TODO: When MATH-282 is resolved, s/3000/10000 below + for (int i = 1; i < 10; i++) { + checkNextPoissonConsistency(randomData.nextUniform(1000, 3000)); + } } /** @@ -367,15 +373,7 @@ public class RandomDataTest extends RetryTestCase { msgBuffer.append(alpha); msgBuffer.append("."); fail(msgBuffer.toString()); - } - - } - - public void testNextPoissonLargeMean() { - for (int i = 0; i < 1000; i++) { - long n = randomData.nextPoisson(1500.0); - assertTrue(0 <= n); - } + } } /** test dispersion and failute modes for nextHex() */