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 * Create an instance.
41 *
42 * @param rng Generator of uniformly distributed random numbers.
43 */
44 public MarsagliaNormalizedGaussianSampler(UniformRandomProvider rng) {
45 this.rng = rng;
46 }
47
48 /** {@inheritDoc} */
49 @Override
50 public double sample() {
51 if (Double.isNaN(nextGaussian)) {
52 // Rejection scheme for selecting a pair that lies within the unit circle.
53 while (true) {
54 // Generate a pair of numbers within [-1 , 1).
55 final double x = 2 * rng.nextDouble() - 1;
56 final double y = 2 * rng.nextDouble() - 1;
57 final double r2 = x * x + y * y;
58
59 if (r2 < 1 && r2 > 0) {
60 // Pair (x, y) is within unit circle.
61 final double alpha = Math.sqrt(-2 * Math.log(r2) / r2);
62
63 // Keep second element of the pair for next invocation.
64 nextGaussian = alpha * y;
65
66 // Return the first element of the generated pair.
67 return alpha * x;
68 }
69
70 // Pair is not within the unit circle: Generate another one.
71 }
72 }
73
74 // Use the second element of the pair (generated at the
75 // previous invocation).
76 final double r = nextGaussian;
77
78 // Both elements of the pair have been used.
79 nextGaussian = Double.NaN;
80
81 return r;
82 }
83
84 /** {@inheritDoc} */
85 @Override
86 public String toString() {
87 return "Box-Muller (with rejection) 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 MarsagliaNormalizedGaussianSampler(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 MarsagliaNormalizedGaussianSampler(rng);
112 }
113 }