Back to all posts
•
draft, sampling, sobol, monte-carlo, pi

Estimating pi with Sobol sequence sampling

Suppose we want to estimate \(\pi\) using a Monte Carlo method. The basic idea is simple: sample many random points in the unit square and count how many fall inside the quarter circle of radius 1.

Here is a quick visual explanation of the idea:

The area of the quarter circle is

\[ \frac{\pi}{4}, \]

while the area of the unit square is 1. So if we sample uniformly in the square, the probability that a point lands inside the quarter circle is

\[ \frac{\pi}{4}. \]

This gives the estimator

\[ \hat{\pi} = 4 \cdot \frac{N_{\text{inside}}}{N_{\text{total}}}. \]

The standard Monte Carlo approach uses random uniform samples. But there is another way: use a quasi-Monte Carlo method based on Sobol sequences.

A Sobol sequence is a low-discrepancy sequence designed to cover the space more evenly than independent random samples. Unlike ordinary random sampling, the points are deterministic and are chosen to fill the domain in a structured way.

This matters because Monte Carlo variance decreases like \(1/\sqrt{n}\) for random sampling, while low-discrepancy sequences can often reduce error for the same number of points, especially in moderate dimensions.

The setup is straightforward. We generate points \((x_i, y_i)\) from a Sobol sequence in \([0, 1]^2\), and then check whether they satisfy

\[ x_i^2 + y_i^2 \le 1. \]

If so, the point lies inside the quarter circle. The estimate is then

\[ \hat{\pi}_{n} = 4 \cdot \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}\{x_i^2 + y_i^2 \le 1\}. \]

This is still a Monte Carlo estimator, but the sampling points are not independent random draws. Instead, they are carefully chosen so that the sample is more evenly distributed in the square.

The benefit is that the estimate can converge faster in practice for many applications, especially when we want a smooth and stable approximation without spending a huge number of points.

What makes Sobol sequences interesting in this context is that they are a more efficient way to explore the unit square. Since the quantity of interest is an area, a well-spread set of points often gives a better estimate than a cloud of purely random points.

This is a natural place to begin when thinking about quasi-Monte Carlo methods: they are not random in the usual sense, but they still produce numerical estimates with strong accuracy.

How does this apply to pipeline integrity

Well, think what the function

\[ x_i^2 + y_i^2 \le 1. \]

can be. For example it can be the limit state for B31G. That means that we can use the same framework to estimate a probability of failure.