ConvCheck function

Testing Convergence of mcmc using package coda

ConvCheck returns an mcmc.list (mcmc) to be used with the coda package and the Potential scale reduction factors (Rhat) of the model parameters computed using the gelman.diag function in the coda package

ConvCheck(mod, startit = 15000, thin = 10)

Arguments

• mod: is a list with $m\ge 1$ elements, one for each chain generated using WrapSp, ProjSp, WrapSpTi or ProjSpTi
• startit: is an integer, the iteration at which the chains start (required to build the mcmc.list)
• thin: is an integer, the thinning applied to chains

Returns

a list of two elements

• mcmc: an mcmc.list (mcmc) to be used with the coda package
• Rhat: the Potential scale reduction factors of the model parameters computed using the gelman.diag function in the coda package

Examples

library(CircSpaceTime)
## auxiliary function
rmnorm<-function(n = 1, mean = rep(0, d), varcov){
  d <- if (is.matrix(varcov))
    ncol(varcov)
  else 1
  z <- matrix(rnorm(n * d), n, d) %*% chol(varcov)
  y <- t(mean + t(z))
  return(y)
}

####
# Simulation with exponential spatial covariance function
####
set.seed(1)
n <- 20
coords <- cbind(runif(n,0,100), runif(n,0,100))
Dist <- as.matrix(dist(coords))

rho     <- 0.05
sigma2  <- 0.3
alpha   <- c(0.5)
SIGMA   <- sigma2*exp(-rho*Dist)

Y <- rmnorm(1,rep(alpha,times=n), SIGMA)
theta <- c()
for(i in 1:n) {
  theta[i] <- Y[i]%%(2*pi)
}
rose_diag(theta)

#validation set
val <- sample(1:n,round(n*0.1))

set.seed(12345)
mod <- WrapSp(
  x       = theta[-val],
  coords    = coords[-val,],
  start   = list("alpha"      = c(.36,0.38),
                 "rho"     = c(0.041,0.052),
                 "sigma2"    = c(0.24,0.32),
                 "k"       = rep(0,(n - length(val)))),
  priors   = list("rho"      = c(0.04,0.08), #few observations require to be more informative
                  "sigma2"    = c(2,1),
                  "alpha" =  c(0,10)
  ),
  sd_prop   = list( "sigma2" = 0.1,  "rho" = 0.1),
  iter    = 1000,
  BurninThin    = c(burnin = 500, thin = 5),
  accept_ratio = 0.234,
  adapt_param = c(start = 40000, end = 45000, exp = 0.5),
  corr_fun = "exponential",
  kappa_matern = .5,
  parallel = FALSE,
  #With doParallel, bigger iter (normally around 1e6) and n_cores>=2 it is a lot faster
  n_chains = 2 ,
  n_cores = 1
)
check <- ConvCheck(mod)
check$Rhat ## close to 1 means convergence has been reached

See Also

ProjKrigSp and WrapKrigSp for posterior spatial estimations, and ProjKrigSpTi and WrapKrigSpTi for posterior spatio-temporal estimations