# x = read.table("zip.train")
# n = nrow(x)
# v1 = var(x[,-1])
# ev1 = eigen(v1)


# i1 = seq(n)[ x[,1]==1 ]
# i4 = seq(n)[ x[,1]==4 ]
# i7 = seq(n)[ x[,1]==7 ]
# i8 = seq(n)[ x[,1]==8 ]

# xb = rbind( apply(x[i1,-1],2,"mean"),
# 	    apply(x[i4,-1],2,"mean"),
# 	    apply(x[i7,-1],2,"mean"),
# 	    apply(x[i8,-1],2,"mean") )

# v2 = var(xb)
# ev2 = eigen(v2)

# ii = c(i1,i4,i7,i8)

plot( as.matrix(x[ii,-1]) %*% ev1$vec[,1:2], type="n",
	xlab="pc1",ylab="pc2",main="Principal Components All Data" )
    points( as.matrix(x[i1,-1]) %*% ev1$vec[,1:2], pch=16, cex=.5 )
    points( as.matrix(x[i4,-1]) %*% ev1$vec[,1:2], pch=16, cex=.5, col=2 )
    points( as.matrix(x[i7,-1]) %*% ev1$vec[,1:2], pch=16, cex=.5, col=3 )
    points( as.matrix(x[i8,-1]) %*% ev1$vec[,1:2], pch=16, cex=.5, col=4 )
	legend(-5,-3,c("1","4","7","8"),col=c(1,2,3,4),lwd=rep(1,4))

  browser()
plot( as.matrix(x[ii,-1]) %*% ev2$vec[,1:2], type="n",
	xlab="pc1",ylab="pc2",main="Principal Components Four Means" )
    points( as.matrix(x[i1,-1]) %*% ev2$vec[,1:2], pch=16, cex=.5 )
    points( as.matrix(x[i4,-1]) %*% ev2$vec[,1:2], pch=16, cex=.5, col=2 )
    points( as.matrix(x[i7,-1]) %*% ev2$vec[,1:2], pch=16, cex=.5, col=3 )
    points( as.matrix(x[i8,-1]) %*% ev2$vec[,1:2], pch=16, cex=.5, col=4 )
	legend(-10,-6.5,c("1","4","7","8"),col=c(1,2,3,4),lwd=rep(1,4))