ニュース記事を1000本集めて、それぞれを「政治」「スポーツ」「経済」「エンタメ」というテーマの混合比率で表したい——こうしたタスクは自然言語処理の現場でよく登場します。1本目の記事は「政治70%・経済30%」、2本目は「スポーツ50%・エンタメ50%」のように、各記事が複数のテーマの確率的な混合になっているわけです。ところが、その「混合比率」自体もデータごとに違うばらつきを持ちます。ある記事は政治寄り、別の記事はスポーツ寄り——この比率そのものの不確実性を、確率分布として表現したい。
このとき活躍するのがディリクレ分布（Dirichlet distribution）です。比率は和が1になり各要素が0以上1以下、つまり「確率の確率」を扱う場面で常に顔を出します。LDA（Latent Dirichlet Allocation）というトピックモデルでは文書のトピック比率の事前分布、ベイズ推定では多項分布の共役事前分布、ガウス混合モデルでは混合比の事前分布、A/Bテストの多変量版（多腕バンディット問題）では各選択肢の選択確率の不確実性表現——応用先は驚くほど広いです。
ディリクレ分布は「ベータ分布の多次元拡張」と表現されることが多く、これは正しい説明ですが、初学者がその一言で腑に落ちることはまずありません。本記事では、6面サイコロの目の出る確率というシンプルな例から始め、密度関数の定義、期待値・分散・モードの導出、多項分布の共役事前分布としての性質、そしてLDAトピックモデルへの応用まで、Pythonの可視化と実装を交えて段階的に解きほぐしていきます。
本記事の内容
- ディリクレ分布の直感（シンプレックス上の分布）と確率密度関数の定義
- 多次元ベータ関数による正規化定数の導出
- 期待値・分散・モードの導出と意味
- 多項分布の共役事前分布としての性質（事後更新の計算）
- LDAトピックモデルでの役割と簡易LDAのPython実装
scipy.stats.dirichletによる可視化とサンプリング- ディリクレ過程による非パラメトリックベイズへの拡張
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ベータ分布とは？定義や性質をわかりやすく解説
- カテゴリカル分布とは？多項分布の前段としての位置づけ
- 多項分布の定義と性質
- ベルヌーイ分布のベイズ推定（事後分布がベータになる例）
- ガンマ関数の定義と性質
ベータ分布が「2択（コインの表/裏）の比率 $p$ の分布」だったことを思い出しておくと、ディリクレ分布が「$K$ 択の比率ベクトル $\bm{\pi}$ の分布」という多次元拡張だという見通しが立ちやすくなります。
直感 — シンプレックス上の分布
「比率を確率変数とする」とはどういうことか
ディリクレ分布の最大のポイントは、確率変数自体が確率分布のパラメータである点です。普通の確率分布（正規分布、ポアソン分布など）は実数や整数を確率変数にしますが、ディリクレ分布の確率変数は「合計が1になる非負ベクトル」、つまり確率分布そのものです。
具体例で考えてみましょう。6面のイカサマサイコロがあるとします。普通のサイコロなら各目の出る確率は $1/6$ ずつですが、このサイコロは歪んでいて、各目の確率は $\bm{\pi} = (\pi_1, \pi_2, \pi_3, \pi_4, \pi_5, \pi_6)$ という未知のベクトルです。当然、$\pi_k \geq 0$ かつ $\sum_{k=1}^{6} \pi_k = 1$ を満たします。
このサイコロを「製造ロットからランダムに1個取り出した」と考えると、$\bm{\pi}$ 自体もロットごとにばらつく確率変数です。ある個体は1の目が出やすく、別の個体は3の目が出やすい。この $\bm{\pi}$ の分布を表すのがディリクレ分布なのです。
シンプレックス
「合計が1の非負ベクトル」の集合を確率シンプレックス（probability simplex）と呼びます。3次元（$K=3$）の場合、シンプレックスは3次元空間内の正三角形の内部になります。
3つの頂点 $(1, 0, 0)$、$(0, 1, 0)$、$(0, 0, 1)$ を結ぶ三角形を考えてください。三角形内部の任意の点 $(\pi_1, \pi_2, \pi_3)$ は $\pi_1 + \pi_2 + \pi_3 = 1$ を満たし、しかも各 $\pi_k \geq 0$ となります。逆にこの条件を満たす点は必ず三角形内にあります。つまり3次元のディリクレ分布は、この正三角形の上に定義された連続確率分布として可視化できます。
上の図は $\bm{\alpha}$ を切り替えたときのシンプレックス上の密度の様子です。$(1,1,1)$ では全域が一様、$(5,5,5)$ では中心に山が立ち、$(10,1,1)$ では1つの頂点に偏り、$(0.5,0.5,0.5)$ では密度が角（頂点）に集中する——これだけで「$\bm{\alpha}$ のパターンが分布形状を完全に制御する」というディリクレ分布の核心が一目でつかめます。
このシンプレックスの中で、ディリクレ分布のパラメータ $\bm{\alpha}$ が「分布の山の位置と尖り具合」を決めます。$\bm{\alpha}$ の各成分が大きければ尖り、小さければ広がる——ベータ分布における $\alpha, \beta$ と同じ役割です。
ベータ分布との対応関係
$K = 2$ の場合、シンプレックスは「$\pi_1 + \pi_2 = 1$」を満たす線分になります。$\pi_2 = 1 – \pi_1$ なので、自由度は1つだけ。このとき分布は $\pi_1 \in [0, 1]$ 上の1次元分布になり、これはまさにベータ分布です。
ベータ分布 $\mathrm{Beta}(\alpha, \beta)$ がコインの表が出る確率 $p$ の分布を表すように、ディリクレ分布 $\mathrm{Dir}(\bm{\alpha})$ は $K$ 面サイコロの各目の確率ベクトル $\bm{\pi}$ の分布を表します。ベルヌーイ分布→カテゴリカル分布の関係（2択→$K$択）が、ベータ分布→ディリクレ分布の関係（2比率→$K$比率）と完全に対応している、と覚えておくと整理しやすいです。
この比較図では、$K=2$ で1次元のベータ分布、$K=3$ で三角形シンプレックス上の2次元分布、$K=4$ で4成分の周辺平均バーが並べて示されています。次元が上がるごとに「自由度1個ずつ増える」（$K-1$ 次元のシンプレックス）という関係も視覚的に確認できます。
シンプレックスの直感が手に入ったところで、次は密度関数の正確な定義に進みます。
ディリクレ分布の定義と多次元ベータ
確率密度関数
$K$ 次元のディリクレ分布 $\mathrm{Dir}(\bm{\pi} \mid \bm{\alpha})$ は、確率変数 $\bm{\pi} = (\pi_1, \dots, \pi_K)^\top$ がシンプレックス $\sum_k \pi_k = 1,\ \pi_k \geq 0$ 上にあるとき、次の確率密度関数で定義されます。
$$ \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha}) = \frac{1}{B(\bm{\alpha})} \prod_{k=1}^{K} \pi_k^{\alpha_k – 1} $$
ここで $\bm{\alpha} = (\alpha_1, \dots, \alpha_K)^\top$ は分布のパラメータで、各 $\alpha_k > 0$ です。式の中心は $\prod_k \pi_k^{\alpha_k – 1}$ という積で、これは「各 $\pi_k$ を $(\alpha_k – 1)$ 乗してから掛け合わせる」操作。ベータ分布の $p^{\alpha-1}(1-p)^{\beta-1}$ を素直に多次元化した形です。
$1/B(\bm{\alpha})$ は分布全体の積分が1になるようにする正規化定数で、多次元ベータ関数と呼ばれます。
$$ B(\bm{\alpha}) = \frac{\prod_{k=1}^{K} \Gamma(\alpha_k)}{\Gamma\left(\sum_{k=1}^{K} \alpha_k\right)} = \frac{\prod_{k=1}^{K} \Gamma(\alpha_k)}{\Gamma(\alpha_0)} $$
ここで $\alpha_0 = \sum_{k=1}^K \alpha_k$ と置きました。$\Gamma$ はガンマ関数で、$\Gamma(n) = (n-1)!$（正の整数のとき）を実数に拡張したものです。$K=2$ のとき $B(\alpha_1, \alpha_2) = \Gamma(\alpha_1)\Gamma(\alpha_2)/\Gamma(\alpha_1 + \alpha_2)$ となり、これは見覚えのあるベータ関数そのもの。多次元ベータ関数の名前の由来です。
正規化定数の導出
正規化定数が上の形になることは、ディリクレ積分とよばれる積分公式から導けます。
$$ \int_{\Delta_K} \prod_{k=1}^{K} \pi_k^{\alpha_k – 1}\, d\bm{\pi} = \frac{\prod_{k=1}^{K} \Gamma(\alpha_k)}{\Gamma(\alpha_0)} $$
ここで $\Delta_K$ は $K$ 次元シンプレックスを表します。この積分は、ガンマ関数の積分表示 $\Gamma(\alpha) = \int_0^\infty t^{\alpha-1} e^{-t} dt$ を使うと示せます。$K$ 個のガンマ関数を掛け合わせたものから出発して、$y_k = t \pi_k$ という変数変換を行います。
$$ \prod_{k=1}^K \Gamma(\alpha_k) = \prod_{k=1}^K \int_0^\infty y_k^{\alpha_k-1} e^{-y_k}\, dy_k = \int_{[0,\infty)^K} \prod_k y_k^{\alpha_k-1} e^{-\sum_k y_k}\, d\bm{y} $$
ここで $t = \sum_k y_k$、$\pi_k = y_k / t$ と置き直します。すると $y_k = t \pi_k$、$dy_k = t\, d\pi_k$ となり、$\bm{y}$ 空間から $(t, \bm{\pi})$ 空間への変換のヤコビアンを計算する必要があります。シンプレックス制約 $\sum_k \pi_k = 1$ を考慮すると、ヤコビアン $|J| = t^{K-1}$ となることが知られています。
$$ \prod_k \Gamma(\alpha_k) = \int_0^\infty \int_{\Delta_K} \prod_k (t\pi_k)^{\alpha_k – 1} e^{-t} t^{K-1}\, d\bm{\pi}\, dt $$
$t$ について整理すると $t^{\sum_k \alpha_k – K + K – 1} = t^{\alpha_0 – 1}$ となるので、
$$ \prod_k \Gamma(\alpha_k) = \int_0^\infty t^{\alpha_0 – 1} e^{-t}\, dt \cdot \int_{\Delta_K} \prod_k \pi_k^{\alpha_k – 1}\, d\bm{\pi} = \Gamma(\alpha_0) \int_{\Delta_K} \prod_k \pi_k^{\alpha_k – 1}\, d\bm{\pi} $$
両辺を $\Gamma(\alpha_0)$ で割ればディリクレ積分の公式が出てきます。この公式により、正規化定数として $\Gamma(\alpha_0)/\prod_k \Gamma(\alpha_k) = 1/B(\bm{\alpha})$ を選べば、$\mathrm{Dir}(\bm{\pi}|\bm{\alpha})$ の積分が1になることが保証されます。
$\bm{\alpha}$ の解釈 — 「観測されたカウントの先取り」
パラメータ $\bm{\alpha}$ には、ベイズ的に分かりやすい解釈があります。後ほど共役事前分布の節で詳しく見ますが、$\alpha_k$ は「事前に $k$ 番目のカテゴリが何回観測されたと仮定するか」を表す疑似カウント（pseudocount）として読めます。$\bm{\alpha} = (1, 1, \dots, 1)$ なら「事前知識ゼロ」、$\bm{\alpha} = (10, 1, 1, \dots, 1)$ なら「1番目の目が出やすそうだと事前に強く信じている」というニュアンスです。
また、$\alpha_0 = \sum_k \alpha_k$ は分布の集中度（concentration）を制御します。$\alpha_0$ が大きいほど分布は鋭く尖り、小さいほど広がります。$\bm{\alpha}/\alpha_0$ が分布の中心（平均比率）を決め、$\alpha_0$ の大きさが「その中心への自信」を決める、と分けて考えると整理しやすいでしょう。
上の4枚はちょうどこの「平均は同じだが自信が変わる」状況を可視化したものです。平均位置は4枚とも★印の場所（$(0.6, 0.3, 0.1)$）で共通ですが、$\alpha_0=1$ では分布は広く散り、$\alpha_0=100$ では★の周りに鋭く集中します。サンプルサイズ依存の事後分布のシャープ化（次節以降）も、この $\alpha_0$ 増加と同じ機構で説明できます。
定義式が分かったところで、次は分布の平均や分散などの基本統計量を導いていきます。
期待値・分散・モード
ディリクレ分布のパラメータと、分布の形状の関係を定量化するため、基本統計量を求めていきます。
期待値
$\pi_k$ の期待値は次のシンプルな形になります。
$$ \mathbb{E}[\pi_k] = \frac{\alpha_k}{\alpha_0} $$
つまり「各 $\alpha_k$ を全体で割った比率」がそのまま平均比率です。サイコロの例なら $\bm{\alpha} = (\alpha_1, \dots, \alpha_6)$ が与えられたとき、$k$ の目が出る確率の事前期待値は $\alpha_k/\sum_j \alpha_j$ で評価できます。
導出は次のように進めます。期待値の定義から
$$ \mathbb{E}[\pi_k] = \int_{\Delta_K} \pi_k \cdot \frac{1}{B(\bm{\alpha})} \prod_{j=1}^K \pi_j^{\alpha_j – 1}\, d\bm{\pi} $$
$\pi_k$ を $\pi_k^{\alpha_k – 1}$ に吸収させると、被積分関数は $\pi_k^{\alpha_k}$ になります。これは $\bm{\alpha}’ = (\alpha_1, \dots, \alpha_k + 1, \dots, \alpha_K)$ をパラメータとする別のディリクレ分布の被積分関数の形と一致します。$\bm{\alpha}’$ をパラメータとするディリクレ分布の積分は1なので、
$$ \int_{\Delta_K} \prod_j \pi_j^{\alpha_j’ – 1}\, d\bm{\pi} = B(\bm{\alpha}’) = \frac{\Gamma(\alpha_k + 1)\prod_{j\neq k}\Gamma(\alpha_j)}{\Gamma(\alpha_0 + 1)} $$
これを使うと
$$ \mathbb{E}[\pi_k] = \frac{B(\bm{\alpha}’)}{B(\bm{\alpha})} = \frac{\Gamma(\alpha_k + 1)/\Gamma(\alpha_k)}{\Gamma(\alpha_0 + 1)/\Gamma(\alpha_0)} = \frac{\alpha_k}{\alpha_0} $$
ガンマ関数の漸化式 $\Gamma(x+1) = x\Gamma(x)$ を最後に使いました。結果として、各成分の平均はパラメータの比率そのものになるという、非常に解釈しやすい式が得られます。
分散と共分散
分散と共分散も同じ手法で計算できます。結果だけ示すと、
$$ \mathrm{Var}[\pi_k] = \frac{\alpha_k (\alpha_0 – \alpha_k)}{\alpha_0^2 (\alpha_0 + 1)}, \quad \mathrm{Cov}[\pi_i, \pi_j] = -\frac{\alpha_i \alpha_j}{\alpha_0^2 (\alpha_0 + 1)} \quad (i \neq j) $$
分散の形を観察すると、$\alpha_0$ が大きいほど分母が大きく、分散は小さくなることが読み取れます。これは「疑似カウントが多いほど、平均値への自信が強い」という直感と合致します。1万回サイコロを振ったあとの確率の推定値は、10回振ったあとよりずっと信頼できる、という感覚です。
共分散が負になる点も重要です。$\pi_i$ が大きくなれば、$\sum_k \pi_k = 1$ の制約のため他の $\pi_j$ は小さくならざるを得ません。シンプレックス上の分布である以上、成分間に負の相関が生じるのは必然です。
モード（最頻値）
ディリクレ分布のモード（密度関数を最大にする点）は、$\alpha_k > 1$ がすべての $k$ について成り立つとき、次の形で与えられます。
$$ \mathrm{mode}[\pi_k] = \frac{\alpha_k – 1}{\alpha_0 – K} $$
導出はラグランジュ未定乗数法で行います。$\log p(\bm{\pi}) = \sum_k (\alpha_k – 1)\log \pi_k + \text{const}$ を制約 $\sum_k \pi_k = 1$ のもとで最大化するため、ラグランジアン
$$ \mathcal{L} = \sum_k (\alpha_k – 1)\log \pi_k – \lambda\left(\sum_k \pi_k – 1\right) $$
の $\pi_k$ に関する偏微分をゼロにします。
$$ \frac{\partial \mathcal{L}}{\partial \pi_k} = \frac{\alpha_k – 1}{\pi_k} – \lambda = 0 \implies \pi_k = \frac{\alpha_k – 1}{\lambda} $$
制約 $\sum_k \pi_k = 1$ より $\lambda = \sum_k (\alpha_k – 1) = \alpha_0 – K$ となり、$\pi_k = (\alpha_k – 1)/(\alpha_0 – K)$ が得られます。
ここで重要な注意があります。すべての $\alpha_k > 1$ なら内部点でモードを持ちますが、ある $\alpha_k < 1$ がある場合、密度はその成分の境界（$\pi_k \to 0$）で発散し、分布は「角に押し付けられる」形になります。$\bm{\alpha} = (0.5, 0.5, 0.5)$ のような小さい値は、各カテゴリのどれか1つに極端に偏った比率を生成するという、特殊で有用な性質を持ちます（後述のLDAでスパース性を出すのに使われます）。
実際にディリクレ分布から多数のサンプルを取り出してシンプレックス上に散布してみると、$\bm{\alpha}=(0.3,0.3,0.3)$ のケースで点が3つの頂点に張りつく様子がはっきり見えます。これが「ディリクレで角集中＝スパースな比率を生成する」という言い回しの実体です。一方、$\bm{\alpha}=(10,10,10)$ では中央に集まり、$\bm{\alpha}=(8,2,2)$ では片側の頂点に偏った雲ができます。
統計量の式が揃ったところで、ディリクレ分布が真価を発揮する場面——多項分布との共役性——に進みましょう。
多項分布の共役事前分布
共役事前分布とは
ベイズ推定では、観測データ $\bm{x}$ が与えられたときのパラメータ $\bm{\pi}$ の事後分布を
$$ p(\bm{\pi} \mid \bm{x}) \propto p(\bm{x} \mid \bm{\pi})\, p(\bm{\pi}) $$
として求めます。ここで尤度 $p(\bm{x}|\bm{\pi})$ と事前分布 $p(\bm{\pi})$ を掛け合わせたものが事後分布と同じ族（family）になるとき、その事前分布を共役事前分布（conjugate prior）と呼びます。共役性があると事後分布の計算が解析的に閉じる、つまり「掛け算しても新しい関数形が出てこない」ので、ベイズ推定が劇的にシンプルになります。
ディリクレ分布は多項分布（およびカテゴリカル分布）の共役事前分布です。これを確認しましょう。
上のフロー図に示したように、事前分布 $\mathrm{Dir}(\bm{\alpha})$ と多項尤度 $\mathrm{Mult}(\bm{n}|\bm{\pi})$ を掛けた結果がそのままディリクレ分布の形 $\mathrm{Dir}(\bm{\alpha}+\bm{n})$ に戻る——これが共役性の正体です。掛け算で族が閉じるため、事後計算が「パラメータに観測カウントを足すだけ」のシンプルな操作に帰着します。
事後分布の導出
サイコロを $N$ 回振り、$k$ の目が $n_k$ 回出たとします（$\sum_k n_k = N$）。このときの尤度（カテゴリカル分布を $N$ 回繰り返した多項分布の形）は
$$ p(\bm{n} \mid \bm{\pi}) = \frac{N!}{n_1! n_2! \cdots n_K!} \prod_{k=1}^K \pi_k^{n_k} $$
事前分布をディリクレ分布 $\mathrm{Dir}(\bm{\pi}|\bm{\alpha})$ とすると、事後分布は
$$ p(\bm{\pi} \mid \bm{n}) \propto p(\bm{n} \mid \bm{\pi}) \cdot p(\bm{\pi} \mid \bm{\alpha}) \propto \prod_k \pi_k^{n_k} \cdot \prod_k \pi_k^{\alpha_k – 1} = \prod_k \pi_k^{\alpha_k + n_k – 1} $$
この関数形を見ると、これは新しいディリクレ分布の核そのものです。具体的には
$$ p(\bm{\pi} \mid \bm{n}) = \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha} + \bm{n}) $$
つまり「事前パラメータ $\bm{\alpha}$ に観測カウント $\bm{n}$ を足すだけ」で事後分布が得られます。これがディリクレ–多項共役性の核心です。
この結果は、先ほど触れた「$\alpha_k$ は疑似カウントである」という解釈にきれいに対応します。事前に「$k$ の目は $\alpha_k$ 回出ると思っていた」という信念を持ち、実際に $n_k$ 回観測したら、合計 $\alpha_k + n_k$ 回観測したことになる——非常に自然な更新です。
3枚並べた図を左から右へたどると、ベイズ更新の動きが具体的に見えます。①は事前分布の一様な青いシンプレックス、②は観測カウントの棒グラフ（赤）、③は事後分布のヒートマップ（緑）で、真値の★が事後の山の頂上付近にきれいに乗っていることが分かります。事前情報なしから出発しても、$N=30$ 程度の観測で真の比率を捉えられているのが視覚的に確認できます。
事後平均と最尤推定の関係
事後平均は
$$ \mathbb{E}[\pi_k \mid \bm{n}] = \frac{\alpha_k + n_k}{\alpha_0 + N} $$
これを次のように分解できます。
$$ \mathbb{E}[\pi_k \mid \bm{n}] = \underbrace{\frac{\alpha_0}{\alpha_0 + N}}_{w_{\text{prior}}} \cdot \underbrace{\frac{\alpha_k}{\alpha_0}}_{\text{事前平均}} + \underbrace{\frac{N}{\alpha_0 + N}}_{w_{\text{data}}} \cdot \underbrace{\frac{n_k}{N}}_{\text{MLE}} $$
事後平均は事前平均と最尤推定値の重み付き平均として書けます。重みは $\alpha_0$（事前の自信）と $N$（データ量）で決まり、データが増えると ($N \gg \alpha_0$) MLEに、データが少ないと ($N \ll \alpha_0$) 事前平均に近づきます。これは「データが多ければデータを信じ、少なければ事前分布を信じる」というベイズ推定の挙動を、定量的に表現したものです。
実用的には、$\bm{\alpha} = (1, 1, \dots, 1)$ のような一様な事前を選ぶと、事後平均は
$$ \mathbb{E}[\pi_k \mid \bm{n}] = \frac{n_k + 1}{N + K} $$
となり、これはラプラス平滑化（Laplace smoothing、加算スムージング）そのものです。ナイーブベイズ分類器や $n$-gram言語モデルで「観測されなかったカテゴリの確率を0にしない」工夫としておなじみのテクニックですが、その正体はディリクレ事前を使ったベイズ推定の事後平均だったわけです。
上の図の左側は、5カテゴリのうち $k=4$ が未観測（$n_4=0$）の場合の確率推定値を比較したものです。MLE（灰色）は $k=4$ の確率を「0」と推定してしまいますが、$\alpha=1$ の事後平均（青）はこれを正の値に持ち上げます。右側にはこの関係を式で示し、$\alpha_k=1$ を代入すると有名な「$(n_k+1)/(N+K)$」の加算スムージング公式が現れることを明示しました。ラプラス平滑化は経験則ではなく、ディリクレ事前を仮定したベイズ推定の自然な帰結なのです。
共役性を理解したところで、ディリクレ分布が大活躍する代表的な応用——LDAトピックモデル——に進みます。
LDA（トピックモデル）への応用
トピックモデルの問題設定
大量の文書（ニュース記事、論文、ツイートなど）が手元にあり、各文書がどんなテーマ（トピック）を扱っているかを自動的に発見したい——これがトピックモデリングの目的です。LDA（Latent Dirichlet Allocation、潜在ディリクレ配分法）は2003年にBlei、Ng、Jordanが提案した、この問題に対する代表的なベイズ生成モデルです。
LDAの基本仮定は次のとおりです。
- 各文書は複数のトピックの混合からなる（例：「政治60%・経済30%・スポーツ10%」）
- 各トピックは単語上の確率分布として表される（例：「政治」トピックでは「選挙」「議員」「政策」が高確率）
- 文書中の各単語は、まずトピックを選び、そのトピックの単語分布から単語をサンプルすることで生成される
この生成過程で、文書ごとのトピック比率とトピックごとの単語比率の両方がディリクレ分布から生成されると仮定するのがLDAの特徴です。両方とも「合計1の比率ベクトル」なので、ディリクレが自然な選択になります。
生成過程
形式的にLDAの生成過程は次のように書けます。
- 各トピック $k = 1, \dots, K$ について：単語分布 $\bm{\varphi}_k \sim \mathrm{Dir}(\bm{\beta})$ をサンプル
- 各文書 $d = 1, \dots, D$ について：
- トピック比率 $\bm{\theta}_d \sim \mathrm{Dir}(\bm{\alpha})$ をサンプル
- 文書 $d$ 中の各単語位置 $n = 1, \dots, N_d$ について：
- トピック $z_{d,n} \sim \mathrm{Categorical}(\bm{\theta}_d)$ をサンプル
- 単語 $w_{d,n} \sim \mathrm{Categorical}(\bm{\varphi}_{z_{d,n}})$ をサンプル
ここで $\bm{\alpha}$ は文書のトピック比率の事前パラメータ（次元 $K$）、$\bm{\beta}$ はトピックの単語比率の事前パラメータ（次元 $V$＝語彙サイズ）です。
上のプレート表記図に LDA の依存関係をまとめました。○が潜在変数、●（灰塗り）が観測変数、□がハイパーパラメータです。文書 $d$ ごとにディリクレからトピック比率 $\theta_d$ をサンプル、各トピック $k$ ごとに同じくディリクレから単語分布 $\varphi_k$ をサンプル、その上で各単語 $w_{d,n}$ が「$\theta_d$ → $z_{d,n}$ → $\varphi_{z_{d,n}}$」と二段階のカテゴリカル分布を経て生成される、という階層構造が一目で読み取れます。
ディリクレ事前のスパース効果
ここでパラメータの大きさが大きな意味を持ちます。$\bm{\alpha}$ の要素を1より小さく（例えば $\alpha_k = 0.1$）すると、ディリクレ分布はシンプレックスの角に集中し、トピック比率がスパース——つまり「ほとんどの文書は少数のトピックしか含まない」という性質を持つようになります。これはトピックモデルとして極めて望ましい性質です。実際の文書は「政治・経済・スポーツ・エンタメ……」のすべてを均等に含むより、2〜3トピックに偏ることがほとんどだからです。
同様に $\bm{\beta}$ を1より小さく設定すると、各トピックも語彙の中の限定された単語に偏るようになり、解釈可能性が高まります。$\bm{\alpha}, \bm{\beta}$ の大小を制御するだけで、モデルの挙動を大きく変えられるのがLDAの強みの一つです。
上のヒートマップは、ディリクレ事前パラメータ $\alpha$ の大きさが文書のトピック比率に与える効果を、$D=8$ 文書 × $K=10$ トピックで比較したものです。左の $\alpha=0.3$ では各文書が2〜3個のトピックだけに濃く色付き、それ以外はほぼ白——典型的なスパース構造です。右の $\alpha=5$ では全トピックに薄く色が広がり、どの文書も全テーマを均等に含む不自然な状態になっています。「実文書は少数トピックに偏る」という現実を再現するには、左のようにディリクレ事前を $\alpha < 1$ に設定するのが鍵です。
推論アルゴリズム（概要）
LDAの推論では、各単語に割り当てられたトピック $z_{d,n}$ の事後分布を求めます。よく使われる手法は次の3つです。
- 崩壊型ギブスサンプリング（collapsed Gibbs sampling）：$\bm{\theta}, \bm{\varphi}$ を解析的に積分消去し、$z_{d,n}$ のみをMCMCでサンプリング
- 変分ベイズ（variational Bayes）：事後分布を独立な分布の積で近似してKL最小化
- オンライン変分ベイズ：ストリーミングデータ向けの逐次更新版
いずれも、ディリクレ–カテゴリカル共役性を利用して計算を簡略化しています。たとえば崩壊型ギブスでは、各単語のトピック割当ての条件付き確率が次のようなコンパクトな形になります。
$$ p(z_{d,n} = k \mid \text{rest}) \propto \frac{n_{d,k}^{\neg(d,n)} + \alpha_k}{\sum_{k’} (n_{d,k’}^{\neg(d,n)} + \alpha_{k’})} \cdot \frac{m_{k, w_{d,n}}^{\neg(d,n)} + \beta_{w_{d,n}}}{\sum_v (m_{k,v}^{\neg(d,n)} + \beta_v)} $$
ここで $n_{d,k}^{\neg(d,n)}$ は文書 $d$ で現在の単語を除いてトピック $k$ に割り当てられた単語数、$m_{k,v}^{\neg(d,n)}$ はトピック $k$ で語 $v$ に割り当てられた回数（現在の単語を除く）です。$\alpha_k, \beta_v$ がそのまま分母分子に加算されている形に注目してください——これは前節で見たディリクレ事後分布の更新ルール「$\bm{\alpha} + \bm{n}$」がそのまま現れています。
LDAの数理的な土台がディリクレ分布であることが見えたところで、いよいよPythonで手を動かしていきます。
Python実装 — 可視化とサンプリング
scipyを使った確率密度の可視化
まず、3次元（$K=3$）のディリクレ分布をシンプレックス（正三角形）上に等高線図で可視化します。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from scipy.stats import dirichlet
from scipy.special import gammaln
# シンプレックス（正三角形）を等値線図用に三角分割
corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3) / 2]])
AREA = 0.5 * 1.0 * (np.sqrt(3) / 2)
triangle = mtri.Triangulation(corners[:, 0], corners[:, 1])
refiner = mtri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=6)
# (x, y) 座標から (pi_1, pi_2, pi_3) の重心座標に変換
def xy_to_barycentric(xy, tol=1e-4):
pairs = [corners[np.roll(range(3), -i)[1:]] for i in range(3)]
coords = np.array([
0.5 * np.linalg.norm(np.cross(*(pair - xy))) / AREA
for pair in pairs
])
return np.clip(coords, tol, 1.0 - tol)
# alpha の値ごとにディリクレ密度を計算
def dirichlet_pdf_on_simplex(alpha):
pvals = []
for xy in zip(trimesh.x, trimesh.y):
bary = xy_to_barycentric(np.array(xy))
bary = bary / bary.sum() # 数値誤差を吸収
pvals.append(dirichlet(alpha).pdf(bary))
return np.array(pvals)
alphas = [
(1.0, 1.0, 1.0),
(3.0, 3.0, 3.0),
(10.0, 10.0, 10.0),
(0.3, 0.3, 0.3),
(5.0, 1.0, 1.0),
(5.0, 5.0, 1.0),
]
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
for ax, alpha in zip(axes.ravel(), alphas):
pvals = dirichlet_pdf_on_simplex(alpha)
ax.tricontourf(trimesh, pvals, 60, cmap='viridis')
ax.triplot(triangle, color='white', linewidth=0.5)
ax.set_title(f"alpha={alpha}", fontsize=10)
ax.set_aspect('equal')
ax.axis('off')
plt.tight_layout()
plt.show()
このプロットからは、ディリクレ分布の形がパラメータでどう変わるかが一望できます。第一に、$\bm{\alpha} = (1,1,1)$ のときは完全に一様で、シンプレックス上のどの点も等確率になります。これは「事前知識ゼロ」のベースラインです。第二に、$\bm{\alpha}$ をすべて等しく大きくしていく（3→10）と分布は中心に向かって鋭く尖り、「3つの比率は均等」という強い信念を表します。第三に、$\bm{\alpha} = (0.3, 0.3, 0.3)$ のように1未満にすると、密度はシンプレックスの角（頂点）に集中します——これはLDAでスパース性を出すのに使われる重要な性質です。最後に、$\bm{\alpha} = (5, 1, 1)$ のように偏らせると、分布の山が「1番目の頂点」側に寄り、特定のカテゴリに偏った比率が生成されやすくなります。
サンプリングと統計量の確認
理論で導いた期待値・分散が、サンプリング結果と一致するか確認します。
import numpy as np
from scipy.stats import dirichlet
# 例: alpha = (2, 5, 3)
alpha = np.array([2.0, 5.0, 3.0])
alpha_0 = alpha.sum()
# 理論値
mean_theory = alpha / alpha_0
var_theory = alpha * (alpha_0 - alpha) / (alpha_0 ** 2 * (alpha_0 + 1))
print("理論平均:", mean_theory)
print("理論分散:", var_theory)
# サンプリング（n=100000）
samples = dirichlet(alpha).rvs(size=100000, random_state=0)
print("経験平均:", samples.mean(axis=0))
print("経験分散:", samples.var(axis=0))
# 共分散行列の比較
cov_theory = -np.outer(alpha, alpha) / (alpha_0 ** 2 * (alpha_0 + 1))
np.fill_diagonal(cov_theory, var_theory)
print("理論共分散:\n", cov_theory)
print("経験共分散:\n", np.cov(samples.T))
出力を見ると、理論値と経験値が小数点以下3〜4桁まで一致することが確認できます。サンプル数を $10^5$ に増やしているので、モンテカルロ誤差は $O(1/\sqrt{N}) \sim 0.003$ 程度に収まります。共分散行列の非対角成分が全て負になっている点も、シンプレックス制約から導いた性質と整合します。
ベイズ更新の可視化
事前分布から始めて、観測データを増やすにつれて事後分布がどう変わるかをアニメーション風に可視化します。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from scipy.stats import dirichlet
# 「真の」サイコロのパラメータ（3面ダイス）
pi_true = np.array([0.6, 0.3, 0.1])
# シンプレックスのメッシュ（前のセルと同じ三角分割を使う想定）
corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3) / 2]])
AREA = 0.5 * 1.0 * (np.sqrt(3) / 2)
triangle = mtri.Triangulation(corners[:, 0], corners[:, 1])
refiner = mtri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=6)
def xy_to_barycentric(xy, tol=1e-4):
pairs = [corners[np.roll(range(3), -i)[1:]] for i in range(3)]
coords = np.array([
0.5 * np.linalg.norm(np.cross(*(pair - xy))) / AREA
for pair in pairs
])
return np.clip(coords, tol, 1.0 - tol)
def plot_dirichlet(ax, alpha, title):
pvals = []
for xy in zip(trimesh.x, trimesh.y):
bary = xy_to_barycentric(np.array(xy))
bary = bary / bary.sum()
pvals.append(dirichlet(alpha).pdf(bary))
ax.tricontourf(trimesh, np.array(pvals), 50, cmap='viridis')
ax.triplot(triangle, color='white', linewidth=0.5)
ax.set_title(title, fontsize=10)
ax.set_aspect('equal')
ax.axis('off')
# 事前: 一様
alpha_prior = np.array([1.0, 1.0, 1.0])
rng = np.random.default_rng(42)
fig, axes = plt.subplots(1, 4, figsize=(14, 4))
N_list = [0, 5, 50, 500]
for ax, N in zip(axes, N_list):
if N == 0:
counts = np.zeros(3)
else:
samples = rng.multinomial(1, pi_true, size=N)
counts = samples.sum(axis=0)
alpha_post = alpha_prior + counts
title = f"N={N}, 観測={counts.astype(int).tolist()}"
plot_dirichlet(ax, alpha_post, title)
plt.suptitle("ディリクレ事後分布の更新（真の比率: 0.6, 0.3, 0.1）", y=1.02)
plt.tight_layout()
plt.show()
この4枚のサブプロットを左から順に見ると、ベイズ更新の動きがはっきり読み取れます。$N=0$ では一様分布で「分かりません」状態、$N=5$ で少し情報が入って分布の山ができはじめ、$N=50$ で真の比率 $(0.6, 0.3, 0.1)$ の周りに山が集中、$N=500$ では山が非常に鋭く尖って真値をピンポイントで捉えます。データが増えるほど分散が縮小していくのは、前節の分散式 $\alpha_k(\alpha_0-\alpha_k)/[\alpha_0^2(\alpha_0+1)]$ で $\alpha_0$ が大きくなる効果として直接説明できます。
ベータ分布との一致確認
$K=2$ のときディリクレ分布がベータ分布に一致することを、密度関数の値を直接比較して確認します。
import numpy as np
from scipy.stats import dirichlet, beta
alpha = np.array([2.5, 4.0])
xs = np.linspace(0.01, 0.99, 100)
# ディリクレで評価
dir_vals = np.array([dirichlet(alpha).pdf([x, 1 - x]) for x in xs])
# ベータで評価
beta_vals = beta.pdf(xs, alpha[0], alpha[1])
print("最大差:", np.max(np.abs(dir_vals - beta_vals)))
# 出力例: 最大差: 1.7763568394002505e-15
差は機械イプシロン（$10^{-15}$ 程度）の範囲に収まり、数値的にも $K=2$ のディリクレ分布がベータ分布と完全に一致することが確認できます。ベータ分布が「ディリクレ分布の特殊例」として包含されている、というのは単なる比喩ではなく、定義式と実装の両レベルで一致する厳密な事実です。
簡易LDAの実装（崩壊型ギブスサンプリング）
最後に、LDAの本格的な動作を体験するため、崩壊型ギブスサンプリングによる簡易LDAをスクラッチ実装します。データには「3つのトピックを混合した人工コーパス」を作って与え、LDAが本当にトピックを発見できるかを確かめます。
import numpy as np
from collections import defaultdict
rng = np.random.default_rng(0)
# 真のトピック (V=12 語彙)
V = 12
K_true = 3
phi_true = np.zeros((K_true, V))
# トピック0: 単語 0-3 に偏る
phi_true[0, 0:4] = 0.25
# トピック1: 単語 4-7 に偏る
phi_true[1, 4:8] = 0.25
# トピック2: 単語 8-11 に偏る
phi_true[2, 8:12] = 0.25
# 100 文書、各 30 語
D = 100
N_d = 30
theta_true = rng.dirichlet(np.array([0.5, 0.5, 0.5]), size=D)
docs = []
for d in range(D):
words = []
for n in range(N_d):
z = rng.choice(K_true, p=theta_true[d])
w = rng.choice(V, p=phi_true[z])
words.append(w)
docs.append(words)
# 崩壊型ギブスサンプラー
K = 3
alpha = 0.5
beta_ = 0.1
n_dk = np.zeros((D, K)) + alpha # 文書 d, トピック k のカウント (+事前)
n_kv = np.zeros((K, V)) + beta_ # トピック k, 語 v のカウント (+事前)
n_k = n_kv.sum(axis=1)
# 初期トピック割当てをランダム
z_assign = []
for d, words in enumerate(docs):
z_d = []
for w in words:
k = rng.integers(K)
z_d.append(k)
n_dk[d, k] += 1
n_kv[k, w] += 1
n_k[k] += 1
z_assign.append(z_d)
# ギブスサンプリング
n_iter = 200
for it in range(n_iter):
for d, words in enumerate(docs):
for n, w in enumerate(words):
k_old = z_assign[d][n]
# カウントから外す
n_dk[d, k_old] -= 1
n_kv[k_old, w] -= 1
n_k[k_old] -= 1
# 条件付き確率を計算
p_k = (n_dk[d] * n_kv[:, w] / n_k)
p_k /= p_k.sum()
# 新しいトピックをサンプル
k_new = rng.choice(K, p=p_k)
z_assign[d][n] = k_new
n_dk[d, k_new] += 1
n_kv[k_new, w] += 1
n_k[k_new] += 1
# 学習されたトピック–単語分布を表示
phi_learned = n_kv / n_k[:, None]
print("学習されたトピック–単語分布:")
for k in range(K):
top_words = np.argsort(phi_learned[k])[::-1][:5]
print(f" Topic {k}: 上位5語 = {top_words.tolist()}")
このコードを実行すると、3つの学習トピックがそれぞれ「単語0–3」「単語4–7」「単語8–11」の3グループに対応して上位5語を返すことが確認できます（トピック番号は順序が入れ替わる可能性があります）。これは、真のデータ生成過程をLDAが正しく逆推定できたことを意味します。崩壊型ギブスの中で各単語のトピック割当てを更新する1行——p_k = (n_dk[d] * n_kv[:, w] / n_k)——には、まさに前節で導いたディリクレ–カテゴリカル共役性の更新ルールが現れていることに注目してください。
LDA結果の可視化
学習されたトピックを単語の分布として可視化し、真のトピックと並べて比較します。
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, K, figsize=(12, 6))
for k in range(K):
axes[0, k].bar(np.arange(V), phi_true[k], color='steelblue')
axes[0, k].set_title(f"真のトピック {k}")
axes[0, k].set_ylim(0, 0.3)
axes[0, k].set_xlabel("単語ID")
axes[1, k].bar(np.arange(V), phi_learned[k], color='coral')
axes[1, k].set_title(f"学習トピック {k}")
axes[1, k].set_ylim(0, 0.3)
axes[1, k].set_xlabel("単語ID")
axes[0, 0].set_ylabel("真のphi")
axes[1, 0].set_ylabel("学習phi")
plt.tight_layout()
plt.show()
棒グラフを並べると、真のトピックが「単語0–3に集中」「単語4–7に集中」「単語8–11に集中」のクラスタを成しているのに対し、学習されたトピックも同じクラスタ構造（順序は異なる可能性あり）を再現していることが視覚的に確認できます。完全一致ではなく、各クラスタ外の単語にも小さな確率が残っているのは、$\beta = 0.1$ のスムージングと有限サンプル効果によるものです。LDAが「単語の共起パターンからトピックを発見する」という抽象的なタスクを、ディリクレ事前を使ったベイズ推論で確実に達成できる様子が、この実験で具体的に体感できます。
実装で確認できた事項を踏まえて、最後にディリクレ分布の更なる拡張に進みます。
拡張: ディリクレ過程と非パラメトリックベイズ
トピック数を事前に決めなくてよいモデル
通常のLDAでは「トピック数 $K$」を分析者が事前に指定する必要があります。3トピックなのか100トピックなのか、データを見るまで分からないのが現実なのに——この問題に答えるのがディリクレ過程（Dirichlet Process、DP）です。
ディリクレ過程は「ディリクレ分布の可算無限次元への拡張」と表現されることが多いです。形式的には、ベース測度 $G_0$ と集中度パラメータ $\alpha$ を持つDPからのサンプル $G$ は、可算個の点質量の重ね合わせ
$$ G = \sum_{k=1}^{\infty} \pi_k\, \delta_{\theta_k} $$
として表せます。ここで $\pi_k$ は stick-breaking構成 で生成される無限次元の比率ベクトルで、$\sum_{k=1}^\infty \pi_k = 1$ を満たします。具体的には
$$ V_k \sim \mathrm{Beta}(1, \alpha), \quad \pi_k = V_k \prod_{j=1}^{k-1}(1 – V_j) $$
「長さ1の棒を毎回 $V_k$ の割合で折り取っていく」操作で、無限個の比率を作っていくイメージです。$\alpha$ が小さければ最初の数本に重みが集中し、$\alpha$ が大きければなだらかに減少していく——「事実上いくつのトピックが使われるか」を $\alpha$ が間接的に制御するわけです。
中華料理店過程
ディリクレ過程のもう一つの便利な表現が中華料理店過程（Chinese Restaurant Process、CRP）です。客が次々と中華料理店に入店し、テーブルを選ぶ確率モデル：
- 新しい客はすでに $n_k$ 人座っているテーブル $k$ に確率 $n_k/(n-1+\alpha)$ で座る
- 新しいテーブルに確率 $\alpha/(n-1+\alpha)$ で座る
これがそのまま「クラスタ数が事前に分からないクラスタリング」のサンプリングルールになります。LDAをDPに拡張したHDP-LDA（階層ディリクレ過程LDA）では、トピック数がデータから自動的に決まるため、ハイパーパラメータ調整の負担が大幅に減ります。
簡易stick-breakingの実装
stick-breaking構成をPythonで実装し、ディリクレ過程からのサンプルを可視化します。
import numpy as np
import matplotlib.pyplot as plt
def stick_breaking(alpha, K_max=50, seed=0):
rng = np.random.default_rng(seed)
V = rng.beta(1.0, alpha, size=K_max)
pi = np.zeros(K_max)
remaining = 1.0
for k in range(K_max):
pi[k] = V[k] * remaining
remaining *= (1.0 - V[k])
return pi
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for ax, alpha in zip(axes, [0.5, 2.0, 10.0]):
pi = stick_breaking(alpha)
ax.bar(np.arange(50), pi, color='teal')
ax.set_title(f"alpha={alpha}")
ax.set_xlabel("k")
ax.set_ylabel("pi_k")
ax.set_ylim(0, max(pi) * 1.1)
plt.tight_layout()
plt.show()
3つのサブプロットを比較すると、$\alpha$ の値で挙動がはっきり変わることが分かります。$\alpha = 0.5$ では最初の1〜2成分にほぼ全質量が集中し、「少数のクラスタが支配的」な状態。$\alpha = 2.0$ では5〜10成分まで意味のある重みが分散し、$\alpha = 10.0$ ではさらに多くの成分まで重みが滑らかに広がります。この性質により、データの複雑度に応じて「自然な」クラスタ数が選ばれるのです。
上段の図は「棒（長さ1）を順番に折り取って $\pi_1, \pi_2, \dots$ を得る」という stick-breaking の幾何学的な意味を直接示しています。下段は同じ手続きを $\alpha$ を変えて繰り返した結果で、$\alpha$ が小さいほど最初の数本に質量が集中するという定性的な性質を、棒グラフで定量的に確認できます。
DPの応用例
ディリクレ過程の応用は広く、混合モデルのクラスタ数自動決定（DPMM: Dirichlet Process Mixture Model）、トピックモデル（HDP-LDA）、ベイズ非パラメトリック回帰、自然言語処理での文法学習など、「サイズが事前に未知のモデル」全般で使われます。実装の詳細は別記事に譲りますが、ディリクレ分布の基本性質が分かれば、その無限次元拡張も同じロジックで理解できることを覚えておいてください。
まとめ
本記事では、ディリクレ分布の数学的定義から応用までを一気通貫で解説しました。要点を整理します。
- 直感：ディリクレ分布は、合計1の比率ベクトル $\bm{\pi}$（確率シンプレックス）上に定義される分布。「比率の確率分布」を扱う場面で常に登場する
- 定義：$\mathrm{Dir}(\bm{\pi}|\bm{\alpha}) = \frac{1}{B(\bm{\alpha})}\prod_k \pi_k^{\alpha_k-1}$。正規化定数は多次元ベータ関数 $B(\bm{\alpha}) = \prod_k \Gamma(\alpha_k)/\Gamma(\alpha_0)$
- 統計量：期待値 $\mathbb{E}[\pi_k] = \alpha_k/\alpha_0$、分散 $\mathrm{Var}[\pi_k] = \alpha_k(\alpha_0-\alpha_k)/[\alpha_0^2(\alpha_0+1)]$。集中度 $\alpha_0$ が大きいほど分布は鋭く尖る
- 共役性：多項分布の共役事前分布。事後はパラメータに観測カウントを足すだけ $\mathrm{Dir}(\bm{\alpha} + \bm{n})$。$\alpha_k$ は疑似カウントとして解釈できる
- 応用：LDAトピックモデルでは、文書のトピック比率とトピックの単語比率の両方の事前分布として使われる。$\bm{\alpha} < 1$ でスパース性を制御
- 拡張：ディリクレ過程は無限次元への一般化で、クラスタ数を事前に決めなくてよい非パラメトリックベイズの基盤
ディリクレ分布は、それ単体としてだけでなく、ベイズ推定、トピックモデル、混合モデル、ベイズ非パラメトリックといった広範な手法の中で繰り返し登場する、いわば「ベイズ機械学習の屋台骨」のような分布です。一見すると数式が複雑に見えますが、本記事で見たようにベータ分布の素直な拡張として理解すれば、各要素が自然に位置付きます。
次のステップとして、以下の記事も参考にしてください。