[Python-bugs-list] [ python-Bugs-527139 ] random.gammavariate hosed
noreply@sourceforge.net
noreply@sourceforge.net
Mon, 13 May 2002 15:49:16 -0700
Bugs item #527139, was opened at 2002-03-07 15: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: Remind
Priority: 5
Submitted By: Tim Peters (tim_one)
Assigned to: Raymond Hettinger (rhettinger)
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: Raymond Hettinger (rhettinger)
Date: 2002-05-13 17:49
Message:
Logged In: YES
user_id=80475
Committed documentation update: librandom.tex 1.26 and
1.25.18.1.
Leaving the bug open until a random.py gets patched to
expand its domain. Attaching test code to validate the
patch.
----------------------------------------------------------------------
Comment By: Tim Peters (tim_one)
Date: 2002-05-13 17:33
Message:
Logged In: YES
user_id=31435
I can't foresee ever making time to look at this, so
reassigned it to you. Have fun <wink>.
----------------------------------------------------------------------
Comment By: Raymond Hettinger (rhettinger)
Date: 2002-05-13 17:26
Message:
Logged In: YES
user_id=80475
Tim, is okay to commit the attached documentation patch and
backport it to release22-maint?
Currently, the documentation mis-advertises the method name
and domain. It intended to say alpha>0; however, as the
poster pointed out, the implementation fails for alpha<0.5.
For an immediate action, I think we ought to bring the docs
into alignment with the implementation's actual name and
domain. See attached doc patch.
For a follow-up action, the implementation should be
modified to handle the full domain, alpha>0. This will
take some research (meaning, I don't immediately see how to
fix it <wink>).
Note 1: I tested the current implementation and found that
it does have the expected mean and variance for alpha in
[.5, .75, 1, 2, 4, 8] and beta in [.1, .2, .4, .8, 1, 2, 4,
8]. So, it looks like the code does function well when
alpha >= 0.5.
Note 2: My CRC handbook defines the gammavariate function
with alpha offet by one, alpha > -1. My college
engineering statistics book defines it with alpha>0. The
executive summary is that the alpha offset is non-
standard. This implementation uses the alpha>0 version.
----------------------------------------------------------------------
Comment By: John Machin (sjmachin)
Date: 2002-03-09 21: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