random.py revision 6c395ba31609eeffce2428280cc5d95e4fb8058a
1ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# R A N D O M V A R I A B L E G E N E R A T O R S 2ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# 3ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# distributions on the real line: 4ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# ------------------------------ 5ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# normal (Gaussian) 6ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# lognormal 7ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# negative exponential 8ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# gamma 995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum# beta 10ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# 11ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# distributions on the circle (angles 0 to 2pi) 12ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# --------------------------------------------- 13ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# circular uniform 14ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# von Mises 15ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 16ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# Translated from anonymously contributed C/C++ source. 17ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 18d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum# Multi-threading note: the random number generator used here is not 19d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum# thread-safe; it is possible that two calls return the same random 20d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum# value. See whrandom.py for more info. 21d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 2233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumimport whrandom 233357561476d3b5570ca9e4936baee444e504c39fGuido van Rossumfrom whrandom import random, uniform, randint, choice, randrange # For export! 2495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumfrom math import log, exp, pi, e, sqrt, acos, cos, sin 25ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 2633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum# Interfaces to replace remaining needs for importing whrandom 2733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum# XXX TO DO: make the distribution functions below into methods. 2833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 2933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef makeseed(a=None): 3033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Turn a hashable value into three seed values for whrandom.seed(). 3133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 3233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum None or no argument returns (0, 0, 0), to seed from current time. 3333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 3433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """ 3533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum if a is None: 3633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return (0, 0, 0) 3733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a = hash(a) 3833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, x = divmod(a, 256) 3933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, y = divmod(a, 256) 4033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, z = divmod(a, 256) 4133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x = (x + a) % 256 or 1 4233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum y = (y + a) % 256 or 1 4333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum z = (z + a) % 256 or 1 4433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return (x, y, z) 4533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 4633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef seed(a=None): 4733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Seed the default generator from any hashable value. 4833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 4933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum None or no argument returns (0, 0, 0) to seed from current time. 5033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """ 5233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x, y, z = makeseed(a) 5333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum whrandom.seed(x, y, z) 5433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumclass generator(whrandom.whrandom): 5633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Random generator class.""" 5733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum def __init__(self, a=None): 5933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Constructor. Seed from current time or hashable value.""" 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum self.seed(a) 6133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 6233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum def seed(self, a=None): 6333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Seed the generator from current time or hashable value.""" 6433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x, y, z = makeseed(a) 6533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum whrandom.whrandom.seed(self, x, y, z) 6633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 6733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef new_generator(a=None): 6833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Return a new random generator instance.""" 6933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return generator(a) 7033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 71ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# Housekeeping function to verify that magic constants have been 72ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# computed correctly 73ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 74ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef verify(name, expected): 75ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum computed = eval(name) 76ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if abs(computed - expected) > 1e-7: 77ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum raise ValueError, \ 78ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 'computed value for %s deviates too much (computed %g, expected %g)' % \ 79ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum (name, computed, expected) 80ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 81ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- normal distribution -------------------- 82ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 83cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumNV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0) 84ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('NV_MAGICCONST', 1.71552776992141) 85ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef normalvariate(mu, sigma): 86ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # mu = mean, sigma = standard deviation 87ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 88ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses Kinderman and Monahan method. Reference: Kinderman, 89ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # A.J. and Monahan, J.F., "Computer generation of random 90ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # variables using the ratio of uniform deviates", ACM Trans 91ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Math Software, 3, (1977), pp257-260. 92ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 93ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 94ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 95ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 96ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = NV_MAGICCONST*(u1-0.5)/u2 97cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum zz = z*z/4.0 98ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if zz <= -log(u2): 99ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 100ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return mu+z*sigma 101ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 102ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- lognormal distribution -------------------- 103ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 104ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef lognormvariate(mu, sigma): 105ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return exp(normalvariate(mu, sigma)) 106ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 107ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- circular uniform -------------------- 108ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 109ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef cunifvariate(mean, arc): 110ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # mean: mean angle (in radians between 0 and pi) 111ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # arc: range of distribution (in radians between 0 and pi) 112ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 113ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return (mean + arc * (random() - 0.5)) % pi 114ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 115ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- exponential distribution -------------------- 116ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 117ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef expovariate(lambd): 118ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # lambd: rate lambd = 1/mean 119ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ('lambda' is a Python reserved word) 120ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 121ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 122ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while u <= 1e-7: 123ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 124ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return -log(u)/lambd 125ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 126ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- von Mises distribution -------------------- 127ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 128cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumTWOPI = 2.0*pi 129ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('TWOPI', 6.28318530718) 130ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 131ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef vonmisesvariate(mu, kappa): 1325810297052003f28788f6790ac799fe8e5373494Guido van Rossum # mu: mean angle (in radians between 0 and 2*pi) 133ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # kappa: concentration parameter kappa (>= 0) 134ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # if kappa = 0 generate uniform random angle 1355810297052003f28788f6790ac799fe8e5373494Guido van Rossum 1365810297052003f28788f6790ac799fe8e5373494Guido van Rossum # Based upon an algorithm published in: Fisher, N.I., 1375810297052003f28788f6790ac799fe8e5373494Guido van Rossum # "Statistical Analysis of Circular Data", Cambridge 1385810297052003f28788f6790ac799fe8e5373494Guido van Rossum # University Press, 1993. 1395810297052003f28788f6790ac799fe8e5373494Guido van Rossum 1405810297052003f28788f6790ac799fe8e5373494Guido van Rossum # Thanks to Magnus Kessler for a correction to the 1415810297052003f28788f6790ac799fe8e5373494Guido van Rossum # implementation of step 4. 1425810297052003f28788f6790ac799fe8e5373494Guido van Rossum 143ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if kappa <= 1e-6: 144ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return TWOPI * random() 145ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 146cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa) 147cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum b = (a - sqrt(2.0 * a))/(2.0 * kappa) 148cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum r = (1.0 + b * b)/(2.0 * b) 149ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 150ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 151ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 152ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 153ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = cos(pi * u1) 154cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum f = (1.0 + r * z)/(r + z) 155ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum c = kappa * (r - f) 156ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 157ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 158ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 159ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)): 160ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 161ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 162ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u3 = random() 163ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if u3 > 0.5: 1645810297052003f28788f6790ac799fe8e5373494Guido van Rossum theta = (mu % TWOPI) + acos(f) 165ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: 1665810297052003f28788f6790ac799fe8e5373494Guido van Rossum theta = (mu % TWOPI) - acos(f) 167ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1685810297052003f28788f6790ac799fe8e5373494Guido van Rossum return theta 169ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 170ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- gamma distribution -------------------- 171ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 172cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumLOG4 = log(4.0) 173ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('LOG4', 1.38629436111989) 174ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 175ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef gammavariate(alpha, beta): 176ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # beta times standard gamma 177cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum ainv = sqrt(2.0 * alpha - 1.0) 178ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) 179ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 180cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumSG_MAGICCONST = 1.0 + log(4.5) 181ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('SG_MAGICCONST', 2.50407739677627) 182ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 183ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef stdgamma(alpha, ainv, bbb, ccc): 184ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ainv = sqrt(2 * alpha - 1) 185ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # bbb = alpha - log(4) 186ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ccc = alpha + ainv 187ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 188ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if alpha <= 0.0: 189ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum raise ValueError, 'stdgamma: alpha must be > 0.0' 190ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 191ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if alpha > 1.0: 192ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 193ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses R.C.H. Cheng, "The generation of Gamma 194ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # variables with non-integral shape parameters", 195ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Applied Statistics, (1977), 26, No. 1, p71-74 196ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 197ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 198ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 199ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 200cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum v = log(u1/(1.0-u1))/ainv 201ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = alpha*exp(v) 202ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = u1*u1*u2 203ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum r = bbb+ccc*v-x 204cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z): 205ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return x 206ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 207ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum elif alpha == 1.0: 208ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # expovariate(1) 209ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 210ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while u <= 1e-7: 211ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 212ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return -log(u) 213ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 214ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: # alpha is between 0 and 1 (exclusive) 215ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 216ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 217ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 218ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 219ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 220ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum b = (e + alpha)/e 221ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum p = b*u 222ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if p <= 1.0: 223ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = pow(p, 1.0/alpha) 224ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: 225ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # p > 1 226ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = -log((b-p)/alpha) 227ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 228ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if not (((p <= 1.0) and (u1 > exp(-x))) or 229ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 230ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 231ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return x 232ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 23395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 23495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum# -------------------- Gauss (faster alternative) -------------------- 23595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 23695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumgauss_next = None 23795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumdef gauss(mu, sigma): 238cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 239cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # When x and y are two variables from [0, 1), uniformly 240cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # distributed, then 241cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # 24272c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # cos(2*pi*x)*sqrt(-2*log(1-y)) 24372c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # sin(2*pi*x)*sqrt(-2*log(1-y)) 244cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # 245cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # are two *independent* variables with normal distribution 246cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # (mu = 0, sigma = 1). 247cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # (Lambert Meertens) 24872c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # (corrected version; bug discovered by Mike Miller, fixed by LM) 249cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 250d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # Multithreading note: When two threads call this function 251d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # simultaneously, it is possible that they will receive the 252d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # same return value. The window is very small though. To 253d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # avoid this, you have to use a lock around all calls. (I 254d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # didn't want to slow this down in the serial case by using a 255d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum # lock here.) 256d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum 25795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum global gauss_next 258cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 259d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum z = gauss_next 260d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum gauss_next = None 261d03e1197cb5052e3f758794e2a7aecf9f5ca5f46Guido van Rossum if z is None: 26295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum x2pi = random() * TWOPI 26372c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum g2rad = sqrt(-2.0 * log(1.0 - random())) 26472c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum z = cos(x2pi) * g2rad 26572c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum gauss_next = sin(x2pi) * g2rad 266cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 26795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum return mu + z*sigma 26895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 26995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum# -------------------- beta -------------------- 27095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 27195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumdef betavariate(alpha, beta): 272cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 273cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # Discrete Event Simulation in C, pp 87-88. 274cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 27595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum y = expovariate(alpha) 27695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum z = expovariate(1.0/beta) 27795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum return z/(y+z) 27895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 2795bdea89c892b1a10281eeae7b60da7f9a0c15ec4Guido van Rossum# -------------------- Pareto -------------------- 280cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 281cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossumdef paretovariate(alpha): 282cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum # Jain, pg. 495 283cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 284cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum u = random() 285cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum return 1.0 / pow(u, 1.0/alpha) 286cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 2875bdea89c892b1a10281eeae7b60da7f9a0c15ec4Guido van Rossum# -------------------- Weibull -------------------- 288cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 289cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossumdef weibullvariate(alpha, beta): 290cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum # Jain, pg. 499; bug fix courtesy Bill Arms 291cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 292cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum u = random() 293cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum return alpha * pow(-log(u), 1.0/beta) 294cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 2956c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum# -------------------- shuffle -------------------- 2966c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum# Not quite a random distribution, but a standard algorithm. 2976c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum# This implementation due to Tim Peters. 2986c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 2996c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossumdef shuffle(x, random=random, int=int): 3006c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum """x, random=random.random -> shuffle list x in place; return None. 3016c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 3026c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum Optional arg random is a 0-argument function returning a random 3036c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum float in [0.0, 1.0); by default, the standard random.random. 3046c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 3056c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum Note that for even rather small len(x), the total number of 3066c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum permutations of x is larger than the period of most random number 3076c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum generators; this implies that "most" permutations of a long 3086c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum sequence can never be generated. 3096c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum """ 3106c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 3116c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum for i in xrange(len(x)-1, 0, -1): 3126c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum # pick an element in x[:i+1] with which to exchange x[i] 3136c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum j = int(random() * (i+1)) 3146c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum x[i], x[j] = x[j], x[i] 3156c395ba31609eeffce2428280cc5d95e4fb8058aGuido van Rossum 316ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- test program -------------------- 317ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 3182922c6dabbd9f8e49975ff3d972644c7212882e9Guido van Rossumdef test(N = 200): 319ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'TWOPI =', TWOPI 320ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'LOG4 =', LOG4 321ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'NV_MAGICCONST =', NV_MAGICCONST 322ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'SG_MAGICCONST =', SG_MAGICCONST 323ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'random()') 324ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'normalvariate(0.0, 1.0)') 325ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'lognormvariate(0.0, 1.0)') 326ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'cunifvariate(0.0, 1.0)') 327ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'expovariate(1.0)') 328ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'vonmisesvariate(0.0, 1.0)') 329ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(0.5, 1.0)') 330ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(0.9, 1.0)') 331ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(1.0, 1.0)') 332ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(2.0, 1.0)') 333ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(20.0, 1.0)') 334ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(200.0, 1.0)') 33595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum test_generator(N, 'gauss(0.0, 1.0)') 33695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum test_generator(N, 'betavariate(3.0, 3.0)') 337cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum test_generator(N, 'paretovariate(1.0)') 338cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum test_generator(N, 'weibullvariate(1.0, 1.0)') 339ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 340ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef test_generator(n, funccall): 34195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum import time 34295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print n, 'times', funccall 343ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum code = compile(funccall, funccall, 'eval') 344ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sum = 0.0 345ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sqsum = 0.0 34695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum smallest = 1e10 347cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum largest = -1e10 34895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum t0 = time.time() 349ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum for i in range(n): 350ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = eval(code) 351ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sum = sum + x 352ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sqsum = sqsum + x*x 35395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum smallest = min(x, smallest) 35495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum largest = max(x, largest) 35595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum t1 = time.time() 35695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print round(t1-t0, 3), 'sec,', 357ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum avg = sum/n 358ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum stddev = sqrt(sqsum/n - avg*avg) 35995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print 'avg %g, stddev %g, min %g, max %g' % \ 36095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum (avg, stddev, smallest, largest) 361ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 362ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 363ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test() 364