Estimating Pi

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

zv7qrnb
Leave a comment

2 Comments

  1. Carl Witthoft

     /  December 12, 2012

    Seems like a lot of SLOC — by comparison, look at a recent post http://giventhedata.blogspot.com/2012/09/estimating-pi-with-r-via-mcs-dart-very.html

    and even that can be tightened up with some vectorization.

    Reply
  1. True Significance of a T Statistic | Statistical Research

Leave a Reply

Your email address will not be published. Required fields are marked *


6 − four =

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

%d bloggers like this: