WilcoxonSRPVs <- function(X)
# WilcoxonSRPVs(X)
# computes for a data vector X (vector of differences)
# Wilcoxon's signed-rank test statistic
#    T  =  sum_{i=1}^n sign(X[i]) R[i] ,
# where R(i) denotes the rank of |X[i]| among the non-zero 
# components of (|X[1]|, |X[2]|, ..., |X[n]|).
# 
# In case of X[i] = 0 we set R[i] = 0, i.e. components
# equal to zero are ignored.
# 
# In case of ties among the absolute values of X's
# non-zero components, we use the usual convention of
# average ranks.
#
# In addition, the left-, right- and two-sided P-values
#    pvl  =  Pr(T_* <= T) ,
#    pvr  =  Pr(T_* >= T) ,
#    pvz  =  2 * min(pvl, pvr)
# are computed. Here T_* = sum_{i=1}^n S[i]*R[i] with
# independent random signs S[1], S[2], ..., S[n], while
# R is viewed temporarily as a fixed vector.
# 
# Additional output arguments are
# - supp : a vector containing all support points of the
#          distribution function F(.) := Pr(T_* <= .),
# - cdf : a vector containing the numbers cdf(i) = F(supp[i]).
# 
# Lutz Duembgen, May 20, 2011
{
	R <- abs(X)
	tmp <- (R > 0)
	R[tmp] <- rank(R[tmp])
	T <- sum(sign(X)*R)
	
	N  <- sum(tmp)
	R2 <- ceiling(2*sort(R[tmp]))
	h2 <- N*(N+1)
	h1 <- ceiling(h2/2)
	# Now we work temporarily with the test statistic
	#    T2  =  T + N(N+1)/2
	#        =  sum_{i=1}^N 1{X[i] > 0} 2R[i]
	# taking values in {0,1,...,N(N+1)}:
	cdf <- rep(1,h2+1)
	m <- 1
	for (i in 1:N)
	{
		m_new <- m+R2[i]
		cdf[(R2[i]+1):m_new] <- (cdf[1:m] + cdf[(R2[i]+1):m_new])/2
		cdf[1:R2[i]] <- cdf[1:R2[i]]/2
		m <- m_new
	}
	
	supp <- (0:h2) - h1
	
	ind <- h1 + T + 1
	pv.left <- cdf[ind]
	if (ind >= 2)
	{
		pv.right <- 1 - cdf[ind-1]
	}
	else
	{
		pv.right <- 1
	}
	pv.twosided <- 2*min(pv.left,pv.right)
	return(list(T=T,
		pv.left=pv.left,pv.right=pv.right,pv.twosided=pv.twosided,
		supp=supp,cdf=cdf))
}