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 df423bf RNG-152: Update samplers to use ZigguratSampler.NormalizedGaussian df423bf is described below commit df423bf1bbd2679177095915242e724a1e625b4a Author: aherbert <aherb...@apache.org> AuthorDate: Thu Jul 8 11:29:40 2021 +0100 RNG-152: Update samplers to use ZigguratSampler.NormalizedGaussian --- .../jmh/sampling/UnitSphereSamplerBenchmark.java | 10 ++--- .../sampling/shape/UnitBallSamplerBenchmark.java | 17 ++++---- .../commons/rng/sampling/UnitSphereSampler.java | 8 ++-- .../AhrensDieterMarsagliaTsangGammaSampler.java | 4 +- .../distribution/LargeMeanPoissonSampler.java | 4 +- .../rng/sampling/distribution/LevySampler.java | 4 +- .../rng/sampling/shape/UnitBallSampler.java | 5 +-- .../sampling/distribution/GaussianSamplerTest.java | 4 +- .../rng/sampling/distribution/LevySamplerTest.java | 50 +++++++++++++--------- .../distribution/LogNormalSamplerTest.java | 6 +-- src/changes/changes.xml | 3 ++ 11 files changed, 62 insertions(+), 53 deletions(-) diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java index 2911155..fce9d31 100644 --- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java +++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/UnitSphereSamplerBenchmark.java @@ -19,7 +19,7 @@ package org.apache.commons.rng.examples.jmh.sampling; import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; -import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; +import org.apache.commons.rng.sampling.distribution.ZigguratSampler; import org.apache.commons.rng.simple.RandomSource; import org.openjdk.jmh.annotations.Benchmark; import org.openjdk.jmh.annotations.BenchmarkMode; @@ -233,7 +233,7 @@ public class UnitSphereSamplerBenchmark { * @param rng the source of randomness */ UnitSphereSampler2D(UniformRandomProvider rng) { - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -291,7 +291,7 @@ public class UnitSphereSamplerBenchmark { * @param rng the source of randomness */ UnitSphereSampler3D(UniformRandomProvider rng) { - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -350,7 +350,7 @@ public class UnitSphereSamplerBenchmark { * @param rng the source of randomness */ UnitSphereSampler4D(UniformRandomProvider rng) { - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -387,7 +387,7 @@ public class UnitSphereSamplerBenchmark { */ ArrayBasedUnitSphereSampler(int dimension, UniformRandomProvider rng) { this.dimension = dimension; - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/UnitBallSamplerBenchmark.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/UnitBallSamplerBenchmark.java index f4523e0..04c098d 100644 --- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/UnitBallSamplerBenchmark.java +++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/shape/UnitBallSamplerBenchmark.java @@ -20,7 +20,6 @@ package org.apache.commons.rng.examples.jmh.sampling.shape; import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.sampling.distribution.ContinuousSampler; import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; -import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; import org.apache.commons.rng.sampling.distribution.ZigguratSampler; import org.apache.commons.rng.simple.RandomSource; import org.openjdk.jmh.annotations.Benchmark; @@ -307,7 +306,7 @@ public class UnitBallSamplerBenchmark { */ HypersphereInternalSampler(UniformRandomProvider rng) { super(rng); - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -340,7 +339,7 @@ public class UnitBallSamplerBenchmark { */ HypersphereDiscardSampler(UniformRandomProvider rng) { super(rng); - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -434,7 +433,7 @@ public class UnitBallSamplerBenchmark { */ BallPointSampler(UniformRandomProvider rng) { super(rng); - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2) // thus is the equivalent of the HypersphereDiscardSampler. // Here we use mean = 1 and scale the output later. @@ -472,7 +471,7 @@ public class UnitBallSamplerBenchmark { */ HypersphereInternalSampler(UniformRandomProvider rng) { super(rng); - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -506,7 +505,7 @@ public class UnitBallSamplerBenchmark { */ HypersphereDiscardSampler(UniformRandomProvider rng) { super(rng); - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -619,7 +618,7 @@ public class UnitBallSamplerBenchmark { BallPointSampler(UniformRandomProvider rng, int dimension) { super(rng); this.dimension = dimension; - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); // Exponential(mean=2) == Chi-squared distribution(degrees freedom=2) // thus is the equivalent of the HypersphereDiscardSampler. // Here we use mean = 1 and scale the output later. @@ -670,7 +669,7 @@ public class UnitBallSamplerBenchmark { super(rng); this.dimension = dimension; power = 1.0 / dimension; - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -714,7 +713,7 @@ public class UnitBallSamplerBenchmark { HypersphereDiscardSampler(UniformRandomProvider rng, int dimension) { super(rng); this.dimension = dimension; - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); } @Override diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java index 76db7b0..53c81f9 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/UnitSphereSampler.java @@ -19,7 +19,7 @@ package org.apache.commons.rng.sampling; import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; -import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; +import org.apache.commons.rng.sampling.distribution.ZigguratSampler; /** * Generate vectors <a href="http://mathworld.wolfram.com/SpherePointPicking.html"> @@ -95,7 +95,7 @@ public class UnitSphereSampler implements SharedStateObjectSampler<double[]> { * @param rng Source of randomness. */ UnitSphereSampler2D(UniformRandomProvider rng) { - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -131,7 +131,7 @@ public class UnitSphereSampler implements SharedStateObjectSampler<double[]> { * @param rng Source of randomness. */ UnitSphereSampler3D(UniformRandomProvider rng) { - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override @@ -171,7 +171,7 @@ public class UnitSphereSampler implements SharedStateObjectSampler<double[]> { */ UnitSphereSamplerND(int dimension, UniformRandomProvider rng) { this.dimension = dimension; - sampler = new ZigguratNormalizedGaussianSampler(rng); + sampler = ZigguratSampler.NormalizedGaussian.of(rng); } @Override diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java index 08c72b8..c3c884d 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java @@ -219,7 +219,7 @@ public class AhrensDieterMarsagliaTsangGammaSampler double alpha, double theta) { super(rng, alpha, theta); - gaussian = new ZigguratNormalizedGaussianSampler(rng); + gaussian = ZigguratSampler.NormalizedGaussian.of(rng); dOptim = alpha - ONE_THIRD; cOptim = ONE_THIRD / Math.sqrt(dOptim); } @@ -231,7 +231,7 @@ public class AhrensDieterMarsagliaTsangGammaSampler MarsagliaTsangGammaSampler(UniformRandomProvider rng, MarsagliaTsangGammaSampler source) { super(rng, source); - gaussian = new ZigguratNormalizedGaussianSampler(rng); + gaussian = ZigguratSampler.NormalizedGaussian.of(rng); dOptim = source.dOptim; cOptim = source.cOptim; } 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 41f9eca..d56974f 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 @@ -132,7 +132,7 @@ public class LargeMeanPoissonSampler } this.rng = rng; - gaussian = new ZigguratNormalizedGaussianSampler(rng); + gaussian = ZigguratSampler.NormalizedGaussian.of(rng); exponential = ZigguratSampler.Exponential.of(rng); // Plain constructor uses the uncached function. factorialLog = NO_CACHE_FACTORIAL_LOG; @@ -177,7 +177,7 @@ public class LargeMeanPoissonSampler } this.rng = rng; - gaussian = new ZigguratNormalizedGaussianSampler(rng); + gaussian = ZigguratSampler.NormalizedGaussian.of(rng); exponential = ZigguratSampler.Exponential.of(rng); // Plain constructor uses the uncached function. factorialLog = NO_CACHE_FACTORIAL_LOG; diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java index f604604..5207755 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LevySampler.java @@ -42,7 +42,7 @@ public final class LevySampler implements SharedStateContinuousSampler { private LevySampler(UniformRandomProvider rng, double location, double scale) { - this.gaussian = new ZigguratNormalizedGaussianSampler(rng); + this.gaussian = ZigguratSampler.NormalizedGaussian.of(rng); this.location = location; this.scale = scale; this.rng = rng; @@ -54,7 +54,7 @@ public final class LevySampler implements SharedStateContinuousSampler { */ private LevySampler(UniformRandomProvider rng, LevySampler source) { - this.gaussian = new ZigguratNormalizedGaussianSampler(rng); + this.gaussian = ZigguratSampler.NormalizedGaussian.of(rng); this.location = source.location; this.scale = source.scale; this.rng = rng; diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/UnitBallSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/UnitBallSampler.java index 80c836c..8e78c0f 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/UnitBallSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/shape/UnitBallSampler.java @@ -21,7 +21,6 @@ import org.apache.commons.rng.UniformRandomProvider; import org.apache.commons.rng.sampling.SharedStateObjectSampler; import org.apache.commons.rng.sampling.distribution.ContinuousSampler; import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler; -import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; import org.apache.commons.rng.sampling.distribution.ZigguratSampler; /** @@ -124,7 +123,7 @@ public abstract class UnitBallSampler implements SharedStateObjectSampler<double * @param rng Source of randomness. */ UnitBallSampler3D(UniformRandomProvider rng) { - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); // Require an Exponential(mean=2). // Here we use mean = 1 and scale the output later. exp = ZigguratSampler.Exponential.of(rng); @@ -169,7 +168,7 @@ public abstract class UnitBallSampler implements SharedStateObjectSampler<double */ UnitBallSamplerND(int dimension, UniformRandomProvider rng) { this.dimension = dimension; - normal = new ZigguratNormalizedGaussianSampler(rng); + normal = ZigguratSampler.NormalizedGaussian.of(rng); // Require an Exponential(mean=2). // Here we use mean = 1 and scale the output later. exp = ZigguratSampler.Exponential.of(rng); diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GaussianSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GaussianSamplerTest.java index 2ad15f8..9222f76 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GaussianSamplerTest.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GaussianSamplerTest.java @@ -34,7 +34,7 @@ public class GaussianSamplerTest { public void testConstructorThrowsWithZeroStandardDeviation() { final RestorableUniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(0L); - final NormalizedGaussianSampler gauss = new ZigguratNormalizedGaussianSampler(rng); + final NormalizedGaussianSampler gauss = ZigguratSampler.NormalizedGaussian.of(rng); final double mean = 1; final double standardDeviation = 0; GaussianSampler.of(gauss, mean, standardDeviation); @@ -99,7 +99,7 @@ public class GaussianSamplerTest { public void testSharedStateSampler() { final UniformRandomProvider rng1 = RandomSource.SPLIT_MIX_64.create(0L); final UniformRandomProvider rng2 = RandomSource.SPLIT_MIX_64.create(0L); - final NormalizedGaussianSampler gauss = new ZigguratNormalizedGaussianSampler(rng1); + final NormalizedGaussianSampler gauss = ZigguratSampler.NormalizedGaussian.of(rng1); final double mean = 1.23; final double standardDeviation = 4.56; final SharedStateContinuousSampler sampler1 = diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java index 4678f36..7430f38 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LevySamplerTest.java @@ -71,7 +71,7 @@ public class LevySamplerTest { public void testSupport() { final double location = 0.0; final double scale = 1.0; - // Force the underlying ZigguratNormalizedGaussianSampler to create 0 + // Force the underlying ZigguratSampler.NormalizedGaussian to create 0 final LevySampler s1 = LevySampler.of( new SplitMix64(0L) { @Override @@ -81,34 +81,42 @@ public class LevySamplerTest { }, location, scale); Assert.assertEquals(Double.POSITIVE_INFINITY, s1.sample(), 0.0); - // Force the underlying ZigguratNormalizedGaussianSampler to create the largest value. - // This is 14.11 + // Force the underlying ZigguratSampler.NormalizedGaussian to create a large + // sample in the tail of the distribution. + // The first three -1,-1,-1 values enters the tail of the distribution. + // Here an exponential is added to 3.6360066255. + // The exponential also requires -1,-1 to recurse. Each recursion adds 7.56927469415 + // to the exponential. A value of 0 stops recursion with a sample of 0. + // Two exponentials are required: x and y. + // The exponential is multiplied by 0.275027001597525 to create x. + // The condition 2y >= x^x must be true to return x. + // Create x = 4 * 7.57 and y = 16 * 7.57 + final long[] sequence = { + // Sample the Gaussian tail + -1, -1, -1, + // Exponential x = 4 * 7.57... * 0.275027001597525 + -1, -1, -1, -1, -1, -1, -1, -1, 0, + // Exponential y = 16 * 7.57... + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, + }; final LevySampler s2 = LevySampler.of( new SplitMix64(0L) { private int i; @Override public long next() { - i++; - if (i == 1) { - // Set the first value to ensure we sample the tail of the ziggurat. - // The lowest 7 bits are zero to select rectangle 0 from the ziggurat. - return (Long.MAX_VALUE << 7) & Long.MAX_VALUE; + if (i++ < sequence.length) { + return sequence[i - 1]; } - if (i == 2) { - // Set the second value to generate y as the largest value possible by - // ensuring Math.log is called with a small value. - return 0L; - } - // The next value generates x which must be set to the largest value x which - // satisfies the condition: - // 2y >= x^2 - return 1377L << 11; + return super.next(); } }, location, scale); - // The tail of the zigguart should be s=12.014118700751192 - // expected is 1/s^2 = 0.006928132149804786 - // This is as close to zero as the sampler can get. - final double expected = 0.006928132149804786; + // The tail of the zigguart should be approximately s=11.963 + final double s = 4 * 7.56927469415 * 0.275027001597525 + 3.6360066255; + // expected is 1/s^2 = 0.006987 + // So the sampler never achieves the lower bound of zero. + // It requires an extreme deviate from the Gaussian. + final double expected = 1 / (s * s); Assert.assertEquals(expected, s2.sample(), 0.0); } } diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LogNormalSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LogNormalSamplerTest.java index 03554ab..a554fd1 100644 --- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LogNormalSamplerTest.java +++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/LogNormalSamplerTest.java @@ -34,7 +34,7 @@ public class LogNormalSamplerTest { public void testConstructorThrowsWithNegativeScale() { final RestorableUniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(0L); - final NormalizedGaussianSampler gauss = new ZigguratNormalizedGaussianSampler(rng); + final NormalizedGaussianSampler gauss = ZigguratSampler.NormalizedGaussian.of(rng); final double scale = -1e-6; final double shape = 1; LogNormalSampler.of(gauss, scale, shape); @@ -47,7 +47,7 @@ public class LogNormalSamplerTest { public void testConstructorThrowsWithZeroShape() { final RestorableUniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(0L); - final NormalizedGaussianSampler gauss = new ZigguratNormalizedGaussianSampler(rng); + final NormalizedGaussianSampler gauss = ZigguratSampler.NormalizedGaussian.of(rng); final double scale = 1; final double shape = 0; LogNormalSampler.of(gauss, scale, shape); @@ -60,7 +60,7 @@ public class LogNormalSamplerTest { public void testSharedStateSampler() { final UniformRandomProvider rng1 = RandomSource.SPLIT_MIX_64.create(0L); final UniformRandomProvider rng2 = RandomSource.SPLIT_MIX_64.create(0L); - final NormalizedGaussianSampler gauss = new ZigguratNormalizedGaussianSampler(rng1); + final NormalizedGaussianSampler gauss = ZigguratSampler.NormalizedGaussian.of(rng1); final double scale = 1.23; final double shape = 4.56; final SharedStateContinuousSampler sampler1 = diff --git a/src/changes/changes.xml b/src/changes/changes.xml index c71270f..38970b7 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -77,6 +77,9 @@ re-run tests that fail, and pass the build if they succeed within the allotted number of reruns (the test will be marked as 'flaky' in the report). "> + <action dev="aherbert" type="update" issue="152"> + Update samplers to use ZigguratSampler.NormalizedGaussian for Gaussian deviates. + </action> <action dev="aherbert" type="fix" issue="146"> "GaussianSampler": Prevent infinite mean and standard deviation. </action>