library(MASS)
#
nn<-c(100,50,50)
# number of observations per class
corrx<-(-.8)
# try 0, -.8, 0.8, -0.5 etc
# play with this number of change the correlation between the x-variables
Sigma<-cbind(c(1,corrx),c(corrx,1))
# 
mu1<-c(0,0)
mu2<-c(1.5,0.75)
mu3<-c(3,3)
# mean location for each class, try different values
x1<-mvrnorm(nn[1],mu1,Sigma)
x2<-mvrnorm(nn[2],mu2,Sigma)
x3<-mvrnorm(nn[3],mu3,Sigma)
# creating the 2D data for 3 classes.
# change the number to change the location of the means
y<-c(rep(0,nn[1]),rep(1,nn[2]),rep(2,nn[3]))
xx<-(rbind(x1,x2,x3))
#
lly<-lda(xx,y)
# run lda on data
print(c("Directions:"))
print(lly$sca[,1])
print(lly$sca[,2])
# the discriminant variables - new coordinate system
par(mfrow=c(1,1))
plot(xx[,1],xx[,2],col="blue",xlab="X1",ylab="X2")
points(x2[,1],x2[,2],col="red")
points(x3[,1],x3[,2],col="green")
points(lly$means[1,1],lly$means[1,2],cex=2)
points(lly$means[2,1],lly$means[2,2],cex=2)
points(lly$means[3,1],lly$means[3,2],cex=2)
# plot the data in the original coordinate system
p<-locator()
abline(lly$scal[1,1],lly$scal[2,1])
abline(abs(lly$scal[2,2]),-1*abs(lly$scal[1,2]))
# these are the vectors to project on for lda classifiction
p<-locator()
discvar<-xx%*%lly$scaling
plot(discvar[y==0,1],discvar[y==0,2],col="blue",xlim=c(min(discvar[,1]),max(discvar[,1])),ylim=c(min(discvar[,2]),max(discvar[,2])),xlab="disc1",ylab="disc2")
points(discvar[y==1,1],discvar[y==1,2],col="red")
points(discvar[y==2,1],discvar[y==2,2],col="green")
centroids<-lly$means%*%lly$scaling
points(centroids[1,1],centroids[1,2],cex=2)
points(centroids[2,1],centroids[2,2],cex=2)
points(centroids[3,1],centroids[3,2],cex=2)
xpl1<-seq(-15,15)
xpl2<-(centroids[2,1]^2+centroids[2,2]^2-centroids[1,1]^2-centroids[1,2]^2+2*(centroids[1,1]-centroids[2,1])*xpl1+2*log(length(y[y==0])/length(y[y==1])))/(2*(centroids[2,2]-centroids[1,2]))
lines(xpl1,xpl2,col="purple")
xpl2<-(centroids[3,1]^2+centroids[3,2]^2-centroids[1,1]^2-centroids[1,2]^2+2*(centroids[1,1]-centroids[3,1])*xpl1+2*log(length(y[y==0])/length(y[y==2])))/(2*(centroids[3,2]-centroids[1,2]))
lines(xpl1,xpl2,col="darkblue")
xpl2<-(centroids[3,1]^2+centroids[3,2]^2-centroids[2,1]^2-centroids[2,2]^2+2*(centroids[2,1]-centroids[3,1])*xpl1+2*log(length(y[y==1])/length(y[y==2])))/(2*(centroids[3,2]-centroids[2,2]))
lines(xpl1,xpl2,col="orange")
# plotting the data in the new coordinate system
#
# testset
p<-locator()
nn<-1000
tx1<-mvrnorm(nn,mu1,Sigma)
tx2<-mvrnorm(nn,mu2,Sigma)
tx3<-mvrnorm(nn,mu3,Sigma)
ty<-c(rep(0,nn),rep(1,nn),rep(2,nn))
txx<-(rbind(tx1,tx2,tx3))
# lda
ply<-predict(lly,newdata=txx)
# predict on new data
plot(xx[,1],xx[,2],col="blue")
points(x2[,1],x2[,2],col="red")
points(x3[,1],x3[,2],col="green")
points(lly$means[1,1],lly$means[1,2],cex=2)
points(lly$means[2,1],lly$means[2,2],cex=2)
points(lly$means[3,1],lly$means[3,2],cex=2)
points(txx[,1],txx[,2],col="lightblue",pch=2)
points(tx2[,1],tx2[,2],col="orange",pch=2)
points(tx3[,1],tx3[,2],col="darkgreen",pch=2)
# plotting the data in the original coordinate system
print(c("LDA misclass",round(length(ply$class[ply$class!=ty])/length(ty),3)))
p<-locator()
# error rate
#
discvart<-txx%*%lly$scaling
plot(discvart[ty==0,1],discvart[ty==0,2],col="lightblue",xlim=c(min(discvart[,1]),max(discvart[,1])),ylim=c(min(discvart[,2]),max(discvart[,2])))
points(discvart[ty==1,1],discvart[ty==1,2],col="orange")
points(discvart[ty==2,1],discvart[ty==2,2],col="darkgreen")
centroids<-lly$means%*%lly$scaling
points(centroids[1,1],centroids[1,2],cex=2)
points(centroids[2,1],centroids[2,2],cex=2)
points(centroids[3,1],centroids[3,2],cex=2)
xpl1<-seq(-15,15)
xpl2<-(centroids[2,1]^2+centroids[2,2]^2-centroids[1,1]^2-centroids[1,2]^2+2*(centroids[1,1]-centroids[2,1])*xpl1+2*log(length(y[y==0])/length(y[y==1])))/(2*(centroids[2,2]-centroids[1,2]))
lines(xpl1,xpl2,col="purple")
xpl2<-(centroids[3,1]^2+centroids[3,2]^2-centroids[1,1]^2-centroids[1,2]^2+2*(centroids[1,1]-centroids[3,1])*xpl1+2*log(length(y[y==0])/length(y[y==2])))/(2*(centroids[3,2]-centroids[1,2]))
lines(xpl1,xpl2,col="darkblue")
xpl2<-(centroids[3,1]^2+centroids[3,2]^2-centroids[2,1]^2-centroids[2,2]^2+2*(centroids[2,1]-centroids[3,1])*xpl1+2*log(length(y[y==1])/length(y[y==2])))/(2*(centroids[3,2]-centroids[2,2]))
lines(xpl1,xpl2,col="orange")
# plotting in new coordinate system
#
# project only onto first
p<-locator()
ply1<-predict(lly,newdata=txx,dimen=1)
print(c("LDA misclass 1 direction",round(length(ply1$class[ply1$class!=ty])/length(ty),3)))
plot(discvart[ty==0,1],discvart[ty==0,2],col="lightblue",xlim=c(min(discvart[,1]),max(discvart[,1])),ylim=c(min(discvart[,2]),max(discvart[,2])))
points(discvart[ty==1,1],discvart[ty==1,2],col="orange")
points(discvart[ty==2,1],discvart[ty==2,2],col="darkgreen")
centroids<-lly$means%*%lly$scaling
points(centroids[1,1],centroids[1,2],cex=2)
points(centroids[2,1],centroids[2,2],cex=2)
points(centroids[3,1],centroids[3,2],cex=2)
cc1<-(centroids[2,1]+centroids[1,1])/2
cc2<-(centroids[3,1]+centroids[1,1])/2
cc3<-(centroids[3,1]+centroids[2,1])/2
abline(v=cc1,col="purple")
abline(v=cc2,col="darkblue")
abline(v=cc3,col="orange")
# onto second 
p<-locator()
plot(discvart[ty==0,1],discvart[ty==0,2],col="lightblue",xlim=c(min(discvart[,1]),max(discvart[,1])),ylim=c(min(discvart[,2]),max(discvart[,2])))
points(discvart[ty==1,1],discvart[ty==1,2],col="orange")
points(discvart[ty==2,1],discvart[ty==2,2],col="darkgreen")
centroids<-lly$means%*%lly$scaling
points(centroids[1,1],centroids[1,2],cex=2)
points(centroids[2,1],centroids[2,2],cex=2)
points(centroids[3,1],centroids[3,2],cex=2)
cc1<-(centroids[2,2]+centroids[1,2])/2
cc2<-(centroids[3,2]+centroids[1,2])/2
cc3<-(centroids[3,2]+centroids[2,2])/2
abline(h=cc1,col="purple")
abline(h=cc2,col="darkblue")
abline(h=cc3,col="orange")