皆さん、こんにちは！
この記事では、Pythonのscilit-learnを用いて、k-means法でクラスタリングを行っていきます。
記事の前半では、コードをコピペするだけで実行可能になるようにしています。
データの生成から実行していきますので、k-means法の基本的なステップを学習することが出来るはずです。
記事の後半では、私が用意したサンプルデータセット（NPB Batters Dataset(2008-2021)）の一部のデータを利用して、クラスタリングを行っています。
データの前処理や、3次元でのグラフ描画も行っているので、少しは実践的な内容になっていると思います。
Contents
k-means法とは？
まずは、k-means法について解説します。
k-means法は、非階層クラスター分析の代表的な手法の1つです。
要するに、与えられたデータセットを、似た者同士（クラスター）で分類するという方法になります。
では、何個のクラスターに分類するのかと言う部分が、「k」に該当し、こちらはユーザー側で定義してやる必要があります。
k-means法のアルゴリズムは以下の通りです。
- クラスターの「核」となるサンプルをk 個選択する。（初期値）
- 全てのサンプルとk個の「核」の距離を測る
- 各サンプルを最も近い「核」と同じクラスターに分類する
- 全てのサンプルをクラスターに分類できたら、再度、クラスターの重心点を求める（これが2イタレーション目の「核」になる）
- ②から⑤を、全てのサンプルとk個の核までの距離が変化しなくなるまで繰り返す。
k-means法を試してみる
それでは早速Pythonでk-means法を試してみましょう。
JupyterNotebookを使っている方は、この記事で使っているコードが全て含まれた.ipynbファイルを使って確認してもらっても大丈夫です。
ファイルアップロードの都合上Zipファイルになっています。
この記事で使うPythonライブラリは以下の通りです。
pip install pandas
pip install numpy
pip install matplotlib
pip install sklearn
pip install
pip install サンプルデータの生成
インストールが完了したら、次はサンプルデータの作成に取り掛かります。
サンプルデータは、sklearn.datasetsのmake_blobs関数を利用します。
make_blobs関数は、主に分類用のサンプルの作成に適したデータセットで、ガウス分布に基づいた形でランダムにデータセットを生成してくれます。
引数は以下の通り（この記事で使うもののみ記載）
| 引数 | 意味 |
|---|---|
| n_samples | データサンプルの総数 |
| centers | 各クラスターの重心座標 |
| cluster_std | クラスターの標準偏差 |
| random_state | ランダムシード（乱数を固定する） |
戻り値はこちらです。
| 戻り値 | 意味 |
|---|---|
| X | 生成された座標（ndarray） |
| y | 生成元のクラスター番号（ndarray） |
戻り値は、作成されたサンプルであるXと、生成元のクラスター番号（正解ラベル）であるyです。
簡単にサンプルを作成できるので、アルゴリズムの挙動を確認するときなどに便利ですね。
ちなみに、クラスター重心のcentersは、2次元以上も可能です。
import numpy as np
import pandas as pd
# クラスタリング用のデータを作成する
from sklearn.datasets import make_blobs
%matplotlib inline
import matplotlib.pyplot as plt
X, labels_true = make_blobs(
n_samples=1000, # サンプル数
centers=[[-5,-5],[0,0],[5,-5]],# クラスタ重心座標設定(3クラスタを2次元で作成)
cluster_std=1.0, # 乱数生成時の標準偏差
random_state=0)
print('生成されたデータ X:\n', X[:5])
print('生成されたデータ y:\n', labels_true[:5])生成されたデータをそれぞれ出力すると以下の通りです。
今回は2次元で作成したので、2次元で出てきています。
yの方は、Xのデータがそれぞれ、0番のクラスター中心から生成されたものなのか、1番目なのかを識別できるようになっています。
ただ、ndarrayの形式のままだと、扱いにくい部分があるので、pandasのデータフレームに成型しておきましょう。
hstackを使って、dataframeに変換できるよう形式にひとまとめにしていきます。
dataframeに成型出来たら、確認のためのに、matplotlibで散布図を作成しておきましょう。
df = pd.DataFrame(np.hstack([X,labels_true.reshape(len(labels_true),1)]),
columns=["X0","X1","label"])
col = df.label.map({0:'b', 1:'g', 2:'r'})
df.plot(x='X0', y='X1', kind='scatter', c=col, colorbar=False, figsize=(5,5))作成すると、このようにしっかりとクラスターに分けられたサンプルデータが作成できていることがわかりますね。
k-means法で学習を行う
サンプルを作成できた後は、k-means法で学習を行っていきましょう。
と言っても、学習だけなら2行で出来てしまいます。
残りのコードは、データを事前に成型するものだったり、結果を可視化したりするものなので、いかに前処理と後処理に工数がかかることがわかりますね。
k-means法は、KMeans関数で、学習のインスタンスを作成して、fitにデータを渡すだけです。
その際は、各変数を標準化しておくことをオススメします。
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
# kmeans法で、3つのクラスタを作成する
n_clusters = 3
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
kmeans.fit(X)
km_label =kmeans.labels_ # クラスタの分類結果ラベルを出力する
df["km_label"] = km_label # データフレームにクラスタの分類結果ラベルを渡す
print(kmeans.cluster_centers_[0], ':Greenクラスターの中心') # クラスターの中心を表す
print(kmeans.cluster_centers_[1], ':Blueクラスターの中心')
print(kmeans.cluster_centers_[2], ':Redクラスターの中心')
print(df)学習結果を表示すると、以下の様になります。
各クラスターの中心と、生成されたデータと、その正解ラベルを1つのデータフレームに収めています。
※label と km_labelの数字が異なっていますが、今回であればlabel = 0の時、km_labelが2であればOKです。
とはいえ、1000サンプルも1行ずつ追っていくのはしんどいので、プロットして確認しましょう。
各サンプルをkm_label基準で色分けしてプロットします。
colors = ['red', 'blue', 'green',"orange","black", "purple"]
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)
for i in range(n_clusters):
ax.scatter(df["X0"][df.km_label==i], df["X1"][df.km_label==i], color=colors[i], alpha=0.5)
#ax.scatter(centers[i, 0], centers[i, 1], marker='x', color=colors[i], s=300)
plt.show()表示するとこんな感じで、それぞれのクラスターごとに色分けがされていることがわかります。
k-means法で適切なクラスタ数を予測する
ただ、k-means法では、ユーザー側がクラスター数であるkをユーザー側が設定しなければならないのが、大きな問題です。
先ほどのような2次元のサンプルであれば、グラフで見ることができますが、実際は3次元、4次元など、複雑なデータを分類する必要があるのが普通でしょう。
これを簡単に見積もることができる方法が「エルボー法」です。
エルボー法とは
k-means法のゴールは、各クラスター内の中心からの距離が変化しなくなる（最短）になることを目指します。
正しくは二乗誤差（inertia）を最小化することが目的ですが、どれくらい小さくすれば良いという絶対的な指標はありません。
あくまでデータによる。と言うことです。
もう1つは、k-means法自体が等方性（例えば分布が円状など、偏りがない）ことが前提になっており、ある次元から見た分布が細長い場合など異方性のデータの分類が苦手です。
とはいえ、上記はすべてk-meansの方法です。
エルボー法は、kをループで変化させ、グラフにX：k、Ｙを二乗誤差でプロットし、どのkが最適かを探索する方法全般を指します。
エルボー法でプロットしてみる
それでは、エルボー法をプロットしていきましょう。
先ほどのサンプルデータを使って、クラスター数を10まで1ずつ増やしていきます。
# エルボー図の作成
distortions = [] # 空のリストを作成
for i in range(1,11): # 10クラスタまで、順に実施
km=KMeans(n_clusters=i, random_state=0).fit(X)
distortions.append(km.inertia_)
plt.plot(range(1,11), distortions, marker='o')プロットするとこんな感じ、定義上、クラスター数を増やせば増やすほど、二乗誤差は小さくなりますが、サンプル生成時に定義したk = 3 までが二乗誤差の減り方が大きく、分類が効果的であることがわかります。
クラスタ数を6にしてみると…？
ちなみに、クラスター数を6にして増やしみると、このように1つのクラスタを更に2分するような分類がされます。
ただ、1つのクラスタを縦に切るのか、横に切るのか、どちらが適切かは微妙なところですし、クラスター数を増やしすぎることで、かえってわかりにくい分類になってしまうことがわかります。
異常で、k-means法の基本的な流れの解説を終わります！
サンプルデータは元から分類してあるようなものだったので、少々物足りなかったかもしれませんね。
次からは実際のデータを使って、分類にチャレンジしてみましょう。
例題：NPBの野手シーズンデータを使って打者タイプを分類しよう！
サンプルデータ
以下の分析では、私が別の記事で紹介している野球のデータを使っていきます。
データ取得に興味がある方は、以下の記事にアクセスしてトレースしてみてください。
ひとまずデータだけ使いたい！
と言う方は、以下のリンクからCSVデータをダウンロードできますので、ダウンロードしてお使いください。
データを開いて操作したい場合、Google Spread Sheet がおすすめです。
データのインポート
上記リンクからCSVをダウンロードしたら、分析を行う、Pythonファイル（私はJupyter Notebookファイル）と同じ階層においてください。
まずはデータをインポートして、どんなデータが入っているか確認するのがおすすめです。
import pandas as pd
%matplotlib inline
import matplotlib.pyplot as plt
# データセットを読み込んでみる
df = pd.read_csv("stats_summary.csv")
df.head()df.head()で、上から5件のみの表示ですが、このように野手のシーズンデータが入っていることがわかります。
（ちなみに4358行あります。）
データの前処理
簡単にですが、データの前処理を行っていきましょう。
今回は、打率、長打率、三振率を使って、打者のタイプ分けを行っていきたいと思います。
前処理では、データフレームを先ほど挙げた変数に絞り込んでいきます。
更に、「三振率」は元のデータセットにないため、三振数を打席数で割って作成することにします。
# データフレームを今回の分析対象に限定
df = df[["選手名","打率",'打席数','三振', '長打率']]
# ”三振率”の列を作成
df["三振率"] = round(df["三振"] / df["打席数"],3)
# 後々の可視化のために、カラム名を半角英数字表記に変更
columns = ("Player Name","AVG", "At Bats", "Strike Outs","SLG", "K%")
df.columns = columns
df.head()
データフレームを小さくすることが出来ました。
データセットを用意したら、describe()を使って、各列の特徴を把握しておきましょう。
上から、データ数、平均値、標準偏差、最小値、下から25％、50％、75％の値、最大値です。
数に着目すると、25％の値がたったの24打席と、打者タイプをカテゴライズするのには小さすぎる値です。
打席数が少ないと、打率10割など、極端な値（ノイズ）が出てしまうことが懸念されます。
規定打席以上でも良いですが、これだと打撃に優れた選手ばかり残ってしまうことも考えられますので、今回はザックリ300打席以上の選手に絞り込みます。（300打席未満は切り捨ててしまいます）
300打席未満を切り捨てて、再度describe()してみましょう。
せっかくなので、このデータを散布図行列として表示してみることにします。
あまり色んなライブラリを使うのもどうかと思いますが、seabornの相関行列が便利すぎるので、使います。
たった2行で、先ほどのデータフレームをバシッと可視化してくれます。
import seaborn as sns
sns.pairplot(df)打率は正規分布のようなきれいな山なりになっており、打席数を除けば、それなりに多様な分布になっていることが分かります。（やはり打率と長打率は正の相関がありそうですが、今回はサンプルなので…）
三次元で表示してみる
また、今回は3つの変数を使ってクラスタリングをしていくため、3次元のグラフで表示してみましょう。
ついでに、先ほどのデータフレームを学習用データセットXに格納しておきます。
# 学習用データセットを作成
X = df[["AVG","SLG","K%"]]
from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d.axes3d import Axes3D
x = X["AVG"]
y = X["SLG"]
z = X["K%"]
# figureを生成する
fig = plt.figure(figsize=(10, 8))
# axをfigureに設定する
ax = Axes3D(fig)
# axesに散布図を設定する
ax.scatter(x, y, z,c='b')
ax.set_xlabel('AVG')
ax.set_ylabel('SLG')
ax.set_zlabel('K%')
# 表示する
plt.show()3次元表示すると、こんな感じでXYZに分布していることがわかります。
最適なkを探索する
早速ですが、エルボー法を使って、最適なクラス数を見積もってみます。
# エルボー図の作成
distortions = [] # 空のリストを作成
# 学習用データをXに格納
X = df[["AVG","SLG","K%"]]
for i in range(1,11): # 10クラスタまで、順に実施
km=KMeans(n_clusters=i, random_state=0).fit(X)
distortions.append(km.inertia_)
plt.plot(range(1,11), distortions, marker='o')
結果は、こちらです。
各クラスタ数の比較は私が裏でやってしまいましたので、皆さんの環境でも色々試してみてください。
ここでは私が一番よさそうと判断したk = 5 を採用して、可視化までを行っていきます。
クラスタごとに色分けもしたいため、分類結果であるlabelsもデータフレームに追記するようにしておきます。
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
n_clusters = 5
# kmeans法
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
kmeans.fit(X)
kmeans =kmeans.labels_ # クラスタのラベル
X["labels"] = kmeans
X
結果をプロット
それでは結果をプロットしていきます。
あらかじめcolors のリストに色を保持しておき、クラスターの数だけループ分でプロットしていきます。
3次元プロットすると、何やら分類できてそうな雰囲気がありますが、斜め視点だとよくわかりません。
colors = ['red', 'blue', 'green',"orange","black", "purple"]
fig = plt.figure(figsize=(10, 8))
# axをfigureに設定する
ax = Axes3D(fig)
for i in range(n_clusters):
ax.scatter(X["AVG"][X.labels==i], X["SLG"][X.labels==i], X["K%"][X.labels==i],color=colors[i], alpha=0.5)
ax.set_xlabel('AVG')
ax.set_ylabel('SLG')
ax.set_zlabel('K%')
plt.show()
それでは、上記にプロットした3次元グラフを、真横や真上から見た視点に変更していきます。
変更方法は、plt.show()の前にax.view_init(elev = 〇〇, azim = 〇〇)を記入しておくと視点を変更できます。
ax.view_init(elev=-90, azim=0)
ax.view_init(elev=0, azim=-90)
ax.view_init(elev=0, azim=0)
分類結果
分類結果をまとめると、以下の様に言語化できそうです。
ラベル0（赤）：高三振・中長打のフリースインガー型
ラベル1（青）：低長打・低打率の下位打線の選手
ラベル2（緑）：低三振・中長打のアベレージヒッター
ラベル3（橙）：低長打・中打率のリードオフマン型
ラベル4（黒）：高打率・高長打の主力打者（クリーンアップなど）
シーズンデータのうち、打率・長打率・三振率という非常に少ない変数ですが、概ね打者の特徴を捉えられてそうです。
規定打席未満の選手を入れたこともあり、入れ替わりが激しいと考えられる下位打線の選手をうまくカテゴライズできているような気がします。
チームとしては、中軸にラベル4（黒）の選手を配置しつつ、如何にラベル1（青）の選手を減らしつつ、ラベル2（緑）やラベル3（橙）の選手を増やすことで厚みのある打線を組めそうです。
打線のアクセントとしてはラベル0（赤）欲しいところですが、6番や7番あたりに1人くらいいれば十分かと思います。
ラベル0（赤）：高三振・中長打のフリースインガー型
df[df.labels == 0].head(20)
ラベル2（緑）：低三振・中長打のアベレージヒッター
df[df.labels == 2].head(20)
ラベル4（黒）：高打率・高長打の主力打者（クリーンアップなど）
df[df.labels == 4].head(20)
ラベル1（青）：低長打・低打率の下位打線の選手
df[df.labels == 1].head(20)
ラベル3（橙）：低長打・中打率のリードオフマン型
df[df.labels == 3].head(20)
なお、今回の分析では、選手の走力や守備力は一切含んでいません。
打撃成績に関しては、四球など選球眼に関する指標も三振に一部考慮されていますが、そこまで重視しているわけではありません。
あくまで打者のカテゴライズと言うことで、それなりにイメージと合致する分類が出来たと思っています。
まとめ
この記事では、PythonのScikit-learnを使って、k-means法をデータの生成から試してみました！
また、後半では、NPBのシーズン成績を利用して、打者のカテゴライズに挑戦しました。
皆さんもこの記事のデータを使って試してみてください！！
参考文献
参考文献を示します。
可視化
k-meansの最適なクラスター数を調べる方法
https://qiita.com/deaikei/items/11a10fde5bb47a2cf2c2
3D scatterplot（Matplotlib ドキュメント）
https://matplotlib.org/3.5.0/gallery/mplot3d/scatter3d.html#sphx-glr-gallery-mplot3d-scatter3d-py
Scikit-learnについて
sklearn.cluster.KMeans（scikit-learn ドキュメント）
https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
K-means Clustering（scikit-learn ドキュメント）
3D scatterplot
https://matplotlib.org/3.5.0/gallery/mplot3d/scatter3d.html#sphx-glr-gallery-mplot3d-scatter3d-py
