Python random Module and Monte Carlo Simulation Examples

By Engr. Henry Akpan

Let the Computer Roll the Dice

What is the chance that two people in a class of 30 share a birthday? How likely is it that a student who guesses every answer passes a 40-question multiple-choice test? Some probability questions are hard to work out with formulas, but easy to answer by simulation: let the computer play out the situation thousands of times and count what happens. This approach is called the Monte Carlo method, named after the famous casino, and it is used in physics, finance, engineering, weather forecasting and artificial intelligence. This tutorial covers Python's random module and Monte Carlo simulation: random numbers, integers and choices, seeds for reproducible results, sampling, shuffling and weighted choices, simulating dice, estimating π, the birthday problem, normally distributed data, and how accuracy improves with more trials. All outputs were produced by running the code on Python 3.13; because the examples use fixed seeds, you will get the same numbers.

Random Numbers in Python

The random module provides pseudo-random numbers: they are produced by a deterministic algorithm (the Mersenne Twister) but look random for all practical purposes. The main functions are:

  • random.random(): a float from 0 up to (but not including) 1.
  • random.randint(a, b): an integer from a to b, including both ends.
  • random.uniform(a, b): a float between a and b.
import random

random.seed(42)
print(random.random())
print(random.randint(1, 6))
print(random.uniform(10, 20))

Output:

0.6394267984578837
1
17.41550499759833

Seeds: Making Randomness Repeatable

A seed sets the starting point of the generator. The same seed always produces the same sequence, which is invaluable for testing, debugging and letting other people reproduce your results:

random.seed(42)
first = [random.randint(1, 6) for _ in range(5)]
random.seed(42)
second = [random.randint(1, 6) for _ in range(5)]
print(first, second, first == second)

Output:

[6, 1, 1, 6, 3] [6, 1, 1, 6, 3] True

Without random.seed(), Python seeds the generator from the operating system, so every run is different. In a game you want that; in a scientific simulation you usually record the seed.

Security warning: the random module is predictable, so never use it for passwords, PINs or tokens. Use the secrets module instead, for example secrets.token_hex(16).

Choosing, Sampling and Shuffling

random.seed(7)
students = ["Ada", "Bayo", "Chidi", "Dami", "Efe", "Fola"]

print(random.choice(students))           # one item
print(random.sample(students, 3))        # 3 different items
print(random.choices(["H", "T"], k=10))  # 10 picks WITH replacement

deck = list(range(1, 11))
random.shuffle(deck)                     # shuffles in place
print(deck)

colours = random.choices(["red", "green", "blue"], weights=[5, 3, 2], k=8)
print(colours)

Output:

Chidi
['Bayo', 'Dami', 'Ada']
['H', 'T', 'H', 'H', 'T', 'H', 'H', 'H', 'H', 'H']
[8, 9, 3, 1, 5, 7, 6, 4, 2, 10]
['red', 'red', 'green', 'red', 'red', 'green', 'green', 'green']
  • sample() picks without replacement, so no item is chosen twice: ideal for selecting a committee or exam questions.
  • choices() picks with replacement, and its weights make some outcomes more likely. Here red has a 50% chance.
  • shuffle() changes the list itself and returns None.

Notice that 8 of the 10 coin tosses were heads. Short runs of random results are often lopsided, which is exactly why simulations need many trials.

The Monte Carlo Method

Every Monte Carlo simulation follows the same pattern:

  1. Write a function that simulates one trial of the random situation.
  2. Run it a large number of times.
  3. Count how often the outcome of interest happens, or average the results.
  4. The proportion estimates the probability, and it gets more accurate as the number of trials grows (the law of large numbers).

Example 1: The Sum of Two Dice

from collections import Counter

random.seed(1)
rolls = [random.randint(1, 6) + random.randint(1, 6) for _ in range(10_000)]
counts = Counter(rolls)

for total in range(2, 13):
    share = counts[total] / len(rolls)
    print(f"{total:2} {share:.3f} {'#' * (counts[total] // 100)}")

Output:

 2 0.027 ##
 3 0.053 #####
 4 0.087 ########
 5 0.114 ###########
 6 0.137 #############
 7 0.168 ################
 8 0.145 ##############
 9 0.108 ##########
10 0.080 ########
11 0.056 #####
12 0.026 ##

The text histogram shows the triangle shape you expect. The exact probability of a 7 is 6/36 ≈ 0.167, and the simulation gave 0.168. The tutorial on probability explains how to get the exact values.

Example 2: Estimating π

Throw random points into a 1 × 1 square. The quarter-circle of radius 1 covers π/4 of the square's area, so the fraction of points that land inside it (where x² + y² ≤ 1) should be close to π/4. Multiply by 4 to estimate π:

def estimate_pi(n):
    inside = 0
    for _ in range(n):
        x, y = random.random(), random.random()
        if x * x + y * y <= 1:
            inside += 1
    return 4 * inside / n

random.seed(2024)
for n in (100, 10_000, 1_000_000):
    print(n, estimate_pi(n))

Output:

100 3.32
10000 3.1152
1000000 3.142668

With 100 points the estimate is rough; with a million it is correct to about 3 significant figures (π = 3.14159…). The error of a Monte Carlo estimate shrinks roughly in proportion to 1/√n: to get one more correct digit you need about 100 times as many trials. Estimating pi by simulation is not the fastest way to compute π, but the same idea estimates areas and volumes that have no simple formula.

Example 3: The Birthday Problem

How many people must be in a room for there to be a better-than-even chance that two share a birthday? Most people guess around 180. Let us simulate it (ignoring leap years):

def shared_birthday(people):
    days = [random.randint(1, 365) for _ in range(people)]
    return len(set(days)) < people      # a repeat shrinks the set

random.seed(3)
trials = 20_000
for people in (10, 23, 30, 50):
    hits = sum(shared_birthday(people) for _ in range(trials))
    print(people, round(hits / trials, 3))

Output:

10 0.118
23 0.506
30 0.704
50 0.972

With just 23 people, the chance is already about 50%, and in a class of 30 it is about 70%. The exact values (0.117, 0.507, 0.706 and 0.970) agree closely. The trick in the code is that a set removes duplicates, so if the set is smaller than the list, at least two birthdays matched.

Example 4: Guessing in a Multiple-Choice Exam

A student guesses every answer in a 40-question test with four options per question. What score should they expect, and how likely are they to reach 20 (50%)?

import statistics

def guess_exam(questions=40, p_correct=0.25):
    return sum(random.random() < p_correct for _ in range(questions))

random.seed(11)
scores = [guess_exam() for _ in range(10_000)]
print("Average score:", round(statistics.mean(scores), 2))
print("Best score:", max(scores))
print("Share scoring 20+:", sum(s >= 20 for s in scores) / len(scores))

Output:

Average score: 10.0
Best score: 23
Share scoring 20+: 0.0009

Guessing gives an average of 10 out of 40, and only about 1 in 1,000 guessers reaches 20. (The exact binomial probability is about 0.0006; with such a rare event, 10,000 trials gives only a rough estimate.) The tutorial on the binomial distribution shows how to compute it exactly.

Example 5: Normally Distributed Data

random.gauss(mu, sigma) generates values from a normal distribution with mean mu and standard deviation sigma. This is useful for creating realistic test data:

random.seed(5)
heights = [random.gauss(170, 8) for _ in range(5)]
print([round(h, 1) for h in heights])

many = [random.gauss(170, 8) for _ in range(10_000)]
print(round(statistics.mean(many), 1), round(statistics.stdev(many), 1))

Output:

[160.6, 160.8, 175.4, 151.6, 168.9]
169.9 8.0

The sample of 10,000 values has a mean and standard deviation very close to the requested 170 and 8. For large simulations, the NumPy library generates millions of random values far faster, as shown in the tutorial on NumPy.

Tips for Good Simulations

  • Test the one-trial function on its own before running thousands of trials.
  • Use enough trials. Start with about 10,000 and increase if the answer is still changing.
  • Fix the seed while developing and for reports, so results are reproducible.
  • Compare with an exact answer when one is available, as a check that the simulation is correct.
  • Remember the 1/√n rule: rare events and high precision need many more trials.

Python random Module: Practice Exercises

  1. Simulate 1,000 tosses of a fair coin and print the proportion of heads.
  2. Write a function that deals a random 5-card hand from a 52-card deck using sample().
  3. Estimate the probability of rolling at least one 6 in four rolls of a die. (The exact answer is 1 − (5/6)⁴ ≈ 0.518.)
  4. Simulate a random walk: start at 0 and take 100 steps of +1 or −1. Repeat 10,000 times and find the average final distance from 0.
  5. Use the π method to estimate the area under y = x² from 0 to 1. (The exact answer is 1/3.)

Conclusion

Python's random module provides random(), randint(), uniform(), choice(), sample(), choices(), shuffle() and gauss(), and seed() makes results repeatable. The Monte Carlo method answers probability questions by simulating one trial, repeating it many times and counting the outcomes. It confirmed that 7 is the most likely total of two dice, estimated π, solved the birthday problem and showed how hopeless guessing is in an exam. More trials give better answers, with the error shrinking like 1/√n.

Related Tutorials