estimatePi.R

library(foreach)
library(doParallel)


# Estimate value of pi using simulation
mcPi <- function(n=1e6) 4 * sum( runif(n)^2 + runif(n)^2 < 1) / n


# Set up cluster
cl <- makeCluster(4)
registerDoParallel(cl)

# Run simulation
est <- foreach(icount(128), .combine=c) %dopar% mcPi()
mean(est)

# Plot results
plot(density(est), main=expression(paste("Simulated value of ", pi)))
abline(v=mean(est), col="red", lwd="3")

abline(v=pi, col="black", lwd="3")
text(pi, 0, label=expression(pi), pos=2, cex=2)
text(mean(est), 50, label="Simulated value", pos=if(mean(est) < pi) 3 else 1, 
     cex=1, col="red", srt=90, adj=c(20, 20))


# Stop cluster
stopCluster(cl)