データをいくつかのまとまり（クラスタ）に分けたいとき、最初に立ちはだかる素朴な、しかし厄介な問いがあります。「クラスタはいくつあるのか？」 です。k-means や通常のガウス混合モデルでは、クラスタ数 $K$ を人間が事前に決めてアルゴリズムに渡さなければなりません。しかし、新しい顧客データや観測データを前にして「これは3グループだ」と断言できる場面はそう多くありません。グループが2つかもしれないし、10個かもしれない。さらに、データが増えるにつれて新しいタイプが現れることもあります。
たとえばニュース記事を話題ごとに分類したいとします。今日のニュースは「政治・経済・スポーツ」の3話題かもしれませんが、来月には「新しい感染症」という話題が突然増えるかもしれません。あるいは天体観測で得たスペクトルを分類するとき、未知の天体種別が紛れ込んでいる可能性があります。こうした「クラスタ数そのものをデータから学びたい」「データが増えたら自然にクラスタが増えてほしい」という要求に応えるのが、本記事で扱うディリクレ過程混合モデル（Dirichlet Process Mixture, DPM） です。
ディリクレ過程混合モデルは、
- 顧客セグメンテーション — 何種類の顧客層がいるか分からないまま購買データから層を発見する
- トピックモデル（自然言語処理） — 文書集合に潜む話題数を事前に決めずに推定する（階層ディリクレ過程, HDP）
といった応用で広く使われています。本記事では、有限のガウス混合モデルの「成分数を無限に飛ばす」極限としてディリクレ過程を導出し、その背後にある美しいメタファー「中華料理店過程（Chinese Restaurant Process, CRP）」を周辺化によって導きます。そして CRP に基づくギブスサンプリングを Python でスクラッチ実装し、真のクラスタ数を一切教えずに2次元データをクラスタリングしてみます。
左の有限GMM ではクラスタ数 K=3 を人間がアルゴリズムに手渡しています。一方、右の DPM ではデータが増えると新しいクラスタ（赤い矢印で示した第4クラスタ）が自動的に生まれるイメージです。この「クラスタ数をデータが決める」という仕組みをどう実現するかが、本記事の核心になります。
本記事の内容
- 有限ディリクレ混合の極限としてのディリクレ過程・スティック折り過程の導出
- 中華料理店過程（CRP）の予測分布を周辺化から導く
- DPM のギブスサンプリングを Python で実装し、クラスタ数を自動推定する
- 集中度 $\alpha$ に対する推定クラスタ数の変化と収束を可視化して読み取る
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が格段に深まります。ディリクレ分布の性質と、ディリクレ・カテゴリカルの共役関係、そしてギブスサンプリングは本記事の土台です。
有限混合モデルからの出発
ディリクレ過程をいきなり「無限次元の確率分布の分布」と定義されても、何のことか分かりません。ここでは慣れ親しんだ有限混合モデルから出発し、成分数を増やしていったらどうなるか、という素朴な思考実験から始めます。これがディリクレ過程の最も直感的な入り口です。
$K$ 個の成分を持つ混合モデルを考えます。各データ点 $x_i$（$i = 1, \dots, N$）は、まずどの成分から生成されるかを表す潜在変数 $z_i \in \{1, \dots, K\}$ がカテゴリカル分布から選ばれ、その成分のパラメータ $\bm{\theta}_{z_i}$ に従って生成されると考えます。生成過程を書き下すと次のようになります。
$$ \begin{align} \bm{\pi} &\sim \mathrm{Dir}(\alpha/K, \dots, \alpha/K) \\ z_i \mid \bm{\pi} &\sim \mathrm{Categorical}(\bm{\pi}) \\ \bm{\theta}_k &\sim H \quad (k = 1, \dots, K) \\ x_i \mid z_i, \{\bm{\theta}_k\} &\sim F(\bm{\theta}_{z_i}) \end{align} $$
ここで $\bm{\pi} = (\pi_1, \dots, \pi_K)$ は各成分の混合比（足して1になる確率ベクトル）、$H$ は成分パラメータの基底分布（事前分布）、$F(\bm{\theta})$ はパラメータ $\bm{\theta}$ の観測分布（たとえばガウス）です。重要なのは混合比 $\bm{\pi}$ の事前分布を、各成分に対して対称な濃度 $\alpha/K$ を持つ対称ディリクレ分布にした点です。総和を $\alpha$ に固定したまま成分数 $K$ を変えられるよう、あえて $\alpha/K$ という形にしてあります。
なぜ $\alpha/K$ という形にするのか、と疑問に思うかもしれません。これは後で $K \to \infty$ の極限を取るための仕掛けです。ディリクレ分布の濃度パラメータの総和が $\alpha$ で一定なら、成分数を増やしても全体の「バラつきの強さ」が一定に保たれます。総和を $\alpha$ に固定しないと、$K$ を増やすにつれて分布が無限に尖ったり平坦になったりして、極限が定義できなくなってしまうのです。
このグラフは、$\alpha=2$ を固定したまま成分数 $K$ を 5・20・200 と増やしたときの混合比 $\bm{\pi}$ の分布です。$K$ が増えるほど各成分の比は小さくなりますが、実質的に使われる成分数（混合比が十分大きい成分）はほぼ一定に保たれています。$\alpha/K$ という設計のおかげで、$K$ をいくら大きくしても情報が散らばり続けるのではなく、有限個のクラスタに集中する性質が維持されます。これが $K\to\infty$ の極限を意味のあるものにする仕掛けです。
ここまでで有限混合モデルの生成過程を確認しました。次に、この混合比 $\bm{\pi}$ を積分消去（周辺化）すると何が起こるかを見ていきます。これがディリクレ過程の隠れた骨格を明らかにします。
混合比を周辺化する — 予測分布の導出
混合モデルで推論する際、混合比 $\bm{\pi}$ そのものに興味があることは稀です。本当に知りたいのは「各データ点がどの成分に属するか」という割り当て $\bm{z} = (z_1, \dots, z_N)$ です。そこで $\bm{\pi}$ を積分で消去してしまい、割り当て $\bm{z}$ の分布だけを残すことを考えます。これは共役性のおかげで解析的に実行できます。
ディリクレ・カテゴリカルの共役関係から、$\bm{\pi}$ を周辺化した割り当て全体の同時分布は次の形になります。$n_k$ を成分 $k$ に割り当てられたデータ点の数とすると、
$$ p(\bm{z} \mid \alpha) = \int p(\bm{z} \mid \bm{\pi})\, p(\bm{\pi} \mid \alpha)\, d\bm{\pi} = \frac{\Gamma(\alpha)}{\Gamma(\alpha + N)} \prod_{k=1}^{K} \frac{\Gamma(n_k + \alpha/K)}{\Gamma(\alpha/K)} $$
となります。ここで $\Gamma(\cdot)$ はガンマ関数です。この式自体はディリクレ・多項分布（ディリクレ・カテゴリカル）の正規化定数から直接得られます。
私たちが本当に欲しいのは、すでに $n_{-i,k}$ 個（自分 $i$ を除いた成分 $k$ の所属数）が決まっているときに、$i$ 番目のデータ点が成分 $k$ に割り当てられる条件付き確率です。ギブスサンプリングではこの形が必要になります。上の同時分布の比を取ると、
$$ p(z_i = k \mid \bm{z}_{-i}, \alpha) = \frac{n_{-i,k} + \alpha/K}{N – 1 + \alpha} $$
が得られます。ここで $\bm{z}_{-i}$ は $i$ 以外の全割り当て、$n_{-i,k}$ は $i$ を除いて成分 $k$ に属する点の数です。
この比の導出を一行ずつ追っておきましょう。同時分布 $p(\bm{z})$ の中で、$z_i = k$ とした場合とそうでない場合の違いは、成分 $k$ のカウント $n_k$ が1だけ増えるかどうかだけに現れます。ガンマ関数の性質 $\Gamma(n+1) = n\,\Gamma(n)$ を使うと、$z_i$ に依存する項だけが $\Gamma(n_{-i,k} + 1 + \alpha/K) / \Gamma(n_{-i,k} + \alpha/K) = n_{-i,k} + \alpha/K$ として残ります。さらに分母の $\Gamma(\alpha + N)/\Gamma(\alpha + N – 1) = \alpha + N – 1$ が規格化に効いてくるため、上の条件付き確率が得られるわけです。
ここまでで「自分を除いた所属数」に比例して所属先が決まるという美しい構造が見えてきました。次に、いよいよこの式で $K \to \infty$ の極限を取り、無限個の成分を持つモデルへと飛躍します。
無限への飛躍 — 中華料理店過程
有限混合モデルの予測分布 $p(z_i = k \mid \bm{z}_{-i}) = (n_{-i,k} + \alpha/K)/(N – 1 + \alpha)$ で、成分数 $K$ を無限大に飛ばすとどうなるでしょうか。直感的には「無限個の引き出しを用意したが、実際にデータが入っているのはごく一部」という状況になります。この極限を取る操作こそが、ディリクレ過程の正体を露わにします。
極限を取る際、成分を2種類に分けて考えます。
(1) すでにデータが入っている成分（$n_{-i,k} > 0$ の成分）に割り当てられる確率は、$K \to \infty$ で $\alpha/K \to 0$ となるため、
$$ p(z_i = k \mid \bm{z}_{-i}) = \frac{n_{-i,k} + \alpha/K}{N – 1 + \alpha} \;\xrightarrow{K \to \infty}\; \frac{n_{-i,k}}{N – 1 + \alpha} $$
に収束します。つまり既存成分には所属人数に比例して割り当てられます。
(2) まだ誰も入っていない成分（$n_{-i,k} = 0$ の成分）は無限個あります。それらをまとめて「新しい成分に割り当てられる確率」を計算します。空の成分1個あたりの確率は $\frac{\alpha/K}{N-1+\alpha}$ で、空の成分は約 $K – K_{+}$ 個（$K_+$ は使用中の成分数、有限）あります。これらを合計すると、
$$ p(z_i = \text{新成分} \mid \bm{z}_{-i}) = \frac{(K – K_+)\cdot \alpha/K}{N – 1 + \alpha} \;\xrightarrow{K \to \infty}\; \frac{\alpha}{N – 1 + \alpha} $$
となります。ここで $K – K_+ \approx K$（$K \to \infty$ で $K_+$ は無視できる）を使うと、$K$ が約分されて $\alpha/(N-1+\alpha)$ という有限の値が残ります。$K$ が消えてくれるのが、$\alpha/K$ という形を選んだご利益です。
この2つを合わせると、中華料理店過程（CRP） の予測分布が得られます。
$$ \begin{equation} p(z_i = k \mid \bm{z}_{-i}, \alpha) = \begin{cases} \dfrac{n_{-i,k}}{N – 1 + \alpha} & (\text{既存の成分 } k) \\[2ex] \dfrac{\alpha}{N – 1 + \alpha} & (\text{新しい成分}) \end{cases} \end{equation} $$
なぜ「中華料理店」なのか
この式には有名なメタファーがあります。無限にテーブルがある中華料理店を想像してください。客が一人ずつ入店します。
- $i$ 番目の客は、すでに人が座っているテーブル $k$ に、そのテーブルの客数 $n_{-i,k}$ に比例する確率で座る（人気のテーブルほど座りやすい — 「金持ちはさらに金持ちに」効果）
- あるいは、確率 $\propto \alpha$ で誰も座っていない新しいテーブルに座る
集中度パラメータ $\alpha$ は「新しいテーブルへの座りやすさ」を表します。$\alpha$ が大きいほど客は新しいテーブルを好み、結果としてテーブル数（=クラスタ数）が増えます。逆に $\alpha$ が小さいと既存テーブルに集まり、少数の大きなクラスタができます。
図は $N=40$ 人、$\alpha=2$ での CRP シミュレーション結果です。各円がテーブルを表し、中の数字がそのテーブルの客数です。最初は1テーブルから始まりますが、確率 $\propto\alpha$ で新テーブルが次々に生まれます。人数の多いテーブルほど新客を引き付けやすく（「金持ちがさらに金持ちに」効果）、結果として少数の大きなテーブルと少数の新テーブルが共存するべき分布が自然に生まれることが分かります。
CRP の重要な性質を2つ挙げておきます。
第一に、交換可能性（exchangeability） です。客の入店順序を入れ替えても、最終的なテーブル分割の確率は変わりません。これにより、ギブスサンプリングで「ある客 $i$ を最後に入店した客とみなして再着席させる」操作が正当化されます。各点を「最後の客」として扱えるからこそ、上の予測分布を全ての点に対して使えるのです。
第二に、テーブル数（クラスタ数）の期待値は $N$ 人の客に対しておよそ $\alpha \log N$ で増えます。データが増えるほどクラスタも対数的にゆっくり増える、という挙動がここから読み取れます。これは「データが増えたら新しいタイプが現れてほしい」という当初の要求に正確に応えています。
3つのパネルは $\alpha=0.5$・$2.0$・$10.0$ での CRP シミュレーション結果（$N=60$）です。$\alpha=0.5$ ではほぼ1〜2テーブルに人数が集中し、大きな少数クラスタが形成されます。$\alpha=2.0$ では中程度のテーブル数に落ち着きます。$\alpha=10.0$ では多数の小テーブルが生まれ、クラスタが細かく分かれます。いずれの場合も人数順に並べると急激に減少する「べき乗則的な」分布を示しており、CRP の富める者はさらに富む（rich-get-richer）性質が一目で分かります。
ここまでで CRP による割り当ての生成を理解しました。次に、混合比そのものを明示的に作る別の見方「スティック折り過程」を見ておきます。CRP が「割り当て」の視点なら、こちらは「混合比」の視点です。
スティック折り過程 — 混合比の構成的定義
CRP は割り当て $\bm{z}$ の周辺分布を与えますが、「無限個の混合比 $\pi_1, \pi_2, \dots$ をどう作るか」という問いには直接答えていません。これに構成的な答えを与えるのがスティック折り過程（stick-breaking process） です。混合比を目で見て作れるので、ディリクレ過程のもう一つの直感的な理解になります。
長さ1の棒（スティック）を用意します。これを次の手順で無限回折っていきます。
$$ \begin{align} \beta_k &\sim \mathrm{Beta}(1, \alpha) \quad (k = 1, 2, \dots) \\ \pi_k &= \beta_k \prod_{j=1}^{k-1}(1 – \beta_j) \end{align} $$
直感的には、棒の左から $\beta_1$ の割合を折り取って $\pi_1$ とし、残りの $(1-\beta_1)$ からさらに $\beta_2$ の割合を折り取って $\pi_2$ とし……を無限に繰り返します。$\prod_{j=1}^{k-1}(1-\beta_j)$ は「$k$ 番目までに残っている棒の長さ」を表します。この構成だと $\sum_{k=1}^\infty \pi_k = 1$ が確率1で成り立ちます（残りが幾何級数的に0へ近づくため）。
集中度 $\alpha$ の役割もここで明快になります。$\beta_k \sim \mathrm{Beta}(1, \alpha)$ の期待値は $1/(1+\alpha)$ です。$\alpha$ が小さいと毎回大きな割合を折り取るので、最初の数本の $\pi_k$ に重みが集中します（=少数の大きなクラスタ）。$\alpha$ が大きいと少しずつしか折り取らないので、重みが多くの成分に広く分散します（=多数の小さなクラスタ）。CRP で見た「$\alpha$ が大きいほどクラスタが増える」挙動と完全に整合しています。
各棒が成分番号 $k$ の混合比 $\pi_k$ を表しています。$\alpha=0.5$ では最初の1〜2成分が棒のほぼ全長を占めており、残りはほぼゼロです。$\alpha=2.0$ では重みが5〜6成分に広がり、$\alpha=8.0$ ではさらに多くの成分に均等に分散します。パネル上部に示した「上位3成分の和」（$\alpha=0.5$ で 0.95 程度、$\alpha=8.0$ で 0.45 程度）を見ると、$\alpha$ が大きいほど先頭の成分への集中が薄まることが定量的に確認できます。
この $\pi_k$ と、基底分布から独立に引いた成分パラメータ $\bm{\theta}_k \sim H$ を組み合わせると、
$$ G = \sum_{k=1}^{\infty} \pi_k\, \delta_{\bm{\theta}_k} $$
という離散確率測度ができあがります。ここで $\delta_{\bm{\theta}_k}$ は点 $\bm{\theta}_k$ に質量を持つディラックのデルタです。この $G$ がディリクレ過程 $\mathrm{DP}(\alpha, H)$ からのサンプルです。ディリクレ過程とは「確率分布 $G$ を生み出す分布」、すなわち分布の分布なのです。基底分布 $H$ が連続であっても、$G$ は離散になる（同じ $\bm{\theta}_k$ が複数のデータに共有される）という点が、クラスタリングを可能にする鍵です。
正式には、ディリクレ過程 $\mathrm{DP}(\alpha, H)$ は次の性質で定義されます。空間 $\Theta$ の任意の有限分割 $A_1, \dots, A_r$ に対して、$G$ がその分割に割り当てる質量がディリクレ分布に従う、というものです。
$$ (G(A_1), \dots, G(A_r)) \sim \mathrm{Dir}(\alpha H(A_1), \dots, \alpha H(A_r)) $$
これがまさに「ディリクレ分布の無限次元への拡張」という呼び名の由来です。どんな細かさで空間を区切っても、その上での質量分配がディリクレ分布になっている、という強い一貫性を持った無限次元オブジェクトなのです。
縦棒が $\mathrm{DP}(\alpha, H)$ のサンプル $G$ を表しており、縦の高さが各点 $\bm{\theta}_k$ の質量 $\pi_k$、横位置がその値 $\bm{\theta}_k$ です。橙の破線は基底分布 $H$（ガウス）です。$H$ は連続な密度関数なのに、$G$ は有限個の点への集中質量（縦棒）になっています。これが「DP のサンプルは離散」という重要な性質で、複数のデータ点が同じ $\bm{\theta}_k$ を共有しクラスタを形成する仕組みの根拠です。$\alpha$ が大きいほど棒の本数が増え、質量が広い範囲に散らばっているのも確認できます。
ここまでで CRP（割り当て）とスティック折り過程（混合比）の2つの視点からディリクレ過程を理解しました。次は、これを使った混合モデル DPM をギブスサンプリングで推論する具体的なアルゴリズムへ進みます。
ディリクレ過程混合モデルとギブスサンプリング
ディリクレ過程をそのまま観測データに当てはめると、各点が必ず別々の $\bm{\theta}$ を持ってしまい使いにくいです。そこで $\bm{\theta}$ からノイズを乗せて観測 $x_i$ を生成する一段階を加えたものが DPM です。生成過程は次の通りです。
$$ \begin{align} G &\sim \mathrm{DP}(\alpha, H) \\ \bm{\theta}_i \mid G &\sim G \\ x_i \mid \bm{\theta}_i &\sim F(\bm{\theta}_i) \end{align} $$
$G$ が離散なので複数の $\bm{\theta}_i$ が同じ値を取り、それが自然にクラスタを形成します。推論の目標は、観測 $\{x_i\}$ から割り当て $\bm{z}$ とクラスタ数を同時に推定することです。
ここで本記事の具体的なモデルを固定します。各クラスタは2次元ガウス分布で、共分散は既知の等方 $\sigma^2 \bm{I}$、平均 $\bm{\mu}_k$ が未知とします。基底分布 $H$ は平均に対するガウス事前 $\mathcal{N}(\bm{\mu}_0, \sigma_0^2 \bm{I})$ とします。これはガウス・ガウス共役なので、クラスタ平均を解析的に周辺化できます（collapsed Gibbs sampling）。
周辺化ギブスサンプラーのアルゴリズム
CRP の交換可能性のおかげで、各点 $i$ を「最後の客」とみなして再着席させられます。各反復で全点 $i$ について次を行います。
- 点 $i$ を現在の所属から外す。所属クラスタの人数 $n_{z_i}$ を1減らす。もしそのクラスタが空になったら削除する。
- 既存クラスタ $k$ への割り当て確率を計算する。CRP 事前 $\times$ 尤度で、 $$ p(z_i = k \mid \cdots) \propto n_{-i,k} \cdot p(x_i \mid \bm{x}_{-i,k}) $$ ここで $p(x_i \mid \bm{x}_{-i,k})$ は、クラスタ $k$ の既存データ $\bm{x}_{-i,k}$ を見たうえでの $x_i$ の事後予測分布です。
- 新クラスタへの割り当て確率を計算する。 $$ p(z_i = \text{新} \mid \cdots) \propto \alpha \cdot p(x_i \mid H) $$ ここで $p(x_i \mid H)$ は、まだ何も観測していない（事前分布 $H$ だけの）ときの $x_i$ の事前予測分布です。
- 上の確率を正規化し、カテゴリカル分布から $z_i$ をサンプリングして再着席させる。
ステップ2・3で尤度をかける点が CRP 単独との違いです。CRP は事前分布、尤度はデータが各クラスタにどれだけ「馴染むか」を表し、両者の積（=事後）でサンプリングするのがベイズの王道です。「人気のテーブル（$n$ 大）」かつ「料理が口に合うテーブル（尤度大）」に座りやすい、というわけです。
予測分布の具体形（ガウス・ガウス共役）
等方ガウス・ガウス共役の事後予測分布は解析的に書けます。クラスタ $k$ に既存点 $\bm{x}_{-i,k}$（$n$ 個、和 $\bm{s} = \sum \bm{x}$）があるとき、平均 $\bm{\mu}_k$ の事後は再びガウスになり、その事後予測分布も2次元ガウスになります。各次元独立なので1次元で書くと、
$$ p(x_i \mid \bm{x}_{-i,k}) = \mathcal{N}\!\left(x_i \;\middle|\; \mu_n,\; \sigma^2 + \sigma_n^2\right) $$
ここで事後平均と事後分散は
$$ \sigma_n^2 = \left(\frac{1}{\sigma_0^2} + \frac{n}{\sigma^2}\right)^{-1}, \qquad \mu_n = \sigma_n^2\left(\frac{\mu_0}{\sigma_0^2} + \frac{s}{\sigma^2}\right) $$
です。式を一つずつ読み解くと、事後分散 $\sigma_n^2$ は「事前の精度 $1/\sigma_0^2$ とデータの精度 $n/\sigma^2$ を足した逆数」、すなわち精度（情報量）が足し合わされる形になっています。事後平均 $\mu_n$ は事前平均とデータ和をそれぞれの精度で重み付けした加重平均です。予測分散に $\sigma^2$ が追加で乗るのは、新しい点 $x_i$ 自体に観測ノイズが乗るからです。新クラスタの場合は $n=0$ とすればよく、予測分布は $\mathcal{N}(\mu_0, \sigma^2 + \sigma_0^2)$ になります。
グレーの破線が「事前予測（データなし）」、緑の実線が「事後予測（データあり）」を表しています。観測数が $n=0$ のとき（新クラスタ）は広い予測分布になっており、$n=5\to20$ と増えるにつれてデータ平均（赤の点線）に引き寄せられて分布が鋭くなります。この「観測が増えるほど確信を深める」挙動が、ギブスサンプリングの各ステップで「データと馴染むクラスタ」を正しく選ぶ原動力になっています。
これでアルゴリズムの全ピースが揃いました。次に、この周辺化ギブスサンプラーを Python でスクラッチ実装し、真のクラスタ数を知らないまま2次元データを分類してみます。
Pythonでの実装
まず、クラスタ数を「知らないふり」をして推定するために、3つの真のクラスタを持つ2次元データを生成します。アルゴリズムにはクラスタ数 $K=3$ を一切教えません。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
# --- 真の3クラスタからデータ生成（クラスタ数はアルゴリズムに教えない）---
true_means = np.array([[0.0, 0.0], [6.0, 6.0], [-5.0, 5.0]])
n_per = [60, 50, 40] # 各クラスタの点数
sigma_true = 1.0 # 真の観測ノイズの標準偏差
X = np.vstack([
m + sigma_true * np.random.randn(n, 2)
for m, n in zip(true_means, n_per)
])
N = X.shape[0]
print(f"総データ点数 N = {N}（真のクラスタ数 = 3 だが推論には使わない）")
plt.figure(figsize=(6, 6))
plt.scatter(X[:, 0], X[:, 1], s=18, c='gray', alpha=0.7)
plt.title("Observed data (true #clusters hidden)")
plt.xlabel("x1"); plt.ylabel("x2")
plt.grid(alpha=0.3)
plt.show()
散布図には3つのまとまりがうっすら見えますが、境界付近の点はどちらに属するか曖昧で、人間の目でも「本当に3つか？」と断言はできません。アルゴリズムにはこの座標だけを渡し、クラスタ数も割り当ても一から推定させます。総点数は150点です。
次に、ガウス・ガウス共役の事後予測分布（対数尤度）を計算する補助関数を用意します。数値安定性のため対数で計算します。
def log_predictive(x, n, s, ss, sigma2, sigma0_2, mu0):
"""クラスタの事後予測分布の対数尤度（2次元・等方ガウス）
x : 評価する点 (2,)
n : クラスタ内の点数（新クラスタなら 0）
s : クラスタ内の点の和 (2,)
sigma2 : 観測ノイズ分散
sigma0_2: 事前分散
mu0 : 事前平均 (2,)
"""
# 事後分散と事後平均（次元独立なのでスカラー精度でOK）
post_var = 1.0 / (1.0 / sigma0_2 + n / sigma2)
post_mean = post_var * (mu0 / sigma0_2 + s / sigma2)
# 予測分散 = 事後分散 + 観測ノイズ
pred_var = post_var + sigma2
# 2次元独立ガウスの対数密度
d = x - post_mean
return np.sum(-0.5 * np.log(2 * np.pi * pred_var) - 0.5 * d**2 / pred_var)
この関数は1つの点 $x$ が、与えられたクラスタの統計量（点数 $n$、和 $s$）から見てどれだけ「馴染むか」を対数尤度として返します。n=0, s=0 を渡せば新クラスタ（事前分布だけ）の予測尤度になり、ステップ3にそのまま流用できます。引数 ss（二乗和）は今回の等方既知分散モデルでは使いませんが、分散も推定する拡張に備えて残しています。
続いて、CRP ギブスサンプリング本体を実装します。クラスタごとに「点数」と「点の和」だけを持てば十分（十分統計量）なので、それらを辞書で管理します。
def crp_gibbs(X, alpha, n_iter=50, sigma2=1.0, sigma0_2=25.0, seed=1):
"""ディリクレ過程混合（等方ガウス）の周辺化ギブスサンプリング"""
rng = np.random.default_rng(seed)
N, D = X.shape
mu0 = np.zeros(D)
# 初期化: 全点を1つのクラスタ0に割り当て
z = np.zeros(N, dtype=int)
counts = {0: N} # クラスタ -> 点数
sums = {0: X.sum(axis=0).copy()} # クラスタ -> 点の和
next_id = 1 # 新クラスタに振る次の番号
k_history = [] # 各反復後の使用クラスタ数を記録
for it in range(n_iter):
for i in range(N):
# --- (1) 点 i を現在の所属から外す ---
ci = z[i]
counts[ci] -= 1
sums[ci] -= X[i]
if counts[ci] == 0: # 空になったクラスタは削除
del counts[ci]; del sums[ci]
# --- (2)(3) 各候補クラスタの対数確率を計算 ---
keys = list(counts.keys())
logp = []
for k in keys: # 既存クラスタ: log(n_k) + log予測尤度
lp = np.log(counts[k]) + log_predictive(
X[i], counts[k], sums[k], None, sigma2, sigma0_2, mu0)
logp.append(lp)
# 新クラスタ: log(alpha) + 事前予測尤度（n=0, s=0）
logp.append(np.log(alpha) + log_predictive(
X[i], 0, np.zeros(D), None, sigma2, sigma0_2, mu0))
# --- (4) 正規化してサンプリング ---
logp = np.array(logp)
p = np.exp(logp - logp.max())
p /= p.sum()
choice = rng.choice(len(p), p=p)
if choice == len(keys): # 新クラスタが選ばれた
z[i] = next_id
counts[next_id] = 1
sums[next_id] = X[i].copy()
next_id += 1
else: # 既存クラスタが選ばれた
k = keys[choice]
z[i] = k
counts[k] += 1
sums[k] += X[i]
k_history.append(len(counts))
return z, counts, sums, k_history
実装のポイントは、クラスタを固定長配列ではなく辞書で持つことです。これにより、空になったクラスタを del で消し、新クラスタを動的に追加できます。クラスタ数が固定でないノンパラメトリックモデルならではの設計で、k-means の実装とは根本的に異なります。対数領域で確率を計算し logp - logp.max() で引き算してから exp するのは、アンダーフロー（極小値の桁落ち）を防ぐ定番テクニックです。
それでは集中度 $\alpha = 1.0$ で実行し、推定されたクラスタリング結果を可視化します。
z, counts, sums, k_hist = crp_gibbs(X, alpha=1.0, n_iter=50, sigma2=1.0)
# クラスタIDを 0,1,2,... に振り直す（表示用）
uniq = list(counts.keys())
remap = {k: i for i, k in enumerate(uniq)}
z_plot = np.array([remap[zi] for zi in z])
print(f"推定クラスタ数 = {len(uniq)} 各クラスタの点数 = {sorted(counts.values(), reverse=True)}")
plt.figure(figsize=(6, 6))
for i, k in enumerate(uniq):
pts = X[z == k]
plt.scatter(pts[:, 0], pts[:, 1], s=18, alpha=0.7, label=f"cluster {i} (n={counts[k]})")
plt.title(r"DPM clustering ($\alpha=1.0$)")
plt.xlabel("x1"); plt.ylabel("x2")
plt.legend(); plt.grid(alpha=0.3)
plt.show()
左が「参考用の真のラベル」、右が「クラスタ数を一切教えずに DPM が推定したラベル」です。推定クラスタ数はおおむね3に落ち着き、3つのまとまりが正しく色分けされています。黒い星印は各クラスタの推定中心（事後平均から計算）であり、真の中心に近い位置に来ています。ときどき境界付近に少数の点からなる小さな余剰クラスタが生じることもありますが、点数の多い上位3クラスタが真の構造を捉えていることが点数の表示（例: [60, 50, 40] に近い値）から読み取れます。クラスタ数 $K$ をハイパーパラメータとして与えずに、データ自身が「3つで十分」と語っている点が DPM の核心です。
次に、本記事の主題である集中度 $\alpha$ が推定クラスタ数をどう左右するかを体系的に調べます。
# --- α を変えて推定クラスタ数の変化を見る ---
alphas = [0.05, 0.2, 1.0, 5.0, 20.0]
est_K = []
for a in alphas:
_, counts_a, _, _ = crp_gibbs(X, alpha=a, n_iter=50, sigma2=1.0, seed=2)
# 点数2以上の「実質的な」クラスタ数を数える
K_a = sum(1 for c in counts_a.values() if c >= 2)
est_K.append(K_a)
print(f"alpha = {a:5.2f} -> 推定クラスタ数 = {K_a}")
plt.figure(figsize=(7, 4.5))
plt.semilogx(alphas, est_K, 'o-', color='C3')
plt.axhline(3, color='gray', ls='--', label='true K = 3')
plt.xlabel(r"concentration $\alpha$ (log scale)")
plt.ylabel("estimated #clusters")
plt.title("Effect of concentration parameter on #clusters")
plt.legend(); plt.grid(alpha=0.3)
plt.show()
赤丸が推定クラスタ数（点数2以上の実質クラスタ）、青四角が理論値 $\alpha\log N$、灰色の点線が真の値 $K=3$ です。このグラフから、集中度 $\alpha$ と推定クラスタ数の関係が明確に読み取れます。$\alpha$ が小さい（0.05程度）と新テーブルが嫌われるため、3つの真のクラスタが1〜2個に過剰に併合される傾向が出ます。$\alpha$ が中程度（0.2〜1.0）では真の値3を安定して当てます。$\alpha$ が大きい（5〜20）と新テーブルへの引力が強まり、4個以上に過分割される傾向が現れます。理論値 $\alpha\log N$ とも傾向が整合しており、CRP のクラスタ数期待値の式が実際の推定挙動をよく説明していることが分かります。$\alpha$ は「クラスタの細かさ」を制御するつまみだと理解できます。
最後に、ギブスサンプリングの反復に伴ってクラスタ数がどう収束していくかを追います。初期状態（全点1クラスタ）からどのくらいで安定するかを見ます。
# --- 反復に伴うクラスタ数の収束（複数の α で比較）---
plt.figure(figsize=(8, 4.5))
for a, col in zip([0.2, 1.0, 5.0], ['C0', 'C1', 'C2']):
_, _, _, k_hist = crp_gibbs(X, alpha=a, n_iter=60, sigma2=1.0, seed=7)
plt.plot(k_hist, color=col, label=fr"$\alpha={a}$")
plt.axhline(3, color='gray', ls='--', label='true K = 3')
plt.xlabel("Gibbs iteration")
plt.ylabel("#active clusters")
plt.title("Convergence of #clusters over Gibbs iterations")
plt.legend(); plt.grid(alpha=0.3)
plt.show()
この収束グラフからは2つのことが読み取れます。第一に、初期は全点が1クラスタ（クラスタ数1）から始まりますが、数回の反復で急速にクラスタが分裂し、$\alpha=1.0$（破線）では真の値3付近で安定します。これはギブスサンプリングがマルコフ連鎖として目標事後分布に近づいていく様子そのものです。第二に、$\alpha=5.0$（点線）では定常状態のクラスタ数が高い水準で揺らぎ、$\alpha=0.2$（実線）では低い水準に収まります。定常後もクラスタ数が完全には一定にならず小刻みに揺れているのは、これがサンプリング（事後分布からのドロー）であって最適化ではないためで、各反復のサンプルは事後分布上の1標本だと解釈すべきことを示しています。実用上は、十分なバーンイン後の複数サンプルを集計して事後要約とします。
まとめ
本記事では、ディリクレ過程混合モデルと中華料理店過程の理論を導出から実装まで一気通貫で解説しました。
- 有限混合からの極限: 対称ディリクレ事前 $\mathrm{Dir}(\alpha/K, \dots)$ を置いた有限混合で混合比を周辺化し、$K \to \infty$ の極限を取ると、$\alpha/K$ という設計のおかげで $K$ が約分されて消え、有限の予測分布が残る
- 中華料理店過程（CRP）: 既存テーブルには人数 $n_{-i,k}$ に比例し、新テーブルには集中度 $\alpha$ に比例して座る予測分布。交換可能性により各点を「最後の客」とみなせ、ギブスサンプリングが正当化される
- スティック折り過程: 棒を $\mathrm{Beta}(1,\alpha)$ で次々に折って混合比 $\pi_k$ を構成し、離散測度 $G = \sum_k \pi_k \delta_{\theta_k}$ としてディリクレ過程のサンプルを得る、混合比視点の構成的定義
- DPM のギブスサンプリング: CRP 事前 $\times$ ガウス・ガウス共役の予測尤度でサンプリングし、クラスタ数を固定せずに辞書で動的管理。Python 実装で真の3クラスタを当てられることを確認
- 集中度 $\alpha$ の役割: $\alpha$ が小さいと少数の大クラスタ、大きいと多数の小クラスタ。推定クラスタ数は $\alpha$ と単調に増減し、理論の $\alpha \log N$ と整合
ディリクレ過程混合モデルは、「クラスタ数を事前に決めなくてよい」という強力な性質によって、データの構造を柔軟に捉えます。次のステップとして、以下の記事も参考にしてください。
- ガウス混合モデル（GMM）とEMアルゴリズム — クラスタ数を固定する有限混合の標準手法。DPM との比較で理解が深まる
- 変分推論（変分ベイズ）入門 — DPM をサンプリングではなく変分法で高速に近似推論する手法
- ノンパラメトリックベイズ入門 — ディリクレ過程を含むノンパラメトリックベイズ全体の地図
参考文献
- C. E. Rasmussen, “The Infinite Gaussian Mixture Model”, NIPS 2000
- R. M. Neal, “Markov Chain Sampling Methods for Dirichlet Process Mixture Models”, Journal of Computational and Graphical Statistics, 2000
- Y. W. Teh et al., “Hierarchical Dirichlet Processes”, JASA, 2006
- パターン認識と機械学習（PRML）