#  This file demonstrates quadratic discriminant function
#  theory.  
#      No data is involved, just the densities.
#      Two and three densities illustrated 
#          with two different covariance matrices
#          with three covariance as follows:
#              two equal, one different
#              all three different covariance matrices
#
#  Authors:
#       H.A. Chipman, 2003
#       R.W. Oldford, 2004
#
#

# First determine the parameters of two multivariate normals (different Sigmas)
mu0 <- c(3,3)
mu1 <- c(6,8)
Sigma0 <- cbind(c(1,2),c(2,20))
Sigma1 <- cbind(c(2,3),c(3,16))


# The following sets up the grid of points on which the densities will
# be evaluated.
x1.grid <-  rep(seq(0,8,l=40),48)
x2.grid <- rep(seq(-5,18,l=48),rep(40,48))
x.grid <- cbind(x1.grid,x2.grid)

# Here are the grid points
plot(x.grid)

# just a handy function to capture the pathname for the files we will load.
web441 <- function(x)
         {paste('http://www.undergrad.math.uwaterloo.ca/~stat441/R-code/',
                 x,
                 sep='')
         }

# get code for creating densities
source(web441('dmvnorm.R')) 

# Create the class 0 density
den0 <- dmvnorm(x.grid,mu0,Sigma0)
# Plot its contours (6 of them)
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
# identify the mean as 0
points(mu0[1],mu0[2], pch="0")

# Repeat for the class 1 density
#
den1 <- dmvnorm(x.grid,mu1,Sigma1)
# Following fools R to believe that the graphics device is new and
# hence doesn't need clearing
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
#
# Now produce the discriminant boundary which in this case will be
# quadratic.  Notice how the boundary falls where the two densities
# are equal.

par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den1),40,48),nlevels=1,lwd=3,col='black')

# now let's do a similar thing with three densities...
# the third density will have a different mu
mu2 <- c(1,10)
# and the same Sigma as the first density
Sigma2 <- Sigma0
# Or, maybe a very different one (try each in turn)
Sigma2 <- cbind(c(1,-2),c(-2,10))
den2 <- dmvnorm(x.grid,mu2,Sigma2)

# draw all three contours
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)

# Could do the classifications pairwise

# den0 versus den1 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den1),40,48),nlevels=1,lwd=4,col='purple',pch="")
  
# den0 versus den2 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den2),40,48),nlevels=1,lwd=4,col='orange',pch="")
  
# den1 versus den2 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den1)-log(den2),40,48),nlevels=1,lwd=4,col='steelblue',pch="")
  
# Or each versus the other two by max likelihood
# redraw all three contours each time

# den1 versus den0 OR den2
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den1)-pmax(log(den2),log(den0)),40,48),nlevels=1,lwd=4,col='blue',pch="1")
  
  
# den2 versus den0 OR den1
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den2)-pmax(log(den1),log(den0)),40,48),nlevels=1,lwd=4,col='green',pch="2")


# den0 versus den1 OR den2 
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-pmax(log(den2),log(den1)),40,48),nlevels=1,lwd=4,col='red',pch="0")


#
# This is a slightly better way of getting the class boundaries.  At
# each point a label for the most probable class, and then used to determine the
# colour of each point in the grid.  
#
# First redraw the contours 
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
# Now find the best.  The value stored in best will be a unique, but arbitrary,
# label for each density.  The numbers chosen here are 2 for den0, 4 for den1
# and 3 for den2.  When used as the value of the col argument in points, 2,4, and 3
# stand for the colours red, blue and green, respectively (which will match the colours
# of the contour plots).  (Why these numbers?  white=0, black=1, red=2, green=3, blue=4.
# Hint: think flatbed 4-colour pen plotters.)

best <- apply(cbind(den0,den1,den2),1,function(x){c(2,4,3)[x==max(x)]})
points(x1.grid,x2.grid,col=best,pch=19,cex=2)
# The border (based on these grid points) can be determined by the contour function
# by giving it g-1 levels where g is the number of groups (here 3).
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(best,40,48),nlevels=2,
        lwd=3,col='black')


######################
#
#
#   Two densities, hyperbolic (still quadratic, just not parabolic)
#   discriminant function
#
#

# The following choices will produce hyperbolic regions 
# Note that the principal axes are as before but that the
# scale parameters are opposite in size.
vecs0 <- eigen(Sigma0)$vectors
Sigma0 <- vecs0 %*% cbind(c(1,0),c(0,10)) %*% t(vecs0)
vecs1 <- eigen(Sigma1)$vectors
Sigma1 <- vecs1 %*% cbind(c(20,0),c(0,1)) %*% t(vecs1)


# Create the class 0 density
den0 <- dmvnorm(x.grid,mu0,Sigma0)
# Plot its contours (6 of them)
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
# identify the mean as 0
points(mu0[1],mu0[2], pch="0")

# Repeat for the class 1 density
#
den1 <- dmvnorm(x.grid,mu1,Sigma1)
# Following fools R to believe that the graphics device is new and
# hence doesn't need clearing
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
#
# Now produce the discriminant boundary which in this case will still be
# quadratic.  However, this time it is hyperbolic in shape rather than
# parabolic (levels=c(0) was added to avoid a glitch in R (or S) which
# would add a level curve at -100; don't know why).
#
# Again, notice how the boundary falls where the two densities
# are equal.

par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-log(den1),40,48),
	        levels=c(0),nlevels=1,lwd=3,col='black')

# now let's do a similar thing with three densities... third one, den2
# is unchanged from before.

# draw all three contours
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)

# Could do the classifications pairwise

# den0 versus den1 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-log(den1),40,48),
        levels=c(0), nlevels=1,lwd=4,col='purple',pch="")
  
# den0 versus den2 

contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)

par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-log(den2),40,48),
        levels=c(0),
        nlevels=1,lwd=4,col='black',pch="")
  
# den1 versus den2 

contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
  
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den1)-log(den2),40,48),
	        levels=c(0),
        nlevels=1,lwd=4,col='black',pch="")
  
# Or each versus the other two by max likelihood
# redraw all three contours each time

# den0 versus den1 OR den2 
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-pmax(log(den2),log(den1)),40,48),
	        levels=c(0),
        nlevels=1,lwd=4,col='black',pch="0")


# den1 versus den0 OR den2
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den1)-pmax(log(den2),log(den0)),40,48),nlevels=1,lwd=4,col='black',pch="1")
  
  
# den2 versus den0 OR den1
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den2)-pmax(log(den1),log(den0)),40,48),nlevels=1,lwd=4,col='black',pch="2")


# The three regions
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den1)-pmax(log(den2),log(den0)),40,48),nlevels=1,
	        lwd=4,col='black',pch="1")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den2)-pmax(log(den1),log(den0)),40,48),
        nlevels=1,lwd=4,col='black',pch="2")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-pmax(log(den2),log(den1)),40,48),
	        levels=c(0),
        nlevels=1,lwd=4,col='black',pch="0")
  


#
# This is a slightly better way of getting the class boundaries.  At
# each point a label for the most probable class, and then used to determine the
# colour of each point in the grid.  
#
# First redraw the contours 
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,40,48),nlevels=6,col='green',lwd=3)
# Now find the best.  The value stored in best will be a unique, but arbitrary,
# label for each density.  The numbers chosen here are 2 for den0, 4 for den1
# and 3 for den2.  When used as the value of the col argument in points, 2,4, and 3
# stand for the colours red, blue and green, respectively (which will match the colours
# of the contour plots).  (Why these numbers?  white=0, black=1, red=2, green=3, blue=4.
# Hint: think flatbed 4-colour pen plotters.)

best <- apply(cbind(den0,den1,den2),1,function(x){c(2,4,3)[x==max(x)]})
points(x1.grid,x2.grid,col=best)
# The border (based on these grid points) can be determined by the contour function
# by giving it g levels where g is the number of groups (here 3).
# (produces a double border everywhere rather than just at one or two borders).
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(best,40,48),nlevels=3,
        lwd=3,col='black')

# Or filling the points and making them larger
points(x1.grid,x2.grid,col=best,pch=19,cex=2)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(best,40,48),nlevels=3,
        lwd=3,col='black')

######################
#
#  Introducing a mixture model (assuming an equal mixture)
#  for the three densities.

contour(unique(x1.grid),unique(x2.grid),
  matrix( (den0 + den1 +den2)/3,40,48),nlevels=8,col='purple',lwd=3)

#  Overlay the discriminant functions based on generalized likelihood ratio
#

# The three regions
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den1)-pmax(log(den2),log(den0)),40,48),nlevels=1,
	        lwd=4,col='black',pch="1")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den2)-pmax(log(den1),log(den0)),40,48),
        nlevels=1,lwd=4,col='black',pch="2")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-pmax(log(den2),log(den1)),40,48),
	        levels=c(0),
        nlevels=1,lwd=4,col='black',pch="0")


#
#  Or if the likelihood ratio has to be at least 10
#

contour(unique(x1.grid),unique(x2.grid),
  matrix( (den0 + den1 +den2)/3,40,48),nlevels=8,col='purple',lwd=3)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den1)-pmax(log(den2),log(den0)),40,48),nlevels=1,
	        levels=c(log(10)),
	        lwd=4,col='black',pch="1")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den2)-pmax(log(den1),log(den0)),40,48),
	        levels=c(log(10)),
        nlevels=1,lwd=4,col='black',pch="2")
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(log(den0)-pmax(log(den2),log(den1)),40,48),
	        levels=c(log(10)),
        nlevels=1,lwd=4,col='black',pch="0")