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))
}