Readers should be familiar with some classic approaches to using monte carlo methods to estimate \(\pi\). However the typical textbook lesson only gives you a good estimate to 1 or 2 decimal points and our goal here is to get accuracy upto 10 decimal places in a reasonable amount of time. Let the quest begin!
Let's recall the typical approach for estimating \(\pi\) via monte carlo: First, sample a large number of points uniformly in the unit square \([0,1]^2\), then count the number of points that land in the quarter circle and divide it by the total number of points that we sampled. Scaling by 4 gives us an estimate for \(\pi\):
$$ \hat{\pi}_N = \frac{4}{N} \sum_{i=1}^N \mathbb{I}_{\{x_i^2 + y_i^2 \leq 1\}} $$ To understand how accurate this estimate of \(\pi\) is, note that the indicator term is Bernoulli, i.e. \(\mathbb{I}_{\{x_i^2 + y_i^2 \leq 1\}} \sim Ber(p = \frac{\pi}{4})\). Also recall that Bernoulli random variables have variacne \(p(1-p)\). So the variance of the estimator \(\hat{\pi}_N\) is given by $$ Var(\hat{\pi}_N) = \frac{16}{N^2} Var \left( \sum_{i=1}^N \mathbb{I}_{\{x_i^2 + y_i^2 \leq 1\}} \right) $$ $$ = \frac{16N}{N^2} p(1-p) $$ $$ = \frac{16}{N} p(1-p) $$
With this we can quickly get the order of magnitude of \(N\) and compare it to the variance:
import math as m
def var_pi(n, p = m.pi/4): return (16/n) * p * (1-p)
magnitueds = [10**i for i in range(1,21)]
for num in magnitueds:
print(int(m.log10(num)), f"{var_pi(num):.20f}")
Output:
1 0.26967662132698144717
2 0.02696766213269814749
3 0.00269676621326981475
4 0.00026967662132698145
5 0.00002696766213269815
6 0.00000269676621326981
7 0.00000026967662132698
8 0.00000002696766213270
9 0.00000000269676621327
10 0.00000000026967662133
11 0.00000000002696766213
12 0.00000000000269676621
13 0.00000000000026967662
14 0.00000000000002696766
15 0.00000000000000269677
16 0.00000000000000026968
17 0.00000000000000002697
18 0.00000000000000000270
19 0.00000000000000000027
20 0.00000000000000000003
# The textbook monte carlo estimator:
def mc_pi(n):
x1 = np.random.uniform(0,1,n)
x2 = np.random.uniform(0,1,n)
#y = x1[np.sqrt(x1**2 + x2**2)<1] # The square root is redundant
y = x1[x1**2 + x2**2 < 1]
pi_est = 4 * len(y) / n
return pi_est
N = 10 ** 4
print(mc_pi(N))
Output:
3.1332 # Not great.
# Quasi monte carlo using a Sobol sampler.
def qmc_pi(n):
sobol = qmc.Sobol(d=2, scramble=True, bits=40) # sobol sequnce
x = sobol.random(n) # generate Quasi-Monte Carlo random numbers
x1 = x[:,0]
x2 = x[:,1]
#y = x1[np.sqrt(x1**2 + x2**2)<1] # square root is redundant
y = x1[x1**2 + x2**2 < 1]
pi_est = 4 * len(y) / (n)
return pi_est
# number of samles needed for 10 decimals accuracy.
N2 = m.log2(2.7 * (10 ** 10))
print(N2)
# remember Sobol only accepts powers of 2
print(qmc_pi(2**20))
Output:
34.65224035614973 # Requiring 2**34 samples is far too much for my computer to handle.
3.141551971435547 # okay we got 3.145 which is just 3 decimal places. So better but not great!
def qmc_pi_integral(n):
sobol = qmc.Sobol(d=1, scramble=True)
x = sobol.random(n).flatten()
fx = np.sqrt(1 - x**2)
return 4 * np.mean(fx)
print(qmc_pi_integral(2**20))
Output:
3.1415926553566957 # Great, now we have 3.14159265. But that's only 8 digits!
import numpy as np
def qmc_pi_antithetic_control(n):
sobol = qmc.Sobol(d=1, scramble=True)
# We are actually using half the number of samples we typically would: 2^n --> 2^(n-1)
x = sobol.random(n // 2).flatten()
x_bar = 1 - x
x_full = np.concatenate([x, x_bar])
# Function and control variate
f = np.sqrt(1 - x_full**2)
h = 1 - x_full
Eh = 0.5
# Optimal beta (from full vector)
beta = np.cov(f, h)[0, 1] / np.var(h)
# Apply control variate
f_cv = f - beta * (h - Eh)
# Antithetic average: pair elements [x_i, 1 - x_i]
f_cv_avg = 0.5 * (f_cv[:n//2] + f_cv[n//2:])
return 4 * np.mean(f_cv_avg)
print(qmc_pi_antithetic_control(2**24))
Output:
3.141592653606968 # Great, now we have 3.141592653. 9 digits now!
def qmc_pi_multiple_runs(n, num_runs=10):
"""Run multiple independent estimates and average them"""
results = []
for i in range(num_runs):
# Use different seeds for each run
sobol = qmc.Sobol(d=1, scramble=True)
x = sobol.random(n // 2).flatten()
x_bar = 1 - x
x_full = np.concatenate([x, x_bar])
f = np.sqrt(1 - x_full**2)
h = 1 - x_full
Eh = 0.5
# Stable covariance calculation
f_mean = np.mean(f)
h_mean = np.mean(h)
cov_fh = np.mean((f - f_mean) * (h - h_mean))
var_h = np.mean((h - h_mean)**2)
beta = cov_fh / var_h if var_h > 0 else 0
f_cv = f - beta * (h - Eh)
f_cv_avg = 0.5 * (f_cv[:n//2] + f_cv[n//2:])
results.append(4 * np.mean(f_cv_avg))
return np.mean(results), np.std(results)
print(qmc_pi_multiple_runs(2**24))
Output:
(np.float64(3.141592653577687), np.float64(5.027322973198329e-11))
# Now we have all 10 digits and better yet our std indicates we can consistantly get it
# even though the process is random!