|
1 """Random variable generators. |
|
2 |
|
3 integers |
|
4 -------- |
|
5 uniform within range |
|
6 |
|
7 sequences |
|
8 --------- |
|
9 pick random element |
|
10 pick random sample |
|
11 generate random permutation |
|
12 |
|
13 distributions on the real line: |
|
14 ------------------------------ |
|
15 uniform |
|
16 normal (Gaussian) |
|
17 lognormal |
|
18 negative exponential |
|
19 gamma |
|
20 beta |
|
21 pareto |
|
22 Weibull |
|
23 |
|
24 distributions on the circle (angles 0 to 2pi) |
|
25 --------------------------------------------- |
|
26 circular uniform |
|
27 von Mises |
|
28 |
|
29 General notes on the underlying Mersenne Twister core generator: |
|
30 |
|
31 * 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 semantics of |
|
34 jumpahead(n) are weakened to simply jump to another distant state and rely |
|
35 on the large period to avoid overlapping sequences. |
|
36 * The random() method is implemented in C, executes in a single Python step, |
|
37 and is, therefore, threadsafe. |
|
38 |
|
39 """ |
|
40 |
|
41 from warnings import warn as _warn |
|
42 from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType |
|
43 from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil |
|
44 from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin |
|
45 from os import urandom as _urandom |
|
46 from binascii import hexlify as _hexlify |
|
47 |
|
48 __all__ = ["Random","seed","random","uniform","randint","choice","sample", |
|
49 "randrange","shuffle","normalvariate","lognormvariate", |
|
50 "expovariate","vonmisesvariate","gammavariate", |
|
51 "gauss","betavariate","paretovariate","weibullvariate", |
|
52 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits", |
|
53 "SystemRandom"] |
|
54 |
|
55 NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) |
|
56 TWOPI = 2.0*_pi |
|
57 LOG4 = _log(4.0) |
|
58 SG_MAGICCONST = 1.0 + _log(4.5) |
|
59 BPF = 53 # Number of bits in a float |
|
60 RECIP_BPF = 2**-BPF |
|
61 |
|
62 |
|
63 # Translated by Guido van Rossum from C source provided by |
|
64 # Adrian Baddeley. Adapted by Raymond Hettinger for use with |
|
65 # the Mersenne Twister and os.urandom() core generators. |
|
66 |
|
67 import _random |
|
68 |
|
69 class Random(_random.Random): |
|
70 """Random number generator base class used by bound module functions. |
|
71 |
|
72 Used to instantiate instances of Random to get generators that don't |
|
73 share state. Especially useful for multi-threaded programs, creating |
|
74 a different instance of Random for each thread, and using the jumpahead() |
|
75 method to ensure that the generated sequences seen by each thread don't |
|
76 overlap. |
|
77 |
|
78 Class Random can also be subclassed if you want to use a different basic |
|
79 generator of your own devising: in that case, override the following |
|
80 methods: random(), seed(), getstate(), setstate() and jumpahead(). |
|
81 Optionally, implement a getrandombits() method so that randrange() |
|
82 can cover arbitrarily large ranges. |
|
83 |
|
84 """ |
|
85 |
|
86 VERSION = 2 # used by getstate/setstate |
|
87 |
|
88 def __init__(self, x=None): |
|
89 """Initialize an instance. |
|
90 |
|
91 Optional argument x controls seeding, as for Random.seed(). |
|
92 """ |
|
93 |
|
94 self.seed(x) |
|
95 self.gauss_next = None |
|
96 |
|
97 def seed(self, a=None): |
|
98 """Initialize internal state from hashable object. |
|
99 |
|
100 None or no argument seeds from current time or from an operating |
|
101 system specific randomness source if available. |
|
102 |
|
103 If a is not None or an int or long, hash(a) is used instead. |
|
104 """ |
|
105 |
|
106 if a is None: |
|
107 try: |
|
108 a = long(_hexlify(_urandom(16)), 16) |
|
109 except NotImplementedError: |
|
110 import time |
|
111 a = long(time.time() * 256) # use fractional seconds |
|
112 |
|
113 super(Random, self).seed(a) |
|
114 self.gauss_next = None |
|
115 |
|
116 def getstate(self): |
|
117 """Return internal state; can be passed to setstate() later.""" |
|
118 return self.VERSION, super(Random, self).getstate(), self.gauss_next |
|
119 |
|
120 def setstate(self, state): |
|
121 """Restore internal state from object returned by getstate().""" |
|
122 version = state[0] |
|
123 if version == 2: |
|
124 version, internalstate, self.gauss_next = state |
|
125 super(Random, self).setstate(internalstate) |
|
126 else: |
|
127 raise ValueError("state with version %s passed to " |
|
128 "Random.setstate() of version %s" % |
|
129 (version, self.VERSION)) |
|
130 |
|
131 ## ---- Methods below this point do not need to be overridden when |
|
132 ## ---- subclassing for the purpose of using a different core generator. |
|
133 |
|
134 ## -------------------- pickle support ------------------- |
|
135 |
|
136 def __getstate__(self): # for pickle |
|
137 return self.getstate() |
|
138 |
|
139 def __setstate__(self, state): # for pickle |
|
140 self.setstate(state) |
|
141 |
|
142 def __reduce__(self): |
|
143 return self.__class__, (), self.getstate() |
|
144 |
|
145 ## -------------------- integer methods ------------------- |
|
146 |
|
147 def randrange(self, start, stop=None, step=1, int=int, default=None, |
|
148 maxwidth=1L<<BPF): |
|
149 """Choose a random item from range(start, stop[, step]). |
|
150 |
|
151 This fixes the problem with randint() which includes the |
|
152 endpoint; in Python this is usually not what you want. |
|
153 Do not supply the 'int', 'default', and 'maxwidth' arguments. |
|
154 """ |
|
155 |
|
156 # This code is a bit messy to make it fast for the |
|
157 # common case while still doing adequate error checking. |
|
158 istart = int(start) |
|
159 if istart != start: |
|
160 raise ValueError, "non-integer arg 1 for randrange()" |
|
161 if stop is default: |
|
162 if istart > 0: |
|
163 if istart >= maxwidth: |
|
164 return self._randbelow(istart) |
|
165 return int(self.random() * istart) |
|
166 raise ValueError, "empty range for randrange()" |
|
167 |
|
168 # stop argument supplied. |
|
169 istop = int(stop) |
|
170 if istop != stop: |
|
171 raise ValueError, "non-integer stop for randrange()" |
|
172 width = istop - istart |
|
173 if step == 1 and width > 0: |
|
174 # Note that |
|
175 # int(istart + self.random()*width) |
|
176 # instead would be incorrect. For example, consider istart |
|
177 # = -2 and istop = 0. Then the guts would be in |
|
178 # -2.0 to 0.0 exclusive on both ends (ignoring that random() |
|
179 # might return 0.0), and because int() truncates toward 0, the |
|
180 # final result would be -1 or 0 (instead of -2 or -1). |
|
181 # istart + int(self.random()*width) |
|
182 # would also be incorrect, for a subtler reason: the RHS |
|
183 # can return a long, and then randrange() would also return |
|
184 # a long, but we're supposed to return an int (for backward |
|
185 # compatibility). |
|
186 |
|
187 if width >= maxwidth: |
|
188 return int(istart + self._randbelow(width)) |
|
189 return int(istart + int(self.random()*width)) |
|
190 if step == 1: |
|
191 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width) |
|
192 |
|
193 # Non-unit step argument supplied. |
|
194 istep = int(step) |
|
195 if istep != step: |
|
196 raise ValueError, "non-integer step for randrange()" |
|
197 if istep > 0: |
|
198 n = (width + istep - 1) // istep |
|
199 elif istep < 0: |
|
200 n = (width + istep + 1) // istep |
|
201 else: |
|
202 raise ValueError, "zero step for randrange()" |
|
203 |
|
204 if n <= 0: |
|
205 raise ValueError, "empty range for randrange()" |
|
206 |
|
207 if n >= maxwidth: |
|
208 return istart + istep*self._randbelow(n) |
|
209 return istart + istep*int(self.random() * n) |
|
210 |
|
211 def randint(self, a, b): |
|
212 """Return random integer in range [a, b], including both end points. |
|
213 """ |
|
214 |
|
215 return self.randrange(a, b+1) |
|
216 |
|
217 def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF, |
|
218 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType): |
|
219 """Return a random int in the range [0,n) |
|
220 |
|
221 Handles the case where n has more bits than returned |
|
222 by a single call to the underlying generator. |
|
223 """ |
|
224 |
|
225 try: |
|
226 getrandbits = self.getrandbits |
|
227 except AttributeError: |
|
228 pass |
|
229 else: |
|
230 # Only call self.getrandbits if the original random() builtin method |
|
231 # has not been overridden or if a new getrandbits() was supplied. |
|
232 # This assures that the two methods correspond. |
|
233 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method: |
|
234 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2) |
|
235 r = getrandbits(k) |
|
236 while r >= n: |
|
237 r = getrandbits(k) |
|
238 return r |
|
239 if n >= _maxwidth: |
|
240 _warn("Underlying random() generator does not supply \n" |
|
241 "enough bits to choose from a population range this large") |
|
242 return int(self.random() * n) |
|
243 |
|
244 ## -------------------- sequence methods ------------------- |
|
245 |
|
246 def choice(self, seq): |
|
247 """Choose a random element from a non-empty sequence.""" |
|
248 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty |
|
249 |
|
250 def shuffle(self, x, random=None, int=int): |
|
251 """x, random=random.random -> shuffle list x in place; return None. |
|
252 |
|
253 Optional arg random is a 0-argument function returning a random |
|
254 float in [0.0, 1.0); by default, the standard random.random. |
|
255 """ |
|
256 |
|
257 if random is None: |
|
258 random = self.random |
|
259 for i in reversed(xrange(1, len(x))): |
|
260 # pick an element in x[:i+1] with which to exchange x[i] |
|
261 j = int(random() * (i+1)) |
|
262 x[i], x[j] = x[j], x[i] |
|
263 |
|
264 def sample(self, population, k): |
|
265 """Chooses k unique random elements from a population sequence. |
|
266 |
|
267 Returns a new list containing elements from the population while |
|
268 leaving the original population unchanged. The resulting list is |
|
269 in selection order so that all sub-slices will also be valid random |
|
270 samples. This allows raffle winners (the sample) to be partitioned |
|
271 into grand prize and second place winners (the subslices). |
|
272 |
|
273 Members of the population need not be hashable or unique. If the |
|
274 population contains repeats, then each occurrence is a possible |
|
275 selection in the sample. |
|
276 |
|
277 To choose a sample in a range of integers, use xrange as an argument. |
|
278 This is especially fast and space efficient for sampling from a |
|
279 large population: sample(xrange(10000000), 60) |
|
280 """ |
|
281 |
|
282 # XXX Although the documentation says `population` is "a sequence", |
|
283 # XXX attempts are made to cater to any iterable with a __len__ |
|
284 # XXX method. This has had mixed success. Examples from both |
|
285 # XXX sides: sets work fine, and should become officially supported; |
|
286 # XXX dicts are much harder, and have failed in various subtle |
|
287 # XXX ways across attempts. Support for mapping types should probably |
|
288 # XXX be dropped (and users should pass mapping.keys() or .values() |
|
289 # XXX explicitly). |
|
290 |
|
291 # Sampling without replacement entails tracking either potential |
|
292 # selections (the pool) in a list or previous selections in a set. |
|
293 |
|
294 # When the number of selections is small compared to the |
|
295 # population, then tracking selections is efficient, requiring |
|
296 # only a small set and an occasional reselection. For |
|
297 # a larger number of selections, the pool tracking method is |
|
298 # preferred since the list takes less space than the |
|
299 # set and it doesn't suffer from frequent reselections. |
|
300 |
|
301 n = len(population) |
|
302 if not 0 <= k <= n: |
|
303 raise ValueError, "sample larger than population" |
|
304 random = self.random |
|
305 _int = int |
|
306 result = [None] * k |
|
307 setsize = 21 # size of a small set minus size of an empty list |
|
308 if k > 5: |
|
309 setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets |
|
310 if n <= setsize or hasattr(population, "keys"): |
|
311 # An n-length list is smaller than a k-length set, or this is a |
|
312 # mapping type so the other algorithm wouldn't work. |
|
313 pool = list(population) |
|
314 for i in xrange(k): # invariant: non-selected at [0,n-i) |
|
315 j = _int(random() * (n-i)) |
|
316 result[i] = pool[j] |
|
317 pool[j] = pool[n-i-1] # move non-selected item into vacancy |
|
318 else: |
|
319 try: |
|
320 selected = set() |
|
321 selected_add = selected.add |
|
322 for i in xrange(k): |
|
323 j = _int(random() * n) |
|
324 while j in selected: |
|
325 j = _int(random() * n) |
|
326 selected_add(j) |
|
327 result[i] = population[j] |
|
328 except (TypeError, KeyError): # handle (at least) sets |
|
329 if isinstance(population, list): |
|
330 raise |
|
331 return self.sample(tuple(population), k) |
|
332 return result |
|
333 |
|
334 ## -------------------- real-valued distributions ------------------- |
|
335 |
|
336 ## -------------------- uniform distribution ------------------- |
|
337 |
|
338 def uniform(self, a, b): |
|
339 """Get a random number in the range [a, b).""" |
|
340 return a + (b-a) * self.random() |
|
341 |
|
342 ## -------------------- normal distribution -------------------- |
|
343 |
|
344 def normalvariate(self, mu, sigma): |
|
345 """Normal distribution. |
|
346 |
|
347 mu is the mean, and sigma is the standard deviation. |
|
348 |
|
349 """ |
|
350 # mu = mean, sigma = standard deviation |
|
351 |
|
352 # Uses Kinderman and Monahan method. Reference: Kinderman, |
|
353 # A.J. and Monahan, J.F., "Computer generation of random |
|
354 # variables using the ratio of uniform deviates", ACM Trans |
|
355 # Math Software, 3, (1977), pp257-260. |
|
356 |
|
357 random = self.random |
|
358 while 1: |
|
359 u1 = random() |
|
360 u2 = 1.0 - random() |
|
361 z = NV_MAGICCONST*(u1-0.5)/u2 |
|
362 zz = z*z/4.0 |
|
363 if zz <= -_log(u2): |
|
364 break |
|
365 return mu + z*sigma |
|
366 |
|
367 ## -------------------- lognormal distribution -------------------- |
|
368 |
|
369 def lognormvariate(self, mu, sigma): |
|
370 """Log normal distribution. |
|
371 |
|
372 If you take the natural logarithm of this distribution, you'll get a |
|
373 normal distribution with mean mu and standard deviation sigma. |
|
374 mu can have any value, and sigma must be greater than zero. |
|
375 |
|
376 """ |
|
377 return _exp(self.normalvariate(mu, sigma)) |
|
378 |
|
379 ## -------------------- exponential distribution -------------------- |
|
380 |
|
381 def expovariate(self, lambd): |
|
382 """Exponential distribution. |
|
383 |
|
384 lambd is 1.0 divided by the desired mean. (The parameter would be |
|
385 called "lambda", but that is a reserved word in Python.) Returned |
|
386 values range from 0 to positive infinity. |
|
387 |
|
388 """ |
|
389 # lambd: rate lambd = 1/mean |
|
390 # ('lambda' is a Python reserved word) |
|
391 |
|
392 random = self.random |
|
393 u = random() |
|
394 while u <= 1e-7: |
|
395 u = random() |
|
396 return -_log(u)/lambd |
|
397 |
|
398 ## -------------------- von Mises distribution -------------------- |
|
399 |
|
400 def vonmisesvariate(self, mu, kappa): |
|
401 """Circular data distribution. |
|
402 |
|
403 mu is the mean angle, expressed in radians between 0 and 2*pi, and |
|
404 kappa is the concentration parameter, which must be greater than or |
|
405 equal to zero. If kappa is equal to zero, this distribution reduces |
|
406 to a uniform random angle over the range 0 to 2*pi. |
|
407 |
|
408 """ |
|
409 # mu: mean angle (in radians between 0 and 2*pi) |
|
410 # kappa: concentration parameter kappa (>= 0) |
|
411 # if kappa = 0 generate uniform random angle |
|
412 |
|
413 # Based upon an algorithm published in: Fisher, N.I., |
|
414 # "Statistical Analysis of Circular Data", Cambridge |
|
415 # University Press, 1993. |
|
416 |
|
417 # Thanks to Magnus Kessler for a correction to the |
|
418 # implementation of step 4. |
|
419 |
|
420 random = self.random |
|
421 if kappa <= 1e-6: |
|
422 return TWOPI * random() |
|
423 |
|
424 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) |
|
425 b = (a - _sqrt(2.0 * a))/(2.0 * kappa) |
|
426 r = (1.0 + b * b)/(2.0 * b) |
|
427 |
|
428 while 1: |
|
429 u1 = random() |
|
430 |
|
431 z = _cos(_pi * u1) |
|
432 f = (1.0 + r * z)/(r + z) |
|
433 c = kappa * (r - f) |
|
434 |
|
435 u2 = random() |
|
436 |
|
437 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c): |
|
438 break |
|
439 |
|
440 u3 = random() |
|
441 if u3 > 0.5: |
|
442 theta = (mu % TWOPI) + _acos(f) |
|
443 else: |
|
444 theta = (mu % TWOPI) - _acos(f) |
|
445 |
|
446 return theta |
|
447 |
|
448 ## -------------------- gamma distribution -------------------- |
|
449 |
|
450 def gammavariate(self, alpha, beta): |
|
451 """Gamma distribution. Not the gamma function! |
|
452 |
|
453 Conditions on the parameters are alpha > 0 and beta > 0. |
|
454 |
|
455 """ |
|
456 |
|
457 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 |
|
458 |
|
459 # Warning: a few older sources define the gamma distribution in terms |
|
460 # of alpha > -1.0 |
|
461 if alpha <= 0.0 or beta <= 0.0: |
|
462 raise ValueError, 'gammavariate: alpha and beta must be > 0.0' |
|
463 |
|
464 random = self.random |
|
465 if alpha > 1.0: |
|
466 |
|
467 # Uses R.C.H. Cheng, "The generation of Gamma |
|
468 # variables with non-integral shape parameters", |
|
469 # Applied Statistics, (1977), 26, No. 1, p71-74 |
|
470 |
|
471 ainv = _sqrt(2.0 * alpha - 1.0) |
|
472 bbb = alpha - LOG4 |
|
473 ccc = alpha + ainv |
|
474 |
|
475 while 1: |
|
476 u1 = random() |
|
477 if not 1e-7 < u1 < .9999999: |
|
478 continue |
|
479 u2 = 1.0 - random() |
|
480 v = _log(u1/(1.0-u1))/ainv |
|
481 x = alpha*_exp(v) |
|
482 z = u1*u1*u2 |
|
483 r = bbb+ccc*v-x |
|
484 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): |
|
485 return x * beta |
|
486 |
|
487 elif alpha == 1.0: |
|
488 # expovariate(1) |
|
489 u = random() |
|
490 while u <= 1e-7: |
|
491 u = random() |
|
492 return -_log(u) * beta |
|
493 |
|
494 else: # alpha is between 0 and 1 (exclusive) |
|
495 |
|
496 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle |
|
497 |
|
498 while 1: |
|
499 u = random() |
|
500 b = (_e + alpha)/_e |
|
501 p = b*u |
|
502 if p <= 1.0: |
|
503 x = p ** (1.0/alpha) |
|
504 else: |
|
505 x = -_log((b-p)/alpha) |
|
506 u1 = random() |
|
507 if p > 1.0: |
|
508 if u1 <= x ** (alpha - 1.0): |
|
509 break |
|
510 elif u1 <= _exp(-x): |
|
511 break |
|
512 return x * beta |
|
513 |
|
514 ## -------------------- Gauss (faster alternative) -------------------- |
|
515 |
|
516 def gauss(self, mu, sigma): |
|
517 """Gaussian distribution. |
|
518 |
|
519 mu is the mean, and sigma is the standard deviation. This is |
|
520 slightly faster than the normalvariate() function. |
|
521 |
|
522 Not thread-safe without a lock around calls. |
|
523 |
|
524 """ |
|
525 |
|
526 # When x and y are two variables from [0, 1), uniformly |
|
527 # distributed, then |
|
528 # |
|
529 # cos(2*pi*x)*sqrt(-2*log(1-y)) |
|
530 # sin(2*pi*x)*sqrt(-2*log(1-y)) |
|
531 # |
|
532 # are two *independent* variables with normal distribution |
|
533 # (mu = 0, sigma = 1). |
|
534 # (Lambert Meertens) |
|
535 # (corrected version; bug discovered by Mike Miller, fixed by LM) |
|
536 |
|
537 # Multithreading note: When two threads call this function |
|
538 # simultaneously, it is possible that they will receive the |
|
539 # same return value. The window is very small though. To |
|
540 # avoid this, you have to use a lock around all calls. (I |
|
541 # didn't want to slow this down in the serial case by using a |
|
542 # lock here.) |
|
543 |
|
544 random = self.random |
|
545 z = self.gauss_next |
|
546 self.gauss_next = None |
|
547 if z is None: |
|
548 x2pi = random() * TWOPI |
|
549 g2rad = _sqrt(-2.0 * _log(1.0 - random())) |
|
550 z = _cos(x2pi) * g2rad |
|
551 self.gauss_next = _sin(x2pi) * g2rad |
|
552 |
|
553 return mu + z*sigma |
|
554 |
|
555 ## -------------------- beta -------------------- |
|
556 ## See |
|
557 ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 |
|
558 ## for Ivan Frohne's insightful analysis of why the original implementation: |
|
559 ## |
|
560 ## def betavariate(self, alpha, beta): |
|
561 ## # Discrete Event Simulation in C, pp 87-88. |
|
562 ## |
|
563 ## y = self.expovariate(alpha) |
|
564 ## z = self.expovariate(1.0/beta) |
|
565 ## return z/(y+z) |
|
566 ## |
|
567 ## was dead wrong, and how it probably got that way. |
|
568 |
|
569 def betavariate(self, alpha, beta): |
|
570 """Beta distribution. |
|
571 |
|
572 Conditions on the parameters are alpha > 0 and beta > 0. |
|
573 Returned values range between 0 and 1. |
|
574 |
|
575 """ |
|
576 |
|
577 # This version due to Janne Sinkkonen, and matches all the std |
|
578 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). |
|
579 y = self.gammavariate(alpha, 1.) |
|
580 if y == 0: |
|
581 return 0.0 |
|
582 else: |
|
583 return y / (y + self.gammavariate(beta, 1.)) |
|
584 |
|
585 ## -------------------- Pareto -------------------- |
|
586 |
|
587 def paretovariate(self, alpha): |
|
588 """Pareto distribution. alpha is the shape parameter.""" |
|
589 # Jain, pg. 495 |
|
590 |
|
591 u = 1.0 - self.random() |
|
592 return 1.0 / pow(u, 1.0/alpha) |
|
593 |
|
594 ## -------------------- Weibull -------------------- |
|
595 |
|
596 def weibullvariate(self, alpha, beta): |
|
597 """Weibull distribution. |
|
598 |
|
599 alpha is the scale parameter and beta is the shape parameter. |
|
600 |
|
601 """ |
|
602 # Jain, pg. 499; bug fix courtesy Bill Arms |
|
603 |
|
604 u = 1.0 - self.random() |
|
605 return alpha * pow(-_log(u), 1.0/beta) |
|
606 |
|
607 ## -------------------- Wichmann-Hill ------------------- |
|
608 |
|
609 class WichmannHill(Random): |
|
610 |
|
611 VERSION = 1 # used by getstate/setstate |
|
612 |
|
613 def seed(self, a=None): |
|
614 """Initialize internal state from hashable object. |
|
615 |
|
616 None or no argument seeds from current time or from an operating |
|
617 system specific randomness source if available. |
|
618 |
|
619 If a is not None or an int or long, hash(a) is used instead. |
|
620 |
|
621 If a is an int or long, a is used directly. Distinct values between |
|
622 0 and 27814431486575L inclusive are guaranteed to yield distinct |
|
623 internal states (this guarantee is specific to the default |
|
624 Wichmann-Hill generator). |
|
625 """ |
|
626 |
|
627 if a is None: |
|
628 try: |
|
629 a = long(_hexlify(_urandom(16)), 16) |
|
630 except NotImplementedError: |
|
631 import time |
|
632 a = long(time.time() * 256) # use fractional seconds |
|
633 |
|
634 if not isinstance(a, (int, long)): |
|
635 a = hash(a) |
|
636 |
|
637 a, x = divmod(a, 30268) |
|
638 a, y = divmod(a, 30306) |
|
639 a, z = divmod(a, 30322) |
|
640 self._seed = int(x)+1, int(y)+1, int(z)+1 |
|
641 |
|
642 self.gauss_next = None |
|
643 |
|
644 def random(self): |
|
645 """Get the next random number in the range [0.0, 1.0).""" |
|
646 |
|
647 # Wichman-Hill random number generator. |
|
648 # |
|
649 # Wichmann, B. A. & Hill, I. D. (1982) |
|
650 # Algorithm AS 183: |
|
651 # An efficient and portable pseudo-random number generator |
|
652 # Applied Statistics 31 (1982) 188-190 |
|
653 # |
|
654 # see also: |
|
655 # Correction to Algorithm AS 183 |
|
656 # Applied Statistics 33 (1984) 123 |
|
657 # |
|
658 # McLeod, A. I. (1985) |
|
659 # A remark on Algorithm AS 183 |
|
660 # Applied Statistics 34 (1985),198-200 |
|
661 |
|
662 # This part is thread-unsafe: |
|
663 # BEGIN CRITICAL SECTION |
|
664 x, y, z = self._seed |
|
665 x = (171 * x) % 30269 |
|
666 y = (172 * y) % 30307 |
|
667 z = (170 * z) % 30323 |
|
668 self._seed = x, y, z |
|
669 # END CRITICAL SECTION |
|
670 |
|
671 # Note: on a platform using IEEE-754 double arithmetic, this can |
|
672 # never return 0.0 (asserted by Tim; proof too long for a comment). |
|
673 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 |
|
674 |
|
675 def getstate(self): |
|
676 """Return internal state; can be passed to setstate() later.""" |
|
677 return self.VERSION, self._seed, self.gauss_next |
|
678 |
|
679 def setstate(self, state): |
|
680 """Restore internal state from object returned by getstate().""" |
|
681 version = state[0] |
|
682 if version == 1: |
|
683 version, self._seed, self.gauss_next = state |
|
684 else: |
|
685 raise ValueError("state with version %s passed to " |
|
686 "Random.setstate() of version %s" % |
|
687 (version, self.VERSION)) |
|
688 |
|
689 def jumpahead(self, n): |
|
690 """Act as if n calls to random() were made, but quickly. |
|
691 |
|
692 n is an int, greater than or equal to 0. |
|
693 |
|
694 Example use: If you have 2 threads and know that each will |
|
695 consume no more than a million random numbers, create two Random |
|
696 objects r1 and r2, then do |
|
697 r2.setstate(r1.getstate()) |
|
698 r2.jumpahead(1000000) |
|
699 Then r1 and r2 will use guaranteed-disjoint segments of the full |
|
700 period. |
|
701 """ |
|
702 |
|
703 if not n >= 0: |
|
704 raise ValueError("n must be >= 0") |
|
705 x, y, z = self._seed |
|
706 x = int(x * pow(171, n, 30269)) % 30269 |
|
707 y = int(y * pow(172, n, 30307)) % 30307 |
|
708 z = int(z * pow(170, n, 30323)) % 30323 |
|
709 self._seed = x, y, z |
|
710 |
|
711 def __whseed(self, x=0, y=0, z=0): |
|
712 """Set the Wichmann-Hill seed from (x, y, z). |
|
713 |
|
714 These must be integers in the range [0, 256). |
|
715 """ |
|
716 |
|
717 if not type(x) == type(y) == type(z) == int: |
|
718 raise TypeError('seeds must be integers') |
|
719 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): |
|
720 raise ValueError('seeds must be in range(0, 256)') |
|
721 if 0 == x == y == z: |
|
722 # Initialize from current time |
|
723 import time |
|
724 t = long(time.time() * 256) |
|
725 t = int((t&0xffffff) ^ (t>>24)) |
|
726 t, x = divmod(t, 256) |
|
727 t, y = divmod(t, 256) |
|
728 t, z = divmod(t, 256) |
|
729 # Zero is a poor seed, so substitute 1 |
|
730 self._seed = (x or 1, y or 1, z or 1) |
|
731 |
|
732 self.gauss_next = None |
|
733 |
|
734 def whseed(self, a=None): |
|
735 """Seed from hashable object's hash code. |
|
736 |
|
737 None or no argument seeds from current time. It is not guaranteed |
|
738 that objects with distinct hash codes lead to distinct internal |
|
739 states. |
|
740 |
|
741 This is obsolete, provided for compatibility with the seed routine |
|
742 used prior to Python 2.1. Use the .seed() method instead. |
|
743 """ |
|
744 |
|
745 if a is None: |
|
746 self.__whseed() |
|
747 return |
|
748 a = hash(a) |
|
749 a, x = divmod(a, 256) |
|
750 a, y = divmod(a, 256) |
|
751 a, z = divmod(a, 256) |
|
752 x = (x + a) % 256 or 1 |
|
753 y = (y + a) % 256 or 1 |
|
754 z = (z + a) % 256 or 1 |
|
755 self.__whseed(x, y, z) |
|
756 |
|
757 ## --------------- Operating System Random Source ------------------ |
|
758 |
|
759 class SystemRandom(Random): |
|
760 """Alternate random number generator using sources provided |
|
761 by the operating system (such as /dev/urandom on Unix or |
|
762 CryptGenRandom on Windows). |
|
763 |
|
764 Not available on all systems (see os.urandom() for details). |
|
765 """ |
|
766 |
|
767 def random(self): |
|
768 """Get the next random number in the range [0.0, 1.0).""" |
|
769 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF |
|
770 |
|
771 def getrandbits(self, k): |
|
772 """getrandbits(k) -> x. Generates a long int with k random bits.""" |
|
773 if k <= 0: |
|
774 raise ValueError('number of bits must be greater than zero') |
|
775 if k != int(k): |
|
776 raise TypeError('number of bits should be an integer') |
|
777 bytes = (k + 7) // 8 # bits / 8 and rounded up |
|
778 x = long(_hexlify(_urandom(bytes)), 16) |
|
779 return x >> (bytes * 8 - k) # trim excess bits |
|
780 |
|
781 def _stub(self, *args, **kwds): |
|
782 "Stub method. Not used for a system random number generator." |
|
783 return None |
|
784 seed = jumpahead = _stub |
|
785 |
|
786 def _notimplemented(self, *args, **kwds): |
|
787 "Method should not be called for a system random number generator." |
|
788 raise NotImplementedError('System entropy source does not have state.') |
|
789 getstate = setstate = _notimplemented |
|
790 |
|
791 ## -------------------- test program -------------------- |
|
792 |
|
793 def _test_generator(n, func, args): |
|
794 import time |
|
795 print n, 'times', func.__name__ |
|
796 total = 0.0 |
|
797 sqsum = 0.0 |
|
798 smallest = 1e10 |
|
799 largest = -1e10 |
|
800 t0 = time.time() |
|
801 for i in range(n): |
|
802 x = func(*args) |
|
803 total += x |
|
804 sqsum = sqsum + x*x |
|
805 smallest = min(x, smallest) |
|
806 largest = max(x, largest) |
|
807 t1 = time.time() |
|
808 print round(t1-t0, 3), 'sec,', |
|
809 avg = total/n |
|
810 stddev = _sqrt(sqsum/n - avg*avg) |
|
811 print 'avg %g, stddev %g, min %g, max %g' % \ |
|
812 (avg, stddev, smallest, largest) |
|
813 |
|
814 |
|
815 def _test(N=2000): |
|
816 _test_generator(N, random, ()) |
|
817 _test_generator(N, normalvariate, (0.0, 1.0)) |
|
818 _test_generator(N, lognormvariate, (0.0, 1.0)) |
|
819 _test_generator(N, vonmisesvariate, (0.0, 1.0)) |
|
820 _test_generator(N, gammavariate, (0.01, 1.0)) |
|
821 _test_generator(N, gammavariate, (0.1, 1.0)) |
|
822 _test_generator(N, gammavariate, (0.1, 2.0)) |
|
823 _test_generator(N, gammavariate, (0.5, 1.0)) |
|
824 _test_generator(N, gammavariate, (0.9, 1.0)) |
|
825 _test_generator(N, gammavariate, (1.0, 1.0)) |
|
826 _test_generator(N, gammavariate, (2.0, 1.0)) |
|
827 _test_generator(N, gammavariate, (20.0, 1.0)) |
|
828 _test_generator(N, gammavariate, (200.0, 1.0)) |
|
829 _test_generator(N, gauss, (0.0, 1.0)) |
|
830 _test_generator(N, betavariate, (3.0, 3.0)) |
|
831 |
|
832 # Create one instance, seeded from current time, and export its methods |
|
833 # as module-level functions. The functions share state across all uses |
|
834 #(both in the user's code and in the Python libraries), but that's fine |
|
835 # for most programs and is easier for the casual user than making them |
|
836 # instantiate their own Random() instance. |
|
837 |
|
838 _inst = Random() |
|
839 seed = _inst.seed |
|
840 random = _inst.random |
|
841 uniform = _inst.uniform |
|
842 randint = _inst.randint |
|
843 choice = _inst.choice |
|
844 randrange = _inst.randrange |
|
845 sample = _inst.sample |
|
846 shuffle = _inst.shuffle |
|
847 normalvariate = _inst.normalvariate |
|
848 lognormvariate = _inst.lognormvariate |
|
849 expovariate = _inst.expovariate |
|
850 vonmisesvariate = _inst.vonmisesvariate |
|
851 gammavariate = _inst.gammavariate |
|
852 gauss = _inst.gauss |
|
853 betavariate = _inst.betavariate |
|
854 paretovariate = _inst.paretovariate |
|
855 weibullvariate = _inst.weibullvariate |
|
856 getstate = _inst.getstate |
|
857 setstate = _inst.setstate |
|
858 jumpahead = _inst.jumpahead |
|
859 getrandbits = _inst.getrandbits |
|
860 |
|
861 if __name__ == '__main__': |
|
862 _test() |