R/rigep.test.R

In GTRT: Graph Theoretic Randomness Tests

#' RIG-EP Test
#'
#' Performs the RIG-EP test of randomness.
#'
#' @param x A numeric vector 
#' @param alpha The level of significance 
#' @return Probability of non-intersection of edges, cutoff for RIG-EP test and adjusted p-values for the RIG-EP test.
#' @examples
#' x <- arima.sim(model = list(ar=0.9), 1000) ## AR(1) model
#' rigep.test(x,0.05)
#' @export
rigep.test <- function(x,alpha) {
  m <- length(x)
  n <- m %/% 4
  r <- m %% 4
  std <- sqrt(2/(9*n))
  c <- qnorm((1-alpha/2), 1/3, std) - 1/3
  p_value <- c()
  p_hat <- c()
  a <- c()
  for (i in 1:(r + 1)) {
    s <- t <- c()
    for (j in 1:(2*n)) {
      if (2*j+i-1 < m) {
        s[j] = x[(2*j-1 + i-1)%%m]
        t[j] = x[(2*j + i-1)%%m]
      }
      else if (2*j+i-1 == m){
        s[j] = x[(2*j-1 + i-1)%%m]
        t[j] = x[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)  ## n pair of intervals out of 2n intervals
      e1[j] = p[1]
      e2[j] = p[2]
      g <- g[-which(g %in% p )]
    }
    p_hat[i] = nip.rig(s,t,e1,e2)
    a[i] <- abs(p_hat[i] - 1/3)
    p_value[i] <- pnorm(a[i] + 1/3, 1/3, std, lower.tail = FALSE) +
      pnorm(-a[i] + 1/3, 1/3, 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     = "RIG-EP Test",
    data.name  = deparse(substitute(x))
  ),
  class = "ep")
}