こーひーをぶんなぐれ！ - CSVデータを加工して主成分分析

本項でもGoogle Colaboratoryを使用する．（もちろんファイル読み込みの部分以外はJupyter Notebookでもよいし，一気にPythonのソースコードをまとめて書いて実行しても問題はない）

主成分分析

これまで「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」という4つの観測値を扱ってきた．この「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」という問題を考える．

ところが，観測値が4つある時，つまり説明変数が4つある時，2次元グラフはそれぞれの軸の組み合わせて6通り存在する．さらに3つ以上の観測値を軸にして見ることはできないため，グラフでぱっと見ではわからない．

そこで「主成分分析」という「次元削減」手法を用いる．

主成分分析は，説明変数が多数存在する時に，元のデータの特徴を十分に表しつつ，かつ次元が少ない説明変数のセットを「合成」する手法である．

実際のには，上図のようにただ足し合わせるだけではなく，変換行列を作成し，それによって次元削減をするが，その変換行列を作成する際に「元のデータの特徴を十分に表現できる」ことを目標として調整をしていく．

これは3DCGのレンダリングの時の3次元座標から2次元座標の変換行列（カメラの位置，アングル，拡大度等から算出される）を「3Dモデルのディティールが上手く掴める」ように調整するのに，イメージ的には近い．

十分に表現できるような「カメラの位置」「カメラの角度」「画角」「拡大」などをうまく調整し，3Dデータの特徴を十分に表現できるような2次元画像にするのと，主成分分析のための変換行列の作成は，似たようなもんである．

使用するデータと加工

今回はデータとして，相関係数を出した時と同じ「東日本の県庁所在地」の「1977年〜2018年」の「8月の平均気温」「8月の日最高気温の月平均」「8月の日最低気温の月平均」「25度を超えた日数」のCSVファイルを使う．

さらにそこから，盛岡，甲府，横浜に絞ってスタック形式のデータを作る．

GoogleスプレッドシートにCSVを読み込む．

例によって不要な行と列を削除．

左下の「＋」マークを押して，新規シートを追加．

シートの名前をわかりやすくつける．

最初のシートから盛岡の全データをコピーして，

新作ったシートに貼り付ける．（一番左の列は地点ラベルが入るので，2番目の列から貼り付ける）

盛岡のデータが貼り付けられた状態．

一番左の列名を「観測地点」にして，

1つのセルに「盛岡」と書き込んでコピーし，

データの他の行に貼り付ける

盛岡のデータは完成．

同じように甲府のデータもコピペして地点ラベルをつける．

横浜のデータも同様にする．

列名を，Rでも読み込める短いものにする．

これまでと同じく「ファイル」→「ウェブに公開」

スタック形式データのシートを選択し，CSVとして公開．公開するシートの選択を間違わないように．

import numpy as np
import pandas as pd

df = pd.read_csv('公開したURL')
df.head()

そのままグラフを描画してみる

まず，説明変数4つをx軸とy軸に取り，2次元のグラフを描画している．説明変数の組み合わせは合計6通りあるが，とりあえず3つの組み合わせでグラフを描画する．

# 生データを散布図でプロットしてみる
# 説明変数が4つで2次元に表示するためには，6通りのxとyの組み合わせがあるが，とりあえず3つだけ
import matplotlib.pyplot as plt

# 都市名ラベルの0番目を除いた説明変数4つ分で3つのグラフを描く
for i in range(1, 4, 1):
  for j, a_city in enumerate(['Morioka', 'Kofu', 'Yokohama']):
    # 42年分 * 3都市分なので，都市ごとに色を変えるため，データを分けてscatterを呼び出す
    x = df.iloc[j*42:(j+1)*42, i].astype(np.float32)
    y = df.iloc[j*42:(j+1)*42, i+1].astype(np.float32)
    plt.scatter(x, y, label=a_city)
  if i==1:
    plt.xlabel('Ave. of Temp in Aug.')
    plt.ylabel('Ave. of Max Temp. in Aug')
  elif i==2:
    plt.xlabel('Ave. of Max Temp. in Aug')
    plt.ylabel('Ave. of Min Temp. in Aug')
  else:
    plt.xlabel('Ave. of Min Temp. in Aug')
    plt.ylabel('Count of over 25 deg in Aug')
  plt.legend()
  plt.show()

X軸が「8月の平均気温」，Y軸が「8月の日最高気温の平均」

X軸が「8月の日最高気温の平均」，Y軸が「8月の日最低気温の平均」

X軸が「8月の日最低気温の平均」，Y軸が「25度を超えた日数」

今は都市ごとに色分けをしているため判別しやすいが，横浜はどのグラフでも分離して見える（線形判別分析が簡単にできそう，すなわち判別の「直線」を引くことが可能に見える，「線形分離可能」と呼ぶ．）のに対して，盛岡と甲府はかなり被っているように見える．（盛岡と甲府は同じ直線上に載っているように見える）

主成分分析は「このグラフが分離して見えるような表示になるような「合成変数」を作成すること」が目的である．合成変数は任意の数（次元）で作成することができる．合成変数を2つ（2次元）にしてしまえば，2次元グラフでも簡単に表示できるようになるはず．

主成分分析の実行

scikit-learnの中に入っているPCAクラスが使いやすい．主成分分析は「説明変数に対して行う」ものなので，PCAクラスでfitする時に食わせるデータは「説明変数のみ」である．

説明変数間で単位が違う場合「標準化」or「正規化」という作業をしなければならない（値のスケールをうまい具合に整えてくれるイメージ）．これは同じくscikit-learnの中に含まれているStandardScalerクラスで簡単に行うことができる．標準化するのは同じく説明変数のみであることに注意．

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 単位が違う説明変数（3つの説明変数の単位は摂氏温度で，1つの単位は日数である）を「標準化」or「正規化」する
scaler = StandardScaler()
df_sc = scaler.fit_transform(df.loc[:, ['ave', 'max.ave', 'min.ave', '25over']])

# 合成変数を2次元で作ると指定
pca_model = PCA(n_components=2)
# データを学習して合成変数への変換行列を作る
pca_model.fit(df_sc)

# 変換行列ができたので，
# 4次元(8月の平均気温，8月の日最高気温の平均，8月の日最低気温の平均，25度を超えた日数）の説明変数を，
# 2次元の合成変数に変換する（合成する）
transformed_variables = pca_model.transform(df.loc[:,  ['ave', 'max.ave', 'min.ave', '25over']])
# 使いやすいようにデータフレーム化しておく
transformed_df = pd.DataFrame(transformed_variables)

# 合成変数2つの散布図を描く
for i, a_city in enumerate( ['Morioka', 'Kofu', 'Yokohama']):
  plt.scatter(transformed_df.iloc[i*42:(i+1)*42, 0], transformed_df.iloc[i*42:(i+1)*42, 1], label=a_city)
plt.legend()
plt.show()

結果は以下のようになる．

このグラフでは，X軸が主成分つまり合成変数0(Principal Component)0 , Y軸が主成分つまり合成変数1である．

合成変数では，盛岡，甲府，横浜が綺麗に分離しており，線形判別分析が容易になっている．この3都市の8月の気温に関しては，主成分分析 - 合成変数が有効であることを示している．

この状態で最初の目的である「4つの観測値から『都市の気候』をうまく表現したい．うまく説明できるようにしたい．」を達成できたかを考えると，「直線の補助線を引くことで，「この線のこちら側が盛岡で，この線のこちら側が甲府で，この線のこちら側が横浜です」とうように明確に説明できるようになる」ので，目的を達成したと言える．

寄与率

主成分分析でもう一つ確認しなければならないのが「寄与率」である．主成分分析をした後に，モデルを保持している変数のexplained_variance_ratio_という変数の中に寄与率が入っている．

# 寄与率の表示
print(pca_model.explained_variance_ratio_)

寄与率の出力結果は以下の通り．

[0.97995861 0.01258219]

寄与率というのは，それぞれの主成分つまり合成変数が「データをどれだけ表現できているか」を表したものである．寄与率の合計が1.0に十分に近ければ，合成変数のみでデータを十分表現できていると言っても良い．

寄与率の出力は配列になっており，0番目が主成分合成変数0の寄与率，1番目が主成分合成変数1の寄与率である．合計が1に十分に近いので，この2次元だけで元データを十分に表現していると言って良い．つまり「4次元あった変数は2次元に落とし込むことが可能」ということである．

逆に合計寄与率が低かった場合，「この次元数では元のデータを表現できているとは言い難い」となる．この場合，合成変数の次元数を増やす必要がある．1次元ずつ増やしていって，寄与率が1に十分近くなる次元数を探す．次元数を増やしても合計寄与率が上がらなくなった場合はそこで止める．

主成分分析の合成変数間の相関係数

# 合成変数間の相関係数の表示
print(transformed_df.corr())

「合成変数間の相関係数」は，普通にデータフレームの内部の相関係数を出力すれば良い．

              0             1
0  1.000000e+00 -1.182588e-15
1 -1.182588e-15  1.000000e+00

合成変数間の相関係数は非常に小さい（ほぼ0）．というのも，主成分分析の内部計算は「合成変数間の相関係数を低くする」ように調整するからである．逆に言えば「相関が高いものをまとめる」という調整とも言える．