Description Usage Arguments Value References See Also Examples
APEcirc
computes the average prediction error (APE),
defined as the average circular distance across pairs
1 |
real |
a vector of the values of the process at the test locations |
sim |
a matrix with |
bycol |
logical. It is TRUE if the columns of sim represent the observations and the rows the posterior samples, the default value is FALSE. |
a list of two elements
ApePoints
a vector of APE, one element for each test point
Ape
the overall mean
G. Jona Lasinio, A. Gelfand, M. Jona-Lasinio, "Spatial analysis of wave direction data using wrapped Gaussian processes", The Annals of Applied Statistics 6 (2013), 1478-1498
ProjKrigSp
and WrapKrigSp
for posterior spatial
estimations,
ProjKrigSpTi
and WrapKrigSpTi
for posterior spatio-temporal
estimations
Other model performance indices: CRPScirc
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 84 85 86 87 88 89 90 | library(CircSpaceTime)
## functions
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 ##
######################################
set.seed(1)
n <- 20
### simulate coordinates from a unifrom distribution
coords <- cbind(runif(n,0,100), runif(n,0,100)) #spatial coordinates
coordsT <- sort(runif(n,0,100)) #time coordinates (ordered)
Dist <- as.matrix(dist(coords))
DistT <- as.matrix(dist(coordsT))
rho <- 0.05 #spatial decay
rhoT <- 0.01 #temporal decay
sep_par <- 0.5 #separability parameter
sigma2 <- 0.3 # variance of the process
alpha <- c(0.5)
#Gneiting covariance
SIGMA <- sigma2 * (rhoT * DistT^2 + 1)^(-1) * exp(-rho * Dist/(rhoT * DistT^2 + 1)^(sep_par/2))
Y <- rmnorm(1,rep(alpha, times = n), SIGMA) #generate the linear variable
theta <- c()
## wrapping step
for(i in 1:n) {
theta[i] <- Y[i] %% (2*pi)
}
### Add plots of the simulated data
rose_diag(theta)
## use this values as references for the definition of initial values and priors
rho_sp.min <- 3/max(Dist)
rho_sp.max <- rho_sp.min+0.5
rho_t.min <- 3/max(DistT)
rho_t.max <- rho_t.min+0.5
val <- sample(1:n,round(n*0.2)) #validation set
set.seed(100)
mod <- WrapSpTi(
x = theta[-val],
coords = coords[-val,],
times = coordsT[-val],
start = list("alpha" = c(.79, .74),
"rho_sp" = c(.33,.52),
"rho_t" = c(.19, .43),
"sigma2" = c(.49, .37),
"sep_par" = c(.47, .56),
"k" = sample(0,length(theta[-val]), replace = TRUE)),
priors = list("rho_sp" = c(0.01,3/4), ### uniform prior on this interval
"rho_t" = c(0.01,3/4), ### uniform prior on this interval
"sep_par" = c(1,1), ### beta prior
"sigma2" = c(5,5),## inverse gamma prior with mode=5/6
"alpha" = c(0,20) ## wrapped gaussian with large variance
) ,
sd_prop = list( "sigma2" = 0.1, "rho_sp" = 0.1, "rho_t" = 0.1,"sep_par"= 0.1),
iter = 7000,
BurninThin = c(burnin = 3000, thin = 10),
accept_ratio = 0.234,
adapt_param = c(start = 1, end = 1000, exp = 0.5),
n_chains = 2 ,
parallel = FALSE,
n_cores = 1
)
check <- ConvCheck(mod,startit = 1 ,thin = 1)
check$Rhat ## convergence has been reached
## when plotting chains remember that alpha is a circular variable
par(mfrow = c(3,2))
coda::traceplot(check$mcmc)
par(mfrow = c(1,1))
############## Prediction on the validation set
Krig <- WrapKrigSpTi(
WrapSpTi_out = mod,
coords_obs = coords[-val,],
coords_nobs = coords[val,],
times_obs = coordsT[-val],
times_nobs = coordsT[val],
x_obs = theta[-val]
)
### checking the prediction
Wrap_Ape <- APEcirc(theta[val], Krig$Prev_out)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.