1 |
m |
|
pts |
|
plotit |
|
cop |
|
center |
|
xlab |
|
ylab |
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 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 | ##---- 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 (m, pts = NA, plotit = TRUE, cop = 3, center = NA, xlab = "VAR 1",
ylab = "VAR 2")
{
library(MASS)
if (cop != 2 && cop != 3 && cop != 4)
stop("Only cop=2, 3 or 4 is allowed")
if (is.list(m))
stop("Store data in a matrix; might use function listm")
m <- as.matrix(m)
pts <- as.matrix(pts)
if (!is.na(pts[1]))
remm <- m
nm <- nrow(m)
nm1 <- nm + 1
if (!is.na(pts[1])) {
if (ncol(m) != ncol(pts))
stop("Number of columns of m is not equal to number of columns for pts")
}
m <- elimna(m)
m <- as.matrix(m)
if (ncol(m) == 1)
dep <- unidepth(as.vector(m[, 1]), pts = pts)
if (ncol(m) > 1) {
if (is.na(center[1])) {
if (cop == 2) {
center <- cov.mcd(m)$center
}
if (cop == 4) {
center <- cov.mve(m)$center
}
if (cop == 3) {
center <- apply(m, 2, median)
}
}
if (is.na(pts[1])) {
mdep <- matrix(NA, nrow = nrow(m), ncol = nrow(m))
}
if (!is.na(pts[1])) {
mdep <- matrix(NA, nrow = nrow(m), ncol = nrow(pts))
}
for (i in 1:nrow(m)) {
B <- m[i, ] - center
dis <- NA
BB <- B^2
bot <- sum(BB)
if (bot != 0) {
if (is.na(pts[1])) {
for (j in 1:nrow(m)) {
A <- m[j, ] - center
temp <- sum(A * B) * B/bot
dis[j] <- sign(sum(A * B)) * sqrt(sum(temp^2))
}
}
if (!is.na(pts[1])) {
m <- rbind(remm, pts)
for (j in 1:nrow(m)) {
A <- m[j, ] - center
temp <- sum(A * B) * B/bot
dis[j] <- sign(sum(A * B)) * sqrt(sum(temp^2))
}
}
if (is.na(pts[1]))
mdep[i, ] <- unidepth(dis)
if (!is.na(pts[1])) {
mdep[i, ] <- unidepth(dis[1:nm], dis[nm1:nrow(m)])
}
}
if (bot == 0)
mdep[i, ] <- rep(0, ncol(mdep))
}
dep <- apply(mdep, 2, min)
if (ncol(m) == 2 && is.na(pts[1])) {
flag <- chull(m)
dep[flag] <- min(dep)
}
}
if (ncol(m) == 2) {
if (is.na(pts[1]) && plotit) {
plot(m, xlab = xlab, ylab = ylab)
points(center[1], center[2], pch = "+")
x <- m
temp <- dep
flag <- (temp >= median(temp))
xx <- x[flag, ]
xord <- order(xx[, 1])
xx <- xx[xord, ]
temp <- chull(xx)
xord <- order(xx[, 1])
xx <- xx[xord, ]
temp <- chull(xx)
lines(xx[temp, ])
lines(xx[c(temp[1], temp[length(temp)]), ])
}
}
dep <- round(dep * nrow(m))/nrow(m)
dep
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.