ProjKrigSp: Kriging using projected normal model.
In CircSpaceTime: Spatial and Spatio-Temporal Bayesian Model for Circular Data

Description Usage Arguments Value References See Also Examples

ProjKrigSp function computes the spatial prediction for circular spatial data using samples from the posterior distribution of the spatial projected normal

1	ProjKrigSp(ProjSp_out, coords_obs, coords_nobs, x_obs)

`ProjSp_out`	the function takes the output of `ProjSp` function
`coords_obs`	coordinates of observed locations (in UTM)
`coords_nobs`	coordinates of unobserved locations (in UTM)
`x_obs`	observed values in [0,2π). If they are not in [0,2π), the function will transform the data in the right interval

a list of 3 elements

M_out: the mean of the associated linear process on the prediction locations coords_nobs (rows) over all the posterior samples (columns) returned by ProjSp
V_out: the variance of the associated linear process on the prediction locations coords_nobs (rows) over all the posterior samples (columns) returned by ProjSp
Prev_out: the posterior predicted values at the unobserved locations.

F. Wang, A. E. Gelfand, "Modeling space and space-time directional data using projected Gaussian processes", Journal of the American Statistical Association,109 (2014), 1565-1580

G. Mastrantonio, G. Jona Lasinio, A. E. Gelfand, "Spatio-temporal circular models with non-separable covariance structure", TEST 25 (2016), 331-350 https://doi.org/10.1007/s11749-015-0458-y

ProjSp for spatial sampling from Projected Normal , WrapSp for spatial sampling from Wrapped Normal and WrapKrigSp for spatial interpolation under the wrapped model

Other spatial interpolations: WrapKrigSp

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[(i-1)*2+2],Y[(i-1)*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))

# move to prediction once convergence is achieved
Krig <- ProjKrigSp(
 ProjSp_out = mod,
 coords_obs =  coords[-val,],
 coords_nobs =  coords[val,],
 x_obs = theta[-val]
)

# The quality of prediction can be checked using APEcirc and CRPScirc
ape  <- APEcirc(theta[val],Krig$Prev_out)
crps <- CRPScirc(theta[val],Krig$Prev_out)