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 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 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}