pymc-modeling

PyMC モデリング

最新の API パターンを使用した PyMC v5+ のベイズモデリングワークフローです。

ノートブックの推奨：プロジェクトですでに Jupyter を使用している場合を除き、対話型モデリングには marimo を使用してください。

モデルの指定

基本構造

import pymc as pm
import arviz as az

with pm.Model(coords=coords) as model:
    # Data containers (for out-of-sample prediction)
    x = pm.Data("x", x_obs, dims="obs")

    # Priors
    beta = pm.Normal("beta", mu=0, sigma=1, dims="features")
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Likelihood
    mu = pm.math.dot(x, beta)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=y_obs, dims="obs")

    # Inference
    idata = pm.sample()

coords と dims

モデルに意味のある構造がある場合は、解釈しやすい InferenceData を得るために coords/dims を使用してください。

coords = {
    "obs": np.arange(n_obs),
    "features": ["intercept", "age", "income"],
    "group": group_labels,
}

オーバーヘッドがメリットを上回る単純なモデルでは省略してください。

パラメーター化

データが少ない階層モデルでは、非中心化パラメーター化を優先してください。

# Non-centered (better for divergences)
offset = pm.Normal("offset", 0, 1, dims="group")
alpha = mu_alpha + sigma_alpha * offset

# Centered (better with strong data)
alpha = pm.Normal("alpha", mu_alpha, sigma_alpha, dims="group")

推論

デフォルトのサンプリング（nutpie）

デフォルトのサンプラーとして nutpie を使用してください。Rust ベースで、通常は 2～5 倍高速です。

with model:
    idata = pm.sample(
        draws=1000, tune=1000, chains=4,
        nuts_sampler="nutpie",
        random_seed=42,
    )

PyMC ネイティブのサンプリング

nutpie を利用できない場合は、PyMC の NUTS に切り替えてください。

with model:
    idata = pm.sample(draws=1000, tune=1000, chains=4, random_seed=42)

代替 MCMC バックエンド

references/inference.md を参照してください。

• NumPyro/JAX：GPU アクセラレーション、ベクトル化チェーン

近似推論

高速（ただし不正確）な事後近似には、次を使用します。

• ADVI/DADVI：ガウス近似を用いた変分推論
• Pathfinder：初期化またはスクリーニングのための準ニュートン最適化

診断と ArviZ ワークフロー

サンプリングの実行後は毎回、次の体系的なワークフローに従ってください。

フェーズ 1：即時チェック（必須）

# 1. Check for divergences (must be 0 or near 0)
n_div = idata.sample_stats["diverging"].sum().item()
print(f"Divergences: {n_div}")

# 2. Summary with convergence diagnostics
summary = az.summary(idata, var_names=["~offset"])  # exclude auxiliary
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "ess_bulk", "ess_tail", "r_hat"]])

# 3. Visual convergence check
az.plot_trace(idata, compact=True)
az.plot_rank(idata, var_names=["beta", "sigma"])

合格基準（次に進む前に、すべてを満たす必要があります）。

• 発散がゼロ（または 0.1% 未満で、ランダムに散在している）
• すべてのパラメーターで r_hat < 1.01
• ess_bulk > 400 かつ ess_tail > 400
• トレースプロットに良好な混合が見られる（密度が重なり、ぼやけたキャタピラー状になる）

フェーズ 2：詳細な収束チェック（フェーズ 1 の結果が境界的な場合）

# ESS evolution (should grow linearly)
az.plot_ess(idata, kind="evolution")

# Energy diagnostic (HMC health)
az.plot_energy(idata)

# Autocorrelation (should decay rapidly)
az.plot_autocorr(idata, var_names=["beta"])

フェーズ 3：モデルの批判的評価（必須）

# Generate posterior predictive
with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True)

# Does the model capture the data?
az.plot_ppc(idata, kind="cumulative")

# Calibration check
az.plot_loo_pit(idata, y="y")

重要なルール：フェーズ 1～3 に合格するまでは、パラメーターを解釈しないでください。

フェーズ 4：パラメーターの解釈

# Posterior summaries
az.plot_posterior(idata, var_names=["beta"], ref_val=0)

# Forest plots for hierarchical parameters
az.plot_forest(idata, var_names=["alpha"], combined=True)

# Parameter correlations (identify non-identifiability)
az.plot_pair(idata, var_names=["alpha", "beta", "sigma"])

ArviZ の包括的な使い方については、references/arviz.md を参照してください。 トラブルシューティングについては、references/diagnostics.md を参照してください。

事前予測チェックと事後予測チェック

事前予測（フィッティング前）

フィッティング前に、必ず事前分布の影響を確認してください。

with model:
    prior_pred = pm.sample_prior_predictive(draws=500)

# Do prior predictions span reasonable outcome range?
az.plot_ppc(prior_pred, group="prior", kind="cumulative")

# Numerical sanity check
prior_y = prior_pred.prior_predictive["y"].values.flatten()
print(f"Prior predictive range: [{prior_y.min():.1f}, {prior_y.max():.1f}]")

警告サイン：事前予測があり得ない値（負のカウント、1 を超える確率）を含む、または極端に広い／狭い。

事後予測（フィッティング後）

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

# Density comparison
az.plot_ppc(idata, kind="kde")

# Cumulative (better for systematic deviations)
az.plot_ppc(idata, kind="cumulative")

# Calibration diagnostic
az.plot_loo_pit(idata, y="y")

解釈：観測データ（濃い線）が事後予測分布（薄い線）の範囲内に収まっている必要があります。詳しい解釈については、references/arviz.md を参照してください。

モデルのデバッグ

モデル構造の確認

# Print model summary (variables, shapes, distributions)
print(model)

# Visualize model as directed graph
pm.model_to_graphviz(model)

モデルの指定ミスの確認

サンプリングの前に、モデルを検証してください。

# Debug model: checks for common issues
model.debug()

# Check initial point log-probabilities
# Identifies which variables have invalid starting values
model.point_logps()

よくある問題

症状	考えられる原因	対処法
ValueError: Shape mismatch	パラメーターと観測値の次元の不一致	インデックスベクトルを使用します：alpha[group_idx]
Initial evaluation failed	分布のサポート外のデータ	境界を確認し、init="adapt_diag" を使用します
Mass matrix contains zeros	スケーリングされていない予測変数、または平坦な事前分布	特徴量を標準化し、弱情報事前分布を使用します
発散の数が多い	ファンネル形状	非中心化パラメーター化を使用します
対数確率に NaN が含まれる	無効なパラメーターの組み合わせ	パラメーター制約を確認し、境界を追加します
対数確率が -inf	尤度のサポート外の観測値	データが分布の定義域と一致することを確認します
離散変数のサンプリングが遅い	NUTS は離散変数と互換性がない	離散変数を周辺化します
サンプル外予測の精度が低い	静的なデータ配列	mutable=True を指定した pm.Data コンテナーを使用します

問題と解決策の包括的なガイドについては、references/troubleshooting.md を参照してください。

発散のデバッグ

# Identify where divergences occur in parameter space
az.plot_pair(idata, var_names=["alpha", "beta", "sigma"], divergences=True)

# Check if divergences cluster in specific regions
# Clustering suggests parameterization or prior issues

遅いモデルのプロファイリング

# Time individual operations in the log-probability computation
profile = model.profile(model.logp())
profile.summary()

# Identify bottlenecks in gradient computation
import pytensor
grad_profile = model.profile(pytensor.grad(model.logp(), model.continuous_value_vars))
grad_profile.summary()

追加のトラブルシューティングについては、references/gotchas.md を参照してください。

モデル比較

LOO-CV（推奨）

# Compute LOO with pointwise diagnostics
loo = az.loo(idata, pointwise=True)
print(f"ELPD: {loo.elpd_loo:.1f} ± {loo.se:.1f}")

# Check Pareto k values (must be < 0.7 for reliable LOO)
print(f"Bad k (>0.7): {(loo.pareto_k > 0.7).sum().item()}")
az.plot_khat(idata)

モデルの比較

comparison = az.compare({
    "model_a": idata_a,
    "model_b": idata_b,
}, ic="loo")

print(comparison[["rank", "elpd_loo", "d_loo", "weight", "dse"]])
az.plot_compare(comparison)

判定ルール：d_loo < 2*dse の場合、モデルは実質的に同等です。

詳細なモデル比較のワークフローについては、references/arviz.md を参照してください。

結果の保存と読み込み

InferenceData の永続化

後で分析したり共有したりできるように、サンプリング結果を保存します。

# Save to NetCDF (recommended format)
idata.to_netcdf("results/model_v1.nc")

# Load
idata = az.from_netcdf("results/model_v1.nc")

圧縮ストレージ

大きな InferenceData オブジェクト（ドロー数が多い場合や、事後予測が大きい場合）には、次を使用します。

# Compress with zlib (reduces file size 50-80%)
idata.to_netcdf(
    "results/model_v1.nc",
    engine="h5netcdf",
    encoding={var: {"zlib": True, "complevel": 4}
              for group in ["posterior", "posterior_predictive"]
              if hasattr(idata, group)
              for var in getattr(idata, group).data_vars}
)

保存される内容

InferenceData は、ベイズワークフロー全体を保持します。

• posterior：MCMC によるパラメーターサンプル
• prior、prior_predictive：事前サンプル（生成した場合）
• posterior_predictive：予測値（生成した場合）
• observed_data、constant_data：フィッティングに使用したデータ
• sample_stats：診断情報（発散、木の深さ、エネルギー）
• log_likelihood：点ごとの対数尤度（LOO-CV 用）
• すべての座標と次元

ワークフローのパターン

# Save after each major step
with model:
    idata = pm.sample(nuts_sampler="nutpie")
idata.to_netcdf("results/step1_posterior.nc")

with model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True)
idata.to_netcdf("results/step2_with_ppc.nc")

# Resume later
idata = az.from_netcdf("results/step2_with_ppc.nc")
az.plot_ppc(idata)  # Continue analysis

事前分布の選択

以下については、references/priors.md を参照してください。

• 分布の種類ごとの弱情報事前分布のデフォルト
• 事前予測チェックのワークフロー
• 分野固有の推奨事項

一般的なパターン

階層型／マルチレベル

with pm.Model(coords={"group": groups, "obs": obs_idx}) as hierarchical:
    # Hyperpriors
    mu_alpha = pm.Normal("mu_alpha", 0, 1)
    sigma_alpha = pm.HalfNormal("sigma_alpha", 1)

    # Group-level (non-centered)
    alpha_offset = pm.Normal("alpha_offset", 0, 1, dims="group")
    alpha = pm.Deterministic("alpha", mu_alpha + sigma_alpha * alpha_offset, dims="group")

    # Likelihood
    y = pm.Normal("y", alpha[group_idx], sigma, observed=y_obs, dims="obs")

GLM

# Logistic regression
with pm.Model() as logistic:
    alpha = pm.Normal("alpha", 0, 2.5)  # intercept
    beta = pm.Normal("beta", 0, 2.5, dims="features")
    
    # Logit link
    logit_p = alpha + pm.math.dot(X, beta)
    p = pm.math.sigmoid(logit_p)
    
    y = pm.Bernoulli("y", p=p, observed=y_obs)

# Poisson regression
with pm.Model() as poisson:
    beta = pm.Normal("beta", 0, 1, dims="features")
    mu = pm.math.exp(pm.math.dot(X, beta))
    y = pm.Poisson("y", mu=mu, observed=y_obs)

ガウス過程

ほとんどの GP 問題（n > 500、1～3D 入力）では、HSGP をデフォルトとして使用してください。計算量は O(n³) ではなく O(nm) です。

with pm.Model() as gp_model:
    # Hyperparameters
    ell = pm.InverseGamma("ell", alpha=5, beta=5)
    eta = pm.HalfNormal("eta", sigma=2)
    sigma = pm.HalfNormal("sigma", sigma=0.5)

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

    # HSGP approximation
    gp = pm.gp.HSGP(m=[20], c=1.5, cov_func=cov)
    f = gp.prior("f", X=X[:, None])  # X must be 2D

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

周期的なパターンには pm.gp.HSGPPeriodic を使用してください。小規模なデータセット（n < 500）には、pm.gp.Marginal または pm.gp.Latent を使用してください。

以下については、references/gp.md を参照してください。

• HSGP パラメーターの選択（m と c の選択、自動ヒューリスティック）
• 季節性／周期性パターン向けの HSGPPeriodic
• 近似品質の診断
• 共分散関数と事前分布
• 一般的なパターン（トレンド＋季節性、分類、異分散）

時系列

with pm.Model(coords={"time": range(T)}) as ar_model:
    rho = pm.Uniform("rho", -1, 1)
    sigma = pm.HalfNormal("sigma", sigma=1)

    y = pm.AR("y", rho=[rho], sigma=sigma, constant=True,
              observed=y_obs, dims="time")

以下については、references/timeseries.md を参照してください。

• 自己回帰モデル（AR、ARMA）
• ランダムウォークモデルと局所水準モデル
• 構造的時系列（トレンド＋季節性）
• 状態空間モデル
• 時系列向け GP
• 複数の季節性への対応
• 予測パターン

BART（ベイズ加法回帰木）

import pymc_bart as pmb

with pm.Model() as bart_model:
    mu = pmb.BART("mu", X=X, Y=y, m=50)
    sigma = pm.HalfNormal("sigma", 1)
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

以下については、references/bart.md を参照してください。

• 回帰と分類
• 変数重要度と部分依存
• BART とパラメトリック成分の組み合わせ
• 設定（木の数、深さの事前分布）

混合モデル

import numpy as np

coords = {"component": range(K)}

with pm.Model(coords=coords) as gmm:
    # Mixture weights
    w = pm.Dirichlet("w", a=np.ones(K), dims="component")

    # Component parameters (with ordering to avoid label switching)
    mu = pm.Normal("mu", mu=0, sigma=10, dims="component",
                   transform=pm.distributions.transforms.ordered)
    sigma = pm.HalfNormal("sigma", sigma=2, dims="component")

    # Mixture likelihood
    y = pm.NormalMixture("y", w=w, mu=mu, sigma=sigma, observed=y_obs)

以下については、references/mixtures.md を参照してください。

• 有限混合モデルと回帰の混合モデル
• ラベルスイッチング問題と解決策（順序制約、ラベルの付け替え）
• 周辺化混合モデル（pymc-extras）
• 混合モデルの診断

特殊な尤度

# Zero-Inflated Poisson (excess zeros)
with pm.Model() as zip_model:
    psi = pm.Beta("psi", alpha=2, beta=2)  # P(structural zero)
    mu = pm.Exponential("mu", lam=1)
    y = pm.ZeroInflatedPoisson("y", psi=psi, mu=mu, observed=y_obs)

# Censored data (e.g., right-censored survival)
with pm.Model() as censored_model:
    mu = pm.Normal("mu", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=5)
    y = pm.Censored("y", dist=pm.Normal.dist(mu=mu, sigma=sigma),
                    lower=None, upper=censoring_time, observed=y_obs)

# Ordinal regression
with pm.Model() as ordinal:
    beta = pm.Normal("beta", mu=0, sigma=2, dims="features")
    cutpoints = pm.Normal("cutpoints", mu=0, sigma=2,
                          transform=pm.distributions.transforms.ordered,
                          shape=n_categories - 1)
    y = pm.OrderedLogistic("y", eta=pm.math.dot(X, beta),
                           cutpoints=cutpoints, observed=y_obs)

以下については、references/specialized_likelihoods.md を参照してください。

• ゼロ過剰モデル（ポアソン、負の二項、二項）
• カウントデータ向けハードルモデル
• 打ち切りデータと切断データ
• 順序回帰
• スチューデントの t 尤度によるロバスト回帰

よくある落とし穴

以下については、references/gotchas.md を参照してください。

• 中心化と非中心化のパラメーター化
• スケールパラメーターの事前分布
• 混合モデルのラベルスイッチング
• パフォーマンスの問題（GP、大規模な Deterministics）
• Python の条件分岐とハードクリッピング
• 階層モデルにおける冗長な切片

以下を網羅した包括的な問題解決ガイドについては、references/troubleshooting.md を参照してください。

• 形状と次元のエラー
• 初期化の失敗
• 質量行列と数値上の問題
• 離散変数に関する課題
• データコンテナーと予測に関する問題

因果推論の操作

pm.do（介入）

変数を固定値に設定する do 計算の介入を適用します。

with pm.Model() as causal_model:
    x = pm.Normal("x", 0, 1)
    y = pm.Normal("y", x, 1)
    z = pm.Normal("z", y, 1)

# Intervene: set x = 2 (breaks incoming edges to x)
with pm.do(causal_model, {"x": 2}) as intervention_model:
    idata = pm.sample_prior_predictive()
    # Samples from P(y, z | do(x=2))

pm.observe（条件付け）

介入せずに、観測値で条件付けします。

# Condition: observe y = 1 (doesn't break causal structure)
with pm.observe(causal_model, {"y": 1}) as conditioned_model:
    idata = pm.sample()
    # Samples from P(x, z | y=1)

do と observe の組み合わせ

# Intervention + observation for causal queries
with pm.do(causal_model, {"x": 2}) as m1:
    with pm.observe(m1, {"z": 0}) as m2:
        idata = pm.sample()
        # P(y | do(x=2), z=0)

pymc-extras

特殊なモデルと推論については、次を参照してください。

import pymc_extras as pmx

# Marginalize discrete parameters from a model
model = pmx.marginalize(model, ["discrete_var"])

# R2D2 prior for regression (requires output_sigma and input_sigma)
residual_sigma, beta = pmx.R2D2M2CP(
    "r2d2",
    output_sigma=y.std(),
    input_sigma=X.std(axis=0),
    dims="features",
    r2=0.5,
)

# Laplace approximation for fast inference
idata = pmx.fit_laplace(model)

カスタム分布とモデルコンポーネント

組み込みの分布を超えて PyMC を拡張する方法については、次をご覧ください。

import pymc as pm
import pytensor.tensor as pt

# Custom likelihood via DensityDist
def custom_logp(value, mu, sigma):
    return pm.logp(pm.Normal.dist(mu=mu, sigma=sigma), value)

with pm.Model() as model:
    mu = pm.Normal("mu", 0, 1)
    y = pm.DensityDist("y", mu, 1.0, logp=custom_logp, observed=y_obs)

# Soft constraints via Potential
with pm.Model() as model:
    alpha = pm.Normal("alpha", 0, 1, dims="group")
    pm.Potential("sum_to_zero", -100 * pt.sqr(alpha.sum()))

詳しくは、references/custom_models.mdをご覧ください。

• カスタム尤度にはpm.DensityDist
• ソフト制約とヤコビアン補正にはpm.Potential
• シミュレーションベースの推論（ABC）にはpm.Simulator
• カスタム事前分布にはpm.CustomDist