Module 5 · Lesson 6

Simulation Case Study: Monte Carlo Estimate of π

Here’s a classic way to estimate π without any trigonometry: scatter random points uniformly inside a 2x2 square centered at the origin, and check what fraction land inside the unit circle inscribed in that square.

The area of the circle is π * r² = π (since r = 1), and the area of the square is 2 * 2 = 4. So:

fraction inside circle ≈ π / 4   =>   π ≈ 4 * fraction inside circle

The whole “check what fraction land inside the circle” step is just boolean indexing from Module 2 — no loop needed: generate x and y coordinate arrays, build the mask x**2 + y**2 <= 1, and take its mean (since a boolean array’s mean IS the fraction of True values).

rng = np.random.default_rng(seed=0)
n = 100_000
x = rng.uniform(-1, 1, size=n)
y = rng.uniform(-1, 1, size=n)

inside = x**2 + y**2 <= 1
pi_estimate = 4 * inside.mean()
print("pi estimate:", pi_estimate)
print("true pi:    ", np.pi)

plt.figure(figsize=(4, 4))
plt.scatter(x[inside][:2000], y[inside][:2000], s=1, label="inside")
plt.scatter(x[~inside][:2000], y[~inside][:2000], s=1, label="outside")
plt.gca().set_aspect("equal")
plt.legend()
plt.title(f"Monte Carlo π ≈ {pi_estimate:.4f}")
plt.show()

Output:

pi estimate: 3.14736
true pi:     3.141592653589793

Write estimate_pi(n_samples, seed) that runs the Monte Carlo procedure above: sample n_samples uniform (x, y) points in [-1, 1] x [-1, 1] using np.random.default_rng(seed), and return 4 * (fraction of points with x**2 + y**2 <= 1) as a float. Use vectorized operations only — no Python loop over the samples.

💡 Hint

Reuse the exact pattern from the example above: rng.uniform(-1, 1, size=n_samples) for x and y, build the inside mask, return 4 * inside.mean().

✅ Solution

def estimate_pi(n_samples, seed):
    rng = np.random.default_rng(seed)
    x = rng.uniform(-1, 1, size=n_samples)
    y = rng.uniform(-1, 1, size=n_samples)
    inside = x**2 + y**2 <= 1
    return 4 * inside.mean()

🎯 Key Takeaways

np.dot/@ on vectors gives a scalar; np.linalg.norm gives length; cosine similarity measures direction regardless of magnitude.
@ on 2-D arrays is true matrix multiplication ((m,n) @ (n,p) -> (m,p)) — completely different from elementwise *, which needs matching/broadcastable shapes.
Solve Ax = b with np.linalg.solve, not by inverting A — it’s faster and more numerically stable. A zero determinant means A is singular and has no unique solution.
np.linalg.eig gives eigenvalues/eigenvectors satisfying A @ v == λ * v; np.linalg.svd generalizes the idea to non-square matrices.
Use np.random.default_rng(seed) for reproducible, independent random streams — never the legacy global np.random.seed/rand functions.
Random sampling + vectorized boolean indexing (Module 2) is enough to build real simulations, like estimating π with Monte Carlo.

Next up: Module 6 — Performance, I/O & Mini-Projects, where you’ll learn how memory layout affects speed, how to save/load arrays, and put everything from this course together in a few capstone exercises.