Estimating Pi by Monte Carlo

Throw random darts at a square, count how many land in the inscribed circle, and read off pi - with a 1/sqrt(N) error.

The idea

A quarter circle of radius 1 has area pi/4; the unit square containing it has area 1. Scatter random points uniformly in the square; the fraction landing inside the circle estimates pi/4, so pi is about 4 times that fraction.

import numpy as np
N=1_000_000
pts=np.random.rand(N,2)
inside=(pts[:,0]**2+pts[:,1]**2)<=1
print(4*inside.mean())   # ~3.1416, wobbles run to run

Accuracy

Monte Carlo error shrinks as 1/sqrt(N), independent of dimension. To gain one more decimal digit you need roughly 100 times more samples - slow. A million points typically gives pi to two or three digits. The estimate is unbiased, and its standard error is about sqrt(p(1-p)/N)*4.

Why bother if it is slow

In one or two dimensions Monte Carlo loses to Simpson's rule badly. Its advantage appears in high dimensions: grid-based quadrature suffers the curse of dimensionality (points grow as n^d), while Monte Carlo error stays 1/sqrt(N) regardless of d. That is why it dominates high-dimensional integrals in physics and finance.

Variance reduction

Practitioners rarely use naive sampling. Importance sampling, stratified sampling, and quasi-random (low-discrepancy) sequences all cut the constant in front of 1/sqrt(N) - or improve the rate - by placing samples more cleverly. The pi demo is the friendly face of a serious toolkit.