ProjSpTi {CircSpaceTime}  R Documentation 
ProjSpTi
produces samples from the posterior distribution of the spatial
projected normal model.
ProjSpTi(x = x, coords = coords, times = c(), start = list(alpha =
c(1, 1, 0.5, 0.5), tau = c(0.1, 0.5), rho_sp = c(0.1, 0.5), rho_t =
c(0.1, 0.5), sep_par = c(0.1, 0.5), sigma2 = c(0.1, 0.5), r = sample(1,
length(x), replace = T)), priors = list(tau = c(8, 14), rho_sp = c(8,
14), rho_t = c(8, 14), sep_par = c(8, 14), sigma2 = c(), alpha_mu = c(1,
1), alpha_sigma = c()), sd_prop = list(sigma2 = 0.5, tau = 0.5, rho_sp
= 0.5, rho_t = 0.5, sep_par = 0.5, sdr = sample(0.05, length(x), replace
= T)), iter = 1000, BurninThin = c(burnin = 20, thin = 10),
accept_ratio = 0.234, adapt_param = c(start = 1, end = 1e+07, exp =
0.9, sdr_update_iter = 50), n_chains = 2, parallel = FALSE,
n_cores = 1)
x 
a vector of n circular data in 
coords 
an nx2 matrix with the sites coordinates 
times 
an n vector with the times of .... 
start 
a list of 4 elements giving initial values for the model parameters. Each elements is a vector with

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

sd_prop 
=list of 4 elements. To run the MCMC for the rho_sp, tau and sigma2 parameters and r vector we use an adaptive metropolis and in sd_prop we build a list of initial guesses for these three parameters and the r vector 
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 4 elements giving the iteration number at which the adaptation must start and end. The third element (exp) must be a number in (0,1) 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. The last element (sdr_update_iter) must be a positive integer defining every how many iterations there is the update of the sd (vector) of (vector) r. 
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
tau
, rho_sp
, rho_t
, sigma2
vectors with the thinned chains
alpha
a matrix with nrow=2
and ncol=
the length of thinned chains
r
a matrix with nrow=length(x)
and ncol=
the length of thinned chains
G. Mastrantonio, G.Jona Lasinio, A. E. Gelfand, "Spatiotemporal circular models with nonseparable covariance structure", TEST 25 (2016), 331–350.
F. Wang, A. E. Gelfand, "Modeling space and spacetime directional data using projected Gaussian processes", Journal of the American Statistical Association,109 (2014), 15651580
T. Gneiting, "Nonseparable, Stationary Covariance Functions for SpaceTime Data", JASA 97 (2002), 590600
ProjKrigSpTi
for spatiotemporal prediction under the spatiotemporal projected normal model,
WrapSpTi
to sample from the posterior distribution of the spatiotemporal
Wrapped Normal model and WrapKrigSpTi
for spatiotemporal prediction under the
same model
Other spatiotemporal models: WrapSpTi
library(CircSpaceTime)
#### simulated example
## auxiliary 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 using a gneiting covariance function
####
set.seed(1)
n < 20
coords < cbind(runif(n, 0, 100), runif(n, 0, 100))
coordsT < cbind(runif(n, 0, 100))
Dist < as.matrix(dist(coords))
DistT < as.matrix(dist(coordsT))
rho < 0.05
rhoT < 0.01
sep_par < 0.1
sigma2 < 1
alpha < c(0.5)
SIGMA < sigma2 * (rhoT * DistT^2 + 1)^(1) * exp(rho * Dist / (rhoT * DistT^2 + 1)^(sep_par / 2))
tau < 0.2
Y < rmnorm(
1, rep(alpha, times = n),
kronecker(SIGMA, matrix(c(sigma2, sqrt(sigma2) * tau, sqrt(sigma2) * tau, 1), nrow = 2))
)
theta < c()
for (i in 1:n) {
theta[i] < atan2(Y[(i  1) * 2 + 2], Y[(i  1) * 2 + 1])
}
theta < theta %% (2 * pi) ## to be sure we have values in (0,2pi)
rose_diag(theta)
################ some useful quantities
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
### validation set 20% of the data
val < sample(1:n, round(n * 0.2))
set.seed(200)
mod < ProjSpTi(
x = theta[val],
coords = coords[val, ],
times = coordsT[val],
start = list(
"alpha" = c(0.66, 0.38, 0.27, 0.13),
"rho_sp" = c(0.29, 0.33),
"rho_t" = c(0.25, 0.13),
"sep_par" = c(0.56, 0.31),
"tau" = c(0.71, 0.65),
"sigma2" = c(0.47, 0.53),
"r" = abs(rnorm(length(theta[val])))
),
priors = list(
"rho_sp" = c(rho_sp.min, rho_sp.max), # Uniform prior in this interval
"rho_t" = c(rho_t.min, rho_t.max), # Uniform prior in this interval
"sep_par" = c(1, 1), # Beta distribution
"tau" = c(1, 1), ## Uniform prior in this interval
"sigma2" = c(10, 3), # inverse gamma
"alpha_mu" = c(0, 0), ## a vector of 2 elements,
## the means of the bivariate Gaussian distribution
"alpha_sigma" = diag(10, 2) # a 2x2 matrix, the covariance matrix of the
# bivariate Gaussian distribution,
),
sd_prop = list(
"sep_par" = 0.1, "sigma2" = 0.1, "tau" = 0.1, "rho_sp" = 0.1, "rho_t" = 0.1,
"sdr" = sample(.05, length(theta), replace = TRUE)
),
iter = 4000,
BurninThin = c(burnin = 2000, thin = 2),
accept_ratio = 0.234,
adapt_param = c(start = 155000, end = 156000, exp = 0.5),
n_chains = 2,
parallel = TRUE,
)
check < ConvCheck(mod)
check$Rhat ### convergence has been reached when the values are close to 1
#### graphical checking
#### recall that it is made of as many lists as the number of chains and it has elements named
#### as the model's parameters
par(mfrow = c(3, 3))
coda::traceplot(check$mcmc)
par(mfrow = c(1, 1))
# move to prediction once convergence is achieved using ProjKrigSpTi