アンケートで5つの選択肢それぞれに何人が投票するか。文章の中に各単語が何回ずつ現れるか。ECサイトの訪問者が「購入・カート放棄・即離脱」のどれを何回選ぶか — これらはすべて「K種類の結果があり得る試行をn回繰り返したとき、各結果が何回ずつ出るか」という同じ構造の問題です。この構造を表現する確率分布が多項分布（multinomial distribution）です。
多項分布は、確率・統計の中でも応用範囲がとくに広い分布です。たとえば次のような場面で主役になります。
- 自然言語処理: 文書を「単語の出現回数ベクトル」とみなすBag-of-Wordsモデルや、トピックモデル（LDA）の文書生成部分は多項分布そのもの
- 多クラス分類の評価・検定: 分類器の予測がK個のクラスに散らばる様子や、適合度検定（カイ二乗検定）の帰無分布の土台
- A/Bテストの多変種版（A/B/nテスト）: 複数パターンへのコンバージョンの割れ方のモデル化
- 遺伝学・アンケート分析: 血液型や回答選択肢など、カテゴリ別の度数データ全般
定義式には階乗や総乗記号が並ぶため一見いかつく見えますが、実体は「歪んだサイコロをn回振る」というそれだけの話です。しかもベルヌーイ分布・二項分布・カテゴリカル分布をすべて特殊ケースとして含む「全部入り」の分布なので、多項分布を1つ理解すれば離散分布の基本ファミリーが一望できます。
本記事の内容
- サイコロのアナロジーによる直感的理解と、PMF（確率質量関数）の定義・多項係数の意味
- ベルヌーイ・二項・カテゴリカル分布との包含関係
- 期待値・分散・共分散の導出（共分散がなぜ負になるかの直感も）
- 周辺分布が二項分布になることの証明
- 最尤推定（ラグランジュの未定乗数法）とディリクレ共役事前分布によるベイズ更新
np.random.multinomialによるサンプリングと確率単体上の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
多項分布とは — 歪んだサイコロをn回振る
多項分布を一言でいうと、「出る目の確率が偏ったサイコロをn回振ったとき、各目が何回ずつ出るかの確率を与える分布」です。
まず具体的な状況を見てしまいましょう。次の図は、6つの目の出る確率が $(0.30, 0.10, 0.20, 0.10, 0.10, 0.20)$ という歪んだサイコロを、20回振ったときの様子です。
左が多項分布への「入力」です。各目の出る確率 $\pi_k$（合計すると1）と、振る回数 $n = 20$ を決めます。右が1回の実験の「出力」で、20回振った結果、各目が出た回数のベクトル $\bm{x} = (8, 1, 5, 2, 0, 4)$ が得られました。回数なのですべて0以上の整数で、合計はちょうど20になっています。確率0.30の目1が8回と最多になっている一方、確率0.10の目5は1回も出ていません。回数ベクトルは振り直すたびに変わる — つまり確率的にばらつくわけで、「どの回数ベクトルがどれくらい出やすいか」を教えてくれるのが多項分布です。
ここで多項分布の登場人物を整理しておきます。
- $K$: カテゴリ（出目）の種類数。普通のサイコロなら $K=6$
- $n$: 試行回数（サイコロを振る回数）
- $\bm{\pi} = (\pi_1, \dots, \pi_K)$: 各カテゴリが出る確率。$\sum_k \pi_k = 1$
- $\bm{x} = (x_1, \dots, x_K)$: 各カテゴリが出た回数。$\sum_k x_k = n$
「サイコロを振る」というたとえを使いましたが、冒頭で挙げたアンケート・単語頻度・A/B/nテストもすべて「カテゴリがK種類ある試行をn回繰り返す」構造なので、まったく同じ枠組みに乗ります。
ところで、「コインをn回投げる」のが二項分布だったことを思い出すと、多項分布はその「サイコロ版」に見えます。実はこの直感は正確で、多項分布はおなじみの離散分布たちのちょうど「親玉」に当たります。先に家系図を見ておくと定義式の理解が一気に楽になるので、次にそれを確認します。
ベルヌーイ・二項・カテゴリカル分布との包含関係
ベルヌーイ分布・二項分布・カテゴリカル分布・多項分布の4つは、「試行回数が1回かn回か」「カテゴリが2種類かK種類か」という2つの軸で整理できます。
この図は次のように読みます。
- ベルヌーイ分布（左上）: コイン（2値）を1回投げる
- 二項分布（右上）: ベルヌーイ試行をn回繰り返し、表の回数を数える
- カテゴリカル分布（左下）: サイコロ（K値）を1回振る
- 多項分布（右下）: カテゴリカル試行をn回繰り返し、各目の回数を数える
つまり多項分布は、二項分布を「2値→K値」に広げたものであり、同時にカテゴリカル分布を「1回→n回」に広げたものでもあります。逆にいうと、多項分布で $K=2$ とおけば二項分布が、$n=1$ とおけばカテゴリカル分布が、$K=2, n=1$ とおけばベルヌーイ分布が、それぞれそのまま出てきます。多項分布の性質（期待値や分散など）を1回導出すれば、残り3つの分布の性質は代入するだけで手に入る — これが「親玉を先に押さえる」うまみです。
家系図で立ち位置がわかったので、いよいよ定義式に進みます。記号が多く見えますが、登場人物はさきほど整理した4つだけです。
多項分布の数学的定義
多項分布のPMF（確率質量関数）は、「回数ベクトル $\bm{x}$ が観測される確率」を次の式で与えます。
$$ \mathrm{Multi}(\bm{x} \mid \bm{\pi}, n) = \frac{n!}{x_1! \, x_2! \cdots x_K!} \prod_{k=1}^{K} \pi_k^{x_k} $$
ただし $\bm{x}$ と $\bm{\pi}$ は次の制約を満たします。
$$ \sum_{k=1}^{K} x_k = n, \qquad \sum_{k=1}^{K} \pi_k = 1, \qquad x_k \in \{0, 1, \dots, n\}, \quad \pi_k \ge 0 $$
この式は「確率の部分」と「並べ方の数の部分」の2つの積になっています。分解して理解しましょう。
後半の積: 1通りの出方の確率
後半の $\prod_k \pi_k^{x_k}$ は、「特定の順番でその出目列が出る確率」です。たとえば $K=3$ のサイコロを5回振って「1, 1, 2, 2, 3」という列が出る確率は、各回が独立なので
$$ \pi_1 \cdot \pi_1 \cdot \pi_2 \cdot \pi_2 \cdot \pi_3 = \pi_1^2 \, \pi_2^2 \, \pi_3^1 $$
と、確率の掛け算で書けます。ここで大事なのは、「2, 1, 3, 1, 2」のように順番を入れ替えた列でも、回数ベクトルが同じ $(2, 2, 1)$ なら確率はまったく同じ $\pi_1^2 \pi_2^2 \pi_3$ になることです。
前半の係数: 順番の入れ替えが何通りあるか
多項分布が知りたいのは「特定の順番」ではなく「回数ベクトル $\bm{x}$」の確率なので、同じ $\bm{x}$ を与える順番のパターンを全部足し合わせる必要があります。各パターンの確率は等しいので、「パターン数 × 1パターンの確率」を計算すればよいわけです。このパターン数が前半の係数
$$ \frac{n!}{x_1! \, x_2! \cdots x_K!} $$
で、多項係数と呼ばれます。「n個のものを、区別のつかない $x_1$ 個、$x_2$ 個、…のグループに分ける並べ方の数」です。$n!$ がn回の試行の並べ方すべて、分母の $x_k!$ が「同じ出目どうしの入れ替えは区別しない」ための割り算です。二項係数 $\binom{n}{x} = \frac{n!}{x!(n-x)!}$ のK種類版だと思えば見慣れた形のはずです。
まとめると、多項分布のPMFは
$$ \underbrace{\frac{n!}{x_1! \cdots x_K!}}_{\text{その回数になる並べ方の数}} \times \underbrace{\pi_1^{x_1} \cdots \pi_K^{x_K}}_{\text{並べ方1通りの確率}} $$
という、ごく自然な「場合の数 × 確率」の構造をしています。記号の多さに反して、やっていることは高校の確率の「反復試行」と同じです。
定義がわかったところで、実際に数値を入れて確率を計算し、目で確かめてみましょう。
具体例 — 手計算とscipyでの確認
$K=3$ のくじ引き（当たり確率 $\bm{\pi} = (0.5, 0.3, 0.2)$）を $n=5$ 回引いて、結果が $\bm{x} = (2, 2, 1)$ になる確率を計算します。定義式に代入すると
$$ P(\bm{x}) = \frac{5!}{2! \, 2! \, 1!} \times 0.5^2 \times 0.3^2 \times 0.2^1 $$
まず多項係数を計算します。$5! = 120$、$2! \, 2! \, 1! = 4$ なので $120 / 4 = 30$ 通りです。次に確率部分は $0.25 \times 0.09 \times 0.2 = 0.0045$。掛け合わせて
$$ P(\bm{x}) = 30 \times 0.0045 = 0.135 $$
つまり13.5%です。scipyで検算してみましょう。
import numpy as np
from scipy.stats import multinomial
n, pi = 5, [0.5, 0.3, 0.2]
dist = multinomial(n, pi)
print(round(dist.pmf([2, 2, 1]), 5)) # 0.135
print(round(dist.pmf([5, 0, 0]), 5)) # 0.03125 (= 0.5^5)
print(round(dist.pmf([0, 0, 5]), 5)) # 0.00032 (= 0.2^5)
手計算どおり0.135が得られました。比較用に出した「全部カテゴリ1」の確率は $0.5^5 = 0.03125$、「全部カテゴリ3」は $0.2^5 = 0.00032$ で、これらは並べ方が1通りしかない（多項係数が1）ため、確率部分がそのまま答えになります。最頻の組み合わせ $(2, 2, 1)$ あたりと比べると、極端な出方がいかに起きにくいかがわかります。
1点ずつの計算だけでなく、PMFの「全体の景色」も見ておきましょう。$K=3$ の場合は $x_3 = n – x_1 – x_2$ と自動的に決まるので、$(x_1, x_2)$ の2次元平面にPMFをヒートマップとして描けます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multinomial
n, pi = 10, [0.5, 0.3, 0.2]
P = np.full((n + 1, n + 1), np.nan)
for x1 in range(n + 1):
for x2 in range(n + 1 - x1):
P[x2, x1] = multinomial(n, pi).pmf([x1, x2, n - x1 - x2])
plt.figure(figsize=(7, 5.5))
plt.imshow(P, origin='lower', cmap='viridis')
plt.colorbar(label='確率')
plt.xlabel('出目1の回数 $x_1$')
plt.ylabel('出目2の回数 $x_2$')
plt.show()
この図から2つのことが読み取れます。第一に、確率が定義されるのは $x_1 + x_2 \le 10$ の三角形の領域だけです（それ以外は回数の合計が10を超えるのであり得ません）。第二に、確率の山のてっぺん（最頻値）は $\bm{x} = (5, 3, 2)$、つまりちょうど $n\bm{\pi} = (10 \times 0.5, 10 \times 0.3, 10 \times 0.2)$ の位置にあります。「確率に比例した回数が一番出やすい」という直感そのままです。
この「$n\pi_k$ が中心」という観察は、偶然ではありません。次のセクションで、期待値が本当に $n\pi_k$ になることを、分散・共分散とまとめて導出します。
期待値・分散・共分散の導出
結論を先にまとめます。
$$ E[x_k] = n\pi_k, \qquad V[x_k] = n\pi_k(1 – \pi_k), \qquad \mathrm{Cov}(x_i, x_j) = -n\pi_i\pi_j \;\; (i \neq j) $$
これらをまとめて導く強力なテクニックが、指示変数（インディケータ）への分解です。
準備: 回数を「1回ごとの試行」に分解する
$t$ 回目の試行でカテゴリ $k$ が出たら1、出なければ0を取る指示変数を $I_k^{(t)}$ と書きます。すると「カテゴリ $k$ が出た回数」は、単に各回の指示変数の合計です。
$$ x_k = \sum_{t=1}^{n} I_k^{(t)} $$
指示変数は0か1しか取らないため、期待値が簡単に計算できます。
$$ E[I_k^{(t)}] = 1 \cdot \pi_k + 0 \cdot (1 – \pi_k) = \pi_k $$
期待値の導出
期待値の線形性（和の期待値は期待値の和）を使うと、一瞬で求まります。
$$ E[x_k] = E\!\left[\sum_{t=1}^{n} I_k^{(t)}\right] = \sum_{t=1}^{n} E[I_k^{(t)}] = n\pi_k $$
確率0.3の目を10回振れば平均3回出る、という当たり前の結果が、きちんと式から出てきました。
分散の導出
各試行は独立なので、和の分散は分散の和になります。指示変数 $I_k^{(t)}$ はベルヌーイ変数なので、その分散は $\pi_k(1 – \pi_k)$ です（$E[I^2] = E[I] = \pi_k$ より $V[I] = \pi_k – \pi_k^2$）。よって
$$ V[x_k] = \sum_{t=1}^{n} V[I_k^{(t)}] = n\pi_k(1 – \pi_k) $$
これは二項分布の分散と同じ形です。後で見るように $x_k$ 単独の分布は二項分布そのものなので、当然の一致です。
共分散の導出 — なぜ負になるのか
多項分布で最も面白いのが、異なるカテゴリの回数 $x_i$ と $x_j$ の共分散です。共分散も試行ごとに分解できます。異なる試行の指示変数は独立で共分散0なので、生き残るのは「同じ試行 $t$ の中での $I_i^{(t)}$ と $I_j^{(t)}$ の共分散」だけです。
$$ \mathrm{Cov}(x_i, x_j) = \sum_{t=1}^{n} \mathrm{Cov}\!\left(I_i^{(t)}, I_j^{(t)}\right) $$
ここで核心の観察をします。1回の試行で出る目は1つだけなので、$I_i^{(t)}$ と $I_j^{(t)}$ が同時に1になることはあり得ません。つまり積は常に0です。
$$ E\!\left[I_i^{(t)} I_j^{(t)}\right] = 0 $$
共分散の定義 $\mathrm{Cov}(A, B) = E[AB] – E[A]E[B]$ に代入すると
$$ \mathrm{Cov}\!\left(I_i^{(t)}, I_j^{(t)}\right) = 0 – \pi_i \pi_j = -\pi_i \pi_j $$
これをn回分合計して、求める共分散が得られます。
$$ \mathrm{Cov}(x_i, x_j) = -n\pi_i\pi_j $$
共分散は必ず負です。直感的にも納得できます。合計回数はnに固定されているので、回数の取り合いはゼロサムゲームです。ある試行で目1が出れば、その回は目2には絶対に回ってきません。目1がたくさん出た実験では、残りの枠が減るぶん目2は少なくなりがち — これが負の相関の正体です。
実際にサンプルを描いて確かめてみましょう。$n=20$、$\bm{\pi}=(0.4, 0.4, 0.2)$ から3000個の回数ベクトルを生成し、$x_1$ と $x_2$ を散布図にします。
点群が右下がりに分布しており、$x_1$ が多い実験ほど $x_2$ が少ない傾向がはっきり見えます。理論値を計算すると $\mathrm{Cov}(x_1, x_2) = -20 \times 0.4 \times 0.4 = -3.2$、相関係数にすると $-3.2 / 4.8 = -2/3 \approx -0.67$ で、図中の標本相関係数 $-0.66$ とよく一致しています。点線の「上限 $x_1 + x_2 = n$」が、回数の取り合いというゼロサム構造を視覚的に示しています。
3つのモーメントが揃ったので、シミュレーションでまとめて検証します。
import numpy as np
rng = np.random.default_rng(0)
n, pi = 20, np.array([0.4, 0.4, 0.2])
S = rng.multinomial(n, pi, size=50000) # shape (50000, 3)
print("期待値: 理論", n * pi, " 実測", S.mean(axis=0).round(3))
print("分散: 理論", n * pi * (1 - pi), " 実測", S.var(axis=0).round(3))
print("共分散: 理論", -n * pi[0] * pi[1], " 実測", np.cov(S[:, 0], S[:, 1])[0, 1].round(3))
5万サンプルの実測値は、期待値 $(8, 8, 4)$、分散 $(4.8, 4.8, 3.2)$、共分散 $-3.2$ という理論値と小数第1位までぴったり一致します。導出した公式が正しいことが数値的にも確認できました。
分散の導出で「$x_k$ 単独の分布は二項分布」とさらっと述べました。これは多項分布の理解を大きく助ける事実なので、次にきちんと確認します。
周辺分布は二項分布になる
多項分布の回数ベクトルから1成分 $x_1$ だけに注目したとき（他の成分は気にしない＝周辺化したとき）、$x_1$ はどんな分布に従うでしょうか。
答えは二項分布 $\mathrm{Bin}(n, \pi_1)$ です。証明は、視点の切り替え1つで済みます。カテゴリ1だけに興味があるなら、各試行の結果を「カテゴリ1が出た（確率 $\pi_1$）／それ以外が出た（確率 $1 – \pi_1$）」の2値に潰して見ても、$x_1$ の値は変わりません。2値の試行をn回繰り返した成功回数は、定義そのものずばり二項分布です。よって
$$ x_1 \sim \mathrm{Bin}(n, \pi_1) $$
形式的にPMFの和（$x_2, \dots, x_K$ についての周辺化）を多項定理で計算しても同じ結果が得られますが、この「2択に潰す」議論のほうが本質を突いています。シミュレーションでも確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
rng = np.random.default_rng(0)
n, pi = 30, [0.5, 0.3, 0.2]
S = rng.multinomial(n, pi, size=20000)
plt.figure(figsize=(8, 4.5))
plt.hist(S[:, 0], bins=np.arange(-0.5, n + 1.5), density=True,
alpha=0.6, label='多項分布からの$x_1$')
vals = np.arange(n + 1)
plt.plot(vals, binom(n, pi[0]).pmf(vals), 'o-', color='red',
label='二項分布 Bin(30, 0.5)')
plt.xlabel('$x_1$'); plt.ylabel('確率'); plt.legend()
plt.show()
多項分布から生成した $x_1$ のヒストグラム（青）の上に、二項分布 $\mathrm{Bin}(30, 0.5)$ のPMF（赤）がぴったり重なっています。多項分布は「K本の二項分布を、合計nという制約（と負の共分散）で束ねたもの」というイメージが持てると、かなり見通しがよくなります。
ここまでは「$\bm{\pi}$ が既知」という前提で分布の性質を調べてきました。しかし実務では話が逆で、観測された回数 $\bm{x}$ から $\bm{\pi}$ を推定したいことがほとんどです。ここからは推定の話に入ります。
最尤推定 — ラグランジュの未定乗数法で π を求める
歪んだサイコロをn回振って回数 $\bm{x}$ を観測したとき、各目の確率 $\bm{\pi}$ をどう見積もるべきでしょうか。直感的には「出た比率 $x_k / n$ をそのまま使えばよい」と思うはずです。この直感が最尤推定の答えと一致することを、きちんと導出します。
ゴールは、対数尤度を $\sum_k \pi_k = 1$ の制約のもとで最大化する $\bm{\pi}$ を求めることです。多項分布の対数尤度は
$$ \log L(\bm{\pi}) = \log n! – \sum_{k=1}^{K} \log x_k! + \sum_{k=1}^{K} x_k \log \pi_k $$
前の2項は $\bm{\pi}$ を含まない定数なので、最大化すべきは $\sum_k x_k \log \pi_k$ だけです。制約付き最大化なので、ラグランジュの未定乗数法を使います。ラグランジュ関数を
$$ \mathcal{L}(\bm{\pi}, \lambda) = \sum_{k=1}^{K} x_k \log \pi_k + \lambda \left(1 – \sum_{k=1}^{K} \pi_k\right) $$
とおき、$\pi_k$ で偏微分して0とおきます。
$$ \frac{\partial \mathcal{L}}{\partial \pi_k} = \frac{x_k}{\pi_k} – \lambda = 0 \quad \Longrightarrow \quad \pi_k = \frac{x_k}{\lambda} $$
未知の $\lambda$ は、制約条件に代入すれば決まります。
$$ \sum_{k=1}^{K} \pi_k = \sum_{k=1}^{K} \frac{x_k}{\lambda} = \frac{n}{\lambda} = 1 \quad \Longrightarrow \quad \lambda = n $$
したがって、最尤推定値は
$$ \hat{\pi}_k = \frac{x_k}{n} $$
つまり「観測された比率そのもの」です。直感どおりの結果が、制約付き最適化からきれいに出てきました。
ところで $\bm{\pi}$ の取り得る値の全体は、「成分が非負で合計1」という条件を満たすベクトルの集合で、確率単体（probability simplex）と呼ばれます。$K=3$ なら正三角形の内部として描けるので、対数尤度の景色を直接眺められます。
左が確率単体の見取り図です。三角形の頂点が「その目しか出ない」極端な $\bm{\pi}$、重心が一様な $(1/3, 1/3, 1/3)$ に対応し、内部のあらゆる点が1つの「あり得るサイコロ」を表します。右は観測 $\bm{x} = (5, 3, 2)$（$n=10$）に対する対数尤度を単体上で色塗りしたもので、最も明るい点（赤い×）がちょうど比率 $\bm{x}/n = (0.5, 0.3, 0.2)$ に一致しています。導出した式の幾何的な姿です。
ただし最尤推定には弱点があります。たとえば10回振って目3が一度も出なかったら $\hat{\pi}_3 = 0$、つまり「目3は永遠に出ない」と断定してしまいます。データが少ないときのこの過信を防ぐのがベイズ的なアプローチで、多項分布には相性抜群の相棒がいます。
ディリクレ共役事前分布とベイズ更新
ベイズ推定では、$\bm{\pi}$ 自体を確率変数とみなして事前分布を置き、データで更新して事後分布を得ます。多項分布の尤度に対しては、ディリクレ分布を事前分布に選ぶと計算が劇的に簡単になります。ディリクレ分布は確率単体上の分布で、密度は
$$ \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha}) \propto \prod_{k=1}^{K} \pi_k^{\alpha_k – 1} $$
という形です（詳細はディリクレ分布の解説記事を参照）。なぜこれが「相性抜群」なのか、事後分布を計算するとわかります。ベイズの定理より、事後分布は尤度と事前分布の積に比例します。
$$ p(\bm{\pi} \mid \bm{x}) \propto \underbrace{\prod_{k} \pi_k^{x_k}}_{\text{多項尤度}} \times \underbrace{\prod_{k} \pi_k^{\alpha_k – 1}}_{\text{ディリクレ事前}} = \prod_{k} \pi_k^{x_k + \alpha_k – 1} $$
最後の式は、パラメータを $\alpha_k + x_k$ に置き換えたディリクレ分布の形そのものです。つまり
$$ p(\bm{\pi} \mid \bm{x}) = \mathrm{Dir}(\bm{\pi} \mid \bm{\alpha} + \bm{x}) $$
事前と事後が同じ分布族に収まる — この性質を共役と呼びます。更新ルールは「事前パラメータ $\bm{\alpha}$ に観測回数 $\bm{x}$ を足すだけ」という驚くべき簡単さです。$\alpha_k$ は「観測前から持っている仮想的な観測回数」と解釈でき、データが少ないうちは推定を滑らかにし（$\hat{\pi}_k = 0$ のような断定を防ぎ）、データが増えるほど影響が薄れます。
更新の様子を確率単体上で見てみましょう。真のサイコロを $\bm{\pi} = (0.5, 0.3, 0.2)$ とし、事前分布 $\mathrm{Dir}(2, 2, 2)$ から始めて、$n=10$ で $\bm{x} = (5, 3, 2)$、$n=100$ で $\bm{x} = (50, 30, 20)$ が観測されたとします。
import numpy as np
alpha = np.array([2.0, 2.0, 2.0]) # 事前分布 Dir(2, 2, 2)
x10 = np.array([5, 3, 2]) # n=10 の観測
x100 = np.array([50, 30, 20]) # n=100 の観測
post10 = alpha + x10 # 足すだけでベイズ更新完了
post100 = alpha + x100
print("n=10 の事後:", post10, " 事後平均:", (post10 / post10.sum()).round(3))
# n=10 の事後: [7. 5. 4.] 事後平均: [0.438 0.312 0.25]
print("n=100の事後:", post100, " 事後平均:", (post100 / post100.sum()).round(3))
# n=100の事後: [52. 32. 22.] 事後平均: [0.491 0.302 0.208]
事後平均は $n=10$ の段階で $(0.438, 0.312, 0.250)$ と、最尤推定値 $(0.5, 0.3, 0.2)$ よりも一様分布側へ引き戻された値になります。これが事前分布による平滑化です。$n=100$ では $(0.491, 0.302, 0.208)$ と、事前分布の影響がほぼ消えて観測比率に近づいています。
左の事前分布は単体の中央付近に広くぼんやりした山で、「どんなサイコロかまだよくわからない」状態を表しています。中央の $n=10$ 観測後は山が真値方向に動き始め、右の $n=100$ 観測後には真の $\bm{\pi}$（赤い×）のすぐ近くに鋭く集中しています。データが増えるほど確信が強まる過程が、パラメータの足し算だけで表現できているのがディリクレ共役の威力です。この枠組みの詳細はディリクレ・カテゴリカルの共役ベイズ解説で深掘りしています。
最後に、ここまで何度も使ってきたサンプリング関数 np.random.multinomial を正面から扱い、推定の収束を可視化して締めくくります。
np.random.multinomialと大数の法則の可視化
NumPyでは np.random.multinomial（または Generator.multinomial）で多項分布からのサンプルを生成できます。使い方は「試行回数n・確率ベクトル・生成数」を渡すだけです。
import numpy as np
rng = np.random.default_rng(42)
pi = [0.5, 0.3, 0.2]
print(rng.multinomial(10, pi)) # [6 3 1] — 1回の実験（n=10の回数ベクトル）
print(rng.multinomial(10, pi, size=3)) # 3回の実験をまとめて生成（各行の合計は10）
# [[7 1 2]
# [3 2 5]
# [6 2 2]]
返り値の各行は合計がちょうどnになる回数ベクトルです。乱数生成の基本はNumPyのランダムサンプリング解説も参考にしてください。
仕上げに、最尤推定値 $\hat{\bm{\pi}} = \bm{x}/n$ がnを増やすとどう振る舞うかを、確率単体上で観察します。$n = 10, 100, 1000$ のそれぞれで実験を300回繰り返し、得られた経験比率を単体に打点します。
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(42)
true_pi = np.array([0.5, 0.3, 0.2])
corners = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3) / 2]]) # 単体→平面の変換
fig, axes = plt.subplots(1, 3, figsize=(13, 4.2))
for ax, n in zip(axes, [10, 100, 1000]):
ratios = rng.multinomial(n, true_pi, size=300) / n # 経験比率 x/n
xy = ratios @ corners
ax.scatter(xy[:, 0], xy[:, 1], s=10, alpha=0.4)
ax.plot(*(true_pi @ corners), 'rx', markersize=11)
tri = np.vstack([corners, corners[0]])
ax.plot(tri[:, 0], tri[:, 1], 'k-', lw=1)
ax.set_title(f'n = {n}'); ax.set_aspect('equal'); ax.axis('off')
plt.show()
$n=10$ では経験比率が単体の広い範囲に散らばっています（取り得る値が0.1刻みなので格子状です）。$n=100$ で真値の周りにまとまり、$n=1000$ では赤い×のごく近傍に集中します。これは大数の法則そのものです。さらに各成分の標準偏差は $\sqrt{\pi_k(1-\pi_k)/n}$、つまり $1/\sqrt{n}$ のオーダーで縮むので、「精度を10倍にしたければデータは100倍必要」という統計の経験則も、この図のばらつきの縮み方（10 → 1000 でばらつき約1/10）から読み取れます。
まとめ
本記事では、多項分布について解説しました。要点をまとめます。
- 多項分布は「歪んだサイコロをn回振ったとき、各目が何回ずつ出るか」の分布。PMFは「並べ方の数（多項係数）× 1通りの確率」という自然な構造を持つ
- ベルヌーイ・二項・カテゴリカル分布は、多項分布の $K=2$ や $n=1$ の特殊ケース。多項分布は離散基本分布の「親玉」
- モーメントは $E[x_k] = n\pi_k$、$V[x_k] = n\pi_k(1-\pi_k)$、$\mathrm{Cov}(x_i, x_j) = -n\pi_i\pi_j$。共分散が必ず負なのは、合計nの枠を取り合うゼロサム構造のため
- 1成分だけ見た周辺分布は二項分布 $\mathrm{Bin}(n, \pi_k)$（「その目か、それ以外か」の2択に潰せる）
- 最尤推定値はラグランジュの未定乗数法から $\hat{\pi}_k = x_k/n$。ベイズではディリクレ事前分布が共役で、更新は $\bm{\alpha} + \bm{x}$ と足すだけ
多項分布が手に馴染むと、トピックモデルや混合モデルなど「カテゴリの世界」のベイズモデリングが一気に読めるようになります。次のステップとして、以下の記事も参考にしてください。