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