# generation of populations predicted from data

library(mva)

# total number of datapoints

N <- 1000

# total number of iterations

k <- 1

# read in data - with class information

c1<-read.table("ex06a.dat")
names(c1)<-c("x1","x2","y")

c2<-read.table("ex06b.dat")
names(c2)<-c("x1","x2","y")

# compute covariance matrices (without class info)

c1.cov<-cov(cbind(c1$x1,c1$x2))
c2.cov<-cov(cbind(c2$x1,c2$x2))

# compute the means

c1.m <- diag(c(mean(c1$x1),mean(c1$x2)))
c2.m <- diag(c(mean(c2$x1),mean(c2$x2)))

# set the covariance and mean
ones <- matrix( c(rep(1,N), rep(1,N)), ncol=2)

# rotate the data
c1.r <- chol(c1.cov)
c2.r <- chol(c2.cov)

# set up contour grid

x<-seq(-6,10,0.1)
y<-seq(-6,10,0.1)

Xcon <- data.frame(rep(x,length(y)), rep(y, each=length(x)))
names(Xcon) <- c("x1","x2")

# set up ploting parameters 

if( k > 6 )
	k<- 1

if( k == 1)
	par(mfcol=c(1,1)) 

if( k == 2 )
	par(mfcol=c(1,2))

if( k == 3 )  
	par(mfcol=c(3,1))

if( k == 4 )
	par(mfcol=c(2,2))

if( k == 5 )
	par(mfcol=c(2,3))

if( k == 6 )
	par(mfcol=c(2,3))

# loop through process i times

for(i in 1:k)
{
	# generate the data: normal distributions

	pop1 <- data.frame(rnorm(N),rnorm(N))
	names(pop1)<-c("x1","x2")
	pop2 <- data.frame(rnorm(N),rnorm(N))
	names(pop2)<-c("x1","x2")

	pop1.r <- as.matrix(pop1) %*% c1.r
	pop2.r <- as.matrix(pop2) %*% c2.r

	# shift the data
	pop1.r.s <- pop1.r + (ones %*% c1.m)
	pop2.r.s <- pop2.r + (ones %*% c2.m)

	# add in the class information

	pop1.r.s.c = data.frame(pop1.r.s, rep(1,length(pop1.r.s[,1])))
	names(pop1.r.s.c)<-c("x1","x2","y")
	pop2.r.s.c = data.frame(pop2.r.s, rep(0,length(pop2.r.s[,1])))
	names(pop2.r.s.c)<-c("x1","x2","y")

	data <- data.frame(rbind(pop1.r.s.c, pop2.r.s.c))

	g <- lm( y ~ x1 + x2 , data = data)
	g.pr <- predict(g, Xcon)

	text.x1=paste("x1=", format(g$coefficients["x1"], digits=3))
	text.x2=paste("x2=", format(g$coefficients["x2"], digits=3))
	text.i=paste("Int=", format(g$coefficients["(Intercept)"], digits=3))


	plot(data$x1,data$x2, xlim=c(-6,10), ylim=c(-6,10), type="n")
	points(data$x1,data$x2,col=(data$y+1))
	contour(x,y,matrix(g.pr,length(x),length(y)), levels=0.5, labex = 0, add=T)
	title("Regression of Response Matrix")
	mtext(paste(text.i,text.x1,text.x2))
}