#  This file demonstrates linear discriminant function
#  theory.  
#      No data is involved, just the densities.
#      Only two densities, common covariance matrix
#
#  Authors:
#       H.A. Chipman, 2003
#       R.W. Oldford, 2004
#
#


# The following sets up the grid of points on which the densities will
# be evaluated.
x1.grid <-  rep(seq(0,19,l=20),24)
x2.grid <- rep(seq(0,23,l=24),rep(20,24))
x.grid <- cbind(x1.grid,x2.grid)

# Here are the grid points
plot(x.grid)
# Create the class 0 density
mu0<-c(3,8)
den0 <- dpois(x.grid[,1],mu0[1]) *  dnorm(x.grid[,2],mu0[2],4) 
# Plot its contours (6 of them)
contour(unique(x1.grid),unique(x2.grid),
        matrix(den0,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red')
# identify the mean as 0
points(mu0[1],mu0[2], pch="0")

# Repeat for the class 1 density
mu1 <- c(8,15)
den1 <-  dpois(x.grid[,1],mu1[1]) *  dchisq(x.grid[,2],mu1[2]) 

par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),
        levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue')
points(mu1[1],mu1[2], pch="1")


par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0), labels=c(""),
  matrix(log(den0)-log(den1),20,24),nlevels=1,lwd=4,col='black')

mu2 <- c(15,10)
den2 <- dpois(x.grid[,1],mu2[1]) *  dnorm(x.grid[,2],mu2[2],4) 
# draw all three contours
contour(unique(x1.grid),unique(x2.grid),matrix(den0,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='green')

# Could do the classifications pairwise

# den0 versus den1 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den1),20,24),nlevels=1,lwd=2,col='purple',pch="")
  
# den0 versus den2 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den2),20,24),nlevels=1,lwd=2,col='orange',pch="")
  
# den1 versus den2 
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den1)-log(den2),20,24),nlevels=1,lwd=2,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,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='green')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den1)-pmax(log(den2),log(den0)),20,24),nlevels=1,lwd=2,col='blue',pch="1")
  
  
# den2 versus den0 OR den1
contour(unique(x1.grid),unique(x2.grid),matrix(den0,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='green')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den2)-pmax(log(den1),log(den0)),20,24),nlevels=1,lwd=2,col='green',pch="2")


# den0 versus den1 OR den2 
contour(unique(x1.grid),unique(x2.grid),matrix(den0,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='green')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-pmax(log(den2),log(den1)),20,24),nlevels=1,lwd=2,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 (a little thicker this time)
contour(unique(x1.grid),unique(x2.grid),matrix(den0,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='red',lwd=2)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='blue',lwd=2)
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den2,20,24),levels=c(0.02,0.01,0.008,0.006,0.004,0.002),col='green',lwd=2)
# 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-1 levels where g is the number of groups (here 3).
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
        matrix(best,20,24),nlevels=2,
        lwd=1,col='black')