R/dptest.R

In Compositional: Compositional Data Analysis

dptest <- function(x1, x2, B = 100) {
  d <- dim(x1)[2]
  u <- Rfast::matrnorm(B, d)
  u / sqrt( Rfast::rowsums(u^2) )
  x1 <- sqrt(x1)  ;   x2 <- sqrt(x2)
  z1 <- tcrossprod(x1, u)
  z2 <- tcrossprod(x2, u)
  p <- numeric(B)
  for (i in 1:B)
  p[i] <- ks.test(z1[, i], z2[, i])$p.value
  p <- sort(p)
  i <- 1:B
  pval <- min(B / i * p )
  list(pvalues = p, pvalue = pval)
}