ProjKrigSp {CircSpaceTime} R Documentation

## Kriging using projected normal model.

### Description

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

### Usage

ProjKrigSp(ProjSp_out, coords_obs, coords_nobs, x_obs)


### Arguments

 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\pi). If they are not in [0,2\pi), the function will transform the data in the right interval

### Value

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.

### References

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

### Examples

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)


[Package CircSpaceTime version 0.9.0 Index]