Monte Carlo Simulation in Python (With NumPy Examples)

A Monte Carlo simulation in Python answers “what’s likely to happen?” by running a random experiment thousands of times and looking at the results. NumPy makes it fast because you can generate millions of random numbers at once:

import numpy as np
rng = np.random.default_rng(2026)
x, y = rng.random(1_000_000), rng.random(1_000_000)
pi_estimate = 4 * np.mean(x**2 + y**2 <= 1)

Every number here comes from a seeded run of Python 3.12.5 with NumPy 2.5.3, so you’ll get the same results if you copy the code.

We’ll estimate pi, check simulated dice odds against the exact math, model 30 years of retirement savings in dollars, forecast a project deadline, and plot the spread of outcomes.

How Monte Carlo simulation works

1. Describe the uncertain inputs with probability distributions (uniform, normal, triangular and so on).
2. Draw random values for those inputs and compute the outcome.
3. Repeat thousands of times.
4. Summarize the outcomes with averages, percentiles and probabilities.

The more runs you do, the closer the answers get to the truth. The error shrinks roughly with the square root of the number of runs, so 100 times more runs gives about 10 times more accuracy.

Estimating pi with a Monte Carlo simulation in Python

Drop random points into a 1 by 1 square. The share that land inside the quarter circle approaches pi/4, so multiplying by 4 estimates pi:

import numpy as np

rng = np.random.default_rng(2026)

for n in (1_000, 100_000, 10_000_000):
    x, y = rng.random(n), rng.random(n)            # random points in a 1 x 1 square
    inside = (x * x + y * y <= 1).sum()            # how many land inside the quarter circle
    estimate = 4 * inside / n
    print(f"{n:>12,} points: pi = {estimate:.5f}  (error {abs(estimate - np.pi):.5f})")

Output:

       1,000 points: pi = 3.17200  (error 0.03041)
     100,000 points: pi = 3.15316  (error 0.01157)
  10,000,000 points: pi = 3.14182  (error 0.00023)

I always seed the generator with default_rng(seed) while developing. It makes runs repeatable, which matters when you’re comparing changes.

Checking a simulation against exact probability

Dice are a good sanity check because you know the exact answers. A million rolls of two dice should land very close to 1 in 6 for each of these:

import numpy as np

rng = np.random.default_rng(7)
n = 1_000_000
rolls = rng.integers(1, 7, size=(n, 2))            # two six-sided dice, a million times
totals = rolls.sum(axis=1)

print(f"P(total = 7)  simulated {np.mean(totals == 7):.4f} | exact {6/36:.4f}")
print(f"P(doubles)    simulated {np.mean(rolls[:, 0] == rolls[:, 1]):.4f} | exact {6/36:.4f}")
print(f"P(total >= 10) simulated {np.mean(totals >= 10):.4f} | exact {6/36:.4f}")

Output:

P(total = 7)  simulated 0.1665 | exact 0.1667
P(doubles)    simulated 0.1666 | exact 0.1667
P(total >= 10) simulated 0.1666 | exact 0.1667

rng.integers(1, 7) excludes the upper bound, so it gives 1 to 6. For other ways to make random numbers, see NumPy random in Python.

Monte Carlo simulation for retirement savings

Here’s where I think Monte Carlo really earns its keep. Start with $50,000, add $6,000 a year, and let yearly returns vary around 7% with 15% volatility. One average return would hide the risk:

import numpy as np

rng = np.random.default_rng(42)
sims, years = 10_000, 30
start, yearly_add = 50_000.00, 6_000.00
returns = rng.normal(loc=0.07, scale=0.15, size=(sims, years))   # 7% average, 15% volatility

balance = np.full(sims, start)
for year in range(years):
    balance = balance * (1 + returns[:, year]) + yearly_add

p10, p50, p90 = np.percentile(balance, [10, 50, 90])
print(f"after {years} years ({sims:,} simulations):")
print(f"  pessimistic (10th pct): ${p10:>13,.2f}")
print(f"  median      (50th pct): ${p50:>13,.2f}")
print(f"  optimistic  (90th pct): ${p90:>13,.2f}")
print(f"  chance of reaching $1M: {np.mean(balance >= 1_000_000):.1%}")

Output:

after 30 years (10,000 simulations):
  pessimistic (10th pct): $   366,542.43
  median      (50th pct): $   769,643.82
  optimistic  (90th pct): $ 1,722,327.50
  chance of reaching $1M: 34.0%

The returns here are illustrative, not investment advice. The point is the method: model the uncertainty, then read the percentiles instead of a single number.

How do you forecast a project deadline with Monte Carlo?

Ask a team for optimistic, likely and pessimistic estimates per task, and model each with a triangular distribution. Adding the “most likely” numbers usually underestimates the total:

import numpy as np

rng = np.random.default_rng(11)
sims = 50_000
# each task: (optimistic, most likely, pessimistic) days -> triangular distribution
tasks = {"design": (3, 5, 9), "build": (8, 12, 20), "test": (2, 4, 8), "launch": (1, 2, 4)}

total = sum(rng.triangular(lo, mode, hi, sims) for lo, mode, hi in tasks.values())
print(f"plan using 'most likely' days: {sum(m for _, m, _ in tasks.values())} days")
print(f"average simulated:             {total.mean():.1f} days")
print(f"chance of finishing in 25 days: {np.mean(total <= 25):.1%}")
print(f"85% confident by:              {np.percentile(total, 85):.1f} days")

Output:

plan using 'most likely' days: 23 days
average simulated:             26.0 days
chance of finishing in 25 days: 39.3%
85% confident by:              29.4 days

I find this one convincing for managers: the simulation shows why plans slip, since the risks mostly push one way. Read more about sum() in Python, which adds the four NumPy arrays here.

Plotting Monte Carlo results

A histogram shows the whole spread of outcomes at a glance, including the long tail of very good and very bad results:

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
returns = rng.normal(0.07, 0.15, size=(10_000, 30))
balance = np.full(10_000, 50_000.0)
for year in range(30):
    balance = balance * (1 + returns[:, year]) + 6_000

plt.figure(figsize=(8, 4.5))
plt.hist(balance / 1e6, bins=80, color="#3b82f6", edgecolor="white")
plt.axvline(np.median(balance) / 1e6, color="black", linestyle="--", label="median")
plt.title("Retirement balance after 30 years (10,000 simulations)")
plt.xlabel("Balance (millions of USD)")
plt.ylabel("Simulations")
plt.xlim(0, 6)
plt.legend()
plt.tight_layout()
plt.show()

For more on histograms, see plotting a histogram with matplotlib.

Tips for faster, better simulations

• Vectorize with NumPy: generate all random numbers in one call instead of looping in Python.
• Use np.random.default_rng(), NumPy’s modern generator, rather than the legacy np.random.seed.
• Run until the answer stops changing much; compare results at 10,000 and 100,000 runs.
• Report percentiles and probabilities, not just the mean.

Explore more Python for data and math:

• Random numbers between 0 and 1
• Principal component analysis in Python
• R vs Python for data science
• Python for loop

Frequently asked questions