blob: 3db1f1a3facd8993e4abef867237909dc1a83c2c [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
Tim Petersd7b5e882001-01-25 03:36:26 +000042from math import log as _log, exp as _exp, pi as _pi, e as _e
43from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
Tim Peters9146f272002-08-16 03:41:39 +000044from math import floor as _floor
Guido van Rossumff03b1a1994-03-09 12:55:02 +000045
Raymond Hettingerf24eb352002-11-12 17:41:57 +000046__all__ = ["Random","seed","random","uniform","randint","choice","sample",
Skip Montanaro0de65802001-02-15 22:15:14 +000047 "randrange","shuffle","normalvariate","lognormvariate",
48 "cunifvariate","expovariate","vonmisesvariate","gammavariate",
49 "stdgamma","gauss","betavariate","paretovariate","weibullvariate",
Raymond Hettinger40f62172002-12-29 23:03:38 +000050 "getstate","setstate","jumpahead"]
Tim Petersd7b5e882001-01-25 03:36:26 +000051
52NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000053TWOPI = 2.0*_pi
Tim Petersd7b5e882001-01-25 03:36:26 +000054LOG4 = _log(4.0)
Tim Petersd7b5e882001-01-25 03:36:26 +000055SG_MAGICCONST = 1.0 + _log(4.5)
Tim Petersd7b5e882001-01-25 03:36:26 +000056
57# Translated by Guido van Rossum from C source provided by
Raymond Hettinger40f62172002-12-29 23:03:38 +000058# Adrian Baddeley. Adapted by Raymond Hettinger for use with
59# the Mersenne Twister core generator.
Tim Petersd7b5e882001-01-25 03:36:26 +000060
Raymond Hettinger145a4a02003-01-07 10:25:55 +000061import _random
Raymond Hettinger40f62172002-12-29 23:03:38 +000062
Raymond Hettinger145a4a02003-01-07 10:25:55 +000063class Random(_random.Random):
Raymond Hettingerc32f0332002-05-23 19:44:49 +000064 """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 Hettingeref4d4bd2002-05-23 23:58:17 +000075
Raymond Hettingerc32f0332002-05-23 19:44:49 +000076 """
Tim Petersd7b5e882001-01-25 03:36:26 +000077
Raymond Hettinger40f62172002-12-29 23:03:38 +000078 VERSION = 2 # used by getstate/setstate
Tim Petersd7b5e882001-01-25 03:36:26 +000079
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 Hettinger40f62172002-12-29 23:03:38 +000087 self.gauss_next = None
Tim Petersd7b5e882001-01-25 03:36:26 +000088
Tim Peters0de88fc2001-02-01 04:59:18 +000089 def seed(self, a=None):
90 """Initialize internal state from hashable object.
Tim Petersd7b5e882001-01-25 03:36:26 +000091
Tim Peters0de88fc2001-02-01 04:59:18 +000092 None or no argument seeds from current time.
93
Tim Petersbcd725f2001-02-01 10:06:53 +000094 If a is not None or an int or long, hash(a) is used instead.
Tim Petersd7b5e882001-01-25 03:36:26 +000095 """
96
Raymond Hettinger145a4a02003-01-07 10:25:55 +000097 super(Random, self).seed(a)
Tim Peters46c04e12002-05-05 20:40:00 +000098 self.gauss_next = None
99
Tim Peterscd804102001-01-25 20:25:57 +0000100 def getstate(self):
101 """Return internal state; can be passed to setstate() later."""
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000102 return self.VERSION, super(Random, self).getstate(), self.gauss_next
Tim Peterscd804102001-01-25 20:25:57 +0000103
104 def setstate(self, state):
105 """Restore internal state from object returned by getstate()."""
106 version = state[0]
Raymond Hettinger40f62172002-12-29 23:03:38 +0000107 if version == 2:
108 version, internalstate, self.gauss_next = state
Raymond Hettinger145a4a02003-01-07 10:25:55 +0000109 super(Random, self).setstate(internalstate)
Tim Peterscd804102001-01-25 20:25:57 +0000110 else:
111 raise ValueError("state with version %s passed to "
112 "Random.setstate() of version %s" %
113 (version, self.VERSION))
114
Tim Peterscd804102001-01-25 20:25:57 +0000115## ---- Methods below this point do not need to be overridden when
116## ---- subclassing for the purpose of using a different core generator.
117
118## -------------------- pickle support -------------------
119
120 def __getstate__(self): # for pickle
121 return self.getstate()
122
123 def __setstate__(self, state): # for pickle
124 self.setstate(state)
125
126## -------------------- integer methods -------------------
127
Tim Petersd7b5e882001-01-25 03:36:26 +0000128 def randrange(self, start, stop=None, step=1, int=int, default=None):
129 """Choose a random item from range(start, stop[, step]).
130
131 This fixes the problem with randint() which includes the
132 endpoint; in Python this is usually not what you want.
133 Do not supply the 'int' and 'default' arguments.
134 """
135
136 # This code is a bit messy to make it fast for the
Tim Peters9146f272002-08-16 03:41:39 +0000137 # common case while still doing adequate error checking.
Tim Petersd7b5e882001-01-25 03:36:26 +0000138 istart = int(start)
139 if istart != start:
140 raise ValueError, "non-integer arg 1 for randrange()"
141 if stop is default:
142 if istart > 0:
143 return int(self.random() * istart)
144 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000145
146 # stop argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000147 istop = int(stop)
148 if istop != stop:
149 raise ValueError, "non-integer stop for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000150 if step == 1 and istart < istop:
151 try:
152 return istart + int(self.random()*(istop - istart))
153 except OverflowError:
154 # This can happen if istop-istart > sys.maxint + 1, and
155 # multiplying by random() doesn't reduce it to something
156 # <= sys.maxint. We know that the overall result fits
157 # in an int, and can still do it correctly via math.floor().
158 # But that adds another function call, so for speed we
159 # avoided that whenever possible.
160 return int(istart + _floor(self.random()*(istop - istart)))
Tim Petersd7b5e882001-01-25 03:36:26 +0000161 if step == 1:
Tim Petersd7b5e882001-01-25 03:36:26 +0000162 raise ValueError, "empty range for randrange()"
Tim Peters9146f272002-08-16 03:41:39 +0000163
164 # Non-unit step argument supplied.
Tim Petersd7b5e882001-01-25 03:36:26 +0000165 istep = int(step)
166 if istep != step:
167 raise ValueError, "non-integer step for randrange()"
168 if istep > 0:
169 n = (istop - istart + istep - 1) / istep
170 elif istep < 0:
171 n = (istop - istart + istep + 1) / istep
172 else:
173 raise ValueError, "zero step for randrange()"
174
175 if n <= 0:
176 raise ValueError, "empty range for randrange()"
177 return istart + istep*int(self.random() * n)
178
179 def randint(self, a, b):
Tim Peterscd804102001-01-25 20:25:57 +0000180 """Return random integer in range [a, b], including both end points.
Tim Petersd7b5e882001-01-25 03:36:26 +0000181 """
182
183 return self.randrange(a, b+1)
184
Tim Peterscd804102001-01-25 20:25:57 +0000185## -------------------- sequence methods -------------------
186
Tim Petersd7b5e882001-01-25 03:36:26 +0000187 def choice(self, seq):
188 """Choose a random element from a non-empty sequence."""
189 return seq[int(self.random() * len(seq))]
190
191 def shuffle(self, x, random=None, int=int):
192 """x, random=random.random -> shuffle list x in place; return None.
193
194 Optional arg random is a 0-argument function returning a random
195 float in [0.0, 1.0); by default, the standard random.random.
196
197 Note that for even rather small len(x), the total number of
198 permutations of x is larger than the period of most random number
199 generators; this implies that "most" permutations of a long
200 sequence can never be generated.
201 """
202
203 if random is None:
204 random = self.random
205 for i in xrange(len(x)-1, 0, -1):
Tim Peterscd804102001-01-25 20:25:57 +0000206 # pick an element in x[:i+1] with which to exchange x[i]
Tim Petersd7b5e882001-01-25 03:36:26 +0000207 j = int(random() * (i+1))
208 x[i], x[j] = x[j], x[i]
209
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000210 def sample(self, population, k, int=int):
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000211 """Chooses k unique random elements from a population sequence.
212
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000213 Returns a new list containing elements from the population while
214 leaving the original population unchanged. The resulting list is
215 in selection order so that all sub-slices will also be valid random
216 samples. This allows raffle winners (the sample) to be partitioned
217 into grand prize and second place winners (the subslices).
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000218
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000219 Members of the population need not be hashable or unique. If the
220 population contains repeats, then each occurrence is a possible
221 selection in the sample.
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000222
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000223 To choose a sample in a range of integers, use xrange as an argument.
224 This is especially fast and space efficient for sampling from a
225 large population: sample(xrange(10000000), 60)
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000226 """
227
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000228 # Sampling without replacement entails tracking either potential
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000229 # selections (the pool) in a list or previous selections in a
230 # dictionary.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000231
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000232 # When the number of selections is small compared to the population,
233 # then tracking selections is efficient, requiring only a small
234 # dictionary and an occasional reselection. For a larger number of
235 # selections, the pool tracking method is preferred since the list takes
236 # less space than the dictionary and it doesn't suffer from frequent
237 # reselections.
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000238
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000239 n = len(population)
240 if not 0 <= k <= n:
241 raise ValueError, "sample larger than population"
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000242 random = self.random
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000243 result = [None] * k
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000244 if n < 6 * k: # if n len list takes less space than a k len dict
Raymond Hettinger311f4192002-11-18 09:01:24 +0000245 pool = list(population)
246 for i in xrange(k): # invariant: non-selected at [0,n-i)
247 j = int(random() * (n-i))
248 result[i] = pool[j]
Raymond Hettinger8b9aa8d2003-01-04 05:20:33 +0000249 pool[j] = pool[n-i-1] # move non-selected item into vacancy
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000250 else:
Raymond Hettinger311f4192002-11-18 09:01:24 +0000251 selected = {}
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000252 for i in xrange(k):
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000253 j = int(random() * n)
Raymond Hettinger311f4192002-11-18 09:01:24 +0000254 while j in selected:
Raymond Hettingerc0b40342002-11-13 15:26:37 +0000255 j = int(random() * n)
256 result[i] = selected[j] = population[j]
Raymond Hettinger311f4192002-11-18 09:01:24 +0000257 return result
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000258
Tim Peterscd804102001-01-25 20:25:57 +0000259## -------------------- real-valued distributions -------------------
260
261## -------------------- uniform distribution -------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000262
263 def uniform(self, a, b):
264 """Get a random number in the range [a, b)."""
265 return a + (b-a) * self.random()
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000266
Tim Peterscd804102001-01-25 20:25:57 +0000267## -------------------- normal distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000268
Tim Petersd7b5e882001-01-25 03:36:26 +0000269 def normalvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000270 """Normal distribution.
271
272 mu is the mean, and sigma is the standard deviation.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000273
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000274 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000275 # mu = mean, sigma = standard deviation
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000276
Tim Petersd7b5e882001-01-25 03:36:26 +0000277 # Uses Kinderman and Monahan method. Reference: Kinderman,
278 # A.J. and Monahan, J.F., "Computer generation of random
279 # variables using the ratio of uniform deviates", ACM Trans
280 # Math Software, 3, (1977), pp257-260.
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000281
Tim Petersd7b5e882001-01-25 03:36:26 +0000282 random = self.random
Raymond Hettinger311f4192002-11-18 09:01:24 +0000283 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000284 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000285 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000286 z = NV_MAGICCONST*(u1-0.5)/u2
287 zz = z*z/4.0
288 if zz <= -_log(u2):
289 break
290 return mu + z*sigma
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000291
Tim Peterscd804102001-01-25 20:25:57 +0000292## -------------------- lognormal distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000293
294 def lognormvariate(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000295 """Log normal distribution.
296
297 If you take the natural logarithm of this distribution, you'll get a
298 normal distribution with mean mu and standard deviation sigma.
299 mu can have any value, and sigma must be greater than zero.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000300
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000301 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000302 return _exp(self.normalvariate(mu, sigma))
303
Tim Peterscd804102001-01-25 20:25:57 +0000304## -------------------- circular uniform --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000305
306 def cunifvariate(self, mean, arc):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000307 """Circular uniform distribution.
308
309 mean is the mean angle, and arc is the range of the distribution,
310 centered around the mean angle. Both values must be expressed in
311 radians. Returned values range between mean - arc/2 and
312 mean + arc/2 and are normalized to between 0 and pi.
313
314 Deprecated in version 2.3. Use:
315 (mean + arc * (Random.random() - 0.5)) % Math.pi
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000316
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000317 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000318 # mean: mean angle (in radians between 0 and pi)
319 # arc: range of distribution (in radians between 0 and pi)
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000320 import warnings
321 warnings.warn("The cunifvariate function is deprecated; Use (mean "
322 "+ arc * (Random.random() - 0.5)) % Math.pi instead",
323 DeprecationWarning)
Tim Petersd7b5e882001-01-25 03:36:26 +0000324
325 return (mean + arc * (self.random() - 0.5)) % _pi
326
Tim Peterscd804102001-01-25 20:25:57 +0000327## -------------------- exponential distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000328
329 def expovariate(self, lambd):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000330 """Exponential distribution.
331
332 lambd is 1.0 divided by the desired mean. (The parameter would be
333 called "lambda", but that is a reserved word in Python.) Returned
334 values range from 0 to positive infinity.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000335
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000336 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000337 # lambd: rate lambd = 1/mean
338 # ('lambda' is a Python reserved word)
339
340 random = self.random
Tim Peters0c9886d2001-01-15 01:18:21 +0000341 u = random()
342 while u <= 1e-7:
343 u = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000344 return -_log(u)/lambd
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000345
Tim Peterscd804102001-01-25 20:25:57 +0000346## -------------------- von Mises distribution --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000347
Tim Petersd7b5e882001-01-25 03:36:26 +0000348 def vonmisesvariate(self, mu, kappa):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000349 """Circular data distribution.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000350
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000351 mu is the mean angle, expressed in radians between 0 and 2*pi, and
352 kappa is the concentration parameter, which must be greater than or
353 equal to zero. If kappa is equal to zero, this distribution reduces
354 to a uniform random angle over the range 0 to 2*pi.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000355
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000356 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000357 # mu: mean angle (in radians between 0 and 2*pi)
358 # kappa: concentration parameter kappa (>= 0)
359 # if kappa = 0 generate uniform random angle
360
361 # Based upon an algorithm published in: Fisher, N.I.,
362 # "Statistical Analysis of Circular Data", Cambridge
363 # University Press, 1993.
364
365 # Thanks to Magnus Kessler for a correction to the
366 # implementation of step 4.
367
368 random = self.random
369 if kappa <= 1e-6:
370 return TWOPI * random()
371
372 a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
373 b = (a - _sqrt(2.0 * a))/(2.0 * kappa)
374 r = (1.0 + b * b)/(2.0 * b)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000375
Raymond Hettinger311f4192002-11-18 09:01:24 +0000376 while True:
Tim Peters0c9886d2001-01-15 01:18:21 +0000377 u1 = random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000378
379 z = _cos(_pi * u1)
380 f = (1.0 + r * z)/(r + z)
381 c = kappa * (r - f)
382
383 u2 = random()
384
385 if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
Tim Peters0c9886d2001-01-15 01:18:21 +0000386 break
Tim Petersd7b5e882001-01-25 03:36:26 +0000387
388 u3 = random()
389 if u3 > 0.5:
390 theta = (mu % TWOPI) + _acos(f)
391 else:
392 theta = (mu % TWOPI) - _acos(f)
393
394 return theta
395
Tim Peterscd804102001-01-25 20:25:57 +0000396## -------------------- gamma distribution --------------------
Tim Petersd7b5e882001-01-25 03:36:26 +0000397
398 def gammavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000399 """Gamma distribution. Not the gamma function!
400
401 Conditions on the parameters are alpha > 0 and beta > 0.
402
403 """
Tim Peters8ac14952002-05-23 15:15:30 +0000404
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000405 # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
Tim Peters8ac14952002-05-23 15:15:30 +0000406
Guido van Rossum570764d2002-05-14 14:08:12 +0000407 # Warning: a few older sources define the gamma distribution in terms
408 # of alpha > -1.0
409 if alpha <= 0.0 or beta <= 0.0:
410 raise ValueError, 'gammavariate: alpha and beta must be > 0.0'
Tim Peters8ac14952002-05-23 15:15:30 +0000411
Tim Petersd7b5e882001-01-25 03:36:26 +0000412 random = self.random
Tim Petersd7b5e882001-01-25 03:36:26 +0000413 if alpha > 1.0:
414
415 # Uses R.C.H. Cheng, "The generation of Gamma
416 # variables with non-integral shape parameters",
417 # Applied Statistics, (1977), 26, No. 1, p71-74
418
Raymond Hettingerca6cdc22002-05-13 23:40:14 +0000419 ainv = _sqrt(2.0 * alpha - 1.0)
420 bbb = alpha - LOG4
421 ccc = alpha + ainv
Tim Peters8ac14952002-05-23 15:15:30 +0000422
Raymond Hettinger311f4192002-11-18 09:01:24 +0000423 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000424 u1 = random()
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000425 if not 1e-7 < u1 < .9999999:
426 continue
427 u2 = 1.0 - random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000428 v = _log(u1/(1.0-u1))/ainv
429 x = alpha*_exp(v)
430 z = u1*u1*u2
431 r = bbb+ccc*v-x
432 if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z):
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000433 return x * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000434
435 elif alpha == 1.0:
436 # expovariate(1)
437 u = random()
438 while u <= 1e-7:
439 u = random()
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000440 return -_log(u) * beta
Tim Petersd7b5e882001-01-25 03:36:26 +0000441
442 else: # alpha is between 0 and 1 (exclusive)
443
444 # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
445
Raymond Hettinger311f4192002-11-18 09:01:24 +0000446 while True:
Tim Petersd7b5e882001-01-25 03:36:26 +0000447 u = random()
448 b = (_e + alpha)/_e
449 p = b*u
450 if p <= 1.0:
451 x = pow(p, 1.0/alpha)
452 else:
453 # p > 1
454 x = -_log((b-p)/alpha)
455 u1 = random()
456 if not (((p <= 1.0) and (u1 > _exp(-x))) or
457 ((p > 1) and (u1 > pow(x, alpha - 1.0)))):
458 break
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000459 return x * beta
460
461
462 def stdgamma(self, alpha, ainv, bbb, ccc):
463 # This method was (and shall remain) undocumented.
464 # This method is deprecated
465 # for the following reasons:
466 # 1. Returns same as .gammavariate(alpha, 1.0)
467 # 2. Requires caller to provide 3 extra arguments
468 # that are functions of alpha anyway
469 # 3. Can't be used for alpha < 0.5
470
471 # ainv = sqrt(2 * alpha - 1)
472 # bbb = alpha - log(4)
473 # ccc = alpha + ainv
474 import warnings
475 warnings.warn("The stdgamma function is deprecated; "
476 "use gammavariate() instead",
477 DeprecationWarning)
478 return self.gammavariate(alpha, 1.0)
479
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000480
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000481
Tim Peterscd804102001-01-25 20:25:57 +0000482## -------------------- Gauss (faster alternative) --------------------
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000483
Tim Petersd7b5e882001-01-25 03:36:26 +0000484 def gauss(self, mu, sigma):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000485 """Gaussian distribution.
486
487 mu is the mean, and sigma is the standard deviation. This is
488 slightly faster than the normalvariate() function.
489
490 Not thread-safe without a lock around calls.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000491
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000492 """
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000493
Tim Petersd7b5e882001-01-25 03:36:26 +0000494 # When x and y are two variables from [0, 1), uniformly
495 # distributed, then
496 #
497 # cos(2*pi*x)*sqrt(-2*log(1-y))
498 # sin(2*pi*x)*sqrt(-2*log(1-y))
499 #
500 # are two *independent* variables with normal distribution
501 # (mu = 0, sigma = 1).
502 # (Lambert Meertens)
503 # (corrected version; bug discovered by Mike Miller, fixed by LM)
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000504
Tim Petersd7b5e882001-01-25 03:36:26 +0000505 # Multithreading note: When two threads call this function
506 # simultaneously, it is possible that they will receive the
507 # same return value. The window is very small though. To
508 # avoid this, you have to use a lock around all calls. (I
509 # didn't want to slow this down in the serial case by using a
510 # lock here.)
Guido van Rossumd03e1191998-05-29 17:51:31 +0000511
Tim Petersd7b5e882001-01-25 03:36:26 +0000512 random = self.random
513 z = self.gauss_next
514 self.gauss_next = None
515 if z is None:
516 x2pi = random() * TWOPI
517 g2rad = _sqrt(-2.0 * _log(1.0 - random()))
518 z = _cos(x2pi) * g2rad
519 self.gauss_next = _sin(x2pi) * g2rad
Guido van Rossumcc32ac91994-03-15 16:10:24 +0000520
Tim Petersd7b5e882001-01-25 03:36:26 +0000521 return mu + z*sigma
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000522
Tim Peterscd804102001-01-25 20:25:57 +0000523## -------------------- beta --------------------
Tim Peters85e2e472001-01-26 06:49:56 +0000524## See
525## http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
526## for Ivan Frohne's insightful analysis of why the original implementation:
527##
528## def betavariate(self, alpha, beta):
529## # Discrete Event Simulation in C, pp 87-88.
530##
531## y = self.expovariate(alpha)
532## z = self.expovariate(1.0/beta)
533## return z/(y+z)
534##
535## was dead wrong, and how it probably got that way.
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000536
Tim Petersd7b5e882001-01-25 03:36:26 +0000537 def betavariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000538 """Beta distribution.
539
540 Conditions on the parameters are alpha > -1 and beta} > -1.
541 Returned values range between 0 and 1.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000542
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000543 """
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000544
Tim Peters85e2e472001-01-26 06:49:56 +0000545 # This version due to Janne Sinkkonen, and matches all the std
546 # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
547 y = self.gammavariate(alpha, 1.)
548 if y == 0:
549 return 0.0
550 else:
551 return y / (y + self.gammavariate(beta, 1.))
Guido van Rossum95bfcda1994-03-09 14:21:05 +0000552
Tim Peterscd804102001-01-25 20:25:57 +0000553## -------------------- Pareto --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000554
Tim Petersd7b5e882001-01-25 03:36:26 +0000555 def paretovariate(self, alpha):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000556 """Pareto distribution. alpha is the shape parameter."""
Tim Petersd7b5e882001-01-25 03:36:26 +0000557 # Jain, pg. 495
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000558
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000559 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000560 return 1.0 / pow(u, 1.0/alpha)
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000561
Tim Peterscd804102001-01-25 20:25:57 +0000562## -------------------- Weibull --------------------
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000563
Tim Petersd7b5e882001-01-25 03:36:26 +0000564 def weibullvariate(self, alpha, beta):
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000565 """Weibull distribution.
566
567 alpha is the scale parameter and beta is the shape parameter.
Raymond Hettingeref4d4bd2002-05-23 23:58:17 +0000568
Raymond Hettingerc32f0332002-05-23 19:44:49 +0000569 """
Tim Petersd7b5e882001-01-25 03:36:26 +0000570 # Jain, pg. 499; bug fix courtesy Bill Arms
Guido van Rossumcf4559a1997-12-02 02:47:39 +0000571
Raymond Hettinger73ced7e2003-01-04 09:26:32 +0000572 u = 1.0 - self.random()
Tim Petersd7b5e882001-01-25 03:36:26 +0000573 return alpha * pow(-_log(u), 1.0/beta)
Guido van Rossum6c395ba1999-08-18 13:53:28 +0000574
Raymond Hettinger40f62172002-12-29 23:03:38 +0000575## -------------------- Wichmann-Hill -------------------
576
577class WichmannHill(Random):
578
579 VERSION = 1 # used by getstate/setstate
580
581 def seed(self, a=None):
582 """Initialize internal state from hashable object.
583
584 None or no argument seeds from current time.
585
586 If a is not None or an int or long, hash(a) is used instead.
587
588 If a is an int or long, a is used directly. Distinct values between
589 0 and 27814431486575L inclusive are guaranteed to yield distinct
590 internal states (this guarantee is specific to the default
591 Wichmann-Hill generator).
592 """
593
594 if a is None:
595 # Initialize from current time
596 import time
597 a = long(time.time() * 256)
598
599 if not isinstance(a, (int, long)):
600 a = hash(a)
601
602 a, x = divmod(a, 30268)
603 a, y = divmod(a, 30306)
604 a, z = divmod(a, 30322)
605 self._seed = int(x)+1, int(y)+1, int(z)+1
606
607 self.gauss_next = None
608
609 def random(self):
610 """Get the next random number in the range [0.0, 1.0)."""
611
612 # Wichman-Hill random number generator.
613 #
614 # Wichmann, B. A. & Hill, I. D. (1982)
615 # Algorithm AS 183:
616 # An efficient and portable pseudo-random number generator
617 # Applied Statistics 31 (1982) 188-190
618 #
619 # see also:
620 # Correction to Algorithm AS 183
621 # Applied Statistics 33 (1984) 123
622 #
623 # McLeod, A. I. (1985)
624 # A remark on Algorithm AS 183
625 # Applied Statistics 34 (1985),198-200
626
627 # This part is thread-unsafe:
628 # BEGIN CRITICAL SECTION
629 x, y, z = self._seed
630 x = (171 * x) % 30269
631 y = (172 * y) % 30307
632 z = (170 * z) % 30323
633 self._seed = x, y, z
634 # END CRITICAL SECTION
635
636 # Note: on a platform using IEEE-754 double arithmetic, this can
637 # never return 0.0 (asserted by Tim; proof too long for a comment).
638 return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
639
640 def getstate(self):
641 """Return internal state; can be passed to setstate() later."""
642 return self.VERSION, self._seed, self.gauss_next
643
644 def setstate(self, state):
645 """Restore internal state from object returned by getstate()."""
646 version = state[0]
647 if version == 1:
648 version, self._seed, self.gauss_next = state
649 else:
650 raise ValueError("state with version %s passed to "
651 "Random.setstate() of version %s" %
652 (version, self.VERSION))
653
654 def jumpahead(self, n):
655 """Act as if n calls to random() were made, but quickly.
656
657 n is an int, greater than or equal to 0.
658
659 Example use: If you have 2 threads and know that each will
660 consume no more than a million random numbers, create two Random
661 objects r1 and r2, then do
662 r2.setstate(r1.getstate())
663 r2.jumpahead(1000000)
664 Then r1 and r2 will use guaranteed-disjoint segments of the full
665 period.
666 """
667
668 if not n >= 0:
669 raise ValueError("n must be >= 0")
670 x, y, z = self._seed
671 x = int(x * pow(171, n, 30269)) % 30269
672 y = int(y * pow(172, n, 30307)) % 30307
673 z = int(z * pow(170, n, 30323)) % 30323
674 self._seed = x, y, z
675
676 def __whseed(self, x=0, y=0, z=0):
677 """Set the Wichmann-Hill seed from (x, y, z).
678
679 These must be integers in the range [0, 256).
680 """
681
682 if not type(x) == type(y) == type(z) == int:
683 raise TypeError('seeds must be integers')
684 if not (0 <= x < 256 and 0 <= y < 256 and 0 <= z < 256):
685 raise ValueError('seeds must be in range(0, 256)')
686 if 0 == x == y == z:
687 # Initialize from current time
688 import time
689 t = long(time.time() * 256)
690 t = int((t&0xffffff) ^ (t>>24))
691 t, x = divmod(t, 256)
692 t, y = divmod(t, 256)
693 t, z = divmod(t, 256)
694 # Zero is a poor seed, so substitute 1
695 self._seed = (x or 1, y or 1, z or 1)
696
697 self.gauss_next = None
698
699 def whseed(self, a=None):
700 """Seed from hashable object's hash code.
701
702 None or no argument seeds from current time. It is not guaranteed
703 that objects with distinct hash codes lead to distinct internal
704 states.
705
706 This is obsolete, provided for compatibility with the seed routine
707 used prior to Python 2.1. Use the .seed() method instead.
708 """
709
710 if a is None:
711 self.__whseed()
712 return
713 a = hash(a)
714 a, x = divmod(a, 256)
715 a, y = divmod(a, 256)
716 a, z = divmod(a, 256)
717 x = (x + a) % 256 or 1
718 y = (y + a) % 256 or 1
719 z = (z + a) % 256 or 1
720 self.__whseed(x, y, z)
721
Tim Peterscd804102001-01-25 20:25:57 +0000722## -------------------- test program --------------------
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000723
Tim Petersd7b5e882001-01-25 03:36:26 +0000724def _test_generator(n, funccall):
Tim Peters0c9886d2001-01-15 01:18:21 +0000725 import time
726 print n, 'times', funccall
727 code = compile(funccall, funccall, 'eval')
728 sum = 0.0
729 sqsum = 0.0
730 smallest = 1e10
731 largest = -1e10
732 t0 = time.time()
733 for i in range(n):
734 x = eval(code)
735 sum = sum + x
736 sqsum = sqsum + x*x
737 smallest = min(x, smallest)
738 largest = max(x, largest)
739 t1 = time.time()
740 print round(t1-t0, 3), 'sec,',
741 avg = sum/n
Tim Petersd7b5e882001-01-25 03:36:26 +0000742 stddev = _sqrt(sqsum/n - avg*avg)
Tim Peters0c9886d2001-01-15 01:18:21 +0000743 print 'avg %g, stddev %g, min %g, max %g' % \
744 (avg, stddev, smallest, largest)
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000745
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000746def _sample_generator(n, k):
747 # Return a fixed element from the sample. Validates random ordering.
748 return sample(xrange(n), k)[k//2]
749
750def _test(N=2000):
Tim Petersd7b5e882001-01-25 03:36:26 +0000751 _test_generator(N, 'random()')
752 _test_generator(N, 'normalvariate(0.0, 1.0)')
753 _test_generator(N, 'lognormvariate(0.0, 1.0)')
754 _test_generator(N, 'cunifvariate(0.0, 1.0)')
Tim Petersd7b5e882001-01-25 03:36:26 +0000755 _test_generator(N, 'vonmisesvariate(0.0, 1.0)')
Raymond Hettingerb760efb2002-05-14 06:40:34 +0000756 _test_generator(N, 'gammavariate(0.01, 1.0)')
757 _test_generator(N, 'gammavariate(0.1, 1.0)')
Tim Peters8ac14952002-05-23 15:15:30 +0000758 _test_generator(N, 'gammavariate(0.1, 2.0)')
Tim Petersd7b5e882001-01-25 03:36:26 +0000759 _test_generator(N, 'gammavariate(0.5, 1.0)')
760 _test_generator(N, 'gammavariate(0.9, 1.0)')
761 _test_generator(N, 'gammavariate(1.0, 1.0)')
762 _test_generator(N, 'gammavariate(2.0, 1.0)')
763 _test_generator(N, 'gammavariate(20.0, 1.0)')
764 _test_generator(N, 'gammavariate(200.0, 1.0)')
765 _test_generator(N, 'gauss(0.0, 1.0)')
766 _test_generator(N, 'betavariate(3.0, 3.0)')
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000767 _test_generator(N, '_sample_generator(50, 5)') # expected s.d.: 14.4
768 _test_generator(N, '_sample_generator(50, 45)') # expected s.d.: 14.4
Tim Peterscd804102001-01-25 20:25:57 +0000769
Tim Peters715c4c42001-01-26 22:56:56 +0000770# Create one instance, seeded from current time, and export its methods
Raymond Hettinger40f62172002-12-29 23:03:38 +0000771# as module-level functions. The functions share state across all uses
772#(both in the user's code and in the Python libraries), but that's fine
773# for most programs and is easier for the casual user than making them
774# instantiate their own Random() instance.
775
Tim Petersd7b5e882001-01-25 03:36:26 +0000776_inst = Random()
777seed = _inst.seed
778random = _inst.random
779uniform = _inst.uniform
780randint = _inst.randint
781choice = _inst.choice
782randrange = _inst.randrange
Raymond Hettingerf24eb352002-11-12 17:41:57 +0000783sample = _inst.sample
Tim Petersd7b5e882001-01-25 03:36:26 +0000784shuffle = _inst.shuffle
785normalvariate = _inst.normalvariate
786lognormvariate = _inst.lognormvariate
787cunifvariate = _inst.cunifvariate
788expovariate = _inst.expovariate
789vonmisesvariate = _inst.vonmisesvariate
790gammavariate = _inst.gammavariate
791stdgamma = _inst.stdgamma
792gauss = _inst.gauss
793betavariate = _inst.betavariate
794paretovariate = _inst.paretovariate
795weibullvariate = _inst.weibullvariate
796getstate = _inst.getstate
797setstate = _inst.setstate
Tim Petersd52269b2001-01-25 06:23:18 +0000798jumpahead = _inst.jumpahead
Tim Petersd7b5e882001-01-25 03:36:26 +0000799
Guido van Rossumff03b1a1994-03-09 12:55:02 +0000800if __name__ == '__main__':
Tim Petersd7b5e882001-01-25 03:36:26 +0000801 _test()