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