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  package org.apache.commons.rng.sampling.distribution;
18  
19  import org.apache.commons.rng.UniformRandomProvider;
20  
21  /**
22   * <a href="https://en.wikipedia.org/wiki/Marsaglia_polar_method">
23   * Marsaglia polar method</a> for sampling from a Gaussian distribution
24   * with mean 0 and standard deviation 1.
25   * This is a variation of the algorithm implemented in
26   * {@link BoxMullerNormalizedGaussianSampler}.
27   *
28   * @since 1.1
29   */
30  public class MarsagliaNormalizedGaussianSampler
31      implements NormalizedGaussianSampler {
32      /** Next gaussian. */
33      private double nextGaussian = Double.NaN;
34      /** Underlying source of randomness. */
35      private final UniformRandomProvider rng;
36  
37      /**
38       * @param rng Generator of uniformly distributed random numbers.
39       */
40      public MarsagliaNormalizedGaussianSampler(UniformRandomProvider rng) {
41          this.rng = rng;
42      }
43  
44      /** {@inheritDoc} */
45      @Override
46      public double sample() {
47          if (Double.isNaN(nextGaussian)) {
48              // Rejection scheme for selecting a pair that lies within the unit circle.
49              while (true) {
50                  // Generate a pair of numbers within [-1 , 1).
51                  final double x = 2 * rng.nextDouble() - 1;
52                  final double y = 2 * rng.nextDouble() - 1;
53                  final double r2 = x * x + y * y;
54  
55                  if (r2 < 1 && r2 > 0) {
56                      // Pair (x, y) is within unit circle.
57                      final double alpha = Math.sqrt(-2 * Math.log(r2) / r2);
58  
59                      // Keep second element of the pair for next invocation.
60                      nextGaussian = alpha * y;
61  
62                      // Return the first element of the generated pair.
63                      return alpha * x;
64                  }
65  
66                  // Pair is not within the unit circle: Generate another one.
67              }
68          } else {
69              // Use the second element of the pair (generated at the
70              // previous invocation).
71              final double r = nextGaussian;
72  
73              // Both elements of the pair have been used.
74              nextGaussian = Double.NaN;
75  
76              return r;
77          }
78      }
79  
80      /** {@inheritDoc} */
81      @Override
82      public String toString() {
83          return "Box-Muller (with rejection) normalized Gaussian deviate [" + rng.toString() + "]";
84      }
85  }