On 6/20/11 2:46 PM, Mikkel Meyer Andersen wrote: > Hi, > > I've now made a commit. I'm now writing to be sure that I did it > correctly. First of all, I of course committed the code. I also added > the following to changes.xml. I did this in the same commit (this is > best, right?):
Yes, its best to include the update to changes.xml in the same commit. > <action dev="mikl" type="fix" issue="MATH-597"> > Implemented faster generation of random exponential > distributed values with > algorithm from Ahrens and Dieter (1972): Computer methods for sampling > from the exponential and normal distributions. > </action> > > As the commit message, I used: > Added fix for MATH-597: Implemented faster generation of random > exponential distributed values with algorithm from Ahrens and Dieter > (1972): Computer methods for sampling from the exponential and normal > distributions. Test case was improved, too. > > Is this okay? Did I miss something or should I done something different? Looks great! Phil > Cheers, Mikkel. > > 2011/6/20 <[email protected]>: >> Author: mikl >> Date: Mon Jun 20 21:42:48 2011 >> New Revision: 1137795 >> >> URL: http://svn.apache.org/viewvc?rev=1137795&view=rev >> Log: >> Added fix for MATH-597: Implemented faster generation of random exponential >> distributed values with algorithm from Ahrens and Dieter (1972): Computer >> methods for sampling from the exponential and normal distributions. Test >> case was improved, too. >> >> 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=1137795&r1=1137794&r2=1137795&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 Jun 20 21:42:48 2011 >> @@ -44,6 +44,7 @@ import org.apache.commons.math.exception >> import org.apache.commons.math.exception.util.LocalizedFormats; >> import org.apache.commons.math.util.FastMath; >> import org.apache.commons.math.util.MathUtils; >> +import org.apache.commons.math.util.ResizableDoubleArray; >> >> /** >> * Implements the {@link RandomData} interface using a {@link >> RandomGenerator} >> @@ -107,6 +108,21 @@ public class RandomDataImpl implements R >> /** Serializable version identifier */ >> private static final long serialVersionUID = -626730818244969716L; >> >> + /** Used when generating Exponential samples >> + * [1] writes: >> + * One table containing the constants >> + * q_i = sum_{j=1}^i (ln 2)^j/j! = ln 2 + (ln 2)^2/2 + ... + (ln 2)^i/i! >> + * until the largest representable fraction below 1 is exceeded. >> + * >> + * Note that >> + * 1 = 2 - 1 = exp(ln 2) - 1 = sum_{n=1}^infty (ln 2)^n / n! >> + * thus q_i -> 1 as i -> infty, >> + * so the higher 1, the closer to one we get (the series is not >> alternating). >> + * >> + * By trying, n = 16 in Java is enough to reach 1.0. >> + */ >> + private static double[] EXPONENTIAL_SA_QI = null; >> + >> /** underlying random number generator */ >> private RandomGenerator rand = null; >> >> @@ -114,6 +130,35 @@ public class RandomDataImpl implements R >> private SecureRandom secRand = null; >> >> /** >> + * Initialize tables >> + */ >> + static { >> + /** >> + * Filling EXPONENTIAL_SA_QI table. >> + * Note that we don't want qi = 0 in the table. >> + */ >> + final double LN2 = FastMath.log(2); >> + double qi = 0; >> + int i = 1; >> + >> + /** >> + * MathUtils provides factorials up to 20, so let's use that limit >> together >> + * with MathUtils.EPSILON to generate the following code (a priori, >> we know that >> + * there will be 16 elements, but instead of hardcoding that, this >> is >> + * prettier): >> + */ >> + final ResizableDoubleArray ra = new ResizableDoubleArray(20); >> + >> + while (qi < 1) { >> + qi += FastMath.pow(LN2, i) / MathUtils.factorial(i); >> + ra.addElement(qi); >> + ++i; >> + } >> + >> + EXPONENTIAL_SA_QI = ra.getElements(); >> + } >> + >> + /** >> * Construct a RandomDataImpl. >> */ >> public RandomDataImpl() { >> @@ -469,10 +514,11 @@ public class RandomDataImpl implements R >> * Returns a random value from an Exponential distribution with the given >> * mean. >> * <p> >> - * <strong>Algorithm Description</strong>: Uses the <a >> - * href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html";> >> Inversion >> - * Method</a> to generate exponentially distributed random values from >> - * uniform deviates. >> + * <strong>Algorithm Description</strong>: Uses the Algorithm SA >> (Ahrens) >> + * from p. 876 in: >> + * [1]: Ahrens, J. H. and Dieter, U. (1972). Computer methods for >> + * sampling from the exponential and normal distributions. >> + * Communications of the ACM, 15, 873-882. >> * </p> >> * >> * @param mean the mean of the distribution >> @@ -483,12 +529,43 @@ public class RandomDataImpl implements R >> if (mean <= 0.0) { >> throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, >> mean); >> } >> - final RandomGenerator generator = getRan(); >> - double unif = generator.nextDouble(); >> - while (unif == 0.0d) { >> - unif = generator.nextDouble(); >> + >> + // Step 1: >> + double a = 0; >> + double u = this.nextUniform(0, 1); >> + >> + // Step 2 and 3: >> + while (u < 0.5) { >> + a += EXPONENTIAL_SA_QI[0]; >> + u *= 2; >> } >> - return -mean * FastMath.log(unif); >> + >> + // Step 4 (now u >= 0.5): >> + u += u - 1; >> + >> + // Step 5: >> + if (u <= EXPONENTIAL_SA_QI[0]) { >> + return mean * (a + u); >> + } >> + >> + // Step 6: >> + int i = 0; // Should be 1, be we iterate before it in while using 0 >> + double u2 = this.nextUniform(0, 1); >> + double umin = u2; >> + >> + // Step 7 and 8: >> + do { >> + ++i; >> + u2 = this.nextUniform(0, 1); >> + >> + if (u2 < umin) { >> + umin = u2; >> + } >> + >> + // Step 8: >> + } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since >> EXPONENTIAL_SA_QI[MAX] = 1 >> + >> + return mean * (a + umin * EXPONENTIAL_SA_QI[0]); >> } >> >> /** >> >> 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=1137795&r1=1137794&r2=1137795&view=diff >> ============================================================================== >> --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) >> +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Jun 20 21:42:48 >> 2011 >> @@ -52,6 +52,11 @@ The <action> type attribute can be add,u >> If the output is not quite correct, check for invisible trailing spaces! >> --> >> <release version="3.0" date="TBD" description="TBD"> >> + <action dev="mikl" type="fix" issue="MATH-597"> >> + Implemented faster generation of random exponential distributed >> values with >> + algorithm from Ahrens and Dieter (1972): Computer methods for >> sampling >> + from the exponential and normal distributions. >> + </action> >> <action dev="luc" type="add" issue="MATH-548"> >> K-means++ clustering can now run multiple trials >> </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=1137795&r1=1137794&r2=1137795&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 Jun 20 21:42:48 2011 >> @@ -30,6 +30,7 @@ import org.apache.commons.math.distribut >> import org.apache.commons.math.distribution.BinomialDistributionTest; >> import org.apache.commons.math.distribution.CauchyDistributionImpl; >> import org.apache.commons.math.distribution.ChiSquaredDistributionImpl; >> +import org.apache.commons.math.distribution.ExponentialDistributionImpl; >> import org.apache.commons.math.distribution.FDistributionImpl; >> import org.apache.commons.math.distribution.GammaDistributionImpl; >> import org.apache.commons.math.distribution.HypergeometricDistributionImpl; >> @@ -245,10 +246,10 @@ public class RandomDataTest { >> >> @Test >> public void testNextPoissonConsistency() throws Exception { >> - >> + >> // Reseed randomGenerator to get fixed sequence >> - randomData.reSeed(1000); >> - >> + randomData.reSeed(1000); >> + >> // Small integral means >> for (int i = 1; i < 100; i++) { >> checkNextPoissonConsistency(i); >> @@ -581,7 +582,7 @@ public class RandomDataTest { >> >> /** test failure modes and distribution of nextExponential() */ >> @Test >> - public void testNextExponential() { >> + public void testNextExponential() throws Exception { >> try { >> randomData.nextExponential(-1); >> Assert.fail("negative mean -- expecting >> MathIllegalArgumentException"); >> @@ -609,6 +610,32 @@ public class RandomDataTest { >> */ >> Assert.assertEquals("exponential cumulative distribution", (double) >> cumFreq >> / (double) largeSampleSize, 0.8646647167633873, .2); >> + >> + /** >> + * Proposal on improving the test of generating exponentials >> + */ >> + double[] quartiles; >> + long[] counts; >> + >> + // Mean 1 >> + quartiles = TestUtils.getDistributionQuartiles(new >> ExponentialDistributionImpl(1)); >> + counts = new long[4]; >> + randomData.reSeed(1000); >> + for (int i = 0; i < 1000; i++) { >> + double value = randomData.nextExponential(1); >> + TestUtils.updateCounts(value, counts, quartiles); >> + } >> + TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> + >> + // Mean 5 >> + quartiles = TestUtils.getDistributionQuartiles(new >> ExponentialDistributionImpl(5)); >> + counts = new long[4]; >> + randomData.reSeed(1000); >> + for (int i = 0; i < 1000; i++) { >> + double value = randomData.nextExponential(5); >> + TestUtils.updateCounts(value, counts, quartiles); >> + } >> + TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> >> /** test reseeding, algorithm/provider games */ >> @@ -810,7 +837,7 @@ public class RandomDataTest { >> Assert.fail("permutation not found"); >> return -1; >> } >> - >> + >> @Test >> public void testNextInversionDeviate() throws Exception { >> // Set the seed for the default random generator >> @@ -830,9 +857,9 @@ public class RandomDataTest { >> for (int i = 0; i < 10; i++) { >> double value = randomData.nextInversionDeviate(betaDistribution); >> >> Assert.assertEquals(betaDistribution.cumulativeProbability(value), >> quantiles[i], 10E-9); >> - } >> + } >> } >> - >> + >> @Test >> public void testNextBeta() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> BetaDistributionImpl(2,5)); >> @@ -844,7 +871,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextCauchy() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> CauchyDistributionImpl(1.2, 2.1)); >> @@ -856,7 +883,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextChiSquare() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> ChiSquaredDistributionImpl(12)); >> @@ -868,7 +895,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextF() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> FDistributionImpl(12, 5)); >> @@ -880,7 +907,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextGamma() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> GammaDistributionImpl(4, 2)); >> @@ -892,7 +919,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextT() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> TDistributionImpl(10)); >> @@ -904,7 +931,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextWeibull() throws Exception { >> double[] quartiles = TestUtils.getDistributionQuartiles(new >> WeibullDistributionImpl(1.2, 2.1)); >> @@ -916,7 +943,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(expected, counts, 0.001); >> } >> - >> + >> @Test >> public void testNextBinomial() throws Exception { >> BinomialDistributionTest testInstance = new >> BinomialDistributionTest(); >> @@ -942,7 +969,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, >> observedCounts, .001); >> } >> - >> + >> @Test >> public void testNextHypergeometric() throws Exception { >> HypergeometricDistributionTest testInstance = new >> HypergeometricDistributionTest(); >> @@ -968,7 +995,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, >> observedCounts, .001); >> } >> - >> + >> @Test >> public void testNextPascal() throws Exception { >> PascalDistributionTest testInstance = new PascalDistributionTest(); >> @@ -993,7 +1020,7 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, >> observedCounts, .001); >> } >> - >> + >> @Test >> public void testNextZipf() throws Exception { >> ZipfDistributionTest testInstance = new ZipfDistributionTest(); >> @@ -1018,5 +1045,5 @@ public class RandomDataTest { >> } >> TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, >> observedCounts, .001); >> } >> - >> + >> } >> >> >> > --------------------------------------------------------------------- > To unsubscribe, e-mail: [email protected] > For additional commands, e-mail: [email protected] > > --------------------------------------------------------------------- To unsubscribe, e-mail: [email protected] For additional commands, e-mail: [email protected]
