ConvCheck: Testing Convergence of mcmc using package coda
In CircSpaceTime: Spatial and Spatio-Temporal Bayesian Model for Circular Data

Description Usage Arguments Value See Also Examples

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

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

`mod`	is a list with m≥ 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

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

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

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