mcest.pi <- function(N) {
#  Function to estimate pi using uniforms on the unit square
#  and comparing their distance from the origin to 1.  
    unifs <- matrix(runif(2*N),nrow=N)
    dists <- unifs[,1]^2+unifs[,2]^2
    pi.hat <- 4*sum(dists<=1)/N
    c(N, pi.hat, sqrt(pi.hat*(4-pi.hat)/N))
}

mcest.pi.1 <- function(N) {
#  Function to estimate pi as the average of (1-u^2) where
#  the u's are N iid U(0,1) variates.
    unifs <- runif(N)
    pi.tilde <- 4*mean(sqrt(1-unifs^2))
    c(N, pi.tilde, sqrt((32-3*pi.tilde^2)/(3*N)))
}

#
#  The code which I used to generate the tables on pages 45 and 48
#  is as follows

N <- c(10,50,100,500,1000,5000,1e4,5e4,1e5,5e5,1e6)
piests.0 <- matrix(NA, length(N), 4)
piests.1 <- piests.0
options(object.size=1e9) # needed for the large simulations.
seed <- .Random.seed
for (i in 1:length(N)) {
  .Random.seed <- seed
  temp <- mcest.pi(N[i])
  piests.0[i,] <- c(temp, abs(temp[2]-pi))
  .Random.seed <- seed
  temp <- mcest.pi.1(N[i])
  piests.1[i,] <- c(temp, abs(temp[2]-pi))
}  



myqqnorm <- function(data, R=999, alpha=0.05, 
                     ylab="data quantiles") {
#  Function to add confidence envelopes to
#  a normal quantile plot.
#
   z <- (data-mean(data))/sqrt(var(data))
   qqs <- qqnorm(z, plot=F)
# Calculate a normal quantile plot but do not plot it.   

   n <- length(z)
   sims <- rbind(z,matrix(rnorm(R*n), nrow=R))
   sims <- apply(sims, 1, sort)
# Sort each row of sims.  Returns an n x (R+1) matrix.

   sims <- apply(sims, 1, sort)
# Sort each of the n rows. Each row corresponds to an
# order statistic. This will return an (R+1) x n matrix

   rows <- round((R+1)*c(alpha/2, 1-alpha/2))
   env <- sims[rows,]
# Select the rows of sims corresponding to the envelope.

   ylim <- c(min(qqs$x,env), max(qqs$x,env))
   plot(qqs, xlab="Quantiles of Standard normal",
        ylim=ylim, ylab=ylab)
   xvals <- sort(qqs$x)
   lines(xvals, env[1,], lty=2)
   lines(xvals, env[2,], lty=2)
# Make the plot and add lines for the envelope bounds.

   invisible(list(x=xvals, orig=sort(z), env=env))
}

#  This is the code used to generate the plots on pages 57 and 59
normdata <- rnorm(15)
normdata <- (normdata-mean(normdata))/sqrt(var(normdata))
cauchydata <- rcauchy(15)
cauchydata <- (cauchydata-mean(cauchydata))/sqrt(var(cauchydata))

motif("-g 600x800") # Open a graphics window in Unix.  
split.screen(c(2,1))
screen(1)
normq <- qqnorm(normdata, ylim=c(-3,3), ylab="Normal Data Quantiles",
                pch=18)
screen(2)
qqnorm(cauchydata, ylim=c(-3,3), ylab="Cauchy Data Quantiles",
                pch=18)
normq <- sort(normq$x)
seed <- .Random.seed
for (i in 1:10)
{  newz <- sort(rnorm(15))
   screen(1, new=F)
   lines(normq, newz, lty=2)
   screen(2, new=F)
   lines(normq, newz, lty=2)
}
# Plot on page 57.


close.screen(all=T)
split.screen(c(2,1))
screen(1)
.Random.seed <- seed
myqqnorm(normdata, ylab="Normal Data Quantiles")
screen(2)
.Random.seed <- seed
myqqnorm(cauchydata, ylab="Cauchy Data Quantiles")
# Plot on page 59.