rmba:

Usage Arguments Examples

Usage

1
rmba(x, csteps = 5, na.rm = TRUE, plotit = FALSE)

Arguments

x
csteps
na.rm
plotit

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
##---- 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, csteps = 5, na.rm = TRUE, plotit = FALSE) 
{
    x = as.matrix(x)
    if (na.rm) 
        x = elimna(x)
    p <- dim(x)[2]
    n <- dim(x)[1]
    covs <- var(x)
    mns <- apply(x, 2, mean)
    for (i in 1:csteps) {
        md2 <- mahalanobis(x, mns, covs)
        medd2 <- median(md2)
        mns <- apply(as.matrix(x[md2 <= medd2, ]), 2, mean)
        covs <- var(x[md2 <= medd2, ])
    }
    covb <- covs
    mnb <- mns
    critb <- prod(diag(chol(covb)))
    covv <- diag(p)
    med <- apply(x, 2, median)
    md2 <- mahalanobis(x, center = med, covv)
    medd2 <- median(md2)
    mns <- apply(as.matrix(x[md2 <= medd2, ]), 2, mean)
    covs <- var(x[md2 <= medd2, ])
    for (i in 1:csteps) {
        md2 <- mahalanobis(x, mns, covs)
        medd2 <- median(md2)
        mns <- apply(as.matrix(x[md2 <= medd2, ]), 2, mean)
        covs <- var(x[md2 <= medd2, ])
    }
    crit <- prod(diag(chol(covs)))
    if (crit < critb) {
        critb <- crit
        covb <- covs
        mnb <- mns
    }
    rd2 <- mahalanobis(x, mnb, covb)
    const <- median(rd2)/(qchisq(0.5, p))
    covb <- const * covb
    rd2 <- mahalanobis(x, mnb, covb)
    up <- qchisq(0.975, p)
    rmnb <- apply(as.matrix(x[rd2 <= up, ]), 2, mean)
    rcovb <- var(x[rd2 <= up, ])
    rd2 <- mahalanobis(x, rmnb, rcovb)
    const <- median(rd2)/(qchisq(0.5, p))
    rcovb <- const * rcovb
    rd2 <- mahalanobis(x, rmnb, rcovb)
    up <- qchisq(0.975, p)
    rmnb <- apply(as.matrix(x[rd2 <= up, ]), 2, mean)
    rcovb <- var(x[rd2 <= up, ])
    rd2 <- mahalanobis(x, rmnb, rcovb)
    const <- median(rd2)/(qchisq(0.5, p))
    rcovb <- const * rcovb
    cor.b = NULL
    temp = outer(sqrt(diag(rcovb)), sqrt(diag(rcovb)), "*")
    if (min(diag(rcovb) > 0)) 
        cor.b = rcovb/temp
    list(center = rmnb, cov = rcovb, cor = cor.b)
  }

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