Skip to content

Instantly share code, notes, and snippets.

@ngopal
Created June 23, 2017 00:07
Show Gist options
  • Save ngopal/9c11ba56b13e374aeeb2a38308bf072f to your computer and use it in GitHub Desktop.
Save ngopal/9c11ba56b13e374aeeb2a38308bf072f to your computer and use it in GitHub Desktop.
####
# Eigens
#####
# How to calculate covariance matrix
# Great video: https://www.youtube.com/watch?v=9B5vEVjH2Pk
dat <- as.matrix(
cbind(c(90, 90, 60, 30, 30),
c(80, 60, 50, 40, 20),
c(40, 80, 70, 70, 70))
)
apply(dat, MARGIN = 2, sum) # sum of scores for each subject
apply(dat, MARGIN = 2, mean) # mean score for each subject
# to calculate covariance matrix, we need to do
# A*At
#
# First, let's calculate the deviation matrix
#
# a = A - [1]*A*(1/5)
# In fact, A represents scores, and ( [1]*A*(1/5) ) represents the average
# scores, since [1]*A gives us the sum of each column in A, and dividing it
# by 1/5 gives us the average for each
# This is the unity matrix
unity <- matrix(1, 5, 5) # 5 x 5 matrix
unity %*% dat / dim(unity)[1] # 5x5 matrix times dat, divided by dat column lengths of 5
a = dat - unity %*% dat / dim(unity)[1] # deviation matrix (scores - avg score)
# these give us different results, one a student covariance matrix,
# and one a subject covariance matrix
a %*% t(a) # student cov
t(a) %*% a # subject cov
subjects = t(a) %*% a
eobj = eigen(subjects)
# eigenvalue % contributions
eobj$values / sum(eobj$values)
eobj$vectors
# to calculate loadings
eobj$vectors %*% sqrt(eobj$values) #referred to as P. P%*%t(P) = correlation matrix
# In fact, X = P L P'
eobj$vectors %*% diag(eobj$values) %*% t(eobj$vectors)
# This equals t(a) %*% a
# calculating scores
# first, "center" and normalize by sd. Then multiply by eigenvectors
apply((((dat - mean(dat)) / sd(dat)) %*% eobj$vectors), 2, var)
# the above should be the same as the eigenvalues
###
svdobj <- svd(a)
svdobj$d*svdobj$d / sum(svdobj$d*svdobj$d) # same as eigenvalues
t(svdobj$v) #identical to eigenvectors
# eigen(a %*% t(a))$vectors is the same as svdobj$u
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment