サイコロを10回振ったら、1の目が4回も出たとします。「このサイコロは1が出やすいように歪んでいるのでは？」と疑いたくなりますが、たった10回ではたまたまかもしれません。では何回振れば「歪んでいる」と自信を持って言えるのでしょうか。そして、まだ20回しか振っていない段階で「次に1が出る確率」を見積もるには、どう計算すればよいのでしょうか。
この「カテゴリが複数ある事象（サイコロの6つの目、文章中の単語、商品のクリック有無など）について、各カテゴリの出現確率を観測データから推定し、しかも観測が少ないうちは不確かさも込みで扱いたい」という問題に、ベイズの枠組みで美しく答えるのがディリクレ分布とカテゴリカル分布の組み合わせです。
この2つの分布は、実用上きわめて広い応用を持っています。
- 自然言語処理のトピックモデル（LDA）: 文書がどのトピックに属するか、トピックがどの単語を生成するかを、すべてディリクレ分布で表現します。
- A/Bテストやレコメンド: 複数の選択肢のクリック率を、観測の少なさによる不確かさを保ったまま推定し、トンプソンサンプリングなどの意思決定に使います。
この記事では、ディリクレ分布をベータ分布の多変量への一般化として導入し、カテゴリカル分布（多項分布）との共役性を一行ずつ丁寧に導出します。そのうえで、確率パラメータ $\bm{\pi}$ を積分で消去した予測分布（周辺化）を導き、最後に偏ったサイコロのベイズ推定をPythonで実装して、観測が増えるにつれて事後分布が真の値に集中していく様子を可視化します。
本記事の内容
- カテゴリカル分布とディリクレ分布の直感的な理解と数学的定義
- 共役性（事後分布が再びディリクレ分布になること）の省略なしの導出
- $\bm{\pi}$ を周辺化した事後予測分布の導出と、MAP・事後平均の比較
- Pythonでの単体上の可視化と、偏ったサイコロのベイズ推定
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ベータ分布をわかりやすく解説 — ディリクレ分布はベータ分布の多変量版です。まずこちらで2カテゴリの場合を押さえておくと理解が一気に楽になります。
- ディリクレ分布をわかりやすく解説 — 分布そのものの形状や性質に焦点を当てた記事です。
- ベイズの定理の基礎 — 事前分布・尤度・事後分布の関係をおさらいできます。
- カテゴリカル分布のベイズ推定 — 本記事と密接に関連する応用記事です。
カテゴリカル分布とは
まず、推定したい「データを生み出す仕組み」の側を整理します。サイコロを1回振ると、$\{1, 2, 3, 4, 5, 6\}$ のいずれか1つの目が出ます。各目が出る確率を $\pi_1, \pi_2, \dots, \pi_6$ とすると、これらは合計して1になります。このように「$K$ 個のカテゴリのうちちょうど1つが選ばれ、カテゴリ $k$ が選ばれる確率が $\pi_k$」という最も基本的な離散分布がカテゴリカル分布です。
コインのように結果が2通り（表・裏）ならベルヌーイ分布ですが、結果が3通り以上になった一般化がカテゴリカル分布だと考えてください。サイコロ、トランプのスート、アンケートの選択肢、文章中の単語など、「いくつかの選択肢から1つ」という状況はすべてカテゴリカル分布で表せます。
カテゴリカル分布を数式で書くために、結果をワンホット表現で表します。カテゴリ $k$ が出たことを、$k$ 番目だけが1で残りが0のベクトル $\bm{x} = (x_1, \dots, x_K)$、ただし $\sum_k x_k = 1$、で表します。このとき1回の試行の確率は
$$ \begin{equation} p(\bm{x} \mid \bm{\pi}) = \prod_{k=1}^{K} \pi_k^{x_k} \end{equation} $$
と書けます。$x_k$ は0か1なので、実際に出たカテゴリ $k$ では $\pi_k^1 = \pi_k$、それ以外では $\pi_j^0 = 1$ となり、積を取ると「出たカテゴリの確率 $\pi_k$」だけが残ります。指数 $x_k$ を使うことで、6本の場合分けを1本の式にまとめている点がポイントです。
ここで確率パラメータ $\bm{\pi} = (\pi_1, \dots, \pi_K)$ は
$$ \pi_k \geq 0, \qquad \sum_{k=1}^{K} \pi_k = 1 $$
を満たします。この制約を満たす $\bm{\pi}$ の集合は、幾何学的には $(K-1)$ 次元の単体（シンプレックス）と呼ばれる図形になります。$K=3$ なら、3頂点を $(1,0,0), (0,1,0), (0,0,1)$ とする正三角形の内部です。
多項分布への拡張
サイコロを1回ではなく $N$ 回振り、カテゴリ $k$ が出た回数を $n_k$ とすると（$\sum_k n_k = N$）、回数の組 $\bm{n} = (n_1, \dots, n_K)$ が従う分布が多項分布です。
$$ \begin{equation} p(\bm{n} \mid \bm{\pi}, N) = \frac{N!}{n_1! \, n_2! \cdots n_K!} \prod_{k=1}^{K} \pi_k^{n_k} \end{equation} $$
前の係数は「どの順番で各目が出たか」の場合の数（多項係数）です。本記事の推定では、$\bm{\pi}$ にかかる部分 $\prod_k \pi_k^{n_k}$ だけが本質的に効いてきます。$N$ 回の独立なカテゴリカル試行の同時確率を取ると、ちょうどこの $\prod_k \pi_k^{n_k}$ が現れるからです。
さて、ここで自然な疑問が生まれます。観測されるのは目の出方 $\bm{n}$ だけで、肝心の確率 $\bm{\pi}$ は隠れています。この $\bm{\pi}$ 自体を「確率変数」とみなして、その分布を考えたい。では、単体の上で定義される確率分布をどう作ればよいのでしょうか。その答えがディリクレ分布です。
ディリクレ分布とは
確率の組 $\bm{\pi}$ は単体の上に住んでいます。「3つの確率 $(\pi_1, \pi_2, \pi_3)$ がどのあたりにありそうか」という信念を、正三角形の上の山の高さ（密度）として表したい。たとえば「どの目も同じくらい出そう」なら三角形の中心が高い分布、「1の目が出やすそう」なら頂点1に寄った分布、という具合です。この「確率の上の確率分布」を与えるのがディリクレ分布です。
コインの表が出る確率 $p$（1次元）に対する信念をベータ分布で表したのと同じことを、多次元の単体の上で行うものだ、と考えてください。実際、後で見るようにディリクレ分布はベータ分布をそのまま多変量に拡張したものです。
ディリクレ分布の密度関数は、パラメータ $\bm{\alpha} = (\alpha_1, \dots, \alpha_K)$（すべて $\alpha_k > 0$）を持ち、次の形をしています。
$$ \begin{equation} \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha}) = \frac{1}{B(\bm{\alpha})} \prod_{k=1}^{K} \pi_k^{\alpha_k – 1} \end{equation} $$
まず注目すべきは右側の $\prod_k \pi_k^{\alpha_k – 1}$ です。これはカテゴリカル尤度 $\prod_k \pi_k^{x_k}$ とそっくりな形をしています。この「形が同じ」ことこそが、あとで共役性を生む鍵になります。指数が $\alpha_k – 1$ と「マイナス1」されているのは、ベータ分布の指数が $\alpha-1, \beta-1$ だったのと同じ理由で、規格化の都合によるものです。
前の係数 $1/B(\bm{\alpha})$ は、単体上で積分して1になるようにする規格化定数です。$B(\bm{\alpha})$ は多変量ベータ関数と呼ばれ、ガンマ関数 $\Gamma$ を使って
$$ \begin{equation} B(\bm{\alpha}) = \frac{\prod_{k=1}^{K} \Gamma(\alpha_k)}{\Gamma\!\left(\sum_{k=1}^{K} \alpha_k\right)} \end{equation} $$
と書けます。これは1次元のベータ関数 $B(\alpha, \beta) = \Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$ の自然な多変量版です。
パラメータ $\bm{\alpha}$ の意味
パラメータ $\alpha_k$ は、直感的には「カテゴリ $k$ を事前に何回くらい見たことにするか」という疑似カウントとして解釈できます。
- すべての $\alpha_k = 1$ のとき、$\prod_k \pi_k^{0} = 1$ となり、密度は単体上で一定（一様分布）になります。「どの確率の組も等しくありうる」という無情報な事前知識を表します。
- $\alpha_k$ を大きくすると、その分布はカテゴリ $k$ の確率 $\pi_k$ が大きい方へ偏ります。
- すべての $\alpha_k$ を一律に大きくする（例：すべて100）と、分布は単体の中心に鋭く集中します。「どの確率もだいたい $1/K$ のはずだ」という強い信念です。
- 逆にすべての $\alpha_k < 1$ にすると、密度は単体の頂点や辺に集中し、「どれか1つのカテゴリだけが支配的」という疎な状況を表します。
総和 $\alpha_0 = \sum_k \alpha_k$ は集中度（concentration）と呼ばれ、大きいほど分布が平均の周りに鋭くなります。ベイズ推定では、この $\alpha_0$ が「事前にどれだけのデータを見たつもりか」という事前知識の強さに対応します。
ベータ分布との関係
$K=2$ の場合を考えると、$\pi_2 = 1 – \pi_1$ なので自由なパラメータは $\pi_1$ だけです。ディリクレ分布の密度は
$$ \mathrm{Dir}(\pi_1, \pi_2 \mid \alpha_1, \alpha_2) = \frac{\Gamma(\alpha_1 + \alpha_2)}{\Gamma(\alpha_1)\Gamma(\alpha_2)} \, \pi_1^{\alpha_1 – 1} (1 – \pi_1)^{\alpha_2 – 1} $$
となり、これはまさにベータ分布 $\mathrm{Beta}(\pi_1 \mid \alpha_1, \alpha_2)$ そのものです。つまりディリクレ分布は、ベータ分布をカテゴリ数 $K$ 個に拡張したものだと正確に言えます。ベータ分布をすでに知っていれば、頭の中で「2 → $K$」と置き換えるだけで理解できます。
ここまでで、データを生むカテゴリカル分布と、確率の信念を表すディリクレ分布が揃いました。次は、この2つを組み合わせたときに何が起きるか——ベイズ更新で事後分布がどうなるかを導出します。
共役性の導出：事後分布が再びディリクレになる
ベイズ推論の心臓部は、ベイズの定理
$$ \begin{equation} p(\bm{\pi} \mid \mathcal{D}) = \frac{p(\mathcal{D} \mid \bm{\pi}) \, p(\bm{\pi})}{p(\mathcal{D})} \end{equation} $$
です。事前分布 $p(\bm{\pi})$ に観測データ $\mathcal{D}$ から得た尤度 $p(\mathcal{D} \mid \bm{\pi})$ を掛け、規格化すると事後分布 $p(\bm{\pi} \mid \mathcal{D})$ が得られます。ここで共役事前分布とは、「事前分布と事後分布が同じ分布族になる」ような事前分布のことです。共役性が成り立つと、ベイズ更新が「分布の形を変えずパラメータを更新するだけ」の簡単な操作になります。
ここで示したいのは次のことです。ディリクレ分布を事前分布、カテゴリカル（多項）分布を尤度に選ぶと、事後分布も再びディリクレ分布になる。 これを一行ずつ導きましょう。
データの尤度
$N$ 個の独立な観測 $\mathcal{D} = \{\bm{x}^{(1)}, \dots, \bm{x}^{(N)}\}$（各 $\bm{x}^{(i)}$ はワンホット）があり、カテゴリ $k$ が出た回数を $n_k = \sum_{i} x_k^{(i)}$ とします。独立性から尤度は各試行の確率の積です。
$$ \begin{align} p(\mathcal{D} \mid \bm{\pi}) &= \prod_{i=1}^{N} \prod_{k=1}^{K} \pi_k^{x_k^{(i)}} \end{align} $$
同じ底 $\pi_k$ の積は指数の和にまとめられます。$i$ についての和を先に取ると、指数は $\sum_i x_k^{(i)} = n_k$ になります。
$$ \begin{align} p(\mathcal{D} \mid \bm{\pi}) &= \prod_{k=1}^{K} \pi_k^{\sum_{i=1}^{N} x_k^{(i)}} = \prod_{k=1}^{K} \pi_k^{n_k} \end{align} $$
つまり尤度は、各カテゴリのカウント $n_k$ を指数に持つ積 $\prod_k \pi_k^{n_k}$ にきれいにまとまりました。多項分布の係数 $N!/\prod_k n_k!$ は $\bm{\pi}$ に依存しない定数なので、$\bm{\pi}$ の関数として見ると無視できます。
事前分布と尤度の積
事前分布をディリクレ分布 $\mathrm{Dir}(\bm{\pi} \mid \bm{\alpha})$ に取ります。事後分布はベイズの定理より、尤度と事前分布の積に比例します（分母 $p(\mathcal{D})$ は $\bm{\pi}$ に依存しない定数なので、いったん比例関係 $\propto$ で進めます）。
$$ \begin{align} p(\bm{\pi} \mid \mathcal{D}) &\propto p(\mathcal{D} \mid \bm{\pi}) \, \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha}) \end{align} $$
それぞれを代入します。尤度に $\prod_k \pi_k^{n_k}$、ディリクレ事前に $\frac{1}{B(\bm{\alpha})}\prod_k \pi_k^{\alpha_k – 1}$ を入れると
$$ \begin{align} p(\bm{\pi} \mid \mathcal{D}) &\propto \left(\prod_{k=1}^{K} \pi_k^{n_k}\right) \cdot \frac{1}{B(\bm{\alpha})}\left(\prod_{k=1}^{K} \pi_k^{\alpha_k – 1}\right) \end{align} $$
規格化定数 $1/B(\bm{\alpha})$ は $\bm{\pi}$ によらない定数なので比例関係に吸収できます。残った2つの積を、同じ底 $\pi_k$ ごとにまとめると、指数が足し算になります。
$$ \begin{align} p(\bm{\pi} \mid \mathcal{D}) &\propto \prod_{k=1}^{K} \pi_k^{n_k} \cdot \pi_k^{\alpha_k – 1} = \prod_{k=1}^{K} \pi_k^{\alpha_k + n_k – 1} \end{align} $$
ここで指数を見比べてみましょう。$\alpha_k + n_k – 1$ という形は、ディリクレ分布の密度 $\prod_k \pi_k^{\alpha_k’ – 1}$ の形と、$\alpha_k’ = \alpha_k + n_k$ とおけば完全に一致します。すなわち事後分布は、パラメータを $\bm{\alpha} + \bm{n}$ に置き換えたディリクレ分布にほかなりません。
$$ \begin{equation} \boxed{\; p(\bm{\pi} \mid \mathcal{D}) = \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha} + \bm{n}) \;} \end{equation} $$
比例関係から等号にできるのは、両辺ともに単体上で積分して1になる確率分布であり、同じ関数形（同じ指数）を持つなら規格化定数まで一致するからです。事後分布の規格化定数は自動的に $1/B(\bm{\alpha} + \bm{n})$ になります。
この結果は驚くほどシンプルです。ベイズ更新は「事前のパラメータ $\alpha_k$ に観測カウント $n_k$ を足すだけ」になります。$\alpha_k$ が「疑似カウント」と呼ばれた理由がここで明確になります。事前知識をあたかも「カテゴリ $k$ をすでに $\alpha_k$ 回見たデータ」のように扱い、そこに実データのカウント $n_k$ を足し込むのです。
ここまでで $\bm{\pi}$ の事後分布が手に入りました。しかし実用では「$\bm{\pi}$ の分布」そのものより、「次にどのカテゴリが出るか」を知りたい場面が多いです。次は $\bm{\pi}$ を積分で消し去った予測分布を導きます。
事後予測分布：$\bm{\pi}$ を周辺化する
サイコロを20回振って結果を観測したあと、「21回目に1の目が出る確率は？」と問うとします。素朴には事後平均 $\hat{\pi}_1$ を答えればよさそうですが、ベイズの立場ではもっと丁寧に、$\bm{\pi}$ の不確かさをすべて積分で平均化した事後予測分布を計算します。
$\bm{\pi}$ が不確かなので、考えうるすべての $\bm{\pi}$ について「その $\bm{\pi}$ のもとで次にカテゴリ $k$ が出る確率」を、事後分布で重み付けして足し合わせます。この「不要なパラメータを積分で消す」操作を周辺化（marginalization）と呼びます。新しい観測を $\bm{x}^{\text{new}}$（カテゴリ $k$ がワンホットで1）とすると
$$ \begin{align} p(x_k^{\text{new}} = 1 \mid \mathcal{D}) &= \int p(x_k^{\text{new}} = 1 \mid \bm{\pi}) \, p(\bm{\pi} \mid \mathcal{D}) \, d\bm{\pi} \end{align} $$
これは「カテゴリカル尤度を、$\bm{\pi}$ の事後分布で期待値を取ったもの」です。$p(x_k^{\text{new}}=1 \mid \bm{\pi}) = \pi_k$ なので、結局これは事後分布のもとでの $\pi_k$ の期待値、つまり事後平均そのものです。
$$ \begin{align} p(x_k^{\text{new}} = 1 \mid \mathcal{D}) = \int \pi_k \, \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha} + \bm{n}) \, d\bm{\pi} = \mathbb{E}[\pi_k \mid \mathcal{D}] \end{align} $$
ディリクレ分布の期待値
そこでディリクレ分布の期待値を計算しておきます。パラメータ $\bm{\alpha}’$（ここでは $\bm{\alpha}’ = \bm{\alpha} + \bm{n}$）を持つディリクレ分布のもとで、$\pi_k$ の期待値は
$$ \begin{equation} \mathbb{E}[\pi_k] = \frac{\alpha_k’}{\sum_{j} \alpha_j’} = \frac{\alpha_k’}{\alpha_0′} \end{equation} $$
となります（$\alpha_0′ = \sum_j \alpha_j’$）。これを導いておきましょう。期待値は密度に $\pi_k$ を掛けて積分したものです。
$$ \begin{align} \mathbb{E}[\pi_k] &= \int \pi_k \cdot \frac{1}{B(\bm{\alpha}’)} \prod_{j} \pi_j^{\alpha_j’ – 1} \, d\bm{\pi} \end{align} $$
被積分関数の中で $\pi_k \cdot \pi_k^{\alpha_k’ – 1} = \pi_k^{\alpha_k’}$ とまとめると、これは $k$ 番目のパラメータだけを $\alpha_k’ + 1$ に増やしたディリクレ分布の「規格化前の形」です。よって積分は新しい多変量ベータ関数 $B(\bm{\alpha}’ + \bm{e}_k)$ に等しくなります（$\bm{e}_k$ は $k$ 番目だけ1のベクトル）。
$$ \begin{align} \mathbb{E}[\pi_k] &= \frac{B(\bm{\alpha}’ + \bm{e}_k)}{B(\bm{\alpha}’)} \end{align} $$
ここで $B(\bm{\alpha}’) = \frac{\prod_j \Gamma(\alpha_j’)}{\Gamma(\alpha_0′)}$ を代入し、比を取ります。分子では $\alpha_k’$ が $\alpha_k’ + 1$ に、総和 $\alpha_0’$ が $\alpha_0′ + 1$ に変わるので
$$ \begin{align} \mathbb{E}[\pi_k] &= \frac{\Gamma(\alpha_k’ + 1)}{\Gamma(\alpha_k’)} \cdot \frac{\Gamma(\alpha_0′)}{\Gamma(\alpha_0′ + 1)} \end{align} $$
ガンマ関数の基本性質 $\Gamma(z+1) = z\,\Gamma(z)$ を使うと、$\Gamma(\alpha_k’+1)/\Gamma(\alpha_k’) = \alpha_k’$、$\Gamma(\alpha_0′)/\Gamma(\alpha_0’+1) = 1/\alpha_0’$ となります。代入すると
$$ \begin{align} \mathbb{E}[\pi_k] &= \alpha_k’ \cdot \frac{1}{\alpha_0′} = \frac{\alpha_k’}{\alpha_0′} \end{align} $$
が得られます。期待値の導出が完了しました。
予測分布の最終形
$\bm{\alpha}’ = \bm{\alpha} + \bm{n}$ を代入すると、求める事後予測分布は
$$ \begin{equation} \boxed{\; p(x_k^{\text{new}} = 1 \mid \mathcal{D}) = \frac{\alpha_k + n_k}{\sum_{j=1}^{K}(\alpha_j + n_j)} = \frac{\alpha_k + n_k}{\alpha_0 + N} \;} \end{equation} $$
となります（$\alpha_0 = \sum_j \alpha_j$、$N = \sum_j n_j$）。これがこの記事の中心的な結果の一つです。「次にカテゴリ $k$ が出る確率」は、$\alpha_k + n_k$ に比例します。
この式は、ラプラスの加法的平滑化（add-$\alpha$ smoothing）そのものです。単純に頻度 $n_k / N$ を使うと、一度も観測されなかったカテゴリ（$n_k = 0$）の確率がゼロになってしまいますが、$\alpha_k$ を足すことで「まだ見ていないだけかもしれない」という不確かさを反映し、確率がゼロにならないようにします。$N$ が大きくなると $\alpha_k$ の影響は薄れ、予測は経験頻度 $n_k/N$ に近づきます。
MAP推定と事後平均の違い
予測分布は事後平均 $\frac{\alpha_k + n_k}{\alpha_0 + N}$ を使いますが、点推定としてよく対比されるのがMAP推定（事後確率最大）です。MAPは事後分布の密度を最大化する $\bm{\pi}$ で、ラグランジュの未定乗数法で制約 $\sum_k \pi_k = 1$ のもとで $\prod_k \pi_k^{\alpha_k + n_k – 1}$ を最大化すると
$$ \begin{equation} \hat{\pi}_k^{\text{MAP}} = \frac{\alpha_k + n_k – 1}{\sum_{j}(\alpha_j + n_j – 1)} = \frac{\alpha_k + n_k – 1}{\alpha_0 + N – K} \end{equation} $$
となります。事後平均と比べると、分子・分母が「$-1$」「$-K$」だけずれています。この差は指数 $\alpha_k + n_k – 1$ の「$-1$」に由来し、分布の最頻値（モード）と平均値の違いを反映しています。観測が少ないうちは両者は大きく異なりますが、$N$ が大きくなるとどちらも経験頻度に収束します。事後平均は予測誤差（二乗損失）を最小化する点推定であり、不確かさの積分を通った値なので、予測には事後平均を使うのが自然です。
これで理論は完成です。共役更新（$\bm{\alpha} \to \bm{\alpha}+\bm{n}$）、予測分布（$\propto \alpha_k + n_k$）、MAPと事後平均の違いがそろいました。次はこれらをPythonで実装し、目で確かめます。
Pythonでの実装
単体上のディリクレ分布の可視化（K=3）
まず $K=3$ のディリクレ分布が単体（正三角形）の上でどんな形をしているかを見ます。パラメータ $\bm{\alpha}$ を変えると分布の形がどう変わるかを確認します。3次元の確率 $(\pi_1, \pi_2, \pi_3)$ を、正三角形の重心座標に変換して2次元平面に描きます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import dirichlet
# 単体（正三角形）の3頂点。重心座標で (1,0,0),(0,1,0),(0,0,1) に対応
corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3) / 2]])
def bary_to_cart(pi):
"""重心座標 pi=(p1,p2,p3) を三角形上の直交座標へ変換"""
return pi @ corners
# 三角形内部を細かくサンプリングするグリッドを作る
res = 200
xs, ys = np.meshgrid(np.linspace(0, 1, res), np.linspace(0, np.sqrt(3) / 2, res))
ここまでで、三角形上の各点を確率の組 $(\pi_1,\pi_2,\pi_3)$ に対応づける準備ができました。次に、各点でディリクレ密度を評価して色で塗ります。
def cart_to_bary(x, y):
"""直交座標 (x,y) を重心座標 (p1,p2,p3) に逆変換"""
T = np.array([[corners[0, 0] - corners[2, 0], corners[1, 0] - corners[2, 0]],
[corners[0, 1] - corners[2, 1], corners[1, 1] - corners[2, 1]]])
lam = np.linalg.solve(T, np.array([x - corners[2, 0], y - corners[2, 1]]))
return np.array([lam[0], lam[1], 1 - lam[0] - lam[1]])
def dirichlet_density_grid(alpha):
"""グリッド上でディリクレ密度を評価（三角形外は NaN）"""
Z = np.full(xs.shape, np.nan)
for i in range(res):
for j in range(res):
pi = cart_to_bary(xs[i, j], ys[i, j])
if np.all(pi > 1e-6): # 三角形の内部のみ
Z[i, j] = dirichlet.pdf(pi, alpha)
return Z
alphas = [np.array([1.0, 1.0, 1.0]), # 一様
np.array([5.0, 5.0, 5.0]), # 中心に集中
np.array([2.0, 5.0, 10.0])] # 頂点3寄り
titles = ["alpha=(1,1,1) uniform", "alpha=(5,5,5) concentrated", "alpha=(2,5,10) skewed"]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, a, t in zip(axes, alphas, titles):
Z = dirichlet_density_grid(a)
ax.contourf(xs, ys, Z, levels=30, cmap="viridis")
ax.plot(*np.append(corners, [corners[0]], axis=0).T, "w-", lw=1)
ax.set_title(t)
ax.axis("equal"); ax.axis("off")
plt.tight_layout()
plt.savefig("dirichlet_simplex.png", dpi=150, bbox_inches="tight")
plt.show()
3枚の図から、パラメータ $\bm{\alpha}$ の役割がはっきり読み取れます。左の $(1,1,1)$ では密度が三角形全体で一定（一様分布）で、「どの確率の組も等しくありうる」という無情報な事前分布を表します。中央の $(5,5,5)$ では密度が三角形の中心 $(1/3,1/3,1/3)$ に集中しており、「3カテゴリはどれもだいたい同じ確率のはず」という強い信念を表します。右の $(2,5,10)$ では密度のピークが頂点3（$\pi_3$ が大きい角）の方へ寄っており、$\alpha_k$ の大小がそのまま確率の偏りに対応することがわかります。
これで事前分布の形を視覚的に理解できました。次は、データを観測したときに事前分布がどう更新され、真の値へ集中していくかを見ます。
観測数の増加に伴う事後分布の集中
偏ったサイコロ（実は3面のサイコロとして $K=3$）を考え、真の確率を $\bm{\pi}^* = (0.2, 0.3, 0.5)$ とします。観測数 $N$ を増やしながら、事後分布 $\mathrm{Dir}(\bm{\alpha} + \bm{n})$ が真の値の周りに集中していく様子を可視化します。事前分布は一様 $\bm{\alpha} = (1,1,1)$ から始めます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import dirichlet
np.random.seed(0)
pi_true = np.array([0.2, 0.3, 0.5]) # 真のカテゴリ確率
alpha_prior = np.array([1.0, 1.0, 1.0]) # 一様事前分布
corners = np.array([[0.0, 0.0], [1.0, 0.0], [0.5, np.sqrt(3) / 2]])
res = 200
xs, ys = np.meshgrid(np.linspace(0, 1, res), np.linspace(0, np.sqrt(3) / 2, res))
def cart_to_bary(x, y):
T = np.array([[corners[0, 0] - corners[2, 0], corners[1, 0] - corners[2, 0]],
[corners[0, 1] - corners[2, 1], corners[1, 1] - corners[2, 1]]])
lam = np.linalg.solve(T, np.array([x - corners[2, 0], y - corners[2, 1]]))
return np.array([lam[0], lam[1], 1 - lam[0] - lam[1]])
def density_grid(alpha):
Z = np.full(xs.shape, np.nan)
for i in range(res):
for j in range(res):
pi = cart_to_bary(xs[i, j], ys[i, j])
if np.all(pi > 1e-6):
Z[i, j] = dirichlet.pdf(pi, alpha)
return Z
ここで、各観測数 $N$ についてサイコロを振ってカウント $\bm{n}$ を作り、事後パラメータ $\bm{\alpha}+\bm{n}$ で密度を描きます。
N_list = [0, 10, 50, 500]
fig, axes = plt.subplots(1, 4, figsize=(18, 4.5))
for ax, N in zip(axes, N_list):
# N 回サイコロを振ってカウント n を作る（多項分布からのサンプル）
counts = np.random.multinomial(N, pi_true) if N > 0 else np.zeros(3)
alpha_post = alpha_prior + counts # 共役更新： alpha + n
Z = density_grid(alpha_post)
ax.contourf(xs, ys, Z, levels=30, cmap="magma")
ax.plot(*np.append(corners, [corners[0]], axis=0).T, "w-", lw=1)
# 真の値の位置に赤丸
true_xy = pi_true @ corners
ax.plot(*true_xy, "ro", ms=8)
ax.set_title(f"N={N}, counts={counts.astype(int)}")
ax.axis("equal"); ax.axis("off")
plt.tight_layout()
plt.savefig("dirichlet_posterior_concentration.png", dpi=150, bbox_inches="tight")
plt.show()
4枚の図は、ベイズ更新の本質を一目で示しています。$N=0$（左）では事前分布の一様分布のままで、確率の組について何も分かっていません。$N=10$ では分布が真の値（赤丸）の方向にゆるく寄りますが、まだ広く広がっており「不確かさが大きい」ことを表しています。$N=50$、$N=500$ と観測が増えるにつれて分布は真の値の周りに鋭く集中していき、推定の確信度が高まる様子がはっきり見えます。重要なのは、データが少ないうちは「広がり（不確かさ）」を分布として保持し、データが増えるとともに自動的に確信を強めていく点です。これが点推定にはないベイズ推論の強みです。
集中していく様子が見えたので、次は点推定（事後平均とMAP）が観測数とともにどう真の値に近づくかを数値で比較します。
事後平均とMAPの比較
偏ったサイコロ（今度は通常の6面、$K=6$）で、真の確率を $\bm{\pi}^* = (0.1, 0.1, 0.1, 0.1, 0.1, 0.5)$（6の目が出やすい）とします。観測数を増やしながら、事後平均と MAP 推定が真の値にどう収束するかを比較します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)
K = 6
pi_true = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.5]) # 6の目が偏っている
alpha_prior = np.ones(K) # 一様事前分布 (add-1 平滑化に対応)
Ns = np.arange(1, 1001)
# 1000回分のサイコロの目を一気に生成
rolls = np.random.choice(K, size=Ns.max(), p=pi_true)
post_mean_6, map_6, freq_6 = [], [], []
for N in Ns:
counts = np.bincount(rolls[:N], minlength=K) # 各目のカウント n_k
alpha_post = alpha_prior + counts # 共役更新
a0 = alpha_post.sum()
post_mean_6.append(alpha_post[5] / a0) # 事後平均
map_6.append((alpha_post[5] - 1) / (a0 - K)) # MAP 推定
freq_6.append(counts[5] / N) # 単純頻度（MLE）
ここまでで各観測数における3種類の推定値（事後平均・MAP・単純頻度）を計算しました。これらを真の値 $0.5$ と比べてプロットします。
plt.figure(figsize=(10, 6))
plt.axhline(0.5, color="k", ls="--", label="true pi_6 = 0.5")
plt.plot(Ns, post_mean_6, label="posterior mean", color="tab:blue")
plt.plot(Ns, map_6, label="MAP", color="tab:orange")
plt.plot(Ns, freq_6, label="frequency (MLE)", color="tab:green", alpha=0.6)
plt.xlabel("number of observations N")
plt.ylabel("estimate of pi_6")
plt.xlim(1, 1000)
plt.ylim(0, 1)
plt.legend()
plt.title("Posterior mean vs MAP vs MLE (true pi_6 = 0.5)")
plt.grid(alpha=0.3)
plt.savefig("mean_vs_map.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから3つのことが読み取れます。第一に、観測数が少ないうち（$N < 50$ あたり）は3つの推定値が大きくばらつき、特に単純頻度（MLE）は激しく上下します。第二に、事後平均は一様事前のため初期は $1/K \approx 0.17$ から出発し、データが増えるにつれて滑らかに真の値 $0.5$ へ近づきます。MAP は事後平均よりやや MLE に近い挙動を示します（$-1$ の補正のため）。第三に、$N$ が大きくなると3つの推定値はすべて真の値 $0.5$ に収束します。少数データでは事前分布による「引き戻し（正則化）」が効いて推定が安定し、データが十分になると事前の影響が消える——これがベイズ推定の望ましい振る舞いです。
最後に、これまでの理論を統合して、偏ったサイコロのベイズ推定を予測分布まで含めて実行します。
偏ったサイコロのベイズ推定（予測分布まで）
実際にサイコロを振った結果から、各目の出る確率の事後平均（＝次に各目が出る予測確率）と、その不確かさ（事後標準偏差）を計算します。ディリクレ分布の分散は $\mathrm{Var}[\pi_k] = \frac{\alpha_k'(\alpha_0′ – \alpha_k’)}{\alpha_0’^2(\alpha_0′ + 1)}$ で与えられます。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(7)
K = 6
pi_true = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.5])
alpha_prior = np.ones(K)
N = 60 # 60回サイコロを振る
rolls = np.random.choice(K, size=N, p=pi_true)
counts = np.bincount(rolls, minlength=K)
alpha_post = alpha_prior + counts # 事後パラメータ alpha + n
a0 = alpha_post.sum()
# 事後予測確率（= 事後平均）と事後標準偏差
pred = alpha_post / a0
var = alpha_post * (a0 - alpha_post) / (a0**2 * (a0 + 1))
std = np.sqrt(var)
print("counts :", counts)
print("naive freq :", np.round(counts / N, 3))
print("predictive :", np.round(pred, 3))
print("posterior std :", np.round(std, 3))
実行すると、たとえば次のような出力が得られます（乱数により多少変動します）。
counts : [ 7 5 6 4 8 30]
naive freq : [0.117 0.083 0.1 0.067 0.133 0.5 ]
predictive : [0.121 0.091 0.106 0.076 0.136 0.47 ]
posterior std : [0.041 0.036 0.039 0.033 0.043 0.063]
この出力から、ベイズ推定の特徴が読み取れます。予測確率（事後平均）は単純頻度に近いものの、すべてのカテゴリに事前の疑似カウント1が足されているため、極端な値がわずかに中心 $1/6$ 側へ引き戻されています。6の目の予測確率は約 $0.47$ で、真の値 $0.5$ にかなり近づいています。さらに重要なのが事後標準偏差で、「この推定がどれくらい信頼できるか」を定量化しています。最後に、この予測確率と不確かさを棒グラフで可視化します。
fig, ax = plt.subplots(figsize=(9, 5))
x = np.arange(1, K + 1)
ax.bar(x, pred, yerr=std, capsize=5, color="tab:blue", alpha=0.7,
label="predictive (posterior mean +/- std)")
ax.plot(x, pi_true, "r*", ms=14, label="true probability")
ax.axhline(1 / K, color="gray", ls=":", label="uniform 1/6")
ax.set_xlabel("die face")
ax.set_ylabel("probability")
ax.set_title(f"Bayesian estimate of a biased die (N={N})")
ax.legend()
ax.grid(alpha=0.3)
plt.savefig("biased_die_bayes.png", dpi=150, bbox_inches="tight")
plt.show()
棒グラフから、ベイズ推定が偏ったサイコロを正しく捉えていることが分かります。6の目（一番右）の予測確率が突出して高く、赤い星で示した真の確率ともよく一致しています。各棒のエラーバー（事後標準偏差）は、$N=60$ という中程度の観測数に対応した不確かさを表しており、6の目はカウントが多いぶん相対的な不確かさが小さいことも見て取れます。点推定だけを返すのではなく「推定値±不確かさ」を返せる点こそ、ディリクレ・カテゴリカル共役モデルの実用上の価値です。
まとめ
本記事では、ディリクレ分布とカテゴリカル分布の共役性を、導出を省略せずに解説しました。
- カテゴリカル分布は「$K$ 個の選択肢から1つ」を表す最も基本的な離散分布で、$N$ 回の試行をまとめると多項分布になり、確率パラメータ $\bm{\pi}$ は単体の上に住みます。
- ディリクレ分布 $\mathrm{Dir}(\bm{\pi} \mid \bm{\alpha}) \propto \prod_k \pi_k^{\alpha_k – 1}$ は、確率の組 $\bm{\pi}$ の上の確率分布であり、ベータ分布の多変量への一般化です。パラメータ $\alpha_k$ は「カテゴリ $k$ の疑似カウント」と解釈できます。
- 共役性：ディリクレ事前にカテゴリカル尤度を掛けると、事後分布は $\mathrm{Dir}(\bm{\pi} \mid \bm{\alpha} + \bm{n})$ になります。ベイズ更新は「$\alpha_k$ に観測カウント $n_k$ を足すだけ」という単純な操作です。
- 事後予測分布：$\bm{\pi}$ を周辺化すると、次にカテゴリ $k$ が出る確率は $\frac{\alpha_k + n_k}{\alpha_0 + N}$ に比例します。これはラプラスの加法的平滑化に対応し、未観測カテゴリの確率がゼロにならないようにします。
- MAPと事後平均：MAP は $\frac{\alpha_k + n_k – 1}{\alpha_0 + N – K}$、事後平均（＝予測確率）は $\frac{\alpha_k + n_k}{\alpha_0 + N}$ で、$N$ が大きくなるとどちらも経験頻度に収束します。
- Python実装：単体上のディリクレ密度の可視化、観測数の増加に伴う事後分布の集中、点推定の収束比較、偏ったサイコロの予測分布までを確認しました。
この共役モデルは、より大きなベイズモデルの構成要素として至るところに登場します。とくに自然言語処理のトピックモデル（LDA）では、文書ごとのトピック分布とトピックごとの単語分布の両方をディリクレ分布で表現し、本記事で導いた共役更新と周辺化を繰り返し使います。
次のステップとして、以下の記事も参考にしてください。
- ベータ分布をわかりやすく解説 — 本記事の2カテゴリ版。共役性の最もシンプルな例です。
- ディリクレ分布をわかりやすく解説 — 分布の形状や性質をさらに詳しく扱います。
- カテゴリカル分布のベイズ推定 — 本記事と密接に関連する応用記事です。
- ベイズの定理の基礎 — 事前・尤度・事後の関係の復習に。