Monte Carlo Simulation for Estimating Pi ¶
Monte Carlo methods are widely used for numerical simulations, optimization problems, and probabilistic modeling. One of the classic problems that showcase the power of Monte Carlo simulations is estimating the value of $\pi$. The idea is simple: randomly generate points inside a square and count how many fall within a quarter-circle. The ratio of points inside the quarter-circle to the total number of points gives an approximation of $\pi$.
This approach leverages the relationship between the areas of a quarter circle (with radius $r$) and its bounding square (with side length r):
- Area of quarter circle = $\frac{\pi r^2}{4}$
- Area of square = $r^2$
Therefore, the ratio of these areas is:
$$ \frac{\text{Area of quarter circle}}{\text{Area of square}} = \frac{\pi r^2/4}{r^2} = \frac{\pi}{4} $$
Thus,
$$ \pi = 4 \cdot \frac{\text{Area of quarter circle}}{\text{Area of square}} $$
As we increase the number of random points, our approximation becomes more accurate, demonstrating the law of large numbers in action. While this isn't the most efficient way to calculate $\pi$, it elegantly demonstrates how randomized algorithms can solve deterministic problems and introduces key concepts in parallel computing.
Our approach ¶
To computationally estimate the ratio using the Monte Carlo method, we need to:
- Generate a large number of random points within a unit square (with coordinates from 0 to 1).
- Determine which points lie within a quarter circle.
- Calculate the ratio of points inside the circle to total points generated.
A point $(x, y)$ lies within a quarter circle of radius 1 if the distance from the point to the origin is less than or equal to 1. Using the Pythagorean theorem, this means that a point is inside the quarter circle if:
$$ x^2 + y^2 \leq 1 $$
First, let's get some imports out of the way.
import time
import random
import numpy as np
import matplotlib.pyplot as plt
from typing import override
import ray
if not ray.is_initialized():
ray.init()
2025-07-14 20:03:46,544 INFO worker.py:1888 -- Started a local Ray instance.
Now let's draw a quarter circle using matplotlib.
def quarter_circle_coords(resolution: int = 100) -> tuple[np.ndarray, np.ndarray]:
"""
Compute the x,y coordinates of a quarter‑circle of radius 1
from θ=0 to θ=π/2.
Args:
resolution: Number of points along the curve.
Returns:
(x_circle, y_circle): Arrays of coordinates.
"""
theta = np.linspace(0, np.pi / 2, resolution)
return np.cos(theta), np.sin(theta)
We can view this circle using matplotlib.
xc, yc = quarter_circle_coords(200)
plt.figure(figsize=(4,4))
plt.plot(xc, yc, lw=2, label='Quarter Circle (r=1)')
plt.axis('equal')
plt.xlabel('x'); plt.ylabel('y')
plt.show()
Here you see the smooth arc from (1,0) up to (0,1). This is just the boundary of the region we’ll be sampling inside.
Our Monte Carlo estimate uses randomly‐placed points in the unit square [0,1]×[0,1].
We’ll compare how many of them land inside that blue arc vs. outside—this ratio gives us π.
def unit_square_coords() -> tuple[list[float], list[float]]:
"""
Return the closed polygon (x,y) for the unit square [0,1]×[0,1].
"""
xs = [0., 1., 1., 0., 0.]
ys = [0., 0., 1., 1., 0.]
return xs, ys
xsq, ysq = unit_square_coords()
plt.figure(figsize=(4,4))
plt.plot(xsq, ysq, lw=2, label='Unit Square')
plt.plot(xc, yc, lw=2, label='Quarter Circle (r=1)')
plt.axis('equal')
plt.xlabel('x'); plt.ylabel('y')
plt.show()
def sample_points(n_samples: int) -> int:
"""
Generate uniformly random (x,y) in [0,1] x [0,1].
Returns:
Number of randomly sampled points located inside the circle.
"""
n_inside = 0
for _ in range(n_samples):
x,y = random.random(), random.random()
if x*x + y*y <= 1:
n_inside += 1
return n_inside
def estimate_pi(n_samples: int, n_inside: int) -> float:
"""
Run one Monte Carlo simulation of num_points and return π estimate.
"""
return 4. * (n_inside / n_samples)
def monte_carlo_pi(n_samples: int) -> tuple[float, float]:
"""
Use Monte Carlo to estimate pi.
Args:
n_samples: total random points to generate.
Returns:
Estimate of pi
Total computation time.
"""
t_start = time.time()
n_inside = sample_points(n_samples)
pi_est = estimate_pi(n_samples, n_inside)
t_stop = time.time()
return pi_est, t_stop - t_start
n_samples = 1_000
pi_est, t_delta = monte_carlo_pi(n_samples)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.17200000 Error = -0.03040735 Time = 0.0001 s
n_samples = 100_000
pi_est, t_delta = monte_carlo_pi(n_samples)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.14308000 Error = -0.00148735 Time = 0.0117 s
n_samples = 1_000_000
pi_est, t_delta = monte_carlo_pi(n_samples)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.14072400 Error = 0.00086865 Time = 0.0832 s
n_samples = 100_000_000
pi_est, t_delta = monte_carlo_pi(n_samples)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.14145316 Error = 0.00013949 Time = 8.3804 s
n_samples = 1_000_000_000
pi_est, t_delta = monte_carlo_pi(n_samples)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.14161175 Error = -0.00001909 Time = 83.3244 s
How Parallelization Helps ¶
The Monte Carlo estimation of $\pi$ is a textbook example of an "embarrassingly parallel" problem.
- Complete Independence : Each random sample is independent of all other samples.
- No Shared State : Workers don't need to communicate or synchronize with each other.
- Minimal Overhead : The computation-to-communication ratio is extremely high.
- Linear Scaling : Doubling the computing resources can nearly double the throughput.
- Simple Aggregation : Results from all workers can be combined with a simple sum.
from raygent import Task
class MonteCarloPiTask(Task[list[int], int]):
@override
def do(self, batch: list[int], *args: object, **kwargs: object) -> int:
"""Simulates a Monte Carlo experiment for estimating Pi.
Args:
batch: The number of random points to generate in this task instance.
Returns:
The count of points that fell inside the quarter-circle.
"""
n_inside = sample_points(batch[0])
return n_inside
Here, the method
do
takes an integer parameter representing the number of points to generate in each task.
It then:
- Creates the specified number of random points with x and y coordinates between 0 and 1;
- Checks each point to see if it falls within the quarter-circle (where x² + y² ≤ 1);
- Counts and returns the total number of points that landed inside the quarter-circle.
Each task instance will process its own batch of points independently, allowing the workload to be distributed efficiently across multiple cores or machines.
Running the Task with
TaskManager
¶
The
raygent
framework's
TaskManager
enables us to distribute this workload efficiently:
- We divide our desired total number of samples into smaller batches.
- Each batch becomes an independent task assigned to worker processes.
- Multiple CPU cores process these tasks simultaneously.
- The framework handles all the task scheduling and result collection.
- Results are aggregated to calculate the final π approximation.
This approach can achieve near-linear speedup relative to the number of CPU cores available.
For example, on a 16-core machine, we can potentially compute results up to 16 times faster than using a single core.
Even better, by leveraging Ray's distributed computing capabilities through
TaskManager
, we can scale beyond a single machine to a cluster of computers with minimal code changes.
Our
MonteCarloPiTask
class encapsulates the core algorithm while
TaskManager
handles all the complexities of distributing the work, making parallelization straightforward and effective.
from raygent import TaskManager
from raygent.results.handlers import SumResultsHandler
def estimate_pi_parallel(n_samples: int, n_workers: int) -> tuple[float, float]:
"""Estimates Pi using Monte Carlo simulation with parallel execution.
Args:
n_samples: Total number of samples.
n_workers: Number of independent Monte Carlo simulations to run.
Returns:
Estimated value of π.
"""
manager: TaskManager[list[int], SumResultsHandler[int]] = TaskManager(
MonteCarloPiTask, handler_cls=SumResultsHandler, n_cores=n_workers, in_parallel=True
)
samples_per_worker = n_samples // n_workers
t_start = time.time()
handler = manager.submit_tasks(data=[samples_per_worker] * n_workers, batch_size=1)
n_inside = handler.get()
pi_estimate = estimate_pi(n_samples, n_inside)
t_stop = time.time()
return pi_estimate, t_stop - t_start
pi_est, t_delta = estimate_pi_parallel(n_samples=1_000_000_000, n_workers=8)
print(f"Pi ≈ {pi_est:.8f}")
print(f"Error = {np.pi - pi_est:.8f}")
print(f"Time = {t_delta:.4f} s")
Pi ≈ 3.14162994 Error = -0.00003729 Time = 15.8711 s