You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/08/10 12:08:47 UTC
[commons-rng] 02/03: Additional tests for ziggurat exponential
sampler
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
commit 68388c23255ae80ddca0e3361fd952cde72323ce
Author: aherbert <ah...@apache.org>
AuthorDate: Tue Aug 10 10:16:44 2021 +0100
Additional tests for ziggurat exponential sampler
---
.../sampling/distribution/ZigguratSamplerTest.java | 128 +++++++++++++++++++--
1 file changed, 117 insertions(+), 11 deletions(-)
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZigguratSamplerTest.java b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZigguratSamplerTest.java
index 3a862c2..4ac616e 100644
--- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZigguratSamplerTest.java
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZigguratSamplerTest.java
@@ -19,10 +19,13 @@ package org.apache.commons.rng.sampling.distribution;
import org.junit.Assert;
import org.junit.Test;
import java.util.Arrays;
+import org.apache.commons.math3.distribution.AbstractRealDistribution;
+import org.apache.commons.math3.distribution.ExponentialDistribution;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.stat.inference.ChiSquareTest;
import org.apache.commons.rng.RestorableUniformRandomProvider;
import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.core.source64.SplitMix64;
import org.apache.commons.rng.sampling.RandomAssert;
import org.apache.commons.rng.simple.RandomSource;
@@ -78,6 +81,65 @@ public class ZigguratSamplerTest {
}
/**
+ * Test the recursion in the exponential distribution.
+ */
+ @Test
+ public void testExponentialRecursion() {
+ // The exponential distribution will enter the edge of the ziggurat if the RNG
+ // outputs -1 (all bits). This performs alias sampling using a long value.
+ // The tail will be selected if the next output is -1.
+ // Thus two -1 values enter recursion where a new exponential sample is added
+ // to the tail value:
+ final double tailValue = 7.569274694148063;
+
+ // Alias sampling assigns the ziggurat layer using the lower 8 bits.
+ // The rest determine if the layer or the alias are used. We do not control this
+ // and leave it to a seeded RNG to select different layers.
+
+ // A value of zero will create a sample of zero (42 is the seed)
+ Assert.assertEquals(0.0, expSample(42, 0), 0.0);
+ Assert.assertEquals(tailValue, expSample(42, -1, -1, 0), 0.0);
+
+ // Use different seeds to test different layers from the edge of the ziggurat.
+ for (final long seed : new long[] {42, -2136612838, 2340923842L, -1263746817818681L}) {
+ // Base value
+ final double x0 = expSample(seed);
+ // Edge value
+ final double x1 = expSample(seed, -1);
+ // Recursion
+ // Note the order of additions is important as the final sample is added to the
+ // tail value multiple times.
+ Assert.assertEquals(x0 + tailValue, expSample(seed, -1, -1), 0.0);
+ Assert.assertEquals(x1 + tailValue, expSample(seed, -1, -1, -1), 0.0);
+ // Double recursion
+ Assert.assertEquals(x0 + tailValue + tailValue, expSample(seed, -1, -1, -1, -1), 0.0);
+ Assert.assertEquals(x1 + tailValue + tailValue, expSample(seed, -1, -1, -1, -1, -1), 0.0);
+ }
+ }
+
+ /**
+ * Create an exponential sample from the sequence of longs, then revert to a seed RNG.
+ *
+ * @param seed the seed
+ * @param longs the longs
+ * @return the sample
+ */
+ private static double expSample(long seed, final long... longs) {
+ final SplitMix64 rng = new SplitMix64(seed) {
+ private int i;
+ @Override
+ public long next() {
+ if (i == longs.length) {
+ // Revert to seeded RNG
+ return super.next();
+ }
+ return longs[i++];
+ }
+ };
+ return ZigguratSampler.Exponential.of(rng).sample();
+ }
+
+ /**
* Test Gaussian samples using a large number of bins based on uniformly spaced quantiles.
* Added for RNG-159.
*/
@@ -89,7 +151,7 @@ public class ZigguratSamplerTest {
for (int i = 0; i < bins; i++) {
quantiles[i] = dist.inverseCumulativeProbability((i + 1.0) / bins);
}
- testGaussianSamples(quantiles);
+ testSamples(quantiles, false);
}
/**
@@ -105,23 +167,60 @@ public class ZigguratSamplerTest {
for (int i = 0; i < bins; i++) {
values[i] = minx + (maxx - minx) * (i + 1.0) / bins;
}
- testGaussianSamples(values);
+ // Ensure upper bound is the support limit
+ values[bins - 1] = Double.POSITIVE_INFINITY;
+ testSamples(values, false);
+ }
+
+ /**
+ * Test exponential samples using a large number of bins based on uniformly spaced quantiles.
+ */
+ @Test
+ public void testExponentialSamplesWithQuantiles() {
+ final int bins = 2000;
+ final ExponentialDistribution dist = new ExponentialDistribution(null, 1.0);
+ final double[] quantiles = new double[bins];
+ for (int i = 0; i < bins; i++) {
+ quantiles[i] = dist.inverseCumulativeProbability((i + 1.0) / bins);
+ }
+ testSamples(quantiles, true);
+ }
+
+ /**
+ * Test exponential samples using a large number of bins uniformly spaced in a range.
+ */
+ @Test
+ public void testExponentialSamplesWithUniformValues() {
+ final int bins = 2000;
+ final double[] values = new double[bins];
+ final double minx = 0;
+ // Enter the tail of the distribution
+ final double maxx = 12;
+ for (int i = 0; i < bins; i++) {
+ values[i] = minx + (maxx - minx) * (i + 1.0) / bins;
+ }
+ // Ensure upper bound is the support limit
+ values[bins - 1] = Double.POSITIVE_INFINITY;
+ testSamples(values, true);
}
/**
- * Test Gaussian samples using the provided bins. Values correspond to the bin
- * upper limit. It is assumed the values span the mean, and ideally most of the
- * normal distribution.
+ * Test samples using the provided bins. Values correspond to the bin upper limit. It
+ * is assumed the values span most of the distribution. Additional tests are performed
+ * using a region of the distribution sampled using the edge of the ziggurat.
*
* @param values Bin upper limits
+ * @param exponential Set the true to use an exponential sampler
*/
- private static void testGaussianSamples(double[] values) {
+ private static void testSamples(double[] values,
+ boolean exponential) {
final int bins = values.length;
final int samples = 10000000;
final long[] observed = new long[bins];
final RestorableUniformRandomProvider rng = RandomSource.XO_SHI_RO_128_PP.create(0xabcdefL);
- final ZigguratSampler.NormalizedGaussian sampler = ZigguratSampler.NormalizedGaussian.of(rng);
+ final ContinuousSampler sampler = exponential ?
+ ZigguratSampler.Exponential.of(rng) : ZigguratSampler.NormalizedGaussian.of(rng);
for (int i = 0; i < samples; i++) {
final double x = sampler.sample();
final int index = findIndex(values, x);
@@ -129,7 +228,8 @@ public class ZigguratSamplerTest {
}
// Compute expected
- final NormalDistribution dist = new NormalDistribution(null, 0.0, 1.0);
+ final AbstractRealDistribution dist = exponential ?
+ new ExponentialDistribution(null, 1.0) : new NormalDistribution(null, 0.0, 1.0);
final double[] expected = new double[bins];
double x0 = Double.NEGATIVE_INFINITY;
for (int i = 0; i < bins; i++) {
@@ -140,13 +240,15 @@ public class ZigguratSamplerTest {
final double significanceLevel = 0.001;
+ final double lowerBound = dist.getSupportLowerBound();
+
final ChiSquareTest chiSquareTest = new ChiSquareTest();
// Pass if we cannot reject null hypothesis that the distributions are the same.
double chi2 = chiSquareTest.chiSquareTest(expected, observed);
Assert.assertFalse(String.format("(%s <= x < %s) Chi-square p-value = %s",
- Double.NEGATIVE_INFINITY, values[bins - 1], chi2), chi2 < 0.001);
+ lowerBound, values[bins - 1], chi2), chi2 < 0.001);
- // Test around the mean
+ // Test bins sampled from the edge of the ziggurat. This is always around zero.
for (final double range : new double[] {0.5, 0.25, 0.1, 0.05}) {
final int min = findIndex(values, -range);
final int max = findIndex(values, range);
@@ -154,7 +256,7 @@ public class ZigguratSamplerTest {
final double[] expected2 = Arrays.copyOfRange(expected, min, max + 1);
chi2 = chiSquareTest.chiSquareTest(expected2, observed2);
Assert.assertFalse(String.format("(%s <= x < %s) Chi-square p-value = %s",
- min == 0 ? Double.NEGATIVE_INFINITY : values[min - 1], values[max], chi2),
+ min == 0 ? lowerBound : values[min - 1], values[max], chi2),
chi2 < significanceLevel);
}
}
@@ -165,6 +267,10 @@ public class ZigguratSamplerTest {
* data[index - 1] <= x < data[index]
* </pre>
*
+ * <p>This is a specialised binary search that assumes the bounds of the data are the
+ * extremes of the support, and the upper support is infinite. Thus an index cannot
+ * be returned as equal to the data length.
+ *
* @param data the data
* @param x the value
* @return the index