library(MASS)

set.seed(770116)

n <- 200
rho <- 0.5

z1 <- rnorm(n)
z2 <- rnorm(n) + rho*z1
f1 <- z1 - 2
f2 <- z2 + 2
c <- "group1"
for (i in 2:n) c <- c(c,"group1") 

z1 <- rnorm(n)
z2 <- rnorm(n) + rho*z1
f1 <- c(f1, z1)
f2 <- c(f2, z2 - 2)
for (i in 1:n) c <- c(c,"group2") 

z1 <- rnorm(n)
z2 <- rnorm(n) + rho*z1
f1 <- c(f1, z1 + 2)
f2 <- c(f2, z2 + 2)
for (i in 1:n) c <- c(c,"group3") 

m.group1 <- cbind(mean(f1[c="group1"]),mean(f2[c="group1"]))
m.group2 <- cbind(mean(f1[c="group2"]),mean(f2[c="group2"]))
m.group3 <- cbind(mean(f1[c="group3"]),mean(f2[c="group3"]))

f  <- cbind(f1,f2)

fit <- lda(f,c) 
print(fit$prior) 

size <- 100

f1.axis <- seq(from=min(f1),to=max(f1),length=100)
f2.axis <- seq(from=min(f2),to=max(f2),length=100)

grid <- mat.or.vec(size*size,2)
for (i in 1:size) {
  for (j in 1:size) {
    grid[size*(j-1)+i,1] <- f1.axis[i];
    grid[size*(j-1)+i,2] <- f2.axis[j];  
  }
}

pred <- predict(fit,grid)

chat <- pred$class

edge <- mat.or.vec(size*size,1)
edge <- as.logical(edge)
for (i in 1:size) { 
  for (j in 1:size) {
    cntr <- size*(j-1)+i
    left <- cntr - size
    down <- cntr - 1
    if ( (left > 0) && (down > 0) ) {
      edge[cntr] <- (chat[left] != chat[cntr]) || (chat[down] != chat[cntr])
    }
    if (i == 1) edge[cntr] <- FALSE
  }
}

source("psopts.r")
postscript(file="lda_rule.eps")
par.default <- par(no.readonly=TRUE) 
plot(grid,type='n',xlab="feature 1",ylab="feature 2")
# points(grid[chat=="group1",],pch='.',col="red")
# points(grid[chat=="group2",],pch='.',col="green")
# points(grid[chat=="group3",],pch='.',col="blue")
# points(m.group1,pch="o",col="red")
# points(m.group2,pch="o",col="green")
# points(m.group3,pch="o",col="blue")
  points(f[c=="group1",],pch='o',col="red")
  points(f[c=="group2",],pch='o',col="green")
  points(f[c=="group3",],pch='o',col="blue")
  points(grid[edge,],pch='+',col="black",cex=0.4)
dev.off()

ps.options(width=7.0,height=3.0);
postscript(file="lda_data.eps")
par(mfrow=c(1,3),mar=c(5.5,4,3,2)+0.1) # mar=c(bottom,left,top,right)
plot(f,type='n',xlab="feature 1",ylab="feature 2")
  points(f[c=="group1",],pch='o',col="red")
plot(f,type='n',xlab="feature 1",ylab="feature 2")
  points(f[c=="group2",],pch='o',col="green")
plot(f,type='n',xlab="feature 1",ylab="feature 2")
  points(f[c=="group3",],pch='o',col="blue")
dev.off()