Data <- read.delim("../../data/iris.txt")

cols=c(rgb(1,0,0,alpha=.5),rgb(0,1,0,alpha=.5),rgb(0,0,1,alpha=.5))
pairs(Data[,1:4],col=cols[as.numeric(Data$Species)],pch=19,gap=0)

require(MASS)
fit <- lda(Species~.,Data)
plot(fit,abbrev=TRUE)
fit

table(predict(fit)$class,Data$Species)

ind <- c(sample(which(Data$Species=="setosa"),5),
         sample(which(Data$Species=="versicolor"),5),
         sample(which(Data$Species=="virginica"),5))
Train <- Data[ind,]
Test <- Data[-ind,]
fit <- lda(Species~.,Train)
pred <- predict(fit,Test)
table(pred$class,Test$Species)
xtable(table(pred$class,Test$Species))

## comparison w/ multinom
require(nnet)
fit.m <- multinom(Species~.,Data)
table(predict(fit.m),Data$Species)

mcr.tab <- function(X){1-sum(diag(X))/sum(X)}
N <- 1000
MC <- matrix(NA,nrow=N,ncol=2)
for (i in 1:N)
  {
    ind <- c(sample(which(Data$Species=="setosa"),5),
             sample(which(Data$Species=="versicolor"),5),
             sample(which(Data$Species=="virginica"),5))
    Train <- Data[ind,]
    Test <- Data[-ind,]
    fit <- lda(Species~.,Train)
    y1 <- predict(fit,Test)$class
    fit.m <- multinom(Species~.,Train,trace=FALSE)
    y2 <- predict(fit.m,Test)
    MC[i,1] <- sum(y1!=Test$Species)
    MC[i,2] <- sum(y2!=Test$Species)
  }
apply(MC,2,sum)/(N*135)