日: 2019年6月25日

主成分分析

物理情報数学Bで説明していている対称行列が出てくる例です．固有ベクトルへの射影も関係します．

#
# 主成分分析
#
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

N = 100
d = 16
x = 50 + np.random.randn(N) * d
y = 0.8 * x + np.random.randn(N) * d * 0.5 + 10
plt.figure(figsize=(5, 5))
plt.xlim([0, 100])
plt.ylim([0, 100])
plt.xlabel('x')
plt.ylabel('y')
plt.scatter(x, y, color='blue', alpha=0.5)

features = np.vstack([x, y]).T
features2 = features.copy()
ave = np.average(features,0)
print("平均")
print(ave)
features = features - ave
cv = np.dot(features.T, features) / (N - 1)
print("共分散行列")
print(cv)
W, v = np.linalg.eig(cv)
print("共分散行列の固有ベクトルを列にもつ行列")
print("{0}".format(v))
print("固有ベクトル")
print("{0}".format(v[:,0].reshape(2,1)))
print("固有ベクトル")
print("{0}".format(v[:,1].reshape(2,1)))
print("共分散行列の固有値：{0}".format(W))
# 寄与率は固有値／（固有値の総和）
p = np.trace(cv)
print("寄与率1: {0}, 2: {1}".format(W[0]/p, W[1]/p))
# 固有値の大きい順にソート
ind = np.argsort(W)[::-1]
sortedW = np.sort(W)[::-1]
sortedv = v[:,ind]
# 固有ベクトル方向の直線
plt.plot(np.array([0, 100]), np.array([sortedv[1,0]/sortedv[0,0]*(0 - ave[0])+ave[1], sortedv[1,0]/sortedv[0,0]*(100 - ave[0])+ave[1]]),  color='blue', alpha=0.5)
plt.plot(np.array([0, 100]), np.array([sortedv[1,1]/sortedv[0,1]*(0 - ave[0])+ave[1], sortedv[1,1]/sortedv[0,1]*(100 - ave[0])+ave[1]]),  color='blue', alpha=0.5)
#
t1 = np.zeros([2, N])
# 固有ベクトル方向に射影する
#p1 = v[:,0].reshape(2,1)
#p2 = v[:,1].reshape(2,1)
#for i in range(N):
#  t1[0, i] = np.dot(p1.T, (features.T[:,i]).reshape(2,1)) # 射影するベクトルは単位ベクトルなので2乗ノルムの2乗で割る必要なし
#  t1[1, i] = np.dot(p2.T, (features.T[:,i]).reshape(2,1))
t1 = np.dot(sortedv.T, features.T)
#
plt.figure(figsize=(5, 5))
plt.xlim([-70, 70])
plt.ylim([-70, 70])
plt.xlabel('p1')
plt.ylabel('p2')
plt.scatter(t1[0,:], t1[1,:], color='red', alpha=0.5)
# 3σの楕円
theta = np.arange(0, 361, 3)
X = np.cos(theta/180*np.pi) * np.sqrt(sortedW[0]) * 3
Y = np.sin(theta/180*np.pi) * np.sqrt(sortedW[1]) * 3
plt.plot(X, Y, color='red', alpha=0.8)
plt.show()
#
#
print("主成分分析のライブラリを使う")
pca = PCA(n_components=2)
pca.fit(features2)
print("components")
print("{0}".format(pca.components_))
print("mean={0}".format(pca.mean_))
print("covariance")
print("{0}".format(pca.get_covariance()))
t2 = pca.fit_transform(features)
print("contribution")
print(pca.explained_variance_ratio_)
plt.figure(figsize=(5, 5))
plt.xlim([-70, 70])
plt.ylim([-70, 70])
plt.xlabel('p1')
plt.ylabel('p2')
plt.scatter(t2[:,0], t2[:,1], color='green', alpha=0.5)
print("固有ベクトルの符号が反転していると左右・上下に反転がある")
plt.show()
平均
[52.53235793 51.86342693]
共分散行列
[[244.04297958 193.42743647]
 [193.42743647 216.02097448]]
共分散行列の固有ベクトルを列にもつ行列
[[ 0.73220426 -0.6810851 ]
 [ 0.6810851   0.73220426]]
固有ベクトル
[[0.73220426]
 [0.6810851 ]]
固有ベクトル
[[-0.6810851 ]
 [ 0.73220426]]
共分散行列の固有値：[423.96619622  36.09775785]
寄与率1: 0.9215375220555805, 2: 0.07846247794441952
直線は固有ベクトルの方向を表す
第1主成分と第2主成分を軸にとったもの．3σの長半径と短半径を持つ楕円を重ねて描いた
主成分分析のライブラリを使う
components
[[ 0.73220426  0.6810851 ]
 [-0.6810851   0.73220426]]
mean=[52.53235793 51.86342693]
covariance
[[244.04297958 193.42743647]
 [193.42743647 216.02097448]]
contribution
[0.92153752 0.07846248]
固有ベクトルの符号が反転していると左右・上下に反転がある
主成分分析のライブラリで計算したもの

第1主成分は，固有値が最も大きい固有ベクトルの方向です．固有値が最も大きい方向は，データの散らばりが最も大きくなる方向です．これは次のように導くことができます．


ベクトル\mathbf x_iの平均は\bar{\mathbf x}=\displaystyle\frac{1}{N}\sum_{i=1}^{N}\mathbf x_iです．
今，データのベクトル\mathbf x_iをデータの散らばりが最も大きくなる方向のベクトル\mathbf u_1に射影することを考えます．このとき，\mathbf u_1を単位ベクトルにとります．\mathbf x_iを\mathbf u_1に射影すると\mathbf x_i^T\mathbf u_1になります．射影の分散は
\sigma^2 = \displaystyle\frac{1}{N}\sum_{i=1}^{N}(\mathbf x_i^T\mathbf u_1 - \bar{\mathbf x}^T\mathbf u_1)^2 = \mathbf u_1^T S \mathbf u_1です．
ここで
S= \displaystyle\frac{1}{N}(\mathbf x_i - \bar{\mathbf x})(\mathbf x_i - \bar{\mathbf x})^Tです．
\mathbf u_1^T S \mathbf u_1を\mathbf u_1^T\mathbf u_1=1の条件の下で最大化します．これはラグランジェの未定乗数法を用いて
\mathbf u_1^T S \mathbf u_1 - \lambda(\mathbf u_1^T\mathbf u_1 - 1)
を最大化して解きます．\mathbf u_1で微分すると
S\mathbf u_1 - \lambda \mathbf u_1=0
になります．つまりS\mathbf u_1 = \lambda \mathbf u_1ですから，\mathbf u_1は行列Sの固有ベクトルです．
この式に左から\mathbf u_1^Tをかけると\mathbf u_1^T S \mathbf u_1 = \lambdaとなり，データ\mathbf x_iの射影の分散は固有値です．また，それを最大にするものは最大の固有値になります．