前回は、決定木を使って、質問を順番に重ねながら分類する方法を見ました。

今回は、クラスタリング を扱います。クラスタリングは、正解ラベルが付いていないデータを見て、「似ているもの同士」を自動的にグループへ分ける方法です。

ファインマン風に言えば、クラスタリングは「名前の付いていない箱に、似ているものを自然に集めていく作業」です。最初から 優良顧客 や 離反しそうな顧客 という答えがあるわけではありません。データの形を見ながら、「このあたりの点はまとまっていそうだ」と探していきます。

目次

クラスタリングで何をしたいのか
この記事で使うデータについて
K-means の基本原理
距離と類似度：何をもって「似ている」と見るか
NumPy で K-means を実装する
Scikit-learn で K-means を使う
階層クラスタリング
階層クラスタリングを距離行列で追う
Scikit-learn で階層クラスタリングを使う
DBSCAN：密度でクラスタを見つける
Scikit-learn で DBSCAN を使う
ガウス混合モデル：確率分布としてクラスタを見る
EM アルゴリズムの考え方
Scikit-learn でガウス混合モデルを使う
ケース：アジアサッカーチーム風の指標をクラスタリングする
まとめ

1. クラスタリングで何をしたいのか

クラスタリングは、教師なし学習 の代表的な方法です。

教師あり学習では、データに正解ラベルがあります。たとえば、住宅価格を予測するなら 価格、ローン審査を分類するなら 承認 や 見送り のようなラベルがあります。

一方、クラスタリングでは、最初から正解ラベルがありません。あるのは特徴量だけです。

よく使われる場面は、次のようなものです。

顧客を購買傾向ごとに分ける
商品を売れ方や価格帯で分ける
Webアクセスログから似た行動パターンを探す
画像や文章を似た特徴ごとに整理する

ここで注意したいのは、クラスタ番号そのものには意味がないことです。cluster=0 が偉い、cluster=2 が悪い、という意味はありません。番号は、アルゴリズムが便宜的に付けた名前です。

大事なのは、分かれたあとに各クラスタの特徴を見て、こちらが意味を読み取ることです。

クラスタリングの目的が見えたところで、この記事で使うデータを用意します。

2. この記事で使うデータについて

この記事では、外部ファイルに依存しないように、説明用の小さな架空データと、Scikit-learn で作れる合成データを使います。

K-means と階層クラスタリングでは、架空の顧客データ customer_jp を使います。
DBSCAN では、曲がった形のクラスタを作れる make_moons を使います。
ガウス混合モデルでは、複数の山を持つデータを作れる make_blobs を使います。
最後のケースでは、架空のチーム指標データを使います。

まず、顧客データを作ります。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler


customer_jp = pd.DataFrame({
    "顧客ID": [
        "C001", "C002", "C003", "C004",
        "C005", "C006", "C007", "C008",
        "C009", "C010", "C011", "C012",
    ],
    "年間購入額_万円": [8, 10, 12, 9, 45, 48, 52, 50, 88, 92, 96, 90],
    "来店回数": [3, 4, 5, 4, 18, 20, 21, 19, 24, 26, 27, 25],
})

print(customer_jp)

各列の意味は次の通りです。

列名	意味
`顧客ID`	顧客を区別するためのID。クラスタリングには使いません。
`年間購入額_万円`	1年間の購入額。単位は万円です。
`来店回数`	1年間の来店回数です。

このデータは説明用に作った架空データです。実在の購買履歴ではありません。

クラスタリングでは、距離を使うことが多いため、特徴量のスケールをそろえておくことが重要です。年間購入額_万円 と 来店回数 は単位が違うので、ここでは標準化してから使います。

たとえば標準化をしないと、年間購入額_万円 は8から96まで動く一方、来店回数 は3から27までしか動きません。ユークリッド距離のような計算では、値の変動幅が大きい特徴量ほど距離への影響が強くなるため、標準化をしないままだと、実質的に 年間購入額_万円 の違いだけでクラスタが決まってしまいます。

feature_cols = ["年間購入額_万円", "来店回数"]

X_customer = customer_jp[feature_cols]
scaler = StandardScaler()
X_customer_scaled = scaler.fit_transform(X_customer)

plt.scatter(X_customer_scaled[:, 0], X_customer_scaled[:, 1])
plt.xlabel("年間購入額_万円（標準化後）")
plt.ylabel("来店回数（標準化後）")
plt.title("顧客データの散布図")
plt.show()

この散布図を見ると、いくつかの点のまとまりがありそうです。次は、このまとまりを自動で見つける代表的な方法である K-means を見ていきます。

3. K-means の基本原理

K-means は、データをあらかじめ決めた $K$ 個のクラスタに分ける方法です。

中心になる考え方は、各クラスタの中心点に近いサンプルを、そのクラスタへ入れる ことです。この中心点を重心または セントロイド と呼びます。

サンプル $x_i$ をどのクラスタへ入れるかは、次のように書けます。

$$
c_i=\arg\min_{k\in{1,\ldots,K}}|x_i-\mu_k|^2
$$

$\mu_k$ は、第 $k$ クラスタのセントロイドです。

クラスタに入ったサンプルが決まったら、セントロイドを平均で更新します。

$$
\mu_k=\frac{1}{|C_k|}\sum_{x_i\in C_k}x_i
$$

K-means が小さくしようとしている量は、クラスタ内平方和です。

$$
J=\sum_{k=1}^{K}\sum_{x_i\in C_k}|x_i-\mu_k|^2
$$

具体的な数字で、1回分の「割り当て→更新」の流れを追ってみます。次の4個の点があるとします。

点	座標
$P_1$	$(1,\ 1)$
$P_2$	$(2,\ 1)$
$P_3$	$(8,\ 8)$
$P_4$	$(9,\ 9)$

$K=2$ とし、初期セントロイドを $\mu_1=(1,1)$（$P_1$ と同じ位置）、$\mu_2=(8,8)$（$P_3$ と同じ位置）とします。

まず割り当てです。各点から $\mu_1$、$\mu_2$ までの距離の2乗を計算すると、次のようになります。

点	$\|x-\mu_1\|^2$	$\|x-\mu_2\|^2$	割り当て
$P_1$	0	98	クラスタ1
$P_2$	1	85	クラスタ1
$P_3$	98	0	クラスタ2
$P_4$	128	2	クラスタ2

次に更新です。クラスタ1は $P_1,P_2$ の平均、クラスタ2は $P_3,P_4$ の平均でセントロイドを動かします。

$$
\mu_1=\left(\frac{1+2}{2},\ \frac{1+1}{2}\right)=(1.5,\ 1),
\qquad
\mu_2=\left(\frac{8+9}{2},\ \frac{8+9}{2}\right)=(8.5,\ 8.5)
$$

新しいセントロイドで距離を計算し直しても、割り当ては変わりません。つまりこの例では、もう1回の「割り当て→更新」で収束したことになります。実際のデータではもっと繰り返しが必要になりますが、やっていることはこの4点の例と同じです。

ファインマン風に言えば、K-means は「仮の集合場所を置く、近い人を集める、集まった人たちの真ん中へ集合場所を動かす」を繰り返す方法です。

K-means は距離を見て動くので、次に「距離」や「似ている」の測り方を整理します。

4. 距離と類似度：何をもって「似ている」と見るか

クラスタリングでは、「似ている」を数字にする必要があります。代表的な見方には、距離と類似度があります。

4.1 Minkowski 距離

Minkowski 距離は、いろいろな距離をまとめて表せる形です。

$$
d(x,y)=\left(\sum_{i=1}^{n}|x_i-y_i|^p\right)^{1/p}
$$

$p=2$ のときは、よく使われるユークリッド距離になります。

$$
d(x,y)=\sqrt{\sum_{i=1}^{n}(x_i-y_i)^2}
$$

4.2 余弦類似度

余弦類似度は、2つのベクトルの向きがどれくらい似ているかを見ます。

$$
s(x,y)=\cos\theta
=\frac{x^Ty}{|x||y|}
$$

値が 1 に近いほど、向きが似ています。

4.3 Pearson 相関係数

Pearson 相関係数は、2つの変数がどれくらい同じように増減するかを見ます。

$$
\rho(x,y)=\frac{\operatorname{cov}(x,y)}{\sigma_x\sigma_y}
$$

NumPy で計算してみます。ここでは、距離と余弦類似度は最初の2人の顧客ベクトルで比べ、Pearson 相関係数は 年間購入額_万円 と 来店回数 が全体として同じ方向に動くかを見るために使います。

def minkowski_distance(x, y, p=2):
    return np.sum(np.abs(x - y) ** p) ** (1 / p)


def cosine_similarity(x, y):
    return np.dot(x, y) / (np.linalg.norm(x) * np.linalg.norm(y))


def pearson_correlation(x, y):
    return np.corrcoef(x, y)[0, 1]


customer_a = X_customer_scaled[0]
customer_b = X_customer_scaled[1]
amount = X_customer["年間購入額_万円"].to_numpy()
visits = X_customer["来店回数"].to_numpy()

print("Euclidean distance:", round(minkowski_distance(customer_a, customer_b, p=2), 3))
print("Cosine similarity:", round(cosine_similarity(customer_a, customer_b), 3))
print("Pearson correlation:", round(pearson_correlation(amount, visits), 3))

距離は「小さいほど近い」、類似度や相関は「大きいほど似ている」と読みます。同じ「似ている」でも、何を似ていると見なすかで結果が変わります。

距離の考え方が見えたところで、K-means を小さなコードで実装してみます。

5. NumPy で K-means を実装する

K-means の流れは、NumPy だけでも書けます。

ここでは、2節で作った標準化済みの X_customer_scaled を使い、顧客を3つのクラスタに分けます。

def kmeans_from_scratch(X, n_clusters, max_iter=100, random_state=42):
    rng = np.random.default_rng(random_state)
    initial_indices = rng.choice(len(X), size=n_clusters, replace=False)
    centroids = X[initial_indices].copy()

    for _ in range(max_iter):
        distances = np.linalg.norm(
            X[:, None, :] - centroids[None, :, :],
            axis=2,
        )
        labels = np.argmin(distances, axis=1)

        new_centroids = []
        for cluster in range(n_clusters):
            members = X[labels == cluster]
            if len(members) == 0:
                new_centroids.append(centroids[cluster])
            else:
                new_centroids.append(members.mean(axis=0))

        new_centroids = np.array(new_centroids)
        if np.allclose(centroids, new_centroids):
            break

        centroids = new_centroids

    return labels, centroids


kmeans_labels, kmeans_centroids = kmeans_from_scratch(
    X_customer_scaled,
    n_clusters=3,
)

customer_clustered = customer_jp.assign(cluster=kmeans_labels)
print(customer_clustered.sort_values("cluster"))

クラスタ番号は、グループ名として仮に付いているだけです。cluster=0、cluster=1、cluster=2 の大小には意味がありません。

散布図でも確認します。

plt.scatter(
    X_customer_scaled[:, 0],
    X_customer_scaled[:, 1],
    c=kmeans_labels,
    cmap="tab10",
)
plt.scatter(
    kmeans_centroids[:, 0],
    kmeans_centroids[:, 1],
    marker="x",
    s=200,
    c="black",
    label="centroid",
)
plt.xlabel("年間購入額_万円（標準化後）")
plt.ylabel("来店回数（標準化後）")
plt.title("NumPy で実装した K-means")
plt.legend()
plt.show()

この実装でやっていることは、3節の手順そのものです。

セントロイドを置く
各点から各セントロイドまでの距離を計算する
一番近いセントロイドへ割り当てる
割り当てられた点の平均でセントロイドを更新する

仕組みが見えたところで、実務でよく使う Scikit-learn の K-means も確認します。

6. Scikit-learn で K-means を使う

Scikit-learn では、K-means は KMeans として用意されています。

必要なライブラリが入っていない場合は、先にインストールします。

pip install numpy pandas matplotlib scikit-learn

K-means を使うには、クラスタ数 $K$ をあらかじめ決めておく必要があります。$K$ の決め方に唯一の正解はありませんが、代表的な手がかりとして、クラスタ内平方和とシルエット係数を比べる方法がよく使われます。まずは $K=2$ から $K=6$ まで試して、指標がどう変わるかを見てみます。

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score


ks = range(2, 7)
inertias = []
silhouette_scores = []

for k in ks:
    model = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = model.fit_predict(X_customer_scaled)
    inertias.append(model.inertia_)
    silhouette_scores.append(silhouette_score(X_customer_scaled, labels))

score_table = pd.DataFrame({
    "K": list(ks),
    "クラスタ内平方和": inertias,
    "silhouette": silhouette_scores,
})

print(score_table.round(3))

クラスタ内平方和（inertia_）は、各点から自分のセントロイドまでの距離の2乗をすべて足したものです。小さいほど、各点がセントロイドの近くにまとまっていることを表します。ただし $K$ を増やすほど、クラスタ内平方和は基本的に小さくなっていきます（極端な話、$K$ をサンプル数と同じにすれば、各点が自分だけのクラスタになり、クラスタ内平方和は0になります）。そのため、「$K$ を1つ増やしても、あまり下がらなくなったところ」を探します。折れ線グラフに描いたときの曲がり方が肘（エルボー）に似ていることから、これは エルボー法 と呼ばれます。

シルエット係数は、「同じクラスタの点とはどれくらい近いか」と「一番近い別のクラスタの点とはどれくらい離れているか」を比べる指標です。同じクラスタの点に近く、ほかのクラスタの点からは離れているほど、値は1に近づきます。反対に、クラスタの境目にあるような点では、値は0に近くなります。ざっくり言えば、シルエット係数は「その点は、今のクラスタにちゃんと馴染めているか」を数字にしたものです。

上の表でシルエット係数が高く、クラスタ内平方和の下がり方も落ち着いてくる $K=3$ 付近を選ぶことにします。

model = KMeans(n_clusters=3, n_init=10, random_state=42)
sklearn_labels = model.fit_predict(X_customer_scaled)

customer_clustered = customer_jp.assign(cluster=sklearn_labels)
print(customer_clustered.sort_values("cluster"))

n_init=10 は、初期セントロイドの置き方を変えて10回試し、その中でクラスタ内平方和が一番小さい結果を採用するための指定です。K-means は初期位置の影響を受けやすいアルゴリズムなので、1回だけで決めない方が安定します。random_state=42 は、その初期セントロイドを選ぶ際の乱数のシードです。同じ random_state を指定しておけば、何度実行しても同じ結果を再現できます。42という数字自体に特別な意味はありません。

K-means はシンプルで強力ですが、「丸いまとまり」を探すのが得意な方法です。次は、クラスタを木のように合併していく階層クラスタリングを見ます。

7. 階層クラスタリング

階層クラスタリングは、サンプル同士やクラスタ同士を少しずつまとめて、階層構造を作る方法です。

大きく分けると、次の2種類があります。

方法	考え方
凝集型	最初は各サンプルを1つずつのクラスタとし、近いものから合併していく。
分割型	最初は全サンプルを1つのクラスタとし、少しずつ分割していく。

よく使われるのは、凝集型です。英語では Agglomerative Clustering と呼ばれます。

凝集型の流れは次の通りです。

すべてのサンプルを、それぞれ1つのクラスタと見る。
クラスタ同士の距離を計算する。
一番近い2つのクラスタを合併する。
目的のクラスタ数になるまで繰り返す。

クラスタ同士の距離には、いくつかの定義があります。

名前	意味
single	2つのクラスタ間で、最も近い点同士の距離を見る。
complete	2つのクラスタ間で、最も遠い点同士の距離を見る。
average	2つのクラスタに属する点の組み合わせ全体の平均距離を見る。
ward	合併したときにクラスタ内のばらつきがどれくらい増えるかを見る。

式で書くと、たとえば single linkage は次のようになります。

$$
d_{\min}(C_i,C_j)=\min_{x\in C_i,\ z\in C_j}d(x,z)
$$

complete linkage は、次のようになります。

$$
d_{\max}(C_i,C_j)=\max_{x\in C_i,\ z\in C_j}d(x,z)
$$

K-means が「中心に集める」方法だとすると、階層クラスタリングは「近いものから順に家系図を作る」方法です。

次は、小さな距離行列を使って、どの順番で合併されるかを追ってみます。

8. 階層クラスタリングを距離行列で追う

5個のサンプル $x_1,\ldots,x_5$ があるとします。初期状態では、それぞれが別々のクラスタです。

$$
G_i={x^{(i)}}
$$

距離行列を次のように置きます。

$$
D=
\begin{bmatrix}
0&7&2&9&3\
7&0&5&4&6\
2&5&0&8&1\
9&4&8&0&5\
3&6&1&5&0
\end{bmatrix}
$$

まず、最小の距離は $D_{35}=1$ なので、$G_3$ と $G_5$ を合併します。

$$
G_6={x^{(3)},x^{(5)}}
$$

single linkage では、クラスタ同士の距離を「一番近い点同士の距離」として見ます。そのため、$G_6$ と $G_1$ の距離は、$x_3$ と $x_1$ の距離 2、$x_5$ と $x_1$ の距離 3 のうち小さい方、つまり 2 になります。

次に $G_1$ と $G_6$ が合併し、最後に $G_2$ と $G_4$ が合併すると、2つのクラスタは次のようになります。

$$
C_1={x^{(1)},x^{(3)},x^{(5)}}
$$

$$
C_2={x^{(2)},x^{(4)}}
$$

Scikit-learn で、同じ距離行列を使って確認してみます。

from sklearn.cluster import AgglomerativeClustering


distance_matrix = np.array([
    [0, 7, 2, 9, 3],
    [7, 0, 5, 4, 6],
    [2, 5, 0, 8, 1],
    [9, 4, 8, 0, 5],
    [3, 6, 1, 5, 0],
])

sample_names = ["x1", "x2", "x3", "x4", "x5"]

hierarchy_from_distance = AgglomerativeClustering(
    n_clusters=2,
    metric="precomputed",
    linkage="single",
)

distance_labels = hierarchy_from_distance.fit_predict(distance_matrix)

print(pd.DataFrame({
    "sample": sample_names,
    "cluster": distance_labels,
}))

クラスタ番号は任意ですが、single linkage の考え方では、$x_1,x_3,x_5$ のまとまりと、$x_2,x_4$ のまとまりができます。

距離行列で合併の流れが見えたので、次は顧客データにも階層クラスタリングを当てはめます。

9. Scikit-learn で階層クラスタリングを使う

2節で作った顧客データを、階層クラスタリングでも分けてみます。

from sklearn.cluster import AgglomerativeClustering


hier_model = AgglomerativeClustering(
    n_clusters=3,
    linkage="ward",
)

hier_labels = hier_model.fit_predict(X_customer_scaled)

customer_hier = customer_jp.assign(cluster=hier_labels)
print(customer_hier.sort_values("cluster"))

ward は、クラスタを合併したときに、クラスタ内のばらつき（各点とクラスタの重心との距離の2乗和）がなるべく増えないような組み合わせを選ぶ方法です。まとまりのよい2つのクラスタ同士を合併してもばらつきはあまり増えませんが、離れた2つのクラスタを無理に合併すると、ばらつきは大きく増えます。ward は、合併したときの増え方が一番小さい組み合わせから順に合併していきます。数値特徴量のクラスタリングでよく使われます。

散布図でも確認します。

plt.scatter(
    X_customer_scaled[:, 0],
    X_customer_scaled[:, 1],
    c=hier_labels,
    cmap="tab10",
)
plt.xlabel("年間購入額_万円（標準化後）")
plt.ylabel("来店回数（標準化後）")
plt.title("階層クラスタリング")
plt.show()

階層クラスタリングは、クラスタを作る過程を追いやすいのが長所です。一方、データ数がとても多い場合は計算量が重くなりやすいので、大規模データでは K-means などと使い分けます。

ここまでの K-means や階層クラスタリングは、距離を使ってまとまりを見ていました。次は、密度を使って「点が濃く集まっている場所」を探す DBSCAN を見ます。

10. DBSCAN：密度でクラスタを見つける

DBSCAN は、点が密に集まっている領域 をクラスタとして見つける方法です。

主に2つのパラメータを使います。

パラメータ	意味
`eps`	近所と見なす半径。数式では $\varepsilon$ と書くことが多いです。
`min_samples`	その半径内に何点あれば密な場所と見るか。

サンプル $x$ の $\varepsilon$ 近傍は、次のように書けます。

$$
N_{\varepsilon}(x)={x_i\mid d(x_i,x)\le\varepsilon}
$$

この近傍の中に十分な数のサンプルがあれば、その点は コア点 です。

$$
|N_{\varepsilon}(x)|\ge\operatorname{MinPts}
$$

DBSCAN では、点をおおまかに次の3種類に分けます。

種類	意味
コア点	近くに十分な点がある、クラスタの中心的な点。
境界点	自分自身はコア点ではないが、どこかのコア点の近くにある点。
ノイズ点	どのクラスタにも入りにくい孤立した点。

「コア点から到達できる点をまとめる」というステップを、もう少し具体的に見ておきます。

ある点 $q$ が、コア点 $p$ から距離 $\varepsilon$ 以内にあるとき、$q$ は $p$ から 直接到達可能 であるといいます。

さらに、$p_1=p,\ p_2,\ldots,p_n=q$ のように、隣り合う点同士が直接到達可能な形でコア点をたどっていけるなら、$q$ は $p$ から 到達可能 であるといいます。途中でたどる点はコア点である必要がありますが、最後の $q$ 自身はコア点でなくてもかまいません。

DBSCAN は、あるコア点から到達可能なすべての点（コア点と境界点の両方を含みます）を、同じクラスタにまとめます。境界点がクラスタの「端」になるのは、境界点自身はコア点ではなく、そこからさらに他の点への到達を広げられないためです。どのコア点からも到達可能でない点だけが、ノイズ点として残ります。

K-means はクラスタ数 $K$ を先に決めますが、DBSCAN はクラスタ数を先に決めません。その代わり、eps と min_samples の選び方が結果に大きく効きます。

DBSCAN の考え方がわかったところで、曲がった形のデータで試してみます。

11. Scikit-learn で DBSCAN を使う

K-means は、丸いまとまりを探すのが得意です。反対に、三日月のように曲がったクラスタは苦手です。

そこで、make_moons で曲がった2つのクラスタを作り、DBSCAN を当ててみます。

from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.preprocessing import StandardScaler


X_moons, _ = make_moons(n_samples=300, noise=0.06, random_state=42)
X_moons = StandardScaler().fit_transform(X_moons)

dbscan_model = DBSCAN(eps=0.25, min_samples=5)
dbscan_labels = dbscan_model.fit_predict(X_moons)

print(pd.Series(dbscan_labels).value_counts().sort_index())

DBSCAN では、ラベル -1 がノイズ点を表します。

散布図で見ます。

plt.scatter(
    X_moons[:, 0],
    X_moons[:, 1],
    c=dbscan_labels,
    cmap="tab10",
    s=20,
)
plt.xlabel("feature 1")
plt.ylabel("feature 2")
plt.title("DBSCAN による密度クラスタリング")
plt.show()

DBSCAN は、形が丸くないクラスタにも対応しやすく、ノイズ点も分けられるのが長所です。ただし、密度が場所によって大きく変わるデータでは、1つの eps で全体をうまく扱うのが難しくなることがあります。

ここまでの方法は、点をはっきりどこかのクラスタへ割り当てていました。次は、「この点はクラスタAに70%、クラスタBに30%くらい属していそう」と確率で見るガウス混合モデルに進みます。

12. ガウス混合モデル：確率分布としてクラスタを見る

ガウス混合モデルは、データが複数のガウス分布から混ざって生まれたと考える方法です。英語では Gaussian Mixture Model、略して GMM と呼ばれます。

まず、1次元のガウス分布は次の形です。

$$
f(x)=\frac{1}{\sqrt{2\pi}\sigma}
\exp\left[-\frac{(x-\mu)^2}{2\sigma^2}\right]
$$

多次元では、平均ベクトル $\mu$ と共分散行列 $\Sigma$ を使います。

$$
f(x)=\frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}}
\exp\left[-\frac12(x-\mu)^T\Sigma^{-1}(x-\mu)\right]
$$

1つのガウス分布だけでは、山が1つのデータしか表しにくいです。そこで、複数のガウス分布を重み付きで足し合わせます。

$$
p(x)=\sum_{k=1}^{K}\pi_k\mathcal{N}(x\mid\mu_k,\Sigma_k)
$$

ここで、混合比 $\pi_k$ は次の条件を満たします。

$$
\pi_k\ge0,\qquad \sum_{k=1}^{K}\pi_k=1
$$

K-means が「一番近い中心を選ぶ」方法だとすると、GMM は「どの山から生まれた可能性が高いか」を確率で見る方法です。

この確率をどう推定するかを見るために、次は EM アルゴリズムを確認します。

13. EM アルゴリズムの考え方

GMM では、どのサンプルがどのガウス成分から生まれたかは見えません。そこで、見えない割り当てを確率として推定しながら、分布のパラメータも更新します。

この反復手順が EM アルゴリズム です。

E ステップでは、サンプル $x_i$ が第 $k$ 成分に属する確率を計算します。この確率は 責任度（負担率、英語では responsibility）と呼ばれます。

$$
\gamma(i,k)=
\frac{\pi_k\mathcal{N}(x_i\mid\mu_k,\Sigma_k)}
{\sum_{j=1}^{K}\pi_j\mathcal{N}(x_i\mid\mu_j,\Sigma_j)}
$$

式の形は、ベイズの定理に近いものです。分子は「第 $k$ 成分の重み $\times$ 第 $k$ 成分がそのあたりでどれくらい山になっているか」、分母はすべての成分について同じ計算をした合計です。つまり $\gamma(i,k)$ は、「もし $x_i$ がどれかの成分から生まれたのだとしたら、それは第 $k$ 成分だった可能性がどれくらいか」を、ほかの成分と比べたうえでの割合として表したものです。$k$ について足し合わせると、必ず1になります。

M ステップでは、その責任度を重みとして、各成分の平均、共分散、混合比を更新します。

$$
N_k=\sum_{i=1}^{N}\gamma(i,k)
$$

$$
\mu_k=\frac{1}{N_k}\sum_{i=1}^{N}\gamma(i,k)x_i
$$

$$
\Sigma_k=
\frac{1}{N_k}\sum_{i=1}^{N}
\gamma(i,k)(x_i-\mu_k)(x_i-\mu_k)^T
$$

$$
\pi_k=\frac{N_k}{N}
$$

流れとしては、次のようになります。

ファインマン風に言えば、EM は「仮にこの山の担当だろうと重みを付ける、その重みに合わせて山の位置と形を直す」を繰り返す方法です。

では、Scikit-learn で GMM を動かしてみます。

14. Scikit-learn でガウス混合モデルを使う

Scikit-learn では、GMM は GaussianMixture として用意されています。

ここでは、make_blobs で3つのまとまりを持つデータを作り、GMM を当てます。

from sklearn.datasets import make_blobs
from sklearn.mixture import GaussianMixture


X_gmm, _ = make_blobs(
    n_samples=400,
    centers=3,
    cluster_std=[0.6, 1.0, 0.8],
    random_state=42,
)

gmm_model = GaussianMixture(
    n_components=3,
    covariance_type="full",
    random_state=42,
)

gmm_labels = gmm_model.fit_predict(X_gmm)
responsibility = gmm_model.predict_proba(X_gmm[:5])

print(pd.DataFrame(responsibility).round(3))

predict_proba で出てくる値が、各サンプルが各成分に属する確率です。K-means のように1つのクラスタへきっぱり割り当てるだけでなく、迷い具合も見られます。

散布図で確認します。

plt.scatter(
    X_gmm[:, 0],
    X_gmm[:, 1],
    c=gmm_labels,
    cmap="tab10",
    s=20,
)
plt.scatter(
    gmm_model.means_[:, 0],
    gmm_model.means_[:, 1],
    marker="x",
    s=200,
    c="black",
    label="mean",
)
plt.xlabel("feature 1")
plt.ylabel("feature 2")
plt.title("GaussianMixture によるクラスタリング")
plt.legend()
plt.show()

GMM は、クラスタを楕円のような分布として表せるのが長所です。K-means より柔らかく、「どのクラスタにどれくらい属しているか」を見たいときに向いています。

最後に、ここまでの流れを、少し実務の分析に近い形でまとめます。

15. ケース：アジアサッカーチーム風の指標をクラスタリングする

最後に、アジアサッカーチームの分析に近い形で、チームの特徴をクラスタリングする流れを見ます。

ここでは、説明用に作った架空のチーム指標データを使います。実在チームの最新成績ではありません。目的は、数値特徴量を標準化し、クラスタリングし、クラスタごとの特徴を読む一連の流れを確認することです。

from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler


team_jp = pd.DataFrame({
    "チーム": [
        "チームA", "チームB", "チームC", "チームD", "チームE", "チームF",
        "チームG", "チームH", "チームI", "チームJ", "チームK", "チームL",
    ],
    "得点力": [82, 85, 80, 84, 50, 48, 52, 46, 35, 38, 33, 40],
    "守備安定度": [76, 79, 74, 78, 68, 70, 66, 71, 42, 45, 40, 48],
    "平均支配率": [56, 58, 55, 57, 42, 44, 40, 43, 34, 36, 32, 37],
    "国際経験": [88, 91, 85, 90, 55, 52, 58, 50, 28, 30, 25, 32],
})

team_features = ["得点力", "守備安定度", "平均支配率", "国際経験"]

各列の意味は次の通りです。

列名	意味
`チーム`	説明用に用意した仮のチーム名です（実在のチームではありません）。
`得点力`	攻撃面の得点力を表す指標です。0から100のスコアで、大きいほど高評価とします。
`守備安定度`	守備の安定度を表す指標です。0から100のスコアで、大きいほど高評価とします。
`平均支配率`	試合におけるボール支配率の平均です。単位は%です。
`国際経験`	国際大会での経験の豊富さを表す指標です。0から100のスコアで、大きいほど高評価とします。

team_scaled = StandardScaler().fit_transform(team_jp[team_features])

team_model = KMeans(n_clusters=3, n_init=10, random_state=42)
team_jp["cluster"] = team_model.fit_predict(team_scaled)

print(team_jp.sort_values("cluster"))

クラスタごとの平均を見ると、どのようなタイプに分かれたかを読みやすくなります。

cluster_summary = team_jp.groupby("cluster")[team_features].mean()
print(cluster_summary.round(1))

たとえば、得点力・守備安定度・平均支配率・国際経験 のすべてが高いクラスタは、実力が高いチーム群と読めます。すべてが中位のクラスタは中堅どころ、すべてがまだ低いクラスタは発展途上のチーム群と読めます。今回のデータでは4つの指標がおおむね同じ方向に動くように作っているため、このようにクラスタ全体の水準の高さで読み分けられます。実際のデータでは、指標ごとに違う方向へ動くこともあるので、その場合は「攻撃は強いが守備は不安定」のように、クラスタごとに違う持ち味を見つけることになります。

クラスタリングの結果は、出しただけではまだ分析になりません。分かれたあとの平均値、分布、元データの行を見ながら、「このクラスタは何を表していそうか」を人間が解釈する必要があります。

16. まとめ

今回の内容をまとめます。

クラスタリングは、正解ラベルのないデータを似ているもの同士に分ける教師なし学習である。
クラスタ番号そのものには大小や善悪の意味はなく、分けたあとに特徴を読むことが大切である。
K-means は、セントロイドに近い点を割り当て、平均でセントロイドを更新する方法である。
距離や類似度の選び方によって、「似ている」の意味は変わる。
K-means では、特徴量のスケールをそろえることが重要である。
階層クラスタリングは、近いクラスタを順に合併して階層構造を作る。
DBSCAN は、密度の高い場所をクラスタとして見つけ、ノイズ点も扱える。
ガウス混合モデルは、データが複数の確率分布から生まれたと考え、所属確率も出せる。
クラスタリングでは、結果を出すだけでなく、クラスタごとの特徴を読み取るところまでが分析である。

クラスタリングは、分類や回帰と違って、最初から正解が与えられていません。そのぶん、結果をどう解釈するかがとても重要です。

まずは「距離で近いものを集める K-means」「階層としてまとめる階層クラスタリング」「密度で見つける DBSCAN」「確率分布で見る GMM」という4つの見方を押さえておくと、教師なし学習の入口がかなり見通しやすくなります。