blob: e2c675c553643a603de9d0a9f3f05a06237e0415 [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
Guido van Rossumff03b1a1994-03-09 12:55:02 +000021
Guido van Rossume7b146f2000-02-04 15:28:42 +000022 distributions on the circle (angles 0 to 2pi)
23 ---------------------------------------------
24 circular uniform
25 von Mises
26
27Translated from anonymously contributed C/C++ source.
28
Tim Peterse360d952001-01-26 10:00:39 +000029Multi-threading note: the random number generator used here is not thread-
30safe; it is possible that two calls return the same random value. However,
31you can instantiate a different instance of Random() in each thread to get
32generators that don't share state, then use .setstate() and .jumpahead() to
33move the generators to disjoint segments of the full period. For example,
34
35def create_generators(num, delta, firstseed=None):
36 ""\"Return list of num distinct generators.
37 Each generator has its own unique segment of delta elements from
38 Random.random()'s full period.
39 Seed the first generator with optional arg firstseed (default is
40 None, to seed from current time).
41 ""\"
42
43 from random import Random
44 g = Random(firstseed)
45 result = [g]
46 for i in range(num - 1):
47 laststate = g.getstate()
48 g = Random()
49 g.setstate(laststate)
50 g.jumpahead(delta)
51 result.append(g)
52 return result
53
54gens = create_generators(10, 1000000)
55
56That creates 10 distinct generators, which can be passed out to 10 distinct
57threads. The generators don't share state so can be called safely in
58parallel. So long as no thread calls its g.random() more than a million
59times (the second argument to create_generators), the sequences seen by
60each thread will not overlap.
61
62The period of the underlying Wichmann-Hill generator is 6,953,607,871,644,
63and that limits how far this technique can be pushed.
64
65Just for fun, note that since we know the period, .jumpahead() can also be
66used to "move backward in time":
67
68>>> g = Random(42) # arbitrary
69>>> g.random()
Tim Peters0de88fc2001-02-01 04:59:18 +0000700.25420336316883324
Tim Peterse360d952001-01-26 10:00:39 +000071>>> g.jumpahead(6953607871644L - 1) # move *back* one
72>>> g.random()
Tim Peters0de88fc2001-02-01 04:59:18 +0000730.25420336316883324
Guido van Rossume7b146f2000-02-04 15:28:42 +000074"""
Tim Petersd7b5e882001-01-25 03:36:26 +000075# XXX The docstring sucks.
Guido van Rossumd03e1191998-05-29 17:51:31 +000076
Tim Petersd7b5e882001-01-25 03:36:26 +000077from math import log as _log, exp as _exp, pi as _pi, e as _e
78from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Tim Peters9146f272002-08-16 03:41:39 +000079from math import floor as _floor
Guido van Rossumff03b1a1994-03-09 12:55:02 +000080
Raymond Hettingerf24eb352002-11-12 17:41:57 +000081__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000082 "randrange","shuffle","normalvariate","lognormvariate",
83 "cunifvariate","expovariate","vonmisesvariate","gammavariate",
84 "stdgamma","gauss","betavariate","paretovariate","weibullvariate",
85 "getstate","setstate","jumpahead","whseed"]
Tim Peters0e6d2132001-02-15 23:56:39 +000086
Tim Petersdc47a892001-11-25 21:12:43 +000087def _verify(name, computed, expected):
Tim Peters0c9886d2001-01-15 01:18:21 +000088 if abs(computed - expected) > 1e-7:
Tim Petersd7b5e882001-01-25 03:36:26 +000089 raise ValueError(
90 "computed value for %s deviates too much "
91 "(computed %g, expected %g)" % (name, computed, expected))
92
93NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersdc47a892001-11-25 21:12:43 +000094_verify('NV_MAGICCONST', NV_MAGICCONST, 1.71552776992141)
Tim Petersd7b5e882001-01-25 03:36:26 +000095
96TWOPI = 2.0*_pi
Tim Petersdc47a892001-11-25 21:12:43 +000097_verify('TWOPI', TWOPI, 6.28318530718)
Tim Petersd7b5e882001-01-25 03:36:26 +000098
99LOG4 = _log(4.0)
Tim Petersdc47a892001-11-25 21:12:43 +0000100_verify('LOG4', LOG4, 1.38629436111989)
Tim Petersd7b5e882001-01-25 03:36:26 +0000101
102SG_MAGICCONST = 1.0 + _log(4.5)
Tim Petersdc47a892001-11-25 21:12:43 +0000103_verify('SG_MAGICCONST', SG_MAGICCONST, 2.50407739677627)
Tim Petersd7b5e882001-01-25 03:36:26 +0000104
105del _verify
106
107# Translated by Guido van Rossum from C source provided by
108# Adrian Baddeley.
109
110class Random:
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000111 """Random number generator base class used by bound module functions.
112
113 Used to instantiate instances of Random to get generators that don't
114 share state. Especially useful for multi-threaded programs, creating
115 a different instance of Random for each thread, and using the jumpahead()
116 method to ensure that the generated sequences seen by each thread don't
117 overlap.
118
119 Class Random can also be subclassed if you want to use a different basic
120 generator of your own devising: in that case, override the following
121 methods: random(), seed(), getstate(), setstate() and jumpahead().
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000122
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000123 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000124
125 VERSION = 1 # used by getstate/setstate
126
127 def __init__(self, x=None):
128 """Initialize an instance.
129
130 Optional argument x controls seeding, as for Random.seed().
131 """
132
133 self.seed(x)
Tim Petersd7b5e882001-01-25 03:36:26 +0000134
Tim Peterscd804102001-01-25 20:25:57 +0000135## -------------------- core generator -------------------
136
Tim Petersd7b5e882001-01-25 03:36:26 +0000137 # Specific to Wichmann-Hill generator. Subclasses wishing to use a
Tim Petersd52269b2001-01-25 06:23:18 +0000138 # different core generator should override the seed(), random(),
Tim Peterscd804102001-01-25 20:25:57 +0000139 # getstate(), setstate() and jumpahead() methods.
Tim Petersd7b5e882001-01-25 03:36:26 +0000140
Tim Peters0de88fc2001-02-01 04:59:18 +0000141 def seed(self, a=None):
142 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +0000143
Tim Peters0de88fc2001-02-01 04:59:18 +0000144 None or no argument seeds from current time.
145
Tim Petersbcd725f2001-02-01 10:06:53 +0000146 If a is not None or an int or long, hash(a) is used instead.
Tim Peters0de88fc2001-02-01 04:59:18 +0000147
148 If a is an int or long, a is used directly. Distinct values between
149 0 and 27814431486575L inclusive are guaranteed to yield distinct
150 internal states (this guarantee is specific to the default
151 Wichmann-Hill generator).
Tim Petersd7b5e882001-01-25 03:36:26 +0000152 """
153
Tim Peters0de88fc2001-02-01 04:59:18 +0000154 if a is None:
Tim Petersd7b5e882001-01-25 03:36:26 +0000155 # Initialize from current time
156 import time
Tim Peters0de88fc2001-02-01 04:59:18 +0000157 a = long(time.time() * 256)
158
159 if type(a) not in (type(3), type(3L)):
160 a = hash(a)
161
162 a, x = divmod(a, 30268)
163 a, y = divmod(a, 30306)
164 a, z = divmod(a, 30322)
165 self._seed = int(x)+1, int(y)+1, int(z)+1
Tim Petersd7b5e882001-01-25 03:36:26 +0000166
Tim Peters46c04e12002-05-05 20:40:00 +0000167 self.gauss_next = None
168
Tim Petersd7b5e882001-01-25 03:36:26 +0000169 def random(self):
170 """Get the next random number in the range [0.0, 1.0)."""
171
172 # Wichman-Hill random number generator.
173 #
174 # Wichmann, B. A. & Hill, I. D. (1982)
175 # Algorithm AS 183:
176 # An efficient and portable pseudo-random number generator
177 # Applied Statistics 31 (1982) 188-190
178 #
179 # see also:
180 # Correction to Algorithm AS 183
181 # Applied Statistics 33 (1984) 123
182 #
183 # McLeod, A. I. (1985)
184 # A remark on Algorithm AS 183
185 # Applied Statistics 34 (1985),198-200
186
187 # This part is thread-unsafe:
188 # BEGIN CRITICAL SECTION
189 x, y, z = self._seed
190 x = (171 * x) % 30269
191 y = (172 * y) % 30307
192 z = (170 * z) % 30323
193 self._seed = x, y, z
194 # END CRITICAL SECTION
195
196 # Note: on a platform using IEEE-754 double arithmetic, this can
197 # never return 0.0 (asserted by Tim; proof too long for a comment).
198 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
199
Tim Peterscd804102001-01-25 20:25:57 +0000200 def getstate(self):
201 """Return internal state; can be passed to setstate() later."""
202 return self.VERSION, self._seed, self.gauss_next
203
204 def setstate(self, state):
205 """Restore internal state from object returned by getstate()."""
206 version = state[0]
207 if version == 1:
208 version, self._seed, self.gauss_next = state
209 else:
210 raise ValueError("state with version %s passed to "
211 "Random.setstate() of version %s" %
212 (version, self.VERSION))
213
214 def jumpahead(self, n):
215 """Act as if n calls to random() were made, but quickly.
216
217 n is an int, greater than or equal to 0.
218
219 Example use: If you have 2 threads and know that each will
220 consume no more than a million random numbers, create two Random
221 objects r1 and r2, then do
222 r2.setstate(r1.getstate())
223 r2.jumpahead(1000000)
224 Then r1 and r2 will use guaranteed-disjoint segments of the full
225 period.
226 """
227
228 if not n >= 0:
229 raise ValueError("n must be >= 0")
230 x, y, z = self._seed
231 x = int(x * pow(171, n, 30269)) % 30269
232 y = int(y * pow(172, n, 30307)) % 30307
233 z = int(z * pow(170, n, 30323)) % 30323
234 self._seed = x, y, z
235
Tim Peters0de88fc2001-02-01 04:59:18 +0000236 def __whseed(self, x=0, y=0, z=0):
237 """Set the Wichmann-Hill seed from (x, y, z).
238
239 These must be integers in the range [0, 256).
240 """
241
242 if not type(x) == type(y) == type(z) == type(0):
243 raise TypeError('seeds must be integers')
244 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
245 raise ValueError('seeds must be in range(0, 256)')
246 if 0 == x == y == z:
247 # Initialize from current time
248 import time
249 t = long(time.time() * 256)
250 t = int((t&0xffffff) ^ (t>>24))
251 t, x = divmod(t, 256)
252 t, y = divmod(t, 256)
253 t, z = divmod(t, 256)
254 # Zero is a poor seed, so substitute 1
255 self._seed = (x or 1, y or 1, z or 1)
256
Tim Peters46c04e12002-05-05 20:40:00 +0000257 self.gauss_next = None
258
Tim Peters0de88fc2001-02-01 04:59:18 +0000259 def whseed(self, a=None):
260 """Seed from hashable object's hash code.
261
262 None or no argument seeds from current time. It is not guaranteed
263 that objects with distinct hash codes lead to distinct internal
264 states.
265
266 This is obsolete, provided for compatibility with the seed routine
267 used prior to Python 2.1. Use the .seed() method instead.
268 """
269
270 if a is None:
271 self.__whseed()
272 return
273 a = hash(a)
274 a, x = divmod(a, 256)
275 a, y = divmod(a, 256)
276 a, z = divmod(a, 256)
277 x = (x + a) % 256 or 1
278 y = (y + a) % 256 or 1
279 z = (z + a) % 256 or 1
280 self.__whseed(x, y, z)
281
Tim Peterscd804102001-01-25 20:25:57 +0000282## ---- Methods below this point do not need to be overridden when
283## ---- subclassing for the purpose of using a different core generator.
284
285## -------------------- pickle support -------------------
286
287 def __getstate__(self): # for pickle
288 return self.getstate()
289
290 def __setstate__(self, state): # for pickle
291 self.setstate(state)
292
293## -------------------- integer methods -------------------
294
Tim Petersd7b5e882001-01-25 03:36:26 +0000295 def randrange(self, start, stop=None, step=1, int=int, default=None):
296 """Choose a random item from range(start, stop[, step]).
297
298 This fixes the problem with randint() which includes the
299 endpoint; in Python this is usually not what you want.
300 Do not supply the 'int' and 'default' arguments.
301 """
302
303 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000304 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000305 istart = int(start)
306 if istart != start:
307 raise ValueError, "non-integer arg 1 for randrange()"
308 if stop is default:
309 if istart > 0:
310 return int(self.random() * istart)
311 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000312
313 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000314 istop = int(stop)
315 if istop != stop:
316 raise ValueError, "non-integer stop for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000317 if step == 1 and istart < istop:
318 try:
319 return istart + int(self.random()*(istop - istart))
320 except OverflowError:
321 # This can happen if istop-istart > sys.maxint + 1, and
322 # multiplying by random() doesn't reduce it to something
323 # <= sys.maxint. We know that the overall result fits
324 # in an int, and can still do it correctly via math.floor().
325 # But that adds another function call, so for speed we
326 # avoided that whenever possible.
327 return int(istart + _floor(self.random()*(istop - istart)))
Tim Petersd7b5e882001-01-25 03:36:26 +0000328 if step == 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000329 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000330
331 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000332 istep = int(step)
333 if istep != step:
334 raise ValueError, "non-integer step for randrange()"
335 if istep > 0:
336 n = (istop - istart + istep - 1) / istep
337 elif istep < 0:
338 n = (istop - istart + istep + 1) / istep
339 else:
340 raise ValueError, "zero step for randrange()"
341
342 if n <= 0:
343 raise ValueError, "empty range for randrange()"
344 return istart + istep*int(self.random() * n)
345
346 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000347 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000348 """
349
350 return self.randrange(a, b+1)
351
Tim Peterscd804102001-01-25 20:25:57 +0000352## -------------------- sequence methods -------------------
353
Tim Petersd7b5e882001-01-25 03:36:26 +0000354 def choice(self, seq):
355 """Choose a random element from a non-empty sequence."""
356 return seq[int(self.random() * len(seq))]
357
358 def shuffle(self, x, random=None, int=int):
359 """x, random=random.random -> shuffle list x in place; return None.
360
361 Optional arg random is a 0-argument function returning a random
362 float in [0.0, 1.0); by default, the standard random.random.
363
364 Note that for even rather small len(x), the total number of
365 permutations of x is larger than the period of most random number
366 generators; this implies that "most" permutations of a long
367 sequence can never be generated.
368 """
369
370 if random is None:
371 random = self.random
372 for i in xrange(len(x)-1, 0, -1):
Tim Peterscd804102001-01-25 20:25:57 +0000373 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000374 j = int(random() * (i+1))
375 x[i], x[j] = x[j], x[i]
376
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000377 def sample(self, population, k, random=None, int=int):
378 """Chooses k unique random elements from a population sequence.
379
380 Returns a new list containing elements from the population. The
381 list itself is in random order so that all sub-slices are also
382 random samples. The original sequence is left undisturbed.
383
384 If the population has repeated elements, then each occurence is
385 a possible selection in the sample.
386
387 If indices are needed for a large population, use xrange as an
388 argument: sample(xrange(10000000), 60)
389
390 Optional arg random is a 0-argument function returning a random
391 float in [0.0, 1.0); by default, the standard random.random.
392 """
393
394 n = len(population)
395 if not 0 <= k <= n:
396 raise ValueError, "sample larger than population"
397 if random is None:
398 random = self.random
399 if n < 6 * k: # if n len list takes less space than a k len dict
400 pool = list(population)
401 for i in xrange(n-1, n-k-1, -1):
402 j = int(random() * (i+1))
403 pool[i], pool[j] = pool[j], pool[i]
404 return pool[-k:]
405 inorder = [None] * k
406 selections = {}
407 for i in xrange(k):
408 j = int(random() * n)
409 while j in selections:
410 j = int(random() * n)
411 selections[j] = inorder[i] = population[j]
412 return inorder # return selections in the order they were picked
413
Tim Peterscd804102001-01-25 20:25:57 +0000414## -------------------- real-valued distributions -------------------
415
416## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000417
418 def uniform(self, a, b):
419 """Get a random number in the range [a, b)."""
420 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000421
Tim Peterscd804102001-01-25 20:25:57 +0000422## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000423
Tim Petersd7b5e882001-01-25 03:36:26 +0000424 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000425 """Normal distribution.
426
427 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000428
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000429 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000430 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000431
Tim Petersd7b5e882001-01-25 03:36:26 +0000432 # Uses Kinderman and Monahan method. Reference: Kinderman,
433 # A.J. and Monahan, J.F., "Computer generation of random
434 # variables using the ratio of uniform deviates", ACM Trans
435 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000436
Tim Petersd7b5e882001-01-25 03:36:26 +0000437 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000438 while 1:
439 u1 = random()
440 u2 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000441 z = NV_MAGICCONST*(u1-0.5)/u2
442 zz = z*z/4.0
443 if zz <= -_log(u2):
444 break
445 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000446
Tim Peterscd804102001-01-25 20:25:57 +0000447## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000448
449 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000450 """Log normal distribution.
451
452 If you take the natural logarithm of this distribution, you'll get a
453 normal distribution with mean mu and standard deviation sigma.
454 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000455
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000456 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000457 return _exp(self.normalvariate(mu, sigma))
458
Tim Peterscd804102001-01-25 20:25:57 +0000459## -------------------- circular uniform --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000460
461 def cunifvariate(self, mean, arc):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000462 """Circular uniform distribution.
463
464 mean is the mean angle, and arc is the range of the distribution,
465 centered around the mean angle. Both values must be expressed in
466 radians. Returned values range between mean - arc/2 and
467 mean + arc/2 and are normalized to between 0 and pi.
468
469 Deprecated in version 2.3. Use:
470 (mean + arc * (Random.random() - 0.5)) % Math.pi
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000471
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000472 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000473 # mean: mean angle (in radians between 0 and pi)
474 # arc: range of distribution (in radians between 0 and pi)
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000475 import warnings
476 warnings.warn("The cunifvariate function is deprecated; Use (mean "
477 "+ arc * (Random.random() - 0.5)) % Math.pi instead",
478 DeprecationWarning)
Tim Petersd7b5e882001-01-25 03:36:26 +0000479
480 return (mean + arc * (self.random() - 0.5)) % _pi
481
Tim Peterscd804102001-01-25 20:25:57 +0000482## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000483
484 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000485 """Exponential distribution.
486
487 lambd is 1.0 divided by the desired mean. (The parameter would be
488 called "lambda", but that is a reserved word in Python.) Returned
489 values range from 0 to positive infinity.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000490
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000491 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000492 # lambd: rate lambd = 1/mean
493 # ('lambda' is a Python reserved word)
494
495 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000496 u = random()
497 while u <= 1e-7:
498 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000499 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000500
Tim Peterscd804102001-01-25 20:25:57 +0000501## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000502
Tim Petersd7b5e882001-01-25 03:36:26 +0000503 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000504 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000505
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000506 mu is the mean angle, expressed in radians between 0 and 2*pi, and
507 kappa is the concentration parameter, which must be greater than or
508 equal to zero. If kappa is equal to zero, this distribution reduces
509 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000510
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000511 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000512 # mu: mean angle (in radians between 0 and 2*pi)
513 # kappa: concentration parameter kappa (>= 0)
514 # if kappa = 0 generate uniform random angle
515
516 # Based upon an algorithm published in: Fisher, N.I.,
517 # "Statistical Analysis of Circular Data", Cambridge
518 # University Press, 1993.
519
520 # Thanks to Magnus Kessler for a correction to the
521 # implementation of step 4.
522
523 random = self.random
524 if kappa <= 1e-6:
525 return TWOPI * random()
526
527 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
528 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
529 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000530
Tim Peters0c9886d2001-01-15 01:18:21 +0000531 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000532 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000533
534 z = _cos(_pi * u1)
535 f = (1.0 + r * z)/(r + z)
536 c = kappa * (r - f)
537
538 u2 = random()
539
540 if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
Tim Peters0c9886d2001-01-15 01:18:21 +0000541 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000542
543 u3 = random()
544 if u3 > 0.5:
545 theta = (mu % TWOPI) + _acos(f)
546 else:
547 theta = (mu % TWOPI) - _acos(f)
548
549 return theta
550
Tim Peterscd804102001-01-25 20:25:57 +0000551## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000552
553 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000554 """Gamma distribution. Not the gamma function!
555
556 Conditions on the parameters are alpha > 0 and beta > 0.
557
558 """
Tim Peters8ac14952002-05-23 15:15:30 +0000559
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000560 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000561
Guido van Rossum570764d2002-05-14 14:08:12 +0000562 # Warning: a few older sources define the gamma distribution in terms
563 # of alpha > -1.0
564 if alpha <= 0.0 or beta <= 0.0:
565 raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters8ac14952002-05-23 15:15:30 +0000566
Tim Petersd7b5e882001-01-25 03:36:26 +0000567 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000568 if alpha > 1.0:
569
570 # Uses R.C.H. Cheng, "The generation of Gamma
571 # variables with non-integral shape parameters",
572 # Applied Statistics, (1977), 26, No. 1, p71-74
573
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000574 ainv = _sqrt(2.0 * alpha - 1.0)
575 bbb = alpha - LOG4
576 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000577
Tim Petersd7b5e882001-01-25 03:36:26 +0000578 while 1:
579 u1 = random()
580 u2 = random()
581 v = _log(u1/(1.0-u1))/ainv
582 x = alpha*_exp(v)
583 z = u1*u1*u2
584 r = bbb+ccc*v-x
585 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000586 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000587
588 elif alpha == 1.0:
589 # expovariate(1)
590 u = random()
591 while u <= 1e-7:
592 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000593 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000594
595 else: # alpha is between 0 and 1 (exclusive)
596
597 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
598
599 while 1:
600 u = random()
601 b = (_e + alpha)/_e
602 p = b*u
603 if p <= 1.0:
604 x = pow(p, 1.0/alpha)
605 else:
606 # p > 1
607 x = -_log((b-p)/alpha)
608 u1 = random()
609 if not (((p <= 1.0) and (u1 > _exp(-x))) or
610 ((p > 1) and (u1 > pow(x, alpha - 1.0)))):
611 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000612 return x * beta
613
614
615 def stdgamma(self, alpha, ainv, bbb, ccc):
616 # This method was (and shall remain) undocumented.
617 # This method is deprecated
618 # for the following reasons:
619 # 1. Returns same as .gammavariate(alpha, 1.0)
620 # 2. Requires caller to provide 3 extra arguments
621 # that are functions of alpha anyway
622 # 3. Can't be used for alpha < 0.5
623
624 # ainv = sqrt(2 * alpha - 1)
625 # bbb = alpha - log(4)
626 # ccc = alpha + ainv
627 import warnings
628 warnings.warn("The stdgamma function is deprecated; "
629 "use gammavariate() instead",
630 DeprecationWarning)
631 return self.gammavariate(alpha, 1.0)
632
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000633
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000634
Tim Peterscd804102001-01-25 20:25:57 +0000635## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000636
Tim Petersd7b5e882001-01-25 03:36:26 +0000637 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000638 """Gaussian distribution.
639
640 mu is the mean, and sigma is the standard deviation. This is
641 slightly faster than the normalvariate() function.
642
643 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000644
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000645 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000646
Tim Petersd7b5e882001-01-25 03:36:26 +0000647 # When x and y are two variables from [0, 1), uniformly
648 # distributed, then
649 #
650 # cos(2*pi*x)*sqrt(-2*log(1-y))
651 # sin(2*pi*x)*sqrt(-2*log(1-y))
652 #
653 # are two *independent* variables with normal distribution
654 # (mu = 0, sigma = 1).
655 # (Lambert Meertens)
656 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000657
Tim Petersd7b5e882001-01-25 03:36:26 +0000658 # Multithreading note: When two threads call this function
659 # simultaneously, it is possible that they will receive the
660 # same return value. The window is very small though. To
661 # avoid this, you have to use a lock around all calls. (I
662 # didn't want to slow this down in the serial case by using a
663 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000664
Tim Petersd7b5e882001-01-25 03:36:26 +0000665 random = self.random
666 z = self.gauss_next
667 self.gauss_next = None
668 if z is None:
669 x2pi = random() * TWOPI
670 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
671 z = _cos(x2pi) * g2rad
672 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000673
Tim Petersd7b5e882001-01-25 03:36:26 +0000674 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000675
Tim Peterscd804102001-01-25 20:25:57 +0000676## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000677## See
678## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
679## for Ivan Frohne's insightful analysis of why the original implementation:
680##
681## def betavariate(self, alpha, beta):
682## # Discrete Event Simulation in C, pp 87-88.
683##
684## y = self.expovariate(alpha)
685## z = self.expovariate(1.0/beta)
686## return z/(y+z)
687##
688## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000689
Tim Petersd7b5e882001-01-25 03:36:26 +0000690 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000691 """Beta distribution.
692
693 Conditions on the parameters are alpha > -1 and beta} > -1.
694 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000695
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000696 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000697
Tim Peters85e2e472001-01-26 06:49:56 +0000698 # This version due to Janne Sinkkonen, and matches all the std
699 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
700 y = self.gammavariate(alpha, 1.)
701 if y == 0:
702 return 0.0
703 else:
704 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000705
Tim Peterscd804102001-01-25 20:25:57 +0000706## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000707
Tim Petersd7b5e882001-01-25 03:36:26 +0000708 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000709 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000710 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000711
Tim Petersd7b5e882001-01-25 03:36:26 +0000712 u = self.random()
713 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000714
Tim Peterscd804102001-01-25 20:25:57 +0000715## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000716
Tim Petersd7b5e882001-01-25 03:36:26 +0000717 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000718 """Weibull distribution.
719
720 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000721
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000722 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000723 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000724
Tim Petersd7b5e882001-01-25 03:36:26 +0000725 u = self.random()
726 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000727
Tim Peterscd804102001-01-25 20:25:57 +0000728## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000729
Tim Petersd7b5e882001-01-25 03:36:26 +0000730def _test_generator(n, funccall):
Tim Peters0c9886d2001-01-15 01:18:21 +0000731 import time
732 print n, 'times', funccall
733 code = compile(funccall, funccall, 'eval')
734 sum = 0.0
735 sqsum = 0.0
736 smallest = 1e10
737 largest = -1e10
738 t0 = time.time()
739 for i in range(n):
740 x = eval(code)
741 sum = sum + x
742 sqsum = sqsum + x*x
743 smallest = min(x, smallest)
744 largest = max(x, largest)
745 t1 = time.time()
746 print round(t1-t0, 3), 'sec,',
747 avg = sum/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000748 stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters0c9886d2001-01-15 01:18:21 +0000749 print 'avg %g, stddev %g, min %g, max %g' % \
750 (avg, stddev, smallest, largest)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000751
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000752def _test_sample(n):
753 # For the entire allowable range of 0 <= k <= n, validate that
754 # the sample is of the correct length and contains only unique items
755 population = xrange(n)
756 for k in xrange(n+1):
757 s = sample(population, k)
758 assert len(dict([(elem,True) for elem in s])) == len(s) == k
759
760def _sample_generator(n, k):
761 # Return a fixed element from the sample. Validates random ordering.
762 return sample(xrange(n), k)[k//2]
763
764def _test(N=2000):
Tim Petersd7b5e882001-01-25 03:36:26 +0000765 print 'TWOPI =', TWOPI
766 print 'LOG4 =', LOG4
767 print 'NV_MAGICCONST =', NV_MAGICCONST
768 print 'SG_MAGICCONST =', SG_MAGICCONST
769 _test_generator(N, 'random()')
770 _test_generator(N, 'normalvariate(0.0, 1.0)')
771 _test_generator(N, 'lognormvariate(0.0, 1.0)')
772 _test_generator(N, 'cunifvariate(0.0, 1.0)')
773 _test_generator(N, 'expovariate(1.0)')
774 _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000775 _test_generator(N, 'gammavariate(0.01, 1.0)')
776 _test_generator(N, 'gammavariate(0.1, 1.0)')
Tim Peters8ac14952002-05-23 15:15:30 +0000777 _test_generator(N, 'gammavariate(0.1, 2.0)')
Tim Petersd7b5e882001-01-25 03:36:26 +0000778 _test_generator(N, 'gammavariate(0.5, 1.0)')
779 _test_generator(N, 'gammavariate(0.9, 1.0)')
780 _test_generator(N, 'gammavariate(1.0, 1.0)')
781 _test_generator(N, 'gammavariate(2.0, 1.0)')
782 _test_generator(N, 'gammavariate(20.0, 1.0)')
783 _test_generator(N, 'gammavariate(200.0, 1.0)')
784 _test_generator(N, 'gauss(0.0, 1.0)')
785 _test_generator(N, 'betavariate(3.0, 3.0)')
786 _test_generator(N, 'paretovariate(1.0)')
787 _test_generator(N, 'weibullvariate(1.0, 1.0)')
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000788 _test_generator(N, '_sample_generator(50, 5)') # expected s.d.: 14.4
789 _test_generator(N, '_sample_generator(50, 45)') # expected s.d.: 14.4
790 _test_sample(1000)
Tim Petersd7b5e882001-01-25 03:36:26 +0000791
Tim Peterscd804102001-01-25 20:25:57 +0000792 # Test jumpahead.
793 s = getstate()
794 jumpahead(N)
795 r1 = random()
796 # now do it the slow way
797 setstate(s)
798 for i in range(N):
799 random()
800 r2 = random()
801 if r1 != r2:
802 raise ValueError("jumpahead test failed " + `(N, r1, r2)`)
803
Tim Peters715c4c42001-01-26 22:56:56 +0000804# Create one instance, seeded from current time, and export its methods
805# as module-level functions. The functions are not threadsafe, and state
806# is shared across all uses (both in the user's code and in the Python
807# libraries), but that's fine for most programs and is easier for the
808# casual user than making them instantiate their own Random() instance.
Tim Petersd7b5e882001-01-25 03:36:26 +0000809_inst = Random()
810seed = _inst.seed
811random = _inst.random
812uniform = _inst.uniform
813randint = _inst.randint
814choice = _inst.choice
815randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000816sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000817shuffle = _inst.shuffle
818normalvariate = _inst.normalvariate
819lognormvariate = _inst.lognormvariate
820cunifvariate = _inst.cunifvariate
821expovariate = _inst.expovariate
822vonmisesvariate = _inst.vonmisesvariate
823gammavariate = _inst.gammavariate
824stdgamma = _inst.stdgamma
825gauss = _inst.gauss
826betavariate = _inst.betavariate
827paretovariate = _inst.paretovariate
828weibullvariate = _inst.weibullvariate
829getstate = _inst.getstate
830setstate = _inst.setstate
Tim Petersd52269b2001-01-25 06:23:18 +0000831jumpahead = _inst.jumpahead
Tim Peters0de88fc2001-02-01 04:59:18 +0000832whseed = _inst.whseed
Tim Petersd7b5e882001-01-25 03:36:26 +0000833
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000834if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000835 _test()