1. EOF解析の実践（基礎編）¶
(A) 特異値分解を用いた方法¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
(A-1) データの読み込み、偏差行列の作成¶
サンプルデータ：
data/sample_data.csvを読み込む（30人の国語、数学、英語の点数）。以下のような偏差行列($X$)を作成する（バーは各科目の平均）。
In [2]:
# データの読み込み
file = "./sample_data.csv"
# --> file = "data/dd/sample_data.csv"に変更のこと
X_org = pd.read_csv(file,encoding='shift_jis',header=0,index_col=None)# headerは一行目、index列はなし、の意味
print(X_org.head()) #最初の5人分
# 偏差行列にする（データは三列あるが、まずは最初の二列（二科目）を解析してみる -> 後の課題1では三科目データを使う)
columns = X_org.columns
data = X_org.values[:,0:2] # .valuesで数値を取り出す（ラベルは剥がす）
mean = data.mean(axis=0) # 平均
anom = data - mean # 平均偏差行列
m,n = anom.shape # 配列の大きさ
(A-2) 特異値分解¶
$m \times n$ $(m >n)$の偏差行列$X$について：
$$ \begin{align} X &= \sum_{k=1}^n \sigma_k \boldsymbol{u}_k \boldsymbol{v}_k^T \\ & = UDV^T \end{align} $$$$ U = (\boldsymbol{u}_1,...,\boldsymbol{u}_m), \ D = \left[\begin{array}{ccc} \sigma_1&\cdots&O \\ \vdots&\ddots&\vdots \\ O&\cdots&\sigma_n \\ &O& \\ \end{array}\right], \ V^T = \left( \begin{array}{ccc} \boldsymbol{v}_1 \\ \vdots \\ \boldsymbol{v}_n \\ \end{array} \right) $$- ### 主成分ベクトル（EOF$k$）: $\boldsymbol{v}_k$
- ### 主成分得点（PC$k$）: $\sigma_k \boldsymbol{u}_k$
元の偏差行列$X$の$v_k$への射影と考えれば良い：
$$ X\boldsymbol{v}_k = \sigma_k \boldsymbol{u}_k $$
特異値分解にはNumpyの関数numpy.linalg.svdを用いる。
[使い方]
U, s, Vh = np.linalg.svd(a)a: （入力）行列U: （出力）直交行列。$k$番目の左特異ベクトルは$U[:,k]$s: （出力）特異値。降順にソートされて出力。Vh:（出力）直交行列。$k$番目の右特異ベクトルは$Vh[k,:]$
- 参考: https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.htmlIn [3]:
# 特異値分解(Numpyの関数を使う)
U, s, Vh = np.linalg.svd(anom)
# 描画
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
# -------（左パネル）------------
## x,y軸をそれぞれ国語、算数の得点として、(1) それぞれの得点をプロット ＋ (2) 固有ベクトルを重ね書き
## (1) オリジナルデータを描画
ax.scatter(anom[:,0],anom[:,1])
## (2) 主成分ベクトル(EOF1, EOF2)を描画
colors = ["red","blue"]
labels = ["EOF1","EOF2"]
for i in range(2):
print("第",i,"モード")
print("特異値 = ", s[i])
print("主成分ベクトル(EOF) = ",Vh[i,:])
print("----")
# 主成分ベクトル．ベクトルの大きさは，主成分得点の標準偏差(std)とする（Vhは大きさ１に規格化されていることに注意）
# std = s/m**0.5 (s:特異値, m: データ数)
x2 = s[i]*Vh[i,0]/m**0.5; y2 = s[i]*Vh[i,1]/m**0.5
ax.quiver([0],[0],[x2],[y2],color=colors[i],label ="EOF"+str(i),angles='xy',scale_units='xy',scale=1)
ax.legend()
ax.set_xlabel("Japanese"); ax.set_ylabel("Math")
ax.set_xlim(-20,20); ax.set_ylim(-20,20)
# -----（左パネル、ここまで）-----
# -------(ここから 右パネル）-------
## x,y軸をそれぞれEOF1, EOF2として（基底を取り直して）、それぞれの主成分得点をプロット
pc1 = s[0]*U[:,0]
pc2 = s[1]*U[:,1]
#元の偏差行列の主成分ベクトルへの射影を計算しても良い
#V = Vh.T#V(=Vhの転置)
#pc1 = anom @ V[:,0]# "@" は行列の積
#pc2 = anom @ V[:,1]
ax2.scatter(pc1,pc2)
ax2.set_xlabel('PC1'); ax2.set_ylabel('PC2')
ax2.set_xlim(-20,20); ax2.set_ylim(-20,20)
# ------（右パネル、ここまで）------
plt.show()
（課題 1-1, -2）三科目データを対象にEOF解析を施して結果を描画せよ。¶
三次元領域でプロットするには、サブプロット領域の作成時に3Dオプションをつけるだけで良い：
ax = fig.add_subplot(121,projection='3d')
%matplotlib notebookと書いておけば（インターラクティブなグラフになり）マウスで視点を動かすことができる。
In [4]:
%matplotlib notebook
data = X_org.values # .valuesで数値を取り出す（ラベルは剥がす）
mean = data.mean(axis=0) # 平均
anom = data - mean # 平均偏差行列
m,n = anom.shape
# 特異値分解(Numpyの関数を使う)
U, s, Vh = np.linalg.svd(anom)
# 描画
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121,projection='3d')
ax2 = fig.add_subplot(122,projection='3d')
# -------（左パネル）------------
## x,y軸をそれぞれ国語、算数の得点として、(1) それぞれの得点をプロット ＋ (2) 固有ベクトルを重ね書き
## (1) オリジナルデータを描画
ax.scatter(anom[:,0],anom[:,1],anom[:,2])
## (2) 主成分ベクトル(EOF1, EOF2)を描画
colors = ["red","blue","green"]
labels = ["EOF1","EOF2","EOF3"]
for i in range(3):
print("第",i,"モード")
print("特異値 = ", s[i])
print("主成分ベクトル(EOF) = ",Vh[i,:])
print("----")
# 主成分ベクトル．ベクトルの大きさは，主成分得点の標準偏差とする
x2 = s[i]*Vh[i,0]/m**0.5; y2 = s[i]*Vh[i,1]/m**0.5; z2 = s[i]*Vh[i,2]/m**0.5
print('EOF',i,' = ', [x2, y2, z2])
ax.quiver([0],[0],[0],[x2],[y2],[z2],color=colors[i],label ="EOF"+str(i))
# 3Dプロットでは、angles や scale_unitなどのオプションは無いようです
ax.legend()
ax.set_xlabel("Japanese"); ax.set_ylabel("Math"); ax.set_zlabel("English")
ax.set_xlim(-20,20); ax.set_ylim(-20,20); ax.set_zlim(-20,20)
# -----（左パネル、ここまで）-----
# -------(ここから 右パネル）-------
## x,y軸をそれぞれEOF1, EOF2として（基底を取り直して）、それぞれの主成分得点をプロット
pc1 = s[0]*U[:,0]
pc2 = s[1]*U[:,1]
pc3 = s[2]*U[:,2]
ax2.scatter(pc1,pc2,pc3)
ax2.set_xlabel('PC1'); ax2.set_ylabel('PC2'); ax2.set_zlabel('PC3')
ax2.set_xlim(-20,20); ax2.set_ylim(-20,20); ax2.set_zlim(-20,20)
# ------（右パネル、ここまで）------
plt.show()
(B) 共分散行列から求める方法¶
(1) 共分散行列を作る: $$ \begin{align} S_x = \frac{1}{m}\tilde X^T \tilde X \end{align} $$
(2) 固有値問題を解くことで、主成分ベクトル($\boldsymbol v_k$)を求める:
$$ \begin{align} S_x \boldsymbol v_k = \lambda \boldsymbol v_k \end{align} $$(3) 偏差行列を$v_k$に射影することで$\sigma_k \boldsymbol u_k$を求める。
（課題 1-3）上記の手順により、（A）と同様の結果が得られることを確認せよ。¶
なお、固有値分解にはNumpyの関数numpy.linalg.eigを用いる。
[使い方]
w, v = np.linalg.eig(a)a: （入力）行列w: （出力）固有値。ソートされていない。v:（出力）直交行列。固有値$w[k]$に対応する固有ベクトルは$v[:,k]$
- 参考: https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.htmlIn [5]:
data = X_org.values# .valuesで数値を取り出す（ラベルは剥がす）
mean = data.mean(axis=0) # 平均
anom = data - mean # 平均偏差行列
m,n = anom.shape
#print(m,n)
# 共分散行列
Cov = 1/m*np.dot(anom.T,anom)
print(Cov)
# 固有値を求める
w, v = np.linalg.eig(Cov)
#固有値iのi番目の固有ベクトルはv[:,i]
print("Singular values", (w*m)**0.5)
print("Eigenvectors:", v[:,0],v[:,1],v[:,2])
In [6]:
var = np.sum( anom.var(axis=0) )
# 分散が固有値の和で表せることを確かめる
print("Variance in original coorinate: ", var)
print("Variance in EOF1-2 coordinate: ", np.sum(w))
l = w.shape
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121,projection='3d')
ax2 = fig.add_subplot(122,projection='3d')
# 元データ
ax.scatter(anom[:,0],anom[:,1],anom[:,2])
# 新しい軸
c = ["red","blue","green"]
labels = ["EOF1","EOF2","EOF3"]
pc = []
for i in range(len(w)):
idx = np.argsort(w)[::-1][i]
#降順（固有値が大きい順）に並べ直したときのi番目のインデックス
#（つまり、固有値が最大値を取るインデックスを取得）
print("index", idx)
ll = w[idx]; print(w[idx])
x2 = ll**0.5*v[0,idx]; y2 = ll**0.5*v[1,idx]; z2 = ll**0.5*v[2,idx]
# 偏差行列を特異値分解して得られる特異値(s)と,共分散行列を固有値分解して得られる固有値（ll)には
# s = sqrt(ll)*sqrt(m)という関係がありました。
# ここではベクトルの長さは s/sqrt(m) = sqrt(ll)としています（標準偏差）
ax.quiver([0],[0],[0],[x2],[y2],[z2], color=colors[i],label ="EOF"+str(i))
#偏差行列の、新しい基底ベクトル(v)への射影を計算。
pci = anom @ v[:,idx]
pc.append(pci)
ax.set_title('Score of Japanese & Math')
ax.set_xlabel("Japanese"); ax.set_ylabel("Math")
ax.set_xlim(-20,20); ax.set_ylim(-20,20)
ax2.scatter(pc[0],pc[1],pc[2])
ax2.set_xlabel('PC1'); ax2.set_ylabel('PC2'); ax2.set_zlabel('PC3')
ax2.set_xlim(-20,20); ax2.set_ylim(-20,20); ax2.set_zlim(-20,20)
ax.legend() #凡例をかく
plt.show()
In [7]:
# データの読み込み
file = "data/Maizuru_1hr/Maizuru_2017_uv.csv"
# -->適宜変更
X_org = pd.read_csv(file,encoding='shift_jis',index_col=[0],parse_dates=[0])
print(X_org.head())
# 偏差行列にする（データは三列あるが、まずは簡単のため最初の二列の数値だけ取り出す)
columns = X_org.columns
data = X_org.values[:,0:2] # .valuesで数値を取り出す（ラベルは剥がす）
In [8]:
mean = data.mean(axis=0) # 平均
anom = data - mean # 平均偏差行列
m,n = anom.shape
In [9]:
# 特異値分解(Numpyの関数を使う)
U, s, Vh = np.linalg.svd(anom)
# 描画
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
# -------（左パネル）------------
## x,y軸をそれぞれ東西風速(U)、南北風速(V)として、(1) 散布図をプロット ＋ (2) 固有ベクトルを重ね書き
## (1) オリジナルデータを描画
ax.scatter(anom[:,0],anom[:,1])
## (2) 主成分ベクトル(EOF1, EOF2)を描画
colors = ["red","blue"]
labels = ["EOF1","EOF2"]
for i in range(2):
# 主成分ベクトル．ベクトルの大きさは，主成分得点の標準偏差とする
x2 = s[i]*Vh[i,0]/m**0.5; y2 = s[i]*Vh[i,1]/m**0.5
print('EOF',i,' = ', [x2, y2])
ax.quiver([0],[0],[x2],[y2],color=colors[i],label ="EOF"+str(i),angles='xy',scale_units='xy',scale=1)
ax.legend()
ax.set_xlabel("U-wind"); ax.set_ylabel("V-wind")
ax.set_xlim(-20,20); ax.set_ylim(-20,20)
# -----（左パネル、ここまで）-----
# -------(ここから 右パネル）-------
## x,y軸をそれぞれEOF1, EOF2として（基底を取り直して）、それぞれの主成分得点をプロット
pc1 = s[0]*U[:,0]
pc2 = s[1]*U[:,1]
ax2.scatter(pc1,pc2,s=3)
ax2.set_xlabel('PC1'); ax2.set_ylabel('PC2')
ax2.set_xlim(-20,20); ax2.set_ylim(-20,20)
# ------（右パネル、ここまで）------
plt.show()
In [ ]: