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 [6]:
# データの読み込み
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 [11]:
# 特異値分解(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と書いておけば（インターラクティブなグラフになり）マウスで視点を動かすことができる。
(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 [4]:
# データの読み込み
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 [ ]: