Playing with Monte Carlo

I use Monte Carlo simulations a lot working on LHCb but I've never really got behind the fancy software we use to try out the basics for myself.  Writing this up as a blog will help me organise my thoughts and who knows, maybe some of you will find this interesting too.

The Monte Carlo method makes use of brute force computational power to estimate answers to complex maths and/or model complex systems using large numbers of random number trials.  At LHCb this means answering questions like: "How many times will two protons colliding in the LHC give me my favourite particle?" Wikipedia suggests a simpler question to answer: "What is the value of π?"

First, draw a circle of radius r, enclosed by a square with sides 2r, so the circle just fits inside: 

The ratio of the areas of this circle and square can be used to calculate π :

So, to find π we have to estimate the areas of the circle and square, actually only the ratio between them.  

Using the Monte Carlo method, I randomly choose N points inside the square and count how many of them, k, are also inside the circle.  If these points are really random (i.e. if the points have an equal chance of landing anywhere inside the square) then the ratio k/N will be an estimate of the ratio of the areas circle/square and 4k/N will be an estimate of π.  

Let's try it with 100 points...
I get π ≈ 3.08.

And again...
3.20, 3.24 and 3.12.

Not bad, but not really close to 3.14159265 either.  We'd all guess that the answer should get better as the number of trials, N, is increased. What quantifiable error should I expect from this estimate?  

For each point, the probability, p, of being inside the circle is just the ratio of the areas circle/square, i.e. π/4.  The probability of not being inside the circle must be 1-p since the probability of falling in the square is 1. This is just binomial statistics from school, so we know the number of points we expect to fall in the circle and the deviation from that number (the error) we can expect in a particular test :

The error we can expect on π is this error on k multiplied by 4/N and so is proportional to 1/√N, in line with what we would all expect.  

With a little re-arranging we can write down how many trials we need to use to get whatever precision, f, we need on our estimate of π :

For a precision f of 10%, for example, we need N = 27 and for 1% we need about N = 2,700.  

To provide the notes for Lu Chao, who claims to have memorised π to 67,890 digits we would need N > 2.7 × 10^(135,780).  This is why Monte Carlo is no substitute for an analytic solution. Let's assume Mr. Chao used a different method...

If any of you are still reading and you want to see some of how Monte Carlo is used at LHCb and in particle physics in general, check out a recent talk from the legendary Rick Field.