1 |
x |
|
tr |
|
grp |
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 | ##---- 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, tr = 0.2, grp = c(1:length(x)))
{
if (!is.list(x) && !is.matrix(x))
stop("Data must be stored in a matrix or in list mode.")
if (is.list(x)) {
J <- length(grp)
m1 <- matrix(x[[grp[1]]], length(x[[grp[1]]]), 1)
for (i in 2:J) {
m2 <- matrix(x[[grp[i]]], length(x[[i]]), 1)
m1 <- cbind(m1, m2)
}
}
if (is.matrix(x)) {
if (length(grp) < ncol(x))
m1 <- as.matrix(x[, grp])
if (length(grp) >= ncol(x))
m1 <- as.matrix(x)
J <- ncol(x)
}
m2 <- matrix(0, nrow(m1), ncol(m1))
xvec <- 1
g <- floor(tr * nrow(m1))
for (j in 1:ncol(m1)) {
m2[, j] <- winval(m1[, j], tr)
xvec[j] <- mean(m1[, j], tr)
}
xbar <- mean(xvec)
qc <- (nrow(m1) - 2 * g) * sum((xvec - xbar)^2)
m3 <- matrix(0, nrow(m1), ncol(m1))
m3 <- sweep(m2, 1, apply(m2, 1, mean))
m3 <- sweep(m3, 2, apply(m2, 2, mean))
m3 <- m3 + mean(m2)
qe <- sum(m3^2)
test <- (qc/(qe/(nrow(m1) - 2 * g - 1)))
v <- winall(m1)$cov
vbar <- mean(v)
vbard <- mean(diag(v))
vbarj <- 1
for (j in 1:J) {
vbarj[j] <- mean(v[j, ])
}
A <- J * J * (vbard - vbar)^2/(J - 1)
B <- sum(v * v) - 2 * J * sum(vbarj^2) + J * J * vbar^2
ehat <- A/B
etil <- (nrow(m2) * (J - 1) * ehat - 2)/((J - 1) * (nrow(m2) -
1 - (J - 1) * ehat))
etil <- min(1, etil)
df1 <- (J - 1) * etil
df2 <- (J - 1) * etil * (nrow(m2) - 2 * g - 1)
siglevel <- 1 - pf(test, df1, df2)
list(test = test, df = c(df1, df2), p.value = siglevel, tmeans = xvec,
ehat = ehat, etil = etil)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.