hodges_lehmann_location_shift_estimate.R

hodges_lehmann_location_shift_estimate <- function (x, y, conf.level) {
  nx <- length(x)
  ny <- length(y)
  # https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_npar1way_a0000000201.htm
  alpha <- 1 - conf.level
  # qnorm in R is the https://en.wikipedia.org/wiki/Probit
  # this is also using the https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test#Normal_approximation_and_tie_correction
  Ca <- floor(nx * ny / 2 - qnorm(1 - alpha / 2) * sqrt(nx * ny * (nx + ny + 1) / 12))
  m <- nx * ny
  U <- vector(mode = "numeric", length = m)
  for (i in 1:nx) {
    for (j in 1:ny) {
      U[(i-1) * ny + j] = x[i] - y[j]
    }
  }
  U <- sort(U)

  lower <- U[Ca]
  median <- if (m %% 2 == 0) (U[m/2] + U[m/2 + 1])/2 else U[(m+1)/2]
  upper <- U[m + 1 - Ca]
  cat("Hodges-Lehmann Estimation of Location Shift\n\n")
  cat("median of all paired differences between observations in the two samples\n")
  cat(paste0("\t", median, "\n"))
  cat(paste0("confidence interval (at least ", 100 * conf.level, " percent)\n"))
  cat(paste0("\t", lower, " ", upper, "\n"))
}

x <- c(330.918433554412, 229.000245801897, 329.965804711059, 300.691712977023, 234.825814560758, 527.963570240824, 308.480228109542, 522.589636888846, 325.180019314907, 259.507753259425, 279.822702423339, 266.629000675101, 348.343075741739, 349.334465676962, 287.134434273056, 255.173593243992, 361.932565928192, 292.828093580835, 315.670479115113, 233.148959408031)
y <- c(311.307077576214, 287.01781844844, 326.260035401324, 354.377742876501, 288.289637886991, 273.101023005202, 274.56113018375, 315.970477279932, 219.797595818019, 347.564349614263, 273.942037169113, 346.293931534072, 245.535433540729, 289.58600732826, 283.287591042244, 329.580201801935, 182.369592581587, 256.401674519353, 260.634916843668, 274.456816887721, 294.182095864991, 273.935688437867, 248.996033572172, 322.630738439236, 250.809890805682, 229.394540671313, 326.040872321786, 235.165610985367, 239.857271090047, 347.272576200616)
conf.level <- 0.95

print(wilcox.test(x, y, conf.int=T, conf.level=conf.level))
hodges_lehmann_location_shift_estimate(x,y,conf.level)
# median of all paired differences between observations in the two samples
# 	19.0410031804765
# confidence interval (at least 95 percent)
# 	-8.46693546365196 51.968823287689