psout <- TRUE # generate postscript file "out.eps"

library(MASS)
library(datasets)
data(iris3)
iris3 <- iris3[,1:2,];
Iris <- data.frame( rbind(iris3[,,1], iris3[,,2], iris3[,,3]),
                    class = c( 	rep("s",nrow(iris3[,,1])),
						 		rep("c",nrow(iris3[,,2])),
								rep("v",nrow(iris3[,,3]))))

# z <- lda(class ~ ., Iris, prior = c(1,1,1)/3 )
# L <- predict(z,Iris)

r1 <- range(Iris[,1])
d <- diff(r1)*0.1
r1 <- c(r1[1]-d,r1[2]+d)
r2 <- range(Iris[,2])
d <- diff(r2)*0.1
r2 <- c(r2[1]-d,r2[2]+d)

x <- seq(r1[1],r1[2],length=100)
y <- seq(r2[1],r2[2],length=100)

# cartesian product !!
xy <- data.frame(as.matrix(expand.grid(x,y)))
colnames(xy) <- colnames(Iris[,1:2]);

bg <- 0.6 # background contrast

#############################################################################
#
#                       Linear Discriminant Analysis
#
#############################################################################
if(psout)
{
	#postscript(file="out.eps",paper="special",width=6,height=6,
	#	onefile=FALSE,horizontal=FALSE)
	pdf(file="qda_iris.pdf")
}

matplot(r1,r2, 
			type= "n", xlab = colnames(Iris)[1], ylab = colnames(Iris)[2],
             main = "QDA: Iris Dataset")

#z <- lda(class ~ ., Iris, prior = c(1,1,1)/3 )
z <- qda(class ~ ., Iris, prior = c(1,1,1)/3 )
			 
zz <- predict(z,xy)
zc <- matrix(NA,ncol=1,nrow=nrow(xy))

zc[which(zz$class=="s")] <- rgb(1.0,bg,bg)
zc[which(zz$class=="c")] <- rgb(bg,bg,1.0)
zc[which(zz$class=="v")] <- rgb(bg,1.0,bg)

dx <- x[2]-x[1];
dy <- y[2]-y[1];
rect(xy[,1],xy[,2],xy[,1]+dx,xy[,2]+dy,col=zc,density=-1,border=FALSE,lwd=0)

Is <- which(Iris[,"class"] == "s")
Ic <- which(Iris[,"class"] == "c")
Iv <- which(Iris[,"class"] == "v")
points(Iris[Is,1], Iris[Is,2], pch = "+", col = rgb(0.5,0.0,0.0))
points(Iris[Ic,1], Iris[Ic,2], pch = "o", col = rgb(0.0,0.0,0.5))
points(Iris[Iv,1], Iris[Iv,2], pch = "x", col = rgb(0.0,0.5,0.0))

if( psout )
{
	dev.off()
}