Random variable generators

diff --git a/Lib/random.py b/Lib/random.py
new file mode 100644 (file)
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()