k-means法クラスター分析
クラスター分析をk-means法で試みた。データは幾つかの都市の平年降水量であり、図１に示すものである。国土交通省のウェブサイト
https://www.mlit.go.jp/river/toukei_chousa/kasen/jiten/toukei/birn38p.html
に掲載されているものから値を読み込んで、Excelファイル（rain.xlsx）とした（2024年6月作成）。
図１ 平年降水量（rain.xlsx）
k-means法では、クラスター数を指定するか、初期値としてクラスターの中心を指定するが、まず、クラスター数を指定する方法を試みる。その後、初期値としてクラスター中心を指定する方法を試みる。
図１のデータに対してWard法によるクラスター分析を行ったものが、このウェブサイトに示されている。そこに示されている樹状図からクラスター数を６個としてk-means法を試みることにした。スクリプトはリスト１に示すものである。スクリプトファイルなどは、まとめて圧縮ファイルkmeansfiles.zipとした。
リスト１のスクリプトを実行すると、クラスター分析の結果を示すための座標軸の向きを決めるため、ホワイト処置後のデータの布置の中心化を行った後、固有分解を行った結果の固有値のスクリープロットが示される（図２）。ホワイト処置については、本ウェブサイトの後半において少し検討を行う。
図２
２つの固有値が他より大きいので、２次元の散布図で布置の大まかな様子を表すことが出来ると考えられる。
図２のウィンドウを閉じると、固有ベクトルのグラフが表示される（図３）。
図３
固有ベクトルのグラフを閉じると、クラスターの中心と各都市の回転後の座標値の散布図が表示される（図４）。
図４
都市名の先頭に付されている数字は、その都市の属するクラスターを表している。図４は２次元に投影された分布であるので、クラスター中心と都市との距離は投影された平面上での見かけの距離であることに注意する必要がある。
図４のウィンドウを閉じると、端末に次に表示する散布図の次元の設定を求めるメッセージが表示される。
Input the dimensions, larger than 0
and smaller than 12.
dim_x = 1
dim_y = 3
上の例では、X軸として第１次元、Y軸として第３次元が設定されている。X軸とY軸に割り当てる次元を設定すると、設定した次元で散布図が描かれる。次元数として、０以下、あるいは存在しない大きい次元数を設定するとスクリプトの実行終了となる。
図４に示された散布図の右下方におけるクラスター中心の分布がやや密であるように見えるので、このことを踏まえてクラスター中心の初期値を設定してk-means法を行ったものを次に示す。
スクリプトをリスト２に示す。スクリプトファイルなどは、まとめて圧縮ファイルkmeansfiles.zipとした。
図４の散布図から、クラスター数を５とし、５つのクラスターの中心の初期値を、新潟、山形、東京、静岡、長崎とすることにした。
コードは、以下のように用意した。
init_c = {'新潟', '山形',
'東京', '静岡', '長崎'}
k = 5
init_cities = []
init_coords = []
for i in range(len(cities)):
if cities[i] in init_c:
init_cities.append(cities[i])
init_coords.append(data_wh[i])
init_coords = np.array(init_coords)
上のクラスター中心の初期値init_coordsを用いて、次のコードでk-means法を実行した。
results = kmeans(data_wh, init_coords)
リスト２のスクリプトを実行すると、まず、散布図を描くときの軸の方向を決める固有分解の固有値のスクリープロットが表示される（図５）。
図５
２つの固有値が他のものより大きいので、２次元の布置で分布の大まかな様子が表されると考えられる。
図５のウィンドウを閉じると、固有ベクトルのグラフが表示される（図６）。
図６
図５のウィンドウを閉じると、クラスター中心と都市の散布図が表示される（図７）。
図７
都市名の先頭に付されている数字は、その都市の属するクラスターを表している。図４は２次元に投影された分布であるので、クラスター中心と都市との距離は投影された平面上での見かけの距離であることに注意する必要がある。
図７のウィンドウを閉じると、端末に次に表示する散布図の次元の設定を求めるメッセージが表示される。
Input the dimensions, larger than 0
and smaller than 12.
dim_x = 1
dim_y = 3
上の例では、X軸として第１次元、Y軸として第３次元が設定されている。X軸とY軸に割り当てる次元を設定すると、設定した次元で散布図が描かれる。次元数として、０以下、あるいは存在しない大きい次元数を設定するとスクリプトの実行終了となる。
ホワイト処理（whitening）において、データ値が標準偏差で割られる。これにより、標準偏差は１になるが、平均値は標準偏差で割られた値になる。ホワイト処理の検討を行うリスト３のスクリプトにおいて、次のコードによりホワイト処理前のデータと処理後のデータの比較を行った。
for i in range(12):
print(f'{i+1}月:')
print('Before: Mean =',
mean_raw[i], ' SD =',
sd_raw[i])
print('After: Mean =', mean_w[i], ' SD =', sd_w[i])
print('mean_raw/sd_raw =',
mean_raw[i]/sd_raw[i])
リスト３を実行すると、以下の出力を得る。
1月:
Before: Mean = 79.27872340425533 SD = 61.486867430821434
After: Mean = 1.2893602604401895 SD = 1.0
mean_raw/sd_raw = 1.2893602604401895
2月:
Before: Mean = 82.45106382978723 SD = 37.47995136319437
After: Mean = 2.1998711532682216 SD = 0.9999999999999998
mean_raw/sd_raw = 2.199871153268221
3月:
Before: Mean = 112.24680851063829 SD = 35.74835720709254
After: Mean = 3.1399151535939205 SD = 0.9999999999999998
mean_raw/sd_raw = 3.13991515359392
4月:
Before: Mean = 129.3127659574468 SD = 44.42636085881883
After: Mean = 2.910721550396304 SD = 1.0
mean_raw/sd_raw = 2.9107215503963033
5月:
Before: Mean = 141.86170212765958 SD = 45.94784533272244
After: Mean = 3.087450588822947 SD = 1.0
mean_raw/sd_raw = 3.087450588822947
6月:
Before: Mean = 209.85744680851067 SD = 89.8300682842293
After: Mean = 2.3361603838984673 SD = 0.9999999999999999
mean_raw/sd_raw = 2.3361603838984677
7月:
Before: Mean = 199.80851063829786 SD = 66.79021724068876
After: Mean = 2.9915834817284295 SD = 0.9999999999999998
mean_raw/sd_raw = 2.9915834817284295
8月:
Before: Mean = 159.55957446808506 SD = 48.0064608221642
After: Mean = 3.3237104284600316 SD = 1.0
mean_raw/sd_raw = 3.3237104284600307
9月:
Before: Mean = 207.54255319148928 SD = 53.4373036676482
After: Mean = 3.8838515221930803 SD = 1.0
mean_raw/sd_raw = 3.883851522193079
10月:
Before: Mean = 123.3957446808511 SD = 33.210354431513174
After: Mean = 3.7155804806395367 SD = 1.0
mean_raw/sd_raw = 3.7155804806395367
11月:
Before: Mean = 97.97446808510637 SD = 48.643131496012295
After: Mean = 2.0141480425276117 SD = 0.9999999999999999
mean_raw/sd_raw = 2.0141480425276113
12月:
Before: Mean = 72.68085106382979 SD = 65.18007996932504
After: Mean = 1.1150776601997843 SD = 0.9999999999999998
mean_raw/sd_raw = 1.1150776601997843
ホワイト処理後、標準偏差が１に揃えられ、平均値は標準偏差で割られた値になっていることが確認できる。
リスト３のスクリプトを実行すると、ホワイト処理前のデータの分布を示す図W1のグラフが表示される。
図W1．ホワイト処理前の分布
３月と10月の標準偏差が小さく、６月の標準偏差が大きい。
図W1のグラフの表示されているウィンドウを閉じると、ホワイト処理後の分布を示すグラフが表示される（図W2）。
図W2．ホワイト処理後の分布
標準偏差が１に揃えられていることが分かる。図W1と図W2を比べると、３月と10月の分布が引き延ばされ、６月の分布が圧縮されて、クラスター分析が行われていることになる。３月と10月の分布が引き延ばされると、これらの月における違いの影響は大きくなり、６月における分布が圧縮されるとこの月における違いの影響は小さくなる。ホワイト処理を行うかどうかは、分析の目的に応じて決められるべきであろう。
リスト１ クラスター数を指定する場合
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.cluster.hierarchy import
linkage
from scipy.cluster.hierarchy import
dendrogram
from scipy.cluster.vq import kmeans,
whiten
# ExcelファイルをPandasのDataFrame型として読み込む
fl_nm = 'rain.xlsx'
with pd.ExcelFile(fl_nm) as xlfile:
df = pd.read_excel(xlfile)
cities = df.values[:,0]
data_v = np.array(df.values[:,1:],
dtype=float)
k = 6
data_wh = whiten(data_v)
#
#
k-meansクラスター分析
#
results = kmeans(data_wh, k)
print('results =\n', results)
#
#
描画の準備のための固有分解
#
#
中心化
Z = data_wh - np.mean(data_wh, axis=0)
# 主軸のウェイト算出のための固有分解
vals, vecs = np.linalg.eigh(Z.T @ Z)
# 固有値の大きさの降順に並べ換える
idx = np.argsort(-vals)
vals = vals[idx]
vecs = vecs[:,idx]
# 固有値のスクリープロット
plt.plot(range(len(vals)), vals)
plt.title('Eigen Values', fontsize=16)
plt.tight_layout()
plt.show()
# 直交回転
data_wh_rot = data_wh @ vecs
k_means_rot = results[0] @ vecs
print('k_means_rot =\n', k_means_rot)
#
# 対象と各クラスター中心までの距離
#
d_obj_clus = np.empty((len(cities),
k))
for i in range(len(cities)):
for j in range(k):
d_obj_clus[i][j] = (np.sum((data_wh_rot[i] - k_means_rot[j])**2))**0.5
#
対称の所属クラスターを求める
obj_clus = np.empty(len(cities))
for i in range(len(cities)):
obj_clus[i] =
np.argmin(d_obj_clus[i])
#
# 回転用ウェイトのプロット
#
for i in range(len(vecs)):
plt.plot(range(len(vecs.T)),
vecs.T[i], label=f'w[{i+1}]')
plt.title('Weights for the Rotation',
fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
#
# 対象とクラスター中心の散布図
#
dim_x = 1
dim_y = 2
while True:
plt.plot(k_means_rot[:,dim_x-1], k_means_rot[:,dim_y-1], linewidth=0,
marker = '*', markersize=25, alpha=0.7, label='Cluster')
plt.plot(data_wh_rot[:,dim_x-1], data_wh_rot[:,dim_y-1], linewidth=0,
marker='o', label='Object')
for i in range(len(cities)):
plt.text(data_wh_rot[i,dim_x-1], data_wh_rot[i,dim_y-1],
f'{int(obj_clus[i])+1}'+ cities[i])
plt.xlabel(f'Axis-{dim_x}')
plt.ylabel(f'Axis-{dim_y}')
plt.legend()
plt.show()
print('Input the dimensions,
larger than 0 and smaller than' +
f' {len(data_wh[0])}.')
dim_x = int(input('dim_x =
'))
if (dim_x < 1) or (dim_x
> len(data_wh[0])):
break
dim_y = int(input('dim_y =
'))
if (dim_y < 1) or (dim_y
> len(data_wh[0])):
break
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.cluster.hierarchy import
linkage
from scipy.cluster.hierarchy import
dendrogram
from scipy.cluster.vq import kmeans,
whiten
# ExcelファイルをPandasのDataFrame型として読み込む
fl_nm = 'rain.xlsx'
with pd.ExcelFile(fl_nm) as xlfile:
df = pd.read_excel(xlfile)
cities = df.values[:,0]
data_v = np.array(df.values[:,1:],
dtype=float)
data_wh = whiten(data_v)
init_c = {'新潟', '山形',
'東京', '静岡', '長崎'}
k = 5
init_cities = []
init_coords = []
for i in range(len(cities)):
if cities[i] in init_c:
init_cities.append(cities[i])
init_coords.append(data_wh[i])
init_coords = np.array(init_coords)
print('init_clusters,,,\n')
for v in zip(init_cities,
init_coords):
print('city =', v[0], ' coord =', v[1])
#
#
k-meansクラスター分析
#
results = kmeans(data_wh, init_coords)
#
#
描画の準備のための固有分解
#
Z = data_wh - np.mean(data_wh, axis=0)
vals, vecs = np.linalg.eigh(Z.T @ Z)
# 固有値の大きさの降順に並べ換える
idx = np.argsort(-vals)
vals = vals[idx]
vecs = vecs[:,idx]
#
# 固有値のスクリープロット
#
plt.plot(range(len(vals)), vals)
plt.title('Eigen Values', fontsize=16)
plt.tight_layout()
plt.show()
#
# 直交回転
#
data_wh_rot = data_wh @ vecs
k_means_rot = results[0] @ vecs
#
# 対象と各クラスター中心までの距離
#
d_obj_clus = np.empty((len(cities),
k))
for i in range(len(cities)):
for j in range(k):
d_obj_clus[i][j] = (np.sum((data_wh_rot[i] - k_means_rot[j])**2))**0.5
# 対象から最小距離のクラスター
obj_clus = np.empty(len(cities))
for i in range(len(cities)):
obj_clus[i] =
np.argmin(d_obj_clus[i])
#
# ウエイトのグラフ
#
for i in range(len(vecs)):
plt.plot(range(len(vecs.T)),
vecs.T[i], label=f'w[{i+1}]')
plt.title('Weights for the Rotation',
fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
#
#
対象とクラスター中心の散布図
#
dim_x = 1
dim_y = 2
while True:
plt.plot(k_means_rot[:,dim_x-1],
k_means_rot[:,dim_y-1], linewidth=0,
marker = '*', markersize=25, alpha=0.7, label='Cluster')
plt.plot(data_wh_rot[:,dim_x-1], data_wh_rot[:,dim_y-1], linewidth=0,
marker='o', label='Object')
for i in range(len(cities)):
plt.text(data_wh_rot[i,dim_x-1], data_wh_rot[i,dim_y-1],
f'{int(obj_clus[i])+1}'+ cities[i])
plt.xlabel(f'Axis-{dim_x}')
plt.ylabel(f'Axis-{dim_y}')
plt.legend()
plt.show()
print('Input the dimensions,
larger than 0 and smaller than' +
f' {len(data_wh[0])}.')
dim_x = int(input('dim_x =
'))
if (dim_x < 1) or (dim_x
> len(data_wh[0])):
break
dim_y = int(input('dim_y =
'))
if (dim_y < 1) or (dim_y >
len(data_wh[0])):
break
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.cluster.vq import kmeans,
whiten
# ExcelファイルをPandasのDataFrame型として読み込む
fl_nm = 'rain.xlsx'
with pd.ExcelFile(fl_nm) as xlfile:
df = pd.read_excel(xlfile)
data_v = np.array(df.values[:,1:],
dtype=float)
# Whitening
data_wh = whiten(data_v)
#
# 各月の平均値と標準偏差
#
mean_raw = np.mean(data_v, axis=0)
sd_raw = np.std(data_v, axis=0)
mean_w = np.mean(data_wh, axis=0)
sd_w = np.std(data_wh, axis=0)
print('Comparison between raw data and
whitened data')
for i in range(12):
print(f'{i+1}月:')
print('Before: Mean =',
mean_raw[i], ' SD =',
sd_raw[i])
print('After: Mean =', mean_w[i], ' SD =', sd_w[i])
print('mean_raw/sd_raw =',
mean_raw[i]/sd_raw[i])
max_v = np.max(data_v)
print('max_v =', max_v)
#
# ヒストグラムの境界値の設定
#
v_bins = np.arange(0, 450, 25)
fig, ax = plt.subplots(3, 4,
sharex=True, sharey=True, figsize=(10,7))
for i in range(12):
r_idx = i // 4
c_idx = i % 4
ax[r_idx,
c_idx].hist(data_v.T[i], bins=v_bins)
ax[r_idx,
c_idx].set_title(f'{i+1}月(SD={sd_raw[i]:.2f})')
ax[r_idx, c_idx].set_xlabel('mm')
ax[r_idx,
c_idx].set_ylabel('度\n\n数', rotation=0, va='center')
fig.suptitle('平年降水量', fontsize=18)
plt.tight_layout()
plt.show()
max_vw = np.max(data_wh)
print('max_vw =', max_vw)
#
# ヒストグラムの境界値の設定
#
v_bins_w = np.arange(0, 8, 0.5)
fig, ax = plt.subplots(3, 4,
sharex=True, sharey=True, figsize=(10,7))
for i in range(12):
r_idx = i // 4
c_idx = i % 4
ax[r_idx,
c_idx].hist(data_wh.T[i], bins=v_bins_w)
ax[r_idx,
c_idx].set_title(f'{i+1}月(SD={sd_w[i]:.2f})')
ax[r_idx,
c_idx].set_xlabel('SD')
ax[r_idx,
c_idx].set_ylabel('度\n\n数', rotation=0, va='center')
fig.suptitle('平年降水量(Whitened)', fontsize=18)
plt.tight_layout()
plt.show()
なお、Pythonデータ分析の入門書＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を上梓した。