1 |
x |
|
csteps |
|
na.rm |
|
plotit |
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)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.