pymc-bayesian-modeling

PyMC ベイズモデリング

概要

PyMC は、ベイズモデリングと確率的プログラミングのための Python ライブラリです。PyMC のモダンな API（バージョン 5.x 以降）を使用して、ベイズモデルの構築、フィッティング、検証、比較を行えます。階層モデル、MCMC サンプリング（NUTS）、変分推論、モデル比較（LOO、WAIC）に対応しています。

このスキルを使用するタイミング

このスキルは次の場合に使用します：

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

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

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

1. データ準備

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

# Load and prepare data
X = ...  # Predictors
y = ...  # Outcomes

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

主要なプラクティス：

• 連続値の説明変数を標準化する（サンプリング効率が向上します）
• 可能な場合は目的変数を中心化する
• 欠損データは明示的に扱う（パラメータとして扱う）
• 分かりやすさのために coords で名前付きの次元を使用する

2. モデルの構築

coords = {
    'predictors': ['var1', 'var2', 'var3'],
    'obs_id': np.arange(len(y))
}

with pm.Model(coords=coords) as model:
    # Priors
    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')

主要なプラクティス：

• 弱情報事前分布を使用する（無情報の一様事前分布ではなく）
• 尺度パラメータには HalfNormal または Exponential を使用する
• 可能な場合は shape の代わりに名前付きの次元（dims）を使用する
• 予測のために値を更新する場合は pm.Data() を使用する

3. 事前予測チェック

フィッティングの前に、必ず事前分布を検証してください：

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

# Visualize
az.plot_ppc(prior_pred, group='prior')

確認事項：

• 事前予測は妥当な値の範囲をカバーしているか？
• ドメイン知識に照らして極端な値は妥当か？
• 事前分布が非現実的なデータを生成する場合は、調整して再確認する

4. モデルのフィッティング

with model:
    # Optional: Quick exploration with ADVI
    # approx = pm.fit(n=20000)

    # Full MCMC inference
    idata = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        target_accept=0.9,
        random_seed=42,
        idata_kwargs={'log_likelihood': True}  # For model comparison
    )

主要なパラメータ：

• draws=2000：チェーンごとのサンプル数
• tune=1000：ウォームアップのサンプル（破棄されます）
• chains=4：収束を確認するために 4 本のチェーンを実行する
• target_accept=0.9：難しい事後分布の場合はより高く（0.95〜0.99）
• モデル比較のために log_likelihood=True を含める

5. 診断の確認

診断スクリプトを使用します：

from scripts.model_diagnostics import check_diagnostics

results = check_diagnostics(idata, var_names=['alpha', 'beta', 'sigma'])

確認事項：

• R-hat < 1.01：チェーンが収束している
• ESS > 400：有効サンプル数が十分である
• ダイバージェンスなし：NUTS が正常にサンプリングできた
• トレースプロット：チェーンがよく混ざっている（ふさふさした毛虫のような形状）

問題が発生した場合：

• ダイバージェンス → target_accept=0.95 に上げ、非中心化パラメータ化を使用する
• ESS が低い → サンプル数を増やす、相関を減らすために再パラメータ化する
• R-hat が高い → より長く実行する、多峰性を確認する

6. 事後予測チェック

モデルの適合度を検証します：

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

# Visualize
az.plot_ppc(idata)

確認事項：

• 事後予測は観測データのパターンを捉えているか？
• 系統的な乖離が見られるか（モデルの誤指定）？
• 適合が悪い場合は代替モデルを検討する

7. 結果の分析

# Summary statistics
print(az.summary(idata, var_names=['alpha', 'beta', 'sigma']))

# Posterior distributions
az.plot_posterior(idata, var_names=['alpha', 'beta', 'sigma'])

# Coefficient estimates
az.plot_forest(idata, var_names=['beta'], combined=True)

8. 予測の実行

X_new = ...  # New predictor values
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
    )

# Extract prediction intervals
y_pred_mean = post_pred.posterior_predictive['y_obs'].mean(dim=['chain', 'draw'])
y_pred_hdi = az.hdi(post_pred.posterior_predictive, var_names=['y_obs'])

一般的なモデルのパターン

線形回帰

連続値の目的変数と線形の関係がある場合：

with pm.Model() as linear_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)
    sigma = pm.HalfNormal('sigma', sigma=1)

    mu = alpha + pm.math.dot(X, beta)
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

テンプレートの使用： assets/linear_regression_template.py

ロジスティック回帰

二値の目的変数の場合：

with pm.Model() as logistic_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

    logit_p = alpha + pm.math.dot(X, beta)
    y = pm.Bernoulli('y', logit_p=logit_p, observed=y_obs)

階層モデル

グループ化されたデータの場合（非中心化パラメータ化を使用）：

with pm.Model(coords={'groups': group_names}) as hierarchical_model:
    # Hyperpriors
    mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=10)
    sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=1)

    # Group-level (non-centered)
    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
    mu = alpha[group_idx]
    sigma = pm.HalfNormal('sigma', sigma=1)
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=y_obs)

テンプレートの使用： assets/hierarchical_model_template.py

重要： ダイバージェンスを避けるため、階層モデルでは必ず非中心化パラメータ化を使用してください。

ポアソン回帰

カウントデータの場合：

with pm.Model() as poisson_model:
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)

    log_lambda = alpha + pm.math.dot(X, beta)
    y = pm.Poisson('y', mu=pm.math.exp(log_lambda), observed=y_obs)

過分散のカウントデータには、代わりに NegativeBinomial を使用してください。

時系列

自己回帰過程の場合：

with pm.Model() as ar_model:
    sigma = pm.HalfNormal('sigma', sigma=1)
    rho = pm.Normal('rho', mu=0, sigma=0.5, shape=ar_order)
    init_dist = pm.Normal.dist(mu=0, sigma=sigma)

    y = pm.AR('y', rho=rho, sigma=sigma, init_dist=init_dist, observed=y_obs)

モデル比較

モデルの比較

モデル比較には LOO または WAIC を使用します：

from scripts.model_comparison import compare_models, check_loo_reliability

# Fit models with log_likelihood
models = {
    'Model1': idata1,
    'Model2': idata2,
    'Model3': idata3
}

# Compare using LOO
comparison = compare_models(models, ic='loo')

# Check reliability
check_loo_reliability(models)

解釈：

• Δloo < 2：モデルはほぼ同等であり、より単純なモデルを選択する
• 2 < Δloo < 4：より良いモデルへの根拠は弱い
• 4 < Δloo < 10：中程度の根拠
• Δloo > 10：より良いモデルへの強い根拠

Pareto-k 値の確認：

• k < 0.7：LOO は信頼できる
• k > 0.7：WAIC または k 分割交差検証を検討する

Model Averaging

モデルが互いに似ている場合は、予測を平均します:

from scripts.model_comparison import model_averaging

averaged_pred, weights = model_averaging(models, var_name='y_obs')

分布の選択ガイド

事前分布

尺度パラメータ (σ, τ):

• 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.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.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
    start = approx.sample(return_inferencedata=False)[0]
    idata = pm.sample(start=start)

トレードオフ:

• 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/)

• distributions.md: PyMC の分布を、連続、離散、多変量、混合、時系列のカテゴリ別にまとめた包括的なカタログです。事前分布や尤度を選択する際に使用します。

• sampling_inference.md: サンプリングアルゴリズム（NUTS、Metropolis、SMC）、変分推論（ADVI、SVGD）、およびサンプリングの問題への対処についての詳細なガイドです。収束の問題に遭遇したとき、または推論手法を選択するときに使用します。

• workflows.md: 一般的なモデルタイプ、データ準備、事前分布の選択、モデル検証のための、完全なワークフロー例とコードパターンです。標準的なベイズ分析のクックブックとして使用してください。

スクリプト (scripts/)

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

• model_comparison.py: LOO/WAIC を用いたモデル比較のユーティリティです。関数: compare_models()、check_loo_reliability()、model_averaging()。

テンプレート (assets/)

• linear_regression_template.py: データ準備、事前チェック、フィッティング、診断、予測までを含む、ベイズ線形回帰の完全なワークフローのテンプレートです。

• hierarchical_model_template.py: 非中心化パラメータ化と group レベルの分析を用いた、階層モデル/多階層モデルの完全なテンプレートです。

クイックリファレンス

モデル構築

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': X_new})
    pred = pm.sample_posterior_predictive(idata.posterior)

補足事項

• PyMC は、可視化と診断のために ArviZ と統合されています
• モデルの構造を可視化するには pm.model_to_graphviz(model) を使用します
• 結果の保存には idata.to_netcdf('results.nc') を使用します
• 読み込みには az.from_netcdf('results.nc') を使用します
• 非常に大規模なモデルの場合は、ミニバッチ ADVI やデータのサブサンプリングを検討してください