アヤメ (iris) データで主成分分析をしてみよう
目次
データの準備
ここでは有名なフィッシャーのアヤメ (iris) データを利用して主成分分析を実行してみます．まず，必要なライブラリをインポートします．
import seaborn as sns
import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler # データの標準化
from pprint import pprint
import matplotlib.pyplot as plt
# from IPython.display import set_matplotlib_formats # バージョンによってはこちらを有効に
from matplotlib_inline.backend_inline import set_matplotlib_formats
set_matplotlib_formats('retina')
Seaborn に登録された iris データをロードします．このデータは萼片 (sepal) の長さと幅，花びら (petal) の長さと幅，アヤメの品種 (species) からなるデータです．
df = sns.load_dataset('iris')
print(df)
sepal_length sepal_width petal_length petal_width species 0 5.1 3.5 1.4 0.2 setosa 1 4.9 3.0 1.4 0.2 setosa 2 4.7 3.2 1.3 0.2 setosa 3 4.6 3.1 1.5 0.2 setosa 4 5.0 3.6 1.4 0.2 setosa .. ... ... ... ... ... 145 6.7 3.0 5.2 2.3 virginica 146 6.3 2.5 5.0 1.9 virginica 147 6.5 3.0 5.2 2.0 virginica 148 6.2 3.4 5.4 2.3 virginica 149 5.9 3.0 5.1 1.8 virginica [150 rows x 5 columns]
このデータを概観するには ydata-profiling のページを参照してください．このような4次元（あるいはそれ以上の高次元）のデータから特徴を発見することは難しくなります．ここでは，高次元のデータについて次元を圧縮する主成分分析を行います．主成分分析によって2次元や3次元に圧縮できれば，散布図を描くことでデータの特徴を視覚的に理解することができるようになります．
データの標準化
主成分分析を行うとき，多くの場合はあらかじめデータを標準化しておきます．これはデータから各特徴量の平均を引いて標準偏差で割るという処理です．この結果，各特徴量は平均が0，標準偏差が1の標準正規分布になります．
あらかじめ変数名のリストを作成しておきます．
variables = df.columns[:4].tolist()
# variables = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width'] # この書き方でも良い
print(variables)
['sepal_length', 'sepal_width', 'petal_length', 'petal_width']
変数のリストを利用してデータフレームから変数の列だけ取り出して Numpy 配列に変換します．そのうえで，先頭の5行だけ表示します．
mat_x = df[variables].values
print(mat_x[:5])
[[5.1 3.5 1.4 0.2]
[4.9 3. 1.4 0.2]
[4.7 3.2 1.3 0.2]
[4.6 3.1 1.5 0.2]
[5. 3.6 1.4 0.2]]
データの平均と分散を確認しておきます．
print(f"平均 : {np.mean(mat_x, axis=0)}")
print(f"分散 : {np.var(mat_x, axis=0)}")
平均 : [5.84333333 3.05733333 3.758 1.19933333] 分散 : [0.68112222 0.18871289 3.09550267 0.57713289]
小数点以下の表示桁数が多いと確認が難しくなることから，表示桁数を4桁に変更します．
np.set_printoptions(precision=4)
もう一度，平均と分散を表示して確認します．もしもここで確認した平均と分散について変数間で大きな差がないようであれば以降の標準化のステップをスキップして，そのまま主成分分析を行なってもよいでしょう．今回はあらかじめ標準化を行うことにします．
print(f"平均 : {np.mean(mat_x, axis=0)}")
print(f"分散 : {np.var(mat_x, axis=0)}")
平均 : [5.8433 3.0573 3.758 1.1993]
分散 : [0.6811 0.1887 3.0955 0.5771]
データの標準化を行なって，先頭5件のデータを表示します．ここで負の値は平均よりも小さな値であることを意味していることに注意してください．例えば，sepal_length の平均は 5.8433 ですが，もとのデータ先頭行の sepal_length は 5.1 であり，平均値よりも小さいために -0.9007 という負の値に標準化されています．一方で，sepal_width の平均値は 3.0573 ですが，もとのデータ先頭行の sepal_width は 3.5 であり，これは平均値よりも大きいことから，正の値である 1.019 という値に標準化されていることになります．さらに，分散の値も利用すればデータが平均値からどの程度離れているかも理解することが可能です．
# データの標準化
scaler = StandardScaler()
mat_x_standardized = scaler.fit_transform(mat_x)
print(mat_x_standardized[:5])
[[-0.9007 1.019 -1.3402 -1.3154]
[-1.143 -0.132 -1.3402 -1.3154]
[-1.3854 0.3284 -1.3971 -1.3154]
[-1.5065 0.0982 -1.2834 -1.3154]
[-1.0218 1.2492 -1.3402 -1.3154]]
標準化されたデータについて，平均と分散を確認します．平均が0，分散が1（つまり標準偏差も1）であることがわかります．
print("平均 :", np.mean(mat_x_standardized, axis=0))
print("分散 :", np.var(mat_x_standardized, axis=0))
平均 : [-4.7370e-16 -7.8160e-16 -4.2633e-16 -4.7370e-16] 分散 : [1. 1. 1. 1.]
上の平均については自動的に指数表示となっていることから，指数表示を行わないように設定を変更します．
np.set_printoptions(suppress=True)
改めて，平均と分散を確認すると結果が読みやすくなりました．
print("平均 :", np.mean(mat_x_standardized, axis=0))
print("分散 :", np.var(mat_x_standardized, axis=0))
平均 : [-0. -0. -0. -0.] 分散 : [1. 1. 1. 1.]
主成分分析
上で標準化したデータを用いて主成分分析を行います．今回は主成分の数を変数の数と同じ 4 に設定します．
pca = PCA(n_components = 4)
pca.fit(mat_x_standardized)
PCA
実行結果にどのようなものが含まれているかを確認するために，得られた結果をすべて表示します．このとき，表示を見やすくするために pprint を利用してます．
pprint(vars(pca))
{'_fit_svd_solver': 'full',
'components_': array([[ 0.5211, -0.2693, 0.5804, 0.5649],
[ 0.3774, 0.9233, 0.0245, 0.0669],
[-0.7196, 0.2444, 0.1421, 0.6343],
[-0.2613, 0.1235, 0.8014, -0.5236]]),
'copy': True,
'explained_variance_': array([2.9381, 0.9202, 0.1477, 0.0209]),
'explained_variance_ratio_': array([0.7296, 0.2285, 0.0367, 0.0052]),
'iterated_power': 'auto',
'mean_': array([-0., -0., -0., -0.]),
'n_components': 4,
'n_components_': 4,
'n_features_in_': 4,
'n_oversamples': 10,
'n_samples_': 150,
'noise_variance_': 0.0,
'power_iteration_normalizer': 'auto',
'random_state': None,
'singular_values_': array([20.9231, 11.7092, 4.6919, 1.7627]),
'svd_solver': 'auto',
'tol': 0.0,
'whiten': False}
実行結果から主成分ベクトルだけを表示します．
print("主成分")
print(pca.components_)
主成分 [[ 0.5211 -0.2693 0.5804 0.5649] [ 0.3774 0.9233 0.0245 0.0669] [-0.7196 0.2444 0.1421 0.6343] [-0.2613 0.1235 0.8014 -0.5236]]
標準化済みのデータについて平均を表示します．これはもちろん0になるはずです．
print("平均")
print(pca.mean_)
平均 [-0. -0. -0. -0.]
主成分の分散，つまりばらつき度合いを取得します．これを見ると，第1主成分のばらつき度合いが最も大きく，第2，第3，第4と徐々に小さくなっています．
print("分散")
print(pca.explained_variance_)
分散 [2.9381 0.9202 0.1477 0.0209]
それぞれの主成分の分散が分散の総和に占める割合を寄与率と呼びます．つまり，次のようにして寄与率を求めることができます．
print(np.sum(pca.explained_variance_).round(4)) # 分散の総和を求めて表示する
print(pca.explained_variance_ / np.sum(pca.explained_variance_)) # 分散を総和で割って表示する
4.0268 [0.7296 0.2285 0.0367 0.0052]
しかしながら，pca には各主成分の寄与率がすでに計算されているので，わざわざ計算することなく，それを表示するとよいでしょう．同じ結果が得られるはずです．
print("寄与率")
print(pca.explained_variance_ratio_)
寄与率 [0.7296 0.2285 0.0367 0.0052]
np.cumsum() を利用すると累積寄与率を求めることができます．この結果，第1主成分だけで 72.96%，第2主成分まで含めると 95.81% を説明可能であることがわかったので，以降では第2主成分まで利用して分析をすることにします．
print("累積寄与率")
print(np.cumsum(pca.explained_variance_ratio_))
累積寄与率 [0.7296 0.9581 0.9948 1. ]
いま求めた寄与率と累積寄与率をグラフにして確認します．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(np.arange(4)+1, pca.explained_variance_ratio_)
ax.scatter(np.arange(4)+1, pca.explained_variance_ratio_)
ax.set_ylim(-0.1,1.1)
ax.set_xticks(np.arange(1, 5, 1))
ax.set_xlabel('principle components')
ax.set_ylabel('contribution rate')
# plt.savefig('pca_contribution.png', dpi=300, facecolor='white')
plt.show()
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(np.arange(4)+1, np.cumsum(pca.explained_variance_ratio_))
ax.scatter(np.arange(4)+1, np.cumsum(pca.explained_variance_ratio_))
ax.set_ylim(-0.1,1.1)
ax.set_xticks(np.arange(1, 5, 1))
ax.set_xlabel('principle components')
ax.set_ylabel('cumurative contribution rate')
# plt.savefig('pca_cumurative_contribution.png', dpi=300, facecolor='white')
plt.show()
因子負荷量
因子負荷量は主成分にどの変数が強く寄与しているかを求めることができる量です．この因子負荷量は主成分の分散（標準偏差）と主成分ベクトルから求めることができます．まず，主成分の分散を改めて表示します．
pca.explained_variance_
array([2.9381, 0.9202, 0.1477, 0.0209])
分散について平方根を計算すると標準偏差が得られます．
np.sqrt(pca.explained_variance_)
array([1.7141, 0.9593, 0.3844, 0.1444])
次に，主成分ベクトルを表示します．
print(pca.components_)
[[ 0.5211 -0.2693 0.5804 0.5649] [ 0.3774 0.9233 0.0245 0.0669] [-0.7196 0.2444 0.1421 0.6343] [-0.2613 0.1235 0.8014 -0.5236]]
第1主成分に関する主成分ベクトルとそれに関する標準偏差を取り出します．
print(pca.components_[0,])
print(np.sqrt(pca.explained_variance_[0]).round(4))
[ 0.5211 -0.2693 0.5804 0.5649] 1.7141
第1主成分の因子負荷量を求めます．これは主成分ベクトルと標準偏差の積であることに注意してください．
fl1 = pca.components_[0,] * np.sqrt(pca.explained_variance_[0])
print(fl1)
[ 0.8932 -0.4617 0.9949 0.9682]
同じように第2主成分の主成分ベクトルと標準偏差を取り出します．
print(pca.components_[1,])
print(np.sqrt(pca.explained_variance_[1]).round(4))
[0.3774 0.9233 0.0245 0.0669] 0.9593
それらの積を計算することで第2主成分の因子負荷量を求めることができます．
fl2 = pca.components_[1,] * np.sqrt(pca.explained_variance_[1])
print(fl2)
[0.362 0.8857 0.0235 0.0642]
ここで，得られた因子負荷量について考察を行います．まず，第1主成分の因子負荷量について散布図を描きます．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(np.arange(4),fl1, alpha=0.5)
ax.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax.set_ylim(-1.1,1.1)
ax.set_title('PC1')
ax.set_xlabel('variables')
ax.set_ylabel('factor loadings')
for i, label in enumerate(variables):
plt.annotate(label, (i, fl1[i]))
# plt.savefig('pca_factor_loadings1.png', dpi=300, facecolor='white')
plt.show()
上の図は変数ごとの因子負荷量をプロットしています．データを標準化して分析していることから，この因子負荷量は -1.0 から +1.0 の範囲の値になります．この因子負荷量が0に近いとその変数が主成分に与える影響は小さいことを意味しています．一方で，+1.0 や -1.0 に近いほどその変数が主成分に与える影響は大きいことを意味してます．したがって，上の図から，第1主成分に影響が強い変数は sepal_length，petal_length，petal_width であり，sepal_width が第1主成分に与える影響は小さいことがわかります．つまり，第1主成分が大きいと計算されたデータは花びら (petal) が長く，幅広く，萼片 (sepal) も大きいものと予想されます．
同様に第2主成分についても因子負荷量をプロットします．
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(np.arange(4), fl2, alpha=0.5)
ax.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax.set_ylim(-1,1)
ax.set_xticks(np.arange(0, 4, 1))
ax.set_title('PC2')
ax.set_xlabel('variables')
ax.set_ylabel('factor loadings')
for i, label in enumerate(variables):
plt.annotate(label, (i, fl2[i]))
# plt.savefig('pca_factor_loadings2.png', dpi=300, facecolor='white')
plt.show()
上の図から，第2主成分には sepal_width が強く影響を与えていることがわかります．つまり，第2主成分が大きい個体は萼片 (sepal) の幅が大きいと判断されるはずです．
第1主成分の因子負荷量を横軸に，第2主成分の因子負荷量を縦軸にして因子負荷量をプロットしてみます．このプロットは後のステップでも利用されることに注意してください．
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.scatter(fl1, fl2, alpha=0.5)
ax.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax.axvline(0, color='gray', linestyle='--', linewidth=0.5)
ax.set_xlim(-1,1)
ax.set_ylim(-.1,1)
ax.set_xlabel('factor loadings for pc1')
ax.set_ylabel('factor loadings for pc2')
for i, label in enumerate(variables):
plt.annotate(label, (fl1[i], fl2[i]))
# plt.savefig('pca_factor_loadings.png', dpi=300, facecolor='white')
plt.show()
散布図を描く
次は主成分分析によって2次元に圧縮したデータを2次元平面上にプロットします．このためには pca.transform を使って標準化されたデータを主成分の座標に変換します．
transformed = pca.transform(mat_x_standardized)
print("先頭の5件")
print(transformed[:5])
先頭の5件 [[-2.2647 0.48 -0.1277 -0.0242] [-2.081 -0.6741 -0.2346 -0.103 ] [-2.3642 -0.3419 0.0442 -0.0284] [-2.2994 -0.5974 0.0913 0.066 ] [-2.3898 0.6468 0.0157 0.0359]]
品種のデータも含めて散布図を描きやすいように，データフレームを作成します．
df_result = pd.DataFrame({
'pc1': transformed[:,0], # 第1主成分
'pc2': transformed[:,1], # 第2主成分
'species': df['species']
})
print(df_result)
pc1 pc2 species 0 -2.264703 0.480027 setosa 1 -2.080961 -0.674134 setosa 2 -2.364229 -0.341908 setosa 3 -2.299384 -0.597395 setosa 4 -2.389842 0.646835 setosa .. ... ... ... 145 1.870503 0.386966 virginica 146 1.564580 -0.896687 virginica 147 1.521170 0.269069 virginica 148 1.372788 1.011254 virginica 149 0.960656 -0.024332 virginica [150 rows x 3 columns]
品種の列から集合 (set) で重複を削除してリスト化し，さらにソートします．
species = sorted(list(set(df['species'].values)))
# species = ['setosa', 'versicolor', 'virginica'] # この書き方でも良いですが，数が多いと大変になります
species
['setosa', 'versicolor', 'virginica']
2次元に圧縮したデータについて散布図を描きます．
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
colors = ['Red', 'Blue', 'Pink']
markers = ['o', 'x', 'v']
for i, s in enumerate(species):
x = df_result.loc[df['species'] == s, 'pc1']
y = df_result.loc[df['species'] == s, 'pc2']
ax.scatter(x, y, alpha=0.5, label=f"{s}", color=colors[i], marker=markers[i])
ax.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax.axvline(0, color='gray', linestyle='--', linewidth=0.5)
ax.set_title("PCA")
ax.set_xlabel('pc1')
ax.set_ylabel('pc2')
ax.legend(loc='lower right')
# plt.savefig('pca.png', dpi=300, facecolor='white')
plt.show()
上の図では，setosa がグラフの左側にプロットされていますことが読み取れます．横軸の第1主成分は右方向に行くほど花びらが大きいことを意味していることから，setosa は他の品種と比べて花びらが小さい特徴を持つことがわかります．一方で，右側にプロットされた virginica は花びらが大きいという特徴がありそうです．しかしながら，versicolor との境界ははっきりせず一部は交差していることから，花びらの大きさだけで versicolor と virginica を完全に識別することできなさそうです．
また，縦軸の第2主成分は主に萼片の幅を意味していました．setosa は萼片の幅が広いものから，狭いものまで個体差があるようです．versicolor について萼片の幅は狭いものから中程度まで，virginica については幅の広い個体から中程度の個体まであることがわかります．
さらに，biplot を描き，この biplot から各個体の特徴が視覚的に把握できることを確認します．上で作成した各個体の散布図について，プロットのx軸，y軸の範囲がいずれも2になるようにあらかじめ値を調整します．そのために必要となる xscale，yscale を次の通り計算します．
xscale = 2.0 / (df_result['pc1'].max() - df_result['pc1'].min())
yscale = 2.0 / (df_result['pc2'].max() - df_result['pc2'].min())
print(xscale)
print(yscale)
0.3286748987214411 0.3745067766597033
いま求めた xscale，yscale を乗じて個体の第1，第2主成分をプロットします．さらに同じグラフ上に因子負荷量のベクトルも矢印でプロットします．なお，この矢印はここで作成した因子負荷量のプロットと一致していることに注意してください．
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
colors = ['Red', 'Blue', 'Pink']
markers = ['o', 'x', 'v']
# 軸の幅が 2.0 になるように調整
xscale = 2.0 / (df_result['pc1'].max() - df_result['pc1'].min())
yscale = 2.0 / (df_result['pc2'].max() - df_result['pc2'].min())
for i, s in enumerate(species):
x = df_result.loc[df['species'] == s, 'pc1'] * xscale
y = df_result.loc[df['species'] == s, 'pc2'] * yscale
ax.scatter(x, y, alpha=0.2, label=f"{s}", color=colors[i], marker=markers[i])
for i in range(len(variables)):
ax.arrow(
x=0, y=0, # 始点の座標
dx=fl1[i], dy=fl2[i], # 終点の座標
width=0.001,
head_width=0.03,
head_length=0.05,
length_includes_head=True,
color='gray'
)
for i, label in enumerate(variables):
plt.annotate(label, (fl1[i], fl2[i]))
for i in range(len(transformed[:,0])):
plt.annotate(i, (transformed[:,0][i]*xscale, transformed[:,1][i]*yscale))
ax.axhline(0, color='gray', linestyle='--', linewidth=0.5)
ax.axvline(0, color='gray', linestyle='--', linewidth=0.5)
ax.set_title("PCA")
ax.set_xlabel('pc1')
ax.set_ylabel('pc2')
ax.legend(loc='lower right')
# plt.savefig('pca_arrow.png', dpi=300, facecolor='white')
plt.show()
上の図について若干の考察を加えてみましょう．図に配置された矢印の向きを確認すると，sepal_width の矢印は11時の方向を向いています．この方向に進むと萼片の幅が大きくなることを意味しています．また他の3本の矢印は2時から3時の方向を向いています．変数の順序を確認しやすいようにもう一度表示しておきます．
print(variables)
['sepal_length', 'sepal_width', 'petal_length', 'petal_width']
図の左上にある 15 のデータは11時の方向に大きな値になっている（第2主成分の値が大きい）ことから，萼片の幅 (sepal_width) が大きいはずです．そのことを確認します．また第1主成分が小さいことから他の変数は小さいはずです．
print(mat_x[15])
[5.7 4.4 1.5 0.4]
図の上下中央，右端にプロットされている 118 のデータは第1主成分の値が大きいことから，花びらの長さと幅 (petal_length, petal_width)，さらに萼片の長さ (sepal_length) も大きいはずです．そのことを確認します．15 のデータと比較するとその事がわかるはずです．
print(mat_x[118])
[7.7 2.6 6.9 2.3]
図の右下（4~5時方向）にある 119 のデータは萼片の幅の矢印と反対方向にあることからその値が小さいはずです．
print(mat_x[119])
[6. 2.2 5. 1.5]
図の左下にプロットされた 41 のデータはすべての値が小さいはずです．
print(mat_x[41])
[4.5 2.3 1.3 0.3]
図の左上にある 18, 14, 33 のデータは11時方向の矢印と概ね平行に並んでいます．つまり，この順番で sepal_width が大きくなっていると予想できます．このことを確認します．
print(mat_x[18])
print(mat_x[14])
print(mat_x[33])
[5.7 3.8 1.7 0.3] [5.8 4. 1.2 0.2] [5.5 4.2 1.4 0.2]
pca 処理の確認
上では pca() 関数を呼び出すだけで簡単に主成分分析を行うことができました．この処理が具体的にどのようなものであるかをここで確認しておきます．主成分分析の中身は特異値分解 (singular value decomposition) そのものです．
まず，上の主成分分析で得られた pca インスタンスから主成分ベクトルの値（行列）である pca.components_ 属性を表示します．
print(pca.components_)
[[ 0.5211 -0.2693 0.5804 0.5649] [ 0.3774 0.9233 0.0245 0.0669] [-0.7196 0.2444 0.1421 0.6343] [-0.2613 0.1235 0.8014 -0.5236]]
続いて，特異値である pca.singular_values_ も表示します．
print(pca.singular_values_)
[20.9231 11.7092 4.6919 1.7627]
特異値分解は scipy ライブラリを利用して実行できます（Numpyでも可能です）．ライブラリをインポートして，特異値分解を実行します．
import scipy
P, S, Qt = scipy.linalg.svd(mat_x_standardized) # 特異値分解
特異値分解の結果，変数 S に特異値が格納されています．
print(S)
[20.9231 11.7092 4.6919 1.7627]
変数 Qt には主成分ベクトルが格納されています．これらの値は符号が異なる部分があるものの，いずれも pca() を実行して得られた結果と同等であることが確認できました．
print(Qt)
[[ 0.5211 -0.2693 0.5804 0.5649] [-0.3774 -0.9233 -0.0245 -0.0669] [ 0.7196 -0.2444 -0.1421 -0.6343] [ 0.2613 -0.1235 -0.8014 0.5236]]
さらに，上の主成分分析では主成分の分散 pca.explained_variance_ や寄与率 pca.explained_variance_ratio_ も得られていました．これらを改めて確認します．
print(pca.explained_variance_)
print(pca.explained_variance_ratio_)
[2.9381 0.9202 0.1477 0.0209] [0.7296 0.2285 0.0367 0.0052]
主成分の分散は特異値 S の分散を意味します．データ数 m も利用して，次のように求めることができ，上の結果を一致していることが確認できました．
m = mat_x_standardized.shape[0] # データ数
explained_variance = (S ** 2) / (m - 1)
print(explained_variance)
[2.9381 0.9202 0.1477 0.0209]
一方で，寄与率は各主成分がどれだけ変動の割合を占めるかということを意味するので，分散の総和で割ることで求めることができます．これも上の結果と一致しているはずです．
total_variance = np.sum(explained_variance)
explained_variance_ratio = explained_variance / total_variance
print(explained_variance_ratio)
[0.7296 0.2285 0.0367 0.0052]
transform 処理の確認
上では，標準化されたデータを主成分に変換するために pca.transform を用いました．より理解を深めるために，この計算がどのように行われているのかを確認します．
まず，標準化されたデータを表示して確認します．
print("先頭の5件")
print(mat_x_standardized[:5])
先頭の5件
[[-0.9007 1.019 -1.3402 -1.3154]
[-1.143 -0.132 -1.3402 -1.3154]
[-1.3854 0.3284 -1.3971 -1.3154]
[-1.5065 0.0982 -1.2834 -1.3154]
[-1.0218 1.2492 -1.3402 -1.3154]]
主成分ベクトルの値（行列）も表示します．
print(pca.components_.T)
[[ 0.5211 0.3774 -0.7196 -0.2613] [-0.2693 0.9233 0.2444 0.1235] [ 0.5804 0.0245 0.1421 0.8014] [ 0.5649 0.0669 0.6343 -0.5236]]
先頭データの第1主成分はデータの先頭行と主成分ベクトル先頭列の値を用いて次のように計算することができます．
print(f"{(-0.9007)*(0.5211) + (1.019)*(-0.2693) + (-1.3402)*(0.5804) + (-1.3154)*(0.5649):.4f}")
-2.2647
先頭データの第2主成分も同様に，データの先頭行と主成分ベクトルの第2列を用いて計算できます．（なお，四捨五入して表示された値を用いて計算していることから若干の誤差は生じることに注意してください．）
print(f"{(-0.9007)*(0.3774) + (1.019)*(0.9233) + (-1.3402)*(0.0245) + (-1.3154)*(0.0669):.4f}")
0.4801
上のような手順で各データに対する主成分を求めることができますが，次の方法によって，すべてのデータと主成分の組み合わせに対してまとめて求めることができます．なお，@ は行列の積を求める演算子です．ハイライトした値が上で求めた値と一致していること（上で求めた方法では僅かな誤差があることも）を確認してください．
mat_tf = mat_x_standardized @ pca.components_.T
# mat_tf = np.dot(mat_x_standardized, pca.components_.T) # この書き方でも良い
print("先頭の5件")
print(mat_tf[:5])
先頭の5件
[[-2.2647 0.48 -0.1277 -0.0242]
[-2.081 -0.6741 -0.2346 -0.103 ]
[-2.3642 -0.3419 0.0442 -0.0284]
[-2.2994 -0.5974 0.0913 0.066 ]
[-2.3898 0.6468 0.0157 0.0359]]
このようにして求められた結果が，transform によって求めたものと一致していることを確認します．
print("先頭の5件")
print(transformed[:5])
先頭の5件 [[-2.2647 0.48 -0.1277 -0.0242] [-2.081 -0.6741 -0.2346 -0.103 ] [-2.3642 -0.3419 0.0442 -0.0284] [-2.2994 -0.5974 0.0913 0.066 ] [-2.3898 0.6468 0.0157 0.0359]]
データの標準化を行わない場合
ここまでの手順では，データを平均0，分散1に標準化してから主成分分析を行いました．データの平均や標準偏差について変数間で大きな差が認められない場合には，もとのデータからそのまま主成分分析を行うことも可能です．ここでは，実際にデータを標準化することなく主成分分析を行なった結果だけを示しておきます．
データを読み込んだあと，標準化せずにそのまま pca.fit 関数に与えます．
pca = PCA(n_components = 4)
pca.fit(mat_x) # データの行列をそのまま与える
pprint(vars(pca))
{'_fit_svd_solver': 'full',
'components_': array([[ 0.3614, -0.0845, 0.8567, 0.3583],
[ 0.6566, 0.7302, -0.1734, -0.0755],
[-0.582 , 0.5979, 0.0762, 0.5458],
[-0.3155, 0.3197, 0.4798, -0.7537]]),
'copy': True,
'explained_variance_': array([4.2282, 0.2427, 0.0782, 0.0238]),
'explained_variance_ratio_': array([0.9246, 0.0531, 0.0171, 0.0052]),
'iterated_power': 'auto',
'mean_': array([5.8433, 3.0573, 3.758 , 1.1993]),
'n_components': 4,
'n_components_': 4,
'n_features_in_': 4,
'n_oversamples': 10,
'n_samples_': 150,
'noise_variance_': 0.0,
'power_iteration_normalizer': 'auto',
'random_state': None,
'singular_values_': array([25.1 , 6.0131, 3.4137, 1.8845]),
'svd_solver': 'auto',
'tol': 0.0,
'whiten': False}
また，pca.transform の処理と同じ処理を行うには，データから列の平均を引いて標準化をしなければなりません．その後主成分ベクトルとの積を求めるとよいでしょう．
smat_x = mat_x - pca.mean_
mat_tf = smat_x @ pca.components_.T
第1主成分と第2主成分の散布図を示します．
標準化しない場合は矢印の長さ（ベクトルの大きさ）が変数ごとに大きく異なる結果が得られました．
