# APEcirc: Average Prediction Error for circular Variables. In CircSpaceTime: Spatial and Spatio-Temporal Bayesian Model for Circular Data

## Description

`APEcirc` computes the average prediction error (APE), defined as the average circular distance across pairs

## Usage

 `1` ```APEcirc(real, sim, bycol = F) ```

## Arguments

 `real` a vector of the values of the process at the test locations `sim` a matrix with `nrow =` the test locations and `ncol =` the number of posterior samples from the posterior distributions by `WrapKrigSp` `WrapKrigSpTi`, `ProjKrigSp`, `ProjKrigSpTi` `bycol` logical. It is TRUE if the columns of sim represent the observations and the rows the posterior samples, the default value is FALSE.

## Value

a list of two elements

`ApePoints`

a vector of APE, one element for each test point

`Ape`

the overall mean

## References

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) ```