classification
Title: Random.paretovariate sometimes raises ZeroDivisionError for small alpha
Type: behavior Stage: resolved
Components: Library (Lib) Versions: Python 3.10, Python 3.9, Python 3.8
process
Status: closed Resolution: works for me
Dependencies: Superseder:
Assigned To: Nosy List: David MacIver, mark.dickinson, rhettinger, serhiy.storchaka, tim.peters
Priority: normal Keywords: patch

Created on 2020-07-28 17:59 by David MacIver, last changed 2020-08-02 00:40 by rhettinger. This issue is now closed.

Pull Requests
URL Status Linked Edit
PR 21695 merged rhettinger, 2020-07-31 09:02
Messages (12)
msg374516 - (view) Author: David MacIver (David MacIver) * Date: 2020-07-28 17:59
The following code raises a ZeroDivisionError:


from random import Random

Random(14064741636871487939).paretovariate(0.01)


This raises:

random.py, line 692, in paretovariate
    return 1.0 / u ** (1.0/alpha)
ZeroDivisionError: float division by zero


That specific stack trace is from 3.8.5 but I've tried this on 3.6 and 3.7 as well.

It's a little hard to tell what the intended correct range of parameter values is for paretovariate, and this may just be lack of validation for the alpha < 1 case - perhaps that was never intended to be supported?

Based on some very informal inspection, what seems to happen is that the probability of getting a ZeroDivisionError approaches one as you approach zero. They rarely occur at this level of alpha (0.01), but for alpha=0.001 they seem to occur just under half the time, and for alpha=0.0001 they occur more than 90% of the time.

(For the interested, this bug was found by Hypothesis as part of its own testing of our integration with the Random API)
msg374517 - (view) Author: David MacIver (David MacIver) * Date: 2020-07-28 18:03
I guess on actual inspection of the code (which I should have done before, sorry) it's obvious why this happens: u ** (1.0 / alpha) will round to 0.0 for small values of u, and the smaller alpha is the higher the probability of that happening, in that it happens with probability c ** alpha for some c.
msg374537 - (view) Author: Tim Peters (tim.peters) * (Python committer) Date: 2020-07-28 23:35
I'm inclined to ignore this. No actual user has complained about this, and I doubt any ever will:  it's got to be rare as hen's teeth to use a parameter outside of, say, [0.1, 10.0], in real life. The division error can't happen for those. For "small" alpha, the distribution simply does contain values too large to represent as binary doubles. So the only other defensible thing to do be done in this case is to return math.inf (or raise OverflowError) - but why slow it down for something no actual user cares about?
msg374538 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2020-07-28 23:36
The point of the "u = 1.0 - self.random()" line was to prevent the case where *u* was exactly equal to zero; however, as you noted, raising a very small number to a small can round down to zero.

We could wrap the calculation in a try/except to catch the ZeroDivisionError but that is probably insufficient.  The codomain can easily spew beyond the range of a C double:

   >>> alpha = 0.01
   >>> u = 0.005
   >>> 1.0 / u ** (1.0 / alpha)
   1.2676506002282268e+230             <== Huge!

We could clamp the values as is done in gammavariate().  Or we can raise a clean Overflow or Underflow exception but that would be unpleasant for users if it were to arise unexpectedly.
msg374546 - (view) Author: Tim Peters (tim.peters) * (Python committer) Date: 2020-07-28 23:54
BTW, if we have to "do something", how about changing

return 1.0 / u ** (1.0/alpha)

to the mathematically equivalent

return (1.0 / u) ** (1.0/alpha)

? Not sure about Linux-y boxes, but on Windows that would raise OverflowError instead of ZeroDivisionError. Which makes more sense, because the Pareto variable we're _trying_ to produce is too large to represent.  If it returns math.inf on some box instead, same thing.

Or, also faster, and suffering fewer rounding errors,

return  u ** (-1.0 / alpha)
msg374548 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2020-07-29 00:05
I don't feel a need to change anything, my post hit a few seconds after yours :-) 

That said, "return u ** (-1.0 / alpha)" would be a slight, but nice improvement.

It might also be reasonable to document a safe range of alpha values.
msg374564 - (view) Author: David MacIver (David MacIver) * Date: 2020-07-29 08:33
I should say, I don't actually care about this bug at all: I only ran into it because of trying to recreate the random API and testing that my recreation worked sensibly. I thought I'd do the good citizen thing and report it, but I'm totally fine with any resolution and given that that code has been unchanged for at least 18 years I doubt anyone else has ever run into this problem.

I do think it would be beneficial to document the range of reasonable alpha on paretovariate - the lack of such documentation is the only reason I hit this at all - but I think leaving the behaviour as is would be fine.
msg374565 - (view) Author: Serhiy Storchaka (serhiy.storchaka) * (Python committer) Date: 2020-07-29 09:01
+1 for "u ** (-1.0 / alpha)"!
msg374659 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2020-08-01 08:18
New changeset 5c3270939c09e4c8e80fd26449b718a998701912 by Raymond Hettinger in branch 'master':
bpo-41421: Algebraic simplification for random.paretovariate() (GH-21695)
https://github.com/python/cpython/commit/5c3270939c09e4c8e80fd26449b718a998701912
msg374666 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2020-08-01 19:36
How about this wording?

    If alpha is smaller than 0.1, an occasional OverflowError can
    arise when the variate exceeds the range of Python float.

I like using 0.1 because it's easy and gives us some wiggle room.  The actual cutoff is lower:

    >>> alpha = 0.05079
    >>> sys.float_info.epsilon ** (-1.0 / alpha)

    >>> alpha = 0.05078
    >>> sys.float_info.epsilon ** (-1.0 / alpha)
    Traceback (most recent call last):
       File "<pyshell#48>", line 1, in <module>
         sys.float_info.epsilon ** (-1.0 / alpha)
    1.590779741838475e+308
msg374672 - (view) Author: Tim Peters (tim.peters) * (Python committer) Date: 2020-08-01 21:07
That text is fine, if you feel something needs to be said at all. I really don't. A Pareto distribution of this kind with parameter <= 1.0 has infinite expected value - VERY long tail. Anyone who knows what they're doing already knows that. The reason the implementation can't "blow up" for parameters >= (roughly) 0.1 isn't that IEEE doubles have such a large dynamic range but rather that they can't represent a number < 1.0 larger than 1 - 2**-53 (so u = 1 - random.random() is always at least 2**-53). The actual distribution has infinite expected value nonetheless, but the implementation is incapable of generating any of its very large values (which, while very rare, are also very large).
msg374674 - (view) Author: Raymond Hettinger (rhettinger) * (Python committer) Date: 2020-08-02 00:40
Okay, marking this as closed.
History
Date User Action Args
2020-08-02 00:40:14rhettingersetstatus: open -> closed
resolution: works for me
messages: + msg374674

stage: patch review -> resolved
2020-08-01 21:07:18tim.peterssetmessages: + msg374672
2020-08-01 19:36:19rhettingersetmessages: + msg374666
2020-08-01 08:18:50rhettingersetmessages: + msg374659
2020-07-31 09:02:10rhettingersetkeywords: + patch
stage: patch review
pull_requests: + pull_request20839
2020-07-29 09:01:57serhiy.storchakasetmessages: + msg374565
2020-07-29 08:34:00David MacIversetmessages: + msg374564
2020-07-29 00:05:17rhettingersetmessages: + msg374548
2020-07-28 23:54:41tim.peterssetmessages: + msg374546
2020-07-28 23:36:31rhettingersettype: behavior
versions: + Python 3.9, Python 3.10, - Python 3.6, Python 3.7
2020-07-28 23:36:08rhettingersetmessages: + msg374538
2020-07-28 23:35:48tim.peterssetmessages: + msg374537
2020-07-28 23:24:26rhettingersetnosy: + tim.peters
2020-07-28 19:19:11serhiy.storchakasetnosy: + serhiy.storchaka
2020-07-28 18:37:42xtreaksetnosy: + rhettinger, mark.dickinson
2020-07-28 18:03:11David MacIversetmessages: + msg374517
2020-07-28 17:59:09David MacIvercreate