MannWhitneyCBs <- function(X,Y,alpha=0.05, display=FALSE) # MannWhitneyCBs(X,Y,alpha) # computes for two independent samples # X in R^m and Y in R^n confidence bounds # for an unknown shift parameter Delta, # assuming that # X(1), ..., X(m), Y(1), ..., Y(n) # are independent, where # X(i) has distribution function F(. - Delta) , # Y(j) has distribution function F , # with an unknown continuous cdf F. # # In case of min(m,n) > 100, the program uses a # normal approximation for the underlying null # distribution. # # Lutz Duembgen, May 2, 2013 { m <- length(X) n <- length(Y) tmp <- LocalMWU_cdf(m,n) c <- tmp\$c Fc <- tmp\$Fc Tmax <- m*n; k1 <- Tmax - min(c[Fc >= 1 - alpha]) l1 <- Tmax + 1 - k1 k2 <- Tmax - min(c[Fc >= 1 - alpha/2]) l2 <- Tmax + 1 - k2 ww <- rep(0,Tmax) for (i in 1:m) { ww[((i-1)*n+1):(i*n)] <- X[i] - Y } ww <- sort(ww) if (k1 > 0) { a1 <- ww[k1] b1 <- ww[l1] } else { a1 <- -Inf b1 <- Inf } if (k2 > 0) { a2 <- ww[k2] b2 <- ww[l2] } else { a2 <- -Inf b2 <- Inf } if (display) { ww2 <- rep(ww,each=2) Tww2 <- c(Tmax,rep((Tmax-1):1,each=2),0) plot(ww2, Tww2, type="l",lwd=2, xlab=expression(italic(Delta)), ylab=expression(italic(T[U](Delta)))) abline(v=a1, col="magenta") abline(v=b1, col="magenta") abline(v=a2, col="red") abline(v=b2, col="red") } return(list( Delta.est=median(ww), lower.bound=a1, upper.bound=b1, conf.int=c(a2,b2), conf.level=1-alpha)) } LocalMWU_cdf <- function(m,n) # computes the support (-> c) and the # distribution function (-> Fc) # of the random variable # sum_{i=1}^m sum_{j=1}^n 1{Pi(i) > Pi(m+j)} , # Pi a random permutation of {1, 2, ..., m+n} . { CMax <- m*n c <- (0:CMax) if (min(m,n) > 100) { mu <- m*n/2 sigma <- sqrt(m*n*(m+n+1)/12) Fc <- pnorm((c + 0.5 - mu)/sigma) return(list(c=c,Fc=Fc)) } if (m <= n) { reverse <- FALSE } else { tmp <- m m <- n n <- tmp reverse <- TRUE } FM <- matrix(1,CMax+1,m) for (s in 1:m) { FM[1:s,s] <- (1:s)/(s+1) } for (t in 2:n) { FM[1:t,1] = (1:t)/(t+1) for (s in 2:m) { c_max <- s*t FM[1:c_max,s] <- FM[1:c_max,s]*(t/(s+t)) FM[(t+1):c_max,s] <- FM[(t+1):c_max,s] + FM[1:(c_max-t),s-1]*(s/(s+t)) } } Fc = FM[,m] if (reverse == 1) { Fc <- 1 - c(Fc[(length(Fc)-1):1],0) } return(list(c=c,Fc=Fc)) }