Home

主成分分析と回転

 

主成分分析の古典的説明では、主成分は１つずつ構成され、それらの組み合わせに対する回転は不適切であるとされる。しかし、主成分分析を低次元空間への正射影として考える場合は、射影先の空間における座標系は数学的には任意である（岡本, 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、mymodule.py、およびサンプルデータファイルDataSubjects.csvを圧縮ファイルpcafiles.zipにまとめた。ファイル名pcafiles.zipをクリックしてダウンロードすれば展開・解凍して利用することができる。

なお、Jupyter Notebook用のファイルは、このウェブサイトに用意した。

メインプログラムMainPrg.pyとPCAを行うクラス型InitPCAのモジュールinit_pca.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スクリプトの実行

 

上の主成分分析を行うPythonスクリプトを用意した。スクリプトファイルとサンプルデータファイルは、圧縮ファイルpcafiles.zipにまとめた。ダウンロード・解凍して利用することができる。

スクリプトの使用法を、表１の仮想データを用いて説明する。

 

表１

通　番

（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を実行すると、入力データファイル名が求められる。

入力データファイル名を設定すると、ファイルが読み込まれて表示され、出力テキストファイル名の設定が求められる。出力テキストファイル名は任意である。

 

 

PS D:\******\pcafiles> python MainPrg.py

Input data file (*.xlsx) = 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]]

Output trxt file (*.txt) = Results.txt

 

 

入力データファイルはExcelファイルとして図２に示す形式で用意する。

 

図２　DataSubjects.xlsx

 

図２のファイルは、表1のデータに対するものである。１行目のセルに変量名を入れ、２行目からデータを入れていく。データの設定が終われば、保存する。保存を行わないと、設定されたデータはプログラムから読み込めない。

入力データフィル名と出力ファイル名を設定すると計算が始まる。まず、固有分解が行われて、固有値が計算され、累積２乗和などが端末に以下のように表示される。

 

 

                        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回転後のパターン行列の表示は終了して、斜交回転が始まる。斜交回転はpromax回転である。斜交回転後は、成分が直交していないので、直交座標による表示はない。パターン行列は出力ファイルに出力されている。

斜交回転後、プログラムは終了し、端末画面には以下のような出力が表示される。

 

 

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

Results.txt was saved.

PS D:\******\pcafiles>

 

 

プログラムの実行終了後、出力ファイル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...

                        Lambd    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.96874   0.04003

       Eng  -0.90061   0.21177

      Hist  -0.94468   0.11008

      Math  -0.30993   0.86587

      Phys   0.02632   0.93005

      Chem  -0.11907   0.95445

 

 

Sorted Pattern for Varimax rotation =

              comp.1    comp.2

       Jpn  -0.96874   0.04003

      Hist  -0.94468   0.11008

       Eng  -0.90061   0.21177

      Chem  -0.11907   0.95445

      Phys   0.02632   0.93005

      Math  -0.30993   0.86587

 

 

Components for Varimax...

        ID    comp.1    comp.2

         1   1.06183  -0.52441

         2   0.54871  -0.67176

         3   1.15054  -0.79973

         4  -1.67903   0.26119

         5   0.91030   0.90467

         6  -1.35463   0.34366

         7   0.22690   0.06211

         8  -1.36359   0.97843

         9   0.85817   0.90026

        10   0.50788   1.35511

        11   1.05295  -1.16132

        12   0.01415   0.28859

        13   1.86754   0.00552

        14  -1.27905  -3.12963

        15  -0.13901   0.40247

        16  -0.18290   1.09781

        17  -1.50110   0.78865

        18  -0.13097   0.09944

        19  -0.01783  -0.40742

        20  -0.55088  -0.79363

 

 

Promax rotation was applied...

 

 

Pattern for the Promax rotated components =

              comp.1    comp.2

       Jpn  -0.98905  -0.08941

       Eng  -0.89541   0.09616

      Hist  -0.95468  -0.01422

      Math  -0.19873   0.84757

      Phys   0.15534   0.95880

      Chem   0.00944   0.96427

 

 

Correlation matrix of the components...

              comp.1    comp.2

    comp.1   1.00000  -0.26097

    comp.2  -0.26097   1.00000

 

 

Sorted Pattern for Promax rotation =

       Jpn  -0.98905  -0.08941

      Hist  -0.95468  -0.01422

       Eng  -0.89541   0.09616

      Chem   0.00944   0.96427

      Phys   0.15534   0.95880

      Math  -0.19873   0.84757

 

 

Components for Promax...

        ID    comp.1    comp.2

         1   1.12102  -0.66116

         2   0.63142  -0.73886

         3   1.24479  -0.94585

         4  -1.69874   0.48249

         5   0.78490   0.77537

         6  -1.38782   0.52102

         7   0.21690   0.03134

         8  -1.47927   1.15133

         9   0.73378   0.77794

        10   0.32730   1.27539

        11   1.19506  -1.29122

        12  -0.02351   0.28413

        13   1.85096  -0.24327

        14  -0.86111  -2.93139

        15  -0.19017   0.41740

        16  -0.32413   1.11239

        17  -1.59093   0.98155

        18  -0.14279   0.11600

        19   0.03532  -0.40142

        20  -0.44297  -0.71319

 

 

 

参考文献

Deisenroth, Faisal, and Ong (2020), “Mathematics for Machine Learning.” Cambridge University Press.

Y. Okamoto (2006) “A Justification of Rotation in Principal Component Analysis: Projective viewpoint of PCA.” Japan Women’s University Journal (Faculty of Integrated Arts and Social Sciences), 17, 59-71.

岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館

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

岡本安晴（2019）「いまさら聞けないPythonでデータ分析：多変量解析、ベイズ統計分析（PyStan、PyMC）」丸善出版

 

 

 

リスト１　メインプログラムMainPrg.py

 

import pandas as pd

import numpy

import math

import mymodule as mm

from init_pca import InitPCA

from rotate import Varimax, Promax

from PlotComp import PlotPattern

from ArrangePat import ArrangePtn

 

"""

           Principal Components Analysis with Rotation of Components

 

                                            Yasuharu Okamoto, 2016.10, 2023.11

 

"""

 

#

#       Prepare the input data file (Excel format)

#

s = input("Input data file (*.xlsx) = ")

with pd.ExcelFile(s) as xlsx:

    data = pd.read_excel(xlsx, header = None)

print(data.values)

 

#

#       Prepare the output file

#

s_out = input("Output trxt file (*.txt) = ")

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(Pat, n_vars, n_comp)

PatV = rot.rotateV()

           

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)

 

   

 

PPP =  PatV * numpy.linalg.inv(numpy.transpose(PatV) * PatV)

PCompV = PComp * numpy.transpose(Pat) * PPP

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 = Promax(PatV, n_vars, n_comp, 4.0)

Patob, Phi = obrot.rotateP();

 

f_out.write("\n\nPromax rotation was applied...\n\n")

f_out.write("\nPattern for the Promax 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")

 

PPP =  Patob * numpy.linalg.inv(numpy.transpose(Patob) * Patob)

PCompob = PComp * numpy.transpose(Pat) * PPP

f_out.write("\n\nComponents for Promax...\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 numpy

import math

import scipy.stats as ss

import matplotlib.pyplot as plt

import mymodule as mm

 

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):  

        data = numpy.array(self.data)

        n_vars = len(data[0]) - 1

   

        #

        #       var_names: List of the variable names

        #

        var_names = data[0]

        n_data = len(data) - 1 

        CaseID = data.T[0][1:]

        print('var_names =\n', var_names)

        print('n_data =', n_data)

        print('Case_ID =\n', CaseID)

 

        X = data[1:][:,1:]

        print('X =\n', X)

 

        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"%X[i][j])

            self.f_out.write("\n")

 

        X = numpy.array(X, dtype = 'float')

        R = numpy.corrcoef(X.T)

        Lmbd2, V = numpy.linalg.eig(R)

        Lmbd = []

        for i in range(0, len(Lmbd2)):

            if Lmbd2[i] > 0.0:

                Lmbd.append(math.sqrt(Lmbd2[i]))

            else:

                Lmbd.append(0.0)

        sort_index = numpy.argsort(numpy.array(Lmbd) * (-1))

        Lmbd = numpy.array(Lmbd)[sort_index]

        n_eign = 0

        for v in Lmbd:

            if v > 0.0:

                n_eign += 1

        ListLV = V.T[sort_index].T

 

        self.f_out.write("\nLambda...\n")

        sum_L2 = 0.0

        for i in range(0, n_eign):

            sum_L2 += mm.sqr(Lmbd[i])

        cum_sum = 0.0

        self.f_out.write("%30s"%"Lambd ")

        self.f_out.write("%10s"%"cum.sqr")

        self.f_out.write("%10s\n"%"%")

        print("\n")

        print(("%30s"%"Lambda") + ("%10s"%"cum.sqr") + ("%10s"%"%") + "\n")

 

        for i in range(0, n_eign):

            self.f_out.write("%20s"%("Lambda-%d = "%(i + 1)))

            self.f_out.write("%10.5f"%Lmbd[i]) 

            cum_sum += mm.sqr(Lmbd[i]) 

            self.f_out.write("%10.5f"%cum_sum)

            self.f_out.write("%10.2f\n"%(100.0 * cum_sum / sum_L2))

            print("%20s"%("Lambda-%d = "%(i + 1)), end="")

            print("%10.5f"%Lmbd[i], end = "") 

            print("%10.5f"%cum_sum, end = "")

            print("%10.2f"%(100.0 * cum_sum / sum_L2))

 

        plt.plot(numpy.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

 

        Z = ss.zscore(X, axis = 0)

 

        L1 = numpy.diag(Lmbd[:n_comp])

 

        V1 = ListLV[:, :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 = numpy.matrix(Z @ V1 @ numpy.linalg.inv(L1))

           

        return var_names, n_data, n_vars, CaseID, n_comp, Pat, PComp

 

 

Home