I’m late for π day

[This article was first published on Back Side Smack » R Stuff, 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.

It is officially no longer pi day, but I didn’t see this Drew Conway post about estimating pi until just a few minutes ago. Because Google Reader doesn’t show github embeds, I also got to try it without seeing Drew’s solution. The estimation method relies on exploiting the area of a circle.
Area of a circle and an inscribed square
We can use R to generate random numbers for our x and y coordinates and count up the number of x,y pairs inside the circle (or quarter of a circle, in our case). Because frac {pi} {4} r^2 is the area of our quarter circle, the ratio of the 4 times the number of random coordinates within the quarter circle to the total number of random coordinates should converge on pi . This is a very simple Monte Carlo integration. So what do we get?

From pi day

Gets pretty close! The final error was 2.54641*10^{-6} , not too shabby! I’m computing the running sample average, so it isn’t a true Monte Carlo, but it converges well enough. Code is below:

pi.ratio<- function(n) {
x<- runif(n, min=-1,max=1)
y<- runif(n, min=-1,max=1)
return(4*sum(x^2 + y^2 <= 1.0)/n)
}
pi.est<- function(iter) {
sum.pi<-stor.rat<-numeric(0)
for (i in 1:iter) {
stor.rat[i]<- pi.ratio(1000)
sum.pi[i]<- sum(stor.rat[1:i])/i
}
return(sum.pi)
}
view raw piday.R hosted with ❤ by GitHub

To leave a comment for the author, please follow the link and comment on their blog: Back Side Smack » R Stuff.

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.

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)