モデル記述から推論まで

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

📎 前提：確率的プログラミング概観　|　関連：ベイズ線形回帰

⚠️ 要最新確認：PyMC 6.0.1 / ArviZ 1.2.0 で検証。sample_prior_predictive・sample_posterior_predictive の引数は版で変わります。

要点（BLUF）

ベイズ分析の標準ワークフローは 事前 → 事前予測チェック → サンプリング → 事後 → 事後予測チェック → 要約。PPL はこの全工程を支えます。
PyMC で線形回帰を宣言（係数に正規事前、尺度に半正規、正規尤度）。pm.sample で事後、sample_posterior_predictive で予測を得ます。
第7章で手実装したベイズ線形回帰（ベイズ線形回帰）が宣言的に書け、真の $a,b,\sigma$ を回収。事後予測の 95% 帯がデータをほぼ覆うことを確かめます。

1. ベイズワークフロー全体

良いベイズ分析は、いきなり結果を出さず、段階を踏んで点検します。

flowchart LR
  M["モデル記述<br/>事前 + 尤度"] --> PR["事前予測チェック<br/>含意が常識的か（第3章）"]
  PR --> S["サンプリング<br/>NUTS で事後を得る"]
  S --> D["収束診断<br/>R-hat・ESS（第4章）"]
  D --> PO["事後の要約<br/>平均・区間・確率（第3章）"]
  PO --> PP["事後予測チェック<br/>当てはまり診断（第8章）"]

第3〜8章で個別に学んだ道具（事前予測・収束診断・事後要約・事後予測チェック）が、PPL では一つの流れにつながります。

2. PyMC で線形回帰を一気通貫に書く

データ $y=a+bx+\varepsilon$ を、PyMC で宣言して推論します。

import numpy as np, pymc as pm, arviz as az

rng = np.random.default_rng(1)
a_true, b_true, sig = 1.0, -2.0, 0.5
x = np.sort(rng.uniform(-2, 2, 30))
y = a_true + b_true*x + rng.normal(0, sig, len(x))

with pm.Model() as model:
    a = pm.Normal("a", 0, 5)                            # 弱情報事前
    b = pm.Normal("b", 0, 5)
    sigma = pm.HalfNormal("sigma", 2)                  # 尺度は正の半正規
    mu = a + b*x                                        # 線形予測子
    pm.Normal("y", mu=mu, sigma=sigma, observed=y)     # 尤度（observed でデータを結ぶ）

    prior = pm.sample_prior_predictive(200, random_seed=0)        # ① 事前予測
    idata = pm.sample(2000, tune=1000, chains=4, cores=1,         # ② サンプリング
                      random_seed=0, progressbar=False)
    pm.sample_posterior_predictive(idata, random_seed=0,         # ③ 事後予測
                                   progressbar=False, extend_inferencedata=True)

for nm, tv in [("a", a_true), ("b", b_true), ("sigma", sig)]:    # ④ 事後の要約
    d = idata.posterior[nm].values.flatten()
    lo, hi = np.percentile(d, [5.5, 94.5])
    print(f"{nm}: 事後平均={d.mean():+.3f} SD={d.std():.3f} 89%=[{lo:+.3f},{hi:+.3f}]  (真={tv})")
print(f"r_hat(b)={float(az.rhat(idata)['b'].values):.3f}  ESS(b)={float(az.ess(idata)['b'].values):.0f}")

pp = idata.posterior_predictive["y"].values.reshape(-1, len(x))  # ⑤ 事後予測チェック
cover = np.mean((y >= np.percentile(pp,2.5,0)) & (y <= np.percentile(pp,97.5,0)))
print(f"事後予測95%帯がデータを覆う割合={cover:.2f}")

出力：

a: 事後平均=+0.988 SD=0.091 89%=[+0.841,+1.132]  (真=1.0)
b: 事後平均=-2.032 SD=0.078 89%=[-2.155,-1.907]  (真=-2.0)
sigma: 事後平均=+0.484 SD=0.069 89%=[+0.388,+0.603]  (真=0.5)
r_hat(b)=1.000  ESS(b)=8580
事後予測95%帯がデータを覆う割合=0.97

出力の意味：宣言したモデルから、真の係数 $a=1.0,b=-2.0,\sigma=0.5$ をすべて 89% 区間内に回収。 $\hat R=1.000$ ・ESS $8580$ で収束も良好。事後予測の 95% 帯がデータ点の 97% を覆い（理想は 95%）、モデルがデータをうまく説明できていると分かります。第7章で行列計算で書いた事後（ベイズ線形回帰）が、ここではモデルを宣言するだけで得られました。

3. 各ステップの役割

① 事前予測チェック：sample_prior_predictive で、データを見る前にモデルが含意する $y$ の範囲を確認（事前の選び方と感度分析）。事前が非常識に広い/狭いと、ここで気づけます。
② サンプリング：NUTS が事後からサンプル。idata（InferenceData）に格納。
③ 事後予測：sample_posterior_predictive で、事後パラメータから複製データを生成（信用区間と事後予測分布）。
④ 事後の要約：平均・区間・確率（事後分布の要約）。収束（ $\hat R$ ・ESS）を必ず併記。
⑤ 事後予測チェック：複製データと観測を比べてモデルの欠陥を探す（事後予測チェック）。

この5段が、PPL を使うときの基本動作です。モデル（事前＋尤度）を差し替えれば、GLM・階層・GP でも同じ流れが回ります。

4. 手実装との対応

第7章のベイズ線形回帰は、計画行列・事後共分散 $S_N$ ・事後平均 $m_N$ を手で計算しました。PyMC では a = pm.Normal(...)・mu = a + b*x・pm.Normal("y", mu, sigma, observed=y) と生成過程を書くだけ。共役でない尤度（ロジスティック・ポアソン、ベイズGLM）でも、尤度の行を変えるだけで同じワークフローが使えます——ここが手実装にない PPL の強みです。

⚠️ よくある誤解

「サンプリングすれば終わり」ではない：事前予測（前）と事後予測（後）の点検が両端に要ります（§3）。
「事後予測がデータを覆えば完璧」ではない：覆う割合は1つの指標。複数の検定統計量で見ます（事後予測チェック）。
「決定的な変数（mu）も確率変数」ではない：mu = a + b*x は確率変数の関数（決定的）。observed= を付けた行だけが尤度（データと結ぶ）です。
「半正規でなく正規を尺度に使ってよい」ではない： $\sigma>0$ の制約があるので HalfNormal 等を使います。正規だと負を許して破綻。