R/L1median.R

In musto101/wilcox_R: What the package does (short line)

L1median <-
function(X, tol = 1e-08, maxit = 200, m.init = apply(X, 2, median),
                     trace = FALSE)
{
  ## L1MEDIAN calculates the multivariate L1 median
  ## I/O: mX=L1median(X,tol);
  ##
  ## X  : the data matrix
  ## tol: the convergence criterium:
  ##      the iterative process stops when ||m_k - m_{k+1}|| < tol.
  ## maxit: maximum number of iterations
  ## init.m: starting value for m; typically coordinatewise median
  ##
  ## Ref: Hossjer and Croux (1995)
  ##  "Generalizing Univariate Signed Rank Statistics for Testing
  ##   and Estimating a Multivariate Location Parameter";
  ##   Non-parametric Statistics, 4, 293-308.
  ##
  ## Implemented by Kristel Joossens
  ## Many thanks to Martin Maechler for improving the program!

  ## slightly faster version of 'sweep(x, 2, m)':
  centr <- function(X,m) X - rep(m, each = n)
  ## computes objective function in m based on X and a:
  mrobj <- function(X,m) sum(sqrt(rowSums(centr(X,m)^2)))

  d <- dim(X); n <- d[1]; p <- d[2]
  m <- m.init
  if(!is.numeric(m) || length(m) != p)
      stop("'m.init' must be numeric of length p =", p)
  k <- 1
  if(trace) nstps <- 0
  while (k <= maxit) {
    mold <- m
    obj.old <- if(k == 1) mrobj(X,mold) else obj
    X. <- centr(X, m)
    Xnorms <- sqrt(rowSums(X. ^ 2))
    inorms <- order(Xnorms)
    dx <- Xnorms[inorms] # smallest first, i.e., 0's if there are
    X  <- X [inorms,]
    X. <- X.[inorms,]
    ## using 1/x weighting {MM: should this be generalized?}
    w <- ## (0 norm -> 0 weight) :
        if (all(dn0 <- dx != 0))  1/dx
        else c(rep.int(0, length(dx)- sum(dn0)), 1/dx[dn0])
    delta <- colSums(X. * rep(w,p)) / sum(w)
    nd <- sqrt(sum(delta^2))

    maxhalf <- if (nd < tol) 0 else ceiling(log2(nd/tol))
    m <- mold + delta    # computation of a new estimate
    ## If step 'delta' is too far, we try halving the stepsize
    nstep <- 0
    while ((obj <- mrobj(X, m)) >= obj.old && nstep <= maxhalf) {
      nstep <- nstep+1
      m <- mold + delta/(2^nstep)
    }
    if(trace) {
        if(trace >= 2)
            cat(sprintf("k=%3d obj=%19.12g m=(",k,obj),
                paste(formatC(m),collapse=","),
                ")", if(nstep) sprintf(" nstep=%2d halvings",nstep) else "",
                "\n", sep="")
        nstps[k] <- nstep
    }
    if (nstep > maxhalf) { ## step halving failed; keep old
        m <- mold
        ## warning("step halving failed in ", maxhalf, " steps")
        break
      }
    k <- k+1
  }
  if (k > maxit) warning("iterations did not converge in ", maxit, " steps")
  if(trace == 1)
      cat("needed", k, "iterations with a total of",
          sum(nstps), "stepsize halvings\n")
  return(m)
}