home *** CD-ROM | disk | FTP | other *** search
/ Chip 2003 January / Chip_2003-01_cd2.bin / convert / eJayMp3Pro / mp3pro_demo.exe / RANDOM.PY < prev    next >
Encoding:
Python Source  |  1998-08-01  |  8.9 KB  |  343 lines
  1. #    R A N D O M   V A R I A B L E   G E N E R A T O R S
  2. #
  3. #    distributions on the real line:
  4. #    ------------------------------
  5. #           normal (Gaussian)
  6. #           lognormal
  7. #           negative exponential
  8. #           gamma
  9. #           beta
  10. #
  11. #    distributions on the circle (angles 0 to 2pi)
  12. #    ---------------------------------------------
  13. #           circular uniform
  14. #           von Mises
  15.  
  16. # Translated from anonymously contributed C/C++ source.
  17.  
  18. # Multi-threading note: the random number generator used here is not
  19. # thread-safe; it is possible that two calls return the same random
  20. # value.  See whrandom.py for more info.
  21.  
  22. import whrandom
  23. from whrandom import random, uniform, randint, choice, randrange # For export!
  24. from math import log, exp, pi, e, sqrt, acos, cos, sin
  25.  
  26. # Interfaces to replace remaining needs for importing whrandom
  27. # XXX TO DO: make the distribution functions below into methods.
  28.  
  29. def makeseed(a=None):
  30.     """Turn a hashable value into three seed values for whrandom.seed().
  31.  
  32.     None or no argument returns (0, 0, 0), to seed from current time.
  33.  
  34.     """
  35.     if a is None:
  36.         return (0, 0, 0)
  37.     a = hash(a)
  38.     a, x = divmod(a, 256)
  39.     a, y = divmod(a, 256)
  40.     a, z = divmod(a, 256)
  41.     x = (x + a) % 256 or 1
  42.     y = (y + a) % 256 or 1
  43.     z = (z + a) % 256 or 1
  44.     return (x, y, z)
  45.  
  46. def seed(a=None):
  47.     """Seed the default generator from any hashable value.
  48.  
  49.     None or no argument returns (0, 0, 0) to seed from current time.
  50.  
  51.     """
  52.     x, y, z = makeseed(a)
  53.     whrandom.seed(x, y, z)
  54.  
  55. class generator(whrandom.whrandom):
  56.     """Random generator class."""
  57.  
  58.     def __init__(self, a=None):
  59.         """Constructor.  Seed from current time or hashable value."""
  60.         self.seed(a)
  61.  
  62.     def seed(self, a=None):
  63.         """Seed the generator from current time or hashable value."""
  64.         x, y, z = makeseed(a)
  65.         whrandom.whrandom.seed(self, x, y, z)
  66.  
  67. def new_generator(a=None):
  68.     """Return a new random generator instance."""
  69.     return generator(a)
  70.  
  71. # Housekeeping function to verify that magic constants have been
  72. # computed correctly
  73.  
  74. def verify(name, expected):
  75.     computed = eval(name)
  76.     if abs(computed - expected) > 1e-7:
  77.         raise ValueError, \
  78.   'computed value for %s deviates too much (computed %g, expected %g)' % \
  79.   (name, computed, expected)
  80.  
  81. # -------------------- normal distribution --------------------
  82.  
  83. NV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0)
  84. verify('NV_MAGICCONST', 1.71552776992141)
  85. def normalvariate(mu, sigma):
  86.     # mu = mean, sigma = standard deviation
  87.  
  88.     # Uses Kinderman and Monahan method. Reference: Kinderman,
  89.     # A.J. and Monahan, J.F., "Computer generation of random
  90.     # variables using the ratio of uniform deviates", ACM Trans
  91.     # Math Software, 3, (1977), pp257-260.
  92.  
  93.     while 1:
  94.         u1 = random()
  95.         u2 = random()
  96.         z = NV_MAGICCONST*(u1-0.5)/u2
  97.         zz = z*z/4.0
  98.         if zz <= -log(u2):
  99.             break
  100.     return mu+z*sigma
  101.  
  102. # -------------------- lognormal distribution --------------------
  103.  
  104. def lognormvariate(mu, sigma):
  105.     return exp(normalvariate(mu, sigma))
  106.  
  107. # -------------------- circular uniform --------------------
  108.  
  109. def cunifvariate(mean, arc):
  110.     # mean: mean angle (in radians between 0 and pi)
  111.     # arc:  range of distribution (in radians between 0 and pi)
  112.  
  113.     return (mean + arc * (random() - 0.5)) % pi
  114.  
  115. # -------------------- exponential distribution --------------------
  116.  
  117. def expovariate(lambd):
  118.     # lambd: rate lambd = 1/mean
  119.     # ('lambda' is a Python reserved word)
  120.  
  121.     u = random()
  122.     while u <= 1e-7:
  123.         u = random()
  124.     return -log(u)/lambd
  125.  
  126. # -------------------- von Mises distribution --------------------
  127.  
  128. TWOPI = 2.0*pi
  129. verify('TWOPI', 6.28318530718)
  130.  
  131. def vonmisesvariate(mu, kappa):
  132.     # mu:    mean angle (in radians between 0 and 2*pi)
  133.     # kappa: concentration parameter kappa (>= 0)
  134.     # if kappa = 0 generate uniform random angle
  135.  
  136.     # Based upon an algorithm published in: Fisher, N.I.,
  137.     # "Statistical Analysis of Circular Data", Cambridge
  138.     # University Press, 1993.
  139.  
  140.     # Thanks to Magnus Kessler for a correction to the
  141.     # implementation of step 4.
  142.  
  143.     if kappa <= 1e-6:
  144.         return TWOPI * random()
  145.  
  146.     a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa)
  147.     b = (a - sqrt(2.0 * a))/(2.0 * kappa)
  148.     r = (1.0 + b * b)/(2.0 * b)
  149.  
  150.     while 1:
  151.         u1 = random()
  152.  
  153.         z = cos(pi * u1)
  154.         f = (1.0 + r * z)/(r + z)
  155.         c = kappa * (r - f)
  156.  
  157.         u2 = random()
  158.  
  159.         if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)):
  160.             break
  161.  
  162.     u3 = random()
  163.     if u3 > 0.5:
  164.         theta = (mu % TWOPI) + acos(f)
  165.     else:
  166.         theta = (mu % TWOPI) - acos(f)
  167.  
  168.     return theta
  169.  
  170. # -------------------- gamma distribution --------------------
  171.  
  172. LOG4 = log(4.0)
  173. verify('LOG4', 1.38629436111989)
  174.  
  175. def gammavariate(alpha, beta):
  176.         # beta times standard gamma
  177.     ainv = sqrt(2.0 * alpha - 1.0)
  178.     return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv)
  179.  
  180. SG_MAGICCONST = 1.0 + log(4.5)
  181. verify('SG_MAGICCONST', 2.50407739677627)
  182.  
  183. def stdgamma(alpha, ainv, bbb, ccc):
  184.     # ainv = sqrt(2 * alpha - 1)
  185.     # bbb = alpha - log(4)
  186.     # ccc = alpha + ainv
  187.  
  188.     if alpha <= 0.0:
  189.         raise ValueError, 'stdgamma: alpha must be > 0.0'
  190.  
  191.     if alpha > 1.0:
  192.  
  193.         # Uses R.C.H. Cheng, "The generation of Gamma
  194.         # variables with non-integral shape parameters",
  195.         # Applied Statistics, (1977), 26, No. 1, p71-74
  196.  
  197.         while 1:
  198.             u1 = random()
  199.             u2 = random()
  200.             v = log(u1/(1.0-u1))/ainv
  201.             x = alpha*exp(v)
  202.             z = u1*u1*u2
  203.             r = bbb+ccc*v-x
  204.             if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z):
  205.                 return x
  206.  
  207.     elif alpha == 1.0:
  208.         # expovariate(1)
  209.         u = random()
  210.         while u <= 1e-7:
  211.             u = random()
  212.         return -log(u)
  213.  
  214.     else:    # alpha is between 0 and 1 (exclusive)
  215.  
  216.         # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
  217.  
  218.         while 1:
  219.             u = random()
  220.             b = (e + alpha)/e
  221.             p = b*u
  222.             if p <= 1.0:
  223.                 x = pow(p, 1.0/alpha)
  224.             else:
  225.                 # p > 1
  226.                 x = -log((b-p)/alpha)
  227.             u1 = random()
  228.             if not (((p <= 1.0) and (u1 > exp(-x))) or
  229.                   ((p > 1)  and  (u1 > pow(x, alpha - 1.0)))):
  230.                 break
  231.         return x
  232.  
  233.  
  234. # -------------------- Gauss (faster alternative) --------------------
  235.  
  236. gauss_next = None
  237. def gauss(mu, sigma):
  238.  
  239.     # When x and y are two variables from [0, 1), uniformly
  240.     # distributed, then
  241.     #
  242.     #    cos(2*pi*x)*sqrt(-2*log(1-y))
  243.     #    sin(2*pi*x)*sqrt(-2*log(1-y))
  244.     #
  245.     # are two *independent* variables with normal distribution
  246.     # (mu = 0, sigma = 1).
  247.     # (Lambert Meertens)
  248.     # (corrected version; bug discovered by Mike Miller, fixed by LM)
  249.  
  250.     # Multithreading note: When two threads call this function
  251.     # simultaneously, it is possible that they will receive the
  252.     # same return value.  The window is very small though.  To
  253.     # avoid this, you have to use a lock around all calls.  (I
  254.     # didn't want to slow this down in the serial case by using a
  255.     # lock here.)
  256.  
  257.     global gauss_next
  258.  
  259.     z = gauss_next
  260.     gauss_next = None
  261.     if z is None:
  262.         x2pi = random() * TWOPI
  263.         g2rad = sqrt(-2.0 * log(1.0 - random()))
  264.         z = cos(x2pi) * g2rad
  265.         gauss_next = sin(x2pi) * g2rad
  266.  
  267.     return mu + z*sigma
  268.  
  269. # -------------------- beta --------------------
  270.  
  271. def betavariate(alpha, beta):
  272.  
  273.     # Discrete Event Simulation in C, pp 87-88.
  274.  
  275.     y = expovariate(alpha)
  276.     z = expovariate(1.0/beta)
  277.     return z/(y+z)
  278.  
  279. # -------------------- Pareto --------------------
  280.  
  281. def paretovariate(alpha):
  282.     # Jain, pg. 495
  283.  
  284.     u = random()
  285.     return 1.0 / pow(u, 1.0/alpha)
  286.  
  287. # -------------------- Weibull --------------------
  288.  
  289. def weibullvariate(alpha, beta):
  290.     # Jain, pg. 499; bug fix courtesy Bill Arms
  291.  
  292.     u = random()
  293.     return alpha * pow(-log(u), 1.0/beta)
  294.  
  295. # -------------------- test program --------------------
  296.  
  297. def test(N = 200):
  298.     print 'TWOPI         =', TWOPI
  299.     print 'LOG4          =', LOG4
  300.     print 'NV_MAGICCONST =', NV_MAGICCONST
  301.     print 'SG_MAGICCONST =', SG_MAGICCONST
  302.     test_generator(N, 'random()')
  303.     test_generator(N, 'normalvariate(0.0, 1.0)')
  304.     test_generator(N, 'lognormvariate(0.0, 1.0)')
  305.     test_generator(N, 'cunifvariate(0.0, 1.0)')
  306.     test_generator(N, 'expovariate(1.0)')
  307.     test_generator(N, 'vonmisesvariate(0.0, 1.0)')
  308.     test_generator(N, 'gammavariate(0.5, 1.0)')
  309.     test_generator(N, 'gammavariate(0.9, 1.0)')
  310.     test_generator(N, 'gammavariate(1.0, 1.0)')
  311.     test_generator(N, 'gammavariate(2.0, 1.0)')
  312.     test_generator(N, 'gammavariate(20.0, 1.0)')
  313.     test_generator(N, 'gammavariate(200.0, 1.0)')
  314.     test_generator(N, 'gauss(0.0, 1.0)')
  315.     test_generator(N, 'betavariate(3.0, 3.0)')
  316.     test_generator(N, 'paretovariate(1.0)')
  317.     test_generator(N, 'weibullvariate(1.0, 1.0)')
  318.  
  319. def test_generator(n, funccall):
  320.     import time
  321.     print n, 'times', funccall
  322.     code = compile(funccall, funccall, 'eval')
  323.     sum = 0.0
  324.     sqsum = 0.0
  325.     smallest = 1e10
  326.     largest = -1e10
  327.     t0 = time.time()
  328.     for i in range(n):
  329.         x = eval(code)
  330.         sum = sum + x
  331.         sqsum = sqsum + x*x
  332.         smallest = min(x, smallest)
  333.         largest = max(x, largest)
  334.     t1 = time.time()
  335.     print round(t1-t0, 3), 'sec,', 
  336.     avg = sum/n
  337.     stddev = sqrt(sqsum/n - avg*avg)
  338.     print 'avg %g, stddev %g, min %g, max %g' % \
  339.           (avg, stddev, smallest, largest)
  340.  
  341. if __name__ == '__main__':
  342.     test()
  343.