主成分分析を計算してみる
主成分分析を計算してみる。
x1 <- c(0.966,0.207,-0.552,-0.931,0.587,1.55,-0.931,0.587,-1.69,0.207)
x2 <- c(0.619,-1.009,-0.358,-0.684,0.944,1.921,0.293,0.293,-1.335,-0.684)
x3 <- c(-0.474,-0.474,-1.684,0.735,-0.474,1.849,-0.474,0.735,-0.474,0.735)
ひっつけて行列にする。
X <- cbind(x1, x2, x3)
x1 x2 x3
[1,] 0.966 0.619 -0.474
[2,] 0.207 -1.009 -0.474
[3,] -0.552 -0.358 -1.684
[4,] -0.931 -0.684 0.735
[5,] 0.587 0.944 -0.474
[6,] 1.550 1.921 1.849
[7,] -0.931 0.293 -0.474
[8,] 0.587 0.293 0.735
[9,] -1.690 -1.335 -0.474
[10,] 0.207 -0.684 0.735
3ベクトルとも平均0、分散1に標準化されている。
mean(x1)
mean(x2)
mean(x3)
[1] -1.664657e-17 [1] 0 [1] 1.110223e-17
var(x1)
var(x2)
var(x3)
[1] 1.000535 [1] 1.000046 [1] 0.9998569
この10行3列の行列Xに主成分分析を適用し、3つの主成分を作成する。
最初から標準化している行列 を用いるので、相関係数行列から実施しようと、共分散行列から実施しようが結果は同じとなる。
行数（n=10）と列数（p=3）を定義しておく。
n <- nrow(X)
p <- ncol(X)
相関係数行列を求める。
cor(X)
x1 x2 x3
x1 1.0000000 0.7518879 0.4557925
x2 0.7518879 1.0000000 0.3995735
x3 0.4557925 0.3995735 1.0000000
X'X / (n-1) でも計算できる。
これは共分散行列の計算だが、標準化データなので相関係数行列と一致。
相関係数行列をRと定義
R <- t(X) %*% X / (n-1)
x1 x2 x3
x1 1.0005353 0.7521066 0.4558819
x2 0.7521066 1.0000464 0.3995542
x3 0.4558819 0.3995542 0.9998569
主成分分析は、合成変数（求める主成分）の分散を最大化する問題を解くに等しい。
これは、行列計算の性質を活かし、行列Rに対する固有値問題を解くことに等しくなる。
Eig <- eigen(R)
固有値が合成変数の分散となる。
固有値を対角要素とした対角行列を作成。
D <- diag(Eig$values)
[,1] [,2] [,3]
[1,] 2.088701 0.0000000 0.0000000
[2,] 0.000000 0.6666368 0.0000000
[3,] 0.000000 0.0000000 0.2451012
順にλ1、λ2、λ3とする。（λ1 > λ2 > λ3）主成分1、2、3の分散に対応する。
各固有値に対応する 固有ベクトル。
これは、合成変数（主成分）の係数となる。Vと定義する。
V <- Eig$vectors
[,1] [,2] [,3]
[1,] 0.6254950 -0.2845735 0.72648048
[2,] 0.6104604 -0.4013624 -0.68282230
[3,] 0.4858951 0.8705895 -0.07732934
各固有ベクトルは長さ1に基準化されている。
sum(V[,1]^2)
sum(V[,2]^2)
sum(V[,3]^2)
[1] 1 [1] 1 [1] 1
主成分スコアを計算する。Tと定義する。T = XV。
T <- X %*% V
[,1] [,2] [,3]
[1,] 0.75178892 -0.93600078 0.31576725
[2,] -0.71679136 -0.06659147 0.87600326
[3,] -1.38206543 -1.16530044 -0.02634424
[4,] -0.64275790 1.17935314 -0.26613994
[5,] 0.71312594 -0.95859021 -0.18148610
[6,] 3.04063179 0.39761389 -0.32863884
[7,] -0.63378525 -0.26532069 -0.83976615
[8,] 0.90316338 0.35523946 0.16954004
[9,] -2.10236553 0.60408863 -0.27953013
[10,] 0.06905545 0.85550848 0.56059485
各主成分の分散が固有値
var(T[,1])
var(T[,2])
var(T[,3])
[1] 2.088701 [1] 0.6666368 [1] 0.2451012
固有値問題を解いたので、RV = VD となる。
R %*% V
V %*% D
[,1] [,2] [,3]
x1 1.306472 -0.1897072 0.17806122
x2 1.275069 -0.2675629 -0.16736054
x3 1.014889 0.5803670 -0.01895351
一致。
固有値問題の性質により、対称行列の固有値（Dの対角要素）はすべて実数となり、各固有ベクトル（V）は直行となる。
直行なので V' = V^-1。
t(V)
solve(V)
[,1] [,2] [,3]
[1,] 0.6254950 0.6104604 0.48589510
[2,] -0.2845735 -0.4013624 0.87058953
[3,] 0.7264805 -0.6828223 -0.07732934
一致。
V'V = VV' = I とも書ける。
t(V) %*% V
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.447410e-17 -2.558988e-16
[2,] -1.447410e-17 1.000000e+00 -1.906841e-17
[3,] -2.558988e-16 -1.906841e-17 1.000000e+00
この性質があるので、各主成分間（Tの各列）の相関が0になる（T' T の非対角要素が0になる）。
Tの相関行列は I。Iは単位行列。
cor(T)
[,1] [,2] [,3]
[1,] 1.000000e+00 1.957105e-17 5.024524e-16
[2,] 1.957105e-17 1.000000e+00 -2.224671e-16
[3,] 5.024524e-16 -2.224671e-16 1.000000e+00
T' T の非対角要素が0になる。
t(T) %*% T
[,1] [,2] [,3]
[1,] 1.879831e+01 2.218583e-16 3.128228e-15
[2,] 2.218583e-16 5.999731e+00 -7.938528e-16
[3,] 3.128228e-15 -7.938528e-16 2.205911e+00
分散とn-1で割るとIとなる。
t(T) %*% T %*% solve(D) / (n-1)
[,1] [,2] [,3]
[1,] 1.000000e+00 3.697804e-17 1.418112e-15
[2,] 1.180203e-17 1.000000e+00 -3.598754e-16
[3,] 1.664101e-16 -1.323147e-16 1.000000e+00
スペクトル分解で考えてみる。
R V = V D から変形（VV'=I）して、 R = VDV' と書くことができる。
V %*% D %*% t(V)
[,1] [,2] [,3]
[1,] 1.0005353 0.7521066 0.4558819
[2,] 0.7521066 1.0000464 0.3995542
[3,] 0.4558819 0.3995542 0.9998569
相関係数行列（R）に一致。
スペクトル分解は行列の近似（元の行列と近似行列の各要素の差の2乗和を最小にする）と考えられる。
主成分分析で言えば、相関係数行列（R）を近似する行列を与えることができる。
V = [v1 v2 v3]とする。 v1、v2、v3、は主成分1、2、3のスコアベクトル。
R = V D V' = λ1 v1' v1 + λ2 v2' v2 + λ3 v3' v3と書けるので、幾つの固有値まで採用するか（幾つの主成分まで取り出すか）で近似の精度が決まる。
すべて取り出すと、先程示したように元の行列と完全に一致する。
第二主成分までだと。R = V D V' = λ1 v1' v1 + λ2 v2' v2。
V[,1:2] %*% D[1:2,1:2] %*% t(V[,1:2])
[,1] [,2] [,3]
[1,] 0.8711773 0.8736907 0.4696512
[2,] 0.8736907 0.8857689 0.3866123
[3,] 0.4696512 0.3866123 0.9983912
Rから遠からず、といった感じである。
Rの関数を使って主成分分析を実施してみる。prcomp()関数を使う。
princomp()という関数もあるようで、これは分散の計算の際nで割るそう。対してprcomp()はn-1で割るらしい。prcomp()の方が一般的ではなかろうか、、。
上の計算でも行列Rを定義した再にn-1で割ったので、prcomp()で実施。
pca_result <- prcomp(X)
summary(pca_result)
Importance of components:
PC1 PC2 PC3
Standard deviation 1.4452 0.8165 0.49508
Proportion of Variance 0.6961 0.2222 0.08169
Cumulative Proportion 0.6961 0.9183 1.00000
固有値。
pca_result$sdev^2
[1] 2.0887007 0.6666368 0.2451012
固有ベクトル。
pca_result$rotation
PC1 PC2 PC3
x1 -0.6254950 0.2845735 0.72648048
x2 -0.6104604 0.4013624 -0.68282230
x3 -0.4858951 -0.8705895 -0.07732934
ベクトルの向きが上の計算と反対になるのは仕方がない。
以下は、JMPによる出力。素晴らしい。
