← ベイズ統計 一覧

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

📎 前提：ベイズファクターとモデル平均化　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計）

⚠️ 要最新確認：PPL の API は速く変わります。本ノートは PyMC 6.0.1 / ArviZ 1.2.0 で実行検証しました。バージョンが違うと関数名・既定値・出力列が変わるので、必ず手元の版で確認してください。

要点（BLUF）

• **確率的プログラミング言語（PPL）**は、生成モデルを宣言的に書けば、推論（NUTS など）をエンジンが自動で行う仕組み。代表は PyMC・Stan・NumPyro。
• モデルを書く＝事前と尤度を宣言するだけ。pm.sample が NUTS（ハミルトニアンモンテカルロとNUTS）で事後サンプルを返します。第4章で手実装した MCMC の自動化です。
• 実装で、PyMC の事後が第2章で手計算した閉形式 $\mathrm{Beta}(28,14)$ と一致することを確かめます。

1. PPL とは：宣言的モデリング

第4章では MH やギブスを numpy で手実装し、第7章ではラプラス近似を自分で書きました。PPL はこれを自動化します。やることは「データがどう生成されたか（事前＋尤度）を宣言する」だけ。あとは：

• 自動微分で対数事後の勾配を計算し、
• NUTS（自動調整 HMC）で事後からサンプリングし、
• 収束診断（R-hat・ESS）まで出してくれる。

推論アルゴリズムを書かずにモデルだけ書けるので、複雑な階層モデル・GLM も宣言的に組めます。

2. 主要 PPL の位置づけ

どれも「モデル宣言 → NUTS で推論 → ArviZ で診断・可視化」という流れは共通（ArviZ は推論結果を扱う共通ライブラリ）。以下は PyMC で示しますが、考え方は移植できます（具体的 API は要最新確認）。

3. コード：最小の PyMC モデルが閉形式と一致する

コイン（ベルヌーイ）の成功確率 $\theta$ を推定します。事前 $\mathrm{Beta}(1,1)$・尤度ベルヌーイなら、第2章の共役で事後は閉形式 $\mathrm{Beta}(1+k,1+n-k)$（ベータ二項モデル）。PyMC の NUTS がこれを再現するか確かめます。

import numpy as np, pymc as pm, arviz as az
from scipy import stats

rng = np.random.default_rng(0)
y = rng.binomial(1, 0.7, size=40); n, k = len(y), int(y.sum())

with pm.Model() as model:
    theta = pm.Beta("theta", alpha=1, beta=1)          # 事前 Beta(1,1)
    pm.Bernoulli("obs", p=theta, observed=y)           # 尤度（データを observed で渡す）
    idata = pm.sample(2000, tune=1000, chains=4, cores=1,   # ← Windowsは cores=1 推奨
                      random_seed=0, progressbar=False)

draws = idata.posterior["theta"].values.flatten()      # 事後サンプルを取り出す
lo, hi = np.percentile(draws, [5.5, 94.5])             # 89%区間
rhat = float(az.rhat(idata)["theta"].values)
ess  = float(az.ess(idata)["theta"].values)
post = stats.beta(1+k, 1+n-k)                          # 閉形式の事後
print(f"データ: n={n}, k={k}")
print(f"PyMC : 平均={draws.mean():.4f} SD={draws.std():.4f} 89%=[{lo:.3f},{hi:.3f}] r_hat={rhat:.3f} ESS={ess:.0f}")
print(f"閉形式: 平均={post.mean():.4f} SD={post.std():.4f} 89%=[{post.ppf(.055):.3f},{post.ppf(.945):.3f}] Beta({1+k},{1+n-k})")

出力：

データ: n=40, k=27
PyMC : 平均=0.6654 SD=0.0722 89%=[0.546,0.776] r_hat=1.000 ESS=3396
閉形式: 平均=0.6667 SD=0.0719 89%=[0.548,0.777] Beta(28,14)

出力の意味：PyMC の事後（平均 $0.665$・89%区間 $[0.546,0.776]$）が、手計算の閉形式 $\mathrm{Beta}(28,14)$（平均 $0.667$・$[0.548,0.777]$）と小数第2位までぴたり一致。$\hat R=1.000$・ESS $3396$ で収束も良好。with pm.Model() の中で事前と尤度を2行宣言しただけで、NUTS が事後を正しくサンプリングしました。第4章で苦労して書いた MCMC が、ここでは自動です。

4. 何が自動化されたのか

• 勾配：対数事後の微分を自動微分で計算（HMC/NUTS に必須、ハミルトニアンモンテカルロとNUTS）。手で導関数を書かなくてよい。
• サンプラー調整：NUTS が軌道長・ステップ幅を自動調整（手調整不要）。
• 診断：idata に複数チェーンが入り、az.rhat・az.ess で収束を即チェック（収束診断）。

宣言と推論が分離されるので、モデルを差し替えるだけで階層モデル・GLM・GP まで同じ流れで扱えます（モデルの書き方から事後予測までは モデル記述から推論まで）。

5. 実行上の注意（要最新確認）

• Windows の並列サンプリング：cores>1 のマルチプロセスは if __name__ == "__main__": ガードが要ります。スクリプトでの手軽さなら cores=1（チェーンを順番に回す）が安全。
• ArviZ の出力列はバージョン依存：本ノートの ArviZ 1.2.0 では区間が eti89（89%等裾）。az.summary の数値は文字列で返るので、計算には idata.posterior[...].values から取り出す。
• PPL は更新が速い：関数名・既定値が変わります。エラーが出たら、まず手元のバージョンの最新ドキュメントを確認。

⚠️ よくある誤解

• 「PPL は中で魔法をしている」ではない：第4章の MCMC・第6章の VI を自動化しただけ。中身は同じ数理です。
• 「宣言すれば必ず正しく推論できる」ではない：漏斗などで NUTS も詰まります（診断とトラブルシュート）。診断は必須。
• 「PyMC のコードはそのまま動く」とは限らない：バージョン差・OS 差（Windows の cores）でつまずきます。要最新確認。
• 「サンプル平均＝真の事後」ではない：収束（$\hat R\approx1$）と十分な ESS を確認してから信じます。

関連ノート

• ベイズファクターとモデル平均化
• モデル記述から推論まで（次のトピック・回帰モデルの完全な流れ）
• ハミルトニアンモンテカルロとNUTS（PPL が使う NUTS の中身）
• ベータ二項モデル（PyMC が再現した閉形式）
• 第9章 確率的プログラミング 目次
• ベイズ統計テキスト 全体目次