## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----install-------------------------------------------------------------
list.of.packages=c("ggfortify")
if(length(which(!list.of.packages %in% installed.packages()))){
  install.packages(list.of.packages[!list.of.packages %in% installed.packages()])
}

## ----Example 5.1---------------------------------------------------------
setwd("C:/Users/dujuan/Documents/workshop_MV/workshop_MV/workshopcode")

 pop.prcomp=function(A)
  {
   eigvv=eigen(A);
   Lam=eigvv$value;
   Ev=eigvv$vectors;
   R=diag(sqrt(1/diag(A)))%*%(Ev)%*%diag(sqrt(Lam))
   ord=order(Lam);
   cat("There are", nrow(A), "PCs\n\n")
   cat("The variances of PCs are\n")
   cat(rev(sort(Lam)),"\n\n")
   cat("The PC coefficients are\n")
   Pcc=Ev[,rev(ord)]
   dimnames(Pcc)=list(NULL,paste("PC",1:nrow(A),sep=""))
   print(Pcc)
   cat("\n\n")
   cat("The Correlation Matrix between Variables and PCs\n")
   dimnames(R)=list(paste("X",1:nrow(A),sep=""),paste("PC",1:nrow(A),sep=""))
   print(R)
  }

 Sigma=matrix(c(1,-2,0,-2,5,0,0,0,2),nrow=3)  
 A<-Sigma
 eigvv=eigen(A);
 Lam=eigvv$value;
 Ev=eigvv$vectors;
 Lam/sum(Lam)
 cumsum(Lam)/sum(Lam)
 pop.prcomp(Sigma)
 

## ----Example52-----------------------------------------------------------
 tract=read.table("T8-5.DAT")
 S=cov(tract)
 tract.pca=prcomp(tract,scale=F)
 tract.R=diag(sqrt(1/diag(S)))%*%(tract.pca$rotation)%*%diag(tract.pca$sdev)
 #dimnames(tract.R)=list(paste("X",1:nrow(S),sep=""),paste("PC",1:nrow(S),sep=""))
 pc=tract.pca$rotation
 na=c("Total population", "Profession", "Employment (%)", "Government employment(%)", "Medium home value")
 dimnames(tract.R)=list(paste(na,sep=""),paste("PC",1:nrow(S),sep=""))
 dimnames(pc)=list(paste(na,sep=""),paste("PC",1:nrow(S),sep=""))
 print(tract.R)
 print(pc)
 cumvar.pca=cumsum(tract.pca$sdev^2)/sum(tract.pca$sdev^2)
 screeplot(tract.pca,type="l")
 library(ggfortify)
 autoplot(tract.pca)
 

## ----Example53-----------------------------------------------------------
 stock=read.table("T8-4.DAT")
 p=dim(stock)[2]
 stock.pca=prcomp(stock,scale=T)
 print(stock.pca)
 stock.R=(stock.pca$rotation)%*%diag(stock.pca$sdev)
 dimnames(stock.R)=list(paste("X",1:p,sep=""),paste("PC",1:p,sep=""))
 print(stock.R)
 nam=c("JP Morgan", "Citibank", "Wells Fargo", "Royal Dutch Shell","ExxonMobil")
 PC<-stock.pca$rotation
 dimnames(PC)=list(paste(nam,sep=""),paste("PC",1:nrow(S),sep=""))
 print(PC)
 cumvar.pca=cumsum(stock.pca$sdev^2)/sum(stock.pca$sdev^2)
 screeplot(stock.pca,type="l")