scientific-bayesian-statistics

Scientific Bayesian Statistics

An analysis skill for Bayesian statistical inference. Performs parameter estimation, model comparison, and decision support leveraging probabilistic programming frameworks (PyMC, Stan, NumPyro).

When to Use

• Posterior distribution estimation of parameters
• Hierarchical Bayesian models (multilevel, mixed effects)
• Bayesian regression and Bayesian classification
• MCMC diagnostics (Rhat, ESS, trace plots)
• Model comparison (WAIC, LOO-CV, Bayes factors)
• Bayesian optimization (Gaussian Process)
• Posterior predictive checks

Quick Start

Bayesian analysis pipeline

Phase 1: Model Specification
  - 尤度関数の選択 (Normal, Poisson, Binomial, etc.)
  - 事前分布の設定 (弱情報事前分布推奨)
  - パラメータ空間の定義
    ↓
Phase 2: Prior Predictive Check
  - 事前予測シミュレーション
  - 事前分布の妥当性検証
  - ドメイン知識との整合性確認
    ↓
Phase 3: MCMC Sampling
  - NUTS (No-U-Turn Sampler) 実行
  - チェーン数・サンプル数設定 (chains≥4, draws≥2000)
  - 収束診断 (Rhat < 1.01, ESS > 400)
    ↓
Phase 4: Posterior Analysis
  - 事後分布可視化 (forest, trace, pair plot)
  - 95% HDI (Highest Density Interval)
  - 事後予測チェック (PPC)
    ↓
Phase 5: Model Comparison
  - WAIC / LOO-CV 算出
  - Bayes Factor 計算
  - モデル重み付け平均 (Stacking)
    ↓
Phase 6: Reporting
  - 推定値サマリーテーブル
  - 事後分布プロット
  - 結果の解釈・レポート生成

Workflow

1. PyMC: Bayesian linear regression

import pymc as pm
import arviz as az
import numpy as np
import matplotlib.pyplot as plt

# === データ生成 (例) ===
np.random.seed(42)
N = 100
X = np.random.randn(N)
true_alpha, true_beta, true_sigma = 1.0, 2.5, 0.5
y = true_alpha + true_beta * X + np.random.normal(0, true_sigma, N)

# === ベイズ線形回帰 ===
with pm.Model() as linear_model:
    # 事前分布
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta = pm.Normal("beta", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=5)

    # 尤度
    mu = alpha + beta * X
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    # MCMC サンプリング (NUTS)
    trace = pm.sample(
        draws=2000,
        tune=1000,
        chains=4,
        cores=4,
        target_accept=0.95,
        return_inferencedata=True,
    )

# === 収束診断 ===
print(az.summary(trace, var_names=["alpha", "beta", "sigma"]))
# Rhat < 1.01 && ESS > 400 を確認

# === トレースプロット ===
az.plot_trace(trace, var_names=["alpha", "beta", "sigma"])
plt.savefig("figures/bayesian_trace.png", dpi=300, bbox_inches="tight")
plt.show()

2. Hierarchical Bayesian model (multilevel)

import pandas as pd

# === 階層ベイズ回帰 (ランダム効果) ===
# 例: 複数学校の生徒の成績予測
# school_id: グループ変数
# x: 予測変数, y: 目的変数

def hierarchical_model(df):
    school_ids = df["school_id"].unique()
    school_idx = df["school_id"].map({s: i for i, s in enumerate(school_ids)}).values
    n_schools = len(school_ids)

    with pm.Model() as hier_model:
        # Hyper-priors (全体平均)
        mu_alpha = pm.Normal("mu_alpha", mu=0, sigma=10)
        sigma_alpha = pm.HalfNormal("sigma_alpha", sigma=5)
        mu_beta = pm.Normal("mu_beta", mu=0, sigma=10)
        sigma_beta = pm.HalfNormal("sigma_beta", sigma=5)

        # 学校レベルのパラメータ (部分プーリング)
        alpha = pm.Normal("alpha", mu=mu_alpha, sigma=sigma_alpha, shape=n_schools)
        beta = pm.Normal("beta", mu=mu_beta, sigma=sigma_beta, shape=n_schools)

        # 残差
        sigma = pm.HalfNormal("sigma", sigma=5)

        # 尤度
        mu = alpha[school_idx] + beta[school_idx] * df["x"].values
        y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=df["y"].values)

        # サンプリング
        trace = pm.sample(2000, tune=1000, chains=4, target_accept=0.95,
                          return_inferencedata=True)

    return trace

# === Forest Plot（学校間パラメータ比較）===
# az.plot_forest(trace, var_names=["alpha", "beta"], combined=True)

3. Posterior predictive check (PPC)

def posterior_predictive_check(model, trace, observed_y):
    """事後予測チェック: モデルの適合度検証"""
    with model:
        ppc = pm.sample_posterior_predictive(trace, random_seed=42)

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))

    # 1. 密度オーバーレイ
    az.plot_ppc(az.from_pymc3(posterior_predictive=ppc, model=model),
                ax=axes[0], kind="kde")
    axes[0].set_title("Posterior Predictive: KDE")

    # 2. 統計量比較 (mean)
    ppc_mean = ppc.posterior_predictive["y_obs"].mean(dim=("chain", "draw")).values
    axes[1].hist(ppc_mean, bins=30, alpha=0.7, label="PPC mean")
    axes[1].axvline(observed_y.mean(), color="red", linestyle="--", label="Observed mean")
    axes[1].set_title("PPC: Mean Comparison")
    axes[1].legend()

    # 3. 残差チェック
    ppc_median = np.median(ppc.posterior_predictive["y_obs"].values, axis=(0, 1))
    residuals = observed_y - ppc_median
    axes[2].scatter(ppc_median, residuals, alpha=0.5, s=10)
    axes[2].axhline(0, color="red", linestyle="--")
    axes[2].set_title("PPC Residuals")

    plt.tight_layout()
    plt.savefig("figures/bayesian_ppc.png", dpi=300, bbox_inches="tight")
    plt.show()

4. Model comparison (WAIC / LOO-CV)

def bayesian_model_comparison(models_dict):
    """
    WAIC / LOO-CV でベイズモデルを比較。
    models_dict: {"model_name": (model, trace)} の辞書
    """
    comparison_data = {}

    for name, (model, trace) in models_dict.items():
        # WAIC
        waic = az.waic(trace)
        # LOO-CV (PSIS-LOO)
        loo = az.loo(trace)

        comparison_data[name] = {
            "waic": round(waic.waic, 2),
            "waic_se": round(waic.waic_se, 2),
            "p_waic": round(waic.p_waic, 2),
            "loo": round(loo.loo, 2),
            "loo_se": round(loo.loo_se, 2),
            "p_loo": round(loo.p_loo, 2),
        }

    # ArviZ モデル比較
    comp_df = az.compare(
        {name: trace for name, (_, trace) in models_dict.items()},
        ic="loo",
    )
    print("Model Comparison (LOO-CV):")
    print(comp_df)

    # 比較プロット
    az.plot_compare(comp_df)
    plt.savefig("figures/bayesian_model_comparison.png", dpi=300, bbox_inches="tight")
    plt.show()

    return comparison_data

5. Bayesian optimization (Gaussian Process)

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
from scipy.optimize import minimize
from scipy.stats import norm

def bayesian_optimization(objective_func, bounds, n_init=5, n_iter=25):
    """
    Gaussian Process ベースのベイズ最適化。
    Expected Improvement (EI) を獲得関数として使用。
    """
    dim = len(bounds)

    # 初期サンプリング (Latin Hypercube)
    X_samples = np.random.uniform(
        [b[0] for b in bounds],
        [b[1] for b in bounds],
        size=(n_init, dim)
    )
    y_samples = np.array([objective_func(x) for x in X_samples])

    # GP モデル
    kernel = Matern(nu=2.5)
    gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-6, normalize_y=True)

    history = {"X": list(X_samples), "y": list(y_samples)}

    for i in range(n_iter):
        gp.fit(np.array(history["X"]), np.array(history["y"]))

        # Expected Improvement
        def neg_ei(x):
            mu, sigma = gp.predict(x.reshape(1, -1), return_std=True)
            best_y = np.min(history["y"])
            z = (best_y - mu) / (sigma + 1e-8)
            ei = (best_y - mu) * norm.cdf(z) + sigma * norm.pdf(z)
            return -ei[0]

        # 次の探索点
        best_x = None
        best_ei = np.inf
        for _ in range(20):
            x0 = np.random.uniform([b[0] for b in bounds], [b[1] for b in bounds])
            result = minimize(neg_ei, x0, bounds=bounds, method="L-BFGS-B")
            if result.fun < best_ei:
                best_ei = result.fun
                best_x = result.x

        new_y = objective_func(best_x)
        history["X"].append(best_x)
        history["y"].append(new_y)

        print(f"Iter {i+1}: x = {best_x}, y = {new_y:.4f}, best = {min(history['y']):.4f}")

    return history

6. Modeling with Stan

# === CmdStanPy インターフェース ===
import cmdstanpy

def fit_stan_model(stan_code, data_dict, **kwargs):
    """Stan モデルのコンパイルとフィット"""
    model = cmdstanpy.CmdStanModel(stan_file=None, stan_file_path=None,
                                     model_code=stan_code)
    fit = model.sample(data=data_dict, chains=4, iter_sampling=2000,
                       iter_warmup=1000, **kwargs)

    # ArviZ 変換
    idata = az.from_cmdstanpy(fit)
    return idata

# Stan モデルコード例: 階層正規モデル
STAN_HIERARCHICAL = """
data {
  int<lower=0> N;
  int<lower=0> J;
  array[N] int<lower=1,upper=J> group;
  vector[N] x;
  vector[N] y;
}
parameters {
  real mu_alpha;
  real mu_beta;
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;
  real<lower=0> sigma;
  vector[J] alpha;
  vector[J] beta;
}
model {
  mu_alpha ~ normal(0, 10);
  mu_beta ~ normal(0, 10);
  sigma_alpha ~ half_normal(5);
  sigma_beta ~ half_normal(5);
  sigma ~ half_normal(5);
  alpha ~ normal(mu_alpha, sigma_alpha);
  beta ~ normal(mu_beta, sigma_beta);
  y ~ normal(alpha[group] + beta[group] .* x, sigma);
}
"""

---

Best Practices

1. Weakly informative priors: Let the data speak. Avoid overly informative priors
2. Convergence diagnostics required: Rhat < 1.01, ESS > 400/chain, no divergences
3. Posterior predictive check: Always verify model fit with PPC
4. LOO-CV preferred: LOO-CV (PSIS-LOO) is more stable than WAIC for model comparison
5. Centered parameterization: Consider non-centered parameterization for hierarchical models
6. Reproducibility: Fix random_seed
7. Sensitivity analysis: Check the stability of results when priors are changed

Completeness Checklist

• Definition of model specification (likelihood and priors)
• Prior predictive check
• MCMC sampling completed
• Convergence diagnostics (Rhat, ESS, divergences)
• Summary of posterior distributions (mean, HDI)
• Posterior predictive check (PPC)
• Model comparison (LOO-CV / WAIC)
• Result report and plot generation

ToolUniverse integration

TU Key	Tool name	Integration content
biotools	bio.tools	Bayesian statistics tool registry search

References

Output Files

File	Format	Generation timing
results/bayesian_summary.json	Posterior distribution summary (JSON)	At sampling completion
figures/bayesian_trace.png	Trace plot	At MCMC completion
figures/bayesian_ppc.png	Posterior predictive check plot	At PPC completion
figures/bayesian_model_comparison.png	Model comparison plot	At comparison completion

Reference skills

Skill	Integration
scientific-statistical-testing	← Foundation for comparison with frequentist tests
scientific-ml-regression	← Bayesian extension of regression models
scientific-doe	← Bayesian optimization of experimental design
scientific-causal-inference	← Bayesian causal model
scientific-quantum-computing	→ Bayesian optimization of VQE parameters
scientific-meta-analysis	← Bayesian meta-analysis