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 * Create an instance. 041 * 042 * @param rng Generator of uniformly distributed random numbers. 043 */ 044 public MarsagliaNormalizedGaussianSampler(UniformRandomProvider rng) { 045 this.rng = rng; 046 } 047 048 /** {@inheritDoc} */ 049 @Override 050 public double sample() { 051 if (Double.isNaN(nextGaussian)) { 052 // Rejection scheme for selecting a pair that lies within the unit circle. 053 while (true) { 054 // Generate a pair of numbers within [-1 , 1). 055 final double x = 2 * rng.nextDouble() - 1; 056 final double y = 2 * rng.nextDouble() - 1; 057 final double r2 = x * x + y * y; 058 059 if (r2 < 1 && r2 > 0) { 060 // Pair (x, y) is within unit circle. 061 final double alpha = Math.sqrt(-2 * Math.log(r2) / r2); 062 063 // Keep second element of the pair for next invocation. 064 nextGaussian = alpha * y; 065 066 // Return the first element of the generated pair. 067 return alpha * x; 068 } 069 070 // Pair is not within the unit circle: Generate another one. 071 } 072 } 073 074 // Use the second element of the pair (generated at the 075 // previous invocation). 076 final double r = nextGaussian; 077 078 // Both elements of the pair have been used. 079 nextGaussian = Double.NaN; 080 081 return r; 082 } 083 084 /** {@inheritDoc} */ 085 @Override 086 public String toString() { 087 return "Box-Muller (with rejection) normalized Gaussian deviate [" + rng.toString() + "]"; 088 } 089 090 /** 091 * {@inheritDoc} 092 * 093 * @since 1.3 094 */ 095 @Override 096 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) { 097 return new MarsagliaNormalizedGaussianSampler(rng); 098 } 099 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}