1 |
x |
|
y |
|
bend |
|
SEED |
|
xout |
|
outfun |
|
... |
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 | ##---- 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, y, bend = 1.345, SEED = FALSE, xout = FALSE, outfun = out,
...)
{
library(MASS)
if (SEED)
set.seed(12)
x <- as.matrix(x)
p <- ncol(x)
m <- elimna(cbind(x, y))
x <- m[, 1:p]
p1 <- p + 1
y <- m[, p1]
if (xout) {
x <- as.matrix(x)
flag <- outfun(x, plotit = plotit, ...)$keep
x <- x[flag, ]
y <- y[flag]
x <- as.matrix(x)
}
x <- as.matrix(x)
cutoff <- bend
mve <- vector("list")
if (ncol(x) == 1) {
mve$center <- median(x)
mve$cov <- mad(x)^2
}
if (ncol(x) >= 2)
mve <- cov.mve(x)
reg0 <- ltsgreg(x, y)
rob <- 1
mx <- mve$center
rob <- mahalanobis(x, mx, mve$cov)
k21 <- qchisq(0.95, p)
c62 <- k21/rob
vecone <- c(rep(1, length(y)))
c30 <- pmin(vecone, c62)
k81 <- median(abs(reg0$residuals))
k72 <- 1.4826 * (1 + (5/(length(y) - p - 1))) * k81
c60 <- reg0$residuals/(k72 * c30)
cvec <- c(rep(cutoff, length(y)))
c27 <- pmin(cvec, c60)
c27 <- pmax(-1 * cutoff, c27)
c66 <- ifelse(abs(c27) <= bend, 1, 0)
m1 <- cbind(1, x)
m2 <- t(m1)
m5 <- diag(c30)
m4 <- diag(c66)
m6 <- m4 %*% m1
m7 <- m2 %*% m6
m8 <- solve(m7)
m9 <- m8 %*% m2
m9 <- m9 %*% m5
m10 <- m9 %*% c27
c20 <- m10 * k72
c21 <- reg0$coef + c20
res <- y - m1 %*% c21
list(coef = t(c21), residuals = res)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.