nfl2014 <- read.csv("http://www.stat.ufl.edu/~winner/sta4702/nfl_combine_2014.csv")
attach(nfl2014); names(nfl2014)

Height_OL <- Height[Position=="OL"]
Weight_OL <- Weight[Position=="OL"]
HandLen_OL <- HandLen[Position=="OL"]
ArmLen_OL <- ArmLen[Position=="OL"]

nfl2014_OL <- data.frame(Height_OL, Weight_OL, HandLen_OL, ArmLen_OL)

detach(nfl2014)
attach(nfl2014_OL)

####  Data now only contains players who are offensive linemen (OL)

(n_OL <- length(Height_OL))     ### n=sample size (number of players)
one_OL <- rep(1,n_OL)           ### nx1 vecor of 1s

### create X matrix by column binding the p=4 variables
X_OL <- cbind(Height_OL, Weight_OL, HandLen_OL, ArmLen_OL)

### Obtain S_n and R_n using built-in R functions cov and cor
(S_n <- ((n_OL-1)/n_OL) * cov(nfl2014_OL))
(R_n <- cor(nfl2014_OL))

### Obtain px1 vector xbar = (1/n) X'1_n   * is scalar multiply, %*% is matrix
(xbar_OL <- (1/n_OL) * t(X_OL) %*% one_OL)

### Repeats each of p=4 means n times
xbar1_OL <- rep(xbar_OL[1,], n_OL)
xbar2_OL <- rep(xbar_OL[2,], n_OL)
xbar3_OL <- rep(xbar_OL[3,], n_OL)
xbar4_OL <- rep(xbar_OL[4,], n_OL)

### Create p=4 deviation vectors, each of length n:  d_ i = x_i - xbar_i1_n
d1_OL <- X_OL[,1] - xbar1_OL
d2_OL <- X_OL[,2] - xbar2_OL
d3_OL <- X_OL[,3] - xbar3_OL
d4_OL <- X_OL[,4] - xbar4_OL

### Shows that xbar_i1_n is perpendicular to d_i
t(xbar1_OL) %*% d1_OL
t(xbar2_OL) %*% d2_OL
t(xbar3_OL) %*% d3_OL
t(xbar4_OL) %*% d4_OL

### J_n is nxn matrix of 1s,   I_n is nxn Identity matrix
J_n <- matrix(rep(1,n_OL^2),ncol=n_OL)
I_n <- diag(n_OL)

### Deviation matrix = (I - (1/n)J)X
dmat_OL <- (I_n - (1/n_OL)*J_n) %*% X_OL 

### Sum of Squares and Crossproducts matrix:  SSCP=X'(I - (1/n)J)X  and S
(SS_dmat_OL <- t(dmat_OL) %*% dmat_OL)
(S_dmat_OL <- (1/(n_OL-1)) * SS_dmat_OL)

### D^(1/2)
(D12_dmat_OL <- diag(sqrt(diag(S_dmat_OL))))

### R = D^(-1/2)SD^(-1/2)
(R_dmat_OL <- solve(D12_dmat_OL) %*% S_dmat_OL %*% solve(D12_dmat_OL))

### Obtain eigenvalues (lambda) and eigenvectors (as matrix P) of S
eigen_OL <- eigen(S_dmat_OL)
(lambda_OL <- eigen_OL$val)
(P_OL <- eigen_OL$vec)      ### P is [e1 ... ep]  (pxp)

### Obtain diagonal matrix of eigenvalues and sqrt(eigvals)
(Lambda_OL <- diag(lambda_OL))
(Lambda12_OL <- diag(sqrt(lambda_OL)))

### Reproduce S, S^{-1} and S^(1/2) from spectral decompositon
(S_sd_OL <- P_OL %*% Lambda_OL %*% t(P_OL))
(Sinv_sd_OL <- P_OL %*% solve(Lambda_OL) %*% t(P_OL))
(S12_sd_OL <- P_OL %*% Lambda12_OL %*% t(P_OL))

S_sd_OL %*% Sinv_sd_OL
S12_sd_OL %*% S12_sd_OL