random.py gives wrong results (+ a solution)
Tim Peters
tim.one at home.com
Fri Jan 26 01:37:27 EST 2001
[Ivan Frohne]
> I'm convinced that Janne Sinkkonen is right: The beta distribution
> generator in module random.py does not return Beta-distributed random
> numbers. Janne's suggested fix should work just fine.
Which is (for reference):
def betavariate(alpha, beta):
y = gammavariate(alpha,1)
if y==0: return 0.0
else: return y/(y+gammavariate(beta,1))
> Here's my guess on how and why this bug bit -- it won't be of interest
> to most but this subject is so obscure sometimes that there needs to
> be a detailed analysis.
> ...
Which I'll skip, since it's immortalized already in c.l.py and the bug
report:
https://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
> ...
> The corrected equation for a beta random variable, using Watkins'
> definition of the gamma density, is
>
> brv(A, B) = grv(A, 1) / ( grv(A, 1) + grv(1/B, 1) ),
>
> which translates to
>
> brv(A, B) = grv(A, 1) / (grv(A, 1) + grv(B, 1)
>
> using the more common gamma density definition (the one used in random.
> py).
Now that seems plain wrong, although I assume it's just typographical
confusions:
1. You don't really mean to generate two independent instances of grv(A, 1),
right?
2. The parens around the summand denominator have magically disappeared,
changing something of the form X/(X+Y) to (X/X)+Y.
> Many standard statistical references give this equation -- two are
> "Non-Uniform random Variate Generation," by Luc Devroye, Springer-Verlag,
> 1986, p. 432, and "Monte Carlo Concepts, Algorithms and Applications,"
> by George S. Fishman, Springer, 1996, p. 200.
So does Knuth, Vol 3 Ed 3 (in the X/(X+Y) form), although he's not careful
(as was Janne) to avoid division by 0. So I'll check in Janne's code as-is.
Thank you for the analysis! It was an enlightening help.
More information about the Python-list
mailing list