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 b3a439bab6a71828b9d3426b37c895492de9a2ca
Author: Alex Herbert <aherb...@apache.org>
AuthorDate: Fri Jul 19 22:08:50 2019 +0100
RNG-110: Provide factory constructors for unreleased samplers.
These samplers have no requirement to maintain an instance constructor.
The constructor has been made private.
---
.../AliasMethodDiscreteSamplerPerformance.java | 4 +-
.../distribution/DiscreteSamplersPerformance.java | 12 +-
.../distribution/GeometricSamplersPerformance.java | 2 +-
.../distribution/PoissonSamplersPerformance.java | 2 +-
.../distribution/AliasMethodDiscreteSampler.java | 18 +-
.../sampling/distribution/GeometricSampler.java | 52 +-
.../distribution/GuideTableDiscreteSampler.java | 158 +--
.../distribution/KempSmallMeanPoissonSampler.java | 62 +-
.../distribution/LargeMeanPoissonSampler.java | 4 +-
.../MarsagliaTsangWangDiscreteSampler.java | 1063 ++++++++++----------
.../AliasMethodDiscreteSamplerTest.java | 14 +-
.../distribution/DiscreteSamplersList.java | 20 +-
.../distribution/GeometricSamplerTest.java | 18 +-
.../GuideTableDiscreteSamplerTest.java | 12 +-
.../KempSmallMeanPoissonSamplerTest.java | 16 +-
.../MarsagliaTsangWangDiscreteSamplerTest.java | 58 +-
src/main/resources/pmd/pmd-ruleset.xml | 3 +-
17 files changed, 778 insertions(+), 740 deletions(-)
diff --git
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/AliasMethodDiscreteSamplerPerformance.java
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/AliasMethodDiscreteSamplerPerformance.java
index 408e1ba..32d8307 100644
---
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/AliasMethodDiscreteSamplerPerformance.java
+++
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/AliasMethodDiscreteSamplerPerformance.java
@@ -96,7 +96,7 @@ public class AliasMethodDiscreteSamplerPerformance {
public void setup() {
probabilities = createProbabilities(size);
UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- sampler = AliasMethodDiscreteSampler.create(rng, probabilities,
alpha);
+ sampler = AliasMethodDiscreteSampler.of(rng, probabilities, alpha);
}
/**
@@ -155,6 +155,6 @@ public class AliasMethodDiscreteSamplerPerformance {
@Benchmark
public Object createSampler(DistributionData dist) {
// For the construction the RNG can be null
- return AliasMethodDiscreteSampler.create(null,
dist.getProbabilities(), dist.getAlpha());
+ return AliasMethodDiscreteSampler.of(null, dist.getProbabilities(),
dist.getAlpha());
}
}
diff --git
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/DiscreteSamplersPerformance.java
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/DiscreteSamplersPerformance.java
index c8cb81d..03bae18 100644
---
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/DiscreteSamplersPerformance.java
+++
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/DiscreteSamplersPerformance.java
@@ -119,17 +119,17 @@ public class DiscreteSamplersPerformance {
// Note: Use with a fractional part to the mean includes a
small mean sample
sampler = new LargeMeanPoissonSampler(rng, 41.7);
} else if ("GeometricSampler".equals(samplerType)) {
- sampler = new GeometricSampler(rng, 0.21);
+ sampler = GeometricSampler.of(rng, 0.21);
} else if
("MarsagliaTsangWangDiscreteSampler".equals(samplerType)) {
- sampler =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng,
DISCRETE_PROBABILITIES);
+ sampler = MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng,
DISCRETE_PROBABILITIES);
} else if ("MarsagliaTsangWangPoissonSampler".equals(samplerType))
{
- sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, 8.9);
+ sampler = MarsagliaTsangWangDiscreteSampler.Poisson.of(rng,
8.9);
} else if
("MarsagliaTsangWangBinomialSampler".equals(samplerType)) {
- sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, 20, 0.33);
+ sampler = MarsagliaTsangWangDiscreteSampler.Binomial.of(rng,
20, 0.33);
} else if ("GuideTableDiscreteSampler".equals(samplerType)) {
- sampler = new GuideTableDiscreteSampler(rng,
DISCRETE_PROBABILITIES);
+ sampler = GuideTableDiscreteSampler.of(rng,
DISCRETE_PROBABILITIES);
} else if ("AliasMethodDiscreteSampler".equals(samplerType)) {
- sampler = AliasMethodDiscreteSampler.create(rng,
DISCRETE_PROBABILITIES);
+ sampler = AliasMethodDiscreteSampler.of(rng,
DISCRETE_PROBABILITIES);
}
}
}
diff --git
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/GeometricSamplersPerformance.java
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/GeometricSamplersPerformance.java
index b87ef8f..7fc3e57 100644
---
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/GeometricSamplersPerformance.java
+++
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/GeometricSamplersPerformance.java
@@ -95,7 +95,7 @@ public class GeometricSamplersPerformance {
final RandomSource randomSource =
RandomSource.valueOf(randomSourceName);
final UniformRandomProvider rng =
RandomSource.create(randomSource);
if ("GeometricSampler".equals(samplerType)) {
- sampler = new GeometricSampler(rng, probabilityOfSuccess);
+ sampler = GeometricSampler.of(rng, probabilityOfSuccess);
} else {
final DiscreteInverseCumulativeProbabilityFunction
geometricFunction =
new
GeometricDiscreteInverseCumulativeProbabilityFunction(probabilityOfSuccess);
diff --git
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/PoissonSamplersPerformance.java
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/PoissonSamplersPerformance.java
index fbce71e..2ef6f8b 100644
---
a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/PoissonSamplersPerformance.java
+++
b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/distribution/PoissonSamplersPerformance.java
@@ -185,7 +185,7 @@ public class PoissonSamplersPerformance {
factory = new DiscreteSamplerFactory() {
@Override
public DiscreteSampler create() {
- return new KempSmallMeanPoissonSampler(generator,
mean);
+ return KempSmallMeanPoissonSampler.of(generator, mean);
}
};
} else if
("BoundedKempSmallMeanPoissonSampler".equals(samplerType)) {
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.java
index 2ad3593..2fa1ef7 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.java
@@ -52,7 +52,7 @@ import java.util.Arrays;
* approximately 12 bytes of storage per input probability, that is {@code n *
12} for size
* {@code n}. Zero-padding only requires 4 bytes of storage per padded value
as the probability is
* known to be zero. A table can be padded to a power of 2 using the utility
function
- * {@link #create(UniformRandomProvider, double[], int)} to construct the
sampler.</p>
+ * {@link #of(UniformRandomProvider, double[], int)} to construct the
sampler.</p>
*
* <p>An optimisation is performed for small table sizes that are a power of
2. In this case the
* sampling uses 1 or 2 calls from {@link UniformRandomProvider#nextInt()} to
generate up to
@@ -293,7 +293,7 @@ public class AliasMethodDiscreteSampler
* power-of-two. Padding is bounded by the upper limit on the size of an
array.</p>
*
* <p>To avoid zero-padding use the
- * {@link #create(UniformRandomProvider, double[], int)} method with a
negative
+ * {@link #of(UniformRandomProvider, double[], int)} method with a negative
* {@code alpha} factor.</p>
*
* @param rng Generator of uniformly distributed random numbers.
@@ -302,11 +302,11 @@ public class AliasMethodDiscreteSampler
* @throws IllegalArgumentException if {@code probabilities} is null or
empty, a
* probability is negative, infinite or {@code NaN}, or the sum of all
* probabilities is not strictly positive.
- * @see #create(UniformRandomProvider, double[], int)
+ * @see #of(UniformRandomProvider, double[], int)
*/
- public static AliasMethodDiscreteSampler create(final
UniformRandomProvider rng,
- final double[]
probabilities) {
- return create(rng, probabilities, DEFAULT_ALPHA);
+ public static SharedStateDiscreteSampler of(final UniformRandomProvider
rng,
+ final double[] probabilities) {
+ return of(rng, probabilities, DEFAULT_ALPHA);
}
/**
@@ -348,9 +348,9 @@ public class AliasMethodDiscreteSampler
* probability is negative, infinite or {@code NaN}, or the sum of all
* probabilities is not strictly positive.
*/
- public static AliasMethodDiscreteSampler create(final
UniformRandomProvider rng,
- final double[]
probabilities,
- int alpha) {
+ public static SharedStateDiscreteSampler of(final UniformRandomProvider
rng,
+ final double[] probabilities,
+ int alpha) {
// The Alias method balances N categories with counts around the mean
into N sections,
// each allocated 'mean' observations.
//
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GeometricSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GeometricSampler.java
index 25a4703..7cbca6b 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GeometricSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GeometricSampler.java
@@ -45,10 +45,7 @@ import org.apache.commons.rng.UniformRandomProvider;
*
* @since 1.3
*/
-public class GeometricSampler implements SharedStateDiscreteSampler {
- /** The appropriate geometric sampler for the parameters. */
- private final SharedStateDiscreteSampler delegate;
-
+public final class GeometricSampler {
/**
* Sample from the geometric distribution when the probability of success
is 1.
*/
@@ -132,48 +129,29 @@ public class GeometricSampler implements
SharedStateDiscreteSampler {
}
}
+ /** Class contains only static methods. */
+ private GeometricSampler() {}
+
/**
- * Creates a new geometric distribution sampler. The samples will be
provided in the set
- * {@code k=[0, 1, 2, ...]} where {@code k} indicates the number of
failures before the first
- * success.
+ * Creates a new geometric distribution sampler. The samples will be
provided in
+ * the set {@code k=[0, 1, 2, ...]} where {@code k} indicates the number of
+ * failures before the first success.
*
- * @param rng Generator of uniformly distributed random numbers
- * @param probabilityOfSuccess The probability of success
- * @throws IllegalArgumentException if {@code probabilityOfSuccess} is not
in the range
- * {@code [0 < probabilityOfSuccess <= 1]})
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param probabilityOfSuccess The probability of success.
+ * @return the sampler
+ * @throws IllegalArgumentException if {@code probabilityOfSuccess} is not
in
+ * the range {@code [0 < probabilityOfSuccess <= 1]})
*/
- public GeometricSampler(UniformRandomProvider rng, double
probabilityOfSuccess) {
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double probabilityOfSuccess) {
if (probabilityOfSuccess <= 0 || probabilityOfSuccess > 1) {
throw new IllegalArgumentException(
"Probability of success (p) must be in the range [0 < p <= 1]:
" +
probabilityOfSuccess);
}
- delegate = probabilityOfSuccess == 1 ?
+ return probabilityOfSuccess == 1 ?
GeometricP1Sampler.INSTANCE :
new GeometricExponentialSampler(rng, probabilityOfSuccess);
}
-
- /**
- * Create a sample from a geometric distribution.
- *
- * <p>The sample will take the values in the set {@code [0, 1, 2, ...]},
equivalent to the
- * number of failures before the first success.
- */
- @Override
- public int sample() {
- return delegate.sample();
- }
-
- /** {@inheritDoc} */
- @Override
- public String toString() {
- return delegate.toString();
- }
-
- /** {@inheritDoc} */
- @Override
- public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- // Direct return of the optimised sampler
- return delegate.withUniformRandomProvider(rng);
- }
}
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSampler.java
index 0ca3315..3ad4218 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSampler.java
@@ -19,27 +19,27 @@ package org.apache.commons.rng.sampling.distribution;
import org.apache.commons.rng.UniformRandomProvider;
/**
- * Compute a sample from a discrete probability distribution. The cumulative
probability
- * distribution is searched using a guide table to set an initial start point.
This implementation
- * is based on:
+ * Compute a sample from {@code n} values each with an associated probability.
If all unique items
+ * are assigned the same probability it is more efficient to use the {@link
DiscreteUniformSampler}.
*
- * <ul>
- * <li>
- * <blockquote>
- * Devroye, Luc (1986). Non-Uniform Random Variate Generation.
- * New York: Springer-Verlag. Chapter 3.2.4 "The method of guide tables" p.
96.
- * </blockquote>
- * </li>
- * </ul>
+ * <p>The cumulative probability distribution is searched using a guide table
to set an
+ * initial start point. This implementation is based on:</p>
+ *
+ * <blockquote>
+ * Devroye, Luc (1986). Non-Uniform Random Variate Generation.
+ * New York: Springer-Verlag. Chapter 3.2.4 "The method of guide tables" p.
96.
+ * </blockquote>
*
* <p>The size of the guide table can be controlled using a parameter. A
larger guide table
* will improve performance at the cost of storage space.</p>
*
* <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
*
+ * @see <a
href="http://en.wikipedia.org/wiki/Probability_distribution#Discrete_probability_distribution";>
+ * Discrete probability distribution (Wikipedia)</a>
* @since 1.3
*/
-public class GuideTableDiscreteSampler
+public final class GuideTableDiscreteSampler
implements SharedStateDiscreteSampler {
/** The default value for {@code alpha}. */
private static final double DEFAULT_ALPHA = 1.0;
@@ -63,38 +63,93 @@ public class GuideTableDiscreteSampler
private final int[] guideTable;
/**
- * Create a new instance using the default guide table size.
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param cumulativeProbabilities The cumulative probability table ({@code
f(x)}).
+ * @param guideTable The inverse cumulative probability guide table.
+ */
+ private GuideTableDiscreteSampler(UniformRandomProvider rng,
+ double[] cumulativeProbabilities,
+ int[] guideTable) {
+ this.rng = rng;
+ this.cumulativeProbabilities = cumulativeProbabilities;
+ this.guideTable = guideTable;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public int sample() {
+ // Compute a probability
+ final double u = rng.nextDouble();
+
+ // Initialise the search using the guide table to find an initial
guess.
+ // The table provides an upper bound on the sample (x+1) for a known
+ // cumulative probability (f(x)).
+ int x = guideTable[getGuideTableIndex(u, guideTable.length)];
+ // Search down.
+ // In the edge case where u is 1.0 then 'x' will be 1 outside the
range of the
+ // cumulative probability table and this will decrement to a valid
range.
+ // In the case where 'u' is mapped to the same guide table index as a
lower
+ // cumulative probability f(x) (due to rounding down) then this will
not decrement
+ // and return the exclusive upper bound (x+1).
+ while (x != 0 && u <= cumulativeProbabilities[x - 1]) {
+ x--;
+ }
+ return x;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return "Guide table deviate [" + rng.toString() + "]";
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return new GuideTableDiscreteSampler(rng, cumulativeProbabilities,
guideTable);
+ }
+
+ /**
+ * Create a new sampler for an enumerated distribution using the given
{@code probabilities}.
+ * The samples corresponding to each probability are assumed to be a
natural sequence
+ * starting at zero.
+ *
+ * <p>The size of the guide table is {@code probabilities.length}.</p>
*
* @param rng Generator of uniformly distributed random numbers.
* @param probabilities The probabilities.
+ * @return the sampler
* @throws IllegalArgumentException if {@code probabilities} is null or
empty, a
* probability is negative, infinite or {@code NaN}, or the sum of all
* probabilities is not strictly positive.
*/
- public GuideTableDiscreteSampler(UniformRandomProvider rng,
- double[] probabilities) {
- this(rng, probabilities, DEFAULT_ALPHA);
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double[] probabilities) {
+ return of(rng, probabilities, DEFAULT_ALPHA);
}
/**
- * Create a new instance.
+ * Create a new sampler for an enumerated distribution using the given
{@code probabilities}.
+ * The samples corresponding to each probability are assumed to be a
natural sequence
+ * starting at zero.
*
- * <p>The size of the guide table is {@code alpha * probabilities.length}.
+ * <p>The size of the guide table is {@code alpha *
probabilities.length}.</p>
*
* @param rng Generator of uniformly distributed random numbers.
* @param probabilities The probabilities.
* @param alpha The alpha factor used to set the guide table size.
+ * @return the sampler
* @throws IllegalArgumentException if {@code probabilities} is null or
empty, a
* probability is negative, infinite or {@code NaN}, the sum of all
* probabilities is not strictly positive, or {@code alpha} is not
strictly positive.
*/
- public GuideTableDiscreteSampler(UniformRandomProvider rng,
- double[] probabilities,
- double alpha) {
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double[] probabilities,
+ double alpha) {
validateParameters(probabilities, alpha);
final int size = probabilities.length;
- cumulativeProbabilities = new double[size];
+ final double[] cumulativeProbabilities = new double[size];
double sumProb = 0;
int count = 0;
@@ -110,12 +165,10 @@ public class GuideTableDiscreteSampler
throw new IllegalArgumentException("Invalid sum of probabilities:
" + sumProb);
}
- this.rng = rng;
-
// Note: The guide table is at least length 1. Compute the size
avoiding overflow
// in case (alpha * size) is too large.
final int guideTableSize = (int) Math.ceil(alpha * size);
- guideTable = new int[Math.max(guideTableSize, guideTableSize + 1)];
+ final int[] guideTable = new int[Math.max(guideTableSize,
guideTableSize + 1)];
// Compute and store cumulative probability.
for (int x = 0; x < size; x++) {
@@ -123,7 +176,8 @@ public class GuideTableDiscreteSampler
cumulativeProbabilities[x] = (norm < 1) ? norm : 1.0;
// Set the guide table value as an exclusive upper bound (x + 1)
- guideTable[getGuideTableIndex(cumulativeProbabilities[x])] = x + 1;
+ final int index = getGuideTableIndex(cumulativeProbabilities[x],
guideTable.length);
+ guideTable[index] = x + 1;
}
// Edge case for round-off
@@ -139,17 +193,8 @@ public class GuideTableDiscreteSampler
for (int i = 1; i < guideTable.length; i++) {
guideTable[i] = Math.max(guideTable[i - 1], guideTable[i]);
}
- }
- /**
- * @param rng Generator of uniformly distributed random numbers.
- * @param source Source to copy.
- */
- private GuideTableDiscreteSampler(UniformRandomProvider rng,
- GuideTableDiscreteSampler source) {
- this.rng = rng;
- cumulativeProbabilities = source.cumulativeProbabilities;
- guideTable = source.guideTable;
+ return new GuideTableDiscreteSampler(rng, cumulativeProbabilities,
guideTable);
}
/**
@@ -171,48 +216,15 @@ public class GuideTableDiscreteSampler
/**
* Gets the guide table index for the probability. This is obtained using
- * {@code p * (guideTable.length - 1)} so is inside the length of the
table.
+ * {@code p * (tableLength - 1)} so is inside the length of the table.
*
* @param p Cumulative probability.
+ * @param tableLength Table length.
* @return the guide table index.
*/
- private int getGuideTableIndex(double p) {
+ private static int getGuideTableIndex(double p, int tableLength) {
// Note: This is only ever called when p is in the range of the
cumulative
// probability table. So assume 0 <= p <= 1.
- return (int) (p * (guideTable.length - 1));
- }
-
- /** {@inheritDoc} */
- @Override
- public int sample() {
- // Compute a probability
- final double u = rng.nextDouble();
-
- // Initialise the search using the guide table to find an initial
guess.
- // The table provides an upper bound on the sample (x+1) for a known
- // cumulative probability (f(x)).
- int x = guideTable[getGuideTableIndex(u)];
- // Search down.
- // In the edge case where u is 1.0 then 'x' will be 1 outside the
range of the
- // cumulative probability table and this will decrement to a valid
range.
- // In the case where 'u' is mapped to the same guide table index as a
lower
- // cumulative probability f(x) (due to rounding down) then this will
not decrement
- // and return the exclusive upper bound (x+1).
- while (x != 0 && u <= cumulativeProbabilities[x - 1]) {
- x--;
- }
- return x;
- }
-
- /** {@inheritDoc} */
- @Override
- public String toString() {
- return "Guide table deviate [" + rng.toString() + "]";
- }
-
- /** {@inheritDoc} */
- @Override
- public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- return new GuideTableDiscreteSampler(rng, this);
+ return (int) (p * (tableLength - 1));
}
}
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSampler.java
index 7c4a09c..353651d 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSampler.java
@@ -41,11 +41,11 @@ import org.apache.commons.rng.UniformRandomProvider;
* <p>Sampling uses 1 call to {@link UniformRandomProvider#nextDouble()}. This
method provides
* an alternative to the {@link SmallMeanPoissonSampler} for slow generators
of {@code double}.</p>
*
- * @since 1.3
* @see <a href="https://www.jstor.org/stable/2346348";>Kemp, A.W. (1981) JRSS
Vol. 30, pp.
* 249-253</a>
+ * @since 1.3
*/
-public class KempSmallMeanPoissonSampler
+public final class KempSmallMeanPoissonSampler
implements SharedStateDiscreteSampler {
/** Underlying source of randomness. */
private final UniformRandomProvider rng;
@@ -61,35 +61,15 @@ public class KempSmallMeanPoissonSampler
/**
* @param rng Generator of uniformly distributed random numbers.
+ * @param p0 Probability of the Poisson sample {@code p(x=0)}.
* @param mean Mean.
- * @throws IllegalArgumentException if {@code mean <= 0} or
- * {@code Math.exp(-mean) == 0}.
- */
- public KempSmallMeanPoissonSampler(UniformRandomProvider rng,
- double mean) {
- if (mean <= 0) {
- throw new IllegalArgumentException("Mean is not strictly positive:
" + mean);
- }
-
- p0 = Math.exp(-mean);
- if (p0 > 0) {
- this.rng = rng;
- this.mean = mean;
- } else {
- // This catches the edge case of a NaN mean
- throw new IllegalArgumentException("No probability for mean " +
mean);
- }
- }
-
- /**
- * @param rng Generator of uniformly distributed random numbers.
- * @param source Source to copy.
*/
private KempSmallMeanPoissonSampler(UniformRandomProvider rng,
- KempSmallMeanPoissonSampler source) {
+ double p0,
+ double mean) {
this.rng = rng;
- p0 = source.p0;
- mean = source.mean;
+ this.p0 = p0;
+ this.mean = mean;
}
/** {@inheritDoc} */
@@ -129,6 +109,32 @@ public class KempSmallMeanPoissonSampler
/** {@inheritDoc} */
@Override
public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- return new KempSmallMeanPoissonSampler(rng, this);
+ return new KempSmallMeanPoissonSampler(rng, p0, mean);
+ }
+
+ /**
+ * Creates a new sampler for the Poisson distribution.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param mean Mean of the distribution.
+ * @return the sampler
+ * @throws IllegalArgumentException if {@code mean <= 0} or
+ * {@code Math.exp(-mean) == 0}.
+ */
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double mean) {
+ if (mean <= 0) {
+ throw new IllegalArgumentException("Mean is not strictly positive:
" + mean);
+ }
+
+ final double p0 = Math.exp(-mean);
+
+ // Probability must be positive. As mean increases then p(0) decreases.
+ if (p0 > 0) {
+ return new KempSmallMeanPoissonSampler(rng, p0, mean);
+ }
+
+ // This catches the edge case of a NaN mean
+ throw new IllegalArgumentException("No probability for mean: " + mean);
}
}
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
index 2a68e32..bc71e5b 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/LargeMeanPoissonSampler.java
@@ -155,7 +155,7 @@ public class LargeMeanPoissonSampler
final double lambdaFractional = mean - lambda;
smallMeanPoissonSampler = (lambdaFractional < Double.MIN_VALUE) ?
NO_SMALL_MEAN_POISSON_SAMPLER : // Not used.
- new KempSmallMeanPoissonSampler(rng, lambdaFractional);
+ KempSmallMeanPoissonSampler.of(rng, lambdaFractional);
}
/**
@@ -196,7 +196,7 @@ public class LargeMeanPoissonSampler
// The algorithm requires a Poisson sample from the remaining lambda
fraction.
smallMeanPoissonSampler = (lambdaFractional < Double.MIN_VALUE) ?
NO_SMALL_MEAN_POISSON_SAMPLER : // Not used.
- new KempSmallMeanPoissonSampler(rng, lambdaFractional);
+ KempSmallMeanPoissonSampler.of(rng, lambdaFractional);
}
/**
diff --git
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java
index 32bc12e..97b9228 100644
---
a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java
+++
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java
@@ -40,17 +40,16 @@ import org.apache.commons.rng.UniformRandomProvider;
* <p>The sampler supports the following distributions:</p>
*
* <ul>
- * <li>Any discrete probability distribution (probabilities must be provided)
+ * <li>Enumerated distribution (probabilities must be provided for each
sample)
* <li>Poisson distribution up to {@code mean = 1024}
* <li>Binomial distribution up to {@code trials = 65535}
* </ul>
*
- * @since 1.3
* @see <a href="http://dx.doi.org/10.18637/jss.v011.i03";>Margsglia, et al
(2004) JSS Vol.
* 11, Issue 3</a>
+ * @since 1.3
*/
-public abstract class MarsagliaTsangWangDiscreteSampler
- implements SharedStateDiscreteSampler {
+public final class MarsagliaTsangWangDiscreteSampler {
/** The value 2<sup>8</sup> as an {@code int}. */
private static final int INT_8 = 1 << 8;
/** The value 2<sup>16</sup> as an {@code int}. */
@@ -60,34 +59,6 @@ public abstract class MarsagliaTsangWangDiscreteSampler
/** The value 2<sup>31</sup> as a {@code double}. */
private static final double DOUBLE_31 = 1L << 31;
- /** The general name of any discrete probability distribution. */
- private static final String DISCRETE_NAME = "discrete";
- /** The name of the Poisson distribution. */
- private static final String POISSON_NAME = "Poisson";
- /** The name of the Binomial distribution. */
- private static final String BINOMIAL_NAME = "Binomial";
-
- /**
- * Upper bound on the mean for the Poisson distribution.
- *
- * <p>The original source code provided in Marsaglia, et al (2004) has no
explicit
- * limit but the code fails at mean >= 1941 as the transform to compute
p(x=mode)
- * produces infinity. Use a conservative limit of 1024.</p>
- */
- private static final double MAX_POISSON_MEAN = 1024;
- /**
- * The threshold for the mean of the Poisson distribution to switch the
method used
- * to compute the probabilities. This is taken from the example software
provided by
- * Marsaglia, et al (2004).
- */
- private static final double POISSON_MEAN_THRESHOLD = 21.4;
-
- /** Underlying source of randomness. */
- protected final UniformRandomProvider rng;
-
- /** The name of the distribution. */
- private final String distributionName;
-
//
=========================================================================
// Implementation note:
//
@@ -109,15 +80,55 @@ public abstract class MarsagliaTsangWangDiscreteSampler
// when provided via an array of probabilities and the Poisson and Binomial
// distributions for a restricted set of parameters. The restrictions are
// imposed by the requirement to compute the entire probability
distribution
- // from the controlling parameter(s) using a recursive method.
+ // from the controlling parameter(s) using a recursive method. Factory
+ // constructors return a SharedStateDiscreteSampler instance. Each
distribution
+ // type is contained in an inner class.
//
=========================================================================
/**
+ * The base class for Marsaglia-Tsang-Wang samplers.
+ */
+ private abstract static class AbstractMarsagliaTsangWangDiscreteSampler
+ implements SharedStateDiscreteSampler {
+ /** Underlying source of randomness. */
+ protected final UniformRandomProvider rng;
+
+ /** The name of the distribution. */
+ private final String distributionName;
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param distributionName Distribution name.
+ */
+ AbstractMarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
+ String distributionName) {
+ this.rng = rng;
+ this.distributionName = distributionName;
+ }
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param source Source to copy.
+ */
+ AbstractMarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
+
AbstractMarsagliaTsangWangDiscreteSampler source) {
+ this.rng = rng;
+ this.distributionName = source.distributionName;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String toString() {
+ return "Marsaglia Tsang Wang " + distributionName + " deviate [" +
rng.toString() + "]";
+ }
+ }
+
+ /**
* An implementation for the sample algorithm based on the decomposition
of the
* index in the range {@code [0,2^30)} into 5 base-64 digits with 8-bit
backing storage.
*/
private static class MarsagliaTsangWangBase64Int8DiscreteSampler
- extends MarsagliaTsangWangDiscreteSampler {
+ extends AbstractMarsagliaTsangWangDiscreteSampler {
/** The mask to convert a {@code byte} to an unsigned 8-bit integer. */
private static final int MASK = 0xff;
@@ -257,7 +268,7 @@ public abstract class MarsagliaTsangWangDiscreteSampler
* index in the range {@code [0,2^30)} into 5 base-64 digits with 16-bit
backing storage.
*/
private static class MarsagliaTsangWangBase64Int16DiscreteSampler
- extends MarsagliaTsangWangDiscreteSampler {
+ extends AbstractMarsagliaTsangWangDiscreteSampler {
/** The mask to convert a {@code byte} to an unsigned 16-bit integer.
*/
private static final int MASK = 0xffff;
@@ -397,7 +408,7 @@ public abstract class MarsagliaTsangWangDiscreteSampler
* index in the range {@code [0,2^30)} into 5 base-64 digits with 32-bit
backing storage.
*/
private static class MarsagliaTsangWangBase64Int32DiscreteSampler
- extends MarsagliaTsangWangDiscreteSampler {
+ extends AbstractMarsagliaTsangWangDiscreteSampler {
/** Limit for look-up table 1. */
private final int t1;
/** Limit for look-up table 2. */
@@ -528,108 +539,10 @@ public abstract class MarsagliaTsangWangDiscreteSampler
}
}
- /**
- * Return a fixed result for the Binomial distribution. This is a special
class to handle
- * an edge case of probability of success equal to 0 or 1.
- */
- private static class MarsagliaTsangWangFixedResultBinomialSampler
- extends MarsagliaTsangWangDiscreteSampler {
- /** The result. */
- private final int result;
-
- /**
- * @param result Result.
- */
- MarsagliaTsangWangFixedResultBinomialSampler(int result) {
- super(null, BINOMIAL_NAME);
- this.result = result;
- }
-
- @Override
- public int sample() {
- return result;
- }
-
- @Override
- public String toString() {
- return BINOMIAL_NAME + " deviate";
- }
-
- @Override
- public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- // No shared state
- return this;
- }
- }
-
- /**
- * Return an inversion result for the Binomial distribution. This assumes
the
- * following:
- *
- * <pre>
- * Binomial(n, p) = 1 - Binomial(n, 1 - p)
- * </pre>
- */
- private static class MarsagliaTsangWangInversionBinomialSampler
- extends MarsagliaTsangWangDiscreteSampler {
- /** The number of trials. */
- private final int trials;
- /** The Binomial distribution sampler. */
- private final SharedStateDiscreteSampler sampler;
-
- /**
- * @param trials Number of trials.
- * @param sampler Binomial distribution sampler.
- */
- MarsagliaTsangWangInversionBinomialSampler(int trials,
- SharedStateDiscreteSampler
sampler) {
- super(null, BINOMIAL_NAME);
- this.trials = trials;
- this.sampler = sampler;
- }
-
- @Override
- public int sample() {
- return trials - sampler.sample();
- }
-
- @Override
- public String toString() {
- return sampler.toString();
- }
-
- @Override
- public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
- return new MarsagliaTsangWangInversionBinomialSampler(this.trials,
- this.sampler.withUniformRandomProvider(rng));
- }
- }
-
- /**
- * @param rng Generator of uniformly distributed random numbers.
- * @param distributionName Distribution name.
- */
- MarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
- String distributionName) {
- this.rng = rng;
- this.distributionName = distributionName;
- }
- /**
- * @param rng Generator of uniformly distributed random numbers.
- * @param source Source to copy.
- */
- MarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
- MarsagliaTsangWangDiscreteSampler
source) {
- this.rng = rng;
- this.distributionName = source.distributionName;
- }
- /** {@inheritDoc} */
- @Override
- public String toString() {
- return "Marsaglia Tsang Wang " + distributionName + " deviate [" +
rng.toString() + "]";
- }
+ /** Class contains only static methods. */
+ private MarsagliaTsangWangDiscreteSampler() {}
/**
* Gets the k<sup>th</sup> base 64 digit of {@code m}.
@@ -643,6 +556,17 @@ public abstract class MarsagliaTsangWangDiscreteSampler
}
/**
+ * Convert the probability to an integer in the range [0,2^30]. This is
the numerator of
+ * a fraction with assumed denominator 2<sup>30</sup>.
+ *
+ * @param p Probability.
+ * @return the fraction numerator
+ */
+ private static int toUnsignedInt30(double p) {
+ return (int) (p * INT_30 + 0.5);
+ }
+
+ /**
* Create a new instance for probabilities {@code p(i)} where the sample
value {@code x} is
* {@code i + offset}.
*
@@ -655,10 +579,10 @@ public abstract class MarsagliaTsangWangDiscreteSampler
* @param offset The offset (must be positive).
* @return Sampler.
*/
- private static MarsagliaTsangWangDiscreteSampler
createSampler(UniformRandomProvider rng,
- String
distributionName,
- int[] prob,
- int offset)
{
+ private static SharedStateDiscreteSampler
createSampler(UniformRandomProvider rng,
+ String
distributionName,
+ int[] prob,
+ int offset) {
// Note: No argument checks for private method.
// Choose implementation based on the maximum index
@@ -673,426 +597,545 @@ public abstract class MarsagliaTsangWangDiscreteSampler
}
//
=========================================================================
- // The following factory methods are the public API to construct a sampler
for:
- // - Any discrete probability distribution (from provided double[]
probabilities)
+ // The following public classes provide factory methods to construct a
sampler for:
+ // - Enumerated probability distribution (from provided double[]
probabilities)
// - Poisson distribution for mean <= 1024
// - Binomial distribution for trials <= 65535
//
=========================================================================
/**
- * Creates a sampler for a given probability distribution.
- *
- * <p>The probabilities will be normalised using their sum. The only
requirement
- * is the sum is positive.</p>
- *
- * <p>The sum of the probabilities is normalised to 2<sup>30</sup>. Note
that
- * probabilities are adjusted to the nearest 1<sup>-30</sup> due to
round-off during
- * the normalisation conversion. Consequently any probability less than
2<sup>-31</sup>
- * will not be observed in samples.</p>
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param probabilities The list of probabilities.
- * @return Sampler.
- * @throws IllegalArgumentException if {@code probabilities} is null or
empty, a
- * probability is negative, infinite or {@code NaN}, or the sum of all
- * probabilities is not strictly positive.
+ * Create a sampler for an enumerated distribution of {@code n} values
each with an
+ * associated probability.
+ * The samples corresponding to each probability are assumed to be a
natural sequence
+ * starting at zero.
*/
- public static MarsagliaTsangWangDiscreteSampler
createDiscreteDistribution(UniformRandomProvider rng,
-
double[] probabilities) {
- return createSampler(rng, DISCRETE_NAME,
normaliseProbabilities(probabilities), 0);
- }
+ public static final class Enumerated {
+ /** The name of the enumerated probability distribution. */
+ private static final String ENUMERATED_NAME = "Enumerated";
- /**
- * Normalise the probabilities to integers that sum to 2<sup>30</sup>.
- *
- * @param probabilities The list of probabilities.
- * @return the normalised probabilities.
- * @throws IllegalArgumentException if {@code probabilities} is null or
empty, a
- * probability is negative, infinite or {@code NaN}, or the sum of all
- * probabilities is not strictly positive.
- */
- private static int[] normaliseProbabilities(double[] probabilities) {
- final double sumProb =
InternalUtils.validateProbabilities(probabilities);
-
- // Compute the normalisation: 2^30 / sum
- final double normalisation = INT_30 / sumProb;
- final int[] prob = new int[probabilities.length];
- int sum = 0;
- int max = 0;
- int mode = 0;
- for (int i = 0; i < prob.length; i++) {
- // Add 0.5 for rounding
- final int p = (int) (probabilities[i] * normalisation + 0.5);
- sum += p;
- // Find the mode (maximum probability)
- if (max < p) {
- max = p;
- mode = i;
- }
- prob[i] = p;
- }
+ /** Class contains only static methods. */
+ private Enumerated() {}
- // The sum must be >= 2^30.
- // Here just compensate the difference onto the highest probability.
- prob[mode] += INT_30 - sum;
+ /**
+ * Creates a sampler for an enumerated distribution of {@code n}
values each with an
+ * associated probability.
+ *
+ * <p>The probabilities will be normalised using their sum. The only
requirement
+ * is the sum is positive.</p>
+ *
+ * <p>The sum of the probabilities is normalised to 2<sup>30</sup>.
Note that
+ * probabilities are adjusted to the nearest 1<sup>-30</sup> due to
round-off during
+ * the normalisation conversion. Consequently any probability less
than 2<sup>-31</sup>
+ * will not be observed in samples.</p>
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param probabilities The list of probabilities.
+ * @return Sampler.
+ * @throws IllegalArgumentException if {@code probabilities} is null
or empty, a
+ * probability is negative, infinite or {@code NaN}, or the sum of all
+ * probabilities is not strictly positive.
+ */
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double[] probabilities) {
+ return createSampler(rng, ENUMERATED_NAME,
normaliseProbabilities(probabilities), 0);
+ }
- return prob;
- }
+ /**
+ * Normalise the probabilities to integers that sum to 2<sup>30</sup>.
+ *
+ * @param probabilities The list of probabilities.
+ * @return the normalised probabilities.
+ * @throws IllegalArgumentException if {@code probabilities} is null
or empty, a
+ * probability is negative, infinite or {@code NaN}, or the sum of all
+ * probabilities is not strictly positive.
+ */
+ private static int[] normaliseProbabilities(double[] probabilities) {
+ final double sumProb =
InternalUtils.validateProbabilities(probabilities);
+
+ // Compute the normalisation: 2^30 / sum
+ final double normalisation = INT_30 / sumProb;
+ final int[] prob = new int[probabilities.length];
+ int sum = 0;
+ int max = 0;
+ int mode = 0;
+ for (int i = 0; i < prob.length; i++) {
+ // Add 0.5 for rounding
+ final int p = (int) (probabilities[i] * normalisation + 0.5);
+ sum += p;
+ // Find the mode (maximum probability)
+ if (max < p) {
+ max = p;
+ mode = i;
+ }
+ prob[i] = p;
+ }
- /**
- * Creates a sampler for the Poisson distribution.
- *
- * <p>Any probability less than 2<sup>-31</sup> will not be observed in
samples.</p>
- *
- * <p>Storage requirements depend on the tabulated probability values.
Example storage
- * requirements are listed below.</p>
- *
- * <pre>
- * mean table size kB
- * 0.25 882 0.88
- * 0.5 1135 1.14
- * 1 1200 1.20
- * 2 1451 1.45
- * 4 1955 1.96
- * 8 2961 2.96
- * 16 4410 4.41
- * 32 6115 6.11
- * 64 8499 8.50
- * 128 11528 11.53
- * 256 15935 31.87
- * 512 20912 41.82
- * 1024 30614 61.23
- * </pre>
- *
- * <p>Note: Storage changes to 2 bytes per index between {@code mean=128}
and {@code mean=256}.</p>
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param mean Mean.
- * @return Sampler.
- * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean >
1024}.
- */
- public static MarsagliaTsangWangDiscreteSampler
createPoissonDistribution(UniformRandomProvider rng,
-
double mean) {
- validatePoissonDistributionParameters(mean);
-
- // Create the distribution either from X=0 or from X=mode when the
mean is high.
- return mean < POISSON_MEAN_THRESHOLD ?
- createPoissonDistributionFromX0(rng, mean) :
- createPoissonDistributionFromXMode(rng, mean);
- }
+ // The sum must be >= 2^30.
+ // Here just compensate the difference onto the highest
probability.
+ prob[mode] += INT_30 - sum;
- /**
- * Validate the Poisson distribution parameters.
- *
- * @param mean Mean.
- * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean >
1024}.
- */
- private static void validatePoissonDistributionParameters(double mean) {
- if (mean <= 0) {
- throw new IllegalArgumentException("mean is not strictly positive:
" + mean);
- }
- if (mean > MAX_POISSON_MEAN) {
- throw new IllegalArgumentException("mean " + mean + " > " +
MAX_POISSON_MEAN);
+ return prob;
}
}
/**
- * Creates the Poisson distribution by computing probabilities recursively
from {@code X=0}.
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param mean Mean.
- * @return Sampler.
+ * Create a sampler for the Poisson distribution.
*/
- private static MarsagliaTsangWangDiscreteSampler
createPoissonDistributionFromX0(
- UniformRandomProvider rng, double mean) {
- final double p0 = Math.exp(-mean);
-
- // Recursive update of Poisson probability until the value is too small
- // p(x + 1) = p(x) * mean / (x + 1)
- double p = p0;
- int i;
- for (i = 1; p * DOUBLE_31 >= 1; i++) {
- p *= mean / i;
- }
+ public static final class Poisson {
+ /** The name of the Poisson distribution. */
+ private static final String POISSON_NAME = "Poisson";
- // Probabilities are 30-bit integers, assumed denominator 2^30
- final int size = i - 1;
- final int[] prob = new int[size];
-
- p = p0;
- prob[0] = toUnsignedInt30(p);
- // The sum must exceed 2^30. In edges cases this is false due to
round-off.
- int sum = prob[0];
- for (i = 1; i < prob.length; i++) {
- p *= mean / i;
- prob[i] = toUnsignedInt30(p);
- sum += prob[i];
- }
+ /**
+ * Upper bound on the mean for the Poisson distribution.
+ *
+ * <p>The original source code provided in Marsaglia, et al (2004) has
no explicit
+ * limit but the code fails at mean >= 1941 as the transform to
compute p(x=mode)
+ * produces infinity. Use a conservative limit of 1024.</p>
+ */
- // If the sum is < 2^30 add the remaining sum to the mode
(floor(mean)).
- prob[(int) mean] += Math.max(0, INT_30 - sum);
+ private static final double MAX_MEAN = 1024;
+ /**
+ * The threshold for the mean of the Poisson distribution to switch
the method used
+ * to compute the probabilities. This is taken from the example
software provided by
+ * Marsaglia, et al (2004).
+ */
+ private static final double MEAN_THRESHOLD = 21.4;
- // Note: offset = 0
- return createSampler(rng, POISSON_NAME, prob, 0);
- }
+ /** Class contains only static methods. */
+ private Poisson() {}
- /**
- * Creates the Poisson distribution by computing probabilities recursively
upward and downward
- * from {@code X=mode}, the location of the largest p-value.
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param mean Mean.
- * @return Sampler.
- */
- private static MarsagliaTsangWangDiscreteSampler
createPoissonDistributionFromXMode(
- UniformRandomProvider rng, double mean) {
- // If mean >= 21.4, generate from largest p-value up, then largest
down.
- // The largest p-value will be at the mode (floor(mean)).
-
- // Find p(x=mode)
- final int mode = (int) mean;
- // This transform is stable until mean >= 1941 where p will result in
Infinity
- // before the divisor i is large enough to start reducing the product
(i.e. i > c).
- final double c = mean * Math.exp(-mean / mode);
- double p = 1.0;
- int i;
- for (i = 1; i <= mode; i++) {
- p *= c / i;
+ /**
+ * Creates a sampler for the Poisson distribution.
+ *
+ * <p>Any probability less than 2<sup>-31</sup> will not be observed
in samples.</p>
+ *
+ * <p>Storage requirements depend on the tabulated probability values.
Example storage
+ * requirements are listed below.</p>
+ *
+ * <pre>
+ * mean table size kB
+ * 0.25 882 0.88
+ * 0.5 1135 1.14
+ * 1 1200 1.20
+ * 2 1451 1.45
+ * 4 1955 1.96
+ * 8 2961 2.96
+ * 16 4410 4.41
+ * 32 6115 6.11
+ * 64 8499 8.50
+ * 128 11528 11.53
+ * 256 15935 31.87
+ * 512 20912 41.82
+ * 1024 30614 61.23
+ * </pre>
+ *
+ * <p>Note: Storage changes to 2 bytes per index between {@code
mean=128} and {@code mean=256}.</p>
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param mean Mean.
+ * @return Sampler.
+ * @throws IllegalArgumentException if {@code mean <= 0} or {@code
mean > 1024}.
+ */
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ double mean) {
+ validatePoissonDistributionParameters(mean);
+
+ // Create the distribution either from X=0 or from X=mode when the
mean is high.
+ return mean < MEAN_THRESHOLD ?
+ createPoissonDistributionFromX0(rng, mean) :
+ createPoissonDistributionFromXMode(rng, mean);
}
- final double pMode = p;
- // Find the upper limit using recursive computation of the p-value.
- // Note this will exit when i overflows to negative so no check on the
range
- for (i = mode + 1; p * DOUBLE_31 >= 1; i++) {
- p *= mean / i;
- }
- final int last = i - 2;
-
- // Find the lower limit using recursive computation of the p-value.
- p = pMode;
- int j = -1;
- for (i = mode - 1; i >= 0; i--) {
- p *= (i + 1) / mean;
- if (p * DOUBLE_31 < 1) {
- j = i;
- break;
+ /**
+ * Validate the Poisson distribution parameters.
+ *
+ * @param mean Mean.
+ * @throws IllegalArgumentException if {@code mean <= 0} or {@code
mean > 1024}.
+ */
+ private static void validatePoissonDistributionParameters(double mean)
{
+ if (mean <= 0) {
+ throw new IllegalArgumentException("mean is not strictly
positive: " + mean);
+ }
+ if (mean > MAX_MEAN) {
+ throw new IllegalArgumentException("mean " + mean + " > " +
MAX_MEAN);
}
}
- // Probabilities are 30-bit integers, assumed denominator 2^30.
- // This is the minimum sample value: prob[x - offset] = p(x)
- final int offset = j + 1;
- final int size = last - offset + 1;
- final int[] prob = new int[size];
-
- p = pMode;
- prob[mode - offset] = toUnsignedInt30(p);
- // The sum must exceed 2^30. In edges cases this is false due to
round-off.
- int sum = prob[mode - offset];
- // From mode to upper limit
- for (i = mode + 1; i <= last; i++) {
- p *= mean / i;
- prob[i - offset] = toUnsignedInt30(p);
- sum += prob[i - offset];
- }
- // From mode to lower limit
- p = pMode;
- for (i = mode - 1; i >= offset; i--) {
- p *= (i + 1) / mean;
- prob[i - offset] = toUnsignedInt30(p);
- sum += prob[i - offset];
+ /**
+ * Creates the Poisson distribution by computing probabilities
recursively from {@code X=0}.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param mean Mean.
+ * @return Sampler.
+ */
+ private static SharedStateDiscreteSampler
createPoissonDistributionFromX0(
+ UniformRandomProvider rng, double mean) {
+ final double p0 = Math.exp(-mean);
+
+ // Recursive update of Poisson probability until the value is too
small
+ // p(x + 1) = p(x) * mean / (x + 1)
+ double p = p0;
+ int i;
+ for (i = 1; p * DOUBLE_31 >= 1; i++) {
+ p *= mean / i;
+ }
+
+ // Probabilities are 30-bit integers, assumed denominator 2^30
+ final int size = i - 1;
+ final int[] prob = new int[size];
+
+ p = p0;
+ prob[0] = toUnsignedInt30(p);
+ // The sum must exceed 2^30. In edges cases this is false due to
round-off.
+ int sum = prob[0];
+ for (i = 1; i < prob.length; i++) {
+ p *= mean / i;
+ prob[i] = toUnsignedInt30(p);
+ sum += prob[i];
+ }
+
+ // If the sum is < 2^30 add the remaining sum to the mode
(floor(mean)).
+ prob[(int) mean] += Math.max(0, INT_30 - sum);
+
+ // Note: offset = 0
+ return createSampler(rng, POISSON_NAME, prob, 0);
}
- // If the sum is < 2^30 add the remaining sum to the mode.
- // If above 2^30 then the effect is truncation of the long tail of the
distribution.
- prob[mode - offset] += Math.max(0, INT_30 - sum);
+ /**
+ * Creates the Poisson distribution by computing probabilities
recursively upward and downward
+ * from {@code X=mode}, the location of the largest p-value.
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param mean Mean.
+ * @return Sampler.
+ */
+ private static SharedStateDiscreteSampler
createPoissonDistributionFromXMode(
+ UniformRandomProvider rng, double mean) {
+ // If mean >= 21.4, generate from largest p-value up, then largest
down.
+ // The largest p-value will be at the mode (floor(mean)).
+
+ // Find p(x=mode)
+ final int mode = (int) mean;
+ // This transform is stable until mean >= 1941 where p will result
in Infinity
+ // before the divisor i is large enough to start reducing the
product (i.e. i > c).
+ final double c = mean * Math.exp(-mean / mode);
+ double p = 1.0;
+ int i;
+ for (i = 1; i <= mode; i++) {
+ p *= c / i;
+ }
+ final double pMode = p;
+
+ // Find the upper limit using recursive computation of the p-value.
+ // Note this will exit when i overflows to negative so no check on
the range
+ for (i = mode + 1; p * DOUBLE_31 >= 1; i++) {
+ p *= mean / i;
+ }
+ final int last = i - 2;
+
+ // Find the lower limit using recursive computation of the p-value.
+ p = pMode;
+ int j = -1;
+ for (i = mode - 1; i >= 0; i--) {
+ p *= (i + 1) / mean;
+ if (p * DOUBLE_31 < 1) {
+ j = i;
+ break;
+ }
+ }
- return createSampler(rng, POISSON_NAME, prob, offset);
+ // Probabilities are 30-bit integers, assumed denominator 2^30.
+ // This is the minimum sample value: prob[x - offset] = p(x)
+ final int offset = j + 1;
+ final int size = last - offset + 1;
+ final int[] prob = new int[size];
+
+ p = pMode;
+ prob[mode - offset] = toUnsignedInt30(p);
+ // The sum must exceed 2^30. In edges cases this is false due to
round-off.
+ int sum = prob[mode - offset];
+ // From mode to upper limit
+ for (i = mode + 1; i <= last; i++) {
+ p *= mean / i;
+ prob[i - offset] = toUnsignedInt30(p);
+ sum += prob[i - offset];
+ }
+ // From mode to lower limit
+ p = pMode;
+ for (i = mode - 1; i >= offset; i--) {
+ p *= (i + 1) / mean;
+ prob[i - offset] = toUnsignedInt30(p);
+ sum += prob[i - offset];
+ }
+
+ // If the sum is < 2^30 add the remaining sum to the mode.
+ // If above 2^30 then the effect is truncation of the long tail of
the distribution.
+ prob[mode - offset] += Math.max(0, INT_30 - sum);
+
+ return createSampler(rng, POISSON_NAME, prob, offset);
+ }
}
/**
- * Creates a sampler for the Binomial distribution.
- *
- * <p>Any probability less than 2<sup>-31</sup> will not be observed in
samples.</p>
- *
- * <p>Storage requirements depend on the tabulated probability values.
Example storage
- * requirements are listed below (in kB).</p>
- *
- * <pre>
- * p
- * trials 0.5 0.1 0.01 0.001
- * 4 0.06 0.63 0.44 0.44
- * 16 0.69 1.14 0.76 0.44
- * 64 4.73 2.40 1.14 0.51
- * 256 8.63 5.17 1.89 0.82
- * 1024 31.12 9.45 3.34 0.89
- * </pre>
- *
- * <p>The method requires that the Binomial distribution probability at
{@code x=0} can be computed.
- * This will fail when {@code (1 - p)^trials == 0} which requires {@code
trials} to be large
- * and/or {@code p} to be small. In this case an exception is raised.</p>
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param trials Number of trials.
- * @param probabilityOfSuccess Probability of success (p).
- * @return Sampler.
- * @throws IllegalArgumentException if {@code trials < 0} or {@code trials
>= 2^16},
- * {@code p} is not in the range {@code [0-1]}, or the probability
distribution cannot
- * be computed.
+ * Create a sampler for the Binomial distribution.
*/
- public static MarsagliaTsangWangDiscreteSampler
createBinomialDistribution(UniformRandomProvider rng,
-
int trials,
-
double probabilityOfSuccess) {
- validateBinomialDistributionParameters(trials, probabilityOfSuccess);
-
- // Handle edge cases
- if (probabilityOfSuccess == 0) {
- return new MarsagliaTsangWangFixedResultBinomialSampler(0);
- }
- if (probabilityOfSuccess == 1) {
- return new MarsagliaTsangWangFixedResultBinomialSampler(trials);
- }
+ public static final class Binomial {
+ /** The name of the Binomial distribution. */
+ private static final String BINOMIAL_NAME = "Binomial";
+
+ /**
+ * Return a fixed result for the Binomial distribution. This is a
special class to handle
+ * an edge case of probability of success equal to 0 or 1.
+ */
+ private static class MarsagliaTsangWangFixedResultBinomialSampler
+ extends AbstractMarsagliaTsangWangDiscreteSampler {
+ /** The result. */
+ private final int result;
+
+ /**
+ * @param result Result.
+ */
+ MarsagliaTsangWangFixedResultBinomialSampler(int result) {
+ super(null, BINOMIAL_NAME);
+ this.result = result;
+ }
- // Check the supported size.
- if (trials >= INT_16) {
- throw new IllegalArgumentException("Unsupported number of trials:
" + trials);
+ @Override
+ public int sample() {
+ return result;
+ }
+
+ @Override
+ public String toString() {
+ return BINOMIAL_NAME + " deviate";
+ }
+
+ @Override
+ public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ // No shared state
+ return this;
+ }
}
- return createBinomialDistributionSampler(rng, trials,
probabilityOfSuccess);
- }
+ /**
+ * Return an inversion result for the Binomial distribution. This
assumes the
+ * following:
+ *
+ * <pre>
+ * Binomial(n, p) = 1 - Binomial(n, 1 - p)
+ * </pre>
+ */
+ private static class MarsagliaTsangWangInversionBinomialSampler
+ extends AbstractMarsagliaTsangWangDiscreteSampler {
+ /** The number of trials. */
+ private final int trials;
+ /** The Binomial distribution sampler. */
+ private final SharedStateDiscreteSampler sampler;
+
+ /**
+ * @param trials Number of trials.
+ * @param sampler Binomial distribution sampler.
+ */
+ MarsagliaTsangWangInversionBinomialSampler(int trials,
+
SharedStateDiscreteSampler sampler) {
+ super(null, BINOMIAL_NAME);
+ this.trials = trials;
+ this.sampler = sampler;
+ }
- /**
- * Validate the Binomial distribution parameters.
- *
- * @param trials Number of trials.
- * @param probabilityOfSuccess Probability of success (p).
- * @throws IllegalArgumentException if {@code trials < 0} or
- * {@code p} is not in the range {@code [0-1]}
- */
- private static void validateBinomialDistributionParameters(int trials,
double probabilityOfSuccess) {
- if (trials < 0) {
- throw new IllegalArgumentException("Trials is not positive: " +
trials);
+ @Override
+ public int sample() {
+ return trials - sampler.sample();
+ }
+
+ @Override
+ public String toString() {
+ return sampler.toString();
+ }
+
+ @Override
+ public SharedStateDiscreteSampler
withUniformRandomProvider(UniformRandomProvider rng) {
+ return new
MarsagliaTsangWangInversionBinomialSampler(this.trials,
+ this.sampler.withUniformRandomProvider(rng));
+ }
}
- if (probabilityOfSuccess < 0 || probabilityOfSuccess > 1) {
- throw new IllegalArgumentException("Probability is not in range
[0,1]: " + probabilityOfSuccess);
+
+ /** Class contains only static methods. */
+ private Binomial() {}
+
+ /**
+ * Creates a sampler for the Binomial distribution.
+ *
+ * <p>Any probability less than 2<sup>-31</sup> will not be observed
in samples.</p>
+ *
+ * <p>Storage requirements depend on the tabulated probability values.
Example storage
+ * requirements are listed below (in kB).</p>
+ *
+ * <pre>
+ * p
+ * trials 0.5 0.1 0.01 0.001
+ * 4 0.06 0.63 0.44 0.44
+ * 16 0.69 1.14 0.76 0.44
+ * 64 4.73 2.40 1.14 0.51
+ * 256 8.63 5.17 1.89 0.82
+ * 1024 31.12 9.45 3.34 0.89
+ * </pre>
+ *
+ * <p>The method requires that the Binomial distribution probability
at {@code x=0} can be computed.
+ * This will fail when {@code (1 - p)^trials == 0} which requires
{@code trials} to be large
+ * and/or {@code p} to be small. In this case an exception is
raised.</p>
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param trials Number of trials.
+ * @param probabilityOfSuccess Probability of success (p).
+ * @return Sampler.
+ * @throws IllegalArgumentException if {@code trials < 0} or {@code
trials >= 2^16},
+ * {@code p} is not in the range {@code [0-1]}, or the probability
distribution cannot
+ * be computed.
+ */
+ public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
+ int trials,
+ double
probabilityOfSuccess) {
+ validateBinomialDistributionParameters(trials,
probabilityOfSuccess);
+
+ // Handle edge cases
+ if (probabilityOfSuccess == 0) {
+ return new MarsagliaTsangWangFixedResultBinomialSampler(0);
+ }
+ if (probabilityOfSuccess == 1) {
+ return new
MarsagliaTsangWangFixedResultBinomialSampler(trials);
+ }
+
+ // Check the supported size.
+ if (trials >= INT_16) {
+ throw new IllegalArgumentException("Unsupported number of
trials: " + trials);
+ }
+
+ return createBinomialDistributionSampler(rng, trials,
probabilityOfSuccess);
}
- }
- /**
- * Creates the Binomial distribution sampler.
- *
- * <p>This assumes the parameters for the distribution are valid. The
method
- * will only fail if the initial probability for {@code X=0} is zero.</p>
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param trials Number of trials.
- * @param probabilityOfSuccess Probability of success (p).
- * @return Sampler.
- * @throws IllegalArgumentException if the probability distribution cannot
be
- * computed.
- */
- private static MarsagliaTsangWangDiscreteSampler
createBinomialDistributionSampler(
- UniformRandomProvider rng, int trials, double
probabilityOfSuccess) {
-
- // The maximum supported value for Math.exp is approximately -744.
- // This occurs when trials is large and p is close to 1.
- // Handle this by using an inversion: generate j=Binomial(n,1-p),
return n-j
- final boolean useInversion = probabilityOfSuccess > 0.5;
- final double p = useInversion ? 1 - probabilityOfSuccess :
probabilityOfSuccess;
-
- // Check if the distribution can be computed
- final double p0 = Math.exp(trials * Math.log(1 - p));
- if (p0 < Double.MIN_VALUE) {
- throw new IllegalArgumentException("Unable to compute
distribution");
+ /**
+ * Validate the Binomial distribution parameters.
+ *
+ * @param trials Number of trials.
+ * @param probabilityOfSuccess Probability of success (p).
+ * @throws IllegalArgumentException if {@code trials < 0} or
+ * {@code p} is not in the range {@code [0-1]}
+ */
+ private static void validateBinomialDistributionParameters(int trials,
double probabilityOfSuccess) {
+ if (trials < 0) {
+ throw new IllegalArgumentException("Trials is not positive: "
+ trials);
+ }
+ if (probabilityOfSuccess < 0 || probabilityOfSuccess > 1) {
+ throw new IllegalArgumentException("Probability is not in
range [0,1]: " + probabilityOfSuccess);
+ }
}
- // First find size of probability array
- double t = p0;
- final double h = p / (1 - p);
- // Find first probability above the threshold of 2^-31
- int begin = 0;
- if (t * DOUBLE_31 < 1) {
- // Somewhere after p(0)
- // Note:
- // If this loop is entered p(0) is < 2^-31.
- // This has been tested at the extreme for p(0)=Double.MIN_VALUE
and either
- // p=0.5 or trials=2^16-1 and does not fail to find the beginning.
- for (int i = 1; i <= trials; i++) {
+ /**
+ * Creates the Binomial distribution sampler.
+ *
+ * <p>This assumes the parameters for the distribution are valid. The
method
+ * will only fail if the initial probability for {@code X=0} is
zero.</p>
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param trials Number of trials.
+ * @param probabilityOfSuccess Probability of success (p).
+ * @return Sampler.
+ * @throws IllegalArgumentException if the probability distribution
cannot be
+ * computed.
+ */
+ private static SharedStateDiscreteSampler
createBinomialDistributionSampler(
+ UniformRandomProvider rng, int trials, double
probabilityOfSuccess) {
+
+ // The maximum supported value for Math.exp is approximately -744.
+ // This occurs when trials is large and p is close to 1.
+ // Handle this by using an inversion: generate j=Binomial(n,1-p),
return n-j
+ final boolean useInversion = probabilityOfSuccess > 0.5;
+ final double p = useInversion ? 1 - probabilityOfSuccess :
probabilityOfSuccess;
+
+ // Check if the distribution can be computed
+ final double p0 = Math.exp(trials * Math.log(1 - p));
+ if (p0 < Double.MIN_VALUE) {
+ throw new IllegalArgumentException("Unable to compute
distribution");
+ }
+
+ // First find size of probability array
+ double t = p0;
+ final double h = p / (1 - p);
+ // Find first probability above the threshold of 2^-31
+ int begin = 0;
+ if (t * DOUBLE_31 < 1) {
+ // Somewhere after p(0)
+ // Note:
+ // If this loop is entered p(0) is < 2^-31.
+ // This has been tested at the extreme for
p(0)=Double.MIN_VALUE and either
+ // p=0.5 or trials=2^16-1 and does not fail to find the
beginning.
+ for (int i = 1; i <= trials; i++) {
+ t *= (trials + 1 - i) * h / i;
+ if (t * DOUBLE_31 >= 1) {
+ begin = i;
+ break;
+ }
+ }
+ }
+ // Find last probability
+ int end = trials;
+ for (int i = begin + 1; i <= trials; i++) {
t *= (trials + 1 - i) * h / i;
- if (t * DOUBLE_31 >= 1) {
- begin = i;
+ if (t * DOUBLE_31 < 1) {
+ end = i - 1;
break;
}
}
- }
- // Find last probability
- int end = trials;
- for (int i = begin + 1; i <= trials; i++) {
- t *= (trials + 1 - i) * h / i;
- if (t * DOUBLE_31 < 1) {
- end = i - 1;
- break;
- }
- }
-
- return createBinomialDistributionSamplerFromRange(rng, trials, p,
useInversion,
- p0, begin, end);
- }
- /**
- * Creates the Binomial distribution sampler using only the probability
values for {@code X}
- * between the begin and the end (inclusive).
- *
- * @param rng Generator of uniformly distributed random numbers.
- * @param trials Number of trials.
- * @param p Probability of success (p).
- * @param useInversion Set to {@code true} if the probability was inverted.
- * @param p0 Probability at {@code X=0}
- * @param begin Begin value {@code X} for the distribution.
- * @param end End value {@code X} for the distribution.
- * @return Sampler.
- */
- private static MarsagliaTsangWangDiscreteSampler
createBinomialDistributionSamplerFromRange(
- UniformRandomProvider rng, int trials, double p,
- boolean useInversion, double p0, int begin, int end) {
-
- // Assign probability values as 30-bit integers
- final int size = end - begin + 1;
- final int[] prob = new int[size];
- double t = p0;
- final double h = p / (1 - p);
- for (int i = 1; i <= begin; i++) {
- t *= (trials + 1 - i) * h / i;
- }
- int sum = toUnsignedInt30(t);
- prob[0] = sum;
- for (int i = begin + 1; i <= end; i++) {
- t *= (trials + 1 - i) * h / i;
- prob[i - begin] = toUnsignedInt30(t);
- sum += prob[i - begin];
+ return createBinomialDistributionSamplerFromRange(rng, trials, p,
useInversion,
+ p0, begin, end);
}
- // If the sum is < 2^30 add the remaining sum to the mode
(floor((n+1)p))).
- // If above 2^30 then the effect is truncation of the long tail of the
distribution.
- final int mode = (int) ((trials + 1) * p) - begin;
- prob[mode] += Math.max(0, INT_30 - sum);
+ /**
+ * Creates the Binomial distribution sampler using only the
probability values for {@code X}
+ * between the begin and the end (inclusive).
+ *
+ * @param rng Generator of uniformly distributed random numbers.
+ * @param trials Number of trials.
+ * @param p Probability of success (p).
+ * @param useInversion Set to {@code true} if the probability was
inverted.
+ * @param p0 Probability at {@code X=0}
+ * @param begin Begin value {@code X} for the distribution.
+ * @param end End value {@code X} for the distribution.
+ * @return Sampler.
+ */
+ private static SharedStateDiscreteSampler
createBinomialDistributionSamplerFromRange(
+ UniformRandomProvider rng, int trials, double p,
+ boolean useInversion, double p0, int begin, int end) {
+
+ // Assign probability values as 30-bit integers
+ final int size = end - begin + 1;
+ final int[] prob = new int[size];
+ double t = p0;
+ final double h = p / (1 - p);
+ for (int i = 1; i <= begin; i++) {
+ t *= (trials + 1 - i) * h / i;
+ }
+ int sum = toUnsignedInt30(t);
+ prob[0] = sum;
+ for (int i = begin + 1; i <= end; i++) {
+ t *= (trials + 1 - i) * h / i;
+ prob[i - begin] = toUnsignedInt30(t);
+ sum += prob[i - begin];
+ }
- final MarsagliaTsangWangDiscreteSampler sampler = createSampler(rng,
BINOMIAL_NAME, prob, begin);
+ // If the sum is < 2^30 add the remaining sum to the mode
(floor((n+1)p))).
+ // If above 2^30 then the effect is truncation of the long tail of
the distribution.
+ final int mode = (int) ((trials + 1) * p) - begin;
+ prob[mode] += Math.max(0, INT_30 - sum);
- // Check if an inversion was made
- return useInversion ?
- new MarsagliaTsangWangInversionBinomialSampler(trials, sampler)
:
- sampler;
- }
+ final SharedStateDiscreteSampler sampler = createSampler(rng,
BINOMIAL_NAME, prob, begin);
- /**
- * Convert the probability to an integer in the range [0,2^30]. This is
the numerator of
- * a fraction with assumed denominator 2<sup>30</sup>.
- *
- * @param p Probability.
- * @return the fraction numerator
- */
- private static int toUnsignedInt30(double p) {
- return (int) (p * INT_30 + 0.5);
+ // Check if an inversion was made
+ return useInversion ?
+ new MarsagliaTsangWangInversionBinomialSampler(trials,
sampler) :
+ sampler;
+ }
}
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSamplerTest.java
index e5f9108..0264c7c 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSamplerTest.java
@@ -68,7 +68,7 @@ public class AliasMethodDiscreteSamplerTest {
@Test
public void testToString() {
- final AliasMethodDiscreteSampler sampler = createSampler(new double[]
{0.5, 0.5});
+ final SharedStateDiscreteSampler sampler = createSampler(new double[]
{0.5, 0.5});
Assert.assertTrue(sampler.toString().toLowerCase().contains("alias
method"));
}
@@ -78,9 +78,9 @@ public class AliasMethodDiscreteSamplerTest {
* @param probabilities the probabilities
* @return the alias method discrete sampler
*/
- private static AliasMethodDiscreteSampler createSampler(double[]
probabilities) {
+ private static SharedStateDiscreteSampler createSampler(double[]
probabilities) {
final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- return AliasMethodDiscreteSampler.create(rng, probabilities, -1);
+ return AliasMethodDiscreteSampler.of(rng, probabilities, -1);
}
/**
@@ -130,7 +130,7 @@ public class AliasMethodDiscreteSamplerTest {
@Test
public void
testNonUniformSamplesWithProbabilitiesWithDefaultFactoryConstructor() {
final double[] expected = {0.1, 0.2, 0.3, 0.1, 0.3 };
-
checkSamples(AliasMethodDiscreteSampler.create(RandomSource.create(RandomSource.SPLIT_MIX_64),
expected), expected);
+
checkSamples(AliasMethodDiscreteSampler.of(RandomSource.create(RandomSource.SPLIT_MIX_64),
expected), expected);
}
/**
@@ -218,7 +218,7 @@ public class AliasMethodDiscreteSamplerTest {
*
* @param expected the expected probabilities
*/
- private static void checkSamples(AliasMethodDiscreteSampler sampler,
double[] probabilies) {
+ private static void checkSamples(SharedStateDiscreteSampler sampler,
double[] probabilies) {
final int numberOfSamples = 10000;
final long[] samples = new long[probabilies.length];
for (int i = 0; i < numberOfSamples; i++) {
@@ -275,8 +275,8 @@ public class AliasMethodDiscreteSamplerTest {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
// Use negative alpha to disable padding
- final AliasMethodDiscreteSampler sampler1 =
- AliasMethodDiscreteSampler.create(rng1, probabilities, -1);
+ final SharedStateDiscreteSampler sampler1 =
+ AliasMethodDiscreteSampler.of(rng1, probabilities, -1);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
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 f6fe40e..73c7cb7 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
@@ -53,12 +53,12 @@ public final class DiscreteSamplersList {
add(LIST, new
org.apache.commons.math3.distribution.BinomialDistribution(unusedRng,
trialsBinomial, probSuccessBinomial),
// range [9,16]
MathArrays.sequence(8, 9, 1),
-
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(RandomSource.create(RandomSource.WELL_19937_A),
trialsBinomial, probSuccessBinomial));
+
MarsagliaTsangWangDiscreteSampler.Binomial.of(RandomSource.create(RandomSource.WELL_19937_A),
trialsBinomial, probSuccessBinomial));
// Inverted
add(LIST, new
org.apache.commons.math3.distribution.BinomialDistribution(unusedRng,
trialsBinomial, 1 - probSuccessBinomial),
// range [4,11] = [20-16, 20-9]
MathArrays.sequence(8, 4, 1),
-
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(RandomSource.create(RandomSource.WELL_19937_C),
trialsBinomial, 1 - probSuccessBinomial));
+
MarsagliaTsangWangDiscreteSampler.Binomial.of(RandomSource.create(RandomSource.WELL_19937_C),
trialsBinomial, 1 - probSuccessBinomial));
// Geometric ("inverse method").
final double probSuccessGeometric = 0.21;
@@ -68,7 +68,7 @@ public final class DiscreteSamplersList {
// Geometric.
add(LIST, new
org.apache.commons.math3.distribution.GeometricDistribution(unusedRng,
probSuccessGeometric),
MathArrays.sequence(10, 0, 1),
- new
GeometricSampler(RandomSource.create(RandomSource.XOR_SHIFT_1024_S),
probSuccessGeometric));
+
GeometricSampler.of(RandomSource.create(RandomSource.XOR_SHIFT_1024_S),
probSuccessGeometric));
// Hypergeometric ("inverse method").
final int popSizeHyper = 34;
@@ -136,10 +136,10 @@ public final class DiscreteSamplersList {
new
SmallMeanPoissonSampler(RandomSource.create(RandomSource.XO_SHI_RO_256_PLUS),
meanPoisson));
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
meanPoisson, epsilonPoisson, maxIterationsPoisson),
MathArrays.sequence(10, 0, 1),
- new
KempSmallMeanPoissonSampler(RandomSource.create(RandomSource.XO_SHI_RO_128_PLUS),
meanPoisson));
+
KempSmallMeanPoissonSampler.of(RandomSource.create(RandomSource.XO_SHI_RO_128_PLUS),
meanPoisson));
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
meanPoisson, epsilonPoisson, maxIterationsPoisson),
MathArrays.sequence(10, 0, 1),
-
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(RandomSource.create(RandomSource.XO_SHI_RO_128_PLUS),
meanPoisson));
+
MarsagliaTsangWangDiscreteSampler.Poisson.of(RandomSource.create(RandomSource.XO_SHI_RO_128_PLUS),
meanPoisson));
// Poisson (40 < mean < 80).
final double largeMeanPoisson = 67.89;
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
largeMeanPoisson, epsilonPoisson, maxIterationsPoisson),
@@ -151,7 +151,7 @@ public final class DiscreteSamplersList {
new
LargeMeanPoissonSampler(RandomSource.create(RandomSource.SPLIT_MIX_64),
largeMeanPoisson));
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
largeMeanPoisson, epsilonPoisson, maxIterationsPoisson),
MathArrays.sequence(50, (int) (largeMeanPoisson - 25), 1),
-
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PLUS),
largeMeanPoisson));
+
MarsagliaTsangWangDiscreteSampler.Poisson.of(RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PLUS),
largeMeanPoisson));
// Poisson (mean >> 40).
final double veryLargeMeanPoisson = 543.21;
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
veryLargeMeanPoisson, epsilonPoisson, maxIterationsPoisson),
@@ -163,16 +163,16 @@ public final class DiscreteSamplersList {
new
LargeMeanPoissonSampler(RandomSource.create(RandomSource.SPLIT_MIX_64),
veryLargeMeanPoisson));
add(LIST, new
org.apache.commons.math3.distribution.PoissonDistribution(unusedRng,
veryLargeMeanPoisson, epsilonPoisson, maxIterationsPoisson),
MathArrays.sequence(100, (int) (veryLargeMeanPoisson - 50), 1),
-
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(RandomSource.create(RandomSource.XO_RO_SHI_RO_64_SS),
veryLargeMeanPoisson));
+
MarsagliaTsangWangDiscreteSampler.Poisson.of(RandomSource.create(RandomSource.XO_RO_SHI_RO_64_SS),
veryLargeMeanPoisson));
// Any discrete distribution
final double[] discreteProbabilities = new double[] {0.1, 0.2,
0.3, 0.4, 0.5};
add(LIST, discreteProbabilities,
-
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(RandomSource.create(RandomSource.XO_SHI_RO_512_PLUS),
discreteProbabilities));
+
MarsagliaTsangWangDiscreteSampler.Enumerated.of(RandomSource.create(RandomSource.XO_SHI_RO_512_PLUS),
discreteProbabilities));
add(LIST, discreteProbabilities,
- new
GuideTableDiscreteSampler(RandomSource.create(RandomSource.XO_SHI_RO_512_SS),
discreteProbabilities));
+
GuideTableDiscreteSampler.of(RandomSource.create(RandomSource.XO_SHI_RO_512_SS),
discreteProbabilities));
add(LIST, discreteProbabilities,
-
AliasMethodDiscreteSampler.create(RandomSource.create(RandomSource.KISS),
discreteProbabilities));
+
AliasMethodDiscreteSampler.of(RandomSource.create(RandomSource.KISS),
discreteProbabilities));
} catch (Exception e) {
// CHECKSTYLE: stop Regexp
System.err.println("Unexpected exception while creating the list
of samplers: " + e);
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GeometricSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GeometricSamplerTest.java
index f384485..baaa77f 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GeometricSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GeometricSamplerTest.java
@@ -33,7 +33,7 @@ public class GeometricSamplerTest {
@Test
public void testProbabilityOfSuccessIsOneGeneratesZeroForSamples() {
final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- final GeometricSampler sampler = new GeometricSampler(rng, 1);
+ final SharedStateDiscreteSampler sampler = GeometricSampler.of(rng, 1);
// All samples should be 0
for (int i = 0; i < 10; i++) {
Assert.assertEquals("p=1 should have 0 for all samples", 0,
sampler.sample());
@@ -56,7 +56,7 @@ public class GeometricSamplerTest {
// The internal exponential sampler validates the mean so demonstrate
creating a
// geometric sampler does not throw.
final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- new GeometricSampler(rng, probabilityOfSuccess);
+ GeometricSampler.of(rng, probabilityOfSuccess);
}
/**
@@ -66,7 +66,7 @@ public class GeometricSamplerTest {
@Test
public void testProbabilityOfSuccessIsOneSamplerToString() {
final UniformRandomProvider unusedRng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- final GeometricSampler sampler = new GeometricSampler(unusedRng, 1);
+ final SharedStateDiscreteSampler sampler =
GeometricSampler.of(unusedRng, 1);
Assert.assertTrue("Missing 'Geometric' from toString",
sampler.toString().contains("Geometric"));
}
@@ -83,7 +83,7 @@ public class GeometricSamplerTest {
@Test
public void
testProbabilityOfSuccessIsAlmostZeroGeneratesMaxValueForSamples() {
final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- final GeometricSampler sampler = new GeometricSampler(rng,
Double.MIN_VALUE);
+ final SharedStateDiscreteSampler sampler = GeometricSampler.of(rng,
Double.MIN_VALUE);
// All samples should be max value
for (int i = 0; i < 10; i++) {
Assert.assertEquals("p=(almost 0) should have Integer.MAX_VALUE
for all samples",
@@ -98,8 +98,7 @@ public class GeometricSamplerTest {
public void testProbabilityOfSuccessAboveOneThrows() {
final UniformRandomProvider unusedRng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
final double probabilityOfSuccess = Math.nextUp(1.0);
- @SuppressWarnings("unused")
- final GeometricSampler sampler = new GeometricSampler(unusedRng,
probabilityOfSuccess);
+ GeometricSampler.of(unusedRng, probabilityOfSuccess);
}
/**
@@ -109,8 +108,7 @@ public class GeometricSamplerTest {
public void testProbabilityOfSuccessIsZeroThrows() {
final UniformRandomProvider unusedRng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
final double probabilityOfSuccess = 0;
- @SuppressWarnings("unused")
- final GeometricSampler sampler = new GeometricSampler(unusedRng,
probabilityOfSuccess);
+ GeometricSampler.of(unusedRng, probabilityOfSuccess);
}
/**
@@ -138,8 +136,8 @@ public class GeometricSamplerTest {
private static void testSharedStateSampler(double probabilityOfSuccess) {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
- final GeometricSampler sampler1 =
- new GeometricSampler(rng1, probabilityOfSuccess);
+ final SharedStateDiscreteSampler sampler1 =
+ GeometricSampler.of(rng1, probabilityOfSuccess);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSamplerTest.java
index e38b3db..fc29e0e 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/GuideTableDiscreteSamplerTest.java
@@ -76,7 +76,7 @@ public class GuideTableDiscreteSamplerTest {
@Test
public void testToString() {
- final GuideTableDiscreteSampler sampler = createSampler(new double[]
{0.5, 0.5}, 1.0);
+ final SharedStateDiscreteSampler sampler = createSampler(new double[]
{0.5, 0.5}, 1.0);
Assert.assertTrue(sampler.toString().toLowerCase().contains("guide
table"));
}
@@ -86,9 +86,9 @@ public class GuideTableDiscreteSamplerTest {
* @param probabilities the probabilities
* @return the alias method discrete sampler
*/
- private static GuideTableDiscreteSampler createSampler(double[]
probabilities, double alpha) {
+ private static SharedStateDiscreteSampler createSampler(double[]
probabilities, double alpha) {
final UniformRandomProvider rng =
RandomSource.create(RandomSource.SPLIT_MIX_64);
- return new GuideTableDiscreteSampler(rng, probabilities, alpha);
+ return GuideTableDiscreteSampler.of(rng, probabilities, alpha);
}
/**
@@ -201,7 +201,7 @@ public class GuideTableDiscreteSamplerTest {
* @param alpha the alpha
*/
private static void checkSamples(double[] probabilies, double alpha) {
- final GuideTableDiscreteSampler sampler = createSampler(probabilies,
alpha);
+ final SharedStateDiscreteSampler sampler = createSampler(probabilies,
alpha);
final int numberOfSamples = 10000;
final long[] samples = new long[probabilies.length];
@@ -244,8 +244,8 @@ public class GuideTableDiscreteSamplerTest {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final double[] probabilities = {0.1, 0, 0.2, 0.3, 0.1, 0.3, 0};
- final GuideTableDiscreteSampler sampler1 =
- new GuideTableDiscreteSampler(rng1, probabilities);
+ final SharedStateDiscreteSampler sampler1 =
+ GuideTableDiscreteSampler.of(rng1, probabilities);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSamplerTest.java
index 2796ff0..cc0b7c8 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/KempSmallMeanPoissonSamplerTest.java
@@ -46,7 +46,7 @@ public class KempSmallMeanPoissonSamplerTest {
public void testConstructorThrowsWithMeanLargerThanUpperBound() {
final double mean = SUPPORTED_UPPER_BOUND + 1;
@SuppressWarnings("unused")
- KempSmallMeanPoissonSampler sampler = new
KempSmallMeanPoissonSampler(dummyRng, mean);
+ SharedStateDiscreteSampler sampler =
KempSmallMeanPoissonSampler.of(dummyRng, mean);
}
/**
@@ -56,7 +56,7 @@ public class KempSmallMeanPoissonSamplerTest {
public void testConstructorThrowsWithZeroMean() {
final double mean = 0;
@SuppressWarnings("unused")
- KempSmallMeanPoissonSampler sampler = new
KempSmallMeanPoissonSampler(dummyRng, mean);
+ SharedStateDiscreteSampler sampler =
KempSmallMeanPoissonSampler.of(dummyRng, mean);
}
/**
@@ -66,7 +66,7 @@ public class KempSmallMeanPoissonSamplerTest {
public void testConstructorThrowsWithNegativeMean() {
final double mean = -1;
@SuppressWarnings("unused")
- KempSmallMeanPoissonSampler sampler = new
KempSmallMeanPoissonSampler(dummyRng, mean);
+ SharedStateDiscreteSampler sampler =
KempSmallMeanPoissonSampler.of(dummyRng, mean);
}
/**
@@ -76,7 +76,7 @@ public class KempSmallMeanPoissonSamplerTest {
public void testConstructorWithNaNMean() {
final double mean = Double.NaN;
@SuppressWarnings("unused")
- KempSmallMeanPoissonSampler sampler = new
KempSmallMeanPoissonSampler(dummyRng, mean);
+ SharedStateDiscreteSampler sampler =
KempSmallMeanPoissonSampler.of(dummyRng, mean);
}
/**
@@ -129,7 +129,7 @@ public class KempSmallMeanPoissonSamplerTest {
PoissonDistribution.DEFAULT_EPSILON,
PoissonDistribution.DEFAULT_MAX_ITERATIONS);
final FixedRNG rng = new FixedRNG(0);
- final KempSmallMeanPoissonSampler sampler = new
KempSmallMeanPoissonSampler(rng, mean);
+ final SharedStateDiscreteSampler sampler =
KempSmallMeanPoissonSampler.of(rng, mean);
// Lower bound should be zero
testSample(rng, sampler, 0, 0, 0);
@@ -154,8 +154,8 @@ public class KempSmallMeanPoissonSamplerTest {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final double mean = 1.23;
- final KempSmallMeanPoissonSampler sampler1 =
- new KempSmallMeanPoissonSampler(rng1, mean);
+ final SharedStateDiscreteSampler sampler1 =
+ KempSmallMeanPoissonSampler.of(rng1, mean);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
@@ -169,7 +169,7 @@ public class KempSmallMeanPoissonSamplerTest {
* @param lower the expected lower limit
* @param upper the expected upper limit
*/
- private static void testSample(FixedRNG rng, KempSmallMeanPoissonSampler
sampler, double cumulativeProbability,
+ private static void testSample(FixedRNG rng, SharedStateDiscreteSampler
sampler, double cumulativeProbability,
int lower, int upper) {
rng.setValue(cumulativeProbability);
final int sample = sampler.sample();
diff --git
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSamplerTest.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSamplerTest.java
index 5f2dd87..4b05b1e 100644
---
a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSamplerTest.java
+++
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSamplerTest.java
@@ -85,9 +85,9 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
* @param probabilities Probabilities.
* @return the sampler
*/
- private static MarsagliaTsangWangDiscreteSampler
createDiscreteDistributionSampler(double[] probabilities) {
+ private static SharedStateDiscreteSampler
createDiscreteDistributionSampler(double[] probabilities) {
final UniformRandomProvider rng = new SplitMix64(0L);
- return
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng,
probabilities);
+ return MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng,
probabilities);
}
// Sampling tests
@@ -146,9 +146,9 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final double[] p2 = createProbabilities(offset2, prob);
final double[] p3 = createProbabilities(offset3, prob);
- final MarsagliaTsangWangDiscreteSampler sampler1 =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng1, p1);
- final MarsagliaTsangWangDiscreteSampler sampler2 =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng2, p2);
- final MarsagliaTsangWangDiscreteSampler sampler3 =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng3, p3);
+ final SharedStateDiscreteSampler sampler1 =
MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng1, p1);
+ final SharedStateDiscreteSampler sampler2 =
MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng2, p2);
+ final SharedStateDiscreteSampler sampler3 =
MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng3, p3);
for (int i = 0; i < values.length; i++) {
// Remove offsets
@@ -189,12 +189,12 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
// First test the table is completely filled to 2^30
final UniformRandomProvider dummyRng = new
FixedSequenceIntProvider(new int[] {0xffffffff});
- final MarsagliaTsangWangDiscreteSampler dummySampler =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(dummyRng,
probabilities);
+ final SharedStateDiscreteSampler dummySampler =
MarsagliaTsangWangDiscreteSampler.Enumerated.of(dummyRng, probabilities);
// This will throw if the table is incomplete as it hits the upper
limit
dummySampler.sample();
// Do a test of the actual sampler
- final MarsagliaTsangWangDiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng,
probabilities);
+ final SharedStateDiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng, probabilities);
final int numberOfSamples = 10000;
final long[] samples = new long[probabilities.length];
@@ -300,7 +300,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final double mean = 1025;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, mean);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Poisson.of(rng, mean);
}
/**
@@ -311,7 +311,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final double mean = 0;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, mean);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Poisson.of(rng, mean);
}
/**
@@ -322,7 +322,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final double mean = 1024;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, mean);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Poisson.of(rng, mean);
}
/**
@@ -333,7 +333,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
public void testCreatePoissonDistributionWithSmallMean() {
final UniformRandomProvider rng = new FixedRNG();
final double mean = 0.25;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, mean);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Poisson.of(rng, mean);
// This will throw if the table does not sum to 2^30
sampler.sample();
}
@@ -347,7 +347,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
public void testCreatePoissonDistributionWithMediumMean() {
final UniformRandomProvider rng = new FixedRNG();
final double mean = 21.4;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createPoissonDistribution(rng, mean);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Poisson.of(rng, mean);
// This will throw if the table does not sum to 2^30
sampler.sample();
}
@@ -361,7 +361,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final int trials = -1;
final double p = 0.5;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
}
/**
@@ -373,7 +373,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final int trials = 1 << 16; // 2^16
final double p = 0.5;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
}
/**
@@ -385,7 +385,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final int trials = 1;
final double p = -0.5;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
}
/**
@@ -397,7 +397,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final int trials = 1;
final double p = 1.5;
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
}
/**
@@ -422,7 +422,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
Assert.assertEquals("Invalid test set-up for p(0)", 0,
getBinomialP0(trials + 1, p), 0);
// This will throw if the table does not sum to 2^30
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
sampler.sample();
}
@@ -439,7 +439,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
// Validate set-up
Assert.assertEquals("Invalid test set-up for p(0)", 0,
getBinomialP0(trials, p), 0);
@SuppressWarnings("unused")
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
}
/**
@@ -482,7 +482,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
Assert.assertEquals("Invalid test set-up for p(0)", Double.MIN_VALUE,
getBinomialP0(trials, p), 0);
Assert.assertEquals("Invalid test set-up for p(0)", 0,
getBinomialP0(trials, Math.nextAfter(p, 1)), 0);
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
// This will throw if the table does not sum to 2^30
sampler.sample();
}
@@ -506,7 +506,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final int trials = 1000000;
final double p = 0;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
for (int i = 0; i < 5; i++) {
Assert.assertEquals(0, sampler.sample());
}
@@ -523,7 +523,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final int trials = 1000000;
final double p = 1;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
for (int i = 0; i < 5; i++) {
Assert.assertEquals(trials, sampler.sample());
}
@@ -541,7 +541,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final int trials = 65000;
final double p = 0.01;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
// This will throw if the table does not sum to 2^30
sampler.sample();
}
@@ -555,7 +555,7 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng = new FixedRNG();
final int trials = 10;
final double p = 0.5;
- final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p);
+ final DiscreteSampler sampler =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p);
// This will throw if the table does not sum to 2^30
sampler.sample();
}
@@ -570,8 +570,8 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final int trials = 10;
final double p1 = 0.4;
final double p2 = 1 - p1;
- final DiscreteSampler sampler1 =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p1);
- final DiscreteSampler sampler2 =
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng, trials, p2);
+ final DiscreteSampler sampler1 =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p1);
+ final DiscreteSampler sampler2 =
MarsagliaTsangWangDiscreteSampler.Binomial.of(rng, trials, p2);
Assert.assertEquals(sampler1.toString(), sampler2.toString());
}
@@ -610,8 +610,8 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
double[] probabilities = createProbabilities(offset, prob);
- final MarsagliaTsangWangDiscreteSampler sampler1 =
-
MarsagliaTsangWangDiscreteSampler.createDiscreteDistribution(rng1,
probabilities);
+ final SharedStateDiscreteSampler sampler1 =
+ MarsagliaTsangWangDiscreteSampler.Enumerated.of(rng1,
probabilities);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
@@ -643,8 +643,8 @@ public class MarsagliaTsangWangDiscreteSamplerTest {
private static void testSharedStateSampler(int trials, double
probabilityOfSuccess) {
final UniformRandomProvider rng1 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
final UniformRandomProvider rng2 =
RandomSource.create(RandomSource.SPLIT_MIX_64, 0L);
- final MarsagliaTsangWangDiscreteSampler sampler1 =
-
MarsagliaTsangWangDiscreteSampler.createBinomialDistribution(rng1, trials,
probabilityOfSuccess);
+ final SharedStateDiscreteSampler sampler1 =
+ MarsagliaTsangWangDiscreteSampler.Binomial.of(rng1, trials,
probabilityOfSuccess);
final SharedStateDiscreteSampler sampler2 =
sampler1.withUniformRandomProvider(rng2);
RandomAssert.assertProduceSameSequence(sampler1, sampler2);
}
diff --git a/src/main/resources/pmd/pmd-ruleset.xml
b/src/main/resources/pmd/pmd-ruleset.xml
index 39d64db..40c1312 100644
--- a/src/main/resources/pmd/pmd-ruleset.xml
+++ b/src/main/resources/pmd/pmd-ruleset.xml
@@ -74,7 +74,8 @@
<properties>
<!-- Array is generated internally in this case. -->
<property name="violationSuppressXPath"
- value="//ClassOrInterfaceDeclaration[@Image='PoissonSamplerCache' or
@Image='AliasMethodDiscreteSampler']"/>
+ value="//ClassOrInterfaceDeclaration[@Image='PoissonSamplerCache' or
@Image='AliasMethodDiscreteSampler'
+ or @Image='GuideTableDiscreteSampler']"/>
</properties>
</rule>
