# STATS 4CI3/6CI3 Winter 2019
# Some further examples of Markov Chain simulations

# Consider sampling from the joint density
# f(x,y) = 8xy 0<x<y<1
N <- 10000
# Accept-Reject sampling from independent uniforms
XY0 <- matrix(runif(8*N*2), ncol=2)
M <- 8 # supremum of f/g
p.accept <- f(XY0[,1], XY0[,2])/M


# First let us consider a Metropolis-Hastings algorithm
# with candidate density g(x,y)=1 0<x<1, 0<y<1

XY <- matrix(NA, N,2)
p <- rep(NA,N)
f <- function(x,y) 
  ifelse((x<y) & (x>0) & (y<1), 8*x*y, 0)

# The rgl package allows for nice 3D plots
library(rgl)
persp3d(f, xlim=c(0,1), ylim=c(0,1))

g <- function(x,y) 
  (x>0)*(x<1)*(y>0)*(y<1)


xy0 <- runif(2)
x0 <- min(xy0)
y0 <- max(xy0)

# Since this is an indpendence MH algorithm we can generate
# all the candidates at once
x <- runif(N)
y <- runif(N)
u <- runif(N)

for (i in 1:N) {
  x1 <- x[i]
  y1 <- y[i]
  p[i] <- f(x1,y1)/f(x0,y0)*g(x0,y0)/g(x1,y1)
  # Note that if x1>y1 then p[i]=0 so we will never
  # accept pairs that are outside the support.
  if (u[i]<p[i]) 
    XY[i,] <- c(x1,y1)
  else
    XY[i,] <- c(x0,y0)
  x0 <- XY[i,1]
  y0 <- XY[i,2]
}

# An alternative guarantees that the accetance probability is >0
# For this we consider the joint distribution given by 
# x~unif(0,1), y|x~unif(x,1)
#
# This corresponds to the joint density
# g(x,y) =1/(1-x) 0<x<y<1
g1 <- function(x,y) 
  ifelse((x>0) & (x<y) & (y<1), 1/(1-x), 0)

mfrow3d(1,2,sharedMouse = TRUE)
persp3d(f, xlim=c(0,1), ylim=c(0,1))
next3d()
persp3d(g1, xlim=c(0,1), ylim=c(0,1))

# this is still an independence MH algorithm but with different
# candidate distribution.  We now need to generate the y from
y.new <- runif(N, x, 1)
XY1 <- matrix(NA, N,2)
p1 <- rep(NA, N)

x0 <- runif(1)
y0 <- runif(1,x,1)

for (i in 1:N) {
  x1 <- x[i]
  y1 <- y.new[i]
  p1[i] <- f(x1,y1)/f(x0,y0)*g1(x0,y0)/g1(x1,y1)
  if (u[i]<p1[i]) 
    XY1[i,] <- c(x1,y1)
  else
    XY1[i,] <- c(x0,y0)
  x0 <- XY1[i,1]
  y0 <- XY1[i,2]
}

# A Gibbs Sampler from the joint density 
# f(xy)=8xy 0<x<y<1

# We can easily generate Y0 from the target marginal
# distribtion of Y since that has cdf
# F(y)=y^4 0<y<1

Y0 <- runif(1)^(1/4)

N <- 100000
U <- matrix(runif(2*N),ncol=2)
XY2 <- matrix(NA, nrow=N, ncol=2)
XY2[1,1] <- Y0*sqrt(U[1,1])
XY2[1,2] <- sqrt(X[1]^2+U[1,2]*(1-X[1]^2))
for (t in 2:N) {
  XY2[t,1] <- XY2[t-1,2]*sqrt(U[t,1])
  XY2[t,2] <- sqrt(XY2[t,1]^2+U[t,2]*(1-XY2[t,1]^2))
}

# Of course in this case we do not actually need to
# use MCMC at all since we can generate from the 
# target directly using the marginal for Y and the
# conditional for X given Y.
XY3 <- cbind(NA,runif(N)^1/4)
XY3[,1] <- XY3[,2]*sqrt(runif(N))

# Simulate from a multinomial distribution conditional 
# using the Gibbs Sampler.

# The full conditionals are all binomial
gibbs.multinomial <- function(N, n, p, N0=0) {
  d <- length(p)
  p <- p/sum(p)
  X <- matrix(NA, nrow=N+N0+1, ncol=d)
  X[1,] <- round(p*n)
  for (i in 1:(N+N0)) {
    X[i+1,] <- X[i,]
    for (j in 1:(d-1))
      X[i+1,j] <- rbinom(1,n-sum(X[i+1,-c(j,d)]),
                         p[j]/(1-sum(p[-c(j,d)])))
    X[i+1,d] <- n-sum(X[i+1,-d])
  }
  X[-(1:(N0+1)),]
}

n <- 20
p <- (1:4)/10
X <- gibbs.multinomial(100000,n,p, N0=1000)
# The theoretical mean of this multinomial is 
n*p

# The empircal estimate from our chain is
colMeans(X)

# The theoretical covariance matrix is 
n*(diag(p)-outer(p,p))

# The empirical estimate from our chain is
cov(X)