事後分布の要約｜ベイズ統計

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

📎 前提：指数型分布族と共役事前　|　関連：点推定と損失関数・信用区間と事後予測分布

要点（BLUF）

事後分布は点・区間・確率の3通りで要約します。点は平均/中央値/最頻値（歪むと3つはずれる）、区間は等裾とHPD（最高密度）、確率は $P(\theta\in A\mid D)$ を直接読む。
HPD区間は同じ確率を最短幅で覆う。歪んだ事後では等裾より狭くなり、サンプルからの「最短区間」アルゴリズムで求められます。
ベイズの強みは仮説の確率をそのまま計算できること。「 $P(\theta>c\mid D)$ 」は事後を積分（またはサンプルの割合）するだけ。このサンプルからの要約が、第4章 MCMC の出口になります。

1. 3通りの要約

第1章では Beta 事後を例に点推定（点推定と損失関数）と信用区間・事後予測（信用区間と事後予測分布）を見ました。本ノートはそれを要約の道具箱として整理し、共役が使えない一般の事後（第4章でサンプルとして得る）にも効く形にします。例として、歪んだ事後 $\mathrm{Gamma}(25,11)$ （ガンマポアソンと指数のポアソン事後）を使います。

要約のタイプ	何を答えるか	代表
点	「いちばんありそうな1値」	事後平均・中央値・最頻値(MAP)
区間	「だいたいこの範囲」	等裾区間・HPD区間
確率	「ある命題が成り立つ確率」	$P(\theta>c\mid D)$ 、 $P(a<\theta<b\mid D)$

2. 点要約：歪むと平均・中央値・最頻値はずれる

対称な事後なら3つは一致しますが、歪んだ事後ではずれます。どれを報告すべきかは損失で決まります（点推定と損失関数：二乗損失→平均、絶対損失→中央値、0-1損失→最頻値）。

import numpy as np
from scipy import stats

# 歪んだ事後（ガンマ–ポアソン事後 Gamma(25,11)）
post = stats.gamma(a=25, scale=1/11)        # scale = 1/rate

mean   = post.mean()
median = post.median()
mode   = (25 - 1) / 11                        # ガンマの最頻値 (α-1)/β
print(f"事後平均={mean:.4f}  中央値={median:.4f}  最頻値(MAP)={mode:.4f}  SD={post.std():.4f}")

出力：

事後平均=2.2727  中央値=2.2425  最頻値(MAP)=2.1818  SD=0.4545

出力の意味：右に裾を引く分布なので 平均 $2.27$ ＞中央値 $2.24$ ＞最頻値 $2.18$ の順。裾が平均を引き上げます。1点だけ報告するときは、この3つのどれかを目的（損失）に応じて選び、必ず散らばり（SD や区間）を添えます。

3. 区間：等裾 vs HPD（最短区間）

同じ 95% でも区間の取り方は一意でありません。

等裾区間：左右の裾を $2.5\%$ ずつ切る。分位点で即座に求まる。
HPD区間（最高密度）：区間内の密度が区間外より高くなるよう取る。同じ確率を最短幅で覆う。

歪んだ事後では HPD が狭くなります。HPD は分位点では出せないので、**サンプルを並べて「全体の95%を含む最短の窓」**を探します。

import numpy as np
from scipy import stats

post = stats.gamma(a=25, scale=1/11)

# 等裾 95%
lo_et, hi_et = post.ppf(0.025), post.ppf(0.975)
print(f"等裾95% = [{lo_et:.4f}, {hi_et:.4f}]  幅={hi_et-lo_et:.4f}")

# HPD 95%：サンプルを並べ、95%を含む最短区間を探す
rng = np.random.default_rng(0)
s = post.rvs(200000, random_state=rng)
def hpd(samples, cred=0.95):
    s = np.sort(samples)
    n = len(s); k = int(np.floor(cred * n))   # 区間に含める点数
    widths = s[k:] - s[:n-k]                   # 幅 k 点ぶんの窓を全位置で
    i = np.argmin(widths)                      # 最短の窓
    return s[i], s[i+k]
lo_h, hi_h = hpd(s, 0.95)
print(f"HPD95%  = [{lo_h:.4f}, {hi_h:.4f}]  幅={hi_h-lo_h:.4f}")

出力：

等裾95% = [1.4708, 3.2464]  幅=1.7756
HPD95%  = [1.4201, 3.1817]  幅=1.7617

出力の意味：HPD の幅 $1.762$ は等裾 $1.776$ より狭い。歪んだぶん、HPD は左側へ少しずれて密度の高い領域を優先的に覆います。対称な事後なら両者は一致します。この「サンプルを並べて最短窓を探す」やり方は、第4章で MCMC サンプルしか手元にない場合にそのまま使えます（arviz.hdi も内部は同じ発想。PyMC/ArviZ の最新 API は要最新確認）。

4. 確率の問い：仮説の確率を直接読む

ベイズのいちばん自然な出力は「ある命題が成り立つ確率」です。頻度論では言えない「 $P(\theta>c\mid D)$ 」を、事後の積分（またはサンプルの割合）でそのまま計算します。

import numpy as np
from scipy import stats

post = stats.gamma(a=25, scale=1/11)
rng = np.random.default_rng(0)
s = post.rvs(200000, random_state=rng)

# 「発生率が 2.5 を超える確率」を 閉形式 と サンプル で
p_closed = post.sf(2.5)              # 生存関数 = P(λ>2.5)
p_sample = np.mean(s > 2.5)          # サンプルの割合
print(f"P(λ>2.5|D): 閉形式={p_closed:.4f}  サンプル={p_sample:.4f}")

出力：

P(λ>2.5|D): 閉形式=0.2910  サンプル=0.2928

出力の意味：「発生率が $2.5$ を超える確率は約 $29\%$ 」と、結論をそのまま確率で言えます。閉形式（ $0.291$ ）とサンプル割合（ $0.293$ ）はモンテカルロ誤差の範囲で一致。サンプルさえあれば不等式の割合を取るだけなので、共役で閉形式が出ない事後でも、MCMC サンプルから同じように確率を読めます。これが「事後分布を持つと意思決定に直結する」というベイズの実利です。

5. すべて「サンプルからの要約」に統一できる

§2〜§4 を見ると、点・区間・確率はどれもサンプルがあれば計算できます——平均は標本平均、分位点は順序統計量、HPD は最短窓、確率は割合。共役なら閉形式とサンプルが一致しますが、共役が崩れる一般のモデルでは閉形式が無く、サンプルだけが頼りになります。そのサンプルを事後から生成する手段が第4章の MCMC です。本章の要約の道具は、MCMC の出力をそのまま受け取れるように作ってあります。

flowchart LR
  P["事後分布 P(θ|D)"] --> S["サンプル θ⁽¹⁾…θ⁽ᴺ⁾"]
  S --> Pt["点：標本平均・中央値・最頻"]
  S --> In["区間：分位点（等裾）・最短窓（HPD）"]
  S --> Pr["確率：割合 #(θ∈A)/N"]

⚠️ よくある誤解

「平均・中央値・最頻値はどれも同じ」ではない：歪んだ事後ではずれます（§2）。対称分布のときだけ一致。報告は目的（損失）に応じて選び、散らばりを必ず添えます。
「等裾区間が唯一の信用区間」ではない：HPD のほうが狭く、歪んだ事後では密度の低い値を等裾が含むこともあります（信用区間と事後予測分布）。
「 $P(\theta>c\mid D)$ は頻度論でも言える」ではない：これはベイズ特有の読み。頻度論では $\theta$ は固定値で確率を割り当てません（信用区間と信頼区間で対比）。
「サンプル要約は閉形式の劣化版」ではない：サンプル数を増やせば任意精度に近づき、閉形式が無いモデルでは唯一の道。MCMC の収束（収束診断）さえ確かめれば信頼できます。