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)
(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))
(incf good-ones)))
(* 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:
(run-sim 1000)
3.112
(run-sim 10000)
3.1584
(run-sim 100000)
3.13652
(run-sim 1000000)
3.144964
(run-sim 10000000)
3.1420372
(run-sim 100000000)
3.1415274
And there you have it.