R/discrete.DRPT.R

In DRPT: Density Ratio Permutation Test

#' A function implementing the discrete version of the DRPT for discrete data with finite support.
#'
#' A function that implements the discrete version of the DRPT for discrete data with finite support
#' as defined in Section 2.1 in \insertCite{BB2025DRPT;textual}{DRPT}.
#'
#' Counts for the permuted samples are drawn using \code{rMFNCHypergeo} from the package \code{BiasedUrn}.
#' When \code{type="U"} the test statistic is the U-statistic (12); when \code{type="V"} the test statistic is the V-statistic (11); setting \code{type="D"}
#' gives the test statistic (56) in Appendix B of the paper.
#'
#' @param X A numeric vector containing the first sample.
#' @param Y A numeric vector containing the second sample.
#' @param r A numeric vector of positive values specifying the hypothesised density ratio in the discrete setting.
#' @param type A character string indicating the test statistic to use. See the Details section for more information.
#' Defaults to \code{"V"}.
#' @param H An integer specifying the number of permutations to use. Defaults to 99.
#'
#' @return The p-value of the DRPT as defined in (2) in \insertCite{BB2025DRPT;textual}{DRPT}.
#' @export
#'
#' @references \insertRef{BB2025DRPT}{DRPT}
#'
#' @importFrom BiasedUrn rMFNCHypergeo
#'
#' @examples
#' n = 100; m = n
#' X = sample(0:3, n, prob = c(1/8, 1/8, 3/8, 3/8), replace = TRUE)
#' Y = sample(0:3, m, prob = c(1/43, 3/43, 16/43, 23/43), replace = TRUE)
#' r = c(1, 3, 3, 10)
#'
#' discrete.DRPT(X,Y,r,H=19)
#' discrete.DRPT(X,Y,r, type = "U", H=19)
#' discrete.DRPT(X,Y,r, type = "D", H=19)

discrete.DRPT = function(X, Y, r, H=99, type = "V") {
  lev = sort(union(X,Y))
  NX = table(factor(X, levels = lev)); NY = table(factor(Y, levels = lev))
  tot = NX + NY
  n = length(X)
  m = length(Y)
  r = r/r[1]

  # Compute T_hat_0
  T_hat_0 = discreteT(NX, NY, r, n, m, type)

  # Compute T_hat
  sum_indicator = 0

  for (b in 1:H) {
    NY_sigma = rMFNCHypergeo(1, tot, m, r)
    T_hat_b = discreteT(tot - NY_sigma, NY_sigma, r, n, m, type)

    sum_indicator = sum_indicator + (T_hat_b >= T_hat_0)
  }

  p_hat = (1 + sum_indicator) / (H + 1)
  return(p_hat)
}