ベイズ統計：MCMCは「歩き方」と「数え方」の合成語だった【第33回】

はじめに

第31章はベイズ法です。前回のモデル選択で「予測が目的か説明が目的か」という区別を扱いましたが、今回はもっと手前の、確率という言葉の使い方そのものが変わる回になります。

正直に書くと、私はベイズ統計を何度も学ぼうとして、何度も同じ場所で止まっていました。ベイズの定理は分かる。事後分布を計算する手順も追える。それでも「結局ベイズ統計って何なのか」が最後まで言われないまま終わる、という感覚がずっとありました。

今回それが解けたのは、計算の話を後回しにして、絵を2枚並べたからです。頻度論とベイズで「何が固定されていて何が動くのか」を描いただけで、話が通りました。難しい数学は1つも出てきません。

そしてこの記事の中盤で、私は一度完全に詰まりました。 MCMC（マルコフ連鎖モンテカルロ法）です。メトロポリス法のアルゴリズムを順番に追うことはできたのに、それが何をしている道具なのか分からない。原因は名前を分解していなかったことでした。MCMC は独立に発明された2つの技術がくっついた言葉で、片方だけ先に理解すれば済む話でした。その回り道もそのまま記録します。

いちばん驚いたのは終盤です。t分布が正規分布の混合分布だったという事実に行き当たりました。第7回で「分散を知らない罰金」として学んだあの分布が、今回の階層ベイズとまったく同じ構造をしていました。第5回の負の二項分布も同じでした。「消すべきパラメータを最尤推定値のような1点に決め打つと、たいてい自信過剰になる」という筋が、連載で別々に学んだ回をまとめて回収してくれました。

いつものように、出てくる数値はすべて自分で計算し、答えが分かっている問題と突き合わせています。MCMC は自分で実装し、2次元グリッドによる数値積分を「答え」として用意してからその結果を検証しました。

この回で扱う用語

ベイズが変えたのは「何が動くか」だけだった

困りごとから始める

第11回の区間推定で、私はひとつ気持ち悪さを抱えたまま先に進んでいました。

95%信頼区間を計算して [0.14, 0.52] という答えが出たとき、「クリック率が 0.14 から 0.52 の間にある確率は 95%」と言ってはいけないと学びました。95% という数字は区間についているのではなく、区間の作り方についている成績だったからです。真のクリック率 $\theta$ は定数なので、確率を持てない。

理屈は分かりました。でも、じゃあ目の前のこの1本の区間について何が言えるのかが分からない。実務で「で、結局クリック率はどのくらいなんですか」と聞かれたときに、「この区間は、同じ手続きを無限回繰り返せば95%の割合で真の値を含むような区間です」と答えるのは、明らかに何かがおかしい。

ベイズ統計はこの気持ち悪さに正面から答えます。

図で見る：同じ問題を見る2つの絵

同じ問題（広告のクリック率）を扱っているのに、動いているものが逆です。

左（頻度論） では、真のクリック率 $\theta = 0.30$ は1点に固定された定数です。神様だけが知っている値で、確率変数ではありません。動くのは標本のほうで、20回試すたびに 7/20、9/20、5/20 と推定値がばらつきます。

右（ベイズ） では、観測した 6/20 はもう見てしまったので固定です。動くのは $\theta$。「今のところ 0.3 あたりが一番ありそうだが、0.2 や 0.45 でもおかしくない」という自分の知識の状態を分布で書く。

これが唯一の本質的な差です。そして $\theta$ が確率変数になった瞬間に、さっき禁じられていた文が書けるようになります。

$\theta$ が 0.146 と 0.522 の間にある確率は 95%

よくある引っかかり：クリック率は1つの値でしょう

ここで多くの人が引っかかります（私も引っかかりました）。「現実のクリック率は1つの値として存在するはずで、それが分布するのはおかしいのでは」という疑問です。

これは正しい疑問で、答えはベイズの分布は $\theta$ の物理的なばらつきではないということです。表しているのは自分の知識の不確かさです。

サイコロの例で言うと、振ってしまって伏せた手のひらの下にある目は、物理的にはもう1つに決まっています。それでも「6である確率は 1/6」と言えます。確率が対象の性質ではなく、自分の情報状態を表しているからです。ベイズの $\theta$ の分布はこれと同じ種類のものです。

この立場の違いには名前があって、確率を「長期的な頻度」と読むのが頻度論、「知識の度合い」と読むのがベイズです。

事後分布は「かけ算して面積を1にした」だけ

計算はこれで全部です。

分母は $\theta$ を含まないただの定数なので、実質はかけ算だけです。だから実務では比例だけ書きます。

①データを見る前の考えに、②実際に見たデータの尤度をかけて、③面積が1になるよう割る。

第10回でやった最尤法は、②のピークを探すだけの作業でした。ベイズは②を捨てずに全体を持ち歩きます。 この差が後で効いてきます。

3つの代表値が出てくる

事前・尤度・事後を並べた図の右端で気づいてほしいのは、ピークは 0.300 なのに平均は 0.318 とズレていることです。分布が右に裾を引いているためです。

MAP が最尤推定と一致するのは重要な接続点です。事前分布を平らにしたベイズの MAP は最尤推定と一致するので、最尤法はベイズの特殊ケースとして読めます。

ただし一致するのは MAP だけです。事後平均 0.3182 と事後中央値 0.3126 は最尤推定 0.3000 と一致しません。さらに「平ら」という性質はパラメータの取り方に依存する（後述）ので、$\psi = \text{logit}(\theta)$ について平らな事前分布を置くと MAP も 0.30 にはなりません。

そして事後平均が $7/22 = (x+1)/(n+2)$ になるのも偶然ではありません。「分子に1、分母に2を足す」というあの謎の操作（機械学習のラプラススムージング）の出どころがこれです。$n=0$ のときに 1/2 を返すので、「何も見ていないなら半々」という自然な振る舞いになります。

なお似た形の Agresti-Coull 法は $z^2 \approx 4$ を使って成功2・失敗2を足す $(x+2)/(n+4)$ で、こちらは $a=b=2$ の事前分布に相当します。$+1/+2$ とは別物なので混同しないよう注意が必要です（自分は一度混同しました）。

共役とは何か：形が掛け算で壊れないこと

ここで用語をひとつ固めます。私はこの回の準備で「共役事前分布は予習済み」と思っていたのですが、記事を書く段階で「二項分布と共役ってどういう意味だっけ」と自分で分からなくなりました。 定義から確認します。

尤度の族を決めたとき、事後分布が事前分布と同じ族に戻ってくるような事前分布の族を、その尤度に対する共役事前分布という。

用語の注意が1つあります。 「二項分布と共役」は言葉を省略した言い方で、正確には「ベータ分布は二項尤度に対する共役事前分布である」。共役は「分布と分布」の関係ではなく、「尤度の族」と「事前分布の族」の関係です。

「ベータと二項が似ている」という話ではありません。掛け算しても形が壊れない相性の話です。

なぜベータ×二項がベータに戻るのか

式を3行並べるだけで分かります。正規化定数は $\theta$ を含まないので無視します。

掛けると指数が足されるだけです。

最後の形は $\text{Beta}(x+a,\; n-x+b)$ そのものです。

核心はこれです。 事前分布と尤度がどちらも「$\theta$ の何乗 × $(1-\theta)$ の何乗」という同じ骨格をしている。だから掛け算すると指数が足されるだけで、骨格が変わりません。

共役の正体は「関数の形が掛け算で保たれること」であって、それ以上の意味はない。

共役でないと本当に戻らないのか

言葉だけだと納得しづらいので、共役でない事前分布と並べました。ロジット正規分布（$\text{logit}(\theta)$ が正規分布に従う）を使います。ベータとよく似た形ですが、二項尤度に共役ではありません。

上段（ベータ） は、掛けた結果がぴったりベータ分布に戻ります。真の事後分布（緑）と $\text{Beta}(3.49, 24.68)$（赤破線）が重なって区別できません。密度の最大差は $4.0 \times 10^{-10}$ で、数値積分の誤差だけです。

下段（ロジット正規） は、②の尤度が上段とまったく同じなのに、③では同じ族で最善の当てはめをしても最大差 0.57 のズレが残ります（ピークの高さが合っていません）。これが「族に戻らない」ということです。

共役な組み合わせの一覧

試験で計算問題として出るのはこの表です。

この表の $\text{Gamma}(a,b)$ は $b$ を率として書いています。後で出てくる逆ガンマは $b$ を尺度として使うので、教科書ごとの流儀の違いに注意してください。

共通の型は1つだけです。

事前分布のパラメータに、データの十分統計量を足す。

理由は第9回の指数型分布族です。指数型分布族の尤度は「$\exp($ 自然パラメータ $\cdot$ 十分統計量 $-$ 正規化項 $)$」の形をしているので、同じ形の事前分布を掛けると指数の中で足し算になる。だから族が保たれます。

ここで注意点が2つあります。私は最初、「共役事前分布が存在するのは指数型分布族のときだけ」と書いたのですが、「だけ」という限定も、そこから引いた理由づけも誤りでした。

まず逆は成り立ちません。 $X \sim U(0, \theta)$ は台が母数に依存するので指数型分布族ではありませんが、パレート分布が共役です（尤度 $\theta^{-n}\mathbf{1}\{\theta \geq x_{\max}\}$ にパレートを掛けるとパレートに戻る）。

そしてロジスティック回帰の尤度は、実は $\beta$ について指数型分布族です。 対数尤度を書き直すと

で、「自然パラメータ $\beta$ ・十分統計量 $\sum y_i x_i$」の形をしています。実際 $\pi(\beta) \propto \exp(\beta^\top\tau - m A(\beta))$ という共役族（Diaconis-Ylvisaker 事前分布）は存在し、更新も閉じます。

後でロジスティック回帰で共役が使えない本当の理由は、その族の正規化定数が閉じた形で書けず、名前のついた分布にならないことです。「共役族が存在しない」のではなく「あっても手で扱えない」が正確でした。

正しい整理はこうです。尤度が指数型分布族なら、自然共役族を必ず構成できる（超パラメータは十分統計量の次元に「見かけの標本サイズ」を1つ足した分）。ただし逆は成り立たず、また構成できても実用的な形になるとは限りません。

信頼区間と信用区間：解釈は違うが数字はほぼ同じ

第11回から引きずっていた気持ち悪さに戻ります。同じデータ（$n=20$ で6クリック）で両方計算して並べました。

Wilson の信頼区間と信用区間が、小数第2位までほぼ一致します。 哲学がまるで違うのに、出てくる数字はほぼ同じです。

なぜ一致するのか

偶然ではありません。事前・尤度・事後を並べた図で見たとおり、事前分布が平らなら事後分布の形は尤度そのものです。そして尤度のピーク周りの形は、頻度論が正規近似で使っているものと同じです。同じ関数を、一方は「$\theta$ の確率分布」と読み、他方は「推定値のばらつき」と読んでいるという状況になります。

$n$ を大きくすると一致はさらに良くなります。ここから実務上の含意が出ます。

通常の条件下でデータが十分あるとき、ベイズにするかどうかで結論の数字はほぼ変わらない。 変わるのは「その数字を何と言い表せるか」だけ。

「通常の条件下で」と付けたのは、これがベルンシュタイン=フォン・ミーゼスの定理という漸近論に頼っているからです。真の値がパラメータ空間の内部にあり、事前分布がその近傍で正の密度を持ち、次元が固定されている必要があります。この記事の後半で見る $x=0$ のケースや不適切事前分布では実際に崩れます。

裏を返すと、$n$ が小さいとき・事前情報が本当にあるときにだけ、ベイズは違う答えを出します。

では信用区間の何が嬉しいのか

解釈です。信用区間なら

$\theta$ が 0.146 と 0.522 の間にある確率は 95%

でこの文が終わります。区間の作り方に言及する必要がありません。

もっとはっきりした差が、任意の質問に答えられることです。ベイズなら

と、聞かれた質問に対して数字が1つ返ります。 頻度論では「$\theta < 0.2$ である確率」という文自体が定義できません（$\theta$ は定数なので、この確率は 0 か 1 のどちらかで、どちらかは分からない、という答えになってしまう）。

用語の整理

日本語が1文字違いなので、試験では英語で覚えたほうが安全です。

そして信用区間には作り方が2種類あります。

単峰かつ対称な事後分布なら両者は一致します（対称でも二峰なら HPD は2つの区間に分かれるので一致しません）。歪んでいるとズレて、HPD のほうが狭くなります。

頻度論の基準で採点したらどうなるか

「ベイズの信用区間は頻度論の基準で見ると成績が悪いのでは」という疑問が湧きます。確かめました。横軸に真の $\theta$ を全部並べ、その区間が $\theta$ を含む確率を厳密に計算します。

$\theta = 0.05$ の行が要点です。教科書に最初に載っている Wald は、95% と名乗りながら実際は 63.89% しか当てていません。 3回に1回以上外している。一方 Clopper-Pearson は 98.41% で、外さないけれど広すぎます。

信用区間は頻度論の被覆率を目標に作っていないのに、頻度論の基準で採点しても Wald よりずっとまともです。 ギザギザしているのは二項分布が離散だからで、これは頻度論側の区間にも同じように出ます。

0クリックだったときが象徴的

$n=20$ で1件もクリックされなかった場合を計算しました。

Wald が潰れる理由は式を見れば明らかで、$\hat p = 0$ なので $\sqrt{\hat p(1-\hat p)/n} = 0$ になるからです。

このとき事後平均は $1/22 = 0.0455$、そして

と直接言えます。「クリック率が10%未満である確率は89%」という、そのまま意思決定に使える形です。

なお HPD 区間で計算すると、事後分布 $\text{Beta}(1,21)$ は単調減少なのでピークが端（$\theta = 0$）にあり、下端が 0 になります。上端は $1 - 0.05^{1/21} = 0.1329$ で、等裾の上端 0.1611 より狭い。片側だけの区間になるという、等裾との違いがはっきり出る例です。

無情報事前分布は本当に「無情報」なのか

事前分布を置くと決めたとき、次に来る疑問はこれです。「何も情報がないときは一様分布を置けばいいのでは」。

私はこれに引っかかりました。一様分布を置くのも、1つの積極的な主張に見えるからです。結論から言うと、その直感が正しいです。

一様分布は尺度を変えると偏った主張になる

$\theta$ に一様分布を置いた後で、同じものを別の尺度で見てみます。

• オッズ $\phi = \theta/(1-\theta)$ で見ると、小さい側に強く偏っています
• 対数オッズ $\psi = \log\bigl(\theta/(1-\theta)\bigr)$ で見ると、0 の近くに集まった山型（ロジスティック分布）になります

数値でも確認しました。$\theta \sim$ 一様 から生成した $\psi$ は、標準偏差 1.81402（理論値 $\pi/\sqrt{3} = 1.81380$）。そして

$\psi$ は全実数を動くので「$\psi$ について一様」は不適切事前分布になってしまいますが、図の範囲 $(-8, 8)$ に切って比べると一様なら $2/16 = 0.125$。実際は3.7倍になっています。

つまり「$\theta$ について無情報」は「$\psi$ について無情報」ではありません。 そしてこれは机上の話ではなく、ロジスティック回帰の係数は $\psi$ の側です。第18回で扱った $\log(\text{オッズ比})$ がまさにこれで、「$\theta$ に一様を置く」と決めた瞬間に、係数については「0 に近い」という主張をしていることになります。

第4回の変数変換（ヤコビアン）が効いているだけの話です。密度は変換で形が変わるので、「平ら」という性質は尺度に依存するのです。

ジェフリーズ事前分布：どの尺度でも同じ答えになるように作る

この問題を解決するために提案されたのがジェフリーズ事前分布です。

二項分布なら $I(\theta) = n/\bigl(\theta(1-\theta)\bigr)$ なので

なぜ $\sqrt{I(\theta)}$ を掛けると解決するのか。 変換不変性を数値で確認しました。

比が全て 1.000000。 「$\theta$ で無情報事前分布を作ってから変換する」と「最初から $\psi$ で無情報事前分布を作る」が、比例定数まで完全に一致します。どの尺度で考えても同じ結論になる――これがジェフリーズ事前分布の設計目標です。一様分布ではこれが成り立ちません。

直感的な読み方はこうです。フィッシャー情報量が大きい領域＝データが敏感に効く領域なので、そこに事前分布の重みを厚く置くと、尺度の取り方に左右されなくなります。

なおパラメータが複数あるときは $\pi(\boldsymbol\theta) \propto \sqrt{\det I(\boldsymbol\theta)}$ になりますが、多母数のジェフリーズ事前分布は必ずしも良い無情報事前分布になりません（そのため参照事前分布などの改良版が提案されています）。この記事では1母数の場合だけを扱います。

事前分布の形と、それが結論をどれだけ動かすか

ジェフリーズ $\text{Beta}(0.5, 0.5)$ は両端で無限に高くなる U 字型です。「0 や 1 に近い可能性を軽視しない」形になっています。

一様とジェフリーズはほぼ同じで、強い事前分布だけが大きくズレます。

$n$ を増やせば事前分布の影響は消えるか

消えます。率を 30% に固定したまま $n$ を増やしました。

$n$ が数百を超えると、事前分布を何にしても結論は実質同じになります。 事前分布の選択に神経を使う必要があるのは小標本のときだけです。

逆に言えば、小標本では事前分布が結論を左右するので、選んだ理由を説明できなければならない。 ここが「ベイズは恣意的だ」という批判の当たっている部分でもあります。

罠：不適切事前分布

上の表のハルデーン事前分布 $\text{Beta}(0,0)$ には注意が必要です。これは積分が発散するので、そもそも確率分布ではありません。こういうものを不適切事前分布（improper prior） と呼びます。

不適切でも事後分布がまともになるなら実用上は使えるのですが、壊れる場合があります。$x=0$ のデータで計算すると、

$a=0$ では事後分布 $\text{Beta}(0, 21)$ 自体が正規化できず、事後分布そのものが不適切になります（数値的には $\theta = 0$ に退化します）。「無情報を極限まで追求する」と使えなくなる、という教訓です。ジェフリーズが $\text{Beta}(0.5, 0.5)$ という「ほどよく弱い」ところに落ち着いているのは、この意味でも合理的です。

ここで一度詰まった：MCMC が何なのか分からない

ここまでは共役事前分布のおかげで、事後分布が式で書けました。$\text{Beta}(a,b) \times$ 二項 $\to \text{Beta}(x+a, n-x+b)$ で、足すだけでした。

問題は共役でないときです。 そこで MCMC という道具が出てきます。

正直に書くと、私はここで完全に止まりました。教科書のメトロポリス法のアルゴリズムを読んで、手順を追うことはできました。「提案して、比を計算して、確率で受容する」。動きも分かる。でもそれが何をしている道具なのか分からない。

原因は後から分かりました。名前を分解していなかったのです。

MCMC は独立に発明された2つの技術がくっついた名前で、それぞれ別々に理解できます。そして先に理解すべきなのは後ろの「モンテカルロ法」のほうでした。こちらにはマルコフ連鎖が一切出てきません。

順番に積み上げます。

Step 1：モンテカルロ法は「計算する代わりに数える」

まず円周率です。統計は一切出てきません。

正方形にダーツをランダムに投げて、円に入った割合を数える。面積比が $\pi/4$ なので、$4 \times (\text{円の中}/\text{全部})$ が $\pi$ になります。1億投で 3.141591（真値との差 $-0.000002$）まで出ました。

積分も公式も使っていません。投げて数えただけです。 これがモンテカルロ法の全部で、それ以上の中身はありません。名前の由来もモナコの Monte Carlo（カジノの街）で、要はサイコロを振って答えを出す方法の総称です。

発想の転換はここです。

難しい「計算」を、たくさんの「数え上げ」に置き換える。

Step 2：統計に持ち込む

いま仮に、事後分布 $\text{Beta}(7,15)$ から出た点が2万個手元にあるとします。

• 事後平均を知りたい → 足して割る
• 95%区間を知りたい → 並べて 2.5% と 97.5% の位置を見る
• $P(\theta < 0.2)$ を知りたい → 0.2 未満を数える

全部「数える」で済み、積分は一度も出てきません。個数を増やすと厳密解に収束します。

さらに強いのが変換してから数えられることです。

「オッズの分布が欲しい」「2群の差が欲しい」と言われても、サンプルを変換して数え直すだけです。積分をやり直す必要がありません。

ちなみに $E[1/\theta] = 3.4987$ と $1/E[\theta] = 3.1404$ が別物なのは第4回のイェンセンの不等式です。サンプルで扱うとこの区別が自動的に正しく処理されます。

Step 3：でも、そのサンプルはどうやって作るのか

ここが本題でした。 Step 2 では「事後分布からサンプルが出せる」ことを前提にしていました。

$\text{Beta}(7,15)$ なら有名な分布なのでライブラリで一発です。しかしベイズで出てくる事後分布は「事前 × 尤度」というその場で作られた変な形の関数で、名前も付いていません。

素朴な方法が棄却法です。分布を囲む箱にダーツを投げ、曲線の下に入ったものだけ採用する（円周率と同じ発想）。

1次元なら採用率 24.3%（図に描いた1500投の試行）で十分実用的です。理論値は 24.8% で、下の表の $d=1$ の行がそれに当たります。ところが次元が上がると壊れます。

理由は棄却法の図の右パネルです。高次元の立方体はほとんど「隅」でできています。 $d=20$ では、立方体に投げたダーツが内接球に入る確率は $2.5 \times 10^{-8}$ しかありません。これが次元の呪い（curse of dimensionality） です。

参考に、真面目にグリッドで数値積分する道も潰れています。各軸100点で刻むと、

階層モデルではパラメータが数十から数千個になります。この道は原理的に閉じています。

棄却法の敗因は「毎回ゼロから独立に候補を出す」ことです。 良い場所を見つけても、その情報を捨てて、また当てもなく投げる。

ならば「今いる良い場所の隣」を次の候補にすればいい。

これがマルコフ連鎖の出番です。

Step 4：マルコフ連鎖の滞在時間は定常分布になる

第22回の内容ですが、必要な部分だけ確認します。

マルコフ連鎖とは「次の状態が今の状態だけで決まる」確率的な移動です。3状態 A・B・C を推移確率に従って渡り歩かせると、歩数を増やすにつれ「各状態にいた時間の割合」が理論値に収束します。

一致します。この一致が MCMC の原理そのものです。

だから話を逆から組み立てられます。

定常分布が「事後分布」になるような歩き方を設計すれば、その足跡が事後分布からのサンプルになる。

全体の組み立て

ここまでを1本に繋げるとこうなります。これが最初から分かっていれば詰まらなかった部分です。

1. 事後平均や区間が知りたい → 積分が必要 → 高次元では不可能
2. でもサンプルがあれば数えるだけで済む（＝モンテカルロ法）
3. じゃあサンプルをどう作るか → 棄却法は高次元で死ぬ
4. 今いる場所の隣を歩くようにすれば高次元でも動ける（＝マルコフ連鎖）
5. 足跡の分布は定常分布になるので、定常分布＝事後分布になるよう歩き方を設計する
6. その足跡を数える（2に戻る）

最後の代償は実測しました。同じ「20000個」でも、独立サンプルなら事後平均のばらつきは 0.000729、MCMC では 0.001395（1.91倍）。逆算すると MCMC の20000個は独立換算で約5458個分です。

つまり MCMC は「独立性を犠牲にして、高次元でもサンプルが引ける」という取引をしています。

メトロポリス法：比を取ると分母が消える

歩き方の設計に入ります。アルゴリズムはこれで全部です。

def log_target(th):
    if th <= 0 or th >= 1: return -inf
    return 6*log(th) + 14*log(1-th)      # 事後 ∝ θ^6 (1-θ)^14

cur = 0.5
for i in range(n_iter):
    prop = cur + normal(0, step)               # ① 隣を提案する
    logr = log_target(prop) - log_target(cur)  # ② 比を取る → 正規化定数が消える
    if log(uniform()) < logr:                  # ③ 確率 min(1, r) で受容
        cur = prop
    chain[i] = cur                             # ④ 棄却でも今の値を記録する

②が全ての仕掛けです。 比を取ると

計算できない分母が約分で消えます。 分子（事前×尤度）は手で書けるので、これは計算できる。つまり正規化定数を知らないまま、その分布からサンプルできてしまう。 ここがベイズ計算の歴史を変えた一点です。

歩き方を言葉にするとこうです。

登るときは必ず行く。下るときはときどき行く。

$r \geq 1$（提案先の密度が高い）なら必ず受容（緑）、$r < 1$ なら確率 $r$ で受容し、残りは棄却（赤）。アニメには比 $r$ と受容確率 $\min(1, r)$ を並べて表示してあるので、$r$ が 1 を超えた歩がすべて受容されていることを確認できます。棄却されるとトレースプロットが横に平らになります。

初期値 0.80 は事後分布の端（密度がほぼ0の場所）なので、中央へ向かう提案は密度が上がる＝坂を登る方向になります。だから最初はほぼすべて受容されて、山の中心へ吸い込まれていきます。実際、最初の40歩のうち15歩が $r > 1$ で、そのすべてが受容されました。

④の「棄却でも記録する」を落とすとバグります。 棄却は「動かない」だけで、その場にもう1回滞在したことになる。記録しないと密度の高い場所の重みが失われます。

足跡を数えるだけで厳密な $\text{Beta}(7,15)$ に一致していきます。標本平均は10回で 0.762（大外れ）→ 200回で 0.344 → 4000回で 0.31807（厳密値 0.31818 との差 $-0.0001$）。

なぜ定常分布が事後分布になるのか：詳細釣り合い

「事後分布が定常分布になるような連鎖」をどう作るのか。詳細釣り合い（detailed balance） という条件を使います。

読み方は単純で、左辺が「$x$ にいて $y$ へ移る量」、右辺が「$y$ にいて $x$ へ移る量」。2点間の行き来が釣り合っているという意味です。全ての点対で流れが釣り合っていれば分布は変化しないので、$\pi$ が定常分布になります。

メトロポリスの受容確率 $\min\bigl(1, \pi(y)/\pi(x)\bigr)$ は、この等式が成り立つように逆算して作られたものです。天から降ってきた式ではありません。確かめました。

差は倍精度の丸め誤差だけです。ぴったり釣り合っています。

$q$ が両辺に出てくることに注意してください。 今回は提案が正規分布で対称 $q(x \to y) = q(y \to x)$ なので $q$ が約分で消え、比が $\pi$ だけで書けました。対称でない提案を使う場合は $q$ の比も残ります。それが「メトロポリス法」と「メトロポリス・ヘイスティングス法」の違いです。

ここで条件の向きを正確にしておきます。 詳細釣り合いは $\pi$ が定常分布であるための十分条件であって、必要条件ではありません。 実際この記事の後半で使うギブスサンプリングは、$\mu \to \tau$ という決まった順番で更新する（決定的スキャン）ので詳細釣り合いを満たしませんが、$\pi$ を不変に保ちます。「MCMC は必ず可逆」と覚えると間違えます。

もうひとつ、詳細釣り合いだけでは収束は保証されません。 第22回でやった既約性（どこからどこへでも行ける）と非周期性が別途必要です。後で出てくる二峰性の失敗例（$\hat R = 5.12$）は、まさに実質的に既約でなくなった状態です。

なお第23回のランダムウォークとの関係も見えます。提案 $y = x + \varepsilon$ はただのランダムウォークで、それに受容判定という重み付けを加えると、目標分布に沿って歩くようになる。だからこの方式はランダムウォーク・メトロポリスと呼ばれます。

バーンイン：最初のほうは捨てる

わざと外した初期値 0.97 から、しかも小さすぎる歩幅 0.01 で出発した例です。745 反復目でようやく 0.45 を下回ります。 この期間のサンプルは事後分布ではないので捨てます。

1500 回捨てると平均 0.31860 で厳密値 0.31818 と一致します。捨てないと 0.32423 とズレたままです。

歩幅のチューニング：受容率は高いほど良いのではない

反復回数は3つとも 60000 で同じなのに、質が全く違います。

歩幅 0.005 は受容率 98.5% でほぼ全部通りますが、一歩が小さすぎて全体を回れません。58000 個記録しているのに、独立サンプル換算では 34 個分しかない。

逆に歩幅 2.0 は提案が遠すぎてほぼ棄却され、トレースが階段状に固まります。

受容率 20〜50% あたりが目安です（1次元の理論的最適は約 44%、高次元では約 23.4% という結果が知られています）。「受容率が高いほど効率がいい」は直感的ですが逆で、ここは引っかかりやすい点です。

なお表を見ると、ESS が小さくても事後平均自体はそれなりに当たっています（歩幅0.005 で誤差 $-0.00096$）。ESS が効くのは「その推定値がどれくらい信用できるか」の側です。

マルコフ連鎖なので隣り合うサンプルは似ています。赤（歩幅が狭い）はラグ 300 でも相関 0.70 が残っている――300回前の値をまだ引きずっています。この「相関が消えるまでの長さ」が ESS を決めます。

まとめて検証する

サンプルの集まりが手に入ると、任意の質問が全部「数えるだけ」で答えられます。 積分をやり直す必要がありません。

共役でない実例に適用する

道具が揃ったので、手では絶対に解けない問題を解きます。ここで大事なのは、変えるのが2箇所だけだという点です。

prop = cur + rng.normal(0, 1, 2) * step        # ← 提案を2次元にした
if log(rng.random()) < log_post(*prop) - log_post(*cur):
    cur = prop
chain[i] = cur

アルゴリズム本体は1文字も変わりません。 差し替えるのは log_post の中身だけです。ここが MCMC が「あらゆるベイズモデルの標準解法」になった理由で、モデルごとに解法を考える必要がなくなりました。

例1：ロジスティック回帰をベイズでやる

第18回のロジスティック回帰を題材にします。記事の長さ $x$（千字）とクリックの有無 $y$（0/1）、$n=20$。

この尤度には、手で扱える（名前のついた）共役事前分布がありません（さきほど見たとおり、共役族そのものは存在しても正規化定数が閉じた形で書けません）。だから事後分布は手では出ません。パラメータは $b_0, b_1$ の2個です。

書くのは事後分布の分子だけです。

def log_post(b0, b1):
    eta = b0 + b1*x
    log_lik = np.sum(y*eta - np.logaddexp(0.0, eta))   # ロジスティックの対数尤度
    log_pri = -(b0**2 + b1**2)/(2*5.0**2)              # 事前 N(0, 5²)
    return log_lik + log_pri

これだけで動きます。

①データと当てはめ。 薄い青の線は事後分布から引いた120本の曲線で、これ自体が不確実性の可視化になっています。

②2次元の事後分布。 $b_0$ と $b_1$ は強く負に相関しています（相関 $-0.9034$）。細長い斜めの形です。直感的には、切片を下げても傾きを上げれば同じデータを説明できるからです。この「斜めさ」は後でギブスの弱点として効いてきます。

③$b_1$ の周辺事後分布。 MCMC のヒストグラムとグリッド積分が完全に重なっています。

なお今回はパラメータが2個なので、グリッド数値積分で答えを別途計算できます。 これを参照解にして MCMC を検証しました。

最後の行がベイズにしか書けない文です。「記事の長さの効果が正である確率は 99.9%」と直接言える。頻度論では p 値（データの珍しさ）しか出てきません。

最尤推定との比較：$n=20$ では結論がズレる

$n=20$ と小さいので一致しません。 $N(0, 5^2)$ の弱い事前分布が極端な値を抑える方向に効き、区間の下限が 0.152 → 0.383 と 0 から遠ざかっています。

さきほど「$n$ が大きければ一致する」と書きましたが、その裏返しがこれです。小標本でこそ差が出ます。

予測：サンプルを変換して数えるだけ

ここが実務で一番効く部分です。

「3千字の記事のクリック率は?」に答えるには、18万個のサンプルそれぞれを変換して数えるだけです。

pr3 = 1/(1+np.exp(-(s[:,0] + s[:,1]*3.0)))   # 全サンプルを変換
pr3.mean()                                    # 0.6989
np.percentile(pr3, [2.5, 97.5])               # [0.4325, 0.9148]
(pr3 > 0.5).mean()                            # 0.9304

参照解（グリッド積分）の事後平均 0.6993 と一致します。

予測の図の右がこれを全ての $x$ でやったものです。頻度論なら誤差伝播の公式（デルタ法）を導出する必要がありますが、サンプルがあると変換して数えるだけです。

この図で自分の思い込みが崩れました。 最初は「データが薄い両端で帯が広がる」と書いたのですが、実測すると逆でした。

計算した $x = 0$ から 6 の範囲では、最大は中央付近（$x \approx 1.8$）で、データが薄い右端がいちばん狭い。 理由は確率スケールだからです。$x$ が大きいと予測が 1 に張り付くので、上にも下にも動く余地がなくなります。

線形予測子（対数オッズ）のスケールで測ると期待通りの形になり、同じ範囲で幅は $x=2.0$ 付近が最小の 2.52、$x=6.0$ が最大の 7.68 と端で広がります。つまり「端で広がる」は対数オッズのスケールでの話で、確率に変換すると天井と床に潰されて逆転する、というのが正しい理解でした。

例2：ギブスサンプリング

もう1つの主要な MCMC 手法です。題材を変えます。身長データ $n=25$ で、平均 $\mu$ と標準偏差 $\sigma$ の両方が未知。

ギブスが使える条件はこれです。

片方を固定すると、もう片方が共役になる。

ここでは $\sigma$ を固定すれば $\mu$ は正規の共役、$\mu$ を固定すれば精度 $1/\sigma^2$ はガンマの共役になります。さきほどの共役の表がそのまま部品として使われます。

for i in range(n_iter):
    # ① σ を固定して μ を引く（正規の共役）
    prec = 1/s0**2 + n*tau
    mu  = rng.normal((m0/s0**2 + tau*data.sum())/prec, sqrt(1/prec))
    # ② μ を固定して τ=1/σ² を引く（ガンマの共役）
    tau = rng.gamma(a0 + n/2, 1/(b0 + ((data-mu)**2).sum()/2))

受容判定がありません。 棄却が0回、受容率 100% です。「提案して受け入れるか決める」のではなく、正しい分布から直接引いているので捨てる必要がありません。

ギブスは「横→縦」の2段階でしか動きません。 斜めには決して動かない――一度に1つのパラメータしか更新しないからです。

軌跡を並べた図の下段のトレースで差がはっきりします。メトロポリス（右）は平らな段が目立ちます（棄却されて動けなかった時間）。ギブスは（連続パラメータなので）毎回必ず動きます。

参照解と一致しています。$\sigma$ の事後分布が右に裾を引いて非対称なのが見どころで、「標準偏差の不確実性は上側に大きい」という、当然だが式では見えにくい事実が図に出ます。

効率の比較：同じ問題・同じ反復数

反復数は両方 195,000 で同じなのに、ギブスの有効サンプル数は 7.8〜10.1倍。しかも ESS が反復数にほぼ等しい＝ほとんど独立なサンプルが得られているということです。

ただし、どちらも正しい答えに収束しています。 差は「正しさ」ではなく「効率」です。ここは混同しやすい点なので強調しておきます。

使い分け

条件付き分布が書けるならギブス、書けないならメトロポリス。 これが判断基準です。

なお両者は別系統の手法ではありません。ギブスは受容確率が恒等的に 1 になるメトロポリス・ヘイスティングス法の特別な場合として導けます。だから「棄却がない」のも偶然ではなく、条件付き分布から直接引いている以上、必ず受け入れられるということです。

さきほど「決定的スキャンのギブスは詳細釣り合いを満たさない」と書いたことと矛盾するように見えますが、両方正しいです。各成分の更新それぞれは詳細釣り合いを満たす（だから MH の特別な場合と言える）のですが、決まった順番に合成したスイープ全体としては満たさない（$\mu \to \tau$ の逆順にはならないため）。それでも各ステップが $\pi$ を不変に保つので、合成も $\pi$ を不変にします。

ギブスの弱点は、さきほどのロジスティック回帰の2次元事後分布に繋がります。パラメータ間の相関が強い（細長い斜めの分布）とギブスは極端に遅くなります。 軸に平行にしか動けないので、斜めの谷を進むのに小刻みなジグザグを繰り返すしかありません。例1のロジスティック回帰（相関 $-0.9034$）でギブスを使わなかったのは、条件付きが共役でないことに加えてこの理由もあります。

なお両者は排他ではなく、一部のパラメータだけメトロポリスで更新する「メトロポリス内ギブス」が実務では多用されます。実際の統計ソフトでは、JAGS や PyMC が変数ごとにサンプラーを割り当てる一方、Stan は一律にハミルトニアン・モンテカルロ法（HMC）とその発展版 NUTS を使います（そのため離散パラメータを直接は扱えません）。

収束判定は主観的なのか

トレースプロットを目で見るだけでは不十分です。なぜ不十分なのかが今回いちばん納得できた部分でした。

二峰性の分布（山が2つ）で、4本の鎖を別々の初期値から走らせました。

上段が失敗例（歩幅 0.5）。4本が2つの山に分断され、互いに行き来できていません。鎖ごとのヒストグラムを見ると、各鎖は片方の山しか見ていない。$\hat R = 5.1247$。

下段が成功例（歩幅 3.0）。山を越えられるので4本とも同じ分布に到達。$\hat R = 1.0010$。

$\hat R$ は何を測っているのか

鎖と鎖の間のばらつき $B$ と 1本の鎖の中のばらつき $W$ を比べています。

• 収束していれば、どの鎖も同じ分布を見ている → $B \approx W$ → $\hat R \approx 1$
• 収束していなければ、鎖が別々の場所に閉じ込められている → $B \gg W$ → $\hat R > 1$

実測では失敗例で $W = 0.4761$ に対して $B = 120{,}272$（25万倍）。「トレースプロットを目で見る」を数値化したものが $\hat R$ です。判定基準は慣習的に 1.01 未満（昔は 1.1 でしたが、近年は厳しめが推奨されます）。

なぜ複数の鎖が必要か

ここが最重要です。 失敗例の第1鎖だけを見ると、

• 平均 $-3.0006$、標準偏差 0.6949
• 前半と後半の平均差 0.0204

完璧に安定して見えます。 ところが第3鎖は平均 $+3.0110$ で真逆の場所にいます。

この例では1本のトレースプロットからは気づけません。 だから「目で見る」だけでは不十分で、複数鎖＋$\hat R$ が必要になります。

そして $\hat R$ の図の右列の通り、失敗例では反復を増やしても $\hat R$ が下がりません（5.0 のまま）。

これが「いくら回しても無駄」のサインです。 収束していない鎖は、時間をかけても収束しません。

さらに厄介な罠

この失敗例、4本まとめた平均は $+0.0047$ で真の値 0 とほぼ一致してしまいます。 山が対称なので偶然打ち消し合ったのです。

しかし分布は完全に間違っており、$P(-1 < \theta < 1)$ は実測 0.00160 対 真値 0.00213 でズレています。

「平均が合っているから大丈夫」は収束の証拠にならない。

これは自分でやってみないと気づけない類の罠でした。$\hat R$ を見ていれば 5.12 で即座に分かります。

階層ベイズ：他の記事の情報を借りる

ここからはブログのアクセス解析に直接使える話です。試験では用語レベルの扱いですが、実務でベイズを選ぶ理由がここに集中しているので数値で追いました。

困りごきから

記事30本のクリック率を知りたい。ただし表示回数がバラバラです（5回のものから2000回のものまで）。総表示 8,749回・総クリック 1,618回・全体率 0.1849。

表示8回でクリック0回の記事は「クリック率0%」なのか。 明らかに違うのに、その記事のデータだけを見ると $0/8 = 0$ しか出てきません。

やり方は3つあります。

「階層」とは何が階層なのか

確率の指定が3段重なっていることを「階層」と呼びます。

下から読むと「データは $\theta$ で決まる」「その $\theta$ は $\text{Beta}(a,b)$ から来ている」「その $a, b$ は…」と積み上がります。推定はこの矢印を逆にたどる作業です。

「記事ごとに違うが、無関係でもない」を階層で書いたのがこのモデルです。①は「全部無関係」、②は「全部同じ」と言っていたので、その中間を表現できるようになったことになります。

$a, b$ は何を表しているのか

ここが私がいちばん詰まった箇所でした。

$\text{Beta}(a, b)$ は「30本の記事の真のクリック率がどう散らばっているか」を表す分布。

記事1本の率ではなく、率の集団を表しています。第2層に置かれているのが目印です。読み方は2通りあります。

読み方A：中心と集中度に分ける（意味を掴むならこちら）

つまり $\text{Beta}(3.49, 16.68)$ は「記事のクリック率は 0.173 ± 0.082 くらいに散らばっている」という主張です。実際の30本の真の率は平均 0.1808・標準偏差 0.0750 なので、データから正しく読み取れています。

読み方B：疑似データとして読む（計算に直結するのはこちら）

この記事の前半で出てきたラプラススムージング（分子に1、分母に2を足す）は $a=b=1$ の場合でした。経験ベイズは「その足す数をデータから決める」版です。

実際の計算を1本の記事で追う

対象： 表示 $n=8$ 回、クリック $x=0$ 回（真の率は 0.0919）

ステップ1：標本比率

これが「クリック率0%」という無意味な答えになる原因です。

ステップ2：事前分布を掛けて事後分布を作る

ベータ分布と二項分布は共役なので、足すだけです。

ステップ3：その平均を取る

真の率 0.0919 にかなり近い値です。

図の中央パネルが要点です。$n$ が小さいと尤度自身の幅が非常に広く、「0 だと断定できない」と言っているので、事前分布に引かれる余地が生まれます。

ステップ4：重み付き平均に書き換える

ステップ3と一致します。自分のデータの重みは 28.4% しかなく、71.6% は全体平均から借りている。 これが「他の記事の情報を借りる」の具体的な意味です。

なぜ2つの式が同じなのかは、分数を書き換えるだけです。

1行目から2行目は分子を分けただけ、2行目から3行目は各項に $n/n$ と $(a+b)/(a+b)$ を掛けただけです。「縮小」という高級そうな話が、ただの分数の書き換えだったのがこの単元の気持ちよいところでした。

ステップ4は必要なのか（自分が引っかかった点）

ここで私はひとつ勘違いをしました。「ステップ4は本来必要な計算だが、今回はベータだから省略できるのだろう」と思ったのです。逆でした。

ステップ4はどんな場合でも計算には不要。ただし共役のときだけ、それが厳密な等式として書ける。

ステップ3で答えはもう出ています。ステップ4は同じ数を別表記にしただけで、新しい情報を1つも生みません。「$w$ だけデータを信じ、残りは全体平均から借りている」という意味を見えるようにするための書き換えです。

つまり共役だと得なことが2つあります。

(B) を実験しました。

ベータでは倍精度の限界で完全一致、ロジット正規では $10^{-3}$ オーダーのずれが残ります。

このきれいな直線の形は共役の副産物です。共役でないと、事後平均は標本比率と事前平均の線形結合では書けません（尤度と事前分布の形が非線形に絡むため）。

これは第19回のリッジ回帰との類似にも前提が付くことを意味します。 「縮小＝重み付き平均」が厳密に成り立つのは、共役や正規モデルの場合だけです。

縮小の効果を測る

①も②も負けて、その中間の③が勝ちます。 精度は 2.27 倍。表示20回以下の記事に絞ると 34.3%（約3倍）まで改善します。逆に表示500回以上では 75.8% で差が小さくなります。

$0/8$ と $0/10$ の記事は標本比率だと 0 になりますが、経験ベイズでは 0.124・0.116 に持ち上げられます（真の率は 0.092・0.143）。

推定された事前分布が「20.2回分の疑似データ」に相当するので、$n = 20.2$ でちょうど半々になります。

同じ式・同じ $a, b$ を使っているのに、$n$ が大きい記事はほとんど動きません。

「$n$ に応じて自動で調整」の中身はこれだけです。足している疑似データが 20.2 回分なので、2000 回の実データの前では無視できる量になります。

これは第19回のリッジ回帰・Lasso と同じ話です。 あれも係数を 0 のほうへ縮めて精度を上げていました。「少し偏らせる代わりに分散を大きく減らす」というバイアス・バリアンストレードオフが、ベイズでは事前分布という形で自然に出てきます。

階層ベイズと経験ベイズは何が違うのか

ここまでの計算は経験ベイズでした。階層ベイズとの関係を整理します。

階層モデルの立場からは、経験ベイズは第3層を1点で近似したものと読める。

ただし「経験ベイズは階層ベイズの劣った版」という意味ではありません。経験ベイズは Robbins や James-Stein の流れで、頻度論的なリスク（平均二乗誤差）の観点から独立に正当化される手法です。事前分布の形を仮定しないノンパラメトリック経験ベイズもあります。

違いは第3層の扱いだけです。

$a, b$ の決め方（経験ベイズ） は、各記事の $\theta_i$ を積分で消してしまい、$a, b$ だけの関数を最大化します。

$\theta_i$ が消えているので2変数の最適化で済みます。これが経験ベイズが軽い理由です。

ただしここに弱点があります。

$a$ を 3.0 や 4.0 にしても周辺尤度は 0.7 しか下がりません。 つまり $(a,b)$ はあまり精度よく決まっていないのに、経験ベイズは 3.493 という1点に決め打ってしまう。この決め打ちの不確実性が結果に反映されません。

階層ベイズを実際に計算する

$a, b$ に超事前分布 $p(a,b) \propto (a+b)^{-5/2}$ を置き、MCMC で推定しました。$\theta_i$ は共役なので解析的に積分でき、$(a,b)$ の2次元だけメトロポリス法で動かせば済みます。 受容率 32.6%、4本の鎖で $\hat R = 1.0010$（$a$）・1.0011（$b$）で収束を確認しました。

アルゴリズムは2段構えです。

for t in range(n_iter):
    # 【前半】a, b をメトロポリス法で1歩動かす（受容/棄却あり）
    prop = (a, b) + 提案のノイズ
    if log(uniform()) < log_marg(prop) - log_marg(cur):
        a, b = prop

    # 【後半】今の a, b で θ_i を1個引く ← これが「ステップ2」の分布
    for i in range(30):
        theta[t, i] = rng.beta(x[i] + a,  n[i] - x[i] + b)

後半の Beta(x+a, n−x+b) が、手計算のステップ2で作った分布そのものです。 共役なので「足すだけ」で書けるため1行で済んでいます。そして前半で $a, b$ が毎回変わるので、後半で使う分布も毎回変わります。

実際の反復を並べるとこうなります（$n=8$, $x=0$ の記事について）。

経験ベイズなら $a=3.4933$, $b=16.6773$ 固定なので、右列が全部 $\text{Beta}(3.493, 24.677)$ の同じ分布になります。

「毎回違う分布から1個ずつ引いて溜める」――これが階層ベイズの MCMC がやっていることの全部です。

結果として事後分布は「ベータ分布の混合」になる

• 経験ベイズの事後分布 = 1つのベータ分布
• 階層ベイズの事後分布 = ベータ分布の混合分布

「混合」なので、成分の平均的な分散よりは必ず広くなります（全分散の公式で $V = E[V] + V[E] \geq E[V]$）。ただし比較相手である経験ベイズの1本は $V(\theta \mid \hat a, \hat b)$ で $E[V]$ とは別物なので、「必ず広い」が自動的に言えるわけではありません。この例では実測で広くなっていることを後で確認します。

差が出るのは点推定ではなく区間幅

超パラメータの事後分布

$a$ の95%区間が [1.39, 6.53] と非常に広い。 記事30本しかないので $a, b$ は精度よく決まりません。

点推定はほぼ同じ

差は最大 0.00495・平均 0.00164 で、実用上区別がつきません。精度を上げたいだけなら経験ベイズで十分です。

違いが出るのは区間の幅

階層ベイズの区間のほうが広い。 $a, b$ 自体の不確実性を織り込んでいるからです。逆に言えば経験ベイズは自信過剰（不確実性を1段階分数え落としている）。

95%区間が真の率を含んだ割合はどちらも 96.7%（29/30本）でした。この例では実害は出ませんでしたが、$n$ がもっと小さい・グループ数がもっと少ない場合は問題になります。

余分な不確実性を分解する

全分散の公式で2つに分けられます。

すべて「データを与えたもとで」の量です。以下では表記を簡単にするため条件付けを省略します。

経験ベイズは第2項を丸ごと落としています。 $a, b$ を1点に決め打った瞬間に $V[E(\theta \mid a,b)] = 0$ になるからです。

検算として、MCMC サンプルでの実測（平均 0.119194・標準偏差 0.064721）と厳密な混合分布の値（平均 0.119450・標準偏差 0.065186）が一致しました。経験ベイズの標準偏差は 0.061024 なので、階層ベイズは 1.068 倍です。

過小評価の度合いは $n$ が小さいほど大きくなります（$n=5$ で 1.094 倍、$n=2000$ で 1.0008 倍）。$n$ が大きい記事は $a, b$ にほとんど頼っていない（重み $w$ が 99%）ので、$a, b$ の不確実性も効かないのです。

ただし単調ではありません。 30本を $n$ の順に並べると9箇所で増加しており（$n=25$ の 1.0361 から $n=30$ の 1.0803 など）、$x$ の値によっても変わります。$n$ と比の順位相関は $-0.9564$ で、「おおむね減るが単調ではない」が正確な表現です。

同じ理由で、この比は第2項の割合だけでは決まりません。$n=8$ なら $\sqrt{1/(1-0.0792)} = 1.042$ のはずですが実測は 1.068 です。第1項 $E[V(\theta \mid a,b)]$ も経験ベイズの $V(\theta \mid \hat a, \hat b)$ と等しくない（$a, b$ を平均する操作が条件付き分散そのものも動かす）ためで、この点は最初に書き間違えて検算で気づきました。

この「二段階」の正体は周辺化だった

階層ベイズの MCMC を見ていて、いちばん面白かったのがここです。

最終的に欲しいのは $a, b$ ではありません。 $a, b$ の分布を求めて、それに従って引いた $\theta$ の分布が欲しい。この二段階の構造に、統計学は名前を持っています。

式で書くとこうです。

MCMC はこの二重積分を「引いて数える」で自動的にやっています。

• アニメの前半（メトロポリス法）＝ $p(a,b \mid \text{データ})$ から引く
• アニメの後半（ベータから抽出）＝ $p(\theta \mid a,b,\text{データ})$ から引く

「二段階」に見えたのは正しく、それが積分の構造をそのまま反映しているからでした。

そして経験ベイズは、この式の後半を1点に集まった分布で近似しています。

デルタ関数で置き換えると積分が消えて、ベータ分布1本だけが残る。 「経験ベイズ＝1本、階層ベイズ＝混合」の正体はこれでした。

同じ構造を第7回で既に見ていた

ここで伏線が回収されます。t分布がまさにこの構造です。

第7回で t分布を「分散を知らない罰金」として学びました。その中身はこうです。

$\sigma^2$ を周辺化すると、ちょうど自由度 $\nu$ の t分布になります。つまり

t分布 = 正規分布の混合分布

今回の「階層ベイズの事後分布 = ベータ分布の混合」とまったく同じ構造です。

実験で確認しました。

左の列は400万個の乱数によるシミュレーション値なので、$10^{-3}$ 程度の差は乱数誤差です。

分散も一致します（$\nu=5$ で混合 1.6641 対 理論値 $\nu/(\nu-2) = 1.6667$）。

そして「決め打つと自信過剰になる」のも同じ

（この列は数表と照合できるよう厳密値を載せています。）

「なぜ $n$ が小さいと $z$ ではなく $t$ を使うのか」の答えがこれです。 $\sigma^2$ を1点（最尤推定値）に決め打つと、区間が $\nu=3$ で 38.4% も短くなります。

そして階層ベイズで見た「経験ベイズは区間が 5% 狭い」とまったく同じ種類の誤りです。一般化するとこうなります。

消すべきパラメータを最尤推定値のような中心的な1点に決め打つと、通常は自信過剰（区間が狭すぎる）になる。 正しく周辺化すると裾が厚くなる。

ただし「必ず」ではありません。代入する点の選び方で向きが変わります。 $\sigma^2$ を最尤推定値ではなく事後平均 $E[\sigma^2] = \nu/(\nu-2)$ で代入すると、$\nu=3$ では $1.96\sqrt{3} = 3.3948$ となり、t の 3.1824 より6.7% 広くなります（自信過小）。$\nu=5$ では 2.5303 対 2.5706 で 1.6% 狭いだけです。

理由は階層ベイズのところで自分で見つけた不一致と同じで、$V(\theta \mid \hat a, \hat b) \neq E[V(\theta \mid a,b)]$ だからです。「1点に潰すと $V[E(\cdot)]$ の項が消える」のは常に正しいが、同時に第1項も動くので、差し引きの符号は代入する点に依存します。

この構造が現れる場所の一覧