主成分分析と回転
2025.06改訂
主成分分析の古典的説明では、主成分は１つずつ構成され、それらの組み合わせに対する回転は不適切であるとされる。しかし、主成分分析を低次元空間への正射影として考える場合は、射影先の空間における座標系は数学的には任意である（岡本, 2006；岡本（2009）「統計学を学ぶための数学入門［下］」培風館、第5.7節「正射影の例」68−70；岡本（2014）「心理学データ分析と測定」勁草書房、第５章「情報の骨子を見る：主成分分析」67−77；岡本（2019）「いまさら聞けないPythonでデータ分析」丸善出版、第７章「主成分分析」115−129）。すなわち、実質的な解釈のための適切な座標系を求めることが可能であり、回転はこの座標系の選択のための数学的方法である。
Deisenroth, Faisal, and Ong (2020), “Mathematics for Machine Learning.” Cambridge University Press, Chapter 10においても、主成分分析が低次元空間への正射影で説明されている。
主成分分析をデータ空間の低次元空間への正射影と考え、初期解としての主成分、および直交回転（バリマックス回転）と斜交回転（コーティミン回転）による単純パターンを求めるPythonプログラムを作成した。メインプログラムファイルMainPrg.pyとモジュールファイルinit_pca.py、PlotComp.py、rotate.py、ArrangePat.py、およびサンプルデータファイルDataSubjects.csvを圧縮ファイルpcafiles.zipにまとめた。ファイル名pcafiles.zipをクリックしてダウンロードすれば展開・解凍して利用することができる。回転のモジュールrotate.pyは独立しているので、他の回転を行う因子分析などのスクリプトで用いることができる。
メインプログラムMainPrg.py、PCAを行うクラス型InitPCAのモジュールinit_pca.py、回転を行うモジュールrotate.pyは、本ウェブサイトの後ろにリスト１、リスト２、およびリスト３として示した。
まず、主成分分析を正射影の観点（Okamoto, 2006）から説明し、その後、プログラムの使用法について説明する。
正射影の観点から見た主成分分析
データは複数の変量の値として与えられる。例えば、身長、体重、血圧の３つの変量の場合は、各人のデータは３次元ベクトル
で表される。このデータを４０人から収集すると、４０人のデータは３次元空間において４０個の点として表される。変量の数が多い場合、データの各点は変量数の次元という高次元空間において表されるが、これを低次元の空間に写してデータの分布が表されるとデータの情報が読み取りやすくなる。例えば、２次元空間でデータの分布が表された場合は、平面上にデータの分布が散布図として表される。データの分布を低次元の空間で表す方法として主成分分析がある。
直感的に説明するために、図1のようにデータが２変量からなる場合を考える。詳しい説明は、Okamoto (2006)を参照されたい。
図１
このとき、データの各点は、平面上の分布として表される。図1では、５個のデータが表されている。これを１次元空間、すなわち直線上の分布として表すことを考える。図1では、直線Line-1と直線Line-2に点の射影をとる場合を表している。射影とは、点に光を当ててその影を求めることと説明される。図1の場合は、直線に対して垂直真上から光を当てて直線上に影を求めているので正射影（orthogonal projection）と呼ぶ。影を落とす直線に対して斜め上から光を当てる場合は、射影（projection）と呼ぶ。
射影は、一般的には、s次元空間からr次元空間への影を求めることと考えられるが、正射影の場合は、r次元空間に対して垂直な方向に光の方向が決まるので、影を投影する空間が決まれば、正射影も決まる。s次元空間の点を
と表したとき、r次元空間への正射影は
と行列
による積で表すことができる（岡本、2009、pp.
65-68）。主成分分析では、元のs次元でのデータの分布を最もよく反映するr次元への正射影を求める。元の分布を最もよく反映するということを、その分散が最大になることと考える。図1の例で言えば、直線Line-1への正射影が、２次元での分布をよく反映するものであると考える。直線Line-1に比べて直線Line-2への正射影は、正射影した点の散らばりが小さいという意味で２次元での分布の様子の反映が劣ると考える。
式（1）の形式のデータが
組与えられているとする。それを行列
で表す。主成分分析では変量間の関係を見るので、各変量を平均０、標準偏差１に揃えて原点と単位の影響を除いて分析することが考えられる。式（2）の変量の値を平均０、標準偏差１の標準得点に直したもの（3）について考える。
この
に対して、射影された
の分散が最大になる正射影Pを求める。r次元空間での分散を次式（4）で表す。
なお、Deisenroch et
al.(2020)は、正射影をととの差の2乗和
を最小にするものとして求めている。
が成り立つので、分散（4）を最大にすることと、誤差の２乗和（D）を最小にすることとは同じである。
式（4）を最大にする正射影Pは次式（5）で与えられる。
ここで、
は次式の特異値分解で与えられるものである。
、
、
、
、
このとき、
である。式（6）の値は、
のとき、
となる。
相関行列をRとおくと、
となり、相関行列の固有分解により
と
を得ることができる。
射影されたの空間において、座標軸をの座標値の分散が最大になるようにお互いに直交するものを選ぶと
になる。このときの座標値は主成分（principal components）と呼ばれている。主成分を
で表すと、次式で与えられる。
空間における基底を以外に取ったときは、の座標値は成分（components）と呼ぶ。成分との関係を表す行列はパターン行列（pattern matrix）と呼ぶ。
主成分のパターン行列
は次式で与えられる。
空間における基底を以外にとると、各座標軸の解釈が容易になる場合がある。いま、基底を
にとったとする。座標間の変換を表す行列を
とおく。すなわち、
である。このを、回転（rotation）と呼ぶ。正規直交基底間の回転のときは、直交回転（orthogonal
rotation）、それ以外のときは斜交回転（oblique rotation）と呼ぶ。
基底における成分
とパターン行列
は次式で与えられる。
２つのパターン行列とが与えられると、回転は次式
で与えられる。
したがって、パターン行列とが与えられたときの成分の変換式は
となる。
以上の式は、岡本（2019）「いまさら聞けないPythonでデータ分析」の第７章「主成分分析」でも説明している。
上の主成分分析を行うPythonスクリプトを用意した。スクリプトファイルとサンプルデータファイルは、圧縮ファイルpcafiles.zipにまとめた。ダウンロード・解凍して利用することができる。
スクリプトの使用法を、表１の仮想データを用いて説明する。
表１ 仮想データDataSubjects.xlsx
|
通 番 （ID） |
国 語 （Jpn） |
英 語 （Eng） |
歴 史 （Hist） |
数 学 （Math） |
物 理 （Phys） |
化 学 （Chem） |
|
1 |
48 |
46 |
60 |
47 |
68 |
44 |
|
2 |
54 |
49 |
63 |
61 |
49 |
45 |
|
3 |
48 |
55 |
47 |
49 |
49 |
47 |
|
4 |
83 |
78 |
85 |
79 |
66 |
63 |
|
5 |
56 |
49 |
56 |
78 |
70 |
73 |
|
6 |
78 |
86 |
76 |
67 |
67 |
75 |
|
7 |
62 |
64 |
59 |
62 |
69 |
57 |
|
8 |
80 |
78 |
82 |
82 |
80 |
75 |
|
9 |
55 |
49 |
58 |
79 |
67 |
75 |
|
10 |
53 |
65 |
63 |
75 |
89 |
77 |
|
11 |
45 |
52 |
53 |
47 |
41 |
42 |
|
12 |
59 |
62 |
72 |
56 |
72 |
72 |
|
13 |
44 |
34 |
49 |
54 |
65 |
59 |
|
14 |
79 |
64 |
75 |
37 |
15 |
9 |
|
15 |
62 |
69 |
66 |
71 |
67 |
66 |
|
16 |
72 |
63 |
66 |
80 |
84 |
72 |
|
17 |
78 |
82 |
83 |
81 |
67 |
81 |
|
18 |
62 |
73 |
63 |
60 |
64 |
67 |
|
19 |
66 |
60 |
61 |
64 |
43 |
63 |
|
20 |
70 |
66 |
71 |
50 |
55 |
52 |
メインスクリプトMainPrg.pyを実行すると、入力データファイル名が求められる。入力データファイルは、Excelファイルとして用意する（図１）。1行目に変数名、１列目にデータのIDが入力されている。
図１
入力データファイル名を設定すると、ファイルが読み込まれて表示され、出力テキストファイル名の設定が求められる。出力テキストファイル名は任意である。
Input data file (*.xlsx) =
DataSubjects.xlsx
Output data file = Results.txt
入力データフィル名と出力ファイル名を設定すると計算が始まる。まず、固有分解が行われて、固有値が計算され、累積２乗和などが端末に以下のように表示される。
Lambda cum.sqr
%
Lambda-1 =
1.83432 3.36473 56.08
Lambda-2 =
1.40444 5.33718 88.95
Lambda-3 =
0.53507 5.62349 93.72
Lambda-4 =
0.46458 5.83932 97.32
Lambda-5 =
0.31091 5.93599 98.93
Lambda-6 =
0.25301 6.00000 100.00
また、スクリープロットも図２のように表示される。
図２
これらから、主成分数は２が適当であると考えられる。
図２のフォームを閉じると、次のように主成分数を聞いてくる。
Number of
components =
いまの場合、２と設定してみる。
Number of
components = 2
主成分数を設定すると計算が始まり、主成分パターンのプロットを行うときの横軸の主成分の設定が求められる。
Number of components = 2
Set the component number (from 1 to 2)
to plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim =
いまの場合、横軸の主成分として１を設定してみる。横軸の主成分を設定すると縦軸の主成分を聞いてくる。
Set the component number (from 1 to 2)
to plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
上の例では、縦軸として主成分２を設定している。縦軸の主成分を指定すると、指定した主成分の組み合わせで描かれたパターンが表示される（図３）。
図３
フォームの右上のX印をクリックしてフォームを閉じると、次の表示成分の設定が求められる。
Number of components = 2
Set the component number (from 1 to 2)
to plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
Set the component number (from 1 to 2)
to plot the pattern.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
成分として１より小さい値を設定すると（上の例では０）、初期値である主成分のパターン行列の表示は終了して、Varimax回転が行われ、Varimax回転後のパータン行列の表示における成分の指定が求められる。
Set the component number (from 1 to 2)
to plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 1
ydim = 2
上の例のように設定すると、図４の散布図が表示される。
図４
フォーム右上のX印をクリックして閉じると、次の成分の指定が求められる。０以下の値を設定すると、Varimax回転後のパターン行列の表示は終了して、斜交回転が始まる。斜交回転はクオーティミン回転である。回転における軸の選択において、足立・山本（2024）では「j < k」の制約が置かれている（p. 140）が、この制約条件では最後の軸は動かないので、制約条件を「j != k」とした（リスト３）。
斜交回転後は、成分が直交していないので、直交座標による表示はない。パターン行列は出力ファイルに出力されている。
斜交回転後、プログラムは終了し、端末画面には以下のような出力が表示される。
Set the component number (from 1 to 2)
to plot the pattern for the Varimax criterion.
If you do not want to plot the
pattern, set the number less than 1.
xdim = 0
ck =
0.06213042525789447 n_try = 1
ck = 0.6850740697132184 n_try = 2
ck = 0.47722207137964623 n_try = 3
ck = 0.0021437810652332264 n_try = 4
ck = 3.2208237454507853e-07 n_try = 5
Results.txt was saved.
「ck = 0.0621304・・・」は、クオーティミン回転における繰り返し計算における変化量である。
プログラムの実行終了後、出力ファイルResults.txtを開くと以下のようになっている。
Input data file = DataSubjects.xlsx
n_vars = 6
Var_names...
ID
Jpn
Eng
Hist
Math
Phys
Chem
1 48.000 46.000 60.000 47.000 68.000 44.000
2 54.000 49.000 63.000 61.000 49.000 45.000
3 48.000 55.000 47.000 49.000 49.000 47.000
4 83.000 78.000 85.000 79.000 66.000 63.000
5 56.000 49.000 56.000 78.000 70.000 73.000
6 78.000 86.000 76.000 67.000 67.000 75.000
7 62.000 64.000 59.000 62.000 69.000 57.000
8 80.000 78.000 82.000 82.000 80.000 75.000
9 55.000 49.000 58.000 79.000 67.000 75.000
10 53.000 65.000 63.000 75.000 89.000 77.000
11 45.000 52.000 53.000 47.000 41.000 42.000
12 59.000 62.000 72.000 56.000 72.000 72.000
13 44.000 34.000 49.000 54.000 65.000 59.000
14 79.000 64.000 75.000 37.000 15.000 9.000
15 62.000 69.000 66.000 71.000 67.000 66.000
16 72.000 63.000 66.000 80.000 84.000 72.000
17 78.000 82.000 83.000 81.000 67.000 81.000
18 62.000 73.000 63.000 60.000 64.000 67.000
19 66.000 60.000 61.000 64.000 43.000 63.000
20 70.000 66.000 71.000 50.000 55.000 52.000
Lambda cum.sqr
%
Lambda-1 1.83432 3.36473 56.08
Lambda-2 1.40444 5.33718 88.95
Lambda-3 0.53507 5.62349 93.72
Lambda-4 0.46458 5.83932 97.32
Lambda-5 0.31091 5.93599 98.93
Lambda-6 0.25301 6.00000 100.00
Pattern matrix for the principal
components =
comp.1 comp.2
Jpn -0.75184 0.61221
Eng -0.81468 0.43846
Hist -0.78026 0.54381
Math -0.80603 -0.44283
Phys -0.59683 -0.71379
Chem -0.72185 -0.63568
Principal components...
ID comp.1 comp.2
1 1.14263 -0.31126
2 0.85614 0.13922
3 1.39156 -0.16393
4 -1.43025 0.91747
5 0.08185 -1.28077
6 -1.24203 0.64068
7 0.12871 -0.19691
8 -1.66954 0.17136
9 0.04574 -1.24292
10 -0.51805 -1.35126
11 1.55819 0.17149
12 -0.18071 -0.22545
13 1.39459 -1.24212
14 1.11698 3.19106
15 -0.37087 -0.20919
16 -0.86467 -0.70070
17 -1.64668 0.40461
18 -0.16398 0.01236
19 0.25673 0.31686
20 0.11365 0.95938
Varimax rotation was applied...
Pattern for the Varimax rotated
components =
comp.1 comp.2
Jpn -0.96872 -0.04061
Eng -0.90049 -0.21231
Hist -0.94461 -0.11065
Math -0.30941 -0.86605
Phys 0.02688 -0.93004
Chem -0.11849 -0.95452
Sorted Pattern for Varimax rotation =
comp.1 comp.2
Jpn -0.96872 -0.04061
Hist -0.94461 -0.11065
Eng -0.90049 -0.21231
Chem -0.11849 -0.95452
Phys 0.02688 -0.93004
Math -0.30941 -0.86605
Components for Varimax...
ID comp.1 comp.2
1 1.06152 0.52504
2 0.54831 0.67209
3 1.15006 0.80042
4 -1.67887 -0.26219
5 0.91085 -0.90412
6 -1.35442 -0.34447
7 0.22694 -0.06197
8 -1.36301 -0.97925
9 0.85871 -0.89975
10 0.50869 -1.35481
11 1.05226 1.16195
12 0.01432 -0.28858
13 1.86754 -0.00440
14 -1.28092 3.12886
15 -0.13876 -0.40256
16 -0.18224 -1.09792
17 -1.50063 -0.78955
18 -0.13091 -0.09952
19 -0.01807 0.40741
20 -0.55136 0.79330
Quartimin rotation was applied...
Pattern for the Quartimin rotated
components =
comp.1 comp.2
Jpn -0.98386 0.07338
Eng -0.89359 -0.10999
Hist -0.95085 -0.00098
Math -0.21176 -0.84744
Phys 0.13917 -0.95262
Chem -0.00627 -0.96038
Correlation matrix of the
components...
comp.1 comp.2
comp.1 1.00000 0.23068
comp.2 0.23068 1.00000
Sorted Pattern for Promax rotation =
Jpn -0.98386 0.07338
Hist -0.95085 -0.00098
Eng -0.89359 -0.10999
Chem -0.00627 -0.96038
Phys 0.13917 -0.95262
Math -0.21176 -0.84744
Components for Quartimin...
ID comp.1 comp.2
1 1.11499 0.64553
2 0.62217 0.73158
3 1.23472 0.92937
4 -1.69791 -0.45665
5 0.80048 -0.79145
6 -1.38511 -0.50044
7 0.21828 -0.03502
8 -1.46686 -1.13187
9 0.74920 -0.79320
10 0.34902 -1.28605
11 1.17926 1.27699
12 -0.01906 -0.28493
13 1.85457 0.21394
14 -0.91147 2.95767
15 -0.18427 -0.41602
16 -0.30767 -1.11169
17 -1.58169 -0.95956
18 -0.14151 -0.11414
19 0.02904 0.40251
20 -0.45617 0.72341
足立浩平・山本倫生（2024）主成分分析と因子分析．共立出版
Deisenroth, Faisal, and Ong (2020), “Mathematics for Machine Learning.” Cambridge University Press.
岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館
岡本安晴（2019）「いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析（PyStan、PyMC）」丸善出版
山本倫生（2019）因子分析モデルにおける因子回転問題．計算機統計学、32、21−44．
#
# Principal
Components Analysis with Rotation of Components
#
Yasuharu Okamoto, 2016.10, 2025.06
#
import numpy as np
import pandas as pd
from init_pca import InitPCA
from rotate import *
from PlotComp import PlotPattern
from ArrangePat import ArrangePtn
s = input("Input data file
(*.xlsx) = ")
data = pd.read_excel(s)
#
# Prepare
the output file
#
s_out = input("Output data file =
")
f_out = open(s_out, "w")
f_out.write("Input data file =
" + s + "\n\n")
#
#
Calculate the principal components
#
calc_init_pca = InitPCA(data, f_out)
var_names, n_data, n_vars, CaseID,
n_comp, Pat, PComp = calc_init_pca.calc_pca()
"""
var_names:
A list of variable names
n_data: The
number of cases in data
n_vars: The
number of variables
CaseID: The
variable for case idenfication (string type)
n_com:
The number of componets
Pat:
The pattern matrix of the principal components
PComp:
The matrix of the principal components
"""
f_out.write("\n\nPrincipal
components...\n")
f_out.write("%10s"%var_names[0])
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PComp[i, j])
f_out.write("\n")
PlotPtn = PlotPattern( Pat, var_names,
n_vars, n_comp, " ...principal components" )
while True:
print("\nSet the
component number (from 1 to %d) to plot the pattern."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
rot = Varimax_Orthomax(Pat)
PatV = rot.rotate()
f_out.write("\n\nVarimax rotation
was applied...\n\n")
f_out.write("\nPattern for the
Varimax rotated components =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
f_out.write("%10.5f"%PatV[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names, PatV,
n_vars, n_comp)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Varimax rotation =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
PlotPtn = PlotPattern( PatV,
var_names, n_vars, n_comp, " ...Varimax rotated components" )
while True:
print("\nSet the
component number (from 1 to %d) to plot the pattern for the Varimax
criterion."%n_comp)
print("If you do not
want to plot the pattern, set the number less than 1.")
xdim = int(input("xdim
= "))
if xdim < 1:
break
ydim = int(input("ydim
= "))
PlotPtn.Plot(xdim, ydim)
PCompV = PComp @ Pat.transpose() @
PatV @ np.linalg.inv(PatV.transpose() @ PatV)
f_out.write("\n\nComponents for
Varimax...\n")
f_out.write("%10s"%var_names[0])
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PCompV[i, j])
f_out.write("\n")
obrot = Oblique_Quartimin(Pat)
Patob, Phi = obrot.rotate()
f_out.write("\n\nQuartimin
rotation was applied...\n\n")
f_out.write("\nPattern for the
Quartimin rotated components =\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_vars):
f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
f_out.write("%10.5f"%Patob[i, j])
f_out.write("\n")
f_out.write("\n\nCorrelation
matrix of the components...\n")
f_out.write("%10s"%"
")
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(i + 1)))
for j in range(0, n_comp):
f_out.write("%10.5f"%Phi[i, j])
f_out.write("\n")
ArgPtn = ArrangePtn(var_names, Patob,
n_vars, n_comp)
PtnArgd = ArgPtn.arrange()
f_out.write("\n\nSorted Pattern
for Promax rotation =\n")
for i in range(0, n_vars):
f_out.write("%10s"%PtnArgd[i][0])
for j in range(0, n_comp):
f_out.write("%10.5f"%PtnArgd[i][1][j])
f_out.write("\n")
PCompob = PComp @ Pat.T @ Patob @
np.linalg.inv(Patob.T @ Patob)
f_out.write("\n\nComponents for
Quartimin...\n")
f_out.write("%10s"%var_names[0])
for j in range(0, n_comp):
f_out.write("%10s"%("comp.%d"%(j + 1)))
f_out.write("\n")
for i in range(0, n_data):
f_out.write("%10s"%CaseID[i])
for j in range(0, n_comp):
f_out.write("%10.5f"%PCompob[i, j])
f_out.write("\n")
f_out.close()
print(s_out + " was saved.")
リスト２ 主成分分析のクラス型InitPCAのモジュールinit_pca.py
import scipy.stats as scst
import numpy as np
import matplotlib.pyplot as plt
class InitPCA:
"""
Calculation of the Principal Components
"""
def __init__(self, data,
f_out):
self.data = data
self.f_out = f_out
def calc_pca(self):
var_names = self.data.keys();
print(var_names)
n_data, n_vars = self.data.shape
n_vars -= 1
print('n_data =', n_data)
print('n_vars =', n_vars)
print('data =\n', self.data)
CaseID = self.data[var_names[0]].values
print('CaseID =\n', CaseID)
print('data.values =\n', self.data.values)
self.f_out.write("\n")
self.f_out.write("n_vars = %d\n"%n_vars)
self.f_out.write("Var_names...\n")
for
vn in var_names:
self.f_out.write("%10s"%vn)
self.f_out.write("\n")
self.f_out.write("\n")
for
i in range(0, n_data):
self.f_out.write("%10s"%CaseID[i])
for j in range(0, n_vars):
self.f_out.write("%10.3f"%self.data.values[i][j+1])
self.f_out.write("\n")
Z =
np.array(self.data.values[:,1:], dtype=float)
print('Z =\n', Z)
Z =
scst.zscore(Z, axis=0)
print('Z =\n', Z)
R =
(Z.transpose() @ Z) / n_data
print('R =\n', R)
Lmbd2, V = np.linalg.eigh(R)
print('Lmbd2 =', Lmbd2)
print('V =\n', V.T)
Lmbdsort = np.argsort(-Lmbd2)
Lmbd2 = Lmbd2[Lmbdsort]
V =
V[:,Lmbdsort]
print('Lmbd2 =', Lmbd2)
print('V =\n', V.T)
Lmbd
= []
for
i in range(0, Lmbd2.shape[0]):
if Lmbd2[i] > 1.0e-12:
Lmbd.append(Lmbd2[i]**0.5)
else:
Lmbd.append(0.0)
n_eign = len(Lmbd)
print('Lmbd =\n', Lmbd)
sum_L2 = np.sum(Lmbd2)
cum_sum = 0.0
print("\n")
print(("%30s"%"Lambda") +
("%10s"%"cum.sqr") + ("%10s"%"%") +
"\n")
self.f_out.write(
("\n\n%30s"%"Lambda")
+ ("%10s"%"cum.sqr") + ("%10s"%"%") +
"\n")
for
i in range(0, n_eign):
self.f_out.write('{0: >20s}'.format("Lambda-{0}".format(i +
1)))
self.f_out.write('{0:>10.5f}'.format(Lmbd[i]))
cum_sum += Lmbd2[i]
self.f_out.write("%10.5f"%cum_sum)
self.f_out.write("%10.2f\n"%(100.0 * cum_sum / sum_L2))
print('{0: >20s}'.format("Lambda-{0}".format(i + 1)), end =
'')
print('{0:>10.5f}'.format(Lmbd[i]), end = '')
print("%10.5f"%cum_sum, end = "")
print("%10.2f"%(100.0 * cum_sum / sum_L2))
plt.plot(np.arange(1, n_eign+1), Lmbd)
plt.ylabel('Lambda', fontsize = 16)
plt.tight_layout()
plt.show()
s0 =
input("\nNumber of components = ")
n_comp = int(s0)
if
(n_comp > n_vars):
n_comp = n_vars
elif
(n_comp < 1):
n_comp = 1
if
(n_comp > n_eign):
n_comp = n_eign
L1 =
np.full((n_comp, n_comp), 0.0)
for
i in range(n_comp):
L1[i, i] = Lmbd[i]
V1 =
V[:,:n_comp]
Pat
= V1 @ L1
self.f_out.write("\n\nPattern matrix for the principal components
=\n")
self.f_out.write("%10s"%" ")
for
i in range(0, n_comp):
self.f_out.write("%10s"%("comp.%d"%(i + 1)))
self.f_out.write("\n")
for
i in range(0, n_vars):
self.f_out.write("%10s"%var_names[i + 1])
for j in range(0, n_comp):
self.f_out.write("%10.5f"%Pat[i, j])
self.f_out.write("\n")
PComp = Z @ V1 @ np.linalg.inv(L1)
return var_names, n_data, n_vars, CaseID, n_comp, Pat, PComp
"""
Y.Okamoto,
2025.06
参考文献
足立浩平・山本倫生（２０２４）
主成分分析と因子分析、共立出版
"""
import numpy as np
import scipy.optimize as so
class Varimax_Orthomax:
def __init__(self, C):
"""
C: 負荷行列
"""
self.C = C
def rotate(self):
"""
足立・山本（2014）、p.137
"""
p =
len(self.C)
Mw =
np.identity(p) - np.ones((p,p))/p # w==1
Gmma
= self.C
while True:
Gmma2 = Gmma * Gmma
S_tilde = Gmma * (Mw @ Gmma2)
U, Delta, Vh = np.linalg.svd(self.C.T @ S_tilde)
T = U @ Vh
Gmma_new = self.C @ T
ck = np.sum((Gmma - Gmma_new)**2)
Gmma = Gmma_new
print(ck)
if ck < 1.0e-6:
break
return Gmma
# 回転後の行列
class Oblique_Quartimin():
"""
C: 負荷行列
"""
def __init__(self, C):
self.C = C
def rotate(self):
"""
足立・山本（2014）、p.140
"""
p, m
= np.shape(self.C)
Gmma
= self.C
Phi
= np.identity(m)
T =
np.identity(m)
def
calc_tjk(gmma, phi, j, k):
class Gammas():
def __init__(self, C, Phi, j, k):
#self.c = Gmma
self.C = C
self.Phi = Phi
self.j = j
self.k = k
def f_gammas(self, t):
"""
Quartimin criterion
"""
c = self.C
Phi = self.Phi
p = np.shape(c)[0]
j = self.j
k = self.k
v = 0.0
for i in range(p):
g2_ij = (1 - 2*t*Phi[j][k] + (t**2)) * (c[i][j]**2)
g2_ik = (t**2) * (c[i][j]**2) + \
2 * t * c[i][j] * c[i][k] + c[i][k]**2
v += g2_ij * g2_ik
return v
calc_g = Gammas(gmma, phi, j, k)
result = so.minimize_scalar(calc_g.f_gammas, [-10, 10])
#print(result)
return result.x
n_try = 0
while True:
Gmma0 = Gmma
for j in range(m):
for k in range(m):
if j != k:
# j < k を j != k に変更
tjk = calc_tjk(Gmma, Phi, j, k)
tjj = (1 - 2 * tjk * Phi[j][k] + (tjk**2)) ** 0.5
T_jk = np.identity(m)
T_jk[j][j] = tjj
T_jk[j][k] = tjk
Gmma = Gmma @ T_jk
inv_T_jk = np.linalg.inv(T_jk)
Phi = inv_T_jk @ Phi @ (inv_T_jk.T)
T = T @ T_jk
n_try += 1
ck = np.sum((Gmma0 - Gmma) ** 2)
print('ck =', ck, '
n_try =', n_try)
if (ck < 1.0e-6) or (n_try > 50):
break
else:
T = np.identity(m)
return Gmma, Phi # 回転後の負荷行列と共分散行列