🙆‍♀️

主成分分析(PCA)の実装_世界の統計

2024/09/02に公開

前書き

機械学習に関する仕事をしながら年に複数回旅行ができれば幸せな人生だと考えているAIエンジニアです。今年は初めて行く都市が多くてなかなか良い年です。人生の幕を閉じる前に行きたいところリストを制覇したいものです。仕事と旅行は計画的に!


今までの記事で
・主成分分析の哲学(考え方)と理論
・データの作成
をまとめました。

今回は実装部分をまとめていきたいと思います。

主成分分析の実装

処理の大きな流れ

主成分分析の理論の結果、分散共分散行列の固有ベクトルと固有値を利用して軸を新しく作ることで特徴的な部分を残しながら情報を削減できる事を学びました。
よってやるべき事をまとめると
・元のデータから分散共分散行列の作成
・分散共分散行列の固有ベクトル・固有値の計算
・固有ベクトルと固有値を利用した元データの変換
となります。

大きく分類するとこれだけですが、細かな部分で注意すべき点もあるので都度書いていきます。

主成分分析の考え方と理論はこちらのページにまとめてあります。
https://zenn.dev/goals_techblog/articles/1b031d8b144ec5

使用するデータ

まずは使用するデータについてです。
総務省統計局が提供してくださっている世界の統計のデータを利用しています。
出典:「世界の統計2024」(総務省統計局) (URL) (2024年7月25日に利用)

https://www.stat.go.jp/data/sekai/0116.html

テーブルデータの作成はこちらの記事にまとめてあります。
よろしければご覧ください。

https://zenn.dev/goals_techblog/articles/178231b1aa829c

データの標準化(Z-score normalization)

まず今回"country"のカラムは使用しないので削除します。またfloat型に変換しておきます。

import os

import numpy as np
import pandas as pd

dir_path = "workspace_dir/input_data"
file_name = "20240725_for_pca.csv"
data_path = os.path.join(dir_path, file_name)
using_df = pd.read_csv(data_path, index_col=0)

pca_df = using_df.iloc[:, 1:]
pca_df = pca_df.astype("float64")
pca_df
land_area forest_area forest_rate artificial_forest_area forest_area_diff forest_area_diff_rate num estimated_num_2020 estimated_num_2022 birth_rate_2015 birth_rate_2020 death_rate_2015 death_rate_2020 gdp_2015 gdp_2017 gdp_2018 gdp_2019 gdp_2020 gdp_2021 life_expectancy male_life_expectancy female_life_expectancy healthy_life_expectancy male_healthy_life_expectancy female_healthy_life_expectancy population_density
0 364560 249350 68.4 101840 -3.1 -0 126146 1.26146e+08 124947 8 6.6 11 12.1 4.4454e+06 4.9313e+06 5.0425e+06 5.1179e+06 5.0483e+06 5.0037e+06 84 81 87 74 73 75 346.023
1 1.62876e+06 107520 6.6 10010 6 0.1 79926 8.729e+07 88551 19.4 14.2 4.8 5.6 417210 503710 526365 519356 543654 594892 77 76 79 66 66 67 53.5929
2 2.97319e+06 721600 24.3 132690 266.4 0.4 1.21086e+06 1.39639e+09 1.41717e+06 18.8 16.6 6.7 7.4 2.14676e+06 2.62433e+06 2.7632e+06 2.85073e+06 2.6722e+06 3.20147e+06 71 70 72 60 60 60 469.66
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
31 263310 98930 37.6 20840 4.4 0.1 4793 5.061e+06 5185 13.1 12.5 6.6 6.4 178064 206624 211953 213435 212214 250451 82 80 84 70 70 71 19.2207

当たり前ですが、作成したデータは列毎に様々なデータが入っていて単位もスケール(値の大きさやバラツキの大きさ)もバラバラな状態です。
このような時、それぞれの値を列毎の平均が0, 標準偏差が1になるような値(Z値)に変換することで単位やスケールの違いを気にせずに済むようになります。
※この変換はZ-score normalizationとか標準化(Standardization)とかよばれます。
(広い意味でのデータの正規化です。)

数式で書くと

z_i = \frac{x_i - \bar{x}} {s_x}

このような計算をします。

以下のコードでZ値に変換後の行列を作成できます。

def generate_z_score_matrix(matrix):
    """Z値の行列を作成する"""
    column_mean = np.mean(matrix, axis=0)
    print("mean")
    print(column_mean)
    column_std = np.std(matrix, axis=0, ddof=1)
    print()
    print("std")
    print(column_std)
    z_score_matrix = ( (matrix - column_mean) /  column_std)
    print()
    print("Z_mean")
    print(np.mean(Z, axis=0))
    print()
    print("Z_std")
    print(np.std(Z, axis=0, ddof=1))
    return z_score_matrix

z_score_matrix = generate_z_score_matrix(pca_df.values)

元の平均と標準偏差はバラバラなのに対し、z値に変換後の平均(Z_mean)と標準偏差(Z_std)は0と1になっていることが確認できます。
(平均(Z_mean)は数値計算の誤差の関係で0ちょうどにはなっていませんが、10のマイナス17乗とかの値になっているのでほとんど0です。ちっちゃいことは気にするな♪)

分散共分散行列の作成

次に分散共分散行列を計算していきます。
列毎の平均(Z値に変換した場合は0)を引いた行列の転置行列と転置前の行列の積を計算しデータ数で割ることで計算できます。

def generate_cov_matrix(matrix):
    """
    分散共分散行列を作成する
    行列の積ではデータ数(不偏標準偏差を利用した場合はデータ数 - 1)倍された値になるのでそれぞれの要素をデータ数で割る
    """
    data_num = matrix.shape[0] - 1
    means = np.mean(matrix, axis=0)
    diff_mean = matrix - means
    cov_sum = (diff_mean.T @ diff_mean)
    cov_matrix = cov_sum / data_num
    return cov_matrix

cov_matrix = generate_cov_matrix(z_score_matrix)
pd.DataFrame(cov_matrix)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
0 1 0.932203 0.0757723 0.569541 0.14188 0.0911632 0.32008 0.296979 0.293587 -0.148606 -0.186378 0.147459 0.187319 0.444149 0.449921 0.449978 0.45003 0.44336 0.456989 0.0284587 -0.0265941 0.0591371 0.00362972 -0.0394569 0.0555815 -0.314272
1 0.932203 1 0.27893 0.404345 -0.117657 0.010283 0.165383 0.148769 0.146134 -0.145854 -0.173204 0.234691 0.289318 0.292974 0.297275 0.292289 0.292896 0.280096 0.286114 -0.00654636 -0.0831449 0.041622 -0.0234347 -0.0788388 0.0469483 -0.277443
2 0.0757723 0.27893 1 0.058358 -0.32963 0.0445849 -0.13779 -0.144617 -0.14722 -0.359105 -0.335283 0.199773 0.270019 0.0461734 0.0476542 0.0404982 0.038841 0.0304043 0.0156736 0.346861 0.288351 0.386606 0.389166 0.340133 0.421597 -0.093562
3 0.569541 0.404345 0.058358 1 0.644057 0.386283 0.762824 0.717163 0.711146 -0.252727 -0.288894 0.0272426 0.0510501 0.680075 0.69382 0.715067 0.712884 0.724514 0.75398 0.128878 0.121591 0.139978 0.167382 0.139716 0.166551 -0.0955482
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
25 -0.314272 -0.277443 -0.093562 -0.0955482 0.110673 0.135335 0.175742 0.179497 0.181717 -0.00831322 0.00106182 -0.162259 -0.216104 -0.0791432 -0.0737726 -0.0705042 -0.0695425 -0.0648732 -0.0681703 -0.0735503 -0.0411716 -0.0885003 -0.0588345 -0.0366955 -0.0907202 1

実対象行列(要素が全て実数で転置しても元の行列と一致する行列)になっていることを確認できます。
ちなみにZ値に変換したデータの分散共分散行列は相関係数行列と一致します。
(相関係数は共分散をそれぞれの項目の標準偏差で割って正規化したものですが、Z値に変換すると標準偏差が1になり分母が1になるためです。)


参考までにnpの関数との差です。

diff_frobenius_norm_np_cov = ((cov_matrix - np.cov(z_score_matrix, y=None, rowvar=False)) ** 2).sum()
print("npの分散共分散行列との差: ", diff_frobenius_norm_np_cov)
# npの分散共分散行列との差:  8.640941423459553e-31

diff_frobenius_norm_np_corrcoef = ((cov_matrix - np.corrcoef(z_score_matrix, y=None, rowvar=False)) ** 2).sum()
print("npの相関係数行列との差: ", diff_frobenius_norm_np_corrcoef)
# npの相関係数行列との差:  5.1418643223629736e-30

↑差がほとんど0になっていることが確認できます。


固有ベクトル・固有値の計算

次に固有ベクトル・固有値を計算していきます。
npのlinalg.eigを使うことで計算できます。

e_values, e_vector_matrix = np.linalg.eig(cov_matrix)
print("固有値")
print(e_values.shape)
# (26,)

print("固有ベクトル")
print(e_vector_matrix.shape)
# (26, 26)

結果の確認

固有値と寄与率の結果

具体的に中身を見ていきます。

先に固有値を見ていきます。

分散共分散行列の固有値は分散を意味します。
全部の分散のうち合成後のそれぞれの軸が占めている分散の割合(寄与率)とその累積値(累積寄与率)を計算してまとめてみます。

def generate_evalue_info_df(e_values):
    """"固有値, 寄与率, 累積寄与率をまとめたdfを作成"""
    e_value_sum = e_values.sum()
    # 寄与率の計算
    contribution_ratio = e_values / e_value_sum
    # 累積寄与率の計算
    cumulative_contribution_ratio = np.cumsum(contribution_ratio)
    # 固有値、寄与率、累積寄与率をまとめる
    e_value_info = np.concatenate([[e_values], [contribution_ratio], [cumulative_contribution_ratio]], 0)

    index = ["固有値", "寄与率", "累積寄与率"]
    columns = ["p" + str(index+1) for index in range(len(e_values))]
    e_value_info_df = pd.DataFrame(e_value_info, index=index, columns=columns)
    return e_value_info_df

e_value_info_df = generate_evalue_info_df(e_values)
e_value_info_df
p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 p13 p14 p15 p16 p17 p18 p19 p20 p21 p22 p23 p24 p25 p26
固有値 9.85899 6.96661 2.98104 1.73924 1.42611 1.01462 0.736971 0.4863 0.370923 0.22565 0.0904013 0.0439865 0.0279827 0.013453 0.00783747 0.00434769 0.00225093 0.00155896 0.000803999 0.000641084 0.000158056 6.85892e-05 4.01978e-05 1.12018e-05 1.0105e-06 5.46367e-06
寄与率 0.379192 0.267947 0.114655 0.0668937 0.0548506 0.0390238 0.028345 0.0187039 0.0142663 0.00867883 0.00347697 0.00169179 0.00107626 0.000517423 0.000301441 0.000167219 8.65741e-05 5.99601e-05 3.0923e-05 2.46571e-05 6.07908e-06 2.63804e-06 1.54607e-06 4.30838e-07 3.88656e-08 2.10141e-07
累積寄与率 0.379192 0.647138 0.761794 0.828687 0.883538 0.922562 0.950907 0.969611 0.983877 0.992556 0.996033 0.997725 0.998801 0.999318 0.99962 0.999787 0.999874 0.999933 0.999964 0.999989 0.999995 0.999998 0.999999 1 1 1

↑累積寄与率の結果から第一主成分(p1)だけで全体の分散の約38%, 第一主成分(p1)と第二主成分(p2)の2軸で約65%、同様に第三主成分(p3)までの3軸で約76%を表現できることが分かります。

また寄与率の結果から第四主成分(p4)以降は全体の分散の約10%以下であることがわかります。

累積寄与率の結果を折れ線グラフで見てみます。

e_value_info_df.T.plot.line(y=["累積寄与率"], ylim=(0, 1), marker="o")

↑グラフからも第一主成分(p1)で40%弱、第一主成分(p1)と第二主成分(p2)の2軸で60%強の分散を表現できていることがわかります。また、第十一主成分(p11)以降の軸はほとんど影響しないことも分かります。

次元を削減する時の基準は様々ありますが、今回は全体の70%以上を表現できているという基準で次元を3つに削減することにします。


固有ベクトルの結果と軸の解釈

次に固有ベクトルの結果を見ていきます。
累積寄与率の結果から次元を3つに絞ることにしたので、第三主成分(p3)までの結果を確認します。

def generate_e_vector_df(e_vector_matrix, pca_df_columns, num_of_dimension=None):
    """固有ベクトルをまとめたdfを作成する"""
    columns = ["p" + str(index+1) for index in range(e_vector_matrix.shape[0])]
    if num_of_dimension:
        columns = columns[:num_of_dimension]
        e_vector_matrix = e_vector_matrix[:, :num_of_dimension]

    e_vector_df = pd.DataFrame(e_vector_matrix, index=pca_df_columns, columns=columns)
    return e_vector_df

e_vector_df = generate_e_vector_df(e_vector_matrix, pca_df.columns, num_of_dimension=3)
e_vector_df
p1 p2 p3
land_area 0.133638 0.13902 -0.253528
forest_area 0.0894894 0.0948706 -0.335279
forest_rate 0.0724927 -0.128071 -0.238961
artificial_forest_area 0.2281 0.18997 0.0363571
forest_area_diff 0.186578 0.0717666 0.271166
forest_area_diff_rate 0.191038 -0.06318 0.219826
num 0.147522 0.249341 0.23795
estimated_num_2020 0.130906 0.255045 0.237969
estimated_num_2022 0.128582 0.256051 0.238496
birth_rate_2015 -0.236983 0.210662 -0.0105548
birth_rate_2020 -0.242097 0.199203 -0.0225353
death_rate_2015 -0.00961641 0.0398703 -0.400759
death_rate_2020 0.0387575 0.000478852 -0.415786
gdp_2015 0.252479 0.177535 -0.107411
gdp_2017 0.254531 0.179619 -0.100807
gdp_2018 0.256215 0.1829 -0.0911535
gdp_2019 0.254451 0.185363 -0.0914418
gdp_2020 0.255112 0.186521 -0.0822235
gdp_2021 0.256547 0.192383 -0.0688899
life_expectancy 0.222152 -0.264555 0.0186876
male_life_expectancy 0.214872 -0.258519 0.0601743
female_life_expectancy 0.221881 -0.268635 -0.00482029
healthy_life_expectancy 0.218476 -0.269892 0.027524
male_healthy_life_expectancy 0.211332 -0.268387 0.0488378
female_healthy_life_expectancy 0.215371 -0.272372 -0.00439609
population_density -0.0197248 0.0129451 0.285948

固有ベクトルは元の軸の重みを意味します。
具体的に、第一主成分(p1)では、陸地の面積(land_area)を約0.13倍, 森林の面積(forest_area)を約0.09倍, ..., 人口密度(population_density)を約-0.02倍していることになります。

しかし、軸(主成分)の説明をする時に"この軸は陸地の面積(land_area)を約0.13倍, 森林の面積(forest_area)を約0.09倍, ..., 人口密度(population_density)を約-0.02倍した軸"ですと言われてもどんな軸なのか理解しにくいと思います。

それぞれの軸(主成分)が"何を意味しているのか"は分析者が考える必要があります。
重みの関係を見ながら軸の意味を考えていきます。

元の軸の絶対値が大きな要素ほどその主成分に大きな影響を与えていることになるのですが、数値だけだと少し分かりにくいと思います。

一度ヒートマップで可視化した結果も見ながらそれぞれの主成分の意味を考えていきます。

import seaborn as sns
import matplotlib.pyplot as plt


def print_heat_map(df, vmax=1, vmin=0, center=0.5, cmap='Blues'):
    """
    ヒートマップを作成する
    デフォルトで最大1, 最小0, 中心0.5で作成
    (デフォルトでは確率の結果を表示することを想定)
    """
    vertical = len(df) / 4
    fig = plt.figure(figsize=(7, vertical))
    ax = sns.heatmap(df, xticklabels=1, yticklabels=1, vmax=vmax, vmin=vmin, center=center, linewidths=0.5, linecolor="Black", cmap=cmap)
    ax.tick_params(labelsize=9)
    plt.show()
    plt.close()
    return

def print_e_vector_heat_map(e_vector_df):
    """
    固有ベクトルごとにヒートマップを作成する
    中心を0にして絶対値の最大値で色のスケールを調整
    プラスは赤、マイナスが青となる表示とする
    """
    for column in e_vector_df:
        print(column)
        v_max = e_vector_df[column].abs().max()
        v_min = (-1) * v_max
        print_heat_map(pd.DataFrame(e_vector_df[column]), vmax=v_max, vmin=v_min, center=0, cmap="bwr")
    return

print_e_vector_heat_map(e_vector_df)

まず、p1(第一主成分)の結果の結果を見てみます。

p1

このヒートマップはプラスに大きいほど赤色が濃く、マイナスに大きいほど青色が濃くなるように作成しています。
赤や青の色が濃い要素ほど、p1(第一主成分)の結果にも大きな影響を及ぼしています。

具体的に見ていくと、gdp, 平均寿命(life_expectancy), 健康寿命(healthy_life_expectancy)は赤色が濃く、これらの値が大きいほどp1(第一主成分)の値も大きくなることを表します。
逆に出生率(birth_rate)は青色が濃くなっているので、出生率の値が小さいほどp1(第一主成分)の値が大きくなることを表します。

まとめるとp1(第一主成分)の値が大きいということは、gdpが大きく寿命も長いが、出生率が小さくなっていることを意味します。

直感的に、経済が発展しきった先進国を彷彿させる軸となっています。
ここでは"現状の経済力の高さ"と解釈することにします。


同様にp2(第二主成分)の結果を見てみます。

p2

この図から
p2(第二主成分)が大きいということは、人口(num)が多く出生率も高いが、寿命が短いことを意味します。
人口増加中の途上国を彷彿させる軸となっています。
ここでは"今後の経済発展力の高さ"と解釈することにします。


p3

p3(第三主成分)が大きくなるほど、死亡率(death_rate)が低く, 森の面積(forest_area)も狭くて人口密度(population_density)が高いということが読み取れます。

解釈が難しいですが、領土に対して人の数が多く死亡率も低いという理由でここでは"人口増加のしやすさ"と解釈することにします。

固有ベクトルと固有値を利用した元データの変換

固有ベクトルと固有値が求まっているので元データを変換します。
変換後のデータを主成分得点(Principal component score)と呼びます。

def compute_principal_component_score_matrix(data, e_values, e_vector_matrix, num_of_dimension=None):
    """主成分得点を計算する"""
    new_e_values = np.zeros_like(e_values)
    new_e_values[:num_of_dimension] = e_values[:num_of_dimension]
    ramda_matrix = np.diag(new_e_values)
    principal_component_score_matrix = data @ e_vector_matrix @ ramda_matrix
    return principal_component_score_matrix

principal_component_score_matrix = compute_principal_component_score_matrix(z_score_matrix, e_values, e_vector_matrix, num_of_dimension=3)
print(principal_component_score_matrix.shape) # (32, 26)

散布図の確認

主成分得点の結果を散布図で確認して見ます。まず散布図を作成するために便利なライブラリを使える状態にします。

# matplotlib日本語表示のインストール
!pip install japanize-matplotlib

# 散布図でテキストが重ならないようにするライブラリのインストール
!pip install adjustText
from adjustText import adjust_text
import japanize_matplotlib

japanize_matplotlib.japanize()

エラーが出なければOKです。
これでグラフ上に日本語が表示できるようになり、散布図上でテキストが被らないようにできます。

それでは実際に散布図を作成していきます。
絞った主成分のそれぞれの組み合わせで散布図を作成する仕様にしています。

import itertools

def print_sccater_plot_from_np_array(data, x_index, y_index, labels=None, data_names=None):
    """散布図を描画"""
    x = data[:, x_index]
    y = data[:, y_index]
    fig = plt.figure(figsize=(7, 5)) # Figureオブジェクトを作成
    ax = fig.add_subplot(1,1,1)
    ax.scatter(x, y)
    if labels is not None:
        ax.set_xlabel(labels[0])
        ax.set_ylabel(labels[1])

    if data_names is not None:
        texts = []
        for i, data_name in enumerate(data_names):
            texts.append(ax.text(x[i], y[i], data_name, fontsize=8))
        adjust_text(texts)
    plt.show()
    plt.close()
    return

def print_pca_result_sccater_plot(data, data_names, vector_name_dict):
    vecror_names = list(vector_name_dict.values())
    for labels in itertools.combinations(vecror_names, 2):
        print(labels)
        x_index = vecror_names.index(labels[0])
        y_index = vecror_names.index(labels[1])
        print_sccater_plot_from_np_array(data, x_index, y_index, labels=labels, data_names=data_names)
    return

vector_name_dict = {
    "p1": "p1_現状の経済力の高さ",
    "p2": "p2_今後の発展力の高さ",
    "p3": "p3_人口増加のしやすさ",
}

data_names = list(using_df.iloc[:, 0])
print_pca_result_sccater_plot(principal_component_score_matrix, data_names, vector_name_dict)

結果を見ていきます。


p1とp2の組み合わせ

解釈が正しければ右上にいればいるほど、
現状の経済力が高く、今後の発展力も見込めることになります。

アメリカや中国の強さがデータからも見えてくると思います。

日本も現状の経済力は3番目にいることが分かります。

アフリカ諸国が左上に多いことからも、アフリカ諸国の今後の発展が期待できそうです。


p1とp3の組み合わせ

p3ではアジア系の国が上に多く、今後も人口が増加していきそうなことが分かります。


p2とp3の組み合わせ

中国、インドは今後の発展も人口増加も期待できそうです。

最も気になる日本は、、、

不愉快なことを言うな。

と言われそうですが、

現実はもっと不愉快です!!

といった状況です。

挽回していきたいものです。

まとめ

基本的な数学の説明から、主成分分析の理論、データの取得と整形、実装までまとめました。
何かと便利な分析手法でおすすめです。

どなたかの参考になれば幸いです。

感想ですが、生データの分析はどんな結果になるか分からないので面白いです。

参考

正規化と標準化
https://qiita.com/oki_kosuke/items/02ec7bee92c85cf10ac2

固有値・固有ベクトルの計算
https://note.nkmk.me/python-numpy-matrix/#nplinalgeig

主成分分析の軸の数の判断
https://www.macromill.com/service/data-analysis/principal-component-analysis/?uuid=20240830_52d08a5e_9605_4477_8aca_e1c8900096b0

ヒートマップ
https://note.nkmk.me/python-seaborn-heatmap/

https://qiita.com/hik0107/items/67ad4cfbc9e84032fc6b

matplitlibの日本語表示
https://qiita.com/uehara1414/items/6286590d2e1ffbf68f6c

散布図で重なる文字の対策
https://qiita.com/mosomoso_1910/items/54b1f8da014676687d2b

itertools
https://qiita.com/Mohrey/items/b0afc8edc5fed742b68a

Goals Tech Blog

Discussion