R/RcppExports.R

Defines functions murakami_stat_perms murakami_stat bws_cdf bws_stat

Documented in bws_cdf bws_stat murakami_stat murakami_stat_perms

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' @title
#' Compute the test statistic of the Baumgartner-Weiss-Schindler test.
#'
#' @description
#'
#' Compute the Baumgartner-Weiss-Schindler test statistic.
#'
#' @details
#'
#' Given vectors \eqn{X} and \eqn{Y}, computes \eqn{B_X} and \eqn{B_Y} as
#' described by Baumgartner \emph{et al.}, returning their average, \eqn{B}.
#' The test statistic approximates the variance-weighted square norm of the
#' difference in CDFs of the two distributions. For sufficiently large sample
#' sizes (more than 20, say), under the null the test statistic approaches the asymptotic
#' value computed in \code{\link{bws_cdf}}.
#'
#' The test value is an approximation of
#' \deqn{\tilde{B} = \frac{mn}{m+n} \int_0^1 \frac{1}{z(1-z)} \left(F_X(z) - F_Y(z)\right)^2 \mathrm{dz},}
#' where \eqn{m} (\eqn{n}) is the number of elements in \eqn{X} (\eqn{Y}), and
#' \eqn{F_X(z)}{F_X(z)} (\eqn{F_Y(z)}{F_Y(z)}) is the CDF of \eqn{X} (\eqn{Y}).
#'
#' The test statistic is based only on the ranks of the input. If the same
#' monotonic transform is applied to both vectors, the result should be unchanged.
#' Moreover, the test is inherently two-sided, so swapping \eqn{X} and \eqn{Y}
#' should also leave the test statistic unchanged.
#'
#' @param x a vector.
#' @param y a vector.
#'
#' @return The BWS test statistic, \eqn{B}.
#' @seealso \code{\link{bws_cdf}}, \code{\link{bws_test}}
#' @examples
#'
#' set.seed(1234)
#' x <- runif(1000)
#' y <- runif(100)
#' bval <- bws_stat(x,y)
#' # check a monotonic transform:
#' ftrans <- function(x) { log(1 + x) }
#' bval2 <- bws_stat(ftrans(x),ftrans(y))
#' stopifnot(all.equal(bval,bval2))
#' # check commutivity
#' bval3 <- bws_stat(y,x)
#' stopifnot(all.equal(bval,bval3))
#'
#' @template etc
#' @template ref-bws
#' @rdname bws_stat
#' @export
bws_stat <- function(x, y) {
    .Call('_BWStest_bws_stat', PACKAGE = 'BWStest', x, y)
}

#' @title
#' CDF of the Baumgartner-Weiss-Schindler test under the null.
#'
#' @description
#'
#' Computes the CDF of the Baumgartner-Weiss-Schindler test statistic under the
#' null hypothesis of equal distributions.
#'
#' @details
#'
#' Given value \eqn{b}, computes the CDF of the BWS statistic under
#' the null, denoted as \eqn{\Psi(b)}{Psi(b)} by 
#' Baumgartner \emph{et al.}  The CDF is computed from 
#' equation (2.5) via numerical quadrature.
#'
#' The expression for the CDF contains the integral
#' \deqn{\int_0^1 \frac{1}{\sqrt{r^3 (1-r)}} \mathrm{exp}\left(\frac{rb}{8} - \frac{\pi^2 (4j+1)^2}{8rb}\right) \mathrm{dr}}
#' By making the change of variables \eqn{x = 2r - 1}, this can
#' be re-expressed as an integral of the form
#' \deqn{\int_{-1}^1 \frac{1}{\sqrt{1-x^2}} f(x) \mathrm{dx},}
#' for some function \eqn{f(x)} involving \eqn{b} and \eqn{j}. 
#' This integral can be approximated
#' via Gaussian quadrature using Chebyshev nodes (of the first kind), which
#' is the approach we take here.
#'
#' @param b a vector of BWS test statistics.
#' @param maxj the maximum value of j to take in the approximate computation
#' of the CDF via equation (2.5). Baumgartner \emph{et. al.} claim that a
#' value of 3 is sufficient.
#' @param lower_tail boolean, when \code{TRUE} returns \eqn{\Psi}{Psi}, otherwise
#' compute the upper tail, \eqn{1-\Psi}{1 - Psi}, which is more useful for hypothesis tests.
#'
#' @return A vector of the CDF of \eqn{b}, \eqn{\Psi(b)}{Psi(b)}.
#' @seealso \code{\link{bws_stat}}, \code{\link{bws_test}}
#' @examples
#'
#' # do it 500 times
#' set.seed(123)
#' bvals <- replicate(500, bws_stat(rnorm(50),rnorm(50)))
#' pvals <- bws_cdf(bvals)
#' # these should be uniform!
#' \donttest{ 
#'   plot(ecdf(pvals)) 
#' }
#' 
#' # compare to Table 1 of Baumgartner et al.
#' bvals <- c(1.933,2.493,3.076,3.880,4.500,5.990)
#' tab1v <- c(0.9,0.95,0.975,0.990,0.995,0.999)
#' pvals <- bws_cdf(bvals,lower_tail=TRUE)
#' show(data.frame(B=bvals,BWS_psi=tab1v,our_psi=pvals))
#'
#' @template etc
#' @template ref-bws
#' @rdname bws_cdf
#' @export
bws_cdf <- function(b, maxj = 5L, lower_tail = TRUE) {
    .Call('_BWStest_bws_cdf', PACKAGE = 'BWStest', b, maxj, lower_tail)
}

#' @title
#' Compute Murakami's test statistic.
#'
#' @description
#'
#' Compute one of the modified Baumgartner-Weiss-Schindler test statistics proposed
#' by Murakami, or Neuhauser.
#'
#' @details
#'
#' Given vectors \eqn{X} and \eqn{Y}, computes \eqn{B_{jX}} and \eqn{B_{jY}} 
#' for some \eqn{j} as described by Murakami and by Neuhauser, returning either their 
#' their average or their average distance.
#' The test statistics approximate the weighted square norm of the
#' difference in CDFs of the two distributions. 
#'
#' The test statistic is based only on the ranks of the input. If the same
#' monotonic transform is applied to both vectors, the result should be unchanged.
#'
#' The various \sQuote{flavor}s of test statistic are:
#' \describe{
#' \item{0}{The statistic of Baumgartner-Weiss-Schindler.}
#' \item{1}{Murakami's \eqn{B_1} statistic, from his 2006 paper.}
#' \item{2}{Neuhauser's difference statistic, denoted by Murakami as \eqn{B_2} in his 
#' 2012 paper.}
#' \item{3}{Murakami's \eqn{B_3} statistic, from his 2012 paper.}
#' \item{4}{Murakami's \eqn{B_4} statistic, from his 2012 paper.}
#' \item{5}{Murakami's \eqn{B_5} statistic, from his 2012 paper, with a log weighting.}
#' }
#'
#' @param x a vector of the first sample.
#' @param y a vector of the second sample.
#' @param flavor which \sQuote{flavor} of test statistic. 
#' @param nx the length of \code{x}, the first sample.
#' @param ny the length of \code{y}, the second sample.
#'
#' @return The BWS test statistic, \eqn{B_j}. For \code{murakami_stat_perms}, a vector of
#' the test statistics for \emph{all} permutations of the input.
#' @note \code{NA} and \code{NaN} are not yet dealt with!
#' @seealso \code{\link{bws_stat}}.
#' @examples
#'
#' set.seed(1234)
#' x <- runif(1000)
#' y <- runif(100)
#' bval <- murakami_stat(x,y,1)
#'
#' \donttest{
#' nx <- 6
#' ny <- 5
#' # monte carlo
#' set.seed(1234)
#' repli <- replicate(3000,murakami_stat(rnorm(nx),rnorm(ny),0L))
#' # under the null, perform the permutation test:
#' allem <- murakami_stat_perms(nx,ny,0L)
#' plot(ecdf(allem)) 
#' lines(ecdf(repli),col='red') 
#' }
#'
#' @template etc
#' @template ref-bws
#' @template ref-modtests
#' @rdname murakami_stat
#' @export
murakami_stat <- function(x, y, flavor = 0L) {
    .Call('_BWStest_murakami_stat', PACKAGE = 'BWStest', x, y, flavor)
}

#' @rdname murakami_stat
#' @export
murakami_stat_perms <- function(nx, ny, flavor = 0L) {
    .Call('_BWStest_murakami_stat_perms', PACKAGE = 'BWStest', nx, ny, flavor)
}

Try the BWStest package in your browser

Any scripts or data that you put into this service are public.

BWStest documentation built on Oct. 11, 2023, 1:07 a.m.