R/rcagep.test.R

In GTRT: Graph Theoretic Randomness Tests

#' RCAG-EP Test
#'
#' Performs the RCAG-EP test of randomness for circular data.
#'
#' @param theta A numeric vector.
#' @param alpha The level of significance 
#' @return Probability of non-intersection of edges, cutoff for RCAG-EP test and adjusted p-values for the RCAG-EP test.
#' @examples
#' x <- arima.sim(model = list(ar=0.9), 1000) ## AR(1) model
#' theta <- ((2*atan(x))%%(2*pi))*(180/pi) ##LAR(1) model
#' rcagep.test(theta,0.05)
#' @export
rcagep.test <- function(theta,alpha) {
  m = length(theta)
  n = m %/% 4
  r = m %% 4
  std = sqrt(5/(36 * n))
  c = qnorm((1-alpha/2), 1/6, std) - 1/6

  adj_p_value <- c()
  p_value <- c()
  p_hat <- c()
  a <- c()

  for (i in 1:(r+1)) {
    s = c()
    t = c()
    for (j in 1:(2*n)) {
      if (2 * j + i - 1 < m) {
        s[j] = theta[(2 * j - 1 + i - 1) %% m]
        t[j] = theta[(2 * j + i - 1) %% m]
      } else if (2 * j + i - 1 == m) {
        s[j] = theta[(2 * j - 1 + i - 1) %% m]
        t[j] = theta[2 * j + i - 1]
      }
    }

    e1 = c()
    e2 = c()
    g <- seq(1, 2 * n)
    p <- c()
    for (j in 1:n) {
      p = sample(g, 2, replace = FALSE)
      e1[j] = p[1]
      e2[j] = p[2]
      g <- g[-which(g %in% p)]
    }

    p_hat[i] = nip.rcag(s, t, e1, e2)
    a[i] = abs(p_hat[i] - 1/6)
    p_value[i] = pnorm(a[i] + 1/6, 1/6, std, lower.tail = FALSE) +
      pnorm(-a[i] + 1/6, 1/6, std, lower.tail = TRUE)
  }

  adj_p_value = p.adjust(p_value, method = "BH")
  structure(list(
    statistic = p_hat,
    cutoff = round(c, 5), 
    p.value = adj_p_value,
    method     = "RCAG-EP Test",
    data.name  = deparse(substitute(theta))
  ),
  class = "ep")
}