random.py gives wrong results (+ a solution)

[Janne Sinkkonen]

At least in Python 2.0 and earlier, the samples returned by the
function betavariate() of random.py are not from a beta distribution
although the function name misleadingly suggests so.

The following would give beta-distributed samples:

def betavariate(alpha, beta):
     y = gammavariate(alpha,1)
     if y==0: return 0.0
     else: return  y/(y+gammavariate(beta,1))

This is from matlab. A comment in the original matlab code refers to
Devroye, L. (1986) Non-Uniform Random Variate Generation, theorem 4.1A
(p. 430). Another reference would be Gelman, A. et al. (1995) Bayesian
data analysis, p. 481, which I have checked and found to agree with
the code above.

-- 
Janne