【Python】主成分分析を感じてみたい

・お題：次元圧縮の手法で、主成分分析というやつがある。軸を引き直して情報の寄与の大きい軸だけとり出すことで、次元を下げるらしい。よく分からないので、体感してみたい。

・例えば10次元の変数セットを主成分分析で軸を取り直し、寄与の大きい軸トップ3をとって来てその後の解析に使うような用途があると思うけれど、正直訳が分からない。3次元ならイメージできるので、3次元から2次元への圧縮を図示してみたい。

・私は素人なので、主成分分析の用語の意味に関しては、大きく間違っていることがあるかもしれない。あくまで私のイメージを記載しているだけなので、正しいことは個別に調べてほしい。

・コードなどは、以下のサイトを参考にさせて頂いた。主成分分析を学びたい方はぜひ元記事をご確認いただきたい。

qiita.com

・とりあえず、3次元の図示をするためのライブラリをインポートする。

%matplotlib notebook#ぐりぐり動かせるようにしておく。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import pandas as pd

・次に、データを作成する。

a = np.linspace(0,1,50)

df = pd.DataFrame({"x":[n*0.5 + np.random.rand()/4 for n in a],
"y":[n*0.7 + np.random.rand() for n in a],
"z":[n*0.6 + np.random.rand()/10 for n in a]})

#標準化しておく
df.apply(lambda x: (x-x.mean())/x.std(), axis=0)

df.head()

・次に、3Dで図示する。

#描画エリアの作成
fig = plt.figure()
ax = fig.add_subplot(projection = '3d')

#散布図を作成
ax.scatter(df.x, df.y, df.z, s = 20, c = a)

#x,y,z軸を作成
u = [1,0,0]
v = [0,1,0]
w = [0,0,1]
ax.quiver(0, 0, 0, u, v, w, arrow_length_ratio = 0.1)

#軸のラベルを付ける
ax.text(1, 0, 0, "X-axis")
ax.text(0, 1, 0, "Y-axis")
ax.text(0, 0, 1, "Z-axis")

#表示範囲と軸を設定
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(0, 2)
ax.axis("off")

plt.show()

・少しずらしてみる。

・ひし形の鋭角の角を原点付近に持ってきたような分布になっている。

・次に、PCAしてみる。

from sklearn.decomposition import PCA

pca = PCA()
pca.fit(df)

・50個のプロットの座標を、新たに作成した軸（主成分）に打ち直す。この新たな軸における各プロットのことを、主成分得点と呼ぶらしい。

feature = pca.transform(df)#主成分空間に写像して、主成分得点を算出
df_pca = pd.DataFrame(feature, columns=["PC1", "PC2", "PC3"])
df_pca.head()

・次に、新たに作成した軸（主成分）のうち、第一主成分と第二主成分を横軸と縦軸にとって、それぞれの点を図示してみる。

plt.scatter(df_pca.PC1, df_pca.PC2, c = a, alpha=0.5)
plt.xlabel("PC1")
plt.ylabel("PC2")

・なんだかひし形っぽくなった。

・PCAは、軸を取り直しただけで、プロットの分布そのものは変化していない、、のであれば、元の3Dプロットをぐりぐり動かせば、どこかで先の二次元散布図が出てくるんじゃないか？

・むやみやたらにぐりぐりするのもしんどいので、元の三次元空間に新たに作成した軸（主成分）を図示してみる。元の三次元空間における新たに作成した軸（主成分）のことを、固有ベクトルと呼ぶらしい。固有ベクトルを呼び出してみる。

pca.components_

を実行すると、以下が返ってきた。

array([[ 0.49532361, 0.66136223, 0.56325352],
[-0.74953172, -0.00240839, 0.66196404],
[-0.43915455, 0.7500628 , -0.49451904]])

・これが何かというと、第一主成分の固有ベクトルのx,y,z成分、第二主成分の固有ベクトルのx,y,z成分、第三主成分の固有ベクトルのx,y,z成分らしい。これを軸として、最初の3Dプロットを描きなおす。。

・次に、これをぐりぐり動かして、横軸にPC1、縦軸にPC2が来るように、PC3がなるべく見えなくなるように調節してみる。

・不器用すぎてPC3がこれ以上短くなる角度にできなかったが、大体やりたかった角度にできた。先の2Dプロットと見比べてみると、確かに同じようになっている。PC1はひし形の長い方の対角線、PC2は短い方の対角線に大体対応している。

・ちなみに、一番寄与が大きいPC1をなるべく消すように、PC2とPC3で軸をとるようにぐりぐり動かしてみると…

・なるほどひし形の厚みのところがPC3になっていたらしい。これだと、ひし形であることがよく分からない。こうやって見ると、PC1とPC2でほとんどの情報を表しており、PC3は正直おまけ程度の位置づけに感じる。次に、先ほどからちょくちょく出てきている寄与の度合いを図示してみる。

import matplotlib.ticker as ticker
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
plt.plot([0] + list( np.cumsum(pca.explained_variance_ratio_)))
plt.xlabel("Number of principal components")
plt.ylabel("Cumulative contribution rate")
plt.grid()
plt.show()

・どうも情報量の90%以上がPC1とPC2で説明できるらしい。PC3がおまけ程度の寄与になっているのもなんとなく分かる。

・元の変数と、PC1とPC2の関係を知りたければ、PC1とPC2におけるx,y,zの寄与度を図示すれば良い。

plt.figure(figsize=(6, 6))
for x, y, name in zip(pca.components_[0], pca.components_[1], df.columns[1:]):
plt.text(x, y, name)
plt.scatter(pca.components_[0], pca.components_[1], alpha=0.8)
plt.grid()
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()

・ちなみに、これはx,y,z軸とPC1,PC2,PC3軸を同一空間中に図示して、PC1とPC2が横軸と縦軸になるようにぐりぐり動かした場合のx,y,z軸の位置みたいな感じ。

#描画エリアの作成
fig = plt.figure()
ax = fig.add_subplot(projection = '3d')

#x,y,z軸
u = [1,0,0]
v = [0,1,0]
w = [0,0,1]

#ax.quiver(0, 0, 0, pca.components_.T[0], pca.components_.T[1], pca.components_.T[1], arrow_length_ratio = 0.1)
ax.quiver(0, 0, 0, pca.components_[0][0],pca.components_[0][1],pca.components_[0][2], arrow_length_ratio = 0.1, color = "red")
ax.quiver(0, 0, 0, pca.components_[1][0],pca.components_[1][1],pca.components_[1][2], arrow_length_ratio = 0.1, color = "red")
ax.quiver(0, 0, 0, pca.components_[2][0],pca.components_[2][1],pca.components_[2][2], arrow_length_ratio = 0.1, color = "red")

ax.text(pca.components_[0][0],pca.components_[0][1],pca.components_[0][2], "PC1")
ax.text(pca.components_[1][0],pca.components_[1][1],pca.components_[1][2], "PC2")
ax.text(pca.components_[2][0],pca.components_[2][1],pca.components_[2][2], "PC3")

#x,y,z軸を作成
u = [1,0,0]
v = [0,1,0]
w = [0,0,1]
ax.quiver(0, 0, 0, u, v, w, arrow_length_ratio = 0.1)

#軸のラベルを付ける
ax.text(1, 0, 0, "X-axis")
ax.text(0, 1, 0, "Y-axis")
ax.text(0, 0, 1, "Z-axis")

#表示範囲
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)
ax.set_zlim(0, 2)

ax.axis("off")

#描画
plt.show()