정준분석 : 코드

자료 생성

Sigma = matrix(rWishart(1, pp, diag(pp)), pp)
R1 = diag(1 / sqrt(diag(Sigma))) %*% Sigma %*% diag(1 / sqrt(diag(Sigma)))
Sigma = matrix(rWishart(1, pp, diag(pp)), pp)
R2 = diag(1 / sqrt(diag(Sigma))) %*% Sigma %*% diag(1 / sqrt(diag(Sigma)))
round(R1, 2)
round(R2, 2)

그런데, 이 자료의 상관성을 보니 X, Y의 상관성이 없는 것 같음.

## 자료 뽑기 & 상관성 확인
nn = 10000
ss = bdiag(R1, R2) # 자료 합치기
pp2 = dim(ss)[1]
mu = rep(0, pp2)
d1 = mvrnorm(nn, mu, ss)
round(cor(d1), 2) # x와 y의 상관성이 없다.

따라서 새 자료를 만듬.

## 행과 열을 교체한 새 자료
d1 = d1[, as.numeric(matrix(1:pp2, 2, byrow = T))]
d1 = data.frame(d1) # 인위적인 변수 조정
names(d1) = paste('X', 1 : pp2, sep='')
round(cor(d1), 2)

정준변수점수 뽑기

library(CCA)

# 행렬 간 상관관계
X = d1[, 1:3]
Y = d1[, -c(1:3)]
print(matcor(X, Y), digits=2)

# 정준상관분석 : cc
d2 = cc(X, Y)
str(d2)

# X에 대한 정준변수점수
score = d2$scores$xscores
round(apply(score, 2, mean), 2)
round(apply(score, 2, sd), 2)
round(cov(score), 2)

손으로 만들기

ss = cov(d1)
Sxx = ss[1:3, 1:3]
Sxy = ss[1:3, 4:6]
Syy = ss[4:6, 4:6]

{\Sigma_{xx}}^{-1/2} ~\Sigma_{xy}~{\Sigma_{xx}}^{-1}~ {\Sigma_{xy}}^{T}~ {\Sigma_{xx}}^{-1/2}~

먼저 공분산행렬의 -1/2를 만든다.

# Sigma의 -1/2승 구하기
tmp = eigen(Sxx)
Sxx2 = tmp$vectors %*% diag(sqrt(tmp$values)) %*% t(tmp$vectors)
tmp = eigen(Syy)
Syy2 = tmp$vectors %*% diag(sqrt(tmp$values)) %*% t(tmp$vectors)

이후 행렬곱을 만들고, 분해해서 하는 작업.

# 연산1
tmp = solve(Sxx2) %*% Sxy %*% solve(Syy) %*% t(Sxy) %*% solve(Sxx2); tmp
tmp = eigen(tmp)
P = tmp$vectors
u = solve(Sxx2) %*% P

# 연산2
tmp = solve(Syy2) %*% t(Sxy) %*% solve(Sxx) %*% Sxy %*% solve(Syy2); tmp
tmp = eigen(tmp)
Q = tmp$vectors
v = solve(Syy2) %*% Q

라이브러리와의 비교

round(d2$xcoef, 4); round(u, 4)
round(d2$ycoef, 4); round(v, 4)

실제로 같은 값이 나옴을 알 수 있다.

응용하기

## 적정한 정준변수쌍 적정 수 산정
install.packages("CCP")
library(CCP)
p.asym(d2$cor, dim(X)[1], dim(X)[2], dim(Y)[2], tstat='Wilks')

## 공헌도
w = (d2$cor)^2; w
## 근사도
cumsum(w) / sum(w) * 100