なぜサンプリングか

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

📎 前提：事前の選び方と感度分析　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計）・事前・尤度・事後・周辺尤度

要点（BLUF）

共役でない事後は正規化定数 $Z=\int p(D\mid\theta)p(\theta)\,d\theta$ が解けず、欲しい期待値もまた積分。高次元では数値積分が次元の呪いで破綻します。
代わりに事後からサンプルを取り、期待値を標本平均で近似するのがモンテカルロ積分。誤差は $\propto 1/\sqrt{N}$ で、しかも次元に依りません。これが高次元で唯一現実的な道。
残る問題は「正規化されていない事後から、どうやってサンプルを取るか」。それを解くのが MCMC（第4章の本体）です。

1. 壁：正規化定数と期待値の積分

ベイズで欲しいのは事後そのものより、そこからの期待値です（事後平均・信用区間・確率、事後分布の要約）。

\mathbb{E}[g(\theta)\mid D]=\int g(\theta)\,p(\theta\mid D)\,d\theta =\frac{\int g(\theta)\,p(D\mid\theta)\,p(\theta)\,d\theta}{\displaystyle\int p(D\mid\theta)\,p(\theta)\,d\theta}

分子も分母も積分です。分母は周辺尤度 $Z=p(D)$ （事前・尤度・事後・周辺尤度）。共役なら閉形式（第2章）ですが、現実のモデルはほとんど非共役で、 $Z$ が解析的に解けません。残された手は数値積分ですが、それも次元が上がると破綻します。

2. 次元の呪い：グリッドは爆発する

数値積分の素朴な方法は、パラメータ空間を格子に切って和を取ることです。しかし $d$ 次元を各軸 $G$ 分割で覆うには $G^d$ 個の格子点が要り、 $d$ とともに指数的に爆発します。

import numpy as np

G = 50
print(f"各軸 {G} 分割で d 次元を覆う総格子点数 G^d:")
for d in [1, 2, 3, 6, 10]:
    print(f"  d={d:>2}: {G**d:.3e} 点")

出力：

各軸 50 分割で d 次元を覆う総格子点数 G^d:
  d= 1: 5.000e+01 点
  d= 2: 2.500e+03 点
  d= 3: 1.250e+05 点
  d= 6: 1.562e+10 点
  d=10: 9.766e+16 点

出力の意味： $d=1$ なら $50$ 点で済むのに、 $d=6$ で $1.6\times10^{10}$ 、 $d=10$ で $9.8\times10^{16}$ 点。階層モデルではパラメータが何十〜何百次元になるので、格子の数値積分は原理的に不可能です。これが次元の呪い。格子の大半は事後密度がほぼゼロの無駄な領域でもあり、二重に効率が悪い。

3. モンテカルロ積分：標本平均・誤差は 1/√N・次元非依存

発想を変えます。「空間を埋め尽くす」のではなく、「事後から $N$ 個サンプル $\theta^{(1)},\dots,\theta^{(N)}$ を取り、期待値を標本平均で近似」。

\mathbb{E}[g(\theta)\mid D]\approx\frac{1}{N}\sum_{t=1}^{N}g(\theta^{(t)})

中心極限定理より、この推定の誤差は標本数 $N$ に対して $\propto 1/\sqrt{N}$ で縮み、次元 $d$ には依りません。確かめます。

import numpy as np
from scipy import stats

# (a) 既知の事後 Beta(16,8) で E[θ] を標本平均近似：誤差は ~1/√N
post = stats.beta(16, 8)
true_mean = post.mean()
rng = np.random.default_rng(0)
print(f"E[θ] のモンテカルロ近似（真値={true_mean:.4f}）:")
for N in [100, 1000, 10000, 100000]:
    est = post.rvs(N, random_state=rng).mean()
    print(f"  N={N:>6}: 推定={est:.4f}  誤差={abs(est-true_mean):.4f}  (1/√N={1/np.sqrt(N):.4f})")

# (b) 誤差は次元に依らない：高次元正規の第1成分の平均を同じ N で
print("\n次元を上げても誤差は同じ（標準正規の第1成分の平均, 真値=0, N=20000）:")
for d in [1, 10, 100]:
    s = rng.standard_normal((20000, d))
    print(f"  次元 d={d:>3}: 推定={s[:,0].mean():+.4f}  (~1/√N={1/np.sqrt(20000):.4f})")

出力：

E[θ] のモンテカルロ近似（真値=0.6667）:
  N=   100: 推定=0.6602  誤差=0.0065  (1/√N=0.1000)
  N=  1000: 推定=0.6739  誤差=0.0072  (1/√N=0.0316)
  N= 10000: 推定=0.6665  誤差=0.0001  (1/√N=0.0100)
  N=100000: 推定=0.6668  誤差=0.0001  (1/√N=0.0032)

次元を上げても誤差は同じ（標準正規の第1成分の平均, 真値=0, N=20000）:
  次元 d=  1: 推定=+0.0041  (~1/√N=0.0071)
  次元 d= 10: 推定=+0.0068  (~1/√N=0.0071)
  次元 d=100: 推定=-0.0033  (~1/√N=0.0071)

出力の意味：(a) サンプル数を増やすと推定誤差はおよそ $1/\sqrt{N}$ のスケールで縮みます（揺らぎはあるが $N$ が10倍で誤差はおよそ $1/\sqrt{10}$ ）。(b) 次元を $1\to100$ に上げても、第1成分の平均の推定誤差は同じ $\sim 1/\sqrt{N}$ のまま——モンテカルロは次元の呪いを受けません。格子が $G^d$ で爆発したのと対照的です。だから高次元の事後では、サンプリングが事実上唯一の道になります。

4. 残る問題：規格化なしの事後からどう引くか

モンテカルロは「事後からサンプルが取れる」前提でした。ところが手元にあるのは正規化されていない事後 $p(D\mid\theta)p(\theta)$ （分子）だけで、 $Z$ が分からない。scipy の .rvs() のように直接引く方法は、規格化された標準分布にしか使えません。

そこで「規格化されていない密度の値しか計算できなくても、その分布からサンプルを取る」工夫が要ります。それが MCMC——目標の事後を定常分布に持つマルコフ連鎖を作り、回して得た点列をサンプルとして使う（理論は MCMC（マルコフ連鎖モンテカルロ））。受容確率に事後が比でしか現れないため $Z$ が約分で消える、という仕掛けが核心です（メトロポリスヘイスティングス）。

flowchart LR
  Z["正規化定数 Z が解けない<br/>＋ 高次元で格子は爆発"] --> MC["モンテカルロ積分<br/>標本平均・誤差~1/√N・次元非依存"]
  MC --> Q["でも規格化なしの事後から<br/>どう引く?"]
  Q --> MCMC["MCMC：事後を定常分布にする<br/>マルコフ連鎖を回す（第4章）"]

⚠️ よくある誤解

「数値積分すればよい」ではない：低次元なら可。高次元は格子が $G^d$ で爆発し不可能（§2）。
「モンテカルロも次元の呪いを受ける」ではない：標本平均の誤差は $1/\sqrt{N}$ で次元非依存（§3）。これがサンプリングを選ぶ最大の理由。
「サンプルが取れれば $Z$ も分かる」ではない：MCMC は $Z$ を計算しません。期待値（比や割合）が欲しいだけなら $Z$ は要らない、というのがミソ（メトロポリスヘイスティングス）。
「 $1/\sqrt{N}$ は遅い」：確かに精度を10倍にするには $N$ を100倍。だが次元に依らない点が決定的で、高次元では他に手がありません（収束を速める HMC はハミルトニアンモンテカルロとNUTS）。