# 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")