環境情報論　第11回 応用編：主成分分析1

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

環境情報論の集大成として，今回から2回に渡って，行列を用いたデータ解析の基礎である主成分分析を学びます。

線型代数の知識をガンガン使うので，正直最初は難しく感じると思います。

学部発展（大学院基礎）レベルなので，「へぇーこんなことで行列使うのかー」と気軽に学んでください。

いつもの準備

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv

宇都宮の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Utsunomiya_temp.csv

（出典は第3回に記述）

In [2]:
tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
utsu_temp = np.genfromtxt("Utsunomiya_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]
utsu = utsu_temp[:, 2]

# 今回は1990年から2019年の30年分のデータを用いる
tokyo = tokyo[(1990 <= y)*(y <= 2019)]
utsu = utsu[(1990 <= y)*(y <= 2019)]
m = m[(1990 <= y)*(y <= 2019)]
y = y[(1990 <= y)*(y <= 2019)]

# 気候値の計算
tokyoc= np.zeros((12))
utsuc= np.zeros((12))
for mm in range(1, 13): 
    tokyoc[mm-1] = np.nanmean(tokyo[m==mm], 0)
    utsuc[mm-1] = np.nanmean(utsu[m==mm], 0)

# 偏差の計算
tokyoa = np.zeros((tokyo.shape))
utsua = np.zeros((utsu.shape))
for yy in range(1990, 2020):
    for mm in range(1, 13):
        tokyoa[(y==yy)*(m==mm)] = tokyo[(y==yy)*(m==mm)] - tokyoc[mm-1]
        utsua[(y==yy)*(m==mm)] = utsu[(y==yy)*(m==mm)] - utsuc[mm-1]

# 今回は温暖化には興味がないのでデトレンド
tokyoa = signal.detrend(tokyoa)
utsua = signal.detrend(utsua)

主成分分析とは

たとえば，札幌・仙台・東京・大阪・福岡・那覇の気温偏差の時系列が与えられたとします。

これらをブレンドして「日本で最も目立つ変動」を抜き出して時系列を一つ作りたいとき，どういう配合にすれば良いのかを知るのが，主成分分析の目的です。

単に6地点の平均を出すだけだと，結果は地点の選び方に大きく依存してしまいますので，望ましくありません。

「日本代表」（＝日本の気温偏差の変動の「主成分」）を客観的に抜き出す方法を目指しましょう。

高校の復習(？): 分散と標準偏差

まず，分散（variance） という量を復習しましょう。分散とは，簡単にいうと変動の大きさを表す指標です。

いま，$N$ヵ月分の月別気温偏差の時系列データを

$$\vec{x} = (x_1, x_2, \cdots, x_N)$$

のように行ベクトルで書くとき，$\vec{x}$の分散 $\mathrm{var}(\vec{x})$の定義は成分の2乗の平均，つまり

$$\mathrm{var}(\vec{x}) := (x_1^2 + x_2^2 + \cdots + x_N^2)/N = |\vec{x}|^2/N$$

です。たとえば東京と宇都宮の気温偏差の分散はそれぞれ

In [3]:
print(np.var(tokyoa))
print(np.var(utsua))

1.0319958222914387
1.0844095840574333

で求められます。単位は℃$^2$です。内陸にある宇都宮の方が分散がわずかに大きく，年々の寒暖差が大きいことがわかります。

ちなみに，単位が℃$^2$ではわかりづらいので，無理矢理ルートをとって単位を合わせたものが 標準偏差(standard deviation) です。

$$\mathrm{std}(\vec{x}) := \sqrt{\mathrm{var}(\vec{x})} = |\vec{x}|/\sqrt{N}$$
In [4]:
print(np.std(tokyoa))
print(np.std(utsua))

1.015871951720018
1.0413498855127576

この結果から，宇都宮の方が平均的に0.03℃くらい気温偏差の振幅が大きいことがわかります。

まずは「関東代表」の時系列を作ってみる

分散が最大となる方向に座標を回転

準備が終わったので，いよいよ本題です。

我々の目的は6地点の気温偏差をブレンドして 「日本代表」の時系列 を一つ作りたいということでした。

しかしいきなり6地点から始まるのはしんどいので，まずは第8回でやったように2変数のときから考えてみましょう。

2変数から 「関東代表」の時系列 を決めるのは簡単そうです（正直どうやってもできそうです）。

ただ，それを6変数に拡張することを見越して，ここではちょっと複雑なことをやります。

In [5]:
def draw_scatter(x, y, xname, yname):
    plt.scatter(x, y)
    plt.xlim(-5, 5)
    plt.ylim(-5, 5)
    plt.xlabel(xname)
    plt.ylabel(yname)
    plt.gca().set_aspect('equal', adjustable='box') # 縦横比を1:1にする

draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')
plt.axhline(y=0, xmin=-10, xmax=10, color='k')
plt.axvline(x=0, ymin=-10, ymax=10, color='k')
plt.show()

東京と宇都宮の気温偏差の散布図には，見覚えがありますね（第8回）。

いま，この散布図を2次元の「データ空間」に点が散らばっているという風に見てみましょう。

そうすると，データの散らばりが大きいのはどう見ても斜め約45˚に傾いた方向の変動で，これこそが今回見つけたい変動の主成分です。

これを少し数学的に言い換えると，

・与えられたデータ空間の$xy$座標を回転することによって新しい$x'y'$座標をつくり，データの$x'$成分の分散が最大になるようにしたい

・今回のデータでは，$x$軸を斜め約45˚方向に回転させて作った$x'$軸の$x'$成分の分散が最大になりそう

ということになります。

In [6]:
draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')

# x成分の分散が最大になるようにx軸を取り直すとだいたいこんな感じになりそう
x = np.array([-10, 10])
plt.plot(x, x, 'k')
plt.plot(x, -x, 'k')
plt.show()

主成分分析は，このように座標軸を回転させる（線型代数の言葉でいうと「基底を取り直す」）ことで，分散が最大になる方向を探して新しく軸を取り直す分析手法のことです。

元の座標では，

・$x$軸：東京の気温偏差の変動

・$y$軸：宇都宮の気温偏差の変動

のセットだったわけでした。

それに対して，座標回転後の$x'y'$座標では，

・$x'$軸：関東で最も目立つ気温偏差の変動（「第1主成分」という）

・$y'$軸：$x'$軸と相関のない（＝直交した）残りの変動（「第2主成分」という）

のセットになっています。

それでは，$x'$軸はどうやって探せば良いのでしょうか？

共分散行列

分散最大方向を客観的に見つけるには，「共分散行列」という行列の固有ベクトルを求めます（！）。

まさかこんなところでも役に立つとは線型代数シリーズPart3くらいでしょうか。

余談：2変数のときは「それって回帰直線の方向じゃないの？」と言いたくなります。実際ほぼその通りなのですが，ここでは多変数に応用することを見越した方法を勉強しています。

いま，$N$ヵ月分の東京の気温偏差時系列を$\vec{x_1}$，宇都宮の気温偏差時系列を$\vec{x_2}$とします。

そして，それらの行ベクトルを縦に並べた2行$N$列の行列をデータ行列$X$と呼ぶことにします。

$$X := \left( \begin{array}{c} \vec{x_1}\\ \vec{x_2} \end{array} \right) = \left( \begin{array}{cccc} x_{1, 1} & x_{1, 2} & \cdots & x_{1, N}\\ x_{2, 1} & x_{2, 2} & \cdots & x_{2, N} \end{array} \right)$$
In [7]:
N = tokyoa.shape[0]
X = np.array([tokyoa, utsua])

次に，共分散行列という行列を計算します。

そもそも共分散 (covariance)とは，時系列$\vec{x_1}$と$\vec{x_2}$がどれだけ一緒に大きく変動しているかを表す指標です。

共分散 $\mathrm{cov}(\vec{x_1}, \vec{x_2})$の定義は各成分同士の積の平均，つまり

$$\mathrm{cov}(\vec{x_1}, \vec{x_2}) := (x_{1,1}x_{2,1} + x_{1,2}x_{2,2} + \cdots + x_{1,N}x_{2,N})/N = (\vec{x_1}\cdot\vec{x_2})/N$$

です。次に共分散行列$C$の定義は，

$$C := \left( \begin{array}{cc} \mathrm{cov}(\vec{x_1}, \vec{x_1}) & \mathrm{cov}(\vec{x_1}, \vec{x_2})\\ \mathrm{cov}(\vec{x_2}, \vec{x_1}) & \mathrm{cov}(\vec{x_2}, \vec{x_2}) \end{array} \right) = XX^T/N$$

です。東京と宇都宮の気温偏差の共分散行列は

In [8]:
C = X@X.T/N #行列の積は@で書ける

C
Out[8]:
array([[1.03199582, 0.99284395],
       [0.99284395, 1.08440958]])

となります。注意して欲しいのは，対角成分には東京と宇都宮の気温偏差の分散の値がそのまま入っているということです。

これは，同じ時系列同士の共分散は分散に等しいからです。

$$\mathrm{cov}(\vec{x}, \vec{x}) = \mathrm{var}(\vec{x})$$

また，非対角成分には二つ同じ値が入っています。これは，

$$\mathrm{cov}(\vec{x1}, \vec{x2}) = \mathrm{cov}(\vec{x2}, \vec{x1})$$

だからです。共分散行列は実対称行列になるというのはものすごく重要な性質です。

分散最大方向を見つける

いよいよ分散最大方向を見つけるには，共分散行列の固有ベクトルを計算します。

復習：行列 $C$の固有ベクトルとは，

$$C\vec{e} = \lambda \vec{e}$$

を満たすような0ベクトルでないベクトル$\vec{e}$のことです。また，スカラー$\lambda$は固有値と呼ばれます。

$n$次正方行列が逆行列を持つとき（rankが$n$のとき），0でない固有値は$n$個，対応する固有ベクトルは$n$組存在します。

固有ベクトルのスカラー倍はすべて固有ベクトルなので，データ解析の文脈では$\vec{e}$は単位ベクトルとします。

いま，固有値を並べた行列を$\Lambda$，固有ベクトルを並べた行列を$E$とします。いわゆる「対角化」です。

$$C\vec{e_1} = \lambda_1 \vec{e_1}, \ \ C\vec{e_2} = \lambda_2 \vec{e_2} \iff CE = E \Lambda \iff E^{-1}CE = \Lambda$$

このとき，共分散行列$C$は実対称行列なので，$E$は直交行列となり（復習：「線型代数学4」），

$$\vec{e_1}\cdot\vec{e_2} =0$$

となります（固有ベクトル同士は直交）。また，実対称行列の固有値は実数になります（復習：「線型代数学4」）。

2次正方行列の場合は，

$$ E = \left( \begin{array}{cc} \vec{e_1} & \vec{e_2} \end{array} \right)$$$$ \Lambda = \left( \begin{array}{cc} \lambda_1 & 0\\ 0 & \lambda_2 \end{array} \right)$$

です。ただし，$\lambda_1 \geq \lambda_2 > 0$とします。

In [9]:
# 固有値，固有ベクトルもnumpyなら楽勝！
[Lam, E] = np.linalg.eig(C)
# 固有値の大きい順にソート
index = Lam.argsort()[::-1]   
Lam = Lam[index]
E = E[:,index]
# 対角成分に固有値を並べる
L = np.diag(Lam) 

print(E)
print(L)

[[-0.69771535 -0.7163751 ]
 [-0.7163751   0.69771535]]
[[2.05139247 0.        ]
 [0.         0.06501294]]

このとき，求めたい分散最大方向は，最大固有値に対応する固有ベクトル$\vec{e_1}$の方向となっています（証明：D問題）。

方向が知りたいだけで，固有ベクトルの符号は任意なので，あとで図示するときの都合でマイナスをつけたい場合は自由につけます。

（ただし，あとで使うので行列$E$の中身も書き換えることを忘れずに。）

In [10]:
e1 = -E[:, 0]
e2 = E[:, 1]

E[:, 0] = e1
E[:, 1] = e2

$\vec{e1}$を赤い矢印，$\vec{e2}$を黄色い矢印で図示してみましょう。

In [11]:
def draw_arrow(vec, color):
    point = {
            'start': [0, 0],
            'end': vec
        }
    plt.annotate('', xy=point['end'], xytext=point['start'],
                arrowprops=dict(shrink=0, width=3, headwidth=8, 
                                headlength=10, connectionstyle='arc3',
                                facecolor=color, edgecolor='black'))

draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')
draw_arrow(e1, 'red')
draw_arrow(e2, 'yellow')

無事に分散最大方向を見つけることができました！

この固有ベクトルを基底として，新しい$x'y'$平面を張り，$x'$成分の時系列を求めれば，欲しかった「関東代表」の時系列となります。

具体的な手続きは，次週のお楽しみ。

課題11（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. 主成分分析の目的と手法を，自分の言葉で簡単に説明してください（答案は短くて良いです）。

B問題. 6次元のデータ行列について上記と同様のことを行い，第一主成分の固有ベクトル（最大の絶対値を持つ固有値に対応する固有ベクトル）を求めてください。ただし使用データは札幌，仙台，東京，大阪，福岡，那覇の気温偏差 http://web.is.ocha.ac.jp/~tsubasa/env_info/jp_tempa.npz とし，

N = tokyoa.shape[0]

X = np.array([sapporoa, sendaia, tokyoa, osakaa, fukuokaa, nahaa])

によって$6\times N$のデータ行列$X$を定義してください。

ヒント：$N$と$X$を定義した後は，ほとんど資料のままのことをやるだけです。意味を噛み締めながら進んでください。ただし，散布図とかベクトルは図示しなくて良いです。そもそも6次元空間なんて描けないし。

C問題. 実対称行列の固有値が実数になること，および固有ベクトルを並べた行列が直交行列になることを証明してください。

ヒント：線型代数学4で履修済みのはずです。集中演習の勉強になりますね！！（）

D問題. 共分散行列（2行2列でも良い）について，固有ベクトルの方向が分散最大の方向と一致していることの証明を調べて読み，証明の流れを簡単に記述してください。

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「環境情報論課題11，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。