set.seed(1231); NN <- 30; coef1 <- 1; coef2 <- 10
hdata <- data.frame(x2 = sort(runif(NN)))
hdata <- transform(hdata, y = rhuber(NN, mu = coef1 + coef2 * x2))
hdata$x2[1] <- 0.0 # Add an outlier
hdata$y[1] <- 10
fit.huber2 <- vglm(y ~ x2, huber2(imethod = 3), hdata, trace = TRUE)
fit.huber1 <- vglm(y ~ x2, huber1(imethod = 3), hdata, trace = TRUE)
coef(fit.huber2, matrix = TRUE)
summary(fit.huber2)
# Plot the results
plot(y ~ x2, hdata, col = "blue", las = 1)
lines(fitted(fit.huber2) ~ x2, hdata, col = "darkgreen", lwd = 2)
fit.lm <- lm(y ~ x2, hdata) # Compare to a LM:
lines(fitted(fit.lm) ~ x2, hdata, col = "lavender", lwd = 3)
# Compare to truth:
lines(coef1 + coef2 * x2 ~ x2, hdata, col = "orange", lwd = 2, lty = "dashed")
legend("bottomright", legend = c("truth", "huber", "lm"),
col = c("orange", "darkgreen", "lavender"),
lty = c("dashed", "solid", "solid"), lwd = c(2, 2, 3))
Run the code above in your browser using DataLab