⚙️ Setup¶
import numpy as np
rng = np.random.default_rng(361)
print("ready")
ready
N=200_000
flips = rng.integers(0,2,size=(N,3)) # 1 = heads
est = (flips.sum(axis=1) >= 1).mean()
exact = 1 - 0.5**3
print(f"estimate = {est:.4f}")
print(f"exact 1 - (1/2)^3 = {exact:.4f}")
estimate = 0.8746 exact 1 - (1/2)^3 = 0.8750
Answer: the estimate matches the exact 7/8 = 0.875. "At least one" is the complement of "none", and P(no heads) = (1/2)³ = 1/8, so the answer is 1 − 1/8 = 7/8. Simulation confirms it without the algebra.
N=500_000
d = rng.integers(1,7,size=(N,2))
est = d.max(axis=1).mean()
print(f"E[max of two dice] estimate = {est:.4f}")
print(f"exact 161/36 = {161/36:.4f}")
E[max of two dice] estimate = 4.4743 exact 161/36 = 4.4722
Answer: E[max] ≈ 4.47 (exactly 161/36 = 4.472). To estimate an expectation you average the values instead of counting hits, the same Monte Carlo loop with .mean() over the quantity itself.
N=200_000; k=23
def has_match():
bdays = rng.integers(0,365,k)
return len(np.unique(bdays)) < k
est = np.mean([has_match() for _ in range(N//20)])
print(f"P(shared birthday, 23 people) estimate = {est:.3f}")
P(shared birthday, 23 people) estimate = 0.506
Answer: about 0.507, just over half, with only 23 people. The count of pairs is C(23,2) = 253, which is what drives the collisions, not the 23 people. The result feels impossible until you simulate it.
N=1_000_000
u = rng.random(N)
est = (u**2).mean() # E[f(U)] approximates the integral over [0,1]
print(f"Monte Carlo integral estimate = {est:.4f}")
print(f"exact 1/3 = {1/3:.4f}")
Monte Carlo integral estimate = 0.3335 exact 1/3 = 0.3333
Answer: the estimate lands near 0.333. For a uniform U on [0,1], E[f(U)] equals the integral of f over [0,1], so averaging f(U) over many samples approximates the integral. This is how Monte Carlo tackles integrals too hard to do analytically, the core trick behind a great deal of computational statistics.
sample = rng.normal(100, 20, size=50)
B=10_000
boot_meds = np.array([np.median(rng.choice(sample, size=len(sample), replace=True)) for _ in range(B)])
lo, hi = np.percentile(boot_meds, [2.5, 97.5])
print(f"sample median = {np.median(sample):.2f}")
print(f"bootstrap 95% CI for median = [{lo:.2f}, {hi:.2f}]")
sample median = 100.83 bootstrap 95% CI for median = [96.45, 111.22]
Answer: the bootstrap resamples the data thousands of times, computes the median of each resample, and reads the 2.5th and 97.5th percentiles as the interval. The median has no tidy standard-error formula, but Monte Carlo does not care, it gives a confidence interval for any statistic you can compute.