13257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel"""Random variable generators. 23257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 33257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel integers 43257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel -------- 53257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel uniform within range 63257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 73257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel sequences 83257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel --------- 93257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pick random element 103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pick random sample 113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel generate random permutation 123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel distributions on the real line: 143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel ------------------------------ 153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel uniform 163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel triangular 173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel normal (Gaussian) 183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel lognormal 193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel negative exponential 203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel gamma 213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel beta 223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pareto 233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Weibull 243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel distributions on the circle (angles 0 to 2pi) 263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel --------------------------------------------- 273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel circular uniform 283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel von Mises 293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielGeneral notes on the underlying Mersenne Twister core generator: 313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel* The period is 2**19937-1. 333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel* It is one of the most extensively tested generators in existence. 343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel* Without a direct way to compute N steps forward, the semantics of 353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel jumpahead(n) are weakened to simply jump to another distant state and rely 363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel on the large period to avoid overlapping sequences. 373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel* The random() method is implemented in C, executes in a single Python step, 383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel and is, therefore, threadsafe. 393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel""" 413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom __future__ import division 433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom warnings import warn as _warn 443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType 453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil 463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin 473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom os import urandom as _urandom 483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielfrom binascii import hexlify as _hexlify 493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielimport hashlib as _hashlib 503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel__all__ = ["Random","seed","random","uniform","randint","choice","sample", 523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "randrange","shuffle","normalvariate","lognormvariate", 533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "expovariate","vonmisesvariate","gammavariate","triangular", 543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "gauss","betavariate","paretovariate","weibullvariate", 553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", 563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "SystemRandom"] 573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielNV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) 593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielTWOPI = 2.0*_pi 603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielLOG4 = _log(4.0) 613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielSG_MAGICCONST = 1.0 + _log(4.5) 623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielBPF = 53 # Number of bits in a float 633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielRECIP_BPF = 2**-BPF 643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# Translated by Guido van Rossum from C source provided by 673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# Adrian Baddeley. Adapted by Raymond Hettinger for use with 683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# the Mersenne Twister and os.urandom() core generators. 693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielimport _random 713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielclass Random(_random.Random): 733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Random number generator base class used by bound module functions. 743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Used to instantiate instances of Random to get generators that don't 763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel share state. Especially useful for multi-threaded programs, creating 773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a different instance of Random for each thread, and using the jumpahead() 783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel method to ensure that the generated sequences seen by each thread don't 793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel overlap. 803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Class Random can also be subclassed if you want to use a different basic 823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel generator of your own devising: in that case, override the following 833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel methods: random(), seed(), getstate(), setstate() and jumpahead(). 843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Optionally, implement a getrandbits() method so that randrange() can cover 853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel arbitrarily large ranges. 863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel VERSION = 3 # used by getstate/setstate 903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def __init__(self, x=None): 923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Initialize an instance. 933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Optional argument x controls seeding, as for Random.seed(). 953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.seed(x) 983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = None 993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def seed(self, a=None): 1013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Initialize internal state from hashable object. 1023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel None or no argument seeds from current time or from an operating 1043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel system specific randomness source if available. 1053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel If a is not None or an int or long, hash(a) is used instead. 1073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 1083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if a is None: 1103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 1113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Seed with enough bytes to span the 19937 bit 1123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # state space for the Mersenne Twister 1133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = long(_hexlify(_urandom(2500)), 16) 1143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except NotImplementedError: 1153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel import time 1163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = long(time.time() * 256) # use fractional seconds 1173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel super(Random, self).seed(a) 1193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = None 1203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def getstate(self): 1223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Return internal state; can be passed to setstate() later.""" 1233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.VERSION, super(Random, self).getstate(), self.gauss_next 1243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def setstate(self, state): 1263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Restore internal state from object returned by getstate().""" 1273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel version = state[0] 1283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if version == 3: 1293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel version, internalstate, self.gauss_next = state 1303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel super(Random, self).setstate(internalstate) 1313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel elif version == 2: 1323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel version, internalstate, self.gauss_next = state 1333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # In version 2, the state was saved as signed ints, which causes 1343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # inconsistencies between 32/64-bit systems. The state is 1353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # really unsigned 32-bit ints, so we convert negative ints from 1363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # version 2 to positive longs for version 3. 1373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 1383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel internalstate = tuple( long(x) % (2**32) for x in internalstate ) 1393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except ValueError, e: 1403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise TypeError, e 1413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel super(Random, self).setstate(internalstate) 1423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 1433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError("state with version %s passed to " 1443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "Random.setstate() of version %s" % 1453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel (version, self.VERSION)) 1463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def jumpahead(self, n): 1483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Change the internal state to one that is likely far away 1493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel from the current state. This method will not be in Py3.x, 1503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel so it is better to simply reseed. 1513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 1523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # The super.jumpahead() method uses shuffling to change state, 1533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # so it needs a large and "interesting" n to work with. Here, 1543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # we use hashing to create a large n for the shuffle. 1553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel s = repr(n) + repr(self.getstate()) 1563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel n = int(_hashlib.new('sha512', s).hexdigest(), 16) 1573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel super(Random, self).jumpahead(n) 1583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## ---- Methods below this point do not need to be overridden when 1603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## ---- subclassing for the purpose of using a different core generator. 1613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- pickle support ------------------- 1633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def __getstate__(self): # for pickle 1653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.getstate() 1663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def __setstate__(self, state): # for pickle 1683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.setstate(state) 1693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def __reduce__(self): 1713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.__class__, (), self.getstate() 1723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- integer methods ------------------- 1743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def randrange(self, start, stop=None, step=1, _int=int, _maxwidth=1L<<BPF): 1763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Choose a random item from range(start, stop[, step]). 1773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel This fixes the problem with randint() which includes the 1793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel endpoint; in Python this is usually not what you want. 1803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 1823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # This code is a bit messy to make it fast for the 1843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # common case while still doing adequate error checking. 1853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel istart = _int(start) 1863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istart != start: 1873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "non-integer arg 1 for randrange()" 1883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if stop is None: 1893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istart > 0: 1903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istart >= _maxwidth: 1913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self._randbelow(istart) 1923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return _int(self.random() * istart) 1933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "empty range for randrange()" 1943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 1953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # stop argument supplied. 1963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel istop = _int(stop) 1973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istop != stop: 1983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "non-integer stop for randrange()" 1993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel width = istop - istart 2003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if step == 1 and width > 0: 2013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Note that 2023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # int(istart + self.random()*width) 2033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # instead would be incorrect. For example, consider istart 2043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # = -2 and istop = 0. Then the guts would be in 2053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # -2.0 to 0.0 exclusive on both ends (ignoring that random() 2063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # might return 0.0), and because int() truncates toward 0, the 2073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # final result would be -1 or 0 (instead of -2 or -1). 2083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # istart + int(self.random()*width) 2093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # would also be incorrect, for a subtler reason: the RHS 2103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # can return a long, and then randrange() would also return 2113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # a long, but we're supposed to return an int (for backward 2123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # compatibility). 2133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if width >= _maxwidth: 2153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return _int(istart + self._randbelow(width)) 2163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return _int(istart + _int(self.random()*width)) 2173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if step == 1: 2183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) 2193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Non-unit step argument supplied. 2213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel istep = _int(step) 2223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istep != step: 2233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "non-integer step for randrange()" 2243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if istep > 0: 2253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel n = (width + istep - 1) // istep 2263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel elif istep < 0: 2273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel n = (width + istep + 1) // istep 2283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 2293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "zero step for randrange()" 2303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if n <= 0: 2323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, "empty range for randrange()" 2333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if n >= _maxwidth: 2353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return istart + istep*self._randbelow(n) 2363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return istart + istep*_int(self.random() * n) 2373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def randint(self, a, b): 2393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Return random integer in range [a, b], including both end points. 2403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 2413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.randrange(a, b+1) 2433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def _randbelow(self, n, _log=_log, _int=int, _maxwidth=1L<<BPF, 2453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): 2463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Return a random int in the range [0,n) 2473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Handles the case where n has more bits than returned 2493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel by a single call to the underlying generator. 2503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 2513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 2533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel getrandbits = self.getrandbits 2543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except AttributeError: 2553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pass 2563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 2573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Only call self.getrandbits if the original random() builtin method 2583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # has not been overridden or if a new getrandbits() was supplied. 2593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # This assures that the two methods correspond. 2603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: 2613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel k = _int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) 2623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r = getrandbits(k) 2633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while r >= n: 2643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r = getrandbits(k) 2653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return r 2663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if n >= _maxwidth: 2673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _warn("Underlying random() generator does not supply \n" 2683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "enough bits to choose from a population range this large") 2693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return _int(self.random() * n) 2703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- sequence methods ------------------- 2723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def choice(self, seq): 2743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Choose a random element from a non-empty sequence.""" 2753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty 2763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def shuffle(self, x, random=None): 2783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """x, random=random.random -> shuffle list x in place; return None. 2793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Optional arg random is a 0-argument function returning a random 2813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel float in [0.0, 1.0); by default, the standard random.random. 2823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 2843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if random is None: 2863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 2873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _int = int 2883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel for i in reversed(xrange(1, len(x))): 2893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # pick an element in x[:i+1] with which to exchange x[i] 2903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel j = _int(random() * (i+1)) 2913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x[i], x[j] = x[j], x[i] 2923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def sample(self, population, k): 2943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Chooses k unique random elements from a population sequence. 2953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 2963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Returns a new list containing elements from the population while 2973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel leaving the original population unchanged. The resulting list is 2983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel in selection order so that all sub-slices will also be valid random 2993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel samples. This allows raffle winners (the sample) to be partitioned 3003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel into grand prize and second place winners (the subslices). 3013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Members of the population need not be hashable or unique. If the 3033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel population contains repeats, then each occurrence is a possible 3043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel selection in the sample. 3053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel To choose a sample in a range of integers, use xrange as an argument. 3073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel This is especially fast and space efficient for sampling from a 3083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel large population: sample(xrange(10000000), 60) 3093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 3103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Sampling without replacement entails tracking either potential 3123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # selections (the pool) in a list or previous selections in a set. 3133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # When the number of selections is small compared to the 3153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # population, then tracking selections is efficient, requiring 3163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # only a small set and an occasional reselection. For 3173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # a larger number of selections, the pool tracking method is 3183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # preferred since the list takes less space than the 3193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # set and it doesn't suffer from frequent reselections. 3203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel n = len(population) 3223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not 0 <= k <= n: 3233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError("sample larger than population") 3243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 3253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _int = int 3263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel result = [None] * k 3273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel setsize = 21 # size of a small set minus size of an empty list 3283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if k > 5: 3293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets 3303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if n <= setsize or hasattr(population, "keys"): 3313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # An n-length list is smaller than a k-length set, or this is a 3323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # mapping type so the other algorithm wouldn't work. 3333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pool = list(population) 3343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel for i in xrange(k): # invariant: non-selected at [0,n-i) 3353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel j = _int(random() * (n-i)) 3363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel result[i] = pool[j] 3373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pool[j] = pool[n-i-1] # move non-selected item into vacancy 3383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 3393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 3403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel selected = set() 3413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel selected_add = selected.add 3423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel for i in xrange(k): 3433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel j = _int(random() * n) 3443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while j in selected: 3453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel j = _int(random() * n) 3463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel selected_add(j) 3473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel result[i] = population[j] 3483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except (TypeError, KeyError): # handle (at least) sets 3493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if isinstance(population, list): 3503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise 3513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.sample(tuple(population), k) 3523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return result 3533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- real-valued distributions ------------------- 3553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- uniform distribution ------------------- 3573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def uniform(self, a, b): 3593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "Get a random number in the range [a, b) or [a, b] depending on rounding." 3603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return a + (b-a) * self.random() 3613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- triangular -------------------- 3633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def triangular(self, low=0.0, high=1.0, mode=None): 3653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Triangular distribution. 3663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Continuous distribution bounded by given lower and upper limits, 3683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel and having a given mode value in-between. 3693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel http://en.wikipedia.org/wiki/Triangular_distribution 3713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 3733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = self.random() 3743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 3753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel c = 0.5 if mode is None else (mode - low) / (high - low) 3763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except ZeroDivisionError: 3773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return low 3783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if u > c: 3793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = 1.0 - u 3803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel c = 1.0 - c 3813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel low, high = high, low 3823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return low + (high - low) * (u * c) ** 0.5 3833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- normal distribution -------------------- 3853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def normalvariate(self, mu, sigma): 3873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Normal distribution. 3883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel mu is the mean, and sigma is the standard deviation. 3903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 3923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # mu = mean, sigma = standard deviation 3933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Uses Kinderman and Monahan method. Reference: Kinderman, 3953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # A.J. and Monahan, J.F., "Computer generation of random 3963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # variables using the ratio of uniform deviates", ACM Trans 3973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Math Software, 3, (1977), pp257-260. 3983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 3993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 4003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while 1: 4013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u1 = random() 4023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u2 = 1.0 - random() 4033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = NV_MAGICCONST*(u1-0.5)/u2 4043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel zz = z*z/4.0 4053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if zz <= -_log(u2): 4063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel break 4073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return mu + z*sigma 4083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- lognormal distribution -------------------- 4103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def lognormvariate(self, mu, sigma): 4123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Log normal distribution. 4133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel If you take the natural logarithm of this distribution, you'll get a 4153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel normal distribution with mean mu and standard deviation sigma. 4163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel mu can have any value, and sigma must be greater than zero. 4173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 4193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return _exp(self.normalvariate(mu, sigma)) 4203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- exponential distribution -------------------- 4223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def expovariate(self, lambd): 4243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Exponential distribution. 4253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel lambd is 1.0 divided by the desired mean. It should be 4273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel nonzero. (The parameter would be called "lambda", but that is 4283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a reserved word in Python.) Returned values range from 0 to 4293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel positive infinity if lambd is positive, and from negative 4303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel infinity to 0 if lambd is negative. 4313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 4333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # lambd: rate lambd = 1/mean 4343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # ('lambda' is a Python reserved word) 4353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # we use 1-random() instead of random() to preclude the 4373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # possibility of taking the log of zero. 4383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return -_log(1.0 - self.random())/lambd 4393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- von Mises distribution -------------------- 4413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def vonmisesvariate(self, mu, kappa): 4433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Circular data distribution. 4443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel mu is the mean angle, expressed in radians between 0 and 2*pi, and 4463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel kappa is the concentration parameter, which must be greater than or 4473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel equal to zero. If kappa is equal to zero, this distribution reduces 4483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel to a uniform random angle over the range 0 to 2*pi. 4493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 4513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # mu: mean angle (in radians between 0 and 2*pi) 4523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # kappa: concentration parameter kappa (>= 0) 4533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # if kappa = 0 generate uniform random angle 4543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Based upon an algorithm published in: Fisher, N.I., 4563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # "Statistical Analysis of Circular Data", Cambridge 4573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # University Press, 1993. 4583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Thanks to Magnus Kessler for a correction to the 4603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # implementation of step 4. 4613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 4633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if kappa <= 1e-6: 4643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return TWOPI * random() 4653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel s = 0.5 / kappa 4673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r = s + _sqrt(1.0 + s * s) 4683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while 1: 4703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u1 = random() 4713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = _cos(_pi * u1) 4723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel d = z / (r + z) 4743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u2 = random() 4753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d): 4763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel break 4773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel q = 1.0 / r 4793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel f = (q + z) / (1.0 + q * z) 4803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u3 = random() 4813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if u3 > 0.5: 4823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel theta = (mu + _acos(f)) % TWOPI 4833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 4843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel theta = (mu - _acos(f)) % TWOPI 4853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return theta 4873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- gamma distribution -------------------- 4893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def gammavariate(self, alpha, beta): 4913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Gamma distribution. Not the gamma function! 4923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Conditions on the parameters are alpha > 0 and beta > 0. 4943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel The probability distribution function is: 4963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 4973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x ** (alpha - 1) * math.exp(-x / beta) 4983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel pdf(x) = -------------------------------------- 4993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel math.gamma(alpha) * beta ** alpha 5003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 5023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 5043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Warning: a few older sources define the gamma distribution in terms 5063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # of alpha > -1.0 5073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if alpha <= 0.0 or beta <= 0.0: 5083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 5093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 5113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if alpha > 1.0: 5123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Uses R.C.H. Cheng, "The generation of Gamma 5143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # variables with non-integral shape parameters", 5153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Applied Statistics, (1977), 26, No. 1, p71-74 5163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel ainv = _sqrt(2.0 * alpha - 1.0) 5183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel bbb = alpha - LOG4 5193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel ccc = alpha + ainv 5203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while 1: 5223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u1 = random() 5233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not 1e-7 < u1 < .9999999: 5243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel continue 5253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u2 = 1.0 - random() 5263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel v = _log(u1/(1.0-u1))/ainv 5273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = alpha*_exp(v) 5283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = u1*u1*u2 5293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r = bbb+ccc*v-x 5303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 5313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return x * beta 5323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel elif alpha == 1.0: 5343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # expovariate(1) 5353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = random() 5363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while u <= 1e-7: 5373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = random() 5383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return -_log(u) * beta 5393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: # alpha is between 0 and 1 (exclusive) 5413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 5433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel while 1: 5453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = random() 5463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel b = (_e + alpha)/_e 5473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel p = b*u 5483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if p <= 1.0: 5493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = p ** (1.0/alpha) 5503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 5513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = -_log((b-p)/alpha) 5523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u1 = random() 5533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if p > 1.0: 5543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if u1 <= x ** (alpha - 1.0): 5553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel break 5563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel elif u1 <= _exp(-x): 5573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel break 5583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return x * beta 5593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- Gauss (faster alternative) -------------------- 5613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def gauss(self, mu, sigma): 5633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Gaussian distribution. 5643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel mu is the mean, and sigma is the standard deviation. This is 5663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel slightly faster than the normalvariate() function. 5673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Not thread-safe without a lock around calls. 5693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 5713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # When x and y are two variables from [0, 1), uniformly 5733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # distributed, then 5743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # 5753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # cos(2*pi*x)*sqrt(-2*log(1-y)) 5763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # sin(2*pi*x)*sqrt(-2*log(1-y)) 5773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # 5783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # are two *independent* variables with normal distribution 5793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # (mu = 0, sigma = 1). 5803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # (Lambert Meertens) 5813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # (corrected version; bug discovered by Mike Miller, fixed by LM) 5823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Multithreading note: When two threads call this function 5843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # simultaneously, it is possible that they will receive the 5853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # same return value. The window is very small though. To 5863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # avoid this, you have to use a lock around all calls. (I 5873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # didn't want to slow this down in the serial case by using a 5883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # lock here.) 5893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel random = self.random 5913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = self.gauss_next 5923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = None 5933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if z is None: 5943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x2pi = random() * TWOPI 5953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel g2rad = _sqrt(-2.0 * _log(1.0 - random())) 5963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = _cos(x2pi) * g2rad 5973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = _sin(x2pi) * g2rad 5983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 5993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return mu + z*sigma 6003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- beta -------------------- 6023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## See 6033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## http://mail.python.org/pipermail/python-bugs-list/2001-January/003752.html 6043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## for Ivan Frohne's insightful analysis of why the original implementation: 6053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## 6063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## def betavariate(self, alpha, beta): 6073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## # Discrete Event Simulation in C, pp 87-88. 6083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## 6093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## y = self.expovariate(alpha) 6103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## z = self.expovariate(1.0/beta) 6113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## return z/(y+z) 6123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## 6133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## was dead wrong, and how it probably got that way. 6143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def betavariate(self, alpha, beta): 6163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Beta distribution. 6173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Conditions on the parameters are alpha > 0 and beta > 0. 6193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Returned values range between 0 and 1. 6203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 6223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # This version due to Janne Sinkkonen, and matches all the std 6243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). 6253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel y = self.gammavariate(alpha, 1.) 6263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if y == 0: 6273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return 0.0 6283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 6293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return y / (y + self.gammavariate(beta, 1.)) 6303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- Pareto -------------------- 6323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def paretovariate(self, alpha): 6343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Pareto distribution. alpha is the shape parameter.""" 6353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Jain, pg. 495 6363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = 1.0 - self.random() 6383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return 1.0 / pow(u, 1.0/alpha) 6393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- Weibull -------------------- 6413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def weibullvariate(self, alpha, beta): 6433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Weibull distribution. 6443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel alpha is the scale parameter and beta is the shape parameter. 6463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 6483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Jain, pg. 499; bug fix courtesy Bill Arms 6493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel u = 1.0 - self.random() 6513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return alpha * pow(-_log(u), 1.0/beta) 6523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- Wichmann-Hill ------------------- 6543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielclass WichmannHill(Random): 6563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel VERSION = 1 # used by getstate/setstate 6583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def seed(self, a=None): 6603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Initialize internal state from hashable object. 6613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel None or no argument seeds from current time or from an operating 6633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel system specific randomness source if available. 6643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel If a is not None or an int or long, hash(a) is used instead. 6663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel If a is an int or long, a is used directly. Distinct values between 6683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 0 and 27814431486575L inclusive are guaranteed to yield distinct 6693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel internal states (this guarantee is specific to the default 6703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Wichmann-Hill generator). 6713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 6723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if a is None: 6743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel try: 6753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = long(_hexlify(_urandom(16)), 16) 6763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel except NotImplementedError: 6773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel import time 6783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = long(time.time() * 256) # use fractional seconds 6793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not isinstance(a, (int, long)): 6813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = hash(a) 6823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, x = divmod(a, 30268) 6843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, y = divmod(a, 30306) 6853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, z = divmod(a, 30322) 6863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self._seed = int(x)+1, int(y)+1, int(z)+1 6873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = None 6893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def random(self): 6913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Get the next random number in the range [0.0, 1.0).""" 6923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 6933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Wichman-Hill random number generator. 6943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # 6953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Wichmann, B. A. & Hill, I. D. (1982) 6963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Algorithm AS 183: 6973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # An efficient and portable pseudo-random number generator 6983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Applied Statistics 31 (1982) 188-190 6993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # 7003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # see also: 7013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Correction to Algorithm AS 183 7023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Applied Statistics 33 (1984) 123 7033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # 7043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # McLeod, A. I. (1985) 7053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # A remark on Algorithm AS 183 7063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Applied Statistics 34 (1985),198-200 7073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # This part is thread-unsafe: 7093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # BEGIN CRITICAL SECTION 7103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x, y, z = self._seed 7113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = (171 * x) % 30269 7123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel y = (172 * y) % 30307 7133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = (170 * z) % 30323 7143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self._seed = x, y, z 7153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # END CRITICAL SECTION 7163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Note: on a platform using IEEE-754 double arithmetic, this can 7183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # never return 0.0 (asserted by Tim; proof too long for a comment). 7193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 7203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def getstate(self): 7223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Return internal state; can be passed to setstate() later.""" 7233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return self.VERSION, self._seed, self.gauss_next 7243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def setstate(self, state): 7263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Restore internal state from object returned by getstate().""" 7273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel version = state[0] 7283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if version == 1: 7293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel version, self._seed, self.gauss_next = state 7303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel else: 7313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError("state with version %s passed to " 7323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "Random.setstate() of version %s" % 7333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel (version, self.VERSION)) 7343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def jumpahead(self, n): 7363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Act as if n calls to random() were made, but quickly. 7373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel n is an int, greater than or equal to 0. 7393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Example use: If you have 2 threads and know that each will 7413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel consume no more than a million random numbers, create two Random 7423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel objects r1 and r2, then do 7433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r2.setstate(r1.getstate()) 7443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel r2.jumpahead(1000000) 7453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Then r1 and r2 will use guaranteed-disjoint segments of the full 7463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel period. 7473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 7483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not n >= 0: 7503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError("n must be >= 0") 7513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x, y, z = self._seed 7523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = int(x * pow(171, n, 30269)) % 30269 7533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel y = int(y * pow(172, n, 30307)) % 30307 7543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = int(z * pow(170, n, 30323)) % 30323 7553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self._seed = x, y, z 7563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def __whseed(self, x=0, y=0, z=0): 7583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Set the Wichmann-Hill seed from (x, y, z). 7593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel These must be integers in the range [0, 256). 7613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 7623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not type(x) == type(y) == type(z) == int: 7643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise TypeError('seeds must be integers') 7653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): 7663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError('seeds must be in range(0, 256)') 7673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if 0 == x == y == z: 7683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Initialize from current time 7693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel import time 7703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t = long(time.time() * 256) 7713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t = int((t&0xffffff) ^ (t>>24)) 7723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t, x = divmod(t, 256) 7733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t, y = divmod(t, 256) 7743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t, z = divmod(t, 256) 7753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel # Zero is a poor seed, so substitute 1 7763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self._seed = (x or 1, y or 1, z or 1) 7773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.gauss_next = None 7793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def whseed(self, a=None): 7813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Seed from hashable object's hash code. 7823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel None or no argument seeds from current time. It is not guaranteed 7843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel that objects with distinct hash codes lead to distinct internal 7853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel states. 7863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel This is obsolete, provided for compatibility with the seed routine 7883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel used prior to Python 2.1. Use the .seed() method instead. 7893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 7903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 7913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if a is None: 7923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.__whseed() 7933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return 7943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a = hash(a) 7953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, x = divmod(a, 256) 7963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, y = divmod(a, 256) 7973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel a, z = divmod(a, 256) 7983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = (x + a) % 256 or 1 7993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel y = (y + a) % 256 or 1 8003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel z = (z + a) % 256 or 1 8013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel self.__whseed(x, y, z) 8023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## --------------- Operating System Random Source ------------------ 8043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielclass SystemRandom(Random): 8063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Alternate random number generator using sources provided 8073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel by the operating system (such as /dev/urandom on Unix or 8083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel CryptGenRandom on Windows). 8093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel Not available on all systems (see os.urandom() for details). 8113257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """ 8123257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8133257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def random(self): 8143257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """Get the next random number in the range [0.0, 1.0).""" 8153257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF 8163257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8173257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def getrandbits(self, k): 8183257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel """getrandbits(k) -> x. Generates a long int with k random bits.""" 8193257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if k <= 0: 8203257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise ValueError('number of bits must be greater than zero') 8213257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel if k != int(k): 8223257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise TypeError('number of bits should be an integer') 8233257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel bytes = (k + 7) // 8 # bits / 8 and rounded up 8243257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = long(_hexlify(_urandom(bytes)), 16) 8253257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return x >> (bytes * 8 - k) # trim excess bits 8263257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8273257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def _stub(self, *args, **kwds): 8283257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "Stub method. Not used for a system random number generator." 8293257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel return None 8303257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel seed = jumpahead = _stub 8313257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8323257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel def _notimplemented(self, *args, **kwds): 8333257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel "Method should not be called for a system random number generator." 8343257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel raise NotImplementedError('System entropy source does not have state.') 8353257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel getstate = setstate = _notimplemented 8363257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8373257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel## -------------------- test program -------------------- 8383257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8393257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanieldef _test_generator(n, func, args): 8403257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel import time 8413257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel print n, 'times', func.__name__ 8423257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel total = 0.0 8433257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel sqsum = 0.0 8443257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel smallest = 1e10 8453257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel largest = -1e10 8463257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t0 = time.time() 8473257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel for i in range(n): 8483257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel x = func(*args) 8493257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel total += x 8503257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel sqsum = sqsum + x*x 8513257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel smallest = min(x, smallest) 8523257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel largest = max(x, largest) 8533257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel t1 = time.time() 8543257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel print round(t1-t0, 3), 'sec,', 8553257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel avg = total/n 8563257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel stddev = _sqrt(sqsum/n - avg*avg) 8573257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel print 'avg %g, stddev %g, min %g, max %g' % \ 8583257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel (avg, stddev, smallest, largest) 8593257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8603257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8613257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanieldef _test(N=2000): 8623257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, random, ()) 8633257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, normalvariate, (0.0, 1.0)) 8643257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, lognormvariate, (0.0, 1.0)) 8653257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, vonmisesvariate, (0.0, 1.0)) 8663257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (0.01, 1.0)) 8673257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (0.1, 1.0)) 8683257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (0.1, 2.0)) 8693257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (0.5, 1.0)) 8703257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (0.9, 1.0)) 8713257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (1.0, 1.0)) 8723257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (2.0, 1.0)) 8733257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (20.0, 1.0)) 8743257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gammavariate, (200.0, 1.0)) 8753257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, gauss, (0.0, 1.0)) 8763257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, betavariate, (3.0, 3.0)) 8773257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test_generator(N, triangular, (0.0, 1.0, 1.0/3.0)) 8783257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8793257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# Create one instance, seeded from current time, and export its methods 8803257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# as module-level functions. The functions share state across all uses 8813257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel#(both in the user's code and in the Python libraries), but that's fine 8823257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# for most programs and is easier for the casual user than making them 8833257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel# instantiate their own Random() instance. 8843257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 8853257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel_inst = Random() 8863257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielseed = _inst.seed 8873257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielrandom = _inst.random 8883257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanieluniform = _inst.uniform 8893257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanieltriangular = _inst.triangular 8903257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielrandint = _inst.randint 8913257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielchoice = _inst.choice 8923257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielrandrange = _inst.randrange 8933257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielsample = _inst.sample 8943257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielshuffle = _inst.shuffle 8953257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielnormalvariate = _inst.normalvariate 8963257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniellognormvariate = _inst.lognormvariate 8973257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielexpovariate = _inst.expovariate 8983257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielvonmisesvariate = _inst.vonmisesvariate 8993257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielgammavariate = _inst.gammavariate 9003257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielgauss = _inst.gauss 9013257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielbetavariate = _inst.betavariate 9023257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielparetovariate = _inst.paretovariate 9033257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielweibullvariate = _inst.weibullvariate 9043257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielgetstate = _inst.getstate 9053257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielsetstate = _inst.setstate 9063257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanieljumpahead = _inst.jumpahead 9073257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielgetrandbits = _inst.getrandbits 9083257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel 9093257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDanielif __name__ == '__main__': 9103257aa99321d745773a6bd1bd4ce7f6fafe74411Daryl McDaniel _test() 911