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 dd8ca0e5 RNG-195: Extend support for alpha parameters in
DirichletSampler
dd8ca0e5 is described below
commit dd8ca0e56574111d4fc069b5a59a4607cd8acb44
Author: Alex Herbert <[email protected]>
AuthorDate: Sat Oct 10 19:04:12 2026 +0100
RNG-195: Extend support for alpha parameters in DirichletSampler
Changes have been made to output log samples from the gamma sampler. The
log samples can be scaled by a constant before normalising by their sum.
This supports Dirichlet samples where the gamma samples are all zero or
overflow due to tiny (<1e-19) or large (>1e33) concentration parameters.
---
.../AhrensDieterMarsagliaTsangGammaSampler.java | 105 ++++-
.../sampling/distribution/DirichletSampler.java | 250 ++++++++++--
...AhrensDieterMarsagliaTsangGammaSamplerTest.java | 81 +++-
.../distribution/DirichletSamplerTest.java | 433 ++++++++++++++++++++-
src/changes/changes.xml | 4 +
5 files changed, 811 insertions(+), 62 deletions(-)
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 42687f3b..e05c96cd 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
@@ -16,6 +16,8 @@
*/
package org.apache.commons.rng.sampling.distribution;
+import java.util.function.DoubleBinaryOperator;
+import java.util.function.DoubleUnaryOperator;
import org.apache.commons.rng.UniformRandomProvider;
/**
@@ -138,6 +140,10 @@ public class AhrensDieterMarsagliaTsangGammaSampler
private final double oneOverAlpha;
/** Optimization (see code). */
private final double bGSOptim;
+ /** Output mapping function for step 2. Accepts u and x. */
+ private final DoubleBinaryOperator out2;
+ /** Output mapping function for step 3. Accepts x. */
+ private final DoubleUnaryOperator out3;
/**
* Constructs a new instance.
@@ -145,15 +151,32 @@ public class AhrensDieterMarsagliaTsangGammaSampler
* @param rng Generator of uniformly distributed random numbers.
* @param alpha Alpha parameter of the distribution.
* @param theta Theta parameter of the distribution.
+ * @param logSamples Return the natural logarithm of the sample.
* @throws IllegalArgumentException if {@code alpha <= 0} or {@code
theta <= 0},
* or if {@code alpha} or {@code theta} are not finite
*/
AhrensDieterGammaSampler(UniformRandomProvider rng,
double alpha,
- double theta) {
+ double theta,
+ boolean logSamples) {
super(rng, alpha, theta);
oneOverAlpha = 1 / alpha;
bGSOptim = 1 + alpha / Math.E;
+ if (logSamples) {
+ final double logTheta = Math.log(theta);
+ // RNG-195: Support log gamma samples for the DirichletSampler
+ // x = ((1 + alpha / Math.E) * u)^1/alpha
+ // log(x) = log((1 + alpha / Math.E) * u) * 1/alpha
+ // = (log(1 + alpha / Math.E) + log(u)) * 1/alpha
+ // alpha < 2^-106 the log1p term (~2^-107) cannot be added to
log(u).
+ final double logBSGOptim = Math.log1p(alpha / Math.E);
+ out2 = (u, x) -> logTheta + (logBSGOptim + Math.log(u)) *
oneOverAlpha;
+ out3 = x -> logTheta + Math.log(x);
+ } else {
+ // RNG-195: Support for the existing gamma sampler behaviour
is unchanged
+ out2 = (u, x) -> theta * x;
+ out3 = x -> theta * x;
+ }
}
/**
@@ -167,6 +190,8 @@ public class AhrensDieterMarsagliaTsangGammaSampler
super(rng, source);
oneOverAlpha = source.oneOverAlpha;
bGSOptim = source.bGSOptim;
+ out2 = source.out2;
+ out3 = source.out3;
}
@Override
@@ -180,6 +205,11 @@ public class AhrensDieterMarsagliaTsangGammaSampler
if (p <= 1) {
// Step 2:
+ // Always taken when alpha < 2^-52 since 1 + alpha / E = 1.
+ // RNG-195: Support for the existing gamma sampler
behaviour is unchanged.
+ // This does not use exp(log(p) * 1/alpha) or log1p
variants to recycle
+ // the log term. Output can differ by a large relative
error.
+ // In the main use case for real samples using pow is
faster than exp+log.
final double x = Math.pow(p, oneOverAlpha);
final double u2 = rng.nextDouble();
@@ -187,15 +217,24 @@ public class AhrensDieterMarsagliaTsangGammaSampler
// Reject.
continue;
}
- return theta * x;
+ return out2.applyAsDouble(u, x);
}
// Step 3:
+ // -log(((1+alpha/E) - p) / alpha)
+ // -log((1+alpha/E) - p) + log(alpha)
+ // At limit alpha ~ 1
+ // max(1 + alpha / E) = 1.368, p > 1
+ // (1 + alpha / Math.E) - p < 0.368 and -> 0 as alpha reduces.
+ // No use for log1p.
+ // alpha ~ 1 sample range:
+ // -log(2^-52) = 36.04
+ // -log(1 / E) = 1
final double x = -Math.log((bGSOptim - p) * oneOverAlpha);
final double u2 = rng.nextDouble();
if (u2 <= Math.pow(x, alpha - 1)) {
- return theta * x;
+ return out3.applyAsDouble(x);
}
// Reject and continue.
}
@@ -228,6 +267,10 @@ public class AhrensDieterMarsagliaTsangGammaSampler
private final double cOptim;
/** Gaussian sampling. */
private final NormalizedGaussianSampler gaussian;
+ /** Output mapping function for step 1. Accepts v and Gaussian sample
x. */
+ private final DoubleBinaryOperator out1;
+ /** Output mapping function for step 2. Accepts v and log(v). */
+ private final DoubleBinaryOperator out2;
/**
* Constructs a new instance.
@@ -235,16 +278,32 @@ public class AhrensDieterMarsagliaTsangGammaSampler
* @param rng Generator of uniformly distributed random numbers.
* @param alpha Alpha parameter of the distribution.
* @param theta Theta parameter of the distribution.
+ * @param logSamples Return the natural logarithm of the sample.
* @throws IllegalArgumentException if {@code alpha <= 0} or {@code
theta <= 0},
* or if {@code alpha} or {@code theta} are not finite
*/
MarsagliaTsangGammaSampler(UniformRandomProvider rng,
double alpha,
- double theta) {
+ double theta,
+ boolean logSamples) {
super(rng, alpha, theta);
gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
+ // Note: As alpha -> large:
+ // dOptim -> alpha
+ // cOptim -> 0
+ // samples -> alpha
dOptim = alpha - ONE_THIRD;
cOptim = ONE_THIRD / Math.sqrt(dOptim);
+ if (logSamples) {
+ final double c = Math.log(theta) + Math.log(dOptim);
+ out1 = (v, x) -> c + 3 * Math.log1p(cOptim * x);
+ out2 = (v, lv) -> c + lv;
+ } else {
+ // RNG-195: Support for the existing gamma sampler behaviour
is unchanged
+ final double c = theta * dOptim;
+ out1 = (v, x) -> c * v;
+ out2 = (v, lv) -> c * v;
+ }
}
/**
@@ -259,6 +318,8 @@ public class AhrensDieterMarsagliaTsangGammaSampler
gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
dOptim = source.dOptim;
cOptim = source.cOptim;
+ out1 = source.out1;
+ out2 = source.out2;
}
@Override
@@ -277,11 +338,15 @@ public class AhrensDieterMarsagliaTsangGammaSampler
// Squeeze.
if (u < 1 - 0.0331 * x2 * x2) {
- return theta * dOptim * v;
+ return out1.applyAsDouble(v, x);
}
- if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + Math.log(v))) {
- return theta * dOptim * v;
+ // Support for large alpha where 1 + cOptim * x is 1.
+ // log(v) = log((1 + cOptim * x)^3)
+ // Ensures log output is consistent on both output pathways.
+ final double lv = 3 * Math.log1p(cOptim * x);
+ if (Math.log(u) < 0.5 * x2 + dOptim * (1 - v + lv)) {
+ return out2.applyAsDouble(v, lv);
}
}
}
@@ -346,9 +411,31 @@ public class AhrensDieterMarsagliaTsangGammaSampler
public static SharedStateContinuousSampler of(UniformRandomProvider rng,
double alpha,
double theta) {
+ return of(rng, alpha, theta, false);
+ }
+
+ /**
+ * Creates a new gamma distribution sampler.
+ *
+ * <p><Strong>This method is intentionally not part of the public
API</strong>.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param alpha Alpha parameter of the distribution (this is a shape
parameter).
+ * @param theta Theta parameter of the distribution (this is a scale
parameter).
+ * @param logSamples Return samples using the natural logarithm of the
sample.
+ * @return the sampler
+ * @throws IllegalArgumentException if {@code alpha <= 0} or {@code theta
<= 0},
+ * or if {@code alpha} or {@code theta} are not finite
+ * @since 1.8
+ */
+ // package-private for use in the DirichletSampler
+ static SharedStateContinuousSampler of(UniformRandomProvider rng,
+ double alpha,
+ double theta,
+ boolean logSamples) {
// Each sampler should check the input arguments.
return alpha < 1 ?
- new AhrensDieterGammaSampler(rng, alpha, theta) :
- new MarsagliaTsangGammaSampler(rng, alpha, theta);
+ new AhrensDieterGammaSampler(rng, alpha, theta, logSamples) :
+ new MarsagliaTsangGammaSampler(rng, alpha, theta, logSamples);
}
}
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DirichletSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DirichletSampler.java
index 0776ca93..a3781156 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DirichletSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/DirichletSampler.java
@@ -35,10 +35,143 @@ import
org.apache.commons.rng.sampling.SharedStateObjectSampler;
public abstract class DirichletSampler implements
SharedStateObjectSampler<double[]> {
/** The minimum number of categories. */
private static final int MIN_CATGEORIES = 2;
+ /** The minimum alpha before using log samples. See RNG-195. */
+ private static final double MIN_ALPHA = 0x1p-5;
+ /** The maximum alpha before using log samples. See RNG-195. */
+ private static final double MAX_ALPHA = 0x1p53;
+ /** The threshold for tiny alpha to use scaling.
+ * A tiny alpha is any value where log1p(alpha/E) + log(u) == log(u) for u
in (0, 1).
+ * The smallest magnitude for log(u) is 2^-53 so a tiny alpha is 2^-53 *
2^-53 = 2^-106.
+ * However this must be lowered by the amount of scaling
+ * applied to alpha. Set using a conservatively low threshold so that
+ * TINY_ALPHA * ALPHA_SCALE is below 2^-106. */
+ private static final double TINY_ALPHA = 0x1p-163;
+ /** Scaling to apply to tiny alpha where all log gamma samplers
+ * will return samples as log(u) / alpha with u in [0, 1).
+ * Any scaled alpha must be below 2^-106 after scaling so the log(u) /
alpha
+ * pathway is always used, and log1p(alpha/E) cannot be added to log(u).
+ * Set so that log(2^-53) / (scale * MIN_VALUE) is finite allowing support
+ * for alpha down to the lowest non-zero double value when the uniform
deviate
+ * is the smallest non-zero value. */
+ private static final double ALPHA_SCALE = 0x1p56;
+ /** Normalizer for the gamma samples. */
+ private final Normalizer normalizer;
/** RNG (used for the toString() method). */
private final UniformRandomProvider rng;
+ /**
+ * Normalize an array of samples to unit length.
+ */
+ // package-private for testing
+ enum Normalizer {
+ /** Normalise real-valued samples.
+ * The sum of the samples is expected to be a finite normal double.
+ * A sub-normal sum cannot be normalised as this uses multiplication
by the
+ * reciprocal of the sum (which overflows when sum is tiny). */
+ REAL {
+ @Override
+ boolean normalize(double[] x) {
+ // Samples are never length zero
+ double norm = x[0];
+ for (int i = 1; i < x.length; i++) {
+ norm += x[i];
+ }
+ // Normalize by dividing by the sum of the samples
+ norm = 1.0 / norm;
+ // Detect an invalid normalization,
+ // e.g. cases of all zero samples, or one or more infinite
samples.
+ if (isNonZeroPositiveFinite(norm)) {
+ // Normalise
+ for (int i = 0; i < x.length; i++) {
+ x[i] *= norm;
+ }
+ return true;
+ }
+ return false;
+ }
+ },
+ /** Normalise log-valued samples.
+ * This will normalize any sample even including infinity (or NaN). */
+ LOG {
+ @Override
+ boolean normalize(double[] x) {
+ return normalize(x, 1.0);
+ }
+ },
+ /** Normalise scaled log-valued samples.
+ * This will normalize any sample even including infinity (or NaN). */
+ SCALED_LOG {
+ @Override
+ boolean normalize(double[] x) {
+ return normalize(x, ALPHA_SCALE);
+ }
+ };
+
+ /**
+ * Normalise the sample to unit length using the sum.
+ *
+ * <p>Returns false if the sample could not be normalized. This occurs
when
+ * all terms underflow and the sum is zero, or one or more terms
+ * overflow and the sum is infinity.
+ *
+ * @param x Sample.
+ * @return true if the sample was normalized
+ */
+ abstract boolean normalize(double[] x);
+
+ /**
+ * Normalise the log sample to unit length using the sum.
+ * Log samples are multiplied by the provided {@code scale}. Used
+ * when log samples originate from a scaled alpha parameter passed
+ * to the gamma sampler.
+ *
+ * @param x Sample.
+ * @param scale Scale.
+ * @return true
+ */
+ static boolean normalize(double[] x, double scale) {
+ // Scale log numbers so we can sum them.
+ // Subtracting the maximum log value stabilizes the sum.
+ // The maximum exp term will be 1, all other terms <=1 and the sum
is finite.
+ double max = x[0];
+ for (int i = 1; i < x.length; i++) {
+ max = Math.max(max, x[i]);
+ }
+ double norm = 0;
+ if (Double.isFinite(max)) {
+ // A finite max will normalize any sample including -infinity.
+ for (int i = 0; i < x.length; i++) {
+ // If scale > 1 then log samples are moved further apart.
+ // It is likely only the max will contribute to the sum
+ // and other sample values are 0.
+ x[i] = Math.exp((x[i] - max) * scale);
+ norm += x[i];
+ }
+ } else {
+ // Any value matching the max is normalised to length 1.
+ // A side effect is that any NaNs will normalize to 1
+ // as Math.max puts NaN above all other numbers.
+ // NaNs or positive infinity are not expected from the log
gamma samples.
+ final long b = Double.doubleToLongBits(max);
+ for (int i = 0; i < x.length; i++) {
+ x[i] = Double.doubleToLongBits(x[i]) == b ? 1 : 0;
+ norm += x[i];
+ }
+ }
+ // Note: In contrast to the real case where samples
+ // are expected across all categories, log samples
+ // for tiny alpha are typically biased to 1 category.
+ // Using divide ensures extreme samples for tiny
+ // alpha have 1.0 in the chosen category and not 1-2^53
+ // which can happen when multiplying by the reciprocal.
+ for (int i = 0; i < x.length; i++) {
+ x[i] /= norm;
+ }
+ return true;
+ }
+ }
+
/**
* Sample from a Dirichlet distribution with different concentration
parameters
* for each category.
@@ -52,10 +185,12 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
*
* @param rng Generator of uniformly distributed random numbers.
* @param samplers Samplers for each category.
+ * @param normalizer Normalizer for the samples.
*/
GeneralDirichletSampler(UniformRandomProvider rng,
- SharedStateContinuousSampler[] samplers) {
- super(rng);
+ SharedStateContinuousSampler[] samplers,
+ Normalizer normalizer) {
+ super(rng, normalizer);
// Array is stored directly as it is generated within the
DirichletSampler class
this.samplers = samplers;
}
@@ -76,7 +211,7 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
for (int i = 0; i < newSamplers.length; i++) {
newSamplers[i] = samplers[i].withUniformRandomProvider(rng);
}
- return new GeneralDirichletSampler(rng, newSamplers);
+ return new GeneralDirichletSampler(rng, newSamplers,
getNormalizer());
}
}
@@ -96,11 +231,13 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
* @param rng Generator of uniformly distributed random numbers.
* @param k Number of categories.
* @param sampler Sampler for the categories.
+ * @param normalizer Normalizer for the samples.
*/
SymmetricDirichletSampler(UniformRandomProvider rng,
int k,
- SharedStateContinuousSampler sampler) {
- super(rng);
+ SharedStateContinuousSampler sampler,
+ Normalizer normalizer) {
+ super(rng, normalizer);
this.k = k;
this.sampler = sampler;
}
@@ -117,7 +254,7 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
@Override
public SymmetricDirichletSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- return new SymmetricDirichletSampler(rng, k,
sampler.withUniformRandomProvider(rng));
+ return new SymmetricDirichletSampler(rng, k,
sampler.withUniformRandomProvider(rng), getNormalizer());
}
}
@@ -125,9 +262,20 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
* Constructs a new instance.
*
* @param rng Generator of uniformly distributed random numbers.
+ * @param normalizer Normalizer for the samples.
*/
- DirichletSampler(UniformRandomProvider rng) {
+ DirichletSampler(UniformRandomProvider rng, Normalizer normalizer) {
this.rng = rng;
+ this.normalizer = normalizer;
+ }
+
+ /**
+ * Gets the normalizer.
+ *
+ * @return the normalizer
+ */
+ Normalizer getNormalizer() {
+ return normalizer;
}
/** {@inheritDoc} */
@@ -141,25 +289,15 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
public double[] sample() {
// Create Gamma(alpha_i, 1) deviates for all alpha
final double[] y = new double[getK()];
- double norm = 0;
for (int i = 0; i < y.length; i++) {
- final double yi = nextGamma(i);
- norm += yi;
- y[i] = yi;
+ y[i] = nextGamma(i);
}
- // Normalize by dividing by the sum of the samples
- norm = 1.0 / norm;
- // Detect an invalid normalization, e.g. cases of all zero samples
- if (!isNonZeroPositiveFinite(norm)) {
+ if (!normalizer.normalize(y)) {
// Sample again using recursion.
// A stack overflow due to a broken RNG will eventually occur
// rather than the alternative which is an infinite loop.
return sample();
}
- // Normalise
- for (int i = 0; i < y.length; i++) {
- y[i] *= norm;
- }
return y;
}
@@ -195,11 +333,22 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
public static DirichletSampler of(UniformRandomProvider rng,
double... alpha) {
validateNumberOfCategories(alpha.length);
+ // Check for the maximum alpha.
+ // The mean of each gamma sample is the concentration parameter alpha.
+ // All samples must be normalised to unit length using the sum.
+ // Use log samples when samples will likely over/underflow on the
largest
+ // sample (create a zero or infinite sample). These can be scaled by
the
+ // largest sample before summing.
+ double maxAlpha = alpha[0];
+ for (int i = 1; i < alpha.length; i++) {
+ maxAlpha = Math.max(maxAlpha, alpha[i]);
+ }
+ final boolean logSamples = requireLogSamples(maxAlpha);
final SharedStateContinuousSampler[] samplers = new
SharedStateContinuousSampler[alpha.length];
for (int i = 0; i < samplers.length; i++) {
- samplers[i] = createSampler(rng, alpha[i]);
+ samplers[i] = createSampler(rng, alpha[i], logSamples, maxAlpha);
}
- return new GeneralDirichletSampler(rng, samplers);
+ return new GeneralDirichletSampler(rng, samplers,
createNormalizer(logSamples, maxAlpha));
}
/**
@@ -217,8 +366,9 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
int k,
double alpha) {
validateNumberOfCategories(k);
- final SharedStateContinuousSampler sampler = createSampler(rng, alpha);
- return new SymmetricDirichletSampler(rng, k, sampler);
+ final boolean logSamples = requireLogSamples(alpha);
+ final SharedStateContinuousSampler sampler = createSampler(rng, alpha,
logSamples, alpha);
+ return new SymmetricDirichletSampler(rng, k, sampler,
createNormalizer(logSamples, alpha));
}
/**
@@ -239,19 +389,41 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
*
* @param rng Generator of uniformly distributed random numbers.
* @param alpha Concentration parameter.
+ * @param logSamples Return samples using the natural logarithm of the
sample.
+ * @param maxAlpha Maximum alpha.
* @return the sampler
* @throws IllegalArgumentException if the concentration parameter is not
strictly positive.
*/
private static SharedStateContinuousSampler
createSampler(UniformRandomProvider rng,
- double alpha) {
+ double alpha,
+ boolean
logSamples,
+ double maxAlpha)
{
InternalUtils.requireStrictlyPositiveFinite(alpha, "alpha
concentration");
// Create a Gamma(shape=alpha, scale=1) sampler.
- if (alpha == 1) {
+ if (alpha == 1 && !logSamples) {
// Special case
// Gamma(shape=1, scale=1) == Exponential(mean=1)
return ZigguratSampler.Exponential.of(rng);
}
- return AhrensDieterMarsagliaTsangGammaSampler.of(rng, alpha, 1);
+
+ // -----
+ // Alpha scaling
+ //
+ // When alpha is very small 1/alpha will overflow.
+ // This invalidates the gamma sampler which creates samples using:
+ // u in [0, 1) : p = u * (1 + alpha / E) ~ u
+ // x = p^(1/alpha)
+ // For any p < 1 the result is zero. Log samples use:
+ // log(p) / alpha
+ // This is stable until alpha approaches MIN_NORMAL where it can
overflow.
+ // Since this sampler normalizes log samples using scaling so they can
be
+ // summed the solution to tiny alpha is to perform some scaling on
alpha
+ // when all samplers use the same gamma sampler method. Log samples
+ // must be rescaled before normalisation:
+ // log(p) / alpha = scale * log(p) / (scale * alpha)
+ // -----
+ final double a = maxAlpha < TINY_ALPHA ? alpha * ALPHA_SCALE : alpha;
+ return AhrensDieterMarsagliaTsangGammaSampler.of(rng, a, 1,
logSamples);
}
/**
@@ -263,4 +435,30 @@ public abstract class DirichletSampler implements
SharedStateObjectSampler<doubl
private static boolean isNonZeroPositiveFinite(double x) {
return x > 0 && x < Double.POSITIVE_INFINITY;
}
+
+ /**
+ * Check if the gamma sampler should use log samples. This extends support
for
+ * small and large concentration parameters where the samples cannot be
normalised
+ * to unit length without underflow or overflow.
+ *
+ * @param alpha Maximum alpha parameter.
+ * @return true if log samples are required
+ */
+ private static boolean requireLogSamples(double alpha) {
+ return alpha < MIN_ALPHA || alpha > MAX_ALPHA;
+ }
+
+ /**
+ * Creates the normalizer for the samples.
+ *
+ * @param logSamples Return samples using the natural logarithm of the
sample.
+ * @param alpha Maximum alpha parameter.
+ * @return the normalizer
+ */
+ private static Normalizer createNormalizer(boolean logSamples, double
alpha) {
+ if (logSamples) {
+ return alpha < TINY_ALPHA ? Normalizer.SCALED_LOG : Normalizer.LOG;
+ }
+ return Normalizer.REAL;
+ }
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSamplerTest.java
index 42e57a28..ed872e3a 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSamplerTest.java
@@ -50,19 +50,21 @@ class AhrensDieterMarsagliaTsangGammaSamplerTest {
}
/**
- * Test the SharedStateSampler implementation.
+ * Test the SharedStateSampler implementation with {@code alpha < 1}.
*/
@Test
void testSharedStateSamplerWithAlphaBelowOne() {
- testSharedStateSampler(0.5, 3.456);
+ testSharedStateSampler(0.5, 3.456, false);
+ testSharedStateSampler(0.5, 3.456, true);
}
/**
- * Test the SharedStateSampler implementation.
+ * Test the SharedStateSampler implementation with {@code alpha > 1}.
*/
@Test
void testSharedStateSamplerWithAlphaAboveOne() {
- testSharedStateSampler(3.5, 3.456);
+ testSharedStateSampler(3.5, 3.456, false);
+ testSharedStateSampler(3.5, 3.456, true);
}
/**
@@ -70,17 +72,82 @@ class AhrensDieterMarsagliaTsangGammaSamplerTest {
*
* @param alpha Alpha.
* @param theta Theta.
+ * @param logSamples Use log samples.
*/
- private static void testSharedStateSampler(double alpha, double theta) {
+ private static void testSharedStateSampler(double alpha, double theta,
boolean logSamples) {
final UniformRandomProvider rng1 = RandomAssert.seededRNG();
final UniformRandomProvider rng2 = RandomAssert.seededRNG();
- // Use instance constructor not factory constructor to exercise 1.X
public API
- final AhrensDieterMarsagliaTsangGammaSampler sampler1 =
+ final SharedStateContinuousSampler sampler1 = logSamples ?
+ AhrensDieterMarsagliaTsangGammaSampler.of(rng1, alpha, theta,
true) :
+ // Use instance constructor not factory constructor to exercise
1.X public API
new AhrensDieterMarsagliaTsangGammaSampler(rng1, alpha, theta);
final SharedStateContinuousSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
+ /**
+ * Test the log samples implementation with {@code alpha < 1}.
+ */
+ @Test
+ void testLogSamplesWithAlphaBelowOne() {
+ // Note: ulp passes with N=1,000,000. N has been lowered for time
savings.
+ final int n = 1000;
+ testLogSamples(0.5, 3.456, n, 25);
+ // alpha = 2^-52 is where (1 + alpha/E) == 1 (See RNG-195)
+ testLogSamples(0x1p-50, 1.0, n, 50);
+ testLogSamples(0x1p-52, 1.0, n, 50);
+ testLogSamples(0x1p-54, 1.0, n, 50);
+ }
+
+ /**
+ * Test the log samples implementation with {@code alpha > 1}.
+ */
+ @Test
+ void testLogSamplesWithAlphaAboveOne() {
+ // Note: ulp passes with N=1,000,000. N has been lowered for time
savings.
+ final int n = 1000;
+ testLogSamples(3.5, 3.456, n, 25);
+ // As alpha -> large then samples -> alpha
+ // Occurs consistently > 3e32 (See RNG-195)
+ testLogSamples(1e30, 1.0, n, 50);
+ testLogSamples(1e32, 1.0, n, 50);
+ testLogSamples(1e34, 1.0, n, 50);
+ }
+
+ /**
+ * Test the log samples implementation matches the real samples output.
+ *
+ * @param alpha Alpha.
+ * @param theta Theta.
+ * @param ulps Units-in-the-last place for matching output.
+ */
+ private static void testLogSamples(double alpha, double theta, int n, int
ulps) {
+ // May spuriously fail. Allow reruns to use a different RNG.
+ final UniformRandomProvider[] rngs = RandomAssert.createRNG(2);
+ final UniformRandomProvider rng1 = rngs[0];
+ final UniformRandomProvider rng2 = rngs[1];
+ final SharedStateContinuousSampler sampler1 =
+ AhrensDieterMarsagliaTsangGammaSampler.of(rng1, alpha, theta,
false);
+ final SharedStateContinuousSampler sampler2 =
+ AhrensDieterMarsagliaTsangGammaSampler.of(rng2, alpha, theta,
true);
+ for (int i = 0; i < n; i++) {
+ // The sequence may not be exactly the same.
+ final double x1 = sampler1.sample();
+ final double x2 = sampler2.sample();
+ // Use exp on the log sample.
+ final double s2 = Math.exp(x2);
+ Assertions.assertEquals(x1, s2, ulps * Math.ulp(x1),
+ () -> String.format("exp %s %s ulp=%d", alpha, theta, (long)
((s2 - x1) / Math.ulp(x1))));
+ // Using log on the real sample has issues mapping
+ // the domains and relative error can be large when x1 -> 1
(log(x1) -> 0)
+ if (x1 < 0.95 && x1 > 1.05) {
+ final double s1 = Math.log(x1);
+ Assertions.assertEquals(x2, s1, ulps * Math.ulp(x2),
+ () -> String.format("log %s %s ulp=%d", alpha, theta,
(long) ((s1 - x2) / Math.ulp(x2))));
+ }
+ }
+ }
+
/**
* Test the toString method. This is added to ensure coverage as the
factory constructor
* used in other tests does not create an instance of the wrapper class.
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DirichletSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DirichletSamplerTest.java
index 1dc6c6c0..222390ea 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DirichletSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/DirichletSamplerTest.java
@@ -16,13 +16,19 @@
*/
package org.apache.commons.rng.sampling.distribution;
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.time.Duration;
import java.util.Arrays;
+import java.util.function.Supplier;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.stat.correlation.Covariance;
import org.apache.commons.math3.stat.descriptive.moment.Mean;
+import org.apache.commons.math3.stat.inference.ChiSquareTest;
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.sampling.distribution.DirichletSampler.Normalizer;
import org.apache.commons.rng.simple.RandomSource;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
@@ -31,6 +37,126 @@ import org.junit.jupiter.api.Test;
* Test for {@link DirichletSampler}.
*/
class DirichletSamplerTest {
+ /** Threshold for tiny alpha.
+ * Set where log1p(alpha/Math.E) + log(u) == log(u) for u in [0, 1). */
+ private static final double TINY_ALPHA = 0x1p-106;
+
+ @Test
+ void testRealNormalizer() {
+ final Normalizer normalizer = Normalizer.REAL;
+ // Fail on all zeros
+ assertNormalizer(normalizer, new double[] {0, 0, 0}, null);
+ // Fail when sum is infinite
+ assertNormalizer(normalizer, new double[] {Double.POSITIVE_INFINITY,
1, 1}, null);
+ assertNormalizer(normalizer, new double[] {Double.MAX_VALUE,
Double.MAX_VALUE}, null);
+ // A sub-normal sum cannot divide using multiplication of the
reciprocal
+ assertNormalizer(normalizer, new double[] {Double.MIN_VALUE,
Double.MIN_VALUE}, null);
+ assertNormalizer(normalizer, new double[] {0, Double.MIN_VALUE}, null);
+ assertNormalizer(normalizer, new double[] {0x1p-1025, 0x1p-1025},
null);
+ // Normalize to unit length
+ assertNormalizer(normalizer, new double[] {1, 1}, new double[] {0.5,
0.5});
+ assertNormalizer(normalizer, new double[] {1, 2, 3, 4}, new double[]
{1.0 / 10, 2.0 / 10, 3.0 / 10, 4.0 / 10});
+ assertNormalizer(normalizer, new double[] {0, Double.MAX_VALUE}, new
double[] {0, 1});
+ assertNormalizer(normalizer, new double[] {0, Double.MIN_NORMAL}, new
double[] {0, 1});
+ }
+
+ @Test
+ void testLogNormalizer() {
+ final Normalizer normalizer = Normalizer.LOG;
+ final double inf = Double.POSITIVE_INFINITY;
+ // Infinite sample
+ assertNormalizer(normalizer, new double[] {inf, 1, 1}, new double[]
{1, 0, 0});
+ assertNormalizer(normalizer, new double[] {inf, 1, inf}, new double[]
{0.5, 0, 0.5});
+ // Finite max normalises even with -infinity
+ assertNormalizer(normalizer, new double[] {-inf,
Math.log(Double.MIN_VALUE)}, new double[] {0, 1});
+ // Sum of exp(x_i) overflows
+ assertNormalizer(normalizer, new double[] {Double.MAX_VALUE,
Double.MAX_VALUE}, new double[] {0.5, 0.5});
+ // Sum of exp(x_i) underflows
+ assertNormalizer(normalizer, new double[] {-Double.MAX_VALUE,
-Double.MAX_VALUE}, new double[] {0.5, 0.5});
+ // Normalize to unit length
+ // Equal log values
+ final double[] e1 = {1.0 / 3, 1.0 / 3, 1.0 / 3};
+ assertNormalizer(normalizer, new double[] {0, 0, 0}, e1);
+ assertNormalizer(normalizer, new double[] {42, 42, 42}, e1);
+ assertNormalizer(normalizer, new double[] {-13, -13, -13}, e1);
+ // Scale independent normalisation of finite real values
+ final double[] e2 = {Math.exp(-1), Math.exp(-2), Math.exp(-3)};
+ Normalizer.REAL.normalize(e2);
+ assertNormalizer(normalizer, new double[] {-1, -2, -3}, e2);
+ assertNormalizer(normalizer, new double[] {-11, -12, -13}, e2);
+ assertNormalizer(normalizer, new double[] {-10001, -10002, -10003},
e2);
+ assertNormalizer(normalizer, new double[] {10003, 10002, 10001}, e2);
+ // Finite max normalises even with -infinity
+ final double[] e3 = {Math.exp(-1), Math.exp(-2), Math.exp(-inf)};
+ Normalizer.REAL.normalize(e3);
+ assertNormalizer(normalizer, new double[] {-1, -2, -inf}, e3);
+
+ // Side effect: NaN is considered the maximum sample.
+ final double nan = Double.NaN;
+ assertNormalizer(normalizer, new double[] {nan, 1, 1}, new double[]
{1, 0, 0});
+ assertNormalizer(normalizer, new double[] {nan, 1, nan}, new double[]
{0.5, 0, 0.5});
+ assertNormalizer(normalizer, new double[] {inf, nan}, new double[] {0,
1});
+ }
+
+ @Test
+ void testScaledLogNormalizer() {
+ final Normalizer normalizer = Normalizer.SCALED_LOG;
+ final double inf = Double.POSITIVE_INFINITY;
+ // Infinite sample
+ assertNormalizer(normalizer, new double[] {inf, 1, 1}, new double[]
{1, 0, 0});
+ assertNormalizer(normalizer, new double[] {inf, 1, inf}, new double[]
{0.5, 0, 0.5});
+ // Finite max normalises even with -infinity
+ assertNormalizer(normalizer, new double[] {-inf,
Math.log(Double.MIN_VALUE)}, new double[] {0, 1});
+ // Sum of exp(x_i) overflows
+ assertNormalizer(normalizer, new double[] {Double.MAX_VALUE,
Double.MAX_VALUE}, new double[] {0.5, 0.5});
+ // Sum of exp(x_i) underflows
+ assertNormalizer(normalizer, new double[] {-Double.MAX_VALUE,
-Double.MAX_VALUE}, new double[] {0.5, 0.5});
+ // Normalize to unit length
+ // Equal log values
+ final double[] e1 = {1.0 / 3, 1.0 / 3, 1.0 / 3};
+ assertNormalizer(normalizer, new double[] {0, 0, 0}, e1);
+ assertNormalizer(normalizer, new double[] {42, 42, 42}, e1);
+ assertNormalizer(normalizer, new double[] {-13, -13, -13}, e1);
+ // Scale independent normalisation of finite real values.
+ // ---
+ // Log values are scaled by 2^56 so we have to downscale the provided
sample.
+ final double downscale = 0x1p-56;
+ // ---
+ final double[] e2 = {Math.exp(-1), Math.exp(-2), Math.exp(-3)};
+ Normalizer.REAL.normalize(e2);
+ assertNormalizer(normalizer, new double[] {-1 * downscale, -2 *
downscale, -3 * downscale}, e2);
+ assertNormalizer(normalizer, new double[] {-4 * downscale, -5 *
downscale, -6 * downscale}, e2);
+ assertNormalizer(normalizer, new double[] {-400 * downscale, -401 *
downscale, -402 * downscale}, e2);
+ // Finite max normalises even with -infinity
+ final double[] e3 = {Math.exp(-1), Math.exp(-2), Math.exp(-inf)};
+ Normalizer.REAL.normalize(e3);
+ assertNormalizer(normalizer, new double[] {-1 * downscale, -2 *
downscale, -inf}, e3);
+
+ // Side effect: NaN is considered the maximum sample.
+ final double nan = Double.NaN;
+ assertNormalizer(normalizer, new double[] {nan, 1, 1}, new double[]
{1, 0, 0});
+ assertNormalizer(normalizer, new double[] {nan, 1, nan}, new double[]
{0.5, 0, 0.5});
+ assertNormalizer(normalizer, new double[] {inf, nan}, new double[] {0,
1});
+ }
+
+ private static void assertNormalizer(Normalizer normalizer, double[] x,
double[] expected) {
+ final double[] y = x.clone();
+ final boolean result = normalizer.normalize(x);
+ final Supplier<String> msg = () -> Arrays.toString(y) + " -> " +
Arrays.toString(x);
+ if (expected == null) {
+ Assertions.assertFalse(result, msg);
+ Assertions.assertArrayEquals(y, x, msg);
+ } else {
+ Assertions.assertTrue(result, msg);
+ final double sum = Arrays.stream(expected).sum();
+ Assertions.assertEquals(1, sum, 0x1p-52, "Expected sum should be
1");
+ // Normalization may not be exact
+ for (int i = 0; i < x.length; i++) {
+ Assertions.assertEquals(expected[i], x[i],
Math.ulp(expected[i]), msg);
+ }
+ }
+ }
+
@Test
void testDistributionThrowsWithInvalidNumberOfCategories() {
final UniformRandomProvider rng = RandomAssert.seededRNG();
@@ -110,42 +236,267 @@ class DirichletSamplerTest {
assertSample(2, sampler.sample());
}
+ /**
+ * Test the assumption that a scaled alpha parameter to a log gamma
sampler will
+ * return the scaled samples. This works when alpha is very small as the
sampler
+ * creates samples using log(p) / alpha == scale * log(p) / (scale *
alpha).
+ */
+ @Test
+ void testScalingAssumptions() {
+ // The threshold where log1p(alpha/E) + log(u) == log(u)
+ // The smallest magnitude for log(u) in [0, 1) is log(1-2^-53) =
-1.1102230246251565E-16.
+ final double l = Math.log(Math.nextDown(1.0));
+ Assertions.assertEquals(l, Math.log1p(TINY_ALPHA / Math.E) + l);
+
+ final RandomSource randomSource = RandomSource.XO_RO_SHI_RO_128_PP;
+ final byte[] seed = randomSource.createSeed();
+ final UniformRandomProvider rng = randomSource.create(seed);
+ final double scale = 0x1p53;
+ final double tiny = TINY_ALPHA / scale;
+ for (int i = 0; i < 10; i++) {
+ final double alpha1 = rng.nextDouble() * tiny;
+ final double alpha2 = alpha1 * scale;
+ final UniformRandomProvider rng1 = randomSource.create(seed);
+ final UniformRandomProvider rng2 = randomSource.create(seed);
+ // log samples
+ final SharedStateContinuousSampler sampler1 =
AhrensDieterMarsagliaTsangGammaSampler.of(rng1, alpha1, 1, true);
+ final SharedStateContinuousSampler sampler2 =
AhrensDieterMarsagliaTsangGammaSampler.of(rng2, alpha2, 1, true);
+ RandomAssert.assertProduceSameSequence(sampler1, () ->
sampler2.sample() * scale);
+ }
+ }
+
+ /**
+ * Test extreme concentration parameters.
+ * This tests concentration parameters that failed in a previous version
+ * due to a sum of gamma samples that was either zero or infinity. This
could
+ * not be used to normalise the sum to 1.
+ * Critical thresholds were 1e-19 and 2^1023 for unbounded recursion.
+ * See RNG-195.
+ */
+ @Test
+ void testExtremeConcentration() {
+ // Thresholds to use log gamma samples are 2^-5 and 2^53.
+ // The assertion checks the sum of the sample is 1.
+ assertExtremeConcentration(false, 0, 1e-5, 1e-6, 1e-7);
+ assertExtremeConcentration(false, 0, 1e18, 1e19, 1e20);
+ assertExtremeConcentration(true, 3, 1e50);
+ assertExtremeConcentration(true, 3, 1e-8);
+
+ // Tiny alpha ~ 2^-106
+ assertExtremeConcentration(false, 0, TINY_ALPHA / 1.5, TINY_ALPHA / 2,
TINY_ALPHA, 2.5);
+ assertExtremeConcentration(true, 2, Double.MIN_VALUE);
+ // One dominant concentration
+ assertExtremeConcentration(false, 0, TINY_ALPHA / 1e200, TINY_ALPHA /
1.5, TINY_ALPHA / 1e100);
+ // Two dominant concentrations
+ assertExtremeConcentration(false, 0, TINY_ALPHA / 1.5, TINY_ALPHA /
1.5, TINY_ALPHA / 1e100);
+
+ // Large alpha
+ assertExtremeConcentration(false, 0, 0x1p1022, 0x1p1023, 0x1.5p1023);
+ assertExtremeConcentration(true, 2, Double.MAX_VALUE);
+
+ // A large alpha and tiny alpha will create some log samples as
-infinity.
+ // Normalisation should handle this and this category will never be
observed.
+ assertExtremeConcentration(false, 0, 0x1p1023, TINY_ALPHA / 2);
+ assertExtremeConcentration(false, 0, Double.MIN_VALUE,
Double.MAX_VALUE);
+ }
+
+ private static void assertExtremeConcentration(boolean symmetric, int k,
double... alpha) {
+ final RandomSource randomSource = RandomSource.XO_RO_SHI_RO_128_PP;
+ final byte[] seed = randomSource.createSeed();
+ final UniformRandomProvider rng1 = randomSource.create(seed);
+ final DirichletSampler sampler1;
+ final double[] alphas;
+ if (symmetric) {
+ sampler1 = DirichletSampler.symmetric(rng1, k, alpha[0]);
+ alphas = new double[k];
+ Arrays.fill(alphas, alpha[0]);
+ } else {
+ sampler1 = DirichletSampler.of(rng1, alpha);
+ alphas = alpha.clone();
+ }
+
+ // Test sample generation does not infinitely recurse.
+ // 10 milliseconds is fine on a fast machine; add some headroom
+ final int n = 100;
+ final Duration timeout = Duration.ofMillis(100);
+ final double[][][] samples = new double[1][][];
+
+ Assertions.assertTimeoutPreemptively(timeout,
+ () -> samples[0] = sampler1.samples(n).toArray(double[][]::new));
+ final double[][] x = samples[0];
+ Assertions.assertNotNull(x);
+
+ if (symmetric) {
+ // Check the general case matches
+ final UniformRandomProvider rng2 = randomSource.create(seed);
+ final DirichletSampler sampler2 = DirichletSampler.of(rng2,
alphas);
+
+ Assertions.assertTimeoutPreemptively(timeout,
+ () -> samples[0] =
sampler2.samples(n).toArray(double[][]::new));
+ final double[][] y = samples[0];
+ // Check nested arrays
+ Assertions.assertArrayEquals(x, y);
+ }
+
+ // Checks on the sample. Correct length and sum to 1.
+ for (int i = 0; i < n; i++) {
+ assertSample(alphas.length, x[i]);
+ }
+
+ // Extreme concentrations have predictable output
+ final double minAlpha = Arrays.stream(alpha).min().getAsDouble();
+ final double maxAlpha = Arrays.stream(alpha).max().getAsDouble();
+
+ // If alpha is above 1e33 then the probability of a gamma sample
p(x=alpha) = 1.
+ // Check if the samples are all the same. The samples should be the
normalised
+ // log alpha values.
+ if (minAlpha > 1e33) {
+ final double[] expected =
Arrays.stream(alphas).map(Math::log).toArray();
+ Normalizer.LOG.normalize(expected);
+ for (int i = 0; i < n; i++) {
+ Assertions.assertArrayEquals(expected, x[i], "All samples
should be the same");
+ }
+ }
+
+ // If alpha is below 2^-53 then the log gamma samples are log(u) /
alpha.
+ // If alpha is all the same then the samples are expected to be 1.0 in
a single
+ // category and 0.0 elsewhere (assuming u output from the RNG is not
identical
+ // within a sample).
+ // If the alpha magnitude difference is very large across categories
some categories
+ // will never be observed.
+ if (maxAlpha < TINY_ALPHA) {
+ // Count observed categories
+ final long[] observed = new long[alphas.length];
+ if (symmetric) {
+ // A symmetric sample for a tiny concentration is to pick one
of the
+ // categories as 1.0 and all others 0.0.
+ // Note a previous check ensured the sum was 1.0. Here we test
for
+ // categories being 1.0 and all others zero, or allow a rare
case of 0.5
+ // in two categories (chance of 2^-53 for two p output from
the RNG to match).
+ for (int i = 0; i < n; i++) {
+ final double[] s = x[i];
+ int cOne = 0;
+ int cHalf = 0;
+ int cZero = 0;
+ for (int j = 0; j < s.length; j++) {
+ if (Math.abs(s[j] - 1.0) < 1e15) {
+ cOne++;
+ observed[j]++;
+ } else if (Math.abs(s[j] - 0.5) < 1e15) {
+ cHalf++;
+ observed[j]++;
+ } else if (Math.abs(s[j]) < 1e15) {
+ cZero++;
+ // not observed
+ }
+ }
+ // Expected 1.0, 0.5 or 0.0. Fail if anything else was
present.
+ Assertions.assertEquals(s.length, cOne + cHalf + cZero,
+ () -> Arrays.toString(s));
+ }
+ // The observed should follow a uniform distribution across
categories
+ final ChiSquareTest chiSquareTest = new ChiSquareTest();
+ final double[] expected = new double[alphas.length];
+ Arrays.fill(expected, 1.0 / alphas.length);
+ Assertions.assertFalse(chiSquareTest.chiSquareTest(expected,
observed, 0.001),
+ () -> "Symetric case was not uniform: " +
Arrays.toString(observed));
+ } else {
+ // Assume all alphas are different as the symmetric case was
not used.
+ // Flag any alpha not expected by comparison to max alpha.
+ // Identify the case where for all p in [0, 1):
+ // log(p1) / alpha1 > log(p2) / alpha2
+ // log(p1) / log(p2) > alpha1 / alpha2
+ // Extreme finite sample when p1 = 2^-53 and p2 = 1-2^-53:
+ // -36.74 / -1.11e-16 = 3.31e17
+ // If an alpha is much smaller than the max then it should not
be observed
+ // assuming p != 0 for the dominant concentration.
+ // This test can fail by chance so check the allowed category
across all samples.
+ final int[] allowed = Arrays.stream(alphas).mapToInt(a ->
maxAlpha > a * 4e17 ? 0 : 1).toArray();
+ for (int i = 0; i < n; i++) {
+ final double[] s = x[i];
+ for (int j = 0; j < s.length; j++) {
+ if (s[j] != 0) {
+ observed[j]++;
+ }
+ }
+ }
+ // Check the allowed categories
+ for (int j = 0; j < alphas.length; j++) {
+ // The disallowed category can occur by chance (frequency
2^-53 when the max
+ // alpha sample had p=0 from the RNG). This is very
unlikely multiple times.
+ if (allowed[j] == 0 && observed[j] > 1) {
+ Assertions.fail(String.format("%s : Unexpected
observation for tiny alpha %s",
+ Arrays.toString(alphas),
Arrays.toString(observed)));
+ }
+ }
+ }
+ }
+ }
+
@Test
void testSharedStateSampler() {
+ assertSharedStateSampler(1, 2, 3);
+ // Tiny alpha
+ assertSharedStateSampler(TINY_ALPHA / 1.5, TINY_ALPHA / 2, TINY_ALPHA
/ 2.5);
+ // Large alpha
+ assertSharedStateSampler(1e100, 2e100, 3e100);
+ // Small alpha (requires log samples) with a category with alpha=1
+ // that can use an exponential sampler
+ assertSharedStateSampler(1, 1e-5, 1e-3);
+ }
+
+ private static void assertSharedStateSampler(double... alpha) {
final RandomSource randomSource = RandomSource.XO_RO_SHI_RO_128_PP;
final byte[] seed = randomSource.createSeed();
final UniformRandomProvider rng1 = randomSource.create(seed);
final UniformRandomProvider rng2 = randomSource.create(seed);
- final DirichletSampler sampler1 = DirichletSampler.of(rng1, 1, 2, 3);
+ final DirichletSampler sampler1 = DirichletSampler.of(rng1, alpha);
final DirichletSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
@Test
void testSharedStateSamplerForSymmetricCase() {
+ assertSharedStateSamplerForSymmetricCase(2, 1.5);
+ // Tiny alpha
+ assertSharedStateSamplerForSymmetricCase(3, TINY_ALPHA / 1.5);
+ // Large alpha
+ assertSharedStateSamplerForSymmetricCase(2, 2.5e50);
+ // Exponential sampler
+ assertSharedStateSamplerForSymmetricCase(4, 1);
+ }
+
+ private static void assertSharedStateSamplerForSymmetricCase(int k, double
alpha) {
final RandomSource randomSource = RandomSource.XO_RO_SHI_RO_128_PP;
final byte[] seed = randomSource.createSeed();
final UniformRandomProvider rng1 = randomSource.create(seed);
final UniformRandomProvider rng2 = randomSource.create(seed);
- final DirichletSampler sampler1 = DirichletSampler.symmetric(rng1, 2,
1.5);
+ final DirichletSampler sampler1 = DirichletSampler.symmetric(rng1, k,
alpha);
final DirichletSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
@Test
void testSymmetricCaseMatchesGeneralCase() {
+ assertSymmetricCaseMatchesGeneralCase(3, 0.5);
+ assertSymmetricCaseMatchesGeneralCase(3, 1.0);
+ assertSymmetricCaseMatchesGeneralCase(3, 1.5);
+ // Tiny alpha
+ assertSymmetricCaseMatchesGeneralCase(4, TINY_ALPHA / 1.5);
+ // Large alpha
+ assertSymmetricCaseMatchesGeneralCase(2, 42e50);
+ }
+
+ private static void assertSymmetricCaseMatchesGeneralCase(int k, double
alpha) {
final RandomSource randomSource = RandomSource.XO_RO_SHI_RO_128_PP;
final byte[] seed = randomSource.createSeed();
final UniformRandomProvider rng1 = randomSource.create(seed);
final UniformRandomProvider rng2 = randomSource.create(seed);
- final int k = 3;
final double[] alphas = new double[k];
- for (final double alpha : new double[] {0.5, 1.0, 1.5}) {
- Arrays.fill(alphas, alpha);
- final DirichletSampler sampler1 = DirichletSampler.symmetric(rng1,
k, alpha);
- final DirichletSampler sampler2 = DirichletSampler.of(rng2,
alphas);
- RandomAssert.assertProduceSameSequence(sampler1, sampler2);
- }
+ Arrays.fill(alphas, alpha);
+ final DirichletSampler sampler1 = DirichletSampler.symmetric(rng1, k,
alpha);
+ final DirichletSampler sampler2 = DirichletSampler.of(rng2, alphas);
+ RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
/**
@@ -185,6 +536,30 @@ class DirichletSamplerTest {
assertSamples(1, 2, 3, 4);
}
+ @Test
+ void testSampling6() {
+ assertSamples(1000, 2000, 3000);
+ }
+
+ @Test
+ void testSampling7() {
+ // Log samples are used when max alpha < 2^-5 or > 2^53
+ assertSamples(0x1p-6, 0x2p-6, 0x3p-6);
+ }
+
+ @Test
+ void testSampling8() {
+ // Log samples are used when max alpha < 2^-5 or > 2^53
+ assertSamples(0x1p53, 0x2p53, 0x3p53);
+ }
+
+ @Test
+ void testSampling9() {
+ // tiny alpha.
+ // BigDecimal handles the small values to compute the mean and
variance.
+ assertSamples(TINY_ALPHA / 1.5, TINY_ALPHA / 2.5, TINY_ALPHA / 2);
+ }
+
/**
* Assert samples from the distribution. The variates are tested against
the expected
* mean and covariance for the given concentration parameters.
@@ -205,38 +580,56 @@ class DirichletSamplerTest {
// Computation of moments:
// https://en.wikipedia.org/wiki/Dirichlet_distribution#Moments
+ // Perform in extended precision to handle extreme concentration
parameters.
+ // 20 digits of precision when dividing is enough for a double result.
+ final MathContext mc = new MathContext(20);
// Compute the sum of the concentration parameters: alpha0
- double alpha0 = 0;
+ BigDecimal alpha0 = BigDecimal.ZERO;
for (int i = 0; i < k; i++) {
- alpha0 += alpha[i];
+ alpha0 = alpha0.add(bd(alpha[i]));
}
// Use a moderate tolerance.
// Differences are usually observed in the 3rd significant figure.
- final double relativeTolerance = 5e-2;
+ // When the alpha is large the covariance is very low and computation
of
+ // actual may be inaccurate. Here differences are in the 1st or
+ // 2nd significant figure.
+ final double meanRelativeTolerance = 5e-2;
+ final double varRelativeTolerance = alpha0.doubleValue() > 0x1p53 ?
0.1 : 5e-2;
// Mean = alpha[i] / alpha0
final double[] means = getColumnMeans(samples);
for (int i = 0; i < k; i++) {
- final double mean = alpha[i] / alpha0;
- Assertions.assertEquals(mean, means[i], mean * relativeTolerance,
"Mean");
+ final double mean = bd(alpha[i]).divide(alpha0, mc).doubleValue();
+ Assertions.assertEquals(mean, means[i], mean *
meanRelativeTolerance, "Mean");
}
// Variance = alpha_i (alpha_0 - alpha_i) / alpha_0^2 (alpha_0 + 1)
// Covariance = -alpha_i * alpha_j / alpha_0^2 (alpha_0 + 1)
final double[][] covars = getCovariance(samples);
- final double denom = alpha0 * alpha0 * (alpha0 + 1);
+ final BigDecimal denom =
alpha0.pow(2).multiply(alpha0.add(BigDecimal.ONE));
for (int i = 0; i < k; i++) {
- final double variance = alpha[i] * (alpha0 - alpha[i]) / denom;
- Assertions.assertEquals(variance, covars[i][i], variance *
relativeTolerance, "Variance");
+ final double variance =
bd(alpha[i]).multiply(alpha0.subtract(bd(alpha[i])))
+ .divide(denom, mc).doubleValue();
+ Assertions.assertEquals(variance, covars[i][i], variance *
varRelativeTolerance, "Variance");
for (int j = i + 1; j < k; j++) {
- final double covar = -alpha[i] * alpha[j] / denom;
- Assertions.assertEquals(covar, covars[i][j], Math.abs(covar) *
relativeTolerance, "Covariance");
+ final double covar =
bd(-alpha[i]).multiply(bd(alpha[j])).divide(denom, mc).doubleValue();
+ Assertions.assertEquals(covar, covars[i][j], Math.abs(covar) *
varRelativeTolerance, "Covariance");
}
}
}
+ /**
+ * Create a BigDecimal from the double value.
+ *
+ * @param x Value.
+ * @return the BigDecimal
+ */
+ private static BigDecimal bd(double x) {
+ return new BigDecimal(x);
+ }
+
/**
* Assert the sample has the correct length and sums to 1.
*
@@ -251,7 +644,7 @@ class DirichletSamplerTest {
for (int i = 2; i < x.length; i++) {
sum += x[i];
}
- Assertions.assertEquals(1.0, sum, 1e-10, "Invalid sum");
+ Assertions.assertEquals(1.0, sum, 2 * Math.ulp(1.0), "Invalid sum");
}
/**
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 606b6f4b..040bb0c6 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -56,6 +56,10 @@ 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="fix" due-to="Security scan, Alex Herbert"
issue="RNG-195">
+ "DirichletSampler": Avoid unbouded recursion for extremely small or
large
+ concentration parameters which created stack overflow errors.
+ </action>
<action dev="aherbert" type="add" due-to="Alex Herbert" issue="RNG-203">
"ZetaSampler": Add a sampler for the zeta distribution.
</action>
