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