View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.rng.sampling.distribution;
19  
20  import org.apache.commons.rng.UniformRandomProvider;
21  
22  /**
23   * <a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm">
24   * Marsaglia and Tsang "Ziggurat" method</a> for sampling from a Gaussian
25   * distribution with mean 0 and standard deviation 1.
26   *
27   * The algorithm is explained in this
28   * <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf">paper</a>
29   * and this implementation has been adapted from the C code provided therein.
30   *
31   * @since 1.1
32   */
33  public class ZigguratNormalizedGaussianSampler
34      implements NormalizedGaussianSampler {
35      /** Start of tail. */
36      private static final double R = 3.442619855899;
37      /** Inverse of R. */
38      private static final double ONE_OVER_R = 1 / R;
39      /** Rectangle area. */
40      private static final double V = 9.91256303526217e-3;
41      /** 2^63 */
42      private static final double MAX = Math.pow(2, 63);
43      /** 2^-63 */
44      private static final double ONE_OVER_MAX = 1d / MAX;
45      /** Number of entries. */
46      private static final int LEN = 128;
47      /** Index of last entry. */
48      private static final int LAST = LEN - 1;
49      /** Auxiliary table. */
50      private static final long[] K = new long[LEN];
51      /** Auxiliary table. */
52      private static final double[] W = new double[LEN];
53      /** Auxiliary table. */
54      private static final double[] F = new double[LEN];
55      /** Underlying source of randomness. */
56      private final UniformRandomProvider rng;
57  
58      static {
59          // Filling the tables.
60  
61          double d = R;
62          double t = d;
63          double fd = gauss(d);
64          final double q = V / fd;
65  
66          K[0] = (long) ((d / q) * MAX);
67          K[1] = 0;
68  
69          W[0] = q * ONE_OVER_MAX;
70          W[LAST] = d * ONE_OVER_MAX;
71  
72          F[0] = 1;
73          F[LAST] = fd;
74  
75          for (int i = LAST - 1; i >= 1; i--) {
76              d = Math.sqrt(-2 * Math.log(V / d + fd));
77              fd = gauss(d);
78  
79              K[i + 1] = (long) ((d / t) * MAX);
80              t = d;
81  
82              F[i] = fd;
83  
84              W[i] = d * ONE_OVER_MAX;
85          }
86      }
87  
88      /**
89       * @param rng Generator of uniformly distributed random numbers.
90       */
91      public ZigguratNormalizedGaussianSampler(UniformRandomProvider rng) {
92          this.rng = rng;
93      }
94  
95      /** {@inheritDoc} */
96      @Override
97      public double sample() {
98          final long j = rng.nextLong();
99          final int i = (int) (j & LAST);
100         if (Math.abs(j) < K[i]) {
101             return j * W[i];
102         } else {
103             return fix(j, i);
104         }
105     }
106 
107     /** {@inheritDoc} */
108     @Override
109     public String toString() {
110         return "Ziggurat normalized Gaussian deviate [" + rng.toString() + "]";
111     }
112 
113     /**
114      * Gets the value from the tail of the distribution.
115      *
116      * @param hz Start random integer.
117      * @param iz Index of cell corresponding to {@code hz}.
118      * @return the requested random value.
119      */
120     private double fix(long hz,
121                        int iz) {
122         double x;
123         double y;
124 
125         x = hz * W[iz];
126         if (iz == 0) {
127             // Base strip.
128             // This branch is called about 5.7624515E-4 times per sample.
129             do {
130                 y = -Math.log(rng.nextDouble());
131                 x = -Math.log(rng.nextDouble()) * ONE_OVER_R;
132             } while (y + y < x * x);
133 
134             final double out = R + x;
135             return hz > 0 ? out : -out;
136         } else {
137             // Wedge of other strips.
138             // This branch is called about 0.027323 times per sample.
139             if (F[iz] + rng.nextDouble() * (F[iz - 1] - F[iz]) < gauss(x)) {
140                 return x;
141             } else {
142                 // Try again.
143                 // This branch is called about 0.012362 times per sample.
144                 return sample();
145             }
146         }
147     }
148 
149     /**
150      * @param x Argument.
151      * @return \( e^{-\frac{x^2}{2}} \)
152      */
153     private static double gauss(double x) {
154         return Math.exp(-0.5 * x * x);
155     }
156 }