marpca.sub:

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

Usage

marpca.sub(x, p = ncol(x) - 1, N1 = 3, N2 = 2, tol = 0.001, B = NULL, a = NULL, LSCALE = TRUE)

Arguments

x
p
N1
N2
tol
B
a
LSCALE

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, p = ncol(x) - 1, N1 = 3, N2 = 2, tol = 0.001, B = NULL, 
    a = NULL, LSCALE = TRUE) 
{
    wt.cov <- NULL
    if (!is.null(B)) {
        B <- as.matrix(B)
        if (ncol(B) == 1) 
            B <- t(B)
    }
    n <- nrow(x)
    m <- ncol(x)
    q <- m - p
    if (q < 0) 
        stop("p and q should have values between 1 and ncol(x)")
    hval <- floor((n + m - q + 2)/2)
    DEL <- Inf
    sig0 <- Inf
    if (is.null(B)) {
        if (p > 0 && p < m) {
            B <- matrix(runif(q * m), nrow = q, ncol = m)
            B <- t(ortho(t(B)))
        }
        if (p == 0) 
            B <- diag(rep(1, m))
    }
    it <- 1
    Bx <- matrix(NA, ncol = q, nrow = n)
    for (i in 1:n) Bx[i, ] <- B %*% as.matrix(x[i, ])
    if (is.null(a)) 
        a <- apply(Bx, 2, FUN = "median")
    while (it < N1 + N2 && DEL > tol) {
        r <- NA
        for (i in 1:n) r[i] <- sum(Bx[i, ] - a)^2
        if (LSCALE) 
            sig <- lscale(r, m, q)
        if (!LSCALE) {
            delta <- delta <- (n - m + q - 1)/(2 * n)
            sig <- mscale(r, delta)
        }
        DEL <- 1 - sig/sig0
        sig0 <- sig
        ord.r <- order(r)
        w <- rep(0, n)
        w[ord.r[1:hval]] <- 1
        xx <- x
        for (i in 1:n) xx[i, ] <- x[i, ] * w[i]
        mu <- apply(xx, 2, FUN = "sum")/sum(w)
        Cmat <- matrix(0, nrow = m, ncol = m)
        for (i in 1:n) {
            temp <- w[i] * as.matrix(x[i, ] - mu) %*% t(as.matrix(x[i, 
                ] - mu))
            Cmat <- Cmat + temp
        }
        wt.cov <- Cmat/sum(w)
        if (it > N1) {
            temp <- eigen(wt.cov)
            ord.eig <- order(temp$values)
            for (iq in 1:q) B[iq, ] <- temp$vectors[, ord.eig[iq]]
        }
        a <- B %*% mu
        it <- it + 1
    }
    list(B = B, a = a, var.op = sig, mu = mu, wt.cov = wt.cov)
  }

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