library(MASS)

rawdata<-matrix(scan("tab1_1.dat"),ncol=3,byrow=T)
group <- rawdata[,1]
X <- 100 * rawdata[,2:3]
Apf <- X[group==1,]
Af <- X[group==0,]
xbar1 <- apply(Af, 2, mean)
S1 <- var(Af)
N1 <- dim(Af)[1]
xbar2 <- apply(Apf, 2, mean)
S2 <- var(Apf)
N2 <- dim(Apf)[1]
S<-((N1-1)*S1+(N2-1)*S2)/(N1+N2-2)
Sinv=solve(S)

d<-xbar1-xbar2
b <- Sinv %*% d
v <- X %*% b

par(mfcol=c(2,2))
plot(data.frame(v,rep(0,length(v[1]))), col=group+1)

d <- data.frame(rawdata)
names(d)<-c("y","x1","x2")

d$x1 = d$x1 * 100
d$x2 = d$x2 * 100

g<-lda( y ~ x1 + x2, data=d)
v2 <- predict(g, d)

par(mfcol=c(2,2))

y<-seq(100,300,2)
x<-seq(100,300,2)
Xcon <- data.frame(rep(x,length(y)), rep(y, each=length(x)))

names(Xcon)<-c("x1","x2")
u<-as.matrix(Xcon) %*% b
plot(X, type="n")
text(X[group==0,], col=1)
text(X[group==1,], col=2)
title("Manual Calculation")
contour(x,y,matrix(u,length(x),length(y)), add=T, levels=5.87, labex = 0)

Z <- predict(g, Xcon)
zp <- Z$post[,1] - Z$post[,2]

plot(X, type="n")
text(X[group==0,], col=1)
text(X[group==1,], col=2)
title("LDA Calculation")
contour(x,y,matrix(zp,length(x),length(y)), add=T, levels=0, labex = 0)

plot(data.frame(v,rep(0,length(v[1]))), type="n")
text(matrix(c(v[group==0], rep(0, length(v[group==0]))), ncol=2), col=1)
text(matrix(c(v[group==1], rep(0, length(v[group==1]))), ncol=2), col=2)
title("Projection onto 1st Canonical (manual)")
mtext(paste("beta_1=",format(b[1],digits=2), "beta_2=",format(b[2],digits=2)))

plot(data.frame(v2$x,rep(0,length(v[1]))), type="n")
title("Projection onto 1st Canonical (LDA)")
mtext(paste("beta_1=",format(g$sc[1],digits=2), "beta_2=",format(g$sc[2],digits=2)))
text(matrix(c(v2$x[group==0], rep(0, length(v2$x[group==0]))), ncol=2), col=1)
text(matrix(c(v2$x[group==1], rep(0, length(v2$x[group==1]))), ncol=2), col=2)