機械学習と情報技術

ベイズ仮説検定とベイズファクターの計算

ベイズ統計

「この薬は効くのか」「コインは公正か」——仮説検定はデータに基づいてこうした問いに答える手法です。頻度論的な検定ではp値を計算して帰無仮説を棄却するか判断しますが、ベイズ的なアプローチでは根本的に異なる問いを立てます。「帰無仮説の事後確率はどれだけか」あるいは「データはどちらの仮説をどれだけ支持しているか」を直接計算するのです。

この「どちらの仮説をどれだけ支持するか」の定量化がベイズファクター（Bayes factor）です。ベイズファクターはp値が抱える多くの問題を回避し、仮説間の証拠の強さを直感的に評価できる道具として注目されています。

ベイズ仮説検定を理解すると、以下のことが可能になります。

  • 仮説の相対的な評価: 帰無仮説だけでなく対立仮説の確率も定量化
  • p値の限界の克服: サンプルサイズに依存しない証拠の強さの評価
  • 逐次的な検定: データが追加されるたびに証拠を更新
  • モデル比較: 2つ以上のモデルの相対的な適合度を評価

本記事の内容

  • ベイズ仮説検定の基本的な枠組み
  • ベイズファクターの定義と計算
  • 頻度論的p値との本質的な違い
  • リンドレーのパラドックス
  • 事前分布の感度分析
  • Pythonによる具体例の実装

前提知識

ベイズ仮説検定の枠組み

頻度論的検定の限界

ベイズ仮説検定の動機を理解するために、まず頻度論的検定の問題点を整理しましょう。

頻度論的検定では、帰無仮説 $H_0$ のもとでのp値——「$H_0$ が正しいと仮定したとき、観測されたデータ以上に極端なデータが得られる確率」——を計算し、それが小さければ（通常は0.05以下）$H_0$ を棄却します。

しかし、p値にはいくつかの根本的な問題があります。

p値は $H_0$ の確率ではない: p値は「$H_0$ が正しいならばデータが極端になる確率」であり、「データが得られたときに $H_0$ が正しい確率」ではありません。多くの研究者がこの違いを混同しています。「p = 0.03 なので $H_0$ が正しい確率は3%」という解釈は誤りです。

サンプルサイズ依存性: サンプルサイズを十分に大きくすれば、実質的に無意味な差でもp値は有意になります。逆に、サンプルサイズが小さければ、大きな効果があってもp値が有意にならないことがあります。

$H_1$ について何も言わない: p値は帰無仮説のもとでのみ計算されるため、対立仮説 $H_1$ がデータをどの程度うまく説明するかについての情報を含んでいません。

オプショナルストッピング問題: 頻度論的検定では、「有意になるまでデータを集め続ける」と偽陽性率が膨張します。実験の停止ルールが結果の解釈に影響するのです。

ベイズ仮説検定は、これらの問題に対して根本的に異なるアプローチを提供します。

仮説に事前確率を置く

ベイズ的アプローチでは、仮説自体に事前確率を割り当てます。

$$ \pi_0 = P(H_0), \quad \pi_1 = P(H_1) = 1 – \pi_0 $$

事前確率 $\pi_0$ は、データを見る前の段階で $H_0$ がどの程度もっともらしいかについての信念を表します。特別な事前情報がなければ $\pi_0 = \pi_1 = 0.5$（五分五分）とすることが多いです。

データ $\bm{x}$ を観測した後、ベイズの定理により仮説の事後確率が計算できます。

$$ \begin{equation} P(H_0 | \bm{x}) = \frac{P(\bm{x} | H_0)\,\pi_0}{P(\bm{x} | H_0)\,\pi_0 + P(\bm{x} | H_1)\,\pi_1} \end{equation} $$

ここで $P(\bm{x} | H_i)$ は仮説 $H_i$ のもとでのデータの周辺尤度（marginal likelihood）です。

この式が意味するのは、仮説の事後確率がデータの尤度と事前確率の両方によって決まるということです。頻度論的検定ではデータの尤度（p値に相当する部分）のみを見ますが、ベイズ検定では事前確率も考慮に入れます。

周辺尤度の計算

$H_i$ がパラメータ $\theta$ を含む場合、周辺尤度は事前分布 $\pi_i(\theta)$ で積分した値です。

$$ \begin{equation} P(\bm{x} | H_i) = \int P(\bm{x} | \theta, H_i)\,\pi_i(\theta)\,d\theta \end{equation} $$

この周辺尤度は「事前分布の予測の良さ」を表しています。事前分布が予測したパラメータの範囲にデータが整合的であれば大きく、事前予測と整合的でなければ小さくなります。

周辺尤度の計算は一般に容易ではありません。解析的に計算できるのは共役事前分布を使った場合に限られ、それ以外はモンテカルロ法などの数値的手法が必要になります。共役事前分布を使った場合の計算例は後ほど具体的に示します。

次に、2つの仮説の相対的な支持度を定量化するベイズファクターを定義しましょう。

ベイズファクターの定義と解釈

定義

ベイズファクター（Bayes factor）は、2つの仮説の周辺尤度の比として定義されます。

$$ \begin{equation} \text{BF}_{10} = \frac{P(\bm{x} | H_1)}{P(\bm{x} | H_0)} \end{equation} $$

$\text{BF}_{10} > 1$ なら $H_1$ がデータによってより支持され、$\text{BF}_{10} < 1$ なら $H_0$ がより支持されます。添字の「10」は「$H_1$ 対 $H_0$ の比」を意味します。逆のベイズファクターは $\text{BF}_{01} = 1/\text{BF}_{10}$ です。

事後オッズとの関係

事後オッズと事前オッズの関係は次のように表されます。

$$ \begin{equation} \underbrace{\frac{P(H_1 | \bm{x})}{P(H_0 | \bm{x})}}_{\text{事後オッズ}} = \underbrace{\text{BF}_{10}}_{\text{ベイズファクター}} \times \underbrace{\frac{\pi_1}{\pi_0}}_{\text{事前オッズ}} \end{equation} $$

ベイズファクターは「データがオッズをどれだけ更新するか」を表す乗数です。事前オッズが1:1（五分五分）なら、事後オッズ = ベイズファクターとなり、ベイズファクターがそのまま仮説間の証拠の比率を示します。

この分解には重要な意味があります。ベイズファクターは事前確率 $\pi_0, \pi_1$ に依存しません。したがって、事前確率について合意できない研究者同士でも、ベイズファクターの値については合意できます。事前確率は主観的ですが、ベイズファクター（データが証拠をどう更新するか）は客観的な量なのです。

ベイズファクターの対数表現

ベイズファクターの対数（ベイズの重み、weight of evidence）を使うと、証拠の蓄積が加法的になるため扱いやすいことがあります。

$$ \ln \text{BF}_{10} = \ln P(\bm{x} | H_1) – \ln P(\bm{x} | H_0) $$

データが独立に $x_1, x_2, \dots, x_n$ と観測される場合、

$$ \ln \text{BF}_{10} = \sum_{i=1}^n \left[\ln P(x_i | H_1) – \ln P(x_i | H_0)\right] $$

と各データ点の寄与の和になります（正確には、逐次更新では事前が更新されるため、$P(x_i | \text{前のデータ}, H_k)$ としての周辺尤度を使います）。

ジェフリーズの基準

ベイズファクターの値を解釈するためのガイドラインとして、ジェフリーズ（Harold Jeffreys）の基準が広く使われています。

$\text{BF}_{10}$ $\ln \text{BF}_{10}$ 証拠の強さ
1 – 3 0 – 1.1 弱い（anecdotal）
3 – 10 1.1 – 2.3 中程度（moderate）
10 – 30 2.3 – 3.4 強い（strong）
30 – 100 3.4 – 4.6 非常に強い（very strong）
> 100 > 4.6 決定的（decisive）

ジェフリーズの基準はあくまで目安であり、分野や文脈に応じて解釈は異なります。物理学（5σ基準）では決定的な証拠を要求しますが、探索的な研究では中程度の証拠でも十分とされることがあります。

重要なのは、ベイズファクターは連続的な量であるということです。頻度論的検定のように「有意 or 非有意」の二分法ではなく、証拠の強さのグラデーションとして結果を解釈できます。

具体例: コインの公正さの検定

ベイズファクターの計算を具体例で見てみましょう。コインを100回投げて65回表が出たとき、このコインは公正（$p = 0.5$）でしょうか。

仮説の設定

  • $H_0$: $p = 0.5$（公正なコイン）— 単純仮説（パラメータが1点で指定される）
  • $H_1$: $p \neq 0.5$（不公正なコイン）— 複合仮説（パラメータに事前分布を置く）

$H_1$ のもとでの $p$ の事前分布として、一様分布 $p \sim \text{Beta}(1, 1)$（= $U(0, 1)$）を使います。これは「$p$ がどんな値でも等確率で起こりうる」という非情報的な事前分布です。

周辺尤度の計算

$H_0$ のもとでの周辺尤度: $p = 0.5$ は固定値なので積分は不要です。

$$ P(\text{data} | H_0) = \binom{100}{65} \cdot 0.5^{65} \cdot 0.5^{35} = \binom{100}{65} \cdot 0.5^{100} $$

$H_1$ のもとでの周辺尤度: Beta-Binomial モデルの周辺尤度が解析的に計算できます。

$$ P(\text{data} | H_1) = \int_0^1 \binom{100}{65} p^{65}(1-p)^{35} \cdot 1\, dp = \binom{100}{65} \cdot B(66, 36) $$

ここで $B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$ はベータ関数です。$\text{Beta}(1, 1)$ の事前分布は $1/B(1,1) = 1$ なので、$\pi(p) = 1$ です。

ベータ関数の性質 $\int_0^1 p^{a-1}(1-p)^{b-1}\,dp = B(a, b)$ を使うと、

$$ P(\text{data} | H_1) = \binom{100}{65} \cdot B(66, 36) = \binom{100}{65} \cdot \frac{65!\cdot 35!}{101!} $$

ベイズファクターの計算

$$ \text{BF}_{10} = \frac{P(\text{data}|H_1)}{P(\text{data}|H_0)} = \frac{B(66, 36)}{0.5^{100}} $$

この値をPythonで数値的に確認し、逐次的な更新も可視化します。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, integrate
from scipy.special import beta as beta_func, comb

np.random.seed(42)

# コインを100回投げて65回表が出た
n = 100
k = 65

# H0: p = 0.5（公正なコイン）
# H1: p ≠ 0.5（p ~ Beta(1,1) = U(0,1)を事前分布に）

# H0のもとでの周辺尤度
loglik_H0 = stats.binom.logpmf(k, n, 0.5)
lik_H0 = np.exp(loglik_H0)

# H1のもとでの周辺尤度（Beta-Binomialモデル）
# P(x|H1) = C(n,k) * B(k+1, n-k+1) / B(1,1) = C(n,k) * B(k+1, n-k+1)
lik_H1 = comb(n, k) * beta_func(k + 1, n - k + 1)

BF_10 = lik_H1 / lik_H0
BF_01 = 1 / BF_10

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# (a) ベイズファクターの概念図
ax = axes[0]
p_range = np.linspace(0, 1, 500)

# 尤度関数
likelihood = stats.binom.pmf(k, n, p_range)
# H0の予測（p=0.5でのピンポイント）
ax.axvline(0.5, color="blue", linewidth=2, linestyle="--",
           label=rf"$H_0$: $p = 0.5$, $P(data|H_0) = {lik_H0:.6f}$")
ax.plot(p_range, likelihood, "r-", linewidth=2, label="Likelihood")
ax.fill_between(p_range, likelihood, alpha=0.2, color="red")
ax.plot(p_range, np.ones_like(p_range) * lik_H1, "g--", linewidth=2,
        label=rf"$P(data|H_1) = {lik_H1:.6f}$ (avg likelihood)")

ax.set_xlabel("p", fontsize=12)
ax.set_ylabel("Likelihood", fontsize=12)
ax.set_title(f"Bayes factor: BF10 = {BF_10:.2f}", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) 逐次的なベイズファクターの更新
ax = axes[1]
data = np.concatenate([np.ones(k), np.zeros(n - k)])
np.random.shuffle(data)

cumulative_bf = []
bf = 1.0  # 初期ベイズファクター

for i in range(1, n + 1):
    k_i = int(data[:i].sum())
    n_i = i
    lik_h0_i = stats.binom.pmf(k_i, n_i, 0.5)
    lik_h1_i = comb(n_i, k_i) * beta_func(k_i + 1, n_i - k_i + 1)
    bf_i = lik_h1_i / lik_h0_i if lik_h0_i > 0 else float('inf')
    cumulative_bf.append(bf_i)

ax.semilogy(range(1, n + 1), cumulative_bf, "b-", linewidth=2)
ax.axhline(1, color="gray", linewidth=1, linestyle="--")
ax.axhline(10, color="green", linewidth=1, linestyle=":", label="Strong evidence (BF=10)")
ax.axhline(1/10, color="red", linewidth=1, linestyle=":", label="Strong for H0 (BF=1/10)")

ax.fill_between(range(1, n+1), 1/3, 3, alpha=0.1, color="gray",
                label="Inconclusive (1/3 < BF < 3)")

ax.set_xlabel("Number of coin flips", fontsize=12)
ax.set_ylabel("Bayes Factor (BF10)", fontsize=12)
ax.set_title("Sequential Bayes factor update", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")

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

# p値との比較
p_value = 2 * (1 - stats.binom.cdf(k - 1, n, 0.5))
print(f"ベイズファクター BF10 = {BF_10:.4f}")
print(f"ベイズファクター BF01 = {BF_01:.4f}")
print(f"p値（両側）= {p_value:.6f}")
print(f"事後確率 P(H0|data) = {1/(1 + BF_10):.4f} (事前オッズ1:1)")

この結果から、以下のことが読み取れます。

  1. 左図: 尤度関数のピークが $p = 0.65$ 付近にあり、$H_0$ の値 $p = 0.5$ から離れている 。$H_1$ のもとでの平均尤度（緑の破線）と $H_0$ のもとでの尤度（青の破線でのピンポイント値）の比がベイズファクターです。$H_1$ は $p$ の値域全体にわたって尤度を平均化するため、尤度のピーク付近では $H_0$ よりも有利ですが、$p$ が0や1に近い領域では不利になります。この「事前分布による自動的なペナルティ」がベイズファクターの重要な特性です。

  2. 右図: データが蓄積されるにつれてベイズファクターが更新される 。逐次的に証拠が蓄積され、最終的にどちらの仮説が支持されるかが明確になります。灰色の帯は「証拠が弱い」領域を示しています。頻度論的検定とは異なり、データを見ながら検定を停止しても（オプショナルストッピング）ベイズファクターの解釈は正当です。これはベイズファクターの尤度原理に基づく性質であり、実験計画の柔軟性を大きく高めます。

p値とベイズファクターの体系的な比較

p値とベイズファクターは常に一致するわけではありません。両者の関係を系統的に比較してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.special import beta as beta_func, comb

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# (a) サンプルサイズに対するp値とBFの比較
ax = axes[0]
n_range = np.arange(20, 501, 5)

# 65%の表が出る場合
proportion = 0.65
p_values = []
bf_values = []

for n in n_range:
    k = int(n * proportion)
    # p値（両側）
    p_val = 2 * min(stats.binom.cdf(k, n, 0.5),
                    1 - stats.binom.cdf(k - 1, n, 0.5))
    p_values.append(min(p_val, 1))

    # ベイズファクター
    lik_h0 = stats.binom.pmf(k, n, 0.5)
    lik_h1 = comb(n, k, exact=True) * beta_func(k + 1, n - k + 1) if k <= n else 0
    bf = lik_h1 / lik_h0 if lik_h0 > 0 else float('inf')
    bf_values.append(bf)

ax.semilogy(n_range, p_values, "r-", linewidth=2, label="p-value")
ax.semilogy(n_range, [1/b if b > 0 else float('inf') for b in bf_values],
            "b-", linewidth=2, label=r"$1/\mathrm{BF}_{10}$ (evidence for $H_0$)")
ax.axhline(0.05, color="red", linewidth=1, linestyle=":", alpha=0.7, label="p = 0.05")
ax.axhline(1/3, color="blue", linewidth=1, linestyle=":", alpha=0.7, label="BF = 1/3")
ax.set_xlabel("Sample size n", fontsize=12)
ax.set_ylabel("Value (log scale)", fontsize=12)
ax.set_title(f"p-value vs BF for {proportion*100:.0f}% heads", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")

# (b) 効果量が小さい場合（51%の表）
ax = axes[1]
proportion_small = 0.51

p_values_small = []
bf_values_small = []

for n in n_range:
    k = int(n * proportion_small)
    p_val = 2 * min(stats.binom.cdf(k, n, 0.5),
                    1 - stats.binom.cdf(k - 1, n, 0.5))
    p_values_small.append(min(p_val, 1))

    lik_h0 = stats.binom.pmf(k, n, 0.5)
    lik_h1 = comb(n, k, exact=True) * beta_func(k + 1, n - k + 1) if k <= n else 0
    bf = lik_h1 / lik_h0 if lik_h0 > 0 else float('inf')
    bf_values_small.append(bf)

ax.semilogy(n_range, p_values_small, "r-", linewidth=2, label="p-value")
ax.semilogy(n_range, bf_values_small, "b-", linewidth=2, label=r"$\mathrm{BF}_{10}$")
ax.axhline(0.05, color="red", linewidth=1, linestyle=":", alpha=0.7, label="p = 0.05")
ax.axhline(3, color="blue", linewidth=1, linestyle=":", alpha=0.7, label="BF = 3")
ax.set_xlabel("Sample size n", fontsize=12)
ax.set_ylabel("Value (log scale)", fontsize=12)
ax.set_title(f"Small effect: {proportion_small*100:.0f}% heads", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")

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

この比較から、p値とベイズファクターの重要な違いが明確になります。

  1. 左図（効果量が大きい場合、65%）: p値もベイズファクターもサンプルサイズの増加とともに $H_0$ に反する証拠が蓄積される 。両者の方向性は一致していますが、ベイズファクターの方が穏やかに変化する傾向があります。これは事前分布による「ペナルティ」の効果です。

  2. 右図（効果量が小さい場合、51%）: 大きなサンプルサイズではp値は有意になるが、ベイズファクターは必ずしも $H_1$ を支持しない 。$n = 400$ 程度で p値が0.05を下回っても、ベイズファクターはまだ証拠が弱い（BF < 3）領域にとどまることがあります。これは、実質的に無意味な差をサンプルサイズの力で「有意」にしてしまうp値の問題と、ベイズファクターがモデルの複雑さを自動的にペナルティするオッカムの剃刀の効果を対比的に示しています。

リンドレーのパラドックス

パラドックスの内容

ベイズファクターとp値は常に同じ結論を導くわけではありません。リンドレーのパラドックス（Lindley’s paradox）は、p値が有意（小さい）であっても、ベイズファクターは $H_0$ を支持する場合があることを示しています。

これは特にサンプルサイズが大きく、効果量が小さい場合に起こります。具体的に見てみましょう。

$X_1, \dots, X_n \overset{\text{i.i.d.}}{\sim} N(\theta, 1)$ として、$H_0: \theta = 0$ vs $H_1: \theta \neq 0$（$\theta \sim N(0, \sigma_0^2)$）を検定します。

標本平均 $\bar{x} = c/\sqrt{n}$（$c$ は定数）のとき、

p値: $p = 2(1 – \Phi(c))$ はサンプルサイズ $n$ に依存せず、$c$ が十分大きければ常に有意です。

ベイズファクター: $n \to \infty$ で

$$ \text{BF}_{01} = \frac{P(\bar{x}|H_0)}{P(\bar{x}|H_1)} \approx \sqrt{1 + n\sigma_0^2} \cdot \exp\left(-\frac{c^2}{2}\cdot\frac{n\sigma_0^2}{1 + n\sigma_0^2}\right) $$

$n \to \infty$ で第一の因子 $\sqrt{1 + n\sigma_0^2}$ が支配的になり、$\text{BF}_{01} \to \infty$、つまり $H_0$ が支持されるのです。

パラドックスの解釈

この一見矛盾した結果は、以下のように理解できます。

$H_1$ の事前分布が広い場合（$\sigma_0^2$ が大きい場合）、$H_1$ は「$\theta$ が0から非常に離れた値をとる」可能性にも事前確率を配分しています。しかし、$\bar{x} = c/\sqrt{n}$ ということは、$n$ が大きいとき $\theta$ の推定値は0に非常に近いということです。$H_1$ の事前分布の大部分が「データと整合しないパラメータ空間」に割り当てられているため、周辺尤度が小さくなるのです。

この現象は「ベイズ統計は自動的にオッカムの剃刀を組み込んでいる」と解釈できます。複雑なモデル（広い事前分布を持つ $H_1$）は、データをうまく説明できるパラメータ値が少ない場合、単純なモデル（$H_0$）に対してペナルティを受けます。

実用的な教訓

リンドレーのパラドックスから得られる実用的な教訓は以下の通りです。

  1. 事前分布の選択が重要: $H_1$ のもとでの事前分布は、実質的に妥当な効果量の範囲を反映すべきです。あまりに広い事前分布は不当に $H_0$ を有利にします。
  2. p値とベイズファクターを相補的に使う: 両者が同じ結論を導く場合は確信を持ちやすく、異なる結論を導く場合は慎重な吟味が必要です。
  3. 効果量を報告する: 統計的有意性だけでなく、推定される効果量とその不確実性も報告すべきです。

事前分布の感度分析

ベイズファクターは事前分布の選択に依存するため、結論のロバスト性を確認する感度分析が重要です。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.special import beta as beta_func, comb

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

n = 100
k = 65

# (a) Beta事前分布のパラメータに対するBFの感度
ax = axes[0]
a_vals = np.linspace(0.5, 10, 100)

for b_prior in [0.5, 1.0, 2.0, 5.0]:
    bf_vals = []
    for a_prior in a_vals:
        lik_h0 = stats.binom.pmf(k, n, 0.5)
        # Beta(a, b)事前分布でのBeta-Binomial周辺尤度
        lik_h1 = comb(n, k) * beta_func(k + a_prior, n - k + b_prior) / beta_func(a_prior, b_prior)
        bf_vals.append(lik_h1 / lik_h0)
    ax.plot(a_vals, bf_vals, linewidth=2,
            label=rf"$\beta = {b_prior}$")

ax.set_xlabel(r"Prior parameter $\alpha$ in Beta($\alpha, \beta$)", fontsize=12)
ax.set_ylabel(r"$\mathrm{BF}_{10}$", fontsize=12)
ax.set_title(f"Sensitivity to prior (n={n}, k={k})", fontsize=13)
ax.legend(fontsize=10)
ax.axhline(1, color="gray", linewidth=1, linestyle="--")
ax.axhline(3, color="green", linewidth=1, linestyle=":", alpha=0.7)
ax.grid(True, alpha=0.3)

# (b) 正規分布の検定での事前分布の幅に対する感度（リンドレーのパラドックス）
ax = axes[1]
sigma0_range = np.linspace(0.1, 10, 200)

for z_stat in [2.0, 2.5, 3.0, 3.5]:
    bf01_vals = []
    for s0 in sigma0_range:
        # BF01 for z-test with N(0, sigma0^2) prior
        n_eff = 100  # サンプルサイズ
        bf01 = np.sqrt(1 + n_eff * s0**2) * np.exp(
            -0.5 * z_stat**2 * n_eff * s0**2 / (1 + n_eff * s0**2))
        bf01_vals.append(bf01)
    ax.semilogy(sigma0_range, bf01_vals, linewidth=2,
                label=rf"$z = {z_stat}$ (p={2*(1-stats.norm.cdf(z_stat)):.4f})")

ax.axhline(1, color="gray", linewidth=1, linestyle="--")
ax.axhline(1/3, color="red", linewidth=1, linestyle=":", alpha=0.7, label=r"BF$_{01}$=1/3")
ax.set_xlabel(r"Prior SD $\sigma_0$", fontsize=12)
ax.set_ylabel(r"$\mathrm{BF}_{01}$ (evidence for $H_0$)", fontsize=12)
ax.set_title(r"Lindley's paradox: BF$_{01}$ vs prior width (n=100)", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")

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

この感度分析から、以下のことが明らかになります。

  1. 左図: Beta事前分布のパラメータによってベイズファクターの値が変わるが、質的な結論は比較的ロバスト 。対称な事前分布（$\alpha = \beta$）は公正なコインの近くに確率質量を集めるため、BFが小さくなる傾向がありますが、一様事前分布（$\alpha = \beta = 1$）から大きく外れない限り、結論は大きく変わりません。

  2. 右図: リンドレーのパラドックスの可視化 。事前分布の幅 $\sigma_0$ が大きくなると、BF$_{01}$（$H_0$ の証拠）が増加します。$z = 2.0$（$p = 0.046$、ぎりぎり有意）の場合、$\sigma_0$ が適度に大きいとBFは $H_0$ を支持します。一方、$z = 3.5$（$p = 0.0005$）のように非常に強い証拠がある場合は、事前分布の幅に対してロバストに $H_1$ が支持されます。

ベイズファクターの実用的な計算法

サヴィッジ・ディッキー比

特定の条件下で、ベイズファクターを事前分布と事後分布の比として簡単に計算できる方法があります。帰無仮説が $H_0: \theta = \theta_0$ という点仮説であり、$H_1$ のもとでの事前分布が $\pi(\theta)$、事後分布が $\pi(\theta|\bm{x})$ であるとき、

$$ \text{BF}_{01} = \frac{\pi(\theta_0 | \bm{x})}{\pi(\theta_0)} $$

つまり、パラメータ値 $\theta_0$ における事後密度と事前密度の比がベイズファクターになります。これはサヴィッジ・ディッキー密度比（Savage-Dickey density ratio）と呼ばれ、MCMCのサンプルから事後密度を推定することでベイズファクターが近似的に計算できます。

この方法の利点は、$H_0$ のもとでのモデルを明示的にフィットする必要がなく、$H_1$ のもとでの事後分布さえ得られればよい点です。

ブリッジサンプリング

より一般的な状況では、ブリッジサンプリング（bridge sampling）やネスティド・サンプリング（nested sampling）といった手法が周辺尤度の数値計算に使われます。これらの方法はMCMCのサンプルから効率的に周辺尤度を推定できるため、Stan や PyMC といった確率プログラミング言語と組み合わせて広く使われています。

まとめ

本記事では、ベイズ仮説検定の枠組みとベイズファクターの理論を解説しました。

  • ベイズ仮説検定は仮説に事前確率を割り当て、データから事後確率を計算する枠組みであり、「仮説が正しい確率」を直接計算できる
  • ベイズファクター $\text{BF}_{10} = P(\bm{x}|H_1)/P(\bm{x}|H_0)$ は証拠の強さを定量化する乗数であり、事前確率に依存しない客観的な量である
  • ベイズファクターは逐次更新が可能であり、オプショナルストッピングの問題がない。これはp値にはない大きな利点である
  • リンドレーのパラドックスはp値とベイズファクターが異なる結論を導く場合があることを示し、特に大標本・小効果量の状況で顕著になる
  • ベイズアプローチはオッカムの剃刀を自動的に組み込んでおり、複雑なモデルはデータで支持されない限りペナルティを受ける
  • 事前分布の感度分析はベイズファクターを使う際に不可欠であり、結論のロバスト性を確認する必要がある

次のステップとして、以下の記事も参考にしてください。

関連記事

前の記事

確率の不等式 — マルコフ・チェビシェフ・ヘフディング・ベルンシュタイン

次の記事

完備十分統計量とレーマン・シェッフェの定理

関連タグ: p値 ベイズファクター ベイズ統計 モデル選択 事後確率 仮説検定