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