###################################################
### chunk number 1: libs
###################################################
options(width = 65) 
library(spBayes)
library(MASS)


###################################################
### chunk number 2: synData
###################################################
set.seed(1)
n <- 50
x <- runif(n, 0, 100)
y <- runif(n, 0, 100)

D <- as.matrix(dist(cbind(x, y)))

phi <- 3/50
sigmasq <- 50
tausq <- 20
mu <- 150

s <- (sigmasq*exp(-phi*D))
w <-  mvrnorm(1, rep(0, n), s)
Y <- mvrnorm(1, rep(mu, n) + w, tausq*diag(n))
X <- as.matrix(rep(1, length(Y)))


###################################################
### chunk number 3: adapMCMC
###################################################
#IG sigma^2 and tau^2
a.sig <- 2 
b.sig <- 100
a.tau <- 2
b.tau <- 100

#U phi
a.phi <- 3/500
b.phi <- 3/3

logit <- function(theta, a, b){
  log((theta-a)/(b-theta))
}

logit.inv <- function(z, a, b){
  b-(b-a)/(1+exp(z))
}

##metrop target
target <- function(theta){
  
  mu.cand <- theta[1]
  sigmasq.cand <- exp(theta[2])
  tausq.cand <- exp(theta[3])
  phi.cand <- logit.inv(theta[4], a.phi, b.phi)

  Sigma <- sigmasq.cand*exp(-phi.cand*D)+tausq.cand*diag(n)
  SigmaInv <- chol2inv(chol(Sigma))
  logDetSigma <- determinant(Sigma, log=TRUE)$modulus[1]
  
  out <- (
          -(a.sig+1)*log(sigmasq.cand) - b.sig/sigmasq.cand 
          -(a.tau+1)*log(tausq.cand) - b.tau/tausq.cand 
          +log(sigmasq.cand) + log(tausq.cand) + log(phi.cand - a.phi) + log(b.phi - phi.cand)
          -0.5*logDetSigma-0.5*(t(Y-X%*%mu.cand)%*%SigmaInv%*%(Y-X%*%mu.cand))
          )
  
  return(out)
}


#starting values
inits <- c(0, log(1000), log(1000), logit(3/10, a.phi, b.phi))

metrop.out <- adaptMetropGibbs(ltd=target, starting=inits, batch=500, batch.length=25, report=100)

summary(mcmc(metrop.out$acceptance))


###################################################
### chunk number 4: accept
###################################################
plot(mcmc(metrop.out$acceptance))


###################################################
### chunk number 5: fullChains
###################################################
p.samples <- metrop.out$p.samples
p.samples[,2] <- exp(metrop.out$p.samples[,2])
p.samples[,3] <- exp(metrop.out$p.samples[,3])
p.samples[,4] <- 3/logit.inv(metrop.out$p.samples[,4], a.phi, b.phi)

colnames(p.samples) <- c("mu", "sigma.sq", "tau.sq", "effective range")

plot(mcmc(p.samples), smooth = FALSE, density=FALSE)


###################################################
### chunk number 6: burninChains
###################################################
##with burn.in
n.samples <- nrow(p.samples)
burn.in <- as.integer(0.25*n.samples)
p.samples <- mcmc(p.samples[burn.in:n.samples,])

summary(p.samples)

plot(p.samples, density=FALSE)