R : Copyright 2002, The R Development Core Team
Version 1.5.1  (2002-06-17)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> invisible(options(echo = TRUE))
> 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) 
   group1    group2    group3 
0.3333333 0.3333333 0.3333333 
> 
> 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()
null device 
          1 
> 
> 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()
null device 
          1 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> proc.time()
[1] 26.46  0.64 27.18  0.00  0.00
>