← ベイズ統計 一覧

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：なぜサンプリングか　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計・詳細釣り合いの証明）

要点（BLUF）

• メトロポリス・ヘイスティングス（MH）は、規格化されていない事後（尤度×事前）の値だけからサンプルを取る最も基本的な MCMC。受容確率に事後が比で入るので、解けない正規化定数 $Z$ が約分で消えます。
• アルゴリズムは「提案 → 受容確率 $\min(1,\text{比})$ で採否」。密度が下がる方向も比の確率で受け入れるから、山頂に貼り付かず分布全体を歩けます。
• 提案幅（step）が要調整：小さすぎると高受容だが動かず（相関大）、大きすぎると低受容。numpy 実装で $\mathrm{Beta}(16,8)$ を再現し、この調整を体験します。

1. アルゴリズム

現在地 $\theta$ から次の点を決める手順です（詳細釣り合い→定常分布の証明は MCMC（マルコフ連鎖モンテカルロ））。

1. 提案：対称な提案分布から候補を出す。例：ランダムウォーク $\theta'=\theta+\varepsilon,\ \varepsilon\sim\mathcal N(0,\text{step}^2)$。
2. 受容確率：対称提案なら $\displaystyle\alpha=\min\!\left(1,\ \frac{\pi(\theta')}{\pi(\theta)}\right)$。
3. 採否：$u\sim\mathrm{Uniform}(0,1)$ を引き、$u\le\alpha$ なら受容して移動、そうでなければ棄却して現在地にとどまる（その点も標本に記録）。

$\pi$ は事後 $p(\theta\mid D)=p(D\mid\theta)p(\theta)/Z$ ですが、受容確率では比なので

と $Z$ が消えます。計算できる分子（尤度×事前）だけで回せる——これが なぜサンプリングか の壁を越える仕掛けです。

2. numpy でフル実装し、Beta(16,8) を再現する

正規化していない事後 $\pi(\theta)\propto\theta^{15}(1-\theta)^{7}$（$\mathrm{Beta}(2,2)$ 事前 × 二項 $n{=}20,k{=}14$）を目標に、MH で標本を取り、真の事後 $\mathrm{Beta}(16,8)$ と一致するか見ます。対数で計算して桁あふれを防ぎます。

import numpy as np
from scipy import stats

# 目標：正規化なしの対数事後 log π(θ) = 15 log θ + 7 log(1-θ)  （[0,1] 外は -inf）
a0, b0, n, k = 2, 2, 20, 14
def log_target(t):
    if t <= 0 or t >= 1:
        return -np.inf
    return (a0 - 1 + k)*np.log(t) + (b0 - 1 + n - k)*np.log(1 - t)

def metropolis(step, N=40000, burn=2000, seed=0):
    rng = np.random.default_rng(seed)
    th = 0.5; lt = log_target(th)
    out = np.empty(N); acc = 0
    for i in range(N):
        cand = th + rng.normal(0, step)            # 対称ランダムウォーク提案
        lc = log_target(cand)
        if np.log(rng.uniform()) < lc - lt:        # u<exp(Δlog) ⇔ 受容
            th, lt = cand, lc; acc += 1
        out[i] = th                                # 棄却時は現在地を記録
    return out[burn:], acc/N

true = stats.beta(16, 8)
print(f"真の事後 Beta(16,8): 平均={true.mean():.4f}, 95%=[{true.interval(0.95)[0]:.3f}, {true.interval(0.95)[1]:.3f}]")
for step in [0.02, 0.1, 0.4]:
    s, ar = metropolis(step)
    lo, hi = np.quantile(s, [0.025, 0.975])
    print(f"  step={step:<4}: 受容率={ar:.2f}  標本平均={s.mean():.4f}  95%=[{lo:.3f}, {hi:.3f}]")

出力：

真の事後 Beta(16,8): 平均=0.6667, 95%=[0.471, 0.836]
  step=0.02: 受容率=0.93  標本平均=0.6699  95%=[0.463, 0.835]
  step=0.1 : 受容率=0.69  標本平均=0.6683  95%=[0.478, 0.837]
  step=0.4 : 受容率=0.28  標本平均=0.6663  95%=[0.469, 0.836]

出力の意味：どの step でも標本平均 $\approx0.667$・95%区間 $\approx[0.47,0.84]$ と、真の事後 $\mathrm{Beta}(16,8)$ を正しく再現しています。正規化定数を一度も計算せず、$\theta^{15}(1-\theta)^7$ の値だけでこれができたのが MH の威力です。一方で受容率は step とともに $0.93\to0.28$ と下がります。次節でこの調整を見ます。

3. 提案幅のトレードオフと、トレースで見る歩き方

• step 小：候補が現在地のすぐ近く → ほぼ受容（$0.93$）。でも一歩が小さく、サンプルが互いに強く相関して分布を回りきるのが遅い。
• step 大：遠くへ跳ぶ → 多くが密度の低い所で棄却（受容率 $0.28$）、同じ点に長くとどまる。
• 経験則として、1次元では受容率 $0.3$〜$0.5$ あたりが効率の良い目安（高次元では $\approx0.234$）。相関の強さは有効サンプルサイズで定量化します（収束診断）。

サンプルの時間変化（トレース）とヒストグラムを描くと、連鎖が事後を埋めていく様子が見えます。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用

a0, b0, n, k = 2, 2, 20, 14
def log_target(t):
    if t <= 0 or t >= 1: return -np.inf
    return (a0-1+k)*np.log(t) + (b0-1+n-k)*np.log(1-t)

rng = np.random.default_rng(0)
th, lt, step = 0.5, log_target(0.5), 0.4
chain = np.empty(40000)
for i in range(40000):
    cand = th + rng.normal(0, step); lc = log_target(cand)
    if np.log(rng.uniform()) < lc - lt: th, lt = cand, lc
    chain[i] = th
chain = chain[2000:]                                  # バーンイン除去

fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].plot(chain[:1500], lw=0.6); ax[0].set_title("トレース（最初の1500ステップ）")
ax[0].set_xlabel("ステップ"); ax[0].set_ylabel("θ")
ax[1].hist(chain, bins=60, density=True, alpha=0.6, label="MH標本")
xx = np.linspace(0, 1, 300)
ax[1].plot(xx, stats.beta(16, 8).pdf(xx), lw=2, label="真の事後 Beta(16,8)")
ax[1].set_title("標本ヒストグラム vs 真の事後"); ax[1].set_xlabel("θ"); ax[1].legend()
plt.tight_layout(); plt.show()

グラフの意味：左のトレースは水平な帯状にばらつき（定常に達したサイン）、上下に揺れながら事後の台 $[0.45,0.85]$ を往復します。右のヒストグラムは真の事後 $\mathrm{Beta}(16,8)$（実線）にぴたりと重なります。トレースがドリフトしたり一点に張り付いたりせず帯状なら「収束していそう」と読めます（厳密な判定は 収束診断）。

⚠️ よくある誤解

• 「棄却したら記録しない」ではない：棄却時は現在地と同じ値を1標本として記録します。滞在時間が事後確率の高さを表すので、これを落とすと分布が歪みます（MCMC（マルコフ連鎖モンテカルロ））。
• 「受容率は高いほど良い」ではない：受容率 $0.99$ は step が小さすぎて動いていないサイン。高すぎず低すぎず（1次元で $0.3$〜$0.5$）が目安。
• 「正規化定数が要る」ではない：受容確率は比なので $Z$ は消えます。計算できない $Z$ を回避できるのが MH の存在理由（§1）。
• 「1本の連鎖が帯状なら必ず収束」ではない：多峰分布では1本だと片方の山に閉じ込められることがあります。複数チェーンと $\hat R$ で確認します（収束診断）。

関連ノート

• なぜサンプリングか
• ギブスサンプリング（次のトピック・条件付き分布で受容判定なし）
• 収束診断（受容率・相関・トレースの定量評価）
• MCMC（マルコフ連鎖モンテカルロ）（統計・受容確率が詳細釣り合いを満たす証明）
• 第4章 MCMC 目次
• ベイズ統計テキスト 全体目次