読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

主成分分析を計算してみる

主成分分析を計算してみる。

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による出力。素晴らしい。

f:id:High_School_Student:20130420101802p:plain