#This R script illustrates how R does PCA via the pcomp function
#and relates the output to the SVD.  The R package rgl is used to 
#give nice 3d graphics.
#
#Please run the code in sections
#Section 1
set.seed(1);#To enable replication of randomly generated data
x=rnorm(100,sd=1); #Generating 100 values of x centered at 0 with standard deviation 1
y=rnorm(100,sd=.75);#Generating 100 values of y centered at 1 with standard deviation .75
z=.1*x-.1*y+rnorm(100,sd=.1);#Generating 100 values of z so that (x,y,z) roughly satisfy z=.1*x-.1*y
#
library(rgl); #Tool that allows 3d plot
#
max=max(x,y,z);#Max absolute value of data to set window size in 3dplot
plot3d(x,y,z,col="yellow",type="s", radius=0.075, xlim=c(-max,max), ylim=c(-max, max), zlim=c(-max,max));#Plots data in R^3
centroid = c(mean(x), mean(y), mean(z));#Determines center of the data
spheres3d(centroid[1], centroid[2], centroid[3], radius=0.1, col="red");#Plots the center of the data as a red dot
#
#Render the plane 
a = -0.1;
b = 0.1;
c = 1;
d = 0;
planes3d(a,b,c,d, alpha=0.4);
#
#
#End of Section 1: You should be able to see a 3Dplot of the 100 data-points, along with the plane that 
#was used to generate these "random" points.
#
#


#Start Section 2. 
#xyz is a 100 by 3 matrix of data representing 100 datapoints, each with 3 measurements
xyz=matrix(c(x,y,z),byrow=FALSE,ncol=3);
#
#Shift data to have means zero
xnew=x-mean(x);
ynew=y-mean(y);
znew=z-mean(z);
#
#N is the 100 by 3 matrix recording the 100 data points after each factor has 
#been shifted to have mean equal to 0, and divided by sqrt{n-1} to
#make into variances.  I have done the scaling because that is want prcomp does.
#
N=(1/99)^{.5}*matrix(c(xnew,ynew,znew),byrow=FALSE,ncol=3);
#
#Compute the svd (singular value decomposiiton of N)
#S$v gives the 3 by 3 matrix of rotations
#S$d gives the 3 singular values or equivalently the sdev's of the prcomp
#We can see this by noting that S$v= pr$rotation, an S$d=pr$sdev
S=svd(N);
pr = prcomp(xyz, center=TRUE,scale.=FALSE)
S$v
pr$rotation
S$d
pr$sdev
#
#
#End Section 2. In this seciton, we've seen that the sdev and rotations for the data-set 
#correspond to the singular values and columns of the svd for the scaled data-set.

#Start section 3
# We can check that each rotation is a eigenvector of N^TN with the  corresponding 
# eigenvalue the  square of the corresponding pr$sdev
#
u=pr$rotation[,1]
t(N)%*%N%*%u
(pr$sdev[1])^(2)*u
#
#
v=pr$rotation[,2]
t(N)%*%N%*%v
(pr$sdev[2])^(2)*v
#
#
w=pr$rotation[,3]
t(N)%*%N%*%w
(pr$sdev[3])^(2)*w
#
#
#End Section 3.


#Section 4. 
# Here we will add the new rotation axes to the data set
# The green axis points in the direction of the most variance among the data
# The blue axis is orthogonal to the green and subject to that points in the
# direction of the next most variance.
# The red axis is orthogonal to the blue and the green and points in the direction
# of the least variance among the data
#render the eigenvectors - move so the vectors start at the centroid of the data
#plot3d(x,y,z,col="yellow",type="s", radius=0.075, xlim=c(-max,max), ylim=c(-max, max), zlim=c(-max,max));#Plots data in R^3
#First eigenvector
segments3d(c(centroid[1], centroid[1] + pr$rotation[1,1]), c(centroid[2], centroid[2] + pr$rotation[2,1]), c(centroid[3], centroid[3] +pr$rotation[3,1]), col="green", lwd=4);
#Second eigenvector
segments3d(c(centroid[1], centroid[1] + pr$rotation[1,2]), c(centroid[2], centroid[2] + pr$rotation[2,2]), c(centroid[3], centroid[3] + pr$rotation[3,2]), col="blue", lwd=4);
#Third eigenvector
segments3d(c(centroid[1], centroid[1] + pr$rotation[1,3]), c(centroid[2], centroid[2] + pr$rotation[2,3]), c(centroid[3], centroid[3] + pr$rotation[3,3]), col="red", lwd=4);
#
#For fun, move the plot around until the red axis points directly into the screen, 
# the green axis points east, and the blue axis points north.
# What you are seing is the projection of the data-points onto the PC1-PC2 plane

#Start Section 5.
#In this section we illustrate how the PC1 vs PC2 plot is obtained
#Multiplying the data-matrix N by the rotation matrix gives the projections
#of the data points onto the various columns of the rotation matrix.
P=N%*%pr$rotation
#X stores coordinates of the projections of the scaled data-points onto the PC1-PC2 plane
X=pr$sdev[1]*P[,1]
#Y stores the coordinates of the projections of the sclaed data-points onto the PC1-PC2 plane  
Y=pr$sdev[2]*P[,2]
#The first plot gives the PC1 vs PC2 using the scaled SVD
#Compare it to what you saw in Section 4 in the 3D-plot
plot(X,Y,xlab="PC1", ylab = "PC2", main = "PC1 / PC2 - plot using the  scaled SVD")
#
#We can use R to do all of this for us. pr$x[,i] simply records the projections of 
#the original (unscalled) data onto the PCi axis. 
#End Section 5

#Start Section 6
#The plot below will give the PC1 vs PC2 plot, but with the scaling on the axis
#reflecting the original data.
plot(pr$x[,1],pr$x[,2], xlab="PC1", ylab = "PC2", main = "PC1 / PC2 - plot using prcomp function")