simpi: Monte Carlo estimate of pi (3.14159265...)

In pscl: Political Science Computational Laboratory

Monte Carlo estimate of pi (3.14159265...)

Description

Monte Carlo estimation of pi

Usage

simpi(n)

Arguments

n	integer, number of Monte Carlo samples, defaults to 1000

Details

A crude Monte Carlo estimate of \pi can be formed as follows. Sample from the unit square many times (i.e., each sample is formed with two independent draws from a uniform density on the unit interval). Compute the proportion p of sampled points that lie inside a unit circle centered on the origin; such points (x,y) have distance from the origin d = \sqrt{x^2 + y^2} less than 1. Four times p is a Monte Carlo estimate of \pi. This function is a wrapper to a simple C function, bringing noticeable speed gains and memory efficiencies over implementations in native R.

Contrast this Monte Carlo method with Buffon's needle and refinements thereof (see the discussion in Ripley (1987, 193ff).

Value

the Monte Carlo estimate of \pi

Author(s)

Simon Jackman [e-mail]

References

Ripley, Brain D. 1987 [2006]. Stochastic Simulation. Wiley: Hoboken, New Jersey.

Examples

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)

pscl documentation built on May 31, 2023, 5:17 p.m.