001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.sampling.distribution;
018
019import org.apache.commons.rng.UniformRandomProvider;
020
021/**
022 * Utility class implementing Cheng's algorithms for beta distribution sampling.
023 *
024 * <blockquote>
025 * <pre>
026 * R. C. H. Cheng,
027 * "Generating beta variates with nonintegral shape parameters",
028 * Communications of the ACM, 21, 317-322, 1978.
029 * </pre>
030 * </blockquote>
031 *
032 * @since 1.0
033 */
034public class ChengBetaSampler
035    extends SamplerBase
036    implements ContinuousSampler {
037    /** First shape parameter. */
038    private final double alphaShape;
039    /** Second shape parameter. */
040    private final double betaShape;
041
042    /**
043     * Creates a sampler instance.
044     *
045     * @param rng Generator of uniformly distributed random numbers.
046     * @param alpha Distribution first shape parameter.
047     * @param beta Distribution second shape parameter.
048     */
049    public ChengBetaSampler(UniformRandomProvider rng,
050                            double alpha,
051                            double beta) {
052        super(rng);
053        alphaShape = alpha;
054        betaShape = beta;
055    }
056
057    /** {@inheritDoc} */
058    @Override
059    public double sample() {
060        final double a = Math.min(alphaShape, betaShape);
061        final double b = Math.max(alphaShape, betaShape);
062
063        if (a > 1) {
064            return algorithmBB(a, b);
065        } else {
066            return algorithmBC(b, a);
067        }
068    }
069
070    /** {@inheritDoc} */
071    @Override
072    public String toString() {
073        return "Cheng Beta deviate [" + super.toString() + "]";
074    }
075
076    /**
077     * Computes one sample using Cheng's BB algorithm, when \( \alpha \) and
078     * \( \beta \) are both larger than 1.
079     *
080     * @param a \( \min(\alpha, \beta) \).
081     * @param b \( \max(\alpha, \beta) \).
082     * @return a random sample.
083     */
084    private double algorithmBB(double a,
085                               double b) {
086        final double alpha = a + b;
087        final double beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
088        final double gamma = a + 1 / beta;
089
090        double r;
091        double w;
092        double t;
093        do {
094            final double u1 = nextDouble();
095            final double u2 = nextDouble();
096            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
097            w = a * Math.exp(v);
098            final double z = u1 * u1 * u2;
099            r = gamma * v - 1.3862944;
100            final double s = a + r - w;
101            if (s + 2.609438 >= 5 * z) {
102                break;
103            }
104
105            t = Math.log(z);
106            if (s >= t) {
107                break;
108            }
109        } while (r + alpha * (Math.log(alpha) - Math.log(b + w)) < t);
110
111        w = Math.min(w, Double.MAX_VALUE);
112
113        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
114    }
115
116    /**
117     * Computes one sample using Cheng's BB algorithm, when at least one of
118     * \( \alpha \) or \( \beta \) is smaller than 1.
119     *
120     * @param a \( \min(\alpha, \beta) \).
121     * @param b \( \max(\alpha, \beta) \).
122     * @return a random sample.
123     */
124    private double algorithmBC(double a,
125                               double b) {
126        final double alpha = a + b;
127        final double beta = 1 / b;
128        final double delta = 1 + a - b;
129        final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 0.777778);
130        final double k2 = 0.25 + (0.5 + 0.25 / delta) * b;
131
132        double w;
133        while (true) {
134            final double u1 = nextDouble();
135            final double u2 = nextDouble();
136            final double y = u1 * u2;
137            final double z = u1 * y;
138            if (u1 < 0.5) {
139                if (0.25 * u2 + z - y >= k1) {
140                    continue;
141                }
142            } else {
143                if (z <= 0.25) {
144                    final double v = beta * (Math.log(u1) - Math.log1p(-u1));
145                    w = a * Math.exp(v);
146                    break;
147                }
148
149                if (z >= k2) {
150                    continue;
151                }
152            }
153
154            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
155            w = a * Math.exp(v);
156            if (alpha * (Math.log(alpha) - Math.log(b + w) + v) - 1.3862944 >= Math.log(z)) {
157                break;
158            }
159        }
160
161        w = Math.min(w, Double.MAX_VALUE);
162
163        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
164    }
165
166    /**
167     * @param a Value.
168     * @param b Value.
169     * @return {@code true} if {@code a} is equal to {@code b}.
170     */
171    private boolean equals(double a,
172                           double b) {
173        return Math.abs(a - b) <= Double.MIN_VALUE;
174    }
175}