ProjKrigSp: Kriging using projected normal model.

In GianlucaMastrantonio/CircSpaceTime: Spatial and Spatio-Temporal Bayesian Model for Circular Data

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π). If they are not in [0,2π), 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

See Also

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)

GianlucaMastrantonio/CircSpaceTime documentation built on July 8, 2020, 2:34 p.m.