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}