PyMC を用いたベイズ推論：実践と応用
ベイズ推論の概要
ベイズ推論は、確率論に基づいた統計的手法であり、未知のパラメータに対する事前の信念（事前分布）を、観測されたデータ（尤度）を用いて更新し、事後分布を得るプロセスです。
この事後分布は、データが観測された後のパラメータに関する最も確からしい情報を提供します。ベイズ推論の核心は、「定理」であるベイズの定理に基づいています。
P(θ|D) = P(D|θ) * P(θ) / P(D)
ここで、
P(θ|D)は事後分布：データDが与えられたときのパラメータθの確率P(D|θ)は尤度：パラメータθのもとでデータDが得られる確率P(θ)は事前分布：データDを観測する前のパラメータθの確率P(D)は周辺尤度（証拠）：データDが得られる確率（正規化定数）
PyMCは、このベイズ推論をPythonで容易に実装するための強力なライブラリです。統計モデルの定義、サンプリング、そして結果の分析まで、一連のプロセスをサポートします。
PyMC の基本構造と機能
PyMCを用いたベイズ推論は、主に以下のステップで構成されます。
1. モデルの定義
PyMCでは、pm.Model()コンテキストマネージャ内で、確率変数（パラメータや観測データ）とその確率分布を定義します。確率分布は、pymc.distributionsモジュールから選択します。例えば、正規分布であればpm.Normal、一様分布であればpm.Uniformなどです。
パラメータは、事前分布を持つ確率変数として定義されます。観測データは、そのパラメータが生成したと仮定される確率変数として定義され、observed引数に実際のデータを指定します。
2. 推論（サンプリング）
モデル定義後、MCMC（Markov Chain Monte Carlo）アルゴリズムを用いて、事後分布からのサンプルを生成します。PyMCは、NUTS（No-U-Turn Sampler）などの効率的なサンプラーを標準で提供しています。
pm.sample()関数を実行することで、指定したステップ数だけマルコフ連鎖を生成し、事後分布の近似を得ます。
3. 結果の分析
サンプリングされたデータは、事後分布からのサンプルとして扱われます。これらのサンプルを用いて、パラメータの平均値、中央値、信用区間（Credible Interval）などの要約統計量を計算できます。また、事後分布の形状を可視化するために、ヒストグラムやカーネル密度推定（KDE） plots を描画することも一般的です。
PyMCは、arvizライブラリと連携することで、これらの分析や可視化を容易に行うことができます。arviz.summary()やarviz.plot_trace()などの関数が便利です。
PyMC を用いた実践例：単純な線形回帰
ここでは、簡単な線形回帰モデルをPyMCで実装する例を示します。モデルは以下のように定義されます。
y = α + βx + ε
ここで、αは切片、βは傾き、εは誤差項であり、ε ~ Normal(0, σ)と仮定します。
モデル定義
import pymc as pm
import numpy as np
# サンプルデータの生成
np.random.seed(0)
x = np.linspace(0, 10, 100)
α_true = 1.0
β_true = 2.0
σ_true = 1.0
y = α_true + β_true * x + np.random.normal(0, σ_true, 100)
with pm.Model() as linear_model:
# 事前分布の設定
α = pm.Normal('α', mu=0, sigma=10) # 切片の事前分布
β = pm.Normal('β', mu=0, sigma=10) # 傾きの事前分布
σ = pm.HalfCauchy('σ', beta=1) # 誤差の標準偏差の事前分布（非負なのでHalfCauchy）
# 尤度の定義
μ = α + β * x
y_obs = pm.Normal('y_obs', mu=μ, sigma=σ, observed=y) # 観測データの尤度
# サンプリングの実行
trace = pm.sample(2000, tune=1000, cores=2) # 2000個のサンプル、1000個のチューニングステップ
結果の分析
サンプリング後、traceオブジェクトに事後分布からのサンプルが格納されます。arvizライブラリを用いて、結果の要約や可視化を行います。
import arviz as az
# 結果の要約
summary = az.summary(trace, hdi_prob=0.94) # 94%信用区間
print(summary)
# トレースプロット（事後分布の可視化）
az.plot_trace(trace)
# plt.show() # 必要に応じて表示
summaryオブジェクトには、各パラメータの事後平均、標準偏差、信用区間などが表示されます。plot_traceは、各パラメータの事後分布とサンプリング過程を可視化し、MCMCの収束性を確認するのに役立ちます。
PyMC の高度なトピックと応用
PyMCは、単純なモデルだけでなく、より複雑な統計モデリングにも対応しています。
階層モデル（Hierarchical Models）
階層モデルは、グループ間でパラメータが共有されるような構造を持つデータをモデリングする際に強力です。例えば、複数の学校の生徒の成績を分析する際に、各学校固有のパラメータと、全学校に共通するパラメータを同時に推定できます。
欠損値のモデリング
PyMCを用いると、観測データに欠損がある場合でも、その欠損値をモデルの一部として定義し、同時に推定することができます。これにより、単純な欠損値補完よりも柔軟で、不確実性を考慮した分析が可能になります。
時系列分析
ARIMAモデルや状態空間モデルなど、時系列データ特有の構造をPyMCで表現し、ベイズ推論を行うことができます。
モデル診断と選択
MCMCサンプリングが適切に収束しているかを確認するための診断手法（例：R-hat統計量、トレースプロット）が重要です。また、異なるモデルの比較には、WAIC（Watanabe-Akaike Information Criterion）やLOO-CV（Leave-One-Out Cross-Validation）などが用いられます。
カスタム分布
PyMCは、標準で提供されている分布以外にも、ユーザーが独自の確率分布を定義してモデルに組み込む機能も備えています。
まとめ
PyMCは、Pythonにおいてベイズ推論を実践するための洗練されたライブラリです。その直感的なモデル定義構文、効率的なサンプリングアルゴリズム、そしてarvizとの連携による強力な分析・可視化機能は、統計モデリングの専門家からデータサイエンティストまで、幅広いユーザーにとって価値のあるツールとなっています。複雑な統計的問題に取り組む上で、PyMCはベイズ的なアプローチを現実のものとしてくれます。
