Random seeds

A footnote toward Rick Wilkin’s recent article “How to Lie with a Simulation”. (Sit in front of a laptop w/o SAS; have to port all SAS/IML codes into R)

Generated 10 seeds randomly to run Stochastic simulation of Buffon's needle experiment by Rick's method. Hardly converge any of them ….

# Replicate Rick Wicklin's SAS/IML codes for Buffon's needle experiment 
simupi <- function(N, seed) {
  set.seed(seed)
  z <- matrix(runif(N*2, 0, 1), N, 2)
  theta <- pi*z[, 1]
  y <- z[, 2] / 2
  P <- sum(y < sin(theta)/2) / N
  piEst <- 2/P
  Trials <- 1:N
  Hits <- (y < sin(theta)/2)
  Pr <- cumsum(Hits)/Trials
  Est <- 2/Pr
  cbind(Est, Trials, seed)
}

# Generated 10 seeds randomly
seed_vector <- floor(runif(1:10)*10000)

# Each simulation with 50000 iterations
N <- 50000

# Run these 10 simulations
rt <- list()
for (i in 1:length(seed_vector)) {
  rt[[i]] <- simupi(N, seed_vector[i])
}
results <- as.data.frame(do.call("rbind", rt))
results$seed <- as.factor(results$seed)

# Visualize individual simulation results
library(ggplot2)
p <- qplot(x = Trials, y = Est, data = results, geom = "line", 
           color = seed, ylim = c(2.9, 3.5))
p + geom_line(aes(x = Trials, y = pi), color = "Black") 
ggsave("c:/plot1.png")

Averaging all results of the 10 simulations out. Then the curve converges easily. The application of this Monte Carlo simulation in Buffon's needle experiment is explained here by Rick Wicklin.

# Visuazlie the average result
rtmean <- aggregate(Est ~ Trials, data = results, mean)
p <- qplot(x = Trials, y = Est, data = rtmean, geom = "line")
p + geom_line(aes(x = Trials, y = pi), color = "Red")
ggsave("c:/plot2.png")