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