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    /** Underlying source of randomness. */
042    private final UniformRandomProvider rng;
043
044    /**
045     * Creates a sampler instance.
046     *
047     * @param rng Generator of uniformly distributed random numbers.
048     * @param alpha Distribution first shape parameter.
049     * @param beta Distribution second shape parameter.
050     */
051    public ChengBetaSampler(UniformRandomProvider rng,
052                            double alpha,
053                            double beta) {
054        super(null);
055        this.rng = rng;
056        alphaShape = alpha;
057        betaShape = beta;
058    }
059
060    /** {@inheritDoc} */
061    @Override
062    public double sample() {
063        final double a = Math.min(alphaShape, betaShape);
064        final double b = Math.max(alphaShape, betaShape);
065
066        if (a > 1) {
067            return algorithmBB(a, b);
068        } else {
069            return algorithmBC(b, a);
070        }
071    }
072
073    /** {@inheritDoc} */
074    @Override
075    public String toString() {
076        return "Cheng Beta deviate [" + rng.toString() + "]";
077    }
078
079    /**
080     * Computes one sample using Cheng's BB algorithm, when \( \alpha \) and
081     * \( \beta \) are both larger than 1.
082     *
083     * @param a \( \min(\alpha, \beta) \).
084     * @param b \( \max(\alpha, \beta) \).
085     * @return a random sample.
086     */
087    private double algorithmBB(double a,
088                               double b) {
089        final double alpha = a + b;
090        final double beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
091        final double gamma = a + 1 / beta;
092
093        double r;
094        double w;
095        double t;
096        do {
097            final double u1 = rng.nextDouble();
098            final double u2 = rng.nextDouble();
099            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
100            w = a * Math.exp(v);
101            final double z = u1 * u1 * u2;
102            r = gamma * v - 1.3862944;
103            final double s = a + r - w;
104            if (s + 2.609438 >= 5 * z) {
105                break;
106            }
107
108            t = Math.log(z);
109            if (s >= t) {
110                break;
111            }
112        } while (r + alpha * (Math.log(alpha) - Math.log(b + w)) < t);
113
114        w = Math.min(w, Double.MAX_VALUE);
115
116        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
117    }
118
119    /**
120     * Computes one sample using Cheng's BB algorithm, when at least one of
121     * \( \alpha \) or \( \beta \) is smaller than 1.
122     *
123     * @param a \( \min(\alpha, \beta) \).
124     * @param b \( \max(\alpha, \beta) \).
125     * @return a random sample.
126     */
127    private double algorithmBC(double a,
128                               double b) {
129        final double alpha = a + b;
130        final double beta = 1 / b;
131        final double delta = 1 + a - b;
132        final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 0.777778);
133        final double k2 = 0.25 + (0.5 + 0.25 / delta) * b;
134
135        double w;
136        while (true) {
137            final double u1 = rng.nextDouble();
138            final double u2 = rng.nextDouble();
139            final double y = u1 * u2;
140            final double z = u1 * y;
141            if (u1 < 0.5) {
142                if (0.25 * u2 + z - y >= k1) {
143                    continue;
144                }
145            } else {
146                if (z <= 0.25) {
147                    final double v = beta * (Math.log(u1) - Math.log1p(-u1));
148                    w = a * Math.exp(v);
149                    break;
150                }
151
152                if (z >= k2) {
153                    continue;
154                }
155            }
156
157            final double v = beta * (Math.log(u1) - Math.log1p(-u1));
158            w = a * Math.exp(v);
159            if (alpha * (Math.log(alpha) - Math.log(b + w) + v) - 1.3862944 >= Math.log(z)) {
160                break;
161            }
162        }
163
164        w = Math.min(w, Double.MAX_VALUE);
165
166        return equals(a, alphaShape) ? w / (b + w) : b / (b + w);
167    }
168
169    /**
170     * @param a Value.
171     * @param b Value.
172     * @return {@code true} if {@code a} is equal to {@code b}.
173     */
174    private boolean equals(double a,
175                           double b) {
176        return Math.abs(a - b) <= Double.MIN_VALUE;
177    }
178}