You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ps...@apache.org on 2009/10/12 04:04:06 UTC
svn commit: r824214 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/random/RandomDataImpl.java
site/xdoc/changes.xml
test/java/org/apache/commons/math/random/RandomDataTest.java
Author: psteitz
Date: Mon Oct 12 02:04:06 2009
New Revision: 824214
URL: http://svn.apache.org/viewvc?rev=824214&view=rev
Log:
Implemented alternative algorithm for generating poisson deviates when the mean is large. JIRA: MATH-294.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/RandomDataImpl.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/RandomDataTest.java
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/RandomDataImpl.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/RandomDataImpl.java?rev=824214&r1=824213&r2=824214&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/RandomDataImpl.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/RandomDataImpl.java Mon Oct 12 02:04:06 2009
@@ -322,30 +322,17 @@
/**
* {@inheritDoc}
* <p>
- * <strong>Algorithm Description</strong>: For small means, uses simulation
- * of a Poisson process using Uniform deviates, as described <a
- * href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
- * </p>
- * <p>
- * The Poisson process (and hence value returned) is bounded by 1000 * mean.
- * </p>
+ * <strong>Algorithm Description</strong>:
+ * <ul><li> For small means, uses simulation of a Poisson process
+ * using Uniform deviates, as described
+ * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
+ * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li>
+ *
+ * <li> For large means, uses the rejection algorithm described in <br/>
+ * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
+ * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p>
*
- * <p>
- * For large means, uses a reject method as described in <a
- * href="http://cg.scs.carleton.ca/~luc/rnbookindex.html">Non-Uniform Random
- * Variate Generation</a>
- * </p>
- *
- * <p>
- * References:
- * <ul>
- * <li>Devroye, Luc. (1986). <i>Non-Uniform Random Variate Generation</i>.
- * New York, NY. Springer-Verlag</li>
- * </ul>
- * </p>
- *
- * @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 @@
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 @@
}
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);
-
- 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;
+ 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 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;
}
}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=824214&r1=824213&r2=824214&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Oct 12 02:04:06 2009
@@ -39,6 +39,10 @@
</properties>
<body>
<release version="2.1" date="TBD" description="TBD">
+ <action dev="psteitz" tyoe="fix" issue="MATH-294">
+ Fixed implementation of RandomDataImpl#nextPoisson by implementing an alternative
+ algorithm for large means.
+ </action>
<action dev="psteitz" tyoe="add" issue="MATH-300" due-to="Gilles Sadowski">
Added support for multidimensional interpolation using the robust microsphere algorithm.
</action>
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/RandomDataTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/RandomDataTest.java?rev=824214&r1=824213&r2=824214&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/RandomDataTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/random/RandomDataTest.java Mon Oct 12 02:04:06 2009
@@ -228,14 +228,20 @@
}
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 @@
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() */