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 }