If
G_{a} and
G_{b} are independent random variables,
G_{a} having a gamma distribution with shape parameter
a > 0 and
G_{b} having a gamma distribution with shape parameter
b > 0, then the random variable
X = G_{a} / ( G_{a} + G_{b} ) has a beta distribution with shape parameters
a and
b. In particular the probability density for
X,
f(x), is:
f(x) =  0  for x < 0  [ 1 / B(a,b) ] x^{ a1}(1  x)^{ b1}  for 0 ≤ x ≤ 1  0  for x > 1 


where
B(a,b) is the beta function.
 double Beta_Random_Variate( double a, double b)
This function returns X where X is a random variable with the density given above. Note that both a and b are positive.
C source code is available for this routine: