# Reading the data from the course webpage
d <- read.table(file=url("http://www.maths.usyd.edu.au/stat2912/r/phonelines.txt"), header = T)
d
diff <- d$test-d$control
median(diff)

# Manual calculation of Hodges-Lehmann estimator
l <- 1
# Initializing the vector that will store the averages of all pairs
p <- 1:(length(diff)*(length(diff)+1)/2)
# Filling in the vector with pairwise averages
for (i in 1:length(diff)) for (j in 1:i) {p[l] <- (diff[i]+diff[j])/2; l <- l+1}
median(p)

# Built-in function to calculate Hodges-Lehmann estimator
library(ICSNP)
hl.loc(diff)


# CI
M <- length(diff)*(length(diff)+1)/2
# for(k in 1:25) print(c(k, psignrank(M-k, length(diff))- psignrank(k-1, length(diff))))
# A better version (thanks go to Katherine Huang)
k = qsignrank(0.025, length(diff))
sort(p)[k]
sort(p)[M-k+1]