1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.rng.sampling.distribution;
18
19 import org.apache.commons.rng.UniformRandomProvider;
20
21 /**
22 * Sampling from a <a href="http://en.wikipedia.org/wiki/Beta_distribution">beta distribution</a>.
23 *
24 * <p>Uses Cheng's algorithms for beta distribution sampling:</p>
25 *
26 * <blockquote>
27 * <pre>
28 * R. C. H. Cheng,
29 * "Generating beta variates with nonintegral shape parameters",
30 * Communications of the ACM, 21, 317-322, 1978.
31 * </pre>
32 * </blockquote>
33 *
34 * <p>Sampling uses {@link UniformRandomProvider#nextDouble()}.</p>
35 *
36 * @since 1.0
37 */
38 public class ChengBetaSampler
39 extends SamplerBase
40 implements SharedStateContinuousSampler {
41 /** Natural logarithm of 4. */
42 private static final double LN_4 = Math.log(4.0);
43
44 /** The appropriate beta sampler for the parameters. */
45 private final SharedStateContinuousSampler delegate;
46
47 /**
48 * Base class to implement Cheng's algorithms for the beta distribution.
49 */
50 private abstract static class BaseChengBetaSampler
51 implements SharedStateContinuousSampler {
52 /**
53 * Flag set to true if {@code a} is the beta distribution alpha shape parameter.
54 * Otherwise {@code a} is the beta distribution beta shape parameter.
55 *
56 * <p>From the original Cheng paper this is equal to the result of: {@code a == a0}.</p>
57 */
58 protected final boolean aIsAlphaShape;
59 /**
60 * First shape parameter.
61 * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
62 */
63 protected final double a;
64 /**
65 * Second shape parameter.
66 * The meaning of this is dependent on the {@code aIsAlphaShape} flag.
67 */
68 protected final double b;
69 /** Underlying source of randomness. */
70 protected final UniformRandomProvider rng;
71 /**
72 * The algorithm alpha factor. This is not the beta distribution alpha shape parameter.
73 * It is the sum of the two shape parameters ({@code a + b}.
74 */
75 protected final double alpha;
76 /** The logarithm of the alpha factor. */
77 protected final double logAlpha;
78
79 /**
80 * @param rng Generator of uniformly distributed random numbers.
81 * @param aIsAlphaShape true if {@code a} is the beta distribution alpha shape parameter.
82 * @param a Distribution first shape parameter.
83 * @param b Distribution second shape parameter.
84 */
85 BaseChengBetaSampler(UniformRandomProvider rng, boolean aIsAlphaShape, double a, double b) {
86 this.rng = rng;
87 this.aIsAlphaShape = aIsAlphaShape;
88 this.a = a;
89 this.b = b;
90 alpha = a + b;
91 logAlpha = Math.log(alpha);
92 }
93
94 /**
95 * @param rng Generator of uniformly distributed random numbers.
96 * @param source Source to copy.
97 */
98 BaseChengBetaSampler(UniformRandomProvider rng,
99 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 final 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 final 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 this(of(rng, alpha, beta));
318 }
319
320 /**
321 * @param delegate Beta sampler.
322 */
323 private ChengBetaSampler(SharedStateContinuousSampler delegate) {
324 super(null);
325 this.delegate = delegate;
326 }
327
328 /** {@inheritDoc} */
329 @Override
330 public double sample() {
331 return delegate.sample();
332 }
333
334 /** {@inheritDoc} */
335 @Override
336 public String toString() {
337 return delegate.toString();
338 }
339
340 /**
341 * {@inheritDoc}
342 *
343 * @since 1.3
344 */
345 @Override
346 public SharedStateContinuousSampler withUniformRandomProvider(UniformRandomProvider rng) {
347 return delegate.withUniformRandomProvider(rng);
348 }
349
350 /**
351 * Creates a new beta distribution sampler.
352 *
353 * @param rng Generator of uniformly distributed random numbers.
354 * @param alpha Distribution first shape parameter.
355 * @param beta Distribution second shape parameter.
356 * @return the sampler
357 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}
358 * @since 1.3
359 */
360 public static SharedStateContinuousSampler of(UniformRandomProvider rng,
361 double alpha,
362 double beta) {
363 InternalUtils.requireStrictlyPositive(alpha, "alpha");
364 InternalUtils.requireStrictlyPositive(beta, "beta");
365
366 // Choose the algorithm.
367 final double a = Math.min(alpha, beta);
368 final double b = Math.max(alpha, beta);
369 final boolean aIsAlphaShape = a == alpha;
370
371 return a > 1 ?
372 // BB algorithm
373 new ChengBBBetaSampler(rng, aIsAlphaShape, a, b) :
374 // The BC algorithm is deliberately invoked with reversed parameters
375 // as the argument order is: max(alpha,beta), min(alpha,beta).
376 // Also invert the 'a is alpha' flag.
377 new ChengBCBetaSampler(rng, !aIsAlphaShape, b, a);
378 }
379 }