blob: f37bc2a111d9540d65e06e2e14c42cb3002212d7 [file] [log] [blame]
Georg Brandl116aa622007-08-15 14:28:22 +00001:mod:`random` --- Generate pseudo-random numbers
2================================================
3
4.. module:: random
5 :synopsis: Generate pseudo-random numbers with various common distributions.
6
Raymond Hettinger10480942011-01-10 03:26:08 +00007**Source code:** :source:`Lib/random.py`
Georg Brandl116aa622007-08-15 14:28:22 +00008
Raymond Hettinger4f707fd2011-01-10 19:54:11 +00009--------------
10
Georg Brandl116aa622007-08-15 14:28:22 +000011This module implements pseudo-random number generators for various
12distributions.
13
Raymond Hettingerb21dac12010-09-07 05:32:49 +000014For integers, there is uniform selection from a range. For sequences, there is
15uniform selection of a random element, a function to generate a random
16permutation of a list in-place, and a function for random sampling without
17replacement.
Georg Brandl116aa622007-08-15 14:28:22 +000018
19On the real line, there are functions to compute uniform, normal (Gaussian),
20lognormal, negative exponential, gamma, and beta distributions. For generating
21distributions of angles, the von Mises distribution is available.
22
Georg Brandl92849d12016-02-19 08:57:38 +010023Almost all module functions depend on the basic function :func:`.random`, which
Georg Brandl116aa622007-08-15 14:28:22 +000024generates a random float uniformly in the semi-open range [0.0, 1.0). Python
25uses the Mersenne Twister as the core generator. It produces 53-bit precision
26floats and has a period of 2\*\*19937-1. The underlying implementation in C is
27both fast and threadsafe. The Mersenne Twister is one of the most extensively
28tested random number generators in existence. However, being completely
29deterministic, it is not suitable for all purposes, and is completely unsuitable
30for cryptographic purposes.
31
32The functions supplied by this module are actually bound methods of a hidden
33instance of the :class:`random.Random` class. You can instantiate your own
Raymond Hettinger28de64f2008-01-13 23:40:30 +000034instances of :class:`Random` to get generators that don't share state.
Georg Brandl116aa622007-08-15 14:28:22 +000035
36Class :class:`Random` can also be subclassed if you want to use a different
Georg Brandl92849d12016-02-19 08:57:38 +010037basic generator of your own devising: in that case, override the :meth:`~Random.random`,
38:meth:`~Random.seed`, :meth:`~Random.getstate`, and :meth:`~Random.setstate` methods.
39Optionally, a new generator can supply a :meth:`~Random.getrandbits` method --- this
Georg Brandl116aa622007-08-15 14:28:22 +000040allows :meth:`randrange` to produce selections over an arbitrarily large range.
41
Benjamin Peterson21896a32010-03-21 22:03:03 +000042The :mod:`random` module also provides the :class:`SystemRandom` class which
43uses the system function :func:`os.urandom` to generate random numbers
44from sources provided by the operating system.
Georg Brandl116aa622007-08-15 14:28:22 +000045
Raymond Hettingerc89a4512014-05-11 02:26:23 -070046.. warning::
47
48 The pseudo-random generators of this module should not be used for
Steven D'Apranob2871fa2016-04-17 01:42:33 +100049 security purposes. For security or cryptographic uses, see the
50 :mod:`secrets` module.
Raymond Hettingerc89a4512014-05-11 02:26:23 -070051
Raymond Hettingere1329102016-11-21 12:33:50 -080052.. seealso::
Georg Brandl116aa622007-08-15 14:28:22 +000053
Raymond Hettingere1329102016-11-21 12:33:50 -080054 M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-dimensionally
55 equidistributed uniform pseudorandom number generator", ACM Transactions on
Serhiy Storchaka0264e462016-11-26 13:49:59 +020056 Modeling and Computer Simulation Vol. 8, No. 1, January pp.3--30 1998.
Raymond Hettingere1329102016-11-21 12:33:50 -080057
58
59 `Complementary-Multiply-with-Carry recipe
60 <https://code.activestate.com/recipes/576707/>`_ for a compatible alternative
61 random number generator with a long period and comparatively simple update
62 operations.
63
64
65Bookkeeping functions
66---------------------
Georg Brandl116aa622007-08-15 14:28:22 +000067
Ezio Melottie0add762012-09-14 06:32:35 +030068.. function:: seed(a=None, version=2)
Georg Brandl116aa622007-08-15 14:28:22 +000069
Raymond Hettingerf763a722010-09-07 00:38:15 +000070 Initialize the random number generator.
Georg Brandl116aa622007-08-15 14:28:22 +000071
Ezio Melottie0add762012-09-14 06:32:35 +030072 If *a* is omitted or ``None``, the current system time is used. If
Raymond Hettingerf763a722010-09-07 00:38:15 +000073 randomness sources are provided by the operating system, they are used
74 instead of the system time (see the :func:`os.urandom` function for details
75 on availability).
Georg Brandl116aa622007-08-15 14:28:22 +000076
Ezio Melottie0add762012-09-14 06:32:35 +030077 If *a* is an int, it is used directly.
Raymond Hettingerf763a722010-09-07 00:38:15 +000078
79 With version 2 (the default), a :class:`str`, :class:`bytes`, or :class:`bytearray`
Raymond Hettinger16eb8272016-09-04 11:17:28 -070080 object gets converted to an :class:`int` and all of its bits are used.
81
82 With version 1 (provided for reproducing random sequences from older versions
83 of Python), the algorithm for :class:`str` and :class:`bytes` generates a
84 narrower range of seeds.
Raymond Hettingerf763a722010-09-07 00:38:15 +000085
86 .. versionchanged:: 3.2
87 Moved to the version 2 scheme which uses all of the bits in a string seed.
Georg Brandl116aa622007-08-15 14:28:22 +000088
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -070089 .. deprecated:: 3.9
90 In the future, the *seed* must be one of the following types:
91 *NoneType*, :class:`int`, :class:`float`, :class:`str`,
92 :class:`bytes`, or :class:`bytearray`.
93
Georg Brandl116aa622007-08-15 14:28:22 +000094.. function:: getstate()
95
96 Return an object capturing the current internal state of the generator. This
97 object can be passed to :func:`setstate` to restore the state.
98
Georg Brandl116aa622007-08-15 14:28:22 +000099
100.. function:: setstate(state)
101
102 *state* should have been obtained from a previous call to :func:`getstate`, and
103 :func:`setstate` restores the internal state of the generator to what it was at
Sandro Tosi985104a2012-08-12 15:12:15 +0200104 the time :func:`getstate` was called.
Georg Brandl116aa622007-08-15 14:28:22 +0000105
Georg Brandl116aa622007-08-15 14:28:22 +0000106
Raymond Hettingerf01d1be2020-05-04 22:52:13 -0700107Functions for bytes
108-------------------
Georg Brandl116aa622007-08-15 14:28:22 +0000109
Victor Stinner9f5fe792020-04-17 19:05:35 +0200110.. function:: randbytes(n)
111
112 Generate *n* random bytes.
113
Raymond Hettingerf01d1be2020-05-04 22:52:13 -0700114 This method should not be used for generating security tokens.
115 Use :func:`secrets.token_bytes` instead.
116
Victor Stinner9f5fe792020-04-17 19:05:35 +0200117 .. versionadded:: 3.9
118
119
Raymond Hettingere1329102016-11-21 12:33:50 -0800120Functions for integers
121----------------------
Georg Brandl116aa622007-08-15 14:28:22 +0000122
Ezio Melottie0add762012-09-14 06:32:35 +0300123.. function:: randrange(stop)
124 randrange(start, stop[, step])
Georg Brandl116aa622007-08-15 14:28:22 +0000125
126 Return a randomly selected element from ``range(start, stop, step)``. This is
127 equivalent to ``choice(range(start, stop, step))``, but doesn't actually build a
128 range object.
129
Raymond Hettinger05156612010-09-07 04:44:52 +0000130 The positional argument pattern matches that of :func:`range`. Keyword arguments
131 should not be used because the function may use them in unexpected ways.
132
133 .. versionchanged:: 3.2
134 :meth:`randrange` is more sophisticated about producing equally distributed
135 values. Formerly it used a style like ``int(random()*n)`` which could produce
136 slightly uneven distributions.
Georg Brandl116aa622007-08-15 14:28:22 +0000137
138.. function:: randint(a, b)
139
Raymond Hettingerafd30452009-02-24 10:57:02 +0000140 Return a random integer *N* such that ``a <= N <= b``. Alias for
141 ``randrange(a, b+1)``.
Georg Brandl116aa622007-08-15 14:28:22 +0000142
Raymond Hettingerf01d1be2020-05-04 22:52:13 -0700143.. function:: getrandbits(k)
144
145 Returns a Python integer with *k* random bits. This method is supplied with
146 the MersenneTwister generator and some other generators may also provide it
147 as an optional part of the API. When available, :meth:`getrandbits` enables
148 :meth:`randrange` to handle arbitrarily large ranges.
149
150 .. versionchanged:: 3.9
151 This method now accepts zero for *k*.
152
Georg Brandl116aa622007-08-15 14:28:22 +0000153
Raymond Hettingere1329102016-11-21 12:33:50 -0800154Functions for sequences
155-----------------------
Georg Brandl116aa622007-08-15 14:28:22 +0000156
157.. function:: choice(seq)
158
159 Return a random element from the non-empty sequence *seq*. If *seq* is empty,
160 raises :exc:`IndexError`.
161
Raymond Hettinger9016f282016-09-26 21:45:57 -0700162.. function:: choices(population, weights=None, *, cum_weights=None, k=1)
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700163
164 Return a *k* sized list of elements chosen from the *population* with replacement.
165 If the *population* is empty, raises :exc:`IndexError`.
166
167 If a *weights* sequence is specified, selections are made according to the
168 relative weights. Alternatively, if a *cum_weights* sequence is given, the
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400169 selections are made according to the cumulative weights (perhaps computed
170 using :func:`itertools.accumulate`). For example, the relative weights
171 ``[10, 5, 30, 5]`` are equivalent to the cumulative weights
172 ``[10, 15, 45, 50]``. Internally, the relative weights are converted to
173 cumulative weights before making selections, so supplying the cumulative
174 weights saves work.
Raymond Hettingere8f1e002016-09-06 17:15:29 -0700175
176 If neither *weights* nor *cum_weights* are specified, selections are made
177 with equal probability. If a weights sequence is supplied, it must be
178 the same length as the *population* sequence. It is a :exc:`TypeError`
179 to specify both *weights* and *cum_weights*.
180
181 The *weights* or *cum_weights* can use any numeric type that interoperates
182 with the :class:`float` values returned by :func:`random` (that includes
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800183 integers, floats, and fractions but excludes decimals). Behavior is
184 undefined if any weight is negative. A :exc:`ValueError` is raised if all
185 weights are zero.
Georg Brandl116aa622007-08-15 14:28:22 +0000186
Raymond Hettinger40ebe942019-01-30 13:30:20 -0800187 For a given seed, the :func:`choices` function with equal weighting
188 typically produces a different sequence than repeated calls to
189 :func:`choice`. The algorithm used by :func:`choices` uses floating
190 point arithmetic for internal consistency and speed. The algorithm used
191 by :func:`choice` defaults to integer arithmetic with repeated selections
192 to avoid small biases from round-off error.
193
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400194 .. versionadded:: 3.6
195
Raymond Hettinger041d8b42019-11-23 02:22:13 -0800196 .. versionchanged:: 3.9
197 Raises a :exc:`ValueError` if all weights are zero.
198
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400199
Georg Brandl116aa622007-08-15 14:28:22 +0000200.. function:: shuffle(x[, random])
201
Raymond Hettingera3950e42016-11-17 01:49:54 -0800202 Shuffle the sequence *x* in place.
Georg Brandl116aa622007-08-15 14:28:22 +0000203
Raymond Hettingera3950e42016-11-17 01:49:54 -0800204 The optional argument *random* is a 0-argument function returning a random
205 float in [0.0, 1.0); by default, this is the function :func:`.random`.
206
207 To shuffle an immutable sequence and return a new shuffled list, use
208 ``sample(x, k=len(x))`` instead.
209
210 Note that even for small ``len(x)``, the total number of permutations of *x*
211 can quickly grow larger than the period of most random number generators.
212 This implies that most permutations of a long sequence can never be
213 generated. For example, a sequence of length 2080 is the largest that
214 can fit within the period of the Mersenne Twister random number generator.
Georg Brandl116aa622007-08-15 14:28:22 +0000215
Raymond Hettinger190fac92020-05-02 16:45:32 -0700216 .. deprecated-removed:: 3.9 3.11
217 The optional parameter *random*.
218
Georg Brandl116aa622007-08-15 14:28:22 +0000219
220.. function:: sample(population, k)
221
Raymond Hettinger1acde192008-01-14 01:00:53 +0000222 Return a *k* length list of unique elements chosen from the population sequence
223 or set. Used for random sampling without replacement.
Georg Brandl116aa622007-08-15 14:28:22 +0000224
Georg Brandl116aa622007-08-15 14:28:22 +0000225 Returns a new list containing elements from the population while leaving the
226 original population unchanged. The resulting list is in selection order so that
227 all sub-slices will also be valid random samples. This allows raffle winners
228 (the sample) to be partitioned into grand prize and second place winners (the
229 subslices).
230
Guido van Rossum2cc30da2007-11-02 23:46:40 +0000231 Members of the population need not be :term:`hashable` or unique. If the population
Georg Brandl116aa622007-08-15 14:28:22 +0000232 contains repeats, then each occurrence is a possible selection in the sample.
233
Raymond Hettingera3950e42016-11-17 01:49:54 -0800234 To choose a sample from a range of integers, use a :func:`range` object as an
Georg Brandl116aa622007-08-15 14:28:22 +0000235 argument. This is especially fast and space efficient for sampling from a large
Raymond Hettingera3950e42016-11-17 01:49:54 -0800236 population: ``sample(range(10000000), k=60)``.
Georg Brandl116aa622007-08-15 14:28:22 +0000237
Raymond Hettingerf07d9492012-07-09 12:43:57 -0700238 If the sample size is larger than the population size, a :exc:`ValueError`
Raymond Hettinger86a20f82012-07-08 16:01:53 -0700239 is raised.
240
Raymond Hettinger4fe00202020-04-19 00:36:42 -0700241 .. deprecated:: 3.9
242 In the future, the *population* must be a sequence. Instances of
243 :class:`set` are no longer supported. The set must first be converted
244 to a :class:`list` or :class:`tuple`, preferably in a deterministic
245 order so that the sample is reproducible.
246
247
Raymond Hettingere1329102016-11-21 12:33:50 -0800248Real-valued distributions
249-------------------------
250
Georg Brandl116aa622007-08-15 14:28:22 +0000251The following functions generate specific real-valued distributions. Function
252parameters are named after the corresponding variables in the distribution's
253equation, as used in common mathematical practice; most of these equations can
254be found in any statistics text.
255
256
257.. function:: random()
258
259 Return the next random floating point number in the range [0.0, 1.0).
260
261
262.. function:: uniform(a, b)
263
Benjamin Petersonb58dda72009-01-18 22:27:04 +0000264 Return a random floating point number *N* such that ``a <= N <= b`` for
265 ``a <= b`` and ``b <= N <= a`` for ``b < a``.
Georg Brandl116aa622007-08-15 14:28:22 +0000266
Raymond Hettingerbe40db02009-06-11 23:12:14 +0000267 The end-point value ``b`` may or may not be included in the range
268 depending on floating-point rounding in the equation ``a + (b-a) * random()``.
Benjamin Peterson35e8c462008-04-24 02:34:53 +0000269
Georg Brandl73dd7c72011-09-17 20:36:28 +0200270
Christian Heimesfe337bf2008-03-23 21:54:12 +0000271.. function:: triangular(low, high, mode)
272
Benjamin Petersonb58dda72009-01-18 22:27:04 +0000273 Return a random floating point number *N* such that ``low <= N <= high`` and
Christian Heimescc47b052008-03-25 14:56:36 +0000274 with the specified *mode* between those bounds. The *low* and *high* bounds
275 default to zero and one. The *mode* argument defaults to the midpoint
276 between the bounds, giving a symmetric distribution.
Christian Heimesfe337bf2008-03-23 21:54:12 +0000277
Georg Brandl116aa622007-08-15 14:28:22 +0000278
279.. function:: betavariate(alpha, beta)
280
Benjamin Petersonb58dda72009-01-18 22:27:04 +0000281 Beta distribution. Conditions on the parameters are ``alpha > 0`` and
282 ``beta > 0``. Returned values range between 0 and 1.
Georg Brandl116aa622007-08-15 14:28:22 +0000283
284
285.. function:: expovariate(lambd)
286
Mark Dickinson2f947362009-01-07 17:54:07 +0000287 Exponential distribution. *lambd* is 1.0 divided by the desired
288 mean. It should be nonzero. (The parameter would be called
289 "lambda", but that is a reserved word in Python.) Returned values
290 range from 0 to positive infinity if *lambd* is positive, and from
291 negative infinity to 0 if *lambd* is negative.
Georg Brandl116aa622007-08-15 14:28:22 +0000292
293
294.. function:: gammavariate(alpha, beta)
295
Benjamin Petersonb58dda72009-01-18 22:27:04 +0000296 Gamma distribution. (*Not* the gamma function!) Conditions on the
297 parameters are ``alpha > 0`` and ``beta > 0``.
Georg Brandl116aa622007-08-15 14:28:22 +0000298
Georg Brandl73dd7c72011-09-17 20:36:28 +0200299 The probability distribution function is::
300
301 x ** (alpha - 1) * math.exp(-x / beta)
302 pdf(x) = --------------------------------------
303 math.gamma(alpha) * beta ** alpha
304
Georg Brandl116aa622007-08-15 14:28:22 +0000305
306.. function:: gauss(mu, sigma)
307
Benjamin Petersonb58dda72009-01-18 22:27:04 +0000308 Gaussian distribution. *mu* is the mean, and *sigma* is the standard
309 deviation. This is slightly faster than the :func:`normalvariate` function
310 defined below.
Georg Brandl116aa622007-08-15 14:28:22 +0000311
312
313.. function:: lognormvariate(mu, sigma)
314
315 Log normal distribution. If you take the natural logarithm of this
316 distribution, you'll get a normal distribution with mean *mu* and standard
317 deviation *sigma*. *mu* can have any value, and *sigma* must be greater than
318 zero.
319
320
321.. function:: normalvariate(mu, sigma)
322
323 Normal distribution. *mu* is the mean, and *sigma* is the standard deviation.
324
325
326.. function:: vonmisesvariate(mu, kappa)
327
328 *mu* is the mean angle, expressed in radians between 0 and 2\*\ *pi*, and *kappa*
329 is the concentration parameter, which must be greater than or equal to zero. If
330 *kappa* is equal to zero, this distribution reduces to a uniform random angle
331 over the range 0 to 2\*\ *pi*.
332
333
334.. function:: paretovariate(alpha)
335
336 Pareto distribution. *alpha* is the shape parameter.
337
338
339.. function:: weibullvariate(alpha, beta)
340
341 Weibull distribution. *alpha* is the scale parameter and *beta* is the shape
342 parameter.
343
344
Raymond Hettingere1329102016-11-21 12:33:50 -0800345Alternative Generator
346---------------------
Georg Brandl116aa622007-08-15 14:28:22 +0000347
Matthias Bussonnier31e8d692019-04-16 09:47:11 -0700348.. class:: Random([seed])
349
350 Class that implements the default pseudo-random number generator used by the
351 :mod:`random` module.
352
Raymond Hettingerd0cdeaa2019-08-22 09:19:36 -0700353 .. deprecated:: 3.9
354 In the future, the *seed* must be one of the following types:
355 :class:`NoneType`, :class:`int`, :class:`float`, :class:`str`,
356 :class:`bytes`, or :class:`bytearray`.
357
Georg Brandl116aa622007-08-15 14:28:22 +0000358.. class:: SystemRandom([seed])
359
360 Class that uses the :func:`os.urandom` function for generating random numbers
361 from sources provided by the operating system. Not available on all systems.
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000362 Does not rely on software state, and sequences are not reproducible. Accordingly,
Raymond Hettingerafd30452009-02-24 10:57:02 +0000363 the :meth:`seed` method has no effect and is ignored.
Georg Brandl116aa622007-08-15 14:28:22 +0000364 The :meth:`getstate` and :meth:`setstate` methods raise
365 :exc:`NotImplementedError` if called.
366
Georg Brandl116aa622007-08-15 14:28:22 +0000367
Raymond Hettinger435cb0f2010-09-06 23:36:31 +0000368Notes on Reproducibility
Antoine Pitroue72b5862010-12-12 20:13:31 +0000369------------------------
Raymond Hettinger435cb0f2010-09-06 23:36:31 +0000370
Julien Palard58a40542020-01-31 10:50:14 +0100371Sometimes it is useful to be able to reproduce the sequences given by a
372pseudo-random number generator. By re-using a seed value, the same sequence should be
Raymond Hettinger435cb0f2010-09-06 23:36:31 +0000373reproducible from run to run as long as multiple threads are not running.
374
375Most of the random module's algorithms and seeding functions are subject to
376change across Python versions, but two aspects are guaranteed not to change:
377
378* If a new seeding method is added, then a backward compatible seeder will be
379 offered.
380
Georg Brandl92849d12016-02-19 08:57:38 +0100381* The generator's :meth:`~Random.random` method will continue to produce the same
Raymond Hettinger435cb0f2010-09-06 23:36:31 +0000382 sequence when the compatible seeder is given the same seed.
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000383
Raymond Hettinger6e353942010-12-04 23:42:12 +0000384.. _random-examples:
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000385
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000386Examples and Recipes
Antoine Pitroue72b5862010-12-12 20:13:31 +0000387--------------------
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000388
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800389Basic examples::
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000390
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800391 >>> random() # Random float: 0.0 <= x < 1.0
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000392 0.37444887175646646
393
Raymond Hettingere1329102016-11-21 12:33:50 -0800394 >>> uniform(2.5, 10.0) # Random float: 2.5 <= x < 10.0
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800395 3.1800146073117523
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000396
Raymond Hettingere1329102016-11-21 12:33:50 -0800397 >>> expovariate(1 / 5) # Interval between arrivals averaging 5 seconds
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800398 5.148957571865031
399
Raymond Hettingere1329102016-11-21 12:33:50 -0800400 >>> randrange(10) # Integer from 0 to 9 inclusive
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000401 7
402
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800403 >>> randrange(0, 101, 2) # Even integer from 0 to 100 inclusive
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000404 26
405
Raymond Hettinger6befb642016-11-21 01:59:39 -0800406 >>> choice(['win', 'lose', 'draw']) # Single random element from a sequence
407 'draw'
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000408
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800409 >>> deck = 'ace two three four'.split()
410 >>> shuffle(deck) # Shuffle a list
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400411 >>> deck
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800412 ['four', 'two', 'ace', 'three']
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000413
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800414 >>> sample([10, 20, 30, 40, 50], k=4) # Four samples without replacement
415 [40, 10, 50, 30]
Raymond Hettinger3cdf8712010-12-02 05:35:35 +0000416
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800417Simulations::
418
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800419 >>> # Six roulette wheel spins (weighted sampling with replacement)
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400420 >>> choices(['red', 'black', 'green'], [18, 18, 2], k=6)
421 ['red', 'green', 'black', 'black', 'red', 'black']
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000422
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800423 >>> # Deal 20 cards without replacement from a deck of 52 playing cards
424 >>> # and determine the proportion of cards with a ten-value
425 >>> # (a ten, jack, queen, or king).
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800426 >>> deck = collections.Counter(tens=16, low_cards=36)
427 >>> seen = sample(list(deck.elements()), k=20)
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800428 >>> seen.count('tens') / 20
Raymond Hettinger0a1a9092016-11-17 00:45:35 -0800429 0.15
430
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800431 >>> # Estimate the probability of getting 5 or more heads from 7 spins
432 >>> # of a biased coin that settles on heads 60% of the time.
Raymond Hettinger9abb7252019-02-15 12:40:18 -0800433 >>> def trial():
434 ... return choices('HT', cum_weights=(0.60, 1.00), k=7).count('H') >= 5
435 ...
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700436 >>> sum(trial() for i in range(10_000)) / 10_000
Raymond Hettinger16ef5d42016-10-31 22:53:52 -0700437 0.4169
438
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800439 >>> # Probability of the median of 5 samples being in middle two quartiles
Raymond Hettinger9abb7252019-02-15 12:40:18 -0800440 >>> def trial():
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700441 ... return 2_500 <= sorted(choices(range(10_000), k=5))[2] < 7_500
Raymond Hettinger9abb7252019-02-15 12:40:18 -0800442 ...
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700443 >>> sum(trial() for i in range(10_000)) / 10_000
Raymond Hettinger71c62e12016-12-04 11:00:34 -0800444 0.7958
445
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400446Example of `statistical bootstrapping
447<https://en.wikipedia.org/wiki/Bootstrapping_(statistics)>`_ using resampling
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700448with replacement to estimate a confidence interval for the mean of a sample::
Raymond Hettinger2fdc7b12010-12-02 02:41:33 +0000449
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400450 # http://statistics.about.com/od/Applications/a/Example-Of-Bootstrapping.htm
Raymond Hettinger47d99872019-02-21 15:06:29 -0800451 from statistics import fmean as mean
Raymond Hettinger1c3a1212016-10-12 01:42:10 -0400452 from random import choices
Raymond Hettingerf5b7c7b2016-09-05 13:15:02 -0700453
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700454 data = [41, 50, 29, 37, 81, 30, 73, 63, 20, 35, 68, 22, 60, 31, 95]
455 means = sorted(mean(choices(data, k=len(data))) for i in range(100))
Raymond Hettinger2589ee32016-11-16 21:34:17 -0800456 print(f'The sample mean of {mean(data):.1f} has a 90% confidence '
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700457 f'interval from {means[5]:.1f} to {means[94]:.1f}')
Raymond Hettinger2589ee32016-11-16 21:34:17 -0800458
Raymond Hettinger00305ad2016-11-16 22:56:11 -0800459Example of a `resampling permutation test
460<https://en.wikipedia.org/wiki/Resampling_(statistics)#Permutation_tests>`_
461to determine the statistical significance or `p-value
462<https://en.wikipedia.org/wiki/P-value>`_ of an observed difference
463between the effects of a drug versus a placebo::
464
465 # Example from "Statistics is Easy" by Dennis Shasha and Manda Wilson
Raymond Hettinger47d99872019-02-21 15:06:29 -0800466 from statistics import fmean as mean
Raymond Hettinger00305ad2016-11-16 22:56:11 -0800467 from random import shuffle
468
469 drug = [54, 73, 53, 70, 73, 68, 52, 65, 65]
470 placebo = [54, 51, 58, 44, 55, 52, 42, 47, 58, 46]
471 observed_diff = mean(drug) - mean(placebo)
472
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700473 n = 10_000
Raymond Hettinger00305ad2016-11-16 22:56:11 -0800474 count = 0
475 combined = drug + placebo
476 for i in range(n):
477 shuffle(combined)
478 new_diff = mean(combined[:len(drug)]) - mean(combined[len(drug):])
479 count += (new_diff >= observed_diff)
480
481 print(f'{n} label reshufflings produced only {count} instances with a difference')
482 print(f'at least as extreme as the observed difference of {observed_diff:.1f}.')
483 print(f'The one-sided p-value of {count / n:.4f} leads us to reject the null')
Raymond Hettinger6befb642016-11-21 01:59:39 -0800484 print(f'hypothesis that there is no difference between the drug and the placebo.')
485
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700486Simulation of arrival times and service deliveries for a multiserver queue::
Raymond Hettinger6befb642016-11-21 01:59:39 -0800487
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700488 from heapq import heappush, heappop
Raymond Hettinger1149d932016-11-21 14:13:07 -0800489 from random import expovariate, gauss
490 from statistics import mean, median, stdev
Raymond Hettinger6befb642016-11-21 01:59:39 -0800491
492 average_arrival_interval = 5.6
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700493 average_service_time = 15.0
494 stdev_service_time = 3.5
495 num_servers = 3
Raymond Hettinger6befb642016-11-21 01:59:39 -0800496
Raymond Hettingerd3a8d612020-04-21 16:11:00 -0700497 waits = []
498 arrival_time = 0.0
499 servers = [0.0] * num_servers # time when each server becomes available
500 for i in range(100_000):
501 arrival_time += expovariate(1.0 / average_arrival_interval)
502 next_server_available = heappop(servers)
503 wait = max(0.0, next_server_available - arrival_time)
504 waits.append(wait)
505 service_duration = gauss(average_service_time, stdev_service_time)
506 service_completed = arrival_time + wait + service_duration
507 heappush(servers, service_completed)
Raymond Hettinger1149d932016-11-21 14:13:07 -0800508
Raymond Hettinger1149d932016-11-21 14:13:07 -0800509 print(f'Mean wait: {mean(waits):.1f}. Stdev wait: {stdev(waits):.1f}.')
510 print(f'Median wait: {median(waits):.1f}. Max wait: {max(waits):.1f}.')
Raymond Hettinger6befb642016-11-21 01:59:39 -0800511
Raymond Hettinger05374052016-11-21 10:52:04 -0800512.. seealso::
513
514 `Statistics for Hackers <https://www.youtube.com/watch?v=Iq9DzN6mvYA>`_
515 a video tutorial by
516 `Jake Vanderplas <https://us.pycon.org/2016/speaker/profile/295/>`_
517 on statistical analysis using just a few fundamental concepts
518 including simulation, sampling, shuffling, and cross-validation.
519
520 `Economics Simulation
521 <http://nbviewer.jupyter.org/url/norvig.com/ipython/Economics.ipynb>`_
522 a simulation of a marketplace by
523 `Peter Norvig <http://norvig.com/bio.html>`_ that shows effective
Raymond Hettinger7f946192016-11-21 15:13:18 -0800524 use of many of the tools and distributions provided by this module
Raymond Hettinger05374052016-11-21 10:52:04 -0800525 (gauss, uniform, sample, betavariate, choice, triangular, and randrange).
526
527 `A Concrete Introduction to Probability (using Python)
528 <http://nbviewer.jupyter.org/url/norvig.com/ipython/Probability.ipynb>`_
529 a tutorial by `Peter Norvig <http://norvig.com/bio.html>`_ covering
530 the basics of probability theory, how to write simulations, and
Raymond Hettinger7f946192016-11-21 15:13:18 -0800531 how to perform data analysis using Python.