You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/08/25 19:54:50 UTC
[commons-rng] branch master updated (b6d5f50 -> 2286528)
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-rng.git.
from b6d5f50 Ensure same samplers are tested in single and sequential
new f420ea0 Add description of sampling method to the ziggurat sampler.
new 2286528 Update observed frequencies with expected where appropriate.
The 2 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../rng/sampling/distribution/ZigguratSampler.java | 463 ++++++++++++++++-----
1 file changed, 351 insertions(+), 112 deletions(-)
[commons-rng] 02/02: Update observed frequencies with expected
where appropriate.
Posted by ah...@apache.org.
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 2286528976920008b523229608e6408101078648
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Aug 25 16:53:44 2021 +0100
Update observed frequencies with expected where appropriate.
The expected frequencies are obtained during the table construction.
They closely match the previously measured observed frequencies.
---
.../rng/sampling/distribution/ZigguratSampler.java | 45 +++++++++++-----------
1 file changed, 22 insertions(+), 23 deletions(-)
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
index ad21c3f..564b274 100644
--- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
@@ -563,15 +563,15 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
if (i < I_MAX) {
// Early exit.
- // This branch is called about 0.984379 times per call into createSample.
- // Note: Frequencies have been empirically measured for the first call to
- // createSample; recursion due to retries have been ignored. Frequencies sum to 1.
+ // Expected frequency = 0.984375
// Drop the sign bit to multiply by [0, 2^63).
return X[i] * (x & MAX_INT64);
}
- // For the first call into createSample:
- // Recursion frequency = 0.000515503
- // Overhang frequency = 0.0151056
+ // Expected frequency = 0.015625
+
+ // Tail frequency = 0.000516062 (recursion)
+ // Overhang frequency = 0.0151089
+
final int j = selectRegion();
// If the tail then exploit the memoryless property of the exponential distribution.
// Add a new sample to the start of the tail (X_0).
@@ -622,17 +622,14 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// _FAST_PRNG_SAMPLE_X(xj, ux)
final double x = sampleX(X, j, u1);
if (uDistance >= E_MAX) {
- // Frequency (per call into createSample): 0.0126230
- // Frequency (per call into expOverhang): 0.823328
// Early Exit: x < y - epsilon
return x;
}
- // Frequency per call into createSample:
- // Return = 0.00248262
- // Recursion = 0.000226013
- // Frequency per call into expOverhang:
- // Return = 0.161930
- // Recursion = 0.0147417
+
+ // Note: Frequencies have been empirically measured per call into expOverhang:
+ // Early Exit = 0.823328
+ // Accept Y = 0.161930
+ // Reject Y = 0.0147417 (recursion)
// _FAST_PRNG_SAMPLE_Y(j, pow(2, 63) - (ux + uDistance))
return sampleY(Y, j, u1 + uDistance) <= Math.exp(-x) ?
@@ -972,7 +969,7 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
if (i < I_MAX) {
// Early exit.
- // Branch frequency: 0.988283
+ // Expected frequency = 0.988281
return X[i] * xx;
}
@@ -989,6 +986,8 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
* @return a sample
*/
private double edgeSample(long xx) {
+ // Expected frequency = 0.0117188
+
// Recycle bits then advance RNG:
// u1 = RANDOM_INT63()
long u1 = xx & MAX_INT64;
@@ -1007,8 +1006,8 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
double x;
if (j > J_INFLECTION) {
// Convex overhang
- // Branch frequency: 0.00892897
- // Loop repeat frequency: 0.389804
+ // Expected frequency: 0.00892899
+ // Observed loop repeat frequency: 0.389804
for (;;) {
x = sampleX(X, j, u1);
// u2 = u1 + (u2 - u1) = u1 + uDistance
@@ -1030,18 +1029,18 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
} else if (j < J_INFLECTION) {
if (j == 0) {
// Tail
- // Branch frequency: 0.000276358
+ // Expected frequency: 0.000276902
// Note: Although less frequent than the next branch, j == 0 is a subset of
// j < J_INFLECTION and must be first.
- // Loop repeat frequency: 0.0634786
+ // Observed loop repeat frequency: 0.0634786
do {
x = ONE_OVER_X_0 * exponential.sample();
} while (exponential.sample() < 0.5 * x * x);
x += X_0;
} else {
// Concave overhang
- // Branch frequency: 0.00249563
- // Loop repeat frequency: 0.0123784
+ // Expected frequency: 0.00249694
+ // Observed loop repeat frequency: 0.0123784
for (;;) {
// u2 = u1 + (u2 - u1) = u1 + uDistance
long uDistance = randomInt63() - u1;
@@ -1061,8 +1060,8 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
}
} else {
// Inflection point
- // Branch frequency: 0.0000159359
- // Loop repeat frequency: 0.500213
+ // Expected frequency: 0.000015914
+ // Observed loop repeat frequency: 0.500213
for (;;) {
x = sampleX(X, j, u1);
if (sampleY(Y, j, randomInt63()) < Math.exp(-0.5 * x * x)) {
[commons-rng] 01/02: Add description of sampling method to the
ziggurat sampler.
Posted by ah...@apache.org.
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 f420ea0d771f118195552cb5cf287e29eba7ed95
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Aug 25 16:36:32 2021 +0100
Add description of sampling method to the ziggurat sampler.
Update the javadoc for all values and methods to be consistent with the
summary description.
---
.../rng/sampling/distribution/ZigguratSampler.java | 418 ++++++++++++++++-----
1 file changed, 329 insertions(+), 89 deletions(-)
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
index 542650e..ad21c3f 100644
--- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratSampler.java
@@ -77,10 +77,12 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// https://github.com/cd-mcfarland/fast_prng
//
// The adaption was faithful to the original as of July-2021.
- // The code uses the naming conventions from the exponential.h and normal.h
- // reference. Comments from the c source have been included.
+ // The code uses similar naming conventions from the exponential.h and normal.h
+ // reference. Naming has been updated to be consistent in the exponential and normal
+ // samplers. Comments from the c source have been included.
// When the c reference was updated to use an underlying RNG
- // available in Commons RNG the c and java code output the same random
+ // available in Commons RNG outputting the same signed and positive long values
+ // the c and java code output the same random
// deviates over 2^30 cycles if identically seeded. Branch frequencies have
// been measured and added as comments.
//
@@ -113,20 +115,161 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// probability of 2^-49 (1.78e-15), or 2^-45 (2.84e-14) respectively. The tables may be
// regenerated in future versions if the reference script receives updates to improve
// accuracy.
+ //
+ // Method Description
+ //
+ // The ziggurat is constructed using layers that fit exactly within the probability density
+ // function. Each layer has the same area. This area is chosen to be a fraction of the total
+ // area under the PDF with the denominator of the fraction a power of 2. These tables
+ // use 1/256 as the volume of each layer. The remaining part of the PDF that is not represented
+ // by the layers is the overhang. There is an overhang above each layer and a final tail.
+ // The following is a ziggurat with 3 layers:
+ //
+ // Y3 |\
+ // | \ j=3
+ // | \
+ // Y2 | \
+ // |----\
+ // | |\
+ // | i=2| \ j=2
+ // | | \
+ // Y1 |--------\
+ // | i=1 | \ j=1
+ // | | \
+ // Y0 |-----------\
+ // | i=0 | \ j=0 (tail)
+ // +--------------
+ // X3 | | X0
+ // | X1
+ // X2
+ //
+ // There are N layers referenced using i in [0, N). The overhangs are referenced using
+ // j in [1, N]; j=0 is the tail. Note that N is < 256.
+ // Information about the ziggurat is pre-computed:
+ // X = The length of each layer (supplemented with zero for Xn)
+ // Y = PDF(X) for each layer (supplemented with PDF(x=0) for Yn)
+ //
+ // Sampling is performed as:
+ // - Pick index i in [0, 256).
+ // - If i is a layer then return a uniform deviate multiplied by the layer length
+ // - If i is not a layer then sample from the overhang or tail
+ //
+ // The overhangs and tail have different volumes. Sampling must pick a region j based the
+ // probability p(j) = vol(j) / sum (vol(j)). This is performed using alias sampling.
+ // (See Walker, AJ (1977) "An Efficient Method for Generating Discrete Random Variables with
+ // General Distributions" ACM Transactions on Mathematical Software 3 (3), 253-256.)
+ // This uses a table that has been constructed to evenly balance A categories with
+ // probabilities around the mean into B sections each allocated the 'mean'. For the 4
+ // regions in the ziggurat shown above balanced into 8 sections:
+ //
+ // 3
+ // 3
+ // 32
+ // 32
+ // 321
+ // 321 => 31133322
+ // 3210 01233322
+ //
+ // section abcdefgh
+ //
+ // A section with an index below the number of categories represents the category j and
+ // optionally an alias. Sections with an index above the number
+ // of categories are entirely filled with the alias. The region is chosen
+ // by selecting a section and then checking if a uniform deviate is above the alias
+ // threshold. If so then the alias is used in place of the original index.
+ //
+ // Alias sampling uses a table size of 256. This allows fast computation of the index
+ // as a power of 2. The probability threshold is stored as the numerator of a fraction
+ // allowing direct comparison with a uniform long deviate.
+ //
+ // MAP = Alias map for j in [0, 256)
+ // IPMF = Alias probability threshold for j
+ //
+ // Note: The IPMF table is larger than the number of regions. Thus the final entries
+ // must represent a probability of zero so that the alias is always used.
+ //
+ // If the selected region j is the tail then sampling uses a sampling method appropriate
+ // for the PDF. If the selected region is an overhang then sampling generates a random
+ // coordinate inside the rectangle covering the overhang using random deviates u1 and u2:
+ //
+ // X[j],Y[j]
+ // |\-->u1
+ // | \ |
+ // | \ |
+ // | \| Overhang j (with hypotenuse not pdf(x))
+ // | \
+ // | |\
+ // | | \
+ // | u2 \
+ // +-------- X[j-1],Y[j-1]
+ //
+ // The random point (x,y) has coordinates:
+ // x = X[j] + u1 * (X[j-1] - X[j])
+ // y = Y[j] + u2 * (Y[j-1] - Y[j])
+ //
+ // The expressions can be evaluated from the opposite direction using (1-u), e.g:
+ // y = Y[j-1] + (1-u2) * (Y[j] - Y[j-1])
+ // This allows the large value to subtract the small value before multiplying by u.
+ // Note that the tables X and Y have been scaled by 2^-63. This allows U to be a uniform
+ // long in [0, 2^63). Thus the value c in 'c + m * x' must be scaled up by 2^63.
+ //
+ // If point (x,y) is below pdf(x) then the sample is accepted.
+ // If u2 > u1 then the point is below the hypotenuse.
+ // If u1 > u2 then the point is above the hypotenuse.
+ // The distance above/below the hypotenuse is the difference u2 - u1: negative is above;
+ // positive is below.
+ //
+ // The pdf(x) may lie completely above or below the hypotenuse. If the region under the pdf
+ // is inside then this is referred to as convex (above) and concave (below). The
+ // exponential function is concave for all regions. The normal function is convex below
+ // x=1, and concave above x=1. x=1 is the point of inflection.
+ //
+ // Concave Convex
+ // |- |----
+ // | - | ---
+ // | - | --
+ // | -- | --
+ // | -- | -
+ // | --- | -
+ // | ---- | -
+ //
+ // Optimisations:
+ //
+ // Regions that are concave can detect a point (x,y) above the hypotenuse and reflect the
+ // point in the hypotenuse by swapping u1 and u2.
+ // Regions that are convex can detect a point (x,y) below the hypotenuse and immediate accept
+ // the sample.
+ // The maximum distance of pdf(x) from the hypotenuse can be precomputed. This can be done for
+ // each region or by taking the largest distance across all regions. This value can be
+ // compared to the distance between u1 and u2 and the point immediately accepted (concave)
+ // or rejected (convex) as it is known to be respectively inside or outside the pdf.
+ // This sampler uses a single value for the maximum distance of pdf(x) from the hypotenuse.
+ // For the normal distribution this is two values to separate the maximum for convex and
+ // concave regions.
+ //
// =========================================================================
/**
* Modified ziggurat method for sampling from an exponential distributions.
*/
public static class Exponential extends ZigguratSampler {
- /** Maximum i value for early exit. */
+ // Ziggurat volumes:
+ // Inside the layers = 98.4375% (252/256)
+ // Fraction outside the layers:
+ // concave overhangs = 96.6972%
+ // tail = 3.3028%
+
+ /** The number of layers in the ziggurat. Maximum i value for early exit. */
private static final int I_MAX = 252;
- /** Maximum distance value for early exit. Equal to approximately 0.0926 scaled by 2^63. */
- private static final long IE_MAX = 853965788476313646L;
- /** Beginning of tail. */
+ /** Maximum deviation of concave pdf(x) below the hypotenuse value for early exit.
+ * Equal to approximately 0.0926 scaled by 2^63. */
+ private static final long E_MAX = 853965788476313646L;
+ /** Beginning of tail. Equal to X[0] * 2^63. */
private static final double X_0 = 7.569274694148063;
- /** The alias map. An integer in [0, 255] stored as a byte to save space. */
+ /** The alias map. An integer in [0, 255] stored as a byte to save space.
+ * Contains the alias j for each index. j=0 is the tail; j in [1, N] is the overhang
+ * for each layer. */
private static final byte[] MAP = toBytes(new int[] {0, 0, 1, 235, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
1, 1, 1, 2, 2, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252,
252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252,
@@ -139,7 +282,10 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
251, 250, 250, 250, 250, 250, 250, 250, 249, 249, 249, 249, 249, 249, 248, 248, 248, 248, 247, 247, 247,
247, 246, 246, 246, 245, 245, 244, 244, 243, 243, 242, 241, 241, 240, 239, 237, 3, 3, 4, 4, 6, 0, 0, 0, 0,
236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 2, 0, 0, 0});
- /** The alias inverse PMF. */
+ /** The alias inverse PMF. This is the probability threshold to use the alias for j in-place of j.
+ * This has been scaled by 2^64 and offset by -2^63. It represents the numerator of a fraction
+ * with denominator 2^64 and can be compared directly to a uniform long deviate.
+ * The value 0 is Long.MIN_VALUE and is used when {@code j > I_MAX}. */
private static final long[] IPMF = {9223372036854774016L, 1623796909450834944L, 2664290944894291200L,
7387971354164060928L, 6515064486552723200L, 8840508362680718848L, 6099647593382936320L,
7673130333659513856L, 6220332867583438080L, 5045979640552813824L, 4075305837223955456L,
@@ -205,8 +351,14 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
3843116882594676224L, 3927231442413903104L, -9223372036854775808L, -9223372036854775808L,
-9223372036854775808L};
/**
- * The precomputed ziggurat lengths, denoted X_i in the main text. X_i = length of
- * ziggurat layer i.
+ * The precomputed ziggurat lengths, denoted X_i in the main text.
+ * <ul>
+ * <li>X_i = length of ziggurat layer i.
+ * <li>X_j is the upper-left X coordinate of overhang j (starting from 1).
+ * <li>X_(j-1) is the lower-right X coordinate of overhang j.
+ * </ul>
+ * <p>Values have been scaled by 2^-63.
+ * Contains {@code I_MAX + 1} entries as the final value is 0.
*/
private static final double[] X = {8.206624067534882E-19, 7.397373235160728E-19, 6.913331337791529E-19,
6.564735882096453E-19, 6.291253995981851E-19, 6.065722412960496E-19, 5.873527610373727E-19,
@@ -272,7 +424,16 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
3.4067836691100565E-20, 3.2128447641564046E-20, 3.0095646916399994E-20, 2.794846945559833E-20,
2.5656913048718645E-20, 2.317520975680391E-20, 2.042669522825129E-20, 1.7261770330213488E-20,
1.3281889259442579E-20, 0.0};
- /** Overhang table. Y_i = f(X_i). */
+ /**
+ * The precomputed ziggurat heights, denoted Y_i in the main text.
+ * <ul>
+ * <li>Y_i = height of ziggurat layer i.
+ * <li>Y_j is the upper-left Y coordinate of overhang j (starting from 1).
+ * <li>Y_(j-1) is the lower-right Y coordinate of overhang j.
+ * </ul>
+ * <p>Values have been scaled by 2^-63.
+ * Contains {@code I_MAX + 1} entries as the final value is pdf(x=0).
+ */
private static final double[] Y = {5.595205495112736E-23, 1.1802509982703313E-22, 1.844442338673583E-22,
2.543903046669831E-22, 3.2737694311509334E-22, 4.0307732132706715E-22, 4.812547831949511E-22,
5.617291489658331E-22, 6.443582054044353E-22, 7.290266234346368E-22, 8.156388845632194E-22,
@@ -405,48 +566,62 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// This branch is called about 0.984379 times per call into createSample.
// Note: Frequencies have been empirically measured for the first call to
// createSample; recursion due to retries have been ignored. Frequencies sum to 1.
+ // Drop the sign bit to multiply by [0, 2^63).
return X[i] * (x & MAX_INT64);
}
// For the first call into createSample:
// Recursion frequency = 0.000515503
// Overhang frequency = 0.0151056
- final int j = expSampleA();
- return j == 0 ? X_0 + createSample() : expOverhang(j);
+ final int j = selectRegion();
+ // If the tail then exploit the memoryless property of the exponential distribution.
+ // Add a new sample to the start of the tail (X_0).
+ return j == 0 ? X_0 + createSample() : sampleOverhang(j);
}
/**
- * Alias sampling.
- * See http://scorevoting.net/WarrenSmithPages/homepage/sampling.abs
+ * Select the overhang region or the tail using alias sampling.
*
- * @return the alias
+ * @return the region
*/
- private int expSampleA() {
+ private int selectRegion() {
final long x = nextLong();
- // j <- I(0, 256)
+ // j in [0, 256)
final int j = ((int) x) & MASK_INT8;
+ // map to j in [0, N] with N the number of layers of the ziggurat
return x >= IPMF[j] ? MAP[j] & MASK_INT8 : j;
}
/**
- * Draws a PRN from overhang.
+ * Sample from overhang region {@code j}.
*
* @param j Index j (must be {@code > 0})
* @return the sample
*/
- private double expOverhang(int j) {
- // To sample a unit right-triangle:
- // U_x <- min(U_1, U_2)
- // distance <- | U_1 - U_2 |
- // U_y <- 1 - (U_x + distance)
- long ux = randomInt63();
- long uDistance = randomInt63() - ux;
+ private double sampleOverhang(int j) {
+ // Sample from the triangle:
+ // X[j],Y[j]
+ // |\-->u1
+ // | \ |
+ // | \ |
+ // | \| Overhang j (with hypotenuse not pdf(x))
+ // | \
+ // | |\
+ // | | \
+ // | u2 \
+ // +-------- X[j-1],Y[j-1]
+ // u2 = u1 + (u2 - u1) = u1 + uDistance
+ // If u2 < u1 then reflect in the hypotenuse by swapping u1 and u2.
+ long u1 = randomInt63();
+ long uDistance = randomInt63() - u1;
if (uDistance < 0) {
+ // Upper-right triangle. Reflect in hypotenuse.
uDistance = -uDistance;
- ux -= uDistance;
+ // Update u1 to be min(u1, u2) by subtracting the distance between them
+ u1 -= uDistance;
}
// _FAST_PRNG_SAMPLE_X(xj, ux)
- final double x = fastPrngSampleX(X, j, ux);
- if (uDistance >= IE_MAX) {
+ final double x = sampleX(X, j, u1);
+ if (uDistance >= E_MAX) {
// Frequency (per call into createSample): 0.0126230
// Frequency (per call into expOverhang): 0.823328
// Early Exit: x < y - epsilon
@@ -460,10 +635,8 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// Recursion = 0.0147417
// _FAST_PRNG_SAMPLE_Y(j, pow(2, 63) - (ux + uDistance))
- // Long.MIN_VALUE is used as an unsigned int with value 2^63:
- // uy = Long.MIN_VALUE - (ux + uDistance)
- return fastPrngSampleY(Y, j, Long.MIN_VALUE - (ux + uDistance)) <= Math.exp(-x) ?
- x : expOverhang(j);
+ return sampleY(Y, j, u1 + uDistance) <= Math.exp(-x) ?
+ x : sampleOverhang(j);
}
/** {@inheritDoc} */
@@ -512,21 +685,35 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
*/
public static final class NormalizedGaussian extends ZigguratSampler
implements NormalizedGaussianSampler, SharedStateContinuousSampler {
- /** Maximum i value for early exit. */
+ // Ziggurat volumes:
+ // Inside the layers = 98.8281% (253/256)
+ // Fraction outside the layers:
+ // concave overhangs = 76.1941%
+ // inflection overhang = 0.1358%
+ // convex overhangs = 21.3072%
+ // tail = 2.3629%
+
+ /** The number of layers in the ziggurat. Maximum i value for early exit. */
private static final int I_MAX = 253;
/** The point where the Gaussian switches from convex to concave.
* This is the largest value of X[j] below 1. */
private static final int J_INFLECTION = 204;
- /** Used for largest deviations of f(x) from y_i. This is negated on purpose. */
- private static final long MAX_IE = -2269182951627976004L;
- /** Used for largest deviations of f(x) from y_i. */
- private static final long MIN_IE = 760463704284035184L;
- /** Beginning of tail. */
+ /** Maximum epsilon distance of convex pdf(x) above the hypotenuse value for early rejection.
+ * Equal to approximately 0.2460 scaled by 2^63. This is negated on purpose as the
+ * distance for a point (x,y) above the hypotenuse is negative:
+ * {@code (|d| < max) == (d >= -max)}. */
+ private static final long CONVEX_E_MAX = -2269182951627976004L;
+ /** Maximum distance of concave pdf(x) below the hypotenuse value for early exit.
+ * Equal to approximately 0.08244 scaled by 2^63. */
+ private static final long CONCAVE_E_MAX = 760463704284035184L;
+ /** Beginning of tail. Equal to X[0] * 2^63. */
private static final double X_0 = 3.6360066255009455861;
- /** 1/X_0. */
- private static final double ONE_OVER_X_0 = 1d / X_0;
+ /** 1/X_0. Used for tail sampling. */
+ private static final double ONE_OVER_X_0 = 1.0 / X_0;
- /** The alias map. An integer in [0, 255] stored as a byte to save space. */
+ /** The alias map. An integer in [0, 255] stored as a byte to save space.
+ * Contains the alias j for each index. j=0 is the tail; j in [1, N] is the overhang
+ * for each layer. */
private static final byte[] MAP = toBytes(
new int[] {0, 0, 239, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253,
253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253,
@@ -540,7 +727,10 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
253, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 251, 251, 251, 251, 251, 251, 251, 250, 250, 250,
250, 250, 249, 249, 249, 248, 248, 248, 247, 247, 247, 246, 246, 245, 244, 244, 243, 242, 240, 2, 2, 3,
3, 0, 0, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 1, 0, 0});
- /** The alias inverse PMF. */
+ /** The alias inverse PMF. This is the probability threshold to use the alias for j in-place of j.
+ * This has been scaled by 2^64 and offset by -2^63. It represents the numerator of a fraction
+ * with denominator 2^64 and can be compared directly to a uniform long deviate.
+ * The value 0 is Long.MIN_VALUE and is used when {@code j > I_MAX}. */
private static final long[] IPMF = {9223372036854775296L, 1100243796534090752L, 7866600928998383104L,
6788754710675124736L, 9022865200181688320L, 6522434035205502208L, 4723064097360024576L,
3360495653216416000L, 2289663232373870848L, 1423968905551920384L, 708364817827798016L, 106102487305601280L,
@@ -606,8 +796,14 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
2889316805630113280L, -2712587580533804032L, 6562498853538167040L, 7975754821147501312L,
-9223372036854775808L, -9223372036854775808L};
/**
- * The precomputed ziggurat lengths, denoted X_i in the main text. X_i = length of
- * ziggurat layer i.
+ * The precomputed ziggurat lengths, denoted X_i in the main text.
+ * <ul>
+ * <li>X_i = length of ziggurat layer i.
+ * <li>X_j is the upper-left X coordinate of overhang j (starting from 1).
+ * <li>X_(j-1) is the lower-right X coordinate of overhang j.
+ * </ul>
+ * <p>Values have been scaled by 2^-63.
+ * Contains {@code I_MAX + 1} entries as the final value is 0.
*/
private static final double[] X = {3.9421662825398133E-19, 3.720494500411901E-19, 3.582702448062868E-19,
3.480747623654025E-19, 3.3990177171882136E-19, 3.330377836034014E-19, 3.270943881761755E-19,
@@ -673,7 +869,16 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
6.092168754828126E-20, 5.885987357557682E-20, 5.666267511609098E-20, 5.430181363089457E-20,
5.173817174449422E-20, 4.8915031722398545E-20, 4.57447418907553E-20, 4.2078802568583416E-20,
3.762598672240476E-20, 3.162858980588188E-20, 0.0};
- /** Overhang table. Y_i = f(X_i). */
+ /**
+ * The precomputed ziggurat heights, denoted Y_i in the main text.
+ * <ul>
+ * <li>Y_i = height of ziggurat layer i.
+ * <li>Y_j is the upper-left Y coordinate of overhang j (starting from 1).
+ * <li>Y_(j-1) is the lower-right Y coordinate of overhang j.
+ * </ul>
+ * <p>Values have been scaled by 2^-63.
+ * Contains {@code I_MAX + 1} entries as the final value is pdf(x=0).
+ */
private static final double[] Y = {1.4598410796619063E-22, 3.0066613427942797E-22, 4.612972881510347E-22,
6.266335004923436E-22, 7.959452476188154E-22, 9.687465502170504E-22, 1.144687700237944E-21,
1.3235036304379167E-21, 1.504985769205313E-21, 1.6889653000719298E-21, 1.8753025382711626E-21,
@@ -791,7 +996,7 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// double sign_bit = u1 & 0x100 ? 1. : -1.
// Use 2 - 1 or 0 - 1
final double signBit = ((u1 >>> 7) & 0x2) - 1.0;
- final int j = normSampleA();
+ final int j = selectRegion();
// Four kinds of overhangs:
// j = 0 : Sample from tail
// 0 < j < J_INFLECTION : Overhang is concave; only sample from Lower-Left triangle
@@ -805,22 +1010,21 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// Branch frequency: 0.00892897
// Loop repeat frequency: 0.389804
for (;;) {
- x = fastPrngSampleX(X, j, u1);
- final long uDiff = randomInt63() - u1;
- if (uDiff >= 0) {
+ x = sampleX(X, j, u1);
+ // u2 = u1 + (u2 - u1) = u1 + uDistance
+ final long uDistance = randomInt63() - u1;
+ if (uDistance >= 0) {
// Lower-left triangle
break;
}
- if (uDiff >= MAX_IE &&
+ if (uDistance >= CONVEX_E_MAX &&
// Within maximum distance of f(x) from the triangle hypotenuse.
// Frequency (per upper-right triangle): 0.431497
// Reject frequency: 0.489630
- // Long.MIN_VALUE is used as an unsigned int with value 2^63:
- // uy = Long.MIN_VALUE - (ux + uDiff)
- fastPrngSampleY(Y, j, Long.MIN_VALUE - (u1 + uDiff)) < Math.exp(-0.5 * x * x)) {
+ sampleY(Y, j, u1 + uDistance) < Math.exp(-0.5 * x * x)) {
break;
}
- // uDiff < MAX_IE (upper-right triangle) or rejected as above the curve
+ // uDistance < MAX_IE (upper-right triangle) or rejected as above the curve
u1 = randomInt63();
}
} else if (j < J_INFLECTION) {
@@ -839,18 +1043,17 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// Branch frequency: 0.00249563
// Loop repeat frequency: 0.0123784
for (;;) {
- // U_x <- min(U_1, U_2)
- // distance <- | U_1 - U_2 |
- // U_y <- 1 - (U_x + distance)
- long uDiff = randomInt63() - u1;
- if (uDiff < 0) {
+ // u2 = u1 + (u2 - u1) = u1 + uDistance
+ long uDistance = randomInt63() - u1;
+ if (uDistance < 0) {
// Upper-right triangle. Reflect in hypotenuse.
- uDiff = -uDiff;
- u1 -= uDiff;
+ uDistance = -uDistance;
+ // Update u1 to be min(u1, u2) by subtracting the distance between them
+ u1 -= uDistance;
}
- x = fastPrngSampleX(X, j, u1);
- if (uDiff > MIN_IE ||
- fastPrngSampleY(Y, j, Long.MIN_VALUE - (u1 + uDiff)) < Math.exp(-0.5 * x * x)) {
+ x = sampleX(X, j, u1);
+ if (uDistance > CONCAVE_E_MAX ||
+ sampleY(Y, j, u1 + uDistance) < Math.exp(-0.5 * x * x)) {
break;
}
u1 = randomInt63();
@@ -861,8 +1064,8 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
// Branch frequency: 0.0000159359
// Loop repeat frequency: 0.500213
for (;;) {
- x = fastPrngSampleX(X, j, u1);
- if (fastPrngSampleY(Y, j, randomInt63()) < Math.exp(-0.5 * x * x)) {
+ x = sampleX(X, j, u1);
+ if (sampleY(Y, j, randomInt63()) < Math.exp(-0.5 * x * x)) {
break;
}
u1 = randomInt63();
@@ -872,15 +1075,15 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
}
/**
- * Alias sampling.
- * See http://scorevoting.net/WarrenSmithPages/homepage/sampling.abs
+ * Select the overhang region or the tail using alias sampling.
*
- * @return the alias
+ * @return the region
*/
- private int normSampleA() {
+ private int selectRegion() {
final long x = nextLong();
- // j <- I(0, 256)
+ // j in [0, 256)
final int j = ((int) x) & MASK_INT8;
+ // map to j in [0, N] with N the number of layers of the ziggurat
return x >= IPMF[j] ? MAP[j] & MASK_INT8 : j;
}
@@ -942,30 +1145,67 @@ public abstract class ZigguratSampler implements SharedStateContinuousSampler {
}
/**
- * Auxiliary function to see if rejection sampling is required in the overhang. See
- * Fig. 2 in the main text.
+ * Compute the x value of the point in the overhang region from the uniform deviate.
+ * <pre>{@code
+ * X[j],Y[j]
+ * |\
+ * | \
+ * | \
+ * | \ Overhang j (with hypotenuse not pdf(x))
+ * | \
+ * | \
+ * | \
+ * |-->u1 \
+ * +-------- X[j-1],Y[j-1]
*
- * @param x Precomputed ziggurat lengths, denoted X_i in the main text.
- * X_i = length of ziggurat layer i.
+ * x = X[j] + u1 * (X[j-1] - X[j])
+ * }</pre>
+ * <p>See Fig. 2 in the main text.
+ *
+ * @param x Precomputed ziggurat lengths, X_j = length of ziggurat layer j.
* @param j j
- * @param ux ux
- * @return the sample
+ * @param u1 u1
+ * @return x
*/
- static double fastPrngSampleX(double[] x, int j, long ux) {
- return x[j] * TWO_POW_63 + (x[j - 1] - x[j]) * ux;
+ static double sampleX(double[] x, int j, long u1) {
+ return x[j] * TWO_POW_63 + u1 * (x[j - 1] - x[j]);
}
/**
- * Auxiliary function to see if rejection sampling is required in the overhang. See
- * Fig. 2 in the main text.
+ * Compute the y value of the point in the overhang region from the uniform deviate.
+ * <pre>{@code
+ * X[j],Y[j]
+ * |\ |
+ * | \|
+ * | \
+ * | |\ Overhang j (with hypotenuse not pdf(x))
+ * | | \
+ * | u2 \
+ * | \
+ * | \
+ * +-------- X[j-1],Y[j-1]
+ *
+ * y = Y[j-1] + (1-u2) * (Y[j] - Y[j-1])
+ * }</pre>
+ * <p>See Fig. 2 in the main text.
*
- * @param y Overhang table. Y_i = f(X_i).
- * @param i i
- * @param uy uy
- * @return the sample
+ * @param y Overhang table. Y_j = f(X_j).
+ * @param j j
+ * @param u2 u2
+ * @return y
*/
- static double fastPrngSampleY(double[] y, int i, long uy) {
- return y[i - 1] * TWO_POW_63 + (y[i] - y[i - 1]) * uy;
+ static double sampleY(double[] y, int j, long u2) {
+ // Note: u2 is in [0, 2^63)
+ // Long.MIN_VALUE is used as an unsigned int with value 2^63:
+ // 1 - u2 = Long.MIN_VALUE - u2
+ //
+ // The subtraction from 1 can be avoided with:
+ // y = Y[j] + u2 * (Y[j-1] - Y[j])
+ // This is effectively sampleX(y, j, u2)
+ // Tests show the alternative is 1 ULP different with approximately 3% frequency.
+ // It is never more than 1 ULP different.
+
+ return y[j - 1] * TWO_POW_63 + (Long.MIN_VALUE - u2) * (y[j] - y[j - 1]);
}
/**