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   * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
29   *
30   * @since 1.1
31   */
32  public class MarsagliaNormalizedGaussianSampler
33      implements NormalizedGaussianSampler, SharedStateContinuousSampler {
34      /** Next gaussian. */
35      private double nextGaussian = Double.NaN;
36      /** Underlying source of randomness. */
37      private final UniformRandomProvider rng;
38  
39      /**
40       * @param rng Generator of uniformly distributed random numbers.
41       */
42      public MarsagliaNormalizedGaussianSampler(UniformRandomProvider rng) {
43          this.rng = rng;
44      }
45  
46      /** {@inheritDoc} */
47      @Override
48      public double sample() {
49          if (Double.isNaN(nextGaussian)) {
50              // Rejection scheme for selecting a pair that lies within the unit circle.
51              while (true) {
52                  // Generate a pair of numbers within [-1 , 1).
53                  final double x = 2 * rng.nextDouble() - 1;
54                  final double y = 2 * rng.nextDouble() - 1;
55                  final double r2 = x * x + y * y;
56  
57                  if (r2 < 1 && r2 > 0) {
58                      // Pair (x, y) is within unit circle.
59                      final double alpha = Math.sqrt(-2 * Math.log(r2) / r2);
60  
61                      // Keep second element of the pair for next invocation.
62                      nextGaussian = alpha * y;
63  
64                      // Return the first element of the generated pair.
65                      return alpha * x;
66                  }
67  
68                  // Pair is not within the unit circle: Generate another one.
69              }
70          }
71  
72          // Use the second element of the pair (generated at the
73          // previous invocation).
74          final double r = nextGaussian;
75  
76          // Both elements of the pair have been used.
77          nextGaussian = Double.NaN;
78  
79          return r;
80      }
81  
82      /** {@inheritDoc} */
83      @Override
84      public String toString() {
85          return "Box-Muller (with rejection) normalized Gaussian deviate [" + rng.toString() + "]";
86      }
87  
88      /**
89       * {@inheritDoc}
90       *
91       * @since 1.3
92       */
93      @Override
94      public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
95          return new MarsagliaNormalizedGaussianSampler(rng);
96      }
97  
98      /**
99       * Create a new normalised Gaussian sampler.
100      *
101      * @param <S> Sampler type.
102      * @param rng Generator of uniformly distributed random numbers.
103      * @return the sampler
104      * @since 1.3
105      */
106     @SuppressWarnings("unchecked")
107     public static <S extends NormalizedGaussianSampler & SharedStateContinuousSampler> S
108             of(UniformRandomProvider rng) {
109         return (S) new MarsagliaNormalizedGaussianSampler(rng);
110     }
111 }