blob: 8462061baae49037885f91f79e8e5689457db70c [file] [log] [blame]
Guido van Rossume7b146f2000-02-04 15:28:42 +00001"""Random variable generators.
Guido van Rossumff03b1a1994-03-09 12:55:02 +00002
Tim Petersd7b5e882001-01-25 03:36:26 +00003 integers
4 --------
5 uniform within range
6
7 sequences
8 ---------
9 pick random element
Raymond Hettingerf24eb352002-11-12 17:41:57 +000010 pick random sample
Tim Petersd7b5e882001-01-25 03:36:26 +000011 generate random permutation
12
Guido van Rossume7b146f2000-02-04 15:28:42 +000013 distributions on the real line:
14 ------------------------------
Tim Petersd7b5e882001-01-25 03:36:26 +000015 uniform
Guido van Rossume7b146f2000-02-04 15:28:42 +000016 normal (Gaussian)
17 lognormal
18 negative exponential
19 gamma
20 beta
Raymond Hettinger40f62172002-12-29 23:03:38 +000021 pareto
22 Weibull
Guido van Rossumff03b1a1994-03-09 12:55:02 +000023
Guido van Rossume7b146f2000-02-04 15:28:42 +000024 distributions on the circle (angles 0 to 2pi)
25 ---------------------------------------------
26 circular uniform
27 von Mises
28
Raymond Hettinger40f62172002-12-29 23:03:38 +000029General notes on the underlying Mersenne Twister core generator:
Guido van Rossume7b146f2000-02-04 15:28:42 +000030
Raymond Hettinger40f62172002-12-29 23:03:38 +000031* The period is 2**19937-1.
32* It is one of the most extensively tested generators in existence
33* Without a direct way to compute N steps forward, the
34 semantics of jumpahead(n) are weakened to simply jump
35 to another distant state and rely on the large period
36 to avoid overlapping sequences.
37* The random() method is implemented in C, executes in
38 a single Python step, and is, therefore, threadsafe.
Tim Peterse360d952001-01-26 10:00:39 +000039
Guido van Rossume7b146f2000-02-04 15:28:42 +000040"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000041
Tim Petersd7b5e882001-01-25 03:36:26 +000042from math import log as _log, exp as _exp, pi as _pi, e as _e
43from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Tim Peters9146f272002-08-16 03:41:39 +000044from math import floor as _floor
Guido van Rossumff03b1a1994-03-09 12:55:02 +000045
Raymond Hettingerf24eb352002-11-12 17:41:57 +000046__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000047 "randrange","shuffle","normalvariate","lognormvariate",
48 "cunifvariate","expovariate","vonmisesvariate","gammavariate",
49 "stdgamma","gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger40f62172002-12-29 23:03:38 +000050 "getstate","setstate","jumpahead"]
Tim Petersd7b5e882001-01-25 03:36:26 +000051
52NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000053TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000054LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000055SG_MAGICCONST = 1.0 + _log(4.5)
Tim Petersd7b5e882001-01-25 03:36:26 +000056
57# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000058# Adrian Baddeley. Adapted by Raymond Hettinger for use with
59# the Mersenne Twister core generator.
Tim Petersd7b5e882001-01-25 03:36:26 +000060
Raymond Hettinger40f62172002-12-29 23:03:38 +000061from _random import Random as CoreGenerator
62
63class Random(CoreGenerator):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000064 """Random number generator base class used by bound module functions.
65
66 Used to instantiate instances of Random to get generators that don't
67 share state. Especially useful for multi-threaded programs, creating
68 a different instance of Random for each thread, and using the jumpahead()
69 method to ensure that the generated sequences seen by each thread don't
70 overlap.
71
72 Class Random can also be subclassed if you want to use a different basic
73 generator of your own devising: in that case, override the following
74 methods: random(), seed(), getstate(), setstate() and jumpahead().
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +000075
Raymond Hettingerc32f0332002-05-23 19:44:49 +000076 """
Tim Petersd7b5e882001-01-25 03:36:26 +000077
Raymond Hettinger40f62172002-12-29 23:03:38 +000078 VERSION = 2 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000079
80 def __init__(self, x=None):
81 """Initialize an instance.
82
83 Optional argument x controls seeding, as for Random.seed().
84 """
85
86 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +000087 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +000088
Tim Peters0de88fc2001-02-01 04:59:18 +000089 def seed(self, a=None):
90 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +000091
Tim Peters0de88fc2001-02-01 04:59:18 +000092 None or no argument seeds from current time.
93
Tim Petersbcd725f2001-02-01 10:06:53 +000094 If a is not None or an int or long, hash(a) is used instead.
Tim Petersd7b5e882001-01-25 03:36:26 +000095 """
96
Raymond Hettinger40f62172002-12-29 23:03:38 +000097 CoreGenerator.seed(self, a)
Tim Peters46c04e12002-05-05 20:40:00 +000098 self.gauss_next = None
99
Tim Peterscd804102001-01-25 20:25:57 +0000100 def getstate(self):
101 """Return internal state; can be passed to setstate() later."""
Raymond Hettinger40f62172002-12-29 23:03:38 +0000102 return self.VERSION, CoreGenerator.getstate(self), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000103
104 def setstate(self, state):
105 """Restore internal state from object returned by getstate()."""
106 version = state[0]
Raymond Hettinger40f62172002-12-29 23:03:38 +0000107 if version == 2:
108 version, internalstate, self.gauss_next = state
109 CoreGenerator.setstate(self, internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000110 else:
111 raise ValueError("state with version %s passed to "
112 "Random.setstate() of version %s" %
113 (version, self.VERSION))
114
Tim Peterscd804102001-01-25 20:25:57 +0000115## ---- Methods below this point do not need to be overridden when
116## ---- subclassing for the purpose of using a different core generator.
117
118## -------------------- pickle support -------------------
119
120 def __getstate__(self): # for pickle
121 return self.getstate()
122
123 def __setstate__(self, state): # for pickle
124 self.setstate(state)
125
126## -------------------- integer methods -------------------
127
Tim Petersd7b5e882001-01-25 03:36:26 +0000128 def randrange(self, start, stop=None, step=1, int=int, default=None):
129 """Choose a random item from range(start, stop[, step]).
130
131 This fixes the problem with randint() which includes the
132 endpoint; in Python this is usually not what you want.
133 Do not supply the 'int' and 'default' arguments.
134 """
135
136 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000137 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000138 istart = int(start)
139 if istart != start:
140 raise ValueError, "non-integer arg 1 for randrange()"
141 if stop is default:
142 if istart > 0:
143 return int(self.random() * istart)
144 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000145
146 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000147 istop = int(stop)
148 if istop != stop:
149 raise ValueError, "non-integer stop for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000150 if step == 1 and istart < istop:
151 try:
152 return istart + int(self.random()*(istop - istart))
153 except OverflowError:
154 # This can happen if istop-istart > sys.maxint + 1, and
155 # multiplying by random() doesn't reduce it to something
156 # <= sys.maxint. We know that the overall result fits
157 # in an int, and can still do it correctly via math.floor().
158 # But that adds another function call, so for speed we
159 # avoided that whenever possible.
160 return int(istart + _floor(self.random()*(istop - istart)))
Tim Petersd7b5e882001-01-25 03:36:26 +0000161 if step == 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000162 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000163
164 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000165 istep = int(step)
166 if istep != step:
167 raise ValueError, "non-integer step for randrange()"
168 if istep > 0:
169 n = (istop - istart + istep - 1) / istep
170 elif istep < 0:
171 n = (istop - istart + istep + 1) / istep
172 else:
173 raise ValueError, "zero step for randrange()"
174
175 if n <= 0:
176 raise ValueError, "empty range for randrange()"
177 return istart + istep*int(self.random() * n)
178
179 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000180 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000181 """
182
183 return self.randrange(a, b+1)
184
Tim Peterscd804102001-01-25 20:25:57 +0000185## -------------------- sequence methods -------------------
186
Tim Petersd7b5e882001-01-25 03:36:26 +0000187 def choice(self, seq):
188 """Choose a random element from a non-empty sequence."""
189 return seq[int(self.random() * len(seq))]
190
191 def shuffle(self, x, random=None, int=int):
192 """x, random=random.random -> shuffle list x in place; return None.
193
194 Optional arg random is a 0-argument function returning a random
195 float in [0.0, 1.0); by default, the standard random.random.
196
197 Note that for even rather small len(x), the total number of
198 permutations of x is larger than the period of most random number
199 generators; this implies that "most" permutations of a long
200 sequence can never be generated.
201 """
202
203 if random is None:
204 random = self.random
205 for i in xrange(len(x)-1, 0, -1):
Tim Peterscd804102001-01-25 20:25:57 +0000206 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000207 j = int(random() * (i+1))
208 x[i], x[j] = x[j], x[i]
209
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000210 def sample(self, population, k, random=None, int=int):
211 """Chooses k unique random elements from a population sequence.
212
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000213 Returns a new list containing elements from the population while
214 leaving the original population unchanged. The resulting list is
215 in selection order so that all sub-slices will also be valid random
216 samples. This allows raffle winners (the sample) to be partitioned
217 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000218
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000219 Members of the population need not be hashable or unique. If the
220 population contains repeats, then each occurrence is a possible
221 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000222
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000223 To choose a sample in a range of integers, use xrange as an argument.
224 This is especially fast and space efficient for sampling from a
225 large population: sample(xrange(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000226
227 Optional arg random is a 0-argument function returning a random
228 float in [0.0, 1.0); by default, the standard random.random.
229 """
230
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000231 # Sampling without replacement entails tracking either potential
232 # selections (the pool) or previous selections.
233
234 # Pools are stored in lists which provide __getitem__ for selection
235 # and provide a way to remove selections. But each list.remove()
236 # rebuilds the entire list, so it is better to rearrange the list,
237 # placing non-selected elements at the head of the list. Tracking
238 # the selection pool is only space efficient with small populations.
239
240 # Previous selections are stored in dictionaries which provide
241 # __contains__ for detecting repeat selections. Discarding repeats
242 # is efficient unless most of the population has already been chosen.
Raymond Hettinger311f4192002-11-18 09:01:24 +0000243 # So, tracking selections is fast only with small sample sizes.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000244
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000245 n = len(population)
246 if not 0 <= k <= n:
247 raise ValueError, "sample larger than population"
248 if random is None:
249 random = self.random
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000250 result = [None] * k
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000251 if n < 6 * k: # if n len list takes less space than a k len dict
Raymond Hettinger311f4192002-11-18 09:01:24 +0000252 pool = list(population)
253 for i in xrange(k): # invariant: non-selected at [0,n-i)
254 j = int(random() * (n-i))
255 result[i] = pool[j]
256 pool[j] = pool[n-i-1]
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000257 else:
Raymond Hettinger311f4192002-11-18 09:01:24 +0000258 selected = {}
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000259 for i in xrange(k):
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000260 j = int(random() * n)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000261 while j in selected:
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000262 j = int(random() * n)
263 result[i] = selected[j] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000264 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000265
Tim Peterscd804102001-01-25 20:25:57 +0000266## -------------------- real-valued distributions -------------------
267
268## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000269
270 def uniform(self, a, b):
271 """Get a random number in the range [a, b)."""
272 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000273
Tim Peterscd804102001-01-25 20:25:57 +0000274## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000275
Tim Petersd7b5e882001-01-25 03:36:26 +0000276 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000277 """Normal distribution.
278
279 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000280
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000281 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000282 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000283
Tim Petersd7b5e882001-01-25 03:36:26 +0000284 # Uses Kinderman and Monahan method. Reference: Kinderman,
285 # A.J. and Monahan, J.F., "Computer generation of random
286 # variables using the ratio of uniform deviates", ACM Trans
287 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000288
Tim Petersd7b5e882001-01-25 03:36:26 +0000289 random = self.random
Raymond Hettinger311f4192002-11-18 09:01:24 +0000290 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000291 u1 = random()
292 u2 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000293 z = NV_MAGICCONST*(u1-0.5)/u2
294 zz = z*z/4.0
295 if zz <= -_log(u2):
296 break
297 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000298
Tim Peterscd804102001-01-25 20:25:57 +0000299## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000300
301 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000302 """Log normal distribution.
303
304 If you take the natural logarithm of this distribution, you'll get a
305 normal distribution with mean mu and standard deviation sigma.
306 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000307
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000308 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000309 return _exp(self.normalvariate(mu, sigma))
310
Tim Peterscd804102001-01-25 20:25:57 +0000311## -------------------- circular uniform --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000312
313 def cunifvariate(self, mean, arc):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000314 """Circular uniform distribution.
315
316 mean is the mean angle, and arc is the range of the distribution,
317 centered around the mean angle. Both values must be expressed in
318 radians. Returned values range between mean - arc/2 and
319 mean + arc/2 and are normalized to between 0 and pi.
320
321 Deprecated in version 2.3. Use:
322 (mean + arc * (Random.random() - 0.5)) % Math.pi
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000323
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000324 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000325 # mean: mean angle (in radians between 0 and pi)
326 # arc: range of distribution (in radians between 0 and pi)
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000327 import warnings
328 warnings.warn("The cunifvariate function is deprecated; Use (mean "
329 "+ arc * (Random.random() - 0.5)) % Math.pi instead",
330 DeprecationWarning)
Tim Petersd7b5e882001-01-25 03:36:26 +0000331
332 return (mean + arc * (self.random() - 0.5)) % _pi
333
Tim Peterscd804102001-01-25 20:25:57 +0000334## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000335
336 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000337 """Exponential distribution.
338
339 lambd is 1.0 divided by the desired mean. (The parameter would be
340 called "lambda", but that is a reserved word in Python.) Returned
341 values range from 0 to positive infinity.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000342
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000343 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000344 # lambd: rate lambd = 1/mean
345 # ('lambda' is a Python reserved word)
346
347 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000348 u = random()
349 while u <= 1e-7:
350 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000351 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000352
Tim Peterscd804102001-01-25 20:25:57 +0000353## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000354
Tim Petersd7b5e882001-01-25 03:36:26 +0000355 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000356 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000357
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000358 mu is the mean angle, expressed in radians between 0 and 2*pi, and
359 kappa is the concentration parameter, which must be greater than or
360 equal to zero. If kappa is equal to zero, this distribution reduces
361 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000362
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000363 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000364 # mu: mean angle (in radians between 0 and 2*pi)
365 # kappa: concentration parameter kappa (>= 0)
366 # if kappa = 0 generate uniform random angle
367
368 # Based upon an algorithm published in: Fisher, N.I.,
369 # "Statistical Analysis of Circular Data", Cambridge
370 # University Press, 1993.
371
372 # Thanks to Magnus Kessler for a correction to the
373 # implementation of step 4.
374
375 random = self.random
376 if kappa <= 1e-6:
377 return TWOPI * random()
378
379 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
380 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
381 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000382
Raymond Hettinger311f4192002-11-18 09:01:24 +0000383 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000384 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000385
386 z = _cos(_pi * u1)
387 f = (1.0 + r * z)/(r + z)
388 c = kappa * (r - f)
389
390 u2 = random()
391
392 if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
Tim Peters0c9886d2001-01-15 01:18:21 +0000393 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000394
395 u3 = random()
396 if u3 > 0.5:
397 theta = (mu % TWOPI) + _acos(f)
398 else:
399 theta = (mu % TWOPI) - _acos(f)
400
401 return theta
402
Tim Peterscd804102001-01-25 20:25:57 +0000403## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000404
405 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000406 """Gamma distribution. Not the gamma function!
407
408 Conditions on the parameters are alpha > 0 and beta > 0.
409
410 """
Tim Peters8ac14952002-05-23 15:15:30 +0000411
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000412 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000413
Guido van Rossum570764d2002-05-14 14:08:12 +0000414 # Warning: a few older sources define the gamma distribution in terms
415 # of alpha > -1.0
416 if alpha <= 0.0 or beta <= 0.0:
417 raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters8ac14952002-05-23 15:15:30 +0000418
Tim Petersd7b5e882001-01-25 03:36:26 +0000419 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000420 if alpha > 1.0:
421
422 # Uses R.C.H. Cheng, "The generation of Gamma
423 # variables with non-integral shape parameters",
424 # Applied Statistics, (1977), 26, No. 1, p71-74
425
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000426 ainv = _sqrt(2.0 * alpha - 1.0)
427 bbb = alpha - LOG4
428 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000429
Raymond Hettinger311f4192002-11-18 09:01:24 +0000430 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000431 u1 = random()
432 u2 = random()
433 v = _log(u1/(1.0-u1))/ainv
434 x = alpha*_exp(v)
435 z = u1*u1*u2
436 r = bbb+ccc*v-x
437 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000438 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000439
440 elif alpha == 1.0:
441 # expovariate(1)
442 u = random()
443 while u <= 1e-7:
444 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000445 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000446
447 else: # alpha is between 0 and 1 (exclusive)
448
449 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
450
Raymond Hettinger311f4192002-11-18 09:01:24 +0000451 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000452 u = random()
453 b = (_e + alpha)/_e
454 p = b*u
455 if p <= 1.0:
456 x = pow(p, 1.0/alpha)
457 else:
458 # p > 1
459 x = -_log((b-p)/alpha)
460 u1 = random()
461 if not (((p <= 1.0) and (u1 > _exp(-x))) or
462 ((p > 1) and (u1 > pow(x, alpha - 1.0)))):
463 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000464 return x * beta
465
466
467 def stdgamma(self, alpha, ainv, bbb, ccc):
468 # This method was (and shall remain) undocumented.
469 # This method is deprecated
470 # for the following reasons:
471 # 1. Returns same as .gammavariate(alpha, 1.0)
472 # 2. Requires caller to provide 3 extra arguments
473 # that are functions of alpha anyway
474 # 3. Can't be used for alpha < 0.5
475
476 # ainv = sqrt(2 * alpha - 1)
477 # bbb = alpha - log(4)
478 # ccc = alpha + ainv
479 import warnings
480 warnings.warn("The stdgamma function is deprecated; "
481 "use gammavariate() instead",
482 DeprecationWarning)
483 return self.gammavariate(alpha, 1.0)
484
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000485
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000486
Tim Peterscd804102001-01-25 20:25:57 +0000487## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000488
Tim Petersd7b5e882001-01-25 03:36:26 +0000489 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000490 """Gaussian distribution.
491
492 mu is the mean, and sigma is the standard deviation. This is
493 slightly faster than the normalvariate() function.
494
495 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000496
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000497 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000498
Tim Petersd7b5e882001-01-25 03:36:26 +0000499 # When x and y are two variables from [0, 1), uniformly
500 # distributed, then
501 #
502 # cos(2*pi*x)*sqrt(-2*log(1-y))
503 # sin(2*pi*x)*sqrt(-2*log(1-y))
504 #
505 # are two *independent* variables with normal distribution
506 # (mu = 0, sigma = 1).
507 # (Lambert Meertens)
508 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000509
Tim Petersd7b5e882001-01-25 03:36:26 +0000510 # Multithreading note: When two threads call this function
511 # simultaneously, it is possible that they will receive the
512 # same return value. The window is very small though. To
513 # avoid this, you have to use a lock around all calls. (I
514 # didn't want to slow this down in the serial case by using a
515 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000516
Tim Petersd7b5e882001-01-25 03:36:26 +0000517 random = self.random
518 z = self.gauss_next
519 self.gauss_next = None
520 if z is None:
521 x2pi = random() * TWOPI
522 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
523 z = _cos(x2pi) * g2rad
524 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000525
Tim Petersd7b5e882001-01-25 03:36:26 +0000526 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000527
Tim Peterscd804102001-01-25 20:25:57 +0000528## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000529## See
530## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
531## for Ivan Frohne's insightful analysis of why the original implementation:
532##
533## def betavariate(self, alpha, beta):
534## # Discrete Event Simulation in C, pp 87-88.
535##
536## y = self.expovariate(alpha)
537## z = self.expovariate(1.0/beta)
538## return z/(y+z)
539##
540## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000541
Tim Petersd7b5e882001-01-25 03:36:26 +0000542 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000543 """Beta distribution.
544
545 Conditions on the parameters are alpha > -1 and beta} > -1.
546 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000547
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000548 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000549
Tim Peters85e2e472001-01-26 06:49:56 +0000550 # This version due to Janne Sinkkonen, and matches all the std
551 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
552 y = self.gammavariate(alpha, 1.)
553 if y == 0:
554 return 0.0
555 else:
556 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000557
Tim Peterscd804102001-01-25 20:25:57 +0000558## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000559
Tim Petersd7b5e882001-01-25 03:36:26 +0000560 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000561 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000562 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000563
Tim Petersd7b5e882001-01-25 03:36:26 +0000564 u = self.random()
565 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000566
Tim Peterscd804102001-01-25 20:25:57 +0000567## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000568
Tim Petersd7b5e882001-01-25 03:36:26 +0000569 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000570 """Weibull distribution.
571
572 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000573
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000574 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000575 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000576
Tim Petersd7b5e882001-01-25 03:36:26 +0000577 u = self.random()
578 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000579
Raymond Hettinger40f62172002-12-29 23:03:38 +0000580## -------------------- Wichmann-Hill -------------------
581
582class WichmannHill(Random):
583
584 VERSION = 1 # used by getstate/setstate
585
586 def seed(self, a=None):
587 """Initialize internal state from hashable object.
588
589 None or no argument seeds from current time.
590
591 If a is not None or an int or long, hash(a) is used instead.
592
593 If a is an int or long, a is used directly. Distinct values between
594 0 and 27814431486575L inclusive are guaranteed to yield distinct
595 internal states (this guarantee is specific to the default
596 Wichmann-Hill generator).
597 """
598
599 if a is None:
600 # Initialize from current time
601 import time
602 a = long(time.time() * 256)
603
604 if not isinstance(a, (int, long)):
605 a = hash(a)
606
607 a, x = divmod(a, 30268)
608 a, y = divmod(a, 30306)
609 a, z = divmod(a, 30322)
610 self._seed = int(x)+1, int(y)+1, int(z)+1
611
612 self.gauss_next = None
613
614 def random(self):
615 """Get the next random number in the range [0.0, 1.0)."""
616
617 # Wichman-Hill random number generator.
618 #
619 # Wichmann, B. A. & Hill, I. D. (1982)
620 # Algorithm AS 183:
621 # An efficient and portable pseudo-random number generator
622 # Applied Statistics 31 (1982) 188-190
623 #
624 # see also:
625 # Correction to Algorithm AS 183
626 # Applied Statistics 33 (1984) 123
627 #
628 # McLeod, A. I. (1985)
629 # A remark on Algorithm AS 183
630 # Applied Statistics 34 (1985),198-200
631
632 # This part is thread-unsafe:
633 # BEGIN CRITICAL SECTION
634 x, y, z = self._seed
635 x = (171 * x) % 30269
636 y = (172 * y) % 30307
637 z = (170 * z) % 30323
638 self._seed = x, y, z
639 # END CRITICAL SECTION
640
641 # Note: on a platform using IEEE-754 double arithmetic, this can
642 # never return 0.0 (asserted by Tim; proof too long for a comment).
643 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
644
645 def getstate(self):
646 """Return internal state; can be passed to setstate() later."""
647 return self.VERSION, self._seed, self.gauss_next
648
649 def setstate(self, state):
650 """Restore internal state from object returned by getstate()."""
651 version = state[0]
652 if version == 1:
653 version, self._seed, self.gauss_next = state
654 else:
655 raise ValueError("state with version %s passed to "
656 "Random.setstate() of version %s" %
657 (version, self.VERSION))
658
659 def jumpahead(self, n):
660 """Act as if n calls to random() were made, but quickly.
661
662 n is an int, greater than or equal to 0.
663
664 Example use: If you have 2 threads and know that each will
665 consume no more than a million random numbers, create two Random
666 objects r1 and r2, then do
667 r2.setstate(r1.getstate())
668 r2.jumpahead(1000000)
669 Then r1 and r2 will use guaranteed-disjoint segments of the full
670 period.
671 """
672
673 if not n >= 0:
674 raise ValueError("n must be >= 0")
675 x, y, z = self._seed
676 x = int(x * pow(171, n, 30269)) % 30269
677 y = int(y * pow(172, n, 30307)) % 30307
678 z = int(z * pow(170, n, 30323)) % 30323
679 self._seed = x, y, z
680
681 def __whseed(self, x=0, y=0, z=0):
682 """Set the Wichmann-Hill seed from (x, y, z).
683
684 These must be integers in the range [0, 256).
685 """
686
687 if not type(x) == type(y) == type(z) == int:
688 raise TypeError('seeds must be integers')
689 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
690 raise ValueError('seeds must be in range(0, 256)')
691 if 0 == x == y == z:
692 # Initialize from current time
693 import time
694 t = long(time.time() * 256)
695 t = int((t&0xffffff) ^ (t>>24))
696 t, x = divmod(t, 256)
697 t, y = divmod(t, 256)
698 t, z = divmod(t, 256)
699 # Zero is a poor seed, so substitute 1
700 self._seed = (x or 1, y or 1, z or 1)
701
702 self.gauss_next = None
703
704 def whseed(self, a=None):
705 """Seed from hashable object's hash code.
706
707 None or no argument seeds from current time. It is not guaranteed
708 that objects with distinct hash codes lead to distinct internal
709 states.
710
711 This is obsolete, provided for compatibility with the seed routine
712 used prior to Python 2.1. Use the .seed() method instead.
713 """
714
715 if a is None:
716 self.__whseed()
717 return
718 a = hash(a)
719 a, x = divmod(a, 256)
720 a, y = divmod(a, 256)
721 a, z = divmod(a, 256)
722 x = (x + a) % 256 or 1
723 y = (y + a) % 256 or 1
724 z = (z + a) % 256 or 1
725 self.__whseed(x, y, z)
726
Tim Peterscd804102001-01-25 20:25:57 +0000727## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000728
Tim Petersd7b5e882001-01-25 03:36:26 +0000729def _test_generator(n, funccall):
Tim Peters0c9886d2001-01-15 01:18:21 +0000730 import time
731 print n, 'times', funccall
732 code = compile(funccall, funccall, 'eval')
733 sum = 0.0
734 sqsum = 0.0
735 smallest = 1e10
736 largest = -1e10
737 t0 = time.time()
738 for i in range(n):
739 x = eval(code)
740 sum = sum + x
741 sqsum = sqsum + x*x
742 smallest = min(x, smallest)
743 largest = max(x, largest)
744 t1 = time.time()
745 print round(t1-t0, 3), 'sec,',
746 avg = sum/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000747 stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters0c9886d2001-01-15 01:18:21 +0000748 print 'avg %g, stddev %g, min %g, max %g' % \
749 (avg, stddev, smallest, largest)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000750
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000751def _sample_generator(n, k):
752 # Return a fixed element from the sample. Validates random ordering.
753 return sample(xrange(n), k)[k//2]
754
755def _test(N=2000):
Tim Petersd7b5e882001-01-25 03:36:26 +0000756 _test_generator(N, 'random()')
757 _test_generator(N, 'normalvariate(0.0, 1.0)')
758 _test_generator(N, 'lognormvariate(0.0, 1.0)')
759 _test_generator(N, 'cunifvariate(0.0, 1.0)')
760 _test_generator(N, 'expovariate(1.0)')
761 _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000762 _test_generator(N, 'gammavariate(0.01, 1.0)')
763 _test_generator(N, 'gammavariate(0.1, 1.0)')
Tim Peters8ac14952002-05-23 15:15:30 +0000764 _test_generator(N, 'gammavariate(0.1, 2.0)')
Tim Petersd7b5e882001-01-25 03:36:26 +0000765 _test_generator(N, 'gammavariate(0.5, 1.0)')
766 _test_generator(N, 'gammavariate(0.9, 1.0)')
767 _test_generator(N, 'gammavariate(1.0, 1.0)')
768 _test_generator(N, 'gammavariate(2.0, 1.0)')
769 _test_generator(N, 'gammavariate(20.0, 1.0)')
770 _test_generator(N, 'gammavariate(200.0, 1.0)')
771 _test_generator(N, 'gauss(0.0, 1.0)')
772 _test_generator(N, 'betavariate(3.0, 3.0)')
773 _test_generator(N, 'paretovariate(1.0)')
774 _test_generator(N, 'weibullvariate(1.0, 1.0)')
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000775 _test_generator(N, '_sample_generator(50, 5)') # expected s.d.: 14.4
776 _test_generator(N, '_sample_generator(50, 45)') # expected s.d.: 14.4
Tim Peterscd804102001-01-25 20:25:57 +0000777
Tim Peters715c4c42001-01-26 22:56:56 +0000778# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000779# as module-level functions. The functions share state across all uses
780#(both in the user's code and in the Python libraries), but that's fine
781# for most programs and is easier for the casual user than making them
782# instantiate their own Random() instance.
783
Tim Petersd7b5e882001-01-25 03:36:26 +0000784_inst = Random()
785seed = _inst.seed
786random = _inst.random
787uniform = _inst.uniform
788randint = _inst.randint
789choice = _inst.choice
790randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000791sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000792shuffle = _inst.shuffle
793normalvariate = _inst.normalvariate
794lognormvariate = _inst.lognormvariate
795cunifvariate = _inst.cunifvariate
796expovariate = _inst.expovariate
797vonmisesvariate = _inst.vonmisesvariate
798gammavariate = _inst.gammavariate
799stdgamma = _inst.stdgamma
800gauss = _inst.gauss
801betavariate = _inst.betavariate
802paretovariate = _inst.paretovariate
803weibullvariate = _inst.weibullvariate
804getstate = _inst.getstate
805setstate = _inst.setstate
Tim Petersd52269b2001-01-25 06:23:18 +0000806jumpahead = _inst.jumpahead
Tim Petersd7b5e882001-01-25 03:36:26 +0000807
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000808if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000809 _test()