PyMCで学ぶ確率的プログラミング — ベイズモデリングの実践入門

import pymc as pm
import numpy as np

# データ
observed_data = np.array([1, 0, 1, 1, 0, 1, 0, 1, 1, 1])

with pm.Model() as model:
    # 事前分布
    theta = pm.Beta("theta", alpha=1, beta=1)

    # 尤度（データの生成モデル）
    y = pm.Bernoulli("y", p=theta, observed=observed_data)

    # MCMCサンプリング
    trace = pm.sample(2000, tune=1000, cores=2, random_seed=42)

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

# データ
n_flips = 20
n_heads = 14
data = np.concatenate([np.ones(n_heads), np.zeros(n_flips - n_heads)])

# PyMCモデル
with pm.Model() as coin_model:
    # 事前分布: Beta(1, 1) = 一様分布
    theta = pm.Beta("theta", alpha=1, beta=1)

    # 尤度: ベルヌーイ分布
    y = pm.Bernoulli("y", p=theta, observed=data)

    # MCMCサンプリング
    trace = pm.sample(3000, tune=1000, cores=2, random_seed=42,
                      return_inferencedata=True)

# 結果の可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: 事後分布 vs 解析解
ax = axes[0]
theta_samples = trace.posterior["theta"].values.flatten()
ax.hist(theta_samples, bins=50, density=True, alpha=0.5, color='blue',
        label='MCMCサンプル')

x = np.linspace(0, 1, 200)
analytic_posterior = stats.beta(15, 7)
ax.plot(x, analytic_posterior.pdf(x), 'r-', linewidth=2,
        label='解析解 Beta(15, 7)')
ax.set_xlabel(r'$\theta$')
ax.set_ylabel('密度')
ax.set_title('事後分布: MCMC vs 解析解')
ax.legend(fontsize=11)

# 右: トレースプロット
ax = axes[1]
for chain in range(trace.posterior.dims['chain']):
    chain_samples = trace.posterior["theta"].values[chain, :]
    ax.plot(chain_samples, alpha=0.5, linewidth=0.5, label=f'Chain {chain}')
ax.set_xlabel('サンプル番号')
ax.set_ylabel(r'$\theta$')
ax.set_title('トレースプロット')
ax.legend(fontsize=11)

plt.tight_layout()
plt.savefig("pymc_coin_flip.png", dpi=150, bbox_inches="tight")
plt.show()

# 要約統計量
print(az.summary(trace, var_names=["theta"]))

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

# データ生成
np.random.seed(42)
N = 50
true_intercept = 1.5
true_slope = 2.3
true_sigma = 0.8

x = np.random.uniform(0, 5, N)
y = true_intercept + true_slope * x + np.random.normal(0, true_sigma, N)

# PyMCモデル
with pm.Model() as linear_model:
    # 事前分布
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=5)

    # 線形モデル
    mu = intercept + slope * x

    # 尤度
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    # サンプリング
    trace = pm.sample(3000, tune=1000, cores=2, random_seed=42,
                      return_inferencedata=True)

# 結果の表示
print(az.summary(trace, var_names=["intercept", "slope", "sigma"]))

# 事後予測分布の可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: 回帰直線の不確実性
ax = axes[0]
ax.scatter(x, y, c='gray', s=20, alpha=0.6, label='データ')

# 事後サンプルから回帰直線を描画
x_plot = np.linspace(0, 5, 100)
intercept_samples = trace.posterior["intercept"].values.flatten()
slope_samples = trace.posterior["slope"].values.flatten()

# ランダムに100本の回帰直線を描画
idx = np.random.choice(len(intercept_samples), 100, replace=False)
for i in idx:
    ax.plot(x_plot, intercept_samples[i] + slope_samples[i] * x_plot,
            'b-', alpha=0.03, linewidth=1)

# 事後平均の回帰直線
ax.plot(x_plot, np.mean(intercept_samples) + np.mean(slope_samples) * x_plot,
        'r-', linewidth=2, label='事後平均')

# 真の直線
ax.plot(x_plot, true_intercept + true_slope * x_plot,
        'k--', linewidth=1.5, label='真の関数')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('ベイズ線形回帰: 事後回帰直線')
ax.legend()

# 右: パラメータの事後分布
ax = axes[1]
from matplotlib.gridspec import GridSpec

fig2, axes2 = plt.subplots(1, 3, figsize=(15, 4))
for ax2, var, true_val in zip(axes2,
                               ["intercept", "slope", "sigma"],
                               [true_intercept, true_slope, true_sigma]):
    samples = trace.posterior[var].values.flatten()
    ax2.hist(samples, bins=50, density=True, alpha=0.6, color='blue')
    ax2.axvline(true_val, color='red', linestyle='--', linewidth=2,
                label=f'真の値 ({true_val})')
    ax2.set_xlabel(var)
    ax2.set_ylabel('密度')
    ax2.set_title(f'{var}の事後分布')
    ax2.legend()

ax.set_visible(False)
plt.tight_layout()
plt.savefig("pymc_linear_regression.png", dpi=150, bbox_inches="tight")
plt.show()

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

# 8学校のデータ
y_obs = np.array([28.39, 7.94, -2.75, 6.82, -0.64, 0.63, 18.01, 12.16])
sigma = np.array([14.9, 10.2, 16.3, 11.0, 9.4, 11.4, 10.4, 17.6])
school_names = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']
J = len(y_obs)

# 非心パラメータ化（ファネル問題の回避）
with pm.Model() as eight_schools:
    # ハイパーパラメータ
    mu = pm.Normal("mu", mu=0, sigma=10)
    tau = pm.HalfNormal("tau", sigma=10)

    # 非心パラメータ化
    theta_raw = pm.Normal("theta_raw", mu=0, sigma=1, shape=J)
    theta = pm.Deterministic("theta", mu + tau * theta_raw)

    # 尤度
    y = pm.Normal("y", mu=theta, sigma=sigma, observed=y_obs)

    # サンプリング
    trace = pm.sample(4000, tune=2000, cores=2, random_seed=42,
                      target_accept=0.95, return_inferencedata=True)

print(az.summary(trace, var_names=["mu", "tau", "theta"]))

# 結果の可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 左上: 各学校の事後分布
ax = axes[0, 0]
theta_samples = trace.posterior["theta"].values.reshape(-1, J)
positions = np.arange(J)

bp = ax.boxplot([theta_samples[:, j] for j in range(J)],
                positions=positions, widths=0.6, patch_artist=True,
                showfliers=False)
for patch in bp['boxes']:
    patch.set_facecolor('lightblue')
    patch.set_alpha(0.7)

ax.scatter(positions, y_obs, c='red', s=80, zorder=5, marker='x',
           linewidth=2, label='観測値')
ax.axhline(np.mean(trace.posterior["mu"].values), color='green',
           linestyle='--', linewidth=1.5, label=r'全体平均 $\mu$')
ax.set_xticks(positions)
ax.set_xticklabels(school_names)
ax.set_xlabel('学校')
ax.set_ylabel('効果量')
ax.set_title('各学校の効果量の事後分布')
ax.legend()

# 右上: 縮小推定の可視化
ax = axes[0, 1]
theta_means = np.mean(theta_samples, axis=0)
mu_mean = np.mean(trace.posterior["mu"].values)

for j in range(J):
    ax.plot([y_obs[j], theta_means[j]], [j, j], 'b-', linewidth=2, alpha=0.7)
    ax.scatter(y_obs[j], j, c='red', s=60, zorder=5, marker='x', linewidth=2)
    ax.scatter(theta_means[j], j, c='blue', s=60, zorder=5)

ax.axvline(mu_mean, color='green', linestyle='--', linewidth=1.5,
           label=r'全体平均 $\mu$')
ax.set_yticks(range(J))
ax.set_yticklabels(school_names)
ax.set_xlabel('効果量')
ax.set_ylabel('学校')
ax.set_title('縮小推定: 観測値(赤x) → 事後平均(青●)')
ax.legend()

# 左下: muとtauの事後分布
ax = axes[1, 0]
mu_samples = trace.posterior["mu"].values.flatten()
ax.hist(mu_samples, bins=50, density=True, alpha=0.6, color='blue',
        label=r'$\mu$ の事後分布')
ax.set_xlabel(r'$\mu$ (全体平均効果)')
ax.set_ylabel('密度')
ax.set_title(r'$\mu$ の事後分布')
ax.legend()

ax = axes[1, 1]
tau_samples = trace.posterior["tau"].values.flatten()
ax.hist(tau_samples, bins=50, density=True, alpha=0.6, color='orange',
        label=r'$\tau$ の事後分布')
ax.set_xlabel(r'$\tau$ (学校間標準偏差)')
ax.set_ylabel('密度')
ax.set_title(r'$\tau$ の事後分布')
ax.legend()

plt.tight_layout()
plt.savefig("pymc_eight_schools.png", dpi=150, bbox_inches="tight")
plt.show()

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

# 8学校モデルの trace を使用（再実行が必要な場合は例題3のコードを実行）
# ここではダミーの簡単なモデルで実演

np.random.seed(42)
data = np.random.normal(5, 2, 100)

with pm.Model() as demo_model:
    mu = pm.Normal("mu", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=5)
    y = pm.Normal("y", mu=mu, sigma=sigma, observed=data)
    trace = pm.sample(3000, tune=1000, cores=4, chains=4,
                      random_seed=42, return_inferencedata=True)

# R-hat と ESS の確認
summary = az.summary(trace, var_names=["mu", "sigma"])
print("=== 収束診断 ===\n")
print(summary[["mean", "sd", "hdi_3%", "hdi_97%", "ess_bulk", "ess_tail", "r_hat"]])

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# トレースプロット
for i, var in enumerate(["mu", "sigma"]):
    # トレース
    ax = axes[i, 0]
    for chain in range(trace.posterior.dims['chain']):
        samples = trace.posterior[var].values[chain, :]
        ax.plot(samples, alpha=0.5, linewidth=0.5)
    ax.set_xlabel('サンプル番号')
    ax.set_ylabel(var)
    ax.set_title(f'{var}: トレースプロット')

    # 自己相関
    ax = axes[i, 1]
    for chain in range(trace.posterior.dims['chain']):
        samples = trace.posterior[var].values[chain, :]
        n = len(samples)
        max_lag = 50
        autocorr = np.array([np.corrcoef(samples[:n-lag], samples[lag:])[0, 1]
                            for lag in range(1, max_lag + 1)])
        ax.plot(range(1, max_lag + 1), autocorr, alpha=0.5)
    ax.axhline(0, color='black', linestyle='-', linewidth=0.5)
    ax.set_xlabel('ラグ')
    ax.set_ylabel('自己相関')
    ax.set_title(f'{var}: 自己相関プロット')

plt.tight_layout()
plt.savefig("pymc_diagnostics.png", dpi=150, bbox_inches="tight")
plt.show()

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

# 簡単なモデルで ArviZ の可視化を実演
np.random.seed(42)
x = np.random.uniform(0, 5, 50)
y = 1.5 + 2.3 * x + np.random.normal(0, 0.8, 50)

with pm.Model() as demo:
    a = pm.Normal("intercept", mu=0, sigma=10)
    b = pm.Normal("slope", mu=0, sigma=10)
    s = pm.HalfNormal("sigma", sigma=5)
    mu = a + b * x
    y_obs = pm.Normal("y_obs", mu=mu, sigma=s, observed=y)
    trace = pm.sample(3000, tune=1000, cores=2, random_seed=42,
                      return_inferencedata=True)

# フォレストプロット（HDI区間の比較）
fig, ax = plt.subplots(figsize=(10, 4))
az.plot_forest(trace, var_names=["intercept", "slope", "sigma"],
               combined=True, hdi_prob=0.94, ax=ax)
ax.set_title('フォレストプロット (94% HDI)')
plt.tight_layout()
plt.savefig("arviz_forest.png", dpi=150, bbox_inches="tight")
plt.show()

# 事後予測チェック（Posterior Predictive Check）
with demo:
    ppc = pm.sample_posterior_predictive(trace, random_seed=42)

fig, ax = plt.subplots(figsize=(10, 5))
ppc_samples = ppc.posterior_predictive["y_obs"].values.reshape(-1, len(y))
for i in range(min(100, len(ppc_samples))):
    ax.hist(ppc_samples[i], bins=20, alpha=0.02, color='blue', density=True)
ax.hist(y, bins=20, alpha=0.5, color='red', density=True, label='観測データ')
ax.set_xlabel('y')
ax.set_ylabel('密度')
ax.set_title('事後予測チェック (PPC)')
ax.legend()
plt.tight_layout()
plt.savefig("arviz_ppc.png", dpi=150, bbox_inches="tight")
plt.show()