WrapSp {CircSpaceTime}R Documentation

Samples from the Wrapped Normal spatial model


The function WrapSp produces samples from the posterior distribution of the wrapped normal spatial model.


WrapSp(x = x, coords = coords, start = list(alpha = c(2, 1), rho =
  c(0.1, 0.5), sigma2 = c(0.1, 0.5), k = sample(0, length(x), replace =
  T)), priors = list(alpha = c(pi, 1, -10, 10), rho = c(8, 14), sigma2 =
  c()), sd_prop = list(sigma2 = 0.5, rho = 0.5), iter = 1000,
  BurninThin = c(burnin = 20, thin = 10), accept_ratio = 0.234,
  adapt_param = c(start = 1, end = 1e+07, exp = 0.9),
  corr_fun = "exponential", kappa_matern = 0.5, n_chains = 1,
  parallel = FALSE, n_cores = 1)



a vector of n circular data in [0,2\pi) If they are not in [0,2\pi), the function will tranform the data in the right interval


an nx2 matrix with the sites coordinates


a list of 4 elements giving initial values for the model parameters. Each elements is a numeric vector with n_chains elements

  • alpha the mean which value is in [0,2\pi).

  • rho the spatial decay parameter

  • sigma2 the process variance

  • k the vector of length(x) winding numbers


a list of 3 elements to define priors for the model parameters:


a vector of 2 elements the mean and the variance of a Wrapped Gaussian distribution, default is mean \pi and variance 1,


a vector of 2 elements defining the minimum and maximum of a uniform distribution,


a vector of 2 elements defining the shape and rate of an inverse-gamma distribution,


list of 3 elements. To run the MCMC for the rho and sigma2 parameters we use an adaptive metropolis and in sd.prop we build a list of initial guesses for these two parameters and the beta parameter


number of iterations


a vector of 2 elements with the burnin and the chain thinning


it is the desired acceptance ratio in the adaptive metropolis


a vector of 3 elements giving the iteration number at which the adaptation must start and end. The third element (exp) must be a number in (0,1) and it is a parameter ruling the speed of changes in the adaptation algorithm, it is recommended to set it close to 1, if it is too small non positive definite matrices may be generated and the program crashes.


characters, the name of the correlation function; currently implemented functions are c("exponential", "matern","gaussian")


numeric, the smoothness parameter of the Matern correlation function, default is kappa_matern = 0.5 (the exponential function)


integer, the number of chains to be launched (default is 1, but we recommend to use at least 2 for model diagnostic)


logical, if the multiple chains must be lunched in parallel (you should install doParallel package). Default is FALSE


integer, required if parallel=TRUE, the number of cores to be used in the implementation. Default value is 1.


It returns a list of n_chains lists each with elements

Implementation Tips

To facilitate the estimations, the observations x are centered around pi, and the prior and starting value of alpha are changed accordingly. After the estimations, posterior samples of alpha are changed back to the original scale


G. Jona Lasinio, A. Gelfand, M. Jona-Lasinio, "Spatial analysis of wave direction data using wrapped Gaussian processes", The Annals of Applied Statistics 6 (2013), 1478-1498

See Also

WrapKrigSp for spatial interpolation, ProjSp for posterior sampling from the Projected Normal model and ProjKrigSp for spatial interpolation under the same model


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

# Simulation with exponential spatial covariance function
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)

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

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
## graphical check
par(mfrow = c(3,1))
par(mfrow = c(1,1))
##### We move to the spatial interpolation see WrapKrigSp

[Package CircSpaceTime version 0.9.0 Index]