I'm pretty sure Tim's seen this already, but just in case... ----- Original Message ----- From: "Ivan Frohne" <frohne@gci.net> Newsgroups: comp.lang.python Sent: Thursday, January 25, 2001 5:20 PM Subject: Re: random.py gives wrong results (+ a solution) > > "Janne Sinkkonen" <janne@oops.nnets.fi> wrote in message > news:m3u26oy1rw.fsf@kinos.nnets.fi... > > > > 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. > > > 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. > > 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. > > The probability density function of the gamma distribution with (positive) > parameters > A and B is usually written > > g(x; A, B) = (x**(A-1) * exp(x/B)) / (Gamma(A) * B**A), where x, A, and > B > 0. > > Here Gamma(A) is the gamma function -- for A a positive integer, Gamma(A) is > the > factorial of A - 1, Gamma(A) = (A-1)!. In fact, this is the definition used > by the authors of random.py in defining gammavariate(alpha, beta), the gamma > distribution random number generator. > > Now it happens that a gamma-distributed random variable with parameters A = > 1 and > B has the (much simpler) exponential distribution with density function > > g(x; 1, B) = exp(-x/B) / B. > > Keep that in mind. > > The reference "Discrete Event Simulation in ," by Kevin Watkins > (McGraw-Hill, 1993) > was consulted by the random.py authors. But this reference defines the > gamma probability distribution a little differently, as > > g1(x; A, B) = (B**A * x**(A-1) * exp(B*x)) / Gamma(A), where x, A, B > > 0. > > (See p. 85). On page 87, Watkins states (incorrectly) that if grv(A, B) is > a function which > returns a gamma random variable with parameters A and B (using his > definition on p. 85), > then the function > > brv(A, B) = grv(1, 1/B) / ( grv(1, 1/B) + grv(1, A) ) [ not > true!] > > will return a random variable which has the beta distribution with > parameters A and B. > > Believing Watkins to be correct, the random.py authors remembered that a > gamma > random variable with parameter A = 1 is just an exponential random variable > and > further simplified their beta generator to > > brv(A, B) = erv(1/B) / (erv(1/B) + erv(A)), where erv(K) is a random > variable > > having the exponential distribution with > > parameter K. > > 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). > 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. > > --Ivan Frohne > > > > > >>> > > > >
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4