Last chance! 50% off unlimited learning
Sale ends in
Monte Carlo estimation of pi
simpi(n)
the Monte Carlo estimate of
integer, number of Monte Carlo samples, defaults to 1000
Simon Jackman simon.jackman@sydney.edu.au
A crude Monte Carlo estimate of
Contrast this Monte Carlo method with Buffon's needle and refinements thereof (see the discussion in Ripley (1987, 193ff).
Ripley, Brain D. 1987 [2006]. Stochastic Simulation. Wiley: Hoboken, New Jersey.
seed <- round(pi*10000) ## hah hah hah
m <- 6
z <- rep(NA,m)
lim <- rep(NA,m)
for(i in 1:m){
cat(paste("simulation for ",i,"\n"))
set.seed(seed)
timings <- system.time(z[i] <- simpi(10^i))
print(timings)
cat("\n")
lim[i] <- qbinom(prob=pi/4,size=10^i,.975)/10^i * 4
}
## convert to squared error
z <-(z - pi)^2
lim <- (lim - pi)^2
plot(x=1:m,
y=z,
type="b",
pch=16,
log="y",
axes=FALSE,
ylim=range(z,lim),
xlab="Monte Carlo Samples",
ylab="Log Squared Error")
lines(1:m,lim,col="blue",type="b",pch=1)
legend(x="topright",
legend=c("95% bound",
"Realized"),
pch=c(1,16),
lty=c(1,1),
col=c("blue","black"),
bty="n")
axis(1,at=1:m,
labels=c(expression(10^{1}),
expression(10^{2}),
expression(10^{3}),
expression(10^{4}),
expression(10^{5}),
expression(10^{6})))
axis(2)
Run the code above in your browser using DataLab