実装関連事項

各種プログラミング言語の基本的な書き方やソフトウェア等の使用方法について.

主成分分析 (principal component analysis) とは多次元のデータを低次元データに縮約する方法のことである.PCA とも呼ばれる.高次元データを2次元か3次元に落とすことで人間が理解しやすい形式に変換するために行う.R で主成分分析を行う関数には princomp() と prcomp() の2種類が存在するが,princomp() にはサンプルサイズと変数の数に制限があるので,prcomp() を使えば,その点を意識せずに済む.以下に示されるような R に付属の11次元からなるデータ mtcars に対して主成分分析を行い次元の縮約を行う.mtcars の次元は,燃費 mpg,シリンダーの数 cyl,排気量 displacement,馬力 hp,重量 wt,ドラッグレースのタイム qsec,V型エンジンか直列型エンジンかの分類 vs,トランスミッション am,ギア数 gear,キャブレーター数 carb の11個で,それに対して,26個の車のデータが格納されている.

                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1

このデータは以下のコマンドで確認できる.

$ mtcars

主成分分析の実行

これに対して PCA を行うには以下のように打つ.これによりオブジェクト rpca に主成分分析の結果が格納される.ここで,scale=T は数値のスケールが合っていないときに用いるオプションであり,スケールが統一されている場合は入力する必要はない.スケールが合っていないとは,次元1の単位が cm であるときに,次元2の単位が kg であるようなときのことをいう.ここで T を指定すると相関行列から,F を指定すると分散共分散行列から主成分分析が行われる.基本的に単位が合っていない解析の方が多いだろうから,分散共分散行列から主成分分析を行うのは限られた場合のみである.

$ rpca=prcomp(x=mtcars,scale=T)

以上のコマンドにて主成分分析が終了し,rpca には主成分スコア $x,標準偏差 $sdev (固有値の正の平方根),$rotation (各主成分軸の固有ベクトル),$center (解析に用いた各項目の値の平均) 等が格納される.多次元データを低次元に縮約したいだけなら,これらの項目の中では主成分スコア,$x だけを考えれば良い.以下のように打つことで,それぞれ PC1軸および PC2軸の値を確認できる.

$ rpca$[,1]
$ rpca$[,2]

主成分得点のプロット

これらを以下のコマンドで散布図上にプロットする.ここで,ラベルをプロットエリア上に描くためのコマンド pointLabel() を使用するための maptools は外部パッケージなので,install.packages() で事前にインストールする必要がある.

$ library(maptools)
$ par(las=1)
$ plot(x=NULL,type="n",xlab="PC1",ylab="PC2",xlim=c(-5,5),ylim=c(-5,5),xaxs="i",yaxs="i",xaxt="n",yaxt="n",bty="n")
$ axis(side=1,at=seq(-5,5,1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-5,-4,-3,-2,-1,0,1,2,3,4,5))
$ axis(side=2,at=seq(-5,5,1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-5,-4,-3,-2,-1,0,1,2,3,4,5))
$ points(x=rpca$x[,1],y=rpca$x[,2],pch=16,col="#ff8c00")
$ pointLabel(x=rpca$x[,1],y=rpca$x[,2],labels=rownames(mtcars),cex=0.8)
$ box(bty="l")

これを実行した結果は以下のようになる.Mazda RX4 と Mazda RX4 Wag のように名前が似ていて,車も似ていると思われる車が近いし,Honda Civic や Toyota Corolla のような類似する大衆車が近い位置にプロットされているのがわかる.

ida_r_pca_01.svg

因子負荷量のプロット

主成分分析において新たに生じた軸,PC1 と PC2 であるが,これが何を意味するのかの解釈は観測者に委ねられる.解釈のヒントは元の軸を散布図上にプロットすることで得られる.以下のように打つ.

$ library(maptools)
$ par(las=1)
$ plot(x=NULL,type="n",xlab="PC1",ylab="PC2",xlim=c(-0.5,0.5),ylim=c(-0.5,0.5),xaxs="i",yaxs="i",xaxt="n",yaxt="n",bty="n")
$ axis(side=1,at=seq(-0.5,0.5,0.1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5))
$ axis(side=2,at=seq(-0.5,0.5,0.1),tck=1.0,lty="dotted",lwd=0.5,col="#dddddd",labels=expression(-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5))
$ for(i in 1:11)
$ {
$ 	arrows(0,0,rpca$rotation[i,1],rpca$rotation[i,2],col="#ff8c00",length=0.1)
$ }
$ pointLabel(x=rpca$rotation[,1],y=rpca$rotation[,2],labels=rownames(rpca$rotation),cex=1)
$ box(bty="l")

これを実行した結果は以下のようになる.これを確認したところ,PC1軸方向には馬力や燃費重さ等が分布している.一方で,PC2軸方向には速さとかギア数が分布している.なので,PC1軸は車のパワー,PC2軸は操作性を示す軸というように解釈できるかもしれない.

ida_r_pca_02.svg

上の散布図と主成分得点の散布図を一度に簡単に描くコマンドが用意されている.以下のように打つと,同じ平面上にそれらふたつの情報をプロットできる.ただし,図の装飾ができないので以下は確認用に留め,見せる用には上のコマンドを用いてちゃんと描画すべき.

$ biplot(x=rpca)
ida_r_pca_03.svg

軸の寄与率

主成分分析を用いるとこのように,高次元のデータを低次元に落とし込むことができる.このとき,上のようなたった2次元で元の11次元のどれほどの割合を説明できるのかどうかを説明する値がある.以下のようにコマンドを打つとそれを確認できる.この値を寄与率という.

$ summary(rpca)$importance

これを実行した結果は以下のようになる.2行目の Proportion of Variance が各軸の寄与率を表しており,PC1軸は,元の11次元の0.600760,すなわち約60%を説明できることとなる.3行目には累積寄与率が示されており,このとき,PC1軸とPC2軸を合わせると,全体の約84%を説明できているということになる.

                            PC1      PC2       PC3       PC4       PC5       PC6       PC7      PC8       PC9      PC10      PC11
Standard deviation     2.570681 1.628026 0.7919579 0.5192277 0.4727061 0.4599958 0.3677798 0.350573 0.2775728 0.2281128 0.1484736
Proportion of Variance 0.600760 0.240950 0.0570200 0.0245100 0.0203100 0.0192400 0.0123000 0.011170 0.0070000 0.0047300 0.0020000
Cumulative Proportion  0.600760 0.841720 0.8987300 0.9232400 0.9435600 0.9627900 0.9750900 0.986260 0.9932700 0.9980000 1.0000000

図示すると以下のようになる.

ida_r_pca_04.svg

主成分得点からデータの逆算

主成分分析では,データを主成分空間上にプロットするが,これとは逆に主成分空間中のあるデータを示す点 (または任意の点) から本来のデータそのものを復元 (または新規導出) することもできる.主成分分析とは,1行a列からなるデータ X の主成分スコア S を (正規化した) X と軸の固有ベクトルA (a行1列) を掛け合わせることで求める.なので,元のデータ X はその主成分スコアにAの逆行列を掛けることで復元することができる.実際にはこの掛け算の結果を,正規化した分だけ補正する.計算は,主成分分析に使用した各変数の標準偏差 (rpca$scale) を各々の要素に掛け (対角行列化した標準偏差を掛けることと同義),さらに各変数の平均値 (rpca$center) を加える.式で表すと以下のようになる.ここで t(A) はAの転置行列を表す.相関行列から主成分分析を行った場合,すなわち,prcomp() において scale=T で計算した場合は rpca$scale を対角行列化したものを使用する必要がある.

X= S * t(A) %*% diag(rpca$scale) + rpca$center

分散共分散行列から行った場合 (scale=F) は以下の式を考える.

X= S * t(A) + rpca$center

上で行なった PCA の結果から Mazda RX4 の値を復元する.最初に以下のコマンドで固有ベクトルAの準備をする.固有ベクトルは主成分分析の結果を格納したオブジェクト rpca の $rotation に格納されているので,それを取り出すには以下のように打つ.

$ eigen01=as.matrix(rpca$rotation[,1])
$ eigen02=as.matrix(rpca$rotation[,2])
$ eigen03=as.matrix(rpca$rotation[,3])
$ eigen04=as.matrix(rpca$rotation[,4])
$ eigen05=as.matrix(rpca$rotation[,5])
$ eigen06=as.matrix(rpca$rotation[,6])
$ eigen07=as.matrix(rpca$rotation[,7])
$ eigen08=as.matrix(rpca$rotation[,8])
$ eigen09=as.matrix(rpca$rotation[,9])
$ eigen10=as.matrix(rpca$rotation[,10])
$ eigen11=as.matrix(rpca$rotation[,11])

この情報を用いて,元のデータの1行目のデータ Mazda Rx4 の値を主成分スコアから逆算するには以下のように打つ.

$ (rpca$x[1,1]*t(eigen01) + rpca$x[1,2]*t(eigen02) + rpca$x[1,3]*t(eigen03) + rpca$x[1,4]*t(eigen04) + rpca$x[1,5]*t(eigen05) + rpca$x[1,6]*t(eigen06) + rpca$x[1,7]*t(eigen07) + rpca$x[1,8]*t(eigen08) + rpca$x[1,9]*t(eigen09) + rpca$x[1,10]*t(eigen10) + rpca$x[1,11]*t(eigen11) ) %*% diag(rpca$scale) + rpca$center

これを実行した結果は以下のようになる.元のデータを完全に再現できていることが分かる.

     [,1] [,2] [,3] [,4] [,5] [,6]  [,7]         [,8] [,9] [,10] [,11]
[1,]   21    6  160  110  3.9 2.62 16.46 3.330669e-16    1     4     4

以上では,全11個の固有ベクトルを用いて逆算したので,完全に元のデータと一致する結果が得られたが,以下のようにすると,主成分1軸および2軸のみからデータを逆算することができる.

$ (rpca$x[1,1]*t(eigen01) + rpca$x[1,2]*t(eigen02)) %*% diag(rpca$scale) + rpca$center

これを実行した結果は以下のようになる.主成分1軸と2軸の累積寄与率が80%を超えており,高い部類なので,再現率は悪くない.

         [,1]     [,2]    [,3]     [,4]     [,5]     [,6]    [,7]      [,8]      [,9]    [,10]    [,11]
[1,] 21.66999 5.888978 190.762 161.1851 3.949175 2.759129 16.6661 0.3380028 0.8480916 4.368929 3.729912

以上では,現実の主成分スコアからデータを逆算したが,全く架空の主成分スコアから元のデータを逆算することもできる.以下の図,赤色の点で示される主成分スコアが PC1=4,PC2=-4 からなる架空のデータを考える.

ida_r_pca_05.svg

これに対して,元のデータを計算するには以下のように打つ.

$ (4*t(eigen01) + (-4)*t(eigen02)) %*% diag(rpca$scale) + rpca$center

これを実行した結果は以下のようになる.このような操作をすることで,素性のよくわかった平面上からデータを逆算できるので,売れ行きデータ等と併せることで新たな車のスペックを設計する際等に活用できるのかもしれない.

         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]      [,8]       [,9]    [,10]    [,11]
[1,] 10.96209 8.546148 437.7041 168.9767 2.379964 5.131667 19.72803 0.2865687 -0.9197895 1.712355 1.523225
Hatena Google+