Estimating Pi with R via MCS-Dart: A very simple example of numerical integration, illustrated and computed in R.
[This article was first published on GivenTheData, 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.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Have you ever played Monte Carlo Dart? If not, read this post and learn how to do it with R and what it can be used for. In fact it is a very easy and prevalent example (which I have come across in a computational economics course last spring semester) that demonstrates the idea behind numerical integration. The aim of this post is not to present a programmatically swell implementation but rather to give an illustrative explanation. Having that said, here’s the problem:
You know of the existence of the number Pi and its role in determining the area of a circle. But, you don’t know its value. How can you estimate the value of Pi using Monte Carlo Simulation, or in other words using random numbers?
Square and circle
First, draw a circle with radius r=1 placed in a square with side length a=2r=2 to illustrate the initial setting. Note that I use the shape package to draw the circle.
library(shape)
r <- 1
par(pty=”s”)
emptyplot(c(0, 0),frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol=”blue”)
lines(x=c(0,1),y=c(0,0), lty=1, col=”red”)
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)
text(x=0.5,y=0.2, labels=expression(“r=1″), col=”red”)
Generally, if a circle with radius r is surrounded by a square with side length a=2r, Pi is equal to the area of the circle divided by the area of the square times 4. Hence, by estimating the fraction of the two areas multiplied by 4 we estimate the value of Pi.
# in case you want to plot the math in R…
emptyplot( xlim=c(0,10), ylim=c(0,11))
text(5,10, expression(paste(A[circle] == r^2 %*% pi, ” and “, A[square] == (2 * r)^{2},”; “, sep=””) ))
text(5,7, expression( paste( frac(A[circle], A[square]) == frac(r^2 %*% pi, (2 * r)^{2}), ” = “, frac(pi, 4), “; thus “, 4 %*% frac(A[circle], A[square]) == pi, sep=””) ))
Playing Monte Carlo Dart
To estimate the fraction of the two areas, we use Monte Carlo Simulation. Imagine it (in this context) as a simulation of randomly throwing darts at the graph above and check for each dart whether it is inside or outside of the circle. The fraction of darts inside the circle is approximately the area of the circle relative to the area of the square (and thus, multiplied by 4, approximately the value of Pi). Intuitively this approximation gets better the greater the number of darts (random numbers) is. But even with a high number of darts the estimated value might in one execution of the experiment be above and in another under the true value of Pi. On average the estimated value is, though, unbiased and pretty close to the real value. Hence, to get a fair estimation we repeat the experiment many times and take the mean of the results.
In fact, the approach described above can be even more simplified. Instead of looking at the whole circle/square it is sufficient to only look in the upper right fourth of the circle/square (as illustrated below). With that target we do exactly the same “dart-throwing” and multiply the resulting fraction by 4.
r <- 1
par(pty=”s”)
emptyplot(frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol=”blue”)
lines(x=c(0,1),y=c(0,0), lty=1, col=”red”)
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)
How to “throw the darts”?
We draw two vectors of random numbers that are ~U(0,1). One for the x coordinates and one for the y coordinates. Each combination of the elements of these two vectors gives us the location where a “dart hit the target”. The remaining question is how to determine whether the dart hit the square inside or outside the circle. The answer comes from the Pythagorean theorem. Imagine the distances from the axes to the point (0.8,0.58) as side lengths a and b of a triangle. If the value of the unknown side length c is greater than the radius r (1), the point is outside the circle; otherwise it is inside (or on) the circle (see illustration below).
points(x=0.8,y=0.58,col=”red”, lwd=2)
# Illustration Pythagorean
par(pty=”s”)
emptyplot(frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol=”blue”)
lines(x=c(0,1),y=c(0,0), lty=1, col=”red”)
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)
points(x=0.8,y=0.58,col=”red”, lwd=2)
lines(x=c(0,0.8), y=c(0,0), col=”red”, lwd=2)
lines(x=c(0.8,0.8), y=c(0,0.58), col=”red”, lwd=2)
lines(x=c(0,0.8), y=c(0,0.58), col=”red”, lty=2,lwd=2)
text(x=0.4,y=0.4,”c”,col=”red”)
title(expression(c == sqrt(a^2 + b^2)))
A function for one round of Monte Carlo Dart
The following function estimates Pi with the approach described above. The only input to the function is the number of “darts” (random coordinates) n. Note that this implementation is programmatically not smart because it uses a loop instead of vectorization. But, I think it is more illustrative for (R-) programming beginners this way.
est.pi <- function(n){
inside <- 0
outside <-0
inside <- 0
outside <-0
pi.est <- 0
# draw random numbers for the coordinates of the “dart-hits”
# draw random numbers for the coordinates of the “dart-hits”
a <- runif(n,0,1)
b <- runif(n,0,1)
b <- runif(n,0,1)
# use the pythagorean theorem
c <- sqrt((a^2) + (b^2) )
# loop: check for each point whether it is inside or outside the circle:
for (i in 1:n) {
if (c[i] > 1) {
# loop: check for each point whether it is inside or outside the circle:
for (i in 1:n) {
if (c[i] > 1) {
outside <- outside + 1
} else {
inside <- inside + 1
}
if (i==n) {pi.est <- (inside/(inside+outside))*4} # compute an estimation of the value of pi
}
return(pi.est)
}
Try that function several times. For the same n you will every time get different estimations of the value of Pi. That is due to the fact, that the whole approach is based on random numbers.
est.pi(n=1000)
est.pi(n=1000)
est.pi(n=1000)
Now, repeat this experiment 1000 times and take the mean of the results as an unbiased estimation of the value of pi. Compare the result with the value of Pi defined in R.
rep <- 1000
pis <- sapply(1:rep, function(x){
est.pi(n=1000)
})
mean(pis)
pi # value of pi defined in R
Finally, get a graphical overview of the estimation.
plot(density(pis), main=expression(paste("Distribution of estimated ",pi,"-values")),cex.main=0.7, cex.axis=0.7, cex.lab=0.7)
# the value of pi defined in R
lines(x=c(pi,pi),y=c(0,10), lty=1,col="blue")
# the estimation
# the estimation
lines(x=c(mean(pis),mean(pis)),y=c(0,10), lty=1,col="red")
What has this to do with numerical integration?
The link between the application of "MCS-Dart" to estimate the value of Pi and the general idea behind numerical integration is straightforward. The core of the approach described above is to draw random numbers (of a certain interval) for coordinates and check, whether the resulting points are under or above a curve. If one is interested in the (approximate) area under a curve of a given function y=f(x) =... in a certain interval and it is very hard to find the integral analytically one can apply exactly the same approach. Draw random numbers for x and y, check if y is greater than f(x), if so, the point (x,y) must be above the curve and otherwise under the curve. Compute the fraction of points under the curve and multiply it by the total area of the rectangle with the x- and y-intervals as side lengths...
To leave a comment for the author, please follow the link and comment on their blog: GivenTheData.
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.