# STATS 4CI3/6CI3 Winter 2019
# R code for Assignment 2

# Question 1(e)
rnorm.ar <- function(n) {
#
# Function to use accept-reject method to generate n
# standard normal random variates.
#
  n1 <- n #n1 will count how many observations we still need.
  Z <- rep(NA, n)
  while (n1>0) {
    X <- -log(runif(n1))
    Y <- -log(runif(n1))
    accept <- 2*X >= (Y-1)^2
    n2 <- sum(accept)
    Y <- Y[accept]
    U3 <- runif(n2)
    Z[(n-n1+1):(n-n1+n2)] <- ifelse(U3<0.5, -Y, Y)
    n1 <- n1-n2
  }
  Z
}

# Question 2(c)
gamma.est <- function(N, alpha) {
  #
  # Function to estimate Gamma(alpha) by Monte Carlo using a
  # simulation size of N.
  #
  # The function requires alpha>0.5 to get a finite variance.
  #
  if (alpha<=0.5)
    stop("This method cannot be used with alpha<=0.5")
  if (N<2) 
    stop("The number of Monte Carlo replicates must be greater than 1")
  if (N!=ceiling(N)) {
    warning("N rounded up to ",ceiling(N))
    N <- ceiling(N)
  }
  U <- runif(N)
  X <- -log(U)
  Xa1 <- X^(alpha-1)
  gamma.hat <- mean(Xa1)
  temp <- mean(Xa1^2)
  se <- sqrt((temp-gamma.hat^2)/N)
  c(gamma.hat, se)
}

# Code to examine the algorithm
N <- 10^(2:5)
al <- c(0.51, 0.6, 0.75, 0.95, 1.1, 1.25, 1.5, 2, 3, 4, 5)
set.seed(240219)
results <- matrix(NA, nrow=length(al), ncol=2*length(N))
for (i in 1:length(al))
  for (j in 1:length(N))
    results[i,(2*j-1):(2*j)] <- gamma.est(N[j], al[i])
results <- cbind(gamma(al), results)
rownames(results) <- paste("alpha=",al, sep="")
colnames(results) <- c("Value", paste(c("est","se"),
                                      rep(N, each=2), sep="."))

# Question 3(a)
N <- 100000
set.seed(24022019)
X <- runif(N)
Ihat.mc <- mean(exp(X^2))
se.Ihat.mc <- sqrt((mean(exp(2*X^2))-Ihat.mc^2)/N)
round(c(Ihat.mc, se.Ihat.mc),6)

# Question 3(b)
C <- mean(X^2)
mu <- 1/3
cov <- cov(exp(X^2), X^2)
beta <- 11.25*cov
beta
Ihat.C <- Ihat.mc-beta*(C-mu)
se.Ihat.C <- sqrt(se.Ihat.mc^2+beta^2/(11.25*N)-2*beta*cov/N)
round(c(Ihat.C, se.Ihat.C),6)

# Question 3(c)
Ihat1 <- Ihat.mc
Ihat2 <- mean(exp((1-X)^2))
Ihat.A <- (Ihat1+Ihat2)/2
vI1 <- var(exp(X^2))/N
vI2 <- var(exp((1-X)^2))/N
covI1I2 <- cov(exp(X^2), exp((1-X)^2))/N
se.Ihat.A <- sqrt(vI1/4+vI2/4+covI1I2/2)
round(c(Ihat.A, se.Ihat.A),6)

# Question 3(d)
xx <- seq(0,1,by=0.01)[-1]
alpha <- seq(1,2, by=0.1)
vars <- rep(NA, 10)
for (i in 1:10) 
  vars[i] <- var(exp(xx^2)/dbeta(xx,alpha[i],1))
a1 <- alpha[which.min(vars)]
alpha <- seq(a1-0.1, a1+0.1, by=0.01)[-1]
vars <- rep(NA, 20)
for (i in 1:20) 
  vars[i] <- var(exp(xx^2)/dbeta(xx,alpha[i],1))
alpha[which.min(vars)]

set.seed(24022019)
X.beta <- rbeta(N, 1.24, 1)
W <- 1/dbeta(X.beta, 1.24, 1)
hX <- exp(X.beta^2)
Ihat.IS <- mean(hX*W)
se.Ihat.IS <- sqrt((mean((hX*W)^2)-Ihat.IS^2)/N)
round(c(Ihat.IS, se.Ihat.IS),6)

# Question 4(a)
n <- c(5*(1:10),10*(6:10),200,500,1000,5000,10000)
maxn <- max(n)
N <- 10000
# I will generate all of the random variables in one go and store them in a matrix.
set.seed(24022019)
Zmat <- matrix(rnorm(N*maxn), nrow=R)
alpha <- 0.05
z.alpha <- qnorm(1-alpha/2)
results.4a <- matrix(NA,ncol=3,nrow=length(n))
colnames(results.4a) <- c("n", "Cover", "se(Cover)")

for (i in 1:length(n)) {
  ni <- n[i]
  # the first n[i] observations in each row can be considered the sample.
  zbar <- rowMeans(Zmat[,1:ni])
  sez <- apply(Zmat[,1:ni],1,sd)/sqrt(ni)
  lower <- zbar-z.alpha*sez
  upper <- zbar+z.alpha*sez
  p <- mean(lower<0 & upper >0)
  sep <- sqrt(p*(1-p)/R)
  results.4a[i,] <- c(n[i],p,sep)
}

# Alternative method
results.4a1 <- results.4a
set.seed(24022019)
for (i in 1:length(n)) {
  ni <- n[i]
  # the first n[i] observations in each row can be considered the sample.
  T <- rt(N, n[i]-1)
  p1 <- mean(T<z.alpha & T > -z.alpha)
  sep1 <- sqrt(p1*(1-p1)/R)
  results.4a1[i,] <- c(n[i],p1,sep1)
}

# Question 4(b)
set.seed(24022019)
Zmat.exp <- matrix(rexp(N*maxn), nrow=N)
alpha <- 0.05
z.alpha <- qnorm(1-alpha/2)
results.4b <- matrix(NA,ncol=5,nrow=length(n))
colnames(results.4b) <- c("n", "P(lo<0)","se(P(lo<0))", "Cover", "se(Cover)")

for (i in 1:length(n)) {
  ni <- n[i]
  zbar <- rowMeans(Zmat.exp[,1:ni])
  sez <- apply(Zmat.exp[,1:ni],1,sd)/sqrt(ni)
  lower <- zbar-z.alpha*sez
  upper <- zbar+z.alpha*sez
  p0 <- mean(lower<0)
  se0 <- sqrt(p0*(1-p0)/R)
  p1 <- mean(lower<1 & upper>1)
  se1 <- sqrt(p1*(1-p1)/R)
  results.4b[i,] <- c(ni,p0,se0,p1,se1)
}

# Question 4(c)
n <- c(5*(1:6),10*(4:10),200,500,1000)
maxn <- max(n)
sigma <- c(0.5,1,2,4)
z.alpha <- qnorm(1-alpha/2)
out <- matrix(NA,ncol=5,nrow=length(n))
colnames(out) <- c("n", "P(lo<0)","se(P(lo<0))", "Cover", "se(Cover)")
results.4c <- list("sigma=0.5"=out, "sigma=1"=out, "sigma=2"=out,
                   "sigma=4"=out)

set.seed(24022019)
Zmat.norm <- matrix(rnorm(N*maxn), nrow=N)
for (j in 1:length(sigma)) {
  Ymat <- exp(sigma[j]*Zmat.norm)
  true.mean <- exp(sigma[j]^2/2)
  for (i in 1:length(n)) {
    ni <- n[i]
    ybar <- rowMeans(Ymat[,1:ni])
    sey <- apply(Ymat[,1:ni],1,sd)/sqrt(ni)
    lower <- ybar-z.alpha*sey
    upper <- ybar+z.alpha*sey
    p0 <- mean(lower<0)
    se0 <- sqrt(p0*(1-p0)/R)
    p1 <- mean(lower<true.mean & upper>true.mean)
    se1 <- sqrt(p1*(1-p1)/R)
    results.4c[[j]][i,] <- c(ni,p0,se0,p1,se1)
  }
}