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 }