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