Random numbers, simulations and statistics

How likely is at least one six in four rolls of a die? Python answers in a second by playing the game 100,000 times and counting, and the same tools pick quiz questions, shuffle a deck and split data for a test. They come from the random module and its pseudo-randompseudo-random: Made by a fixed algorithm that only looks random. The random module's Mersenne Twister turns a seed into such a stream, so anyone who learns the generator's state can work out every number that follows. generator: a fixed algorithm, the Mersenne Twister, turns a starting value called the seedseed: The starting value of a pseudo-random generator. The same seed gives the same stream of numbers, so random.seed(42) or random.Random(42) makes a simulation or a test repeatable. into a long stream of numbers that pass for random, and the same seed always gives the same stream.

You will make runs repeatable, pick with weights and without repeats, measure how far a simulated answer can be from the truth, and describe the numbers you collect with the statistics module. One rule comes before all of it: random is made for games, simulations and tests, never for passwords or tokens, which come from the secrets module.

Seeds make a run repeatable

random.random() returns a float from 0.0 up to, but not including, 1.0, and every other function of the module draws on the same hidden generator. Python seeds that generator from the operating system’s source of randomness when the module is imported, so two runs of a program normally see different numbers. Call random.seed() with a value of your own and the stream starts again from a known point:

import random

# The same seed gives the same stream, run after run.
random.seed(42)
first = [random.randint(1, 6) for _ in range(10)]
random.seed(42)
again = [random.randint(1, 6) for _ in range(10)]
print(first)
print("same again:", again == first)

# Any other draw from the shared generator moves the stream along.
random.seed(42)
random.random()
print([random.randint(1, 6) for _ in range(10)])

# A Random object has a generator of its own.
dice = random.Random(42)
coins = random.Random(7)
print([dice.randint(1, 6) for _ in range(10)])
print("".join(coins.choice("HT") for _ in range(10)))

[6, 1, 1, 6, 3, 2, 2, 2, 6, 1]
same again: True
[1, 6, 3, 2, 2, 2, 6, 1, 6, 6]
[6, 1, 1, 6, 3, 2, 2, 2, 6, 1]
THTHHHTHHH

• The first two lines are the whole point of a seed: ten rolls after seed(42), then exactly the same ten after a second seed(42).
• In the third line, one call to random() came first. It used up the two internal numbers behind the first two rolls of the first line, so the rolls start where its third roll was. Any code that draws from the module’s shared generator moves the stream along for everyone else.
• A random.Random object carries a generator of its own: Random(42) rolls the same dice as seed(42), and nothing outside the object can disturb its stream. Give each simulation or test its own Random instead of seeding the module.

Step through the same run as requests to a generator and the rolls it hands back, and watch where each line of rolls begins:

The stream never changed; only the point at which the rolls began to read it did. What sits behind the stream is drawn below. The state is 624 numbers of 32 bits, and each output is one of them, scrambled:

That design gives a period of 2¹⁹⁹³⁷ − 1 outputs before the stream repeats, and floats with 53 bits of precision. Seeded streams do not depend on the computer: the example above printed the same rolls on CPython 3.14.8 and on the Python that runs in your browser.

Across Python versions the promise is narrower. The documentation guarantees only that random() keeps producing the same stream for the same seed; the algorithms behind randint(), choice(), shuffle() and the rest may change in a new release. A test that compares a seeded shuffle with a stored list can therefore break after an upgrade, while a test that checks properties of the result (the same items, the right length, a share near the expected one) keeps passing. Threads are the other catch: the module’s functions share one generator, so in a program with several threads the order of draws, and with it what each thread sees, can change from run to run.

choices() can repeat an item, sample() never does

A quiz app makes a good test bench: it needs questions picked at random, the hard ones less often, and no question twice in one quiz.

import random
from collections import Counter

rng = random.Random(5)
bank = {
    "easy": ["E1", "E2", "E3", "E4"],
    "medium": ["M1", "M2", "M3"],
    "hard": ["H1", "H2"],
}
levels = list(bank)
weights = [5, 3, 2]  # out of every ten picks: five easy, three medium, two hard

# choices(): k picks with replacement, in proportion to the weights
print(rng.choices(levels, weights=weights, k=10))
counts = Counter(rng.choices(levels, weights=weights, k=10_000))
print({level: f"{counts[level] / 10_000:.1%}" for level in levels})

# sample(): k different items, never the same one twice
every_question = [q for level in levels for q in bank[level]]
print(rng.sample(every_question, k=4))

# shuffle(): reorders the list itself and returns None
order = every_question.copy()
result = rng.shuffle(order)
print(order, result)

# randint() includes both ends; randrange() stops before its end
print(rng.randint(1, 3), rng.randrange(1, 3), round(rng.uniform(0, 10), 3))

['medium', 'medium', 'medium', 'hard', 'medium', 'hard', 'easy', 'easy', 'hard', 'medium']
{'easy': '50.0%', 'medium': '30.1%', 'hard': '19.9%'}
['E3', 'M3', 'M2', 'E2']
['M3', 'E4', 'M1', 'E3', 'E2', 'E1', 'M2', 'H1', 'H2'] None
3 2 5.338

• choices(population, weights=…, k=…) makes k independent picks with replacementwith replacement: Drawing so that every pick can take any item again, as if each one were put back first. random.choices() draws with replacement; random.sample() draws without, using each position at most once., each in proportion to its weight. Ten picks can look nothing like 5 : 3 : 2 (only two of them came out easy), but 10,000 picks land within a few tenths of a percent of 50, 30 and 20 %. cum_weights takes running totals instead (5, 8, 10), which saves work when you pick from the same weights many times.
• sample(population, k) picks k different positions, so no question comes up twice. It raises ValueError when k is larger than the population, and since Python 3.11 it accepts only sequences: pass sorted(a_set), not a set. Its counts argument repeats items without building the longer list: sample(["gold", "silver"], counts=[1, 9], k=3) draws three of ten coins, so at most one of the three is gold.
• shuffle(x) reorders the list in placein place: Changing an object itself instead of making a new one. list.sort() sorts a list in place and returns None, while sorted() leaves its input as it was and returns a new list. More in “Lists: create, index, slice and modify” and returns None, the None the fourth line printed. For a shuffled copy, use sample(x, k=len(x)). For a list of more than 2,080 items the number of possible orders is larger than the generator’s period, so most orders can never come up.
• randint(1, 3) can return 3, because it includes both ends. randrange(1, 3) cannot: like range(), it stops before its end. uniform(0, 10) returns a float between the two.

Which call fits depends on three things about your draw. Answer them in the table: it rules out the calls that do not fit, with the reason from the run above.

In short: choices() for independent picks, with weights if you need them; sample() for different items; shuffle() for a list you no longer need in its old order; sample(x, k=len(x)) for a shuffled copy. choice(seq) is the simple case of one item, for one winner or one move in a game.

The module also draws from distributions: uniform(a, b) gives any value in a range, gauss(mu, sigma) a value on a bell curve and expovariate() a waiting time, and since Python 3.12 binomialvariate(n, p) gives a whole number, the successes in n tries. Python 3.13 also added a command line: python -m random 6 prints a whole number from 1 to 6, and python -m random tea coffee picks one of the words.

A simulation estimates a probability by counting

A simulation answers a probability question by repeating the experiment many times and counting how often the event happened. Take an old gambling question: roll a die four times; how likely is at least one six? Each trial below is one game of four rolls, and the program prints the estimates of four runs, of 100, 1,000, 10,000 and 100,000 games:

import math
import random


def at_least_one_six(rng, rolls=4):
    """One experiment: roll a die four times. Did a six come up?"""
    return any(rng.randint(1, 6) == 6 for _ in range(rolls))


exact = 1 - (5 / 6) ** 4
print(f"exact: {exact:.4f}")
rng = random.Random(2)
for trials in (100, 1_000, 10_000, 100_000):
    hits = sum(at_least_one_six(rng) for _ in range(trials))
    estimate = hits / trials
    typical = math.sqrt(exact * (1 - exact) / trials)
    print(f"{trials:>7,} trials: {estimate:.4f}  off by {abs(estimate - exact):.4f}  typical error {typical:.4f}")

exact: 0.5177
    100 trials: 0.5100  off by 0.0077  typical error 0.0500
  1,000 trials: 0.5320  off by 0.0143  typical error 0.0158
 10,000 trials: 0.5158  off by 0.0019  typical error 0.0050
100,000 trials: 0.5176  off by 0.0002  typical error 0.0016

The estimates close in on 0.5177, the exact answer the program prints first. It comes from counting the other way: a game has no six only if all four rolls miss, each with probability 5/6, so the chance of at least one six is 1 − (5/6)⁴. The number of sixes in four rolls follows the binomial distributionbinomial distribution: The number of successes in n independent tries that each succeed with probability p. The chance of exactly k successes is C(n, k) × pᵏ × (1 − p)ⁿ⁻ᵏ; the sixes in four rolls of a die follow it with n = 4 and p = 1/6., which the next figure draws for n tries that each succeed with probability p. It starts at the dice game: 4 tries, a chance of 0.166667 (1/6) and at least 1 success.

The highlighted bars add up to 0.5177, the answer the simulation is aiming at, and the bar at 0 holds the rest, (5/6)⁴ ≈ 0.4823. Where both exist, the exact answer is the best test of the simulation’s code. The simulation earns its place when the exact answer is hard to work out: the same few lines count any rule you can write as code, while a formula has to be found again for every new rule.

Four times the trials halve the error

The estimates in the table are never exact, but their error shrinks in a predictable way. For a probability p estimated from n independent trials, the typical error, called the standard errorstandard error: The typical distance between an estimate and the true value. For a probability p estimated from n independent trials it is √(p(1 − p) / n), so four times as many trials halve it., is:

In words: the error is largest when p is near one half and shrinks with the square root of the number of trials. For the dice game and 1,000 trials it is √(0.5177 × 0.4823 ÷ 1,000) ≈ 0.0158, the typical error the program printed in that row.

An estimate’s error behaves like a draw from a bell curve centred on 0 whose spread is the standard error, so in about two runs out of three the estimate lands within one standard error of the truth, and in about 19 out of 20 within two. The figure draws the errors of 1,000-trial runs in percentage points, with a spread of 1.58 points and a band of ±1.58 shaded:

The band holds 0.6827 of all runs; widen it to ±3.16 points, two standard errors, and it holds 0.9545. Every row of the program’s table landed within one standard error of its own size: the 1,000-trial run, for one, was off by 1.43 points, inside the band of ±1.58 drawn here.

Because n sits under a square root, halving the error takes four times as many trials. Seeding the Random object makes the whole table repeatable: run it again and every estimate is the same. Another seed gives other estimates, spread around the exact answer by about the same typical error.

statistics finds the middle, the spread and the quartiles

The statistics module describes a list of numbers: where its middle is, how spread out it is and how it divides into quarters. Here are the marks of a class of twenty:

import statistics as st

marks = [72, 85, 64, 91, 58, 77, 85, 69, 94, 81, 66, 85, 73, 88, 47, 79, 70, 83, 61, 90]

print("students ", len(marks))
print("mean     ", st.mean(marks))
print("median   ", st.median(marks))
print("mode     ", st.mode(marks))
print("stdev    ", round(st.stdev(marks), 2), "  pstdev", round(st.pstdev(marks), 2))
print("quartiles", st.quantiles(marks, n=4))
print("inclusive", st.quantiles(marks, n=4, method="inclusive"))
print("deciles  ", st.quantiles(marks, n=10))

grades = ["B", "A", "C", "A", "B"]
print("tied     ", st.mode(grades), st.multimode(grades))

# Two tests taken by the same eight students
mid = [55, 62, 70, 48, 81, 66, 74, 90]
final = [58, 70, 72, 52, 85, 63, 80, 94]
print("correlation", round(st.correlation(mid, final), 3))
fit = st.linear_regression(mid, final)
print(f"final = {fit.slope:.3f} × mid + {fit.intercept:.2f}")

for summary, data in ((st.mean, []), (st.stdev, [42])):
    try:
        summary(data)
    except st.StatisticsError as error:
        print("StatisticsError:", error)

students  20
mean      75.9
median    78.0
mode      85
stdev     12.48   pstdev 12.17
quartiles [66.75, 78.0, 85.0]
inclusive [68.25, 78.0, 85.0]
deciles   [58.3, 64.4, 69.3, 72.4, 78.0, 82.2, 85.0, 87.4, 90.9]
tied      B ['B', 'A']
correlation 0.974
final = 1.009 × mid + 2.87
StatisticsError: mean requires at least one data point
StatisticsError: stdev requires at least two data points

The figure draws the same twenty marks as dots on one line, with the box between the first and third quartilesquartile: One of the three cut points that split sorted data into four groups of the same size; the middle one is the median. statistics.quantiles(data, n=4) returns all three. that the program printed, the median across it and the mean below:

• mean() is the average and median() the middle value, here 75.9 and 78.0. One low mark, 47, pulls the mean down more than the median, which is why the mean’s diamond sits left of the median’s line, and why medians describe skewed data such as incomes better.
• mode() is the most common value. With a tie it returns the first of the tied values in the data, so "B" wins here; multimode() returns all of them.
• stdev() is the sample standard deviation. Its sum of squares is divided by n − 1 rather than n, a correction for the group being only part of everyone you care about; pstdev() divides by n and describes exactly the data you have. Use pstdev() when the class is all you want to describe, and stdev() when it stands for every student who sat the test.
• quantiles(data, n=4) returns the three cut points that split the data into quarters; n=10 gives deciles. The default method, "exclusive", assumes the real population can go beyond your lowest and highest values; "inclusive" treats them as the true minimum and maximum, so its outer quartiles are never further out: here the lower one moves from 66.75 to 68.25.
• correlation() measures how closely two lists follow a straight line together, from −1 to 1, and linear_regression() fits that line: here a student’s final mark is about 1.009 times their mid-term mark plus 2.87.
• With too few values there is nothing to compute: mean([]) and stdev([42]) raise StatisticsError, a kind of ValueError.

The functions accept ints, floats, Decimal and Fraction values, and they keep exact types exact: the mean of Fraction values is a Fraction. Mixing types in one list is not supported, so convert first, for example with map(float, data). The module aims at the level of a graphing or scientific calculator and leaves large-scale work to libraries such as NumPy and SciPy.

random is not for secrets

The documentation says it plainly: the generators of random should not be used for security. The stream is predictable, and not only to someone who knows the seed. Each output is a state number put through a scrambling step that can be undone, so 624 outputs in a row hand over the whole state. This program watches a generator seeded by the operating system, whose seed nobody knows, and predicts it; then it runs the same check on the generator that the secrets module uses:

import random
import secrets
import string


def undo_right(value, shift):
    """Undo value ^= value >> shift on a 32-bit number."""
    result = value
    for _ in range(32 // shift):
        result = value ^ (result >> shift)
    return result


def undo_left(value, shift, mask):
    """Undo value ^= (value << shift) & mask on a 32-bit number."""
    result = value
    for _ in range(32 // shift):
        result = value ^ ((result << shift) & mask)
    return result & 0xFFFFFFFF


def unscramble(output):
    """Turn one 32-bit output back into the state number it was made from."""
    y = undo_right(output, 18)
    y = undo_left(y, 15, 0xEFC60000)
    y = undo_left(y, 7, 0x9D2C5680)
    return undo_right(y, 11)


def predicts(generator):
    """Watch 624 outputs, copy the state they reveal, and check the copy's next five numbers."""
    seen = [generator.getrandbits(32) for _ in range(624)]
    copy = random.Random()
    copy.setstate((3, tuple(unscramble(n) for n in seen) + (624,), None))
    upcoming = [generator.getrandbits(32) for _ in range(5)]
    return upcoming == [copy.getrandbits(32) for _ in range(5)], copy


server = random.Random()  # seeded by the operating system: nobody knows its seed
hit, copy = predicts(server)
print("next 5 numbers predicted:", hit)

alphabet = string.ascii_letters + string.digits
token = "".join(server.choice(alphabet) for _ in range(20))
guess = "".join(copy.choice(alphabet) for _ in range(20))
print("20-character token guessed:", token == guess)

# The same check on the generator behind the secrets module: its outputs reveal nothing
hit, _ = predicts(secrets.SystemRandom())
print("secrets.SystemRandom predicted:", hit)

next 5 numbers predicted: True
20-character token guessed: True
secrets.SystemRandom predicted: False

After seeing 624 numbers, the copy predicts the next five exactly, and with them a 20-character “reset token” built with choice(). The same attempt on secrets.SystemRandom fails, because its numbers come from the operating system’s most secure source and have no state to copy. The difference in a table:

For anything that guards an account or money (passwords, reset links, session ids, one-time codes, API keys) use the secrets module:

import secrets
import string

# The values change on every run, so this prints only their shape.
token = secrets.token_urlsafe(32)  # 32 random bytes, Base64-encoded for URLs
print(len(token), set(token) <= set(string.ascii_letters + string.digits + "-_"))

alphabet = string.ascii_letters + string.digits
password = "".join(secrets.choice(alphabet) for _ in range(16))
print(len(password), password.isalnum())

code = f"{secrets.randbelow(1_000_000):06d}"  # a six-digit one-time code
print(len(code), code.isdigit())

# Compare a submitted token with the stored one in constant time
print(secrets.compare_digest(token, token), secrets.compare_digest(token, token[:-1] + "?"))

43 True
16 True
6 True
True False

• secrets.token_urlsafe(32) gives 32 random bytes as 43 characters that are safe in a URL. The documentation names 32 bytes, 256 bits, as enough for typical uses; token_hex() and token_bytes() work the same way.
• secrets.choice() picks a character, and secrets.randbelow(n) a whole number from 0 to n − 1, here padded to a six-digit code.
• secrets.compare_digest() compares two secrets in constant time, so the time an answer takes does not reveal how many leading characters were right.

The program prints only lengths and checks because the tokens change on every run, which is exactly what they should do.

Practise: shared birthdays and a marks summary

How many people must be in a room before two of them probably share a birthday? Fewer than most people guess: 23 is enough for a chance just above one half. Simulate it with a generator of your own, then summarise a list of marks; the sample tests allow for the random error of a simulation.

import random
import statistics


def simulate_birthday(n, trials, seed):
    """Estimate the chance that at least two of n people share a birthday (365 equally likely days)."""
    # Replace this line with your code.
    return 0.0


def summarise(marks):
    """Return the count, mean, median, sample stdev, lowest, highest and quartiles of the marks."""
    # Replace this line with your code.
    return {}

import math

from simulate import simulate_birthday, summarise


def test_twenty_three_people():
    """estimates about a 50.7 % chance for 23 people"""
    estimate = simulate_birthday(23, 4000, seed=1)
    assert math.isclose(estimate, 0.5073, abs_tol=0.03)


def test_other_group_sizes():
    """gets close to the exact answers for 10 and 50 people"""
    assert math.isclose(simulate_birthday(10, 4000, seed=4), 0.1169, abs_tol=0.03)
    assert math.isclose(simulate_birthday(50, 2000, seed=3), 0.9704, abs_tol=0.03)


def test_certain_cases():
    """gives 0.0 for one person and 1.0 for 366 people"""
    assert simulate_birthday(1, 100, seed=1) == 0.0
    assert simulate_birthday(366, 10, seed=1) == 1.0


def test_repeatable():
    """returns the same estimate for the same seed"""
    assert simulate_birthday(30, 500, seed=7) == simulate_birthday(30, 500, seed=7)


def test_summary():
    """describes ten marks"""
    summary = summarise([72, 85, 64, 91, 58, 77, 85, 69, 94, 81])
    assert summary["count"] == 10
    assert math.isclose(summary["mean"], 77.6)
    assert summary["median"] == 79.0
    assert math.isclose(summary["stdev"], 11.7587, abs_tol=0.0001)
    assert summary["lowest"] == 58
    assert summary["highest"] == 94
    assert summary["quartiles"] == [67.75, 79.0, 86.5]


def test_too_few_marks():
    """raises a ValueError for a single mark"""
    try:
        summarise([50])
    except ValueError:
        return
    assert False, "expected a ValueError for one mark"

Eighteen questions on the six ideas of this lesson follow. A session opens with the idea your answers show you know least, moves up a level after each right answer, and says why every choice is right or wrong the moment you pick it:

Interview questions

Fresher · What is the difference between random.sample() and random.choices()?

sample(population, k) draws without replacement: it returns k items from k different positions, so an item that appears once in the population appears at most once in the result, and k may not be larger than the population. choices(population, weights=None, k=1) draws with replacement: each of the k picks is independent, so items can repeat, k can be any size, and weights or cumulative weights make some items more likely than others. Use sample() for a raffle, a quiz without repeated questions or splitting data into a test set; use choices() for simulations in which the same outcome can happen again, such as customers choosing products. For a shuffled copy of a whole list, sample(x, k=len(x)) does the job.

Key takeaways

• random produces a repeatable stream from a seed, shared by every function of the module; give each simulation or test its own random.Random(seed), and test seeded code by its properties, because only random() is promised to keep its stream across versions.
• choices() picks with replacement and weights, sample() without repeats, and shuffle() reorders a list in place and returns None; randint() includes both ends and randrange() stops before its end.
• A simulation estimates a probability as the share of trials in which the event happened, and an exact answer such as 1 − (5/6)⁴ is the best check of its code.
• A simulation’s typical error is √(p(1 − p) / n): ten times the trials shrink it about 3.2 times, and halving it costs four times the trials.
• statistics gives the mean, median, mode, sample and population standard deviations, quantiles, correlation and a regression line, and raises StatisticsError when there are too few values.
• Never use random for secrets, because 624 outputs reveal its state: use secrets.token_urlsafe(), secrets.choice(), secrets.randbelow() and secrets.compare_digest().