32
32
circular uniform
33
33
von Mises
34
34
35
+ discrete distributions
36
+ ----------------------
37
+ binomial
38
+
39
+
35
40
General notes on the underlying Mersenne Twister core generator:
36
41
37
42
* The period is 2**19937-1.
49
54
from math import log as _log , exp as _exp , pi as _pi , e as _e , ceil as _ceil
50
55
from math import sqrt as _sqrt , acos as _acos , cos as _cos , sin as _sin
51
56
from math import tau as TWOPI , floor as _floor , isfinite as _isfinite
57
+ from math import lgamma as _lgamma , fabs as _fabs , log2 as _log2
52
58
try :
53
59
from os import urandom as _urandom
54
60
except ImportError :
@@ -72,7 +78,7 @@ def _urandom(*args, **kwargs):
72
78
73
79
try :
74
80
# hashlib is pretty heavy to load, try lean internal module first
75
- from _sha512 import sha512 as _sha512
81
+ from _sha2 import sha512 as _sha512
76
82
except ImportError :
77
83
# fallback to official implementation
78
84
from hashlib import sha512 as _sha512
@@ -81,6 +87,7 @@ def _urandom(*args, **kwargs):
81
87
"Random" ,
82
88
"SystemRandom" ,
83
89
"betavariate" ,
90
+ "binomialvariate" ,
84
91
"choice" ,
85
92
"choices" ,
86
93
"expovariate" ,
@@ -249,7 +256,7 @@ def _randbelow_with_getrandbits(self, n):
249
256
"Return a random int in the range [0,n). Defined for n > 0."
250
257
251
258
getrandbits = self .getrandbits
252
- k = n .bit_length () # don't use (n-1) here because n can be 1
259
+ k = n .bit_length ()
253
260
r = getrandbits (k ) # 0 <= r < 2**k
254
261
while r >= n :
255
262
r = getrandbits (k )
@@ -304,58 +311,25 @@ def randrange(self, start, stop=None, step=_ONE):
304
311
305
312
# This code is a bit messy to make it fast for the
306
313
# common case while still doing adequate error checking.
307
- try :
308
- istart = _index (start )
309
- except TypeError :
310
- istart = int (start )
311
- if istart != start :
312
- _warn ('randrange() will raise TypeError in the future' ,
313
- DeprecationWarning , 2 )
314
- raise ValueError ("non-integer arg 1 for randrange()" )
315
- _warn ('non-integer arguments to randrange() have been deprecated '
316
- 'since Python 3.10 and will be removed in a subsequent '
317
- 'version' ,
318
- DeprecationWarning , 2 )
314
+ istart = _index (start )
319
315
if stop is None :
320
316
# We don't check for "step != 1" because it hasn't been
321
317
# type checked and converted to an integer yet.
322
318
if step is not _ONE :
323
- raise TypeError (' Missing a non-None stop argument' )
319
+ raise TypeError (" Missing a non-None stop argument" )
324
320
if istart > 0 :
325
321
return self ._randbelow (istart )
326
322
raise ValueError ("empty range for randrange()" )
327
323
328
- # stop argument supplied.
329
- try :
330
- istop = _index (stop )
331
- except TypeError :
332
- istop = int (stop )
333
- if istop != stop :
334
- _warn ('randrange() will raise TypeError in the future' ,
335
- DeprecationWarning , 2 )
336
- raise ValueError ("non-integer stop for randrange()" )
337
- _warn ('non-integer arguments to randrange() have been deprecated '
338
- 'since Python 3.10 and will be removed in a subsequent '
339
- 'version' ,
340
- DeprecationWarning , 2 )
324
+ # Stop argument supplied.
325
+ istop = _index (stop )
341
326
width = istop - istart
342
- try :
343
- istep = _index (step )
344
- except TypeError :
345
- istep = int (step )
346
- if istep != step :
347
- _warn ('randrange() will raise TypeError in the future' ,
348
- DeprecationWarning , 2 )
349
- raise ValueError ("non-integer step for randrange()" )
350
- _warn ('non-integer arguments to randrange() have been deprecated '
351
- 'since Python 3.10 and will be removed in a subsequent '
352
- 'version' ,
353
- DeprecationWarning , 2 )
327
+ istep = _index (step )
354
328
# Fast path.
355
329
if istep == 1 :
356
330
if width > 0 :
357
331
return istart + self ._randbelow (width )
358
- raise ValueError ("empty range for randrange() (%d, %d, %d)" % ( istart , istop , width ) )
332
+ raise ValueError (f "empty range in randrange({ start } , { stop } )" )
359
333
360
334
# Non-unit step argument supplied.
361
335
if istep > 0 :
@@ -365,7 +339,7 @@ def randrange(self, start, stop=None, step=_ONE):
365
339
else :
366
340
raise ValueError ("zero step for randrange()" )
367
341
if n <= 0 :
368
- raise ValueError ("empty range for randrange()" )
342
+ raise ValueError (f "empty range in randrange({ start } , { stop } , { step } )" )
369
343
return istart + istep * self ._randbelow (n )
370
344
371
345
def randint (self , a , b ):
@@ -531,7 +505,14 @@ def choices(self, population, weights=None, *, cum_weights=None, k=1):
531
505
## -------------------- real-valued distributions -------------------
532
506
533
507
def uniform (self , a , b ):
534
- "Get a random number in the range [a, b) or [a, b] depending on rounding."
508
+ """Get a random number in the range [a, b) or [a, b] depending on rounding.
509
+
510
+ The mean (expected value) and variance of the random variable are:
511
+
512
+ E[X] = (a + b) / 2
513
+ Var[X] = (b - a) ** 2 / 12
514
+
515
+ """
535
516
return a + (b - a ) * self .random ()
536
517
537
518
def triangular (self , low = 0.0 , high = 1.0 , mode = None ):
@@ -542,6 +523,11 @@ def triangular(self, low=0.0, high=1.0, mode=None):
542
523
543
524
http://en.wikipedia.org/wiki/Triangular_distribution
544
525
526
+ The mean (expected value) and variance of the random variable are:
527
+
528
+ E[X] = (low + high + mode) / 3
529
+ Var[X] = (low**2 + high**2 + mode**2 - low*high - low*mode - high*mode) / 18
530
+
545
531
"""
546
532
u = self .random ()
547
533
try :
@@ -623,7 +609,7 @@ def lognormvariate(self, mu, sigma):
623
609
"""
624
610
return _exp (self .normalvariate (mu , sigma ))
625
611
626
- def expovariate (self , lambd ):
612
+ def expovariate (self , lambd = 1.0 ):
627
613
"""Exponential distribution.
628
614
629
615
lambd is 1.0 divided by the desired mean. It should be
@@ -632,12 +618,15 @@ def expovariate(self, lambd):
632
618
positive infinity if lambd is positive, and from negative
633
619
infinity to 0 if lambd is negative.
634
620
635
- """
636
- # lambd: rate lambd = 1/mean
637
- # ('lambda' is a Python reserved word)
621
+ The mean (expected value) and variance of the random variable are:
638
622
623
+ E[X] = 1 / lambd
624
+ Var[X] = 1 / lambd ** 2
625
+
626
+ """
639
627
# we use 1-random() instead of random() to preclude the
640
628
# possibility of taking the log of zero.
629
+
641
630
return - _log (1.0 - self .random ()) / lambd
642
631
643
632
def vonmisesvariate (self , mu , kappa ):
@@ -693,8 +682,12 @@ def gammavariate(self, alpha, beta):
693
682
pdf(x) = --------------------------------------
694
683
math.gamma(alpha) * beta ** alpha
695
684
685
+ The mean (expected value) and variance of the random variable are:
686
+
687
+ E[X] = alpha * beta
688
+ Var[X] = alpha * beta ** 2
689
+
696
690
"""
697
- # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
698
691
699
692
# Warning: a few older sources define the gamma distribution in terms
700
693
# of alpha > -1.0
@@ -753,6 +746,11 @@ def betavariate(self, alpha, beta):
753
746
Conditions on the parameters are alpha > 0 and beta > 0.
754
747
Returned values range between 0 and 1.
755
748
749
+ The mean (expected value) and variance of the random variable are:
750
+
751
+ E[X] = alpha / (alpha + beta)
752
+ Var[X] = alpha * beta / ((alpha + beta)**2 * (alpha + beta + 1))
753
+
756
754
"""
757
755
## See
758
756
## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html
@@ -793,6 +791,97 @@ def weibullvariate(self, alpha, beta):
793
791
return alpha * (- _log (u )) ** (1.0 / beta )
794
792
795
793
794
+ ## -------------------- discrete distributions ---------------------
795
+
796
+ def binomialvariate (self , n = 1 , p = 0.5 ):
797
+ """Binomial random variable.
798
+
799
+ Gives the number of successes for *n* independent trials
800
+ with the probability of success in each trial being *p*:
801
+
802
+ sum(random() < p for i in range(n))
803
+
804
+ Returns an integer in the range: 0 <= X <= n
805
+
806
+ The mean (expected value) and variance of the random variable are:
807
+
808
+ E[X] = n * p
809
+ Var[x] = n * p * (1 - p)
810
+
811
+ """
812
+ # Error check inputs and handle edge cases
813
+ if n < 0 :
814
+ raise ValueError ("n must be non-negative" )
815
+ if p <= 0.0 or p >= 1.0 :
816
+ if p == 0.0 :
817
+ return 0
818
+ if p == 1.0 :
819
+ return n
820
+ raise ValueError ("p must be in the range 0.0 <= p <= 1.0" )
821
+
822
+ random = self .random
823
+
824
+ # Fast path for a common case
825
+ if n == 1 :
826
+ return _index (random () < p )
827
+
828
+ # Exploit symmetry to establish: p <= 0.5
829
+ if p > 0.5 :
830
+ return n - self .binomialvariate (n , 1.0 - p )
831
+
832
+ if n * p < 10.0 :
833
+ # BG: Geometric method by Devroye with running time of O(np).
834
+ # https://dl.acm.org/doi/pdf/10.1145/42372.42381
835
+ x = y = 0
836
+ c = _log2 (1.0 - p )
837
+ if not c :
838
+ return x
839
+ while True :
840
+ y += _floor (_log2 (random ()) / c ) + 1
841
+ if y > n :
842
+ return x
843
+ x += 1
844
+
845
+ # BTRS: Transformed rejection with squeeze method by Wolfgang Hörmann
846
+ # https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.47.8407&rep=rep1&type=pdf
847
+ assert n * p >= 10.0 and p <= 0.5
848
+ setup_complete = False
849
+
850
+ spq = _sqrt (n * p * (1.0 - p )) # Standard deviation of the distribution
851
+ b = 1.15 + 2.53 * spq
852
+ a = - 0.0873 + 0.0248 * b + 0.01 * p
853
+ c = n * p + 0.5
854
+ vr = 0.92 - 4.2 / b
855
+
856
+ while True :
857
+
858
+ u = random ()
859
+ u -= 0.5
860
+ us = 0.5 - _fabs (u )
861
+ k = _floor ((2.0 * a / us + b ) * u + c )
862
+ if k < 0 or k > n :
863
+ continue
864
+
865
+ # The early-out "squeeze" test substantially reduces
866
+ # the number of acceptance condition evaluations.
867
+ v = random ()
868
+ if us >= 0.07 and v <= vr :
869
+ return k
870
+
871
+ # Acceptance-rejection test.
872
+ # Note, the original paper erroneously omits the call to log(v)
873
+ # when comparing to the log of the rescaled binomial distribution.
874
+ if not setup_complete :
875
+ alpha = (2.83 + 5.1 / b ) * spq
876
+ lpq = _log (p / (1.0 - p ))
877
+ m = _floor ((n + 1 ) * p ) # Mode of the distribution
878
+ h = _lgamma (m + 1 ) + _lgamma (n - m + 1 )
879
+ setup_complete = True # Only needs to be done once
880
+ v *= alpha / (a / (us * us ) + b )
881
+ if _log (v ) <= h - _lgamma (k + 1 ) - _lgamma (n - k + 1 ) + (k - m ) * lpq :
882
+ return k
883
+
884
+
796
885
## ------------------------------------------------------------------
797
886
## --------------- Operating System Random Source ------------------
798
887
@@ -859,6 +948,7 @@ def _notimplemented(self, *args, **kwds):
859
948
gammavariate = _inst .gammavariate
860
949
gauss = _inst .gauss
861
950
betavariate = _inst .betavariate
951
+ binomialvariate = _inst .binomialvariate
862
952
paretovariate = _inst .paretovariate
863
953
weibullvariate = _inst .weibullvariate
864
954
getstate = _inst .getstate
@@ -883,15 +973,17 @@ def _test_generator(n, func, args):
883
973
low = min (data )
884
974
high = max (data )
885
975
886
- print (f'{ t1 - t0 :.3f} sec, { n } times { func .__name__ } ' )
976
+ print (f'{ t1 - t0 :.3f} sec, { n } times { func .__name__ } { args !r } ' )
887
977
print ('avg %g, stddev %g, min %g, max %g\n ' % (xbar , sigma , low , high ))
888
978
889
979
890
- def _test (N = 2000 ):
980
+ def _test (N = 10_000 ):
891
981
_test_generator (N , random , ())
892
982
_test_generator (N , normalvariate , (0.0 , 1.0 ))
893
983
_test_generator (N , lognormvariate , (0.0 , 1.0 ))
894
984
_test_generator (N , vonmisesvariate , (0.0 , 1.0 ))
985
+ _test_generator (N , binomialvariate , (15 , 0.60 ))
986
+ _test_generator (N , binomialvariate , (100 , 0.75 ))
895
987
_test_generator (N , gammavariate , (0.01 , 1.0 ))
896
988
_test_generator (N , gammavariate , (0.1 , 1.0 ))
897
989
_test_generator (N , gammavariate , (0.1 , 2.0 ))
0 commit comments