The Monte Carlo is a probability-based method, very popular in Physics circles, to perform numerical estimations of quantities. It relies on the very simple idea of repeated random sampling and is often used to estimate integrals. In practice, what you do is drawing a lot of random numbers, and observing the number of them which respect a certain property, property that will give you the estimate of the quantity you're looking for and that is difficult to calculate analytically.
The method is a very simple and elegant one, and this is why it's one of the silver bullets of Physics. It was proposed by S Ulam while working on military-related projects at the Los Alamos labs in 1940 and became a big contributor to the work of the Project Manhattan, see this historical review paper for a dive into those times. The name is a clear reference to the casino in Monaco.
Monte Carlo estimations are regularly used in problems of
- optimisation
- numerical integration
- drawing from a probability distribution
The basic idea can be visualised by imagining that we can estimate the surface area of a lake by throwing stones across and counting the number of those which fall in the lake with respect to those that fall outside.
The Monte Carlo method is fundamentally based on the law of large numbers (see page).
{% page-ref page="the-law-of-large-numbers.md" %}
A notebook with the code presented here can be seen here.
This is a pedagogical example, cited in many places, for instance Wikipedia. Because of the relation linking
we can easily estimate$$\pi$$via considering a quarter-circle inscribed in a square. The area of the square would simply be$$r^2$$, hence the ratio of the two areas comes down to$$\pi$$if the radius is 1, which is how we would estimate
For other examples, see the reference below.
def f_quartercircle(x):
"""Quarter of a cirle in the first quadrant."""
return np.sqrt(1-x**2)
# get 1000 numbers between 0 and 1, equally spaced, and the quarter circle function on them
x = np.linspace(0, 1, num=1000)
y = [f_quartercircle(item) for item in x]
# loop over the number of extracted points, and extract them uniformly between 0 and 1, both x and y
for n in [100, 1000, 5000, 10000, 20000, 30000, 50000]:
points = np.random.uniform(0, 1, size=(n, 2)) # n points in the plane, randomly (uniformly) extracted in [0,1]
under_points = []
over_points = []
# select if point is below or above the circle
for point in points:
if point[1] <= f_quartercircle(point[0]):
under_points.append(point)
else:
over_points.append(point)
# estimate pi as the ratio of number of points below circle to total
est_pi = float(len(under_points)) / n * 4
# compute the relative error to the real pi, in percentage
perc_err = abs(float(est_pi - np.pi))/np.pi
print('Estimated pi at %d points: %f, with relative error %f' %(n, est_pi, perc_err))This will print the estimation of
plt.figure(figsize=(10, 10))
plt.title('Final estimation plot')
plt.plot(x, y, color='r', lw=3)
plt.plot([point[0] for point in under_points], [point[1] for point in under_points], 'x', alpha=0.5)
plt.plot([point[0] for point in over_points], [point[1] for point in over_points], 'x', alpha=0.5)
plt.show();
Using Monte Carlo for distribution sampling means doing this:
- generate some independent datasets under the condition of interest
- computer the numerical value of the estimator for each dataset, that is, the test statistic
- compute summary statistics over all the dataset test statistic computed
If for instance we want to estimate the mean and standard deviation of a random variable
- Stan Ulam, John Von Neumann and the Monte Carlo method, Los Alamos Science special issue, 1987
- Wikipedia on the Monte Carlo method
- Some useful slides on the method, with examples, from the University of Geneva