@@ -844,8 +844,8 @@ def binomialvariate(self, n=1, p=0.5):
844
844
# BTRS: Transformed rejection with squeeze method by Wolfgang Hörmann
845
845
# https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.47.8407&rep=rep1&type=pdf
846
846
assert n * p >= 10.0 and p <= 0.5
847
- setup_complete = False
848
847
848
+ setup_complete = False
849
849
spq = _sqrt (n * p * (1.0 - p )) # Standard deviation of the distribution
850
850
b = 1.15 + 2.53 * spq
851
851
a = - 0.0873 + 0.0248 * b + 0.01 * p
@@ -860,22 +860,23 @@ def binomialvariate(self, n=1, p=0.5):
860
860
k = _floor ((2.0 * a / us + b ) * u + c )
861
861
if k < 0 or k > n :
862
862
continue
863
+ v = random ()
863
864
864
865
# The early-out "squeeze" test substantially reduces
865
866
# the number of acceptance condition evaluations.
866
- v = random ()
867
867
if us >= 0.07 and v <= vr :
868
868
return k
869
869
870
- # Acceptance-rejection test.
871
- # Note, the original paper erroneously omits the call to log(v)
872
- # when comparing to the log of the rescaled binomial distribution.
873
870
if not setup_complete :
874
871
alpha = (2.83 + 5.1 / b ) * spq
875
872
lpq = _log (p / (1.0 - p ))
876
873
m = _floor ((n + 1 ) * p ) # Mode of the distribution
877
874
h = _lgamma (m + 1 ) + _lgamma (n - m + 1 )
878
875
setup_complete = True # Only needs to be done once
876
+
877
+ # Acceptance-rejection test.
878
+ # Note, the original paper erroneously omits the call to log(v)
879
+ # when comparing to the log of the rescaled binomial distribution.
879
880
v *= alpha / (a / (us * us ) + b )
880
881
if _log (v ) <= h - _lgamma (k + 1 ) - _lgamma (n - k + 1 ) + (k - m ) * lpq :
881
882
return k
0 commit comments