[This article was first published on Wiekvoet, 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.

What irritates me, is the low efficiency of this estimate. What do you get for 10 000 simulations? Probably, but not even certain, the first two digits.

In the past years I have been thinking how to get that more efficient, but that is not obvious. For instance, it is possible to use the three dimensional equivalent, a ball:

This is even worse, the variation is higher.

At some point I thought this is due to the limited information in such a calculation, it is binomial and one simulation gives one bit of information. And it could be more simple. If the first random number is known, say y, then all second random numbers over sqrt(1-y²) give distance larger than 1, while the remainder gives distance less than 1. Thus should pi be equal to the mean of random numbers transformed like sqrt(1-y²)?
pin <- function(N=1000) {
  4*sum(sqrt(1-runif(N)^2))/N
}
summary(sapply(1:1000,function(x) pin(10000)))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.113   3.135   3.142   3.141   3.147   3.171 
These numbers are closer, but there are additional calculations. Hence the number of simulations should be adapted to reflect the work done. Luckily we have microbenchmark() to calibrate this. After a bit of experimenting, these are the number of simulations giving roughly equivalent computation times.
microbenchmark(pi2d(10000),pi3d(6666),pin(22000))

It could obviously be thought that the random numbers are not needed. An integration can be done. But that is much less fun.

---

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.