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 ….

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.

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