DBSCAN密度ベースクラスタリングの理論と実装

都市の地図を見て「繁華街」を見つけるとき、私たちは建物が密集している地域を自然に見分けています。建物がまばらな地域は「郊外」、完全に孤立した建物は「外れ値」です。この「密集度」に基づいてグループを見つける考え方を機械学習に持ち込んだのがDBSCAN(Density-Based Spatial Clustering of Applications with Noise)です。

k-meansのような中心ベースのクラスタリングは、クラスタが球状で均等なサイズであることを暗黙に仮定しています。しかし現実のデータでは、三日月形やリング状のクラスタ、密度の異なるクラスタ、ノイズ(外れ値)が混在することが少なくありません。DBSCANはこれらの問題をすべて解決できる、実用上非常に重要なアルゴリズムです。

DBSCANを理解すると、以下のような場面で活用できます。

  • 地理データ分析: GPSデータからの滞在地・移動パターンの検出
  • 異常検知: ノイズ点の自動検出による外れ値の発見
  • 画像処理: 任意形状の領域分割
  • 科学データ分析: 粒子衝突データや天体観測データのクラスタ検出

本記事の内容

  • DBSCANの基本概念(ε-近傍、コア点、到達可能性)
  • アルゴリズムの詳細と計算量
  • パラメータ選択の方法(k-距離グラフ)
  • k-meansとの比較実験
  • HDBSCANへの拡張
  • Pythonでのスクラッチ実装

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

DBSCANの基本概念

密度ベースクラスタリングの直感

DBSCANの中心的なアイデアは「データ点が密集している領域を1つのクラスタとし、密度の低い領域でクラスタを分離する」というものです。

この密度を定義するために、DBSCANは2つのパラメータを使います。

  • ε(イプシロン): 近傍の半径。各点の周囲にこの半径の「円」(高次元では超球)を描く
  • MinPts: コア点と見なすために必要な近傍内の最小点数

この2つのパラメータにより、データ中の各点は3種類に分類されます。

3種類の点

コア点(Core Point): ε-近傍内に自分自身を含めてMinPts個以上の点を持つ点。クラスタの「内部」を構成する点です。

数式で表すと、点 $\bm{p}$ がコア点であるとは

$$ |N_\epsilon(\bm{p})| \geq \text{MinPts} $$

を満たすことです。ここで $N_\epsilon(\bm{p}) = \{\bm{q} \in \mathcal{D} : d(\bm{p}, \bm{q}) \leq \epsilon\}$ は点 $\bm{p}$ のε-近傍(距離 $\epsilon$ 以内の全点の集合)です。

ボーダー点(Border Point): 自身はコア点ではないが、あるコア点のε-近傍に含まれる点。クラスタの「境界」に位置する点です。

$$ |N_\epsilon(\bm{p})| < \text{MinPts} \quad \text{かつ} \quad \exists \bm{q} \in \mathcal{D} : \bm{p} \in N_\epsilon(\bm{q}) \text{ and } |N_\epsilon(\bm{q})| \geq \text{MinPts} $$

ノイズ点(Noise Point): コア点でもボーダー点でもない点。どのクラスタにも属さない外れ値です。

密度到達可能性と密度接続

クラスタを形成するために、DBSCANは2つの概念を定義します。

直接密度到達可能(Directly Density-Reachable): 点 $\bm{q}$ が点 $\bm{p}$ から直接密度到達可能であるとは、$\bm{p}$ がコア点で、かつ $\bm{q} \in N_\epsilon(\bm{p})$ であることです。

この関係は非対称であることに注意してください。$\bm{p}$ がコア点で $\bm{q}$ がボーダー点の場合、$\bm{q}$ は $\bm{p}$ から直接密度到達可能ですが、$\bm{p}$ は $\bm{q}$ から直接密度到達可能とは限りません($\bm{q}$ はコア点ではないため)。

密度到達可能(Density-Reachable): 点 $\bm{q}$ が点 $\bm{p}$ から密度到達可能であるとは、点の列 $\bm{p}_1 = \bm{p}, \bm{p}_2, \dots, \bm{p}_k = \bm{q}$ が存在して、各 $\bm{p}_{i+1}$ が $\bm{p}_i$ から直接密度到達可能であることです。つまり、コア点を辿る「チェーン」で接続できるということです。

$$ \bm{p}_1 \xrightarrow{\text{直接}} \bm{p}_2 \xrightarrow{\text{直接}} \cdots \xrightarrow{\text{直接}} \bm{p}_k $$

密度接続(Density-Connected): 2つの点 $\bm{p}$ と $\bm{q}$ が密度接続されているとは、ある点 $\bm{o}$ が存在して、$\bm{p}$ も $\bm{q}$ も $\bm{o}$ から密度到達可能であることです。

$$ \exists \bm{o} : \bm{p} \xleftarrow{\text{到達可能}} \bm{o} \xrightarrow{\text{到達可能}} \bm{q} $$

密度接続は対称な関係です。DBSCANにおけるクラスタは、互いに密度接続された点の極大集合として定義されます。

クラスタの形式的定義

DBSCANのクラスタ $C$ は、以下の2つの条件を満たす空でない部分集合です。

最大性: $\bm{p} \in C$ かつ $\bm{q}$ が $\bm{p}$ から密度到達可能ならば、$\bm{q} \in C$

接続性: 任意の $\bm{p}, \bm{q} \in C$ について、$\bm{p}$ と $\bm{q}$ は密度接続されている

この形式的定義により、DBSCANのクラスタは任意の形状を取ることができます。k-meansのクラスタがボロノイ分割(各点が最寄りの中心に属する)に制約されるのとは対照的です。

これらの概念を踏まえて、次にDBSCANのアルゴリズムの具体的な手順を見ていきましょう。

DBSCANのアルゴリズム

手順

DBSCANのアルゴリズムは驚くほどシンプルです。

入力: データセット $\mathcal{D}$、ε(近傍半径)、MinPts(最小点数)

  1. すべての点を「未訪問」とマークする
  2. 各未訪問点 $\bm{p}$ について: – $\bm{p}$ を「訪問済み」にマーク – $N_\epsilon(\bm{p})$ を計算(ε-近傍の全点を取得) – $|N_\epsilon(\bm{p})| < \text{MinPts}$ ならば、$\bm{p}$ をノイズとしてマーク(後でボーダー点に変更される可能性あり) - $|N_\epsilon(\bm{p})| \geq \text{MinPts}$ ならば:
    • 新しいクラスタ $C$ を作成し、$\bm{p}$ を $C$ に追加
    • $N_\epsilon(\bm{p})$ 内の全点を「シード集合」$S$ に追加
    • $S$ の各点 $\bm{q}$ について:
    • $\bm{q}$ が未訪問ならば、「訪問済み」にマークし $N_\epsilon(\bm{q})$ を計算
    • $|N_\epsilon(\bm{q})| \geq \text{MinPts}$ ならば、$N_\epsilon(\bm{q})$ を $S$ に追加
    • $\bm{q}$ がまだどのクラスタにも属さないならば、$\bm{q}$ を $C$ に追加

計算量

ε-近傍の探索が全データ点に対して行われるため、ナイーブな実装では

$$ O(n^2) $$

の計算量がかかります。しかし、空間インデックス(kd-tree や ball-tree)を使用すると、近傍探索が $O(\log n)$ で行えるため、全体の計算量は

$$ O(n \log n) $$

に改善されます。ただし、高次元データでは空間インデックスの効率が低下し(次元の呪い)、$O(n^2)$ に近づく場合があります。

k-meansとの比較

DBSCANとk-meansの根本的な違いを整理しましょう。

特性 DBSCAN k-means
クラスタ形状 任意の形状 球状(凸型)
クラスタ数 自動決定 事前指定が必要
ノイズ処理 ノイズ点を自動検出 全点をいずれかのクラスタに割当
パラメータ ε, MinPts k(クラスタ数)
密度の違い 均一密度を仮定 均一密度を仮定
決定論的 はい(ボーダー点を除く) いいえ(初期値依存)
計算量 $O(n \log n)$〜$O(n^2)$ $O(nk \cdot \text{iter})$

ここで重要な注意点があります。DBSCANも「均一密度」を暗黙に仮定しています。ε と MinPts がデータ全体で固定されているため、密度の大きく異なるクラスタが共存する場合、適切なパラメータの設定が困難になります。この問題はHDBSCANで解決されますが、後のセクションで詳しく解説します。

パラメータの選び方を知る前に、まずDBSCANをPythonでスクラッチ実装してアルゴリズムの動きを確認しましょう。

Pythonによるスクラッチ実装

DBSCANの実装

まず、DBSCANを一から実装します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

class DBSCANScratch:
    """DBSCANのスクラッチ実装"""

    NOISE = -1
    UNVISITED = -2

    def __init__(self, eps=0.5, min_pts=5):
        self.eps = eps
        self.min_pts = min_pts
        self.labels_ = None

    def fit(self, X):
        n_samples = X.shape[0]
        self.labels_ = np.full(n_samples, self.UNVISITED)

        # 全点間の距離行列を事前計算(小規模データ向け)
        dist_matrix = cdist(X, X, metric="euclidean")

        cluster_id = 0

        for i in range(n_samples):
            if self.labels_[i] != self.UNVISITED:
                continue

            # ε-近傍を取得
            neighbors = self._get_neighbors(dist_matrix, i)

            if len(neighbors) < self.min_pts:
                # コア点の条件を満たさない → 暫定的にノイズ
                self.labels_[i] = self.NOISE
            else:
                # 新しいクラスタを形成
                self._expand_cluster(dist_matrix, i, neighbors, cluster_id)
                cluster_id += 1

        self.n_clusters_ = cluster_id
        self.n_noise_ = np.sum(self.labels_ == self.NOISE)
        return self

    def _get_neighbors(self, dist_matrix, point_idx):
        """点のε-近傍を取得"""
        return np.where(dist_matrix[point_idx] <= self.eps)[0]

    def _expand_cluster(self, dist_matrix, point_idx, neighbors, cluster_id):
        """コア点からクラスタを拡張"""
        self.labels_[point_idx] = cluster_id

        # シード集合(探索待ちの近傍点)
        seed_set = list(neighbors)
        idx = 0

        while idx < len(seed_set):
            q = seed_set[idx]

            if self.labels_[q] == self.NOISE:
                # ノイズ → ボーダー点としてクラスタに追加
                self.labels_[q] = cluster_id
            elif self.labels_[q] == self.UNVISITED:
                # 未訪問 → クラスタに追加
                self.labels_[q] = cluster_id

                # qの近傍を確認
                q_neighbors = self._get_neighbors(dist_matrix, q)
                if len(q_neighbors) >= self.min_pts:
                    # qもコア点 → 近傍をシード集合に追加
                    for n in q_neighbors:
                        if n not in seed_set:
                            seed_set.append(n)

            idx += 1

    def fit_predict(self, X):
        self.fit(X)
        return self.labels_


# テストデータ: 三日月形データ
from sklearn.datasets import make_moons, make_blobs
np.random.seed(42)

X_moons, y_moons = make_moons(n_samples=300, noise=0.08, random_state=42)

# DBSCANの実行
dbscan = DBSCANScratch(eps=0.2, min_pts=5)
labels = dbscan.fit_predict(X_moons)

print(f"検出されたクラスタ数: {dbscan.n_clusters_}")
print(f"ノイズ点の数: {dbscan.n_noise_}")

# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# DBSCANの結果
colors = plt.cm.tab10(np.linspace(0, 1, 10))
for cluster_id in range(dbscan.n_clusters_):
    mask = labels == cluster_id
    axes[0].scatter(X_moons[mask, 0], X_moons[mask, 1],
                    color=colors[cluster_id], label=f"クラスタ {cluster_id}",
                    s=30, edgecolors="k", linewidths=0.5)
noise_mask = labels == -1
if noise_mask.sum() > 0:
    axes[0].scatter(X_moons[noise_mask, 0], X_moons[noise_mask, 1],
                    color="gray", marker="x", s=30, label="ノイズ")
axes[0].set_title("DBSCAN (ε=0.2, MinPts=5)", fontsize=13)
axes[0].legend(fontsize=10)
axes[0].set_xlabel("$x_1$")
axes[0].set_ylabel("$x_2$")

# 真のラベル
axes[1].scatter(X_moons[:, 0], X_moons[:, 1], c=y_moons, cmap="tab10",
                s=30, edgecolors="k", linewidths=0.5)
axes[1].set_title("真のラベル", fontsize=13)
axes[1].set_xlabel("$x_1$")
axes[1].set_ylabel("$x_2$")

plt.tight_layout()
plt.savefig("dbscan_moons.png", dpi=150, bbox_inches="tight")
plt.show()

三日月形のデータに対するDBSCANの結果を見ると、2つの三日月を正しく分離できていることが分かります。k-meansではこのような非凸形状のクラスタを分離することができません(中心点ベースのため、直線的な分離しかできない)。また、ノイズ点が自動的に検出されている点も重要です。

k-meansとの比較実験

DBSCANとk-meansの違いを、さまざまな形状のデータで比較してみましょう。

from sklearn.cluster import KMeans
from sklearn.datasets import make_circles

# 3種類のテストデータ
datasets = []

# 1. 三日月形
X1, y1 = make_moons(n_samples=300, noise=0.08, random_state=42)
datasets.append(("三日月形", X1, y1, 0.2, 5))

# 2. 同心円
X2, y2 = make_circles(n_samples=300, noise=0.05, factor=0.5, random_state=42)
datasets.append(("同心円", X2, y2, 0.15, 5))

# 3. 不均等な密度のブロブ + ノイズ
X3_blobs, y3_blobs = make_blobs(
    n_samples=[100, 100, 50], centers=[[-3, 0], [3, 0], [0, 4]],
    cluster_std=[0.5, 0.5, 0.3], random_state=42
)
noise = np.random.RandomState(42).uniform(-5, 5, (30, 2))
X3 = np.vstack([X3_blobs, noise])
y3 = np.hstack([y3_blobs, np.full(30, -1)])
datasets.append(("ブロブ+ノイズ", X3, y3, 0.8, 5))

fig, axes = plt.subplots(3, 3, figsize=(15, 15))

for row, (name, X, y_true, eps, min_pts) in enumerate(datasets):
    # 真のラベル
    axes[row, 0].scatter(X[:, 0], X[:, 1], c=y_true, cmap="tab10",
                         s=20, edgecolors="k", linewidths=0.3)
    axes[row, 0].set_title(f"{name}: 真のラベル", fontsize=12)

    # k-means
    n_clusters_true = len(set(y_true)) - (1 if -1 in y_true else 0)
    kmeans = KMeans(n_clusters=n_clusters_true, random_state=42, n_init=10)
    km_labels = kmeans.fit_predict(X)
    axes[row, 1].scatter(X[:, 0], X[:, 1], c=km_labels, cmap="tab10",
                         s=20, edgecolors="k", linewidths=0.3)
    axes[row, 1].set_title(f"k-means (k={n_clusters_true})", fontsize=12)

    # DBSCAN
    dbscan = DBSCANScratch(eps=eps, min_pts=min_pts)
    db_labels = dbscan.fit_predict(X)
    scatter_colors = np.where(db_labels == -1, -1, db_labels)
    axes[row, 2].scatter(X[:, 0], X[:, 1], c=scatter_colors, cmap="tab10",
                         s=20, edgecolors="k", linewidths=0.3)
    noise_count = np.sum(db_labels == -1)
    axes[row, 2].set_title(
        f"DBSCAN (ε={eps}, MinPts={min_pts}, ノイズ={noise_count}点)",
        fontsize=11
    )

for ax in axes.flat:
    ax.set_xlabel("$x_1$", fontsize=10)
    ax.set_ylabel("$x_2$", fontsize=10)

plt.suptitle("k-means vs DBSCAN: 形状への適応力", fontsize=15, y=1.01)
plt.tight_layout()
plt.savefig("kmeans_vs_dbscan.png", dpi=150, bbox_inches="tight")
plt.show()

この比較実験から、DBSCANの強みが明確に分かります。三日月形データではk-meansが中央で直線的に分割してしまうのに対し、DBSCANは曲線的なクラスタを正しく検出します。同心円データでもk-meansは内円と外円を分離できませんが、DBSCANは密度の違いを捉えて正確に分離します。ノイズを含むデータでは、k-meansが全点をいずれかのクラスタに無理に割り当てるのに対し、DBSCANはノイズ点を自動的に検出・分離します。

パラメータ選択

εの選び方: k-距離グラフ

DBSCANの性能はε と MinPts の選択に大きく依存します。特にεの選択が重要で、小さすぎるとほぼ全点がノイズになり、大きすぎると全データが1つのクラスタに統合されてしまいます。

εの選択にはk-距離グラフ(k-distance graph)が有効です。各点について、k番目に近い点までの距離(k-距離)を計算し、降順にソートしてプロットします。このグラフの「肘」(エルボー)の位置が適切なεの値を示唆します。

ここで $k = \text{MinPts} – 1$ とするのが一般的です(自分自身を含めてMinPts個なので、k番目の近傍は自分以外でMinPts – 1番目)。

from sklearn.neighbors import NearestNeighbors

def plot_k_distance(X, k=4, ax=None):
    """k-距離グラフを描画"""
    nn = NearestNeighbors(n_neighbors=k + 1)  # 自分自身を含む
    nn.fit(X)
    distances, _ = nn.kneighbors(X)

    # k番目の近傍までの距離(自分自身を除く)
    k_distances = distances[:, k]
    k_distances = np.sort(k_distances)[::-1]  # 降順ソート

    if ax is None:
        fig, ax = plt.subplots(figsize=(8, 5))

    ax.plot(range(len(k_distances)), k_distances, linewidth=2)
    ax.set_xlabel("データ点(距離の降順)", fontsize=12)
    ax.set_ylabel(f"{k}-距離", fontsize=12)
    ax.set_title(f"k-距離グラフ (k={k})", fontsize=14)
    ax.grid(True, alpha=0.3)

    return k_distances


# 三日月形データでのk-距離グラフ
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

k_dist = plot_k_distance(X_moons, k=4, ax=axes[0])
axes[0].axhline(y=0.2, color="red", linestyle="--", label="ε = 0.2")
axes[0].legend(fontsize=11)

# 異なるεでのDBSCAN結果
eps_values = [0.1, 0.2, 0.3, 0.5]
results_text = []
for eps in eps_values:
    db = DBSCANScratch(eps=eps, min_pts=5)
    db.fit(X_moons)
    results_text.append(f"ε={eps}: {db.n_clusters_}クラスタ, {db.n_noise_}ノイズ")

axes[1].text(0.1, 0.5, "\n".join(results_text), fontsize=13,
             transform=axes[1].transAxes, verticalalignment="center",
             family="monospace",
             bbox=dict(boxstyle="round", facecolor="lightyellow", alpha=0.8))
axes[1].set_title("ε値によるクラスタリング結果の変化", fontsize=14)
axes[1].axis("off")

plt.tight_layout()
plt.savefig("k_distance_graph.png", dpi=150, bbox_inches="tight")
plt.show()

k-距離グラフでは、クラスタ内の点は近傍までの距離が小さく、ノイズ点は距離が大きくなります。グラフの急な変化点(エルボー)が、クラスタ内の点とノイズ点の境界を示しています。この例ではε = 0.2付近にエルボーがあり、実際にε = 0.2で良好なクラスタリング結果が得られています。

MinPtsの選び方

MinPtsについては、以下の経験則が知られています。

  • 一般的な目安: $\text{MinPts} \geq d + 1$($d$ は次元数)。2次元データなら MinPts ≥ 3
  • 推奨値: $\text{MinPts} = 2d$(Ester et al., 1996の原著論文)。2次元なら MinPts = 4
  • ノイズが多い場合: MinPtsを大きくすることでロバスト性が向上。ただし大きすぎると小さなクラスタを見逃す
# MinPtsの影響を可視化
fig, axes = plt.subplots(1, 4, figsize=(20, 4))
min_pts_values = [3, 5, 10, 20]

for ax, min_pts in zip(axes, min_pts_values):
    db = DBSCANScratch(eps=0.2, min_pts=min_pts)
    labels = db.fit_predict(X_moons)

    for cid in range(db.n_clusters_):
        mask = labels == cid
        ax.scatter(X_moons[mask, 0], X_moons[mask, 1], s=20,
                   edgecolors="k", linewidths=0.3)
    noise = labels == -1
    if noise.sum() > 0:
        ax.scatter(X_moons[noise, 0], X_moons[noise, 1],
                   color="gray", marker="x", s=20)
    ax.set_title(f"MinPts={min_pts}\n{db.n_clusters_}クラスタ, {db.n_noise_}ノイズ",
                 fontsize=11)
    ax.set_xlabel("$x_1$")
    ax.set_ylabel("$x_2$")

plt.suptitle("MinPtsの影響 (ε=0.2固定)", fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("minpts_effect.png", dpi=150, bbox_inches="tight")
plt.show()

MinPtsを変化させた結果を見ると、MinPtsが小さすぎると(例: 3)ノイズが少なくなる反面、ノイズ点がクラスタに取り込まれやすくなります。MinPtsが大きすぎると(例: 20)多くの点がノイズと判定され、クラスタが断片化する可能性があります。このデータでは MinPts = 5 が適切なバランスを示しています。

パラメータ選択の方法を理解したところで、次にDBSCANの限界と、それを克服するHDBSCANについて見ていきましょう。

DBSCANの限界

密度の異なるクラスタの問題

DBSCANの最大の弱点は、εがグローバルに固定されているため、密度の異なるクラスタを同時に検出するのが困難なことです。

# 密度の異なるクラスタ
np.random.seed(42)

# 高密度クラスタ
X_dense = np.random.randn(200, 2) * 0.3 + np.array([-2, 0])

# 低密度クラスタ
X_sparse = np.random.randn(100, 2) * 1.5 + np.array([4, 0])

# ノイズ
X_noise = np.random.uniform(-6, 8, (30, 2))

X_mixed = np.vstack([X_dense, X_sparse, X_noise])
y_mixed = np.hstack([np.zeros(200), np.ones(100), np.full(30, -1)])

# 異なるεでの結果
fig, axes = plt.subplots(1, 4, figsize=(20, 4))

# 真のラベル
axes[0].scatter(X_mixed[:200, 0], X_mixed[:200, 1], c="blue", s=15, label="高密度")
axes[0].scatter(X_mixed[200:300, 0], X_mixed[200:300, 1], c="red", s=15, label="低密度")
axes[0].scatter(X_mixed[300:, 0], X_mixed[300:, 1], c="gray", marker="x", s=15, label="ノイズ")
axes[0].set_title("真のラベル", fontsize=12)
axes[0].legend(fontsize=9)

for ax, eps in zip(axes[1:], [0.3, 0.8, 1.5]):
    db = DBSCANScratch(eps=eps, min_pts=5)
    labels = db.fit_predict(X_mixed)

    for cid in range(db.n_clusters_):
        mask = labels == cid
        ax.scatter(X_mixed[mask, 0], X_mixed[mask, 1], s=15,
                   edgecolors="k", linewidths=0.2)
    noise = labels == -1
    if noise.sum() > 0:
        ax.scatter(X_mixed[noise, 0], X_mixed[noise, 1],
                   color="gray", marker="x", s=15)
    ax.set_title(f"ε={eps}: {db.n_clusters_}クラスタ, {db.n_noise_}ノイズ", fontsize=11)

for ax in axes:
    ax.set_xlabel("$x_1$")
    ax.set_ylabel("$x_2$")

plt.suptitle("DBSCANの限界: 密度の異なるクラスタ", fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("dbscan_limitation.png", dpi=150, bbox_inches="tight")
plt.show()

実験結果を見ると、ε = 0.3 では高密度クラスタは正しく検出されますが、低密度クラスタは細かく分割されるかノイズ扱いになります。ε = 1.5 では低密度クラスタは検出されますが、高密度クラスタとノイズが統合されてしまいます。ε = 0.8 は中間的ですが、どちらのクラスタも完全には正しく検出できません。これがDBSCANの本質的な限界です。

高次元データでの課題

高次元空間では「次元の呪い」により、全ての点間の距離が似通ってきます。そのため、εの設定が極めて困難になります。一般的に、DBSCANは20次元程度までのデータに適しており、それ以上の高次元データでは次元削減を前処理として行うのが望ましいです。

HDBSCANへの拡張

DBSCANの密度パラメータεの問題を解決するのがHDBSCAN(Hierarchical DBSCAN)です。HDBSCANはεを固定する代わりに、すべての可能なε値に対するDBSCANの結果を階層的に統合します。

HDBSCANの直感

HDBSCANの中心的なアイデアは次のとおりです。

  1. εを $\infty$ から $0$ まで徐々に小さくしていくことを想像する
  2. ε が大きいときは全データが1つのクラスタ
  3. εを小さくするにつれて、密度の低い領域でクラスタが分裂していく
  4. この分裂過程を木構造(樹形図)で表現する
  5. 木構造から「最も安定した」クラスタを抽出する

相互到達距離

HDBSCANは相互到達距離(mutual reachability distance)を導入します。点 $\bm{p}$ のコア距離 $d_{\text{core}}(\bm{p})$ を、$k$ 番目に近い点までの距離とします($k = \text{MinPts}$)。

$$ d_{\text{core}_k}(\bm{p}) = d(\bm{p}, \text{NN}_k(\bm{p})) $$

相互到達距離は

$$ d_{\text{mreach}}(\bm{p}, \bm{q}) = \max\{d_{\text{core}}(\bm{p}),\ d_{\text{core}}(\bm{q}),\ d(\bm{p}, \bm{q})\} $$

と定義されます。これにより、密集領域にある点同士は実際の距離で評価されますが、疎な領域にある点は「膨らまされた」距離で評価されるため、密度の違いに適応できるようになります。

scikit-learnでのHDBSCAN

HDBSCANはscikit-learn 1.3以降で利用可能です。

from sklearn.cluster import HDBSCAN

# 密度の異なるデータに対してHDBSCANを適用
hdbscan = HDBSCAN(min_cluster_size=10, min_samples=5)
hdbscan_labels = hdbscan.fit_predict(X_mixed)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 真のラベル
axes[0].scatter(X_mixed[:200, 0], X_mixed[:200, 1], c="blue", s=15, label="高密度")
axes[0].scatter(X_mixed[200:300, 0], X_mixed[200:300, 1], c="red", s=15, label="低密度")
axes[0].scatter(X_mixed[300:, 0], X_mixed[300:, 1], c="gray", marker="x", s=15, label="ノイズ")
axes[0].set_title("真のラベル", fontsize=13)
axes[0].legend(fontsize=10)

# DBSCAN (最良のε)
db_best = DBSCANScratch(eps=0.8, min_pts=5)
db_labels = db_best.fit_predict(X_mixed)
for cid in range(db_best.n_clusters_):
    mask = db_labels == cid
    axes[1].scatter(X_mixed[mask, 0], X_mixed[mask, 1], s=15,
                    edgecolors="k", linewidths=0.2)
noise = db_labels == -1
if noise.sum() > 0:
    axes[1].scatter(X_mixed[noise, 0], X_mixed[noise, 1],
                    c="gray", marker="x", s=15)
axes[1].set_title(f"DBSCAN (ε=0.8): {db_best.n_clusters_}クラスタ", fontsize=13)

# HDBSCAN
n_hdbscan_clusters = len(set(hdbscan_labels)) - (1 if -1 in hdbscan_labels else 0)
for cid in range(n_hdbscan_clusters):
    mask = hdbscan_labels == cid
    axes[2].scatter(X_mixed[mask, 0], X_mixed[mask, 1], s=15,
                    edgecolors="k", linewidths=0.2)
noise = hdbscan_labels == -1
if noise.sum() > 0:
    axes[2].scatter(X_mixed[noise, 0], X_mixed[noise, 1],
                    c="gray", marker="x", s=15)
axes[2].set_title(f"HDBSCAN: {n_hdbscan_clusters}クラスタ", fontsize=13)

for ax in axes:
    ax.set_xlabel("$x_1$", fontsize=11)
    ax.set_ylabel("$x_2$", fontsize=11)

plt.suptitle("DBSCAN vs HDBSCAN: 密度の異なるクラスタの検出", fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig("dbscan_vs_hdbscan.png", dpi=150, bbox_inches="tight")
plt.show()

HDBSCANの結果を見ると、εを固定するDBSCANでは困難だった密度の異なるクラスタの検出が、HDBSCANでは適切に行われています。高密度クラスタと低密度クラスタの両方が正しく検出され、ノイズ点も適切に分離されています。HDBSCANはεのチューニングが不要で、min_cluster_size(クラスタの最小サイズ)のみを指定すればよいという実用上の利点もあります。

OPTICSによる可視化

DBSCANのもう一つの拡張としてOPTICS(Ordering Points To Identify the Clustering Structure)があります。OPTICSは到達可能距離のプロットを生成することで、データの密度構造を視覚的に把握できます。

from sklearn.cluster import OPTICS

# OPTICSの実行
optics = OPTICS(min_samples=5, xi=0.05, min_cluster_size=0.05)
optics.fit(X_mixed)

# 到達可能距離プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 到達可能距離プロット
reachability = optics.reachability_[optics.ordering_]
axes[0].bar(range(len(reachability)), reachability, width=1.0,
            color="steelblue", alpha=0.7)
axes[0].set_xlabel("データ点(OPTICSの順序)", fontsize=12)
axes[0].set_ylabel("到達可能距離", fontsize=12)
axes[0].set_title("OPTICSの到達可能距離プロット", fontsize=14)
axes[0].set_ylim(0, np.percentile(reachability[np.isfinite(reachability)], 99) * 1.1)
axes[0].grid(True, alpha=0.3)

# OPTICSのクラスタリング結果
optics_labels = optics.labels_
n_optics_clusters = len(set(optics_labels)) - (1 if -1 in optics_labels else 0)
for cid in range(n_optics_clusters):
    mask = optics_labels == cid
    axes[1].scatter(X_mixed[mask, 0], X_mixed[mask, 1], s=15,
                    edgecolors="k", linewidths=0.2, label=f"クラスタ {cid}")
noise = optics_labels == -1
if noise.sum() > 0:
    axes[1].scatter(X_mixed[noise, 0], X_mixed[noise, 1],
                    c="gray", marker="x", s=15, label="ノイズ")
axes[1].set_title(f"OPTICS: {n_optics_clusters}クラスタ", fontsize=13)
axes[1].set_xlabel("$x_1$", fontsize=11)
axes[1].set_ylabel("$x_2$", fontsize=11)
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.savefig("optics_reachability.png", dpi=150, bbox_inches="tight")
plt.show()

OPTICSの到達可能距離プロットでは、クラスタは「谷」として現れ、クラスタ間の境界は「山」として現れます。谷の深さは密度の高さに対応し、谷の幅はクラスタのサイズに対応します。このプロットにより、データの階層的な密度構造を1つの図で把握でき、異なるεの値でDBSCANを実行した場合の結果を推測することもできます。

クラスタリング評価指標

クラスタリング結果の評価には、真のラベルが利用可能な場合と利用不可能な場合で異なる指標を使います。

from sklearn.metrics import (
    adjusted_rand_score, normalized_mutual_info_score,
    silhouette_score
)

# 三日月形データでの評価
X_eval, y_eval = make_moons(n_samples=500, noise=0.08, random_state=42)

# DBSCAN
db_eval = DBSCANScratch(eps=0.2, min_pts=5)
db_labels = db_eval.fit_predict(X_eval)

# k-means
km_eval = KMeans(n_clusters=2, random_state=42, n_init=10)
km_labels = km_eval.fit_predict(X_eval)

print("=== クラスタリング評価指標 ===\n")
print("【外部指標(真のラベルが必要)】")

# ノイズ点を除いた評価(DBSCANはノイズ点を除外)
non_noise = db_labels != -1
print(f"  DBSCAN - ARI: {adjusted_rand_score(y_eval[non_noise], db_labels[non_noise]):.4f}")
print(f"  DBSCAN - NMI: {normalized_mutual_info_score(y_eval[non_noise], db_labels[non_noise]):.4f}")
print(f"  k-means - ARI: {adjusted_rand_score(y_eval, km_labels):.4f}")
print(f"  k-means - NMI: {normalized_mutual_info_score(y_eval, km_labels):.4f}")

print("\n【内部指標(真のラベル不要)】")
if len(set(db_labels[non_noise])) > 1:
    print(f"  DBSCAN - シルエットスコア: {silhouette_score(X_eval[non_noise], db_labels[non_noise]):.4f}")
print(f"  k-means - シルエットスコア: {silhouette_score(X_eval, km_labels):.4f}")

print(f"\n【補足情報】")
print(f"  DBSCAN - 検出クラスタ数: {db_eval.n_clusters_}")
print(f"  DBSCAN - ノイズ点数: {db_eval.n_noise_} ({db_eval.n_noise_/len(X_eval)*100:.1f}%)")

評価指標の結果を見ると、ARI(Adjusted Rand Index)とNMI(Normalized Mutual Information)の両方でDBSCANがk-meansを大幅に上回っていることが分かります。三日月形データでは、k-meansは本質的に正しいクラスタリングができないためです。一方、シルエットスコアはk-meansの方が高く出ることがあります。これはシルエットスコアが球状のクラスタを前提としているためで、非凸形状のクラスタではシルエットスコアが適切な指標にならないことを示しています。

実応用: 地理データのクラスタリング

最後に、実用的な例として地理データ(緯度・経度)のクラスタリングを行います。

# 模擬的な地理データ(東京周辺の店舗位置データ)
np.random.seed(42)

# 渋谷エリア(高密度)
shibuya = np.random.randn(100, 2) * 0.005 + np.array([139.7005, 35.6580])
# 新宿エリア(高密度)
shinjuku = np.random.randn(80, 2) * 0.006 + np.array([139.6990, 35.6895])
# 池袋エリア(中密度)
ikebukuro = np.random.randn(50, 2) * 0.008 + np.array([139.7107, 35.7295])
# 郊外の孤立した店舗(ノイズ)
outliers = np.random.uniform(
    [139.60, 35.60], [139.80, 35.80], (20, 2)
)

locations = np.vstack([shibuya, shinjuku, ikebukuro, outliers])

# DBSCAN
from sklearn.cluster import DBSCAN as SklearnDBSCAN
geo_dbscan = SklearnDBSCAN(eps=0.012, min_samples=5)
geo_labels = geo_dbscan.fit_predict(locations)

n_geo_clusters = len(set(geo_labels)) - (1 if -1 in geo_labels else 0)

fig, ax = plt.subplots(figsize=(10, 8))
colors_map = plt.cm.tab10(np.linspace(0, 1, 10))
area_names = {0: "エリアA", 1: "エリアB", 2: "エリアC"}

for cid in range(n_geo_clusters):
    mask = geo_labels == cid
    name = area_names.get(cid, f"エリア{cid}")
    ax.scatter(locations[mask, 0], locations[mask, 1],
               color=colors_map[cid], s=30, label=name,
               edgecolors="k", linewidths=0.3)

noise = geo_labels == -1
if noise.sum() > 0:
    ax.scatter(locations[noise, 0], locations[noise, 1],
               color="gray", marker="x", s=30, label="孤立店舗")

ax.set_xlabel("経度", fontsize=12)
ax.set_ylabel("緯度", fontsize=12)
ax.set_title(f"地理データのDBSCANクラスタリング ({n_geo_clusters}エリア検出)", fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("geo_dbscan.png", dpi=150, bbox_inches="tight")
plt.show()

print(f"検出エリア数: {n_geo_clusters}")
print(f"孤立店舗数: {noise.sum()}")

地理データの例では、DBSCANが繁華街(高密度エリア)を自動的にクラスタとして検出し、郊外の孤立した店舗をノイズとして分離しています。k-meansではクラスタ数を事前に指定する必要があり、孤立店舗も無理やりいずれかのクラスタに割り当てられてしまいます。DBSCANは地理データ分析において非常に実用的な手法です。

まとめ

本記事では、DBSCAN(密度ベースクラスタリング)の理論と実装を解説しました。

重要なポイントを振り返ります。

  • DBSCANの基本概念: ε-近傍とMinPtsにより、コア点・ボーダー点・ノイズ点を分類し、密度接続された点の極大集合をクラスタとする
  • k-meansとの違い: クラスタ数の自動決定、任意形状のクラスタ検出、ノイズの自動分離
  • パラメータ選択: k-距離グラフのエルボーでεを推定、MinPtsは $2d$ を目安
  • DBSCANの限界: グローバルなεでは密度の異なるクラスタを同時に検出できない
  • HDBSCANへの拡張: 相互到達距離と階層的クラスタリングにより、密度の異なるクラスタに対応
  • OPTICSによる可視化: 到達可能距離プロットでデータの密度構造を把握

DBSCANは表形式データの探索的分析、地理データ分析、異常検知など、幅広い場面で活用できる強力なアルゴリズムです。データの特性に応じて、k-means、DBSCAN、HDBSCANを使い分けることが重要です。