diff options
Diffstat (limited to 'Lib/random.py')
| -rw-r--r-- | Lib/random.py | 212 |
1 files changed, 212 insertions, 0 deletions
diff --git a/Lib/random.py b/Lib/random.py new file mode 100644 index 0000000..51ecb32 --- /dev/null +++ b/Lib/random.py @@ -0,0 +1,212 @@ +# R A N D O M V A R I A B L E G E N E R A T O R S +# +# distributions on the real line: +# ------------------------------ +# normal (Gaussian) +# lognormal +# negative exponential +# gamma +# +# distributions on the circle (angles 0 to 2pi) +# --------------------------------------------- +# circular uniform +# von Mises + +# Translated from anonymously contributed C/C++ source. + +from whrandom import random, uniform, randint, choice # Also for export! +from math import log, exp, pi, e, sqrt, acos, cos + +# Housekeeping function to verify that magic constants have been +# computed correctly + +def verify(name, expected): + computed = eval(name) + if abs(computed - expected) > 1e-7: + raise ValueError, \ + 'computed value for %s deviates too much (computed %g, expected %g)' % \ + (name, computed, expected) + +# -------------------- normal distribution -------------------- + +NV_MAGICCONST = 4*exp(-0.5)/sqrt(2) +verify('NV_MAGICCONST', 1.71552776992141) +def normalvariate(mu, sigma): + # mu = mean, sigma = standard deviation + + # Uses Kinderman and Monahan method. Reference: Kinderman, + # A.J. and Monahan, J.F., "Computer generation of random + # variables using the ratio of uniform deviates", ACM Trans + # Math Software, 3, (1977), pp257-260. + + while 1: + u1 = random() + u2 = random() + z = NV_MAGICCONST*(u1-0.5)/u2 + zz = z*z/4 + if zz <= -log(u2): + break + return mu+z*sigma + +# -------------------- lognormal distribution -------------------- + +def lognormvariate(mu, sigma): + return exp(normalvariate(mu, sigma)) + +# -------------------- circular uniform -------------------- + +def cunifvariate(mean, arc): + # mean: mean angle (in radians between 0 and pi) + # arc: range of distribution (in radians between 0 and pi) + + return (mean + arc * (random() - 0.5)) % pi + +# -------------------- exponential distribution -------------------- + +def expovariate(lambd): + # lambd: rate lambd = 1/mean + # ('lambda' is a Python reserved word) + + u = random() + while u <= 1e-7: + u = random() + return -log(u)/lambd + +# -------------------- von Mises distribution -------------------- + +TWOPI = 2*pi +verify('TWOPI', 6.28318530718) + +def vonmisesvariate(mu, kappa): + # mu: mean angle (in radians between 0 and 180 degrees) + # kappa: concentration parameter kappa (>= 0) + + # if kappa = 0 generate uniform random angle + if kappa <= 1e-6: + return TWOPI * random() + + a = 1.0 + sqrt(1 + 4 * kappa * kappa) + b = (a - sqrt(2 * a))/(2 * kappa) + r = (1 + b * b)/(2 * b) + + while 1: + u1 = random() + + z = cos(pi * u1) + f = (1 + r * z)/(r + z) + c = kappa * (r - f) + + u2 = random() + + if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)): + break + + u3 = random() + if u3 > 0.5: + theta = mu + 0.5*acos(f) + else: + theta = mu - 0.5*acos(f) + + return theta % pi + +# -------------------- gamma distribution -------------------- + +LOG4 = log(4) +verify('LOG4', 1.38629436111989) + +def gammavariate(alpha, beta): + # beta times standard gamma + ainv = sqrt(2 * alpha - 1) + return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) + +SG_MAGICCONST = 1+log(4.5) +verify('SG_MAGICCONST', 2.50407739677627) + +def stdgamma(alpha, ainv, bbb, ccc): + # ainv = sqrt(2 * alpha - 1) + # bbb = alpha - log(4) + # ccc = alpha + ainv + + if alpha <= 0.0: + raise ValueError, 'stdgamma: alpha must be > 0.0' + + if alpha > 1.0: + + # Uses R.C.H. Cheng, "The generation of Gamma + # variables with non-integral shape parameters", + # Applied Statistics, (1977), 26, No. 1, p71-74 + + while 1: + u1 = random() + u2 = random() + v = log(u1/(1-u1))/ainv + x = alpha*exp(v) + z = u1*u1*u2 + r = bbb+ccc*v-x + if r + SG_MAGICCONST - 4.5*z >= 0 or r >= log(z): + return x + + elif alpha == 1.0: + # expovariate(1) + u = random() + while u <= 1e-7: + u = random() + return -log(u) + + else: # alpha is between 0 and 1 (exclusive) + + # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle + + while 1: + u = random() + b = (e + alpha)/e + p = b*u + if p <= 1.0: + x = pow(p, 1.0/alpha) + else: + # p > 1 + x = -log((b-p)/alpha) + u1 = random() + if not (((p <= 1.0) and (u1 > exp(-x))) or + ((p > 1) and (u1 > pow(x, alpha - 1.0)))): + break + return x + +# -------------------- test program -------------------- + +def test(): + print 'TWOPI =', TWOPI + print 'LOG4 =', LOG4 + print 'NV_MAGICCONST =', NV_MAGICCONST + print 'SG_MAGICCONST =', SG_MAGICCONST + N = 100 + test_generator(N, 'random()') + test_generator(N, 'normalvariate(0.0, 1.0)') + test_generator(N, 'lognormvariate(0.0, 1.0)') + test_generator(N, 'cunifvariate(0.0, 1.0)') + test_generator(N, 'expovariate(1.0)') + test_generator(N, 'vonmisesvariate(0.0, 1.0)') + test_generator(N, 'gammavariate(0.5, 1.0)') + test_generator(N, 'gammavariate(0.9, 1.0)') + test_generator(N, 'gammavariate(1.0, 1.0)') + test_generator(N, 'gammavariate(2.0, 1.0)') + test_generator(N, 'gammavariate(20.0, 1.0)') + test_generator(N, 'gammavariate(200.0, 1.0)') + +def test_generator(n, funccall): + import sys + print '%d calls to %s:' % (n, funccall), + sys.stdout.flush() + code = compile(funccall, funccall, 'eval') + sum = 0.0 + sqsum = 0.0 + for i in range(n): + x = eval(code) + sum = sum + x + sqsum = sqsum + x*x + avg = sum/n + stddev = sqrt(sqsum/n - avg*avg) + print 'avg %g, stddev %g' % (avg, stddev) + +if __name__ == '__main__': + test() |