# ================
# R Example on PCA
# ================

dat <- iris[,-5];
n <- nrow(dat)
pca <- prcomp(dat, scale.= TRUE, center=T, retx=T)
pca; summary(pca)

names(pca)
# The rotated variables - values of PCs
pca$x

# Scree Plot
plot(pca)
biplot(pca)

prcomp(~  Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
        data = dat, scale = TRUE)


# A Detailed look at PCA?
# -----------------------

dat
X0 <- apply(dat, 2, FUN=function(x){x-mean(x)}) 
cov0 <- (t(X0)%*%X0)/(n-1); cov0
cov(dat)

eigen.decomp <-  eigen(cor(dat))
names(eigen.decomp)
eigen.decomp$values
cbind(sqrt.eigenvalue=sqrt(eigen.decomp$values), sd.pc=pca$sdev)

eigen.decomp$vectors
pca$rotation









# -------------------------
# ANOTHER FUNCTION princomp
# -------------------------

princomp(USArrests, cor = TRUE) # =^= prcomp(USArrests, scale=TRUE)
## Similar, but different:
## The standard deviations differ by a factor of sqrt(49/50)

summary(pc.cr <- princomp(USArrests, cor = TRUE))
loadings(pc.cr)  ## note that blank entries are small but not zero
plot(pc.cr) # shows a screeplot.
biplot(pc.cr)

## Formula interface
princomp(~ ., data = USArrests, cor = TRUE)
# NA-handling
USArrests[1, 2] <- NA
pc.cr <- princomp(~ Murder + Assault + UrbanPop,
                  data = USArrests, na.action=na.exclude, cor = TRUE)
pc.cr$scores