Bayesian PI

Introduction

While I usually prefer the Day/Month/Year calendar format, I have to admit that the Month/Day/Year format offers something unique: it gives us 3/14, officially recognized as $\pi$ day by the United States and Belize. A common example in sampling/Monte Carlo methods is the estimation of $\pi$. These examples typically involve running many simulations, generating random points, and calculating how often they fall within a circle. It is then shown that, on average, the result approximates the real value of $\pi$.

However, these approaches often stop at point estimation and, at best, showing that multiple independent simulations concentrate around the real $\pi$ rather than providing a formal quantification of uncertainty. Naturally, this brings us to the question of uncertainty quantification, how confident are we in our estimate of at any given point?

This is where Bayesian statistics comes in. Instead of treating our estimate as a single number, Bayesian inference allows us to compute a probability distribution over possible values of $\pi$. This approach lets us quantify how confident we are of our estimate This post extends that approach by applying Bayesian reasoning to the estimation of $\pi$. Just as I previously used Bayesian inference to assess bias in coin flips, here we will construct confidence intervals for that update dynamically as more data becomes available, providing a structured way to quantify uncertainty in our estimates.

Before starting, let us import relevant libraries and set a reproducible seed. In order to celebrate $\pi$ day, the ideal seed would be $\pi$ itself. Unfortunately in numpy seeds need to be integers. Let us settle with a number written from all $\pi$ digits that I can remember off the top of my head.

import numpy as np
from matplotlib import pyplot as plt
np.random.seed(314159)

  

Sampling $\pi$

We sample random points uniformly inside a square and check whether they fall inside a circle of radius 1.

x = np.random.uniform(-1,1,10000)
y = np.random.uniform(-1,1,10000)
rho = np.sqrt(x**2+y**2)
success = rho < 1

Here, success is a boolean array where True indicates that a sampled point falls inside the circle.

The reasoning is that area of the unit circle (radius = 1) is $\pi$, while the area of the square $[−1,1]×[−1,1]$ is 4. Thus, the probability of a uniformly sampled point falling inside the circle is ratio of these areas, i.e., $\pi/4$.

$$\frac{\pi}{4} \approx \frac{\text{\textcolor{green}{Points Inside Circle}}}{\text{Total Points}}$$

Let us test this. That would mean that the mean of success would be approximately $\pi$. In numpy, True is treated as 1, and False is treated as 0 in numerical operations. This means that when we take the mean of the success array, we are computing the fraction of points where success is True, or the fraction of points that landed inside the circle:

4*success.mean()

3.1332

At first glance, this seems promising: the estimate lands reasonably close to $\pi$, correct to one decimal place. But how much should we trust this number? What if we ran this experiment again? Would we get the same result?

Beyond Point Estimates: Bayes

Estimating $\pi$ by sampling random numubers naturally involves randomness. Each run of our simulation produces a slightly different estimate. If we increase the number of samples, our estimates tend to cluster around the true value of $\pi$ but by how much? And how do we express this uncertainty mathematically?

This is where Bayesian inference comes in. Instead of treating our estimate as a single number, we seek to determine a probability distribution over possible values of $\pi$. This approach allows us to:

• Quantify uncertainty – how confident are we in our estimate?

• Update our belief dynamically – as more data becomes available, our estimate of π should become more refined.

Likelihood

Each sampled point either falls inside the circle (success, 1) or outside (failure, 0). If we assume each sample is independent, the number of successes $s$ out of $N$ trials follows a Binomial distribution

$$p(s|q,N) = \binom{N}{s} q^s (1 – q)^{N – s}$$

where $q$, the probability of success is our stand for $\frac{\pi}{4}$.

Posterior

Now that we have our likelihood function, we can apply Bayes’ theorem to update our belief about qq after observing data. The logic follows the same done for our previous blog on the subject

$$p(q|s,N) = p(q) \frac{p(s|q,N)}{p(s|N)} \propto p(q) {p(s|q,N)}$$

Since we work in a limited interval ($q$ is between 0 and 1 and represents $\pi/4$) we started with a uniform prior, meaning $p(q) = 1$ for all $q$. Thus the posterior simplifies to

$$p(q∣s,N) = C q^s (1 – q)^{N – s}$$

where $C$ is a normalization constant ensuring the total probability sums to 1.

def  likelihood(success,q):
	return q**(success.sum())* (1-q)**((1-success).sum())

def  posterior(data,precision=1000000):
	q = np.linspace(0,1,precision)
	dq = q[1]-q[0]

	#non-marginalized posterior,
	nm_posterior = likelihood(data,q) #*prior(q) 

	#denominator integral for marginalization
	integral = np.sum(nm_posterior)*dq
	posterior = nm_posterior/integral

	return q, posterior

Transforming the Posterior into $\pi$

Since $$\pi =4q$$

we can transform samples from $q$‘s posterior to a posterior for $\pi$

$$p(\pi|s,N) = \frac{1}{4}p(q|s,N)$$

To visualize how Bayesian inference improves our estimate of qq, let’s plot its posterior distribution after 10, 100, and 1000 samples

for j in [10,100,1000]:
	
	q, p = posterior(success[:j])
	plt.plot(4*q,p/4,label='{} datapoints'.format(j))
	plt.axvline(np.pi,color = 'k',label=r'Target: $\pi$ ')
	plt.xlabel(r'Estimated $\pi$')
	plt.ylabel('Posterior')
	plt.legend()
	plt.savefig('pi_post')

Interestingly, in our first batch of 10 samples, all points happened to fall inside the circle. This is possible when dealing with randomness. This highlights why uncertainty quantification is essential: if we only relied on a point estimate, we might wrongly assume $\pi = 4$. This is the power of Bayesian inference: as data accumulates, our uncertainty decreases, and the posterior indeed concentrates around the true value of $\pi$.

Visualizing the Bayesian update

To further illustrate this, let’s look at an animation that dynamically shows how our estimate of $\pi$ evolves as more data is collected. (Full code on my GitHub)

Left Panel: The sampled points appear within the circle. We mark in green when fall inside the circle (successes), while red points fall outside (failures).The fraction of points inside the circle is computed at each step.

Right Panel: Posterior distribution for π updated in real-time, with the 95% credible interval shown as a red shaded region. (A 95% credible interval means there’s a 95% probability, given our data, that $\pi$ lies within that shaded area.)

Initially, with only a small number of samples, this interval is broad and skewed towards larger values. However, as more points are collected, the posterior distribution narrows and concentrates tightly around the true value of $\pi$. This visualization highlights the strength of Bayesian methods in quantifying uncertainty, allowing our confidence in the estimate to grow as more data is collected.

Happy $\pi$ day