commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [05/50] commons-rng git commit: RNG-30: Continuous distributions.
Date Mon, 21 Nov 2016 17:03:55 GMT
RNG-30: Continuous distributions.

Code copied and adapted from classes in "o.a.c.math4.distribution" package.


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

Branch: refs/heads/master
Commit: 3ece8fd401e0ed1a37db5267c32b3f7b5f43d219
Parents: 0996034
Author: Gilles <erans@apache.org>
Authored: Sat Nov 12 16:33:59 2016 +0100
Committer: Gilles <erans@apache.org>
Committed: Sat Nov 12 16:33:59 2016 +0100

----------------------------------------------------------------------
 .../AhrensDieterExponentialSampler.java         | 117 +++++++++++++
 .../AhrensDieterMarsagliaTsangGammaSampler.java | 139 +++++++++++++++
 .../distribution/BoxMullerLogNormalSampler.java |  59 +++++++
 .../sampling/distribution/ChengBetaSampler.java | 175 +++++++++++++++++++
 .../distribution/ContinuousUniformSampler.java  |  56 ++++++
 5 files changed, 546 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-rng/blob/3ece8fd4/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.java
new file mode 100644
index 0000000..8807ce8
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.java
@@ -0,0 +1,117 @@
+/*
+ * 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;
+import org.apache.commons.rng.sampling.AbstractContinuousSampler;
+
+/**
+ * Sampling from an exponential distribution.
+ *
+ * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html">Exponential
distribution (MathWorld)</a>
+ */
+public class AhrensDieterExponentialSampler extends AbstractContinuousSampler {
+    /**
+     * Table containing the constants
+     * \( q_i = sum_{j=1}^i (\ln 2)^j / j! = \ln 2 + (\ln 2)^2 / 2 + ... + (\ln 2)^i / i!
\)
+     * until the largest representable fraction below 1 is exceeded.
+     *
+     * Note that
+     * \( 1 = 2 - 1 = \exp(\ln 2) - 1 = sum_{n=1}^\infinity (\ln 2)^n / n! \)
+     * thus \( q_i \rightarrow 1 as i \rightarrow +\infinity \),
+     * so the higher \( i \), the closer we get to 1 (the series is not alternating).
+     *
+     * By trying, n = 16 in Java is enough to reach 1.
+     */
+    private static final double[] EXPONENTIAL_SA_QI = new double[16];
+    /** The mean of this distribution. */
+    private final double mean;
+
+    /**
+     * Initialize tables.
+     */
+    static {
+        /**
+         * Filling EXPONENTIAL_SA_QI table.
+         * Note that we don't want qi = 0 in the table.
+         */
+        final double ln2 = Math.log(2);
+        double qi = 0;
+
+        for (int i = 0; i < EXPONENTIAL_SA_QI.length; i++) {
+            qi += Math.pow(ln2, i + 1) / InternalUtils.factorial(i + 1);
+            EXPONENTIAL_SA_QI[i] = qi;
+        }
+    }
+
+    /**
+     * @param rng Generator of uniformly distributed random numbers.
+     * @param mean Mean of this distribution.
+     */
+    public AhrensDieterExponentialSampler(UniformRandomProvider rng,
+                                          double mean) {
+        super(rng);
+        this.mean = mean;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        // Step 1:
+        double a = 0;
+        double u = nextUniform();
+
+        // Step 2 and 3:
+        while (u < 0.5) {
+            a += EXPONENTIAL_SA_QI[0];
+            u *= 2;
+        }
+
+        // Step 4 (now u >= 0.5):
+        u += u - 1;
+
+        // Step 5:
+        if (u <= EXPONENTIAL_SA_QI[0]) {
+            return mean * (a + u);
+        }
+
+        // Step 6:
+        int i = 0; // Should be 1, be we iterate before it in while using 0.
+        double u2 = nextUniform();
+        double umin = u2;
+
+        // Step 7 and 8:
+        do {
+            ++i;
+            u2 = nextUniform();
+
+            if (u2 < umin) {
+                umin = u2;
+            }
+
+            // Step 8:
+        } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX]
= 1.
+
+        return mean * (a + umin * EXPONENTIAL_SA_QI[0]);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "[Ahrens-Dieter method for exponential distribution " + super.toString() +
"]";
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-rng/blob/3ece8fd4/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java
new file mode 100644
index 0000000..14e6382
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.java
@@ -0,0 +1,139 @@
+/*
+ * 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;
+import org.apache.commons.rng.sampling.AbstractContinuousSampler;
+
+/**
+ * <p>
+ * Sampling from the <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma
distribution</a>.
+ * <ul>
+ *  <li>
+ *   For {@code 0 < shape < 1}:
+ *   <blockquote>
+ *    Ahrens, J. H. and Dieter, U.,
+ *    <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
+ *    Computing, 12, 223-246, 1974.
+ *   </blockquote>
+ *  </li>
+ *  <li>
+ *  For {@code shape >= 1}:
+ *   <blockquote>
+ *   Marsaglia and Tsang, <i>A Simple Method for Generating
+ *   Gamma Variables.</i> ACM Transactions on Mathematical Software,
+ *   Volume 26 Issue 3, September, 2000.
+ *   </blockquote>
+ *  </li>
+ * </ul>
+ * </p>
+ */
+public class AhrensDieterMarsagliaTsangGammaSampler extends AbstractContinuousSampler {
+    /** The shape parameter. */
+    private final double theta;
+    /** The alpha parameter. */
+    private final double alpha;
+    /** Gaussian sampling. */
+    private final BoxMullerGaussianSampler gaussian;
+
+    /**
+     * @param rng Generator of uniformly distributed random numbers.
+     * @param alpha Alpha parameter of the distribution.
+     * @param theta Theta parameter of the distribution.
+     */
+    public AhrensDieterMarsagliaTsangGammaSampler(UniformRandomProvider rng,
+                                                  double alpha,
+                                                  double theta) {
+        super(rng);
+        this.alpha = alpha;
+        this.theta = theta;
+        gaussian = new BoxMullerGaussianSampler(rng, 0, 1);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        if (theta < 1) {
+            // [1]: p. 228, Algorithm GS.
+
+            while (true) {
+                // Step 1:
+                final double u = nextUniform();
+                final double bGS = 1 + theta / Math.E;
+                final double p = bGS * u;
+
+                if (p <= 1) {
+                    // Step 2:
+
+                    final double x = Math.pow(p, 1 / theta);
+                    final double u2 = nextUniform();
+
+                    if (u2 > Math.exp(-x)) {
+                        // Reject.
+                        continue;
+                    } else {
+                        return alpha * x;
+                    }
+                } else {
+                    // Step 3:
+
+                    final double x = -1 * Math.log((bGS - p) / theta);
+                    final double u2 = nextUniform();
+
+                    if (u2 > Math.pow(x, theta - 1)) {
+                        // Reject.
+                        continue;
+                    } else {
+                        return alpha * x;
+                    }
+                }
+            }
+        }
+
+        // Now theta >= 1.
+
+        final double d = theta - 0.333333333333333333;
+        final double c = 1 / (3 * Math.sqrt(d));
+
+        while (true) {
+            final double x = gaussian.sample();
+            final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);
+
+            if (v <= 0) {
+                continue;
+            }
+
+            final double x2 = x * x;
+            final double u = nextUniform();
+
+            // Squeeze.
+            if (u < 1 - 0.0331 * x2 * x2) {
+                return alpha * d * v;
+            }
+
+            if (Math.log(u) < 0.5 * x2 + d * (1 - v + Math.log(v))) {
+                return alpha * d * v;
+            }
+        }
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "[Ahrens-Dieter-Marsaglia-Tsang method for Gamma distribution " + super.toString()
+ "]";
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-rng/blob/3ece8fd4/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerLogNormalSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerLogNormalSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerLogNormalSampler.java
new file mode 100644
index 0000000..4cda8d4
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/BoxMullerLogNormalSampler.java
@@ -0,0 +1,59 @@
+/*
+ * 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;
+import org.apache.commons.rng.sampling.AbstractContinuousSampler;
+
+/**
+ * <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform">
+ * Box-Muller algorithm</a> for sampling from a Log normal distribution.
+ */
+public class BoxMullerLogNormalSampler extends AbstractContinuousSampler {
+    /** Scale. */
+    private final double scale;
+    /** Shape. */
+    private final double shape;
+    /** Gaussian sampling. */
+    private final BoxMullerGaussianSampler gaussian;
+
+    /**
+     * @param rng Generator of uniformly distributed random numbers.
+     * @param scale Scale of the Log normal distribution.
+     * @param shape Shape of the Log normal distribution.
+     */
+    public BoxMullerLogNormalSampler(UniformRandomProvider rng,
+                                     double scale,
+                                     double shape) {
+        super(null); // Not used.
+        this.scale = scale;
+        this.shape = shape;
+        gaussian = new BoxMullerGaussianSampler(rng, 0, 1);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        return Math.exp(scale + shape * gaussian.sample());
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "[Box-Muller method for Log normal distribution " + gaussian.toString() +
"]";
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-rng/blob/3ece8fd4/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
new file mode 100644
index 0000000..d0bfd41
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ChengBetaSampler.java
@@ -0,0 +1,175 @@
+/*
+ * 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;
+import org.apache.commons.rng.sampling.AbstractContinuousSampler;
+
+
+/**
+ * Utility class implementing Cheng's algorithms for beta distribution sampling.
+ *
+ * <blockquote>
+ * <pre>
+ * R. C. H. Cheng,
+ * "Generating beta variates with nonintegral shape parameters",
+ * Communications of the ACM, 21, 317-322, 1978.
+ * </pre>
+ * </blockquote>
+ *
+ * @since 1.0
+ */
+public class ChengBetaSampler extends AbstractContinuousSampler {
+    /** First shape parameter. */
+    private final double alphaShape;
+    /** Second shape parameter. */
+    private final double betaShape;
+
+    /**
+     * Creates a sampler instance.
+     *
+     * @param rng Generator of uniformly distributed random numbers.
+     * @param alpha Distribution first shape parameter.
+     * @param beta Distribution second shape parameter.
+     */
+    public ChengBetaSampler(UniformRandomProvider rng,
+                            double alpha,
+                            double beta) {
+        super(rng);
+        alphaShape = alpha;
+        betaShape = beta;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        final double a = Math.min(alphaShape, betaShape);
+        final double b = Math.max(alphaShape, betaShape);
+
+        if (a > 1) {
+            return algorithmBB(a, b);
+        } else {
+            return algorithmBC(b, a);
+        }
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "[Cheng method for Beta distribution " + super.toString() + "]";
+    }
+
+    /**
+     * Computes one sample using Cheng's BB algorithm, when \( \alpha \) and
+     * \( \beta \) are both larger than 1.
+     *
+     * @param a \( \min(\alpha, \beta) \).
+     * @param b \( \max(\alpha, \beta) \).
+     * @return a random sample.
+     */
+    private double algorithmBB(double a,
+                               double b) {
+        final double alpha = a + b;
+        final double beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
+        final double gamma = a + 1 / beta;
+
+        double r;
+        double w;
+        double t;
+        do {
+            final double u1 = nextUniform();
+            final double u2 = nextUniform();
+            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
+            w = a * Math.exp(v);
+            final double z = u1 * u1 * u2;
+            r = gamma * v - 1.3862944;
+            final double s = a + r - w;
+            if (s + 2.609438 >= 5 * z) {
+                break;
+            }
+
+            t = Math.log(z);
+            if (s >= t) {
+                break;
+            }
+        } while (r + alpha * (Math.log(alpha) - Math.log(b + w)) < t);
+
+        w = Math.min(w, Double.MAX_VALUE);
+
+        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
+    }
+
+    /**
+     * Computes one sample using Cheng's BB algorithm, when at least one of
+     * \( \alpha \) or \( \beta \) is smaller than 1.
+     *
+     * @param a \( \min(\alpha, \beta) \).
+     * @param b \( \max(\alpha, \beta) \).
+     * @return a random sample.
+     */
+    private double algorithmBC(double a,
+                               double b) {
+        final double alpha = a + b;
+        final double beta = 1 / b;
+        final double delta = 1 + a - b;
+        final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 0.777778);
+        final double k2 = 0.25 + (0.5 + 0.25 / delta) * b;
+
+        double w;
+        while (true) {
+            final double u1 = nextUniform();
+            final double u2 = nextUniform();
+            final double y = u1 * u2;
+            final double z = u1 * y;
+            if (u1 < 0.5) {
+                if (0.25 * u2 + z - y >= k1) {
+                    continue;
+                }
+            } else {
+                if (z <= 0.25) {
+                    final double v = beta * (Math.log(u1) - Math.log1p(-u1));
+                    w = a * Math.exp(v);
+                    break;
+                }
+
+                if (z >= k2) {
+                    continue;
+                }
+            }
+
+            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
+            w = a * Math.exp(v);
+            if (alpha * (Math.log(alpha) - Math.log(b + w) + v) - 1.3862944 >= Math.log(z))
{
+                break;
+            }
+        }
+
+        w = Math.min(w, Double.MAX_VALUE);
+
+        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
+    }
+
+    /**
+     * @param a Value.
+     * @param b Value.
+     * @return {@code true} if {@code a} is equal to {@code b}.
+     */
+    private boolean equals(double a,
+                           double b) {
+        return Math.abs(a - b) <= Double.MIN_VALUE;
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-rng/blob/3ece8fd4/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ContinuousUniformSampler.java
----------------------------------------------------------------------
diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ContinuousUniformSampler.java
b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ContinuousUniformSampler.java
new file mode 100644
index 0000000..a77cddc
--- /dev/null
+++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ContinuousUniformSampler.java
@@ -0,0 +1,56 @@
+/*
+ * 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;
+import org.apache.commons.rng.sampling.AbstractContinuousSampler;
+
+/**
+ * Sampling from a uniform distribution.
+ */
+public class ContinuousUniformSampler extends AbstractContinuousSampler {
+    /** Lower bound. */
+    private final double lo;
+    /** Higher bound. */
+    private final double hi;
+
+    /**
+     * @param rng Generator of uniformly distributed random numbers.
+     * @param lo Lower bound.
+     * @param hi Higher bound.
+     */
+    public ContinuousUniformSampler(UniformRandomProvider rng,
+                                    double lo,
+                                    double hi) {
+        super(rng);
+        this.lo = lo;
+        this.hi = hi;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public double sample() {
+        final double u = nextUniform();
+        return u * hi + (1 - u) * lo;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public String toString() {
+        return "[Uniform distribution " + super.toString() + "]";
+    }
+}


Mime
View raw message