WrapSp {CircSpaceTime}  R Documentation 
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)
x 
a vector of n circular data in 
coords 
an nx2 matrix with the sites coordinates 
start 
a list of 4 elements giving initial values for the model parameters. Each elements is a numeric vector with

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

sd_prop 
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 
iter 
number of iterations 
BurninThin 
a vector of 2 elements with the burnin and the chain thinning 
accept_ratio 
it is the desired acceptance ratio in the adaptive metropolis 
adapt_param 
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. 
corr_fun 
characters, the name of the correlation function; currently implemented functions are c("exponential", "matern","gaussian") 
kappa_matern 
numeric, the smoothness parameter of the Matern
correlation function, default is 
n_chains 
integer, the number of chains to be launched (default is 1, but we recommend to use at least 2 for model diagnostic) 
parallel 
logical, if the multiple chains must be lunched in parallel (you should install doParallel package). Default is FALSE 
n_cores 
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
alpha
, rho
,sigma2
vectors with the thinned chains,
k
a matrix with nrow = length(x)
and ncol =
the length of thinned chains
corr_fun
characters with the type of spatial correlation chosen.
distribution
characters, always "WrapSp"
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. JonaLasinio, "Spatial analysis of wave direction data using wrapped Gaussian processes", The Annals of Applied Statistics 6 (2013), 14781498
WrapKrigSp
for spatial interpolation,
ProjSp
for posterior sampling from the
Projected Normal model and ProjKrigSp
for
spatial interpolation under the same model
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
## graphical check
par(mfrow = c(3,1))
coda::traceplot(check$mcmc)
par(mfrow = c(1,1))
##### We move to the spatial interpolation see WrapKrigSp