This is an automated email from the ASF dual-hosted git repository.
asf-gitbox-commits 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 05b4f366 RNG-203: Add a zeta distribution sampler
05b4f366 is described below
commit 05b4f3661349b0de1559a79d822cefef460d3252
Author: Alex Herbert <[email protected]>
AuthorDate: Fri Sep 18 12:51:49 2026 +0100
RNG-203: Add a zeta distribution sampler
---
.../rng/sampling/distribution/ZetaSampler.java | 222 +++++++++++++++++++++
.../DiscreteSamplerParametricTest.java | 23 ++-
.../distribution/DiscreteSamplerTestData.java | 55 ++++-
.../distribution/DiscreteSamplersList.java | 72 ++++++-
.../rng/sampling/distribution/ZetaSamplerTest.java | 180 +++++++++++++++++
src/changes/changes.xml | 3 +
6 files changed, 537 insertions(+), 18 deletions(-)
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZetaSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZetaSampler.java
new file mode 100644
index 00000000..8e30c8cc
--- /dev/null
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZetaSampler.java
@@ -0,0 +1,222 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.rng.sampling.distribution;
+
+import java.util.function.LongToDoubleFunction;
+import org.apache.commons.rng.UniformRandomProvider;
+
+/**
+ * Sampling from a zeta distribution.
+ *
+ * <p>Note that the zeta distribution has an upper limit of positive infinity.
This
+ * implementation is clipped to {@link Integer#MAX_VALUE}. As the exponent \(
s \to 1 \)
+ * sampling will become biased to 2<sup>31</sup> - 1 as the distribution
density is
+ * truncated by the {@code int} limit. When the exponent is large the density
+ * of the zeta distribution is concentrated at 1. In {@code double} precision
this
+ * limits the rejection method to {@code s <= 54}; larger {@code s} will
always sample 1.
+ *
+ * <p>Rejection sampling method for a zeta distribution adapted from:
+ * <blockquote>
+ * Luc Devroye (1986)
+ * <i>"Non-uniform random variate generation",</i><br>
+ * <strong>Springer New York, NY</strong> pp 550-552.
+ * </blockquote>
+ *
+ * <p>Sampling uses {@link UniformRandomProvider#nextLong()}.</p>
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Zeta_distribution";>Zeta
distribution (Wikipedia)</a>
+ * @since 1.4
+ */
+public final class ZetaSampler {
+ /** The minimum exponent. The exponent must be above 1. */
+ private static final double MIN_EXPONENT = Math.nextUp(1.0);
+ /**
+ * The maximum exponent for the rejection sampler.
+ * Note that if the sampler uses u in (0, 1] to bias towards 1 then the
sample
+ * x = floor ( U^{-1/(a-1)} ) is always 1 when a is large.
+ * The threshold for a is 54 if u=2^-53.
+ * Note the zeta distribution CDF(x=1; s=54) = 1.0 in {@code double}
precision.
+ */
+ private static final double REJECTION_EXPONENT = 54;
+
+ /**
+ * Sample from the zeta distribution using a rejection method.
+ * Package-private for testing.
+ */
+ static final class RejectionZetaSampler implements
SharedStateDiscreteSampler {
+ /**
+ * The threshold to bias the extreme sample to 1 or infinity. Change
the
+ * extreme sample of the zeta distribution using the midpoint of the
support
+ * domain, i.e. x = 2^31 / 2; cdf(x; a) = sf(x; a) ~ 0.5.
+ */
+ private static final double THRESHOLD = 1.0324376395045163;
+ /** ln(2). */
+ private static final double LN2 = Math.log(2.0);
+
+ /** Source of randomness. */
+ private final UniformRandomProvider rng;
+ /** a - 1. */
+ private final double am1;
+ /** Reciprocal of (a - 1) = 1 / (a - 1). */
+ private final double ram1;
+ /** 2^(a-1) / (2^(a-1) - 1). */
+ private final double bObm1;
+ /** Function to compute u in [0, 1]. */
+ private final LongToDoubleFunction nextU;
+
+ /**
+ * Create an instance.
+ *
+ * @param rng Source of randomness.
+ * @param a Exponent of the zeta distribution ({@code a > 1}).
+ */
+ RejectionZetaSampler(UniformRandomProvider rng, double a) {
+ this.rng = rng;
+ am1 = a - 1;
+ ram1 = 1 / am1;
+ // b = 2^(a-1)
+ // This will not overflow within the usable range of the algorithm.
+ // We never expect infinity / infinity = NaN.
+ final double b = Math.exp(LN2 * am1);
+ final double bm1 = Math.expm1(LN2 * am1);
+ bObm1 = b / bm1;
+ // Note:
+ // u in [0, 1]
+ // u == 0 : x == inf
+ // u == 1 : x == 1
+ // When a -> 1 then bias to infinity; otherwise bias to 1.
+ nextU = a <= THRESHOLD ?
+ // u in [0, 1)
+ InternalUtils::makeDouble :
+ // u in (0, 1]
+ InternalUtils::makeNonZeroDouble;
+ }
+
+ /**
+ * Copy constructor.
+ *
+ * @param rng Source of randomness.
+ * @param source Source to copy.
+ */
+ private RejectionZetaSampler(UniformRandomProvider rng,
RejectionZetaSampler source) {
+ this.rng = rng;
+ am1 = source.am1;
+ ram1 = source.ram1;
+ bObm1 = source.bObm1;
+ nextU = source.nextU;
+ }
+
+ @Override
+ public int sample() {
+ double u;
+ double v;
+ double x;
+ double t;
+ double tm1;
+ for (;;) {
+ // Generate iid uniform [0, 1] random variate U, V.
+ // Sampling only uses nextLong.
+ u = nextU.applyAsDouble(rng.nextLong());
+ v = InternalUtils.makeDouble(rng.nextLong());
+
+ // X = floor ( U^{-1/(a-1)} ) , X in [1, inf]
+ x = Math.floor(Math.pow(u, -ram1));
+
+ // T = (1 + 1/x)^(a-1)
+ // This can create T ~ 1 for large X so
+ // avoid precision loss using exp(log(1 + 1/x) * (a-1)).
+ t = Math.log1p(1 / x) * am1;
+ tm1 = Math.expm1(t);
+ t = Math.exp(t);
+
+ // Until:
+ // T-1 T
+ // VX --- <= -
+ // b-1 b
+
+ // Note: If X==infinity then T==1.
+ // Avoid infinity * (T - 1) == NaN by rearrangement:
+ // b T
+ // VX --- <= ---
+ // b-1 T-1
+
+ if (v * x * bObm1 <= t / tm1) {
+ // Truncates x >= 2^31 to integer max
+ return (int) x;
+ }
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return "Zeta deviate [" + rng.toString() + "]";
+ }
+
+ @Override
+ public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return new RejectionZetaSampler(rng, this);
+ }
+ }
+
+ /**
+ * Sample from the zeta distribution when the density is entirely
concentrated at x=1.
+ */
+ private static final class LargeExponentZetaSampler implements
SharedStateDiscreteSampler {
+ /** The single instance. */
+ static final LargeExponentZetaSampler INSTANCE = new
LargeExponentZetaSampler();
+
+ @Override
+ public int sample() {
+ return 1;
+ }
+
+ @Override
+ public String toString() {
+ return "Zeta(x=1) deviate";
+ }
+
+ @Override
+ public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ // No requirement for a new instance
+ return this;
+ }
+ }
+
+ /** Class contains only static methods. */
+ private ZetaSampler() {}
+
+ /**
+ * Creates a new zeta distribution sampler.
+ *
+ * <p>When {@code s > 54} the sample will always be 1. See the {@linkplain
ZetaSampler
+ * class-level} documentation for details.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param exponent Exponent.
+ * @return the sampler
+ * @throws IllegalArgumentException if {@code exponent <= 1} or is
non-finite.
+ */
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double exponent) {
+ InternalUtils.requireRangeClosed(MIN_EXPONENT,
Double.POSITIVE_INFINITY, exponent, "exponent");
+ return exponent <= REJECTION_EXPONENT ?
+ new RejectionZetaSampler(rng, exponent) :
+ LargeExponentZetaSampler.INSTANCE;
+ }
+}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerParametricTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerParametricTest.java
index a972464a..0cbc488e 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerParametricTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerParametricTest.java
@@ -28,6 +28,13 @@ import org.junit.jupiter.params.provider.MethodSource;
* Tests for random deviates generators.
*/
class DiscreteSamplerParametricTest {
+
+ /** Method to test two integer values. */
+ @FunctionalInterface
+ interface IntIntBiPredicate {
+ boolean test(int a, int b);
+ }
+
private static Iterable<DiscreteSamplerTestData> getSamplerTestData() {
return DiscreteSamplersList.list();
}
@@ -40,7 +47,8 @@ class DiscreteSamplerParametricTest {
check(sampleSize,
data.getSampler(),
data.getPoints(),
- data.getProbabilities());
+ data.getProbabilities(),
+ data.isRange());
}
/**
@@ -53,11 +61,13 @@ class DiscreteSamplerParametricTest {
* @param sampleSize Number of random values to generate.
* @param points Outcomes.
* @param expected Expected counts of the given outcomes.
+ * @param range True if the probabilities are for a range.
*/
private static void check(long sampleSize,
DiscreteSampler sampler,
int[] points,
- double[] expected) {
+ double[] expected,
+ boolean range) {
final ChiSquareTest chiSquareTest = new ChiSquareTest();
final int numTests = 50;
@@ -67,6 +77,13 @@ class DiscreteSamplerParametricTest {
final int numBins = points.length;
final long[] observed = new long[numBins];
+ // Support testing a probability mass function PMF(x)
+ // or a range p(x_i < X < x_i+1). If points are sorted
+ // we only require comparing using <= points[i].
+ final IntIntBiPredicate test = range ?
+ (x, y) -> x <= y :
+ (x, y) -> x == y;
+
// For storing chi2 larger than the critical value.
final List<Double> failedStat = new ArrayList<>();
try {
@@ -76,7 +93,7 @@ class DiscreteSamplerParametricTest {
final int value = sampler.sample();
for (int k = 0; k < numBins; k++) {
- if (value == points[k]) {
+ if (test.test(value, points[k])) {
++observed[k];
continue SAMPLE;
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerTestData.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerTestData.java
index 3a7ef85c..11bfbf54 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerTestData.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplerTestData.java
@@ -16,22 +16,42 @@
*/
package org.apache.commons.rng.sampling.distribution;
-import java.util.Arrays;
+import org.junit.jupiter.api.Assertions;
/**
* Data store for {@link DiscreteSamplerParametricTest}.
+ *
+ * <p>By default the probabilities are for {@code p_i(x_i)} where {@code x_i}
is the point.
+ * The data should contain probabilities for the entire range of {@code x}
expected to
+ * be produced by the sampler.
+ *
+ * <p>The data can be used to represent a cumulative probability distribution.
Points
+ * must be sorted and the probabilities correspond to:
+ * <pre>{@code
+ * P_0(X <= x_0)
+ * P_i+1(x_i < X <= x_i+1) , i in [0, n)
+ * }</pre>
*/
class DiscreteSamplerTestData {
private final DiscreteSampler sampler;
private final int[] points;
private final double[] probabilities;
+ private final boolean range;
DiscreteSamplerTestData(DiscreteSampler sampler,
int[] points,
- double[] probabilities) {
+ double[] probabilities,
+ boolean range) {
this.sampler = sampler;
this.points = points.clone();
this.probabilities = probabilities.clone();
+ this.range = range;
+ if (range) {
+ for (int k = 1; k < points.length; k++) {
+ Assertions.assertTrue(points[k - 1] < points[k],
+ "Points must be sorted for range probabilities");
+ }
+ }
}
public DiscreteSampler getSampler() {
@@ -46,13 +66,36 @@ class DiscreteSamplerTestData {
return probabilities.clone();
}
+ /**
+ * Check if the probabilities are for a range. Points must be sorted.
+ *
+ * @return true if the probabilities are for a range
+ * @since 1.8
+ */
+ public boolean isRange() {
+ return range;
+ }
+
@Override
public String toString() {
+ final StringBuilder sb = new StringBuilder(2048)
+ .append(sampler.toString()).append(':');
final int len = points.length;
- final String[] p = new String[len];
- for (int i = 0; i < len; i++) {
- p[i] = "p(" + points[i] + ")=" + probabilities[i];
+ if (range) {
+ sb.append(" p(<=").append(points[0]).append(")=")
+ .append(probabilities[0]);
+ for (int i = 1; i < len; i++) {
+ // Use a half-open interval, e.g. p((4,5])=
+ sb.append(" p((").append(points[i - 1]).append(',')
+ .append(points[i]).append("])=")
+ .append(probabilities[i]);
+ }
+ } else {
+ for (int i = 0; i < len; i++) {
+ sb.append(" p(").append(points[i]).append(")=")
+ .append(probabilities[i]);
+ }
}
- return sampler.toString() + ": " + Arrays.toString(p);
+ return sb.toString();
}
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplersList.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplersList.java
index 0db0e500..1128520d 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplersList.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DiscreteSamplersList.java
@@ -122,7 +122,7 @@ public final class DiscreteSamplersList {
MathArrays.sequence(5, 1, 1),
RejectionInversionZipfSampler.of(RandomAssert.createRNG(),
numElementsZipf, exponentCloseToOneZipf));
// Zipf (exponent = 0).
- add(LIST, MathArrays.sequence(5, 1, 1), new double[] {0.2, 0.2,
0.2, 0.2, 0.2},
+ add(LIST, MathArrays.sequence(5, 1, 1), new double[] {0.2, 0.2,
0.2, 0.2, 0.2}, false,
RejectionInversionZipfSampler.of(RandomAssert.createRNG(),
numElementsZipf, 0.0));
// Poisson ("inverse method").
@@ -180,16 +180,66 @@ public final class DiscreteSamplersList {
final int[] discretePoints = {0, 1, 2, 3, 4};
final double[] discreteProbabilities = {0.1, 0.2, 0.3, 0.4, 0.5};
final long[] discreteFrequencies = {1, 2, 3, 4, 5};
- add(LIST, discretePoints, discreteProbabilities,
+ add(LIST, discretePoints, discreteProbabilities, false,
MarsagliaTsangWangDiscreteSampler.Enumerated.of(RandomAssert.createRNG(),
discreteProbabilities));
- add(LIST, discretePoints, discreteProbabilities,
+ add(LIST, discretePoints, discreteProbabilities, false,
GuideTableDiscreteSampler.of(RandomAssert.createRNG(),
discreteProbabilities));
- add(LIST, discretePoints, discreteProbabilities,
+ add(LIST, discretePoints, discreteProbabilities, false,
AliasMethodDiscreteSampler.of(RandomAssert.createRNG(),
discreteProbabilities));
- add(LIST, discretePoints, discreteProbabilities,
+ add(LIST, discretePoints, discreteProbabilities, false,
FastLoadedDiceRollerDiscreteSampler.of(RandomAssert.createRNG(),
discreteFrequencies));
- add(LIST, discretePoints, discreteProbabilities,
+ add(LIST, discretePoints, discreteProbabilities, false,
FastLoadedDiceRollerDiscreteSampler.of(RandomAssert.createRNG(),
discreteProbabilities));
+
+ // Zeta distribution: x in [1, infinity].
+ // Range points generated with scipy.stats (1.18.1) zipf(s).
+ // from scipy.stats import zipf
+ // import numpy as np; np.set_printoptions(precision=17)
+ // x = np.unique(zipf(s).ppf([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
0.8, 0.9])).astype(int)
+ // p = np.append(zipf(s).cdf(x), [1])
+ // Convert cumulative to range probabilities:
+ // p[1:] -= p[:-1]
+
+ // s = 1.15
+ add(LIST,
+ new int[] {1, 3, 6, 17, 58, 256, 1742, 26005, 2641987,
Integer.MAX_VALUE},
+ new double[] {
+ 0.13784177793759045, 0.10108148525222727,
0.06720522534524914,
+ 0.09570277096807939, 0.09903685404346751,
0.09925485973116599,
+ 0.09988225352737812, 0.09999485138279085,
0.09999992334440055,
+ 0.09999999846765073},
+ true,
+ ZetaSampler.of(RandomAssert.createRNG(), 1.15));
+ // s = 1.35. Skewed to 1
+ add(LIST,
+ new int[] {1, 2, 4, 8, 18, 58, 417, Integer.MAX_VALUE},
+ new double[] {
+ 0.28908106624157415, 0.11340420378808813,
0.11008788688054949,
+ 0.09700240404611993, 0.09297268629333943,
0.09863885643999604,
+ 0.09887693680911303, 0.0999359595012198},
+ true,
+ ZetaSampler.of(RandomAssert.createRNG(), 1.35));
+ // s = 2.35. Very skewed to 1
+ add(LIST,
+ new int[] {1, 2, 3, Integer.MAX_VALUE},
+ new double[] {
+ 0.7107946035897981, 0.13941953571184507,
0.0537661789726942,
+ 0.09601968172566266},
+ true,
+ ZetaSampler.of(RandomAssert.createRNG(), 2.35));
+ // s = 1.02. Skewed to infinity. Tests sampler with truncation of
the distribution.
+ // scipy zipf is too slow to explore small s as it sums the pmf
for the cdf.
+ // Computed using CDF=(zeta(s, x+1) - zeta(s)) / zeta(s) using the
+ // the zeta function from mpmath (1.14.1) with x=[2**3, 2**7,
2**11, 2**15, 2**20, 2**25]
+ // and 50 digits of precision.
+ add(LIST,
+ new int[] {8, 128, 2048, 32768, 1048576, 33554432,
Integer.MAX_VALUE},
+ new double[] {0.05287102382096685, 0.0500627338177098,
0.04832778811226474,
+ 0.04577927420430394, 0.05377155271772868,
0.05017084772170899,
+ // Most of the density is above 2^31
+ 0.699016779605317},
+ true,
+ ZetaSampler.of(RandomAssert.createRNG(), 1.02));
} catch (Exception e) {
// CHECKSTYLE: stop Regexp
System.err.println("Unexpected exception while creating the list
of samplers: " + e);
@@ -228,7 +278,8 @@ public final class DiscreteSamplersList {
});
list.add(new DiscreteSamplerTestData(inverseMethodSampler,
points,
- getProbabilities(dist, points)));
+ getProbabilities(dist, points),
+ false));
}
/**
@@ -243,7 +294,8 @@ public final class DiscreteSamplersList {
final DiscreteSampler sampler) {
list.add(new DiscreteSamplerTestData(sampler,
points,
- getProbabilities(dist, points)));
+ getProbabilities(dist, points),
+ false));
}
/**
@@ -255,10 +307,12 @@ public final class DiscreteSamplersList {
private static void add(List<DiscreteSamplerTestData> list,
int[] points,
final double[] probabilities,
+ boolean range,
final DiscreteSampler sampler) {
list.add(new DiscreteSamplerTestData(sampler,
points,
- probabilities));
+ probabilities,
+ range));
}
/**
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZetaSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZetaSamplerTest.java
new file mode 100644
index 00000000..3c91a3a9
--- /dev/null
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ZetaSamplerTest.java
@@ -0,0 +1,180 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.rng.sampling.distribution;
+
+import java.time.Duration;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.RandomAssert;
+import
org.apache.commons.rng.sampling.distribution.ZetaSampler.RejectionZetaSampler;
+import org.apache.commons.rng.simple.RandomSource;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.ValueSource;
+
+/**
+ * Test for the {@link ZetaSampler}. The tests hit edge cases for the sampler.
+ */
+class ZetaSamplerTest {
+ /**
+ * Test the constructor with an {@code exponent <= 1}.
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {-1, 0, 1, Double.NaN})
+ void testConstructorThrowsWithBadExponent(double s) {
+ final UniformRandomProvider rng = RandomAssert.seededRNG();
+ Assertions.assertThrows(IllegalArgumentException.class,
+ () -> ZetaSampler.of(rng, s));
+ }
+
+ /**
+ * Test the SharedStateSampler implementation with exponents that are
above and below
+ * the threshold to switch the sampler bias from 1 to infinity.
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {1.01, 1.1, 1.3, 2, 100})
+ void testSharedStateSampler(double s) {
+ final UniformRandomProvider rng1 = RandomAssert.seededRNG();
+ final UniformRandomProvider rng2 = RandomAssert.seededRNG();
+ final SharedStateDiscreteSampler sampler1 = ZetaSampler.of(rng1, s);
+ final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
+ RandomAssert.assertProduceSameSequence(sampler1, sampler2);
+ }
+
+ /**
+ * Test a large exponent biases towards 1.
+ *
+ * <p>Note: This uses a package-private constructor for the rejection
sampler.
+ * The test validates that an exponent above 54 does not sample values
other
+ * than 1. If the exponent is too large (e.g. above 1025) then the sampler
will
+ * overflow to infinity when computing 2^(a-1). This will result in inf /
inf
+ * division producing a NaN and infinite recursion. This limits the
exponent used
+ * by this test.
+ *
+ * <p>In practice the factory method knows the rejection method will not
work
+ * and returns a sampler that always returns 1.
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {54.1, 100})
+ void testLargeExponent(double exponent) {
+ // x = floor ( U^{-1/(a-1)} )
+ // Bias towards 1 has U in (0, 1] to avoid U=0 where x=infinity
+ // Threshold at 1.9999... = pow(2^-53, -(1 / (a - 1))
+ // a = 1 - Math.log(0x1p-53) / Math.log(Math.nextDown(2.0)) =
54.00000000000001
+ // The sample should always be 1
+ final int expected = 1;
+ UniformRandomProvider rng;
+ DiscreteSampler s;
+ // u from first long : Can be [0, 1) or (0, 1]
+ // v from second long should control rejection
+ rng = createRNG(0, -1, 0, -1, 0, -1);
+ s = new RejectionZetaSampler(rng, exponent);
+ for (int i = 0; i < 3; i++) {
+ Assertions.assertEquals(expected, s.sample());
+ }
+ rng = createRNG(-1, -1, -1, -1, -1, -1);
+ s = new RejectionZetaSampler(rng, exponent);
+ for (int i = 0; i < 3; i++) {
+ Assertions.assertEquals(expected, s.sample());
+ }
+ // Any RNG
+ s = new RejectionZetaSampler(RandomAssert.createRNG(), exponent);
+ for (int i = 0; i < 100; i++) {
+ Assertions.assertEquals(expected, s.sample());
+ }
+ }
+
+ /**
+ * Test a tiny exponent cannot avoid small samples even though it should
+ * bias towards infinity.
+ */
+ @Test
+ void testTinyExponent() {
+ final double exponent = Math.nextUp(1.0);
+ // x = floor ( U^{-1/(a-1)} )
+ // Bias towards infinity has U in [0, 1) to avoid U=1
+ // But this cannot avoid small samples when U is close to 1.
+ // Math.pow(1.0, -(1/0x1p-52)) = 1
+ // Math.pow(Math.nextDown(1.0), -(1/0x1p-52)) = 1.6487212707001282
+ // Math.pow(Math.nextDown(Math.nextDown(1.0)), -(1/0x1p-52)) =
2.7182818284590455
+ // Math.pow(Math.nextDown(Math.nextDown(Math.nextDown(1.0))),
-(1/0x1p-52)) = 4.481689070338066
+ final UniformRandomProvider rng = createRNG(
+ 0, -1, // u=0.0; v~1.0
+ -1, -1, // u=1.0 - 2^-53; v~1.0
+
+ // Note: These lead to rejection in the current implementation
when v is high.
+ // So here we lower v to allow the rejection test to pass.
+ -2L << 11, -1L >>> 1, // u=1.0 - 2 * 2^-53; v~0.5
+ -3L << 11, -1L >>> 1 // u=1.0 - 3 * 2^-53; v~0.5
+ );
+ final DiscreteSampler s = ZetaSampler.of(rng, exponent);
+ Assertions.assertEquals(Integer.MAX_VALUE, s.sample());
+ // Expected samples from values above: 1.645; 2.718; 4.481.
+ Assertions.assertEquals(1, s.sample());
+ Assertions.assertEquals(2, s.sample());
+ Assertions.assertEquals(4, s.sample());
+
+ // Test the sampler can output samples (i.e. do not reject forever
with other u in [0, 1).
+ final int[][] sample = new int[1][];
+ UniformRandomProvider rng2 = RandomSource.L128_X256_MIX.create();
+ Assertions.assertTimeoutPreemptively(Duration.ofMillis(100), () ->
+ sample[0] = ZetaSampler.of(rng2, exponent).samples(1000).toArray()
+ );
+ // The survival function of the zeta distribution with tiny s:
+ // sf(x=2147483647; s=1 + 2^-52) = 0.9999999999999951
+ // We expect the sample to be at the upper bound.
+ // Use of a fresh RNG for repeat invocation if the test fails should
ensure
+ // this passes. Otherwise consider a fixed seed.
+ for (final int x : sample[0]) {
+ Assertions.assertEquals(Integer.MAX_VALUE, x, "Consider using a
fixed RNG seed");
+ }
+ }
+
+ /**
+ * Test the toString method for cases not hit in the rest of the test
suite.
+ * This test asserts the toString method always contains the string 'zeta'.
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {1.23, 100})
+ void testToString(double exponent) {
+ final UniformRandomProvider rng = RandomAssert.seededRNG();
+ final String s = ZetaSampler.of(rng,
exponent).toString().toLowerCase();
+ Assertions.assertTrue(s.contains("zeta"));
+ }
+
+ /**
+ * Creates the RNG to return the given values from the nextLong() method.
+ *
+ * @param values Long values
+ * @return the RNG
+ */
+ private static UniformRandomProvider createRNG(long... values) {
+ return new UniformRandomProvider() {
+ private int i;
+
+ @Override
+ public long nextLong() {
+ return values[i++];
+ }
+
+ @Override
+ public double nextDouble() {
+ throw new IllegalStateException("nextDouble cannot be trusted
to be in [0, 1) and should be ignored");
+ }
+ };
+ }
+}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 591e7e54..606b6f4b 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -56,6 +56,9 @@ If the output is not quite correct, check for invisible
trailing spaces!
<release version="1.8" date="TBD" description="
New features, updates and bug fixes (requires Java 8).
">
+ <action dev="aherbert" type="add" due-to="Alex Herbert" issue="RNG-203">
+ "ZetaSampler": Add a sampler for the zeta distribution.
+ </action>
<action dev="aherbert" type="update" due-to="Security scan, Alex
Herbert" issue="RNG-202">
"AliasMethodDiscreteSampler" and "GuideTableDiscreteSampler": Improve
documentation of the effects of the scaling factor on the internal
table
