#-*- R -*-

# generation of populations predicted from data

library(MASS)
library(class)
library(nnet)

# total number of datapoints

N <- 200

# total number of iterations

j <- 10

# read in data - with class information

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

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

# 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)

# generate the training 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("a",length(pop1.r.s[,1])))
names(pop1.r.s.c)<-c("x1","x2","y")
pop2.r.s.c = data.frame(pop2.r.s, rep("b",length(pop2.r.s[,1])))
names(pop2.r.s.c)<-c("x1","x2","y")

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

p <- as.matrix(Pop[, -3])
tp <- Pop$y

ki = seq(100,1)
ke = rep(0,100)

for(m in 1:100)
{
    total=0
	for(i in 1:j)
	{
		# 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("a",length(pop1.r.s[,1])))
		names(pop1.r.s.c)<-c("x1","x2","y")
		pop2.r.s.c = data.frame(pop2.r.s, rep("b",length(pop2.r.s[,1])))
		names(pop2.r.s.c)<-c("x1","x2","y")

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

		pt <- as.matrix(Pop[, -3])

		Z <- knn(p, pt, tp, k = ki[m])
		zp = sum(class.ind(Z)[1:N,2]) + sum(class.ind(Z)[(N+1):(2*N),1])
		total = total + zp
	}
	ke[m]<-total/j
}
plot(ki,ke,type="b")