頻度論的な仮説検定では、p値を用いて帰無仮説を棄却するかどうかを判定します。しかし、p値は「帰無仮説が正しい場合に今回のデータ以上に極端な結果が得られる確率」であり、「帰無仮説が正しい確率」を直接教えてくれるわけではありません。
ベイズファクター（Bayes Factor） は、この限界を克服するベイズ統計学に基づく仮説比較の枠組みです。2つの仮説（モデル）の 周辺尤度（marginal likelihood） の比として定義され、データがどちらの仮説をどれだけ支持するかを直接的に定量化します。ベイズファクターは、帰無仮説を「支持する」方向の証拠も評価できるという点で、p値にはない対称性を持っています。
本記事の内容
- ベイズファクターの定義（周辺尤度の比）
- 事後オッズと事前オッズの関係
- ジェフリーズのスケールによる解釈
- 周辺尤度の計算方法（解析的計算とモンテカルロ法）
- 正規分布の平均検定でのベイズファクター
- p値との比較
- Pythonでの実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ベイズファクターとは
直感的な理解
2つのモデル（仮説）$M_1$ と $M_2$ があり、データ $\bm{x}$ がどちらのモデルをより強く支持するかを比較したいとします。
各モデルのもとで「データ $\bm{x}$ がどれくらいもっともらしいか」を測る指標が 周辺尤度（marginal likelihood） です。ベイズファクターは、この周辺尤度の比です。
$$ BF_{12} = \frac{P(\bm{x} \mid M_1)}{P(\bm{x} \mid M_2)} $$
$BF_{12} > 1$ であれば $M_1$ が $M_2$ よりもデータによって支持されており、$BF_{12} < 1$ であれば $M_2$ が支持されています。
周辺尤度の定義
モデル $M_j$ がパラメータ $\bm{\theta}_j$ を持つとき、周辺尤度は尤度を事前分布で積分したものです。
$$ P(\bm{x} \mid M_j) = \int P(\bm{x} \mid \bm{\theta}_j, M_j) \cdot \pi(\bm{\theta}_j \mid M_j) \, d\bm{\theta}_j $$
ここで、
- $P(\bm{x} \mid \bm{\theta}_j, M_j)$ はパラメータ $\bm{\theta}_j$ のもとでの尤度
- $\pi(\bm{\theta}_j \mid M_j)$ はパラメータの事前分布
周辺尤度は「パラメータの不確実性を事前分布で平均化した尤度」と解釈できます。これが「周辺」と呼ばれる理由は、パラメータを 周辺化（marginalize out） しているためです。
ベイズファクターの正式な定義
$$ BF_{12} = \frac{P(\bm{x} \mid M_1)}{P(\bm{x} \mid M_2)} = \frac{\displaystyle\int P(\bm{x} \mid \bm{\theta}_1, M_1) \pi(\bm{\theta}_1 \mid M_1) \, d\bm{\theta}_1}{\displaystyle\int P(\bm{x} \mid \bm{\theta}_2, M_2) \pi(\bm{\theta}_2 \mid M_2) \, d\bm{\theta}_2} $$
事後オッズと事前オッズの関係
ベイズの定理のモデル版
各モデルの事後確率はベイズの定理から、
$$ P(M_j \mid \bm{x}) = \frac{P(\bm{x} \mid M_j) P(M_j)}{P(\bm{x})} $$
2つのモデルの事後確率の比（事後オッズ）は、
$$ \frac{P(M_1 \mid \bm{x})}{P(M_2 \mid \bm{x})} = \frac{P(\bm{x} \mid M_1)}{P(\bm{x} \mid M_2)} \cdot \frac{P(M_1)}{P(M_2)} $$
オッズの分解
$$ \underbrace{\frac{P(M_1 \mid \bm{x})}{P(M_2 \mid \bm{x})}}_{\text{事後オッズ}} = \underbrace{BF_{12}}_{\text{ベイズファクター}} \times \underbrace{\frac{P(M_1)}{P(M_2)}}_{\text{事前オッズ}} $$
この関係式は非常に重要です。
- 事前オッズ: データを見る前の2つのモデルへの確信の比
- ベイズファクター: データによる証拠の強さ
- 事後オッズ: データを見た後のモデルへの確信の比
ベイズファクターは、事前の信念をデータによってどれだけ更新すべきかを示す「証拠の重み」です。事前オッズが等しい（$P(M_1) = P(M_2) = 0.5$）場合、事後オッズはベイズファクターそのものになります。
ジェフリーズのスケール
Harold Jeffreysは、ベイズファクターの大きさを解釈するための目安を提案しました。
| $BF_{12}$ | $\log_{10} BF_{12}$ | $M_1$ を支持する証拠の強さ |
|---|---|---|
| 1 〜 3.2 | 0 〜 0.5 | わずかな証拠（Barely worth mentioning） |
| 3.2 〜 10 | 0.5 〜 1 | 実質的な証拠（Substantial） |
| 10 〜 32 | 1 〜 1.5 | 強い証拠（Strong） |
| 32 〜 100 | 1.5 〜 2 | 非常に強い証拠（Very strong） |
| > 100 | > 2 | 決定的な証拠（Decisive） |
$BF_{12} < 1$ の場合は $M_2$ を支持しており、$1/BF_{12}$ をジェフリーズのスケールに当てはめて解釈します。
周辺尤度の計算方法
解析的計算（共役事前分布）
事前分布と尤度が共役関係にある場合、周辺尤度は閉じた形で計算できます。
正規分布・正規事前分布の場合
データ $x_1, \dots, x_n \sim N(\mu, \sigma^2)$（$\sigma^2$ 既知）で、$\mu$ の事前分布を $\mu \sim N(\mu_0, \tau^2)$ とします。
尤度は
$$ P(\bm{x} \mid \mu) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i – \mu)^2}{2\sigma^2}\right) = (2\pi\sigma^2)^{-n/2} \exp\left(-\frac{\sum(x_i – \mu)^2}{2\sigma^2}\right) $$
事前分布は
$$ \pi(\mu) = \frac{1}{\sqrt{2\pi\tau^2}} \exp\left(-\frac{(\mu – \mu_0)^2}{2\tau^2}\right) $$
周辺尤度は
$$ P(\bm{x}) = \int_{-\infty}^{\infty} P(\bm{x} \mid \mu) \pi(\mu) \, d\mu $$
この積分を計算します。指数部分をまとめると、
$$ -\frac{\sum(x_i – \mu)^2}{2\sigma^2} – \frac{(\mu – \mu_0)^2}{2\tau^2} $$
$\sum(x_i – \mu)^2 = \sum(x_i – \bar{x})^2 + n(\bar{x} – \mu)^2$ を代入すると、
$$ -\frac{\sum(x_i – \bar{x})^2}{2\sigma^2} – \frac{n(\bar{x} – \mu)^2}{2\sigma^2} – \frac{(\mu – \mu_0)^2}{2\tau^2} $$
$\mu$ に依存する部分のみ抜き出すと、
$$ -\frac{1}{2}\left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)\mu^2 + \left(\frac{n\bar{x}}{\sigma^2} + \frac{\mu_0}{\tau^2}\right)\mu + \text{const} $$
$A = \frac{n}{\sigma^2} + \frac{1}{\tau^2}$ とおくと、これは $\mu$ に関する正規分布の指数部です。
ガウス積分 $\int_{-\infty}^{\infty} \exp(-\frac{A}{2}\mu^2 + B\mu) \, d\mu = \sqrt{\frac{2\pi}{A}} \exp\left(\frac{B^2}{2A}\right)$ を用いると、
$$ P(\bm{x}) = (2\pi\sigma^2)^{-n/2} \cdot \frac{1}{\sqrt{2\pi\tau^2}} \cdot \sqrt{\frac{2\pi}{A}} \cdot \exp\left(-\frac{\sum(x_i – \bar{x})^2}{2\sigma^2} + \frac{B^2}{2A}\right) $$
ここで $B = \frac{n\bar{x}}{\sigma^2} + \frac{\mu_0}{\tau^2}$ です。
整理すると、$1/A = \frac{\sigma^2\tau^2}{n\tau^2 + \sigma^2}$ なので、
$$ P(\bm{x}) = (2\pi\sigma^2)^{-n/2} \cdot \frac{1}{\sqrt{n\tau^2/\sigma^2 + 1}} \cdot \exp\left(-\frac{\sum(x_i – \bar{x})^2}{2\sigma^2} – \frac{(\bar{x} – \mu_0)^2}{2(\sigma^2/n + \tau^2)}\right) $$
最後の指数部分は、$\bar{x}$ の周辺分布が $\bar{x} \sim N(\mu_0, \sigma^2/n + \tau^2)$ であることを反映しています。
モンテカルロ法
共役事前分布でない場合は、モンテカルロ法で周辺尤度を近似します。
単純モンテカルロ法
事前分布 $\pi(\theta)$ から $M$ 個の標本 $\theta^{(1)}, \dots, \theta^{(M)}$ を生成し、
$$ P(\bm{x} \mid M_j) \approx \frac{1}{M}\sum_{m=1}^{M} P(\bm{x} \mid \theta^{(m)}, M_j) $$
これは $E_{\pi}[P(\bm{x} \mid \theta)]$ のモンテカルロ推定です。
問題点
単純モンテカルロ法は、事前分布と事後分布が大きく異なる場合に効率が悪くなります。事前分布からサンプリングした $\theta$ の多くは尤度が低い領域にあるため、多数の標本が必要になります。
ハーモニック平均法
事後分布からのサンプル $\theta^{(1)}, \dots, \theta^{(M)}$（MCMCの出力）を使って周辺尤度を推定する方法です。
$$ \frac{1}{P(\bm{x})} = E_{\text{post}}\left[\frac{1}{P(\bm{x} \mid \theta)}\right] $$
したがって、
$$ P(\bm{x}) \approx \left[\frac{1}{M}\sum_{m=1}^{M} \frac{1}{P(\bm{x} \mid \theta^{(m)})}\right]^{-1} $$
ただし、ハーモニック平均法は分散が非常に大きくなりやすく、実用上は推奨されません。
ブリッジサンプリング
より精度の高い方法として ブリッジサンプリング があります。事前分布と事後分布の両方からのサンプルを使い、最適なブリッジ関数を介して周辺尤度を推定します。
正規分布の平均検定でのベイズファクター
問題設定
データ $x_1, \dots, x_n \sim N(\mu, \sigma^2)$ に対して、
$$ M_0: \mu = \mu_0 \quad \text{vs} \quad M_1: \mu \neq \mu_0 $$
$\sigma^2$ は既知とします。$M_1$ のもとでは $\mu$ の事前分布を $\mu \sim N(\mu_0, \tau^2)$ とします（$M_0$ のもとでは $\mu = \mu_0$ で固定なのでパラメータなし）。
$M_0$ の周辺尤度
$M_0$ のもとではパラメータが固定されているので、
$$ P(\bm{x} \mid M_0) = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x_i – \mu_0)^2}{2\sigma^2}\right) $$
$$ = (2\pi\sigma^2)^{-n/2}\exp\left(-\frac{\sum(x_i – \mu_0)^2}{2\sigma^2}\right) $$
$M_1$ の周辺尤度
前節で導出した結果を使います。
$$ P(\bm{x} \mid M_1) = (2\pi\sigma^2)^{-n/2} \cdot \frac{1}{\sqrt{n\tau^2/\sigma^2 + 1}} \cdot \exp\left(-\frac{\sum(x_i – \bar{x})^2}{2\sigma^2} – \frac{(\bar{x} – \mu_0)^2}{2(\sigma^2/n + \tau^2)}\right) $$
ベイズファクターの導出
$$ BF_{01} = \frac{P(\bm{x} \mid M_0)}{P(\bm{x} \mid M_1)} $$
分子と分母で共通の $(2\pi\sigma^2)^{-n/2}$ は約分されます。
$$ \begin{align} BF_{01} &= \sqrt{n\tau^2/\sigma^2 + 1} \cdot \exp\left(-\frac{\sum(x_i – \mu_0)^2}{2\sigma^2} + \frac{\sum(x_i – \bar{x})^2}{2\sigma^2} + \frac{(\bar{x} – \mu_0)^2}{2(\sigma^2/n + \tau^2)}\right) \end{align} $$
$\sum(x_i – \mu_0)^2 = \sum(x_i – \bar{x})^2 + n(\bar{x} – \mu_0)^2$ を使うと、
$$ \begin{align} BF_{01} &= \sqrt{1 + n\tau^2/\sigma^2} \cdot \exp\left(-\frac{n(\bar{x} – \mu_0)^2}{2\sigma^2} + \frac{(\bar{x} – \mu_0)^2}{2(\sigma^2/n + \tau^2)}\right) \end{align} $$
指数部分を整理します。$r = n\tau^2/\sigma^2$ とおくと、$\sigma^2/n + \tau^2 = \sigma^2(1 + r)/n$ ですから、
$$ \frac{n(\bar{x} – \mu_0)^2}{2\sigma^2} – \frac{(\bar{x} – \mu_0)^2}{2\sigma^2(1+r)/n} = \frac{n(\bar{x} – \mu_0)^2}{2\sigma^2}\left(1 – \frac{1}{1+r}\right) = \frac{n(\bar{x} – \mu_0)^2}{2\sigma^2} \cdot \frac{r}{1+r} $$
z統計量 $z = \frac{\bar{x} – \mu_0}{\sigma/\sqrt{n}}$ を使うと、$n(\bar{x} – \mu_0)^2/\sigma^2 = z^2$ なので、
$$ \boxed{BF_{01} = \sqrt{1 + r} \cdot \exp\left(-\frac{z^2}{2} \cdot \frac{r}{1 + r}\right)} $$
ここで $r = n\tau^2/\sigma^2$ です。
ベイズファクターの解釈
$BF_{01} > 1$: $M_0$（帰無仮説）を支持する証拠
$BF_{01} < 1$: $M_1$（対立仮説）を支持する証拠
$BF_{10} = 1/BF_{01}$ で、対立仮説を支持するベイズファクターに変換できます。
事前分布の幅 $\tau$ の影響
$r = n\tau^2/\sigma^2$ を通じて、事前分布の幅 $\tau$ がベイズファクターに影響します。
- $\tau$ が大きい（広い事前分布）: $M_1$ は広い範囲の $\mu$ に事前確率を分散させるため、データがない限り $M_0$ に有利になる（リンドレーのパラドックス）
- $\tau$ が小さい（狭い事前分布）: $M_1$ の事前分布が $M_0$ と似てくるため、判別が難しくなる
p値との比較
根本的な違い
| p値 | ベイズファクター | |
|---|---|---|
| 定義 | $P(\text{data以上に極端} \mid H_0)$ | $P(\text{data} \mid M_1) / P(\text{data} \mid M_0)$ |
| 解釈 | $H_0$ の否定の証拠 | 2つのモデルの相対的な支持度 |
| $H_0$ の支持 | 不可（棄却しないだけ） | 可能（$BF_{01} > 1$） |
| 事前分布 | 不要 | 必要 |
| 標本サイズ依存 | $n \to \infty$ で $p \to 0$ | モデルに応じて有限値に収束しうる |
Lindley’s paradox
大標本では、p値とベイズファクターが矛盾する場合があります。
$n$ が非常に大きいとき、わずかな効果でも $p < 0.05$ で統計的に有意になります。しかし、ベイズファクターは対立仮説の事前分布の広がりにペナルティを課すため、$BF_{01} > 1$（帰無仮説を支持）となることがあります。
これを リンドレーのパラドックス と呼びます。p値が「有意」を示しても、ベイズファクターが「帰無仮説を支持」する場合、データは帰無仮説からのわずかな逸脱を示しているだけで、実質的に重要な差ではない可能性があります。
具体例
数値例
工場の製品の重量が設計値 $\mu_0 = 100$ g であるかを検定します。$\sigma = 5$ g（既知）、$n = 25$ 個の測定で $\bar{x} = 101.5$ g でした。
$M_1$ の事前分布: $\mu \sim N(100, 10^2)$（$\tau = 10$）
z統計量:
$$ z = \frac{101.5 – 100}{5/\sqrt{25}} = \frac{1.5}{1.0} = 1.5 $$
パラメータ $r$:
$$ r = \frac{n\tau^2}{\sigma^2} = \frac{25 \times 100}{25} = 100 $$
ベイズファクター:
$$ BF_{01} = \sqrt{1 + 100} \cdot \exp\left(-\frac{1.5^2}{2} \cdot \frac{100}{101}\right) = \sqrt{101} \cdot \exp\left(-\frac{2.25}{2} \cdot 0.990\right) $$
$$ = 10.05 \times \exp(-1.114) = 10.05 \times 0.329 = 3.31 $$
$BF_{01} = 3.31 > 1$ であり、ジェフリーズのスケールでは「実質的な証拠」の境界付近で、帰無仮説を弱く支持しています。
一方、p値は $p = 2(1 – \Phi(1.5)) = 2 \times 0.0668 = 0.134$ であり、有意水準5%では帰無仮説を棄却しません。この場合はp値とベイズファクターの結論は一致しています。
Pythonでの実装
ベイズファクターの計算（正規分布・分散既知）
import numpy as np
from scipy import stats
def bayes_factor_normal_known_var(data, mu_0, sigma, tau):
"""
正規分布の平均の検定に対するベイズファクター（分散既知）
H0: μ = μ_0 vs H1: μ ~ N(μ_0, τ²)
Parameters
----------
data : array-like
標本データ
mu_0 : float
帰無仮説のもとでの平均
sigma : float
母標準偏差（既知）
tau : float
対立仮説での事前分布の標準偏差
Returns
-------
BF01 : float
帰無仮説を支持するベイズファクター
BF10 : float
対立仮説を支持するベイズファクター
"""
data = np.array(data, dtype=float)
n = len(data)
x_bar = np.mean(data)
# z統計量
z = (x_bar - mu_0) / (sigma / np.sqrt(n))
# r パラメータ
r = n * tau**2 / sigma**2
# ベイズファクター BF01
BF01 = np.sqrt(1 + r) * np.exp(-z**2 / 2 * r / (1 + r))
BF10 = 1 / BF01
return BF01, BF10
# --- 具体例 ---
np.random.seed(42)
mu_0 = 100 # 帰無仮説の母平均
sigma = 5.0 # 母標準偏差（既知）
tau = 10.0 # 事前分布の標準偏差
n = 25
# データ生成
data = np.random.normal(101.5, sigma, n)
x_bar = np.mean(data)
z = (x_bar - mu_0) / (sigma / np.sqrt(n))
BF01, BF10 = bayes_factor_normal_known_var(data, mu_0, sigma, tau)
p_value = 2 * (1 - stats.norm.cdf(np.abs(z)))
print("=== ベイズファクター vs p値 ===")
print(f"標本サイズ: n = {n}")
print(f"標本平均: x̄ = {x_bar:.4f}")
print(f"z統計量: z = {z:.4f}")
print(f"p値: {p_value:.4f}")
print()
print(f"BF₀₁ = {BF01:.4f} (H₀を支持)")
print(f"BF₁₀ = {BF10:.4f} (H₁を支持)")
print()
# ジェフリーズのスケール
if BF01 > 1:
bf_for_scale = BF01
direction = "H₀を支持"
else:
bf_for_scale = BF10
direction = "H₁を支持"
if bf_for_scale > 100:
strength = "決定的"
elif bf_for_scale > 30:
strength = "非常に強い"
elif bf_for_scale > 10:
strength = "強い"
elif bf_for_scale > 3:
strength = "実質的"
else:
strength = "わずか"
print(f"ジェフリーズのスケール: {strength}な証拠で{direction}")
事前分布の幅がベイズファクターに与える影響
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# 固定パラメータ
mu_0 = 0
sigma = 1.0
n = 30
# z値を変えてBFを計算
z_values = np.linspace(0, 4, 100)
tau_values = [0.5, 1.0, 2.0, 5.0, 10.0]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
# BF10 vs z
for tau in tau_values:
r = n * tau**2 / sigma**2
BF01 = np.sqrt(1 + r) * np.exp(-z_values**2 / 2 * r / (1 + r))
BF10 = 1 / BF01
ax1.plot(z_values, BF10, linewidth=2, label=f'$\\tau$ = {tau}')
ax1.axhline(y=1, color='gray', linestyle='--', linewidth=1)
ax1.axhline(y=3, color='orange', linestyle=':', linewidth=1, label='BF = 3')
ax1.axhline(y=10, color='red', linestyle=':', linewidth=1, label='BF = 10')
ax1.set_xlabel('|z| statistic', fontsize=12)
ax1.set_ylabel('$BF_{10}$', fontsize=12)
ax1.set_title('Bayes Factor vs z-statistic', fontsize=13)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')
ax1.set_ylim(0.01, 1000)
# p値 vs BF10（固定 τ = 1）
tau = 1.0
r = n * tau**2 / sigma**2
BF01 = np.sqrt(1 + r) * np.exp(-z_values**2 / 2 * r / (1 + r))
BF10 = 1 / BF01
p_values = 2 * (1 - stats.norm.cdf(z_values))
ax2.plot(p_values, BF10, 'b-', linewidth=2)
ax2.axvline(x=0.05, color='red', linestyle='--', linewidth=1.5, label='p = 0.05')
ax2.axhline(y=1, color='gray', linestyle='--', linewidth=1)
ax2.axhline(y=3, color='orange', linestyle=':', linewidth=1, label='BF = 3')
# p = 0.05のときのBFを表示
z_005 = stats.norm.ppf(0.975)
BF01_005 = np.sqrt(1 + r) * np.exp(-z_005**2 / 2 * r / (1 + r))
BF10_005 = 1 / BF01_005
ax2.plot(0.05, BF10_005, 'ro', markersize=10)
ax2.annotate(f'p=0.05: BF₁₀={BF10_005:.2f}', xy=(0.05, BF10_005),
xytext=(0.15, BF10_005*2), fontsize=10,
arrowprops=dict(arrowstyle='->', color='red'))
ax2.set_xlabel('p-value', fontsize=12)
ax2.set_ylabel('$BF_{10}$', fontsize=12)
ax2.set_title(f'p-value vs Bayes Factor ($\\tau$ = {tau}, n = {n})', fontsize=13)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_yscale('log')
ax2.set_xlim(0, 0.5)
plt.tight_layout()
plt.show()
モンテカルロ法による周辺尤度の推定
import numpy as np
from scipy import stats
def marginal_likelihood_mc(data, mu_0, sigma, tau, M=100000, random_state=42):
"""
単純モンテカルロ法による周辺尤度の推定
Parameters
----------
data : array-like
標本データ
mu_0 : float
事前分布の平均
sigma : float
母標準偏差（既知）
tau : float
事前分布の標準偏差
M : int
モンテカルロサンプル数
Returns
-------
ml_mc : float
周辺尤度のモンテカルロ推定値
ml_se : float
推定値の標準誤差
"""
rng = np.random.RandomState(random_state)
data = np.array(data, dtype=float)
n = len(data)
# 事前分布からサンプリング
mu_samples = rng.normal(mu_0, tau, M)
# 各サンプルに対する対数尤度
log_likes = np.zeros(M)
for m in range(M):
log_likes[m] = np.sum(stats.norm.logpdf(data, loc=mu_samples[m], scale=sigma))
# 数値安定性のためにlog-sum-expトリック
max_ll = np.max(log_likes)
ml_mc = np.exp(max_ll) * np.mean(np.exp(log_likes - max_ll))
# 標準誤差
likes = np.exp(log_likes)
ml_se = np.std(likes) / np.sqrt(M)
return ml_mc, ml_se
# --- 具体例 ---
np.random.seed(42)
mu_0 = 100
sigma = 5.0
tau = 10.0
n = 25
data = np.random.normal(101.5, sigma, n)
# M0の周辺尤度（パラメータなし、直接計算）
log_ml_m0 = np.sum(stats.norm.logpdf(data, loc=mu_0, scale=sigma))
ml_m0 = np.exp(log_ml_m0)
# M1の周辺尤度（モンテカルロ推定）
ml_m1_mc, ml_m1_se = marginal_likelihood_mc(data, mu_0, sigma, tau, M=500000)
# M1の周辺尤度（解析解）
x_bar = np.mean(data)
r = n * tau**2 / sigma**2
marginal_var = sigma**2 / n + tau**2
log_ml_m1_exact = (np.sum(stats.norm.logpdf(data, loc=x_bar, scale=sigma))
+ stats.norm.logpdf(x_bar, loc=mu_0, scale=np.sqrt(marginal_var))
- stats.norm.logpdf(x_bar, loc=x_bar, scale=sigma/np.sqrt(n)))
# ベイズファクター
BF01_mc = ml_m0 / ml_m1_mc
z = (x_bar - mu_0) / (sigma / np.sqrt(n))
BF01_exact = np.sqrt(1 + r) * np.exp(-z**2 / 2 * r / (1 + r))
print("=== 周辺尤度の計算方法の比較 ===")
print(f"\nM₀の周辺尤度（解析）: log P(x|M₀) = {log_ml_m0:.4f}")
print(f"M₁の周辺尤度（MC）: P(x|M₁) = {ml_m1_mc:.4e} (SE = {ml_m1_se:.4e})")
print(f"\nBF₀₁（MC法）: {BF01_mc:.4f}")
print(f"BF₀₁（解析解）: {BF01_exact:.4f}")
リンドレーのパラドックスの可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 標本サイズを変えたときのp値とBFの挙動
sample_sizes = np.arange(10, 2010, 10)
sigma = 1.0
tau = 1.0
# 固定された効果量: δ = 0.1（非常に小さな効果）
delta = 0.1
mu_0 = 0
mu_true = delta
p_values = []
BF01_values = []
for n in sample_sizes:
# z統計量（母集団の効果から直接計算）
se = sigma / np.sqrt(n)
z = delta / se
# p値
p = 2 * (1 - stats.norm.cdf(z))
p_values.append(p)
# ベイズファクター
r = n * tau**2 / sigma**2
BF01 = np.sqrt(1 + r) * np.exp(-z**2 / 2 * r / (1 + r))
BF01_values.append(BF01)
p_values = np.array(p_values)
BF01_values = np.array(BF01_values)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# p値の推移
ax1.plot(sample_sizes, p_values, 'b-', linewidth=2)
ax1.axhline(y=0.05, color='red', linestyle='--', linewidth=1.5,
label='$\\alpha = 0.05$')
ax1.set_xlabel('Sample Size $n$', fontsize=12)
ax1.set_ylabel('p-value', fontsize=12)
ax1.set_title(f"Lindley's Paradox: p-value vs n (true $\\delta$ = {delta})", fontsize=14)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_yscale('log')
# n*のマーク（p < 0.05になる点）
n_star = sample_sizes[np.argmax(p_values < 0.05)]
ax1.axvline(x=n_star, color='green', linestyle=':', linewidth=1.5,
label=f'p < 0.05 at n = {n_star}')
ax1.legend(fontsize=11)
# ベイズファクターの推移
ax2.plot(sample_sizes, BF01_values, 'r-', linewidth=2, label='$BF_{01}$')
ax2.axhline(y=1, color='gray', linestyle='--', linewidth=1)
ax2.axhline(y=3, color='orange', linestyle=':', linewidth=1,
label='$BF_{01} = 3$ (substantial for $H_0$)')
ax2.axvline(x=n_star, color='green', linestyle=':', linewidth=1.5,
label=f'p < 0.05 at n = {n_star}')
ax2.set_xlabel('Sample Size $n$', fontsize=12)
ax2.set_ylabel('$BF_{01}$', fontsize=12)
ax2.set_title(f"Lindley's Paradox: Bayes Factor vs n (true $\\delta$ = {delta}, $\\tau$ = {tau})",
fontsize=14)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_yscale('log')
plt.tight_layout()
plt.show()
print(f"p値が0.05を下回るn: {n_star}")
print(f"そのときのBF01: {BF01_values[np.argmax(p_values < 0.05)]:.4f}")
print("→ p値は有意だが、BFは帰無仮説を支持する方向（リンドレーのパラドックス）")
異なる事前分布でのベイズファクターの感度分析
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# データ生成
n = 50
sigma = 1.0
mu_0 = 0
data = np.random.normal(0.3, sigma, n) # 小さな効果
x_bar = np.mean(data)
z = x_bar / (sigma / np.sqrt(n))
p_value = 2 * (1 - stats.norm.cdf(np.abs(z)))
print(f"x̄ = {x_bar:.4f}, z = {z:.4f}, p = {p_value:.4f}")
# 様々なτでベイズファクターを計算
tau_range = np.logspace(-2, 2, 200)
BF10_values = []
for tau in tau_range:
r = n * tau**2 / sigma**2
BF01 = np.sqrt(1 + r) * np.exp(-z**2 / 2 * r / (1 + r))
BF10_values.append(1 / BF01)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(tau_range, BF10_values, 'b-', linewidth=2.5)
ax.axhline(y=1, color='gray', linestyle='--', linewidth=1)
ax.axhline(y=3, color='orange', linestyle=':', linewidth=1, label='BF = 3')
ax.axhline(y=10, color='red', linestyle=':', linewidth=1, label='BF = 10')
ax.axhline(y=1/3, color='green', linestyle=':', linewidth=1, label='BF = 1/3')
# JZSベイズファクターに使われるτ=1の位置
ax.axvline(x=1, color='purple', linestyle='--', linewidth=1.5,
label='$\\tau = 1$ (JZS default)')
ax.set_xlabel('Prior Width $\\tau$', fontsize=13)
ax.set_ylabel('$BF_{10}$', fontsize=13)
ax.set_title(f'Sensitivity Analysis: BF vs Prior Width (z = {z:.2f}, p = {p_value:.4f})',
fontsize=14)
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0.01, 100)
plt.tight_layout()
plt.show()
まとめ
本記事では、ベイズファクターとベイズ検定の理論について解説しました。
- ベイズファクター: 2つのモデルの周辺尤度の比 $BF_{12} = P(\bm{x} \mid M_1) / P(\bm{x} \mid M_2)$
- 事後オッズ: 事後オッズ = ベイズファクター $\times$ 事前オッズ
- ジェフリーズのスケール: $BF > 3$（実質的）、$> 10$（強い）、$> 100$（決定的）
- 周辺尤度: 尤度を事前分布で積分。共役事前分布なら解析的に計算可能、そうでなければモンテカルロ法
- 正規分布の検定: $BF_{01} = \sqrt{1+r} \cdot \exp(-z^2 r / [2(1+r)])$（$r = n\tau^2/\sigma^2$）
- p値との違い: ベイズファクターは帰無仮説の支持も定量化でき、リンドレーのパラドックスでp値と矛盾することがある
- 事前分布の影響: 事前分布の幅 $\tau$ が結果に影響するため、感度分析が重要
次のステップとして、以下の記事も参考にしてください。