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