Python/Matplotlibでグラフを描いてみよう03(主成分分析結果をグラフで可視化)

Python/Matplotlibでグラフを描いてみよう03(主成分分析結果をグラフで可視化)

データを要約(次元削減)して、多次元データの特徴を可視化してみよう




PCA(主成分分析)の簡単な使い方とグラフ表現方法

PythonでMatplotlibを使ったグラフ作成練習その3をやるよ。ナノネ、今回は主成分分析(PCA、Principal Component Analysis)を使ったデータの可視化をやってみよう。

主成分分析を使うとグラフを描くときに何かいいことあるの?

たとえばIrisデータセットだと特徴量の変数が4つ(SepalLength、SepalWidth、PetalLength、PetalWidth)あったから、一つのグラフではプロットできなかったよね。変数が多い場合は複数のグラフに分け描いたり、ちょっと頑張っても3軸の3Dグラフにして普通はもうお手上げ。そんなときでも、主成分分析を使ってうまくデータが要約ができると、変数が沢山ある多次元データでも可視化できたりするんだ。

モモノキ、PCAを利用すると4次元以上のグラフが描けるの?

そんなグラフは描けないんだけど...。データを要約することで、情報量をあまり落とさず、より少ない次元(次元削減)でデータの特徴を表現できるってことだよ。

2次元までのデータなら、普通にグラフを描けるから扱いやすいなぁ。3Dグラフは向きの調整とか、けっこう面倒そうだし。

じゃぁ早速だけど、PCAを使ったデータ処理でグラフを作ってみよう。

その前に今回使うパッケージを一式インポートしておくよ。データは、またいつものIrisで。

In [0]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
In [0]:
# ********************************
# 動作テスト環境 @2019/05
# ********************************
# Google Colaboratory
# Ubuntu 18.04.2 LTS
# Python 3.6.7
# :: matplotlib 3.0.3
# :: numpy 1.16.3
# :: scikit-learn 0.21.1
# :: pandas 0.24.2
# :: seaborn 0.9.0
In [3]:
from sklearn import datasets
iris = datasets.load_iris()

print(iris.data.shape) # データ 150x4(がく片長さ、がく片幅、花弁長さ、花弁幅)
print(iris.feature_names) # sepal length, sepal width, petal length, petal width
print(iris.target.shape) # アヤメ品種(0=setosa、1=versicolor、2=virginica)
print(iris.target_names)
(150, 4)
['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
(150,)
['setosa' 'versicolor' 'virginica']

主成分分析(PCA)を使ったデータ処理でグラフに可視化

今回はグラフの作成練習がメインなので、PCAは使い方だけ記載するね(PCAの詳しい解説はネットにたくさん情報があるから調べてみてね)。

手抜きだね...PCAについては自分で調べるよ。

PCAを実行するまえにデータの前処理をやっておくよ。PCAに使うデータは、標準化しておいた方よいケースが多いよ。特に単位が違う(重さ、長さ、kg、tonなど)データを扱うときは必須になるね。

データの標準化は、たしか...平均ゼロ、分散1にするんだっけ?平均値を引いて、標準偏差で割って...けっこう面倒かも。

普通の計算で処理でもできるけど、skeleanのStandardScalerを使うと一発だよ。

In [0]:
# データ前処理(標準化)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
iris_data_scaler = scaler.fit_transform(iris.data)

標準化したデータをチェックしておこう。データが見やすくなるように、pandasのDataFrameに格納するよ。

In [5]:
df_iris_data_scaler = pd.DataFrame(
    iris_data_scaler,
    columns=iris.feature_names)
df_iris_data_scaler.head()
Out[5]:
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm)
0 -0.900681 1.019004 -1.340227 -1.315444
1 -1.143017 -0.131979 -1.340227 -1.315444
2 -1.385353 0.328414 -1.397064 -1.315444
3 -1.506521 0.098217 -1.283389 -1.315444
4 -1.021849 1.249201 -1.340227 -1.315444
In [0]:
 
In [6]:
df_iris_data_scaler.describe()
Out[6]:
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm)
count 1.500000e+02 1.500000e+02 1.500000e+02 1.500000e+02
mean -1.690315e-15 -1.842970e-15 -1.698641e-15 -1.409243e-15
std 1.003350e+00 1.003350e+00 1.003350e+00 1.003350e+00
min -1.870024e+00 -2.433947e+00 -1.567576e+00 -1.447076e+00
25% -9.006812e-01 -5.923730e-01 -1.226552e+00 -1.183812e+00
50% -5.250608e-02 -1.319795e-01 3.364776e-01 1.325097e-01
75% 6.745011e-01 5.586108e-01 7.627583e-01 7.906707e-01
max 2.492019e+00 3.090775e+00 1.785832e+00 1.712096e+00

ほぼ平均ゼロ、標準偏差1(分散1)になるっているから、うまく標準化できたみたいだね。

続いてPCAの実行。sklearnのPCAを使うよ。PCA()でインスタンスを作って、fit_taransform(データ)を実行するだけでOK。

In [0]:
from sklearn.decomposition import PCA
pca = PCA() # PCAインスタンス作成
X_pca = pca.fit_transform(iris_data_scaler) # PCA実行して、結果をX_pcaに格納

PCAの処理結果もDataFrameに入れて、中身をチェックしてみよう。

In [0]:
# PCAの結果をDataFrameに格納にてデータチェック
df_pca = pd.DataFrame(
    X_pca,
    columns=['第1主成分(PC1)', '第2主成分(PC2)', '第3主成分(PC3)', '第4主成分(PC4)'])
In [9]:
df_pca.head()
Out[9]:
第1主成分(PC1) 第2主成分(PC2) 第3主成分(PC3) 第4主成分(PC4)
0 -2.264703 0.480027 -0.127706 -0.024168
1 -2.080961 -0.674134 -0.234609 -0.103007
2 -2.364229 -0.341908 0.044201 -0.028377
3 -2.299384 -0.597395 0.091290 0.065956
4 -2.389842 0.646835 0.015738 0.035923
In [10]:
df_pca.tail()
Out[10]:
第1主成分(PC1) 第2主成分(PC2) 第3主成分(PC3) 第4主成分(PC4)
145 1.870503 0.386966 0.256274 -0.389257
146 1.564580 -0.896687 -0.026371 -0.220192
147 1.521170 0.269069 0.180178 -0.119171
148 1.372788 1.011254 0.933395 -0.026129
149 0.960656 -0.024332 0.528249 0.163078

第1主成分、第2主成分...とかって何??

データを要約したあとの新しい合成変数で、第1主成分に最も多くの情報が詰まっているんだ。各主成分がどれくらいの情報をもっているかは、explained_variance_ratio_で確認できるよ。この値は寄与率ともいうよ。

In [11]:
pca.explained_variance_ratio_
Out[11]:
array([0.72962445, 0.22850762, 0.03668922, 0.00517871])

おおよそ第1主成分73%、第2主成分が23%、第3主成分が4%、第4主成分が0.5%って感じだね。第1主成分が一番大きくて、第2、第3と、だんだん小さくな値になるんだね。

寄与率を順に足し合わせた値を累積寄与率というよ。第〇主成分までで、元の情報をどれ位説明できるかがわかるよ。グラフでプロットしてみよう。

In [12]:
# 累積寄与率の計算
pca.explained_variance_ratio_.cumsum()
Out[12]:
array([0.72962445, 0.95813207, 0.99482129, 1.        ])
In [13]:
# 累積寄与率プロット
plt.figure()
plt.plot(np.hstack([0, pca.explained_variance_ratio_]).cumsum(), 'o-') # 0からプロット
plt.xticks(range(5))
plt.xlabel('components')
plt.ylabel('explained variance ratio')
plt.grid()

今回の場合だと、第2主成分までで元の情報の95%以上を説明できているんだね。

それじゃぁ、第1主成分と第2主成分のデータを使って、散布図を描いてみよう。

In [14]:
# PCA処理結果をプロット1
import seaborn as sns
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.scatter(X_pca[:, 0], X_pca[:, 1])
ax.set_xlabel('Principal Component 1', size=14)
ax.set_ylabel('Principal Component 2', size=14)
ax.set_title('Iirs PCA', size=16)
plt.show()

さらにアヤメの種類で、マーカーを色分けするとこんな感じ。

In [15]:
# PCA処理結果をプロット2(アヤメの種類で色分け)
import seaborn as sns
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.scatter(X_pca[0:50:, 0], X_pca[0:50, 1], label='setosa')
ax.scatter(X_pca[50:100, 0], X_pca[50:100, 1], label='versicolor')
ax.scatter(X_pca[100:150, 0], X_pca[100:150, 1], label='virginica')
ax.set_xlabel('Principal Component 1', size=14)
ax.set_ylabel('Principal Component 2', size=14)
ax.set_title('Iirs PCA', size=16)
ax.legend(loc='best', fontsize=14)
plt.show()

Irisデータの特徴を2次元散布図で可視化できたってこと?

今回の場合は、ほとんど情報を失わずに4次元あったデータを2次元で可視化できたね。PCAの用途はいくつかあるんけど、多次元データの特徴を可視化する手段として、とても有用なんだよ。

2次元散布図でデータを表現できたのは良いことなんだけど。もともとあった変数との関係がよく分かんないなぁ。名前が変わって第1主成分、第2主成分とか言われてもチンプンカンプンだよ。SepalLengthとSepalWidthとか、一体どこに行っちゃったの?

新しく作成された主成分(PC〇)と元の変数の関係は、PCAで算出された固有ベクトルの情報で確認できるよ。値の確認はPCAインスタンスのcomponents_を参照するといいよ。

In [16]:
# 固有ベクトル
pca.components_
Out[16]:
array([[ 0.52106591, -0.26934744,  0.5804131 ,  0.56485654],
       [ 0.37741762,  0.92329566,  0.02449161,  0.06694199],
       [-0.71956635,  0.24438178,  0.14212637,  0.63427274],
       [-0.26128628,  0.12350962,  0.80144925, -0.52359713]])

このままだと見にくいので、名前をつけてDataFrameに入れたものを表示するね。

In [17]:
feature_names = ['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth']
pd.DataFrame(pca.components_,
             columns=feature_names,
             index=["PC{}".format(x + 1) for x in range(len(feature_names))])
Out[17]:
SepalLength SepalWidth PetalLength PetalWidth
PC1 0.521066 -0.269347 0.580413 0.564857
PC2 0.377418 0.923296 0.024492 0.066942
PC3 -0.719566 0.244382 0.142126 0.634273
PC4 -0.261286 0.123510 0.801449 -0.523597

もっと分かりやすいように、ヒートマップで表示してみるね。

In [18]:
# ヒートマップ
import seaborn as sns
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
sns.heatmap(pca.components_,
           cmap='Blues',
           annot=True,
           annot_kws={'size': 14},
           fmt='.2f',
           xticklabels=['SepalLength', 'SepalWidth', 'PetalLength', 'PetalLength'],
           yticklabels=['PC1', 'PC2', 'PC3', 'PC4'],
           ax=ax)
plt.show()

これで第〇主成分と元の変数との関係性が分かるんだ。表示されている数値は各主成分と元の変数との相関関係(-1~1)を意味しているよ。第1主成分(PC)だとSepalLength、PetalLength、PetalWidthの3つが平均的に大きな値だね。また第2主成分はSepalWidthが突出して強く影響しているね。

なるほど。主成分がどんな内容の軸であるかは、PCAの結果値を見て判断するんだね。

うん、主成分軸の解釈やネーミングなどは分析者自身で考える必要があるんだ。今回はグラフの作成練習だから、プロットできただけでヨシとしよう。

PCAの結果の見かたは、何度か練習して慣れる必要がありそうだね。

グラフ作成練習とは少し脱線するけど、ついでだからPCAについてもう少しだけ説明するね。ざっくりだけどPCAは、データの分散共分散行列もしくは相関行列(標準化した場合)を使って、固有値と固有ベクトルを求める処理なんだ。

たしか、固有ベクトルは元の変数と新しい主成分との関係、影響の度合いを意味していたね。固有値はどんな関係があるの?

固有値は分散にあたるもので、その主成分がどれくらいの情報を持っているかが分かるよ。固有値はPCAインスタンスのexplained_variance_で確認できるよ。

In [19]:
# 固有値(分散)
print(pca.explained_variance_)
[2.93808505 0.9201649  0.14774182 0.02085386]

小数の誤差はあるけど、今回の場合、固有値(分散)を全部足すとだいたい4になるはずだよ。標準化(平均ゼロ、分散1)したデータで変数が4個あったから、分散の合計は4になるんだ。ちなみにPCA実行前後も分散の総和は変化しないからね。

もし固有値が1より大きければ(標準化データの場合)、元の変数1個よりも多くの情報を持っていることになるんだ。次元を削減したいときは、固有値の小さい(情報量の少ない)主成分は無視して、固有値が大きい主成分だけを扱うよにすれば目的が達成できるね。

各固有値を固有値の合計で割った値は、pca.explained_variance_ratio_(寄与率)の出力と同じになるよ。検算してみよう。

In [20]:
# 各固有値 / 固有値の合計 -> 寄与率
pca.explained_variance_/ pca.explained_variance_.sum()
Out[20]:
array([0.72962445, 0.22850762, 0.03668922, 0.00517871])
In [21]:
# PCAの寄与率
pca.explained_variance_ratio_
Out[21]:
array([0.72962445, 0.22850762, 0.03668922, 0.00517871])

検算結果と等しいね。

固有値は少し理解できたけど、PCAの結果としてグラフにプロットしていた数値は一体何なの?、どうやって求めたの?

グラフにプロットしていた数値は、元データを新しく作られた主成分軸に射影したときの値で、主成分得点と言われるよ。主成分得点は元のデータと対応する固有ベクトルを掛けて足し合わせることで、求めることができるよ。

たとえば、一番最初のデータを使って第1主成分の主成分得点を検算してみよう。

In [22]:
iris_data_scaler[0,] # 一番最初の元データ(標準化済み)
Out[22]:
array([-0.90068117,  1.01900435, -1.34022653, -1.3154443 ])
In [23]:
pca.components_[0,] # 第1主成分の固有ベクトル
Out[23]:
array([ 0.52106591, -0.26934744,  0.5804131 ,  0.56485654])
In [24]:
(iris_data_scaler[0,] * pca.components_[0,]).sum() # 主成分得点の計算(データ×固有ベクトルの総和)
Out[24]:
-2.264702808807589
In [25]:
X_pca[0, 0] # PCAの結果(一番最初のデータ、第1主成分の主成分得点)
Out[25]:
-2.2647028088075962

小数の誤差は少しあるけど、検算結果でも同じ値が求められたね。

最後はおまけで、R言語のbiplotっぽいグラフを描いてみるね。

In [26]:
import seaborn as sns
import matplotlib.patches as patches

# **************************
# 主成分得点プロット *******
# **************************
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
ax.scatter(X_pca[0:50:, 0], X_pca[0:50, 1], label='setosa')
ax.scatter(X_pca[50:100, 0], X_pca[50:100, 1], label='versicolor')
ax.scatter(X_pca[100:150, 0], X_pca[100:150, 1], label='virginica')
ax.set_xlabel('Principal Component 1', size=14)
ax.set_ylabel('Principal Component 2', size=14)
ax.set_title('Iirs PCA', size=16)
ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=14)

# **************************
# 固有ベクトルを図示 *******
# **************************

pc1 = pca.components_[0]
pc2 = pca.components_[1]
arrow_magnification = 2 # 矢印の倍率
feature_names = ['SepalLength', 'SepalWidth', 'PetalLength', 'PetalWidth']

# ガイド円
patch_circle = patches.Circle(
    xy=(0, 0), radius=1 * arrow_magnification,
    fc='white',
    ec='darkred',
    fill=False)

ax.add_patch(patch_circle)

# 矢印と変数ラベル
for i in range(len(feature_names)):
  # 矢印
  ax.arrow(0, 0,
           pc1[i] * arrow_magnification,
           pc2[i] * arrow_magnification,
           head_width=0.1,
           head_length=0.1,
          color='darkred')
  # 変数ラベル
  ax.text(pc1[i] * arrow_magnification * 1.2,
          pc2[i] * arrow_magnification * 1.2,
          feature_names[i])
  
plt.axis('scaled')
ax.set_aspect('equal')

plt.show()

矢印を付けると、主成分軸に対してどの変数が強く影響しているか、少し分かりやすくなるね。

今回のPCAを使ったグラフの作成練習は以上で。またね!

お疲れさまでした。またね!







スポンサーリンク

0 件のコメント :

コメントを投稿