主成分分析:正射影によるデータ次元の圧縮
簡単な解説とプログラム(実行形式ファイルとソースコードファイル)
なお、主成分分析の説明を<岡本安晴「心理学データ分析と測定」>で行っている。
データがp個の変量からなるとき、その第i番目のデータを
で表わす。N個のデータ、・・・、はp次元空間内のN個の点で表わされる。主成分分析(Principal Component Analysis; PCA)は、これらのp次元空間におけるデータに含まれる情報をできるだけ解り易くするために低次元空間にを正射影するものである(Okamoto (2006), “A Justification of Rotation in Principal Component Analysis: Projective viewpoint of PCA”)。図1は、2変量のデータを2次元空間(平面)において表わしたものである。点A、点B、点Cが正射影Pによって直線(1次元空間)上に射影されて点AP、点BP、点CPとなっている。
図1
一般に、p次元空間の点のq次元空間への正射影は階数がqであるp次の正方行列P(正射影行列)を右から掛けることによってと表わすことができる(射影行列については、<岡本安晴(2009)「統計学を学ぶための数学入門[下]」培風館;Pp. 65-70>参照)。
N個のデータをまとめて
と表わすとき、データを正射影したものは
と表わされる。Xに含まれる情報をできるだけ含む正射影PをXPの多変量の分散
(1)
を最大にするものとして選ぶ。ここで、データXは各変量が平均0、分散1であるように標準化されているものとし、は相関行列
である。
式(1)を最大にする正射影Pは次の特異値分解
(2)
から構成することができる。q次元空間への正射影の場合は、特異値を大きいものからq個選び、対応するVの部分行列を
と置くとき、式(1)を最大にする正射影は
で与えられる。は、で張られるq次元空間への正射影である。
の座標系に関する座標は
で与えられる。とはに対応する特異値と特異ベクトルである。これらの座標を標準化したものは主成分(principal component)と呼ばれ、次式で与えられる。
この主成分Cに対するパターン行列は
である。ただし、である。
いま、射影された空間XPにおける座標系をXPの他の座標系に変えたときのXPの座標に対するパターン行列をとおき、パターン行列の変換を
(3)
とおくと、
が成り立つ。このパターン行列の変換(3)を回転という。
ここでは、主成分分析をデータに含まれる情報をできるだけ保持した低次元空間への正射影として説明した。射影された先の低次元空間における座標軸は任意であるので、正射影したデータの解釈に合わせて座標軸を設定することができる。この座標軸の再設定は、式(3)の回転によって表される。これに対して、主成分を1つ1つ独自の意味をもつものとして求める方法がある。この1つ1つ順番に求める方法により得られた主成分は、それぞれ独自の意味をもつものとして固定される。正射影による方法の場合は、式(1)を最大にする正射影すなわちによる正射影の射影の空間が意味をもち、その中の座標系は任意である。
なお、因子分析のウェブサイトで説明した次式
(3)
を最小にするものとして与えられる解は、主成分分析における式(1)を最大にする解と数学的には同じ式で与えられるが、数学的には同じ形式であっても実質科学的には異なる概念に基づくものであることに注意する必要がある。一般に、数式の表面的な類似性から実質科学的解釈を行うことは危険である。物理学において、数式の表面的類似に基づいて議論を展開することの危険が指摘されている(サイフェ(著)・林大(訳)「宇宙を復号する」早川書房,2007、p. 100)。主成分分析における主成分は射影によって得られる(別の観点からいえば重み付き)合成量であるが、因子分析の場合の因子は潜在変数として想定されているものであるという実質科学的に異なる性質の量であることに注意しなければならない。
プログラムpca.exeは、正射影によってデータを低次元に圧縮して主成分(principal components)を求め、さらに射影された空間内において座標軸の回転を行い、データに含まれる情報の読み取りを容易にすることを試みるものである。プログラムは.NET Framework3.5対応の実行形式のファイルであるので、ファイル名pca.exeをクリックしてダウンロード・実行することができる。実行すると図1のフォームが表示される。
図1
ラジオボタンは、「PrintData」の方にチェックが入っていると、読み込んだデータがファイルに出力される。出力されたファイルを見ることによって、データが正しく読み込まれているかどうか確認することができる。射影された空間における座標値(主成分/成分値)は「PrintCom」にチェックが入っているとファイルに出力される。初期解は主成分を与える。主成分の数はリッチテキストボックスに設定することができる。座標軸(主成分)の回転は、直交回転としてバリマックス回転、斜交回転としてプロマックス回転を行っている。rはプロマックス回転におけるパラメータである。通常はでよいが、他の値を試すこともできる。
「GO」ボタンをクリックすると、図1のフォームに設定された条件で主成分分析(データの最適な正射影の計算)が行われる。まず、入力データファイル名の設定を求めるダイアログウィンドウが表示される(図2)。
図2
入力データファイルは、以下のような形式で用意する。
======================= 入力データファイルの例 ===========================
/
6
ID
P-1
P-2
P-3
P-4
P-5
P-6
1
46
57
53 65
54 67
2
19
16
38
34
26
40
3
28
30
23
20
35
18
4
40
44
55
48
43
45
5 64
62 67
57
71 70
6
49
43
61
61
52
59
7
27
36
57
64
30
50
8
52
48
46
52
44
50
9
66
65
57
57
63
48
10
54
52
26
34
55
26
11
62
58
55
65
54
65
12
81
83
54
51 84 54
13
32
36
19
33
36
26
14
81
80
34
27
69
31
15
32
41
71
56
41
58
16
70
70
29
24 66 30
17
48
45
35
34
46
38
18
56
49
46
46
49
35
19
66
65
71
72
66
57
20
55
64
60 69
63 47
/
======================== 入力データファイルの例 =========================
データは、スラッシュ「/」で始まる行で挟む。最初のスラッシュ「/」で始まる行の次の行に、主成分分析を行うデータの変量の数pを書く。次の行に各変量のラベルを空白で区切って書く。データはその次の行から、1行に1ケースのデータを書いていく。各行の第1変量はデータのケース番号を書く。このケース番号は分析には用いられないが、データの入力を解りやすくするための番号である。また、ケース番号は文字列として読み込まれているので、英字などを含んだ「abc123」のようなものでもよい。ケース番号は区切り文字を含まないものであるが、区切り文字はプログラムにおいて次のコードで指定されている。
array<wchar_t> ^ separator = { ' ', '\t', ',',' ','/' };
詳しい説明は<岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房、p. 209 >を参照されたい。
各ケースのデータの変量の値は2番目から書く。データの最終行の次の行は、スラッシュ「/」で始まる行を置く。この2つのスラッシュで始まる行で挟まれた部分以外は、プログラムでは無視されるので、コメントなどを書くことができる。
図2の入力ファイル名を設定するダイアログウィンドウでファイル名を設定して「開く」ボタンをクリックすると、次に出力ファイル名の設定を求めるダイアログウィンドウが表示される(図3)。
図3
適当なテキストファイル名を設定して「保存」ボタンをクリックする。「保存」ボタンをクリックすると計算が始まり、計算が終了すると計算結果を出力ファイルに保存したことを知らせるダイアログウィンドウが表示される(図4)。
図4
図4の「OK」ボタンのクリックでプログラムの実行が終了する。実行後、出力ファイルを開くと以下のようになっている。
まず、入力ファイル名が出力されている。
Data File =
F:\yasuharu\www\LaCoocan\VCpp\pca\pca\pca\DataV6.txt
入力ファイル名の出力に続いて、スラッシュで始まる行まで、読み込んだ行がそのまま出力される。
次に、読み込んだ変量名と読み込んだケース数nが出力される。
Vars =
P-1
P-2
P-3
P-4 P-5 P-6
n = 20
続いて、各変量の平均値と標準偏差が出力される。
mean sd
P-1
51.400 17.261
P-2
52.200 16.345
P-3
47.850 15.551
P-4
48.450 15.628
P-5
52.350 14.728
P-6
45.700 14.488
Eigen valueとして、特異値の2乗の値と、その累積和が出力される。
Eigen values...
Squared lambda(1) = 3.45665 Cum. sum =
3.45665
Squared lambda(2) = 2.18338 Cum. sum =
5.64003
Squared lambda(3) = 0.15401 Cum. sum =
5.79404
Squared lambda(4) = 0.11771 Cum. sum =
5.91175
Squared lambda(5) = 0.05112 Cum. sum =
5.96287
Squared lambda(6) = 0.03713 Cum. sum =
6.00000
続いてパターン行列(と主成分:主成分の出力を指定した場合)が出力される。
パターン行列
P.C.-1
P.C.-2
P-1
0.76686 -0.61656
P-2 0.81292 -0.55673
P-3
0.72398
0.63147
P-4
0.68197
0.66907
P-5
0.83447 -0.51828
P-6
0.72260
0.61503
特異値分解とパターン行列および主成分の計算は以下のコードによって行っている。
SVD( x, rnk, U, V, lambda ); // 特異値分解
fout->WriteLine();
fout->WriteLine("Eigen values...");
double
sum_lambda2 = 0.0;
for (int i = 0; i < rnk; i++){
sum_lambda2 += sqr(lambda[i][0]);
fout->WriteLine("Squared lambda(" + (i+1).ToString() + ") = " + (sqr(lambda[i][0])).ToString("F5")
+ " Cum. sum =
" + sum_lambda2.ToString("F5"));
}
myMat Pat(p,q); // パターン行列
myMat C(n,q); // 主成分行列
for (int i = 0; i < p; i++)
for (int j = 0; j < q; j++)
Pat[i][j] = V[i][j] * lambda[j][0]; // パターン行列
for (int i = 0; i < n; i++)
for (int j = 0; j < q; j++)
C[i][j] = U[i][j] * Math::Sqrt(n); // 主成分得点(座標値)
これらの初期値としての主成分の出力後、直交回転としてバリマックス(VM)回転、斜交回転としてプロマックス(PM)回転を行ったときのパターン行列(と座標値/成分値)が算出・出力される。
バリマックス回転後のパターン行列
P.C.(VM)-1 P.C.(VM)-2
P-1
0.04715
0.98285
P-2
0.12248
0.97765
P-3
0.95277
0.12304
P-4
0.95305
0.06666
P-5
0.16555
0.96827
P-6
0.93955
0.13291
プロマックス回転後のパターン行列
P.C.(PM)-1 P.C.(PM)-2
P-1 -0.06582 0.99665
P-2
0.01159
0.98263
P-3
0.95714
0.01527
P-4
0.96396 -0.04225
P-5
0.05659
0.96809
P-6
0.94252
0.02686
成分間相関行列...
P.C.(PM)-1
P.C.(PM)-2
P.C.(PM)-1
1.00000
0.22419
P.C.(PM)-2
0.22419
1.00000
プログラムについての詳しい内容はソースコードを参照されたい。ソースコードファイルなどは圧縮ファイルpca.zipとしてまとめた。圧縮ファイル名pca.zipをクリックしてダウンロードすることができる。ダウンロードした圧縮ファイルは、マウスの右ボタンクリックで表示されるメニュから「展開」あるいは「解凍」を選んで解凍することができる。解凍したファイルはVisual C++2015で開くことができる。
岡本安晴 作