Movatterモバイル変換
[0]ホーム
[Python-Dev] Fw: random.py gives wrong results (+ a solution)
Fredrik Lundhfredrik@effbot.org
Thu, 25 Jan 2001 20:23:50 +0100
I'm pretty sure Tim's seen this already, but justin case...----- Original Message ----- From: "Ivan Frohne" <frohne@gci.net>Newsgroups: comp.lang.pythonSent: Thursday, January 25, 2001 5:20 PMSubject: 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>>>>> >>>>>>>
[8]ページ先頭