なんとなくベイズ推定の学習と推論を理解したものの、もう少し具体例とコードを動かして確実に理解したい、という人向けに、カテゴリカル分布に対するベイズ推定の学習と推論を数式と具体例、そしてPythonの実装コードを交えてやってみます。
ベイズ推定について、まだ理解できている自信がない人は、下の記事から読んでみることをお勧めします。
本記事の内容
- カテゴリカル分布の定義
- ディリクレ分布を事前分布としたベイズ推定の導出
- サイコロの出目確率の推定問題
- Python での実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
サイコロの問題でベイズ推定を行う
まず問題設定から入ります。出目がまばらなサイコロ（6面体）を考えます。
通常のサイコロであれば、6つの面がそれぞれ均等に $\frac{1}{6}$ ずつ目が出ますが、今回はサイコロの内部に鉛が仕込んであり、特定の出目が出やすくなるような仕掛けがしてあるとします。
このようなサイコロを何回か投げる試行のデータを通して、それぞれの出目が出る確率をベイズ推定で予測します。
カテゴリカル分布の定義
$K$ 個の値をとる離散確率変数を考えるので、カテゴリカル分布を用います。
$$ \begin{equation} p(\bm{s}|\bm{\pi}) = \text{Cat}(\bm{s}|\bm{\pi}) = \prod_{k=1}^{K} \pi_k^{s_k} \end{equation} $$
ここで $\bm{s}$ は1-of-K表現（one-hot表現）で、$K$ 次元のベクトルのうち $k$ 番目の要素だけが1、他が0のベクトルです。パラメータ $\bm{\pi} = (\pi_1, \dots, \pi_K)$ は、
$$ \pi_k \in (0, 1), \quad \sum_{k=1}^{K} \pi_k = 1 $$
を満たします。
事前分布：ディリクレ分布
パラメータ $\bm{\pi}$ の事前分布としてディリクレ分布を使います。
$$ \begin{equation} p(\bm{\pi}) = \text{Dir}(\bm{\pi}|\bm{\alpha}) = \frac{\Gamma(\sum_{k=1}^{K}\alpha_k)}{\prod_{k=1}^{K}\Gamma(\alpha_k)} \prod_{k=1}^{K} \pi_k^{\alpha_k – 1} \end{equation} $$
$\bm{\alpha} = (\alpha_1, \dots, \alpha_K)$ はハイパーパラメータで、$\alpha_k > 0$ です。
尤度関数
$N$ 個のデータ $\mathcal{S} = \{\bm{s}_1, \bm{s}_2, \dots, \bm{s}_N\}$ の尤度関数は、
$$ \begin{align} p(\mathcal{S}|\bm{\pi}) &= \prod_{n=1}^{N} \text{Cat}(\bm{s}_n|\bm{\pi}) = \prod_{n=1}^{N}\prod_{k=1}^{K} \pi_k^{s_{nk}} = \prod_{k=1}^{K} \pi_k^{\sum_{n=1}^{N} s_{nk}} = \prod_{k=1}^{K} \pi_k^{N_k} \end{align} $$
ここで $N_k = \sum_{n=1}^{N} s_{nk}$ はカテゴリ $k$ が観測された回数です。
事後分布の導出
ベイズの定理より、
$$ \begin{align} p(\bm{\pi}|\mathcal{S}) &\propto p(\mathcal{S}|\bm{\pi}) \cdot p(\bm{\pi}) \\ &= \prod_{k=1}^{K} \pi_k^{N_k} \cdot \prod_{k=1}^{K} \pi_k^{\alpha_k – 1} \\ &= \prod_{k=1}^{K} \pi_k^{N_k + \alpha_k – 1} \end{align} $$
これはディリクレ分布のカーネルに一致するので、
$$ \begin{equation} p(\bm{\pi}|\mathcal{S}) = \text{Dir}(\bm{\pi}|\hat{\bm{\alpha}}) \end{equation} $$
ここで $\hat{\alpha}_k = N_k + \alpha_k$ です。ハイパーパラメータ $\alpha_k$ は「仮想的な観測回数」として解釈できます。
予測分布
新しいデータ $\bm{s}_*$ の予測分布は、
$$ \begin{align} p(\bm{s}_*|\mathcal{S}) &= \int p(\bm{s}_*|\bm{\pi}) \cdot p(\bm{\pi}|\mathcal{S}) \, d\bm{\pi} \\ &= E_{\bm{\pi}|\mathcal{S}}[\pi_k] = \frac{\hat{\alpha}_k}{\sum_{j=1}^{K}\hat{\alpha}_j} = \frac{N_k + \alpha_k}{N + \sum_{j=1}^{K}\alpha_j} \end{align} $$
Python での実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import dirichlet
# 偏りのあるサイコロの真の確率
K = 6
true_pi = np.array([0.1, 0.1, 0.15, 0.15, 0.2, 0.3])
true_pi /= true_pi.sum()
# 事前分布のハイパーパラメータ（一様ディリクレ）
alpha_prior = np.ones(K)
# データ生成
np.random.seed(42)
N_total = 200
data = np.random.choice(K, size=N_total, p=true_pi)
# 逐次ベイズ更新
fig, axes = plt.subplots(2, 3, figsize=(16, 9))
axes = axes.flatten()
observations = [0, 5, 20, 50, 100, 200]
x = np.arange(1, K+1)
for idx, N in enumerate(observations):
ax = axes[idx]
obs = data[:N]
# カテゴリごとの観測回数
N_k = np.array([np.sum(obs == k) for k in range(K)])
alpha_post = alpha_prior + N_k
# 事後分布の期待値（予測確率）
pred_pi = alpha_post / alpha_post.sum()
# 棒グラフで比較
width = 0.25
ax.bar(x - width, true_pi, width, color='red', alpha=0.7, label='True')
ax.bar(x, pred_pi, width, color='blue', alpha=0.7, label='Posterior mean')
if N > 0:
freq = N_k / N
ax.bar(x + width, freq, width, color='green', alpha=0.7, label='Frequency')
ax.set_xlabel('Face', fontsize=10)
ax.set_ylabel('Probability', fontsize=10)
ax.set_title(f'N = {N}', fontsize=12)
ax.set_xticks(x)
ax.legend(fontsize=8)
ax.set_ylim(0, 0.45)
ax.grid(True, alpha=0.3, axis='y')
plt.suptitle('Bayesian Estimation of Loaded Dice (Categorical + Dirichlet)', fontsize=14)
plt.tight_layout()
plt.show()
データ数が増えるにつれて、事後平均（青）が真の確率（赤）に近づいていく様子が確認できます。
import numpy as np
import matplotlib.pyplot as plt
# ディリクレ分布からのサンプリングによる不確かさの可視化
np.random.seed(42)
N_total = 200
data = np.random.choice(K, size=N_total, p=true_pi)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
ns = [10, 50, 200]
for idx, N in enumerate(ns):
ax = axes[idx]
obs = data[:N]
N_k = np.array([np.sum(obs == k) for k in range(K)])
alpha_post = alpha_prior + N_k
# ディリクレ分布からサンプリング
samples = np.random.dirichlet(alpha_post, 500)
# 箱ひげ図で不確かさを表示
bp = ax.boxplot([samples[:, k] for k in range(K)], positions=range(1, K+1), widths=0.5)
ax.scatter(range(1, K+1), true_pi, color='red', s=100, zorder=5, marker='*', label='True')
ax.set_xlabel('Face', fontsize=11)
ax.set_ylabel('Probability', fontsize=11)
ax.set_title(f'Posterior Uncertainty (N={N})', fontsize=12)
ax.legend(fontsize=10)
ax.set_ylim(0, 0.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
データが増えるにつれて箱ひげ図の幅が狭くなり、パラメータの不確かさが減少していくことが確認できます。
まとめ
本記事では、カテゴリカル分布に対するベイズ推定について解説しました。
- カテゴリカル分布の共役事前分布はディリクレ分布であり、事後分布もディリクレ分布になる
- 事後分布のパラメータは $\hat{\alpha}_k = N_k + \alpha_k$ で、ハイパーパラメータは仮想的な観測回数と解釈できる
- 予測確率は事後ディリクレ分布の期待値で計算できる
- データが増えるにつれて事後分布が真のパラメータに収束し、不確かさも減少する
次のステップとして、以下の記事も参考にしてください。