Monte Carlo Estimate of Pi

Estimate pi by sampling points in a square, count the fraction inside the inscribed circle, and quantify the error scaling.

Problem

Monte Carlo integration estimates quantities by random sampling. To estimate pi we drop random points into a unit square and count how many fall inside the quarter circle of radius 1; the fraction inside times four approximates pi. It is a clean illustration of how sampling error shrinks.

Sampling

For N uniform points in the unit square, the fraction with x^2 + y^2 <= 1 estimates the area ratio pi/4. The standard error of the estimate falls as 1/sqrt(N), independent of dimension, which is why Monte Carlo wins for high-dimensional integrals.

import numpy as np
rng=np.random.default_rng(7)
for N in [1000,10000,100000,1000000]:
    p=rng.random((N,2))
    inside=(p[:,0]**2+p[:,1]**2<=1).mean()
    est=4*inside
    err=abs(est-np.pi)
    print('N',N,'estimate',round(est,4),'error',round(err,4))

Result

Each tenfold increase in samples cuts the typical error by about a factor of three (sqrt(10)), matching the 1/sqrt(N) law. To gain one more decimal digit of accuracy you need roughly a hundred times more samples, which is the fundamental limitation of plain Monte Carlo. Variance-reduction tricks like stratified or importance sampling improve the constant, not the exponent.

• Convergence is 1/sqrt(N) regardless of dimension, unlike grid methods that suffer the curse of dimensionality.
• The result is a random estimate; report a confidence interval, not just the point value.
• Kronos uses Monte Carlo for neutron transport in breeder-blanket studies, where high dimensionality rules out grids.