alterlab-pymc

PyMC ベイズモデリング

概要

PyMC は、ベイズモデリングと確率的プログラミングのための Python ライブラリです。階層モデル、MCMC サンプリング（NUTS）、変分推論、PSIS-LOO によるモデル比較など、現行 API（PyMC ≥ 6、ArviZ ≥ 1.1）を使用して、ベイズモデルの構築、適合、検証、比較を行えます。

uv pip install "pymc[nutpie]" "arviz[matplotlib,h5netcdf]"   # nutpie sampler, plotting backend, NetCDF I/O

PyMC 6 では、多くのチュートリアルで今も使われているいくつかのデフォルト値と名称が変更されています。PyMC 5／ArviZ 0.x のコードを再利用する前に、最後にあるバージョンに関する注記をお読みください。

このスキルを使用する場面

このスキルは、次の場合に使用してください。

• ベイズモデル（線形回帰／ロジスティック回帰、階層モデル、時系列など）を構築する
• MCMC サンプリングまたは変分推論を行う
• 事前予測チェック／事後予測チェックを行う
• サンプリングの問題（発散、収束、ESS）を診断する
• PSIS-LOO 交差検証を用いて複数のモデルを比較する
• ベイズ手法を用いて不確実性を定量化する
• 階層／マルチレベルのデータ構造を扱う
• 欠測データや測定誤差を適切に扱う

該当しない場面

シナリオ	代わりに使用
p 値、標準誤差、残差診断を用いる頻度論的回帰	alterlab-statsmodels
ICC と分散成分の報告を伴う、入れ子データ向けの報告重視の混合効果モデル（lme4/glmmTMB/brms、statsmodels MixedLM、または bambi の数式）	alterlab-multilevel-models
複数の研究間で効果量を統合する（固定効果／ランダム効果、異質性、ファンネルプロット）	alterlab-meta-analysis
点予測を行う機械学習（ランダムフォレスト、勾配ブースティング、CV チューニング）	alterlab-scikit-learn

標準的なベイズワークフロー

ベイズモデルを構築、検証する際は、次の 8 ステップのワークフローに従ってください。

1. データの準備 — 予測子を標準化し、欠測データを処理し、coords を設定します
2. モデルの構築 — 弱情報事前分布、名前付き dims、予測可能な変数には pm.Data() を使用します
3. 事前予測チェック — pm.sample_prior_predictive を実行し、適合前に事前分布を検証します
4. 適合 — pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9) を実行し、モデルを比較する場合は続けて pm.compute_log_likelihood(idata) を実行します
5. 診断 — R-hat < 1.01、ESS > 400、発散なし、良好なトレースの混合を確認します
6. 事後予測チェック — pm.sample_posterior_predictive を実行し、適合度と観測データを確認します
7. 分析 — az.summary、az.plot_dist（以前の plot_posterior）、az.plot_forest を使用します
8. 予測 — pm.set_data を実行してから pm.sample_posterior_predictive を実行し、HDI 区間を抽出します

ステップごとのコード全体：references/workflow_examples.md。

一般的なモデルパターン

PyMC は、線形／ロジスティック／ポアソン回帰、階層（マルチレベル）モデル、 時系列（AR）に対応しています。それぞれすぐに応用できるコードは references/model_patterns.md にあります。

グループごとの観測数が少ない場合、階層モデルには非中心化パラメータ化を使用してください。中心化形式ではファンネル形状の幾何構造が生じ、NUTS で探索しにくくなり、発散として現れます。 テンプレート：assets/linear_regression_template.py、assets/hierarchical_model_template.py。

分布の選択

事前分布と尤度の選択は、モデリングにおいて最も大きな影響をもたらす判断です。簡易選択ガイド （スケールパラメータ、非有界、正値、確率、相関行列、連続／カウント／二値／ カテゴリカル尤度）は references/distribution_selection.md にあります。包括的な一覧は references/distributions.md にあります。

モデルの比較とサンプリング

• モデルを比較するには、PSIS-LOO（scripts/model_comparison.py）を使用します。elpd_diffをdseと照らし合わせて解釈し、Pareto-kを確認します。ArviZ 1.xにはWAICがありません。
• サンプリングするには、デフォルトでNUTSを使用します（インストールされている場合はnutpie）。発散がある場合はtarget_acceptを引き上げます。高速な近似にはADVIを使用します。
• 診断するには、scripts/model_diagnostics.py（check_diagnostics、create_diagnostic_report）を使用します。
• 発散、低いESS、高いR-hat、サンプリングの遅さをトラブルシューティングします。

コードと判断基準については、references/model_comparison.mdを参照してください。サンプリングアルゴリズムの詳細なガイド： references/sampling_inference.md。

ベストプラクティス

モデルの構築

1. サンプリングを改善するため、予測変数は必ず標準化します
2. 弱情報事前分布を使用します（一様事前分布ではなく）
3. 明確さのため、名前付き次元（dims）を使用します
4. 階層モデルには非中心化パラメータ化を使用します
5. フィッティングの前に事前予測を確認します

サンプリング

1. 収束のため、複数のチェーン（少なくとも4つ）を実行します
2. target_accept=0.9を基準として使用します（必要に応じて引き上げます）
3. モデル比較の前に**pm.compute_log_likelihood(idata)を呼び出します**
4. 再現性のため、乱数シードを設定します

検証

1. 解釈の前に診断結果を確認します（R-hat、ESS、発散）
2. モデル検証のために事後予測チェックを行います
3. 適切な場合は複数のモデルを比較します
4. 不確実性を報告します（点推定だけでなく、HDI区間も報告します）

まずはシンプルなモデルから始め、複雑さを徐々に加えながら、各予測チェックに基づいてモデルを繰り返し調整してください（上記の8ステップのワークフローを参照）。

リソース

このスキルには次の内容が含まれます。

リファレンス（references/）

• workflow_examples.md：ベイズワークフローの8ステップについて、データ準備から予測までを網羅したコードを順を追って説明します。
• model_patterns.md：線形回帰、ロジスティック回帰、ポアソン回帰、階層モデル、AR時系列モデルにすぐに応用できるコードです。
• distribution_selection.md：パラメータや結果の種類に応じた事前分布と尤度のクイックガイドです。
• model_comparison.md：PSIS-LOOによる比較、診断スクリプト、トラブルシューティング（発散、ESS、R-hat、サンプリングの遅さ）について説明します。
• distributions.md：カテゴリ別（連続、離散、多変量、混合、時系列）に整理された、PyMCの分布の包括的なカタログです。事前分布や尤度を選ぶ際に使用します。
• sampling_inference.md：サンプリングアルゴリズム（NUTS、Metropolis、SMC）、変分推論（ADVI、SVGD）、サンプリングの問題への対処について詳しく説明します。収束の問題が発生した場合や、推論方法を選ぶ際に使用します。
• workflows.md：データ準備から保存までを網羅した、最初から最後まで実行可能な単一スクリプトに加え、workflow_examples.mdにはない追加のレシピ（欠損データの補完、QR再パラメータ化、混合モデル、モデル平均化ヘルパー）を掲載しています。

スクリプト（scripts/）

• model_diagnostics.py：診断チェックとレポート生成を自動化します（PyMC-NUTSとnutpieの両方のサンプル統計に対応）。関数：簡易チェック用のcheck_diagnostics()、プロットを含む包括的な分析用のcreate_diagnostic_report()。

• model_comparison.py：PSIS-LOOによるモデル比較ユーティリティです。関数：compare_models()、check_loo_reliability()、model_averaging()。

テンプレート（assets/）

• linear_regression_template.py：データ準備、事前分布のチェック、フィッティング、診断、予測までを網羅した、ベイズ線形回帰の完全なテンプレートです。

• hierarchical_model_template.py：非中心化パラメータ化とグループレベル分析を含む、階層モデル／マルチレベルモデルの完全なテンプレートです。

クイックリファレンス

モデルの構築

with pm.Model(coords={'var': names}) as model:
    # Priors
    param = pm.Normal('param', mu=0, sigma=1, dims='var')
    # Likelihood
    y = pm.Normal('y', mu=..., sigma=..., observed=data)

サンプリング

idata = pm.sample(draws=2000, tune=1000, chains=4, target_accept=0.9)  # returns xarray.DataTree
pm.compute_log_likelihood(idata)  # needed for LOO / az.compare

診断

from scripts.model_diagnostics import check_diagnostics
check_diagnostics(idata)

モデルの比較

from scripts.model_comparison import compare_models
compare_models({'m1': idata1, 'm2': idata2})  # PSIS-LOO; higher elpd is better

予測

# X must have been wrapped at build time: pm.Data('X', X, dims=('obs', 'predictors'))
with model:
    pm.set_data({'X': X_new}, coords={'obs': range(len(X_new))})
    pm.sample_posterior_predictive(idata, predictions=True, extend_inferencedata=True)
# predictions land in idata.predictions

追加の注意事項

• PyMC は、可視化と診断に ArviZ を統合しています
• pm.model_to_graphviz(model) を使うと、モデルの構造を可視化できます（graphviz パッケージが必要です）
• idata.to_netcdf('results.nc') で結果を保存できます（h5netcdf または netCDF4 が必要です。例：arviz[h5netcdf]）。az.from_netcdf('results.nc') で読み込みます
• 非常に大規模なモデルでは、ミニバッチ ADVI またはデータのサブサンプリングを検討してください

バージョンに関する注意事項（PyMC 6 / ArviZ 1.x）

PyMC 6.0（2026 年 5 月）では、ArviZ 1.x と PyTensor 3 に移行しました。PyMC 5 / ArviZ 0.x 向けに書かれたコードは、以下の点で動作しなくなります。

• 結果は az.InferenceData ではなく、xarray.DataTree です。グループへのアクセス（idata.posterior、idata["posterior"]）は引き続き可能です。idata.extend(other) は idata.update(other) に変更されました。グループの一覧は list(idata.children) で取得できます。
• 対数尤度：サンプリング後に pm.compute_log_likelihood(idata) を呼び出してください。pm.sample(idata_kwargs={"log_likelihood": True}) は非推奨です。
• サンプラー：nutpie がインストールされている場合、デフォルトの NUTS サンプラーになります（デフォルトのチューニングステップ数は 400 です。PyMC 独自の NUTS では 1000 です）。nuts_sampler="pymc" | "nutpie" | "numpyro" | "blackjax" を使って明示的に選択してください。サンプル統計の名前（tree_depth と depth）は異なるため、診断コードではどちらが存在するか確認してください。初期値は initvals= に指定します（start= は失敗します）。
• バックエンド：PyTensor 3 はデフォルトで Numba を使ってコンパイルするため、初回の呼び出しには JIT コンパイル時間が含まれます。
• 要約：az.summary は、94% HDI（hdi_3%、hdi_97%）の代わりに 89% 等裾区間（eti89_lb、eti89_ub）を報告します。az.summary(idata, ci_kind="hdi", ci_prob=0.94) を使うか、az.rcParams["stats.ci_kind"] / ["stats.ci_prob"] を設定してください。az.hdi(..., prob=0.95) は hdi_prob= に代わる引数です。
• プロットは PlotCollection を返します（pc.savefig("file.png") で保存し、pc.show() で表示します）。ax= / axes= 引数は廃止されました。名前の変更：plot_posterior → plot_dist、plot_ppc → plot_ppc_dist（事前予測チェックでは group="prior_predictive" を指定し、num_pp_samples= の代わりに num_samples= を使います）、plot_trace → plot_trace_dist（トレースと密度）、plot_dist_comparison → plot_prior_posterior。バックエンドは arviz[matplotlib] でインストールしてください。
• モデル比較：az.waic は廃止され、az.compare(dict) は ic= / scale= を受け付けなくなりました。スタッキング重みを使って PSIS-LOO elpd の高い順に並べ替え（値が高いほど良好）、elpd_diff、dse、p_worse、diag_* 列を報告します。
• 新しいデータで予測するには、モデルの構築時に予測子を pm.Data('X', X, dims=...) でラップする必要があります。その場合に限り、pm.set_data({'X': X_new}, coords={...}) で値を入れ替えられます。グラフに直接組み込まれた通常の NumPy 配列は置き換えられません。
• pm.sample_prior_predictive では draws= を指定します（以前の samples= キーワードは削除されました）。
• 標本外予測を行うには、pm.sample_posterior_predictive(idata, predictions=True, extend_inferencedata=True, ...) を呼び出してください。結果は idata.posterior_predictive ではなく idata.predictions に格納されます。