[Python-bugs-list] [ python-Bugs-527139 ] random.gammavariate hosed

[noreply@sourceforge.net]

Bugs item #527139, was opened at 2002-03-08 07:44
You can respond by visiting: 
http://sourceforge.net/tracker/?func=detail&atid=105470&aid=527139&group_id=5470

Category: Python Library
Group: None
Status: Open
Resolution: None
Priority: 5
Submitted By: Tim Peters (tim_one)
Assigned to: Nobody/Anonymous (nobody)
Summary: random.gammavariate hosed

Initial Comment:
Report from c.l.py; something is certainly wrong here, 
but will take some digging to figure out whether docs 
or code:

"""
From: Weet Vanniks
Sent: Thursday, March 07, 2002 11:14 AM
To: python-list@python.org
Subject: Bug in the standard module random ?


Hi,

The gammavariate function of the standard module is 
documented as taking two parameters, the first one 
alpha is required to be > -1 and the second beta is 
required to be >0.  However, examining the 
implementation, it seems that the requirement for 
alpha is to be >0.

In spite of this, I still have a problem since I 
called the gammavariate function with a parameter 
alpha equal to 0.2 and it fails logically on the 
following line:

ainv=_sqrt(2.0 * alpha - 1.0)

Apparently, the implementation requires alpha to be > 
0.5.  Am i missing something or is it a bug?
"""

----------------------------------------------------------------------

Comment By: John Machin (sjmachin)
Date: 2002-03-10 13:09

Message:
Logged In: YES 
user_id=480138

Edited & revised version of my posts on c.l.py:

=== Point 1 =================
There is a doc problem. Some definitions of the gamma 
distribution define it (a) such that alpha > -1 and has the 
exponential distribution as a special case when alpha == 0; 
others take the tack (b) that the lower bound is zero and 
the exponential distribution is when alpha == 1.

Option (b) seems to be the most popular recently.
Example of (a): my very old probability textbook.
Example of (b):
http://www.nag.com/numeric/cl/manual/pdf/G05/g05ffc_cl05.pdf

Here we have the doc taking option (a) and the 
implementation taking
option (b).

N.B. the following is VERY relevant:
http://mail.python.org/pipermail/python-dev/2001-
January/012181.html

=== Point 2 =============
gammavariate() calls stdgamma() with arguments that include 
ainv, bbb, and ccc. stdgamma() uses those arguments only if 
alpha > 1.0. As the OP found, ainv is off the planet for 
alpha < 0.5.

stdgamma() is called (inside random.py) only once, from 
gammavariate()

The module exposes both gammavariate() and stdgamma() --- 
stdgamma() with those three extra args seems useless to me. 
It should not be an exposed function IMO.

=== Point 3 ========

The documentation says there is a gamma() but doesn't 
mention gammavariate() and stdgamma(). The exposed function 
should be called "gammavariate" (a) for consistency with 
other functions (b) to avoid the perceived need to explain 
that this "gamma" is not the gamma function!!!

[post-posting points]

=== Point 4 ==== 
Serendipitous coincidence: http://www.jstatsoft.org/v03/i03/

=== Point 5 ====

The algorithm that Python uses starts falling apart from 
underflow for small alpha. So does the algorithm in GSL 
(GNU Scientific Library). So does the algorithm in the 
Monty Python paper. The problem in all of them is that they 
do pow(random(), 1.0 / alpha) or exp(log(random) / alpha) 
or suchlike. For alpha == 0.00001 (say), the result is 
vanishingly small and most generated variates are == 0.0.  
Note that the Monty Python paper says "the rare but 
difficult case alpha < 1". Given the bug we currently have 
with alpha < 0.5, it would appear to be even rarer than 
using floats as dict keys :-) -- as we don't purport to be 
NAG or even GSL, I propose we do no more about this than a 
one-line warning in the doc. Perhaps the whole random 
module needs a "caveat emptor" up the front.

=== Point 6 ====
I shall prepare a patch for random.py and upload it real 
soon now.

----------------------------------------------------------------------

You can respond by visiting: 
http://sourceforge.net/tracker/?func=detail&atid=105470&aid=527139&group_id=5470