2つのグループのデータがあるとき、「グループ間のばらつきは、グループ内のばらつきに比べて本当に大きいのだろうか？」という問いに答えたいことがあります。たとえば、3種類の肥料を使って育てたトマトの収量に有意な差があるかどうか、あるいは2つの測定装置の精度が同じかどうか。これらの問いに答えるには、分散の比がどのような確率分布に従うかを知る必要があります。
この「分散の比の分布」こそがF分布（F-distribution）です。F分布は、統計学者 R. A. フィッシャー（Fisher）の名にちなんで、ジョージ・スネデカー（Snedecor）が命名しました。
F分布を理解すると、以下のような応用が開けます。
- 分散分析（ANOVA）: 3つ以上の母集団の平均値に差があるかを検定する際の中核的な分布
- 等分散性の検定（F検定）: 2つの母集団の分散が等しいかどうかを検定する
- 回帰分析: 回帰モデル全体の有意性を検定するF検定統計量がF分布に従う
- 実験計画法: 要因効果の有意性判定に不可欠
本記事では、F分布が「2つの独立なχ²分布の比」として自然に現れることを直感的に理解し、確率密度関数を厳密に導出します。さらに、モーメントや他の分布との関係を整理し、分散分析への応用をPythonで実装します。
本記事の内容
- F分布の直感的理解と数学的定義
- 確率密度関数の導出
- モーメント（期待値・分散）の計算
- F分布の性質（逆数の関係、t分布との関係、正規近似）
- Pythonによる可視化とモンテカルロシミュレーション
- 分散分析（一元配置ANOVA）への応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
F分布とは — 直感的な理解
分散の比としてのイメージ
F分布を直感的に理解するために、具体的な場面を考えましょう。
ある学校で、A組（30人）とB組（25人）が同じ数学のテストを受けました。両方のクラスが同じ母集団から来ているとすれば、A組の標本分散 $s_A^2$ とB組の標本分散 $s_B^2$ はだいたい同じ値になるはずです。しかし、たまたまの揺らぎで $s_A^2 / s_B^2$ が2になったり0.5になったりすることもあり得ます。
問題は、この比がどの程度大きければ「たまたまではなく、本当に分散が異なる」と判断できるかです。この判断基準を与えるのがF分布です。
$s_A^2 / s_B^2$ の比がどのような分布に従うかを追いかけてみましょう。χ²分布の記事で学んだように、正規母集団の標本分散は
$$ \frac{(n_A – 1)s_A^2}{\sigma^2} \sim \chi^2(n_A – 1), \quad \frac{(n_B – 1)s_B^2}{\sigma^2} \sim \chi^2(n_B – 1) $$
に従います。帰無仮説 $\sigma_A^2 = \sigma_B^2 = \sigma^2$ のもとでは、
$$ F = \frac{s_A^2}{s_B^2} = \frac{\chi^2(n_A – 1)/(n_A – 1)}{\chi^2(n_B – 1)/(n_B – 1)} $$
となります。つまり、F統計量は2つの独立なχ²確率変数をそれぞれの自由度で割った比なのです。
自由度の2つのパラメータ
χ²分布のパラメータが自由度1つだけだったのに対して、F分布は2つの自由度 $d_1$（分子の自由度）と $d_2$（分母の自由度）をパラメータとして持ちます。これは、分子と分母にそれぞれ異なるχ²分布が入るためです。
直感的には、$d_1$ と $d_2$ が両方とも大きくなるほど、$F$ の値は1の周りに集中していきます。これは、サンプルサイズが大きければ標本分散が母分散に近づくため、比が1に近くなるからです。逆に、自由度が小さいと分散の推定が不安定になり、$F$ の分布の裾が重くなります。
この直感を数学的に定式化しましょう。
F分布の数学的定義
$Q_1 \sim \chi^2(d_1)$ と $Q_2 \sim \chi^2(d_2)$ が独立であるとき、
$$ \begin{equation} F = \frac{Q_1 / d_1}{Q_2 / d_2} \end{equation} $$
が従う分布を自由度 $(d_1, d_2)$ のF分布といい、$F \sim F(d_1, d_2)$ と書きます。
この分布の確率密度関数は次の式で与えられます。
$$ \begin{equation} f(x;\, d_1, d_2) = \frac{1}{B\!\left(\frac{d_1}{2},\, \frac{d_2}{2}\right)} \left(\frac{d_1}{d_2}\right)^{d_1/2} \frac{x^{d_1/2 – 1}}{\left(1 + \frac{d_1}{d_2}x\right)^{(d_1 + d_2)/2}}, \quad x > 0 \end{equation} $$
ここで $B(\cdot, \cdot)$ はベータ関数です。
この式は一見複雑ですが、構造を読み解くことができます。$x^{d_1/2 – 1}$ が $x = 0$ 付近の振る舞いを制御し、$(1 + d_1 x/d_2)^{-(d_1+d_2)/2}$ が右裾の減衰速度を制御しています。べき乗の減衰であるため、χ²分布の指数的減衰よりもゆっくり減衰し、裾が重い分布になっています。
次に、この確率密度関数の導出を見ていきましょう。
確率密度関数の導出
導出の方針
$F = \frac{Q_1/d_1}{Q_2/d_2} = \frac{d_2 Q_1}{d_1 Q_2}$ の確率密度関数を求めます。$Q_1$ と $Q_2$ は独立なχ²確率変数です。
導出は次の3ステップで進めます。
- $Q_1$ と $Q_2$ の同時確率密度関数を書き下す
- 変数変換 $(Q_1, Q_2) \to (F, Q_2)$ を行う
- $Q_2$ について積分（周辺化）して $F$ の確率密度関数を得る
ステップ1: 同時確率密度関数
$Q_1 \sim \chi^2(d_1)$, $Q_2 \sim \chi^2(d_2)$ が独立なので、同時確率密度関数はそれぞれの確率密度関数の積です。
$$ f_{Q_1, Q_2}(q_1, q_2) = \frac{q_1^{d_1/2-1}\,e^{-q_1/2}}{2^{d_1/2}\,\Gamma(d_1/2)} \cdot \frac{q_2^{d_2/2-1}\,e^{-q_2/2}}{2^{d_2/2}\,\Gamma(d_2/2)} $$
ステップ2: 変数変換
$F = \frac{d_2 q_1}{d_1 q_2}$ と定義します。$q_2$ はそのまま残し、$q_1 = \frac{d_1}{d_2} F \cdot q_2$ と表します。
ヤコビアンを計算します。
$$ \frac{\partial q_1}{\partial F} = \frac{d_1}{d_2}\,q_2, \quad \frac{\partial q_1}{\partial q_2} = \frac{d_1}{d_2}\,F $$
$$ \frac{\partial q_2}{\partial F} = 0, \quad \frac{\partial q_2}{\partial q_2} = 1 $$
したがって、ヤコビアンの絶対値は
$$ |J| = \left|\frac{d_1}{d_2}\,q_2 \cdot 1 – \frac{d_1}{d_2}\,F \cdot 0\right| = \frac{d_1}{d_2}\,q_2 $$
変換後の同時確率密度関数は $f_{Q_1, Q_2}(q_1(F, q_2),\, q_2) \cdot |J|$ です。$q_1 = \frac{d_1}{d_2}F q_2$ を代入すると、
$$ f_{F, Q_2}(F, q_2) = \frac{1}{2^{(d_1+d_2)/2}\,\Gamma(d_1/2)\,\Gamma(d_2/2)} \left(\frac{d_1}{d_2}F q_2\right)^{d_1/2-1} q_2^{d_2/2-1} \exp\left(-\frac{q_2}{2}\left(\frac{d_1}{d_2}F + 1\right)\right) \cdot \frac{d_1}{d_2}q_2 $$
$q_2$ のべきを整理します。$q_2$ の因子は $q_2^{d_1/2-1} \cdot q_2^{d_2/2-1} \cdot q_2 = q_2^{(d_1+d_2)/2 – 1}$ です。
ステップ3: 周辺化
$F$ の確率密度関数を得るために $q_2$ について $0$ から $\infty$ まで積分します。
$$ f_F(F) = \int_0^\infty f_{F, Q_2}(F, q_2)\,dq_2 $$
$F$ に依存しない因子をまとめ、$q_2$ の積分部分を抽出します。$u = \frac{q_2}{2}\left(\frac{d_1}{d_2}F + 1\right)$ と置換すると、この積分はガンマ関数の積分表現に帰着します。
$$ \int_0^\infty q_2^{(d_1+d_2)/2 – 1} \exp\left(-\frac{q_2}{2}\left(\frac{d_1 F}{d_2} + 1\right)\right) dq_2 = \frac{\Gamma\left(\frac{d_1+d_2}{2}\right)}{\left(\frac{1}{2}\left(\frac{d_1 F}{d_2} + 1\right)\right)^{(d_1+d_2)/2}} $$
すべての因子を集めて整理すると、次の確率密度関数が得られます。
$$ f_F(x) = \frac{\Gamma\left(\frac{d_1+d_2}{2}\right)}{\Gamma\left(\frac{d_1}{2}\right)\Gamma\left(\frac{d_2}{2}\right)} \left(\frac{d_1}{d_2}\right)^{d_1/2} \frac{x^{d_1/2 – 1}}{\left(1 + \frac{d_1}{d_2}x\right)^{(d_1 + d_2)/2}} $$
ベータ関数 $B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$ を使うと、先ほどの定義式と一致することが確認できます。
$$ \begin{equation} f(x;\, d_1, d_2) = \frac{1}{B\!\left(\frac{d_1}{2},\, \frac{d_2}{2}\right)} \left(\frac{d_1}{d_2}\right)^{d_1/2} \frac{x^{d_1/2 – 1}}{\left(1 + \frac{d_1}{d_2}x\right)^{(d_1 + d_2)/2}} \end{equation} $$
導出は完了です。この確率密度関数がχ²分布の比の定義から自然に導かれることを確認しました。
次に、F分布の基本的な統計量を計算しましょう。
期待値と分散
期待値
$F \sim F(d_1, d_2)$ の期待値は、$d_2 > 2$ のとき存在し、
$$ \begin{equation} E[F] = \frac{d_2}{d_2 – 2} \quad (d_2 > 2) \end{equation} $$
で与えられます。
導出: $F = \frac{Q_1/d_1}{Q_2/d_2}$ であり、$Q_1$ と $Q_2$ が独立なので、
$$ E[F] = \frac{d_2}{d_1} E[Q_1] \cdot E\left[\frac{1}{Q_2}\right] $$
$E[Q_1] = d_1$ はχ²分布の期待値です。$Q_2 \sim \chi^2(d_2)$（= $\text{Gamma}(d_2/2, 2)$）の逆数の期待値は、
$$ E\left[\frac{1}{Q_2}\right] = \frac{1}{2^{d_2/2}\Gamma(d_2/2)}\int_0^\infty q_2^{d_2/2-2}\,e^{-q_2/2}\,dq_2 = \frac{1}{d_2 – 2} \quad (d_2 > 2) $$
と計算されます（ガンマ関数の性質 $\Gamma(n) = (n-1)\Gamma(n-1)$ を利用）。したがって、
$$ E[F] = \frac{d_2}{d_1} \cdot d_1 \cdot \frac{1}{d_2 – 2} = \frac{d_2}{d_2 – 2} $$
注目すべき点は、期待値が分子の自由度 $d_1$ に依存しないことです。分母の自由度 $d_2$ だけで決まります。$d_2$ が大きくなると $E[F] \to 1$ に近づき、$d_2$ が小さいと $E[F] > 1$ で、期待値が1より大きくなります。
分散
$d_2 > 4$ のとき、分散は次のように求まります。
$$ \begin{equation} \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) \end{equation} $$
$d_2 \leq 4$ のとき分散は存在しない（発散する）ことに注意が必要です。F分布は裾が重い分布であり、分母の自由度が小さいと高次モーメントが発散します。
最頻値
$d_1 > 2$ のとき、最頻値は
$$ \begin{equation} \text{mode} = \frac{d_1 – 2}{d_1} \cdot \frac{d_2}{d_2 + 2} \end{equation} $$
で与えられます。最頻値は常に期待値 $d_2/(d_2-2)$ よりも小さく、分布が右に裾を引いている（正の歪度を持つ）ことが反映されています。
次に、F分布が持つ重要な理論的性質を見ていきましょう。
F分布の重要な性質
逆数の性質
F分布の最も重要な性質の一つは、逆数に関する次の関係です。
$$ \begin{equation} F \sim F(d_1, d_2) \quad \Longrightarrow \quad \frac{1}{F} \sim F(d_2, d_1) \end{equation} $$
すなわち、F分布に従う確率変数の逆数は、自由度の順序を入れ替えたF分布に従います。
証明: $F = \frac{Q_1/d_1}{Q_2/d_2}$ ならば $1/F = \frac{Q_2/d_2}{Q_1/d_1}$ であり、これは定義より $F(d_2, d_1)$ に従います。$\square$
この性質により、F分布のパーセント点の表は片側（上側）だけあれば十分です。下側のパーセント点は逆数の関係から求められます。
$$ F_{\alpha}(d_1, d_2) = \frac{1}{F_{1-\alpha}(d_2, d_1)} $$
t分布との関係
$T \sim t(d)$ のとき、
$$ \begin{equation} T^2 \sim F(1, d) \end{equation} $$
が成り立ちます。これは、t検定（両側）とF検定（$d_1 = 1$ の場合）が同値であることを意味します。
証明: $T = Z/\sqrt{Q/d}$（$Z \sim N(0,1)$, $Q \sim \chi^2(d)$ が独立）のとき、
$$ T^2 = \frac{Z^2}{Q/d} = \frac{Z^2/1}{Q/d} = \frac{\chi^2(1)/1}{\chi^2(d)/d} \sim F(1, d) $$
ベータ分布との関係
$F \sim F(d_1, d_2)$ のとき、次の変換
$$ W = \frac{d_1 F / d_2}{1 + d_1 F / d_2} = \frac{d_1 F}{d_2 + d_1 F} $$
により得られる $W$ は、$\text{Beta}(d_1/2, d_2/2)$ に従います。この変換は $F$ を $[0, \infty)$ から $[0, 1)$ に写し、確率密度関数の導出でも利用されます。
大標本での正規近似
$d_1, d_2$ がともに大きいとき、F分布は正規分布に近づきます。特に、$d_2 \to \infty$ のとき、
$$ d_1 F \xrightarrow{d} \chi^2(d_1) $$
が成り立ちます。これは、分母が $Q_2/d_2 \to E[Q_2]/d_2 = 1$（大数の法則）に確率収束するためです。
これらの性質を可視化で確認してみましょう。
Pythonによる確率密度関数の可視化
さまざまな自由度の組み合わせに対するF分布の確率密度関数を描画します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) d1を固定し、d2を変化
ax = axes[0]
x = np.linspace(0.01, 5, 500)
d1 = 5
for d2 in [5, 10, 20, 50, 100]:
pdf = stats.f.pdf(x, dfn=d1, dfd=d2)
ax.plot(x, pdf, linewidth=2, label=rf"$F({d1}, {d2})$")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("f(x)", fontsize=12)
ax.set_title(rf"Varying $d_2$ with fixed $d_1 = {d1}$", fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1.0)
ax.grid(True, alpha=0.3)
# (b) さまざまな(d1, d2)の組み合わせ
ax = axes[1]
params = [(1, 1), (2, 5), (5, 5), (10, 10), (50, 50), (100, 100)]
for d1, d2 in params:
pdf = stats.f.pdf(x, dfn=d1, dfd=d2)
ax.plot(x, pdf, linewidth=2, label=rf"$F({d1}, {d2})$")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("f(x)", fontsize=12)
ax.set_title("Various F distributions", fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 2.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("f_distribution_pdf.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、F分布の形状について以下の特徴が読み取れます。
-
左図: $d_2$ が大きくなるにつれて分布が1の周辺に集中していく 。$d_2 = 5$ のときは右裾が長く広がっていますが、$d_2 = 100$ では鋭いピークを持つ分布になっています。これは、分母の自由度が大きいほどχ²/自由度の推定精度が上がり、F比のばらつきが小さくなることを反映しています
-
右図: $d_1 = d_2 = 1$ の場合は非常に裾の重い分布 であり、大きな値が出やすいことがわかります。両方の自由度が増加すると、分布は対称的な釣り鐘型に近づいていきます
-
すべてのF分布は正の値のみをとる 。分散の比は常に正であるため、これは自然な性質です
-
分布のピーク（最頻値）は常に1未満 。$d_1 \geq 3$ のとき最頻値は $(d_1-2)/d_1 \cdot d_2/(d_2+2) < 1$ であり、分布は1を中心に左に偏った形をしています
次に、モンテカルロシミュレーションで理論を検証しましょう。
モンテカルロシミュレーションによる検証
χ²分布の比を実際に生成し、F分布の理論と一致するかを確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
n_samples = 200000
param_sets = [(5, 10), (10, 10), (2, 30), (30, 30)]
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
for idx, (d1, d2) in enumerate(param_sets):
ax = axes[idx // 2, idx % 2]
# χ²分布からサンプリング
Q1 = np.random.chisquare(df=d1, size=n_samples)
Q2 = np.random.chisquare(df=d2, size=n_samples)
# F統計量の計算
F_samples = (Q1 / d1) / (Q2 / d2)
# ヒストグラム
upper = np.percentile(F_samples, 99)
ax.hist(F_samples, bins=200, density=True, alpha=0.6,
color="steelblue", edgecolor="white", range=(0, upper),
label="Monte Carlo")
# 理論的な確率密度関数
x = np.linspace(0.01, upper, 500)
pdf = stats.f.pdf(x, dfn=d1, dfd=d2)
ax.plot(x, pdf, "r-", linewidth=2.5, label=rf"$F({d1}, {d2})$ PDF")
# 統計量の比較
mean_theory = d2 / (d2 - 2) if d2 > 2 else float('inf')
ax.set_title(rf"$F({d1}, {d2})$: mean = {F_samples.mean():.3f} "
rf"(theory: {mean_theory:.3f})", fontsize=11)
ax.set_xlabel("x", fontsize=11)
ax.set_ylabel("Density", fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("f_distribution_monte_carlo.png", dpi=150, bbox_inches="tight")
plt.show()
このシミュレーション結果から、以下のことが確認できます。
-
すべてのパラメータの組み合わせで、シミュレーションのヒストグラムが理論曲線と一致している 。200,000サンプルの経験分布は、理論的なF分布の確率密度関数を忠実に再現しています
-
経験的な平均が理論値 $d_2/(d_2-2)$ とよく一致している 。特に $d_2$ が大きい場合、平均は1に近く、$d_2 = 10$ のときは $10/8 = 1.25$ となります
-
$F(2, 30)$ のように $d_1$ が小さい場合、分布は原点近くでのピークが高く、右裾が長い 。これは、分子の自由度が小さいときにχ²/自由度の推定が不安定であることを反映しています
理論と数値の一致を確認しました。次に、F分布の最も重要な応用である分散分析を実装してみましょう。
分散分析（一元配置ANOVA）への応用
ANOVAの理論
分散分析（Analysis of Variance, ANOVA）は、3つ以上の母集団の平均値に差があるかどうかを検定する手法です。一元配置ANOVAでは、$a$ 個のグループがあり、各グループ $i$ に $n_i$ 個の観測値があるとします。
帰無仮説は「すべてのグループの母平均が等しい」、すなわち $H_0: \mu_1 = \mu_2 = \cdots = \mu_a$ です。
ANOVAの基本的なアイデアは、データ全体のばらつき（全変動）を群間変動（between-group variation）と群内変動（within-group variation）に分解することです。
全変動 $SS_T$ は次のように分解されます。
$$ \underbrace{\sum_{i=1}^a \sum_{j=1}^{n_i} (X_{ij} – \bar{X})^2}_{SS_T} = \underbrace{\sum_{i=1}^a n_i(\bar{X}_i – \bar{X})^2}_{SS_B} + \underbrace{\sum_{i=1}^a \sum_{j=1}^{n_i} (X_{ij} – \bar{X}_i)^2}_{SS_W} $$
ここで $\bar{X}$ は全体平均、$\bar{X}_i$ はグループ $i$ の平均です。
$SS_B$（群間変動）は「グループ平均のばらつき」を測り、$SS_W$（群内変動）は「各グループ内部のばらつき」を測ります。
帰無仮説のもとで、以下が成り立ちます。
$$ \frac{SS_B}{\sigma^2} \sim \chi^2(a – 1), \quad \frac{SS_W}{\sigma^2} \sim \chi^2(N – a) $$
ここで $N = \sum_{i=1}^a n_i$ は全サンプル数です。したがって、F統計量
$$ \begin{equation} F = \frac{SS_B/(a-1)}{SS_W/(N-a)} = \frac{MS_B}{MS_W} \end{equation} $$
は $H_0$ のもとで $F(a-1, N-a)$ に従います。$MS_B$（群間平均平方）と $MS_W$（群内平均平方）はそれぞれの変動を自由度で割った値です。
帰無仮説が偽であれば（つまり、グループ間に本当に差があれば）、$SS_B$ は大きくなり、$F$ 統計量が大きくなります。F統計量がF分布の臨界値を超えれば、帰無仮説を棄却します。
PythonによるANOVAの実装
3種類の肥料によるトマトの収量を比較する例をスクラッチで実装します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# 3種類の肥料によるトマトの収量データ（シミュレーション）
# 帰無仮説が偽: グループ間に差がある
n_per_group = 20
group_A = np.random.normal(loc=50, scale=8, size=n_per_group) # 肥料A
group_B = np.random.normal(loc=55, scale=8, size=n_per_group) # 肥料B
group_C = np.random.normal(loc=58, scale=8, size=n_per_group) # 肥料C
groups = [group_A, group_B, group_C]
group_labels = ["Fertilizer A", "Fertilizer B", "Fertilizer C"]
a = len(groups)
N = sum(len(g) for g in groups)
# 全体平均
all_data = np.concatenate(groups)
grand_mean = np.mean(all_data)
# 群間変動 SS_B と群内変動 SS_W の計算
SS_B = sum(len(g) * (np.mean(g) - grand_mean)**2 for g in groups)
SS_W = sum(np.sum((g - np.mean(g))**2) for g in groups)
SS_T = np.sum((all_data - grand_mean)**2)
# 自由度
df_between = a - 1
df_within = N - a
# 平均平方
MS_B = SS_B / df_between
MS_W = SS_W / df_within
# F統計量
F_stat = MS_B / MS_W
# p値
p_value = 1 - stats.f.cdf(F_stat, dfn=df_between, dfd=df_within)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# 箱ひげ図
ax = axes[0]
bp = ax.boxplot(groups, labels=group_labels, patch_artist=True,
boxprops=dict(facecolor="lightblue", alpha=0.7))
colors = ["#5B9BD5", "#ED7D31", "#70AD47"]
for patch, color in zip(bp["boxes"], colors):
patch.set_facecolor(color)
patch.set_alpha(0.6)
# 各グループの平均を表示
for i, g in enumerate(groups):
ax.plot(i + 1, np.mean(g), "D", color="red", markersize=8, zorder=5)
ax.axhline(grand_mean, color="gray", linestyle="--", linewidth=1,
label=f"Grand mean = {grand_mean:.1f}")
ax.set_ylabel("Yield", fontsize=12)
ax.set_title("Tomato yield by fertilizer type", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")
# F分布と検定統計量
ax = axes[1]
x = np.linspace(0, 10, 500)
pdf = stats.f.pdf(x, dfn=df_between, dfd=df_within)
ax.plot(x, pdf, "b-", linewidth=2, label=rf"$F({df_between}, {df_within})$")
ax.fill_between(x, pdf, where=(x >= F_stat), alpha=0.3, color="red",
label=f"p-value = {p_value:.4f}")
ax.axvline(F_stat, color="red", linewidth=2, linestyle="--",
label=rf"$F$ = {F_stat:.2f}")
critical_value = stats.f.ppf(0.95, dfn=df_between, dfd=df_within)
ax.axvline(critical_value, color="green", linewidth=1.5, linestyle=":",
label=f"Critical value (5%) = {critical_value:.2f}")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("F-test for one-way ANOVA", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("f_distribution_anova.png", dpi=150, bbox_inches="tight")
plt.show()
# ANOVA表
print("=" * 60)
print("ANOVA Table")
print("=" * 60)
print(f"{'Source':<12} {'SS':>10} {'df':>5} {'MS':>10} {'F':>8} {'p-value':>10}")
print("-" * 60)
print(f"{'Between':<12} {SS_B:>10.2f} {df_between:>5d} {MS_B:>10.2f} {F_stat:>8.3f} {p_value:>10.6f}")
print(f"{'Within':<12} {SS_W:>10.2f} {df_within:>5d} {MS_W:>10.2f}")
print(f"{'Total':<12} {SS_T:>10.2f} {N-1:>5d}")
print("=" * 60)
print(f"\n判定（有意水準5%）: {'帰無仮説を棄却 → 肥料間に有意な差あり' if p_value < 0.05 else '棄却できない'}")
この分析結果から、以下のことが読み取れます。
-
左図の箱ひげ図: 肥料A、B、Cの収量には視覚的に差が見えます。赤いダイヤモンドで示された各グループの平均値は、肥料Cが最も高く、肥料Aが最も低くなっています。ただし、各グループ内のばらつき（箱の高さ）もあるため、この差が統計的に有意かどうかはF検定で判定します
-
右図のF分布: 計算されたF統計量（赤の破線）が臨界値（緑の点線）を超えており、p値が有意水準0.05を大きく下回っています。これにより帰無仮説「3種類の肥料の効果に差はない」が棄却されます
-
ANOVA表: 群間変動が群内変動に比べて相対的に大きいことが、F統計量が大きくなる原因です。群間平均平方（$MS_B$）と群内平均平方（$MS_W$）の比が大きいほど、グループ間の差が顕著であることを意味します
F検定の検出力
検出力とは
F検定の検出力（power）は、帰無仮説が偽であるとき（つまり、グループ間に本当に差があるとき）に正しく帰無仮説を棄却できる確率です。検出力が高いほど、存在する差を見逃しにくい検定であると言えます。
検出力は、効果の大きさ（グループ間の差の程度）、サンプルサイズ、有意水準、そして群内のばらつきに依存します。
以下のPythonコードで、サンプルサイズと効果量に対する検出力の変化を可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
def power_one_way_anova(n_per_group, a, effect_sizes, sigma, alpha=0.05):
"""一元配置ANOVAの検出力をシミュレーションで推定"""
n_sim = 10000
powers = []
for delta in effect_sizes:
# グループ平均: [0, delta, 2*delta, ...] をcenteredにする
means = np.linspace(0, delta, a)
means -= means.mean()
rejections = 0
for _ in range(n_sim):
groups = [np.random.normal(loc=mu, scale=sigma, size=n_per_group)
for mu in means]
all_data = np.concatenate(groups)
grand_mean = all_data.mean()
N = a * n_per_group
SS_B = sum(n_per_group * (np.mean(g) - grand_mean)**2 for g in groups)
SS_W = sum(np.sum((g - np.mean(g))**2) for g in groups)
F_stat = (SS_B / (a - 1)) / (SS_W / (N - a))
p_val = 1 - stats.f.cdf(F_stat, dfn=a-1, dfd=N-a)
if p_val < alpha:
rejections += 1
powers.append(rejections / n_sim)
return powers
a = 3 # グループ数
sigma = 8 # 群内標準偏差
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 効果量を固定し、サンプルサイズを変化
ax = axes[0]
delta_fixed = 10
n_range = [5, 10, 15, 20, 30, 50]
effect_sizes = [delta_fixed]
powers_by_n = []
for n in n_range:
p = power_one_way_anova(n, a, effect_sizes, sigma)
powers_by_n.append(p[0])
ax.plot(n_range, powers_by_n, "bo-", linewidth=2, markersize=8)
ax.axhline(0.8, color="red", linestyle="--", linewidth=1, label="Power = 0.80")
ax.set_xlabel("Sample size per group", fontsize=12)
ax.set_ylabel("Power", fontsize=12)
ax.set_title(rf"Power vs sample size ($\Delta\mu = {delta_fixed}$, $\sigma = {sigma}$)",
fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1.05)
ax.grid(True, alpha=0.3)
# (b) サンプルサイズを固定し、効果量を変化
ax = axes[1]
n_fixed = 20
effect_range = np.arange(0, 21, 2)
powers_by_effect = power_one_way_anova(n_fixed, a, effect_range, sigma)
ax.plot(effect_range, powers_by_effect, "ro-", linewidth=2, markersize=8)
ax.axhline(0.8, color="blue", linestyle="--", linewidth=1, label="Power = 0.80")
ax.set_xlabel(r"Effect size $\Delta\mu$ (max - min of group means)", fontsize=12)
ax.set_ylabel("Power", fontsize=12)
ax.set_title(rf"Power vs effect size ($n = {n_fixed}$, $\sigma = {sigma}$)",
fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1.05)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("f_distribution_power.png", dpi=150, bbox_inches="tight")
plt.show()
この検出力分析から、以下の実用的な知見が得られます。
-
左図: サンプルサイズが増加するにつれて検出力が上昇する 。効果量 $\Delta\mu = 10$、群内標準偏差 $\sigma = 8$ の場合、1グループあたり約15サンプルで検出力0.80（一般的に望ましいとされる水準）に達します。サンプルサイズが小さすぎると、差が存在しても検出できないリスクがあります
-
右図: 効果量が大きいほど検出力が高い 。これは当然ですが、グループ間の差（$\Delta\mu$）が群内のばらつき（$\sigma$）に対して十分に大きくなければ検定は無力であることを示しています。効果量が0のとき検出力は約0.05（= 有意水準）であり、正しく帰無仮説を棄却しない（第一種の過誤率）ことが確認できます
-
実験計画への示唆: 実験を開始する前に、想定される効果量と望む検出力から必要なサンプルサイズを決定すべきです。この事前分析を怠ると、差があるのに検出できない無意味な実験になってしまいます
まとめ
本記事では、F分布の定義から導出、性質、分散分析への応用までを解説しました。
- F分布は2つの独立なχ²分布の比の分布 であり、$F = \frac{Q_1/d_1}{Q_2/d_2} \sim F(d_1, d_2)$ と定義されます。2つの自由度 $(d_1, d_2)$ をパラメータとします
- 確率密度関数は変数変換とガンマ関数の積分により導出され、ベータ関数を用いたコンパクトな表現を持ちます
- 基本統計量: 期待値 $d_2/(d_2-2)$（$d_2 > 2$）、分散は $d_2 > 4$ のとき存在します。注目すべきことに、期待値は分子の自由度 $d_1$ に依存しません
- 逆数の性質 $1/F(d_1, d_2) \sim F(d_2, d_1)$ とt分布との関係 $t^2(d) \sim F(1, d)$ は特に重要です
- 分散分析（ANOVA） では、群間変動と群内変動の比としてF統計量が構成され、グループ間の平均値に差があるかを検定します
- 検出力分析により、実験計画におけるサンプルサイズの決定ができます
次のステップとして、以下の記事も参考にしてください。