Suppose we want to measure the area of a circle. Yeah, yeah, you already know this, it’s \( \pi r^2 \), but let’s say we want to measure a circle of radius 1 with something called Monte Carlo integration:
- We generate pairs \( (x,y) \) where \( x \) and \( y \) are independent and uniformly distributed between 0 and 1.
- We count the fraction of pairs \( (x,y) \) that are inside our circle, by testing each one: the point \( (x,y) \) is inside the circle if \( x^2 + y^2 < 1 \).
This will give us the area of the upper right quadrant of the circle, and we can multiply by four to get something that should be close to \( \pi \):
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def monte_carlo_integration(f,R,N,x_func=None,y_func=None,plotter=False):
x = R.rand(N)
x = x if x_func is None else x_func(x)
y = R.rand(N)
y = y if y_func is None else y_func(y)
inside = f(x,y)
S = sum(inside)
if plotter:
def plotter(ax):
ax.axis('equal')
ax.set_xlim(0,1)
ax.set_ylim(0,1)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
msz = 0.9
ax.plot(x[inside], y[inside], '.', markersize=msz, color='#4090ff')
ax.plot(x[~inside], y[~inside], '.', markersize=msz, color='#e0e0e0')
return S, plotter
else:
return S
def f_circle(x,y):
return x*x + y*y < 1
R = np.random.RandomState(seed=12345)
N = 20000
fig, ax = plt.subplots(figsize=(9,9))
S, plotter = monte_carlo_integration(f_circle, R, N=N, plotter=True)
plotter(ax)
pi_approx = 4 * S / N
print('%f = 4 * %d / %d' % (pi_approx,S,N))3.130800 = 4 * 15654 / 20000
In case you haven’t figured it out, you won’t get that many digits of accuracy from this method. We can get a sense of how accurate it is in two ways, first by running a few iterations, and then by doing some thinking.
import re
def approx_pi(R, N):
x = R.rand(N)
y = R.rand(N)
inside = x*x + y*y < 1
return 4*sum(inside)/N
def print_to_column(something, col, separator=','):
s = str(something)
pattern = re.compile(separator)
line_start = 0
for match in pattern.finditer(s):
this_pos = match.end()
if this_pos - line_start > col:
print(s[line_start:last_pos])
line_start = last_pos
last_pos = this_pos
else:
last_pos = this_pos
print(s[line_start:])
R = np.random.RandomState(seed=12345)
Niter = 16
for N in [20000, 200000, 2000000]:
pi_iterations = [approx_pi(R, N=N) for _ in range(Niter)]
print("%d iterations of N=%d; mean=%.6f, std=%.6f" %
(Niter, N, np.mean(pi_iterations), np.std(pi_iterations)))
print_to_column(pi_iterations, col=80)16 iterations of N=20000; mean=3.139088, std=0.008702 [3.1308, 3.1362, 3.1402, 3.1384, 3.1392, 3.1364, 3.1302, 3.126, 3.154, 3.144, 3.1236, 3.1462, 3.15, 3.133, 3.1494, 3.1478] 16 iterations of N=200000; mean=3.140106, std=0.002850 [3.13944, 3.14074, 3.1382, 3.13858, 3.1451, 3.14092, 3.14186, 3.14538, 3.14142, 3.13566, 3.14234, 3.14196, 3.13778, 3.13588, 3.13644, 3.14] 16 iterations of N=2000000; mean=3.141958, std=0.001117 [3.140182, 3.141652, 3.142454, 3.141808, 3.143592, 3.143868, 3.140966, 3.139914, 3.141036, 3.14187, 3.141084, 3.142752, 3.142714, 3.141974, 3.143414, 3.142046]
Increasing the number of samples by a factor of K decreases the standard deviation of the result by a factor of roughly \( \sqrt{K} \)… so if we want one more digit of \( \pi \), we have to run 100× more trials. Ugh.
But for some kinds of problems, the accuracy requirements aren’t high, and it’s easier to use Monte Carlo integration than some other technique. The usefulness of Monte Carlo increases with higher dimensions: for one-dimensional integration there’s always the trapezoidal rule and Simpson’s Rule, but 2-D or 3-D integration (or higher dimensions in theoretical applications) requires a combinatorial explosion of calculations. Wojciech Jarosz explained it this way in his doctoral dissertation, Efficient Monte Carlo Methods for Light Transport in Scattering Media:
Standard integration techniques exist which converge much faster in one dimension; however, these techniques suffer from the curse of dimensionality, where the convergence rate becomes exponentially worse with increased dimensions. The basic Monte Carlo estimator above can easily be extended to multiple dimensions, and, in contrast to deterministic quadrature techniques, the convergence rate for Monte Carlo is independent of the number of dimensions in the integral. This makes Monte Carlo integration the only practical technique for many high dimensional integration problems, such as those encountered when computing global illumination.
Aside from just running several iterations of the same problem and computing the variance numerically, the other way to figure out how much variation is present in Monte Carlo integration is to consider each iteration as a “coin flip” with a certain probability \( p \); in our example of the area of a circle, it’s \( p = \pi/4 \approx \) 0.7854. If the result of each experiment of \( N \) points is a total count \( S \), then \( S \) has a binomial distribution \( B(N,p) \) which has a mean of \( Np \) and standard deviation \( \sqrt{Np(1-p)} \), so the value \( 4S/N \) has a mean of \( 4p = \pi \) and a standard deviation of \( 4\sqrt\frac{p(1-p)}{N} = \sqrt\frac{\pi(4-\pi)}{N} \approx \frac{1.642}{\sqrt{N}}. \)
p = np.pi/4
print(" N expected standard deviation")
for N in [20000, 200000, 2000000]:
print('%8d %.8f' % (N, 4*np.sqrt(p*(1-p)/N)))N expected standard deviation 20000 0.01161199 200000 0.00367203 2000000 0.00116120
We could use the same technique to compute the area of a superellipse of degree 4 (also known as a “squircle”), shown as the red area in the diagram below: sample random points within unit square A, count the number of points \( S \) which lie within the region of interest, and compute the area of the region as \( 4S/N \) since the full region has reflection symmetry in the x and y axes.
This has a more complicated expression for its area \( A \), involving the gamma function:
$$A = \frac{8\Gamma\left(\frac{5}{4}\right)^2}{\sqrt{\pi}} \approx 3.708149$$
But the Monte Carlo integration experiment here isn’t any more complicated than the one for the circle:
def f_superellipse4(x,y):
return x**4 + y**4 < 1
R = np.random.RandomState(seed=12345)
N = 20000
fig, ax = plt.subplots(figsize=(9,9))
S, plotter = monte_carlo_integration(f_superellipse4, R, N=N, plotter=True)
plotter(ax)
A_approx = 4 * S / N
print('%f = 4 * %d / %d' % (A_approx,S,N))3.700000 = 4 * 18500 / 20000
One way to improve the circle and superellipse area computations is to use stratified sampling techniques which partition the sampling regions so that the resulting subregions are each better matched to their own calculation than the overall problem as a whole. For example, we’ve been generating \( (x,y) \) uniformly within the square region \( 0 < x,y < 1 \), corresponding to the square A below:
But we could just as easily run Monte Carlo experiments for the rectangle B \( (x_1 < x < x_2, 0.5 < y < x_1) \) and rectangle C \( (x_2 < x < 1, 0 < y < 0.5) \) where for the circle, \( x_1 = \frac{\sqrt{2}}{2} \approx 0.7071, x_2 = \frac{\sqrt{3}}{2} \approx 0.866 \), and for the superellipse with degree 4, \( x_1 = \sqrt[4]{\frac{1}{2}} \approx 0.8409, x_2 = \sqrt[4]{1 - (\frac{1}{2})^4} \approx 0.984 \). This excludes rectangular regions in square A which have easy-to-calculate areas — there’s not much point in wasting some of our randomly-generated points when they’re far within the boundary of interest, and we can just compute those separately.
Then the area of each region is equal to
$$\begin{aligned} A &= 4(0.25 + (x_2 - 0.5) + (x_1 - 0.5)^2 + 2A_B + 2A_C) \\ &= 4x_2 + 4x_1{}^2 - 4x_1 + 8A_B + 8A_C \\ &= 4x_2 + 4x_1{}^2 - 4x_1 + 8(x_2-x_1)(x_1-0.5)s_B + 4(1-x_2)s_C \end{aligned}$$
where \( A_B \) and \( A_C \) are the areas of the desired region within regions B and C — both regions get counted twice because their reflections B’ and C’ near the y-axis have the same area — and \( s_B \) and \( s_C \) are the fraction of randomly-generated points within regions B and C.
R = np.random.RandomState(seed=12345)
N = 20000
def scale_coordinate(cmin, cmax):
def f(c):
# c is uniformly distributed between 0 and 1
return cmin + (cmax-cmin)*c
return f
for region, func, exponent in [('circle', f_circle, 2),
('superellipse', f_superellipse4, 4)]:
def f_x(x):
return (1 - x**exponent)**(1/exponent)
x_1 = 0.5**(1/exponent)
x_2 = f_x(0.5)
print('%s (x_0 = %.6f)' % (region, x_0))
fig, ax = plt.subplots(figsize=(9,9))
S_B, plot_B = monte_carlo_integration(func,R,N,
x_func=scale_coordinate(x_1,x_2),
y_func=scale_coordinate(0.5,x_1),
plotter=True)
S_C, plot_C = monte_carlo_integration(func,R,N,
x_func=scale_coordinate(x_2,1),
y_func=scale_coordinate(0,0.5),
plotter=True)
A_B = (x_2 - x_1)*(x_1 - 0.5) * S_B/N
A_C = (1-x_2)/2*S_C/N
area = 4*x_2 + 4*x_1*x_1 - 4*x_1 + 8*A_C + 8*A_B
# Calculate exact areas, for reference
if region == 'circle':
theta_2 = np.arcsin(0.5)
A_B_exact = ((np.pi/4 - theta_2)/2 # area of the circular sector with arc in region B
- (x_2 - x_1)*0.5/2
- (x_1 - 0.5)*x_1/2) # base * height / 2 = area of two triangles outside region B
A_C_exact = (theta_2/2 # area of the circular sector with arc in region C
- 0.5*x_2*0.5) # 0.5 * base * height = area of triangle outside region C
else:
Nsegments = 1000
u = np.arange(Nsegments+1)/Nsegments
x_B = u*(x_2-x_1) + x_1
A_B_exact = np.trapz(f_x(x_B)-0.5,x_B)
y_C = (1-u)*0.5
A_C_exact = np.trapz(y_C,f_x(y_C)-x_2)
print('s_B = %.6f, A_B = %.6f (exact=%.6f)' % (S_B/N, A_B, A_B_exact))
print('s_C = %.6f, A_C = %.6f (exact=%.6f)' % (S_C/N, A_C, A_C_exact))
print('area = %.6f (exact=%.6f)' % (area, 4*x_2 +4*x_1*x_1 - 4*x_1 + 8*A_C_exact + 8*A_B_exact))
plot_C(ax)
plot_B(ax)
theta = np.arange(0,1.001,0.002)*np.pi/2
ax.plot(np.cos(theta)**(2/exponent),
np.sin(theta)**(2/exponent),
color='black',
linewidth=1)
ax.set_title('Stratified sampling, %s' % region)
myblue = '#4090ff20'
ax.add_patch(plt.Rectangle((0, 0), 0.5, 0.5, edgecolor='black', facecolor=myblue))
ax.add_patch(plt.Rectangle((0.5, 0.5), x_1-0.5, x_1-0.5, edgecolor='black', facecolor=myblue))
ax.add_patch(plt.Rectangle((0.5, 0), x_2-0.5, 0.5, edgecolor='black', facecolor=myblue))
ax.add_patch(plt.Rectangle((0, 0.5), 0.5, x_2-0.5, edgecolor='black', facecolor=myblue))
ax.annotate('B',(x_1,0.5),(1,1),textcoords='offset points',fontsize=15)
ax.annotate('C',(x_2,0),(1,1),textcoords='offset points',fontsize=15)
ax.set_ylim(0,1.01)
ax.set_xlim(0,1.01)
circle (x_0 = 0.983995) s_B = 0.543200, A_B = 0.017878 (exact=0.017947) s_C = 0.674500, A_C = 0.045183 (exact=0.045293) area = 3.140165 (exact=3.141593) superellipse (x_0 = 0.983995) s_B = 0.659200, A_B = 0.032157 (exact=0.031997) s_C = 0.799900, A_C = 0.006401 (exact=0.006419) area = 3.709286 (exact=3.708149)
These are more accurate because the randomness makes up a smaller part of the calculation.
An analogy: suppose you are estimating the population of Arizona, which has an area of approximately 295,000 km², and you have 1000 workers. It doesn’t make sense to divide up the geography equally and assign them each 295 km² to cover. Instead, you send the bulk of your workers to estimate the centers of population, where population density is high — Phoenix, Tucson, Yuma, and Flagstaff — and spread out the remaining time among the rural areas, where population density is low. (Unfortunately the political version of this strategy leads to resentment among rural voters whose views are drowned out by the city-dwellers.)
Another analogy: suppose you are trying to determine the quantitative composition of minestrone soup, by taking random samples of the soup. This gets complicated because the soup is liquid and all the damn stuff in it keeps moving around. So you filter the soup with a strainer, and separate it into broth (liquid and small particles), and the solid ingredients. Then you can find the composition of the broth, using techniques for liquids, and the composition of the solid ingredients, using techniques for solids that don’t move. I don’t know what these techniques are, but maybe you get the idea.
Okay, now that we’ve got that straightened out....
Monte Carlo integration isn’t just used for calculation of areas or volumes… we’ve shown examples of integrating the constant 1 over a region, but you can just as easily use it to integrate some function that has spatial variation, by weighting each sample with the integrand at its location.
In other words, suppose I wanted to integrate \( g(x,y) = x^8 + y^4 \) over the unit circle \( x^2 + y^2 < 1 \). Instead of counting the number of points generated randomly within the unit circle, we need to evaluate \( g(x,y) \) for each random point and add up the total, divided by the area “allocated” for each point:
$$\begin{aligned} G &= \int\limits_{x^2+y^2 < 1} g(x,y)\, dx\, dy \\ \\ &\approx \sum\limits_{x^2+y^2 < 1} g(x,y)\, \Delta A \end{aligned}$$
for \( N \) random points \( x,y \) within the square \( |x| < 1, |y| < 1 \) of area 4, in which case \( \Delta A = 4/N \).
It’s also possible to do a better job selecting sample points. Here I’ve done a poor job dividing up my sampling points for the superellipse area: area C is much smaller, so more sampling points should be allocated for area B. And of course there are tighter bounds on the region that could be used; I could have drawn triangles, or used a chain of smaller rectangles along the region boundary. But it’s time to move on.
History of Monte Carlo Methods
Shortly after the end of the Second World War, mathematicians and physicists at Los Alamos used ENIAC, the first general-purpose digital computer, to answer computational questions involving the processes of atomic fission and fusion. Los Alamos physicists Nicholas Metropolis and Stan Frankel had visited the ENIAC team at the University of Pennsylvania, and briefed other staff at Los Alamos, including mathematician Stanislaw Ulam. Ulam later had the idea to utilize statistical techniques:
The first thoughts and attempts I made to practice [the Monte Carlo method] were suggested by a question which occurred to me in 1946 as I was convalescing from an illness and playing solitaires. The question was what are the chances that a Canfield solitaire laid out with 52 cards will come out successfully? After spending a lot of time trying to estimate them by pure combinatorial calculations, I wondered whether a more practical method than “abstract thinking” might not be to lay it out say one hundred times and simply observe and count the number of successful plays. This was already possible to envisage with the beginning of the new era of fast computers, and I immediately thought of problems of neutron diffusion and other questions of mathematical physics, and more generally how to change processes described by certain differential equations into an equivalent form interpretable as a succession of random operations. Later…[in 1946, I] described the idea to John von Neumann and we began to plan actual calculations.
Von Neumann, then at the Institute for Advanced Study at Princeton, took this idea and ran with it, writing to Robert Richtmeyer, the Theoretical Division Leader at Los Alamos, in March 1947:
I have been thinking a good deal about the possibility of using statistical methods to solve neutron diffusion and multiplication problems, in accordance with the principle suggested by Stan Ulam. The more I think about this, the more I become convinced that the idea has great merit.
Von Neumann’s applications of Monte Carlo were used to estimate the explosive yields of fission bombs under development.
Although at the time these were military applications, Monte Carlo methods spread quickly among the scientific community and have been used for solving many different problems ever since.
What Does This Have to Do With Embedded Systems?
What does this have to do with embedded systems? Admittedly, not much, although I did mention Monte Carlo analysis briefly in an article on tolerance analysis, where I used the technique to estimate the tolerance limits of the output of a resistor voltage divider, given the two resistors in the divider have resistances described by a normal distribution. Monte Carlo analysis is used frequently for all sorts of statistical estimates in circuit design.
Stay tuned for more uses of randomness!
References
Los Alamos National Laboratory published a magazine called Los Alamos Science, with a special issue in 1987 dedicated to the memory of Stanislaw Ulam. There are several articles in this issue about Monte Carlo methods, most notably:
-
Nicholas Metropolis, The Beginning of the Monte Carlo Method
-
Roger Eckhardt, Stan Ulam, John von Neumann, and the Monte Carlo Method. (The quotes above from Ulam and von Neumann are from this article.)
-
Gary D. Doolen and John Hendricks, Monte Carlo at Work, which outlines some applications of Monte Carlo in physics.
Metropolis and Ulam published an article in 1949 on the Monte Carlo Method:
- Nicholas Metropolis and Stanislav Ulam, The Monte Carlo Method, Journal of the American Statistical Association, September 1949. (Search for this online and you can find freely available copies.)
Other sources:
-
Wojciech Jarosz, Appendix A: Monte Carlo Integration, part of doctoral dissertation, September 2008.
-
Owen Summerscales, Hitting the Jackpot: The Birth of the Monte Carlo Method, Actinide Resesarch Quarterly, November 2023.
-
Anthony T. Patera and Masayuki Yano, Monte Carlo Integration… in a Nutshell, 2014.
© 2026 Jason M. Sachs, all rights reserved.
