Description Usage Arguments Value Implementation Tips References See Also Examples
WrapKrigSp
function computes the spatial prediction
for circular spatial data using samples from the posterior distribution
of the spatial wrapped normal
1 | WrapKrigSp(WrapSp_out, coords_obs, coords_nobs, x_obs)
|
WrapSp_out |
the functions takes the output of |
coords_obs |
coordinates of observed locations (in UTM) |
coords_nobs |
coordinates of unobserved locations (in UTM) |
x_obs |
observed values |
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 WrapSp
V_out
the variance of the associated linear process on the prediction locations coords_nobs (rows) over all the posterior samples (columns) returned by WrapSp
Prev_out
the posterior predicted values at the unobserved locations.
To facilitate the estimations, the observations x are centered around pi, and the posterior samples of x and posterior mean are changed back to the original scale
G. Jona-Lasinio, A .E. Gelfand, M. Jona-Lasinio, "Spatial analysis of wave direction data using wrapped Gaussian processes", The Annals of Applied Statistics, 6 (2012), 1478-1498
WrapSp
for spatial sampling from
Wrapped Normal ,
ProjSp
for spatial sampling from
Projected Normal and ProjKrigSp
for
Kriging estimation
Other spatial interpolations: ProjKrigSp
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 | 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 with 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
sigma2 <- 0.3
alpha <- c(0.5)
SIGMA <- sigma2*exp(-rho*Dist)
Y <- rmnorm(1,rep(alpha,times=n), SIGMA)
theta <- c()
for(i in 1:n) {
theta[i] <- Y[i]%%(2*pi)
}
rose_diag(theta)
#validation set
val <- sample(1:n,round(n*0.1))
set.seed(12345)
mod <- WrapSp(
x = theta[-val],
coords = coords[-val,],
start = list("alpha" = c(.36,0.38),
"rho" = c(0.041,0.052),
"sigma2" = c(0.24,0.32),
"k" = rep(0,(n - length(val)))),
priors = list("rho" = c(0.04,0.08), #few observations require to be more informative
"sigma2" = c(2,1),
"alpha" = c(0,10)
),
sd_prop = list( "sigma2" = 0.1, "rho" = 0.1),
iter = 1000,
BurninThin = c(burnin = 500, thin = 5),
accept_ratio = 0.234,
adapt_param = c(start = 40000, end = 45000, exp = 0.5),
corr_fun = "exponential",
kappa_matern = .5,
parallel = FALSE,
#With doParallel, bigger iter (normally around 1e6) and n_cores>=2 it is a lot faster
n_chains = 2 ,
n_cores = 1
)
check <- ConvCheck(mod)
check$Rhat ## close to 1 means convergence has been reached
## graphical check
par(mfrow = c(3,1))
coda::traceplot(check$mcmc)
par(mfrow = c(1,1))
##### We move to the spatial interpolation
Krig <- WrapKrigSp(
WrapSp_out = mod,
coords_obs = coords[-val,],
coords_nobs = coords[val,],
x_obs = theta[-val]
)
#### check the quality of the prediction using APE and CRPS
ApeCheck <- APEcirc(theta[val],Krig$Prev_out)
CrpsCheck <- CRPScirc(theta[val],Krig$Prev_out)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.