Writing a Simple Probability Simulator in Python

What the script will do

A probability simulator needs only three ingredients: a way to generate random events, a function that plays one trial, and a loop that repeats the trial and counts successes. Python’s standard library provides everything required, so there is nothing to install. This outline builds a small tool in four pieces and checks each piece against an answer you can work out by hand.

The code is a teaching sketch. It favors clarity over speed, and you should read it, run it and change it.

Piece 1: a generic estimator

The core is a function that takes any one-trial function and repeats it.

import random

def estimate(trial, n=10000, seed=None):
    """Run trial() n times and return the fraction that returned True."""
    if seed is not None:
        random.seed(seed)
    successes = 0
    for _ in range(n):
        if trial():
            successes += 1
    return successes / n

Passing a function as an argument keeps the loop separate from the rules of the game. You will reuse estimate for every question below. Setting the seed makes a run repeatable, which helps when debugging; the generator is a pseudo-random one, so the same seed gives the same sequence.

Piece 2: dice rolls

Define one trial as “roll four dice and see if any shows a six.”

def roll_die():
    return random.randint(1, 6)

def at_least_one_six_in_four():
    return any(roll_die() == 6 for _ in range(4))

print(estimate(at_least_one_six_in_four, n=100000, seed=1))

random.randint(1, 6) includes both ends, so it returns 1 through 6 with equal chance. The exact answer is 1 - (5/6)^4 = 671/1296, about 0.518. Your output should be near that, for example 0.51 to 0.53 for 100,000 trials, though the exact digits depend on the seed. If it is far off, there is a bug.

To explore sums instead, define a trial as roll_die() + roll_die() == 7. The exact probability is 6/36, about 0.167.

Piece 3: card draws

For cards, build a deck once and draw from it without putting cards back.

suits = ["hearts", "diamonds", "clubs", "spades"]
ranks = list(range(1, 14))
deck = [(rank, suit) for suit in suits for rank in ranks]

def two_hearts():
    hand = random.sample(deck, 2)
    return all(card[1] == "hearts" for card in hand)

print(estimate(two_hearts, n=100000, seed=2))

random.sample(deck, 2) picks two different cards, which models drawing without replacement. The deck has 52 cards, 13 of them hearts. The exact probability of two hearts is (13/52) × (12/51) = 156/2652, which simplifies to 1/17, about 0.059. If you used random.choice twice by mistake, the second draw could repeat the first card and the estimate would drift toward (1/4) × (1/4) = 0.0625. Seeing that gap is a nice way to understand why each card changes the odds.

Piece 4: adding the uncertainty

An estimate without an error range is incomplete. Add a helper that computes the approximate standard error.

import math

def estimate_with_error(trial, n=10000, seed=None):
    p = estimate(trial, n, seed)
    se = math.sqrt(p * (1 - p) / n)
    return p, 1.96 * se   # estimate and approximate 95% margin

p, margin = estimate_with_error(two_hearts, n=100000, seed=3)
print(f"{p:.4f} +/- {margin:.4f}")

For p near 0.059 and 100,000 trials, the margin is about 1.96 × the square root of (0.059 × 0.941 / 100,000), which is roughly 0.0015. So you would expect results between about 0.057 and 0.060. The relation between trials and precision is the same square-root pattern that applies to every simulation.

A multi-step game: waiting for a six

The same structure handles games with variable length. Return a number instead of True or False and average it.

def turns_until_six():
    turns = 0
    while True:
        turns += 1
        if roll_die() == 6:
            return turns

random.seed(4)
lengths = [turns_until_six() for _ in range(100000)]
print(sum(lengths) / len(lengths))

The average should come out close to 6, since the expected number of rolls until the first six is 1/p = 6. Storing all the lengths, not just the mean, lets you also compute the median or the longest game.

Habits that prevent mistakes

• Check against a known answer first. Every example above has an exact answer you can verify before moving to a question that does not.
• Keep the trial function self-contained. Each call should start fresh; do not carry state over between trials.
• Use enough trials. Fifty is far too few; start with tens of thousands for common events and more for rare ones.
• Print the seed. When you report a result, record the seed so the run can be reproduced.
• Do not overread the output. A simulation estimates the model you coded. It does not prove anything about real-world gambling, where games with a house edge lose money over time and past results do not change future odds.

Where to go next

Once these pieces work, try modeling a simple board track, comparing two stopping rules in a dice game, or counting how often a particular five-card hand appears. Each is a variation on the same pattern: write one trial, repeat it, count, and attach a margin of error. The skill of building small, checkable pieces matters more than any particular library.