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/Marsaglia_polar_method"> 023 * Marsaglia polar method</a> for sampling from a Gaussian distribution 024 * with mean 0 and standard deviation 1. 025 * This is a variation of the algorithm implemented in 026 * {@link BoxMullerNormalizedGaussianSampler}. 027 * 028 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p> 029 * 030 * @since 1.1 031 */ 032public class MarsagliaNormalizedGaussianSampler 033 implements NormalizedGaussianSampler, SharedStateContinuousSampler { 034 /** Next gaussian. */ 035 private double nextGaussian = Double.NaN; 036 /** Underlying source of randomness. */ 037 private final UniformRandomProvider rng; 038 039 /** 040 * @param rng Generator of uniformly distributed random numbers. 041 */ 042 public MarsagliaNormalizedGaussianSampler(UniformRandomProvider rng) { 043 this.rng = rng; 044 } 045 046 /** {@inheritDoc} */ 047 @Override 048 public double sample() { 049 if (Double.isNaN(nextGaussian)) { 050 // Rejection scheme for selecting a pair that lies within the unit circle. 051 while (true) { 052 // Generate a pair of numbers within [-1 , 1). 053 final double x = 2 * rng.nextDouble() - 1; 054 final double y = 2 * rng.nextDouble() - 1; 055 final double r2 = x * x + y * y; 056 057 if (r2 < 1 && r2 > 0) { 058 // Pair (x, y) is within unit circle. 059 final double alpha = Math.sqrt(-2 * Math.log(r2) / r2); 060 061 // Keep second element of the pair for next invocation. 062 nextGaussian = alpha * y; 063 064 // Return the first element of the generated pair. 065 return alpha * x; 066 } 067 068 // Pair is not within the unit circle: Generate another one. 069 } 070 } 071 072 // Use the second element of the pair (generated at the 073 // previous invocation). 074 final double r = nextGaussian; 075 076 // Both elements of the pair have been used. 077 nextGaussian = Double.NaN; 078 079 return r; 080 } 081 082 /** {@inheritDoc} */ 083 @Override 084 public String toString() { 085 return "Box-Muller (with rejection) normalized Gaussian deviate [" + rng.toString() + "]"; 086 } 087 088 /** 089 * {@inheritDoc} 090 * 091 * @since 1.3 092 */ 093 @Override 094 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) { 095 return new MarsagliaNormalizedGaussianSampler(rng); 096 } 097 098 /** 099 * 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}