Guido van Rossum | e7b146f | 2000-02-04 15:28:42 +0000 | [diff] [blame] | 1 | """Random variable generators. |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 2 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 3 | integers |
| 4 | -------- |
| 5 | uniform within range |
| 6 | |
| 7 | sequences |
| 8 | --------- |
| 9 | pick random element |
| 10 | generate random permutation |
| 11 | |
Guido van Rossum | e7b146f | 2000-02-04 15:28:42 +0000 | [diff] [blame] | 12 | distributions on the real line: |
| 13 | ------------------------------ |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 14 | uniform |
Guido van Rossum | e7b146f | 2000-02-04 15:28:42 +0000 | [diff] [blame] | 15 | normal (Gaussian) |
| 16 | lognormal |
| 17 | negative exponential |
| 18 | gamma |
| 19 | beta |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 20 | |
Guido van Rossum | e7b146f | 2000-02-04 15:28:42 +0000 | [diff] [blame] | 21 | distributions on the circle (angles 0 to 2pi) |
| 22 | --------------------------------------------- |
| 23 | circular uniform |
| 24 | von Mises |
| 25 | |
| 26 | Translated from anonymously contributed C/C++ source. |
| 27 | |
Tim Peters | e360d95 | 2001-01-26 10:00:39 +0000 | [diff] [blame^] | 28 | Multi-threading note: the random number generator used here is not thread- |
| 29 | safe; it is possible that two calls return the same random value. However, |
| 30 | you can instantiate a different instance of Random() in each thread to get |
| 31 | generators that don't share state, then use .setstate() and .jumpahead() to |
| 32 | move the generators to disjoint segments of the full period. For example, |
| 33 | |
| 34 | def create_generators(num, delta, firstseed=None): |
| 35 | ""\"Return list of num distinct generators. |
| 36 | Each generator has its own unique segment of delta elements from |
| 37 | Random.random()'s full period. |
| 38 | Seed the first generator with optional arg firstseed (default is |
| 39 | None, to seed from current time). |
| 40 | ""\" |
| 41 | |
| 42 | from random import Random |
| 43 | g = Random(firstseed) |
| 44 | result = [g] |
| 45 | for i in range(num - 1): |
| 46 | laststate = g.getstate() |
| 47 | g = Random() |
| 48 | g.setstate(laststate) |
| 49 | g.jumpahead(delta) |
| 50 | result.append(g) |
| 51 | return result |
| 52 | |
| 53 | gens = create_generators(10, 1000000) |
| 54 | |
| 55 | That creates 10 distinct generators, which can be passed out to 10 distinct |
| 56 | threads. The generators don't share state so can be called safely in |
| 57 | parallel. So long as no thread calls its g.random() more than a million |
| 58 | times (the second argument to create_generators), the sequences seen by |
| 59 | each thread will not overlap. |
| 60 | |
| 61 | The period of the underlying Wichmann-Hill generator is 6,953,607,871,644, |
| 62 | and that limits how far this technique can be pushed. |
| 63 | |
| 64 | Just for fun, note that since we know the period, .jumpahead() can also be |
| 65 | used to "move backward in time": |
| 66 | |
| 67 | >>> g = Random(42) # arbitrary |
| 68 | >>> g.random() |
| 69 | 0.24855401895528142 |
| 70 | >>> g.jumpahead(6953607871644L - 1) # move *back* one |
| 71 | >>> g.random() |
| 72 | 0.24855401895528142 |
Guido van Rossum | e7b146f | 2000-02-04 15:28:42 +0000 | [diff] [blame] | 73 | """ |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 74 | # XXX The docstring sucks. |
Guido van Rossum | d03e119 | 1998-05-29 17:51:31 +0000 | [diff] [blame] | 75 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 76 | from math import log as _log, exp as _exp, pi as _pi, e as _e |
| 77 | from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 78 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 79 | def _verify(name, expected): |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 80 | computed = eval(name) |
| 81 | if abs(computed - expected) > 1e-7: |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 82 | raise ValueError( |
| 83 | "computed value for %s deviates too much " |
| 84 | "(computed %g, expected %g)" % (name, computed, expected)) |
| 85 | |
| 86 | NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0) |
| 87 | _verify('NV_MAGICCONST', 1.71552776992141) |
| 88 | |
| 89 | TWOPI = 2.0*_pi |
| 90 | _verify('TWOPI', 6.28318530718) |
| 91 | |
| 92 | LOG4 = _log(4.0) |
| 93 | _verify('LOG4', 1.38629436111989) |
| 94 | |
| 95 | SG_MAGICCONST = 1.0 + _log(4.5) |
| 96 | _verify('SG_MAGICCONST', 2.50407739677627) |
| 97 | |
| 98 | del _verify |
| 99 | |
| 100 | # Translated by Guido van Rossum from C source provided by |
| 101 | # Adrian Baddeley. |
| 102 | |
| 103 | class Random: |
| 104 | |
| 105 | VERSION = 1 # used by getstate/setstate |
| 106 | |
| 107 | def __init__(self, x=None): |
| 108 | """Initialize an instance. |
| 109 | |
| 110 | Optional argument x controls seeding, as for Random.seed(). |
| 111 | """ |
| 112 | |
| 113 | self.seed(x) |
| 114 | self.gauss_next = None |
| 115 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 116 | ## -------------------- core generator ------------------- |
| 117 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 118 | # Specific to Wichmann-Hill generator. Subclasses wishing to use a |
Tim Peters | d52269b | 2001-01-25 06:23:18 +0000 | [diff] [blame] | 119 | # different core generator should override the seed(), random(), |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 120 | # getstate(), setstate() and jumpahead() methods. |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 121 | |
| 122 | def __whseed(self, x=0, y=0, z=0): |
| 123 | """Set the Wichmann-Hill seed from (x, y, z). |
| 124 | |
| 125 | These must be integers in the range [0, 256). |
| 126 | """ |
| 127 | |
| 128 | if not type(x) == type(y) == type(z) == type(0): |
| 129 | raise TypeError('seeds must be integers') |
| 130 | if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256): |
| 131 | raise ValueError('seeds must be in range(0, 256)') |
| 132 | if 0 == x == y == z: |
| 133 | # Initialize from current time |
| 134 | import time |
| 135 | t = long(time.time()) * 256 |
| 136 | t = int((t&0xffffff) ^ (t>>24)) |
| 137 | t, x = divmod(t, 256) |
| 138 | t, y = divmod(t, 256) |
| 139 | t, z = divmod(t, 256) |
| 140 | # Zero is a poor seed, so substitute 1 |
| 141 | self._seed = (x or 1, y or 1, z or 1) |
| 142 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 143 | def random(self): |
| 144 | """Get the next random number in the range [0.0, 1.0).""" |
| 145 | |
| 146 | # Wichman-Hill random number generator. |
| 147 | # |
| 148 | # Wichmann, B. A. & Hill, I. D. (1982) |
| 149 | # Algorithm AS 183: |
| 150 | # An efficient and portable pseudo-random number generator |
| 151 | # Applied Statistics 31 (1982) 188-190 |
| 152 | # |
| 153 | # see also: |
| 154 | # Correction to Algorithm AS 183 |
| 155 | # Applied Statistics 33 (1984) 123 |
| 156 | # |
| 157 | # McLeod, A. I. (1985) |
| 158 | # A remark on Algorithm AS 183 |
| 159 | # Applied Statistics 34 (1985),198-200 |
| 160 | |
| 161 | # This part is thread-unsafe: |
| 162 | # BEGIN CRITICAL SECTION |
| 163 | x, y, z = self._seed |
| 164 | x = (171 * x) % 30269 |
| 165 | y = (172 * y) % 30307 |
| 166 | z = (170 * z) % 30323 |
| 167 | self._seed = x, y, z |
| 168 | # END CRITICAL SECTION |
| 169 | |
| 170 | # Note: on a platform using IEEE-754 double arithmetic, this can |
| 171 | # never return 0.0 (asserted by Tim; proof too long for a comment). |
| 172 | return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0 |
| 173 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 174 | def seed(self, a=None): |
| 175 | """Seed from hashable object's hash code. |
| 176 | |
| 177 | None or no argument seeds from current time. It is not guaranteed |
| 178 | that objects with distinct hash codes lead to distinct internal |
| 179 | states. |
| 180 | """ |
| 181 | |
| 182 | if a is None: |
| 183 | self.__whseed() |
| 184 | return |
| 185 | a = hash(a) |
| 186 | a, x = divmod(a, 256) |
| 187 | a, y = divmod(a, 256) |
| 188 | a, z = divmod(a, 256) |
| 189 | x = (x + a) % 256 or 1 |
| 190 | y = (y + a) % 256 or 1 |
| 191 | z = (z + a) % 256 or 1 |
| 192 | self.__whseed(x, y, z) |
| 193 | |
| 194 | def getstate(self): |
| 195 | """Return internal state; can be passed to setstate() later.""" |
| 196 | return self.VERSION, self._seed, self.gauss_next |
| 197 | |
| 198 | def setstate(self, state): |
| 199 | """Restore internal state from object returned by getstate().""" |
| 200 | version = state[0] |
| 201 | if version == 1: |
| 202 | version, self._seed, self.gauss_next = state |
| 203 | else: |
| 204 | raise ValueError("state with version %s passed to " |
| 205 | "Random.setstate() of version %s" % |
| 206 | (version, self.VERSION)) |
| 207 | |
| 208 | def jumpahead(self, n): |
| 209 | """Act as if n calls to random() were made, but quickly. |
| 210 | |
| 211 | n is an int, greater than or equal to 0. |
| 212 | |
| 213 | Example use: If you have 2 threads and know that each will |
| 214 | consume no more than a million random numbers, create two Random |
| 215 | objects r1 and r2, then do |
| 216 | r2.setstate(r1.getstate()) |
| 217 | r2.jumpahead(1000000) |
| 218 | Then r1 and r2 will use guaranteed-disjoint segments of the full |
| 219 | period. |
| 220 | """ |
| 221 | |
| 222 | if not n >= 0: |
| 223 | raise ValueError("n must be >= 0") |
| 224 | x, y, z = self._seed |
| 225 | x = int(x * pow(171, n, 30269)) % 30269 |
| 226 | y = int(y * pow(172, n, 30307)) % 30307 |
| 227 | z = int(z * pow(170, n, 30323)) % 30323 |
| 228 | self._seed = x, y, z |
| 229 | |
| 230 | ## ---- Methods below this point do not need to be overridden when |
| 231 | ## ---- subclassing for the purpose of using a different core generator. |
| 232 | |
| 233 | ## -------------------- pickle support ------------------- |
| 234 | |
| 235 | def __getstate__(self): # for pickle |
| 236 | return self.getstate() |
| 237 | |
| 238 | def __setstate__(self, state): # for pickle |
| 239 | self.setstate(state) |
| 240 | |
| 241 | ## -------------------- integer methods ------------------- |
| 242 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 243 | def randrange(self, start, stop=None, step=1, int=int, default=None): |
| 244 | """Choose a random item from range(start, stop[, step]). |
| 245 | |
| 246 | This fixes the problem with randint() which includes the |
| 247 | endpoint; in Python this is usually not what you want. |
| 248 | Do not supply the 'int' and 'default' arguments. |
| 249 | """ |
| 250 | |
| 251 | # This code is a bit messy to make it fast for the |
| 252 | # common case while still doing adequate error checking |
| 253 | istart = int(start) |
| 254 | if istart != start: |
| 255 | raise ValueError, "non-integer arg 1 for randrange()" |
| 256 | if stop is default: |
| 257 | if istart > 0: |
| 258 | return int(self.random() * istart) |
| 259 | raise ValueError, "empty range for randrange()" |
| 260 | istop = int(stop) |
| 261 | if istop != stop: |
| 262 | raise ValueError, "non-integer stop for randrange()" |
| 263 | if step == 1: |
| 264 | if istart < istop: |
| 265 | return istart + int(self.random() * |
| 266 | (istop - istart)) |
| 267 | raise ValueError, "empty range for randrange()" |
| 268 | istep = int(step) |
| 269 | if istep != step: |
| 270 | raise ValueError, "non-integer step for randrange()" |
| 271 | if istep > 0: |
| 272 | n = (istop - istart + istep - 1) / istep |
| 273 | elif istep < 0: |
| 274 | n = (istop - istart + istep + 1) / istep |
| 275 | else: |
| 276 | raise ValueError, "zero step for randrange()" |
| 277 | |
| 278 | if n <= 0: |
| 279 | raise ValueError, "empty range for randrange()" |
| 280 | return istart + istep*int(self.random() * n) |
| 281 | |
| 282 | def randint(self, a, b): |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 283 | """Return random integer in range [a, b], including both end points. |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 284 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 285 | (Deprecated; use randrange(a, b+1).) |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 286 | """ |
| 287 | |
| 288 | return self.randrange(a, b+1) |
| 289 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 290 | ## -------------------- sequence methods ------------------- |
| 291 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 292 | def choice(self, seq): |
| 293 | """Choose a random element from a non-empty sequence.""" |
| 294 | return seq[int(self.random() * len(seq))] |
| 295 | |
| 296 | def shuffle(self, x, random=None, int=int): |
| 297 | """x, random=random.random -> shuffle list x in place; return None. |
| 298 | |
| 299 | Optional arg random is a 0-argument function returning a random |
| 300 | float in [0.0, 1.0); by default, the standard random.random. |
| 301 | |
| 302 | Note that for even rather small len(x), the total number of |
| 303 | permutations of x is larger than the period of most random number |
| 304 | generators; this implies that "most" permutations of a long |
| 305 | sequence can never be generated. |
| 306 | """ |
| 307 | |
| 308 | if random is None: |
| 309 | random = self.random |
| 310 | for i in xrange(len(x)-1, 0, -1): |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 311 | # pick an element in x[:i+1] with which to exchange x[i] |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 312 | j = int(random() * (i+1)) |
| 313 | x[i], x[j] = x[j], x[i] |
| 314 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 315 | ## -------------------- real-valued distributions ------------------- |
| 316 | |
| 317 | ## -------------------- uniform distribution ------------------- |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 318 | |
| 319 | def uniform(self, a, b): |
| 320 | """Get a random number in the range [a, b).""" |
| 321 | return a + (b-a) * self.random() |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 322 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 323 | ## -------------------- normal distribution -------------------- |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 324 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 325 | def normalvariate(self, mu, sigma): |
| 326 | # mu = mean, sigma = standard deviation |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 327 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 328 | # Uses Kinderman and Monahan method. Reference: Kinderman, |
| 329 | # A.J. and Monahan, J.F., "Computer generation of random |
| 330 | # variables using the ratio of uniform deviates", ACM Trans |
| 331 | # Math Software, 3, (1977), pp257-260. |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 332 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 333 | random = self.random |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 334 | while 1: |
| 335 | u1 = random() |
| 336 | u2 = random() |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 337 | z = NV_MAGICCONST*(u1-0.5)/u2 |
| 338 | zz = z*z/4.0 |
| 339 | if zz <= -_log(u2): |
| 340 | break |
| 341 | return mu + z*sigma |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 342 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 343 | ## -------------------- lognormal distribution -------------------- |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 344 | |
| 345 | def lognormvariate(self, mu, sigma): |
| 346 | return _exp(self.normalvariate(mu, sigma)) |
| 347 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 348 | ## -------------------- circular uniform -------------------- |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 349 | |
| 350 | def cunifvariate(self, mean, arc): |
| 351 | # mean: mean angle (in radians between 0 and pi) |
| 352 | # arc: range of distribution (in radians between 0 and pi) |
| 353 | |
| 354 | return (mean + arc * (self.random() - 0.5)) % _pi |
| 355 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 356 | ## -------------------- exponential distribution -------------------- |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 357 | |
| 358 | def expovariate(self, lambd): |
| 359 | # lambd: rate lambd = 1/mean |
| 360 | # ('lambda' is a Python reserved word) |
| 361 | |
| 362 | random = self.random |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 363 | u = random() |
| 364 | while u <= 1e-7: |
| 365 | u = random() |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 366 | return -_log(u)/lambd |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 367 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 368 | ## -------------------- von Mises distribution -------------------- |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 369 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 370 | def vonmisesvariate(self, mu, kappa): |
| 371 | # mu: mean angle (in radians between 0 and 2*pi) |
| 372 | # kappa: concentration parameter kappa (>= 0) |
| 373 | # if kappa = 0 generate uniform random angle |
| 374 | |
| 375 | # Based upon an algorithm published in: Fisher, N.I., |
| 376 | # "Statistical Analysis of Circular Data", Cambridge |
| 377 | # University Press, 1993. |
| 378 | |
| 379 | # Thanks to Magnus Kessler for a correction to the |
| 380 | # implementation of step 4. |
| 381 | |
| 382 | random = self.random |
| 383 | if kappa <= 1e-6: |
| 384 | return TWOPI * random() |
| 385 | |
| 386 | a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) |
| 387 | b = (a - _sqrt(2.0 * a))/(2.0 * kappa) |
| 388 | r = (1.0 + b * b)/(2.0 * b) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 389 | |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 390 | while 1: |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 391 | u1 = random() |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 392 | |
| 393 | z = _cos(_pi * u1) |
| 394 | f = (1.0 + r * z)/(r + z) |
| 395 | c = kappa * (r - f) |
| 396 | |
| 397 | u2 = random() |
| 398 | |
| 399 | if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)): |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 400 | break |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 401 | |
| 402 | u3 = random() |
| 403 | if u3 > 0.5: |
| 404 | theta = (mu % TWOPI) + _acos(f) |
| 405 | else: |
| 406 | theta = (mu % TWOPI) - _acos(f) |
| 407 | |
| 408 | return theta |
| 409 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 410 | ## -------------------- gamma distribution -------------------- |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 411 | |
| 412 | def gammavariate(self, alpha, beta): |
| 413 | # beta times standard gamma |
| 414 | ainv = _sqrt(2.0 * alpha - 1.0) |
| 415 | return beta * self.stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) |
| 416 | |
| 417 | def stdgamma(self, alpha, ainv, bbb, ccc): |
| 418 | # ainv = sqrt(2 * alpha - 1) |
| 419 | # bbb = alpha - log(4) |
| 420 | # ccc = alpha + ainv |
| 421 | |
| 422 | random = self.random |
| 423 | if alpha <= 0.0: |
| 424 | raise ValueError, 'stdgamma: alpha must be > 0.0' |
| 425 | |
| 426 | if alpha > 1.0: |
| 427 | |
| 428 | # Uses R.C.H. Cheng, "The generation of Gamma |
| 429 | # variables with non-integral shape parameters", |
| 430 | # Applied Statistics, (1977), 26, No. 1, p71-74 |
| 431 | |
| 432 | while 1: |
| 433 | u1 = random() |
| 434 | u2 = random() |
| 435 | v = _log(u1/(1.0-u1))/ainv |
| 436 | x = alpha*_exp(v) |
| 437 | z = u1*u1*u2 |
| 438 | r = bbb+ccc*v-x |
| 439 | if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): |
| 440 | return x |
| 441 | |
| 442 | elif alpha == 1.0: |
| 443 | # expovariate(1) |
| 444 | u = random() |
| 445 | while u <= 1e-7: |
| 446 | u = random() |
| 447 | return -_log(u) |
| 448 | |
| 449 | else: # alpha is between 0 and 1 (exclusive) |
| 450 | |
| 451 | # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle |
| 452 | |
| 453 | while 1: |
| 454 | u = random() |
| 455 | b = (_e + alpha)/_e |
| 456 | p = b*u |
| 457 | if p <= 1.0: |
| 458 | x = pow(p, 1.0/alpha) |
| 459 | else: |
| 460 | # p > 1 |
| 461 | x = -_log((b-p)/alpha) |
| 462 | u1 = random() |
| 463 | if not (((p <= 1.0) and (u1 > _exp(-x))) or |
| 464 | ((p > 1) and (u1 > pow(x, alpha - 1.0)))): |
| 465 | break |
| 466 | return x |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 467 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 468 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 469 | ## -------------------- Gauss (faster alternative) -------------------- |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 470 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 471 | def gauss(self, mu, sigma): |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 472 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 473 | # When x and y are two variables from [0, 1), uniformly |
| 474 | # distributed, then |
| 475 | # |
| 476 | # cos(2*pi*x)*sqrt(-2*log(1-y)) |
| 477 | # sin(2*pi*x)*sqrt(-2*log(1-y)) |
| 478 | # |
| 479 | # are two *independent* variables with normal distribution |
| 480 | # (mu = 0, sigma = 1). |
| 481 | # (Lambert Meertens) |
| 482 | # (corrected version; bug discovered by Mike Miller, fixed by LM) |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 483 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 484 | # Multithreading note: When two threads call this function |
| 485 | # simultaneously, it is possible that they will receive the |
| 486 | # same return value. The window is very small though. To |
| 487 | # avoid this, you have to use a lock around all calls. (I |
| 488 | # didn't want to slow this down in the serial case by using a |
| 489 | # lock here.) |
Guido van Rossum | d03e119 | 1998-05-29 17:51:31 +0000 | [diff] [blame] | 490 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 491 | random = self.random |
| 492 | z = self.gauss_next |
| 493 | self.gauss_next = None |
| 494 | if z is None: |
| 495 | x2pi = random() * TWOPI |
| 496 | g2rad = _sqrt(-2.0 * _log(1.0 - random())) |
| 497 | z = _cos(x2pi) * g2rad |
| 498 | self.gauss_next = _sin(x2pi) * g2rad |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 499 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 500 | return mu + z*sigma |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 501 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 502 | ## -------------------- beta -------------------- |
Tim Peters | 85e2e47 | 2001-01-26 06:49:56 +0000 | [diff] [blame] | 503 | ## See |
| 504 | ## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470 |
| 505 | ## for Ivan Frohne's insightful analysis of why the original implementation: |
| 506 | ## |
| 507 | ## def betavariate(self, alpha, beta): |
| 508 | ## # Discrete Event Simulation in C, pp 87-88. |
| 509 | ## |
| 510 | ## y = self.expovariate(alpha) |
| 511 | ## z = self.expovariate(1.0/beta) |
| 512 | ## return z/(y+z) |
| 513 | ## |
| 514 | ## was dead wrong, and how it probably got that way. |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 515 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 516 | def betavariate(self, alpha, beta): |
Tim Peters | 85e2e47 | 2001-01-26 06:49:56 +0000 | [diff] [blame] | 517 | # This version due to Janne Sinkkonen, and matches all the std |
| 518 | # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution"). |
| 519 | y = self.gammavariate(alpha, 1.) |
| 520 | if y == 0: |
| 521 | return 0.0 |
| 522 | else: |
| 523 | return y / (y + self.gammavariate(beta, 1.)) |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 524 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 525 | ## -------------------- Pareto -------------------- |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 526 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 527 | def paretovariate(self, alpha): |
| 528 | # Jain, pg. 495 |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 529 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 530 | u = self.random() |
| 531 | return 1.0 / pow(u, 1.0/alpha) |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 532 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 533 | ## -------------------- Weibull -------------------- |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 534 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 535 | def weibullvariate(self, alpha, beta): |
| 536 | # Jain, pg. 499; bug fix courtesy Bill Arms |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 537 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 538 | u = self.random() |
| 539 | return alpha * pow(-_log(u), 1.0/beta) |
Guido van Rossum | 6c395ba | 1999-08-18 13:53:28 +0000 | [diff] [blame] | 540 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 541 | ## -------------------- test program -------------------- |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 542 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 543 | def _test_generator(n, funccall): |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 544 | import time |
| 545 | print n, 'times', funccall |
| 546 | code = compile(funccall, funccall, 'eval') |
| 547 | sum = 0.0 |
| 548 | sqsum = 0.0 |
| 549 | smallest = 1e10 |
| 550 | largest = -1e10 |
| 551 | t0 = time.time() |
| 552 | for i in range(n): |
| 553 | x = eval(code) |
| 554 | sum = sum + x |
| 555 | sqsum = sqsum + x*x |
| 556 | smallest = min(x, smallest) |
| 557 | largest = max(x, largest) |
| 558 | t1 = time.time() |
| 559 | print round(t1-t0, 3), 'sec,', |
| 560 | avg = sum/n |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 561 | stddev = _sqrt(sqsum/n - avg*avg) |
Tim Peters | 0c9886d | 2001-01-15 01:18:21 +0000 | [diff] [blame] | 562 | print 'avg %g, stddev %g, min %g, max %g' % \ |
| 563 | (avg, stddev, smallest, largest) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 564 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 565 | def _test(N=200): |
| 566 | print 'TWOPI =', TWOPI |
| 567 | print 'LOG4 =', LOG4 |
| 568 | print 'NV_MAGICCONST =', NV_MAGICCONST |
| 569 | print 'SG_MAGICCONST =', SG_MAGICCONST |
| 570 | _test_generator(N, 'random()') |
| 571 | _test_generator(N, 'normalvariate(0.0, 1.0)') |
| 572 | _test_generator(N, 'lognormvariate(0.0, 1.0)') |
| 573 | _test_generator(N, 'cunifvariate(0.0, 1.0)') |
| 574 | _test_generator(N, 'expovariate(1.0)') |
| 575 | _test_generator(N, 'vonmisesvariate(0.0, 1.0)') |
| 576 | _test_generator(N, 'gammavariate(0.5, 1.0)') |
| 577 | _test_generator(N, 'gammavariate(0.9, 1.0)') |
| 578 | _test_generator(N, 'gammavariate(1.0, 1.0)') |
| 579 | _test_generator(N, 'gammavariate(2.0, 1.0)') |
| 580 | _test_generator(N, 'gammavariate(20.0, 1.0)') |
| 581 | _test_generator(N, 'gammavariate(200.0, 1.0)') |
| 582 | _test_generator(N, 'gauss(0.0, 1.0)') |
| 583 | _test_generator(N, 'betavariate(3.0, 3.0)') |
| 584 | _test_generator(N, 'paretovariate(1.0)') |
| 585 | _test_generator(N, 'weibullvariate(1.0, 1.0)') |
| 586 | |
Tim Peters | cd80410 | 2001-01-25 20:25:57 +0000 | [diff] [blame] | 587 | # Test jumpahead. |
| 588 | s = getstate() |
| 589 | jumpahead(N) |
| 590 | r1 = random() |
| 591 | # now do it the slow way |
| 592 | setstate(s) |
| 593 | for i in range(N): |
| 594 | random() |
| 595 | r2 = random() |
| 596 | if r1 != r2: |
| 597 | raise ValueError("jumpahead test failed " + `(N, r1, r2)`) |
| 598 | |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 599 | # Initialize from current time. |
| 600 | _inst = Random() |
| 601 | seed = _inst.seed |
| 602 | random = _inst.random |
| 603 | uniform = _inst.uniform |
| 604 | randint = _inst.randint |
| 605 | choice = _inst.choice |
| 606 | randrange = _inst.randrange |
| 607 | shuffle = _inst.shuffle |
| 608 | normalvariate = _inst.normalvariate |
| 609 | lognormvariate = _inst.lognormvariate |
| 610 | cunifvariate = _inst.cunifvariate |
| 611 | expovariate = _inst.expovariate |
| 612 | vonmisesvariate = _inst.vonmisesvariate |
| 613 | gammavariate = _inst.gammavariate |
| 614 | stdgamma = _inst.stdgamma |
| 615 | gauss = _inst.gauss |
| 616 | betavariate = _inst.betavariate |
| 617 | paretovariate = _inst.paretovariate |
| 618 | weibullvariate = _inst.weibullvariate |
| 619 | getstate = _inst.getstate |
| 620 | setstate = _inst.setstate |
Tim Peters | d52269b | 2001-01-25 06:23:18 +0000 | [diff] [blame] | 621 | jumpahead = _inst.jumpahead |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 622 | |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 623 | if __name__ == '__main__': |
Tim Peters | d7b5e88 | 2001-01-25 03:36:26 +0000 | [diff] [blame] | 624 | _test() |