1 | mulwmwcrit(mm1, mm2, plotit = TRUE, cop = 3, iter = 1000, alpha = 0.05, SEED = NA, pr = FALSE)
|
mm1 |
|
mm2 |
|
plotit |
|
cop |
|
iter |
|
alpha |
|
SEED |
|
pr |
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 | ##---- 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 (mm1, mm2, plotit = TRUE, cop = 3, iter = 1000, alpha = 0.05,
SEED = NA, pr = FALSE)
{
if (!is.matrix(mm1))
stop("Data are assumed to be stored in a matrix having two or more columns. For univariate data, use the function outbox or out")
if (is.na(SEED))
set.seed(2)
if (!is.na(SEED))
set.seed(SEED)
val <- NA
n1 <- nrow(mm1)
n2 <- nrow(mm2)
for (it in 1:iter) {
ivec1 <- sample(c(1:n1), replace = TRUE)
ivec2 <- sample(c(1:n2), replace = TRUE)
m1 <- mm1[ivec1, ]
m2 <- mm2[ivec2, ]
if (cop == 1) {
if (ncol(m1) > 2) {
center1 <- dmean(m1, tr = 0.5)
center2 <- dmean(m2, tr = 0.5)
}
if (ncol(m1) == 2) {
tempd <- NA
for (i in 1:nrow(m1)) tempd[i] <- depth(m1[i,
1], m1[i, 2], m1)
mdep <- max(tempd)
flag <- (tempd == mdep)
if (sum(flag) == 1)
center1 <- m1[flag, ]
if (sum(flag) > 1)
center1 <- apply(m1[flag, ], 2, mean)
for (i in 1:nrow(m2)) tempd[i] <- depth(m2[i,
1], m2[i, 2], m2)
mdep <- max(tempd)
flag <- (tempd == mdep)
if (sum(flag) == 1)
center2 <- m2[flag, ]
if (sum(flag) > 1)
center2 <- apply(m2[flag, ], 2, mean)
}
}
if (cop == 2) {
center1 <- cov.mcd(m1)$center
center2 <- cov.mcd(m2)$center
}
if (cop == 3) {
center1 <- apply(m1, 2, median)
center2 <- apply(m2, 2, median)
}
center <- (center1 + center2)/2
B <- center1 - center2
if (sum(center1^2) > sum(center2^2))
B <- (0 - 1) * B
BB <- B^2
bot <- sum(BB)
disx <- NA
disy <- NA
if (bot != 0) {
for (j in 1:nrow(m1)) {
AX <- m1[j, ] - center
tempx <- sum(AX * B) * B/bot
disx[j] <- sign(sum(AX * B)) * sqrt(sum(tempx^2))
}
for (j in 1:nrow(m2)) {
AY <- m2[j, ] - center
tempy <- sum(AY * B) * B/bot
disy[j] <- sign(sum(AY * B)) * sqrt(sum(tempy^2))
}
}
m <- outer(disx, disy, FUN = "-")
m <- sign(m)
val[it] <- (1 - mean(m))/2
if (bot == 0)
val[it] <- 0.5
if (pr)
print(paste("Iteration ", it, " of ", iter, " complete"))
}
val <- sort(val)
low <- round(alpha * iter/2) + 1
up <- iter - low
crit <- NA
crit[1] <- val[low]
crit[2] <- val[up]
crit
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.