Random Numbers & Simulation

Random numbers sound like the opposite of careful analysis, but data scientists use them every day: to pick a fair sample, to split data for machine learning, to create test data, and to simulate situations that are hard to calculate by hand.

In this lesson you will use NumPy's modern random number generator, learn why a seed makes your results repeatable, and answer three questions by simulation, from rolling dice to planning stock for a shop.

The modern way: default_rng

Create a generator object with np.random.default_rng(), then call its methods. Passing a number, called the seed, makes the "random" sequence the same every time you run the program.

import numpy as np

rng = np.random.default_rng(42)

print(rng.integers(1, 7, size=10))          # dice rolls: 1 to 6 (7 is excluded)
print(rng.random(3))                        # decimals from 0 up to 1
print(rng.uniform(25, 35, size=3).round(1)) # decimals between 25 and 35
print(rng.normal(loc=60, scale=12, size=5).round(1))  # bell-curve marks
print(rng.choice(["Rohtak", "Pune", "Kochi"], size=4))

[1 5 4 3 3 6 1 5 2 1]
[0.97562235 0.7611397  0.78606431]
[26.3 29.5 28.7]
[69.3 60.8 73.5 65.6 49.7]
['Rohtak' 'Pune' 'Kochi' 'Rohtak']

Run it twice: you get exactly the same numbers, because the seed is fixed. Remove 42 and the numbers change on every run. These numbers are pseudo-random: a formula produces them, and the seed is its starting point. That is good enough for statistics and simulation (but not for passwords; Python's secrets module is for that).

Choosing with probabilities and shuffling

rng.choice() accepts a p argument with one probability per item, adding up to 1. This lets you simulate realistic data, such as payment methods at an online shop.

import numpy as np

rng = np.random.default_rng(7)

modes = ["UPI", "Card", "Cash on Delivery", "Wallet"]
payments = rng.choice(modes, size=1000, p=[0.55, 0.2, 0.15, 0.1])

values, counts = np.unique(payments, return_counts=True)
for v, c in zip(values, counts):
    print(f"{v:<17}{c:>5}")

students = np.array(["Aarav", "Priya", "Kabir", "Meera", "Rohan", "Isha"])
print(rng.permutation(students))     # random order for viva exams

Card               195
Cash on Delivery   158
UPI                552
Wallet              95
['Kabir' 'Isha' 'Rohan' 'Meera' 'Priya' 'Aarav']

With 1,000 simulated payments, the counts come close to the probabilities you set (about 550 UPI, 200 card and so on), but not exactly. That small wobble is random variation, and it shrinks as you simulate more.

Simulation: estimating a probability

Monte Carlo simulation means answering a question by repeating a random experiment thousands of times and counting how often something happens. It is named after the famous casino. Here is a classic: when you roll two dice, how often is the total 7?

import numpy as np

rng = np.random.default_rng(2025)
n = 100_000

die1 = rng.integers(1, 7, size=n)
die2 = rng.integers(1, 7, size=n)
total = die1 + die2

print("P(total = 7)  simulated:", (total == 7).mean())
print("P(total = 7)  exact:    ", round(6 / 36, 5))
print("P(total >= 10) simulated:", (total >= 10).mean())

P(total = 7)  simulated: 0.16606
P(total = 7)  exact:     0.16667
P(total >= 10) simulated: 0.16308

There are no loops: 100,000 games are two arrays and one addition. (total == 7).mean() works because the mean of a boolean array is the fraction of True values. The simulated answer is within a fraction of a percent of the exact 6/36.

Simulation: the birthday problem

Simulation really helps when the maths is awkward. In a class of 40 students, what is the chance that at least two share a birthday? Most people guess "low". Let us test it.

import numpy as np

rng = np.random.default_rng(1)
classes = 20_000
size = 40

birthdays = rng.integers(1, 366, size=(classes, size))   # 20,000 classes, 40 students each
sorted_days = np.sort(birthdays, axis=1)
has_shared = (np.diff(sorted_days, axis=1) == 0).any(axis=1)

print("Classes simulated:", classes)
print("Chance of a shared birthday:", has_shared.mean().round(3))

Classes simulated: 20000
Chance of a shared birthday: 0.889

Each row is one class. Sorting each row puts equal birthdays next to each other, np.diff() subtracts neighbours, and a difference of 0 means a shared birthday. .any(axis=1) asks "does this class have at least one?". The answer is close to 9 in 10, far higher than most people guess. (The simulation ignores leap years and assumes every day is equally likely.)

A business simulation: how much stock?

A sweet shop in Lucknow sells on average 30 boxes of kaju katli a day. Daily demand varies, and a common model for "number of customers or orders in a period" is the Poisson distribution. How many boxes should the shop prepare so it runs out on no more than 1 day in 20?

import numpy as np

rng = np.random.default_rng(10)
demand = rng.poisson(lam=30, size=100_000)     # 100,000 simulated days

print("Average demand:", demand.mean().round(2))
for stock in [30, 35, 40]:
    sold_out = (demand > stock).mean()
    print(f"Stock {stock}: sells out on {sold_out:.1%} of days")

print("Stock needed for 95% of days:", np.percentile(demand, 95))

Average demand: 29.99
Stock 30: sells out on 45.0% of days
Stock 35: sells out on 15.8% of days
Stock 40: sells out on 3.2% of days
Stock needed for 95% of days: 39.0

Preparing exactly the average (30) leaves the shop short on almost half of all days, and even 40 boxes run out on about 3 days in 100. The 95th percentile, the value that 95% of simulated days stay at or below, tells the owner how many boxes cover 19 days out of 20. You will meet the Poisson distribution properly in Probability Distributions.

Common mistakes

• Expecting integers(1, 6) to include 6. The upper limit is excluded. Use integers(1, 7) for a die.
• Creating a new generator inside a loop with the same seed. Every loop gives the same numbers. Create one generator at the top and reuse it.
• Trusting a tiny simulation. 100 runs wobble a lot; use thousands or more and check that the answer settles.
• Mixing old and new styles. np.random.seed() does not affect a default_rng() generator.

Key points

• Create a generator with rng = np.random.default_rng(seed) and call methods such as integers, normal, poisson and choice.
• A seed makes pseudo-random results repeatable, which is essential for reports and experiments.
• The mean of a boolean array is the fraction of True values, a handy way to estimate probabilities.
• Monte Carlo simulation repeats a random experiment many times to estimate answers that are hard to calculate.
• The older np.random.seed() and np.random.randint() style still works but is not recommended for new code.

import numpy as np

rng = np.random.default_rng(5)
balls = rng.random(size=(100_000, 6)) < 0.30     # True = boundary
boundaries = balls.sum(axis=1)
print((boundaries >= 2).mean().round(3))

0.581