commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [07/14] commons-rng git commit: RNG-36: Variation of the Box-Muller algorithm for Gaussian sampling.
Date Thu, 09 Feb 2017 00:22:16 GMT
RNG-36: Variation of the Box-Muller algorithm for Gaussian sampling.


Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/c14e09ca
Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/c14e09ca
Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/c14e09ca

Branch: refs/heads/master
Commit: c14e09caa4161dc7cdbaa9fae3dbbd3cf02b70e3
Parents: 94f0843
Author: Gilles <erans@apache.org>
Authored: Thu Jan 26 23:54:52 2017 +0100
Committer: Gilles <erans@apache.org>
Committed: Thu Jan 26 23:54:52 2017 +0100

----------------------------------------------------------------------
 ...rWithRejectionNormalizedGaussianSampler.java | 88 ++++++++++++++++++++
 .../distribution/ContinuousSamplersList.java    |  4 +
 2 files changed, 92 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-rng/blob/c14e09ca/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerWithRejectionNormalizedGaussianSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerWithRejectionNormalizedGaussianSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerWithRejectionNormalizedGaussianSampler.java
new file mode 100644
index 0000000..58d1b1f
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerWithRejectionNormalizedGaussianSampler.java
@@ -0,0 +1,88 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.rng.sampling.distribution;
+
+import org.apache.commons.rng.UniformRandomProvider;
+
+/**
+ * <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform">
+ * Box-Muller algorithm</a> for sampling from Gaussian distribution with
+ * mean 0 and standard deviation 1.
+ * This is a variation, suggested in <a href="http://haifux.org/lectures/79/random.pdf">
+ * this presentation</a> (page 39), of the algorithm implemented in
+ * {@link BoxMullerNormalizedGaussianSampler}, 
+ *
+ * @since 1.1
+ */
+public class BoxMullerWithRejectionNormalizedGaussianSampler
+    extends SamplerBase
+    implements NormalizedGaussianSampler {
+    /** Next gaussian. */
+    private double nextGaussian = Double.NaN;
+
+    /**
+     * @param rng Generator of uniformly distributed random numbers.
+     */
+    public BoxMullerWithRejectionNormalizedGaussianSampler(UniformRandomProvider rng) {
+        super(rng);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        final double random;
+        if (Double.isNaN(nextGaussian)) {
+            // Rejection scheme for selecting a pair that lies within the unit circle.
+            SAMPLE: while (true) {
+                // Generate a pair of numbers within [-1 , 1).
+                final double x = 2 * nextDouble() - 1;
+                final double y = 2 * nextDouble() - 1;
+                final double r2 = x * x + y * y;
+
+                if (r2 < 1) {
+                    // Pair (x, y) is within unit circle.
+
+                    final double r = Math.sqrt(r2);
+                    final double alpha = 2 * Math.sqrt(-Math.log(r)) / r;
+
+                    // Return the first element of the generated pair.
+                    random = alpha * x;
+
+                    // Keep second element of the pair for next invocation.
+                    nextGaussian = alpha * y;
+                    break SAMPLE;
+                }
+                // Pair is not within the unit circle: Generate another one.
+            }
+        } else {
+            // Use the second element of the pair (generated at the
+            // previous invocation).
+            random = nextGaussian;
+
+            // Both elements of the pair have been used.
+            nextGaussian = Double.NaN;
+        }
+
+        return random;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "Box-Muller (with rejection) normalized Gaussian deviate [" + super.toString()
+ "]";
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-rng/blob/c14e09ca/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
index c24e786..aa0e458 100644
--- a/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
+++ b/commons-rng-sampling/src/test/java/org/apache/commons/rng/sampling/distribution/ContinuousSamplersList.java
@@ -48,6 +48,10 @@ public class ContinuousSamplersList {
             add(LIST, new org.apache.commons.math3.distribution.NormalDistribution(meanNormal,
sigmaNormal),
                 new GaussianSampler(new BoxMullerNormalizedGaussianSampler(RandomSource.create(RandomSource.MT)),
                                     meanNormal, sigmaNormal));
+            // Gaussian ("Box-Muller" with rejection).
+            add(LIST, new org.apache.commons.math3.distribution.NormalDistribution(meanNormal,
sigmaNormal),
+                new GaussianSampler(new BoxMullerWithRejectionNormalizedGaussianSampler(RandomSource.create(RandomSource.MT)),
+                                    meanNormal, sigmaNormal));
 
             // Beta ("inverse method").
             final double alphaBeta = 4.3;


Mime
View raw message