rm(list = ls())  # clear the memory
###########################

## Example 10.1 in text

# Enter or compute covariances:
# Example 10.1:

S11 = cbind(c(1,.4), c(.4,1))
S12 = cbind(c(.5, .3), c(.6, .4))
S21 = t(S12)
S22 = cbind(c(1,.2), c(.2, 1))

p = nrow(S11)
q = nrow(S22)

# Get symmetric square roots: 
s1 = eigen(S11)
rootS11 = s1$vectors%*%diag(sqrt(s1$values))%*%t(s1$vectors)

s2 = eigen(S22)
rootS22 = s2$vectors%*%diag(sqrt(s2$values))%*%t(s2$vectors)

S12tilde = solve(rootS11, S12)%*%solve(rootS22)

s = svd(S12tilde, nu = p, nv = q)
G = s$u
H = s$v
R = cbind(diag(s$d),matrix(nrow=p, data = rep(0, p*(q-p))))

#check: 
round(S12tilde - G%*%R%*%t(H),5)


# Compute matrices Atilde and Btilde of (unnormalized) canonical coefficients
Atilde = t(solve(rootS11,G))
Btilde = t(solve(rootS22,H))
Btilde = Btilde[1:p,1:q]

# Normalize these so their rows have unit norm:
Anorms = sqrt(diag(Atilde%*%t(Atilde)))
A = solve(diag(Anorms),Atilde)
Bnorms = sqrt(diag(Btilde%*%t(Btilde)))
B = solve(diag(Bnorms),Btilde)

round(cbind(A, B),2)  # coefficients of the canonical pairs; unit norms

round(1/diag(Atilde%*%t(Atilde)),2)
round(1/diag(Btilde%*%t(Btilde)),2)

# Canonical correlations:
round(diag(R),2)



################### job satisfaction ############
# Triangular data in E10-5.DAT

L = matrix(0, 12, 12)
L[outer(1:12, 1:12, '<=')] <- scan("http://www.stat.ualberta.ca/~wiens/stat575/datasets/E10-5.DAT") 
S = L + t(L) - diag(diag(L))

rownames(S) = NULL
colnames(S) = NULL

p = 5
q = 7
S11 = S[1:p, 1:p]
S12 = S[1:p, (p+1):(p+q)]
S21 = t(S12)
S22 = S[(p+1):(p+q),(p+1):(p+q)]

# Get symmetric square roots: 
s1 = eigen(S11)
rootS11 = s1$vectors%*%diag(sqrt(s1$values))%*%t(s1$vectors)

s2 = eigen(S22)
rootS22 = s2$vectors%*%diag(sqrt(s2$values))%*%t(s2$vectors)

S12tilde = solve(rootS11, S12)%*%solve(rootS22)

s = svd(S12tilde, nu = p, nv = q)
G = -s$u
H = -s$v
R = cbind(diag(s$d),matrix(nrow=p, data = rep(0, p*(q-p))))


# Compute matrices Atilde and Btilde of (unnormalized) canonical coefficients
Atilde = t(solve(rootS11,G))
Btilde = t(solve(rootS22,H))
Btilde = Btilde[1:p,1:q]

# Normalize:
Anorms = sqrt(diag(Atilde%*%t(Atilde)))
A = solve(diag(Anorms),Atilde)
Bnorms = sqrt(diag(Btilde%*%t(Btilde)))
B = solve(diag(Bnorms),Btilde)

round(cbind(A, B),2)  # coefficients; unit norms

# Canonical correlations:
round(diag(R),2)

corUX = Atilde%*%S11/sqrt(diag(S11))
corUY = Atilde%*%S12/sqrt(diag(S22))
corVX = Btilde%*%S21/sqrt(diag(S11))
corVY = Btilde%*%S22/sqrt(diag(S22))

corU1X = corUX[1,]
corV1X = corVX[1,]
corU1Y = corUY[1,]
corV1Y = corVY[1,]

round(corU1X,2)
round(corV1X,2)
round(corU1Y,2)
round(corV1Y,2)

# LR tests
pval = vector(length = p)
chisq = vector(length = p)
df = vector(length = p)
rho = diag(R)
for (k in 0:(p-1)) {
n = 784
deg = (p-k)*(q-k)
chisquare = -(n-1-(p+q+1)/2)*log(prod((1-rho^2)[(k+1):p]))
pvalue = 1 - pchisq(chisquare, df = deg)

pval[k+1] = pvalue
df[k+1] = deg
chisq[k+1] = chisquare
}
num.sig = 0:(p-1)
round(cbind(num.sig,chisq, df, pval),2)