環境情報論　第12回 応用編：主成分分析2

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

環境情報論もいよいよフィナーレです！

前回に引き続き，行列を用いたデータ解析の基礎である主成分分析を学びます。

今日は，探した基底の使い方についてです。

いつもの準備（思えばいっぱい勉強したものです）

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)

関数も用意しておきます。

In [3]:
# 散布図を描く関数
def draw_scatter(x, y, xname, yname):
    plt.scatter(x, y)
    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    plt.xlabel(xname)
    plt.ylabel(yname)
    plt.gca().set_aspect('equal', adjustable='box') # 縦横比を1:1にする

# ベクトルを描く関数
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'))
    
# 月別の時系列を描画する関数
def plot_mon_time(time_series, lower = -6, upper = 6, init_year=1990, fin_year=2020):
    mon = np.arange(1990, 2020, 1/12)
    plt.figure;
    plt.plot(mon, time_series)
    plt.plot(mon, 0*time_series, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()

-----いつもの準備ここまで-----

主成分分析の流れ（前回の復習）

まず，先週やったことを一気にやりつつおさらいします。

In [4]:
# データ行列の定義
N = tokyoa.shape[0] #データ行列の大きさ
X = np.array([tokyoa, utsua]) #データ行列

# 共分散行列の計算
C = X@X.T/N 

# 共分散行列の固有値固有ベクトルを求める
[Lam, E] = np.linalg.eig(C)
# 固有値の大きい順にソート
index = Lam.argsort()[::-1]   
Lam = Lam[index]
E = E[:,index]
# 対角成分に固有値を並べる
L = np.diag(Lam)

# 固有ベクトルを抜き出すと，それが分散最大方向とその直交方向の基底となる
#（固有ベクトルの符号は自由につけてよいがEの方も直す）
e1 = -E[:, 0] #第1主成分の基底
e2 = E[:, 1] #第2主成分の基底
E[:, 0] = e1
E[:, 1] = e2

これによって，分散最大方向の単位ベクトル$\vec{e_1}$（赤）とその直交方向の単位ベクトル$\vec{e_2}$（黄色）を見つけたのが前回までの内容でした。

In [5]:
draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)') # 前回の図よりもちょっと拡大して表示してます
draw_arrow(e1, 'red')
draw_arrow(e2, 'yellow')

$x'$座標（e1への射影「PC1」）を計算する

いま，基底$\vec{e_1}$の向きを前回同様$x'$軸と呼ぶことにします。

たとえば，1990年2月の点を表すベクトル$\vec{a}$（以下ではピンク矢印で描画することにします）について，$x'$座標を求めるにはどうすれば良いでしょうか。

In [6]:
draw_scatter(tokyoa, utsua, 'Tokyo (℃)', 'Utsunomiya (℃)')
x = np.array([-10, 10])
plt.plot(x, x, 'k')
draw_arrow(e1, 'red') # 基底e1を赤で描画

a = X[:, 1] # 1990年2月の点
draw_arrow(a, 'pink') # ベクトルaをピンクで描画

plt.plot(x, -(x-a[0])+a[1], 'k--') # x'座標がわかりやすいように点線を描画しているだけ
plt.show()

$x'$座標を求めるには，基底とデータ点を表すベクトルとの内積をとります（ここで，$\vec{e_1}$を単位ベクトルにしておいたことが効いてきます）。

このことを，$\vec{e_1}$に射影すると言います。$\vec{e_1}$（赤いベクトル）と$\vec{a}$（ピンク色のベクトル）のなす角を$\theta$とすると，

$$ x' = |\vec{a}|\cos\theta = \vec{e_1} \cdot \vec{a} = \vec{e_1}^T \vec{a}= \left( \begin{array}{cc} e_{1,x} \ \ e_{1,y} \end{array} \right) \left( \begin{array}{c} {a_1}\\ {a_2} \end{array} \right) $$
In [7]:
x_prime = np.dot(e1, a) # e1とaの内積を計算する

x_prime
Out[7]:
2.1209285142499463

こうすると，1990年2月の$x'$座標が求まりました！

では，同様に他の時刻についても計算したいときは，各月についてのfor文を回さなければいけないのでしょうか？

実はそんなことはありません。一発で計算する方法があります。行列の積です！

$$\vec{x'}= \vec{e_1}^T X = \left( \begin{array}{cc} e_{1,x} \ \ e_{1,y} \end{array} \right) \left( \begin{array}{c} \vec{x_1}\\ \vec{x_2} \end{array} \right) = \left( \begin{array}{cc} e_{1,x} \ \ e_{1,y} \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) = \left( \begin{array}{cccc} x'_{1} & x'_{2} & \cdots & x'_{N} \end{array} \right) $$

こうすることで，全ての月において$x'$座標を求めた時系列が求まります。

この時系列のことを，第一主成分（First Principal Component）の時系列，または単にPC1（時系列）といいます。

In [8]:
pc1 = e1.T@X

# PC1を描画する
plot_mon_time(pc1)

これが，求めたかった「関東代表」の時系列です！$\vec{e_1}$の中身は

In [9]:
print(e1)

[0.69771535 0.7163751 ]

ですから，PC1（すなわち$x'$座標）が正のときは東京も宇都宮も同程度に暖かく，逆にPC1が負のときは東京も宇都宮も同程度に寒かったことになります。

$y'$座標（e2への射影「PC2」）を計算する

同様に$\vec{e_2}$にデータ行列$X$を射影することで，$x'$座標と直交する$y'$座標も求めてみます。これが，第二主成分の時系列またはPC2（時系列）です。

In [10]:
pc2 = e2.T@X

# PC1を描画する
plot_mon_time(pc2)

あえて縦軸はPC1のときと同じにしてありますが，明らかにPC1よりも分散が小さいです。

これは，データからPC1の変動を取り除いてしまったあとの「出がらし」であるPC2には，あまり情報が残っていないことを意味します。

$\vec{e_2}$の中身は

In [11]:
print(e2)

[-0.7163751   0.69771535]

ですから，PC2（すなわち$y'$座標）が正のときは東京は寒く宇都宮が暖かいこと，逆にPC2が負のときは東京が暖かく宇都宮が寒いことを意味しています。

実際の$i$ヶ月目の気温（$x_{1, i}$が東京，$x_{2, i}$が宇都宮）は，各月のPC1とPC2の重ね合わせで戻すことができます。

$$\left( \begin{array}{c} x_{1, i}\\ x_{2, i} \end{array} \right) = PC1_i \ \vec{e_1} + PC2_i \ \vec{e_2} \ \ \ \ \ (i = 1, 2, ..., N) \iff X = \vec{e_1}\vec{PC1} + \vec{e_2} \vec{PC2}$$

（右の式の右辺は，行列としての積です）

つまり主成分分析をすることによって，元の気温偏差時系列を，第1項の「東京と宇都宮が同位相で変化する第1主成分」（振幅が大きい），第2項の「東京と宇都宮が逆位相で変化する第2主成分」（振幅が小さい）に分解できたということになります。

In [12]:
plot_mon_time(tokyoa, -4, 4) # 元の東京の気温偏差時系列
plot_mon_time(pc1*e1[0] + pc2*e2[0], -4, 4) # 二つの主成分を重ね合わせて戻した時系列
plt.show() # 完 全 に 一 致

そしてなんと，PC1，PC2，...を一発で同時に求める方法があります！！

上記では射影の意味を理解してもらうために回りくどく計算しましたが，慣れたらいきなり以下の行列計算で終わらせましょう（データ行列の次元が増えれば増えるほど有効です）。

$$ Z = E^{T}X \\ (\iff \vec{PC1} = \vec{e_1}^T X, \vec{PC2} = \vec{e_2}^T X)$$

ただし，ここで$Z$はPC時系列を縦に並べた行列

$$Z := \left( \begin{array}{c} \vec{PC1}\\ \vec{PC2} \end{array} \right)$$

です。つまり，固有ベクトルを並べた行列自体で一気に射影を取れば良いのですね。

逆に，分解した主成分から元のデータに戻すのも簡単（EZ，イーズィー）で（僕の師匠の渾身のギャグです），

$$X=EZ$$

です（$E$は直交行列なので，$E^{-1} = E^{T}$であり，これを$Z = E^{T}X$に左から掛ければ良い）。

In [13]:
Z = E.T@X #行列の積は@で書ける
pc1 = Z[0, :] # せやかて工藤，実はpc1もpc2もこれだけであっさり求められるんや
pc2 = Z[1, :] # 今まで黙ってて悪かった，許したってーな

主成分分析の性質1: 二重直交性

主成分分析の大事な性質として，二重直交性（double orthogonality）があります。

すなわち，固有ベクトル同士・PC時系列同士はそれぞれ直交しているということです。

In [14]:
np.dot(e1, e2) # 固有ベクトル同士の内積はゼロ（前回のC問題）
Out[14]:
0.0
In [15]:
np.dot(pc1, pc2) # PC時系列同士の内積もゼロ
Out[15]:
-7.149836278586008e-14

すなわち，PC時系列同士は無相関であるということです。分散の大きい順に，互いに相関のない成分に分解するのが主成分分析であるともいえます。

In [16]:
np.corrcoef(pc1, pc2)[1, 0]
Out[16]:
-5.50593451171466e-16

PC時系列同士が直交することの証明：$Z$の共分散行列が対角行列になることを示せば良い。 $$ZZ^{T}/N = (E^{T}XX^{T}E)/N = E^{T}CE = E^{-1}CE = \Lambda \ \ \ ▪️$$

主成分分析の性質2: 固有値は対応する固有ベクトル方向の分散

次に，PC1とPC2の分散を計算してみましょう。

In [17]:
print(np.var(pc1))
print(np.var(pc2))

2.0513924657032288
0.06501294064564263

PC1の分散が最大になるように基底を回転したのですから，PC1の方が圧倒的に大きくなっています。

ここで，固有値を並べた行列$\Lambda$の中身をもう一度見てみましょう。

In [18]:
print(L)

[[2.05139247 0.        ]
 [0.         0.06501294]]

このように，実は$\Lambda$の対角成分には，各PC時系列の分散が入っています。

主成分分析の性質3: 分散の和が保存される

さらに，上記の行列$\Lambda$と，元データの共分散行列$C$を比べて，何か気づくことはありますでしょうか？

In [19]:
print(C)

[[1.03199582 0.99284395]
 [0.99284395 1.08440958]]

実は，$C$と$\Lambda$の対角成分の和（トレース）は一致しています。

復習：基底変換に関してトレースは不変で，固有値の和と等しい。

これを，データの視点から解釈し直すと，

東京の気温偏差の分散 + 宇都宮の気温偏差の分散 = PC1の分散 + PC2の分散

という関係が成り立っていることがわかります。つまり主成分分析は，

・分散の和は不変とする

・二重直交性を持たせる

という制約のもとで，出来る限り上位の主成分に分散を押し付ける解析方法であると言えます。

「各PCが，データ内の分散の和のうちどのくらいの割合を説明できているのか」のことを，寄与率（explained variance）といいます。

つまり一般に，$M$個の時系列に対して主成分分析を行った場合，$i = 1, 2, ..., M$に対して，PC$i$の寄与率EV$i$は

$$EVi:= \frac{\mathrm{var}(\vec{PCi})}{\sum_{j = 1}^M \mathrm{var}(\vec{x_j})} = \frac{\lambda_i}{\sum_{j=1}^{M}\lambda_j}$$

今回の例でも，PC1の寄与率EV1を求めてみると，

In [20]:
ev1 = L[0, 0]/np.trace(L) # トレースの計算はnp.traceを用いる

ev1
Out[20]:
0.9692814332969409

つまり，東京と宇都宮における気温偏差の分散の和のうち，97%が第一主成分のみによって説明されるということになります。

（要は「東京と宇都宮の気温変動はほとんど同位相」という見たまんまの結論ですが，この方法をそのまま$M$次元に拡張できるというのがパワフルな点です。）

寄与率が大きな主成分が得られたならば，データ内の大事な変動成分はほぼその主成分に集約できたことになります（「次元削減」と「特徴抽出」）。

主成分分析は大成功です！

余談：気象学や海洋物理学では，主成分分析のことを経験的直交関数(empirical orthogonal function; EOF)解析と呼びます。「経験的」はデータから基底を求めるという意味で，直交は基底が直交しているからです。PC1, PC2, ...に対応する固有ベクトルのことを，それぞれEOF1，EOF2, ...などと呼びます。大気海洋データ解析業界の方言のようなものです。

課題12（締切：1月末までということにしてみますが，他の課題とかテストが忙しかったら多少遅れても良いです。ただし，神山研の配属を考えてくれている方は，1月中旬の配属順位決定までに提出された課題（他の回のC問題・D問題も含む）のみを考慮しますので注意してください。）

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

A問題. 以下の二つについて答えてください。

  1. PC時系列，および寄与率とは何か，自分の言葉で簡単に説明してください（答案は短くて良い）。

  2. 環境情報論全体を通して，感想を自由にお願いします。

B問題. 前回のB問題の続きを上記資料にしたがって行い，次の問いに答えてください。

  1. e1 = E[:, 0] の中身を見ると，各都市はどのような「ブレンド」になっているでしょうか。
  2. 「日本代表」の時系列（PC1時系列）を求めてください。

  3. 第一主成分（$\vec{e_1}$とPC1時系列のセット）は，概ねどのような変動を意味しているでしょうか（ヒント：上記の例では，「東京と宇都宮が同位相に変動する成分」を表していました）。

  4. また，得られたPC1時系列は札幌，仙台，東京，大阪，福岡，那覇の気温偏差の分散の和のうち何パーセントを説明できているでしょうか。

  5. 第二主成分についても同様に解釈を加えてください。

C問題. B問題において，元のデータ行列は

$$X = \vec{e_1}\vec{PC1} + \vec{e_2}\vec{PC2} + ... + \vec{e_6}\vec{PC6}$$

のように，寄与率の大きい順に6成分に分解できたことになります。いま，寄与率の大きい上位3成分のみを残し，残りは捨てることにすると，元のデータ行列$X$に含まれる変動の情報量はどの程度残っていると言えるでしょうか（ヒント：上位3成分の寄与率の合計を求めれば良く，これを累積寄与率といいます）。

また，上位3成分のみから復元した各都市の時系列と，元データの各都市の時系列を比較し，次元削減がうまくいっているかを確かめてください。

D問題. 環境情報論で学んだ解析方法を用いて，好きなデータについて何か面白いレポートを書いてください。

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