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/Box%E2%80%93Muller_transform">
23   * Box-Muller algorithm</a> for sampling from Gaussian distribution with
24   * mean 0 and standard deviation 1.
25   *
26   * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
27   *
28   * @since 1.1
29   */
30  public class BoxMullerNormalizedGaussianSampler
31      implements NormalizedGaussianSampler, SharedStateContinuousSampler {
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 BoxMullerNormalizedGaussianSampler(UniformRandomProvider rng) {
41          this.rng = rng;
42      }
43  
44      /** {@inheritDoc} */
45      @Override
46      public double sample() {
47          double random;
48          if (Double.isNaN(nextGaussian)) {
49              // Generate a pair of Gaussian numbers.
50  
51              // Avoid zero for the uniform deviate y.
52              // The extreme tail of the sample is:
53              // y = 2^-53
54              // r = 8.57167
55              final double x = rng.nextDouble();
56              final double y = InternalUtils.makeNonZeroDouble(rng.nextLong());
57              final double alpha = 2 * Math.PI * x;
58              final double r = Math.sqrt(-2 * Math.log(y));
59  
60              // Return the first element of the generated pair.
61              random = r * Math.cos(alpha);
62  
63              // Keep second element of the pair for next invocation.
64              nextGaussian = r * Math.sin(alpha);
65          } else {
66              // Use the second element of the pair (generated at the
67              // previous invocation).
68              random = nextGaussian;
69  
70              // Both elements of the pair have been used.
71              nextGaussian = Double.NaN;
72          }
73  
74          return random;
75      }
76  
77      /** {@inheritDoc} */
78      @Override
79      public String toString() {
80          return "Box-Muller normalized Gaussian deviate [" + rng.toString() + "]";
81      }
82  
83      /**
84       * {@inheritDoc}
85       *
86       * @since 1.3
87       */
88      @Override
89      public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
90          return new BoxMullerNormalizedGaussianSampler(rng);
91      }
92  
93      /**
94       * Create a new normalised Gaussian sampler.
95       *
96       * @param <S> Sampler type.
97       * @param rng Generator of uniformly distributed random numbers.
98       * @return the sampler
99       * @since 1.3
100      */
101     @SuppressWarnings("unchecked")
102     public static <S extends NormalizedGaussianSampler & SharedStateContinuousSampler> S
103             of(UniformRandomProvider rng) {
104         return (S) new BoxMullerNormalizedGaussianSampler(rng);
105     }
106 }