#  This file demonstrates linear discriminant function
#  theory when using a mixture with known probabilities.
# 
#      No data is involved, just the densities.
#      Only two densities:
#          Linear discriminant function: common covariance matrix
#          Quadratic discriminant function: different covariance matrices
#
#  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,8,l=40),48)
x2.grid <- rep(seq(-5,18,l=48),rep(40,48))
x.grid <- cbind(x1.grid,x2.grid)

# 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  
# (in class we just loaded this file since the 
#  machine wasn't connected to the internet).
source(web441('dmvnorm.R'))

# First specify the multivariate normals parameters for each 
# density

mu0 <- c(3,3)
mu1 <- c(6,8)
Sigma0 <- cbind(c(1,2),c(2,20))
# To demonstrate effect on linear discriminant we use
Sigma1 <- Sigma0
# Or for quadratic we use:  (TRY EACH ONE ON ALL THAT FOLLOWS)
Sigma1 <- cbind(c(2,3),c(3,16))

# Create the class 0 density
den0 <- dmvnorm(x.grid,mu0,Sigma0)
# and the class 1 density
#
den1 <- dmvnorm(x.grid,mu1,Sigma1)

#  Show the density assuming an equal mixture
#   First, overlaid the others
contour(unique(x1.grid),unique(x2.grid),matrix(den0,40,48),nlevels=6,col='red')
# identify the mean as 0
points(mu0[1],mu0[2], pch="0")
# Repeat for the class 1 density
par(new=T)
contour(unique(x1.grid),unique(x2.grid),matrix(den1,40,48),nlevels=6,col='blue')
#  Now overlay the mixture
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(0.5 * (den0 + den1),40,48),nlevels=8,lwd=2, col='brown')

# Now have the mixture on its own with the discriminant function overlaid.
#
contour(unique(x1.grid),unique(x2.grid),
  matrix(0.5 * (den0 + den1),40,48),nlevels=8,lwd=1, col='brown')
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den1),40,48),nlevels=1,lwd=3,col='black')

# What if the data are more likely to come from the class 0
# than from the class 1?
#
# If they are twice as likely from class 0 than class 1
# then this could be described as a different mixture.
# say as

contour(unique(x1.grid),unique(x2.grid),
  matrix(((2/3) * den0 + (1/3)* den1),40,48),nlevels=10,lwd=1, col='brown')

# Overlaying the discriminant function as before
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(log(den0)-log(den1),40,48),nlevels=1,lwd=3,col='black')

#
# This can be accommodated by making use of these probabilities
# as in

par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0),
  matrix(log(den0)-log(den1)+log(2),40,48),nlevels=1,lwd=2,col='green')
  
#
# Similarly, had class 0 been half as likely to occur as class 1
#

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

#
# This can be accommodated by making use of these probabilities
# as in

par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0),
  matrix(log(den0)-log(den1)+log(1/2),40,48),nlevels=1,lwd=2,col='green')
  
# If they are 10 times more likely from class 0 than class 1
# then this could be described as a different mixture.
# say overlaid as follows:
par(new=T)
contour(unique(x1.grid),unique(x2.grid),
  matrix(((10/11) * den0 + (1/11)* den1),40,48),nlevels=30,lwd=1, col='brown')

# note that you can also modify the prior probabilities, which has the
# has the effect of shifting the line one way or the other.  We add to
# the log of the ratio the log of the ratio of prior probs.  For
# example, if class 0 is twice as likely as class 1, we do:
par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0),
  matrix(log(den0)-log(den1)+log(2),40,48),nlevels=1,lwd=2,col='green')


# Or, if class 0 is half as likely as class 1, we do:
par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0),
  matrix(log(den0)-log(den1)+log(1/2),40,48),nlevels=1,lwd=2,col='orange')
  
# Or, if class 0 is 10 times as likely as class 1, we do:
par(new=T)
contour(unique(x1.grid),unique(x2.grid),levels=c(0),
  matrix(log(den0)-log(den1)+log(10),40,48),nlevels=1,lwd=2,col='blue')