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 7d0a4aa42d65acdc28d818ed4e6c8d44d81425f3
Author: Alex Herbert <aherb...@apache.org>
AuthorDate: Sun Oct 6 14:57:12 2019 +0100
RNG-121: ChengBetaSampler algorithms delegated to specialised classes.
Allows precomputation of factors including one Math.log() constant used
each loop iteration.
---
.../sampling/distribution/ChengBetaSampler.java | 387 +++++++++++++--------
.../distribution/ChengBetaSamplerTest.java | 53 ++-
src/changes/changes.xml | 4 +
3 files changed, 299 insertions(+), 145 deletions(-)
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
index f96e198..0e82f79 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
@@ -19,7 +19,7 @@ package org.apache.commons.rng.sampling.distribution;
import org.apache.commons.rng.UniformRandomProvider;
/**
- * Sampling from a <a
href="http://en.wikipedia.org/wiki/Beta_distribution";>Beta distribution</a>.
+ * Sampling from a <a
href="http://en.wikipedia.org/wiki/Beta_distribution";>beta distribution</a>.
*
* <p>Uses Cheng's algorithms for beta distribution sampling:</p>
*
@@ -38,181 +38,268 @@ import org.apache.commons.rng.UniformRandomProvider;
public class ChengBetaSampler
extends SamplerBase
implements SharedStateContinuousSampler {
- /** 1/2. */
- private static final double ONE_HALF = 1d / 2;
- /** 1/4. */
- private static final double ONE_QUARTER = 1d / 4;
-
- /** First shape parameter. */
- private final double alphaShape;
- /** Second shape parameter. */
- private final double betaShape;
- /** Underlying source of randomness. */
- private final UniformRandomProvider rng;
+ /** The appropriate beta sampler for the parameters. */
+ private final SharedStateContinuousSampler delegate;
/**
- * Creates a sampler instance.
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param alpha Distribution first shape parameter.
- * @param beta Distribution second shape parameter.
- * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta
<= 0}
+ * Base class to implement Cheng's algorithms for the beta distribution.
*/
- public ChengBetaSampler(UniformRandomProvider rng,
- double alpha,
- double beta) {
- super(null);
- if (alpha <= 0) {
- throw new IllegalArgumentException("alpha is not strictly
positive: " + alpha);
+ private abstract static class BaseChengBetaSampler
+ implements SharedStateContinuousSampler {
+ /**
+ * Flag set to true if {@code a} is the beta distribution alpha shape
parameter.
+ * Otherwise {@code a} is the beta distribution beta shape parameter.
+ */
+ protected final boolean aIsAlphaShape;
+ /**
+ * First shape parameter.
+ * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
+ */
+ protected final double a;
+ /**
+ * Second shape parameter.
+ * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
+ */
+ protected final double b;
+ /** Underlying source of randomness. */
+ protected final UniformRandomProvider rng;
+ /**
+ * The algorithm alpha factor. This is not the beta distribution alpha
shape parameter.
+ * It is the sum of the two shape parameters ({@code a + b}.
+ */
+ protected final double alpha;
+ /** The logarithm of the alpha factor. */
+ protected final double logAlpha;
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param aIsAlphaShape true if {@code a} is the beta distribution
alpha shape parameter.
+ * @param a Distribution first shape parameter.
+ * @param b Distribution second shape parameter.
+ */
+ BaseChengBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape,
double a, double b) {
+ this.rng = rng;
+ this.aIsAlphaShape = aIsAlphaShape;
+ this.a = a;
+ this.b = b;
+ alpha = a + b;
+ logAlpha = Math.log(alpha);
}
- if (beta <= 0) {
- throw new IllegalArgumentException("beta is not strictly positive:
" + beta);
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param source Source to copy.
+ */
+ private BaseChengBetaSampler(UniformRandomProvider rng,
+ BaseChengBetaSampler source) {
+ this.rng = rng;
+ aIsAlphaShape = source.aIsAlphaShape;
+ a = source.a;
+ b = source.b;
+ // Compute algorithm factors.
+ alpha = source.alpha;
+ logAlpha = source.logAlpha;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return "Cheng Beta deviate [" + rng.toString() + "]";
}
- this.rng = rng;
- alphaShape = alpha;
- betaShape = beta;
}
/**
- * @param rng Generator of uniformly distributed random numbers.
- * @param source Source to copy.
+ * Computes one sample using Cheng's BB algorithm, when beta distribution
{@code alpha} and
+ * {@code beta} shape parameters are both larger than 1.
*/
- private ChengBetaSampler(UniformRandomProvider rng,
- ChengBetaSampler source) {
- super(null);
- this.rng = rng;
- alphaShape = source.alphaShape;
- betaShape = source.betaShape;
- }
+ private static class ChengBBBetaSampler extends BaseChengBetaSampler {
+ /** The algorithm beta factor. This is not the beta distribution beta
shape parameter. */
+ private final double beta;
+ /** The algorithm gamma factor. */
+ private final double gamma;
- /** {@inheritDoc} */
- @Override
- public double sample() {
- final double a = Math.min(alphaShape, betaShape);
- final double b = Math.max(alphaShape, betaShape);
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param aIsAlphaShape true if {@code a} is the beta distribution
alpha shape parameter.
+ * @param a min(alpha, beta) shape parameter.
+ * @param b max(alpha, beta) shape parameter.
+ */
+ ChengBBBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape,
double a, double b) {
+ super(rng, aIsAlphaShape, a, b);
+ beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
+ gamma = a + 1 / beta;
+ }
- if (a > 1) {
- return algorithmBB(a, b);
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param source Source to copy.
+ */
+ private ChengBBBetaSampler(UniformRandomProvider rng,
+ ChengBBBetaSampler source) {
+ super(rng, source);
+ // Compute algorithm factors.
+ beta = source.beta;
+ gamma = source.gamma;
}
- // The algorithm is deliberately invoked with reversed parameters
- // as the argument order is: max(alpha,beta), min(alpha,beta).
- return algorithmBC(b, a);
- }
- /** {@inheritDoc} */
- @Override
- public String toString() {
- return "Cheng Beta deviate [" + rng.toString() + "]";
- }
+ @Override
+ public double sample() {
+ double r;
+ double w;
+ double t;
+ do {
+ final double u1 = rng.nextDouble();
+ final double u2 = rng.nextDouble();
+ final double v = beta * (Math.log(u1) - Math.log1p(-u1));
+ w = a * Math.exp(v);
+ final double z = u1 * u1 * u2;
+ r = gamma * v - 1.3862944;
+ final double s = a + r - w;
+ if (s + 2.609438 >= 5 * z) {
+ break;
+ }
- /** {@inheritDoc} */
- @Override
- public SharedStateContinuousSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- return new ChengBetaSampler(rng, this);
+ t = Math.log(z);
+ if (s >= t) {
+ break;
+ }
+ } while (r + alpha * (logAlpha - Math.log(b + w)) < t);
+
+ w = Math.min(w, Double.MAX_VALUE);
+
+ return aIsAlphaShape ? w / (b + w) : b / (b + w);
+ }
+
+ @Override
+ public SharedStateContinuousSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return new ChengBBBetaSampler(rng, this);
+ }
}
/**
- * Computes one sample using Cheng's BB algorithm, when \( \alpha \) and
- * \( \beta \) are both larger than 1.
- *
- * @param a \( \min(\alpha, \beta) \).
- * @param b \( \max(\alpha, \beta) \).
- * @return a random sample.
+ * Computes one sample using Cheng's BC algorithm, when at least one of
beta distribution
+ * {@code alpha} or {@code beta} shape parameters is smaller than 1.
*/
- private double algorithmBB(double a,
- double b) {
- final double alpha = a + b;
- final double beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
- final double gamma = a + 1 / beta;
-
- double r;
- double w;
- double t;
- do {
- final double u1 = rng.nextDouble();
- final double u2 = rng.nextDouble();
- final double v = beta * (Math.log(u1) - Math.log1p(-u1));
- w = a * Math.exp(v);
- final double z = u1 * u1 * u2;
- r = gamma * v - 1.3862944;
- final double s = a + r - w;
- if (s + 2.609438 >= 5 * z) {
- break;
- }
+ private static class ChengBCBetaSampler extends BaseChengBetaSampler {
+ /** 1/2. */
+ private static final double ONE_HALF = 1d / 2;
+ /** 1/4. */
+ private static final double ONE_QUARTER = 1d / 4;
- t = Math.log(z);
- if (s >= t) {
- break;
- }
- } while (r + alpha * (Math.log(alpha) - Math.log(b + w)) < t);
+ /** The algorithm beta factor. This is not the beta distribution beta
shape parameter. */
+ private final double beta;
+ /** The algorithm delta factor. */
+ private final double delta;
+ /** The algorithm k1 factor. */
+ private final double k1;
+ /** The algorithm k2 factor. */
+ private final double k2;
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param aIsAlphaShape true if {@code a} is the beta distribution
alpha shape parameter.
+ * @param a max(alpha, beta) shape parameter.
+ * @param b min(alpha, beta) shape parameter.
+ */
+ ChengBCBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape,
double a, double b) {
+ super(rng, aIsAlphaShape, a, b);
+ // Compute algorithm factors.
+ beta = 1 / b;
+ delta = 1 + a - b;
+ k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 0.777778);
+ k2 = 0.25 + (0.5 + 0.25 / delta) * b;
+ }
- w = Math.min(w, Double.MAX_VALUE);
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param source Source to copy.
+ */
+ private ChengBCBetaSampler(UniformRandomProvider rng,
+ ChengBCBetaSampler source) {
+ super(rng, source);
+ beta = source.beta;
+ delta = source.delta;
+ k1 = source.k1;
+ k2 = source.k2;
+ }
- return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
- }
+ @Override
+ public double sample() {
+ double w;
+ while (true) {
+ final double u1 = rng.nextDouble();
+ final double u2 = rng.nextDouble();
+ final double y = u1 * u2;
+ final double z = u1 * y;
+ if (u1 < ONE_HALF) {
+ if (0.25 * u2 + z - y >= k1) {
+ continue;
+ }
+ } else {
+ if (z <= ONE_QUARTER) {
+ final double v = beta * (Math.log(u1) -
Math.log1p(-u1));
+ w = a * Math.exp(v);
+ break;
+ }
- /**
- * Computes one sample using Cheng's BC algorithm, when at least one of
- * \( \alpha \) or \( \beta \) is smaller than 1.
- *
- * @param a \( \max(\alpha, \beta) \).
- * @param b \( \min(\alpha, \beta) \).
- * @return a random sample.
- */
- private double algorithmBC(double a,
- double b) {
- final double alpha = a + b;
- final double beta = 1 / b;
- final double delta = 1 + a - b;
- final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta -
0.777778);
- final double k2 = 0.25 + (0.5 + 0.25 / delta) * b;
-
- double w;
- while (true) {
- final double u1 = rng.nextDouble();
- final double u2 = rng.nextDouble();
- final double y = u1 * u2;
- final double z = u1 * y;
- if (u1 < ONE_HALF) {
- if (0.25 * u2 + z - y >= k1) {
- continue;
- }
- } else {
- if (z <= ONE_QUARTER) {
- final double v = beta * (Math.log(u1) - Math.log1p(-u1));
- w = a * Math.exp(v);
- break;
+ if (z >= k2) {
+ continue;
+ }
}
- if (z >= k2) {
- continue;
+ final double v = beta * (Math.log(u1) - Math.log1p(-u1));
+ w = a * Math.exp(v);
+ if (alpha * (logAlpha - Math.log(b + w) + v) - 1.3862944 >=
Math.log(z)) {
+ break;
}
}
- final double v = beta * (Math.log(u1) - Math.log1p(-u1));
- w = a * Math.exp(v);
- if (alpha * (Math.log(alpha) - Math.log(b + w) + v) - 1.3862944 >=
Math.log(z)) {
- break;
- }
- }
+ w = Math.min(w, Double.MAX_VALUE);
- w = Math.min(w, Double.MAX_VALUE);
+ return aIsAlphaShape ? w / (b + w) : b / (b + w);
+ }
- return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
+ @Override
+ public SharedStateContinuousSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return new ChengBCBetaSampler(rng, this);
+ }
}
/**
- * @param a Value.
- * @param b Value.
- * @return {@code true} if {@code a} is equal to {@code b}.
+ * Creates a sampler instance.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param alpha Distribution first shape parameter.
+ * @param beta Distribution second shape parameter.
+ * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta
<= 0}
*/
- private boolean equals(double a,
- double b) {
- return Math.abs(a - b) <= Double.MIN_VALUE;
+ public ChengBetaSampler(UniformRandomProvider rng,
+ double alpha,
+ double beta) {
+ super(null);
+ delegate = of(rng, alpha, beta);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ return delegate.sample();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return delegate.toString();
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public SharedStateContinuousSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return delegate.withUniformRandomProvider(rng);
}
/**
- * Creates a new Beta distribution sampler.
+ * Creates a new beta distribution sampler.
*
* @param rng Generator of uniformly distributed random numbers.
* @param alpha Distribution first shape parameter.
@@ -223,6 +310,24 @@ public class ChengBetaSampler
public static SharedStateContinuousSampler of(UniformRandomProvider rng,
double alpha,
double beta) {
- return new ChengBetaSampler(rng, alpha, beta);
+ if (alpha <= 0) {
+ throw new IllegalArgumentException("alpha is not strictly
positive: " + alpha);
+ }
+ if (beta <= 0) {
+ throw new IllegalArgumentException("beta is not strictly positive:
" + beta);
+ }
+
+ // Choose the algorithm.
+ final double a = Math.min(alpha, beta);
+ final double b = Math.max(alpha, beta);
+ final boolean aIsAlphaShape = a == alpha;
+
+ return a > 1 ?
+ // BB algorithm
+ new ChengBBBetaSampler(rng, aIsAlphaShape, a, b) :
+ // The BC algorithm is deliberately invoked with reversed
parameters
+ // as the argument order is: max(alpha,beta), min(alpha,beta).
+ // Also invert the 'a is alpha' flag.
+ new ChengBCBetaSampler(rng, !aIsAlphaShape, b, a);
}
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ChengBetaSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ChengBetaSamplerTest.java
index 8325244..0af2887 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ChengBetaSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ChengBetaSamplerTest.java
@@ -20,6 +20,7 @@ import org.apache.commons.rng.RestorableUniformRandomProvider;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.RandomAssert;
import org.apache.commons.rng.simple.RandomSource;
+import org.junit.Assert;
import org.junit.Test;
/**
@@ -54,14 +55,58 @@ public class ChengBetaSamplerTest {
* Test the SharedStateSampler implementation.
*/
@Test
- public void testSharedStateSampler() {
+ public void
testSharedStateSamplerWithAlphaAndBetaAbove1AndAlphaBelowBeta() {
+ testSharedStateSampler(1.23, 4.56);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation.
+ */
+ @Test
+ public void
testSharedStateSamplerWithAlphaAndBetaAbove1AndAlphaAboveBeta() {
+ testSharedStateSampler(4.56, 1.23);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation.
+ */
+ @Test
+ public void testSharedStateSamplerWithAlphaOrBetaBelow1AndAlphaBelowBeta()
{
+ testSharedStateSampler(0.23, 4.56);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation.
+ */
+ @Test
+ public void testSharedStateSamplerWithAlphaOrBetaBelow1AndAlphaAboveBeta()
{
+ testSharedStateSampler(4.56, 0.23);
+ }
+
+ /**
+ * Test the SharedStateSampler implementation.
+ *
+ * @param alpha Alpha.
+ * @param beta Beta.
+ */
+ private static void testSharedStateSampler(double alpha, double beta) {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
- final double alpha = 1.23;
- final double beta = 4.56;
+ // Use instance constructor not factory constructor to exercise 1.X
public API
final SharedStateContinuousSampler sampler1 =
- ChengBetaSampler.of(rng1, alpha, beta);
+ new ChengBetaSampler(rng1, alpha, beta);
final SharedStateContinuousSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
+
+ /**
+ * 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.
+ */
+ @Test
+ public void testToString() {
+ final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
+ Assert.assertTrue(new ChengBetaSampler(rng, 1.0, 2.0).toString()
+ .toLowerCase().contains("beta"));
+ }
}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 795e256..ae40c29 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -75,6 +75,10 @@ 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="RNG-121">
+ "ChengBetaSampler": Algorithms for different distribution parameters
have
+ been delegated to specialised classes.
+ </action>
<action dev="aherbert" type="update" issue="RNG-120">
Update security of serialization code for java.util.Random instances.
Implement
look-ahead deserialization or remove the use of
ObjectInputStream.readObject().
