require(MASS)
require(ggplot2)
data(iris)

newData <- iris[51:150,]
newData[,5] <- as.factor(as.vector(newData[,5]))

irisGlm <- glm(Species~.,data=newData,family=binomial)
irisLda <- lda(Species~.,data=newData)

newObs <- data.frame(Sepal.Length=5.9,Sepal.Width=2.7,Petal.Length=4.6,Petal.Width=1.3)
predict(irisGlm,newObs,type="response")
###The automatic encoding will result in virginica being 1. The probability above
###is the estimated probability of virginica. It is close to 0, hence we classify as
###versicolor
predict(irisLda,newObs)

summary(irisGlm)
summary(glm(Species~.,data=newData[,-1],family=binomial))
###Etc.

###Q3.2

trainPosterior <- data.frame(logis.virginica=predict(irisGlm,type="response"),lda=predict(irisLda)$posterior[,2])
format(trainPosterior,scientific=FALSE)
ggplot(data=melt(cbind(data.frame(index=1:100,diff=trainPosterior[,1] - trainPosterior[,2]+0.5),trainPosterior),id="index"),aes(x=index,y=value,col=variable)) + geom_point()


###Q3.3

mu1Hat <- irisLda$means[2,]
mu0Hat <- irisLda$means[1,]

resid <- as.matrix(newData[,-5] - cbind(matrix(c(rep(1,50),rep(0,50))),matrix(c(rep(0,50),rep(1,50)))) %*% irisLda$means)

SigmaHat <- t(resid) %*% resid/98

tauHatBeta0 <- (t(mu0Hat) %*% solve(SigmaHat,mu0Hat) - t(mu1Hat) %*% solve(SigmaHat,mu1Hat)   )/2
tauHatBeta <- solve(SigmaHat,mu1Hat - mu0Hat)
hat <- data.frame(tauHat = c(tauHatBeta0,tauHatBeta),Hat = irisGlm$coef)
row.names(hat) <- c("beta_0","beta_1","beta_2","beta_3","beta_4")
print(hat)