RNG-37: fixes
Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/83ab92fc Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/83ab92fc Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/83ab92fc Branch: refs/heads/master Commit: 83ab92fcc765a15cb6de07f5706032aa55bc354e Parents: 60cd84f Author: Olga Kirillova <o...@revunit.com> Authored: Tue Oct 3 15:39:25 2017 -0700 Committer: Olga Kirillova <o...@revunit.com> Committed: Tue Oct 3 15:40:30 2017 -0700 ---------------------------------------------------------------------- .../distribution/ZigguratGaussianSampler.java | 70 +++++++++++++------- 1 file changed, 46 insertions(+), 24 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/83ab92fc/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java index 2673f73..9fae54e 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratGaussianSampler.java @@ -26,6 +26,8 @@ import org.apache.commons.rng.UniformRandomProvider; * The Ziggurat Method for Generating Random Variables * by George Marsaglia and Wai Wan Tsang * + * @see <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf";>Ziggurat Method for Generating Random Variables</a> + * * @since 1.0 */ @@ -35,6 +37,10 @@ public class ZigguratGaussianSampler /** * Generates values from Gaussian (normal) probability distribution + * It uses two tables, integers KN and reals WN. Some 99% of the time, + * the required x is produced by: + * generate a random 32-bit integer j and let i be the index formed from + * the rightmost 8 bits of j. If j < k_i return x = j * w_i. */ private static final int[] KN = new int[128]; @@ -47,46 +53,61 @@ public class ZigguratGaussianSampler static { /** * Filling the tables. + * k_0 = 2^32 * r * f(dn) / vn + * k_i = 2^32 * ( x_{i-1} / x_i ) + * w_0 = .5^32 * vn / f(dn) + * w_i = .5^32 * x_i + * where dn - the rightmost x_i + * vn - the area of the rectangle + * f(dn) = exp(-.5 * dn * dn) */ - final double m = 2147483648.0; + final double m = 2147483648.0; // 2^31 + + /* provides z(r) = 0, where z(r): x_255 = r, vn = r*f(r)+integral_r^inf f(x)dx */ + final double vn = 9.91256303526217e-3; - double dn=3.442619855899; - double tn=dn; - double vn=9.91256303526217e-3; - double q=vn/Math.exp(-.5*dn*dn); + double dn = 3.442619855899; + double tn = dn; + double e = Math.exp(-.5 * dn * dn); + final double q = vn / e; - KN[0]= (int) ((dn/q)*m); - KN[1]= 0; + KN[0] = (int) ((dn / q) * m); + KN[1] = 0; - WN[0]= q/m; - WN[127]= dn/m; + WN[0] = q / m; + WN[127] = dn / m; - FN[0]= 1.0; - FN[127]= Math.exp(-.5*dn*dn); + FN[0] = 1.0; + FN[127] = e; - for (int i=126; i>=1; i--){ - dn=Math.sqrt(-2.*Math.log(vn/dn+Math.exp(-.5*dn*dn))); - KN[i+1]= (int) ((dn/tn)*m); - tn=dn; - FN[i]= Math.exp(-.5*dn*dn); - WN[i]= (dn/m); + for (int i = 126; i >= 1; i--){ + e = Math.exp(-.5 *dn * dn); + dn = Math.sqrt(-2. * Math.log(vn / dn + e)); + KN[i+1] = (int) ((dn / tn) * m); + tn = dn; + FN[i] = e; + WN[i] = dn / m; } } + /** + * @param rng Generator of uniformly distributed random numbers. + */ public ZigguratGaussianSampler(UniformRandomProvider rng) { super(rng); } + /** {@inheritDoc} */ @Override public double sample() { - int j=nextInt(); - int i=j&127; - return (j < KN[i]) ? j*WN[i] : nfix(j,i); + int j = nextInt(); + int i = j & 127; + return (j < KN[i]) ? j * WN[i] : nfix(j,i); } private double nfix(int hz, int iz) { /* The start of the right tail */ - double r = 3.442619855899; + final double r = 3.442619855899; double uni; double x; @@ -99,8 +120,8 @@ public class ZigguratGaussianSampler if (iz == 0) { /* return x from the tail */ do { - x = -Math.log(uni) * 0.2904764; y = -Math.log(uni); + x = y * 0.2904764; uni = .5 + nextInt() * .2328306e-9; } while (y + y < x * x); return (hz > 0) ? r + x : -r - x; @@ -109,14 +130,15 @@ public class ZigguratGaussianSampler if (FN[iz] + uni * (FN[iz - 1] - FN[iz]) < Math.exp(-.5 * x * x)) { return x; } - hz=nextInt(); - iz=hz&127; + hz = nextInt(); + iz = hz & 127; if (Math.abs(hz) < KN[iz]) { return (hz * WN[iz]); } } } + /** {@inheritDoc} */ @Override public String toString() { return "Ziggurat Gaussian deviate [" + super.toString() + "]";
