주성분분석 : 코드 (1)

자료 만들기

library(MASS)
# 공분산행렬 제작
D = diag(c(1, 2, 0.4))
data = c(1, 0.6, 0.1, 0.6, 1, 0.2, 0.1, 0.2, 1)
rr = matrix(data, nrow=3, byrow=T)
Sigma = D^{1/2} %*% rr %*% D^{1/2}

# 평균벡터
mu = rep(0, 3)

# 자료행렬 d1 제작
n = 100
d1 = mvrnorm(n, mu, Sigma)

# 표준화된 자료행렬 d2 제작
d2 = scale(d1)

# d1의 성분 확인
cov(d1); cor(d1)

이후, 다음 작업을 실시.

eeR = eigen(cov(d2))

cov(d2) = 표준화된 데이터의 공분산행렬 = $\Sigma$
eeR = 공분산행렬를 분해시킨 뒤, 각 요소를 담고 있는 배열
- eeR$values : 고유값을 의미
- eeR$vectors : 고유벡터를 의미

왜 이렇게 하나?

f(\bold{X}) \approx \exp(Y^\prime \cdot S^{-1} \cdot Y) = \exp(Y^\prime \cdot U \cdot \Lambda^{-1}\cdot U^{\prime} \cdot Y)

X가 자료행렬이라고 하자. 다변량 정규분포를 따르므로, 그 구조는 근사 시 위의 식과 같다.
- exp 속 $U^{\prime} \cdot Y$ = 성분들을 담은 행렬
- 여기서 몇 개를 골라내야 하며, 그 기준은 설명력이 큰 것 = 분산이 큰 것
따라서 $U^{\prime} \cdot Y$ 의 분산을 알아야 한다!

\begin{aligned} Cov(U^{\prime} \cdot Y) &= Y^{\prime} \cdot Cov(Y) \cdot Y \\ &= Y^{\prime} \cdot U \cdot \Lambda \cdot U^{\prime} \cdot Y \end{aligned}

여기에서 각 성분의 설명력을 결정짓는 것은 \Lambda 이다.
- 대각선에 있는 각 람다 = 오름차순으로 있음.
따라서 표준화자료의 공분산을 찾고, 그것의 고유값을 알아야 한다.
- $Y = X - \underline{\mu}$ = 표준화된 자료 = d2
- 코드에서 d2의 공분산행렬을 구하고, 그것을 분해시킨 것 = eeR을 구한 이유.

핵심 아이디어 : 성분을 직접 구하지 않고, 성분의 설명력만을 가지고 차원 축소를 진행한다.

성분과 그 점수

이렇게 구해진 고유값을 이용해서 각 주성분의 분산을 생각해볼 수 있다

> ## 성분의 분산
> l1 = eeR$values; l1
[1] 1.7943800 0.8887445 0.3168755

> ## 각 성분의 점수 계산
> l1/sum(l1) * 100
[1] 59.81267 29.62482 10.56252

> ## 각 성분의 누적 점수 계산
> cumsum(l1/sum(l1) * 100)
[1]  59.81267  89.43748 100.00000

예를 들어, 두 번째 성분의 경우

첫번째 성분의 분산 : 0.8887445
첫번째 성분의 점수 : 29.62482
첫번째 성분의 누적점수 : 89.43748

eeR 분석하기

eeR들의 각 요소들을 모아봤을 때 cov(d2)가 만들어지지 보자.

> U %*% L %*% t(U)
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6830738 0.2062988
[2,] 0.6830738 1.0000000 0.2142176
[3,] 0.2062988 0.2142176 1.0000000

> cov(d2)
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6830738 0.2062988
[2,] 0.6830738 1.0000000 0.2142176
[3,] 0.2062988 0.2142176 1.0000000

모두 같다는 것을 알 수 있다.

또햔, cov(d2)의 고유벡터가 정규직교기저를 이루는지 확인해보자.

> round(U %*% t(U), 2)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

단위행렬이 나오는 것을 보고, 잘 나오는 것을 확인할 수 있다.

주성분의 점수 (1)

d2를 행렬로 바꿔서 Y라고 하고, 여기서 첫번째 행만을 뽑아내서 이것을 Y11라고 하자.

첫번째 subject의 첫번째 주성분 점수를 뽑아 낼 수 있다.

Y = as.matrix(d2)
Y11 = matrix(Y[1, ], ncol = 1) # 첫번째 subject
p11 = t(U[, 1]) %*% Y11; p11

주성분의 점수 (2)

이번에는 모든 subject의 첫번째 주성분 점수들을 한번에 도출해보자.

p1 = Y %*% U[, 1]
p1[1:3]
mean(p1)
sd(p1)

시각화 (1) : 주성분 점수

첫번째 주성분의 시각화

p1 = 첫번째 주성분 점수들을 모은 벡터.

이때 p1의 분포를 그리면 거의 정규분포에 가깝다.

par(mfrow = c(1, 2))
hist(p1)
qqplot(qnorm(ppoints(length(p1)), sd = sqrt(L[1, 1])), p1,
       xlab = 'Normal dist', ylab = 'emphirical dist')
qqline(p1, distribution = function(p) qnorm(p, sd = sqrt(L[1, 1])))

정규분포를 따르는 것을 확인 가능!

1st 주성분 & 2nd 주성분의 결합분포

3차원 그래프를 이용해서 두 주성분 점수벡터(p1, p2)를 동시에 시각화하자.

p1 = Y %*% U[, 1]
p2 = Y %*% U[, 2]
require(KernSmooth)
par(mfrow = c(1, 2))
z1 = bkde2D(cbind(p1, p2), 1)
persp(z1$fhat, theta = 0, phi = 90)
z2 = bkde2D(Y[, 1:2], 1)
persp(z2$fhat, theta = 0, phi = 90)

1st 주성분 & 2nd 주성분의 biplot

par(mfrow = c(1, 1))
plot(p1, p2, cex = .5, pch = 19, col = "gray", 
     ylim = c(-4, 4), xlim = c(4, -4), bty = 'l')
abline(h = 0, v = 0, lty = 3)

시각화 (2) : 주성분 점수

4x4 공분산행렬

D <- diag(c(1, 2, 0.4, 0.5))
data <- c(1, 0.6, 0.1, 0.3,
          0.6, 1, 0.2, 0.4,
          0.1, 0.2, 1, 0.5,
          0.3, 0.4, 0.5, 1)
rr <- matrix(data, nrow=4, byrow=TRUE)
Sigma = D^{1/2} %*% rr %*% D^{1/2}
mu = rep(0, 4)
n = 100
d3 = mvrnorm(n, mu, Sigma)
d4 = scale(d3)
eeR = eigen(cov(d4))

biplot 모음

P = d3 %*% eeR$vectors[, 1:2]
P = data.frame(P)

display = function(d, x, y)
{
  t1 = paste(x, "번째 var")
  t2 = paste(y, "번째 var")
  plot(d[, x], d[, y], pch = 19, col='grey', 
       cex = 1, xlab = t1, ylab = t2, bty = 'l')
}

par(mfrow = c(2, 2))
display(d3, 1, 2)
display(d3, 1, 3)
display(d3, 3, 4)
plot(P$X1, P$X2, pch = 19, col = "grey", 
     cex = 1, xlab = "1st pc",
     ylab = "2nd pc", bty = 'l')

3차원 그래프

library(plotly)
P = d3 %*% eeR$vectors[, 1:3]
P = data.frame(P)
plot_ly(x = P$X1, y = P$X2, z = P$X3, type = "scatter3d",
        width = 8, height = 9, mode = 'markers', 
        marker = list(size = 5))

고유값의 합

> p = 4
> sum(eeR$values > 0) == p
[1] TRUE
> eeR$values
[1] 1.8308114 1.1107489 0.6377428 0.4206968

고유값 = 데이터의 상관 행렬에서 주성분에 대한 분산의 크기
따라서 이를 체크하기 위한 코드