Description Usage Arguments Value References See Also Examples
ProjSp
produces samples from the posterior distribtion
of the spatial projected normal model.
1 2 3 4 5 6 7 8 9  ProjSp(x = x, coords = coords, start = list(alpha = c(1, 1, 0.5,
0.5), tau = c(0.1, 0.5), rho = c(0.1, 0.5), sigma2 = c(0.1, 0.5), r =
rep(1, times = length(x))), priors = list(tau = c(8, 14), rho = c(8,
14), sigma2 = c(), alpha_mu = c(1, 1), alpha_sigma = c()),
sd_prop = list(sigma2 = 0.5, tau = 0.5, rho = 0.5, sdr = sample(0.05,
length(x), replace = TRUE)), 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), corr_fun = "exponential",
kappa_matern = 0.5, n_chains = 2, parallel = FALSE, n_cores = 1)

x 
a vector of n circular data in [0,2π). If they are not in [0,2π), the function will transform the data in the right interval 
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 vector with

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

sd_prop 
list of 4 elements. To run the MCMC for the rho, 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 
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. 
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
rho
,tau
, 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
corr_fun
characters with the type of spatial correlation chosen
distribution
characters, always "ProjSp"
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
ProjKrigSp
for spatial interpolation under the projected normal model,
WrapSp
for spatial sampling from
Wrapped Normal and WrapKrigSp
for
Kriging estimation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  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 using 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
tau < 0.2
sigma2 < 1
alpha < c(0.5,0.5)
SIGMA < sigma2*exp(rho*Dist)
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[(i1)*2+2],Y[(i1)*2+1])
}
theta < theta %% (2*pi) #to be sure to have values in (0,2pi)
hist(theta)
rose_diag(theta)
val < sample(1:n,round(n*0.1))
################some useful quantities
rho.min < 3/max(Dist)
rho.max < rho.min+0.5
set.seed(100)
mod < ProjSp(
x = theta[val],
coords = coords[val,],
start = list("alpha" = c(0.92, 0.18, 0.56, 0.35),
"rho" = c(0.51,0.15),
"tau" = c(0.46, 0.66),
"sigma2" = c(0.27, 0.3),
"r" = abs(rnorm( length(theta)) )),
priors = list("rho" = c(rho.min,rho.max),
"tau" = c(1,1),
"sigma2" = c(10,3),
"alpha_mu" = c(0, 0),
"alpha_sigma" = diag(10,2)
) ,
sd_prop = list("sigma2" = 0.1, "tau" = 0.1, "rho" = 0.1,
"sdr" = sample(.05,length(theta), replace = TRUE)),
iter = 10000,
BurninThin = c(burnin = 7000, thin = 10),
accept_ratio = 0.234,
adapt_param = c(start = 130000, end = 120000, exp = 0.5),#no adaptation
corr_fun = "exponential",
kappa_matern = .5,
n_chains = 2 ,
parallel = TRUE ,
n_cores = 2
)
# If you don't want to install/use DoParallel
# please set parallel = FALSE. Keep in mind that it can be substantially slower
# How much it takes?
check < ConvCheck(mod)
check$Rhat #close to 1 we have convergence
#### graphical check
par(mfrow=c(3,2))
coda::traceplot(check$mcmc)
par(mfrow=c(1,1))
# once convergence is achieved move to prediction using ProjKrigSp

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.