5番目の資料

主成分分析(PCA, Principal Component Analysis)
多くの変量の値を、できるだけ情報の損失を少なくして少数個の総合的指標(主成分)で代表させる。
１．考え方
ｎ人について、科目毎のｐ個の得点が得られたとき、全科目の総合得点として合理的な指標を求
めたい。まず、２変量(p=2)の場合を考えると、重み係数a 1 , a 2 を持ちこんで次のｚを考える。
z = a1 ( x1 − x1m ) + a2 ( x2 − x2 m ) （ただし、 a12 + a 22 = 1 とする）・・・・（１）
係数 a1 , a 2 はｚが x1 ， x 2 をよく代表するように決めなければならない。
「よく代表する」ことを客観的に表現するために、ｚの分散を最大化するように a1 , a 2 を決める。
そこで変量 x1 ， x 2 の特徴をみるために平均、分散、共分散を用いる。
平均：
x1m =
分散：
V
11
=
x2
1 n
1 n
，
= ∑ x2 i
∑
x
x
1i
2m
n i =1
n i =1
2
1 n
，
−
(
x
x
)
∑
1
i
1
m
n i =1
共分散：
V 12 = V 21 =
V
22
=
・・(2)
2
1 n
−
(
)
x
x
∑
2
2
i
m
n i =1
θ2
1 n
∑ ( − )( − )
n i =1 x1i x1m x2i x2 m
V ( z) ≡
=
H
θ1
x1 cos θ1
O
このとき、合成変量ｚの分散 V(z)は(3)式となる。
x2 cosθ 2
( x1 , x2 )
z
x1
2
1 n
−
(
)
z
z
∑
i
m
n i =1
2
1 n
∑{ ( − ) + ( − )}
n i =1 a1 x1i x1m a 2 x 2i x 2m
= a1 V 11 + 2 a1 a2V 12 + a2V 22
2
2
主成分分析のイメージ図
・・・・(3)
図のように各点をプロットしたときに、ちらばりの最も大きい方向 OZ の方向を見つけ、その
方向の成分を総合得点とすれば良い。そこで、任意の点 P の OZ 軸上の座標ｚ（P から OZ に下
ろした垂線の足 H までの原点 O からの距離 OH）は、 z = x1 cos θ1 + x 2 cos θ 2 となり、(1)式と比
較して a1 , a 2 は方向余弦ともなっている。すなわち、 a12 + a2 2 = 1
・・・・(4) の条件の制約の
もとで、(3)式の V(z)を最大にするような係数 a1 , a 2 を求める最適化問題に帰着する。
２．Lagrange の未定乗数法による解法
このような制約条件付き最大化問題は、Lagrange の未定乗数λを用いて次のように定式化する。
F (a1 , a 2 , λ ) = a1 V11 + 2a1 a 2V12 + a 2 V22 − λ (a1 + a 2 − 1) ・・・・・・・(5)
2
2
2
2
(5)式を a1 , a 2 , λ に関して最大にする問題となる。そこで、F を a1 , a 2 , λ で偏微分する。
∂F
= 2a1V11 + 2a2V12 − 2λa1 = 0
∂a1
，
∂F
2
2
= −( a1 + a 2 − 1) = 0
∂λ
∂F
= 2a1V12 + 2a 2V22 − 2λa 2 = 0
∂a 2
・・・(6)
51
(6)式を変形させて
(V11 − λ )a1 + V12 a 2 = 0
V V
a
a
V −λ
V12  a1 
または、  11 12  1  = λ  1  、または  11
  = 0 ・・(7)
V V  a 
a 
 V
−
V
λ
12
22
2
2
22
12

 
 
 a 2 

V12 a1 + (V22 − λ )a 2 = 0
(7)式は、分散共分散行列の固有値問題となっている。
これよりλに関する２次方程式を得る。
(V11 − λ )(V22 − λ ) − V12 = 0 ，すなわち、 λ2 − (V11 + V22 )λ + (V11V22 − V12 ) = 0
2
2
・・・(8)
2
2
(8)式の解 λ は、 λ = V11 + V22 ± (V11 − V22 ) + 4V12 ・・・(9) の２つの実数解である。
2
λ1 + λ 2 = V11 + V22 ・・・(9)' が成り立つ。
この解を λ1 ， λ2 （ λ1 ≧ λ2 ≧0）とし、(7)式に代入して a1 , a 2 の比を求め、(4)式の条件から
解と係数の関係より
、
（ a12 , a 22 ）が例えば次のように求まる。
２つのλに対応して、２組の（ a11 , a 21 ）
a
11
=
V
V + (λ1−V 11)
，
12
a
2
2
21
=
λ1−V
V + (λ1−V 11)
11
12
2
2
12
a
=
λ 2 −V
V + (λ 2 −V 22)
2
12
，
22
2
a
=
22
V
V + (λ 2 −V 22)
・・(10) を得る。
12
2
2
12
12
(7)式、(10)式より
V11a1i + 2V12 a1i a2i + V22 a2i = λi (a1i + a2i ) = λi
2
2
2
2
(i=1,2)・・(11) となる（下記［参考］参照）。
３．主成分と寄与率
(11)式の左辺は(3)式よりV(z)の分散であり、(4)式を使えばV(z i )=λ i が成り立つ。(９)式の２つの
解λ i (i=1,2) のうちの大きい固有解λ 1 に対応する固有ベクトル ( a11 , a 21 ) が求める係数になり、
これを用いた合成変量： z1 = a11 ( x1 − x1m ) + a21 ( x2 − x2 m ) を第１主成分と言う。λ 1 が十分大きく
ない場合は第２主成分も考える。これは２番目に大きい解λ 2 に対する固有ベクトル ( a12 , a 22 ) を
用いた合成変量
z2 = a12 ( x1 − x1m ) + a22 ( x2 − x2 m )
を採用する。z 2 の分散はλ 2 である。
第１主成分の寄与率は、 λ1 /(V11 + V22 ) となる。変量の数が２以上でも同様に考える。
合成変量ｚを用いた総合得点の順位付けの場合、最大固有値の場合の固有ベクトルを用いて主成
分得点を算出し評価する。
主成分得点： z1i = a11 ( x1i − x1m ) + a21 ( x2i − x2 m )
・・・(12)
○単位固有ベクトルに固有値の平方根を掛けたものを主成分負荷量(または因子負荷量)と呼ぶ。
［注意］変量x 1 、x 2 の変域、単位が異なる場合などは、その影響が出ないように各変数を基準化
し、その分散共分散行列は相関行列となるので、これをもとに解析を進める。
［参考］
V11a1 + 2V12 a1a2 + V22 a2
2
2
V11a1 + 2V12 a1a2 + V22 a2 = λ1 (a1 + a2 ) = λ の導出
= V11a1 ⋅ a1 + 2V12 a1a2 + V22 a2 ⋅ a2
(7)式より V11a1 = λa1 − V12 a2 ， V22 a2 = λa2 − V12 a1
= λa1 − V12 a1a2 + 2V12 a1a2 + λa2 − V12 a1a2
2
2
2
2
= (λa1 − V12 a2 )a1 + 2V12 a1a2 + (λa2 − V12 a1 )a2
2
2
= λ (a1 + a2 ) = λ
これらを代入すると、右のとおりとなる。
2
52
2
【演習 21】
（下記ファイル：http://www.coins.tsukuba.ac.jp/~fukui/da/PCA.xls ）
８人について、英語(x1)、数学(x2)の成績が表１のように出た。両者を総合して順位を決めるのに、
それぞれの点数に重みをつけた以下の式によって合理的な順位付けを行いたい。
主成分分析によって総合得点指標ｚ＝a1*(x1-x1m)＋a2*(x2-x2m) を求め、それによって順位付けよ。
そのときの主成分の寄与率も求めよ。　また、Ｒを使っても分析せよ。
表．並び替え用リスト
氏名
A
B
C
D
E
F
G
H
平均
x1i
77
75
68
55
46
62
51
48
x1m
60.3
分散 V11
146.8
共分散 V12
7.75
x2i 合計
90 167
79 154
77 145
88 143
89 135
72 134
74 125
75 123
x2m
80.5
n=
a=
V22
54.0
b=
c=
λ＝
zi
8
1
a11=V12/sqrt(V12^2+(λ-V11)^2)
=
a21=(λ-V11)/sqrt(V12^2+(λ-V11)^2)
=
R=λ/(V11+V22)
＝
手順
氏名
Ｚｉ
A
１．英語、数学それぞれの得点の分散 V11，V22 を求める。
（分散は各得点から平均値を引いた偏差平方和を自由度で割る） B
C
２．英語と数学の得点の共分散 V12 を求める
（共分散は各得点から平均値を引いた偏差積和を自由度で割る） D
３．分散共分散行列の固有方程式（次式）の解を求める
E
F
λ2－(V11+V22)λ＋(V11V22-V12*V12)＝0
この方程式の係数をa=1, b=-(V11+V22), c=V11*V22-V12^2 とおくと G
大きいほうの解λは、λ＝（-b + sqrt(b^2-4*a*c))/(2*a)
H
４．この解の大きい方のλを次の式に代入して重み係数ａ11、ａ21を求める
a11=V12/sqrt(V12^2+(λ-V11)^2)
a21=(λ-V11)/sqrt(V12^2+(λ-V11)^2)
５．寄与率Rは、R=λ/(V11+V22)となる
６．大きい方のλ値から導いたa11, a21を用いて、各人の主成分得点値を
ｚi＝a11*(x1i-x1m)＋a21*(x2i-x2m) から求めて左表の右端の欄に記入する
７．ｚi値を、右表の並び替え用リストに、値のみコピーする
（値のみコピーは、貼り付けのとき、［形式を選択して貼り付け］から「値」だけを選択）
８．右表内の氏名、データ全体を選択し、メニューの[データ］から
［並び替え］を用いて、得点の高い順にソートする
単なる合計の得点順と、順位が異なっているところに注目する。
【解法】ソルバーを使って、分散が最大となる係数（固有ベクトル、主成分負荷量）を求める
「データ」メニューから「分析」に「ソルバー」が表示されないときは、次のステップでアクティブにすること
１．エクセル画面の左上の Office マーク
をクリックして、
［Excel のオプション］を選択
２．
「アドイン」を選択して出るリストの中から「ソルバーアドイン」を選択して「設定」を選択
エクセル表（左下図）で係数（固有ベクトル、主成分負荷量に比例）a 1 , a 2 を適当な初期値を与え、
その平方和を求めるセルを設定する。また、各学生の総合指標となるｚの値を初期値のa 1 , a 2 を使
って求めて、ｚの分散を求めるセルを設定する。
その後、ソルバーの設定ウィンドウ（右下図）で、目的セル（左下図ではD19）が最大値をとる
ように目標値に設定し、制約条件として、平方和（D6 セル）が 1 とするように設定し、変動させ
る変数を（B6:B7）に設定を行った後「実行」で、係数a 1 , a 2 の値が求まる。
初期値
=sumsq(B6:C6) 平方和
=sumproduct(B$6:C$6,B11:C11) の積和計算
=varp(B11:B18)
=varp(D11:D18) で(母集団の)分散計算
53
R で主成分分析を使う場合
> English <- c(77,75,68,55,46,62,51,48)
> Math <- c(90,79,77,88,89,72,74,75)
> pcompo <- cbind(English,Math)
> rownames(pcompo) <- c("A","B","C","D","E","F","G","H")
> result <- princomp(pcompo,cor=F)
> summary(result)
Importance of components:
Comp.1
Comp.2
Standard deviation
11.3653279 6.8203241
Proportion of Variance 0.7352298 0.2647702
Cumulative Proportion 0.7352298 1.0000000
> loadings(result)
Loadings:
Comp.1 Comp.2
English -0.996
Math
-0.996
SS loadings
Proportion Var
Cumulative Var
> plot(result)
> biplot(result)
>
Comp.1 Comp.2
1.0
1.0
0.5
0.5
0.5
1.0
【演習 22】
ソルバーを用いて第一主成分を求める
20 名による各野菜の好みを、主成分負荷量（固有ベクトル）
a3
a4
a2
a1
1.000
1.000
1.000
1.000
1(最低)~10(最高)段階で評価
した結果の表がある。
茄子
人参
大根キャベツ
y
u
v
x
第一主成分得点を求め、その
4
4
5
A
3
寄与率を求めよ。
8
4
6
7
B
C
D
E
F
G
H
I
J
K
L
M
N
O
P
Q
R
S
T
3
7
6
9
6
10
7
7
10
6
10
5
7
7
7
2
8
5
3
9
8
7
8
10
9
4
7
5
8
8
8
5
5
2
4
9
5
7
5
6
6
10
8
8
6
5
10
5
8
8
9
1
8
5
2
9
7
6
6
8
5
4
7
7
10
10
7
9
5
3
6
7
分散（固有値）
全分散（全固有値の和）　＝　寄与率　＝(主成分pの分散）/（全分散）＝
54
a5
1.000
牛蒡
w
5
6
1
7
8
8
9
10
5
5
7
7
8
7
8
7
6
3
4
9
平方和
5.00
主成分得点
p
R の主成分分析結果の説明（例)
column bind の意味
> x1 <- c(86,71,42,62,96,39,50,78,51,89)
rbind (row bind)だと
> x2 <- c(79,75,43,58,97,33,53,66,44,92)
       
> x3 <- c(67,78,39,98,61,45,64,52,76,93)
       
 ( x1 ) 
 のように配置

  x1  x 2   x3   x 4  
> x4 <- c(68,84,44,95,63,50,72,47,72,91)
 ( x 2 )
       










 ( x3 ) 
> pc <- cbind(x1,x2,x3,x4)


 ( x 4 )
> pc


のように配置
x1 x2 x3 x4
[1,] 86 79 67 68
[2,] 71 75 78 84
principal component
[3,] 42 43 39 44
[4,] 62 58 98 95
[5,] 96 97 61 63
T では、相関行列に変換して
[6,] 39 33 45 50
から主成分分析する
[7,] 50 53 64 72
F では、分散共分散行列を用
[8,] 78 66 52 47
いて主成分分析する
[9,] 51 44 76 72
[10,] 89 92 93 91
> pr <- princomp(pc,cor=T)
> summary(pr)
Importance of components:
Comp.1
Comp.2
Comp.3
Comp.4
Standard deviation
1.6494644 1.1053504 0.22893548 0.071105876
Proportion of Variance
0.6801832 0.3054499 0.01310286 0.001264011
Cumulative Proportion 0.6801832 0.9856331 0.99873599 1.000000000
> loadings(pr)
主成分得点
Loadings: （主成分負荷量）
Comp.1 Comp.2 Comp.3 Comp.4
x1 -0.487 0.527 -0.499 0.485
主成分負荷量
x2 -0.511 0.474 0.539 -0.474
x3 -0.508 -0.481 -0.504 -0.506
Sum of Square loadings
x4 -0.493 -0.516 0.455 0.533 主成分負荷量の平方和
Comp.1 Comp.2 Comp.3 Comp.4
SS loadings
1.00 1.00 1.00 1.00
Proportion Var 0.25 0.25 0.25 0.25
Cumulative Var 0.25 0.50 0.75 1.00
> plot(pr)
寄与率
> biplot(pr)
主成分得点
主成分負荷量
累積寄与率
主成分負荷量(固有ベクトル、矢印)と主成分得点(数字)の散布図
オプション課題 (主成分分析)
３０名の相撲力士の身長(Height)、体重(Weight)、胸囲(Chest)、座高(Seating height)、体重／身
長比(Weight Height ratio)のデータが下記サイトに PCAsumo.txt のファイル名で R で使える形
式で参照できるように置いてある。このデータを使って主成分分析を行い、第一主成分は何を表
すかを適切な言葉で表現し、その寄与率を求めよ。
http://www.coins.tsukuba.ac.jp/~fukui/da/PCAsumo.txt
55

Download Report