L1medcen:

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

Usage

L1medcen(X, tol = 1e-08, maxit = 200, m.init = apply(X, 2, median), trace = FALSE)

Arguments

X
tol
maxit
m.init
trace

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (X, tol = 1e-08, maxit = 200, m.init = apply(X, 2, median), 
    trace = FALSE) 
{
    centr <- function(X, m) X - rep(m, each = n)
    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]
        X <- X[inorms, ]
        X. <- X.[inorms, ]
        w <- 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
        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) {
            m <- mold
            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")
    list(center = m)
  }

musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.