Python で 主成分分析 をやってみた
2018-09-04

今回は多変量解析手法の1つである 主成分分析を行います。

Principal Component Analysis の頭文字を取って PCA と呼ばれます。

主成分 (Principal component)

主成分とは データのばらつきを表す ベクトルです。

射影したデータの分散が大きい順に 第一主成分、第二主成分 ... とよび、 最大でデータの次元数かデータ数(のいずれか小さい方)と同じ数の主成分を持ちます。

Note

  • 次元数が主成分の最大となるのは 計算の過程で求める 共分散行列 のサイズが次元数となるからですが、 データ数によっても主成分が制限されるのは、データ数以上の主成分(直行するベクトル)は分散がほぼ0で意味のある値にならないためだそうです。 (nishio さん から教えていただきました。もし誤解してたら訂正ください!)

  • 第一主成分の傾きは主成分(ベクトル)とデータの距離合計が最小になるような直線です。

    • このように言うと 最小二乗法で求めた回帰直線とどう違うのかと思うんですが、 回帰直線の場合は Y方向の距離で、主成分の直線は 直交方向の距離を最小にするようです。

それぞれの主成分がどの程度データの特徴を表しているかを 主成分負荷量 といい、その比率を 寄与率 といいます。

寄与率の小さな主成分を排除することを次元削減といい、これが主成分分析を行う大きなメリットです。

固有値 (Eigen value)

まずは仕組みを理解するためにライブラリを使わずに固有値・固有ベクトルから主成分を求め、描画してみます。

流れとしては以下のようになります

  • 分散共分散行列 を求める

  • 求めた分散共分散行列の固有値と固有ベクトルを求める

    • 主成分ベクトルは 列方向に 並ぶので、今回は転置してから使う

In [1]:
%autosave 0
%matplotlib notebook
Autosave disabled
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import Matrix, init_printing
init_printing()
In [3]:
data = np.array([
    [172, 58, 63],
    [180, 80, 85],
    [175, 60, 70],
    [162, 50, 65],
    [167, 63, 70],
    [151, 45, 59],
    [188, 70, 90],
    [170, 70, 70],
    [166, 61, 66],
])
In [4]:
avg = np.array([np.mean(data[:, 0]), np.mean(data[:, 1]), np.mean(data[:, 2])])
avg
Out[4]:
array([170.11111111,  61.88888889,  70.88888889])
In [5]:
cov = np.cov(data.T)
Matrix(cov)
Out[5]:
$$\left[\begin{matrix}112.861111111111 & 90.7638888888889 & 95.0138888888889\\90.7638888888889 & 113.361111111111 & 88.7361111111111\\95.0138888888889 & 88.7361111111111 & 103.611111111111\end{matrix}\right]$$
In [6]:
# 固有値, 固有ベクトル
ev, evec = np.linalg.eig(cov)
evec = evec.T
Matrix(evec)
Out[6]:
$$\left[\begin{matrix}-0.58904357698053 & -0.576998859387866 & -0.565773789321389\\-0.634318860412278 & -0.103630316639308 & 0.766094211437801\\0.500666803116854 & -0.810143859902881 & 0.304958486549253\end{matrix}\right]$$
In [7]:
Matrix(ev)
Out[7]:
$$\left[\begin{matrix}293.029493323743\\12.9370765367566\\23.8667634728339\end{matrix}\right]$$
In [8]:
# 寄与率
crate = ev / np.sum(ev)
Matrix(crate)
Out[8]:
$$\left[\begin{matrix}0.888416856969407\\0.0392230718648507\\0.0723600711657421\end{matrix}\right]$$
In [9]:
# 累積寄与率
(crate[0], crate[0] + crate[2], crate[0] + crate[1] + crate[2])
Out[9]:
$$\left ( 0.8884168569694072, \quad 0.9607769281351493, \quad 0.9999999999999999\right )$$
In [10]:
# 描画のために寄与率順にソートする
evec = np.array([evec[0], evec[2], evec[1]])
Matrix(evec)
Out[10]:
$$\left[\begin{matrix}-0.58904357698053 & -0.576998859387866 & -0.565773789321389\\0.500666803116854 & -0.810143859902881 & 0.304958486549253\\-0.634318860412278 & -0.103630316639308 & 0.766094211437801\end{matrix}\right]$$
In [11]:
# それぞれの内積はほぼ0
evec[0].dot(evec[1]), evec[0].dot(evec[2]), evec[1].dot(evec[2])
Out[11]:
$$\left ( -2.498001805406602e-16, \quad 5.551115123125783e-17, \quad -5.551115123125783e-17\right )$$
In [12]:
fig = plt.figure()
ax = fig.gca(projection='3d', adjustable='box')
ax.axis('equal')
ax.set_xlabel('Height')
ax.set_ylabel('Weight')
ax.set_zlabel('Waist circumference')

distance = np.sqrt(((data.max(axis=0) - data.min(axis=0)) ** 2).sum())
right = distance / 2
left = -right

xs0 = np.array([evec[0][0] * left, evec[0][0] * right]) + avg[0]
ys0 = np.array([evec[0][1] * left, evec[0][1] * right]) + avg[1]
zs0 = np.array([evec[0][2] * left, evec[0][2] * right]) + avg[2]

xs1 = np.array([evec[1][0] * left, evec[1][0] * right]) + avg[0]
ys1 = np.array([evec[1][1] * left, evec[1][1] * right]) + avg[1]
zs1 = np.array([evec[1][2] * left, evec[1][2] * right]) + avg[2]

xs2 = np.array([evec[2][0] * left, evec[2][0] * right]) + avg[0]
ys2 = np.array([evec[2][1] * left, evec[2][1] * right]) + avg[1]
zs2 = np.array([evec[2][2] * left, evec[2][2] * right]) + avg[2]

ax.plot(xs0, ys0, zs0, 'r')
ax.plot(xs1, ys1, zs1, 'g')
ax.plot(xs2, ys2, zs2, 'b')

ax.scatter(data[:,0], data[:,1], data[:,2])
Out[12]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f25a847d198>

特異値分解することで共分散行列を用いずに求める方法もあるみたい ですが、それはまたの機会ということで。

scikit-learn

続いて、 scikit-learn の PCA クラスを使います。

pip でインストールできます

$ pip install scikit-learn
  • PCA のインスタンス作成時 n_components 引数に次元数を指定する

    • 今回は先程と同じ結果を求めたいので 次元数は 3 とする

  • 作成した PCA インスタンスに fit メソッドでデータを適用する

  • 求めるべき主成分ベクトルは components_ 属性に格納されている

これらの前提を踏まえ先程と同じグラフを描画します。

In [1]:
%autosave 0
%matplotlib notebook
Autosave disabled
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import Matrix, init_printing
from sklearn.decomposition import PCA
init_printing()
In [3]:
data = np.array([
    [172, 58, 63],
    [180, 80, 85],
    [175, 60, 70],
    [162, 50, 65],
    [167, 63, 70],
    [151, 45, 59],
    [188, 70, 90],
    [170, 70, 70],
    [166, 61, 66],
])
In [4]:
avg = np.array([np.mean(data[:, 0]), np.mean(data[:, 1]), np.mean(data[:, 2])])
avg
Out[4]:
array([170.11111111,  61.88888889,  70.88888889])
In [5]:
pca = PCA(n_components=3)
pca.fit(data)
Out[5]:
PCA(copy=True, iterated_power='auto', n_components=3, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [6]:
Matrix(pca.components_)
Out[6]:
$$\left[\begin{matrix}-0.58904357698053 & -0.576998859387866 & -0.565773789321389\\0.500666803116852 & -0.810143859902881 & 0.304958486549255\\0.63431886041228 & 0.103630316639305 & -0.7660942114378\end{matrix}\right]$$
In [7]:
# それぞれの内積はほぼ0
(
    pca.components_[0].dot(pca.components_[1]),
    pca.components_[0].dot(pca.components_[2]),
    pca.components_[1].dot(pca.components_[2])
)
Out[7]:
$$\left ( 0.0, \quad -5.551115123125783e-17, \quad -2.7755575615628914e-17\right )$$
In [8]:
fig = plt.figure()
ax = fig.gca(projection='3d', adjustable='box')
ax.axis('equal')
ax.set_xlabel('Height')
ax.set_ylabel('Weight')
ax.set_zlabel('Waist circumference')

distance = np.sqrt(((data.max(axis=0) - data.min(axis=0)) ** 2).sum())
right = distance / 2
left = -right

xs0 = np.array([pca.components_[0][0] * left, pca.components_[0][0] * right]) + avg[0]
ys0 = np.array([pca.components_[0][1] * left, pca.components_[0][1] * right]) + avg[1]
zs0 = np.array([pca.components_[0][2] * left, pca.components_[0][2] * right]) + avg[2]

xs1 = np.array([pca.components_[1][0] * left, pca.components_[1][0] * right]) + avg[0]
ys1 = np.array([pca.components_[1][1] * left, pca.components_[1][1] * right]) + avg[1]
zs1 = np.array([pca.components_[1][2] * left, pca.components_[1][2] * right]) + avg[2]

xs2 = np.array([pca.components_[2][0] * left, pca.components_[2][0] * right]) + avg[0]
ys2 = np.array([pca.components_[2][1] * left, pca.components_[2][1] * right]) + avg[1]
zs2 = np.array([pca.components_[2][2] * left, pca.components_[2][2] * right]) + avg[2]

ax.plot(xs0, ys0, zs0, 'r')
ax.plot(xs1, ys1, zs1, 'g')
ax.plot(xs2, ys2, zs2, 'b')

ax.scatter(data[:,0], data[:,1], data[:,2])
Out[8]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f7ba84fea20>

主成分得点 (Scores)

求めた主成分を新たな軸として、データの新たな座標を求めます。

主成分得点を求めるためには pca.transform メソッドを使います。

これは元データを主成分で一次変換しているのとほぼ同じです。

ほぼ というのは、transform では 0 (原点) がデータの中心になるので、 一次変換したデータから平均値をマイナスすることで同一の結果となりました。

今回は 第二主成分まで求めます。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import Matrix, init_printing
from sklearn.decomposition import PCA
init_printing()
In [3]:
data = np.array([
    [172, 58, 63],
    [180, 80, 85],
    [175, 60, 70],
    [162, 50, 65],
    [167, 63, 70],
    [151, 45, 59],
    [188, 70, 90],
    [170, 70, 70],
    [166, 61, 66],
])
In [4]:
pca = PCA(n_components=2)
pca.fit(data)
Out[4]:
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [5]:
Matrix(pca.components_)
Out[5]:
$$\left[\begin{matrix}-0.58904357698053 & -0.576998859387866 & -0.565773789321389\\0.500666803116852 & -0.810143859902881 & 0.304958486549255\end{matrix}\right]$$
In [6]:
transformed = pca.transform(data)  # fit_transform でもいい
Matrix(transformed)
Out[6]:
$$\left[\begin{matrix}5.59457314574721 & 1.69047980051002\\-24.2587737417006 & -5.41826398833493\\-1.28697182921984 & 3.7069018958996\\14.9694522120127 & 3.77487962166306\\1.69438020846081 & -2.72886410874386\\27.7285685916662 & 0.488513167596558\\-26.0300027102732 & 8.21330146837497\\-4.11174253819584 & -6.89787071871347\\5.70051666150263 & -2.82907713825197\end{matrix}\right]$$
In [7]:
plt.xlabel('pca1')
plt.ylabel('pca2')
plt.grid(True)
plt.plot(transformed[:,0], transformed[:,1], 'ro')
plt.show()
In [8]:
transformed2 = data.dot(pca.components_.T)
Matrix(transformed2)
Out[8]:
$$\left[\begin{matrix}-170.425177812395 & 58.3387309143345\\-200.278524699843 & 51.2299871254896\\-177.306722787362 & 60.3551530097241\\-161.050298746129 & 60.4231307354876\\-174.325370749681 & 53.9193870050807\\-148.291182366476 & 57.1367642814211\\-202.049753668415 & 64.8615525821995\\-180.131493496338 & 49.7503803951111\\-170.31923429664 & 53.8191739755726\end{matrix}\right]$$
In [9]:
# グラフの形は同じだが座標は違う
plt.xlabel('pca1')
plt.ylabel('pca2')
plt.grid(True)
plt.plot(transformed2[:,0], transformed2[:,1], 'ro')
plt.show()
In [10]:
transformed3 = data.dot(pca.components_.T)
transformed3 -= transformed3.mean(axis=0)
Matrix(transformed3)
Out[10]:
$$\left[\begin{matrix}5.59457314574723 & 1.69047980051002\\-24.2587737417006 & -5.41826398833494\\-1.28697182921982 & 3.70690189589961\\14.9694522120127 & 3.77487962166307\\1.69438020846079 & -2.72886410874386\\27.7285685916662 & 0.488513167596565\\-26.0300027102732 & 8.21330146837497\\-4.11174253819584 & -6.89787071871346\\5.70051666150263 & -2.82907713825197\end{matrix}\right]$$
In [11]:
plt.xlabel('pca1')
plt.ylabel('pca2')
plt.grid(True)
plt.plot(transformed3[:,0], transformed3[:,1], 'ro')
plt.show()
In [12]:
pca.explained_variance_ratio_
Out[12]:
array([0.88841686, 0.07236007])

画像の主成分分析 (Images)

最後に画像の主成分を求めそれを描画してみます。

画像における寄与率の低い主成分を削除した(後に次元を戻した)らどのように描画されるのかを見てみます。

In [1]:
%autosave 0
%matplotlib inline
Autosave disabled
In [2]:
from PIL import Image
from IPython.display import display
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import fetch_olivetti_faces
In [3]:
faces = fetch_olivetti_faces()
In [4]:
faces.images.shape
Out[4]:
(400, 64, 64)
In [5]:
# flat にしたものが data に入っている
faces.data.shape
Out[5]:
(400, 4096)
In [6]:
fig = plt.figure(figsize=(16, 6))
for i in range(30):
    ax = fig.add_subplot(3, 10, i + 1, xticks=[], yticks=[])
    ax.imshow(faces.images[i], cmap=plt.cm.bone)
In [7]:
# 1主成分
pca1 = PCA(n_components=1)
pca1.fit(faces.data)
Out[7]:
PCA(copy=True, iterated_power='auto', n_components=1, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [8]:
# 2主成分
pca2 = PCA(n_components=4)
pca2.fit(faces.data)
Out[8]:
PCA(copy=True, iterated_power='auto', n_components=4, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [9]:
# 4主成分
pca4 = PCA(n_components=4)
pca4.fit(faces.data)
Out[9]:
PCA(copy=True, iterated_power='auto', n_components=4, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [10]:
# 16主成分
pca16 = PCA(n_components=16)
pca16.fit(faces.data)
Out[10]:
PCA(copy=True, iterated_power='auto', n_components=16, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [11]:
# 32主成分 
pca32 = PCA(n_components=32)
pca32.fit(faces.data)
Out[11]:
PCA(copy=True, iterated_power='auto', n_components=32, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [12]:
# 64主成分
pca64 = PCA(n_components=64)
pca64.fit(faces.data)
Out[12]:
PCA(copy=True, iterated_power='auto', n_components=64, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [13]:
# 400主成分(データ数により、これがMAX)
pca400 = PCA(n_components=400)
pca400.fit(faces.data)
Out[13]:
PCA(copy=True, iterated_power='auto', n_components=400, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [14]:
# 主成分とそれぞれの寄与率を見ていく
display(pd.DataFrame(np.array([
    pca400.explained_variance_ratio_ * 100,
    pca400.explained_variance_ratio_.cumsum() * 100,
]).T, columns=['寄与率', '累積寄与率']))
寄与率 累積寄与率
0 2.381274e+01 23.812737
1 1.399396e+01 37.806702
2 7.968611e+00 45.775311
3 4.998332e+00 50.773643
4 3.609843e+00 54.383488
5 3.156932e+00 57.540417
6 2.426832e+00 59.967251
7 2.036401e+00 62.003647
8 1.958113e+00 63.961761
9 1.672123e+00 65.633888
10 1.595223e+00 67.229111
11 1.436979e+00 68.666092
12 1.246742e+00 69.912834
13 1.147133e+00 71.059967
14 1.062877e+00 72.122841
15 9.777201e-01 73.100563
16 9.190605e-01 74.019623
17 8.155730e-01 74.835190
18 7.538748e-01 75.589066
19 7.469854e-01 76.336052
20 6.985782e-01 77.034630
21 6.146284e-01 77.649261
22 5.839278e-01 78.233185
23 5.697224e-01 78.802910
24 5.461911e-01 79.349106
25 5.318553e-01 79.880966
26 5.138166e-01 80.394783
27 4.958412e-01 80.890617
28 4.576644e-01 81.348282
29 4.411840e-01 81.789467
... ... ...
370 4.089816e-03 99.913818
371 4.022978e-03 99.917839
372 3.962825e-03 99.921806
373 3.902175e-03 99.925705
374 3.869481e-03 99.929581
375 3.802802e-03 99.933380
376 3.788238e-03 99.937172
377 3.648663e-03 99.940819
378 3.608338e-03 99.944427
379 3.552887e-03 99.947975
380 3.426009e-03 99.951401
381 3.378371e-03 99.954781
382 3.306739e-03 99.958092
383 3.286914e-03 99.961380
384 3.212103e-03 99.964592
385 3.179530e-03 99.967766
386 3.107459e-03 99.970871
387 3.040177e-03 99.973907
388 2.934402e-03 99.976845
389 2.803527e-03 99.979645
390 2.706007e-03 99.982353
391 2.600349e-03 99.984947
392 2.490905e-03 99.987442
393 2.400202e-03 99.989845
394 2.204693e-03 99.992050
395 2.176392e-03 99.994225
396 2.141459e-03 99.996361
397 2.002386e-03 99.998367
398 1.639758e-03 100.000015
399 1.229380e-10 100.000015

400 rows × 2 columns

In [15]:
# 平均顔
plt.imshow(pca64.mean_.reshape(faces.images[0].shape), cmap=plt.cm.bone)
# ただの平均をとったのと同じ
# plt.imshow(faces.data.mean(axis=0).reshape(faces.images[0].shape), cmap=plt.cm.bone)
Out[15]:
<matplotlib.image.AxesImage at 0x7f77e51d2898>
In [16]:
# 以下では 基本的に 9割程度説明できる64の主成分を使っていく
pca64.components_.shape
Out[16]:
(64, 4096)
In [17]:
# 主成分の顔を描画してみる
fig = plt.figure(figsize=(16, 6))
for i in range(30):
    ax = fig.add_subplot(3, 10, i + 1, xticks=[], yticks=[])
    ax.imshow(pca64.components_[i].reshape(faces.images[0].shape),
              cmap=plt.cm.bone)
In [18]:
# 主成分の平均顔(どうでもいいけど)
plt.imshow(pca64.components_.sum(axis=0).reshape(faces.images[0].shape), cmap=plt.cm.bone)
Out[18]:
<matplotlib.image.AxesImage at 0x7f77e5265080>
In [19]:
# 戦闘データを64次元に削減
transformed = pca64.transform(faces.data[:1])
transformed.shape
Out[19]:
(1, 64)
In [20]:
# このままでは描画できないので主成分(転地)で射影して次元を戻す
face1 = transformed.dot(pca64.components_)
In [21]:
# 描画 (だいぶ元の顔に近いかんじ)
plt.imshow(face1.reshape(faces.images[0].shape), cmap=plt.cm.gray)
Out[21]:
<matplotlib.image.AxesImage at 0x7f77e57fb668>