commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [02/10] commons-rng git commit: RNG-37: fixes
Date Wed, 11 Oct 2017 08:45:29 GMT
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 <olga@revunit.com>
Authored: Tue Oct 3 15:39:25 2017 -0700
Committer: Olga Kirillova <olga@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() + "]";


Mime
View raw message