WrapSpTi {CircSpaceTime}

Samples from the posterior distribution of the Wrapped Normal spatial temporal model

Description

The WrapSpTi function returns samples from the posterior distribution of the spatio-temporal Wrapped Gaussian Model

Usage

WrapSpTi(x = x, coords = coords, times, start = list(alpha = c(2, 1),
  rho_sp = c(0.1, 0.5), rho_t = c(0.1, 1), sep_par = c(0.01, 0.1), k =
  sample(0, length(x), replace = T)), priors = list(alpha = c(pi, 1, -10,
  10), rho_sp = c(8, 14), rho_t = c(1, 2), sep_par = c(0.001, 1), sigma2 =
  c()), sd_prop = list(rho_sp = 0.5, rho_t = 0.5, sep_par = 0.5, sigma2 =
  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), n_chains = 1, parallel = FALSE, n_cores = 1)

Arguments

x

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

coords

an nx2 matrix with the sites coordinates

times

an n vector with the times of the observations x

start

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

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

• rho_sp the spatial decay parameter,

• rho_t the temporal decay parameter,

• sigma2 the process variance,

• sep_par the separation parameter,

• k the vector of length(x) winding numbers

priors

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

alpha

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

rho_sp

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

rho_t

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

sep_par

a vector of 2 elements defining the two parameters of a beta distribution,

sigma2

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

sd_prop

list of 3 elements. To run the MCMC for the rho_sp 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

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.

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.

Value

it returns a list of n_chains lists each with elements

alpha, rho_sp, rho_t, sep_par, sigma2

vectors with the thinned chains

k

a matrix with nrow = length(x) and ncol = the length of thinned chains

distribution

characters, always "WrapSpTi"

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

References

G. Mastrantonio, G. Jona Lasinio, A. E. Gelfand, "Spatio-temporal circular models with non-separable covariance structure", TEST 25 (2016), 331–350.

T. Gneiting, "Nonseparable, Stationary Covariance Functions for Space-Time Data", JASA 97 (2002), 590-600

See Also

WrapKrigSpTi for spatio-temporal prediction, ProjSpTi to sample from the posterior distribution of the spatio-temporal Projected Normal model and ProjKrigSpTi for spatio-temporal prediction under the same model

Other spatio-temporal models: ProjSpTi

Examples

library(CircSpaceTime)
## functions
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                       ##
######################################
set.seed(1)
n <- 20
### simulate coordinates from a unifrom distribution
coords  <- cbind(runif(n,0,100), runif(n,0,100)) #spatial coordinates
coordsT <- sort(runif(n,0,100)) #time coordinates (ordered)
Dist <- as.matrix(dist(coords))
DistT <- as.matrix(dist(coordsT))

rho     <- 0.05 #spatial decay
rhoT    <- 0.01 #temporal decay
sep_par <- 0.5 #separability parameter
sigma2  <- 0.3 # variance of the process
alpha   <- c(0.5)
#Gneiting covariance
SIGMA <- sigma2 * (rhoT * DistT^2 + 1)^(-1) * exp(-rho * Dist/(rhoT * DistT^2 + 1)^(sep_par/2))

Y <- rmnorm(1,rep(alpha, times = n), SIGMA) #generate the linear variable
theta <- c()
## wrapping step
for(i in 1:n) {
  theta[i] <- Y[i] %% (2*pi)
}
### Add plots of the simulated data

rose_diag(theta)
## use this values as references for the definition of initial values and priors
rho_sp.min <- 3/max(Dist)
rho_sp.max <- rho_sp.min+0.5
rho_t.min  <- 3/max(DistT)
rho_t.max  <- rho_t.min+0.5
val <- sample(1:n,round(n*0.2)) #validation set
set.seed(100)
mod <- WrapSpTi(
  x       = theta[-val],
  coords    = coords[-val,],
  times    = coordsT[-val],
  start   = list("alpha"      = c(.79, .74),
                 "rho_sp"     = c(.33,.52),
                 "rho_t"     = c(.19, .43),
                 "sigma2"    = c(.49, .37),
                 "sep_par"  = c(.47, .56),
                 "k"       = sample(0,length(theta[-val]), replace = TRUE)),
  priors   = list("rho_sp"      = c(0.01,3/4), ### uniform prior on this interval
                  "rho_t"      = c(0.01,3/4), ### uniform prior on this interval
                  "sep_par"  = c(1,1), ### beta prior
                  "sigma2"    = c(5,5),## inverse gamma prior with mode=5/6
                  "alpha" =  c(0,20) ## wrapped gaussian with large variance
  )  ,
  sd_prop   = list( "sigma2" = 0.1,  "rho_sp" = 0.1,  "rho_t" = 0.1,"sep_par"= 0.1),
  iter    = 7000,
  BurninThin    = c(burnin = 3000, thin = 10),
  accept_ratio = 0.234,
  adapt_param = c(start = 1, end = 1000, exp = 0.5),
  n_chains = 2 ,
  parallel = FALSE,
  n_cores = 1
)
check <- ConvCheck(mod,startit = 1 ,thin = 1)
check$Rhat ## convergence has been reached
## when plotting chains remember that alpha is a circular variable
par(mfrow = c(3,2))
coda::traceplot(check$mcmc)
par(mfrow = c(1,1))

#### move to the prediction step with WrapKrigSpTi