random.py revision 76ca1d428f96284ed58f4523b698ed95c6fdbdb2
1e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum"""Random variable generators. 2e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 3d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters integers 4d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters -------- 5d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform within range 6d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 7d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequences 8d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters --------- 9d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters pick random element 10f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger pick random sample 11d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generate random permutation 12d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 13e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the real line: 14e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum ------------------------------ 15d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters uniform 16e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum normal (Gaussian) 17e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum lognormal 18e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum negative exponential 19e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum gamma 20e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum beta 2140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger pareto 2240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Weibull 23e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 24e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum distributions on the circle (angles 0 to 2pi) 25e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum --------------------------------------------- 26e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum circular uniform 27e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum von Mises 28e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum 2940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond HettingerGeneral notes on the underlying Mersenne Twister core generator: 3040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 3140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The period is 2**19937-1. 3240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* It is one of the most extensively tested generators in existence 3340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* Without a direct way to compute N steps forward, the 3440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger semantics of jumpahead(n) are weakened to simply jump 3540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger to another distant state and rely on the large period 3640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger to avoid overlapping sequences. 3740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger* The random() method is implemented in C, executes in 3840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a single Python step, and is, therefore, threadsafe. 3940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 40e7b146fb3bdca62a0d5ecc06dbf3348e5a4fe757Guido van Rossum""" 41d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 42d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import log as _log, exp as _exp, pi as _pi, e as _e 43d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 449146f27b7799dab231083f194a14c6157b57549fTim Petersfrom math import floor as _floor 45d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 46f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger__all__ = ["Random","seed","random","uniform","randint","choice","sample", 470de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "randrange","shuffle","normalvariate","lognormvariate", 480de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "cunifvariate","expovariate","vonmisesvariate","gammavariate", 490de65807e6bdc5254f5a7e99b2f39adeea6b883bSkip Montanaro "stdgamma","gauss","betavariate","paretovariate","weibullvariate", 5040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "getstate","setstate","jumpahead"] 51ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 52d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 53d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersTWOPI = 2.0*_pi 54d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersLOG4 = _log(4.0) 55d7b5e88e8e40b77813ceb25dc28b87d672538403Tim PetersSG_MAGICCONST = 1.0 + _log(4.5) 5633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 57d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters# Translated by Guido van Rossum from C source provided by 5840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# Adrian Baddeley. Adapted by Raymond Hettinger for use with 5940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# the Mersenne Twister core generator. 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 61145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerimport _random 6240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettingerclass Random(_random.Random): 64c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Random number generator base class used by bound module functions. 65c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 66c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Used to instantiate instances of Random to get generators that don't 67c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger share state. Especially useful for multi-threaded programs, creating 68c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger a different instance of Random for each thread, and using the jumpahead() 69c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger method to ensure that the generated sequences seen by each thread don't 70c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger overlap. 71c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 72c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Class Random can also be subclassed if you want to use a different basic 73c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger generator of your own devising: in that case, override the following 74c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger methods: random(), seed(), getstate(), setstate() and jumpahead(). 75ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 76c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 7733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 7840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 2 # used by getstate/setstate 7933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 80d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def __init__(self, x=None): 81d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Initialize an instance. 8233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 83d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional argument x controls seeding, as for Random.seed(). 84d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 8533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 86d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.seed(x) 8740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 88ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 890de88fc4b108751b86443852b6741680d704168fTim Peters def seed(self, a=None): 900de88fc4b108751b86443852b6741680d704168fTim Peters """Initialize internal state from hashable object. 91d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 920de88fc4b108751b86443852b6741680d704168fTim Peters None or no argument seeds from current time. 930de88fc4b108751b86443852b6741680d704168fTim Peters 94bcd725fc456faca13f4598f87c0517f917711cdaTim Peters If a is not None or an int or long, hash(a) is used instead. 95d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 96d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 97145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).seed(a) 9846c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters self.gauss_next = None 9946c04e140cf26d1b44935c28c6f15ea467400d22Tim Peters 100d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def getstate(self): 101d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Return internal state; can be passed to setstate() later.""" 102145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger return self.VERSION, super(Random, self).getstate(), self.gauss_next 103d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 104d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def setstate(self, state): 105d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Restore internal state from object returned by getstate().""" 106d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters version = state[0] 10740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 2: 10840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, internalstate, self.gauss_next = state 109145a4a0f10009f7ce2644465ccd359938b034ac4Raymond Hettinger super(Random, self).setstate(internalstate) 110d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 111d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError("state with version %s passed to " 112d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters "Random.setstate() of version %s" % 113d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters (version, self.VERSION)) 114d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 115cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- Methods below this point do not need to be overridden when 116cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## ---- subclassing for the purpose of using a different core generator. 117d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 118cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- pickle support ------------------- 119d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 120cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __getstate__(self): # for pickle 121cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters return self.getstate() 122d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 123cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters def __setstate__(self, state): # for pickle 124cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters self.setstate(state) 125cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 126cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- integer methods ------------------- 127d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 128d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randrange(self, start, stop=None, step=1, int=int, default=None): 129d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random item from range(start, stop[, step]). 130d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 131d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters This fixes the problem with randint() which includes the 132d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters endpoint; in Python this is usually not what you want. 133d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Do not supply the 'int' and 'default' arguments. 134d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 135d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 136d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # This code is a bit messy to make it fast for the 1379146f27b7799dab231083f194a14c6157b57549fTim Peters # common case while still doing adequate error checking. 138d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istart = int(start) 139d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart != start: 140d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer arg 1 for randrange()" 141d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if stop is default: 142d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istart > 0: 143d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return int(self.random() * istart) 144d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1459146f27b7799dab231083f194a14c6157b57549fTim Peters 1469146f27b7799dab231083f194a14c6157b57549fTim Peters # stop argument supplied. 147d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istop = int(stop) 148d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istop != stop: 149d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer stop for randrange()" 1509146f27b7799dab231083f194a14c6157b57549fTim Peters if step == 1 and istart < istop: 15176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # Note that 15276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # int(istart + self.random()*(istop - istart)) 15376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # instead would be incorrect. For example, consider istart 15476ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # = -2 and istop = 0. Then the guts would be in 15576ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # -2.0 to 0.0 exclusive on both ends (ignoring that random() 15676ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # might return 0.0), and because int() truncates toward 0, the 15776ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # final result would be -1 or 0 (instead of -2 or -1). 15876ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # istart + int(self.random()*(istop - istart)) 15976ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # would also be incorrect, for a subtler reason: the RHS 16076ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # can return a long, and then randrange() would also return 16176ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # a long, but we're supposed to return an int (for backward 16276ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters # compatibility). 16376ca1d428f96284ed58f4523b698ed95c6fdbdb2Tim Peters return int(istart + int(self.random()*(istop - istart))) 164d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if step == 1: 165d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 1669146f27b7799dab231083f194a14c6157b57549fTim Peters 1679146f27b7799dab231083f194a14c6157b57549fTim Peters # Non-unit step argument supplied. 168d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters istep = int(step) 169d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep != step: 170d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "non-integer step for randrange()" 171d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if istep > 0: 172d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep - 1) / istep 173d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif istep < 0: 174d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters n = (istop - istart + istep + 1) / istep 175d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 176d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "zero step for randrange()" 177d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 178d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if n <= 0: 179d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters raise ValueError, "empty range for randrange()" 180d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return istart + istep*int(self.random() * n) 181d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 182d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def randint(self, a, b): 183cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters """Return random integer in range [a, b], including both end points. 184d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 185d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 186d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return self.randrange(a, b+1) 187d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 188cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- sequence methods ------------------- 189cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 190d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def choice(self, seq): 191d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Choose a random element from a non-empty sequence.""" 192d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return seq[int(self.random() * len(seq))] 193d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 194d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def shuffle(self, x, random=None, int=int): 195d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """x, random=random.random -> shuffle list x in place; return None. 196d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 197d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Optional arg random is a 0-argument function returning a random 198d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters float in [0.0, 1.0); by default, the standard random.random. 199d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 200d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters Note that for even rather small len(x), the total number of 201d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters permutations of x is larger than the period of most random number 202d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters generators; this implies that "most" permutations of a long 203d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters sequence can never be generated. 204d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """ 205d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 206d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if random is None: 207d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 208d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters for i in xrange(len(x)-1, 0, -1): 209cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters # pick an element in x[:i+1] with which to exchange x[i] 210d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters j = int(random() * (i+1)) 211d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x[i], x[j] = x[j], x[i] 212d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 213fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger def sample(self, population, k): 214f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """Chooses k unique random elements from a population sequence. 215f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 216c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Returns a new list containing elements from the population while 217c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger leaving the original population unchanged. The resulting list is 218c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger in selection order so that all sub-slices will also be valid random 219c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger samples. This allows raffle winners (the sample) to be partitioned 220c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger into grand prize and second place winners (the subslices). 221f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 222c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger Members of the population need not be hashable or unique. If the 223c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger population contains repeats, then each occurrence is a possible 224c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger selection in the sample. 225f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 226c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger To choose a sample in a range of integers, use xrange as an argument. 227c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger This is especially fast and space efficient for sampling from a 228c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger large population: sample(xrange(10000000), 60) 229f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger """ 230f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 231c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger # Sampling without replacement entails tracking either potential 2328b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections (the pool) in a list or previous selections in a 2338b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary. 234c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 2358b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # When the number of selections is small compared to the population, 2368b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # then tracking selections is efficient, requiring only a small 2378b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # dictionary and an occasional reselection. For a larger number of 2388b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # selections, the pool tracking method is preferred since the list takes 2398b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # less space than the dictionary and it doesn't suffer from frequent 2408b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger # reselections. 241c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger 242f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger n = len(population) 243f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if not 0 <= k <= n: 244f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger raise ValueError, "sample larger than population" 2458b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger random = self.random 246fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger _int = int 247c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result = [None] * k 248f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger if n < 6 * k: # if n len list takes less space than a k len dict 249311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger pool = list(population) 250311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger for i in xrange(k): # invariant: non-selected at [0,n-i) 251fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * (n-i)) 252311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger result[i] = pool[j] 2538b9aa8dbba644da23ce8417f2d30b218392b3282Raymond Hettinger pool[j] = pool[n-i-1] # move non-selected item into vacancy 254c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger else: 255311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger selected = {} 256c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger for i in xrange(k): 257fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 258311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while j in selected: 259fdbe5223b7402ee34c4f0c06caa6faabd9e73e35Raymond Hettinger j = _int(random() * n) 260c0b4034b8165ce958a23f2c865b51ae0f52040f5Raymond Hettinger result[i] = selected[j] = population[j] 261311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger return result 262f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 263cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- real-valued distributions ------------------- 264cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 265cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- uniform distribution ------------------- 266d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 267d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def uniform(self, a, b): 268d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters """Get a random number in the range [a, b).""" 269d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return a + (b-a) * self.random() 270ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 271cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- normal distribution -------------------- 272ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 273d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def normalvariate(self, mu, sigma): 274c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Normal distribution. 275c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 276c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. 277ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 278c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 279d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu = mean, sigma = standard deviation 280d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 281d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses Kinderman and Monahan method. Reference: Kinderman, 282d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # A.J. and Monahan, J.F., "Computer generation of random 283d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables using the ratio of uniform deviates", ACM Trans 284d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Math Software, 3, (1977), pp257-260. 285d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 286d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 287311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 288d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 28973ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 290d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = NV_MAGICCONST*(u1-0.5)/u2 291d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters zz = z*z/4.0 292d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if zz <= -_log(u2): 293d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 294d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 295ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 296cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- lognormal distribution -------------------- 297ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 298d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def lognormvariate(self, mu, sigma): 299c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Log normal distribution. 300c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 301c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger If you take the natural logarithm of this distribution, you'll get a 302c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger normal distribution with mean mu and standard deviation sigma. 303c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu can have any value, and sigma must be greater than zero. 304ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 305c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 306d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return _exp(self.normalvariate(mu, sigma)) 307ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 308cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- circular uniform -------------------- 309ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 310d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def cunifvariate(self, mean, arc): 311c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular uniform distribution. 312c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 313c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mean is the mean angle, and arc is the range of the distribution, 314c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger centered around the mean angle. Both values must be expressed in 315c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger radians. Returned values range between mean - arc/2 and 316c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mean + arc/2 and are normalized to between 0 and pi. 317c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 318c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Deprecated in version 2.3. Use: 319c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger (mean + arc * (Random.random() - 0.5)) % Math.pi 320ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 321c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 322d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mean: mean angle (in radians between 0 and pi) 323d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # arc: range of distribution (in radians between 0 and pi) 324c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger import warnings 325c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger warnings.warn("The cunifvariate function is deprecated; Use (mean " 326c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger "+ arc * (Random.random() - 0.5)) % Math.pi instead", 327c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger DeprecationWarning) 328ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 329d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return (mean + arc * (self.random() - 0.5)) % _pi 330ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 331cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- exponential distribution -------------------- 332ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 333d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def expovariate(self, lambd): 334c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Exponential distribution. 335c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 336c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger lambd is 1.0 divided by the desired mean. (The parameter would be 337c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger called "lambda", but that is a reserved word in Python.) Returned 338c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger values range from 0 to positive infinity. 339ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 340c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 341d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lambd: rate lambd = 1/mean 342d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # ('lambda' is a Python reserved word) 343ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 344d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 3450c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 346d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 347d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 348d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return -_log(u)/lambd 349ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 350cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- von Mises distribution -------------------- 351ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 352d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def vonmisesvariate(self, mu, kappa): 353c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Circular data distribution. 354ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 355c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean angle, expressed in radians between 0 and 2*pi, and 356c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger kappa is the concentration parameter, which must be greater than or 357c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger equal to zero. If kappa is equal to zero, this distribution reduces 358c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger to a uniform random angle over the range 0 to 2*pi. 359ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 360c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 361d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # mu: mean angle (in radians between 0 and 2*pi) 362d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # kappa: concentration parameter kappa (>= 0) 363d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # if kappa = 0 generate uniform random angle 3645810297052003f28788f6790ac799fe8e5373494Guido van Rossum 365d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Based upon an algorithm published in: Fisher, N.I., 366d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # "Statistical Analysis of Circular Data", Cambridge 367d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # University Press, 1993. 3685810297052003f28788f6790ac799fe8e5373494Guido van Rossum 369d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Thanks to Magnus Kessler for a correction to the 370d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # implementation of step 4. 3715810297052003f28788f6790ac799fe8e5373494Guido van Rossum 372d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 373d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if kappa <= 1e-6: 374d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return TWOPI * random() 375ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 376d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) 377d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (a - _sqrt(2.0 * a))/(2.0 * kappa) 378d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = (1.0 + b * b)/(2.0 * b) 379ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 380311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 381d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 382ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 383d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(_pi * u1) 384d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters f = (1.0 + r * z)/(r + z) 385d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters c = kappa * (r - f) 386ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 387d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u2 = random() 388ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 389d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): 390d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 391ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 392d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u3 = random() 393d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if u3 > 0.5: 394d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) + _acos(f) 395d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 396d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters theta = (mu % TWOPI) - _acos(f) 397ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 398d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return theta 399ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 400cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- gamma distribution -------------------- 401ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 402d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gammavariate(self, alpha, beta): 403c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gamma distribution. Not the gamma function! 404c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 405c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > 0 and beta > 0. 406c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 407c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 4088ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 409b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 4108ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 411570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # Warning: a few older sources define the gamma distribution in terms 412570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum # of alpha > -1.0 413570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum if alpha <= 0.0 or beta <= 0.0: 414570764ddce285afc32e6bd4bce031e421376b382Guido van Rossum raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 4158ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 416d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 417d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if alpha > 1.0: 418d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 419d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses R.C.H. Cheng, "The generation of Gamma 420d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # variables with non-integral shape parameters", 421d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Applied Statistics, (1977), 26, No. 1, p71-74 422d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 423ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ainv = _sqrt(2.0 * alpha - 1.0) 424ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger bbb = alpha - LOG4 425ca6cdc2c0259b1b74a3f4c2d29a35e76617a3019Raymond Hettinger ccc = alpha + ainv 4268ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters 427311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 428d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 42973ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger if not 1e-7 < u1 < .9999999: 43073ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger continue 43173ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u2 = 1.0 - random() 432d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters v = _log(u1/(1.0-u1))/ainv 433d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = alpha*_exp(v) 434d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = u1*u1*u2 435d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters r = bbb+ccc*v-x 436d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 437b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 438d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 439d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters elif alpha == 1.0: 440d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # expovariate(1) 4410c9886d589ddebf32de0ca3f027a173222ed383aTim Peters u = random() 442d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters while u <= 1e-7: 443d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 444b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return -_log(u) * beta 445d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 446d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: # alpha is between 0 and 1 (exclusive) 447d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 448d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 449d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 450311f4196284b894f86f56c287c71a0e59c4a72a2Raymond Hettinger while True: 451d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u = random() 452d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters b = (_e + alpha)/_e 453d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters p = b*u 454d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if p <= 1.0: 455d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = pow(p, 1.0/alpha) 456d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters else: 457d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # p > 1 458d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x = -_log((b-p)/alpha) 459d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters u1 = random() 460d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if not (((p <= 1.0) and (u1 > _exp(-x))) or 461d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 462d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters break 463b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return x * beta 464b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 465b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 466b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger def stdgamma(self, alpha, ainv, bbb, ccc): 467b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # This method was (and shall remain) undocumented. 468b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # This method is deprecated 469b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # for the following reasons: 470b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 1. Returns same as .gammavariate(alpha, 1.0) 471b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 2. Requires caller to provide 3 extra arguments 472b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # that are functions of alpha anyway 473b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # 3. Can't be used for alpha < 0.5 474b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 475b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # ainv = sqrt(2 * alpha - 1) 476b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # bbb = alpha - log(4) 477b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger # ccc = alpha + ainv 478b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger import warnings 479b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger warnings.warn("The stdgamma function is deprecated; " 480b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger "use gammavariate() instead", 481b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger DeprecationWarning) 482b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger return self.gammavariate(alpha, 1.0) 483b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger 484ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 48595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 486cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Gauss (faster alternative) -------------------- 48795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 488d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def gauss(self, mu, sigma): 489c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Gaussian distribution. 490c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 491c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger mu is the mean, and sigma is the standard deviation. This is 492c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger slightly faster than the normalvariate() function. 493c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 494c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Not thread-safe without a lock around calls. 495ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 496c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 497d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 498d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # When x and y are two variables from [0, 1), uniformly 499d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # distributed, then 500d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 501d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # cos(2*pi*x)*sqrt(-2*log(1-y)) 502d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # sin(2*pi*x)*sqrt(-2*log(1-y)) 503d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # 504d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # are two *independent* variables with normal distribution 505d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (mu = 0, sigma = 1). 506d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (Lambert Meertens) 507d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # (corrected version; bug discovered by Mike Miller, fixed by LM) 508d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 509d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Multithreading note: When two threads call this function 510d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # simultaneously, it is possible that they will receive the 511d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # same return value. The window is very small though. To 512d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # avoid this, you have to use a lock around all calls. (I 513d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # didn't want to slow this down in the serial case by using a 514d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # lock here.) 515d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 516d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters random = self.random 517d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = self.gauss_next 518d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = None 519d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters if z is None: 520d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters x2pi = random() * TWOPI 521d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters g2rad = _sqrt(-2.0 * _log(1.0 - random())) 522d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters z = _cos(x2pi) * g2rad 523d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters self.gauss_next = _sin(x2pi) * g2rad 524d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 525d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return mu + z*sigma 52695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 527cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- beta -------------------- 52885e2e4742d0a1accecd02058a7907df36308297eTim Peters## See 52985e2e4742d0a1accecd02058a7907df36308297eTim Peters## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 53085e2e4742d0a1accecd02058a7907df36308297eTim Peters## for Ivan Frohne's insightful analysis of why the original implementation: 53185e2e4742d0a1accecd02058a7907df36308297eTim Peters## 53285e2e4742d0a1accecd02058a7907df36308297eTim Peters## def betavariate(self, alpha, beta): 53385e2e4742d0a1accecd02058a7907df36308297eTim Peters## # Discrete Event Simulation in C, pp 87-88. 53485e2e4742d0a1accecd02058a7907df36308297eTim Peters## 53585e2e4742d0a1accecd02058a7907df36308297eTim Peters## y = self.expovariate(alpha) 53685e2e4742d0a1accecd02058a7907df36308297eTim Peters## z = self.expovariate(1.0/beta) 53785e2e4742d0a1accecd02058a7907df36308297eTim Peters## return z/(y+z) 53885e2e4742d0a1accecd02058a7907df36308297eTim Peters## 53985e2e4742d0a1accecd02058a7907df36308297eTim Peters## was dead wrong, and how it probably got that way. 54095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 541d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def betavariate(self, alpha, beta): 542c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Beta distribution. 543c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 544c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Conditions on the parameters are alpha > -1 and beta} > -1. 545c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger Returned values range between 0 and 1. 546ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 547c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 548ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 54985e2e4742d0a1accecd02058a7907df36308297eTim Peters # This version due to Janne Sinkkonen, and matches all the std 55085e2e4742d0a1accecd02058a7907df36308297eTim Peters # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 55185e2e4742d0a1accecd02058a7907df36308297eTim Peters y = self.gammavariate(alpha, 1.) 55285e2e4742d0a1accecd02058a7907df36308297eTim Peters if y == 0: 55385e2e4742d0a1accecd02058a7907df36308297eTim Peters return 0.0 55485e2e4742d0a1accecd02058a7907df36308297eTim Peters else: 55585e2e4742d0a1accecd02058a7907df36308297eTim Peters return y / (y + self.gammavariate(beta, 1.)) 55695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 557cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Pareto -------------------- 558cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 559d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def paretovariate(self, alpha): 560c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Pareto distribution. alpha is the shape parameter.""" 561d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 495 562cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 56373ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 564d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return 1.0 / pow(u, 1.0/alpha) 565cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 566cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- Weibull -------------------- 567cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 568d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters def weibullvariate(self, alpha, beta): 569c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """Weibull distribution. 570c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger 571c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger alpha is the scale parameter and beta is the shape parameter. 572ef4d4bdc3c99d3120a401b7af6e06610716d2e47Raymond Hettinger 573c32f0336e062dd36f82fb236c66ac25e2bac217bRaymond Hettinger """ 574d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters # Jain, pg. 499; bug fix courtesy Bill Arms 575cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 57673ced7ee995180c0bd8d96ff7d7fb614a744ad7dRaymond Hettinger u = 1.0 - self.random() 577d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters return alpha * pow(-_log(u), 1.0/beta) 5786c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 57940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger## -------------------- Wichmann-Hill ------------------- 58040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 58140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettingerclass WichmannHill(Random): 58240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 58340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger VERSION = 1 # used by getstate/setstate 58440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 58540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def seed(self, a=None): 58640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Initialize internal state from hashable object. 58740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 58840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. 58940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is not None or an int or long, hash(a) is used instead. 59140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger If a is an int or long, a is used directly. Distinct values between 59340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 0 and 27814431486575L inclusive are guaranteed to yield distinct 59440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger internal states (this guarantee is specific to the default 59540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Wichmann-Hill generator). 59640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 59740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 59840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 59940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 60040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 60140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = long(time.time() * 256) 60240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not isinstance(a, (int, long)): 60440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 60540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 60640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 30268) 60740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 30306) 60840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 30322) 60940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = int(x)+1, int(y)+1, int(z)+1 61040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 61240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def random(self): 61440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Get the next random number in the range [0.0, 1.0).""" 61540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 61640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichman-Hill random number generator. 61740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 61840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Wichmann, B. A. & Hill, I. D. (1982) 61940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Algorithm AS 183: 62040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # An efficient and portable pseudo-random number generator 62140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 31 (1982) 188-190 62240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 62340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # see also: 62440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Correction to Algorithm AS 183 62540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 33 (1984) 123 62640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # 62740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # McLeod, A. I. (1985) 62840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # A remark on Algorithm AS 183 62940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Applied Statistics 34 (1985),198-200 63040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 63140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # This part is thread-unsafe: 63240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # BEGIN CRITICAL SECTION 63340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 63440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (171 * x) % 30269 63540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (172 * y) % 30307 63640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (170 * z) % 30323 63740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 63840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # END CRITICAL SECTION 63940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Note: on a platform using IEEE-754 double arithmetic, this can 64140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # never return 0.0 (asserted by Tim; proof too long for a comment). 64240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 64340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def getstate(self): 64540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Return internal state; can be passed to setstate() later.""" 64640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return self.VERSION, self._seed, self.gauss_next 64740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 64840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def setstate(self, state): 64940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Restore internal state from object returned by getstate().""" 65040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version = state[0] 65140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if version == 1: 65240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger version, self._seed, self.gauss_next = state 65340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger else: 65440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("state with version %s passed to " 65540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger "Random.setstate() of version %s" % 65640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger (version, self.VERSION)) 65740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 65840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def jumpahead(self, n): 65940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Act as if n calls to random() were made, but quickly. 66040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger n is an int, greater than or equal to 0. 66240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 66340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Example use: If you have 2 threads and know that each will 66440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger consume no more than a million random numbers, create two Random 66540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger objects r1 and r2, then do 66640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.setstate(r1.getstate()) 66740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger r2.jumpahead(1000000) 66840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger Then r1 and r2 will use guaranteed-disjoint segments of the full 66940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger period. 67040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 67140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 67240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not n >= 0: 67340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError("n must be >= 0") 67440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x, y, z = self._seed 67540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = int(x * pow(171, n, 30269)) % 30269 67640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = int(y * pow(172, n, 30307)) % 30307 67740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = int(z * pow(170, n, 30323)) % 30323 67840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = x, y, z 67940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def __whseed(self, x=0, y=0, z=0): 68140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Set the Wichmann-Hill seed from (x, y, z). 68240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger These must be integers in the range [0, 256). 68440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 68540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 68640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not type(x) == type(y) == type(z) == int: 68740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise TypeError('seeds must be integers') 68840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 68940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger raise ValueError('seeds must be in range(0, 256)') 69040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if 0 == x == y == z: 69140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Initialize from current time 69240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger import time 69340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = long(time.time() * 256) 69440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t = int((t&0xffffff) ^ (t>>24)) 69540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, x = divmod(t, 256) 69640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, y = divmod(t, 256) 69740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger t, z = divmod(t, 256) 69840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger # Zero is a poor seed, so substitute 1 69940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self._seed = (x or 1, y or 1, z or 1) 70040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.gauss_next = None 70240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger def whseed(self, a=None): 70440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """Seed from hashable object's hash code. 70540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 70640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger None or no argument seeds from current time. It is not guaranteed 70740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger that objects with distinct hash codes lead to distinct internal 70840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger states. 70940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger This is obsolete, provided for compatibility with the seed routine 71140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger used prior to Python 2.1. Use the .seed() method instead. 71240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger """ 71340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 71440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger if a is None: 71540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed() 71640f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger return 71740f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a = hash(a) 71840f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, x = divmod(a, 256) 71940f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, y = divmod(a, 256) 72040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger a, z = divmod(a, 256) 72140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger x = (x + a) % 256 or 1 72240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger y = (y + a) % 256 or 1 72340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger z = (z + a) % 256 or 1 72440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger self.__whseed(x, y, z) 72540f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 726cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters## -------------------- test program -------------------- 727ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 728d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersdef _test_generator(n, funccall): 7290c9886d589ddebf32de0ca3f027a173222ed383aTim Peters import time 7300c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print n, 'times', funccall 7310c9886d589ddebf32de0ca3f027a173222ed383aTim Peters code = compile(funccall, funccall, 'eval') 732b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total = 0.0 7330c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = 0.0 7340c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = 1e10 7350c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = -1e10 7360c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t0 = time.time() 7370c9886d589ddebf32de0ca3f027a173222ed383aTim Peters for i in range(n): 7380c9886d589ddebf32de0ca3f027a173222ed383aTim Peters x = eval(code) 739b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger total += x 7400c9886d589ddebf32de0ca3f027a173222ed383aTim Peters sqsum = sqsum + x*x 7410c9886d589ddebf32de0ca3f027a173222ed383aTim Peters smallest = min(x, smallest) 7420c9886d589ddebf32de0ca3f027a173222ed383aTim Peters largest = max(x, largest) 7430c9886d589ddebf32de0ca3f027a173222ed383aTim Peters t1 = time.time() 7440c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print round(t1-t0, 3), 'sec,', 745b98154e4243a8d73f758dfee9a81bbe36ddc05cbRaymond Hettinger avg = total/n 746d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters stddev = _sqrt(sqsum/n - avg*avg) 7470c9886d589ddebf32de0ca3f027a173222ed383aTim Peters print 'avg %g, stddev %g, min %g, max %g' % \ 7480c9886d589ddebf32de0ca3f027a173222ed383aTim Peters (avg, stddev, smallest, largest) 749ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 750f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettinger 751f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingerdef _test(N=2000): 752d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'random()') 753d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'normalvariate(0.0, 1.0)') 754d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'lognormvariate(0.0, 1.0)') 755d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'cunifvariate(0.0, 1.0)') 756d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'vonmisesvariate(0.0, 1.0)') 757b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.01, 1.0)') 758b760efb08d509bb2acfdef8f4e59dfafa20ca57fRaymond Hettinger _test_generator(N, 'gammavariate(0.1, 1.0)') 7598ac1495a6a1d18111a626cec0c7f2eb67df3edb3Tim Peters _test_generator(N, 'gammavariate(0.1, 2.0)') 760d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.5, 1.0)') 761d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(0.9, 1.0)') 762d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(1.0, 1.0)') 763d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(2.0, 1.0)') 764d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(20.0, 1.0)') 765d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gammavariate(200.0, 1.0)') 766d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'gauss(0.0, 1.0)') 767d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test_generator(N, 'betavariate(3.0, 3.0)') 768cd804108548e1939bd8646634ed52ef388ee9f44Tim Peters 769715c4c412b21f68ad59773698d06eea8eb0c5a44Tim Peters# Create one instance, seeded from current time, and export its methods 77040f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# as module-level functions. The functions share state across all uses 77140f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger#(both in the user's code and in the Python libraries), but that's fine 77240f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# for most programs and is easier for the casual user than making them 77340f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger# instantiate their own Random() instance. 77440f621709286a7a0f7e6f260c0fd020d0fac0de0Raymond Hettinger 775d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters_inst = Random() 776d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersseed = _inst.seed 777d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandom = _inst.random 778d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersuniform = _inst.uniform 779d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandint = _inst.randint 780d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterschoice = _inst.choice 781d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersrandrange = _inst.randrange 782f24eb35d185c0623315cfbd9977d37c509860dcfRaymond Hettingersample = _inst.sample 783d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersshuffle = _inst.shuffle 784d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersnormalvariate = _inst.normalvariate 785d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterslognormvariate = _inst.lognormvariate 786d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterscunifvariate = _inst.cunifvariate 787d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersexpovariate = _inst.expovariate 788d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersvonmisesvariate = _inst.vonmisesvariate 789d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgammavariate = _inst.gammavariate 790d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersstdgamma = _inst.stdgamma 791d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgauss = _inst.gauss 792d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersbetavariate = _inst.betavariate 793d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersparetovariate = _inst.paretovariate 794d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersweibullvariate = _inst.weibullvariate 795d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Petersgetstate = _inst.getstate 796d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peterssetstate = _inst.setstate 797d52269bfd029c4a517ea74c17edd5c3a250c366cTim Petersjumpahead = _inst.jumpahead 798d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters 799ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 800d7b5e88e8e40b77813ceb25dc28b87d672538403Tim Peters _test() 801