ベイズ推定をPythonで実装する医療データ分析の完全ガイド

ベイズ推定をPythonで実装する医療データ分析の基礎と応用

「頻度論で十分」と思っているあなたは、検査の偽陽性率を3倍以上過大評価しているかもしれません。

ベイズ推定の基本概念とPythonで実装する前に理解すべき数式

ベイズ推定は「新しいデータが得られるたびに、自分の確信度を更新していく」という考え方を数式で表したものです。頻度論的な統計が「データの反復試行における長期的な確率」を扱うのに対し、ベイズ推定は「今持っている情報に基づく信念の度合い」を扱います。この違いは医療現場においてきわめて重要です。

中心となる数式はベイズの定理で、次のように表されます。

$$P(H|D) = \frac{P(D|H) \cdot P(H)}{P(D)}$$

各項の意味は以下のとおりです。

• P(H)：事前確率（Prior）。データを見る前に、ある仮説Hがどれくらい確からしいかを表します。医療でいえば「検査前の疾患有病率」に相当します。
• P(D|H)：尤度（Likelihood）。仮説Hが正しいとしたとき、観測データDが得られる確率です。「感度（真陽性率）」がこれにあたります。
• P(H|D)：事後確率（Posterior）。データDを観測した後に、仮説Hがどれくらい確からしいかを表します。「陽性的中率（PPV）」の考え方そのものです。
• P(D)：周辺尤度（Evidence）。正規化定数として機能し、すべての仮説を考慮した際のデータの観測確率です。

つまり「事後確率 ∝ 尤度 × 事前確率」が基本です。

Pythonで実装する前に、この関係をしっかり理解しておくことが、後のコードを読み解く上での土台になります。特に「事前分布（Prior Distribution）」の選択が推定結果に大きく影響するという点は、医療データ分析において見落とされがちな落とし穴です。

例えば、ある感染症の有病率が人口全体で0.1%（事前確率0.001）であるとします。感度95%、特異度99%の検査で陽性が出た場合の陽性的中率をベイズの定理で計算すると、実は約8.7%にしかなりません。「陽性＝ほぼ確定」と直感的に考えがちですが、計算すると約91%が偽陽性という結果です。意外ですね。

この事実は医療現場でも十分に意識されているとは言えず、ベイズ推定をPythonで実装して可視化することで、このような直感との乖離を明確に示せます。

ベイズ推定のPython実装に使うライブラリ：PyMC・SciPy・NumPyの選び方

Pythonでベイズ推定を実装する際に使用する主なライブラリは3系統に分かれています。これらは目的と習熟度に応じて使い分けることが実践上の効率を左右します。

まず最も直感的な実装に向いているのがNumPy + SciPyによる解析的アプローチです。共役事前分布が使える場合（例：二項分布とベータ分布の組み合わせ）は、後述の数式を直接コードに落とし込めます。

```python
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# 例：コインの表が出る確率をベイズ推定する
# 事前分布：Beta(1, 1)（一様分布 = 無情報事前分布）
alpha_prior = 1
beta_prior = 1

# 観測データ：20回中14回表
n_heads = 14
n_tails = 6

# 事後分布：Beta(alpha_prior + n_heads, beta_prior + n_tails)
alpha_post = alpha_prior + n_heads
beta_post = beta_prior + n_tails

# 事後分布の可視化
x = np.linspace(0, 1, 200)
prior = stats.beta.pdf(x, alpha_prior, beta_prior)
posterior = stats.beta.pdf(x, alpha_post, beta_post)

plt.figure(figsize=(8, 4))
plt.plot(x, prior, label='事前分布 Beta(1,1)', linestyle='--')
plt.plot(x, posterior, label=f'事後分布 Beta({alpha_post},{beta_post})')
plt.xlabel('表が出る確率 θ')
plt.ylabel('確率密度')
plt.legend()
plt.title('ベイズ更新の可視化')
plt.tight_layout()
plt.show()

print(f"事後分布の平均（MAP推定）: {alpha_post / (alpha_post + beta_post):.3f}")
```

これは使えそうです。

次に、より複雑なモデルを扱うならPyMC（旧PyMC3）が強力な選択肢です。PyMCはMCMC（マルコフ連鎖モンテカルロ法）やVI（変分推論）を自動化してくれるため、共役事前分布が存在しない問題にも対応できます。2025年時点ではPyMC v5が主流で、NumPyroやTensorFlow Probabilityと比較しても医療・製薬分野での採用実績が豊富です。

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

# 観測データ（例：ある処置後の回復人数）
n_patients = 50
n_recovered = 31

with pm.Model() as model:
# 事前分布：回復率は0〜1の一様分布（無情報事前分布）
theta = pm.Beta('theta', alpha=1, beta=1)

# 尤度：二項分布
obs = pm.Binomial('obs', n=n_patients, p=theta, observed=n_recovered)

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

# 結果の確認
az.plot_posterior(trace, var_names='theta')
print(az.summary(trace, var_names='theta'))
```

そしてNumPyroはJAXバックエンドを使うことで高速化が図れ、大規模な医療データセット（数万件以上の電子カルテなど）を扱う際に優位性があります。ただし習熟コストが高いため、最初はPyMCから入るのが原則です。

ライブラリ選択の基準をまとめると以下のようになります。

ライブラリ	適した場面	難易度	速度
NumPy + SciPy	共役モデル・教育目的	★☆☆	速い
PyMC	中〜大規模・実務分析	★★☆	中程度
NumPyro	大規模・高速化が必要	★★★	非常に速い

ベイズ推定のPython実装で医療現場の陽性的中率を計算する具体例

医療従事者にとって最も身近なベイズ推定の応用は、検査の陽性的中率（PPV：Positive Predictive Value）と陰性的中率（NPV）の計算です。頻度論的な感度・特異度の数字だけでは、実際の患者における意味が変わります。これが原則です。

以下のコードは、有病率・感度・特異度を入力すれば自動的にPPVとNPVを計算し、有病率の変化に応じた推移をグラフ化するものです。

```python
import numpy as np
import matplotlib.pyplot as plt

def calc_ppv_npv(prevalence, sensitivity, specificity):
"""
有病率・感度・特異度からPPVとNPVを計算
"""
tp = prevalence * sensitivity
fp = (1 - prevalence) * (1 - specificity)
fn = prevalence * (1 - sensitivity)
tn = (1 - prevalence) * specificity

ppv = tp / (tp + fp) if (tp + fp) > 0 else 0
npv = tn / (tn + fn) if (tn + fn) > 0 else 0
return ppv, npv

# 感度95%、特異度99%の検査を想定
sensitivity = 0.95
specificity = 0.99

# 有病率を0.1%〜50%の範囲で変化させてPPVを計算
prevalence_range = np.linspace(0.001, 0.5, 200)
ppv_list =
npv_list =

for prev in prevalence_range:
ppv, npv = calc_ppv_npv(prev, sensitivity, specificity)
ppv_list.append(ppv)
npv_list.append(npv)

# グラフ化
plt.figure(figsize=(9, 5))
plt.plot(prevalence_range * 100, ppv_list, label='PPV（陽性的中率）', color='crimson')
plt.plot(prevalence_range * 100, npv_list, label='NPV（陰性的中率）', color='steelblue')
plt.axvline(x=0.1, color='gray', linestyle=':', label='有病率0.1%（低リスク集団）')
plt.axvline(x=10, color='orange', linestyle=':', label='有病率10%（中リスク集団）')
plt.xlabel('有病率（%）')
plt.ylabel('的中率')
plt.title('有病率の変化によるPPV・NPVの推移（感度95%, 特異度99%）')
plt.legend()
plt.tight_layout()
plt.show()

# 具体的な数値も出力
for prev_pct in 0.1, 1.0, 5.0, 10.0, 20.0:
ppv, npv = calc_ppv_npv(prev_pct / 100, sensitivity, specificity)
print(f"有病率{prev_pct:5.1f}% → PPV: {ppv*100:.1f}% NPV: {npv*100:.2f}%")
```

このコードを実行すると、有病率0.1%のスクリーニング検査では感度95%・特異度99%でもPPVが8.7%しかない一方、有病率10%の精密検査場面ではPPVが91.4%まで上昇することが一目で確認できます。

数字で言うと、スクリーニング陽性100人のうち実際に疾患があるのは9人程度ということです。この情報を患者へのインフォームドコンセントに活かすためにも、Pythonで動的に計算できる環境を整えておく価値があります。

有病率が検査の解釈を根本から変えるということですね。このような計算ツールをJupyter Notebookで院内共有しておくだけで、部署全体の検査解釈のリテラシーを底上げできます。

参考として、検査性能と有病率に関する基礎的な統計的背景はこちらをご覧ください。

ベイズ推定でPython実装するMCMCの収束診断と失敗しないチューニング方法

MCMCベースのベイズ推定でよく起きる問題が「サンプルが収束していない」状態です。収束していないトレースからは信頼できる推論が得られません。厳しいところですね。

収束を診断するための主な指標として、以下を確認するのが標準的な手順です。

• R-hat（Gelman-Rubin統計量）：複数のチェーンが同じ分布に収束しているかを測る指標。1.01以下が目安で、1.1を超えていたら収束失敗と判断します。PyMCではaz.summary()で自動表示されます。
• ESS（有効サンプルサイズ）：MCMCの自己相関を考慮した実質的なサンプル数。最低でも400以上、可能なら1000以上を目安にします。
• トレースプロット：チェーンが特定の値に貼り付いていたり、周期的なパターンを示す場合は収束失敗のサインです。

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

# 収束診断のデモ（2パラメータの正規モデル）
np.random.seed(42)
data = np.random.normal(loc=5.0, scale=1.5, size=80)

with pm.Model() as normal_model:
mu = pm.Normal('mu', mu=0, sigma=10)
sigma = pm.HalfNormal('sigma', sigma=5)
obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=data)

# chains=4 で複数チェーンを走らせてR-hatを計算可能にする
trace = pm.sample(
draws=2000,
tune=1000,
chains=4,
target_accept=0.9, # NUTSのステップサイズ調整。デフォルト0.8より高めに設定
return_inferencedata=True,
random_seed=42
)

# 収束診断
summary = az.summary(trace, var_names='mu', 'sigma')
print(summary['mean', 'sd', 'hdi_3%', 'hdi_97%', 'r_hat', 'ess_bulk'])

# トレースプロット
az.plot_trace(trace, var_names='mu', 'sigma')
```

収束しない原因の多くは事前分布が非常に広い（無情報すぎる）か、パラメータ空間のスケールが変数間で大きく異なる場合です。医療データでは測定値の単位が変数によって血圧（mmHg）・体重（kg）・酵素活性（U/L）などと異なるため、モデル構築前に標準化（z-score化）を行うことが収束改善の近道です。

```python
# データの標準化（収束改善のため）
from scipy.stats import zscore

data_scaled = zscore(data) # 平均0、標準偏差1に変換
# ← このデータをobservedに渡すとサンプリングが安定しやすい
```

また、`target_accept`パラメータを0.9〜0.95程度に設定するとNUTSアルゴリズムのステップサイズが小さくなり、急峻な事後分布面での発散エラー（divergence）を減らせます。発散の数は0が理想で、全サンプルの1%を超えるとモデルの再検討が必要です。

事前分布の選択とスケーリングが収束のカギということです。

ベイズ推定Python実装で医療従事者が見落としがちな「階層モデル」の活用法

一般的なベイズ推定の記事ではあまり取り上げられませんが、医療データ分析で特に効果を発揮するのが階層ベイズモデル（Hierarchical Bayesian Model）です。

病院Aと病院Bで同じ処置をしたとき、それぞれの成功率を独立に推定するのは実は非効率です。なぜなら、患者数が少ない病院では推定の不確実性が大きくなり、結果が不安定になりやすいからです。階層モデルは「各病院の成功率は、全病院共通の母集団分布から引き出されている」と仮定することで、データが少ない施設の推定を「全体の平均に引き寄せる（shrinkage）」効果をもたらします。

これはまるで、研修医の評価を「その病院だけで完結させる」のではなく「全国平均も参照しながら調整する」ようなものです。

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

# 5つの施設での処置成功データ（件数が少ない施設も含む）
n_hospitals = 5
n_treated = np.array(120, 45, 200, 18, 90) # 各施設の処置件数
n_success = np.array( 89, 30, 161, 10, 72) # 各施設の成功件数

with pm.Model() as hierarchical_model:
# 超事前分布（全施設共通の母集団パラメータ）
mu_logit = pm.Normal('mu_logit', mu=0, sigma=1)
sigma_logit = pm.HalfNormal('sigma_logit', sigma=1)

# 各施設のlogitスケールの成功率
theta_logit = pm.Normal('theta_logit', mu=mu_logit, sigma=sigma_logit,
shape=n_hospitals)

# 確率スケールに変換
theta = pm.Deterministic('theta', pm.math.sigmoid(theta_logit))

# 尤度
obs = pm.Binomial('obs', n=n_treated, p=theta, observed=n_success)

trace_hier = pm.sample(2000, tune=1000, chains=4,
return_inferencedata=True, random_seed=42)

# 各施設の成功率と不確実性の可視化
az.plot_forest(trace_hier, var_names='theta',
combined=True, hdi_prob=0.94)
```

このモデルで特に注目すべきは、処置件数が18件と少ない施設4の推定値が、全体の傾向に引き寄せられて安定する点です。独立推定なら成功率55.6%（10/18）という不安定な値になりますが、階層モデルでは全施設の情報を借用して適切に補正されます。

これが階層ベイズの本質です。多施設共同研究や地域間比較分析、あるいはコホート研究での施設効果の調整など、実臨床研究で直面する多くの場面でこのアプローチは威力を発揮します。

階層モデルの詳細な理論背景と実装については、以下の資料が参考になります。

PyMC公式ドキュメント：Hierarchical Partial Pooling（階層的部分プーリングの実装例）

また、ベイズ統計の医療応用に関する日本語の学術的解説としては以下も有用です。

ベイズ推定をPythonで実装する際の事前分布の選択と医療データへの応用まとめ

ここまでの内容を踏まえると、ベイズ推定のPython実装で最も意思決定に影響するのが「事前分布（Prior Distribution）」の選択です。事前分布は単なる数学的な設定ではなく、「この分析に組み込む専門知識の量」を決める設計判断です。

事前分布の主な種類と医療データへの対応を整理します。

事前分布の種類	特徴	医療データでの使用例
無情報事前分布（一様分布など）	特定の値を優遇しない。データに完全に依存	新薬の有効率（既存知識なし）
弱情報事前分布（Half-Normal, Normal）	パラメータの範囲に緩い制約。実務で最も推奨	副作用発現率、検査値の基準範囲推定
情報的事前分布（過去のメタ解析結果など）	既存のエビデンスを反映。サンプルが少ない場面で強力	希少疾患の臨床試験、小規模コホート

弱情報事前分布が実務では最も安全な選択です。

特に医療従事者が注意すべきなのは、「事前分布を意識せずに無情報事前分布を使うと、サンプル数が少ない場面で非常識な推定値（例：有効率150%など）が事後分布の裾に入ることがある」という点です。これを防ぐためのPyMCコードの書き方を示します。

```python
import pymc as pm

# ❌ 避けるべき例：完全な無情報事前分布
with pm.Model() as bad_model:
# Uniform(0, 1)は問題ないが、連続値の場合は注意
rate = pm.Uniform('rate', lower=0, upper=1)

# ✅ 推奨例：弱情報事前分布を使う
with pm.Model() as good_model:
# 有効率が0〜1の範囲に収まり、かつ極端な値に対してペナルティを与える
rate = pm.Beta('rate', alpha=2, beta=2)
# Beta(2,2)は「0.5付近が多少優遇されるが、0や1も除外しない」弱情報事前分布

# 別の場面：治療効果量（連続値）に弱情報事前分布
effect_size = pm.Normal('effect_size', mu=0, sigma=1)
# 「効果ゼロ付近が事前に最も妥当」というドメイン知識を反映
```

また、事前分布の感度分析（Prior Sensitivity Analysis）は臨床研究での報告に際して推奨されます。使用した事前分布を変えても結論が変わらないことを示すことで、研究の頑健性を担保できます。PyMCでは複数のモデルを定義してaz.compare()でWAIC（広義赤池情報量規準）を比較する方法が標準的です。

```python
import arviz as az

# 2つのモデルのトレースを比較（model_trace_1, model_trace_2 が事前に定義済みの場合）
# comparison = az.compare({'モデルA': trace_a, 'モデルB': trace_b})
# print(comparison)
# WAICが小さい方が予測精度の高いモデルと判断する
```

ベイズ推定のPython実装は、コードを書くことが目的ではありません。最終的な目的は「限られた医療データから、より適切な臨床判断を引き出すこと」です。PyMCやSciPyはその手段であり、事前分布という形でドメイン知識をモデルに組み込める点こそが、頻度論的手法にはないベイズ推定の最大の強みです。

実装を始める際は、まずSciPyによる解析的ベイズ更新から試し、次にPyMCで階層モデルへと段階的に進むルートが挫折しにくい方法です。ベイズ推定の学習を体系的に進めるには以下の書籍・オンラインリソースも活用できます。

岩波書店：「ベイズ統計の理論と方法」久保拓弥 著（日本語ベイズ統計の定番参考書）