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