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.huber <- vglm(y ~ x2, huber(meth = 3), hdata, trace = TRUE)
coef(fit.huber, matrix = TRUE)
summary(fit.huber)
# Plot the results
plot(y ~ x2, hdata, col = "blue", las = 1)
lines(fitted(fit.huber) ~ 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 = "red", lwd = 2, lty = "dashed")
legend("bottomright", legend = c("truth", "huber", "lm"),
col = c("red", "darkgreen", "lavender"),
lty = c("dashed", "solid", "solid"), lwd = c(2, 2, 3))
Run the code above in your browser using DataLab