pymc-bayesian-modeling

PyMC ベイズモデリング

概要

PyMC は、ベイズ統計モデリングと確率的プログラミングのための Python ライブラリです。確率モデルを定義するための表現力豊かな構文と、MCMC（NUTS）および変分法（ADVI）による効率的な推論を提供します。このスキルでは、モデルの指定から診断、比較、予測まで、ベイズモデリングのサイクル全体を扱います。

使用する場面

• 点推定だけでなく、信用区間などで不確実性を完全に定量化してパラメーターを推定する場合
• グループ化データや入れ子構造のデータに階層モデル／マルチレベルモデルを適合させる場合
• モデルの仮定を検証するために、事前予測チェックや事後予測チェックを行う場合
• 情報量規準（LOO-CV、WAIC）を使って候補モデルを比較する場合
• ベイズフレームワークで回帰モデル（線形、ロジスティック、ポアソン）を構築する場合
• 欠損データや測定誤差を潜在パラメーターとして扱う場合
• 自己回帰事前分布やランダムウォーク事前分布を使って時系列をモデル化する場合
• 不確実性の範囲を含む、新しい観測値の事後予測を生成する場合
• 大規模モデルでコンパイル済みの、よりスケーラブルなベイズ推論を行う場合は、代わりに Stan/PyStan を使用します。頻度論的な統計検定には statsmodels を使用します

前提条件

• Python パッケージ：pymc >= 5.0、arviz、numpy、matplotlib
• データ：数値列を含む NumPy 配列または pandas DataFrame
• 環境：ほとんどのモデルでは CPU で十分です。大規模モデルでは JAX バックエンド経由で GPU を利用できます

pip install pymc arviz numpy matplotlib
# Optional: JAX backend for GPU acceleration
pip install pymc[jax]

クイックスタート

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

# Simulate data
np.random.seed(42)
X = np.random.randn(100)
y = 2.5 + 1.3 * X + np.random.randn(100) * 0.5

# Build and fit model
with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0, sigma=5)
    beta = pm.Normal("beta", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=1)
    mu = alpha + beta * X
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)
    idata = pm.sample(1000, tune=1000, chains=4, random_seed=42)

print(az.summary(idata, var_names=["alpha", "beta", "sigma"]))
# Expected: alpha ~ 2.5, beta ~ 1.3, sigma ~ 0.5

ワークフロー

ステップ 1：データを準備する

サンプリング効率を高めるため、連続予測変数を標準化します。読みやすいモデルと ArviZ との統合のために、名前付き座標を使用します。

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

# Load data
X = np.random.randn(200, 3)  # 200 obs, 3 predictors
y = X @ np.array([1.0, -0.5, 0.3]) + np.random.randn(200) * 0.8

# Standardize predictors
X_mean, X_std = X.mean(axis=0), X.std(axis=0)
X_scaled = (X - X_mean) / X_std

# Define coordinates for named dimensions
coords = {
    "predictors": ["var1", "var2", "var3"],
    "obs_id": np.arange(len(y)),
}
print(f"Data shape: X={X_scaled.shape}, y={y.shape}")

ステップ 2：モデルを定義し、事前分布を設定する

pm.Model() コンテキスト内でモデル構造を指定します。弱情報事前分布、名前付き次元を指定する dims、およびスケールパラメーターには HalfNormal または Exponential を使用します。

with pm.Model(coords=coords) as model:
    # Priors — weakly informative, not flat
    alpha = pm.Normal("alpha", mu=0, sigma=1)
    beta = pm.Normal("beta", mu=0, sigma=1, dims="predictors")
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Linear predictor
    mu = alpha + pm.math.dot(X_scaled, beta)

    # Likelihood
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y, dims="obs_id")

# Inspect model variables
print(model.basic_RVs)  # Lists: [alpha, beta, sigma, y_obs]

ステップ 3：事前予測チェック

適合前に、事前分布から妥当な範囲のデータが生成されることを検証します。シミュレーションデータが不合理な場合は、事前分布を調整します。

with model:
    prior_pred = pm.sample_prior_predictive(samples=1000, random_seed=42)

# Check prior-implied data range
prior_y = prior_pred.prior_predictive["y_obs"].values.flatten()
print(f"Prior predictive range: [{prior_y.min():.1f}, {prior_y.max():.1f}]")
print(f"Observed data range:    [{y.min():.1f}, {y.max():.1f}]")

az.plot_ppc(prior_pred, group="prior", num_pp_samples=100)

ステップ 4：事後分布をサンプリングする（MCMC）

複数のチェーンを使って NUTS サンプリングを実行します。後でモデルを比較する予定がある場合は、log_likelihood=True を指定します。

with model:
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={"log_likelihood": True},
    )

print(f"Posterior shape: {idata.posterior['beta'].shape}")
# Expected: (4 chains, 2000 draws, 3 predictors)

ステップ 5：サンプリングを診断する

結果を解釈する前に、収束を確認します。3 つの診断項目（R-hat、ESS、発散）がすべて基準を満たす必要があります。

# Summary with convergence diagnostics
summary = az.summary(idata, var_names=["alpha", "beta", "sigma"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "r_hat", "ess_bulk"]])

# R-hat convergence check
bad_rhat = summary[summary["r_hat"] > 1.01]
if len(bad_rhat) > 0:
    print(f"WARNING: {len(bad_rhat)} parameters with R-hat > 1.01")
    print(bad_rhat[["r_hat"]])

# Effective sample size check
low_ess = summary[summary["ess_bulk"] < 400]
if len(low_ess) > 0:
    print(f"WARNING: {len(low_ess)} parameters with ESS < 400")

# Divergence check
n_div = idata.sample_stats.diverging.sum().item()
total = len(idata.posterior.draw) * len(idata.posterior.chain)
print(f"Divergences: {n_div}/{total} ({n_div / total * 100:.2f}%)")

# Visual diagnostics — trace plots and rank plots
az.plot_trace(idata, var_names=["alpha", "beta", "sigma"])
az.plot_rank(idata, var_names=["alpha", "beta", "sigma"])

ステップ 6：事後予測チェック

事後分布からシミュレーションしたデータと観測データを比較して、モデルの適合度を検証します。

with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True, random_seed=42)

az.plot_ppc(idata, num_pp_samples=100)
# Blue = observed data, grey = posterior simulations
# Systematic deviations indicate model misspecification

ステップ 7：モデルを比較する

LOO-CV または WAIC を使って候補モデルを比較します。情報量規準は値が小さいほど良好です。

# Fit multiple models with log_likelihood=True, then compare
# Example: compare linear vs a second model
idatas = {"linear": idata}  # add more fitted models here

comparison = az.compare(idatas, ic="loo")
print(comparison[["rank", "elpd_loo", "p_loo", "d_loo", "weight"]])

# Check LOO reliability via Pareto-k diagnostics
loo_result = az.loo(idata, pointwise=True)
high_k = (loo_result.pareto_k > 0.7).sum().item()
print(f"Observations with Pareto-k > 0.7: {high_k}")
# Interpretation: Dloo < 2 = similar models; Dloo > 10 = strong evidence

az.plot_compare(comparison)

ステップ 8：予測を生成する

不確実性を完全に伝播させて、新しいデータの事後予測を生成します。

X_new = np.array([[0.5, -1.0, 0.2]])
X_new_scaled = (X_new - X_mean) / X_std

with model:
    pm.set_data({"X_scaled": X_new_scaled})
    post_pred = pm.sample_posterior_predictive(
        idata.posterior, var_names=["y_obs"], random_seed=42
    )

y_pred = post_pred.posterior_predictive["y_obs"]
print(f"Predicted mean: {y_pred.mean().item():.3f}")
print(f"94% HDI: {az.hdi(y_pred, hdi_prob=0.94).values}")

主なパラメーター

パラメーター	デフォルト	範囲／オプション	効果
draws	1000	500-10000	チェーンごとの事後サンプル数
tune	1000	500-5000	ウォームアップ反復数（破棄されます）。複雑な事後分布では増やします
chains	4	2-8	独立したチェーンの数。信頼できる R-hat を得るには最低 4 が必要です
cores	すべての CPU	1-N	並列チェーン数。完全な並列処理には chains と同じ値を設定します
target_accept	0.8	0.8-0.99	NUTS の受容率。発散を減らすには値を上げます
init	"auto"	"adapt_diag"、"jitter+adapt_diag"、"advi"	サンプラーの初期化方法
random_seed	None	任意の int	再現性のためのシード
idata_kwargs	{}	{"log_likelihood": True}	LOO／WAIC によるモデル比較のために対数尤度を保存します
method (pm.fit)	"advi"	"advi"、"fullrank_advi"、"svgd"	変分推論アルゴリズム
n (pm.fit)	10000	5000-100000	VI の最適化反復数
samples (prior pred)	500	100-5000	検証用の事前予測サンプル数

主な概念

事前分布の選択ガイド

分布	使用する場面	主なパラメーター
Normal(mu, sigma)	有界でない実数値パラメーター（標準化データ）	mu：中心、sigma：ばらつき
HalfNormal(sigma)	スケール／標準偏差パラメーター（正の値）	sigma：正の半分のばらつき
Exponential(lam)	スケールパラメーター。HalfNormal の代替	lam：率（1／平均）
StudentT(nu, mu, sigma)	Normal のロバストな代替（外れ値に強い）	nu：自由度（10 未満では裾がより重い）
Beta(alpha, beta)	[0,1] の確率または比率	alpha=beta=2：弱情報事前分布
Gamma(alpha, beta)	正のパラメーター（率、濃度）	alpha：形状、beta：率
LogNormal(mu, sigma)	乗法効果を伴う正のパラメーター	mu、sigma：基礎となる Normal の値
LKJCorr(n, eta)	相関行列の事前分布	eta=1：一様、eta>1：単位行列を優先
Dirichlet(a)	確率ベクトル（合計が 1）	a：濃度。すべて同じ場合は一様
Bernoulli(p / logit_p)	二値結果の尤度	数値安定性のため logit_p を使用
Poisson(mu)	カウントデータ（等分散）	mu：率。過分散の場合は NegBinomial を使用
NegativeBinomial(mu, alpha)	過分散カウントデータ	alpha：分散（小さいほど過分散が大きい）

診断のしきい値

指標	しきい値	解釈	基準を満たさない場合の対処
R-hat	< 1.01	チェーンが収束している	チェーンを長く実行する。多峰性を確認する
ESS bulk	> 400	独立サンプルが十分にある	draws を増やす。再パラメーター化する
ESS tail	> 400	裾の推定値が信頼できる	draws を増やす
発散	0	NUTS による探索が正常に行われた	target_accept を増やす。非中心化パラメーター化を行う
Pareto-k (LOO)	< 0.7	LOO 推定値が信頼できる	WAIC または k 分割交差検証を使用
最大木の深さ	< 10	軌跡が途中で打ち切られていない	再パラメーター化するか、max_treedepth を増やす

モデルのバリエーション概要

問題の種類	レシピ	尤度	主な特徴
グループ化／入れ子データ	階層モデル	Normal（変動あり）	非中心化パラメーター化、部分プーリング
二値結果	ロジスティック回帰	Bernoulli	logit_p リンク関数
非線形／空間データ	ガウス過程	Normal	カーネルベースの共分散、柔軟な形状
カウントデータ	（ワークフローで Poisson を使用）	Poisson / NegBinomial	対数リンク。過分散には NegBinomial
時系列	（参考資料を参照）	AR / GaussianRandomWalk	自己回帰係数
混合／クラスタリング	（参考資料を参照）	Mixture / NormalMixture	Dirichlet による成分の重み

よく使うレシピ

レシピ：階層モデル

使用する場面：データに自然なグループ構造がある場合（病院内の患者、学校内の生徒など）。非中心化パラメーター化により、ファネル形状に起因する発散を回避します。

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

n_groups, n_per_group = 5, 30
group_idx = np.repeat(np.arange(n_groups), n_per_group)
group_names = [f"group_{i}" for i in range(n_groups)]

# Simulated grouped data
true_alphas = np.random.normal(3.0, 1.5, n_groups)
y_obs = np.random.normal(true_alphas[group_idx], 0.5)

with pm.Model(coords={"groups": group_names}) as hierarchical_model:
    # Hyperpriors (population level)
    mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=5)
    sigma_alpha = pm.HalfNormal("sigma_alpha", sigma=2)

    # Non-centered parameterization
    alpha_offset = pm.Normal("alpha_offset", mu=0, sigma=1, dims="groups")
    alpha = pm.Deterministic(
        "alpha", mu_alpha + sigma_alpha * alpha_offset, dims="groups"
    )

    # Observation level
    sigma = pm.HalfNormal("sigma", sigma=1)
    y = pm.Normal("y", mu=alpha[group_idx], sigma=sigma, observed=y_obs)

    idata_hier = pm.sample(2000, tune=1000, target_accept=0.95, random_seed=42)

print(az.summary(idata_hier, var_names=["mu_alpha", "sigma_alpha", "alpha", "sigma"]))

レシピ：ロジスティック回帰

使用する場面：二値の結果変数（成功／失敗、病気／健康など）。確率を直接計算する代わりに、数値安定性のため logit_p を使用します。

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

# Simulated binary outcome
np.random.seed(42)
n, p = 200, 3
X_lr = np.random.randn(n, p)
true_beta = np.array([1.0, -0.5, 0.3])
prob = 1 / (1 + np.exp(-(0.2 + X_lr @ true_beta)))
y_binary = np.random.binomial(1, prob)

with pm.Model() as logistic_model:
    alpha = pm.Normal("alpha", mu=0, sigma=2)
    beta = pm.Normal("beta", mu=0, sigma=2, shape=p)
    logit_p = alpha + pm.math.dot(X_lr, beta)
    y = pm.Bernoulli("y", logit_p=logit_p, observed=y_binary)

    idata_logit = pm.sample(2000, tune=1000, target_accept=0.9, random_seed=42)

summary = az.summary(idata_logit, var_names=["alpha", "beta"])
print(summary[["mean", "sd", "hdi_3%", "hdi_97%"]])
# Odds ratios
beta_samples = idata_logit.posterior["beta"].values.reshape(-1, p)
print(f"Odds ratios (median): {np.median(np.exp(beta_samples), axis=0)}")

レシピ：ガウス過程

使用する場面：未知の非線形関数、空間データ、または滑らかな潜在過程をモデル化する場合。n > 1000 では計算コストが高いため、より大きなデータセットではスパース近似を検討してください。

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

# Simulated nonlinear data
np.random.seed(42)
X_gp = np.sort(np.random.uniform(0, 10, 60))[:, None]
y_gp = np.sin(X_gp[:, 0]) + np.random.randn(60) * 0.3

with pm.Model() as gp_model:
    # GP hyperparameters
    ls = pm.Gamma("lengthscale", alpha=2, beta=1)
    eta = pm.HalfNormal("amplitude", sigma=2)
    sigma_noise = pm.HalfNormal("noise", sigma=0.5)

    # Covariance function
    cov = eta**2 * pm.gp.cov.Matern52(1, ls=ls)
    gp = pm.gp.Marginal(cov_func=cov)

    # Marginal likelihood
    y_ = gp.marginal_likelihood("y", X=X_gp, y=y_gp, sigma=sigma_noise)
    idata_gp = pm.sample(1000, tune=1000, target_accept=0.9, random_seed=42)

# Predict on new points
X_new_gp = np.linspace(0, 10, 100)[:, None]
with gp_model:
    f_pred = gp.conditional("f_pred", X_new_gp)
    pred_samples = pm.sample_posterior_predictive(idata_gp, var_names=["f_pred"])

print(f"GP predictions shape: {pred_samples.posterior_predictive['f_pred'].shape}")

期待される出力

• トレースプロット（az.plot_trace）――チェーンの混合状態を可視化します。良好なチェーンは「ふわふわした毛虫」のように見えます
• 事後分布の要約（az.summary）――パラメーターごとの平均、SD、HDI、R-hat、ESS
• 事前・事後予測プロット（az.plot_ppc）――シミュレーションデータと観測データの分布
• フォレストプロット（az.plot_forest）――信用区間付きの係数推定値
• モデル比較テーブル（az.compare）――LOO/WAIC、重み、警告を含むモデルの順位
• ランクプロット（az.plot_rank）――一様なヒストグラムは良好な混合を示します
• エネルギープロット（az.plot_energy）――HMC のエネルギー遷移の診断
• InferenceData ファイル（idata.to_netcdf("results.nc")）――後で分析できるようにシリアライズした結果

トラブルシューティング

問題	原因	解決策
発散遷移	NUTS にとって事後分布の幾何構造が難しい（ファンネル、リッジ）	target_accept=0.95 に増やす。非中心化パラメーター化を使う。より強い事前分布を設定する
ESS が低い（< 400）	サンプル間の自己相関が高い	draws=5000 に増やす。相関を減らすよう再パラメーター化する。相関のある予測変数には QR 分解を試す
R-hat > 1.01	チェーンが収束していない	tune=2000, draws=5000 に増やす。多峰性を確認する。ADVI で初期化する
サンプリングが遅い	事後分布が複雑、またはデータセットが大きい	初期化に ADVI を使う。モデルの複雑さを減らす。cores を増やす。JAX バックエンドを試す
LOO で Pareto-k > 0.7	LOO 推定値に影響する観測値がある	代わりに WAIC を使う。影響の大きい点を調べる。k-fold CV を検討する
最大ツリー深度に到達	モデルの幾何構造に長い軌道が必要	モデルを再パラメーター化する。max_treedepth パラメーターを増やす
開始時に SamplingError が発生する	初期化が不適切	init="adapt_diag" または init="jitter+adapt_diag" を試す。MAP 経由で initvals を指定する
事後分布に偏りがある（事前分布が支配的）	データに対して事前分布が強すぎる	事前分布を弱める（sigma を大きくする）。事前予測を確認する。データが正しく渡されていることを確認する
ValueError: logp = -inf	パラメーターが不可能な領域に達した	データに NaN/Inf がないか確認する。正のパラメーターには HalfNormal/Gamma を使う。尤度がデータ型と一致していることを確認する

ベストプラクティス

1. フラットな事前分布ではなく、弱情報事前分布を使ってください。 フラットな事前分布（広い範囲の Uniform）はサンプリングを難しくし、適切であることはほとんどありません。
2. フィッティングの前に、必ず事前予測チェックを実行してください。 事前分布から得られるデータがもっともらしくない場合は、サンプリング前に事前分布を調整してください。
3. 連続予測変数を平均=0、sd=1 に標準化すると、NUTS の幾何構造が改善し、サンプリングが速くなります。
4. 階層モデルには非中心化パラメーター化を使ってください。 param ~ Normal(mu, sigma) の代わりに、offset ~ Normal(0,1); param = mu + sigma * offset とします。
5. 信頼性の高い R-hat と ESS の診断のため、少なくとも 4 つのチェーンを実行してください。 単一のチェーンを決して信用しないでください。
6. 事後分布を解釈する前に、3 つの診断すべてを確認してください（R-hat < 1.01、ESS > 400、発散 0）。
7. LOO または WAIC でモデルを比較する予定がある場合は、フィッティング時に log_likelihood=True を含めてください。
8. 名前付き次元には shape ではなく dims を使ってください。 ラベル付きの要約やサブセット化で ArviZ と連携できます。
9. 点推定値だけでなく、信用区間（HDI）を報告してください。 結果は事後分布であり、単一の数値ではありません。
10. 再現性の確保と後での再分析のため、idata.to_netcdf("results.nc") で結果を保存してください。

同梱リソース

オンデマンドで詳しく確認できるよう、2 つのリファレンスファイルを用意しています。

• references/distributions_inference.md -- カテゴリ別（連続、離散、多変量、混合、時系列、特殊／修飾）の分布カタログです。すべてのパラメーターシグネチャ、サポート状況、一般的な用途を記載しています。サンプリングと推論の手法（NUTS、Metropolis、Slice、CompoundStep、SMC、ADVI、fullrank ADVI、SVGD、MAP）については、コード例と手法選択ガイドの表を掲載しています。再パラメーター化のテクニック（非中心化、QR 分解）についてもコードを記載しています。 内容：元の references/distributions.md（320 行）と references/sampling_inference.md（424 行）を統合したものです。 本文から移動した内容：分布選択ガイドの基本表（主要概念）、診断しきい値の表（主要概念）、ADVI の基本的な使い方（元のレシピ「変分推論」。こちらのリファレンスで解説）。 省略した内容：形状ブロードキャストの例（単純な NumPy 形式の shape=5）、基本的な dims の使用方法（ワークフローのステップ 2 で解説）、事前予測サンプリング（ワークフローのステップ 3 で解説）。

• references/advanced_workflows.md -- モデル比較のワークフロー（LOO、WAIC、Pareto-k の信頼性チェック、重み付き予測によるモデル平均化）、包括的な診断レポートの生成（トレースプロット、ランクプロット、自己相関、エネルギー、ESS の推移）、事前分布と事後分布の比較パターン、データ準備のベストプラクティス（標準化、中心化、欠損データの補完）、事前分布の選択ガイドライン（弱情報事前分布と、ドメイン知識に基づく情報事前分布）、名前付き次元（dims）の使用方法と xarray によるサブセット化、保存／読み込みパターン（NetCDF、pickle）を記載しています。混合モデルのパターンも含みます。 内容：元の references/workflows.md（526 行）を統合したものです。 本文から移動した内容：診断チェックの基本ロジック（ワークフローのステップ 5）、モデル比較の基本（ワークフローのステップ 7）、事前予測チェック（ワークフローのステップ 3）、線形／ロジスティック／階層モデルのコード（レシピ）。 省略した内容：包括的なモノリシックワークフローのテンプレート（8 ステップのワークフローセクションと重複するため）。

リファレンスファイルごとの処理（元ファイルすべて）

元の references/（3 ファイル）：

1. distributions.md（320 行） -> sampling_inference.md と統合し、新しい references/distributions_inference.md に集約
2. sampling_inference.md（424 行） -> distributions.md と統合し、新しい references/distributions_inference.md に集約
3. workflows.md（526 行） -> 新しい references/advanced_workflows.md に移行

元の scripts/（2 ファイル）：

1. model_diagnostics.py（350 行） -> check_diagnostics() のロジック（R-hat、ESS、発散、ツリー深度のチェック）はワークフローのステップ 5 に記載。create_diagnostic_report() のプロット生成パターンは「期待される出力」で説明し、references/advanced_workflows.md を参照。compare_prior_posterior() ユーティリティは references/advanced_workflows.md で解説
2. model_comparison.py（387 行） -> compare_models() と check_loo_reliability() のパターンは、az.compare と az.loo を直接使ってワークフローのステップ 7 に記載。model_averaging() のパターンと cross_validation_comparison() のガイダンスは references/advanced_workflows.md で解説。plot_model_comparison() は az.plot_compare の薄いラッパーのため、ステップ 7 に記載

関連スキル

• arviz（予定） -- 事後分布の可視化と診断。PyMC は ArviZ InferenceData オブジェクトを返します
• bambi（予定） -- PyMC の数式ベースのインターフェース（R スタイルの y ~ x1 + x2 構文）。標準的な GLM に適した、よりシンプルな API を提供します
• numpyro（予定） -- JAX ベースの確率的プログラミング。ハードウェアアクセラレーションによる、より高速な NUTS を提供します
• statsmodels-statistical-modeling -- 頻度論的回帰と GLM。ベイズフレームワークが不要な場合に使用します
• scikit-learn-machine-learning -- ML による分類／回帰。不確実性の定量化よりも予測精度が重要な場合に使用します

参考資料

• PyMC documentation -- 公式 API リファレンスとチュートリアル
• PyMC GitHub -- ソースコードと課題トラッカー
• ArviZ documentation -- 事後分布の分析と可視化
• Bayesian Modeling and Computation in Python -- PyMC を使用した教科書
• PyMC examples gallery -- 厳選されたノートブック集