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/RNG-30__sampling Commit: 3ece8fd401e0ed1a37db5267c32b3f7b5f43d219 Parents: 0996034 Author: Gilles <er...@apache.org> Authored: Sat Nov 12 16:33:59 2016 +0100 Committer: Gilles <er...@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() + "]"; + } +}
