random.py revision 33d7f1a76c3544d2901492cfb6fc9db85f2dfbd6
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 1833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumimport whrandom 19ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumfrom whrandom import random, uniform, randint, choice # Also for export! 2095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumfrom math import log, exp, pi, e, sqrt, acos, cos, sin 21ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 2233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum# Interfaces to replace remaining needs for importing whrandom 2333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum# XXX TO DO: make the distribution functions below into methods. 2433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 2533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef makeseed(a=None): 2633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Turn a hashable value into three seed values for whrandom.seed(). 2733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 2833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum None or no argument returns (0, 0, 0), to seed from current time. 2933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 3033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """ 3133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum if a is None: 3233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return (0, 0, 0) 3333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a = hash(a) 3433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, x = divmod(a, 256) 3533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, y = divmod(a, 256) 3633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum a, z = divmod(a, 256) 3733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x = (x + a) % 256 or 1 3833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum y = (y + a) % 256 or 1 3933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum z = (z + a) % 256 or 1 4033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return (x, y, z) 4133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 4233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef seed(a=None): 4333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Seed the default generator from any hashable value. 4433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 4533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum None or no argument returns (0, 0, 0) to seed from current time. 4633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 4733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """ 4833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x, y, z = makeseed(a) 4933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum whrandom.seed(x, y, z) 5033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumclass generator(whrandom.whrandom): 5233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Random generator class.""" 5333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum def __init__(self, a=None): 5533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Constructor. Seed from current time or hashable value.""" 5633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum self.seed(a) 5733d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 5833d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum def seed(self, a=None): 5933d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Seed the generator from current time or hashable value.""" 6033d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum x, y, z = makeseed(a) 6133d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum whrandom.whrandom.seed(self, x, y, z) 6233d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 6333d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossumdef new_generator(a=None): 6433d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum """Return a new random generator instance.""" 6533d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum return generator(a) 6633d7f1a76c3544d2901492cfb6fc9db85f2dfbd6Guido van Rossum 67ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# Housekeeping function to verify that magic constants have been 68ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# computed correctly 69ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 70ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef verify(name, expected): 71ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum computed = eval(name) 72ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if abs(computed - expected) > 1e-7: 73ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum raise ValueError, \ 74ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 'computed value for %s deviates too much (computed %g, expected %g)' % \ 75ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum (name, computed, expected) 76ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 77ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- normal distribution -------------------- 78ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 79cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumNV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0) 80ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('NV_MAGICCONST', 1.71552776992141) 81ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef normalvariate(mu, sigma): 82ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # mu = mean, sigma = standard deviation 83ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 84ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses Kinderman and Monahan method. Reference: Kinderman, 85ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # A.J. and Monahan, J.F., "Computer generation of random 86ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # variables using the ratio of uniform deviates", ACM Trans 87ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Math Software, 3, (1977), pp257-260. 88ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 89ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 90ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 91ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 92ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = NV_MAGICCONST*(u1-0.5)/u2 93cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum zz = z*z/4.0 94ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if zz <= -log(u2): 95ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 96ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return mu+z*sigma 97ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 98ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- lognormal distribution -------------------- 99ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 100ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef lognormvariate(mu, sigma): 101ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return exp(normalvariate(mu, sigma)) 102ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 103ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- circular uniform -------------------- 104ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 105ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef cunifvariate(mean, arc): 106ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # mean: mean angle (in radians between 0 and pi) 107ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # arc: range of distribution (in radians between 0 and pi) 108ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 109ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return (mean + arc * (random() - 0.5)) % pi 110ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 111ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- exponential distribution -------------------- 112ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 113ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef expovariate(lambd): 114ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # lambd: rate lambd = 1/mean 115ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ('lambda' is a Python reserved word) 116ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 117ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 118ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while u <= 1e-7: 119ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 120ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return -log(u)/lambd 121ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 122ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- von Mises distribution -------------------- 123ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 124cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumTWOPI = 2.0*pi 125ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('TWOPI', 6.28318530718) 126ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 127ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef vonmisesvariate(mu, kappa): 1285810297052003f28788f6790ac799fe8e5373494Guido van Rossum # mu: mean angle (in radians between 0 and 2*pi) 129ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # kappa: concentration parameter kappa (>= 0) 130ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # if kappa = 0 generate uniform random angle 1315810297052003f28788f6790ac799fe8e5373494Guido van Rossum 1325810297052003f28788f6790ac799fe8e5373494Guido van Rossum # Based upon an algorithm published in: Fisher, N.I., 1335810297052003f28788f6790ac799fe8e5373494Guido van Rossum # "Statistical Analysis of Circular Data", Cambridge 1345810297052003f28788f6790ac799fe8e5373494Guido van Rossum # University Press, 1993. 1355810297052003f28788f6790ac799fe8e5373494Guido van Rossum 1365810297052003f28788f6790ac799fe8e5373494Guido van Rossum # Thanks to Magnus Kessler for a correction to the 1375810297052003f28788f6790ac799fe8e5373494Guido van Rossum # implementation of step 4. 1385810297052003f28788f6790ac799fe8e5373494Guido van Rossum 139ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if kappa <= 1e-6: 140ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return TWOPI * random() 141ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 142cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa) 143cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum b = (a - sqrt(2.0 * a))/(2.0 * kappa) 144cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum r = (1.0 + b * b)/(2.0 * b) 145ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 146ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 147ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 148ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 149ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = cos(pi * u1) 150cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum f = (1.0 + r * z)/(r + z) 151ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum c = kappa * (r - f) 152ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 153ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 154ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 155ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)): 156ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 157ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 158ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u3 = random() 159ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if u3 > 0.5: 1605810297052003f28788f6790ac799fe8e5373494Guido van Rossum theta = (mu % TWOPI) + acos(f) 161ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: 1625810297052003f28788f6790ac799fe8e5373494Guido van Rossum theta = (mu % TWOPI) - acos(f) 163ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 1645810297052003f28788f6790ac799fe8e5373494Guido van Rossum return theta 165ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 166ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- gamma distribution -------------------- 167ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 168cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumLOG4 = log(4.0) 169ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('LOG4', 1.38629436111989) 170ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 171ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef gammavariate(alpha, beta): 172ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # beta times standard gamma 173cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum ainv = sqrt(2.0 * alpha - 1.0) 174ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) 175ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 176cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van RossumSG_MAGICCONST = 1.0 + log(4.5) 177ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumverify('SG_MAGICCONST', 2.50407739677627) 178ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 179ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef stdgamma(alpha, ainv, bbb, ccc): 180ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ainv = sqrt(2 * alpha - 1) 181ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # bbb = alpha - log(4) 182ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # ccc = alpha + ainv 183ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 184ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if alpha <= 0.0: 185ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum raise ValueError, 'stdgamma: alpha must be > 0.0' 186ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 187ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if alpha > 1.0: 188ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 189ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses R.C.H. Cheng, "The generation of Gamma 190ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # variables with non-integral shape parameters", 191ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Applied Statistics, (1977), 26, No. 1, p71-74 192ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 193ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 194ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 195ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u2 = random() 196cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum v = log(u1/(1.0-u1))/ainv 197ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = alpha*exp(v) 198ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum z = u1*u1*u2 199ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum r = bbb+ccc*v-x 200cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z): 201ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return x 202ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 203ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum elif alpha == 1.0: 204ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # expovariate(1) 205ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 206ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while u <= 1e-7: 207ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 208ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return -log(u) 209ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 210ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: # alpha is between 0 and 1 (exclusive) 211ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 212ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 213ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 214ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum while 1: 215ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u = random() 216ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum b = (e + alpha)/e 217ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum p = b*u 218ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if p <= 1.0: 219ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = pow(p, 1.0/alpha) 220ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum else: 221ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum # p > 1 222ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = -log((b-p)/alpha) 223ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum u1 = random() 224ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum if not (((p <= 1.0) and (u1 > exp(-x))) or 225ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum ((p > 1) and (u1 > pow(x, alpha - 1.0)))): 226ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum break 227ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum return x 228ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 22995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 23095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum# -------------------- Gauss (faster alternative) -------------------- 23195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 23295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumgauss_next = None 23395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumdef gauss(mu, sigma): 234cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 235cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # When x and y are two variables from [0, 1), uniformly 236cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # distributed, then 237cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # 23872c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # cos(2*pi*x)*sqrt(-2*log(1-y)) 23972c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # sin(2*pi*x)*sqrt(-2*log(1-y)) 240cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # 241cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # are two *independent* variables with normal distribution 242cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # (mu = 0, sigma = 1). 243cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # (Lambert Meertens) 24472c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum # (corrected version; bug discovered by Mike Miller, fixed by LM) 245cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 24695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum global gauss_next 247cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 24895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum if gauss_next != None: 24995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum z = gauss_next 25095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum gauss_next = None 25195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum else: 25295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum x2pi = random() * TWOPI 25372c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum g2rad = sqrt(-2.0 * log(1.0 - random())) 25472c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum z = cos(x2pi) * g2rad 25572c2e1b56e35c7fc4a80e90b14541494426e3cd0Guido van Rossum gauss_next = sin(x2pi) * g2rad 256cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 25795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum return mu + z*sigma 25895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 25995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum# -------------------- beta -------------------- 26095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 26195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossumdef betavariate(alpha, beta): 262cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 263cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum # Discrete Event Simulation in C, pp 87-88. 264cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum 26595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum y = expovariate(alpha) 26695bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum z = expovariate(1.0/beta) 26795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum return z/(y+z) 26895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum 2695bdea89c892b1a10281eeae7b60da7f9a0c15ec4Guido van Rossum# -------------------- Pareto -------------------- 270cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 271cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossumdef paretovariate(alpha): 272cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum # Jain, pg. 495 273cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 274cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum u = random() 275cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum return 1.0 / pow(u, 1.0/alpha) 276cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 2775bdea89c892b1a10281eeae7b60da7f9a0c15ec4Guido van Rossum# -------------------- Weibull -------------------- 278cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 279cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossumdef weibullvariate(alpha, beta): 280cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum # Jain, pg. 499; bug fix courtesy Bill Arms 281cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 282cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum u = random() 283cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum return alpha * pow(-log(u), 1.0/beta) 284cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum 285ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum# -------------------- test program -------------------- 286ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 2872922c6dabbd9f8e49975ff3d972644c7212882e9Guido van Rossumdef test(N = 200): 288ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'TWOPI =', TWOPI 289ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'LOG4 =', LOG4 290ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'NV_MAGICCONST =', NV_MAGICCONST 291ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum print 'SG_MAGICCONST =', SG_MAGICCONST 292ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'random()') 293ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'normalvariate(0.0, 1.0)') 294ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'lognormvariate(0.0, 1.0)') 295ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'cunifvariate(0.0, 1.0)') 296ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'expovariate(1.0)') 297ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'vonmisesvariate(0.0, 1.0)') 298ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(0.5, 1.0)') 299ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(0.9, 1.0)') 300ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(1.0, 1.0)') 301ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(2.0, 1.0)') 302ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(20.0, 1.0)') 303ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test_generator(N, 'gammavariate(200.0, 1.0)') 30495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum test_generator(N, 'gauss(0.0, 1.0)') 30595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum test_generator(N, 'betavariate(3.0, 3.0)') 306cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum test_generator(N, 'paretovariate(1.0)') 307cf4559a62ec9316a3bb55a67c6fca81ec1ad0d18Guido van Rossum test_generator(N, 'weibullvariate(1.0, 1.0)') 308ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 309ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumdef test_generator(n, funccall): 31095bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum import time 31195bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print n, 'times', funccall 312ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum code = compile(funccall, funccall, 'eval') 313ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sum = 0.0 314ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sqsum = 0.0 31595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum smallest = 1e10 316cc32ac9704128c799170e1cd7bdbfb3a90da43c1Guido van Rossum largest = -1e10 31795bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum t0 = time.time() 318ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum for i in range(n): 319ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum x = eval(code) 320ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sum = sum + x 321ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum sqsum = sqsum + x*x 32295bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum smallest = min(x, smallest) 32395bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum largest = max(x, largest) 32495bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum t1 = time.time() 32595bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print round(t1-t0, 3), 'sec,', 326ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum avg = sum/n 327ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum stddev = sqrt(sqsum/n - avg*avg) 32895bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum print 'avg %g, stddev %g, min %g, max %g' % \ 32995bfcda3e0be2ace895e021296328a383eafb273Guido van Rossum (avg, stddev, smallest, largest) 330ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum 331ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossumif __name__ == '__main__': 332ff03b1ae5bba4d6712563efb7c77ace57dbe6788Guido van Rossum test() 333