「コインを5回投げて表が3回出た。このコインの表が出る真の確率はいくつだろうか？」
頻度論的に答えれば $3/5 = 0.6$ となります。しかしこの答えは少し奇妙です。たった5回の試行で、本当に「0.6」と断言できるのでしょうか。50回投げて30回表が出たときも答えは同じ0.6ですが、こちらのほうがずっと自信が持てるはずです。さらに、5回投げて全部表だったら確率は1.0なのでしょうか？「裏が出ない無敵のコイン」と本気で信じる人はいないでしょう。
この違和感の正体は、頻度論的推定は「点」しか返さないことにあります。「真値はおそらく0.6だが、0.4や0.8の可能性も十分ある」「データが少ないほど不確かである」という、人間が直感的に持つ確信度の濃淡が表現できないのです。ベイズ推論は、この欠落を埋めます。点推定の代わりに確率分布そのものを推定することで、「どの値がどれくらい確からしいか」をまるごと表現できる枠組みです。
この考え方は、現代の機械学習やデータ分析のあちこちで活躍しています。WebサイトのA/Bテストでは、「新デザインのクリック率は本当に旧デザインより高いか」を、限られたサンプルから事後分布として推定し、リスクを許容しながら意思決定します。スパムフィルタは、メールに含まれる単語の事後確率を逐次更新してスパム判定を更新していきます。ガウス過程回帰やベイズ最適化はハイパーパラメータ探索やロボット制御の標準ツールであり、いずれもベイズ推論を骨組みとしています。さらに、近年話題の確率的深層学習（Bayesian Deep Learning）は、ニューラルネットの重みに事後分布を入れて予測の不確実性を定量化する試みです。
本記事の内容
- ベイズの定理と、事前分布・尤度・事後分布・予測分布の役割
- ベルヌーイ-ベータ共役モデルでの完全な解析的導出
- MAP推定と最尤推定との対応関係、事前分布の正則化的解釈
- 共役でない場合のMCMC・変分推論のさわり
- 信頼区間と信用区間の違い — 似て非なる2つの区間
- PythonとPyMCによる実装、解析解とMCMCの突き合わせ
- 階層ベイズとベイズ最適化への入口
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
直感 — 「信念をデータで更新する」
数式に入る前に、ベイズ推論が何をしている儀式なのかをイメージで掴んでおきましょう。
ある朝、あなたは知らないコインを手渡されたとします。投げる前、このコインが「公平なコイン」だという気がしているかもしれません。あるいは、相手が手品師だったら「イカサマの確率がそこそこある」と疑うかもしれません。この「投げる前の信念」を確率分布で表したものが、事前分布（prior）です。何の情報もなく「公平っぽい」と思っている状態は、$\mu = 0.5$ あたりに山がある分布として表せますし、まったく見当がつかないなら「$\mu$ は0から1の間にどれも等しくありそう」という一様分布で表せます。
次にあなたはコインを何回か投げます。表が多く出れば「このコインは表に偏っているのでは？」と疑い、半々に近ければ「やっぱり公平か」と確信を強めます。この観測データに照らして信念を更新する操作こそ、ベイズ推論の中心にある仕組みです。更新された信念は事後分布（posterior）と呼ばれます。
- 事前分布 = データを見る前の信念
- データを取得する
- 事後分布 = データを見た後の信念
ベイズ推論は、この「信念の更新」を数学的にきちんと書き下したものです。そして驚くべきことに、その更新ルールはたった一つの式 — ベイズの定理 — に集約されます。次の節でその式と中身を確認しましょう。
ベイズの定理 — 事後分布を組み立てる方程式
条件付き確率からベイズの定理へ
確率論の基本道具である条件付き確率から出発します。事象 $A, B$ について、同時確率 $P(A, B)$ は二通りに書けます。
$$ P(A, B) = P(A \mid B)\, P(B) = P(B \mid A)\, P(A) $$
中辺と右辺を等号で結び、$P(B)$ で割ると、
$$ P(A \mid B) = \frac{P(B \mid A)\, P(A)}{P(B)} $$
これがベイズの定理です。式自体は1行ですが、$A$ をパラメータ $\theta$、$B$ を観測データ $\mathcal{D}$ に置き換えると、ベイズ推論の中心方程式が現れます。
$$ \boxed{\; p(\theta \mid \mathcal{D}) = \frac{p(\mathcal{D} \mid \theta)\, p(\theta)}{p(\mathcal{D})} \;} $$
各項には固有名詞が付いています。
| 記号 | 名前 | 意味 |
|---|---|---|
| $p(\theta)$ | 事前分布（prior） | データを見る前のパラメータへの信念 |
| $p(\mathcal{D} \mid \theta)$ | 尤度（likelihood） | パラメータが $\theta$ のときデータ $\mathcal{D}$ が観測される確率（の関数） |
| $p(\theta \mid \mathcal{D})$ | 事後分布（posterior） | データを見た後のパラメータへの信念 |
| $p(\mathcal{D})$ | エビデンス（evidence）/ 周辺尤度 | データ自体が観測される確率（規格化定数） |
エビデンスは「規格化」のために存在する
エビデンス $p(\mathcal{D})$ は、$\theta$ を積分消去したものです。
$$ p(\mathcal{D}) = \int p(\mathcal{D} \mid \theta)\, p(\theta)\, d\theta $$
これは $\theta$ に依らない定数なので、$\theta$ の関数として事後分布の形を見るときは無視できます。
$$ p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta)\, p(\theta) $$
この比例関係は実用上きわめて重要です。多くの場合、エビデンスの積分は解析的に解けず数値計算も困難ですが、形が分かれば事後分布の特定には十分なケースが大半だからです。MCMCのようなサンプリング法も、この「規格化定数を計算せずに事後分布から標本を取る」という発想を活用しています。
ベイズの定理の数値的な感触
抽象的な式だけだと掴みにくいので、具体的な数で味わってみましょう。「ある病気の罹患率は $1/1000$、検査の感度（罹患者を陽性と判定する確率）は $0.99$、特異度（健常者を陰性と判定する確率）は $0.95$」というよく出る例題です。陽性反応が出たとき、本当に罹患している確率は？
$$ P(\text{罹患} \mid \text{陽性}) = \frac{P(\text{陽性} \mid \text{罹患})\, P(\text{罹患})}{P(\text{陽性})} $$
分母 $P(\text{陽性})$ は罹患者の陽性 + 健常者の偽陽性の和で、
$$ P(\text{陽性}) = 0.99 \times 0.001 + 0.05 \times 0.999 \approx 0.0509 $$
これを代入すると、
$$ P(\text{罹患} \mid \text{陽性}) = \frac{0.99 \times 0.001}{0.0509} \approx 0.0194 $$
陽性反応にもかかわらず、本当に病気である確率は約 1.9% にすぎません。直感的にはずいぶん低く感じますが、これは事前確率 $P(\text{罹患}) = 0.001$ が小さいために偽陽性が支配するからです。事前分布を無視するとひどい誤推論をやらかすという、ベイズ的思考の威力（と落とし穴）が同時に見える例題です。
10万人を仮想的に並べた2×2の図で見ると、なぜ陽性予測値がこれほど低くなるのかが一目で分かります。罹患者100人中99人が真陽性になる一方、健常者99,900人のうち5%にあたる4,995人が偽陽性になります。同じ「陽性」のグループに、罹患者99人と健常者4,995人が混じっているため、陽性者全体に占める罹患者比率は約 1.9% に留まります。事前確率の影響を視覚化する典型例です。
ここまでで、ベイズの定理という「事後分布を組み立てる工具」が手に入りました。では、機械学習の文脈ではこの工具をどう使うのか、もう少し体系立てて見ていきましょう。
学習と予測 — ベイズ流の二段ロケット
機械学習における学習と予測のサイクルを、ベイズ推論はこう書きます。
ステップ1: 学習（事後分布を求める）
データ $\mathcal{D} = \{x_1, \dots, x_N\}$ が独立同分布（i.i.d.）でモデル $p(x \mid \theta)$ から生成されたと仮定すると、尤度は積で書けます。
$$ p(\mathcal{D} \mid \theta) = \prod_{n=1}^{N} p(x_n \mid \theta) $$
これと事前分布 $p(\theta)$ を組み合わせて、事後分布が得られます。
$$ p(\theta \mid \mathcal{D}) \propto \left[\prod_{n=1}^{N} p(x_n \mid \theta)\right] p(\theta) $$
頻度論なら $\hat{\theta}_{\text{MLE}} = \arg\max_\theta p(\mathcal{D} \mid \theta)$ という点を返すところを、ベイズ流では分布を返す。これが本質的な違いです。
ステップ2: 予測（予測分布で次の観測を予想する）
学んだ事後分布を使って、未来の観測 $\hat{x}$ を予測するには次の式を使います。
$$ \boxed{\; p(\hat{x} \mid \mathcal{D}) = \int p(\hat{x} \mid \theta)\, p(\theta \mid \mathcal{D})\, d\theta \;} $$
これを予測分布（posterior predictive distribution）と呼びます。導出は次の通り。同時分布 $p(\hat{x}, \theta \mid \mathcal{D})$ を $\theta$ について積分消去すれば
$$ p(\hat{x} \mid \mathcal{D}) = \int p(\hat{x}, \theta \mid \mathcal{D})\, d\theta = \int p(\hat{x} \mid \theta, \mathcal{D})\, p(\theta \mid \mathcal{D})\, d\theta $$
ここで、未来観測 $\hat{x}$ がパラメータ $\theta$ のもとで過去データ $\mathcal{D}$ と条件独立だと仮定すると（i.i.d.の自然な帰結）$p(\hat{x} \mid \theta, \mathcal{D}) = p(\hat{x} \mid \theta)$ となり、目的の式に着きます。
頻度論的な「プラグイン予測」が単に $p(\hat{x} \mid \hat{\theta}_{\text{MLE}})$ で済ますのに対し、ベイズ予測は「$\theta$ がどの値であるかも分からないのだから、可能性のある全ての $\theta$ について加重平均を取ろう」と考えます。これによりパラメータ不確実性が予測の不確実性に自然に伝播するのが、ベイズ予測の最大の強みです。
ここまで一般論を見てきましたが、抽象的すぎてピンと来ない部分もあるはずです。次の節からは、最も透明な具体例 — ベルヌーイ・ベータ共役モデル — で全てを手計算してみましょう。
ベルヌーイ-ベータ共役モデル — 完全解析
問題設定
裏表のあるコインを $N$ 回投げて、結果 $\bm{X} = \{x_1, \dots, x_N\}$（$x_n \in \{0, 1\}$、$1$ が表）が得られました。次の試行で表が出る確率を推定したい。
- モデル: 各試行は独立にベルヌーイ分布 $p(x \mid \mu) = \mu^x (1-\mu)^{1-x}$ から生成。$\mu \in [0, 1]$ が表確率。
- 事前分布: $\mu$ にはベータ分布 $p(\mu) = \mathrm{Beta}(\mu \mid a, b)$ を置きます。
ベータ分布を選ぶのは偶然ではありません。ベルヌーイ尤度との掛け算が再びベータ分布の形になるという美しい性質 — 共役性（conjugacy） — があるからです。これからその仕組みを目撃します。
事前分布: ベータ分布
ベータ分布の確率密度関数は
$$ \mathrm{Beta}(\mu \mid a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\, \mu^{a-1}(1-\mu)^{b-1} $$
ここで $\Gamma$ はガンマ関数で、$a, b > 0$ が形状パラメータです。$a = b = 1$ なら一様分布、$a = b = 0.5$ なら両端が立ち上がるU字、$a = b$ で大きければ $\mu = 0.5$ を中心とする山型になります。
尤度
$N$ 個の独立観測の尤度は、ベルヌーイの確率質量関数の積です。
$$ p(\bm{X} \mid \mu) = \prod_{n=1}^{N} \mu^{x_n}(1-\mu)^{1-x_n} = \mu^{\sum_n x_n}(1-\mu)^{N – \sum_n x_n} $$
ここで $k = \sum_{n=1}^{N} x_n$（表の回数）と置くと
$$ p(\bm{X} \mid \mu) = \mu^k (1-\mu)^{N-k} $$
簡潔な形にまとまりました。これは表の回数 $k$ という十分統計量だけで尤度が決まることを示しています。
事後分布の導出
ベイズの定理より、事後分布は尤度と事前分布の積に比例します。
$$ p(\mu \mid \bm{X}) \propto p(\bm{X} \mid \mu)\, p(\mu) = \mu^k (1-\mu)^{N-k} \cdot \mu^{a-1}(1-\mu)^{b-1} $$
指数法則で同じ底の冪をまとめると
$$ p(\mu \mid \bm{X}) \propto \mu^{k + a – 1}(1-\mu)^{N – k + b – 1} $$
ここで重要な観察があります。この式は再びベータ分布 $\mathrm{Beta}(\mu \mid \hat{a}, \hat{b})$ の形をしています。よって規格化定数を補えば、
$$ \boxed{\; p(\mu \mid \bm{X}) = \mathrm{Beta}(\mu \mid \hat{a}, \hat{b}), \quad \hat{a} = a + k, \;\; \hat{b} = b + (N – k) \;} $$
これが共役性の威力です。事前分布のパラメータ $(a, b)$ に「表の回数 $k$」と「裏の回数 $N – k$」をそれぞれ足すだけで事後分布が得られます。難しい積分や数値計算は一切要りません。
直感的解釈もシンプルです。事前分布の $(a, b)$ は「過去に $a-1$ 回表、$b-1$ 回裏を見たような仮想的な経験」を表していると考えられます。データを観測すると、その仮想経験に実観測が加算されていく。$a, b$ は「擬似カウント（pseudo-count）」と呼ばれることもあります。
事前→データ→事後の3段階を密度プロットで並べると、共役性が「数を足すだけ」という性質であることが視覚的に確認できます。事前 $\mathrm{Beta}(2, 2)$ に $N=10, k=7$ のデータを観測すると、形式的には事前と同じベータ族のまま、パラメータが $(2+7, 2+3) = (9, 5)$ に更新されるだけ。難しい積分は一切要りません。これが解析的に閉じる共役モデルの恩恵で、複雑なモデルに進む前の「教科書的に綺麗な特殊例」として位置づけられます。
予測分布の導出
事後分布が手に入ったので、次のコイン投げの結果 $\hat{x}$ を予測します。予測分布の定義式に代入し、
$$ p(\hat{x} \mid \bm{X}) = \int_0^1 p(\hat{x} \mid \mu)\, p(\mu \mid \bm{X})\, d\mu = \int_0^1 \mu^{\hat{x}}(1-\mu)^{1-\hat{x}}\, \mathrm{Beta}(\mu \mid \hat{a}, \hat{b})\, d\mu $$
$\hat{x} = 1$（表）の場合を計算します。
$$ p(\hat{x} = 1 \mid \bm{X}) = \int_0^1 \mu \cdot \frac{\Gamma(\hat{a}+\hat{b})}{\Gamma(\hat{a})\Gamma(\hat{b})}\, \mu^{\hat{a}-1}(1-\mu)^{\hat{b}-1}\, d\mu $$
これは $\mu$ に関する $\mathrm{Beta}(\hat{a}, \hat{b})$ のもとでの $\mu$ の期待値 $\mathbb{E}[\mu]$ にほかなりません。ベータ分布の期待値の公式 $\mathbb{E}[\mu] = \hat{a}/(\hat{a}+\hat{b})$ を使うと、
$$ \boxed{\; p(\hat{x} = 1 \mid \bm{X}) = \frac{\hat{a}}{\hat{a} + \hat{b}} = \frac{a + k}{a + b + N} \;} $$
具体例で確かめましょう。$N = 5$、$k = 3$、一様事前 $a = b = 1$ なら、
$$ p(\hat{x} = 1 \mid \bm{X}) = \frac{1 + 3}{1 + 1 + 5} = \frac{4}{7} \approx 0.571 $$
頻度論の点推定 $k/N = 0.6$ より少しだけ $0.5$ に引き寄せられた値が得られます。これはラプラスの継承則（rule of succession）として古典的に知られた結果です。5回投げて全部表（$k = 5$）でも、
$$ p(\hat{x} = 1 \mid \bm{X}) = \frac{1 + 5}{1 + 1 + 5} = \frac{6}{7} \approx 0.857 $$
頻度論なら $5/5 = 1.0$ ですが、ベイズなら「もしかしたら裏も出るかも」というニュアンスが残ります。少数データへの過剰な確信を防ぐ自然な仕組みです。
ここまでで「事後分布」と「予測分布」が手に入ったので、ベイズ推論の学習・予測サイクルが一周しました。次は、頻度論との接続点 — MAP推定と最尤推定 — を見ておきましょう。
MAP推定と最尤推定の対応
ベイズ事後分布から最頻値（mode）を取り出した推定量をMAP推定（Maximum A Posteriori）と呼びます。
$$ \hat{\theta}_{\text{MAP}} = \arg\max_\theta p(\theta \mid \mathcal{D}) = \arg\max_\theta \left[ \log p(\mathcal{D} \mid \theta) + \log p(\theta) \right] $$
一方、最尤推定（MLE）は事前分布を考えず尤度だけを最大化する手法です。
$$ \hat{\theta}_{\text{MLE}} = \arg\max_\theta \log p(\mathcal{D} \mid \theta) $$
つまりMAP = MLE + 事前分布による正則化項という関係です。ベルヌーイ-ベータの例で確認しましょう。事後分布 $\mathrm{Beta}(\hat{a}, \hat{b})$ の最頻値は（$\hat{a}, \hat{b} > 1$ なら）
$$ \hat{\mu}_{\text{MAP}} = \frac{\hat{a} – 1}{\hat{a} + \hat{b} – 2} = \frac{a + k – 1}{a + b + N – 2} $$
$a = b = 1$（一様事前）なら $\hat{\mu}_{\text{MAP}} = k/N$ となり、MLEと完全一致します。事前分布が無情報なら、MAPはMLEに退化するわけです。$a = b > 1$ にすれば「$\mu$ は $0.5$ あたりに引っ張られる」効果が入り、L2正則化や重み減衰と本質的に同じ役割を担います。
この対応関係はベイズ推論の解釈に欠かせません。よく聞く「リッジ回帰はガウス事前のMAP推定」「LASSO はラプラス事前のMAP推定」もこの図式の応用例で、機械学習の正則化はベイズ推論として書き直せるものが多いのです。
ただし、MAPはあくまで事後分布の代表点に過ぎず、不確実性の情報を捨ててしまいます。予測分布のように分布全体を使うほうが、ベイズの旨味は大きいのです。
同じデータ（$N = 5, k = 4$）に対して、事前分布だけを変えた2パターンを並べてみると、MAPとMLEの関係が手で触れるように見えてきます。左の一様事前 Beta(1, 1) では事前が水平線なので、事後分布の最頻値（MAP）は尤度の最大点（MLE = 0.8）と完全に一致します。右の対称事前 Beta(5, 5) では「$\mu = 0.5$ あたり」という前知識が効いて、事後のピークは MLE よりも 0.5 寄りに引き戻されます。これがまさに「事前分布による正則化」が機械学習の点推定に与える影響の縮図です。
ここまでは尤度と事前の組が共役であり、解析解がきれいに書けるという恵まれた状況でした。現実にはそうではないケースが大半です。次の節では共役を外れた世界を覗きます。
共役でない場合 — MCMCと変分推論
何が困るのか
例えば、ロジスティック回帰の重み $\bm{w}$ にガウス事前 $\mathcal{N}(0, \sigma^2 I)$ を置くと、尤度 $\prod_n \sigma(\bm{w}^\top \bm{x}_n)^{y_n}(1 – \sigma(\bm{w}^\top \bm{x}_n))^{1 – y_n}$ はシグモイドの積で、事前と掛け合わせても標準分布の形にはなりません。事後分布は解析的に書けず、規格化定数 $p(\mathcal{D}) = \int p(\mathcal{D} \mid \bm{w}) p(\bm{w}) d\bm{w}$ も高次元積分で評価不能です。
このような場合、二つの大きな処方箋があります。
MCMC（Markov Chain Monte Carlo）
事後分布から「サンプル」を取り出し、それを使って期待値や分位点を近似する戦略です。代表的なアルゴリズムにメトロポリス・ヘイスティングス、ギブスサンプリング、ハミルトニアンモンテカルロ（HMC）、NUTS（HMCの自動チューニング版）などがあります。
メトロポリス・ヘイスティングスの心臓は単純です。現在の点 $\theta$ から提案分布で $\theta’$ を生成し、受容率
$$ \alpha = \min\left(1, \frac{p(\theta’ \mid \mathcal{D})}{p(\theta \mid \mathcal{D})}\right) = \min\left(1, \frac{p(\mathcal{D} \mid \theta’)\, p(\theta’)}{p(\mathcal{D} \mid \theta)\, p(\theta)}\right) $$
で受容/棄却を決めます。比を取るのでエビデンス $p(\mathcal{D})$ は消え、規格化なしの事後分布だけ計算できれば回ります。これがMCMCを実用的にしている肝です。
変分推論（Variational Inference, VI）
事後分布 $p(\theta \mid \mathcal{D})$ を扱いやすい分布族 $q_\phi(\theta)$（例: ガウス）で近似し、両者のKLダイバージェンスを最小化します。
$$ \phi^* = \arg\min_\phi \mathrm{KL}(q_\phi(\theta) \,\|\, p(\theta \mid \mathcal{D})) $$
これはELBO（Evidence Lower Bound）
$$ \mathcal{L}(\phi) = \mathbb{E}_{q_\phi}[\log p(\mathcal{D}, \theta)] – \mathbb{E}_{q_\phi}[\log q_\phi(\theta)] $$
の最大化と等価で、勾配法で解けます。VAE（Variational Autoencoder）はこの枠組みを深層学習に拡張したものです。MCMCに比べて高速ですが、$q_\phi$ の表現力に近似精度が縛られます。
左パネルは「複雑な真の事後（二峰性）」を「扱いやすいガウス $q$」で近似する様子を示しています。$q$ は表現力に縛られて片方の山に張り付き、もう一方を取り逃がします。これがVIの典型的なバイアスです。右パネルでは横軸を変分パラメータ $\phi$ にとってELBOの形を可視化しました。エビデンス $\log p(\mathcal{D})$ は $\phi$ に依存しない上限で、ELBOがどこまで上がっても残る間隔がそのまま KL($q \Vert p$) です。最適化アルゴリズムが ELBO を持ち上げると、それに比例して KL が縮む — このシーソー関係こそ変分推論の心臓です。
どちらを使うか
| 観点 | MCMC | 変分推論 |
|---|---|---|
| 精度 | 理論的には正解に収束 | $q$ の表現力に依存 |
| 速度 | 低速 | 高速 |
| 高次元 | 苦手（特にメトロポリス） | 得意 |
| 主な用途 | 統計モデリング、小〜中規模 | 深層学習、大規模 |
実務では「PyMCやStanで小さなモデルをMCMC、TensorFlow ProbabilityやPyroで大規模VI」という棲み分けが定着しています。
ここまで「事後分布を求めるアルゴリズム」の議論をしましたが、ベイズ推論で得られる区間が、頻度論のいう信頼区間とどう違うのかも見落とせないポイントです。次節でそこを整理します。
信頼区間 vs 信用区間 — 似て非なる2つの区間
統計学者から教えてもらった「信頼区間（confidence interval）」と、ベイズの「信用区間（credible interval）」は、見た目こそ似ていても解釈がまったく違います。
頻度論の信頼区間
母比率 $\mu$ の 95% 信頼区間が $[0.45, 0.65]$ と出たとき、これが意味するのは
同じ実験を何度も繰り返してそのたびに信頼区間を計算すると、そのうち約 95% の区間が真の $\mu$ を含む
であって、「真の $\mu$ がこの区間に入る確率が 95%」ではない点に注意です。頻度論では $\mu$ は確率変数ではなく固定値なので、「$\mu$ が区間に入る確率」という言い方自体が不適切です。
ベイズの信用区間
ベイズの 95% 信用区間が $[0.45, 0.65]$ なら、その意味はまさしく
事後分布のもとで、$\mu$ がこの区間に入る確率は 95%
です。直感に最も近いのはこちらのほうで、しばしばニュースや科学コミュニケーションで「信頼区間」と呼ばれているもののイメージは、実はベイズの信用区間に近いことが多いのです。
どちらが正しいか？
「正しい」というより前提が違うだけです。ベイズは「パラメータも確率分布で表現する」立場、頻度論は「パラメータは固定値、データが確率変数」という立場。前提を選んだら、それに沿った解釈をしなければなりません。実務では「直感的に伝わりやすい」「事前情報を取り入れやすい」という理由でベイズ流の信用区間が好まれる場面が増えていますが、規制業界（医薬品など）ではいまだに頻度論的検定が支配的です。
左パネルは、同じ真値 $\mu = 0.5$（黒破線）から20回の実験を繰り返した頻度論的シミュレーションです。実験ごとにサンプルが揺れ、その都度95%信頼区間が引かれます。20本のうち1〜2本は真値を含まず（赤）、残り（青）が真値を捕まえます。「区間自体が確率変数で、$\mu$ は固定」という頻度論の立場が一望できます。右パネルはベイズ流で、データから1本だけ事後分布を作り、その中央95%領域を塗ったもの。$\mu$ そのものが確率変数として描かれ、「この区間内に $\mu$ が入っている確率は95%」と直接読めるのが信用区間の強みです。
理論を一巡したので、いよいよPythonで手を動かして確かめます。
Pythonで実装する — 解析解とMCMCの照合
環境準備
pip install numpy matplotlib scipy pymc arviz
PyMCはMCMCの定番ライブラリで、NUTSサンプラーを標準装備しています。arvizは結果可視化用です。
1. 解析解で事後分布と予測分布を可視化
ベルヌーイ-ベータ共役モデルの事後分布を、データ件数を増やしながらアニメーション風に描画して、信念が更新される様子を見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
# 真のコインの表確率
mu_true = 0.7
rng = np.random.default_rng(seed=42)
# 事前分布のハイパーパラメータ
a_prior, b_prior = 1.0, 1.0 # 一様事前
# データ件数のスナップショット
N_list = [0, 1, 5, 20, 100, 500]
fig, axes = plt.subplots(2, 3, figsize=(13, 7))
mu_grid = np.linspace(0, 1, 500)
# 累積データ
data = rng.binomial(1, mu_true, size=max(N_list))
for ax, N in zip(axes.flat, N_list):
k = int(data[:N].sum()) if N > 0 else 0
a_post = a_prior + k
b_post = b_prior + (N - k)
pdf = beta.pdf(mu_grid, a_post, b_post)
ax.plot(mu_grid, pdf, lw=2, color="steelblue")
ax.fill_between(mu_grid, pdf, alpha=0.25, color="steelblue")
ax.axvline(mu_true, color="crimson", linestyle="--", label=f"true μ = {mu_true}")
ax.set_title(f"N = {N}, k = {k}\nBeta({a_post:.0f}, {b_post:.0f})")
ax.set_xlabel("μ")
ax.set_ylabel("posterior density")
ax.legend(loc="upper left", fontsize=9)
plt.suptitle("Sequential Bayesian updating (Beta-Bernoulli)", y=1.02, fontsize=13)
plt.tight_layout()
plt.show()
上のサブプロット群から、データ件数の増加に伴って事後分布がどう変化するかが一目で読み取れます。$N = 0$（事前分布のみ）では一様分布ですが、$N = 1$ で早くも片側に傾き、$N = 20$ では真値 $\mu = 0.7$ 周辺に山がはっきり立ち上がります。$N = 500$ では分布の幅がぐっと狭まり、真値の周りに針のような事後分布が形成されます。データが増えるほど確信度（事後分布の集中度）が高まるというベイズの基本動作が視覚的に確認できます。
各パネルにはMAP（オレンジ点線）と95%信用区間（軸下の帯）も重ねてあります。$N$ が小さいうちは区間が極端に広く、$\mu$ が0.2から0.9まで「ありえる」と評価されますが、$N = 50$ あたりで区間幅が0.1程度まで収縮し、$N = 500$ になると針のような幅で真値を捕まえます。MAP値も $N$ が増えるにつれて真値0.7に近づいていきます。点推定（MAP）と区間推定（信用区間）が同じ事後分布から自然に取り出せるのは、分布として全情報を保持しているベイズ的アプローチの利点です。
次に、予測分布が表である確率の収束を見ます。
# 予測分布: p(x_hat = 1 | X) = a_post / (a_post + b_post)
N_range = np.arange(1, 501)
data_long = rng.binomial(1, mu_true, size=N_range[-1])
k_cum = np.cumsum(data_long)
pred_mean = (a_prior + k_cum) / (a_prior + b_prior + N_range)
# 事後分布の信用区間 (2.5% / 97.5%)
ci_low = beta.ppf(0.025, a_prior + k_cum, b_prior + (N_range - k_cum))
ci_high = beta.ppf(0.975, a_prior + k_cum, b_prior + (N_range - k_cum))
plt.figure(figsize=(10, 5))
plt.fill_between(N_range, ci_low, ci_high, alpha=0.25, color="steelblue", label="95% credible interval of μ")
plt.plot(N_range, pred_mean, color="navy", lw=2, label="posterior predictive p(x̂=1|X)")
plt.axhline(mu_true, color="crimson", linestyle="--", label=f"true μ = {mu_true}")
plt.xscale("log")
plt.xlabel("number of observations N")
plt.ylabel("probability")
plt.title("Convergence of posterior predictive and credible interval")
plt.legend(loc="lower right")
plt.grid(alpha=0.3)
plt.show()
ここから二つの大事な傾向が読み取れます。第一に、予測平均（紺線）が真値 $0.7$（赤破線）に収束していくこと。これは大数の法則のベイズ版と呼べる現象です。第二に、95% 信用区間（水色帯）が $N$ の増加とともに収縮していくこと。これはパラメータの不確実性が観測によって減っていく様子を直接表現しています。少数データでは区間が広く、ベイズが「自信のなさ」をちゃんと正直に出していることが分かります。
2. PyMCで同じ問題をMCMCで解いて比較する
共役なので解析解で十分ですが、検算のためMCMCで近似事後分布を作り、解析解と重ね描きしてみます。
import pymc as pm
import arviz as az
# データを固定して再現性を確保
mu_true = 0.7
rng = np.random.default_rng(seed=0)
N = 100
data_obs = rng.binomial(1, mu_true, size=N)
k = int(data_obs.sum())
# PyMCモデル
with pm.Model() as coin_model:
mu = pm.Beta("mu", alpha=1.0, beta=1.0) # 事前
x = pm.Bernoulli("x", p=mu, observed=data_obs) # 尤度
idata = pm.sample(draws=4000, tune=1000, chains=4,
random_seed=0, progressbar=False)
posterior_samples = idata.posterior["mu"].values.flatten()
# 解析解
a_post = 1.0 + k
b_post = 1.0 + (N - k)
mu_grid = np.linspace(0, 1, 500)
analytic = beta.pdf(mu_grid, a_post, b_post)
plt.figure(figsize=(9, 5))
plt.hist(posterior_samples, bins=80, density=True, alpha=0.45,
color="steelblue", label="PyMC MCMC samples")
plt.plot(mu_grid, analytic, color="crimson", lw=2,
label=f"Analytic Beta({a_post:.0f}, {b_post:.0f})")
plt.axvline(mu_true, color="black", linestyle="--", label="true μ = 0.7")
plt.xlabel("μ")
plt.ylabel("density")
plt.title(f"MCMC vs analytic posterior (N={N}, k={k})")
plt.legend()
plt.tight_layout()
plt.show()
print(az.summary(idata, var_names=["mu"], round_to=4))
ヒストグラム（青）と解析解の曲線（赤）はほぼ完全に重なります。MCMCサンプラーは事後分布から正しくサンプリングできていることが視覚的に確認でき、要約統計（az.summary）が出す平均・標準偏差・HDI（最高密度区間）も解析解と一致します。これは「MCMCが期待通り動いている」という基本テストでもあり、複雑なモデルに進む前の信頼性確認として常に意識しておくとよいでしょう。
左パネルでは、メトロポリス・ヘイスティングス法で得た16,000サンプルのヒストグラムと、解析解 $\mathrm{Beta}(\hat a, \hat b)$ の曲線がほぼぴったり重なる様子が確認できます。右パネルはMCMCのトレースプロットで、サンプル列が解析平均（赤破線）の周りを定常的に振動しており、収束していることが視覚的に判定できます。受容率が0.2〜0.5あたりに収まっているのもメトロポリス・ヘイスティングスの健全な兆候です。PyMCはこれをより洗練されたNUTSサンプラーで自動化していますが、根本のアイデアは同じです。
PyMCの強みは、ここから尤度や事前を差し替えるだけで非共役なモデルにも即対応できる点です。例えばベルヌーイをポアソンに、ベータをガンマに変える、階層構造を組む、回帰係数に正則化的な事前を入れる、といった拡張が pm.Model ブロック内の数行で実現します。
3. 異なる事前分布の影響を比較する
事前の置き方で事後がどう変わるかも見ておきます。
priors = [(1, 1, "uniform: a=b=1"),
(5, 5, "weakly biased: a=b=5"),
(50, 50, "strongly biased: a=b=50"),
(1, 9, "asymmetric: a=1, b=9")]
mu_true = 0.7
rng = np.random.default_rng(seed=1)
N = 20
data_obs = rng.binomial(1, mu_true, size=N)
k = int(data_obs.sum())
fig, ax = plt.subplots(figsize=(10, 5))
mu_grid = np.linspace(0, 1, 500)
for a0, b0, label in priors:
a_post = a0 + k
b_post = b0 + (N - k)
ax.plot(mu_grid, beta.pdf(mu_grid, a_post, b_post), lw=2, label=label)
ax.axvline(mu_true, color="black", linestyle="--", label="true μ = 0.7")
ax.set_xlabel("μ")
ax.set_ylabel("posterior density")
ax.set_title(f"Effect of prior on posterior (N={N}, k={k})")
ax.legend()
plt.tight_layout()
plt.show()
この比較から、事前分布が事後分布に与える影響の構造が読み取れます。一様事前（青）は素直にデータに従い、$\mu = 0.7$ 付近にピークを形成します。一方、強い事前（$a = b = 50$、緑など）はデータの $20$ サンプルでは簡単には動かず、ピークは依然として $0.5$ 寄りに留まります。非対称な事前（$a = 1, b = 9$、紫）は「裏が出やすい」と信じ込んでいるため、$20$ サンプルでは事後ピークが真値より低い側に引きずられます。事前の強さが「データに対するすり合わせのしにくさ」を支配するわけで、データ量とのバランスを意識して事前を設計する必要があります。
ただし、$N \to \infty$ ではどの事前から始めても事後分布は同じ点に収束します（これをベイズの一致性と呼びます）。事前の影響は「データが少ない領域での補正項」と考えるとイメージしやすいでしょう。
左パネルは4種類の事前分布そのものの形状です。一様 Beta(1,1) は無情報、Jeffreys Beta(0.5,0.5) は両端で発散するU字、強い対称 Beta(10,10) は $\mu=0.5$ に山、非対称 Beta(2,8) は「裏が出やすい」という事前知識を表します。右パネルは同じデータ ($N = 20, k = 14$) を観測した後の事後分布で、一様事前はデータに素直に従って $k/N = 0.7$ 付近にピークを作る一方、強い対称事前は依然として $0.5$ 寄りに留まり、非対称事前は真値より低い側に引きずられます。データ20件では事前の強さが事後の位置を左右することが、視覚的にはっきり確認できます。
最後に、ここまで学んだベイズ推論の枠組みが、どこへ拡張されていくのか展望を述べておきます。
実用例 — ベイズA/Bテスト
最後に、ベイズ推論が実務で最も多用される場面のひとつ — A/Bテスト — に触れておきます。Webサイトの旧デザイン（A）と新デザイン（B）でクリック率を比較したいとしましょう。頻度論なら $p$ 値や信頼区間で「有意差あり/なし」を判定しますが、ベイズ流の問いはもっと直接的です。「新デザインが旧デザインより優れている確率はいくつか？」
モデル設定
- A: $n_A$ 人に提示、$k_A$ 人クリック。クリック率 $\mu_A \sim \mathrm{Beta}(1, 1)$、観測は $k_A \sim \mathrm{Binomial}(n_A, \mu_A)$。
- B: 同様に $n_B, k_B$、$\mu_B \sim \mathrm{Beta}(1, 1)$。
事後分布は独立にベータで、
$$ \mu_A \mid \mathcal{D} \sim \mathrm{Beta}(1 + k_A, 1 + n_A – k_A), \quad \mu_B \mid \mathcal{D} \sim \mathrm{Beta}(1 + k_B, 1 + n_B – k_B) $$
知りたいのは
$$ P(\mu_B > \mu_A \mid \mathcal{D}) = \int \int \mathbb{1}[\mu_B > \mu_A]\, p(\mu_A \mid \mathcal{D})\, p(\mu_B \mid \mathcal{D})\, d\mu_A d\mu_B $$
この積分はモンテカルロで簡単に評価できます。事後分布から大量にサンプルを取って、$\mu_B > \mu_A$ となった割合を数えるだけ。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta as beta_dist
# 観測データ
n_A, k_A = 1000, 120 # 旧デザイン: 1000人中120クリック
n_B, k_B = 1000, 145 # 新デザイン: 1000人中145クリック
# 事後分布
a_A, b_A = 1 + k_A, 1 + n_A - k_A
a_B, b_B = 1 + k_B, 1 + n_B - k_B
# モンテカルロサンプリング
N_mc = 200_000
rng = np.random.default_rng(seed=42)
samples_A = rng.beta(a_A, b_A, size=N_mc)
samples_B = rng.beta(a_B, b_B, size=N_mc)
# 「Bが優れている確率」と「リフト」
prob_B_better = np.mean(samples_B > samples_A)
lift = (samples_B - samples_A) / samples_A
lift_mean = lift.mean()
lift_hdi = np.percentile(lift, [2.5, 97.5])
print(f"P(μ_B > μ_A | data) = {prob_B_better:.4f}")
print(f"相対リフト 期待値 = {lift_mean*100:.2f}%, 95% CI = [{lift_hdi[0]*100:.2f}%, {lift_hdi[1]*100:.2f}%]")
# 可視化: 事後分布の重ね描き
mu_grid = np.linspace(0.08, 0.20, 500)
plt.figure(figsize=(10, 5))
plt.plot(mu_grid, beta_dist.pdf(mu_grid, a_A, b_A), lw=2, label=f"A: Beta({a_A}, {b_A})")
plt.plot(mu_grid, beta_dist.pdf(mu_grid, a_B, b_B), lw=2, label=f"B: Beta({a_B}, {b_B})")
plt.fill_between(mu_grid, beta_dist.pdf(mu_grid, a_A, b_A), alpha=0.25)
plt.fill_between(mu_grid, beta_dist.pdf(mu_grid, a_B, b_B), alpha=0.25)
plt.title(f"Posterior of click-through rates (P(B > A) = {prob_B_better:.3f})")
plt.xlabel("click-through rate μ")
plt.ylabel("posterior density")
plt.legend()
plt.tight_layout()
plt.show()
出力結果から、Bの事後分布がAより右にずれていることが視覚的に分かります。P(μ_B > μ_A | data) ≈ 0.94 程度の数値が得られ、「Bのほうが優れている確率は約94%」とそのまま意思決定者に伝えられます。さらに相対リフトの期待値と信用区間も同時に出るため、「期待リフト約20%、95%の確信で5%〜37%」のように利益見積もりまで一気通貫で語れるのがベイズA/Bテストの強みです。
頻度論的なA/Bテストとの最大の違いは、サンプルサイズを事前に決めなくてよい点と、結果がそのまま意思決定の言葉になる点です。途中で結果を覗いて止めても多重比較の問題が起きにくく（停止規則を事前に決めれば）、ビジネスサイドにも説明しやすい。これがネット業界で「ベイズA/Bテスト」が広く採用されている理由です。
拡張への入口 — 階層ベイズ、ベイズ最適化、確率的プログラミング
階層ベイズ
「ハイパーパラメータ自身にも事前を置く」という階層モデリングは、データのグループ構造を扱う強力な道具です。例えば、複数のWebサイトでクリック率 $\mu_i$ を推定したいとき、$\mu_i \sim \mathrm{Beta}(\alpha, \beta)$ とおき、さらに $\alpha, \beta$ にもハイパー事前を置きます。すると「サイト間で情報を共有しつつ個別の異質性も保つ」部分プーリング（partial pooling）が自然に実現します。完全プーリング（全サイトを束ねる）と独立推定（サイトごとに別々）の中間的な振る舞いを示し、データの少ないサイトでも安定した推定が得られるのです。
3層構造のグラフィカルモデルで階層ベイズを描くと、情報の流れが明確になります。最上層のハイパー事前 $p(\alpha, \beta)$ から中層のハイパーパラメータ $\alpha, \beta$ が決まり、それを共通の親としてサイトごとの $\mu_i$ が生成され、最終的に各サイトで $N_i$ 件の観測が得られる。このグラフを伝って情報が伝播するため、観測数の少ないサイト（例：サイト3, 5）は他サイトの情報を間接的に借り、推定が安定します。逆に観測数の多いサイト（サイト1）は自分のデータが優勢で、独自の値を保ちます。これが「部分プーリング」と呼ばれる現代ベイズ統計の代表的な構造です。
ベイズ最適化
「評価コストの高い関数の最大値を、最少回数の評価で見つけたい」というニーズに応えるのがベイズ最適化です。関数にガウス過程の事前を置き、観測点から事後分布を更新、次の評価点を獲得関数（expected improvement など）で決めるサイクルを回します。ハイパーパラメータチューニング、材料探索、化学反応の条件最適化など、現代の科学的探索の標準ツールになっています。
確率的プログラミング
PyMC、Stan、Pyro、NumPyro といった確率的プログラミング言語（PPL）を使うと、複雑な階層モデルや確率モデルを数十行で記述し、推論はライブラリが自動で行ってくれます。MCMCも変分推論も同じインターフェースから呼べることが多く、モデルの試作が劇的に楽になりました。ベイズ推論を実務に持ち込むハードルは、ここ10年で大幅に下がったと言えます。
まとめ
本記事では、ベイズ推論の枠組みを一気通貫で見てきました。
- ベイズの定理: $p(\theta \mid \mathcal{D}) \propto p(\mathcal{D} \mid \theta)\, p(\theta)$。事前 × 尤度 = 事後の比例関係が、ベイズ推論の心臓です。
- 共役性: ベルヌーイ-ベータの例で見たように、共役な事前を選ぶと事後が同じ族に留まり、解析的に閉じます。事前のパラメータは「擬似カウント」と解釈できます。
- 予測分布: 学んだ事後分布で $\theta$ を積分消去することで、パラメータ不確実性を予測へ伝播できます。
- MAP推定: 事後分布の最頻値はMLEに事前項を足したもので、L2正則化など機械学習の正則化と対応します。
- MCMCと変分推論: 共役でないモデルの事後分布を、サンプリング/最適化で扱う2大処方箋です。
- 信用区間 vs 信頼区間: 直感に近いのはベイズの信用区間。ただし前提が違うだけで「どちらが正しい」議論ではありません。
- 拡張: 階層ベイズ、ベイズ最適化、確率的プログラミングが、現代ベイズ推論の主戦場です。
ベイズ推論の真の魅力は、不確実性を「定数」ではなく「分布」として扱う発想にあります。データが少ないときの控えめな予測、事前知識の組み込み、A/Bテストでの確率的意思決定 — どれもベイズの自然な帰結です。
次のステップとして、以下の記事も参考にしてください。