サンプルサイズの決め方6章 2つの母分散の比の検定

はじめに

このブログ記事では、永田靖著『サンプルサイズの決め方』（朝倉書店）の6章を参考に、2つの母分散の比の検定、検出力の計算方法、適切なサンプルサイズの設計方法について解説します。これら統計学の基本的な概念をPythonを用いて実装しながら確認していきます。

2つの母分散の比の検定

母分散を検定する手順は書籍では以下のように紹介されています。

• 手順1 仮説の設定

  • 帰無仮説 $H_0: \sigma^2_1 = \sigma_2^2$
  • 対立仮説:
    • (1)$H_1: \sigma^2_1 \neq \sigma_2^2$
    • (2)$H_1: \sigma^2_1 > \sigma_2^2$
    • (3)$H_1: \sigma^2_1 < \sigma_2^2$

• 手順2 有意水準の設定

  • 有意水準 $\alpha$ の設定 (通常は $\alpha = 0.05$ とする)

• 手順3 棄却域の設定

  • 棄却域 $R$ の設定:

    • 対立仮説 (1)の場合:
    $$\begin{aligned} R: &F_0 \leq F(\phi_1, \phi_2, 1 - \alpha/2) = \frac{1}{F(\phi_2,\phi_1;\alpha/2)}\\ &F_0 \geq F(\phi_1, \phi_2; \alpha/2) \end{aligned}$$
    • 対立仮説 (2)の場合:
    $$\begin{aligned} R: F_0 \geq F(\phi_1, \phi_2, \alpha) \end{aligned}$$
    • 対立仮説 (3)の場合:
    $$\begin{aligned} R: F_0 \leq F(\phi_1, \phi_2, 1 - \alpha) = \frac{1}{F(\phi_2,\phi_1;\alpha)} \end{aligned}$$

• 手順4 検定統計量の計算

  • データ $x_{11}, x_{12}, \ldots, x_{1n_1}$および$x_{21}, x_{22}, \ldots, x_{2n_2}$をとり、検定統計量 $F_0$ の値を計算する。
  $$F_0 = \frac{V_1}{V_2} \\[1em] V_1 = \frac{S_1}{n_1-1} = \frac{\Sigma^{n_1}_{i=1}(x_{1i}-\bar{x}_1)^2}{n_1-1} = \frac{\Sigma^{n_1}_{i=1}x_{1i}^2-(\Sigma x_{1i})^2/n_1}{n_1-1} \\[1em] V_2 = \frac{S_2}{n_2-1} = \frac{\Sigma^{n_1}_{i=1}(x_{2i}-\bar{x}_2)^2}{n_2-1} = \frac{\Sigma^{n_1}_{i=1}x_{2i}^2-(\Sigma x_{2i})^2/n_2}{n_2-1} \\[1em] \phi_1 = n_1 - 1 \\[1em] \phi_2 = n_2 - 1 \\[1em]$$

• 手順5 判定

  • $F_0$ が手順3で設定した棄却域 $R$ にあれば有意と判定し、$H_0$ を棄却する。

この検定手順に則って検定を行う関数をpythonで実装しました。

import numpy as np
from scipy import stats

def f_test(data1, data2, alpha=0.05, test_type='two-sided'):
    """
    パラメータ:
    data1, data2 : 配列のようなオブジェクト
        比較する2つのサンプル
    alpha : float
        有意水準（デフォルトは0.05）
    test_type : {'two-sided', 'greater', 'less'}
    """
    
    # 検定統計量の計算
    n1, n2 = len(data1), len(data2)
    df1, df2 = n1 - 1, n2 - 1
    v1 = np.var(data1, ddof=1)
    v2 = np.var(data2, ddof=1)
    f_statistic = v1 / v2
    
    # 棄却域の設定と判定
    if test_type == 'two-sided':
        critical_value_lower = stats.f.ppf(alpha/2, df1, df2)
        critical_value_upper = stats.f.ppf(1 - alpha/2, df1, df2)
        reject = f_statistic <= critical_value_lower or f_statistic >= critical_value_upper
        critical_values = (critical_value_lower, critical_value_upper)
    elif test_type == 'greater':
        critical_value = stats.f.ppf(1 - alpha, df1, df2)
        reject = f_statistic >= critical_value
        critical_values = (critical_value,)
    elif test_type == 'less':
        critical_value = stats.f.ppf(alpha, df1, df2)
        reject = f_statistic <= critical_value
        critical_values = (critical_value,)
    else:
        raise ValueError("test_typeは'two-sided'、'greater'または'less'でなければなりません")
    
    return reject, f_statistic, critical_values

この関数を用いて練習問題を解いていきます。

$$\begin{aligned} &\text{問題6.1} \\ &n_1 = 10\text{、} n_2 = 11 \text{ の下記のデータに基づいて、} H_0：\sigma^2_1 = \sigma^2_2\text{、} H_1：\sigma^2_1 \neq \sigma^2_2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ1：} 4.2, 2.8, 5.3, 3.5, 4.5, 2.9, 4.8, 5.9, 4.6, 5.3 \\ &\text{データ2：} 7.3, 6.9, 7.3, 8.4, 8.5, 8.1, 5.3, 7.6, 6.2, 5.9, 6.7 \end{aligned}$$

# 問題6.1
print("問題6.1の解答:")
data1_6_1 = [4.2, 2.8, 5.3, 3.5, 4.5, 2.9, 4.8, 5.9, 4.6, 5.3]
data2_6_1 = [7.3, 6.9, 7.3, 8.4, 8.5, 8.1, 5.3, 7.6, 6.2, 5.9, 6.7]

result_6_1 = f_test(data1_6_1, data2_6_1, alpha=0.05, test_type='two-sided')
print(f"帰無仮説の棄却: {result_6_1[0]}")
print(f"F統計量: {result_6_1[1]}")
print(f"臨界値: {result_6_1[2]}")
print()

結果

問題6.1の解答:
帰無仮説の棄却: False
F統計量: 1.0139362266151635
臨界値: (0.2522790156161077, 3.7789626340915725)

F統計量(1.0139)が両側の臨界値(0.2523, 3.7790)の間にあるため、帰無仮説を棄却できません。つまり、2つのグループの分散に有意な差があるとは言えません。

$$\begin{aligned} &\text{問題6.2} \\ &n_1 = 9\text{、} n_2 = 8 \text{ の下記のデータに基づいて、} H_0：\sigma^2_1 = \sigma^2_2\text{、} H_1：\sigma^2_1 > \sigma^2_2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ1：} 13.8, 10.2, 8.7, 7.9, 13.0, 11.6, 10.5, 8.7, 7.1 \\ &\text{データ2：} 10.3, 10.0, 9.8, 8.8, 10.5, 9.4, 10.6, 10.1 \end{aligned}$$

# 問題6.2
print("問題6.2の解答:")
data1_6_2 = [13.8, 10.2, 8.7, 7.9, 13.0, 11.6, 10.5, 8.7, 7.1]
data2_6_2 = [10.3, 10.0, 9.8, 8.8, 10.5, 9.4, 10.6, 10.1]

result_6_2 = f_test(data1_6_2, data2_6_2, alpha=0.05, test_type='greater')
print(f"帰無仮説の棄却: {result_6_2[0]}")
print(f"F統計量: {result_6_2[1]}")
print(f"臨界値: {result_6_2[2]}")
print()

結果

問題6.2の解答:
帰無仮説の棄却: True
F統計量: 14.673945409429296
臨界値: (3.725725317122702,)

F統計量(14.6739)が臨界値(3.7257)よりも大きいため、帰無仮説を棄却します。つまり、データ1の分散がデータ2の分散より有意に大きいと結論づけられます。

$$\begin{aligned} &\text{問題6.3} \\ &n_1 = 8\text{、} n_2 = 10 \text{ の下記のデータに基づいて、} H_0：\sigma^2_1 = \sigma^2_2\text{、} H_1：\sigma^2_1 < \sigma^2_2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ1：} 21, 19, 16, 19, 22, 18, 20, 21 \\ &\text{データ2：} 19, 22, 21, 22, 25, 19, 24, 23, 19, 22 \end{aligned}$$

# 問題6.3
print("問題6.3の解答:")
data1_6_3 = [21, 19, 16, 19, 22, 18, 20, 21]
data2_6_3 = [19, 22, 21, 22, 25, 19, 24, 23, 19, 22]

result_6_3 = f_test(data1_6_3, data2_6_3, alpha=0.05, test_type='less')
print(f"帰無仮説の棄却: {result_6_3[0]}")
print(f"F統計量: {result_6_3[1]}")
print(f"臨界値: {result_6_3[2]}")

結果

問題6.3の解答:
帰無仮説の棄却: False
F統計量: 0.8274398868458276
臨界値: (0.27198489991198765,)

F統計量(0.8274)が臨界値(0.2720)よりも大きいため、帰無仮説を棄却できません。つまり、データ1の分散がデータ2の分散より小さいとは言えません。

検出力の計算

検出力を計算する手順は書籍では以下のように紹介されています。

$$\begin{aligned} \text{両側検定 }H_1: \sigma^2_1 \neq \sigma_2^2 \text{}\\[1em] &Pr\left(F \leq \frac{1}{\Delta^2 F(\phi_2, \phi_1;\alpha/2)}\right) + Pr\left(F \geq \frac{F(\phi_1,\phi_2;\alpha/2)}{\Delta^2}\right)\\ \text{片側検定 }H_1: \sigma^2_1 > \sigma_2^2\text{} \\[1em] &Pr\left(F \geq \frac{F(\phi_1,\phi_2;\alpha)}{\Delta^2}\right) \\ \text{片側検定 }H_1: \sigma^2_1 < \sigma_2^2\text{} \\ &Pr\left(F \leq \frac{1}{\Delta^2F(\phi_2, \phi_1;\alpha)}\right) \end{aligned}$$

以下のプログラムは検出力の計算をpythonで実装したものとなっています。
このcalculate_power関数を用いて、練習問題を解いてきます。

def calculate_power(n, alpha, delta, test_type='two-sided'):
    """
    n: サンプルサイズ
    alpha: 有意水準
    delta: Δ 
    test_type: 'two-sided', 'greater' or 'less'
    """
    df = n - 1
    
    if test_type == 'two-sided':
        chi2_lower = stats.chi2.ppf(alpha/2, df)
        chi2_upper = stats.chi2.ppf(1 - alpha/2, df)
        power = stats.chi2.cdf(chi2_lower / delta**2, df) + \
                (1 - stats.chi2.cdf(chi2_upper / delta**2, df))
    elif test_type == 'greater':
        chi2_crit = stats.chi2.ppf(1 - alpha, df)
        power = 1 - stats.chi2.cdf(chi2_crit / delta**2, df)
    elif test_type == 'less':
        chi2_crit = stats.chi2.ppf(alpha, df)
        power = stats.chi2.cdf(chi2_crit / delta**2, df)
    else:
        raise ValueError("test_type must be 'two-sided', 'greater' or 'less'")
    
    return power

$$\begin{aligned} &\text{問題6.4} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma_1^2 = \sigma_2^2\text{、} H_1：\sigma_1^2 \neq \sigma_2^2 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n_1 = 10\text{、} n_2 = 11\text{、}\Delta = \sigma_1/\sigma_2 = 0.5 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n_1 = 10\text{、} n_2 = 11\text{、}\Delta = \sigma_1/\sigma_2 = 2.0 \text{ の場合に検出力を求めよ。} \end{aligned}$$

# 問題6.4
print("問題6.4の解答:")
n1_6_4, n2_6_4 = 10, 11
alpha_6_4 = 0.05

# (1) Δ = 0.5 の場合
delta_6_4_1 = 0.5
power_6_4_1 = calculate_power(n1_6_4, n2_6_4, delta_6_4_1, alpha_6_4, 'two-sided')
print(f"(1) Δ = 0.5 の場合の検出力: {power_6_4_1:.4f}")

# (2) Δ = 2.0 の場合
delta_6_4_2 = 2.0
power_6_4_2 = calculate_power(n1_6_4, n2_6_4, delta_6_4_2, alpha_6_4, 'two-sided')
print(f"(2) Δ = 2.0 の場合の検出力: {power_6_4_2:.4f}")

結果

問題6.4の解答:
(1) Δ = 0.5 の場合の検出力: 0.5101
(2) Δ = 2.0 の場合の検出力: 0.5295

$Δ$の値を0.5から2.0に変更しても、検出力はそれほど変わっていません（0.5101から0.5295）。

$$\begin{aligned} &\text{問題6.5} \\ &\text{帰無仮説と対立仮説を } H_0：\mu_1 = \mu_2\text{、} H_1：\mu_1 > \mu_2 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n_1 = 9\text{、} n_2 = 8\text{、}\Delta = \sigma_1/\sigma_2 = 2.0 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n_1 = 9\text{、} n_2 = 8\text{、}\Delta = \sigma_1/\sigma_2 = 2.5 \text{ の場合に検出力を求めよ。} \end{aligned}$$

# 問題6.5
print("問題6.5の解答:")
n1_6_5, n2_6_5 = 9, 8
alpha_6_5 = 0.05

# (1) Δ = 2.0 の場合
delta_6_5_1 = 2.0
power_6_5_1 = calculate_power(n1_6_5, n2_6_5, delta_6_5_1, alpha_6_5, 'greater')
print(f"(1) Δ = 2.0 の場合の検出力: {power_6_5_1:.4f}")

# (2) Δ = 2.5 の場合
delta_6_5_2 = 2.5
power_6_5_2 = calculate_power(n1_6_5, n2_6_5, delta_6_5_2, alpha_6_5, 'greater')
print(f"(2) Δ = 2.5 の場合の検出力: {power_6_5_2:.4f}")

結果

問題6.5の解答:
(1) Δ = 2.0 の場合の検出力: 0.5440
(2) Δ = 2.5 の場合の検出力: 0.7582

$Δ$の値を2.0から2.5に増やすことで、検出力が0.5440から0.7582に大幅に上昇しました。これは、$Δ$が大きいほど、その差を統計的に検出しやすくなることを示しています。

$$\begin{aligned} &\text{問題6.6} \\ &\text{帰無仮説と対立仮説を } H_0：\mu_1 = \mu_2\text{、} H_1：\mu_1 < \mu_2 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n_1 = 8\text{、} n_2 = 10\text{、}\Delta = \sigma_1/\sigma_2 = 1/2 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n_1 = 8\text{、} n_2 = 10\text{、}\Delta = \sigma_1/\sigma_2 = 1/3 \text{ の場合に検出力を求めよ。} \end{aligned}$$

# 問題6.6
print("問題6.6の解答:")
n1_6_6, n2_6_6 = 8, 10
alpha_6_6 = 0.05

# (1) Δ = 1/2 の場合
delta_6_6_1 = 1/2
power_6_6_1 = calculate_power(n1_6_6, n2_6_6, delta_6_6_1, alpha_6_6, 'less')
print(f"(1) Δ = 1/2 の場合の検出力: {power_6_6_1:.4f}")

# (2) Δ = 1/3 の場合
delta_6_6_2 = 1/3
power_6_6_2 = calculate_power(n1_6_6, n2_6_6, delta_6_6_2, alpha_6_6, 'less')
print(f"(2) Δ = 1/3 の場合の検出力: {power_6_6_2:.4f}")

結果

問題6.6の解答:
(1) Δ = 1/2 の場合の検出力: 0.5576
(2) Δ = 1/3 の場合の検出力: 0.8944

$Δ$の値を1/2から1/3に減少させることで、検出力が0.5576から0.8944に大きく上昇しました。

サンプルサイズの設計方法

サンプルサイズの設計方法は書籍では以下のように紹介されています。

$$\begin{align*} & \text{両側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 \neq \sigma_0^2) \\[0.5em] & \quad \text{(a) } \Delta_{01} > 1 \text{ の場合:} \quad n \approx 1+\left(\frac{z_{\alpha/2} - z_{1-\beta}}{ln\Delta_{01}}\right)^2 \\[1em] & \quad \text{(b) } \Delta_{02} < 1 \text{ の場合:} \quad n \approx 1+\left(\frac{z_{\alpha/2} - z_{1-\beta}}{ln\Delta_{02}}\right)^2 \\[2em] & \text{片側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 > \sigma_0^2) \\[0.5em] & \quad n \approx 1+\left(\frac{z_{\alpha} - z_{1-\beta}}{ln\Delta_0}\right)^2 \\[2em] & \text{片側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 < \sigma_0^2) \\[0.5em] & \quad n \approx 1+\left(\frac{z_{\alpha} - z_{1-\beta}}{ln\Delta_0}\right)^2 \end{align*}$$

上記の式通りにサンプルサイズを計算する関数をpythonで実装しました。
これと検出力の計算の部分で実装したcalculate_powerを利用して、練習問題を解いていきます。

def calculate_sample_size(alpha, beta, delta, test_type='two-sided'):
    z_alpha = stats.norm.ppf(1 - alpha/2) if test_type == 'two-sided' else stats.norm.ppf(1 - alpha)
    z_beta = stats.norm.ppf(beta)
    
    n = 1 + ((z_alpha - z_beta) / np.log(delta))**2
    return np.ceil(n)

def find_optimal_sample_size(alpha, beta, delta, test_type='two-sided'):
    if test_type == 'two-sided':
        n1 = int(calculate_sample_size(alpha, 1-beta, delta, test_type))
        n2 = int(calculate_sample_size(alpha, 1-beta, 1/delta, test_type))
        n = max(n1, n2)
        
        power = calculate_power(n, n, delta, alpha, test_type)
        while power < beta:
            n += 1
            power = calculate_power(n, n, delta, alpha, test_type)

        return n, power
    else:
        n = int(calculate_sample_size(alpha, 1-beta, delta, test_type))
        
        power = calculate_power(n, n, delta, alpha, test_type)
        while power < beta:
            n += 1
            power = calculate_power(n, n, delta, alpha, test_type)

        return n, power

$$\begin{aligned} &\text{問題6.7} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma_1^2 = \sigma_2^2\text{、} H_1：\sigma_1^2 \neq \sigma_2^2 \text{ と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma_1/\sigma_2 \geq 2.0 \text{ または } \Delta = \sigma_1/\sigma_2 \leq 2/5 \text{ のとき検出力が } 0.90 \text{ となるように} \\ &\quad\text{サンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma_1/\sigma_2 \geq 3.0 \text{ または } \Delta = \sigma_1/\sigma_2 \leq 1/4 \text{ のとき検出力が } 0.95 \text{ となるように} \\ &\quad\text{サンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題6.7
print("問題6.7:")
n_6_7_1, power_6_7_1 = find_optimal_sample_size(0.05, 0.90, 2.0, 'two-sided')
print(f"(1) サンプルサイズ: {n_6_7_1}, 検出力: {power_6_7_1:.4f}")
n_6_7_2, power_6_7_2 = find_optimal_sample_size(0.05, 0.95, 3.0, 'two-sided')
print(f"(2) サンプルサイズ: {n_6_7_2}, 検出力: {power_6_7_2:.4f}")

結果

問題6.7:
(1) サンプルサイズ: 24, 検出力: 0.9019
(2) サンプルサイズ: 13, 検出力: 0.9535

両側検定では、$Δ$が1より大きい場合と小さい場合の両方を考慮する必要があります。ここでは、$Δ=2.0（または1/2）$の場合より$Δ=3.0（または1/3）$の場合の方が、同じ検出力を得るために必要なサンプルサイズが小さくなっています。

$$\begin{aligned} &\text{問題6.8} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma_1^2 = \sigma_2^2\text{、} H_1：\sigma_1^2 > \sigma_2^2 \text{ と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma_1/\sigma_2 \geq 2.0 \text{ のとき検出力が } 0.90 \text{ となるようにサンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma_1/\sigma_2 \geq 2.5 \text{ のとき検出力が } 0.90 \text{ となるようにサンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題6.8
print("\n問題6.8:")
n_6_8_1, power_6_8_1 = find_optimal_sample_size(0.05, 0.90, 2.0, 'greater')
print(f"(1) サンプルサイズ: {n_6_8_1}, 検出力: {power_6_8_1:.4f}")
n_6_8_2, power_6_8_2 = find_optimal_sample_size(0.05, 0.90, 2.5, 'greater')
print(f"(2) サンプルサイズ: {n_6_8_2}, 検出力: {power_6_8_2:.4f}")

結果

問題6.8:
(1) サンプルサイズ: 20, 検出力: 0.9044
(2) サンプルサイズ: 13, 検出力: 0.9211

$Δ$の値を2.0から2.5に増やすことで、同じ検出力（0.90）を得るために必要なサンプルサイズが20から13に減少しました。

$$\begin{aligned} &\text{問題6.9} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma_1^2 = \sigma_2^2\text{、} H_1：\sigma_1^2 < \sigma_2^2 \text{ と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma_1/\sigma_2 \leq 1/2 \text{ のとき検出力が } 0.90 \text{ となるようにサンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma_1/\sigma_2 \leq 1/4 \text{ のとき検出力が } 0.95 \text{ となるようにサンプルサイズ } n(= n_1 = n_2) \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題6.9
print("\n問題6.9:")
n_6_9_1, power_6_9_1 = find_optimal_sample_size(0.05, 0.90, 1/2, 'less')
print(f"(1) サンプルサイズ: {n_6_9_1}, 検出力: {power_6_9_1:.4f}")
n_6_9_2, power_6_9_2 = find_optimal_sample_size(0.05, 0.95, 1/4, 'less')
print(f"(2) サンプルサイズ: {n_6_9_2}, 検出力: {power_6_9_2:.4f}")

結果

問題6.9:
(1) サンプルサイズ: 20, 検出力: 0.9044
(2) サンプルサイズ: 8, 検出力: 0.9617

$Δ$の値を1/2から1/4に減少させ、目標とする検出力を0.90から0.95に上げたにもかかわらず、必要なサンプルサイズは20から8に減少しました。これは、$Δ$が極めて大きい場合、高い検出力を少ないサンプル数で達成できることを示しています。

最後に

本ブログでは、永田靖著『サンプルサイズの決め方』の第4章「１つの母分散の検定」の内容を、Pythonを用いて実践的に解説しました。ここで学んだ主要なポイントは以下の通りです：

1. 母分散の検定の基本概念と手順
2. 検出力の計算方法とその重要性
3. 適切なサンプルサイズの設計方法

次回のブログでは、『サンプルサイズの決め方』の第7章の内容を取り上げる予定です。このブログを通して少しでもサンプルサイズの設計方法について理解が進めば幸いです。