# Computing pi by simulation (only for illustration)

#Method 1: Monte Carlo Integration of f(x,y)=1_{x^2+y^2 < 1} on the unit square 
pi.sd <- rep(0,20)
pi.hut <- rep(0,20)
for (i in (1:20)) {
  pi.hut[i] <- 4*mean((runif(10000)^2+runif(10000)^2) < 1)
}
pi.sd <- 0.01*sqrt(pi.hut*(4-pi.hut))
plot(1:20,pi.hut,ylim=c(min(pi.hut-2*pi.sd),max(pi.hut+2*pi.sd)),
     xlab="replication",ylab="confidence interval",
     main="confidence intervals for pi")
segments((1:20),pi.hut-1.96*pi.sd,(1:20),pi.hut+1.96*pi.sd)
abline(h=pi)
mean(pi.sd)

#Same integration with a regular grid
x <- 0.01*((1:100)-0.5)
4*mean(outer(x,x,FUN=function(x,y) (x^2+y^2 < 1)))-pi


#Method 2: Taking r(x,y)=Indicator of a polygon as control variable.
par(mfrow=c(2,2))
m.h <- rep(0,20)
m.r <- rep(0,20)
a <- sqrt(2)-1
x <- 0.001*(0:1000)
plot(x,sqrt(1-x^2),type="l",xaxs="i",yaxs="i")
segments(0,1,1/sqrt(2),1/sqrt(2),col=2)
segments(1,0,1/sqrt(2),1/sqrt(2),col=2)
for (i in (1:20)) {
  u <- runif(10000)
  v <- runif(10000)
  h <- ((u^2+v^2) < 1)
  r <- (v < pmin(1-a*u,(1-u)/a))
  m.r[i] <- mean(r) 
  m.h[i] <- mean(h)
}
plot(m.r,m.h,main="relation between target and control variable")
abline(v=0.5*sqrt(2),h=0.25*pi)
copt <- (1-m.h)/(1-m.r) # estimated regression coefficient.
boxplot(copt,main="regression coefficient")
abline(h=0.7327,col=2) #exact value would be (1-pi/4)/(1-0.5*sqrt(2))=0.7327
pi.hut <- 4*(m.h-copt*(m.r-0.5*sqrt(2)))
pi.sd <- 0.04*sqrt(m.h*(1-m.h) - copt^2*m.r*(1-m.r))
plot(1:20,pi.hut,ylim=c(min(pi.hut-2*pi.sd),max(pi.hut+2*pi.sd)),
     xlab="replication",ylab="confidence interval",
     main="confidence intervals for pi")
segments((1:20),pi.hut-1.96*pi.sd,(1:20),pi.hut+1.96*pi.sd)
abline(h=pi)
mean(pi.sd)

#Method 3: Monte Carlo integration of f(x)=sqrt(1-x^2) on the unit
# interval (obtained from method 1 by conditioning on U)
pi.hut <- rep(0,20)
pi.sd <- rep(0,20)
for (i in (1:20)) {
  x <- sqrt(1-runif(10000)^2)
  pi.hut[i] <- 4*mean(x)
  pi.sd[i] <- 0.04*sqrt(var(x))
}
plot(1:20,pi.hut,ylim=c(min(pi.hut-2*pi.sd),max(pi.hut+2*pi.sd)),
     xlab="replication",ylab="confidence interval",
     main="confidence intervals for pi")
segments((1:20),pi.hut-1.96*pi.sd,(1:20),pi.hut+1.96*pi.sd)
abline(h=pi)
mean(pi.sd)

#same with a regular grid of points
4*mean(sqrt(1 - (0.0001*((1:10000)-0.5))^2))-pi

#Method 4 = Method 3 with antithetic variables
par(mfrow=c(1,1))
pi.hut <- rep(0,20)
pi.sd <- rep(0,20)
for (i in (1:20)) {
  u <- runif(5000)
  x <- 0.5*(sqrt(1-u^2)+sqrt(1-(1-u)^2))
  pi.hut[i] <- 4*mean(x)
  pi.sd[i] <- 4*sqrt(var(x))/sqrt(5000)
}
plot(1:20,pi.hut,ylim=c(min(pi.hut-2*pi.sd),max(pi.hut+2*pi.sd)),
     xlab="replication",ylab="confidence interval",
     main="confidence intervals for pi")
segments((1:20),pi.hut-1.96*pi.sd,(1:20),pi.hut+1.96*pi.sd)
abline(h=pi)
mean(pi.sd)