# Project 1, problem 2
#
# F(x,y) = exp( - 1/x - 1/y - 1/(x*y))
#
r2 <- function(n) {
  # Simulate x
  u <- runif(n)
  x <- -1/log(u)
  # Compute w_1 (w_2 = 1 - w_1 not needed)
  w1 <- 1/(1+x)
  alpha <- c(1,2)
  beta <- 1 + 1/x
  u <- runif(n)
  k <- ifelse(u < w1, 1, 2)
  y <- 1/rgamma(n, shape=alpha[k], rate=beta)
  cbind(x,y)
}
XY <- r2(1000)
par(mfrow=c(2,2))
# Scatter plot of realisations from joint
plot(XY, pch=".", xlim=c(0,20), ylim=c(0,20))
# cdf vs ecdf of marginal distribution of x
plot(ecdf(XY[,1]), xlim=c(0, 20), ylim=c(0,1), pch=".")
curve(exp(-1/x), add=TRUE, col="green")
# Joint cdf vs ecdf
n <- 201
x <- seq(0, 20, len=n)
y <- seq(0, 20, len=n)
cdf <- outer(x,y, function(x,y) exp(- 1/x - 1/y - 1/(x*y)))
image(x, y, cdf, col=hcl.colors(20), breaks=seq(0,1, length=21))
ecdf <- outer(x, y, Vectorize(function(x,y) mean(XY[,1]<=x & XY[,2]<=y),vectorize.args = c("x","y")))
image(x, y, ecdf, col=hcl.colors(20), breaks=seq(0,1,length=21))


# Monte-carlo integration
#
w <- .3
mu <- c(1,1)
Sigma <- toeplitz(c(1,.6))
rx <- function(u, mu, Sigma) {
  r <- sqrt(-log(u[,1]))
  theta <- 2*pi*u[,2]
  z <- r*cbind(cos(theta), sin(theta))
  eig <- eigen(Sigma)
  x <- mu + eig$vectors %*% diag(sqrt(eig$values)) %*% t(z)
  t(x)
}
h <- function(x) log(w*exp(x[1] + (1-w)*exp(x[2])))
par(mfcol=c(3,2))
n <- 1024
u <- matrix(runif(n*2), ncol=2)
plot(u, pch=".")
plot(ecdf(u[,1]))
x <- rx(u, mu, Sigma)
plot(x, pch=".")
mean(apply(x, 1, h)) # final MC-estimate

# Randomized quasi-Monte Carlo
u <- qrng::sobol(n, d=2, randomize = "digital.shift")
plot(u, pch=".")
plot(ecdf(u[,1]))
x <- rx(u, mu, Sigma)
plot(x, pch=".")
mean(apply(x, 1, h)) # QMC-estimate


# Harmonic mean estimator of model evidence f(x), toy example
#
x <- 1
# samples from posterior
theta <- rnorm(1e+4, 3/4, sqrt(3/4))
reciproclikelihoods <- 1/dnorm(x, mean = theta, 1)
cummean <- function(x) cumsum(x)/(1:length(x))
cumharmonicmean <- function(x) 1/cummean(x)
par(mfrow=c(2,1))
plot(cummean(reciproclikelihoods), pch=".")
plot(cumharmonicmean(reciproclikelihoods), pch=".")
abline(h=dnorm(x, 0, sqrt(1+3)))


# Project 2, problem 2e
library(coda)
gibbs <- function(y, t, method="gibbs", sigma=.1, mcmc=1e+4) {
  n <- length(y)
  y <- c(NA,y,NA)
  
  chain <- matrix(NA, mcmc, 4)
  accepted <- 0
  
  lambda <- 1
  alpha <- 1
  for (i in 1:mcmc) {
    y[1] <- y[2] - qgamma(runif(1, min=pgamma(y[2], alpha, lambda)), alpha, lambda)
    y[n+2] <- y[n+1] + qgamma(runif(1, min=pgamma(t-y[n+1], alpha, lambda)), alpha, lambda)
    switch(method,
           "gibbs"={
             lambda <- rgamma(1, (n+1)*alpha + 1, y[n+2] - y[1])
             alphap <- rnorm(1, alpha, sigma)
             if (alphap<=0)
               accept <- 0
             else
               accept <- min(1, exp(
                 -(n+1)*(lgamma(alphap) - lgamma(alpha)) 
                 + (alphap - alpha)*((n+1)*log(lambda) - 1 + sum(log(diff(y))))))
             if (runif(1)<accept) {
               alpha <- alphap
               accepted <- accepted + 1               
             }
           },
           "block"={
             alphap <- rnorm(1, alpha, sigma)
             if (alphap<=0)
               accept <- 0
             else {
               lambdap <- rgamma(1, (n+1)*alphap + 1, y[n+2] - y[1])
               accept <- min(1, exp(
                 (alphap - alpha)*(- 1 + sum(log(diff(y))) - (n+1)*log(y[n+2] - y[1]))
                 - (n+1)*(lgamma(alphap) - lgamma(alpha))
                 + lgamma(alphap*(n + 1) + 1) - lgamma(alpha*(n + 1) + 1)
               ))
             }
             if (runif(1)<accept) {
               lambda <- lambdap
               alpha <- alphap
               accepted <- accepted + 1               
             }
           }
    )
    chain[i, ] <- c(alpha, lambda, y[1], y[n+2])
  }
  cat("Acceptance rate =", accepted/mcmc,"\n")
  colnames(chain) <- c("alpha", "lambda", "y0", "ynp1")
  mcmc(chain)
}
y <- c(5.14, 10.15, 18.4, 31.72, 35.4, 36.5, 37.95, 40.82, 64.84, 
       96.8)
t <- 100
chain <- gibbs(y, t, sigma=.7, mcmc=10000, method="gibbs")
summary(chain)$stat
chain <- gibbs(y, t, sigma=1.1, mcmc=10000, method="block")
summary(chain)$stat
pairs(unclass(chain), pch=".")

sigma <- seq(.1, 2, len=20)
effsize <- sapply(sigma, function(sigma) effectiveSize(gibbs(y, t, sigma=sigma, mcmc=50000, method="block")))
effsizegibbs <- sapply(sigma, function(sigma) effectiveSize(gibbs(y, t, sigma=sigma, mcmc=50000, method="gibbs")))
plot(sigma, effsize[1,])
points(sigma, effsizegibbs[1,], col="red")