R/hsmode.R

In twolodzko/twextras: twextras

#' Half-sample mode
#' 
#' @param x numeric vector.
#' 
#' @references 
#' 
#' Bickel, D.R. (2002). Robust estimators of the mode and skewness of continuous data.
#' Computational Statistics & Data Analysis, 39, 153-163.
#' 
#' Bickel, D.R. and R. Frühwirth. (2006). On a fast, robust estimator of the mode: comparisons to other estimators with applications.
#' Computational Statistics & Data Analysis, 50, 3500-3530.
#' 
#' Dalenius, T. (1965). The mode - A neglected statistical parameter.
#' Journal, Royal Statistical Society A, 128, 110-117.
#' 
#' Robertson, T. and J.D. Cryer. (1974). An iterative procedure for estimating the mode.
#' Journal, American Statistical Association, 69, 1012-1016.
#' 
#' @export

hsmode <- function(x) {
  stopifnot(is.numeric(x) && is.vector(x))
  n <- length(x)
  x <- sort(x)
  
  if (n == 1) {
    return(x)
  } else if (n == 3) {
    dx <- diff(x)
    if (dx[1] == dx[2])
      return(x[2])
    if (dx[1] < dx[2])
      x <- x[1:2]
    else
      x <- x[2:3]
  } else {
    while (n > 2) {
      h <- floor(n/2)
      i <- which.min(x[1:(n-h)+h] - x[1:(n-h)])
      x <- x[i:(i+h)]
      n <- h+1
    }
  }
  
  return(sum(x)/2)
}


#' Shortest half
#' 
#' @param x numeric vector.
#' 
#' @export

shorth <- function(x) {
  stopifnot(is.numeric(x))
  x <- sort(na.omit(x))
  n <- length(x)
  h <- floor(n/2)
  res <- x[(n-h+1):n] - x[1:(n-h)]
  i <- which.min(res)
  c(mean  = mean(x[i:(i+h)]),
    mid   = (x[i] + x[i+h])/2,
    low   = x[i],
    high  = x[i+h],
    range = res[i])
}