blob: 80af32b0f1f116178206a60a09e80cd616983ee6 [file] [log] [blame]
Guido van Rossume7b146f2000-02-04 15:28:42 +00001"""Random variable generators.
Guido van Rossumff03b1a1994-03-09 12:55:02 +00002
Tim Petersd7b5e882001-01-25 03:36:26 +00003 integers
4 --------
5 uniform within range
6
7 sequences
8 ---------
9 pick random element
Raymond Hettingerf24eb352002-11-12 17:41:57 +000010 pick random sample
Tim Petersd7b5e882001-01-25 03:36:26 +000011 generate random permutation
12
Guido van Rossume7b146f2000-02-04 15:28:42 +000013 distributions on the real line:
14 ------------------------------
Tim Petersd7b5e882001-01-25 03:36:26 +000015 uniform
Guido van Rossume7b146f2000-02-04 15:28:42 +000016 normal (Gaussian)
17 lognormal
18 negative exponential
19 gamma
20 beta
Raymond Hettinger40f62172002-12-29 23:03:38 +000021 pareto
22 Weibull
Guido van Rossumff03b1a1994-03-09 12:55:02 +000023
Guido van Rossume7b146f2000-02-04 15:28:42 +000024 distributions on the circle (angles 0 to 2pi)
25 ---------------------------------------------
26 circular uniform
27 von Mises
28
Raymond Hettinger40f62172002-12-29 23:03:38 +000029General notes on the underlying Mersenne Twister core generator:
Guido van Rossume7b146f2000-02-04 15:28:42 +000030
Raymond Hettinger40f62172002-12-29 23:03:38 +000031* 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 Peterse360d952001-01-26 10:00:39 +000039
Guido van Rossume7b146f2000-02-04 15:28:42 +000040"""
Guido van Rossumd03e1191998-05-29 17:51:31 +000041
Raymond Hettinger2f726e92003-10-05 09:09:15 +000042from warnings import warn as _warn
43from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
Tim Petersd7b5e882001-01-25 03:36:26 +000044from math import log as _log, exp as _exp, pi as _pi, e as _e
45from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +000046from os import urandom as _urandom
47from binascii import hexlify as _hexlify
Guido van Rossumff03b1a1994-03-09 12:55:02 +000048
Raymond Hettingerf24eb352002-11-12 17:41:57 +000049__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000050 "randrange","shuffle","normalvariate","lognormvariate",
Raymond Hettingerf8a52d32003-08-05 12:23:19 +000051 "expovariate","vonmisesvariate","gammavariate",
52 "gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger356a4592004-08-30 06:14:31 +000053 "getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
Raymond Hettinger23f12412004-09-13 22:23:21 +000054 "SystemRandom"]
Tim Petersd7b5e882001-01-25 03:36:26 +000055
56NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000057TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000058LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000059SG_MAGICCONST = 1.0 + _log(4.5)
Raymond Hettinger2f726e92003-10-05 09:09:15 +000060BPF = 53 # Number of bits in a float
Tim Peters7c2a85b2004-08-31 02:19:55 +000061RECIP_BPF = 2**-BPF
Tim Petersd7b5e882001-01-25 03:36:26 +000062
Raymond Hettinger356a4592004-08-30 06:14:31 +000063
Tim Petersd7b5e882001-01-25 03:36:26 +000064# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000065# Adrian Baddeley. Adapted by Raymond Hettinger for use with
Raymond Hettinger3fa19d72004-08-31 01:05:15 +000066# the Mersenne Twister and os.urandom() core generators.
Tim Petersd7b5e882001-01-25 03:36:26 +000067
Raymond Hettinger145a4a02003-01-07 10:25:55 +000068import _random
Raymond Hettinger40f62172002-12-29 23:03:38 +000069
Raymond Hettinger145a4a02003-01-07 10:25:55 +000070class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000071 """Random number generator base class used by bound module functions.
72
73 Used to instantiate instances of Random to get generators that don't
74 share state. Especially useful for multi-threaded programs, creating
75 a different instance of Random for each thread, and using the jumpahead()
76 method to ensure that the generated sequences seen by each thread don't
77 overlap.
78
79 Class Random can also be subclassed if you want to use a different basic
80 generator of your own devising: in that case, override the following
81 methods: random(), seed(), getstate(), setstate() and jumpahead().
Raymond Hettinger2f726e92003-10-05 09:09:15 +000082 Optionally, implement a getrandombits() method so that randrange()
83 can cover arbitrarily large ranges.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +000084
Raymond Hettingerc32f0332002-05-23 19:44:49 +000085 """
Tim Petersd7b5e882001-01-25 03:36:26 +000086
Raymond Hettinger40f62172002-12-29 23:03:38 +000087 VERSION = 2 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000088
89 def __init__(self, x=None):
90 """Initialize an instance.
91
92 Optional argument x controls seeding, as for Random.seed().
93 """
94
95 self.seed(x)
Raymond Hettinger40f62172002-12-29 23:03:38 +000096 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +000097
Tim Peters0de88fc2001-02-01 04:59:18 +000098 def seed(self, a=None):
99 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +0000100
Raymond Hettinger23f12412004-09-13 22:23:21 +0000101 None or no argument seeds from current time or from an operating
102 system specific randomness source if available.
Tim Peters0de88fc2001-02-01 04:59:18 +0000103
Tim Petersbcd725f2001-02-01 10:06:53 +0000104 If a is not None or an int or long, hash(a) is used instead.
Tim Petersd7b5e882001-01-25 03:36:26 +0000105 """
106
Raymond Hettinger3081d592003-08-09 18:30:57 +0000107 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000108 try:
109 a = long(_hexlify(_urandom(16)), 16)
110 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000111 import time
112 a = long(time.time() * 256) # use fractional seconds
Raymond Hettinger356a4592004-08-30 06:14:31 +0000113
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000114 super(Random, self).seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +0000115 self.gauss_next = None
116
Tim Peterscd804102001-01-25 20:25:57 +0000117 def getstate(self):
118 """Return internal state; can be passed to setstate() later."""
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000119 return self.VERSION, super(Random, self).getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000120
121 def setstate(self, state):
122 """Restore internal state from object returned by getstate()."""
123 version = state[0]
Raymond Hettinger40f62172002-12-29 23:03:38 +0000124 if version == 2:
125 version, internalstate, self.gauss_next = state
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000126 super(Random, self).setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000127 else:
128 raise ValueError("state with version %s passed to "
129 "Random.setstate() of version %s" %
130 (version, self.VERSION))
131
Tim Peterscd804102001-01-25 20:25:57 +0000132## ---- Methods below this point do not need to be overridden when
133## ---- subclassing for the purpose of using a different core generator.
134
135## -------------------- pickle support -------------------
136
137 def __getstate__(self): # for pickle
138 return self.getstate()
139
140 def __setstate__(self, state): # for pickle
141 self.setstate(state)
142
Raymond Hettinger5f078ff2003-06-24 20:29:04 +0000143 def __reduce__(self):
144 return self.__class__, (), self.getstate()
145
Tim Peterscd804102001-01-25 20:25:57 +0000146## -------------------- integer methods -------------------
147
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000148 def randrange(self, start, stop=None, step=1, int=int, default=None,
149 maxwidth=1L<<BPF):
Tim Petersd7b5e882001-01-25 03:36:26 +0000150 """Choose a random item from range(start, stop[, step]).
151
152 This fixes the problem with randint() which includes the
153 endpoint; in Python this is usually not what you want.
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000154 Do not supply the 'int', 'default', and 'maxwidth' arguments.
Tim Petersd7b5e882001-01-25 03:36:26 +0000155 """
156
157 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000158 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000159 istart = int(start)
160 if istart != start:
161 raise ValueError, "non-integer arg 1 for randrange()"
162 if stop is default:
163 if istart > 0:
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000164 if istart >= maxwidth:
165 return self._randbelow(istart)
Tim Petersd7b5e882001-01-25 03:36:26 +0000166 return int(self.random() * istart)
167 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000168
169 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000170 istop = int(stop)
171 if istop != stop:
172 raise ValueError, "non-integer stop for randrange()"
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000173 width = istop - istart
174 if step == 1 and width > 0:
Tim Peters76ca1d42003-06-19 03:46:46 +0000175 # Note that
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000176 # int(istart + self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000177 # instead would be incorrect. For example, consider istart
178 # = -2 and istop = 0. Then the guts would be in
179 # -2.0 to 0.0 exclusive on both ends (ignoring that random()
180 # might return 0.0), and because int() truncates toward 0, the
181 # final result would be -1 or 0 (instead of -2 or -1).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000182 # istart + int(self.random()*width)
Tim Peters76ca1d42003-06-19 03:46:46 +0000183 # would also be incorrect, for a subtler reason: the RHS
184 # can return a long, and then randrange() would also return
185 # a long, but we're supposed to return an int (for backward
186 # compatibility).
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000187
188 if width >= maxwidth:
Tim Peters58eb11c2004-01-18 20:29:55 +0000189 return int(istart + self._randbelow(width))
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000190 return int(istart + int(self.random()*width))
Tim Petersd7b5e882001-01-25 03:36:26 +0000191 if step == 1:
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000192 raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)
Tim Peters9146f272002-08-16 03:41:39 +0000193
194 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000195 istep = int(step)
196 if istep != step:
197 raise ValueError, "non-integer step for randrange()"
198 if istep > 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000199 n = (width + istep - 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000200 elif istep < 0:
Raymond Hettingerffdb8bb2004-09-27 15:29:05 +0000201 n = (width + istep + 1) // istep
Tim Petersd7b5e882001-01-25 03:36:26 +0000202 else:
203 raise ValueError, "zero step for randrange()"
204
205 if n <= 0:
206 raise ValueError, "empty range for randrange()"
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000207
208 if n >= maxwidth:
209 return istart + self._randbelow(n)
Tim Petersd7b5e882001-01-25 03:36:26 +0000210 return istart + istep*int(self.random() * n)
211
212 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000213 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000214 """
215
216 return self.randrange(a, b+1)
217
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000218 def _randbelow(self, n, _log=_log, int=int, _maxwidth=1L<<BPF,
219 _Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
220 """Return a random int in the range [0,n)
221
222 Handles the case where n has more bits than returned
223 by a single call to the underlying generator.
224 """
225
226 try:
227 getrandbits = self.getrandbits
228 except AttributeError:
229 pass
230 else:
231 # Only call self.getrandbits if the original random() builtin method
232 # has not been overridden or if a new getrandbits() was supplied.
233 # This assures that the two methods correspond.
234 if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
235 k = int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2)
236 r = getrandbits(k)
237 while r >= n:
238 r = getrandbits(k)
239 return r
240 if n >= _maxwidth:
241 _warn("Underlying random() generator does not supply \n"
242 "enough bits to choose from a population range this large")
243 return int(self.random() * n)
244
Tim Peterscd804102001-01-25 20:25:57 +0000245## -------------------- sequence methods -------------------
246
Tim Petersd7b5e882001-01-25 03:36:26 +0000247 def choice(self, seq):
248 """Choose a random element from a non-empty sequence."""
Raymond Hettinger5dae5052004-06-07 02:07:15 +0000249 return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty
Tim Petersd7b5e882001-01-25 03:36:26 +0000250
251 def shuffle(self, x, random=None, int=int):
252 """x, random=random.random -> shuffle list x in place; return None.
253
254 Optional arg random is a 0-argument function returning a random
255 float in [0.0, 1.0); by default, the standard random.random.
256
257 Note that for even rather small len(x), the total number of
258 permutations of x is larger than the period of most random number
259 generators; this implies that "most" permutations of a long
260 sequence can never be generated.
261 """
262
263 if random is None:
264 random = self.random
Raymond Hettinger85c20a42003-11-06 14:06:48 +0000265 for i in reversed(xrange(1, len(x))):
Tim Peterscd804102001-01-25 20:25:57 +0000266 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000267 j = int(random() * (i+1))
268 x[i], x[j] = x[j], x[i]
269
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000270 def sample(self, population, k):
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000271 """Chooses k unique random elements from a population sequence.
272
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000273 Returns a new list containing elements from the population while
274 leaving the original population unchanged. The resulting list is
275 in selection order so that all sub-slices will also be valid random
276 samples. This allows raffle winners (the sample) to be partitioned
277 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000278
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000279 Members of the population need not be hashable or unique. If the
280 population contains repeats, then each occurrence is a possible
281 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000282
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000283 To choose a sample in a range of integers, use xrange as an argument.
284 This is especially fast and space efficient for sampling from a
285 large population: sample(xrange(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000286 """
287
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000288 # Sampling without replacement entails tracking either potential
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000289 # selections (the pool) in a list or previous selections in a
290 # dictionary.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000291
Jeremy Hylton2b55d352004-02-23 17:27:57 +0000292 # When the number of selections is small compared to the
293 # population, then tracking selections is efficient, requiring
294 # only a small dictionary and an occasional reselection. For
295 # a larger number of selections, the pool tracking method is
296 # preferred since the list takes less space than the
297 # dictionary and it doesn't suffer from frequent reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000298
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000299 n = len(population)
300 if not 0 <= k <= n:
301 raise ValueError, "sample larger than population"
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000302 random = self.random
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000303 _int = int
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000304 result = [None] * k
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000305 if n < 6 * k: # if n len list takes less space than a k len dict
Raymond Hettinger311f4192002-11-18 09:01:24 +0000306 pool = list(population)
307 for i in xrange(k): # invariant: non-selected at [0,n-i)
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000308 j = _int(random() * (n-i))
Raymond Hettinger311f4192002-11-18 09:01:24 +0000309 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000310 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000311 else:
Raymond Hettinger66d09f12003-09-06 04:25:54 +0000312 try:
313 n > 0 and (population[0], population[n//2], population[n-1])
314 except (TypeError, KeyError): # handle sets and dictionaries
315 population = tuple(population)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000316 selected = {}
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000317 for i in xrange(k):
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000318 j = _int(random() * n)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000319 while j in selected:
Raymond Hettingerfdbe5222003-06-13 07:01:51 +0000320 j = _int(random() * n)
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000321 result[i] = selected[j] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000322 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000323
Tim Peterscd804102001-01-25 20:25:57 +0000324## -------------------- real-valued distributions -------------------
325
326## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000327
328 def uniform(self, a, b):
329 """Get a random number in the range [a, b)."""
330 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000331
Tim Peterscd804102001-01-25 20:25:57 +0000332## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000333
Tim Petersd7b5e882001-01-25 03:36:26 +0000334 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000335 """Normal distribution.
336
337 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000338
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000339 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000340 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000341
Tim Petersd7b5e882001-01-25 03:36:26 +0000342 # Uses Kinderman and Monahan method. Reference: Kinderman,
343 # A.J. and Monahan, J.F., "Computer generation of random
344 # variables using the ratio of uniform deviates", ACM Trans
345 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000346
Tim Petersd7b5e882001-01-25 03:36:26 +0000347 random = self.random
Raymond Hettinger42406e62005-04-30 09:02:51 +0000348 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000349 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000350 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000351 z = NV_MAGICCONST*(u1-0.5)/u2
352 zz = z*z/4.0
353 if zz <= -_log(u2):
354 break
355 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000356
Tim Peterscd804102001-01-25 20:25:57 +0000357## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000358
359 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000360 """Log normal distribution.
361
362 If you take the natural logarithm of this distribution, you'll get a
363 normal distribution with mean mu and standard deviation sigma.
364 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000365
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000366 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000367 return _exp(self.normalvariate(mu, sigma))
368
Tim Peterscd804102001-01-25 20:25:57 +0000369## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000370
371 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000372 """Exponential distribution.
373
374 lambd is 1.0 divided by the desired mean. (The parameter would be
375 called "lambda", but that is a reserved word in Python.) Returned
376 values range from 0 to positive infinity.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000377
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000378 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000379 # lambd: rate lambd = 1/mean
380 # ('lambda' is a Python reserved word)
381
382 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000383 u = random()
384 while u <= 1e-7:
385 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000386 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000387
Tim Peterscd804102001-01-25 20:25:57 +0000388## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000389
Tim Petersd7b5e882001-01-25 03:36:26 +0000390 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000391 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000392
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000393 mu is the mean angle, expressed in radians between 0 and 2*pi, and
394 kappa is the concentration parameter, which must be greater than or
395 equal to zero. If kappa is equal to zero, this distribution reduces
396 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000397
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000398 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000399 # mu: mean angle (in radians between 0 and 2*pi)
400 # kappa: concentration parameter kappa (>= 0)
401 # if kappa = 0 generate uniform random angle
402
403 # Based upon an algorithm published in: Fisher, N.I.,
404 # "Statistical Analysis of Circular Data", Cambridge
405 # University Press, 1993.
406
407 # Thanks to Magnus Kessler for a correction to the
408 # implementation of step 4.
409
410 random = self.random
411 if kappa <= 1e-6:
412 return TWOPI * random()
413
414 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
415 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
416 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000417
Raymond Hettinger42406e62005-04-30 09:02:51 +0000418 while 1:
Tim Peters0c9886d2001-01-15 01:18:21 +0000419 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000420
421 z = _cos(_pi * u1)
422 f = (1.0 + r * z)/(r + z)
423 c = kappa * (r - f)
424
425 u2 = random()
426
Raymond Hettinger42406e62005-04-30 09:02:51 +0000427 if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
Tim Peters0c9886d2001-01-15 01:18:21 +0000428 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000429
430 u3 = random()
431 if u3 > 0.5:
432 theta = (mu % TWOPI) + _acos(f)
433 else:
434 theta = (mu % TWOPI) - _acos(f)
435
436 return theta
437
Tim Peterscd804102001-01-25 20:25:57 +0000438## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000439
440 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000441 """Gamma distribution. Not the gamma function!
442
443 Conditions on the parameters are alpha > 0 and beta > 0.
444
445 """
Tim Peters8ac14952002-05-23 15:15:30 +0000446
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000447 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000448
Guido van Rossum570764d2002-05-14 14:08:12 +0000449 # Warning: a few older sources define the gamma distribution in terms
450 # of alpha > -1.0
451 if alpha <= 0.0 or beta <= 0.0:
452 raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters8ac14952002-05-23 15:15:30 +0000453
Tim Petersd7b5e882001-01-25 03:36:26 +0000454 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000455 if alpha > 1.0:
456
457 # Uses R.C.H. Cheng, "The generation of Gamma
458 # variables with non-integral shape parameters",
459 # Applied Statistics, (1977), 26, No. 1, p71-74
460
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000461 ainv = _sqrt(2.0 * alpha - 1.0)
462 bbb = alpha - LOG4
463 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000464
Raymond Hettinger42406e62005-04-30 09:02:51 +0000465 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000466 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000467 if not 1e-7 < u1 < .9999999:
468 continue
469 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000470 v = _log(u1/(1.0-u1))/ainv
471 x = alpha*_exp(v)
472 z = u1*u1*u2
473 r = bbb+ccc*v-x
474 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000475 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000476
477 elif alpha == 1.0:
478 # expovariate(1)
479 u = random()
480 while u <= 1e-7:
481 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000482 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000483
484 else: # alpha is between 0 and 1 (exclusive)
485
486 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
487
Raymond Hettinger42406e62005-04-30 09:02:51 +0000488 while 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000489 u = random()
490 b = (_e + alpha)/_e
491 p = b*u
492 if p <= 1.0:
Raymond Hettinger42406e62005-04-30 09:02:51 +0000493 x = p ** (1.0/alpha)
Tim Petersd7b5e882001-01-25 03:36:26 +0000494 else:
Tim Petersd7b5e882001-01-25 03:36:26 +0000495 x = -_log((b-p)/alpha)
496 u1 = random()
Raymond Hettinger42406e62005-04-30 09:02:51 +0000497 if p > 1.0:
498 if u1 <= x ** (alpha - 1.0):
499 break
500 elif u1 <= _exp(-x):
Tim Petersd7b5e882001-01-25 03:36:26 +0000501 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000502 return x * beta
503
Tim Peterscd804102001-01-25 20:25:57 +0000504## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000505
Tim Petersd7b5e882001-01-25 03:36:26 +0000506 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000507 """Gaussian distribution.
508
509 mu is the mean, and sigma is the standard deviation. This is
510 slightly faster than the normalvariate() function.
511
512 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000513
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000514 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000515
Tim Petersd7b5e882001-01-25 03:36:26 +0000516 # When x and y are two variables from [0, 1), uniformly
517 # distributed, then
518 #
519 # cos(2*pi*x)*sqrt(-2*log(1-y))
520 # sin(2*pi*x)*sqrt(-2*log(1-y))
521 #
522 # are two *independent* variables with normal distribution
523 # (mu = 0, sigma = 1).
524 # (Lambert Meertens)
525 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000526
Tim Petersd7b5e882001-01-25 03:36:26 +0000527 # Multithreading note: When two threads call this function
528 # simultaneously, it is possible that they will receive the
529 # same return value. The window is very small though. To
530 # avoid this, you have to use a lock around all calls. (I
531 # didn't want to slow this down in the serial case by using a
532 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000533
Tim Petersd7b5e882001-01-25 03:36:26 +0000534 random = self.random
535 z = self.gauss_next
536 self.gauss_next = None
537 if z is None:
538 x2pi = random() * TWOPI
539 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
540 z = _cos(x2pi) * g2rad
541 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000542
Tim Petersd7b5e882001-01-25 03:36:26 +0000543 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000544
Tim Peterscd804102001-01-25 20:25:57 +0000545## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000546## See
547## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
548## for Ivan Frohne's insightful analysis of why the original implementation:
549##
550## def betavariate(self, alpha, beta):
551## # Discrete Event Simulation in C, pp 87-88.
552##
553## y = self.expovariate(alpha)
554## z = self.expovariate(1.0/beta)
555## return z/(y+z)
556##
557## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000558
Tim Petersd7b5e882001-01-25 03:36:26 +0000559 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000560 """Beta distribution.
561
562 Conditions on the parameters are alpha > -1 and beta} > -1.
563 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000564
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000565 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000566
Tim Peters85e2e472001-01-26 06:49:56 +0000567 # This version due to Janne Sinkkonen, and matches all the std
568 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
569 y = self.gammavariate(alpha, 1.)
570 if y == 0:
571 return 0.0
572 else:
573 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000574
Tim Peterscd804102001-01-25 20:25:57 +0000575## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000576
Tim Petersd7b5e882001-01-25 03:36:26 +0000577 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000578 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000579 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000580
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000581 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000582 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000583
Tim Peterscd804102001-01-25 20:25:57 +0000584## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000585
Tim Petersd7b5e882001-01-25 03:36:26 +0000586 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000587 """Weibull distribution.
588
589 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000590
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000591 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000592 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000593
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000594 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000595 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000596
Raymond Hettinger40f62172002-12-29 23:03:38 +0000597## -------------------- Wichmann-Hill -------------------
598
599class WichmannHill(Random):
600
601 VERSION = 1 # used by getstate/setstate
602
603 def seed(self, a=None):
604 """Initialize internal state from hashable object.
605
Raymond Hettinger23f12412004-09-13 22:23:21 +0000606 None or no argument seeds from current time or from an operating
607 system specific randomness source if available.
Raymond Hettinger40f62172002-12-29 23:03:38 +0000608
609 If a is not None or an int or long, hash(a) is used instead.
610
611 If a is an int or long, a is used directly. Distinct values between
612 0 and 27814431486575L inclusive are guaranteed to yield distinct
613 internal states (this guarantee is specific to the default
614 Wichmann-Hill generator).
615 """
616
617 if a is None:
Raymond Hettingerc1c43ca2004-09-05 00:00:42 +0000618 try:
619 a = long(_hexlify(_urandom(16)), 16)
620 except NotImplementedError:
Raymond Hettinger356a4592004-08-30 06:14:31 +0000621 import time
622 a = long(time.time() * 256) # use fractional seconds
Raymond Hettinger40f62172002-12-29 23:03:38 +0000623
624 if not isinstance(a, (int, long)):
625 a = hash(a)
626
627 a, x = divmod(a, 30268)
628 a, y = divmod(a, 30306)
629 a, z = divmod(a, 30322)
630 self._seed = int(x)+1, int(y)+1, int(z)+1
631
632 self.gauss_next = None
633
634 def random(self):
635 """Get the next random number in the range [0.0, 1.0)."""
636
637 # Wichman-Hill random number generator.
638 #
639 # Wichmann, B. A. & Hill, I. D. (1982)
640 # Algorithm AS 183:
641 # An efficient and portable pseudo-random number generator
642 # Applied Statistics 31 (1982) 188-190
643 #
644 # see also:
645 # Correction to Algorithm AS 183
646 # Applied Statistics 33 (1984) 123
647 #
648 # McLeod, A. I. (1985)
649 # A remark on Algorithm AS 183
650 # Applied Statistics 34 (1985),198-200
651
652 # This part is thread-unsafe:
653 # BEGIN CRITICAL SECTION
654 x, y, z = self._seed
655 x = (171 * x) % 30269
656 y = (172 * y) % 30307
657 z = (170 * z) % 30323
658 self._seed = x, y, z
659 # END CRITICAL SECTION
660
661 # Note: on a platform using IEEE-754 double arithmetic, this can
662 # never return 0.0 (asserted by Tim; proof too long for a comment).
663 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
664
665 def getstate(self):
666 """Return internal state; can be passed to setstate() later."""
667 return self.VERSION, self._seed, self.gauss_next
668
669 def setstate(self, state):
670 """Restore internal state from object returned by getstate()."""
671 version = state[0]
672 if version == 1:
673 version, self._seed, self.gauss_next = state
674 else:
675 raise ValueError("state with version %s passed to "
676 "Random.setstate() of version %s" %
677 (version, self.VERSION))
678
679 def jumpahead(self, n):
680 """Act as if n calls to random() were made, but quickly.
681
682 n is an int, greater than or equal to 0.
683
684 Example use: If you have 2 threads and know that each will
685 consume no more than a million random numbers, create two Random
686 objects r1 and r2, then do
687 r2.setstate(r1.getstate())
688 r2.jumpahead(1000000)
689 Then r1 and r2 will use guaranteed-disjoint segments of the full
690 period.
691 """
692
693 if not n >= 0:
694 raise ValueError("n must be >= 0")
695 x, y, z = self._seed
696 x = int(x * pow(171, n, 30269)) % 30269
697 y = int(y * pow(172, n, 30307)) % 30307
698 z = int(z * pow(170, n, 30323)) % 30323
699 self._seed = x, y, z
700
701 def __whseed(self, x=0, y=0, z=0):
702 """Set the Wichmann-Hill seed from (x, y, z).
703
704 These must be integers in the range [0, 256).
705 """
706
707 if not type(x) == type(y) == type(z) == int:
708 raise TypeError('seeds must be integers')
709 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
710 raise ValueError('seeds must be in range(0, 256)')
711 if 0 == x == y == z:
712 # Initialize from current time
713 import time
714 t = long(time.time() * 256)
715 t = int((t&0xffffff) ^ (t>>24))
716 t, x = divmod(t, 256)
717 t, y = divmod(t, 256)
718 t, z = divmod(t, 256)
719 # Zero is a poor seed, so substitute 1
720 self._seed = (x or 1, y or 1, z or 1)
721
722 self.gauss_next = None
723
724 def whseed(self, a=None):
725 """Seed from hashable object's hash code.
726
727 None or no argument seeds from current time. It is not guaranteed
728 that objects with distinct hash codes lead to distinct internal
729 states.
730
731 This is obsolete, provided for compatibility with the seed routine
732 used prior to Python 2.1. Use the .seed() method instead.
733 """
734
735 if a is None:
736 self.__whseed()
737 return
738 a = hash(a)
739 a, x = divmod(a, 256)
740 a, y = divmod(a, 256)
741 a, z = divmod(a, 256)
742 x = (x + a) % 256 or 1
743 y = (y + a) % 256 or 1
744 z = (z + a) % 256 or 1
745 self.__whseed(x, y, z)
746
Raymond Hettinger23f12412004-09-13 22:23:21 +0000747## --------------- Operating System Random Source ------------------
Raymond Hettinger356a4592004-08-30 06:14:31 +0000748
Raymond Hettinger23f12412004-09-13 22:23:21 +0000749class SystemRandom(Random):
750 """Alternate random number generator using sources provided
751 by the operating system (such as /dev/urandom on Unix or
752 CryptGenRandom on Windows).
Raymond Hettinger356a4592004-08-30 06:14:31 +0000753
754 Not available on all systems (see os.urandom() for details).
755 """
756
757 def random(self):
758 """Get the next random number in the range [0.0, 1.0)."""
Tim Peters7c2a85b2004-08-31 02:19:55 +0000759 return (long(_hexlify(_urandom(7)), 16) >> 3) * RECIP_BPF
Raymond Hettinger356a4592004-08-30 06:14:31 +0000760
761 def getrandbits(self, k):
762 """getrandbits(k) -> x. Generates a long int with k random bits."""
Raymond Hettinger356a4592004-08-30 06:14:31 +0000763 if k <= 0:
764 raise ValueError('number of bits must be greater than zero')
765 if k != int(k):
766 raise TypeError('number of bits should be an integer')
767 bytes = (k + 7) // 8 # bits / 8 and rounded up
768 x = long(_hexlify(_urandom(bytes)), 16)
769 return x >> (bytes * 8 - k) # trim excess bits
770
771 def _stub(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000772 "Stub method. Not used for a system random number generator."
Raymond Hettinger356a4592004-08-30 06:14:31 +0000773 return None
774 seed = jumpahead = _stub
775
776 def _notimplemented(self, *args, **kwds):
Raymond Hettinger23f12412004-09-13 22:23:21 +0000777 "Method should not be called for a system random number generator."
778 raise NotImplementedError('System entropy source does not have state.')
Raymond Hettinger356a4592004-08-30 06:14:31 +0000779 getstate = setstate = _notimplemented
780
Tim Peterscd804102001-01-25 20:25:57 +0000781## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000782
Raymond Hettinger62297132003-08-30 01:24:19 +0000783def _test_generator(n, func, args):
Tim Peters0c9886d2001-01-15 01:18:21 +0000784 import time
Raymond Hettinger62297132003-08-30 01:24:19 +0000785 print n, 'times', func.__name__
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000786 total = 0.0
Tim Peters0c9886d2001-01-15 01:18:21 +0000787 sqsum = 0.0
788 smallest = 1e10
789 largest = -1e10
790 t0 = time.time()
791 for i in range(n):
Raymond Hettinger62297132003-08-30 01:24:19 +0000792 x = func(*args)
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000793 total += x
Tim Peters0c9886d2001-01-15 01:18:21 +0000794 sqsum = sqsum + x*x
795 smallest = min(x, smallest)
796 largest = max(x, largest)
797 t1 = time.time()
798 print round(t1-t0, 3), 'sec,',
Raymond Hettingerb98154e2003-05-24 17:26:02 +0000799 avg = total/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000800 stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters0c9886d2001-01-15 01:18:21 +0000801 print 'avg %g, stddev %g, min %g, max %g' % \
802 (avg, stddev, smallest, largest)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000803
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000804
805def _test(N=2000):
Raymond Hettinger62297132003-08-30 01:24:19 +0000806 _test_generator(N, random, ())
807 _test_generator(N, normalvariate, (0.0, 1.0))
808 _test_generator(N, lognormvariate, (0.0, 1.0))
809 _test_generator(N, vonmisesvariate, (0.0, 1.0))
810 _test_generator(N, gammavariate, (0.01, 1.0))
811 _test_generator(N, gammavariate, (0.1, 1.0))
812 _test_generator(N, gammavariate, (0.1, 2.0))
813 _test_generator(N, gammavariate, (0.5, 1.0))
814 _test_generator(N, gammavariate, (0.9, 1.0))
815 _test_generator(N, gammavariate, (1.0, 1.0))
816 _test_generator(N, gammavariate, (2.0, 1.0))
817 _test_generator(N, gammavariate, (20.0, 1.0))
818 _test_generator(N, gammavariate, (200.0, 1.0))
819 _test_generator(N, gauss, (0.0, 1.0))
820 _test_generator(N, betavariate, (3.0, 3.0))
Tim Peterscd804102001-01-25 20:25:57 +0000821
Tim Peters715c4c42001-01-26 22:56:56 +0000822# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000823# as module-level functions. The functions share state across all uses
824#(both in the user's code and in the Python libraries), but that's fine
825# for most programs and is easier for the casual user than making them
826# instantiate their own Random() instance.
827
Tim Petersd7b5e882001-01-25 03:36:26 +0000828_inst = Random()
829seed = _inst.seed
830random = _inst.random
831uniform = _inst.uniform
832randint = _inst.randint
833choice = _inst.choice
834randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000835sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000836shuffle = _inst.shuffle
837normalvariate = _inst.normalvariate
838lognormvariate = _inst.lognormvariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000839expovariate = _inst.expovariate
840vonmisesvariate = _inst.vonmisesvariate
841gammavariate = _inst.gammavariate
Tim Petersd7b5e882001-01-25 03:36:26 +0000842gauss = _inst.gauss
843betavariate = _inst.betavariate
844paretovariate = _inst.paretovariate
845weibullvariate = _inst.weibullvariate
846getstate = _inst.getstate
847setstate = _inst.setstate
Tim Petersd52269b2001-01-25 06:23:18 +0000848jumpahead = _inst.jumpahead
Raymond Hettinger2f726e92003-10-05 09:09:15 +0000849getrandbits = _inst.getrandbits
Tim Petersd7b5e882001-01-25 03:36:26 +0000850
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000851if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000852 _test()