x<-rnorm(100, 2,1)
y<-2*x+rnorm(100)
plot(x,y, xlim=c(-1,5), ylim=c(-2,12))
abline(lm(y~x))
xy=cbind(x,y)

plot.pca(xy, col="red")


plot.pca <- function(x, ...) {

   xe <- eigen(cov(x[,1:2]))

   xev <- xe$vectors

   xew <- xe$values

   lines(mean(x[,1])+c(-xev[1,1],xev[1,1])*xew[1],
         mean(x[,2])+c(-xev[2,1],xev[2,1])*xew[1], ...)

   lines(mean(x[,1])+c(-xev[1,2],xev[1,2])*xew[2],
         mean(x[,2])+c(-xev[2,2],xev[2,2])*xew[2], ...)
}