library(MASS)


sigma.test<-matrix(c(8,3,5,8),2,2)
biv.norm<-mvrnorm(200,mu=rep(0, 2),Sigma=sigma.test)
biv.pca<-prcomp(biv.norm,center=FALSE)

#An array is just a high dimensional matrix.  So this is 2x2x50.  There is nothing special about the number 50, I just picked that to smooth the
#animation.
#The numbers in the sequence make transition.ar[,,1] the 2x2 identity matrix and transition.ar[,,50] the final transformed rotation matrix from
#the principal components.

transition.ar<-array(0,c(2,2,50))
transition.ar[1,1,]<-seq(1,0.7040705,length.out=50)
transition.ar[1,2,]<-seq(0,-0.7101301,length.out=50)
transition.ar[2,1,]<-seq(0,0.7101301,length.out=50)
transition.ar[2,2,]<-seq(1,0.7040705,length.out=50)

#Unlike transforming the segments by hand, we can use matrix multiplication to transform the data points.
#pre-multiplying biv.norm by the identity matrix gives you biv.norm again.  Pre-multiplying it by the
#rotation matrix gives you the same points arranged on the axes of their variation.  We interpolate between
#those two.
##This part requires that you have Imagemagick installed.
library(animation)
saveMovie({
for (i in 1:50) {
	plot(1,xlim=c(-10,10),ylim=c(-10,10),xlab="PC1",ylab="PC2",main="Transition to PCA",type="n")
	abline(coef(lm(biv.norm[,2]~biv.norm[,1]))[1],coef(lm(biv.norm[,2]~biv.norm[,1]))[2],col=2,lty=3)
	points((biv.norm%*%transition.ar[,,i])[,1],(biv.norm%*%transition.ar[,,i])[,2],pch=20)
	abline(coef(lm((biv.norm%*%transition.ar[,,i])[,2]~(biv.norm%*%transition.ar[,,i])[,1]))[1],coef(lm((biv.norm%*%transition.ar[,,i])[,2]~(biv.norm%*%transition.ar[,,i])[,1]))[2],col=2)
	Sys.sleep(ani.options("interval"))
	}
},interval = 0.03, moviename = "pcaanim.gif",outdir = getwd(), width = 600, height = 600);
#That saves the movie.  If you just want to see the plot in R you can type in the console (w/o the comment marks of course!)
#for (i in 1:50) {
#	plot(1,xlim=c(-10,10),ylim=c(-10,10),xlab="PC1",ylab="PC2",main="Transition to PCA",type="n")
#	abline(coef(lm(biv.norm[,2]~biv.norm[,1]))[1],coef(lm(biv.norm[,2]~biv.norm[,1]))[2],col=2,lty=3)
#	points((biv.norm%*%transition.ar[,,i])[,1],(biv.norm%*%transition.ar[,,i])[,2],pch=20)
#	abline(coef(lm((biv.norm%*%transition.ar[,,i])[,2]~(biv.norm%*%transition.ar[,,i])[,1]))[1],coef(lm((biv.norm%*%transition.ar[,,i])[,2]~(biv.norm%*%transition.ar[,,i])[,1]))[2],col=2)
#	Sys.sleep(0.03)
#	}