X <- as.matrix(subset(iris, select = -Species))
W <- crossprod(resid(aov(X ~ Species, iris)))
B <- crossprod(resid(aov(X ~ 1, iris))) - W
n <- nrow(iris)
p <- length(levels(iris$Species))
(ev <- eigen2(B/(p - 1), W/(n - p))) ## hand-made discriminant analysis
DF <- X %*% ev$vectors[, 1:2]
with(iris, {
plot(DF, col = Species, pch = 20,
xlab = expression(DF[1]), ylab = expression(DF[2]))
legend("topleft", levels(Species), pch = 20, col = 1:3)
})
Run the code above in your browser using DataLab