pymc

PyMC によるベイズモデリング

概要

PyMC は、ベイズモデリングと確率的プログラミングのための Python ライブラリです。PyMC の最新 API（バージョン 6.x 以降）を使って、階層モデル、MCMC サンプリング（NUTS）、変分推論、事後予測チェック、モデル比較（LOO、WAIC）などのベイズモデルを構築、適合、検証、比較できます。

現行バージョンとセットアップ

2026 年 6 月時点で、PyMC 6.0.1 が最新の安定版です。Python 3.12 以降が必要で、計算グラフのバックエンドとして PyTensor 3 を使用し、デフォルトでは Numba などのコンパイル済みバックエンドを使用します。再現可能なローカル環境を用意するには、バージョンを固定します。

uv pip install "pymc[nutpie]==6.0.1"

nutpie エクストラを使うと、より高速な Rust/Numba NUTS 実装を利用できます。NumPyro または BlackJAX を使う場合は、同じ環境にそれらのオプションのサンプラー依存パッケージをインストールし、プロジェクトのロックファイルにもバージョンを固定してください。

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

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

• ベイズモデル（線形回帰／ロジスティック回帰、階層モデル、時系列など）を構築する
• MCMC サンプリングまたは変分推論を行う
• 事前予測チェック／事後予測チェックを実施する
• サンプリングの問題（発散、収束、ESS）を診断する
• 情報量規準（LOO、WAIC）を使って複数のモデルを比較する
• ベイズ手法によって不確実性を定量化する
• 階層型／マルチレベルのデータ構造を扱う
• 欠損データや測定誤差を、原則に沿った方法で扱う

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

まずサンプリングして後から確認するのは避けてください。references/standard_workflow.md にコードとともに記載されている 8 段階のワークフローは次のとおりです。

1. データの準備 — 事前分布を解釈しやすくするための予測変数の標準化などを含みます。
2. モデルの構築 — pm.Model コンテキスト内で事前分布と尤度を定義します。
3. 事前予測チェック — モデルを適合する前に、事前分布がもっともらしいデータを示唆することを確認します。
4. モデルの適合 — 明示的なシードを指定して pm.sample() を実行します。
5. 診断の確認 — R-hat、ESS、発散を確認します。発散があると適合結果は無効です。target_accept を上げて期待するのではなく、モデルを修正するか、再パラメーター化してください。
6. 事後予測チェック — 適合したモデルが観測データを再現するか確認します。
7. 結果の分析 — 事後分布から要約と区間を求めます。
8. 予測の実施 — pm.set_data と事後予測サンプリングを使って新しいデータを予測します。

再利用可能なモデル構造とモデル比較については、 references/model_patterns.md を参照してください。

分布選択ガイド

事前分布

スケールパラメーター（σ、τ）：

• pm.HalfNormal('sigma', sigma=1) - デフォルトの選択肢
• pm.Exponential('sigma', lam=1) - 代替案
• pm.Gamma('sigma', alpha=2, beta=1) - より情報量の多い選択肢

有界でないパラメーター：

• pm.Normal('theta', mu=0, sigma=1) - 標準化済みデータ用
• pm.StudentT('theta', nu=3, mu=0, sigma=1) - 外れ値に対してロバスト

正のパラメーター：

• pm.LogNormal('theta', mu=0, sigma=1)
• pm.Gamma('theta', alpha=2, beta=1)

確率：

• pm.Beta('p', alpha=2, beta=2) - 弱情報事前分布
• pm.Uniform('p', lower=0, upper=1) - 無情報事前分布（使用は控えめに）

相関行列：

• pm.LKJCholeskyCov('chol', n=n_vars, eta=2, sd_dist=pm.HalfNormal.dist(1)) - 推奨される共分散事前分布
• pm.LKJCorr('corr', n=n_vars, eta=2) - 相関のみの事前分布。eta=1 は一様、eta>1 は単位行列を優先

尤度

連続アウトカム：

• pm.Normal('y', mu=mu, sigma=sigma) - 連続データのデフォルト
• pm.StudentT('y', nu=nu, mu=mu, sigma=sigma) - 外れ値に対して頑健

カウントデータ：

• pm.Poisson('y', mu=lambda) - 分散が平均に等しいカウントデータ
• pm.NegativeBinomial('y', mu=mu, alpha=alpha) - 過分散のカウントデータ
• pm.ZeroInflatedPoisson('y', psi=psi, mu=mu) - ゼロが過剰なデータ
• pm.HurdleNegativeBinomial('y', psi=psi, mu=mu, alpha=alpha) - ゼロが過剰で、過分散もあるデータ

二値アウトカム：

• pm.Bernoulli('y', p=p) または pm.Bernoulli('y', logit_p=logit_p)

カテゴリカルアウトカム：

• pm.Categorical('y', p=probs)

参照先： 包括的な分布リファレンスについては references/distributions.md を参照してください。

サンプリングと推論

NUTS による MCMC

ほとんどのモデルでデフォルトとして推奨されています。

idata = pm.sample(
    draws=2000,
    tune=1000,
    chains=4,
    target_accept=0.9,
    random_seed=42
)

必要に応じて調整します：

• 発散が発生する場合 → target_accept=0.95 以上に設定します
• サンプリングが遅い場合 → 初期化に ADVI を使用します
• 離散パラメーターがある場合 → 離散変数には pm.Metropolis() を使用します

変分推論

探索や初期化に適した高速な近似です。

with model:
    approx = pm.fit(n=20000, method='advi')

    # Use for initialization
    initvals = approx.sample(return_inferencedata=False)[0]
    idata = pm.sample(initvals=initvals)

トレードオフ：

• MCMC より大幅に高速です
• 近似であるため、不確実性を過小評価することがあります
• 大規模モデルや迅速な探索に適しています

参照先： サンプリングの詳細なガイドについては references/sampling_inference.md を参照してください。

診断スクリプト

包括的な診断

from scripts.model_diagnostics import create_diagnostic_report

create_diagnostic_report(
    idata,
    var_names=['alpha', 'beta', 'sigma'],
    output_dir='diagnostics/'
)

次の項目を作成します：

• トレースプロット
• ランクプロット（混合の確認）
• 自己相関プロット
• エネルギープロット
• 局所 ESS プロット
• 要約統計量の CSV

簡易診断チェック

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata)

R-hat、ESS、発散、ツリー深度を確認します。

よくある問題と解決策

発散

症状： idata.sample_stats.diverging.sum() > 0

解決策：

1. target_accept=0.95 または 0.99 に引き上げます
2. 非中心化パラメーター化を使用します（階層モデル）
3. パラメーターを制約するために、より強い事前分布を追加します
4. モデルの誤指定がないか確認します

有効サンプルサイズが小さい

症状： ESS < 400

解決策：

1. サンプル数を増やします：draws=5000
2. 事後相関を減らすために再パラメーター化します
3. 相関のある予測変数を用いた回帰には QR 分解を使用します

R-hat が大きい

症状： R-hat > 1.01

解決策：

1. チェーンを長く実行します：tune=2000, draws=5000
2. 多峰性がないか確認します
3. ADVI を使用して初期化を改善します

サンプリングが遅い

解決策：

1. ADVI による初期化を使用します
2. モデルの複雑さを軽減します
3. 並列化を増やします：cores=8, chains=8
4. 適切な場合は変分推論を使用します

ベストプラクティス

モデル構築

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

サンプリング

1. 収束を確認するため、複数のチェーン（少なくとも 4 本）を実行します
2. 基準として target_accept=0.9 を使用します（必要に応じて引き上げます）
3. モデル比較のために log_likelihood=True を含めます
4. 再現性のために乱数シードを設定します

検証

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

ワークフロー

1. シンプルなモデルから始め、徐々に複雑さを加えます
2. 事前予測チェック → フィッティング → 診断 → 事後予測チェック
3. チェックの結果に基づいてモデルの仕様を反復的に見直します
4. 仮定と事前分布の選択を記録します

リソース

このスキルには次のものが含まれます：

References (references/)

• distributions.md：カテゴリ別（連続、離散、多変量、混合、時系列）に整理された、PyMC の分布の包括的なカタログです。事前分布や尤度を選ぶ際に使用します。

• sampling_inference.md：サンプリングアルゴリズム（NUTS、Metropolis、SMC）、変分推論（ADVI、SVGD）、サンプリングの問題への対処方法を詳しく解説します。収束の問題が発生した場合や、推論方法を選ぶ際に使用します。

• workflows.md：一般的なモデルの種類、データの準備、事前分布の選択、モデルの検証に関する、完全なワークフローの例とコードパターンです。標準的なベイズ解析のレシピ集として使用します。

Scripts (scripts/)

• model_diagnostics.py：診断チェックとレポート生成を自動化します。関数：簡易チェック用の check_diagnostics()、プロットを含む包括的な分析用の create_diagnostic_report()。

• model_comparison.py：ArviZ 1.x の compare() が順位付けに使用する唯一の基準である PSIS-LOO ELPD に基づく、モデル比較ユーティリティです。関数：compare_models()、check_loo_reliability()、model_averaging()。

Templates (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)

診断

from scripts.model_diagnostics import check_diagnostics
check_diagnostics(idata)

モデルの比較

from scripts.model_comparison import compare_models
compare_models({'m1': idata1, 'm2': idata2}, ic='loo')

予測

with model:
    pm.set_data({'X_data': X_new})
    pred = pm.sample_posterior_predictive(idata, predictions=True)

補足事項

• PyMC は可視化と診断のために ArviZ と連携します。PyMC 6／ArviZ 1 では、.posterior や .posterior_predictive などの使い慣れたグループを維持しながら、xarray の DataTree を使用します。
• pm.model_to_graphviz(model) を使ってモデル構造を可視化します。
• idata.to_netcdf('results.nc') で結果を保存します。
• az.from_netcdf('results.nc') で読み込みます。
• 非常に大規模なモデルでは、ミニバッチ ADVI またはデータのサブサンプリングを検討してください。

エージェントの作業手順

1. 環境を確認します。 Python 環境とライブラリのバージョン（python -c "import pkg; print(pkg.__version__)"）を確認し、データの形状、型、欠損値を調べます。
2. 入力を明確にします。 データまたはユーザーから、形式、識別子、パラメータを確認します。結果に影響する値は推測せず、質問してください。
3. まず小規模に実行します。 最初にサンプルまたは単一のフォールドで実行し、実行時間とメモリを確認します。
4. 上記の指示と参考資料を使って、タスク全体を実行します。
5. 結果を検証します。 ホールドアウトデータ、固定した乱数シード、適切な指標を使用し、データリークがないか確認します。不確実性（CI、シード間の標準偏差）を報告します。
6. 報告します。 実行内容（バージョン、コマンド、パラメータ）、確認内容、未確認事項を明記します。

発生した状況	対応
メモリ不足または実行が非常に遅い	サブサンプリング、チャンク処理または遅延計算を使うか、モデルの規模を縮小し、変更内容をユーザーに伝えます。
これらの指示に記載された関数、フラグ、エンドポイントがインストール済みバージョンに存在しない	インストール済みバージョンの公式ドキュメント（help()、--help、公式ドキュメント）を確認し、適宜対応して、その旨をユーザーに伝えます。API を捏造してはいけません。
必須の入力、識別子、パラメータが曖昧	ユーザーに質問するか、実行前に仮定を明示します。

誠実性に関するルール

• 結果、パラメータ、識別子、引用、統計値を捏造してはいけません。実行または検証できないことがあれば、明確に伝えます。
• このセッションで計算していない指標を報告してはいけません。各数値を生成したコードパスを示します。
• ここに記載されたバージョン固有の詳細は古くなっている可能性があるため、インストール済みバージョンの公式ドキュメントで確認します。
• 費用が発生する操作、共有 GPU やクラウドクォータを消費する操作、個人データや患者データに触れる操作、取り消せない操作を行う前に、ユーザーに確認します。

関連スキル

• statistical-analysis：研究データの統計分析をガイドします。検定の選択、仮定の確認、効果量、検出力分析、ベイズ的な代替手法などに対応します…