Source code

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 * Sampling from a <a href="http://en.wikipedia.org/wiki/Beta_distribution">beta distribution</a>.
023 *
024 * <p>Uses Cheng's algorithms for beta distribution sampling:</p>
025 *
026 * <blockquote>
027 * <pre>
028 * R. C. H. Cheng,
029 * "Generating beta variates with nonintegral shape parameters",
030 * Communications of the ACM, 21, 317-322, 1978.
031 * </pre>
032 * </blockquote>
033 *
034 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
035 *
036 * @since 1.0
037 */
038public class ChengBetaSampler
039    extends SamplerBase
040    implements SharedStateContinuousSampler {
041    /** Natural logarithm of 4. */
042    private static final double LN_4 = Math.log(4.0);
043
044    /** The appropriate beta sampler for the parameters. */
045    private final SharedStateContinuousSampler delegate;
046
047    /**
048     * Base class to implement Cheng's algorithms for the beta distribution.
049     */
050    private abstract static class BaseChengBetaSampler
051            implements SharedStateContinuousSampler {
052        /**
053         * Flag set to true if {@code a} is the beta distribution alpha shape parameter.
054         * Otherwise {@code a} is the beta distribution beta shape parameter.
055         *
056         * <p>From the original Cheng paper this is equal to the result of: {@code a == a0}.</p>
057         */
058        protected final boolean aIsAlphaShape;
059        /**
060         * First shape parameter.
061         * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
062         */
063        protected final double a;
064        /**
065         * Second shape parameter.
066         * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
067         */
068        protected final double b;
069        /** Underlying source of randomness. */
070        protected final UniformRandomProvider rng;
071        /**
072         * The algorithm alpha factor. This is not the beta distribution alpha shape parameter.
073         * It is the sum of the two shape parameters ({@code a + b}.
074         */
075        protected final double alpha;
076        /** The logarithm of the alpha factor. */
077        protected final double logAlpha;
078
079        /**
080         * @param rng Generator of uniformly distributed random numbers.
081         * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
082         * @param a Distribution first shape parameter.
083         * @param b Distribution second shape parameter.
084         */
085        BaseChengBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
086            this.rng = rng;
087            this.aIsAlphaShape = aIsAlphaShape;
088            this.a = a;
089            this.b = b;
090            alpha = a + b;
091            logAlpha = Math.log(alpha);
092        }
093
094        /**
095         * @param rng Generator of uniformly distributed random numbers.
096         * @param source Source to copy.
097         */
098        private BaseChengBetaSampler(UniformRandomProvider rng,
099                                     BaseChengBetaSampler source) {
100            this.rng = rng;
101            aIsAlphaShape = source.aIsAlphaShape;
102            a = source.a;
103            b = source.b;
104            // Compute algorithm factors.
105            alpha = source.alpha;
106            logAlpha = source.logAlpha;
107        }
108
109        /** {@inheritDoc} */
110        @Override
111        public String toString() {
112            return "Cheng Beta deviate [" + rng.toString() + "]";
113        }
114
115        /**
116         * Compute the sample result X.
117         *
118         * <blockquote>
119         * If a == a0 deliver X = W/(b + W); otherwise deliver X = b/(b + W).
120         * </blockquote>
121         *
122         * <p>The finalisation step is shared between the BB and BC algorithm (as step 5 of the
123         * BB algorithm and step 6 of the BC algorithm).</p>
124         *
125         * @param w Algorithm value W.
126         * @return the sample value
127         */
128        protected double computeX(double w) {
129            // Avoid (infinity / infinity) producing NaN
130            final double tmp = Math.min(w, Double.MAX_VALUE);
131            return aIsAlphaShape ? tmp / (b + tmp) : b / (b + tmp);
132        }
133    }
134
135    /**
136     * Computes one sample using Cheng's BB algorithm, when beta distribution {@code alpha} and
137     * {@code beta} shape parameters are both larger than 1.
138     */
139    private static class ChengBBBetaSampler extends BaseChengBetaSampler {
140        /** 1 + natural logarithm of 5. */
141        private static final double LN_5_P1 = 1 + Math.log(5.0);
142
143        /** The algorithm beta factor. This is not the beta distribution beta shape parameter. */
144        private final double beta;
145        /** The algorithm gamma factor. */
146        private final double gamma;
147
148        /**
149         * @param rng Generator of uniformly distributed random numbers.
150         * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
151         * @param a min(alpha, beta) shape parameter.
152         * @param b max(alpha, beta) shape parameter.
153         */
154        ChengBBBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
155            super(rng, aIsAlphaShape, a, b);
156            beta = Math.sqrt((alpha - 2) / (2 * a * b - alpha));
157            gamma = a + 1 / beta;
158        }
159
160        /**
161         * @param rng Generator of uniformly distributed random numbers.
162         * @param source Source to copy.
163         */
164        private ChengBBBetaSampler(UniformRandomProvider rng,
165                                   ChengBBBetaSampler source) {
166            super(rng, source);
167            // Compute algorithm factors.
168            beta = source.beta;
169            gamma = source.gamma;
170        }
171
172        @Override
173        public double sample() {
174            double r;
175            double w;
176            double t;
177            do {
178                // Step 1:
179                final double u1 = rng.nextDouble();
180                final double u2 = rng.nextDouble();
181                final double v = beta * (Math.log(u1) - Math.log1p(-u1));
182                w = a * Math.exp(v);
183                final double z = u1 * u1 * u2;
184                r = gamma * v - LN_4;
185                final double s = a + r - w;
186                // Step 2:
187                if (s + LN_5_P1 >= 5 * z) {
188                    break;
189                }
190
191                // Step 3:
192                t = Math.log(z);
193                if (s >= t) {
194                    break;
195                }
196                // Step 4:
197            } while (r + alpha * (logAlpha - Math.log(b + w)) < t);
198
199            // Step 5:
200            return computeX(w);
201        }
202
203        @Override
204        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
205            return new ChengBBBetaSampler(rng, this);
206        }
207    }
208
209    /**
210     * Computes one sample using Cheng's BC algorithm, when at least one of beta distribution
211     * {@code alpha} or {@code beta} shape parameters is smaller than 1.
212     */
213    private static class ChengBCBetaSampler extends BaseChengBetaSampler {
214        /** 1/2. */
215        private static final double ONE_HALF = 1d / 2;
216        /** 1/4. */
217        private static final double ONE_QUARTER = 1d / 4;
218
219        /** The algorithm beta factor. This is not the beta distribution beta shape parameter. */
220        private final double beta;
221        /** The algorithm delta factor. */
222        private final double delta;
223        /** The algorithm k1 factor. */
224        private final double k1;
225        /** The algorithm k2 factor. */
226        private final double k2;
227
228        /**
229         * @param rng Generator of uniformly distributed random numbers.
230         * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
231         * @param a max(alpha, beta) shape parameter.
232         * @param b min(alpha, beta) shape parameter.
233         */
234        ChengBCBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
235            super(rng, aIsAlphaShape, a, b);
236            // Compute algorithm factors.
237            beta = 1 / b;
238            delta = 1 + a - b;
239            // The following factors are divided by 4:
240            // k1 = delta * (1 + 3b) / (18a/b - 14)
241            // k2 = 1 + (2 + 1/delta) * b
242            k1 = delta * (1.0 / 72.0 + 3.0 / 72.0 * b) / (a * beta - 7.0 / 9.0);
243            k2 = 0.25 + (0.5 + 0.25 / delta) * b;
244        }
245
246        /**
247         * @param rng Generator of uniformly distributed random numbers.
248         * @param source Source to copy.
249         */
250        private ChengBCBetaSampler(UniformRandomProvider rng,
251                                   ChengBCBetaSampler source) {
252            super(rng, source);
253            beta = source.beta;
254            delta = source.delta;
255            k1 = source.k1;
256            k2 = source.k2;
257        }
258
259        @Override
260        public double sample() {
261            double w;
262            while (true) {
263                // Step 1:
264                final double u1 = rng.nextDouble();
265                final double u2 = rng.nextDouble();
266                // Compute Y and Z
267                final double y = u1 * u2;
268                final double z = u1 * y;
269                if (u1 < ONE_HALF) {
270                    // Step 2:
271                    if (ONE_QUARTER * u2 + z - y >= k1) {
272                        continue;
273                    }
274                } else {
275                    // Step 3:
276                    if (z <= ONE_QUARTER) {
277                        final double v = beta * (Math.log(u1) - Math.log1p(-u1));
278                        w = a * Math.exp(v);
279                        break;
280                    }
281
282                    // Step 4:
283                    if (z >= k2) {
284                        continue;
285                    }
286                }
287
288                // Step 5:
289                final double v = beta * (Math.log(u1) - Math.log1p(-u1));
290                w = a * Math.exp(v);
291                if (alpha * (logAlpha - Math.log(b + w) + v) - LN_4 >= Math.log(z)) {
292                    break;
293                }
294            }
295
296            // Step 6:
297            return computeX(w);
298        }
299
300        @Override
301        public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
302            return new ChengBCBetaSampler(rng, this);
303        }
304    }
305
306    /**
307     * Creates a sampler instance.
308     *
309     * @param rng Generator of uniformly distributed random numbers.
310     * @param alpha Distribution first shape parameter.
311     * @param beta Distribution second shape parameter.
312     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}
313     */
314    public ChengBetaSampler(UniformRandomProvider rng,
315                            double alpha,
316                            double beta) {
317        super(null);
318        delegate = of(rng, alpha, beta);
319    }
320
321    /** {@inheritDoc} */
322    @Override
323    public double sample() {
324        return delegate.sample();
325    }
326
327    /** {@inheritDoc} */
328    @Override
329    public String toString() {
330        return delegate.toString();
331    }
332
333    /**
334     * {@inheritDoc}
335     *
336     * @since 1.3
337     */
338    @Override
339    public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
340        return delegate.withUniformRandomProvider(rng);
341    }
342
343    /**
344     * Creates a new beta distribution sampler.
345     *
346     * @param rng Generator of uniformly distributed random numbers.
347     * @param alpha Distribution first shape parameter.
348     * @param beta Distribution second shape parameter.
349     * @return the sampler
350     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}
351     * @since 1.3
352     */
353    public static SharedStateContinuousSampler of(UniformRandomProvider rng,
354                                                  double alpha,
355                                                  double beta) {
356        if (alpha <= 0) {
357            throw new IllegalArgumentException("alpha is not strictly positive: " + alpha);
358        }
359        if (beta <= 0) {
360            throw new IllegalArgumentException("beta is not strictly positive: " + beta);
361        }
362
363        // Choose the algorithm.
364        final double a = Math.min(alpha, beta);
365        final double b = Math.max(alpha, beta);
366        final boolean aIsAlphaShape = a == alpha;
367
368        return a > 1 ?
369            // BB algorithm
370            new ChengBBBetaSampler(rng, aIsAlphaShape, a, b) :
371            // The BC algorithm is deliberately invoked with reversed parameters
372            // as the argument order is: max(alpha,beta), min(alpha,beta).
373            // Also invert the 'a is alpha' flag.
374            new ChengBCBetaSampler(rng, !aIsAlphaShape, b, a);
375    }
376}