commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <phil.ste...@gmail.com>
Subject Re: [MATH] Commit review of MATH-597
Date Mon, 20 Jun 2011 22:12:33 GMT
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  <mikl@apache.org>:
>> 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: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message