気候データマネゴト解析：実行編

この記事の続きになります！

記事の概要

自然豊かな北海道。この記事は北海道の夏季（6,7,8月）の気温の変動傾向がどんな要因によって支配されているかを主成分分析を利用して確認してみます。（その実行編！）

利用データとツール

DL_Link：気象庁過去のデータダウンロード
場所：気象庁の北海道道内22地点の観測所
項目：気温（℃）
年数：1970~2020の51年間の6,7,8月
地点：北海道内22地点
解析ツール：Jupyter notebook（Python）

前回は気象庁からDLしたデータの整形を行い、行成分に年代情報、列成分に地点情報を含むデータフレームXを作成するところまで行いました。

今回は主成分分析を行い、その結果を検討まで行えればなと思います。

#念の為...以下で使用するモジュールたちです
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

まずはXを確認します。

#Xの先頭5行を表示
X.head()

#Xの行列数を確認
print(X.shape)

#出力(Xは51行22列の行列)
#(51, 22)

主成分分析の実行

ここから主成分分析を実行していきます。

#sklearnから標準化用のモジュールをインポート
from sklearn.preprocessing import StandardScaler

#クラスを作成sc = StandardScaler()

#データフレームXを標準化
X_std = sc.fit_transform(X)

以上でX_stdという変数に標準化処理を行ったXを代入します。標準化はX全体の平均を0、標準偏差を1に統一する処理になります。
※この後に、Xの分散共分散行列を作成するのでカラム間のばらつきを打ち消すために無次元にする必要があるためです。

それでは次に標準化したX（X_std）の分散共分散行列を求めます。
分散共分散行列についてはこちら

#np.cov()は共分散を求める関数になります。
#変数cov_matはX_stdの分散共分散行列。
cov_mat = np.cov(X_std.T)

#eigen vals, eigen vecsにそれぞれ、固有値、固有ベクトルを代入にします。
#linalg.eig()関数は固有値と固有ベクトルを求めます。
eigen_vals, eigen_vecs=np.linalg.eig(cov_mat)

#固有値を出力してみます。
#%記法ですね。%sの部分に、engen_valsを代入する形で出力する書き方です。
print('\nEigenvalues \n%s' % eigen_vals)

##出力##
Eigenvalues 
[2.04901297e+01 6.04838941e-01 4.97032553e-01 2.37624774e-01
1.41416078e-01 1.02522916e-01 7.20535035e-02 6.38128191e-02
5.19063694e-02 4.04408288e-02 2.73766628e-02 2.06771175e-02
1.96468479e-02 1.66726837e-02 1.33778032e-02 2.14823263e-03
9.39770554e-03 8.36786308e-03 6.88359836e-03 5.77602124e-03
3.60561844e-03 4.29136700e-03]

np.cov()は列ベクトルを計算するので、X_stdを.Tで転置したものを計算します。（これによって、列ベクトルが年代ごとのデータから、地点ごとのデータになります。今回は地点ごとの気候の違いに着目するので意外と注意が必要な処理だと思います。）

linalg.eig()関数についてはこちら

ここまでで、固有値と固有ベクトルが計算できてしまいました、、関数サマサマですほんとに。

分散説明率の可視化

次に分散説明率を計算しましょう。これで説明に使う主成分の数に当たりをつけることができます。

#固有値の分散説明率をプロット
tot = sum(eigen_vals)

#分散説明率の説明var_exp = [(i/tot) for i in sorted (eigen_vals, reverse=True)]

#分散説明率の累積和の取得
cum_var_exp = np.cumsum(var_exp)

import matplotlib.pyplot as plt

#分散説明率の棒グラフを作成
plt.bar(range(1,23), var_exp, alpha=0.5, align="center", label="Individual explained variance")

#分散説明率の累積和の階段グラフを作成
plt.step(range(1,23), cum_var_exp, where="mid", label="Cumulative explained variance")

plt.ylabel("Explained variance ratio")
plt.xlabel("principal component index")

plt.legend(loc="best")
plt.tight_layout()
plt.show()

結果は第一主成分で分散の90%、第二主成分までで95%が説明可能であることがわかりました。よって解釈に利用するデータは二次元で十分そうです。

次元削減

それでは説明率の高い２つの固有値固有ベクトルを抽出する処理を行い、特徴量変換（新しい軸）を作成します。

#(固有値, 固有ベクトル)のタプルのリストをリスト内包表記で作成
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(key=lambda k: k[0], reverse=True)

eigen_pairsに固有値, 固有ベクトルを作成します。
リスト内包表記でのリスト作成はこちらが参考になります。

その後、sort関数を用いて、固有値の降順で並び替えます。
sort関数のoptionはk : k[0]で一番目の要素をkey選択させてreverseで降順を指定しています。

#出力してみる
eigen_pairs

固有値と固有ベクトルがセットになっていることがわかります。22次元（地点）分存在します。

今回は分散説明率の高い上位2成分の固有ベクトルをここから抽出して新たな2次元配列を作成します。

#特徴量変換。新しい主成分軸に変換する#固有ベクトルの配列の作成
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)

#特徴量変換（22次元のデータに2列の行列を掛け算し、分散を最大化する条件で2次元のデータに変換）
X_pca = X_std.dot(w)

hstack関数を用いて、作成した新たな配列を結合します。
ここでwは第一、二主成分の固有ベクトルになります。

.dot()関数でX_stdとwのドット積（行列の積）を取ることで、特徴量変換を行います。

固有ベクトルの可視化

ここから、地点ごとの第一、第二主成分の固有ベクトルを可視化して、北海道の気候について考えてみます。

#固有ベクトルを分割してそれぞれのリストを作成
pc1=[]
pc2=[]
for i in range(len(w)):
   pc1.append(w[i][0])
   pc2.append(w[i][1])

%matplotlib inline
left = np.arange(len(w))  # numpyで横軸を設定
labels = ['函館','江差','倶知安','寿都','小樽','札幌','苫小牧','室蘭','岩見沢','羽幌','留萌','稚内','北見枝幸','旭川','浦河','雄武','紋別','網走','帯広','広尾','根室','釧路']
height = 0.3
fig = plt.figure(figsize=(12,9))

#第一主成分の固有ベクトルを出力するときはこちらを実行
#plt.barh(left, pc1, color='r', height=height, align='center')

#第二主成分の固有ベクトルを出力するときはこちらを実行
#plt.barh(left+height, pc2, color='b', height=height, align='center')

plt.yticks(left + height/2, labels)
plt.show()

出力されたグラフは以下のようになります。

第一主成分の固有ベクトルは全ての地点で負の値となり、値の変動も小さいです。各地の気温の平均値と第一主成分の経年変化はほぼ同じ変動を示します（図略）。

よって主成分は各地の気温の経年変動を示していると解釈できるでしょう。

一方、第二主成分の固有ベクトルは正負も異なり、変動が大きくなっています。正負を地点ごとにプロットすると以下のようになります。

道東部・道南部は同西部と正負が明確に異なることがわかりました。
よって第二主成分は北海道内の気候の地域差を示す指標と解釈することができそうです。

（ホントは正負の地域ごとに気温の平均値をとって差が有意であることを示したいのですが長くなるので社会人になって時間があるときに取り組みたいと思います。）

まとめ

北海道の約50年間の気温データに対して主成分分析を行いました。

第一主成分は経年変動
第二主成分は地域特性
を示すことが考えられそうです！