データをクラスタリングしたいとき、「クラスタの数をいくつに設定すべきか」は厄介な問題です。k-meansやガウス混合モデルでは、分析者がクラスタ数 $K$ をあらかじめ指定する必要があります。しかし実際のデータでは、クラスタ数が事前に分かっていることは稀です。しかも、データが増えるにつれて新しいクラスタが現れる可能性もあります。
もし「クラスタ数をデータに語らせる」ことができたら——この願望を叶えるのがベイズノンパラメトリクス（Bayesian nonparametrics）であり、その中核にあるのがディリクレ過程（Dirichlet process, DP）です。
ベイズノンパラメトリクスを理解すると、以下のような応用が開けます。
- クラスタリング: クラスタ数を自動推定する無限混合モデル（DP混合モデル）
- 密度推定: パラメトリックな仮定を置かない柔軟な確率密度の推定
- 自然言語処理: トピックモデルの拡張（HDP-LDA）、言語モデルのスムージング
- コンピュータビジョン: 画像セグメンテーション、物体認識における無限モデル
本記事の内容
- ベイズノンパラメトリクスの動機 — なぜ無限次元が必要か
- ディリクレ分布からディリクレ過程へ — 有限から無限への自然な拡張
- ディリクレ過程の3つの構成法（棒折り・中華料理店・ポリアの壺）
- ディリクレ過程混合モデル（DPMM）の理論
- Pythonによる棒折り過程とDP混合モデルの実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ベイズノンパラメトリクスとは
パラメトリックモデルの限界
統計モデリングにおいて、私たちは通常「パラメトリックモデル」を使います。正規分布なら平均 $\mu$ と分散 $\sigma^2$ の2つのパラメータ、混合正規分布なら各成分の平均・分散・混合重みという有限個のパラメータでモデルを特定します。
しかし、パラメトリックモデルには本質的な制約があります。モデルの「複雑さ」がパラメータの数で固定されるため、データの複雑さが事前の想定を超えると適切にフィットできません。ガウス混合モデルで $K = 3$ と指定したのにデータが実際には5つのクラスタを含んでいれば、推定結果は歪みます。
BICやAICなどの情報量基準で $K$ を選ぶこともできますが、これは「$K = 1, 2, 3, \ldots$ のモデルをそれぞれ推定して比較する」という手順であり、本質的にはベイズの枠組みの外でモデルの複雑さを扱っています。
ノンパラメトリックの意味
「ノンパラメトリック」という名前はやや誤解を招きます。パラメータがないわけではなく、パラメータの数が有限に固定されていない（事実上無限のパラメータを持つ）という意味です。
より正確に言えば、ベイズノンパラメトリクスでは関数空間や無限次元の対象に対して事前分布を置くことがポイントです。たとえば、「確率分布そのもの」に対して事前分布を定義します。通常のベイズでは $\theta$ という有限次元パラメータに事前分布 $p(\theta)$ を置きますが、ベイズノンパラメトリクスでは確率分布 $G$ そのものに事前分布 $p(G)$ を置くのです。
このアイデアは最初は奇妙に感じるかもしれません。「分布の分布」とはどういう意味なのでしょうか。確率分布 $G$ は無限次元の対象ですから、有限個のパラメータでは指定できません。ここで登場するのがディリクレ過程です。
ベイズノンパラメトリクスの概念がつかめたところで、次にその中核となるディリクレ過程を定義するために、まず有限次元のディリクレ分布から出発しましょう。
ディリクレ分布の復習
有限次元のディリクレ分布
ディリクレ過程を理解するための自然な出発点は、有限次元のディリクレ分布です。
$K$ 次元のディリクレ分布 $\text{Dir}(\alpha_1, \alpha_2, \ldots, \alpha_K)$ は、$K$ 個の非負の値 $(p_1, p_2, \ldots, p_K)$ で $\sum_{k=1}^{K} p_k = 1$ を満たすベクトル——すなわち確率ベクトル——上の分布です。
確率密度関数は次のように書けます。
$$ \begin{equation} f(p_1, \ldots, p_K) = \frac{\Gamma(\alpha_0)}{\prod_{k=1}^{K}\Gamma(\alpha_k)} \prod_{k=1}^{K} p_k^{\alpha_k – 1} \end{equation} $$
ここで $\alpha_0 = \sum_{k=1}^{K}\alpha_k$ は集中度パラメータの合計です。
ディリクレ分布の直感
ディリクレ分布が生成する確率ベクトルの性質は、パラメータ $\alpha_k$ の値によって大きく変わります。
- $\alpha_k$ が全て1のとき: 単体上の一様分布。どんな確率ベクトルも同様に確からしい
- $\alpha_k$ が全て1より大きいとき: 均等に近い確率ベクトルが出やすい（集中）
- $\alpha_k$ が全て1より小さいとき: 頂点付近の確率ベクトルが出やすい（疎、スパース）
とくに $\alpha_k = \alpha/K$（全て等しく、合計が $\alpha$）として $K \to \infty$ とするとどうなるでしょうか。この極限操作こそがディリクレ過程への道を開きます。
ディリクレ過程の定義
確率分布上の分布
ディリクレ過程 $\text{DP}(\alpha, G_0)$ は2つのパラメータで特徴づけられます。
- 集中度パラメータ $\alpha > 0$: 生成される分布がベース分布 $G_0$ にどれだけ近いかを制御する
- ベース分布 $G_0$: 生成される分布の「中心」となる確率分布
$G \sim \text{DP}(\alpha, G_0)$ のとき、$G$ は確率分布です。つまり、ディリクレ過程は「確率分布を生成する確率過程」です。
形式的な定義
ディリクレ過程の定義は、測度空間上の可測分割に対する条件として与えられます。
確率分布 $G$ がディリクレ過程 $\text{DP}(\alpha, G_0)$ に従うとは、任意の有限分割 $\{A_1, A_2, \ldots, A_K\}$（$A_k$ は標本空間の互いに素な可測集合で和集合が全体）に対して、次が成り立つことです。
$$ \begin{equation} (G(A_1), G(A_2), \ldots, G(A_K)) \sim \text{Dir}(\alpha G_0(A_1), \alpha G_0(A_2), \ldots, \alpha G_0(A_K)) \end{equation} $$
この定義は「確率分布 $G$ の有限次元射影がディリクレ分布に従う」ことを意味しています。どんな分割を取っても対応する確率がディリクレ分布に従うという、無限次元への自然な拡張になっています。
定義の直感的な意味
この定義が何を言っているか、具体的に考えてみましょう。
たとえば、実数直線上のベース分布 $G_0 = N(0, 1)$ としましょう。区間 $A_1 = (-\infty, 0)$ と $A_2 = [0, \infty)$ で分割すると、$G_0(A_1) = G_0(A_2) = 0.5$ なので
$$ (G(A_1), G(A_2)) \sim \text{Dir}(\alpha/2, \alpha/2) $$
$\alpha$ が大きいほど $G(A_1) \approx 0.5$ に集中し（ベース分布に近い）、$\alpha$ が小さいほど $(G(A_1), G(A_2))$ は $(0, 1)$ や $(1, 0)$ に近い極端な値をとりやすくなります。
同じ議論がどんな細かい分割に対しても成り立つため、$\alpha$ が小さいとき $G$ は少数の点に質量が集中した「スパイク状」の分布になり、$\alpha$ が大きいとき $G$ はベース分布 $G_0$ に近い滑らかな分布になります。
ここでの重要なポイントは、ディリクレ過程から生成される分布 $G$ がほぼ確実に離散分布になる（確率1で離散）ということです。これは一見すると制限に思えますが、実はクラスタリングへの応用で中心的な役割を果たします。
では、この離散性はどのように現れるのでしょうか。それを明示的に見せてくれるのが、次に説明する棒折り構成法です。
棒折り構成法（Stick-Breaking Construction）
Sethuraman の構成
1994年にSethuramanが示した棒折り構成法は、ディリクレ過程からのサンプル $G$ を具体的に構成する方法です。
$G \sim \text{DP}(\alpha, G_0)$ は次のように表現できます。
$$ \begin{equation} G = \sum_{k=1}^{\infty} \pi_k \delta_{\theta_k} \end{equation} $$
ここで $\delta_{\theta_k}$ は点 $\theta_k$ に質量1を置くディラックのデルタ測度です。パラメータ $\theta_k$ と重み $\pi_k$ は次のように生成されます。
アトムの生成:
$$ \theta_k \stackrel{\text{iid}}{\sim} G_0, \quad k = 1, 2, 3, \ldots $$
重みの生成（棒折り）:
$$ V_k \stackrel{\text{iid}}{\sim} \text{Beta}(1, \alpha), \quad k = 1, 2, 3, \ldots $$
$$ \begin{equation} \pi_k = V_k \prod_{j=1}^{k-1}(1 – V_j) \end{equation} $$
棒折りの直感
「棒折り」という名前の由来を理解するために、長さ1の棒を折る操作を想像しましょう。
- 長さ1の棒がある。まず $V_1 \sim \text{Beta}(1, \alpha)$ の割合で折り、左側の部分を $\pi_1 = V_1$ とする
- 残りの棒（長さ $1 – V_1$）から $V_2$ の割合で折り、$\pi_2 = V_2(1 – V_1)$ とする
- 残りの棒（長さ $(1-V_1)(1-V_2)$）から $V_3$ の割合で折り、$\pi_3 = V_3(1-V_1)(1-V_2)$ とする
- これを無限に繰り返す
$V_k \sim \text{Beta}(1, \alpha)$ の期待値は $1/(1+\alpha)$ なので、$\alpha$ が小さいと各ステップで大きな割合を折り取る（少数の成分に質量が集中）、$\alpha$ が大きいと小さな割合しか折り取らない（多数の成分に質量が分散）ことになります。
重要な性質として、$\sum_{k=1}^{\infty} \pi_k = 1$（確率1で）が保証されています。これは、棒の全体がちょうど使い切られることを意味します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
def stick_breaking(alpha, K_max=50):
"""棒折り過程で重みを生成する"""
V = np.random.beta(1, alpha, size=K_max)
pi = np.zeros(K_max)
remaining = 1.0
for k in range(K_max):
pi[k] = remaining * V[k]
remaining *= (1 - V[k])
return pi
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
alphas = [0.1, 1.0, 5.0, 50.0]
for ax, alpha in zip(axes.ravel(), alphas):
pi = stick_breaking(alpha, K_max=30)
ax.bar(range(len(pi)), pi, color="steelblue", alpha=0.8)
ax.set_title(rf"$\alpha = {alpha}$", fontsize=13)
ax.set_xlabel("Component k", fontsize=11)
ax.set_ylabel(r"$\pi_k$", fontsize=11)
ax.set_xlim(-0.5, 29.5)
top5_sum = np.sum(np.sort(pi)[::-1][:5])
ax.text(0.97, 0.95, f"Top-5 sum: {top5_sum:.3f}",
transform=ax.transAxes, ha="right", va="top", fontsize=10,
bbox=dict(boxstyle="round,pad=0.3", facecolor="lightyellow"))
plt.suptitle("Stick-Breaking Weights for Different Concentration Parameters",
fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("stick_breaking.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、集中度パラメータ $\alpha$ の効果が明確に読み取れます。
-
$\alpha = 0.1$（左上）: 最初の数個の成分にほぼ全ての質量が集中しています。Top-5の重みの合計がほぼ1に達し、実質的には2〜3個の成分しか必要ありません。これは「データが少数のクラスタに分かれる」という事前信念を表現しています
-
$\alpha = 1.0$（右上）: 先頭の成分が最大の重みを持ちますが、その後も緩やかに減衰し、10個程度の成分が無視できない重みを持ちます
-
$\alpha = 5.0$（左下）: 重みがより均等に分散し始めます。多くの成分が小さいながらも有意な重みを持ち、20個以上の成分が必要です
-
$\alpha = 50.0$（右下）: 重みはほぼ均等に分散し、各成分の重みが非常に小さくなります。これは多数の小さなクラスタを想定する事前信念に対応します
このように、$\alpha$ は「事前に想定するクラスタ数の多さ」を制御するパラメータとして解釈できます。
棒折り構成は、ディリクレ過程の離散性を明示的に示してくれます。分布 $G$ は可算無限個の点に質量を持つ離散分布ですが、質量のほとんどは最初の数十個の点に集中するため、実用的には有限の切断（truncation）で近似できます。
では、この離散性がデータ生成の観点からどのような意味を持つのでしょうか。それを明らかにするのが、次の中華料理店過程です。
中華料理店過程（Chinese Restaurant Process）
逐次的なクラスタリング
中華料理店過程（CRP）は、ディリクレ過程の予測分布（predictive distribution）を表現する、極めて直感的な比喩です。
無限個のテーブルがある中華料理店を想像してください。客が一人ずつ入店し、テーブルに座ります。
- 最初の客: 必ず最初のテーブル（テーブル1）に座る
- $n+1$ 番目の客: 以下の確率でテーブルを選ぶ
$$ \begin{equation} P(\text{テーブル } k \text{ に座る}) = \begin{cases} \frac{n_k}{n + \alpha} & \text{既存のテーブル } k \text{ に } n_k \text{ 人がいる場合} \\ \frac{\alpha}{n + \alpha} & \text{新しいテーブルを開く場合} \end{cases} \end{equation} $$
ここで $n_k$ はテーブル $k$ に既に座っている客数、$n$ は入店済みの客の総数、$\alpha$ は集中度パラメータです。
CRPの性質
この確率割り当てには、いくつかの重要な性質があります。
「金持ちはますます金持ちに」効果: 既に多くの客が座っているテーブルほど、新しい客も座る確率が高くなります。$n_k/( n + \alpha)$ は $n_k$ に比例するため、人気のテーブルはますます人気になります。これは「優先的選択」（preferential attachment）や「ポリアの壺」と同じメカニズムです。
新しいテーブルの確率: $\alpha/(n + \alpha)$ は $n$ が増えるにつれて減少しますが、決して0にはなりません。つまり、データが増えても新しいクラスタが常に現れる可能性があります。
テーブル数の期待値: $n$ 人の客が座った後のテーブル数の期待値は
$$ \begin{equation} E[K_n] = \alpha \sum_{i=1}^{n} \frac{1}{\alpha + i – 1} \approx \alpha \ln\left(\frac{n + \alpha}{\alpha}\right) \approx \alpha \ln n \quad (n \gg \alpha) \end{equation} $$
テーブル数は $n$ に対して対数的に増加します。1000人の客がいても、$\alpha = 1$ なら平均的なテーブル数はわずか $\ln(1000) \approx 7$ 個程度です。
CRPとディリクレ過程の関係
中華料理店過程は、ディリクレ過程の「積分表現」に対応します。
$G \sim \text{DP}(\alpha, G_0)$ から独立に $\theta_1, \theta_2, \ldots$ を生成する場合、$\theta_i$ の条件付き分布（既に $\theta_1, \ldots, \theta_n$ が観測されたとき）は次のようになります。
$$ \begin{equation} \theta_{n+1} | \theta_1, \ldots, \theta_n \sim \frac{\alpha}{\alpha + n} G_0 + \frac{1}{\alpha + n}\sum_{i=1}^{n} \delta_{\theta_i} \end{equation} $$
$G$ を積分消去（marginalize out）して得られるこの予測分布はブラックウェル＝マッケイの式（Blackwell-MacQueen Polya urn scheme）と呼ばれます。
右辺の第1項はベース分布 $G_0$ から新しい値を生成する（= 新しいテーブルを開く）確率 $\alpha/(\alpha + n)$ に対応し、第2項は既存の値 $\theta_i$ の一つと同じ値をとる（= 既存のテーブルに座る）確率 $1/(\alpha + n)$ に対応します。同じ値 $\theta^*$ を持つ $\theta_i$ が $n_k$ 個あれば、その値が選ばれる確率は $n_k/(\alpha + n)$ となり、まさにCRPの規則が再現されます。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
def chinese_restaurant_process(n, alpha):
"""中華料理店過程をシミュレーションする"""
tables = [1] # 最初の客はテーブル1に座る
assignments = [0] # 各客のテーブル番号
for i in range(1, n):
probs = np.array(tables + [alpha], dtype=float)
probs /= probs.sum()
chosen = np.random.choice(len(probs), p=probs)
if chosen == len(tables):
# 新しいテーブル
tables.append(1)
assignments.append(len(tables) - 1)
else:
tables[chosen] += 1
assignments.append(chosen)
return tables, assignments
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
n = 500
alphas = [0.5, 1.0, 5.0, 20.0]
for ax, alpha in zip(axes.ravel(), alphas):
tables, assignments = chinese_restaurant_process(n, alpha)
K = len(tables)
sorted_tables = sorted(tables, reverse=True)
ax.bar(range(K), sorted_tables, color="coral", alpha=0.8, edgecolor="darkred",
linewidth=0.5)
ax.set_title(rf"$\alpha = {alpha}$, Tables = {K}", fontsize=13)
ax.set_xlabel("Table (sorted by size)", fontsize=11)
ax.set_ylabel("Number of customers", fontsize=11)
ax.grid(True, alpha=0.3, axis="y")
plt.suptitle(f"Chinese Restaurant Process (n = {n} customers)",
fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("crp.png", dpi=150, bbox_inches="tight")
plt.show()
このシミュレーション結果から、CRPの特徴が読み取れます。
-
$\alpha = 0.5$: テーブル数が非常に少なく（3〜5程度）、最大のテーブルに客の大半が集中しています。「金持ちはますます金持ちに」の効果が強く表れています
-
$\alpha = 1.0$: テーブル数はやや増えますが（5〜8程度）、依然として少数の大きなテーブルが支配的です。対数的な増加 $\alpha \ln n \approx 6.2$ とおおよそ一致します
-
$\alpha = 5.0$: テーブル数がかなり増え（20〜35程度）、サイズの分布もより均等になります
-
$\alpha = 20.0$: 多数のテーブルに比較的均等に客が分散し、各テーブルのサイズは小さくなります
CRPが「テーブルへの客の割り当て」としてクラスタリングを定義することを見てきました。しかし、これまでは離散的な構造（クラスタ番号）の話でした。実際のデータ分析では、各クラスタにパラメータを割り当て、連続値のデータをモデル化する必要があります。それを実現するのが、次に紹介するディリクレ過程混合モデルです。
ポリアの壺図式（Polya Urn Scheme）
色付きの玉による解釈
ディリクレ過程のもう一つの直感的な構成法として、ポリアの壺図式があります。
壺に色付きの玉が入っています。以下の操作を繰り返します。
- 壺から玉を1つ取り出す
- その玉と同じ色の玉を新たに1つ加えて、元の玉と合わせて壺に戻す（つまり同じ色が1つ増える）
- ただし、確率 $\alpha/(\alpha + n)$ で「新しい色」の玉を追加する
ここで「色」が $\theta$ の値に対応し、新しい色はベース分布 $G_0$ から生成されます。
数学的には、CRPと全く同じ予測分布を持ちますが、ポリアの壺図式では「値の生成」の側面が強調されます。同じ色の玉が増えるほど次にその色が引かれやすくなるため、出現頻度の高い値がさらに出やすくなるという「正のフィードバック」が自然に組み込まれています。
交換可能性
ポリアの壺図式（そしてCRP）が持つ重要な性質は交換可能性（exchangeability）です。客がレストランに入る順番や、壺から玉を取り出す順番を入れ替えても、最終的なクラスタリングの分布は変わりません。
この性質はデ・フィネッティの表現定理と深く関係しています。交換可能な列 $\theta_1, \theta_2, \ldots$ に対して、ある確率分布 $G$ が存在して
$$ \theta_i | G \stackrel{\text{iid}}{\sim} G $$
と書けます。ディリクレ過程の場合、この混合分布 $G$ 自体が $\text{DP}(\alpha, G_0)$ に従うのです。
交換可能性は実用上も重要です。データの到着順序に結果が依存しないため、推論アルゴリズムの設計が簡潔になります。
ここまでで、ディリクレ過程の3つの顔——棒折り、中華料理店、ポリアの壺——を見てきました。これらは全て同じ確率過程の異なる表現です。次に、これらをデータ分析の道具として使うために、ディリクレ過程混合モデルを定式化しましょう。
ディリクレ過程の性質
期待値と分散
ディリクレ過程 $G \sim \text{DP}(\alpha, G_0)$ の基本的な性質を整理します。
任意の可測集合 $A$ に対して、$G(A)$ の期待値と分散は次のように計算できます。
$$ \begin{equation} E[G(A)] = G_0(A) \end{equation} $$
$$ \begin{equation} \text{Var}[G(A)] = \frac{G_0(A)(1 – G_0(A))}{\alpha + 1} \end{equation} $$
期待値がベース分布 $G_0(A)$ に等しいことは、ディリクレ過程がベース分布を「中心」として生成することを意味します。分散は $\alpha$ が大きいほど小さくなり、$G$ が $G_0$ に集中することを定量的に示しています。
これらの式の導出はディリクレ分布の性質から直接得られます。分割 $\{A, A^c\}$ に対して $(G(A), G(A^c)) \sim \text{Dir}(\alpha G_0(A), \alpha G_0(A^c))$ ですから、ベータ分布の期待値・分散の公式を使えば
$$ E[G(A)] = \frac{\alpha G_0(A)}{\alpha G_0(A) + \alpha G_0(A^c)} = \frac{\alpha G_0(A)}{\alpha} = G_0(A) $$
$$ \text{Var}[G(A)] = \frac{\alpha G_0(A) \cdot \alpha G_0(A^c)}{(\alpha)^2(\alpha + 1)} = \frac{G_0(A)(1 – G_0(A))}{\alpha + 1} $$
と得られます。
事後分布の共役性
ディリクレ過程の最も美しい性質の一つは、自分自身と共役であることです。
$G \sim \text{DP}(\alpha, G_0)$ から $\theta_1, \ldots, \theta_n \stackrel{\text{iid}}{\sim} G$ を観測したとき、$G$ の事後分布は次のようになります。
$$ \begin{equation} G | \theta_1, \ldots, \theta_n \sim \text{DP}\left(\alpha + n, \; \frac{\alpha}{\alpha + n}G_0 + \frac{n}{\alpha + n}\hat{G}_n\right) \end{equation} $$
ここで $\hat{G}_n = \frac{1}{n}\sum_{i=1}^{n}\delta_{\theta_i}$ は経験分布です。
この式の意味を読み解きましょう。事後分布は再びディリクレ過程になり、集中度パラメータは $\alpha + n$ に増加し、ベース分布は事前のベース分布 $G_0$ と経験分布 $\hat{G}_n$ の重み付き平均になります。データが増えるにつれて（$n \to \infty$）、事後のベース分布は経験分布に近づき、$G$ はデータの経験分布に収束します。
この共役性は、ディリクレ過程が理論的に扱いやすい理由の一つであり、推論アルゴリズムの設計にも活用されます。
ディリクレ過程混合モデル（DPMM）
モデルの定式化
ディリクレ過程混合モデル（Dirichlet Process Mixture Model, DPMM）は、ディリクレ過程を混合モデルの事前分布として使うモデルです。有限混合モデルのクラスタ数 $K$ を固定する代わりに、ディリクレ過程を用いて事実上無限のクラスタを許容します。
モデルは次のように定式化されます。
$$ \begin{align} G &\sim \text{DP}(\alpha, G_0) \\ \theta_i &\sim G, \quad i = 1, \ldots, n \\ x_i &\sim F(\cdot | \theta_i), \quad i = 1, \ldots, n \end{align} $$
ここで $F(\cdot | \theta_i)$ はパラメータ $\theta_i$ を持つ観測モデル（たとえば正規分布 $N(\mu_i, \sigma_i^2)$）です。
$G$ が離散分布であるため、$\theta_i$ のうちのいくつかは同じ値を取ります。$\theta_i = \theta_j$ であれば、データ点 $x_i$ と $x_j$ は同じクラスタに属します。
棒折り表現
棒折り構成法を使うと、DPMMは次のように書けます。
$$ \begin{align} V_k &\sim \text{Beta}(1, \alpha), \quad k = 1, 2, \ldots \\ \pi_k &= V_k \prod_{j < k}(1 - V_j) \\ \theta_k^* &\sim G_0, \quad k = 1, 2, \ldots \\ z_i &\sim \text{Categorical}(\pi_1, \pi_2, \ldots), \quad i = 1, \ldots, n \\ x_i &\sim F(\cdot | \theta_{z_i}^*), \quad i = 1, \ldots, n \end{align} $$
ここで $z_i \in \{1, 2, \ldots\}$ はデータ点 $x_i$ のクラスタ割り当て変数です。実際には有限の $K$ で切断して近似します。
ガウス混合への適用
最も一般的なDPMMの例として、ガウスDPMMを考えます。
$$ \begin{align} G &\sim \text{DP}(\alpha, \text{NIW}(\mu_0, \kappa_0, \nu_0, \Psi_0)) \\ (\mu_k, \Sigma_k) &\sim G \\ x_i | z_i = k &\sim N(\mu_k, \Sigma_k) \end{align} $$
ここでNIW（Normal-Inverse-Wishart）は正規分布の平均と共分散行列に対する共役事前分布です。
有限ガウス混合モデルとの対比は明確です。
| 有限ガウス混合 | DP混合 | |
|---|---|---|
| クラスタ数 | $K$（固定） | $\infty$（データに適応） |
| 重み | $(\pi_1, \ldots, \pi_K) \sim \text{Dir}(\alpha/K, \ldots, \alpha/K)$ | 棒折り過程で生成 |
| 推論 | EM or ギブスサンプリング | ギブスサンプリング or 変分推論 |
| モデル選択 | BIC/AICで $K$ を選ぶ | $\alpha$ を推定（クラスタ数は自動決定） |
推論アルゴリズム
DPMMの推論には主に2つのアプローチがあります。
ギブスサンプリング（CRPベース）: $G$ を積分消去し、CRPの予測分布を利用してクラスタ割り当て $z_i$ をギブスサンプリングで更新します。各データ点 $x_i$ に対して
$$ P(z_i = k | z_{-i}, x, \alpha) \propto \begin{cases} n_{-i,k} \cdot F(x_i | \theta_k^*) & \text{既存クラスタ } k \\ \alpha \cdot \int F(x_i | \theta) \, dG_0(\theta) & \text{新規クラスタ} \end{cases} $$
ここで $z_{-i}$ は $z_i$ を除いた割り当て、$n_{-i,k}$ はクラスタ $k$ の（$x_i$ を除いた）データ数です。
変分推論（棒折りベース）: 棒折り表現を有限次元で切断し、変分推論を適用します。計算量が少なく大規模データに適しています。
Pythonによるディリクレ過程混合モデルの実装
シンプルな1次元DP混合モデル
ここでは、CRPベースのギブスサンプリングによるDP混合モデルを1次元ガウス混合データに対して実装します。推論アルゴリズムの各ステップが理論とどう対応するかを確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, t as t_dist
np.random.seed(42)
# --- 合成データの生成 ---
# 4つのクラスタから生成（モデルにはクラスタ数を教えない）
true_means = [-5.0, 0.0, 4.0, 8.0]
true_stds = [0.8, 1.2, 0.6, 1.0]
true_weights = [0.2, 0.35, 0.25, 0.2]
n = 300
data = []
true_labels = []
for i, (mu, std, w) in enumerate(zip(true_means, true_stds, true_weights)):
n_k = int(n * w)
data.extend(np.random.normal(mu, std, n_k))
true_labels.extend([i] * n_k)
data = np.array(data)
true_labels = np.array(true_labels)
n = len(data)
# --- CRPベースのギブスサンプリング ---
# ハイパーパラメータ
alpha = 1.0 # 集中度パラメータ
mu_0 = 0.0 # 事前平均
kappa_0 = 0.1 # 事前精度の重み
nu_0 = 3.0 # 事前の自由度
sigma2_0 = 4.0 # 事前の分散スケール
def marginal_likelihood(x_cluster, mu_0, kappa_0, nu_0, sigma2_0):
"""クラスタのデータの周辺尤度（NIW事前で積分消去後）"""
n_k = len(x_cluster)
if n_k == 0:
return 0.0
x_bar = np.mean(x_cluster)
kappa_n = kappa_0 + n_k
nu_n = nu_0 + n_k
mu_n = (kappa_0 * mu_0 + n_k * x_bar) / kappa_n
sigma2_n = (nu_0 * sigma2_0
+ np.sum((x_cluster - x_bar)**2)
+ kappa_0 * n_k * (x_bar - mu_0)**2 / kappa_n) / nu_n
# t分布の尤度
return np.sum(t_dist.logpdf(x_cluster, df=nu_n, loc=mu_n,
scale=np.sqrt(sigma2_n * (1 + 1/kappa_n))))
def prior_predictive(x, mu_0, kappa_0, nu_0, sigma2_0):
"""新規クラスタの事前予測分布"""
return t_dist.logpdf(x, df=nu_0, loc=mu_0,
scale=np.sqrt(sigma2_0 * (1 + 1/kappa_0)))
# 初期化: 全データを1つのクラスタに
z = np.zeros(n, dtype=int)
n_iter = 200
# クラスタ数の履歴
K_history = []
for iteration in range(n_iter):
for i in range(n):
# データ点iを現在のクラスタから除去
old_k = z[i]
z_temp = np.delete(z, i)
x_temp = np.delete(data, i)
# アクティブなクラスタを取得
unique_clusters = np.unique(z_temp)
log_probs = []
for k in unique_clusters:
mask = z_temp == k
n_k = np.sum(mask)
# CRP確率 × 尤度
log_prior = np.log(n_k)
x_cluster = x_temp[mask]
log_lik = t_dist.logpdf(data[i], df=nu_0 + n_k,
loc=(kappa_0*mu_0 + n_k*np.mean(x_cluster)) / (kappa_0 + n_k),
scale=np.sqrt(
((nu_0*sigma2_0 + np.sum((x_cluster - np.mean(x_cluster))**2)
+ kappa_0*n_k*(np.mean(x_cluster) - mu_0)**2/(kappa_0+n_k))
/ (nu_0 + n_k))
* (1 + 1/(kappa_0 + n_k))
))
log_probs.append(log_prior + log_lik)
# 新規クラスタ
log_prior_new = np.log(alpha)
log_lik_new = prior_predictive(data[i], mu_0, kappa_0, nu_0, sigma2_0)
log_probs.append(log_prior_new + log_lik_new)
# 正規化してサンプリング
log_probs = np.array(log_probs)
log_probs -= np.max(log_probs)
probs = np.exp(log_probs)
probs /= probs.sum()
choices = list(unique_clusters) + [max(unique_clusters) + 1 if len(unique_clusters) > 0 else 0]
z[i] = np.random.choice(choices, p=probs)
# クラスタ番号をリラベル
unique_k = np.unique(z)
mapping = {old: new for new, old in enumerate(unique_k)}
z = np.array([mapping[k] for k in z])
K_history.append(len(unique_k))
print(f"Final number of clusters: {K_history[-1]}")
print(f"Cluster sizes: {[np.sum(z == k) for k in range(K_history[-1])]}")
# --- 結果の可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
# (a) クラスタリング結果
ax = axes[0]
colors = plt.cm.Set2(np.linspace(0, 1, max(K_history[-1], 4)))
for k in range(K_history[-1]):
mask = z == k
ax.scatter(data[mask], np.zeros(np.sum(mask)) + np.random.normal(0, 0.05, np.sum(mask)),
color=colors[k], alpha=0.6, s=20, label=f"Cluster {k+1} (n={np.sum(mask)})")
ax.set_xlabel("x", fontsize=12)
ax.set_title("DPMM Clustering Result", fontsize=13)
ax.legend(fontsize=8, loc="upper right")
ax.grid(True, alpha=0.3)
# (b) 真のクラスタとの比較
ax = axes[1]
for k in range(4):
mask = true_labels == k
ax.scatter(data[mask], np.zeros(np.sum(mask)) + np.random.normal(0, 0.05, np.sum(mask)),
color=colors[k], alpha=0.6, s=20, label=f"True {k+1} (n={np.sum(mask)})")
ax.set_xlabel("x", fontsize=12)
ax.set_title("True Clusters", fontsize=13)
ax.legend(fontsize=8, loc="upper right")
ax.grid(True, alpha=0.3)
# (c) クラスタ数の推移
ax = axes[2]
ax.plot(K_history, color="steelblue", linewidth=1.5)
ax.axhline(4, color="red", linestyle="--", linewidth=1.5, label="True K = 4")
ax.set_xlabel("Iteration", fontsize=12)
ax.set_ylabel("Number of clusters", fontsize=12)
ax.set_title("Convergence of Cluster Count", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("dpmm_result.png", dpi=150, bbox_inches="tight")
plt.show()
上のグラフから、DP混合モデルの振る舞いが確認できます。
-
左図（クラスタリング結果）: ギブスサンプリングにより、クラスタ数を事前に指定することなく、データ中の構造が自動的に検出されています。最終的なクラスタ数は真のクラスタ数4に近い値に収束していることが多いですが、集中度パラメータ $\alpha$ の設定やデータの重なり具合によって変動します
-
中央図（真のクラスタ）: 真のクラスタ構造との比較です。明確に分離した4つのクラスタが存在し、DPMM がこれをおおよそ正しく復元できていることが確認できます
-
右図（クラスタ数の推移）: ギブスサンプリングの反復に伴い、クラスタ数が真の値 $K = 4$（赤い破線）の近辺で安定化しています。初期の不安定な挙動の後、定常分布に収束していく様子が観察できます
ディリクレ過程の発展と関連モデル
階層ディリクレ過程（HDP）
複数のグループがあり、各グループ内でクラスタリングを行いたいが、グループ間でクラスタを共有したい場合に使われるのが階層ディリクレ過程（Hierarchical Dirichlet Process, HDP）です。
$$ \begin{align} G_0 &\sim \text{DP}(\gamma, H) \\ G_j &\sim \text{DP}(\alpha, G_0), \quad j = 1, \ldots, J \\ \theta_{ji} &\sim G_j \end{align} $$
HDPの代表的な応用はトピックモデルです。各文書 $j$ がDP $G_j$ を持ち、文書のトピック分布はDP混合モデルで表現されます。$G_0$ がグローバルなトピック集合を定義するため、異なる文書間でトピックが共有されます。
Pitman-Yor過程
ディリクレ過程のべき乗則的な一般化がPitman-Yor過程（PYP）です。DPに追加のパラメータ $d \in [0, 1)$（割引パラメータ）を導入します。
$$ V_k \sim \text{Beta}(1 – d, \alpha + kd) $$
$d = 0$ のときDPに帰着します。$d > 0$ のとき、テーブルサイズの分布がべき乗則に従うようになり、自然言語のジップ則のような現象をモデル化できます。
Indian Buffet Process（IBP）
ディリクレ過程がクラスタリング（各データ点が1つのクラスタに属する）のためのモデルであるのに対し、Indian Buffet Processは特徴選択（各データ点が複数の特徴を持つ）のためのモデルです。DPの「中華料理店過程」に対応する「インド料理ビュッフェ過程」として位置づけられます。
これらの発展モデルは、ベイズノンパラメトリクスの柔軟性が様々な問題設定に適応できることを示しています。
まとめ
本記事では、ベイズノンパラメトリクスの中核であるディリクレ過程について解説しました。
- ベイズノンパラメトリクスは、パラメータの数を固定せず、データに応じてモデルの複雑さを自動調整する枠組みである
- ディリクレ過程 $\text{DP}(\alpha, G_0)$ は確率分布上の分布であり、集中度パラメータ $\alpha$ とベース分布 $G_0$ で特徴づけられる
- 棒折り構成法は、DPからのサンプルが無限個の点に質量を持つ離散分布であることを明示的に示す。$\alpha$ が小さいほど質量が少数の点に集中する
- 中華料理店過程は、DPの予測分布を逐次的なクラスタリングとして表現し、「金持ちはますます金持ちに」効果と新規クラスタの生成を自然に組み込む
- DP混合モデルは、有限混合モデルのクラスタ数を無限に拡張し、クラスタ数をデータから自動推定するモデルである
- DPの事後分布は再びDPとなる（共役性）ため、理論的にも計算的にも扱いやすい
次のステップとして、以下の記事も参考にしてください。