#
# first simulated data set in R
#

library(mva)
library(lattice)

# generate the data: normal distributions

grpA <- cbind(rnorm(1000),rnorm(1000))
grpB <- cbind(rnorm(1000),rnorm(1000))

ones <- matrix(c(rep(1,1000), rep(1,1000)), ncol=2 )

# set the covariance 
sigma0 <- cbind(c(4,2),c(4,9))
sigma1 <- cbind(c(4,2),c(4,9))

# rotate the data
r <- chol(sigma0)
grpA.r <- grpA %*% r

r <- chol(sigma1)
grpB.r <- grpB %*% r

# shift the data by the means

m0 <- c(-5,5)
m1 <- c(5,-5)

grpA.r.s <- grpA.r + (ones %*% diag(m0))
grpB.r.s <- grpB.r + (ones %*% diag(m1))


trellis.device("windows", bg="white")

a <- xyplot(grpA[,2] ~ grpA[,1], xlim=c(-20,20), ylim=c(-20,20))
b <- xyplot(grpA.r[,2] ~ grpA.r[,1], xlim=c(-20,20), ylim=c(-20,20))
c <- xyplot(grpA.r.s[,2] ~ grpA.r.s[,1], xlim=c(-20,20), ylim=c(-20,20))

# add the class information

grpA.r.s.pca <- prcomp(grpA.r.s)

g <- rbind(grpA.r.s, grpB.r.s)
g.c <- rbind(cbind(grpA.r.s, rep(0,1000)), cbind(grpB.r.s, rep(1,1000)))

d <- xyplot(grpA.r.s.pca$x[,2] ~ grpA.r.s.pca$x[,1], xlim=c(-20,20), ylim=c(-20,20))

e <- xyplot(g.c[,2] ~ g.c[,1], groups=g.c[,3], panel=panel.superpose, xlim=c(-20,20), ylim=c(-20,20))

g.pca <- prcomp(g)

f <- xyplot(g.pca$x[,2] ~ g.pca$x[,1], xlim=c(-20,20), ylim=c(-20,20))

print.trellis(a, split=c(1,1,3,3), more=T)
print.trellis(b, split=c(2,1,3,3), more=T)
print.trellis(c, split=c(1,2,3,3), more=T)
print.trellis(d, split=c(2,2,3,3), more=T)
print.trellis(e, split=c(1,3,3,3), more=T)
print.trellis(f, split=c(2,3,3,3), more=F)