001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.rng.sampling.distribution; 018 019import org.apache.commons.rng.UniformRandomProvider; 020 021/** 022 * <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform"> 023 * Box-Muller algorithm</a> for sampling from Gaussian distribution with 024 * mean 0 and standard deviation 1. 025 * 026 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p> 027 * 028 * @since 1.1 029 */ 030public class BoxMullerNormalizedGaussianSampler 031 implements NormalizedGaussianSampler, SharedStateContinuousSampler { 032 /** Next gaussian. */ 033 private double nextGaussian = Double.NaN; 034 /** Underlying source of randomness. */ 035 private final UniformRandomProvider rng; 036 037 /** 038 * @param rng Generator of uniformly distributed random numbers. 039 */ 040 public BoxMullerNormalizedGaussianSampler(UniformRandomProvider rng) { 041 this.rng = rng; 042 } 043 044 /** {@inheritDoc} */ 045 @Override 046 public double sample() { 047 double random; 048 if (Double.isNaN(nextGaussian)) { 049 // Generate a pair of Gaussian numbers. 050 051 // Avoid zero for the uniform deviate y. 052 // The extreme tail of the sample is: 053 // y = 2^-53 054 // r = 8.57167 055 final double x = rng.nextDouble(); 056 final double y = InternalUtils.makeNonZeroDouble(rng.nextLong()); 057 final double alpha = 2 * Math.PI * x; 058 final double r = Math.sqrt(-2 * Math.log(y)); 059 060 // Return the first element of the generated pair. 061 random = r * Math.cos(alpha); 062 063 // Keep second element of the pair for next invocation. 064 nextGaussian = r * Math.sin(alpha); 065 } else { 066 // Use the second element of the pair (generated at the 067 // previous invocation). 068 random = nextGaussian; 069 070 // Both elements of the pair have been used. 071 nextGaussian = Double.NaN; 072 } 073 074 return random; 075 } 076 077 /** {@inheritDoc} */ 078 @Override 079 public String toString() { 080 return "Box-Muller normalized Gaussian deviate [" + rng.toString() + "]"; 081 } 082 083 /** 084 * {@inheritDoc} 085 * 086 * @since 1.3 087 */ 088 @Override 089 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) { 090 return new BoxMullerNormalizedGaussianSampler(rng); 091 } 092 093 /** 094 * Create a new normalised Gaussian sampler. 095 * 096 * @param <S> Sampler type. 097 * @param rng Generator of uniformly distributed random numbers. 098 * @return the sampler 099 * @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}