# ProjSpTi: Samples from the posterior distribution of the Projected... In CircSpaceTime: Spatial and Spatio-Temporal Bayesian Model for Circular Data

## Description

`ProjSpTi` produces samples from the posterior distribution of the spatial projected normal model.

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```ProjSpTi(x = x, coords = coords, times = c(), start = list(alpha = c(1, 1, 0.5, 0.5), tau = c(0.1, 0.5), rho_sp = c(0.1, 0.5), rho_t = c(0.1, 0.5), sep_par = c(0.1, 0.5), sigma2 = c(0.1, 0.5), r = sample(1, length(x), replace = T)), priors = list(tau = c(8, 14), rho_sp = c(8, 14), rho_t = c(8, 14), sep_par = c(8, 14), sigma2 = c(), alpha_mu = c(1, 1), alpha_sigma = c()), sd_prop = list(sigma2 = 0.5, tau = 0.5, rho_sp = 0.5, rho_t = 0.5, sep_par = 0.5, sdr = sample(0.05, length(x), replace = T)), iter = 1000, BurninThin = c(burnin = 20, thin = 10), accept_ratio = 0.234, adapt_param = c(start = 1, end = 1e+07, exp = 0.9, sdr_update_iter = 50), n_chains = 2, parallel = FALSE, n_cores = 1) ```

## Arguments

 `x` a vector of n circular data in [0,2π) If they are not in [0,2π), the function will tranform the data in the right interval `coords` an nx2 matrix with the sites coordinates `times` an n vector with the times of .... `start` a list of 4 elements giving initial values for the model parameters. Each elements is a vector with `n_chains` elements alpha the 2-d vector of the means of the Gaussian bi-variate distribution, tau the correlation of the two components of the linear representation, rho_sp the spatial decay parameter, rho_t the temporal decay parameter, sigma2 the process variance, sep_par the separation parameter, r the vector of `length(x)`, latent variable `priors` a list of 7 elements to define priors for the model parameters: alpha_mua vector of 2 elements, the means of the bivariate Gaussian distribution, alpha_sigmaa 2x2 matrix, the covariance matrix of the bivariate Gaussian distribution, rho_spa vector of 2 elements defining the minimum and maximum of a uniform distribution, rho_ta vector of 2 elements defining the minimum and maximum of a uniform distribution, tau vector of 2 elements defining the minimum and maximum of a uniform distribution with the constraints min(tau) >= -1 and max(tau) <= 1, sep_para vector of 2 elements defining the two parameters of a beta distribution, sigma2a vector of 2 elements defining the shape and rate of an inverse-gamma distribution, `sd_prop` =list of 4 elements. To run the MCMC for the rho_sp, tau and sigma2 parameters and r vector we use an adaptive metropolis and in sd_prop we build a list of initial guesses for these three parameters and the r vector `iter` iter number of iterations `BurninThin` a vector of 2 elements with the burnin and the chain thinning `accept_ratio` it is the desired acceptance ratio in the adaptive metropolis `adapt_param` a vector of 4 elements giving the iteration number at which the adaptation must start and end. The third element (exp) must be a number in (0,1) is a parameter ruling the speed of changes in the adaptation algorithm, it is recommended to set it close to 1, if it is too small non positive definite matrices may be generated and the program crashes. The last element (sdr_update_iter) must be a positive integer defining every how many iterations there is the update of the sd (vector) of (vector) r. `n_chains` integer, the number of chains to be launched (default is 1, but we recommend to use at least 2 for model diagnostic) `parallel` logical, if the multiple chains must be lunched in parallel (you should install doParallel package). Default is FALSE `n_cores` integer, required if parallel=TRUE, the number of cores to be used in the implementation. Default value is 1.

## Value

it returns a list of `n_chains` lists each with elements

`tau`, `rho_sp`, `rho_t`, `sigma2`

vectors with the thinned chains

`alpha`

a matrix with `nrow=2` and `ncol=` the length of thinned chains

`r`

a matrix with `nrow=length(x)` and `ncol=` the length of thinned chains

## References

G. Mastrantonio, G.Jona Lasinio, A. E. Gelfand, "Spatio-temporal circular models with non-separable covariance structure", TEST 25 (2016), 331–350.

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

T. Gneiting, "Nonseparable, Stationary Covariance Functions for Space-Time Data", JASA 97 (2002), 590-600

`ProjKrigSpTi` for spatio-temporal prediction under the spatio-temporal projected normal model, `WrapSpTi` to sample from the posterior distribution of the spatio-temporal Wrapped Normal model and `WrapKrigSpTi` for spatio-temporal prediction under the same model
Other spatio-temporal models: `WrapSpTi`
 ``` 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 91 92 93 94 95 96``` ```library(CircSpaceTime) #### simulated example ## auxiliary 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 using a gneiting covariance function #### set.seed(1) n <- 20 coords <- cbind(runif(n, 0, 100), runif(n, 0, 100)) coordsT <- cbind(runif(n, 0, 100)) Dist <- as.matrix(dist(coords)) DistT <- as.matrix(dist(coordsT)) rho <- 0.05 rhoT <- 0.01 sep_par <- 0.1 sigma2 <- 1 alpha <- c(0.5) SIGMA <- sigma2 * (rhoT * DistT^2 + 1)^(-1) * exp(-rho * Dist / (rhoT * DistT^2 + 1)^(sep_par / 2)) tau <- 0.2 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 we have values in (0,2pi) rose_diag(theta) ################ some useful quantities 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 ### validation set 20% of the data val <- sample(1:n, round(n * 0.2)) set.seed(200) mod <- ProjSpTi( x = theta[-val], coords = coords[-val, ], times = coordsT[-val], start = list( "alpha" = c(0.66, 0.38, 0.27, 0.13), "rho_sp" = c(0.29, 0.33), "rho_t" = c(0.25, 0.13), "sep_par" = c(0.56, 0.31), "tau" = c(0.71, 0.65), "sigma2" = c(0.47, 0.53), "r" = abs(rnorm(length(theta[-val]))) ), priors = list( "rho_sp" = c(rho_sp.min, rho_sp.max), # Uniform prior in this interval "rho_t" = c(rho_t.min, rho_t.max), # Uniform prior in this interval "sep_par" = c(1, 1), # Beta distribution "tau" = c(-1, 1), ## Uniform prior in this interval "sigma2" = c(10, 3), # inverse gamma "alpha_mu" = c(0, 0), ## a vector of 2 elements, ## the means of the bivariate Gaussian distribution "alpha_sigma" = diag(10, 2) # a 2x2 matrix, the covariance matrix of the # bivariate Gaussian distribution, ), sd_prop = list( "sep_par" = 0.1, "sigma2" = 0.1, "tau" = 0.1, "rho_sp" = 0.1, "rho_t" = 0.1, "sdr" = sample(.05, length(theta), replace = TRUE) ), iter = 4000, BurninThin = c(burnin = 2000, thin = 2), accept_ratio = 0.234, adapt_param = c(start = 155000, end = 156000, exp = 0.5), n_chains = 2, parallel = TRUE, ) check <- ConvCheck(mod) check\$Rhat ### convergence has been reached when the values are close to 1 #### graphical checking #### recall that it is made of as many lists as the number of chains and it has elements named #### as the model's parameters par(mfrow = c(3, 3)) coda::traceplot(check\$mcmc) par(mfrow = c(1, 1)) # move to prediction once convergence is achieved using ProjKrigSpTi ```