プロ野球のシーズン開幕から1ヶ月が経ち、各選手の打率が出揃ったとします。ある選手が20打数10安打で打率.500を記録しています。しかし、シーズン終了時の打率が.500のまま維持されるとは誰も思わないでしょう。なぜなら、リーグ全体の打率が.260程度であることを「事前情報」として知っているからです。
では、その選手のシーズン打率をどう予測すべきでしょうか。打率.500をそのまま使うのは楽観的すぎますが、リーグ平均の.260を使うのは個人のデータを無視しすぎます。最適な予測は両者の間のどこかにあるはずです。
この「データから事前分布を推定し、個別の推定を改善する」というアイデアが経験ベイズ法（empirical Bayes method）です。通常のベイズ統計では事前分布は分析者が主観的に設定しますが、経験ベイズではデータそのものから事前分布（のハイパーパラメータ）を推定します。
経験ベイズ法を理解すると、以下のような応用が開けます。
- ゲノム科学: 数千の遺伝子の発現差を同時に検定するとき、偽発見率（FDR）の制御に使われる
- スポーツ分析: 選手の能力推定（打率、勝率、レーティング）の改善
- 保険数理: 多数の契約者の事故率推定
- 自然言語処理: 言語モデルにおける確率のスムージング
本記事の内容
- 経験ベイズ法の動機と基本的な考え方
- ジェームズ・スタイン推定量 — 縮小推定の原点
- パラメトリック経験ベイズ法の理論
- ノンパラメトリック経験ベイズ法（Robbinsの方法）
- 多重検定における経験ベイズ法の応用
- Pythonによる実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
経験ベイズ法の基本的な考え方
完全ベイズ vs 経験ベイズ
ベイズ統計の枠組みでは、データ $y$ が与えられたときのパラメータ $\theta$ の推定は次の手順で行われます。
- 事前分布 $p(\theta | \eta)$ を指定する（$\eta$ はハイパーパラメータ）
- 尤度 $p(y | \theta)$ とベイズの定理から事後分布 $p(\theta | y, \eta) \propto p(y | \theta) p(\theta | \eta)$ を計算する
- 事後分布に基づいて推定・予測を行う
ここで問題になるのが、ステップ1のハイパーパラメータ $\eta$ をどう設定するかです。
完全ベイズ（full Bayes）: $\eta$ にも事前分布 $p(\eta)$ を置き、$\theta$ と $\eta$ の同時事後分布 $p(\theta, \eta | y)$ を計算する（階層ベイズ）
経験ベイズ（empirical Bayes）: データ $y$ から $\eta$ の点推定値 $\hat{\eta}$ を求め、それを使って事後分布 $p(\theta | y, \hat{\eta})$ を計算する
経験ベイズは階層ベイズの「近似」と見ることができます。$\eta$ の事後分布が鋭いピークを持つとき（データが多いとき）、点推定 $\hat{\eta}$ を使っても $\eta$ の不確実性を無視する影響は小さくなります。
経験ベイズの特長
経験ベイズの魅力は「ベイズの良いとこ取り」にあります。
頻度論的な保証: 事前分布をデータから推定するため、事前分布の主観的な設定に関する批判を回避できます。さらに、多くの場合、経験ベイズ推定量は頻度論的にも優れた性質（一致性、漸近効率性）を持ちます。
計算の簡便性: 完全ベイズ（階層ベイズ）のように $\eta$ 上の積分やMCMCを実行する必要がなく、$\hat{\eta}$ の計算は通常、最尤推定やモーメント法で効率的に行えます。
多重パラメータ問題への適合性: 多数のパラメータ $\theta_1, \ldots, \theta_p$ を同時に推定する問題で特に威力を発揮します。個々のパラメータのデータは少なくても、全体のデータから事前分布を推定することで、「情報の借り入れ」（borrowing strength）が実現されます。
ただし、$\eta$ の不確実性を無視するため、推定の不確実性（信頼区間）をやや過小評価する傾向があります。
経験ベイズの基本的な考え方を理解したところで、次にその歴史的な原点であるジェームズ・スタイン推定量から話を始めましょう。
ジェームズ・スタイン推定量 — 縮小推定の発見
スタインのパラドクス
1961年、チャールズ・スタインとウィラード・ジェームズは統計学に衝撃を与える結果を発表しました。
問題設定: $p \geq 3$ 個の独立な正規分布からの観測値 $y_i \sim N(\theta_i, 1)$（$i = 1, \ldots, p$）が与えられたとき、$\theta = (\theta_1, \ldots, \theta_p)$ を推定したい。
直感的には、各 $\theta_i$ の最良の推定値は $y_i$ そのもの（最尤推定値）だと思えます。事実、各 $\theta_i$ を個別に考えれば、$y_i$ は $\theta_i$ の最良不偏推定量です。
ところが、スタインは次の驚くべき結果を示しました。
スタインのパラドクス: $p \geq 3$ のとき、最尤推定量 $\hat{\theta}_{\text{MLE}} = y$ は許容的でない（inadmissible）。すなわち、あらゆる $\theta$ に対してより小さい（あるいは等しい）二乗誤差を達成する推定量が存在する。
ジェームズ・スタイン推定量の定義
最尤推定量を支配（dominate）する推定量として提案されたのが、ジェームズ・スタイン推定量です。
$$ \begin{equation} \hat{\theta}_{\text{JS}} = \left(1 – \frac{p – 2}{\|y\|^2}\right) y \end{equation} $$
ここで $\|y\|^2 = \sum_{i=1}^{p} y_i^2$ は観測ベクトルのノルムの二乗です。
この推定量は、全ての観測値を原点に向かって縮小します。縮小の度合いは $1 – (p-2)/\|y\|^2$ で、$\|y\|^2$ が大きいほど（データが原点から離れているほど）縮小が小さくなります。
リスクの比較
ジェームズ・スタイン推定量の二乗誤差リスクは次のように計算できます。
$$ \begin{equation} R(\theta, \hat{\theta}_{\text{JS}}) = E\left[\|\hat{\theta}_{\text{JS}} – \theta\|^2\right] = p – (p – 2)^2 E\left[\frac{1}{\|y\|^2}\right] \end{equation} $$
ここで $\|y\|^2 \sim \chi^2_p(\|\theta\|^2)$（非心χ²分布）なので、$E[1/\|y\|^2] > 0$ が保証されます。
最尤推定量のリスク $R(\theta, \hat{\theta}_{\text{MLE}}) = p$ と比較すると
$$ R(\theta, \hat{\theta}_{\text{JS}}) < p = R(\theta, \hat{\theta}_{\text{MLE}}) \quad \text{for all } \theta $$
つまり、あらゆる $\theta$ に対してジェームズ・スタイン推定量は最尤推定量よりリスクが低い（$p \geq 3$ のとき）のです。
経験ベイズの解釈
ジェームズ・スタイン推定量は、実は経験ベイズ法として自然に導出できます。
$\theta_i \stackrel{\text{iid}}{\sim} N(0, A)$（$A > 0$ は未知）、$y_i | \theta_i \sim N(\theta_i, 1)$ という階層モデルを考えます。$A$ が既知のとき、$\theta_i$ の事後平均は
$$ E[\theta_i | y_i, A] = \frac{A}{A + 1} y_i = \left(1 – \frac{1}{A + 1}\right) y_i $$
$A$ が未知なので、データから推定します。$y_i$ の周辺分布は $y_i \sim N(0, A + 1)$ であり、$\sum_{i} y_i^2 \sim (A+1) \chi^2_p$ なので、$A + 1$ の不偏推定量は $\|y\|^2 / p$ です。よって $1/(A+1)$ の推定量として $(p-2)/\|y\|^2$ を代入すると（$p-2$ は不偏性のための補正）、まさにジェームズ・スタイン推定量が得られます。
$$ \hat{\theta}_i^{\text{JS}} = \left(1 – \frac{p-2}{\|y\|^2}\right) y_i $$
この導出は、経験ベイズの核心を端的に示しています。個々のデータからだけでなく、全体のデータから事前分布のパラメータを推定することで、全ての推定が改善されるのです。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# --- ジェームズ・スタイン推定量のリスク比較 ---
p_values = [3, 5, 10, 20]
n_sim = 10000
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) リスクの理論曲線
ax = axes[0]
theta_norm_sq = np.linspace(0, 50, 200)
for p in p_values:
# JS推定量のリスク（近似）
# R(theta, JS) = p - (p-2)^2 * E[1/||y||^2]
# ||y||^2 ~ chi^2_p(||theta||^2) = noncentral chi-squared
risk_js = np.zeros_like(theta_norm_sq)
for j, ncp in enumerate(theta_norm_sq):
# E[1/X] for X ~ chi2(p, ncp) をモンテカルロ推定
samples = np.random.noncentral_chisquare(p, ncp, 5000)
risk_js[j] = p - (p - 2)**2 * np.mean(1.0 / samples)
risk_mle = np.full_like(theta_norm_sq, p)
ax.plot(theta_norm_sq, risk_js / risk_mle, linewidth=2, label=f"p = {p}")
ax.axhline(1.0, color="black", linestyle="--", linewidth=1, alpha=0.5)
ax.set_xlabel(r"$\|\theta\|^2$", fontsize=12)
ax.set_ylabel(r"$R(\theta, \hat{\theta}_{JS}) / R(\theta, \hat{\theta}_{MLE})$", fontsize=12)
ax.set_title("James-Stein Risk Ratio (JS / MLE)", fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1.1)
ax.grid(True, alpha=0.3)
# (b) 縮小推定の可視化（p=10の例）
ax = axes[1]
p = 10
theta_true = np.random.normal(3, 2, p)
y = np.random.normal(theta_true, 1)
# ジェームズ・スタイン推定量
shrinkage = max(0, 1 - (p - 2) / np.sum(y**2))
theta_js = shrinkage * y
for i in range(p):
ax.plot([i, i], [y[i], theta_js[i]], "gray", linewidth=1, alpha=0.7)
ax.scatter(range(p), y, s=80, color="red", marker="^", zorder=5,
label=r"MLE $\hat{\theta}_i = y_i$")
ax.scatter(range(p), theta_js, s=80, color="blue", marker="o", zorder=5,
label=rf"JS (shrinkage={shrinkage:.3f})")
ax.scatter(range(p), theta_true, s=40, color="green", marker="D", zorder=5,
label=r"True $\theta_i$")
ax.axhline(0, color="purple", linestyle="--", linewidth=1.5, alpha=0.5,
label="Shrinkage target (0)")
ax.set_xlabel("Parameter index i", fontsize=12)
ax.set_ylabel("Value", fontsize=12)
ax.set_title(f"Shrinkage Estimation (p = {p})", fontsize=13)
ax.legend(fontsize=9, loc="upper right")
ax.grid(True, alpha=0.3)
mse_mle = np.mean((y - theta_true)**2)
mse_js = np.mean((theta_js - theta_true)**2)
ax.text(0.02, 0.02, f"MSE(MLE)={mse_mle:.2f}, MSE(JS)={mse_js:.2f}",
transform=ax.transAxes, fontsize=10,
bbox=dict(boxstyle="round,pad=0.3", facecolor="lightyellow"))
plt.tight_layout()
plt.savefig("james_stein.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、ジェームズ・スタイン推定量の性質が確認できます。
-
左図（リスク比）: ジェームズ・スタイン推定量のリスクは、あらゆる $\|\theta\|^2$ に対して最尤推定量のリスクを下回っています（比が1未満）。$\|\theta\|^2$ が小さいとき（パラメータが原点に近いとき）に改善が最も大きく、$\|\theta\|^2$ が大きくなると改善は小さくなりますが、依然として1未満を保っています。また、次元 $p$ が大きいほど改善が顕著です
-
右図（縮小推定の可視化）: MLE（赤い三角）が原点に向かって縮小されてJS推定量（青い丸）になっています。灰色の線が縮小の量を示しています。真の値（緑のダイヤ）とのMSEを比較すると、JS推定量の方が改善されていることが確認できます
ジェームズ・スタイン推定量は原点への縮小しか行いませんでしたが、一般の中心への縮小や、より柔軟な事前分布を考えることで、パラメトリック経験ベイズ法の一般的な理論が構築されます。
パラメトリック経験ベイズ法
一般的な定式化
パラメトリック経験ベイズ法では、事前分布のパラメトリックな形を仮定し、そのハイパーパラメータをデータから推定します。
モデル:
$$ \begin{align} \theta_i &\sim g(\theta | \eta), \quad i = 1, \ldots, p \\ y_i | \theta_i &\sim f(y | \theta_i) \end{align} $$
周辺尤度（marginal likelihood、evidence）:
$$ \begin{equation} m(y | \eta) = \prod_{i=1}^{p} m_i(y_i | \eta), \quad m_i(y_i | \eta) = \int f(y_i | \theta) g(\theta | \eta) d\theta \end{equation} $$
ハイパーパラメータの推定: 周辺尤度を最大化して $\hat{\eta}$ を求める
$$ \begin{equation} \hat{\eta} = \arg\max_{\eta} \sum_{i=1}^{p} \ln m_i(y_i | \eta) \end{equation} $$
事後推定: $\hat{\eta}$ を代入して各 $\theta_i$ の事後分布を計算
$$ p(\theta_i | y_i, \hat{\eta}) = \frac{f(y_i | \theta_i) g(\theta_i | \hat{\eta})}{m_i(y_i | \hat{\eta})} $$
正規-正規モデルの例
最もよく使われるパラメトリック経験ベイズの例は、正規-正規モデルです。
$$ \theta_i \sim N(\mu, A), \quad y_i | \theta_i \sim N(\theta_i, \sigma_i^2) $$
ここで $\sigma_i^2$ は既知（あるいは十分なデータから正確に推定可能）とします。
$y_i$ の周辺分布は $y_i \sim N(\mu, A + \sigma_i^2)$ なので、周辺尤度の最大化は
$$ \hat{\mu}, \hat{A} = \arg\max_{\mu, A} \sum_{i=1}^{p} \left[-\frac{1}{2}\ln(A + \sigma_i^2) – \frac{(y_i – \mu)^2}{2(A + \sigma_i^2)}\right] $$
$\sigma_i^2$ が全て等しい場合（$\sigma_i^2 = \sigma^2$）:
$$ \hat{\mu} = \bar{y}, \quad \hat{A} = \max\left(0, \frac{1}{p}\sum_{i=1}^{p}(y_i – \bar{y})^2 – \sigma^2\right) $$
事後平均は
$$ \begin{equation} \hat{\theta}_i^{\text{EB}} = \hat{B}_i \hat{\mu} + (1 – \hat{B}_i) y_i, \quad \hat{B}_i = \frac{\sigma_i^2}{\hat{A} + \sigma_i^2} \end{equation} $$
$\hat{B}_i$ は縮小係数（shrinkage factor）で、データの精度が低いほど（$\sigma_i^2$ が大きいほど）全体平均 $\hat{\mu}$ への縮小が強くなります。
階層ベイズとの違い
経験ベイズと階層ベイズの本質的な違いは、ハイパーパラメータ $\eta$ の扱いにあります。
| 経験ベイズ | 階層ベイズ | |
|---|---|---|
| $\eta$ の扱い | 点推定 $\hat{\eta}$（周辺尤度最大化） | 事前分布 $p(\eta)$ を置いて事後分布を計算 |
| $\eta$ の不確実性 | 無視 | 考慮される |
| 計算コスト | 低い（最適化 + 事後計算） | 高い（MCMCなど） |
| 信頼区間 | $\eta$ の不確実性を無視するため、やや過小評価 | 正確 |
| 適用場面 | $p$ が大きい（パラメータが多い）とき特に有効 | $p$ が小さくても有効 |
$p$ が大きいとき、$\hat{\eta}$ の推定精度が高くなるため、経験ベイズと階層ベイズの差は小さくなります。
ノンパラメトリック経験ベイズ法
Robbinsの方法
ハーバート・ロビンズ（Robbins, 1956）は、事前分布 $g(\theta)$ の形を仮定せずにベイズ推定を行う方法を提案しました。これがノンパラメトリック経験ベイズ法の原点です。
ポアソンモデルを例に考えましょう。
$$ \theta_i \sim g(\theta), \quad y_i | \theta_i \sim \text{Poisson}(\theta_i) $$
$\theta_i$ の事後平均（二乗損失のもとでのベイズ推定量）は
$$ E[\theta_i | y_i] = (y_i + 1) \frac{f(y_i + 1)}{f(y_i)} $$
ここで $f(y) = \int \frac{e^{-\theta}\theta^y}{y!} g(\theta) d\theta$ は $y$ の周辺確率質量関数です。
驚くべきことに、事後平均を計算するのに事前分布 $g(\theta)$ を推定する必要がないのです。必要なのは周辺分布 $f(y)$ だけであり、これはデータの経験度数 $f(y) \approx n_y / p$（$y$ という値が観測された回数 $n_y$ を全体数 $p$ で割ったもの）で推定できます。
$$ \begin{equation} \hat{\theta}_i^{\text{EB}} = (y_i + 1) \frac{\hat{f}(y_i + 1)}{\hat{f}(y_i)} = (y_i + 1) \frac{n_{y_i + 1}}{n_{y_i}} \end{equation} $$
Robbins推定量の直感
この推定量の直感的な意味を考えてみましょう。
たとえば、$y_i = 3$ を観測したとき、ノンパラメトリック経験ベイズ推定量は $4 \times n_4 / n_3$ です。
- $n_4 / n_3 \approx 1$ ならば $\hat{\theta}_i \approx 4$: 観測値 $y_i = 3$ はほぼそのまま使われる
- $n_4 / n_3 < 1$ ならば $\hat{\theta}_i < 4$: $y = 4$ のケースが相対的に少ないことから、$\theta$ はやや小さめと推定される
- $n_4 / n_3 > 1$ ならば $\hat{\theta}_i > 4$: $y = 4$ のケースが多いことから、$\theta$ はやや大きめと推定される
周辺分布の「傾き」が事後推定に影響を与えるのです。
正規モデルでのノンパラメトリック経験ベイズ
正規モデル $y_i | \theta_i \sim N(\theta_i, 1)$ の場合、ツイーディの公式（Tweedie’s formula）が対応する結果を与えます。
$$ \begin{equation} E[\theta | y] = y + \frac{d}{dy}\ln f(y) = y + \frac{f'(y)}{f(y)} \end{equation} $$
ここで $f(y) = \int \phi(y – \theta) g(\theta) d\theta$ は $y$ の周辺密度（$\phi$ は標準正規密度）です。
事後平均は、観測値 $y$ に周辺密度の対数微分を足したものです。$f'(y)/f(y)$ は周辺分布のスコア関数であり、密度が増加する方向（$f'(y) > 0$）に推定値をシフトさせます。
実用的には、$f(y)$ をカーネル密度推定やスプラインで推定し、その微分を計算します。
多重検定における経験ベイズ
大規模同時推論
経験ベイズ法が現代統計学で最も影響力を発揮しているのが、大規模同時推論（large-scale simultaneous inference）の分野です。
ゲノム研究では、数千から数万の遺伝子について同時に仮説検定を行います。各遺伝子 $i$ について「この遺伝子は疾患と関連があるか」を検定し、$p$ 値 $p_i$ を計算します。
個別の検定なら有意水準 $\alpha = 0.05$ で判断すればよいのですが、$10{,}000$ 個の検定を同時に行えば、全ての帰無仮説が正しくても約500個が「有意」になってしまいます。これが多重検定問題です。
局所偽発見率（local fdr）
ブラッドリー・エフロンは、経験ベイズの枠組みで多重検定問題を次のように定式化しました。
各遺伝子 $i$ について、帰無仮説 $H_0$ のもとでは検定統計量 $z_i$ が $N(0, 1)$ に従うとします。帰無仮説が正しい確率を $\pi_0$、対立仮説が正しいときの $z_i$ の分布を $f_1(z)$ とすると
$$ f(z) = \pi_0 \phi(z) + (1 – \pi_0) f_1(z) $$
各遺伝子の局所偽発見率（local false discovery rate, lfdr）は、観測された $z_i$ のもとで帰無仮説が正しい事後確率として定義されます。
$$ \begin{equation} \text{lfdr}(z) = P(H_0 | z) = \frac{\pi_0 \phi(z)}{f(z)} \end{equation} $$
$f(z)$ は全遺伝子のデータのヒストグラムから推定できます。$\pi_0$ も $f(z)$ の裾の振る舞いから推定できます。つまり、事前分布（$\pi_0$ と $f_1$）をデータから推定しているのであり、これはまさに経験ベイズの考え方です。
Benjamini-Hochberg法との関係
有名なBenjamini-Hochberg（BH）法による偽発見率（FDR）制御も、実は経験ベイズ的な解釈が可能です。
BH法の棄却規則は、$p$ 値を昇順に並べ $p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}$ として
$$ \text{Reject } H_{(i)} \text{ for } i \leq k^*, \quad k^* = \max\left\{i : p_{(i)} \leq \frac{i}{m}\alpha\right\} $$
この手続きにおいて、$p_{(i)}$ に対応する閾値 $i\alpha/m$ は、帰無仮説の割合 $\pi_0$ が1であるという保守的な仮定の下でのFDR制御を保証します。
エフロンの経験ベイズ法では $\pi_0$ をデータから推定するため、BH法よりも検出力（power）を高めることができます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# --- 大規模検定のシミュレーション ---
m = 5000 # 遺伝子数
pi0 = 0.9 # 帰無仮説が正しい割合（90%はノイズ）
n_null = int(m * pi0)
n_alt = m - n_null
# 検定統計量の生成
z_null = np.random.normal(0, 1, n_null) # 帰無仮説（効果なし）
z_alt = np.random.normal(2.5, 1.0, n_alt) # 対立仮説（効果あり）
z = np.concatenate([z_null, z_alt])
is_null = np.concatenate([np.ones(n_null), np.zeros(n_alt)]).astype(bool)
# シャッフル
idx = np.random.permutation(m)
z = z[idx]
is_null = is_null[idx]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) z値のヒストグラムと混合分布の推定
ax = axes[0]
z_grid = np.linspace(-4, 6, 300)
# 理論的な混合密度
f_true = pi0 * stats.norm.pdf(z_grid, 0, 1) + (1-pi0) * stats.norm.pdf(z_grid, 2.5, 1)
f_null = stats.norm.pdf(z_grid, 0, 1)
ax.hist(z, bins=80, density=True, color="lightblue", alpha=0.7,
edgecolor="gray", linewidth=0.5)
ax.plot(z_grid, f_true, "b-", linewidth=2, label="Mixture $f(z)$")
ax.plot(z_grid, pi0 * f_null, "r--", linewidth=2, label=r"$\pi_0 \phi(z)$")
ax.set_xlabel("z-score", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("Distribution of Test Statistics", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) 局所偽発見率（lfdr）
ax = axes[1]
# カーネル密度推定で f(z) を推定
from scipy.stats import gaussian_kde
kde = gaussian_kde(z, bw_method=0.15)
f_hat = kde(z_grid)
# pi0 の推定（中央部分の密度から）
# z ~ 0 付近では f(z) ≈ pi0 * phi(0) なので
pi0_hat = min(1.0, kde(0)[0] / stats.norm.pdf(0))
# lfdr の計算
lfdr = pi0_hat * stats.norm.pdf(z_grid) / np.maximum(f_hat, 1e-10)
lfdr = np.clip(lfdr, 0, 1)
ax.plot(z_grid, lfdr, "b-", linewidth=2)
ax.axhline(0.2, color="red", linestyle="--", linewidth=1.5, alpha=0.7,
label="lfdr = 0.2 threshold")
ax.set_xlabel("z-score", fontsize=12)
ax.set_ylabel("Local FDR", fontsize=12)
ax.set_title(rf"Empirical Bayes lfdr ($\hat{{\pi}}_0 = {pi0_hat:.3f}$)", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(-0.05, 1.05)
# (c) BH法との比較
ax = axes[2]
# p値の計算
p_values = 2 * (1 - stats.norm.cdf(np.abs(z)))
# BH法
alpha_fdr = 0.1
sorted_p = np.sort(p_values)
bh_threshold = np.arange(1, m+1) * alpha_fdr / m
bh_reject = sorted_p <= bh_threshold
k_star = np.max(np.where(bh_reject)[0]) + 1 if np.any(bh_reject) else 0
bh_cutoff = sorted_p[k_star - 1] if k_star > 0 else 0
bh_rejected = p_values <= bh_cutoff
# 経験ベイズ（lfdr < 0.2 で棄却）
lfdr_per_point = pi0_hat * stats.norm.pdf(z) / np.maximum(kde(z), 1e-10)
lfdr_per_point = np.clip(lfdr_per_point, 0, 1)
eb_rejected = lfdr_per_point < 0.2
# 結果の比較
categories = ["BH method", "Empirical Bayes"]
tp = [np.sum(bh_rejected & ~is_null), np.sum(eb_rejected & ~is_null)]
fp = [np.sum(bh_rejected & is_null), np.sum(eb_rejected & is_null)]
fn = [np.sum(~bh_rejected & ~is_null), np.sum(~eb_rejected & ~is_null)]
x_pos = np.arange(len(categories))
width = 0.25
ax.bar(x_pos - width, tp, width, color="green", alpha=0.7, label="True Positive")
ax.bar(x_pos, fp, width, color="red", alpha=0.7, label="False Positive")
ax.bar(x_pos + width, fn, width, color="orange", alpha=0.7, label="False Negative")
for j in range(len(categories)):
fdr_actual = fp[j] / max(1, tp[j] + fp[j])
ax.text(j, max(tp[j], fp[j], fn[j]) + 10, f"FDR={fdr_actual:.3f}",
ha="center", fontsize=10, fontweight="bold")
ax.set_ylabel("Count", fontsize=12)
ax.set_title("BH vs Empirical Bayes", fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels(categories, fontsize=11)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("empirical_bayes_testing.png", dpi=150, bbox_inches="tight")
plt.show()
このシミュレーション結果から、経験ベイズ法の多重検定への適用効果が読み取れます。
-
左図（検定統計量の分布）: 5,000個の遺伝子のz値のヒストグラムは、帰無仮説の標準正規分布（赤い破線、$\pi_0\phi(z)$）と対立仮説の分布が混合した形（青い実線）になっています。中央付近では帰無仮説の分布が支配的ですが、右側（$z > 2$ 付近）では対立仮説からのシグナルが見えています
-
中央図（局所偽発見率）: 経験ベイズで推定した局所偽発見率（lfdr）は、z値が大きくなるほど小さくなります。$z \approx 2$ 付近でlfdrが急速に減少し、閾値0.2を下回ります。$\hat{\pi}_0$ が真の値0.9に近い値に推定されていることも確認できます
-
右図（BH法との比較）: 経験ベイズ法はBH法と比較して、同程度のFDR制御を保ちながら、より多くの真の陽性（True Positive）を検出し、偽陰性（False Negative）を減らしています。これは $\pi_0$ の推定による検出力の改善効果です
経験ベイズ法の応用例：打率の推定
野球の打率推定
冒頭で触れた打率推定問題を、パラメトリック経験ベイズで解いてみましょう。
各打者 $i$ の真の打率を $\theta_i$、打数を $n_i$、安打数を $y_i$ とします。二項-ベータモデルを設定します。
$$ \theta_i \sim \text{Beta}(\alpha, \beta), \quad y_i | \theta_i \sim \text{Binomial}(n_i, \theta_i) $$
$y_i$ の周辺分布はベータ-二項分布 $\text{BetaBinomial}(n_i, \alpha, \beta)$ であり、$\alpha$ と $\beta$ を周辺尤度の最大化で推定します。
事後分布は $\theta_i | y_i \sim \text{Beta}(\alpha + y_i, \beta + n_i – y_i)$ であり、事後平均は
$$ \begin{equation} \hat{\theta}_i^{\text{EB}} = \frac{\hat{\alpha} + y_i}{\hat{\alpha} + \hat{\beta} + n_i} \end{equation} $$
打数 $n_i$ が少ないとき、$\hat{\alpha}/(\hat{\alpha} + \hat{\beta})$（事前平均、すなわちリーグ平均打率）に強く引き寄せられ、$n_i$ が多いとき、観測打率 $y_i/n_i$ に近づきます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize_scalar, minimize
np.random.seed(42)
# --- 模擬的な打率データ ---
n_players = 200
# 真の打率はベータ分布から生成（リーグ平均 .265、ばらつきあり）
alpha_true, beta_true = 40, 111 # mean = 40/151 ≈ .265
theta_true = np.random.beta(alpha_true, beta_true, n_players)
# 打数（シーズン前半なので50-300打席）
n_at_bats = np.random.randint(50, 300, n_players)
# 安打数
hits = np.array([np.random.binomial(n, theta) for n, theta in zip(n_at_bats, theta_true)])
# 観測打率
batting_avg = hits / n_at_bats
# --- ハイパーパラメータの推定（周辺尤度最大化） ---
def neg_log_marginal_likelihood(params):
a, b = params
if a <= 0 or b <= 0:
return 1e10
log_ml = 0
for y, n in zip(hits, n_at_bats):
# log Beta-Binomial PMF
from scipy.special import gammaln, comb
log_ml += (gammaln(n+1) - gammaln(y+1) - gammaln(n-y+1)
+ gammaln(a+b) - gammaln(a) - gammaln(b)
+ gammaln(a+y) + gammaln(b+n-y) - gammaln(a+b+n))
return -log_ml
result = minimize(neg_log_marginal_likelihood, [30, 100], method="Nelder-Mead")
alpha_hat, beta_hat = result.x
print(f"Estimated hyperparameters: alpha={alpha_hat:.2f}, beta={beta_hat:.2f}")
print(f"Prior mean (league avg): {alpha_hat/(alpha_hat+beta_hat):.3f}")
print(f"True prior mean: {alpha_true/(alpha_true+beta_true):.3f}")
# --- 経験ベイズ推定 ---
eb_estimate = (alpha_hat + hits) / (alpha_hat + beta_hat + n_at_bats)
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 縮小推定の効果
ax = axes[0]
ax.scatter(batting_avg, eb_estimate, c=n_at_bats, cmap="viridis",
alpha=0.6, s=30, edgecolors="gray", linewidth=0.3)
cbar = plt.colorbar(ax.collections[0], ax=ax)
cbar.set_label("At-bats", fontsize=11)
# 45度線
lim = [0.1, 0.45]
ax.plot(lim, lim, "r--", linewidth=1.5, alpha=0.5, label="No shrinkage")
league_avg = alpha_hat / (alpha_hat + beta_hat)
ax.axhline(league_avg, color="green", linestyle=":", linewidth=1.5, alpha=0.5,
label=f"League avg = {league_avg:.3f}")
ax.set_xlabel("Observed batting avg (MLE)", fontsize=12)
ax.set_ylabel("Empirical Bayes estimate", fontsize=12)
ax.set_title("Shrinkage toward League Average", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_aspect("equal")
# (b) MSEの比較
ax = axes[1]
mse_mle = (batting_avg - theta_true)**2
mse_eb = (eb_estimate - theta_true)**2
# 打数でビン分け
bins = [50, 100, 150, 200, 250, 300]
bin_labels = ["50-99", "100-149", "150-199", "200-249", "250-300"]
mse_mle_binned = []
mse_eb_binned = []
for low, high in zip(bins[:-1], bins[1:]):
mask = (n_at_bats >= low) & (n_at_bats < high)
mse_mle_binned.append(np.mean(mse_mle[mask]))
mse_eb_binned.append(np.mean(mse_eb[mask]))
x_pos = np.arange(len(bin_labels))
width = 0.35
ax.bar(x_pos - width/2, mse_mle_binned, width, color="salmon", alpha=0.8,
label="MLE")
ax.bar(x_pos + width/2, mse_eb_binned, width, color="steelblue", alpha=0.8,
label="Empirical Bayes")
ax.set_xlabel("At-bats range", fontsize=12)
ax.set_ylabel("Mean Squared Error", fontsize=12)
ax.set_title("MSE by Number of At-bats", fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels(bin_labels, fontsize=9)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")
# (c) 事前分布の推定
ax = axes[2]
theta_grid = np.linspace(0.1, 0.45, 200)
prior_true = stats.beta.pdf(theta_grid, alpha_true, beta_true)
prior_estimated = stats.beta.pdf(theta_grid, alpha_hat, beta_hat)
ax.hist(theta_true, bins=30, density=True, color="lightblue", alpha=0.7,
edgecolor="gray", linewidth=0.5, label="True abilities (histogram)")
ax.plot(theta_grid, prior_true, "r-", linewidth=2, label="True prior")
ax.plot(theta_grid, prior_estimated, "b--", linewidth=2, label="Estimated prior (EB)")
ax.set_xlabel(r"True batting average $\theta$", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("Prior Distribution Estimation", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("empirical_bayes_batting.png", dpi=150, bbox_inches="tight")
plt.show()
# 全体のMSE
print(f"\nOverall MSE(MLE): {np.mean(mse_mle):.6f}")
print(f"Overall MSE(EB): {np.mean(mse_eb):.6f}")
print(f"Improvement: {(1 - np.mean(mse_eb)/np.mean(mse_mle))*100:.1f}%")
このグラフから、打率推定における経験ベイズ法の効果が読み取れます。
-
左図（縮小推定）: 各点が1人の打者を表しています。観測打率（横軸）に対して経験ベイズ推定値（縦軸）が45度線より内側（リーグ平均に向かって）に縮小されています。色はデータ量（打数）を表しており、打数が少ない選手（暗い色）ほど縮小が大きく、打数が多い選手（明るい色）ほど観測値に近い推定値を取っています。これは「データが少ないほど全体の情報を借りる」という経験ベイズの本質を視覚的に示しています
-
中央図（MSEの比較）: 打数の範囲ごとにMSEを比較すると、全ての打数範囲で経験ベイズ（青）がMLE（赤）を上回っています。特に打数が少ないグループ（50-99打数）で改善が最も大きく、打数が増えるにつれて差は小さくなります。これはデータが多い場合にはMLE自体が十分正確であるためです
-
右図（事前分布の推定）: 経験ベイズで推定した事前分布（青い破線）が真の事前分布（赤い実線）に非常に近いことが確認できます。真の打率の分布（ヒストグラム）と推定された事前分布が一致していることから、周辺尤度最大化によるハイパーパラメータ推定が適切に機能していることがわかります
経験ベイズの理論的保証
漸近的最適性
経験ベイズ推定量の理論的保証は、パラメータの数 $p \to \infty$ のもとで得られます。
パラメトリック経験ベイズの場合、事前分布のパラメトリックな形が正しいとき（well-specified case）、$\hat{\eta}$ は $\eta$ の一致推定量となり、経験ベイズ推定量のリスクはオラクルベイズ推定量（$\eta$ が既知のときのベイズ推定量）のリスクに収束します。
$$ \frac{1}{p}\sum_{i=1}^{p} R(\theta_i, \hat{\theta}_i^{\text{EB}}) \to \frac{1}{p}\sum_{i=1}^{p} R(\theta_i, \hat{\theta}_i^{\text{Oracle}}) \quad (p \to \infty) $$
ノンパラメトリック経験ベイズの場合でも、Robbinsの推定量は事前分布の形によらず漸近的にベイズ最適に近づくことが示されています。
有限サンプルでの注意点
経験ベイズの有限サンプルでの問題点も知っておく必要があります。
$p$ が小さい場合: ハイパーパラメータの推定が不安定になり、経験ベイズの利点が薄れます。目安として $p \geq 20$ 程度は欲しいところです。$p$ が小さい場合は完全ベイズ（階層ベイズ）の方が適切です。
不確実性の過小評価: ハイパーパラメータを点推定で固定するため、信頼区間が実際よりも狭くなる傾向があります。これを修正するためのパラメトリックブートストラップなどの方法が提案されています。
事前分布の誤特定: パラメトリック経験ベイズでは、事前分布の形を仮定するため、この仮定が間違っていると推定にバイアスが生じます。心配な場合はノンパラメトリック経験ベイズや完全ベイズを検討すべきです。
まとめ
本記事では、経験ベイズ法の理論と応用について解説しました。
- 経験ベイズ法はデータからハイパーパラメータ（事前分布）を推定し、個別の推定を改善する方法である
- ジェームズ・スタイン推定量は経験ベイズの原点であり、3次元以上で最尤推定量を一様に上回る「スタインのパラドクス」を実現する
- パラメトリック経験ベイズは事前分布の形を仮定し、周辺尤度最大化でハイパーパラメータを推定する。正規-正規モデルやベータ-二項モデルが代表例
- ノンパラメトリック経験ベイズ（Robbins, 1956）は事前分布の形を仮定せず、周辺分布の推定だけでベイズ推定を行う
- 多重検定において、経験ベイズは局所偽発見率（lfdr）の推定を通じて、効率的なFDR制御を実現する
- 縮小推定の効果はデータが少ないパラメータほど大きく、データが多いパラメータではMLEに近づく（適応的な情報の借り入れ）
次のステップとして、以下の記事も参考にしてください。
- 階層ベイズモデル — 経験ベイズの完全ベイズ版
- ベイズモデル比較 — ハイパーパラメータ選択の別の視点
- 偽発見率の制御 — BH法と経験ベイズの関係
- ラオ・ブラックウェルの定理 — 推定量改善の理論的基盤