Estimating Pi with Monte Carlo Methods

This article demonstrates how to use a Monte Carlo simulation to calculate a value for Pi. An example program is given, written in Erlang, which can be run in a parallel environment. As well as the theory behind the method, we will look at how random numbers can be generated in Erlang and how to time Erlang programs. We will also perform an informal analysis of the results.

Monte Carlo Simulations

A Monte Carlo method is the general name of an algorithm that uses a stochastic process to perform repeated random sampling on a problem. A statistical approach is used such that the numerical result to the problem can then be given by considering the number of ‘hits’ that are achieved with the sampling.

The accuracy of the result, and the confidence in it, will increase as the number of samples we take increases.

More details about Monte Carlo methods can be found on the Wikipedia article and in the article The Beginning of the Monte Carlo Method.

Estimating Pi

A nice example of the use of Monte Carlo methods involves estimating a value for Pi. Let us consider a circle (with a radius of one) which is within a square (with edges of length two). Refer to the following diagram:

Now, by choosing random points within the square, we are able to calculate whether or not each point is within the circle or not. We can consider the analogy of throwing darts ‘randomly’ at a dart board.

Since we know the area of a circle with radius r is given by:

And the area of a square with sides of length x is given by:

And in our example, r=1 and x=2 and we know that the ratio, R, of the area of the circle to the total area of the square is:

Therefore, if we can come up with this ratio by experimentation, we can use it to derive a value for Pi:

Generating Random Numbers in Erlang

Our method requires that we generate pairs of random numbers to identify the x- and y-co-ordinates of random points within the square. We need to make sure that these numbers have a uniform distribution—so that there is an equal chance of hitting any point on the dart board.

Fortunately, Erlang has just such a function: random:uniform/0. It is essentially as simple as calling random:uniform() to generate a random float between 0.0 and 1.0. Details can be found in the manual page.

Seeding

A note at the bottom of the manual page states:

If a process calls uniform/0 or uniform/1 without setting a seed first, seed/0 is called automatically.

Interestingly, if a number of processes are spawned together, they will all produce the same list of ‘random’ numbers. This is because each process is seeding the random number generator with the same default values.

To avoid this, we can seed the random number generator ourselves. The manual recommends using the built-in-function now/0. We can achieve this by calling the following line from the newly spawned process:

apply(random, seed, tuple_to_list(now()))

The function now/0 will always return a unique value, so we can be confident that each process has a uniquely seeded random number generator. Incidentally, I am led to believe that now/0 will actually move ahead of the system time if called frequently enough.

Code Listing

-module(pi).
 
-export([start/2]).
 
%% Start processing on the specified number of processes for the specified
%% number of points. The number of points will be split equally amongst the
%% number of processes.
 
start(Processes, Points) ->
	PointsPerProcess = Points div Processes,
	start_processes(Processes, PointsPerProcess),
	wait_for_results(Processes, Points).
 
 
%% Spawn the required number of processes. Each process is passed the Pid of
%% the current process so that it can send the result back.
 
start_processes(0, _PointsPerProcess) ->
	io:format("Started all processes.~n");
 
start_processes(Processes, PointsPerProcess) ->
	Pid = self(),
	spawn(fun() -> process_points(PointsPerProcess, Pid) end),
	start_processes(Processes - 1, PointsPerProcess).
 
 
%% Wait for the results to come back from each of the processeses. The result
%% from each process comes back containing the count of 'hits'. A total sum
%% of hits is calculated, and then when all the results have come back, we
%% give a calculated value for Pi.
%% 
%% Note: this is heavily reliant on the spawned processes not crashing, and
%% returning their respective counts. (Otherwise we will be waiting forever).
 
wait_for_results(Processes, Points) ->
	wait_for_results(Processes, Points, 0).
 
wait_for_results(0, Points, Sum) ->
	Pi = (Sum / Points) * 4,
	io:format("Pi: ~w~n", [Pi]);
 
wait_for_results(Processes, Points, Sum) ->
	receive
		{result, Count} ->
			wait_for_results(Processes - 1, Points, Sum + Count)
	end.
 
 
%% Setup processing of the points. We seed the random number generator, count
%% the number of 'hits' in the simulation, and then send the result back to
%% the parent process.
 
process_points(Points, Pid) ->
	apply(random, seed, tuple_to_list(now())),
	Count = test_points(Points, 0),
	Pid ! {result, Count}.
 
 
%% Each point is chosen randomly and then tested to see whether it is within
%% the circle. If it is, it is added to the count. Finally, the count is
%% returned.
 
test_points(0, Count) -> Count;
 
test_points(Points, Count) ->
	X = (random:uniform() * 2) - 1,
	Y = (random:uniform() * 2) - 1,
	I = X*X + Y*Y,
	if
		I < 1 ->
			test_points(Points - 1, Count + 1);
		true ->
			test_points(Points - 1, Count)
	end.

Note: I’ve only just started learning Erlang. Any feedback on my (simple) code will be welcomed.

Timing an Erlang Program

In order to experiment with running the program on multiple processors, it would be useful to be able to time the execution of the program. There is a module in Erlang for just this purpose. The function to use is: timer:tc/3. The three parameters that it takes are the module, the function and a list of arguments (these three parameters are commonly referred to as the MFA—module, function, arguments).

The following two lines allow us to time the program for ten thousand points running on one process and two processes, respectively:

timer:tc(pi, start, [1, 10000]).
timer:tc(pi, start, [2, 10000]).

The timer:tc/3 function returns the tuple {Time, Value}, where Time is the time that the function took to execute, in microseconds, and Value is the return value from the function we are timing.

Analysis of Results

I’m currently running this code on my MacBook, so as a result of having two CPU cores, running with any more than two processes doesn’t give any speed benefits. It is nice to see, however, that running with two processes does give a significant speed increase over just a single process:

Times above (and below) are averaged from three runs of each scenario and given in seconds.

Let’s now look at the estimations given for Pi:

We can observe the convergence to the true value of Pi on a graph:

The blue line shows our estimation of Pi, and the red line shows the true value of Pi. Note that the x-axis has a logarithmic scale.

Extensions

A nice extension to my program would be to be able to leave it running, and be able to call a function to see what the current progress is. For example, the number of points tested, and the current estimation of Pi would be given.

It would also be interesting to see how the program behaved on more than two processors. Specifically, when the processors are distributed over a network. I don’t currently have any other machines set up to try this out on.