R/davies.R

In snpsettest: A Set-Based Association Test using GWAS Summary Statistics

##' A `davies` function from `CompQuadForm` package (version 1.4.3)
##'
##' Computes the distrbiution of quadratic forms in normal variables using
##' Davies' algorithm (Davies 1980).
##'
##' @noRd
davies <- function(q, lambda, h = rep(1, length(lambda)),
                   delta = rep(0, length(lambda)),
                   sigma = 0, lim = 10000, acc = 0.0001) {
  r <- length(lambda)
  if (any(delta < 0)) {
    stop("All non centrality parameters in 'delta' should be positive!")
  }
  if (length(h) != r) {
    stop("lambda and h should have the same length!")
  }
  if (length(delta) != r) {
    stop("lambda and delta should have the same length!")
  }

  out <- .C("qfc", lambdas = as.double(lambda), noncentral = as.double(delta),
            df = as.integer(h), r = as.integer(r), sigma = as.double(sigma),
            q = as.double(q), lim = as.integer(lim), acc = as.double(acc),
            trace = as.double(rep(0, 7)), ifault = as.integer(0),
            res = as.double(0), PACKAGE = "snpsettest")

  out$res <- 1 - out$res

  ## Warning unmeaningful results by additionally checking out$res <= 0
  if (out$res > 1 || out$res <= 0) {
    warning("Consider playing with 'lim' or 'acc'.")
  }

  return(list(trace = out$trace, ifault = out$ifault, Qq = out$res))
}