F分布の導出と分散分析(ANOVA)への応用

「3つの薬の効果に差があるか」「5つの工場の品質に違いがあるか」— このような「複数グループの平均の差の検定」は、t検定を繰り返すのではなく分散分析(ANOVA)で行います。ANOVAの検定統計量が従う分布がF分布(F-distribution)です。

F分布はロナルド・フィッシャー(Ronald Fisher)にちなんで名付けられ、統計的推測の中核をなす分布です。

  • 分散分析(ANOVA): 複数グループの平均の比較
  • 等分散検定: 2つの母集団の分散の比の検定
  • 回帰分析: モデル全体の有意性のF検定
  • 実験計画法: 要因効果の検定

本記事の内容

  • F分布の定義(カイ二乗分布の比)
  • 密度関数と期待値・分散
  • t分布との関係
  • 分散分析(一元配置ANOVA)への応用
  • Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

F分布の定義

カイ二乗分布の比としての導出

$U \sim \chi^2(d_1)$ と $V \sim \chi^2(d_2)$ が独立のとき

$$ \begin{equation} F = \frac{U/d_1}{V/d_2} \sim F(d_1, d_2) \end{equation} $$

$d_1$ は分子の自由度、$d_2$ は分母の自由度です。

直感的には、F分布は「2つの分散の比」の分布です。分子も分母もカイ二乗分布(正規分布の二乗和)をその自由度で割ったものであり、もし2つの母分散が等しければ $F \approx 1$ になるはずです。

密度関数

$$ f(x) = \frac{1}{B(d_1/2, d_2/2)} \left(\frac{d_1}{d_2}\right)^{d_1/2} \frac{x^{d_1/2 – 1}}{(1 + d_1 x/d_2)^{(d_1+d_2)/2}}, \quad x > 0 $$

期待値と分散

$$ E[F] = \frac{d_2}{d_2 – 2} \quad (d_2 > 2) $$

$$ \text{Var}(F) = \frac{2d_2^2(d_1 + d_2 – 2)}{d_1(d_2-2)^2(d_2-4)} \quad (d_2 > 4) $$

$d_2 > 2$ のとき $E[F] > 1$ であり、$d_2 \to \infty$ で $E[F] \to 1$ です。

t分布との関係

$T \sim t(\nu)$ のとき $T^2 \sim F(1, \nu)$ です。つまり、2つのグループの平均の比較(t検定)はF検定($d_1 = 1$ のANOVA)と等価です。

分散分析(一元配置ANOVA)

問題設定

$k$ 個のグループがあり、グループ $i$ から $n_i$ 個のデータ $X_{i1}, \ldots, X_{in_i}$ が得られています。

$$ X_{ij} = \mu_i + \varepsilon_{ij}, \quad \varepsilon_{ij} \sim N(0, \sigma^2) $$

帰無仮説 $H_0: \mu_1 = \mu_2 = \cdots = \mu_k$ を検定します。

F統計量

グループ間変動(SSB)とグループ内変動(SSW)を計算します。

$$ \text{SSB} = \sum_{i=1}^{k} n_i (\bar{X}_{i\cdot} – \bar{X}_{\cdot\cdot})^2, \quad \text{SSW} = \sum_{i=1}^{k} \sum_{j=1}^{n_i} (X_{ij} – \bar{X}_{i\cdot})^2 $$

F統計量は

$$ F = \frac{\text{SSB}/(k-1)}{\text{SSW}/(N-k)} = \frac{\text{MSB}}{\text{MSW}} $$

$H_0$ のもとで $F \sim F(k-1, N-k)$ に従います。

Pythonでの実装と可視化

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

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

# (a) 異なる自由度でのPDF
ax = axes[0, 0]
x = np.linspace(0.01, 5, 500)
params_f = [(1, 1), (2, 5), (5, 5), (10, 10), (10, 50), (50, 50)]

for d1, d2 in params_f:
    pdf = stats.f.pdf(x, d1, d2)
    ax.plot(x, pdf, linewidth=2, label=f'F({d1},{d2})')

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('f(x)', fontsize=12)
ax.set_title('F-Distribution PDF', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.2)

# (b) T^2 = F(1,ν) の確認
ax = axes[0, 1]
np.random.seed(42)
nu = 10
n_sim = 100000
t_samples = np.random.standard_t(nu, n_sim)
f_from_t = t_samples**2

x_range = np.linspace(0, 8, 300)
ax.hist(f_from_t, bins=100, density=True, range=(0, 8),
        color='steelblue', alpha=0.6, label='$T^2$ (simulation)')
ax.plot(x_range, stats.f.pdf(x_range, 1, nu), 'r-', linewidth=2.5,
        label=f'F(1, {nu}) theory')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'$T^2 \\sim F(1, \\nu)$ verification ($\\nu$={nu})', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# (c) ANOVAのシミュレーション
ax = axes[1, 0]
np.random.seed(42)
n_anova_sim = 10000

# H0のもとでのF統計量
f_null = []
k = 3
n_per_group = 20
for _ in range(n_anova_sim):
    groups = [np.random.normal(0, 1, n_per_group) for _ in range(k)]
    all_data = np.concatenate(groups)
    grand_mean = np.mean(all_data)
    ssb = sum(n_per_group * (np.mean(g) - grand_mean)**2 for g in groups)
    ssw = sum(np.sum((g - np.mean(g))**2) for g in groups)
    f_stat = (ssb/(k-1)) / (ssw/(k*n_per_group-k))
    f_null.append(f_stat)

x_range = np.linspace(0, 6, 300)
ax.hist(f_null, bins=80, density=True, color='steelblue', alpha=0.6,
        label='Simulation (H0 true)')
ax.plot(x_range, stats.f.pdf(x_range, k-1, k*n_per_group-k), 'r-',
        linewidth=2.5, label=f'F({k-1}, {k*n_per_group-k}) theory')

# 棄却域
f_crit = stats.f.ppf(0.95, k-1, k*n_per_group-k)
ax.axvline(f_crit, color='green', linestyle='--', linewidth=2,
           label=f'Critical value (α=0.05): {f_crit:.2f}')
ax.fill_between(x_range[x_range >= f_crit],
                stats.f.pdf(x_range[x_range >= f_crit], k-1, k*n_per_group-k),
                alpha=0.3, color='red', label='Rejection region')

ax.set_xlabel('F statistic', fontsize=12)
ax.set_ylabel('Density', fontsize=12)
ax.set_title(f'One-way ANOVA (k={k}, n={n_per_group})', fontsize=12)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# (d) 検出力の比較
ax = axes[1, 1]
effect_sizes = np.linspace(0, 2, 50)
powers = []

for delta in effect_sizes:
    rejections = 0
    for _ in range(2000):
        groups = [np.random.normal(i*delta, 1, n_per_group) for i in range(k)]
        all_data = np.concatenate(groups)
        grand_mean = np.mean(all_data)
        ssb = sum(n_per_group * (np.mean(g) - grand_mean)**2 for g in groups)
        ssw = sum(np.sum((g - np.mean(g))**2) for g in groups)
        f_stat = (ssb/(k-1)) / (ssw/(k*n_per_group-k))
        if f_stat > f_crit:
            rejections += 1
    powers.append(rejections / 2000)

ax.plot(effect_sizes, powers, 'b-', linewidth=2.5)
ax.axhline(0.05, color='red', linestyle='--', linewidth=1.5,
           label='$\\alpha$ = 0.05')
ax.axhline(0.80, color='green', linestyle=':', linewidth=1.5,
           label='Power = 0.80')
ax.set_xlabel('Effect size ($\\delta$)', fontsize=12)
ax.set_ylabel('Power', fontsize=12)
ax.set_title('ANOVA Power Curve', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('f_distribution.png', dpi=150, bbox_inches='tight')
plt.show()

この可視化から、F分布とANOVAの関係が確認できます。

  1. 左上(PDF): 自由度が大きくなるにつれてピークが1に近づき、分布が集中します。$F(1,1)$ は右に長い裾を持ちます。

  2. 右上($T^2 = F(1,\nu)$): t分布の二乗がF分布に一致することが数値的に確認されています。

  3. 左下(ANOVA): $H_0$ のもとでのF統計量のヒストグラムが理論的なF分布(赤い実線)と一致しています。臨界値(緑の破線)を超えた領域が棄却域です。

  4. 右下(検出力): 効果量が大きくなるにつれて検出力が上昇します。一般的な目安である検出力80%を達成するのに必要な効果量が読み取れます。

まとめ

本記事では、F分布の導出から分散分析への応用まで解説しました。

  • F分布は2つの独立なカイ二乗分布の比(各自由度で割ったもの)として定義される
  • t分布との関係: $T^2 \sim F(1, \nu)$
  • 分散分析(ANOVA)のF統計量は「グループ間変動/グループ内変動」の比であり、帰無仮説のもとでF分布に従う
  • 自由度が大きくなると分布は1の周りに集中する

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