⚙️ Setup¶
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(36)
plt.rcParams.update({"figure.dpi":110,"font.size":11,"axes.spines.top":False,"axes.spines.right":False})
INDIGO="#4f46e5"; PINK="#db2777"; GRAY="#94a3b8"
print("ready")
ready
N = 500_000
rolls = rng.integers(1, 7, size=(N, 4)) # N experiments of 4 dice each
at_least_one_six = (rolls == 6).any(axis=1) # success in each experiment
estimate = at_least_one_six.mean()
exact = 1 - (5/6)**4
print(f"Monte Carlo estimate = {estimate:.4f}")
print(f"exact 1 - (5/6)^4 = {exact:.4f}")
Monte Carlo estimate = 0.5168 exact 1 - (5/6)^4 = 0.5177
The estimate lands right on the exact value 0.5177. That is the entire idea: a probability is a long-run frequency (Chapter 30), so simulating the long run and counting is a measurement of the probability. The same loop estimates an expectation, just average the values instead of counting hits.
N = 100_000
x = rng.random(N); y = rng.random(N)
inside = x**2 + y**2 <= 1.0
pi_est = 4 * inside.mean()
print(f"fraction inside quarter circle = {inside.mean():.4f} (theory pi/4 = {np.pi/4:.4f})")
print(f"pi estimate = 4 x fraction = {pi_est:.4f} (true pi = {np.pi:.4f})")
fraction inside quarter circle = 0.7854 (theory pi/4 = 0.7854) pi estimate = 4 x fraction = 3.1414 (true pi = 3.1416)
fig,ax=plt.subplots(figsize=(4.6,4.6))
samp=2500
ax.scatter(x[:samp][inside[:samp]], y[:samp][inside[:samp]], s=4, color=INDIGO, label="inside")
ax.scatter(x[:samp][~inside[:samp]], y[:samp][~inside[:samp]], s=4, color=GRAY, label="outside")
th=np.linspace(0,np.pi/2,200); ax.plot(np.cos(th), np.sin(th), color=PINK, lw=2)
ax.set_aspect("equal"); ax.set_title(f"{samp} darts: pi approx {pi_est:.3f}"); ax.legend(loc="upper right", fontsize=8)
plt.tight_layout(); plt.show()
No formula for pi was used, only random points and an area ratio. This is Monte Carlo integration in miniature: a hard quantity (an area, an integral, an expectation) becomes an average over random samples.
Ns = np.array([100, 300, 1000, 3000, 10_000, 30_000, 100_000, 300_000])
errs = []
for n in Ns:
xs = rng.random(n); ys = rng.random(n)
est = 4*((xs**2+ys**2)<=1).mean()
errs.append(abs(est - np.pi))
for n,e in zip(Ns,errs):
print(f"N={n:7d} |error| = {e:.4f}")
N= 100 |error| = 0.1016 N= 300 |error| = 0.2584 N= 1000 |error| = 0.0936 N= 3000 |error| = 0.0803 N= 10000 |error| = 0.0064 N= 30000 |error| = 0.0129 N= 100000 |error| = 0.0067 N= 300000 |error| = 0.0008
fig,ax=plt.subplots(figsize=(7,3.2))
ax.loglog(Ns, errs, "o-", color=INDIGO, label="observed error")
ax.loglog(Ns, 1.6/np.sqrt(Ns), "--", color=PINK, label="reference 1/sqrt(N)")
ax.set_xlabel("number of samples N"); ax.set_ylabel("|estimate - pi|")
ax.set_title("Monte Carlo error shrinks like 1/sqrt(N)"); ax.legend()
plt.tight_layout(); plt.show()
On log-log axes the error tracks a straight line of slope -1/2: the signature of 1/sqrt(N) convergence. The practical lesson is sobering, getting one more decimal place of accuracy costs about 100 times more samples. Simulation is easy to write but can be expensive to make precise.
N = 200_000
car = rng.integers(0, 3, N) # where the car is
pick = rng.integers(0, 3, N) # the contestant's first guess
# Key fact: after the host reveals a goat, SWITCHING wins exactly when the first pick was wrong.
stay_wins = (pick == car).mean()
switch_wins = (pick != car).mean()
print(f"P(win | stay) = {stay_wins:.4f} (theory 1/3 = {1/3:.4f})")
print(f"P(win | switch) = {switch_wins:.4f} (theory 2/3 = {2/3:.4f})")
P(win | stay) = 0.3362 (theory 1/3 = 0.3333) P(win | switch) = 0.6638 (theory 2/3 = 0.6667)
Switching wins about 2/3 of the time, double the odds of staying. The intuition trap is treating the two remaining doors as 50/50; in fact the host's choice is not random, it leaks information, and switching captures it. A skeptic can argue forever, but 200,000 simulated games end the debate.
sample = rng.normal(50, 12, size=40) # one sample of 40 observations
B = 10_000
boot_means = np.array([rng.choice(sample, size=len(sample), replace=True).mean() for _ in range(B)])
lo, hi = np.percentile(boot_means, [2.5, 97.5])
print(f"sample mean = {sample.mean():.2f}")
print(f"bootstrap 95% CI for mean = [{lo:.2f}, {hi:.2f}]")
sample mean = 49.83 bootstrap 95% CI for mean = [46.94, 52.72]
By resampling the data itself, the bootstrap builds a whole distribution of plausible means and reads a confidence interval straight off the percentiles, no algebra required. It is one of the most quietly powerful tools in data science, and it is pure Monte Carlo.