A Basic Monte-Carlo Simulation in Common Lisp
Monte-Carlo simulations are algorithms that approximate numerical
results by repeatedly sampling some space. Generally, the more samples
you collect, the higher the accuracy of the result. In this blog post
I am going to demonstrate how to use such a technique to approximate
the value of the mathematical constant pi.
The high level idea is as follows. We know the area of a circle can
be found using the formula A = pi * r ^ 2, where r is the radius of
the circle. By using a unit circle (a circle whose radius is 1) the
above equation simplifies itself to A = pi. Therefore if we can come
up with a simple method for approximating the area of a unit circle,
we have a method for approximating pi.
The technique we will use is the following: Start with a 1x1 square,
and draw a quarter-circle inside of it - starting from the top-left
corner down to the bottom-right corner. If we want to approximate the
area of the quarter-circle, we can sprinkle some sand onto the square.
The approximate area is thus the number of grains of sand that lie within
the quarter-circle divided by the total number of grains dropped.
Therefore, to approximate the area of a unit circle, it suffices to
approximate the area of this quarter-circle, and then multiply the
result by 4.
In order to simulate this process on a computer, we first need a
way of generating random grains of sand. The following procedure
will generate two points (one for x, the other for y) between
0 and 1:
(defun random-point ()
"Generate a random point in the range [0, 1)."
(cons (random 1.0)
Running this a few times, you will see values such as
(0.9045708 . 0.368289)
(0.42225754 . 0.63779473)
(0.15361631 . 0.61039114)
and so on.
We next need a routine to tell us whether a given random placement
of a grain of sand lies within the quarter-circle. We will use
the standard euclidean distance formula to determine the distance
from the point to the origin (0, 0). If this distance is less than
the radius (1) of the circle, then we proclaim that that particular
point lies within the boundary:
(defun in-circle (point)
"Return true if the given point is within the circle boundary."
(let ((x (car point))
(y (cdr point)))
(< (sqrt (+ (* x x) (* y y))) 1)))
And finally, in order to generate a good approximation, we need to
run this test many times with many samples. The following function
will run the test for a given number of times, counting the number
of times the random points lie within the quarter-circle. Once done,
it multiplies the ratio by 4, giving us our approximation of pi.
(defun run-sim (iterations)
"Run the simulation for so many iterations."
(declare (optimize (speed 3) (debug 0)))
(let ((good-ones 0))
(dotimes (i iterations)
(when (in-circle (random-point))
(* 4.0 (/ good-ones iterations))))
Running this with larger and larger numbers, you can see our approximation
get closer and closer to the true value of pi:
And there you have it.