Repository: commons-math Updated Branches: refs/heads/MATH-1153 9ba0a1cad -> 84a642266
Use static internal class and method for Cheng's sampling. The static components allow an implementation closer to the original patch contribution, avoiding weird variables renamings. Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/84a64226 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/84a64226 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/84a64226 Branch: refs/heads/MATH-1153 Commit: 84a642266dcebf24ed02670e26d0e41153af4f9c Parents: 9ba0a1c Author: Luc Maisonobe <l...@apache.org> Authored: Wed Dec 17 11:05:56 2014 +0100 Committer: Luc Maisonobe <l...@apache.org> Committed: Wed Dec 17 11:05:56 2014 +0100 ---------------------------------------------------------------------- .../math3/distribution/BetaDistribution.java | 170 +++++++++++-------- 1 file changed, 101 insertions(+), 69 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/84a64226/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java b/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java index 458fe23..d7a3f28 100644 --- a/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java +++ b/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java @@ -270,91 +270,123 @@ public class BetaDistribution extends AbstractRealDistribution { */ @Override public double sample() { - if (FastMath.min(alpha, beta) > 1) { - return algorithmBB(); - } else { - return algorithmBC(); - } + return ChengBetaSampler.sample(random, alpha, beta); } - /** Returns one sample using Cheng's BB algorithm, when both α and β are greater than 1. - * @return sampled value + /** Utility class implementing Cheng's algorithms for beta distribution sampling. + * <p> + * R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.". + * Communications of the ACM, 21, 317â322, 1978. + * </p> + * @since 3.4 */ - private double algorithmBB() { - final double a = FastMath.min(alpha, beta); - final double b = FastMath.max(alpha, beta); - final double newAlpha = a + b; - final double newBeta = FastMath.sqrt((newAlpha - 2.) / (2. * a * b - newAlpha)); - final double gamma = a + 1. / newBeta; - - double r; - double w; - double t; - do { - final double u1 = random.nextDouble(); - final double u2 = random.nextDouble(); - final double v = newBeta * FastMath.log(u1 / (1. - u1)); - w = a * FastMath.exp(v); - final double newZ = u1 * u1 * u2; - r = gamma * v - 1.3862944; - final double s = a + r - w; - if (s + 2.609438 >= 5 * newZ) { - break; - } + private static class ChengBetaSampler { - t = FastMath.log(newZ); - if (s >= t) { - break; - } - } while (r + newAlpha * FastMath.log(newAlpha / (b + w)) < t); + /** Private constructor for a utility method. + */ + private ChengBetaSampler() { + } - w = FastMath.min(w, Double.MAX_VALUE); - return Precision.equals(a, alpha) ? w / (b + w) : b / (b + w); - } + /** Returns one sample using Cheng's BB algorithm, when both α and β are greater than 1. + * @param random random generator to use + * @param a0 distribution first shape parameter + * @param b0 distribution second shape parameter + * @return sampled value + */ + public static double algorithmBB(final RandomGenerator random, final double a0, final double b0) { + final double a = FastMath.min(a0, b0); + final double b = FastMath.max(a0, b0); + final double alpha = a + b; + final double beta = FastMath.sqrt((alpha - 2.) / (2. * a * b - alpha)); + final double gamma = a + 1. / beta; - /** Returns one sample using Cheng's BC algorithm, when at least one of α and β is smaller than 1. - * @return sampled value - */ - private double algorithmBC() { - final double a = FastMath.max(alpha, beta); - final double b = FastMath.min(alpha, beta); - final double newAlpha = a + b; - final double newBeta = 1. / b; - final double delta = 1. + a - b; - final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * newBeta - 0.777778); - final double k2 = 0.25 + (0.5 + 0.25 / delta) * b; - - double w; - for (;;) { - final double u1 = random.nextDouble(); - final double u2 = random.nextDouble(); - final double y = u1 * u2; - final double newZ = u1 * y; - if (u1 < 0.5) { - if (0.25 * u2 + newZ - y >= k1) { - continue; + double r; + double w; + double t; + do { + final double u1 = random.nextDouble(); + final double u2 = random.nextDouble(); + final double v = beta * FastMath.log(u1 / (1. - u1)); + w = a * FastMath.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; } - } else { - if (newZ <= 0.25) { - final double v = newBeta * FastMath.log(u1 / (1. - u1)); - w = a * FastMath.exp(v); + + t = FastMath.log(z); + if (s >= t) { break; } + } while (r + alpha * FastMath.log(alpha / (b + w)) < t); + + w = FastMath.min(w, Double.MAX_VALUE); + return Precision.equals(a, a0) ? w / (b + w) : b / (b + w); + } + + /** Returns one sample using Cheng's BC algorithm, when at least one of α and β is smaller than 1. + * @param random random generator to use + * @param a0 distribution first shape parameter + * @param b0 distribution second shape parameter + * @return sampled value + */ + public static double algorithmBC(final RandomGenerator random, final double a0, final double b0) { + final double a = FastMath.max(a0, b0); + final double b = FastMath.min(a0, b0); + 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; + for (;;) { + final double u1 = random.nextDouble(); + final double u2 = random.nextDouble(); + 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 * FastMath.log(u1 / (1. - u1)); + w = a * FastMath.exp(v); + break; + } + + if (z >= k2) { + continue; + } + } - if (newZ >= k2) { - continue; + final double v = beta * FastMath.log(u1 / (1. - u1)); + w = a * FastMath.exp(v); + if (alpha * (FastMath.log(alpha / (b + w)) + v) - 1.3862944 >= FastMath.log(z)) { + break; } } - final double v = newBeta * FastMath.log(u1 / (1. - u1)); - w = a * FastMath.exp(v); - if (newAlpha * (FastMath.log(newAlpha / (b + w)) + v) - 1.3862944 >= FastMath.log(newZ)) { - break; + w = FastMath.min(w, Double.MAX_VALUE); + return Precision.equals(a, a0) ? w / (b + w) : b / (b + w); + } + + /** Returns one sample using Cheng's algorithms. + * @param random random generator to use + * @param alpha distribution first shape parameter + * @param beta distribution second shape parameter + * @return sampled value + */ + public static double sample(final RandomGenerator random, final double alpha, final double beta) { + if (FastMath.min(alpha, beta) > 1) { + return algorithmBB(random, alpha, beta); + } else { + return algorithmBC(random, alpha, beta); } } - w = FastMath.min(w, Double.MAX_VALUE); - return Precision.equals(a, alpha) ? w / (b + w) : b / (b + w); } }
