dataset <- read.table("http://azzalini.stat.unipd.it/Book-DM/classes.dat", header = TRUE)
x1 <- dataset[ , 1]
x2 <- dataset[ , 2]
gr <- (2 - dataset[ , 3])
K  <- 3
n <- nrow(dataset)
x1 <- z1 <- dataset[ , 1]
x2 <- z2 <- dataset[ , 2]
gr <- dataset[,3]
group <- factor(gr)
n.grid <- 800

library(MASS)
d1 <- lda(group ~ x1 + x2)
p <- pred.square(d1, x1, x2, n.grid)
pred <- array(p$pred$posterior, c(n.grid, n.grid, 3))
ind <- apply(pred, c(1, 2), order)[3, , ] 

# figure 5.11a
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]),
     ylab = expression(italic(x)[2]))
for(k in 1:3) points(dataset[gr == k, 1:2], pch = 15 + k, col = 1 + k)   
contour(p$x, p$y, ind, add = TRUE, drawlabels = FALSE, nlevels = 2, lty = 2)

d2 <- lda(group ~ x1 + x2 + I(x1^2) + I(x1*x2) + I(x2^2))
p <- pred.square(d2, x1, x2, n.grid)
pred <- array(p$pred$posterior, c(n.grid, n.grid, 3))
ind <- apply(pred, c(1, 2), order)[3, , ] 

# figure 5.11b
plot(dataset[, 1:2], type = "n", xlab = expression(italic(x)[1]),
     ylab = expression(italic(x)[2]))
for(k in 1:3) points(dataset[gr == k, 1:2], pch = 15 + k, col = 1 + k)   
contour(p$x, p$y, ind, add = TRUE, drawlabels = FALSE, nlevels = 2, lty = 2)


### QDA ###

d3 <- qda(group ~ x1 + x2)
p <- pred.square(d3, x1, x2, n.grid)
pred <- array(p$pred$posterior, c(n.grid, n.grid, 3))
ind <- apply(pred, c(1, 2), order)[3, , ] 

# figure 5.12a
plot(dataset[, 1:2], type = "n", xlab = expression(italic(x)[1]),
     ylab = expression(italic(x)[2]))
for(k in 1:3) points(dataset[gr == k, 1:2], pch = 15 + k, col = 1  + k)   
contour(p$x, p$y, ind, add = TRUE, drawlabels = FALSE, nlevels = 2)

d4 <- qda(group ~ x1 + x2 + I(x1^2) + I(x1*x2) + I(x2^2))
p <- pred.square(d4, x1, x2, n.grid)
pred <- array(p$pred$posterior, c(n.grid, n.grid, 3))
ind <- apply(pred, c(1, 2), order)[3, , ] 

#  figure 5.12b
plot(dataset[ , 1:2], type = "n", xlab = expression(italic(x)[1]),
     ylab = expression(italic(x)[2]))
for(k in 1:3) points(dataset[gr == k, 1:2], pch = 15 + k, col = 1 + k)   
contour(p$x, p$y, ind, add = TRUE, drawlabels = FALSE, nlevels = 2, lty = 2)




f1 <- as.formula(paste("choice ~ ", paste(names(juice1)[-1], collapse = "+"),
                       collapse = NULL))
d1 <- lda(f1, data = juice1)
d2 <- qda(f1, data = juice1)
p1 <- predict(d1, newdata = juice2)
p2 <- predict(d2, newdata = juice2)

# table 5.7
confusion.matrix(p1$class, juice2[ , "choice"])
confusion.matrix(p2$class, juice2[ , "choice"])

g <- as.numeric(juice2[ , 1] == "MM")
lr1 <- lift.roc(p1$post[ , 2], g, type = "bin", plot.it = FALSE)
lr2 <- lift.roc(p2$post[ , 2], g, type = "bin", plot.it = FALSE)

# figure 5.13b
plot(lr1[[3]], lr1[[4]], type = "b", xlim = c(0,1), pch = 16,
     ylim = c(0, 1), cex = 1.5, xlab = "1-specificity",
     ylab = "Sensibility", col = 2)
lines(lr2[[3]], lr2[[4]], type = "b", cex = 1.5, col = 3, pch = 17)
abline(0, 1, lty = 2, col = 4)
legend('bottomright', legend = c("LDA", "QDA"), col = 2:3, pch = 16:17)