Estimating Pi

December 11, 2012
By

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

Recently I’ve been working on some jackknife and bootstrapping problems.  While working on those projects I figured it would be a fun distraction to take the process and estimate pi.  I’m sure this problem has been tackled countless times but I have never bothered to try it using a Monte Carlo approach.  Here is the code that can be used to estimate pi.  The idea is to generate a 2 by 2 box and then draw a circle inside of it.  Then start randomly selecting points inside the box.  This example uses 10,000 iterations using a sample size of 500.  An interesting follow-up exercise would be to change the sample size, mean, and standard deviation to and see what it does to the jackknife variance.

library(bootstrap);
nsims <- 10000;
size <- 500;
alpha <- .05;
theta <- function(x){sd(x)}
pi.estimate <- function(nsims,size,alpha){
out <- as.data.frame(NA);
call <- match.call();
pi.sim <- matrix(NA, nrow=nsims);
for(i in 1:nsims){
x <- runif(size,-1,1);
y <- runif(size,-1,1);
radius <- -1;
x.y <- cbind(x,y);
d.origin <- sqrt(x^2+y^2);
x.y.d <- cbind(x.y,d.origin);
est.pi <- 4*length(x.y.d[,3][x.y.d[,3]<1])/
length(x.y.d[,3])
pi.sim[i] <- est.pi;
}
pi.mean <- mean(pi.sim);
jk <- jackknife(c(pi.sim), theta);
jk.me <- jk$jack.se*qnorm(1-alpha/2);
pi.mean <- pi.mean;
pi.se <- jk$jack.se;
jk.me <- jk.me;
pi.me.u <- pi.mean+jk.me;
pi.me.l <- pi.mean-jk.me;
return(list(pi.mean=pi.mean, pi.se=pi.se,pi.me.u=pi.me.u,pi.me.l=pi.me.l, call=call));
}
pi.estimate(nsims,size,alpha);

To leave a comment for the author, please follow the link and comment on his blog: Statistical Research » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: 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.