# STAT 4CI3/6CI3 Winter 2019
# R code for Assignment 3

# Question 1(a)
rgamma.1a <- function(N,k,beta) {
#
# Function to generate N gamma random variates as a sum of
# independent exponentials which are generated from uniforms.
#
# Arguments
#   N: Sample size required
#   k: The shape parameter of the gamma which must be a
#      positive integer
#   beta: The scale parameter which must be positive
#
  if (k%%1 != 0) stop("k must be an integer")
  if (k<0 | beta<0) stop("Both parameters must be positive")
  U <- matrix(runif(N*k), ncol=k)
  Y <- -beta*log(U)
  rowSums(Y)
}


# Question 1(b)

alpha <- 3.2
beta <- 2
k <- floor(alpha)
b <- alpha*beta/k
x.tilde <- alpha*beta
M <- dgamma(x.tilde,alpha,scale=beta)/dgamma(x.tilde,k,scale=b)
N <- 10000

set.seed(21032019)
Y.1b <- rgamma.1a(N,k,b)
U.1b <- runif(N)
p.1b <- dgamma(Y.1b,alpha,scale=beta)/(M*dgamma(Y.1b,k,scale=b))
accept.1b <- U.1b < p.1b
X.1b <- Y.1b[accept.1b]

# Question 1(c)

x0 <- rgamma(1,alpha,scale=beta)
X.1c <- rep(NA, N)
accept.1c <- rep(NA, N)
p.1c <- rep(NA,N)
for (i in 1:N) {
  if (i ==1) x <- x0
  else x <- X.1c[i-1]
  y <- Y.1b[i]
  p.1c[i] <- dgamma(y,alpha,scale=beta)/dgamma(y,k,scale=b)*
               dgamma(x,k,scale=b)/dgamma(x,alpha,scale=beta)
  p.1c[i] <- min(p.1c[i],1)
  accept.1c[i] <- U.1b[i]< p.1c[i]
  if (accept.1c[i]) X.1c[i] <- y
  else X.1c[i] <- x
}

# Question 2(a)
rnorm.hastings <- function(N, x0, delta) {
#
# An implementation of the Hastings random walk algorithm
# to generate a chain of N standard normal random variates
# starting at x0 and using uniform (-delta, delta) step sizes
# for the candidate at each iteration.  The algorithm will
# return the observed chain as well as the corresponding
# acceptance chain.
#
   eps <- runif(N, -delta, delta)
   U <- runif(N)
   X <- rep(x0, N+1)
   accept <- rep(NA, N)
   for (i in 1:N) {
     y <- X[i]+eps[i]
     accept[i] <- U[i] < dnorm(y)/dnorm(X[i])
     if (accept[i]) X[i+1] <- y
     else X[i+1] <- X[i]
   }
   list(sample=X[-1], accept=accept)
}

# Question 2(b)
N <- 10000
N0 <- 5000
set.seed(18032019)
out.2b <- rnorm.hastings(N+N0, 0, 1)
X.2b <- out.2b$sample[-(1:N0)]

# Code to create a pdf file with the first plot
pdf("A3Q2b1.pdf", width=6, height=9)
par(mfcol=c(3,1))
plot(X.2b, xlab="Iteration Number")
hist(X.2b, breaks=100, prob=TRUE)
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)), col="red", lwd=2)
acf(X.2b)
dev.off()

# Code for the second plot of the auto-correlation function.
pdf("A3Q2b2.pdf", width=6, height=4)
acf(X.2b, lag.max=1000)
dev.off()

# Question 3(c)
set.seed(180320191)
out.2c.1 <- rnorm.hastings(N+N0, 0, 0.1)

set.seed(180320191)
out.2c.2 <- rnorm.hastings(N+N0, 0, 0.5)

set.seed(180320191)
out.2c.3 <- rnorm.hastings(N+N0, 0, 2)

set.seed(180320191)
out.2c.4 <- rnorm.hastings(N+N0, 0, 10)

# Acceptance rates
cbind(c(0.1,0.5,1,2,10),
      c(mean(out.2c.1$accept[-(1:N0)]),
        mean(out.2c.2$accept[-(1:N0)]),
        mean(out.2b$accept[-(1:N0)]),
        mean(out.2c.3$accept[-(1:N0)]),
        mean(out.2c.4$accept[-(1:N0)])))

# Some summary statistics
stats <- function(x) 
    c(mean=mean(x), sd=sd(x), median=median(x))
rbind(stats(out.2c.1$sample[-(1:N0)]),
            stats(out.2c.2$sample[-(1:N0)]),
      stats(out.2b$sample[-(1:N0)]),
      stats(out.2c.3$sample[-(1:N0)]),
      stats(out.2c.4$sample[-(1:N0)]))

pdf("A3Q2c1.pdf", width=6, height=9)
par(mfcol=c(3,1))
plot(out.2c.1$sample[-(1:N0)], main="delta=0.1", 
     ylab="X", xlab="Iteration Number")
hist(out.2c.1$sample[-(1:N0)], breaks=100, prob=TRUE, main="Histogram")
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)), col="red", lwd=2)
acf(out.2c.1$sample[-(1:N0)], lag.max=200, main="Auto-correlation Function")
dev.off()

pdf("A3Q2c2.pdf", width=6, height=9)
par(mfcol=c(3,1))
plot(out.2c.2$sample[-(1:N0)], main="delta=0.5", 
     ylab="X", xlab="Iteration Number")
hist(out.2c.2$sample[-(1:N0)], breaks=100, prob=TRUE, main="Histogram")
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)), col="red", lwd=2)
acf(out.2c.2$sample[-(1:N0)], lag.max=200, main="Auto-correlation Function")
dev.off()

pdf("A3Q2c3.pdf", width=6, height=9)
par(mfcol=c(3,1))
plot(out.2c.3$sample[-(1:N0)], main="delta=2", 
     ylab="X", xlab="Iteration Number")
hist(out.2c.3$sample[-(1:N0)], breaks=100, prob=TRUE, main="Histogram")
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)), col="red", lwd=2)
acf(out.2c.3$sample[-(1:N0)], lag.max=200, main="Auto-correlation Function")
dev.off()

pdf("A3Q2c4.pdf", width=6, height=9)
par(mfcol=c(3,1))
plot(out.2c.4$sample[-(1:N0)], main="delta=10", 
     ylab="X", xlab="Iteration Number")
hist(out.2c.4$sample[-(1:N0)], breaks=100, prob=TRUE, main="Histogram")
lines(seq(-4,4,by=0.01), dnorm(seq(-4,4,by=0.01)), col="red", lwd=2)
acf(out.2c.4$sample[-(1:N0)], lag.max=200, main="Auto-correlation Function")
dev.off()

# Question 3(a)
rbinorm.mh <- function(N, mu1, mu2, s1, s2, rho) {
#
# Function to run a Metropolis-Hastings algorithm to
# generate bivariate normals using a random walk with
# independent standard normal innovations.
#
# Arguments:
#  N - length of the chain
#  mu1, mu2 - marginal means
#  s1, s2 - marginal standard deviations
#  rho - correlation
#
   dbinorm <- function(x, mu1, mu2, s1, s2, rho) {
   # Local function to evaluate the bivariate normal
   # density up to the constant which will cancel in
   # the acceptance probabilities anyway.
       z1 <- (x[1]-mu1)/s1
       z2 <- (x[2]-mu2)/s2
       exp(-1/(2*(1-rho^2))*(z1^2-2*rho*z1*z2+z2^2))
   }

   x0 <- rnorm(2, c(mu1, mu2), c(s1,s2))
   eps <- cbind(rnorm(N,0,s1), rnorm(N,0,s2))
   U <- runif(N)
   X <- matrix(NA, nrow=N+1, ncol=2)
   X[1,] <- x0
   for (i in 1:N) {
      y <- X[i,]+eps[i,]
      p <- dbinorm(y,mu1,mu2,s1,s2,rho)/
             dbinorm(X[i,],mu1,mu2,s1,s2,rho)
      if (U[i] < p) X[i+1,] <- y
      else X[i+1,] <- X[i,]
   }
   X[-1,]
}

# Question 3(c)
rbinorm.gibbs <- function(N, mu1, mu2, s1, s2, rho) {
#
# Function to run a Gibbs sampling algorithm to
# generate bivariate normals.
#
#
# Arguments:
#  N - length of the chain
#  mu1, mu2 - marginal means
#  s1, s2 - marginal standard deviations
#  rho - correlation
#
   x0 <- rnorm(2, c(mu1, mu2), c(s1,s2))
   # calculate the conditional standard deviations which do not
   # depend on the current state of the chain.
   s1c <- s1*sqrt(1-rho^2)
   s2c <- s2*sqrt(1-rho^2)
   X <- matrix(NA, nrow=N+1, ncol=2)
   X[1,] <- x0
   for (i in 1:N) {
      X[i+1,1] <- rnorm(1, mu1+rho*s1/s2*(X[i,2]-mu2), s1c)
      X[i+1,2] <- rnorm(1, mu2+rho*s2/s1*(X[i+1,1]-mu1), s2c)
   }
   X[-1,]
}

# Question 3(d)
set.seed(19032019)
X.q3.mh <- rbinorm.mh(10000, -1, 1, 1, 2, -0.5)
set.seed(19032019)
X.q3.gibbs <- rbinorm.gibbs(10000, -1, 1, 1, 2, -0.5)
colMeans(X.q3.mh)
apply(X.q3.mh,2,sd)
cor(X.q3.mh[,1], X.q3.mh[,2])
colMeans(X.q3.gibbs)
apply(X.q3.gibbs,2,sd)
cor(X.q3.gibbs[,1], X.q3.gibbs[,2])

# Plot of the 4 histograms
pdf("A3Q3d.pdf")
par(mfrow=c(2,2))
hist(X.q3.mh[,1], breaks=100, prob=TRUE, main="Metropolis", xlab="X")
curve(dnorm(x,-1,1), c(-5,3), col="red", lwd=2, add=TRUE)
hist(X.q3.mh[,2], breaks=100, prob=TRUE, main="Metropolis", xlab="Y")
curve(dnorm(x,1,2), c(-7,9), col="red", lwd=2, add=TRUE)
hist(X.q3.gibbs[,1], breaks=100, prob=TRUE, main="Gibbs", xlab="X")
curve(dnorm(x,-1,1), c(-5,3), col="red", lwd=2, add=TRUE)
hist(X.q3.gibbs[,2], breaks=100, prob=TRUE, main="Gibbs", xlab="Y")
curve(dnorm(x,1,2), c(-7,9), col="red", lwd=2, add=TRUE)
dev.off()

par(mfrow=c(2,2))
acf(X.q3.mh[,1], main="M-H chain 1")
acf(X.q3.mh[,2], main="M-H chain 2")
acf(X.q3.gibbs[,1], main="GS chain 1")
acf(X.q3.gibbs[,2], main="GS chain 2")

# Question 4

# Question 4(d)
N <- 10000
N0 <- 2000
mu <- 16
alpha <- 2
beta <- 4
x <- 2

Gibbs.Q4d <- matrix(NA, nrow=N+N0, ncol=2)
set.seed(21032019)
n <- rpois(1,mu)
th <- rbeta(1, alpha, beta)
for (i in 1:(N+N0)) {
  n <- x+rpois(1,mu*(1-th))
  th <- rbeta(1, x+alpha, n-x+beta)
  Gibbs.Q4d[i,] <- c(n, th)
}
Gibbs.Q4d <- Gibbs.Q4[-(1:N0),]
pdf(file="A3Q4d.pdf", width=6, height=6)
par(mfrow=c(2,1), mar=c(3,2,2,1)+0.1)
plot(Gibbs.Q4d[,1], main="Gibbs Sampler for n")
plot(Gibbs.Q4d[,2], main="Gibbs Sampler for theta")
dev.off()

colMeans(Gibbs.Q4d)
colVars(Gibbs.Q4d)
cor(Gibbs.Q4d[,1], Gibbs.Q4d[,2])

# Question 4(e)
Gibbs.Q4e <- matrix(NA, nrow=N+N0, ncol=3)
set.seed(21032019)
n <- rpois(1,mu)
th <- rbeta(1, alpha, beta)
for (i in 1:(N+N0)) {
  x <- rbinom(1, n, th)
  n <- x+rpois(1,mu*(1-th))
  th <- rbeta(1, x+alpha, n-x+beta)
  Gibbs.Q4e[i,] <- c(x, n, th)
}
Gibbs.Q4e <- Gibbs.Q4e[-(1:N0),]

colMeans(Gibbs.Q4e)
colVars(Gibbs.Q4e)
cor(Gibbs.Q4e[,2], Gibbs.Q4e[,3])