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