WilcoxonSRCBs <- function(X,alpha=0.05,display=FALSE)
# Suppose that X is a random vector such that for some
# unknown real parameter mu, the vector X' := X - mu
# (with i-th component X[i] - mu) satisfies two conditions:
# - The distribution of X' remains unchanged if we
#   multiply its components with random signs.
# - The absolute values |X[1]|, |X[2]|, ..., |X[n]| are
#   almost surely nonzero and pairwise different.
# (For instance, let the components of X be i.i.d. with
# c.d.f. F_o(. - mu), where F_o is a continuous c.d.f.
# such that F_o(-r) = 1 - F_o(r)  for all real r.)
# Under this assmption, this procedure computes one- and
# two-sided (1 - alpha)-confidence bounds for the
# parameter mu.
# The one-sided bounds a1, b1 are given by W[k1], W[l1],
# and the confidence interval equals [W[k2], W[l2]], where
#    W[1] <= W[2] <= ... <= W[nt]
# are the ordered nt = n(n+1)/2 Walsh-means
#    (X[i] + X[j])/2,  1 <= i <= j <= n ,
# and k1, l1, k2, l2 are indices in {1,2,...,nt} depending
# on n and alpha and satisfying
#    lj = nt+1 - kj . 
# 
# Lutz Duembgen, June 2, 2016.
{
	n <- length(X)
	Tmax <- n*(n+1)/2
	cdf <- rep(1,Tmax+1)
	nt <- 0
	for (i in 1:n)
	{
		nt_new <- nt+i
		cdf[(i+1):nt_new] <-
			(cdf[(i+1):nt_new] + cdf[1:nt])/2
		cdf[1:i] <- cdf[1:i]/2
		nt <- nt_new
	}
	
	k1 <- Tmax - min((0:Tmax)[cdf >= 1 - alpha])
	l1 <- Tmax + 1 - k1
	k2 <- Tmax - min((0:Tmax)[cdf >= 1 - alpha/2])
	l2 <- Tmax + 1 - k2
 
	ww <- rep(0,Tmax)
	for (j in 1:n)
	{
		ww[(j-1)*j/2 + (1:j)] <- (X[1:j] + X[j])/2
	}
	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)
	{
		ww <- rep(ww,each=2)
		Tww <- rep(0,length(ww))
		Tww[seq(1,length(ww)-1,by=2)] <- seq(Tmax,  2-Tmax,by=-2)
		Tww[seq(2,length(ww),  by=2)] <- seq(Tmax-2,-Tmax, by=-2)
		ww <- c(1.05*ww[1] - 0.05*ww[length(ww)],
			ww,
			1.05*ww[length(ww)] - 0.05*ww[1])
		Tww <- c(Tww[1], Tww, Tww[length(Tww)])
		plot(ww, Tww, type="l",
			lwd=2,col="black",
			xlab=expression(italic(m)),
			ylab=expression(italic(T(X-m))),
			main="")
		abline(v=a1,col="magenta")
		abline(v=b1,col="magenta")
		abline(v=a2,col="red")
		abline(v=b2,col="red")
	}
	
	return(list(lower.bound=a1,
		upper.bound=b1,
		conf.int=c(a2,b2),
		k1=k1,l1=l1,k2=k2,l2=l2,nt=nt,
		n=n,alpha=alpha))
}