@@ -844,8 +844,8 @@ def binomialvariate(self, n=1, p=0.5):
844844 # BTRS: Transformed rejection with squeeze method by Wolfgang Hörmann
845845 # https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.47.8407&rep=rep1&type=pdf
846846 assert n * p >= 10.0 and p <= 0.5
847- setup_complete = False
848847
848+ setup_complete = False
849849 spq = _sqrt (n * p * (1.0 - p )) # Standard deviation of the distribution
850850 b = 1.15 + 2.53 * spq
851851 a = - 0.0873 + 0.0248 * b + 0.01 * p
@@ -860,22 +860,23 @@ def binomialvariate(self, n=1, p=0.5):
860860 k = _floor ((2.0 * a / us + b ) * u + c )
861861 if k < 0 or k > n :
862862 continue
863+ v = random ()
863864
864865 # The early-out "squeeze" test substantially reduces
865866 # the number of acceptance condition evaluations.
866- v = random ()
867867 if us >= 0.07 and v <= vr :
868868 return k
869869
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.
873870 if not setup_complete :
874871 alpha = (2.83 + 5.1 / b ) * spq
875872 lpq = _log (p / (1.0 - p ))
876873 m = _floor ((n + 1 ) * p ) # Mode of the distribution
877874 h = _lgamma (m + 1 ) + _lgamma (n - m + 1 )
878875 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.
879880 v *= alpha / (a / (us * us ) + b )
880881 if _log (v ) <= h - _lgamma (k + 1 ) - _lgamma (n - k + 1 ) + (k - m ) * lpq :
881882 return k
0 commit comments