## Compute number of Canonical Variates
k <- min(dim(X)[2],dim(Y)[2])
## Get number of observations in X and Y
N.x <- dim(X)[1]
N.y <- dim(Y)[1]
##
S.xy <- cov(X,Y)
S.xx <- var(X)
S.yx <- cov(Y,X)
S.yy <- var(Y)
A <- eigen(solve(S.xx) %*% S.xy %*% solve(S.yy) %*% S.yx)$vectors[,1:k]
B <- eigen(solve(S.yy) %*% S.yx %*% solve(S.xx) %*% S.xy)$vectors[,1:k]
R <- sqrt(eigen(solve(S.yy) %*% S.yx %*% solve(S.xx) %*% S.xy)$values[1:k])
## Singly standardized weights (SAS "raw")
A.single <- A %*% solve(sqrt(diag(diag(var(X %*% A)))))
B.single <- B %*% solve(sqrt(diag(diag(var(Y %*% B)))))
rownames(A.single) <- colnames(X)
rownames(B.single) <- colnames(Y)
## fully standardized weights
##Deviation score X,Y
X.dev <- Q(UnitVector(N.x)) %*% X
Y.dev <- Q(UnitVector(N.y)) %*% Y
##Z-score x,y
D.x <- solve(sqrt(diag(diag(var(X)))))
D.y <- solve(sqrt(diag(diag(var(Y)))))
Z.x <- X.dev %*% D.x
Z.y <- Y.dev %*% D.y
##Correlation Matrices
R.xy <- cor(X,Y)
R.xx <- cor(X)
R.yx <- cor(Y,X)
R.yy <- cor(Y)
A.s <- eigen(solve(R.xx) %*% R.xy %*% solve(R.yy) %*% R.yx)$vectors[,1:k]
B.s <- eigen(solve(R.yy) %*% R.yx %*% solve(R.xx) %*% R.xy)$vectors[,1:k]
A.fully <- A.s %*% solve(sqrt(diag(diag(var(Z.x %*% A.s)))))
B.fully <- B.s %*% solve(sqrt(diag(diag(var(Z.y %*% B.s)))))
rownames(A.fully) <- colnames(X)
rownames(B.fully) <- colnames(Y)

##grab UCLA data

mm <- read.table("http://www.ats.ucla.edu/stat/R/dae/mmreg.csv", sep = ",", header = TRUE)
attach(mm)