library(MASS)

##Generate training data
set.seed(401)

nTrain=1000
y=runif(nTrain,0,1)>0.3
x=matrix(NA,nTrain,2)
x[y==1,]=t(replicate(sum(y==1),mvrnorm(1,mu = c(0,1),Sigma = matrix(c(1.9,0.5,0.5,0.3),2,2))))
x[y==0,]=t(replicate(sum(y==0),mvrnorm(1,mu = c(-2,1.5),Sigma = matrix(c(1,-0.25,-0.25,0.1),2,2))))

##Generate test data
nTest=500
yTest=runif(nTest,0,1)>0.3
xTest=matrix(NA,nTest,2)
xTest[yTest==1,]=t(replicate(sum(yTest==1),mvrnorm(1,mu = c(0,1),Sigma = matrix(c(1.9,0.5,0.5,0.3),2,2))))
xTest[yTest==0,]=t(replicate(sum(yTest==0),mvrnorm(1,mu = c(-2,1.5),Sigma = matrix(c(1,-0.25,-0.25,0.1),2,2))))

##fit LDA model
ajusteLinear=lda(x=x,grouping = y)

##Plot results with the separting hyperplane line
plot(x[y==1,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)
points(x[y==0,],pch=18,col=4)
coeff=coefficients(ajusteLinear)
abline(a=coeff[1]+2,b=coeff[1],lwd=4)

predictedY=predict(ajusteLinear,xTest)

##plot training data
plot(x[y==1,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7,main="Conjunto de Treinamento")
points(x[y==0,],pch=18,col=4)

##plot test data
plot(xTest,pch=18,col=1,xlim=c(-5,3),ylim=c(-1,4),xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7,main="Conjunto de Teste")

##plot results: test different cutoff points
par(mfrow=c(1,2))
plot(xTest[predictedY$posterior[,1]<0.02,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),main="Predito (Corte: 0.02)",xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)
points(xTest[predictedY$posterior[,1]>0.02,],pch=18,col=4)

plot(xTest[yTest==1,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),main="Real",xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)
points(xTest[yTest==0,],pch=18,col=4)

##plot ROC curve
par(mfrow=c(1,1))
library(ROCR)
predLinear <- prediction(predictedY$posterior[,2], yTest)
perfLinear <- performance(predLinear, measure = "tpr", x.measure = "fpr") 
plot(perfLinear, col=2,lwd=4,xlab="1-Especificidade",ylab="Sensibilidade",cex.lab=1.7)
abline(a=0,b=1,lwd=3)

##fit QDA model and predict in the tet data
ajusteQuadratico=qda(x=x,grouping = y)
predictedY=predict(ajusteQuadratico,xTest)

##plot results: test different cutoff points
par(mfrow=c(1,2))
plot(xTest[predictedY$posterior[,1]<0.03,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),main="Predito (Corte: 0.03)",xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)
points(xTest[predictedY$posterior[,1]>0.03,],pch=18,col=4,xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)

plot(xTest[yTest==1,],pch=18,col=2,xlim=c(-5,3),ylim=c(-1,4),main="Real",xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)
points(xTest[yTest==0,],pch=18,col=4,xlab=expression(x[1]),ylab=expression(x[2]),cex.lab=1.7)

##plot ROC curve
par(mfrow=c(1,1))
predQuad <- prediction(predictedY$posterior[,2], yTest)
perfQuad <- performance(predQuad, measure = "tpr", x.measure = "fpr") 
plot(perfLinear, col=2,lwd=4,xlab="1-Especificidade",ylab="Sensibilidade",cex.lab=1.7)
abline(a=0,b=1,lwd=3)
plot(perfQuad, add=TRUE,col=4,lwd=3)
legend("bottomright",c("Linear","Quadrática"),col=c(2,4),lwd=3,cex=1.7,bty = "n")