カテゴリカルデータの主成分分析：単純２段法

カテゴリカルデータの主成分分析

単純２段法

　質問項目がカテゴリカル変量である場合は多くみられるが、その分析法として村上（1999）は非計量的主成分分析を提案している。この分析法では、k番目（）のカテゴリカル変量に対するデータを指標行列、に対する重み（数量化）行列を、の因子モデルにおける得点行列を、負荷行列をで表すとき、次式

を制約条件

のもとで最小化する解が求められる。しかし、ここではの情報を最大限反映する数量化に基づいた主成分分析を行うため、次の単純２段主成分分析を考える。この単純２段主成分分析を行うプログラムについては、このウェブサイトの後半で説明する。

2　単純２段主成分分析

各カテゴリカル変量（）の数量化を行い（第１段）、その後、数量化されたカテゴリカルデータと他の数量データがあればそれを加えて主成分分析を行う。

（１）カテゴリカルデータの数量化

数量化のための重みをの分散が最大（極大）になるように選ぶ。カテゴリ値を重みの値に置き換えたとき、データに含まれる情報をよく反映する重み付けは、そのときのデータの分散が大きくなると考える。

、　

とおき、の制約のもとでを最大にするを求めるため、目的関数を

とおく。

なので、求めるは（２）式を満たす。

（２）式に対して、例えば　　は１組の解である。

（２）式を満たす2組の解、と、およびと、に対して

が成り立つ。したがって、

である。また、

より、

となる。これより、標準化得点

を得る。

とおく。

（２）主成分分析

　すなわちの数量化による標準得点に対して

とおく。

（２a）データがカテゴリカル変量のみからなる場合

の主成分分析を行う。

（２b）データがカテゴリカル変量と数量データとからなる場合

に対して主成分分析を行う。

主成分分析については、岡本安晴（2014）、pp. 73－77で説明している。

以下のプログラム例では、（2a）カテゴリカル変量のみからなる場合について扱う。

参考文献

村上隆（1999）．科学研究費補助金（09610114）研究成果報告書

岡本安晴（2013）「カテゴリカルデータの単純２段主成分分析」日本行動計量学会第41回大会抄録集、192-195．

岡本安晴（2014）「心理学データ分析と測定」勁草書房

プログラム

まず、VC++/CLRによるプログラムを説明する。その後、Pythonスクリプト例について説明する。

VC++/CLRプログラム例

　（２a）のカテゴリカル変量のみからなる場合に対する上の方法によるプログラムPCACatVar.exeを用意した。プログラムファイル名PCACatVar.exeをクリックしてダウンロード・実行することができる。実行開始時に「***.dllがインストールされていません。」というようなエラーメッセージが表示された場合は、無料開発ツールVisual Studio Express 2012 for Windows Desktopを

http://www.microsoft.com/visualstudio/jpn/downloads

からダウンロード・インストールすると必要なファイルがインストールされ、プログラムが実行できるようになる。プログラムのソースファイルなどは圧縮ファイルPCACatVar.zipとしてまとめたので、圧縮ファイル名PCACatVar.zipをクリックしてダウンロード・解凍すれば、Visual Studio 2012で開くことができる。

プログラムPCACatVar.exeを実行すると図１のフォームが表示される。

図１

「N of Comp. = 」の右側のテキストボックスに求める主成分の個数を設定して、「GO」ボタンをクリックする。図２の入力ファイル名設定のダイアログウィンドウが表示される。

図２

入力ファイル名（ファイルの形式は後で説明する）を設定して、「開く」ボタンをクリックすると、計算結果が出力される出力用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。

図３

出力用のテキストファイル名を設定して「保存」ボタンをクリックすると計算が始まる。計算が終了すると、まず、各データ（ケース）に対する主成分の値を座標とする散布図が表示される（図４）。

図４

図４の左上は、図５のようになっている。

図５

横軸と縦軸に主成分の軸を選び直して、「Draw」ボタンをクリックすると、指定した軸による散布図が描かれる。表示されている散布図は、「Save」ボタンのクリックで保存することができる。保存した散布図は、後でWord文書などに挿入したりすることができる。「Close」ボタンをクリックすると、数量化得点と主成分との相関係数を座標とする散布図が表示される（図６）。

図６

この数量化得点と主成分との関係（構造）を表す散布図により、主成分の解釈を行うことができる。数量化得点の解釈は、出力ファイルの説明のところで扱う。図６の左上のボタンなどの配置は、図５と同じである。「Close」ボタンのクリックでフォームは閉じて、図７のダイアログボックスが表示される。

図７

「OK」ボタンのクリックでプログラムの実行は終了する。

入力ファイルと出力ファイル

　入力ファイルは、図２で設定したDataMurakami.txtの場合、以下のようになっている。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

村上　隆　平成１１年（１９９９）科研報告書、Table 1

Var1: 生国；　ドイツ・オーストリア１、フランス２、ロシア３

Var2: 作品の性格；　絶対音楽１、中間的２、表題・劇音楽３

Var3: 時期；　１９世紀前半１、１９世紀後半２、２０世紀前半３

/ カテゴリー変数の総数

/ カテゴリー変数名

Country Charactr period

/ 変数ごとのカテゴリー数

3 3 3

/ データ：　ＩＤ文字列カテゴリー値・・・

Beethoven 1 1 1

Schubert 1 2 1

Wagner 1 3 2

Schoenberg 1 1 3

Berlioz 2 3 1

Bizet 2 3 2

Debussey 2 2 3

Tchaikovsky 3 2 2

Musorgsky 3 3 2

Stravinsky 3 2 3

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

最初の「スラッシュ / 」で始まる行まではコメントである。この入力ファイルのデータが、BeethovenからStravinskyまでの１０人の作曲家についての、生国、作品の性格、時期に関するそれぞれ３カテゴリーによる評定であることが示されている。

この最初の「スラッシュ / 」で始まる行の次の行に、カテゴリー変数の個数を書く。上のデータの場合は、生国、作品の性格、時期の３変数である。

次の「スラッシュ / 」で始まる行の次の行に、各カテゴリー変数のラベル、「Country」、「Charactr」、「period」、を順番に並べて書く。これはプログラムの出力に用いられるので、分かり易いものを用いる。

３番目の「スラッシュ / 」で始まる行の次の行には、各カテゴリー変数における選択カテゴリー数を並べて書く。上のデータの場合、すべてカテゴリー数は３個である。カテゴリー値は、カテゴリー数がK個の場合は、１からKまでの整数値を用いる。カテゴリー値は名義尺度であるので、どの選択カテゴリーにどの数値を対応させるかは任意である。

４番目と５番目の「スラッシュ / 」で始まる行の間にデータ値を１行に1ケース単位で書いていく。まず、各データの表すケースのラベルを書く。これは出力に用いられるので、分かり易ものを用いる。続いて、各カテゴリー変数の選択カテゴリーを１からKまでの整数で表したものを順番に書いていく。

上のデータを入力ファイルとしたときの出力ファイルは以下のようになる。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Data file...D:\yasuharu\www\LaCoocan\booksetc\pm2010\PCACatVars\PCACatVar\PCACatVar\DataMurakami.txt

村上　隆　平成１１年（１９９９）科研報告書、Table 1

Var1: 生国；　ドイツ・オーストリア１、フランス２、ロシア３

Var2: 作品の性格；　絶対音楽１、中間的２、表題・劇音楽３

Var3: 時期；　１９世紀前半１、１９世紀後半２、２０世紀前半３

/ カテゴリー変数の総数

/ カテゴリー変数名

/ 変数ごとのカテゴリー数

CatVar[0] --> Country n_categories = 3

CatVar[1] --> Charactr n_categories = 3

CatVar[2] --> period n_categories = 3

/ データ：　ＩＤ文字列カテゴリー値・・・

Number of cases = 10

Data...

Beethoven 1 1 1

Schubert 1 2 1

Wagner 1 3 2

Schoenberg 1 1 3

Berlioz 2 3 1

Bizet 2 3 2

Debussey 2 2 3

Tchaikovsky 3 2 2

Musorgsky 3 3 2

Stravinsky 3 2 3

Indicator matrices...

Beethoven 1 0 0 1 0 0 1 0 0

Schubert 1 0 0 0 1 0 1 0 0

Wagner 1 0 0 0 0 1 0 1 0

Schoenberg 1 0 0 1 0 0 0 0 1

Berlioz 0 1 0 0 0 1 1 0 0

Bizet 0 1 0 0 0 1 0 1 0

Debussey 0 1 0 0 1 0 0 0 1

Tchaikovsky 0 0 1 0 1 0 0 1 0

Musorgsky 0 0 1 0 0 1 0 1 0

Stravinsky 0 0 1 0 1 0 0 0 1

CatVar-Country

Lambda...

3.6 3

Omega...

0.816 0.000

-0.408 0.707

-0.408 -0.707

CatVar-Charactr

Lambda...

4 2.4

Omega...

0.000 0.816

0.707 -0.408

-0.707 -0.408

CatVar-period

Lambda...

3.6 3

Omega...

-0.408 0.707

0.816 0.000

-0.408 -0.707

Variable Labels...

Country-0

Country-1

Charactr-0

Charactr-1

period-0

period-1

Z...

Country-0 Country-1 Charactr-0 Charactr-1 period-0 period-1

Beethoven 1.225 0.000 0.000 2.000 -0.816 1.291

Schubert 1.225 0.000 1.118 -0.500 -0.816 1.291

Wagner 1.225 0.000 -1.118 -0.500 1.225 0.000

Schoenberg 1.225 0.000 0.000 2.000 -0.816 -1.291

Berlioz -0.816 1.291 -1.118 -0.500 -0.816 1.291

Bizet -0.816 1.291 -1.118 -0.500 1.225 0.000

Debussey -0.816 1.291 1.118 -0.500 -0.816 -1.291

Tchaikovsky -0.816 -1.291 1.118 -0.500 1.225 0.000

Musorgsky -0.816 -1.291 -1.118 -0.500 1.225 0.000

Stravinsky -0.816 -1.291 1.118 -0.500 -0.816 -1.291

Singular values...

1.39311 1.25179 1.05663 0.89511 0.61414 0.44430

Rank of Z = 6

The principal components...

Beethoven: 1.760 -0.678

Schubert: 0.846 0.009

Wagner: -0.331 -1.006

Schoenberg: 1.538 0.453

Berlioz: -0.164 -1.266

Bizet: -1.046 -1.086

Debussey: -0.091 1.061

Tchaikovsky: -0.959 0.919

Musorgsky: -1.254 -0.277

Stravinsky: -0.299 1.870

Structure matrix...

comp. 1 comp. 2

Country-0 0.778 -0.249

Country-1 0.156 -0.491

Charactr-0 0.256 0.838

Charactr-1 0.825 -0.056

period-0 -0.733 -0.296

period-1 0.167 -0.687

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

まず、読みこまれたデータと、カテゴリー変数の値を指標行列で表したものが出力されている。

続いて、各カテゴリー変数ごとに、数量化のときに算出されたとの値が出力されている。これは、との対が、の下にの値が縦ベクトルの形で示されている。例えば、カテゴリー変数Countryの場合、に対してはであり、に対してはである。

の値のパターンから、対応する数量化の解釈ができる。Countryにおけるに対するはは、「ドイツ・オーストリア」に対して0.816を与え、「フランス」と「ロシア」には-0.408を与えるものであるので、「ドイツ・オーストリア」と（「フランス」および「ロシア」）を区別する数量化である。

次に、各カテゴリー変量を数量化した標準得点のラベルが示される。このラベルはカテゴリー変数のラベルに０からの整数値が順番に付けられた形で作成される。順番は、上のの値の大きさの順番に対応している。

これらの数量化標準得点に対して主成分分析が行われる。主成分分析を行ったときの特異値が「Singular values…」の後に示されている。

各ケースの主成分の値は、「The principal components…」の後に示されている。これらの値は、図４の散布図を描くのに用いられている。

主成分値とカテゴリー変数の標準化得点との相関係数は、「Structure matrix…」の後に示されている。この値は、図６の散布図の作成に用いられているが、主成分の解釈を行うのに有用である。例えば、図６を見ると、第１主成分とCharactr-1（絶対音楽vs＜中間、表題＞）およびCountry-0（ドイツ・オーストリアvs＜フランス、ロシア＞）の間に正の相関関係があり、このことは図４においてBeethovenとSchoenbergが右端に位置していることに対応している。図６において第１主成分はperiod-0（１９ｃ後半vs＜１９ｃ前半、２０ｃ前半＞）との間に負の相関関係があることが示されているが、これは図４ではTchaikovsky、Musorgsky、Bizet、Wagnerが左側に位置していることに対応している。第１主成分は、（ドイツ・オーストリア、絶対音楽）に対して正の値、１９ｃ後半に対して負の値をとる傾向の成分であるといえる。（生国、作品の性格）の特性と時期の特性が第１主成分の正の値と負の値によって表されている。

参考文献

岡本安晴（2010）大学生のための心理学VC++プログラミング入門．勁草書房．

Pythonスクリプト例

　（２a）のカテゴリカル変量のみからなる場合についてのPythonスクリプト例をリストP1のように用意した。スクリプトファイルなどは、まとめて圧縮ファイルPCACatVarPy.zipとした。ダウンロード・展開すれば、自由に利用できる。

リストP1のスクリプトでは、入力データはExcelファイルとして用意されたものを読み込んでいる。Excelファイルを読み込むために、ライブラリPandasを利用しているが、このときライブラリopenpyxlが必要である。コマンド

pip install openpyxl

を実行してインストールしておく必要がある。

リストP1のスクリプトを実行すると、入力データファイル名が求められる。

Excel Data File(*.xlsx) =

入力データファイルは、図P1のような形式で用意する。

図P1　村上隆　平成11年（1999）科研報告書、Table 1より

図P1のデータは、村上　隆　平成１１年（１９９９）科研報告書、Table 1より作成したものである。

Country:

１＜ドイツ・オーストリア＞、２＜フランス＞、３＜ロシア＞

Character:

１＜絶対音楽＞、２＜中間的＞、３＜表題・劇音楽＞

Period:

１＜１９世紀前半＞、２＜１９世紀後半＞、３＜２０世紀前半＞

である。

１行目に、変数のラベルを置く。

２行目に、各カテゴリ変数のカテゴリ数を設定する。１列目の列は、各データの識別名の変数であるので、カテゴリ数は書かずに空セルとする。

３行目以降に、データ値を書く。１列目はデータの識別名である。２列目以降に、各カテゴリー変数における選択カテゴリーの値を並べて書く。上のデータの場合、すべてカテゴリー数は３個である。カテゴリー値は、カテゴリー数がK個の場合は、１からKまでの整数値を用いる。カテゴリー値は名義尺度（もちろん、順序尺度でもよい）のものとして扱うので、どの選択カテゴリーにどの数値を対応させるかは任意である。

図P1のファイルDataMurakami.xlsxを入力ファイルとするときは、次のようにファイル名を設定する。

Excel Data File(*.xlsx) = DataMurakami.xlsx

入力データファイル名を設定すると、計算が始まる。

計算が進むと、まず主成分分析の特異値のグラフが表示される（図P2）

図P2

図P2のフォームを閉じると、主成分分析の主成分の個数の設定が求めらっる。

Number of components =

図P2のグラフから、主成分の個数として２を設定してみる。

Number of components = 2

主成分の個数２で計算が始まる。

まず、主成分の値を座標値とする散布図が表示される（図P3）。

図P3

図P2は、第１主成分を横座標、第２主成分を縦座標とする散布図である。

図P3のフォームを閉じると、次に描く散布図の横軸の主成分の設定が求められる。

X-Axis =

希望する主成分を設定すると、次に縦軸の主成分の設定が求められる。

主成分値の散布図の描画は終えて、次に進むときは、横軸の主成分として０以下の値を設定する。

X-Axis = 0

０以下の値を入力すると、次に構造行列の値（変数と主成分の相関係数）を座標値とする散布図が表示される（P4）。

図P4

図P4のフォームを閉じると、次の表示における主成分の設定が求められる。

X-Axis =

横軸（X-Axis）の主成分を設定すると、次に縦軸の主成分の設定が求められる。

横軸（X-Axis）の設定において、０以下の値を設定すると構造行列の描画が終了する。

構造行列の描画が終了すると、プログラムの実行終了である。

端末には、次のメッセージが表示される。

Results.txt was saved.

ファイルResults.txtを開くと、以下のようになっている。

Data file = DataMurakami.xlsx

data...

Data_ID Country Character period

0 NaN 3 3 3

1 Beethoven 1 1 1

2 Schubert 1 2 1

3 Wagner 1 3 2

4 Schoenberg 1 1 3

5 Berlioz 2 3 1

6 Bizet 2 3 2

7 Debussey 2 2 3

8 Tchaikovsky 3 2 2

9 Musorgsky 3 3 2

10 Stravinsky 3 2 3

Country

G =

[[1. 0. 0.]

[1. 0. 0.]

[0. 1. 0.]

[0. 0. 1.]

[0. 0. 1.]]

f =

[[4. 3. 3.]]

Lambda...

3.600 3.000

Omega...

0.816 0.000

-0.408 -0.707

-0.408 0.707

Country_0 Country_1

Character

G =

[[1. 0. 0.]

[0. 1. 0.]

[0. 0. 1.]

[1. 0. 0.]

[0. 0. 1.]

[0. 1. 0.]

[0. 0. 1.]

[0. 1. 0.]]

f =

[[2. 4. 4.]]

Lambda...

4.000 2.400

Omega...

-0.000 0.816

-0.707 -0.408

0.707 -0.408

Character_0 Character_1

period

G =

[[1. 0. 0.]

[1. 0. 0.]

[0. 1. 0.]

[0. 0. 1.]

[1. 0. 0.]

[0. 1. 0.]

[0. 0. 1.]

[0. 1. 0.]

[0. 0. 1.]]

f =

[[3. 4. 3.]]

Lambda...

3.600 3.000

Omega...

0.408 0.707

-0.816 0.000

0.408 -0.707

period_0 period_1

Quant. Vars...

Country_0 Country_1 Character_0 Character_1 period_0 period_1

Beethoven 1.22474 0.00000 -0.00000 2.00000 0.81650 1.29099

Schubert 1.22474 0.00000 -1.11803 -0.50000 0.81650 1.29099

Wagner 1.22474 0.00000 1.11803 -0.50000 -1.22474 0.00000

Schoenberg 1.22474 0.00000 -0.00000 2.00000 0.81650 -1.29099

Berlioz -0.81650 -1.29099 1.11803 -0.50000 0.81650 1.29099

Bizet -0.81650 -1.29099 1.11803 -0.50000 -1.22474 0.00000

Debussey -0.81650 -1.29099 -1.11803 -0.50000 0.81650 -1.29099

Tchaikovsky -0.81650 1.29099 -1.11803 -0.50000 -1.22474 0.00000

Musorgsky -0.81650 1.29099 1.11803 -0.50000 -1.22474 0.00000

Stravinsky -0.81650 1.29099 -1.11803 -0.50000 0.81650 -1.29099

Singular values...

[1.39311266 1.25179426 1.05662675 0.89510525 0.61414187 0.44430229]

Rank of genX = 6

n_comp = 2

Principal components...

comp_1 comp_2

Beethoven -1.76034 0.67795

Schubert -0.84579 -0.00930

Wagner 0.33143 1.00592

Schoenberg -1.53814 -0.45344

Berlioz 0.16427 1.26576

Bizet 1.04611 1.08558

Debussey 0.09110 -1.06104

Tchaikovsky 0.95855 -0.91854

Musorgsky 1.25392 0.27686

Stravinsky 0.29890 -1.86975

Structure matrix...

comp_1 comp_2

Country_0 -0.77829 0.24926

Country_1 0.15620 -0.49080

Character_0 0.25636 0.83772

Character_1 -0.82462 0.05613

period_0 -0.73281 -0.29594

period_1 -0.16702 0.68664

重みや主成分の値は、上のVC++/CLRのものと比べて、符号が反転しているものがあるが、これは式（２）の解、主成分の解（岡本、2014、式（5.2.2））のベクトルが符号を反転したものも解であることによる。

出力ファイルには、まず、読みこまれたデータファイル名とデータが出力されている。

続いて、変数ごとに、カテゴリー変数の値を指標行列で表したものと各カテゴリの頻度が出力されている。続いて、数量化のときに算出されたとの値が出力されている。これは、との対が、の下にの値が縦ベクトルの形で示されている。例えば、カテゴリー変数Countryの場合、λ=3.600に対してはω=(0.816,-0.408,-0.408)'であり、λ=3に対してはω=(0.000,-0.707,0.707)'である。

ωの値のパターンから、対応する数量化の解釈ができる。Countryにおけるλ=3.6に対するはω=(0.816,-0.408,-0.408)'は、「ドイツ・オーストリア」に対して0.816を与え、「フランス」と「ロシア」には-0.408を与えるものであるので、「ドイツ・オーストリア」と（「フランス」および「ロシア」）を区別する数量化である。

次に、各カテゴリー変量を数量化した標準得点のラベルが示される。このラベルはカテゴリー変数のラベルに０からの整数値が順番に付けられた形で作成される。順番は、上のλの値の大きさの順番に対応している。

これらの数量化標準得点の値が、「Quant. Vars…」の下に示されている。

これらの数量化標準得点に対して主成分分析が行われる。主成分分析を行ったときの特異値が「Singular values…」の後に示されている。

各ケースの主成分の値は、「Principal components…」の後に示されている。これらの値は、図P3の散布図を描くのに用いられている。

主成分値とカテゴリー変数の数量化標準得点との相関係数は、「Structure matrix…」の後に示されている。この値は、図P4の散布図の作成に用いられているが、主成分の解釈を行うのに有用である。例えば、図P4を見ると、第１主成分とCharactr-1（絶対音楽vs＜中間、表題＞）およびCountry-0（ドイツ・オーストリアvs＜フランス、ロシア＞）の間に負の相関関係があり、このことは図P3においてBeethovenとSchoenbergが左端に位置していることに対応している。図P4において第１主成分はperiod-0（１９ｃ後半vs＜１９ｃ前半、２０ｃ前半＞）との間に負の相関関係があることが示されているが、これは図P3ではTchaikovsky、Musorgsky、Bizet、Wagnerが右側に位置していることに対応している。第１主成分は、（ドイツ・オーストリア、絶対音楽）に対して負の値、１９ｃ後半に対して正の値をとる傾向の成分であるといえる。（生国、作品の性格）の特性と時期の特性が第１主成分の正の値と負の値によって表されている。

参考文献

岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版．

リストP1　カテゴリカル変量の主成分分析

"""

Principal Components Analysis of Categorical Data

Yasuharu Okamoto, 2023.08

"""

import pandas as pd

import numpy as np

import scipy.stats as ss

import matplotlib.pyplot as plt

fl_nm = input('Excel Data File(*.xlsx) = ')

xlfile = pd.ExcelFile(fl_nm)

df = pd.read_excel(xlfile)

fout_nm = 'Results.txt'

fout = open(fout_nm, 'w')

fout.write(f'Data file = {fl_nm}\n')

print('DataFrame =\n', df)

fout.write(f'data...\n{df}\n')

case_id = df.values.T[0][1:]

var_names = np.array(df.keys()[1:])

var_types = df.values[0][1:]

X = df.values[1:,1:]

# Quantification

genvars = []

genX = []

n = len(case_id)

for i in range(len(var_names)):

k = var_types[i]

G = np.zeros((n, k))

for j in range(n):

G[j][X[j][i] - 1] = 1

print('\n', var_names[i])

print('G =\n', G)

fout.write(f'\n{var_names[i]}')

fout.write(f'\nG =\n')

fout.write(f'{G}\n')

f = np.reshape(np.sum(G, axis=0), (1, k))

print('f =\n', f)

fout.write(f'f =\n{f}\n')

H = G - np.ones((n,1)) @ f / n

vals, vecs = np.linalg.eig(H.T @ H)

idxsort = np.argsort(-vals)

vals = vals[idxsort]

vecs = vecs.T[idxsort].T

ck_k = np.sum([vals > 1.0e-9])

print('ck_k =', ck_k)

print(f'\n{var_names[i]}:')

print('\nLambda...')

for t in range(ck_k):

print(f'{vals[t]:>20.3f}', end = '')

print()

fout.write('Lambda...\n')

for t in range(ck_k):

fout.write(f'{vals[t]:>20.3f}')

fout.write('\n')

print('Omega...\n')

fout.write('Omega...\n')

for t in range(k):

for j in range(ck_k):

print(f'{vecs[t][j]:>20.3f}', end = '')

fout.write(f'{vecs[t][j]:>20.3f}')

print()

fout.write('\n')

for t in range(ck_k):

v = f'{var_names[i]}_{t}'

genvars.append(v)

fout.write(f"{v:>20s}")

lmbd = vals[t]

omg = np.reshape(vecs[:, t], (k, 1))

z = H @ omg / ((lmbd/n)**0.5)

z = np.reshape(z, len(z))

genX.append(z)

fout.write('\n')

genX = np.array(genX).T

print('\nQuant. Vars...\n')

fout.write('\nQuant. Vars...\n')

print(f'{" ":>15s}', end = '')

fout.write(f'{" ":>15s}')

for i in range(len(genvars)):

print(f'{genvars[i]:>15s}', end = '')

fout.write(f'{genvars[i]:>15s}')

print()

fout.write('\n')

for i in range(len(genX)):

v = f'{case_id[i]:>15s}'

print(v, end = '')

fout.write(v)

for v in genX[i]:

print(f'{v:>15.5f}', end = '')

fout.write(f'{v:>15.5f}')

print()

fout.write('\n')

# Principal Components Analysis

U, S, Vh = np.linalg.svd(genX/(n**0.5))

print('\nSingular values...\n', S)

fout.write(f'\nSingular values...\n{S}\n')

plt.plot(range(len(S)), S)

plt.ylabel('Singular value')

plt.xticks(range(len(S)), np.arange(len(S))+1)

plt.show()

n_comp = int(input('Number of components = '))

rnk_S = np.sum([S > 1.0e-9])

print('\nrnk_S =', rnk_S)

print('\nRank of genX =', rnk_S)

fout.write(f'\nRank of genX = {rnk_S}\n')

if rnk_S < n_comp:

n_comp = rnk_S

print('\nn_comp =', n_comp)

fout.write(f'\nn_comp = {n_comp}\n')

C = U * (n**0.5) # Principal Components

print('\nPrincipal components...')

fout.write('\nPrincipal components...\n')

print(f'{" ":>15s}', end = '')

fout.write(f'{" ":>15s}')

for i in range(n_comp):

v = f'comp_{i+1}'

print(f'{v:>15s}', end = '')

fout.write(f'{v:>15s}')

print()

fout.write('\n')

for i in range(len(case_id)):

print(f'{case_id[i]:>15s}', end = '')

fout.write(f'{case_id[i]:>15s}')

for j in range(n_comp):

print(f'{C[i][j]:>15.5f}', end = '')

fout.write(f'{C[i][j]:>15.5f}')

print()

fout.write('\n')

max_c = np.max(C)

min_c = np.min(C)

if max_c < -min_c:

max_c = -min_c

if min_c > -max_c:

min_c = -max_c

dimx = 1

dimy = 2 if n_comp > 1 else 1

while True:

for i in range(len(case_id)):

plt.plot(C[i][dimx-1], C[i][dimy-1], 'o', c = 'b')

plt.text(C[i][dimx-1], C[i][dimy-1], case_id[i])

plt.xlabel(f'comp_{dimx}')

plt.ylabel(f'comp_{dimy}')

plt.plot([min_c, max_c], [0,0], lw = 1, c ='k', ls = ':')

plt.plot([0,0], [min_c, max_c], lw = 1, c = 'k', ls = ':')

plt.show()

dimx = int(input('X-Axis = '))

if dimx < 1:

break

if dimx > n_comp:

dimx = n_comp

dimy = int(input('Y-Axis = '))

if dimy < 1:

break

if dimy > n_comp:

dimy = n_comp

Strctr = genX.T @ C / n # Structure Matrix

dimx = 1

dimy = 2 if n_comp > 1 else 1

while True:

for i in range(len(genX[0])):

plt.plot(Strctr[i][dimx-1], Strctr[i][dimy-1], 'o', c = 'b')

plt.text(Strctr[i][dimx-1], Strctr[i][dimy-1], genvars[i])

plt.plot([-1, 1], [0,0], lw = 1, c ='k', ls = ':')

plt.plot([0,0], [-1, 1], lw = 1, c = 'k', ls = ':')

plt.xlabel(f'comp_{dimx}')

plt.ylabel(f'comp_{dimy}')

plt.show()

dimx = int(input('X-Axis = '))

if dimx < 1:

break

if dimx > n_comp:

dimx = n_comp

dimy = int(input('Y-Axis = '))

if dimy < 1:

break

if dimy > n_comp:

dimy = n_comp

print('\nStructure matrix...')

fout.write('\nStructure matrix...\n')

print(f'{" ":>15s}', end = '')

fout.write(f'{" ":>15s}')

for i in range(n_comp):

v = f'comp_{i+1}'

print(f'{v:>15s}', end = '')

fout.write(f'{v:>15s}')

print()

fout.write('\n')

for i in range(len(genvars)):

print(f'{genvars[i]:>15s}', end = '')

fout.write(f'{genvars[i]:>15s}')

for j in range(n_comp):

print(f'{Strctr[i][j]:>15.5f}', end = '')

fout.write(f'{Strctr[i][j]:15.5f}')

print()

fout.write('\n')

fout.close()

print(f'\n{fout_nm} was saved.\n')