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