Estimation of the number PI – A Monte Carlo simulation

October 15, 2012
By

(This article was first published on ProbaPerception, and kindly contributed to R-bloggers)


How to estimate PI when we only have R and the formula for the surface of a circle (Surface = PI * r * r)?
The estimation of this number has been one of the greatest challenge in the history of mathematics. PI is the ratio between a circle’s circumference and diameter. The first estimation I have heard about has been done by an Egyptian 3700 years ago. His estimation was around 3.16. Which is fairly accurate for a 3700-year-old estimation.
In this post I propose an easy way to estimate PI. Though it is far from being the most efficient estimation ever made, it is a good opportunity to introduce Monte Carlo simulation in the blog.
Monte Carlo simulations are methods to estimate results by repeating a random process. The limit of this method is the source of randomness in the results. By construction of these methods, it cannot be mathematically proved, but only “confidence interval” results. Besides, most of these methods are very slow if we want to obtain a fairly accurate result. The optimization of Monte Carlo is an amazing area which is really useful in many disciplines. 
Again, the method I program here, is far from being efficient, it is only a basic of the Monte Carlo theory and I think one of the easiest illustrations.
Let’s explain how it works.
As we know, if r is the radius of a circle, then the surface of this circle is PI * r * r. Therefore, if we know the radius of a circle, and if we can estimate the surface of this circle, we can estimate PI. How can we estimate the surface of a circle with a random method. The answer is found in the square! Yes, in the square. Especially in the following figure. 
A circle in a square is our solution. Indeed, this is very easy to know the surface of a square. It is the square of the length of one of its sides. Do you see where this is leading us? What if we were able to estimate the rate between the surface of the circle and the surface of the square? Then we would know the surface of the circle.
Indeed, 
surfaceCircle = (surfaceCircle / surfaceSquare) * surfaceSquare. 
Finally, we just need to estimate (surfaceCircle/surfaceSquare). And this is where Monte Carlo is useful. If we uniformly distribute k random variables in the square and count how many of them are in the circle, then we have a proxy of the quantity surfaceCircle/surfaceSquare.
We would obtain something like the next graph. Then the proportion between the number of blue spots and the total number of spots is equal to surfaceCircle/surfaceSquare.
I have done this experience with R. The next plot is the estimation of PI we can obtain with this method. The red line is the real value of PI. The black line is my estimator of PI according to the number of simulated random variables. As we can see, the accuracy tends to increase with the size of our sample. After 10,000 simulations, we obtain an estimation of PI of 3.1444 which is not so good. Indeed, you can already find in the Old Testament a better estimation of PI!


Even though it is a weak estimation, Monte Carlo simulation is a very powerful method. Indeed, if the circle is now a well known mathematical object, all the mathematical concepts are not as well understood. The advantage of Monte Carlo simulation is that we do not need any analytic study of the object. This is why Monte Carlo is used in many applied areas.
 The code (R):
simulation = function(long){
  c = rep(0,long)
  numberIn = 0
  for(i in 1:long){
    x = runif(2,-1,1)
    if(sqrt(x[1]*x[1] + x[2]*x[2]) <= 1){
      numberIn = numberIn + 1
    }
    prop = numberIn / i
    piHat = prop *4
    c[i] = piHat
  }
  return(c)
}
size = 1000
res = simulation(size)
ini = 1
plot(res[ini:size], type = ‘l’)
lines(rep(pi, size)[ini:size], col = ‘red’)

To leave a comment for the author, please follow the link and comment on their blog: ProbaPerception.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)