BoxMullerNormalizedGaussianSampler.java

  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. import org.apache.commons.rng.UniformRandomProvider;

  19. /**
  20.  * <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform">
  21.  * Box-Muller algorithm</a> for sampling from Gaussian distribution with
  22.  * mean 0 and standard deviation 1.
  23.  *
  24.  * <p>Sampling uses:</p>
  25.  *
  26.  * <ul>
  27.  *   <li>{@link UniformRandomProvider#nextDouble()}
  28.  *   <li>{@link UniformRandomProvider#nextLong()}
  29.  * </ul>
  30.  *
  31.  * @since 1.1
  32.  */
  33. public class BoxMullerNormalizedGaussianSampler
  34.     implements NormalizedGaussianSampler, SharedStateContinuousSampler {
  35.     /** Next gaussian. */
  36.     private double nextGaussian = Double.NaN;
  37.     /** Underlying source of randomness. */
  38.     private final UniformRandomProvider rng;

  39.     /**
  40.      * Create an instance.
  41.      *
  42.      * @param rng Generator of uniformly distributed random numbers.
  43.      */
  44.     public BoxMullerNormalizedGaussianSampler(UniformRandomProvider rng) {
  45.         this.rng = rng;
  46.     }

  47.     /** {@inheritDoc} */
  48.     @Override
  49.     public double sample() {
  50.         final double random;
  51.         if (Double.isNaN(nextGaussian)) {
  52.             // Generate a pair of Gaussian numbers.

  53.             // Avoid zero for the uniform deviate y.
  54.             // The extreme tail of the sample is:
  55.             // y = 2^-53
  56.             // r = 8.57167
  57.             final double x = rng.nextDouble();
  58.             final double y = InternalUtils.makeNonZeroDouble(rng.nextLong());
  59.             final double alpha = 2 * Math.PI * x;
  60.             final double r = Math.sqrt(-2 * Math.log(y));

  61.             // Return the first element of the generated pair.
  62.             random = r * Math.cos(alpha);

  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.             // Both elements of the pair have been used.
  70.             nextGaussian = Double.NaN;
  71.         }

  72.         return random;
  73.     }

  74.     /** {@inheritDoc} */
  75.     @Override
  76.     public String toString() {
  77.         return "Box-Muller normalized Gaussian deviate [" + rng.toString() + "]";
  78.     }

  79.     /**
  80.      * {@inheritDoc}
  81.      *
  82.      * @since 1.3
  83.      */
  84.     @Override
  85.     public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
  86.         return new BoxMullerNormalizedGaussianSampler(rng);
  87.     }

  88.     /**
  89.      * Create a new normalised Gaussian sampler.
  90.      *
  91.      * @param <S> Sampler type.
  92.      * @param rng Generator of uniformly distributed random numbers.
  93.      * @return the sampler
  94.      * @since 1.3
  95.      */
  96.     @SuppressWarnings("unchecked")
  97.     public static <S extends NormalizedGaussianSampler & SharedStateContinuousSampler> S
  98.             of(UniformRandomProvider rng) {
  99.         return (S) new BoxMullerNormalizedGaussianSampler(rng);
  100.     }
  101. }