[This article was first published on ProbaPerception, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.

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’)

---

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.