主成分分析 (1)

12. 主成分分析 (1)#

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation
from IPython.display import HTML

\(\def\bm{\boldsymbol}\)以下は、教育用標準データセット（SSDSE: Standardized Statistical Data Set for Education）の「C. 都道府県庁所在市別、家計消費データ」から麺類に関する部分を抜き出したものである（多いので冒頭の15行だけを表示した）。SSDSE-Cは、総務省統計局「家計調査」から、二人以上の世帯の１世帯当たりの年間支出金額を都道府県庁所在市別・品目別に表形式で収録したものである。

Hide code cell source 
df = pd.read_excel('https://www.nstac.go.jp/SSDSE/data/2021/SSDSE-C-2021.xlsx', skiprows=1)
df = df[['都道府県', '生うどん・そば', "乾うどん・そば", "パスタ", "中華麺", "カップ麺", "即席麺", "他の麺類"]]
df = df.iloc[1:]
df[:15]
都道府県 生うどん・そば 乾うどん・そば パスタ 中華麺 カップ麺 即席麺 他の麺類
1 北海道 3162 2082 1266 4152 5189 1609 726
2 青森県 2964 2224 1114 5271 6088 2039 554
3 岩手県 3349 2475 1305 5991 5985 1889 898
4 宮城県 3068 2407 1339 5017 5295 1709 842
5 秋田県 3231 3409 1019 5172 5013 1680 554
6 山形県 4478 3084 1288 5236 5875 1745 754
7 福島県 2963 2705 1064 4397 5862 1687 919
8 茨城県 3353 2477 1248 4034 4562 1440 653
9 栃木県 3908 2218 1391 4534 4945 1860 742
10 群馬県 4563 1948 1203 4153 5049 1544 546
11 埼玉県 4016 2256 1487 4512 4584 1568 949
12 千葉県 3389 2277 1441 4582 4513 1840 835
13 東京都 3088 2385 1595 4592 3953 1734 974
14 神奈川県 3410 2684 1472 4484 4035 1741 1049
15 新潟県 2697 3409 1308 4192 6095 2201 612

このデータを分析することで、日本における麺類の消費行動の全体的な傾向や、地域における特徴が明らかになるかもしれない。また、カップ麺と即席麺は消費行動が似ているかもしれないのでどちらか一方だけ調査すれば十分である、ということがデータから明らかになるかもしれない。

そこで、麺類の消費行動の傾向を掴むために、このデータを可視化することを考える。仮に、このデータが2次元（2種類の麺類によるデータ）であるならば、2次元平面上に各都道府県庁所在地をプロットすればよい。ところが、このデータは各都道府県庁所在地が7次元のベクトルで表されているため、2次元平面上にプロットすることができない。そこで、このデータの特徴空間を7次元から2次元に圧縮することを考えたい。

本章では、このような分析に役立つ手法として、主成分分析（principal component analysis）を紹介する。主成分分析は、データの情報をできるだけ豊富に表現できるように配慮しながら、もとの特徴量の線形結合で表現される新しい「軸」を導き出す手法である。元データの情報量を最も豊富に表現できる軸は第１主成分（first principal component）と呼ばれる。その第１主成分と直交し、かつ２番目に情報量が豊富な軸を第２主成分、第１主成分と第２主成分に直交し、かつ３番目に情報量が豊富な軸を第３主成分、･･･という形で、元データを表現する情報量の多い順にお互いに直交する複数の軸を求めることができる。

元のデータの各事例は、主成分分析で求められた新しい軸に射影される。主成分分析で求められた全ての軸を用いることで、元のデータを完全に再現することができる。また、元データを第１主成分、第２主成分といった主要な軸だけで表現することで、元データよりも少ない軸に要約することができる。このため、主成分分析はデータの次元圧縮やノイズ除去、可視化等に用いられる。

12.1. ２次元データの例#

いきなり3次元以上のデータで考えるのは難しいので、まずは2次元のデータの場合を考える。簡単な例として、以下のグラフで示されたデータ\(\mathcal{D}_s\)を考える。このデータは\(xy\)平面上の\(4\)点\((-7,-2),(-3,-3),(4,1),(6,4)\)で構成されている。

Hide code cell source 
D = np.array([[-7, -2], [-3, -3], [4, 1], [6, 4]])
C = ['red', 'blue', 'green', 'orange']

fig = plt.figure(dpi=100, figsize=(6, 6))
ax = fig.add_subplot(1,1,1)
ax.scatter(D[:,0], D[:,1], marker='o', color=C)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
ax.set_xlim(-8, 8)
ax.set_ylim(-8, 8)
ax.grid()
../_images/b4a41716e49a12f044d3ae7c922925f48b1539f6e4dc522111ea2b8ddd400ff4.png

このデータ\(\mathcal{D}_s\)を次のように表すことにする（\(i \in \{1, 2, 3, 4\}\)は事例に対するインデックスである）。

(12.1)#\[\begin{align} \mathcal{D}_s = \{(x_i, y_i)\}_{i=1}^{4} = \{(-7,-2),(-3,-3),(4,1),(6,4)\} \end{align}\]

なお、

(12.2)#\[\begin{align} \frac{1}{4}\left(-7 -3 + 4 + 6\right) = 0,\; \frac{1}{4}\left(-2 -3 + 1 + 4\right) = 0 \\ \end{align}\]

であるから、この４点の重心は原点\((0,0)\)である（データの重心が原点にない場合の取り扱いは次章で説明する）。

12.2. データ点を直線へ射影#

さて、各点\((x_i, y_i) \in \mathbb{R}^2\)を\(xy\)平面上で原点を通る直線に射影し、その直線上の点\((\hat{x}_i, \hat{y}_i)\)で近似的に表現したい。原点を通る直線は、長さ\(1\)の方向ベクトル\(\bm{u} \in \mathbb{R}^2\)と媒介変数\(a \in \mathbb{R}\)を用いて、\(a \bm{u}\)で表される。したがって、点\((x_i, y_i)\)を射影した点\((\hat{x}_i, \hat{y}_i)\)は、\(a_i \in \mathbb{R}\)を用いて次式で表される。

(12.3)#\[\begin{align} \begin{pmatrix}\hat{x}_i \\ \hat{y}_i\end{pmatrix} = a_i \bm{u} \end{align}\]

なお、\(|a_i|\)は射影された点\((\hat{x}_i, \hat{y}_i)\)の原点からの距離を表す。射影された点が方向ベクトルと同じ向きにあるときは\(a_i\)は正の値、反対の向きにあるときは\(a_i\)は負の値をとることとし、\(a_i\)を直線上の位置と呼ぶ。以上の定式化により、元のデータ点\((x_i, y_i)\)は\(\bm{u}\)が定める（１次元の）数直線上の位置\(a_i\)として近似的に表現される。

例えば、直線\(y=\frac{1}{2}x\)を考えると、元のデータ点（丸印）は射影された点（データ点から降ろした垂線の足、四角印）で近似的に表現される。

Hide code cell source 
def draw_projections(ax, q, fill=False, show_a = False, show_e = False, show_xy=False):
    artist = []
    ya = ['bottom', 'top', 'top', 'bottom']
    yoffs = np.array([0.2, -0.2, -0.2, 0.2])
    axoffs = np.array([-0.2, -0.3, -0.4, -0.2])
    ayoffs = np.array([-0.7, 0.5, 0.5, -0.7])
    exoffs = np.array([-0.8, -0.3, 0.3, -0.1])
    eyoffs = np.array([0.35, -0.8, -0.4, 0.35])
    
    P = []
    for i in range(D.shape[0]):
        x, y = D[i,0], D[i,1]
        # Project (x, y) onto the line q.
        d = np.dot(np.array([x, y]), q)
        p = d * q
        artist.append(ax.scatter(p[0], p[1], marker='s', color=C[i]))
        if show_a:
            artist.append(ax.text(p[0] + axoffs[i], p[1] + ayoffs[i], '$a_{}$'.format(i+1)))
        if show_e:
            artist.append(ax.text(p[0] + exoffs[i], p[1] + eyoffs[i], '$\epsilon_{}$'.format(i+1)))            
        if show_xy:
            artist.append(ax.text(x, y + yoffs[i], '$(x_{0}, y_{0})$'.format(i+1), ha='center', va=ya[i]))
        artist += ax.plot([x, p[0]], [y, p[1]], ls='--', color=C[i])
        P.append((i, p)) 

    # Sort the projected points so that we can draw lines in an appropriate order.
    P.sort(key=lambda x: np.dot(x[1], x[1]), reverse=True)
    for i, p in P:
        artist += ax.plot([0, p[0]], [0, p[1]], ls='-', color=C[i])
        if fill:
            tri = plt.Polygon(((0, 0), D[i], p), fc=C[i], alpha=0.2)
            artist.append(ax.add_patch(tri))

    return artist

u = np.array([2, 1])
uy = 8 * u[1] / u[0]

fig = plt.figure(dpi=100, figsize=(6, 6))
ax = fig.add_subplot(1,1,1)
ax.scatter(D[:,0], D[:,1], marker='o', color=C)
ax.plot([-8, 8], [-uy, uy], ls='-', color='black', lw=1)
draw_projections(ax, u / np.linalg.norm(u), show_a=True, show_xy=True)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
ax.set_xlim(-8, 8)
ax.set_ylim(-8, 8)
ax.grid()
../_images/8f30e237f508f0d38b5407c0366fb6d657fef50e63012bb3ee45e75f1d8eee3c.png

\(a_i\)を求めるため、直線\(y=\frac{1}{2}x\)に対応する長さ\(1\)の方向ベクトル\(\bm{u}\)を求める。

(12.4)#\[\begin{align} \bm{u} = \frac{1}{\sqrt{2^2 + 1^2}}\begin{pmatrix}2 \\ 1\end{pmatrix} = \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} \end{align}\]

ある点\((x_i, y_i)\)を方向ベクトル\(\bm{u}\)上に射影したとき、その位置\(a_i\)は、

(12.5)#\[\begin{split} \begin{align} a_i = \left\|\begin{pmatrix}x_i \\ y_i\end{pmatrix}\right\| \cos \theta = \left\|\begin{pmatrix}x_i \\ y_i\end{pmatrix}\right\| \frac{\begin{pmatrix}x_i \\ y_i\end{pmatrix}^\top \bm{u}}{\left\|\begin{pmatrix}x_i \\ y_i\end{pmatrix}\right\| \|\bm{u}\|} = \begin{pmatrix}x_i \\ y_i\end{pmatrix}^\top \bm{u} \end{align} \end{split}\]

と求めることができる。ここで、\(\theta\)は位置ベクトル\((x_i, y_i)\)と方向ベクトル\(\bm{u}\)のなす角である。ゆえに、

(12.6)#\[\begin{align} a_1 &= \begin{pmatrix}-7 & -2\end{pmatrix} \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} \approx -7.16 \\ a_2 &= \begin{pmatrix}-3 & -3\end{pmatrix} \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} \approx -4.02 \\ a_3 &= \begin{pmatrix}4 & 1\end{pmatrix} \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} \approx 4.02 \\ a_4 &= \begin{pmatrix}6 & 4\end{pmatrix} \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} \approx 7.16 \\ \end{align}\]

以上により、データ\(\mathcal{D}_s\)の各事例\((x_i, y_i)\)を方向ベクトル\(\bm{u}\)と、位置\(a_i\)で近似的に表現したことになる。各事例\((x_i, y_i)\)を射影した点\((\hat{x}_i, \hat{y}_i)\)は、

(12.7)#\[\begin{align} \begin{pmatrix}\hat{x}_1 \\ \hat{y}_1\end{pmatrix} &= -7.16 \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} = \begin{pmatrix}-6.4 \\ -3.2 \end{pmatrix}\\ \begin{pmatrix}\hat{x}_2 \\ \hat{y}_2\end{pmatrix} &= -4.02 \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} = \begin{pmatrix}-3.6 \\ -1.8 \end{pmatrix}\\ \begin{pmatrix}\hat{x}_3 \\ \hat{y}_3\end{pmatrix} &= 4.02 \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} = \begin{pmatrix}3.6 \\ 1.8 \end{pmatrix} \\ \begin{pmatrix}\hat{x}_4 \\ \hat{y}_4\end{pmatrix} &= 7.16 \begin{pmatrix}\frac{2}{\sqrt{5}} \\ \frac{1}{\sqrt{5}}\end{pmatrix} = \begin{pmatrix}6.4 \\ 3.2 \end{pmatrix} \\ \end{align}\]

12.3. 近似の良さの定量化#

直線への射影による近似の良さを定量化するため、元データの各点\((x_i, y_i)\)から直線上に垂線を降ろしたときの長さを\(\epsilon_i\)とし、これを残差と呼ぶ。三平方の定理（ピタゴラスの定理）より、

(12.8)#\[\begin{align} x_i^2 + y_i^2 = a_i^2 + \epsilon_i^2 \end{align}\]

が成り立つ。以下のグラフでは、各事例に関して直角三角形を色付けして描いている。直線\(y=\frac{1}{2}x\)上の面を底辺、垂線を高さとすると、

  • \(x_i^2 + y_i^2\): 斜面の長さの二乗（原点から丸印までの距離の二乗）

  • \(a_i^2\): 底辺の長さの二乗（原点から四角印までの距離の二乗）

  • \(\epsilon_i^2\): 高さ（点線）の二乗（残差の二乗）

と整理できる。

Hide code cell source 
u = np.array([2, 1])
uy = 8 * u[1] / u[0]

fig = plt.figure(dpi=100, figsize=(6, 6))
ax = fig.add_subplot(1,1,1)
ax.scatter(D[:,0], D[:,1], marker='o', color=C)
ax.plot([-8, 8], [-uy, uy], ls='-', color='black', lw=1)
draw_projections(ax, u / np.linalg.norm(u), True, show_a=True, show_e=True, show_xy=True)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
ax.set_xlim(-8, 8)
ax.set_ylim(-8, 8)
ax.grid()
../_images/096ab10fd3d55f7a6af7e6dead9eca4eda7b219124a8b47979038546f27de87e.png

これまでの議論では、データ点を射影する直線（方向ベクトル\(\bm{u}\)）を固定していたが、以降は射影された点の良さを考慮しながら様々な直線を検討したい。今回のデータ\(\mathcal{D}_s\)の重心は原点であるため、直線の切片は\(0\)に固定し、直線の傾きだけを変更する。回帰分析と同様に、方向ベクトル\(\bm{u}\)で元のデータ点を射影したときの「良さ」を、全てのデータ点の残差の二乗和で定量化する。

(12.9)#\[\begin{align} \sum_{i=1}^{N} \epsilon_i^2 \end{align}\]

なお、2次元のデータにおいて、線形回帰では\(y\)軸に関する実際の値と近似値の値の差\(\epsilon = y - \hat{y}\)を残差として定義したのに対し、主成分分析では事例を直線に射影したときの垂線の長さ\(\epsilon = \sqrt{(x - \hat{x})^2 + (y - \hat{y})^2}\)を残差として定義する所が異なる。

さて、直線の傾きを変えながら射影された点をプロットするアニメーションを示す。

Hide code cell source 
K = 50
X = np.hstack([np.linspace(-8, 8, K, endpoint=False), np.full(K-1, 8)])
Y = np.hstack([np.full(K-1, -8), np.linspace(-8, 8, K, endpoint=False)])

fig = plt.figure(dpi=100, figsize=(6, 6))
ax = fig.add_subplot(1,1,1)
ax.scatter(D[:,0], D[:,1], marker='o', color=C)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_aspect('equal')
ax.set_xlim(-8, 8)
ax.set_ylim(-8, 8)
ax.grid()

artists = []
for x, y in zip(X, Y):
    q = np.array([-x, -y]) - np.array([x, y])
    artist = []
    artist += ax.plot([x, -x], [y, -y], ls='-', color='black', lw=1)
    artist += draw_projections(ax, q / np.linalg.norm(q), True)
    artists.append(artist)

ani = matplotlib.animation.ArtistAnimation(fig, artists, interval=100)
#ani.save('pca.mp4', writer="ffmpeg")
html = ani.to_jshtml()
plt.close(fig)
HTML(html)